EstimatePlotter.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2003 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
// psrchive/Util/units/EstimatePlotter.h
 
#ifndef __EstimatePlotter_h
#define __EstimatePlotter_h
 
#include "Estimate.h"
#include "Error.h"
#include "myfinite.h"
 
#include <vector>
#include <iostream>
#include <cmath>
 
class EstimatePlotter {
 
 public:
 
  static bool report_mean;
 
  bool report_mean_on_single_line;
 
  EstimatePlotter ();
 
  void clear ();
 
  void set_control_viewport (bool flag);
 
  void set_border (float fraction_x, float fraction_y);
 
  void set_xrange (float x_min, float x_max);
 
  void set_minimum_error (float error);
 
  void set_maximum_error (float error);
 
  void set_graph_marker (int symbol);
 
  void separate_viewports (bool scaled = true, bool vertical = true);
 
  void set_viewport (unsigned index);
 
  void set_world (float x1, float x2, float y1, float y2);
 
  float get_x_min () const { return x_min; }
  float get_x_max () const { return x_max; }
 
  float get_y_min () const { return y_min; }
  float get_y_max () const { return y_max; }
 
  void restore_viewport ();
 
  template<class T> void add_plot (const std::vector< Estimate<T> >& data);
 
  template<class Xt, class Yt>
  void add_plot (const std::vector<Xt>&, const std::vector< Estimate<Yt> >&);
 
  template<class Xt, class Yt> 
  void add_plot (const std::vector<Xt>& xdata,
                 const std::vector< std::complex< Estimate<Yt> > >& ydata);
 
  template<class T> void plot (const std::vector< Estimate<T> >& data);
 
  unsigned plot (unsigned index);
 
  unsigned excise (unsigned index);
 
  unsigned plot ();
 
  void set_minmax_range (unsigned i_min, unsigned i_max);
 
  void minmax (bool& xrange_set, float& xmin, float& xmax,
               bool& yrange_set, float& ymin, float& ymax,
               const std::vector<float>& x,
               const std::vector<float>& y,
               const std::vector<float>& yerr);
 
  void set_plot_error_bars (bool flag)
  { plot_error_bars = flag; }
  
 protected:
 
  bool control_viewport;
 
  float x_border, y_border;
 
  float xrange_min, xrange_max;
 
  float x_min, x_max;
 
  float y_min, y_max;
 
  unsigned i_min, i_max;
 
  bool xrange_set;
  bool yrange_set;
 
  float minimum_error;
 
  float maximum_error;
 
  int graph_marker;
 
  bool plot_error_bars;
 
  float log_err_madm_threshold;
  bool log_err_madm_applied;
  
  std::vector< std::vector<float> > xval;
  std::vector< std::vector<float> > yval;
  std::vector< std::vector<float> > yerr;
  
 private:
 
  std::vector<float> data_xmin;
  std::vector<float> data_xmax;
  std::vector<float> data_ymin;
  std::vector<float> data_ymax;
 
  float vp_x1, vp_x2, vp_y1, vp_y2;
 
  bool viewports_set;
  bool viewports_vertical;
  bool viewports_scaled;
 
};
 
template<class T> 
void EstimatePlotter::plot (const std::vector< Estimate<T> >& data)
{
  clear ();
  add_plot (data);
  plot (0);
}
 
template<class T> 
void EstimatePlotter::add_plot (const std::vector< Estimate<T> >& data)
{
  unsigned ipt = 0, npt = data.size();
  if (npt == 0)
    return;
 
  log_err_madm_applied = false;
 
  xval.push_back ( std::vector<float>(npt) );
  yval.push_back ( std::vector<float>(npt) );
  yerr.push_back ( std::vector<float>(npt) );
 
  std::vector<float>& x = xval.back();
  std::vector<float>& y = yval.back();
  std::vector<float>& ye = yerr.back();
 
  double xscale = 0.0;
  if (npt > 1)
    xscale = double(xrange_max - xrange_min) / double(npt-1);
 
  MeanEstimate<T> mean;
 
  for (ipt=0; ipt<npt; ipt++)
  {
    x[ipt] = xrange_min + xscale * double(ipt);
 
    if (!myfinite( data[ipt].val ))
      ye[ipt] = 0;
    else 
    {
      y[ipt] = data[ipt].val;
      ye[ipt] = sqrt (data[ipt].var);
 
      if (data[ipt].var)
        mean += data[ipt];
    }
 
  }
 
  // useful at times, excessive at others ...
  if (report_mean)
  {
    std::cerr << "Mean = " << mean << std::endl;
  }
  else if (report_mean_on_single_line)
  {
    std::cout << mean.get_Estimate().get_value() << " ";
    std::cout << mean.get_Estimate().get_error() << " ";
  }
 
  minmax (xrange_set, x_min, x_max, yrange_set, y_min, y_max, x, y, ye);
}
 
template<class Xt, class Yt> 
void EstimatePlotter::add_plot (const std::vector<Xt>& xdata,
                                const std::vector< Estimate<Yt> >& ydata)
{
  unsigned ipt = 0, npt = xdata.size();
  if (npt == 0)
    return;
  
  if (ydata.size() != npt)
    throw Error (InvalidParam, "EstimatePlotter::add_plot (Xt, Yt)",
                 "xdata.size=%d != ydata.size=%d", npt, ydata.size());
 
  log_err_madm_applied = false;
 
  xval.push_back ( std::vector<float>(npt) );
  yval.push_back ( std::vector<float>(npt) );
  yerr.push_back ( std::vector<float>(npt) );
 
  std::vector<float>& x = xval.back();
  std::vector<float>& y = yval.back();
  std::vector<float>& ye = yerr.back();
 
  for (ipt=0; ipt<npt; ipt++) {
    ye[ipt] = sqrt (ydata[ipt].var);
    x[ipt] = xdata[ipt];
    y[ipt] = ydata[ipt].val;
  }
 
  minmax (xrange_set, x_min, x_max, yrange_set, y_min, y_max, x, y, ye);
}
 
template<class Xt, class Yt> 
void EstimatePlotter::add_plot
(const std::vector<Xt>& xdata,
 const std::vector< std::complex< Estimate<Yt> > >& ydata)
{
  unsigned ipt = 0, npt = xdata.size();
  if (npt == 0)
    return;
  
  if (ydata.size() != npt)
    throw Error (InvalidParam, "EstimatePlotter::add_plot (Xt, Yt)",
                 "xdata.size=%d != ydata.size=%d", npt, ydata.size());
 
  log_err_madm_applied = false;
 
  xval.push_back ( std::vector<float>(npt) );
  std::vector<float>& x = xval.back();
 
  std::vector<float> re_y (npt);
  std::vector<float> re_ye (npt);
 
  std::vector<float> im_y (npt);
  std::vector<float> im_ye (npt);
 
  for (ipt=0; ipt<npt; ipt++)
  {
    x[ipt] = xdata[ipt];
 
    re_y[ipt] = ydata[ipt].real().val;
    re_ye[ipt] = sqrt (ydata[ipt].real().var);
 
    im_y[ipt] = ydata[ipt].imag().val;
    im_ye[ipt] = sqrt (ydata[ipt].imag().var);
  }
 
  minmax (xrange_set, x_min, x_max, yrange_set, y_min, y_max, x, re_y, re_ye);
  minmax (xrange_set, x_min, x_max, yrange_set, y_min, y_max, x, im_y, im_ye);
 
  yval.push_back ( re_y );
  yval.push_back ( im_y );
 
  yerr.push_back ( re_ye );
  yerr.push_back ( im_ye );
}
 
#endif