alternaria_model.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 <cassert>
#include <cmath>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <list>
#include <vector>
#include <algorithm>
#include <stdexcept>

#include <inttypes.h>

#include <jwsc/config.h>
#include <jwsc/base/limits.h>
#include <jwsc/prob/pmf.h>
#include <jwsc/prob/pdf.h>
#include <jwsc/matblock/matblock.h>
#include <jwsc/matblock/matblock_conv.h>
#include <jwsc/imgblock/imgblock.h>
#include <jwsc/imgblock/imgblock_io.h>

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

#include "density.h"
#include "printable.h"
#include "structure.h"
#include "spore.h"
#include "hypha.h"
#include "psf_model.h"
#include "imaging_model.h"
#include "alternaria_model.h"


#define  NUM_LEVELS_P  0.7f

#define  SPORE_P            0.3f
#define  APICAL_HYPHA_P     0.6f
#define  LATERAL_HYPHA_1_P  0.1f
#define  LATERAL_HYPHA_N_P  0.1f

#define  C_1_MIN          0.0f
#define  C_1_MAX          0.0f


using std::isinf;
using std::isnan;
using std::log;
using std::sqrt;
using std::fabs;
using std::list;
using std::vector;
using std::ifstream;
using std::ofstream;
using std::ostringstream;
using std::istringstream;
using namespace jwsc;
using jwscxx::base::IO_error;
using jwscxx::base::Arg_error;
using jwscxx::base::Result_error;


/* -------------------------  Alternaria_density  --------------------------- */

Alternaria_density::Alternaria_density() throw (Arg_error)
{
    set_num_levels_p(NUM_LEVELS_P);
    set_structure_p(SPORE_P, APICAL_HYPHA_P, LATERAL_HYPHA_1_P, 
            LATERAL_HYPHA_N_P);
    set_c_1(C_1_MIN, C_1_MAX);
}


void Alternaria_density::set_num_levels_p(float p) throw (Arg_error)
{
    if (p <= 0 || p >= 1)
    {
        throw Arg_error("Must be in (0,1)");
    }
    num_levels_p = p;
}


void Alternaria_density::set_structure_p
(
    float spore_p,
    float apical_hypha_p,
    float lateral_hypha_1_p,
    float lateral_hypha_n_p
) 
throw (Arg_error)
{
    if (spore_p < 0 || spore_p > 1 || 
        apical_hypha_p < 0 || apical_hypha_p > 1 || 
        lateral_hypha_1_p < 0 || lateral_hypha_1_p > 1 ||
        lateral_hypha_n_p < 0 || lateral_hypha_n_p > 1)
    {
        throw Arg_error("Must be in [0,1]");
    }
    this->spore_p = spore_p;
    this->apical_hypha_p = apical_hypha_p;
    this->lateral_hypha_1_p = lateral_hypha_1_p;
    this->lateral_hypha_n_p = lateral_hypha_n_p;
}


void Alternaria_density::set_c_1(float min, float max) throw (Arg_error)
{
    if (min < 0)
    {
        throw Arg_error("Minimum must be >= 0");
    }
    c_1.set(min, max);
}


float Alternaria_density::sample_c_1() const
{
    return (float)sample_uniform_pdf_d(c_1.get_min(), c_1.get_max());
}


Alternaria* Alternaria_density::sample
(
    size_t                       N,
    float                        centroid_x, 
    float                        centroid_y, 
    float                        centroid_z,
    float                        theta,
    float                        psi,
    size_t                       base_level,
    const Spore_density*         spore_d,
    const Apical_hypha_density*  apical_hypha_d,
    const Lateral_hypha_density* lateral_hypha_1_d,
    const Lateral_hypha_density* lateral_hypha_n_d
)
const throw (Arg_error)
{
    Alternaria* alternaria = new Alternaria(sample_c_1(), theta,
            psi, base_level, this, spore_d, apical_hypha_d, lateral_hypha_1_d, 
            lateral_hypha_n_d);

    try
    {
        generate_structure(alternaria, N, centroid_x, centroid_y, centroid_z);
    }
    catch (Arg_error e)
    {
        delete alternaria;
        throw e;
    }

    return alternaria;
}


void Alternaria_density::generate_structure
(
    Alternaria* alternaria, 
    size_t      N, 
    float       centroid_x, 
    float       centroid_y, 
    float       centroid_z
) 
const throw (Arg_error)
{
    if (N == 0)
        return;

    const Spore_density* spore_d = alternaria->get_spore_d();
    const Apical_hypha_density* apical_hypha_d 
        = alternaria->get_apical_hypha_d();
    const Lateral_hypha_density* lateral_hypha_1_d 
        = alternaria->get_lateral_hypha_1_d();
    const Lateral_hypha_density* lateral_hypha_n_d 
        = alternaria->get_lateral_hypha_n_d();

    Apical_hypha* apical = apical_hypha_d->sample(centroid_x, centroid_y, 
            centroid_z, alternaria->get_theta(), alternaria->get_psi(),
            alternaria->get_base_level());
    alternaria->add_apical_hypha(apical);

    Structure* s;

    for (size_t n = 1; n < N; n++)
    {
        size_t level = alternaria->get_uniform_random_level();

        double u = sample_uniform_pdf_d(0, 1);

        if (level == 0)
        {
            //double lat_prob = lateral_hypha_1_p / (lat_prob + apical_hypha_p);
            double lat_prob = lateral_hypha_1_p*geometric_pmf_d(num_levels_p,level);

            if (u <= lat_prob)
            {
                s = alternaria->get_random_structure(level);
                alternaria->add_lateral_hypha(lateral_hypha_1_d->sample(s));
            }
            else
            {
                s = alternaria->get_random_terminal(level);
                alternaria->add_apical_hypha(apical_hypha_d->sample(s));
            }
        }
        else
        {
            //double lat_prob = lateral_hypha_n_p*geometric_pmf_d(num_levels_p,level-1);
            double lat_prob = lateral_hypha_n_p*geometric_pmf_d(num_levels_p,level);
            double spore_prob = spore_p;

            lat_prob   /= (lat_prob + spore_p + apical_hypha_p);
            spore_prob /= (lat_prob + spore_p + apical_hypha_p);

            if (u <= lat_prob)
            {
                s = alternaria->get_random_structure(level);
                alternaria->add_lateral_hypha(lateral_hypha_n_d->sample(s));
            }
            else if (u <= spore_prob)
            {
                s = alternaria->get_random_terminal(level);
                alternaria->add_spore(spore_d->sample(s));
            }
            else
            {
                s = alternaria->get_random_terminal(level);
                alternaria->add_apical_hypha(apical_hypha_d->sample(s));
            }
        }
    }
}

/* -------------------------------------------------------------------------- */











/* -----------------------------  Alternaria  ------------------------------- */

Alternaria::Alternaria
(
    float                        c_1,
    float                        theta,
    float                        psi,
    size_t                       base_level,
    const Alternaria_density*    alt_d,
    const Spore_density*         spore_d,
    const Apical_hypha_density*  apical_hypha_d,
    const Lateral_hypha_density* lateral_hypha_1_d,
    const Lateral_hypha_density* lateral_hypha_n_d
)
throw (Arg_error)
{
    this->alt_d             = alt_d;
    this->spore_d           = spore_d;
    this->apical_hypha_d    = apical_hypha_d;
    this->lateral_hypha_1_d = lateral_hypha_1_d;
    this->lateral_hypha_n_d = lateral_hypha_n_d;
    this->root              = 0;
    this->log_prob          = 0;
    this->c_1               = c_1;
    this->theta             = theta;
    this->psi               = psi;
    this->base_level        = base_level;

    check_c_1();

    level_structs.reserve(10);
    level_apical_hypha.reserve(10);
    level_lateral_hypha.reserve(10);
    level_spores.reserve(10);
    level_terminals.reserve(10);
}


Alternaria::Alternaria(const Alternaria& alt) 
{
    alt_d             = alt.alt_d;
    spore_d           = alt.spore_d;
    apical_hypha_d    = alt.apical_hypha_d;
    lateral_hypha_1_d = alt.lateral_hypha_1_d;
    lateral_hypha_n_d = alt.lateral_hypha_n_d;
    log_prob          = alt.log_prob;
    c_1               = alt.c_1;
    theta             = alt.theta;
    psi               = alt.psi;
    base_level        = alt.base_level;

    level_structs.reserve(10);
    level_apical_hypha.reserve(10);
    level_lateral_hypha.reserve(10);
    level_spores.reserve(10);
    level_terminals.reserve(10);

    root = 0;
    if (alt.root)
    {
        root = alt.root->recursively_clone();
        assert(base_level == root->get_level());
        build_levels(root);
    }
}


