sampler.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 <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <list>

#include <inttypes.h>
#include <sys/time.h>

#include <jwsc/prob/pdf.h>
#include <jwsc/matblock/matblock.h>

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

#include "alternaria_model.h"
#include "psf_model.h"
#include "imaging_model.h"
#include "sampler.h"


using std::list;
using std::ostream;
using std::ofstream;
using std::ostringstream;
using namespace jwsc;
using jwscxx::base::Exception;
using jwscxx::base::Arg_error;
using jwscxx::base::IO_error;


Sampler_move_parameters::Sampler_move_parameters
(
    Alternaria_model*       alternaria,
    PSF_model*              psf,
    Imaging_model*          imaging,
    const jwsc::Matblock_f* data
)
throw (Arg_error)
{
    this->alternaria = alternaria;
    this->psf        = psf;
    this->imaging    = imaging;
    this->data       = data;

    data_avg = 0;

    uint32_t num_elts = data->num_mats*data->num_rows*data->num_cols;
    float*   elts     = **(data->elts);

    for (uint32_t i = 0; i < num_elts; i++)
    {
        data_avg += elts[i];
    }
    data_avg /= (float)num_elts;

    ll = alternaria->calc_log_likelihood(psf, imaging, data, data_avg);

    alternaria_proposal = alternaria->clone();
    psf_proposal        = psf->clone();
    imaging_proposal    = imaging->clone();
    ll_proposal         = ll;
}


Sampler_move_parameters::~Sampler_move_parameters()
{
    delete alternaria;
    delete psf;
    delete imaging;
    delete alternaria_proposal;
    delete psf_proposal;
    delete imaging_proposal;
}


void Sampler_move_parameters::init_proposals()
{
    *alternaria_proposal = *alternaria;
    *psf_proposal        = *psf;
    *imaging_proposal    = *imaging;
    ll_proposal          = ll;
}


Sampler_diffusion_move::Sampler_diffusion_move(const char* name, double prob) 
    throw (Arg_error)
{
    num_accepted = 0;
    num_attempts = 0;
    this->name   = name;
    set_prob(prob);
}


void Sampler_diffusion_move::set_prob(double prob) throw (Arg_error)
{
    if (0 > prob || prob > 1.0)
    {
        throw Arg_error("Sampler move probability not in [0,1]");
    }

    this->prob = prob;
}


Sampler_jump_move::Sampler_jump_move
(
    const char* name_1, 
    const char* name_2,
    double      prob_1,
    double      prob_2
) 
throw (Arg_error)
{
    num_accepted_1 = 0;
    num_attempts_1 = 0;
    this->name_1   = name_1;
    set_prob_1(prob_1);

    num_accepted_2 = 0;
    num_attempts_2 = 0;
    this->name_2   = name_2;
    set_prob_2(prob_2);

    name = name_1;
}


bool Sampler_jump_move::run(Sampler_move_parameters* params) throw (Exception)
{
    if (sample_uniform_pdf_d(0, prob_1+prob_2) <= prob_1)
    {
        name = name_1;
        return run_1(params);
    }
    name = name_2;
    return run_2(params);
}


const std::ostringstream& Sampler_jump_move::get_info
(
    Sampler_move_parameters* params
) 
{
    assert(params->alternaria);
    info.str("");
    info << std::left << std::setw(4) 
         << params->alternaria->get_num_levels();
    info << std::left << std::setw(4) 
         << params->alternaria->get_num_apical_hypha();
    info << std::left << std::setw(4) 
         << params->alternaria->get_num_lateral_hypha();
    info << std::left << std::setw(4) 
         << params->alternaria->get_num_spores();
    return info;
}


const std::ostringstream& Sampler_jump_move::get_proposal_info
(
    Sampler_move_parameters* params
)
{ 
    assert(params->alternaria_proposal);
    proposal_info.str("");
    proposal_info << std::left << std::setw(4)
                  << params->alternaria_proposal->get_num_levels();
    proposal_info << std::left << std::setw(4) 
                  << params->alternaria_proposal->get_num_apical_hypha();
    proposal_info << std::left << std::setw(4)
                  << params->alternaria_proposal->get_num_lateral_hypha();
    proposal_info << std::left << std::setw(4)
                  << params->alternaria_proposal->get_num_spores();
    return proposal_info; 
}


