00001 #ifndef Impala_Core_Feature_IntWeibullNgbPnLoop_h
00002 #define Impala_Core_Feature_IntWeibullNgbPnLoop_h
00003
00004 #include "Core/Array/Pattern/Categories.h"
00005 #include "Core/Array/Pattern/Cnum.h"
00006 #include "Core/Array/Element/E1Cast.h"
00007 #include "Core/Array/Element/Vec3Real64.h"
00008 #include "Core/Array/Arrays.h"
00009
00010 #include "limits.h"
00011 #include "float.h"
00012 #include <errno.h>
00013
00014 namespace Impala
00015 {
00016 namespace Core
00017 {
00018 namespace Feature
00019 {
00020
00021
00022 #ifndef M_PI
00023 #define M_PI 3.1415926535897932384626433832795029
00024 #endif
00025 #ifndef MMGAMMA
00026 #define MMGAMMA 0.5772156649015328606065120900824024
00027 #endif
00028 #ifndef M_LN2
00029 #define M_LN2 0.6931471805599453094172321214581766
00030 #endif
00031
00032
00035 template<class DstArrayT, class SrcArrayT>
00036 class IntWeibullNgbPnLoop
00037 {
00038 public:
00039
00041 typedef Array::Pattern::TagLoop IteratorCategory;
00042
00044 typedef Array::Pattern::TagNPhase PhaseCategory;
00045
00046 typedef typename DstArrayT::StorType DstStorType;
00047 typedef typename DstArrayT::ArithType DstArithType;
00048 typedef typename SrcArrayT::StorType SrcStorType;
00049 typedef typename SrcArrayT::ArithType SrcArithType;
00050
00051
00053 IntWeibullNgbPnLoop(bool verbose, int ngbWidth, int ngbHeight)
00054 {
00055 mNgbWidth = ngbWidth;
00056 mNgbHeight = ngbHeight;
00057 mVerbose = verbose;
00058 mFrac = mNgbWidth*mNgbHeight;
00059 mGamma = 0.01;
00060 mGammaNext = 0.0;
00061 mLast = false;
00062 mData = new SrcArrayT(mFrac, 1, 0, 0);
00063 mEcdf = new SrcArrayT(mFrac, 1, 0, 0);
00064 }
00065
00067 ~IntWeibullNgbPnLoop()
00068 {
00069 delete mData;
00070 delete mEcdf;
00071 }
00072
00074 int
00075 Width()
00076 {
00077 return mNgbWidth;
00078 }
00079
00081 int
00082 Height()
00083 {
00084 return mNgbHeight;
00085 }
00086
00088 void
00089 Init(int phase, int x, int y, SrcArithType value)
00090 {
00091 if(mLast) return;
00092 int ph = phase%2;
00093
00094
00095
00096
00097 if (phase == 1)
00098 {
00099
00100
00101 mGamma = 0.01;
00102 mGammaNext = 0.0;
00103 mPrecision = 0.005;
00104 mMaxNumIters = 20;
00105 mIter = 1;
00106 mMaxGamma = 5;
00107 }
00108 mPhase = phase;
00109 if (ph)
00110 {
00111 mSumDataPowerGamma = 0.0;
00112 mSumDataLog = 0.0;
00113 mF = 0.0;
00114 mFprime = 0.0;
00115 }
00116
00117 }
00118
00120 void
00121 NextEl(int x, int y, SrcArithType value)
00122 {
00123 SrcArithType dataPowerGamma;
00124 SrcArithType sumItem;
00125 SrcArithType data = fabs(value) + 0.000001;
00126 dataPowerGamma = pow(data, mGamma);
00127 int phase = mPhase%2;
00128
00129 if(mVerbose)
00130 if (mGamma == 0.01 && phase == 0)
00131 std::cout << value << std::endl;
00132
00133
00134
00135 if(mLast)
00136 {
00137 mData->SetValue(fabs(value), x*mNgbWidth + y, 0);
00138 return;
00139 }
00140 switch (phase)
00141 {
00142 case 1 :
00143 mSumDataPowerGamma += dataPowerGamma;
00144 mSumDataLog += dataPowerGamma*log(data);
00145 break;
00146
00147 case 0 :
00148 sumItem = dataPowerGamma/mSumDataPowerGamma;
00149 mF += (sumItem)*log(mFrac*sumItem);
00150 mFprime += (log(dataPowerGamma/mSumDataPowerGamma) + 1)*
00151 (dataPowerGamma*(log(data)*mSumDataPowerGamma - mSumDataLog))/
00152 (mSumDataPowerGamma*mSumDataPowerGamma);
00153 break;
00154 }
00155
00156 }
00157
00159 void
00160 Done(int phase)
00161 {
00162 int ph = phase%2;
00163 if (mVerbose)
00164 std::cout << "done: phase = " << phase << std::endl;
00165 if(!ph)
00166 {
00167 if (mVerbose)
00168 std::cout << "mGamma = " << mGamma << std::endl;
00169
00170 mFprime = (1/pow(mGamma,2))*(1 - log(mGamma) - DiGamma(1/mGamma) -
00171 (1/mGamma)*TriGamma(1/mGamma) + mF) -(1/mGamma)*mFprime;
00172 mF = 1 + (1/mGamma)*(log(mGamma) + DiGamma(1/mGamma) - mF);
00173 mGammaNext = mGamma - mF/mFprime;
00174 if (mVerbose)
00175 std::cout << "mGammaNext = " << mGammaNext << std::endl;
00176 }
00177
00178 }
00179
00181 bool
00182 HasNextPhase(int lastPhase)
00183 {
00184 if(mLast)
00185 {
00186 mLast = false;
00187 return false;
00188 }
00189 int ph = lastPhase%2;
00190 if (mVerbose)
00191 std::cout << " WeibullNgbPnLoop::hasNext(" << lastPhase << ") "
00192 << std::endl;
00193
00194 if (ph)
00195 return true;
00196
00197 else
00198 {
00199 if (fabs(mGamma - mGammaNext) < mPrecision || mIter > mMaxNumIters ||
00200 mGammaNext > mMaxGamma || isnan(mGammaNext) || isinf(mGammaNext) || mGammaNext < 0)
00201 {
00202 if (mVerbose)
00203 {
00204 std::cout << "no next phase! " << "fabs(mGamma - mGammaNext) = " << fabs(mGamma - mGammaNext) << std::endl;
00205 std::cout << "mIter = " << mIter << " maxIter = " << mMaxNumIters << std::endl;
00206 std::cout << "maxGamma = " << mMaxGamma << std::endl;
00207 std::cout << "mGamma = " << mGamma << " mGammaNext = " << mGammaNext << std::endl;
00208 }
00209
00210 mIter = 1;
00211 mLast = true;
00212 return true;
00213 }
00214 else
00215 {
00216 if (mVerbose)
00217 std::cout << "there is next phase: iteration = " << mIter << std::endl;
00218 mGamma = mGammaNext;
00219 mIter++;
00220 return true;
00221 }
00222 }
00223 }
00224
00226 DstArithType
00227 Result()
00228 {
00229 double Beta = pow(mSumDataPowerGamma/mFrac, (1/mGamma));
00230 if(Beta == 0)
00231 Beta = 0.000001;
00232 Core::Array::Sort(mEcdf, mData);
00233 Core::Array::Set(mData, mEcdf);
00234 Core::Array::MulVal(mEcdf, mEcdf, 1/mData->Value(mFrac-1, 0));
00235 double AD = 0.0;
00236 double F;
00237 double dF;
00238 for (int i=0; i<mFrac; i++)
00239 {
00240 F = 1 - exp(-(1/mGamma)*pow(mData->Value(i, 0)/Beta, mGamma));
00241
00242 if(F !=0 && F != 1)
00243 {
00244 dF = (1/Beta)*pow( mData->Value(i, 0)/Beta, mGamma-1 )*F;
00245 AD += pow((mEcdf->Value(i, 0) - F), 2)*dF/(F*(1-F));
00246 }
00247 }
00248 AD = AD/mFrac;
00249 if (mVerbose)
00250 std::cout << "result: Gamma = " << mGamma << " Beta = " << Beta << " AD = " << AD << std::endl;
00251 return Core::Array::Element::Vec3Real64(Beta, mGamma, AD);
00252 }
00253
00254 double
00255 gammaFunc(double xx) const
00256 {
00257 double x,y,tmp,ser;
00258 static double cof[6] =
00259 {
00260 76.18009172947146, -86.50532032941677,
00261 24.01409824083091, -1.231739572450155,
00262 0.1208650973866179e-2, -0.5395239384953e-5
00263 };
00264
00265 int j;
00266 y=x=xx;
00267 tmp=x+5.5;
00268 tmp -= (x+0.5)*log(tmp);
00269 ser=1.000000000190015;
00270 for (j=0;j<=5;j++)
00271 ser += cof[j]/++y;
00272 xx = -tmp+log(2.5066282746310005*ser/x);
00273 return exp(xx);
00274 }
00275
00276 double
00277 DiGamma(double x) const
00278 {
00279
00280 if( x < 0.0 )
00281 return DiGamma(1.0-x)+M_PI/tan(M_PI*(1.0-x)) ;
00282 else if( x < 1.0 )
00283 return DiGamma(1.0+x)-1.0/x ;
00284 else if ( x == 1.0)
00285 return -MMGAMMA ;
00286 else if ( x == 2.0)
00287 return 1.0-MMGAMMA ;
00288 else if ( x == 3.0)
00289 return 1.5-MMGAMMA ;
00290 else if ( x > 3.0)
00291
00292 return 0.5*(DiGamma(x/2.0)+DiGamma((x+1.0)/2.0))+M_LN2 ;
00293 else
00294 {
00295 static double Kncoe[] =
00296 { .30459198558715155634315638246624251,
00297 .72037977439182833573548891941219706, -.12454959243861367729528855995001087,
00298 .27769457331927827002810119567456810e-1, -.67762371439822456447373550186163070e-2,
00299 .17238755142247705209823876688592170e-2, -.44817699064252933515310345718960928e-3,
00300 .11793660000155572716272710617753373e-3, -.31253894280980134452125172274246963e-4,
00301 .83173997012173283398932708991137488e-5, -.22191427643780045431149221890172210e-5,
00302 .59302266729329346291029599913617915e-6, -.15863051191470655433559920279603632e-6,
00303 .42459203983193603241777510648681429e-7, -.11369129616951114238848106591780146e-7,
00304 .304502217295931698401459168423403510e-8, -.81568455080753152802915013641723686e-9,
00305 .21852324749975455125936715817306383e-9, -.58546491441689515680751900276454407e-10,
00306 .15686348450871204869813586459513648e-10, -.42029496273143231373796179302482033e-11,
00307 .11261435719264907097227520956710754e-11, -.30174353636860279765375177200637590e-12,
00308 .80850955256389526647406571868193768e-13, -.21663779809421233144009565199997351e-13,
00309 .58047634271339391495076374966835526e-14, -.15553767189204733561108869588173845e-14,
00310 .41676108598040807753707828039353330e-15, -.11167065064221317094734023242188463e-15 } ;
00311
00312 double Tn_1 = 1.0 ;
00313 double Tn = x-2.0 ;
00314 double resul = Kncoe[0] + Kncoe[1]*Tn ;
00315
00316 x -= 2.0 ;
00317
00318 for(int n = 2 ; n < sizeof(Kncoe)/sizeof(double) ;n++)
00319 {
00320 const double Tn1 = 2.0 * x * Tn - Tn_1 ;
00321 resul += Kncoe[n]*Tn1 ;
00322 Tn_1 = Tn ;
00323 Tn = Tn1 ;
00324 }
00325 return resul ;
00326 }
00327 }
00328
00331 #define ISNAN(x) (isnan(x) !=0)
00332 #define ML_POSINF (1.0 / 0.0)
00333 #define ML_NEGINF ((-1.0) / 0.0)
00334 #define ML_NAN (0.0 / 0.0)
00335 #define n_max (100)
00336 #define M_LOG10_2 0.30102999566398119521
00337
00338 int Rf_i1mach(int i)
00339 {
00340 switch(i)
00341 {
00342
00343 case 1: return 5;
00344 case 2: return 6;
00345 case 3: return 0;
00346 case 4: return 0;
00347
00348 case 5: return CHAR_BIT * sizeof(int);
00349 case 6: return sizeof(int)/sizeof(char);
00350
00351 case 7: return 2;
00352 case 8: return CHAR_BIT * sizeof(int) - 1;
00353 case 9: return INT_MAX;
00354
00355 case 10: return FLT_RADIX;
00356
00357 case 11: return FLT_MANT_DIG;
00358 case 12: return FLT_MIN_EXP;
00359
00360 case 13: return FLT_MAX_EXP;
00361
00362 case 14: return DBL_MANT_DIG;
00363 case 15: return DBL_MIN_EXP;
00364 case 16: return DBL_MAX_EXP;
00365
00366 default: return 0;
00367 }
00368 }
00369
00370 double Rf_d1mach(int i)
00371 {
00372 switch(i)
00373 {
00374
00375 case 1: return DBL_MIN;
00376 case 2: return DBL_MAX;
00377
00378 case 3:
00379
00380 return 0.5*DBL_EPSILON;
00381
00382 case 4:
00383
00384 return DBL_EPSILON;
00385
00386 case 5: return M_LOG10_2;
00387
00388 default: return 0.0;
00389 }
00390 }
00391
00392 double fmax2(double x, double y)
00393 {
00394 #ifdef IEEE_754
00395 if (ISNAN(x) || ISNAN(y))
00396 return x + y;
00397 #endif
00398 return (x < y) ? y : x;
00399 }
00400
00401 double fmin2(double x, double y)
00402 {
00403 #ifdef IEEE_754
00404 if (ISNAN(x) || ISNAN(y))
00405 return x + y;
00406 #endif
00407 return (x < y) ? x : y;
00408 }
00409
00410 int imin2(int x, int y)
00411 {
00412 return (x < y) ? x : y;
00413 }
00414
00415 void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr)
00416 {
00417 const static double bvalues[] = {
00418 1.00000000000000000e+00,
00419 -5.00000000000000000e-01,
00420 1.66666666666666667e-01,
00421 -3.33333333333333333e-02,
00422 2.38095238095238095e-02,
00423 -3.33333333333333333e-02,
00424 7.57575757575757576e-02,
00425 -2.53113553113553114e-01,
00426 1.16666666666666667e+00,
00427 -7.09215686274509804e+00,
00428 5.49711779448621554e+01,
00429 -5.29124242424242424e+02,
00430 6.19212318840579710e+03,
00431 -8.65802531135531136e+04,
00432 1.42551716666666667e+06,
00433 -2.72982310678160920e+07,
00434 6.01580873900642368e+08,
00435 -1.51163157670921569e+10,
00436 4.29614643061166667e+11,
00437 -1.37116552050883328e+13,
00438 4.88332318973593167e+14,
00439 -1.92965793419400681e+16
00440 };
00441
00442 int i, j, k, mm, mx, nn, np, nx, fn;
00443 double arg, den, elim, eps, fln, fx, rln, rxsq,
00444 r1m4, r1m5, s, slope, t, ta, tk, tol, tols, tss, tst,
00445 tt, t1, t2, wdtol, xdmln, xdmy, xinc, xln = 0.0 ,
00446 xm, xmin, xq, yint;
00447 double trm[23], trmr[n_max + 1];
00448
00449 *ierr = 0;
00450 if (n < 0 || kode < 1 || kode > 2 || m < 1) {
00451 *ierr = 1;
00452 return;
00453 }
00454 if (x <= 0.) {
00455
00456
00457
00458 if (x == (long)x) {
00459
00460 for(j=0; j < m; j++)
00461 ans[j] = ((j+n) % 2) ? ML_POSINF : ML_NAN;
00462 return;
00463 }
00464 dpsifn(1. - x, n, 1, m, ans, nz, ierr);
00465
00466
00467
00468
00469
00470 if(m > 1 || n > 3) {
00471
00472 *ierr = 4; return;
00473 }
00474 x *= M_PI;
00475 if (n == 0)
00476 tt = cos(x)/sin(x);
00477 else if (n == 1)
00478 tt = -1/pow(sin(x),2);
00479 else if (n == 2)
00480 tt = 2*cos(x)/pow(sin(x),3);
00481 else if (n == 3)
00482 tt = -2*(2*pow(cos(x),2) + 1)/pow(sin(x),4);
00483 else
00484 tt = ML_NAN;
00485
00486
00487 s = (n % 2) ? -1. : 1.;
00488
00489
00490 t1 = t2 = s = 1.;
00491 for(k=0, j=k-n; j < m; k++, j++, s = -s) {
00492
00493 t1 *= M_PI;
00494 if(k >= 2)
00495 t2 *= k;
00496 if(j >= 0)
00497 ans[j] = s*(ans[j] + t1/t2 * tt);
00498 }
00499 if (n == 0 && kode == 2)
00500 ans[0] += xln;
00501 return;
00502 }
00503
00504 *nz = 0;
00505 mm = m;
00506 nx = imin2(-Rf_i1mach(15), Rf_i1mach(16));
00507 r1m5 = Rf_d1mach(5);
00508 r1m4 = Rf_d1mach(4) * 0.5;
00509 wdtol = fmax2(r1m4, 0.5e-18);
00510
00511
00512
00513 elim = 2.302 * (nx * r1m5 - 3.0);
00514 xln = log(x);
00515 for(;;) {
00516 nn = n + mm - 1;
00517 fn = nn;
00518 t = (fn + 1) * xln;
00519
00520
00521
00522 if (fabs(t) > elim) {
00523 if (t <= 0.0) {
00524 *nz = 0;
00525 *ierr = 2;
00526 return;
00527 }
00528 }
00529 else {
00530 if (x < wdtol) {
00531 ans[0] = pow(x, -n-1.0);
00532 if (mm != 1) {
00533 for(k = 1; k < mm ; k++)
00534 ans[k] = ans[k-1] / x;
00535 }
00536 if (n == 0 && kode == 2)
00537 ans[0] += xln;
00538 return;
00539 }
00540
00541
00542
00543 rln = r1m5 * Rf_i1mach(14);
00544 rln = fmin2(rln, 18.06);
00545 fln = fmax2(rln, 3.0) - 3.0;
00546 yint = 3.50 + 0.40 * fln;
00547 slope = 0.21 + fln * (0.0006038 * fln + 0.008677);
00548 xm = yint + slope * fn;
00549 mx = (int)xm + 1;
00550 xmin = mx;
00551 if (n != 0) {
00552 xm = -2.302 * rln - fmin2(0.0, xln);
00553 arg = xm / n;
00554 arg = fmin2(0.0, arg);
00555 eps = exp(arg);
00556 xm = 1.0 - eps;
00557 if (fabs(arg) < 1.0e-3)
00558 xm = -arg;
00559 fln = x * xm / eps;
00560 xm = xmin - x;
00561 if (xm > 7.0 && fln < 15.0)
00562 break;
00563 }
00564 xdmy = x;
00565 xdmln = xln;
00566 xinc = 0.0;
00567 if (x < xmin) {
00568 nx = (int)x;
00569 xinc = xmin - nx;
00570 xdmy = x + xinc;
00571 xdmln = log(xdmy);
00572 }
00573
00574
00575
00576 t = fn * xdmln;
00577 t1 = xdmln + xdmln;
00578 t2 = t + xdmln;
00579 tk = fmax2(fabs(t), fmax2(fabs(t1), fabs(t2)));
00580 if (tk <= elim)
00581 goto L10;
00582 }
00583 nz++;
00584 mm--;
00585 ans[mm] = 0.;
00586 if (mm == 0)
00587 return;
00588 }
00589 nn = (int)fln + 1;
00590 np = n + 1;
00591 t1 = (n + 1) * xln;
00592 t = exp(-t1);
00593 s = t;
00594 den = x;
00595 for(i=1; i <= nn; i++) {
00596 den += 1.;
00597 trm[i] = pow(den, (double)-np);
00598 s += trm[i];
00599 }
00600 ans[0] = s;
00601 if (n == 0 && kode == 2)
00602 ans[0] = s + xln;
00603
00604 if (mm != 1) {
00605
00606 tol = wdtol / 5.0;
00607 for(j = 1; j < mm; j++) {
00608 t /= x;
00609 s = t;
00610 tols = t * tol;
00611 den = x;
00612 for(i=1; i <= nn; i++) {
00613 den += 1.;
00614 trm[i] /= den;
00615 s += trm[i];
00616 if (trm[i] < tols)
00617 break;
00618 }
00619 ans[j] = s;
00620 }
00621 }
00622 return;
00623
00624 L10:
00625 tss = exp(-t);
00626 tt = 0.5 / xdmy;
00627 t1 = tt;
00628 tst = wdtol * tt;
00629 if (nn != 0)
00630 t1 = tt + 1.0 / fn;
00631 rxsq = 1.0 / (xdmy * xdmy);
00632 ta = 0.5 * rxsq;
00633 t = (fn + 1) * ta;
00634 s = t * bvalues[2];
00635 if (fabs(s) >= tst) {
00636 tk = 2.0;
00637 for(k = 4; k <= 22; k++) {
00638 t = t * ((tk + fn + 1)/(tk + 1.0))*((tk + fn)/(tk + 2.0)) * rxsq;
00639 trm[k] = t * bvalues[k-1];
00640 if (fabs(trm[k]) < tst)
00641 break;
00642 s += trm[k];
00643 tk += 2.;
00644 }
00645 }
00646 s = (s + t1) * tss;
00647 if (xinc != 0.0) {
00648
00649
00650
00651 nx = (int)xinc;
00652 np = nn + 1;
00653 if (nx > n_max) {
00654 *nz = 0;
00655 *ierr = 3;
00656 return;
00657 }
00658 else {
00659 if (nn==0)
00660 goto L20;
00661 xm = xinc - 1.0;
00662 fx = x + xm;
00663
00664
00665 for(i = 1; i <= nx; i++) {
00666 trmr[i] = pow(fx, (double)-np);
00667 s += trmr[i];
00668 xm -= 1.;
00669 fx = x + xm;
00670 }
00671 }
00672 }
00673 ans[mm-1] = s;
00674 if (fn == 0)
00675 goto L30;
00676
00677
00678
00679 for(j = 2; j <= mm; j++) {
00680 fn--;
00681 tss *= xdmy;
00682 t1 = tt;
00683 if (fn!=0)
00684 t1 = tt + 1.0 / fn;
00685 t = (fn + 1) * ta;
00686 s = t * bvalues[2];
00687 if (fabs(s) >= tst) {
00688 tk = 4 + fn;
00689 for(k=4; k <= 22; k++) {
00690 trm[k] = trm[k] * (fn + 1) / tk;
00691 if (fabs(trm[k]) < tst)
00692 break;
00693 s += trm[k];
00694 tk += 2.;
00695 }
00696 }
00697 s = (s + t1) * tss;
00698 if (xinc != 0.0) {
00699 if (fn == 0)
00700 goto L20;
00701 xm = xinc - 1.0;
00702 fx = x + xm;
00703 for(i=1 ; i<=nx ; i++) {
00704 trmr[i] = trmr[i] * fx;
00705 s += trmr[i];
00706 xm -= 1.;
00707 fx = x + xm;
00708 }
00709 }
00710 ans[mm - j] = s;
00711 if (fn == 0)
00712 goto L30;
00713 }
00714 return;
00715
00716 L20:
00717 for(i = 1; i <= nx; i++)
00718 s += 1. / (x + nx - i);
00719
00720 L30:
00721 if (kode!=2)
00722 ans[0] = s - xdmln;
00723 else if (xdmy != x) {
00724 xq = xdmy / x;
00725 ans[0] = s - log(xq);
00726 }
00727 return;
00728 }
00729
00730 double TriGamma(double x)
00731 {
00732 double ans;
00733 int nz, ierr;
00734 if(ISNAN(x)) return x;
00735 dpsifn(x, 1, 1, 1, &ans, &nz, &ierr);
00736 if(ierr != 0)
00737 {
00738 errno = EDOM;
00739 return ML_NAN;
00740 }
00741 return ans;
00742 }
00743
00744
00745
00746 private:
00747 int mNgbWidth;
00748 int mNgbHeight;
00749 int mFrac;
00750 SrcArithType mFactor;
00751 SrcArithType mGamma;
00752 SrcArithType mGammaNext;
00753 SrcArithType mMaxGamma;
00754 SrcArithType mSumDataPowerGamma;
00755 SrcArithType mSumDataLog;
00756 SrcArrayT* mData;
00757 SrcArrayT* mEcdf;
00758 bool mLast;
00759 SrcArithType mF;
00760 SrcArithType mFprime;
00761 SrcArithType mPrecision;
00762 int mIter;
00763 int mMaxNumIters;
00764 int mPhase;
00765 bool mVerbose;
00766 };
00767
00768 }
00769 }
00770 }
00771
00772 #endif