NegativeFreeEnergy Class Reference

Class for computing the negitve free energy using thermodynamic integration. More...

Inheritance diagram for NegativeFreeEnergy:
DataForMCMC

List of all members.

Public Member Functions

 NegativeFreeEnergy (const DataForMCMC &d)
double variance_from_sample (double sample, double samples, double base_sigma, double scale_pow) const
double compute_with_mask (const Vector<> &spots, const vector< vector< int > > &spot_pixels) const
double operator() (const Vector<> &spots) const
MT19937get_rng () const

Protected Attributes

const vector< ImageRef > & pixels
const vector< vector< double > > & pixel_intensities
const double mu_brightness
const double sigma_brightness
const double mu_blur
const double sigma_blur
const double variance
const int samples
const int sample_iterations
const Matrix< 3 > A
const Vector< 3 > pi
MT19937rng

Detailed Description

Class for computing the negitve free energy using thermodynamic integration.

Definition at line 258 of file multispot5.cc.


Constructor & Destructor Documentation

NegativeFreeEnergy::NegativeFreeEnergy ( const DataForMCMC d  )  [inline]
Parameters:
d Necessary data

Definition at line 263 of file multispot5.cc.

00264     :DataForMCMC(d)
00265     {
00266     }


Member Function Documentation

double NegativeFreeEnergy::variance_from_sample ( double  sample,
double  samples,
double  base_sigma,
double  scale_pow 
) const [inline]

Give the noise variance given a sample number and growth parameters.

Parameters:
sample Sample number
samples Total number of samples
base_sigma Starting value of sigme
scale_pow Exponent scaling

Definition at line 273 of file multispot5.cc.

References sq().

Referenced by compute_with_mask(), and operator()().

00274     {
00275         double scale = pow(1.25, sample * 1. / samples * scale_pow);
00276         double sigma = base_sigma * scale;
00277         double new_variance = sq(sigma);
00278 
00279         return new_variance;
00280     }

double NegativeFreeEnergy::compute_with_mask ( const Vector<> &  spots,
const vector< vector< int > > &  spot_pixels 
) const [inline]

Estimate free energy using the Slow Growth Thermodynamic Integration method given in Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford.

M. Neal, 1993 Except using a 5 point stencil instead of forward differenceing in Eq 6.30

Parameters:
spots list of spots
spot_pixels Mask around each spot to use to save on computation

Definition at line 287 of file multispot5.cc.

References DataForMCMC::A, Kahan::add(), assert_same_size(), SampledMultispot::compute_spot_intensity(), DataForMCMC::get_rng(), log_log_normal(), DataForMCMC::mu_blur, DataForMCMC::mu_brightness, DataForMCMC::pi, DataForMCMC::pixel_intensities, DataForMCMC::pixels, DataForMCMC::sample_iterations, DataForMCMC::samples, DataForMCMC::sigma_blur, DataForMCMC::sigma_brightness, spots_to_vector(), sq(), DataForMCMC::variance, and variance_from_sample().

