ThreeB 1.1

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

- Class List
- Class Hierarchy
- Class Members

Public Member Functions | Protected Attributes
NegativeFreeEnergy Class Reference

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

Inheritance diagram for NegativeFreeEnergy:

List of all members.

## 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 280 of file multispot5.cc.

## Constructor & Destructor Documentation

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

Parameters:

  d | Necessary data

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

References 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 309 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 391 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 208 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 198 of file multispot5.cc.

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

Definition at line 199 of file multispot5.cc.

const double DataForMCMC::mu_brightness`[protected, inherited]`

Definition at line 200 of file multispot5.cc.

const double DataForMCMC::sigma_brightness`[protected, inherited]`

Definition at line 200 of file multispot5.cc.

const double DataForMCMC::mu_blur`[protected, inherited]`

Definition at line 200 of file multispot5.cc.

const double DataForMCMC::sigma_blur`[protected, inherited]`

Definition at line 200 of file multispot5.cc.

const double DataForMCMC::variance`[protected, inherited]`

Definition at line 201 of file multispot5.cc.

const int DataForMCMC::samples`[protected, inherited]`

Definition at line 202 of file multispot5.cc.

const int DataForMCMC::sample_iterations`[protected, inherited]`

Definition at line 202 of file multispot5.cc.

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

Definition at line 203 of file multispot5.cc.

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

Definition at line 204 of file multispot5.cc.

MT19937& DataForMCMC::rng`[protected, inherited]`

Definition at line 205 of file multispot5.cc.

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

- multispot5.cc