|| || || || || || || ||

---

IntWeibullNgbPnLoop.h

Go to the documentation of this file.

#ifndef Impala_Core_Feature_IntWeibullNgbPnLoop_h
#define Impala_Core_Feature_IntWeibullNgbPnLoop_h

#include "Core/Array/Pattern/Categories.h"
#include "Core/Array/Pattern/Cnum.h"
#include "Core/Array/Element/E1Cast.h"
#include "Core/Array/Element/Vec3Real64.h"
#include "Core/Array/Arrays.h"
//#include <cmath>
#include "limits.h"
#include "float.h"
#include <errno.h>

namespace Impala
{
namespace Core
{
namespace Feature
{

/* constants: pi, Euler's constant, and log(2) */
#ifndef M_PI
#define M_PI 3.1415926535897932384626433832795029
#endif
#ifndef MMGAMMA
#define MMGAMMA 0.5772156649015328606065120900824024
#endif
#ifndef M_LN2
#define M_LN2 0.6931471805599453094172321214581766
#endif


template<class DstArrayT, class SrcArrayT>
class IntWeibullNgbPnLoop
{
public:

    typedef Array::Pattern::TagLoop IteratorCategory;

    typedef Array::Pattern::TagNPhase PhaseCategory;

    typedef typename DstArrayT::StorType DstStorType;
    typedef typename DstArrayT::ArithType DstArithType;
    typedef typename SrcArrayT::StorType SrcStorType;
    typedef typename SrcArrayT::ArithType SrcArithType;


    IntWeibullNgbPnLoop(bool verbose, int ngbWidth, int ngbHeight)
    {
        mNgbWidth = ngbWidth;
        mNgbHeight = ngbHeight;
        mVerbose = verbose;
        mFrac = mNgbWidth*mNgbHeight;
        mGamma = 0.01;
        mGammaNext = 0.0;
        mLast = false;
        mData = new SrcArrayT(mFrac, 1, 0, 0);
        mEcdf = new SrcArrayT(mFrac, 1, 0, 0);
    }

    ~IntWeibullNgbPnLoop()
    {
      delete mData;
      delete mEcdf;
    }

    int
    Width()
    {
        return mNgbWidth;
    }

    int
    Height()
    {
        return mNgbHeight;
    }

    void
    Init(int phase, int x, int y, SrcArithType value)
    {
      if(mLast) return;
      int ph = phase%2;
      //        if (mVerbose)
      //  std::cout << "  WeibullNgbPnLoop::init(" << phase << ","
      //            << x << "," << y << "," << value << ") " << std::endl;

      if (phase == 1)
      {
        //      std::cout << "int: x = " << x << " y = " << y << std::endl; 
        //      if(mGamma > 5) mGamma = 0.01;
        mGamma = 0.01; 
        mGammaNext = 0.0; 
        mPrecision = 0.005; 
        mMaxNumIters = 20;
        mIter = 1;
        mMaxGamma = 5; 
      }
      mPhase = phase;
      if (ph)
      {
        mSumDataPowerGamma = 0.0;
        mSumDataLog = 0.0;
        mF = 0.0;
        mFprime = 0.0;
      }
      
    }

    void
    NextEl(int x, int y, SrcArithType value)
    {
      SrcArithType dataPowerGamma;
      SrcArithType sumItem;
      SrcArithType data = fabs(value) + 0.000001; 
      dataPowerGamma = pow(data, mGamma);
      int phase = mPhase%2;
        
      if(mVerbose)
        if (mGamma == 0.01 && phase == 0)
          std::cout << value << std::endl;
      //      if (mVerbose)
      //        std::cout << "      WeibullNgbPnLoop::next(" << x << "," << y << ","<< value << ") " << std::endl;

      if(mLast)
      {
        mData->SetValue(fabs(value), x*mNgbWidth + y, 0);
        return;
      }
      switch (phase)
      {
      case 1 : 
        mSumDataPowerGamma += dataPowerGamma;
        mSumDataLog += dataPowerGamma*log(data);
        break;
        
      case 0 :
        sumItem = dataPowerGamma/mSumDataPowerGamma; 
        mF += (sumItem)*log(mFrac*sumItem);
        mFprime += (log(dataPowerGamma/mSumDataPowerGamma) + 1)*
          (dataPowerGamma*(log(data)*mSumDataPowerGamma - mSumDataLog))/
          (mSumDataPowerGamma*mSumDataPowerGamma);
        break;
      }
        
    }