00288     {
00289         //Estimate free energy using the Slow Growth Thermodynamic Integration method given in
00290         //Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford. M. Neal, 1993
00291         //Except using a 5 point stencil instead of forward differenceing in Eq 6.30
00292         double base_sigma = sqrt(variance); // should be 1
00293         double scale_pow = 100.0;
00294 
00295         const unsigned int nspots  = spots.size()/4;
00296         const unsigned int nframes = pixel_intensities.size();
00297         const unsigned int npixels = pixels.size();
00298         assert(spots.size() %4 == 0);
00299         assert(spot_pixels.size() == nspots);
00300         assert_same_size(spot_pixels);
00301 
00302         //Compute the intensities and derivatives for all spot
00303         vector<vector<double> > spot_intensity; //[spot][pixel]
00304         for(unsigned int i=0; i < nspots; i++)
00305             spot_intensity.push_back(compute_spot_intensity(pixels, spots.slice<Dynamic,4>(i*4,4)));
00306         
00307         GibbsSampler2 sampler(pixel_intensities, spot_intensity, spots_to_vector(spots), spot_pixels, A, pi, variance, sample_iterations);
00308 
00309         double sum = 0;
00310         Kahan ksum;
00311         for(int sample=0; sample < samples; sample++)
00312         {
00313             //Compute the positions of the surrounding steps
00314             double var1 = variance_from_sample(sample-2, samples, base_sigma, scale_pow);
00315             double var2 = variance_from_sample(sample-1, samples, base_sigma, scale_pow);
00316             double var3 = variance_from_sample(sample+1, samples, base_sigma, scale_pow);
00317             double var4 = variance_from_sample(sample+2, samples, base_sigma, scale_pow);
00318 
00319             //Take a sample
00320             sampler.set_variance(var2);
00321             sampler.next(DataForMCMC::get_rng());
00322             
00323             //Compute the SSD. This is fixed regardless of sigma.   
00324             double err_sum=0;
00325             for(unsigned int frame=0; frame < nframes; frame++)
00326                 for(unsigned int pixel=0; pixel < npixels; pixel++)
00327                     err_sum -= sq(pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel]);
00328             
00329             //Compute the derivative using a five point stencil
00330             //This could be done in a better but less clear way
00331             double e1 = err_sum / (2*var1) - npixels*nframes*::log(2*M_PI*var1)/2;
00332             double e2 = err_sum / (2*var2) - npixels*nframes*::log(2*M_PI*var2)/2;
00333             double e3 = err_sum / (2*var3) - npixels*nframes*::log(2*M_PI*var3)/2;
00334             double e4 = err_sum / (2*var4) - npixels*nframes*::log(2*M_PI*var4)/2;
00335             sum += (-e1 + 8*e2 - 8*e3 + e4)/12;
00336 
00337             ksum.add(-e1/12);
00338             ksum.add(8*e2/12);
00339             ksum.add(-8*e3/12);
00340             ksum.add(e4/12);
00341         }
00342 
00343         double log_final = (log(variance_from_sample(samples, samples, base_sigma, scale_pow)*2*M_PI)/2) * npixels * nframes;
00344         
00345         double priors=0;
00346         for(unsigned int i=0; i < nspots; i++)
00347         {
00348             priors += log_log_normal(spots[i*4+0], mu_brightness, sigma_brightness);
00349             priors += log_log_normal(spots[i*4+1], mu_blur, sigma_blur);
00350         }
00351 
00352         /*cout << "Thermo:\n";
00353         cout << "sum: " << sum -log_final << endl;
00354         cout << "kah: " << ksum.sum -log_final << endl;
00355         cout << "priors: " << priors + sum -log_final << endl;
00356         cout << "   kah: " << priors + ksum.sum -log_final << endl;
00357 */
00358         //cout << log_final << endl;
00359         //cout << sum + log_final << endl;
00360         
00361         
00362         sampler.set_variance(variance);
00363         return -(sum+priors - log_final);
00364     }

double NegativeFreeEnergy::operator() ( const Vector<> &  spots  )  const [inline]

Estimate free energy using the Slow Growth Thermodynamic Integration method given in Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford.

M. Neal, 1993 Except using a 5 point stencil instead of forward differenceing in Eq 6.30

Parameters:
spots list of spots

Definition at line 370 of file multispot5.cc.

References DataForMCMC::A, Kahan::add(), SampledMultispot::compute_spot_intensity(), DataForMCMC::get_rng(), log_log_normal(), DataForMCMC::mu_blur, DataForMCMC::mu_brightness, DataForMCMC::pi, DataForMCMC::pixel_intensities, DataForMCMC::pixels, DataForMCMC::sample_iterations, DataForMCMC::samples, DataForMCMC::sigma_blur, DataForMCMC::sigma_brightness, spots_to_vector(), sq(), DataForMCMC::variance, and variance_from_sample().