Alternaria::Alternaria
(
    const char*                  fname,
    const Alternaria_density*    alt_d,
    const Spore_density*         spore_d,
    const Apical_hypha_density*  apical_hypha_d,
    const Lateral_hypha_density* lateral_hypha_1_d,
    const Lateral_hypha_density* lateral_hypha_n_d
)
throw (IO_error, Arg_error)
{
    istringstream ist;
    ostringstream ost;
    ifstream in;

    this->alt_d             = alt_d;
    this->spore_d           = spore_d;
    this->apical_hypha_d    = apical_hypha_d;
    this->lateral_hypha_1_d = lateral_hypha_1_d;
    this->lateral_hypha_n_d = lateral_hypha_n_d;
    this->root              = 0;

    in.open(fname);
    if (in.fail())
    {
        ost << fname << ": Could not open file";
        throw IO_error(ost.str());
    }

    const char* member_value;

    if (!(member_value = get_member_value("C_1", in)))
    {
        ost << fname << ": Missing C_1";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> c_1;
    if (ist.fail())
    {
        ost << fname << ": Invalid C_1";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    check_c_1();

    if (!(member_value = get_member_value("Theta", in)))
    {
        ost << fname << ": Missing theta";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> theta;
    if (ist.fail())
    {
        ost << fname << ": Invalid theta";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    if (!(member_value = get_member_value("Psi", in)))
    {
        ost << fname << ": Missing psi";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> psi;
    if (ist.fail())
    {
        ost << fname << ": Invalid psi";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    if (!(member_value = get_member_value("Level", in)))
    {
        ost << fname << ": Missing base level";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> base_level;
    if (ist.fail())
    {
        ost << fname << ": Invalid base level";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    if (!(member_value = get_member_value("Log Prob", in)))
    {
        ost << fname << ": Missing log probability";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> log_prob;
    if (ist.fail())
    {
        ost << fname << ": Invalid log probability";
        throw Arg_error(ost.str());
    }

    try
    {
        read(in);
    }
    catch (Arg_error e)
    {
        ost << fname << ": " << e.get_msg();
        throw Arg_error(ost.str());
    }
    catch (IO_error e)
    {
        ost << fname << ": " << e.get_msg();
        throw IO_error(ost.str());
    }

    in.close();
    if (in.fail())
    {
        ost << fname << ": Could not close file";
        throw IO_error(ost.str());
    }

    level_structs.reserve(10);
    level_apical_hypha.reserve(10);
    level_lateral_hypha.reserve(10);
    level_spores.reserve(10);
    level_terminals.reserve(10);
}


Alternaria::~Alternaria()
{
    delete root;
    for (size_t i = 0; i < level_structs.size(); i++)
    {
        delete level_structs[ i ];
        delete level_apical_hypha[ i ];
        delete level_lateral_hypha[ i ];
        delete level_spores[ i ];
        delete level_terminals[ i ];
    }
}


Alternaria* Alternaria::clone() const
{
    return new Alternaria(*this);
}


Alternaria& Alternaria::operator= (const Alternaria& alternaria)
{
    alt_d             = alternaria.alt_d;
    spore_d           = alternaria.spore_d;
    apical_hypha_d    = alternaria.apical_hypha_d;
    lateral_hypha_1_d = alternaria.lateral_hypha_1_d;
    lateral_hypha_n_d = alternaria.lateral_hypha_n_d;
    log_prob          = alternaria.log_prob;
    c_1               = alternaria.c_1;
    theta             = alternaria.theta;
    psi               = alternaria.psi;
    base_level        = alternaria.base_level;

    for (size_t i = 0; i < level_structs.size(); i++)
    {
        delete level_structs[ i ];
        delete level_apical_hypha[ i ];
        delete level_lateral_hypha[ i ];
        delete level_spores[ i ];
        delete level_terminals[ i ];
    }

    level_structs.clear();
    level_apical_hypha.clear();
    level_lateral_hypha.clear();
    level_spores.clear();
    level_terminals.clear();

    level_structs.reserve(10);
    level_apical_hypha.reserve(10);
    level_lateral_hypha.reserve(10);
    level_spores.reserve(10);
    level_terminals.reserve(10);

    delete root;
    root = 0;
    if (alternaria.root)
    {
        root = alternaria.root->recursively_clone();
        assert(base_level == root->get_level());
        build_levels(root);
    }

    return *this;
}


double Alternaria::get_log_prob() const
{
#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_prob) && !isnan(log_prob));
#endif

    assert(log_prob <= 0);

    if (root)
    {
        assert(root->get_log_prob_descend() <= 0);
        return log_prob + root->get_log_prob_descend();
    }

    return log_prob;
}


void Alternaria::set_root(Structure* root) throw (Arg_error)
{ 
    assert(this->root==0); 
    assert(base_level == root->get_level());
    root->set_base_orientation(theta, psi);
    this->root = root;
}


void Alternaria::set_c_1(float c_1) throw (Arg_error)
{
    float old_c_1 = this->c_1;

    try
    {
        this->c_1 = c_1;
        check_c_1();
    }
    catch (Arg_error e)
    {
        set_c_1(old_c_1);
        throw e;
    }
}


void Alternaria::print(const char* fname) const throw (IO_error)
{
    ofstream out;

    out.open(fname);
    if (out.fail())
    {
        ostringstream ost;
        ost << fname << ": Could not open file";
        throw IO_error(ost.str());
    }

    print(out);

    out.close();
    if (out.fail())
    {
        ostringstream ost;
        ost << fname << ": Could not close file";
        throw IO_error(ost.str());
    }
}


void Alternaria::print(std::ostream& out) const
{
    out << "        C_1: " << c_1 << '\n'
        << "      Theta: " << theta << '\n'
        << "        Psi: " << psi << '\n'
        << "      Level: " << base_level << '\n'
        << "   Log Prob: " << log_prob << '\n';

    if (root)
    {
        root->recursively_print(out);
    }
}


size_t Alternaria::get_num_structures() const
{
    size_t n = 0;

    for (size_t i = 0; i < level_structs.size(); i++)
    {
        n += level_structs.at(i)->size();
    }
    
    return n;
}


size_t Alternaria::get_num_apical_hypha() const
{
    size_t n = 0;

    for (size_t i = 0; i < level_apical_hypha.size(); i++)
    {
        n += level_apical_hypha.at(i)->size();
    }
    
    return n;
}


size_t Alternaria::get_num_lateral_hypha() const
{
    size_t n = 0;

    for (size_t i = 0; i < level_lateral_hypha.size(); i++)
    {
        n += level_lateral_hypha.at(i)->size();
    }
    
    return n;
}


size_t Alternaria::get_num_spores() const
{
    size_t n = 0;

    for (size_t i = 0; i < level_spores.size(); i++)
    {
        n += level_spores.at(i)->size();
    }
    
    return n;
}


size_t Alternaria::get_uniform_random_level() const throw (Result_error)
{
    if (level_structs.size() == 0)
    {
        throw Result_error("No levels in Alternaria");
    }

    size_t level = static_cast<size_t>(floor(sample_uniform_pdf_d(
                    base_level, 0.9999*level_structs.size())));

    return level;
}


size_t Alternaria::get_geometric_random_level() const throw (Result_error)
{
    if (level_structs.size() == 0)
    {
        throw Result_error("No levels in Alternaria");
    }

    size_t level = level_structs.size();
    while (level >= level_structs.size())
    {
        level = sample_geometric_pmf_d(alt_d->get_num_levels_p());
    }

    return level;
}


Structure* Alternaria::get_random_structure(size_t level) const 
    throw (Arg_error, Result_error)
{
    try
    {
        if (level_structs.at(level)->size() == 0)
        {
            ostringstream ost;
            ost << "No structures at level " << level;
            throw Result_error(ost.str());
        }
    }
    catch (std::out_of_range e)
    {
        ostringstream ost;
        ost << "No level " << level << " in Alternaria";
        throw Arg_error(ost.str());
    }

    size_t i = static_cast<size_t>(floor(sample_uniform_pdf_d(0, 
                    0.9999 * level_structs.at(level)->size())));
    return level_structs.at(level)->at(i);
}


Apical_hypha* Alternaria::get_random_apical_hypha(size_t level) const
    throw (Arg_error, Result_error)
{
    try
    {
        if (level_apical_hypha.at(level)->size() == 0)
        {
            ostringstream ost;
            ost << "No apical hypha at level " << level;
            throw Result_error(ost.str());
        }
    }
    catch (std::out_of_range e)
    {
        ostringstream ost;
        ost << "No level " << level << " in Alternaria";
        throw Arg_error(ost.str());
    }

    size_t i = static_cast<size_t>(floor(sample_uniform_pdf_d(0, 
                    0.9999 * level_apical_hypha.at(level)->size())));
    return level_apical_hypha.at(level)->at(i);
}


Lateral_hypha* Alternaria::get_random_lateral_hypha(size_t level) const
    throw (Arg_error, Result_error)
{
    try
    {
        if (level_lateral_hypha.at(level)->size() == 0)
        {
            ostringstream ost;
            ost << "No lateral hypha at level " << level;
            throw Result_error(ost.str());
        }
    }
    catch (std::out_of_range e)
    {
        ostringstream ost;
        ost << "No level " << level << " in Alternaria";
        throw Arg_error(ost.str());
    }

    size_t i = static_cast<size_t>(floor(sample_uniform_pdf_d(0, 
                    0.9999 * level_lateral_hypha.at(level)->size())));
    return level_lateral_hypha.at(level)->at(i);
}


Spore* Alternaria::get_random_spore(size_t level) const
    throw (Arg_error, Result_error)
{
    try
    {
        if (level_spores.at(level)->size() == 0)
        {
            ostringstream ost;
            ost << "No spores at level " << level;
            throw Result_error(ost.str());
        }
    }
    catch (std::out_of_range e)
    {
        ostringstream ost;
        ost << "No level " << level << " in Alternaria";
        throw Arg_error(ost.str());
    }

    size_t i = static_cast<size_t>(floor(sample_uniform_pdf_d(0, 
                    0.9999 * level_spores.at(level)->size())));
    return level_spores.at(level)->at(i);
}


Structure* Alternaria::get_random_terminal(size_t level) const
    throw (Arg_error, Result_error)
{
    try
    {
        if (level_terminals.at(level)->size() == 0)
        {
            ostringstream ost;
            ost << "No terminals at level " << level;
            throw Result_error(ost.str());
        }
    }
    catch (std::out_of_range e)
    {
        ostringstream ost;
        ost << "No level " << level << " in Alternaria";
        throw Arg_error(ost.str());
    }

    size_t i = static_cast<size_t>(floor(sample_uniform_pdf_d(0, 
                    0.9999 * level_terminals.at(level)->size())));
    return level_terminals.at(level)->at(i);
}


Apical_hypha* Alternaria::get_random_terminal_apical_hypha(size_t level) const
    throw (Arg_error, Result_error)
{
    try
    {
        if (level_apical_hypha.at(level)->size() == 0)
        {
            ostringstream ost;
            ost << "No apical hypha at level " << level;
            throw Result_error(ost.str());
        }
        size_t n = 0;
        for (size_t i = 0; i < level_apical_hypha.at(level)->size(); i++)
        {
            if (level_apical_hypha.at(level)->at(i)->get_apical() == 0)
            {
                n++;
            }
        }
        if (n == 0)
        {
            ostringstream ost;
            ost << "No terminal apical hypha at level " << level;
            throw Result_error(ost.str());
        }
    }
    catch (std::out_of_range e)
    {
        ostringstream ost;
        ost << "No level " << level << " in Alternaria";
        throw Arg_error(ost.str());
    }

    Apical_hypha* hypha = 0;
    while (!hypha)
    {
        size_t i = static_cast<size_t>(floor(sample_uniform_pdf_d(0, 
                        0.9999 * level_apical_hypha.at(level)->size())));
        if (level_apical_hypha.at(level)->at(i)->get_apical() == 0)
        {
            hypha = level_apical_hypha.at(level)->at(i);
        }
    }

    return hypha;
}


Spore* Alternaria::get_random_terminal_spore(size_t level) const
    throw (Arg_error, Result_error)
{
    try
    {
        if (level_spores.at(level)->size() == 0)
        {
            ostringstream ost;
            ost << "No spores at level " << level;
            throw Result_error(ost.str());
        }
        size_t n = 0;
        for (size_t i = 0; i < level_spores.at(level)->size(); i++)
        {
            if (level_spores.at(level)->at(i)->get_apical() == 0)
            {
                n++;
            }
        }
        if (n == 0)
        {
            ostringstream ost;
            ost << "No terminal spores at level " << level;
            throw Result_error(ost.str());
        }
    }
    catch (std::out_of_range e)
    {
        ostringstream ost;
        ost << "No level " << level << " in Alternaria";
        throw Arg_error(ost.str());
    }

    Spore* spore = 0;
    while (!spore)
    {
        size_t i = static_cast<size_t>(floor(sample_uniform_pdf_d(0, 
                        0.9999 * level_spores.at(level)->size())));
        if (level_spores.at(level)->at(i)->get_apical() == 0)
        {
            spore = level_spores.at(level)->at(i);
        }
    }

    return spore;
}


Lateral_hypha* Alternaria::get_random_terminal_lateral_hypha(size_t level) const
    throw (Arg_error, Result_error)
{
    try
    {
        if (level_lateral_hypha.at(level)->size() == 0)
        {
            ostringstream ost;
            ost << "No lateral hypha at level " << level;
            throw Result_error(ost.str());
        }
        size_t n = 0;
        for (size_t i = 0; i < level_lateral_hypha.at(level)->size(); i++)
        {
            if (level_lateral_hypha.at(level)->at(i)->get_apical() == 0)
            {
                n++;
            }
        }
        if (n == 0)
        {
            ostringstream ost;
            ost << "No terminal apical hypha at level " << level;
            throw Result_error(ost.str());
        }
    }
    catch (std::out_of_range e)
    {
        ostringstream ost;
        ost << "No level " << level << " in Alternaria";
        throw Arg_error(ost.str());
    }

    Lateral_hypha* hypha = 0;
    while (!hypha)
    {
        size_t i = static_cast<size_t>(floor(sample_uniform_pdf_d(0, 
                        0.9999 * level_lateral_hypha.at(level)->size())));
        if (level_lateral_hypha.at(level)->at(i)->get_apical() == 0)
        {
            hypha = level_lateral_hypha.at(level)->at(i);
        }
    }

    return hypha;
}


void Alternaria::add_lateral_hypha(Lateral_hypha* h) throw (Arg_error)
{
    size_t level = h->get_level();

    ensure_level_space(level);

    if (!h->get_parent())
    {
        set_root(h);
    }

    level_terminals[ level ]->push_back(h);
    level_structs[ level ]->push_back(h);
    level_lateral_hypha[ level ]->push_back(h);

    double p;
    double lateral_hypha_1_p = alt_d->get_lateral_hypha_1_p();
    double lateral_hypha_n_p = alt_d->get_lateral_hypha_n_p();
    double spore_p = alt_d->get_spore_p();
    double apical_hypha_p = alt_d->get_apical_hypha_p();
    double num_levels_p = alt_d->get_num_levels_p();

    // Note that the level must be > 0.
    if (level == 1)
    {
        //p = lateral_hypha_1_p / (lateral_hypha_1_p + apical_hypha_p);
        p = lateral_hypha_1_p*geometric_pmf_d(num_levels_p, level-1);
        p /= p + apical_hypha_p;
    }
    else
    {
        //p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, level-2);
        p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, level-1);
        p /= p + spore_p + apical_hypha_p;
    }

    log_prob += (p < JWSC_MIN_LOG_ARG) ? log(JWSC_MIN_LOG_ARG) : log(p);

    check_levels();
}


void Alternaria::remove_lateral_hypha(Lateral_hypha* h) throw (Arg_error)
{
    size_t level = h->get_level();

    ensure_level_space(level);

    if (find(level_terminals[ level ]->begin(), 
             level_terminals[ level ]->end(), h) 
        == level_terminals[ level ]->end())
    {
        throw Arg_error("Cannot remove a non-terminal lateral hypha");
    }

    Structure* parent = h->get_parent();

    if (parent)
    {
        parent->remove_lateral(h);
    }
    else
    {
        root = 0;
    }

    level_structs[ level ]->erase(
            find(level_structs[ level ]->begin(),
                 level_structs[ level ]->end(), (Structure*)h));
    level_terminals[ level ]->erase(
            find(level_terminals[ level ]->begin(),
                 level_terminals[ level ]->end(), (Structure*)h));
    level_lateral_hypha[ level ]->erase(
            find(level_lateral_hypha[ level ]->begin(),
                 level_lateral_hypha[ level ]->end(), h));

    double p;
    double lateral_hypha_1_p = alt_d->get_lateral_hypha_1_p();
    double lateral_hypha_n_p = alt_d->get_lateral_hypha_n_p();
    double spore_p = alt_d->get_spore_p();
    double apical_hypha_p = alt_d->get_apical_hypha_p();
    double num_levels_p = alt_d->get_num_levels_p();

    // Note that the level must be > 0.
    if (level == 1)
    {
        //p = lateral_hypha_1_p / (lateral_hypha_1_p + apical_hypha_p);
        p = lateral_hypha_1_p*geometric_pmf_d(num_levels_p, level-1);
        p /= p + apical_hypha_p;
    }
    else
    {
        //p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, level-2);
        p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, level-1);
        p /= p + spore_p + apical_hypha_p;
    }

    log_prob -= (p < JWSC_MIN_LOG_ARG) ? log(JWSC_MIN_LOG_ARG) : log(p);

    delete h;

    check_levels();
}


void Alternaria::add_apical_hypha(Apical_hypha* h) throw (Arg_error)
{
    size_t level = h->get_level();

    ensure_level_space(level);

    if (h->get_parent())
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                h->get_parent(), (Structure*)h);
    }
    else
    {
        set_root(h);
        level_terminals[ level ]->push_back(h);
    }

    level_structs[ level ]->push_back(h);
    level_apical_hypha[ level ]->push_back(h);

    log_prob += log(alt_d->get_apical_hypha_p());

    check_levels();
}


void Alternaria::remove_apical_hypha(Apical_hypha* h) throw (Arg_error)
{
    size_t level = h->get_level();

    ensure_level_space(level);

    if (find(level_terminals[ level ]->begin(), 
             level_terminals[ level ]->end(), h) 
        == level_terminals[ level ]->end())
    {
        throw Arg_error("Cannot remove a non-terminal apical hypha");
    }

    Structure* parent = h->get_parent();

    if (parent)
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                (Structure*)h, parent);
        parent->set_apical(0);
    }
    else
    {
        level_terminals[ level ]->erase(
                find(level_terminals[ level ]->begin(),
                     level_terminals[ level ]->end(), (Structure*)h));
        root = 0;
    }

    level_structs[ level ]->erase(
            find(level_structs[ level ]->begin(),
                 level_structs[ level ]->end(), (Structure*)h));
    level_apical_hypha[ level ]->erase(
            find(level_apical_hypha[ level ]->begin(),
                 level_apical_hypha[ level ]->end(), h));

    log_prob -= log(alt_d->get_apical_hypha_p());

    delete h;

    check_levels();
}


void Alternaria::merge_2_apical_with_spore
(
    Apical_hypha* hypha_1,
    const Spore*  rvals
) 
throw (Arg_error)
{
    size_t level = hypha_1->get_level();

    ensure_level_space(level);

    if (find(level_structs[ level ]->begin(), 
             level_structs[ level ]->end(), hypha_1) 
        == level_structs[ level ]->end())
    {
        throw Arg_error("Cannot find apical hypha to merge");
    }

    Spore* merged = 0;

    Apical_hypha* hypha_2 = hypha_1->merge_with_apical(&merged, rvals);

    if (!(merged->get_apical()))
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                (Structure*)merged, (Structure*)hypha_2);
    }
    level_structs[ level ]->erase(
            find(level_structs[ level ]->begin(),
                 level_structs[ level ]->end(), (Structure*)hypha_1));
    level_structs[ level ]->erase(
            find(level_structs[ level ]->begin(),
                 level_structs[ level ]->end(), (Structure*)hypha_2));
    level_apical_hypha[ level ]->erase(
            find(level_apical_hypha[ level ]->begin(),
                 level_apical_hypha[ level ]->end(), (Structure*)hypha_1));
    level_apical_hypha[ level ]->erase(
            find(level_apical_hypha[ level ]->begin(),
                 level_apical_hypha[ level ]->end(), (Structure*)hypha_2));

    log_prob -= log(alt_d->get_apical_hypha_p());
    log_prob -= log(alt_d->get_apical_hypha_p());

    delete hypha_1;
    delete hypha_2;

    if (!(merged->get_parent()))
    {
        root = 0;
        set_root(merged);
    }

    level_structs[ level ]->push_back(merged);
    level_spores[ level ]->push_back(merged);

    log_prob += log(alt_d->get_spore_p());

    check_levels();
}


void Alternaria::split_apical_hypha_into_apical
(
    Apical_hypha*       h,
    const Apical_hypha* rvals_1,
    const Apical_hypha* rvals_2
) 
throw (Arg_error)
{
    size_t level = h->get_level();

    ensure_level_space(level);

    if (find(level_structs[ level ]->begin(), 
             level_structs[ level ]->end(), h) 
        == level_structs[ level ]->end())
    {
        throw Arg_error("Cannot find apical hypha to split");
    }

    Apical_hypha* split = h->split_into_apical(rvals_1, rvals_2);
    if (!(split->get_apical()))
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                (Structure*)h, (Structure*)split);
    }

    level_structs[ level ]->push_back(split);
    level_apical_hypha[ level ]->push_back(split);

    log_prob += log(alt_d->get_apical_hypha_p());

    check_levels();
}


