alternaria_moves.cpp

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <config.h>

#include <cstdlib>
#include <cmath>
#include <cassert>
#include <iostream>
#include <iomanip>

#include <inttypes.h>

#include <jwsc/base/limits.h>
#include <jwsc/math/constants.h>
#include <jwsc/prob/pdf.h>
#include <jwsc/prob/pmf.h>

#include <jwsc++/base/exception.h>

#include "sampler.h"
#include "structure.h"
#include "hypha.h"
#include "spore.h"
#include "alternaria_model.h"
#include "alternaria_moves.h"


using std::isinf;
using std::isnan;
using namespace jwsc;
using jwscxx::base::Exception;
using jwscxx::base::Arg_error;


bool Apical_hypha_birth_death_move::run_1(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Structure* parent;
    Apical_hypha* hypha;
    const Apical_hypha_density* density;

    density = params->alternaria_proposal->get_apical_hypha_d();

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        parent = params->alternaria_proposal->get_random_terminal(level);

        try
        {
            hypha = density->sample(parent);
            params->alternaria_proposal->add_apical_hypha(hypha);
        }
        catch (Exception e)
        {
            return false;
        }
    }
    catch (Exception e)
    {
        try
        {
            hypha = density->sample(params->alternaria->get_theta(),
                                    params->alternaria->get_psi(),
                                    params->alternaria->get_base_level());
            params->alternaria_proposal->add_apical_hypha(hypha);
        }
        catch (Exception ee)
        {
            return false;
        }
    }

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
                params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double hypha_log_prob = hypha->get_log_prob();

    double birth_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double death_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha -= hypha_log_prob;
    log_alpha -= birth_move_log_prob;
    log_alpha += death_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


bool Apical_hypha_birth_death_move::run_2(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Apical_hypha* hypha;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        hypha = params->alternaria_proposal
                      ->get_random_terminal_apical_hypha(level);

        // Check that hypha doesn't have any laterals.
        if (hypha->get_laterals().size() != 0)
        {
            return false;
        }

        // TEMPORARY Check that this is not the root.
        if (!(hypha->get_parent()))
        {
            return false;
        }
    }
    catch (Exception e)
    {
        return false;
    }

    double hypha_log_prob = hypha->get_log_prob();

    try
    {
        params->alternaria_proposal->remove_apical_hypha(hypha);
    }
    catch (Arg_error e)
    {
        return false;
    }

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
                params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double birth_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double death_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += hypha_log_prob;
    log_alpha += birth_move_log_prob;
    log_alpha -= death_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}




bool Spore_birth_death_move::run_1(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Structure* parent;
    Spore* spore;
    const Spore_density* density;

    density = params->alternaria_proposal->get_spore_d();

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        if (level == 0)
        {
            return false;
        }

        parent = params->alternaria_proposal->get_random_terminal(level);

        try
        {
            spore = density->sample(parent);
            params->alternaria_proposal->add_spore(spore);
        }
        catch (Exception e)
        {
            return false;
        }
    }
    catch (Exception e)
    {
        try
        {
            spore = density->sample(params->alternaria->get_theta(),
                                    params->alternaria->get_psi(),
                                    params->alternaria->get_base_level());
            params->alternaria_proposal->add_spore(spore);
        }
        catch (Exception ee)
        {
            return false;
        }
    }

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
                params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double spore_log_prob = spore->get_log_prob();

    double birth_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double death_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha -= spore_log_prob;
    log_alpha -= birth_move_log_prob;
    log_alpha += death_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


bool Spore_birth_death_move::run_2(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Spore* spore;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        spore = params->alternaria_proposal->get_random_terminal_spore(level);

        // Check that spore doesn't have any laterals.
        if (spore->get_laterals().size() != 0)
        {
            return false;
        }

        // TEMPORARY Check that this is not the root.
        if (!(spore->get_parent()))
        {
            return false;
        }
    }
    catch (Exception e)
    {
        return false;
    }

    double spore_log_prob = spore->get_log_prob();

    try
    {
        params->alternaria_proposal->remove_spore(spore);
    }
    catch (Arg_error e)
    {
        return false;
    }

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
                params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double birth_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double death_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += spore_log_prob;
    log_alpha += birth_move_log_prob;
    log_alpha -= death_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}

bool Spore_transition_apical_move::run_1(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Apical_hypha* hypha;
    Spore* spore = 0;
    const Spore_density* density;
    double apical_to_spore_proposal_log_prob;
    double spore_to_apical_proposal_log_prob;

    density = params->alternaria_proposal->get_spore_d();

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        if (level == 0)
        {
            return false;
        }

        hypha = params->alternaria_proposal->get_random_apical_hypha(level);
        spore_to_apical_proposal_log_prob = hypha->get_log_prob();

        spore = new Spore(*hypha, density->sample_width(), density);
    }
    catch (Exception e)
    {
        return false;
    }

    apical_to_spore_proposal_log_prob = spore->get_log_prob();

    try
    {
        params->alternaria_proposal->replace_apical_hypha_with_spore(hypha,
                spore);
    }
    catch (Exception e)
    {
        // We are going to try and handle this now in the sampler by catching
        // the exception and rejecting the proposal.
        throw e;
    }

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
                params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double apical_to_spore_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) 
        ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double spore_to_apical_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) 
        ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha -= apical_to_spore_proposal_log_prob;
    log_alpha += spore_to_apical_proposal_log_prob;
    log_alpha -= apical_to_spore_move_log_prob;
    log_alpha += spore_to_apical_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


bool Spore_transition_apical_move::run_2(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Apical_hypha* hypha = 0;
    Spore* spore;
    const Apical_hypha_density* density;
    double spore_to_apical_proposal_log_prob;
    double apical_to_spore_proposal_log_prob;

    density = params->alternaria_proposal->get_apical_hypha_d();

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        spore = params->alternaria_proposal->get_random_spore(level);
        apical_to_spore_proposal_log_prob = spore->get_log_prob();

        hypha = new Apical_hypha(*spore, density->sample_width(), density);
    }
    catch (Exception e)
    {
        return false;
    }

    spore_to_apical_proposal_log_prob = hypha->get_log_prob();

    try
    {
        params->alternaria_proposal->replace_spore_with_apical_hypha(spore,
                hypha);
    }
    catch (Exception e)
    {
        // We are going to try and handle this now in the sampler by catching
        // the exception and rejecting the proposal.
        throw e;
    }

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
                params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double apical_to_spore_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) 
        ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double spore_to_apical_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) 
        ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += apical_to_spore_proposal_log_prob;
    log_alpha -= spore_to_apical_proposal_log_prob;
    log_alpha += apical_to_spore_move_log_prob;
    log_alpha -= spore_to_apical_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}

