FitSpots Class Reference
[Storm classes]

Mega class which actually does the meat of the spot fitting. More...

List of all members.

Public Member Functions

 FitSpots (const vector< Image< float > > &ims_, FitSpotsGraphics &graphics_, StateParameters &params, ofstream &save_spots_)
void optimize_each_spot_in_turn_for_several_passes ()
void try_modifying_model ()
void run ()

Private Attributes

const vector< Image< float > > & ims
FitSpotsGraphicsgraphics
const vector< ImageRef > pixels
vector< Vector< 4 > > spots
const int start_iteration
int start_pass
MT19937rng
const double variance
const double intensity_mu
const double intensity_sigma
const double blur_mu
const double blur_sigma
const double area_extra_radius
set< ImageRef > allowed_area
const int use_position_prior
const double position_prior
const double max_motion
const int sample_iterations
const int main_cg_max_iterations
const int main_samples
const int main_passes
const int outer_loop_iterations
const int add_remove_tries
const int add_remove_opt_samples
const int add_remove_opt_retries
const int add_remove_opt_hess_inner_samples
const int h_outer_samples
const int h_inner_samples
const int tsamples
const Image< float > ave
ofstream & save_spots
double time_gibbs
double time_cg
const bool scale_brightness_limit
const Vector< 4 > limit
const Matrix< 3 > A
const Vector< 3 > pi
vector< vector< double > > pixel_intensities
DataForMCMC data_for_t_mcmc
DataForMCMC data_for_h_mcmc
int iteration

Detailed Description

Mega class which actually does the meat of the spot fitting.

probably could be refactored a bit.

Definition at line 1346 of file multispot5.cc.


Constructor & Destructor Documentation

FitSpots::FitSpots ( const vector< Image< float > > &  ims_,
FitSpotsGraphics graphics_,
StateParameters params,
ofstream &  save_spots_ 
) [inline]

Definition at line 1433 of file multispot5.cc.

References assert_same_size(), ims, pixel_intensities, and pixels.

01434     :ims(ims_), graphics(graphics_),
01435      
01436      //Set up the main parameters
01437      pixels(params.pixels),
01438      spots(params.spots),
01439      start_iteration(params.iteration),
01440      start_pass(params.pass),
01441      rng(*params.rng),
01442 
01443 
01444      
01445      //Set up all the system parameters
01446      variance(1), // it should be
01447 
01448      //To scale the X axis of a log-normal distribution, only
01449      //the mu parameter needs to be changed...
01450      intensity_mu(GV3::get<double>("intensity.rel_mu", 0., -1)  + log(sqrt(variance))),
01451      intensity_sigma(GV3::get<double>("intensity.rel_sigma", 0., -1)),
01452      blur_mu(GV3::get<double>("blur.mu", 0., -1)),
01453      blur_sigma(GV3::get<double>("blur.sigma", 0., -1)),
01454  
01455      //Spot position prior
01456      area_extra_radius(GV3::get<double>("position.extra_radius", 0., -1)),
01457      allowed_area(dilate_mask(pixels, area_extra_radius)),
01458      use_position_prior(GV3::get<bool>("position.use_prior", true, -1)),
01459      position_prior(1.0 / allowed_area.size()),
01460         
01461      //General optimizing and sampling parameters
01462      max_motion(GV3::get<double>("cg.max_motion", 0., -1)),
01463      sample_iterations(GV3::get<int>("gibbs.mixing_iterations", 0, -1)),
01464      
01465      //Task specific optimizing and sampling parameters:
01466  
01467      //Main optimization loop
01468      main_cg_max_iterations(GV3::get<double>("main.cg.max_iterations", 0., -1)),
01469      main_samples(GV3::get<int>("main.gibbs.samples", 0, -1)),
01470      main_passes(GV3::get<int>("main.passes", 0, -1)),
01471      outer_loop_iterations(GV3::get<int>("main.total_iterations", 100000000, 1)),
01472      
01473      //Spot selection loop
01474      add_remove_tries(GV3::get<int>("add_remove.tries", 0, -1)),
01475      add_remove_opt_samples(GV3::get<int>("add_remove.optimizer.samples", 0, -1)),
01476      add_remove_opt_retries(GV3::get<int>("add_remove.optimizer.attempts", 0, -1)),
01477      add_remove_opt_hess_inner_samples(GV3::get<int>("add_remove.optimizer.hessian_inner_samples", 0, -1)),
01478      h_outer_samples(GV3::get<int>("add_remove.hessian.outer_samples", 0, -1)),
01479      h_inner_samples(GV3::get<int>("add_remove.hessian.inner_samples", 0, -1)),
01480      tsamples(GV3::get<int>("add_remove.thermo.samples", 0, -1)),
01481 
01482      ave(average_image(ims_)),
01483 
01484      save_spots(save_spots_),
01485      
01486      time_gibbs(0),
01487      time_cg(0),
01488     
01489      scale_brightness_limit(GV3::get<bool>("max_motion.use_brightness_std", 0, -1)),
01490      limit(makeVector(brightness_motion_limit(intensity_mu, intensity_sigma, scale_brightness_limit), 1, 1, 1)*max_motion),
01491 
01492      A(GV3::get<Matrix<3> >("A", Zeros, 1)),
01493      pi(GV3::get<Vector<3> >("pi", Zeros, 1)),
01494 
01495 
01496      data_for_t_mcmc(pixels, pixel_intensities, intensity_mu, intensity_sigma, blur_mu, blur_sigma, variance, tsamples, sample_iterations, A, pi, rng),
01497      data_for_h_mcmc(pixels, pixel_intensities, intensity_mu, intensity_sigma, blur_mu, blur_sigma, variance, h_outer_samples, sample_iterations, A, pi, rng)
01498 
01499     {
01500         assert_same_size(ims);
01501 
01502         //Pixel intensities for all images [frame][pixel]
01503         pixel_intensities.resize(ims.size(), vector<double>(pixels.size()));
01504         for(unsigned int frame=0; frame < ims.size(); frame++)
01505             for(unsigned int p=0; p < pixels.size(); p++)
01506                 pixel_intensities[frame][p] = ims[frame][pixels[p]];
01507     
01508     }


