ThreeB 1.1

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

- File List
- File Members

multispot5.cc
Go to the documentation of this file.

    00001 /*
    00002     This file is part of B-cubed.
    00003
    00004     Copyright (C) 2009, 2010, 2011, Edward Rosten and Susan Cox
    00005
    00006     B-cubed is free software; you can redistribute it and/or
    00007     modify it under the terms of the GNU Lesser General Public
    00008     License as published by the Free Software Foundation; either
    00009     version 3.0 of the License, or (at your option) any later version.
    00010
    00011     This library is distributed in the hope that it will be useful,
    00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
    00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    00014     Lesser General Public License for more details.
    00015
    00016     You should have received a copy of the GNU General Public License
    00017     along with this program.  If not, see <http://www.gnu.org/licenses/>
    00018 */
    00019
    00020 #include <cstdlib>
    00021 #include <cerrno>
    00022 #include <cstring>
    00023 #include <stack>
    00024 #include <algorithm>
    00025 #include <climits>
    00026 #include <iomanip>
    00027 #include <map>
    00028 #include <tr1/memory>
    00029 #include <cvd/image_io.h>
    00030 #include <cvd/image_convert.h>
    00031 #include <cvd/morphology.h>
    00032 #include <cvd/connected_components.h>
    00033 #include <cvd/draw.h>
    00034 #include <cvd/vector_image_ref.h>
    00035
    00036 #include <cvd/random.h>
    00037 #include <cvd/timer.h>
    00038 #include <gvars3/instances.h>
    00039
    00040 #include <TooN/functions/derivatives.h>
    00041 #include <TooN/determinant.h>
    00042 #include <TooN/SymEigen.h>
    00043 #include <TooN/optimization/conjugate_gradient.h>
    00044
    00045 #include "conjugate_gradient_only.h"
    00046 #include "forward_algorithm.h"
    00047 #include "numerical_derivatives.h"
    00048 #include "storm.h"
    00049 #include "storm_imagery.h"
    00050 #include "debug.h"
    00051 #include "sampled_multispot.h"
    00052 #include "mt19937.h"
    00053 #include "utility.h"
    00054 #include "multispot5.h"
    00055
    00056
    00057 //For benchmarking...
    00058 #define TIME(X)
    00059 //#define TIME(X) X
    00060
    00061 using namespace std;
    00062 using namespace CVD;
    00063 using namespace GVars3;
    00064 using namespace TooN;
    00065 using namespace std::tr1;
    00066
    00067 ///Empty destructor
    00068 UserInterfaceCallback::~UserInterfaceCallback(){}
    00069
    00070 ///User interface callback class which does nothing.
    00071 class NullUICallback: public UserInterfaceCallback
    00072 {
    00073     void per_spot(int, int, int, int){};
    00074     void per_modification(int, int, int){};
    00075     void per_pass(int, int, const std::vector<TooN::Vector<4> >&){};
    00076     void perhaps_stop(){};
    00077 };
    00078
    00079 ///Factory function to generate an instance of NullGraphics
    00080 ///@ingroup gStorm
    00081 auto_ptr<UserInterfaceCallback> null_ui()
    00082 {
    00083     return auto_ptr<UserInterfaceCallback>(new NullUICallback);
    00084 }
    00085
    00086
    00087 //Declare the graphics classes.
    00088 //These provide all the debug drawing operations so that the code can run in GUI and
    00089 //headless mode easily and without macros.
    00090 FitSpotsGraphics::~FitSpotsGraphics(){}
    00091
    00092 ///Graphics class which does absoloutely nothing
    00093 ///@ingroup gStorm
    00094 class NullGraphics: public  FitSpotsGraphics
    00095 {
    00096     public:
    00097         virtual void init(CVD::ImageRef){}
    00098         virtual void draw_krap(const std::vector<TooN::Vector<4> >&, const CVD::Image<CVD::byte>&, const BBox&, int, TooN::Vector<4>){}
    00099         virtual void swap(){}
    00100         virtual void draw_pixels(const std::vector<CVD::ImageRef>&, float, float, float, float){}
    00101         virtual void draw_bbox(const BBox&){}
    00102         virtual void glDrawCross(const TooN::Vector<2>&, int){}
    00103         virtual ~NullGraphics(){}
    00104 };
    00105
    00106 ///Factory function to generate an instance of NullGraphics
    00107 ///@ingroup gStorm
    00108 auto_ptr<FitSpotsGraphics> null_graphics()
    00109 {
    00110     return auto_ptr<FitSpotsGraphics>(new NullGraphics);
    00111 }
    00112
    00113
    00114 ///There are two sensible ways of storing the state vector of spot positions.
    00115 ///This function converts between them. See also spots_to_vector.
    00116 ///@param s list of spots to convert
    00117 ///@ingroup gUtility
    00118 Vector<> spots_to_Vector(const vector<Vector<4> >& s)
    00119 {
    00120     Vector<> r(s.size()*4);
    00121     for(unsigned int i=0; i < s.size(); i++)
    00122     {
    00123         r[i*4+0] = s[i][0];
    00124         r[i*4+1] = s[i][1];
    00125         r[i*4+2] = s[i][2];
    00126         r[i*4+3] = s[i][3];
    00127     }
    00128     return r;
    00129 }
    00130
    00131 ///There are two sensible ways of storing the state vector of spot positions.
    00132 ///This function converts between them. See also spots_to_Vector.
    00133 ///@param s list of spots to convert
    00134 ///@ingroup gUtility
    00135 vector<Vector<4> > spots_to_vector(const Vector<>& s)
    00136 {
    00137     vector<Vector<4> > r(s.size()/4);
    00138     for(unsigned int i=0; i < r.size(); i++)
    00139         r[i] = s.slice<Dynamic, 4>(i*4, 4);
    00140     return r;
    00141 }
    00142
    00143 ///Normalize an image for display purposes.
    00144 ///@ingroup gUtility
    00145 Image<byte> scale_to_bytes(const Image<float>& im, float lo, float hi)
    00146 {
    00147         Image<byte> out(im.size());
    00148         for(int r=0; r < out.size().y-0; r++)
    00149                 for(int c=0; c < out.size().x-0; c++)
    00150                         out[r][c] = (int)floor((im[r][c]-lo)*255/(hi-lo));
    00151         return out;
    00152 }
    00153
    00154 ///Normalize an image for display purposes.
    00155 ///@ingroup gUtility
    00156 Image<byte> scale_to_bytes(const Image<float>& im)
    00157 {
    00158         float lo = *min_element(im.begin(), im.end());
    00159         float hi = *max_element(im.begin(), im.end());
    00160         Image<byte> out(im.size());
    00161         for(int r=0; r < out.size().y-0; r++)
    00162                 for(int c=0; c < out.size().x-0; c++)
    00163                         out[r][c] = (int)floor((im[r][c]-lo)*255/(hi-lo));
    00164
    00165         return out;
    00166 }
    00167
    00168 ///Average the input image stack for display purposes
    00169 ///@ingroup gUtility
    00170 Image<float> average_image(const vector<Image<float> >& ims)
    00171 {
    00172     assert_same_size(ims);
    00173     Image<float> r(ims[0].size(), 0);
    00174
    00175     for(unsigned int i=0; i < ims.size(); i++)
    00176         transform(r.begin(), r.end(), ims[i].begin(), r.begin(), plus<float>());
    00177
    00178     transform(r.begin(), r.end(), r.begin(), bind2nd(multiplies<float>(), 1./ims.size()));
    00179     return r;
    00180 }
    00181
    00182
    00183 ////////////////////////////////////////////////////////////////////////////////
    00184 ////////////////////////////////////////////////////////////////////////////////
    00185 ////////////////////////////////////////////////////////////////////////////////
    00186 ////////////////////////////////////////////////////////////////////////////////
    00187 ////////////////////////////////////////////////////////////////////////////////
    00188 ////////////////////////////////////////////////////////////////////////////////
    00189 ////////////////////////////////////////////////////////////////////////////////
    00190 ////////////////////////////////////////////////////////////////////////////////
    00191
    00192 ///Closure hoding the data required do use GibbsSampler2
    00193 ///See FitSpots for naming of variables.
    00194 ///@ingroup gStorm
    00195 class DataForMCMC
    00196 {
    00197     protected:
    00198     const vector<ImageRef>& pixels;
    00199     const vector<vector<double> >& pixel_intensities;//[frame][pixel]
    00200     const double mu_brightness, sigma_brightness, mu_blur, sigma_blur;
    00201     const double variance;
    00202     const int samples, sample_iterations;
    00203     const Matrix<3> A;
    00204     const Vector<3> pi;
    00205     MT19937& rng;
    00206     public:
    00207
    00208     MT19937& get_rng()const
    00209     {
    00210         return rng;
    00211     }
    00212
    00213     public:
    00214     DataForMCMC(const vector<ImageRef>& pixels_,
    00215                 const vector<vector<double> >& pixel_intensities_,
    00216                 double mu_brightness_,
    00217                 double sigma_brightness_,
    00218                 double mu_blur_,
    00219                 double sigma_blur_,
    00220                 double variance_,
    00221                 int samples_,
    00222                 int sample_iterations_,
    00223                 Matrix<3> A_,
    00224                 Vector<3> pi_,
    00225                 MT19937& rng_)
    00226     :pixels(pixels_),
    00227      pixel_intensities(pixel_intensities_),
    00228      mu_brightness(mu_brightness_),
    00229      sigma_brightness(sigma_brightness_),
    00230      mu_blur(mu_blur_),
    00231      sigma_blur(sigma_blur_),
    00232      variance(variance_),
    00233      samples(samples_),
    00234      sample_iterations(sample_iterations_),
    00235      A(A_),
    00236      pi(pi_),
    00237      rng(rng_)
    00238     {}
    00239 };
    00240
    00241 ///Class implementing the Kahan summation algorithm
    00242 ///to allow accurate summation of very large numbers of
    00243 ///doubles.
    00244 ///@ingroup gUtility
    00245 class Kahan{
    00246     private:
    00247         double y; ///< Input with the compensation removed. Temporary working space.
    00248         double c; ///< Running compenstation for low-order bits.
    00249         double t; ///< y + sum, which loses low order bits of y. Temporary working space.
    00250     public:
    00251         double sum; ///< running sum
    00252
    00253         Kahan()
    00254         :c(0),sum(0)
    00255         {}
    00256
    00257         /// Add a number to the running sum
    00258         /// @param i Number to add
    00259         void add(double i)
    00260         {
    00261             //y = input -c
    00262             y = i;
    00263             y-= c;
    00264
    00265             //t = sum + y
    00266             t = sum;
    00267             t += y;
    00268
    00269             //c = (t - sum) - y
    00270             //c = ((sum + y) - sum) - y)
    00271             c = t;
    00272             c -= sum;
    00273             c -= y;
    00274             sum = t;
    00275         }
    00276 };
    00277
    00278
    00279 ///Class for computing the negitve free energy using thermodynamic integration.
    00280 class NegativeFreeEnergy: public DataForMCMC
    00281 {
    00282     public:
    00283
    00284     ///@param d Necessary data
    00285     NegativeFreeEnergy(const DataForMCMC& d)
    00286     :DataForMCMC(d)
    00287     {
    00288     }
    00289
    00290     ///Give the noise variance given a sample number and growth parameters
    00291     ///@param sample Sample number
    00292     ///@param samples Total number of samples
    00293     ///@param base_sigma Starting value of sigme
    00294     ///@param scale_pow Exponent scaling
    00295     double variance_from_sample(double sample, double samples, double base_sigma, double scale_pow) const
    00296     {
    00297         double scale = pow(1.25, sample * 1. / samples * scale_pow);
    00298         double sigma = base_sigma * scale;
    00299         double new_variance = sq(sigma);
    00300
    00301         return new_variance;
    00302     }
    00303
    00304     ///Estimate free energy using the Slow Growth Thermodynamic Integration method given in
    00305     ///Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford. M. Neal, 1993
    00306     ///Except using a 5 point stencil instead of forward differenceing in Eq 6.30
    00307     ///@param spots list of spots
    00308     ///@param spot_pixels Mask around each spot to use to save on computation
    00309     double compute_with_mask(const Vector<>& spots, const vector<vector<int> >& spot_pixels) const
    00310     {
    00311         //Estimate free energy using the Slow Growth Thermodynamic Integration method given in
    00312         //Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford. M. Neal, 1993
    00313         //Except using a 5 point stencil instead of forward differenceing in Eq 6.30
    00314         double base_sigma = sqrt(variance); // should be 1
    00315         double scale_pow = 100.0;
    00316
    00317         const unsigned int nspots  = spots.size()/4;
    00318         const unsigned int nframes = pixel_intensities.size();
    00319         const unsigned int npixels = pixels.size();
    00320         assert(spots.size() %4 == 0);
    00321         assert(spot_pixels.size() == nspots);
    00322
    00323         //Compute the intensities and derivatives for all spot
    00324         vector<vector<double> > spot_intensity; //[spot][pixel]
    00325         for(unsigned int i=0; i < nspots; i++)
    00326             spot_intensity.push_back(compute_spot_intensity(pixels, spots.slice<Dynamic,4>(i*4,4)));
    00327
    00328         GibbsSampler2 sampler(pixel_intensities, spot_intensity, spots_to_vector(spots), spot_pixels, A, pi, variance, sample_iterations);
    00329
    00330         double sum = 0;
    00331         Kahan ksum;
    00332         for(int sample=0; sample < samples; sample++)
    00333         {
    00334             //Compute the positions of the surrounding steps
    00335             double var1 = variance_from_sample(sample-2, samples, base_sigma, scale_pow);
    00336             double var2 = variance_from_sample(sample-1, samples, base_sigma, scale_pow);
    00337             double var3 = variance_from_sample(sample+1, samples, base_sigma, scale_pow);
    00338             double var4 = variance_from_sample(sample+2, samples, base_sigma, scale_pow);
    00339
    00340             //Take a sample
    00341             sampler.set_variance(var2);
    00342             sampler.next(DataForMCMC::get_rng());
    00343
    00344             //Compute the SSD. This is fixed regardless of sigma.
    00345             double err_sum=0;
    00346             for(unsigned int frame=0; frame < nframes; frame++)
    00347                 for(unsigned int pixel=0; pixel < npixels; pixel++)
    00348                     err_sum -= sq(pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel]);
    00349
    00350             //Compute the derivative using a five point stencil
    00351             //This could be done in a better but less clear way
    00352             double e1 = err_sum / (2*var1) - npixels*nframes*::log(2*M_PI*var1)/2;
    00353             double e2 = err_sum / (2*var2) - npixels*nframes*::log(2*M_PI*var2)/2;
    00354             double e3 = err_sum / (2*var3) - npixels*nframes*::log(2*M_PI*var3)/2;
    00355             double e4 = err_sum / (2*var4) - npixels*nframes*::log(2*M_PI*var4)/2;
    00356             sum += (-e1 + 8*e2 - 8*e3 + e4)/12;
    00357
    00358             ksum.add(-e1/12);
    00359             ksum.add(8*e2/12);
    00360             ksum.add(-8*e3/12);
    00361             ksum.add(e4/12);
    00362         }
    00363
    00364         double log_final = (log(variance_from_sample(samples, samples, base_sigma, scale_pow)*2*M_PI)/2) * npixels * nframes;
    00365
    00366         double priors=0;
    00367         for(unsigned int i=0; i < nspots; i++)
    00368         {
    00369             priors += log_log_normal(spots[i*4+0], mu_brightness, sigma_brightness);
    00370             priors += log_log_normal(spots[i*4+1], mu_blur, sigma_blur);
    00371         }
    00372
    00373         /*cout << "Thermo:\n";
    00374         cout << "sum: " << sum -log_final << endl;
    00375         cout << "kah: " << ksum.sum -log_final << endl;
    00376         cout << "priors: " << priors + sum -log_final << endl;
    00377         cout << "   kah: " << priors + ksum.sum -log_final << endl;
    00378 */
    00379         //cout << log_final << endl;
    00380         //cout << sum + log_final << endl;
    00381
    00382
    00383         sampler.set_variance(variance);
    00384         return -(sum+priors - log_final);
    00385     }
    00386
    00387     ///Estimate free energy using the Slow Growth Thermodynamic Integration method given in
    00388     ///Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford. M. Neal, 1993
    00389     ///Except using a 5 point stencil instead of forward differenceing in Eq 6.30
    00390     ///@param spots list of spots
    00391     double operator()(const Vector<>& spots) const
    00392     {
    00393         double base_sigma = sqrt(variance); // should be 1
    00394         double scale_pow = 100.0;
    00395
    00396         const unsigned int nspots  = spots.size()/4;
    00397         const unsigned int nframes = pixel_intensities.size();
    00398         const unsigned int npixels = pixels.size();
    00399         assert(spots.size() %4 == 0);
    00400
    00401         //Compute the intensities and derivatives for all spot
    00402         vector<vector<double> > spot_intensity; //[spot][pixel]
    00403         for(unsigned int i=0; i < nspots; i++)
    00404             spot_intensity.push_back(compute_spot_intensity(pixels, spots.slice<Dynamic,4>(i*4,4)));
    00405
    00406         GibbsSampler sampler(pixel_intensities, spot_intensity, spots_to_vector(spots), A, pi, variance, sample_iterations);
    00407
    00408         double sum = 0;
    00409         Kahan ksum;
    00410         for(int sample=0; sample < samples; sample++)
    00411         {
    00412             //Compute the positions of the surrounding steps
    00413             double var1 = variance_from_sample(sample-2, samples, base_sigma, scale_pow);
    00414             double var2 = variance_from_sample(sample-1, samples, base_sigma, scale_pow);
    00415             double var3 = variance_from_sample(sample+1, samples, base_sigma, scale_pow);
    00416             double var4 = variance_from_sample(sample+2, samples, base_sigma, scale_pow);
    00417
    00418             //Take a sample
    00419             sampler.set_variance(var2);
    00420             sampler.next(DataForMCMC::get_rng());
    00421
    00422             //Compute the SSD. This is fixed regardless of sigma.
    00423             double err_sum=0;
    00424             for(unsigned int frame=0; frame < nframes; frame++)
    00425                 for(unsigned int pixel=0; pixel < npixels; pixel++)
    00426                     err_sum -= sq(pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel]);
    00427
    00428             //Compute the derivative using a five point stencil
    00429             //This could be done in a better but less clear way
    00430             double e1 = err_sum / (2*var1) - npixels*nframes*::log(2*M_PI*var1)/2;
    00431             double e2 = err_sum / (2*var2) - npixels*nframes*::log(2*M_PI*var2)/2;
    00432             double e3 = err_sum / (2*var3) - npixels*nframes*::log(2*M_PI*var3)/2;
    00433             double e4 = err_sum / (2*var4) - npixels*nframes*::log(2*M_PI*var4)/2;
    00434             sum += (-e1 + 8*e2 - 8*e3 + e4)/12;
    00435
    00436             ksum.add(-e1/12);
    00437             ksum.add(8*e2/12);
    00438             ksum.add(-8*e3/12);
    00439             ksum.add(e4/12);
    00440         }
    00441
    00442         double log_final = (log(variance_from_sample(samples, samples, base_sigma, scale_pow)*2*M_PI)/2) * npixels * nframes;
    00443
    00444         double priors=0;
    00445         for(unsigned int i=0; i < nspots; i++)
    00446         {
    00447             priors += log_log_normal(spots[i*4+0], mu_brightness, sigma_brightness);
    00448             priors += log_log_normal(spots[i*4+1], mu_blur, sigma_blur);
    00449         }
    00450
    00451         /*cout << "Thermo:\n";
    00452         cout << "sum: " << sum -log_final << endl;
    00453         cout << "kah: " << ksum.sum -log_final << endl;
    00454         cout << "priors: " << priors + sum -log_final << endl;
    00455         cout << "   kah: " << priors + ksum.sum -log_final << endl;
    00456 */
    00457         //cout << log_final << endl;
    00458         //cout << sum + log_final << endl;
    00459
    00460
    00461         sampler.set_variance(variance);
    00462         return -(sum+priors - log_final);
    00463     }
    00464 };
    00465
    00466 /// Class for sorting a list of indexes to an array of spots lexicographically
    00467 /// according to the 2D positions of the spots.
    00468 ///
    00469 ///@param Cmp comparator function to specify less or greater
    00470 ///@param First most significant position index
    00471 ///@ingroup gMultiSpot
    00472 template<class Cmp, int First> struct IndexLexicographicPosition{
    00473     const vector<Vector<4> >& spots;  ///< Keep around the array of spots for later comprison
    00474
    00475     /// @param s Vector to sort indices of
    00476     IndexLexicographicPosition(const vector<Vector<4> >& s)
    00477     :spots(s)
    00478     {}
    00479
    00480
    00481     static const int Second = First==2?3:2; ///< Second most siginifcant position index for sorting
    00482
    00483     /// Compare two indexes into the array of spots
    00484     bool operator()(int a, int b)
    00485     {
    00486         Cmp cmp;
    00487
    00488         if(cmp(spots[a][First], spots[b][First]))
    00489             return true;
    00490         else if(spots[a][First] == spots[b][First])
    00491             return cmp(spots[a][Second], spots[b][Second]);
    00492         else
    00493             return false;
    00494     }
    00495 };
    00496
    00497
    00498 ///Closure holding image data generated using samples drawn from the model.
    00499 ///NB this is used with one spot removed (i.e. set to dark). As a result, the data
    00500 ///is treated as the background when considering that spot in isolation.
    00501 ///See FitSpots for naming of variables.
    00502 ///@ingroup gStorm
    00503 struct SampledBackgroundData
    00504 {
    00505     const vector<vector<vector<double> > >& sample_intensities_without_spot;    //[sample][frame][pixel]
    00506     const vector<vector<double> >& pixel_intensities;    //[frame][pixel]
    00507     const vector<ImageRef> pixels;
    00508
    00509     double mu_brightness, sigma_brightness, mu_blur, sigma_blur;
    00510     const Matrix<3> A;
    00511     const Vector<3> pi;
    00512     double variance;
    00513
    00514     const vector<int> O;
    00515
    00516     SampledBackgroundData(
    00517         const vector<vector<vector<double> > >& sample_intensities_without_spot_,
    00518         const vector<vector<double> >& pixel_intensities_,
    00519         const vector<ImageRef> pixels_,
    00520         double mu_brightness_,
    00521         double sigma_brightness_,
    00522         double mu_blur_,
    00523         double sigma_blur_,
    00524         const Matrix<3> A_,
    00525         const Vector<3> pi_,
    00526         double variance_)
    00527     :sample_intensities_without_spot(sample_intensities_without_spot_),
    00528         pixel_intensities(pixel_intensities_),
    00529         pixels(pixels_),
    00530         mu_brightness(mu_brightness_),
    00531         sigma_brightness(sigma_brightness_),
    00532         mu_blur(mu_blur_),
    00533         sigma_blur(sigma_blur_),
    00534         A(A_),
    00535         pi(pi_),
    00536         variance(variance_),
    00537         O(sequence(pixel_intensities.size()))
    00538     {
    00539     }
    00540 };
    00541
    00542 //!@cond Doxygen_Suppress
    00543 //Do not use.
    00544 struct SpotProbabilityWithSampledBackgroundFAKE: public SampledBackgroundData
    00545 {
    00546     SpotProbabilityWithSampledBackgroundFAKE(const SampledBackgroundData& d)
    00547     :SampledBackgroundData(d)
    00548     {
    00549     }
    00550
    00551     //Compute the probability of spot
    00552     double operator()(const Vector<4>& spot) const
    00553     {
    00554         vector<double> spot_intensities = compute_spot_intensity(pixels, spot);
    00555
    00556         double sum_log_prob = 0;
    00557
    00558         for(unsigned int s=0; s < sample_intensities_without_spot.size(); s++)
    00559         {
    00560             SpotWithBackground B(sample_intensities_without_spot[s], spot_intensities, pixel_intensities, variance);
    00561             double log_prob = forward_algorithm(A, pi, B, O);
    00562
    00563             double logprior = log_log_normal(spot[0], mu_brightness, sigma_brightness) +
    00564                           log_log_normal(spot[1], mu_blur, sigma_blur);
    00565
    00566             //cout << setprecision(10) <<fixed <<  "Prob : " << log_prob + logprior << endl;
    00567
    00568             sum_log_prob += log_prob + logprior;
    00569         }
    00570
    00571         double average_log_prob = sum_log_prob / sample_intensities_without_spot.size();
    00572
    00573     //  cout << "Ave Prob : " << average_log_prob << endl;
    00574
    00575         if(spot[0] < 0 || spot[1] < 0)
    00576             return 1e100;
    00577
    00578         return average_log_prob;
    00579     }
    00580 };
    00581 //@endcond
    00582
    00583 ///Compute the derivative of the negative log probability with respect
    00584 ///to the parameters of one spot, given some samples of the other spots.
    00585 ///@ingroup gStorm
    00586 struct SpotNegProbabilityDiffWithSampledBackground: public SampledBackgroundData
    00587 {
    00588     ///@param d Necessary data for construction
    00589     SpotNegProbabilityDiffWithSampledBackground(const SampledBackgroundData& d)
    00590     :SampledBackgroundData(d)
    00591     {
    00592     }
    00593
    00594     ///Compute the probability of spot
    00595     ///@param spot Spot position
    00596     Vector<4> operator()(const Vector<4>& spot) const
    00597     {
    00598         if(spot[0] <= 0 || spot[1] <= 0)
    00599             return Ones * std::numeric_limits<double>::quiet_NaN();
    00600
    00601         vector<pair<double, Vector<4> > > spot_intensities = compute_spot_intensity_derivatives(pixels, spot);
    00602
    00603         Vector<4> sum_diff_log = Zeros;
    00604
    00605         for(unsigned int s=0; s < sample_intensities_without_spot.size(); s++)
    00606         {
    00607             SpotWithBackground B(sample_intensities_without_spot[s], spot_intensities, pixel_intensities, variance);
    00608
    00609             pair<double, Vector<4> > r = forward_algorithm_deriv(A, pi, B, O);
    00610
    00611             sum_diff_log += r.second;
    00612         }
    00613
    00614         Vector<4> diff_log = sum_diff_log / sample_intensities_without_spot.size();
    00615
    00616         //Compute the log probability of the prior
    00617         Vector<4> logprior_deriv = makeVector(diff_log_log_normal(spot[0], mu_brightness, sigma_brightness),
    00618                                               diff_log_log_normal(spot[1], mu_blur, sigma_blur), 0, 0);
    00619
    00620         return -(diff_log + logprior_deriv);
    00621     }
    00622 };
    00623
    00624
    00625 ///Class for computing the Hessian of the negative free energy.
    00626 ///@ingroup gStorm
    00627 class FreeEnergyHessian: public DataForMCMC
    00628 {
    00629     public:
    00630
    00631     ///Constructor
    00632     ///@param d All data required
    00633     FreeEnergyHessian(const DataForMCMC& d)
    00634     :DataForMCMC(d)
    00635     {
    00636     }
    00637
    00638     ///Compute the Hessian
    00639     ///@param spots All spot positions
    00640     ///@param spot spot to compute hessian for
    00641     Matrix<4> hessian(const vector<Vector<4> >& spots, int spot) const
    00642     {
    00643         cout << "----Computing pure MCMC hessian\n";
    00644         const unsigned int nspots  = spots.size();
    00645         const unsigned int nframes = pixel_intensities.size();
    00646         const unsigned int npixels = pixels.size();
    00647         cout << spot << " " << nspots << " " << nframes << " " << npixels << endl;
    00648
    00649         vector<vector<double> > spot_intensity; //[spot][pixel]
    00650         for(unsigned int i=0; i < nspots; i++)
    00651             spot_intensity.push_back(compute_spot_intensity(pixels, spots[i]));
    00652
    00653         vector<tuple<double, Vector<4>, Matrix<4> > > spot_hess_etc = compute_spot_intensity_hessian(pixels, spots[spot]);
    00654
    00655         GibbsSampler sampler(pixel_intensities, spot_intensity, spots, A, pi, variance, sample_iterations);
    00656
    00657         //Compute derivative of log probability, summed (ie integrated)
    00658         Matrix<4> sum_hess1 = Zeros(spots.size());
    00659         Matrix<4> sum_hess2 = Zeros(spots.size());
    00660         Vector<4> sum_deriv = Zeros(spots.size());
    00661
    00662         for(int sample=0; sample < samples; sample++)
    00663         {
    00664             sampler.next(rng);
    00665
    00666             //Compute d log P(data | x, model) / d model, for a given sample
    00667             //And the hessian
    00668             Matrix<4> hess = Zeros(spots.size());
    00669             Vector<4> deriv = Zeros(spots.size());
    00670             for(unsigned int frame=0; frame < nframes; frame++)
    00671             {
    00672                 for(unsigned int pixel=0; pixel < npixels; pixel++)
    00673                 {
    00674                     double e = pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel];
    00675                     //Build up the derivative
    00676                     if(sampler.sample()[spot][frame] == 0)
    00677                     {
    00678                         hess += e * get<2>(spot_hess_etc[pixel]) - get<1>(spot_hess_etc[pixel]).as_col() * get<1>(spot_hess_etc[pixel]).as_row();
    00679                         deriv += e * get<1>(spot_hess_etc[pixel]);
    00680
    00681                     }
    00682                 }
    00683             }
    00684
    00685             hess[0][0] += hess_log_log_normal(spots[spot][0], mu_brightness, sigma_brightness);
    00686             hess[1][1] += hess_log_log_normal(spots[spot][1], mu_blur, sigma_blur);
    00687             sum_hess1 += hess;
    00688
    00689             deriv[0] += diff_log_log_normal(spots[spot][0], mu_brightness, sigma_brightness);
    00690             deriv[1] += diff_log_log_normal(spots[spot][1], mu_blur, sigma_blur);
    00691
    00692             sum_hess2 += deriv.as_col() * deriv.as_row();
    00693
    00694             sum_deriv = sum_deriv + deriv;
    00695         }
    00696
    00697         sum_hess1 /= (samples * variance);
    00698         sum_hess2 /= (samples * variance);
    00699         sum_deriv /= (samples * variance);
    00700
    00701
    00702         cout << sum_hess1 << endl;
    00703         cout << sum_hess2 << endl;
    00704         cout << sum_deriv.as_col() * sum_deriv.as_row() <<  endl;
    00705
    00706         cout << "......." << sum_deriv << endl;
    00707         //Put in the prior
    00708         //The derivative prior parts cancel out.
    00709         //Rather sensibly this means that the second derivatives can be
    00710         //computed without reference to the prior, and then the prior can
    00711         //be added in later, i.e.: P(data, parameters) = P(data|parameter) P(parameters)
    00712
    00713         //The second derivatives have been constructed to be diagonal
    00714         DiagonalMatrix<4> hess_prior = Zeros(spots.size());
    00715
    00716         cout << "sum of parts = \n" << sum_hess1 + sum_hess2 - sum_deriv.as_col() * sum_deriv.as_row() << endl;
    00717         //TooN cannot currently add DiagonalMatrix to Matrix!!
    00718         //sum_hess.diagonal_slice() += hess_prior.diagonal_slice();
    00719
    00720         cout << "++++Done Computing pure MCMC hessian\n";
    00721         return sum_hess1 + sum_hess2 - sum_deriv.as_col() * sum_deriv.as_row();
    00722     }
    00723 };
    00724
    00725 ///Comparator functor for the first element of a std::pair
    00726 struct LessSecond
    00727 {
    00728     ///Comparison function
    00729     ///@param a
    00730     ///@param b
    00731     template<class A, class B> bool operator()(const pair<A,B>& a, const pair<A,B>& b) const
    00732     {
    00733         return a.second < b.second;
    00734     }
    00735 };
    00736
    00737 //!@cond Doxygen_Suppress
    00738 struct NthDeriv{
    00739
    00740     const SpotNegProbabilityDiffWithSampledBackground& compute_deriv;
    00741     int i;
    00742
    00743     NthDeriv(const SpotNegProbabilityDiffWithSampledBackground& c, int ii)
    00744     :compute_deriv(c),i(ii)
    00745     {}
    00746
    00747     double  operator()(const Vector<4>& f) const
    00748     {
    00749         return compute_deriv(f)[i];
    00750     }
    00751 };
    00752 //!@endcond
    00753
    00754 ///Compute the Hessian of the log probability. The background is sampled rather sparsely, and the
    00755 ///spot in question is sampled much more densely using FFBS.
    00756 ///@param spot Spot parameters
    00757 ///@param d Background brightness (from other spots)
    00758 ///@param bs_iterations Exter backward sampling iterations
    00759 ///@param rng Random number generator
    00760 ///@returns the Hessian of the log probability around the spot
    00761 Matrix<4> sampled_background_spot_hessian_ffbs(const Vector<4>& spot, const SampledBackgroundData& d, int bs_iterations, MT19937& rng)
    00762 {
    00763     vector<tuple<double, Vector<4>, Matrix<4> > > spot_hess_etc = compute_spot_intensity_hessian(d.pixels, spot);
    00764     vector<double> spot_intensities = compute_spot_intensity(d.pixels, spot);
    00765
    00766     Matrix<4> sum_hess_log = Zeros;
    00767     Matrix<4> sum_diff2_log = Zeros;
    00768
    00769     vector<State> current_sample;
    00770
    00771     const unsigned int nframes = d.pixel_intensities.size();
    00772     const unsigned int npixels = d.pixels.size();
    00773
    00774     Matrix<4> sum_hess  = Zeros;
    00775     Vector<4> sum_deriv = Zeros;
    00776
    00777     vector<pair<Matrix<4>, Vector<4> > > hess_and_deriv_part(nframes);
    00778
    00779     for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++)
    00780     {
    00781         SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensities, d.pixel_intensities, d.variance);
    00782
    00783         //Compute what the per-frame hess and deriv parts are
    00784         //if the spot is on in a frame.
    00785         for(unsigned int frame=0; frame < nframes; frame++)
    00786         {
    00787             Matrix<4> hess = Zeros;
    00788             Vector<4> deriv = Zeros;
    00789
    00790             for(unsigned int pixel=0; pixel < npixels; pixel++)
    00791             {
    00792                 double e = d.pixel_intensities[frame][pixel] - (d.sample_intensities_without_spot[s][frame][pixel] + spot_intensities[pixel]);
    00793                 //Build up the derivative
    00794                 hess += e * get<2>(spot_hess_etc[pixel]) - get<1>(spot_hess_etc[pixel]).as_col() * get<1>(spot_hess_etc[pixel]).as_row();
    00795                 deriv += e * get<1>(spot_hess_etc[pixel]);
    00796             }
    00797             hess_and_deriv_part[frame] = make_pair(hess, deriv);
    00798         }
    00799
    00800         //Forward filtering
    00801         std::vector<array<double, 3> > delta = forward_algorithm_delta(d.A, d.pi, B, d.O);
    00802
    00803         for(int i=0; i < bs_iterations; i++)
    00804         {
    00805             current_sample = backward_sampling<3,State>(d.A, delta, rng);
    00806
    00807             Matrix<4> hess = Zeros;
    00808             Vector<4> deriv = Zeros;
    00809             for(unsigned int frame=0; frame < nframes; frame++)
    00810                 if(current_sample[frame] == 0)
    00811                 {
    00812                         hess += hess_and_deriv_part[frame].first;
    00813                         deriv += hess_and_deriv_part[frame].second;
    00814                 }
    00815
    00816             sum_hess += hess + deriv.as_col() * deriv.as_row();
    00817             sum_deriv += deriv;
    00818         }
    00819     }
    00820
    00821     sum_hess  /= (bs_iterations * d.sample_intensities_without_spot.size() * d.variance);
    00822     sum_deriv /= (bs_iterations * d.sample_intensities_without_spot.size() * d.variance);
    00823
    00824     sum_hess -= sum_deriv.as_col() * sum_deriv.as_row();
    00825
    00826     sum_hess[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness);
    00827     sum_hess[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);
    00828     sum_deriv[0] += diff_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness);
    00829     sum_deriv[1] += diff_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);
    00830
    00831     //cout << "Turboderiv:" << sum_deriv << endl;
    00832     //cout << "Turbohess:\n" << sum_hess << endl;
    00833
    00834     return sum_hess;
    00835 }
    00836
    00837 ///Debugging function. Not mathematically correct. Do not use.
    00838 Matrix<4> sampled_background_spot_hessian2(const Vector<4>& spot, const SampledBackgroundData& d)
    00839 {
    00840     vector<tuple<double, Vector<4>, Matrix<4> > > spot_intensities = compute_spot_intensity_hessian(d.pixels, spot);
    00841
    00842     Matrix<4> sum_hess_log = Zeros;
    00843     Matrix<4> sum_diff2_log = Zeros;
    00844
    00845     for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++)
    00846     {
    00847         SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensities, d.pixel_intensities, d.variance);
    00848
    00849         double prob;
    00850         Vector<4> diff;
    00851         Matrix<4> hess;
    00852
    00853         tie(prob, diff, hess) = forward_algorithm_hessian(d.A, d.pi, B, d.O);
    00854
    00855         sum_hess_log += hess;
    00856
    00857         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);
    00858         sum_diff2_log += diff.as_col() * diff.as_row();
    00859     }
    00860
    00861     Matrix<4> hess_log = sum_hess_log / d.sample_intensities_without_spot.size();
    00862     Matrix<4> diff2_log = sum_diff2_log / d.sample_intensities_without_spot.size();
    00863
    00864     //Add in the prior
    00865
    00866     hess_log[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness);
    00867     hess_log[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);
    00868
    00869     return hess_log + diff2_log;
    00870 }
    00871
    00872 ///Debugging function. Not mathematically correct. Do not use.
    00873 Matrix<4> sampled_background_spot_hessian_FAKE(const Vector<4>& spot, const SampledBackgroundData& d)
    00874 {
    00875     vector<tuple<double, Vector<4>, Matrix<4> > > spot_intensities = compute_spot_intensity_hessian(d.pixels, spot);
    00876
    00877     Matrix<4> sum_hess_log = Zeros;
    00878
    00879     for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++)
    00880     {
    00881         SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensities, d.pixel_intensities, d.variance);
    00882
    00883         double prob;
    00884         Vector<4> diff;
    00885         Matrix<4> hess;
    00886
    00887         tie(prob, diff, hess) = forward_algorithm_hessian(d.A, d.pi, B, d.O);
    00888
    00889         sum_hess_log += hess;
    00890     }
    00891
    00892     Matrix<4> hess_log = sum_hess_log / d.sample_intensities_without_spot.size();
    00893
    00894     //Add in the prior
    00895
    00896     hess_log[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness);
    00897     hess_log[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);
    00898
    00899     return hess_log;
    00900 }
    00901
    00902 ///Which pixels belong to a given spot?
    00903 ///Find the indices of those pixels
    00904 ///@ingroup gStorm
    00905 void get_spot_pixels(const vector<ImageRef>& pixels, const Vector<4>& spot, vector<int>& out)
    00906 {
    00907     //Go out to three sigma
    00908
    00909     vector<ImageRef> pix = getDisc(spot[1]*6 + 1);
    00910     out.resize(0);
    00911     ImageRef offset = ir_rounded(spot.slice<2,2>());
    00912     for(unsigned int j=0; j < pix.size(); j++)
    00913     {
    00914         int pos = lower_bound(pixels.begin(), pixels.end(), pix[j] + offset) - pixels.begin();
    00915         if(pos != (int)pixels.size() && pixels[pos] == pix[j] + offset)
    00916             out.push_back(pos);
    00917     }
    00918
    00919     if(out.size() == 0)
    00920     {
    00921         cout << "********************************\n";
    00922         cout << "********************************\n";
    00923         cout << "********************************\n";
    00924         cout << "********************************\n";
    00925         cout << "********************************\n";
    00926         cout << "Oe noes!11one\n";
    00927         cout << pix.size() << endl;
    00928     }
    00929 }
    00930
    00931
    00932 ///Tokenize a line
    00933 ///@ingroup gUtility
    00934 vector<string> split(const string& line)
    00935 {
    00936     vector<string> v;
    00937     istringstream i(line);
    00938     string s;
    00939
    00940     while(!i.eof())
    00941     {
    00942         i >> s;
    00943         if(i.fail())
    00944             break;
    00945         v.push_back(s);
    00946     }
    00947     return v;
    00948 }
    00949
    00950 ///Generic version of itoa.
    00951 ///How many times has this been reimplemented??
    00952 ///@param x Value to convert to string
    00953 ///@ingroup gUtility
    00954 template<class C> inline string xtoa(const C& x)
    00955 {
    00956     ostringstream os;
    00957     os << x;
    00958     return os.str();
    00959 }
    00960
    00961 ///Inverse of xtoa()
    00962 ///How many times has this been reimplemented??
    00963 ///@param s String to convert
    00964 ///@param msg Mesage to print on error
    00965 ///@ingroup gUtility
    00966 template<class C> inline C atox(const string& s, const string& msg)
    00967 {
    00968     C c;
    00969     istringstream i(s);
    00970     i >> c;
    00971     if(i.fail())
    00972         throw LogFileParseError("Error parsing " + msg + ". Text is `" + s + "'.");
    00973     return c;
    00974 }
    00975
    00976 /**Parser for multispot 5 log files.
    00977
    00978 Log files are mostly line oriented and contain various records
    00979
    00980 The main records are:
    00981
    00982 Iteraton: \#ITNUM
    00983 MAIN: <list of spot parameters>
    00984
    00985 Pass: \#PASSNUM
    00986 MT19937 <random number generator state>
    00987 PASS\#PASSNUM: <list of spot parameters>
    00988 ENDCHECKPOINT
    00989
    00990 Note that MAIN is redundant since it will be the same as the following PASS 1 (or the first pass
    00991 computed if restoring from a checkpoint).
    00992
    00993 Data should only be considered valid after ENDCHECKPOINT has been read
    00994
    00995 Iteration is written once per iteration, not once per pass. (FIXME)
    00996
    00997 Which moron invented this file format???
    00998
    00999 Note that the file format hasn't beren fixed, so that the output can easily be compared to the
    01000 output of the historic version which is known to be good.
    01001
    01002 @param in Stream to parse file from
    01003 @ingroup gStorm
    01004 */
    01005 StateParameters parse_log_file(istream& in)
    01006 {
    01007     //A line read from the file
    01008     string line;
    01009
    01010     //State lines known to be OK
    01011     string rngline, passline, iterationline;
    01012     bool state_ok=0;
    01013
    01014     //State lines read in, with flags of goodness
    01015     string new_rngline, new_passline, new_iterationline;
    01016     bool new_rngline_ok=0, new_passline_ok=0, new_iterationline_ok=0;
    01017
    01018     unsigned int lineno=0;
    01019     bool doing_gvars = 0;
    01020
    01021     vector<ImageRef> pixels;
    01022
    01023     while(!in.eof())
    01024     {
    01025         getline(in, line);
    01026         if(in.fail())
    01027             break;
    01028
    01029         lineno++;
    01030
    01031         if(line == "ENDGVARLIST")
    01032         {
    01033             if(!doing_gvars)
    01034                 throw LogFileParseError("Spurious end of GVars");
    01035             doing_gvars = 0;
    01036         }
    01037         else if(doing_gvars)
    01038         {
    01039             GUI.ParseLine(line);
    01040         }
    01041         else if(line == "BEGINGVARLIST")
    01042         {
    01043             doing_gvars = 1;
    01044         }
    01045         if(line.substr(0, 11)  == "Iteration: ")
    01046         {
    01047             new_iterationline = line;
    01048             new_iterationline_ok = true;
    01049         }
    01050         else if(line.substr(0, 4) == "PASS")
    01051         {
    01052             new_passline = line;
    01053             if(new_passline_ok)
    01054                 throw LogFileParseError("Duplicate PASS on line " + xtoa(lineno));
    01055             new_passline_ok = true;
    01056         }
    01057         else if(line.substr(0, 8) == "MT19937 ")
    01058         {
    01059             new_rngline = line;
    01060             if(new_rngline_ok)
    01061                 throw LogFileParseError("Duplicate MT19937 on line " + xtoa(lineno));
    01062
    01063             new_rngline_ok = true;
    01064
    01065         }
    01066         else if(line == "ENDCHECKPOINT")
    01067         {
    01068             if(new_passline_ok && new_rngline_ok && new_iterationline_ok)
    01069             {
    01070                 iterationline = new_iterationline;
    01071                 rngline = new_rngline;
    01072                 passline = new_passline;
    01073             }
    01074             else
    01075                 throw LogFileParseError("Reached checkpoint with missing data: "
    01076                          "it=" + xtoa(new_iterationline_ok) +
    01077                         " pa=" + xtoa(new_passline_ok) +
    01078                         " rg=" + xtoa(new_rngline_ok) + " on line " + xtoa(lineno));
    01079
    01080             //Don't reset iteration since it only appears once for each entire
    01081             //set of passes.
    01082             new_rngline_ok = 0;
    01083             new_passline_ok = 0;
    01084
    01085             state_ok = true;
    01086         }
    01087         else if(line.substr(0, 7) == "PIXELS ")
    01088         {
    01089             vector<string> l = split(line);
    01090             if( (l.size() - 1)%2 == 0)
    01091             {
    01092                 int n = (l.size()-1)/2;
    01093                 pixels.resize(n);
    01094                 for(int i=0; i < n; i++)
    01095                 {
    01096                     pixels[i].x = atox<int>(l[i*2+1+0], "pixels");
    01097                     pixels[i].y = atox<int>(l[i*2+1+1], "pixels");
    01098                 }
    01099             }
    01100             else
    01101                 throw LogFileParseError("Bad PIXELS line");
    01102         }
    01103     }
    01104
    01105     if(!state_ok)
    01106         throw LogFileParseError("No state found");
    01107
    01108     if(pixels.size() == 0)
    01109         throw LogFileParseError("No pixels, or pixels is empty");
    01110
    01111     //Now parse the lines
    01112     StateParameters p;
    01113     vector<string> l;
    01114
    01115     //Parse the iterations
    01116     l = split(iterationline);
    01117     p.iteration = atox<int>(l[1], "iteration");
    01118
    01119     //Parse the random number generator
    01120     p.rng =shared_ptr<MT19937>(new MT19937);
    01121     {
    01122         istringstream rng_s(rngline);
    01123         try{
    01124             p.rng->read(rng_s);
    01125         }
    01126         catch(MT19937::ParseError p)
    01127         {
    01128             throw LogFileParseError("Error parsing MT19937");
    01129         }
    01130     }
    01131
    01132     //Parse PASS and the listing of spots
    01133     l = split(passline);
    01134     if( (l.size() - 1 ) % 4 == 0)
    01135     {
    01136         p.pass = atox<int>(l[0].substr(4), "pass");
    01137
    01138         for(unsigned int i=0; i < (l.size()-1)/4; i++)
    01139         {
    01140             cout << l[i*4+1+0] << endl;
    01141             cout << l[i*4+1+1] << endl;
    01142             cout << l[i*4+1+2] << endl;
    01143             cout << l[i*4+1+3] << endl;
    01144             p.spots.push_back(makeVector(
    01145                         atox<double>(l[i*4+1+0], "spot"),
    01146                         atox<double>(l[i*4+1+1], "spot"),
    01147                         atox<double>(l[i*4+1+2], "spot"),
    01148                         atox<double>(l[i*4+1+3], "spot")));
    01149         }
    01150
    01151     }
    01152     else
    01153         throw LogFileParseError("Wrong number of elements in PASS line");
    01154
    01155     //Set up the pixels (oh god the pixels)
    01156     p.pixels = pixels;
    01157
    01158     return p;
    01159 }
    01160
    01161
    01162
    01163 ///Setup the parameters for a run using the old and deeply peculiar method.
    01164 ///This includes the unpleasant and diffucult so use de-checkpointing code.
    01165 ///wtf. The use of this function is very strongly deprecated.
    01166 ///@param log_ratios Image from which region is selected.
    01167 ///@param ims  Input data
    01168 ///@param pixels Region for spot fitting to run in
    01169 StateParameters generate_state_parameters_ye_olde(const BasicImage<double>& log_ratios, const vector<Image<float> >& ims, vector<ImageRef> pixels){
    01170     sort(pixels.begin(), pixels.end());
    01171
    01172     const double variance = 1; // it should be
    01173
    01174     //To scale the X axis of a log-normal distribution, only
    01175     //the mu parameter needs to be changed...
    01176     const double intensity_mu = GV3::get<double>("intensity.rel_mu", 0., -1)  + log(sqrt(variance));
    01177     const double intensity_sigma = GV3::get<double>("intensity.rel_sigma", 0., -1);
    01178     const double blur_mu = GV3::get<double>("blur.mu", 0., -1);
    01179     const double blur_sigma = GV3::get<double>("blur.sigma", 0., -1);
    01180
    01181     //The region was extracted at a certain threshold.
    01182     //These regions may be too small, so some post-region finding dilation
    01183     //may be performed. New spots are only placed at pixels which exceed the threshold.
    01184     //post_dilate.threshold is (if set) used as the placing threshold so the placing threshold
    01185     //can be different from the region-finding threshold.
    01186
    01187     //Note that as a result of dliation, regions of <pixels> may be below the threshold.
    01188     //In the historic version, this could affect new spot placement. This feature is not supported
    01189     //in this version.
    01190     double threshold = GV3::get<double>("threshold", 0, -1);
    01191     const double post_threshold = GV3::get<double>("post_dilate.threshold", -1, 1);
    01192     if(post_threshold != -1)
    01193         threshold = post_threshold;
    01194
    01195
    01196     //If dilation after region finding is to be performed, then do it here.
    01197     const double post_dilate_radius = GV3::get<double>("post_dilate.radius", 0, -1);
    01198     if(post_dilate_radius != 0)
    01199     {
    01200         Image<byte> pix(ims[0].size());
    01201         pix.fill(0);
    01202
    01203         for(unsigned int i=0; i < pixels.size(); i++)
    01204             pix[pixels[i]] = 255;
    01205
    01206         Image<byte> dilated = morphology(pix, getDisc(post_dilate_radius), Morphology::BinaryDilate<byte>());
    01207
    01208         pixels.clear();
    01209
    01210         ImageRef p(0,0);
    01211         do
    01212             if(dilated[p])
    01213                 pixels.push_back(p);
    01214         while(p.next(dilated.size()));
    01215     }
    01216
    01217
    01218     assert_same_size(ims);
    01219     if(log_ratios.size() != ims[0].size())
    01220     {
    01221         cerr << "Bad log ratios size\n";
    01222         exit(1);
    01223     }
    01224
    01225     vector<Vector<4> > spots;
    01226     //Spots can be either put down automatically, or specified
    01227     //The auto-initialization is very strange.
    01228     if(GV3::get<bool>("spots.auto_initialise", 1, 1))
    01229     {
    01230         //You never get two spots in the same disc in the second stage of the algorithm
    01231         vector<ImageRef> disc = getDisc(GV3::get<double>("spot_spread", 3.1, 1));
    01232
    01233         //Record all the pixels
    01234         map<ImageRef, double> valid_pixels;
    01235         for(unsigned int i=0; i < pixels.size(); i++)
    01236             if(log_ratios[pixels[i]] > threshold)
    01237                 valid_pixels.insert(make_pair(pixels[i], log_ratios[pixels[i]]));
    01238
    01239
    01240         //Get some initial spots by finding the local maxima
    01241         ImageRef neighbours[8] = {
    01242             ImageRef(-1, -1),
    01243             ImageRef( 0, -1),
    01244             ImageRef( 1, -1),
    01245
    01246             ImageRef(-1,  0),
    01247             ImageRef( 1,  0),
    01248
    01249             ImageRef(-1,  1),
    01250             ImageRef( 0,  1),
    01251             ImageRef( 1,  1),
    01252         };
    01253         for(unsigned int i=0; i < pixels.size(); i++)
    01254         {
    01255             if(!(log_ratios[pixels[i]] > threshold))
    01256                 goto not_a_maximum;
    01257
    01258             for(int j=0; j < 8; j++)
    01259                 if(!log_ratios.in_image(pixels[i] + neighbours[j]) || ! (log_ratios[pixels[i]] > log_ratios[pixels[i] + neighbours[j]]))
    01260                     goto not_a_maximum;
    01261
    01262             spots.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), log_normal_mode(blur_mu, blur_sigma), pixels[i].x, pixels[i].y));
    01263
    01264             //Remove the pixels around the initial spots
    01265             for(unsigned int j=0; j < disc.size(); j++)
    01266                 valid_pixels.erase(pixels[i] + disc[j]);
    01267
    01268             not_a_maximum:;
    01269         }
    01270
    01271         for(unsigned int i=0; i < spots.size(); i++)
    01272             cout << spots[i] << endl;
    01273
    01274         //Now place down extra spots in the remaining space.
    01275         while(!valid_pixels.empty())
    01276         {
    01277             ImageRef p = max_element(valid_pixels.begin(), valid_pixels.end(), LessSecond())->first;
    01278             spots.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), log_normal_mode(blur_mu, blur_sigma), p.x, p.y));
    01279
    01280             for(unsigned int j=0; j < disc.size(); j++)
    01281                 valid_pixels.erase(p + disc[j]);
    01282         }
    01283
    01284         //This line allows extra spots to be placed down around each spot already put down.
    01285         //This is a shocking hack and jenerally very unpleasant.
    01286         double extra_r = GV3::get<double>("extra_spots", 0, 1);
    01287         vector<ImageRef> extra = getDisc(extra_r);
    01288         vector<Vector<4> > more_spots;
    01289         for(unsigned int i=0; i < extra.size(); i++)
    01290             if(extra[i] != ImageRef_zero)
    01291                 for(unsigned int j=0; j < spots.size(); j++)
    01292                     more_spots.push_back(spots[j] + makeVector(0, 0, extra[i].x, extra[i].y) / (2*extra_r+1));
    01293
    01294         copy(more_spots.begin(), more_spots.end(), back_inserter(spots));
    01295     }
    01296     else
    01297     {
    01298         Vector<>  loaded_spots = GV3::get<Vector<> >("spots.manual_spots", "", -1);
    01299
    01300         if(loaded_spots.size()%4 != 0)
    01301         {
    01302             cerr << "Loaded spot size is not a multiple of 4\n";
    01303             exit(1);
    01304         }
    01305
    01306         else
    01307             spots = spots_to_vector(loaded_spots);
    01308     }
    01309
    01310     //Initialize the MT19937 RNG from a seed.
    01311     shared_ptr<MT19937> rng(new MT19937);
    01312     rng->simple_seed(GV3::get<int>("seed", 0, 1));
    01313
    01314     //Load in a checkpoint (precise RNG state, iteration and pass).
    01315     int start_iteration=0;
    01316     int start_pass=0;
    01317     if(GV3::get<bool>("checkpoint", 0, 1))
    01318     {
    01319         string rng_state = GV3::get<string>("checkpoint.rng.state", "", -1);
    01320         istringstream rs(rng_state);
    01321         rng->read(rs);
    01322         start_iteration=GV3::get<int>("checkpoint.iteration", 0, -1);
    01323         start_pass=GV3::get<int>("checkpoint.pass", 0, -1);
    01324     }
    01325
    01326     StateParameters p;
    01327     p.spots = spots;
    01328     p.rng = rng;
    01329     p.pass = start_pass;
    01330     p.iteration = start_iteration;
    01331     p.pixels = pixels;
    01332
    01333     return p;
    01334 }
    01335
    01336 ///Very simple and inefficient dilation function.
    01337 ///@ingroup gUtility
    01338 set<ImageRef>  dilate_mask(const vector<ImageRef>& v, double r)
    01339 {
    01340     vector<ImageRef> m = getDisc(r);
    01341
    01342     set<ImageRef> ret;
    01343
    01344     for(unsigned int i=0; i < v.size(); i++)
    01345         for(unsigned int j=0; j < m.size(); j++)
    01346             ret.insert(v[i] + m[j]);
    01347
    01348     return ret;
    01349 }
    01350
    01351 ///How far should steps in brightness be limited to?
    01352 ///@ingroup gStorm
    01353 double brightness_motion_limit(double mu, double sigma, bool not_one)
    01354 {
    01355     if(not_one)
    01356         return log_normal_std(mu, sigma);
    01357     else
    01358         return 1;
    01359 }
    01360
    01361
    01362
    01363
    01364 ///Mega class which actually does the meat of the spot fitting.
    01365 ///probably could be refactored a bit.
    01366 ///@ingroup gStorm
    01367 class FitSpots
    01368 {
    01369     const vector<Image<float> >& ims; ///< Input data
    01370     FitSpotsGraphics& graphics;       ///< Graphics class.
    01371     UserInterfaceCallback& ui;        ///< Callbacks to provide user interface
    01372     const vector<ImageRef> pixels;    ///< Area in which to perform model fitting
    01373
    01374     //Spot positions
    01375     vector<Vector<4> > spots;         ///< State in terms of current spot positions
    01376
    01377     //Starting point
    01378     const int start_iteration;        ///< Starting iteration number (for restarting from checkpoint)
    01379     int start_pass;                   ///< Starting pass (for restarting from checkpoint)
    01380
    01381     MT19937& rng;                     ///< Random numbewr generator
    01382
    01383     const double variance;            ///< Variance of noise in data. Must be 1.
    01384     const double intensity_mu;        ///< Prior for spot intensity
    01385     const double intensity_sigma;     ///< Prior for spot intensity
    01386     const double blur_mu;             ///< Prior for spot shape
    01387     const double blur_sigma;          ///< Prior for spt shape
    01388
    01389     //pixels is dilated slightly, by a disk of size area_extra_radius.
    01390     //The result is stored in allowed_area. This is used to implement a uniform
    01391     //prior over position, which is uniform within allowed_area, and zero
    01392     //everywhere else.
    01393     //
    01394     //This is implemented as a set, because the mask may extend beyond the borders of the
    01395     //image slightly
    01396     const double area_extra_radius;   ///< Extra size beyone marked region in which spots are allowed to exist
    01397     set<ImageRef> allowed_area;       ///< Total allowed region, pixels dilated by area_extra_radius
    01398     const int use_position_prior;     ///< Should a proper prior over position be uesd? A clue: yes.
    01399     const double position_prior;      ///< Value for the posision prior, i.e. reciprocal of area
    01400
    01401     //General optimizing and sampling parameters
    01402     const double max_motion;          ///< Maximum motion on any axes for optimization. See ConjugateGradientOnly.
    01403     const int sample_iterations;      ///< Number of mixing samples to use in Gibbs sampling
    01404
    01405     //Task specific optimizing and sampling parameters:
    01406
    01407     //Main optimization loop
    01408     const int main_cg_max_iterations;  ///< Maximum iterations allowed for ConjugateGradientOnly in main optimization loop
    01409     const int main_samples;            ///< Number of samples to use in main loop
    01410     const int main_passes;             ///< Number of passes to perform per iteration in main loop
    01411     const int outer_loop_iterations;   ///< Total number of iterations to perform
    01412
    01413     //Spot selection loop
    01414     const int add_remove_tries;        ///< Number of attemts to add/remove a spot
    01415     const int add_remove_opt_samples;  ///< Number of samples to use in model modification phase
    01416     const int add_remove_opt_retries;  ///< Number of attempts restarting the optimization to escape saddle points
    01417     const int add_remove_opt_hess_inner_samples; ///< Number of extra FFBS samples to use for computing Hessian when testing for convergence to non-saddle point
    01418     const int h_outer_samples;         ///< Number of samples used for computing Hessian as part of Laplace's approximation
    01419     const int h_inner_samples;         ///< Number of additional FFBS samples to use for computing Hessian as part of Laplace's approximation
    01420     const int tsamples;                ///< Number of samples to use in thermodynamic integration
    01421
    01422     // Average of all inputs (used for debugging)
    01423     const Image<float> ave;            ///< Average of input data: used for
    01424
    01425
    01426     //File to save output to
    01427     ofstream& save_spots;               ///< Output stream for log file
    01428
    01429     //Time accumulators for benchmarking
    01430     double time_gibbs;                  ///< Benchmarking data
    01431     double time_cg;                     ///< Benchmarking data
    01432
    01433     ///Motion limit for ConjugateGradientOnly
    01434     ///The x distance, y distance and spot size are all approximately the same scale
    01435     ///which is of order 1. The brigntness is not. By default, the brightness limit is also 1.
    01436     ///This flag controls whether is should be set to the standard deviation of the brightness prior
    01437     ///distribution. This setting will put the motion limit on the same scale as the other three
    01438     ///parameters.
    01439     const bool scale_brightness_limit;
    01440     const Vector<4> limit;              ///< Limit vector for ConjugateGradientOnly
    01441
    01442     const Matrix<3> A;                  ///< Transition matrix
    01443     const Vector<3> pi;                 ///< Initial probabilities
    01444
    01445
    01446     vector<vector<double> > pixel_intensities; ///<Pixel intensities for all images [frame][pixel]
    01447
    01448     DataForMCMC data_for_t_mcmc;  ///<Aggergated data for thermodynamic integration
    01449     DataForMCMC data_for_h_mcmc;  ///<Aggergated data for finding hessian
    01450
    01451     int iteration;                     ///< Iteration number
    01452
    01453
    01454     TIME(cvd_timer timer;)
    01455     public:
    01456
    01457     FitSpots(const vector<Image<float> >& ims_, FitSpotsGraphics& graphics_, UserInterfaceCallback& ui_, StateParameters& params, ofstream& save_spots_)
    01458     :ims(ims_), graphics(graphics_),ui(ui_),
    01459
    01460      //Set up the main parameters
    01461      pixels(params.pixels),
    01462      spots(params.spots),
    01463      start_iteration(params.iteration),
    01464      start_pass(params.pass),
    01465      rng(*params.rng),
    01466
    01467
    01468
    01469      //Set up all the system parameters
    01470      variance(1), // it should be
    01471
    01472      //To scale the X axis of a log-normal distribution, only
    01473      //the mu parameter needs to be changed...
    01474      intensity_mu(GV3::get<double>("intensity.rel_mu", 0., -1)  + log(sqrt(variance))),
    01475      intensity_sigma(GV3::get<double>("intensity.rel_sigma", 0., -1)),
    01476      blur_mu(GV3::get<double>("blur.mu", 0., -1)),
    01477      blur_sigma(GV3::get<double>("blur.sigma", 0., -1)),
    01478
    01479      //Spot position prior
    01480      area_extra_radius(GV3::get<double>("position.extra_radius", 0., -1)),
    01481      allowed_area(dilate_mask(pixels, area_extra_radius)),
    01482      use_position_prior(GV3::get<bool>("position.use_prior", true, -1)),
    01483      position_prior(1.0 / allowed_area.size()),
    01484
    01485      //General optimizing and sampling parameters
    01486      max_motion(GV3::get<double>("cg.max_motion", 0., -1)),
    01487      sample_iterations(GV3::get<int>("gibbs.mixing_iterations", 0, -1)),
    01488
    01489      //Task specific optimizing and sampling parameters:
    01490
    01491      //Main optimization loop
    01492      main_cg_max_iterations(GV3::get<double>("main.cg.max_iterations", 0., -1)),
    01493      main_samples(GV3::get<int>("main.gibbs.samples", 0, -1)),
    01494      main_passes(GV3::get<int>("main.passes", 0, -1)),
    01495      outer_loop_iterations(GV3::get<int>("main.total_iterations", 100000000, 1)),
    01496
    01497      //Spot selection loop
    01498      add_remove_tries(GV3::get<int>("add_remove.tries", 0, -1)),
    01499      add_remove_opt_samples(GV3::get<int>("add_remove.optimizer.samples", 0, -1)),
    01500      add_remove_opt_retries(GV3::get<int>("add_remove.optimizer.attempts", 0, -1)),
    01501      add_remove_opt_hess_inner_samples(GV3::get<int>("add_remove.optimizer.hessian_inner_samples", 0, -1)),
    01502      h_outer_samples(GV3::get<int>("add_remove.hessian.outer_samples", 0, -1)),
    01503      h_inner_samples(GV3::get<int>("add_remove.hessian.inner_samples", 0, -1)),
    01504      tsamples(GV3::get<int>("add_remove.thermo.samples", 0, -1)),
    01505
    01506      ave(average_image(ims_)),
    01507
    01508      save_spots(save_spots_),
    01509
    01510      time_gibbs(0),
    01511      time_cg(0),
    01512
    01513      scale_brightness_limit(GV3::get<bool>("max_motion.use_brightness_std", 0, -1)),
    01514      limit(makeVector(brightness_motion_limit(intensity_mu, intensity_sigma, scale_brightness_limit), 1, 1, 1)*max_motion),
    01515
    01516      A(GV3::get<Matrix<3> >("A", Zeros, 1)),
    01517      pi(GV3::get<Vector<3> >("pi", Zeros, 1)),
    01518
    01519
    01520      data_for_t_mcmc(pixels, pixel_intensities, intensity_mu, intensity_sigma, blur_mu, blur_sigma, variance, tsamples, sample_iterations, A, pi, rng),
    01521      data_for_h_mcmc(pixels, pixel_intensities, intensity_mu, intensity_sigma, blur_mu, blur_sigma, variance, h_outer_samples, sample_iterations, A, pi, rng)
    01522
    01523     {
    01524         assert_same_size(ims);
    01525
    01526         //Pixel intensities for all images [frame][pixel]
    01527         pixel_intensities.resize(ims.size(), vector<double>(pixels.size()));
    01528         for(unsigned int frame=0; frame < ims.size(); frame++)
    01529             for(unsigned int p=0; p < pixels.size(); p++)
    01530                 pixel_intensities[frame][p] = ims[frame][pixels[p]];
    01531
    01532     }
    01533
    01534     ///Perform a complete iteration of the optimzation stage of the spot firrint algorithm
    01535     void optimize_each_spot_in_turn_for_several_passes()
    01536     {
    01537         //Precompute the intensities for all spot pixels
    01538         vector<vector<double> > spot_intensities; //[spot][pixel]
    01539         for(unsigned int i=0; i < spots.size(); i++)
    01540             spot_intensities.push_back(compute_spot_intensity(pixels, spots[i]));
    01541
    01542         //Which pixels does each spot have?
    01543         vector<vector<int> > spot_pixels; //[spot][pixel]
    01544         spot_pixels.resize(spots.size());
    01545         for(unsigned int s=0; s < spots.size(); s++)
    01546             get_spot_pixels(pixels, spots[s], spot_pixels[s]);
    01547
    01548
    01549         //Optimize the model, N spots at a time.
    01550         //
    01551         for(int pass=start_pass; pass < main_passes; pass++)
    01552         {
    01553             save_spots << "Pass: " << pass << endl;
    01554             rng.write(save_spots);
    01555             save_spots << endl;
    01556
    01557             start_pass=0; // This is only nonzero first time, since we might chekpoint midway through an iteration
    01558             save_spots << "PASS" << pass << ": " << setprecision(20) << scientific << spots_to_Vector(spots) << endl;
    01559             save_spots << "ENDCHECKPOINT" << endl << flush;
    01560
    01561             ui.per_pass(iteration, pass, spots);
    01562             //Sort spots according to pass%4
    01563
    01564             //Sort the spots so that the optimization runs in sweeps
    01565             //This heiristic seems to increase the rate of propagation of information
    01566             //about spot positions.
    01567
    01568             //Create a list of indices
    01569             vector<int> index = sequence(spots.size());
    01570
    01571             int passs = pass + iteration;
    01572             //Sort the indices according to the position of the spot that they index
    01573             if(passs%4 == 0)
    01574                 sort(index.begin(), index.end(), IndexLexicographicPosition<less<double>, 2>(spots));
    01575             else if(passs%4==1)
    01576                 sort(index.begin(), index.end(), IndexLexicographicPosition<greater<double>, 2>(spots));
    01577             else if(passs%4==1)
    01578                 sort(index.begin(), index.end(), IndexLexicographicPosition<less<double>, 3>(spots));
    01579             else
    01580                 sort(index.begin(), index.end(), IndexLexicographicPosition<greater<double>, 3>(spots));
    01581
    01582             //Reorder the spots and their intensities and their pixel lists
    01583             {
    01584                 vector<Vector<4> > tmp_spot(index.size());
    01585                 vector<vector<double> > tmp_spot_intensities(index.size());
    01586                 vector<vector<int> > tmp_spot_pixels(index.size());
    01587                 for(unsigned int i=0; i < index.size(); i++)
    01588                 {
    01589                     tmp_spot[i] = spots[index[i]];
    01590                     swap(tmp_spot_intensities[i], spot_intensities[index[i]]);
    01591                     swap(tmp_spot_pixels[i], spot_pixels[i]);
    01592                 }
    01593
    01594                 swap(tmp_spot, spots);
    01595                 swap(tmp_spot_intensities, spot_intensities);
    01596                 swap(tmp_spot_pixels, spot_pixels);
    01597             }
    01598
    01599             //Sweep through and optimize each spot in turn
    01600             for(int s=0; s < (int)spots.size(); s++)
    01601             {
    01602                 ui.per_spot(iteration, pass, s, spots.size());
    01603                 ui.perhaps_stop();
    01604
    01605                 TIME(timer.reset();)
    01606                 //Get some samples with Gibbs sampling
    01607                 vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
    01608                 vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
    01609
    01610                 GibbsSampler2 sampler(pixel_intensities, spot_intensities, spots, spot_pixels, A, pi, variance, sample_iterations);
    01611                 for(int i=0; i < main_samples; i++)
    01612                 {
    01613                     sampler.next(rng);
    01614                     sample_list.push_back(sampler.sample());
    01615                     sample_intensities.push_back(sampler.sample_intensities());
    01616
    01617                     ui.perhaps_stop();
    01618                 }
    01619
    01620                 //First, remove the spot from all the samples.
    01621                 for(unsigned int i=0; i < sample_list.size(); i++)
    01622                     remove_spot(sample_intensities[i], spot_intensities[s], sample_list[i][s]);
    01623
    01624                 //cout << timer.get_time() << endl;
    01625                 TIME(time_gibbs += timer.reset();)
    01626
    01627                 //Package up all the data
    01628                 SampledBackgroundData data(sample_intensities, pixel_intensities, pixels,
    01629                                            intensity_mu, intensity_sigma, blur_mu, blur_sigma,
    01630                                            A, pi, variance);
    01631
    01632                 //Derivative computer:
    01633                 SpotNegProbabilityDiffWithSampledBackground compute_deriv(data);
    01634
    01635
    01636                 graphics.draw_pixels(pixels, 0, 0, 1, 1);
    01637                 graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), s);
    01638                 graphics.swap();
    01639
    01640                 //Optimize spot "s"
    01641                 //Optimize with the derivatives only since the actual probability
    01642                 //is much harder to compute
    01643                 ConjugateGradientOnly<4> cg(spots[s], compute_deriv, limit);
    01644
    01645
    01646                 cg.max_iterations = main_cg_max_iterations;
    01647
    01648
    01649                 #if 0
    01650                     cout << setprecision(10);
    01651                     cout << spots_to_Vector(spots) << endl;
    01652                     Matrix<4> hess, hess_errors;
    01653                     cout << "Hello, my name is Inigo Montoya\n";
    01654                     /*for(int i=0; i < 4; i++)
    01655                     {
    01656                         Matrix<4, 2> d = numerical_gradient_with_errors(NthDeriv(compute_deriv, i), cg.x);
    01657                         hess[i] = d.T()[0];
    01658                         hess_errors[i] = d.T()[1];
    01659                     }
    01660                     */
    01661                     //cout << "Errors:\n" << hess_errors << endl;
    01662                     //cout << "NHess:\n" << hess<< endl;
    01663                     Matrix<4> rhess =  -sampled_background_spot_hessian(cg.x, data);
    01664                     cout << "Hess:\n" << rhess << endl;
    01665                     //cout << "Err:\n" << hess - rhess << endl;
    01666
    01667                     //Vector<4> grad = compute_deriv(cg.x);
    01668
    01669                     //Matrix<4> e = hess - rhess;
    01670
    01671                     //for(int i=0; i < 4; i++)
    01672                     //  for(int j=0; j < 4; j++)
    01673                     //      e[i][j] /= hess_errors[i][j];
    01674
    01675                     //cout << "Err:\n" << e << endl;
    01676                     cout << "Deriv:" <<  compute_deriv(cg.x) << endl;
    01677                     //cout << "Full:\n" << sampled_background_spot_hessian2(cg.x, data) - grad.as_col()*grad.as_row() << endl;
    01678
    01679                     FreeEnergyHessian hesscomputer(data_for_h_mcmc);
    01680
    01681                     Matrix<4> nhess = hesscomputer.hessian(spots, 0);
    01682                     cout << "NHess:\n" << nhess << endl;
    01683
    01684                     cout << "Turbo-N Hess:\n" << sampled_background_spot_hessian_ffbs(cg.x, data, 10000) << endl;
    01685
    01686                     cout << "TI energy: " << NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(spots)) << endl;
    01687                     cout << "FA energy: " <<  SpotProbabilityWithSampledBackgroundFAKE(data)(cg.x) << endl;
    01688
    01689
    01690
    01691                     //cout << "Numerical hessian from FD:\n" << numerical_hessian(SpotProbabilityWithSampledBackgroundFAKE(data), cg.x) << endl;
    01692                     exit(0);
    01693                 #endif
    01694                 //cout << "Starting opt... " << cg.x << endl;
    01695                 while(cg.iterate(compute_deriv))
    01696                 {
    01697                     graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), s, cg.x);
    01698                     graphics.draw_pixels(pixels, 0, 0, 1, .2);
    01699                     graphics.swap();
    01700                     ui.perhaps_stop();
    01701                     //cout << cg.x << endl;
    01702                 }
    01703
    01704                 //Update to use the result of the optimization
    01705                 spots[s] = cg.x;
    01706                 //cout << "End: " << cg.x << endl;
    01707
    01708                 graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), -1);
    01709                 graphics.swap();
    01710
    01711                 //Recompute the new spot intensity, since the spot has changed
    01712                 spot_intensities[s] = compute_spot_intensity(pixels, spots[s]);
    01713
    01714                 //Recompute which are the useful pixels
    01715                 get_spot_pixels(pixels, spots[s], spot_pixels[s]);
    01716
    01717                 //Is the spot within the allowed area, i.e. is it's prior 0?
    01718                 //The prior is sero only if it we are using it and we're in an invalid area
    01719
    01720                 ImageRef quantized_spot_position = ir_rounded(spots[s].slice<2,2>());
    01721                 bool zero_prior = use_position_prior && (allowed_area.count(quantized_spot_position)==0);
    01722
    01723                 //Check to see if spot has been ejected. If spot_pixels is empty, then it has certainly been ejected.
    01724                 if(spot_pixels[s].empty() || zero_prior)
    01725                 {
    01726                     //Spot has been ejected. Erase it.
    01727                     cout  << " Erasing ejected spot: " << spot_pixels[s].empty() << " " << zero_prior << endl;
    01728                     cout << spots[s] << endl;
    01729                     //GUI_Pause(1);
    01730
    01731                     spot_intensities.erase(spot_intensities.begin() + s);
    01732                     spot_pixels.erase(spot_pixels.begin() + s);
    01733                     spots.erase(spots.begin() + s);
    01734                     s--;
    01735                     //exit(0);
    01736                 }
    01737
    01738                 //cout << timer.get_time() << endl;
    01739                 TIME(time_cg += timer.reset();)
    01740
    01741                 //cout << "Times: " << time_gibbs << " " << time_cg << endl;
    01742                 //save_spots << "INTERMEDIATE: " << setprecision(15) << scientific << spots_to_Vector(spots) << endl;
    01743             }
    01744         }
    01745     }
    01746
    01747     ///Perform a complete iteration of the model size modification stage of the spot fitting algorithm
    01748     void try_modifying_model()
    01749     {
    01750
    01751         for(int i=0; i < add_remove_tries; i++)
    01752         {
    01753             ui.per_modification(iteration, i, add_remove_tries);
    01754             ui.perhaps_stop();
    01755
    01756             cout << endl << endl << "Modifying the model" << endl << "======================\n";
    01757             cout << "Hello\n";
    01758             bool add_spot = (rng() > 0.5) || (spots.size() == 1);
    01759             cout << "World\n";
    01760
    01761             vector<Vector<4> > model_1, model_2;
    01762
    01763
    01764             int original; //What is the original model? Model 1 or Model 2?
    01765
    01766             if(add_spot)
    01767             {
    01768                 model_1 = spots;
    01769                 model_2 = model_1;
    01770
    01771                 //Pick a pixel within the thresholded ones as a starting point
    01772                 int r;
    01773                 do
    01774                 {
    01775                     r = (int)(rng() * pixels.size());
    01776                     //cout << r << " " << log_ratios[pixels[r]] << " " << pixels[r] << " " << threshold << endl;
    01777                 }
    01778                 while(0);
    01779
    01780                 //This version does not (yet?) suppotrt thresholding on log_ratios
    01781                 //for placing spots, since the purpose of log_ratios has diminished.
    01782                 //while(log_ratios[pixels[r]] < threshold);
    01783
    01784
    01785                 model_2.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma),
    01786                                              log_normal_mode(blur_mu, blur_sigma),
    01787                                              pixels[r].x + rng()-.5, pixels[r].y + rng() - .5));
    01788                 cout << "Adding a spot\n";
    01789
    01790                 original = 1;
    01791             }
    01792             else
    01793             {
    01794                 //Pick a random point to optimize/remove
    01795                 int a_random_spot = static_cast<int>(rng() * spots.size());
    01796                 model_1 = spots;
    01797                 swap(model_1[model_1.size()-1], model_1[a_random_spot]);
    01798
    01799                 model_2 = model_1;
    01800
    01801                 model_1.pop_back();
    01802                 cout << "Removing a spot\n";
    01803                 original = 2;
    01804             }
    01805
    01806             //The mobile spot is always the last spot of model_2
    01807             const int spot = model_2.size() - 1;
    01808
    01809             cout << "Original model: " << original << endl;
    01810
    01811             //Precompute the intensities for all spot pixels
    01812             //Model 2 is always one longer than model 1 and
    01813             //differs only on the extra element
    01814             vector<vector<double> > model2_spot_intensities; //[spot][pixel]
    01815             for(unsigned int i=0; i < model_2.size(); i++)
    01816                 model2_spot_intensities.push_back(compute_spot_intensity(pixels, model_2[i]));
    01817
    01818             //Which pixels does each spot have?
    01819             vector<vector<int> > model2_spot_pixels(model_2.size()); //[spot][pixel]
    01820             for(unsigned int s=0; s < model_2.size(); s++)
    01821                 get_spot_pixels(pixels, model_2[s], model2_spot_pixels[s]);
    01822
    01823             //Optimize spot:
    01824             {
    01825                 cout << "Optimizing spot for model selection\n";
    01826
    01827
    01828                 //Get some samples with Gibbs sampling
    01829                 vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
    01830                 vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
    01831
    01832                 GibbsSampler2 sampler(pixel_intensities, model2_spot_intensities, model_2, model2_spot_pixels, A, pi, variance, sample_iterations);
    01833                 for(int i=0; i < add_remove_opt_samples; i++)
    01834                 {
    01835                     sampler.next(rng);
    01836                     sample_list.push_back(sampler.sample());
    01837                     sample_intensities.push_back(sampler.sample_intensities());
    01838                     ui.perhaps_stop();
    01839                 }
    01840
    01841                 //First, remove the spot from all the samples.
    01842                 for(unsigned int i=0; i < sample_list.size(); i++)
    01843                     remove_spot(sample_intensities[i], model2_spot_intensities[spot], sample_list[i][spot]);
    01844
    01845                 //Package up all the data
    01846                 SampledBackgroundData data(sample_intensities, pixel_intensities, pixels,
    01847                                            intensity_mu, intensity_sigma, blur_mu, blur_sigma,
    01848                                            A, pi, variance);
    01849
    01850                 //Derivative computer:
    01851                 SpotNegProbabilityDiffWithSampledBackground compute_deriv(data);
    01852
    01853                 graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), spot);
    01854                 graphics.swap();
    01855
    01856                 //Optimize spot "s"
    01857                 //Optimize with the derivatives only since the actual probability
    01858                 //is much harder to compute
    01859                 ConjugateGradientOnly<4> cg(model_2[spot], compute_deriv, limit);
    01860
    01861
    01862                 for(int attempt=0; attempt < add_remove_opt_retries; attempt++)
    01863                 {
    01864                     cout << "Attempt " << attempt << " of " << add_remove_opt_retries << endl;
    01865                     ui.perhaps_stop();
    01866
    01867                     //Optimize with conjugate gradient
    01868                     while(cg.iterate(compute_deriv))
    01869                     {
    01870                         ui.perhaps_stop();
    01871                         graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), spot, cg.x);
    01872                         graphics.swap();
    01873                     }
    01874
    01875                     //Check for being at a saddle point (no point checking on the last try)
    01876                     //All eigenvectors should be negative at a maximum.
    01877                     //WTF: is this a bug? WTF?
    01878                     //It was this:
    01879                     //if(attempt < add_remove_tries - 1)
    01880                     if(attempt < add_remove_opt_retries - 1)
    01881                     {
    01882                         Matrix<4> hessian = sampled_background_spot_hessian_ffbs(cg.x, data, add_remove_opt_hess_inner_samples, rng);
    01883                         SymEigen<4> hess_decomp(hessian);
    01884
    01885                         //cout << "What the fuck:" << endl << hessian << endl << hessian3<< endl << hessian2 << endl;
    01886
    01887                         cout << "Eigenvalues are: " << hess_decomp.get_evalues() << endl;
    01888
    01889                         if(hess_decomp.is_negdef())
    01890                             break;
    01891                         else
    01892                         {
    01893                             //Restart in the direction of the best uphill part
    01894                             cg.init(cg.x + 0.1 * hess_decomp.get_evectors()[3], (hess_decomp.get_evectors()[3]));
    01895
    01896
    01897                             cout << "Grad = " << compute_deriv(cg.x) << endl;
    01898                             for(int i=0; i < 4; i++)
    01899                             {
    01900                                 cout << "Direction: " << i << endl;
    01901                                 cout << unit(compute_deriv(cg.x + 0.1*hess_decomp.get_evectors()[i])) * hess_decomp.get_evectors()[i] << endl;
    01902                             }
    01903
    01904                             for(int i=0; i < 4; i++)
    01905                             {
    01906                                 cout << "Direction: " << i << endl;
    01907                                 Vector<4> d = Zeros;
    01908                                 d[i] = 1;
    01909                                 cout << unit(compute_deriv(cg.x + d)) * hess_decomp.get_evectors()[i] << endl;
    01910                             }
    01911                         }
    01912                     }
    01913                 }
    01914
    01915
    01916                 //Update to use the result of the optimization
    01917                 model_2[spot] = cg.x;
    01918
    01919                 graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), -1);
    01920                 graphics.swap();
    01921
    01922
    01923                 //Update cached data based on the new spot position
    01924                 model2_spot_intensities[spot] = compute_spot_intensity(pixels, model_2[spot]);
    01925                 get_spot_pixels(pixels, model_2[spot], model2_spot_pixels[spot]);
    01926
    01927                 cout << "Done optimizing for model selection\n";
    01928             }
    01929
    01930
    01931             //Which model to keep?
    01932             int keep=original;
    01933
    01934             //Compute position prior (and we might be able to reject it really quickly here)
    01935             bool zero_prior = use_position_prior && (allowed_area.count(ir_rounded(model_2[spot].slice<2,2>()))==0);
    01936
    01937             if(zero_prior)
    01938             {
    01939                 //Model 2 went bad, so we clearly keep model 1
    01940                 keep = 1;
    01941             }
    01942             else
    01943             {
    01944
    01945                 //The position prior is independent
    01946                 //Compute the difference  model2 - model1
    01947                 //This is only valid if model2 is in the valid region
    01948                 double position_log_prior_model2_minus_model1;
    01949                 if(use_position_prior)
    01950                     position_log_prior_model2_minus_model1 = (model_2.size() - model_1.size()) * ln(position_prior);
    01951                 else
    01952                     position_log_prior_model2_minus_model1 = 0;
    01953
    01954
    01955                 //Integrate model_2
    01956                 //First compute the Hessian since this might go wrong.
    01957
    01958                 //FreeEnergyHessian hesscomputer(data_for_h_mcmc);
    01959                 Matrix<4> hess;// = hesscomputer.hessian(model_2, spot);
    01960
    01961                 //Use turbohess here since it is much faster, as the backwards sampling step is fast
    01962                 //We expect this hessian to be really quite different, if the spot has moved from
    01963                 //a long way away, since the sampling will have changed dramatically
    01964                 {
    01965                     //Get some samples with Gibbs sampling
    01966                     vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling
    01967                     vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel]
    01968
    01969                     GibbsSampler sampler(pixel_intensities, model2_spot_intensities, model_2, A, pi, variance, sample_iterations);
    01970                     for(int i=0; i < h_outer_samples; i++)
    01971                     {
    01972                         ui.perhaps_stop();
    01973                         sampler.next(rng);
    01974                         sample_list.push_back(sampler.sample());
    01975                         sample_intensities.push_back(sampler.sample_intensities());
    01976                     }
    01977
    01978                     //First, remove the spot from all the samples.
    01979                     for(unsigned int i=0; i < sample_list.size(); i++)
    01980                         remove_spot(sample_intensities[i], model2_spot_intensities[spot], sample_list[i][spot]);
    01981
    01982                     //Package up all the data
    01983                     SampledBackgroundData data(sample_intensities, pixel_intensities, pixels,
    01984                                                intensity_mu, intensity_sigma, blur_mu, blur_sigma,
    01985                                                A, pi, variance);
    01986
    01987                     hess = sampled_background_spot_hessian_ffbs(model_2[spot], data, h_inner_samples, rng);
    01988                 }
    01989
    01990
    01991                 double det = determinant(-hess / (M_PI*2));
    01992                 SymEigen<4> hess_decomp(-hess);
    01993                 cout << "Hessien Eigenvalues are: " << hess_decomp.get_evalues() << endl;
    01994                 const double smallest_evalue = 1e-6;
    01995
    01996                 //If the determinant is negative, then we are still at a saddle point
    01997                 //despite the attempts above, so abandon the attempt.
    01998                 if(hess_decomp.get_evalues()[0] >  smallest_evalue)
    01999                 {
    02000
    02001                     //Compute the free energy of model 2 at the MLE estimate
    02002                     cout << "Model 2:\n";
    02003             //      double model_2_energy = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_2));
    02004                     double model_2_energy = -NegativeFreeEnergy(data_for_t_mcmc).compute_with_mask(spots_to_Vector(model_2), model2_spot_pixels);
    02005                     cout << "Energy: " << model_2_energy << endl;
    02006
    02007                     //Combine the MLE energy and Hessian using Laplace's approximation
    02008                     double model_2_occam =  -0.5 * log(det);
    02009                     double model_2_prob = model_2_energy + model_2_occam + position_log_prior_model2_minus_model1;
    02010
    02011                     cout << "Occam: " << model_2_occam << endl;
    02012                     cout << "Position: " << position_log_prior_model2_minus_model1 << endl;
    02013                     cout << "Prob: " << model_2_prob << endl;
    02014
    02015                     //Integrate model_1
    02016                     //It has no parameters, in this formulation.
    02017                     //double model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_1));
    02018
    02019                     //Note that model_1 always has one fewer spots, and the last spot is always the one
    02020                     //missing, so we can make the correct mask very easily:
    02021                     model2_spot_pixels.pop_back();
    02022                     double model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc).compute_with_mask(spots_to_Vector(model_1), model2_spot_pixels);
    02023                     cout << "Prob: " << model_1_prob << endl;
    02024                     //model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_1));
    02025
    02026                     if(model_2_prob  >  model_1_prob)
    02027                         keep=2;
    02028                     else
    02029                         keep=1;
    02030
    02031                     cout << "Models evaluated\n";
    02032                 }
    02033                 else
    02034                 {
    02035                     cout << "Determinant has bad eigenvalues!\n";
    02036                     keep = original;
    02037                     cout << hess_decomp.get_evalues() << endl;
    02038                 }
    02039             }
    02040
    02041             if(keep == 2)
    02042             {
    02043                 spots = model_2;
    02044                 cout << "Keeping model 2\n";
    02045
    02046             }
    02047             else
    02048             {
    02049                 spots = model_1;
    02050                 cout << "Keeping model 1\n";
    02051             }
    02052
    02053             if(original != keep)
    02054             {
    02055                 cout << "Model changed!\n";
    02056                 //break;
    02057             }
    02058         }
    02059     }
    02060
    02061     ///Run the complete optimization algorithm.
    02062     void run()
    02063     {
    02064         graphics.init(ims[0].size());
    02065         save_spots << "LOGVERSION 2" << endl;
    02066         save_spots << "PIXELS";
    02067         for(unsigned int i=0; i < pixels.size(); i++)
    02068             save_spots << " " << pixels[i].x << " " << pixels[i].y;
    02069         save_spots << endl;
    02070
    02071         save_spots << "BEGINGVARLIST" << endl;
    02072         GV3::print_var_list(save_spots, "", 1);
    02073         save_spots << "ENDGVARLIST" << endl;
    02074
    02075         //TODO all GVARS are set, so dump out gvars.
    02076
    02077         cout << "Limit vector: " << limit << endl;
    02078
    02079         for(iteration=start_iteration; iteration < outer_loop_iterations ;iteration++)
    02080         {
    02081             save_spots << "Iteration: " << iteration << " (" << iteration *  main_passes << ")" << endl;
    02082             save_spots << "MAIN: " << setprecision(20) << scientific << spots_to_Vector(spots) << endl;
    02083
    02084             cout << endl << endl << "----------------------" << endl << "Optimizing:\n";
    02085             cout << spots.size() << endl;
    02086
    02087
    02088             optimize_each_spot_in_turn_for_several_passes();
    02089
    02090             //spot_intensities is be correct here!
    02091             try_modifying_model();
    02092         }
    02093         save_spots << "FINAL: " << setprecision(15) << scientific << spots_to_Vector(spots) << endl;
    02094     }
    02095
    02096 };
    02097
    02098 ///Wrapper function for using FitSpots
    02099 ///@ingroup gStorm
    02100 void fit_spots_new(const vector<Image<float> >& ims, StateParameters& p, ofstream& save_spots, FitSpotsGraphics& gr)
    02101 {
    02102     auto_ptr<UserInterfaceCallback> ui = null_ui();
    02103     FitSpots fit(ims, gr, *ui, p, save_spots);
    02104     fit.run();
    02105 }
    02106
    02107 ///Wrapper function for using FitSpots
    02108 ///@ingroup gStorm
    02109 void fit_spots_new(const vector<Image<float> >& ims, StateParameters& p, ofstream& save_spots, FitSpotsGraphics& gr, UserInterfaceCallback& ui)
    02110 {
    02111     try{
    02112         FitSpots fit(ims, gr, ui, p, save_spots);
    02113         fit.run();
    02114     }
    02115     catch(UserInterfaceCallback::UserIssuedStop)
    02116     {
    02117     }
    02118 }