text:
/*
 * particle_filters.c
 *
 *  Created on: 18 Sep 2020
 *      Author: heine
 */
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_eigen.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <Accelerate/Accelerate.h>

#include "particle_filters.h"

void MLbootstrapfilter(double* y, double sig_std, double obs_std, int data_length, double x0,
                       int* sample_sizes, double *worst_case_sign_ratio, double* x_hats, double *filtering_time, int* mesh_sizes,  int kl, int ku, double* ab0, double* ab1, int ldab, int* ipiv0, int* ipiv1, int ldb0, int ldb1, int info, double w, double d, double E, double I, double h0, double h1, double L , double* meass, int N_meas) {

  int make_correction = 1; // If one ML telescoping correction is made, if 0 not.
  int N_levels = 2; // hard coded for now
  int N = 0;

  const char prog_bar[51] = "-------------------------------------------------";
  int* minds0 = (int*)malloc(2 * N_meas * sizeof(int));
  int* minds1 = minds0 + N_meas;
  double* obss = (double*)malloc(N_meas*sizeof(double));
  double lower_bound = 0.25, upper_bound = L - 0.25;

  clock_t start = clock();
  int bar_segment_count = 0;

  for (short i = 0; i < N_levels; i++) {
    N += sample_sizes[i];
    printf("N%i = %i, ", i, sample_sizes[i]);
  }
  printf("N = %i, L0 = %i, L1 = %i\n", N, mesh_sizes[0],mesh_sizes[1]);

  double *X = (double*) malloc(2 * N * sizeof(double));
  double W1_tmp;
  double* B0 = (double*) malloc(sample_sizes[0]*mesh_sizes[0]*sizeof(double));
  double* B1 = (double*) malloc(sample_sizes[1]*mesh_sizes[1]*sizeof(double));
  double* pred_meass0 = (double*) malloc(2 * sample_sizes[1] * N_meas * sizeof(double));
  double* pred_meass1 = pred_meass0 + sample_sizes[1] * N_meas;
  double* W = (double*) malloc(2 * N * sizeof(double));
  double* absW = W + N;
  short* signs = (short*) malloc(2 * N * sizeof(short));
  short* signs_res = signs + N;
  int* ind = (int*) malloc(2 * N * sizeof(int));
  int* permutation = ind + N; // needed for uniform resampling

  double x_hat = 0;
  double p_hat = 0;

  double* X0 = X;
  double* X1 = X + sample_sizes[0];
  double* W0 = W;
  double* W1 = W + sample_sizes[0];
  double* X_res = X + N;

  printf(
         "%s\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b",
         prog_bar);
  fflush(stdout);

  /* Make copies of the band matrices as DGBSV destroys them */
  double* a0_copy = (double*) malloc(ldab * mesh_sizes[0] * sizeof(double));
  double* a1_copy = (double*) malloc(ldab * mesh_sizes[1] * sizeof(double));

  /*
   * Initial sample and level indicator
   */
  gsl_rng * rng = gsl_rng_alloc(gsl_rng_taus);
  gsl_rng_set(rng, clock());
  for (int i = 0; i < N; i++) {
    X[i] = gsl_ran_gaussian(rng, sig_std) + x0; // init sample
    W[i] = (double) 1.0 / (double) N; // init weights
    signs_res[i] = 1;
  }

  double normaliser = 0;
  double abs_normaliser = 0;
  //  double diff;
  double ESS, ESSABS;
  //  double tmp_x_hat;
  double positive_mass, negative_mass;

  //  FILE* bpf_out = fopen("mlbpf_out.txt", "w");

  clock_t iteration_timer;
  clock_t resampling_timer;
  clock_t level_timer_start;
  double level_0_particle_cost = 0;
  double level_1_particle_cost = 0;
  double rsamp_time_per_prcl = 0;
  worst_case_sign_ratio[0] = 100;
  double sign_balance;
  double obs_std_scale = 0.8;
  double intercepts[2];
  double slopes[2];

  int res_pos, res_neg;
  int dodb = 0;

  // Mesh point indices for the measurement locations
  for(int i = 0; i < N_meas; i++) {
    minds0[i] = (int)(meass[i] / L * (double) mesh_sizes[0]);
    minds1[i] = (int)(meass[i] / L * (double) mesh_sizes[1]);
  }

  /*
   * FILTER MAIN LOOP
   */
  for (int n = 0; n < data_length; n++) {

    // Restore the destroyed band matrix
    array_copy(ab0, a0_copy, ldab * mesh_sizes[0]);
    array_copy(ab1, a1_copy, ldab * mesh_sizes[1]);

    iteration_timer = clock(); // start the stopwatch
    level_timer_start = clock();

    // Extract the observations for the current time step
    for(int i = 0; i < N_meas; i++)
      obss[i] = y[n + i * data_length];

    /*
     * Level 1 weight calculation
     */

    /* Calculate Level 1 first, so we can estimate the bias of the observation mapping */

    /* Solve the beam for level 0 mesh and sample N1 */
    ebb_solve(sample_sizes[1], X1, mesh_sizes[0], B1, w, d, h0, E, I, L, kl, ku, a0_copy, ldab, ipiv0, ldb0, info, N_meas, minds0, pred_meass0);

    /* Solve the beam for level 1 mesh and sample N1 */
    ebb_solve(sample_sizes[1], X1, mesh_sizes[1], B1, w, d, h1, E, I, L, kl, ku, a1_copy, ldab, ipiv1, ldb1, info, N_meas, minds1, pred_meass1);

    /* Fit a first order LS model for the error between level 0 and level 1 */
    if(sample_sizes[1] > 0){ // do only if N1 > 0
      ls_fit(N_meas, sample_sizes[1], X1, pred_meass0, pred_meass1, intercepts, slopes);
    }

    /* Now that we know the bias, we can evaluate the weights for level 1 */
    double diff;
    for (int j = 0; j < sample_sizes[1]; j++) {
      W1[j] = 1;
      W1_tmp = 1;
      for(int k = 0; k < N_meas; k++){
        /* Level 0 */
        diff = obss[k] - pred_meass0[sample_sizes[1] * k + j] - intercepts[k] - slopes[k] * X1[j];
        W1_tmp *= exp( - diff * diff / (double) 2.0 / (obs_std * obs_std_scale) / (obs_std * obs_std_scale));
        /* Level 1 */
        diff = obss[k] - pred_meass1[sample_sizes[1] * k + j];
        W1[j] *= exp( - diff * diff / (double) 2.0 / obs_std / obs_std );
      }

      W1[j] -= W1_tmp; // take the difference and...
      W1[j] /= (double) sample_sizes[1]; // ... normalise by sample size N1
//      W1[j] = 0; /* uncomment to disable level 1 correction*/
    }

    /*
     Level 0 weight calculation
     */

    array_copy(ab0, a0_copy, ldab * mesh_sizes[0]);

    likelihood(X0, sample_sizes[0], obss, obs_std_scale*obs_std, mesh_sizes[0], kl, ku, a0_copy, ldab, ipiv0, ldb0, info, B0, w, d, h0, E, I, W0, L, minds0, N_meas, dodb, intercepts, slopes);

    // normalise by the sample size (N0)
    for (int i = 0; i < sample_sizes[0]; i++)
      W[i] /= (double)sample_sizes[0];

    if(!make_correction) {
      for(int i = 0; i < sample_sizes[1]; i++) {
        W1[i] = 0;
      }
    }

    /* Joint normalisation of the weights */
    normaliser = 0;
    abs_normaliser = 0;
    for (long i = 0; i < N; i++) {

      W[i] *= (double) signs_res[i]; // Take the sign into account

      absW[i] = fabs(W[i]);
      signs[i] = W[i] > 0 ? 1 : -1;
      normaliser += W[i];
      abs_normaliser += absW[i];
    }

    /* Normalise */
    positive_mass = 0;
    negative_mass = 0;
    for (long i = 0; i < N; i++) {
      W[i] /= normaliser;
      if (W[i] > 0) {
        positive_mass += W[i];
      } else {
        negative_mass += fabs(W[i]);
      }
      absW[i] /= abs_normaliser;
    }
    sign_balance = positive_mass / negative_mass;
    if (sign_balance < worst_case_sign_ratio[0]) {
      worst_case_sign_ratio[0] = sign_balance;
    }

    ESS = 0;
    ESSABS = 0;
    for (long i = 0; i < N; i++) {
      ESS += W[i] * W[i];
      ESSABS += absW[i] * absW[i];
    }
    ESS = (double) 1.0 / ESS;
    ESSABS = (double) 1.0 / ESSABS;

    /* Resample */
    resampling_timer = clock();
    resample(N, absW, ind, rng);
    random_permuter(permutation, N, rng);
    rsamp_time_per_prcl += (double) (clock() - resampling_timer) / CLOCKS_PER_SEC / (double) N;

    normaliser = 0;
    res_pos = 0;
    res_neg = 0;
    for (long i = 0; i < N; i++) {
      X_res[permutation[i]] = X[ind[i]];
      signs_res[permutation[i]] = signs[ind[i]];
      normaliser += signs[ind[i]];
      if(signs[ind[i]]>0) {
        res_pos++;
      }else{
        res_neg++;
      }
    }
    for (long i = 0; i < N; i++) {
      W[i] = (double) signs_res[i] / normaliser;
    }

    /* Calculate the output: posterior mean */
    if(x_hats!=NULL){
      x_hat = 0;
      for (int i = 0; i < N; i++) {
        x_hat += X_res[i] * W[i];
        if(isnan(x_hat)){
          printf("NaN detected!\n X_res[i] = %e, W[i] = %e, i = %i, normaliser = %e + = %e, - = %e\nESS=%e, res_+ =%i, res_- = %i\n",
                 X_res[i], W[i], i,normaliser,positive_mass, negative_mass,ESS,res_pos,res_neg);
          fflush(stdout);
          getchar();
        }
      }
      x_hats[n] = x_hat;

      /* Calculate the output: posterior variance */
      p_hat = 0;
      for (long i = 0; i < N; i++) {
        diff = X_res[i] - x_hat;
        p_hat += diff * diff * W[i];
      }
    }

    /* Mutate */
    for (long i = 0; i < N; i++) {
      X[i] = X_res[i] + gsl_ran_gaussian(rng, sig_std);

      if(X[i] < lower_bound)
        X[i] = lower_bound + lower_bound-X[i];
      if(X[i] > upper_bound)
        X[i] = upper_bound - (X[i]-upper_bound);
    }

    if (floor(n * (double) 50 / data_length) > bar_segment_count) {
      printf("█");
      fflush(stdout);
      bar_segment_count++;
    }

  }

  filtering_time[0] = (double) (clock() - start) / CLOCKS_PER_SEC / (double) data_length;
  filtering_time[1] = level_0_particle_cost / (double) data_length / (double) sample_sizes[0];
  filtering_time[2] = level_1_particle_cost / (double) data_length / (double) sample_sizes[1];
  filtering_time[3] = rsamp_time_per_prcl / (double) data_length / (double) N;
  printf(" %5.5f sec\n", filtering_time[0]);

  free(X);
  free(W);
  free(signs);
  free(ind);
  free(a0_copy);
  free(a1_copy);
  free(B0);
  free(B1);
  free(minds0);
  free(obss);
  free(pred_meass0);
  gsl_rng_free(rng);

}

