|
ThreeB 1.1
|
Mega class which actually does the meat of the spot fitting. More...
Public Member Functions | |
| FitSpots (const vector< Image< float > > &ims_, FitSpotsGraphics &graphics_, UserInterfaceCallback &ui_, StateParameters ¶ms, 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 |
| FitSpotsGraphics & | graphics |
| UserInterfaceCallback & | ui |
| const vector< ImageRef > | pixels |
| vector< Vector< 4 > > | spots |
| const int | start_iteration |
| int | start_pass |
| MT19937 & | rng |
| 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 |
Mega class which actually does the meat of the spot fitting.
probably could be refactored a bit.
Definition at line 1348 of file multispot5.cc.
| FitSpots::FitSpots | ( | const vector< Image< float > > & | ims_, |
| FitSpotsGraphics & | graphics_, | ||
| UserInterfaceCallback & | ui_, | ||
| StateParameters & | params, | ||
| ofstream & | save_spots_ | ||
| ) | [inline] |
Definition at line 1438 of file multispot5.cc.
References assert_same_size().
:ims(ims_), graphics(graphics_),ui(ui_), //Set up the main parameters pixels(params.pixels), spots(params.spots), start_iteration(params.iteration), start_pass(params.pass), rng(*params.rng), //Set up all the system parameters variance(1), // it should be //To scale the X axis of a log-normal distribution, only //the mu parameter needs to be changed... intensity_mu(GV3::get<double>("intensity.rel_mu", 0., -1) + log(sqrt(variance))), intensity_sigma(GV3::get<double>("intensity.rel_sigma", 0., -1)), blur_mu(GV3::get<double>("blur.mu", 0., -1)), blur_sigma(GV3::get<double>("blur.sigma", 0., -1)), //Spot position prior area_extra_radius(GV3::get<double>("position.extra_radius", 0., -1)), allowed_area(dilate_mask(pixels, area_extra_radius)), use_position_prior(GV3::get<bool>("position.use_prior", true, -1)), position_prior(1.0 / allowed_area.size()), //General optimizing and sampling parameters max_motion(GV3::get<double>("cg.max_motion", 0., -1)), sample_iterations(GV3::get<int>("gibbs.mixing_iterations", 0, -1)), //Task specific optimizing and sampling parameters: //Main optimization loop main_cg_max_iterations(GV3::get<double>("main.cg.max_iterations", 0., -1)), main_samples(GV3::get<int>("main.gibbs.samples", 0, -1)), main_passes(GV3::get<int>("main.passes", 0, -1)), outer_loop_iterations(GV3::get<int>("main.total_iterations", 100000000, 1)), //Spot selection loop add_remove_tries(GV3::get<int>("add_remove.tries", 0, -1)), add_remove_opt_samples(GV3::get<int>("add_remove.optimizer.samples", 0, -1)), add_remove_opt_retries(GV3::get<int>("add_remove.optimizer.attempts", 0, -1)), add_remove_opt_hess_inner_samples(GV3::get<int>("add_remove.optimizer.hessian_inner_samples", 0, -1)), h_outer_samples(GV3::get<int>("add_remove.hessian.outer_samples", 0, -1)), h_inner_samples(GV3::get<int>("add_remove.hessian.inner_samples", 0, -1)), tsamples(GV3::get<int>("add_remove.thermo.samples", 0, -1)), ave(average_image(ims_)), save_spots(save_spots_), time_gibbs(0), time_cg(0), scale_brightness_limit(GV3::get<bool>("max_motion.use_brightness_std", 0, -1)), limit(makeVector(brightness_motion_limit(intensity_mu, intensity_sigma, scale_brightness_limit), 1, 1, 1)*max_motion), A(GV3::get<Matrix<3> >("A", Zeros, 1)), pi(GV3::get<Vector<3> >("pi", Zeros, 1)), data_for_t_mcmc(pixels, pixel_intensities, intensity_mu, intensity_sigma, blur_mu, blur_sigma, variance, tsamples, sample_iterations, A, pi, rng), data_for_h_mcmc(pixels, pixel_intensities, intensity_mu, intensity_sigma, blur_mu, blur_sigma, variance, h_outer_samples, sample_iterations, A, pi, rng) { assert_same_size(ims); //Pixel intensities for all images [frame][pixel] pixel_intensities.resize(ims.size(), vector<double>(pixels.size())); for(unsigned int frame=0; frame < ims.size(); frame++) for(unsigned int p=0; p < pixels.size(); p++) pixel_intensities[frame][p] = ims[frame][pixels[p]]; }
| 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 1516 of file multispot5.cc.
References boundingbox(), SampledMultispot::compute_spot_intensity(), get_spot_pixels(), FreeEnergyHessian::hessian(), ConjugateGradientOnly< Size, Precision >::iterate(), ConjugateGradientOnly< Size, Precision >::max_iterations, SampledMultispot::GibbsSampler2::next(), SampledMultispot::remove_spot(), SampledMultispot::GibbsSampler2::sample(), SampledMultispot::GibbsSampler2::sample_intensities(), sampled_background_spot_hessian_ffbs(), scale_to_bytes(), SampledMultispot::sequence(), spots_to_Vector(), TIME, and ConjugateGradientOnly< Size, Precision >::x.
{
//Precompute the intensities for all spot pixels
vector<vector<double> > spot_intensities; //[spot][pixel]
for(unsigned int i=0; i < spots.size(); i++)
spot_intensities.push_back(compute_spot_intensity(pixels, spots[i]));
//Which pixels does each spot have?
vector<vector<int> > spot_pixels; //[spot][pixel]
spot_pixels.resize(spots.size());
for(unsigned int s=0; s < spots.size(); s++)
get_spot_pixels(pixels, spots[s], spot_pixels[s]);
//Optimize the model, N spots at a time.
//
for(int pass=start_pass; pass < main_passes; pass++)
{
save_spots << "Pass: " << pass << endl;
rng.write(save_spots);
save_spots << endl;
start_pass=0; // This is only nonzero first time, since we might chekpoint midway through an iteration
save_spots << "PASS" << pass << ": " << setprecision(20) << scientific << spots_to_Vector(spots) << endl;
save_spots << "ENDCHECKPOINT" << endl << flush;
ui.per_pass(iteration, pass, spots);
//Sort spots according to pass%4
//Sort the spots so that the optimization runs in sweeps
//This heiristic seems to increase the rate of propagation of information
//about spot positions.
//Create a list of indices
vector<int> index = sequence(spots.size());
int passs = pass + iteration;
//Sort the indices according to the position of the spot that they index
if(passs%4 == 0)
sort(index.begin(), index.end(), IndexLexicographicPosition<less<double>, 2>(spots));
else if(passs%4==1)
sort(index.begin(), index.end(), IndexLexicographicPosition<greater<double>, 2>(spots));
else if(passs%4==1)
sort(index.begin(), index.end(), IndexLexicographicPosition<less<double>, 3>(spots));
else
sort(index.begin(), index.end(), IndexLexicographicPosition<greater<double>, 3>(spots));
//Reorder the spots and their intensities and their pixel lists
{
vector<Vector<4> > tmp_spot(index.size());
vector<vector<double> > tmp_spot_intensities(index.size());
vector<vector<int> > tmp_spot_pixels(index.size());
for(unsigned int i=0; i < index.size(); i++)
{
tmp_spot[i] = spots[index[i]];
swap(tmp_spot_intensities[i], spot_intensities[index[i]]);
swap(tmp_spot_pixels[i], spot_pixels[i]);
}
swap(tmp_spot, spots);
swap(tmp_spot_intensities, spot_intensities);
swap(tmp_spot_pixels, spot_pixels);
}
//Sweep through and optimize each spot in turn
for(int s=0; s < (int)spots.size(); s++)
{
ui.per_spot(iteration, pass, s, spots.size());
ui.perhaps_stop();
TIME(timer.reset();)
//Get some samples with Gibbs sampling
vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
GibbsSampler2 sampler(pixel_intensities, spot_intensities, spots, spot_pixels, A, pi, variance, sample_iterations);
for(int i=0; i < main_samples; i++)
{
sampler.next(rng);
sample_list.push_back(sampler.sample());
sample_intensities.push_back(sampler.sample_intensities());
ui.perhaps_stop();
}
//First, remove the spot from all the samples.
for(unsigned int i=0; i < sample_list.size(); i++)
remove_spot(sample_intensities[i], spot_intensities[s], sample_list[i][s]);
//cout << timer.get_time() << endl;
TIME(time_gibbs += timer.reset();)
//Package up all the data
SampledBackgroundData data(sample_intensities, pixel_intensities, pixels,
intensity_mu, intensity_sigma, blur_mu, blur_sigma,
A, pi, variance);
//Derivative computer:
SpotNegProbabilityDiffWithSampledBackground compute_deriv(data);
graphics.draw_pixels(pixels, 0, 0, 1, 1);
graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), s);
graphics.swap();
//Optimize spot "s"
//Optimize with the derivatives only since the actual probability
//is much harder to compute
ConjugateGradientOnly<4> cg(spots[s], compute_deriv, limit);
cg.max_iterations = main_cg_max_iterations;
#if 0
cout << setprecision(10);
cout << spots_to_Vector(spots) << endl;
Matrix<4> hess, hess_errors;
cout << "Hello, my name is Inigo Montoya\n";
/*for(int i=0; i < 4; i++)
{
Matrix<4, 2> d = numerical_gradient_with_errors(NthDeriv(compute_deriv, i), cg.x);
hess[i] = d.T()[0];
hess_errors[i] = d.T()[1];
}
*/
//cout << "Errors:\n" << hess_errors << endl;
//cout << "NHess:\n" << hess<< endl;
Matrix<4> rhess = -sampled_background_spot_hessian(cg.x, data);
cout << "Hess:\n" << rhess << endl;
//cout << "Err:\n" << hess - rhess << endl;
//Vector<4> grad = compute_deriv(cg.x);
//Matrix<4> e = hess - rhess;
//for(int i=0; i < 4; i++)
// for(int j=0; j < 4; j++)
// e[i][j] /= hess_errors[i][j];
//cout << "Err:\n" << e << endl;
cout << "Deriv:" << compute_deriv(cg.x) << endl;
//cout << "Full:\n" << sampled_background_spot_hessian2(cg.x, data) - grad.as_col()*grad.as_row() << endl;
FreeEnergyHessian hesscomputer(data_for_h_mcmc);
Matrix<4> nhess = hesscomputer.hessian(spots, 0);
cout << "NHess:\n" << nhess << endl;
cout << "Turbo-N Hess:\n" << sampled_background_spot_hessian_ffbs(cg.x, data, 10000) << endl;
cout << "TI energy: " << NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(spots)) << endl;
cout << "FA energy: " << SpotProbabilityWithSampledBackgroundFAKE(data)(cg.x) << endl;
//cout << "Numerical hessian from FD:\n" << numerical_hessian(SpotProbabilityWithSampledBackgroundFAKE(data), cg.x) << endl;
exit(0);
#endif
//cout << "Starting opt... " << cg.x << endl;
while(cg.iterate(compute_deriv))
{
graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), s, cg.x);
graphics.draw_pixels(pixels, 0, 0, 1, .2);
graphics.swap();
ui.perhaps_stop();
//cout << cg.x << endl;
}
//Update to use the result of the optimization
spots[s] = cg.x;
//cout << "End: " << cg.x << endl;
graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), -1);
graphics.swap();
//Recompute the new spot intensity, since the spot has changed
spot_intensities[s] = compute_spot_intensity(pixels, spots[s]);
//Recompute which are the useful pixels
get_spot_pixels(pixels, spots[s], spot_pixels[s]);
//Is the spot within the allowed area, i.e. is it's prior 0?
//The prior is sero only if it we are using it and we're in an invalid area
ImageRef quantized_spot_position = ir_rounded(spots[s].slice<2,2>());
bool zero_prior = use_position_prior && (allowed_area.count(quantized_spot_position)==0);
//Check to see if spot has been ejected. If spot_pixels is empty, then it has certainly been ejected.
if(spot_pixels[s].empty() || zero_prior)
{
//Spot has been ejected. Erase it.
cout << " Erasing ejected spot: " << spot_pixels[s].empty() << " " << zero_prior << endl;
cout << spots[s] << endl;
//GUI_Pause(1);
spot_intensities.erase(spot_intensities.begin() + s);
spot_pixels.erase(spot_pixels.begin() + s);
spots.erase(spots.begin() + s);
s--;
//exit(0);
}
//cout << timer.get_time() << endl;
TIME(time_cg += timer.reset();)
//cout << "Times: " << time_gibbs << " " << time_cg << endl;
//save_spots << "INTERMEDIATE: " << setprecision(15) << scientific << spots_to_Vector(spots) << endl;
}
}
}
| void FitSpots::try_modifying_model | ( | ) | [inline] |
Perform a complete iteration of the model size modification stage of the spot fitting algorithm.
Definition at line 1729 of file multispot5.cc.
References SampledMultispot::add_spot(), boundingbox(), SampledMultispot::compute_spot_intensity(), NegativeFreeEnergy::compute_with_mask(), get_spot_pixels(), ConjugateGradientOnly< Size, Precision >::init(), ConjugateGradientOnly< Size, Precision >::iterate(), ln(), log_normal_mode(), SampledMultispot::GibbsSampler::next(), SampledMultispot::GibbsSampler2::next(), SampledMultispot::remove_spot(), SampledMultispot::GibbsSampler::sample(), SampledMultispot::GibbsSampler2::sample(), SampledMultispot::GibbsSampler::sample_intensities(), SampledMultispot::GibbsSampler2::sample_intensities(), sampled_background_spot_hessian_ffbs(), scale_to_bytes(), spots_to_Vector(), and ConjugateGradientOnly< Size, Precision >::x.
{
for(int i=0; i < add_remove_tries; i++)
{
ui.per_modification(iteration, i, add_remove_tries);
ui.perhaps_stop();
cout << endl << endl << "Modifying the model" << endl << "======================\n";
cout << "Hello\n";
bool add_spot = (rng() > 0.5) || (spots.size() == 1);
cout << "World\n";
vector<Vector<4> > model_1, model_2;
int original; //What is the original model? Model 1 or Model 2?
if(add_spot)
{
model_1 = spots;
model_2 = model_1;
//Pick a pixel within the thresholded ones as a starting point
int r;
do
{
r = (int)(rng() * pixels.size());
//cout << r << " " << log_ratios[pixels[r]] << " " << pixels[r] << " " << threshold << endl;
}
while(0);
//This version does not (yet?) suppotrt thresholding on log_ratios
//for placing spots, since the purpose of log_ratios has diminished.
//while(log_ratios[pixels[r]] < threshold);
model_2.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma),
log_normal_mode(blur_mu, blur_sigma),
pixels[r].x + rng()-.5, pixels[r].y + rng() - .5));
cout << "Adding a spot\n";
original = 1;
}
else
{
//Pick a random point to optimize/remove
int a_random_spot = static_cast<int>(rng() * spots.size());
model_1 = spots;
swap(model_1[model_1.size()-1], model_1[a_random_spot]);
model_2 = model_1;
model_1.pop_back();
cout << "Removing a spot\n";
original = 2;
}
//The mobile spot is always the last spot of model_2
const int spot = model_2.size() - 1;
cout << "Original model: " << original << endl;
//Precompute the intensities for all spot pixels
//Model 2 is always one longer than model 1 and
//differs only on the extra element
vector<vector<double> > model2_spot_intensities; //[spot][pixel]
for(unsigned int i=0; i < model_2.size(); i++)
model2_spot_intensities.push_back(compute_spot_intensity(pixels, model_2[i]));
//Which pixels does each spot have?
vector<vector<int> > model2_spot_pixels(model_2.size()); //[spot][pixel]
for(unsigned int s=0; s < model_2.size(); s++)
get_spot_pixels(pixels, model_2[s], model2_spot_pixels[s]);
//Optimize spot:
{
cout << "Optimizing spot for model selection\n";
//Get some samples with Gibbs sampling
vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
GibbsSampler2 sampler(pixel_intensities, model2_spot_intensities, model_2, model2_spot_pixels, A, pi, variance, sample_iterations);
for(int i=0; i < add_remove_opt_samples; i++)
{
sampler.next(rng);
sample_list.push_back(sampler.sample());
sample_intensities.push_back(sampler.sample_intensities());
ui.perhaps_stop();
}
//First, remove the spot from all the samples.
for(unsigned int i=0; i < sample_list.size(); i++)
remove_spot(sample_intensities[i], model2_spot_intensities[spot], sample_list[i][spot]);
//Package up all the data
SampledBackgroundData data(sample_intensities, pixel_intensities, pixels,
intensity_mu, intensity_sigma, blur_mu, blur_sigma,
A, pi, variance);
//Derivative computer:
SpotNegProbabilityDiffWithSampledBackground compute_deriv(data);
graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), spot);
graphics.swap();
//Optimize spot "s"
//Optimize with the derivatives only since the actual probability
//is much harder to compute
ConjugateGradientOnly<4> cg(model_2[spot], compute_deriv, limit);
for(int attempt=0; attempt < add_remove_opt_retries; attempt++)
{
cout << "Attempt " << attempt << " of " << add_remove_opt_retries << endl;
ui.perhaps_stop();
//Optimize with conjugate gradient
while(cg.iterate(compute_deriv))
{
ui.perhaps_stop();
graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), spot, cg.x);
graphics.swap();
}
//Check for being at a saddle point (no point checking on the last try)
//All eigenvectors should be negative at a maximum.
//WTF: is this a bug? WTF?
//It was this:
//if(attempt < add_remove_tries - 1)
if(attempt < add_remove_opt_retries - 1)
{
Matrix<4> hessian = sampled_background_spot_hessian_ffbs(cg.x, data, add_remove_opt_hess_inner_samples, rng);
SymEigen<4> hess_decomp(hessian);
//cout << "What the fuck:" << endl << hessian << endl << hessian3<< endl << hessian2 << endl;
cout << "Eigenvalues are: " << hess_decomp.get_evalues() << endl;
if(hess_decomp.is_negdef())
break;
else
{
//Restart in the direction of the best uphill part
cg.init(cg.x + 0.1 * hess_decomp.get_evectors()[3], (hess_decomp.get_evectors()[3]));
cout << "Grad = " << compute_deriv(cg.x) << endl;
for(int i=0; i < 4; i++)
{
cout << "Direction: " << i << endl;
cout << unit(compute_deriv(cg.x + 0.1*hess_decomp.get_evectors()[i])) * hess_decomp.get_evectors()[i] << endl;
}
for(int i=0; i < 4; i++)
{
cout << "Direction: " << i << endl;
Vector<4> d = Zeros;
d[i] = 1;
cout << unit(compute_deriv(cg.x + d)) * hess_decomp.get_evectors()[i] << endl;
}
}
}
}
//Update to use the result of the optimization
model_2[spot] = cg.x;
graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), -1);
graphics.swap();
//Update cached data based on the new spot position
model2_spot_intensities[spot] = compute_spot_intensity(pixels, model_2[spot]);
get_spot_pixels(pixels, model_2[spot], model2_spot_pixels[spot]);
cout << "Done optimizing for model selection\n";
}
//Which model to keep?
int keep=original;
//Compute position prior (and we might be able to reject it really quickly here)
bool zero_prior = use_position_prior && (allowed_area.count(ir_rounded(model_2[spot].slice<2,2>()))==0);
if(zero_prior)
{
//Model 2 went bad, so we clearly keep model 1
keep = 1;
}
else
{
//The position prior is independent
//Compute the difference model2 - model1
//This is only valid if model2 is in the valid region
double position_log_prior_model2_minus_model1;
if(use_position_prior)
position_log_prior_model2_minus_model1 = (model_2.size() - model_1.size()) * ln(position_prior);
else
position_log_prior_model2_minus_model1 = 0;
//Integrate model_2
//First compute the Hessian since this might go wrong.
//FreeEnergyHessian hesscomputer(data_for_h_mcmc);
Matrix<4> hess;// = hesscomputer.hessian(model_2, spot);
//Use turbohess here since it is much faster, as the backwards sampling step is fast
//We expect this hessian to be really quite different, if the spot has moved from
//a long way away, since the sampling will have changed dramatically
{
//Get some samples with Gibbs sampling
vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
GibbsSampler sampler(pixel_intensities, model2_spot_intensities, model_2, A, pi, variance, sample_iterations);
for(int i=0; i < h_outer_samples; i++)
{
ui.perhaps_stop();
sampler.next(rng);
sample_list.push_back(sampler.sample());
sample_intensities.push_back(sampler.sample_intensities());
}
//First, remove the spot from all the samples.
for(unsigned int i=0; i < sample_list.size(); i++)
remove_spot(sample_intensities[i], model2_spot_intensities[spot], sample_list[i][spot]);
//Package up all the data
SampledBackgroundData data(sample_intensities, pixel_intensities, pixels,
intensity_mu, intensity_sigma, blur_mu, blur_sigma,
A, pi, variance);
hess = sampled_background_spot_hessian_ffbs(model_2[spot], data, h_inner_samples, rng);
}
double det = determinant(-hess / (M_PI*2));
SymEigen<4> hess_decomp(-hess);
cout << "Hessien Eigenvalues are: " << hess_decomp.get_evalues() << endl;
const double smallest_evalue = 1e-6;
//If the determinant is negative, then we are still at a saddle point
//despite the attempts above, so abandon the attempt.
if(hess_decomp.get_evalues()[0] > smallest_evalue)
{
//Compute the free energy of model 2 at the MLE estimate
cout << "Model 2:\n";
// double model_2_energy = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_2));
double model_2_energy = -NegativeFreeEnergy(data_for_t_mcmc).compute_with_mask(spots_to_Vector(model_2), model2_spot_pixels);
cout << "Energy: " << model_2_energy << endl;
//Combine the MLE energy and Hessian using Laplace's approximation
double model_2_occam = -0.5 * log(det);
double model_2_prob = model_2_energy + model_2_occam + position_log_prior_model2_minus_model1;
cout << "Occam: " << model_2_occam << endl;
cout << "Position: " << position_log_prior_model2_minus_model1 << endl;
cout << "Prob: " << model_2_prob << endl;
//Integrate model_1
//It has no parameters, in this formulation.
//double model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_1));
//Note that model_1 always has one fewer spots, and the last spot is always the one
//missing, so we can make the correct mask very easily:
model2_spot_pixels.pop_back();
double model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc).compute_with_mask(spots_to_Vector(model_1), model2_spot_pixels);
cout << "Prob: " << model_1_prob << endl;
//model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_1));
if(model_2_prob > model_1_prob)
keep=2;
else
keep=1;
cout << "Models evaluated\n";
}
else
{
cout << "Determinant has bad eigenvalues!\n";
keep = original;
cout << hess_decomp.get_evalues() << endl;
}
}
if(keep == 2)
{
spots = model_2;
cout << "Keeping model 2\n";
}
else
{
spots = model_1;
cout << "Keeping model 1\n";
}
if(original != keep)
{
cout << "Model changed!\n";
//break;
}
}
}
| void FitSpots::run | ( | ) | [inline] |
Run the complete optimization algorithm.
Definition at line 2043 of file multispot5.cc.
References spots_to_Vector().
Referenced by fit_spots_new().
{
graphics.init(ims[0].size());
save_spots << "LOGVERSION 2" << endl;
save_spots << "PIXELS";
for(unsigned int i=0; i < pixels.size(); i++)
save_spots << " " << pixels[i].x << " " << pixels[i].y;
save_spots << endl;
save_spots << "BEGINGVARLIST" << endl;
GV3::print_var_list(save_spots, "", 1);
save_spots << "ENDGVARLIST" << endl;
//TODO all GVARS are set, so dump out gvars.
cout << "Limit vector: " << limit << endl;
for(iteration=start_iteration; iteration < outer_loop_iterations ;iteration++)
{
save_spots << "Iteration: " << iteration << " (" << iteration * main_passes << ")" << endl;
save_spots << "MAIN: " << setprecision(20) << scientific << spots_to_Vector(spots) << endl;
cout << endl << endl << "----------------------" << endl << "Optimizing:\n";
cout << spots.size() << endl;
optimize_each_spot_in_turn_for_several_passes();
//spot_intensities is be correct here!
try_modifying_model();
}
save_spots << "FINAL: " << setprecision(15) << scientific << spots_to_Vector(spots) << endl;
}
const vector<Image<float> >& FitSpots::ims [private] |
Input data.
Definition at line 1350 of file multispot5.cc.
FitSpotsGraphics& FitSpots::graphics [private] |
Graphics class.
Definition at line 1351 of file multispot5.cc.
UserInterfaceCallback& FitSpots::ui [private] |
Callbacks to provide user interface.
Definition at line 1352 of file multispot5.cc.
const vector<ImageRef> FitSpots::pixels [private] |
Area in which to perform model fitting.
Definition at line 1353 of file multispot5.cc.
vector<Vector<4> > FitSpots::spots [private] |
State in terms of current spot positions.
Definition at line 1356 of file multispot5.cc.
const int FitSpots::start_iteration [private] |
Starting iteration number (for restarting from checkpoint)
Definition at line 1359 of file multispot5.cc.
int FitSpots::start_pass [private] |
Starting pass (for restarting from checkpoint)
Definition at line 1360 of file multispot5.cc.
MT19937& FitSpots::rng [private] |
Random numbewr generator.
Definition at line 1362 of file multispot5.cc.
const double FitSpots::variance [private] |
Variance of noise in data. Must be 1.
Definition at line 1364 of file multispot5.cc.
const double FitSpots::intensity_mu [private] |
Prior for spot intensity.
Definition at line 1365 of file multispot5.cc.
const double FitSpots::intensity_sigma [private] |
Prior for spot intensity.
Definition at line 1366 of file multispot5.cc.
const double FitSpots::blur_mu [private] |
Prior for spot shape.
Definition at line 1367 of file multispot5.cc.
const double FitSpots::blur_sigma [private] |
Prior for spt shape.
Definition at line 1368 of file multispot5.cc.
const double FitSpots::area_extra_radius [private] |
Extra size beyone marked region in which spots are allowed to exist.
Definition at line 1377 of file multispot5.cc.
set<ImageRef> FitSpots::allowed_area [private] |
Total allowed region, pixels dilated by area_extra_radius.
Definition at line 1378 of file multispot5.cc.
const int FitSpots::use_position_prior [private] |
Should a proper prior over position be uesd? A clue: yes.
Definition at line 1379 of file multispot5.cc.
const double FitSpots::position_prior [private] |
Value for the posision prior, i.e. reciprocal of area.
Definition at line 1380 of file multispot5.cc.
const double FitSpots::max_motion [private] |
Maximum motion on any axes for optimization. See ConjugateGradientOnly.
Definition at line 1383 of file multispot5.cc.
const int FitSpots::sample_iterations [private] |
Number of mixing samples to use in Gibbs sampling.
Definition at line 1384 of file multispot5.cc.
const int FitSpots::main_cg_max_iterations [private] |
Maximum iterations allowed for ConjugateGradientOnly in main optimization loop.
Definition at line 1389 of file multispot5.cc.
const int FitSpots::main_samples [private] |
Number of samples to use in main loop.
Definition at line 1390 of file multispot5.cc.
const int FitSpots::main_passes [private] |
Number of passes to perform per iteration in main loop.
Definition at line 1391 of file multispot5.cc.
const int FitSpots::outer_loop_iterations [private] |
Total number of iterations to perform.
Definition at line 1392 of file multispot5.cc.
const int FitSpots::add_remove_tries [private] |
Number of attemts to add/remove a spot.
Definition at line 1395 of file multispot5.cc.
const int FitSpots::add_remove_opt_samples [private] |
Number of samples to use in model modification phase.
Definition at line 1396 of file multispot5.cc.
const int FitSpots::add_remove_opt_retries [private] |
Number of attempts restarting the optimization to escape saddle points.
Definition at line 1397 of file multispot5.cc.
const int FitSpots::add_remove_opt_hess_inner_samples [private] |
Number of extra FFBS samples to use for computing Hessian when testing for convergence to non-saddle point.
Definition at line 1398 of file multispot5.cc.
const int FitSpots::h_outer_samples [private] |
Number of samples used for computing Hessian as part of Laplace's approximation.
Definition at line 1399 of file multispot5.cc.
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 1400 of file multispot5.cc.
const int FitSpots::tsamples [private] |
Number of samples to use in thermodynamic integration.
Definition at line 1401 of file multispot5.cc.
const Image<float> FitSpots::ave [private] |
Average of input data: used for.
Definition at line 1404 of file multispot5.cc.
ofstream& FitSpots::save_spots [private] |
Output stream for log file.
Definition at line 1408 of file multispot5.cc.
double FitSpots::time_gibbs [private] |
Benchmarking data.
Definition at line 1411 of file multispot5.cc.
double FitSpots::time_cg [private] |
Benchmarking data.
Definition at line 1412 of file multispot5.cc.
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 1420 of file multispot5.cc.
const Vector<4> FitSpots::limit [private] |
Limit vector for ConjugateGradientOnly.
Definition at line 1421 of file multispot5.cc.
const Matrix<3> FitSpots::A [private] |
Transition matrix.
Definition at line 1423 of file multispot5.cc.
const Vector<3> FitSpots::pi [private] |
Initial probabilities.
Definition at line 1424 of file multispot5.cc.
vector<vector<double> > FitSpots::pixel_intensities [private] |
Pixel intensities for all images [frame][pixel].
Definition at line 1427 of file multispot5.cc.
DataForMCMC FitSpots::data_for_t_mcmc [private] |
Aggergated data for thermodynamic integration.
Definition at line 1429 of file multispot5.cc.
DataForMCMC FitSpots::data_for_h_mcmc [private] |
Aggergated data for finding hessian.
Definition at line 1430 of file multispot5.cc.
int FitSpots::iteration [private] |
Iteration number.
Definition at line 1432 of file multispot5.cc.
1.7.4