ThreeB 1.1

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- Class List
- Class Hierarchy
- Class Members

Public Member Functions | Private Attributes
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_, UserInterfaceCallback &ui_, 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
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

## Detailed Description

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

probably could be refactored a bit.

Definition at line 1367 of file multispot5.cc.

## Constructor & Destructor Documentation

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

Definition at line 1457 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]];

        }

## 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 1535 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 1748 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 2062 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;
        }

## Member Data Documentation

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

Input data.

Definition at line 1369 of file multispot5.cc.

FitSpotsGraphics& FitSpots::graphics`[private]`

Graphics class.

Definition at line 1370 of file multispot5.cc.

UserInterfaceCallback& FitSpots::ui`[private]`

Callbacks to provide user interface.

Definition at line 1371 of file multispot5.cc.

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

Area in which to perform model fitting.

Definition at line 1372 of file multispot5.cc.

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

State in terms of current spot positions.

Definition at line 1375 of file multispot5.cc.

const int FitSpots::start_iteration`[private]`

Starting iteration number (for restarting from checkpoint)

Definition at line 1378 of file multispot5.cc.

int FitSpots::start_pass`[private]`

Starting pass (for restarting from checkpoint)

Definition at line 1379 of file multispot5.cc.

MT19937& FitSpots::rng`[private]`

Random numbewr generator.

Definition at line 1381 of file multispot5.cc.

const double FitSpots::variance`[private]`

Variance of noise in data. Must be 1.

Definition at line 1383 of file multispot5.cc.

const double FitSpots::intensity_mu`[private]`

Prior for spot intensity.

Definition at line 1384 of file multispot5.cc.

const double FitSpots::intensity_sigma`[private]`

Prior for spot intensity.

Definition at line 1385 of file multispot5.cc.

const double FitSpots::blur_mu`[private]`

Prior for spot shape.

Definition at line 1386 of file multispot5.cc.

const double FitSpots::blur_sigma`[private]`

Prior for spt shape.

Definition at line 1387 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 1396 of file multispot5.cc.

set<ImageRef> FitSpots::allowed_area`[private]`

Total allowed region, pixels dilated by area_extra_radius.

Definition at line 1397 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 1398 of file multispot5.cc.

const double FitSpots::position_prior`[private]`

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

Definition at line 1399 of file multispot5.cc.

const double FitSpots::max_motion`[private]`

Maximum motion on any axes for optimization. See ConjugateGradientOnly.

Definition at line 1402 of file multispot5.cc.

const int FitSpots::sample_iterations`[private]`

Number of mixing samples to use in Gibbs sampling.

Definition at line 1403 of file multispot5.cc.

const int FitSpots::main_cg_max_iterations`[private]`

Maximum iterations allowed for ConjugateGradientOnly in main optimization loop.

Definition at line 1408 of file multispot5.cc.

const int FitSpots::main_samples`[private]`

Number of samples to use in main loop.

Definition at line 1409 of file multispot5.cc.

const int FitSpots::main_passes`[private]`

Number of passes to perform per iteration in main loop.

Definition at line 1410 of file multispot5.cc.

const int FitSpots::outer_loop_iterations`[private]`

Total number of iterations to perform.

Definition at line 1411 of file multispot5.cc.

const int FitSpots::add_remove_tries`[private]`

Number of attemts to add/remove a spot.

Definition at line 1414 of file multispot5.cc.

const int FitSpots::add_remove_opt_samples`[private]`

Number of samples to use in model modification phase.

Definition at line 1415 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 1416 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 1417 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 1418 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 1419 of file multispot5.cc.

const int FitSpots::tsamples`[private]`

Number of samples to use in thermodynamic integration.

Definition at line 1420 of file multispot5.cc.

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

Average of input data: used for.

Definition at line 1423 of file multispot5.cc.

ofstream& FitSpots::save_spots`[private]`

Output stream for log file.

Definition at line 1427 of file multispot5.cc.

double FitSpots::time_gibbs`[private]`

Benchmarking data.

Definition at line 1430 of file multispot5.cc.

double FitSpots::time_cg`[private]`

Benchmarking data.

Definition at line 1431 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 1439 of file multispot5.cc.

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

Limit vector for ConjugateGradientOnly.

Definition at line 1440 of file multispot5.cc.

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

Transition matrix.

Definition at line 1442 of file multispot5.cc.

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

Initial probabilities.

Definition at line 1443 of file multispot5.cc.

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

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

Definition at line 1446 of file multispot5.cc.

DataForMCMC FitSpots::data_for_t_mcmc`[private]`

Aggergated data for thermodynamic integration.

Definition at line 1448 of file multispot5.cc.

DataForMCMC FitSpots::data_for_h_mcmc`[private]`

Aggergated data for finding hessian.

Definition at line 1449 of file multispot5.cc.

int FitSpots::iteration`[private]`

Iteration number.

Definition at line 1451 of file multispot5.cc.

The documentation for this class was generated from the following file:

- multispot5.cc