text:
/*
** Task-related edge density, main algorithm
**
** Ref:  Lohmann et al (2016) PLoS One, 11(6):e0158185
**
** G.Lohmann, MPI-KYB, Jan 2015
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_sort_float.h>

#include "viaio/Vlib.h"
#include "viaio/VImage.h"
#include "viaio/mu.h"
#include "viaio/option.h"

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

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

extern VImage *VImagePointer(VAttrList list,int *nt);
extern void   VReadImagePointer(VAttrList list,VImage *src);
extern VImage VoxelMap(VImage mask,size_t *nvoxels);
extern void   VDataMatrix(VImage *src,int first,int len,VImage map,gsl_matrix_float *X);
extern void   VCheckMatrix(gsl_matrix_float *X);
extern long   VMaskCoverage(VAttrList list,VImage mask);
extern float  ZMatrix(gsl_histogram *histogram,gsl_matrix_float *SNR1,gsl_matrix_float *SNR2,int n1,int n2,
		      VImage roi,VImage map,VImage,int,float elength,float quantile,int,int);
extern void   GetSNR(gsl_matrix_float **X1,gsl_matrix_float **X2,int *,int n,gsl_matrix_float *SNR,int);
extern void   GetMedian(gsl_matrix_float **X1,gsl_matrix_float **X2,int *,int n,gsl_matrix_float *SNR,int);
extern void   VPrintHistogram(gsl_histogram *histogram,int numperm,VString filename);
extern void   HistoUpdate(float *A,size_t nvox,gsl_histogram *hist);

extern size_t EdgeDensity(float *C,int *I,int *J,size_t nedges,
			  gsl_histogram *histogram,gsl_matrix_float *SNR1,gsl_matrix_float *SNR2,int n1,int n2,
			  VImage roi,VImage map,VImage,int,float,float,int,int,int,float,int);

extern size_t EstimateEdges(float *E,int *I,int *J,size_t fulldim,
			    gsl_histogram *TedHist,gsl_matrix_float *SNR1,gsl_matrix_float *SNR2,int n1,int n2,
			    VImage roi,VImage map,VImage mapimage,int adjdef,float elength,float zthreshold,float,int);

/* generate permutation table */
int genperm(gsl_rng *rx,int *table,int n,int numperm)
{
  int s=0;
  int ks = 0;
  int kc=3;
  if (numperm < 10) {
    if (n%2==0) kc = n/2;
    else kc = n/2-1;
  }
  int iter=0;
  while ((ks < kc || ks > n-kc)) {
    if (iter > 100) VError(" not possible to generate random permutations, perhaps not enough input data ?");
    ks=0;
    for (s=0; s<n; s++) {
      table[s] = 0;
      if (gsl_ran_bernoulli (rx,(double)0.5) == 1) {
	table[s] = 1;
	ks++;
      }
    }
    iter++;
  }
  return ks;
}

void VReleaseStorage(VAttrList list)
{
  VImage src;
  VAttrListPosn posn;

  for(VFirstAttr(list, & posn); VAttrExists(& posn); VNextAttr(& posn)) {
    if(VGetAttrRepn(& posn) != VImageRepn) continue;
    VGetAttrValue(& posn, NULL, VImageRepn, & src);
    VDestroyImage(src);
    src = NULL;
  }
}

/*  Check mask coverage */
VAttrList VCheckMaskCoverage(VStringConst *in_filenames,size_t n,VImage mask)
{
  size_t i;
  VAttrList list=NULL;
  VAttrList geolist = NULL;
  FILE *fp=NULL;

  for (i=0; i<n; i++) {
    fp = VOpenInputFile (in_filenames[i], TRUE);
    list = VReadFile (fp, NULL);
    if (! list)  VError("Error reading image");
    fclose(fp);

    /* read geometry info, if unknown */
    if (geolist == NULL) geolist = VGetGeoInfo(list);

    long count = VMaskCoverage(list,mask);
    if (count > 100) fprintf(stderr," incomplete mask_coverage  %s:  %ld\n",in_filenames[i],count);
    VReleaseStorage(list);
    list = NULL;
  }
  return geolist;
}

gsl_matrix_float **VReadImageData(VStringConst *in_filenames,VImage *src,VImage map,size_t n,size_t nvox,size_t first,size_t len)
{
  size_t i;
  VAttrList list=NULL;
  FILE *fp=NULL;

  /* allocate data struct */
  gsl_matrix_float **X = VCalloc(n,sizeof(gsl_matrix_float));
  for (i=0; i<n; i++) {
    X[i] = gsl_matrix_float_calloc(nvox,len);
    if (!X[i]) VError(" err allocating data matrix");
  }
  /* read data */
  for (i=0; i<n; i++) {
    fp = VOpenInputFile (in_filenames[i], TRUE);
    list = VReadFile (fp, NULL);
    if (! list)  VError("Error reading image");
    fclose(fp);
    VReadImagePointer(list,src);
    VDataMatrix(src,(int)first,(int)len,map,X[i]);
    VCheckMatrix(X[i]);
    VReleaseStorage(list);
    list = NULL;
  }
  return X;
}

VDictEntry ADJDict[] = {
  { "6", 0 },
  { "18", 1 },
  { "26", 2 },
  { NULL }
};

VDictEntry TypeDict[] = {
  { "SNR", 0 },
  { "median", 1 },
  { NULL }
};

VDictEntry MetricDict[] = {
  { "pearson", 0 },
  { "spearman", 1 },
  { NULL }
};

