bcedistr.hpp

// This file is part of the BCESolve library for games of incomplete
// information
// Copyright (C) 2022 Benjamin A. Brooks
// 
// BCESolve 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.
// 
// BCESolve 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
// <http://www.gnu.org/licenses/>.
// 
// Benjamin A. Brooks
// ben@benjaminbrooks.net
// Chicago, IL
 
#ifndef BCEDISTR_HPP
#define BCEDISTR_HPP
 
#include <cmath>
#include <vector>
 
using namespace std;
 
 
class BCEDistr
{
public:
  virtual double CDF(double v0, double v1) const = 0;
 
 
  double PDF(double v0, double incr0,
         double v1, double incr1) const
  {
    return (CDF(v0,v1) + CDF(v0-incr0,v1-incr1)
        - CDF(v0-incr0,v1) - CDF(v0,v1-incr1));
  }
 
 
  double PDF(double v0, double v1, double incr) const
  {
    return PDF(v0,incr,v1,incr);
  }
 
  virtual ~BCEDistr() {}
};
 
 
class BCEDistrArray : public BCEDistr
{
public:
  vector<BCEDistr*> distributions;
  vector<double> weights;
 
 
  BCEDistrArray() {}
  
 
  ~BCEDistrArray()
  {
    vector<BCEDistr*>::iterator distr;
    for (distr = distributions.begin();
     distr != distributions.end();
     distr++)
      delete *distr;
  }
 
 
  bool push_back(BCEDistr* newDistr, double newWeight)
  {
    // Avoid recursive catastrophes!
    if (!contains(newDistr))
      {
    distributions.push_back(newDistr);
    weights.push_back(newWeight);
    return false;
      }
    return true;
  }
 
 
  bool contains(BCEDistr * distr) const
  {
    // First check if this object is itself distr.
    if (this==distr)
      return true;
    
    // Now check if any of the elements of distributions are distr. If
    // any of the elements of distributions are themselves arrays,
    // initiate a recursive check.
    for (vector<BCEDistr*>::const_iterator it = distributions.begin();
     it != distributions.end();
     ++it)
      {
    BCEDistrArray * distrArray = dynamic_cast<BCEDistrArray *>(distr);
    
    if (distr == *it 
        || (distrArray != NULL
        && distrArray->contains(distr) ) )
      return true;
      }
 
    return false;
  }
 
  void clear()
  {
    vector<BCEDistr*>::iterator distr;
    for (distr = distributions.begin();
     distr != distributions.end();
     distr++)
      delete *distr;
    distributions.clear();
    weights.clear();
  }
 
  double CDF(double v0, double v1) const
  {
    vector<BCEDistr*>::const_iterator distr;
    vector<double>::const_iterator weight;
 
    double p = 0;
    for (distr = distributions.begin(),
        weight = weights.begin();
     distr != distributions.end();
     distr++, weight++)
      p += *weight * (*distr)->CDF(v0,v1);
    return p;
  }
};
 
 
class unimodal: public BCEDistr
{
  double m0;
  double c0;
  double d0;
  double e0;
  int k0;
 
  double m1;
  double c1;
  double d1;
  double e1;
  int k1;
 
public:
  unimodal(double _c0, double _m0, int _k0,
       double _c1, double _m1, int _k1):
    c0(_c0), k0(_k0), m0(_m0),
    c1(_c1), k1(_k1), m1(_m1)
  {
    e0 = pow(-m0,k0);
    d0 = 1.0/c0+pow(1-m0,k0)-e0;
    e1 = pow(-m1,k1);
    d1 = 1.0/c1+pow(1-m1,k1)-e1;
  }
 
  unimodal(double _c, double _m, int _k):
    c0(_c), k0(_k), m0(_m),
    c1(_c), k1(_k), m1(_m)
  {
    e0 = pow(-m0,k0);
    d0 = 1.0/c0+pow(1-m0,k0)-e0;
    e1 = e0; d1 = d0;
  }
 
  double CDF (double v0, double v1) const
  {
    double p = 1.0;
    if (v0<0)
      return 0.0;
    else if (v0<1.0)
      p *= c0*( d0*v0 + pow(v0-m0,k0) + e0 );
 
    if (v1<0)
      return 0.0;
    else if (v1 < 1)
      p *= c1*( d1*v1 + pow(v1-m1,k1) + e1 );
    
    return p;
  }
 
}; // unimodal
 
class additiveCorrelated: public BCEDistr
{
public:
  additiveCorrelated() {}
 
