text:
//**************************************************************************\
//* This file is property of and copyright by the ALICE Project            *\
//* ALICE Experiment at CERN, All rights reserved.                         *\
//*                                                                        *\
//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *\
//*                  for The ALICE HLT Project.                            *\
//*                                                                        *\
//* Permission to use, copy, modify and distribute this software and its   *\
//* documentation strictly for non-commercial purposes is hereby granted   *\
//* without fee, provided that the above copyright notice appears in all   *\
//* copies and that both the copyright notice and this permission notice   *\
//* appear in the supporting documentation. The authors make no claims     *\
//* about the suitability of this software for any purpose. It is          *\
//* provided "as is" without express or implied warranty.                  *\
//**************************************************************************

/// \file GPUO2InterfaceConfiguration.h
/// \author David Rohr

#ifndef GPUO2INTERFACECONFIGURATION_H
#define GPUO2INTERFACECONFIGURATION_H

#ifndef HAVE_O2HEADERS
#define HAVE_O2HEADERS
#endif
#ifndef GPUCA_TPC_GEOMETRY_O2
#define GPUCA_TPC_GEOMETRY_O2
#endif
#ifndef GPUCA_O2_INTERFACE
#define GPUCA_O2_INTERFACE
#endif

#include <memory>
#include <array>
#include <vector>
#include <functional>
#include <gsl/gsl>
#include "GPUSettings.h"
#include "GPUDataTypes.h"
#include "GPUHostDataTypes.h"
#include "GPUOutputControl.h"
#include "DataFormatsTPC/Constants.h"

class TH1F;
class TH1D;
class TH2F;

namespace o2
{
namespace tpc
{
class TrackTPC;
class Digit;
} // namespace tpc
namespace gpu
{
class TPCFastTransform;
class GPUReconstruction;
struct GPUSettingsO2;

struct GPUInterfaceQAOutputs {
  const std::vector<TH1F>* hist1;
  const std::vector<TH2F>* hist2;
  const std::vector<TH1D>* hist3;
};

struct GPUInterfaceOutputs : public GPUTrackingOutputs {
  GPUInterfaceQAOutputs qa;
};

// Full configuration structure with all available settings of GPU...
struct GPUO2InterfaceConfiguration {
  GPUO2InterfaceConfiguration() = default;
  ~GPUO2InterfaceConfiguration() = default;
  GPUO2InterfaceConfiguration(const GPUO2InterfaceConfiguration&) = default;

  // Settings for the Interface class
  struct GPUInterfaceSettings {
    int dumpEvents = 0;
    bool outputToExternalBuffers = false;
    float memoryBufferScaleFactor = 1.f;
    // These constants affect GPU memory allocation only and do not limit the CPU processing
    unsigned long maxTPCZS = 8192ul * 1024 * 1024;
    unsigned int maxTPCHits = 1024 * 1024 * 1024;
    unsigned int maxTRDTracklets = 128 * 1024;
    unsigned int maxITSTracks = 96 * 1024;
  };

  GPUSettingsDeviceBackend configDeviceBackend;
  GPUSettingsProcessing configProcessing;
  GPUSettingsGRP configGRP;
  GPUSettingsRec configReconstruction;
  GPUSettingsDisplay configDisplay;
  GPUSettingsQA configQA;
  GPUInterfaceSettings configInterface;
  GPURecoStepConfiguration configWorkflow;
  GPUCalibObjects configCalib;

  GPUSettingsO2 ReadConfigurableParam();

 private:
  friend class GPUReconstruction;
  GPUSettingsO2 ReadConfigurableParam_internal();
};

} // namespace gpu
} // namespace o2

#endif


text:
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "chaincolln.h"
#include "cokus.h"
#include "chain.h"
#include "config.h"
#include "parameters.h"
#include "irmutils.h"
#include "relation.h"
#include "domain.h"

#ifdef GSL
#include <gsl/gsl_randist.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_rng.h>
#endif

/*****************************************************************************
A collection of chains. We maintain several chains in case we want to run
Metropolis-coupled MCMC (MC^3). We also keep one spare chain so that
split-merge proposals can be tried without losing the current state.
*****************************************************************************/

struct chaincolln_str
{ int nchains;
  chain *chains;
  chain sparechain;
  char prefix[MAXSTRING];	/* unused right now */
  int itercount;		/* iteration count */
} chaincolln_str;

/* nchains already includes one extra for sparechain */
chaincolln chaincolln_create(int nchains, int ndomains, int nrelns, char *pref)
{
  int i;
  chaincolln cc;
  double temp;
  cc  = (chaincolln) my_malloc(sizeof(struct chaincolln_str));
  cc->chains = (chain *) my_malloc(nchains*sizeof(chain ));
  /* include spare chain */
  for (i = 0; i < nchains; i++) {
    temp = 1.0/(1+ps.temp*i);
    cc->chains[i] = chain_create(ndomains, nrelns, temp);
  }
  cc->nchains = nchains-1;
  cc->sparechain = cc->chains[nchains-1];
  strcpy(cc->prefix, pref);
  cc->itercount = 0;
  return cc;
}

