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

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_spblas.h>

/*
create_random_sparse()
  Create a random sparse matrix with approximately
M*N*density non-zero entries

Inputs: M       - number of rows
        N       - number of columns
        density - sparse density \in [0,1]
                  0 = no non-zero entries
                  1 = all m*n entries are filled
        r       - random number generator

Return: pointer to sparse matrix in triplet format (must be freed by caller)

Notes:
1) non-zero matrix entries are uniformly distributed in [0,1]
*/

static gsl_spmatrix *
create_random_sparse(const size_t M, const size_t N, const double density,
                     const gsl_rng *r)
{
  size_t nnzwanted = (size_t) floor(M * N * GSL_MIN(density, 1.0));
  gsl_spmatrix *m = gsl_spmatrix_alloc_nzmax(M, N,
                                             nnzwanted,
                                             GSL_SPMATRIX_TRIPLET);

  while (gsl_spmatrix_nnz(m) < nnzwanted)
    {
      /* generate a random row and column */
      size_t i = gsl_rng_uniform(r) * M;
      size_t j = gsl_rng_uniform(r) * N;

      /* generate random m_{ij} and add it */
      double x = gsl_rng_uniform(r);
      gsl_spmatrix_set(m, i, j, x);
    }

  return m;
} /* create_random_sparse() */

static void
create_random_vector(gsl_vector *v, const gsl_rng *r)
{
  size_t i;

  for (i = 0; i < v->size; ++i)
    {
      double x = gsl_rng_uniform(r);
      gsl_vector_set(v, i, x);
    }
} /* create_random_vector() */

static int
test_vectors(gsl_vector *observed, gsl_vector *expected, const double tol,
             const char *str)
{
  int s = 0;
  size_t N = observed->size;
  size_t i;

  for (i = 0; i < N; ++i)
    {
      double x_obs = gsl_vector_get(observed, i);
      double x_exp = gsl_vector_get(expected, i);

      gsl_test_rel(x_obs, x_exp, tol, "N=%zu i=%zu %s", N, i, str);
    }

  return s;
} /* test_vectors() */

static void
test_dgemv(const size_t N, const size_t M, const double alpha,
           const double beta, const CBLAS_TRANSPOSE_t TransA,
           const gsl_rng *r)
{
  gsl_spmatrix *A = create_random_sparse(M, N, 0.2, r);
  gsl_spmatrix *B, *C;
  gsl_matrix *A_dense = gsl_matrix_alloc(M, N);
  gsl_vector *x, *y, *y_gsl, *y_sp;
  size_t lenX, lenY;

  if (TransA == CblasNoTrans)
    {
      lenX = N;
      lenY = M;
    }
  else
    {
      lenX = M;
      lenY = N;
    }

  x = gsl_vector_alloc(lenX);
  y = gsl_vector_alloc(lenY);
  y_gsl = gsl_vector_alloc(lenY);
  y_sp = gsl_vector_alloc(lenY);

  /* create random dense vectors */
  create_random_vector(x, r);
  create_random_vector(y, r);

  /* copy A into A_dense */
  gsl_spmatrix_sp2d(A_dense, A);

  gsl_vector_memcpy(y_gsl, y);
  gsl_vector_memcpy(y_sp, y);

  /* compute y = alpha*op(A)*x + beta*y0 with gsl */
  gsl_blas_dgemv(TransA, alpha, A_dense, x, beta, y_gsl);

  /* compute y = alpha*op(A)*x + beta*y0 with spblas/triplet */
  gsl_spblas_dgemv(TransA, alpha, A, x, beta, y_sp);

  /* test y_sp = y_gsl */
  test_vectors(y_sp, y_gsl, 1.0e-10, "test_dgemv: triplet format");

  /* compute y = alpha*op(A)*x + beta*y0 with spblas/CCS */
  B = gsl_spmatrix_ccs(A);
  gsl_vector_memcpy(y_sp, y);
  gsl_spblas_dgemv(TransA, alpha, B, x, beta, y_sp);

  /* test y_sp = y_gsl */
  test_vectors(y_sp, y_gsl, 1.0e-10, "test_dgemv: CCS format");

  /* compute y = alpha*op(A)*x + beta*y0 with spblas/CRS */
  C = gsl_spmatrix_crs(A);
  gsl_vector_memcpy(y_sp, y);
  gsl_spblas_dgemv(TransA, alpha, C, x, beta, y_sp);

  /* test y_sp = y_gsl */
  test_vectors(y_sp, y_gsl, 1.0e-10, "test_dgemv: CRS format");

  gsl_spmatrix_free(A);
  gsl_spmatrix_free(B);
  gsl_spmatrix_free(C);
  gsl_matrix_free(A_dense);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(y_gsl);
  gsl_vector_free(y_sp);
} /* test_dgemv() */

static void
test_dgemm(const double alpha, const size_t M, const size_t N,
           const gsl_rng *r)
{
  const size_t max = GSL_MAX(M, N);
  size_t i, j, k;
  gsl_matrix *A_dense = gsl_matrix_alloc(M, max);
  gsl_matrix *B_dense = gsl_matrix_alloc(max, N);
  gsl_matrix *C_dense = gsl_matrix_alloc(M, N);
  gsl_spmatrix *C = gsl_spmatrix_alloc_nzmax(M, N, 1, GSL_SPMATRIX_CCS);

  for (k = 1; k <= max; ++k)
    {
      gsl_matrix_view Ad = gsl_matrix_submatrix(A_dense, 0, 0, M, k);
      gsl_matrix_view Bd = gsl_matrix_submatrix(B_dense, 0, 0, k, N);
      gsl_spmatrix *TA = create_random_sparse(M, k, 0.2, r);
      gsl_spmatrix *TB = create_random_sparse(k, N, 0.2, r);
      gsl_spmatrix *A = gsl_spmatrix_ccs(TA);
      gsl_spmatrix *B = gsl_spmatrix_ccs(TB);

      gsl_spmatrix_set_zero(C);
      gsl_spblas_dgemm(alpha, A, B, C);

      /* make dense matrices and use standard dgemm to multiply them */
      gsl_spmatrix_sp2d(&Ad.matrix, TA);
      gsl_spmatrix_sp2d(&Bd.matrix, TB);
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, alpha, &Ad.matrix,
                     &Bd.matrix, 0.0, C_dense);

      /* compare C and C_dense */
      for (i = 0; i < M; ++i)
        {
          for (j = 0; j < N; ++j)
            {
              double Cij = gsl_spmatrix_get(C, i, j);
              double Dij = gsl_matrix_get(C_dense, i, j);

              gsl_test_rel(Cij, Dij, 1.0e-12, "test_dgemm: _dgemm");
            }
        }

      gsl_spmatrix_free(TA);
      gsl_spmatrix_free(TB);
      gsl_spmatrix_free(A);
      gsl_spmatrix_free(B);
    }

  gsl_spmatrix_free(C);
  gsl_matrix_free(A_dense);
  gsl_matrix_free(B_dense);
  gsl_matrix_free(C_dense);
} /* test_dgemm() */

int
main()
{
  const size_t N_max = 40;
  size_t m, n;
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  for (m = 1; m <= N_max; ++m)
    {
      for (n = 1; n <= N_max; ++n)
        {
          test_dgemv(m, n, 1.0, 0.0, CblasNoTrans, r);
          test_dgemv(m, n, 1.0, 0.0, CblasTrans, r);

          test_dgemv(m, n, 2.4, -0.5, CblasNoTrans, r);
          test_dgemv(m, n, 2.4, -0.5, CblasTrans, r);

          test_dgemv(m, n, 0.1, 10.0, CblasNoTrans, r);
          test_dgemv(m, n, 0.1, 10.0, CblasTrans, r);
        }
    }

  test_dgemm(1.0, 10, 10, r);
  test_dgemm(2.3, 20, 15, r);
  test_dgemm(1.8, 12, 30, r);
  test_dgemm(0.4, 45, 35, r);

  gsl_rng_free(r);

  exit (gsl_test_summary());
} /* main() */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/** @file */
#ifndef __CCL_UTILS_H_INCLUDED__
#define __CCL_UTILS_H_INCLUDED__

#include <gsl/gsl_spline.h>

#define CCL_MIN(a, b)  (((a) < (b)) ? (a) : (b))
#define CCL_MAX(a, b)  (((a) > (b)) ? (a) : (b))

CCL_BEGIN_DECLS

/**
 * Compute bin edges of N-1 linearly spaced bins on the interval [xmin,xmax]
 * @param xmin minimum value of spacing
 * @param xmax maximum value of spacing
 * @param N number of bins plus one (number of bin edges)
 * @return x, bin edges in range [xmin, xmax]
 */
double * ccl_linear_spacing(double xmin, double xmax, int N);

/**
 * Compute bin edges of N-1 logarithmically and then linearly spaced bins on the interval [xmin,xmax]
 * @param xminlog minimum value of spacing
 * @param xmin value when logarithmical ends and linear spacing begins
 * @param xmax maximum value of spacing
 * @param Nlin number of linear bins plus one (number of bin edges)
 * @param Nlog number of logarithmic bins plus one (number of bin edges)
 * @return x, bin edges in range [xminlog, xmax]
 */
double * ccl_linlog_spacing(double xminlog, double xmin, double xmax, int Nlin, int Nlog);

/**
 * Compute bin edges of N-1 logarithmically spaced bins on the interval [xmin,xmax]
 * @param xmin minimum value of spacing
 * @param xmax maximum value of spacing
 * @param N number of bins plus one (number of bin edges)
 * @return x, bin edges in range [xmin, xmax]
 */
double * ccl_log_spacing(double xmin, double xmax, int N);
//Returns array of N logarithmically-spaced values between xmin and xmax

double ccl_j_bessel(int l,double x);
//Spherical Bessel function of order l (adapted from CAMB)

/**
 * Spline wrapper
 * Used to take care of evaluations outside the supported range
 */
typedef struct {
  gsl_interp_accel *intacc; //GSL spline
  gsl_spline *spline;
  double x0,xf; //Interpolation limits
  double y0,yf; //Constant values to use beyond interpolation limit
} SplPar;

SplPar *ccl_spline_init(int n,double *x,double *y,double y0,double yf);

double ccl_spline_eval(double x,SplPar *spl);

void ccl_spline_free(SplPar *spl);

CCL_END_DECLS

#endif

ext: h

max_stars_repo_name: Russell-Jones-OxPhys/CCL

max_issues_repo_name: Russell-Jones-OxPhys/CCL

max_forks_repo_name: Russell-Jones-OxPhys/CCL


text:
// Petsc initialization
#pragma once

#include <hollow/petsc/config.h>
#include <geode/utility/debug.h>
#include <vector>
#include <string>
// Work around configuration issue in old PETSc versions
#define PETSC_RESTRICT PETSC_CXX_RESTRICT
#include <petsc.h>
namespace hollow {

using std::string;
using std::vector;

// Initialize petsc, and automatically finalize on application exit
HOLLOW_EXPORT void petsc_initialize(const string& help, const vector<string>& args);

// Check if petsc is initialized
HOLLOW_EXPORT bool petsc_initialized();

// Ensure that petsc is initialized.  Can be called multiple times, but no help or args support.
HOLLOW_EXPORT void petsc_reinitialize();

// Explicitly finalize petsc.  This happens automatically at program exit if petsc_initialize
// is called, but sometime it is useful to do it sooner.
HOLLOW_EXPORT void petsc_finalize();

// Reset petsc options to the given list.  Any petsc objects already allocated will not change.
HOLLOW_EXPORT void petsc_set_options(const vector<string>& args);

// Add new petsc options without removing old options.  Like argv[0], the first option is ignored.
HOLLOW_EXPORT void petsc_add_options(const vector<string>& args);

}

max_stars_repo_name: otherlab/hollow

max_issues_repo_name: otherlab/hollow

max_forks_repo_name: otherlab/hollow


text:
/*
** Implementation of LISA algorithm
** for statistical inference of fMRI images
**
** 2nd level inference using a GLM design matrix,
**   e.g. for anova designs
**
** G.Lohmann, 2018
*/
#include <viaio/Vlib.h>
#include <viaio/file.h>
#include <viaio/mu.h>
#include <viaio/option.h>
#include <viaio/os.h>
#include <viaio/VImage.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_histogram.h>

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

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

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