void Alternaria::merge_apical_hypha_with_apical
(
    Apical_hypha*       h, 
    const Apical_hypha* rvals
) 
throw (Arg_error)
{
    size_t level = h->get_level();

    ensure_level_space(level);

    if (find(level_structs[ level ]->begin(), 
             level_structs[ level ]->end(), h) 
        == level_structs[ level ]->end())
    {
        throw Arg_error("Cannot find apical hypha to merge");
    }

    Apical_hypha* merged = h->merge_with_apical(rvals);

    if (!(h->get_apical()))
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                (Structure*)merged, (Structure*)h);
    }
    level_structs[ level ]->erase(
            find(level_structs[ level ]->begin(),
                 level_structs[ level ]->end(), (Structure*)merged));
    level_apical_hypha[ level ]->erase(
            find(level_apical_hypha[ level ]->begin(),
                 level_apical_hypha[ level ]->end(), (Structure*)merged));

    log_prob -= log(alt_d->get_apical_hypha_p());

    delete merged;

    check_levels();
}


void Alternaria::split_apical_hypha_into_lateral
(
    Apical_hypha*        h,
    const Apical_hypha*  rvals_1,
    const Lateral_hypha* rvals_2
) 
throw (Arg_error)
{
    size_t level = h->get_level();

    ensure_level_space(level+1);

    if (find(level_structs[ level ]->begin(), 
             level_structs[ level ]->end(), h) 
        == level_structs[ level ]->end())
    {
        throw Arg_error("Cannot find apical hypha to split");
    }

    Lateral_hypha* split = h->split_into_lateral(rvals_1, rvals_2);
    if (split->get_apical())
    {
        level_terminals[ level ]->push_back(h);
        recursively_update_level(split->get_apical(), level);
    }
    else
    {
        level_terminals[ level+1 ]->push_back(split);
    }

    for (list<Lateral_structure*>::iterator it = split->get_laterals().begin();
            it != split->get_laterals().end(); it++)
    {
        recursively_update_level(*it, level+1);
    }

    level_structs[ level+1 ]->push_back(split);
    level_lateral_hypha[ level+1 ]->push_back(split);

    double p;
    double lateral_hypha_1_p = alt_d->get_lateral_hypha_1_p();
    double lateral_hypha_n_p = alt_d->get_lateral_hypha_n_p();
    double spore_p = alt_d->get_spore_p();
    double apical_hypha_p = alt_d->get_apical_hypha_p();
    double num_levels_p = alt_d->get_num_levels_p();

    // Note that the level must be >= 0.
    if (level == 0)
    {
        //p = lateral_hypha_1_p / (lateral_hypha_1_p + apical_hypha_p);
        p = lateral_hypha_1_p*geometric_pmf_d(num_levels_p, level);
        p /= p + apical_hypha_p;
    }
    else
    {
        //p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, level-1);
        p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, level);
        p /= p + spore_p + apical_hypha_p;
    }

    log_prob += (p < JWSC_MIN_LOG_ARG) ? log(JWSC_MIN_LOG_ARG) : log(p);

    check_levels();
}