bool Spore_split_merge_2_apical_move::run_1(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Spore* spore;

    Apical_hypha* split_1 = 0;
    double split_1_log_prob;

    Apical_hypha* split_2 = 0;
    double split_2_log_prob;

    Spore* merge = 0;
    double merge_log_prob;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        spore = params->alternaria_proposal->get_random_spore(level);

        const Apical_hypha_density* hypha_d = params->alternaria_proposal
                                                    ->get_apical_hypha_d();
        const Spore_density* spore_d = params->alternaria_proposal
                                             ->get_spore_d();

        propose_1st_split(&split_1, &split_1_log_prob, hypha_d, spore);
        propose_2nd_split(&split_2, &split_2_log_prob, split_1);

        propose_merge(&merge, &merge_log_prob, spore->get_parent(), spore_d,
                split_1, split_2);

        params->alternaria_proposal->split_spore_into_2_apical(spore, 
                split_1, split_2);
    }
    catch (Exception e)
    {
        delete merge;
        delete split_1;
        delete split_2;
        return false;
    }

    delete merge;
    delete split_1;
    delete split_2;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double split_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double merge_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += merge_log_prob;
    log_alpha -= split_1_log_prob;
    log_alpha -= split_2_log_prob;
    log_alpha -= split_move_log_prob;
    log_alpha += merge_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


bool Spore_split_merge_2_apical_move::run_2(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Apical_hypha* hypha_1;
    Apical_hypha* hypha_2;

    Spore* merge = 0;
    double merge_log_prob;

    Apical_hypha* split_1 = 0;
    double split_1_log_prob;

    Apical_hypha* split_2 = 0;
    double split_2_log_prob;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        hypha_1 = params->alternaria_proposal->get_random_apical_hypha(level);
        if (!(hypha_2 = dynamic_cast<Apical_hypha*>(hypha_1->get_apical())))
        {
            return false;
        }

        const Apical_hypha_density* hypha_d = params->alternaria_proposal
                                                    ->get_apical_hypha_d();
        const Spore_density* spore_d = params->alternaria_proposal
                                             ->get_spore_d();

        propose_merge(&merge, &merge_log_prob, hypha_1->get_parent(), spore_d,
                hypha_1, hypha_2);
        
        propose_1st_split(&split_1, &split_1_log_prob, hypha_d, merge);
        propose_2nd_split(&split_2, &split_2_log_prob, split_1);

        params->alternaria_proposal->merge_2_apical_with_spore(hypha_1, merge);
    }
    catch (Exception e)
    {
        delete split_1;
        delete split_2;
        delete merge;
        return false;
    }

    delete split_1;
    delete split_2;
    delete merge;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double split_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double merge_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += split_1_log_prob;
    log_alpha += split_2_log_prob;
    log_alpha -= merge_log_prob;
    log_alpha += split_move_log_prob;
    log_alpha -= merge_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


void Spore_split_merge_2_apical_move::propose_1st_split
(
    Apical_hypha**              proposal_out,
    double*                     log_prob_out,
    const Apical_hypha_density* density,
    const Spore*                spore
)
throw (Arg_error)
{
    Apical_hypha* proposal;
    double        log_prob;

    assert(*proposal_out == 0);

    *proposal_out = density->sample();
    proposal = *proposal_out;
    log_prob = proposal->get_log_prob();

    float  mu, sigma;
    float  min, max;
    double p_1, p_2;
    double delta;

    const Structure* parent = spore->get_parent();

    // Set a parent-dependant width.
    if (parent && (dynamic_cast<const Apical_hypha*>(parent) || 
                   dynamic_cast<const Lateral_hypha*>(parent)))
    {
        float  width = proposal->get_width();

        // Take out the width log prob.
        mu        = density->get_width().get_mu();
        sigma     = density->get_width().get_sigma();
        min       = density->get_width().get_min();
        max       = density->get_width().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

        float dwidth = density->sample_dwidth();
        proposal->set_width(parent->get_width() + dwidth);

        // Put in the differential width log prob.
        mu        = density->get_dwidth().get_mu();
        sigma     = density->get_dwidth().get_sigma();
        min       = density->get_dwidth().get_min();
        max       = density->get_dwidth().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, dwidth, dwidth+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);
    }

    float theta = proposal->get_theta();

    // Take out the theta log prob.
    mu        = density->get_theta().get_mu();
    sigma     = density->get_theta().get_sigma();
    min       = density->get_theta().get_min();
    max       = density->get_theta().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, theta, theta+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a theta close to the original.
    mu    = spore->get_theta();
    sigma = density->get_theta().get_sigma();
    min   = density->get_theta().get_min();
    max   = density->get_theta().get_max();
    theta = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_theta(theta);

    // Put in the re-sampled theta log prob.
    mu        = density->get_theta().get_mu();
    sigma     = density->get_theta().get_sigma();
    min       = density->get_theta().get_min();
    max       = density->get_theta().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, theta, theta+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float psi;

    // Take out the psi log prob.
    min       = density->get_psi().get_min();
    max       = density->get_psi().get_max();
    p_1       = max - min;
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);

    // Re-sample a psi close to the original (using a Gaussian).
    mu    = spore->get_psi();
    sigma = density->get_theta().get_sigma(); // Not a bug.
    min   = density->get_psi().get_min();
    max   = density->get_psi().get_max();
    psi   = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_psi(psi);

    // Put in the re-sampled psi log prob.
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, psi, psi+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float opacity = proposal->get_opacity();

    // Take out the opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a opacity close to the original.
    mu    = spore->get_opacity();
    sigma = density->get_opacity().get_sigma();
    min   = density->get_opacity().get_min();
    max   = density->get_opacity().get_max();
    opacity = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_opacity(opacity);

    // Put in the re-sampled opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *log_prob_out = log_prob;
}


