structure.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>

#if defined ALTERNARIA_HAVE_OPENGL_FRAMEWORK
#include <OpenGL/gl.h>
#include <OpenGL/glu.h>
#elif defined ALTERNARIA_HAVE_OPENGL
#include <GL/gl.h>
#include <GL/glu.h>
#endif

#include <cstring>
#include <cmath>
#include <cctype>
#include <cassert>
#include <list>
#include <vector>
#include <iostream>
#include <iomanip>
#include <sstream>

#include <inttypes.h>

#include <jwsc/base/limits.h>
#include <jwsc/math/constants.h>
#include <jwsc/prob/pdf.h>
#include <jwsc/vector/vector.h>
#include <jwsc/vector/vector_math.h>
#include <jwsc/matrix/matrix.h>
#include <jwsc/matrix/matrix_math.h>
#include <jwsc/matblock/matblock.h>

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

#include "density.h"
#include "printable.h"
#include "structure.h"


using std::list;
using std::vector;
using std::ios_base;
using std::istream;
using std::ostringstream;
using std::istringstream;
using std::streampos;
using std::isspace;
using namespace jwsc;
using jwscxx::base::Arg_error;
using jwscxx::base::Result_error;
using jwscxx::base::IO_error;


/* -------------------------  Structure_density  ---------------------------- */

void Structure_density::set_centroid_x(float min, float max) throw (Arg_error)
{
    centroid_x.set(min, max);
}


void Structure_density::set_centroid_y(float min, float max) throw (Arg_error)
{
    centroid_y.set(min, max);
}


void Structure_density::set_centroid_z(float min, float max) throw (Arg_error)
{
    centroid_z.set(min, max);
}


void Structure_density::set_length
(
    float mu, 
    float sigma, 
    float min, 
    float max
) 
throw (Arg_error)
{
    if (min < 0)
    {
        throw Arg_error("Minimum must be >= 0");
    }
    length.set(mu, sigma, min, max);
}


void Structure_density::set_width
(
    float mu, 
    float sigma, 
    float min, 
    float max
) 
throw (Arg_error)
{
    if (min < 0)
    {
        throw Arg_error("Minimum must be >= 0");
    }
    width.set(mu, sigma, min, max);
}


void Structure_density::set_theta
(
    float mu, 
    float sigma, 
    float min, 
    float max
) 
throw (Arg_error)
{
    if (min < 0)
    {
        throw Arg_error("Minimum must be >= 0");
    }
    if (max > JWSC_PI)
    {
        throw Arg_error("Maximum must be <= PI");
    }
    theta.set(mu, sigma, min, max);
}


void Structure_density::set_psi(float min, float max) throw (Arg_error)
{
    if (min < -JWSC_PI)
    {
        throw Arg_error("Minimum must be >= PI");
    }
    if (max > JWSC_PI)
    {
        throw Arg_error("Maximum must be <= PI");
    }
    psi.set(min, max);
}


void Structure_density::set_opacity
(
    float mu, 
    float sigma, 
    float min, 
    float max
) 
throw (Arg_error)
{
    if (min < 0 || max > 1)
    {
        throw Arg_error("Parameter range must be in [0,1]");
    }
    opacity.set(mu, sigma, min, max);
}


float Structure_density::sample_centroid_x() const
{
    return (float)sample_uniform_pdf_d(centroid_x.get_min(), 
            centroid_x.get_max());
}


float Structure_density::sample_centroid_y() const
{
    return (float)sample_uniform_pdf_d(centroid_y.get_min(), 
            centroid_y.get_max());
}


float Structure_density::sample_centroid_z() const
{
    return (float)sample_uniform_pdf_d(centroid_z.get_min(), 
            centroid_z.get_max());
}


float Structure_density::sample_length() const
{
    return (float)sample_gaussian_pdf_d(length.get_mu(), length.get_sigma(), 
            length.get_min(), length.get_max());
}


float Structure_density::sample_width() const
{
    return (float)sample_gaussian_pdf_d(width.get_mu(), width.get_sigma(), 
            width.get_min(), width.get_max());
}


float Structure_density::sample_theta() const
{
    return (float)sample_gaussian_pdf_d(theta.get_mu(), 
            theta.get_sigma(), theta.get_min(), theta.get_max());
}


float Structure_density::sample_psi() const
{
    return (float)sample_uniform_pdf_d(psi.get_min(), psi.get_max());
}


float Structure_density::sample_opacity() const
{
    return (float)sample_gaussian_pdf_d(opacity.get_mu(), 
            opacity.get_sigma(), opacity.get_min(), opacity.get_max());
}

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









/* -------------------------  Lateral_structure_density  -------------------- */

void Lateral_structure_density::set_lat_dist
(
    float mu, 
    float sigma, 
    float min, 
    float max
) 
throw (Arg_error)
{
    if (min < 0 || max > 1)
    {
        throw Arg_error("Parameter range must be in [0,1]");
    }
    lat_dist.set(mu, sigma, min, max);
}


float Lateral_structure_density::sample_lat_dist() const
{
    return (float)sample_gaussian_pdf_d(lat_dist.get_mu(), 
            lat_dist.get_sigma(), lat_dist.get_min(), lat_dist.get_max());
}

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







/* ------------------------------  Structure  ------------------------------- */

Structure::Structure
(
    float length,
    float width,
    float base_theta,
    float base_psi,
    float theta,
    float psi,
    float opacity,
    size_t level,
    const class Structure_density* density
)
throw (Arg_error) 
    : laterals()
{
    this->length     = length;
    this->width      = width;
    this->base_theta = base_theta;
    this->base_psi   = base_psi;
    this->theta      = theta;
    this->psi        = psi;
    this->opacity    = opacity;
    this->level      = level;
    this->density    = density;

    check_length();
    check_width();
    check_base_theta();
    check_base_psi();
    check_theta();
    check_psi();
    check_opacity();

    parent       = 0;
    apical       = 0;
    begin_pt     = 0;
    end_pt       = 0;
    growth_dir   = 0;
    centroid     = 0;
    R            = 0;
    log_prob     = 0;
    log_prob_descend = 0;

    id = generate_id();
}


Structure::Structure(const Structure& s) : laterals(s.laterals)
{
    id         = s.id;
    level      = s.level;
    length     = s.length;
    width      = s.width;
    base_theta = s.base_theta;
    base_psi   = s.base_psi;
    theta      = s.theta;
    psi        = s.psi;
    opacity    = s.opacity;
    density    = s.density;
    parent     = s.parent;
    apical     = s.apical;

    begin_pt = 0;
    copy_vector_f(&begin_pt, s.begin_pt);

    end_pt = 0;
    copy_vector_f(&end_pt, s.end_pt);

    growth_dir = 0;
    copy_vector_f(&growth_dir, s.growth_dir);

    centroid = 0;
    copy_vector_f(&centroid, s.centroid);

    R = 0;
    copy_matrix_f(&R, s.R);

    log_prob = s.log_prob;
    log_prob_descend = s.log_prob_descend;
}


Structure::~Structure()
{
    free_vector_f(begin_pt);
    free_vector_f(end_pt);
    free_vector_f(centroid);
    free_matrix_f(R);
    free_vector_f(growth_dir);

    delete apical;

    for (list<Lateral_structure*>::iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        delete *it;
    }
}


Structure& Structure::operator= (const Structure& s)
{
    id         = s.id;
    level      = s.level;
    length     = s.length;
    width      = s.width;
    base_theta = s.base_theta;
    base_psi   = s.base_psi;
    theta      = s.theta;
    psi        = s.psi;
    opacity    = s.opacity;
    density    = s.density;
    parent     = s.parent;
    apical     = s.apical;

    laterals.clear();
    laterals = s.laterals;

    copy_vector_f(&begin_pt, s.begin_pt);

    copy_vector_f(&end_pt, s.end_pt);

    copy_vector_f(&growth_dir, s.growth_dir);

    copy_vector_f(&centroid, s.centroid);

    copy_matrix_f(&R, s.R);

    log_prob = s.log_prob;
    log_prob_descend = s.log_prob_descend;

    return *this;
}


Structure* Structure::get_root()
{
    if (parent)
        return parent->get_root();

    return this;
}


