text:
/*
 * Module that implements in c a Bayesian network trained with probabilistic
 * backpropagation (PBP).
 *
 * Author: Jose Miguel Hernandez Lobato
 * Date: 10 March 2015
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <malloc.h>
#include <cblas.h>

#include "network.h"
#include "pnorm.h"

/* Prototipes of private functions */

void refine_prior(NETWORK *network);
void do_ADF_update(NETWORK * network, double *x, double y);
double randn(double mu, double sigma);

/**
 * Generates a sample from a Gaussian distribution
 * @param mu    Mean of the Gaussian.
 * @param sigma Standard deviation of the Gaussian
 *
 */

double randn(double mu, double sigma) {

  double U1, U2, W, mult;
  static double X1, X2;
  static int call = 0;

  if (call == 1)
    {
      call = !call;
      return (mu + sigma * (double) X2);
    }

  do
    {
      U1 = -1 + ((double) rand () / RAND_MAX) * 2;
      U2 = -1 + ((double) rand () / RAND_MAX) * 2;
      W = pow (U1, 2) + pow (U2, 2);
    }
  while (W >= 1 || W == 0);

  mult = sqrt ((-2 * log (W)) / W);
  X1 = U1 * mult;
  X2 = U2 * mult;

  call = !call;

  return (mu + sigma * (double) X1);
}

/**
 * Constructor for a network.
 *
 * @param size_hidden_layers    The size of each hidden layer.
 * @param n_hidden_layers       The number of hidden layers.
 * @param d                     The input dimensionality.
 * @param random_noise          Random noise for the initialization of the
 *                              posterior means.
 *
 */

NETWORK *init_network(int *size_hidden_layers, int n_hidden_layers, int d,double *random_noise) {

    NETWORK *ret;
    int n, i, j, k;

    /* We save memory for the structure */

    ret = malloc(sizeof(NETWORK));

    /* We initialize the number of hidden layers */

    ret->n_hidden_layers = n_hidden_layers;

    /* We refine the prior on the noise variance and on the prior variance */

    ret->a_noise = 2.0 * 3;
    ret->b_noise = 2.0 * 3;
    ret->a_prior = 2.0 * 3;
    ret->b_prior = 2.0 * 3;

    /* We save memory for the array with the number of neurons per layer,
     * the activation type, the scaling of the activation per layer,
       and he starting position for the different arrays */

    ret->neurons_per_layer = malloc(sizeof(int) * (n_hidden_layers + 3));
    ret->linear = malloc(sizeof(int) * (n_hidden_layers + 3));
    ret->scaling_a = malloc(sizeof(double) * (n_hidden_layers + 3));
    ret->start_w = malloc(sizeof(int) * (n_hidden_layers + 3));
    ret->start_z = malloc(sizeof(int) * (n_hidden_layers + 3));
    ret->start_a = malloc(sizeof(int) * (n_hidden_layers + 3));

    /* We initialize the number of neurons per layer, the activation type and
     * compute the size of the auxiliary variables */

    /* The input layer */

    ret->start_w[ 0 ] = -1;
    ret->start_z[ 0 ] = 0;
    ret->start_a[ 0 ] = 0;

    ret->neurons_per_layer[ 0 ] = d;
    ret->linear[ 0 ] = 0;
    ret->scaling_a[ 0 ] = 1.0;

    ret->size_w = 0;
    ret->size_a = d;
    ret->size_z = d + 1;

    /* The hidden layers */

    for (i = 1 ; i < n_hidden_layers + 1 ; i ++) {
        ret->start_w[ i ] = ret->size_w;
        ret->start_z[ i ] = ret->size_z;
        ret->start_a[ i ] = ret->size_a;

        ret->neurons_per_layer[ i ] = size_hidden_layers[ i - 1 ];
        ret->linear[ i ] = 0;
        ret->scaling_a[ i ] = 1.0 / (ret->neurons_per_layer[ i - 1 ] + 1);

        ret->size_w += (ret->neurons_per_layer[ i - 1 ] + 1) * ret->neurons_per_layer[ i ];
        ret->size_a += ret->neurons_per_layer[ i ];
        ret->size_z += ret->neurons_per_layer[ i ] + 1;
    }

    /* The output layer */

    ret->start_w[ i ] = ret->size_w;
    ret->start_z[ i ] = ret->size_z;
    ret->start_a[ i ] = ret->size_a;

    ret->neurons_per_layer[ i ] = 1;
    ret->linear[ i ] = 1;
    ret->scaling_a[ i ] = 1.0 / (ret->neurons_per_layer[ i - 1 ] + 1);

    ret->size_w += (ret->neurons_per_layer[ i - 1 ] + 1) * ret->neurons_per_layer[ i ];
    ret->size_a += ret->neurons_per_layer[ i ];
    ret->size_z += ret->neurons_per_layer[ i ] + 1;

    /* The fake output layer */

    i++;
    ret->start_w[ i ] = ret->size_w;
    ret->start_z[ i ] = ret->size_z;
    ret->start_a[ i ] = ret->size_a;

    ret->neurons_per_layer[ i ] = 1;
    ret->linear[ i ] = 1;
    ret->scaling_a[ i ] = 1.0;

    ret->size_w += (ret->neurons_per_layer[ i - 1 ] + 1) * ret->neurons_per_layer[ i ];
    ret->size_a += ret->neurons_per_layer[ i ];
    ret->size_z += ret->neurons_per_layer[ i ] + 1;

    /* We save memory for the auxiliary variables */

    ret->sample_w = (double *) malloc(sizeof(double) * ret->size_w);
    ret->m_w = (double *) malloc(sizeof(double) * ret->size_w);
    ret->m_w_new = (double *) malloc(sizeof(double) * ret->size_w);
    ret->m_w_squared = (double *)  malloc(sizeof(double) * ret->size_w);
    ret->v_w = (double *) malloc(sizeof(double) * ret->size_w);
    ret->v_w_new = (double *) malloc(sizeof(double) * ret->size_w);
    ret->m_z = (double *) malloc(sizeof(double) * ret->size_z);
    ret->m_z_scaled = (double *) malloc(sizeof(double) * ret->size_z);
    ret->m_z_squared = (double *) malloc(sizeof(double) * ret->size_z);
    ret->v_z = (double *) malloc(sizeof(double) * ret->size_z);
    ret->v_z_scaled = (double *) malloc(sizeof(double) * ret->size_z);
    ret->m_a = (double *) malloc(sizeof(double) * ret->size_a);
    ret->v_a = (double *) malloc(sizeof(double) * ret->size_a);
    ret->alpha = (double *) malloc(sizeof(double) * ret->size_a);
    ret->gamma = (double *) malloc(sizeof(double) * ret->size_a);
    ret->delta_m = (double *) malloc(sizeof(double) * ret->size_z);
    ret->delta_v = (double *) malloc(sizeof(double) * ret->size_z);

    ret->dm_z_d_m_a = (double *) malloc(sizeof(double) * ret->size_z);
    ret->dv_z_d_m_a = (double *) malloc(sizeof(double) * ret->size_z);
    ret->dm_z_d_v_a = (double *) malloc(sizeof(double) * ret->size_z);
    ret->dv_z_d_v_a = (double *) malloc(sizeof(double) * ret->size_z);
    ret->dm_a_d_m_a = (double *) malloc(sizeof(double) * ret->size_w);
    ret->dv_a_d_m_a = (double *) malloc(sizeof(double) * ret->size_w);
    ret->dm_a_d_v_a = (double *) malloc(sizeof(double) * ret->size_w);
    ret->dv_a_d_v_a = (double *) malloc(sizeof(double) * ret->size_w);
    ret->dm_a_d_m_z = (double *) malloc(sizeof(double) * ret->size_w);
    ret->dv_a_d_m_z = (double *) malloc(sizeof(double) * ret->size_w);
    ret->dm_a_d_v_z = (double *) malloc(sizeof(double) * ret->size_w);
    ret->dv_a_d_v_z = (double *) malloc(sizeof(double) * ret->size_w);
    ret->dm_a_d_m_w = (double *) malloc(sizeof(double) * ret->size_w);
    ret->dv_a_d_m_w = (double *) malloc(sizeof(double) * ret->size_w);
    ret->dm_a_d_v_w = (double *) malloc(sizeof(double) * ret->size_w);
    ret->dv_a_d_v_w = (double *) malloc(sizeof(double) * ret->size_w);
    ret->grad_m_w = (double *) malloc(sizeof(double) * ret->size_w);
    ret->grad_v_w = (double *) malloc(sizeof(double) * ret->size_w);

    ret->m_w_hat_nat = (double *) malloc(sizeof(double) * ret->size_w);
    ret->v_w_hat_nat = (double *) malloc(sizeof(double) * ret->size_w);
    ret->a_w_hat_nat = (double *) malloc(sizeof(double) * ret->size_w);
    ret->b_w_hat_nat = (double *) malloc(sizeof(double) * ret->size_w);

    ret->m_w_old = (double *) malloc(sizeof(double) * ret->size_w);
    ret->v_w_old = (double *) malloc(sizeof(double) * ret->size_w);

    /* We initialize the posterior approximation */

    k = 0;
    for (i = 1 ; i < n_hidden_layers + 3 ; i ++) {
        n = (ret->neurons_per_layer[ i - 1 ] + 1) * ret->neurons_per_layer[ i ];
        for (j = 0 ; j < n ; j++) {
            ret->m_w[ k ] = 1.0 / sqrt(ret->neurons_per_layer[ i - 1 ] + 1) * random_noise[ k ];
            ret->v_w[ k ] = ret->b_prior / (ret->a_prior - 1);
            ret->m_w_hat_nat[ k ] = 0;
            ret->v_w_hat_nat[ k ] = (ret->a_prior - 1) / ret->b_prior;
            ret->a_w_hat_nat[ k ] = 0;
            ret->b_w_hat_nat[ k ] = 0;
            k++;
        }
    }

    /* We initialize the weights of the fake ouput layer */

    ret->m_w[ ret->start_w[ n_hidden_layers + 2 ] ] = 1.0;
    ret->m_w[ ret->start_w[ n_hidden_layers + 2 ] + 1 ] = 0.0;
    ret->v_w[ ret->start_w[ n_hidden_layers + 2 ] ] = 0.0;
    ret->v_w[ ret->start_w[ n_hidden_layers + 2 ] + 1 ] = 0.0;

    return ret;
}

/**
 * Constructor for a network.
 *
 * @param network   Pointer to the network to destroy.
 *
 */

void destroy_network(NETWORK *n) {

    free(n->neurons_per_layer);
    free(n->linear);
    free(n->scaling_a);
    free(n->start_w);
    free(n->start_a);
    free(n->start_z);
    free(n->sample_w);
    free(n->m_w);
    free(n->m_w_new);
    free(n->m_w_squared);
    free(n->v_w);
    free(n->v_w_new);
    free(n->m_z);
    free(n->m_z_squared);
    free(n->m_z_scaled);
    free(n->v_z);
    free(n->m_a);
    free(n->v_a);
    free(n->alpha);
    free(n->gamma);
    free(n->delta_m);
    free(n->delta_v);
    free(n->dm_z_d_m_a);
    free(n->dv_z_d_m_a);
    free(n->dm_z_d_v_a);
    free(n->dv_z_d_v_a);
    free(n->dm_a_d_m_a);
    free(n->dv_a_d_m_a);
    free(n->dm_a_d_v_a);
    free(n->dv_a_d_v_a);
    free(n->dm_a_d_m_z);
    free(n->dv_a_d_m_z);
    free(n->dm_a_d_v_z);
    free(n->dv_a_d_v_z);
    free(n->dm_a_d_m_w);
    free(n->dv_a_d_m_w);
    free(n->dm_a_d_v_w);
    free(n->dv_a_d_v_w);
    free(n->grad_m_w);
    free(n->grad_v_w);
    free(n->m_w_hat_nat);
    free(n->v_w_hat_nat);
    free(n->a_w_hat_nat);
    free(n->b_w_hat_nat);
    free(n->m_w_old);
    free(n->v_w_old);

    return;
}

/**
 * Deterministic forward propagation of the data.
 *
 */

void deterministc_forward_PBP(NETWORK *n, int index_layer) {

    int i;
    double sqrt_aux;

    /* Auxiliary pointers */

    double * __restrict__ sample_w;
    double * __restrict__ m_z;
    double * __restrict__ m_z_previous_layer;
    double * __restrict__ m_a;

    /* Auxiliary variables */

    double scaling_a;
    int n_neurons;
    int input_size;
    int linear;

    /* We initialize the pointers */

    sample_w = n->sample_w + n->start_w[ index_layer ];
    m_z = n->m_z + n->start_z[ index_layer ];
    m_z_previous_layer = n->m_z + n->start_z[ index_layer - 1 ];
    m_a = n->m_a + n->start_a[ index_layer ];

    /* We update the number of neurons, input size, type of activation and scaling constant */

    n_neurons = n->neurons_per_layer[ index_layer ];
    input_size = n->neurons_per_layer[ index_layer - 1 ] + 1;
    linear = n->linear[ index_layer ];
    scaling_a = n->scaling_a[ index_layer ];

    /* We compute the scaled activation mean: sqrt_aux * sample_w x m_z */

    sqrt_aux = sqrt(scaling_a);
    cblas_dgemv(CblasRowMajor, CblasNoTrans, n_neurons, input_size, sqrt_aux, sample_w, input_size, m_z_previous_layer, 1, 0.0, m_a, 1);

    /* We apply the non-linearity or compute a linear activation function */

    if (linear) {
        for (i = 0 ; i < n_neurons ; i++) {
            m_z[ i ] = m_a[ i ];
        }
    } else {
        for (i = 0 ; i < n_neurons ; i++) {
            m_z[ i ] = (m_a[ i ] > 0) ? m_a[ i ] : 0;
        }
    }

    /* We add the bias if we are not in the two last layers */

    if (index_layer < n->n_hidden_layers + 1) {
        m_z[ n_neurons ] = 1.0;
    }
}

/**
 * Forward propagation of probabilities.
 *
 */