void bootstrapfilter(double* y, double sig_std, double obs_std, int data_length, double x0,
                     int N, double* x_hats, double* p_hats, double *filtering_time, int mesh_size,	int kl, int ku, double* ab, int ldab, int* ipiv, int ldb, int info, double w, double d, double E, double I, double h , double L, double *meass, int N_meas) {

  const char prog_bar[51] = "-------------------------------------------------";
  clock_t start;
  int bar_segment_count = 0;

  printf("Basic BPF with N = %i\n", N);

  double *X = (double*) malloc(2 * N * sizeof(double));
  double *X_res = X + N;
  double *W = (double*) malloc(N * sizeof(double));
  int *ind = (int*) malloc(N * sizeof(int));
  double* payloads = (double*) malloc(mesh_size * N * sizeof(double));
  // Because DGBSV destroys the AB matrix, we need to make a copy of it
  double* ab_copy = (double*) malloc(ldab * mesh_size * sizeof(double));
  double x_hat = 0;
  double p_hat = 0;
  double lower_bound = 0.25, upper_bound = L - 0.25;
  double intercepts[2] = {0,0};
  double slopes[2] = {0,0};

  printf(
         "%s\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b",
         prog_bar);
  fflush(stdout);

  start = clock();

  /*
   * Initial sample and level indicator
   */
  gsl_rng * rng = gsl_rng_alloc(gsl_rng_taus);
  gsl_rng_set(rng, clock());
  for (int i = 0; i < N; i++) {
    X[i] = gsl_ran_gaussian(rng, sig_std) + x0;
    W[i] = (double) 1.0 / (double) N;
  }

  double normaliser = 0;
  double diff;
  double ESS;
  double cost_per_particle = 0;

  int* minds = (int*)malloc(N_meas * sizeof(int));
  double* obss = (double*)malloc(N_meas * sizeof(double));
  for(int i = 0; i < N_meas; i++){
    minds[i] = (int)(meass[i] / L * (double) mesh_size);
  }

  /* Classic BPF main loop */
  /* --------------------- */
  for (int n = 0; n < data_length; n++) {

    /* Weight calculation */
    normaliser = 0;

    array_copy(ab, ab_copy, ldab * mesh_size);

    /*
     * Weight calculation
     */
    clock_t stopwatch_start = clock();
    for(int i = 0; i < N_meas; i++)
      obss[i] = y[n + i*data_length];

    likelihood(X, N, obss, obs_std, mesh_size, kl, ku, ab_copy, ldab, ipiv, ldb, info, payloads, w, d, h, E, I, W, L, minds, N_meas, 0, intercepts, slopes);
    clock_t stopwatch_stop = clock();

    cost_per_particle += (double)(stopwatch_stop-stopwatch_start)/(double) CLOCKS_PER_SEC;

    // Normalise the weights
    normaliser = 0;
    for (int i = 0; i < N; i++)
      normaliser += W[i];
    for (int i = 0; i < N; i++)
      W[i] /= normaliser;

    ESS = 0;
    for (int i = 0; i < N; i++) {
      ESS += W[i] * W[i];
    }
    ESS = (double) 1.0 / ESS;

    /* Resample */
    resample(N, W, ind, rng);
    for (int i = 0; i < N; i++) {
      X_res[i] = X[ind[i]];
    }

    /* Calculate the output: posterior mean */
    if(x_hats != NULL) {
      x_hat = 0;
      for (int i = 0; i < N; i++) {
        x_hat += X_res[i] / (double) N;
      }
      x_hats[n] = x_hat;
      /* Calculate the output: posterior variance */
      p_hat = 0;
      for (int i = 0; i < N; i++) {
        diff = X_res[i] - x_hat;
        p_hat += diff * diff / (double) (N-1);
      }
      if(p_hats!=NULL)
        p_hats[n] = p_hat;
    }
    /* Mutate */
    for (int i = 0; i < N; i++) {
      X[i] = X_res[i] + gsl_ran_gaussian(rng, sig_std);

      if(X[i] < lower_bound)
        X[i] = lower_bound + lower_bound-X[i];
      if(X[i] > upper_bound)
        X[i] = upper_bound - (X[i]-upper_bound);
    }

    if (floor(n * (double) 50 / data_length) > bar_segment_count) {
      printf("█");
      fflush(stdout);
      bar_segment_count++;
    }
  }
  filtering_time[0] = (double) (clock() - start) / CLOCKS_PER_SEC / (double) data_length;
  filtering_time[1] = cost_per_particle / (double) data_length;
  filtering_time[2] = 0;
  filtering_time[3] = 0;
  printf(" %5.5f sec\n", filtering_time[0]);

  free(X);
  free(W);
  free(ind);
  free(ab_copy);
  free(payloads);
  free(minds);
  free(obss);
  gsl_rng_free(rng);
}

/* Computes a first order LS fit for error between level 0 and level 1 */
void ls_fit(int N_meas, int N1, double* X1, double* pred_meass0, double* pred_meass1, double* intercepts, double* slopes) {

  double a,b,c,d,e,f;

  for(int j = 0; j < N_meas; j++) {

    /* First order polynomial fit to the bias */
    a = 0;
    b = 0;
    d = 0;
    e = 0;
    f = 0;
    for(int i = 0; i < N1; i++) {
      a++;
      b += X1[i];
      d += X1[i] * X1[i];
      e += pred_meass1[j * N1 + i] - pred_meass0[j * N1 + i];
      f += X1[i] * (pred_meass1[j * N1 + i] - pred_meass0[j * N1 + i]);
    }
    c = b;
    intercepts[j] = (e * d - b * f) / (a * d - b * c);
    slopes[j] = (-c * e  + a * f) / (a * d - b * c);
  }
}

void likelihood(double* x, int N, double* y, double std, int mesh_size, int kl, int ku, double* ab, int ldab, int* ipiv, int ldb, int info, double* b, double w, double d, double h, double E, double I, double* likes, double L, int* minds, int N_meas, int dodb, double* intercepts, double* slopes) {

  // Create the instantaneous payload
  create_payload(N, mesh_size, b, w, d, x, h, E, I, L);

  // Solve the Euler-Bernoulli beam for each particle
  dgbsv_(&mesh_size, &kl, &ku, &N, ab, &ldab, ipiv, b, &ldb, &info);

  double diff;
  for (int j = 0; j < N; j++) {
    diff = 1;
    likes[j] = 1;
    for(int k = 0; k < N_meas; k++){
      diff = y[k] - b[j * mesh_size + minds[k]] - x[j]*slopes[k] - intercepts[k];
      likes[j] *= exp( - diff * diff / (double) 2.0 / std / std) ;
    }
  }
}

void random_permuter(int *permutation, int N, gsl_rng *r) {

  for (int i = 0; i < N; i++)
    permutation[i] = i;

  int j;
  int tmp;
  for (int i = N - 1; i > 0; i--) {
    j = (int)gsl_rng_uniform_int(r, i + 1);
    tmp = permutation[j];
    permutation[j] = permutation[i];
    permutation[i] = tmp;
  }

}

void resample(int size, double *w, int *ind, gsl_rng *r) {

  /* Generate the exponentials */
  double *e = (double*) malloc((size + 1) * sizeof(double));
  double g = 0;
  for (int i = 0; i <= size; i++) {
    e[i] = gsl_ran_exponential(r, 1.0);
    g += e[i];
  }
  /* Generate the uniform order statistics */
  double *u = (double *) malloc((size + 1) * sizeof(double));
  u[0] = 0;
  for (int i = 1; i <= size; i++)
    u[i] = u[i - 1] + e[i - 1] / g;

  /* Do the actual sampling with inverse cdf */
  double cdf = w[0];
  int j = 0;
  for (int i = 0; i < size; i++) {
    while (cdf < u[i + 1]) {
      j++;
      cdf += w[j];
    }
    ind[i] = j;
  }

  free(e);
  free(u);
}

void create_payload(int N, int n, double *b, double w, double d, double* x, double h, double E, double I, double L) {

  double g = -9.81; // Gravity acceleration
  double scaler = h * h * h * h / (E * I);
  double slab_thickness = 0.20; // meters
  double density = 7874; // kg/ m^3
  double slab_width = 0.72; // m
  int location_index_centre = 0, location_index_right = 0, location_index_left = 0;
  double deviation = slab_thickness / (double) 8;
  double diff;

  for (long j = 0; j < N;j++) { // iterate over particles

    for(int i = 0; i < n; i++)
      b[j * n + i] = g * w * d * (double)480;

    // Find the index corresponding to the location x
    location_index_centre = (int)(x[j] / L * (double) n);

    // Left end of the slab
    location_index_left = location_index_centre -  ceil(slab_thickness / h / (double) 2);
    // make sure index is non-negative
    location_index_left = location_index_left < 0 ? 0 : location_index_left;

    // Right end of the slab
    location_index_right = location_index_centre +  ceil(slab_thickness / h / (double) 2);
    // make sure index is less than mesh size
    location_index_right = location_index_right > n - 1 ? n - 1 : location_index_right;

    if(location_index_centre < 0) {
      printf("NEGATIVE INDEX!\n");
      fflush(stdout);
    }

    // Add weight to the location
    for(int i = location_index_left; i <= location_index_right ; i++){
      diff = x[j] - i*h;

      b[j * n + i] -=  slab_width * slab_width * density * 9.81 * exp(- diff * diff / (double) 2 / deviation / deviation);
    }
    // Scale appropriately (there is some physical justification to this)
    for(int i = 0; i < n;i++) {
      b[j * n + i] *= scaler;
    }
  }
}

void ebb_solve(int N, double* x, int mesh_size, double* b, double w, double d, double h, double E, double I, double L, int kl, int ku, double *ab, int ldab, int *ipiv, int ldb, int info, int N_meas, int* minds, double* pred_measurements) {

  create_payload(N, mesh_size, b, w, d, x, h, E, I, L);
  dgbsv_(&mesh_size, &kl, &ku, &N, ab, &ldab, ipiv, b, &ldb, &info);

  for (int j = 0; j < N; j++) {
      for(int k = 0; k < N_meas; k++){
        pred_measurements[N * k + j] = b[j * mesh_size + minds[k]];
      }
    }
}

void array_copy(double *src, double* dst, int length) {
  for(int i = 0; i < length; i++)
    dst[i] = src[i];
}

max_stars_repo_name: heinekmp/MLBPF

max_issues_repo_name: heinekmp/MLBPF

max_forks_repo_name: heinekmp/MLBPF