    void
    Done(int phase)
    {
      int ph = phase%2;
      if (mVerbose)
        std::cout << "done: phase = " << phase <<  std::endl;
      if(!ph)
      {
        if (mVerbose)
          std::cout << "mGamma = " << mGamma << std::endl;
        
        mFprime = (1/pow(mGamma,2))*(1 - log(mGamma) - DiGamma(1/mGamma) - 
                                     (1/mGamma)*TriGamma(1/mGamma) + mF) -(1/mGamma)*mFprime;
        mF = 1 + (1/mGamma)*(log(mGamma) + DiGamma(1/mGamma) - mF);
        mGammaNext = mGamma - mF/mFprime;
        if (mVerbose)
          std::cout << "mGammaNext = " << mGammaNext << std::endl;
      }
      
    }

    bool
    HasNextPhase(int lastPhase)
    {
      if(mLast)
      {
        mLast = false;
        return false;
      }
      int ph = lastPhase%2;  
      if (mVerbose)
        std::cout << "  WeibullNgbPnLoop::hasNext(" << lastPhase << ") "
                  << std::endl;

      if (ph)
        return true;

      else
      {
        if (fabs(mGamma - mGammaNext) < mPrecision || mIter > mMaxNumIters || 
            mGammaNext > mMaxGamma || isnan(mGammaNext) || isinf(mGammaNext) || mGammaNext < 0)
        {         
        if (mVerbose)
          {
            std::cout << "no next phase! " << "fabs(mGamma - mGammaNext) = " << fabs(mGamma - mGammaNext)  << std::endl;
            std::cout << "mIter = " << mIter << " maxIter =  " << mMaxNumIters << std::endl; 
            std::cout << "maxGamma = " << mMaxGamma << std::endl;
            std::cout << "mGamma = " << mGamma << " mGammaNext = " << mGammaNext << std::endl;
          }
        //      std::cout << "mIter = " << mIter << std::endl; 
          mIter = 1;
          mLast = true;
          return true;
        }
        else 
        {
          if (mVerbose)
            std::cout << "there is next phase: iteration = " << mIter << std::endl;
          mGamma = mGammaNext;
          mIter++;
          return true;
        }
      }
    }

    DstArithType
    Result()
    {
      double Beta = pow(mSumDataPowerGamma/mFrac, (1/mGamma));
      if(Beta == 0)
        Beta = 0.000001;
      Core::Array::Sort(mEcdf, mData);
      Core::Array::Set(mData, mEcdf);
      Core::Array::MulVal(mEcdf, mEcdf, 1/mData->Value(mFrac-1, 0));
      double AD = 0.0;
      double F;
      double dF;
      for (int i=0; i<mFrac; i++)
      {
        F = 1 - exp(-(1/mGamma)*pow(mData->Value(i, 0)/Beta, mGamma));
        
        if(F !=0 && F != 1)
        {
          dF = (1/Beta)*pow( mData->Value(i, 0)/Beta, mGamma-1 )*F;
          AD += pow((mEcdf->Value(i, 0) - F), 2)*dF/(F*(1-F));
        }
      }
      AD = AD/mFrac;
      if (mVerbose)
        std::cout << "result: Gamma = " << mGamma  << " Beta = " << Beta << " AD = " << AD  << std::endl;
      return Core::Array::Element::Vec3Real64(Beta, mGamma, AD);
   }

    double
    gammaFunc(double xx) const
    {
        double x,y,tmp,ser;
        static double cof[6] =
            {
                76.18009172947146, -86.50532032941677,
                24.01409824083091, -1.231739572450155,
                0.1208650973866179e-2, -0.5395239384953e-5
            };
        
        int j;
        y=x=xx;
        tmp=x+5.5;
        tmp -= (x+0.5)*log(tmp);
        ser=1.000000000190015;
        for (j=0;j<=5;j++)
            ser += cof[j]/++y;
        xx = -tmp+log(2.5066282746310005*ser/x);
        return exp(xx);
    }