void forward_PBP(NETWORK *n, int index_layer) {

    int i, limit;
    double a, a2, p, g, sqrt_aux;

    /* Auxiliary pointers */

    double * __restrict__ m_w;
    double * __restrict__ m_w_squared;
    double * __restrict__ v_w;
    double * __restrict__ m_z;
    double * __restrict__ v_z;
    double * __restrict__ m_a;
    double * __restrict__ v_a;
    double * __restrict__ m_z_previous_layer;
    double * __restrict__ m_z_previous_layer_squared;
    double * __restrict__ v_z_previous_layer;
    double * __restrict__ alpha;
    double * __restrict__ gamma;

    /* Auxiliary variables */

    double scaling_a;
    int n_neurons;
    int input_size;
    int linear;

    /* We initialize the pointers */

    m_w = n->m_w + n->start_w[ index_layer ];
    m_w_squared = n->m_w_squared + n->start_w[ index_layer ];
    v_w = n->v_w + n->start_w[ index_layer ];
    m_z = n->m_z + n->start_z[ index_layer ];
    v_z = n->v_z + n->start_z[ index_layer ];
    m_a = n->m_a + n->start_a[ index_layer ];
    v_a = n->v_a + n->start_a[ index_layer ];
    m_z_previous_layer = n->m_z + n->start_z[ index_layer - 1 ];
    v_z_previous_layer = n->v_z + n->start_z[ index_layer - 1 ];
    m_z_previous_layer_squared = n->m_z_squared + n->start_z[ index_layer - 1 ];
    alpha = n->alpha + n->start_a[ index_layer ];
    gamma = n->gamma + n->start_a[ index_layer ];

    /* We update the number of neurons, input size, type of activation and scaling constant */

    n_neurons = n->neurons_per_layer[ index_layer ];
    input_size = n->neurons_per_layer[ index_layer - 1 ] + 1;
    linear = n->linear[ index_layer ];
    scaling_a = n->scaling_a[ index_layer ];

    /* We compute the m_z_previous_layer_squared and the m_w_squared */

    limit = n_neurons * input_size;
    for (i = 0 ; i < limit ; i++)
        m_w_squared[ i ] = m_w[ i ] * m_w[ i ];
    for (i = 0 ; i < input_size ; i++)
        m_z_previous_layer_squared[ i ] = m_z_previous_layer[ i ] * m_z_previous_layer[ i ];

    /* We compute the scaled activation mean: sqrt_aux * m_w x m_z */

    sqrt_aux = sqrt(scaling_a);
    cblas_dgemv(CblasRowMajor, CblasNoTrans, n_neurons, input_size, sqrt_aux, m_w, input_size, m_z_previous_layer, 1, 0.0, m_a, 1);

    /* We compute the scaled activation variance: (m_w_squared x v_z_previous_layer +
       v_w x m_z_previous_layer_squared + v_w x v_z_previous_layer) */

    cblas_dgemv(CblasRowMajor, CblasNoTrans, n_neurons, input_size, scaling_a, m_w_squared, input_size, v_z_previous_layer, 1, 0, v_a, 1);
    cblas_dgemv(CblasRowMajor, CblasNoTrans, n_neurons, input_size, scaling_a, v_w, input_size, m_z_previous_layer_squared, 1, 1.0, v_a, 1);
    cblas_dgemv(CblasRowMajor, CblasNoTrans, n_neurons, input_size, scaling_a, v_w, input_size, v_z_previous_layer, 1, 1.0, v_a, 1);

    /* We apply the non-linearity or compute a linear activation function */

    if (linear) {
        for (i = 0 ; i < n_neurons ; i++) {
            m_z[ i ] = m_a[ i ];
            v_z[ i ] = v_a[ i ];
        }
    } else {
        for (i = 0 ; i < n_neurons ; i++) {

            sqrt_aux = sqrt(v_a[ i ]);
            a = m_a[ i ] / sqrt_aux;
            a2 = a * a;
            p = pnorm(a);
            g = (a < -30) ? -a - 1.0 / a + 2.0 / (a * a2) : 0.398942280401432677939946059 * exp(-0.5 * a2) / p;

            /* We compute the output mean and variance */

            m_z[ i ] = p * (m_a[ i ] + sqrt_aux * g);
            v_z[ i ] = m_z[ i ] * (m_a[ i ] + sqrt_aux * g) * (1 - p) + p * v_a[ i ] * (1 - g * g - g * a);

            alpha[ i ] = a;
            gamma[ i ] = g;
        }
    }

    /* We add the bias if we are not in the two last layers */

    if (index_layer < n->n_hidden_layers + 1) {
        m_z[ n_neurons ] = 1.0;
        v_z[ n_neurons ] = 0.0;
    }
}

/**
 * Backward computation of gradients.
 *
 */

void backward_PBP(NETWORK *n, int index_layer) {

    double g, a, da_dm_a, da_dv_a, dg_dm_a, dg_dv_a, d, p;
    int i, j, k, limit_1, limit_2;

    /* Auxiliary pointers */

    double * __restrict__ m_w;
    double * __restrict__ m_w_layer_above;
    double * __restrict__ v_w_layer_above;
    double * __restrict__ m_z;
    double * __restrict__ m_z_scaled;
    double * __restrict__ m_a;
    double * __restrict__ v_a;
    double * __restrict__ m_z_layer_below;
    double * __restrict__ m_z_layer_below_scaled;
    double * __restrict__ v_z_layer_below;
    double * __restrict__ v_z_layer_below_scaled;
    double * __restrict__ alpha;
    double * __restrict__ gamma;
    double * __restrict__ delta_m_layer_above;
    double * __restrict__ delta_v_layer_above;
    double * __restrict__ delta_m;
    double * __restrict__ delta_v;
    double * __restrict__ dm_z_d_m_a;
    double * __restrict__ dv_z_d_m_a;
    double * __restrict__ dm_z_d_v_a;
    double * __restrict__ dv_z_d_v_a;
    double * __restrict__ dm_a_d_m_a;
    double * __restrict__ dv_a_d_m_a;
    double * __restrict__ dm_a_d_v_a;
    double * __restrict__ dv_a_d_v_a;
    double * __restrict__ dm_a_d_m_z;
    double * __restrict__ dv_a_d_m_z;
    double * __restrict__ dm_a_d_v_z;
    double * __restrict__ dv_a_d_v_z;
    double * __restrict__ dm_a_d_m_w;
    double * __restrict__ dv_a_d_m_w;
    double * __restrict__ dm_a_d_v_w;
    double * __restrict__ dv_a_d_v_w;
    double * __restrict__ grad_m_w;
    double * __restrict__ grad_v_w;

    /* Auxiliary variables */

    double scaling_a, scaling_a_layer_above, sqrt_aux, sqrt_v_a;
    double *p_aux_1, *p_aux_2, *p_aux_3, *p_aux_4;
    double *p_aux_1_2, *p_aux_2_2, *p_aux_3_2, *p_aux_4_2;
    int n_neurons;
    int linear;

    /* We initialize the auxiliary pointers */

    m_w = n->m_w + n->start_w[ index_layer ];
    m_w_layer_above = n->m_w + n->start_w[ index_layer + 1 ];
    v_w_layer_above = n->v_w + n->start_w[ index_layer + 1 ];
    m_z = n->m_z + n->start_z[ index_layer ];
    m_z_scaled = n->m_z_scaled + n->start_z[ index_layer ];
    m_a = n->m_a + n->start_a[ index_layer ];
    v_a = n->v_a + n->start_a[ index_layer ];
    m_z_layer_below = n->m_z + n->start_z[ index_layer - 1 ];
    v_z_layer_below = n->v_z + n->start_z[ index_layer - 1 ];
    m_z_layer_below_scaled = n->m_z_scaled + n->start_z[ index_layer - 1 ];
    v_z_layer_below_scaled = n->v_z_scaled + n->start_z[ index_layer - 1 ];
    alpha = n->alpha + n->start_a[ index_layer ];
    gamma = n->gamma + n->start_a[ index_layer ];
    delta_m_layer_above = n->delta_m + n->start_z[ index_layer + 1 ];
    delta_v_layer_above = n->delta_v + n->start_z[ index_layer + 1 ];
    delta_m = n->delta_m + n->start_z[ index_layer ];
    delta_v = n->delta_v + n->start_z[ index_layer ];

    dm_z_d_m_a = n->dm_z_d_m_a + n->start_z[ index_layer ];
    dv_z_d_m_a = n->dv_z_d_m_a + n->start_z[ index_layer ];
    dm_z_d_v_a = n->dm_z_d_v_a + n->start_z[ index_layer ];
    dv_z_d_v_a = n->dv_z_d_v_a + n->start_z[ index_layer ];
    dm_a_d_m_a = n->dm_a_d_m_a + n->start_w[ index_layer ];
    dv_a_d_m_a = n->dv_a_d_m_a + n->start_w[ index_layer ];
    dm_a_d_v_a = n->dm_a_d_v_a + n->start_w[ index_layer ];
    dv_a_d_v_a = n->dv_a_d_v_a + n->start_w[ index_layer ];
    dm_a_d_m_z = n->dm_a_d_m_z + n->start_w[ index_layer ];
    dv_a_d_m_z = n->dv_a_d_m_z + n->start_w[ index_layer ];
    dm_a_d_v_z = n->dm_a_d_v_z + n->start_w[ index_layer ];
    dv_a_d_v_z = n->dv_a_d_v_z + n->start_w[ index_layer ];
    dm_a_d_m_w = n->dm_a_d_m_w + n->start_w[ index_layer ];
    dv_a_d_m_w = n->dv_a_d_m_w + n->start_w[ index_layer ];
    dm_a_d_v_w = n->dm_a_d_v_w + n->start_w[ index_layer ];
    dv_a_d_v_w = n->dv_a_d_v_w + n->start_w[ index_layer ];
    grad_m_w = n->grad_m_w + n->start_w[ index_layer ];
    grad_v_w = n->grad_v_w + n->start_w[ index_layer ];

    /* We update the number of neurons, input size, type of activation and scaling constant */

    n_neurons = n->neurons_per_layer[ index_layer ];
    linear = n->linear[ index_layer ];
    scaling_a = n->scaling_a[ index_layer ];
    scaling_a_layer_above = n->scaling_a[ index_layer + 1 ];

    /* We compute the gradient of the non-linear activations with respect to the activations */

    if (linear) {
        for (i = 0 ; i < n_neurons ; i++) {
            dm_z_d_m_a[ i ] = 1.0;
            dm_z_d_v_a[ i ] = 0.0;
            dv_z_d_m_a[ i ] = 0.0;
            dv_z_d_v_a[ i ] = 1.0;

        }
    } else {
        for (i = 0 ; i < n_neurons ; i++) {
            g = gamma[ i ];
            a = alpha[ i ];
            da_dm_a = 1.0 / sqrt(v_a[ i ]);
            da_dv_a = m_a[ i ] / (2 * v_a[ i ]) * da_dm_a;
            if (a < -30) {
                g = -a - 1.0 / a + 2.0 / (a * a * a);
                dg_dm_a = -da_dm_a + 1.0 / (a * a) * da_dm_a - 6.0 / (a * a * a * a) * da_dm_a;
                dg_dv_a = -da_dv_a + 1.0 / (a * a) * da_dv_a - 6.0 / (a * a * a * a) * da_dv_a;
            } else {
                dg_dm_a = -(g * a + g * g) * da_dm_a;
                dg_dv_a = -(g * a + g * g) * da_dv_a;
            }
            d = 0.398942280401432677939946059 * exp(-0.5 * a * a);
            p = pnorm(a);
            sqrt_v_a = sqrt(v_a[ i ]);
            dm_z_d_m_a[ i ] = da_dm_a * d * (m_a[ i ] + sqrt_v_a * g) + p * (1 + sqrt_v_a * dg_dm_a);
            dm_z_d_v_a[ i ] = da_dv_a * d * (m_a[ i ] + sqrt_v_a * g) + p * (g / (2 * sqrt_v_a) + sqrt_v_a * dg_dv_a);
            dv_z_d_m_a[ i ] = dm_z_d_m_a[ i ] * (m_a[ i ] + sqrt_v_a * g) * (1 - p) + m_z[ i ] * ((1 + sqrt_v_a * dg_dm_a) * (1 - p) -
                (m_a[ i ] + sqrt_v_a * g) * d * da_dm_a) + d * da_dm_a * v_a[ i ] * (1 - g * g - g * a) - p * v_a[ i ] *
                (2 * g * dg_dm_a + dg_dm_a * a + g * da_dm_a);
            dv_z_d_v_a[ i ] = dm_z_d_v_a[ i ] * (m_a[ i ] + sqrt_v_a * g) * (1 - p) + m_z[ i ] * ((0.5 / sqrt_v_a * g + dg_dv_a * sqrt_v_a) * (1 - p) -
                (m_a[ i ]+ sqrt_v_a * g) * d * da_dv_a) + d * da_dv_a * v_a[ i ] * (1 - g * g - g * a) + p * ((1 - g * g - g * a) + v_a[ i ] *
                (-2 * g * dg_dv_a - dg_dv_a * a - g * da_dv_a));
        }
    }

    /* We scale m_z */

    limit_1 = n_neurons + 1;
    for (i = 0 ; i < limit_1 ; i++) {
        m_z_scaled[ i ] = m_z[ i ] * scaling_a_layer_above;
    }

    /* We initialize the rows of dv_a_dm_z to m_z_scaled */

    limit_1 = n->neurons_per_layer[ index_layer + 1 ];
    limit_2 = n_neurons + 1;
    for (i = 0 ; i < limit_1 ; i++) {
        p_aux_1 = dv_a_d_m_z + i * limit_2;
        for (j = 0 ; j < limit_2 ; j++) {
            p_aux_1[ j ] = m_z_scaled[ j ];
        }
    }

    /* We compute the gradient of the activation of the layer above with respect to the non-linear activations in this layer */

    sqrt_aux = sqrt(scaling_a_layer_above);
    limit_1 = n->neurons_per_layer[ index_layer + 1 ] * (n_neurons + 1);
    for (i = 0 ; i < limit_1 ; i++) {
        dm_a_d_m_z[ i ] = m_w_layer_above[ i ] * sqrt_aux;
        dm_a_d_v_z[ i ] = 0;
        dv_a_d_m_z[ i ] = dv_a_d_m_z[ i ] * 2 * v_w_layer_above[ i ];
        dv_a_d_v_z[ i ] = (m_w_layer_above[ i ] * m_w_layer_above[ i ] + v_w_layer_above[ i ]) * scaling_a_layer_above;
    }

    /* We scale the activatons from the layer below */

    sqrt_aux = sqrt(scaling_a);
    limit_1 = n->neurons_per_layer[ index_layer - 1 ] + 1;
    for (i = 0 ; i < limit_1 ; i++) {
        m_z_layer_below_scaled[ i ] = m_z_layer_below[ i ] * sqrt_aux;
        v_z_layer_below_scaled[ i ] = v_z_layer_below[ i ] * scaling_a;
    }

    /* We initialize dm_a_dm_w, dv_a_dm_w and dv_a_dv_w */

    limit_1 = n_neurons;
    limit_2 = n->neurons_per_layer[ index_layer - 1 ] + 1;
    for (i = 0 ; i < limit_1 ; i++) {
        p_aux_1 = dm_a_d_m_w + i * limit_2;
        p_aux_2 = dv_a_d_m_w + i * limit_2;
        p_aux_3 = dv_a_d_v_w + i * limit_2;
        for (j = 0 ; j < limit_2 ; j++) {
            p_aux_1[ j ] = m_z_layer_below_scaled[ j ];
            p_aux_2[ j ] = v_z_layer_below_scaled[ j ];
            p_aux_3[ j ] = (m_z_layer_below_scaled[ j ] * m_z_layer_below_scaled[ j ] + v_z_layer_below_scaled[ j ]);
        }
    }

    limit_1 = n_neurons * (n->neurons_per_layer[ index_layer - 1 ] + 1);
    for (i = 0 ; i < limit_1 ; i++) {
         dm_a_d_v_w[ i ] = 0;
         dv_a_d_m_w[ i ] = 2 * dv_a_d_m_w[ i ] * m_w[ i ];
    }

    /* We compute the gradient of the activations of the top layer with respect to the activations of the current layer */

    limit_1 = n->neurons_per_layer[ index_layer + 1 ];
    limit_2 = n_neurons + 1;
    k = 0;
    for (i = 0 ; i < limit_1 ; i++) {
        p_aux_1 = dm_a_d_m_a + i * limit_2;
        p_aux_2 = dv_a_d_m_a + i * limit_2;
        p_aux_3 = dm_a_d_v_a + i * limit_2;
        p_aux_4 = dv_a_d_v_a + i * limit_2;
        p_aux_1_2 = dm_a_d_m_z + i * limit_2;
        p_aux_2_2 = dm_a_d_v_z + i * limit_2;
        p_aux_3_2 = dv_a_d_m_z + i * limit_2;
        p_aux_4_2 = dv_a_d_v_z + i * limit_2;
        for (j = 0 ; j < limit_2 - 1 ; j++) {
            p_aux_1[ j ] = p_aux_1_2[ j ] * dm_z_d_m_a[ j ] + p_aux_2_2[ j ] * dv_z_d_m_a[ j ];
            p_aux_2[ j ] = p_aux_3_2[ j ] * dm_z_d_m_a[ j ] + p_aux_4_2[ j ] * dv_z_d_m_a[ j ];
            p_aux_3[ j ] = p_aux_1_2[ j ] * dm_z_d_v_a[ j ] + p_aux_2_2[ j ] * dv_z_d_v_a[ j ];
            p_aux_4[ j ] = p_aux_3_2[ j ] * dm_z_d_v_a[ j ] + p_aux_4_2[ j ] * dv_z_d_v_a[ j ];
        }
    }

    /* We compute the deltas */

    limit_1 = n->neurons_per_layer[ index_layer + 1 ];
    cblas_dgemv(CblasRowMajor, CblasTrans, limit_1, n_neurons + 1, 1.0, dm_a_d_m_a, n_neurons + 1, delta_m_layer_above, 1, 0.0, delta_m, 1);
    cblas_dgemv(CblasRowMajor, CblasTrans, limit_1, n_neurons + 1, 1.0, dv_a_d_m_a, n_neurons + 1, delta_v_layer_above, 1, 1.0, delta_m, 1);
    cblas_dgemv(CblasRowMajor, CblasTrans, limit_1, n_neurons + 1, 1.0, dm_a_d_v_a, n_neurons + 1, delta_m_layer_above, 1, 0.0, delta_v, 1);
    cblas_dgemv(CblasRowMajor, CblasTrans, limit_1, n_neurons + 1, 1.0, dv_a_d_v_a, n_neurons + 1, delta_v_layer_above, 1, 1.0, delta_v, 1);

    /* We compute the gradients */

    limit_1 = n_neurons;
    limit_2 = n->neurons_per_layer[ index_layer - 1 ] + 1;
    k = 0;
    for (i = 0 ; i < limit_1 ; i++) {
        for (j = 0 ; j < limit_2 ; j++) {
            grad_m_w[ k ] = delta_m[ i ] * dm_a_d_m_w[ k ] + delta_v[ i ] * dv_a_d_m_w[ k ];
            grad_v_w[ k ] = delta_m[ i ] * dm_a_d_v_w[ k ] + delta_v[ i ] * dv_a_d_v_w[ k ];
            k++;
        }
    }
}

