text:
#ifndef L_1_LASSO_D_H
#define L_1_LASSO_D_H

#include "Matrix.h"
#include "APPROX2.h"
#include "ALM_APG.h"
#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <math.h>

//This class solves problem of the form f(x)+g(x)+ h(Mx) by ASGARD_DL;
// where f(x)= 0
//and g(x)=frac{lambda2}{2}|x|_2+lambda1|x|_1,
// h(x)= lambda3|x- b|_1.

template<typename L, typename D>
class L_1_Lasso_d: public ALM_APG<L, D>
{
private:

  D lambda1;

  D lambda2;

  D lambda3;

  Matrix<L,D> my_M;

  D val_lambda_f;

protected:

public:

  L_1_Lasso_d(const char* Matrix_file,D val_lambda1, D val_lambda2, D val_lambda3)
  :ALM_APG<L,D>(),my_M(Matrix_file)
  {
    lambda1=val_lambda1;
    lambda2=val_lambda2;
    lambda3=val_lambda3;
    val_lambda_f=0;
    this->val_mu_f=0;
    this->val_mu_g=lambda2;
  }

  L get_n(){return my_M.nfeatures;}
  L get_m(){return my_M.nsamples;}

  inline D gradient_of_phi_j(D x1, L i){
    return 0;
  }

  inline D value_of_g_i(D x, L i){
    return lambda2*x*x/2+lambda1*fabs(x);
  }

  inline D value_of_phi_j(D x1, L i){return 0;}

  inline D value_of_h_j(D x, L i){
    return lambda3*fabs(x- my_M.b[i]);
  }

  inline D value_of_h_star_j(D x, L j){
    if (x<= lambda3 && x>= -lambda3){
      return my_M.b[j]*x;
    }
    else{
      cout<< "error in h*(x)"<< endl;
      return std::numeric_limits<double>::max();
    }
  }

  inline D prox_of_g_i(D x1,D x2,D x3, L i){
    return compute_one_step(1.0/x2, x1, x3)-x3;
  }

  D compute_one_step(D tau, D u, D x){
    D new_x;
    if(x>tau*(lambda1+u))
    new_x=(x-tau*(lambda1+u))/(1+lambda2*tau);
    else if(x<tau*(u-lambda1))
    new_x=(x-tau*(u-lambda1))/(1+lambda2*tau);
    else
    new_x=0;
    return new_x;
  }

  // compute min_x h^*(x)+ x2/2*(x- x1)*(x- x1)
  inline D prox_of_h_star_j(D x1, D x2, L j){
    if (x1- my_M.b[j]/x2> lambda3){
      return lambda3;
    }
    else if (x1- my_M.b[j]/x2< -lambda3){
      return -lambda3;
    }
    else{
      return x1- my_M.b[j]/x2;
    }
  }

  inline D feasible_dual(vector<D> & x, vector<D> & y){

    if(lambda2>0)
    {
      return 1;
    }
    else
    {
      D scal=1;
      L l=x.size();
      for(L i=0;i<l;i++)
      if(fabs(x[i]+y[i])>lambda1)
      scal=min(lambda1/fabs(x[i]+y[i]),scal);
      return scal;
    }
  }

  inline D value_of_phistar_i(D x,L i) {return 0;}

  inline D value_of_g_tilde_star(D scal,vector<D> & x, vector<D> & y){

    return 0;
  }
  inline void set_matrix_M(){
    this->data_M=my_M;
  }

  inline void set_matrix_A(){
    this->data_A.nsamples=0;
    this->data_A.nfeatures=this->data_M.nfeatures;
    this->data_A.nnz=0;
    this->data_A.ptr.resize(1,0);
    this->data_A.ptr_t.resize(this->data_M.nfeatures+1,0);
  }

  inline D distance_to_subgradient_of_g(){
    D res=0;
    D tmp;
    D xv;
    for(L i=0;i<this->n;i++){
      tmp=-this->gradient_of_f[i]-this->M_tlambda_s[i];
      xv=this->x_s[i];
      if(xv>0) res+=(tmp-lambda2*xv-lambda1)*(tmp-lambda2*xv-lambda1);
      else if(xv<0) res+=(tmp-lambda2*xv+lambda1)*(tmp-lambda2*xv+lambda1);
      else if(tmp-lambda2*xv>lambda1) res+=(tmp-lambda2*xv-lambda1)*(tmp-lambda2*xv-lambda1);
      else if(tmp-lambda2*xv<-lambda1) res+=(tmp-lambda2*xv+lambda1)*(tmp-lambda2*xv+lambda1);
      else res+=0;
    }
    return sqrt(res);
  }

  void ALM_APG_solver(D beta_0, D epsilon_0,  D omega, vector<D> & x0,vector<D> & y0,L val_tau, L max_nb_outer, L p_N_1, L p_N_2,string filename1, string filename2, D time){
    this->ALM_APG_solve_with_APPROX(beta_0, epsilon_0, omega,x0,y0, val_tau, max_nb_outer,  p_N_1,  p_N_2,  val_lambda_f, filename1,  filename2, time);
  }

};
#endif


text:
#include <mpi.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>
#include <signal.h>
#include <gsl/gsl_rng.h>

#include "allvars.h"
#include "proto.h"

#define TAG_DATAIN  100
#define TAG_DATAOUT 101

static void serial_sort(char *base, size_t nmemb, size_t size, int (*compar) (const void *, const void *));
static void msort_serial_with_tmp(char *base, size_t n, size_t s, int (*compar) (const void *, const void *),
				  char *t);
static void parallel_merge(int master, int ncpu, size_t * nlist, size_t nmax, char *base, size_t nmemb,
			   size_t size, int (*compar) (const void *, const void *));

void parallel_sort(void *base, size_t nmemb, size_t size, int (*compar) (const void *, const void *))
{
  int i, ncpu_in_group, master, groupnr;
  size_t nmax, *nlist;

  serial_sort((char *) base, nmemb, size, compar);

  nlist = (size_t *) mymalloc(NTask * sizeof(size_t));

  MPI_Allgather(&nmemb, sizeof(size_t), MPI_BYTE, nlist, sizeof(size_t), MPI_BYTE, MPI_COMM_WORLD);

  for(i = 0, nmax = 0; i < NTask; i++)
    if(nlist[i] > nmax)
      nmax = nlist[i];

  for(ncpu_in_group = 2; ncpu_in_group <= (1 << PTask); ncpu_in_group *= 2)
    {
      groupnr = ThisTask / ncpu_in_group;

      master = ncpu_in_group * groupnr;

      parallel_merge(master, ncpu_in_group, nlist, nmax, (char *) base, nmemb, size, compar);
    }

  myfree(nlist);
}

void parallel_merge(int master, int ncpu, size_t * nlist, size_t nmax,
		    char *base, size_t nmemb, size_t size, int (*compar) (const void *, const void *))
{
  size_t na, nb, nr;
  int cpua, cpub, cpur;
  char *list_a, *list_b, *list_r;

  if(master + ncpu / 2 >= NTask)	/* nothing to do */
    return;

  if(ThisTask != master)
    {
      if(nmemb)
	{
	  list_r = (char *) mymalloc(nmemb * size);

	  MPI_Request requests[2];

	  MPI_Isend(base, nmemb * size, MPI_BYTE, master, TAG_DATAIN, MPI_COMM_WORLD, &requests[0]);
	  MPI_Irecv(list_r, nmemb * size, MPI_BYTE, master, TAG_DATAOUT, MPI_COMM_WORLD, &requests[1]);
	  MPI_Waitall(2, requests, MPI_STATUSES_IGNORE);

	  memcpy(base, list_r, nmemb * size);
	  myfree(list_r);
	}
    }
  else
    {
      list_a = (char *) mymalloc(nmax * size);
      list_b = (char *) mymalloc(nmax * size);
      list_r = (char *) mymalloc(nmax * size);

      cpua = master;
      cpub = master + ncpu / 2;
      cpur = master;

      na = 0;
      nb = 0;
      nr = 0;

      memcpy(list_a, base, nmemb * size);
      if(nlist[cpub])
	MPI_Recv(list_b, nlist[cpub] * size, MPI_BYTE, cpub, TAG_DATAIN, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

      while(cpur < master + ncpu && cpur < NTask)
	{
	  while(na >= nlist[cpua] && cpua < master + ncpu / 2 - 1)
	    {
	      cpua++;
	      if(nlist[cpua])
		MPI_Recv(list_a, nlist[cpua] * size, MPI_BYTE, cpua, TAG_DATAIN, MPI_COMM_WORLD,
			 MPI_STATUS_IGNORE);
	      na = 0;
	    }
	  while(nb >= nlist[cpub] && cpub < master + ncpu - 1 && cpub < NTask - 1)
	    {
	      cpub++;
	      if(nlist[cpub])
		MPI_Recv(list_b, nlist[cpub] * size, MPI_BYTE, cpub, TAG_DATAIN, MPI_COMM_WORLD,
			 MPI_STATUS_IGNORE);
	      nb = 0;
	    }

	  while(nr >= nlist[cpur])
	    {
	      if(cpur == master)
		memcpy(base, list_r, nr * size);
	      else
		{
		  if(nlist[cpur])
		    MPI_Send(list_r, nlist[cpur] * size, MPI_BYTE, cpur, TAG_DATAOUT, MPI_COMM_WORLD);
		}
	      nr = 0;
	      cpur++;
	    }

	  if(na < nlist[cpua] && nb < nlist[cpub])
	    {
	      if(compar(list_a + na * size, list_b + nb * size) < 0)
		{
		  memcpy(list_r + nr * size, list_a + na * size, size);
		  na++;
		  nr++;
		}
	      else
		{
		  memcpy(list_r + nr * size, list_b + nb * size, size);
		  nb++;
		  nr++;
		}
	    }
	  else if(na < nlist[cpua])
	    {
	      memcpy(list_r + nr * size, list_a + na * size, size);
	      na++;
	      nr++;
	    }
	  else if(nb < nlist[cpub])
	    {
	      memcpy(list_r + nr * size, list_b + nb * size, size);
	      nb++;
	      nr++;
	    }

	}

      myfree(list_r);
      myfree(list_b);
      myfree(list_a);
    }
}

static void serial_sort(char *base, size_t nmemb, size_t size, int (*compar) (const void *, const void *))
{
  const size_t storage = nmemb * size;

  char *tmp = (char *) mymalloc(storage);

  msort_serial_with_tmp(base, nmemb, size, compar, tmp);

  myfree(tmp);
}

static void msort_serial_with_tmp(char *base, size_t n, size_t s, int (*compar) (const void *, const void *),
				  char *t)
{
  char *tmp;
  char *b1, *b2;
  size_t n1, n2;

  if(n <= 1)
    return;

  n1 = n / 2;
  n2 = n - n1;
  b1 = base;
  b2 = base + n1 * s;

  msort_serial_with_tmp(b1, n1, s, compar, t);
  msort_serial_with_tmp(b2, n2, s, compar, t);

  tmp = t;

  while(n1 > 0 && n2 > 0)
    {
      if(compar(b1, b2) < 0)
	{
	  --n1;
	  memcpy(tmp, b1, s);
	  tmp += s;
	  b1 += s;
	}
      else
	{
	  --n2;
	  memcpy(tmp, b2, s);
	  tmp += s;
	  b2 += s;
	}
    }

  if(n1 > 0)
    memcpy(tmp, b1, n1 * s);

  memcpy(base, t, (n - n2) * s);
}

ext: c

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
/* Initialize topic and POV assignments uniformly at random, and set
   various hyper-parameters needed for inference. Assignments must be
   initialized before doing inference for the first time. */

#include <assert.h>
#include <gsl/gsl_rng.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#include "index.h"
#include "parse_mmaps.h"
#include "sample.h"

int main(int argc, char **argv) {
  if (argc != 11) {
    printf("Usage: %s mmap_directory num_topics pov_per_topic num_threads psi_alpha "
	   "psi_beta gamma_alpha gamma_beta beta alpha\n",
           argv[0]);
    exit(1);
  }
  int num_threads = atoi(argv[4]);
  char *endptr;
  create_indexes(argv[1],
		 strtod(argv[5], &endptr),
		 strtod(argv[6], &endptr),
		 strtod(argv[7], &endptr),
		 strtod(argv[8], &endptr),
		 strtod(argv[9], &endptr),
		 strtod(argv[10], &endptr),
		 atoi(argv[2]),
		 atoi(argv[3]),
		 num_threads);
  struct mmap_info mmap_info = open_mmaps_mmap(argv[1]);

  // If revisions are numbered from zero, make sure the first
  // revision gets assigned a null topic/POV
  struct revision_assignment* first_assignment
    = get_revision_assignment(&mmap_info, 0);
  first_assignment->topic = -1;
  first_assignment->pov = -1;

  struct sample_threads sample_threads;
  initialize_threads(&sample_threads, num_threads, &mmap_info);
  resample_null(&sample_threads);
  resample_uniform(&sample_threads);

  destroy_threads(&sample_threads);
  close_mmaps(mmap_info);
}

max_stars_repo_name: allenlavoie/topic-pov

max_issues_repo_name: allenlavoie/topic-pov

max_forks_repo_name: allenlavoie/topic-pov


text:
/*
 * GaussianMixtureModel.h
 * Provides implementation of Gaussian Mixture Model return type.
 *
   Copyright 2017 Michal Gallus

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
*/

#pragma once
#include <vector>
#include <gsl.h>

namespace spectre::unsupervised::gmm
{
/// <summary>
/// Represents a single component of a Gaussian Mixture.
/// </summary>
struct GaussianComponent
{
    /// <summary>
    /// Mean value of the component or peak's location.
    /// </summary>
    double mean;

    /// <summary>
    /// Standard deviation of the component or peak's width.
    /// </summary>
    double deviation;

    /// <summary>
    /// Weight of the component or peak's height.
    /// </summary>
    double weight;
};

/// <summary>
/// Represents a composite of gaussian distribution components
/// that build up a Gaussian Mixture Model.
/// </summary>
struct GaussianMixtureModel
{
    /// <summary>
    /// Constructor used for initialization of m/z, intesities and gaussian components
    /// collections.
    /// </summary>
    /// <param name="mzArray">M/z data shared by all spectra.</param>
    /// <param name="intensities">Mean intensities at each point.</param>
    /// <param name="numberOfComponents">Number of Gaussian Components to be set.</param>
    GaussianMixtureModel(const gsl::span<double> &mzArray,
                         const gsl::span<double> &intensities,
                         const std::vector<GaussianComponent> &&components) :
        originalMzArray(mzArray.begin(), mzArray.end()),
        originalMeanSpectrum(intensities.begin(), intensities.end()),
        components(std::move(components)),
        isMerged(), isNoiseReduced(), mzMergingThreshold() // will be used in future.
    { }

    /// <summary>
    /// Collection of Gaussian components.
    /// </summary>
    const std::vector<GaussianComponent> components;

    /// <summary>
    /// Collection of initially supplied mz values.
    /// </summary>
    const std::vector<double> originalMzArray;

    /// <summary>
    /// Collection of average supplied spectra.
    /// </summary>
    const std::vector<double> originalMeanSpectrum;

    /// <summary>
    /// M/z threshold used in components merging.
    /// </summary>
    const double mzMergingThreshold;

    /// <summary>
    /// Value indicating whether this instance is merged.
    /// </summary>
    const bool isMerged;

    /// <summary>
    /// Gets a value indicating whether this instance is noise components reduced.
    /// </summary>
    const bool isNoiseReduced;
};
}

ext: h

max_stars_repo_name: spectre-team/native-algorithms

max_issues_repo_name: spectre-team/native-algorithms

max_forks_repo_name: spectre-team/native-algorithms


text:
#ifndef FIT_GSL_H
#define FIT_GSL_H

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

class Fit;

//! Structure for fitting data
struct FitData {
  size_t n;// number of points to be fitted (size of X, Y and sigma arrays)
  size_t p;// number of fit parameters
  double * X;// the data to be fitted (abscissae)
  double * Y; // the data to be fitted (ordinates)
  double * sigma; // the weighting data
  Fit *fitter; //pointer to the fitter object (used only for the NonLinearFit class)
};

int expd3_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J);
int expd3_df (const gsl_vector * x, void *params, gsl_matrix * J);
int expd3_f (const gsl_vector * x, void *params, gsl_vector * f);
double expd3_d (const gsl_vector * x, void *params);