void Spore_split_merge_2_apical_move::propose_2nd_split
(
    Apical_hypha**      proposal_out,
    double*             log_prob_out,
    const Apical_hypha* split_1
)
throw (Arg_error)
{
    Apical_hypha* proposal;
    double        log_prob;

    assert(*proposal_out == 0);

    const Apical_hypha_density* density = split_1->get_density();

    *proposal_out = density->sample();
    proposal = *proposal_out;
    log_prob = proposal->get_log_prob();

    float  mu, sigma;
    float  min, max;
    double p_1, p_2;
    double delta;
    float  width = proposal->get_width();

    // Take out the width log prob.
    mu        = density->get_width().get_mu();
    sigma     = density->get_width().get_sigma();
    min       = density->get_width().get_min();
    max       = density->get_width().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float dwidth = density->sample_dwidth();
    proposal->set_width(split_1->get_width() + dwidth);

    // Put in the differential width log prob.
    mu        = density->get_dwidth().get_mu();
    sigma     = density->get_dwidth().get_sigma();
    min       = density->get_dwidth().get_min();
    max       = density->get_dwidth().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, dwidth, dwidth+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float opacity = proposal->get_opacity();

    // Take out the opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a opacity close to the original.
    mu    = split_1->get_opacity();
    sigma = density->get_opacity().get_sigma();
    min   = density->get_opacity().get_min();
    max   = density->get_opacity().get_max();
    opacity = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_opacity(opacity);

    // Put in the re-sampled opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *log_prob_out = log_prob;
}


void Spore_split_merge_2_apical_move::propose_merge
(
    Spore**              proposal_out,
    double*              log_prob_out,
    const Structure*     parent,
    const Spore_density* density,
    const Apical_hypha*  hypha_1,
    const Apical_hypha*  hypha_2
)
throw (Arg_error)
{
    Spore* proposal;
    double log_prob;

    assert(*proposal_out == 0);

    *proposal_out = density->sample();
    proposal = *proposal_out;
    log_prob = proposal->get_log_prob();

    float  mu, sigma;
    float  min, max;
    double p_1, p_2;
    double delta;

    float opacity = proposal->get_opacity();

    // Take out the opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a opacity close to the original.
    opacity = 0.5f*hypha_1->get_opacity() + 0.5f*hypha_2->get_opacity();
    sigma = density->get_opacity().get_sigma();
    min   = -6.0f*sigma;
    max   = 6.0f*sigma;
    opacity += sample_gaussian_pdf_d(0, sigma, min, max);
    proposal->set_opacity(opacity);

    // Put in the re-sampled opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *proposal_out = proposal;
    *log_prob_out = log_prob;
}

bool Spore_split_merge_spore_move::run_1
(
    Sampler_move_parameters* params
) 
throw (Exception)
{
    params->init_proposals();

    Spore* spore;

    Spore* split_1 = 0;
    double split_1_log_prob;

    Spore* split_2 = 0;
    double split_2_log_prob;

    Spore* merge = 0;
    double merge_log_prob;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        spore = params->alternaria_proposal->get_random_spore(level);

        propose_1st_split(&split_1, &split_1_log_prob, spore);
        propose_2nd_split(&split_2, &split_2_log_prob, split_1);

        propose_merge(&merge, &merge_log_prob, spore->get_parent(), split_1, 
                split_2);

        params->alternaria_proposal->split_spore_into_spore(spore, 
                split_1, split_2);
    }
    catch (Exception e)
    {
        delete merge;
        delete split_1;
        delete split_2;
        return false;
    }

    delete merge;
    delete split_1;
    delete split_2;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double split_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double merge_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += merge_log_prob;
    log_alpha -= split_1_log_prob;
    log_alpha -= split_2_log_prob;
    log_alpha -= split_move_log_prob;
    log_alpha += merge_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


bool Spore_split_merge_spore_move::run_2
(
    Sampler_move_parameters* params
) 
throw (Exception)
{
    params->init_proposals();

    Spore* spore;
    Apical_structure* apical;

    Spore* merge = 0;
    double merge_log_prob;

    Spore* split_1 = 0;
    double split_1_log_prob;

    Spore* split_2 = 0;
    double split_2_log_prob;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        spore = params->alternaria_proposal->get_random_spore(level);
        apical = spore->get_apical();

        if (!apical || !dynamic_cast<Spore*>(apical))
        {
            return false;
        }

        propose_merge(&merge, &merge_log_prob, spore->get_parent(), spore, 
                (Spore*)apical);
        
        propose_1st_split(&split_1, &split_1_log_prob, merge);
        propose_2nd_split(&split_2, &split_2_log_prob, split_1);

        params->alternaria_proposal->merge_spore_with_spore(spore, 
                merge);
    }
    catch (Exception e)
    {
        delete split_1;
        delete split_2;
        delete merge;
        return false;
    }

    delete split_1;
    delete split_2;
    delete merge;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double split_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double merge_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += split_1_log_prob;
    log_alpha += split_2_log_prob;
    log_alpha -= merge_log_prob;
    log_alpha += split_move_log_prob;
    log_alpha -= merge_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


void Spore_split_merge_spore_move::propose_1st_split
(
    Spore**      proposal_out,
    double*      log_prob_out,
    const Spore* spore
)
throw (Arg_error)
{
    Spore* proposal;
    double log_prob;

    assert(*proposal_out == 0);

    const Spore_density* density = spore->get_density();

    *proposal_out = density->sample();
    proposal = *proposal_out;
    log_prob = proposal->get_log_prob();

    float  mu, sigma;
    float  min, max;
    double p_1, p_2;
    double delta;

    float width = proposal->get_width();

    // Take out the width log prob.
    mu        = density->get_width().get_mu();
    sigma     = density->get_width().get_sigma();
    min       = density->get_width().get_min();
    max       = density->get_width().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a width close to the original.
    width = 0.65f*spore->get_width();
    sigma = density->get_width().get_sigma();
    min   = -6.0f*sigma;
    max   = 6.0f*sigma;
    width += sample_gaussian_pdf_d(0, sigma, min, max);
    proposal->set_width(width);

    // Put in the re-sampled width log prob.
    mu        = density->get_width().get_mu();
    sigma     = density->get_width().get_sigma();
    min       = density->get_width().get_min();
    max       = density->get_width().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float theta = proposal->get_theta();

    // Take out the theta log prob.
    mu        = density->get_theta().get_mu();
    sigma     = density->get_theta().get_sigma();
    min       = density->get_theta().get_min();
    max       = density->get_theta().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, theta, theta+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a theta close to the original.
    mu    = spore->get_theta();
    sigma = density->get_theta().get_sigma();
    min   = density->get_theta().get_min();
    max   = density->get_theta().get_max();
    theta = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_theta(theta);

    // Put in the re-sampled theta log prob.
    mu        = density->get_theta().get_mu();
    sigma     = density->get_theta().get_sigma();
    min       = density->get_theta().get_min();
    max       = density->get_theta().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, theta, theta+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float psi;

    // Take out the psi log prob.
    min       = density->get_psi().get_min();
    max       = density->get_psi().get_max();
    p_1       = max - min;
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);

    // Re-sample a psi close to the original (using a Gaussian).
    mu    = spore->get_psi();
    sigma = density->get_theta().get_sigma(); // Not a bug.
    min   = density->get_psi().get_min();
    max   = density->get_psi().get_max();
    psi   = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_psi(psi);

    // Put in the re-sampled psi log prob.
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, psi, psi+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float opacity = proposal->get_opacity();

    // Take out the opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a opacity close to the original.
    mu    = spore->get_opacity();
    sigma = density->get_opacity().get_sigma();
    min   = density->get_opacity().get_min();
    max   = density->get_opacity().get_max();
    opacity = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_opacity(opacity);

    // Put in the re-sampled opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *log_prob_out = log_prob;
}