/**
 * Function that performs an ADF update.
 *
 * @param network   Pointer to the network.
 * @param x         Pointer to the input features.
 * @param y         Pointer to the target.
 *
 */

void do_ADF_update(NETWORK * n, double *x, double y) {

    int i, limit_1;

    double * __restrict__ m_z;
    double * __restrict__ v_z;
    double * __restrict__ delta_m;
    double * __restrict__ delta_v;

    double * __restrict__ m_w;
    double * __restrict__ v_w;
    double * __restrict__ grad_m_w;
    double * __restrict__ grad_v_w;

    double m, v, v1, v2, logZ, logZ1, logZ2, d_logZ_d_m, d_logZ_d_v, a_new, b_new, m_aux, v_aux;

    /* We initialize the non-linear activations of the input layer with the data */

    limit_1 = n->start_a[ 1 ];
    m_z = n->m_z;
    v_z = n->v_z;
    for (i = 0 ; i < limit_1 ; i++) {
        m_z[ i ] = x[ i ];
        v_z[ i ] = 0.0;
    }
    m_z[ i ] = 1.0;
    v_z[ i ] = 0.0;

    /* We do a forward pass */

    limit_1 = n->n_hidden_layers + 2;
    for (i = 1 ; i < limit_1 ; i++)
        forward_PBP(n, i);

    /* We obtain logZ, logZ1 and logZ2 and the updates for a and b */

    m_z = n->m_z + n->start_z[ n->n_hidden_layers + 1 ];
    v_z = n->v_z + n->start_z[ n->n_hidden_layers + 1 ];

    m = m_z[ 0 ];
    v = v_z[ 0 ] + n->b_noise / (n->a_noise - 1);
    v1 = v_z[ 0 ] + n->b_noise / (n->a_noise - 0);
    v2 = v_z[ 0 ] + n->b_noise / (n->a_noise + 1);
    logZ = -0.5 * (log(v) + (m - y) * (m - y) / v);
    logZ1 = -0.5 * (log(v1) + (m - y) * (m - y) / v1);
    logZ2 = -0.5 * (log(v2) + (m - y) * (m - y) / v2);

    a_new = 1.0 / (exp(logZ2 - 2 * logZ1 + logZ) * (n->a_noise + 1.0) / n->a_noise - 1.0);
    b_new = 1.0 / (exp(logZ2 - logZ1) * (n->a_noise + 1) / n->b_noise - exp(logZ1 - logZ) * n->a_noise / n->b_noise);

    /* We initialize the deltas for the output layer */

    d_logZ_d_m = -(m - y) / v;
    d_logZ_d_v = -0.5 / v + 0.5 * (m - y) * (m - y) / (v * v);
    delta_m = n->delta_m + n->start_z[ n->n_hidden_layers + 2 ];
    delta_v = n->delta_v + n->start_z[ n->n_hidden_layers + 2 ];
    delta_m[ 0 ] = d_logZ_d_m;
    delta_v[ 0 ] = d_logZ_d_v;

    /* We do a backward pass */

    limit_1 = n->n_hidden_layers + 1;
    for (i = limit_1 ; i >= 1 ; i--)
        backward_PBP(n, i);

    /* We update the mean and variance parameters, except those in the fake output layer */

    m_w = n->m_w;
    v_w = n->v_w;
    grad_m_w = n->grad_m_w;
    grad_v_w = n->grad_v_w;

    limit_1 = n->size_w - 2;
    for (i = 0 ; i < limit_1 ; i++) {
        v_aux = v_w[ i ] - v_w[ i ] * v_w[ i ] * (grad_m_w[ i ] * grad_m_w[ i ] - 2 * grad_v_w[ i ]);
        m_aux = m_w[ i ] + v_w[ i ] * grad_m_w[ i ];
        if (v_aux > 1e-100 && m_aux != NAN && v_aux != NAN) {
            m_w[ i ] = m_aux;
            v_w[ i ] = v_aux;
        }
    }

    /* We update the noise variables */

    n->b_noise = b_new;
    n->a_noise = a_new;
}

/**
 * Function that refines the approximate factor for the prior.
 *
 * @param network   Pointer to the network.
 *
 */

void refine_prior(NETWORK *n) {

    int i;

    double v_w_nat, m_w_nat, v_w_cav_nat, m_w_cav_nat, v_w_cav, m_w_cav,
        a_w_nat, b_w_nat, a_w_cav_nat, b_w_cav_nat, a_w_cav, b_w_cav,
        v, v1, v2, logZ, logZ1, logZ2, d_logZ_d_m_w_cav, d_logZ_d_v_w_cav,
        m_w_new, v_w_new, a_w_new, b_w_new, v_w_new_nat, m_w_new_nat,
        a_w_new_nat, b_w_new_nat;

    /* We iterate over the layers refining the prior */

    for (i = 0 ; i < n->size_w - 2 ; i++) {
        v_w_nat = 1.0 / n->v_w[ i ];
        m_w_nat = n->m_w[ i ] / n->v_w[ i ];
        v_w_cav_nat = v_w_nat - n->v_w_hat_nat[ i ];
        m_w_cav_nat = m_w_nat - n->m_w_hat_nat[ i ];

        v_w_cav = 1.0 / v_w_cav_nat;
        m_w_cav = m_w_cav_nat / v_w_cav_nat;
        a_w_nat = n->a_prior - 1;
        b_w_nat = -n->b_prior;
        a_w_cav_nat = a_w_nat - n->a_w_hat_nat[ i ];
        b_w_cav_nat = b_w_nat - n->b_w_hat_nat[ i ];
        a_w_cav = a_w_cav_nat + 1;
        b_w_cav = -b_w_cav_nat;

        if (v_w_cav > 0 && b_w_cav > 0 && a_w_cav > 1 && v_w_cav < 1e6) {

            v = v_w_cav + b_w_cav / (a_w_cav - 1);
            v1  = v_w_cav + b_w_cav / a_w_cav;
            v2  = v_w_cav + b_w_cav / (a_w_cav + 1);
            logZ = -0.5 * log(v) - 0.5 * m_w_cav * m_w_cav / v;
            logZ1 = -0.5 * log(v1) - 0.5 * m_w_cav * m_w_cav / v1;
            logZ2 = -0.5 * log(v2) - 0.5 * m_w_cav * m_w_cav / v2;
            d_logZ_d_m_w_cav = -m_w_cav / v;
            d_logZ_d_v_w_cav = -0.5 / v + 0.5 * m_w_cav * m_w_cav / (v * v);
            m_w_new = m_w_cav + v_w_cav * d_logZ_d_m_w_cav;
            v_w_new = v_w_cav - v_w_cav * v_w_cav * (d_logZ_d_m_w_cav * d_logZ_d_m_w_cav - 2 * d_logZ_d_v_w_cav) ;
            a_w_new = 1.0 / (exp(logZ2 - 2 * logZ1 + logZ) * (a_w_cav + 1) / a_w_cav - 1.0);
            b_w_new = 1.0 / (exp(logZ2 - logZ1) * (a_w_cav + 1) / (b_w_cav) - exp(logZ1 - logZ) * a_w_cav / b_w_cav);
            v_w_new_nat = 1.0 / v_w_new;
            m_w_new_nat = m_w_new / v_w_new;
            a_w_new_nat = a_w_new - 1;
            b_w_new_nat = -b_w_new;

            n->m_w_hat_nat[ i ] = m_w_new_nat - m_w_cav_nat;
            n->v_w_hat_nat[ i ] = v_w_new_nat - v_w_cav_nat;
            n->a_w_hat_nat[ i ] = a_w_new_nat - a_w_cav_nat;
            n->b_w_hat_nat[ i ] = b_w_new_nat - b_w_cav_nat;

            n->m_w[ i ] = m_w_new;
            n->v_w[ i ] = v_w_new;

            n->a_prior = a_w_new;
            n->b_prior = b_w_new;
        }
    }
}

/**
 * Function that performs one learning epoch.
 *
 * @param network       Pointer to the network.
 * @param x             Pointer to the input feature matrix.
 * @param y             Pointer to the target vector.
 * @param n_datapoints  Number of data points in the training set.
 * @param d             Dimensionality of the training data.
 * @param permutation   A random permuation of the input data.
 *
 */

NETWORK *one_learning_epoch(NETWORK *n, double *x, double *y, int n_datapoints, int d, int *permutation) {

    int i, index;

    /* We do one ADF upte for each datapoint */

    for (i = 0 ; i < n_datapoints ; i++) {
        index = permutation[ i ] * d;
        do_ADF_update(n, x + index, y[ permutation[ i ] ]);

        if (i % 1000 == 0) {
            printf(".");
            fflush(stdout);
        }
    }
    printf("\n");
    fflush(stdout);

    /* We refine the prior */

    refine_prior(n);

    /* We are done */

    return n;
}

/**
 * Function that predicts deterministically.
 *
 * @param network       Pointer to the network.
 * @param x_test        Pointer to the input feature matrix.
 * @param y_test        Pointer to store the predictions.
 * @param n_datapoints  Number of data points in the test set.
 * @param d             Dimensionality of the test data.
 *
 */

NETWORK *predict_deterministic(NETWORK *n, double *x_test, double *y_test, int n_datapoints, int d) {

    int i, j, index, limit_1;

    double *m_z;

    for (i = 0 ; i < n_datapoints ; i++) {

        /* We initialize the non-linear activations of the input layer with the data */

        limit_1 = n->start_a[ 1 ];
        m_z = n->m_z;
        for (j = 0 ; j < limit_1 ; j++) {
            index = i * d + j;
            m_z[ j ] = x_test[ index ];
        }
        m_z[ j ] = 1.0;

        /* We do a forward pass */

        limit_1 = n->n_hidden_layers + 2;
        for (j = 1 ; j < limit_1 ; j++)
            deterministc_forward_PBP(n, j);

        /* We store the results */

        m_z = n->m_z + n->start_z[ n->n_hidden_layers + 1 ];

        y_test[ i ] = m_z[ 0 ];
    }

    /* We are done */

    return n;
}