text:
/***************************************************************************
    File                 : FrequencyCountDialog.h
    Project              : QtiPlot
    --------------------------------------------------------------------
    Copyright            : (C) 2008 by Ion Vasilief
    Email (use @ for *)  : ion_vasilief*yahoo.fr
    Description          : Frequency count options dialog

 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *  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 FREQUENCYCOUNTDIALOG_H
#define FREQUENCYCOUNTDIALOG_H

#include <QDialog>
#include <gsl/gsl_vector.h>

class QPushButton;
class DoubleSpinBox;
class Table;

//! Filter options dialog
class FrequencyCountDialog : public QDialog
{
    Q_OBJECT

public:
    FrequencyCountDialog(Table *t, QWidget* parent = 0, Qt::WFlags fl = 0 );
    ~FrequencyCountDialog();

private slots:
    bool apply();
    void accept();

private:
    Table *d_source_table;
    Table *d_result_table;
    QString d_col_name;
    gsl_vector *d_col_values;
	int d_bins;

    QPushButton* buttonApply;
	QPushButton* buttonCancel;
	QPushButton* buttonOk;
	DoubleSpinBox* boxStart;
	DoubleSpinBox* boxEnd;
	DoubleSpinBox* boxStep;
};

#endif

ext: h

max_stars_repo_name: hoehnp/SpaceDesignTool

max_issues_repo_name: hoehnp/SpaceDesignTool

max_forks_repo_name: hoehnp/SpaceDesignTool


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <array>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"
#include "PointLight.h"

namespace Library
{
	class PointLight;
	class ProxyModel;
}

namespace Rendering
{
	class MultiplePointLightsMaterial;

	class MultiplePointLightsDemo final : public Library::DrawableGameComponent
	{
	public:
		MultiplePointLightsDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		MultiplePointLightsDemo(const MultiplePointLightsDemo&) = delete;
		MultiplePointLightsDemo(MultiplePointLightsDemo&&) = default;
		MultiplePointLightsDemo& operator=(const MultiplePointLightsDemo&) = default;
		MultiplePointLightsDemo& operator=(MultiplePointLightsDemo&&) = default;
		~MultiplePointLightsDemo();

		bool AnimationEnabled() const;
		void SetAnimationEnabled(bool enabled);
		void ToggleAnimation();

		float AmbientLightIntensity() const;
		void SetAmbientLightIntensity(float intensity);

		const std::array<Library::PointLight, 4>& PointLights() const;
		void SetPointLight(const Library::PointLight& light, size_t index);

		const size_t SelectedLightIndex() const;
		const Library::PointLight& SelectedLight() const;
		void UpdateSelectedLight(const Library::PointLight& light);
		void SelectLight(size_t index);

		float SpecularPower() const;
		void SetSpecularPower(float power);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		inline static const float RotationRate{ DirectX::XM_PI };

		std::shared_ptr<MultiplePointLightsMaterial> mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::uint32_t mIndexCount{ 0 };
		std::array<std::unique_ptr<Library::ProxyModel>, 4> mProxyModels;
		size_t mSelectedLightIndex{ 0 };
		float mModelRotationAngle{ 0.0f };
		bool mAnimationEnabled{ true };
		bool mUpdateMaterial{ true };
	};
}


text:
//This does CELL (~soma) stage of LSTM (long short-term memory) model.
//This requires each neuron to have 4 input time series, Xc, Xi, Xf, Xo,
//where Xc is the usual (or "cellular") input and Xi, Xf, Xo the inputs for the input, forget, output gates.
//Xc, Xi, Xf, Xo are the output of separate linear IN stages (weights and baises).
//In this version, these are stacked into one matrix
//For dim=0: X = [Xc; Xi; Xf; Xo]; for dim=1: X = [Xc Xi Xf Xo].

//For dim=0, C[:,t] = tanh{Xc[:,t] + Uc*Y[:,t-1]}
//           I[:,t] = sig{Xi[:,t] + Ui*Y[:,t-1]}
//           F[:,t] = sig{Xf[:,t] + Uf*Y[:,t-1]}
//           O[:,t] = sig{Xo[:,t] + Uo*Y[:,t-1]}
//           H[:,t] = I[:,t].*C[:,t] + F[:,t].*H[:,t-1]
//           Y[:,t] = O[:,t].*tanh{H[:,t]}
//with sizes Xc, Xi, Xf, Xo: N x T
//           Uc, Ui, Uf, Uo: N x N
//           Y             : N x T
//
//For dim=1, C[t,:] = tanh{Xc[t,:] + Y[t-1,:]*Uc}
//           I[t,:] = sig{Xi[t,:] + Y[t-1,:]*Ui}
//           F[t,:] = sig{Xf[t,:] + Y[t-1,:]*Uf}
//           O[t,:] = sig{Xo[t,:] + Y[t-1,:]*Uo}
//           H[t,:] = I[t,:].*C[t,:] + F[t,:].*H[t-1,:]
//           Y[t,:] = O[t,:].*tanh{H[t,:]}
//with sizes Xc, Xi, Xf, Xo: T x N
//           Uc, Ui, Uf, Uo: N x N
//           Y             : T x N
//
//where sig is the logistic (sigmoid) nonlinearity = 1/(1+exp(-x)),
//I is the input gate, F is the forget gate, O is the output gate,
//C is the "cell input activation vector",
//H is an intermediate (hidden) vector (sometimes called the "cell state vector"),
//Uc, Ui, Uf, Uo are NxN matrices, and Y is the final output (sometimes called the "hidden state vector").

//Note that, the neurons of a layer are independent only if Uc, Ui, Uf, Uo are diagonal matrices.
//This is only really a CELL (~soma) stage in that case.

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

#ifdef I
    #undef I
#endif

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int lstm_s (float *Y, const float *X, const float *Uc, const float *Ui, const float *Uf, const float *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int lstm_d (double *Y, const double *X, const double *Uc, const double *Ui, const double *Uf, const double *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int lstm_inplace_s (float *X, const float *Uc, const float *Ui, const float *Uf, const float *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int lstm_inplace_d (double *X, const double *Uc, const double *Ui, const double *Uf, const double *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int lstm_s (float *Y, const float *X, const float *Uc, const float *Ui, const float *Uf, const float *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;
    const size_t N2 = 2*N, N3 = 3*N, NT = N*T, NT2 = 2*N*T, NT3 = 3*N*T;
    size_t nT, tN, tN4;

    float *C, *I, *F, *O, *H;
    if (!(C=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(I=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(F=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(O=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm_s: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanhf(X[n]) / (1.0f+expf(-X[N+n]));
                Y[n] = tanhf(H[n]) / (1.0f+expf(-X[N3+n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN4 = 4*tN;
                cblas_scopy((int)N,&X[tN4],1,C,1); cblas_scopy((int)N,&X[tN4+N],1,I,1);
                cblas_scopy((int)N,&X[tN4+N2],1,F,1); cblas_scopy((int)N,&X[tN4+N3],1,O,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Y[tN-N],1,o,C,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Y[tN-N],1,o,I,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Y[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Y[tN+n] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanhf(X[nT]) / (1.0f+expf(-X[NT+nT]));
                Y[nT] = tanhf(H[n]) / (1.0f+expf(-X[NT3+nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[t],(int)T,C,1); cblas_scopy((int)N,&X[NT+t],(int)T,I,1);
                cblas_scopy((int)N,&X[NT2+t],(int)T,F,1); cblas_scopy((int)N,&X[NT3+t],(int)T,O,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Y[t-1],(int)T,o,C,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Y[t-1],(int)T,o,I,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Y[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Y[t+n*T] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanhf(X[nT]) / (1.0f+expf(-X[NT+nT]));
                Y[nT] = tanhf(H[n]) / (1.0f+expf(-X[NT3+nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[t],(int)T,C,1); cblas_scopy((int)N,&X[NT+t],(int)T,I,1);
                cblas_scopy((int)N,&X[NT2+t],(int)T,F,1); cblas_scopy((int)N,&X[NT3+t],(int)T,O,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Y[t-1],(int)T,o,C,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Y[t-1],(int)T,o,I,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Y[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Y[t+n*T] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanhf(X[n]) / (1.0f+expf(-X[N+n]));
                Y[n] = tanhf(H[n]) / (1.0f+expf(-X[N3+n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN4 = 4*tN;
                cblas_scopy((int)N,&X[tN4],1,C,1); cblas_scopy((int)N,&X[tN4+N],1,I,1);
                cblas_scopy((int)N,&X[tN4+N2],1,F,1); cblas_scopy((int)N,&X[tN4+N3],1,O,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Y[tN-N],1,o,C,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Y[tN-N],1,o,I,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Y[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Y[tN+n] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in lstm_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int lstm_d (double *Y, const double *X, const double *Uc, const double *Ui, const double *Uf, const double *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;
    const size_t N2 = 2*N, N3 = 3*N, NT = N*T, NT2 = 2*N*T, NT3 = 3*N*T;
    size_t nT, tN, tN4;

    double *C, *I, *F, *O, *H;
    if (!(C=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(I=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(F=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(O=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm_d: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanh(X[n]) / (1.0+exp(-X[N+n]));
                Y[n] = tanh(H[n]) / (1.0+exp(-X[N3+n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN4 = 4*tN;
                cblas_dcopy((int)N,&X[tN4],1,C,1); cblas_dcopy((int)N,&X[tN4+N],1,I,1);
                cblas_dcopy((int)N,&X[tN4+N2],1,F,1); cblas_dcopy((int)N,&X[tN4+N3],1,O,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Y[tN-N],1,o,C,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Y[tN-N],1,o,I,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Y[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Y[tN+n] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanh(X[nT]) / (1.0+exp(-X[NT+nT]));
                Y[nT] = tanh(H[n]) / (1.0+exp(-X[NT3+nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[t],(int)T,C,1); cblas_dcopy((int)N,&X[NT+t],(int)T,I,1);
                cblas_dcopy((int)N,&X[NT2+t],(int)T,F,1); cblas_dcopy((int)N,&X[NT3+t],(int)T,O,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Y[t-1],(int)T,o,C,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Y[t-1],(int)T,o,I,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Y[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Y[t+n*T] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanh(X[nT]) / (1.0+exp(-X[NT+nT]));
                Y[nT] = tanh(H[n]) / (1.0+exp(-X[NT3+nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[t],(int)T,C,1); cblas_dcopy((int)N,&X[NT+t],(int)T,I,1);
                cblas_dcopy((int)N,&X[NT2+t],(int)T,F,1); cblas_dcopy((int)N,&X[NT3+t],(int)T,O,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Y[t-1],(int)T,o,C,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Y[t-1],(int)T,o,I,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Y[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Y[t+n*T] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanh(X[n]) / (1.0+exp(-X[N+n]));
                Y[n] = tanh(H[n]) / (1.0+exp(-X[N3+n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN4 = 4*tN;
                cblas_dcopy((int)N,&X[tN4],1,C,1); cblas_dcopy((int)N,&X[tN4+N],1,I,1);
                cblas_dcopy((int)N,&X[tN4+N2],1,F,1); cblas_dcopy((int)N,&X[tN4+N3],1,O,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Y[tN-N],1,o,C,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Y[tN-N],1,o,I,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Y[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Y[tN+n] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in lstm_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int lstm_inplace_s (float *X, const float *Uc, const float *Ui, const float *Uf, const float *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;
    const size_t N2 = 2*N, N3 = 3*N, N4 = 4*N, NT = N*T, NT2 = 2*N*T, NT3 = 3*N*T;
    size_t nT, tN, tN4;

    float *C, *I, *F, *O, *H;
    if (!(C=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(I=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(F=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(O=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm_inplace_s: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanhf(X[n]) / (1.0f+expf(-X[N+n]));
                X[n] = tanhf(H[n]) / (1.0f+expf(-X[N3+n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN4 = 4*tN;
                cblas_scopy((int)N,&X[tN4],1,C,1); cblas_scopy((int)N,&X[tN4+N],1,I,1);
                cblas_scopy((int)N,&X[tN4+N2],1,F,1); cblas_scopy((int)N,&X[tN4+N3],1,O,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&X[tN4-N4],1,o,C,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&X[tN4-N4],1,o,I,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&X[tN4-N4],1,o,F,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&X[tN4-N4],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    X[tN4+n] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanhf(X[nT]) / (1.0f+expf(-X[NT+nT]));
                X[nT] = tanhf(H[n]) / (1.0f+expf(-X[NT3+nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[t],(int)T,C,1); cblas_scopy((int)N,&X[NT+t],(int)T,I,1);
                cblas_scopy((int)N,&X[NT2+t],(int)T,F,1); cblas_scopy((int)N,&X[NT3+t],(int)T,O,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&X[t-1],(int)T,o,C,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&X[t-1],(int)T,o,I,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&X[t-1],(int)T,o,F,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&X[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    X[t+n*T] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanhf(X[nT]) / (1.0f+expf(-X[NT+nT]));
                X[nT] = tanhf(H[n]) / (1.0f+expf(-X[NT3+nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[t],(int)T,C,1); cblas_scopy((int)N,&X[NT+t],(int)T,I,1);
                cblas_scopy((int)N,&X[NT2+t],(int)T,F,1); cblas_scopy((int)N,&X[NT3+t],(int)T,O,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&X[t-1],(int)T,o,C,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&X[t-1],(int)T,o,I,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&X[t-1],(int)T,o,F,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&X[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    X[t+n*T] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanhf(X[n]) / (1.0f+expf(-X[N+n]));
                X[n] = tanhf(H[n]) / (1.0f+expf(-X[N3+n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN4 = 4*tN;
                cblas_scopy((int)N,&X[tN4],1,C,1); cblas_scopy((int)N,&X[tN4+N],1,I,1);
                cblas_scopy((int)N,&X[tN4+N2],1,F,1); cblas_scopy((int)N,&X[tN4+N3],1,O,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&X[tN4-N4],1,o,C,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&X[tN4-N4],1,o,I,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&X[tN4-N4],1,o,F,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&X[tN4-N4],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    X[tN4+n] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in lstm_inplace_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int lstm_inplace_d (double *X, const double *Uc, const double *Ui, const double *Uf, const double *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;
    const size_t N2 = 2*N, N3 = 3*N, N4 = 4*N, NT = N*T, NT2 = 2*N*T, NT3 = 3*N*T;
    size_t nT, tN, tN4;

    double *C, *I, *F, *O, *H;
    if (!(C=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(I=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(F=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(O=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm_inplace_d: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanh(X[n]) / (1.0+exp(-X[N+n]));
                X[n] = tanh(H[n]) / (1.0+exp(-X[N3+n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN4 = 4*tN;
                cblas_dcopy((int)N,&X[tN4],1,C,1); cblas_dcopy((int)N,&X[tN4+N],1,I,1);
                cblas_dcopy((int)N,&X[tN4+N2],1,F,1); cblas_dcopy((int)N,&X[tN4+N3],1,O,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&X[tN4-N4],1,o,C,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&X[tN4-N4],1,o,I,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&X[tN4-N4],1,o,F,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&X[tN4-N4],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    X[tN4+n] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanh(X[nT]) / (1.0+exp(-X[NT+nT]));
                X[nT] = tanh(H[n]) / (1.0+exp(-X[NT3+nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[t],(int)T,C,1); cblas_dcopy((int)N,&X[NT+t],(int)T,I,1);
                cblas_dcopy((int)N,&X[NT2+t],(int)T,F,1); cblas_dcopy((int)N,&X[NT3+t],(int)T,O,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&X[t-1],(int)T,o,C,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&X[t-1],(int)T,o,I,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&X[t-1],(int)T,o,F,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&X[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    X[t+n*T] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanh(X[nT]) / (1.0+exp(-X[NT+nT]));
                X[nT] = tanh(H[n]) / (1.0+exp(-X[NT3+nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[t],(int)T,C,1); cblas_dcopy((int)N,&X[NT+t],(int)T,I,1);
                cblas_dcopy((int)N,&X[NT2+t],(int)T,F,1); cblas_dcopy((int)N,&X[NT3+t],(int)T,O,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&X[t-1],(int)T,o,C,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&X[t-1],(int)T,o,I,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&X[t-1],(int)T,o,F,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&X[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    X[t+n*T] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanh(X[n]) / (1.0+exp(-X[N+n]));
                X[n] = tanh(H[n]) / (1.0+exp(-X[N3+n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN4 = 4*tN;
                cblas_dcopy((int)N,&X[tN4],1,C,1); cblas_dcopy((int)N,&X[tN4+N],1,I,1);
                cblas_dcopy((int)N,&X[tN4+N2],1,F,1); cblas_dcopy((int)N,&X[tN4+N3],1,O,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&X[tN4-N4],1,o,C,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&X[tN4-N4],1,o,I,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&X[tN4-N4],1,o,F,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&X[tN4-N4],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    X[tN4+n] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in lstm_inplace_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif

ext: c


text:
// Copyright 2018 Jeremy Mason
//
// Licensed under the Apache License, Version 2.0, <LICENSE-APACHE or
// http://apache.org/licenses/LICENSE-2.0> or the MIT license <LICENSE-MIT or
// http://opensource.org/licenses/MIT>, at your option. This file may not be
// copied, modified, or distributed except according to those terms.

//! \file sb_desc.c
//! Contains functions that actually calculate the spherical Bessel descriptors.

#include <cblas.h>      // dscal
#include <math.h>       // pow
#include <stdint.h>     // uint32_t
#include <stdlib.h>     // abort
#include "sbessel.h"    // _sbessel
#include "sb_desc.h"
#include "sb_matrix.h"  // sb_mat_malloc
#include "sb_structs.h" // sb_vec
#include "sb_utility.h" // SB_CHK_ERR
#include "sb_vector.h"  // sb_vec_calloc
#include "safety.h"

// Lookup tables used in get_radial_basis. Built using function in `tables.c`.
static const double _u_data[10152] = {
  #include "unl.tbl"
};
static const size_t _u_n_max = 140;

static const double _c1_data[10011] = {
  #include "c1.tbl"
};
static const size_t _c1_n_max = 140;

static const double _c2_data[10011] = {
  #include "c2.tbl"
};
static const size_t _c2_n_max = 140;

/// Calculates the radial basis functions for the spherical Bessel descriptors.
/// Numerical efficiency depends heavily on the lookup tables defined above,
/// allowing the function to be reduced to evaluating the recursion relations.
///
/// # Parameters
/// - `gnl`: pointer to a matrix to hold the result
/// - `r_data`: radial coordinates of the atoms
/// - `n_max`: defines the number of descriptors calculated
/// - `l`: order of the spherical Bessel functions
/// - `n_atom`: number of atoms in the environment
/// - `rc`: cutoff radius for the environment
///
/// # Returns
/// A copy of `gnl`
static sb_mat * get_radial_basis(
    sb_mat * gnl,
    double * r_data,
    uint32_t l,
    uint32_t n_atom,
    double rc) {
  // access lookup tables directly
  const double * u_data  = _u_data  + l * (2 * _u_n_max  - l + 5) / 2;
  const double * c1_data = _c1_data + l * (2 * _c1_n_max - l + 3) / 2;
  const double * c2_data = _c2_data + l * (2 * _c2_n_max - l + 3) / 2;

  // gnl->n_cols
  const size_t n_cols = gnl->n_cols;

  // forward declaration of variables without initializations
  size_t n, a;
  double u0, u1, u2, d0, d1, e;
  double * g_data;

  // fnl built in gnl
  g_data = gnl->data;
  for (n = 0; n < n_cols; ++n) {
    for (a = 0; a < n_atom; ++a) {
      g_data[a] = c1_data[n] * _sbessel(l, r_data[a] * u_data[n])
        - c2_data[n] * _sbessel(l, r_data[a] * u_data[n + 1]);
    }
    g_data += n_atom;
  }
  sb_mat_smul(gnl, pow(rc, -1.5));

  // initialize quantities used for recursion
  u1 = SB_SQR(u_data[0]);
  u2 = SB_SQR(u_data[1]);
  d1 = 1.;

  // convert to gnl
  g_data = gnl->data;
  for (n = 1; n < n_cols; ++n) {
    u0 = u1;
    u1 = u2;
    u2 = SB_SQR(u_data[n + 1]);

    e = (u0 * u2) / ((u0 + u1) * (u1 + u2));
    d0 = d1;
    d1 = 1. - e / d0;

    g_data += n_atom;
    cblas_dscal(n_atom, 1. / sqrt(d1), g_data, 1);
    cblas_daxpy(n_atom, sqrt(e / (d1 * d0)), g_data - n_atom, 1, g_data, 1);
  }

  return gnl;
}

/// Calculates spherical Bessel descriptors for the given atomic environment.
/// `desc` should contain space for the descriptors, labelled by (n, l) and
/// ordered lexicographically. `disp` should contain the relative Cartesian
/// coordinates of the surrouding atoms in the format `[x_1, y_1, z_1, ...]`.
/// `weights` should contain the weights used in the construction of the
/// neighbor density function (e.g., `[1., ...]`). `restrict` is not used to
/// help with portability.
///
/// # Parameters
/// - `desc`: pointer to an array to hold the result
/// - `disp`: pointer to an array of the relative displacements
/// - `weights`: pointer to an array of the atomic weights
/// - `rc`: cutoff radius for the environment
/// - `n_atom`: number of atoms in the environment
/// - `n_max`: defines the number of descriptors calculated
///
/// # Returns
/// A copy of `desc`
///
/// # Performance
/// The following preprocessor definitions (usually in `safety.h`) enable
/// various safety checks:
/// - `SAFE_MEMORY`: `desc`, `disp`, and `weights` are not `NULL`
/// - `SAFE_FINITE`: `rc` is nonnegative
///
/// # Warning
/// You are reponsible for ensuring that enough memory is allocated for the
/// relevant arrays, and should expect undefined behavior otherwise. The
/// lengths should be:
/// - `desc`: `(n_max + 1) * (n_max + 2) / 2`
/// - `disp`: at least `3 * n_atom`
/// - `weights`: at least `n_atom`
///
/// # Examples
/// ```
/// #include <stddef.h>     // size_t
/// #include <stdio.h>      // printf
/// #include <stdint.h>     // uint32_t
/// #include "sb_desc.h"    // sb_descriptors
///
/// /// An example where `sb_descriptors()` is used to calculate the spherical
/// /// Bessel descriptors for an atomic environment containing four atoms; the
/// /// first and second descriptors should be `0.031870` and `0.138078`.
/// int main(void) {
///   // Sets the number of descriptors returned.
///   uint32_t n_max = 4;
///
///   // Allocate memory for the result.
///   double desc[15] = { 0. };
///
///   // Number of atoms in the environment.
///   uint32_t n_atom = 4;
///
///   // Displacements to the surrounding atoms in Angstroms.
///   // [x_1, y_1, z_1, ...]
///   double disp[12] = {
///      1.3681827, -1.3103517, -1.3131874,
///     -1.5151760,  1.3360077, -1.3477119,
///     -1.3989598, -1.2973683,  1.3679189,
///      1.2279369,  1.3400378,  1.4797429
///   };
///
///   // Weights for the surrounding atoms.
///   double weights[4] = { 1., 1., 1., 1. };
///
///   // Cutoff radius in Angstroms.
///   double rc = 3.7711;
///
///   sb_descriptors(desc, disp, weights, rc, n_atom, n_max);
///
///   // Output the result
///   // Labelled by (n, l), ordered lexicographically
///   printf("SB descriptors:\n");
///   for (size_t a = 0; a < (n_max + 1) * (n_max + 2) / 2; ++a)
///     printf("%.6f\n", desc[a]);
///
///   printf("Completed successfully!\n");
/// }
/// ```
double * sb_descriptors(
    double * desc_arr,
    double * disp_arr,
    const double * weights_arr,
    const double rc,
    const uint32_t n_atom,
    const uint32_t n_max) {
#ifdef SAFE_MEMORY
  SB_CHK_ERR(!desc_arr, abort(), "sb_descriptors: desc cannot be NULL");
  SB_CHK_ERR(!disp_arr, abort(), "sb_descriptors: disp cannot be NULL");
  SB_CHK_ERR(!weights_arr, abort(), "sb_descriptors: weights cannot be NULL");
#endif
#ifdef SAFE_FINITE
  SB_CHK_ERR(rc < 0., abort(), "sb_descriptors: rc cannot be negative");
#endif
  // Check that n_max is within limit defined by tables
  SB_CHK_ERR(n_max > _u_n_max || n_max > _c1_n_max || n_max > _c2_n_max,
      abort(), "sb_descriptors: n_max above limit defined by lookup tables");

  // Convert raw pointers to sb_vec and sb_mat
  sb_mat * disp = malloc(sizeof(sb_mat));
  SB_CHK_ERR(!disp, abort(), "sb_descriptors: failed to allocate disp");

  disp->n_rows = 3;
  disp->n_cols = n_atom;
  disp->n_elem = 3 * n_atom;
  disp->data   = disp_arr;

  double * data1, * data2, * data3;
  size_t a, b;

  // Calculate radial coordinates
  sb_vec * radius = sb_vec_calloc(n_atom, 'r');
  data1 = disp->data;
  data2 = radius->data;
  for (a = 0; a < n_atom; ++a) {
    for (b = 0; b < 3; ++b) {
      *data2 += SB_SQR(data1[b]);
    }
    data1 += 3;
    data2 += 1;
  }
  sb_vec_sqrt(radius);

  // Normalize displacement vectors
  sb_mat_vdiv(disp, radius, 'c');
  sb_vec_smul(radius, 1. / rc);

  // Calculate angle cosines
  sb_mat * gamma = sb_mat_malloc(n_atom, n_atom);
  sb_mat_mm_mul(gamma, disp, disp, "tn");

  // Legendre polynomials
  sb_mat * * lp = malloc((n_max + 1) * sizeof(sb_mat *));
  SB_CHK_ERR(!lp, abort(), "sb_descriptors: failed to allocate lp");
  for (a = 0; a <= n_max; ++a) {
    lp[a] = sb_mat_malloc(n_atom, n_atom);
  }

  sb_mat_set_all(lp[0], 1.);
  sb_mat_memcpy(lp[1], gamma);
  for (a = 2; a <= n_max; ++a) { // l = a
    sb_mat_memcpy(lp[a], gamma);
    sb_mat_pmul(lp[a], lp[a - 1]);
    sb_mat_smul(lp[a], (2. * a - 1.) / (a - 1.));
    sb_mat_psub(lp[a], lp[a - 2]);
    sb_mat_smul(lp[a], (a - 1.) / a);
  }

  // Include multiplier here to simplify calculation below
  for (a = 0; a <= n_max; ++a) { // l = a
    sb_mat_smul(lp[a], (2. * a + 1.) / 12.566370614359172);
  }

  // Radial basis functions
  sb_mat * * gnl = malloc((n_max + 1) * sizeof(sb_mat *));
  SB_CHK_ERR(!gnl, abort(), "sb_descriptors: failed to allocate gnl");
  for (a = 0; a <= n_max; ++a) { // l = a
    gnl[a] = sb_mat_malloc(n_atom, n_max - a + 1);
    get_radial_basis(gnl[a], radius->data, a, n_atom, rc);
    // Scale by the weights
    for (b = 0; b < n_atom; ++b) {
      cblas_dscal(gnl[a]->n_cols, weights_arr[b], gnl[a]->data + b, n_atom);
    }
  }

  // radius can be used for workspace
  data1 = radius->data;
  for (a = 0; a <= n_max; ++a) { // l = a
    data2 = lp[a]->data;
    data3 = gnl[a]->data;
    for (b = a; b <= n_max; ++b) { // n = b
      cblas_dgemv(CblasColMajor, CblasNoTrans, n_atom, n_atom,
          1., data2, n_atom, data3 + (b - a) * n_atom, 1, 0., data1, 1);
      desc_arr[b * (b + 1) / 2 + a] = cblas_ddot(n_atom, data3 + (b - a) * n_atom, 1, data1, 1);
    }
  }

  // Free memory
  for (a = 0; a <= n_max; ++a) {
    SB_MAT_FREE_ALL(lp[a], gnl[a]);
  }
  SB_FREE_ALL(disp, lp, gnl);
  SB_VEC_FREE_ALL(radius);
  SB_MAT_FREE_ALL(gamma);

  return desc_arr;
}


text:
/**
 * author: Jochen K"upper
 * created: Jan 2002
 * file: pygsl/src/statisticsmodule.c
 * $Id: floatmodule.c,v 1.8 2004/03/24 08:40:45 schnizer Exp $
 *
 * "
 */