void Spore_split_merge_spore_move::propose_2nd_split
(
    Spore**      proposal_out,
    double*      log_prob_out,
    const Spore* split_1
)
throw (Arg_error)
{
    Spore* proposal;
    double log_prob;

    assert(*proposal_out == 0);

    const Spore_density* density = split_1->get_density();

    *proposal_out = density->sample();
    proposal = *proposal_out;
    log_prob = proposal->get_log_prob();

    float  mu, sigma;
    float  min, max;
    double p_1, p_2;
    double delta;

    float opacity = proposal->get_opacity();

    // Take out the opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a opacity close to the original.
    mu    = split_1->get_opacity();
    sigma = density->get_opacity().get_sigma();
    min   = density->get_opacity().get_min();
    max   = density->get_opacity().get_max();
    opacity = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_opacity(opacity);

    // Put in the re-sampled opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *log_prob_out = log_prob;
}


void Spore_split_merge_spore_move::propose_merge
(
    Spore**         proposal_out,
    double*         log_prob_out,
    const Structure*parent,
    const Spore*    spore_1,
    const Spore*    spore_2
)
throw (Arg_error)
{
    Spore* proposal;
    double log_prob;

    assert(*proposal_out == 0);

    const Spore_density* density = spore_1->get_density();

    *proposal_out = density->sample();
    proposal = *proposal_out;
    log_prob = proposal->get_log_prob();

    float  mu, sigma;
    float  min, max;
    double p_1, p_2;
    double delta;

    float width = proposal->get_width();

    // Take out the width log prob.
    mu        = density->get_width().get_mu();
    sigma     = density->get_width().get_sigma();
    min       = density->get_width().get_min();
    max       = density->get_width().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a width close to the original.
    width = 0.65f*spore_1->get_width() + 0.65f*spore_2->get_width();
    sigma = density->get_width().get_sigma();
    min   = -6.0f*sigma;
    max   = 6.0f*sigma;
    width += sample_gaussian_pdf_d(0, sigma, min, max);
    proposal->set_width(width);

    // Put in the re-sampled width log prob.
    mu        = density->get_width().get_mu();
    sigma     = density->get_width().get_sigma();
    min       = density->get_width().get_min();
    max       = density->get_width().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float opacity = proposal->get_opacity();

    // Take out the opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a opacity close to the original.
    opacity = 0.5f*spore_1->get_opacity() + 0.5f*spore_2->get_opacity();
    sigma = density->get_opacity().get_sigma();
    min   = -6.0f*sigma;
    max   = 6.0f*sigma;
    opacity += sample_gaussian_pdf_d(0, sigma, min, max);
    proposal->set_opacity(opacity);

    // Put in the re-sampled opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);
#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *proposal_out = proposal;
    *log_prob_out = log_prob;
}

bool Lateral_hypha_birth_death_move::run_1(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Apical_structure* parent;
    Lateral_hypha* hypha;
    const Lateral_hypha_density* density;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        //parent = params->alternaria_proposal->get_random_apical_hypha(level);
        parent = dynamic_cast<Apical_structure*>(
                 params->alternaria_proposal->get_random_structure(level));

        if (!parent)
        {
            return false;
        }

        // HACK Do not allow terminal parents.
        if (!parent->get_apical())
        {
            return false;
        }

        density = (level == 0) 
                ? params->alternaria_proposal->get_lateral_hypha_1_d()
                : params->alternaria_proposal->get_lateral_hypha_n_d();

        hypha = density->sample(parent);
        params->alternaria_proposal->add_lateral_hypha(hypha);
    }
    catch (Exception e)
    {
        return false;
    }

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double hypha_log_prob = hypha->get_log_prob();

    double birth_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double death_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha -= hypha_log_prob;
    log_alpha -= birth_move_log_prob;
    log_alpha += death_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


bool Lateral_hypha_birth_death_move::run_2(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Lateral_hypha* hypha;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        hypha = params->alternaria_proposal
                      ->get_random_terminal_lateral_hypha(level);

        // Check that hypha doesn't have any laterals.
        if (hypha->get_laterals().size() != 0)
        {
            return false;
        }
    }
    catch (Exception e)
    {
        return false;
    }

    double hypha_log_prob = hypha->get_log_prob();

    try
    {
        params->alternaria_proposal->remove_lateral_hypha(hypha);
    }
    catch (Arg_error e)
    {
        return false;
    }

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double birth_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double death_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += hypha_log_prob;
    log_alpha += birth_move_log_prob;
    log_alpha -= death_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


bool Apical_hypha_split_merge_apical_move::run_1
(
    Sampler_move_parameters* params
) 
throw (Exception)
{
    params->init_proposals();

    Apical_hypha* hypha;

    Apical_hypha* split_1 = 0;
    double split_1_log_prob;

    Apical_hypha* split_2 = 0;
    double split_2_log_prob;

    Apical_hypha* merge = 0;
    double merge_log_prob;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        hypha = params->alternaria_proposal->get_random_apical_hypha(level);

        propose_1st_split(&split_1, &split_1_log_prob, hypha);
        propose_2nd_split(&split_2, &split_2_log_prob, split_1);

        propose_merge(&merge, &merge_log_prob, hypha->get_parent(), split_1, 
                split_2);

        params->alternaria_proposal->split_apical_hypha_into_apical(hypha, 
                split_1, split_2);
    }
    catch (Exception e)
    {
        delete merge;
        delete split_1;
        delete split_2;
        return false;
    }

    delete merge;
    delete split_1;
    delete split_2;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double split_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double merge_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += merge_log_prob;
    log_alpha -= split_1_log_prob;
    log_alpha -= split_2_log_prob;
    log_alpha -= split_move_log_prob;
    log_alpha += merge_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