void Alternaria::merge_apical_hypha_with_lateral
(
    Apical_hypha*       h, 
    Lateral_hypha*      lat,
    const Apical_hypha* rvals
) 
throw (Arg_error)
{
    size_t level = h->get_level();

    ensure_level_space(level+1);

    if (find(level_terminals[ level ]->begin(), 
             level_terminals[ level ]->end(), h) 
        == level_terminals[ level ]->end())
    {
        throw Arg_error("Cannot find terminal apical hypha to merge");
    }

    Lateral_hypha* merged = h->merge_with_lateral(lat, rvals);

    if (h->get_apical())
    {
        recursively_update_level(h->get_apical(), level+1);

        level_terminals[ level ]->erase(
                find(level_terminals[ level ]->begin(),
                     level_terminals[ level ]->end(), (Structure*)h));
    }
    else
    {
        level_terminals[ level+1 ]->erase(
                find(level_terminals[ level+1 ]->begin(),
                     level_terminals[ level+1 ]->end(), (Structure*)merged));
    }

    for (list<Lateral_structure*>::iterator it = h->get_laterals().begin();
            it != h->get_laterals().end(); it++)
    {
        try
        {
            recursively_update_level(*it, level+2);
        }
        catch (Arg_error e)
        {
            // Ignore, didn't need updating.
        }
    }

    level_structs[ level+1 ]->erase(
            find(level_structs[ level+1 ]->begin(),
                 level_structs[ level+1 ]->end(), (Structure*)merged));
    level_lateral_hypha[ level+1 ]->erase(
            find(level_lateral_hypha[ level+1 ]->begin(),
                 level_lateral_hypha[ level+1 ]->end(), (Structure*)merged));

    double p;
    double lateral_hypha_1_p = alt_d->get_lateral_hypha_1_p();
    double lateral_hypha_n_p = alt_d->get_lateral_hypha_n_p();
    double spore_p = alt_d->get_spore_p();
    double apical_hypha_p = alt_d->get_apical_hypha_p();
    double num_levels_p = alt_d->get_num_levels_p();

    // Note that the level must be >= 0.
    if (level == 0)
    {
        //p = lateral_hypha_1_p / (lateral_hypha_1_p + apical_hypha_p);
        p = lateral_hypha_1_p*geometric_pmf_d(num_levels_p, level);
        p /= p + apical_hypha_p;
    }
    else
    {
        //p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, level-1);
        p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, level);
        p /= p + spore_p + apical_hypha_p;
    }

    log_prob -= (p < JWSC_MIN_LOG_ARG) ? log(JWSC_MIN_LOG_ARG) : log(p);

    delete merged;

    check_levels();
}