void Sampler_jump_move::set_prob(double prob) throw (Arg_error)
{
    if (0 > prob || prob > 1.0)
    {
        throw Arg_error("Sampler move probability not in [0,1]");
    }

    if ((prob_1 + prob_2) == 0 || prob == 0)
    {
        set_prob_1(0);
        set_prob_2(0);
    }

    double d = (prob_1 + prob_2) / prob;

    set_prob_1(prob_1/d);
    set_prob_2(prob_2/d);
}


void Sampler_jump_move::set_prob_1(double prob_1) throw (Arg_error)
{
    if (0 > prob_1 || prob_1 > 1.0)
    {
        throw Arg_error("Sampler move probability not in [0,1]");
    }

    this->prob_1 = prob_1;
}


void Sampler_jump_move::set_prob_2(double prob_2) throw (Arg_error)
{
    if (0 > prob_2 || prob_2 > 1.0)
    {
        throw Arg_error("Sampler move probability not in [0,1]");
    }

    this->prob_2 = prob_2;
}


Sampler::~Sampler()
{
    for (list<Sampler_move*>::iterator it = moves.begin(); 
            it != moves.end(); it++)
    {
        delete *it;
    }
}


void Sampler::add_move(Sampler_move* move)
{
    moves.push_back(move);
}


void Sampler::normalize_move_probs() throw (Arg_error)
{
    double sum = 0;

    for (list<Sampler_move*>::iterator it = moves.begin(); 
            it != moves.end(); it++)
    {
        sum += (*it)->get_prob();
    }

    if (sum <= 1.0e-10)
    {
        throw Arg_error("Sum of move probabilities is zero");
    }

    for (list<Sampler_move*>::iterator it = moves.begin(); 
            it != moves.end(); it++)
    {
        (*it)->set_prob((*it)->get_prob() / sum);
    }
}