vector<Structure*>* Structure::get_descendants()
{
    list<Structure*>* descendants = new list<Structure*>();
    descendants->push_back(this);

    if (apical)
    {
        vector<Structure*>* apical_descend = apical->get_descendants();
        descendants->insert(descendants->end(), apical_descend->begin(), 
                apical_descend->end());
        delete apical_descend;
    }
    for (list<Lateral_structure*>::const_iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        vector<Structure*>* lateral_descend = (*it)->get_descendants();
        descendants->insert(descendants->end(), lateral_descend->begin(),
                lateral_descend->end());
        delete lateral_descend;
    }

    vector<Structure*>* descendants_v = new vector<Structure*>(
            descendants->begin(), descendants->end());
    delete descendants;

    return descendants_v;
}


void Structure::remove_lateral(class Lateral_structure* s) throw (Arg_error)
{
    for (list<Lateral_structure*>::iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        if ((*it) == s)
        {
            laterals.erase(it);
            update_ancestry_log_prob();
            return;
        }
    }

    throw Arg_error("Could not find lateral to remove");
}


void Structure::set_length(float length) throw (Arg_error)
{
    set_size(length, width);
}


void Structure::set_width(float width) throw (Arg_error)
{
    set_size(length, width);
}


void Structure::set_size(float length, float width) throw (Arg_error)
{
    static int entry = 0;

    entry++;

    float old_length = this->length;
    float old_width  = this->width;

    try
    {
        this->length = length;
        check_length();

        this->width  = width;
        check_width();

        update_centroid();
        update_end_pt();
        update_apical();
        update_laterals();
        update_ancestry_log_prob();
    }
    catch (Arg_error e)
    {
        if (entry < 2)
        {
            set_size(old_length, old_width);
        }
        else
        {
            // TODO the length has been coming up wrong because of precision.
            // Probably could throw something like Result_error exception
            // here.
            entry = 0;
        }
        throw e;
    }
}


void Structure::set_theta(float theta) throw (Arg_error)
{
    set_orientation(theta, psi);
}


void Structure::set_base_theta(float base_theta) throw (Arg_error)
{
    set_base_orientation(base_theta, base_psi);
}


void Structure::set_psi(float psi) throw (Arg_error)
{
    set_orientation(theta, psi);
}


void Structure::set_base_psi(float base_psi) throw (Arg_error)
{
    set_base_orientation(base_theta, base_psi);
}


void Structure::set_orientation(float theta, float psi) throw (Arg_error)
{
    static int entry = 0;

    entry++;

    float old_theta = this->theta;
    float old_psi   = this->psi;

    try
    {
        this->theta = theta;
        check_theta();

        this->psi = psi;
        check_psi();

        update_R();
        update_growth_dir_from_R();
        update_centroid();
        update_end_pt();
        update_apical();
        update_laterals();
        update_ancestry_log_prob();
    }
    catch (Arg_error e)
    {
        if (entry < 2)
        {
            set_orientation(old_theta, old_psi);
        }
        else
        {
            // TODO the length has been coming up wrong because of precision.
            // Probably could throw something like Result_error exception
            // here.
            entry = 0;
        }
        throw e;
    }
}


void Structure::set_base_orientation(float base_theta, float base_psi) 
    throw (Arg_error)
{
    static int entry = 0;

    entry++;

    float old_base_theta = this->base_theta;
    float old_base_psi   = this->base_psi;

    try
    {
        this->base_theta = base_theta;
        check_base_theta();

        this->base_psi = base_psi;
        check_base_psi();

        update_R();
        update_growth_dir_from_R();
        update_centroid();
        update_end_pt();
        update_apical();
        update_laterals();
        update_ancestry_log_prob();
    }
    catch (Arg_error e)
    {
        if (entry < 2)
        {
            set_base_orientation(old_base_theta, old_base_psi);
        }
        else
        {
            // TODO the length has been coming up wrong because of precision.
            // Probably could throw something like Result_error exception
            // here.
            entry = 0;
        }
        throw e;
    }
}


void Structure::set_opacity(float opacity) throw (Arg_error)
{
    static int entry = 0;

    entry++;

    float old_opacity = this->opacity;

    try
    {
        this->opacity = opacity;
        check_opacity();

        update_ancestry_log_prob();
    }
    catch (Arg_error e)
    {
        if (entry < 2)
        {
            set_opacity(old_opacity);
        }
        else
        {
            // TODO the length has been coming up wrong because of precision.
            // Probably could throw something like Result_error exception
            // here.
            entry = 0;
        }
        throw e;
    }
}


bool Structure::recursively_draw_in_matblock_f
(
     jwsc::Matblock_f* M_blk, 
     float             x, 
     float             y,
     float             z,
     float             x_scale,
     float             y_scale,
     float             z_scale,
     bool              fill
)
const
{
    bool intersect = false;

    if (apical)
    {
        if (apical->recursively_draw_in_matblock_f(M_blk, x, y, z, x_scale, 
                    y_scale, z_scale, fill))
        {
            intersect = true;
        }
    }
    for (list<Lateral_structure*>::const_iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        if ((*it)->recursively_draw_in_matblock_f(M_blk, x, y, z, x_scale, 
                    y_scale, z_scale, fill))
        {
            intersect = true;
        }
    }

    if (draw_in_matblock_f(M_blk, x, y, z, x_scale, y_scale, z_scale, fill))
    {
        intersect = true;
    }

    return intersect;
}


#if defined ALTERNARIA_HAVE_OPENGL_FRAMEWORK || defined ALTERNARIA_HAVE_OPENGL
void Structure::recursively_draw_in_opengl(GLUquadric* quad, float scale)
const
{
    draw_in_opengl(quad, scale);

    if (apical)
    {
        apical->recursively_draw_in_opengl(quad, scale);
    }
    for (list<Lateral_structure*>::const_iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        (*it)->recursively_draw_in_opengl(quad, scale);
    }
}
#endif


Structure* Structure::get_terminal_apical()
{
    if (apical)
    {
        return apical->get_terminal_apical();
    }
    else
    {
        return this;
    }
}


Lateral_structure* Structure::get_random_lateral()
    throw (Result_error)
{
    if (laterals.empty())
    {
        throw Result_error("Structure contains no laterals");
    }

    int u = static_cast<int>(
            floor(sample_uniform_pdf_d(0, 0.9999*laterals.size())));

    list<Lateral_structure*>::const_iterator it = laterals.begin();

    for (int i = 0; i < u; i++)
    {
        it++;
    }

    return *it;
}


void Structure::print(ostream& out) const
{
    out << "  Structure: " << get_type_str() << '\n'
        << "         ID: " << id << '\n'
        << "      Level: " << level << '\n'
        << "   Centroid: " << centroid->elts[0] << ' '
                          << centroid->elts[1] << ' '
                          << centroid->elts[2] << '\n'
        << "     Length: " << length << '\n'
        << "      Width: " << width << '\n'
        << " Base Theta: " << base_theta << '\n'
        << "   Base Psi: " << base_psi << '\n'
        << "      Theta: " << theta << '\n'
        << "        Psi: " << psi << '\n'
        << "    Opacity: " << opacity << '\n';

    out << "     Parent: " ;
    if (parent)
    {
        out << parent->id;
    }
    out << '\n';

    out << "     Apical: " ;
    if (apical)
    {
        out << apical->id;
    }
    out << '\n';

    out << "   Laterals: ";
    for (list<Lateral_structure*>::const_iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        out << (*it)->id << ' ';
    }
    out << '\n';

    out << "   Log Prob: " << std::setprecision(10) << log_prob << '\n';
}


void Structure::recursively_print(ostream& out) const
{
    print(out);

    if (apical)
    {
        apical->recursively_print(out);
    }
    for (list<Lateral_structure*>::const_iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        (*it)->recursively_print(out);
    }
}


void Structure::recursively_check_structure() const
{
    check_width();
    check_length();
    check_theta();
    check_psi();
    check_base_theta();
    check_base_psi();
    check_centroid();
    check_opacity();

    if (apical)
    {
        assert(fabs(end_pt->elts[0] - apical->begin_pt->elts[0]) < 1.0e-6);
        assert(fabs(end_pt->elts[1] - apical->begin_pt->elts[1]) < 1.0e-6);
        assert(fabs(end_pt->elts[2] - apical->begin_pt->elts[2]) < 1.0e-6);
        apical->recursively_check_structure();
    }
    for (list<Lateral_structure*>::const_iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        (*it)->recursively_check_structure();
    }
}