#include <Python.h>
#include <gsl/gsl_statistics.h>
#include <pygsl/error_helpers.h>
#include <pygsl/block_helpers.h>

/* include real functions for default data-types (double in C) */

#define STATMOD_WEIGHTED
#define STATMOD_APPEND_PY_TYPE(X) X ## Float32
#define STATMOD_APPEND_PYC_TYPE(X) X ## FLOAT
#define STATMOD_FUNC_EXT(X, Y) X ## _float ## Y
#define STATMOD_PY_AS_C PyFloat_AsDouble
#define STATMOD_C_TYPE float
#include "functions.c"

/* initialization */

PyGSL_STATISTICS_INIT(float, "float")

/*
 * Local Variables:
 * mode: c
 * c-file-style: "Stroustrup"
 * End:
 */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#include <stdio.h>
#include <sys/mman.h>
#include <fcntl.h>
#include <unistd.h>
#include <string.h>
#include <semaphore.h>
#include <signal.h>
#include <sys/stat.h>
#include <sys/times.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "circular_buffer.h"
#include "colorprint.h"

bool running = true;
pid_t pid;
circular_buffer *addr;
double total_wait_time = 0.0;
double total_blocked_time = 0.0;
int total_consumed_messages = 0;

sem_t *sem_mem_id;
sem_t *sem_pro_id;
sem_t *sem_con_id;