/**
 * Function that predicts.
 *
 * @param network       Pointer to the network.
 * @param x_test        Pointer to the input feature matrix.
 * @param m_test        Pointer to store the predictive mean.
 * @param v_test        Pointer to store the predictive variance.
 * @param v_noise       Pointer to store the noise variance.
 * @param n_datapoints  Number of data points in the test set.
 * @param d             Dimensionality of the test data.
 *
 */

NETWORK *predict(NETWORK *n, double *x_test, double *m_test, double *v_test, double *v_noise, int n_datapoints, int d) {

    int i, j, index, limit_1;

    double *m_z, *v_z;

    for (i = 0 ; i < n_datapoints ; i++) {

        /* We initialize the non-linear activations of the input layer with the data */

        limit_1 = n->start_a[ 1 ];
        m_z = n->m_z;
        v_z = n->v_z;
        for (j = 0 ; j < limit_1 ; j++) {
            index = i * d + j;
            m_z[ j ] = x_test[ index ];
            v_z[ j ] = 0.0;
        }
        m_z[ j ] = 1.0;
        v_z[ j ] = 0.0;

        /* We do a forward pass */

        limit_1 = n->n_hidden_layers + 2;
        for (j = 1 ; j < limit_1 ; j++)
            forward_PBP(n, j);

        /* We store the results */

        m_z = n->m_z + n->start_z[ n->n_hidden_layers + 1 ];
        v_z = n->v_z + n->start_z[ n->n_hidden_layers + 1 ];

        m_test[ i ] = m_z[ 0 ];
        v_test[ i ] = v_z[ 0 ];
    }

    /* We store the noise value */

    *v_noise = n->b_noise / (n->a_noise - 1);

    /* We are done */

    return n;
}

/* Function that returns the parameters of the network */

void get_params(NETWORK *n, double *sample_w_out, double *m_w_out, double *v_w_out, double *m_w_hat_nat_out, double *v_w_hat_nat_out,
    double *a_w_hat_nat_out, double *b_w_hat_nat_out, double *a_noise_out, double *b_noise_out, double *a_prior_out,
    double *b_prior_out, int *neurons_per_layer_out) {

    int i, limit_1;

    /* We copy the parameters */

    limit_1 = n->size_w - 2;
    for (i = 0 ; i < limit_1 ; i++) {
        sample_w_out[ i ] = n->sample_w[ i ];
        m_w_out[ i ] = n->m_w[ i ];
        v_w_out[ i ] = n->v_w[ i ];
        m_w_hat_nat_out[ i ] = n->m_w_hat_nat[ i ];
        v_w_hat_nat_out[ i ] = n->v_w_hat_nat[ i ];
        a_w_hat_nat_out[ i ] = n->a_w_hat_nat[ i ];
        b_w_hat_nat_out[ i ] = n->b_w_hat_nat[ i ];
    }

    *a_noise_out = n->a_noise;
    *b_noise_out = n->b_noise;
    *a_prior_out = n->a_prior;
    *b_prior_out = n->b_prior;

    limit_1 = n->n_hidden_layers + 2;
    for (i = 0 ; i < limit_1 ; i++) {
        neurons_per_layer_out[ i ] = n->neurons_per_layer[ i ];
    }
}

/* Function that returns the sizes of the parameters of the network */

void get_size_params(NETWORK *n, int *size_weihts, int *n_layers) {

    *size_weihts = n->size_w - 2;
    *n_layers = n->n_hidden_layers + 2;
}

/* Function that sets the parameters of the network */

void set_params(NETWORK *n, double *sample_w_in, double *m_w_in, double *v_w_in, double *m_w_hat_nat_in, double *v_w_hat_nat_in,
    double *a_w_hat_nat_in, double *b_w_hat_nat_in, double a_noise_in, double b_noise_in, double a_prior_in, double b_prior_in) {

    int i, limit_1;

    /* We copy the parameters */

    limit_1 = n->size_w - 2;
    for (i = 0 ; i < limit_1 ; i++) {
        n->sample_w[ i ] = sample_w_in[ i ];
        n->m_w[ i ] = m_w_in[ i ];
        n->v_w[ i ] = v_w_in[ i ];
        n->m_w_hat_nat[ i ] = m_w_hat_nat_in[ i ];
        n->v_w_hat_nat[ i ] = v_w_hat_nat_in[ i ];
        n->a_w_hat_nat[ i ] = a_w_hat_nat_in[ i ];
        n->b_w_hat_nat[ i ] = b_w_hat_nat_in[ i ];
    }

    n->a_noise = a_noise_in;
    n->b_noise = b_noise_in;
    n->a_prior = a_prior_in;
    n->b_prior = b_prior_in;
}


text:
/* histogram/gsl_histogram.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
 *
 * This program 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 2 of the License, or (at
 * your option) any later version.
 *
 * This program 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_HISTOGRAM_H__
#define __GSL_HISTOGRAM_H__

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_types.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct {
  size_t n ;
  double * range ;
  double * bin ;
} gsl_histogram ;

typedef struct {
  size_t n ;
  double * range ;
  double * sum ;
} gsl_histogram_pdf ;

GSL_EXPORT gsl_histogram * gsl_histogram_alloc (size_t n);

GSL_EXPORT gsl_histogram * gsl_histogram_calloc (size_t n);
GSL_EXPORT gsl_histogram * gsl_histogram_calloc_uniform (const size_t n, const double xmin, const double xmax);
GSL_EXPORT void gsl_histogram_free (gsl_histogram * h);
GSL_EXPORT int gsl_histogram_increment (gsl_histogram * h, double x);
GSL_EXPORT int gsl_histogram_accumulate (gsl_histogram * h, double x, double weight);
GSL_EXPORT int gsl_histogram_find (const gsl_histogram * h,
                                   const double x, size_t * i);

GSL_EXPORT double gsl_histogram_get (const gsl_histogram * h, size_t i);
GSL_EXPORT int gsl_histogram_get_range (const gsl_histogram * h, size_t i,
                                        double * lower, double * upper);

GSL_EXPORT double gsl_histogram_max (const gsl_histogram * h);
GSL_EXPORT double gsl_histogram_min (const gsl_histogram * h);
GSL_EXPORT size_t gsl_histogram_bins (const gsl_histogram * h);

GSL_EXPORT void gsl_histogram_reset (gsl_histogram * h);

GSL_EXPORT gsl_histogram * gsl_histogram_calloc_range(size_t n, double * range);

GSL_EXPORT
int
gsl_histogram_set_ranges (gsl_histogram * h, const double range[], size_t size);
GSL_EXPORT
int
gsl_histogram_set_ranges_uniform (gsl_histogram * h, double xmin, double xmax);

GSL_EXPORT
int
gsl_histogram_memcpy(gsl_histogram * dest, const gsl_histogram * source);

GSL_EXPORT
gsl_histogram *
gsl_histogram_clone(const gsl_histogram * source);

GSL_EXPORT double gsl_histogram_max_val (const gsl_histogram * h);

GSL_EXPORT size_t gsl_histogram_max_bin (const gsl_histogram * h);

GSL_EXPORT double gsl_histogram_min_val (const gsl_histogram * h);

GSL_EXPORT size_t gsl_histogram_min_bin (const gsl_histogram * h);

GSL_EXPORT
int
gsl_histogram_equal_bins_p(const gsl_histogram *h1, const gsl_histogram *h2);

GSL_EXPORT
int
gsl_histogram_add(gsl_histogram *h1, const gsl_histogram *h2);

GSL_EXPORT
int
gsl_histogram_sub(gsl_histogram *h1, const gsl_histogram *h2);

GSL_EXPORT
int
gsl_histogram_mul(gsl_histogram *h1, const gsl_histogram *h2);

GSL_EXPORT
int
gsl_histogram_div(gsl_histogram *h1, const gsl_histogram *h2);

GSL_EXPORT
int
gsl_histogram_scale(gsl_histogram *h, double scale);

GSL_EXPORT
int
gsl_histogram_shift (gsl_histogram * h, double shift);

GSL_EXPORT double gsl_histogram_sigma (const gsl_histogram * h);

GSL_EXPORT double gsl_histogram_mean (const gsl_histogram * h);

GSL_EXPORT double gsl_histogram_sum (const gsl_histogram * h);

GSL_EXPORT int gsl_histogram_fwrite (FILE * stream, const gsl_histogram * h) ;
GSL_EXPORT int gsl_histogram_fread (FILE * stream, gsl_histogram * h);
GSL_EXPORT int gsl_histogram_fprintf (FILE * stream, const gsl_histogram * h,
                                      const char * range_format, const char * bin_format);
GSL_EXPORT int gsl_histogram_fscanf (FILE * stream, gsl_histogram * h);

GSL_EXPORT gsl_histogram_pdf * gsl_histogram_pdf_alloc (const size_t n);
GSL_EXPORT int gsl_histogram_pdf_init (gsl_histogram_pdf * p, const gsl_histogram * h);
GSL_EXPORT void gsl_histogram_pdf_free (gsl_histogram_pdf * p);
GSL_EXPORT double gsl_histogram_pdf_sample (const gsl_histogram_pdf * p, double r);

__END_DECLS

#endif /* __GSL_HISTOGRAM_H__ */

ext: h

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#include <stdio.h>
#include <string.h>
#include <string>
#include <time.h>
#include <math.h>
#include <vector>
#include <algorithm>
#include <iterator>
#include <sys/stat.h>
//#define GSL_DLL
//#include <gsl/gsl_cdf.h> //use gaussian p-value function
using namespace std;
void printProgress(double percentage);


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

typedef struct {
    double l, g;
} Params;

int func (double t, const double y[], double f[], void *params) {
    Params *p = (Params *) params;
    double l = p->l;
    double g = p->g;
    f[0] = y[1];
    f[1] = -(g/l)*y[0];
    return GSL_SUCCESS;
}

int jac (double t, const double y[], double *dfdy, double dfdt[],
         void *params) {
    Params *p = (Params *) params;
    double l = p->l;
    double g = p->g;
    gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);
    gsl_matrix * m = &(dfdy_mat.matrix);
    gsl_matrix_set (m, 0, 0, 0.0);
    gsl_matrix_set (m, 0, 1, 1.0);
    gsl_matrix_set (m, 1, 0, -(g/l));
    gsl_matrix_set (m, 1, 1, 0.0);
    dfdt[0] = 0.0;
    dfdt[1] = 0.0;
    return GSL_SUCCESS;
}

int main (void) {
    Params params;
    params.l = 1.0;
    params.g = 9.81;
    gsl_odeiv2_system sys = {func, jac, 2, &params};

    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (
            &sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0
    );
    int i;
    double t = 0.0, t1 = 100.0;
    double y[2] = { 0.1, 0.0 };

    for (i = 1; i <= 1000; i++) {
        double ti = i*t1/1000.0;
        int status = gsl_odeiv2_driver_apply (d, &t, ti, y);

        if (status != GSL_SUCCESS) {
            fprintf(stderr, "error, return value=%d\n", status);
            break;
        }

        printf ("%.5le %.5le %.5le\n", t, y[0], y[1]);
    }

    gsl_odeiv2_driver_free (d);
    return EXIT_SUCCESS;
}

ext: c

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material


text:
/* GPC: A library for the solution of General Point Correspondence problems.
  Copyright (C) 2006 Andrea Censi (andrea at censi dot org)

  This program 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 2
  of the License, or (at your option) any later version.

  This program 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, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

#include <math.h>
#include <gsl/gsl_math.h>
#include "gpc.h"
#include "gpc_utils.h"

/* Note that we use static values here so we don't need to evaluate that all the time */
#define M(matrix, rows, col) static gsl_matrix*matrix=0; if(!matrix) matrix = gsl_matrix_alloc(rows,col);
#define MF(matrix) /*gsl_matrix_free(matrix)*/