bool Apical_hypha_split_merge_apical_move::run_2
(
    Sampler_move_parameters* params
) 
throw (Exception)
{
    params->init_proposals();

    Apical_hypha* hypha;
    Apical_structure* apical;

    Apical_hypha* merge = 0;
    double merge_log_prob;

    Apical_hypha* split_1 = 0;
    double split_1_log_prob;

    Apical_hypha* split_2 = 0;
    double split_2_log_prob;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        hypha = params->alternaria_proposal->get_random_apical_hypha(level);
        apical = hypha->get_apical();

        if (!apical || !dynamic_cast<Apical_hypha*>(apical))
        {
            return false;
        }

        propose_merge(&merge, &merge_log_prob, hypha->get_parent(), hypha, 
                (Apical_hypha*)apical);
        
        propose_1st_split(&split_1, &split_1_log_prob, merge);
        propose_2nd_split(&split_2, &split_2_log_prob, split_1);

        params->alternaria_proposal->merge_apical_hypha_with_apical(hypha, 
                merge);
    }
    catch (Exception e)
    {
        delete split_1;
        delete split_2;
        delete merge;
        return false;
    }

    delete split_1;
    delete split_2;
    delete merge;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double split_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double merge_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += split_1_log_prob;
    log_alpha += split_2_log_prob;
    log_alpha -= merge_log_prob;
    log_alpha += split_move_log_prob;
    log_alpha -= merge_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


void Apical_hypha_split_merge_apical_move::propose_1st_split
(
    Apical_hypha**      proposal_out,
    double*             log_prob_out,
    const Apical_hypha* hypha
)
throw (Arg_error)
{
    Apical_hypha* proposal;
    double        log_prob;

    assert(*proposal_out == 0);

    const Apical_hypha_density* density = hypha->get_density();

    *proposal_out = density->sample();
    proposal = *proposal_out;
    log_prob = proposal->get_log_prob();

    float  mu, sigma;
    float  min, max;
    double p_1, p_2;
    double delta;

    const Structure* parent = hypha->get_parent();

    // Set a parent-dependant width.
    if (parent && (dynamic_cast<const Apical_hypha*>(parent) || 
                   dynamic_cast<const Lateral_hypha*>(parent)))
    {
        float  width = proposal->get_width();

        // Take out the width log prob.
        mu        = density->get_width().get_mu();
        sigma     = density->get_width().get_sigma();
        min       = density->get_width().get_min();
        max       = density->get_width().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

        float dwidth = density->sample_dwidth();
        proposal->set_width(parent->get_width() + dwidth);

        // Put in the differential width log prob.
        mu        = density->get_dwidth().get_mu();
        sigma     = density->get_dwidth().get_sigma();
        min       = density->get_dwidth().get_min();
        max       = density->get_dwidth().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, dwidth, dwidth+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);
    }

    float theta = proposal->get_theta();

    // Take out the theta log prob.
    mu        = density->get_theta().get_mu();
    sigma     = density->get_theta().get_sigma();
    min       = density->get_theta().get_min();
    max       = density->get_theta().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, theta, theta+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a theta close to the original.
    mu    = hypha->get_theta();
    sigma = density->get_theta().get_sigma();
    min   = density->get_theta().get_min();
    max   = density->get_theta().get_max();
    theta = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_theta(theta);

    // Put in the re-sampled theta log prob.
    mu        = density->get_theta().get_mu();
    sigma     = density->get_theta().get_sigma();
    min       = density->get_theta().get_min();
    max       = density->get_theta().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, theta, theta+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float psi;

    // Take out the psi log prob.
    min       = density->get_psi().get_min();
    max       = density->get_psi().get_max();
    p_1       = max - min;
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);

    // Re-sample a psi close to the original (using a Gaussian).
    mu    = hypha->get_psi();
    sigma = density->get_theta().get_sigma(); // Not a bug.
    min   = density->get_psi().get_min();
    max   = density->get_psi().get_max();
    psi   = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_psi(psi);

    // Put in the re-sampled psi log prob.
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, psi, psi+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float opacity = proposal->get_opacity();

    // Take out the opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a opacity close to the original.
    mu    = hypha->get_opacity();
    sigma = density->get_opacity().get_sigma();
    min   = density->get_opacity().get_min();
    max   = density->get_opacity().get_max();
    opacity = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_opacity(opacity);

    // Put in the re-sampled opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *log_prob_out = log_prob;
}


void Apical_hypha_split_merge_apical_move::propose_2nd_split
(
    Apical_hypha**      proposal_out,
    double*             log_prob_out,
    const Apical_hypha* split_1
)
throw (Arg_error)
{
    Apical_hypha* proposal;
    double        log_prob;

    assert(*proposal_out == 0);

    const Apical_hypha_density* density = split_1->get_density();

    *proposal_out = density->sample();
    proposal = *proposal_out;
    log_prob = proposal->get_log_prob();

    float  mu, sigma;
    float  min, max;
    double p_1, p_2;
    double delta;
    float  width = proposal->get_width();

    // Take out the width log prob.
    mu        = density->get_width().get_mu();
    sigma     = density->get_width().get_sigma();
    min       = density->get_width().get_min();
    max       = density->get_width().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float dwidth = density->sample_dwidth();
    proposal->set_width(split_1->get_width() + dwidth);

    // Put in the differential width log prob.
    mu        = density->get_dwidth().get_mu();
    sigma     = density->get_dwidth().get_sigma();
    min       = density->get_dwidth().get_min();
    max       = density->get_dwidth().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, dwidth, dwidth+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float opacity = proposal->get_opacity();

    // Take out the opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a opacity close to the original.
    mu    = split_1->get_opacity();
    sigma = density->get_opacity().get_sigma();
    min   = density->get_opacity().get_min();
    max   = density->get_opacity().get_max();
    opacity = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_opacity(opacity);

    // Put in the re-sampled opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *log_prob_out = log_prob;
}