bool Structure::has_child_id(uint32_t id) const
{
    if (apical && apical->id == id)
    {
        return true;
    }

    for (list<Lateral_structure*>::const_iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        if ((*it)->id == id)
        {
            return true;
        }
    }

    return false;
}


void Structure::update_angles() throw (Arg_error)
{
    Matrix_f* M = 0;
    Vector_f* v = 0;

    if (parent)
    {
        transpose_matrix_f(&M, parent->R);
        multiply_matrix_and_vector_f(&v, M, growth_dir);
    }
    else
    {
        copy_vector_f(&v, growth_dir);
    }

    if (base_theta != 0 || base_psi != 0)
    {
        create_euler_rotation_matrix_f(&M, 0, base_theta, base_psi);
        transpose_matrix_f(&M, M);
        multiply_matrix_and_vector_f(&v, M, v);
    }

    double x = v->elts[0];
    double y = v->elts[1];
    double z = v->elts[2];

    bool zero_theta = false;

    // Calculate the euler angle psi: rotation about z-axis to put the apical 
    // growth_dir on the z,(y+) plane.
    //
    // Check some special cases.
    if (fabs(x) < 1.0e-6 && y > -1.0e-6)
    {
        psi = 0;
        if (fabs(y) < 1.0e-6)
        { 
            zero_theta = true; 
        }
    }
    else if (fabs(x) < 1.0e-6)
    {
        psi = JWSC_PI;
    }
    else
    {
        double s = std::sqrt(x*x + y*y);
        double d = y / s;

        if (d > 1)
            d = 1;
        else if (d < -1)
            d = -1;

        psi = (x > 0) ? std::acos(d) : -std::acos(d);

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

    float min = density->get_psi().get_min();
    float 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;
    }

    // Undo the psi rotation to put v on the z,(y+) plane.
    if (zero_theta)
    {
        theta = 0;
    }
    else
    {
        create_euler_rotation_matrix_f(&M, 0, 0, psi);
        transpose_matrix_f(&M, M);
        multiply_matrix_and_vector_f(&v, M, v);

        y = v->elts[1];
        z = v->elts[2];

        // Calculate the euler angle theta: rotation about the x-axis to put
        // apical growth_dir on the z-axis.
        assert(y > 0);

        if (z > 1)
            z = 1;

        theta = std::acos(z);

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

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

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

    free_matrix_f(M);
    free_vector_f(v);

    check_theta();
    check_psi();
}


void Structure::update_R()
{
    create_euler_rotation_matrix_f(&R, 0, theta, psi);

    if (base_theta != 0 || base_psi != 0)
    {
        Matrix_f* base_R = 0;
        create_euler_rotation_matrix_f(&base_R, 0, base_theta, base_psi);
        multiply_matrices_f(&R, base_R, R);
        free_matrix_f(base_R);
    }

    if (parent)
    {
        multiply_matrices_f(&R, parent->R, R);
    }
}


void Structure::update_growth_dir_from_R()
{
    create_zero_vector_f(&growth_dir, 3);
    growth_dir->elts[2] = 1;

    multiply_matrix_and_vector_f(&growth_dir, R, growth_dir);
}


void Structure::update_growth_dir_and_length_from_pts() throw (Arg_error)
{
    subtract_vectors_f(&growth_dir, end_pt, begin_pt);
    calc_vector_mag_f(&length, growth_dir);
    normalize_vector_mag_f(&growth_dir, growth_dir);

    check_length();
}


void Structure::update_centroid() throw (Arg_error)
{
    multiply_vector_by_scalar_f(&centroid, growth_dir, 0.5f*length);
    add_vectors_f(&centroid, centroid, begin_pt);

    check_centroid();
}


void Structure::update_begin_pt()
{
    multiply_vector_by_scalar_f(&begin_pt, growth_dir, -0.5*length);
    add_vectors_f(&begin_pt, begin_pt, centroid);
}


void Structure::update_end_pt()
{
    multiply_vector_by_scalar_f(&end_pt, growth_dir, 0.5*length);
    add_vectors_f(&end_pt, end_pt, centroid);
}


void Structure::update_apical() throw (Arg_error)
{
    if (!apical)
        return;

    apical->update_from_parent_end_pt();
}


void Structure::update_laterals() throw (Arg_error)
{
    for (list<Lateral_structure*>::iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        (*it)->update_from_parent_end_pt();
    }
}


void Structure::update_log_prob()
{
    float    mu, sigma;
    float    min, max;
    double   delta;
    double   p_1, p_2;

    log_prob = 0;
    log_prob_descend = 0;

    if (apical)
    {
        log_prob_descend += apical->log_prob + apical->log_prob_descend;
    }

    for (list<Lateral_structure*>::iterator it = laterals.begin();
            it != laterals.end(); it++)
    {
        log_prob_descend += (*it)->log_prob + (*it)->log_prob_descend;
    }

    /* Length */
    mu        = density->get_length().get_mu();
    sigma     = density->get_length().get_sigma();
    min       = density->get_length().get_min();
    max       = density->get_length().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, length, length+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);

    /* Width */
    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);

    /* Theta */
    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);

    /* Psi */
    min       = density->get_psi().get_min();
    max       = density->get_psi().get_max();
    log_prob -= log(max - min);

    /* Opacity */
    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));
    assert(!std::isinf(log_prob_descend) && !std::isnan(log_prob_descend));
#endif
}


void Structure::update_ancestry_log_prob()
{
    update_log_prob();

    if (parent)
    {
        parent->update_ancestry_log_prob();
    }
}


void Structure::recursively_update_levels(int level_delta) throw (Result_error)
{
    if (level_delta == 0)
    {
        return;
    }
    else if (level_delta < 0 && level < static_cast<size_t>(abs(level)))
    {
        throw Result_error("Structure level < 0");
    }
    else
    {
        if (apical)
        {
            apical->recursively_update_levels(level_delta);
        }
        for (list<Lateral_structure*>::iterator it = laterals.begin();
                it != laterals.end(); it++)
        {
            (*it)->recursively_update_levels(level_delta);
        }
        level += level_delta;
    }
}


void Structure::check_length() const throw (Arg_error)
{
    if (density->get_length().get_min() > length || 
        density->get_length().get_max() < length)
    {
        ostringstream ost;
        ost << "Length " << length << " outside density range ["
            << density->get_length().get_min() << ','
            << density->get_length().get_max() << ']';
        throw Arg_error(ost.str());
    }
}


void Structure::check_width() const throw (Arg_error)
{
    if (density->get_width().get_min() > width || 
        density->get_width().get_max() < width)
    {
        ostringstream ost;
        ost << "Width " << width << " outside density range ["
            << density->get_width().get_min() << ','
            << density->get_width().get_max() << ']';
        throw Arg_error(ost.str());
    }
}


void Structure::check_theta() const throw (Arg_error)
{
    float min = density->get_theta().get_min();
    float max = density->get_theta().get_max();

    if (min > theta || max < theta)
    {
        ostringstream ost;
        ost << "Theta " << theta << " outside density range ["
            << min << ',' << max << ']';
        throw Arg_error(ost.str());
    }
}


void Structure::check_base_theta() const throw (Arg_error)
{
    float min = 0;
    float max = JWSC_PI;

    if (min > base_theta || max < base_theta)
    {
        ostringstream ost;
        ost << "Base theta " << base_theta << " outside density range ["
            << min << ',' << max << ']';
        throw Arg_error(ost.str());
    }
}


void Structure::check_psi() const throw (Arg_error)
{
    float min = density->get_psi().get_min();
    float max = density->get_psi().get_max();

    if (min > psi || max < psi)
    {
        ostringstream ost;
        ost << "Psi " << psi << " outside density range ["
            << min << ',' << max << ']';
        throw Arg_error(ost.str());
    }
}


void Structure::check_base_psi() const throw (Arg_error)
{
    float min = -JWSC_PI;
    float max = JWSC_PI;

    if (min > base_psi || max < base_psi)
    {
        ostringstream ost;
        ost << "Base psi " << base_psi << " outside density range ["
            << min << ',' << max << ']';
        throw Arg_error(ost.str());
    }
}


