AlgAdapt.cxx

Go to the documentation of this file.

// $Id: AlgAdapt.cxx,v 1.11 2009/10/15 18:30:13 jyuko Exp $

// C/C++
#include <cassert>
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <iostream>

// MINOS
#include "Registry/Registry.h"
#include "Util/UtilString.h"

// Local
#include "PhysicsNtuple/AlgAdapt.h"

using namespace std;

//---------------------------------------------------------------------------------
Anp::AlgAdapt::AlgAdapt()
{
}

//---------------------------------------------------------------------------------
Anp::AlgAdapt::~AlgAdapt()
{
}

//------------------------------------------------------------------------------------------
long double Anp::AlgAdapt::GetPT() const
{
   if(fCurPar.size() != fDefPar.size() || fCurPar.size() != fDefErr.size())
   {
      cerr << "AlgAdapt::GetPT - parameter and/or error vectors have different size" << endl;
      return 0.0;
   }

   long double pterm = 0.0;

   for(unsigned int ipar = 0; ipar < fCurPar.size(); ++ipar)
   {
      const double cur_par = fCurPar[ipar];
      const double def_par = fDefPar[ipar];
      const double def_err = fDefErr[ipar];

      if(def_err > 0.0)
      {
         pterm += (cur_par - def_par)*(cur_par - def_par)/(def_err*def_err);
      }
   }

   return pterm;
}

//---------------------------------------------------------------------------------
const vector<double>& Anp::AlgAdapt::Get(const string &key) const
{
   if(key == "CurErr")
   {
      return fCurErr;
   }
   else if(key == "CurPar")
   {
      return fCurPar;
   }
   else if(key == "DefErr")
   {
      return fDefErr;
   }
   else if(key == "DefPar")
   {
      return fDefPar;
   }
   
   cerr << "AlgAdapt::Get - unknown key: " << key << endl;
   
   assert(false && "unknown key");

   return fCurPar;
}

//---------------------------------------------------------------------------------
void Anp::AlgAdapt::Set(const vector<double> &dvec, const string &key)
{
   if(dvec.size() != GetNPar())
   {
      cerr << "AlgAdapt::Set - number of parameters and size of input vector do not match: " 
           << " GetNPar() = " << GetNPar() << " and dvec.size() = " << dvec.size() << endl;
      return;
   }

   if(key == "CurErr")
   {
      fCurErr = dvec;
      return;
   }
   else if(key == "CurPar")
   {
      fCurPar = dvec;
      return;
   }
   else if(key == "DefErr")
   {
      fDefErr = dvec;
      return;
   }
   else if(key == "DefPar")
   {
      fDefPar = dvec;
      return;
   }
   
   cerr << "AlgAdapt::Set - unknown key: " << key << endl;

   assert(false && "unknown key");
}

//---------------------------------------------------------------------------------
void Anp::AlgAdapt::ReadVec(const string &prefix, const Registry &reg)
{
   //
   // Read in default and and current error and parameter vectors from
   // Registry, if these keys exist in Registry. Use virtual functions Set()
   // to set parameter values.
   //

   const vector<double> pvec = GetVec(prefix + "Par", reg);
   if(!pvec.empty())
   {
      Set(pvec, "CurPar");
   }

   const vector<double> evec = GetVec(prefix + "Err", reg);
   if(!evec.empty())
   {
      Set(evec, "CurErr");
   }

   const vector<double> evec_def = GetVec(prefix + "DefErr", reg);
   if(!evec_def.empty())
   {
      Set(evec_def, "DefErr");
   }

   const vector<double> pvec_def = GetVec(prefix + "DefPar", reg);
   if(!pvec_def.empty())
   {
      Set(pvec_def, "DefPar");
   }

   const vector<double> pvec_cur = GetVec(prefix + "CurPar", reg);
   if(!pvec_cur.empty())
   {
      Set(pvec_cur, "CurPar");
   }

   const vector<double> evec_cur = GetVec(prefix + "CurErr", reg);
   if(!evec_cur.empty())
   {
      Set(evec_cur, "CurErr");
   }
}

//---------------------------------------------------------------------------------
const vector<double> Anp::AlgAdapt::GetVec(const string &key, const Registry &reg) const
{
   //
   // Extract Registry value stored at key and convert string numbers separated by
   // " " or "," or both to vector of doubles.
   //
 
   vector<double> dvec;
   
   const char *value_char = 0;
   if(!reg.Get(key.c_str(), value_char) || !value_char)
   {
      return dvec;
   }

   vector<string> svec;
   UtilString::StringTok(svec, string(value_char), ", ");
   
   dvec = vector<double>(svec.size(), 0.0);

   for(unsigned int i = 0; i < svec.size(); ++i)
   {
      if(!svec[i].empty())
      {
         dvec[i] = std::atof(svec[i].c_str());
      }
   }

   return dvec;
}