void Apical_hypha_split_merge_apical_move::propose_merge
(
    Apical_hypha**      proposal_out,
    double*             log_prob_out,
    const Structure*    parent,
    const Apical_hypha* hypha_1,
    const Apical_hypha* hypha_2
)
throw (Arg_error)
{
    Apical_hypha* proposal;
    double        log_prob;

    assert(*proposal_out == 0);

    const Apical_hypha_density* density = hypha_1->get_density();

    *proposal_out = density->sample();
    proposal = *proposal_out;
    log_prob = proposal->get_log_prob();

    float  mu, sigma;
    float  min, max;
    double p_1, p_2;
    double delta;

    // Set a parent-dependant width.
    if (parent && (dynamic_cast<const Apical_hypha*>(parent) || 
                   dynamic_cast<const Lateral_hypha*>(parent)))
    {
        float  width = proposal->get_width();

        // Take out the width log prob.
        mu        = density->get_width().get_mu();
        sigma     = density->get_width().get_sigma();
        min       = density->get_width().get_min();
        max       = density->get_width().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

        float dwidth = density->sample_dwidth();
        proposal->set_width(parent->get_width() + dwidth);

        // Put in the differential width log prob.
        mu        = density->get_dwidth().get_mu();
        sigma     = density->get_dwidth().get_sigma();
        min       = density->get_dwidth().get_min();
        max       = density->get_dwidth().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, dwidth, dwidth+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);
    }

    float opacity = proposal->get_opacity();

    // Take out the opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a opacity close to the original.
    opacity = 0.5f*hypha_1->get_opacity() + 0.5f*hypha_2->get_opacity();
    sigma = density->get_opacity().get_sigma();
    min   = -6.0f*sigma;
    max   = 6.0f*sigma;
    opacity += sample_gaussian_pdf_d(0, sigma, min, max);
    proposal->set_opacity(opacity);

    // Put in the re-sampled opacity log prob.
    mu        = density->get_opacity().get_mu();
    sigma     = density->get_opacity().get_sigma();
    min       = density->get_opacity().get_min();
    max       = density->get_opacity().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, opacity, opacity+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *proposal_out = proposal;
    *log_prob_out = log_prob;
}


bool Apical_hypha_split_merge_lateral_move::run_1
(
    Sampler_move_parameters* params
) 
throw (Exception)
{
    params->init_proposals();

    Apical_hypha* hypha;

    Apical_hypha* split_1 = 0;
    double split_1_log_prob;

    Lateral_hypha* split_2 = 0;
    double split_2_log_prob;

    Apical_hypha* merge = 0;
    double merge_log_prob;

    const Lateral_hypha_density* lateral_d;

    // This move sucks!
    return false;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        hypha = params->alternaria_proposal->get_random_apical_hypha(level);

        // HACK Do not split a terminal.
        if (!hypha->get_apical())
        {
            return false;
        }

        lateral_d = (level == 0) 
                ? params->alternaria_proposal->get_lateral_hypha_1_d()
                : params->alternaria_proposal->get_lateral_hypha_n_d();

        propose_1st_split(&split_1, &split_1_log_prob, hypha);
        propose_2nd_split(&split_2, &split_2_log_prob, split_1, lateral_d);

        propose_merge(&merge, &merge_log_prob, hypha->get_parent(), split_1, 
                split_2);

        params->alternaria_proposal->split_apical_hypha_into_lateral(hypha, 
                split_1, split_2);
    }
    catch (Exception e)
    {
        delete merge;
        delete split_1;
        delete split_2;
        return false;
    }

    delete split_1;
    delete split_2;
    delete merge;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double split_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double merge_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += merge_log_prob;
    log_alpha -= split_1_log_prob;
    log_alpha -= split_2_log_prob;
    log_alpha -= split_move_log_prob;
    log_alpha += merge_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


bool Apical_hypha_split_merge_lateral_move::run_2
(
    Sampler_move_parameters* params
) 
throw (Exception)
{
    params->init_proposals();

    Apical_hypha* hypha;
    Lateral_structure* lateral;

    Apical_hypha* merge = 0;
    double merge_log_prob;

    Apical_hypha* split_1 = 0;
    double split_1_log_prob;

    Lateral_hypha* split_2 = 0;
    double split_2_log_prob;

    const Lateral_hypha_density* lateral_d;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        hypha = params->alternaria_proposal
                      ->get_random_terminal_apical_hypha(level);

        lateral = hypha->get_random_lateral();

        if (!dynamic_cast<Lateral_hypha*>(lateral))
        {
            return false;
        }

        propose_merge(&merge, &merge_log_prob, hypha->get_parent(), hypha, 
                (Lateral_hypha*)lateral);

        lateral_d = (level == 0) 
                ? params->alternaria_proposal->get_lateral_hypha_1_d()
                : params->alternaria_proposal->get_lateral_hypha_n_d();
        
        propose_1st_split(&split_1, &split_1_log_prob, merge);
        propose_2nd_split(&split_2, &split_2_log_prob, split_1, lateral_d);

        params->alternaria_proposal->merge_apical_hypha_with_lateral(hypha,
                (Lateral_hypha*)lateral, merge);
    }
    catch (Exception e)
    {
        delete split_1;
        delete split_2;
        delete merge;
        return false;
    }

    delete split_1;
    delete split_2;
    delete merge;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double split_move_log_prob = (prob_1 >= JWSC_MIN_LOG_ARG) ? log(prob_1) 
        : log(JWSC_MIN_LOG_ARG);

    double merge_move_log_prob = (prob_2 >= JWSC_MIN_LOG_ARG) ? log(prob_2) 
        : log(JWSC_MIN_LOG_ARG);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();
    log_alpha += split_1_log_prob;
    log_alpha += split_2_log_prob;
    log_alpha -= merge_log_prob;
    log_alpha += split_move_log_prob;
    log_alpha -= merge_move_log_prob;

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