void Structure::check_opacity() const throw (Arg_error)
{
    if (density->get_opacity().get_min() > opacity || 
        density->get_opacity().get_max() < opacity)
    {
        ostringstream ost;
        ost << "Opacity " << opacity << " outside density range ["
            << density->get_opacity().get_min() << ','
            << density->get_opacity().get_max() << ']';
        throw Arg_error(ost.str());
    }
}


uint32_t Structure::generate_id()
{
    static uint32_t next_id = 1;
    return next_id++;
}


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







/* -------------------------  Printed_structure  ---------------------------- */

void Printed_structure::read_id(std::istream& in) throw (IO_error, Arg_error)
{
    const char* member_value;

    if (!(member_value = Structure::get_member_value("ID", in)))
    {
        throw Arg_error("Structure missing ID");
    }
    istringstream ist(member_value);
    ist >> id;
    if (ist.fail())
    {
        throw Arg_error("Structure has invalid ID");
    }
}


void Printed_structure::read_level(std::istream& in) throw (IO_error, Arg_error)
{
    const char* member_value;

    if (!(member_value = Structure::get_member_value("Level", in)))
    {
        throw Arg_error("Structure missing level");
    }
    istringstream ist(member_value);
    ist >> level;
    if (ist.fail())
    {
        throw Arg_error("Structure has invalid level");
    }
}


void Printed_structure::read_centroid(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    // Centroid X
    if (!(member_value = Structure::get_member_value("Centroid", in)))
    {
        throw Arg_error("Missing centroid x");
    }
    istringstream ist(member_value);
    ist >> centroid_x;
    if (ist.fail())
    {
        throw Arg_error("Invalid centroid x");
    }
    ist.clear(std::ios_base::goodbit);

    // Centroid Y
    if (!(member_value = Structure::get_next_member_value(member_value)))
    {
        throw Arg_error("Missing centroid y");
    }
    ist.str(member_value);
    ist >> centroid_y;
    if (ist.fail())
    {
        throw Arg_error("Invalid centroid y");
    }
    ist.clear(std::ios_base::goodbit);

    // Centroid Z
    if (!(member_value = Structure::get_next_member_value(member_value)))
    {
        throw Arg_error("Missing centroid z");
    }
    ist.str(member_value);
    ist >> centroid_z;
    if (ist.fail())
    {
        throw Arg_error("Invalid centroid z");
    }
    ist.clear(std::ios_base::goodbit);
}


void Printed_structure::read_length(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    if (!(member_value = Structure::get_member_value("Length", in)))
    {
        throw Arg_error("Missing length");
    }
    istringstream ist(member_value);
    ist >> length;
    if (ist.fail())
    {
        throw Arg_error("Invalid length");
    }
}


void Printed_structure::read_width(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    if (!(member_value = Structure::get_member_value("Width", in)))
    {
        throw Arg_error("Missing width");
    }
    istringstream ist(member_value);
    ist >> width;
    if (ist.fail())
    {
        throw Arg_error("Invalid width");
    }
}


void Printed_structure::read_theta(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    if (!(member_value = Structure::get_member_value("Theta", in)))
    {
        throw Arg_error("Missing theta");
    }
    istringstream ist(member_value);
    ist >> theta;
    if (ist.fail())
    {
        throw Arg_error("Invalid theta");
    }
}


void Printed_structure::read_base_theta(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    if (!(member_value = Structure::get_member_value("Base Theta", in)))
    {
        throw Arg_error("Missing base theta");
    }
    istringstream ist(member_value);
    ist >> base_theta;
    if (ist.fail())
    {
        throw Arg_error("Invalid base theta");
    }
}


void Printed_structure::read_psi(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    if (!(member_value = Structure::get_member_value("Psi", in)))
    {
        throw Arg_error("Missing psi");
    }
    istringstream ist(member_value);
    ist >> psi;
    if (ist.fail())
    {
        throw Arg_error("Invalid psi");
    }
}


void Printed_structure::read_base_psi(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    if (!(member_value = Structure::get_member_value("Base Psi", in)))
    {
        throw Arg_error("Missing base psi");
    }
    istringstream ist(member_value);
    ist >> base_psi;
    if (ist.fail())
    {
        throw Arg_error("Invalid base psi");
    }
}


void Printed_structure::read_opacity(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    if (!(member_value = Structure::get_member_value("Opacity", in)))
    {
        throw Arg_error("Missing opacity");
    }
    istringstream ist(member_value);
    ist >> opacity;
    if (ist.fail())
    {
        throw Arg_error("Invalid opacity");
    }
}


void Printed_structure::read_log_prob(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

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


void Printed_structure::read_parent(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    parent = 0;
    if ((member_value = Structure::get_member_value("Parent", in)))
    {
        istringstream ist(member_value);
        ist >> parent;
        if (ist.fail())
        {
            throw Arg_error("Invalid parent");
        }
    }
}


void Printed_structure::read_apical(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    apical = 0;
    if ((member_value = Structure::get_member_value("Apical", in)))
    {
        istringstream ist(member_value);
        ist >> apical;
        if (ist.fail())
        {
            throw Arg_error("Invalid apical");
        }
    }
}


void Printed_structure::read_laterals(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    laterals.clear();
    if ((member_value = Structure::get_member_value("Laterals", in)))
    {
        uint32_t lateral;
        istringstream ist(member_value);
        ist >> lateral;

        laterals.push_back(lateral);
        while ((member_value = Structure::get_next_member_value(member_value)))
        {
            ist.str(member_value);
            ist >> lateral;
            if (ist.fail())
            {
                throw Arg_error("Invalid lateral");
            }
            ist.clear(std::ios_base::goodbit);
            laterals.push_back(lateral);
        }
    }
}


bool Printed_structure::is_on_istream(std::istream& in, const char* type_str) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    streampos in_posn = in.tellg();

    if (!(member_value = Structure::get_member_value("Structure", in)))
    {
        throw Arg_error("Missing structure type");
    }

    in.seekg(in_posn);

    if (strncmp(member_value, type_str, strlen(type_str)) != 0)
    {
        return false;
    }

    return true;
}

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







/* -------------------------  Apical_structure  ----------------------------- */

Apical_structure::Apical_structure
(
    float centroid_x, 
    float centroid_y,
    float centroid_z,
    float length,
    float width,
    float base_theta,
    float base_psi,
    float theta,
    float psi,
    float opacity,
    size_t level,
    const class Apical_structure_density* density
) 
throw (Arg_error)
    : Structure(length, width, base_theta, base_psi, theta, psi, opacity, 
            level, density)
{
    this->density = density;

    create_vector_f(&centroid, 3);
    centroid->elts[0] = centroid_x;
    centroid->elts[1] = centroid_y;
    centroid->elts[2] = centroid_z;

    check_centroid();

    update_R();
    update_growth_dir_from_R();
    update_begin_pt();
    update_end_pt();
}


Apical_structure::Apical_structure
(
    Structure* parent,
    float      length,
    float      width,
    float      base_theta,
    float      base_psi,
    float      theta,
    float      psi,
    float      opacity,
    const class Apical_structure_density* density
)
throw (Arg_error, Apical_error)
    : Structure(length, width, base_theta, base_psi, theta, psi, opacity, 
            parent->level, density)
{
    if (parent->apical)
    {
        throw Apical_error("Structure already contains an apical growth");
    }

    this->parent  = parent;
    this->density = density;

    copy_vector_f(&begin_pt, parent->end_pt);

    update_R();
    update_growth_dir_from_R();
    update_centroid();
    update_end_pt();

    parent->apical = this;
}


Apical_structure::Apical_structure
(
    Structure* parent,
    float      centroid_x, 
    float      centroid_y,
    float      centroid_z,
    float      length,
    float      width,
    float      base_theta,
    float      base_psi,
    float      opacity,
    const class Apical_structure_density* density
)
throw (Arg_error, Apical_error)
    : Structure(length, width, base_theta, base_psi, parent->theta, 
            parent->psi, opacity, parent->level, density)
{
    if (parent->apical)
    {
        throw Apical_error("Structure already contains an apical growth");
    }

    this->parent  = parent;
    this->density = density;

    create_vector_f(&centroid, 3);
    centroid->elts[0] = centroid_x;
    centroid->elts[1] = centroid_y;
    centroid->elts[2] = centroid_z;
    check_centroid();

    copy_vector_f(&begin_pt, parent->end_pt);

    subtract_vectors_f(&growth_dir, centroid, begin_pt);
    normalize_vector_mag_f(&growth_dir, growth_dir);

    update_end_pt();
    update_angles();
    update_R();

    parent->apical = this;
}


