multispot5.cc File Reference

Fit spots to the data. More...

#include <cstdlib>
#include <cerrno>
#include <cstring>
#include <stack>
#include <algorithm>
#include <climits>
#include <iomanip>
#include <map>
#include <tr1/memory>
#include <cvd/image_io.h>
#include <cvd/image_convert.h>
#include <cvd/morphology.h>
#include <cvd/connected_components.h>
#include <cvd/draw.h>
#include <cvd/vector_image_ref.h>
#include <cvd/random.h>
#include <cvd/timer.h>
#include <gvars3/instances.h>
#include <gvars3/GUI_readline.h>
#include <gvars3/GStringUtil.h>
#include <TooN/functions/derivatives.h>
#include <TooN/determinant.h>
#include <TooN/SymEigen.h>
#include <TooN/optimization/conjugate_gradient.h>
#include "conjugate_gradient_only.h"
#include <TooN/TooN.h>
#include <utility>
#include <tr1/tuple>
#include <tr1/array>
#include <vector>
#include <cmath>
#include <cvd/image.h>
#include <string>
#include <cvd/byte.h>
#include <cassert>
#include <cvd/image_ref.h>
#include "spot_with_background.hh"
#include "randomc.h"
#include <iostream>
#include <sstream>
#include "utility.h"
#include <fstream>
#include <TooN/so2.h>
#include "mt19937.h"

Go to the source code of this file.

Classes

class  NullGraphics
 Graphics class which does absoloutely nothing. More...
class  DataForMCMC
 Closure hoding the data required do use GibbsSampler2 See FitSpots for naming of variables. More...
class  Kahan
 Class implementing the Kahan summation algorithm to allow accurate summation of very large numbers of doubles. More...
class  NegativeFreeEnergy
 Class for computing the negitve free energy using thermodynamic integration. More...
struct  IndexLexicographicPosition< Cmp, First >
 Class for sorting a list of indexes to an array of spots lexicographically according to the 2D positions of the spots. More...
struct  SampledBackgroundData
 Closure holding image data generated using samples drawn from the model. More...
struct  SpotNegProbabilityDiffWithSampledBackground
 Compute the derivative of the negative log probability with respect to the parameters of one spot, given some samples of the other spots. More...
class  FreeEnergyHessian
 Class for computing the Hessian of the negative free energy. More...
struct  LessSecond
 Comparator functor for the first element of a std::pair. More...
class  FitSpots
 Mega class which actually does the meat of the spot fitting. More...

Functions

auto_ptr< FitSpotsGraphicsnull_graphics ()
Vector spots_to_Vector (const vector< Vector< 4 > > &s)
vector< Vector< 4 > > spots_to_vector (const Vector<> &s)
Image< byte > scale_to_bytes (const Image< float > &im, float lo, float hi)
Image< byte > scale_to_bytes (const Image< float > &im)
Image< float > average_image (const vector< Image< float > > &ims)
Matrix< 4 > sampled_background_spot_hessian_ffbs (const Vector< 4 > &spot, const SampledBackgroundData &d, int bs_iterations, MT19937 &rng)
Matrix< 4 > sampled_background_spot_hessian2 (const Vector< 4 > &spot, const SampledBackgroundData &d)
Matrix< 4 > sampled_background_spot_hessian_FAKE (const Vector< 4 > &spot, const SampledBackgroundData &d)
void get_spot_pixels (const vector< ImageRef > &pixels, const Vector< 4 > &spot, vector< int > &out)
vector< string > split (const string &line)
template<class C >
string xtoa (const C &x)
template<class C >
atox (const string &s, const string &msg)
StateParameters parse_log_file (istream &in)
StateParameters generate_state_parameters_ye_olde (const BasicImage< double > &log_ratios, const vector< Image< float > > &ims, vector< ImageRef > pixels)
set< ImageRef > dilate_mask (const vector< ImageRef > &v, double r)
double brightness_motion_limit (double mu, double sigma, bool not_one)
void fit_spots_historic (const BasicImage< double > &log_ratios, const vector< Image< float > > &ims, vector< ImageRef > pixels, FitSpotsGraphics &graphics)
void fit_spots_new (const vector< Image< float > > &ims, StateParameters &p, ofstream &save_spots, FitSpotsGraphics &gr)

Detailed Description

Fit spots to the data.

Definition in file multispot5.cc.


Function Documentation

Matrix<4> sampled_background_spot_hessian_ffbs ( const Vector< 4 > &  spot,
const SampledBackgroundData d,
int  bs_iterations,
MT19937 rng 
)

Compute the Hessian of the log probability.

The background is sampled rather sparsely, and the spot in question is sampled much more densely using FFBS.

Parameters:
spot Spot parameters
d Background brightness (from other spots)
bs_iterations Exter backward sampling iterations
rng Random number generator
Returns:
the Hessian of the log probability around the spot

Definition at line 740 of file multispot5.cc.

References SampledBackgroundData::A, SampledMultispot::compute_spot_intensity(), SampledMultispot::compute_spot_intensity_hessian(), diff_log_log_normal(), forward_algorithm_delta(), hess_log_log_normal(), SampledBackgroundData::mu_blur, SampledBackgroundData::mu_brightness, SampledBackgroundData::O, SampledBackgroundData::pi, SampledBackgroundData::pixel_intensities, SampledBackgroundData::pixels, SampledBackgroundData::sample_intensities_without_spot, SampledBackgroundData::sigma_blur, SampledBackgroundData::sigma_brightness, and SampledBackgroundData::variance.

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