void Apical_hypha_split_merge_lateral_move::propose_1st_split
(
    Apical_hypha**      proposal_out,
    double*             log_prob_out,
    const Apical_hypha* hypha
)
throw (Arg_error)
{
    Apical_hypha* proposal;
    double        log_prob;

    assert(*proposal_out == 0);

    const Apical_hypha_density* density = hypha->get_density();

    *proposal_out = density->sample();
    proposal = *proposal_out;
    log_prob = proposal->get_log_prob();

    float  mu, sigma;
    float  min, max;
    double p_1, p_2;
    double delta;

    const Structure* parent = hypha->get_parent();

    // Set a parent-dependant width.
    if (parent && (dynamic_cast<const Apical_hypha*>(parent) || 
                   dynamic_cast<const Lateral_hypha*>(parent)))
    {
        float  width = proposal->get_width();

        // Take out the width log prob.
        mu        = density->get_width().get_mu();
        sigma     = density->get_width().get_sigma();
        min       = density->get_width().get_min();
        max       = density->get_width().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

        float dwidth = density->sample_dwidth();
        proposal->set_width(parent->get_width() + dwidth);

        // Put in the differential width log prob.
        mu        = density->get_dwidth().get_mu();
        sigma     = density->get_dwidth().get_sigma();
        min       = density->get_dwidth().get_min();
        max       = density->get_dwidth().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, dwidth, dwidth+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);
    }

    float theta = proposal->get_theta();

    // Take out the theta log prob.
    mu        = density->get_theta().get_mu();
    sigma     = density->get_theta().get_sigma();
    min       = density->get_theta().get_min();
    max       = density->get_theta().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, theta, theta+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    // Re-sample a theta close to the original.
    mu    = hypha->get_theta();
    sigma = density->get_theta().get_sigma();
    min   = density->get_theta().get_min();
    max   = density->get_theta().get_max();
    theta = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_theta(theta);

    // Put in the re-sampled theta log prob.
    mu        = density->get_theta().get_mu();
    sigma     = density->get_theta().get_sigma();
    min       = density->get_theta().get_min();
    max       = density->get_theta().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, theta, theta+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    float psi;

    // Take out the psi log prob.
    min       = density->get_psi().get_min();
    max       = density->get_psi().get_max();
    p_1       = max - min;
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);

    // Re-sample a psi close to the original (using a Gaussian).
    mu    = hypha->get_psi();
    sigma = density->get_theta().get_sigma(); // Not a bug.
    min   = density->get_psi().get_min();
    max   = density->get_psi().get_max();
    psi   = sample_gaussian_pdf_d(mu, sigma, min, max);
    proposal->set_psi(psi);

    // Put in the re-sampled psi log prob.
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, psi, psi+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *log_prob_out = log_prob;
}


void Apical_hypha_split_merge_lateral_move::propose_2nd_split
(
    Lateral_hypha**              proposal_out,
    double*                      log_prob_out,
    const Apical_hypha*          split_1,
    const Lateral_hypha_density* density
)
{
    assert(*proposal_out == 0);

    Lateral_hypha* proposal = density->sample();
    double log_prob = proposal->get_log_prob();

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *proposal_out = proposal;
    *log_prob_out = log_prob;
}


void Apical_hypha_split_merge_lateral_move::propose_merge
(
    Apical_hypha**       proposal_out,
    double*              log_prob_out,
    const Structure*     parent,
    const Apical_hypha*  hypha_1,
    const Lateral_hypha* hypha_2
)
throw (Arg_error)
{
    Apical_hypha* proposal;
    double        log_prob;

    assert(*proposal_out == 0);

    const Apical_hypha_density* density = hypha_1->get_density();

    *proposal_out = density->sample();
    proposal = *proposal_out;
    log_prob = proposal->get_log_prob();

    // Set a parent-dependant width.
    if (parent && (dynamic_cast<const Apical_hypha*>(parent) || 
                   dynamic_cast<const Lateral_hypha*>(parent)))
    {
        float  mu, sigma;
        float  min, max;
        double p_1, p_2;
        double delta;
        float  width = proposal->get_width();

        // Take out the width log prob.
        mu        = density->get_width().get_mu();
        sigma     = density->get_width().get_sigma();
        min       = density->get_width().get_min();
        max       = density->get_width().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, width, width+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob -= log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob += log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

        float dwidth = density->sample_dwidth();
        proposal->set_width(parent->get_width() + dwidth);

        // Put in the differential width log prob.
        mu        = density->get_dwidth().get_mu();
        sigma     = density->get_dwidth().get_sigma();
        min       = density->get_dwidth().get_min();
        max       = density->get_dwidth().get_max();
        delta     = (max - min) / 1000;
        p_1       = integrate_gaussian_pdf_d(mu, sigma, dwidth, dwidth+delta); 
        p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
        log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
        log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);
    }

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!std::isinf(log_prob) && !std::isnan(log_prob));
#endif

    *log_prob_out = log_prob;
}


bool Structure_resize_move::run(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Structure* s;
    Structure* parent;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        s = params->alternaria_proposal->get_random_structure(level);
        parent = s->get_parent();
    }
    catch (Exception e)
    {
        return false;
    }

    Structure_density* density = s->get_density()->clone();

    float mu    = s->get_length();
    float sigma = density->get_length().get_sigma();
    float min   = density->get_length().get_min();
    float max   = density->get_length().get_max();
    density->set_length(mu, sigma, min, max);

    float width;

    if (parent && (dynamic_cast<Apical_hypha*>(parent) || 
                dynamic_cast<Lateral_hypha*>(parent)) && 
            dynamic_cast<Apical_hypha_density*>(density))
    {
        Apical_hypha_density* hypha_density = (Apical_hypha_density*)density;

        width = parent->get_width() + hypha_density->sample_dwidth();
    }
    else
    {
        mu    = s->get_width();
        sigma = density->get_width().get_sigma();
        min   = density->get_width().get_min();
        max   = density->get_width().get_max();
        density->set_width(mu, sigma, min, max);

        width = density->sample_width();
    }

    try
    {
        s->set_size(density->sample_length(), width);
    }
    catch (Arg_error e)
    {
        delete density;
        return false;
    }

    delete density;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


const std::ostringstream& Structure_resize_move::get_info
(
    Sampler_move_parameters* params
) 
{
    info.str("");
    return info; 
}


const std::ostringstream& Structure_resize_move::get_proposal_info
(
    Sampler_move_parameters* params
)
{ 
    proposal_info.str("");
    return proposal_info; 
}


bool Structure_rotate_move::run(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Structure* s;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        s = params->alternaria_proposal->get_random_structure(level);
    }
    catch (Exception e)
    {
        return false;
    }

    Structure_density* density = s->get_density()->clone();

    float mu    = s->get_theta();
    float sigma = density->get_theta().get_sigma();
    float min   = density->get_theta().get_min();
    float max   = density->get_theta().get_max();
    density->set_theta(mu, sigma, min, max);
    float theta = density->sample_theta();

    /*
    mu    = s->get_psi();
    sigma = density->get_theta().get_sigma(); // Not a bug, but dumb.
    min   = density->get_psi().get_min();
    max   = density->get_psi().get_max();
    float psi = sample_gaussian_pdf_d(mu, sigma, min, max);
    */
    mu    = 0;
    sigma = density->get_theta().get_sigma(); // Not a bug, but dumb.
    min   = -6.0f*sigma;
    max   = 6.0f*sigma;
    float psi = s->get_psi() + sample_gaussian_pdf_d(mu, sigma, min, max);

    min = density->get_psi().get_min();
    max = density->get_psi().get_max();

    if (psi > max && (psi - JWSC_2_PI) >= min)
    {
        psi -= JWSC_2_PI;
    }
    else if (psi < min && (psi + JWSC_2_PI) <= max)
    {
        psi += JWSC_2_PI;
    }

    try
    {
        s->set_orientation(theta, psi);
    }
    catch (Arg_error e)
    {
        delete density;
        return false;
    }

    delete density;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


