GSL.cpp

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include "bias_config.h"
#include <Base/Common/W32Compat.hh>
#include <MathAlgo/GSL.hh>
#include <Base/Debug/Debug.hh>
#include <Base/Debug/Error.hh>
#include <Base/Debug/Exception.hh>
#include <Base/Common/CompareFloatingPoint.hh>

#include <cmath>

#ifdef BIAS_HAVE_GSL

#  include <gsl/gsl_cdf.h>
#  include <gsl/gsl_randist.h>
#  include <gsl/gsl_sf_gamma.h>

#endif

using namespace BIAS;
using namespace std;


double 
BIAS::ChiSquareProbDensFun(double x, int deg_freedom)
{
#ifdef BIAS_HAVE_GSL
  if (x<0.) 
    BEXCEPTION("ChiSquareProbDensFun: invalid argument, x must be > 0: "<<x);
  if (deg_freedom<1) 
    BEXCEPTION("ChiSquareProbDensFun: invalid argument, deg_fred must be > 0 "
               <<deg_freedom);
  double res = gsl_ran_chisq_pdf(x, deg_freedom);
  BIASASSERT(!isnan(res) && !isinf(res));
  return res;
#else
  BIASERR("recompile BIAS with gsl support enabled");
  BIASABORT;
  return 0.0;
#endif
}

double 
BIAS::ChiSquareCulmProbFun(double x, int deg_freedom)
{
#ifdef BIAS_HAVE_GSL
  double res = gsl_cdf_chisq_P(x, deg_freedom);
  BIASASSERT(!isnan(res) && !isinf(res));
  return res;
#else 
  BIASERR("recompile BIAS with gsl support enabled");
  BIASABORT;
  return 0.0;
#endif
}

double 
BIAS::InvChiSquareCulmProbFun(double x, int deg_freedom)
{
#ifdef BIAS_HAVE_GSL
  double res = gsl_cdf_chisq_Pinv(x, deg_freedom);
  BIASASSERT(!isnan(res) && !isinf(res));
  return res;
#else 
  BIASERR("recompile BIAS with gsl support enabled");
  BIASABORT;
  return 0.0;
#endif
}

double 
BIAS::Gamma(double x)
{
#ifdef BIAS_HAVE_GSL
  double res = gsl_sf_gamma (x);
  BIASASSERT(!isnan(res) && !isinf(res));
  return res;
#else 
  BIASERR("recompile BIAS with gsl support enabled");
  BIASABORT;
  return 0.0;
#endif
}

double 
BIAS::InvGamma(double x)
{
#ifdef BIAS_HAVE_GSL
  double res = gsl_sf_gammainv (x);
  BIASASSERT(!isnan(res) && !isinf(res));
  return res;
#else 
  BIASERR("recompile BIAS with gsl support enabled");
  BIASABORT;
  return 0.0;
#endif
}

double 
BIAS::BetaDistribution(double a, double b, double x)
{
#ifdef BIAS_HAVE_GSL
  if ((a+b)<=GSL_SF_GAMMA_XMAX){// check numerically stablility 
    double res = gsl_ran_beta_pdf(x, a, b);
    BIASASSERT(!isnan(res) && !isinf(res));
    return res;
  } else {
    BEXCEPTION("a or b out of range: a="<<a<<" b="<<b<<endl);
  }
#else 
  BIASERR("recompile BIAS with gsl support enabled");
  BIASABORT;
  return 0.0;
#endif
}


double 
BIAS::CulmBetaDistribution(double a, double b, double x)
{
#ifdef BIAS_HAVE_GSL
  double res = gsl_cdf_beta_P(x, a, b);
#ifdef BIAS_DEBUG
  if (isnan(res) || isinf(res)){
    cout << "a: "<<a<<"  b: "<<b<<"  x: "<<x<<endl;
  }
  cout << "a: "<<a<<"  b: "<<b<<"  x: "<<x<<endl;
#endif
  BIASASSERT(!isnan(res) && !isinf(res));
  return res;
#else 
  BIASERR("recompile BIAS with gsl support enabled");
  BIASABORT;
  return 0.0;
#endif
}


double 
BIAS::CulmulativeNormalDistribution(double x, double sigma)
{
#ifdef BIAS_HAVE_GSL
  double res = gsl_cdf_gaussian_P(x, sigma);
#ifdef BIAS_DEBUG
  if (isnan(res) || isinf(res)){
    cout << "sigma: "<<sigma<<"  x: "<<x<<"  res: "<<res<<endl;
  }
  cout << "sigma: "<<sigma<<"  x: "<<x<<"  res: "<<res<<endl;
#endif
  BIASASSERT(!isnan(res) && !isinf(res));
  return res;  
#else 
  BIASERR("recompile BIAS with gsl support enabled");
  BIASABORT;
  return 0.0;
#endif
}


double 
BIAS::BetaDistributionApprox(double a, double b, double x)
{
#ifdef BIAS_HAVE_GSL
  if ((a+b)<=GSL_SF_GAMMA_XMAX){// check numerically stablility 
    double res = gsl_ran_beta_pdf(x, a, b);
    BIASASSERT(!isnan(res) && !isinf(res));
    return res;
  } else if (Equal(a,b)){ 
    // the beta distribution is very steep and narrow, 
    // approximate beta distribution by gaussian
    double dist = x-0.5;
    // The variance of the beta distribution beta(a, b, x) is given by
    // sigma^2 = (a*b)/(a+b)^2/(a+b+1) 
    // (http://en.wikipedia.org/wiki/Beta_distribution)
    double sigma_sq = (a*b)/(a+b)/(a+b)/(a+b+1);
    double res = 1.0/sqrt(sigma_sq * 2.0 * M_PI) * 
      exp(-(dist*dist)/(2.0 * sigma_sq));
    bool debug = false;
    if (debug){
      cout <<"approximating beta distribution by gaussian: alpha = beta = "<<a
           <<"  x = "<<x<<"  sigma: "<<sqrt(sigma_sq)<<"  result: "<<res<<endl;
    }
    return res;
  } else {
    BIASERR("a or b to big, cannot compute beta distribution\n");
    BIASABORT;
    return -1.0;
  }
#else 
  BIASERR("recompile BIAS with gsl support enabled");
  BIASABORT;
  return 0.0;
#endif
}