00741 {
00742     vector<tuple<double, Vector<4>, Matrix<4> > > spot_hess_etc = compute_spot_intensity_hessian(d.pixels, spot);
00743     vector<double> spot_intensities = compute_spot_intensity(d.pixels, spot);
00744 
00745     Matrix<4> sum_hess_log = Zeros;
00746     Matrix<4> sum_diff2_log = Zeros;
00747 
00748     vector<State> current_sample;
00749 
00750     const unsigned int nframes = d.pixel_intensities.size();
00751     const unsigned int npixels = d.pixels.size();
00752 
00753     Matrix<4> sum_hess  = Zeros;
00754     Vector<4> sum_deriv = Zeros;
00755     
00756     vector<pair<Matrix<4>, Vector<4> > > hess_and_deriv_part(nframes);
00757     
00758     for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++)
00759     {
00760         SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensities, d.pixel_intensities, d.variance);
00761         
00762         //Compute what the per-frame hess and deriv parts are
00763         //if the spot is on in a frame.
00764         for(unsigned int frame=0; frame < nframes; frame++)
00765         {
00766             Matrix<4> hess = Zeros;
00767             Vector<4> deriv = Zeros;
00768 
00769             for(unsigned int pixel=0; pixel < npixels; pixel++)
00770             {
00771                 double e = d.pixel_intensities[frame][pixel] - (d.sample_intensities_without_spot[s][frame][pixel] + spot_intensities[pixel]);
00772                 //Build up the derivative
00773                 hess += e * get<2>(spot_hess_etc[pixel]) - get<1>(spot_hess_etc[pixel]).as_col() * get<1>(spot_hess_etc[pixel]).as_row();
00774                 deriv += e * get<1>(spot_hess_etc[pixel]);
00775             }
00776             hess_and_deriv_part[frame] = make_pair(hess, deriv);
00777         }
00778         
00779         //Forward filtering
00780         std::vector<array<double, 3> > delta = forward_algorithm_delta(d.A, d.pi, B, d.O);
00781         
00782         for(int i=0; i < bs_iterations; i++)
00783         {
00784             current_sample = backward_sampling<3,State>(d.A, delta, rng);
00785 
00786             Matrix<4> hess = Zeros;
00787             Vector<4> deriv = Zeros;
00788             for(unsigned int frame=0; frame < nframes; frame++)
00789                 if(current_sample[frame] == 0)
00790                 {
00791                         hess += hess_and_deriv_part[frame].first;
00792                         deriv += hess_and_deriv_part[frame].second;
00793                 }
00794 
00795             sum_hess += hess + deriv.as_col() * deriv.as_row();
00796             sum_deriv += deriv;
00797         }
00798     }
00799 
00800     sum_hess  /= (bs_iterations * d.sample_intensities_without_spot.size() * d.variance);
00801     sum_deriv /= (bs_iterations * d.sample_intensities_without_spot.size() * d.variance);
00802 
00803     sum_hess -= sum_deriv.as_col() * sum_deriv.as_row();
00804 
00805     sum_hess[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness);
00806     sum_hess[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);
00807     sum_deriv[0] += diff_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness);
00808     sum_deriv[1] += diff_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);
00809 
00810     //cout << "Turboderiv:" << sum_deriv << endl;
00811     //cout << "Turbohess:\n" << sum_hess << endl;
00812 
00813     return sum_hess;
00814 }

Matrix<4> sampled_background_spot_hessian2 ( const Vector< 4 > &  spot,
const SampledBackgroundData d 
)

Debugging function. Not mathematically correct. Do not use.

Definition at line 817 of file multispot5.cc.

References SampledBackgroundData::A, SampledMultispot::compute_spot_intensity_hessian(), diff_log_log_normal(), forward_algorithm_hessian(), hess_log_log_normal(), SampledBackgroundData::mu_blur, SampledBackgroundData::mu_brightness, SampledBackgroundData::O, SampledBackgroundData::pi, SampledBackgroundData::pixel_intensities, SampledBackgroundData::pixels, SampledBackgroundData::sample_intensities_without_spot, SampledBackgroundData::sigma_blur, SampledBackgroundData::sigma_brightness, and SampledBackgroundData::variance.

00818 {
00819     vector<tuple<double, Vector<4>, Matrix<4> > > spot_intensities = compute_spot_intensity_hessian(d.pixels, spot);
00820 
00821     Matrix<4> sum_hess_log = Zeros;
00822     Matrix<4> sum_diff2_log = Zeros;
00823 
00824     for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++)
00825     {
00826         SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensities, d.pixel_intensities, d.variance);
00827 
00828         double prob;
00829         Vector<4> diff;
00830         Matrix<4> hess;
00831 
00832         tie(prob, diff, hess) = forward_algorithm_hessian(d.A, d.pi, B, d.O);
00833         
00834         sum_hess_log += hess;
00835 
00836         diff += makeVector(diff_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness), diff_log_log_normal(spot[1], d.mu_blur, d.sigma_blur), 0, 0);
00837         sum_diff2_log += diff.as_col() * diff.as_row();
00838     }
00839 
00840     Matrix<4> hess_log = sum_hess_log / d.sample_intensities_without_spot.size();
00841     Matrix<4> diff2_log = sum_diff2_log / d.sample_intensities_without_spot.size();
00842 
00843     //Add in the prior
00844 
00845     hess_log[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness);
00846     hess_log[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);
00847 
00848     return hess_log + diff2_log;
00849 }

Matrix<4> sampled_background_spot_hessian_FAKE ( const Vector< 4 > &  spot,
const SampledBackgroundData d 
)

Debugging function. Not mathematically correct. Do not use.

Definition at line 852 of file multispot5.cc.

References SampledBackgroundData::A, SampledMultispot::compute_spot_intensity_hessian(), forward_algorithm_hessian(), hess_log_log_normal(), SampledBackgroundData::mu_blur, SampledBackgroundData::mu_brightness, SampledBackgroundData::O, SampledBackgroundData::pi, SampledBackgroundData::pixel_intensities, SampledBackgroundData::pixels, SampledBackgroundData::sample_intensities_without_spot, SampledBackgroundData::sigma_blur, SampledBackgroundData::sigma_brightness, and SampledBackgroundData::variance.

00853 {
00854     vector<tuple<double, Vector<4>, Matrix<4> > > spot_intensities = compute_spot_intensity_hessian(d.pixels, spot);
00855 
00856     Matrix<4> sum_hess_log = Zeros;
00857 
00858     for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++)
00859     {
00860         SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensities, d.pixel_intensities, d.variance);
00861 
00862         double prob;
00863         Vector<4> diff;
00864         Matrix<4> hess;
00865 
00866         tie(prob, diff, hess) = forward_algorithm_hessian(d.A, d.pi, B, d.O);
00867         
00868         sum_hess_log += hess;
00869     }
00870 
00871     Matrix<4> hess_log = sum_hess_log / d.sample_intensities_without_spot.size();
00872 
00873     //Add in the prior
00874 
00875     hess_log[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness);
00876     hess_log[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);
00877     
00878     return hess_log;
00879 }

StateParameters generate_state_parameters_ye_olde ( const BasicImage< double > &  log_ratios,
const vector< Image< float > > &  ims,
vector< ImageRef >  pixels 
)

Setup the parameters for a run using the old and deeply peculiar method.

This includes the unpleasant and diffucult so use de-checkpointing code. wtf. The use of this function is very strongly deprecated.

Parameters:
log_ratios Image from which region is selected.
ims Input data
pixels Region for spot fitting to run in

Definition at line 1148 of file multispot5.cc.

References assert_same_size(), StateParameters::iteration, log_normal_mode(), StateParameters::pass, StateParameters::pixels, StateParameters::rng, StateParameters::spots, and spots_to_vector().