/* read the configuration file and the graph */
chaincolln chaincolln_readdata(void) {
  FILE *fileptr, *initzsfile;
  int i, j, k, ndom, nreln, d, r, nitem, dim, maxclass, initclass, relcl, ndim,
	domlabel, clusterflag, itemind, nchains, cind, zind;
  int *domlabels, *participants, participant;
  double val;
  double nig[DISTSIZE];
  domain *doms;
  relation rn;
  int *initclasses, ***edgecounts, *relsizes;
  char prefix[MAXSTRING];

  chaincolln cc;
  chain c, c0;
#ifdef GSL
  gsl_rng *rng;
  const gsl_rng_type *T;
  gsl_permutation *perm ;
  size_t N;

  gsl_rng_env_setup();
  T = gsl_rng_default;
  rng = gsl_rng_alloc(T);
#endif

  fprintf(stdout,"A\n");
  nchains = ps.nchains+1;
  nig[0] = ps.m; nig[1] = ps.v; nig[2] = ps.a; nig[3] = ps.b;

  fileptr = fopen(ps.configfile,"r");
  if (fileptr == NULL) {
    fprintf(stderr, "couldn't read config file\n"); exit(1);
  }

  /* initial read of ps.configfile to get ps.maxdim, ps.maxrel, ps.maxitem,
     ps.maxclass */
  fscanf(fileptr, "%s", prefix);
  fscanf(fileptr, "%d %d", &ndom, &nreln);
  relsizes=  (int *) my_malloc(nreln*sizeof(int));
  ps.maxrel = nreln;
  ps.maxitem = 0; ps.maxclass = 0;
  for (d = 0; d < ndom; d++) {
    fscanf(fileptr, "%d %d %d %d", &nitem, &maxclass, &initclass, &clusterflag);
    if (nitem > ps.maxitem) {
      ps.maxitem = nitem;
    }
    if (maxclass > ps.maxclass) {
      ps.maxclass= maxclass;
    }
  }
  fprintf(stdout,"B\n");
  ps.maxdim = 0;
  for (r = 0; r < nreln; r++) {
    fscanf(fileptr, "%d", &ndim);
    relsizes[r] = ndim;
    if (ndim > ps.maxdim) {
      ps.maxdim = ndim;
    }
    for (dim=0; dim < ndim; dim++) {
      fscanf(fileptr, "%d", &domlabel);
    }
  }
  fclose(fileptr);

  fprintf(stdout,"C\n");
  domlabels=	 (int *) my_malloc(ps.maxdim*sizeof(int));
  participants=  (int *) my_malloc(ps.maxdim*sizeof(int));
  initclasses =  (int *) my_malloc(ps.maxitem*sizeof(int));

  fprintf(stdout,"D \n");
  /* initial read of ps.graphname to get ps.maxobjtuples */
  edgecounts =  (int ***) my_malloc(ps.maxrel*sizeof(int **));
  for (i = 0; i < ps.maxrel; i++) {
    edgecounts[i] =  (int **) my_malloc(ps.maxdim*sizeof(int *));
    for (j = 0; j < ps.maxdim; j++) {
      edgecounts[i][j] =  (int *) my_malloc(ps.maxitem*sizeof(int));
      for (k = 0; k < ps.maxitem; k++) {
        edgecounts[i][j][k] = 0;
      }
    }
  }
  ps.maxobjtuples = 0;

  fprintf(stdout,"D2 \n");
  fileptr = fopen(ps.graphname,"r");
  if (fileptr == NULL) {
    fprintf(stderr, "couldn't read graph\n"); exit(1);
  }
  while( fscanf( fileptr, " %d", &r)!=EOF ) {
    fprintf(stdout,"%s %d %d\n",__FILE__,__LINE__,r);
    ndim = relsizes[r];
    fprintf(stdout,"%s %d %d\n",__FILE__,__LINE__,ndim);
    for (dim = 0; dim < ndim; dim++) {
      fscanf(fileptr, "%d", &participant);
      participants[dim] = participant;
    }
    fscanf(fileptr, "%lf", &val);

    for (dim = 0; dim < ndim; dim++) {
      fprintf(stdout,"D2 %d %d %d \n",r,dim,participants[dim]);
        edgecounts[r][dim][participants[dim]]++;
      fprintf(stdout,"D2 %d %d %d \n",r,dim,participants[dim]);
    }
  }
  fprintf(stdout,"E\n");
  fclose(fileptr);
  for (i = 0; i < ps.maxrel; i++) {
    for (j = 0; j < ps.maxdim; j++) {
      for (k = 0; k < ps.maxitem; k++) {
        if (edgecounts[i][j][k] > ps.maxobjtuples) {
          ps.maxobjtuples = edgecounts[i][j][k];
        }
        edgecounts[i][j][k]= 0;
      }
    }
  }

  fprintf(stdout,"F\n");
  free(relsizes);
  for (i = 0; i < ps.maxrel; i++) {
    for (j = 0; j < ps.maxdim; j++) {
      free(edgecounts[i][j]);
    }
    free(edgecounts[i]);
  }
  free(edgecounts);

  fprintf(stdout,"G\n");
  /* second read of ps.configfile where we set up datastructures */

  fileptr = fopen(ps.configfile,"r");
  if (ps.outsideinit) {
    initzsfile= fopen(ps.initfile,"r");
    if (initzsfile == NULL) {
      fprintf(stderr, "couldn't read initzsfile\n"); exit(1);
    }
  } else {
    initzsfile = NULL;
  }

  fprintf(stdout,"H\n");
  fscanf(fileptr, "%s", prefix);
  fscanf(fileptr, "%d %d", &ndom, &nreln);

  cc = chaincolln_create(nchains, ndom, nreln, prefix);
  c0 = chaincolln_getchain(cc, 0);

  fprintf(stdout,"I\n");
  /* read domains */
  /* input file: nitem maxclass initclass clusterflag*/
  for (d = 0; d < ndom; d++) {
    fscanf(fileptr, "%d %d %d %d", &nitem, &maxclass, &initclass, &clusterflag);
#ifdef GSL
    N = nitem;
#endif
    if (ps.outsideinit) {
      for (zind = 0; zind < nitem; zind++) {
        fscanf(initzsfile, "%d", &initclasses[zind]);
      }
    }
  fprintf(stdout,"J\n");

    /* add domains and items to chains */
    for (cind = 0; cind < nchains; cind++) {
      c = chaincolln_getchain(cc, cind);
      chain_adddomain(c, d, nitem, maxclass, clusterflag, ps.alpha,
		      ps.alphahyp, initclasses);
#ifdef GSL
      perm =  gsl_permutation_alloc(N);
      gsl_permutation_init(perm);
      gsl_ran_shuffle(rng, perm->data, N, sizeof(size_t));
#endif
      /* assign items to classes */
      relcl = 0;
      for (i = 0; i < nitem; i++) {
        if (ps.outsideinit) {
	  chain_additemtoclass(c, d, i, initclasses[i]);
	} else {
          if (relcl == initclass) relcl = 0;

	  /* without the GNUSL, each chain gets initialized the same way. This
	   * is suboptimal */
	  itemind = i;
#ifdef GSL
          itemind = gsl_permutation_get(perm, i);
#endif
          chain_additemtoclass(c, d, itemind, relcl);
          relcl++;
        }
      }
#ifdef GSL
      gsl_permutation_free(perm);
#endif
    }
  }
#ifdef GSL
  gsl_rng_free(rng);
#endif

  fprintf(stdout,"K\n");
  /* read relations*/
  /* input file: ndim d0 ... dn */

  for (r = 0; r < nreln; r++) {
    fscanf(fileptr, "%d", &ndim);
    for (dim=0; dim < ndim; dim++) {
      fscanf(fileptr, "%d", &domlabel);
      domlabels[dim] = domlabel;
    }
    for (cind = 0; cind < nchains; cind++) {
      c = chaincolln_getchain(cc, cind);
      chain_addrelation(c, r, ndim, ps.betaprop, ps.betamag, nig, domlabels);
    }
  }
  if (ps.outsideinit) {
    fclose(initzsfile);
  }

  fprintf(stdout,"L\n");
  fclose(fileptr);
  /* second read of ps.graphname: store edges*/
  fileptr = fopen(ps.graphname,"r");
  /* input file: relind p0 p1 p2 .. pn val */
  while( fscanf( fileptr, " %d", &r)!= EOF ) {
    ndim = relation_getdim( chain_getrelation(c0, r) );
    doms = relation_getdoms( chain_getrelation(c0, r) );
    for (dim = 0; dim < ndim; dim++) {
      fscanf(fileptr, "%d", &participant);
      fprintf(stdout,"M %d %d\n",dim,participant);
      participants[dim] = participant;
      domlabels[dim] = domain_getlabel(doms[dim]);
    }

    for (i = 0; i < ndim; i++) {
      for (j = 0; j < i; j++) {
        if (participants[i] == participants[j] &&
	    domlabels[i] == domlabels[j]) {
	  fprintf(stderr, "Self links not allowed.\n"); exit(1);
	}
      }
    }

    fscanf(fileptr, "%lf", &val);
      fprintf(stderr,"%d\n",nchains);
    for (cind = 0; cind < nchains; cind++) {
      c = chaincolln_getchain(cc, cind);
      chain_addedge(c, r, val, participants);

      rn = chain_getrelation(c, r);

      if (doubleeq(val, 0)) {
	relation_setmissing(rn, 1);
      }

      if (val > 1.5 && relation_getdtype(rn) != CONT) {
	relation_setdtype(rn, FREQ);
      }

      if (!doubleeq(val, (int) val)) {
	relation_setdtype(rn, CONT);
	relation_setmissing(rn, 1); /* XXX: no sparse continuous matrices */
      }

    }
  }

  fprintf(stderr,"N\n");
  fclose(fileptr);

  for (cind = 0; cind < nchains; cind++) {
    c = chaincolln_getchain(cc, cind);
    for (i = 0; i < chain_getndomains(c); i++) {
      chain_updatedomprobs(c, i);
    }
  }

  fprintf(stderr,"O\n");
  free(domlabels); free(participants); free(initclasses);

  return cc;
}

int chaincolln_getnchain(chaincolln cc) {
  return cc->nchains;
}

chain chaincolln_getchain(chaincolln cc, int chainind) {
  return cc->chains[chainind];
}

char *chaincolln_getprefix(chaincolln cc) {
  return cc->prefix;
}

int chaincolln_getitercount(chaincolln cc) {
  return cc->itercount;
}

void chaincolln_print(chaincolln cc) {
  chaincolln_printassignments(cc);
  chaincolln_printstatus(cc, NULL);
  chaincolln_printstatus(cc, stderr);
}

/* run hillclimbing with random restarts */
void chaincolln_climb(chaincolln cc, int maxloops) {
  int i, j, changeflag, nrepeats;
  chain chn, schn;
  double oldscore, newscore, currscore;

  chn = chaincolln_getchain(cc, 0);
  schn = cc->sparechain;
  changeflag = 1;
  nrepeats = 0;
  oldscore = -MAXDOUBLE;
  currscore= -MAXDOUBLE;

  chaincolln_print(cc);
  for (i = 0; i < maxloops; i++) {
    if (changeflag > 0) {
      /* maximizing scan: move each object to the best class for it */
      changeflag = chain_climbscan(chn);
    }
    /*chain_print(chn); fprintf(stderr, "splitting\n");*/
    changeflag = changeflag + chain_climbsplitfast(chn, schn);

    /*chain_print(chn); fprintf(stderr, "merging\n");*/
    changeflag = changeflag + chain_climbmergefast(chn, schn);

    for (j = 0; j < ps.hypupdates; j++) {
      changeflag = changeflag + chaincolln_hypupdates(cc, &proposal_best);
    }

    cc->itercount++;
    newscore = chain_getprob(chn);
    if ( fabs(newscore - oldscore) < 0.00001 * fabs(newscore) ) {
      nrepeats++;
      if (nrepeats == 8) {
        chaincolln_print(cc);
        fprintf(stderr, "*************** randomize\n");
        nrepeats = 0;
	chain_randomize(cc->chains[0]);
        currscore = chain_getprob(chn);
        newscore = chain_getprob(chn);
      }
    } else {
      nrepeats = 0;
    }
    oldscore = newscore;
    chaincolln_printstatus(cc, stderr);
  }
}

/* run gibbs sampler with Metropolis updates (MC^3, split-merge) */
void chaincolln_mcmc(chaincolln cc, int maxloops) {
  int i, j;
  chaincolln_print(cc);
  for (i = 0; i < maxloops; i++) {
    fprintf(stderr,"loop %d",i);
    /* gibbs scans and split merge updates */
    chaincolln_sample(cc, 1);
    /* MC^3: try swaps between chains at different temperatures */
    chaincolln_tryswaps(cc);
    /* hyperparameter updates? */
    for (j = 0; j < ps.hypupdates; j++) {
      chaincolln_hypupdates(cc, &proposal_mcmc);
    }
    chaincolln_print(cc);
  }
}

/* update hyperparameters */
int chaincolln_hypupdates(chaincolln cc,
	int (*proposal_choose) (double, double, double)) {
  int i, changeflag;
  changeflag = 0;
  for (i = 0; i < cc->nchains; i++) {
    chain_hypupdate(chaincolln_getchain(cc, i), proposal_choose, &changeflag);
  }
  return changeflag;
}

void chaincolln_itergibbsm(chaincolln cc) {
  int i;
  if (1) { /* regular gibbs scans */
    for (i = 0; i < cc->nchains; i++) {
      chain_itergibbs(chaincolln_getchain(cc, i));
    }
  }
  if (1) { /* split-merge updates */
    for (i = 0; i < cc->nchains; i++) {
      chain_itersm(chaincolln_getchain(cc, i), cc->sparechain);
    }
  }
  cc->itercount++;
}

void chaincolln_sample(chaincolln cc, int niter) {
  int i;
  for (i = 0; i < niter; i++) {
    chaincolln_itergibbsm(cc);
  }
}