Apical_structure::Apical_structure(const Apical_structure& s) : Structure(s)
{
    density  = s.density;
}


Apical_structure* Apical_structure::recursively_clone() const
{
    Apical_structure* s = clone();

    s->apical = 0;
    s->laterals.clear();

    if (apical)
    {
        s->apical = apical->recursively_clone();
        s->apical->parent = s;
    }
    for (list<Lateral_structure*>::const_iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        Lateral_structure* lateral = (*it)->recursively_clone();
        s->laterals.push_back(lateral);
        lateral->parent = s;
    }

    return s;
}


Apical_structure& Apical_structure::operator= (const Apical_structure& s)
{
    Structure::operator=(s);

    density = s.density;

    return *this;
}


void Apical_structure::set_position
(
    float centroid_x, 
    float centroid_y, 
    float centroid_z
) 
throw (jwscxx::base::Arg_error)
{
    static int entry = 0;

    entry++;

    Vector_f* old_centroid = 0;
    copy_vector_f(&old_centroid, centroid);

    Vector_f* T = 0;
    copy_vector_f(&T, centroid);
    T->elts[0] = centroid_x - T->elts[0];
    T->elts[1] = centroid_y - T->elts[1];
    T->elts[2] = centroid_z - T->elts[2];

    float mag;
    calc_vector_mag_f(&mag, T);
    if (mag < 1.0e-2)
    {
        free_vector_f(old_centroid);
        free_vector_f(T);
        return;
    }

    try
    {
        add_vectors_f(&begin_pt, begin_pt, T);
        add_vectors_f(&end_pt, end_pt, T);
        update_growth_dir_and_length_from_pts();
        update_centroid();

        if (parent)
        {
            copy_vector_f(&(parent->end_pt), begin_pt);
            parent->update_growth_dir_and_length_from_pts();
            parent->update_centroid();
            parent->update_angles();
            parent->update_R();
            parent->update_laterals();
        }


        update_angles();
        update_R();
        update_laterals();

        if (apical)
        {
            copy_vector_f(&(apical->begin_pt), end_pt);
            apical->update_growth_dir_and_length_from_pts();
            apical->update_centroid();
            apical->update_angles();
            apical->update_R();
            apical->update_apical();
            apical->update_laterals();
            apical->update_ancestry_log_prob();
        }
        else
        {
            update_ancestry_log_prob();
        }

    }
    catch (Arg_error e)
    {
        if (entry < 2)
        {
            set_position(old_centroid->elts[0], old_centroid->elts[1], 
                    old_centroid->elts[2]);
        }
        else
        {
            // TODO the length has been coming up wrong because of precision.
            // Probably could throw something like Result_error exception
            // here.
            entry = 0;
        }
        free_vector_f(old_centroid);
        free_vector_f(T);
        throw e;
    }

    free_vector_f(old_centroid);
    free_vector_f(T);

    entry = 0;
}


Apical_structure* Apical_structure::split_into_apical
(
    const Apical_structure* rvals_1,
    const Apical_structure* rvals_2
) 
throw (Arg_error)
{
    Apical_structure* split = clone();
    split->id = generate_id();
    split->base_theta = rvals_2->get_base_theta();
    split->base_psi = rvals_2->get_base_psi();

    if (split->apical)
    {
        split->apical->parent = split;
    }
    apical = split;
    split->parent = this;

    // Back-up some of the structure in case of failure.
    float old_theta  = theta;
    float old_psi    = psi;
    float old_length = length;
    float old_width  = width;
    float old_opacity = opacity;

    Vector_f* old_growth_dir = 0;
    copy_vector_f(&old_growth_dir, growth_dir);

    list<Lateral_structure*> old_laterals = laterals;

    list<float> old_lat_dists;
    list<Lateral_structure*>::iterator lat_it;
    for (lat_it = laterals.begin(); lat_it != laterals.end(); lat_it++)
    {
        old_lat_dists.push_back((*lat_it)->lat_dist);
    }

    // Try the split.
    try
    {
        theta  = rvals_1->get_theta();
        psi    = rvals_1->get_psi();
        width  = rvals_1->get_width();
        length = 0.5f*(length + 
                 (rvals_1->get_length() - density->get_length().get_mu()));
        opacity = rvals_1->get_opacity();

        check_length();

        update_R();
        update_growth_dir_from_R();
        update_centroid();
        update_end_pt();

        split->width = rvals_2->get_width();
        split->opacity = rvals_2->get_opacity();

        float t;
        calc_vector_dot_prod_f(&t, old_growth_dir, growth_dir);
        t *= length / old_length;

        laterals.clear();
        split->laterals.clear();

        for (lat_it = old_laterals.begin(); lat_it != old_laterals.end(); 
                lat_it++)
        {
            if ((*lat_it)->lat_dist < t)
            {
                (*lat_it)->lat_dist /= t;
                laterals.push_back((*lat_it));
            }
            else
            {
                (*lat_it)->lat_dist = ((*lat_it)->lat_dist - t) / (1 - t);
                (*lat_it)->parent = split;
                split->laterals.push_back((*lat_it));
            }

            (*lat_it)->check_lat_dist();
        }

        update_apical();
        update_laterals();

        update_ancestry_log_prob();
    }
    catch (Arg_error e)
    {
        if ((apical = split->apical))
        {
            apical->parent = this;
        }

        split->apical = 0;
        split->laterals.clear();
        delete split;

        theta    = old_theta;
        psi      = old_psi;
        length   = old_length;
        width    = old_width;
        opacity  = old_opacity;
        laterals = old_laterals;

        free_vector_f(old_growth_dir);

        list<float>::iterator dists_it = old_lat_dists.begin();
        for (lat_it = laterals.begin(); lat_it != laterals.end(); lat_it++)
        {
            (*lat_it)->lat_dist = *dists_it;
            (*lat_it)->parent = this;
            dists_it++;
        }

        update_R();
        update_growth_dir_from_R();
        update_centroid();
        update_end_pt();
        update_apical();
        update_laterals();

        throw e;
    }

    free_vector_f(old_growth_dir);

    return split;
}


Apical_structure* Apical_structure::merge_with_apical
(
    const Apical_structure* rvals
) 
throw (Arg_error)
{
    Apical_structure* merged = apical;

    if (!merged)
    {
        throw Arg_error("No apical to merge with");
    }

    list<Lateral_structure*>::iterator it;

    // Back-up some of the structures in case of failure.
    Apical_structure* this_clone = this->clone();
    Apical_structure* merged_clone = merged->clone();

    list<float> this_lat_dists;
    for (it = laterals.begin(); it != laterals.end(); it++)
    {
        this_lat_dists.push_back((*it)->lat_dist);
    }
    list<float> merged_lat_dists;
    for (it = merged->laterals.begin(); it != merged->laterals.end(); it++)
    {
        merged_lat_dists.push_back((*it)->lat_dist);
    }

    // Try the merge.
    try
    {
        width  = rvals->get_width();
        opacity = rvals->get_opacity();

        copy_vector_f(&end_pt, merged->end_pt);

        update_growth_dir_and_length_from_pts();
        update_angles();
        update_R();
        update_centroid();

        if ((apical = merged->apical))
        {
            apical->parent = this;
        }
        update_apical();

        float this_dp;
        calc_vector_dot_prod_f(&this_dp, this_clone->growth_dir, growth_dir);

        float merged_dp;
        calc_vector_dot_prod_f(&merged_dp, merged_clone->growth_dir,growth_dir);

        // update lat_dist
        for (it = laterals.begin(); it != laterals.end(); it++)
        {
            (*it)->lat_dist *= this_dp * this_clone->length / length;

            (*it)->check_lat_dist();
        }
        for (it = merged->laterals.begin(); it != merged->laterals.end(); it++)
        {
            (*it)->lat_dist = 
                ((*it)->lat_dist * merged_dp * merged_clone->length
                         + this_dp * this_clone->length) / length;

            (*it)->check_lat_dist();

            (*it)->parent = this;
            laterals.push_back(*it);
        }
        update_laterals();

        update_ancestry_log_prob();
    }
    catch (Arg_error e)
    {
        // Undo the changes to the merged Apical_structure.
        *merged = *merged_clone;
        if (merged->apical)
        {
            merged->apical->parent = merged;
        }
        merged->update_apical();

        list<float>::iterator lat_dist_it = merged_lat_dists.begin(); 

        for (it = merged->laterals.begin(); it != merged->laterals.end(); it++)
        {
            (*it)->parent = merged;
            (*it)->lat_dist = *lat_dist_it;
            lat_dist_it++;
        }
        merged->update_laterals();

        // Undo effects to this Apical_structure.
        *this = *this_clone;

        lat_dist_it = this_lat_dists.begin(); 

        for (it = laterals.begin(); it != laterals.end(); it++)
        {
            (*it)->parent = this;
            (*it)->lat_dist = *lat_dist_it;
            lat_dist_it++;
        }
        update_laterals();

        this_clone->apical = 0;
        this_clone->laterals.clear();
        merged_clone->apical = 0;
        merged_clone->laterals.clear();
        delete this_clone;
        delete merged_clone;

        throw e;
    }

    // Clear the laterals and the apical so the merged apical can be deleted.
    merged->apical = 0;
    merged->laterals.clear();

    // Clear the laterals and the apical so the backups can be deleted.
    this_clone->apical = 0;
    this_clone->laterals.clear();
    merged_clone->apical = 0;
    merged_clone->laterals.clear();
    delete this_clone;
    delete merged_clone;

    return merged;
}