01148                                                                                                                                                   {
01149     sort(pixels.begin(), pixels.end());
01150 
01151     const double variance = 1; // it should be
01152 
01153     //To scale the X axis of a log-normal distribution, only
01154     //the mu parameter needs to be changed...
01155     const double intensity_mu = GV3::get<double>("intensity.rel_mu", 0., -1)  + log(sqrt(variance));
01156     const double intensity_sigma = GV3::get<double>("intensity.rel_sigma", 0., -1);
01157     const double blur_mu = GV3::get<double>("blur.mu", 0., -1);
01158     const double blur_sigma = GV3::get<double>("blur.sigma", 0., -1);
01159 
01160     //The region was extracted at a certain threshold.
01161     //These regions may be too small, so some post-region finding dilation 
01162     //may be performed. New spots are only placed at pixels which exceed the threshold.
01163     //post_dilate.threshold is (if set) used as the placing threshold so the placing threshold
01164     //can be different from the region-finding threshold.
01165 
01166     //Note that as a result of dliation, regions of <pixels> may be below the threshold.
01167     //In the historic version, this could affect new spot placement. This feature is not supported
01168     //in this version.
01169     double threshold = GV3::get<double>("threshold", 0, -1);
01170     const double post_threshold = GV3::get<double>("post_dilate.threshold", -1, 1);
01171     if(post_threshold != -1)
01172         threshold = post_threshold;
01173 
01174     
01175     //If dilation after region finding is to be performed, then do it here.
01176     const double post_dilate_radius = GV3::get<double>("post_dilate.radius", 0, -1);
01177     if(post_dilate_radius != 0)
01178     {
01179         Image<byte> pix(ims[0].size());
01180         pix.fill(0);
01181         
01182         for(unsigned int i=0; i < pixels.size(); i++)
01183             pix[pixels[i]] = 255;
01184 
01185         Image<byte> dilated = morphology(pix, getDisc(post_dilate_radius), Morphology::BinaryDilate<byte>());
01186 
01187         pixels.clear();
01188 
01189         ImageRef p(0,0);
01190         do
01191             if(dilated[p])
01192                 pixels.push_back(p);
01193         while(p.next(dilated.size()));
01194     }
01195 
01196 
01197     assert_same_size(ims);
01198     if(log_ratios.size() != ims[0].size())
01199     {
01200         cerr << "Bad log ratios size\n";
01201         exit(1);
01202     }
01203 
01204     vector<Vector<4> > spots;
01205     //Spots can be either put down automatically, or specified 
01206     //The auto-initialization is very strange.
01207     if(GV3::get<bool>("spots.auto_initialise", 1, 1))
01208     {
01209         //You never get two spots in the same disc in the second stage of the algorithm
01210         vector<ImageRef> disc = getDisc(GV3::get<double>("spot_spread", 3.1, 1));
01211         
01212         //Record all the pixels
01213         map<ImageRef, double> valid_pixels;
01214         for(unsigned int i=0; i < pixels.size(); i++)
01215             if(log_ratios[pixels[i]] > threshold)
01216                 valid_pixels.insert(make_pair(pixels[i], log_ratios[pixels[i]]));
01217 
01218 
01219         //Get some initial spots by finding the local maxima
01220         ImageRef neighbours[8] = {
01221             ImageRef(-1, -1),
01222             ImageRef( 0, -1),
01223             ImageRef( 1, -1),
01224 
01225             ImageRef(-1,  0),
01226             ImageRef( 1,  0),
01227 
01228             ImageRef(-1,  1),
01229             ImageRef( 0,  1),
01230             ImageRef( 1,  1),
01231         };
01232         for(unsigned int i=0; i < pixels.size(); i++)
01233         {
01234             if(!(log_ratios[pixels[i]] > threshold))
01235                 goto not_a_maximum;
01236 
01237             for(int j=0; j < 8; j++)
01238                 if(!log_ratios.in_image(pixels[i] + neighbours[j]) || ! (log_ratios[pixels[i]] > log_ratios[pixels[i] + neighbours[j]]))
01239                     goto not_a_maximum;
01240             
01241             spots.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), log_normal_mode(blur_mu, blur_sigma), pixels[i].x, pixels[i].y));
01242 
01243             //Remove the pixels around the initial spots
01244             for(unsigned int j=0; j < disc.size(); j++)
01245                 valid_pixels.erase(pixels[i] + disc[j]);
01246 
01247             not_a_maximum:;
01248         }
01249 
01250         for(unsigned int i=0; i < spots.size(); i++)
01251             cout << spots[i] << endl;
01252         
01253         //Now place down extra spots in the remaining space.
01254         while(!valid_pixels.empty())
01255         {
01256             ImageRef p = max_element(valid_pixels.begin(), valid_pixels.end(), LessSecond())->first;
01257             spots.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), log_normal_mode(blur_mu, blur_sigma), p.x, p.y));
01258 
01259             for(unsigned int j=0; j < disc.size(); j++)
01260                 valid_pixels.erase(p + disc[j]);
01261         }
01262         
01263         //This line allows extra spots to be placed down around each spot already put down.
01264         //This is a shocking hack and jenerally very unpleasant.
01265         double extra_r = GV3::get<double>("extra_spots", 0, 1);
01266         vector<ImageRef> extra = getDisc(extra_r);
01267         vector<Vector<4> > more_spots;
01268         for(unsigned int i=0; i < extra.size(); i++)
01269             if(extra[i] != ImageRef_zero)
01270                 for(unsigned int j=0; j < spots.size(); j++)
01271                     more_spots.push_back(spots[j] + makeVector(0, 0, extra[i].x, extra[i].y) / (2*extra_r+1));
01272 
01273         copy(more_spots.begin(), more_spots.end(), back_inserter(spots));
01274     }
01275     else
01276     {
01277         Vector<>  loaded_spots = GV3::get<Vector<> >("spots.manual_spots", "", -1);
01278 
01279         if(loaded_spots.size()%4 != 0)
01280         {
01281             cerr << "Loaded spot size is not a multiple of 4\n";
01282             exit(1);
01283         }
01284 
01285         else    
01286             spots = spots_to_vector(loaded_spots);
01287     }
01288 
01289     //Initialize the MT19937 RNG from a seed.
01290     shared_ptr<MT19937> rng(new MT19937);
01291     rng->simple_seed(GV3::get<int>("seed", 0, 1));
01292         
01293     //Load in a checkpoint (precise RNG state, iteration and pass).
01294     int start_iteration=0;
01295     int start_pass=0;
01296     if(GV3::get<bool>("checkpoint", 0, 1))
01297     {
01298         string rng_state = GV3::get<string>("checkpoint.rng.state", "", -1);
01299         istringstream rs(rng_state);
01300         rng->read(rs);
01301         start_iteration=GV3::get<int>("checkpoint.iteration", 0, -1);
01302         start_pass=GV3::get<int>("checkpoint.pass", 0, -1);
01303     }
01304     
01305     StateParameters p;
01306     p.spots = spots;
01307     p.rng = rng;
01308     p.pass = start_pass;
01309     p.iteration = start_iteration;
01310     p.pixels = pixels;
01311 
01312     return p;
01313 }