void print_stats()
{
  struct tms cpu_times;
  times(&cpu_times);
  printf_color(1, "[cyan]Tiempo de espera total:[/cyan] [yellow]%f s[/yellow]\n", total_wait_time);
  printf_color(1, "[cyan]Tiempo bloqueado total:[/cyan] [yellow]%f s[/yellow]\n", total_blocked_time);
  printf_color(1, "[cyan]Mensajes consumidos:[/cyan] [yellow]%i[/yellow]\n", total_consumed_messages);
  printf_color(1, "[cyan]Tiempo de usuario:[/cyan] [yellow]%ld[/yellow]\n", (long)cpu_times.tms_utime);
  printf_color(1, "[cyan]Tiempo de kernel:[/cyan] [yellow]%ld[/yellow]\n", (long)cpu_times.tms_stime);
}

void exit_by_id()
{
  printf_color(1, "[bb][lw][info][/lw][/bb] Consumidor finalizado por ID (aleatorio). PID: %d.\n", pid);
  print_stats();
}

void exit_by_finalizer()
{
  printf_color(1, "[bb][lw][info][/lw][/bb] Consumidor cerrado por finalizador. PID: %d.\n", pid);
  print_stats();
}

void intHandler(int dummy)
{
  printf("Finalización forzada.\n");
  addr->current_consumers--;
  print_stats();
  exit(0);
}

int main(int argc, char *argv[])
{
  signal(SIGINT, intHandler);

  const gsl_rng_type *T;
  gsl_rng *r;
  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc(T);
  gsl_ran_poisson(r, 5);

  int fd;
  time_t rawtime;
  struct tm *timeinfo;

  // Get argv
  if (argc < 3)
  {
    printf_color(1, "[red][br][lw][error][/lw][/br] No se han dado todos los argumentos necesarios.[/red]\n");
    return -1;
  }
  char buffer_name[strlen(argv[1])];
  strcpy(buffer_name, argv[1]);
  int wait_time = atoi(argv[2]);

  // Initialize strings for semaphore names
  char *sem_mem_name_base = "semaphore_memory_";
  char *sem_con_name_base = "semaphore_consumers_";
  char *sem_prod_name_base = "semaphore_producers_";

  char sem_mem_name[strlen(buffer_name) + strlen(sem_mem_name_base)];
  char sem_con_name[strlen(buffer_name) + strlen(sem_con_name_base)];
  char sem_prod_name[strlen(buffer_name) + strlen(sem_prod_name_base)];

  strcpy(sem_mem_name, sem_mem_name_base);
  strcat(sem_mem_name, buffer_name);
  strcpy(sem_con_name, sem_con_name_base);
  strcat(sem_con_name, buffer_name);
  strcpy(sem_prod_name, sem_prod_name_base);
  strcat(sem_prod_name, buffer_name);

  pid = getpid();

  printf_color(1, "[bb][lw][info][/lw][/bb] Inicializando consumidor. PID: %d.\n", pid);
  printf_color(1, "[bb][lw][info][/lw][/bb] Nombre de buffer dado: %s.\n", buffer_name);
  printf_color(1, "[bb][lw][info][/lw][/bb] Tiempo de espera promedio: %i ms.\n\n", wait_time);

  // get shared memory file descriptor (NOT a file)
  fd = shm_open(buffer_name, O_RDWR, S_IRUSR | S_IWUSR);
  if (fd == -1)
  {
    printf_color(1, "[br][lw][error][/lw][/br] [red]El buffer no se ha inicializado.[/red]\n");
    return 10;
  }

  // Get shared memory size from file descriptor
  struct stat finfo;
  fstat(fd, &finfo);
  off_t shared_memory_size = finfo.st_size;

  // map shared memory to process address space
  addr = mmap(NULL, shared_memory_size, PROT_WRITE, MAP_SHARED, fd, 0);
  if (addr == MAP_FAILED)
  {
    printf_color(1, "[br][lw][error][/lw][/br] [red]No se ha podido asignar suficiente memoria al proceso.[/red]\n");
    return 30;
  }

  // Initialize semaphores
  sem_mem_id = sem_open(sem_mem_name, O_CREAT, 0600, 1);
  if (sem_mem_id == SEM_FAILED)
  {
    perror("SEMAPHORE_MEMORY_SYNC  : [sem_open] Failed\n");
  }

  sem_pro_id = sem_open(sem_prod_name, O_CREAT, 0600, addr->buffer_size);
  if (sem_pro_id == SEM_FAILED)
  {
    perror("SEMAPHORE_MEMORY_SYNC  : [sem_open] Failed\n");
  }

  sem_con_id = sem_open(sem_con_name, O_CREAT, 0600, 0);
  if (sem_con_id == SEM_FAILED)
  {
    perror("SEMAPHORE_MEMORY_SYNC  : [sem_open] Failed\n");
  }

  sem_wait(sem_mem_id);
  addr->current_consumers++;
  addr->total_consumers++;
  sem_post(sem_mem_id);

  int current_slot = 0;
  long long t = 0;

  while (running)
  {
    t = current_timestamp();

    sem_wait(sem_con_id);
    sem_wait(sem_mem_id);

    // Add to total blocked time
    t = current_timestamp() - t;
    total_blocked_time = total_blocked_time + t / 1000.0;

    current_slot = addr->next_message_to_consume;

    cbuffer_message message = consume_message(addr);
    total_consumed_messages++;

    printf_color(1, "\n[bb]--------------------------------------------[/bb]\n");
    printf_color(1, "[cyan]Consumidor PID:[/cyan] [yellow]%i[/yellow]\n", pid);
    printf_color(1, "[cyan]Recibido mensaje del productor PID:[/cyan] [yellow]%i[/yellow]\n", message.producer_id);
    printf_color(1, "[cyan]Leído del slot [yellow]%i[/yellow] del buffer[/cyan]\n", current_slot);
    printf_color(1, "[cyan]Contenido del mensaje:[/cyan] [yellow]%i[/yellow]\n", message.content);
    printf_color(1, "[cyan]Número mágico:[/cyan] [yellow]%i[/yellow]\n", message.random);

    time(&rawtime);
    timeinfo = localtime(&rawtime);
    printf_color(1, "[cyan]Hora actual:[/cyan] [yellow]%s[/yellow]", asctime(timeinfo));
    printf_color(1, "[cyan]Consumidores conectados:[/cyan] [yellow]%i[/yellow]\n", addr->current_consumers);
    printf_color(1, "[cyan]Productores conectados:[/cyan] [yellow]%i[/yellow]\n", addr->current_producers);

    printf_color(1, "[bb]--------------------------------------------[/bb]\n\n");

    sem_post(sem_mem_id);
    sem_post(sem_pro_id); //+1 al semaforo para que consuman

    if (message.type == KILL_CONSUMER)
    {
      running = false;
      sem_wait(sem_con_id);
      addr->current_consumers--;
      sem_post(sem_mem_id);
      exit_by_finalizer();
    }
    else if (message.random == pid % 6)
    {
      running = false;
      sem_wait(sem_con_id);
      addr->current_consumers--;
      addr->consumers_killed_by_id++;
      sem_post(sem_mem_id);
      exit_by_id();
    }
    else
    {
      int random_wait = gsl_ran_poisson(r, wait_time);
      //       int random_wait = gsl_ran_exponential(r,2000);
      printf_color(1, "[green]*********** Esperando %d ms ***********[/green]\n", random_wait);
      total_wait_time = total_wait_time + random_wait / 1000.0;
      //       sleep(random_wait/1000);
      usleep(random_wait * 1000);
    }
  }

  gsl_rng_free(r);
  return 0;
}