/* try swapping chains at different temperatures */
void chaincolln_tryswaps(chaincolln cc) {
  int nchains;
  int sw1, sw2, swcount;
  chain c1, c2;
  double t1, t2, aswap, rn;

  nchains = cc->nchains;
  if (nchains > 1) {
    for (swcount = 0; swcount < 1; swcount++) {
      sw1 = randomitem(nchains);
      sw2 = sw1;
      while (sw2 == sw1) { sw2 = randomitem(nchains); }

      c1 = cc->chains[sw1]; c2 = cc->chains[sw2];
      t1 = chain_gettemp(c1); t2 = chain_gettemp(c2);
      aswap = exp( (t1 - t2)* (chain_getprob(c2) - chain_getprob(c1)) );
      if (aswap > 1) { aswap=1; }
      rn = myrand();
      if (rn < aswap) {
        chaincolln_swapchains(cc, sw1, sw2);
        fprintf(stderr, "****CS: %f\n", aswap);
      }
    }
  }
}

/* swap chains SW1 and SW2 */
void chaincolln_swapchains(chaincolln cc, int sw1, int sw2) {
  chain c1, c2;
  double t1, t2;
  c1 = chaincolln_getchain(cc, sw1);
  c2 = chaincolln_getchain(cc, sw2);
  t1 = chain_gettemp(c1); t2 = chain_gettemp(c2);
  chain_settemp(c1, t2); chain_settemp(c2, t1);
  cc->chains[sw1] = c2;
  cc->chains[sw2] = c1;
}

void chaincolln_printassignments(chaincolln cc) {
  int i;
  chain chn;
  char filename[MAXSTRING];
  FILE * domfileptr;

  chn = chaincolln_getchain(cc, 0);
  for (i = 0; i < chain_getndomains(chn); i++) {
    sprintf(filename, "%s_dom%d", ps.outroot, i);
    domfileptr = fopen(filename, "a");
    if (domfileptr == NULL) {
      fprintf(stderr, "couldn't read domain file\n");
      exit(1);
    }
    domain_printassignments(domfileptr, chain_getdomain(chn, i));
    fclose(domfileptr);
  }
}

void chaincolln_printstatus(chaincolln cc, FILE *fptr) {
  int i, nchains;
  chain chn;
  char filename[MAXSTRING];
  FILE * domfileptr;

  nchains = chaincolln_getnchain(cc);
  if (fptr == NULL) {
    sprintf(filename, "%s_status", ps.outroot);
    domfileptr = fopen(filename, "a");
    if (domfileptr == NULL) {
      fprintf(stderr, "couldn't read statusfile\n");
      exit(1);
    }
  } else {
    domfileptr = fptr;
  }
  fprintf(domfileptr, "%4d: ", chaincolln_getitercount(cc));
  chn = chaincolln_getchain(cc, 0);
  chain_printdetailedstatus(chn, domfileptr);
  for (i = 1; i < nchains; i++) {
    chn = chaincolln_getchain(cc, i);
    chain_printbriefstatus(chn, domfileptr);
  }
  chn = chaincolln_getchain(cc,0);
  chain_printhyps(chn, domfileptr);
  fprintf(domfileptr, "\n");

  if (fptr == NULL) {
    fclose(domfileptr);
  }
}

void chaincolln_free(chaincolln cc)
{
  int i;
  /* include spare chain */
  for (i = 0; i < cc->nchains+1; i++) {
    chain_free(cc->chains[i]);
  }
  free(cc->chains);
  free(cc);
}

ext: c

max_stars_repo_name: adamsumm/CausalMario

max_issues_repo_name: adamsumm/CausalMario

max_forks_repo_name: adamsumm/CausalMario


text:
#ifndef _GSLRNG_H
#define _GSLRNG_H

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <ctime>
#include <vector>
using namespace std;

namespace myGSL {
  class Rng;
  class PoissonGenerator;
};

typedef myGSL::Rng GSLRng;

class myGSL::Rng {
public:
  Rng(const gsl_rng_type * T = gsl_rng_taus2, unsigned long seed = time(NULL))
    : rng(gsl_rng_alloc(T)), discrete_t(NULL)
  { gsl_rng_set(rng,seed); }

  Rng(const Rng& R)
    : rng(gsl_rng_clone(R.constptr())), discrete_t(NULL)
  {}

  ~Rng() {
    gsl_rng_free(rng);
    if (discrete_t != NULL) discrete_free();
  }

  inline void set(unsigned long seed) { gsl_rng_set(rng,seed); }
  inline unsigned long min() { return gsl_rng_min(rng); }
  inline unsigned long max() { return gsl_rng_max(rng); }

  // uniform distributions
  inline unsigned long get() { return gsl_rng_get(rng); }
  inline double uniform() { return gsl_rng_uniform(rng); }
  inline double uniform_pos() { return gsl_rng_uniform_pos(rng); }
  inline double uniform(double a, double b) { return a + (b-a)*uniform(); }
  inline unsigned long uniform_int(unsigned long n) { return gsl_rng_uniform_int(rng,n); }

  // other distributions
  inline double gaussian(double mu, double s) { return gsl_ran_gaussian(rng,s) + mu; }
  inline unsigned long binomial(double r, unsigned n) { return gsl_ran_binomial(rng,r,n); }
  inline unsigned int poisson(double mu) { return gsl_ran_poisson(rng,mu); }
  inline double exponential(double mu) { return gsl_ran_exponential(rng,mu); }

  // discrete distributions
  inline void discrete_preproc(int K, const double * P) { discrete_t = gsl_ran_discrete_preproc(K,P); }
  inline int  get_discrete() { return gsl_ran_discrete(rng,discrete_t); }
  inline void discrete_free() { gsl_ran_discrete_free(discrete_t); discrete_t = NULL; }

  template<typename T> inline void shuffle(vector<T>& src, size_t first = 0,
      size_t last = -1) {
    if (first >= src.size()) return;
    if (last > src.size()) last = src.size();
    gsl_ran_shuffle(rng,&src[first],last-first,sizeof(T));
  }

  template<typename T> inline void choose(vector<T>& src, vector<T>& dest) {
    gsl_ran_choose(rng,&dest[0],dest.size(),&src[0],src.size(),sizeof(T));
  }

  inline gsl_rng * ptr() { return rng; }
  inline const gsl_rng * constptr() const { return rng; }

private:
  gsl_rng * rng;
  gsl_ran_discrete_t * discrete_t;
};

class myGSL::PoissonGenerator {
  public:
    double mu;
    Rng* rng;
    PoissonGenerator(double m, Rng* r) { mu = m; rng = r; }
    double operator()() { return rng->poisson(mu); }
};

#endif

ext: h


text:
#include <gsl/gsl_errno.h>
#include <gsl/matrix/gsl_matrix.h>
#include <string.h>

#define BASE_DOUBLE
#include <gsl/templates_on.h>
#include <gsl/matrix/copy_source.c>
#include <gsl/templates_off.h>
#undef BASE_DOUBLE

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:
#ifndef MATRIX_H
#define MATRIX_H

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>

#include <config.h>

#include "collections.h"
#include "file_utils.h"
#include "vector.h"
#include "vector_math.h"

#ifdef HAVE_CBLAS
#include <cblas.h>
#else
#warning "No CLBAS"
#endif

typedef enum {
    MATRIX_DENSE,
    MATRIX_SPARSE
} matrix_type_t;