int gpc_solve(int K, const struct gpc_corr*c,
	const double*x0, const double *cov_x0,
	double *x_out)
{
	M(bigM,    4,4); M(g,  4,1); M(bigM_k,2,4);
	M(bigM_k_t,4,2); M(C_k,2,2); M(q_k,   2,1);
	M(temp41,  4,1); M(temp22,2,2);	M(temp22b,2,2);
	M(temp42,  4,2); M(temp44,4,4);	M(temp21, 2,1);
	M(temp22c, 2,2); M(temp12,1,2);

	gsl_matrix_set_zero(bigM);
	gsl_matrix_set_zero(g);
	gsl_matrix_set_zero(temp42);

	double d_bigM[4][4] = { {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}};
	double d_g[4] = {0, 0, 0, 0};
	int k;
	for(k=0;k<K;k++) {
		if(!c[k].valid) continue;

		double C00 = c[k].C[0][0];
		double C01 = c[k].C[0][1];
		double C10 = c[k].C[1][0];
		double C11 = c[k].C[1][1];

		if(C01 != C10) {
			fprintf(stderr, "k=%d; I expect C to be a symmetric matrix.\n", k);
			return 0;
		}
#if GPC_CHECK_SEMIDEF
		double det = C00*C11 - C01*C10;
		double trace = C00 + C11;
		int is_semidef_pos = (det >= 0) && (trace>0);
		if(!is_semidef_pos) {
			fprintf(stderr, "k=%d; I expect the matrix to be semidef positive (det>=0 and trace>0), det = %.20f trace= %.10f C = [%.15f,%.15f;%.15f %.15f]\n",
				k, det, trace,	C00,C01,C10,C11);
			return 0;
		}
#endif

		double qx = c[k].q[0];
		double qy = c[k].q[1];
		double px = c[k].p[0];
		double py = c[k].p[1];

		/* [ C00,  c01,  px C00 + c01 py , c01 px - py C00 ] */
		d_bigM[0][0] += C00;
		d_bigM[0][1] += C01;
		d_bigM[0][2] += +px * C00 + py * C01;
		d_bigM[0][3] += -py * C00 + px * C01;

		/*  [ C10 ,  C11 , py C11 + px C10 , px C11 - py C10 ] */
		d_bigM[1][0] += C10;
		d_bigM[1][1] += C11;
		d_bigM[1][2] += +px * C10 + py * C11;
		d_bigM[1][3] += +px * C11 - py * C10;

		/*Col 1 = [ py C10 + px C00 ]
		 Col 2 = [ py C11 + c01 px ]
		 Col 3 = [ py (py C11 + px C10) + px (px C00 + c01 py) ]
		 Col 4 = [ py (px C11 - py C10) + px (c01 px - py C00) ]
		*/
		d_bigM[2][0] += px * C00 + py * C10;
		d_bigM[2][1] += px * C01 + py * C11;
		d_bigM[2][2] += (px*px)*(+C00) + (px*py)*(+C10+C01) + (py*py)*(+C11);
		d_bigM[2][3] += (px*px)*(+C01) + (px*py)*(-C00+C11) + (py*py)*(-C10);

		/*Col 1 = [ px C10 - py C00 ]
		  Col 2 = [ px C11 - c01 py ]
		 Col 3 = [ px (py C11 + px C10) - py (px C00 + c01 py) ]
		 Col 4 = [ px (px C11 - py C10) - py (c01 px - py C00) ]*/
		d_bigM[3][0] += -py * C00 + px * C10;
		d_bigM[3][1] += -py * C01 + px * C11;
		d_bigM[3][2] += (px*px)*(+C10) + (px*py)*(-C00+C11) + (py*py)*(-C01);
		d_bigM[3][3] += (px*px)*(+C11) + (px*py)*(-C10-C01) + (py*py)*(+C00);

		d_g[0] += C00*qx+C10*qy;
		d_g[1] += C01*qx+C11*qy;
		d_g[2] += qx * (C00*px+C01*py) + qy * (C10*px+C11*py);
		d_g[3] += qx * (C00*(-py)+C01*px) + qy * (C10*(-py)+C11*px);

	}

	{
	unsigned int a,b;
	 for(a=0;a<4;a++)
			*gsl_matrix_ptr(g, a, 0) = -2 * d_g[a];
	 for(a=0;a<4;a++)
		for(b=0;b<4;b++)
		gsl_matrix_set(bigM, a, b,  2 * d_bigM[a][b]);
	}

	M(mA,2,2); gms(mA,0,0, gmg(bigM,0,0)); gms(mA,0,1, gmg(bigM,0,1));
	           gms(mA,1,0, gmg(bigM,1,0)); gms(mA,1,1, gmg(bigM,1,1));
	M(mB,2,2); gms(mB,0,0, gmg(bigM,0,2)); gms(mB,0,1, gmg(bigM,0,3));
	           gms(mB,1,0, gmg(bigM,1,2)); gms(mB,1,1, gmg(bigM,1,3));
	M(mD,2,2); gms(mD,0,0, gmg(bigM,2,2)); gms(mD,0,1, gmg(bigM,2,3));
	           gms(mD,1,0, gmg(bigM,3,2)); gms(mD,1,1, gmg(bigM,3,3));

	M(mS,2,2); M(mSa,2,2);

	/*	mS = mD - mB.trans * mA.inv * mB;
	  temp22b = inv(A) */
	m_inv(mA, temp22b);
	/* temp22c = inv(A) * mB           */
	m_mult(temp22b, mB, temp22c);
	/* temp22 = mB'               */
	m_trans(mB, temp22);
	m_mult(temp22, temp22c, temp22b);
	m_scale(-1.0,temp22b);
	m_add(mD,temp22b,mS);

	/* mSa = mS.inv * mS.det; */
	m_inv(mS, mSa);
	m_scale(m_det(mS), mSa);

	if(TRACE_ALGO) {
		m_display("mA",mA);
		m_display("mB",mB);
		m_display("mD",mD);
		m_display("mS",mS);
		m_display("mSa",mSa);
	}

	M(g1,2,1);	M(g2,2,1);
	M(g1t,1,2);	M(g2t,1,2);
	M(mAi,2,2);	M(mBt,2,2);

	gms(g1,0,0,gmg(g,0,0));
	gms(g1,1,0,gmg(g,1,0));
	gms(g2,0,0,gmg(g,2,0));
	gms(g2,1,0,gmg(g,3,0));
	m_trans(g1, g1t);
	m_trans(g2, g2t);
	m_trans(mB, mBt);
	m_inv(mA, mAi);

	M(m1t,1,2);	M(m1,2,1);
	M(m2t,1,2);	M(m2,2,1);
	M(m3t,1,2);	M(m3,2,1);

	/* m1t = g1t*mAi*mB */
	m_mult(g1t,mAi,temp12);
	m_mult(temp12,mB,m1t);

	m_trans(m1t,m1);
	/*     m2t = m1t*mSa    */
	m_mult(m1t,mSa,m2t);
	m_trans(m2t,m2);
	/* m3t = g2t*mSa     */
	m_mult(g2t,mSa,m3t);
	m_trans(m3t,m3);

	double p[3] = {
		  m_dot(m2t,m2) - 2*m_dot(m2t,m3) +   m_dot(m3t,m3),
		4*m_dot(m2t,m1) - 8*m_dot(m2t,g2) + 4*m_dot(g2t,m3),
		4*m_dot(m1t,m1) - 8*m_dot(m1t,g2) + 4*m_dot(g2t,g2)};

	double l[3] = {m_det(mS), 2*gmg(mS,0,0)+2*gmg(mS,1,1), 4};

	/* q = p - l^2        */
	double q[5] = {p[0]-(l[0]*l[0]), p[1]-(2*l[1]*l[0]),
		p[2]-(l[1]*l[1]+2*l[0]*l[2]), -(2*l[2]*l[1]), -(l[2]*l[2])};

	if(TRACE_ALGO) {
		fprintf(stderr, "p = %f %f %f \n", p[2], p[1], p[0]);
		fprintf(stderr, "l = %f %f %f \n", l[2], l[1], l[0]);
		fprintf(stderr, "q = %f %f %f %f %f \n", q[4],  q[3],  q[2], q[1], q[0]);
	}

	/*
	double lambdas[4];
	if(!poly_real_roots(5, q, lambdas)) {
		fprintf(stderr, "Cannot solve polynomial.\n");
		return 0;
	}

	double lambdas_error[4];
	double lambdas_pose[4][3];

	for(int i=0;i<4;i++) {
		double lambda = lambdas[i];

		if(TRACE_ALGO) {
			fprintf(stderr, "lambda = %f \n", lambda);
		}

		M(W,4,4); gsl_matrix_set_zero(W); gms(W,2,2,1.0); gms(W,3,3,1.0);
		M(x,4,1);

		m_scale(2*lambda, W);
		gsl_matrix_add(bigM,W);
		m_inv(bigM, temp44);
		m_mult(temp44, g, x);
		m_scale(-1.0, x);

		lambdas_pose[i][0] = gmg(x,0,0);
		lambdas_pose[i][1] = gmg(x,1,0);
		lambdas_pose[i][2] = atan2(gmg(x,3,0),gmg(x,2,0));

		lambdas_error[i] = gpc_total_error(c, K, lambdas_pose[i]);
	}

	int lowest_error = 0;
	for(int i=0;i<4;i++) {
		printf("#%d lambda=%lf error=%lf\n",i,lambdas[i],lambdas_error[i]);
		if(lambdas_error[i] < lambdas_error[lowest_error])
			lowest_error = i;
	}

	double lr;
	poly_greatest_real_root(5,q,&lr);
	printf("Choose %d: lambda = %lf   bigger real root = %lf \n",lowest_error,lambdas[lowest_error],lr);

	x_out[0]=lambdas_pose[lowest_error][0];
	x_out[1]=lambdas_pose[lowest_error][1];
	x_out[2]=lambdas_pose[lowest_error][2];
	*/

	double lambda;
	if(!poly_greatest_real_root(5, q, &lambda)) return 0;

	M(W,4,4); gsl_matrix_set_zero(W); gms(W,2,2,1.0); gms(W,3,3,1.0);
	M(x,4,1);

	m_scale(2*lambda, W);
	gsl_matrix_add(bigM,W);
	m_inv(bigM, temp44);
	m_mult(temp44, g, x);
	m_scale(-1.0, x);

	x_out[0] = gmg(x,0,0);
	x_out[1] = gmg(x,1,0);
	x_out[2] = atan2(gmg(x,3,0),gmg(x,2,0));

	if(TRACE_ALGO) {
		fprintf(stderr, "x =  %f  %f %f deg\n", x_out[0], x_out[1],x_out[2]*180/M_PI);
	}

	MF(mA); MF(mB); MF(mD); MF(mS); MF(mSa);
	MF(m1t);	MF(m1);	MF(m2t);	MF(m2);
	MF(m3t);	MF(m3);	MF(W); 	MF(x);
	MF(bigM); MF(g); MF(bigM_k);
	MF(bigM_k_t); MF(C_k); MF(q_k);
	MF(temp42); MF(temp44);	MF(temp21);
	MF(temp41); MF(temp22);	MF(temp22b);
	MF(temp22c); MF(temp12);
	MF(g1);	MF(g2);
	MF(g1t);	MF(g2t);
	MF(mAi);	MF(mBt);
	return 1;
}

double gpc_error(const struct gpc_corr*co, const double*x) {
	double c = cos(x[2]);
	double s = sin(x[2]);
	double e[2];
	e[0] = c*(co->p[0]) -s*(co->p[1]) + x[0] - co->q[0];
	e[1] = s*(co->p[0]) +c*(co->p[1]) + x[1] - co->q[1];
	double this_error = e[0]*e[0]*co->C[0][0]+2*e[0]*e[1]*co->C[0][1]+e[1]*e[1]*co->C[1][1];

	if(0) /* due to limited numerical precision, error might be negative */
	if(this_error < 0) {
		fprintf(stderr, "Something fishy: error = %lf e = [%lf %lf]  C = [%lf,%lf;%lf,%lf]\n", this_error,
			e[0],e[1],co->C[0][0],co->C[0][1],co->C[1][0],co->C[1][1]);
	}
	return this_error;
}

double gpc_total_error(const struct gpc_corr*co, int n, const double*x){
	double error = 0;
	for(int i=0;i<n;i++) {
		if(!co[i].valid) continue;
		error += gpc_error(co+i,x);
	}
	if(0) /* due to limited numerical precision, error might be negative */
	if(error<0) {
		fprintf(stderr, "Something fishy!\n");
	}
	return error;
}

/*
               [       C00       ]         [       c01       ]
                [                 ]         [                 ]
                [       C10       ]         [       C11       ]
(%o10)  Col 1 = [                 ] Col 2 = [                 ]
                [ py C10 + px C00 ]         [ py C11 + c01 px ]
                [                 ]         [                 ]
                [ px C10 - py C00 ]         [ px C11 - c01 py ]
         [              px C00 + c01 py              ]
         [                                           ]
         [              py C11 + px C10              ]
         [                                           ]
 Col 3 = [   2                     2                 ]
         [ py  C11 + px py C10 + px  C00 + c01 px py ]
         [                                           ]
         [               2                         2 ]
         [ px py C11 + px  C10 - px py C00 - c01 py  ]
         [              c01 px - py C00              ]
         [                                           ]
         [              px C11 - py C10              ]
         [                                           ]
 Col 4 = [               2                         2 ]
         [ px py C11 - py  C10 - px py C00 + c01 px  ]
         [                                           ]
         [   2                     2                 ]
         [ px  C11 - px py C10 + py  C00 - c01 px py ]*/


text:
/* randist/pascal.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
 *
 * This program 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.
 *
 * This program 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* The Pascal distribution is a negative binomial with valued integer n

   prob(k) =  (n - 1 + k)!/(n!(k - 1)!) *  p^n (1-p)^k for k = 0, 1, ..., n

   */

unsigned int
gsl_ran_pascal (const gsl_rng * r, double p, unsigned int n)
{
  /* This is a separate interface for the pascal distribution so that
     it can be optimized differently from the negative binomial in
     future.

     e.g. if n < 10 it might be faster to generate the Pascal
     distributions as the sum of geometric variates directly.  */

  unsigned int k = gsl_ran_negative_binomial (r, p, (double) n);
  return k;
}

