DataPointSet.h

// -*- C++ -*-
//
// DataPointSet.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef LWH_DataPointSet_H
#define LWH_DataPointSet_H
//
// This is the declaration of the DataPointSet class representing
//


#include <vector>
#include <limits>
#include <cmath>
#include <algorithm>
#include "AIDataPointSet.h"
#include "ManagedObject.h"
#include "DataPoint.h"

namespace LWH {

using namespace AIDA;

class DataPointSet: public IDataPointSet, public ManagedObject {

public:

  DataPointSet(int D): dim(D) {}

  virtual ~DataPointSet() {}

  IAnnotation & annotation() {
    throw std::runtime_error("LWH cannot handle annotations");
  }

  const IAnnotation & annotation() const {
    throw std::runtime_error("LWH cannot handle annotations");
  }

  std::string title() const {
    return theTitle;
  }

  std::string name() const {
    return theTitle;
  }

  bool setTitle(const std::string & title) {
    theTitle = title;
    return true;
  }

  int dimension() const {
    return dim;
  }

  void clear() {
    dset.clear();
  }

  int size() const {
    return dset.size();
  }

  IDataPoint * point(int index) {
    return &(dset[index]);
  }

  bool setCoordinate(int coord,
                     const std::vector<double>  & val,
                     const std::vector<double>  & err) {
    return setCoordinate(coord, val, err, err);
  }

  bool setCoordinate(int coord,
                     const std::vector<double>  & val,
                     const std::vector<double>  & errp,
                     const std::vector<double>  & errm) {
    if ( coord < 0 || coord >= dimension() ) return false;
    if ( val.size() != dset.size() || errp.size() != dset.size() ||
         errm.size() != dset.size() ) return false;
    for ( int i = 0, N = val.size(); i < N; ++i ) {
      dset[i].coordinate(coord)->setValue(val[i]);
      dset[i].coordinate(coord)->setErrorPlus(errp[i]);
      dset[i].coordinate(coord)->setErrorMinus(errm[i]);
    }
    return true;
  }

  const IDataPoint * point(int index) const {
    if ( index < 0 || unsigned(index) >= dset.size() ) return 0;
    return &(dset[index]);
  }

  IDataPoint * addPoint() {
    dset.push_back(DataPoint(dimension()));
    return &(dset.back());
  }

  bool addPoint(const IDataPoint & point) {
    if ( dimension() && dimension() != point.dimension() ) return false;
    dset.push_back(DataPoint(point));
    return true;
  }

  bool removePoint(int index) {
    if ( index < 0 || unsigned(index) >= dset.size() ) return false;
    dset.erase(dset.begin() + index);
    return true;
  }

  double lowerExtent(int coord) const {
    if ( dset.empty() ) return std::numeric_limits<double>::quiet_NaN();
    if ( coord < 0 || coord >= dimension() )
      return std::numeric_limits<double>::quiet_NaN();
    double low = dset[0].coordinate(coord)->value();
    for ( int i = 1, N = dset.size(); i < N; ++i )
      low = std::min(low, dset[i].coordinate(coord)->value());
    return low;
  }

  double upperExtent(int coord) const {
    if ( dset.empty() ) return std::numeric_limits<double>::quiet_NaN();
    if ( coord < 0 || coord >= dimension() )
      return std::numeric_limits<double>::quiet_NaN();
    double upp = dset[0].coordinate(coord)->value();
    for ( int i = 1, N = dset.size(); i < N; ++i )
      upp = std::max(upp, dset[i].coordinate(coord)->value());
    return upp;
  }

  bool scale(double scale) {
    for ( int i = 0, N = dset.size(); i < N; ++i )
      for ( int j = 0, M = dset[i].dimension(); j < M; ++j ) {
        IMeasurement * m = dset[i].coordinate(j);
        m->setValue(m->value()*scale);
        m->setErrorPlus(m->errorPlus()*scale);
        m->setErrorMinus(m->errorPlus()*scale);
      }
    return true;
  }
        
  bool scaleValues(double scale) {
    for ( int i = 0, N = dset.size(); i < N; ++i )
      for ( int j = 0, M = dset[i].dimension(); j < M; ++j ) {
        IMeasurement * m = dset[i].coordinate(j);
        m->setValue(m->value()*scale);
      }
    return true;
  }

  bool scaleErrors(double scale) {
    for ( int i = 0, N = dset.size(); i < N; ++i )
      for ( int j = 0, M = dset[i].dimension(); j < M; ++j ) {
        IMeasurement * m = dset[i].coordinate(j);
        m->setErrorPlus(m->errorPlus()*scale);
        m->setErrorMinus(m->errorPlus()*scale);
      }
    return true;
  }

  void * cast(const std::string &) const {
    return 0;
  }

  bool writeXML(std::ostream & os, std::string path, std::string name) {
    os << "  <dataPointSet name=\"" << name
       << "\"\n    title=\"" << title()
       << "\" path=\"" << path
       << "\" dimension=\"" << dimension() << "\">\n";
    for ( int d = 0; d < dimension(); ++d )
      os << "    <dimension dim=\"" << d << "\" title=\"unknown\" />\n";
    for ( int i = 0, N = size(); i < N; ++i ) {
      os << "    <dataPoint>\n";
      for ( int j = 0, M = dimension(); j < M; ++j )
        os << "      <measurement value=\""
           << point(i)->coordinate(j)->value()
           << "\" errorPlus=\""
           << point(i)->coordinate(j)->errorPlus()
           << "\" errorMinus=\""
           << point(i)->coordinate(j)->errorMinus()
           << "\"/>\n";
      os << "    </dataPoint>\n";
    }
    os << "  </dataPointSet>" << std::endl;
    return true;
  }

  bool writeFLAT(std::ostream & os, std::string path, std::string name) {
    os << "# " << path << "/" << name << " " << size()
       << " \"" << title() << " \" dimension " << dimension() << std::endl;
    for ( int i = 0, N = size(); i < N; ++i ) {
      for ( int j = 0, M = dimension(); j < M; ++j )
        os << point(i)->coordinate(j)->value() << " ";
      for ( int j = 0, M = dimension(); j < M; ++j )
        os << point(i)->coordinate(j)->errorPlus() << " ";
      for ( int j = 0, M = dimension(); j < M; ++j )
        os << point(i)->coordinate(j)->errorMinus() << " ";
      os << std::endl;
    }
    os << std::endl;
    return true;
  }

private:

  std::string theTitle;

  std::vector<DataPoint> dset;

  unsigned int dim;


};

}

#endif /* LWH_DataPointSet_H */