void Alternaria::replace_apical_hypha_with_spore
(
    Apical_hypha* hypha,
    Spore*        spore
) 
throw (Arg_error)
{
    size_t level = hypha->get_level();

    ensure_level_space(level);

    if (find(level_structs[ level ]->begin(), 
             level_structs[ level ]->end(), hypha) 
        == level_structs[ level ]->end())
    {
        throw Arg_error("Cannot find apical hypha to replace");
    }

    hypha->replace(spore);

    if (!(spore->get_apical()))
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                (Structure*)hypha, (Structure*)spore);
    }
    level_structs[ level ]->erase(
            find(level_structs[ level ]->begin(),
                 level_structs[ level ]->end(), (Structure*)hypha));
    level_apical_hypha[ level ]->erase(
            find(level_apical_hypha[ level ]->begin(),
                 level_apical_hypha[ level ]->end(), (Structure*)hypha));

    log_prob -= log(alt_d->get_apical_hypha_p());

    delete hypha;

    if (!(spore->get_parent()))
    {
        root = 0;
        set_root(spore);
    }

    level_structs[ level ]->push_back(spore);
    level_spores[ level ]->push_back(spore);

    log_prob += log(alt_d->get_spore_p());

    check_levels();
}


void Alternaria::add_spore(Spore* s) throw (Arg_error)
{
    size_t level = s->get_level();

    ensure_level_space(level);

    if (s->get_parent())
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                s->get_parent(), (Structure*)s);
    }
    else
    {
        set_root(s);
        level_terminals[ level ]->push_back(s);
    }

    level_structs[ level ]->push_back(s);
    level_spores[ level ]->push_back(s);

    log_prob += log(alt_d->get_spore_p());

    check_levels();
}


void Alternaria::remove_spore(Spore* s) throw (Arg_error)
{
    size_t level = s->get_level();

    ensure_level_space(level);

    if (find(level_terminals[ level ]->begin(), 
             level_terminals[ level ]->end(), s) 
        == level_terminals[ level ]->end())
    {
        throw Arg_error("Cannot remove a non-terminal spore");
    }

    Structure* parent = s->get_parent();

    if (parent)
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                (Structure*)s, parent);
        parent->set_apical(0);
    }
    else
    {
        level_terminals[ level ]->erase(
                find(level_terminals[ level ]->begin(),
                     level_terminals[ level ]->end(), (Structure*)s));
        root = 0;
    }

    level_structs[ level ]->erase(
            find(level_structs[ level ]->begin(),
                 level_structs[ level ]->end(), (Structure*)s));
    level_spores[ level ]->erase(
            find(level_spores[ level ]->begin(),
                 level_spores[ level ]->end(), s));

    log_prob -= log(alt_d->get_spore_p());

    delete s;

    check_levels();
}


void Alternaria::replace_spore_with_apical_hypha
(
    Spore*        spore,
    Apical_hypha* hypha
) 
throw (Arg_error)
{
    size_t level = spore->get_level();

    ensure_level_space(level);

    if (find(level_structs[ level ]->begin(), 
             level_structs[ level ]->end(), spore) 
        == level_structs[ level ]->end())
    {
        throw Arg_error("Cannot find spore to replace");
    }

    spore->replace(hypha);

    if (!(hypha->get_apical()))
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                (Structure*)spore, (Structure*)hypha);
    }
    level_structs[ level ]->erase(
            find(level_structs[ level ]->begin(),
                 level_structs[ level ]->end(), (Structure*)spore));
    level_spores[ level ]->erase(
            find(level_spores[ level ]->begin(),
                 level_spores[ level ]->end(), (Structure*)spore));

    log_prob -= log(alt_d->get_spore_p());

    delete spore;

    if (!(hypha->get_parent()))
    {
        root = 0;
        set_root(hypha);
    }

    level_structs[ level ]->push_back(hypha);
    level_apical_hypha[ level ]->push_back(hypha);

    log_prob += log(alt_d->get_apical_hypha_p());

    check_levels();
}


void Alternaria::split_spore_into_2_apical
(
    Spore*              spore,
    const Apical_hypha* rvals_1,
    const Apical_hypha* rvals_2
) 
throw (Arg_error)
{
    size_t level = spore->get_level();

    ensure_level_space(level);

    if (find(level_structs[ level ]->begin(), 
             level_structs[ level ]->end(), spore) 
        == level_structs[ level ]->end())
    {
        throw Arg_error("Cannot find spore to split");
    }

    Apical_hypha* split_1 = spore->split_into_apical(rvals_1, rvals_2);
    Apical_hypha* split_2 = dynamic_cast<Apical_hypha*>(split_1->get_apical());
    assert(split_2);

    if (!(split_2->get_apical()))
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                (Structure*)spore, (Structure*)split_2);
    }
    level_structs[ level ]->erase(
            find(level_structs[ level ]->begin(),
                 level_structs[ level ]->end(), (Structure*)spore));
    level_spores[ level ]->erase(
            find(level_spores[ level ]->begin(),
                 level_spores[ level ]->end(), (Structure*)spore));

    log_prob -= log(alt_d->get_spore_p());

    delete spore;

    if (!(split_1->get_parent()))
    {
        root = 0;
        set_root(split_1);
    }

    level_structs[ level ]->push_back(split_1);
    level_structs[ level ]->push_back(split_2);
    level_apical_hypha[ level ]->push_back(split_1);
    level_apical_hypha[ level ]->push_back(split_2);

    log_prob += log(alt_d->get_apical_hypha_p());
    log_prob += log(alt_d->get_apical_hypha_p());

    check_levels();
}