extern void   VIsolatedVoxels(VImage src,float threshold);
extern void   VHistogram(gsl_histogram *histogram,VString filename);
extern void   VCheckImage(VImage src);
extern void   FDR(VImage src,VImage dest,gsl_histogram *nullhist,gsl_histogram *realhist,double);
extern double ttest1(double *data1,int n);
extern void   ImageStats(VImage src,double *,double *,double *hmin,double *hmax);
extern void   VBilateralFilter(VImage src,VImage dest,int radius,double var1,double var2,int);
extern double VImageVar(VImage src);
extern void   VImageCount(VImage src);
extern void   VGetHistRange(VImage src,double *hmin,double *hmax);
extern void   VZScale(VImage src,float mode,float stddev);
extern float  VGetMode(VImage src);

extern void XReadExchange(VString filename,int *exchange,int n);
extern void GLM2(VImage *src,gsl_matrix *X,gsl_vector *contrast,int *,int *permtable,int *signtable,int signswitch,VImage dest);
extern gsl_matrix *XRead2ndLevel(VString);
extern gsl_matrix *PseudoInv(gsl_matrix *A,gsl_matrix *B);
extern int SignSwitch (gsl_matrix *X,gsl_vector *contrast,int *);
extern void genperm(long seed,int *exchange,int,int **,int **,int,int,gsl_vector *);
extern void HistoUpdate(VImage,gsl_histogram *);
extern gsl_matrix *VReadCovariates(VString,VBoolean);

int main (int argc, char *argv[])
{
  static VArgVector in_files1;
  static VString  out_filename="";
  static VString  design_filename="";
  static VString  exchange_filename="";
  static VString  nuisance_filename="";
  static VString  mask_filename="";
  static VArgVector cont;
  static VFloat   alpha = 0.05;
  static VShort   radius = 2;
  static VFloat   rvar = 2.0;
  static VFloat   svar = 2.0;
  static VShort   numiter = 2;
  static VShort   numperm = 5000;
  static VLong    seed = 99402622;
  static VBoolean centering = FALSE;
  static VBoolean demean  = TRUE;
  static VBoolean cleanup = TRUE;
  static VShort   nproc = 0;
  static VOptionDescRec options[] = {
    {"in", VStringRepn, 0, & in_files1, VRequiredOpt, NULL,"Input files" },
    {"out", VStringRepn, 1, & out_filename, VRequiredOpt, NULL,"Output file" },
    {"design",VStringRepn,1,(VPointer) &design_filename,VRequiredOpt,NULL,"Design file (2nd level)"},
    {"grp", VStringRepn, 1, & exchange_filename, VOptionalOpt, NULL,"Group Ids needed for exchangeability" },
    {"contrast", VFloatRepn, 0, (VPointer) &cont, VRequiredOpt, NULL, "contrast vector"},
    {"nuisance", VStringRepn, 1, & nuisance_filename, VOptionalOpt, NULL,"Nuisance regressors" },
    {"demean",VBooleanRepn,1,(VPointer) &demean,VOptionalOpt,NULL,"Whether to subtract mean in nuisance regressors"},
    {"alpha",VFloatRepn,1,(VPointer) &alpha,VOptionalOpt,NULL,"FDR significance level"},
    {"perm",VShortRepn,1,(VPointer) &numperm,VOptionalOpt,NULL,"Number of permutations"},
    {"mask", VStringRepn, 1, (VPointer) &mask_filename, VRequiredOpt, NULL, "Mask"},
    {"seed",VLongRepn,1,(VPointer) &seed,VOptionalOpt,NULL,"Seed for random number generation"},
    {"radius",VShortRepn,1,(VPointer) &radius,VOptionalOpt,NULL,"Bilateral parameter (radius in voxels)"},
    {"rvar",VFloatRepn,1,(VPointer) &rvar,VOptionalOpt,NULL,"Bilateral parameter (radiometric)"},
    {"svar",VFloatRepn,1,(VPointer) &svar,VOptionalOpt,NULL,"Bilateral parameter (spatial)"},
    {"filteriterations",VShortRepn,1,(VPointer) &numiter,VOptionalOpt,NULL,"Bilateral parameter (number of iterations)"},
    {"cleanup",VBooleanRepn,1,(VPointer) &cleanup,VOptionalOpt,NULL,"Whether to remove isolated voxels"},
    {"j",VShortRepn,1,(VPointer) &nproc,VOptionalOpt,NULL,"Number of processors to use, '0' to use all"},
  };
  FILE *fp=NULL;
  VString in_filename;
  VAttrList list1=NULL,out_list=NULL,geolist=NULL;
  VImage *src1;
  int i,j,nimages,npix=0;
  double u=0;
  char *prg_name=GetLipsiaName("vlisa_2ndlevel");
  fprintf (stderr, "%s\n", prg_name);

  /* parse command line */
  if (! VParseCommand (VNumber (options), options, & argc, argv)) {
    VReportUsage (argv[0], VNumber (options), options, NULL);
    exit (EXIT_FAILURE);
  }
  if (argc > 1) {
    VReportBadArgs (argc, argv);
    exit (EXIT_FAILURE);
  }

  /* 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 */

  /* images  */
  VImage mask = VReadImageFile(mask_filename);
  if (mask==NULL) VError("Error reading mask file %s",mask_filename);

  nimages = in_files1.number;
  src1 = (VImage *) VCalloc(nimages,sizeof(VImage));
  for (i = 0; i < nimages; i++) {
    in_filename = ((VString *) in_files1.vector)[i];
    list1   = VReadAttrList(in_filename,0L,TRUE,FALSE);
    VMaskMinval(list1,mask,0.0);
    src1[i] = VReadImage(list1);
    if (src1[i] == NULL) VError(" no input image found");
    if (VPixelRepn(src1[i]) != VFloatRepn) VError(" input pixel repn must be float");
    if (i == 0) npix = VImageNPixels(src1[i]);
    else if (npix != VImageNPixels(src1[i])) VError(" inconsistent image dimensions");
    fprintf(stderr," %3d:  %s\n",i,in_filename);

    /* use geometry info from 1st file */
    if (geolist == NULL) geolist = VGetGeoInfo(list1);
  }

  /* get nuisance file */
  gsl_matrix *XN = NULL;
  if (strlen(nuisance_filename) > 1) {
    XN = VReadCovariates(nuisance_filename,demean);
    fprintf(stderr," nuisance covariates dimensions:  %d x %d\n",(int)XN->size1,(int)XN->size2);
    if (XN->size1 != nimages)
      VError(" number of input images (%d) does not match number of rows in nuisance covariates file (%d)",nimages,(int)XN->size1);
  }

  /* get design file and its inverse */
  gsl_matrix *X0 = XRead2ndLevel(design_filename);
  fprintf(stderr," design file dimensions:  %d x %d\n",(int)X0->size1,(int)X0->size2);
  if (X0->size1 != nimages)
    VError(" number of input images (%d) does not match number of rows in design file (%d)",nimages,(int)X0->size1);

  /* read contrast vector */
  if (cont.number != X0->size2)
    VError(" Dimension of contrast vector (%d) does not match number of columns in design file (%d)",cont.number,X0->size2);
  gsl_vector *contrast0 = gsl_vector_alloc(cont.number);
  for (i=0; i < cont.number; i++) {
    double u = ((VFloat *)cont.vector)[i];
    gsl_vector_set(contrast0, i, u);
  }

  /* concatenate design matrix and nuisance covariates */
  gsl_matrix *X = NULL;
  gsl_vector *contrast = NULL;
  int *permflag = NULL;

  /* include nuisance covariates, if present */
  if (XN != NULL) {
    permflag = (int *) VCalloc(XN->size2 + X0->size2,sizeof(int));
    for (j=0; j<XN->size2 + X0->size2; j++) permflag[j] = 0;
    for (j=0; j<X0->size2; j++) permflag[j] = 1;

    X = gsl_matrix_calloc(X0->size1, X0->size2 + XN->size2);
    for (i=0; i<X->size1; i++) {
      for (j=0; j<X0->size2; j++) {
	u = gsl_matrix_get(X0,i,j);
	gsl_matrix_set(X,i,j,u);
      }
      for (j=0; j<XN->size2; j++) {
	u = gsl_matrix_get(XN,i,j);
	gsl_matrix_set(X,i,j+X0->size2,u);
      }
    }
    contrast = gsl_vector_calloc(cont.number + XN->size2);
    gsl_vector_set_zero(contrast);
    for (j=0; j<contrast0->size; j++) gsl_vector_set(contrast,j,contrast0->data[j]);
  }

  /* no nuisance covariates */
  else {
    if (X0->size2 == 1) VError(" Please use 'vlisa_onesample' for onesample tests");
    permflag = (int *) VCalloc(X0->size2,sizeof(int));
    for (j=0; j<X0->size2; j++) permflag[j] = 1;
    X = gsl_matrix_calloc(X0->size1,X0->size2);
    gsl_matrix_memcpy(X,X0);
    contrast = gsl_vector_calloc(cont.number);
    gsl_vector_memcpy(contrast,contrast0);
  }

  /* read exchangeability information */
  int *exchange = (int *) VCalloc(X->size1,sizeof(int));
  for (i=0; i<X->size1; i++) exchange[i] = 1;  /* no restrictions w.r.t. exchangeability */
  if (strlen(exchange_filename) > 1) {
    XReadExchange(exchange_filename,exchange,(int)X->size1);
  }

  /* ini random permutations and sign switching */
  int signswitch = SignSwitch(X,contrast,permflag);
  int nperm=0;
  int **permtable = (int **) VCalloc(numperm,sizeof(int *));
  int **signtable = (int **) VCalloc(numperm,sizeof(int *));
  genperm((long)seed,exchange,signswitch,permtable,signtable,(int)nimages,(int)numperm,contrast);

  /* estimate null variance based on first 30 permutations */
  double hmin=0,hmax=0;
  float stddev=1.0;
  if (numperm > 0) {
    int tstperm = 30;
    if (tstperm > numperm) tstperm = numperm;
    double varsum=0,nx=0;

#pragma omp parallel for shared(src1,permtable) schedule(dynamic)
    for (nperm = 0; nperm < tstperm; nperm++) {
      VImage zmap = VCreateImageLike(src1[0]);
      GLM2(src1,X,contrast,permflag,permtable[nperm],signtable[nperm],signswitch,zmap);
#pragma omp critical
      {
	varsum += VImageVar(zmap);
	nx++;
      }
      VDestroyImage(zmap);
    }
    double meanvar = varsum/nx;
    stddev = sqrt(meanvar);
  }

  /* no permutation */
  int *nopermtable = (int *) VCalloc(nimages,sizeof(int));
  int *nosigntable = (int *) VCalloc(nimages,sizeof(int));
  for (i=0; i<nimages; i++) nopermtable[i] = i;
  for (i=0; i<nimages; i++) nosigntable[i] = 1;
  VImage dst1  = VCreateImageLike (src1[0]);
  VImage zmap1 = VCreateImageLike(src1[0]);
  GLM2(src1,X,contrast,permflag,nopermtable,nosigntable,signswitch,zmap1);

  if (numperm == 0) {
    double z = VImageVar(zmap1);
    stddev = (float)(sqrt(z)); /* update stddev */
  }
  float mode=0;
  if (centering) mode = VGetMode(zmap1);
  if (numperm > 0) VZScale(zmap1,mode,stddev);
  VBilateralFilter(zmap1,dst1,(int)radius,(double)(rvar),(double)svar,(int)numiter);

  /* ini histograms */
  VGetHistRange(dst1,&hmin,&hmax);
  size_t nbins = 20000;
  gsl_histogram *hist0 = gsl_histogram_alloc (nbins);
  gsl_histogram_set_ranges_uniform (hist0,hmin,hmax);
  gsl_histogram *histz = gsl_histogram_alloc (nbins);
  gsl_histogram_set_ranges_uniform (histz,hmin,hmax);
  HistoUpdate(dst1,histz);

  /* random permutations */
#pragma omp parallel for shared(src1,permtable) schedule(dynamic)
  for (nperm = 0; nperm < numperm; nperm++) {
    if (nperm%20 == 0) fprintf(stderr," perm  %4d  of  %d\r",nperm,(int)numperm);

    VImage zmap = VCreateImageLike(src1[0]);
    VImage dst  = VCreateImageLike (zmap);
    GLM2(src1,X,contrast,permflag,permtable[nperm],signtable[nperm],signswitch,zmap);
    float mode=0;
    if (centering) mode = VGetMode(zmap);
    VZScale(zmap,mode,stddev);
    VBilateralFilter(zmap,dst,(int)radius,(double)(rvar),(double)svar,(int)numiter);

#pragma omp critical
    {
      HistoUpdate(dst,hist0);
    }
    VDestroyImage(dst);
    VDestroyImage(zmap);
  }

  /* apply fdr */
  VImage fdrimage = VCopyImage (dst1,NULL,VAllBands);
  if (numperm > 0) {
    FDR(dst1,fdrimage,hist0,histz,(double)alpha);

    if (cleanup && alpha < 1.0) {
      VIsolatedVoxels(fdrimage,(float)(1.0-alpha));
    }
  }

  /* write output to disk */
  out_list = VCreateAttrList ();
  VHistory(VNumber(options),options,prg_name,&list1,&out_list);
  VSetGeoInfo(geolist,out_list);
  VAppendAttr (out_list,"image",NULL,VImageRepn,fdrimage);
  fp = VOpenOutputFile (out_filename, TRUE);
  if (! VWriteFile (fp, out_list)) exit (1);
  fclose(fp);
  fprintf (stderr, "\n%s: done.\n", argv[0]);
  exit(0);
}