int expd2_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J);
int expd2_df (const gsl_vector * x, void *params, gsl_matrix * J);
int expd2_f (const gsl_vector * x, void *params, gsl_vector * f);
double expd2_d (const gsl_vector * x, void *params);

int exp_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J);
int exp_df (const gsl_vector * x, void *params, gsl_matrix * J);
int exp_f (const gsl_vector * x, void *params, gsl_vector * f);
double exp_d (const gsl_vector * x, void *params);

int boltzmann_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J);
int boltzmann_df (const gsl_vector * x, void *params, gsl_matrix * J);
int boltzmann_f (const gsl_vector * x, void *params, gsl_vector * f);
double boltzmann_d (const gsl_vector * x, void *params);

int logistic_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J);
int logistic_df (const gsl_vector * x, void *params, gsl_matrix * J);
int logistic_f (const gsl_vector * x, void *params, gsl_vector * f);
double logistic_d (const gsl_vector * x, void *params);

int gauss_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J);
int gauss_df (const gsl_vector * x, void *params, gsl_matrix * J);
int gauss_f (const gsl_vector * x, void *params,gsl_vector * f);
double gauss_d (const gsl_vector * x, void *params);

int gauss_multi_peak_f (const gsl_vector * x, void *params, gsl_vector * f);
double gauss_multi_peak_d (const gsl_vector * x, void *params);
int gauss_multi_peak_df (const gsl_vector * x, void *params, gsl_matrix * J);
int gauss_multi_peak_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J);

int lorentz_multi_peak_f (const gsl_vector * x, void *params, gsl_vector * f);
double lorentz_multi_peak_d (const gsl_vector * x, void *params);
int lorentz_multi_peak_df (const gsl_vector * x, void *params, gsl_matrix * J);
int lorentz_multi_peak_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J);

int user_f(const gsl_vector * x, void *params, gsl_vector * f);
double user_d(const gsl_vector * x, void *params);
int user_df(const gsl_vector * x, void *params,gsl_matrix * J);
int user_fdf(const gsl_vector * x, void *params,gsl_vector * f, gsl_matrix * J);

#endif

max_stars_repo_name: hoehnp/SpaceDesignTool

max_issues_repo_name: hoehnp/SpaceDesignTool

max_forks_repo_name: hoehnp/SpaceDesignTool


text:
/*
** 1st level statistical inference using LISA
** GLM (general linear modeling) using precoloring
**
** G.Lohmann, MPI-KYB, 2018
*/
#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

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

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>

extern gsl_matrix *PseudoInv(gsl_matrix *A,gsl_matrix *B);
extern void printmat(gsl_matrix *R,char *str);
extern void printvec(gsl_vector *x,char *str);