Member Function Documentation

void FitSpots::optimize_each_spot_in_turn_for_several_passes (  )  [inline]

Perform a complete iteration of the optimzation stage of the spot firrint algorithm.

Definition at line 1511 of file multispot5.cc.

References A, allowed_area, ave, blur_mu, blur_sigma, boundingbox(), SampledMultispot::compute_spot_intensity(), data_for_h_mcmc, data_for_t_mcmc, FitSpotsGraphics::draw_krap(), FitSpotsGraphics::draw_pixels(), get_spot_pixels(), graphics, FreeEnergyHessian::hessian(), intensity_mu, intensity_sigma, ConjugateGradientOnly< Size, Precision >::iterate(), iteration, limit, main_cg_max_iterations, main_passes, main_samples, ConjugateGradientOnly< Size, Precision >::max_iterations, pi, pixel_intensities, pixels, SampledMultispot::remove_spot(), rng, sample_iterations, sampled_background_spot_hessian_ffbs(), save_spots, scale_to_bytes(), SampledMultispot::sequence(), spots, spots_to_Vector(), start_pass, FitSpotsGraphics::swap(), time_cg, time_gibbs, use_position_prior, variance, MT19937::write(), and ConjugateGradientOnly< Size, Precision >::x.

Referenced by run().

01512     {
01513         //Precompute the intensities for all spot pixels
01514         vector<vector<double> > spot_intensities; //[spot][pixel]
01515         for(unsigned int i=0; i < spots.size(); i++)
01516             spot_intensities.push_back(compute_spot_intensity(pixels, spots[i]));
01517 
01518         //Which pixels does each spot have?
01519         vector<vector<int> > spot_pixels; //[spot][pixel]
01520         spot_pixels.resize(spots.size());
01521         for(unsigned int s=0; s < spots.size(); s++)
01522             get_spot_pixels(pixels, spots[s], spot_pixels[s]);
01523 
01524         
01525         //Optimize the model, N spots at a time.
01526         //
01527         for(int pass=start_pass; pass < main_passes; pass++)
01528         {
01529             save_spots << "Pass: " << pass << endl;
01530             rng.write(save_spots);
01531             save_spots << endl;
01532 
01533             start_pass=0; // This is only nonzero first time, since we might chekpoint midway through an iteration
01534             save_spots << "PASS" << pass << ": " << setprecision(20) << scientific << spots_to_Vector(spots) << endl;
01535             save_spots << "ENDCHECKPOINT" << endl << flush;
01536             //Sort spots according to pass%4
01537 
01538             //Sort the spots so that the optimization runs in sweeps
01539             //This heiristic seems to increase the rate of propagation of information
01540             //about spot positions.
01541 
01542             //Create a list of indices
01543             vector<int> index = sequence(spots.size());
01544             
01545             int passs = pass + iteration;
01546             //Sort the indices according to the position of the spot that they index
01547             if(passs%4 == 0)
01548                 sort(index.begin(), index.end(), IndexLexicographicPosition<less<double>, 2>(spots));
01549             else if(passs%4==1)
01550                 sort(index.begin(), index.end(), IndexLexicographicPosition<greater<double>, 2>(spots));
01551             else if(passs%4==1)
01552                 sort(index.begin(), index.end(), IndexLexicographicPosition<less<double>, 3>(spots));
01553             else
01554                 sort(index.begin(), index.end(), IndexLexicographicPosition<greater<double>, 3>(spots));
01555 
01556             //Reorder the spots and their intensities and their pixel lists
01557             {
01558                 vector<Vector<4> > tmp_spot(index.size());
01559                 vector<vector<double> > tmp_spot_intensities(index.size());
01560                 vector<vector<int> > tmp_spot_pixels(index.size());
01561                 for(unsigned int i=0; i < index.size(); i++)
01562                 {
01563                     tmp_spot[i] = spots[index[i]];
01564                     swap(tmp_spot_intensities[i], spot_intensities[index[i]]);
01565                     swap(tmp_spot_pixels[i], spot_pixels[i]);
01566                 }
01567 
01568                 swap(tmp_spot, spots);
01569                 swap(tmp_spot_intensities, spot_intensities);
01570                 swap(tmp_spot_pixels, spot_pixels);
01571             }
01572 
01573             //Sweep through and optimize each spot in turn
01574             for(int s=0; s < (int)spots.size(); s++)
01575             {
01576                 timer.reset();
01577                 //Get some samples with Gibbs sampling
01578                 vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
01579                 vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
01580 
01581                 GibbsSampler2 sampler(pixel_intensities, spot_intensities, spots, spot_pixels, A, pi, variance, sample_iterations);
01582                 for(int i=0; i < main_samples; i++)
01583                 {
01584                     sampler.next(rng);
01585                     sample_list.push_back(sampler.sample());
01586                     sample_intensities.push_back(sampler.sample_intensities());
01587                 }
01588 
01589                 //First, remove the spot from all the samples.
01590                 for(unsigned int i=0; i < sample_list.size(); i++)
01591                     remove_spot(sample_intensities[i], spot_intensities[s], sample_list[i][s]);
01592                 
01593                 //cout << timer.get_time() << endl;
01594                 time_gibbs += timer.reset();
01595 
01596                 //Package up all the data
01597                 SampledBackgroundData data(sample_intensities, pixel_intensities, pixels, 
01598                                            intensity_mu, intensity_sigma, blur_mu, blur_sigma, 
01599                                            A, pi, variance);
01600                 
01601                 //Derivative computer:
01602                 SpotNegProbabilityDiffWithSampledBackground compute_deriv(data);
01603                 
01604 
01605                 graphics.draw_pixels(pixels, 0, 0, 1, 1);
01606                 graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), s);
01607                 graphics.swap();
01608 
01609                 //Optimize spot "s"
01610                 //Optimize with the derivatives only since the actual probability
01611                 //is much harder to compute
01612                 ConjugateGradientOnly<4> cg(spots[s], compute_deriv, limit);
01613 
01614 
01615                 cg.max_iterations = main_cg_max_iterations;
01616 
01617 
01618                 #if 0
01619                     cout << setprecision(10);
01620                     cout << spots_to_Vector(spots) << endl;
01621                     Matrix<4> hess, hess_errors;
01622                     cout << "Hello, my name is Inigo Montoya\n";
01623                     /*for(int i=0; i < 4; i++)
01624                     {
01625                         Matrix<4, 2> d = numerical_gradient_with_errors(NthDeriv(compute_deriv, i), cg.x);
01626                         hess[i] = d.T()[0];
01627                         hess_errors[i] = d.T()[1];
01628                     }
01629                     */
01630                     //cout << "Errors:\n" << hess_errors << endl;
01631                     //cout << "NHess:\n" << hess<< endl;
01632                     Matrix<4> rhess =  -sampled_background_spot_hessian(cg.x, data);
01633                     cout << "Hess:\n" << rhess << endl;
01634                     //cout << "Err:\n" << hess - rhess << endl;
01635 
01636                     //Vector<4> grad = compute_deriv(cg.x);
01637 
01638                     //Matrix<4> e = hess - rhess;
01639 
01640                     //for(int i=0; i < 4; i++)
01641                     //  for(int j=0; j < 4; j++)
01642                     //      e[i][j] /= hess_errors[i][j];
01643 
01644                     //cout << "Err:\n" << e << endl;
01645                     cout << "Deriv:" <<  compute_deriv(cg.x) << endl;
01646                     //cout << "Full:\n" << sampled_background_spot_hessian2(cg.x, data) - grad.as_col()*grad.as_row() << endl;
01647 
01648                     FreeEnergyHessian hesscomputer(data_for_h_mcmc);
01649 
01650                     Matrix<4> nhess = hesscomputer.hessian(spots, 0);
01651                     cout << "NHess:\n" << nhess << endl;
01652 
01653                     cout << "Turbo-N Hess:\n" << sampled_background_spot_hessian_ffbs(cg.x, data, 10000) << endl;
01654 
01655                     cout << "TI energy: " << NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(spots)) << endl;
01656                     cout << "FA energy: " <<  SpotProbabilityWithSampledBackgroundFAKE(data)(cg.x) << endl;
01657 
01658 
01659 
01660                     //cout << "Numerical hessian from FD:\n" << numerical_hessian(SpotProbabilityWithSampledBackgroundFAKE(data), cg.x) << endl;
01661                     exit(0);
01662                 #endif
01663                 //cout << "Starting opt... " << cg.x << endl;
01664                 while(cg.iterate(compute_deriv))
01665                 {
01666                     graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), s, cg.x);
01667                     graphics.draw_pixels(pixels, 0, 0, 1, .2);
01668                     graphics.swap();
01669                     //cout << cg.x << endl;
01670                 }
01671 
01672                 //Update to use the result of the optimization
01673                 spots[s] = cg.x;
01674                 //cout << "End: " << cg.x << endl;
01675                 
01676                 graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), -1);
01677                 graphics.swap();
01678                     
01679                 //Recompute the new spot intensity, since the spot has changed
01680                 spot_intensities[s] = compute_spot_intensity(pixels, spots[s]);
01681 
01682                 //Recompute which are the useful pixels
01683                 get_spot_pixels(pixels, spots[s], spot_pixels[s]);
01684                 
01685                 //Is the spot within the allowed area, i.e. is it's prior 0?
01686                 //The prior is sero only if it we are using it and we're in an invalid area
01687 
01688                 ImageRef quantized_spot_position = ir_rounded(spots[s].slice<2,2>());
01689                 bool zero_prior = use_position_prior && (allowed_area.count(quantized_spot_position)==0);
01690 
01691                 //Check to see if spot has been ejected. If spot_pixels is empty, then it has certainly been ejected.
01692                 if(spot_pixels[s].empty() || zero_prior)
01693                 {
01694                     //Spot has been ejected. Erase it.
01695                     cout  << " Erasing ejected spot: " << spot_pixels[s].empty() << " " << zero_prior << endl;
01696                     cout << spots[s] << endl;
01697                     //GUI_Pause(1);
01698 
01699                     spot_intensities.erase(spot_intensities.begin() + s);
01700                     spot_pixels.erase(spot_pixels.begin() + s);
01701                     spots.erase(spots.begin() + s);
01702                     s--;
01703                     //exit(0);
01704                 }
01705                 
01706                 //cout << timer.get_time() << endl;
01707                 time_cg += timer.reset();
01708 
01709                 //cout << "Times: " << time_gibbs << " " << time_cg << endl;
01710                 //save_spots << "INTERMEDIATE: " << setprecision(15) << scientific << spots_to_Vector(spots) << endl;
01711             }
01712         }
01713     }

