Class for computing the negitve free energy using thermodynamic integration. More...
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 |
MT19937 & | get_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 |
MT19937 & | rng |
Class for computing the negitve free energy using thermodynamic integration.
Definition at line 258 of file multispot5.cc.
NegativeFreeEnergy::NegativeFreeEnergy | ( | const DataForMCMC & | d | ) | [inline] |
d | Necessary data |
Definition at line 263 of file multispot5.cc.
00264 :DataForMCMC(d) 00265 { 00266 }
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.
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()().
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
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
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 }
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().