Lateral_structure* Apical_structure::split_into_lateral
(
    const Apical_structure*  rvals_1,
    const Lateral_structure* rvals_2
) 
throw (Arg_error)
{
    Lateral_structure* split = rvals_2->clone();
    split->parent = this;
    split->level = level;
    split->id = generate_id();
    split->apical = apical;
    copy_vector_f(&(split->end_pt), end_pt);

    if (apical)
    {
        apical->parent = split;
        apical = 0;
    }

    // Back-up some of the structure in case of failure.
    float old_theta  = theta;
    float old_psi    = psi;
    float old_length = length;
    float old_width  = width;

    Vector_f* old_growth_dir = 0;
    copy_vector_f(&old_growth_dir, growth_dir);

    list<Lateral_structure*> old_laterals = laterals;

    list<float> old_lat_dists;
    list<Lateral_structure*>::iterator lat_it;
    for (lat_it = laterals.begin(); lat_it != laterals.end(); lat_it++)
    {
        old_lat_dists.push_back((*lat_it)->lat_dist);
    }

    // Try the split.
    try
    {
        theta  = rvals_1->get_theta();
        psi    = rvals_1->get_psi();
        width  = rvals_1->get_width();
        length = 0.5f*(length + 
                 (rvals_1->get_length() - density->get_length().get_mu()));

        check_length();

        update_R();
        update_growth_dir_from_R();
        update_centroid();
        update_end_pt();

        float dp;
        calc_vector_dot_prod_f(&dp, old_growth_dir, growth_dir);
        float t1 = dp*length / old_length;
        float t2 = dp*split->lat_dist*length / old_length;

        laterals.clear();
        split->laterals.clear();

        for (lat_it = old_laterals.begin(); lat_it != old_laterals.end(); 
                lat_it++)
        {
            if ((*lat_it)->lat_dist < t1)
            {
                (*lat_it)->lat_dist /= t1;
                laterals.push_back((*lat_it));
            }
            else
            {
                (*lat_it)->lat_dist = ((*lat_it)->lat_dist - t2) / (1 - t2);
                (*lat_it)->parent = split;
                split->laterals.push_back((*lat_it));
            }

            (*lat_it)->check_lat_dist();
        }

        laterals.push_back(split);
        update_laterals();

        split->recursively_update_levels(1);

        update_ancestry_log_prob();
    }
    catch (Arg_error e)
    {
        if ((apical = split->apical))
        {
            apical->parent = this;
        }

        split->apical = 0;
        split->laterals.clear();
        delete split;

        theta    = old_theta;
        psi      = old_psi;
        length   = old_length;
        width    = old_width;
        laterals = old_laterals;

        free_vector_f(old_growth_dir);

        list<float>::iterator dists_it = old_lat_dists.begin();
        for (lat_it = laterals.begin(); lat_it != laterals.end(); lat_it++)
        {
            (*lat_it)->lat_dist = *dists_it;
            (*lat_it)->parent = this;
            dists_it++;
        }

        update_R();
        update_growth_dir_from_R();
        update_centroid();
        update_end_pt();
        update_apical();
        update_laterals();

        throw e;
    }

    free_vector_f(old_growth_dir);

    return split;
}


Lateral_structure* Apical_structure::merge_with_lateral
(
    Lateral_structure*      lateral,
    const Apical_structure* rvals
) 
throw (Arg_error)
{
    list<Lateral_structure*>::iterator it;
    bool found = 0;

    for (it = laterals.begin(); !found && it != laterals.end(); it++)
    {
        if (*it == lateral)
            found = 1;
    }

    if (!found)
    {
        throw Arg_error("Lateral to merge with not found in apical");
    }

    Lateral_structure* merged = lateral;

    // Back-up some of the structures in case of failure.
    Apical_structure* this_clone = this->clone();
    Lateral_structure* merged_clone = merged->clone();

    list<float> this_lat_dists;
    for (it = laterals.begin(); it != laterals.end(); it++)
    {
        this_lat_dists.push_back((*it)->lat_dist);
    }
    list<float> merged_lat_dists;
    for (it = merged->laterals.begin(); it != merged->laterals.end(); it++)
    {
        merged_lat_dists.push_back((*it)->lat_dist);
    }

    // Remove merged from this laterals.
    for (it = laterals.begin(); it != laterals.end(); it++)
    {
        if ((*it) == merged)
        {
            laterals.erase(it);
            break;
        }
    }

    // Try the merge.
    try
    {
        width  = rvals->get_width();

        copy_vector_f(&end_pt, merged->end_pt);

        update_growth_dir_and_length_from_pts();
        update_angles();
        update_R();
        update_centroid();

        if ((apical = merged->apical))
        {
            apical->parent = this;
        }
        update_apical();

        float this_dp;
        calc_vector_dot_prod_f(&this_dp, this_clone->growth_dir, growth_dir);

        float merged_dp;
        calc_vector_dot_prod_f(&merged_dp, merged_clone->growth_dir,growth_dir);

        // update lat_dist
        for (it = laterals.begin(); it != laterals.end(); it++)
        {
            (*it)->lat_dist *= this_dp * this_clone->length / length;

            (*it)->check_lat_dist();
        }
        for (it = merged->laterals.begin(); it != merged->laterals.end(); it++)
        {
            (*it)->lat_dist = 
                ((*it)->lat_dist * merged_dp * merged_clone->length
                       + this_dp * merged_clone->lat_dist * this_clone->length) 
                / length;

            (*it)->check_lat_dist();

            (*it)->parent = this;
            laterals.push_back(*it);
        }
        update_laterals();

        merged->recursively_update_levels(-1);

        update_ancestry_log_prob();
    }
    catch (Arg_error e)
    {
        // Undo the changes to the merged Lateral_structure.
        *merged = *merged_clone;
        if (merged->apical)
        {
            merged->apical->parent = merged;
        }
        merged->update_apical();

        list<float>::iterator lat_dist_it = merged_lat_dists.begin(); 

        for (it = merged->laterals.begin(); it != merged->laterals.end(); it++)
        {
            (*it)->parent = merged;
            (*it)->lat_dist = *lat_dist_it;
            lat_dist_it++;
        }
        merged->update_laterals();

        // Undo effects to this Apical_structure.
        *this = *this_clone;

        lat_dist_it = this_lat_dists.begin(); 

        for (it = laterals.begin(); it != laterals.end(); it++)
        {
            (*it)->parent = this;
            (*it)->lat_dist = *lat_dist_it;
            lat_dist_it++;
        }
        update_laterals();

        this_clone->apical = 0;
        this_clone->laterals.clear();
        merged_clone->apical = 0;
        merged_clone->laterals.clear();
        delete this_clone;
        delete merged_clone;

        throw e;
    }

    // Clear the laterals and the apical so the merged apical can be deleted.
    merged->apical = 0;
    merged->laterals.clear();

    // Clear the laterals and the apical so the backups can be deleted.
    this_clone->apical = 0;
    this_clone->laterals.clear();
    merged_clone->apical = 0;
    merged_clone->laterals.clear();
    delete this_clone;
    delete merged_clone;

    return merged;
}