ext: c

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <gsl/blas/cblas.h>

/**
	All the methods below are obtained directly from the ATLAS implementation of BLAS, available at: http://www.netlib.org/clapack/cblas/
*/

#include "f2c.h"

doublereal ddot(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy){

    /* System generated locals */
    integer i__1;
    doublereal ret_val;

    /* Local variables */
    static integer i, m;
    static doublereal dtemp;
    static integer ix, iy, mp1;

/*     forms the dot product of two vectors.
       uses unrolled loops for increments equal to one.
       jack dongarra, linpack, 3/11/78.
       modified 12/3/93, array(1) declarations changed to array(*)

   Parameter adjustments
       Function Body */
#define DY(I) dy[(I)-1]
#define DX(I) dx[(I)-1]

    ret_val = 0.;
    dtemp = 0.;
    if (*n <= 0) {
	return ret_val;
    }
    if (*incx == 1 && *incy == 1) {
	goto L20;
    }

/*        code for unequal increments or equal increments
            not equal to 1 */

    ix = 1;
    iy = 1;
    if (*incx < 0) {
	ix = (-(*n) + 1) * *incx + 1;
    }
    if (*incy < 0) {
	iy = (-(*n) + 1) * *incy + 1;
    }
    i__1 = *n;
    for (i = 1; i <= *n; ++i) {
	dtemp += DX(ix) * DY(iy);
	ix += *incx;
	iy += *incy;
/* L10: */
    }
    ret_val = dtemp;
    return ret_val;

/*        code for both increments equal to 1

          clean-up loop */

L20:
    m = *n % 5;
    if (m == 0) {
	goto L40;
    }
    i__1 = m;
    for (i = 1; i <= m; ++i) {
	dtemp += DX(i) * DY(i);
/* L30: */
    }
    if (*n < 5) {
	goto L60;
    }
L40:
    mp1 = m + 1;
    i__1 = *n;
    for (i = mp1; i <= *n; i += 5) {
	dtemp = dtemp + DX(i) * DY(i) + DX(i + 1) * DY(i + 1) + DX(i + 2) *
		DY(i + 2) + DX(i + 3) * DY(i + 3) + DX(i + 4) * DY(i + 4);
/* L50: */
    }
L60:
    ret_val = dtemp;
    return ret_val;
} /* ddot_ */

doublereal dnrm2(integer *n, doublereal *x, integer *incx)
{
    /* System generated locals */
    integer i__1, i__2;
    doublereal ret_val, d__1;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static doublereal norm, scale, absxi;
    static integer ix;
    static doublereal ssq;

/*  DNRM2 returns the euclidean norm of a vector via the function
    name, so that

       DNRM2 := sqrt( x'*x )

    -- This version written on 25-October-1982.
       Modified on 14-October-1993 to inline the call to DLASSQ.
       Sven Hammarling, Nag Ltd.

   Parameter adjustments
       Function Body */
#define X(I) x[(I)-1]

    if (*n < 1 || *incx < 1) {
	norm = 0.;
    } else if (*n == 1) {
	norm = abs(X(1));
    } else {
	scale = 0.;
	ssq = 1.;
/*        The following loop is equivalent to this call to the LAPACK

          auxiliary routine:
          CALL DLASSQ( N, X, INCX, SCALE, SSQ ) */

	i__1 = (*n - 1) * *incx + 1;
	i__2 = *incx;
	for (ix = 1; *incx < 0 ? ix >= (*n-1)**incx+1 : ix <= (*n-1)**incx+1; ix += *incx) {
	    if (X(ix) != 0.) {
		absxi = (d__1 = X(ix), abs(d__1));
		if (scale < absxi) {
/* Computing 2nd power */
		    d__1 = scale / absxi;
		    ssq = ssq * (d__1 * d__1) + 1.;
		    scale = absxi;
		} else {
/* Computing 2nd power */
		    d__1 = absxi / scale;
		    ssq += d__1 * d__1;
		}
	    }
/* L10: */
	}
	norm = scale * sqrt(ssq);
    }

    ret_val = norm;
    return ret_val;

/*     End of DNRM2. */

} /* dnrm2_ */