max_stars_repo_name: JosephTico/semaforos-bozon-circular

max_issues_repo_name: JosephTico/semaforos-bozon-circular

max_forks_repo_name: JosephTico/semaforos-bozon-circular


text:
// ************************************************
// rand_helpers.h
// authors: Lee Howes and David B. Thomas
//
// Contains support code for the random number
// generation necessary for initialising the
// cuda simulations correctly.
//
// Ziggurat code taken from Marsaglia's
// paper.
// ************************************************

#ifndef __rand_helpers_h
#define __rand_helpers_h

// RNG choices
#define USE_ZIGGURAT 0
#define USE_GSL 0

#include <math.h>
#include <assert.h>

#if USE_GSL
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#endif

unsigned Kiss()
{
	static unsigned z=362436069, w=521288629, jsr=123456789, jcong=380116160;

	z=36969*(z&65535)+(z>>16);
	w=18000*(w&65535)+(w>>16);
	unsigned mwc=(z<<16)+w;
	jsr^=(jsr<<17);
	jsr^=(jsr>>13);
	jsr^=(jsr<<5);
	jcong=69069*jcong+1234567;
	return (mwc^jcong)+jsr;
}

double Rand()
{
	unsigned long long x=Kiss();
	x=(x<<32)|Kiss();
	return x*5.4210108624275221703311375920553e-20;
}

#if USE_GSL
const gsl_rng_type **t, **t0;
gsl_rng *rng;
bool initialised = false;

void initRand()
{
    gsl_rng_env_setup();

    t0 = gsl_rng_types_setup ();

    printf ("Available generators:\n");

    for (t = t0; *t != 0; t++)
    {
        printf ("%s\n", (*t)->name);
        if( strcmp("mt19937_1999", (*t)->name) == 0 ) break;
    }

    rng = gsl_rng_alloc (*t);
}

double RandN()
{
    return gsl_ran_gaussian_ziggurat  (rng, 1.0);

}
#else

#if USE_ZIGGURAT

/* Period parameters */
#define CPU_MT_N 624
#define CPU_MT_M 397
#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
#define CPU_MT_UPPER_MASK 0x80000000UL /* most significant w-r bits */
#define CPU_MT_LOWER_MASK 0x7fffffffUL /* least significant r bits */

// MT
unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;

static unsigned long mt[CPU_MT_N]; /* the array for the state vector  */
static int mti=CPU_MT_N+1; /* mti==N+1 means mt[N] is not initialized */

