yap-6.3/packages/ProbLog/simplecudd/problogmath.c

/******************************************************************************\
*                                                                              *
*    SimpleCUDD library (www.cs.kuleuven.be/~theo/tools/simplecudd.html)       *
*  SimpleCUDD was developed at Katholieke Universiteit Leuven(www.kuleuven.be) *
*                                                                              *
*  Copyright Katholieke Universiteit Leuven 2008, 2009, 2010                   *
*                                                                              *
*  Author: Bernd Gutmann                                                       *
*  File: problogmath.c                                                         *
*  $Date:: 2010-10-06 13:20:59 +0200 (Wed, 06 Oct 2010)                      $ *
*  $Revision:: 4880                                                          $ *
*                                                                              *
********************************************************************************
*                                                                              *
* Artistic License 2.0                                                         *
*                                                                              *
* Copyright (c) 2000-2006, The Perl Foundation.                                *
*                                                                              *
* Everyone is permitted to copy and distribute verbatim copies of this license *
* document, but changing it is not allowed.                                    *
*                                                                              *
* Preamble                                                                     *
*                                                                              *
* This license establishes the terms under which a given free software Package *
* may be copied, modified, distributed, and/or redistributed. The intent is    *
* that the Copyright Holder maintains some artistic control over the           *
* development of that Package while still keeping the Package available as     *
* open source and free software.                                               *
*                                                                              *
* You are always permitted to make arrangements wholly outside of this license *
* directly with the Copyright Holder of a given Package. If the terms of this  *
* license do not permit the full use that you propose to make of the Package,  *
* you should contact the Copyright Holder and seek a different licensing       *
* arrangement.                                                                 *
* Definitions                                                                  *
*                                                                              *
* "Copyright Holder" means the individual(s) or organization(s) named in the   *
* copyright notice for the entire Package.                                     *
*                                                                              *
* "Contributor" means any party that has contributed code or other material to *
* the Package, in accordance with the Copyright Holder's procedures.           *
*                                                                              *
* "You" and "your" means any person who would like to copy, distribute, or     *
* modify the Package.                                                          *
*                                                                              *
* "Package" means the collection of files distributed by the Copyright Holder, *
* and derivatives of that collection and/or of those files. A given Package    *
* may consist of either the Standard Version, or a Modified Version.           *
*                                                                              *
* "Distribute" means providing a copy of the Package or making it accessible   *
* to anyone else, or in the case of a company or organization, to others       *
* outside of your company or organization.                                     *
*                                                                              *
* "Distributor Fee" means any fee that you charge for Distributing this        *
* Package or providing support for this Package to another party. It does not  *
* mean licensing fees.                                                         *
*                                                                              *
* "Standard Version" refers to the Package if it has not been modified, or has *
* been modified only in ways explicitly requested by the Copyright Holder.     *
*                                                                              *
* "Modified Version" means the Package, if it has been changed, and such       *
* changes were not explicitly requested by the Copyright Holder.               *
*                                                                              *
* "Original License" means this Artistic License as Distributed with the       *
* Standard Version of the Package, in its current version or as it may be      *
* modified by The Perl Foundation in the future.                               *
*                                                                              *
* "Source" form means the source code, documentation source, and configuration *
* files for the Package.                                                       *
*                                                                              *
* "Compiled" form means the compiled bytecode, object code, binary, or any     *
* other form resulting from mechanical transformation or translation of the    *
* Source form.                                                                 *
* Permission for Use and Modification Without Distribution                     *
*                                                                              *
* (1) You are permitted to use the Standard Version and create and use         *
* Modified Versions for any purpose without restriction, provided that you do  *
* not Distribute the Modified Version.                                         *
* Permissions for Redistribution of the Standard Version                       *
*                                                                              *
* (2) You may Distribute verbatim copies of the Source form of the Standard    *
* Version of this Package in any medium without restriction, either gratis or  *
* for a Distributor Fee, provided that you duplicate all of the original       *
* copyright notices and associated disclaimers. At your discretion, such       *
* verbatim copies may or may not include a Compiled form of the Package.       *
*                                                                              *
* (3) You may apply any bug fixes, portability changes, and other              *
* modifications made available from the Copyright Holder. The resulting        *
* Package will still be considered the Standard Version, and as such will be   *
* subject to the Original License.                                             *
* Distribution of Modified Versions of the Package as Source                   *
*                                                                              *
* (4) You may Distribute your Modified Version as Source (either gratis or for *
* a Distributor Fee, and with or without a Compiled form of the Modified       *
* Version) provided that you clearly document how it differs from the Standard *
* Version, including, but not limited to, documenting any non-standard         *
* features, executables, or modules, and provided that you do at least ONE of  *
* the following:                                                               *
*                                                                              *
* (a) make the Modified Version available to the Copyright Holder of the       *
* Standard Version, under the Original License, so that the Copyright Holder   *
* may include your modifications in the Standard Version.                      *
* (b) ensure that installation of your Modified Version does not prevent the   *
* user installing or running the Standard Version. In addition, the Modified   *
* Version must bear a name that is different from the name of the Standard     *
* Version.                                                                     *
* (c) allow anyone who receives a copy of the Modified Version to make the     *
* Source form of the Modified Version available to others under                *
* (i) the Original License or                                                  *
* (ii) a license that permits the licensee to freely copy, modify and          *
* redistribute the Modified Version using the same licensing terms that apply  *
* to the copy that the licensee received, and requires that the Source form of *
* the Modified Version, and of any works derived from it, be made freely       *
* available in that license fees are prohibited but Distributor Fees are       *
* allowed.                                                                     *
* Distribution of Compiled Forms of the Standard Version or Modified Versions  *
* without the Source                                                           *
*                                                                              *
* (5) You may Distribute Compiled forms of the Standard Version without the    *
* Source, provided that you include complete instructions on how to get the    *
* Source of the Standard Version. Such instructions must be valid at the time  *
* of your distribution. If these instructions, at any time while you are       *
* carrying out such distribution, become invalid, you must provide new         *
* instructions on demand or cease further distribution. If you provide valid   *
* instructions or cease distribution within thirty days after you become aware *
* that the instructions are invalid, then you do not forfeit any of your       *
* rights under this license.                                                   *
*                                                                              *
* (6) You may Distribute a Modified Version in Compiled form without the       *
* Source, provided that you comply with Section 4 with respect to the Source   *
* of the Modified Version.                                                     *
* Aggregating or Linking the Package                                           *
*                                                                              *
* (7) You may aggregate the Package (either the Standard Version or Modified   *
* Version) with other packages and Distribute the resulting aggregation        *
* provided that you do not charge a licensing fee for the Package. Distributor *
* Fees are permitted, and licensing fees for other components in the           *
* aggregation are permitted. The terms of this license apply to the use and    *
* Distribution of the Standard or Modified Versions as included in the         *
* aggregation.                                                                 *
*                                                                              *
* (8) You are permitted to link Modified and Standard Versions with other      *
* works, to embed the Package in a larger work of your own, or to build        *
* stand-alone binary or bytecode versions of applications that include the     *
* Package, and Distribute the result without restriction, provided the result  *
* does not expose a direct interface to the Package.                           *
* Items That are Not Considered Part of a Modified Version                     *
*                                                                              *
* (9) Works (including, but not limited to, modules and scripts) that merely   *
* extend or make use of the Package, do not, by themselves, cause the Package  *
* to be a Modified Version. In addition, such works are not considered parts   *
* of the Package itself, and are not subject to the terms of this license.     *
* General Provisions                                                           *
*                                                                              *
* (10) Any use, modification, and distribution of the Standard or Modified     *
* Versions is governed by this Artistic License. By using, modifying or        *
* distributing the Package, you accept this license. Do not use, modify, or    *
* distribute the Package, if you do not accept this license.                   *
*                                                                              *
* (11) If your Modified Version has been derived from a Modified Version made  *
* by someone other than you, you are nevertheless required to ensure that your *
* Modified Version complies with the requirements of this license.             *
*                                                                              *
* (12) This license does not grant you the right to use any trademark, service *
* mark, tradename, or logo of the Copyright Holder.                            *
*                                                                              *
* (13) This license includes the non-exclusive, worldwide, free-of-charge      *
* patent license to make, have made, use, offer to sell, sell, import and      *
* otherwise transfer the Package with respect to any patent claims licensable  *
* by the Copyright Holder that are necessarily infringed by the Package. If    *
* you institute patent litigation (including a cross-claim or counterclaim)    *
* against any party alleging that the Package constitutes direct or            *
* contributory patent infringement, then this Artistic License to you shall    *
* terminate on the date that such litigation is filed.                         *
*                                                                              *
* (14) Disclaimer of Warranty: THE PACKAGE IS PROVIDED BY THE COPYRIGHT HOLDER *
* AND CONTRIBUTORS "AS IS' AND WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES. THE  *
* IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR  *
* NON-INFRINGEMENT ARE DISCLAIMED TO THE EXTENT PERMITTED BY YOUR LOCAL LAW.   *
* UNLESS REQUIRED BY LAW, NO COPYRIGHT HOLDER OR CONTRIBUTOR WILL BE LIABLE    *
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING IN    *
* ANY WAY OUT OF THE USE OF THE PACKAGE, EVEN IF ADVISED OF THE POSSIBILITY OF *
* SUCH DAMAGE.                                                                 *
*                                                                              *
*         The End                                                              *
*                                                                              *
\******************************************************************************/