double
gsl_ran_pascal_pdf (const unsigned int k, const double p, unsigned int n)
{
  double P = gsl_ran_negative_binomial_pdf (k, p, (double) n);
  return P;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
    #ifndef LVVMODEL_H
    #define LVVMODEL_H

    #include <iostream>
    #include <iomanip>
    using std::fixed;
    using std::showpos;

    #include <string>
    #include <gsl/gsl_multifit.h>
    #include <gsl/gsl_poly.h>
    #include "lvv.h"
    #include "math.h"
    using lvv::powi;

    namespace lvv {

    int  const  poly_order = 4;

//	template<int M>
class   Model   {     public: ///////////////////////////////////////////////////////////// Model Data - μ² polinom fitting

Model () : n(0), m(0), x(NULL), y(NULL), X(NULL),  cov(NULL),  c(NULL)              {}

Model (double _x[], double _y[], int _n, int _m=-1) : n(_n), m(_m), x(_x), y(_y)    { ////////////////////////////////   CTOR

        // N - data points [0..N) ( 5 for NF )
 	// M - number of model parameters [0..M)
 	// if polynomial used in model:
 	//      M-1 -  poly order  C[0]..C[m-1]

	assert(_n >= 3 && "not implemented");
	if (_m == -1) {      //  then need autodetect
	    if      (n < 6 ) m = n-2;
	    else if (n < 10) m = n-3;
	    else if (n < 20) m = n-4;
	}

	assert(n >= m);

	X   = gsl_matrix_alloc (n, m) ;		assert(X);
	c   = gsl_vector_alloc (m   ) ;		assert(c);
	cov = gsl_matrix_alloc (m, m) ;		assert(cov);

	xv = gsl_vector_view_array (_x, n);
	yv = gsl_vector_view_array (_y, n);

	/////    INIT - PREP MATRIX for POLIFIT
	for (i = 0; i<n; i++)    for (int p=0; p<m; ++p)     gsl_matrix_set (X, i, p, powi(_x[i], p));

	/////    FIT
	gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (n, m);
	//  gsl_multifit_linear     (X, y,                         c, cov, &chisq, work); // int gsl_multifit_linear     (const gsl_matrix *X, const gsl_vector *y,                           gsl_vector *c, gsl_matrix *cov, double *chisq, gsl_multifit_linear_workspace *ws)
	    gsl_multifit_linear_svd (X, &yv.vector, svd_tol, &rank, c, cov, &chisq, work); // int gsl_multifit_linear_svd (const gsl_matrix *X, const gsl_vector *y, double tol, size t* rank, gsl_vector *c, gsl_matrix *cov, double *chisq, gsl_multifit_linear_workspace *ws)
	gsl_multifit_linear_free (work);
 };

~Model() { ////////////////////////////////////////////////////////////////// DTOR
	if(cov)    gsl_matrix_free (cov);
	if(c  )    gsl_vector_free (c);
	if(X  )    gsl_matrix_free (X);
 };

double  estimate (const double x)   const {
	double y_est = gsl_vector_get(c,0) ;
	double pow_x = 1.0                 ; // x[0,0]^0

	for (int p=1;  p<m;  ++p) {
	    pow_x *= x;
	    y_est += gsl_vector_get(c,p)*pow_x;
	}

	return y_est;
 };

double  inverse_estimate (const double yy)    const {

	double xx = -1.0;

	double *C=c->data;

	if        (m>2) {

		double c0_save = C[0];
		C[0] =  c0_save-yy;

		struct { double real; double imag; } z[poly_order-1];

		int order = poly_order;
		//if ( abs(C[poly_order-1]) < 0.000001 )	// reduce polinomial order of last elem == 0
		if ( C[poly_order-1] == 0 )	// reduce polinomial order of last elem == 0
			order = poly_order-1;
		    			assert (C[poly_order-1] != 0);	// hope we don't need to do it twice

		//  solve
		gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (order);
		gsl_poly_complex_solve (gsl_vector_ptr(c,0), order, w, (double *)z);
		gsl_poly_complex_workspace_free (w);

		//  select real root from list of complex roots
		double oorr_distance         = numeric_limits<double>::max(); // out_of_range root distance to normal interval
		bool   got_root              = false;
		bool   got_out_of_range_root = false;
		double out_of_range_root     = -100;
		double real_root             = -100000;

		//double x_step                = x[1]-x[0];                    // with step p-rmse go from 0.9904215  to 0.9907456
		double x_lower_bound         = x[0]    /* -x_step/2.0 */ ;
		double x_upper_bound         = x[n-1]  /* +x_step/2.0 */ ;

		for  (int i = 0;  i < order-1;  i++)    {				// for all roots
												//cerr << " \t("<<i<<") "<<  z[2*i] << "+i" << z[2*i+1];
			if  (abs(z[i].imag)  >  0.00000001)  continue;

			if   ( x_lower_bound <=  z[i].real   &&  z[i].real <=  x_upper_bound)   {  // in x valid range
				real_root = z[i].real;
				got_root  = true;
															//cerr << " rr=" << real_root;
				break;	// found it
			}

			double below_lower_distance =  x_lower_bound - z[i].real;
			if  (  below_lower_distance > 0 )  {
				if  (( got_out_of_range_root  &&  below_lower_distance < oorr_distance )  ||  ! got_out_of_range_root)  {
					oorr_distance         =   below_lower_distance; 	// new shortest oorr distance
					got_out_of_range_root = true;
					out_of_range_root     = x_lower_bound;
													//cerr  << " oor=" << out_of_range_root ;
					continue;
				}
			}

			double above_upper_distance =  z[i].real - x_upper_bound ;
			if   ( above_upper_distance > 0 )  {
				if (( got_out_of_range_root  && above_upper_distance  < oorr_distance )  ||  ! got_out_of_range_root)  {
					oorr_distance         = above_upper_distance;
					got_out_of_range_root = true;
					out_of_range_root     = x_upper_bound;
													//cerr  << " oor=" << out_of_range_root ;
					continue;
				}
			}

		};

	    if	     (got_root)				xx = real_root;
	    else  if (got_out_of_range_root)		xx = out_of_range_root;
	    else					{ xx=3.66; cerr << "#model::inverse_estimate(): root not found\n";  ; }
														//cerr << "  ROOT=" << xx << endl;

	    gsl_vector_set(c, 0, c0_save);

	    /*
	    double c0 =  gsl_vector_get(c,0);
	    double c1 =  gsl_vector_get(c,1);
	    double c2 =  gsl_vector_get(c,2);
	    double c3 =  gsl_vector_get(c,3);
	    if ( 0 > (27*pow2(c3)*pow2(yy)+(-54*c0*pow2(c3)+18*c1*c2*c3-4*pow3(c2))*yy+27*pow2(c0)*pow2(c3)+(4*pow3(c1)-18*c0*c1*c2)*c3+4*c0*pow3(c2)-pow2(c1)*pow2(c2)))  {
		cout << "inverse_estimate-4-1: cannot solve: c0-2" <<  c0 << " " << c1 << " " <<  c2 << " " << c3 << " y: " << yy << endl;
		return 3.6;
	    }
	    if ( 0 > 27*pow2(c3)*pow2(yy)+(-54*c0*pow2(c3)+18*c1*c2*c3-4*pow3(c2))*yy+27*pow2(c0)*pow2(c3)+(4*pow3(c1)-18*c0*c1*c2)*c3+4*c0*pow3(c2)-pow2(c1)*pow2(c2))  {
		cout << "inverse_estimate-4-2: cannot solve: c0-2" <<  c0 << " " << c1 << " " <<  c2 << " " << c3 << " y: " << yy << endl;
		return 3.6;
	    }
	    xx = cbrt((sqrt(27*pow2(c3)*pow2(yy)+(-54*c0*pow2(c3)+18*c1*c2*c3-4*pow3(c2))*yy+27*pow2(c0)*pow2(c3)+(4*pow3(c1)-18*c0*c1*c2)*c3+4*c0*pow3(c2)-pow2(c1)*pow2(c2))/(6*sqrt(3)*pow2(c3))-(-27*pow2(c3)*yy+27*c0*pow2(c3)-9*c1*c2*c3+2*pow3(c2))/(54*pow3(c3))))-cbrt((3*c1*c3-pow2(c2))/(9*pow2(c3)*(sqrt(27*pow2(c3)*pow2(yy)+(-54*c0*pow2(c3)+18*c1*c2*c3-4*pow3(c2))*yy+27*pow2(c0)*pow2(c3)+(4*pow3(c1)-18*c0*c1*c2)*c3+4*c0*pow3(c2)-pow2(c1)*pow2(c2))/(6*sqrt(3)*pow2(c3))-(-27*pow2(c3)*yy+27*c0*pow2(c3)-9*c1*c2*c3+2*pow3(c2))/(54*pow3(c3)))))-c2/(3*c3);

	} else if (m==3){
	    double c0 =  gsl_vector_get(c,0);
	    double c1 =  gsl_vector_get(c,1);
	    double c2 =  gsl_vector_get(c,2);

	    if(4*c2*yy-4*c0*c2+pow2(c1) < 0) {
		cout << "inverse_estimate-3: cannot solve: c0-2" <<  c0 << " " << c1 << " " <<  c2 << " y: " << yy << endl;
		return 3.6;
	    }

	    xx = (sqrt(4*c2*yy-4*c0*c2+pow2(c1))-c1)/(2*c2);
	    //xx = -(sqrt(4*c2*yy-4*c0*c2+pow2(c1))+c1)/(2*c2);

	   */
	} else if (m==2){
	    double c0 =  gsl_vector_get(c,0);
	    double c1 =  gsl_vector_get(c,1);
	    xx = (yy-c0)/c1;
	} else {
	    assert(false);
	}

	return xx;
 };

void	print (string title="")  const {
	cout << "# :gnuplot: set key bottom; set ytics 0,0.5; set yrange [-0.1:1.1]; set xrange [0:6]; set grid; set pointsize 1.5;"
		<< (title.length()>0  ?  "set title \" " +  title + "\";"  :  "") <<

		    "plot "
			" \"pipe\" using 1:3 title \"model\" with lines 5, "
			" \"pipe\" using 4:2 title \"inv x(y)\" with point 4,"
			" \"pipe\" using 1:2 title \"data\" with points 11 "
		";\n";

	format		table(" % 9f    %20t % 9f   %40t % 9f   %60t % 9f\n");
	string		_("\" \"");
							cout << table	% "# X" % "  Y"		% "Y-mod"     	% "X-mod-inv" ;

	for (int i=0;  i<n; ++i                )	cout << table	% x[i] 	% y[i] 		% _    		% _ ;
	for (double xx=0.0;  xx<=6.0; xx+=1./3.)	cout << table	% xx   	% _		% estimate(xx)	% _  ;
	for (double yy=0.0;  yy<=1.0; yy+=0.05 )	cout << table	% _ 	% yy   		% _    		% inverse_estimate(yy) ;
 };

private: // member vars
        int                  n         , m;
        const static double  svd_tol = 0.000001;
        size_t               rank;
        int                  i;
        double               xi        , yi, ei, chisq;
        gsl_vector_view      xv        , yv;
        double    	   * x;
        double             * y;
        gsl_matrix          *X         , *cov;

    public:
        gsl_vector          *c;
 };

    }      // namespace lvv
    #endif // LVVMODEL_H

max_stars_repo_name: lvv/lvvlib

max_issues_repo_name: lvv/lvvlib

max_forks_repo_name: lvv/lvvlib


text:
/*
 This file originates from LibSVM <http://www.csie.ntu.edu.tw/~cjlin/libsvm/>

 Copyright (c) 2000-2013 Chih-Chung Chang and Chih-Jen Lin
 All rights reserved.

 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions
 are met:

 1. Redistributions of source code must retain the above copyright
 notice, this list of conditions and the following disclaimer.

 2. Redistributions in binary form must reproduce the above copyright
 notice, this list of conditions and the following disclaimer in the
 documentation and/or other materials provided with the distribution.

 3. Neither name of copyright holders nor the names of its contributors
 may be used to endorse or promote products derived from this software
 without specific prior written permission.

 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR
 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

 */

/*
 This file is part of the Princeton FCMA Toolbox-- namely
 the "svm_cross_validation_no_shuffle" function, an edited
 version of libsvm's "svm_cross_validation"
*/

#ifndef _SVM_H
#define _SVM_H

#define LIBSVM_VERSION 310
#define MAXSVMITERATION 20000000

#if defined(__GNUC__)
// gcc supports this syntax
#define ALIGNED(x) __attribute__((aligned(x)))
#else
// visual c++, clang, icc
#define ALIGNED(x) __declspec(align(x))
#endif

//#include "LibSVM.h"
#include <cassert>

// BLAS dependency
#ifdef USE_MKL
#include <mkl.h>
#else
typedef int MKL_INT;
#if defined __APPLE__
#include <Accelerate/Accelerate.h>
#else
extern "C" {
#include <cblas.h>
}
#endif
#endif

// SSE/AVX instrinsics
#if defined(_MSC_VER)
// Microsoft C/C++-compatible compiler */
#include <intrin.h>
#elif defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
// GCC-compatible compiler (gcc,clang,icc) targeting x86/x86-64
#include <x86intrin.h>
#endif

#ifdef __cplusplus
extern "C" {
#endif

// unsigned long long counter;
typedef float Qfloat;
typedef signed char schar;

extern int libsvm_version;

// use SOA instead of AOS --later!!
struct svm_node {
  int index;
  float value;
};

// change y from double to schar, x is one-step pointer now
struct svm_problem {
  int l;
  schar *y;
  struct svm_node **x;
};

struct svm_parameter {
  int svm_type;
  int kernel_type;
  int degree;   /* for poly */
  double gamma; /* for poly/rbf/sigmoid */
  double coef0; /* for poly/sigmoid */

  /* these are for training only */
  double cache_size; /* in MB */
  double eps;        /* stopping criteria */
  double C;          /* for C_SVC, EPSILON_SVR and NU_SVR */
  int nr_weight;     /* for C_SVC */
  int *weight_label; /* for C_SVC */
  double *weight;    /* for C_SVC */
  double nu;         /* for NU_SVC, ONE_CLASS, and NU_SVR */
  double p;          /* for EPSILON_SVR */
  int shrinking;     /* use the shrinking heuristics */
  int probability;   /* do probability estimates */
};

//
// svm_model
//
struct svm_model {
  struct svm_parameter param; /* parameter */
  int nr_class;         /* number of classes, = 2 in regression/one class svm */
  int l;                /* total #SV */
  struct svm_node **SV; /* SVs (SV[l]) */
  double **sv_coef;     /* coefficients for SVs in decision functions
                       (sv_coef[k-1][l]) */
  double *rho; /* constants in decision functions (rho[k*(k-1)/2]) */

  /* for classification only */

  int *label; /* label of each class (label[k]) */
  int *nSV;   /* number of SVs for each class (nSV[k]) */
              /* nSV[0] + nSV[1] + ... + nSV[k-1] = l */
  /* XXX */
  int free_sv; /* 1 if svm_model is created by svm_load_model*/
               /* 0 if svm_model is created by svm_train */
};

enum {
  C_SVC
};
    /* svm_type */  // only C_SVC is needed for now
enum {
  LINEAR,
  PRECOMPUTED
}; /* kernel_type */

typedef struct svm_problem SVMProblem;
typedef struct svm_parameter SVMParameter;
typedef struct svm_node SVMNode;

struct svm_model *svm_train(const struct svm_problem *prob,
                            const struct svm_parameter *param);
void svm_cross_validation(const struct svm_problem *prob,
                          const struct svm_parameter *param,
                          int nr_fold, double *target);

void svm_cross_validation_no_shuffle(const struct svm_problem *prob,
                                     const struct svm_parameter *param,
                                     int nr_fold, double *target);

int svm_get_svm_type(const struct svm_model *model);
int svm_get_nr_class(const struct svm_model *model);
void svm_get_labels(const struct svm_model *model, int *label);

double svm_predict_values(const struct svm_model *model,
                          const struct svm_node *x, double *dec_values);
double svm_predict(const struct svm_model *model, const struct svm_node *x);
double svm_predict_distance(const struct svm_model *model,
                            const struct svm_node *x);

void svm_free_model_content(struct svm_model *model_ptr);
void svm_free_and_destroy_model(struct svm_model **model_ptr_ptr);
void svm_destroy_param(struct svm_parameter *param);

const char *svm_check_parameter(const struct svm_problem *prob,
                                const struct svm_parameter *param);
int svm_check_probability_model(const struct svm_model *model);

void svm_set_print_string_function(void (*print_func)(const char *));

int svm_save_model(const char *model_file_name, const struct svm_model *model);
struct svm_model *svm_load_model(const char *model_file_name);

#ifdef __cplusplus
}
#endif

#endif /* _SVM_H */

max_stars_repo_name: PrincetonUniversity/fcma-toolbox

max_issues_repo_name: PrincetonUniversity/fcma-toolbox

max_forks_repo_name: PrincetonUniversity/fcma-toolbox


text:
/* rng/rng.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 *
 * This program 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 2 of the License, or (at
 * your option) any later version.
 *
 * This program 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>

gsl_rng *
gsl_rng_alloc (const gsl_rng_type * T)
{

  gsl_rng *r = (gsl_rng *) malloc (sizeof (gsl_rng));

  if (r == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for rng struct",
                        GSL_ENOMEM, 0);
    };

  r->state = malloc (T->size);

  if (r->state == 0)
    {
      free (r);         /* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for rng state",
                        GSL_ENOMEM, 0);
    };

  r->type = T;

  gsl_rng_set (r, gsl_rng_default_seed);        /* seed the generator */

  return r;
}