void Alternaria::split_spore_into_spore
(
    Spore*       s,
    const Spore* rvals_1,
    const Spore* rvals_2
) 
throw (Arg_error)
{
    size_t level = s->get_level();

    ensure_level_space(level);

    if (find(level_structs[ level ]->begin(), 
             level_structs[ level ]->end(), s) 
        == level_structs[ level ]->end())
    {
        throw Arg_error("Cannot find spore to split");
    }

    Spore* split = s->split_into_apical(rvals_1, rvals_2);
    if (!(split->get_apical()))
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                (Structure*)s, (Structure*)split);
    }

    level_structs[ level ]->push_back(split);
    level_spores[ level ]->push_back(split);

    log_prob += log(alt_d->get_spore_p());

    check_levels();
}


void Alternaria::merge_spore_with_spore
(
    Spore*       s, 
    const Spore* rvals
) 
throw (Arg_error)
{
    size_t level = s->get_level();

    ensure_level_space(level);

    if (find(level_structs[ level ]->begin(), 
             level_structs[ level ]->end(), s) 
        == level_structs[ level ]->end())
    {
        throw Arg_error("Cannot find spore to merge");
    }

    Spore* merged = s->merge_with_apical(rvals);

    if (!(s->get_apical()))
    {
        replace(level_terminals[ level ]->begin(),
                level_terminals[ level ]->end(), 
                (Structure*)merged, (Structure*)s);
    }
    level_structs[ level ]->erase(
            find(level_structs[ level ]->begin(),
                 level_structs[ level ]->end(), (Structure*)merged));
    level_spores[ level ]->erase(
            find(level_spores[ level ]->begin(),
                 level_spores[ level ]->end(), (Structure*)merged));

    log_prob -= log(alt_d->get_spore_p());

    delete merged;

    check_levels();
}


void Alternaria::read(std::istream& in) throw (IO_error, Arg_error)
{
    list<Printed_structure*>         printed_structs;
    list<Printed_apical_structure*>  printed_apicals;
    list<Printed_lateral_structure*> printed_laterals;

    log_prob = 0;

    while (in.peek() != EOF)
    {
        Printed_apical_structure*  apical;
        Printed_lateral_structure* lateral;

        if (Printed_structure::is_on_istream(in, Apical_hypha::type_str()))
        {
            apical = new Printed_apical_hypha(in, apical_hypha_d);
            printed_structs.push_back(apical);
            printed_apicals.push_back(apical);

            log_prob += log(alt_d->get_apical_hypha_p());
        }
        else if (Printed_structure::is_on_istream(in, 
                    Lateral_hypha::type_str()))
        {
            lateral = new Printed_lateral_hypha(in, lateral_hypha_1_d, 
                    lateral_hypha_n_d);
            printed_structs.push_back(lateral);
            printed_laterals.push_back(lateral);

            double p;
            double lateral_hypha_1_p = alt_d->get_lateral_hypha_1_p();
            double lateral_hypha_n_p = alt_d->get_lateral_hypha_n_p();
            double spore_p = alt_d->get_spore_p();
            double apical_hypha_p = alt_d->get_apical_hypha_p();
            double num_levels_p = alt_d->get_num_levels_p();

            size_t level = lateral->get_level();

            // Note that the level must be > 0.
            if (level == 1)
            {
                //p = lateral_hypha_1_p / (lateral_hypha_1_p + apical_hypha_p);
                p = lateral_hypha_1_p*geometric_pmf_d(num_levels_p, level-1);
                p /= p + apical_hypha_p;
            }
            else
            {
                //p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, level-2);
                p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, level-1);
                p /= p + spore_p + apical_hypha_p;
            }

            log_prob += (p < JWSC_MIN_LOG_ARG) ? log(JWSC_MIN_LOG_ARG) : log(p);
        }
        else if (Printed_structure::is_on_istream(in, Spore::type_str()))
        {
            apical = new Printed_spore(in, spore_d);
            printed_structs.push_back(apical);
            printed_apicals.push_back(apical);

            log_prob += log(alt_d->get_spore_p());
        }
        else
        {
            throw Arg_error("Contains an unknown structure type");
        }
    }

    Printed_structure* printed_root = 0;
    list<Printed_structure*>::const_iterator it;

    for (it = printed_structs.begin(); it != printed_structs.end(); it++)
    {
        if ((*it)->get_parent() == 0)
        {
            printed_root = (*it);
            break;
        }
    }

    if (!printed_root)
    {
        throw Arg_error("No root in structures to read");
    }

    try
    {
        root = printed_root->recursively_convert(0, printed_apicals.begin(), 
                printed_apicals.end(), printed_laterals.begin(), 
                printed_laterals.end());
        if (root)
        {
            build_levels(root);
        }
    }
    catch (Arg_error e)
    {
        delete root;
        for (it = printed_structs.begin(); it != printed_structs.end(); it++)
        {
            delete *it;
        }
        throw e;
    }

    for (it = printed_structs.begin(); it != printed_structs.end(); it++)
    {
        delete *it;
    }
}


void Alternaria::build_levels(Structure* s)
{
    size_t level = s->get_level();

    ensure_level_space(level);

    level_structs[ level ]->push_back(s);

    if (dynamic_cast<Apical_hypha*>(s))
    {
        level_apical_hypha[ level ]->push_back((Apical_hypha*)s);
    }
    else if (dynamic_cast<Lateral_hypha*>(s))
    {
        level_lateral_hypha[ level ]->push_back((Lateral_hypha*)s);
    }
    else if (dynamic_cast<Spore*>(s))
    {
        level_spores[ level ]->push_back((Spore*)s);
    }

    Structure* apical = s->get_apical();
    if (apical)
    {
        assert(level == s->get_apical()->get_level());
        build_levels(s->get_apical());
    }
    else
    {
        level_terminals[ level ]->push_back(s);
    }

    list<Lateral_structure*> laterals = s->get_laterals();
    for (list<Lateral_structure*>::iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        assert((level+1) == (*it)->get_level());
        build_levels(*it);
    }
}


void Alternaria::check_levels()
{
    for (size_t i = 0; i < level_structs.size(); i++)
    {
        for (size_t j = 0; j < level_structs.at(i)->size(); j++)
        {
            assert(level_structs.at(i)->at(j)->get_level() == i);
        }
    }
    for (size_t i = 0; i < level_lateral_hypha.size(); i++)
    {
        for (size_t j = 0; j < level_lateral_hypha.at(i)->size(); j++)
        {
            assert(level_lateral_hypha.at(i)->at(j)->get_level() == i);
        }
    }
    for (size_t i = 0; i < level_apical_hypha.size(); i++)
    {
        for (size_t j = 0; j < level_apical_hypha.at(i)->size(); j++)
        {
            assert(level_apical_hypha.at(i)->at(j)->get_level() == i);
        }
    }
}


void Alternaria::check_c_1() const throw (Arg_error)
{
    if (c_1 < 0)
    {
        throw Arg_error("C_1 < 0");
    }
}


void Alternaria::ensure_level_space(size_t level)
{
    while (level_structs.size() <= level)
    {
        vector<Structure*>* structs = new vector<Structure*>();
        structs->reserve(200);
        level_structs.push_back(structs);

        vector<Apical_hypha*>* a_hypha = new vector<Apical_hypha*>();
        a_hypha->reserve(200);
        level_apical_hypha.push_back(a_hypha);

        vector<Lateral_hypha*>* l_hypha = new vector<Lateral_hypha*>();
        l_hypha->reserve(200);
        level_lateral_hypha.push_back(l_hypha);

        vector<Spore*>* spores = new vector<Spore*>();
        spores->reserve(200);
        level_spores.push_back(spores);

        structs = new vector<Structure*>();
        structs->reserve(200);
        level_terminals.push_back(structs);
    }
}