int dswap(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
{

    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer i, m;
    static doublereal dtemp;
    static integer ix, iy, mp1;

/*     interchanges two vectors.
       uses unrolled loops for increments equal one.
       jack dongarra, linpack, 3/11/78.
       modified 12/3/93, array(1) declarations changed to array(*)

   Parameter adjustments
       Function Body */
#define DY(I) dy[(I)-1]
#define DX(I) dx[(I)-1]

    if (*n <= 0) {
	return 0;
    }
    if (*incx == 1 && *incy == 1) {
	goto L20;
    }

/*       code for unequal increments or equal increments not equal
           to 1 */

    ix = 1;
    iy = 1;
    if (*incx < 0) {
	ix = (-(*n) + 1) * *incx + 1;
    }
    if (*incy < 0) {
	iy = (-(*n) + 1) * *incy + 1;
    }
    i__1 = *n;
    for (i = 1; i <= *n; ++i) {
	dtemp = DX(ix);
	DX(ix) = DY(iy);
	DY(iy) = dtemp;
	ix += *incx;
	iy += *incy;
/* L10: */
    }
    return 0;

/*       code for both increments equal to 1

         clean-up loop */

L20:
    m = *n % 3;
    if (m == 0) {
	goto L40;
    }
    i__1 = m;
    for (i = 1; i <= m; ++i) {
	dtemp = DX(i);
	DX(i) = DY(i);
	DY(i) = dtemp;
/* L30: */
    }
    if (*n < 3) {
	return 0;
    }
L40:
    mp1 = m + 1;
    i__1 = *n;
    for (i = mp1; i <= *n; i += 3) {
	dtemp = DX(i);
	DX(i) = DY(i);
	DY(i) = dtemp;
	dtemp = DX(i + 1);
	DX(i + 1) = DY(i + 1);
	DY(i + 1) = dtemp;
	dtemp = DX(i + 2);
	DX(i + 2) = DY(i + 2);
	DY(i + 2) = dtemp;
/* L50: */
    }
    return 0;
} /* dswap_ */

int dcopy(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
{

    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer i, m, ix, iy, mp1;

/*     copies a vector, x, to a vector, y.
       uses unrolled loops for increments equal to one.
       jack dongarra, linpack, 3/11/78.
       modified 12/3/93, array(1) declarations changed to array(*)

   Parameter adjustments
       Function Body */
#define DY(I) dy[(I)-1]
#define DX(I) dx[(I)-1]

    if (*n <= 0) {
	return 0;
    }
    if (*incx == 1 && *incy == 1) {
	goto L20;
    }

/*        code for unequal increments or equal increments
            not equal to 1 */

    ix = 1;
    iy = 1;
    if (*incx < 0) {
	ix = (-(*n) + 1) * *incx + 1;
    }
    if (*incy < 0) {
	iy = (-(*n) + 1) * *incy + 1;
    }
    i__1 = *n;
    for (i = 1; i <= *n; ++i) {
	DY(iy) = DX(ix);
	ix += *incx;
	iy += *incy;
/* L10: */
    }
    return 0;

/*        code for both increments equal to 1

          clean-up loop */

L20:
    m = *n % 7;
    if (m == 0) {
	goto L40;
    }
    i__1 = m;
    for (i = 1; i <= m; ++i) {
	DY(i) = DX(i);
/* L30: */
    }
    if (*n < 7) {
	return 0;
    }
L40:
    mp1 = m + 1;
    i__1 = *n;
    for (i = mp1; i <= *n; i += 7) {
	DY(i) = DX(i);
	DY(i + 1) = DX(i + 1);
	DY(i + 2) = DX(i + 2);
	DY(i + 3) = DX(i + 3);
	DY(i + 4) = DX(i + 4);
	DY(i + 5) = DX(i + 5);
	DY(i + 6) = DX(i + 6);
/* L50: */
    }
    return 0;
} /* dcopy_ */

int dgemm(char *transa, char *transb, integer *m, integer *	n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc)
{

    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
	    i__3;

    /* Local variables */
    static integer info;
    static logical nota, notb;
    static doublereal temp;
    static integer i, j, l, ncola;
    extern logical lsame_(char *, char *);
    static integer nrowa, nrowb;
    extern /* Subroutine */ int xerbla_(char *, integer *);

/*  Purpose
    =======

    DGEMM  performs one of the matrix-matrix operations

       C := alpha*op( A )*op( B ) + beta*C,

    where  op( X ) is one of

       op( X ) = X   or   op( X ) = X',

    alpha and beta are scalars, and A, B and C are matrices, with op( A )

    an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.

    Parameters
    ==========

    TRANSA - CHARACTER*1.
             On entry, TRANSA specifies the form of op( A ) to be used in

             the matrix multiplication as follows:

                TRANSA = 'N' or 'n',  op( A ) = A.

                TRANSA = 'T' or 't',  op( A ) = A'.

                TRANSA = 'C' or 'c',  op( A ) = A'.

             Unchanged on exit.

    TRANSB - CHARACTER*1.
             On entry, TRANSB specifies the form of op( B ) to be used in

             the matrix multiplication as follows:

                TRANSB = 'N' or 'n',  op( B ) = B.

                TRANSB = 'T' or 't',  op( B ) = B'.

                TRANSB = 'C' or 'c',  op( B ) = B'.

             Unchanged on exit.

    M      - INTEGER.
             On entry,  M  specifies  the number  of rows  of the  matrix

             op( A )  and of the  matrix  C.  M  must  be at least  zero.

             Unchanged on exit.

    N      - INTEGER.
             On entry,  N  specifies the number  of columns of the matrix

             op( B ) and the number of columns of the matrix C. N must be

             at least zero.
             Unchanged on exit.

    K      - INTEGER.
             On entry,  K  specifies  the number of columns of the matrix

             op( A ) and the number of rows of the matrix op( B ). K must

             be at least  zero.
             Unchanged on exit.

    ALPHA  - DOUBLE PRECISION.
             On entry, ALPHA specifies the scalar alpha.
             Unchanged on exit.

    A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is

             k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
             Before entry with  TRANSA = 'N' or 'n',  the leading  m by k

             part of the array  A  must contain the matrix  A,  otherwise

             the leading  k by m  part of the array  A  must contain  the

             matrix A.
             Unchanged on exit.

    LDA    - INTEGER.
             On entry, LDA specifies the first dimension of A as declared

             in the calling (sub) program. When  TRANSA = 'N' or 'n' then

             LDA must be at least  max( 1, m ), otherwise  LDA must be at

             least  max( 1, k ).
             Unchanged on exit.

    B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is

             n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
             Before entry with  TRANSB = 'N' or 'n',  the leading  k by n

             part of the array  B  must contain the matrix  B,  otherwise

             the leading  n by k  part of the array  B  must contain  the

             matrix B.
             Unchanged on exit.

    LDB    - INTEGER.
             On entry, LDB specifies the first dimension of B as declared

             in the calling (sub) program. When  TRANSB = 'N' or 'n' then

             LDB must be at least  max( 1, k ), otherwise  LDB must be at

             least  max( 1, n ).
             Unchanged on exit.

    BETA   - DOUBLE PRECISION.
             On entry,  BETA  specifies the scalar  beta.  When  BETA  is

             supplied as zero then C need not be set on input.
             Unchanged on exit.

    C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
             Before entry, the leading  m by n  part of the array  C must

             contain the matrix  C,  except when  beta  is zero, in which

             case C need not be set on entry.
             On exit, the array  C  is overwritten by the  m by n  matrix

             ( alpha*op( A )*op( B ) + beta*C ).

    LDC    - INTEGER.
             On entry, LDC specifies the first dimension of C as declared

             in  the  calling  (sub)  program.   LDC  must  be  at  least

             max( 1, m ).
             Unchanged on exit.

    Level 3 Blas routine.

    -- Written on 8-February-1989.
       Jack Dongarra, Argonne National Laboratory.
       Iain Duff, AERE Harwell.
       Jeremy Du Croz, Numerical Algorithms Group Ltd.
       Sven Hammarling, Numerical Algorithms Group Ltd.

       Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not

       transposed and set  NROWA, NCOLA and  NROWB  as the number of rows

       and  columns of  A  and the  number of  rows  of  B  respectively.

   Parameter adjustments
       Function Body */

#define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]
#define B(I,J) b[(I)-1 + ((J)-1)* ( *ldb)]
#define C(I,J) c[(I)-1 + ((J)-1)* ( *ldc)]

    nota = lsame_(transa, "N");
    notb = lsame_(transb, "N");
    if (nota) {
	nrowa = *m;
	ncola = *k;
    } else {
	nrowa = *k;
	ncola = *m;
    }
    if (notb) {
	nrowb = *k;
    } else {
	nrowb = *n;
    }

/*     Test the input parameters. */

    info = 0;
    if (! nota && ! lsame_(transa, "C") && ! lsame_(transa, "T")) {
	info = 1;
    } else if (! notb && ! lsame_(transb, "C") && ! lsame_(transb,
	    "T")) {
	info = 2;
    } else if (*m < 0) {
	info = 3;
    } else if (*n < 0) {
	info = 4;
    } else if (*k < 0) {
	info = 5;
    } else if (*lda < max(1,nrowa)) {
	info = 8;
    } else if (*ldb < max(1,nrowb)) {
	info = 10;
    } else if (*ldc < max(1,*m)) {
	info = 13;
    }
    if (info != 0) {
	xerbla_("DGEMM ", &info);
	return 0;
    }

/*     Quick return if possible. */

    if (*m == 0 || *n == 0 || (*alpha == 0. || *k == 0) && *beta == 1.) {
	return 0;
    }

/*     And if  alpha.eq.zero. */

    if (*alpha == 0.) {
	if (*beta == 0.) {
	    i__1 = *n;
	    for (j = 1; j <= *n; ++j) {
		i__2 = *m;
		for (i = 1; i <= *m; ++i) {
		    C(i,j) = 0.;
/* L10: */
		}
/* L20: */
	    }
	} else {
	    i__1 = *n;
	    for (j = 1; j <= *n; ++j) {
		i__2 = *m;
		for (i = 1; i <= *m; ++i) {
		    C(i,j) = *beta * C(i,j);
/* L30: */
		}
/* L40: */
	    }
	}
	return 0;
    }

/*     Start the operations. */

    if (notb) {
	if (nota) {

/*           Form  C := alpha*A*B + beta*C. */

	    i__1 = *n;
	    for (j = 1; j <= *n; ++j) {
		if (*beta == 0.) {
		    i__2 = *m;
		    for (i = 1; i <= *m; ++i) {
			C(i,j) = 0.;
/* L50: */
		    }
		} else if (*beta != 1.) {
		    i__2 = *m;
		    for (i = 1; i <= *m; ++i) {
			C(i,j) = *beta * C(i,j);
/* L60: */
		    }
		}
		i__2 = *k;
		for (l = 1; l <= *k; ++l) {
		    if (B(l,j) != 0.) {
			temp = *alpha * B(l,j);
			i__3 = *m;
			for (i = 1; i <= *m; ++i) {
			    C(i,j) += temp * A(i,l);
/* L70: */
			}
		    }
/* L80: */
		}
/* L90: */
	    }
	} else {

/*           Form  C := alpha*A'*B + beta*C */

	    i__1 = *n;
	    for (j = 1; j <= *n; ++j) {
		i__2 = *m;
		for (i = 1; i <= *m; ++i) {
		    temp = 0.;
		    i__3 = *k;
		    for (l = 1; l <= *k; ++l) {
			temp += A(l,i) * B(l,j);
/* L100: */
		    }
		    if (*beta == 0.) {
			C(i,j) = *alpha * temp;
		    } else {
			C(i,j) = *alpha * temp + *beta * C(i,j);
		    }
/* L110: */
		}
/* L120: */
	    }
	}
    } else {
	if (nota) {

/*           Form  C := alpha*A*B' + beta*C */

	    i__1 = *n;
	    for (j = 1; j <= *n; ++j) {
		if (*beta == 0.) {
		    i__2 = *m;
		    for (i = 1; i <= *m; ++i) {
			C(i,j) = 0.;
/* L130: */
		    }
		} else if (*beta != 1.) {
		    i__2 = *m;
		    for (i = 1; i <= *m; ++i) {
			C(i,j) = *beta * C(i,j);
/* L140: */
		    }
		}
		i__2 = *k;
		for (l = 1; l <= *k; ++l) {
		    if (B(j,l) != 0.) {
			temp = *alpha * B(j,l);
			i__3 = *m;
			for (i = 1; i <= *m; ++i) {
			    C(i,j) += temp * A(i,l);
/* L150: */
			}
		    }
/* L160: */
		}
/* L170: */
	    }
	} else {

/*           Form  C := alpha*A'*B' + beta*C */

	    i__1 = *n;
	    for (j = 1; j <= *n; ++j) {
		i__2 = *m;
		for (i = 1; i <= *m; ++i) {
		    temp = 0.;
		    i__3 = *k;
		    for (l = 1; l <= *k; ++l) {
			temp += A(l,i) * B(j,l);
/* L180: */
		    }
		    if (*beta == 0.) {
			C(i,j) = *alpha * temp;
		    } else {
			C(i,j) = *alpha * temp + *beta * C(i,j);
		    }
/* L190: */
		}
/* L200: */
	    }
	}
    }

    return 0;

/*     End of DGEMM . */

} /* dgemm_ */

int dgemv(char *trans, integer *m, integer *n, doublereal *	alpha, doublereal *a, integer *lda, doublereal *x, integer *incx, doublereal *beta, doublereal *y, integer *incy)
{

    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2;

    /* Local variables */
    static integer info;
    static doublereal temp;
    static integer lenx, leny, i, j;
    extern logical lsame_(char *, char *);
    static integer ix, iy, jx, jy, kx, ky;
    extern /* Subroutine */ int xerbla_(char *, integer *);

/*  Purpose
    =======

    DGEMV  performs one of the matrix-vector operations

       y := alpha*A*x + beta*y,   or   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.

    Parameters
    ==========

    TRANS  - CHARACTER*1.
             On entry, TRANS specifies the operation to be performed as
             follows:

                TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.

                TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.

                TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.

             Unchanged on exit.

    M      - INTEGER.
             On entry, M specifies the number of rows of the matrix A.
             M must be at least zero.
             Unchanged on exit.

    N      - INTEGER.
             On entry, N specifies the number of columns of the matrix A.

             N must be at least zero.
             Unchanged on exit.

    ALPHA  - DOUBLE PRECISION.
             On entry, ALPHA specifies the scalar alpha.
             Unchanged on exit.

    A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
             Before entry, the leading m by n part of the array A must
             contain the matrix of coefficients.
             Unchanged on exit.

    LDA    - INTEGER.
             On entry, LDA specifies the first dimension of A as declared

             in the calling (sub) program. LDA must be at least
             max( 1, m ).
             Unchanged on exit.

    X      - DOUBLE PRECISION array of DIMENSION at least
             ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
             and at least
             ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
             Before entry, the incremented array X must contain the
             vector x.
             Unchanged on exit.

    INCX   - INTEGER.
             On entry, INCX specifies the increment for the elements of
             X. INCX must not be zero.
             Unchanged on exit.

    BETA   - DOUBLE PRECISION.
             On entry, BETA specifies the scalar beta. When BETA is
             supplied as zero then Y need not be set on input.
             Unchanged on exit.

    Y      - DOUBLE PRECISION array of DIMENSION at least
             ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
             and at least
             ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
             Before entry with BETA non-zero, the incremented array Y
             must contain the vector y. On exit, Y is overwritten by the

             updated vector y.

    INCY   - INTEGER.
             On entry, INCY specifies the increment for the elements of
             Y. INCY must not be zero.
             Unchanged on exit.

    Level 2 Blas routine.

    -- Written on 22-October-1986.
       Jack Dongarra, Argonne National Lab.
       Jeremy Du Croz, Nag Central Office.
       Sven Hammarling, Nag Central Office.
       Richard Hanson, Sandia National Labs.

       Test the input parameters.

   Parameter adjustments
       Function Body */
#define X(I) x[(I)-1]
#define Y(I) y[(I)-1]

#define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]

    info = 0;
    if (! lsame_(trans, "N") && ! lsame_(trans, "T") && !
	    lsame_(trans, "C")) {
	info = 1;
    } else if (*m < 0) {
	info = 2;
    } else if (*n < 0) {
	info = 3;
    } else if (*lda < max(1,*m)) {
	info = 6;
    } else if (*incx == 0) {
	info = 8;
    } else if (*incy == 0) {
	info = 11;
    }
    if (info != 0) {
	xerbla_("DGEMV ", &info);
	return 0;
    }

/*     Quick return if possible. */

    if (*m == 0 || *n == 0 || *alpha == 0. && *beta == 1.) {
	return 0;
    }

/*     Set  LENX  and  LENY, the lengths of the vectors x and y, and set

       up the start points in  X  and  Y. */

    if (lsame_(trans, "N")) {
	lenx = *n;
	leny = *m;
    } else {
	lenx = *m;
	leny = *n;
    }
    if (*incx > 0) {
	kx = 1;
    } else {
	kx = 1 - (lenx - 1) * *incx;
    }
    if (*incy > 0) {
	ky = 1;
    } else {
	ky = 1 - (leny - 1) * *incy;
    }

/*     Start the operations. In this version the elements of A are
       accessed sequentially with one pass through A.

       First form  y := beta*y. */

    if (*beta != 1.) {
	if (*incy == 1) {
	    if (*beta == 0.) {
		i__1 = leny;
		for (i = 1; i <= leny; ++i) {
		    Y(i) = 0.;
/* L10: */
		}
	    } else {
		i__1 = leny;
		for (i = 1; i <= leny; ++i) {
		    Y(i) = *beta * Y(i);
/* L20: */
		}
	    }
	} else {
	    iy = ky;
	    if (*beta == 0.) {
		i__1 = leny;
		for (i = 1; i <= leny; ++i) {
		    Y(iy) = 0.;
		    iy += *incy;
/* L30: */
		}
	    } else {
		i__1 = leny;
		for (i = 1; i <= leny; ++i) {
		    Y(iy) = *beta * Y(iy);
		    iy += *incy;
/* L40: */
		}
	    }
	}
    }
    if (*alpha == 0.) {
	return 0;
    }
    if (lsame_(trans, "N")) {

/*        Form  y := alpha*A*x + y. */

	jx = kx;
	if (*incy == 1) {
	    i__1 = *n;
	    for (j = 1; j <= *n; ++j) {
		if (X(jx) != 0.) {
		    temp = *alpha * X(jx);
		    i__2 = *m;
		    for (i = 1; i <= *m; ++i) {
			Y(i) += temp * A(i,j);
/* L50: */
		    }
		}
		jx += *incx;
/* L60: */
	    }
	} else {
	    i__1 = *n;
	    for (j = 1; j <= *n; ++j) {
		if (X(jx) != 0.) {
		    temp = *alpha * X(jx);
		    iy = ky;
		    i__2 = *m;
		    for (i = 1; i <= *m; ++i) {
			Y(iy) += temp * A(i,j);
			iy += *incy;
/* L70: */
		    }
		}
		jx += *incx;
/* L80: */
	    }
	}
    } else {

/*        Form  y := alpha*A'*x + y. */

	jy = ky;
	if (*incx == 1) {
	    i__1 = *n;
	    for (j = 1; j <= *n; ++j) {
		temp = 0.;
		i__2 = *m;
		for (i = 1; i <= *m; ++i) {
		    temp += A(i,j) * X(i);
/* L90: */
		}
		Y(jy) += *alpha * temp;
		jy += *incy;
/* L100: */
	    }
	} else {
	    i__1 = *n;
	    for (j = 1; j <= *n; ++j) {
		temp = 0.;
		ix = kx;
		i__2 = *m;
		for (i = 1; i <= *m; ++i) {
		    temp += A(i,j) * X(ix);
		    ix += *incx;
/* L110: */
		}
		Y(jy) += *alpha * temp;
		jy += *incy;
/* L120: */
	    }
	}
    }

    return 0;

/*     End of DGEMV . */

} /* dgemv_ */

logical lsame_(char *ca, char *cb)
{

    /* System generated locals */
    logical ret_val;

    /* Local variables */
    static integer inta, intb, zcode;

/*  -- LAPACK auxiliary routine (version 2.0) --
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
       Courant Institute, Argonne National Lab, and Rice University
       January 31, 1994

    Purpose
    =======

    LSAME returns .TRUE. if CA is the same letter as CB regardless of
    case.

    Arguments
    =========

    CA      (input) CHARACTER*1
    CB      (input) CHARACTER*1
            CA and CB specify the single characters to be compared.

   =====================================================================

       Test if the characters are equal */

    ret_val = *(unsigned char *)ca == *(unsigned char *)cb;
    if (ret_val) {
	return ret_val;
    }

/*     Now test for equivalence if both characters are alphabetic. */

    zcode = 'Z';

/*     Use 'Z' rather than 'A' so that ASCII can be detected on Prime
       machines, on which ICHAR returns a value with bit 8 set.
       ICHAR('A') on Prime machines returns 193 which is the same as
       ICHAR('A') on an EBCDIC machine. */

    inta = *(unsigned char *)ca;
    intb = *(unsigned char *)cb;

    if (zcode == 90 || zcode == 122) {

/*        ASCII is assumed - ZCODE is the ASCII code of either lower o
r
          upper case 'Z'. */

	if (inta >= 97 && inta <= 122) {
	    inta += -32;
	}
	if (intb >= 97 && intb <= 122) {
	    intb += -32;
	}

    } else if (zcode == 233 || zcode == 169) {

/*        EBCDIC is assumed - ZCODE is the EBCDIC code of either lower
 or
          upper case 'Z'. */

	if (inta >= 129 && inta <= 137 || inta >= 145 && inta <= 153 || inta
		>= 162 && inta <= 169) {
	    inta += 64;
	}
	if (intb >= 129 && intb <= 137 || intb >= 145 && intb <= 153 || intb
		>= 162 && intb <= 169) {
	    intb += 64;
	}

    } else if (zcode == 218 || zcode == 250) {

/*        ASCII is assumed, on Prime machines - ZCODE is the ASCII cod
e
          plus 128 of either lower or upper case 'Z'. */

	if (inta >= 225 && inta <= 250) {
	    inta += -32;
	}
	if (intb >= 225 && intb <= 250) {
	    intb += -32;
	}
    }
    ret_val = inta == intb;

/*     RETURN

       End of LSAME */

    return ret_val;
} /* lsame_ */
int xerbla_(char *srname, integer *info)
{
/*  -- LAPACK auxiliary routine (version 2.0) --
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
       Courant Institute, Argonne National Lab, and Rice University
       September 30, 1994

    Purpose
    =======

    XERBLA  is an error handler for the LAPACK routines.
    It is called by an LAPACK routine if an input parameter has an
    invalid value.  A message is printed and execution stops.

    Installers may consider modifying the STOP statement in order to
    call system-specific exception-handling facilities.

    Arguments
    =========

    SRNAME  (input) CHARACTER*6
            The name of the routine which called XERBLA.

    INFO    (input) INTEGER
            The position of the invalid parameter in the parameter list

            of the calling routine.

   =====================================================================
*/

    printf("** On entry to %6s, parameter number %2i had an illegal value\n",
		srname, *info);

/*     End of XERBLA */

    return 0;
} /* xerbla_ */

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
/*
 * CircularDB implementation for time series data.
 *
 * Copyright (c) 2007-2009 Powerset, Inc
 * Copyright (c) Dan Grillo, Manish Dubey, Dan Sully
 *
 * All rights reserved.
 */

#include "config.h"

#include <ctype.h>
#include <errno.h>
#include <fcntl.h>
#include <float.h>
#include <inttypes.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <time.h>
#include <unistd.h>

/* For the aggregation interface */
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>

#include <circulardb_interface.h>

/* Future win32 support */
#ifndef O_BINARY
#define O_BINARY 0
#endif

/* Save away errno so it doesn't get overwritten */
int cdb_error(void) {
    int cdb_error = errno;
    return cdb_error;
}

static void _print_record(FILE *fh, cdb_time_t time, double value, const char *date_format) {

    if (date_format == NULL || strcmp(date_format, "") == 0) {

        fprintf(fh, "%d %.8g\n", (int)(time), value);

    } else {

        char formatted[256];
        time_t stime = (time_t)time;

        strftime(formatted, sizeof(formatted), date_format, localtime(&stime));

        fprintf(fh, "%d [%s] %.8g\n", (int)(time), formatted, value);
    }
}

static uint64_t _physical_record_for_logical_record(cdb_header_t *header, int64_t logical_record) {

    uint64_t physical_record = 0;

    /* -ve indicates nth record from the end
       +ve indicates nth record from the beginning (zero based counting for
          this case)

       -3 = seek to 3rd record from the end
        5 = seek to 6th record from the beginning
        0 = seek to first record (start) from the beginning

       change the request for -nth record (nth record from the end)
       to a request for mth record from the beginning, where:

       m = N + n (n is -ve)

       N : total num of records. */

    if (logical_record < 0) {
        int64_t rec_num_from_start = header->num_records + logical_record;

        /* if asking for records more than there are in the db, just give as many as we can */
        logical_record = rec_num_from_start >= 0 ? rec_num_from_start : 0;
    }

    if (logical_record >= header->num_records) {
#ifdef DEBUG
        printf("Can't seek to record [%"PRIu64"] with only [%"PRIu64"] in db\n", logical_record, header->num_records);
#endif
        return 0;
    }

    /* DRK unclear behavior here. Looks like if I specify -N where N >
       num_records, I get pointed to #0, but if I specify +N where N >
       num_records, it's an error */

    /* Nth logical record (from the beginning) maps to Mth physical record
       in the file, where:

       m = (n+s) % N

       s = physical record that is 0th logical record
       N = total number of records in the db.  */

    physical_record = logical_record + header->start_record;

    if (header->num_records > 0) {
        physical_record = (physical_record % header->num_records);
    }

    return physical_record;
}

static int64_t _seek_to_logical_record(cdb_t *cdb, int64_t logical_record) {

    uint64_t physical_record = _physical_record_for_logical_record(cdb->header, logical_record);
    uint64_t offset = HEADER_SIZE + (physical_record * RECORD_SIZE);

    if (lseek(cdb->fd, offset, SEEK_SET) != offset) {
        return -1;
    }

    return physical_record;
}

static cdb_time_t _time_for_logical_record(cdb_t *cdb, int64_t logical_record) {

    cdb_record_t record[RECORD_SIZE];
    cdb_time_t time = 0;

    /* skip over any record that has NULL time or bad time values
       Such datapoints in cdb indicate a corrupted cdb. */
    while (!time || time <= 0) {

        if (_seek_to_logical_record(cdb, logical_record) < 0) {
            time = 0;
            break;
        }

        logical_record += 1;

        if (read(cdb->fd, &record, RECORD_SIZE) < 0) {
            time = 0;
            break;
        }

        time = record->time;
    }

    return time;
}

/* note - if no exact match, will return a record with a time greater than the requested value */
static int64_t _logical_record_for_time(cdb_t *cdb, cdb_time_t req_time, int64_t start_logical_record, int64_t end_logical_record) {

    bool first_time = false;
    cdb_time_t start_time, next_time, center_time;
    int64_t delta, center_logical_record, next_logical_record;

    uint64_t num_recs = cdb->header->num_records;

    /* for the very first time find start and end */
    if (start_logical_record == 0 && end_logical_record == 0) {

        start_logical_record = 0;
        end_logical_record   = num_recs - 1;

        first_time = true;
    }

    /* if no particular time was requested, just return the first one. */
    if (req_time == 0) {
        return start_logical_record;
    }

    /* if there are only 2 or 1 records in the search space, return the last record */
    if (end_logical_record - start_logical_record <= 1) {
        return end_logical_record;
    }

    start_time = _time_for_logical_record(cdb, start_logical_record);

    /* if requested time is less than start_time, start_time is the best we can do. */
    if (req_time <= start_time) {
        return start_logical_record;
    }

    /* don't go out of bounds */
    if (start_logical_record + 1 >= num_recs) {
        return start_logical_record;
    }

    next_logical_record = start_logical_record;
    next_time = start_time;

    /* if _seek_to_logical_record encounters bad data it fabricates and
       returns next good value. Try to get the *real* next rec. */
    while ((next_time - start_time) == 0) {

        next_logical_record += 1;

        next_time = _time_for_logical_record(cdb, next_logical_record);

        if (next_logical_record >= num_recs) {
            break;
        }
    }

    delta = next_time - start_time;

    /* delta = 0 means that _seek_to_logical_record fabricated data on the fly */
    if (delta == 0) {
        return start_logical_record;
    /* if we have wrapped over, or if the requested time is in the range of start and next, return start */
    } else if (delta < 0) {
        return start_logical_record;
    } else if (req_time <= next_time) {
        return next_logical_record;
    }

    /* for the very first binary search pivot point, take an aggressive guess to make search converge faster */
    if (first_time) {
        center_logical_record = (req_time - start_time) / (next_time - start_time) - 1;
    } else {
        center_logical_record = start_logical_record + (end_logical_record - start_logical_record) / 2;
    }

    if (num_recs > 0) {
        center_logical_record = (center_logical_record % num_recs);
    }

    center_time = _time_for_logical_record(cdb, center_logical_record);

    if (req_time >= center_time) {
        start_logical_record = center_logical_record;
    } else {
        end_logical_record = center_logical_record;
    }

    return _logical_record_for_time(cdb, req_time, start_logical_record, end_logical_record);
}

bool _cdb_is_writable(cdb_t *cdb) {

    /* We can't check for the O_RDONLY bit being because its defined value is zero.
     * i.e. there are no bits set to look for. We therefore assume
     * O_RDONLY if neither O_WRONLY nor O_RDWR are set. */
    if (cdb->flags & O_RDWR) {
        return true;
    }

    return false;
}

int cdb_read_header(cdb_t *cdb) {
    struct stat st;

    /* If the header has already been read from backing store do not read again */
    if (cdb->synced == true) {
        return CDB_SUCCESS;
    }

    if (cdb_open(cdb) != 0) {
        return cdb_error();
    }

    if (pread(cdb->fd, cdb->header, HEADER_SIZE, 0) != HEADER_SIZE) {
        return cdb_error();
    }

    if (strncmp(cdb->header->token, CDB_TOKEN, sizeof(CDB_TOKEN)) != 0) {
        return CDB_EBADTOK;
    }

    if (strncmp(cdb->header->version, CDB_VERSION, sizeof(CDB_VERSION)) != 0) {
        return CDB_EBADVER;
    }

    cdb->synced = true;

    /* Calculate the number of records */
    if (fstat(cdb->fd, &st) == 0) {
        cdb->header->num_records = (st.st_size - HEADER_SIZE) / RECORD_SIZE;
    } else {
        cdb->header->num_records = 0;
    }

    return CDB_SUCCESS;
}

int cdb_write_header(cdb_t *cdb) {

    if (cdb->synced) {
        return CDB_SUCCESS;
    }

    if (_cdb_is_writable(cdb) == false) {
        return CDB_ERDONLY;
    }

    cdb->synced = false;

    if (cdb_open(cdb) > 0) {
        return cdb_error();
    }

    if (pwrite(cdb->fd, cdb->header, HEADER_SIZE, 0) != HEADER_SIZE) {
        return cdb_error();
    }

    cdb->synced = true;

    return CDB_SUCCESS;
}

void cdb_print_header(cdb_t * cdb) {

    printf("version: [%s]\n", cdb->header->version);
    printf("name: [%s]\n", cdb->header->name);
    printf("desc: [%s]\n", cdb->header->desc);
    printf("units: [%s]\n", cdb->header->units);

    if (cdb->header->type == CDB_TYPE_COUNTER) {
        printf("type: [COUNTER]\n");
    }

    if (cdb->header->type == CDB_TYPE_GAUGE) {
        printf("type: [GAUGE]\n");
    }

    printf("min_value: [%g]\n", cdb->header->min_value);
    printf("max_value: [%g]\n", cdb->header->max_value);
    printf("max_records: [%"PRIu64"]\n", cdb->header->max_records);
    printf("num_records: [%"PRIu64"]\n", cdb->header->num_records);
    printf("start_record: [%"PRIu64"]\n", cdb->header->start_record);
}

int cdb_write_records(cdb_t *cdb, cdb_record_t *records, uint64_t len, uint64_t *num_recs) {

    /* read old header if it exists.
       write a header out, since the db may not have existed.
    */
    uint64_t i   = 0;
    uint64_t j   = 0;
    off_t offset = 0;
    *num_recs    = 0;

    if (cdb_read_header(cdb) != CDB_SUCCESS) {
        return cdb_error();
    }

    if (_cdb_is_writable(cdb) == false) {
        return CDB_ERDONLY;
    }

    if (cdb->header->max_records <= 0) {
        return CDB_EINVMAX;
    }

    /* Logic for writes:
    cdb is 5 records.
        try to write 7 records
        len = 7
        len + 0 > 5
        j = (7 + 0) - 5
        j = 2
        i = 7 - 2
        i = 5

        need to write j records at offset: HEADER_SIZE + (cdb->header->start_record * RECORD_SIZE)
        index into records is i

        need to write i records at offset: HEADER_SIZE + (cdb->header->num_records * RECORD_SIZE)
        index into records is 0;
    */

    /* Calculate our indicies into the records array. */
    if (len + cdb->header->num_records >= cdb->header->max_records) {
        j = (len + cdb->header->num_records) - cdb->header->max_records;
        i = (len - j);
    } else {
        i = len;
    }

    /* If we need to wrap around */
    if (j > 0) {

        offset = HEADER_SIZE + (cdb->header->start_record * RECORD_SIZE);

        cdb->header->start_record += j;
        cdb->header->start_record %= cdb->header->max_records;

        if (pwrite(cdb->fd, &records[i], (RECORD_SIZE * j), offset) != (RECORD_SIZE * j)) {
            return cdb_error();
        }

        cdb->synced = false;

        /* start_record is no longer 0, so update the header */
        if (cdb_write_header(cdb) != CDB_SUCCESS) {
            return cdb_error();
        }
    }

    /* Normal case */
    offset = HEADER_SIZE + (cdb->header->num_records * RECORD_SIZE);
    cdb->header->num_records += i;

    if (pwrite(cdb->fd, &records[0], (RECORD_SIZE * i), offset) != (RECORD_SIZE * i)) {
        return cdb_error();
    }

    *num_recs += len;

#ifdef DEBUG
    printf("write_records: wrote [%"PRIu64"] records\n", *num_recs);
#endif

    return CDB_SUCCESS;
}

bool cdb_write_record(cdb_t *cdb, cdb_time_t time, double value) {

    cdb_record_t record[RECORD_SIZE];
    uint64_t num_recs = 0;

    record->time  = time;
    record->value = value;

    if (cdb_write_records(cdb, record, 1, &num_recs) != CDB_SUCCESS) {
        return false;
    }

    return true;
}

int cdb_update_records(cdb_t *cdb, cdb_record_t *records, uint64_t len, uint64_t *num_recs) {

    int ret    = CDB_SUCCESS;
    *num_recs  = 0;
    uint64_t i = 0;

    if (cdb_read_header(cdb) != CDB_SUCCESS) {
        return cdb_error();
    }

    if (_cdb_is_writable(cdb) == false) {
        return CDB_ERDONLY;
    }

#ifdef DEBUG
    printf("in update_records with [%"PRIu64"] num_recs\n", cdb->header->num_records);
#endif

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

        cdb_time_t time = records[i].time;
        cdb_time_t rtime;
        uint64_t lrec;

        lrec = _logical_record_for_time(cdb, time, 0, 0);

        if (lrec >= 1) {
            lrec -= 1;
            /* DRK things like this are extremely confusing. some logical
               functions are zero based? some are one? or you're trying to back up? */
        }

        rtime = _time_for_logical_record(cdb, lrec);

        while (rtime < time && lrec < cdb->header->num_records - 1) {

            /* DRK is this true? or is the condition (lrec - start_lrec) <
               num_recs -- seems you can't wrap here (i.e. what if the initial
               lrec was num_recs - 1) */

            lrec += 1;
            rtime = _time_for_logical_record(cdb, lrec);
        }

#ifdef DEBUG
        printf("update_records: value: lrec [%"PRIu64"] time [%d] rtime [%d]\n", lrec, (int)time, (int)rtime);
#endif

        while (time == rtime && lrec < cdb->header->num_records - 1) {

            _seek_to_logical_record(cdb, lrec);

            if (write(cdb->fd, &records[0], RECORD_SIZE) != RECORD_SIZE) {
                ret = cdb_error();
                break;
            }

            lrec += 1;

            rtime = _time_for_logical_record(cdb, lrec);
        }
    }

    if (ret == CDB_SUCCESS) {

        if (i > 0) {
            cdb->synced = false;
            *num_recs = i;
        }

        if (cdb_write_header(cdb) != CDB_SUCCESS) {
            ret = cdb_error();
        }
    }

    return ret;
}

bool cdb_update_record(cdb_t *cdb, cdb_time_t time, double value) {

    cdb_record_t record[RECORD_SIZE];
    uint64_t num_recs = 0;

    record->time  = time;
    record->value = value;

    if (cdb_update_records(cdb, record, 1, &num_recs) != 0) {
        return false;
    }

    return true;
}

int cdb_discard_records_in_time_range(cdb_t *cdb, cdb_request_t *request, uint64_t *num_recs) {

    uint64_t i = 0;
    int64_t lrec;
    off_t offset = RECORD_SIZE;
    *num_recs = 0;

    if (cdb_read_header(cdb) != CDB_SUCCESS) {
        return cdb_error();
    }

    if (_cdb_is_writable(cdb) == false) {
        return CDB_ERDONLY;
    }

    lrec = _logical_record_for_time(cdb, request->start, 0, 0);

    if (lrec >= 1) {
        lrec -= 1;
    }

    for (i = lrec; i < cdb->header->num_records; i++) {

        cdb_time_t rtime = _time_for_logical_record(cdb, i);

        if (rtime >= request->start && rtime <= request->end) {

            cdb_record_t record[RECORD_SIZE];

            record->time  = rtime;
            record->value = CDB_NAN;

            if (pwrite(cdb->fd, &record, RECORD_SIZE, offset) != RECORD_SIZE) {
                return cdb_error();
            }

            *num_recs += 1;
        }
    }

    if (*num_recs > 0) {
        cdb->synced = false;
    }

    if (cdb_write_header(cdb) != CDB_SUCCESS) {
        return cdb_error();
    }

    return CDB_SUCCESS;
}

static int _compute_scale_factor_and_num_records(cdb_t *cdb, int64_t *num_records, int32_t *factor) {

    if (cdb->header->type == CDB_TYPE_COUNTER) {
        if (*num_records != 0) {

            if (*num_records > 0) {
                *num_records += 1;
            } else {
                *num_records -= 1;
            }
        }
    }

    if (strlen(cdb->header->units) > 0) {

        int32_t multiplier = 1;
        char *frequency;

        if ((frequency = calloc(strlen(cdb->header->units), sizeof(char))) == NULL) {
            return CDB_ENOMEM;
        }

        if ((sscanf(cdb->header->units, "per %d %s", &multiplier, frequency) == 2) ||
            (sscanf(cdb->header->units, "per %s", frequency) == 1) ||
            (sscanf(cdb->header->units, "%*s per %s", frequency) == 1)) {

            if (strcmp(frequency, "min") == 0) {
                *factor = 60;
            } else if (strcmp(frequency, "hour") == 0) {
                *factor = 60 * 60;
            } else if (strcmp(frequency, "sec") == 0 || strcmp(frequency, "second") == 0) {
                *factor = 1;
            } else if (strcmp(frequency, "day") == 0) {
                *factor = 60 * 60 * 24;
            } else if (strcmp(frequency, "week") == 0) {
                *factor = 60 * 60 * 24 * 7;
            } else if (strcmp(frequency, "month") == 0) {
                *factor = 60 * 60 * 24 * 30;
            } else if (strcmp(frequency, "quarter") == 0) {
                *factor = 60 * 60 * 24 * 90;
            } else if (strcmp(frequency, "year") == 0) {
                *factor = 60 * 60 * 24 * 365;
            }

            if (*factor != 0) {
                *factor *= multiplier;
            }
        }

        free(frequency);
    }

    return CDB_SUCCESS;
}

/* Statistics code
 * Make only one call to reading for a particular time range and compute all our stats
 */
void _compute_statistics(cdb_range_t *range, uint64_t *num_recs, cdb_record_t *records) {

    uint64_t i     = 0;
    uint64_t valid = 0;
    double sum     = 0.0;
    double *values = calloc(*num_recs, sizeof(double));

    for (i = 0; i < *num_recs; i++) {

        if (!isnan(records[i].value)) {

            sum += values[valid] = records[i].value;
            valid++;
        }
    }

    range->num_recs = valid;
    range->mean     = gsl_stats_mean(values, 1, valid);
    range->max      = gsl_stats_max(values, 1, valid);
    range->min      = gsl_stats_min(values, 1, valid);
    range->sum      = sum;
    range->stddev   = gsl_stats_sd(values, 1, valid);
    range->absdev   = gsl_stats_absdev(values, 1, valid);

    /* The rest need sorted data */
    gsl_sort(values, 1, valid);

    range->median   = gsl_stats_median_from_sorted_data(values, 1, valid);
    range->pct95th  = gsl_stats_quantile_from_sorted_data(values, 1, valid, 0.95);
    range->pct75th  = gsl_stats_quantile_from_sorted_data(values, 1, valid, 0.75);
    range->pct50th  = gsl_stats_quantile_from_sorted_data(values, 1, valid, 0.50);
    range->pct25th  = gsl_stats_quantile_from_sorted_data(values, 1, valid, 0.25);

    /* MAD must come last because it alters the values array
     * http://en.wikipedia.org/wiki/Median_absolute_deviation */
    for (i = 0; i < valid; i++) {
        values[i] = fabs(values[i] - range->median);

        if (values[i] < 0.0) {
            values[i] *= -1.0;
        }
    }

    /* Final sort is required MAD */
    gsl_sort(values, 1, valid);
    range->mad = gsl_stats_median_from_sorted_data(values, 1, valid);

    free(values);
}

double cdb_get_statistic(cdb_range_t *range, cdb_statistics_enum_t type) {

    switch (type) {
        case CDB_MEDIAN:
            return range->median;
        case CDB_MAD:
            return range->mad;
        case CDB_95TH:
            return range->pct95th;
        case CDB_75TH:
            return range->pct75th;
        case CDB_50TH:
            return range->pct50th;
        case CDB_25TH:
            return range->pct25th;
        case CDB_MEAN:
            return range->mean;
        case CDB_SUM:
            return range->sum;
        case CDB_MAX:
            return range->max;
        case CDB_MIN:
            return range->min;
        case CDB_STDDEV:
            return range->stddev;
        case CDB_ABSDEV:
            return range->absdev;
        default:
            fprintf(stderr, "aggregate_using_function_for_records() function: [%d] not supported\n", type);
            return CDB_FAILURE;
    }
}

static int _cdb_read_records(cdb_t *cdb, cdb_request_t *request, uint64_t *num_recs, cdb_record_t **records) {

    int64_t first_requested_logical_record;
    int64_t last_requested_logical_record;
    uint64_t last_requested_physical_record;
    int64_t seek_physical_record;

    cdb_record_t *buffer = NULL;

    if (cdb_read_header(cdb) != CDB_SUCCESS) {
        return cdb_error();
    }

    if (request->start != 0 && request->end != 0 && request->end < request->start) {
        return CDB_ETMRANGE;
    }

    if (cdb->header == NULL || cdb->synced == false) {
        return CDB_ESANITY;
    }

    /* bail out if there are no records */
    if (cdb->header->num_records <= 0) {
        return CDB_ENORECS;
    }

    /*
      get the number of requested records:

      -ve indicates n records from the beginning
      +ve indicates n records off of the end.
      0 or undef means the whole thing.

      switch the meaning of -ve/+ve to be more array like
    */
    if (request->count != 0) {
        request->count = -request->count;
    }

#ifdef DEBUG
    printf("read_records start: [%ld]\n", request->start);
    printf("read_records end: [%ld]\n", request->end);
    printf("read_records num_requested: [%"PRIu64"]\n", request->count);
#endif

    if (request->count != 0 && request->count < 0 && request->start == 0) {
        /* if reading only few records from the end, just set -ve offset to seek to */
        first_requested_logical_record = request->count;

    } else {
        /* compute which record to start reading from the beginning, based on start time specified. */
        first_requested_logical_record = _logical_record_for_time(cdb, request->start, 0, 0);
    }

    /* if end is not defined, read all the records or only read uptill the specified record. */
    if (request->end == 0) {

        last_requested_logical_record = cdb->header->num_records - 1;

    } else {

        last_requested_logical_record = _logical_record_for_time(cdb, request->end, 0, 0);

        /* this can return something > end, check for that */
        if (_time_for_logical_record(cdb, last_requested_logical_record) > request->end) {
            last_requested_logical_record -= 1;
        }
    }

    last_requested_physical_record = (last_requested_logical_record + cdb->header->start_record) % cdb->header->num_records;

    /* After _seek_to_logical_record(), we're at the offset to read from. */
    seek_physical_record = _seek_to_logical_record(cdb, first_requested_logical_record);

    if (last_requested_physical_record >= seek_physical_record) {

        uint64_t nrec = (last_requested_physical_record - seek_physical_record + 1);
        uint64_t rlen = RECORD_SIZE * nrec;

        if ((buffer = calloc(1, rlen)) == NULL) {
            free(buffer);
            return CDB_ENOMEM;
        }

        if (read(cdb->fd, buffer, rlen) != rlen) {
            free(buffer);
            return cdb_error();
        }

        *num_recs = nrec;

    } else {

        /* We've wrapped around the end of the file */
        uint64_t nrec1 = (cdb->header->num_records - seek_physical_record);
        uint64_t nrec2 = (last_requested_physical_record + 1);

        uint64_t rlen1 = RECORD_SIZE * nrec1;
        uint64_t rlen2 = RECORD_SIZE * nrec2;

        if ((buffer = calloc(1, rlen1 + rlen2)) == NULL) {
            free(buffer);
            return CDB_ENOMEM;
        }

        /* Read at the offset set by _seek_to_logical_record() */
        if (read(cdb->fd, buffer, rlen1) != rlen1) {
            free(buffer);
            return cdb_error();
        }

        /* And then the wrap around portion past the header. */
        if (pread(cdb->fd, &buffer[nrec1], rlen2, HEADER_SIZE) != rlen2) {
            free(buffer);
            return cdb_error();
        }

        *num_recs = nrec1 + nrec2;
    }

    /* Deal with cooking the output */
    if (request->cooked) {

        bool check_min_max = true;
        int32_t factor = 0;
        uint64_t i = 0;
        uint64_t cooked_recs = 0;
        double prev_value = 0.0;
        cdb_time_t prev_date = 0;
        cdb_record_t *crecords;

        if ((crecords = calloc(*num_recs, RECORD_SIZE)) == NULL) {
            free(crecords);
            free(buffer);
            return CDB_ENOMEM;
        }

        if (_compute_scale_factor_and_num_records(cdb, &request->count, &factor)) {
            free(crecords);
            free(buffer);
            return cdb_error();
        }

        if (cdb->header->min_value == 0 && cdb->header->max_value == 0) {
            check_min_max = false;
        }

        for (i = 0; i < *num_recs; i++) {

            cdb_time_t date = buffer[i].time;
            double value    = buffer[i].value;

            if (cdb->header->type == CDB_TYPE_COUNTER) {
                double new_value = value;
                value = CDB_NAN;

                if (!isnan(prev_value) && !isnan(new_value)) {

                    double val_delta = new_value - prev_value;

                    if (val_delta >= 0) {
                        value = val_delta;
                    }
                }

                prev_value = new_value;
            }

            if (factor != 0 && cdb->header->type == CDB_TYPE_COUNTER) {

                /* Skip the first entry, since it's absolute and is needed
                 * to calculate the second */
                if (prev_date == 0) {
                    prev_date = date;
                    continue;
                }

                cdb_time_t time_delta = date - prev_date;

                if (time_delta > 0 && !isnan(value)) {
                    value = factor * (value / time_delta);
                }

                prev_date = date;
            }

            /* Check for min/max boundaries */
            /* Should this be done on write instead of read? */
            if (check_min_max && !isnan(value)) {
                if (value > cdb->header->max_value || value < cdb->header->min_value) {
                    value = CDB_NAN;
                }
            }

            /* Copy the munged data to our new array, since we might skip
             * elements. Also keep in mind mmap for the future */
            crecords[cooked_recs].time  = date;
            crecords[cooked_recs].value = value;
            cooked_recs += 1;
        }

        /* Now swap our cooked records for the buffer, so we can slice it as needed. */
        free(buffer);
        buffer = crecords;
        *num_recs = cooked_recs;
    }

    /* If we've been requested to average the records & timestamps */
    if (request->step > 1) {

        cdb_record_t *arecords;
        uint32_t step      = request->step;
        uint64_t step_recs = 0;
        uint64_t leftover  = (*num_recs % step);
        uint64_t walkend   = (*num_recs - leftover);
        uint64_t i = 0;

        if ((arecords = calloc(((*num_recs / step) + leftover), RECORD_SIZE)) == NULL) {
            free(arecords);
            free(buffer);
            return CDB_ENOMEM;
        }

        /* Walk our list of cooked records, jumping ahead by the given step.
           For each set of records within that step, we want to get the average
           for those records and place them into a new array.
         */
        for (i = 0; i < walkend; i += step) {

            uint64_t j = 0;
            double xi[step];
            double yi[step];

            for (j = 0; j < step; j++) {

                /* No NaNs on average - they cause bogus graphs. Is there a
                 * better value than 0 to use here? */
                if (isnan(buffer[i+j].value)) {
                    buffer[i+j].value = 0;
                }

                xi[j] = (double)buffer[i+j].time;
                yi[j] = buffer[i+j].value;
            }

            arecords[step_recs].time  = (cdb_time_t)gsl_stats_mean(xi, 1, step);
            arecords[step_recs].value = gsl_stats_mean(yi, 1, step);
            step_recs += 1;
        }

        /* Now collect from the last step point to the end & average. */
        if (leftover > 0) {
            uint64_t leftover_start = *num_recs - leftover;

            uint64_t j = 0;
            double xi[leftover];
            double yi[leftover];

            for (i = leftover_start; i < *num_recs; i++) {

                /* No NaNs on average - they cause bogus graphs. Is there a
                 * better value than 0 to use here? */
                if (isnan(buffer[i].value)) {
                    buffer[i].value = 0;
                }

                xi[j] = (double)buffer[i].time;
                yi[j] = buffer[i].value;
                j++;
            }

            arecords[step_recs].time  = (cdb_time_t)gsl_stats_mean(xi, 1, j);
            arecords[step_recs].value = gsl_stats_mean(yi, 1, j);
            step_recs += 1;
        }

        free(buffer);
        buffer = arecords;
        *num_recs = step_recs;
    }

    /* now pull out the number of requested records if asked */
    if (request->count != 0 && *num_recs >= abs(request->count)) {

        uint64_t start_index = 0;

        if (request->count <= 0) {
            start_index = *num_recs - abs(request->count);
        }

        *num_recs = abs(request->count);

        if ((*records  = calloc(*num_recs, RECORD_SIZE)) == NULL) {
            free(buffer);
            return CDB_ENOMEM;
        }

        memcpy(*records, &buffer[start_index], RECORD_SIZE * *num_recs);

        free(buffer);

    } else {

        *records = buffer;
    }

    return CDB_SUCCESS;
}

int cdb_read_records(cdb_t *cdb, cdb_request_t *request,
    uint64_t *num_recs, cdb_record_t **records, cdb_range_t *range) {

    int ret   = CDB_SUCCESS;

    ret = _cdb_read_records(cdb, request, num_recs, records);

    if (ret == CDB_SUCCESS) {

        if (*num_recs > 0) {
            range->start_time = request->start;
            range->end_time   = request->end;

            _compute_statistics(range, num_recs, *records);
        }
    }

    return ret;
}

void cdb_print_records(cdb_t *cdb, cdb_request_t *request, FILE *fh, const char *date_format) {

    uint64_t i = 0;
    uint64_t num_recs = 0;

    cdb_record_t *records = NULL;

    if (_cdb_read_records(cdb, request, &num_recs, &records) == CDB_SUCCESS) {

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

            _print_record(fh, records[i].time, records[i].value, date_format);
        }
    }

    free(records);
}

