Class for computing the negitve free energy using thermodynamic integration. More...

Inheritance diagram for NegativeFreeEnergy:

Public Member Functions
	NegativeFreeEnergy (const DataForMCMC &d)
double	variance_from_sample (double sample, double samples, double base_sigma, double scale_pow) const
double	compute_with_mask (const Vector<> &spots, const vector< vector< int > > &spot_pixels) const
double	operator() (const Vector<> &spots) const
MT19937 &	get_rng () const
Protected Attributes
const vector< ImageRef > &	pixels
const vector< vector< double > > &	pixel_intensities
const double	mu_brightness
const double	sigma_brightness
const double	mu_blur
const double	sigma_blur
const double	variance
const int	samples
const int	sample_iterations
const Matrix< 3 >	A
const Vector< 3 >	pi
MT19937 &	rng

Detailed Description

Class for computing the negitve free energy using thermodynamic integration.

Definition at line 261 of file multispot5.cc.

Constructor & Destructor Documentation

NegativeFreeEnergy::NegativeFreeEnergy ( const DataForMCMC & d ) [inline]

Parameters:

d	Necessary data

Definition at line 266 of file multispot5.cc.

    :DataForMCMC(d)
    {
    }

Member Function Documentation

double NegativeFreeEnergy::variance_from_sample	(	double	sample,
		double	samples,
		double	base_sigma,
		double	scale_pow
	)		const `[inline]`

Give the noise variance given a sample number and growth parameters.

Parameters:

sample	Sample number
samples	Total number of samples
base_sigma	Starting value of sigme
scale_pow	Exponent scaling

Definition at line 276 of file multispot5.cc.

References scale(), and sq().

    {
        double scale = pow(1.25, sample * 1. / samples * scale_pow);
        double sigma = base_sigma * scale;
        double new_variance = sq(sigma);

        return new_variance;
    }

double NegativeFreeEnergy::compute_with_mask	(	const Vector<> &	spots,
		const vector< vector< int > > &	spot_pixels
	)		const `[inline]`

Estimate free energy using the Slow Growth Thermodynamic Integration method given in Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford.

M. Neal, 1993 Except using a 5 point stencil instead of forward differenceing in Eq 6.30

Parameters:

spots	list of spots
spot_pixels	Mask around each spot to use to save on computation

Definition at line 290 of file multispot5.cc.

References Kahan::add(), SampledMultispot::compute_spot_intensity(), DataForMCMC::get_rng(), log_log_normal(), SampledMultispot::GibbsSampler2::next(), SampledMultispot::GibbsSampler2::sample_intensities(), SampledMultispot::GibbsSampler2::set_variance(), spots_to_vector(), and sq().