void Alternaria::recursively_update_level(Structure* s, size_t old_level)
    throw (Arg_error)
{
    size_t new_level = s->get_level();

    if (old_level == new_level)
        return;

    ensure_level_space((old_level > new_level) ? old_level : new_level);

    if (find(level_structs[ old_level ]->begin(), 
             level_structs[ old_level ]->end(), s) 
        == level_structs[ old_level ]->end())
    {
        ostringstream ost;
        ost << "Structure " << s->get_id() << " not in level " << old_level;
        throw Arg_error(ost.str());
    }

    level_structs[ old_level ]->erase(
            find(level_structs[ old_level ]->begin(),
                 level_structs[ old_level ]->end(), s));
    level_structs[ new_level ]->push_back(s);

    if (level_terminals[ old_level ]->size() > 0 &&
        find(level_terminals[ old_level ]->begin(), 
             level_terminals[ old_level ]->end(), s) 
        != level_terminals[ old_level ]->end())
    {
        level_terminals[ old_level ]->erase(
                find(level_terminals[ old_level ]->begin(),
                     level_terminals[ old_level ]->end(), s));
        level_terminals[ new_level ]->push_back(s);
    }

    if (dynamic_cast<Apical_hypha*>(s))
    {
        level_apical_hypha[ old_level ]->erase(
                find(level_apical_hypha[ old_level ]->begin(),
                     level_apical_hypha[ old_level ]->end(), (Apical_hypha*)s));
        level_apical_hypha[ new_level ]->push_back((Apical_hypha*)s);
    }
    else if (dynamic_cast<Lateral_hypha*>(s))
    {
        level_lateral_hypha[ old_level ]->erase(
                find(level_lateral_hypha[ old_level ]->begin(),
                   level_lateral_hypha[ old_level ]->end(), (Lateral_hypha*)s));
        level_lateral_hypha[ new_level ]->push_back((Lateral_hypha*)s);

        double p;
        double lateral_hypha_1_p = alt_d->get_lateral_hypha_1_p();
        double lateral_hypha_n_p = alt_d->get_lateral_hypha_n_p();
        double spore_p = alt_d->get_spore_p();
        double apical_hypha_p = alt_d->get_apical_hypha_p();
        double num_levels_p = alt_d->get_num_levels_p();

        // Note that the old_level must be > 0.
        if (old_level == 1)
        {
            //p = lateral_hypha_1_p / (lateral_hypha_1_p + apical_hypha_p);
            p = lateral_hypha_1_p*geometric_pmf_d(num_levels_p, old_level-1);
            p /= p + apical_hypha_p;
        }
        else
        {
            //p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, old_level-2);
            p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, old_level-1);
            p /= p + spore_p + apical_hypha_p;
        }
        log_prob -= (p < JWSC_MIN_LOG_ARG) ? log(JWSC_MIN_LOG_ARG) : log(p);

        // Note that the new_level must be > 0.
        if (new_level == 1)
        {
            //p = lateral_hypha_1_p / (lateral_hypha_1_p + apical_hypha_p);
            p = lateral_hypha_1_p*geometric_pmf_d(num_levels_p, new_level-1);
            p /= p + apical_hypha_p;
        }
        else
        {
            //p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, new_level-2);
            p  = lateral_hypha_n_p*geometric_pmf_d(num_levels_p, new_level-1);
            p /= p + spore_p + apical_hypha_p;
        }
        log_prob += (p < JWSC_MIN_LOG_ARG) ? log(JWSC_MIN_LOG_ARG) : log(p);
    }
    else if (dynamic_cast<Spore*>(s))
    {
        level_spores[ old_level ]->erase(
                find(level_spores[ old_level ]->begin(),
                     level_spores[ old_level ]->end(), (Spore*)s));
        level_spores[ new_level ]->push_back((Spore*)s);
    }
    else
    {
        ostringstream ost;
        ost << "Structure " << s->get_id() << " is unknown type";
        std::cerr << "BUG: Alternaria::recursively_update_level: " 
                  << ost.str() << '\n';
        assert(0);
    }

    if (s->get_apical())
    {
        recursively_update_level(s->get_apical(), old_level);
    }
    for (list<Lateral_structure*>::iterator it = s->get_laterals().begin();
            it != s->get_laterals().end(); it++)
    {
        recursively_update_level(*it, old_level+1);
    }
}

/* -------------------------------------------------------------------------- */











/* -----------------------  Alternaria_model_density  ----------------------- */

#if defined ALTERNARIA_HAVE_FFTW3F
Alternaria_model* Alternaria_model_density::sample
(
    size_t                       N,
    float                        centroid_x, 
    float                        centroid_y, 
    float                        centroid_z,
    float                        theta,
    float                        psi,
    size_t                       base_level,
    const Spore_density*         spore_d,
    const Apical_hypha_density*  apical_hypha_d,
    const Lateral_hypha_density* lateral_hypha_1_d,
    const Lateral_hypha_density* lateral_hypha_n_d,
    const PSF_model*             psf,
    const Imaging_model*         imaging
)
const throw (Arg_error)
{
    Alternaria_model* alternaria = new Alternaria_model(sample_c_1(), theta,
            psi, base_level, this, spore_d, apical_hypha_d, lateral_hypha_1_d, 
            lateral_hypha_n_d, psf, imaging);

    try
    {
        generate_structure(alternaria, N, centroid_x, centroid_y, centroid_z);
    }
    catch (Arg_error e)
    {
        delete alternaria;
        throw e;
    }

    alternaria->update_scene(psf, imaging);

    return alternaria;
}
#endif

/* -------------------------------------------------------------------------- */










/* ----------------------------  Alternaria_model  -------------------------- */

#if defined ALTERNARIA_HAVE_FFTW3F
Alternaria_model::Alternaria_model
(
    float                        c_1,
    float                        theta,
    float                        psi,
    size_t                       base_level,
    const Alternaria_density*    alt_d,
    const Spore_density*         spore_d,
    const Apical_hypha_density*  apical_hypha_d,
    const Lateral_hypha_density* lateral_hypha_1_d,
    const Lateral_hypha_density* lateral_hypha_n_d,
    const PSF_model*             psf,
    const Imaging_model*         imaging
)
throw (Arg_error) 
: Alternaria(c_1, theta, psi, base_level, alt_d, spore_d, apical_hypha_d, 
        lateral_hypha_1_d, lateral_hypha_n_d)
{
    scene = 0;
    s_1   = 0;
    s_2   = 0;
    s_3   = 0;
    s_4   = 0;

    update_scene(psf, imaging);
}
#endif


#if defined ALTERNARIA_HAVE_FFTW3F
Alternaria_model::Alternaria_model(const Alternaria_model& model)
    : Alternaria(model)
{
    scene = 0;
    if (model.scene)
    {
        copy_matblock_f(&scene, model.scene);
    }
    s_1 = 0;
    if (model.s_1)
    {
        copy_matblock_cf(&s_1, model.s_1);
    }
    s_2 = 0;
    if (model.s_2)
    {
        copy_matblock_cf(&s_2, model.s_2);
    }
    s_3 = 0;
    if (model.s_3)
    {
        copy_matblock_f(&s_3, model.s_3);
    }
    s_4 = 0;
    if (model.s_4)
    {
        copy_matblock_cf(&s_4, model.s_4);
    }
}
#endif


#if defined ALTERNARIA_HAVE_FFTW3F
Alternaria_model::Alternaria_model
(
    const Alternaria&    alt,
    const PSF_model*     psf,
    const Imaging_model* imaging
) 
: Alternaria(alt)
{
    scene = 0;
    s_1   = 0;
    s_2   = 0;
    s_3   = 0;
    s_4   = 0;

    update_scene(psf, imaging);
}
#endif