void FitSpots::try_modifying_model (  )  [inline]

Perform a complete iteration of the model size modification stage of the spot fitting algorithm.

Definition at line 1716 of file multispot5.cc.

References A, add_remove_opt_hess_inner_samples, add_remove_opt_retries, add_remove_opt_samples, add_remove_tries, SampledMultispot::add_spot(), allowed_area, ave, blur_mu, blur_sigma, boundingbox(), SampledMultispot::compute_spot_intensity(), data_for_t_mcmc, FitSpotsGraphics::draw_krap(), get_spot_pixels(), graphics, h_inner_samples, h_outer_samples, ConjugateGradientOnly< Size, Precision >::init(), intensity_mu, intensity_sigma, ConjugateGradientOnly< Size, Precision >::iterate(), limit, ln(), log_normal_mode(), pi, pixel_intensities, pixels, position_prior, SampledMultispot::remove_spot(), rng, sample_iterations, sampled_background_spot_hessian_ffbs(), scale_to_bytes(), spots, spots_to_Vector(), FitSpotsGraphics::swap(), use_position_prior, variance, and ConjugateGradientOnly< Size, Precision >::x.

Referenced by run().

01717     {
01718 
01719         for(int i=0; i < add_remove_tries; i++)
01720         {
01721             cout << endl << endl << "Modifying the model" << endl << "======================\n";
01722             cout << "Hello\n";
01723             bool add_spot = (rng() > 0.5) || (spots.size() == 1);
01724             cout << "World\n";
01725 
01726             vector<Vector<4> > model_1, model_2;
01727 
01728             
01729             int original; //What is the original model? Model 1 or Model 2?
01730 
01731             if(add_spot)
01732             {
01733                 model_1 = spots;
01734                 model_2 = model_1;
01735                 
01736                 //Pick a pixel within the thresholded ones as a starting point
01737                 int r;
01738                 do
01739                 {
01740                     r = (int)(rng() * pixels.size());
01741                     //cout << r << " " << log_ratios[pixels[r]] << " " << pixels[r] << " " << threshold << endl;
01742                 }
01743                 while(0);
01744 
01745                 //This version does not (yet?) suppotrt thresholding on log_ratios
01746                 //for placing spots, since the purpose of log_ratios has diminished.
01747                 //while(log_ratios[pixels[r]] < threshold);
01748 
01749 
01750                 model_2.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), 
01751                                              log_normal_mode(blur_mu, blur_sigma), 
01752                                              pixels[r].x + rng()-.5, pixels[r].y + rng() - .5));
01753                 cout << "Adding a spot\n";
01754 
01755                 original = 1;
01756             }
01757             else
01758             {   
01759                 //Pick a random point to optimize/remove
01760                 int a_random_spot = static_cast<int>(rng() * spots.size());
01761                 model_1 = spots;
01762                 swap(model_1[model_1.size()-1], model_1[a_random_spot]);
01763 
01764                 model_2 = model_1;
01765 
01766                 model_1.pop_back();
01767                 cout << "Removing a spot\n";
01768                 original = 2;
01769             }
01770             
01771             //The mobile spot is always the last spot of model_2
01772             const int spot = model_2.size() - 1;
01773 
01774             cout << "Original model: " << original << endl;
01775 
01776             //Precompute the intensities for all spot pixels
01777             //Model 2 is always one longer than model 1 and 
01778             //differs only on the extra element
01779             vector<vector<double> > model2_spot_intensities; //[spot][pixel]
01780             for(unsigned int i=0; i < model_2.size(); i++)
01781                 model2_spot_intensities.push_back(compute_spot_intensity(pixels, model_2[i]));
01782 
01783             //Which pixels does each spot have?
01784             vector<vector<int> > model2_spot_pixels(model_2.size()); //[spot][pixel]
01785             for(unsigned int s=0; s < model_2.size(); s++)
01786                 get_spot_pixels(pixels, model_2[s], model2_spot_pixels[s]);
01787 
01788             //Optimize spot:
01789             {
01790                 cout << "Optimizing spot for model selection\n";
01791 
01792 
01793                 //Get some samples with Gibbs sampling
01794                 vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
01795                 vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
01796 
01797                 GibbsSampler2 sampler(pixel_intensities, model2_spot_intensities, model_2, model2_spot_pixels, A, pi, variance, sample_iterations);
01798                 for(int i=0; i < add_remove_opt_samples; i++)
01799                 {
01800                     sampler.next(rng);
01801                     sample_list.push_back(sampler.sample());
01802                     sample_intensities.push_back(sampler.sample_intensities());
01803                 }
01804 
01805                 //First, remove the spot from all the samples.
01806                 for(unsigned int i=0; i < sample_list.size(); i++)
01807                     remove_spot(sample_intensities[i], model2_spot_intensities[spot], sample_list[i][spot]);
01808 
01809                 //Package up all the data
01810                 SampledBackgroundData data(sample_intensities, pixel_intensities, pixels, 
01811                                            intensity_mu, intensity_sigma, blur_mu, blur_sigma, 
01812                                            A, pi, variance);
01813                 
01814                 //Derivative computer:
01815                 SpotNegProbabilityDiffWithSampledBackground compute_deriv(data);
01816                 
01817                 graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), spot);
01818                 graphics.swap();
01819 
01820                 //Optimize spot "s"
01821                 //Optimize with the derivatives only since the actual probability
01822                 //is much harder to compute
01823                 ConjugateGradientOnly<4> cg(model_2[spot], compute_deriv, limit);
01824 
01825 
01826                 for(int attempt=0; attempt < add_remove_opt_retries; attempt++)
01827                 {
01828                     cout << "Attempt " << attempt << " of " << add_remove_opt_retries << endl;
01829 
01830                     //Optimize with conjugate gradient
01831                     while(cg.iterate(compute_deriv))
01832                     {
01833                         graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), spot, cg.x);
01834                         graphics.swap();
01835                     }
01836                     
01837                     //Check for being at a saddle point (no point checking on the last try)
01838                     //All eigenvectors should be negative at a maximum.
01839                     //WTF: is this a bug? WTF?
01840                     //It was this:
01841                     //if(attempt < add_remove_tries - 1)
01842                     if(attempt < add_remove_opt_retries - 1)
01843                     {
01844                         Matrix<4> hessian = sampled_background_spot_hessian_ffbs(cg.x, data, add_remove_opt_hess_inner_samples, rng);
01845                         SymEigen<4> hess_decomp(hessian);
01846 
01847                         //cout << "What the fuck:" << endl << hessian << endl << hessian3<< endl << hessian2 << endl;
01848                         
01849                         cout << "Eigenvalues are: " << hess_decomp.get_evalues() << endl;
01850 
01851                         if(hess_decomp.is_negdef())
01852                             break;
01853                         else
01854                         {
01855                             //Restart in the direction of the best uphill part
01856                             cg.init(cg.x + 0.1 * hess_decomp.get_evectors()[3], (hess_decomp.get_evectors()[3]));
01857 
01858 
01859                             cout << "Grad = " << compute_deriv(cg.x) << endl;
01860                             for(int i=0; i < 4; i++)
01861                             {
01862                                 cout << "Direction: " << i << endl;
01863                                 cout << unit(compute_deriv(cg.x + 0.1*hess_decomp.get_evectors()[i])) * hess_decomp.get_evectors()[i] << endl;
01864                             }
01865 
01866                             for(int i=0; i < 4; i++)
01867                             {
01868                                 cout << "Direction: " << i << endl;
01869                                 Vector<4> d = Zeros;
01870                                 d[i] = 1;
01871                                 cout << unit(compute_deriv(cg.x + d)) * hess_decomp.get_evectors()[i] << endl;
01872                             }
01873                         }
01874                     }
01875                 }
01876 
01877 
01878                 //Update to use the result of the optimization
01879                 model_2[spot] = cg.x;
01880                 
01881                 graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), -1);
01882                 graphics.swap();
01883     
01884 
01885                 //Update cached data based on the new spot position
01886                 model2_spot_intensities[spot] = compute_spot_intensity(pixels, model_2[spot]);
01887                 get_spot_pixels(pixels, model_2[spot], model2_spot_pixels[spot]);
01888 
01889                 cout << "Done optimizing for model selection\n";
01890             }
01891 
01892 
01893             //Which model to keep?
01894             int keep=original;
01895         
01896             //Compute position prior (and we might be able to reject it really quickly here)
01897             bool zero_prior = use_position_prior && (allowed_area.count(ir_rounded(model_2[spot].slice<2,2>()))==0); 
01898             
01899             if(zero_prior)
01900             {
01901                 //Model 2 went bad, so we clearly keep model 1
01902                 keep = 1;
01903             }
01904             else
01905             {
01906                 
01907                 //The position prior is independent
01908                 //Compute the difference  model2 - model1
01909                 //This is only valid if model2 is in the valid region
01910                 double position_log_prior_model2_minus_model1;
01911                 if(use_position_prior)
01912                     position_log_prior_model2_minus_model1 = (model_2.size() - model_1.size()) * ln(position_prior);
01913                 else
01914                     position_log_prior_model2_minus_model1 = 0;
01915 
01916 
01917                 //Integrate model_2
01918                 //First compute the Hessian since this might go wrong.
01919 
01920                 //FreeEnergyHessian hesscomputer(data_for_h_mcmc);
01921                 Matrix<4> hess;// = hesscomputer.hessian(model_2, spot);
01922 
01923                 //Use turbohess here since it is much faster, as the backwards sampling step is fast
01924                 //We expect this hessian to be really quite different, if the spot has moved from
01925                 //a long way away, since the sampling will have changed dramatically
01926                 {
01927                     //Get some samples with Gibbs sampling
01928                     vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
01929                     vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
01930 
01931                     GibbsSampler sampler(pixel_intensities, model2_spot_intensities, model_2, A, pi, variance, sample_iterations);
01932                     for(int i=0; i < h_outer_samples; i++)
01933                     {
01934                         sampler.next(rng);
01935                         sample_list.push_back(sampler.sample());
01936                         sample_intensities.push_back(sampler.sample_intensities());
01937                     }
01938                     
01939                     //First, remove the spot from all the samples.
01940                     for(unsigned int i=0; i < sample_list.size(); i++)
01941                         remove_spot(sample_intensities[i], model2_spot_intensities[spot], sample_list[i][spot]);
01942 
01943                     //Package up all the data
01944                     SampledBackgroundData data(sample_intensities, pixel_intensities, pixels, 
01945                                                intensity_mu, intensity_sigma, blur_mu, blur_sigma, 
01946                                                A, pi, variance);
01947 
01948                     hess = sampled_background_spot_hessian_ffbs(model_2[spot], data, h_inner_samples, rng);
01949                 }
01950 
01951 
01952                 double det = determinant(-hess / (M_PI*2));
01953                 SymEigen<4> hess_decomp(-hess);
01954                 cout << "Hessien Eigenvalues are: " << hess_decomp.get_evalues() << endl;
01955                 const double smallest_evalue = 1e-6;
01956 
01957                 //If the determinant is negative, then we are still at a saddle point
01958                 //despite the attempts above, so abandon the attempt.
01959                 if(hess_decomp.get_evalues()[0] >  smallest_evalue)
01960                 {
01961 
01962                     //Compute the free energy of model 2 at the MLE estimate
01963                     cout << "Model 2:\n";
01964             //      double model_2_energy = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_2));
01965                     double model_2_energy = -NegativeFreeEnergy(data_for_t_mcmc).compute_with_mask(spots_to_Vector(model_2), model2_spot_pixels);
01966                     cout << "Energy: " << model_2_energy << endl;
01967                 
01968                     //Combine the MLE energy and Hessian using Laplace's approximation
01969                     double model_2_occam =  -0.5 * log(det);
01970                     double model_2_prob = model_2_energy + model_2_occam + position_log_prior_model2_minus_model1;
01971 
01972                     cout << "Occam: " << model_2_occam << endl;
01973                     cout << "Position: " << position_log_prior_model2_minus_model1 << endl;
01974                     cout << "Prob: " << model_2_prob << endl;
01975 
01976                     //Integrate model_1
01977                     //It has no parameters, in this formulation.
01978                     //double model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_1));
01979                     
01980                     //Note that model_1 always has one fewer spots, and the last spot is always the one
01981                     //missing, so we can make the correct mask very easily:
01982                     model2_spot_pixels.pop_back();
01983                     double model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc).compute_with_mask(spots_to_Vector(model_1), model2_spot_pixels);
01984                     cout << "Prob: " << model_1_prob << endl;
01985                     //model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_1));
01986 
01987                     if(model_2_prob  >  model_1_prob)
01988                         keep=2;
01989                     else
01990                         keep=1;
01991 
01992                     cout << "Models evaluated\n";
01993                 }
01994                 else
01995                 {
01996                     cout << "Determinant has bad eigenvalues!\n";
01997                     keep = original;
01998                     cout << hess_decomp.get_evalues() << endl;
01999                 }
02000             }
02001 
02002             if(keep == 2)
02003             {
02004                 spots = model_2;
02005                 cout << "Keeping model 2\n";
02006 
02007             }
02008             else
02009             {
02010                 spots = model_1;
02011                 cout << "Keeping model 1\n";
02012             }
02013 
02014             if(original != keep)
02015             {
02016                 cout << "Model changed!\n";
02017                 //break;
02018             }
02019         }
02020     }