/*
** general linear model (GLM) using pecoloring
*/
void VGLM(gsl_matrix *Data,gsl_matrix *X,gsl_matrix *XInv,gsl_vector *con,VImage map,VImage zmap)
{
  int i;
  int m = Data->size2;
  int n = con->size;
  gsl_set_error_handler_off();
  gsl_vector *y = gsl_vector_calloc (m);
  gsl_vector *beta = gsl_vector_calloc (n);

  /* compute pseudoinverse */
  XInv = PseudoInv(X,XInv);

  /* main loop */
  VFillImage(zmap,VAllBands,0);
  size_t nvox=0;
  for (nvox=0; nvox < Data->size1; nvox++) {

    y->data = gsl_matrix_ptr(Data,nvox,0);
    gsl_blas_dgemv(CblasNoTrans,1.0,XInv,y,0.0,beta);

    /* contrast image */
    double sum=0;
    for (i=0; i<beta->size; i++) {
      sum += (beta->data[i]*con->data[i]);
    }
    if (gsl_isnan(sum) || gsl_isinf(sum)) continue;
    int b = VPixel(map,0,0,nvox,VShort);
    int r = VPixel(map,0,1,nvox,VShort);
    int c = VPixel(map,0,2,nvox,VShort);
    VPixel(zmap,b,r,c,VFloat) = sum;
  }

  gsl_vector_free(beta);
  gsl_vector_free(y);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#ifndef STRUCTS_H
#define STRUCTS_H

#include <gsl/gsl_randist.h>
#include <stdint.h>

#define N_FEATURES 32

#define BOARD_WIDTH 10
#define BOARD_HEIGHT 20

struct genotype {
    float * feature_weights;
    int * feature_enabled;
};

struct phenotype {
    int fitness;
    int has_fitness;
    struct genotype* genotype;
};

struct population {
    int size;
    struct phenotype** individuals;
};

struct board {
    uint16_t lines[BOARD_HEIGHT];
};

struct tetromino {
    int p_top;
    int p_left;
    int p_right;
    int p_bottom;

    uint16_t lines[4];
};

struct tetromino tetrominos[19];

struct t_last_placement {
    struct tetromino * tetromino;
    int x;
    int y;
    int n_lines_removed;
    int * lines_removed;
};

enum selection {
    TOURNAMENT,
    SUS,
    SIGMA,
};

struct feature {
    char * name;
    int weights;
    float (* function) (struct board *, struct board *, struct t_last_placement *);
};

struct feature features[N_FEATURES];

struct options {
    int feature_enabled[N_FEATURES];
    int enabled_f_indices[N_FEATURES];
    int n_features_enabled;
    int n_weights_enabled;

    int verbose;
    int population_size;
    int tournament_group_size;
    int max_n_generations;
    int crossover_points;
    int elitism;
    int no_log;
    int no_change_duration;
    int reset_volume;
    int n_trials;
    int print_board;
    int n_piece_lookahead;
    int randomization_range;
    int mutation_range;

    float feature_enable_rate;
    float mutation_rate;
    float crossover_rate;
    float tournament_group_random_selection;

    enum selection selection;

    char* log_directory;

    gsl_rng * rng;
};

#endif

max_stars_repo_name: badeball/tetris-ea

max_issues_repo_name: badeball/tetris-ea

max_forks_repo_name: badeball/tetris-ea


text:
#ifndef __HCORE__
#define __HCORE__

#ifdef MKL
  #include <mkl.h>
  #include <mkl_lapack.h>
  //#pragma message("MKL is used")
#else
  #ifdef ARMPL
    #include <armpl.h>
  #else
    #include <cblas.h>
  #endif
  #ifdef LAPACKE_UTILS
    #include <lapacke_utils.h>
  #endif
  #include <lapacke.h>
  //#pragma message("MKL is NOT used")
#endif

#ifndef hcore_min
#define hcore_min(a, b) ((a) < (b) ? (a) : (b))
#endif
#ifndef hcore_max
#define hcore_max(a, b) ((a) < (b) ? (b) : (a))
#endif

#define HCORE_NoTrans         111
#define HCORE_Trans           112
#define HCORE_ConjTrans       113

#define HCORE_Upper           121
#define HCORE_Lower           122
#define HCORE_UpperLower      123

#define HCORE_NonUnit         131
#define HCORE_Unit            132

#define HCORE_Left            141
#define HCORE_Right           142

typedef int  HCORE_enum;
#endif

max_stars_repo_name: ecrc/hcore

max_issues_repo_name: ecrc/hcore

max_forks_repo_name: ecrc/hcore


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#define COMPEARTH_PRIVATE_DET3X3 1
#define COMPEARTH_PRIVATE_CROSS3 1
#define COMPEARTH_PRIVATE_NORM3 1
#define COMPEARTH_PRIVATE_GEM3 1
#define COMPEARTH_PRIVATE_GEMT3 1
#include "compearth.h"
#ifdef COMPEARTH_USE_MKL
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#pragma clang diagnostic ignored "-Wstrict-prototypes"
#endif
#include <mkl_lapacke.h>
//#include <mkl_cblas.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <lapacke.h>
//#include <cblas.h>
#endif

#define LWORK 15

/*!
 * @brief Converts a nearly orthogonal matrix into a numerically orthogonal
 *        matrix.
 *
 * @param[in] n        Number of bases.
 * @param[in] itype    Orthogonalization strategy. \n
 *                       CE_ORTH_SVD orthgonalizes with the SVD. \n
 *                       CE_ORTH_TAPE2012 orthogalizes with Tape 2012c
 *                       Appendix E. \n
 *                       CE_ORTH_QUAT orthogonalizes with quaternions.
 * @param[in] Uin      [3 x 3 x n] set of bases where each [3 x 3] basis is
 *                     in column major order.
 *
 * @param[out] Uout    [3 x 3 x n] set of re-orthgonalized bases where
 *                     each [3 x 3] bais is in column major order.
 * @param[out] dtUin   If NULL then this will be ignored.  \n
 *                     Otherwise, this is an array of dimension [n] holding
 *                     the determinants of the input basis.
 * @param[out] dtUout  If NULL then this will be ignored.  \n
 *                     Otherwise, this is an array of dimension [n] holding
 *                     the determinants of the output basis.
 *
 * @result 0 indicates success.
 *
 * @author Carl Tape and converted to C by Ben Baker.
 *
 */
int compearth_Uorth(const int n,
                    const enum ceOrthoType_enum itype,
                    const double *__restrict__ Uin,
                    double *__restrict__ Uout,
                    double *__restrict__ dtUin,
                    double *__restrict__ dtUout)
{
    double Ut[9] __attribute__((aligned(64)));
    double U[9] __attribute__((aligned(64)));
    double Vt[9] __attribute__((aligned(64)));
    double work[LWORK], s[3], p[3], det, normb, normp;
    int i, ierr;
    bool lwantDetIn;
    const double tol = DBL_EPSILON*100.0;
    // Check inputs
    if (n < 1 || Uin == NULL || Uout == NULL)
    {
        if (n < 1){fprintf(stderr, "%s: No matrices\n", __func__);}
        if (Uin == NULL){fprintf(stderr, "%s: Uin is NULL\n", __func__);}
        if (Uout == NULL){fprintf(stderr, "%s: Uout is NULL\n", __func__);}
        return -1;
    }
    lwantDetIn = false;
    if (dtUin != NULL){lwantDetIn = true;}
    // Orthogonalize basis with SVD
    if (itype == CE_ORTH_SVD)
    {
        for (i=0; i<n; i++)
        {
            // This will frequently be called with matrices that may already
            // be sufficiently orthonormal.
            det = det3x3ColumnMajor(&Uin[9*i]);
            if (fabs(det - 1.0) > 1.e-12)
            {
                // Compute SVD
                memcpy(Ut, &Uin[9*i], 9*sizeof(double));
                ierr = LAPACKE_dgesvd_work(LAPACK_COL_MAJOR, 'A', 'A',
                                           3, 3, Ut, 3, s, U, 3,
                                           Vt, 3, work, LWORK);
                if (ierr != 0)
                {
                    fprintf(stderr, "%s: Error computing Ut\n", __func__);
                    return -1;
                }
                // Compute U*V' - note Vt is already computed by SVD
                gemm3_colMajorNoTransNoTrans(U, Vt, &Uout[9*i]);
                // N.B. I could get the determinant here from the singular
                // values but i'll wait
            }
            else
            {
                memcpy(&Uout[9*i], &Uin[9*i], 9*sizeof(double));
            }
            if (lwantDetIn){dtUin[i] = det;}
        }
    }
    // Orthgonalize with suggestion in TapeTape2012c Appendix E
    else if (itype == CE_ORTH_TAPE2012)
    {
        for (i=0; i<n; i++)
        {
            // This will frequently be called with matrices that may already
            // be sufficiently orthonormal.
            det = det3x3ColumnMajor(&Uin[9*i]);
            if (fabs(det - 1.0) > tol)
            {
                cross3(&Uin[9*i], &Uin[9*i+3], p);
                normb = norm3(&Uin[9*i+3]);
                normp = norm3(p);
                Uout[9*i]   = Uin[9*i];
                Uout[9*i+1] = Uin[9*i+1];
                Uout[9*i+2] = Uin[9*i+2];
                Uout[9*i+3] = Uin[9*i+3]/normb;
                Uout[9*i+4] = Uin[9*i+4]/normb;
                Uout[9*i+5] = Uin[9*i+5]/normb;
                Uout[9*i+6] = Uin[9*i+6]/normp;
                Uout[9*i+7] = Uin[9*i+7]/normp;
                Uout[9*i+8] = Uin[9*i+8]/normp;
            }
            else
            {
                memcpy(&Uout[9*i], &Uin[9*i], 9*sizeof(double));
            }
            if (lwantDetIn){dtUin[i] = det;}
        }
    }
    else if (itype == CE_ORTH_QUAT)
    {
        fprintf(stderr, "%s: Error quaternions not yet programmed\n",
                 __func__);
        return -1;
    }
    else if (itype == CE_NO_ORTH)
    {
        memcpy(Uout, Uin, 9*(size_t) n*sizeof(double));
        //cblas_dcopy(9*n, Uin, 1, Uout, 1);
    }
    else
    {
        fprintf(stderr, "%s: Only itype=1 is programmed\n", __func__);
        return -1;
    }
    //if (dtUin != NULL)
    //{
    //    for (i=0; i<n; i++)
    //    {
    //        dtUin[i] = det3x3ColumnMajor(&Uin[9*i]);
    //    }
    //}
    if (dtUout != NULL)
    {
        for (i=0; i<n; i++)
        {
            dtUout[i] = det3x3ColumnMajor(&Uout[9*i]);
        }
    }
    return 0;
}

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach


text:
/* ieee-utils/make_rep.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_ieee_utils.h>

#include "endian.c"
#include "standardize.c"

static void sprint_nybble(int i, char *s) ;
static void sprint_byte(int i, char *s) ;
static int determine_ieee_type (int non_zero, int exponent, int max_exponent);

/* For the IEEE float format the bits are found from the following
   masks,

   sign      = 0x80000000
   exponent  = 0x7f800000
   mantisssa = 0x007fffff

   For the IEEE double format the masks are,

   sign      = 0x8000000000000000
   exponent  = 0x7ff0000000000000
   mantissa  = 0x000fffffffffffff

   */

void
gsl_ieee_float_to_rep (const float * x, gsl_ieee_float_rep * r)
{
  int e, non_zero;

  union {
    float f;
    struct  {
      unsigned char byte[4] ;
    } ieee ;
  } u;

  u.f = *x ;

  if (little_endian_p())
    make_float_bigendian(&(u.f)) ;

  /* note that r->sign is signed, u.ieee.byte is unsigned */

  if (u.ieee.byte[3]>>7)
    {
      r->sign = 1 ;
    }
  else
    {
      r->sign = 0 ;
    }

  e = (u.ieee.byte[3] & 0x7f) << 1 | (u.ieee.byte[2] & 0x80)>>7 ;

  r->exponent = e - 127 ;

  sprint_byte((u.ieee.byte[2] & 0x7f) << 1,r->mantissa) ;
  sprint_byte(u.ieee.byte[1],r->mantissa + 7) ;
  sprint_byte(u.ieee.byte[0],r->mantissa + 15) ;

  r->mantissa[23] = '\0' ;

  non_zero = u.ieee.byte[0] || u.ieee.byte[1] || (u.ieee.byte[2] & 0x7f);

  r->type = determine_ieee_type (non_zero, e, 255) ;
}

void
gsl_ieee_double_to_rep (const double * x, gsl_ieee_double_rep * r)
{

  int e, non_zero;

  union
  {
    double d;
    struct  {
      unsigned char byte[8];
    } ieee ;
  } u;

  u.d= *x ;

  if (little_endian_p())
    make_double_bigendian(&(u.d)) ;

  /* note that r->sign is signed, u.ieee.byte is unsigned */

  if (u.ieee.byte[7]>>7)
    {
      r->sign = 1 ;
    }
  else
    {
      r->sign = 0 ;
    }

  e =(u.ieee.byte[7] & 0x7f)<<4 ^ (u.ieee.byte[6] & 0xf0)>>4 ;

  r->exponent = e - 1023 ;

  sprint_nybble(u.ieee.byte[6],r->mantissa) ;
  sprint_byte(u.ieee.byte[5],r->mantissa + 4) ;
  sprint_byte(u.ieee.byte[4],r->mantissa + 12) ;
  sprint_byte(u.ieee.byte[3],r->mantissa + 20) ;
  sprint_byte(u.ieee.byte[2],r->mantissa + 28) ;
  sprint_byte(u.ieee.byte[1],r->mantissa + 36) ;
  sprint_byte(u.ieee.byte[0],r->mantissa + 44) ;

  r->mantissa[52] = '\0' ;

  non_zero = (u.ieee.byte[0] || u.ieee.byte[1] || u.ieee.byte[2]
              || u.ieee.byte[3] || u.ieee.byte[4] || u.ieee.byte[5]
              || (u.ieee.byte[6] & 0x0f)) ;

  r->type = determine_ieee_type (non_zero, e, 2047) ;
}

/* A table of character representations of nybbles */

static char nybble[16][5]={ /* include space for the \0 */
  "0000", "0001", "0010", "0011",
  "0100", "0101", "0110", "0111",
  "1000", "1001", "1010", "1011",
  "1100", "1101", "1110", "1111"
}  ;

static void
sprint_nybble(int i, char *s)
{
  char *c ;
  c=nybble[i & 0x0f ];
  *s=c[0] ;  *(s+1)=c[1] ;  *(s+2)=c[2] ;  *(s+3)=c[3] ;
}

static void
sprint_byte(int i, char *s)
{
  char *c ;
  c=nybble[(i & 0xf0)>>4];
  *s=c[0] ;  *(s+1)=c[1] ;  *(s+2)=c[2] ;  *(s+3)=c[3] ;
  c=nybble[i & 0x0f];
  *(s+4)=c[0] ;  *(s+5)=c[1] ;  *(s+6)=c[2] ;  *(s+7)=c[3] ;
}

static int
determine_ieee_type (int non_zero, int exponent, int max_exponent)
{
  if (exponent == max_exponent)
    {
      if (non_zero)
        {
          return GSL_IEEE_TYPE_NAN ;
        }
      else
        {
          return GSL_IEEE_TYPE_INF ;
        }
    }
  else if (exponent == 0)
    {
      if (non_zero)
        {
          return GSL_IEEE_TYPE_DENORMAL ;
        }
      else
        {
          return GSL_IEEE_TYPE_ZERO ;
        }
    }
  else
    {
      return GSL_IEEE_TYPE_NORMAL ;
    }
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*
** centrality maps using spectral coherence
**
** get cross-periodogram via multiplication of FFTs
**
** Ref: T. Subba Rao,
** On the Cross Periodogram of a Stationary Gaussian Vector Process,
** The Annals of Mathematical Statistics, Vol. 38, No. 2 (Apr., 1967),
**  pp. 593-597.  (2239172.pdf)
**
** Implementation follows SAS, see
**  http://www.tau.ac.il/cc/pages/docs/sas8/ets/chap17/
**
** G.Lohmann, September 2009
*/
#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_math.h>

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

#define PI 3.14159265

#define SQR(x) ((x) * (x))
#define ABS(x) ((x) > 0 ? (x) : -(x))

#ifdef _OPENMP
#include <omp.h>
#endif /*_OPENMP*/

/* re-implementation of cblas_sspmv,  cblas_sspmv causes problems */
void my_sspmv(float *A,float *x,float *y,int n)
{
  size_t i,j,k,kk;
  float tmp1=0,tmp2=0;

  for (j=0; j<n; j++) y[j] = 0;
  kk = k = 0;
  for (j=0; j<n; j++) {
    tmp1 = x[j];
    tmp2 = 0;
    k = kk;
    for (i=0; i<j; i++) {
      y[i] += tmp1 * A[k];
      tmp2 += A[k] * x[i];
      k++;
    }
    y[j] += tmp1*A[kk+j] + tmp2;
    kk += (j+1);
  }
}

double Normalize(double *data,int n)
{
  int i;
  double sum,nx=(double)n;

  sum = 0;
  for (i=0; i<n; i++) sum += data[i];
  sum /= nx;

  for (i=0; i<n; i++) data[i] -= sum;
  return sum;
}

double Tukey(double x,double wx)
{
  if (ABS(x) >= wx) return 0;
  else return 0.5*(1.0 + cos((PI*x)/wx));
}

int GetFFtIndex(int n,double tr,double wavelength)
{
  int k0;
  double nx,cycle;

  nx = (double)n;
  cycle = wavelength / tr;
  k0 = (int)(nx/cycle + 0.5);
  return k0;
}

double *Weights(int n,int k0,int p)
{
  int k;
  double *w=NULL;
  double kx,wx,sum;

  w  = (double *) VCalloc(n,sizeof(double));
  wx = (double) p;
  kx = 0;
  for (k=0; k<p; k++) {
    w[k] = Tukey(kx,wx);
    kx++;
  }
  sum = 0;
  for (k=0; k<p; k++) sum += w[k];
  sum *= (2.0*PI);
  for (k=0; k<p; k++) w[k] /= sum;
  return w;
}

void Tables(gsl_matrix *tablecos,gsl_matrix *tablesin)
{
  int i,k,n;
  double kx,ix,nx,w;

  n = tablecos->size1;
  nx = (double)n;
  kx = 0;
  for (k=0; k<n; k++) {
    w = 2.0*PI*kx/nx;
    ix = 0;
    for (i=1; i<n; i++) {
      gsl_matrix_set(tablecos,k,i,cos(w*ix));
      gsl_matrix_set(tablesin,k,i,sin(w*ix));
      ix++;
    }
    kx++;
  }
}

void FFT(double *in,double *xcos,double *xsin,gsl_matrix *tablecos,gsl_matrix *tablesin)
{
  int i,k,n;
  double sum1,sum2,kx,ix,nx;

  n  = tablecos->size1;
  nx = 1;
  kx = 0;
  for (k=0; k<n; k++) {
    if (k >= tablecos->size1) VError(" k= %d",k);
    if (k >= tablesin->size1) VError(" k= %d",k);
    sum1 = sum2 = 0;
    /* w = 2.0*PI*kx/nx; */
    ix = 0;
    for (i=0; i<n; i++) {
      sum1 += in[i]*gsl_matrix_get(tablecos,k,i);
      sum2 += in[i]*gsl_matrix_get(tablesin,k,i);
      ix++;
    }
    xcos[k] = 2.0*sum1/nx;
    xsin[k] = 2.0*sum2/nx;
    kx++;
  }
}

void DegreeCentrality(float *A,float *ev,size_t n)
{
  size_t i,j,k;
  double sum=0,nx=(double)n;
  fprintf(stderr," degree centrality, n= %lu\n",n);

  for (i=0; i<n; i++) {
    sum=0;
    for (j=0; j<n; j++) {
      if (i==j) continue;
      if (j < i) k=j+i*(i+1)/2;
      else k=i+j*(j+1)/2;
      sum += A[k];
    }
    ev[i] = (float)(100.0*sum/nx);
  }
}

void EigenvectorCentrality(float *A,float *ev,size_t n)
{
  size_t i,iter,maxiter;
  float sum,d,nx;
  static float *y=NULL;

  fprintf(stderr," power iteration, n= %lu\n",n);
  if (y==NULL) y = (float *) VCalloc(n,sizeof(float));
  if (!y) VError(" err allocating vector ");
  nx = sqrt((double)n);
  for (i=0; i<n; i++) ev[i] = y[i] = 1.0/nx;

  maxiter=50;
  for (iter=0; iter < maxiter; iter++) {

    /* y = Ax,  A symmetric and lower-triangular */
    /* cblas_sspmv(CblasRowMajor,CblasLower,(int)n,1.0f,A,ev,1,1.0f,y,1); */
    my_sspmv(A,ev,y,n);

    sum = 0;
    for (i=0; i<n; i++) sum += y[i]*y[i];
    sum = sqrt(sum);

    d = 0;
    for (i=0; i<n; i++) {
      d += SQR(ev[i] - y[i]/sum);
      ev[i] = y[i]/sum;
    }
    fprintf(stderr," %5lu  %f\n",iter,d);
    if (d < 1.0e-6) break;
  }
  for (i=0; i<n; i++) ev[i] *= 100.0;
  fprintf(stderr," ecm done..\n");
}

VImage WriteOutput(VImage src,VImage map,int nslices,int nrows, int ncols, float *ev, size_t n)
{
  size_t i,b,r,c;

  VImage dest = VCreateImage(nslices,nrows,ncols,VFloatRepn);
  VFillImage(dest,VAllBands,0);
  VCopyImageAttrs (src, dest);

  for (i=0; i<n; i++) {
    b = VPixel(map,0,0,i,VInteger);
    r = VPixel(map,0,1,i,VInteger);
    c = VPixel(map,0,2,i,VInteger);
    /* if (ev[i] > 0) fprintf(stderr," %6d    %3d %3d %3d  %f\n",i,c,r,b,ev[i]); */
    VPixel(dest,b,r,c,VFloat) = ev[i];
  }
  return dest;
}

VAttrList VSpectralECM(VAttrList list,VImage mask,VDouble wavelength,VShort nlags,
		       VShort first,VShort length,VShort type)
{
  int b,r,c,j,k,kk,k0,n,p;
  size_t i;
  size_t m,nvox;
  int last;
  gsl_matrix *matcos=NULL,*matsin=NULL;
  float *A=NULL,*ev=NULL;
  float tr=-1,xtr=-1;
  double *in,*xsin=NULL,*xcos=NULL,*w=NULL;
  double nx,u,sumx,freq=0;
  gsl_matrix *tablecos=NULL,*tablesin=NULL;
  double pi = 2.0*acos(0);

  /* get image dimensions */
  int nrows=0,ncols=0,nt=0;
  int nslices = VAttrListNumImages(list);
  VImage *src = VAttrListGetImages(list,nslices);
  VImageDimensions(src,nslices,&nt,&nrows,&ncols);
  size_t ntimesteps = nt;

  if (VImageNRows(mask) != nrows)
    VError(" inconsistent image dims, mask has %d rows, image has %d rows",VImageNRows(mask),nrows);
  if (VImageNColumns(mask) != ncols)
    VError(" inconsistent image dims, mask has %d columns, image has %d columns",VImageNColumns(mask),ncols);
  if (VImageNBands(mask) != nslices)
    VError(" inconsistent image dims, mask has %d slices, image has %d slices",VImageNBands(mask),nslices);

  /* get time steps to include */
  if (first < 0) first = 0;
  if (length <= 0) length = ntimesteps;
  last = first + length -1;
  if (last >= ntimesteps) last = ntimesteps-1;

  n = last - first + 1;
  if (n < 2) VError(" not enough timesteps, nt= %d",n);
  nx = (double)n;

  /*
  ** get freq range
  */
  if (VGetAttr (VImageAttrList (src[0]), "repetition_time", NULL,
		VFloatRepn, (VPointer) & xtr) == VAttrFound) {
    tr = (xtr/1000.0);
  }
  if (tr <= 0) VError(" illegal tr: %f",tr);

  k0   = GetFFtIndex((int)n,tr,wavelength);
  freq = 1.0/wavelength;
  fprintf(stderr," TR= %g secs,  wavelength= %g secs,  freq= %.3f Hz,  k0= %d\n",
	  tr,wavelength,freq,k0);
  if (k0 < 1 || k0 >= n/2) VError(" illegal freq, index must be in [1,%d]\n",n/2);

  /* weights */
  if (nlags < 0)
    p = n/15;   /* default */
  else
    p = (int)nlags;
  if (p < 3) VError(" nlags too small, %d",p);
  w = Weights((int)n,k0,p);

  /* count number of voxels */
  nvox = 0;
  for (b=0; b<nslices; b++) {
    if (VImageNRows(src[b]) < 2) continue;
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {
	if (fabs(VGetPixel(mask,b,r,c)) < TINY) continue;
	nvox++;
      }
    }
  }
  fprintf(stderr," nvoxels: %ld,  n= %d,  nlags= %d\n",(long)nvox, (int)n, (int)p);

  /*
  ** voxel maps
  */
  VImage map = VCreateImage(1,4,nvox,VIntegerRepn);
  if (map == NULL) VError(" error allocating addr map");
  VFillImage(map,VAllBands,0);
  VPixel(map,0,3,0,VInteger) = nslices;
  VPixel(map,0,3,1,VInteger) = nrows;
  VPixel(map,0,3,2,VInteger) = ncols;

  double *xmap = (double *)VCalloc(nvox,sizeof(double));

  i = 0;
  for (b=0; b<nslices; b++) {
    if (VImageNRows(src[b]) < 2) continue;
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {
	if (fabs(VGetPixel(mask,b,r,c)) < TINY) continue;
	VPixel(map,0,0,i,VInteger) = b;
        VPixel(map,0,1,i,VInteger) = r;
        VPixel(map,0,2,i,VInteger) = c;
	i++;
      }
    }
  }

  /* alloc memory */
  in   = (double *)VCalloc(n,sizeof(double));
  xsin = (double *)VCalloc(n,sizeof(double));
  xcos = (double *)VCalloc(n,sizeof(double));
  tablecos = gsl_matrix_calloc((int)n,(int)n);
  if (tablecos == NULL) VError(" err alloc cos");
  tablesin = gsl_matrix_calloc((int)n,(int)n);
  if (tablesin == NULL) VError(" err alloc sin");
  Tables(tablecos,tablesin);

  /*
  ** precompute FFT
  */
  matcos = gsl_matrix_calloc(nvox,(int)n);
  if (!matcos) VError(" err alloc matcos");
  matsin = gsl_matrix_calloc(nvox,(int)n);
  if (!matsin) VError(" err alloc matsin");
  Tables(tablecos,tablesin);

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

    b = VPixel(map,0,0,i,VInteger);
    r = VPixel(map,0,1,i,VInteger);
    c = VPixel(map,0,2,i,VInteger);

    j = 0;
    for (k=first; k<=last; k++) {
      in[j++] = (double) VGetPixel(src[b],k-first,r,c);
    }
    double q = Normalize(in,n);
    if (fabs(q) < TINY) continue;

    FFT(in,xcos,xsin,tablecos,tablesin);
    nx = (double)n;

    sumx = 0;
    for (k=-p; k<=p; k++) {
      kk = k+k0;
      if (kk < 0) kk = -k;
      if (kk >= n) kk = n-k;
      u = 2.0*(xcos[kk]*xcos[kk] + xsin[kk]*xsin[kk])/nx;
      sumx += w[ABS(k)]*u;
    }
    xmap[i] = sumx;

    for (k=0; k<n; k++) {
      gsl_matrix_set(matcos,i,k,xcos[k]);
      gsl_matrix_set(matsin,i,k,xsin[k]);
    }
  }

  /*
  ** compute spectral coherence matrix
  */
  m = (nvox*(nvox+1))/2;
  fprintf(stderr," matrix computation, n= %ld...\n",(long)nvox);
  A = (float *) VCalloc(m,sizeof(float));
  if (!A) VError(" err allocating matrix");
  memset(A,0,m*sizeof(float));
   size_t progress=0;

