/***************************************************************************
* GModelTemporalLightCurve.cpp - Temporal light curve model class *
* ----------------------------------------------------------------------- *
* copyright (C) 2017 by Juergen Knoedlseder *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see . *
* *
***************************************************************************/
/**
* @file GModelTemporalLightCurve.cpp
* @brief Light curve model class interface implementation
* @author Juergen Knoedlseder
*/
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include
#endif
#include "GTools.hpp"
#include "GException.hpp"
#include "GModelTemporalLightCurve.hpp"
#include "GModelTemporalRegistry.hpp"
#include "GFits.hpp"
/* __ Constants __________________________________________________________ */
/* __ Globals ____________________________________________________________ */
const GModelTemporalLightCurve g_temporal_lightcurve_seed;
const GModelTemporalRegistry g_temporal_lightcurve_registry(&g_temporal_lightcurve_seed);
/* __ Method name definitions ____________________________________________ */
#define G_READ "GModelTemporalLightCurve::read(GXmlElement&)"
#define G_WRITE "GModelTemporalLightCurve::write(GXmlElement&)"
#define G_LOAD_NODES "GModelTemporalLightCurve::load_nodes(GFilename&)"
/* __ Macros _____________________________________________________________ */
/* __ Coding definitions _________________________________________________ */
/* __ Debug definitions __________________________________________________ */
/*==========================================================================
= =
= Constructors/destructors =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Void constructor
***************************************************************************/
GModelTemporalLightCurve::GModelTemporalLightCurve(void) : GModelTemporal()
{
// Initialise members
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief File constructor
*
* @param[in] filename File name of nodes.
* @param[in] norm Normalization factor.
*
* Constructs light curve model from a list of nodes that is found in the
* specified FITS file. See the load_nodes() method for more information
* about the expected structure of the file.
***************************************************************************/
GModelTemporalLightCurve::GModelTemporalLightCurve(const GFilename& filename,
const double& norm) :
GModelTemporal()
{
// Initialise members
init_members();
// Load nodes
load_nodes(filename);
// Set normalization
m_norm.value(norm);
// Return
return;
}
/***********************************************************************//**
* @brief XML constructor
*
* @param[in] xml XML element.
*
* Constructs light curve model by extracting information from an XML
* element. See the read() method for more information about the expected
* structure of the XML element.
***************************************************************************/
GModelTemporalLightCurve::GModelTemporalLightCurve(const GXmlElement& xml) :
GModelTemporal()
{
// Initialise members
init_members();
// Read information from XML element
read(xml);
// Return
return;
}
/***********************************************************************//**
* @brief Copy constructor
*
* @param[in] model Light curve model.
***************************************************************************/
GModelTemporalLightCurve::GModelTemporalLightCurve(const GModelTemporalLightCurve& model) :
GModelTemporal(model)
{
// Initialise members
init_members();
// Copy members
copy_members(model);
// Return
return;
}
/***********************************************************************//**
* @brief Destructor
***************************************************************************/
GModelTemporalLightCurve::~GModelTemporalLightCurve(void)
{
// Free members
free_members();
// Return
return;
}
/*==========================================================================
= =
= Operators =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Assignment operator
*
* @param[in] model Light curve model.
* @return Light curve model.
***************************************************************************/
GModelTemporalLightCurve& GModelTemporalLightCurve::operator=(const GModelTemporalLightCurve& model)
{
// Execute only if object is not identical
if (this != &model) {
// Copy base class members
this->GModelTemporal::operator=(model);
// Free members
free_members();
// Initialise members
init_members();
// Copy members
copy_members(model);
} // endif: object was not identical
// Return
return *this;
}
/*==========================================================================
= =
= Public methods =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Clear light curve model
***************************************************************************/
void GModelTemporalLightCurve::clear(void)
{
// Free class members (base and derived classes, derived class first)
free_members();
this->GModelTemporal::free_members();
// Initialise members
this->GModelTemporal::init_members();
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief Clone light curve model
*
* @return Pointer to deep copy of light curve model.
***************************************************************************/
GModelTemporalLightCurve* GModelTemporalLightCurve::clone(void) const
{
// Clone constant temporal model
return new GModelTemporalLightCurve(*this);
}
/***********************************************************************//**
* @brief Evaluate function
*
* @param[in] srcTime True photon arrival time (not used).
* @param[in] gradients Compute gradients?
* @return Value of light curve model.
*
* Computes
*
* \f[
* S_{\rm t}(t) = r(t) \times {\tt m\_norm}
* \f]
*
* where
* \f$r(t)\f$ is the light curve, computed by linear interpolation between
* the nodes in a FITS file, and \f${\tt m\_norm}\f$ is the normalization
* constant.
*
* If the @p gradients flag is true the method will also evaluate the partial
* derivatives of the model with respect to the normalization parameter using
*
* \f[
* \frac{\delta S_{\rm t}(t)}{\delta {\tt m\_norm}} = r(t)
* \f]
***************************************************************************/
double GModelTemporalLightCurve::eval(const GTime& srcTime,
const bool& gradients) const
{
// Initialise value
double value = 0.0;
// Optionally initialise gradient
if (gradients) {
m_norm.factor_gradient(0.0);
}
// Set value if the time is within the validity range
if (srcTime >= m_tmin && srcTime <= m_tmax) {
// Interpolate file function
double func = m_nodes.interpolate(srcTime.convert(m_timeref), m_values);
// Compute function value
value = m_norm.value() * func;
// Optionally compute gradients
if (gradients) {
// Compute partial derivatives of the parameter values
double g_norm = (m_norm.is_free()) ? m_norm.scale() * func : 0.0;
// Set gradients
m_norm.factor_gradient(g_norm);
} // endif: gradient computation was requested
} // endif: time was within the validity range
// Return
return value;
}
/***********************************************************************//**
* @brief Returns vector of random event times
*
* @param[in] rate Mean event rate (events per second).
* @param[in] tmin Minimum event time.
* @param[in] tmax Maximum event time.
* @param[in,out] ran Random number generator.
*
* This method returns a vector of random event times between @p tmin and
* @p tmax assuming a light curve specified in a FITS file.
***************************************************************************/
GTimes GModelTemporalLightCurve::mc(const double& rate, const GTime& tmin,
const GTime& tmax, GRan& ran) const
{
// Allocates empty vector of times
GTimes times;
// Update Monte Carlo cache
mc_update(tmin, tmax);
// Continue only if effective duration is positive and cache is not empty
if (m_mc_eff_duration > 0.0 && m_mc_cum.size() > 0) {
// Compute mean number of times by multiplying the rate with the
// effective duration and the normalization factor
double lambda = rate * norm() * m_mc_eff_duration;
// Compute number of times to be sampled
int ntimes = int(ran.poisson(lambda)+0.5);
// Loop over number of times
for (int i = 0; i < ntimes; ++i) {
// Determine in which bin we reside
int inx = 0;
if (m_mc_cum.size() > 1) {
double u = ran.uniform();
for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
if (m_mc_cum[inx-1] <= u) {
break;
}
}
}
// Get random time
double seconds;
if (m_mc_slope[inx] == 0.0) {
seconds = m_mc_dt[inx] * ran.uniform() + m_mc_time[inx];
}
else {
seconds = (std::sqrt(m_mc_offset[inx]*m_mc_offset[inx] +
2.0 * m_mc_slope[inx] * ran.uniform()) -
m_mc_offset[inx]) / m_mc_slope[inx] + m_mc_time[inx];
}
// Append random time
GTime time(seconds, m_timeref);
times.append(time);
} // endfor: looped over times
} // endif: cache was valid
// Return vector of times
return times;
}
/***********************************************************************//**
* @brief Read model from XML element
*
* @param[in] xml XML element.
*
* Reads the temporal information from an XML element. The XML element should
* have the format
*
*
*
*
***************************************************************************/
void GModelTemporalLightCurve::read(const GXmlElement& xml)
{
// Get parameter pointers
const GXmlElement* norm = gammalib::xml_get_par(G_READ, xml, m_norm.name());
// Read parameters
m_norm.read(*norm);
// Load nodes from file
load_nodes(gammalib::xml_file_expand(xml, xml.attribute("file")));
// Return
return;
}
/***********************************************************************//**
* @brief Write model into XML element
*
* @param[in] xml XML element.
*
* @exception GExpection::invalid_value
* Invalid XML format encountered.
*
* Writes the temporal information into an XML element in the format
*
*
*
*
***************************************************************************/
void GModelTemporalLightCurve::write(GXmlElement& xml) const
{
// Set model type
if (xml.attribute("type") == "") {
xml.attribute("type", type());
}
// Verify model type
if (xml.attribute("type") != type()) {
std::string msg = "Temporal model of type "+xml.attribute("type")+
" encountered while \""+type()+"\" was expected.";
throw GException::invalid_value(G_WRITE, msg);
}
// Get XML parameters
GXmlElement* norm = gammalib::xml_need_par(G_WRITE, xml, m_norm.name());
// Write parameters
m_norm.write(*norm);
// Set file attribute
xml.attribute("file", gammalib::xml_file_reduce(xml, m_filename));
// Return
return;
}
/***********************************************************************//**
* @brief Print light curve information
*
* @param[in] chatter Chattiness.
* @return String containing light curve information.
***************************************************************************/
std::string GModelTemporalLightCurve::print(const GChatter& chatter) const
{
// Initialise result string
std::string result;
// Continue only if chatter is not silent
if (chatter != SILENT) {
// Append header
result.append("=== GModelTemporalLightCurve ===");
// Append information
result.append("\n"+gammalib::parformat("Light curve file"));
result.append(m_filename.url());
result.append("\n"+gammalib::parformat("Number of parameters"));
result.append(gammalib::str(size()));
for (int i = 0; i < size(); ++i) {
result.append("\n"+m_pars[i]->print(chatter));
}
// Append reference MJD and time range
result.append("\n"+gammalib::parformat("Reference MDJ"));
result.append(gammalib::str(m_timeref.mjdref()));
result.append("\n"+gammalib::parformat("Time range"));
result.append(gammalib::str(m_tmin.convert(m_timeref)));
result.append(" - ");
result.append(gammalib::str(m_tmax.convert(m_timeref)));
result.append(" "+m_timeref.timeunit());
result.append(" ("+m_timeref.timesys()+")");
// Append node information
result.append("\n"+gammalib::parformat("Number of nodes"));
result.append(gammalib::str(m_nodes.size()));
for (int i = 0; i < m_nodes.size(); ++i) {
result.append("\n");
result.append(gammalib::parformat(" Time "+gammalib::str(m_nodes[i])+" s"));
result.append(gammalib::str(m_values[i]));
}
} // endif: chatter was not silent
// Return result
return result;
}
/*==========================================================================
= =
= Private methods =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Initialise class members
***************************************************************************/
void GModelTemporalLightCurve::init_members(void)
{
// Initialise normalisation parameter
m_norm.clear();
m_norm.name("Normalization");
m_norm.unit("(relative value)");
m_norm.scale(1.0);
m_norm.value(1.0);
m_norm.range(0.0,1000.0);
m_norm.fix();
m_norm.gradient(0.0);
m_norm.has_grad(true);
// Set parameter pointer(s)
m_pars.clear();
m_pars.push_back(&m_norm);
// Initialise members
m_nodes.clear();
m_values.clear();
m_filename.clear();
m_timeref.clear();
m_tmin.clear();
m_tmax.clear();
// Initialise cache
m_mc_tmin.clear();
m_mc_tmax.clear();
m_mc_eff_duration = 0.0;
m_mc_cum.clear();
m_mc_slope.clear();
m_mc_offset.clear();
m_mc_time.clear();
m_mc_dt.clear();
// Return
return;
}
/***********************************************************************//**
* @brief Copy class members
*
* @param[in] model Light curve model.
***************************************************************************/
void GModelTemporalLightCurve::copy_members(const GModelTemporalLightCurve& model)
{
// Copy members
m_norm = model.m_norm;
m_nodes = model.m_nodes;
m_values = model.m_values;
m_filename = model.m_filename;
m_timeref = model.m_timeref;
m_tmin = model.m_tmin;
m_tmax = model.m_tmax;
// Copy cache
m_mc_tmin = model.m_mc_tmin;
m_mc_tmax = model.m_mc_tmax;
m_mc_eff_duration = model.m_mc_eff_duration;
m_mc_cum = model.m_mc_cum;
m_mc_slope = model.m_mc_slope;
m_mc_offset = model.m_mc_offset;
m_mc_time = model.m_mc_time;
m_mc_dt = model.m_mc_dt;
// Set parameter pointer(s)
m_pars.clear();
m_pars.push_back(&m_norm);
// Return
return;
}
/***********************************************************************//**
* @brief Delete class members
***************************************************************************/
void GModelTemporalLightCurve::free_members(void)
{
// Return
return;
}
/***********************************************************************//**
* @brief Load nodes from file
*
* @param[in] filename File name.
*
* @exception GException::invalid_value
* File function FITS file is invalid
*
* Load the light curve nodes from a FITS file. The light curve nodes are
* expected in the first extension of the FITS file, containing two mandatory
* columns with names "TIME" and "NORM".
***************************************************************************/
void GModelTemporalLightCurve::load_nodes(const GFilename& filename)
{
// Set maximum light curve normalization value, including a small margin
const double max_norm = 1.0 + 1.0e-8;
// Clear nodes and values
m_nodes.clear();
m_values.clear();
// Set filename
m_filename = filename;
// Load FITS file
GFits fits = GFits(filename);
// Extract binary table (so far always load extension 1 as table)
GFitsTable* table = fits.table(1);
// Read time reference from binary table
m_timeref.read(*table);
// Extract columns
GFitsTableCol* time_col = (*table)["TIME"];
GFitsTableCol* norm_col = (*table)["NORM"];
// Check that there are at least two nodes in table
if (time_col->nrows() < 2) {
std::string msg = "\"TIME\" column contains "+
gammalib::str(time_col->nrows())+" rows but at "
"least two rows are required. Please specify a valid "
"temporal file function.";
throw GException::invalid_value(G_LOAD_NODES, msg);
}
// Check that both columns are consistent
if (time_col->nrows() != norm_col->nrows()) {
std::string msg = "\"TIME\" and \"NORM\" columns have inconsistent "
"number of rows ("+
gammalib::str(time_col->nrows())+", "+
gammalib::str(norm_col->nrows())+"). Please "
"specify a valid temporal file function.";
throw GException::invalid_value(G_LOAD_NODES, msg);
}
// Set number of nodes
int nodes = time_col->nrows();
// Check that time values are in ascending order and that no node is
// larger than 1
double last_time = -1.0;
for (int i = 0; i < nodes; ++i) {
// Check if time has increased
if (last_time >= 0.0 && time_col->real(i) <= last_time) {
std::string msg = "Time "+gammalib::str(time_col->real(i))+
" in row "+gammalib::str(i+1)+" of \"TIME\" "
"column is equal to or smaller than preceeding "
"value. Please provide a light curve file with "
"monotonically increasing times.";
throw GException::invalid_value(G_LOAD_NODES, msg);
}
// Check if value is smaller than maximum allowed normalisation
if (norm_col->real(i) > max_norm) {
std::string msg = "Value "+gammalib::str(norm_col->real(i))+" at "
"time "+gammalib::str(time_col->real(i))+" is "
"larger than 1. Please provide a light curve file "
"with normalizations not exceeding 1.";
throw GException::invalid_value(G_LOAD_NODES, msg);
}
} // endfor: looped over nodes
// Extract nodes
for (int i = 0; i < nodes; ++i) {
m_nodes.append(time_col->real(i));
m_values.push_back(norm_col->real(i));
}
// Make sure that no node exceeds 1
for (int i = 0; i < m_values.size(); ++i) {
if (m_values[i] > 1.0) {
m_values[i] = 1.0;
}
}
// Set minimum and maximum times (assumes that times are ordered)
m_tmin.set(time_col->real(0), m_timeref);
m_tmax.set(time_col->real(time_col->nrows()-1), m_timeref);
// Close FITS file
fits.close();
// Return
return;
}
/***********************************************************************//**
* @brief Set MC pre-computation cache
*
* @param[in] tmin Minimum time of time interval.
* @param[in] tmax Maximum time of time interval.
*
* This method sets up an array of indices and the cumulative distribution
* function needed for MC simulations.
***************************************************************************/
void GModelTemporalLightCurve::mc_update(const GTime& tmin,
const GTime& tmax) const
{
// Update the cache only if the time interval has changed
if (tmin != m_mc_tmin || tmax != m_mc_tmax) {
// Store new time interval
m_mc_tmin = tmin;
m_mc_tmax = tmax;
// Initialise cache
m_mc_eff_duration = 0.0;
m_mc_cum.clear();
m_mc_slope.clear();
m_mc_offset.clear();
m_mc_time.clear();
m_mc_dt.clear();
// Initialise duration
double duration = 0.0;
// Continue only if time interval overlaps with temporal file function
// and if time interval is valid
if (tmax > m_tmin && tmin < m_tmax && tmax > tmin) {
// Loop over all intervals between the nodes
for (int i = 0; i < m_nodes.size()-1; ++i) {
// Get start and stop time of interval
GTime tstart(m_nodes[i], m_timeref);
GTime tstop(m_nodes[i+1], m_timeref);
// Make sure that we only consider the interval that overlaps
// with the requested time interval
if (tmin > tstart) {
tstart = tmin;
}
if (tmax < tstop) {
tstop = tmax;
}
// If time interval is empty then skip this node
if (tstart >= tstop) {
continue;
}
// Evaluate rate at start and stop time
double rstart = eval(tstart);
double rstop = eval(tstop);
// Compute mean normalization for this interval
double norm = 0.5 * (rstart + rstop);
// Compute integral for this interval
double dt = tstop - tstart;
double cum = norm * dt;
// Compute slope, offset and start time of interval
double slope = (rstop - rstart) / dt;
double offset = rstart;
double renorm = (0.5 * slope * dt + offset) * dt;
slope /= renorm;
offset /= renorm;
// Update effective duration and real duration
m_mc_eff_duration += cum;
duration += dt;
// Put mean normalization and integral on stack
m_mc_cum.push_back(cum);
m_mc_slope.push_back(slope);
m_mc_offset.push_back(offset);
m_mc_time.push_back(tstart.convert(m_timeref));
m_mc_dt.push_back(dt);
} // endfor: next interval
// Build cumulative distribution
for (int i = 1; i < m_mc_cum.size(); ++i) {
m_mc_cum[i] += m_mc_cum[i-1];
}
double norm = m_mc_cum[m_mc_cum.size()-1];
for (int i = 0; i < m_mc_cum.size(); ++i) {
m_mc_cum[i] /= norm;
}
} // endif: time interval overlaps
} // endif: Update was required
// Return
return;
}