void FitSpots::run (  )  [inline]

Run the complete optimization algorithm.

Definition at line 2023 of file multispot5.cc.

References graphics, ims, FitSpotsGraphics::init(), iteration, limit, optimize_each_spot_in_turn_for_several_passes(), outer_loop_iterations, pixels, save_spots, spots, spots_to_Vector(), start_iteration, and try_modifying_model().

Referenced by fit_spots_new().

02024     {
02025         graphics.init(ims[0].size());
02026         save_spots << "LOGVERSION 2" << endl;
02027         save_spots << "PIXELS";
02028         for(unsigned int i=0; i < pixels.size(); i++)
02029             save_spots << " " << pixels[i].x << " " << pixels[i].y;
02030         save_spots << endl;
02031 
02032         save_spots << "BEGINGVARLIST" << endl;
02033         GV3::print_var_list(save_spots, "", 1);
02034         save_spots << "ENDGVARLIST" << endl;
02035 
02036         //TODO all GVARS are set, so dump out gvars.
02037 
02038         cout << "Limit vector: " << limit << endl;
02039 
02040         for(iteration=start_iteration; iteration < outer_loop_iterations ;iteration++)
02041         {
02042             save_spots << "Iteration: " << iteration << endl;
02043             save_spots << "MAIN: " << setprecision(20) << scientific << spots_to_Vector(spots) << endl;
02044 
02045             cout << endl << endl << "----------------------" << endl << "Optimizing:\n";
02046             cout << spots.size() << endl;
02047             
02048 
02049             optimize_each_spot_in_turn_for_several_passes();
02050 
02051             //spot_intensities is be correct here!
02052             try_modifying_model();
02053         }
02054         save_spots << "FINAL: " << setprecision(15) << scientific << spots_to_Vector(spots) << endl;
02055     }