  double CDF(double v0, double v1) const
  {
    v0 = min(max(v0,0.0),1.0);
    v1 = min(max(v1,0.0),1.0);
    return (v0+v1)*(v0+v1)*v0*v1/4;
  }
}; // additiveCorrelated
 
class uniform: public BCEDistr
{
public:
  uniform() {};
 
  double CDF(double v0, double v1) const
  {
    double p = 1.0;
    if (v0<0)
      return 0.0;
    else if (v0<1)
      p *= v0;
 
    if (v1<0)
      return 0.0;
    else if (v1 < 1)
      p *= v1;
    
    return p;
  }
 
}; // uniform
 
class bernoulli: public BCEDistr
{
private:
  double lowProb;
public:
  bernoulli(double _lowProb): lowProb(_lowProb) {};
 
  double CDF(double v0, double v1) const
  {
    double p = 1.0;
    if (v0<0)
      return 0.0;
    else if (v0<1)
      p *= lowProb;
 
    if (v1<0)
      return 0.0;
    else if (v1 < 1)
      p *= lowProb;
    
    return p;
  }
 
}; // bernoulli
 
class uniformWithMassPnt: public BCEDistr
{
  double alpha;
 
public:
  uniformWithMassPnt() {}
  uniformWithMassPnt(double _alpha): alpha(_alpha) {}
 
  double CDF(double v0, double v1) const
  {
    double p = 1.0;
    if (v0<-1e-10)
      return 0.0;
    else if (v0<1e-10)
      p *= alpha;
    else if (v0<1)
      p *= (1-alpha)*v0;
 
    if (v1<-1e-10)
      return 0.0;
    else if (v1<1e-10)
      p *= alpha;
    else if (v1 < 1)
      p *= (1-alpha)*v1;
    
    return p;
  }
 
}; // uniform
 
class uniformNoTie: public BCEDistr
{
  int k;
  
public:
  uniformNoTie() {}
  uniformNoTie(int _k): k(_k) {}
 
  double CDF(double v0, double v1) const
  {
    double p = 1.0;
    if (v0<0)
      return 0.0;
    else if (v0<1)
      p *= v0;
 
    if (v1<0)
      return 0.0;
    else if (v1 < 1)
      p *= v1;
    
    if (abs(v0-v1)<(0.5/k))
      p -= 1.0/(k*k);
 
    return p;
  }
 
}; // uniformNoTie
 
class vToTheAlpha: public BCEDistr
{
  double alpha0;
  double alpha1;
 
public:
  vToTheAlpha() {}
 
  vToTheAlpha(double a0, double a1): 
    alpha0(a0), 
    alpha1(a1)
  {}
 
  vToTheAlpha(double a): 
    alpha0(a), 
    alpha1(a)
  {}
 
  double CDF(double v0, double v1) const
  {
    double p = 1.0;
    if (v0<0)
      return 0.0;
    else if (v0<1)
      p *= pow(v0,alpha0);
 
    if (v1<0)
      return 0.0;
    else if (v1 < 1)
      p *= pow(v1,alpha1);
    
    return p;
  } // CDF
 
}; // vToTheAlpha
 
class truncated: public BCEDistr
{
public:
  BCEDistr * dist;
  double min0;
  double min1;
  double max0;
  double max1;
 
 
  truncated(BCEDistr * _dist):
    dist(_dist),
    min0(0),
    min1(0),
    max0(1),
    max1(1)
  {}
 
 
  truncated(BCEDistr * _dist, 
        double _min, double _max): 
    dist(_dist),
    min0(_min),
    min1(_min),
    max0(_max),
    max1(_max)
  {}
 
 
  truncated(BCEDistr * _dist, 
        double _min0, double _min1,
        double _max0, double _max1): 
    dist(_dist),
    min0(_min0),
    min1(_min1),
    max0(_max0),
    max1(_max1)
  {
    assert(max0 >= min0);
    assert(max1 >= min1);
  }
 
  double CDF(double v0, double v1) const
  {
    // if (v0 < min0)
    //   v0 = min0;
    // else if (v0 > max0)
    //   v0 = max0;
 
    // if (v1 < min1)
    //   v1 = min1;
    // else if (v1 > max1)
    //   v1 = max1;
 
    if (v0 < min0-1e-13 || v1 < min1-1e-13)
      return 0;
 
    return dist->CDF(v0,v1)-dist->CDF(min0,min1);
  } // CDF
};
 
#endif