void cdb_print(cdb_t *cdb) {

    const char *date_format = "%Y-%m-%d %H:%M:%S";
    cdb_request_t request;

    request.start  = 0;
    request.end    = 0;
    request.count  = 0;
    request.step   = 0;
    request.cooked = false;

    printf("============== Header ================\n");

    if (cdb_read_header(cdb) == CDB_SUCCESS) {
        cdb_print_header(cdb);
    }

    if (cdb->header->type == CDB_TYPE_COUNTER) {

        printf("============= Raw Counter Records =============\n");
        cdb_print_records(cdb, &request, stdout, date_format);
        printf("============== End Raw Counter Records ==============\n");

        printf("============== Cooked Records ================\n");

    } else {
        printf("============== Records ================\n");
    }

    request.cooked = true;
    cdb_print_records(cdb, &request, stdout, date_format);

    printf("============== End ================\n");
}

/* Take in an array of cdbs */
int cdb_read_aggregate_records(cdb_t **cdbs, int num_cdbs, cdb_request_t *request,
    uint64_t *driver_num_recs, cdb_record_t **records, cdb_range_t *range) {

    uint64_t i = 0;
    int ret    = CDB_SUCCESS;
    cdb_record_t *driver_records = NULL;
    *driver_num_recs = 0;

    if (cdbs[0] == NULL) {
        return CDB_ESANITY;
    }

    /* The first cdb is the driver */
    ret = _cdb_read_records(cdbs[0], request, driver_num_recs, &driver_records);

    if (ret != CDB_SUCCESS) {
        fprintf(stderr, "Bailed on: %s\n", cdbs[0]->filename);
        free(driver_records);
        return ret;
    }

    if ((*records = calloc(*driver_num_recs, RECORD_SIZE)) == NULL) {
        free(driver_records);
        return CDB_ENOMEM;
    }

    if (*driver_num_recs <= 1) {
        free(driver_records);
        return CDB_EINTERPD;
    }

    double *driver_x_values   = calloc(*driver_num_recs, sizeof(double));
    double *driver_y_values   = calloc(*driver_num_recs, sizeof(double));
    double *follower_x_values = calloc(*driver_num_recs, sizeof(double));
    double *follower_y_values = calloc(*driver_num_recs, sizeof(double));

    for (i = 0; i < *driver_num_recs; i++) {
        (*records)[i].time  = driver_x_values[i] = driver_records[i].time;
        (*records)[i].value = driver_y_values[i] = driver_records[i].value;
    }

    /* initialize and allocate the gsl objects  */
    gsl_interp_accel *accel = gsl_interp_accel_alloc();
    gsl_interp *interp = gsl_interp_alloc(gsl_interp_linear, *driver_num_recs);

    /* Allows 0.0 to be returned as a valid yi */
    gsl_set_error_handler_off();

    gsl_interp_init(interp, driver_x_values, driver_y_values, *driver_num_recs);

    for (i = 1; i < num_cdbs; i++) {

        uint64_t j = 0;
        uint64_t follower_num_recs = 0;

        cdb_record_t *follower_records = NULL;

        ret = _cdb_read_records(cdbs[i], request, &follower_num_recs, &follower_records);

        /* Just bail, free all allocations below and let the error bubble up */
        if (ret == CDB_SUCCESS && follower_num_recs != 0) {

            for (j = 0; j < *driver_num_recs; j++) {

                /* Check for out of bounds */
                if (j >= follower_num_recs) {
                    break;
                }

                follower_x_values[j] = follower_records[j].time;
                follower_y_values[j] = follower_records[j].value;
            }

            for (j = 0; j < *driver_num_recs; j++) {

                double yi = gsl_interp_eval(interp, follower_x_values, follower_y_values, driver_x_values[j], accel);

                if (isnormal(yi)) {
                    (*records)[j].value += yi;
                }
            }
        }

        free(follower_records);

        if (ret != CDB_SUCCESS) {
            break;
        }
    }

    if (ret == CDB_SUCCESS && *driver_num_recs > 0) {
        /* Compute all the statistics for this range */
        range->start_time = request->start;
        range->end_time   = request->end;

        _compute_statistics(range, driver_num_recs, *records);
    }

    free(driver_x_values);
    free(driver_y_values);
    free(follower_x_values);
    free(follower_y_values);

    gsl_interp_free(interp);
    gsl_interp_accel_free(accel);
    free(driver_records);

    return ret;
}