int main (int argc,char *argv[])
{
  static VArgVector in_files1;
  static VArgVector in_files2;
  static VString  out_filename = "";
  static VString  mask_filename = "";
  static VString  roi_filename = "";
  static VString  hist_filename= "";
  static VShort   first = 0;
  static VShort   len = 0;
  static VLong    seed = 99402622;
  static VFloat   qthreshold = 0.99;
  static VFloat   elength = 5;
  static VFloat   noise_cutoff = 0.01;
  static VShort   adjdef = 2;
  static VShort   numperm = 0;
  static VShort   type = 0;
  static VShort   metric = 0;
  static VShort   nproc = 10;
  static VShort   step = 2;
  static VOptionDescRec  options[] = {
    {"in1", VStringRepn, 0, & in_files1, VRequiredOpt, NULL,"Input files 1" },
    {"in2", VStringRepn, 0, & in_files2, VRequiredOpt, NULL,"Input files 2" },
    {"out", VStringRepn, 1, & out_filename, VOptionalOpt, NULL,"Output file" },
    {"mask", VStringRepn, 1, & mask_filename, VRequiredOpt, NULL,"Mask file" },
    {"roi", VStringRepn, 1, & roi_filename, VOptionalOpt, NULL,"ROI file" },
    {"histogram",VStringRepn,1,(VPointer) &hist_filename,VRequiredOpt,NULL,"Output histogram filename"},
    {"permutations",VShortRepn,1,(VPointer) &numperm,VOptionalOpt,NULL,"Number of permutations"},
    {"qthreshold",VFloatRepn,1,(VPointer) &qthreshold,VOptionalOpt,NULL,"Initial quantile threshold"},
    {"adj",VShortRepn,1,(VPointer) &adjdef,VOptionalOpt,ADJDict,"Definition of adjacency"},
    {"seed",VLongRepn,1,(VPointer) &seed,VOptionalOpt,NULL,"Seed for random number generator"},
    {"first",VShortRepn,1,(VPointer) &first,VOptionalOpt,NULL,"First timepoint to use within trial"},
    {"len",VShortRepn,1,(VPointer) &len,VOptionalOpt,NULL,"Number of timepoints to use, '0' to use all"},
    {"type",VShortRepn,1,(VPointer) &type,VOptionalOpt,TypeDict,"Type of trial average"},
    {"metric",VShortRepn,1,(VPointer) &metric,VOptionalOpt,MetricDict,"Correlation metric"},
    {"edgelength",VFloatRepn,1,(VPointer) &elength,VOptionalOpt,NULL,"Minimal edge length in voxels"},
    {"j",VShortRepn,1,(VPointer) &nproc,VOptionalOpt,NULL,"Number of processors to use, '0' to use all"},
    /*     {"noisecutoff",VFloatRepn,1,(VPointer) &noise_cutoff,VOptionalOpt,NULL,"estimate of noisy edges"}, */
    /*     {"step",VShortRepn,1,(VPointer) &step,VOptionalOpt,NULL,"first iteration ZMatrix step size"} */
  };
  FILE *fp=NULL;
  VAttrList list=NULL,list1=NULL,geolist=NULL;
  VAttrListPosn posn;
  VImage mask=NULL;
  size_t nvox=0,i=0;

  /* 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);
  }

  /* whether to produce output file, or only compute histogram */
  VBoolean histonly = FALSE;
  if (strlen(out_filename) < 2) histonly = TRUE;

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

  /* input filenames */
  size_t n1 = (size_t)in_files1.number;
  size_t n2 = (size_t)in_files2.number;
  if (n1 != n2) VError(" n1 != n2, %d %d",n1,n2);

  VStringConst *in_filenames1 = (VStringConst *) VCalloc(n1,sizeof(VStringConst));
  for (i=0; i<n1; i++) {
    in_filenames1[i] = ((VStringConst *) in_files1.vector)[i];
  }
  VStringConst *in_filenames2 = (VStringConst *) VCalloc(n2,sizeof(VStringConst));
  for (i=0; i<n2; i++) {
    in_filenames2[i] = ((VStringConst *) in_files2.vector)[i];
  }

  /* read mask */
  fp = VOpenInputFile (mask_filename, TRUE);
  list1 = VReadFile (fp, NULL);
  if (! list1) VError("Error reading mask file");
  fclose(fp);

  for (VFirstAttr (list1, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
    if (VGetAttrRepn (& posn) != VImageRepn) continue;
    VGetAttrValue (& posn, NULL,VImageRepn, & mask);
    if (VPixelRepn(mask) == VFloatRepn || VPixelRepn(mask) == VDoubleRepn) {
      mask = NULL;
      continue;
    }
  }
  if (mask == NULL) VError(" no mask found");
  int nslices = VImageNBands(mask);
  int nrows = VImageNRows(mask);
  int ncols = VImageNColumns(mask);

  /* read ROI */
  VAttrList list2 = NULL;
  VImage roi = NULL;
  if (strlen(roi_filename) > 1) {
    fp = VOpenInputFile (roi_filename, TRUE);
    list2 = VReadFile (fp, NULL);
    if (! list2) VError("Error reading roi file");
    fclose(fp);

    for (VFirstAttr (list2, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
      if (VGetAttrRepn (& posn) != VImageRepn) continue;
      VGetAttrValue (& posn, NULL,VImageRepn, & roi);
      if (VPixelRepn(roi) == VFloatRepn || VPixelRepn(roi) == VDoubleRepn) {
	roi = NULL;
	continue;
      }
    }
    if (roi == NULL) VError(" no roi found");
  }

  /* Check mask coverage */
  VAttrList geo=NULL;
  geo = VCheckMaskCoverage(in_filenames1,n1,mask);
  geo = VCheckMaskCoverage(in_filenames2,n2,mask);
  if (geolist == NULL) geolist = geo;

  /* voxel map */
  VImage map = VoxelMap(mask,&nvox);

  /* get image dimensions */
  fp = VOpenInputFile (in_filenames1[0], TRUE);
  list = VReadFile (fp, NULL);
  if (! list)  VError("Error reading image");
  fclose(fp);

  int nt=0;
  VImage *src = VImagePointer(list,&nt);

  if (first >= nt || first < 0) VError(" illegal value, first= %d, nt= %d",first,nt);
  if (len <= 0) len = nt-first;
  if (first + len >= nt) len = nt-first;
  if (len < 2) VError(" len= %d",len);
  fprintf(stderr," image: %d x %d x %d, nt: %d, nvox: %ld\n",(int)nslices,(int)nrows,(int)ncols,nt,nvox);

  /* read image data */
  gsl_matrix_float **X1 = VReadImageData(in_filenames1,src,map,n1,nvox,first,len);
  gsl_matrix_float **X2 = VReadImageData(in_filenames2,src,map,n2,nvox,first,len);

  /* voxel addresses */
  VImage mapimage = VCreateImage(nslices,nrows,ncols,VIntegerRepn);
  VFillImage(mapimage,VAllBands,(VInteger)(-1));
  for (i=0; i<nvox; i++) {
    int b = VPixel(map,0,0,i,VShort);
    int r = VPixel(map,0,1,i,VShort);
    int c = VPixel(map,0,2,i,VShort);
    VPixel(mapimage,b,r,c,VInteger) = (VInteger)i;
  }

  /* tmp storage for SNR time courses */
  int dim=0;
  if (type < 3) dim = len;
  else if (type == 3) dim = len*n1;
  else VError(" illegal type %d",type);
  gsl_matrix_float *SNR1 = gsl_matrix_float_calloc(nvox,dim);
  gsl_matrix_float *SNR2 = gsl_matrix_float_calloc(nvox,dim);

  /* ini zhist */
  size_t hbins = 10000;
  double hmin = 0.0,hmax = 1.001;
  gsl_histogram *TedHist = gsl_histogram_alloc (hbins);
  gsl_histogram_set_ranges_uniform (TedHist,hmin,hmax);
  gsl_histogram_reset(TedHist);

  /* fuer z-threshold */
  size_t nbins = 25000;
  hmin = -1.001,hmax = 1.001;
  gsl_histogram *ZvalHist = gsl_histogram_alloc (nbins);
  gsl_histogram_set_ranges_uniform (ZvalHist,hmin,hmax);
  gsl_histogram_reset(ZvalHist);

  /* ini random */
  size_t n = n1;
  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *rx = gsl_rng_alloc(T);
  gsl_rng_set(rx,(unsigned long int)seed);
  int *table = (int *) VCalloc(n,sizeof(int));
  for (i=0; i<n; i++) table[i] = 0;

  /* allocate output structure */
  size_t nedges=0;
  size_t fulldim = (nvox*(nvox+1)/2);
  double qx = 1.0-qthreshold;
  qx += 0.001;
  size_t nedges_estimated = (size_t) (qx*(double)fulldim);

  /* ini matrices */
  float *E=NULL;   /* matrix of edge densities */
  int *I=NULL;     /* row indices in matrix E */
  int *J=NULL;     /* col indices in matrix E */

  if (histonly == FALSE && noise_cutoff < 0.00001) {
    E = VCalloc(nedges_estimated,sizeof(float));  /* matrix of edge densities */
    I = VCalloc(nedges_estimated,sizeof(int));    /* row indices in matrix E */
    J = VCalloc(nedges_estimated,sizeof(int));    /* col indices in matrix E */
  }

  /* main loop */
  size_t s;
  int ks=0,nperm=0;
  int startperm=0;
  if (numperm > 0) startperm = 1;

  for (nperm = startperm; nperm <= numperm; nperm++) {

    fprintf(stderr,"\n perm %3d: ",nperm);
    if (nperm > 0) {
      ks = genperm(rx,table,n,(int)numperm);
    }
    for (s=0; s<n; s++) {
      fprintf(stderr,"%d",table[s]);
    }
    fprintf(stderr,"  %3d\n",ks);

    /* get SNR data */
    if (type == 0) {
      GetSNR(X1,X2,table,n,SNR1,metric);
      GetSNR(X2,X1,table,n,SNR2,metric);
    }
    if (type == 1) {
      GetMedian(X1,X2,table,n,SNR1,metric);
      GetMedian(X2,X1,table,n,SNR2,metric);
    }

    /* corr matrices */
    float zthr = ZMatrix(ZvalHist,SNR1,SNR2,n1,n2,roi,map,mapimage,(int)adjdef,
			 (float)elength,(float)qthreshold,(int)step,(int)metric);

    /* estimate number of truly needed edges, estimate noise */
    if (noise_cutoff > 0.0 && nperm == numperm && histonly == FALSE) {
      size_t old_estimate = nedges_estimated;
      nedges_estimated = EstimateEdges(E,I,J,old_estimate,TedHist,SNR1,SNR2,n1,n2,roi,map,mapimage,
				       (int)adjdef,elength,zthr,noise_cutoff,(int)metric);

      E = VCalloc(nedges_estimated,sizeof(float));  /* matrix of edge densities */
      I = VCalloc(nedges_estimated,sizeof(int));    /* row indices in matrix E */
      J = VCalloc(nedges_estimated,sizeof(int));    /* col indices in matrix E */
    }

    /* edge densities */
    fprintf(stderr," Computing edge densities...\n");
    nedges = EdgeDensity(E,I,J,nedges_estimated,TedHist,SNR1,SNR2,n1,n2,roi,map,mapimage,(int)adjdef,
			 elength,zthr,nperm,numperm,(int)1,noise_cutoff,(int)metric);
  }

  /* free storage */
  for (i=0; i<n; i++) {
    gsl_matrix_float_free(X1[i]);
    gsl_matrix_float_free(X2[i]);
  }
  gsl_matrix_float_free(SNR1);
  gsl_matrix_float_free(SNR2);
  VFree(table);

  /* print histogram */
  VPrintHistogram(TedHist,(int)numperm,hist_filename);

  if (histonly == TRUE) {
    fprintf (stderr," %s: done.\n", argv[0]);
    exit(0);
  }

  /* write edge density matrix in triplet sparse format */
  size_t edim = nedges*3;
  float filesize1 = (float)(edim*sizeof(float))/(1000.0*1000.0);
  float filesize2 = (float)(edim*sizeof(float))/(1000.0*1000.0*1000.0);
  fprintf(stderr," nedges= %ld,  filesize= %.3f MByte  (%.3f GByte)\n",nedges,filesize1,filesize2);

  /* output */
  VAttrList elist = VCreateAttrList();
  VAttrList ilist = VCreateAttrList();
  VAttrList jlist = VCreateAttrList();
  VBundle edgedens = VCreateBundle ("data",elist,nedges*sizeof(float),(VPointer)E);
  VBundle rowindex = VCreateBundle ("data",ilist,nedges*sizeof(int),(VPointer)I);
  VBundle colindex = VCreateBundle ("data",jlist,nedges*sizeof(int),(VPointer)J);

  VAttrList out_list = VCreateAttrList();
  VSetGeoInfo(geolist,out_list);
  VAppendAttr(out_list,"EdgeDensity",NULL,VBundleRepn,(VBundle)edgedens);
  VAppendAttr(out_list,"RowIndex",NULL,VBundleRepn,(VBundle)rowindex);
  VAppendAttr(out_list,"ColIndex",NULL,VBundleRepn,(VBundle)colindex);
  VAppendAttr(out_list,"nedges",NULL,VLongRepn,(VLong)nedges);
  VAppendAttr(out_list,"map",NULL,VImageRepn,map);
  VAppendAttr(out_list,"qthreshold",NULL,VFloatRepn,(VFloat)qthreshold);

  FILE *fp_out = VOpenOutputFile (out_filename, TRUE);
  if (! VWriteFile (fp_out, out_list)) exit (1);
  fclose(fp_out);
  fprintf (stderr," %s: done.\n\n", argv[0]);
  exit(0);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#pragma once

#include <string>
#include <gsl/span>

namespace crypto {

	/*
		Use the Windows CNG API to compute a MD5 hash of the
		provided data view and return it as a lowercase hex
		string.
	*/
	std::string MD5AsString(gsl::span<uint8_t> data);

}

ext: h

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
// MIT License
//
// Copyright (c) 2021 Daniel Robertson
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.

#ifndef HX711_DISCOVERY_H_82654AF1_3AA3_4885_9BE7_FA38117DE328
#define HX711_DISCOVERY_H_82654AF1_3AA3_4885_9BE7_FA38117DE328

#include "HX711.h"
#include <chrono>
#include <vector>
#include <pthread.h>
#include <sched.h>
#include <thread>
#include <algorithm>
#include <cmath>
#include <gsl/gsl_statistics.h>

namespace HX711 {

using namespace std::chrono;

struct TimingResult {
public:
    Value v;
    high_resolution_clock::time_point start;
    high_resolution_clock::time_point waitStart;
    high_resolution_clock::time_point waitEnd;
    high_resolution_clock::time_point convertStart;
    high_resolution_clock::time_point convertEnd;
    high_resolution_clock::time_point end;

    std::chrono::microseconds getWaitTime() const noexcept {
        return std::chrono::duration_cast<std::chrono::microseconds>(
            this->waitEnd - this->waitStart);
    }

    std::chrono::microseconds getConversionTime() const noexcept {
        return std::chrono::duration_cast<std::chrono::microseconds>(
            this->convertEnd - this->convertStart);
    }

    std::chrono::microseconds getTotalTime() const noexcept {
        return std::chrono::duration_cast<std::chrono::microseconds>(
            this->end - this->start);
    }

};

class TimingCollection : public std::vector<TimingResult> {
public:
    struct Stats {
        double min;
        double max;
        double med;
        double std;

        bool inRange(const double n) const noexcept {
            const double minDev = this->med - this->std;
            const double maxDev = this->med + this->std;
            return n >= minDev && n <= maxDev;
        }

    };

protected:

    Stats _statsFromVector(std::vector<double>& vec) const noexcept {

        Stats s;

        std::sort(vec.begin(), vec.end());

        s.min = gsl_stats_min(vec.data(), 1, vec.size());
        s.max = gsl_stats_max(vec.data(), 1, vec.size());
        s.med = gsl_stats_median_from_sorted_data(vec.data(), 1, vec.size());

        //double work[vec.size()];
        s.std = gsl_stats_sd(vec.data(), 1, vec.size());

        return s;

    }

public:

    Stats getWaitTimeStats() const noexcept {

        std::vector<double> vec;
        vec.reserve(this->size());

        for(auto it = this->cbegin(); it != this->cend(); ++it) {
            vec.push_back(static_cast<double>(it->getWaitTime().count()));
        }

        return _statsFromVector(vec);

    }

    Stats getConversionTimeStats() const noexcept {

        std::vector<double> vec;
        vec.reserve(this->size());

        for(auto it = this->cbegin(); it != this->cend(); ++it) {
            vec.push_back(static_cast<double>(it->getConversionTime().count()));
        }

        return _statsFromVector(vec);

    }

    Stats getTotalTimeStats() const noexcept {

        std::vector<double> vec;
        vec.reserve(this->size());

        for(auto it = this->cbegin(); it != this->cend(); ++it) {
            vec.push_back(static_cast<double>(it->getTotalTime().count()));
        }

        return _statsFromVector(vec);

    }

};

class Discovery : public HX711 {
public:

    Discovery(const int dataPin, const int clockPin, const Rate rate) :
        HX711(dataPin, clockPin, rate) {
            this->connect();
    }

    TimingCollection getTimings(const std::size_t samples) {

        using namespace std::chrono;

        TimingCollection vec;
        vec.reserve(samples);

        for(size_t i = 0; i < samples; ++i) {

            TimingResult tr;

            tr.start = high_resolution_clock::now();

            tr.waitStart = high_resolution_clock::now();
            while(!this->isReady()) ;
            tr.waitEnd = high_resolution_clock::now();

            tr.convertStart = high_resolution_clock::now();
            tr.v = this->readValue();
            tr.convertEnd = high_resolution_clock::now();

            tr.end = high_resolution_clock::now();

            vec.push_back(tr);

        }

        return vec;

    }

    std::vector<std::chrono::nanoseconds> getTimeToReady(const std::size_t samples) {

        using namespace std::chrono;

        std::vector<nanoseconds> timings;
        timings.reserve(samples);

        //do an initial read
        while(!this->isReady()) ;
        this->readValue();

        for(std::size_t i = 0; i < samples; ++i) {

            high_resolution_clock::time_point start = high_resolution_clock::now();
            while(!this->isReady()) ;
            high_resolution_clock::time_point end = high_resolution_clock::now();

            timings.push_back(end - start);

            this->readValue();

        }

        return timings;

    }

};
};
#endif


text:
//This does CELL (~soma) stage of GRU (gated recurrent unit) model.
//This requires each neuron to have 3 input time series, X, Xr, Xz,
//where X is the usual input and Xr, Xz the inputs for the reset and update gates.
//X, Xr, Xz are the output of a separate linear IN stage (weights and baises).

//For dim=0, R[:,t] = sig{Xr[:,t] + Ur*Y[:,t-1]} \n";
//           Z[:,t] = sig{Xz[:,t] + Uz*Y[:,t-1]} \n";
//           H[:,t] = R[:,t].*Y[:,t-1] \n";
//           Y[:,t] = Z[:,t].*Y[:,t-1] + (1-Z[:,t]).*tanh{X[:,t] + U*H[:,t]} \n";
//with sizes X, Xr, Xz: N x T \n";
//           U, Ur, Uz: N x N \n";
//           Y        : N x T \n";
//
//For dim=1, R[t,:] = sig{Xr[t,:] + Y[t-1,:]*Ur} \n";
//           Z[t,:] = sig{Xz[t,:] + Y[t-1,:]*Uz} \n";
//           H[t,:] = R[t,:].*Y[t-1,:] \n";
//           Y[t,:] = Z[t,:].*Y[t-1,:] + (1-Z[t,:]).*tanh{X[t,:] + H[t,:]*U} \n";
//with sizes X, Xr, Xz: T x N \n";
//           U, Ur, Uz: N x N \n";
//           Y        : T x N \n";
//
//where sig is the logistic (sigmoid) nonlinearity = 1/(1+exp(-x)),
//R is the reset gate, Z is the update gate, H is an intermediate (hidden) vector,
//U, Ur, Uz are NxN matrices, and Y is the output.

//Note that, the neurons of a layer are independent only if U, Ur, Uz are diagonal matrices.
//This is only really a CELL (~soma) stage in that case.

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

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

int gru3_s (float *Y, const float *X, const float *Xr, const float *Xz, const float *U, const float *Ur, const float *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int gru3_d (double *Y, const double *X, const double *Xr, const double *Xz, const double *U, const double *Ur, const double *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int gru3_inplace_s (float *X, const float *Xr, const float *Xz, const float *U, const float *Ur, const float *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int gru3_inplace_d (double *X, const double *Xr, const double *Xz, const double *U, const double *Ur, const double *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int gru3_s (float *Y, const float *X, const float *Xr, const float *Xz, const float *U, const float *Ur, const float *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;
    size_t nT, tN;

    float *R, *Z, *H;
    if (!(R=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru3_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(Z=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru3_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru3_s: problem with malloc. "); perror("malloc"); return 1; }

    if (N==1u)
    {
        //R[0] = 1.0f / (1.0f+expf(-Xr[0]));
        Z[0] = 1.0f / (1.0f+expf(-Xz[0]));
        Y[0] = (1.0f-Z[0]) * tanhf(X[0]);
        for (size_t t=1; t<T; ++t)
        {
            //R[0] = 1.0f / (1.0f+expf(-Xr[t]-Ur[0]*Y[t-1]));
            Z[0] = 1.0f / (1.0f+expf(-Xz[t]-Uz[0]*Y[t-1]));
            //H[0] = R[0] * Y[t-1];
            H[0] = Y[t-1] / (1.0f+expf(-Xr[t]-Ur[0]*Y[t-1]));
            Y[t] = Z[0]*Y[t-1] + (1.0f-Z[0])*tanhf(X[t]+U[0]*H[0]);
        }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                //R[n] = 1.0f / (1.0f+expf(-Xr[n]));
                Z[n] = 1.0f / (1.0f+expf(-Xz[n]));
                Y[n] = (1.0f-Z[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_scopy((int)N,&Xr[tN],1,R,1); cblas_scopy((int)N,&Xz[tN],1,Z,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&Y[tN-N],1,o,R,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&Y[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    //R[n] = 1.0f / (1.0f+expf(-R[n]));
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    //H[n] = R[n] * Y[tN-N+n];
                    H[n] = Y[tN-N+n] / (1.0f+expf(-R[n]));
                }
                cblas_scopy((int)N,&X[tN],1,&Y[tN],1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = Z[n]*Y[tN-N+n] + (1.0f-Z[n])*tanhf(Y[tN+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                //R[n] = 1.0f / (1.0f+expf(-Xr[nT]));
                Z[n] = 1.0f / (1.0f+expf(-Xz[nT]));
                Y[nT] = (1.0f-Z[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&Xr[t],(int)T,R,1); cblas_scopy((int)N,&Xz[t],(int)T,Z,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&Y[t-1],(int)T,o,R,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&Y[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    //R[n] = 1.0f / (1.0f+expf(-R[n]));
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    //H[n] = R[n] * Y[t-1+n*T];
                    H[n] = Y[t-1+n*T] / (1.0f+expf(-R[n]));
                }
                cblas_scopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = Z[n]*Y[t-1+nT] + (1.0f-Z[n])*tanhf(Y[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0f / (1.0f+expf(-Xz[nT]));
                Y[nT] = (1.0f-Z[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&Xr[t],(int)T,R,1); cblas_scopy((int)N,&Xz[t],(int)T,Z,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&Y[t-1],(int)T,o,R,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&Y[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = Y[t-1+n*T] / (1.0f+expf(-R[n]));
                }
                cblas_scopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = Z[n]*Y[t-1+nT] + (1.0f-Z[n])*tanhf(Y[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0f / (1.0f+expf(-Xz[n]));
                Y[n] = (1.0f-Z[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_scopy((int)N,&Xr[tN],1,R,1); cblas_scopy((int)N,&Xz[tN],1,Z,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&Y[tN-N],1,o,R,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&Y[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = Y[tN-N+n] / (1.0f+expf(-R[n]));
                }
                cblas_scopy((int)N,&X[tN],1,&Y[tN],1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = Z[n]*Y[tN-N+n] + (1.0f-Z[n])*tanhf(Y[tN+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru3_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int gru3_d (double *Y, const double *X, const double *Xr, const double *Xz, const double *U, const double *Ur, const double *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;
    size_t nT, tN;

    double *R, *Z, *H;
    if (!(R=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru3_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(Z=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru3_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru3_d: problem with malloc. "); perror("malloc"); return 1; }

    if (N==1u)
    {
        Z[0] = 1.0 / (1.0+exp(-Xz[0]));
        Y[0] = (1.0-Z[0]) * tanh(X[0]);
        for (size_t t=1; t<T; ++t)
        {
            Z[0] = 1.0 / (1.0+exp(-Xz[t]-Uz[0]*Y[t-1]));
            H[0] = Y[t-1] / (1.0+exp(-Xr[t]-Ur[0]*Y[t-1]));
            Y[t] = Z[0]*Y[t-1] + (1.0-Z[0])*tanh(X[t]+U[0]*H[0]);
        }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0 / (1.0+exp(-Xz[n]));
                Y[n] = (1.0-Z[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_dcopy((int)N,&Xr[tN],1,R,1); cblas_dcopy((int)N,&Xz[tN],1,Z,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&Y[tN-N],1,o,R,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&Y[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = Y[tN-N+n] / (1.0+exp(-R[n]));
                }
                cblas_dcopy((int)N,&X[tN],1,&Y[tN],1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = Z[n]*Y[tN-N+n] + (1.0-Z[n])*tanh(Y[tN+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0 / (1.0+exp(-Xz[nT]));
                Y[nT] = (1.0-Z[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&Xr[t],(int)T,R,1); cblas_dcopy((int)N,&Xz[t],(int)T,Z,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&Y[t-1],(int)T,o,R,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&Y[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = Y[t-1+n*T] / (1.0+exp(-R[n]));
                }
                cblas_dcopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = Z[n]*Y[t-1+nT] + (1.0-Z[n])*tanh(Y[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0 / (1.0+exp(-Xz[nT]));
                Y[nT] = (1.0-Z[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&Xr[t],(int)T,R,1); cblas_dcopy((int)N,&Xz[t],(int)T,Z,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&Y[t-1],(int)T,o,R,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&Y[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = Y[t-1+n*T] / (1.0+exp(-R[n]));
                }
                cblas_dcopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = Z[n]*Y[t-1+nT] + (1.0-Z[n])*tanh(Y[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0 / (1.0+exp(-Xz[n]));
                Y[n] = (1.0-Z[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_dcopy((int)N,&Xr[tN],1,R,1); cblas_dcopy((int)N,&Xz[tN],1,Z,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&Y[tN-N],1,o,R,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&Y[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = Y[tN-N+n] / (1.0+exp(-R[n]));
                }
                cblas_dcopy((int)N,&X[tN],1,&Y[tN],1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = Z[n]*Y[tN-N+n] + (1.0-Z[n])*tanh(Y[tN+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru3_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int gru3_inplace_s (float *X, const float *Xr, const float *Xz, const float *U, const float *Ur, const float *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;
    size_t nT, tN;

    float *R, *Z, *H;
    if (!(R=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru3_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(Z=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru3_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru3_inplace_s: problem with malloc. "); perror("malloc"); return 1; }

    if (N==1u)
    {
        //R[0] = 1.0f / (1.0f+expf(-Xr[0]));
        Z[0] = 1.0f / (1.0f+expf(-Xz[0]));
        X[0] = (1.0f-Z[0]) * tanhf(X[0]);
        for (size_t t=1; t<T; ++t)
        {
            //R[0] = 1.0f / (1.0f+expf(-Xr[t]-Ur[0]*X[t-1]));
            Z[0] = 1.0f / (1.0f+expf(-Xz[t]-Uz[0]*X[t-1]));
            //H[0] = R[0] * X[t-1];
            H[0] = X[t-1] / (1.0f+expf(-Xr[t]-Ur[0]*X[t-1]));
            X[t] = Z[0]*X[t-1] + (1.0f-Z[0])*tanhf(X[t]+U[0]*H[0]);
        }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                //R[n] = 1.0f / (1.0f+expf(-Xr[n]));
                Z[n] = 1.0f / (1.0f+expf(-Xz[n]));
                X[n] = (1.0f-Z[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_scopy((int)N,&Xr[tN],1,R,1); cblas_scopy((int)N,&Xz[tN],1,Z,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&X[tN-N],1,o,R,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&X[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    //R[n] = 1.0f / (1.0f+expf(-R[n]));
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    //H[n] = R[n] * X[tN-N+n];
                    H[n] = X[tN-N+n] / (1.0f+expf(-R[n]));
                }
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN+n] = Z[n]*X[tN-N+n] + (1.0f-Z[n])*tanhf(X[tN+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0f / (1.0f+expf(-Xz[nT]));
                X[nT] = (1.0f-Z[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&Xr[t],(int)T,R,1); cblas_scopy((int)N,&Xz[t],(int)T,Z,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,R,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = X[t-1+n*T] / (1.0f+expf(-R[n]));
                }
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = Z[n]*X[t-1+nT] + (1.0f-Z[n])*tanhf(X[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0f / (1.0f+expf(-Xz[nT]));
                X[nT] = (1.0f-Z[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&Xr[t],(int)T,R,1); cblas_scopy((int)N,&Xz[t],(int)T,Z,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,R,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = X[t-1+n*T] / (1.0f+expf(-R[n]));
                }
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = Z[n]*X[t-1+nT] + (1.0f-Z[n])*tanhf(X[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0f / (1.0f+expf(-Xz[n]));
                X[n] = (1.0f-Z[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_scopy((int)N,&Xr[tN],1,R,1); cblas_scopy((int)N,&Xz[tN],1,Z,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&X[tN-N],1,o,R,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&X[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = X[tN-N+n] / (1.0f+expf(-R[n]));
                }
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN+n] = Z[n]*X[tN-N+n] + (1.0f-Z[n])*tanhf(X[tN+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru3_inplace_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int gru3_inplace_d (double *X, const double *Xr, const double *Xz, const double *U, const double *Ur, const double *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;
    size_t nT, tN;

    double *R, *Z, *H;
    if (!(R=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru3_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(Z=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru3_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru3_inplace_d: problem with malloc. "); perror("malloc"); return 1; }

    if (N==1u)
    {
        Z[0] = 1.0 / (1.0+exp(-Xz[0]));
        X[0] = (1.0-Z[0]) * tanh(X[0]);
        for (size_t t=1; t<T; ++t)
        {
            Z[0] = 1.0 / (1.0+exp(-Xz[t]-Uz[0]*X[t-1]));
            H[0] = X[t-1] / (1.0+exp(-Xr[t]-Ur[0]*X[t-1]));
            X[t] = Z[0]*X[t-1] + (1.0-Z[0])*tanh(X[t]+U[0]*H[0]);
        }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0 / (1.0+exp(-Xz[n]));
                X[n] = (1.0-Z[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_dcopy((int)N,&Xr[tN],1,R,1); cblas_dcopy((int)N,&Xz[tN],1,Z,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&X[tN-N],1,o,R,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&X[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = X[tN-N+n] / (1.0+exp(-R[n]));
                }
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN+n] = Z[n]*X[tN-N+n] + (1.0-Z[n])*tanh(X[tN+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0 / (1.0+exp(-Xz[nT]));
                X[nT] = (1.0-Z[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&Xr[t],(int)T,R,1); cblas_dcopy((int)N,&Xz[t],(int)T,Z,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,R,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = X[t-1+n*T] / (1.0+exp(-R[n]));
                }
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = Z[n]*X[t-1+nT] + (1.0-Z[n])*tanh(X[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0 / (1.0+exp(-Xz[nT]));
                X[nT] = (1.0-Z[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&Xr[t],(int)T,R,1); cblas_dcopy((int)N,&Xz[t],(int)T,Z,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,R,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = X[t-1+n*T] / (1.0+exp(-R[n]));
                }
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = Z[n]*X[t-1+nT] + (1.0-Z[n])*tanh(X[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0 / (1.0+exp(-Xz[n]));
                X[n] = (1.0-Z[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_dcopy((int)N,&Xr[tN],1,R,1); cblas_dcopy((int)N,&Xz[tN],1,Z,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&X[tN-N],1,o,R,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&X[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = X[tN-N+n] / (1.0+exp(-R[n]));
                }
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN+n] = Z[n]*X[tN-N+n] + (1.0-Z[n])*tanh(X[tN+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru3_inplace_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif

ext: c


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <omp.h>
#include <float.h>
#include <math.h>
#include "sacio.h"
#include "parmt_utils.h"
#ifdef PARMT_USE_INTEL
#include <mkl.h>
#include <mkl_lapacke.h>
#include <mkl_cblas.h>
#else
#include <lapacke.h>
#include <cblas.h>
#endif
#include "parmt_mtsearch.h"
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"
#include "iscl/random/random.h"
#include "iscl/signal/convolve.h"
#include "iscl/time/time.h"

int testBlockSize(void);

int main()
{
    char *ffList[10] = {"rotation_data/B00103ZDS.sac\0",
                        "rotation_data/B00106ZSS.sac\0",
                        "rotation_data/B00101ZDD.sac\0",
                        "rotation_data/B00109ZEX.sac\0",
                        "rotation_data/B00104RDS.sac\0",
                        "rotation_data/B00107RSS.sac\0",
                        "rotation_data/B00102RDD.sac\0",
                        "rotation_data/B00110REX.sac\0",
                        "rotation_data/B00105TDS.sac\0",
                        "rotation_data/B00108TSS.sac\0"};
    struct parmtNoiseBasis_struct basis;
    struct sacData_struct sacZDS, sacZSS, sacZDD, sacZEX,
                          sacRDS, sacRSS, sacRDD, sacREX,
                          sacTDS, sacTSS, sac;
    double *CeInv, *CeInvDiag, *d, *dn, *G, *R, *X,
           *mts, *betas, *gammas, *kappas, *M0s,
           *phi, *sigmas, *thetas, *xrand, *var, *eye;
    double baz, betaMax, betaMin, gammaMax, gammaMin, kappaMin, kappaMax,
           m0Min, m0Max, sigmaMax, sigmaMin, thetaMax, thetaMin;
    double cmpaz, cmpinc, sigma, xnorm;
    int *lags, i, icomp, ierr, imt, imtopt, j, ldc, ldz, ng, nb, nk,
        nlags, ns, nt, nm, nmt, npgrns, npts, rank;
    bool lrescale = false;
    const double oneDeg = M_PI/180.0;  // one degree for padding
    const double az = 25.0;
    const int ldm = 8;
    const int blockSize = 32;
    iscl_init();
    // Discretize the moment tensor space
    ng = 11; // longitudes
    nb = 11; // colatitudes
    nk = 11;  // strike angles
    ns = 11;  // rake
    nt = 11;  // dip
    nm = 4;  // magnitudes
    betaMin = 0.0 + oneDeg;
    betaMax = M_PI - oneDeg;
    gammaMin =-M_PI/6.0 + oneDeg;
    gammaMax = M_PI/6.0 - oneDeg;
    kappaMin = 0.0 + oneDeg;
    kappaMax = 2.0*M_PI - oneDeg;
    thetaMin = 0.0 + oneDeg;
    thetaMax = 0.5*M_PI - oneDeg;
    sigmaMin =-0.5*M_PI + oneDeg;
    sigmaMax = 0.5*M_PI - oneDeg;
    m0Min = 1.0/sqrt(2.0);
    m0Max = 2.0/sqrt(2.0);
    //m0Max = 1.e8/sqrt(2.0) + (double) (nm - 1)*1.e8/sqrt(2.0);
    //omp_set_num_threads(nthreads);
    mkl_set_num_threads(1);
    // Discretize the MT space
    M0s    = array_linspace64f(m0Min,    m0Max,    nm, &ierr);
    betas  = array_linspace64f(betaMin,  betaMax,  nb, &ierr);
    gammas = array_linspace64f(gammaMin, gammaMax, ng, &ierr);
    kappas = array_linspace64f(kappaMin, kappaMax, nk, &ierr);
    thetas = array_linspace64f(thetaMin, thetaMax, nt, &ierr);
    sigmas = array_linspace64f(sigmaMin, sigmaMax, ns, &ierr);
    nmt = ng*nb*nk*ns*nt*nm;
    mts = memory_calloc64f(nmt*ldm);
    ierr = parmt_discretizeMT64f(ng, gammas,
                                 nb, betas,
                                 nm, M0s,
                                 nk, kappas,
                                 nt, thetas,
                                 ns, sigmas,
                                 ldm, nmt, mts);
    if (ierr != 0)
    {
        printf("Failed to discretize MTs\n");
        return EXIT_FAILURE;
    }
    memory_free64f(&M0s);
    memory_free64f(&betas);
    memory_free64f(&gammas);
    memory_free64f(&kappas);
    memory_free64f(&thetas);
    memory_free64f(&sigmas);
    // Read the green's functions generated by CPS
    memset(&sacZDS, 0, sizeof(struct sacData_struct));
    memset(&sacZSS, 0, sizeof(struct sacData_struct));
    memset(&sacZDD, 0, sizeof(struct sacData_struct));
    memset(&sacZEX, 0, sizeof(struct sacData_struct));
    memset(&sacRDS, 0, sizeof(struct sacData_struct));
    memset(&sacRSS, 0, sizeof(struct sacData_struct));
    memset(&sacRDD, 0, sizeof(struct sacData_struct));
    memset(&sacREX, 0, sizeof(struct sacData_struct));
    memset(&sacTDS, 0, sizeof(struct sacData_struct));
    memset(&sacTSS, 0, sizeof(struct sacData_struct));
    ierr  = sacio_readTimeSeriesFile(ffList[0], &sacZDS);
    ierr += sacio_readTimeSeriesFile(ffList[1], &sacZSS);
    ierr += sacio_readTimeSeriesFile(ffList[2], &sacZDD);
    ierr += sacio_readTimeSeriesFile(ffList[3], &sacZEX);
    ierr += sacio_readTimeSeriesFile(ffList[4], &sacRDS);
    ierr += sacio_readTimeSeriesFile(ffList[5], &sacRSS);
    ierr += sacio_readTimeSeriesFile(ffList[6], &sacRDD);
    ierr += sacio_readTimeSeriesFile(ffList[7], &sacREX);
    ierr += sacio_readTimeSeriesFile(ffList[8], &sacTDS);
    ierr += sacio_readTimeSeriesFile(ffList[9], &sacTSS);
    if (ierr != 0)
    {
        printf("Failed to load fundmaental faults\n");
        return EXIT_FAILURE;
    }
    // Set the green's functions
    npgrns = MIN(64, sacZDS.npts);
    npts = npgrns;
    icomp = 1;     // vertical
    cmpaz = 0.0;   // channel azimuth
    cmpinc =-90.0; // channel orientation (-90+ up 0 is horizontal)
    G = memory_calloc64f(6*npgrns);
    printf("Setting Green's functions...\n");
    baz = fmod(az + 180, 360.0);
    ierr = parmt_utils_ff2mtGreens64f(npgrns, icomp,
                                      az, baz,
                                      cmpaz, cmpinc,
                                      sacZDS.data, sacZSS.data,
                                      sacZDD.data, sacZEX.data,
                                      sacRDS.data, sacRSS.data,
                                      sacRDD.data, sacREX.data,
                                      sacTDS.data, sacTSS.data,
                                      &G[0*npgrns], &G[1*npgrns], &G[2*npgrns],
                                      &G[3*npgrns], &G[4*npgrns], &G[5*npgrns]);
    if (ierr != 0)
    {
        printf("Failed to set Green's functions\n");
        return EXIT_FAILURE;
    }
    // Free intermediate memory
    sacio_free(&sacZDS);
    sacio_free(&sacZSS);
    sacio_free(&sacZDD);
    sacio_free(&sacZEX);
    sacio_free(&sacRDS);
    sacio_free(&sacRSS);
    sacio_free(&sacRDD);
    sacio_free(&sacREX);
    sacio_free(&sacTDS);
    sacio_free(&sacTSS);
    // make a noise matrix
    ldc = npts;
    CeInv = memory_calloc64f(npts*ldc);
    R = memory_calloc64f(npts*ldc);
    X = memory_calloc64f((2*npts-1)*ldc);
    //xrand = random_randn64f(npts, 0.0, 0.01, &ierr);
    xrand = random_rand64f(npts, &ierr);
    xrand[0] = 0.0;
    xrand[npts-1] = 0.0;
    xnorm = cblas_dnrm2(npts, xrand, 1);
    cblas_dscal(npts, 1.0/xnorm, xrand, 1);
    ierr = convolve_corrmtx64f_work(npts, xrand,
                                    npts-1, ldc, X,
                                    CORRMTX_AUTOCORRELATION,
                                    true, ldc, R);
FILE *fmat = fopen("xcmat.txt", "w");
for (i=0; i<npts; i++)
{
  for (j=0; j<npts; j++)
  {
   fprintf(fmat, "%d %d %e\n", i+1, npts-j, R[i*ldc+j]);
  }
fprintf(fmat,"\n");
}
fclose(fmat);
    memory_free64f(&X);
    cblas_dcopy(npts*ldc, R, 1, CeInv, 1);
    // Make a noise basis
    ierr = parmt_utils_getNoiseBasisFromCtC64f(npts, ldc, R, &basis);
    if (ierr != 0)
    {
        printf("Failed to compute eigenbasis\n");
        return EXIT_FAILURE;
    }
    // Compute inverse of R
    ierr = LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'U', npts, CeInv, ldc);
    if (ierr != 0)
    {
        printf("Failed to compute chol(CeInv)\n");
        return EXIT_FAILURE;
    }
    ierr = LAPACKE_dpotri(LAPACK_ROW_MAJOR, 'U', npts, CeInv, ldc);
    if (ierr != 0)
    {
        printf("Failed to compute CeInv %d\n", ierr);
        return EXIT_FAILURE;
    }
    // The KL basis uses sqrt(eig) but for consistency with inv(Ce) need
    // just the regular old eigenvalues
    for (i=0; i<basis.nvals; i++)
    {
        basis.sqrtEvals[i] = pow(basis.sqrtEvals[i], 2);
    }
    // Copy upper to lower
    for (i=0; i<npts; i++)
    {
        for (j=0; j<i; j++)
        {
            CeInv[i*ldc+j] = CeInv[j*ldc+i];
        }
    }
    CeInvDiag = memory_calloc64f(npts);
    for (i=0; i<npts; i++){CeInvDiag[i] = CeInv[ldc*i+i];}
    // Make a synthetic
    imtopt = MIN(nmt-1, MAX(0, (3*nmt)/8));
    d = memory_calloc64f(npts);
    cblas_dgemv(CblasColMajor, CblasNoTrans,
                npgrns, 6, 1.0, G, npgrns,
                &mts[ldm*imtopt], 1, 0.0, d, 1);
    // make a noisy synthetic so Ei doesn't blow up
double relError = 0.2;
rank = parmt_utils_getEffectiveRank(relError, basis);
printf("%d\n", rank);
//    rank = basis.npts/4;
//printf("%d\n", rank);
//getchar();
printf("Effective rank used: %d\n", rank);
    sigma = 1.0;
    double snr = 15.0;
    dn = memory_calloc64f(npts);
    parmt_utils_makeKLNoise(rank, npts, basis, dn);
    double xsum = array_sum64f(npts, dn, &ierr);
    for (i=0; i<npts; i++){dn[i] = dn[i] - xsum/(double) npts;}
    double signalEnergy = cblas_dnrm2(npts, d, 1); // signal variance
    double noiseEnergy = cblas_dnrm2(npts, dn, 1); // noise variance
printf("%f %f\n", signalEnergy, noiseEnergy);
    double xvar = 0.0;
    double xfact = (1.0/noiseEnergy)*signalEnergy/pow(10.0, 0.05*snr);
    for (i=0; i<npts; i++)
    {
        dn[i] = xfact*dn[i] + d[i];
    }
    sigma = 0.0;
    for (i=0; i<npts; i++)
    {
        sigma = sigma + pow(dn[i] - d[i], 2); // noise variance
    }
    sigma = sqrt(sigma);
    //sigma = pow(cblas_dnrm2(npts, dn, 1), 2);
printf("variance: %f\n", sigma);
//sigma = 0.002;
printf("%e %e %e %e\n", xfact, sigma, d[14], dn[14]);
    //------------------------------------------------------------------------//
    //                        unlagged conventional                           //
    //------------------------------------------------------------------------//
    // run the inversion for a diagonal matrix
    phi = memory_calloc64f(nmt);
    time_tic();
    nlags = 0;
    ierr = parmt_mtSearch64f(ldm, nlags, nmt,
                             npts, ldc,
                             true,
                             &G[0*npgrns], &G[1*npgrns], &G[2*npgrns],
                             &G[3*npgrns], &G[4*npgrns], &G[5*npgrns],
                             CeInvDiag, mts, d, phi);
    if (ierr != 0)
    {
        printf("Failed calling mtSearch64f 1\n");
        return EXIT_FAILURE;
    }
    imt = array_argmin64f(nmt, phi, &ierr);
    if (imt != imtopt)
    {
        printf("Failed to recover optimal index in test 1\n");
        return EXIT_FAILURE;
    }
    printf("Unlagged diagonal matrix time: %f (s)\n", time_toc());
    time_tic();
 printf("est + true optimal indices: %d %d %e\n", imt, imtopt, phi[imt]);

    ierr = parmt_mtSearch64f(ldm, nlags, nmt,
                             npts, ldc,
                             false,
                             &G[0*npgrns], &G[1*npgrns], &G[2*npgrns],
                             &G[3*npgrns], &G[4*npgrns], &G[5*npgrns],
                             CeInv, mts, dn, phi);
    if (ierr != 0)
    {
        printf("Failed calling mtSearch64f 2\n");
        return EXIT_FAILURE;
    }
    imt = array_argmin64f(nmt, phi, &ierr);
    if (imt != imtopt)
    {
        printf("Failed to recover optimal index in test 2 %d %d\n", imt, imtopt);
        //return EXIT_FAILURE;
    }
imt = array_argmin64f(nmt, phi, &ierr);
    printf("est + true optimal indices: %d %d %e\n", imt, imtopt, phi[imt]);
    printf("Unlagged time: %f (s)\n", time_toc());

    lags = memory_calloc32i(nmt);
    var = memory_calloc64f(npts);
    nlags = 0;
    eye = array_set64f(npts, 1.0, &ierr);
    time_tic();
    ierr = parmt_mtSearchL164f(ldm, nmt, npts, blockSize,
                               nlags, false,
                               lrescale,
                               &G[0*npgrns], &G[1*npgrns], &G[2*npgrns],
                               &G[3*npgrns], &G[4*npgrns], &G[5*npgrns],
                               eye, mts, d, phi, var, lags);
    imt = array_argmax64f(nmt, phi, &ierr);
    printf("L1 est + true optimal indices: %d %d %e %e\n",
           imt, imtopt, phi[imt], phi[imtopt]);
    printf("Unlagged general matrix time: %f (s)\n", time_toc());
    nlags = 4;
omp_set_num_threads(1);
    ierr = parmt_mtSearchL164f(ldm, nmt, npts, blockSize,
                               nlags, true,
                               lrescale,
                               &G[0*npgrns], &G[1*npgrns], &G[2*npgrns],
                               &G[3*npgrns], &G[4*npgrns], &G[5*npgrns],
                               eye, mts, d, phi, var, lags);
    imt = array_argmax64f(nmt, phi, &ierr);
    printf("L1 lagged est + true optimal indices: %d %d %d %e %e\n",
           imt, imtopt, lags[imt], phi[imt], phi[imtopt]);
    printf("Lagged time: %f (s)\n", time_toc());

    time_tic();
    ierr = parmt_mtSearchL164f(ldm, nmt, npts, blockSize,
                               nlags, true,
                               true, //lrescale,
                               &G[0*npgrns], &G[1*npgrns], &G[2*npgrns],
                               &G[3*npgrns], &G[4*npgrns], &G[5*npgrns],
                               eye, mts, d, phi, var, lags);
    int imtNew = array_argmax64f(nmt, phi, &ierr);
    if (imtNew != imt)
    {
        if (fabs(phi[imtNew] - phi[imt]) > 1.e-14)
        {
            printf("Failed to locate lrescaled script %e\n",
                   fabs(phi[imtNew] - phi[imt]));
            return EXIT_FAILURE;
        }
    }
    printf("L1 rescaled lagged est + true optimal indices: %d %d %d %e %e\n",
           imt, imtopt, lags[imt], phi[imt], phi[imtopt]);
    printf("Rescaled lagged time: %f (s)\n", time_toc());

    // repeat with lags
    ierr = testBlockSize();
/*
    nlags = 10;
    eye = array_set64f(npts, 1.0, &ierr);
    ierr = parmt_mtSearchL164f(ldm, nmt, npts, blockSize,
                               nlags, true,
                               &G[0*npgrns], &G[1*npgrns], &G[2*npgrns],
                               &G[3*npgrns], &G[4*npgrns], &G[5*npgrns],
                               eye, mts, d, phi, var, lags);
    imt = array_argmax64f(nmt, phi, &ierr);
    printf("L1 est + true optimal indices: %d %d %e %d\n", imt, imtopt, phi[imt], lags[imt]);
    printf("Unlagged general matrix time: %f (s)\n", time_toc());
*/
/*
    //------------------------------------------------------------------------//
    //                             unlagged kl                                //
    //------------------------------------------------------------------------//
    time_tic();
    ldz = basis.lde;
    ierr = parmt_mtSearchKL64f(nlags, nmt,
                               npts, ldz,
                               rank,
                               sigma,
                               &G[0*npgrns], &G[1*npgrns], &G[2*npgrns],
                               &G[3*npgrns], &G[4*npgrns], &G[5*npgrns],
                               basis.sqrtEvals,
                               basis.evects,
                               mts, dn, phi);
    if (ierr != 0)
    {
        printf("Failed calling mtSearchKL64f unlagged\n");
        return EXIT_FAILURE;
    }
    imt = array_argmax64f(nmt, phi, &ierr);
*/
/*
    if (imt != imtopt)
    {
        printf("Failed to locate unlagged kl optimum\n");
        return EXIT_FAILURE;
    }
*/
/*
    printf("est + true optimal indices: %d %d %e %e %e\n", imt, imtopt,
           phi[imt], phi[imtopt], phi[array_argmin64f(nmt, phi)]);
    printf("%e\n", array_sum64f(nmt, phi, &ierr));
    printf("Unlagged general kl time %f (s)\n", time_toc());
*/
//TEST:;
/*
double *est = memory_calloc64f(npts);
    cblas_dgemv(CblasColMajor, CblasNoTrans,
                npgrns, 6, 1.0, G, npgrns,
                &mts[ldm*imt], 1, 0.0, est, 1);
FILE *fn = fopen("noise.txt","w");
for (i=0; i<npts; i++)
{
 fprintf(fn, "%d %e %e %e %e\n", i, d[i], dn[i], est[i], d[i] - dn[i]);
}
fclose(fn);
*/
    // Free memory
    parmt_utils_freeNoiseBasis(&basis);
    memory_free64f(&G);
    memory_free64f(&CeInv);
    memory_free64f(&CeInvDiag);
    memory_free64f(&mts);
    memory_free64f(&d);
    memory_free64f(&phi);
    memory_free64f(&var);
    memory_free64f(&eye);
    iscl_finalize();
    return EXIT_SUCCESS;
}
//============================================================================//
int testPolarityMTSearch(void)
{
    const char *fcnm = "testPolarityMTSearch\0";
    const int blockSizes[4] = {1, 2, 32, 64};
    int ierr, nb, ng, nk, nm, ns, nt;
    double *betas, *gammas, *kappas, *mts, *M0s, *thetas, *sigmas;
    double betaMin, betaMax, gammaMin, gammaMax, kappaMin, kappaMax,
           m0Min, m0Max, sigmaMin, sigmaMax, thetaMin, thetaMax;
    const double oneDeg = M_PI/180.0;  // one degree for padding
    // Discretize the moment tensor space
    ng = 5; // longitudes
    nb = 6; // colatitudes
    nk = 7; // strike angles
    ns = 8; // rake
    nt = 9; // dip
    nm = 1; // magnitudes
    betaMin = 0.0 + oneDeg;
    betaMax = M_PI - oneDeg;
    gammaMin =-M_PI/6.0 + oneDeg;
    gammaMax = M_PI/6.0 - oneDeg;
    kappaMin = 0.0 + oneDeg;
    kappaMax = 2.0*M_PI - oneDeg;
    thetaMin = 0.0 + oneDeg;
    thetaMax = 0.5*M_PI - oneDeg;
    sigmaMin =-0.5*M_PI + oneDeg;
    sigmaMax = 0.5*M_PI - oneDeg;
    m0Min = 1.0/sqrt(2.0);
    m0Max = 2.0/sqrt(2.0);
    ierr = parmt_discretizeCells64f(nb, betaMin, betaMax,
                                    ng, gammaMin, gammaMax,
                                    nk, kappaMin, kappaMax,
                                    ns, sigmaMin, sigmaMax,
                                    nt, thetaMin, thetaMax,
                                    nm, m0Min, m0Max, false,
                                    &betas, &gammas, &kappas,
                                    &sigmas, &thetas, &M0s);
    if (ierr != 0)
    {
        printf("%s: Error making cell discretization\n", fcnm);
        return EXIT_FAILURE;
    }
    memory_free64f(&betas);
    memory_free64f(&gammas);
    memory_free64f(&kappas);
    memory_free64f(&sigmas);
    memory_free64f(&thetas);
    memory_free64f(&M0s);
    // Investigate different block sizes
    return EXIT_SUCCESS;
}
//============================================================================//
/*!
 * @brief Tests the effect of blocksize in gemm on performance
 */
int testBlockSize(void)
{
    const char *fcnm = "testBlockSize\0";
    FILE *fout;
    char fname[PATH_MAX];
    double *G, *d, *eye, *phi, *phiOri, *var, time, dnorm, xmin;
    int *lags = NULL;
    int i, ierr, imt, nmt, nt0, npmax, npts, npgrns, pSize;
    const int nblocks = 10;
    const int nps = 6;
    int pSizes[6] = {32, 64, 128, 256, 512, 1024};
    int blockSizes[10] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512};
    const int nlags = 0;
    const int ldm = 6;
    const bool lwantLags = false;
    int imtopt, ng, nb, nk, ns, nt, nm;
    const bool lrescale = false;
    double *betas, *gammas, *kappas, *mts, *M0s, *thetas, *sigmas;
    double betaMin, betaMax, gammaMin, gammaMax, kappaMin, kappaMax,
           m0Min, m0Max, sigmaMin, sigmaMax, thetaMin, thetaMax;
    const double oneDeg = M_PI/180.0;  // one degree for padding
    // disretize the moment tensor
    // Discretize the moment tensor space
    ng = 5; // longitudes
    nb = 5; // colatitudes
    nk = 7; // strike angles
    ns = 7; // rake
    nt = 7; // dip
    nm = 1; // magnitudes
    betaMin = 0.0 + oneDeg;
    betaMax = M_PI - oneDeg;
    gammaMin =-M_PI/6.0 + oneDeg;
    gammaMax = M_PI/6.0 - oneDeg;
    kappaMin = 0.0 + oneDeg;
    kappaMax = 2.0*M_PI - oneDeg;
    thetaMin = 0.0 + oneDeg;
    thetaMax = 0.5*M_PI - oneDeg;
    sigmaMin =-0.5*M_PI + oneDeg;
    sigmaMax = 0.5*M_PI - oneDeg;
    m0Min = 1.0/sqrt(2.0);
    m0Max = 2.0/sqrt(2.0);
    // Discretize the MT space
    M0s    = array_linspace64f(m0Min,    m0Max,    nm, &ierr);
    betas  = array_linspace64f(betaMin,  betaMax,  nb, &ierr);
    gammas = array_linspace64f(gammaMin, gammaMax, ng, &ierr);
    kappas = array_linspace64f(kappaMin, kappaMax, nk, &ierr);
    thetas = array_linspace64f(thetaMin, thetaMax, nt, &ierr);
    sigmas = array_linspace64f(sigmaMin, sigmaMax, ns, &ierr);
    nmt = ng*nb*nk*ns*nt*nm;
    imtopt = nmt/2;
    mts = memory_calloc64f(nmt*ldm);
    ierr = parmt_discretizeMT64f(ng, gammas,
                                 nb, betas,
                                 nm, M0s,
                                 nk, kappas,
                                 nt, thetas,
                                 ns, sigmas,
                                 ldm, nmt, mts);
    if (ierr != 0)
    {
        printf("Failed to discretize MTs\n");
        return EXIT_FAILURE;
    }
    memory_free64f(&M0s);
    memory_free64f(&betas);
    memory_free64f(&gammas);
    memory_free64f(&kappas);
    memory_free64f(&thetas);
    memory_free64f(&sigmas);
    // make some random data for testing
    npmax = array_max32i(nps, pSizes, &ierr);
    G = random_rand64f(6*npmax, &ierr);
    // set space
    phiOri = memory_calloc64f(nmt);
    phi = memory_calloc64f(nmt);
    var = memory_calloc64f(npmax);
    d = memory_calloc64f(npmax);
    eye = array_set64f(npmax, 1.0, &ierr);
    // fix the number of threads
    nt0 = omp_get_num_threads();
    omp_set_num_threads(1);
    // loop on point sizes
    for (pSize=0; pSize<nps; pSize++)
    {
        // make synthetic data
        memset(fname, 0, PATH_MAX*sizeof(char));
        sprintf(fname, "blockSize_%d.txt", pSizes[pSize]);
        fout = fopen(fname, "w");
        npts = pSizes[pSize];
        npgrns = npts;
        cblas_dgemv(CblasColMajor, CblasNoTrans,
                    npgrns, 6, 1.0, G, npgrns, &mts[ldm*imtopt], 1, 0.0, d, 1);
        for (i=0; i<nblocks; i++)
        {
            time_tic();
            ierr = parmt_mtSearchL164f(ldm, nmt, npts, blockSizes[i],
                                       nlags, lwantLags,
                                       lrescale,
                                       &G[0*npgrns], &G[1*npgrns], &G[2*npgrns],
                                       &G[3*npgrns], &G[4*npgrns], &G[5*npgrns],
                                       eye, mts, d, phi, var, lags);
            time = time_toc();
            fprintf(fout, "%d %d %e %e %d\n",
                    npts, blockSizes[i], 1.e6*time/(double) nmt, time/(double) nmt, nmt);
            if (i == 0)
            {
                ierr = array_copy64f_work(nmt, phi, phiOri);
            }
            dnorm = array_normDiff64f(nmt, phi, phiOri, ONE_NORM, 1.0, &ierr);
            xmin = array_min64f(nmt, phi, &ierr);
            if (dnorm > DBL_EPSILON*(double) nmt)
            {
                printf("Error computing blocksize\n");
                ierr = 1;
                goto ERROR;
            }
            imt = array_argmax64f(nmt, phi, &ierr);
            if (imt != imtopt)
            {
                printf("Failed to locate optimum\n");
                ierr = 1;
                goto ERROR;
            }
        }
        fclose(fout);
    }
ERROR:;
    // restore number of threads
    omp_set_num_threads(nt0);
    // clear memory
    memory_free64f(&d);
    memory_free64f(&phiOri);
    memory_free64f(&phi);
    memory_free64f(&eye);
    memory_free64f(&var);
    memory_free64f(&mts);
    memory_free64f(&G);
    return 0;
}


text:
#pragma once

#include <gsl/span>
#include "simple_math.h"

constexpr auto VECTOR_SIZE = sizeof(float) * 4;

class CBufferBase;

struct CBufferAlign
{
	size_t size;
	explicit CBufferAlign(size_t size_ = VECTOR_SIZE) : size(size_) {}
};

__forceinline CBufferAlign cbuff_align(size_t size = VECTOR_SIZE)
{
	return CBufferAlign(size);
}

class ICBuffer
{
public:
	virtual ~ICBuffer() = default;
	virtual void write(CBufferBase& cbuf) const = 0;

	size_t cbuffer_size() const;

	template <typename T>
	static size_t cbuffer_size()
	{
		T t;
		return t.cbuffer_size();
	}
};

class CBufferBase
{
protected:
	size_t offset_ = 0;
	size_t alignment_ = 0;

public:
	virtual ~CBufferBase() = default;
	bool align(size_t size = VECTOR_SIZE);
	void add(size_t size);
	void reset();

	size_t offset() const { return offset_; }
	size_t alignment() const { return alignment_; }

	template <typename T>
	CBufferBase& operator<<(const T& data) = delete;

	template <typename T>
	CBufferBase& operator<<(const gsl::span<T>& data) = delete;

	template <typename T>
	CBufferBase& operator<<(const gsl::span<const T>& data) = delete;

	CBufferBase& operator<<(const CBufferAlign& align_of);

	template <typename T, size_t size>
	__forceinline CBufferBase& operator<<(const std::array<T, size>& array)
	{
		return *this << gsl::span<const T>(array);
	}

	template <typename T, size_t size>
	__forceinline CBufferBase& operator<<(const T(&array)[size])
	{
		return *this << gsl::span<const T>(array);
	}

	template <typename T>
	__forceinline CBufferBase& operator<<(const dirty_t<T>& value)
	{
		return *this << value.data();
	}

	virtual void write(const void* data, size_t size)
	{
		align(size);
		add(size);
	}
};

template <>
CBufferBase& CBufferBase::operator<<(const int32_t& data);
template <>
CBufferBase& CBufferBase::operator<<(const uint32_t& data);
template <>
CBufferBase& CBufferBase::operator<<(const float& data);
template <>
CBufferBase& CBufferBase::operator<<(const DirectX::SimpleMath::Matrix& data);
template <>
CBufferBase& CBufferBase::operator<<(const DirectX::SimpleMath::Vector2& data);
template <>
CBufferBase& CBufferBase::operator<<(const DirectX::SimpleMath::Vector3& data);
template <>
CBufferBase& CBufferBase::operator<<(const DirectX::SimpleMath::Vector4& data);
template <>
CBufferBase& CBufferBase::operator<<(const gsl::span<float>& data);
template <>
CBufferBase& CBufferBase::operator<<(const gsl::span<const float>& data);

template <>
__forceinline CBufferBase& CBufferBase::operator<<(const DWORD& data)
{
	return *this << static_cast<uint32_t>(data);
}

template <>
__forceinline CBufferBase& CBufferBase::operator<<(const bool& data)
{
	return *this << (data ? 1 : 0);
}

class CBufferWriter : public CBufferBase
{
	uint8_t* ptr = nullptr;

public:
	CBufferWriter(uint8_t* ptr_);

private:
	void write(const void* data, size_t size) override;
};


text:
#include <stdio.h>
#include <cblas.h>

double A[] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 };

double x[] = { 1.0, -1.0, 2.0 };
double y[] = { -1.0, 0.0, 2.0 };

int main() {
  int i;
  cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 2.5, A, 3, x, 1, 1.2, y, 1);
  for (i = 0; i < 3; ++i) { printf("%6.1f\n", y[i]); }

  return 0;
}

max_stars_repo_name: julnamoo/practice-linux

max_issues_repo_name: julnamoo/practice-linux

max_forks_repo_name: julnamoo/practice-linux


text:
/* Copyright (C) 2010-2019, The Regents of The University of Michigan.
 All rights reserved.

 This software was developed as part of the The Vulcan project in the Intelligent Robotics Lab
 under the direction of Benjamin Kuipers, kuipers@umich.edu. Use of this code is governed by an
 MIT-style License that can be found at "https://github.com/h2ssh/Vulcan".
*/

/**
 * \file     mpepc_optimizer.h
 * \author   Collin Johnson and Jong Jin Park
 *
 * Declaration of MPEPCOptimizer.
 */

#ifndef MPEPC_OPTIMIZER_H
#define MPEPC_OPTIMIZER_H

#include "core/pose.h"
#include "mpepc/trajectory/params.h"
#include "mpepc/trajectory/robot_trajectory_info.h"
#include <memory>
#include <nlopt.h>
#include <vector>

namespace vulcan
{
namespace mpepc
{

const int kOptimizationDimension = 5;

const int kIdxR = 0;
const int kIdxTheta = 1;
const int kIdxDelta = 2;
const int kIdxVGain = 3;
const int kIdxK1 = 4;

/**
 * optimizer_coords_t is a coordinate of optimization space for the MPEPCOptimizer,
 * i.e. the argument to the objective function to be minimized.
 */
struct optimizer_coords_t
{
    double r;
    double theta;
    double delta;
    double velocityGain;
    double k1;
    double k2;

    explicit optimizer_coords_t(double r = 0,
                                double theta = 0,
                                double delta = 0,
                                double velocityGain = 0,
                                double k1 = 1.0,
                                double k2 = 1.0)
    : r(r)
    , theta(theta)
    , delta(delta)
    , velocityGain(velocityGain)
    , k1(k1)
    , k2(k2)
    {
    }

    explicit optimizer_coords_t(const double x[kOptimizationDimension])
    : r(x[kIdxR])
    , theta(x[kIdxTheta])
    , delta(x[kIdxDelta])
    , velocityGain(x[kIdxVGain])
    , k1(x[kIdxK1])
    {
    }

    void toArray(double x[kOptimizationDimension]) const
    {
        x[kIdxR] = r;
        x[kIdxTheta] = theta;
        x[kIdxDelta] = delta;
        x[kIdxVGain] = velocityGain;
        x[kIdxK1] = k1;
    }

    bool operator==(const optimizer_coords_t& rhs) const
    {
        return (r == rhs.r) && (theta == rhs.theta) && (delta == rhs.delta) && (velocityGain == rhs.velocityGain)
          && (k1 == rhs.k1);
    }

    bool operator!=(const optimizer_coords_t& rhs) const { return !(*this == rhs); }
};

class RobotSimulator;
class TrajectoryEvaluator;
class TaskManifold;
struct motion_target_t;
struct trajectory_planner_debug_info_t;

/**
 * MPEPCTrajectoryOptimizer handles the optimization task within the MPEPC planner.
 * The optimizer recieves the simulator and the cost evaluator to use, along with
 * the state for the current planning problem, i.e. planner_data_t, and finds the
 * time-optimal trajectory in the form of a motion target.
 *
 * The optimizer handles the simulate/evaluate loop internally and finds the optimal
 * target along with the intermediate trajectories generated by the optimization
 * algorithm. NLopt package is used for the numerical optimzation, and a number of
 * different algorithms can be implemented.
 */
class MPEPCOptimizer
{

public:
    /**
     * Constructor for MPEPCOptimizer.
     *
     * \param    params          Parameters for the optimization
     */
    MPEPCOptimizer(const mpepc_optimizer_params_t& params);

    /**
     * Destructor for MPEPCOptimizer.
     */
    ~MPEPCOptimizer(void);

    /**
     * Setup sets up the optmizer with the robot simulator and trajectory evaluator to be used.
     *
     * \param    simulator       RobotSimulator to use for the trajectory generation.
     * \param    evaluator       TrajectoryEvaluator for calculating the cost of a trajectory.
     */
    void setup(RobotSimulator& simulator,
               TrajectoryEvaluator& evaluator);   // NOTE: This feels very unsafe. A better way?

    /**
     * setInitialGuesses adds a goal pose (if navigation task) and the previous optimum
     * to the initial guesses for local optimizer.
     *
     * \param    task                 TaskManifold that contains the goal.
     * \param    previousMotionTarget Previous output of the metric planner, i.e. the previous optimum.
     */
    void setInitialGuesses(const TaskManifold& task, const motion_target_t& previousMotionTarget);

    /**
     * run runs the optimizer and finds the best motion target by evaluating the expected
     * costs of an estimated trajectory generated from a candidate motion target.
     *
     * It returns the optimal target found, while the intermediate information including
     * all trajectories considered are stored in the provided debug info instance.
     *
     * \param    debugInfo        Place to store relevant information
     * \return   Optimal target found
     */
    motion_target_t runOptimizer(trajectory_planner_debug_info_t& debugInfo);

    /**
     * evaluateCost calculates the cost of a fixed-length trajectory represented by the given optimizer_target_t.
     * The length of the trajectory is specified in the given instance of the planner_data_t.
     *
     * \param    coords          Optimization variable, currently encoding a target pose in the robot frame and a
     * velocity gain for the kinematic control law. \return   Cost of the trajectory.
     */
    double evaluateCost(const optimizer_coords_t& coords);

    // get getSurvivability of the trajectroy associated with the optimizer coords. Will assert failure if the given
    // coordinate is not the same as the one used for cost evaluation
    double getChanceConstraint(const optimizer_coords_t& coords);

    /**
     * evaluateGradient calculates the gradient around a particular point, by a simple forward differencing.
     *
     * \param    coords          Optimization variable, currently encoding a target pose in the robot frame and a
     * velocity gain for the kinematic control law. \param    cost            Evaluated cost at the given coordinate.
     * \param    gradient        Calculated gradient for each dimension (r, theta, delta, velocity) [output]
     */
    void evaluateGradient(const optimizer_coords_t& coords, double cost, double gradient[kOptimizationDimension]);

    bool haveFoundSolution(void) const { return haveFoundSolution_; };

    //     /**
    //     * evaluateSingleMotionTarget is a method which resets the data and info struct and evaluate a single
    //     optimizer target.
    //     * The optimization_info_t will hold a single trajectory considered. TODO: Do i need this?
    //     *
    //     * \param candidateMotionTarget
    //     * \param data
    //     * \param info
    //     */
    //     double evaluateSingleMotionTarget(const motion_target_t& candidateMotionTarget, const planner_data_t& data,
    //     optimization_info_t& info);

private:
    // setting up the optimizer
    void setupGlobalOptimizer(void);
    void setupLocalOptimizer(void);

    // loading a predetermined set of initial points for local optimizer
    void loadGuessesFromFile(const std::string& filename);

    // checking if the target is within a specified bounds
    bool isTargetValid(const optimizer_coords_t& target) const;

    // selecting the best inital point for the local optimizer
    // assumes there are valid trajectories stored in optimizerInfo_
    void selectBestInitialPoint(const std::vector<optimizer_coords_t>& initialGuesses,
                                double x[kOptimizationDimension]);
    std::size_t bestTrajectoryIdx(void);

    // optimizers
    nlopt_opt globalOptimizer;
    nlopt_opt localOptimizer;

    optimizer_coords_t coordEvaluated_;

    double globalOptimizerSeed_[kOptimizationDimension];

    // simulator and evaluator
    RobotSimulator* simulator_;
    TrajectoryEvaluator* evaluator_;

    robot_trajectory_info_t trajectory;   // Trajectory to use for everything -- just reset it before each use

    // debug info
    trajectory_planner_debug_info_t* optimizerInfo_;

    // a set of initial guesses for optmization
    std::vector<optimizer_coords_t> defaultGuesses_;
    std::vector<optimizer_coords_t> newGuesses_;

    // inital robot state
    pose_t robotPose_;

    // indicator for finding a right solution
    bool haveFoundSolution_;

    // parameters
    mpepc_optimizer_params_t params_;
};

}   // namespace mpepc
}   // namespace vulcan

#endif   // MPEPC_OPTIMIZER_H

max_stars_repo_name: anuranbaka/Vulcan

max_issues_repo_name: anuranbaka/Vulcan

max_forks_repo_name: anuranbaka/Vulcan


text:
/* specfunc/gsl_sf_erf.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Author:  G. Jungman */

#ifndef __GSL_SF_ERF_H__
#define __GSL_SF_ERF_H__

#include <gsl/gsl_sf_result.h>

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

__BEGIN_DECLS

/* Complementary Error Function
 * erfc(x) := 2/Sqrt[Pi] Integrate[Exp[-t^2], {t,x,Infinity}]
 *
 * exceptions: none
 */
int gsl_sf_erfc_e(double x, gsl_sf_result * result);
double gsl_sf_erfc(double x);

/* Log Complementary Error Function
 *
 * exceptions: none
 */
int gsl_sf_log_erfc_e(double x, gsl_sf_result * result);
double gsl_sf_log_erfc(double x);

/* Error Function
 * erf(x) := 2/Sqrt[Pi] Integrate[Exp[-t^2], {t,0,x}]
 *
 * exceptions: none
 */
int gsl_sf_erf_e(double x, gsl_sf_result * result);
double gsl_sf_erf(double x);

/* Probability functions:
 * Z(x) :  Abramowitz+Stegun 26.2.1
 * Q(x) :  Abramowitz+Stegun 26.2.3
 *
 * exceptions: none
 */
int gsl_sf_erf_Z_e(double x, gsl_sf_result * result);
int gsl_sf_erf_Q_e(double x, gsl_sf_result * result);
double gsl_sf_erf_Z(double x);
double gsl_sf_erf_Q(double x);

/* Hazard function, also known as the inverse Mill's ratio.
 *
 *   H(x) := Z(x)/Q(x)
 *         = Sqrt[2/Pi] Exp[-x^2 / 2] / Erfc[x/Sqrt[2]]
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_hazard_e(double x, gsl_sf_result * result);
double gsl_sf_hazard(double x);

__END_DECLS

#endif /* __GSL_SF_ERF_H__ */

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


text:
#ifndef CWANNIER_HTIGHTBINDING_H
#define CWANNIER_HTIGHTBINDING_H

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>
#include "bstrlib/bstrlib.h"

typedef struct {
    int num_rs, num_bands;
    double *ras;
    double *rbs;
    double *rcs;
    double *degens;
    gsl_matrix_complex **Hrs;
} HTightBinding;

HTightBinding* ExtractHTightBinding(char *filePath);

void FreeHTightBinding(HTightBinding *Hrs);

void HkRecip(HTightBinding *Hrs, double k[3], gsl_matrix_complex *Hk);

gsl_matrix_complex* HrAtR(HTightBinding *Hrs, double R[3], double *degen);

#endif // CWANNIER_HTIGHTBINDING_H

max_stars_repo_name: tflovorn/cwannier

max_issues_repo_name: tflovorn/cwannier

max_forks_repo_name: tflovorn/cwannier