#define MATRIX_INIT(name, type, type_name, array_type)                                                          \
    typedef struct {                                                                                \
        size_t m, n;                                                                                \
        type *values;                                                                                \
    } name##_t;                                                                                \
                                                                                \
    static name##_t *name##_new(size_t m, size_t n) {                                                           \
        name##_t *matrix = malloc(sizeof(name##_t));                                                            \
                                                                                \
        if (matrix == NULL) {                                                                                \
            return NULL;                                                                                \
        }                                                                                \
                                                                                \
        matrix->m = m;                                                                                \
        matrix->n = n;                                                                                \
                                                                                \
        matrix->values = malloc(sizeof(type) * m * n);                                                          \
        if (matrix->values == NULL) {                                                                           \
            free(matrix);                                                                                \
            return NULL;                                                                                \
        }                                                                                \
                                                                                \
        return matrix;                                                                                \
                                                                                \
    }                                                                                \
                                                                                \
    static name##_t *name##_new_aligned(size_t m, size_t n, size_t alignment) {                                 \
        name##_t *matrix = malloc(sizeof(name##_t));                                                            \
                                                                                \
        if (matrix == NULL) {                                                                                \
            return NULL;                                                                                \
        }                                                                                \
                                                                                \
        matrix->m = m;                                                                                \
        matrix->n = n;                                                                                \
                                                                                \
        matrix->values = _aligned_malloc(sizeof(type) * m * n, alignment);                                      \
        if (matrix->values == NULL) {                                                                           \
            free(matrix);                                                                                \
            return NULL;                                                                                \
        }                                                                                \
                                                                                \
        return matrix;                                                                                \
                                                                                \
    }                                                                                \
                                                                                \
    static void name##_destroy(name##_t *self) {                                                                \
        if (self == NULL) return;                                                                               \
                                                                                \
        if (self->values != NULL) {                                                                             \
            free(self->values);                                                                                \
        }                                                                                \
                                                                                \
        free(self);                                                                                \
    }                                                                                \
                                                                                \
    static void name##_destroy_aligned(name##_t *self) {                                                        \
        if (self == NULL) return;                                                                               \
                                                                                \
        if (self->values != NULL) {                                                                             \
            _aligned_free(self->values);                                                                        \
        }                                                                                \
                                                                                \
        free(self);                                                                                \
    }                                                                                \
                                                                                \
    static inline void name##_zero(name##_t *self) {                                                            \
        memset(self->values, 0, self->m * self->n * sizeof(type));                                              \
    }                                                                                \
                                                                                \
                                                                                \
    static inline bool name##_resize(name##_t *self, size_t m, size_t n) {                                      \
        if (self == NULL) return false;                                                                         \
                                                                                \
        if (m * n > (self->m * self->n)) {                                                                      \
            type *ptr = realloc(self->values, sizeof(type) * m * n);                                            \
            if (ptr == NULL) {                                                                                \
                return false;                                                                                \
            }                                                                                \
            self->values = ptr;                                                                                \
        }                                                                                \
                                                                                \
        self->m = m;                                                                                \
        self->n = n;                                                                                \
                                                                                \
        return true;                                                                                \
    }                                                                                \
                                                                                \
    static inline bool name##_resize_aligned(name##_t *self, size_t m, size_t n, size_t alignment) {            \
        if (self == NULL) return false;                                                                         \
                                                                                \
        if (m * n > (self->m * self->n)) {                                                                      \
            type *ptr = _aligned_realloc(self->values, sizeof(type) * m * n, alignment);                        \
            if (ptr == NULL) {                                                                                \
                return false;                                                                                \
            }                                                                                \
            self->values = ptr;                                                                                \
        }                                                                                \
                                                                                \
        self->m = m;                                                                                \
        self->n = n;                                                                                \
                                                                                \
        return true;                                                                                \
    }                                                                                \
                                                                                \
    static inline bool name##_resize_fill_zeros(name##_t *self, size_t m, size_t n) {                           \
        size_t old_m = self->m;                                                                                \
        bool ret = name##_resize(self, m, n);                                                                   \
        if (ret && m > old_m) {                                                                                \
            memset(self->values + old_m, 0, (m - old_m) * self->n * sizeof(type));                              \
        }                                                                                \
        return ret;                                                                                \
    }                                                                                \
                                                                                \
    static inline bool name##_resize_aligned_fill_zeros(name##_t *self, size_t m, size_t n, size_t alignment) { \
        size_t old_m = self->m;                                                                                \
        bool ret = name##_resize_aligned(self, m, n, alignment);                                                \
        if (ret && m > old_m) {                                                                                \
            memset(self->values + old_m, 0, (m - old_m) * self->n * sizeof(type));                              \
        }                                                                                \
        return ret;                                                                                \
    }                                                                                \
                                                                                \
    static inline name##_t *name##_new_copy(name##_t *self) {                                                   \
        name##_t *cpy = name##_new(self->m, self->n);                                                           \
        size_t num_values = self->m * self->n;                                                                  \
        memcpy(cpy->values, self->values, num_values * sizeof(type));                                           \
                                                                                \
        return cpy;                                                                                \
    }                                                                                \
                                                                                \
    static inline bool name##_copy(name##_t *self, name##_t *other) {                                           \
        if (self->m != other->m || self->n != other->n) {                                                       \
            return false;                                                                                \
        }                                                                                \
        size_t num_values = self->m * self->n;                                                                  \
                                                                                \
        memcpy(other->values, self->values, num_values * sizeof(type));                                         \
        return true;                                                                                \
    }                                                                                \
                                                                                \
    static inline void name##_init_values(name##_t *self, type *values) {                                       \
        size_t num_values = self->m * self->n;                                                                  \
        memcpy(self->values, values, num_values * sizeof(type));                                                \
    }                                                                                \
                                                                                \
    static inline void name##_set(name##_t *self, type value) {                                                 \
        array_type##_set(self->values, value, self->m * self->n);                                               \
    }                                                                                \
                                                                                \
    static inline void name##_set_row(name##_t *self, size_t index, type *row) {                                \
        size_t offset = index * self->n;                                                                        \
        type *values = self->values;                                                                            \
        size_t n = self->n;                                                                                \
        memcpy(values + offset, row, n * sizeof(type));                                                         \
    }                                                                                \
                                                                                \
    static inline void name##_set_scalar(name##_t *self, size_t row_index, size_t col_index, type value) {      \
        size_t offset = row_index * self->n + col_index;                                                        \
        self->values[offset] = value;                                                                           \
    }                                                                                \
                                                                                \
    static inline void name##_add_scalar(name##_t *self, size_t row_index, size_t col_index, type value) {      \
        size_t offset = row_index * self->n + col_index;                                                        \
        self->values[offset] += value;                                                                          \
    }                                                                                \
                                                                                \
    static inline void name##_sub_scalar(name##_t *self, size_t row_index, size_t col_index, type value) {      \
        size_t offset = row_index * self->n + col_index;                                                        \
        self->values[offset] -= value;                                                                          \
    }                                                                                \
                                                                                \
    static inline void name##_mul_scalar(name##_t *self, size_t row_index, size_t col_index, type value) {      \
        size_t offset = row_index * self->n + col_index;                                                        \
        self->values[offset] *= value;                                                                          \
    }                                                                                \
                                                                                \
    static inline void name##_div_scalar(name##_t *self, size_t row_index, size_t col_index, type value) {      \
        size_t offset = row_index * self->n + col_index;                                                        \
        self->values[offset] /= value;                                                                          \
    }                                                                                \
                                                                                \
    static inline type name##_get(name##_t *self, size_t row_index, size_t col_index) {                         \
        size_t index = row_index * self->n + col_index;                                                         \
        return self->values[index];                                                                             \
    }                                                                                \
                                                                                \
    static inline type *name##_get_row(name##_t *self, size_t row_index) {                                      \
        size_t index = row_index * self->n;                                                                     \
        return self->values + index;                                                                            \
    }                                                                                \
                                                                                \
    static inline name##_t *name##_new_value(size_t m, size_t n, type value) {                                  \
        name##_t *matrix = name##_new(m, n);                                                                    \
        name##_set(matrix, value);                                                                              \
        return matrix;                                                                                \
    }                                                                                \
                                                                                \
    static inline name##_t *name##_new_zeros(size_t m, size_t n) {                                              \
        name##_t *matrix = name##_new(m, n);                                                                    \
        name##_zero(matrix);                                                                                \
        return matrix;                                                                                \
    }                                                                                \
                                                                                \
    static inline name##_t *name##_new_ones(size_t m, size_t n) {                                               \
        return name##_new_value(m, n, (type)1);                                                                 \
    }                                                                                \
                                                                                \
    static inline name##_t *name##_new_values(size_t m, size_t n, type *values) {                               \
        name##_t *matrix = name##_new(m, n);                                                                    \
        memcpy(matrix->values, values, m * n * sizeof(type));                                                   \
        return matrix;                                                                                \
    }                                                                                \
                                                                                \
    static inline void name##_div(name##_t *self, type value) {                                                 \
        array_type##_div(self->values, value, self->m * self->n);                                               \
    }                                                                                \
                                                                                \
    static inline bool name##_div_matrix(name##_t *self, name##_t *other) {                                     \
        if (self->m != other->m || self->n != other->n) return false;                                           \
        array_type##_div_array(self->values, other->values, self->m * self->n);                                 \
        return true;                                                                                \
    }                                                                                \
                                                                                \
    static inline bool name##_div_matrix_times_scalar(name##_t *self, name##_t *other, type v) {                \
        if (self->m != other->m || self->n != other->n) return false;                                           \
        array_type##_div_array_times_scalar(self->values, other->values, v, self->m * self->n);                 \
        return true;                                                                                \
    }                                                                                \
                                                                                \
    static inline void name##_mul(name##_t *self, type value) {                                                 \
        array_type##_mul(self->values, value, self->m * self->n);                                               \
    }                                                                                \
                                                                                \
    static inline bool name##_mul_matrix(name##_t *self, name##_t *other) {                                     \
        if (self->m != other->m || self->n != other->n) return false;                                           \
        array_type##_mul_array(self->values, other->values, self->m * self->n);                                 \
        return true;                                                                                \
    }                                                                                \
                                                                                \
    static inline bool name##_mul_matrix_times_scalar(name##_t *self, name##_t *other, type v) {                \
        if (self->m != other->m || self->n != other->n) return false;                                           \
        array_type##_mul_array_times_scalar(self->values, other->values, v, self->m * self->n);                 \
        return true;                                                                                \
    }                                                                                \
                                                                                \
    static inline void name##_add(name##_t *self, type value) {                                                 \
        array_type##_add(self->values, self->m * self->n, value);                                               \
    }                                                                                \
                                                                                \
                                                                                \
    static inline bool name##_add_matrix(name##_t *self, name##_t *other) {                                     \
        if (self->m != other->m || self->n != other->n) return false;                                           \
        array_type##_add_array(self->values, other->values, self->m * self->n);                                 \
        return true;                                                                                \
    }                                                                                \
                                                                                \
    static inline bool name##_add_matrix_times_scalar(name##_t *self, name##_t *other, type v) {                \
        if (self->m != other->m || self->n != other->n) return false;                                           \
        array_type##_add_array_times_scalar(self->values, other->values, v, self->m * self->n);                 \
        return true;                                                                                \
    }                                                                                \
                                                                                \
    static inline void name##_sub(name##_t *self, type value) {                                                 \
        array_type##_sub(self->values, value, self->m * self->n);                                               \
    }                                                                                \
                                                                                \
    static inline bool name##_sub_matrix(name##_t *self, name##_t *other) {                                     \
        if (self->m != other->m || self->n != other->n) return false;                                           \
        array_type##_sub_array(self->values, other->values, self->m * self->n);                                 \
        return true;                                                                                \
    }                                                                                \
                                                                                \
    static inline bool name##_sub_matrix_times_scalar(name##_t *self, name##_t *other, type v) {                \
        if (self->m != other->m || self->n != other->n) return false;                                           \
        array_type##_sub_array_times_scalar(self->values, other->values, v, self->m * self->n);                 \
        return true;                                                                                \
    }                                                                                \
                                                                                \
    static name##_t *name##_read(FILE *f) {                                                                     \
        name##_t *mat = malloc(sizeof(name##_t));                                                               \
        if (mat == NULL) return NULL;                                                                           \
                                                                                \
        mat->values = NULL;                                                                                \
                                                                                \
        uint64_t m = 0;                                                                                \
        uint64_t n = 0;                                                                                \
                                                                                \
        if (!file_read_uint64(f, &m) ||                                                                         \
            !file_read_uint64(f, &n)) {                                                                         \
            goto exit_##name##_allocated;                                                                       \
        }                                                                                \
                                                                                \
        mat->m = (size_t)m;                                                                                \
        mat->n = (size_t)n;                                                                                \
                                                                                \
        size_t len_data = mat->m * mat->n;                                                                      \
                                                                                \
        type *data = malloc(len_data * sizeof(type));                                                           \
        if (data == NULL) {                                                                                \
            log_error("error in data malloc\n");                                                                \
            goto exit_##name##_allocated;                                                                       \
        }                                                                                \
                                                                                \
        if (!file_read_##array_type(f, data, len_data)) {                                                       \
            free(data);                                                                                \
            goto exit_##name##_allocated;                                                                       \
        }                                                                                \
                                                                                \
        mat->values = data;                                                                                \
                                                                                \
        return mat;                                                                                \
                                                                                \
    exit_##name##_allocated:                                                                                \
        name##_destroy(mat);                                                                                \
        return NULL;                                                                                \
    }                                                                                \
                                                                                \
    static bool name##_write(name##_t *self, FILE *f) {                                                         \
        if (self == NULL || self->values == NULL) {                                                             \
            return false;                                                                                \
        }                                                                                \
                                                                                \
        if (!file_write_uint64(f, (uint64_t)self->m) ||                                                         \
            !file_write_uint64(f, (uint64_t)self->n)) {                                                         \
            return false;                                                                                \
        }                                                                                \
                                                                                \
        uint64_t len_data = (uint64_t)self->m * (uint64_t)self->n;                                              \
                                                                                \
        for (uint64_t i = 0; i < len_data; i++) {                                                               \
            if (!file_write_##type_name(f, self->values[i])) {                                                  \
                return false;                                                                                \
            }                                                                                \
        }                                                                                \
                                                                                \
        return true;                                                                                \
    }