void cdb_print_aggregate_records(cdb_t **cdbs, int32_t num_cdbs, cdb_request_t *request, FILE *fh, const char *date_format) {

    uint64_t i = 0;
    uint64_t num_recs = 0;

    cdb_record_t *records = NULL;
    cdb_range_t *range = calloc(1, sizeof(cdb_range_t));

    cdb_read_aggregate_records(cdbs, num_cdbs, request, &num_recs, &records, range);

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

        _print_record(fh, records[i].time, records[i].value, date_format);
    }

    free(range);
    free(records);
}

void cdb_generate_header(cdb_t *cdb, char* name, char* desc, uint64_t max_records, int32_t type,
    char* units, uint64_t min_value, uint64_t max_value) {

    if (max_records == 0) {
        max_records = CDB_DEFAULT_RECORDS;
    }

    if (type == 0) {
        cdb->header->type = CDB_DEFAULT_DATA_TYPE;
    }

    if (units == NULL || (strcmp(units, "") == 0)) {
        units = (char*)CDB_DEFAULT_DATA_UNIT;
    }

    if (desc == NULL) {
        desc = (char*)"";
    }

    memset(cdb->header->name, 0, sizeof(cdb->header->name));
    memset(cdb->header->desc, 0, sizeof(cdb->header->desc));
    memset(cdb->header->units, 0, sizeof(cdb->header->units));
    memset(cdb->header->version, 0, sizeof(cdb->header->version));
    memset(cdb->header->token, 0, sizeof(cdb->header->token));

    strncpy(cdb->header->name, name, sizeof(cdb->header->name));
    strncpy(cdb->header->desc, desc, sizeof(cdb->header->desc));
    strncpy(cdb->header->units, units, sizeof(cdb->header->units));
    strncpy(cdb->header->version, CDB_VERSION, sizeof(cdb->header->version));
    strncpy(cdb->header->token, CDB_TOKEN, sizeof(cdb->header->token));

    cdb->header->type         = type;
    cdb->header->max_records  = max_records;
    cdb->header->min_value    = min_value;
    cdb->header->max_value    = max_value;
    cdb->header->num_records  = 0;
    cdb->header->start_record = 0;
}

