Class for computing the Hessian of the negative free energy. More...
Public Member Functions | |
FreeEnergyHessian (const DataForMCMC &d) | |
Matrix< 4 > | hessian (const vector< Vector< 4 > > &spots, int spot) 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 Hessian of the negative free energy.
Definition at line 606 of file multispot5.cc.
FreeEnergyHessian::FreeEnergyHessian | ( | const DataForMCMC & | d | ) | [inline] |
Constructor.
d | All data required |
Definition at line 612 of file multispot5.cc.
00613 :DataForMCMC(d) 00614 { 00615 }
Matrix<4> FreeEnergyHessian::hessian | ( | const vector< Vector< 4 > > & | spots, | |
int | spot | |||
) | const [inline] |
Compute the Hessian.
spots | All spot positions | |
spot | spot to compute hessian for |
Definition at line 620 of file multispot5.cc.
References DataForMCMC::A, SampledMultispot::compute_spot_intensity(), SampledMultispot::compute_spot_intensity_hessian(), diff_log_log_normal(), hess_log_log_normal(), DataForMCMC::mu_blur, DataForMCMC::mu_brightness, DataForMCMC::pi, DataForMCMC::pixel_intensities, DataForMCMC::pixels, DataForMCMC::rng, DataForMCMC::sample_iterations, DataForMCMC::samples, DataForMCMC::sigma_blur, DataForMCMC::sigma_brightness, and DataForMCMC::variance.
Referenced by fit_spots_historic(), and FitSpots::optimize_each_spot_in_turn_for_several_passes().
00621 { 00622 cout << "----Computing pure MCMC hessian\n"; 00623 const unsigned int nspots = spots.size(); 00624 const unsigned int nframes = pixel_intensities.size(); 00625 const unsigned int npixels = pixels.size(); 00626 cout << spot << " " << nspots << " " << nframes << " " << npixels << endl; 00627 00628 vector<vector<double> > spot_intensity; //[spot][pixel] 00629 for(unsigned int i=0; i < nspots; i++) 00630 spot_intensity.push_back(compute_spot_intensity(pixels, spots[i])); 00631 00632 vector<tuple<double, Vector<4>, Matrix<4> > > spot_hess_etc = compute_spot_intensity_hessian(pixels, spots[spot]); 00633 00634 GibbsSampler sampler(pixel_intensities, spot_intensity, spots, A, pi, variance, sample_iterations); 00635 00636 //Compute derivative of log probability, summed (ie integrated) 00637 Matrix<4> sum_hess1 = Zeros(spots.size()); 00638 Matrix<4> sum_hess2 = Zeros(spots.size()); 00639 Vector<4> sum_deriv = Zeros(spots.size()); 00640 00641 for(int sample=0; sample < samples; sample++) 00642 { 00643 sampler.next(rng); 00644 00645 //Compute d log P(data | x, model) / d model, for a given sample 00646 //And the hessian 00647 Matrix<4> hess = Zeros(spots.size()); 00648 Vector<4> deriv = Zeros(spots.size()); 00649 for(unsigned int frame=0; frame < nframes; frame++) 00650 { 00651 for(unsigned int pixel=0; pixel < npixels; pixel++) 00652 { 00653 double e = pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel]; 00654 //Build up the derivative 00655 if(sampler.sample()[spot][frame] == 0) 00656 { 00657 hess += e * get<2>(spot_hess_etc[pixel]) - get<1>(spot_hess_etc[pixel]).as_col() * get<1>(spot_hess_etc[pixel]).as_row(); 00658 deriv += e * get<1>(spot_hess_etc[pixel]); 00659 00660 } 00661 } 00662 } 00663 00664 hess[0][0] += hess_log_log_normal(spots[spot][0], mu_brightness, sigma_brightness); 00665 hess[1][1] += hess_log_log_normal(spots[spot][1], mu_blur, sigma_blur); 00666 sum_hess1 += hess; 00667 00668 deriv[0] += diff_log_log_normal(spots[spot][0], mu_brightness, sigma_brightness); 00669 deriv[1] += diff_log_log_normal(spots[spot][1], mu_blur, sigma_blur); 00670 00671 sum_hess2 += deriv.as_col() * deriv.as_row(); 00672 00673 sum_deriv = sum_deriv + deriv; 00674 } 00675 00676 sum_hess1 /= (samples * variance); 00677 sum_hess2 /= (samples * variance); 00678 sum_deriv /= (samples * variance); 00679 00680 00681 cout << sum_hess1 << endl; 00682 cout << sum_hess2 << endl; 00683 cout << sum_deriv.as_col() * sum_deriv.as_row() << endl; 00684 00685 cout << "......." << sum_deriv << endl; 00686 //Put in the prior 00687 //The derivative prior parts cancel out. 00688 //Rather sensibly this means that the second derivatives can be 00689 //computed without reference to the prior, and then the prior can 00690 //be added in later, i.e.: P(data, parameters) = P(data|parameter) P(parameters) 00691 00692 //The second derivatives have been constructed to be diagonal 00693 DiagonalMatrix<4> hess_prior = Zeros(spots.size()); 00694 00695 cout << "sum of parts = \n" << sum_hess1 + sum_hess2 - sum_deriv.as_col() * sum_deriv.as_row() << endl; 00696 //TooN cannot currently add DiagonalMatrix to Matrix!! 00697 //sum_hess.diagonal_slice() += hess_prior.diagonal_slice(); 00698 00699 cout << "++++Done Computing pure MCMC hessian\n"; 00700 return sum_hess1 + sum_hess2 - sum_deriv.as_col() * sum_deriv.as_row(); 00701 }
MT19937& DataForMCMC::get_rng | ( | ) | const [inline, inherited] |
Definition at line 186 of file multispot5.cc.
References DataForMCMC::rng.
Referenced by NegativeFreeEnergy::compute_with_mask(), and NegativeFreeEnergy::operator()().
00187 { 00188 return rng; 00189 }
const vector<ImageRef>& DataForMCMC::pixels [protected, inherited] |
Definition at line 176 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), hessian(), and NegativeFreeEnergy::operator()().
const vector<vector<double> >& DataForMCMC::pixel_intensities [protected, inherited] |
Definition at line 177 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), hessian(), and NegativeFreeEnergy::operator()().
const double DataForMCMC::mu_brightness [protected, inherited] |
Definition at line 178 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), hessian(), and NegativeFreeEnergy::operator()().
const double DataForMCMC::sigma_brightness [protected, inherited] |
Definition at line 178 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), hessian(), and NegativeFreeEnergy::operator()().
const double DataForMCMC::mu_blur [protected, inherited] |
Definition at line 178 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), hessian(), and NegativeFreeEnergy::operator()().
const double DataForMCMC::sigma_blur [protected, inherited] |
Definition at line 178 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), hessian(), and NegativeFreeEnergy::operator()().
const double DataForMCMC::variance [protected, inherited] |
Definition at line 179 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), hessian(), and NegativeFreeEnergy::operator()().
const int DataForMCMC::samples [protected, inherited] |
Definition at line 180 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), hessian(), and NegativeFreeEnergy::operator()().
const int DataForMCMC::sample_iterations [protected, inherited] |
Definition at line 180 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), hessian(), and NegativeFreeEnergy::operator()().
const Matrix<3> DataForMCMC::A [protected, inherited] |
Definition at line 181 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), hessian(), and NegativeFreeEnergy::operator()().
const Vector<3> DataForMCMC::pi [protected, inherited] |
Definition at line 182 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), hessian(), and NegativeFreeEnergy::operator()().
MT19937& DataForMCMC::rng [mutable, protected, inherited] |
Definition at line 183 of file multispot5.cc.
Referenced by DataForMCMC::get_rng(), and hessian().