#define MATRIX_INIT_FLOAT_BASE(name, type, type_name, array_type)                       \
    MATRIX_INIT(name, type, type_name, array_type)                                      \
                                                                                \
    static inline void name##_log(name##_t *self) {                                     \
        array_type##_log(self->values, self->m * self->n);                              \
    }                                                                                \
                                                                                \
    static inline void name##_exp(name##_t *self) {                                     \
        array_type##_exp(self->values, self->m * self->n);                              \
    }                                                                                \
                                                                                \
    static inline void name##_dot_vector(name##_t *self, type *vec, type *result) {     \
        type *values = self->values;                                                    \
        size_t m = self->m;                                                             \
        size_t n = self->n;                                                             \
        for (size_t i = 0; i < m; i++) {                                                \
            for (size_t j = 0; j < n; j++) {                                            \
                result[i] += values[n * i + j] * vec[j];                                \
            }                                                                           \
        }                                                                               \
    }

#ifdef HAVE_CBLAS
#define MATRIX_INIT_FLOAT(name, type, type_name, array_type, blas_prefix)                                   \
    MATRIX_INIT_FLOAT_BASE(name, type, type_name, array_type)                                               \
                                                                                \
    static inline bool name##_dot_matrix(name##_t *m1, name##_t *m2, name##_t *result) {                    \
        if (m1->n != m2->m || m1->m != result->m || m2->n != result->n) {                                   \
            return false;                                                                                \
        }                                                                                \
                                                                                \
        log_debug("doing CBLAS\n");                                                                         \
        cblas_##blas_prefix##gemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,                                \
                    m1->m, m2->n, m1->n, 1.0,                                                               \
                    m1->values, m1->n,                                                                      \
                    m2->values, m2->n, 0.0,                                                                 \
                    result->values, result->n                                                               \
                    );                                                                                \
                                                                                \
        return true;                                                                                \
    }

#else
#define MATRIX_INIT_FLOAT(name, type, type_name, array_type, blas_prefix)                                       \
    MATRIX_INIT_FLOAT_BASE(name, type, type_name, array_type)                                                   \
                                                                                \
    static inline bool name##_dot_matrix(name##_t *m1, name##_t *m2, name##_t *result) {                        \
        if (m1->n != m2->m || m1->m != result->m || m2->n != result->n) {                                       \
            return false;                                                                                \
        }                                                                                \
                                                                                \
        size_t m1_rows = m1->m;                                                                                \
        size_t m1_cols = m1->n;                                                                                \
        size_t m2_rows = m2->m;                                                                                \
        size_t m2_cols = m2->n;                                                                                \
                                                                                \
        type *m1_values = m1->values;                                                                           \
        type *m2_values = m2->values;                                                                           \
        type *result_values = result->values;                                                                   \
                                                                                \
        for (size_t i = 0; i < m1_rows; i++) {                                                                  \
            for (size_t j = 0; j < m2_cols; j++) {                                                              \
                size_t result_index = m2_cols * i + j;                                                          \
                result_values[result_index] = 0.0;                                                              \
                for (size_t k = 0; k < m2_rows; k++) {                                                          \
                    result_values[result_index] += m1_values[m1_cols * i + k] * m2_values[m2_cols * k + j];     \
                }                                                                                \
            }                                                                                \
        }                                                                                \
                                                                                \
        return true;                                                                                \
    }
#endif

MATRIX_INIT(uint32_matrix, uint32_t, uint32, uint32_array)

MATRIX_INIT_FLOAT(float_matrix, float, float, float_array,s)
MATRIX_INIT_FLOAT(double_matrix, double, double, double_array,d)

#endif

max_stars_repo_name: rinigus/deb-libpostal

max_issues_repo_name: rinigus/deb-libpostal

max_forks_repo_name: rinigus/deb-libpostal


text:
#include <gsl/gsl_sf_bessel.h>
#include <stdio.h>

/**
 * The example program from the GSL documentation
 * http://www.gnu.org/software/gsl/doc/html/usage.html#an-example-program
 */

int main() {
    double x = 5.0;
    double y = gsl_sf_bessel_J0 (x);
    printf ("J0(%g) = %.18e\n", x, y);
    return 0;
}

max_stars_repo_name: rockandsalt/conan-center-index

max_issues_repo_name: rockandsalt/conan-center-index

max_forks_repo_name: rockandsalt/conan-center-index


text:
/* Copyright 2019, Alistair Boyle, 3-clause BSD License */
#include "config.h"
#include <stdlib.h> /* malloc, free */
#include <assert.h> /* assert */
#include <stdio.h> /* printf */
#include <string.h> /* memset */
#include <float.h> /* DBL_MIN */
#include <lapacke.h>

#include "matrix.h"
#include "fwd.h"
#include "inv.h"

#ifdef UNIT_TESTING
extern void * _mock_test_malloc(const size_t size, const char * file, const int line);
extern void _test_free(void * const ptr, const char * file, const int line);
#define malloc(size) _mock_test_malloc(size, __FILE__, __LINE__)
#define free(ptr) _test_free(ptr, __FILE__, __LINE__)
#endif

#define SUCCESS 1
#define FAILURE 0

// #define DEBUG_SVD_MATRICES 1

#ifdef DEBUG_SVD_MATRICES
static void set_matrix(matrix * X, char symbol[], char name[], char units[], int m, int n, double * x)
{
    memset(X, 0, sizeof(matrix));
    X->m = m;
    X->n = n;
    X->type = DENSE;
    X->symmetric = false;
    X->scale = 1.0;
    X->symbol = &(symbol[0]);
    X->name = &(name[0]);
    X->units = &(units[0]);
    X->x.dense = x;
}
#endif

/*
static int svd(matrix * J, matrix * U, matrix * sv, matrix * Vt)
{
    assert(J != NULL);
    const int m = J->m;
    const int n = J->n;
    return m*n;
}
*/

/* max/min: watch out for double evaluation of arguments e.g. a++ == a+=2 */
#define min(a,b) (((a) < (b)) ? (a) : (b) )
#define max(a,b) (((a) > (b)) ? (a) : (b) )
/* dense matrix indexing */
#define colmaj(i,j,ld) ((i) + (j)*(ld)) /* column major indexing */

/* solve Δx=(JᵀJ+λ²I)⁻¹Jb for b=b₂-b₁
 * where J=UΣVᵀ (by the SVD)
 * and
 * Δx=VΣ⁺Uᵀb=VϕΣ⁻¹Uᵀb with Σ⁺=(Σ²+λ²I)⁻¹Σ or ϕₖ=σₖ²/(σₖ²+λ²)
 * for filter factor ϕ
 */
