$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/* Copyright 2008-2011 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express is 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.                                    */

/* QUB Express 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,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

#include <math.h>
#include <string.h>
#include <iostream>
#include "sampled_gaussian.h"

using namespace std;

#ifndef NULL
#define NULL 0
#endif

#ifdef _WIN32
#include <boost/math/special_functions/erf.hpp>
using boost::math::erf;
#endif

SampledGaussian::SampledGaussian()
  : A(NULL)
{
  setup(0.0, 1.0, 1e-3, 6.0, 256, false);
}

SampledGaussian::SampledGaussian(const SampledGaussian& other)
  : mean(other.mean), std(other.std), resolution(other.resolution), cutoff_std(other.cutoff_std),
    s2s2i(other.s2s2i), r2(other.r2), max_for_complete(other.max_for_complete), N(other.N),
    A(other.A ? new double[other.N] : NULL), A_out(other.A_out), y2i_m(other.y2i_m), y2i_b(other.y2i_b),
    use_log(other.use_log)
{
  if ( A )
    memcpy(A, other.A, N*sizeof(double));
}

SampledGaussian& SampledGaussian::operator=(const SampledGaussian& other)
{
  delete [] A;
  mean = other.mean;
  std = other.std;
  resolution = other.resolution;
  cutoff_std = other.cutoff_std;
  s2s2i = other.s2s2i;
  r2 = other.r2;
  max_for_complete = other.max_for_complete;
  N = other.N;
  A = other.A ? new double[other.N] : NULL;
  A_out = other.A_out;
  y2i_m = other.y2i_m;
  y2i_b = other.y2i_b;
  use_log = other.use_log;
  return *this;
}

void SampledGaussian::setup(double mean_, double std_, double resolution_,
			    double cutoff_std_, int max_for_complete_, bool sample, bool use_log_)
{
  delete [] A;
  A = NULL;
  
  mean = mean_;
  std = std_;
  resolution = resolution_;
  cutoff_std = cutoff_std_;
  max_for_complete = max_for_complete_;
  use_log = use_log_;
  
  s2s2i = 1.0 / (std * sqrt(2.0));
  r2 = resolution / 2.0;
  N = (int) (2 * std * cutoff_std / resolution + 0.5);
  if ( (! sample) || (N > max_for_complete) ) {
    //cerr << "non-sampled " << mean << " +/- " << std << " -> " << N << " at " << resolution << endl;
    N = 0;
  } else {
    make_A();
  }
}

void SampledGaussian::make_A()
{
  A = new double[N];
  double cdf_b = 0.0;
  double b = mean - std * cutoff_std;
  for (int i=0; i<N; ++i) {
    double a = b;
    double cdf_a = cdf_b ? cdf_b : erf((a - mean) * s2s2i);
    b += resolution;
    cdf_b = erf((b - mean) * s2s2i);
    A[i] = 0.5 * (cdf_b - cdf_a);
    if ( use_log )
      A[i] = log(A[i]);
  }
  A_out = G_area(mean+(cutoff_std+1)*std, mean+(cutoff_std+1)*std + resolution);
  y2i_m = N / (2 * std * cutoff_std);
  y2i_b = N/2.0 - mean * y2i_m;
}

SampledGaussian::~SampledGaussian()
{
  delete [] A;
}

double SampledGaussian::operator()(double x)
{
  return N ? G_sampled(x) : G_centered(x);
}

double SampledGaussian::G_area(double a, double b)
{
  if ( use_log )
    return log(0.5 * (erf((b - mean) * s2s2i) - erf((a - mean) * s2s2i)));
  else
    return 0.5 * (erf((b - mean) * s2s2i) - erf((a - mean) * s2s2i));
}

double SampledGaussian::G_centered(double x)
{
  return G_area(x - r2, x + r2);
}

double SampledGaussian::G_sampled(double x)
{
  int i = (int) (y2i_m * x + y2i_b);
  if ( (0 <= i) && (i < N) )
    return A[i];
  return A_out;
}