#pragma omp parallel for shared(progress) schedule(guided) firstprivate(A,xmap)
  for (i=0; i<nvox; i++) {
    if (i%1000 == 0) fprintf(stderr," %d000  of  %ld\r",(int)(++progress),(long)nvox);

    size_t j=0,kij=0;
    int k=0,kk=0;
    double v=0,x=0,y=0,cs=0,qs=0,c1,s1,c2,s2,re,im;

    x = xmap[i];
    for (j=0; j<i; j++) {

      y = xmap[j];

      cs = qs = 0;
      for (k=-p; k<=p; k++) {

	kk = k+k0;
	if (kk  < 0) kk = -k;
	if (kk >= n) kk = n-k;

	c1 = gsl_matrix_get(matcos,i,kk);
	s1 = gsl_matrix_get(matsin,i,kk);

	c2 = gsl_matrix_get(matcos,j,kk);
	s2 = gsl_matrix_get(matsin,j,kk);

	re = 2.0*(c1*c2 + s1*s2)/nx;
	im = 2.0*(c1*s2 - s1*c2)/nx;

	cs += w[ABS(k)] * re;
	qs += w[ABS(k)] * im;
      }

      if (type == 0) {  /* spectral coherence */
	if (x*y > TINY) {
	  v = (cs*cs + qs*qs) / (x*y);
	  v = sqrt(v);
	}
	else v = -1;
      }
      else {            /* phase sync */
	v = atan2(qs,cs)/pi;
	/* v = 1.0 - ABS(v); */
	v = cos(v);
      }
      if (v < TINY) v = TINY;  /* make matrix irreducible */

      kij=j+i*(i+1)/2;
      A[kij] = v;
    }
  }

  /* free space */
  gsl_matrix_free(matcos);
  gsl_matrix_free(matsin);

  /*
  ** eigenvector centrality
  */
  ev = (float *) VCalloc(nvox,sizeof(float));
  EigenvectorCentrality(A,ev,nvox);
  /* DegreeCentrality(A,ev,nvox); */
  VImage dest0 = WriteOutput(src[0],map,nslices,nrows,ncols,ev,nvox);
  VSetAttr(VImageAttrList(dest0),"name",NULL,VStringRepn,"eigenvector_centrality");
  VAttrList out_list = VCreateAttrList();
  VAppendAttr(out_list,"image",NULL,VImageRepn,dest0);
  return out_list;
}

VDictEntry TYPDict[] = {
  { "freq", 0 },
  { "sync", 1 },
  { NULL }
};