#if defined ALTERNARIA_HAVE_FFTW3F
Alternaria_model::Alternaria_model
(
    const char*                  fname,
    const Alternaria_density*    alt_d,
    const Spore_density*         spore_d,
    const Apical_hypha_density*  apical_hypha_d,
    const Lateral_hypha_density* lateral_hypha_1_d,
    const Lateral_hypha_density* lateral_hypha_n_d,
    const PSF_model*             psf,
    const Imaging_model*         imaging
)
throw (IO_error, Arg_error) 
: Alternaria(0.0f, 0, 0, 0, alt_d, spore_d, apical_hypha_d, lateral_hypha_1_d,
        lateral_hypha_n_d) 
{
    ostringstream ost;
    istringstream ist;
    ifstream in;

    scene = 0;
    s_1   = 0;
    s_2   = 0;
    s_3   = 0;
    s_4   = 0;

    in.open(fname);
    if (in.fail())
    {
        ost << fname << ": Could not open file";
        throw IO_error(ost.str());
    }

    const char* member_value;

    if (!(member_value = get_member_value("C_1", in)))
    {
        ost << fname << ": Missing C_1";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> c_1;
    if (ist.fail())
    {
        ost << fname << ": Invalid C_1";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    check_c_1();

    if (!(member_value = get_member_value("Theta", in)))
    {
        ost << fname << ": Missing theta";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> theta;
    if (ist.fail())
    {
        ost << fname << ": Invalid theta";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    if (!(member_value = get_member_value("Psi", in)))
    {
        ost << fname << ": Missing psi";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> psi;
    if (ist.fail())
    {
        ost << fname << ": Invalid psi";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    if (!(member_value = get_member_value("Level", in)))
    {
        ost << fname << ": Missing base level";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> base_level;
    if (ist.fail())
    {
        ost << fname << ": Invalid base level";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    if (!(member_value = get_member_value("Log Prob", in)))
    {
        ost << fname << ": Missing log probability";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> log_prob;
    if (ist.fail())
    {
        ost << fname << ": Invalid log probability";
        throw Arg_error(ost.str());
    }

    try
    {
        read(in);
    }
    catch (Arg_error e)
    {
        ost << fname << ": " << e.get_msg();
        throw Arg_error(ost.str());
    }
    catch (IO_error e)
    {
        ost << fname << ": " << e.get_msg();
        throw IO_error(ost.str());
    }

    in.close();
    if (in.fail())
    {
        ost << fname << ": Could not close file";
        throw IO_error(ost.str());
    }

    update_scene(psf, imaging);
}
#endif


#if defined ALTERNARIA_HAVE_FFTW3F
Alternaria_model::~Alternaria_model()
{
    free_matblock_f(scene);
    free_matblock_cf(s_1);
    free_matblock_cf(s_2);
    free_matblock_f(s_3);
    free_matblock_cf(s_4);
}
#endif


#if defined ALTERNARIA_HAVE_FFTW3F
Alternaria_model* Alternaria_model::clone() const
{
    return new Alternaria_model(*this);
}
#endif


#if defined ALTERNARIA_HAVE_FFTW3F
Alternaria_model& Alternaria_model::operator= 
(
    const Alternaria_model& alternaria
)
{
    Alternaria::operator=(alternaria);

    if (alternaria.scene)
    {
        copy_matblock_f(&scene, alternaria.scene);
    }
    if (alternaria.s_1)
    {
        copy_matblock_cf(&s_1, alternaria.s_1);
    }
    if (alternaria.s_2)
    {
        copy_matblock_cf(&s_2, alternaria.s_2);
    }
    if (alternaria.s_3)
    {
        copy_matblock_f(&s_3, alternaria.s_3);
    }
    if (alternaria.s_4)
    {
        copy_matblock_cf(&s_4, alternaria.s_4);
    }

    return *this;
}
#endif


#if defined ALTERNARIA_HAVE_FFTW3F
double Alternaria_model::calc_log_likelihood
(
    const PSF_model*        psf,
    const Imaging_model*    imaging,
    const jwsc::Matblock_f* data,
    float                   data_avg
) 
const throw (Arg_error)
{
    static const double log_c = -0.91893853321;

    uint32_t num_mats = imaging->get_window()->get_num_imgs();
    uint32_t num_rows = imaging->get_window()->get_num_rows();
    uint32_t num_cols = imaging->get_window()->get_num_cols();

    uint32_t fft_padding  = psf->get_padding()->fft;
    uint32_t conv_padding = psf->get_padding()->conv;

    uint32_t scene_start_mat = fft_padding + conv_padding;
    uint32_t scene_start_row = fft_padding + conv_padding;
    uint32_t scene_start_col = fft_padding + conv_padding;

    uint32_t data_start_mat = imaging->get_window()->get_img();
    uint32_t data_start_row = imaging->get_window()->get_row();
    uint32_t data_start_col = imaging->get_window()->get_col();

    double ll  = 0;
    double c_2 = imaging->get_c_2();

    uint32_t scene_mat = scene_start_mat; 
    uint32_t data_mat  = data_start_mat;
    while (scene_mat < (num_mats+scene_start_mat))
    {
        uint32_t scene_row = scene_start_row; 
        uint32_t data_row  = data_start_row;
        while (scene_row < (num_rows+scene_start_row))
        {
            uint32_t scene_col = scene_start_col; 
            uint32_t data_col  = data_start_col;
            while (scene_col < (num_cols+scene_start_col))
            {
                if (data_mat >= data->num_mats ||
                    data_row >= data->num_rows ||
                    data_col >= data->num_cols)
                {
                    throw Arg_error("Data not large enough for window");
                }

                double x     = data->elts[data_mat][data_row][data_col];
                double mu    = scene->elts[scene_mat][scene_row][scene_col];
                double sigma = sqrt(c_1*fabs(x - data_avg) + c_2*x);

                if (sigma < 1.0e-10) 
                {
                    std::cerr << "calc_log_likelihood: Sigma too small: " 
                              << sigma << '\n';
                    sigma = 1.0e-10;
                }

                ll += log_c - log(sigma) - 0.5*(x - mu)*(x - mu)/(sigma*sigma);

                scene_col++; 
                data_col++;
            }

            scene_row++; 
            data_row++;
        }

        scene_mat++; 
        data_mat++;
    }

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

    return ll;
}
#endif


#if defined ALTERNARIA_HAVE_FFTW3F
bool Alternaria_model::update_scene
(
    const PSF_model*     psf, 
    const Imaging_model* imaging
)
{
    float x_scale = imaging->get_scale()->get_x();
    float y_scale = imaging->get_scale()->get_y();
    float z_scale = imaging->get_scale()->get_z();

    const Imaging_window* window = imaging->get_window();
    const PSF_padding* padding = psf->get_padding();

    uint32_t num_mats = window->get_num_imgs() 
        + 2*(padding->conv + padding->fft);
    uint32_t num_rows = window->get_num_rows() 
        + 2*(padding->conv + padding->fft);
    uint32_t num_cols = window->get_num_cols() 
        + 2*(padding->conv + padding->fft);

    create_zero_matblock_f(&scene, num_mats, num_rows, num_cols);

    float x = x_scale*((float)window->get_col() - padding->conv - padding->fft);
    float y = y_scale*((float)window->get_row() - padding->conv - padding->fft);
    float z = z_scale*((float)window->get_img() - padding->conv - padding->fft);

    bool intersect = false;

    if (root)
    {
        intersect = root->recursively_draw_in_matblock_f(scene, x, y, z, 
                x_scale, y_scale, z_scale);
    }

    float bg = imaging->get_bg();

    for (uint32_t mat = 0; mat < num_mats; mat++)
    {
        for (uint32_t row = 0; row < num_rows; row++)
        {
            for (uint32_t col = 0; col < num_cols; col++)
            {
                if ((mat < padding->fft) || (mat > num_mats-padding->fft-1) ||
                    (row < padding->fft) || (row > num_rows-padding->fft-1) ||
                    (col < padding->fft) || (col > num_cols-padding->fft-1))
                {
                    scene->elts[ mat ][ row ][ col ] = bg;
                }
                else
                {
                    float id = floor(scene->elts[ mat ][ row ][ col ]);
                    float opacity = scene->elts[ mat ][ row ][ col ] - id;

                    scene->elts[ mat ][ row ][ col ] = bg * (1 - opacity);
                }
            }
        }
    }

    if (root && psf->get_h())
    {
        assert(convolve_matblock_fft_scratch_f(&scene, scene, psf->get_h(), 
                    &s_1, &s_2, &s_3, &s_4) == 0);
    }

    return intersect;
}
#endif


#if defined ALTERNARIA_HAVE_FFTW3F
void Alternaria_model::write_scene
(
    const char*          fname_fmt,
    const PSF_model*     psf
) 
const throw (IO_error)
{
    if (scene == 0)
        return;

    uint32_t scene_num_mats = scene->num_mats;
    uint32_t scene_num_rows = scene->num_rows;
    uint32_t scene_num_cols = scene->num_cols;

    const PSF_padding* padding = psf->get_padding();

    uint32_t num_mats = scene_num_mats - 2*(padding->fft + padding->conv);
    uint32_t num_rows = scene_num_rows - 2*(padding->fft + padding->conv);
    uint32_t num_cols = scene_num_cols - 2*(padding->fft + padding->conv);

    uint32_t start_mat = padding->fft + padding->conv;
    uint32_t start_row = padding->fft + padding->conv;
    uint32_t start_col = padding->fft + padding->conv;

    Matblock_f* scene_no_pad = 0;
    copy_matblock_block_f(&scene_no_pad, scene, start_mat, start_row, start_col,
            num_mats, num_rows, num_cols);

    Imgblock_f* imgs = 0;
    create_imgblock_from_matblock_f(&imgs, scene_no_pad);
    free_matblock_f(scene_no_pad);

    Error* e;
    if ((e = write_imgblock_numbered_f(imgs, fname_fmt)) != 0)
    {
        IO_error ex(e->msg);
        free_error(e);
        throw ex;
    }
    free_imgblock_f(imgs);
}
#endif
/* -------------------------------------------------------------------------- */