void Sampler::run
(
    Alternaria_model**       best_alternaria,
    PSF_model**              best_psf,
    Imaging_model**          best_imaging,
    Sampler_move_parameters* params,
    uint32_t                 iterations,
    const char*              move_fname,
    const char*              best_fname,
    const char*              alt_pro_fmt
) 
throw (Exception)
{
    ostringstream ost;
    ofstream move_out;
    ofstream best_out;

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

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

    run(best_alternaria, best_psf, best_imaging, params, iterations, 
            move_out, best_out, alt_pro_fmt);

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

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


void Sampler::run
(
    Alternaria_model**       best_alternaria,
    PSF_model**              best_psf,
    Imaging_model**          best_imaging,
    Sampler_move_parameters* params,
    uint32_t                 iterations,
    std::ostream&            move_out,
    std::ostream&            best_out,
    const char*              alt_pro_fmt
) 
throw (Exception)
{
    assert(*best_alternaria == 0);
    assert(*best_psf == 0);
    assert(*best_imaging == 0);

    *best_alternaria = params->alternaria->clone();
    *best_psf        = params->psf->clone();
    *best_imaging    = params->imaging->clone();

    double ll_best = params->ll;
    double lp_best = ll_best;
    lp_best += params->alternaria->get_log_prob();
    lp_best += params->psf->get_log_prob();
    lp_best += params->imaging->get_log_prob();

    double lp = params->ll;
    lp += params->alternaria->get_log_prob();
    lp += params->psf->get_log_prob();
    lp += params->imaging->get_log_prob();

    list<Sampler_move*>::iterator it;
    list<Sampler_move*>::iterator itt;

    // Print some inital state information for the moves.
    best_out << std::left << std::setw(7)  << 0
             << std::left << std::setw(14) << ll_best
             << std::left << std::setw(14) << lp_best;
    for (it = moves.begin(); it != moves.end(); it++)
    {
        Sampler_jump_move* jump_move;

        if ((jump_move = dynamic_cast<Sampler_jump_move*>(*it)))
        {
            move_out << std::left << std::setw(7)  << 0
                     << std::left << std::setw(12)  << jump_move->get_name_1()
                     << std::left << std::setw(8)  << "Init"
                     << std::left << std::setw(14) << params->ll
                     << std::left << std::setw(14) << lp
                     << jump_move->get_info(params).str() << '\n';
            move_out << std::left << std::setw(7)  << 0
                     << std::left << std::setw(12)  << jump_move->get_name_2()
                     << std::left << std::setw(8)  << "Init"
                     << std::left << std::setw(14) << params->ll
                     << std::left << std::setw(14) << lp
                     << jump_move->get_info(params).str() << '\n';
        }
        else
        {
            move_out << std::left << std::setw(7)  << 0
                     << std::left << std::setw(12)  << (*it)->get_name()
                     << std::left << std::setw(8)  << "Init"
                     << std::left << std::setw(14) << params->ll
                     << std::left << std::setw(14) << lp
                     << (*it)->get_info(params).str() << '\n';
        }
        move_out.flush();

        best_out << (*it)->get_info(params).str();
    }
    best_out << '\n';
    best_out.flush();

    struct timeval tv1;
    struct timeval tv2;

    // Iteratively generate model samples, keeping track of the best models.
    for (uint32_t i = 1; i <= iterations; i++)
    {
        double u = sample_uniform_pdf_d(0, 1);
        double p = 0;

        move_out << std::left << std::setw(7) << i;

        // Execute one of the sampler moves.
        for (it = moves.begin(); it != moves.end(); it++)
        {
            Sampler_move* move = *it;

            p += move->get_prob();
            if (u <= p)
            {
                bool accept;

                try
                {
                    assert(gettimeofday(&tv1, NULL) == 0);
                    accept = move->run(params);
                    assert(gettimeofday(&tv2, NULL) == 0);
                }
                catch (Exception e)
                {
                    try
                    {
                        params->alternaria->get_root()->recursively_check_structure();
                    }
                    catch (Exception ee)
                    {
                        // This should never happen!
                        std::cerr << "alternaria: sampler.cpp: " 
                                  << ee.get_msg() << '\n';
                        abort();
                    }

                    try
                    {
                        params->alternaria_proposal->get_root()->recursively_check_structure();
                    }
                    catch (Exception ee)
                    {
                        accept = false;
                    }
                }

                double lp_proposal = params->ll_proposal;
                lp_proposal += params->alternaria_proposal->get_log_prob();
                lp_proposal += params->psf_proposal->get_log_prob();
                lp_proposal += params->imaging_proposal->get_log_prob();

                move_out << std::left << std::setw(12) << move->get_name();

                if (accept)
                {
                    move_out << std::left << std::setw(8) << "Accept";

                    *(params->alternaria) = *(params->alternaria_proposal);
                    *(params->psf)        = *(params->psf_proposal);
                    *(params->imaging)    = *(params->imaging_proposal);

                    params->ll = params->ll_proposal;

                    // Update the best models if necessary.
                    if (lp_proposal > lp_best)
                    {
                        **best_alternaria = *(params->alternaria_proposal);
                        **best_psf        = *(params->psf_proposal);
                        **best_imaging    = *(params->imaging_proposal);
                        ll_best           = params->ll_proposal;
                        lp_best           = lp_proposal;

                        best_out << std::left << std::setw(7) << i
                                 << std::left << std::setw(14) << ll_best
                                 << std::left << std::setw(14) << lp_best;
                        for (itt = moves.begin(); itt != moves.end(); itt++)
                        {
                            best_out << (*itt)->get_info(params).str();
                        }
                        best_out << '\n';
                        best_out.flush();
                    }
                }
                else
                {
                    move_out << std::left << std::setw(8) << "Reject";
                }

                double seconds = (tv2.tv_sec - tv1.tv_sec) + 
                         (tv2.tv_usec - tv1.tv_usec)*1.0e-6;

                move_out << std::left << std::setw(14) << params->ll_proposal 
                         << std::left << std::setw(14) << lp_proposal
                         << std::left << std::setprecision(4) << std::setw(10)
                         << seconds
                         << move->get_proposal_info(params).str() << '\n';
                move_out.flush();

                try
                {
                  params->alternaria->get_root()->recursively_check_structure();
                }
                catch (Exception e)
                {
                    std::cerr << "alternaria: sampler.cpp: " 
                              << e.get_msg() << '\n';
                    abort();
                }

                break;
            }
        }

        if (alt_pro_fmt)
        {
            char buf[256] = {0};
            sprintf(buf, alt_pro_fmt, i);
            params->alternaria_proposal->print(buf);
        }
    }
}