int main (int argc,char *argv[])
{
  static VDouble  wavelength = 10;
  static VShort   nlags      = 10;
  static VShort   first      = 0;
  static VShort   length     = 0;
  static VShort   type       = 0;
  static VString  mask_filename   = "";
  static VShort   nproc  = 0;
  static VOptionDescRec  options[] = {
    {"wavelength", VDoubleRepn,1,(VPointer) &wavelength,VRequiredOpt,NULL,"wavelength in secs"},
    {"first",VShortRepn,1,(VPointer) &first,VOptionalOpt,NULL,"first timestep to use"},
    {"length",VShortRepn,1,(VPointer) &length,VOptionalOpt,NULL,"length of time series to use"},
    {"mask",VStringRepn,1,(VPointer) &mask_filename,VRequiredOpt,NULL,"mask"},
    {"nlags",VShortRepn,1,(VPointer) &nlags,VOptionalOpt,NULL,
     "num lags in cross-corr, use 0 to get default value"},
    {"type",VShortRepn,1,(VPointer) &type,VOptionalOpt,TYPDict,"frequency or phase coherence"},
    {"j",VShortRepn,1,(VPointer) &nproc,VOptionalOpt,NULL,"Number of processors to use, '0' to use all"}
  };
  FILE *in_file,*out_file;
  VString in_filename=NULL;
  char *prg = "vspectralecm";

  /* parse command line */
  VParseFilterCmdZ (VNumber (options),options,argc,argv,&in_file,&out_file,&in_filename);
  if (type > 2) VError(" illegal type");

  /* omp-stuff */
#ifdef _OPENMP
  int num_procs=omp_get_num_procs();
  if (nproc > 0 && nproc < num_procs) num_procs = nproc;
  fprintf(stderr,"using %d cores\n",(int)num_procs);
  omp_set_num_threads(num_procs);
#endif /* _OPENMP */

  /* read mask */
  VAttrList listm = VReadAttrList(mask_filename,0L,TRUE,FALSE);
  VImage mask = VReadImage(listm);
  if (mask == NULL) VError(" no ROI mask found");

  /* read functional data */
  VAttrList list = VReadAttrListZ(in_file,in_filename,0L,TRUE,FALSE);
  if (list == NULL) VError(" error reading input file %s",in_file);
  VAttrList geolist = VGetGeoInfo(list);

  /*
  ** process
  */
  VAttrList out_list = VSpectralECM(list,mask,wavelength,nlags,first,length,type);

  /* update geoinfo, 4D to 3D */
  if (geolist != NULL) {
    double *D = VGetGeoDim(geolist,NULL);
    D[0] = 3;  /* 3D */
    D[4] = 1;  /* just one timestep */
    VSetGeoDim(geolist,D);
    VSetGeoInfo(geolist,out_list);
  }

  VHistory(VNumber(options),options,prg,&list,&out_list);
  if (! VWriteFile (out_file, out_list)) exit (1);
  fprintf (stderr, "%s: done.\n", argv[0]);
  return 0;
}


text:
#pragma once
#include <mylibs/container.hpp> // pl::Vector
#include <initializer_list> // std::initializer_list
#include "FoodImpl.h" // fa::FoodWithAmount
#include "NoteImpl.h" // fa::NoteImpl
#include "PlanEntry.h" // fa::PlanEntryDynamicDelegator
#include "Date.h" // fa::Date
#include "GlobalHeader.h" // fa::OStream
#include <utility> // std::forward
#include <mylibs/utility.hpp> // pl::forEachArgument
#include <gsl.h> // gsl::not_null

namespace fa {
    namespace daily_plan_detail {
        template <class Type>
        using Container = pl::Vector<Type>;
    } // END of namespace daily_plan_detail

    class DailyPlanImpl final {
    public:
        using this_type = DailyPlanImpl;
        using DailyPlanFood = FoodWithAmount;
        using DailyPlanNote = NoteImpl;
        using Interface = PlanEntryInterface;
        using FoodUniqueOwner = pl::PolymorphicUniqueOwner<Interface, PlanEntryDynamicDelegator>;
        using NoteUniqueOwner = pl::PolymorphicUniqueOwner<Interface, PlanEntryDynamicDelegator>;
        using FoodCont = daily_plan_detail::Container<FoodUniqueOwner>;
        using NoteCont = daily_plan_detail::Container<NoteUniqueOwner>;

        //! put FoodWithAmounts and NoteImpls in here.
        template <class ...Args>
        explicit DailyPlanImpl(Args &&...args) : date_{ Date::today() } {
            pl::forEachArgument([this](auto &&e) {
                                    initCont(std::forward<decltype(e)>(e));
                                }, std::forward<Args>(args)...);
        }

        Date getDate() const;
        void setDate(Date const &);
        void setDate(Date &&);
        NutrientTable getAllNutrients() const;
        NutrientTable getMacroNutrients() const;
        NutrientTable getMicroNutrients() const;
        double getMacroPercentage() const;
        double getMicroPercentage() const;
        double getKcalPercentage() const;
        double getOverallPercentage() const;
        void addEntry(FoodWithAmount const &);
        void addEntry(FoodWithAmount &&);
        void addEntry(NoteImpl const &);
        void addEntry(NoteImpl &&);

        pl::View<FoodCont const> getFoods() const;
        pl::View<FoodCont> getFoods();
        pl::View<NoteCont const> getNotes() const;
        pl::View<NoteCont> getNotes();
        daily_plan_detail::Container<gsl::not_null<Interface const *>> getEntries() const;
        daily_plan_detail::Container<gsl::not_null<Interface *>> getEntries();
        daily_plan_detail::Container<gsl::not_null<DailyPlanFood *>> getFoodsDownCasted();
        daily_plan_detail::Container<gsl::not_null<DailyPlanFood const *>> getFoodsDownCasted() const;
        daily_plan_detail::Container<gsl::not_null<DailyPlanNote *>> getNotesDownCasted();
        daily_plan_detail::Container<gsl::not_null<DailyPlanNote const *>> getNotesDownCasted() const;

        friend OStream &operator<<(OStream &, this_type const &);

    private:
        void initCont(DailyPlanFood const &dpf);

        void initCont(DailyPlanFood &&dpf);

        void initCont(DailyPlanNote const &note);

        void initCont(DailyPlanNote &&note);

        NutrientTable getNutrients(FoodCont::const_iterator b,
                                   FoodCont::const_iterator e,
                                   NutrientTable accu) const;

        FoodCont foodCont_;
        NoteCont noteCont_;
        Date date_;
    }; // END of class DailyPlanImpl
} // END of namespace fa


text:
//This gets polynomial coeffs for each row/col of X using the Burg method.

//An opt mean0 is added to zero the mean of each row/col of X first.
//In this case, this is mean0 -> autocov_fft -> sig2ar_burg.

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

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

int sig2ar_burg_s (float *Y, float *V, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t P, const int mean0);
int sig2ar_burg_d (double *Y, double *V, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t P, const int mean0);