#define phi(s,hp) ((s*s)/((s*s)+(hp*hp)))
static int inv_solve_svd(model * const mdl, double const * b, double * x)
{
    int ret = FAILURE;
    assert(mdl != NULL);
    assert(b != NULL);
    assert(x != NULL);
    const int m = mdl->n_stimmeas;
    const int n = mdl->n_params[0];
    /* init J */
    // printf("J is %d x %d = %d --> %d B\n", m, n, m * n, m * n * (int) sizeof(double));
    double * Ux = NULL, * sx = NULL, * Vtx = NULL, * superb = NULL;
    double * Jx = malloc(sizeof(double) * m * n);
    if(Jx == NULL) {
        printf("error: %s: out of memory\n", "jacobian");
        return FAILURE;
    }
    if(!calc_jacobian(mdl, &Jx[0])) {
        printf("error: inv_solve/calc_jacobian failed\n");
        free(Jx);
        return FAILURE;
    }
#ifdef DEBUG_SVD_MATRICES
    matrix J;
    set_matrix(&J, "J", "Jacobian", "Vm/S", m, n, &Jx[0]);
    printf_matrix(&J);
#endif
    /* calculate SVD of J; J is destroyed in the process */
    /* SVD config */
    const int layout = LAPACK_COL_MAJOR;
    /* A= m cols of U --> u; S= min(m,n) cols U --> u; O=min(m,n) cols U --> a (overwrite); 'N'=no U */
    const char jobu = 'A';
    const char jobvt = 'A';
    const int ldJ = (layout == LAPACK_ROW_MAJOR) ? n : m;
    const int ldU = ((jobu == 'S') && (layout == LAPACK_ROW_MAJOR)) ? min(m, n) : m;
    const int ldVt = ((jobu == 'S') && (layout == LAPACK_COL_MAJOR)) ? min(m, n) : n;
    /* SVD outputs */
    const int Usz = ((jobu == 'S') && (layout == LAPACK_COL_MAJOR)) ? ldU * min(m, n) : ldU * m;
    const int ssz = min(m, n);
    const int Vtsz = ((jobu == 'S') && (layout == LAPACK_COL_MAJOR)) ? ldVt * min(m, n) : ldVt * n;
    const int sbsz = min(m, n) - 1;
    Ux = malloc(sizeof(double) * Usz);
    sx = malloc(sizeof(double) * ssz);
    Vtx = malloc(sizeof(double) * Vtsz);
    superb = malloc(sizeof(double) * sbsz);
    assert(Ux != NULL);
    assert(sx != NULL);
    assert(Vtx != NULL);
    assert(superb != NULL);
    int info = LAPACKE_dgesvd(layout, jobu, jobvt,
                              m, n, Jx, ldJ,
                              sx,
                              Ux, ldU,
                              Vtx, ldVt,
                              superb);
    //printf("LAPACKE_dgesvd info = %d\n", info);
    assert(info == 0);
    free(Jx);
    free(superb);
#if DEBUG_SVD_MATRICES
    matrix U, s, Vt;
    set_matrix(&U, "U", "left singular vectors", NULL, m, Usz / ldU, &Ux[0]);
    set_matrix(&s, "Σ", "singular values", NULL, ssz, 1, &sx[0]);
    set_matrix(&Vt, "Vᵀ", "right singular vectors", NULL, n, Vtsz / ldVt, &Vtx[0]);
    printf_matrix(&U);
    printf_matrix(&s);
    printf_matrix(&Vt);
#endif
    /* calculate Utb = Uᵀb */
    double Utb[ssz];
    memset(Utb, 0, sizeof(double)*ssz);
    for(int i = 0; i < m; i++) { /* i: b row = Uᵀ col = U row */
        for(int j = 0; j < ssz; j++) { /* j: Uᵀb row = Uᵀ row = U col */
            Utb[j] += Ux[colmaj(i, j, ldU)] * b[i];
        }
    }
#if DEBUG_SVD_MATRICES
    matrix Utbm;
    set_matrix(&Utbm, "Uᵀb", NULL, NULL, ssz, 1, &Utb[0]);
    printf_matrix(&Utbm);
#endif
    /* calculate Utb = ϕΣ⁻¹(Uᵀb) */
    const double hp = mdl->hp;
    for(int j = 0; j < ssz; j++) {
        Utb[j] = (phi(sx[j], hp) * Utb[j]) / (sx[j] + DBL_EPSILON); /* filtered Uᵀb */
    }
#ifdef DEBUG_SVD_MATRICES
    for(int j = 0; j < ssz; j++) {
        printf("ϕ[%d]=%8g Utb[%d]=%8g\n", j, phi(sx[j], hp), j, Utb[j]);
    }
#endif
    /* complete the solution as Δx=V(ϕΣ⁻¹Uᵀb) */
    memset(x, 0, sizeof(double) * n);
    for(int i = 0; i < ssz; i++) { /* i: Uᵀb row = V col = Vᵀ row */
        const double Utbi = Utb[i];
        for(int j = 0; j < n; j++) { /* j: x row = V row = Vᵀ col */
            x[j] += Vtx[colmaj(i, j, ldVt)] * Utbi;
        }
    }
#ifdef DEBUG_SVD_MATRICES
    matrix xm;
    set_matrix(&xm, "x", NULL, NULL, n, 1, &x[0]);
    printf_matrix(&xm);
#endif
    free(Ux);
    free(sx);
    free(Vtx);
    ret = SUCCESS;
    return ret;
}

static int inv_solve_diff(model * const mdl, double const * b1, double const * b2, double * x)
{
    assert(mdl != NULL);
    assert(b1 != NULL);
    assert(b2 != NULL);
    const int m = mdl->n_stimmeas;
    double b[m];
    for(int i = 0; i < m; i++) {
        b[i] = b2[i] - b1[i];
    }
#ifdef DEBUG_SVD_MATRICES
    printf("b = ");
    for(int i = 0; i < m; i++) {
        printf(" %g", b[i]);
    }
    printf("\n");
#endif
    return inv_solve_svd(mdl, &b[0], x);
}

int inv_solve(model * const mdl, double * x)
{
    assert(mdl != NULL);
    assert(x != NULL);
    const int m = mdl->n_stimmeas;
    const int n = mdl->n_params[0];
    for(int frame_idx = 1; frame_idx < mdl->n_data[1]; frame_idx++) {
        int ret = inv_solve_diff(mdl, &(mdl->data[0]), &(mdl->data[frame_idx * m]), &x[(frame_idx - 1) * n]);
        if(!ret) {
            return FAILURE;
        }
        if(mdl->n_params[1] > 1) {
            break;
        }
    }
    // printf_model(mdl);
    return SUCCESS;
}

max_stars_repo_name: boyle/zedhat

max_issues_repo_name: boyle/zedhat

max_forks_repo_name: boyle/zedhat


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>

static void
initialise (double * ri, double * rj, double * rg, double * rh,
            double alpha, double beta);

gsl_integration_qaws_table *
gsl_integration_qaws_table_alloc (double alpha, double beta, int mu, int nu)
{
  gsl_integration_qaws_table * t;

  if (alpha < -1.0)
    {
      GSL_ERROR_VAL ("alpha must be greater than -1.0", GSL_EINVAL, 0);
    }

  if (beta < -1.0)
    {
      GSL_ERROR_VAL ("beta must be greater than -1.0", GSL_EINVAL, 0);
    }

  if (mu != 0 && mu != 1)
    {
      GSL_ERROR_VAL ("mu must be 0 or 1", GSL_EINVAL, 0);
    }

  if (nu != 0 && nu != 1)
    {
      GSL_ERROR_VAL ("nu must be 0 or 1", GSL_EINVAL, 0);
    }

  t = (gsl_integration_qaws_table *)
    malloc(sizeof(gsl_integration_qaws_table));

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

  t->alpha = alpha;
  t->beta = beta;
  t->mu = mu;
  t->nu = nu;

  initialise (t->ri, t->rj, t->rg, t->rh, alpha, beta);

  return t;
}

int
gsl_integration_qaws_table_set (gsl_integration_qaws_table * t,
                                double alpha, double beta, int mu, int nu)
{
  if (alpha < -1.0)
    {
      GSL_ERROR ("alpha must be greater than -1.0", GSL_EINVAL);
    }

  if (beta < -1.0)
    {
      GSL_ERROR ("beta must be greater than -1.0", GSL_EINVAL);
    }

  if (mu != 0 && mu != 1)
    {
      GSL_ERROR ("mu must be 0 or 1", GSL_EINVAL);
    }

  if (nu != 0 && nu != 1)
    {
      GSL_ERROR ("nu must be 0 or 1", GSL_EINVAL);
    }

  t->alpha = alpha;
  t->beta = beta;
  t->mu = mu;
  t->nu = nu;

  initialise (t->ri, t->rj, t->rg, t->rh, alpha, beta);

  return GSL_SUCCESS;
}

void
gsl_integration_qaws_table_free (gsl_integration_qaws_table * t)
{
  RETURN_IF_NULL (t);
  free (t);
}

