psf_moves.cpp

Go to the documentation of this file.

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


#include <config.h>

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

#include <inttypes.h>

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

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

#include "sampler.h"
#include "psf_model.h"
#include "psf_moves.h"


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


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

    PSF_model_density* density = params->psf->get_density()->clone();

    float mu    = params->psf->get_alpha();
    float sigma = density->get_alpha().get_sigma();
    float min   = density->get_alpha().get_min();
    float max   = density->get_alpha().get_max();
    density->set_alpha(mu, sigma, min, max);

    mu    = params->psf->get_beta();
    sigma = density->get_beta().get_sigma();
    min   = density->get_beta().get_min();
    max   = density->get_beta().get_max();
    density->set_beta(mu, sigma, min, max);

    mu    = params->psf->get_gamma();
    sigma = density->get_gamma().get_sigma();
    min   = density->get_gamma().get_min();
    max   = density->get_gamma().get_max();
    density->set_gamma(mu, sigma, min, max);

    params->psf_proposal->set(density->sample_alpha(), density->sample_beta(),
            density->sample_gamma(), params->imaging_proposal);

    delete density;

    params->psf_proposal->update_h(params->imaging_proposal);

    params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal);

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

    double log_alpha = params->ll_proposal;
    log_alpha -= params->ll;
    log_alpha += params->psf_proposal->get_log_prob();
    log_alpha -= params->psf->get_log_prob();

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

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

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


const std::ostringstream& PSF_move::get_info(Sampler_move_parameters* params)
{ 
    assert(params->psf);
    info.str("");
    info << std::left << std::setprecision(4) << std::setw(9)
         << params->psf->get_alpha()
         << std::left << std::setprecision(4) << std::setw(9)
         << params->psf->get_beta() 
         << std::left << std::setprecision(4) << std::setw(9)
         << params->psf->get_gamma(); 
    return info;
}


const std::ostringstream& PSF_move::get_proposal_info
(
    Sampler_move_parameters* params
)
{ 
    assert(params->psf_proposal);
    proposal_info.str("");
    proposal_info << std::left << std::setprecision(4) << std::setw(9)
                  << params->psf_proposal->get_alpha()
                  << std::left << std::setprecision(4) << std::setw(9)
                  << params->psf_proposal->get_beta() 
                  << std::left << std::setprecision(4) << std::setw(9)
                  << params->psf_proposal->get_gamma(); 
    return proposal_info;
}


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

    PSF_model_density* density = params->psf->get_density()->clone();
    float mu    = params->psf->get_alpha();
    float sigma = density->get_alpha().get_sigma();
    float min   = density->get_alpha().get_min();
    float max   = density->get_alpha().get_max();
    density->set_alpha(mu, sigma, min, max);

    params->psf_proposal->set_alpha(density->sample_alpha(), 
            params->imaging_proposal);

    delete density;

    params->psf_proposal->update_h(params->imaging_proposal);

    params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal);

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

    double log_alpha = params->ll_proposal;
    log_alpha -= params->ll;
    log_alpha += params->psf_proposal->get_log_prob();
    log_alpha -= params->psf->get_log_prob();

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

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

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


const std::ostringstream& PSF_alpha_move::get_info
(
    Sampler_move_parameters* params
)
{ 
    assert(params->psf);
    info.str("");
    info << std::left << std::setprecision(4) << std::setw(9)
         << params->psf->get_alpha(); 
    return info;
}


const std::ostringstream& PSF_alpha_move::get_proposal_info
(
    Sampler_move_parameters* params
) 
{ 
    assert(params->psf_proposal);
    proposal_info.str("");
    proposal_info << std::left << std::setprecision(4) << std::setw(9)
                  << params->psf_proposal->get_alpha(); 
    return proposal_info;
}


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

    PSF_model_density* density = params->psf->get_density()->clone();
    float mu    = params->psf->get_beta();
    float sigma = density->get_beta().get_sigma();
    float min   = density->get_beta().get_min();
    float max   = density->get_beta().get_max();
    density->set_beta(mu, sigma, min, max);

    params->psf_proposal->set_beta(density->sample_beta(),
            params->imaging_proposal);

    delete density;

    params->psf_proposal->update_h(params->imaging_proposal);

    params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal);

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

    double log_alpha = params->ll_proposal;
    log_alpha -= params->ll;
    log_alpha += params->psf_proposal->get_log_prob();
    log_alpha -= params->psf->get_log_prob();

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

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

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


const std::ostringstream& PSF_beta_move::get_info
(
    Sampler_move_parameters* params
)
{ 
    assert(params->psf);
    info.str("");
    info << std::left << std::setprecision(4) << std::setw(9)
         << params->psf->get_beta(); 
    return info;
}


const std::ostringstream& PSF_beta_move::get_proposal_info
(
    Sampler_move_parameters* params
) 
{ 
    assert(params->psf_proposal);
    proposal_info.str("");
    proposal_info << std::left << std::setprecision(4) << std::setw(9)
                  << params->psf_proposal->get_beta(); 
    return proposal_info;
}


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

    PSF_model_density* density = params->psf->get_density()->clone();
    float mu    = params->psf->get_gamma();
    float sigma = density->get_gamma().get_sigma();
    float min   = density->get_gamma().get_min();
    float max   = density->get_gamma().get_max();
    density->set_gamma(mu, sigma, min, max);

    params->psf_proposal->set_gamma(density->sample_gamma(),
            params->imaging_proposal);

    delete density;

    params->psf_proposal->update_h(params->imaging_proposal);

    params->alternaria_proposal->update_scene(params->psf_proposal, 
            params->imaging_proposal);

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

    double log_alpha = params->ll_proposal;
    log_alpha -= params->ll;
    log_alpha += params->psf_proposal->get_log_prob();
    log_alpha -= params->psf->get_log_prob();

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

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

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


const std::ostringstream& PSF_gamma_move::get_info
(
    Sampler_move_parameters* params
)
{ 
    assert(params->psf);
    info.str("");
    info << std::left << std::setprecision(4) << std::setw(9)
         << params->psf->get_gamma(); 
    return info;
}


const std::ostringstream& PSF_gamma_move::get_proposal_info
(
    Sampler_move_parameters* params
) 
{ 
    assert(params->psf_proposal);
    proposal_info.str("");
    proposal_info << std::left << std::setprecision(4) << std::setw(9)
                  << params->psf_proposal->get_gamma(); 
    return proposal_info;
}