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< FitSpotsGraphics > | null_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 > | |
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) |
Fit spots to the data.
Definition in file multispot5.cc.
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.
spot | Spot parameters | |
d | Background brightness (from other spots) | |
bs_iterations | Exter backward sampling iterations | |
rng | Random number generator |
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.
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.
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 }