TransformedParameter.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1

#ifndef BPP_NUMERIC_TRANSFORMEDPARAMETER_H
#define BPP_NUMERIC_TRANSFORMEDPARAMETER_H

#include <cmath>

#include "NumConstants.h"
#include "Parameter.h"

namespace bpp
{
class TransformedParameter :
  public Parameter
{
public:
  TransformedParameter(const std::string& name, double value) :
    Parameter(name, value) {}

  TransformedParameter* clone() const = 0;

public:
  virtual void setOriginalValue(double value) = 0;

  virtual double getOriginalValue() const = 0;

  virtual double getFirstOrderDerivative() const = 0;

  virtual double getSecondOrderDerivative() const = 0;
};

class RTransformedParameter :
  public TransformedParameter
{
private:
  double scale_;
  double bound_;
  bool positive_;

public:
  RTransformedParameter(const std::string& name, double value, double bound = 0, bool positive = true, double scale = 1) :
    TransformedParameter(name, 1.),
    scale_(scale),
    bound_(bound),
    positive_(positive)
  {
    setOriginalValue(value);
  }

  RTransformedParameter* clone() const { return new RTransformedParameter(*this); }

public:
  void setOriginalValue(double value)
  {
    if (positive_ ? value <= bound_ : value >= bound_) throw ConstraintException("RTransformedParameter::setValue", this, value);
    if (positive_  & (value < 1 + bound_)) setValue(log(scale_ * (value - bound_)));
    if (positive_  & (value >= 1 + bound_)) setValue(scale_ * (value - 1. - bound_));
    if (!positive_ & (value > -1 + bound_)) setValue(log(-scale_ * (value - bound_)));
    if (!positive_ & (value <= -1 + bound_)) setValue(-scale_ * (value - 1. - bound_));
  }

  double getOriginalValue() const
  {
    double x = getValue();
    if (positive_)
      if (x < 0) return exp(x) / scale_ + bound_;
      else return x / scale_ + 1. + bound_;
    else if (x < 0) return -exp(-x) / scale_ + bound_;
    else return -x / scale_ - 1. + bound_;
  }

  double getFirstOrderDerivative() const
  {
    double x = getValue();
    if (positive_)
      if (x < 0) return exp(x) / scale_;
      else return 1. / scale_;
    else if (x < 0) return exp(-x) / scale_;
    else return -1. / scale_;
  }

  double getSecondOrderDerivative() const
  {
    double x = getValue();
    if (positive_)
      if (x < 0) return exp(x) / scale_;
      else return 0;
    else if (x < 0) return -exp(-x) / scale_;
    else return 0;
  }
};

class IntervalTransformedParameter :
  public TransformedParameter
{
private:
  double scale_;
  double lowerBound_;
  double upperBound_;
  bool hyper_;
  double tiny_;

public:
  IntervalTransformedParameter(const std::string& name, double value, double lowerBound = 0, double upperBound = 1, double scale = 1, bool hyper = true) :
    TransformedParameter(name, hyper ?
        scale * atanh(2. * (value - lowerBound) / (upperBound - lowerBound) - 1.) :
        scale * tan(NumConstants::PI() * (value - lowerBound) / (upperBound - lowerBound) - NumConstants::PI() / 2.)),
    scale_(scale),
    lowerBound_(lowerBound),
    upperBound_(upperBound),
    hyper_(hyper),
    tiny_(NumConstants::TINY())
  {}

  IntervalTransformedParameter* clone() const { return new IntervalTransformedParameter(*this); }

public:
  void setOriginalValue(double value)
  {
    if (value <= lowerBound_ || value >= upperBound_) throw ConstraintException("IntervalTransformedParameter::setValue", this, value);
    setValue(hyper_ ?
        scale_ * atanh(2. * (value - lowerBound_) / (upperBound_ - lowerBound_) - 1.) :
        scale_ * std::tan(NumConstants::PI() * (value - lowerBound_) / (upperBound_ - lowerBound_) - NumConstants::PI() / 2.));
  }

  double getOriginalValue() const
  {
    double x = getValue();
    double x2 = hyper_ ?
        (tanh(x / scale_) + 1.) * (upperBound_ - lowerBound_) / 2. + lowerBound_ :
        (atan(x / scale_) + NumConstants::PI() / 2.) * (upperBound_ - lowerBound_) / NumConstants::PI() + lowerBound_;
    return x2;
  }


  double getFirstOrderDerivative() const
  {
    double x = getValue();
    double x2 = hyper_ ?
        1. / (std::pow(cosh(x / scale_), 2)) * (upperBound_ - lowerBound_) / (2. * scale_) :
        (upperBound_ - lowerBound_) / (NumConstants::PI() * scale_ * (std::pow(x / scale_, 2) + 1.));
    return x2;
  }
  double getSecondOrderDerivative() const
  {
    double x = getValue();
    double x2 = hyper_ ?
        -1. / (std::pow(cosh(x / scale_), 2)) * tanh(x / scale_) * (upperBound_ - lowerBound_) / (scale_ * scale_) :
        -2. * x * (upperBound_ - lowerBound_) / (NumConstants::PI() * std::pow(scale_, 3) * std::pow((std::pow(x / scale_, 2) + 1.), 2));
    return x2;
  }
};

class PlaceboTransformedParameter :
  public TransformedParameter
{
public:
  PlaceboTransformedParameter(const std::string& name, double value) :
    TransformedParameter(name, value)
  {}

  PlaceboTransformedParameter* clone() const { return new PlaceboTransformedParameter(*this); }

public:
  void setOriginalValue(double value)
  {
    setValue(value);
  }

  double getOriginalValue() const
  {
    return getValue();
  }

  double getFirstOrderDerivative() const { return 1.; }

  double getSecondOrderDerivative() const { return 0.; }
};
} // end of namespace bpp.
#endif // BPP_NUMERIC_TRANSFORMEDPARAMETER_H