Referenced by FitSpots::try_modifying_model().

    {
        //Estimate free energy using the Slow Growth Thermodynamic Integration method given in
        //Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford. M. Neal, 1993
        //Except using a 5 point stencil instead of forward differenceing in Eq 6.30
        double base_sigma = sqrt(variance); // should be 1
        double scale_pow = 100.0;

        const unsigned int nspots  = spots.size()/4;
        const unsigned int nframes = pixel_intensities.size();
        const unsigned int npixels = pixels.size();
        assert(spots.size() %4 == 0);
        assert(spot_pixels.size() == nspots);

        //Compute the intensities and derivatives for all spot
        vector<vector<double> > spot_intensity; //[spot][pixel]
        for(unsigned int i=0; i < nspots; i++)
            spot_intensity.push_back(compute_spot_intensity(pixels, spots.slice<Dynamic,4>(i*4,4)));
        
        GibbsSampler2 sampler(pixel_intensities, spot_intensity, spots_to_vector(spots), spot_pixels, A, pi, variance, sample_iterations);

        double sum = 0;
        Kahan ksum;
        for(int sample=0; sample < samples; sample++)
        {
            //Compute the positions of the surrounding steps
            double var1 = variance_from_sample(sample-2, samples, base_sigma, scale_pow);
            double var2 = variance_from_sample(sample-1, samples, base_sigma, scale_pow);
            double var3 = variance_from_sample(sample+1, samples, base_sigma, scale_pow);
            double var4 = variance_from_sample(sample+2, samples, base_sigma, scale_pow);

            //Take a sample
            sampler.set_variance(var2);
            sampler.next(DataForMCMC::get_rng());
            
            //Compute the SSD. This is fixed regardless of sigma.   
            double err_sum=0;
            for(unsigned int frame=0; frame < nframes; frame++)
                for(unsigned int pixel=0; pixel < npixels; pixel++)
                    err_sum -= sq(pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel]);
            
            //Compute the derivative using a five point stencil
            //This could be done in a better but less clear way
            double e1 = err_sum / (2*var1) - npixels*nframes*::log(2*M_PI*var1)/2;
            double e2 = err_sum / (2*var2) - npixels*nframes*::log(2*M_PI*var2)/2;
            double e3 = err_sum / (2*var3) - npixels*nframes*::log(2*M_PI*var3)/2;
            double e4 = err_sum / (2*var4) - npixels*nframes*::log(2*M_PI*var4)/2;
            sum += (-e1 + 8*e2 - 8*e3 + e4)/12;

            ksum.add(-e1/12);
            ksum.add(8*e2/12);
            ksum.add(-8*e3/12);
            ksum.add(e4/12);
        }

        double log_final = (log(variance_from_sample(samples, samples, base_sigma, scale_pow)*2*M_PI)/2) * npixels * nframes;
        
        double priors=0;
        for(unsigned int i=0; i < nspots; i++)
        {
            priors += log_log_normal(spots[i*4+0], mu_brightness, sigma_brightness);
            priors += log_log_normal(spots[i*4+1], mu_blur, sigma_blur);
        }

        /*cout << "Thermo:\n";
        cout << "sum: " << sum -log_final << endl;
        cout << "kah: " << ksum.sum -log_final << endl;
        cout << "priors: " << priors + sum -log_final << endl;
        cout << "   kah: " << priors + ksum.sum -log_final << endl;
*/
        //cout << log_final << endl;
        //cout << sum + log_final << endl;
        
        
        sampler.set_variance(variance);
        return -(sum+priors - log_final);
    }

double NegativeFreeEnergy::operator() ( const Vector<> & spots ) const [inline]

Estimate free energy using the Slow Growth Thermodynamic Integration method given in Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford.

M. Neal, 1993 Except using a 5 point stencil instead of forward differenceing in Eq 6.30

Parameters:

spots list of spots

Definition at line 372 of file multispot5.cc.