Member Data Documentation

const vector<Image<float> >& FitSpots::ims [private]

Input data.

Definition at line 1348 of file multispot5.cc.

Referenced by FitSpots(), and run().

Graphics class.

Definition at line 1349 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), run(), and try_modifying_model().

const vector<ImageRef> FitSpots::pixels [private]

Area in which to perform model fitting.

Definition at line 1350 of file multispot5.cc.

Referenced by FitSpots(), optimize_each_spot_in_turn_for_several_passes(), run(), and try_modifying_model().

vector<Vector<4> > FitSpots::spots [private]

State in terms of current spot positions.

Definition at line 1353 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), run(), and try_modifying_model().

const int FitSpots::start_iteration [private]

Starting iteration number (for restarting from checkpoint).

Definition at line 1356 of file multispot5.cc.

Referenced by run().

int FitSpots::start_pass [private]

Starting pass (for restarting from checkpoint).

Definition at line 1357 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes().

MT19937& FitSpots::rng [private]

Random numbewr generator.

Definition at line 1359 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

const double FitSpots::variance [private]

Variance of noise in data. Must be 1.

Definition at line 1361 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

const double FitSpots::intensity_mu [private]

Prior for spot intensity.

Definition at line 1362 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

const double FitSpots::intensity_sigma [private]

Prior for spot intensity.

Definition at line 1363 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

const double FitSpots::blur_mu [private]

Prior for spot shape.

Definition at line 1364 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

const double FitSpots::blur_sigma [private]

Prior for spt shape.

Definition at line 1365 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

const double FitSpots::area_extra_radius [private]

Extra size beyone marked region in which spots are allowed to exist.

Definition at line 1374 of file multispot5.cc.

set<ImageRef> FitSpots::allowed_area [private]

Total allowed region, pixels dilated by area_extra_radius.

Definition at line 1375 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

const int FitSpots::use_position_prior [private]

Should a proper prior over position be uesd? A clue: yes.

Definition at line 1376 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

const double FitSpots::position_prior [private]

Value for the posision prior, i.e. reciprocal of area.

Definition at line 1377 of file multispot5.cc.

Referenced by try_modifying_model().

const double FitSpots::max_motion [private]

Maximum motion on any axes for optimization. See ConjugateGradientOnly.

Definition at line 1380 of file multispot5.cc.

const int FitSpots::sample_iterations [private]

Number of mixing samples to use in Gibbs sampling.

Definition at line 1381 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

const int FitSpots::main_cg_max_iterations [private]

Maximum iterations allowed for ConjugateGradientOnly in main optimization loop.

Definition at line 1386 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes().

const int FitSpots::main_samples [private]

Number of samples to use in main loop.

Definition at line 1387 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes().

const int FitSpots::main_passes [private]

Number of passes to perform per iteration in main loop.

Definition at line 1388 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes().

const int FitSpots::outer_loop_iterations [private]

Total number of iterations to perform.

Definition at line 1389 of file multispot5.cc.

Referenced by run().

const int FitSpots::add_remove_tries [private]

Number of attemts to add/remove a spot.

Definition at line 1392 of file multispot5.cc.

Referenced by try_modifying_model().

const int FitSpots::add_remove_opt_samples [private]

Number of samples to use in model modification phase.

Definition at line 1393 of file multispot5.cc.

Referenced by try_modifying_model().

const int FitSpots::add_remove_opt_retries [private]

Number of attempts restarting the optimization to escape saddle points.

Definition at line 1394 of file multispot5.cc.

Referenced by try_modifying_model().

Number of extra FFBS samples to use for computing Hessian when testing for convergence to non-saddle point.

Definition at line 1395 of file multispot5.cc.

Referenced by try_modifying_model().

const int FitSpots::h_outer_samples [private]

Number of samples used for computing Hessian as part of Laplace's approximation.

Definition at line 1396 of file multispot5.cc.

Referenced by try_modifying_model().