cdb_t* cdb_new(void) {

    cdb_t *cdb  = calloc(1, sizeof(cdb_t));
    cdb->header = calloc(1, HEADER_SIZE);

    cdb->fd = -1;
    cdb->synced = false;
    cdb->flags = -1;
    cdb->mode = -1;

    return cdb;
}

cdb_request_t cdb_new_request(void) {
    cdb_request_t request;
    memset (&request, 0, sizeof (request));
    request.start  = 0;
    request.end    = 0;
    request.count  = 0;
    request.cooked = true;
    request.step   = 0;
    return request;
}

int cdb_open(cdb_t *cdb) {

    if (cdb->fd >= 0) {
        return CDB_SUCCESS;
    }

    /* Default flags if none were set */
    if (cdb->flags == -1) {
        cdb->flags = O_RDONLY|O_BINARY;
    }

    /* A cdb can't be write only - we need to read the header */
    if (cdb->flags & O_WRONLY) {
        cdb->flags = O_RDWR;
    }

    if (cdb->mode == -1) {
        cdb->mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH;
    }

    cdb->fd = open(cdb->filename, cdb->flags, cdb->mode);

    if (cdb->fd < 0) {
        return cdb_error();
    }

    return CDB_SUCCESS;
}

int cdb_close(cdb_t *cdb) {

    if (cdb != NULL) {

        if (cdb->fd > 0) {
            if (close(cdb->fd) != 0) {
                return cdb_error();
            }
            cdb->fd = -1;
        }
    }

    return CDB_SUCCESS;
}