int sig2ar_burg_s (float *Y, float *V, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t P, const int mean0)
{
    if (dim>3) { fprintf(stderr,"error in sig2ar_burg_s: dim must be in [0 3]\n"); return 1; }

    const float o = 1.0f;
    const int N = (dim==0) ? R : C;
    int r, c, p, q;
    float m, g, b2, f2;
    float *b, *f, *oldf, *AStmp;

    //Checks
    if (R<1) { fprintf(stderr,"error in sig2ar_burg_s: nrows X must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in sig2ar_burg_s: ncols X must be positive\n"); return 1; }
    if (dim==0 && P>=R) { fprintf(stderr,"error in sig2ar_burg_s: P must be < nrows X for dim==0\n"); return 1; }
    if (dim==1 && P>=C) { fprintf(stderr,"error in sig2ar_burg_s: P must be < ncols X for dim==1\n"); return 1; }

    //Allocate
    if (!(b=(float *)malloc((size_t)(N)*sizeof(float)))) { fprintf(stderr,"error in sig2ar_burg_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(f=(float *)malloc((size_t)(N)*sizeof(float)))) { fprintf(stderr,"error in sig2ar_burg_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(oldf=(float *)malloc((size_t)(N)*sizeof(float)))) { fprintf(stderr,"error in sig2ar_burg_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(AStmp=(float *)malloc((size_t)(P)*sizeof(float)))) { fprintf(stderr,"error in sig2ar_burg_s: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0)
    {
        if (iscolmajor)
        {
            for (size_t c=0; c<C; c++)
            {
                cblas_scopy((int)R,&X[c*R],1,b,1);
                if (mean0)
                {
                    m = cblas_sdot((int)R,b,1,&o,0) / R;
                    cblas_saxpy((int)R,-m,&o,0,b,1);
                }
                cblas_scopy((int)R,b,1,f,1);
                V[c] = cblas_sdot((int)R,b,1,f,1) / R;
                for (p=1; p<=P; p++)
                {
                    b2 = cblas_sdot(R-p,b,1,b,1);
                    f2 = cblas_sdot(R-p,&f[p],1,&f[p],1);
                    Y[c*P+p-1] = g = (-2.0f/(b2+f2)) * cblas_sdot(R-p,b,1,&f[p],1);
                    V[c] *= fmaf(g,-g,1.0f);
                    if (p>1) { for (q=0; q<p; q++) { Y[c*P+q] += g * AStmp[p-q-2]; } }
                    if (p<P)
                    {
                        cblas_scopy(p,&Y[c*P],1,AStmp,1);
                        cblas_scopy(N-p,&f[p],1,&oldf[p],1);
                        cblas_saxpy(N-p,g,b,1,&f[p],1);
                        cblas_saxpy(N-p,g,&oldf[p],1,b,1);
                    }
                }
            }
        }
        else
        {
            for (size_t c=0; c<C; c++)
            {
                cblas_scopy((int)R,&X[c],(int)C,b,1);
                if (mean0)
                {
                    m = cblas_sdot((int)R,b,1,&o,0) / R;
                    cblas_saxpy((int)R,-m,&o,0,b,1);
                }
                cblas_scopy((int)R,b,1,f,1);
                V[c] = cblas_sdot((int)R,b,1,f,1) / R;
                for (p=1; p<=P; p++)
                {
                    b2 = cblas_sdot(R-p,b,1,b,1);
                    f2 = cblas_sdot(R-p,&f[p],1,&f[p],1);
                    Y[c+(p-1)*C] = g = (-2.0f/(b2+f2)) * cblas_sdot(R-p,b,1,&f[p],1);
                    V[c] *= fmaf(g,-g,1.0f);
                    if (p>1) { for (q=0; q<p; q++) { Y[c+q*C] += g * AStmp[p-q-2]; } }
                    if (p<P)
                    {
                        cblas_scopy(p,&Y[c],(int)C,AStmp,1);
                        cblas_scopy(N-p,&f[p],1,&oldf[p],1);
                        cblas_saxpy(N-p,g,b,1,&f[p],1);
                        cblas_saxpy(N-p,g,&oldf[p],1,b,1);
                    }
                }
            }
        }
        cblas_sscal(P*C,-1.0f,Y,1);
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (size_t r=0; r<R; r++)
            {
                cblas_scopy((int)C,&X[r],(int)R,b,1);
                if (mean0)
                {
                    m = cblas_sdot((int)C,b,1,&o,0) / C;
                    cblas_saxpy((int)C,-m,&o,0,b,1);
                }
                cblas_scopy((int)C,b,1,f,1);
                V[r] = cblas_sdot((int)C,b,1,f,1) / C;
                for (p=1; p<=P; p++)
                {
                    b2 = cblas_sdot(C-p,b,1,b,1);
                    f2 = cblas_sdot(C-p,&f[p],1,&f[p],1);
                    Y[r+(p-1)*R] = g = (-2.0f/(b2+f2)) * cblas_sdot(C-p,b,1,&f[p],1);
                    V[r] *= fmaf(g,-g,1.0f);
                    if (p>1) { for (q=0; q<p; q++) { Y[r+q*R] += g * AStmp[p-q-2]; } }
                    if (p<P)
                    {
                        cblas_scopy(p,&Y[r],(int)R,AStmp,1);
                        cblas_scopy(N-p,&f[p],1,&oldf[p],1);
                        cblas_saxpy(N-p,g,b,1,&f[p],1);
                        cblas_saxpy(N-p,g,&oldf[p],1,b,1);
                    }
                }
            }
        }
        else
        {
            for (size_t r=0; r<R; r++)
            {
                cblas_scopy((int)C,&X[r*C],1,b,1);
                if (mean0)
                {
                    m = cblas_sdot((int)C,b,1,&o,0) / C;
                    cblas_saxpy((int)C,-m,&o,0,b,1);
                }
                cblas_scopy((int)C,b,1,f,1);
                V[r] = cblas_sdot((int)C,b,1,f,1) / C;
                for (p=1; p<=P; p++)
                {
                    b2 = cblas_sdot(C-p,b,1,b,1);
                    f2 = cblas_sdot(C-p,&f[p],1,&f[p],1);
                    Y[r*P+p-1] = g = (-2.0f/(b2+f2)) * cblas_sdot(C-p,b,1,&f[p],1);
                    V[r] *= fmaf(g,-g,1.0f);
                    if (p>1) { for (q=0; q<p; q++) { Y[r*P+q] += g * AStmp[p-q-2]; } }
                    if (p<P)
                    {
                        cblas_scopy(p,&Y[r*P],1,AStmp,1);
                        cblas_scopy(N-p,&f[p],1,&oldf[p],1);
                        cblas_saxpy(N-p,g,b,1,&f[p],1);
                        cblas_saxpy(N-p,g,&oldf[p],1,b,1);
                    }
                }
            }
        }
        cblas_sscal(R*P,-1.0f,Y,1);
    }
    else
    {
        fprintf(stderr,"error in sig2ar_burg_s: dim must be 0 or 1.\n"); return 1;
    }

    //Exit
    free(b); free(f); free(oldf); free(AStmp);
    return 0;
}

int sig2ar_burg_d (double *Y, double *V, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t P, const int mean0)
{
    if (dim>3) { fprintf(stderr,"error in sig2ar_burg_d: dim must be in [0 3]\n"); return 1; }

    const double o = 1.0;
    const int N = (dim==0) ? R : C;
    int r, c, p, q;
    double m, g, b2, f2;
    double *b, *f, *oldf, *AStmp;

    //Checks
    if (R<1) { fprintf(stderr,"error in sig2ar_burg_d: nrows X must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in sig2ar_burg_d: ncols X must be positive\n"); return 1; }
    if (dim==0 && P>=R) { fprintf(stderr,"error in sig2ar_burg_d: P must be < nrows X for dim==0\n"); return 1; }
    if (dim==1 && P>=C) { fprintf(stderr,"error in sig2ar_burg_d: P must be < ncols X for dim==1\n"); return 1; }

    //Allocate
    if (!(b=(double *)malloc((size_t)(N)*sizeof(double)))) { fprintf(stderr,"error in sig2ar_burg_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(f=(double *)malloc((size_t)(N)*sizeof(double)))) { fprintf(stderr,"error in sig2ar_burg_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(oldf=(double *)malloc((size_t)(N)*sizeof(double)))) { fprintf(stderr,"error in sig2ar_burg_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(AStmp=(double *)malloc((size_t)(P)*sizeof(double)))) { fprintf(stderr,"error in sig2ar_burg_d: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0)
    {
        if (iscolmajor)
        {
            for (size_t c=0; c<C; c++)
            {
                cblas_dcopy((int)R,&X[c*R],1,b,1);
                if (mean0)
                {
                    m = cblas_ddot((int)R,b,1,&o,0) / R;
                    cblas_daxpy((int)R,-m,&o,0,b,1);
                }
                cblas_dcopy((int)R,b,1,f,1);
                V[c] = cblas_ddot((int)R,b,1,f,1) / R;
                for (p=1; p<=P; p++)
                {
                    b2 = cblas_ddot(R-p,b,1,b,1);
                    f2 = cblas_ddot(R-p,&f[p],1,&f[p],1);
                    Y[c*P+p-1] = g = (-2.0/(b2+f2)) * cblas_ddot(R-p,b,1,&f[p],1);
                    V[c] *= fma(g,-g,1.0);
                    if (p>1) { for (q=0; q<p; q++) { Y[c*P+q] += g * AStmp[p-q-2]; } }
                    if (p<P)
                    {
                        cblas_dcopy(p,&Y[c*P],1,AStmp,1);
                        cblas_dcopy(N-p,&f[p],1,&oldf[p],1);
                        cblas_daxpy(N-p,g,b,1,&f[p],1);
                        cblas_daxpy(N-p,g,&oldf[p],1,b,1);
                    }
                }
            }
        }
        else
        {
            for (size_t c=0; c<C; c++)
            {
                cblas_dcopy((int)R,&X[c],(int)C,b,1);
                if (mean0)
                {
                    m = cblas_ddot((int)R,b,1,&o,0) / R;
                    cblas_daxpy((int)R,-m,&o,0,b,1);
                }
                cblas_dcopy((int)R,b,1,f,1);
                V[c] = cblas_ddot((int)R,b,1,f,1) / R;
                for (p=1; p<=P; p++)
                {
                    b2 = cblas_ddot(R-p,b,1,b,1);
                    f2 = cblas_ddot(R-p,&f[p],1,&f[p],1);
                    Y[c+(p-1)*C] = g = (-2.0/(b2+f2)) * cblas_ddot(R-p,b,1,&f[p],1);
                    V[c] *= fma(g,-g,1.0);
                    if (p>1) { for (q=0; q<p; q++) { Y[c+q*C] += g * AStmp[p-q-2]; } }
                    if (p<P)
                    {
                        cblas_dcopy(p,&Y[c],(int)C,AStmp,1);
                        cblas_dcopy(N-p,&f[p],1,&oldf[p],1);
                        cblas_daxpy(N-p,g,b,1,&f[p],1);
                        cblas_daxpy(N-p,g,&oldf[p],1,b,1);
                    }
                }
            }
        }
        cblas_dscal(P*C,-1.0,Y,1);
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (size_t r=0; r<R; r++)
            {
                cblas_dcopy((int)C,&X[r],(int)R,b,1);
                if (mean0)
                {
                    m = cblas_ddot((int)C,b,1,&o,0) / C;
                    cblas_daxpy((int)C,-m,&o,0,b,1);
                }
                cblas_dcopy((int)C,b,1,f,1);
                V[r] = cblas_ddot((int)C,b,1,f,1) / C;
                for (p=1; p<=P; p++)
                {
                    b2 = cblas_ddot(C-p,b,1,b,1);
                    f2 = cblas_ddot(C-p,&f[p],1,&f[p],1);
                    Y[r+(p-1)*R] = g = (-2.0/(b2+f2)) * cblas_ddot(C-p,b,1,&f[p],1);
                    V[r] *= fma(g,-g,1.0);
                    if (p>1) { for (q=0; q<p; q++) { Y[r+q*R] += g * AStmp[p-q-2]; } }
                    if (p<P)
                    {
                        cblas_dcopy(p,&Y[r],(int)R,AStmp,1);
                        cblas_dcopy(N-p,&f[p],1,&oldf[p],1);
                        cblas_daxpy(N-p,g,b,1,&f[p],1);
                        cblas_daxpy(N-p,g,&oldf[p],1,b,1);
                    }
                }
            }
        }
        else
        {
            for (size_t r=0; r<R; r++)
            {
                cblas_dcopy((int)C,&X[r*C],1,b,1);
                if (mean0)
                {
                    m = cblas_ddot((int)C,b,1,&o,0) / C;
                    cblas_daxpy((int)C,-m,&o,0,b,1);
                }
                cblas_dcopy((int)C,b,1,f,1);
                V[r] = cblas_ddot((int)C,b,1,f,1) / C;
                for (p=1; p<=P; p++)
                {
                    b2 = cblas_ddot(C-p,b,1,b,1);
                    f2 = cblas_ddot(C-p,&f[p],1,&f[p],1);
                    Y[r*P+p-1] = g = (-2.0/(b2+f2)) * cblas_ddot(C-p,b,1,&f[p],1);
                    V[r] *= fma(g,-g,1.0);
                    if (p>1) { for (q=0; q<p; q++) { Y[r*P+q] += g * AStmp[p-q-2]; } }
                    if (p<P)
                    {
                        cblas_dcopy(p,&Y[r*P],1,AStmp,1);
                        cblas_dcopy(N-p,&f[p],1,&oldf[p],1);
                        cblas_daxpy(N-p,g,b,1,&f[p],1);
                        cblas_daxpy(N-p,g,&oldf[p],1,b,1);
                    }
                }
            }
        }
        cblas_dscal(R*P,-1.0,Y,1);
    }
    else
    {
        fprintf(stderr,"error in sig2ar_burg_d: dim must be 0 or 1.\n"); return 1;
    }

    //Exit
    free(b); free(f); free(oldf); free(AStmp);
    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
#include <cblas.h>
#include <stdio.h>
#include <assert.h>
/*

 * y = alpha*A*x + beta*y,
 * where alpha and beta are scalars, x and y are vectors and A is an
 * m by n matrix.
 */

void cblas_dgemv(const enum CBLAS_ORDER order,
                 const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
                 const double alpha, const double *A, const int lda,
                 const double *X, const int incX, const double beta,
		 double *Y, const int incY)

{
  assert(order == CblasColMajor);
  assert(TransA == CblasNoTrans);
  int i;
  for (i = 0; i < N; ++i)
    {
      Y[i*incY] += beta* Y[i*incY] + alpha* cblas_ddot(N, A[i] , lda , X , incX);
    }

}

max_stars_repo_name: ethiery/HPCassignments

max_issues_repo_name: ethiery/HPCassignments

max_forks_repo_name: ethiery/HPCassignments


text:
#ifndef OPENMC_TALLIES_FILTER_COLLISIONS_H
#define OPENMC_TALLIES_FILTER_COLLISIONS_H

#include <vector>
#include <unordered_map>
#include <gsl/gsl>

#include "openmc/tallies/filter.h"

namespace openmc {

//==============================================================================
//! Bins the incident neutron energy.
//==============================================================================

class CollisionFilter : public Filter {
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors

  ~CollisionFilter() = default;

  //----------------------------------------------------------------------------
  // Methods

  std::string type() const override { return "collision"; }

  void from_xml(pugi::xml_node node) override;

  void get_all_bins(const Particle& p, TallyEstimator estimator,
    FilterMatch& match) const override;

  void to_statepoint(hid_t filter_group) const override;

  std::string text_label(int bin) const override;

  //----------------------------------------------------------------------------
  // Accessors

  const std::vector<int>& bins() const { return bins_; }
  void set_bins(gsl::span<const int> bins);

protected:
  //----------------------------------------------------------------------------
  // Data members

  std::vector<int> bins_;

  std::unordered_map<int,int> map_;

};

} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_COLLISIONS_H

max_stars_repo_name: HunterBelanger/openmc

max_issues_repo_name: HunterBelanger/openmc

max_forks_repo_name: HunterBelanger/openmc


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2018 Couchbase, Inc.
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */

#pragma once

#include <nlohmann/json_fwd.hpp>
#include <gsl/gsl>
#include <list>
#include <string>

class Module;
class Event;

/**
 * Is this build for enterprise edition?
 *
 * @return true when building EE, false for CE
 */
bool is_enterprise_edition();

/**
 * In order to allow making unit tests we want to be able to mock the
 * enterprise edition settings dynamically
 */
void set_enterprise_edition(bool enable);

/**
 * Load the requested file and parse it as JSON
 *
 * @param fname the name of the file
 * @return the json representation of the file
 * @throws std::system_error if we fail to read the file
 *         std::runtime_error if we fail to parse the content of the file
 */
nlohmann::json load_file(const std::string& fname);

/**
 * Iterate over the module descriptor json and populate each entry
 * in the modules array into the provided modules list.
 *
 * @param ptr The JSON representation of the module description. See
 *            ../README.md for a description of the syntax
 * @param modules Where to store the list of all of the entries found
 * @param srcroot The source root to prepend to all of the paths in the spec
 * @param objroot The object root to prepend to all of the paths in the spec
 * @throws std::invalid_argument if the provided JSON is of an unexpected format
 */
void parse_module_descriptors(const nlohmann::json&,
                              std::list<std::unique_ptr<Module>>& modules,
                              const std::string& srcroot,
                              const std::string& objroot);

/**
 * Build the master event file
 *
 * @param modules The modules to include
 * @param output_file Where to store the result
 * @throws std::system_error if we fail to write the file
 */
void create_master_file(const std::list<std::unique_ptr<Module>>& modules,
                        const std::string& output_file);


text:
//  This matrix class is a C++ wrapper for the GNU Scientific Library
//  Copyright (C)  ULP-IPB Strasbourg

//  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., 675 Mass Ave, Cambridge, MA 02139, USA.

#ifndef _vector_double_h
#define _vector_double_h

#ifdef __HP_aCC
#include <iostream.h>
#else
#include <iostream>
#endif

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gslwrap/vector_double.h>

//#define NDEBUG 0

#include <assert.h>
namespace gsl
{

#ifndef __HP_aCC
	using std::ostream;
//using std::string;
//using std::runtime_error;
#endif

class vector_view;

class vector
{
protected:
	gsl_vector *gsldata;
	void free(){if(gsldata) gsl_vector_free(gsldata);gsldata=NULL;}
	void alloc(size_t n) {gsldata=gsl_vector_alloc(n);}
	void calloc(size_t n){gsldata=gsl_vector_calloc(n);}
public:
	typedef double value_type;
	vector() : gsldata(NULL) {;}
	vector( const vector &other ):gsldata(NULL) {copy(other);}
	template<class oclass>
	vector( const oclass &other ):gsldata(NULL) {copy(other);}
	~vector(){free();}
	vector(const size_t& n,bool clear=true)
	{
		if(clear){this->calloc(n);}
		else     {this->alloc(n);}
	}
	vector(const int& n,bool clear=true)
	{
		if(clear){this->calloc(n);}
		else     {this->alloc(n);}
	}

	void resize(size_t n);

	template <class oclass>
		void copy(const oclass &other)
		{
			if ( static_cast<const void *>( this ) == static_cast<const void *>( &other ) )
				return;

			if (!other.is_set())
			{
				gsldata=NULL;
				return;
			}
			resize(other.size());
			for (size_t i=0;i<size();i++)
			{
				gsl_vector_set(gsldata, i, (double)other[i]);
			}
		}
	void copy(const vector& other);
	bool is_set() const{if (gsldata) return true; else return false;}
//	void clone(vector& other);

//	size_t size() const {if (!gsldata) {cout << "vector::size vector not initialized" << endl; exit(-1);}return gsldata->size;}
	size_t size() const {assert (gsldata); return gsldata->size;}

	/** for interfacing with gsl c */
/*  	gsl_vector       *gslobj()       {if (!gsldata){cout << "vector::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return gsldata;} */
/*  	const gsl_vector *gslobj() const {if (!gsldata){cout << "vector::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return gsldata;} */
	gsl_vector       *gslobj()       {assert(gsldata);return gsldata;}
	const gsl_vector *gslobj() const {assert(gsldata);return gsldata;}

	static vector_view create_vector_view( const gsl_vector_view &other );

// ********Accessing vector elements

//  Unlike FORTRAN compilers, C compilers do not usually provide support for range checking of vectors and matrices (2). However, the functions gsl_vector_get and gsl_vector_set can perform range checking for you and report an error if you attempt to access elements outside the allowed range.

//  The functions for accessing the elements of a vector or matrix are defined in `gsl_vector.h' and declared extern inline to eliminate function-call overhead. If necessary you can turn off range checking completely without modifying any source files by recompiling your program with the preprocessor definition GSL_RANGE_CHECK_OFF. Provided your compiler supports inline functions the effect of turning off range checking is to replace calls to gsl_vector_get(v,i) by v->data[i*v->stride] and and calls to gsl_vector_set(v,i,x) by v->data[i*v->stride]=x. Thus there should be no performance penalty for using the range checking functions when range checking is turned off.

//      This function returns the i-th element of a vector v. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked and 0 is returned.
	double get(size_t i) const {return gsl_vector_get(gsldata,i);}

//      This function sets the value of the i-th element of a vector v to x. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked.
	void  set(size_t i,double x){gsl_vector_set(gsldata,i,x);}

//      These functions return a pointer to the i-th element of a vector v. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked
	double       &operator[](size_t i)       { return *gsl_vector_ptr(gsldata,i);}
	const double &operator[](size_t i) const { return *gsl_vector_ptr(gsldata,i);}

	double       &operator()(size_t i)       { return *gsl_vector_ptr(gsldata,i);}
	const double &operator()(size_t i) const { return *gsl_vector_ptr(gsldata,i);}

//  ***** Initializing vector elements

//      This function sets all the elements of the vector v to the value x.
	void set_all(double x){gsl_vector_set_all (gsldata,x);}
//      This function sets all the elements of the vector v to zero.
	void set_zero(){gsl_vector_set_zero (gsldata);}

//      This function makes a basis vector by setting all the elements of the vector v to zero except for the i-th element which is set to one.
	int set_basis (size_t i) {return gsl_vector_set_basis (gsldata,i);}

//  **** Reading and writing vectors

//  The library provides functions for reading and writing vectors to a file as binary data or formatted text.

//      This function writes the elements of the vector v to the stream stream in binary format. The return value is 0 for success and GSL_EFAILED if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
	int fwrite (FILE * stream) const {return gsl_vector_fwrite (stream, gsldata);}

//      This function reads into the vector v from the open stream stream in binary format. The vector v must be preallocated with the correct length since the function uses the size of v to determine how many bytes to read. The return value is 0 for success and GSL_EFAILED if there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture.
	int fread (FILE * stream) {return gsl_vector_fread (stream, gsldata);}

	void load( const char *filename );
	///
	void save( const char *filename ) const;

//      This function writes the elements of the vector v line-by-line to the stream stream using the format specifier format, which should be one of the %g, %e or %f formats for floating point numbers and %d for integers. The function returns 0 for success and GSL_EFAILED if there was a problem writing to the file.
	int fprintf (FILE * stream, const char * format) const {return gsl_vector_fprintf (stream, gsldata,format) ;}

//      This function reads formatted data from the stream stream into the vector v. The vector v must be preallocated with the correct length since the function uses the size of v to determine how many numbers to read. The function returns 0 for success and GSL_EFAILED if there was a problem reading from the file.
	int fscanf (FILE * stream)  {return gsl_vector_fscanf (stream, gsldata); }

//  ******* Vector views

//  In addition to creating vectors from slices of blocks it is also possible to slice vectors and create vector views. For example, a subvector of another vector can be described with a view, or two views can be made which provide access to the even and odd elements of a vector.

//  A vector view is a temporary object, stored on the stack, which can be used to operate on a subset of vector elements. Vector views can be defined for both constant and non-constant vectors, using separate types that preserve constness. A vector view has the type gsl_vector_view and a constant vector view has the type gsl_vector_const_view. In both cases the elements of the view can be accessed as a gsl_vector using the vector component of the view object. A pointer to a vector of type gsl_vector * or const gsl_vector * can be obtained by taking the address of this component with the & operator.

//      These functions return a vector view of a subvector of another vector v. The start of the new vector is offset by offset elements from the start of the original
//      vector. The new vector has n elements. Mathematically, the i-th element of the new vector v' is given by,

//      v'(i) = v->data[(offset + i)*v->stride]

//      where the index i runs from 0 to n-1.

//      The data pointer of the returned vector struct is set to null if the combined parameters (offset,n) overrun the end of the original vector.

//      The new vector is only a view of the block underlying the original vector, v. The block containing the elements of v is not owned by the new vector. When the
//      new vector goes out of scope the original vector v and its block will continue to exist. The original memory can only be deallocated by freeing the original vector.
//      Of course, the original vector should not be deallocated while the new vector is still in use.

//      The function gsl_vector_const_subvector is equivalent to gsl_vector_subvector but can be used for vectors which are declared const.
	vector_view subvector (size_t offset, size_t n);
	const vector_view subvector (size_t offset, size_t n) const;
//	vector_const_view subvector (size_t offset, size_t n) const;

//  	class view
//  	{
//  		gsl_vector_view *gsldata;
//  	public:
//  		view();
//  	};
//  	view subvector(size_t offset, size_t n)
//  	{
//  		return view(gsl_vector_subvector(gsldata,offset,n);
//  	}
//  	const view subvector(size_t offset, size_t n) const
//  	{
//  		return view(gsl_vector_const_subvector(gsldata,offset,n);
//  	}

//  Function: gsl_vector gsl_vector_subvector_with_stride (gsl_vector *v, size_t offset, size_t stride, size_t n)
//  Function: gsl_vector_const_view gsl_vector_const_subvector_with_stride (const gsl_vector * v, size_t offset, size_t stride, size_t n)
//      These functions return a vector view of a subvector of another vector v with an additional stride argument. The subvector is formed in the same way as for
//      gsl_vector_subvector but the new vector has n elements with a step-size of stride from one element to the next in the original vector. Mathematically,
//      the i-th element of the new vector v' is given by,

//      v'(i) = v->data[(offset + i*stride)*v->stride]

//      where the index i runs from 0 to n-1.

//      Note that subvector views give direct access to the underlying elements of the original vector. For example, the following code will zero the even elements of the
//      vector v of length n, while leaving the odd elements untouched,

//      gsl_vector_view v_even = gsl_vector_subvector_with_stride (v, 0, 2, n/2);
//      gsl_vector_set_zero (&v_even.vector);

//      A vector view can be passed to any subroutine which takes a vector argument just as a directly allocated vector would be, using &view.vector. For example, the
//      following code computes the norm of odd elements of v using the BLAS routine DNRM2,

//      gsl_vector_view v_odd = gsl_vector_subvector_with_stride (v, 1, 2, n/2);
//      double r = gsl_blas_dnrm2 (&v_odd.vector);

//      The function gsl_vector_const_subvector_with_stride is equivalent to gsl_vector_subvector_with_stride but can be used for
//      vectors which are declared const.

//  Function: gsl_vector_view gsl_vector_complex_real (gsl_vector_complex *v)
//  Function: gsl_vector_const_view gsl_vector_complex_const_real (const gsl_vector_complex *v)
//      These functions return a vector view of the real parts of the complex vector v.

//      The function gsl_vector_complex_const_real is equivalent to gsl_vector_complex_real but can be used for vectors which are declared
//      const.

//  Function: gsl_vector_view gsl_vector_complex_imag (gsl_vector_complex *v)
//  Function: gsl_vector_const_view gsl_vector_complex_const_imag (const gsl_vector_complex *v)
//      These functions return a vector view of the imaginary parts of the complex vector v.

//      The function gsl_vector_complex_const_imag is equivalent to gsl_vector_complex_imag but can be used for vectors which are declared
//      const.

//  Function: gsl_vector_view gsl_vector_view_array (double *base, size_t n)
//  Function: gsl_vector_const_view gsl_vector_const_view_array (const double *base, size_t n)
//      These functions return a vector view of an array. The start of the new vector is given by base and has n elements. Mathematically, the i-th element of the new
//      vector v' is given by,

//      v'(i) = base[i]

//      where the index i runs from 0 to n-1.

//      The array containing the elements of v is not owned by the new vector view. When the view goes out of scope the original array will continue to exist. The
//      original memory can only be deallocated by freeing the original pointer base. Of course, the original array should not be deallocated while the view is still in use.

//      The function gsl_vector_const_view_array is equivalent to gsl_vector_view_array but can be used for vectors which are declared const.

//  Function: gsl_vector_view gsl_vector_view_array_with_stride (double * base, size_t stride, size_t n)
//  Function: gsl_vector_const_view gsl_vector_const_view_array_with_stride (const double * base, size_t stride, size_t n)
//      These functions return a vector view of an array base with an additional stride argument. The subvector is formed in the same way as for
//      gsl_vector_view_array but the new vector has n elements with a step-size of stride from one element to the next in the original array. Mathematically,
//      the i-th element of the new vector v' is given by,

//      v'(i) = base[i*stride]

//      where the index i runs from 0 to n-1.

//      Note that the view gives direct access to the underlying elements of the original array. A vector view can be passed to any subroutine which takes a vector
//      argument just as a directly allocated vector would be, using &view.vector.

//      The function gsl_vector_const_view_array_with_stride is equivalent to gsl_vector_view_array_with_stride but can be used for
//      arrays which are declared const.

//  ************* Copying vectors

//  Common operations on vectors such as addition and multiplication are available in the BLAS part of the library (see section BLAS Support). However, it is useful to have a small number of utility functions which do not require the full BLAS code. The following functions fall into this category.

//      This function copies the elements of the vector src into the vector dest.
	vector& operator=(const vector& other){copy(other);return (*this);}

//  Function: int gsl_vector_swap (gsl_vector * v, gsl_vector * w)
//      This function exchanges the elements of the vectors v and w by copying. The two vectors must have the same length.

//  ***** Exchanging elements

//  The following function can be used to exchange, or permute, the elements of a vector.

//  Function: int gsl_vector_swap_elements (gsl_vector * v, size_t i, size_t j)
//      This function exchanges the i-th and j-th elements of the vector v in-place.
	int swap_elements (size_t i, size_t j) {return gsl_vector_swap_elements (gsldata, i,j);}

//  Function: int gsl_vector_reverse (gsl_vector * v)
//      This function reverses the order of the elements of the vector v.
	int reverse () {return  gsl_vector_reverse (gsldata) ;}

// ******* Vector operations

//  The following operations are only defined for real vectors.

//      This function adds the elements of vector b to the elements of vector a, a'_i = a_i + b_i. The two vectors must have the same length.
	int operator+=(const vector &other) {return gsl_vector_add (gsldata, other.gsldata);}

//      This function subtracts the elements of vector b from the elements of vector a, a'_i = a_i - b_i. The two vectors must have the same length.
	int operator-=(const vector &other) {return gsl_vector_sub (gsldata, other.gsldata);}

//  Function: int gsl_vector_mul (gsl_vector * a, const gsl_vector * b)
//      This function multiplies the elements of vector a by the elements of vector b, a'_i = a_i * b_i. The two vectors must have the same length.
	int operator*=(const vector &other) {return gsl_vector_mul (gsldata, other.gsldata);}

//      This function divides the elements of vector a by the elements of vector b, a'_i = a_i / b_i. The two vectors must have the same length.
	int operator/=(const vector &other) {return gsl_vector_div (gsldata, other.gsldata);}

//      This function multiplies the elements of vector a by the constant factor x, a'_i = x a_i.
	int operator*=(double x) {return gsl_vector_scale (gsldata, x);}

//  Function: int gsl_vector_add_constant (gsl_vector * a, const double x)
//      This function adds the constant value x to the elements of the vector a, a'_i = a_i + x.
	int operator+=(double x) {return gsl_vector_add_constant (gsldata,x);}

//      This function multiplies the elements of vector a by the constant factor x, a'_i = x a_i.
	int operator/=(double x) {return gsl_vector_scale (gsldata, 1/x);}

// bool operators:
	bool operator==(const vector& other) const;
	bool operator!=(const vector& other) const { return (!((*this)==other));}

// stream output:
//	friend ostream& operator<< ( ostream& os, const vector& vect );
	/** returns sum of all the vector elements. */
    double sum() const;
	// returns sqrt(v.t*v);
    double norm2() const;

// **** Finding maximum and minimum elements of vectors

//      This function returns the maximum value in the vector v.
    double max() const{return gsl_vector_max (gsldata) ;}

//  Function: double gsl_vector_min (const gsl_vector * v)
//      This function returns the minimum value in the vector v.
    double min() const{return gsl_vector_min (gsldata) ;}

//  Function: void gsl_vector_minmax (const gsl_vector * v, double * min_out, double * max_out)
//      This function returns the minimum and maximum values in the vector v, storing them in min_out and max_out.

//      This function returns the index of the maximum value in the vector v. When there are several equal maximum elements then the lowest index is returned.
	size_t max_index(){return gsl_vector_max_index (gsldata);}

//  Function: size_t gsl_vector_min_index (const gsl_vector * v)
//      This function returns the index of the minimum value in the vector v. When there are several equal minimum elements then the lowest index is returned.
	size_t min_index(){return gsl_vector_min_index (gsldata);}

//  Function: void gsl_vector_minmax_index (const gsl_vector * v, size_t * imin, size_t * imax)
//      This function returns the indices of the minimum and maximum values in the vector v, storing them in imin and imax. When there are several equal minimum
//      or maximum elements then the lowest indices are returned.

//  Vector properties

//  Function: int gsl_vector_isnull (const gsl_vector * v)
//      This function returns 1 if all the elements of the vector v are zero, and 0 otherwise. };
	bool isnull(){return gsl_vector_isnull (gsldata);}
};

// When you add create a view it will stick to its with the view until you call change_view
// ex:
// matrix_float m(5,5);
// vector_float v(5);
// // ...
// m.column(3) = v; //the 3rd column of the matrix m will equal v.
class vector_view : public vector
{
 public:
	vector_view(const vector&     other) :vector(){init(other);}
	vector_view(const vector_view& other):vector(){init(other);}
	vector_view(const gsl_vector& gsl_other) : vector() {init_with_gsl_vector(gsl_other);}

	void init(const vector& other);
	void init_with_gsl_vector(const gsl_vector& gsl_other);
	void change_view(const vector& other){init(other);}
 private:
};

ostream& operator<< ( ostream& os, const vector & vect );

// vector_type<>::type is a template interface to vector_?
// it is usefull for in templated situations for getting the correct vector type
#define tmp_type_is
#ifdef tmp_type_is
typedef vector vector_double;
template<class T>
struct vector_type  {typedef vector_double   type;};

template<class T>
struct value_type  {typedef double   type;};

#else
template<> struct vector_type<double> {typedef vector type;};
#endif
#undef tmp_type_is

}
#endif// _vector_double_h

max_stars_repo_name: boomsbloom/dtm-fmri

max_issues_repo_name: boomsbloom/dtm-fmri

max_forks_repo_name: boomsbloom/dtm-fmri


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT license.

#pragma once

#if SEAL_COMPILER == SEAL_COMPILER_GCC

// We require GCC >= 6
#if (__GNUC__ < 6) || not defined(__cplusplus)
#pragma GCC error "SEAL requires __GNUC__ >= 6"
#endif

// Read in config.h
#include "seal/util/config.h"

#if (__GNUC__ == 6) && defined(SEAL_USE_IF_CONSTEXPR)
#pragma GCC error "g++-6 cannot compile Microsoft SEAL as C++17; set CMake build option `SEAL_USE_CXX17' to OFF"
#endif

// Are we using MSGSL?
#ifdef SEAL_USE_MSGSL
#include <gsl/gsl>
#endif

// Are intrinsics enabled?
#ifdef SEAL_USE_INTRIN
#include <x86intrin.h>

#ifdef SEAL_USE___BUILTIN_CLZLL
#define SEAL_MSB_INDEX_UINT64(result, value) {                                      \
    *result = 63UL - static_cast<unsigned long>(__builtin_clzll(value));            \
}
#endif

#ifdef SEAL_USE___INT128
#define SEAL_MULTIPLY_UINT64_HW64(operand1, operand2, hw64) {                       \
    *hw64 = static_cast<unsigned long long>(                                        \
            ((static_cast<unsigned __int128>(operand1)                              \
            * static_cast<unsigned __int128>(operand2)) >> 64));                    \
}

#define SEAL_MULTIPLY_UINT64(operand1, operand2, result128) {                       \
    unsigned __int128 product = static_cast<unsigned __int128>(operand1) * operand2;\
    result128[0] = static_cast<unsigned long long>(product);                        \
    result128[1] = static_cast<unsigned long long>(product >> 64);                  \
}