#include "problogmath.h"


double sigmoid(double x, double slope) {
  return 1.0 / (1.0 + exp(-x * slope));
}

// This function calculates the accumulated density of the normal distribution
// For details see G. Marsaglia, Evaluating the Normal Distribution, Journal of Statistical Software, 2004:11(4).
double Phi(double x)  {
  double s=x;
  double t=0.0;
  double b=x;
  double q=x*x;
  double i=1;

  // if the value is too small or too big, return
  // 0/1 to avoid long computations
  if (x < -10.0) {
    return 0.0;
  }

  if (x > 10.0) {
    return 1.0;
  }

  // t is the value from last iteration
  // s is the value from the current iteration
  // iterate until they are equal
  while(fabs(s-t) >= DBL_MIN) {
    t=s;
    i+=2;
    b*=q/i;
    s+=b;
  }

  return 0.5+s*exp(-0.5*q-0.91893853320467274178);
}

// integrates the normal distribution over [low,high]
double cumulative_normal(double low, double high, double mu, double sigma) {
  return Phi((high-mu)/sigma) - Phi((low-mu)/sigma);
}

// integrates the normal distribution over [-oo,high]
double cumulative_normal_upper(double high, double mu, double sigma) {
  return Phi((high-mu)/sigma);
}


// evaluates the density of the normal distribution
double normal(double x, double mu,double sigma) {
  double inner=(x-mu)/sigma;
  double denom=sigma*sqrt(2*3.14159265358979323846);
  return exp(-inner*inner/2)/denom;
}