void Apical_structure::replace(Apical_structure* s) throw (Arg_error)
{
    static int entry = 0;

    entry++;

    try
    {
        s->apical = apical;
        if (apical)
        {
            apical->parent = s;
            apical = 0;
        }

        s->laterals.clear();
        for (list<Lateral_structure*>::iterator it = laterals.begin(); 
                it != laterals.end(); it++)
        {
            s->laterals.push_back(*it);
            (*it)->parent = s;
        }
        laterals.clear();

        s->parent = parent;
        if (parent)
        {
            parent->apical = s;
            parent = 0;
        }

        s->update_apical();
        s->update_laterals();
        s->update_ancestry_log_prob();
    }
    catch (Arg_error e)
    {
        if (entry < 2)
        {
            s->replace(this);
        }
        else
        {
            // TODO 
            // Probably could throw something like Result_error exception
            // here.
            entry = 0;
        }
        throw e;
    }

    entry = 0;
}



/*
 * @throw  jwscxx::base::Arg_error  Apical_structure values are not in the 
 *                                  density parameter range.
 */
void Apical_structure::update_from_parent_end_pt() throw (Arg_error)
{
    assert(parent);

    copy_vector_f(&begin_pt, parent->end_pt);

    update_growth_dir_and_length_from_pts();
    update_angles();
    update_R();
    update_centroid();
    if (apical)
    {
        apical->update_from_parent_end_pt();
    }
    update_laterals();

    update_log_prob();
}