int
gsl_rng_memcpy (gsl_rng * dest, const gsl_rng * src)
{
  if (dest->type != src->type)
    {
      GSL_ERROR ("generators must be of the same type", GSL_EINVAL);
    }

  memcpy (dest->state, src->state, src->type->size);

  return GSL_SUCCESS;
}

gsl_rng *
gsl_rng_clone (const gsl_rng * q)
{
  gsl_rng *r = (gsl_rng *) malloc (sizeof (gsl_rng));

  if (r == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for rng struct",
                        GSL_ENOMEM, 0);
    };

  r->state = malloc (q->type->size);

  if (r->state == 0)
    {
      free (r);         /* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for rng state",
                        GSL_ENOMEM, 0);
    };

  r->type = q->type;

  memcpy (r->state, q->state, q->type->size);

  return r;
}

void
gsl_rng_set (const gsl_rng * r, unsigned long int seed)
{
  (r->type->set) (r->state, seed);
}

#ifndef HIDE_INLINE_STATIC
unsigned long int
gsl_rng_get (const gsl_rng * r)
{
  return (r->type->get) (r->state);
}

double
gsl_rng_uniform (const gsl_rng * r)
{
  return (r->type->get_double) (r->state);
}

double
gsl_rng_uniform_pos (const gsl_rng * r)
{
  double x ;
  do
    {
      x = (r->type->get_double) (r->state) ;
    }
  while (x == 0) ;

  return x ;
}

/* Note: to avoid integer overflow in (range+1) we work with scale =
   range/n = (max-min)/n rather than scale=(max-min+1)/n, this reduces
   efficiency slightly but avoids having to check for the out of range
   value.  Note that range is typically O(2^32) so the addition of 1
   is negligible in most usage. */

unsigned long int
gsl_rng_uniform_int (const gsl_rng * r, unsigned long int n)
{
  unsigned long int offset = r->type->min;
  unsigned long int range = r->type->max - offset;
  unsigned long int scale;
  unsigned long int k;

  if (n > range || n == 0)
    {
      GSL_ERROR_VAL ("invalid n, either 0 or exceeds maximum value of generator",
                     GSL_EINVAL, 0) ;
    }

  scale = range / n;

  do
    {
      k = (((r->type->get) (r->state)) - offset) / scale;
    }
  while (k >= n);

  return k;
}
#endif

unsigned long int
gsl_rng_max (const gsl_rng * r)
{
  return r->type->max;
}

unsigned long int
gsl_rng_min (const gsl_rng * r)
{
  return r->type->min;
}

const char *
gsl_rng_name (const gsl_rng * r)
{
  return r->type->name;
}

size_t
gsl_rng_size (const gsl_rng * r)
{
  return r->type->size;
}

void *
gsl_rng_state (const gsl_rng * r)
{
  return r->state;
}

void
gsl_rng_print_state (const gsl_rng * r)
{
  size_t i;
  unsigned char *p = (unsigned char *) (r->state);
  const size_t n = r->type->size;

  for (i = 0; i < n; i++)
    {
      /* FIXME: we're assuming that a char is 8 bits */
      printf ("%.2x", *(p + i));
    }

}

void
gsl_rng_free (gsl_rng * r)
{
  free (r->state);
  free (r);
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include "parameters.h" // Includes all defined parameters
#include "axi.h" // Axisymmetric coordinates
#include "navier-stokes/centered.h" // To solve the Navier-Stokes
#include "two-phase.h" // Implements two-phase flow
// #include <omp.h> // For openMP parallel
#include <gsl/gsl_fit.h>

double * forces_array;
double * times_array;

double gradient = 0.5;
double intercept = 4.9;

int main () {

    forces_array = malloc(INTERP_NO * sizeof(double));
    times_array = malloc(INTERP_NO * sizeof(double));

    DT = 1e-2;

    run();

    free(forces_array);
    free(times_array);
}

event force(i++) {
    if (i < INTERP_NO) {
        forces_array[i] = gradient * t + intercept;
        times_array[i] = t;
    } else {
        double c0, c1, cov00, cov01, cov11, sumsq;
        gsl_fit_linear ( times_array, 1, forces_array, 1, INTERP_NO, &c0, &c1, \
            &cov00, &cov01, &cov11, &sumsq);

        fprintf(stderr, "t = %g, c0 = %g, c1 = %g\n", t, c0, c1);

        double current_force = gradient * t + intercept;
        double interp_force = c0 + c1 * t;
        fprintf(stderr, "f = %g, f_interp = %g", current_force, interp_force);
        fprintf(stderr, "\n");

        // #pragma omp critical
        for (int j = 0; j < INTERP_NO - 1; j++) {
            forces_array[j] = forces_array[j + 1];
            times_array[j] = times_array[j + 1];
        }
        forces_array[INTERP_NO - 1] = current_force;
        times_array[INTERP_NO - 1] = t;
    }
}

event end(t = 1) {
    fprintf(stderr, "Ended\n");
}

max_stars_repo_name: MNegus/plate-impact

max_issues_repo_name: MNegus/plate-impact

max_forks_repo_name: MNegus/plate-impact


text:
/* filter/median.c
 *
 * Contains routines related to the standard median filter
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * This program 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.
 *
 * This program 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_filter.h>
#include <gsl/gsl_errno.h>

gsl_filter_median_workspace *
gsl_filter_median_alloc(const size_t K)
{
  gsl_filter_median_workspace *w;
  size_t H = K / 2;

  w = calloc(1, sizeof(gsl_filter_median_workspace));
  if (w == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
    }

  w->movstat_workspace_p = gsl_movstat_alloc(2*H + 1);
  if (w->movstat_workspace_p == NULL)
    {
      gsl_filter_median_free(w);
      GSL_ERROR_NULL ("failed to allocate space for movstat workspace", GSL_ENOMEM);
    }

  return w;
}

void
gsl_filter_median_free(gsl_filter_median_workspace * w)
{
  if (w->movstat_workspace_p)
    gsl_movstat_free(w->movstat_workspace_p);

  free(w);
}

/*
gsl_filter_median()
  Standard median filter

Inputs: endtype - end point handling
        x       - input vector
        y       - output vector
        w       - workspace
*/

int
gsl_filter_median(const gsl_filter_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_filter_median_workspace * w)
{
  int status = gsl_movstat_median(endtype, x, y, w->movstat_workspace_p);
  return status;
}

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include "example_z.h"
#include <stdlib.h>
#include <lapacke.h>

void random_zvec(int len, BBLAS_Complex64_t *vec)
{
//	BBLAS_Complex64_t *work =
//		(BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*len);
//	int info = 0;
	double scalar = 1.0;
	int seed = 2;
	int zero = 0;
	int one = 1;
	int lenm1 = len-1;
	int colmaj = BblasColMajor;
	LAPACKE_zlagge(colmaj, len, one, lenm1, zero, &scalar, vec, len, &seed);
}

void random_mat(int m, int n, BBLAS_Complex64_t *mat)
{
//	BBLAS_Complex64_t *work = malloc(sizeof(BBLAS_Complex64_t)*m*n);
//	int info = 0;
	double scalar = 1.0;
	int seed = 2;
	int mm1 = m-1;
	int nn1 = n-1;
	int colmaj = BblasColMajor;
	LAPACKE_zlagge(colmaj, m, n, mm1, nn1, &scalar, mat, m, &seed);
}

BBLAS_Complex64_t randz()
{
	BBLAS_Complex64_t val;
	random_mat(1, 1, &val);
	return val;
}

void set_params_fixed_zgemv(struct zgemv_batchf_example *zgemv_example)
{
	int batch_count = rand() % 50;
	zgemv_example->batch_count = batch_count;
	zgemv_example->trans = rand() % 3+111;
	int m = rand() % 200;
	zgemv_example->m = m;
	int n = rand() % 200;
	zgemv_example->n = n;
	zgemv_example->alpha = randz();
	/* arrayA */
	BBLAS_Complex64_t **arrayA =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	for (int i = 0; i < batch_count; i++)
	{
		arrayA[i] = (BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*m*n);
		random_mat(m, n, arrayA[i]);
	}
	zgemv_example->arrayA = arrayA;
	zgemv_example->lda = zgemv_example->m;
	/* arrayx */
	BBLAS_Complex64_t **arrayx =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	int curdim;
	if (zgemv_example->trans == BblasNoTrans)
	{
		curdim = n;
	}
	else
	{
		curdim = m;
	}
	for (int i = 0; i < batch_count; i++)
	{
		arrayx[i] =
			(BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*curdim);
		random_zvec(curdim, arrayx[i]);
	}
	zgemv_example->arrayx = arrayx;
	zgemv_example->incx = 1;
	zgemv_example->beta = randz();
	/* arrayy */
	BBLAS_Complex64_t **arrayy =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	if (zgemv_example->trans == BblasNoTrans)
	{
		curdim = m;
	}
	else
	{
		curdim = n;
	}
	for (int i = 0; i < batch_count; i++)
	{
		arrayy[i] =
			(BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*curdim);
		random_zvec(curdim, arrayy[i]);
	}
	zgemv_example->arrayy = arrayy;
	zgemv_example->incy = 1;
}

void set_params_variable_zgemv(struct zgemv_batchv_example *zgemv_example)
{
	int batch_count = rand() % 50;
	zgemv_example->batch_count = batch_count;
	int *trans = (int*) malloc(sizeof(int)*batch_count);
	for(int i = 0; i < batch_count; i++)
	{
		trans[i] = rand() % 3 + 111;
	}
	zgemv_example->trans = (enum BBLAS_TRANS*) trans;
	int *m = (int*) malloc(sizeof(int)*batch_count);
	for (int i = 0; i < batch_count; i++)
	{
		m[i] = rand() % 200;
	}
	zgemv_example->m = m;
	int *n = (int*) malloc(sizeof(int)*batch_count);
	for (int i = 0; i < batch_count; i++)
	{
		n[i] = rand() % 200;
	}
	zgemv_example->n = n;
	BBLAS_Complex64_t *alpha = malloc(sizeof(BBLAS_Complex64_t)*batch_count);
	for (int i = 0; i < batch_count; i++)
	{
		alpha[i] = randz();
	}
	zgemv_example->alpha = alpha;
	/* arrayA */
	BBLAS_Complex64_t **arrayA =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	for (int i = 0; i < batch_count; i++)
	{
		arrayA[i] = (BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*m[i]*n[i]);
		random_mat(m[i], n[i], arrayA[i]);
	}
	zgemv_example->arrayA = arrayA;
	int *lda = (int*) malloc(sizeof(int)*batch_count);
	for (int i = 0; i < batch_count; i++)
	{
		lda[i] = m[i];
	}
	zgemv_example->lda = lda;
	/* arrayx */
	BBLAS_Complex64_t **arrayx =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	int curdim;
	for (int i = 0; i < batch_count; i++)
	{
		if (trans[i] == BblasNoTrans)
		{
			curdim = n[i];
		}
		else
		{
			curdim = m[i];
		}
		arrayx[i] =
			(BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*curdim);
		random_zvec(curdim, arrayx[i]);
	}
	zgemv_example->arrayx = arrayx;
	int *incx = (int*) malloc(sizeof(int)*batch_count);
	for (int i = 0; i < batch_count; i++)
	{
		incx[i] = 1;
	}
	zgemv_example->incx = incx;
	BBLAS_Complex64_t *beta = malloc(sizeof(BBLAS_Complex64_t)*batch_count);
	for (int i = 0; i < batch_count; i++)
	{
		beta[i] = randz();
	}
	zgemv_example->beta = beta;
	/* arrayy */
	BBLAS_Complex64_t **arrayy =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	for (int i = 0; i < batch_count; i++)
	{
		if (trans[i] == BblasNoTrans)
		{
			curdim = m[i];
		}
		else
		{
			curdim = n[i];
		}
		arrayy[i] =
			(BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*curdim);
		random_zvec(curdim, arrayy[i]);
	}
	zgemv_example->arrayy = arrayy;
	int *incy = (int*) malloc(sizeof(int)*batch_count);
	for (int i = 0; i < batch_count; i++)
	{
		incy[i] = 1;
	}
	zgemv_example->incy = incy;
}

void set_params_fixed_zgemm(struct zgemm_batchf_example *zgemm_example)
{
    int i;
    int batch_count;
    int m;
    int n;
    int k;
    int lda;
    int ldb;

	batch_count = rand() % 5;
	zgemm_example->batch_count = batch_count;
	zgemm_example->transA = rand() % 3+111;
	zgemm_example->transB = rand() % 3+111;
	m = rand() % 20;
	zgemm_example->m = m;
	n = rand() % 20;
	zgemm_example->n = n;
	k = rand() % 20;
	zgemm_example->k = k;
	zgemm_example->alpha = randz();
	zgemm_example->beta = randz();

    /* arrayA */
	BBLAS_Complex64_t **arrayA =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	for (i = 0; i < batch_count; i++)
	{
		arrayA[i] = (BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*m*k);
		random_mat(m, k, arrayA[i]);
	}
	zgemm_example->arrayA = arrayA;
	if (zgemm_example->transA == BblasNoTrans)
	{
		lda = m;
	}
	else
	{
		lda = k;
	}
	zgemm_example->lda = lda;

    /* arrayB */
	BBLAS_Complex64_t **arrayB =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	for (i = 0; i < batch_count; i++)
	{
		arrayB[i] = (BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*k*n);
		random_mat(k, n, arrayB[i]);
	}
	zgemm_example->arrayB = arrayB;
	if (zgemm_example->transB == BblasNoTrans)
	{
		ldb = k;
	}
	else
	{
		ldb = n;
	}
	zgemm_example->ldb = ldb;

    /* arrayC */
	BBLAS_Complex64_t **arrayC =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	for (i = 0; i < batch_count; i++)
	{
		arrayC[i] = (BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*m*n);
		random_mat(m, n, arrayC[i]);
	}
	zgemm_example->arrayC = arrayC;
	zgemm_example->ldc = m;
}

// One pointer (op) approach

void set_params_fixed_zgemm_stride(struct zgemm_batchf_example_stride *zgemm_example)
{
    int i;
    int batch_count;
    int m;
    int n;
    int k;
    int lda;
    int ldb;
    enum BBLAS_TRANS transA;
    enum BBLAS_TRANS transB;

	batch_count = rand() % 5;
	zgemm_example->batch_count = batch_count;
    transA = BblasNoTrans;
    zgemm_example->transA = transA;
	transB = BblasNoTrans;
    zgemm_example->transB = transB;
	m = rand() % 20;
	zgemm_example->m = m;
	n = rand() % 20;
	zgemm_example->n = n;
	k = rand() % 20;
	zgemm_example->k = k;
	zgemm_example->alpha = randz();
	zgemm_example->beta = randz();

    /* arrayA */
	BBLAS_Complex64_t *arrayA =
		(BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*batch_count*m*k);
	for (i = 0; i < batch_count; i++)
	{
		random_mat(m, k, &arrayA[i*m*k]);
	}
	zgemm_example->arrayA = arrayA;
	if (zgemm_example->transA == BblasNoTrans)
	{
		lda = m;
	}
	else
	{
		lda = k;
	}
	zgemm_example->lda = lda;
	zgemm_example->strideA = m*k;

    /* arrayB */
	BBLAS_Complex64_t *arrayB =
		(BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t*)*batch_count*k*n);
	for (i = 0; i < batch_count; i++)
	{
		random_mat(k, n, &arrayB[i*k*n]);
	}
	zgemm_example->arrayB = arrayB;
	if (zgemm_example->transB == BblasNoTrans)
	{
		ldb = k;
	}
	else
	{
		ldb = n;
	}
	zgemm_example->ldb = ldb;
	zgemm_example->strideB = k*n;

    /* arrayC */
	BBLAS_Complex64_t *arrayC =
		(BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t*)*batch_count*m*n);
	for (i = 0; i < batch_count; i++)
	{
		random_mat(m, n, &arrayC[i*m*n]);
	}
	zgemm_example->arrayC = arrayC;
	zgemm_example->ldc = m;
	zgemm_example->strideC = m*n;
}