00371     {
00372         double base_sigma = sqrt(variance); // should be 1
00373         double scale_pow = 100.0;
00374 
00375         const unsigned int nspots  = spots.size()/4;
00376         const unsigned int nframes = pixel_intensities.size();
00377         const unsigned int npixels = pixels.size();
00378         assert(spots.size() %4 == 0);
00379 
00380         //Compute the intensities and derivatives for all spot
00381         vector<vector<double> > spot_intensity; //[spot][pixel]
00382         for(unsigned int i=0; i < nspots; i++)
00383             spot_intensity.push_back(compute_spot_intensity(pixels, spots.slice<Dynamic,4>(i*4,4)));
00384         
00385         GibbsSampler sampler(pixel_intensities, spot_intensity, spots_to_vector(spots), A, pi, variance, sample_iterations);
00386 
00387         double sum = 0;
00388         Kahan ksum;
00389         for(int sample=0; sample < samples; sample++)
00390         {
00391             //Compute the positions of the surrounding steps
00392             double var1 = variance_from_sample(sample-2, samples, base_sigma, scale_pow);
00393             double var2 = variance_from_sample(sample-1, samples, base_sigma, scale_pow);
00394             double var3 = variance_from_sample(sample+1, samples, base_sigma, scale_pow);
00395             double var4 = variance_from_sample(sample+2, samples, base_sigma, scale_pow);
00396 
00397             //Take a sample
00398             sampler.set_variance(var2);
00399             sampler.next(DataForMCMC::get_rng());
00400             
00401             //Compute the SSD. This is fixed regardless of sigma.   
00402             double err_sum=0;
00403             for(unsigned int frame=0; frame < nframes; frame++)
00404                 for(unsigned int pixel=0; pixel < npixels; pixel++)
00405                     err_sum -= sq(pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel]);
00406             
00407             //Compute the derivative using a five point stencil
00408             //This could be done in a better but less clear way
00409             double e1 = err_sum / (2*var1) - npixels*nframes*::log(2*M_PI*var1)/2;
00410             double e2 = err_sum / (2*var2) - npixels*nframes*::log(2*M_PI*var2)/2;
00411             double e3 = err_sum / (2*var3) - npixels*nframes*::log(2*M_PI*var3)/2;
00412             double e4 = err_sum / (2*var4) - npixels*nframes*::log(2*M_PI*var4)/2;
00413             sum += (-e1 + 8*e2 - 8*e3 + e4)/12;
00414 
00415             ksum.add(-e1/12);
00416             ksum.add(8*e2/12);
00417             ksum.add(-8*e3/12);
00418             ksum.add(e4/12);
00419         }
00420 
00421         double log_final = (log(variance_from_sample(samples, samples, base_sigma, scale_pow)*2*M_PI)/2) * npixels * nframes;
00422         
00423         double priors=0;
00424         for(unsigned int i=0; i < nspots; i++)
00425         {
00426             priors += log_log_normal(spots[i*4+0], mu_brightness, sigma_brightness);
00427             priors += log_log_normal(spots[i*4+1], mu_blur, sigma_blur);
00428         }
00429 
00430         /*cout << "Thermo:\n";
00431         cout << "sum: " << sum -log_final << endl;
00432         cout << "kah: " << ksum.sum -log_final << endl;
00433         cout << "priors: " << priors + sum -log_final << endl;
00434         cout << "   kah: " << priors + ksum.sum -log_final << endl;
00435 */
00436         //cout << log_final << endl;
00437         //cout << sum + log_final << endl;
00438         
00439         
00440         sampler.set_variance(variance);
00441         return -(sum+priors - log_final);
00442     }

MT19937& DataForMCMC::get_rng (  )  const [inline, inherited]

Definition at line 186 of file multispot5.cc.

References DataForMCMC::rng.

Referenced by compute_with_mask(), and operator()().

00187     {
00188         return rng;
00189     }


Member Data Documentation

const vector<ImageRef>& DataForMCMC::pixels [protected, inherited]

Definition at line 176 of file multispot5.cc.

Referenced by compute_with_mask(), FreeEnergyHessian::hessian(), and operator()().

const vector<vector<double> >& DataForMCMC::pixel_intensities [protected, inherited]

Definition at line 177 of file multispot5.cc.

Referenced by compute_with_mask(), FreeEnergyHessian::hessian(), and operator()().

const double DataForMCMC::mu_brightness [protected, inherited]

Definition at line 178 of file multispot5.cc.

Referenced by compute_with_mask(), FreeEnergyHessian::hessian(), and operator()().

const double DataForMCMC::sigma_brightness [protected, inherited]

Definition at line 178 of file multispot5.cc.

Referenced by compute_with_mask(), FreeEnergyHessian::hessian(), and operator()().

const double DataForMCMC::mu_blur [protected, inherited]

Definition at line 178 of file multispot5.cc.

Referenced by compute_with_mask(), FreeEnergyHessian::hessian(), and operator()().

const double DataForMCMC::sigma_blur [protected, inherited]

Definition at line 178 of file multispot5.cc.

Referenced by compute_with_mask(), FreeEnergyHessian::hessian(), and operator()().

const double DataForMCMC::variance [protected, inherited]

Definition at line 179 of file multispot5.cc.

Referenced by compute_with_mask(), FreeEnergyHessian::hessian(), and operator()().

const int DataForMCMC::samples [protected, inherited]

Definition at line 180 of file multispot5.cc.

Referenced by compute_with_mask(), FreeEnergyHessian::hessian(), and operator()().

const int DataForMCMC::sample_iterations [protected, inherited]

Definition at line 180 of file multispot5.cc.

Referenced by compute_with_mask(), FreeEnergyHessian::hessian(), and operator()().

const Matrix<3> DataForMCMC::A [protected, inherited]

Definition at line 181 of file multispot5.cc.

Referenced by compute_with_mask(), FreeEnergyHessian::hessian(), and operator()().

const Vector<3> DataForMCMC::pi [protected, inherited]

Definition at line 182 of file multispot5.cc.

Referenced by compute_with_mask(), FreeEnergyHessian::hessian(), and operator()().

MT19937& DataForMCMC::rng [mutable, protected, inherited]

Definition at line 183 of file multispot5.cc.

Referenced by DataForMCMC::get_rng(), and FreeEnergyHessian::hessian().


The documentation for this class was generated from the following file:
Generated on Wed Nov 2 18:00:03 2011 for BCUBED by  doxygen 1.6.3