void Apical_structure::update_log_prob()
{
    float min, max;

    Structure::update_log_prob();

    /* Centroid_x */
    min       = density->get_centroid_x().get_min();
    max       = density->get_centroid_x().get_max();
    log_prob -= log(max - min);

    /* Centroid_y */
    min       = density->get_centroid_y().get_min();
    max       = density->get_centroid_y().get_max();
    log_prob -= log(max - min);

    /* Centroid_x */
    min       = density->get_centroid_z().get_min();
    max       = density->get_centroid_z().get_max();
    log_prob -= log(max - min);

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


void Apical_structure::check_centroid() const throw (Arg_error)
{
    if (density->get_centroid_x().get_min() > centroid->elts[0] || 
        density->get_centroid_x().get_max() < centroid->elts[0])
    {
        ostringstream ost;
        ost << "Centroid-x " << centroid->elts[0] << " outside density range ["
            << density->get_centroid_x().get_min() << ','
            << density->get_centroid_x().get_max() << ']';
        throw Arg_error(ost.str());
    }

    if (density->get_centroid_y().get_min() > centroid->elts[1] || 
        density->get_centroid_y().get_max() < centroid->elts[1])
    {
        ostringstream ost;
        ost << "Centroid-y " << centroid->elts[1] << " outside density range ["
            << density->get_centroid_y().get_min() << ','
            << density->get_centroid_y().get_max() << ']';
        throw Arg_error(ost.str());
    }

    if (density->get_centroid_z().get_min() > centroid->elts[2] || 
        density->get_centroid_z().get_max() < centroid->elts[2])
    {
        ostringstream ost;
        ost << "Centroid-z " << centroid->elts[2] << " outside density range ["
            << density->get_centroid_z().get_min() << ','
            << density->get_centroid_z().get_max() << ']';
        throw Arg_error(ost.str());
    }
}

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










/* --------------------  Printed_apical_structure  -------------------------- */

Apical_structure* Printed_apical_structure::recursively_convert
(
    Structure* parent,
    std::list<class Printed_apical_structure*>::const_iterator apical_begin, 
    std::list<class Printed_apical_structure*>::const_iterator apical_end, 
    std::list<class Printed_lateral_structure*>::const_iterator lateral_begin, 
    std::list<class Printed_lateral_structure*>::const_iterator lateral_end
) 
const throw (Arg_error)
{
    Apical_structure* s = convert(parent);

    try
    {
        if (apical)
        {
            bool found = false;
            for (list<Printed_apical_structure*>::const_iterator 
                    it = apical_begin; !found && it != apical_end; it++)
            {
                if ((*it)->get_id() == apical)
                {
                    found = true;
                    (*it)->recursively_convert(s, apical_begin, apical_end, 
                                               lateral_begin, lateral_end);
                }
            }
            if (!found)
            {
                ostringstream ost;
                ost << "Failed converting printed spore " << id 
                    << ": Missing printed apical structure " << apical;
                throw Arg_error(ost.str());
            }
        }

        for (list<uint32_t>::const_iterator lateral = laterals.begin(); 
                lateral != laterals.end(); lateral++)
        {
            bool found = false;
            for (list<Printed_lateral_structure*>::const_iterator 
                    it = lateral_begin; !found && it != lateral_end; it++)
            {
                if ((*it)->get_id() == *lateral)
                {
                    found = true;
                    (*it)->recursively_convert(s, apical_begin, apical_end, 
                                               lateral_begin, lateral_end);
                }
            }
            if (!found)
            {
                ostringstream ost;
                ost << "Failed converting printed spore " << id 
                    << ": Missing printed lateral structure " << apical;
                throw Arg_error(ost.str());
            }
        }
    }
    catch (Arg_error e)
    {
        if (parent)
        {
            parent->set_apical(0);
        }
        delete s;
        throw e;
    }

    return s;
}

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










/* ----------------------------  Lateral_structure  ------------------------- */

Lateral_structure::Lateral_structure
(
    float centroid_x, 
    float centroid_y,
    float centroid_z,
    float lat_dist,
    float length,
    float width,
    float base_theta,
    float base_psi,
    float theta,
    float psi,
    float opacity,
    size_t level,
    const class Lateral_structure_density* density
) 
throw (Arg_error)
    : Structure(length, width, base_theta, base_psi, theta, psi, opacity, 
            level, density)
{
    this->lat_dist = lat_dist;
    this->density  = density;

    create_vector_f(&centroid, 3);
    centroid->elts[0] = centroid_x;
    centroid->elts[1] = centroid_y;
    centroid->elts[2] = centroid_z;

    check_centroid();
    check_lat_dist();

    update_R();
    update_growth_dir_from_R();
    update_begin_pt();
    update_end_pt();
}


Lateral_structure::Lateral_structure
(
    Structure* parent,
    float      lat_dist,
    float      length,
    float      width,
    float      base_theta,
    float      base_psi,
    float      theta,
    float      psi,
    float      opacity,
    const class Lateral_structure_density* density
)
throw (Arg_error)
    : Structure(length, width, base_theta, base_psi, theta, psi, opacity, 
            parent->level+1, density)
{
    this->parent   = parent;
    this->lat_dist = lat_dist;
    this->density  = density;

    check_lat_dist();

    multiply_vector_by_scalar_f(&begin_pt, parent->growth_dir, 
            parent->length*lat_dist);
    add_vectors_f(&begin_pt, begin_pt, parent->begin_pt);

    update_R();
    update_growth_dir_from_R();
    update_centroid();
    update_end_pt();

    parent->laterals.push_back(this);
}


Lateral_structure::Lateral_structure(const Lateral_structure& s) : Structure(s)
{
    density  = s.density;
    lat_dist = s.lat_dist;
}


Lateral_structure* Lateral_structure::recursively_clone() const
{
    Lateral_structure* s = clone();

    s->apical = 0;
    s->laterals.clear();

    if (apical)
    {
        s->apical = apical->recursively_clone();
        s->apical->parent = s;
    }
    for (list<Lateral_structure*>::const_iterator it = laterals.begin(); 
            it != laterals.end(); it++)
    {
        Lateral_structure* lateral = (*it)->recursively_clone();
        s->laterals.push_back(lateral);
        lateral->parent = s;
    }

    return s;
}


Lateral_structure& Lateral_structure::operator= (const Lateral_structure& s)
{
    Structure::operator=(s);

    density  = s.density;
    lat_dist = s.lat_dist;

    return *this;
}


void Lateral_structure::set_parent(Structure* parent, float lat_dist) 
    throw (Arg_error)
{
    static int entry = 0;

    entry++;

    Structure* old_parent = this->parent;
    float old_lat_dist    = this->lat_dist;

    try
    {
        if (old_parent)
        {
            old_parent->remove_lateral(this);
        }
    }
    catch (Arg_error e)
    {
        std::cerr << "BUG: Lateral_structure::set_parent: " << e.get_msg() 
                  << '\n';
        assert(0);
    }

    try
    {
        this->parent   = parent;
        this->lat_dist = lat_dist;

        if (!parent)
        {
            check_lat_dist();
            update_ancestry_log_prob();
            return;
        }

        parent->add_lateral(this);

        check_lat_dist();

        multiply_vector_by_scalar_f(&begin_pt, parent->growth_dir, 
                parent->length*lat_dist);
        add_vectors_f(&begin_pt, begin_pt, parent->begin_pt);

        update_growth_dir_and_length_from_pts();
        update_angles();
        update_R();
        update_centroid();
        if (apical)
        {
            apical->update_from_parent_end_pt();
        }
        update_laterals();
        update_ancestry_log_prob();
    }
    catch (Arg_error e)
    {
        if (entry < 2)
        {
            set_parent(old_parent, old_lat_dist);
        }
        else
        {
            // TODO the length has been coming up wrong because of precision.
            // Probably could throw something like Result_error exception
            // here.
            entry = 0;
        }
        throw e;
    }

    entry = 0;
}


void Lateral_structure::set_lat_dist(float lat_dist) throw (Arg_error)
{
    static int entry = 0;

    entry++;

    float old_lat_dist = this->lat_dist;

    try
    {
        this->lat_dist = lat_dist;

        check_lat_dist();

        if (!parent)
            return;

        multiply_vector_by_scalar_f(&begin_pt, parent->growth_dir, 
                parent->length*lat_dist);
        add_vectors_f(&begin_pt, begin_pt, parent->begin_pt);

        update_growth_dir_and_length_from_pts();
        update_angles();
        update_R();
        update_centroid();
        if (apical)
        {
            apical->update_from_parent_end_pt();
        }
        update_laterals();
        update_ancestry_log_prob();
    }
    catch (Arg_error e)
    {
        if (entry < 2)
        {
            set_lat_dist(old_lat_dist);
        }
        else
        {
            // TODO the length has been coming up wrong because of precision.
            // Probably could throw something like Result_error exception
            // here.
            entry = 0;
        }
        throw e;
    }

    entry = 0;
}


void Lateral_structure::print(ostream& out) const
{
    Structure::print(out);

    out << "   Lat Dist: " << lat_dist << '\n';
}


/*
 * @throw  jwscxx::base::Arg_error  Apical_structure values are not in the 
 *                                  density parameter range.
 */
void Lateral_structure::update_from_parent_end_pt() throw (Arg_error)
{
    assert(parent);

    multiply_vector_by_scalar_f(&begin_pt, parent->growth_dir, 
            parent->length*lat_dist);
    add_vectors_f(&begin_pt, begin_pt, parent->begin_pt);

    update_growth_dir_and_length_from_pts();
    update_angles();
    update_R();
    update_centroid();
    if (apical)
    {
        apical->update_from_parent_end_pt();
    }
    update_laterals();

    update_log_prob();
}


void Lateral_structure::update_log_prob()
{
    float    mu, sigma;
    float    min, max;
    double   delta;
    double   p_1, p_2;

    Structure::update_log_prob();

    /* Centroid_x */
    min       = density->get_centroid_x().get_min();
    max       = density->get_centroid_x().get_max();
    log_prob -= log(max - min);

    /* Centroid_y */
    min       = density->get_centroid_y().get_min();
    max       = density->get_centroid_y().get_max();
    log_prob -= log(max - min);

    /* Centroid_x */
    min       = density->get_centroid_z().get_min();
    max       = density->get_centroid_z().get_max();
    log_prob -= log(max - min);

    /* Lateral distance */
    mu        = density->get_lat_dist().get_mu();
    sigma     = density->get_lat_dist().get_sigma();
    min       = density->get_lat_dist().get_min();
    max       = density->get_lat_dist().get_max();
    delta     = (max - min) / 1000;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, lat_dist, lat_dist+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
}


void Lateral_structure::check_centroid() const throw (Arg_error)
{
    if (density->get_centroid_x().get_min() > centroid->elts[0] || 
        density->get_centroid_x().get_max() < centroid->elts[0])
    {
        ostringstream ost;
        ost << "Centroid-x " << centroid->elts[0] << " outside density range ["
            << density->get_centroid_x().get_min() << ','
            << density->get_centroid_x().get_max() << ']';
        throw Arg_error(ost.str());
    }

    if (density->get_centroid_y().get_min() > centroid->elts[1] || 
        density->get_centroid_y().get_max() < centroid->elts[1])
    {
        ostringstream ost;
        ost << "Centroid-y " << centroid->elts[1] << " outside density range ["
            << density->get_centroid_y().get_min() << ','
            << density->get_centroid_y().get_max() << ']';
        throw Arg_error(ost.str());
    }

    if (density->get_centroid_z().get_min() > centroid->elts[2] || 
        density->get_centroid_z().get_max() < centroid->elts[2])
    {
        ostringstream ost;
        ost << "Centroid-z " << centroid->elts[2] << " outside density range ["
            << density->get_centroid_z().get_min() << ','
            << density->get_centroid_z().get_max() << ']';
        throw Arg_error(ost.str());
    }
}


void Lateral_structure::check_lat_dist() const throw (Arg_error)
{
    if (density->get_lat_dist().get_min() > lat_dist || 
        density->get_lat_dist().get_max() < lat_dist)
    {
        ostringstream ost;
        ost << "Lateral distance " << lat_dist << " outside density range ["
            << density->get_lat_dist().get_min() << ','
            << density->get_lat_dist().get_max() << ']';
        throw Arg_error(ost.str());
    }
}

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







/* --------------------  Printed_lateral_structure  ------------------------- */

void Printed_lateral_structure::read_lat_dist(std::istream& in) 
    throw (IO_error, Arg_error)
{
    const char* member_value;

    if (!(member_value = Structure::get_member_value("Lat Dist", in)))
    {
        throw Arg_error("Missing lateral distance");
    }
    istringstream ist(member_value);
    ist >> lat_dist;
    if (ist.fail())
    {
        throw Arg_error("Invalid lateral distance");
    }
}


Lateral_structure* Printed_lateral_structure::recursively_convert
(
    Structure* parent,
    std::list<class Printed_apical_structure*>::const_iterator apical_begin, 
    std::list<class Printed_apical_structure*>::const_iterator apical_end, 
    std::list<class Printed_lateral_structure*>::const_iterator lateral_begin, 
    std::list<class Printed_lateral_structure*>::const_iterator lateral_end
) 
const throw (Arg_error)
{
    Lateral_structure* s = convert(parent);

    try
    {
        if (apical)
        {
            bool found = false;
            for (list<Printed_apical_structure*>::const_iterator 
                    it = apical_begin; !found && it != apical_end; it++)
            {
                if ((*it)->get_id() == apical)
                {
                    found = true;
                    (*it)->recursively_convert(s, apical_begin, apical_end, 
                                               lateral_begin, lateral_end);
                }
            }
            if (!found)
            {
                ostringstream ost;
                ost << "Failed converting printed spore " << id 
                    << ": Missing printed apical structure " << apical;
                throw Arg_error(ost.str());
            }
        }

        for (list<uint32_t>::const_iterator lateral = laterals.begin(); 
                lateral != laterals.end(); lateral++)
        {
            bool found = false;
            for (list<Printed_lateral_structure*>::const_iterator 
                    it = lateral_begin; !found && it != lateral_end; it++)
            {
                if ((*it)->get_id() == *lateral)
                {
                    found = true;
                    (*it)->recursively_convert(s, apical_begin, apical_end, 
                                               lateral_begin, lateral_end);
                }
            }
            if (!found)
            {
                ostringstream ost;
                ost << "Failed converting printed spore " << id 
                    << ": Missing printed lateral structure " << apical;
                throw Arg_error(ost.str());
            }
        }
    }
    catch (Arg_error e)
    {
        if (parent)
        {
            parent->remove_lateral(s);
        }
        delete s;
        throw e;
    }

    return s;
}

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