References Kahan::add(), SampledMultispot::compute_spot_intensity(), DataForMCMC::get_rng(), log_log_normal(), SampledMultispot::GibbsSampler::next(), SampledMultispot::GibbsSampler::sample_intensities(), SampledMultispot::GibbsSampler::set_variance(), spots_to_vector(), and sq().

    {
        double base_sigma = sqrt(variance); // should be 1
        double scale_pow = 100.0;

        const unsigned int nspots  = spots.size()/4;
        const unsigned int nframes = pixel_intensities.size();
        const unsigned int npixels = pixels.size();
        assert(spots.size() %4 == 0);

        //Compute the intensities and derivatives for all spot
        vector<vector<double> > spot_intensity; //[spot][pixel]
        for(unsigned int i=0; i < nspots; i++)
            spot_intensity.push_back(compute_spot_intensity(pixels, spots.slice<Dynamic,4>(i*4,4)));
        
        GibbsSampler sampler(pixel_intensities, spot_intensity, spots_to_vector(spots), A, pi, variance, sample_iterations);

        double sum = 0;
        Kahan ksum;
        for(int sample=0; sample < samples; sample++)
        {
            //Compute the positions of the surrounding steps
            double var1 = variance_from_sample(sample-2, samples, base_sigma, scale_pow);
            double var2 = variance_from_sample(sample-1, samples, base_sigma, scale_pow);
            double var3 = variance_from_sample(sample+1, samples, base_sigma, scale_pow);
            double var4 = variance_from_sample(sample+2, samples, base_sigma, scale_pow);

            //Take a sample
            sampler.set_variance(var2);
            sampler.next(DataForMCMC::get_rng());
            
            //Compute the SSD. This is fixed regardless of sigma.   
            double err_sum=0;
            for(unsigned int frame=0; frame < nframes; frame++)
                for(unsigned int pixel=0; pixel < npixels; pixel++)
                    err_sum -= sq(pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel]);
            
            //Compute the derivative using a five point stencil
            //This could be done in a better but less clear way
            double e1 = err_sum / (2*var1) - npixels*nframes*::log(2*M_PI*var1)/2;
            double e2 = err_sum / (2*var2) - npixels*nframes*::log(2*M_PI*var2)/2;
            double e3 = err_sum / (2*var3) - npixels*nframes*::log(2*M_PI*var3)/2;
            double e4 = err_sum / (2*var4) - npixels*nframes*::log(2*M_PI*var4)/2;
            sum += (-e1 + 8*e2 - 8*e3 + e4)/12;

            ksum.add(-e1/12);
            ksum.add(8*e2/12);
            ksum.add(-8*e3/12);
            ksum.add(e4/12);
        }

        double log_final = (log(variance_from_sample(samples, samples, base_sigma, scale_pow)*2*M_PI)/2) * npixels * nframes;
        
        double priors=0;
        for(unsigned int i=0; i < nspots; i++)
        {
            priors += log_log_normal(spots[i*4+0], mu_brightness, sigma_brightness);
            priors += log_log_normal(spots[i*4+1], mu_blur, sigma_blur);
        }

        /*cout << "Thermo:\n";
        cout << "sum: " << sum -log_final << endl;
        cout << "kah: " << ksum.sum -log_final << endl;
        cout << "priors: " << priors + sum -log_final << endl;
        cout << "   kah: " << priors + ksum.sum -log_final << endl;
*/
        //cout << log_final << endl;
        //cout << sum + log_final << endl;
        
        
        sampler.set_variance(variance);
        return -(sum+priors - log_final);
    }

MT19937& DataForMCMC::get_rng ( ) const [inline, inherited]

Definition at line 189 of file multispot5.cc.

Referenced by compute_with_mask(), and operator()().

    {
        return rng;
    }

Member Data Documentation

const vector<ImageRef>& DataForMCMC::pixels [protected, inherited]

Definition at line 179 of file multispot5.cc.

const vector<vector<double> >& DataForMCMC::pixel_intensities [protected, inherited]

Definition at line 180 of file multispot5.cc.

const double DataForMCMC::mu_brightness [protected, inherited]

Definition at line 181 of file multispot5.cc.

const double DataForMCMC::sigma_brightness [protected, inherited]

Definition at line 181 of file multispot5.cc.

const double DataForMCMC::mu_blur [protected, inherited]

Definition at line 181 of file multispot5.cc.

const double DataForMCMC::sigma_blur [protected, inherited]

Definition at line 181 of file multispot5.cc.

const double DataForMCMC::variance [protected, inherited]

Definition at line 182 of file multispot5.cc.

const int DataForMCMC::samples [protected, inherited]

Definition at line 183 of file multispot5.cc.

const int DataForMCMC::sample_iterations [protected, inherited]

Definition at line 183 of file multispot5.cc.

const Matrix<3> DataForMCMC::A [protected, inherited]

Definition at line 184 of file multispot5.cc.

const Vector<3> DataForMCMC::pi [protected, inherited]

Definition at line 185 of file multispot5.cc.

MT19937& DataForMCMC::rng [protected, inherited]

Definition at line 186 of file multispot5.cc.

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

multispot5.cc

Public Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation