imaging_model.cpp

Go to the documentation of this file.

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


#include <config.h>

#include <cassert>
#include <cmath>
#include <iostream>
#include <fstream>
#include <sstream>

#include <inttypes.h>

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

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

#include "density.h"
#include "printable.h"
#include "imaging_model.h"


#define  WINDOW_COL         0
#define  WINDOW_ROW         0
#define  WINDOW_IMG         0
#define  WINDOW_NUM_COLS    1
#define  WINDOW_NUM_ROWS    1
#define  WINDOW_NUM_IMGS    1

#define  SCALE_X            1.0f
#define  SCALE_Y            1.0f
#define  SCALE_Z            2.0f

#define  BG_MU              0.8f
#define  BG_SIGMA           0.1f
#define  BG_MIN             0
#define  BG_MAX             1.0f

#define  C_2_MIN            2.0e-4f
#define  C_2_MAX            2.0e-4f


using std::isinf;
using std::isnan;
using std::log;
using std::ostream;
using std::istream;
using std::ofstream;
using std::ifstream;
using std::ostringstream;
using std::istringstream;
using namespace jwsc;
using jwscxx::base::IO_error;
using jwscxx::base::Arg_error;


/* --------------------------  Imaging_window  ------------------------------ */

Imaging_window::Imaging_window() throw (Arg_error)
{
    set(WINDOW_COL, WINDOW_ROW, WINDOW_IMG, WINDOW_NUM_COLS, WINDOW_NUM_ROWS,
            WINDOW_NUM_IMGS);
}


void Imaging_window::set
(
    uint32_t col, 
    uint32_t row, 
    uint32_t img, 
    uint32_t num_cols,
    uint32_t num_rows,
    uint32_t num_imgs
)
throw (Arg_error)
{
    if (num_cols == 0 || num_rows == 0 || num_imgs == 0)
    {
        throw Arg_error("Window size is zero");
    }

    this->col      = col;
    this->row      = row;
    this->img      = img;
    this->num_cols = num_cols;
    this->num_rows = num_rows;
    this->num_imgs = num_imgs;
}


ostream& operator<<(ostream& out, const Imaging_window& window)
{
    return out << window.col << ' ' << window.row << ' ' << window.img << ' '
               << window.num_cols << ' ' << window.num_rows << ' '
               << window.num_imgs;
}


istream& operator>>(istream& in, Imaging_window& window)
{
    in >> window.col;
    in >> window.row;
    in >> window.img;
    in >> window.num_cols;
    in >> window.num_rows;
    in >> window.num_imgs;
    return in;
}

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









/* --------------------------  Imaging_scale  ------------------------------- */

Imaging_scale::Imaging_scale() throw (Arg_error)
{
    set(SCALE_X, SCALE_Y, SCALE_Z);
}


void Imaging_scale::set(float x, float y, float z) throw (Arg_error)
{
    if (x <= 0 || y <= 0 || z <= 0)
    {
        throw Arg_error("Imaging scale <= 0");
    }

    this->x = x;
    this->y = y;
    this->z = z;
}


ostream& operator<<
(
     ostream&             out, 
     const Imaging_scale& scale
)
{
    return out << scale.x << ' ' << scale.y << ' ' << scale.z;
}


istream& operator>>
(
     istream&       in, 
     Imaging_scale& scale
)
{
    in >> scale.x;
    in >> scale.y;
    in >> scale.z;
    return in;
}

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









/* --------------------  Imaging_model_density  ----------------------------- */

Imaging_model_density::Imaging_model_density() throw (Arg_error)
{
    float mu, sigma;
    float min, max;

    mu    = BG_MU;
    sigma = BG_SIGMA;
    min   = BG_MIN;
    max   = BG_MAX;
    set_bg(mu, sigma, min, max);

    min   = C_2_MIN;
    max   = C_2_MAX;
    set_c_2(min, max);
}


Imaging_model_density* Imaging_model_density::clone() const
{
    float mu, sigma;
    float min, max;

    Imaging_model_density* d = new Imaging_model_density();

    mu    = bg.get_mu();
    sigma = bg.get_sigma();
    min   = bg.get_min();
    max   = bg.get_max();
    d->set_bg(mu, sigma, min, max);

    min   = c_2.get_min();
    max   = c_2.get_max();
    d->set_c_2(min, max);

    return d;
}


void Imaging_model_density::set_bg(float mu, float sigma, float min, float max) 
    throw (Arg_error)
{
    if (min < 0)
    {
        throw Arg_error("Minimum must be >= 0");
    }
    if (max > 1)
    {
        throw Arg_error("Maximum must be <= 1");
    }
    bg.set(mu, sigma, min, max);
}


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


float Imaging_model_density::sample_bg() const
{
    return (float)sample_gaussian_pdf_d(bg.get_mu(), bg.get_sigma(), 
            bg.get_min(), bg.get_max());
}


float Imaging_model_density::sample_c_2() const
{
    return (float)sample_uniform_pdf_d(c_2.get_min(), c_2.get_max());
}


Imaging_model* Imaging_model_density::sample
(
    const Imaging_window* window,
    const Imaging_scale*  scale
)
const
{
    return new Imaging_model(sample_bg(), sample_c_2(), window, scale, this);
}

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









/* ----------------------------  Imaging_model  ----------------------------- */

Imaging_model::Imaging_model
(
    float                        bg,
    float                        c_2,
    const Imaging_window*        window,
    const Imaging_scale*         scale,
    const Imaging_model_density* density
)
throw (Arg_error)
{
    this->bg      = bg;
    this->c_2     = c_2;
    this->window  = window;
    this->scale   = scale;
    this->density = density;

    check_bg();
    check_c_2();

    update_log_prob();
}


Imaging_model::Imaging_model(const Imaging_model& model)
{
    bg       = model.bg;
    c_2      = model.c_2;
    window   = model.window;
    scale    = model.scale;
    density  = model.density;
    log_prob = model.log_prob;
}


Imaging_model::Imaging_model
(
    const char*                  fname, 
    const Imaging_window*        window,
    const Imaging_scale*         scale,
    const Imaging_model_density* density
)
throw (IO_error, Arg_error)
{
    ostringstream ost;
    istringstream ist;
    ifstream      in;

    this->window  = window;
    this->scale   = scale;
    this->density = density;

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

    const char* member_value;

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

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

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

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

    check_bg();
    check_c_2();

    update_log_prob();
}


Imaging_model* Imaging_model::clone()
{
    return new Imaging_model(*this);
}


void Imaging_model::set_bg(float bg) throw (Arg_error)
{
    float old_bg = this->bg;

    try
    {
        this->bg = bg;
        check_bg();

        update_log_prob();
    }
    catch (Arg_error e)
    {
        set_bg(old_bg);
        throw e;
    }
}


void Imaging_model::set_c_2(float c_2) throw (Arg_error)
{
    float old_c_2 = this->c_2;

    try
    {
        this->c_2 = c_2;
        check_c_2();

        update_log_prob();
    }
    catch (Arg_error e)
    {
        set_c_2(old_c_2);
        throw e;
    }
}


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

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

    print(out);

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


void Imaging_model::print(std::ostream& out) const
{
    out << "       Bg: " << bg << '\n'
        << "      C_2: " << c_2 << '\n'
        << " Log Prob: " << log_prob << '\n';
}


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

    log_prob = 0;

    mu        = density->get_bg().get_mu();
    sigma     = density->get_bg().get_sigma();
    min       = density->get_bg().get_min();
    max       = density->get_bg().get_max();
    delta     = (max - min) / 1000.0;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, bg, bg+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(!isinf(log_prob) && !isnan(log_prob));
#endif
}


void Imaging_model::check_bg() const throw (Arg_error)
{
    if (density->get_bg().get_min() > bg || 
        density->get_bg().get_max() < bg)
    {
        ostringstream ost;
        ost << "Background " << bg << " outside density range ["
            << density->get_bg().get_min() << ','
            << density->get_bg().get_max() << ']';
        throw Arg_error(ost.str());
    }
}


void Imaging_model::check_c_2() const throw (Arg_error)
{
    if (density->get_c_2().get_min() > c_2 || 
        density->get_c_2().get_max() < c_2)
    {
        ostringstream ost;
        ost << "C_2 " << c_2 << " outside density range ["
            << density->get_c_2().get_min() << ','
            << density->get_c_2().get_max() << ']';
        throw Arg_error(ost.str());
    }
}

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