#define SEAL_DIVIDE_UINT128_UINT64(numerator, denominator, result) {                \
    unsigned __int128 n, q;                                                         \
    n = (static_cast<unsigned __int128>(numerator[1]) << 64) |                      \
        (static_cast<unsigned __int128>(numerator[0]));                             \
    q = n / denominator;                                                            \
    n -= q * denominator;                                                           \
    numerator[0] = static_cast<std::uint64_t>(n);                                   \
    numerator[1] = static_cast<std::uint64_t>(n >> 64);                             \
    quotient[0] = static_cast<std::uint64_t>(q);                                    \
    quotient[1] = static_cast<std::uint64_t>(q >> 64);                              \
}
#endif

#ifdef SEAL_USE__ADDCARRY_U64
#define SEAL_ADD_CARRY_UINT64(operand1, operand2, carry, result) _addcarry_u64(     \
    carry, operand1, operand2, result)
#endif

#ifdef SEAL_USE__SUBBORROW_U64
#if ((__GNUC__ == 7) && (__GNUC_MINOR__ >= 2)) || (__GNUC__  >= 8)
// The inverted arguments problem was fixed in GCC-7.2
// (https://patchwork.ozlabs.org/patch/784309/)
#define SEAL_SUB_BORROW_UINT64(operand1, operand2, borrow, result) _subborrow_u64(  \
    borrow, operand1, operand2, result)
