CongealGroupPhased.h

Go to the documentation of this file.
00001 #ifndef __CONGEALGROUPPHASED_H__
00002 #define __CONGEALGROUPPHASED_H__
00003 
00004 #include "libOptimize_NelderMead.h"
00005 #include "libOptimize_GreedySwarm.h"
00006 #include "libOptimize_RandomWalk.h"
00007 #include "libErrorFunctions.h"
00008 #include <string.h> //for memcpy()
00009 #include "Pixels.h"
00010 #include "libBaseTypes.h"
00011 #include "Registry.h"
00012 #include "libConfiguration.h"
00013 #include "Recipie.h"
00014 
00015 #undef DEBUG
00016 #define DEBUG 0
00017 #include "libDebug.h"
00018 
00019 #define MULTITHREADED 1
00020 
00021 #if MULTITHREADED
00022 #  include <pthread.h>
00023 #  include "libBarrier.h"
00024 #  include <unistd.h>
00025 #endif
00026 
00027 #ifndef FEATURE_LBFGS
00028 #  define FEATURE_LBFGS 0
00029 #endif
00030 
00031 #define LBFGS 0
00032 #define GRADIENTDESCENT 2
00033 #define NELDERMEAD 3
00034 
00035 #if FEATURE_LBFGS
00036 #  include <lbfgs.h>
00037 #endif
00038 
00039 class NVar
00040 {
00041 public:
00042   int m_n;
00043   Real m_rVar;
00044   static int SortByVar(const void* pnvar1, const void *pnvar2)
00045   {
00046     const NVar &nvar1 = *((NVar*) pnvar1);
00047     const NVar &nvar2 = *((NVar*) pnvar2);
00048 
00049     if (nvar1.m_rVar < nvar2.m_rVar)
00050       {
00051       return -1;
00052       }
00053     return 1;
00054   }
00055 };
00056 
00057 //============================================================================
00058 //=================================================p===========================
00059 template<class RECIPIE>
00060 class CongealGroupPhasedOf
00061 {
00062   static long int m_cLookups;
00063 
00064   typedef CongealGroupPhasedOf<RECIPIE> THIS_TYPE;
00065   static const int DIMENSIONALITY = RECIPIE::type_dimensionality;
00066   typedef typename RECIPIE::type_precision PRECISION;
00067   typedef typename RECIPIE::type_data DATA;
00068   typedef PointOf<DIMENSIONALITY, PRECISION> POINT;
00069 
00070   RECIPIE* m_precipie;
00071 
00072   RegistryOfParameters* m_vregParams;
00073   RegistryOfInitialSteps m_regSteps;
00074   int m_cSamples;
00075   int m_cSources;RANGEPROMOTION(DATA)* m_vdataSum;
00076   RANGEPROMOTION(DATA)* m_vdataSumSq;
00077   Real m_rStep;
00078 
00079   DATA* m_mdata;
00080   POINT* m_vpt;
00081 
00082   Parameter* m_vparam;
00083 
00084 #if MULTITHREADED
00085   bool m_bDie;
00086   pthread_t* m_vthread;
00087   pthread_mutex_t m_mutexReady;
00088   pthread_cond_t m_condReady;
00089   pthread_barrier_t m_barrierDone;
00090 #endif
00091 
00092   struct AlignInfo // ai
00093   {
00094     bool m_bReady;
00095     int m_nSource;
00096     Real m_rErr;
00097     Real* m_vparam;
00098     CongealGroupPhasedOf<RECIPIE>* m_pcongeal;
00099   };
00100 
00101   AlignInfo* m_vai;
00102 
00103 public:
00104   //--------------------------------------------------------------------------
00105   //--------------------------------------------------------------------------
00106   CongealGroupPhasedOf(RECIPIE* precipie)
00107   {
00108     static Real rSigma = ConfigValue("congeal.error[parzen].sigma", 10.);
00109     InitFastGaussian(rSigma);
00110 
00111     m_precipie = precipie;
00112     ClaimPointer(m_precipie);
00113 
00114     m_cSources = m_precipie->CSources();
00115     m_vregParams = new RegistryOfParameters[m_cSources];
00116 
00117     m_cSamples = -1;
00118     m_mdata = NULL;
00119     m_vpt = NULL;
00120     m_vdataSum = NULL;
00121     m_vdataSumSq = NULL;
00122 
00123     // create multiple AlignInfo's so we can eventually launch
00124     // threads for each source
00125     m_vai = new AlignInfo[m_cSources];
00126     for (int nSource = 0; nSource < m_cSources; nSource++)
00127       {
00128       m_vai[nSource].m_bReady = false;
00129       m_vai[nSource].m_nSource = nSource;
00130       m_vai[nSource].m_rErr = 0;
00131       m_vai[nSource].m_pcongeal = this;
00132       m_vai[nSource].m_vparam = NULL;
00133       }
00134 
00135 #if MULTITHREADED
00136     // set up multithreading sync devices
00137     pthread_mutex_init(&m_mutexReady, NULL);
00138     pthread_cond_init(&m_condReady, NULL);
00139     pthread_barrier_init(&m_barrierDone, NULL, m_cSources + 1);
00140 
00141     m_bDie = false;
00142     // launch threads
00143     m_vthread = new pthread_t[m_cSources];
00144     for (int nSource = 0; nSource < m_cSources; nSource++)
00145       {
00146       pthread_create(&(m_vthread[nSource]), NULL, OptimizeSource,
00147           &m_vai[nSource]);
00148       }
00149 #endif
00150   }
00151 
00152   //--------------------------------------------------------------------------
00153   //--------------------------------------------------------------------------
00154   ~CongealGroupPhasedOf()
00155   {
00156 #if MULTITHREADED
00157     pthread_mutex_lock(&m_mutexReady);
00158     m_bDie = true;
00159     pthread_cond_broadcast(&m_condReady);
00160     pthread_mutex_unlock(&m_mutexReady);
00161     pthread_barrier_wait(&m_barrierDone);
00162 
00163     pthread_mutex_destroy(&m_mutexReady);
00164     pthread_cond_destroy(&m_condReady);
00165     pthread_barrier_destroy(&m_barrierDone);
00166     for (int nSource = 0; nSource < m_cSources; nSource++)
00167       {
00168       pthread_detach(m_vthread[nSource]);
00169       }
00170     delete[] m_vthread;
00171 #endif
00172 
00173     delete[] m_vregParams;
00174 
00175     for (int nSource = 0; nSource < m_cSources; nSource++)
00176       {
00177       SafeDeleteArray(m_vai[nSource].m_vparam);
00178       }
00179     delete[] m_vai;
00180 
00181     SafeDeleteArray(m_mdata);
00182     SafeDeleteArray(m_vpt);
00183     SafeDeleteArray( m_vdataSum);
00184     SafeDeleteArray( m_vdataSumSq);
00185 
00186     ReleasePointer(m_precipie);
00187   }
00188 
00189 private:
00190   void AllocateSampleSpace(int cSamples)
00191   {
00192     if (m_cSamples != cSamples)
00193       {
00194       m_cSamples = cSamples;
00195 
00196       SafeDeleteArray(m_mdata);
00197       SafeDeleteArray(m_vpt);
00198       SafeDeleteArray( m_vdataSum);
00199       SafeDeleteArray( m_vdataSumSq);
00200 
00201       m_mdata = new DATA[m_cSamples * m_cSources];
00202       m_vpt = new POINT[m_cSamples];
00203 m_vdataSum      =new RANGEPROMOTION(DATA)[m_cSamples];
00204       m_vdataSumSq=new RANGEPROMOTION(DATA)[m_cSamples];
00205       }
00206     }
00207 
00208   void SetUpRegistries()
00209     {
00210     // register parameters and steps
00211     // NOTE: this assumes all recipies have the same number,
00212     // same ordering and same types of of parameters
00213     m_regSteps.Clear();
00214     m_precipie->RegisterInitialSteps(0,m_regSteps);
00215     for (int nSource=0; nSource<m_cSources; nSource++)
00216       {
00217       m_vregParams[nSource].Clear();
00218       m_precipie->RegisterParameters(nSource,m_vregParams[nSource]);
00219       ASSERTf(m_vregParams[nSource].C() == m_regSteps.C(),"%d != %d",m_vregParams[nSource].C(), m_regSteps.C());
00220       }
00221 
00222     // allocate scratch transfer space for using generalize optimization
00223     // methods which require a packed array of parameters
00224     int cSources=m_cSources;
00225     for (int nSource=0; nSource<cSources; nSource++)
00226       {
00227       SafeDeleteArray(m_vai[nSource].m_vparam);
00228       m_vai[nSource].m_vparam = new Parameter[m_regSteps.C()];
00229       }
00230     }
00231 
00232 public:
00233 
00234   //--------------------------------------------------------------------------
00235   //--------------------------------------------------------------------------
00236   static int SortPoints(const void* ppt1 ,const void *ppt2)
00237     {
00238     const POINT &pt1 = *((POINT*)ppt1);
00239     const POINT &pt2 = *((POINT*)ppt2);
00240 
00241     if (pt1.Z()<pt2.Z())
00242       {
00243       return -1;
00244       }
00245     if (pt1.Z()>pt2.Z())
00246       {
00247       return 1;
00248       }
00249 
00250     if (pt1.Y()<pt2.Y())
00251       {
00252       return -1;
00253       }
00254     if (pt1.Y()>pt2.Y())
00255       {
00256       return 1;
00257       }
00258 
00259     if (pt1.X()<pt2.X())
00260       {
00261       return -1;
00262       }
00263     if (pt1.X()>pt2.X())
00264       {
00265       return 1;
00266       }
00267 
00268     return 0;
00269     }
00270 
00271   //--------------------------------------------------------------------------
00272   //--------------------------------------------------------------------------
00273   void StatisticalPhasedCongeal(
00274       int cSamples,
00275       Real rThreshold=.001,
00276       int cIterations=1000,
00277       int nFirstIteration=0,
00278       int nLastIteration=-1
00279   )
00280     {
00281     if (nLastIteration<nFirstIteration)
00282       {
00283       nLastIteration=cIterations-1;
00284       }
00285 
00286     AllocateSampleSpace(cSamples);
00287     SetUpRegistries();
00288 
00289     m_precipie->PrepareForAccess();
00290     UD("OPTIMIZING %d PARAMETERS PER VOLUME", m_regSteps.C());
00291     UD("WITH %d MEASUREMENTS (%2.4g%%)", cSamples, Real(100*cSamples)/m_precipie->CSize());
00292 
00293     Real rErr=REAL_MAX;
00294     for (int nIteration=nFirstIteration; nIteration<=nLastIteration && rErr > rThreshold; nIteration++)
00295       {
00296       m_cLookups=0;
00297 
00298       m_rStep=(Real(cIterations-nIteration)/cIterations);
00299 
00300       UD( "---------------- ITERATION %d ---------------- %g, %g", nIteration,rErr,m_rStep);
00301 
00302       if (1||nIteration==0)
00303         {
00304         // compute a random collection of samples, retrying if the point is not within a good pixel stack
00305         m_precipie->PrepareForAccess();
00306         for (int nSample=0; nSample<cSamples; nSample++)
00307           {
00308           bool bGoodPoint;
00309 
00310           int cTries=0;
00311           do
00312             {
00313             bGoodPoint=true;
00314             m_vpt[nSample]=NewRandomPoint(m_precipie->Size());
00315             DATA *vdataGroup=&m_mdata[nSample*m_cSources];
00316             for (int nSource=0; nSource<m_cSources; nSource++)
00317               {
00318               m_precipie->PSource(nSource)->GetPoint(
00319                   vdataGroup[nSource],m_vpt[nSample]
00320               );
00321               cTries++;
00322               m_cLookups++;
00323               if (vdataGroup[nSource]==0)
00324                 {
00325                 bGoodPoint=false;
00326                 //                              UD(String("point failed for source ") + nSource + " :: " +m_vpt[nSample]);
00327 
00328                 break;
00329                 }
00330               }
00331             }
00332           while(!bGoodPoint);
00333           //                    UD("Got point %d of %d in %d tries", nSample, cSamples,cTries);
00334           }
00335 
00336         //              qsort(m_vpt,cSamples,sizeof(POINT),CongealGroupPhasedOf<RECIPIE>::SortPoints);
00337 
00338         // set up group statistics
00339         for (int nSample=0; nSample<cSamples; nSample++)
00340           {
00341           DATA *vdataGroup=&m_mdata[nSample*m_cSources];
00342           m_vdataSum[nSample]=0;
00343           m_vdataSumSq[nSample]=0;
00344           for (int nSource=0; nSource<m_cSources; nSource++)
00345             {
00346             m_vdataSum[nSample] += vdataGroup[nSource];
00347             m_vdataSumSq[nSample] += sqr(vdataGroup[nSource]);
00348             }
00349           }
00350 
00351         NVar* vnvar=new NVar[cSamples];
00352         // compute most unique points
00353         for (int nSample=0; nSample<cSamples; nSample++)
00354           {
00355           vnvar[nSample].m_n=nSample;
00356           vnvar[nSample].m_rVar=m_vdataSumSq[nSample]/cSamples - sqr(m_vdataSum[nSample]/cSamples);
00357           }
00358         qsort(vnvar,cSamples,sizeof(NVar),NVar::SortByVar);
00359 
00360         int cPoints=ConfigValue("congeal.optimize.bestpoints",100);
00361         POINT* vpt=new POINT[cPoints];
00362         for (int nPoint=0; nPoint<cPoints; nPoint++)
00363           {
00364           vpt[nPoint]=m_vpt[vnvar[nPoint].m_n];
00365           }
00366 
00367         cSamples=cPoints;
00368 
00369         AllocateSampleSpace(cSamples);
00370         SetUpRegistries();
00371 
00372         m_precipie->PrepareForAccess();
00373 
00374         for (int nPoint=0; nPoint<cPoints; nPoint++)
00375           {
00376           m_vpt[nPoint]=vpt[nPoint];
00377           //                    D(String("POINT ") + nPoint + " = "  + m_vpt[nPoint]);
00378 
00379           }
00380         delete [] vpt;
00381         delete [] vnvar;
00382 
00383         qsort(m_vpt,cSamples,sizeof(POINT),CongealGroupPhasedOf<RECIPIE>::SortPoints);
00384 
00385         for (int nSample=0; nSample<cSamples; nSample++)
00386           {
00387           DATA *vdataGroup=&m_mdata[nSample*m_cSources];
00388           for (int nSource=0; nSource<m_cSources; nSource++)
00389             {
00390             m_precipie->PSource(nSource)->GetPoint(
00391                 vdataGroup[nSource],m_vpt[nSample]
00392             );
00393             m_cLookups++;
00394             //                      D(String("POINT ") + nSample + "/" + nSource + " = "  + m_vpt[nSample] + " ==> " + vdataGroup[nSource]);
00395             }
00396           }
00397         }
00398 
00399       // optimize each source against the collection
00400       TICK(dtmOptimize);
00401 #if MULTITHREADED
00402       pthread_mutex_lock(&m_mutexReady);
00403       for (int nSource=0; nSource<m_cSources; nSource++)
00404         {
00405         m_vai[nSource].m_bReady=true;
00406         }
00407       pthread_cond_broadcast(&m_condReady);
00408       pthread_mutex_unlock(&m_mutexReady);
00409 
00410       pthread_barrier_wait(&m_barrierDone);
00411 #else 
00412       // sequentially
00413       for (int nSource=0; nSource<m_cSources; nSource++)
00414         {
00415         D("optimizing source %d", nSource);
00416         OptimizeSource(&m_vai[nSource]);
00417         }
00418 #endif
00419       UDTOCK("Time to optimize group %g", dtmOptimize);
00420       UD("TOTAL LOOKUPS %ld (%g lookups/s)", m_cLookups, (Real)m_cLookups/TOCK(dtmOptimize));
00421 
00422       // Normalize the data
00423       m_precipie->Normalize();
00424 
00425       rErr=0;
00426       for (int nSource=0; nSource<m_cSources; nSource++)
00427         {
00428         rErr+=m_vai[nSource].m_rErr;
00429         }
00430       rErr/=m_cSources;
00431       }
00432     }
00433 
00434   //--------------------------------------------------------------------------
00435   //--------------------------------------------------------------------------
00436   static void* OptimizeSource(void*p)
00437     {
00438     AlignInfo* pai = (AlignInfo*)p;
00439     pai->m_pcongeal->OptimizeSource(pai);
00440     return NULL;
00441     }
00442 
00443   //--------------------------------------------------------------------------
00444   //--------------------------------------------------------------------------
00445   void OptimizeSource(AlignInfo *pai)
00446     {
00447     int nSource=pai->m_nSource;
00448     String sAlgorithm = ConfigValue("congeal.optimize.algorithm","gradientdescent");
00449     String sError = ConfigValue("congeal.optimize.error","parzen");
00450 
00451     while(1)
00452       {
00453 
00454 #if MULTITHREADED
00455       pthread_mutex_lock(&m_mutexReady);
00456       while(!pai->m_bReady && !m_bDie)
00457         {
00458         pthread_cond_wait(&m_condReady,&m_mutexReady);
00459         }
00460       pai->m_bReady=false;
00461       pthread_mutex_unlock(&m_mutexReady);
00462 
00463       if (m_bDie)
00464         {
00465         pthread_barrier_wait(&m_barrierDone);
00466         return;
00467         }
00468 #endif
00469 
00470       // gradientdescent -------------------------------------
00471       if (sAlgorithm == "gradientdescent")
00472         {
00473         pai->m_rErr = Optimize_OneGradientStep
00474         (
00475             m_vregParams[nSource].C(),
00476             m_vregParams[nSource].V(),
00477             m_regSteps.V(),
00478             ( (sError == "parzen")
00479                 ? MeasureInverseParzenProbabilityInPlace
00480                 : MeasureVarianceInPlace),
00481             pai,
00482             m_rStep
00483         );
00484         }
00485 
00486       // gradientdescent -------------------------------------
00487 
00488       else if (sAlgorithm == "randomwalk")
00489         {
00490         pai->m_rErr = Optimize_GreedyRandomWalk
00491         (
00492             m_vregParams[nSource].C(),
00493             m_vregParams[nSource].V(),
00494             m_regSteps.V(),
00495             ( (sError == "parzen")
00496                 ? MeasureInverseParzenProbabilityInPlace
00497                 : MeasureVarianceInPlace),
00498             pai,
00499             m_rStep
00500         );
00501         }
00502 
00503       // brute force -------------------------------------
00504 
00505       else if(sAlgorithm == "bruteforce")
00506         {
00507         pai->m_rErr = Optimize_BruteForce
00508         (
00509             m_vregParams[nSource].C(),
00510             m_vregParams[nSource].V(),
00511             m_regSteps.V(),
00512             MeasureInverseParzenProbabilityInPlace,
00513             pai,
00514             m_rStep
00515         );
00516         }
00517 #if FEATURE_LBFGS
00518       // L-BFGS -------------------------------------------
00519 
00520       else if (sAlgorithm == "lbfgs")
00521         {
00522         ParamsToArray(pai->m_vparam, m_vregParams[nSource]);
00523 
00524         lbfgs(
00525             m_vregParams[nSource].C(),
00526             pai->m_vparam,
00527             &pai->m_rErr,
00528             EvaluateForLBFGS,
00529             NULL,
00530             pai,
00531             NULL
00532         );
00533 
00534         ParamsFromArray(m_vregParams[nSource],pai->m_vparam);
00535         }
00536 #endif
00537       // Nelder-Mead -------------------------------------------
00538 
00539       else if (sAlgorithm == "neldermead" || sAlgorithm == "simplex")
00540         {
00541         ParamsToArray(pai->m_vparam, m_vregParams[nSource]);
00542 
00543         // TODO: think about m_rStep interacting with m_regSteps.V()
00544         pai->m_rErr = NelderMeadOptimize
00545         (
00546             m_vregParams[nSource].C(),
00547             pai->m_vparam,
00548             m_regSteps.V(),
00549             ( (sError == "parzen")
00550                 ? MeasureInverseParzenProbability
00551                 : MeasureVariance),
00552             pai,
00553             .0001,
00554             10
00555         );
00556 
00557         ParamsFromArray(m_vregParams[nSource],pai->m_vparam);
00558         }
00559 
00560       // Unknown algorithm -------------------------------------------
00561 
00562       else
00563         {
00564         P(DescribeConfiguration());
00565         ERROR(
00566             "Unknown algorithm."
00567             "Parameter congeal.optimize.algorithm must be one of {lbfgs|gradientdescent|neldermead}."
00568             "Algorithm given: '%s'", sAlgorithm.VCH()
00569         );
00570         }
00571 
00572 #if MULTITHREADED
00573       pthread_barrier_wait(&m_barrierDone);
00574 #else
00575       return;
00576 #endif
00577       }
00578     }
00579 
00580   static void ParamsFromArray(RegistryOfParameters& reg, const Parameter* vparam)
00581     {
00582     int c=reg.C();
00583     for (int n=0; n<c; n++)
00584       {
00585       (*(reg.Get(n))) = vparam[n];
00586       }
00587     }
00588 
00589   static void ParamsToArray(Parameter* vparam, const RegistryOfParameters& reg)
00590     {
00591     int c=reg.C();
00592     for (int n=0; n<c; n++)
00593       {
00594       vparam[n] = (*(reg.Get(n)));
00595       }
00596     }
00597 
00598   //--------------------------------------------------------------------------
00599   //--------------------------------------------------------------------------
00600 #if FEATURE_LBFGS
00601   static lbfgsfloatval_t EvaluateForLBFGS(
00602       void *pObjectiveInfo,
00603       const lbfgsfloatval_t *vparam,
00604       lbfgsfloatval_t *vparamDerivative,
00605       const int cParams,
00606       const lbfgsfloatval_t nStep
00607   )
00608     {
00609     AlignInfo* pai=(AlignInfo*)pObjectiveInfo;
00610     const Parameter* vSteps = pai->m_pcongeal->m_regSteps.V();
00611 
00612     static Real rKernelFactor=ConfigValue("congeal.optimize[lbfgs].kernel",30.);
00613     Real rErr=MeasureInverseParzenProbability(vparam, REAL_MAX, pObjectiveInfo);
00614 
00615     for (int n=0; n<cParams; n++)
00616       {
00617       Parameter dKernel=vSteps[n] * rKernelFactor;
00618       Real rErrP1;
00619       // WARNING: making a const non-const for a moment
00620 
00621         {
00622         // change the value temporarily
00623         ((Parameter*)vparam)[n] += dKernel;
00624         rErrP1=MeasureInverseParzenProbability(vparam, REAL_MAX, pObjectiveInfo);
00625         // get back to original value
00626         ((Parameter*)vparam)[n] -= dKernel;
00627         }
00628 
00629       // compute gradient
00630       vparamDerivative[n]=(rErrP1-rErr)/dKernel;
00631 
00632       /*
00633        Using central difference is too expensive,
00634        requiring 2x function evaluations
00635        vparam[n] -= 2*dKernel;
00636        Real rErrM1=MeasureInverseParzenProbability(vparam, REAL_MAX, pObjectiveInfo);
00637 
00638        // get back to original value
00639        vparam[n] += dKernel;
00640 
00641        vparamDerivative[n]=(rErrP1-rErrM1)/dKernel/2;
00642        */
00643       }
00644     return rErr;
00645     }
00646 #endif
00647 
00648   //--------------------------------------------------------------------------
00649   //--------------------------------------------------------------------------
00650   static Real Optimize_BruteForce(
00651       int cParams,
00652       Parameter** vpParams,
00653       Parameter* vSteps,
00654       Real fxnObjective(const void*),
00655       void* pObjectiveInfo,
00656       Real rStep
00657   )
00658     {
00659     static Real rRangeFactor=ConfigValue("congeal.optimize[bruteforce].range",30.);
00660 
00661     Real rErrMin=REAL_MAX;
00662 
00663     for (int n=0; n<cParams; n++)
00664       {
00665       int dMin=0;
00666       rErrMin=REAL_MAX;
00667       for (int d=-20; d<20; d++)
00668         {
00669 
00670         Parameter dRange=vSteps[n] * rRangeFactor * d;
00671         (*(vpParams[n])) += dRange;
00672         Real rErrP1=fxnObjective(pObjectiveInfo);
00673 
00674         if (rErrP1<rErrMin)
00675           {
00676           rErrMin=rErrP1;
00677           dMin=d;
00678           }
00679         //              UD("ERROR ALONG %5d: delta of %5.5f --> %5.5f == %5.5f", n, dRange, (*(vpParams[n])), rErrP1);
00680 
00681         (*(vpParams[n])) -= dRange;
00682         }
00683 
00684       Parameter dRange=vSteps[n] * rRangeFactor * dMin;
00685       (*(vpParams[n])) += dRange;
00686       }
00687     return rErrMin;
00688     }
00689 
00690   //--------------------------------------------------------------------------
00691   //--------------------------------------------------------------------------
00692   static Real Optimize_GreedyRandomWalk(
00693       int cParams,
00694       Parameter** vpParams,
00695       Parameter* vSteps,
00696       Real fxnObjective(const void*),
00697       void* pObjectiveInfo,
00698       Real rStep
00699   )
00700     {
00701     static Real rKernelFactor=ConfigValue("congeal.optimize[randomwalk].kernel",30.);
00702     static int cSteps=ConfigValue("congeal.optimize[randomwalk].steps",1000);
00703     static int cDirections=ConfigValue("congeal.optimize[randomwalk].directions",100);
00704 
00705     /*
00706      Real rErrTotalIn;
00707      UDTIMEEXECUTION(
00708      "Time to measure full error %g",
00709      rErrTotalIn=MeasureFullInverseParzenProbabilityInPlace(pObjectiveInfo)
00710      );
00711      */
00712 
00713     Real rErrIn=fxnObjective(pObjectiveInfo);
00714     Real rErrMin=rErrIn;
00715 
00716     for (int nDir=0; nDir<cDirections; nDir++)
00717       {
00718       int nDim=rand()%cParams;
00719       Parameter dStep=vSteps[nDim] * rStep * rKernelFactor;
00720       dStep *= Random(-1.,1.);
00721 
00722       Parameter dStepThreshold = vSteps[nDim] / 100000;
00723 
00724       for (int nStep=0; nStep<cSteps; nStep++)
00725         {
00726         //              Real paramFrom=(*(vpParams[nDim]));
00727         (*(vpParams[nDim])) += dStep;
00728         //              Real paramTo=(*(vpParams[nDim]));
00729 
00730         Real rErrNew=fxnObjective(pObjectiveInfo);
00731 
00732         if (rErrNew>=rErrMin)
00733           {
00734           //                    UD("     wrong %g + (%g) == %g err %g to %g", paramFrom, dStep, paramTo, rErrMin, rErrNew);
00735           (*(vpParams[nDim])) -= dStep;
00736           dStep=-dStep*.9;
00737           if (dStep > dStepThreshold || dStep < -dStepThreshold)
00738             {
00739             nStep--;
00740             }
00741           }
00742         else
00743           {
00744           rErrMin=rErrNew;
00745           //                    UD("[%d] +++RIGHT %g + (%g) == %g ERR %g TO %g", nDim,paramFrom, dStep, paramTo, rErrMin, rErrNew);
00746           }
00747         }
00748       }
00749 
00750     /*
00751      Real rErrTotalNew=MeasureFullInverseParzenProbabilityInPlace(pObjectiveInfo);
00752 
00753      UD("MOVED FROM STOCERROR(%g-->%g) TOTALERROR(%g-->%g) in %d directions",
00754      rErrIn,rErrMin,rErrTotalIn,rErrTotalNew,cDirections
00755      );
00756      */
00757 
00758     return rErrMin;
00759     }
00760 
00761   //--------------------------------------------------------------------------
00762   //--------------------------------------------------------------------------
00763   static Real Optimize_OneGradientStep(
00764       int cParams,
00765       Parameter** vpParams,
00766       Parameter* vSteps,
00767       Real fxnObjective(const void*),
00768       void* pObjectiveInfo,
00769       Real rStep
00770   )
00771     {
00772     static Real rKernelFactor=ConfigValue("congeal.optimize[gradientdescent].kernel",30.);
00773     Real rErr=fxnObjective(pObjectiveInfo);
00774 
00775     // initialize static buffer for holding derivative
00776     Parameter* vparamDerivative=new Parameter[cParams];
00777 
00778     for (int n=0; n<cParams; n++)
00779       {
00780       Parameter dKernel=vSteps[n] * rStep * rKernelFactor;
00781 
00782       (*(vpParams[n])) += dKernel;
00783       Real rErrP1=fxnObjective(pObjectiveInfo);
00784 
00785       UD("LOOKING AT GRADIENT ALONG %d SHIFTED BY %g GIVES %g-->%g", n, dKernel,rErr,rErrP1);
00786       // get back to original value
00787       (*(vpParams[n])) -= dKernel;
00788 
00789       // compute gradient
00790       vparamDerivative[n]=(rErrP1-rErr)/dKernel;
00791 
00792       /*
00793        Using central difference is too expensive,
00794        requiring 2x function evaluations
00795        (*(vpParams[n])) -= 2*dKernel;
00796        Real rErrM1=fxnObjective(pObjectiveInfo);
00797 
00798        // get back to original value
00799        (*(vpParams[n])) += dKernel;
00800 
00801        UD("LOOKING AT GRADIENT ALONG %d SHIFTED BY %g GIVES %g-->%g", n, -dKernel,rErr,rErrM1);
00802 
00803        vparamDerivative[n]=(rErrP1-rErrM1)/dKernel/2;
00804        */
00805       }
00806     //      AlignInfo* pai=(AlignInfo*)pObjectiveInfo;
00807 
00808     // normalize the gradient
00809     /*
00810      Real rT=0;
00811      for (int n=0; n<cParams; n++){
00812      rT += sqr(vparamDerivative[n]);
00813      }
00814      rT=sqrt(rT);
00815      for (int n=0; n<cParams; n++){
00816      vparamDerivative[n]/=rT;
00817      }
00818      */
00819 
00820     static Real rStepFactor=ConfigValue("congeal.optimize[gradientdescent].step",30.);
00821     static bool bWalkAlongGradient=ConfigValue("congeal.optimize[gradientdescent].walk",true);
00822 
00823     if (!bWalkAlongGradient)
00824       {
00825       for (int n=0; n<cParams; n++)
00826         {
00827         (*(vpParams[n])) += -vparamDerivative[n] * rStep * rStepFactor;
00828         }
00829       Real rErrNew=fxnObjective(pObjectiveInfo);
00830       delete [] vparamDerivative;
00831       return rErrNew;
00832       }
00833 
00834     static int cWalkSteps=ConfigValue("congeal.optimize[gradientdescent].walk.steps",30);
00835 
00836     Real rdStep=1.5 * rStep * rStepFactor / cWalkSteps;
00837     Real rErrMin=rErr;
00838     int nWalkStepMin=-1;
00839 
00840     // note we're taking the negative of the derivative here
00841     for (int n=0; n<cParams; n++)
00842       {
00843       vparamDerivative[n] = -vparamDerivative[n] * rdStep;
00844       }
00845     for (int nWalkStep=0; nWalkStep<cWalkSteps; nWalkStep++)
00846       {
00847       for (int n=0; n<cParams; n++)
00848         {
00849         (*(vpParams[n])) += vparamDerivative[n];
00850         }
00851       Real rErrNew=fxnObjective(pObjectiveInfo);
00852       if (rErrNew<rErrMin)
00853         {
00854         UD("%d] %g < %g ***", nWalkStep, rErrNew, rErrMin);
00855         rErrMin=rErrNew;
00856         nWalkStepMin=nWalkStep;
00857         }
00858       else
00859         {
00860         UD("%d] %g", nWalkStep, rErrNew);
00861         }
00862       }
00863 
00864     //      UD("min error was %g compared to %g", rErrMin, rErr);
00865     //      UD("shifting by: %d of %d",  (cWalkSteps - (nWalkStepMin+1)), cWalkSteps);
00866     for (int n=0; n<cParams; n++)
00867       {
00868       (*(vpParams[n])) -= vparamDerivative[n] * (cWalkSteps - (nWalkStepMin+1));
00869       }
00870     Real rErrNew=fxnObjective(pObjectiveInfo);
00871 
00872     UD("%g --> %g %s", rErr, rErrNew, (rErr==rErrNew?"":"*"));
00873     delete [] vparamDerivative;
00874     return rErrNew;
00875     }
00876 
00877   //--------------------------------------------------------------------------
00878   //--------------------------------------------------------------------------
00879   // note this function copies the parameters from the AlignInfo into
00880   // the recipie and then measures the variance
00881   static Real MeasureVariance(const Parameter* vparam, const Real& rErrMin, const void* p)
00882     {
00883     AlignInfo* pai = (AlignInfo*)p;
00884 
00885     ParamsFromArray(
00886         pai->m_pcongeal->m_vregParams[pai->m_nSource],
00887         vparam
00888     );
00889     return MeasureVarianceInPlace(p);
00890     }
00891 
00892   //--------------------------------------------------------------------------
00893   //--------------------------------------------------------------------------
00894   // note this function assumes parameter values are already inserted into
00895   // the recipie
00896   static Real MeasureVarianceInPlace(const void* p)
00897     {
00898     AlignInfo* pai = (AlignInfo*)p;
00899 
00900     const int& nSource = pai->m_nSource;
00901     const int& cSamples=pai->m_pcongeal->m_cSamples;
00902     const int& cSources=pai->m_pcongeal->m_cSources;
00903 
00904     // get ready to access
00905     pai->m_pcongeal->m_precipie->PrepareForAccess(nSource);
00906 
00907     Real rErrT=0;
00908 
00909     for (int nSample=0; nSample<cSamples; nSample++)
00910       {
00911 
00912       DATA dataReplacement;
00913       ::GetPoint(
00914           *(pai->m_pcongeal->m_precipie->PSource(nSource)),
00915           dataReplacement,
00916           pai->m_pcongeal->m_vpt[nSample]
00917       );
00918 
00919       // compute a new sum and sumsq by taking out the excluded data
00920       // (which is presumably the source we're modifying)
00921       // and adding in the new data
00922       const DATA& dataExcluded = pai->m_pcongeal->m_mdata[nSample*cSources+nSource];
00923 
00924       const RANGEPROMOTION(DATA)& dataSum = (
00925           pai->m_pcongeal->m_vdataSum[nSample]
00926           - dataExcluded
00927           + dataReplacement
00928       );
00929 
00930       const RANGEPROMOTION(DATA)& dataSumSq = (
00931           pai->m_pcongeal->m_vdataSumSq[nSample]
00932           - sqr(dataExcluded)
00933           + sqr(dataReplacement)
00934       );
00935 
00936       // TODO: apply Length() to this for general purpose computation
00937       Real rErr = dataSumSq/cSources - sqr(dataSum/cSources);
00938 
00939       rErrT+=rErr;
00940       }
00941     rErrT /= cSamples;
00942 
00943     return rErrT;
00944     }
00945 
00946   //--------------------------------------------------------------------------
00947   //--------------------------------------------------------------------------
00948   // note this function copies the parameters from the AlignInfo into
00949   // the recipie and then measures the variance
00950   static Real MeasureInverseParzenProbability(const Parameter* vparam, const Real& rErrMin, const void* p)
00951     {
00952     AlignInfo* pai = (AlignInfo*)p;
00953 
00954     ParamsFromArray(
00955         pai->m_pcongeal->m_vregParams[pai->m_nSource],
00956         vparam
00957     );
00958     return MeasureInverseParzenProbabilityInPlace(p);
00959     }
00960 
00961 #define FASTGAUSSIAN_TABLESIZE 1024
00962 #define FASTGAUSSIAN_TABLESIZE_DIV_2 512
00963 
00964   static Real m_vrFastGaussian[FASTGAUSSIAN_TABLESIZE];
00965   //--------------------------------------------------------------------------
00966   //--------------------------------------------------------------------------
00967   static void InitFastGaussian(Real &rSigma)
00968     {
00969     D("InitFastGaussian");
00970     for (int n=0; n<FASTGAUSSIAN_TABLESIZE; n++)
00971       {
00972       m_vrFastGaussian[n]=Gaussian(n-FASTGAUSSIAN_TABLESIZE_DIV_2, 0, rSigma);
00973       }
00974     }
00975   //--------------------------------------------------------------------------
00976   //--------------------------------------------------------------------------
00977   static Real FastGaussian(const DATA &r, const DATA &rMean, Real &rSigma)
00978     {
00979     int n=Bound(rMean-r + FASTGAUSSIAN_TABLESIZE_DIV_2,0,FASTGAUSSIAN_TABLESIZE-1);
00980     return m_vrFastGaussian[n];
00981     }
00982 
00983   //--------------------------------------------------------------------------
00984   //--------------------------------------------------------------------------
00985   // note this function assumes parameter values are already inserted into
00986   // the recipie
00987   static Real MeasureInverseParzenProbabilityInPlace(const void* p)
00988     {
00989     static Real rSigma=ConfigValue("congeal.error[parzen].sigma",10.);
00990     AlignInfo* pai = (AlignInfo*)p;
00991 
00992     const int& nSourceTest = pai->m_nSource;
00993     const int& cSamples=pai->m_pcongeal->m_cSamples;
00994     const int& cSources=pai->m_pcongeal->m_cSources;
00995 
00996     // get ready to access
00997     pai->m_pcongeal->m_precipie->PrepareForAccess(nSourceTest);
00998 
00999     Real rLogProbTotal=1.;
01000 
01001     int cConsideredSamples=0;
01002     //      UD(String("COMPARING SAMPLES: ") + cSamples);
01003     for (int nSample=0; nSample<cSamples; nSample++)
01004       {
01005       DATA dataTest;
01006       ::GetPoint(
01007           *(pai->m_pcongeal->m_precipie->PSource(nSourceTest)),
01008           dataTest,
01009           pai->m_pcongeal->m_vpt[nSample]
01010       );
01011       m_cLookups++;
01012 
01013       /*
01014 
01015        if (dataTest==0){
01016        continue;
01017        }
01018 
01019 
01020        while (dataTest==0){
01021        ::GetPoint(
01022        *(pai->m_pcongeal->m_precipie->PSource(nSourceTest)),
01023        dataTest,
01024        NewRandomPoint(pai->m_pcongeal->m_precipie->Size())
01025        );
01026        m_cLookups++;
01027        }
01028        */
01029 
01030       cConsideredSamples++;
01031       Real rLogProbSample=0;
01032       int cConsideredSources=0;
01033 
01034       for (int nSource=0; nSource<cSources; nSource++)
01035         {
01036         if (nSource==nSourceTest)
01037           {
01038           continue;
01039           }
01040 
01041         DATA dataPrior = pai->m_pcongeal->m_mdata[nSample*cSources+nSource];
01042 
01043         if (dataPrior==0)
01044           {
01045           UD(String("Hmm? dataPrior is 0?") + pai->m_pcongeal->m_vpt[nSample]);
01046           continue;
01047           }
01048 
01049         /*
01050          int cRelookups=0;
01051          while (dataPrior==0){
01052          ::GetPoint(
01053          *(pai->m_pcongeal->m_precipie->PSource(nSource)),
01054          dataPrior,
01055          NewRandomPoint(pai->m_pcongeal->m_precipie->Size())
01056          );
01057          m_cLookups++;
01058          cRelookups++;
01059          D("relookups: %d", cRelookups);
01060          }
01061          */
01062 
01063         cConsideredSources++;
01064 
01065         Real rProbSource = Gaussian(dataTest, dataPrior, rSigma);
01066 
01067         ASSERTf(rProbSource>=0, "NEGATIVE ERROR %g", rProbSource);
01068 
01069         static Real rApriori=ConfigValue("congeal.error[parzen].apriori",.2);
01070         Real rLogProbSource = -log(rProbSource + rApriori);
01071         rLogProbSample += rLogProbSource;
01072 
01073         //              UD(String("COMPARING ") + dataTest + " with " + dataPrior + " getting " + rProbSource + " logs to " + rLogProbSource + " and combines to " + rLogProbSample);
01074         }
01075 
01076       if (cConsideredSources>0)
01077         {
01078         rLogProbTotal +=rLogProbSample/cConsideredSources;
01079         }
01080 
01081       //            UD(String("SOURCES COMPARED: ") + cConsideredSources);
01082       }
01083 
01084     if (cConsideredSamples>0)
01085       {
01086       rLogProbTotal /= cConsideredSamples;
01087       }
01088     //      UD(String("SAMPLES COMPARED: ") + cConsideredSamples + " yields " + rLogProbTotal);
01089 
01090     Real rError = rLogProbTotal;
01091 
01092     return rError;
01093     }
01094 
01095   //--------------------------------------------------------------------------
01096   //--------------------------------------------------------------------------
01097   // note this function assumes parameter values are already inserted into
01098   // the recipie
01099   static Real MeasureFullInverseParzenProbabilityInPlace(const void* p)
01100     {
01101     static Real rSigma=ConfigValue("congeal.error[parzen].sigma",10.);
01102     AlignInfo* pai = (AlignInfo*)p;
01103 
01104     const int& nSourceTest = pai->m_nSource;
01105     const int& cSources=pai->m_pcongeal->m_cSources;
01106 
01107     // get ready to access
01108     pai->m_pcongeal->m_precipie->PrepareForAccess(nSourceTest);
01109 
01110     Real rLogProbTotal=1.;
01111 
01112     // loop over all the points and measure the error
01113     forpoint(TypeOfPoint(RECIPIE),pt,0,(pai->m_pcongeal->m_precipie->PSource(nSourceTest))->Size())
01114       {
01115       DATA dataTest;
01116       ::GetPoint(
01117           *(pai->m_pcongeal->m_precipie->PSource(nSourceTest)),
01118           dataTest,
01119           pt
01120       );
01121 
01122       Real rLogProbSample=0;
01123       for (int nSource=0; nSource<cSources; nSource++)
01124         {
01125         if (nSource==nSourceTest)
01126           {
01127           continue;
01128           }
01129         DATA dataPrior;
01130         ::GetPoint(
01131             *(pai->m_pcongeal->m_precipie->PSource(nSource)),
01132             dataPrior,
01133             pt
01134         );
01135 
01136         Real rProbSource = FastGaussian(dataTest, dataPrior, rSigma);
01137 
01138         ASSERTf(rProbSource>=0, "NEGATIVE ERROR %g", rProbSource);
01139 
01140         static Real rApriori=ConfigValue("congeal.error[parzen].apriori",.2);
01141         rLogProbSample += -log(rProbSource + rApriori);
01142         }
01143 
01144       rLogProbTotal +=rLogProbSample;
01145       }
01146 
01147     rLogProbTotal /= (pai->m_pcongeal->m_precipie->PSource(nSourceTest))->CSize();
01148     rLogProbTotal /= cSources;
01149 
01150     Real rError = rLogProbTotal;
01151 
01152     return rError;
01153     }
01154 
01155   };
01156 
01157 template<class RECIPIE>
01158 long int CongealGroupPhasedOf<RECIPIE>::m_cLookups;
01159 
01160 #endif //__CONGEALGROUPPHASED_H__
01161 

Generated on 6 Apr 2011 for Slicer3 by  doxygen 1.6.1