//------------------------------------------------------------------------------------------
void Anp::AlgAdapt::PrintVec(std::ostream &o, const int w, const int p) const
{
   if(!fCurPar.empty())
   {
      o << "   CurPar = ";
      PrintVec(fCurPar, o, w, p);
   }
   
   if(!fCurErr.empty())
   {
      o << "   CurErr = ";
      PrintVec(fCurErr, o, w, p);
   }
   
   if(!fDefPar.empty())
   {
      o << "   DefPar = ";
      PrintVec(fDefPar, o, w, p);
   }
   
   if(!fDefErr.empty())
   {
      o << "   DefErr = ";
      PrintVec(fDefErr, o, w, p);
   }
}

//---------------------------------------------------------------------------------
void Anp::AlgAdapt::PrintVec(const vector<double> &v, ostream &o, const int w, const int p) const
{
   for(unsigned int i = 0; i < v.size(); ++i)
   {
      if(i == 0)
      {
         o << "(";
      }
      else
      {
         o << ", ";
      }
      
      if(p > 0)
      {
         o << setprecision(p);
      }
      if(w > 0)
      {
         o << setw(w);
      }

      o << v[i];
      
      if(i + 1 == v.size())
      {
         o << ")";
      }
   }

   o << endl;
}

//---------------------------------------------------------------------------------
void Anp::AlgAdapt::Print(std::ostream &os) const
{
   os << "AlgAdapt::Print - this function is not implemented for this class" << endl;
}

//---------------------------------------------------------------------------------
void Anp::AlgAdapt::Save(const string &fname)
{
   if(fCurErr.size() != fCurPar.size() || 
      fCurErr.size() != fDefErr.size() || 
      fCurErr.size() != fDefPar.size() || fCurErr.empty())
   {
      cerr << "AlgAdapt::Save - invalid internal state" << endl;
      return;
   }

   ofstream ofile(fname.c_str());
   if(!ofile.is_open())
   {
      cerr << "AlgAdapt::Save - failed to open output file:\n   " << fname << endl;
      return;
   }

   ofile << "CurErr = " << Anp::ConvertVector2String(fCurErr) << endl
         << "CurPar = " << Anp::ConvertVector2String(fCurPar) << endl
         << "DefErr = " << Anp::ConvertVector2String(fDefErr) << endl
         << "DefPar = " << Anp::ConvertVector2String(fDefPar) << endl;
}

//---------------------------------------------------------------------------------
void Anp::AlgAdapt::Read(const string &fname)
{
   std::ifstream infile(fname.c_str());
   if(!infile.is_open())
   {
      cerr << "AlgAdapt::Read - failed to open input file:\n   " << fname << endl;
      return;
   }

   vector<double> cur_err, cur_par, def_err, def_par;

   while(!infile.eof())
   {
      std::string line;
      std::getline(infile, line);
      
      if(line.empty())
      {
         continue;
      }

      const std::size_t ipos = line.find_last_of("=");
      if(ipos == string::npos)
      {
         cerr << "AlgAdapt::Read - failed to find equal sign:" << endl << line << endl;
         continue;
      }

      const string name = line.substr(0, ipos - 1);
      const string data = line.substr(ipos + 1);

      if     (name == "CurErr") cur_err = Anp::ConvertString2Vector(data);
      else if(name == "CurPar") cur_par = Anp::ConvertString2Vector(data);
      else if(name == "DefErr") def_err = Anp::ConvertString2Vector(data);
      else if(name == "DefPar") def_par = Anp::ConvertString2Vector(data);
   }
   
   if(cur_err.size() != cur_par.size() || 
      cur_err.size() != def_err.size() || 
      cur_err.size() != def_par.size() || cur_err.empty())
   {
      cerr << "AlgAdapt::Read - input vectors do not match" << endl;
      return;
   }

   Set(cur_err, "CurErr");
   Set(cur_par, "CurPar");
   Set(def_par, "DefErr");
   Set(def_par, "DefPar");
}

//---------------------------------------------------------------------------------
const string Anp::ConvertVector2String(const vector<double> &dvec)
{
   stringstream dstr;
   for(vector<double>::const_iterator dit = dvec.begin(); dit != dvec.end(); ++dit)
   {
      dstr << setprecision(9) << fixed << *dit;
      if(dit + 1 != dvec.end()) dstr << ", ";
   }
   
   return dstr.str();
}

//---------------------------------------------------------------------------------
const vector<double> Anp::ConvertString2Vector(const string &dstr)
{
   vector<double> dvec;

   //
   // Parse into individual strings
   //
   vector<string> svec;
   UtilString::StringTok(svec, dstr, ", ");

   if(svec.empty())
   {
      cerr << " Anp::ConvertString2Vector - input string has no data:\n   " << dstr << endl;
      return dvec;
   }

   for(unsigned int i = 0; i < svec.size(); ++i)
   {
      stringstream temp;
      temp << svec[i];

      double value = -1.0e9;
      temp >> value;

      if(!temp.fail())
      {
         dvec.push_back(value);
      }
      else
      {
         cerr << " Anp::ConvertString2Vector - stringstream >> failed for " << svec[i] << endl;
      }
   }

   return dvec;
}

---