void fit_spots_historic ( const BasicImage< double > &  log_ratios,
const vector< Image< float > > &  ims,
vector< ImageRef >  pixels,
FitSpotsGraphics graphics 
)

The new function has started to diverge from this one.

1. There is no placement code to perform post-dilation and re-thresholding. 2. The new function has support for position priors. by setting use_position_prior=1 3. limit is Ones * max_motion. The new function scales the brightness prior.

Parameters:
log_ratios Image from which region is selected.
ims Input data
pixels Region for spot fitting to run in
graphics Graphics class for initialisation

Definition at line 2068 of file multispot5.cc.

References SampledMultispot::add_spot(), assert_same_size(), average_image(), boundingbox(), SampledMultispot::compute_spot_intensity(), FitSpotsGraphics::draw_krap(), FitSpotsGraphics::draw_pixels(), get_spot_pixels(), FreeEnergyHessian::hessian(), ConjugateGradientOnly< Size, Precision >::init(), FitSpotsGraphics::init(), ConjugateGradientOnly< Size, Precision >::iterate(), log_normal_mode(), ConjugateGradientOnly< Size, Precision >::max_iterations, MT19937::read(), SampledMultispot::remove_spot(), sampled_background_spot_hessian_ffbs(), scale_to_bytes(), SampledMultispot::sequence(), MT19937::simple_seed(), spots_to_Vector(), spots_to_vector(), FitSpotsGraphics::swap(), MT19937::write(), and ConjugateGradientOnly< Size, Precision >::x.