static void
initialise (double * ri, double * rj, double * rg, double * rh,
            double alpha, double beta)
{
  const double alpha_p1 = alpha + 1.0;
  const double beta_p1 = beta + 1.0;

  const double alpha_p2 = alpha + 2.0;
  const double beta_p2 = beta + 2.0;

  const double r_alpha = pow (2.0, alpha_p1);
  const double r_beta = pow (2.0, beta_p1);

  size_t i;

  double an, anm1;

  ri[0] = r_alpha / alpha_p1;
  ri[1] = ri[0] * alpha / alpha_p2;

  an = 2.0;
  anm1 = 1.0;

  for (i = 2; i < 25; i++)
    {
      ri[i] = -(r_alpha + an * (an - alpha_p2) * ri[i - 1])
        / (anm1 * (an + alpha_p1));
      anm1 = an;
      an = an + 1.0;
    }

  rj[0] = r_beta / beta_p1;
  rj[1] = rj[0] * beta / beta_p2;

  an = 2.0;
  anm1 = 1.0;

  for (i = 2; i < 25; i++)
    {
      rj[i] = -(r_beta + an * (an - beta_p2) * rj[i - 1])
        / (anm1 * (an + beta_p1));
      anm1 = an;
      an = an + 1.0;
    }

  rg[0] = -ri[0] / alpha_p1;
  rg[1] = -rg[0] - 2.0 * r_alpha / (alpha_p2 * alpha_p2);

  an = 2.0;
  anm1 = 1.0;

  for (i = 2; i < 25; i++)
    {
      rg[i] = -(an * (an - alpha_p2) * rg[i - 1] - an * ri[i - 1]
                + anm1 * ri[i]) / (anm1 * (an + alpha_p1));
      anm1 = an;
      an = an + 1.0;
    }

  rh[0] = -rj[0] / beta_p1;
  rh[1] = -rh[0] - 2.0 * r_beta / (beta_p2 * beta_p2);

  an = 2.0;
  anm1 = 1.0;

  for (i = 2; i < 25; i++)
    {
      rh[i] = -(an * (an - beta_p2) * rh[i - 1] - an * rj[i - 1]
                + anm1 * rj[i]) / (anm1 * (an + beta_p1));
      anm1 = an;
      an = an + 1.0;
    }

  for (i = 1; i < 25; i += 2)
    {
      rj[i] *= -1;
      rh[i] *= -1;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <gsl.h>

#include <mitkBaseData.h>

namespace crimson
{
class QtPropertyStorage;

/*! \brief   A solver parameters data interface. */
class ISolverParametersData : public mitk::BaseData
{
public:
    mitkClassMacro(ISolverParametersData, BaseData);

    /*!
    * \brief   Gets property storage for the SolverParametersData properties.
    */
    virtual gsl::not_null<QtPropertyStorage*> getPropertyStorage() = 0;

protected:
    ISolverParametersData() { mitk::BaseData::InitializeTimeGeometry(1); }
    virtual ~ISolverParametersData() {}

    ISolverParametersData(const Self&) = default;

    void SetRequestedRegion(const itk::DataObject*) override {}
    void SetRequestedRegionToLargestPossibleRegion() override {}
    bool RequestedRegionIsOutsideOfTheBufferedRegion() override { return true; }
    bool VerifyRequestedRegion() override { return true; }
};
}

max_stars_repo_name: carthurs/CRIMSONGUI

max_issues_repo_name: carthurs/CRIMSONGUI

max_forks_repo_name: carthurs/CRIMSONGUI


text:
/*
 * Copyright 2014 RWTH Aachen University. All rights reserved.
 *
 * Licensed under the RWTH LM License (the "License");
 * you may not use this file except in compliance with the License.
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
#include <gsl/gsl_cblas.h>
#include <cassert>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <numeric>

typedef double Real;

inline void FastZero(const int size, float x[]) {
  std::fill(x, x + size, 0.0f);
}

inline void FastZero(const int size, double x[]) {
  std::fill(x, x + size, 0.0);
}

inline void FastCopy(const float source[],
                     const int size,
                     float destination[]) {
  memcpy(destination, source, size * sizeof(float));
}

inline void FastCopy(const double source[],
                     const int size,
                     double destination[]) {
  memcpy(destination, source, size * sizeof(double));
}

inline float FastMax(const float x[], const int size) {
  return *std::max_element(x, x + size);
}

inline double FastMax(const double x[], const int size) {
  return *std::max_element(x, x + size);
}

inline float FastComputeSum(const float x[], const int size) {
  return std::accumulate(x, x + size, 0.0f);
}

inline double FastComputeSum(const double x[], const int size) {
  return std::accumulate(x, x + size, 0.0);
}

inline void FastAddConstant(const float source[],
                            const int size,
                            const float value,
                            float destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [value](const float x) { return x + value; });
}

inline void FastAddConstant(const double source[],
                            const int size,
                            const double value,
                            double destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [value](const double x) { return x + value; });
}

inline void FastSubtractConstant(const float source[],
                                 const int size,
                                 const float value,
                                 float destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [value](const double x) { return x - value; });
}

inline void FastSubtractConstant(const double source[],
                                 const int size,
                                 const double value,
                                 double destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [value](const double x) { return x - value; });
}

inline void FastReverseSubtractConstant(
    const float source[],
    const int size,
    const float value,
    float destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [value](const float x) { return value - x; });
}

inline void FastReverseSubtractConstant(
    const double source[],
    const int size,
    const double value,
    double destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [value](const double x) { return value - x; });
}

inline void FastMultiplyByConstant(const float source[],
                                   const int size,
                                   const float value,
                                   float destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [value](const float x) { return value * x; });
}

inline void FastMultiplyByConstant(const double source[],
                                   const int size,
                                   const double value,
                                   double destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [value](const double x) { return value * x; });
}

inline void FastDivideByConstant(const float source[],
                                 const int size,
                                 const float value,
                                 float destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [value](const float x) { return x / value; });
}

inline void FastDivideByConstant(const double source[],
                                 const int size,
                                 const double value,
                                 double destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [value](const double x) { return x / value; });
}

inline void FastInvert(const float source[],
                       const int size,
                       float destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [](const double x) { return 1.0f / x; });
}

inline void FastInvert(const double source[],
                       const int size,
                       double destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [](const double x) { return 1.0 / x; });
}

inline void FastAdd(const float a[],
                    const int size,
                    const float b[],
                    float c[]) {
  std::transform(a,
                 a + size,
                 b,
                 c,
                 [](const float x, const float y) { return x + y; });
}

inline void FastAdd(const double a[],
                    const int size,
                    const double b[],
                    double c[]) {
  std::transform(a,
                 a + size,
                 b,
                 c,
                 [](const double x, const double y) { return x + y; });
}

inline void FastMultiplyByConstantAdd(const float alpha,
                                      const float a[],
                                      const int size,
                                      float b[]) {
  cblas_saxpy(size, alpha, a, 1, b, 1);
}

inline void FastMultiplyByConstantAdd(const double alpha,
                                      const double a[],
                                      const int size,
                                      double b[]) {
  cblas_daxpy(size, alpha, a, 1, b, 1);
}

inline void FastSub(const float a[],
                    const int size,
                    const float b[],
                    float c[]) {
  std::transform(a,
                 a + size,
                 b,
                 c,
                 [](const float x, const float y) { return x - y; });
}

inline void FastSub(const double a[],
                    const int size,
                    const double b[],
                    double c[]) {
  std::transform(a,
                 a + size,
                 b,
                 c,
                 [](const double x, const double y) { return x - y; });
}

inline void FastMultiply(const float a[],
                         const int size,
                         const float b[],
                         float c[]) {
  std::transform(a,
                 a + size,
                 b,
                 c,
                 [](const float x, const float y) { return x * y; });
}

inline void FastMultiply(const double a[],
                         const int size,
                         const double b[],
                         double c[]) {
  std::transform(a,
                 a + size,
                 b,
                 c,
                 [](const double x, const double y) { return x * y; });
}

inline void FastMultiplyAdd(const float a[],
                            const int size,
                            const float b[],
                            float c[]) {
  for (int i = 0; i < size; ++i)
    c[i] += a[i] * b[i];
}

inline void FastMultiplyAdd(const double a[],
                            const int size,
                            const double b[],
                            double c[]) {
  for (int i = 0; i < size; ++i)
    c[i] += a[i] * b[i];
}

inline void FastTanh(const float source[],
                     const int size,
                     float destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [](const float x) { return std::tanh(x); });
}

inline void FastTanh(const double source[],
                     const int size,
                     double destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [](const double x) { return std::tanh(x); });
}

inline void FastExponential(const float source[],
                            const int size,
                            float destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [](const float x) { return std::exp(x); });
}

inline void FastExponential(const double source[],
                            const int size,
                            double destination[]) {
  std::transform(source,
                 source + size,
                 destination,
                 [](const double x) { return std::exp(x); });
}

inline void FastMatrixVectorMultiply(const float a[],
                                     const bool transpose_a,
                                     const int rows_a,
                                     const int columns_a,
                                     const float x[],
                                     float y[]) {
  cblas_sgemv(CblasColMajor,
              transpose_a ? CblasTrans : CblasNoTrans,
              rows_a,
              columns_a,
              1.0f,
              a,
              rows_a,
              x,
              1,
              1.0f,
              y,
              1);
}

inline void FastMatrixVectorMultiply(const double a[],
                                     const bool transpose_a,
                                     const int rows_a,
                                     const int columns_a,
                                     const double x[],
                                     double y[]) {
  cblas_dgemv(CblasColMajor,
              transpose_a ? CblasTrans : CblasNoTrans,
              rows_a,
              columns_a,
              1.0,
              a,
              rows_a,
              x,
              1,
              1.0,
              y,
              1);
}

inline void FastOuterProduct(const float alpha,
                             const float x[],
                             const int size_x,
                             const float y[],
                             const int size_y,
                             float a[]) {
  cblas_sger(CblasColMajor,
             size_x,
             size_y,
             alpha,
             x,
             1,
             y,
             1,
             a,
             size_x);
}

inline void FastOuterProduct(const double alpha,
                             const double x[],
                             const int size_x,
                             const double y[],
                             const int size_y,
                             double a[]) {
  cblas_dger(CblasColMajor,
             size_x,
             size_y,
             alpha,
             x,
             1,
             y,
             1,
             a,
             size_x);
}

inline void FastMatrixMatrixMultiply(
      const float alpha,
      const float a[],
      const bool transpose_a,
      const int rows_a,
      const int columns_a,
      const float b[],
      const bool transpose_b,
      const int columns_b,
      float c[]) {
  cblas_sgemm(CblasColMajor,
              transpose_a ? CblasTrans : CblasNoTrans,
              transpose_b ? CblasTrans : CblasNoTrans,
              rows_a,
              columns_b,
              columns_a,
              alpha,
              a,
              transpose_a ? columns_a : rows_a,
              b,
              transpose_b ? columns_b : columns_a,
              1.0f,
              c,
              rows_a);
}

inline void FastMatrixMatrixMultiply(
      const double alpha,
      const double a[],
      const bool transpose_a,
      const int rows_a,  // m
      const int columns_a,  // k
      const double b[],
      const bool transpose_b,
      const int columns_b,  // n
      double c[]) {
  cblas_dgemm(CblasColMajor,
              transpose_a ? CblasTrans : CblasNoTrans,
              transpose_b ? CblasTrans : CblasNoTrans,
              rows_a,  // m
              columns_b,  // n
              columns_a,  // k
              alpha,
              a,
              transpose_a ? columns_a : rows_a,
              b,
              transpose_b ? columns_b : columns_a,
              1.0,
              c,
              rows_a);
}

inline Real *FastMalloc(const int size) {
  Real *result = new Real[size];
  assert(result != nullptr || size == 0);
  return result;
}

inline void FastFree(Real *x) {
  delete [] x;
}


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_errno.h>
#include <time.h>

int main (int argc, char* argv[])
{
  /* Parse the command line arguments. */
  const char* usage = "Usage: zeroes <xmin> <xmax> <ymin> <ymax> <xres> <coeff> ... <coeff> - <degree> ... <degree>\n";

  if (argc < 8) {
    fprintf (stderr, "%s", usage);
    exit(EXIT_FAILURE);
  }

  /* The window in the plane. */
  const double xmin = atof(argv[1]);
  const double xmax = atof(argv[2]);
  const double ymin = atof(argv[3]);
  const double ymax = atof(argv[4]);

  /* Image size, width is given, height is computed. */
  const int xres = atoi(argv[5]);
  const int yres = (int)((xres * (ymax-ymin)) / (xmax-xmin));

  /* The coefficients */
  int ci_max = 0 ;
  while (6 + ci_max < argc && (argv[6+ci_max][0] != '-' || argv[6+ci_max][1] != 0)) { ci_max++; }
  if (ci_max == 0) {
    fprintf (stderr, "%s\nError: specify at least one coefficient.", usage) ;
    exit(EXIT_FAILURE);
  }
  double coeff[ci_max];
  for (int ci = 0; ci < ci_max; ci++) { coeff[ci] = atof(argv[6+ci]); }

  /* The degrees */
  const int di_max = argc - 7 - ci_max;
  if (di_max <= 0) {
    fprintf (stderr, "%s\nError: specify at least one degree.", usage) ;
    exit(EXIT_FAILURE);
  }
  int degree[di_max];
  for (int di = 0; di < di_max; di++) { degree[di] = atoi(argv[7+ci_max+di]); }

  /* Allocate the image and initialize it */
  unsigned int *image = (unsigned int *)calloc(sizeof(unsigned int), yres * xres) ;
  const unsigned int count_bound = 0xffffffff;
  unsigned int max_count = 0 ;

  /* Turn off the GSL error handler that aborts on error (some polynomials have tricky zeroes) */
  gsl_set_error_handler_off ();

  /* We're going to skip zeroes whose imaginary part is almost zero. */
  const double epsilon = 1.0e-20 ;

  time_t start = time(NULL);

  /* Run through degrees and compute the zeroes for each one. */
  for (int di = 0; di < di_max; di++) {
    const int d = degree[di]; /* The degree we're working on. */
    time_t now = time(NULL);
    fprintf (stderr, "Degree %d (%.2lf seconds)\n", d, difftime(now, start));
    /* Initialize the polynomial. */
    double poly[d+1];
    for (int j = 0; j <= d; j++) { poly[j] = coeff[0]; }
    /* Initialize the counters. */
    int counter[d+1];
    for (int j = 0; j <= d; j++) { counter[j] = 0; }
    int j = 0;
    /* Allocate the workspace. */
    gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc(d+1);
    do {
      /* Compute the zeroes if head coefficient is non-zero. */
      if (poly[d] > epsilon || poly[d] < -epsilon) {
        double z[2*d]; /* the zeroes are stored here */
        int status = gsl_poly_complex_solve (poly, d+1, w, z);
        /* Only use zeroes if GSL reported success */
        if (status == 0) {
          /* Draw zeroes, skipping the real ones that are away from the origin */
          for (int i=0; i < d; i++) {
            if (-epsilon < z[2*i+1] && z[2*i+1] < epsilon && (z[2*i] < -0.5 || z[2*i] > 0.5)) {
              continue;
            }
            int x = (int)((xres * (z[2*i] - xmin)) / (xmax - xmin));
            int y = yres - (int)((yres * (z[2*i+1] - ymin)) / (ymax - ymin));
            if (0 <= x && x < xres && 0 <= y && y < yres && image[xres * y + x] < count_bound) {
              int c = ++image[xres * y + x];
              if (max_count < c) { max_count = c; }
            }
          }
        }
      }
      /* calculate the next polynomial */
      for (j = 0; j <= d && counter[j] == ci_max-1; j++) { counter[j] = 0; poly[j] = coeff[0]; }
      if (j <= d) { counter[j]++; poly[j] = coeff[counter[j]]; }
    } while (j <= d);
    /* Deallocate workspace. */
    gsl_poly_complex_workspace_free (w);
  }

  /* Output the computed result. */
  /* ASCII header to the file*/
  printf("P6\n# Zeroes, xmin=%lf, xmax=%lf, ymin=%lf, ymax=%lf, max_count=%d, coeffs=[", xmin, xmax, ymin, ymax, max_count);
  for (int ci=0; ci < ci_max; ci++) { printf("%lf%s", coeff[ci], (ci < ci_max-1 ? ", " : "],")); }
  printf("degrees = [");
  for (int di=0; di < di_max; di++) { printf("%d%s", degree[di], (di < di_max-1 ? ", " : "],")); }
  printf("\n");
  printf("%d\n%d\n%d\n", xres, yres, 0xffff);
  /* Image */
  double r = 65535.0 / log((double)max_count) ;
  for (int y = 0; y < yres; y++) {
    for (int x = 0; x < xres; x++) {
        unsigned char color[6];
        unsigned int k = (unsigned int)(r * log((float)image[xres * y + x])) ;
        color[0] = k >> 8;
        color[1] = k & 255;
        color[2] = k >> 8;
        color[3] = k & 255;
        color[4] = k >> 8;
        color[5] = k & 255;
        fwrite(color, 6, 1, stdout);
    }
  }
  free(image);
  return(0);
}

max_stars_repo_name: andrejbauer/zeroes

max_issues_repo_name: andrejbauer/zeroes

max_forks_repo_name: andrejbauer/zeroes


text:
#include <math.h>
#include <gsl/gsl_sf_gamma.h>
#include <galpy_potentials.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
//PowerSphericalPotentialwCutoff
//3  arguments: amp, alpha, rc
double mass(double r2,double alpha, double rc){
  return 2. * M_PI * pow ( rc , 3. - alpha ) * ( gsl_sf_gamma ( 1.5 - 0.5 * alpha ) - gsl_sf_gamma_inc ( 1.5 - 0.5 * alpha , r2 / rc / rc ) );
}
double PowerSphericalPotentialwCutoffEval(double R,double Z, double phi,
					  double t,
					  struct potentialArg * potentialArgs){
  double * args= potentialArgs->args;
  //Get args
  double amp= *args++;
  double alpha= *args++;
  double rc= *args;
  //Radius
  double r2= R*R+Z*Z;
  double r= sqrt(r2);
  return amp * 2. * M_PI * pow(rc,3.-alpha) / r * ( r / rc * ( gsl_sf_gamma ( 1. - 0.5 * alpha ) - gsl_sf_gamma_inc ( 1. - 0.5 * alpha , r2 / rc / rc ) ) - ( gsl_sf_gamma ( 1.5 - 0.5 * alpha ) - gsl_sf_gamma_inc ( 1.5 - 0.5 * alpha , r2 / rc / rc) ) );
}
double PowerSphericalPotentialwCutoffRforce(double R,double Z, double phi,
					    double t,
					    struct potentialArg * potentialArgs){
  double * args= potentialArgs->args;
  //Get args
  double amp= *args++;
  double alpha= *args++;
  double rc= *args;
  //Radius
  double r2= R*R+Z*Z;
  //Calculate Rforce
  return - amp * mass (r2,alpha,rc) * R / pow(r2,1.5);
}
double PowerSphericalPotentialwCutoffPlanarRforce(double R,double phi,
						  double t,
						  struct potentialArg * potentialArgs){
  double * args= potentialArgs->args;
  //Get args
  double amp= *args++;
  double alpha= *args++;
  double rc= *args;
  //Radius
  double r2= R*R;
  //Calculate Rforce
  return - amp * mass (r2,alpha,rc) / r2;
}
double PowerSphericalPotentialwCutoffzforce(double R,double Z,double phi,
					    double t,
					    struct potentialArg * potentialArgs){
  double * args= potentialArgs->args;
  //Get args
  double amp= *args++;
  double alpha= *args++;
  double rc= *args;
  //Radius
  double r2= R*R+Z*Z;
  //Calculate Rforce
  return - amp * mass (r2,alpha,rc) * Z / pow(r2,1.5);
}
double PowerSphericalPotentialwCutoffPlanarR2deriv(double R,double phi,
						   double t,
						   struct potentialArg * potentialArgs){
  double * args= potentialArgs->args;
  //Get args
  double amp= *args++;
  double alpha= *args++;
  double rc= *args;
  //Radius
  double r2= R*R;
  //Calculate R2deriv
  return amp * ( 4. * M_PI * pow(r2,- 0.5 * alpha) * exp(-r2/rc/rc) - 2. * mass(r2,alpha,rc)/pow(r2,1.5) );
}

max_stars_repo_name: turnergarrow/galpy

max_issues_repo_name: BurcuAkbulut/galpy

max_forks_repo_name: BurcuAkbulut/galpy