const std::ostringstream& Structure_rotate_move::get_info
(
    Sampler_move_parameters* params
) 
{ 
    info.str("");
    return info; 
}


const std::ostringstream& Structure_rotate_move::get_proposal_info
(
    Sampler_move_parameters* params
)
{ 
    proposal_info.str("");
    return proposal_info;
}


bool Structure_opacity_move::run(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Structure* s;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        s = params->alternaria_proposal->get_random_structure(level);
    }
    catch (Exception e)
    {
        return false;
    }

    Structure_density* density = s->get_density()->clone();

    float mu    = s->get_opacity();
    float sigma = density->get_opacity().get_sigma();
    float min   = density->get_opacity().get_min();
    float max   = density->get_opacity().get_max();
    density->set_opacity(mu, sigma, min, max);

    s->set_opacity(density->sample_opacity());

    delete density;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


const std::ostringstream& Structure_opacity_move::get_info
(
    Sampler_move_parameters* params
) 
{ 
    info.str("");
    return info; 
}


const std::ostringstream& Structure_opacity_move::get_proposal_info
(
    Sampler_move_parameters* params
)
{ 
    proposal_info.str("");
    return proposal_info; 
}


bool Apical_structure_position_move::run(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Apical_structure* s;

    try
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        s = dynamic_cast<Apical_structure*>(
                params->alternaria_proposal->get_random_structure(level));

        if (!s)
        {
            return false;
        }
    }
    catch (Exception e)
    {
        return false;
    }

    const Vector_f* centroid = s->get_centroid();
    float x = centroid->elts[0] + sample_gaussian_pdf_d(0, 1, -10, 10);
    float y = centroid->elts[1] + sample_gaussian_pdf_d(0, 1, -10, 10);
    float z = centroid->elts[2] + sample_gaussian_pdf_d(0, 1, -10, 10);

    try
    {
        s->set_position(x, y, z);
    }
    catch (Arg_error e)
    {
        return false;
    }

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


const std::ostringstream& Apical_structure_position_move::get_info
(
    Sampler_move_parameters* params
) 
{ 
    info.str("");
    return info; 
}


const std::ostringstream& Apical_structure_position_move::get_proposal_info
(
    Sampler_move_parameters* params
)
{ 
    proposal_info.str("");
    return proposal_info; 
}


bool Lateral_structure_lat_dist_move::run(Sampler_move_parameters* params) 
    throw (Exception)
{
    params->init_proposals();

    Lateral_structure* lat;

    try 
    {
        size_t level = params->alternaria_proposal->get_uniform_random_level();
        lat = params->alternaria_proposal->get_random_lateral_hypha(level);
    }
    catch (Exception e)
    {
        return false;
    }

    Structure* parent = 0;
    Apical_structure* apical = 0;

    if (lat->get_parent())
    {
        parent = lat->get_parent()->get_parent();
        apical = lat->get_parent()->get_apical();
    }

    Lateral_structure_density* density = lat->get_density()->clone();

    try
    {
        float mu    = lat->get_lat_dist();
        float sigma = density->get_lat_dist().get_sigma();
        float min   = density->get_lat_dist().get_min();
        float max   = density->get_lat_dist().get_max();
        float lat_dist;

        if (dynamic_cast<Apical_structure*>(lat->get_parent()) 
                && parent && apical)
        {
            float t_mu  = mu - min;
            float t_max = 2.0f*(max - min);
            float t_min = -max + min;

            lat_dist = sample_gaussian_pdf_d(t_mu, sigma, t_min, t_max);

            if (lat_dist < 0)
            {
                lat->set_parent(parent, -lat_dist + min);
            }
            else if (lat_dist > (max - min))
            {
                lat->set_parent(apical, (lat_dist - (max - min)) + min);
            }
            else
            {
                lat->set_lat_dist(lat_dist + min);
            }
        }
        else if (dynamic_cast<Apical_structure*>(lat->get_parent()) 
                && parent && !apical)
        {
            float t_mu  = mu - min;
            float t_max = max - min;
            float t_min = -max + min;

            lat_dist = sample_gaussian_pdf_d(t_mu, sigma, t_min, t_max);

            if (lat_dist < 0)
            {
                lat->set_parent(parent, -lat_dist + min);
            }
            else
            {
                lat->set_lat_dist(lat_dist + min);
            }
        }
        else if ((dynamic_cast<Lateral_structure*>(lat->get_parent()) || 
                    !parent) && apical)
        {
            float t_mu  = mu - min;
            float t_max = 2.0f*(max - min);
            float t_min = 0;

            lat_dist = sample_gaussian_pdf_d(t_mu, sigma, t_min, t_max);

            if (lat_dist > (max - min))
            {
                lat->set_parent(apical, (lat_dist - (max - min)) + min);
            }
            else
            {
                lat->set_lat_dist(lat_dist + min);
            }
        }
        else if ((dynamic_cast<Lateral_structure*>(lat->get_parent()) || 
                    !parent) && !apical)
        {
            density->set_lat_dist(mu, sigma, min, max);
            lat->set_lat_dist(density->sample_lat_dist());
        }
        else
        {
            std::cerr << "BUG: Lateral_structure_lat_dist_move::run: invalid case\n";
            assert(0);
        }
    }
    catch (Arg_error e)
    {
        delete density;
        return false;
    }

    delete density;

    if (params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal))
    {
        return false;
    }

    params->ll_proposal = params->alternaria_proposal->calc_log_likelihood(
            params->psf_proposal, params->imaging_proposal, params->data,
            params->data_avg);

    double log_alpha = params->ll_proposal - params->ll;
    log_alpha += params->alternaria_proposal->get_log_prob();
    log_alpha -= params->alternaria->get_log_prob();

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_alpha) && !isnan(log_alpha));
#endif

    double log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

    if (log_u <= log_alpha)
    {
        return true;
    }
    return false;
}


const std::ostringstream& Lateral_structure_lat_dist_move::get_info
(
    Sampler_move_parameters* params
) 
{ 
    info.str("");
    return info;
}


const std::ostringstream& Lateral_structure_lat_dist_move::get_proposal_info
(
    Sampler_move_parameters* params
)
{ 
    proposal_info.str("");
    return proposal_info; 
}