int cdb_free(cdb_t *cdb) {
    int ret = CDB_SUCCESS;

    if (cdb != NULL) {

        if (cdb->header != NULL) {
            ret = cdb_close(cdb);
            free(cdb->header);
            cdb->header = NULL;
        }

        free(cdb);
        cdb = NULL;
    }

    return ret;
}

/* -*- Mode: C; tab-width: 4 -*- */
/* vim: set tabstop=4 expandtab shiftwidth=4: */

max_stars_repo_name: dsully/circulardb

max_issues_repo_name: dsully/circulardb

max_forks_repo_name: dsully/circulardb


text:
/*
 * Gsl_blas_dgemm_matrix_view: test gsl_blas_dgemm (matrix . matrix)
 *
 * Copyright (c) 2012 Jérémie Decock
 *
 * Required: GSL library (libgsl0-dev)
 * Usage: gcc gsl_blas_dgemm_matrix_view.c -lgsl -lgslcblas -lm
 *    or: gcc gsl_blas_dgemm_matrix_view.c $(pkg-config --libs gsl)
 */

#include <stdio.h>
#include <gsl/gsl_blas.h>

int main (void)
{
    double a[] = { 0.11, 0.12, 0.13,
                   0.21, 0.22, 0.23 };

    double b[] = { 1011, 1012,
                   1021, 1022,
                   1031, 1032 };

    double c[] = { 0.00, 0.00,
                   0.00, 0.00 };

    gsl_matrix_view A = gsl_matrix_view_array(a, 2, 3);
    gsl_matrix_view B = gsl_matrix_view_array(b, 3, 2);
    gsl_matrix_view C = gsl_matrix_view_array(c, 2, 2);

    /* Compute C = A B */

    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
                   1.0, &A.matrix, &B.matrix,
                   0.0, &C.matrix);

    printf("[ %g, %g\n", c[0], c[1]);
    printf("  %g, %g ]\n", c[2], c[3]);

    return 0;
}

max_stars_repo_name: jeremiedecock/snippets

max_issues_repo_name: jeremiedecock/snippets

max_forks_repo_name: jeremiedecock/snippets


text:
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_chebyshev.h>

double f1(double x, void *p) {
    (void) (p); /* avoid unused parameter warning */
    return exp(pow(x, 2));
}

double f2(double x, void *p) {
    (void) (p); /* avoid unused parameter warning */
    return fabs(x + pow(x, 3));
}

double f3(double x, void *p) {
    (void) (p); /* avoid unused parameter warning */
    return x <= 0 ? -1 : 1;
}

double (*f[])(double, void *) = {f1, f2, f3};

char *a[] = {"pierwsze1.dat", "pierwsze2.dat", "pierwsze3.dat"};

int main(void) {
    int i, n = 10000;
    gsl_cheb_series *cs = gsl_cheb_alloc(40);
    gsl_function F;
    FILE *file;
    for (int j = 0; j < 3; j++) {
        file = fopen(a[j], "w");
        F.function = f[j];
        F.params = 0;
        gsl_cheb_init(cs, &F, -1.0, 1.0);
        for (i = -n; i < n; i++) {
            double x = i / (double) n;
            double r10 = gsl_cheb_eval_n(cs, 10, x);
            double r40 = gsl_cheb_eval(cs, x);
            fprintf(file, "%g %g %g %g\n",
                    x, GSL_FN_EVAL (&F, x), r10, r40);
        }
        fclose(file);
    }
    gsl_cheb_free(cs);
    return EXIT_SUCCESS;
}

max_stars_repo_name: mistyfiky/agh-mownit

max_issues_repo_name: mistyfiky/agh-mownit

max_forks_repo_name: mistyfiky/agh-mownit


text:
/*
 * Copyright 2014 Marc Normandin
 *
 * 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.
 */

/*
 * pso.h
 *
 *  Created on: Jun 30, 2013
 *      Author: marc
 */

#ifndef PSO_H_
#define PSO_H_

#include <iostream>
#include <vector>
#include <algorithm>
#include <functional>
#include <cmath>
#include <cstdlib>
#include <gsl/gsl_rng.h>
#include <limits>
#include <fstream>
#include <cassert>

#include "rng.h"
#include "dim.h"
#include "particle.h"

template<class FitnessFunction>
class PSO
{
public:
    // This routine is used by the PSO unit test
    PSO(const unsigned int numParticles, const std::vector<Dim>& dim,
        const gslseed_t seed, FitnessFunction& fitnessFunction, const unsigned int maxIterations)
        : mNumParticles(numParticles), mGBest(dim), mDim(dim), mRng(seed), mFitnessFunction(fitnessFunction),
        mMaxIterations(maxIterations)
    {
        mParticles.reserve(mNumParticles);
    }

    unsigned int getNumParticles() const
    {
        return mNumParticles;
    }

    const Particle& iterate()
    {
        createRandomParticles();

        unsigned int numIterations = 0;
        std::vector<double> particleFitnesses(mParticles.size());
        do
        {
            // Evaluate the fitness/objective function
        	mFitnessFunction(mParticles, &particleFitnesses);

            //For each particle
            for (unsigned int i = 0; i < mParticles.size(); i++)
            {
                // Calculate fitness value
                const prob_t fitness = particleFitnesses[i]; //

                // If the fitness value is better than the best fitness value (pBest) in history
                // set current value as the new pBest
                mParticles[i].updateFitness( fitness );
            }

            // Choose the particle with the best fitness value of all the particles as the gBest
            std::vector<Particle>::const_iterator best = std::max_element(mParticles.begin(), mParticles.end());
            mGBest = Particle( best->getBestPosition(), best->getBestFitness() );

            // Update the inertia weight
            const double inertiaWeight = computeInertiaWeight(numIterations, mMaxIterations);

            // For each particle
            for (unsigned int i = 0; i < mParticles.size(); i++)
            {
                mParticles[i].updatePosition( mGBest, mDim, mCognitiveWeight, mSocialWeight, mRng, inertiaWeight );
            }

            numIterations++;
        }
        while(numIterations < mMaxIterations);

        return mGBest;
    }

protected:
    PSO(const PSO&);
    void operator=(const PSO&);

    void createRandomParticles() {
        mParticles.clear(); // Remove previous particles in the container

        // For each particle
        for (unsigned int i = 0; i < mNumParticles; i++)
        {
            std::vector<dim_t> pos;

            // For each dimension
            for (unsigned int d = 0; d < mDim.size(); d++)
            {
                dim_t posd = mRng.uniform( mDim[d].min(), mDim[d].max() );
                pos.push_back(posd);
            }

            // Add the particle to the collection
            mParticles.push_back( Particle(pos) );
        }
    }

    // Compute inertia. This is based on equation 4.1 from:
    // http://www.hindawi.com/journals/ddns/2010/462145/
    double computeInertiaWeight(const unsigned int iteration, const unsigned int maxInterations) const
    {
        // Note. We add +1 because our iterations go from 0 to max-1.
        double inertiaWeight = (mOmega1 - mOmega2) * ( (maxInterations - (iteration+1.0)) / (1.0 * (iteration+1.0) ) ) + mOmega2;
        return inertiaWeight;
    }

private:
    unsigned int            mNumParticles;
    std::vector<Particle> 	mParticles; // Particle positions
    Particle 				mGBest;
    std::vector<Dim>		mDim;

    // These values control how random the particle velocities are
    static const double 	mCognitiveWeight;
    static const double 	mSocialWeight;

    // These values control the rate of convergence
    static const double     mOmega1;
    static const double     mOmega2;

    RandomNumberGenerator 	mRng;

    FitnessFunction&        mFitnessFunction;
    unsigned int            mMaxIterations;
};

#endif /* PSO_H_ */

max_stars_repo_name: marcnormandin/ParticleSwarmOptimization

max_issues_repo_name: marcnormandin/ParticleSwarmOptimization

max_forks_repo_name: marcnormandin/ParticleSwarmOptimization


text:
/***
 * Copyright 2019 The Katla Authors
 *
 * 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.
 */

#ifndef KATLA_PIPE_H
#define KATLA_PIPE_H

#include "katla/core/core.h"

#include <gsl/span>

#include <optional>

namespace outcome = OUTCOME_V2_NAMESPACE;

namespace katla {

class PosixPipe {
public:
    PosixPipe();
    ~PosixPipe();

    outcome::result<void> open();

    outcome::result<ssize_t> read(gsl::span<std::byte>& buffer);
    outcome::result<ssize_t> write(gsl::span<std::byte>& buffer);

    outcome::result<void> redirectToRead(int fd_src);
    outcome::result<void> redirectToWrite(int fd_src);

    outcome::result<void> close();
    outcome::result<void> closeRead();
    outcome::result<void> closeWrite();
private:
    int _fd[2];
};

}

#endif // KATLA_PIPE_H

max_stars_repo_name: plok/katla

max_issues_repo_name: plok/katla

max_forks_repo_name: plok/katla