#else
// Warning: Note the inverted order of operand1 and operand2
#define SEAL_SUB_BORROW_UINT64(operand1, operand2, borrow, result) _subborrow_u64(  \
    borrow, operand2, operand1, result)
#endif //(__GNUC__ == 7) && (__GNUC_MINOR__ >= 2)
#endif

#endif //SEAL_USE_INTRIN

#endif

max_issues_repo_name: KoizumiRuby/SEAL

max_forks_repo_name: KoizumiRuby/SEAL


text:
// Copyright 2014 The Chromium Authors. All rights reserved.
// Inspired by Chromium video capture interface
// Simplified and stripped from internal base code

#ifndef S3D_VIDEO_CAPTURE_VIDEO_CAPTURE_DEVICE_DECKLINK_H
#define S3D_VIDEO_CAPTURE_VIDEO_CAPTURE_DEVICE_DECKLINK_H

#include <s3d/video/capture/video_capture_device.h>
#include <s3d/video/capture/video_capture_device_factory.h>

#include <gsl/gsl>

#include <chrono>

namespace s3d {

class DecklinkCaptureDelegate;

class VideoCaptureDeviceDecklink : public VideoCaptureDevice {
 public:
  explicit VideoCaptureDeviceDecklink(const VideoCaptureDeviceDescriptor& deviceDescriptor);

  gsl::owner<VideoCaptureDevice*> clone() const override;

  ~VideoCaptureDeviceDecklink() override;

  void AllocateAndStart(const VideoCaptureFormat& format,
                        VideoCaptureDevice::Client* client) override;

  void StopAndDeAllocate() override;

  // called from delegate
  void OnIncomingCapturedData(const VideoCaptureDevice::Client::Images& images,
                              const VideoCaptureFormat& frameFormat,
                              std::chrono::microseconds /*timestamp*/);

  VideoCaptureFormat DefaultFormat() override;

 private:
  VideoCaptureDevice::Client* client_{};
  std::unique_ptr<DecklinkCaptureDelegate> captureDelegate_;
};

}  // namespace s3d

#endif  // S3D_VIDEO_CAPTURE_VIDEO_CAPTURE_DEVICE_DECKLINK_H


text:
/* linalg/test_ql.c
 *
 * Copyright (C) 2019 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>

static int
test_QL_decomp_eps(const gsl_matrix * m, const double eps, const char * desc)
{
  int s = 0;
  const size_t M = m->size1;
  const size_t N = m->size2;
  size_t i, j;

  gsl_matrix * QL = gsl_matrix_alloc(M, N);
  gsl_vector * tau = gsl_vector_alloc(N);
  gsl_matrix * A  = gsl_matrix_alloc(M, N);
  gsl_matrix * L  = gsl_matrix_alloc(M, N);
  gsl_matrix * Q  = gsl_matrix_alloc(M, M);

  gsl_matrix_memcpy(QL, m);

  s += gsl_linalg_QL_decomp(QL, tau);
  s += gsl_linalg_QL_unpack(QL, tau, Q, L);

  /* compute A = Q L */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, Q, L, 0.0, A);

  for (i = 0; i < M; i++)
    {
      for (j = 0; j < N; j++)
        {
          double aij = gsl_matrix_get(A, i, j);
          double mij = gsl_matrix_get(m, i, j);

          gsl_test_rel(aij, mij, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i,j, aij, mij);
        }
    }

  gsl_matrix_free(QL);
  gsl_vector_free(tau);
  gsl_matrix_free(A);
  gsl_matrix_free(Q);
  gsl_matrix_free(L);

  return s;
}

static int
test_QL_decomp(gsl_rng * r)
{
  int s = 0;
  size_t M, N;

  for (M = 1; M <= 30; ++M)
    {
      for (N = 1; N <= M; ++N)
        {
          gsl_matrix * A = gsl_matrix_alloc(M, N);

          create_random_matrix(A, r);
          s += test_QL_decomp_eps(A, 1.0e5 * GSL_MAX(M, N) * GSL_DBL_EPSILON, "QL_decomp random");

          gsl_matrix_free(A);
        }
    }

  s += test_QL_decomp_eps(m53,   1.0e2 * GSL_DBL_EPSILON, "QL_decomp m(5,3)");
  s += test_QL_decomp_eps(hilb2, 1.0e2 * GSL_DBL_EPSILON, "QL_decomp hilbert(2)");
  s += test_QL_decomp_eps(hilb3, 1.0e2 * GSL_DBL_EPSILON, "QL_decomp hilbert(3)");
  s += test_QL_decomp_eps(hilb4, 1.0e2 * GSL_DBL_EPSILON, "QL_decomp hilbert(4)");
  s += test_QL_decomp_eps(hilb12, 1.0e2 * GSL_DBL_EPSILON, "QL_decomp hilbert(12)");
  s += test_QL_decomp_eps(vander2, 1.0e1 * GSL_DBL_EPSILON, "QL_decomp vander(2)");
  s += test_QL_decomp_eps(vander3, 1.0e1 * GSL_DBL_EPSILON, "QL_decomp vander(3)");
  s += test_QL_decomp_eps(vander4, 1.0e2 * GSL_DBL_EPSILON, "QL_decomp vander(4)");

  return s;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
//cmc_stateSpace.c --
// this program creates and outputs a coalescent markov chain's complete state space
// given a particular model

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include "AFS.h"

int n1, n2;

void getParameters(int argc,const char **argv);
void usage();

int main(int argc, const char * argv[]){
	int i;
	afsObject *A;
	afsStateSpace *stateSpace;

	getParameters(argc,argv);

	stateSpace = afsStateSpaceNew();

	A = afsObjectNew(n1,n2);
 	afsObjectInit(A,n1,n2);
//
// 	printf("%s\n",A->matString );
// //	exit(1);
	stateSpace->states[0] = A;
	stateSpace->nstates=1;
//	stateSpace->tmpStorage[A->nalleles][n1][0] = A;
//	stateSpace->nstateArray[A->nalleles][n1] =1;

	g_hash_table_insert(stateSpace->stateHash, A->matString, "TRUE");

	preorderEvolveState(stateSpace, A, n1, n2);

	//qsort(stateSpace->states,stateSpace->nstates, sizeof(stateSpace->states[0]), afsObjectQsortCompare);

	printf("n1: %d n2: %d nstates: %d\n",n1,n2,stateSpace->nstates);
	for(i=0;i<stateSpace->nstates;i++){
		afsObjectPrint(stateSpace->states[i]);
		printf("-----------\n");

	}

	afsStateSpaceFree(stateSpace);
	return(0);
}

void getParameters(int argc,const char **argv){
//	int args;
//	int i;

	if( argc < 2){
		usage();
	}
	n1 = atoi(argv[1]);
	n2 = atoi(argv[2]);
}

void usage(){
	fprintf(stderr,"usage: cmc_stateSpace n1 n2\n");
//	fprintf(stderr,"parameters: \n");
	exit(1);
}

max_stars_repo_name: dortegadelv/IMaDNA

max_issues_repo_name: dortegadelv/IMaDNA