    double
    DiGamma(double x) const
    {
        /* force into the interval 1..3 */
        if( x < 0.0 )
            return DiGamma(1.0-x)+M_PI/tan(M_PI*(1.0-x)) ;  /* refection formula */
        else if( x < 1.0 )
            return DiGamma(1.0+x)-1.0/x ;
        else if ( x == 1.0)
            return -MMGAMMA ;
        else if ( x == 2.0)
            return 1.0-MMGAMMA ;
        else if ( x == 3.0)
            return 1.5-MMGAMMA ;
        else if ( x > 3.0)
            /* duplication formula */
            return 0.5*(DiGamma(x/2.0)+DiGamma((x+1.0)/2.0))+M_LN2 ;
        else
        {
            static double Kncoe[] =
                { .30459198558715155634315638246624251,
                  .72037977439182833573548891941219706, -.12454959243861367729528855995001087,
                  .27769457331927827002810119567456810e-1, -.67762371439822456447373550186163070e-2,
                  .17238755142247705209823876688592170e-2, -.44817699064252933515310345718960928e-3,
                  .11793660000155572716272710617753373e-3, -.31253894280980134452125172274246963e-4,
                  .83173997012173283398932708991137488e-5, -.22191427643780045431149221890172210e-5,
                  .59302266729329346291029599913617915e-6, -.15863051191470655433559920279603632e-6,
                  .42459203983193603241777510648681429e-7, -.11369129616951114238848106591780146e-7,
                  .304502217295931698401459168423403510e-8, -.81568455080753152802915013641723686e-9,
                  .21852324749975455125936715817306383e-9, -.58546491441689515680751900276454407e-10,
                  .15686348450871204869813586459513648e-10, -.42029496273143231373796179302482033e-11,
                  .11261435719264907097227520956710754e-11, -.30174353636860279765375177200637590e-12,
                  .80850955256389526647406571868193768e-13, -.21663779809421233144009565199997351e-13,
                  .58047634271339391495076374966835526e-14, -.15553767189204733561108869588173845e-14,
                  .41676108598040807753707828039353330e-15, -.11167065064221317094734023242188463e-15 } ;

            double Tn_1 = 1.0 ; /* T_{n-1}(x), started at n=1 */
            double Tn = x-2.0 ; /* T_{n}(x) , started at n=1 */
            double resul = Kncoe[0] + Kncoe[1]*Tn ;

            x -= 2.0 ;

            for(int n = 2 ; n < sizeof(Kncoe)/sizeof(double) ;n++)
            {
                const double Tn1 = 2.0 * x * Tn - Tn_1 ;    /* Chebyshev recursion, Eq. 22.7.4 Abramowitz-Stegun */
                resul += Kncoe[n]*Tn1 ;
                Tn_1 = Tn ;
                Tn = Tn1 ;
            }
            return resul ;
        }
    } //  DiGamma(double x)

#define ISNAN(x) (isnan(x) !=0)
#define ML_POSINF (1.0 / 0.0)
#define ML_NEGINF ((-1.0) / 0.0)
#define ML_NAN (0.0 / 0.0)
#define n_max (100)
#define M_LOG10_2 0.30102999566398119521

int Rf_i1mach(int i)
{
  switch(i) 
  {
    
  case  1: return 5;
  case  2: return 6;
  case  3: return 0;
  case  4: return 0;
    
  case  5: return CHAR_BIT * sizeof(int);
  case  6: return sizeof(int)/sizeof(char);
    
  case  7: return 2;
  case  8: return CHAR_BIT * sizeof(int) - 1;
  case  9: return INT_MAX;
    
  case 10: return FLT_RADIX;
    
  case 11: return FLT_MANT_DIG;
  case 12: return FLT_MIN_EXP;

  case 13: return FLT_MAX_EXP;
    
  case 14: return DBL_MANT_DIG;
        case 15: return DBL_MIN_EXP;
        case 16: return DBL_MAX_EXP;

  default: return 0;
  }
}
 
double Rf_d1mach(int i)
{
  switch(i) 
  {
    
  case 1: return DBL_MIN;
  case 2: return DBL_MAX;
    
  case 3: /* = FLT_RADIX  ^ - DBL_MANT_DIG
             for IEEE:  = 2^-53 = 1.110223e-16 = .5*DBL_EPSILON */
    return 0.5*DBL_EPSILON;
    
  case 4: /* = FLT_RADIX  ^ (1- DBL_MANT_DIG) =
             for IEEE:  = 2^-52 = DBL_EPSILON */
    return DBL_EPSILON;
    
  case 5: return M_LOG10_2;
    
  default: return 0.0;
  }
}
 
double fmax2(double x, double y)
{
#ifdef IEEE_754
  if (ISNAN(x) || ISNAN(y))
    return x + y;
#endif
  return (x < y) ? y : x;
}

double fmin2(double x, double y)
{
#ifdef IEEE_754
  if (ISNAN(x) || ISNAN(y))
    return x + y;
#endif
  return (x < y) ? x : y;
}

int imin2(int x, int y)
{
  return (x < y) ? x : y;
}

void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr)
{
  const static double bvalues[] = {/* Bernoulli Numbers */
          1.00000000000000000e+00,
          -5.00000000000000000e-01,
          1.66666666666666667e-01,
          -3.33333333333333333e-02,
          2.38095238095238095e-02,
          -3.33333333333333333e-02,
          7.57575757575757576e-02,
          -2.53113553113553114e-01,
          1.16666666666666667e+00,
          -7.09215686274509804e+00,
          5.49711779448621554e+01,
          -5.29124242424242424e+02,
          6.19212318840579710e+03,
          -8.65802531135531136e+04,
          1.42551716666666667e+06,
          -2.72982310678160920e+07,
          6.01580873900642368e+08,
          -1.51163157670921569e+10,
          4.29614643061166667e+11,
          -1.37116552050883328e+13,
          4.88332318973593167e+14,
          -1.92965793419400681e+16
        };

        int i, j, k, mm, mx, nn, np, nx, fn;
        double arg, den, elim, eps, fln, fx, rln, rxsq,
          r1m4, r1m5, s, slope, t, ta, tk, tol, tols, tss, tst,
          tt, t1, t2, wdtol, xdmln, xdmy, xinc, xln = 0.0 /* -Wall */, 
          xm, xmin, xq, yint;
        double trm[23], trmr[n_max + 1];

        *ierr = 0;
        if (n < 0 || kode < 1 || kode > 2 || m < 1) {
          *ierr = 1;
          return;
        }
        if (x <= 0.) {
          /* useAbramowitz & Stegun 6.4.7 "Reflection Formula"
           *psi(k, x) = (-1)^k psi(k, 1-x)-  pi^{n+1} (d/dx)^n cot(x)
           */
          if (x == (long)x) {
            /* non-positive integer : +Inf or NaN depends on n */
            for(j=0; j < m; j++) /* k = j + n : */
              ans[j] = ((j+n) % 2) ? ML_POSINF : ML_NAN;
            return;
          }
          dpsifn(1. - x, n, /*kode = */ 1, m, ans, nz, ierr);
          /* ans[j] == (-1)^(k+1) / gamma(k+1) * psi(k, 1 - x)
           *     for j = 0:(m-1) ,k = n + j
           */

          /* Cheat for now: only work for m = 1, n in {0,1,2,3} : */
          if(m > 1 || n > 3) {/* doesn't happen for digamma() .. pentagamma() */
            /* not yet implemented */
            *ierr = 4; return;
          }
          x *= M_PI; /* pi * x */
          if (n == 0)
            tt = cos(x)/sin(x);
          else if (n == 1)
            tt = -1/pow(sin(x),2);
          else if (n == 2)
            tt = 2*cos(x)/pow(sin(x),3);
          else if (n == 3)
            tt = -2*(2*pow(cos(x),2) + 1)/pow(sin(x),4);
          else /* can not happen! */
            tt = ML_NAN;
          /* end cheat */

          s = (n % 2) ? -1. : 1.;/* s = (-1)^n */
          /* t := pi^(n+1) * d_n(x) / gamma(n+1), where
           *   d_n(x) := (d/dx)^n cot(x)*/
          t1 = t2 = s = 1.;
          for(k=0, j=k-n; j < m; k++, j++, s = -s) {
            /* k == n+j , s = (-1)^k */
            t1 *= M_PI;/* t1 == pi^(k+1) */
            if(k >= 2)
              t2 *= k;/* t2 == k! == gamma(k+1) */
            if(j >= 0) /* by cheat above,  tt === d_k(x) */
              ans[j] = s*(ans[j] + t1/t2 * tt);
          }
          if (n == 0 && kode == 2)
            ans[0] += xln;
          return;
        } /* x <= 0 */

        *nz = 0;
        mm = m;
        nx = imin2(-Rf_i1mach(15), Rf_i1mach(16));/* = 1021 */
        r1m5 = Rf_d1mach(5);
        r1m4 = Rf_d1mach(4) * 0.5;
        wdtol = fmax2(r1m4, 0.5e-18); /* 1.11e-16 */

        /* elim = approximate exponential over and underflow limit */

        elim = 2.302 * (nx * r1m5 - 3.0);/* = 700.6174... */
        xln = log(x);
        for(;;) {
          nn = n + mm - 1;
          fn = nn;
          t = (fn + 1) * xln;

          /* overflow and underflow test for small and large x */

          if (fabs(t) > elim) {
            if (t <= 0.0) {
              *nz = 0;
              *ierr = 2;
              return;
            }
          }
          else {
            if (x < wdtol) {
              ans[0] = pow(x, -n-1.0);
              if (mm != 1) {
                for(k = 1; k < mm ; k++)
                  ans[k] = ans[k-1] / x;
              }
              if (n == 0 && kode == 2)
                ans[0] += xln;
              return;
            }

            /* compute xmin and the number of terms of the series,  fln+1 */

            rln = r1m5 * Rf_i1mach(14);
            rln = fmin2(rln, 18.06);
            fln = fmax2(rln, 3.0) - 3.0;
            yint = 3.50 + 0.40 * fln;
            slope = 0.21 + fln * (0.0006038 * fln + 0.008677);
            xm = yint + slope * fn;
            mx = (int)xm + 1;
            xmin = mx;
            if (n != 0) {
              xm = -2.302 * rln - fmin2(0.0, xln);
              arg = xm / n;
              arg = fmin2(0.0, arg);
              eps = exp(arg);
              xm = 1.0 - eps;
              if (fabs(arg) < 1.0e-3)
                xm = -arg;
              fln = x * xm / eps;
              xm = xmin - x;
              if (xm > 7.0 && fln < 15.0)
                break;
            }
            xdmy = x;
            xdmln = xln;
            xinc = 0.0;
            if (x < xmin) {
              nx = (int)x;
              xinc = xmin - nx;
              xdmy = x + xinc;
              xdmln = log(xdmy);
            }

            /* generate w(n+mm-1, x) by the asymptotic expansion */

            t = fn * xdmln;
            t1 = xdmln + xdmln;
            t2 = t + xdmln;
            tk = fmax2(fabs(t), fmax2(fabs(t1), fabs(t2)));
            if (tk <= elim)
              goto L10;
          }
          nz++;
          mm--;
          ans[mm] = 0.;
          if (mm == 0)
            return;
        }
        nn = (int)fln + 1;
        np = n + 1;
        t1 = (n + 1) * xln;
        t = exp(-t1);
        s = t;
        den = x;
        for(i=1; i <= nn; i++) {
          den += 1.;
          trm[i] = pow(den, (double)-np);
          s += trm[i];
        }
        ans[0] = s;
        if (n == 0 && kode == 2)
          ans[0] = s + xln;

        if (mm != 1) { /* generate higher derivatives, j > n */

          tol = wdtol / 5.0;
          for(j = 1; j < mm; j++) {
            t /= x;
            s = t;
            tols = t * tol;
            den = x;
            for(i=1; i <= nn; i++) {
              den += 1.;
              trm[i] /= den;
              s += trm[i];
              if (trm[i] < tols)
                break;
            }
            ans[j] = s;
          }
        }
        return;

      L10:
        tss = exp(-t);
        tt = 0.5 / xdmy;
        t1 = tt;
        tst = wdtol * tt;
        if (nn != 0)
          t1 = tt + 1.0 / fn;
        rxsq = 1.0 / (xdmy * xdmy);
        ta = 0.5 * rxsq;
        t = (fn + 1) * ta;
        s = t * bvalues[2];
        if (fabs(s) >= tst) {
          tk = 2.0;
          for(k = 4; k <= 22; k++) {
            t = t * ((tk + fn + 1)/(tk + 1.0))*((tk + fn)/(tk + 2.0)) * rxsq;
            trm[k] = t * bvalues[k-1];
            if (fabs(trm[k]) < tst)
              break;
            s += trm[k];
            tk += 2.;
          }
        }
        s = (s + t1) * tss;
        if (xinc != 0.0) {

          /* backward recur from xdmy to x */

          nx = (int)xinc;
          np = nn + 1;
          if (nx > n_max) {
            *nz = 0;
            *ierr = 3;
            return;
          }
          else {
            if (nn==0)
              goto L20;
            xm = xinc - 1.0;
            fx = x + xm;

            /* this loop should not be changed. fx is accurate when x is small */
            for(i = 1; i <= nx; i++) {
              trmr[i] = pow(fx, (double)-np);
              s += trmr[i];
              xm -= 1.;
              fx = x + xm;
            }
          }
        }
        ans[mm-1] = s;
        if (fn == 0)
          goto L30;

        /* generate lower derivatives,  j < n+mm-1 */

        for(j = 2; j <= mm; j++) {
          fn--;
          tss *= xdmy;
          t1 = tt;
          if (fn!=0)
            t1 = tt + 1.0 / fn;
          t = (fn + 1) * ta;
          s = t * bvalues[2];
          if (fabs(s) >= tst) {
            tk = 4 + fn;
            for(k=4; k <= 22; k++) {
              trm[k] = trm[k] * (fn + 1) / tk;
              if (fabs(trm[k]) < tst)
                break;
              s += trm[k];
              tk += 2.;
            }
          }
          s = (s + t1) * tss;
          if (xinc != 0.0) {
            if (fn == 0)
              goto L20;
            xm = xinc - 1.0;
            fx = x + xm;
            for(i=1 ; i<=nx ; i++) {
              trmr[i] = trmr[i] * fx;
              s += trmr[i];
              xm -= 1.;
              fx = x + xm;
            }
          }
          ans[mm - j] = s;
          if (fn == 0)
            goto L30;
        }
        return;

      L20:
        for(i = 1; i <= nx; i++)
          s += 1. / (x + nx - i);

      L30:
        if (kode!=2)
          ans[0] = s - xdmln;
        else if (xdmy != x) {
          xq = xdmy / x;
          ans[0] = s - log(xq);
        }
        return;
      } /* dpsifn() */

    double TriGamma(double x)
    {
      double ans;
      int nz, ierr;
      if(ISNAN(x)) return x;
      dpsifn(x, 1, 1, 1, &ans, &nz, &ierr);
      if(ierr != 0) 
      {
        errno = EDOM;
        return ML_NAN;
      }
      return ans;
    } // double triGamma()
// !!!!!!!!!!!!!!!!!!!!!!! end of trigamma routine !!!!!!!!!!!!!!!!!


private:
    int          mNgbWidth;
    int          mNgbHeight;
    int mFrac;
    SrcArithType mFactor;
    SrcArithType mGamma;
    SrcArithType mGammaNext;
    SrcArithType mMaxGamma;
    SrcArithType mSumDataPowerGamma;
    SrcArithType mSumDataLog;
    SrcArrayT* mData;
    SrcArrayT* mEcdf;
    bool mLast;
    SrcArithType mF;
    SrcArithType mFprime;
    SrcArithType mPrecision;
    int mIter;
    int mMaxNumIters;
    int mPhase;
    bool         mVerbose;
};

} // namespace Feature
} // namespace Core
} // namespace Impala

#endif

---