02069 {
02070     sort(pixels.begin(), pixels.end());
02071 
02072 
02073     const double post_dilate_radius = GV3::get<double>("post_dilate.radius", 0, -1);
02074 
02075     if(post_dilate_radius != 0)
02076     {
02077         Image<byte> pix(ims[0].size());
02078         pix.fill(0);
02079         
02080         for(unsigned int i=0; i < pixels.size(); i++)
02081             pix[pixels[i]] = 255;
02082 
02083         Image<byte> dilated = morphology(pix, getDisc(post_dilate_radius), Morphology::BinaryDilate<byte>());
02084 
02085         pixels.clear();
02086 
02087         ImageRef p(0,0);
02088         do
02089             if(dilated[p])
02090                 pixels.push_back(p);
02091         while(p.next(dilated.size()));
02092     }
02093 
02094 
02095     const double variance = 1; // it should be
02096 
02097     //To scale the X axis of a log-normal distribution, only
02098     //the mu parameter needs to be changed...
02099     const double intensity_mu = GV3::get<double>("intensity.rel_mu", 0., -1)  + log(sqrt(variance));
02100     const double intensity_sigma = GV3::get<double>("intensity.rel_sigma", 0., -1);
02101     const double blur_mu = GV3::get<double>("blur.mu", 0., -1);
02102     const double blur_sigma = GV3::get<double>("blur.sigma", 0., -1);
02103 
02104     double threshold = GV3::get<double>("threshold", 0, -1);
02105 
02106 
02107     const double post_threshold = GV3::get<double>("post_dilate.threshold", -1, 1);
02108     if(post_threshold != -1)
02109         threshold = post_threshold;
02110 
02111     //General optimizing and sampling parameters
02112     const double max_motion = GV3::get<double>("cg.max_motion", 0., -1);
02113     const int sample_iterations = GV3::get<int>("gibbs.mixing_iterations", 0, -1);
02114     
02115     //Task specific optimizing and sampling parameters:
02116 
02117     //Main optimization loop
02118     const int main_cg_max_iterations = GV3::get<double>("main.cg.max_iterations", 0., -1);
02119     const int main_samples = GV3::get<int>("main.gibbs.samples", 0, -1);
02120     const int main_passes = GV3::get<int>("main.passes", 0, -1);
02121     const int outer_loop_iterations = GV3::get<int>("main.total_iterations", 100000000, 1);
02122     
02123     //Spot selection loop
02124     const int add_remove_tries = GV3::get<int>("add_remove.tries", 0, -1);
02125     const int add_remove_opt_samples = GV3::get<int>("add_remove.optimizer.samples", 0, -1);
02126     const int add_remove_opt_retries = GV3::get<int>("add_remove.optimizer.attempts", 0, -1);
02127     const int add_remove_opt_hess_inner_samples = GV3::get<int>("add_remove.optimizer.hessian_inner_samples", 0, -1);
02128     const int h_outer_samples = GV3::get<int>("add_remove.hessian.outer_samples", 0, -1);
02129     const int h_inner_samples = GV3::get<int>("add_remove.hessian.inner_samples", 0, -1);
02130     const int tsamples = GV3::get<int>("add_remove.thermo.samples", 0, -1);
02131 
02132     ofstream save_spots;
02133     save_spots.open(GV3::get<string>("save_spots", "", -1).c_str());
02134     int err = errno;
02135 
02136     if(!save_spots.good())
02137     {
02138         cerr << "***********************************************************\n";
02139         cerr << "WARNING: failed to open " << GV3::get<string>("save_spots") << ": " <<strerror(err) << endl;
02140         cerr << "***********************************************************\n";
02141         return;
02142     }
02143 
02144     const Matrix<3> A  = GV3::get<Matrix<3> >("A", Zeros, 1);
02145     const Vector<3> pi = GV3::get<Vector<3> >("pi", Zeros, 1);
02146 
02147 /*  
02148     SGRAPHICS(
02149         int debug_window_zoom = GV3::get<int>("debug.zoom", 3, 1);
02150         GLWindowWrapper win(log_ratios.size()*debug_window_zoom);
02151     );
02152 
02153     GRAPHICS(
02154         set_GL_zoom_size(log_ratios.size(), debug_window_zoom);
02155 
02156         glEnable(GL_BLEND);
02157         glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
02158         glEnable(GL_LINE_SMOOTH);
02159     );
02160 */
02161     graphics.init(log_ratios.size());
02162 
02163     assert_same_size(ims);
02164     if(log_ratios.size() != ims[0].size())
02165     {
02166         cerr << "Bad log ratios size\n";
02167         exit(1);
02168     }
02169 
02170     const Image<float> ave = average_image(ims);
02171     
02172     vector<Vector<4> > spots;
02173     
02174     if(GV3::get<bool>("spots.auto_initialise", 1, 1))
02175     {
02176         vector<ImageRef> disc = getDisc(GV3::get<double>("spot_spread", 3.1, 1));
02177         
02178         //Record all the pixels
02179         map<ImageRef, double> valid_pixels;
02180         for(unsigned int i=0; i < pixels.size(); i++)
02181             if(log_ratios[pixels[i]] > threshold)
02182                 valid_pixels.insert(make_pair(pixels[i], log_ratios[pixels[i]]));
02183 
02184 
02185         //Get some initial spots by finding the local maxima
02186         ImageRef neighbours[8] = {
02187             ImageRef(-1, -1),
02188             ImageRef( 0, -1),
02189             ImageRef( 1, -1),
02190 
02191             ImageRef(-1,  0),
02192             ImageRef( 1,  0),
02193 
02194             ImageRef(-1,  1),
02195             ImageRef( 0,  1),
02196             ImageRef( 1,  1),
02197         };
02198         for(unsigned int i=0; i < pixels.size(); i++)
02199         {
02200             if(!(log_ratios[pixels[i]] > threshold))
02201                 goto not_a_maximum;
02202 
02203             for(int j=0; j < 8; j++)
02204                 if(!log_ratios.in_image(pixels[i] + neighbours[j]) || ! (log_ratios[pixels[i]] > log_ratios[pixels[i] + neighbours[j]]))
02205                     goto not_a_maximum;
02206             
02207             spots.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), log_normal_mode(blur_mu, blur_sigma), pixels[i].x, pixels[i].y));
02208 
02209             //Remove the pixels around the initial spots
02210             for(unsigned int j=0; j < disc.size(); j++)
02211                 valid_pixels.erase(pixels[i] + disc[j]);
02212 
02213             not_a_maximum:;
02214         }
02215 
02216         //Are they OK?
02217         /*GRAPHICS(
02218             glDrawPixels(scale_to_bytes(ave));
02219             glColor3f(1, 0, 0);
02220             draw_bbox(boundingbox(pixels));
02221             
02222             draw_pixels(pixels, 0, 0, 1, .5);
02223             glBegin(GL_LINES);
02224             glColor3f(1, 0, 0);
02225         );*/
02226         for(unsigned int i=0; i < spots.size(); i++)
02227         {
02228             cout << spots[i] << endl;
02229             //GRAPHICS(glDrawCross(spots[i].slice<2, 2>()););
02230         }
02231         /*GRAPHICS(
02232             glEnd();
02233             
02234             glFlush();
02235             win.swap_buffers();
02236 
02237             glPrintErrors();
02238             cout << "Hello\n";
02239             GUI_Pause();
02240         );*/
02241         
02242         while(!valid_pixels.empty())
02243         {
02244             ImageRef p = max_element(valid_pixels.begin(), valid_pixels.end(), LessSecond())->first;
02245             spots.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), log_normal_mode(blur_mu, blur_sigma), p.x, p.y));
02246 
02247             for(unsigned int j=0; j < disc.size(); j++)
02248                 valid_pixels.erase(p + disc[j]);
02249             
02250             /*GRAPHICS(
02251                 glDrawPixels(scale_to_bytes(ave));
02252                 glColor3f(1, 0, 0);
02253                 draw_bbox(boundingbox(pixels));
02254                 
02255                 glBegin(GL_LINES);
02256                 for(unsigned int i=0; i < spots.size(); i++)
02257                     glDrawCross(spots[i].slice<2, 2>());
02258                 glEnd();
02259                 glFlush();
02260                 
02261                 win.swap_buffers();
02262 
02263                 glPrintErrors();
02264                 GUI_Pause();
02265             );*/
02266         }
02267 
02268         //spots.resize(2);
02269         //spots[0] = makeVector(4.050783236168317e+00, 6.893690887878990e-01, 4.165886930921706e+01, 1.184866845082972e+02);
02270         //spots[1] = makeVector(4.050783236168317e+00, 6.893690887878990e-01, 4.165886930921706e+01, 1.184866845082972e+02+2.5);
02271 
02272         
02273         double extra_r = GV3::get<double>("extra_spots", 0, 1);
02274         vector<ImageRef> extra = getDisc(extra_r);
02275         vector<Vector<4> > more_spots;
02276         for(unsigned int i=0; i < extra.size(); i++)
02277             if(extra[i] != ImageRef_zero)
02278                 for(unsigned int j=0; j < spots.size(); j++)
02279                     more_spots.push_back(spots[j] + makeVector(0, 0, extra[i].x, extra[i].y) / (2*extra_r+1));
02280 
02281         copy(more_spots.begin(), more_spots.end(), back_inserter(spots));
02282             
02283     }
02284     else
02285     {
02286         Vector<>  loaded_spots = GV3::get<Vector<> >("spots.manual_spots", "", -1);
02287 
02288         if(loaded_spots.size()%4 != 0)
02289         {
02290             cerr << "Loaded spot size is not a multiple of 4\n";
02291             exit(1);
02292         }
02293 
02294         else    
02295             spots = spots_to_vector(loaded_spots);
02296     
02297             /*GRAPHICS(
02298                 glDrawPixels(scale_to_bytes(ave));
02299                 glColor3f(1, 0, 0);
02300                 draw_bbox(boundingbox(pixels));
02301                 
02302                 glBegin(GL_LINES);
02303                 for(unsigned int i=0; i < spots.size(); i++)
02304                     glDrawCross(spots[i].slice<2, 2>());
02305                 glEnd();
02306                 glFlush();
02307                 
02308                 win.swap_buffers();
02309 
02310                 glPrintErrors();
02311                 GUI_Pause();
02312             );*/
02313 
02314     }
02315 
02316 
02317     
02318     //Pixel intensities for all images [frame][pixel]
02319     vector<vector<double> > pixel_intensities(ims.size(), vector<double>(pixels.size()));
02320     for(unsigned int frame=0; frame < ims.size(); frame++)
02321         for(unsigned int p=0; p < pixels.size(); p++)
02322             pixel_intensities[frame][p] = ims[frame][pixels[p]];
02323 
02324 
02325 
02326     MT19937 rng;
02327     rng.simple_seed(GV3::get<int>("seed", 0, 1));
02328         
02329     int start_iteration=0;
02330     int start_pass=0;
02331     if(GV3::get<bool>("checkpoint", 0, 1))
02332     {
02333         string rng_state = GV3::get<string>("checkpoint.rng.state", "", -1);
02334         istringstream rs(rng_state);
02335         rng.read(rs);
02336         start_iteration=GV3::get<int>("checkpoint.iteration", 0, -1);
02337         start_pass=GV3::get<int>("checkpoint.pass", 0, -1);
02338     }
02339 
02340     DataForMCMC data_for_t_mcmc(pixels, pixel_intensities, intensity_mu, intensity_sigma, blur_mu, blur_sigma, variance, tsamples, sample_iterations, A, pi, rng);
02341     DataForMCMC data_for_h_mcmc(pixels, pixel_intensities, intensity_mu, intensity_sigma, blur_mu, blur_sigma, variance, h_outer_samples, sample_iterations, A, pi, rng);
02342 
02343     //Outer loop like Multispot2:
02344     //
02345     
02346     double time_gibbs=0;
02347     double time_cg = 0;
02348 
02349     const Vector<4> limit = Ones * max_motion;
02350 
02351     for(int iteration=start_iteration; iteration < outer_loop_iterations ;iteration++)
02352     {
02353         save_spots << "Iteration: " << iteration << endl;
02354         save_spots << "MAIN: " << setprecision(20) << scientific << spots_to_Vector(spots) << endl;
02355 
02356         cout << endl << endl << "----------------------" << endl << "Optimizing:\n";
02357         cout << spots.size() << endl;
02358 
02359         {
02360             //Precompute the intensities for all spot pixels
02361             vector<vector<double> > spot_intensities; //[spot][pixel]
02362             for(unsigned int i=0; i < spots.size(); i++)
02363                 spot_intensities.push_back(compute_spot_intensity(pixels, spots[i]));
02364 
02365             //Which pixels does each spot have?
02366             vector<vector<int> > spot_pixels; //[spot][pixel]
02367             spot_pixels.resize(spots.size());
02368             for(unsigned int s=0; s < spots.size(); s++)
02369                 get_spot_pixels(pixels, spots[s], spot_pixels[s]);
02370 
02371             
02372             //Optimize the model, N spots at a time.
02373             //
02374             for(int pass=start_pass; pass < main_passes; pass++)
02375             {
02376                 save_spots << "Pass: " << pass << endl;
02377                 rng.write(save_spots);
02378                 save_spots << endl;
02379 
02380                 start_pass=0; // This is only nonzero first time, since we might chekpoint midway through an iteration
02381                 save_spots << "PASS" << pass << ": " << setprecision(20) << scientific << spots_to_Vector(spots) << endl;
02382                 save_spots << "ENDCHECKPOINT" << endl << flush;
02383                 //Sort spots according to pass%4
02384 
02385                 //Sort the spots so that the optimization runs in sweeps
02386                 //This heiristic seems to increase the rate of propagation of information
02387                 //about spot positions.
02388 
02389                 //Create a list of indices
02390                 vector<int> index = sequence(spots.size());
02391                 
02392                 int passs = pass + iteration;
02393                 //Sort the indices according to the position of the spot that they index
02394                 if(passs%4 == 0)
02395                     sort(index.begin(), index.end(), IndexLexicographicPosition<less<double>, 2>(spots));
02396                 else if(passs%4==1)
02397                     sort(index.begin(), index.end(), IndexLexicographicPosition<greater<double>, 2>(spots));
02398                 else if(passs%4==1)
02399                     sort(index.begin(), index.end(), IndexLexicographicPosition<less<double>, 3>(spots));
02400                 else
02401                     sort(index.begin(), index.end(), IndexLexicographicPosition<greater<double>, 3>(spots));
02402 
02403                 //Reorder the spots and their intensities and their pixel lists
02404                 {
02405                     vector<Vector<4> > tmp_spot(index.size());
02406                     vector<vector<double> > tmp_spot_intensities(index.size());
02407                     vector<vector<int> > tmp_spot_pixels(index.size());
02408                     for(unsigned int i=0; i < index.size(); i++)
02409                     {
02410                         tmp_spot[i] = spots[index[i]];
02411                         swap(tmp_spot_intensities[i], spot_intensities[index[i]]);
02412                         swap(tmp_spot_pixels[i], spot_pixels[i]);
02413                     }
02414 
02415                     swap(tmp_spot, spots);
02416                     swap(tmp_spot_intensities, spot_intensities);
02417                     swap(tmp_spot_pixels, spot_pixels);
02418                 }
02419 
02420                 //Sweep through and optimize each spot in turn
02421                 for(int s=0; s < (int)spots.size(); s++)
02422                 {
02423                     timer.reset();
02424                     //Get some samples with Gibbs sampling
02425                     vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
02426                     vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
02427 
02428                     GibbsSampler2 sampler(pixel_intensities, spot_intensities, spots, spot_pixels, A, pi, variance, sample_iterations);
02429                     for(int i=0; i < main_samples; i++)
02430                     {
02431                         sampler.next(rng);
02432                         sample_list.push_back(sampler.sample());
02433                         sample_intensities.push_back(sampler.sample_intensities());
02434                     }
02435 
02436                     //First, remove the spot from all the samples.
02437                     for(unsigned int i=0; i < sample_list.size(); i++)
02438                         remove_spot(sample_intensities[i], spot_intensities[s], sample_list[i][s]);
02439                     
02440                     //cout << timer.get_time() << endl;
02441                     time_gibbs += timer.reset();
02442 
02443                     //Package up all the data
02444                     SampledBackgroundData data(sample_intensities, pixel_intensities, pixels, 
02445                                                intensity_mu, intensity_sigma, blur_mu, blur_sigma, 
02446                                                A, pi, variance);
02447                     
02448                     //Derivative computer:
02449                     SpotNegProbabilityDiffWithSampledBackground compute_deriv(data);
02450                     
02451 
02452                     graphics.draw_pixels(pixels, 0, 0, 1, 1);
02453                     graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), s);
02454                     graphics.swap();
02455 
02456                     //Optimize spot "s"
02457                     //Optimize with the derivatives only since the actual probability
02458                     //is much harder to compute
02459                     ConjugateGradientOnly<4> cg(spots[s], compute_deriv, limit);
02460 
02461 
02462                     cg.max_iterations = main_cg_max_iterations;
02463 
02464 
02465                     #if 0
02466                     cout << setprecision(10);
02467                     cout << spots_to_Vector(spots) << endl;
02468                     Matrix<4> hess, hess_errors;
02469                     cout << "Hello, my name is Inigo Montoya\n";
02470                     /*for(int i=0; i < 4; i++)
02471                     {
02472                         Matrix<4, 2> d = numerical_gradient_with_errors(NthDeriv(compute_deriv, i), cg.x);
02473                         hess[i] = d.T()[0];
02474                         hess_errors[i] = d.T()[1];
02475                     }
02476                     */
02477                     //cout << "Errors:\n" << hess_errors << endl;
02478                     //cout << "NHess:\n" << hess<< endl;
02479                     Matrix<4> rhess =  -sampled_background_spot_hessian(cg.x, data);
02480                     cout << "Hess:\n" << rhess << endl;
02481                     //cout << "Err:\n" << hess - rhess << endl;
02482 
02483                     //Vector<4> grad = compute_deriv(cg.x);
02484 
02485                     //Matrix<4> e = hess - rhess;
02486 
02487                     //for(int i=0; i < 4; i++)
02488                     //  for(int j=0; j < 4; j++)
02489                     //      e[i][j] /= hess_errors[i][j];
02490 
02491                     //cout << "Err:\n" << e << endl;
02492                     cout << "Deriv:" <<  compute_deriv(cg.x) << endl;
02493                     //cout << "Full:\n" << sampled_background_spot_hessian2(cg.x, data) - grad.as_col()*grad.as_row() << endl;
02494 
02495                     FreeEnergyHessian hesscomputer(data_for_h_mcmc);
02496 
02497                     Matrix<4> nhess = hesscomputer.hessian(spots, 0);
02498                     cout << "NHess:\n" << nhess << endl;
02499 
02500                     cout << "Turbo-N Hess:\n" << sampled_background_spot_hessian_ffbs(cg.x, data, 10000) << endl;
02501 
02502                     cout << "TI energy: " << NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(spots)) << endl;
02503                     cout << "FA energy: " <<  SpotProbabilityWithSampledBackgroundFAKE(data)(cg.x) << endl;
02504 
02505 
02506 
02507                     //cout << "Numerical hessian from FD:\n" << numerical_hessian(SpotProbabilityWithSampledBackgroundFAKE(data), cg.x) << endl;
02508                     exit(0);
02509                     #endif
02510                     //cout << "Starting opt... " << cg.x << endl;
02511                     while(cg.iterate(compute_deriv))
02512                     {
02513                         graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), s, cg.x);
02514                         graphics.draw_pixels(pixels, 0, 0, 1, .2);
02515                         graphics.swap();
02516                         //cout << cg.x << endl;
02517                     }
02518 
02519                     //Update to use the result of the optimization
02520                     spots[s] = cg.x;
02521                     //cout << "End: " << cg.x << endl;
02522                     
02523                     graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), -1);
02524                     graphics.swap();
02525                         
02526                     //Recompute the new spot intensity, since the spot has changed
02527                     spot_intensities[s] = compute_spot_intensity(pixels, spots[s]);
02528 
02529                     //Recompute which are the useful pixels
02530                     get_spot_pixels(pixels, spots[s], spot_pixels[s]);
02531 
02532                     if(spot_pixels[s].empty())
02533                     {
02534                         //Spot has been ejected. Erase it.
02535                         cout  << " Erasing ejected spot\n";
02536                         cout << spots[s] << endl;
02537                         //GUI_Pause(1);
02538 
02539                         spot_intensities.erase(spot_intensities.begin() + s);
02540                         spot_pixels.erase(spot_pixels.begin() + s);
02541                         spots.erase(spots.begin() + s);
02542                         s--;
02543                         //exit(0);
02544                     }
02545                     
02546                     //cout << timer.get_time() << endl;
02547                     time_cg += timer.reset();
02548 
02549                     //cout << "Times: " << time_gibbs << " " << time_cg << endl;
02550                     //save_spots << "INTERMEDIATE: " << setprecision(15) << scientific << spots_to_Vector(spots) << endl;
02551                 }
02552             }
02553         }
02554 
02555         //spot_intensities is be correct here!
02556         
02557         for(int i=0; i < add_remove_tries; i++)
02558         {
02559             cout << endl << endl << "Modifying the model" << endl << "======================\n";
02560             cout << "Hello\n";
02561             bool add_spot = (rng() > 0.5) || (spots.size() == 1);
02562             cout << "World\n";
02563 
02564             vector<Vector<4> > model_1, model_2;
02565 
02566             
02567             int original; //What is the original model? Model 1 or Model 2?
02568 
02569             if(add_spot)
02570             {
02571                 model_1 = spots;
02572                 model_2 = model_1;
02573                 
02574                 //Pick a pixel within the thresholded ones as a starting point
02575                 int r;
02576                 do
02577                 {
02578                     r = (int)(rng() * pixels.size());
02579                     cout << r << " " << log_ratios[pixels[r]] << " " << pixels[r] << " " << threshold << endl;
02580                 }
02581                 while(log_ratios[pixels[r]] < threshold);
02582 
02583 
02584                 model_2.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), 
02585                                              log_normal_mode(blur_mu, blur_sigma), 
02586                                              pixels[r].x + rng()-.5, pixels[r].y + rng() - .5));
02587                 cout << "Adding a spot\n";
02588 
02589                 original = 1;
02590             }
02591             else
02592             {   
02593                 //Pick a random point to optimize/remove
02594                 int a_random_spot = static_cast<int>(rng() * spots.size());
02595                 model_1 = spots;
02596                 swap(model_1[model_1.size()-1], model_1[a_random_spot]);
02597 
02598                 model_2 = model_1;
02599 
02600                 model_1.pop_back();
02601                 cout << "Removing a spot\n";
02602                 original = 2;
02603             }
02604             
02605             //The mobile spot is always the last spot of model_2
02606             int spot = model_2.size() - 1;
02607 
02608             cout << "Original model: " << original << endl;
02609 
02610             //Precompute the intensities for all spot pixels
02611             //Model 2 is always one longer than model 1 and 
02612             //differs only on the extra element
02613             vector<vector<double> > model2_spot_intensities; //[spot][pixel]
02614             for(unsigned int i=0; i < model_2.size(); i++)
02615                 model2_spot_intensities.push_back(compute_spot_intensity(pixels, model_2[i]));
02616 
02617             //Which pixels does each spot have?
02618             vector<vector<int> > model2_spot_pixels(model_2.size()); //[spot][pixel]
02619             for(unsigned int s=0; s < model_2.size(); s++)
02620                 get_spot_pixels(pixels, model_2[s], model2_spot_pixels[s]);
02621 
02622             //Optimize spot:
02623             {
02624                 cout << "Optimizing spot for model selection\n";
02625 
02626 
02627                 //Get some samples with Gibbs sampling
02628                 vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
02629                 vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
02630 
02631                 GibbsSampler2 sampler(pixel_intensities, model2_spot_intensities, model_2, model2_spot_pixels, A, pi, variance, sample_iterations);
02632                 for(int i=0; i < add_remove_opt_samples; i++)
02633                 {
02634                     sampler.next(rng);
02635                     sample_list.push_back(sampler.sample());
02636                     sample_intensities.push_back(sampler.sample_intensities());
02637                 }
02638 
02639                 //First, remove the spot from all the samples.
02640                 for(unsigned int i=0; i < sample_list.size(); i++)
02641                     remove_spot(sample_intensities[i], model2_spot_intensities[spot], sample_list[i][spot]);
02642 
02643                 //Package up all the data
02644                 SampledBackgroundData data(sample_intensities, pixel_intensities, pixels, 
02645                                            intensity_mu, intensity_sigma, blur_mu, blur_sigma, 
02646                                            A, pi, variance);
02647                 
02648                 //Derivative computer:
02649                 SpotNegProbabilityDiffWithSampledBackground compute_deriv(data);
02650                 
02651                 graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), spot);
02652                 graphics.swap();
02653 
02654                 //Optimize spot "s"
02655                 //Optimize with the derivatives only since the actual probability
02656                 //is much harder to compute
02657                 ConjugateGradientOnly<4> cg(model_2[spot], compute_deriv, limit);
02658 
02659 
02660                 for(int attempt=0; attempt < add_remove_opt_retries; attempt++)
02661                 {
02662                     cout << "Attempt " << attempt << " of " << add_remove_opt_retries << endl;
02663 
02664                     //Optimize with conjugate gradient
02665                     while(cg.iterate(compute_deriv))
02666                     {
02667                         graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), spot, cg.x);
02668                         graphics.swap();
02669                     }
02670                     
02671                     //Check for being at a saddle point (no point checking on the last try)
02672                     //All eigenvectors should be negative at a maximum.
02673                     //WTF: is this a bug? WTF?
02674                     if(attempt < add_remove_tries - 1)
02675                     {
02676                         Matrix<4> hessian = sampled_background_spot_hessian_ffbs(cg.x, data, add_remove_opt_hess_inner_samples, rng);
02677                         SymEigen<4> hess_decomp(hessian);
02678 
02679                         //cout << "What the fuck:" << endl << hessian << endl << hessian3<< endl << hessian2 << endl;
02680                         
02681                         cout << "Eigenvalues are: " << hess_decomp.get_evalues() << endl;
02682 
02683                         if(hess_decomp.is_negdef())
02684                             break;
02685                         else
02686                         {
02687                             //Restart in the direction of the best uphill part
02688                             cg.init(cg.x + 0.1 * hess_decomp.get_evectors()[3], (hess_decomp.get_evectors()[3]));
02689 
02690 
02691                             cout << "Grad = " << compute_deriv(cg.x) << endl;
02692                             for(int i=0; i < 4; i++)
02693                             {
02694                                 cout << "Direction: " << i << endl;
02695                                 cout << unit(compute_deriv(cg.x + 0.1*hess_decomp.get_evectors()[i])) * hess_decomp.get_evectors()[i] << endl;
02696                             }
02697 
02698                             for(int i=0; i < 4; i++)
02699                             {
02700                                 cout << "Direction: " << i << endl;
02701                                 Vector<4> d = Zeros;
02702                                 d[i] = 1;
02703                                 cout << unit(compute_deriv(cg.x + d)) * hess_decomp.get_evectors()[i] << endl;
02704                             }
02705                         }
02706                     }
02707                 }
02708 
02709 
02710                 //Update to use the result of the optimization
02711                 model_2[spot] = cg.x;
02712                 
02713                 graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), -1);
02714                 graphics.swap();
02715     
02716 
02717                 //Update cached data based on the new spot position
02718                 model2_spot_intensities[spot] = compute_spot_intensity(pixels, model_2[spot]);
02719                 get_spot_pixels(pixels, model_2[spot], model2_spot_pixels[spot]);
02720 
02721                 cout << "Done optimizing for model selection\n";
02722             }
02723             
02724 
02725             //Which model to keep?
02726             int keep=original;
02727 
02728             //Integrate model_2
02729             //First compute the Hessian since this might go wrong.
02730 
02731             //FreeEnergyHessian hesscomputer(data_for_h_mcmc);
02732             Matrix<4> hess;// = hesscomputer.hessian(model_2, spot);
02733 
02734             //Use turbohess here since it is much faster, as the backwards sampling step is fast
02735             //We expect this hessian to be really quite different, if the spot has moved from
02736             //a long way away, since the sampling will have changed dramatically
02737             {
02738                 //Get some samples with Gibbs sampling
02739                 vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
02740                 vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
02741 
02742                 GibbsSampler sampler(pixel_intensities, model2_spot_intensities, model_2, A, pi, variance, sample_iterations);
02743                 for(int i=0; i < h_outer_samples; i++)
02744                 {
02745                     sampler.next(rng);
02746                     sample_list.push_back(sampler.sample());
02747                     sample_intensities.push_back(sampler.sample_intensities());
02748                 }
02749                 
02750                 //First, remove the spot from all the samples.
02751                 for(unsigned int i=0; i < sample_list.size(); i++)
02752                     remove_spot(sample_intensities[i], model2_spot_intensities[spot], sample_list[i][spot]);
02753 
02754                 //Package up all the data
02755                 SampledBackgroundData data(sample_intensities, pixel_intensities, pixels, 
02756                                            intensity_mu, intensity_sigma, blur_mu, blur_sigma, 
02757                                            A, pi, variance);
02758 
02759                 hess = sampled_background_spot_hessian_ffbs(model_2[spot], data, h_inner_samples, rng);
02760             }
02761 
02762 
02763             double det = determinant(-hess / (M_PI*2));
02764             SymEigen<4> hess_decomp(-hess);
02765             cout << "Hessien Eigenvalues are: " << hess_decomp.get_evalues() << endl;
02766             const double smallest_evalue = 1e-6;
02767 
02768             //If the determinant is negative, then we are still at a saddle point
02769             //despite the attempts above, so abandon the attempt.
02770             if(hess_decomp.get_evalues()[0] >  smallest_evalue)
02771             {
02772 
02773                 //Compute the free energy of model 2 at the MLE estimate
02774                 cout << "Model 2:\n";
02775         //      double model_2_energy = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_2));
02776                 double model_2_energy = -NegativeFreeEnergy(data_for_t_mcmc).compute_with_mask(spots_to_Vector(model_2), model2_spot_pixels);
02777                 cout << "Energy: " << model_2_energy << endl;
02778             
02779                 //Combine the MLE energy and Hessian using Laplace's approximation
02780                 double model_2_occam =  -0.5 * log(det);
02781                 double model_2_prob = model_2_energy + model_2_occam;
02782 
02783                 cout << "Occam: " << model_2_occam << endl;
02784                 cout << "Prob: " << model_2_prob << endl;
02785 
02786                 //Integrate model_1
02787                 //It has no parameters, in this formulation.
02788                 //double model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_1));
02789                 
02790                 //Note that model_1 always has one fewer spots, and the last spot is always the one
02791                 //missing, so we can make the correct mask very easily:
02792                 model2_spot_pixels.pop_back();
02793                 double model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc).compute_with_mask(spots_to_Vector(model_1), model2_spot_pixels);
02794                 cout << "Prob: " << model_1_prob << endl;
02795                 //model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_1));
02796 
02797                 if(model_2_prob > model_1_prob)
02798                     keep=2;
02799                 else
02800                     keep=1;
02801 
02802                 cout << "Models evaluated\n";
02803             }
02804             else
02805             {
02806                 cout << "Determinant has bad eigenvalues!\n";
02807                 keep = original;
02808                 cout << hess_decomp.get_evalues() << endl;
02809             }
02810 
02811             if(keep == 2)
02812             {
02813                 spots = model_2;
02814                 cout << "Keeping model 2\n";
02815 
02816             }
02817             else
02818             {
02819                 spots = model_1;
02820                 cout << "Keeping model 1\n";
02821             }
02822 
02823             if(original != keep)
02824             {
02825                 cout << "Model changed!\n";
02826                 //break;
02827             }
02828         }
02829     }
02830     save_spots << "FINAL: " << setprecision(15) << scientific << spots_to_Vector(spots) << endl;
02831 }

Generated on Wed Nov 2 18:00:00 2011 for BCUBED by  doxygen 1.6.3