/* initializes mt[N] with a seed */
void init_genrand(unsigned long s)
{
    mt[0]= s & 0xffffffffUL;
    for (mti=1; mti<CPU_MT_N; mti++) {
        mt[mti] =
	    (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
        /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
        /* In the previous versions, MSBs of the seed affect   */
        /* only MSBs of the array mt[].                        */
        /* 2002/01/09 modified by Makoto Matsumoto             */
        mt[mti] &= 0xffffffffUL;
        /* for >32 bit machines */
    }
}

/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
/* slight change for C++, 2004/2/26 */
void init_by_array(unsigned long init_key[], int key_length)
{
    int i, j, k;
    init_genrand(19650218UL);
    i=1; j=0;
    k = (CPU_MT_N>key_length ? CPU_MT_N : key_length);
    for (; k; k--) {
        mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
          + init_key[j] + j; /* non linear */
        mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
        i++; j++;
        if (i>=CPU_MT_N) { mt[0] = mt[CPU_MT_N-1]; i=1; }
        if (j>=key_length) j=0;
    }
    for (k=CPU_MT_N-1; k; k--) {
        mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
          - i; /* non linear */
        mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
        i++;
        if (i>=CPU_MT_N) { mt[0] = mt[CPU_MT_N-1]; i=1; }
    }

    mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
}

/* generates a random number on [0,0xffffffff]-interval */
unsigned long genrand_int32(void)
{
    unsigned long y;
    static unsigned long mag01[2]={0x0UL, MATRIX_A};
    /* mag01[x] = x * MATRIX_A  for x=0,1 */

    if (mti >= CPU_MT_N) { /* generate N words at one time */
        int kk;

        if (mti == CPU_MT_N+1)   /* if init_genrand() has not been called, */
            init_genrand(5489UL); /* a default initial seed is used */

        for (kk=0;kk<CPU_MT_N-CPU_MT_M;kk++) {
            y = (mt[kk]&CPU_MT_UPPER_MASK)|(mt[kk+1]&CPU_MT_LOWER_MASK);
            mt[kk] = mt[kk+CPU_MT_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }
        for (;kk<CPU_MT_N-1;kk++) {
            y = (mt[kk]&CPU_MT_UPPER_MASK)|(mt[kk+1]&CPU_MT_LOWER_MASK);
            mt[kk] = mt[kk+(CPU_MT_M-CPU_MT_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }
        y = (mt[CPU_MT_N-1]&CPU_MT_UPPER_MASK)|(mt[0]&CPU_MT_LOWER_MASK);
        mt[CPU_MT_N-1] = mt[CPU_MT_M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];

        mti = 0;
    }

    y = mt[mti++];

    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return y;
}

double genrand_real2(void)
{
    return genrand_int32()*(1.0/4294967296.0);
    /* divided by 2^32 */
}

double genrand_real3(void)
{
    return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0);
    /* divided by 2^32 */
}

// ZIGGURAT
#define abs(X) abs((int)X)
//#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
#define SHR3 genrand_int32()
//#define UNI (.5 + (signed) SHR3 * .2328306e-9)
#define UNI genrand_real3()
#define RandN() (hz=SHR3, iz=hz&127, (abs(hz)<kn[iz])? hz*wn[iz] : nfix())

static unsigned long iz,jz,jsr=123456789,kn[128],ke[256];
static long hz; static float wn[128],fn[128], we[256],fe[256];

float nfix(void) { /*provides RNOR if #define cannot */
    const float r = 3.442620f; static float x, y;
    for(;;){ x=hz*wn[iz];
        if(iz==0){
            do{
                x=-log(UNI)*0.2904764;
                y=-log(UNI);
            } while(y+y<x*x);
            return (hz>0)? r+x : -r-x;
        }
        if( fn[iz]+UNI*(fn[iz-1]-fn[iz]) < exp(-.5*x*x) ) return x;
        hz=SHR3; iz=hz&127;if(abs(hz)<kn[iz]) return (hz*wn[iz]);
    }
}

/*--------This procedure sets the seed and creates the tables------*/
void initRand() {
    unsigned long jsrseed = 123456789;
    const double m1 = 2147483648.0, m2 = 4294967296.;
    double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
    double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
    int i; jsr=jsrseed;

    /* Tables for RNOR: */ q=vn/exp(-.5*dn*dn);
    kn[0]=(dn/q)*m1; kn[1]=0;
    wn[0]=q/m1; wn[127]=dn/m1;
    fn[0]=1.; fn[127]=exp(-.5*dn*dn);
    for(i=126;i>=1;i--) {
        dn=sqrt(-2.*log(vn/dn+exp(-.5*dn*dn)));
        kn[i+1]=(dn/tn)*m1; tn=dn;
        fn[i]=exp(-.5*dn*dn); wn[i]=dn/m1;
    }

    /* Tables for REXP */ q = ve/exp(-de);
    ke[0]=(de/q)*m2; ke[1]=0;
    we[0]=q/m2; we[255]=de/m2;
    fe[0]=1.; fe[255]=exp(-de);
    for(i=254;i>=1;i--) {
        de=-log(ve/de+exp(-de));
        ke[i+1]= (de/te)*m2; te=de;
        fe[i]=exp(-de); we[i]=de/m2;
    }

    init_by_array(init, length);
}

#else

void initRand()
{

}

double RandN()
{
	static bool cached=false;
	static double cn;

	if(cached){
		cached=false;
		return cn;
	}

	double a=std::sqrt(-2*std::log(Rand()));
	double b=6.283185307179586476925286766559*Rand();
	cn=std::sin(b)*a;
	cached=true;
	return std::cos(b)*a;
}
#endif // USE_ZIGGURAT
#endif // USE_GSL

double MakeChi2Scale(unsigned N)
{
	const double chic1 = std::sqrt ( std::sqrt (1.0 - 1.0 / N));
	const double chic2 = std::sqrt (1.0 - chic1 * chic1);
	return chic1+chic2*RandN();
}

#endif

max_stars_repo_name: fodinabor/HeCBench

max_issues_repo_name: fodinabor/HeCBench

max_forks_repo_name: fodinabor/HeCBench


text:
#ifdef __cplusplus
extern "C" {
#endif

#ifndef GSL_GAUSS
#define GLS_GAUSS

#include <math.h>
#include <gsl_math.h>
#include <gsl_cdf.h>

double gsl_cdf_ugaussian_P (const double x);
double gsl_cdf_ugaussian_Q (const double x);
double gsl_cdf_gaussian_P (const double x, const double sigma);
double gsl_cdf_gaussian_Q (const double x, const double sigma);

#endif // GLS_GAUSS

#ifdef __cplusplus
}
#endif

max_stars_repo_name: zeeev/wham

max_issues_repo_name: JamesLinus/wham

max_forks_repo_name: JamesLinus/wham


text:
/**
 * @file batchf_zdotu_sub.c
 *
 * Part of API test for Batched BLAS routines.
 *
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date
 *
 * @precisions normal z -> c
 *
 **/

#include <cblas.h>
#include "bblas.h"

#define COMPLEX

void batchf_zdotu_sub(
		  const int n,
		  BBLAS_Complex64_t const * const * x,
		  const int incx,
		  BBLAS_Complex64_t const * const * y,
		  const int incy,
		  BBLAS_Complex64_t *dotu,
		  const int batch_count, int *info)
{
	/* Local variables */
	int first_index = 0;
	int batch_iter = 0;
	char func_name[15] = "batchf_zdotu";

	/*initialize the result */
	for (batch_iter = 0; batch_iter < batch_count; batch_iter++)
	  {
	    dotu[batch_iter] = (BBLAS_Complex64_t)0.0;
	  }

	/* Check input arguments */
	if (batch_count < 0)
	{
		xerbla_batch(func_name, BBLAS_ERR_BATCH_COUNT, -1);
	}

	if (n < 0)
	{
		xerbla_batch(func_name, BBLAS_ERR_N, first_index);
		info[first_index] = BBLAS_ERR_N;

	}

	if (incx < 1)
	{
		xerbla_batch(func_name, BBLAS_ERR_INCX, first_index);
		info[first_index] = BBLAS_ERR_INCX;
	}
	if (incy < 1)
	{
		xerbla_batch(func_name, BBLAS_ERR_INCY, first_index);
		info[first_index] = BBLAS_ERR_INCY;
	}
	/* Call CBLAS */
	for (batch_iter = 0; batch_iter < batch_count; batch_iter++)
	{
	  cblas_zdotu_sub(
		      n,
		      (void *)x[batch_iter],
		      incx,
		      (void *)y[batch_iter],
		      incy,
		      &dotu[batch_iter]);
	  /* Successful */
	} /* End fixed size for loop */
	info[first_index] = BBLAS_SUCCESS;
}
#undef COMPLEX

max_stars_repo_name: mawussi/BBLAS-group


text:
/*
/
/ adkGSL.c
/
/ homebrewed addons to gsl */

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_statistics.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include "../extern/include/f2c.h"
#include "../extern/include/blaswrap.h"
#include "../extern/include/clapack.h"
//#include <vecLib/cblas.h>
//#include <accelerate/accelerate.h>
#include "adkGSL.h"

#define SUM_LOG_THRESHOLD  -10

//gsl_matrix_covariance-- computes sample covariance matrix and stores it in matrix cov
void gsl_matrix_covariance(gsl_matrix *data, gsl_matrix *cov){
	gsl_vector_view a, b;
	size_t i, j;
	double v;
	for (i = 0; i < data->size2; i++) {
 		for (j = 0; j < data->size2; j++) {
 			a = gsl_matrix_column (data, i);
 			b = gsl_matrix_column (data, j);
 			v = gsl_stats_covariance (a.vector.data, a.vector.stride,
 				b.vector.data, b.vector.stride, a.vector.size);
 			gsl_matrix_set (cov, i, j, v);
		}
	}
}

/*gsl_vector_includes() - returns 1 or 0 depending on result */
int gsl_vector_includes(gsl_vector *aVec, double aValue){
  int i;

  for(i = 0; i < aVec->size; i++){
    if (gsl_vector_get(aVec, i) == aValue){
      return(1);
    }
  }
  return(0);
}

double gsl_vector_sum(gsl_vector *aVec, int aVecSize){
	double sum;
	int i;

	sum = 0;
	for (i = 0; i < aVecSize; i++){
		sum += gsl_vector_get(aVec, i);
		}
	return sum;
}

double gsl_vector_dot_product(gsl_vector *vec1, gsl_vector *vec2){
  int i;
  double dot = 0.0;

  assert(vec1->size == vec2->size);
  for(i = 0; i < vec1->size; i++){
    dot += gsl_vector_get(vec1, i) *  gsl_vector_get(vec2, i);
  }
  return(dot);
}

void gsl_vector_outer_product(gsl_vector *vec1, gsl_vector *vec2, gsl_matrix *result){
	int i,j;

	assert(vec1->size == result->size1 || vec2->size == result->size2);
	for(i = 0; i < vec1->size; i++){
		for(j = 0; j < vec2->size; j++){
			gsl_matrix_set(result,i,j,gsl_vector_get(vec1,i) * gsl_vector_get(vec2,j));
		}
	}
}

double gsl_matrix_row_sum(gsl_matrix *mat, int iRow, int rowSize){
	gsl_vector_view aRow;
	double sum = 0;
	int i;

	aRow = gsl_matrix_row(mat, iRow);
	for (i = 0; i< rowSize; i++){
		sum += gsl_vector_get(&aRow.vector, i);
		}
	return sum;
}

/* efficiently compute log of sum of values, which themselves are
   stored as logs: that is, return log(sum_i exp(l_i)).  The largest
   of the elements of l (call it maxval) is factored out, so that
   log(sum_i(exp(l_i))) = maxval + log(1 + sum_i(exp(l_i-maxval))),
   where the new sum is taken over 2 <= i < n.  All of the quantities
   in the exp must be negative, and those smaller than some reasonable
   threshold can be ignored. [Thanks to David Haussler for showing me
   this trick].  This code was adapted from code kindly given to me by Adam Siepel */

double log_sum(gsl_vector *vec) {
  double maxval, expsum;
  int k;

  if (vec->size > 1){
  	gsl_sort_vector(vec);
	gsl_vector_reverse(vec);
	}

  maxval = gsl_vector_get(vec, 0);
  expsum = 1;
  k = 1;

  while (k < vec->size && gsl_vector_get(vec, k) - maxval > SUM_LOG_THRESHOLD){
    expsum += exp(gsl_vector_get(vec, k++) - maxval);
	}
  return maxval + log(expsum);
}

/*the following methods are wrappers for use with LAPACK */

void set_lapack_entry(double *A, int i, int j, int nrows, double val){
  A[j*nrows+i] = val;
}

double get_lapack_entry(double *a, int i, int j,int nrows){
  printf("ha");
  return a[j*nrows+i];
}

/*this is meant to be a conversion utility */
double *gsl_matrix_2_lapack(gsl_matrix *m){
  double *lapack_mat;
  int i, j;

  lapack_mat = malloc(m->size1 * m->size2 * sizeof(double));
  for(i = 0; i < m->size1; i++){
    for(j = 0; j < m->size2; j++){
      set_lapack_entry(lapack_mat, i, j,m->size1, gsl_matrix_get(m, i, j));
    }
  }
  return lapack_mat;
}

gsl_matrix *lapack_2_gsl_matrix(double *A, int nrows, int ncols){
  int i, j;
  gsl_matrix *new;

  new = gsl_matrix_alloc(nrows,ncols);
  for(i=0;i<nrows;i++){
    for(j=0;j<ncols;j++){
      gsl_matrix_set(new,i,j,get_lapack_entry(A,i,j,nrows));
    }
  }
  return new;
}

gsl_matrix *gsl_matrix_power_logs(gsl_matrix *aMatrix, int power){
  gsl_matrix *powerMat, *logMat;
  int i, j;
  double x;

  powerMat = gsl_matrix_power(aMatrix, power);
  logMat = gsl_matrix_alloc(aMatrix->size1, aMatrix->size2);
  for(i = 0; i < aMatrix->size1; i++){
    for(j = 0; j < aMatrix->size2; j++){
      x = gsl_matrix_get(powerMat, i, j);
      gsl_matrix_set(logMat, i, j, log(x));
    }
  }
  gsl_matrix_free(powerMat);
  return(logMat);
}

gsl_matrix *gsl_matrix_power(gsl_matrix *aMatrix, int power){
  gsl_matrix *revect, *levect;
  double *dataMat, *wr, *wi, *vl, *vr, *aMat, *bMat, *cMat, *dMat, *work;
  int lwork, i, j, n ,info;
  assert(aMatrix->size1 == aMatrix->size2);
  n = aMatrix->size1;

  if (power == 1){
    gsl_matrix * solMat = gsl_matrix_alloc(n, n);
    for(i=0; i < n; i++){
      for(j = 0; j< n; j++){
	gsl_matrix_set(solMat, i, j, gsl_matrix_get(aMatrix,i,j));
      }
    }
    return(solMat);
  }

  //allocate space for eigenvector matrices, and other result holders
  levect = gsl_matrix_alloc(n, n);
  revect = gsl_matrix_alloc(n, n);
  wr = malloc(n * sizeof(double));
  wi =  malloc(n * sizeof(double));
  vl = malloc(n * n * sizeof(double));
  vr = malloc(n * n * sizeof(double));
  lwork = 20 * n;
  work = malloc(lwork * sizeof(double));

  //move data to LAPACK (col major) format
  dataMat =  gsl_matrix_2_lapack(aMatrix);

  //run LAPACK routine dgeev
  #ifdef OSX
  dgeev('V','V',n, dataMat, n, wr, wi, vl, n, vr, n, work, lwork);
  #else
  dgeev_('V','V',n, dataMat, n, wr, wi, vl, n, vr, n, work, lwork, info);
  #endif
  //collect eigenvalues into matrices
  for(i=0; i < n; i++){
    for(j = 0; j< n; j++){
      gsl_matrix_set(revect, i, j, vr[j*2+i]);
      gsl_matrix_set(levect, j, i, vl[j*2+i]);
    }
  }

  /* rescale such that revect and levect are inverses */
  /* see Press et al. (Numerical Recipes) for a very clear explanation
     of the relationship between the left and right eigenvectors  */
   for(i = 0; i < n; i++){
    double dotProd = 0.0;
    /* compute dot product of row i in levect and col i in revect */
    for(j = 0; j< n; j++){
      dotProd += gsl_matrix_get(levect, i, j) * gsl_matrix_get(revect,j,i);
      }

      /* now rescale levect */
    for(j = 0; j< n; j++){
      double old = gsl_matrix_get(levect, i, j);
      double scaled = old / dotProd;
      gsl_matrix_set(levect, i, j, scaled);
      }
  }

   //prepare for matrix multiplication
   gsl_matrix * solMat = gsl_matrix_alloc(n, n);
   gsl_matrix * diag = gsl_matrix_alloc(n, n);
   gsl_matrix_set_all(diag, 0);
   for(j = 0; j< n; j++){
     gsl_matrix_set(diag, j,j, pow(wr[j],power));
   }

   aMat = gsl_matrix_2_lapack(diag);
   bMat = gsl_matrix_2_lapack(revect);
   cMat = malloc(n * n * sizeof(double));

   //compute a^t = s^-1 * diag^t * s
   cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, n, n, n, 1.0, aMat, n, bMat, n, 0.0, cMat, n);
   free(bMat);

   bMat = gsl_matrix_2_lapack(levect);
   dMat =  malloc(n * n * sizeof(double));
   cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, n, n, n, 1, cMat, n, bMat, n, 0.0, dMat, n);

   //put results into gsl_matrix
   for(i = 0; i < n; i++){
     for(j = 0; j < n; j++){
       gsl_matrix_set(solMat, i, j, dMat[j*n+i]);
     }
   }

   //cleanup
   gsl_matrix_free(levect);
   gsl_matrix_free(revect);
   gsl_matrix_free(diag);
   free(work);
   free(vl);
   free(vr);
   free(wi);
   free(wr);
   free(aMat);
   free(bMat);
   free(cMat);
   free(dMat);
   free(dataMat);
   return(solMat);
}

//returns the trace of a matrix-- sum of diag
double gsl_matrix_trace(gsl_matrix *x){
	int i;
	double sum = 0.0;
	for(i=0;i<x->size1;i++){
		sum += gsl_matrix_get(x,i,i);
	}
	return(sum);
}

//copies the lower triangle of src in to dest, sets others to zero
void gsl_matrix_lower_tri_copy(gsl_matrix *src, gsl_matrix *dest){
	int i, j;
	assert(src->size1 == dest->size1);
	assert(src->size2 == dest->size2);
	for(i = 0; i < src->size1; i++){
		for(j = 0; j < src->size2; j++){
			if(i<=j){
				gsl_matrix_set(dest,i,j,gsl_matrix_get(src,i,j));
			}
			else{
				gsl_matrix_set(dest,i,j,0.0);
			}
		}
	}
}

//copies the upper triangle of src in to dest, sets others to zero
void gsl_matrix_upper_tri_copy(gsl_matrix *src, gsl_matrix *dest){
	int i, j;
	assert(src->size1 == dest->size1);
	assert(src->size2 == dest->size2);
	for(i = 0; i < src->size1 ;i++){
		for(j = 0; j < src->size2; j++){
			if(i>=j){
				gsl_matrix_set(dest,i,j,gsl_matrix_get(src,i,j));
			}
			else{
				gsl_matrix_set(dest,i,j,0.0);
			}
		}
	}
}

//zeros out everything but lower tri-- in place
void gsl_matrix_lower_tri(gsl_matrix *x){
	int i, j;
	for(i = 0; i < x->size1; i++){
		for(j = 0; j < x->size2; j++){
			if(i>j)	gsl_matrix_set(x,i,j,0.0);
		}
	}
}

//zeros out everything but upper tri-- in place
void gsl_matrix_upper_tri(gsl_matrix *x){
		int i, j;
		for(i = 0; i < x->size1; i++){
			for(j = 0; j < x->size2; j++){
				if(i<j)	gsl_matrix_set(x,i,j,0.0);
			}
		}
}

//fillMatrixFromArray-- fills up a matrix based on nrow and ncol given *double
void fillMatrixFromArray(double *numbers, gsl_matrix *dest, int nrow, int ncol){
	int i, j, count = 0;
	assert(nrow == dest->size1);
	for(i = 0; i < nrow; i++){
		for(j=0;j<ncol;j++){
			gsl_matrix_set(dest,i,j,numbers[count++]);
		}
	}
}

//fillMatrixFromArray-- fills up a matrix based on nrow and ncol given gsl_vector
void fillMatrixFromVector(gsl_vector *numbers, gsl_matrix *dest, int nrow, int ncol){
	int i, j, count = 0;
	assert(nrow == dest->size1);
	for(i = 0; i < nrow; i++){
		for(j=0;j<ncol;j++){
			gsl_matrix_set(dest,i,j,gsl_vector_get(numbers,count++));
		}
	}
}

//fillArrayFromMatrix-- fills up a vector  based gsl_matrix
void fillArrayFromMatrix(gsl_matrix *src, gsl_vector *dest){
	int i, j, count = 0;

	for(i = 0; i < src->size1; i++){
		for(j=0;j<src->size2;j++){
			gsl_vector_set(dest,count++, gsl_matrix_get(src,i,j));
		}
	}
}

//fillMatrixFromCholArray-- fills up a matrix based on nrow and ncol given arraay
// representing the Cholesky decomp
void fillMatrixFromCholArray(double *numbers, gsl_matrix *dest, int nrow, int ncol){
	int i, j, count = 0;

	for(i = 0; i < nrow; i++){
		for(j=0;j<ncol;j++){
			if(i<=j)gsl_matrix_set(dest,i,j,numbers[count++]);
		}
	}
}

void fillCholArrayFromMatrix(gsl_matrix *src, gsl_vector *dest){
	int i, j, count = 0;

	for(i = 0; i < src->size1; i++){
		for(j=0;j<src->size2;j++){
			if(i<=j){
				 gsl_vector_set(dest,count++, gsl_matrix_get(src,i,j));
				//printf("i:%d j: %d val: %f\n",i,j,gsl_matrix_get(src,i,j) );
			}

		}
	}

}

//////////////////////////////////////////////////////////////////////////////////////////////////////////
#ifdef OSX

/*some wrappers for use with the OS X vecLib LAPACK implementation */
/*wrapper for the lapack eigensystem routine dgeev_*/

int dgeev(char jobvl, char jobvr, int n,  double *a, int lda, double *wr, double *wi, double *vl,
	  int ldvl,double *vr, int ldvr, double *work, int lwork){

  __CLPK_integer _n, _lda, _ldvl, _ldvr, _lwork, info;

  _n = (__CLPK_integer) n;
  _lda =  (__CLPK_integer) lda;
  _ldvl =(__CLPK_integer) ldvl;
  _ldvr =(__CLPK_integer) ldvr;
  _lwork =(__CLPK_integer) lwork;

  dgeev_(&jobvl, &jobvr, &_n, a, &_lda, wr, wi, vl, &_ldvl, vr, &_ldvr, work, &_lwork, &info);
  return info;
}

/*wrapper for the lapack eigensystem routine dsyev_
int dsyev(char jobvl, char uplo, int n,  double *a, int lda, double *w, double *work, int lwork){

  __CLPK_integer _n, _lda, _lwork, info;

  _n = (__CLPK_integer) n;
  _lda =  (__CLPK_integer) lda;
  _lwork =(__CLPK_integer) lwork;

  dsyev_(&jobvl, &uplo, &_n, a, &_lda, w, work, &_lwork, &info);
  return info;
}

int dgehrd(int n, int ilo, int ihi, double *a, int lda, double *tau, double *work, int lwork){
  __CLPK_integer _n, _ilo, _ihi, _lda,  _lwork, info;

   _n = (__CLPK_integer) n;
  _lda =  (__CLPK_integer) lda;
  _lwork =(__CLPK_integer) lwork;
  _ilo =(__CLPK_integer) ilo;
  _ihi =(__CLPK_integer) ihi;

  dgehrd_(&_n, &_ilo, &_ihi, a, &_lda, tau, work, &_lwork, &info);
  return info;
}

int dgebal(char job, int n,  double *a, int lda, int ilo, int ihi, double *scale){
  __CLPK_integer _n, _ilo, _ihi, _lda, info;

   _n = (__CLPK_integer) n;
  _lda =  (__CLPK_integer) lda;
  _ilo =(__CLPK_integer) ilo;
  _ihi =(__CLPK_integer) ihi;

  dgebal_(&job, &_n, a, &_lda, &_ilo, &_ihi, scale, &info);
  return info;
}

int dhseqr(char job, char compz, int n, int ilo, int ihi,double *h, int ldh, double *wr, double *wi, double *z, int ldz, double *work, int lwork){
  __CLPK_integer _n, _ilo, _ihi, _ldh, _ldz, _lwork, info;

   _n = (__CLPK_integer) n;
  _ldh =  (__CLPK_integer) ldh;
  _ilo =(__CLPK_integer) ilo;
  _ihi =(__CLPK_integer) ihi;
  _ldz =  (__CLPK_integer) ldz;
  _lwork = (__CLPK_integer) lwork;
  dhseqr_(&job, &compz, &_n, &_ilo, &_ihi, h, &_ldh, wr, wi, z,&_ldz,  work, &_lwork, &info);
  return info;
}
#endif

int dgeev(char jobvl, char jobvr, int n,  double *a, int lda, double *wr, double *wi, double *vl,
          int ldvl,double *vr, int ldvr, double *work, int lwork){

  long int _n, _lda, _ldvl, _ldvr, _lwork, info;

  _n =  n;
  _lda =   lda;
  _ldvl = ldvl;
  _ldvr = ldvr;
  _lwork = lwork;

  dgeev_(&jobvl, &jobvr, &_n, a, &_lda, wr, wi, vl, &_ldvl, vr, &_ldvr, work, &_lwork, &info);
  return info;
}

int dsyev(char jobvl, char uplo, int n,  double *a, int lda, double *w, double *work, int lwork){

  long int _n, _lda, _lwork, info;

  _n =  n;
  _lda =   lda;
  _lwork = lwork;

  dsyev_(&jobvl, &uplo, &_n, a, &_lda, w, work, &_lwork, &info);
  return info;
}

int dgehrd(int n, int ilo, int ihi, double *a, int lda, double *tau, double *work, int lwork){
  long int _n, _ilo, _ihi, _lda,  _lwork, info;

   _n =  n;
  _lda =   lda;
  _lwork = lwork;
  _ilo = ilo;
  _ihi = ihi;

  dgehrd_(&_n, &_ilo, &_ihi, a, &_lda, tau, work, &_lwork, &info);
  return info;
}

int dgebal(char job, int n,  double *a, int lda, int ilo, int ihi, double *scale){
  long int _n, _ilo, _ihi, _lda, info;

   _n =  n;
  _lda =   lda;
  _ilo = ilo;
  _ihi = ihi;

  dgebal_(&job, &_n, a, &_lda, &_ilo, &_ihi, scale, &info);
  return info;
}

int dhseqr(char job, char compz, int n, int ilo, int ihi,double *h, int ldh, double *wr, double *wi, double *z, int ldz, double *work, int lwork){
  long int _n, _ilo, _ihi, _ldh, _ldz, _lwork, info;

   _n =  n;
  _ldh =   ldh;
  _ilo = ilo;
  _ihi = ihi;
  _ldz =   ldz;
  _lwork =  lwork;
  dhseqr_(&job, &compz, &_n, &_ilo, &_ihi, h, &_ldh, wr, wi, z,&_ldz,  work, &_lwork, &info);
  return info;
}
*/
#endif

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


text:
/*****************************************
Emitting C Generated Code
*******************************************/
#include <cblas.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
/**************** Snippet ****************/
void Snippet(int x0) {
  int x1[6] = { 1, 2, 3, 4, 5, 6 };
  int x2[6] = { 6, 5, 4, 3, 2, 1 };
  int* x3 = (int*)malloc(6 * sizeof(int));
  int x4 = 0;
  while (x4 != 6) {
    int x5 = x4;
    x3[x5] = x1[x5] + x2[x5];
    x4 = x4 + 1;
  }
  int x6 = 0;
  while (x6 != 6) {
    printf("%d ", x3[x6]);
    x6 = x6 + 1;
  }
  int* x7 = (int*)malloc(6 * sizeof(int));
  int x8 = 0;
  while (x8 != 6) {
    int x9 = x8;
    x7[x9] = x1[x9] - x2[x9];
    x8 = x8 + 1;
  }
  int x10 = 0;
  while (x10 != 6) {
    printf("%d ", x7[x10]);
    x10 = x10 + 1;
  }
  int* x11 = (int*)malloc(6 * sizeof(int));
  int x12 = 0;
  while (x12 != 6) {
    int x13 = x12;
    x11[x13] = x1[x13] * x2[x13];
    x12 = x12 + 1;
  }
  int x14 = 0;
  while (x14 != 6) {
    printf("%d ", x11[x14]);
    x14 = x14 + 1;
  }
  int* x15 = (int*)malloc(6 * sizeof(int));
  int x16 = 0;
  while (x16 != 6) {
    int x17 = x16;
    x15[x17] = x1[x17] / x2[x17];
    x16 = x16 + 1;
  }
  int x18 = 0;
  while (x18 != 6) {
    printf("%d ", x15[x18]);
    x18 = x18 + 1;
  }
}
/*****************************************
End of C Generated Code
*******************************************/
int main(int argc, char *argv[]) {
  if (argc != 2) {
    printf("usage: %s <arg>\n", argv[0]);
    return 0;
  }
  Snippet(atoi(argv[1]));
  return 0;
}

max_stars_repo_name: Shangyint/lms-clean

max_issues_repo_name: Shangyint/lms-clean

max_forks_repo_name: Shangyint/lms-clean