double cumulative_normal_dmu(double low, double high,double mu,double sigma) {
  return normal(low,mu,sigma) - normal(high,mu,sigma);
}

double cumulative_normal_upper_dmu(double high,double mu,double sigma) {
  return  - normal(high,mu,sigma);
}


double cumulative_normal_dsigma(double low, double high,double mu,double sigma) {
  return (((mu-high)*normal(high,mu,sigma) - (mu-low)*normal(low,mu,sigma))/sigma);
}

double cumulative_normal_upper_dsigma(double high,double mu,double sigma) {
  return (mu-high)*normal(high,mu,sigma);
}


// this function parses two strings "$a;$b" and "???_???l$ch$d" where $a-$d are (real) numbers
// it is used to parse in the parameters of continues variables from the input file
density_integral parse_density_integral_string(char *input, char *variablename) {
  density_integral result;
  double sigma;
  int i;
  char garbage[64], s1[64],s2[64],s3[64],s4[64];

  if(sscanf(input, "%64[^;];%64[^;]", s1,s2) != 2) {
    fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string\n",input);
    fprintf(stderr, "The string should contain 2 fields seperated by ; characters.\n");
    exit(EXIT_FAILURE);
  }

  if (!getRealNumber(s1, &result.mu)) {
    fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string\n",input);
    fprintf(stderr, "%s is not a number\n",s1);
    exit(EXIT_FAILURE);
  }

  if (!getRealNumber(s2, &sigma) || sigma<=0.0) {
    fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string\n",input);
    fprintf(stderr, "%s is not a number\n",s2);
    exit(EXIT_FAILURE);
  }
  result.log_sigma=log(sigma);

/*  if (result.sigma<=0) { */
/*     fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string",input); */
/*     fprintf(stderr, "The value for sigma has to be larger than 0.\n"); */

/*     exit(EXIT_FAILURE); */
/*   } */

  if (sscanf(variablename,"%64[^lh]l%64[^lh]h%64[^lh]",garbage,s3,s4) != 3) {
    fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string\n",variablename);
    fprintf(stderr, "The string should contain 2 fields seperated by ; characters.\n");
    exit(EXIT_FAILURE);
  }

  //  replace the d by . in s1 and s2
  for(i=0; s3[i]!='\0' ; i++) {
    if (s3[i]=='d') {
      s3[i]='.';
    }
    if (s3[i]=='m') {
      s3[i]='-';
    }
  }
  for(i=0; s4[i]!='\0' ; i++) {
    if (s4[i]=='d') {
      s4[i]='.';
    }
    if (s4[i]=='m') {
      s4[i]='-';
    }
  }

  if (!getRealNumber(s3, &result.low)) {
    fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string\n",input);
    fprintf(stderr, "%s is not a number\n",s1);
    exit(EXIT_FAILURE);
  }

 if (!getRealNumber(s4, &result.high)) {
    fprintf(stderr, "Error ar parsing the string %s in the function parse_density_integral_string\n",input);
    fprintf(stderr, "%s is not a number\n",s1);
    exit(EXIT_FAILURE);
  }

  
  if (result.low>result.high) {
    fprintf(stderr, "Error ar parsing the string %s in the function parse_density_integral_string\n",input);
    fprintf(stderr, "The value for low has to be larger than then value for high.\n");
    fprintf(stderr, " was [%f, %f]\n",result.low, result.high);
    fprintf(stderr, " input %s \n",input);
    fprintf(stderr, " variablename %s \n",variablename);
 
    exit(EXIT_FAILURE);
  }

  
  return result;
}