const int FitSpots::h_inner_samples [private]

Number of additional FFBS samples to use for computing Hessian as part of Laplace's approximation.

Definition at line 1397 of file multispot5.cc.

Referenced by try_modifying_model().

const int FitSpots::tsamples [private]

Number of samples to use in thermodynamic integration.

Definition at line 1398 of file multispot5.cc.

const Image<float> FitSpots::ave [private]

Average of input data: used for.

Definition at line 1401 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

ofstream& FitSpots::save_spots [private]

Output stream for log file.

Definition at line 1405 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and run().

double FitSpots::time_gibbs [private]

Benchmarking data.

Definition at line 1408 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes().

double FitSpots::time_cg [private]

Benchmarking data.

Definition at line 1409 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes().

const bool FitSpots::scale_brightness_limit [private]

Motion limit for ConjugateGradientOnly The x distance, y distance and spot size are all approximately the same scale which is of order 1.

The brigntness is not. By default, the brightness limit is also 1. This flag controls whether is should be set to the standard deviation of the brightness prior distribution. This setting will put the motion limit on the same scale as the other three parameters.

Definition at line 1417 of file multispot5.cc.

const Vector<4> FitSpots::limit [private]

Limit vector for ConjugateGradientOnly.

Definition at line 1418 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), run(), and try_modifying_model().

const Matrix<3> FitSpots::A [private]

Transition matrix.

Definition at line 1420 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

const Vector<3> FitSpots::pi [private]

Initial probabilities.

Definition at line 1421 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

vector<vector<double> > FitSpots::pixel_intensities [private]

Pixel intensities for all images [frame][pixel].

Definition at line 1424 of file multispot5.cc.

Referenced by FitSpots(), optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

Aggergated data for thermodynamic integration.

Definition at line 1426 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and try_modifying_model().

Aggergated data for finding hessian.

Definition at line 1427 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes().

int FitSpots::iteration [private]

Iteration number.

Definition at line 1429 of file multispot5.cc.

Referenced by optimize_each_spot_in_turn_for_several_passes(), and run().


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