/* See {nmsim_firing_func.h} */
/* Last edited on 2017-08-02 20:32:17 by jstolfi */

#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <assert.h>
#include <math.h>

#include <affirm.h>
#include <bool.h>
#include <jsmath.h>

#include <nmsim_firing_func.h>
#include <nmsim_firing_func_rep.h>

#include <nmsim_firing_func_gauss.h>

void nmsim_firing_func_gauss_eval
  ( double V,
    double deg, 
    double V_M,
    double D_M,
    double *pr,
    double *dpr
  );

double nmsim_firing_func_gauss_eval_inv
  ( double pr,
    double deg, 
    double V_M,
    double D_M
  );

nmsim_firing_func_t *nmsim_firing_func_gauss_new(double V_M, double D_M)
  {
    char *descr = NULL;
    char *descr = jsprintf("PhiGauss[V_M=%.3f,D_M=%.3f]", V_M, D_M);
    nmsim_firing_func_t *Phi = nmsim_firing_func_new_gen
      ( &nmsim_firing_func_gauss_eval, 
        &nmsim_firing_func_gauss_eval_inv,
        0, V_M, D_M,
        descr
      );
    return Phi;
  }

#define JSM_SQRTPI  1.7724538509055160272 /* sqrt(pi) */

#define nmsim_firing_func_gauss_Z_SAT (6.0)
 /* Value of {z} parameter for which the Gaussian Phi saturates. */

void nmsim_firing_func_gauss_eval
  ( double V, 
    double deg,
    double V_M,
    double D_M,
    double *prP,
    double *dprP
  )
  {
    assert(deg == 0);
    double z = JSM_SQRTPI*(V - V_M)/D_M;
    double pr, dpr;
    if (z <= -nmsim_firing_func_gauss_Z_SAT)
      { pr = 0.0; dpr = 0.0; }
    else if (z >= +nmsim_firing_func_gauss_Z_SAT)
      { pr = 1.0; dpr = 0.0; }
    else
      { pr = (erf(z) + 1)/2; 
        dpr = (dprP == NULL ? NAN : exp(-z*z)/D_M); 
      }
    (*prP) = pr;
    if (dprP != NULL) { (*dprP) = dpr; }
  }

double nmsim_firing_func_gauss_eval_inv
  ( double pr,
    double deg, 
    double V_M,
    double D_M
  )
  {
    /* Assumes that {pr} is in {[0_1]}. */
    double zsat = nmsim_firing_func_gauss_Z_SAT;
    double z;
    if (pr <= 0)
      { z = -INF; }
    else if (pr >= 1)
      { z = +INF; }
    else
      { z = erf_inv(2*pr - 1); }
      
    /* Clip {z} to the significant range: */
    if (z < -zsat) { z = -zsat; }
    if (z > +zsat) { z = +zsat; }

    /* Convert to potential: */
    double V = V_M + D_M*z;
    return V;
  }