void set_params_variable_zgemm(struct zgemm_batchv_example *zgemm_example)
{
    int i;
    int batch_count;
    int lda;
    int ldb;

    batch_count = rand() % 5;
    zgemm_example->batch_count = batch_count;
	for (i = 0; i < batch_count; i++)
	{
	    zgemm_example->transA[i] = rand() % 3+111;
	    zgemm_example->transB[i] = rand() % 3+111;
	    zgemm_example->m[i] = rand() % 20;
	    zgemm_example->n[i] = rand() % 20;
	    zgemm_example->k[i] = rand() % 20;
	    zgemm_example->alpha[i] = randz();
	    zgemm_example->beta[i] = randz();
    }

    /* arrayA */
	BBLAS_Complex64_t **arrayA =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	for (i = 0; i < batch_count; i++)
	{
		arrayA[i] = (BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*
                zgemm_example->m[i]*zgemm_example->k[i]);
		random_mat(zgemm_example->m[i], zgemm_example->k[i], arrayA[i]);
	    if (zgemm_example->transA[i] == BblasNoTrans)
	    {
		    lda = zgemm_example->m[i];
	    }
	    else
	    {
		    lda = zgemm_example->k[i];
	    }
	    zgemm_example->lda[i] = lda;
	}
	zgemm_example->arrayA = arrayA;

    /* arrayB */
	BBLAS_Complex64_t **arrayB =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	for (i = 0; i < batch_count; i++)
	{
		arrayB[i] = (BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*
                zgemm_example->k[i]*zgemm_example->n[i]);
		random_mat(zgemm_example->k[i], zgemm_example->n[i], arrayB[i]);
        if (zgemm_example->transB[i] == BblasNoTrans)
	    {
		    ldb = zgemm_example->k[i];
	    }
	    else
	    {
		    ldb = zgemm_example->n[i];
	    }
	    zgemm_example->ldb[i] = ldb;
	}
	zgemm_example->arrayB = arrayB;

    /* arrayC */
	BBLAS_Complex64_t **arrayC =
		(BBLAS_Complex64_t**) malloc(sizeof(BBLAS_Complex64_t*)*batch_count);
	for (i = 0; i < batch_count; i++)
	{
		arrayC[i] = (BBLAS_Complex64_t*) malloc(sizeof(BBLAS_Complex64_t)*
                zgemm_example->m[i]*zgemm_example->n[i]);
		random_mat(zgemm_example->m[i], zgemm_example->n[i], arrayC[i]);
	    zgemm_example->ldc[i] = zgemm_example->m[i];
	}
	zgemm_example->arrayC = arrayC;

}

max_stars_repo_name: mawussi/BBLAS-group


text:
/*

gcc-mp-4.8 -O3 cudaddottest.c common.c -o cudaddottest -I../../../../netlib/CBLAS -I/usr/local/cuda/include/ -L/usr/local/cuda/lib -lcublas
export DYLD_LIBRARY_PATH=/usr/local/cuda/lib
./cudaddottest  > ../../../results/mac_os_x-x86_64-ddot-cuda_nooh.csv

*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <cblas.h>
#include "common.h"
#include <cublas.h>

void checkStatus(char* message, cublasStatus status) {
    if (status != CUBLAS_STATUS_SUCCESS) {
    	fprintf (stderr, "!!!! %s fail %d\n", message, status);
    	exit(EXIT_FAILURE);
    }
}

long benchmark(int size) {
	long requestStart, requestEnd;
	int incx = 1, incy = 1, n = size;
	double *cuA, *cuB;
	cublasStatus status;

    double* a = random_array(size);
    double* b = random_array(size);

	status = cublasAlloc(n, sizeof(double),(void**)&cuA);
	checkStatus("A", status);
	status = cublasAlloc(n, sizeof(double),(void**)&cuB);
	checkStatus("B", status);

	status = cublasSetVector(n, sizeof(double), a, incx, cuA, incx);
	checkStatus("setA", status);

	status = cublasSetVector(n, sizeof(double), b, incy, cuB, incy);
	checkStatus("setB", status);

	requestStart = currentTimeNanos();

	cublasDdot(n, cuA, incx, cuB, incy);

	requestEnd = currentTimeNanos();

	status = cublasFree(cuA);
	checkStatus("freeA", status);
	status = cublasFree(cuB);
	checkStatus("freeB", status);

    free(a);
    free(b);

    return (requestEnd - requestStart);
  }

main()
{
	cublasStatus status;

	status = cublasInit();
	checkStatus("init",  status);

	srand(time(NULL));

    double factor = 6.0 / 100.0;
    int i, j;
    for (i = 0 ; i < 10 ; i++) {
        for (j = 1 ; j <= 100 ; j++) {
            int size = (int) pow(10.0, factor * j);
            if (size < 10) continue;
            long took = benchmark(size);
            printf("\"%d\",\"%lu\"\n", size, took);
        }
    }
}


text:
#pragma once

#include <src/common/Span.h>
#include <gsl/gsl>

#include <type_traits>

namespace miner {

    //simple replacement for std::vector<char> in situations where huge
    //allocations happen, since the standard doesn't define how
    //much more memory std::vector allocates on reserve or resize.

    //basically an owning gsl::span
    template<class T = uint8_t>
    class DynamicBuffer {

        gsl::owner<T*> owner = nullptr;
        span<T> buffer;

    public:

        DynamicBuffer() = default;
        // only use alignment if a stronger power of two alignment than the natural alignment of T is needed
        DynamicBuffer(size_t size, size_t alignment = alignof(T))
            : owner(reinterpret_cast<T*>(::operator new(size * sizeof(T) + alignment - 1)))
            , buffer(reinterpret_cast<T*>((uintptr_t(owner) + alignment - 1) & ~uintptr_t(alignment - 1))
                , ptrdiff_t(size)) {
            if (!std::is_trivially_constructible<T>::value) {
                new(buffer.data()) T[size];
            }
        }

        ~DynamicBuffer() {
            if (owner) { //if not moved-from
                if (!std::is_trivially_destructible<T>::value) {
                    for (auto &element : buffer) {
                        element.~T();
                    }
                }
                ::operator delete(owner);
            }
        }

        DynamicBuffer(DynamicBuffer &&o) noexcept
            : owner(o.owner)
            , buffer(o.buffer) {
            o.owner = nullptr; //remove ownership
        }

        DynamicBuffer &operator=(DynamicBuffer &&o) noexcept {
            owner = o.owner;
            buffer = o.buffer;
            o.owner = nullptr; //remove ownership
            return *this;
        }

        //copy
        DynamicBuffer(const DynamicBuffer &) = delete;
        DynamicBuffer &operator=(const DynamicBuffer &) = delete;

        T *data() {
            return buffer.data();
        }

        uint8_t *bytes() {
            return reinterpret_cast<uint8_t*>(buffer.data());
        }

        const T *data() const {
            return buffer.data();
        }

        const uint8_t *bytes() const {
            return reinterpret_cast<const uint8_t*>(buffer.data());
        }

        size_t size_bytes() const {
            return static_cast<size_t>(buffer.size_bytes());
        }

        size_t size() const {
            return static_cast<size_t>(buffer.size());
        }

        span<const T> getSpan() const {
            return buffer;
        }

        span<const uint8_t> getByteSpan() const {
            return {bytes(), static_cast<ptrdiff_t>(size_bytes())};
        }

        operator bool() const {
            return owner != nullptr;
        }

    };

}


text:
/* randist/nbinomial.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
 *
 * This program 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.
 *
 * This program 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>

/* The negative binomial distribution has the form,

   prob(k) =  Gamma(n + k)/(Gamma(n) Gamma(k + 1))  p^n (1-p)^k

   for k = 0, 1, ... . Note that n does not have to be an integer.

   This is the Leger's algorithm (given in the answers in Knuth) */

unsigned int
gsl_ran_negative_binomial (const gsl_rng * r, double p, double n)
{
  double X = gsl_ran_gamma (r, n, 1.0) ;
  unsigned int k = gsl_ran_poisson (r, X*(1-p)/p) ;
  return k ;
}

double
gsl_ran_negative_binomial_pdf (const unsigned int k, const double p, double n)
{
  double P;

  double f = gsl_sf_lngamma (k + n) ;
  double a = gsl_sf_lngamma (n) ;
  double b = gsl_sf_lngamma (k + 1.0) ;

  P = exp(f-a-b) * pow (p, n) * pow (1 - p, (double)k);

  return P;
}


text:
#ifndef _PK_SLOPE_H_
#define _PK_SLOPE_H_

#include <fstream>
#include <string>
#include <vector>
#include <cmath>
#include <vector_types.h>
#include <gsl/gsl_deriv.h>
#include "cspline.h"
#include "file_check.h"

template <typename T> class pkSlope{
    int N;
    std::vector<T> k, Pk, n;
    cspline<T> slope;

    void get_pk_spline();

    void calculate_n();

    void get_n_spline();

    public:
        cspline<T> power;

        pkSlope();

        pkSlope(std::string in_pk_file);

        void initialize(std::string in_pk_file);

        void calculate();

        void get_device_spline(std::vector<float4> &d_spline);

        void get_device_spline(std::vector<double4> &d_spline);

};

template <typename T> double f(double x, void *params) {
    cspline<T> Pow = *(cspline<T> *)params;
    return Pow.evaluate(x);
}

template <typename T> void pkSlope<T>::get_pk_spline() {
    pkSlope<T>::power.initialize(pkSlope<T>::k, pkSlope<T>::Pk);
}

template <typename T> void pkSlope<T>::calculate_n() {
    gsl_function F;
    double result, err;

    F.function = &f<T>;
    F.params = &power;

    for (int i = 0; i < pkSlope<T>::N; ++i) {
        if (i == 0) {
            gsl_deriv_forward(&F, pkSlope<T>::k[i], 1e-8, &result, &err);
        } else if (i == pkSlope<T>::N - 1) {
            gsl_deriv_backward(&F, pkSlope<T>::k[i], 1e-8, &result, &err);
        } else {
            gsl_deriv_central(&F, pkSlope<T>::k[i], 1e-8, &result, &err);
        }
        pkSlope<T>::n.push_back(result);
    }
}

template <typename T> void pkSlope<T>::get_n_spline() {
    pkSlope<T>::slope.initialize(pkSlope<T>::k, pkSlope<T>::n);
}

template <typename T> pkSlope<T>::pkSlope() {
    pkSlope<T>::N = 0;
}

template <typename T> pkSlope<T>::pkSlope(std::string in_pk_file) {
    pkSlope<T>::initialize(in_pk_file);
}

template <typename T> void pkSlope<T>::initialize(std::string in_pk_file) {
    if (check_file_exists(in_pk_file)) {
        std::ifstream fin(in_pk_file);
        while(!fin.eof()) {
            T kt, pt;
            fin >> kt >> pt;
            if (!fin.eof()) {
                pkSlope<T>::k.push_back(log10(kt));
                pkSlope<T>::Pk.push_back(log10(pt));
            }
        }
        fin.close();

        pkSlope<T>::N = pkSlope<T>::k.size();
    }
}

template <typename T> void pkSlope<T>::calculate() {
    pkSlope<T>::get_pk_spline();
    pkSlope<T>::calculate_n();
    pkSlope<T>::get_n_spline();
}

template <typename T> void pkSlope<T>::get_device_spline(std::vector<float4> &d_spline) {
    pkSlope<T>::slope.set_pointer_for_device(d_spline);
}

template <typename T> void pkSlope<T>::get_device_spline(std::vector<double4> &d_spline) {
    pkSlope<T>::slope.set_pointer_for_device(d_spline);
}

#endif