text:
//changes to version glm_popdyn_1.0_test2:
//lookuptable for G(x)=exp(-0.5*DT*c*exp(x/deltaV)) instead of lambda and Plam

#ifndef GLM
#define GLM 0
#endif

#ifndef GLM2
#define GLM2 2
#endif

#ifndef GLIF
#define GLIF 10
#endif

#ifndef GLIF1
#define GLIF1 11
#endif

#ifndef GLIF2
#define GLIF2 12
#endif

#ifndef GLIF4
#define GLIF4 14
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <complex.h>
#include <fftw3.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include "util.c"

#define maxisi 2.
#define PFIRE_MAX 0.99999
#define PFIRE_MIN 0.00001

double extime;
unsigned int ntrial;

struct PopParametersGLM{
  double tref;
  double taum;
  double *taus1; //incoming synaptic time constants, decay time of 1st filter
  double *taus2; //incoming synaptic time constants, decay time of 2nd filter
  double *taur1; //incoming synaptic time constants, rise time of 1st filter
  double *taur2; //incoming synaptic time constants, rise time of 2nd filter
  double *a1;    //relative weight of fast synaptic current
  double *a2;    //relative weight of slow synaptic current
  double mu;
  double c;
  double deltaV;
  double delay;
  double vth;
  double u_reset;      //reset potential for u in glif mode
  int N;         //number of neurons in population
  double *J;     //incoming synaptic weights
  double *Iext;  //external input (in mV), i.e. Iext is actually R*I
  double *tau_theta;
  double *J_theta; //amplitudes of the N_theta exponential kernels in mV
  int N_theta;  //number of exponentials for dynamic threshold

  //internal constants
  int K;
  int end;
  int indx;
  double Em;
  double Em2;
  double *E_theta;  //Jtheta*exp(-t/tau_theta)};
  double *E_theta2;  //exp(-dt/tau_theta)};
  double *gamma_QR;  //deltaV*(1-exp(-gamma/deltaV))/N; gamma=sum_k Jtheta_k*exp(-t/tautheta_k)
  double *w1;          // effective weight
  double *w2;
  double g;
  double g2;
  double dV;  // resolution of voltage grid for lookup table of hazard/firing prob
  double Vmin;       //min value of voltage grid
  double Vmax;       //max value of voltage grid (this value is not included, i.e. max is Vmax-dv)
  double *Gtable;      //lookup table for firing prob P=1-exp(-lam*DT)
};

struct PopVariablesGLM{
  double h;
  double *u;
  double *Xsyn1; //incoming synaptic currents (filtered with fast decay time, e.g. AMPA, GABA)
  double *Xsyn2; //incoming synaptic currents (filtered with slow decay time, e.g. NMDA)
  double *Isyn1; //total incoming synaptic current (sum of Xsyn1+Xsyn2 filtered with short rise-time)
  double *Isyn2; //total incoming synaptic current (sum of Xsyn1+Xsyn2 filtered with short rise-time)
  unsigned int *n;
  double *m0;
  double *theta;  //threshold kernel
  double *g; //QR threshold variables
  double *v;
  double z;
  double x;
  double Gfree;
  double *G;
};

double DT;
double DTBIN;

////////////////////////////////////////////////////////////////////////////////

double get_mean(double *x,int n);
void print_matrix1(int N,int M, double A[][M]);
void print_matrix2(double **A,int N,int M);
void print_matrix3(int N,int M, double *A[N]);
void print_pop_parameters(struct PopParametersGLM p[], int Npop);
unsigned long int random_seed();

////////////////////////////////////////////////////////////////////////////////

void init_synaptic_filters(double **Es1, double **Es2, double **Er1, double **Er2, struct PopParametersGLM p[], int Npop)
{
  int i,j;
  for (i=0;i<Npop;i++)
    for (j=0;j<Npop;j++)
      {
	if (p[i].taur1[j]>0) Er1[i][j]=exp(-DT/p[i].taur1[j]);
	else Er1[i][j]=0;
	if (p[i].taus1[j]>0) Es1[i][j]=exp(-DT/p[i].taus1[j]);
	else Es1[i][j]=0;

	if (p[i].taur2[j]>0) Er2[i][j]=exp(-DT/p[i].taur2[j]);
	else Er2[i][j]=0;
	if (p[i].taus2[j]>0) Es2[i][j]=exp(-DT/p[i].taus2[j]);
	else Es2[i][j]=0;
      }
}

void get_inputs(double *input, int k,double **Es1,double **Es2,double **Er1,double **Er2, struct PopParametersGLM p[], struct PopVariablesGLM pop[],int Npop)
{
  int i,j,ndelay;

  double A[Npop];
  for (j=0;j<Npop;j++)
    {
      ndelay=p[j].delay/DT;
//      A[j]= pop[j].n[(p[j].end-1-ndelay+p[j].K)%p[j].K] / (p[j].N * DT);
      A[j]= pop[j].n[(p[j].end-ndelay+p[j].K)%p[j].K] / (p[j].N * DT);
    }

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

      //base line input
      input[i]=p[i].mu;

      //external input, assumed to be slow compared to DT
      if (p[i].Iext != NULL)
	input[i] += p[i].Iext[k] * (1-p[i].Em);  //add external input if provided

      //synaptic input
      for (j=0;j<Npop;j++)
	{
	  if (p[i].taur1[j]>0)
	    {
	      printf("Warning: rise time not correctly implemented\n");
	      input[i]+=p[i].w1[j]*pop[i].Isyn1[j];

	      if (p[i].taus1[j] > 0)
		{
		  pop[i].Isyn1[j] = pop[i].Xsyn1[j] + (pop[i].Isyn1[j] - pop[i].Xsyn1[j]) * Er1[i][j];
		  pop[i].Xsyn1[j] = A[j] + (pop[i].Xsyn1[j] - A[j]) * Es1[i][j];
		}
	      else
		  pop[i].Isyn1[j] = A[j] + (pop[i].Isyn1[j] - A[j]) * Es1[i][j];
	    }
	  else //zero-rise time
	    if (p[i].taus1[j] > 0)
	      {
		input[i] += p[i].w1[j] * (A[j] + (p[i].taus1[j] * (pop[i].Xsyn1[j] - A[j]) * Es1[i][j] - (pop[i].Xsyn1[j] * p[i].taus1[j] - A[j] * p[i].taum) * p[i].Em) / (p[i].taus1[j] - p[i].taum));

		pop[i].Xsyn1[j] = A[j] + (pop[i].Xsyn1[j] - A[j]) * Es1[i][j];
	      }
	    else //delta current
	      {
		input[i]+=p[i].w1[j]*A[j] * (1 - p[i].Em);
	      }

	  //same for 2nd filter
	  if (p[i].taur2[j]>0)
	    {
	      input[i]+=p[i].w2[j]*pop[i].Isyn2[j];

	      if (p[i].taus2[j] > 0)
		{
		  pop[i].Isyn2[j] = pop[i].Xsyn2[j] + (pop[i].Isyn2[j] - pop[i].Xsyn2[j]) * Er2[i][j];
		  pop[i].Xsyn2[j] = A[j] + (pop[i].Xsyn2[j] - A[j]) * Es2[i][j];
		}
	      else
		  pop[i].Isyn2[j] = A[j] + (pop[i].Isyn2[j] - A[j]) * Es2[i][j];
	    }
	  else //zero-rise time
	    if (p[i].taus2[j] > 0)
	      {
		input[i]+=p[i].w2[j]*pop[i].Xsyn2[j];
		pop[i].Xsyn2[j] = A[j] + (pop[i].Xsyn2[j] - A[j]) * Es2[i][j];
	      }
	    else //delta current
	      input[i]+=p[i].w2[j]*A[j];
	}

    }
}

double cond_rate(double x,struct PopParametersGLM *p)
{
  if (p->deltaV > 0) return p->c * exp(x / p->deltaV);
  else //p->D==0, no noise
    {
      if (x < 0) return 0;
      else return 1./DT;
    }
}

static double get_Plam(double V, struct PopParametersGLM *p, double *G_old)
{
  if (V<p->Vmax)
    {
      if (V>=p->Vmin)
	{
	  //int x=(int)((V - p->Vmin) / p->dV);
	  double G = p->Gtable[(int)((V - p->Vmin) / p->dV)];
	  double Plam = 1 - *G_old * G;
	  //printf("V=%g x=%d lam=%g\n",V,x,lam);
	  *G_old = G;
	  return Plam;
	}
      else
	{
	  *G_old=1;
	  return 0.;
	}
    }
  else
    {
      *G_old = 0;
      return 1.;
    }
}

void update(struct PopVariablesGLM *pop, double *nmean, double input,gsl_rng *rng,struct PopParametersGLM *p, int mode)
{

  int i,k,r;
  int kref=(int)(p->tref/DT);
  double W=0, X=0, Y=0, Z=0, theta=p->vth;
  int end=p->end, K=p->K;
  double threshold;

  // first, determine non-refractory spiking prob.
  for (r=0;r < p->N_theta; r++)
    {
      pop->g[r] = pop->g[r]*p->E_theta2[r] + (1 - p->E_theta2[r]) * pop->n[end]/(p->N * DT);
      theta += p->E_theta[r] * pop->g[r];
    }

  double Pfree;
  pop->h = input + (pop->h - p->mu) * p->Em;
  Pfree=get_Plam(pop->h-theta,p,&(pop->Gfree));

  for (k=0; k<K; k++) X += pop->m0[k];

  double ztot,Plam,PLAM, mm;
  theta -= pop->n[end] * p->gamma_QR[0];//substracts effect of n[end] above, added again in 1st loop below, but ensures that threshold theta(t+DT|k=end) is correct

  switch(mode)
    {
    case GLM:
      for (k=end;k<=end-kref+K;k++)
	{
	  i=k%K;
	  threshold = pop->theta[k-end] + theta;
	  theta += pop->n[i] * p->gamma_QR[k-end];
	  Plam=get_Plam(pop->h-threshold,p,pop->G+i);

	  Z += pop->v[i];
	  mm = Plam * pop->m0[i];
	  Y += Plam * pop->v[i];
	  W += mm;
	  pop->v[i] = (1-Plam)*(1-Plam)*pop->v[i] + mm;
	  pop->m0[i] -= mm;
	}
      break;

    case GLIF:
      {
      for (k=end;k<=end-kref+K;k++)
	{
	  i=k%K;
	  threshold = pop->theta[k-end] + theta;
	  theta += pop->n[i] * p->gamma_QR[k-end];
	  pop->u[i] = input + (pop->u[i] - p->mu) * p->Em;//GLIF
	  Plam=get_Plam(pop->u[i]-threshold,p,pop->G+i);

	  Z += pop->v[i];
	  mm = Plam * pop->m0[i];
	  Y += Plam * pop->v[i];
	  W += mm;
	  pop->v[i] = (1-Plam)*(1-Plam)*pop->v[i] + mm;
	  pop->m0[i] -= mm;

	}
      }
      pop->u[end] = p->u_reset;//GLIF
      break;
    }

  ztot = Z + pop->z;
  if (ztot>0) PLAM=(Y+Pfree*pop->z)/ztot;
  else PLAM=0;
  //  printf("W=%g Z=%g z=%g Y=%g y=%g\n",W, Z,pop->z,Y,Pfree*pop->z);

  *nmean = W + Pfree * pop->x + PLAM * (p->N - X - pop->x);
  pop->n[end] = gsl_ran_binomial(rng, *nmean/p->N, p->N);

  pop->z=(1-Pfree)*(1-Pfree)*pop->z+Pfree*pop->x+pop->v[end];
  pop->x=(1-Pfree)*pop->x+pop->m0[end];

  pop->m0[end] = pop->n[end];
  pop->v[end] = 0.;
  pop->G[end] = 1.;

}

void simulate(int Nbin,double **A,double **rate,struct PopVariablesGLM pop[],gsl_rng *rng,struct PopParametersGLM p[],int Npop, int mode, int dispprog)
{
  double input[Npop],normfact[Npop], nmean[Npop], nmean_bin[Npop];
  int n_bin[Npop];
  int L,i;
  int l=0,k=0;
  L=(int)(DTBIN/DT+0.5);
  for (i=0;i<Npop;i++)
    {
      n_bin[i]=0;
      nmean_bin[i]=0;
      normfact[i]=1./(L*DT*p[i].N);
    }

  double **Es1=dmatrix(Npop,Npop);
  double **Es2=dmatrix(Npop,Npop);
  double **Er1=dmatrix(Npop,Npop);
  double **Er2=dmatrix(Npop,Npop);
  init_synaptic_filters(Es1,Es2,Er1,Er2,p,Npop);

  int dispcount=(int)(Nbin/100);

  //  clock_t start=clock();
  while (k<Nbin)
    {
      //      printf("Nbin=%d, k=%d l=%d\n",Nbin,k,l);
      get_inputs(input,k,Es1,Es2,Er1,Er2,p,pop,Npop);

      for (i=0;i<Npop;i++)
	{
	  update(pop+i,nmean+i,input[i],rng,p+i,mode);
      	  n_bin[i]+=pop[i].n[p[i].end];
      	  nmean_bin[i]+=nmean[i];
   	}
      l++;
      //           printf("nmeanbin=%g\n",nmean_bin[0]);
	   //      printf("\n");
      if (l>=L)
      	{
      	  for (i=0;i<Npop;i++)
	    {
	      A[i][k] = n_bin[i] * normfact[i];
	      rate[i][k] = nmean_bin[i] * normfact[i];

	      //if (i==0) printf("k=%d rate=%g\n",k,rate[i][k]);
	      //	      aa[i][k]=a[i];
      	      n_bin[i]=0;
	      nmean_bin[i]=0;
	    }
      	  l=0;
      	  k++;
      	}

      for (i=0;i<Npop;i++) p[i].end=(p[i].end + 1) % p[i].K;

      if (dispprog)
      	if ((k+1)%dispcount==0)
      	  {
      	    printf("%d%%  \r",(int)((k+1) * 100. / Nbin));
      	    fflush(stdout);
      	  }

    }

  //  printf("Execution time of mesoscopic sim: %g seconds, K=%d\n",(double)(clock()-start)/CLOCKS_PER_SEC,p[0].K);

  free_dmatrix(Es1);
  free_dmatrix(Es2);
  free_dmatrix(Er1);
  free_dmatrix(Er2);

}

double threshold_kernel(double t, struct PopParametersGLM *p)
{
  double theta=0;
  int r;
  for (r=0; r<p->N_theta;r++)
    if (t>=0)  theta += p->J_theta[r] / p->tau_theta[r] * exp(-t / p->tau_theta[r]);
  return theta;
}

void init_glm(struct PopParametersGLM p[],struct PopVariablesGLM pop[],int Npop, gsl_rng *rng, int mode)
//memory allocation and initialization of variables,
{
  int i,j,k,r, K;
  for (i=0;i<Npop;i++)
    {
      K=p[i].K;
      //      printf("i=%d K=%d\n",i,K);
      pop[i].n=uivector(K);
      pop[i].m0=dvector(K);
      pop[i].v=dvector(K);
      if (mode>=GLIF) pop[i].u=dvector(K);
      else pop[i].u=NULL;
      pop[i].Isyn1=dvector(Npop);
      pop[i].Isyn2=dvector(Npop);
      pop[i].Xsyn1=dvector(Npop);
      pop[i].Xsyn2=dvector(Npop);

      p[i].Em= exp(-DT / p[i].taum);
      p[i].Em2= exp(-0.5*DT / p[i].taum);
      p[i].end = 0;

      pop[i].theta=dvector(K);
      p[i].gamma_QR=dvector(K);
      pop[i].G=dvector(K);

      for (k=0;k<K;k++)
	{
	  /* if (k<K-1)  */
	  /*   { */
	  /*     pop[i].m0[k]=(double)(p[i].N)/(K-1); //normalization (the K-1 element is set zero below) */
	  /*     pop[i].n[k]=(int)(pop[i].m0[k]+1); //initializing with spike count corresponding to zero hazard */

	  /* pop[i].n[k]=(int)(100*DT*p[i].N+1); //initializing with spike count corresponding to 100 Hz */
	  /* pop[i].m0[k]=exp(-(K-k-1)*DT*100)*pop[i].n[k]; */
	  /* n+=pop[i].m0[k]; */
	  //if (i==0) printf("m0=%g n=%d\n", pop[i].m0[k], pop[i].n[k]);
	  /*   } */

	  pop[i].m0[k]=0;
	  pop[i].n[k]=0;
	  pop[i].v[k]=0;
	  if (mode>=GLIF) pop[i].u[k]=p[i].vth-20*gsl_rng_uniform(rng);

	  pop[i].theta[k]=threshold_kernel( (K-k-1) * DT, p+i);//shifted by DT to get filter at theta(s+DT)
	  p[i].gamma_QR[k] = p[i].deltaV * (1 - exp(-pop[i].theta[k] / p[i].deltaV)) / p[i].N;
	  //	  if (i==0) printf("i=%d k=%d gamma_QR=%g\n",i,k,p[i].gamma_QR[k]);

	  pop[i].G[k]=1.;
	}
      pop[i].Gfree=1.;

      pop[i].n[K-1]=p[i].N;
      pop[i].m0[K-1]=p[i].N;
      pop[i].z=0.;
      pop[i].x=0.;
      //pop[i].x=(p[i].N - n) * DT;
      /* printf("i=%d n=%g x/dt=%g\n",i,n,pop[i].x/DT); */

      for (j=0;j<Npop;j++) pop[i].Isyn1[j]=0;
      for (j=0;j<Npop;j++) pop[i].Isyn2[j]=0;
      for (j=0;j<Npop;j++) pop[i].Xsyn1[j]=0;
      for (j=0;j<Npop;j++) pop[i].Xsyn2[j]=0;

      pop[i].h=p[i].mu;

      //pre-calculating constant J_k*exp(-tau_rel/tau_k)
      pop[i].g=dvector(p[i].N_theta);
      for (r=0; r < p[i].N_theta; r++)
	pop[i].g[r] = pop[i].n[0]/(p[i].N * DT)/p[i].tau_theta[r];

      p[i].E_theta=dvector(p[i].N_theta);
      p[i].E_theta2=dvector(p[i].N_theta);
      //      printf("E_theta= ");
      for (r=0; r < p[i].N_theta; r++)
	{
	  p[i].E_theta[r] = p[i].J_theta[r] * exp(- K * DT / p[i].tau_theta[r]);
	  p[i].E_theta2[r] = exp(- DT / p[i].tau_theta[r]);
	  //	  printf("%g tautheta=%g",p[i].E_theta[r],p[i].tau_theta[r]);
	}
      //      printf("\n");

            //pre-computing firing probability as lookup table
      p[i].dV=p[i].deltaV/100;
      p[i].Vmin = p[i].deltaV * log(-log(1-PFIRE_MIN)/DT / p[i].c);
      p[i].Vmax = p[i].deltaV * log(-log(1-PFIRE_MAX)/DT / p[i].c);
      int L=(p[i].Vmax - p[i].Vmin)/p[i].dV;
      //      printf("Lookup tbl (popul %d): Use L=%d, Vmin=%g Vmax=%g dV=%g\n",i+1,L,p[i].Vmin,p[i].Vmax,p[i].dV);
      p[i].Gtable=dvector(L);
      for (k=0;k<L;k++)
	{
	  double V=p[i].Vmin + k*p[i].dV;
	  p[i].Gtable[k]= exp(-0.5 * p[i].c * exp(V/p[i].deltaV) * DT);
	}
    }
}

void get_history_size(struct PopParametersGLM p[],int Npop)
{
  int i;
  double tmax=20.;//sec

  for (i=0;i<Npop;i++)
    {
      int k=tmax/DT;
      int kmin=5*p[i].taum/DT;
      while ((threshold_kernel(k * DT, p+i) / p[i].deltaV < 0.1) && (k>kmin)) k--;
      //      printf("K=%d\n",k);
      if (k*DT<=p[i].delay)
	k=(int)(p[i].delay/DT)+1;
      //      printf("K=%d\n",k);
      if (k*DT<=p[i].tref)
	k=(int)(p[i].tref/DT)+1;
      p[i].K=k;
      //      p[i].K=20000;
      printf("Use K=%d bins for history of population %d (T=%g sec)\n",p[i].K,i+1,p[i].K*DT);
    }
}

void free_pop(struct PopVariablesGLM pop[],struct PopParametersGLM p[],int Npop, int mode)
{
  int i;
  for (i=0;i<Npop;i++)
    {
      free_uivector(pop[i].n);
      free_dvector(pop[i].m0);
      free_dvector(pop[i].v);
      if (mode>=GLIF) free_dvector(pop[i].u);
      free_dvector(pop[i].Isyn1);
      free_dvector(pop[i].Isyn2);
      free_dvector(pop[i].Xsyn1);
      free_dvector(pop[i].Xsyn2);

      free_dvector(pop[i].g);
      free_dvector(pop[i].theta);
      free_dvector(p[i].E_theta);
      free_dvector(p[i].E_theta2);
      free_dvector(p[i].gamma_QR);
      free_dvector(pop[i].G);
      free_dvector(p[i].Gtable);
    }
}

void get_psd_pop(double **SA,int Nbin,int Ntrials,struct PopParametersGLM neuron[],int Npop, int mode)
{
  gsl_rng *rng=gsl_rng_alloc(gsl_rng_taus2);
  unsigned long int seed;
  seed = random_seed();
  gsl_rng_set(rng,seed);
  //  gsl_rng_set (rng,(long)time(NULL));

  get_history_size(neuron,Npop);

  int i,j,n;

  struct PopVariablesGLM pop[Npop];
  init_glm(neuron,pop,Npop, rng, mode);
  double sum=0;
  for (i=0;i<neuron[0].K;i++)
    {
      sum+=pop[0].m0[i];
    }

  double **A, **rate;
  A=dmatrix(Npop,Nbin);
  rate=dmatrix(Npop,Nbin);

  simulate(Nbin,A,rate,pop,rng,neuron,Npop, mode, 0);

  double complex *AF[Npop];
  for (j=0;j<Npop;j++)
    {
      AF[j]=(double complex *)malloc(sizeof(double complex)*Nbin);
      SA[j]=(double *)malloc(sizeof(double)*(Nbin/2));//psd of A
      for (i=0;i<Nbin/2;i++) SA[j][i]=0;
    }

  fftw_plan plan=fftw_plan_dft_r2c_1d(Nbin,A[0],AF[0],FFTW_ESTIMATE);

  for (n=0;n<Ntrials;n++)
    {
      //      if (n%10==9)
	{
           printf("trial %d ",n+1);
           for (j=0;j<Npop;j++) printf("%g ",get_mean(A[j],Nbin));
           printf("\r");
	   fflush(stdout);
	}
	simulate(Nbin,A,rate,pop,rng,neuron,Npop,mode, 0);

      for (j=0;j<Npop;j++)
	{
	  fftw_execute_dft_r2c(plan,A[j],AF[j]);
	  for (i=1;i<Nbin/2+1;i++) SA[j][i-1]+=creal(AF[j][i]*conj(AF[j][i]))*DTBIN/Nbin;
	}
    }
  printf("trial %d ",n);
  for (j=0;j<Npop;j++) printf("%g ",get_mean(A[j],Nbin));
  printf("\n");

  for (j=0;j<Npop;j++)
    for (i=0;i<Nbin/2;i++)  SA[j][i]/=Ntrials;

  gsl_rng_free (rng);
  free_pop(pop,neuron,Npop, mode);

  free_dmatrix(A);
  free_dmatrix(rate);

  for (j=0;j<Npop;j++) free(AF[j]);
  fftw_destroy_plan(plan);
}

void get_trajectory(double **A,double **rate,int Nbin, int Npop,struct PopParametersGLM neuron[],double dt, double dtbin, int mode, int seed)
{

  gsl_rng *rng=gsl_rng_alloc(gsl_rng_taus2);
  gsl_rng_set (rng,(long)seed);
  /* printf("seed=%d\n",seed); */
  /* printf("rndnr=%g\n",gsl_rng_uniform(rng)); */
  DT=dt;
  DTBIN=dtbin;

  get_history_size(neuron,Npop);

  struct PopVariablesGLM pop[Npop];
  init_glm(neuron,pop,Npop, rng, mode);

  clock_t start=clock();
  simulate(Nbin,A,rate,pop,rng,neuron,Npop, mode, 1);
  double sim_t=(double)(clock()-start)/CLOCKS_PER_SEC;
  printf("Execution time of mesoscopic dynamics: %g seconds, %g s per biosecond\n",sim_t,sim_t/Nbin/DTBIN);

  int i;
  //for (i=0;i<20;i++) printf("rate=%g\n",rate[0][i]);
  for (i=0;i<Npop;i++) printf("mean rate of population %d: %g \n",i+1,get_mean(rate[i],Nbin));
  //  for (i=0;i<Nbin;i++) printf("%g\n",rate[0][i]);

  gsl_rng_free (rng);
  free_pop(pop,neuron,Npop, mode);
}

void init_population(struct PopParametersGLM p[], int Npop,double tref[], double taum[], double taus1[][Npop], double taus2[][Npop], double taur1[][Npop], double taur2[][Npop], double a1[][Npop], double a2[][Npop], double mu[], double c[], double deltaV[], double delay[], double vth[], double vreset[], int N[], double J[][Npop], double p_conn[][Npop], double **signal, int N_theta[], double J_ref[], double *J_theta[], double *tau_theta[], double sigma[])
{

  int k,l;
  for (k=0;k<Npop;k++)
    {
      p[k].indx=k;
      p[k].tref=tref[k];
      p[k].taum=taum[k];
      p[k].mu=mu[k];
      p[k].c=c[k];
      p[k].deltaV=deltaV[k];
      p[k].delay=delay[k];
      p[k].vth=vth[k];
      p[k].u_reset = vreset[k];
      p[k].N=N[k];
      p[k].g=sqrt(1-exp(-2*DT/p[k].taum)) * sigma[k];
      p[k].g2=sqrt(1-exp(-DT/p[k].taum)) * sigma[k];
      p[k].Iext=signal[k];
      p[k].taus1=taus1[k];
      p[k].taus2=taus2[k];
      p[k].taur1=taur1[k];
      p[k].taur2=taur2[k];
      if (J_ref[k]==0)
	p[k].N_theta=N_theta[k];
      else
	p[k].N_theta=N_theta[k]+1; //putting exponential refractory kernel into threshold (theta kernel)

      p[k].J_theta=dvector(p[k].N_theta);
      p[k].tau_theta=dvector(p[k].N_theta);
      for (l=0;l<N_theta[k];l++)
	{
	  p[k].J_theta[l]=J_theta[k][l];
	  p[k].tau_theta[l]=tau_theta[k][l];
	  //	  printf("k=%d l=%d J=%g tau=%g\n",k,l,J_theta[k][l],tau_theta[k][l]);
	}
      if (J_ref[k]!=0.)
	{
	  p[k].J_theta[N_theta[k]] = J_ref[k]*taum[k];
	  p[k].tau_theta[N_theta[k]] = taum[k];
	}

      //effective synaptic weights
      p[k].w1=dvector(Npop);
      p[k].w2=dvector(Npop);
      for (l=0;l<Npop;l++)
	{
	  p[k].w1[l] = J[k][l] * p_conn[k][l] * N[l] * taum[k] * a1[k][l] / (a1[k][l]+a2[k][l]);
	  p[k].w2[l] = J[k][l] * p_conn[k][l] * N[l] * taum[k] * a2[k][l] / (a1[k][l]+a2[k][l]);
	  //printf("k=%d l=%d a1=%g w1=%g\n",k,l,a1[k][l],p[k].w1[l]);
	}
    }

}

void clean_population(struct PopParametersGLM p[], int Npop)
{
  int k;
  for (k=0;k<Npop;k++) free_dvector(p[k].w1);
  for (k=0;k<Npop;k++) free_dvector(p[k].w2);
  for (k=0;k<Npop;k++) free_dvector(p[k].J_theta);
  for (k=0;k<Npop;k++) free_dvector(p[k].tau_theta);
}

void get_psd_with_fullparameterlist(double **SA, int Nbin, int Ntrials, int Npop, double *tref, double taum[], double taus1[][Npop], double taus2[][Npop], double taur1[][Npop], double taur2[][Npop], double a1[][Npop], double a2[][Npop], double mu[], double c[], double deltaV[], double delay[], double vth[], double vreset[], int N[], double J[][Npop], double p_conn[][Npop], int N_theta[], double Jref[], double *J_theta[], double *tau_theta[], double sigma[], double dt,double dtbin, int mode)
{
  double **signal=dmatrix(Npop,1);
  int i;
  for (i=0;i<Npop;i++)
    signal[i]=NULL;

  DT=dt;
  DTBIN=dtbin;

  struct PopParametersGLM p[Npop];
  init_population(p, Npop, tref, taum, taus1, taus2, taur1, taur2, a1, a2, mu, c, deltaV, delay, vth, vreset, N, J, p_conn, signal, N_theta, Jref, J_theta, tau_theta, sigma);

  get_psd_pop(SA, Nbin, Ntrials, p, Npop, mode);
  free_dmatrix(signal);
  clean_population(p,Npop);

}

 void get_psd_with_2D_arrays(int Nf, double SA[][Nf], int Ntrials, int Npop, double *tref, double taum[], double taus1[][Npop], double taus2[][Npop], double taur1[][Npop], double taur2[][Npop], double a1[][Npop], double a2[][Npop], double mu[], double c[], double deltaV[], double delay[], double vth[], double vreset[], int N[], double J[][Npop], double p_conn[][Npop], int N_theta[], double Jref[], double J_theta[], double tau_theta[], double sigma[], double dt,double dtbin, int mode)
{
  //convert 2D array to double**
  double *SA_tmp[Npop], *J_theta_ptr[Npop], *tau_theta_ptr[Npop];
  int i,j;
  int Nbin=2*Nf;

  int indx=0;
  for (i=0;i<Npop;i++)
    {
      J_theta_ptr[i]=&(J_theta[indx]);
      tau_theta_ptr[i]=&(tau_theta[indx]);
      indx+=N_theta[i];
      /* if (N_theta[i]>0) */
      /* 	{ */
      /* 	  J_theta_ptr[i]=&(J_theta[indx]); */
      /* 	  tau_theta_ptr[i]=&(tau_theta[indx]); */
      /* 	  indx+=N_theta[i]; */
      /* 	} */
      /* else */
      /* 	{ */
      /* 	  J_theta_ptr[i]=NULL; */
      /* 	  tau_theta_ptr[i]=NULL; */
      /* 	} */
    }

  get_psd_with_fullparameterlist(SA_tmp, Nbin, Ntrials, Npop, tref, taum, taus1, taus2, taur1, taur2, a1, a2, mu, c, deltaV, delay, vth, vreset, N, J, p_conn, N_theta, Jref, J_theta_ptr, tau_theta_ptr, sigma, dt, dtbin, mode);

  for (j=0;j<Npop;j++)
    for (i=0;i<Nf;i++)
      SA[j][i]=SA_tmp[j][i];

  for (j=0;j<Npop;j++)
    free(SA_tmp[j]);

}

void get_trajectory_with_fullparameterlist(double **A, double **rate, int Nbin, int Npop, double *tref, double taum[], double taus1[][Npop], double taus2[][Npop], double taur1[][Npop], double taur2[][Npop], double a1[][Npop], double a2[][Npop], double mu[], double c[], double deltaV[], double delay[], double vth[], double vreset[], int N[], double J[][Npop], double p_conn[][Npop], double **signal, int N_theta[], double Jref[], double *J_theta[], double *tau_theta[], double sigma[], double dt,double dtbin, int mode, int seed)
{
  DT=dt;
  DTBIN=dtbin;

  struct PopParametersGLM p[Npop];
  init_population(p, Npop, tref, taum, taus1, taus2, taur1, taur2, a1, a2, mu, c, deltaV, delay, vth, vreset, N, J, p_conn, signal, N_theta, Jref, J_theta, tau_theta, sigma);
  //  print_pop_parameters(p,Npop);

  get_trajectory(A, rate, Nbin, Npop, p, dt, dtbin, mode, seed);

  clean_population(p,Npop);
}

void get_trajectory_with_2D_arrays(int Nbin, double A[][Nbin], double rate[][Nbin], int Npop, double *tref, double taum[], double taus1[][Npop], double taus2[][Npop], double taur1[][Npop], double taur2[][Npop], double a1[][Npop], double a2[][Npop], double mu[], double c[], double deltaV[], double delay[], double vth[], double vreset[], int N[], double J[][Npop], double p_conn[][Npop], double s[][Nbin], int N_theta[], double Jref[], double J_theta[], double tau_theta[], double sigma[], double dt,double dtbin, int mode, int seed)
{

  //convert 2D array to double**
  double *rate_ptr[Npop],*AA[Npop],*signal[Npop], *J_theta_ptr[Npop], *tau_theta_ptr[Npop];
  int i;

  for (i=0;i<Npop;i++)
    AA[i]=A[i];
  for (i=0;i<Npop;i++)
    rate_ptr[i]=rate[i];
  for (i=0;i<Npop;i++)
    signal[i]=s[i];

  int indx=0;
  for (i=0;i<Npop;i++)
    {
      if (N_theta[i]>0)
  	{
  	  J_theta_ptr[i]=&(J_theta[indx]);
  	  tau_theta_ptr[i]=&(tau_theta[indx]);
  	  indx+=N_theta[i];
  	}
      else
  	{
  	  J_theta_ptr[i]=NULL;
  	  tau_theta_ptr[i]=NULL;
  	}
    }

  get_trajectory_with_fullparameterlist(AA, rate_ptr, Nbin, Npop, tref, taum, taus1, taus2, taur1, taur2, a1, a2, mu, c, deltaV, delay, vth, vreset, N, J, p_conn, signal, N_theta, Jref, J_theta_ptr, tau_theta_ptr, sigma, dt, dtbin, mode, seed);

}

//===================================================================
// SRM Model
//===================================================================

void print_pop_parameters(struct PopParametersGLM p[], int Npop)
{

  int k,l;
  printf("\n");
  for (k=0;k<Npop;k++)
    {
      printf("POPULATION %d\n",k+1);
      printf("tref=%g\n",p[k].tref);
      printf("taum=%g\n",p[k].taum);
      printf("mu=%g\n",p[k].mu);
      printf("c=%g\n",p[k].c);
      printf("DeltaV=%g\n",p[k].deltaV);
      printf("delay=%g\n",p[k].delay);
      printf("vth=%g\n",p[k].vth);
      printf("vreset=%g\n",p[k].u_reset);
      printf("N=%d\n",p[k].N);
      printf("Ntheta=%d\n",p[k].N_theta);
      for (l=0;l<p[k].N_theta;l++)
	{
	  printf("Jtheta=%g\n",p[k].J_theta[l]);
	  printf("tautheta=%g\n",p[k].tau_theta[l]);
	}
      for (l=0;l<Npop;l++)
	{
	  printf("pop%d to pop%d:\n",l+1,k+1);
	  printf("   w1=%g\n",p[k].w1[l]);
	  printf("   w2=%g\n",p[k].w2[l]);
	  printf("   taus1=%g taur1=%g\n",p[k].taus1[l],p[k].taur1[l]);
	  printf("   taus2=%g taur2=%g\n",p[k].taus2[l],p[k].taur2[l]);
	}
      printf("\n");
    }
  printf("DT=%g DTBIN=%g\n",DT,DTBIN);
}

double get_mean(double *x,int n)
{
  int i;
  double m=0;
  if (n>200) //discard initial transient
    {
      for (i=200;i<n;i++) m+=x[i];
      return m/(n-200);
    }
  else
    {
      for (i=0;i<n;i++) m+=x[i];
      return m/n;
    }
}

void print_matrix1(int N,int M, double A[][M])
{
  int i,j;
  for (i=0;i<N;i++)
    {
      for (j=0;j<M;j++)
	printf("%g ",A[i][j]);
      printf("\n");
    }
}

void print_matrix2(double **A,int N,int M)
{
  int i,j;
  for (i=0;i<N;i++)
    {
      for (j=0;j<M;j++)
	printf("%g ",A[i][j]);
      printf("\n");
    }
}

void print_matrix3(int N,int M, double *A[N])
{
  int i,j;
  for (i=0;i<N;i++)
    {
      for (j=0;j<M;j++)
	printf("%g ",A[i][j]);
      printf("\n");
    }
}

/* unsigned long int random_seed() */
/* { */

/*  unsigned int seed; */
/*  FILE *devrandom; */

/*  if ((devrandom = fopen("/dev/random","r")) == NULL) { */
/*    fprintf(stderr,"Cannot open /dev/random, setting seed to 0\n"); */
/*    seed = 0; */
/*  } else { */
/*    fread(&seed,sizeof(seed),1,devrandom); */
/*    if(verbose == D_SEED) printf("Got seed %u from /dev/random\n",seed); */
/*    fclose(devrandom); */
/*  } */

/*  return(seed); */

/* } */

#include <sys/time.h>
unsigned long int random_seed()
{

 unsigned int seed;
 struct timeval tv;
 FILE *devrandom;

 if ((devrandom = fopen("/dev/urandom","r")) == NULL) {
   gettimeofday(&tv,0);
   seed = tv.tv_sec + tv.tv_usec;
 } else {
   fread(&seed,sizeof(seed),1,devrandom);
   fclose(devrandom);
 }

 return(seed);

}

max_stars_repo_name: dekamps/miind

max_issues_repo_name: dekamps/miind

max_forks_repo_name: dekamps/miind


text:
#pragma once

#include "mexObjectiveFunction.h"

#include <nlopt.h>
#include <mex.h>

struct mexConstraintFunction
{
  static double con_fun(unsigned n, const double *x, double *gradient, void *d_);
  static void vcon_fun(unsigned m, double *result, unsigned n, const double* x, double* grad, void* f_data);

  bool isvector;
  nlopt_func fun;   // (unsigned n, const double *x, double *fc, void *func_data)
  nlopt_mfunc mfun; // (unsigned m, double *result, unsigned n, const double *x, double *fc, void *func_data)

  mxArray *prhs[2];          // feval mexMatlabCall input arguments for objective function evaluation

  nlopt_opt &opt;
  mxArray *&lasterror;        // mexObjectiveFunction's to store trapped MException
  bool &stop;                 // mexObjectiveFunction's to flag if user issued a stop

  mexConstraintFunction(mxArray * mxFun, mexObjectiveFunction &data);
  ~mexConstraintFunction();
private:
  double evalFun(unsigned n, const double *x, double *gradient);
  void evalVecFun(unsigned m, unsigned n, const double* x, double *c, double* Cgrad);
  bool call_matlab_feval_with_trap(int nlhs, mxArray *plhs[], const int n, const double *x);
};

ext: h

max_stars_repo_name: hokiedsp/matlab-nlopt

max_issues_repo_name: hokiedsp/matlab-nlopt

max_forks_repo_name: hokiedsp/matlab-nlopt


text:
#ifndef __INTERPOLATOR_H__
#define __INTERPOLATOR_H__

#include <vector>
#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>

namespace RHEHpt{

class Interpolator
{
	public:
		Interpolator (unsigned choice = 0);
		virtual ~Interpolator ();

		double operator()(unsigned order, double xp) const;

		void set_choice(unsigned choice){ _choice = choice; }

	private:
		unsigned _choice;
		gsl_interp_accel* _gsl_acc;

		/* matrix of gsl interpolators one row per _choice and one column per order */
		std::vector < std::vector < gsl_spline *> > _interpolators;

};

} /* close RHEHpt namespace* */

#endif /* __INTERPOLATOR_H__ */

max_stars_repo_name: gvita/RHEHpt

max_issues_repo_name: gvita/RHEHpt

max_forks_repo_name: gvita/RHEHpt


text:
#include <petsc.h>
#include <petscmath.h>
#include "compressibleFlow.h"
#include "mesh.h"
#include "petscdmplex.h"
#include "petscts.h"

typedef struct {
    PetscInt dim;
    PetscReal rc;
    PetscReal xc;
    PetscReal yc;
    PetscReal L;
    PetscReal Tinf;
    PetscReal rhoInf;
    PetscReal Mac;
    PetscReal MaxInf;
    PetscReal MayInf;
    PetscReal Rgas;
    PetscReal gamma;
} Constants;

typedef struct {
    Constants constants;
    FlowData flowData;
} ProblemSetup;

static PetscLogStage monitorStage;

static PetscErrorCode EulerExact(PetscInt dim, PetscReal time, const PetscReal xyz[], PetscInt Nf, PetscScalar *node, void *ctx) {
    PetscFunctionBeginUser;

    Constants *constants = (Constants *)ctx;

    // define the speed of sound and other variables in the infinity
    PetscReal aInf = PetscSqrtReal(constants->gamma*constants->Rgas*constants->Tinf);
    PetscReal pInf = constants->rhoInf*constants->Rgas*constants->Tinf;
    PetscReal uInf = constants->MaxInf*aInf;
    PetscReal vInf = constants->MayInf*aInf;

    // compute the values for the vortex
    PetscReal Tvort = constants->Tinf/(1.0 + 0.5*(constants->gamma - 1.0)*constants->Mac*constants->Mac);
    PetscReal aVort = PetscSqrtReal(constants->gamma*constants->Rgas*Tvort);

    // compute the values at each location based upon the time and the var field velocity
    PetscReal distTraveledX = time*uInf;
    PetscReal distTraveledY = time*vInf;
    // compute the updated center
    PetscReal xc = constants->xc + distTraveledX;
    PetscReal yc = constants->yc + distTraveledY;

    while(xc > constants->L*2){
        xc -= constants->L*2;
    }
    while(yc > constants->L){
        yc -= constants->L;
    }

    PetscReal x = xyz[0];
    PetscReal y = xyz[1];

    PetscReal xStar = (x- xc)/constants->rc;
    PetscReal yStar = (y - yc)/constants->rc;
    PetscReal rStar = PetscSqrtReal(PetscSqr((x - xc))+ PetscSqr((y - yc)))/constants->rc;

    PetscReal u = uInf - constants->Mac*aVort*yStar* PetscExpReal(0.5*(1.0 - PetscSqr(rStar)));
    PetscReal v = vInf + constants->Mac*aVort*xStar* PetscExpReal(0.5*(1.0 - PetscSqr(rStar)));
    PetscReal machTerm = 0.5*(constants->gamma-1.0)*PetscSqr(constants->Mac);
    PetscReal T = constants->Tinf*(1.0 - machTerm/(1+machTerm)*PetscSqr(rStar)*PetscExpReal(1-PetscSqr(rStar)) );
    PetscReal p = pInf* PetscPowReal(T/constants->Tinf, constants->gamma/(constants->gamma - 1.0));
    PetscReal rho = p/constants->Rgas/T;
    PetscReal e = p/((constants->gamma - 1.0)*rho);
    PetscReal eT = e + 0.5*(u*u + v*v);

    // Store the values
    node[RHO] = rho;
    node[RHOE] = rho*eT;
    node[RHOU + 0] = rho*u;
    node[RHOU + 1] = rho*v;

    PetscFunctionReturn(0);
}

static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx) {
    PetscFunctionBeginUser;
    PetscErrorCode     ierr;
    PetscLogStagePush(monitorStage);

    PetscInt interval = 1;
    if(step % interval == 0) {
        // Get the DM
        DM dm;
        ierr = TSGetDM(ts, &dm);
        CHKERRQ(ierr);

        ierr = VecViewFromOptions(u, NULL, "-sol_view");
        CHKERRQ(ierr);

        // Open a vtk viewer
        //    PetscViewer viewer;
        //    char        filename[PETSC_MAX_PATH_LEN];
        //    ierr = PetscSNPrintf(filename,sizeof(filename),"/Users/mcgurn/chrestScratch/results/vortex/flow%.4D.vtu",step);CHKERRQ(ierr);
        //    ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)dm),filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
        //    ierr = VecView(u,viewer);CHKERRQ(ierr);
        //    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);

        ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "TS at %f\n", time);
        CHKERRQ(ierr);

        // Compute the error
        void *exactCtxs[1];
        PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
        PetscDS ds;
        ierr = DMGetDS(dm, &ds);
        CHKERRQ(ierr);

        // Get the exact solution
        ierr = PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &exactCtxs[0]);
        CHKERRQ(ierr);

        // Create an vector to hold the exact solution
        Vec exactVec;
        ierr = VecDuplicate(u, &exactVec);
        CHKERRQ(ierr);
        ierr = DMProjectFunction(dm, time, exactFuncs, exactCtxs, INSERT_ALL_VALUES, exactVec);
        CHKERRQ(ierr);

        ierr = PetscObjectSetName((PetscObject)exactVec, "exact");
        CHKERRQ(ierr);
        ierr = VecViewFromOptions(exactVec, NULL, "-exact_view");
        CHKERRQ(ierr);

        // For each component, compute the l2 norms
        ierr = VecAXPY(exactVec, -1.0, u);
        CHKERRQ(ierr);

        PetscReal ferrors[4];
        ierr = VecSetBlockSize(exactVec, 4);
        CHKERRQ(ierr);

        ierr = PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t\n", (int)step, (double)time);
        CHKERRQ(ierr);
        ierr = VecStrideNormAll(exactVec, NORM_2, ferrors);
        CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_WORLD, "\tL_2 Error: [%2.3g, %2.3g, %2.3g, %2.3g]\n", (double)(ferrors[0]), (double)(ferrors[1]), (double)(ferrors[2]), (double)(ferrors[3]));
        CHKERRQ(ierr);

        // And the infinity error
        ierr = VecStrideNormAll(exactVec, NORM_INFINITY, ferrors);
        CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_WORLD, "\tL_Inf Error: [%2.3g, %2.3g, %2.3g, %2.3g]\n", (double)ferrors[0], (double)ferrors[1], (double)ferrors[2], (double)ferrors[3]);
        CHKERRQ(ierr);

        ierr = PetscObjectSetName((PetscObject)exactVec, "error");
        CHKERRQ(ierr);
        ierr = VecViewFromOptions(exactVec, NULL, "-exact_view");
        CHKERRQ(ierr);
        ierr = VecDestroy(&exactVec);
        CHKERRQ(ierr);
    }
    PetscLogStagePop();
    PetscFunctionReturn(0);
}

static PetscErrorCode PhysicsBoundary_Euler(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx) {
    PetscFunctionBeginUser;
    Constants *constants = (Constants *)ctx;

    // Offset the calc assuming the cells are square
    PetscReal x[3];

    for(PetscInt i =0; i < constants->dim; i++){
        x[i] = c[i] + n[i]*0.5;
    }

    EulerExact(constants->dim, time, x, 0, a_xG, ctx);
    PetscFunctionReturn(0);
}

int main(int argc, char **argv)
{
    // Setup the problem
    Constants constants;

    // sub sonic
    constants.dim = 2;
    constants.L = 1.0;
    constants.xc = constants.L/2.0;
    constants.yc = constants.L/2.0;
    constants.Tinf =298.0;
    constants.rhoInf = 1.0;
    constants.Rgas = 287.0;
    constants.gamma = 1.4;
    constants.rc = .1*constants.L;
    constants.gamma = 1.4;
    constants.Mac = 0.3;
    constants.MaxInf = 0.3;
    constants.MayInf = 0.0;

    PetscErrorCode ierr;
    // create the mesh
    // setup the ts
    DM dm;                 /* problem definition */
    TS ts;                 /* timestepper */

    // initialize petsc and mpi
    PetscInitialize(&argc, &argv, NULL, "HELP");

    // Create a ts
    ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);
    ierr = TSSetProblemType(ts, TS_NONLINEAR);CHKERRQ(ierr);
    ierr = TSSetType(ts, TSRK);CHKERRQ(ierr);
    ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);

    PetscInt lengthFactor = 2;
    ierr =  PetscOptionsGetInt(NULL, NULL, "-lengthFactor", &lengthFactor, NULL);CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD, "LengthFactor %d\n", lengthFactor);

    // Create a mesh
    // hard code the problem setup
    PetscReal start[] = {0.0, 0.0};
    PetscReal end[] = {constants.L*lengthFactor, constants.L};
    PetscReal nxHeight = 10;
    PetscInt nx[] = {lengthFactor*nxHeight, nxHeight};
    ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, constants.dim, PETSC_FALSE, nx, start, end, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr);

    // Setup the flow data
    FlowData flowData;     /* store some of the flow data*/
    ierr = FlowCreate(&flowData);CHKERRQ(ierr);

    // Combine the flow data
    ProblemSetup problemSetup;
    problemSetup.flowData = flowData;
    problemSetup.constants = constants;

    //Setup
    CompressibleFlow_SetupDiscretization(flowData, &dm);

    // Add in the flow parameters
    PetscScalar params[TOTAL_COMPRESSIBLE_FLOW_PARAMETERS];
    params[CFL] = 0.5;
    params[GAMMA] = constants.gamma;

    // set up the finite volume fluxes
    CompressibleFlow_StartProblemSetup(flowData, TOTAL_COMPRESSIBLE_FLOW_PARAMETERS, params);
    DMView(flowData->dm, PETSC_VIEWER_STDOUT_WORLD);
    // Add in any boundary conditions
    PetscDS prob;
    ierr = DMGetDS(flowData->dm, &prob);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    const PetscInt idsLeft[]= {1, 2, 3, 4};
    ierr = PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall left", "Face Sets", 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler, NULL, 4, idsLeft, &constants);CHKERRQ(ierr);

    // Complete the problem setup
    ierr = CompressibleFlow_CompleteProblemSetup(flowData, ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Name the flow field
    ierr = PetscObjectSetName(((PetscObject)flowData->flowField), "Numerical Solution");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Setup the TS
    ierr = TSSetFromOptions(ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
//    ierr = TSMonitorSet(ts, MonitorError, &constants, NULL);CHKERRQ(ierr);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // set the initial conditions
    PetscErrorCode     (*func[2]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {EulerExact};
    void* ctxs[1] ={&constants};
    ierr    = DMProjectFunction(flowData->dm,0.0,func,ctxs,INSERT_ALL_VALUES,flowData->flowField);CHKERRQ(ierr);

    // for the mms, add the exact solution
    ierr = PetscDSSetExactSolution(prob, 0, EulerExact, &constants);CHKERRQ(ierr);

    // Output the mesh
    ierr = DMViewFromOptions(flowData->dm, NULL, "-dm_view");CHKERRQ(ierr);

    // Compute the end time so it goes around once
    PetscReal aInf = PetscSqrtReal(constants.gamma*constants.Rgas*constants.Tinf);
    PetscReal u_x = constants.MaxInf*aInf;
    PetscReal endTime = constants.L/u_x;

    TSSetMaxTime(ts, endTime);
    TSSetMaxSteps(ts, 50);
    PetscDSView(prob, PETSC_VIEWER_STDOUT_WORLD);
    PetscLogStage solveStage;
    PetscLogStageRegister("TSSolve",&solveStage);
    PetscLogStageRegister("Monitor", &monitorStage);

    PetscLogStagePush(solveStage);
    ierr = TSSolve(ts,flowData->flowField);CHKERRQ(ierr);
    PetscLogStagePop();
    PetscPrintf(PETSC_COMM_WORLD, "Done");
    FlowDestroy(&flowData);
    TSDestroy(&ts);

    return PetscFinalize();

}

ext: c

max_stars_repo_name: mmcgurn/MattFlowCases

max_issues_repo_name: mmcgurn/MattFlowCases

max_forks_repo_name: mmcgurn/MattFlowCases


text:
#include "imd.h"
#include <sys/time.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>

// #define USEFLOAT  // hauptsächlich in der funktion genexptint. Profiling zeigte, dass
                  // hier die meiste zeit verbraucht wird -> float verdoppelt performance

#ifdef USEFLOAT
typedef float Real;
#define REALTYPE MPI_FLOAT
#else
typedef double Real;
#define REALTYPE MPI_DOUBLE
#endif

#ifdef USEFLOAT
  #define EXPR     expf //exp zu floaten ist eine ganz mieeese idee
  #define SQRTR    sqrtf
  #define POWR     powf
  #define LOGR     logf
#else
  #define EXPR     exp
  #define SQRTR    sqrt
  #define POWR     pow
#define LOGR       log
#endif

// *********************************************************
// PHYSICAL CONSTANTS
// *********************************************************
// const double eV2J=1.6021766E-19;
const Real eV2H=0.03674932; //eV to Hartree
const Real colrad_reltol=1e-5;
const Real colrad_abstol=10.0;

// const Real J2eV=6.2415091E18;
const Real planck=6.62607004E-34; // J/s
const Real bohr_radius=0.52917721067E-10; // m
const Real bohr_radius_sq=2.800285202924816e-21;
const Real hbar_cub=1.172812163789953e-102; //hbar^3
const Real double_emass_pow_3_2 = 2.459112949719466e-45; // (2*emass)^3/2
const int MAXLINE = 255;
const Real  pi=3.141592653589793;
const Real  pi_sq=9.869604401089358;
const Real E_ion_H=13.6; // eV
const Real E_ion_H_J=2.178960176000000e-18; // J
const Real E_ion_H_sq_J=4.747867448593952e-36;

const Real colrad_tequi=1e-12;//TEST// 1e-12; //bei initial equi ohne Temperatur-variation erst einmal
                                //die Saha-besetzungsdichten equilibrieren

//const double  LIGHTSPEED=2.997925458e8; // m/s
Real  LASERFREQ;

int colrad_ydot(double t, N_Vector u, N_Vector udot, void *user_data);
void do_Saha(Real Te,Real totalc,Real ne,N_Vector y);
int colrad_GetCoeffs(N_Vector y,Real It, void * user_data);

// Die Zwei müssen nach Prototypes.h
// void do_colrad(double dt);
// void colrad_init(void);

void colrad_read_states(void);
void colrad_Saha_init(int i,int j,int k);

// ******************************************************************************
// *                CROSS SECTION INTEGRATION STUFF
// ******************************************************************************

gsl_integration_workspace * winteg_inner=NULL;
gsl_integration_workspace * winteg_outer=NULL;
gsl_integration_workspace * winteg_fermi=NULL;
gsl_integration_workspace * winteg_exc=NULL; //excitation

gsl_integration_romberg_workspace * winteg_rb_inner=NULL;
gsl_integration_romberg_workspace * winteg_rb_outer=NULL;

struct my_f_params { Real ne; Real T;Real mu; Real E;Real DeltaE; int allowed;};
// struct my_f_params fparams_inner; //For inner integrand
// struct my_f_params fparams_outer; //outer integrand
// struct my_f_params fparams_fermi;
// struct my_f_params fparams_exc;

double inner_integrand_ionization(double x, void *p); // integrate along E'
double outer_integrand_ionization(double x,void *p);  // integrate along E

Real double_integral_ionization(Real ne,Real T, Real mu, Real DeltaE); //evaluates double integral

double inner_integrand_recombination(double x, void *p);
double outer_integrand_recombination(double x,void *p);
Real double_integral_recombination(Real ne,Real T, Real mu, Real DeltaE);

double integrand_excitation(double x,void *p);
Real eval_excitation_integral(Real ne,Real T,Real mu, Real DeltaE, int allowed);

Real eval_dexcitation_integral(Real ne,Real T,Real mu, Real DeltaE, int allowed);
double integrand_deexcitation(double x,void *p);

double fermi_integrand(double x, void *p);
Real eval_fermi_integrand(Real ne,Real T, Real mu);

double integrand_excitation_debug(double x,void *p);

double outer_integrand_ionization2(double x,struct my_f_params* p);
Real double_integral_ionization2(Real ne,Real T, Real mu, Real DeltaE); //evaluates double integral
double inner_integrand_ionization2(double x, struct my_f_params* p);

// **********************************************************************************************
// *          PAR INTEGRAL STUFF
// **********************************************************************************************
int terminate_gkq;
int terminate_gkq_outer;
int terminate_gkq_inner;
int terminate_serial;
int gkq_iter_serial; // nr of iterations

const double gkq_alpha=0.816496580927726;
const double gkq_beta=0.447213595499958;
static const double xgkq[12] =
{
  0.0,
  -0.942882415695480,
  -0.816496580927726,
  -0.641853342345781,
  -0.447213595499958,
  -0.236383199662150,
  0.0,
  0.236383199662150,
  0.447213595499958,
  0.641853342345781,
  0.816496580927726,
  0.942882415695480
};

Real integral_simpson(Real (*f)(Real, void*), Real a, Real b,int n,void* p);

int simpson_error;
const Real tolmax=1e-20;
const Real simpson_itermax=120;
#define INITIAL_STACK_SIZE 128   /* initial size of new stacks */

/* the stack structure */
struct stack_s{
  int el_count;            /* count of elements on stack */
  int el_size;             /* size of an element */
  int mem_reserve;         /* allocated memory for stack */
  void* elements;          /* pointer to begin of stack */
};

typedef struct _work_t{
  double a;
  double b;
  double tol;
  double S;
  double fa;
  double fb;
  double fm;
  double rec;
  int iter;
  struct my_f_params * p; //pointer auf params
} work_t;

typedef struct _work_t_gkq{
  double a;
  double b;
  double toler;
  double I_13;
  double I_prev;
  double fa;
  double fb;
  struct my_f_params * p; //pointer auf params
  shortint iter;
} work_gkq;

typedef struct stack_s* stack_t;
double integral_simpson_par(double (*f)(double, struct my_f_params*), stack_t stack);

double gkq_adapt_OMP(double (*f)(double, struct my_f_params*), stack_t stack);
double gkq_OMP(double (*f)(double, struct my_f_params*), double a, double b, double TOL, struct my_f_params* p,stack_t stack);

double gkq_serial(double (*f)(double, struct my_f_params*), double a, double b, double TOL, struct my_f_params* p);
double gkq_adapt_serial(double (*f)(double, struct my_f_params*), double a, double b,
                        double fa,double fb, double toler,double I_13, struct my_f_params* p);

// void create_stack(stack_t* stack, int element_size);
// int  empty_stack(stack_t stack);
// void push_stack(stack_t stack, void* element);
// void pop_stack(stack_t stack, void* element);

/******************************************
 * create new stack
 ******************************************/
void create_stack(
             stack_t* stack,     /* stack to create */
             int element_size)   /* size of a stack element */
{
  int initial_size = INITIAL_STACK_SIZE;

  /* allocate memory for new stack struct */
  (*stack) = (stack_t) malloc(sizeof(struct stack_s));
  if (!(*stack)){
    char errstr[255];
    sprintf(errstr, "error: could not allocate memory for stack.. Abort.\n");
    error(errstr);
    // exit(1);
  }

  /* allocate memory for stack elements */
  (*stack)->elements = (void*) malloc(element_size * initial_size);
  (*stack)->mem_reserve = initial_size;
  if (!(*stack)->elements){
    char errstr[255];
    sprintf(errstr, "error: could not allocate memory for stack.. Abort.\n");
    error(errstr);
  }

  (*stack)->el_size = element_size;
  (*stack)->el_count = 0;

}

/*****************************************
 * check if the stack is empty
 *****************************************/
int empty_stack(stack_t stack)
{
  return stack->el_count <= 0;
}

/*****************************************
 * push a element on stack
 *****************************************/
void push_stack(stack_t stack,    /* target stack */
           void* element)    /* element to push */
{
  int i, new_reserve;
  int log2_count;

  /* check if we need more memory for stack */
  if (stack->el_count >= stack->mem_reserve)
  {

      /* calculate new size for the stack
         it should be a power of two */
      for (i = stack->el_count, log2_count = 0;
           i > 0;
           i>>1, log2_count++);
      new_reserve = 1 << log2_count;

      /* reallocate memory for phase thread tables
         and nullify new values */
      stack->elements = (void *) realloc(stack->elements,
                      stack->el_size * new_reserve);
      if (!stack->elements){
        char errstr [255];
        sprintf(errstr, "error: can't reallocate stack.. Aborting\n");
        error(errstr);
        // exit(1);
      }

      stack->mem_reserve = new_reserve;
  }

  /* now push the element on top of the stack */
  memcpy((char*)stack->elements + stack->el_count*stack->el_size,
            element, stack->el_size);
  stack->el_count++;

}

/*****************************************
 * pop an element from stack
 *****************************************/
void pop_stack(
          stack_t stack,    /* target stack */
          void* element)    /* where poped el. should be stored */
{
  if (stack->el_count <= 0){
    char errstr[255];
    sprintf(errstr, "error: trying to pop from empty stack.\n");
    error(errstr);
    // exit(2);
  }

  stack->el_count--;
  memcpy(element,
          (char*)stack->elements + stack->el_count*stack->el_size,
          stack->el_size);

}

// ***************************************************************************
// *      Gauss-kronard quadrature, parallel
// ***************************************************************************
double gkq_OMP(double (*f)(double, struct my_f_params*), double a, double b, double TOL, struct my_f_params* p,stack_t stack)
{
  //1st integration
  double result=0.0;
// *********************************************

  double m=0.5*(a+b);
  double h=0.5*(b-a);
  double y[13];
  double fa=y[0]=f(a,p);
  double fb=y[12]=f(b,p);
  int i;
  for(i=1;i<12;i++)
    y[i]=f(m+xgkq[i]*h,p);

  double I_4= (h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));  // 4-point gauss-lobatto
  double I_7= (h/1470.0)*(77.0*(y[0]+y[12])+432.0*(y[2]+y[10])+  // 7-point kronrod
               625.0*(y[4]+y[8])+672.0*y[6]);

  double I_13= h*(0.0158271919734802*(y[0]+y[12])+0.0942738402188500*(y[1]+y[11])+0.155071987336585*(y[2]+y[10])+
                  0.188821573960182*(y[3]+y[9])+0.199773405226859*(y[4]+y[8])+0.224926465333340*(y[5]+y[7])+
                  0.242611071901408*y[6]);   //13-point Kronrod

  double Err1=fabs(I_7-I_13);
  double Err2=fabs(I_4-I_13);

  double r=(Err2 != 0.0) ? Err1/Err2 : 1.0;
  double toler=(r > 0.0 && r < 1.0) ? TOL/r : TOL;

  if(I_13 == 0)
    I_13=b-a;
  I_13=fabs(I_13);

  //Prepare work and push onto stack
  work_gkq work;

  work.a = a;
  work.b = b;
  work.toler = toler;
  work.I_13=I_13;
  work.fa=fa;
  work.fb=fb;
  work.p=p;
  work.I_prev=I_7;

  //ANTI-FOLGENDES:
  //OUT OF TOLERANCE !!!, mll:3.0162e-18, a:3.0162e-18, b:3.0162e-18, mrr:3.0162e-18,I_7-I_4:0.0000e+00, tol:1.6002e-315,I_13:7.0585e-313
  if(I_13 < 1e-150)
    return 0;

  push_stack(stack, &work);
  result=gkq_adapt(f,stack);
  return result;
}

double gkq_serial(double (*f)(double, struct my_f_params*), double a, double b, double TOL, struct my_f_params* p)
{
  //1st integration

  double result=0.0;
  gkq_iter_serial=0;
// *********************************************

  double m=0.5*(a+b);
  double h=0.5*(b-a);
  double y[13];
  double fa=y[0]=f(a,p);
  double fb=y[12]=f(b,p);
  int i;
  for(i=1;i<12;i++)
    y[i]=f(m+xgkq[i]*h,p);

  double I_4= (h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));  // 4-point gauss-lobatto
  double I_7= (h/1470.0)*(77.0*(y[0]+y[12])+432.0*(y[2]+y[10])+  // 7-point kronrod
               625.0*(y[4]+y[8])+672.0*y[6]);

  double I_13= h*(0.0158271919734802*(y[0]+y[12])+0.0942738402188500*(y[1]+y[11])+0.155071987336585*(y[2]+y[10])+
                  0.188821573960182*(y[3]+y[9])+0.199773405226859*(y[4]+y[8])+0.224926465333340*(y[5]+y[7])+
                  0.242611071901408*y[6]);   //13-point Kronrod

  double Err1=fabs(I_7-I_13);
  double Err2=fabs(I_4-I_13);

  double r=(Err2 != 0.0) ? Err1/Err2 : 1.0;
  double toler=(r > 0.0 && r < 1.0) ? TOL/r : TOL;

  if(I_13 == 0)
    I_13=b-a;
  I_13=fabs(I_13);

  result=gkq_adapt_serial(f,a,b,fa,fb,toler,I_13, p);

  return result;
}

// ***********************************************
// * RECURSIVE ADAPTION ROUTINE FOR PARALLEL-GK-QUADRATURE
// **********************************************
double gkq_adapt_OMP(double (*f)(double, struct my_f_params*), stack_t stack)
{
  work_gkq work;
  work.iter=0;
  int ready, idle, busy;
  double integral_result = 0.0;

  busy = 0;
  terminate_gkq=0;

#pragma omp parallel default(none) \
    shared(stack, integral_result,f,busy,terminate_gkq,myid) \
    private(work, idle, ready)
  {
// printf("me:%d, err:%d\n",omp_get_thread_num(),simpson_error);

    ready = 0;
    idle = 1;

    while(!ready) // && !terminate_gkq)//  && !simpson_error) //<-- so NICHT!
    {
      #pragma omp critical (stack)
      {
        if (!empty_stack(stack))
        {
          /* we have new work */
          pop_stack(stack, &work);
          if (idle)
          {
            /* say others i'm busy */
            busy += 1;
            idle = 0;
          }
        }
        else
        {
          /* no new work on stack */
          if (!idle){
            busy -= 1;
            idle = 1;
          }

          /* nobody has anything to do; let us leave the loop */
          if (busy == 0)
          {
            ready = 1;
          }
        }
      } /* end critical(stack) */

      if (idle)
        continue; //if ready==1 --> leave loop
double I_prev=work.I_prev;

      double a = work.a;
      double b = work.b;
      double toler = work.toler;
      double I_13=work.I_13;
      double fa=work.fa;
      double fb=work.fb;
      int iter=work.iter;
      // double *y= work.y; // brauch ich nicht!
      struct my_f_params * p = work.p;

      double m = (a+b)/2;
      double h = (b -a)/2;
      double mll=m-gkq_alpha*h;
      double ml=m-gkq_beta*h;
      double mr=m+gkq_beta*h;
      double mrr=m+gkq_alpha*h;

      double fmll=f(mll,p);
      double fml=f(ml,p);
      double fm=f(m,p);
      double fmr=f(mr,p);
      double fmrr=f(mrr,p);
      double I_4=h/6.0*(fa+fb+5.0*(fml+fmr));   // 4-point Gauss-Lobatto formula.
      double I_7=h/1470.0*(77.0*(fa+fb)+432.0*(fmll+fmrr)+625.0*(fml+fmr)+672.0*fm);

// if(myid==1)
//   printf("I_7:%.4e, I_13:%.4e,I_4:%.4e, minus:%.4e, to:%.4e\n",I_7,I_13,I_4,I_7-I_4, toler*I_13);
int maxiter=50; //max. subdivisions
double abstol=1e-30;
work.I_prev=I_7; // für abstolcheck in nächster recursion

      if (fabs(I_7-I_4) <= toler*I_13 || mll <= a || b <= mrr || iter > maxiter || fabs(I_7-I_prev) < abstol )
      {
        if ((mll <= a || b <= mrr)) //Error
         {
           // out_of_tolerance=true; // Interval contains no more machine numbers
           // printf("OUT OF TOLERANCE !!!, mll:%.4e, a:%.4e, b:%.4e, mrr:%.4e,I_7-I_4:%.4e, tol:%.4e,I_13:%.4e\n",
           //        mll,b,b,mrr,I_7-I_4, toler*I_13,I_13);
           terminate_gkq=1;
         }
        #pragma omp critical (integral_result)
        {
          integral_result += I_7;      //Terminate recursion.
        }
// printf("me ok:%d, a:%f,b:%f, tler:%.5e,I_4:%f,I_7:%f,ubteg;%.4e\n", omp_get_thread_num(), a,b,toler,I_4,I_7,integral_result);

      }
      else  //subdivide interval and push new work on stack
      {
        #pragma omp critical (stack)
        {

          // printf("me NOOOO:%d, a:%f,b:%f, tler:%.5e,I_4:%f,I_7:%f\n", omp_get_thread_num(), a,b,toler,I_4,I_7);

          work.iter=iter+1;

          work.a=a;
          work.b=mll;
          work.fa=fa;
          work.fb=fmll;
          push_stack(stack, &work);

          work.a=mll;
          work.b=ml;
          work.fa=fmll;
          work.fb=fml;
          push_stack(stack, &work);

          work.a=ml;
          work.b=m;
          work.fa=fml;
          work.fb=fm;
          push_stack(stack, &work);

          work.a=m;
          work.b=mr;
          work.fa=fm;
          work.fb=fmr;
          push_stack(stack, &work);

          work.a=mr;
          work.b=mrr;
          work.fa=fmr;
          work.fb=fmrr;
          push_stack(stack, &work);

          work.a=mrr;
          work.b=b;
          work.fa=fmrr;
          work.fb=fb;
          push_stack(stack, &work);

        } // pragma critical stack
      }   // else ..non-acceptable error
    } // while
  } /* end omp parallel */
  return integral_result;
}

double gkq_adapt_serial(double (*f)(double, struct my_f_params*), double a, double b, double fa,
                        double fb, double toler,double I_13, struct my_f_params* p)
{
  double m = (a+b)/2;
  double h = (b -a)/2;
  double mll=m-gkq_alpha*h;
  double ml=m-gkq_beta*h;
  double mr=m+gkq_beta*h;
  double mrr=m+gkq_alpha*h;

  double fmll=f(mll,p);
  double fml=f(ml,p);
  double fm=f(m,p);
  double fmr=f(mr,p);
  double fmrr=f(mrr,p);
  double I_4=h/6.0*(fa+fb+5.0*(fml+fmr));   // 4-point Gauss-Lobatto formula.
  double I_7=h/1470.0*(77.0*(fa+fb)+432.0*(fmll+fmrr)+625.0*(fml+fmr)+672.0*fm);

  gkq_iter_serial++;

  if ( (fabs(I_7-I_4) <= toler*I_13 || mll <= a || b <= mrr) && gkq_iter_serial)
  {
    if ((mll <= a || b <= mrr) && !terminate_serial) //Error
     {
      // out_of_tolerance=true; // Interval contains no more machine numbers
      printf("OUT OF TOLERANCE !!!, mll:%.4e, a:%.4e, b:%.4e, mrr:%.4e\n", mll,b,b,mrr);
      terminate_serial=1;
     }

// printf("me ok:%d, a:%f,b:%f, tler:%.5e,I_4:%f,I_7:%f\n", omp_get_thread_num(), a,b,toler,I_4,I_7);
    return I_7;
  }
  else
  {

    // printf("me NOOOO:%d, a:%f,b:%f, tler:%.5e,I_4:%f,I_7:%f\n", omp_get_thread_num(), a,b,toler,I_4,I_7);

    return  gkq_adapt_serial(f, a,mll,fa,fmll,toler,I_13,p) +
            gkq_adapt_serial(f, mll,ml,fmll,fml,toler,I_13,p) +
            gkq_adapt_serial(f, ml,m,fml,fm,toler,I_13,p) +
            gkq_adapt_serial(f, m,mr,fm,fmr,toler,I_13,p) +
            gkq_adapt_serial(f, mr,mrr,fmr,fmrr,toler,I_13,p) +
            gkq_adapt_serial(f, mrr,b,fmrr,fb,toler,I_13,p);
  }

}

max_stars_repo_name: fmqeisfeld/IMD

max_issues_repo_name: fmqeisfeld/IMD

max_forks_repo_name: fmqeisfeld/IMD


text:
static char help[] =
"Solve a hyperbolic system in one space dimension (1D):\n"
"    q_t + F(t,x,q)_x = g(t,x,q)\n"
"where solution q(t,x), flux F(t,x,q), and source g(t,x,q) are column vectors\n"
"of length n.  The domain is (t,x) in [0,T] x [a,b].  The initial condition is\n"
"q(0,x) = f(x).  The flux may be of the form\n"
"    F(t,x,q) = A(t,x,q) q\n"
"but this is not required.  Flux boundary conditions are used in all cases,\n"
"    F = bdryflux_a(t,qright)\n"
"    F = bdryflux_b(t,qleft)\n"
"including nonreflecting, reflecting, outflow, and periodic boundary\n"
"conditions.\n\n"
"Uses a finite volume discretization so the grid values represent cell\n"
"averages.  Each case implements a Riemann solver\n"
"    F = faceflux(t,x,qleft,qright)\n"
"at cell faces.  (The Riemann solver computes the value of the solution on\n"
"the cell face going forward in time.)  Implements the following\n"
"slope-limiters when computing fluxes:\n"
"    -limiter none       Godunov's method, i.e. first-order upwinding\n"
"    -limiter fromm      formula (6.14) in LeVeque 2002\n"
"    -limiter mc         formula (6.29)\n"
"    -limiter minmod     formula (6.26)\n"
"Control the spatial grid by PETSc option\n"
"    -da_grid_x M                             [grid of M cells/points]\n\n"
"Time stepping is by semi-discretization in space (method of lines) and then\n"
"application of PETSc's (generally) adaptive and higher-order TS solvers.\n"
"Control time stepping and solution information by these PETSc options, among\n"
"others:\n"
"    -ts_monitor                              [shows time steps]\n"
"    -ts_monitor_solution draw                [generate simple movie]\n"
"        -draw_pause 0.1 -draw_size 2000,200  [control the movie]\n"
"    -ts_type                                 [default is rk]\n"
"        -ts_rk_type X                        [default is 3bs]\n"
"    -ts_dt 0.01 -ts_adapt_type none          [turn off adaptive]\n"
"Note that TS solver type SSP (-ts_type ssp) is recommended for these\n"
"hyperbolic problems.\n\n"
"Use option -problem selects the problem case:\n"
"    -problem acoustic   wave equation in system form (n=2) [default]\n"
"    -problem advection  scalar advection equation (n=1)\n"
"    -problem swater     shallow water equations (n=2)\n"
"    -problem traffic    scalar, nonlinear traffic equation (n=1)\n"
"To see possible initial conditions for problem X see the corresponding\n"
".h file.  Based on the problem-specific possiblities, set\n"
"    -initial Y\n\n"
"See the makefile for test examples, and do 'make test' to test.\n\n"
"This program is documented by the slides in the fvolume/ directory at\n"
"https://github.com/bueler/slide-teach\n\n";

#include <petsc.h>

/* The struct "ProblemCtx" is defined in cases.h.  The comments in cases.h
show how to add new problems. */
#include "cases.h"

// minmod(a,b) as define on LeVeque page 111
static PetscReal minmod(PetscReal a, PetscReal b) {
    if (a*b > 0) // both signs agree
        return (PetscAbs(a) < PetscAbs(b)) ? a : b;
    else
        return 0.0;
}

typedef enum {NONE,FROMM,MC,MINMOD} LimiterType;
static const char* LimiterTypes[] = {"none","fromm","mc","minmod",
                                     "LimiterType", "", NULL};

static LimiterType limiter = NONE;     // slope-limiter

extern PetscErrorCode FormInitial(DMDALocalInfo*, Vec, PetscReal, ProblemCtx*);
extern PetscErrorCode GetMaxSpeed(DMDALocalInfo*, Vec, PetscReal, PetscReal*, ProblemCtx*);
extern PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo*, PetscReal,
        PetscReal*, PetscReal*, void*);

int main(int argc,char **argv) {
    PetscErrorCode   ierr;
    TS               ts;                 // ODE solver for method-of-lines (MOL)
    DM               da;                 // structured grid
    Vec              q;                  // the solution
    DMDALocalInfo    info;               // structured grid info
    ProblemType      problem = ACOUSTIC; // which problem we are solving
    ProblemCtx       user;               // problem-specific information
    PetscInt         swidth, k, steps;
    PetscBool        flg;
    PetscReal        hx, qmin, qmax, t0, tf, dt, c;

    PetscInitialize(&argc,&argv,(char*)0,help);

    // get which problem we are solving
    // (ProblemType, ProblemTypes, InitializerPtrs are defined in cases.h)
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"",
               "riemann (hyperbolic system solver) options",""); CHKERRQ(ierr);
    ierr = PetscOptionsEnum("-problem", "problem type",
               "riemann.c",ProblemTypes,(PetscEnum)(problem),(PetscEnum*)&problem,
               NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnum("-limiter", "limiter type",
               "riemann.c",LimiterTypes,(PetscEnum)(limiter),(PetscEnum*)&limiter,
               NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    // call the initializer for the given case
    // (it allocates list of strings in user->field_names thus PetscFree below)
    ierr = (*InitializerPtrs[problem])(&user); CHKERRQ(ierr);

    // create grid
    swidth = (limiter == NONE) ? 1 : 2;
    ierr = DMDACreate1d(PETSC_COMM_WORLD,
                        user.periodic_bcs ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE,
                        4,          // default resolution
                        user.n_dim, // system dimension (d.o.f.)
                        swidth,     // stencil (half) width
                        NULL,&da); CHKERRQ(ierr);
    ierr = DMSetFromOptions(da); CHKERRQ(ierr);
    ierr = DMSetUp(da); CHKERRQ(ierr);
    ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    hx = (user.b_right - user.a_left) / info.mx;
    ierr = DMDASetUniformCoordinates(da,user.a_left+hx/2.0,user.b_right-hx/2.0,
                                     0.0,1.0,0.0,1.0); CHKERRQ(ierr);

    // set field names so that visualization makes sense
    for (k = 0; k < info.dof; k++) {
        ierr = DMDASetFieldName(da,k,(user.field_names)[k]); CHKERRQ(ierr);
    }

    // create TS:  dq/dt = G(t,q)  form
    ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
    ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
    ierr = TSSetDM(ts,da); CHKERRQ(ierr);
    ierr = DMDATSSetRHSFunctionLocal(da,INSERT_VALUES,
           (DMDATSRHSFunctionLocal)FormRHSFunctionLocal,&user); CHKERRQ(ierr);
    ierr = TSSetType(ts,TSRK); CHKERRQ(ierr);  // defaults to -ts_rk_type 3bs

    // set up time axis
    ierr = TSSetTime(ts,user.t0_default); CHKERRQ(ierr);
    ierr = TSSetMaxTime(ts,user.tf_default); CHKERRQ(ierr);
    dt = user.tf_default - user.t0_default;
    ierr = TSSetTimeStep(ts,dt); CHKERRQ(ierr);  // usually reset below
    ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
    ierr = TSSetFromOptions(ts); CHKERRQ(ierr);

    // get initial values
    ierr = DMCreateGlobalVector(da,&q); CHKERRQ(ierr);
    ierr = TSGetTime(ts,&t0); CHKERRQ(ierr);
    ierr = FormInitial(&info,q,t0,&user); CHKERRQ(ierr);
    //ierr = VecView(q,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

    // use CFL to reset initial time-step dt (unless user sets)
    ierr = PetscOptionsHasName(NULL,NULL,"-ts_dt",&flg); CHKERRQ(ierr);
    ierr = GetMaxSpeed(&info,q,t0,&c,&user); CHKERRQ(ierr);
    if (!flg && c > 0.0) {
        ierr = TSGetMaxTime(ts,&tf); CHKERRQ(ierr);
        dt = PetscMin(hx / c, tf-t0);
        ierr = TSSetTimeStep(ts,dt); CHKERRQ(ierr);
    } else {
        ierr = TSGetTimeStep(ts,&dt); CHKERRQ(ierr);
    }

    // solve
    ierr = PetscPrintf(PETSC_COMM_WORLD,
               "solving problem %s, a system of %d equations,\n"
               "  on %d-point grid with dx=%.6f and initial dt=%.6f...\n",
               ProblemTypes[problem],info.dof,info.mx,hx,dt); CHKERRQ(ierr);
    ierr = TSSolve(ts,q); CHKERRQ(ierr);

    // report on solution
    ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr);
    ierr = TSGetTime(ts,&tf); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,
               "  ... completed %d steps for %.4f <= t <= %.4f\n",
               steps,t0,tf); CHKERRQ(ierr);
    for (k = 0; k < info.dof; k++) {
        ierr = VecStrideMin(q,k,NULL,&qmin); CHKERRQ(ierr);
        ierr = VecStrideMax(q,k,NULL,&qmax); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_WORLD,
               "  range of %d component: %8.5f <= %s <= %8.5f\n",
               k,qmin,(user.field_names)[k],qmax); CHKERRQ(ierr);
    }

    // free memory
    VecDestroy(&q);  TSDestroy(&ts);  DMDestroy(&da);
    ierr = PetscFree(user.field_names); CHKERRQ(ierr);
    return PetscFinalize();
}

PetscErrorCode FormInitial(DMDALocalInfo *info, Vec q, PetscReal t0, ProblemCtx *user) {
    PetscErrorCode   ierr;
    const PetscReal  hx = (user->b_right - user->a_left) / info->mx;
    PetscInt         j;
    PetscReal        x, *aq;

    ierr = DMDAVecGetArray(info->da, q, &aq); CHKERRQ(ierr);
    for (j=info->xs; j<info->xs+info->xm; j++) {
        x = user->a_left + (j+0.5) * hx;
        ierr = user->f_initial(t0,x,&aq[(info->dof)*j]); CHKERRQ(ierr);
    }
    ierr = DMDAVecRestoreArray(info->da, q, &aq); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode GetMaxSpeed(DMDALocalInfo *info, Vec q, PetscReal t,
                           PetscReal *maxspeed, ProblemCtx *user) {
    PetscErrorCode   ierr;
    const PetscReal  hx = (user->b_right - user->a_left) / info->mx;
    PetscInt         j;
    PetscReal        x, cj, locmax, *aq;
    MPI_Comm         comm;

    ierr = DMDAVecGetArray(info->da, q, &aq); CHKERRQ(ierr);
    locmax = 0.0;
    for (j=info->xs; j<info->xs+info->xm; j++) {
        x = user->a_left + (j+0.5) * hx;
        ierr = user->maxspeed(t,x,&aq[(info->dof)*j],&cj); CHKERRQ(ierr);
        locmax = PetscMax(locmax,cj);
    }
    ierr = DMDAVecRestoreArray(info->da, q, &aq); CHKERRQ(ierr);
    ierr = PetscObjectGetComm((PetscObject)info->da,&comm); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&locmax,maxspeed,1,MPIU_REAL,MPIU_MAX,comm); CHKERRQ(ierr);
    return 0;
}

static inline void ncopy(PetscInt n, PetscReal *src, PetscReal *tgt) {
    PetscInt k;
    for (k = 0; k < n; k++)
        tgt[k] = src[k];
}

static inline void slopemodify(PetscInt n, PetscReal C,
                               PetscReal *sig, PetscReal *Q_in,
                               PetscReal *Q_out) {
    PetscInt k;
    for (k = 0; k < n; k++)
        Q_out[k] = Q_in[k] + C * sig[k];
}

// Right-hand-side of method-of-lines discretization form of PDE.  Implements
// Gudonov (i.e. Riemann-solver upwind) method with a slope limiter.
PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo *info, PetscReal t,
        PetscReal *aq, PetscReal *aG, void *ctx) {
    PetscErrorCode   ierr;
    ProblemCtx       *user = (ProblemCtx*)ctx;
    const PetscInt   n = info->dof;
    const PetscReal  hx = (user->b_right - user->a_left) / info->mx;
    Vec              sig;
    PetscInt         j, k;
    PetscReal        x, *asig, sl, sr,
                     *Ql, *Qr,  // slope-limited values of solution at either
                                //     side of current face
                     *Fl, *Fr;  // face fluxes on either end of current cell

    // for each owned cell get slope using the slope-limiter
    ierr = DMCreateLocalVector(info->da,&sig); CHKERRQ(ierr);
    ierr = VecSet(sig,0.0); CHKERRQ(ierr);  // implements limiter == NONE
    ierr = DMDAVecGetArray(info->da,sig,&asig); CHKERRQ(ierr);
    if (limiter != NONE) {  // following block assumes swidth >= 2
        for (j = info->xs-1; j < info->xs + info->xm+1; j++) {   // x_j is cell center
            for (k = 0; k < n; k++) {
                if ((j < 0 || j > info->mx-1) && user->periodic_bcs == PETSC_FALSE)
                    continue;
                if (j == 0 && user->periodic_bcs == PETSC_FALSE)
                    asig[n*j+k] = (aq[n*(j+1) + k] - aq[n*j + k]) / hx;
                else if (j == info->mx-1 && user->periodic_bcs == PETSC_FALSE)
                    asig[n*j+k] = (aq[n*j + k] - aq[n*(j-1) + k]) / hx;
                else {
                    if (limiter == FROMM) {
                        asig[n*j+k] = (aq[n*(j+1) + k] - aq[n*(j-1) + k]) / (2.0 * hx);
                    } else {
                        sr = (aq[n*(j+1) + k] - aq[n*j + k]) / hx;
                        sl = (aq[n*j + k] - aq[n*(j-1) + k]) / hx;
                        if (limiter == MINMOD)
                            asig[n*j+k] = minmod(sl,sr);
                        else if (limiter == MC)
                            asig[n*j+k] = minmod(0.5*(sl+sr),2.0*minmod(sl,sr));
                        else {
                            SETERRQ(PETSC_COMM_SELF,1,"how did I get here?\n");
                        }
                    }
                }
            }
        }
    }

    // get left-face flux Fl for first cell owned by process; may be at x=a
    ierr = PetscMalloc4(n,&Ql,n,&Qr,n,&Fl,n,&Fr); CHKERRQ(ierr);
    if (info->xs == 0 && user->periodic_bcs == PETSC_FALSE) {
        // use right slope
        slopemodify(n,-hx/2.0,&asig[0],&aq[0],Qr);
        ierr = user->bdryflux_a(t,hx,Qr,Fl); CHKERRQ(ierr);
    } else {
        // use left and right (limited) slope [left is owned by other process]
        x = user->a_left + (info->xs+0.5) * hx;
        slopemodify(n,hx/2.0,&asig[n*(info->xs-1)],&aq[n*(info->xs-1)],Ql);
        slopemodify(n,-hx/2.0,&asig[n*(info->xs)],&aq[n*(info->xs)],Qr);
        ierr = user->faceflux(t,x-hx/2.0,Ql,Qr,Fl); CHKERRQ(ierr);
    }

    // for each owned cell, compute RHS  G(t,x,q)
    for (j = info->xs; j < info->xs + info->xm; j++) {   // x_j is cell center
        x = user->a_left + (j+0.5) * hx;
        // set aG[n j + k] = g(t,x_j,u)_k
        ierr = user->g_source(t,x,&aq[n*j],&aG[n*j]); CHKERRQ(ierr);
        // get right-face flux Fr for cell; may be at x=b
        if (j == info->mx-1 && user->periodic_bcs == PETSC_FALSE) {
            // user left slope
            slopemodify(n,hx/2.0,&asig[n*j],&aq[n*j],Ql);
            ierr = user->bdryflux_b(t,hx,Ql,Fr); CHKERRQ(ierr);
        } else {
            // use left and right (limited) slope; see formulas LeVeque page 193
            slopemodify(n,hx/2.0,&asig[n*j],&aq[n*j],Ql);
            slopemodify(n,-hx/2.0,&asig[n*(j+1)],&aq[n*(j+1)],Qr);
            ierr = user->faceflux(t,x+hx/2.0,Ql,Qr,Fr); CHKERRQ(ierr);
        }
        // complete the RHS:
        //   aG[n j + k] = g(t,x_j,u)_k + (F_{j-1/2}_k - F_{j+1/2}_k) / hx
        for (k = 0; k < n; k++)
            aG[n*j+k] += (Fl[k] - Fr[k]) / hx;
        ncopy(n,Fr,Fl); // transfer Fr to Fl, for next loop
    }

    // clean up
    ierr = PetscFree4(Ql,Qr,Fl,Fr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(info->da,sig,&asig); CHKERRQ(ierr);
    ierr = VecDestroy(&sig); CHKERRQ(ierr);
    return 0;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <QuEST.h>
#include "QuEST_debug.h"
#include "QuEST_tn.h"
#include <cblas.h>

#define DEBUG 0

#if DEBUG
# define DEBUG_PRINT(x) printf x
#else
# define DEBUG_PRINT(x) do {} while (0)
#endif

int* getTensorIndexPermutation(int* contractionIndices, int numContractions,
        int* freeIndices, int numFreeIndices, int tensor){

    int *perm = (int*) malloc((numContractions+numFreeIndices)*sizeof(int));

    int firstSize, secondSize;
    int *firstArray, *secondArray;
    if (tensor==1){
        firstSize=numFreeIndices;
        firstArray=freeIndices;
        secondSize=numContractions;
        secondArray=contractionIndices;
    } else {
        secondSize=numFreeIndices;
        secondArray=freeIndices;
        firstSize=numContractions;
        firstArray=contractionIndices;
    }

    for (int i=0; i<firstSize; i++){
        perm[i]=firstArray[i];
    }
    for (int i=0; i<secondSize; i++){
        perm[i+firstSize]=secondArray[i];
    }
    return perm;
}

/*
 * All indices are of size 2
 */
int* getTensorSizes(int numQubits){
    int *sizes = (int*) malloc(numQubits*sizeof(int));
    for (int i=0; i<numQubits; i++) sizes[i]=2;
    return sizes;
}

Tensor createTensor(int numPq, int numVq, QuESTEnv env){
    Tensor tensor;

    //! Probably don't need to store these any more
    tensor.numPq = numPq;
    tensor.numVq = numVq;

    tensor.qureg = createQureg(tensor.numPq + tensor.numVq, env);
    //! might need to store next available vq

    return tensor;
}

/*
 * Swap corresponding qubit IDs in perm and return new index
 */
int swapBits(long long int b, int *perm, int nQubits) {
  long long int out = 0;
  for(int j=0;j<nQubits; j++) {
    out += (b >> j & 1) * (1 << perm[j]);
  }
  return out;
}

qreal* permuteArray(Qureg qureg, int *perm) {
  qreal* outArr = (qreal *) malloc(sizeof(qreal)*2*qureg.numAmpsPerChunk);
  for (long long int i = 0; i < qureg.numAmpsPerChunk; i++) {
    long long int newIndex = 2*swapBits(i, perm, qureg.numQubitsRepresented);
    outArr[newIndex] = qureg.stateVec.real[i];
    outArr[newIndex+1] = qureg.stateVec.imag[i];
  }
  return outArr;
}

/*
 * Expects tensor1Contractions to be in order from smallest to largest.
 * Expects tensor2Contractions to be ordered to match indices in tensor1.
 * Ie if tensor1 and tensor2 each have 3 indices and (tensor1, index 1) is
 * contracted with (tensor2, index 2) and (tensor1, index 2) is contracted with
 * (tensor2, index 0), use:
 * tensor1Contractions = [1, 2]
 * tensor2Contractions = [2, 0]
 */
Tensor contractIndices(Tensor tensor1, Tensor tensor2,
        int *tensor1Contractions, int *tensor2Contractions, int numContractions,
        int *tensor1FreeIndices, int numTensor1FreeIndices,
        int *tensor2FreeIndices, int numTensor2FreeIndices,
        QuESTEnv env){

    printf("Begin contracting\n");
    int numTensor1Qubits, numTensor2Qubits;
    numTensor1Qubits = tensor1.qureg.numQubitsRepresented;
    numTensor2Qubits = tensor2.qureg.numQubitsRepresented;
    int totalNumQ = numTensor1FreeIndices + numTensor2FreeIndices;

    printf("Permute indices to transform tensor contraction into the form of a matrix matrix multiply\n");
    int *tensor1Perm = getTensorIndexPermutation(tensor1Contractions, numContractions,
						 tensor1FreeIndices, numTensor1FreeIndices, 1);
    int *tensor2Perm = getTensorIndexPermutation(tensor2Contractions, numContractions,
						 tensor2FreeIndices, numTensor2FreeIndices, 2);

    qreal* tensor1StateVecPermuted = permuteArray(tensor1.qureg, tensor1Perm);
    qreal* tensor2StateVecPermuted = permuteArray(tensor2.qureg, tensor2Perm);

    // Working output array for BLAS routine. This needs to be separate to contractedQureg as
    // that represents complex numbers as a struct of arrays where BLAS will output an array of complex types
    qreal* outputQureg = (qreal *) malloc(sizeof(qreal)*(1LL << (totalNumQ+1LL)));

    // Dimensions of MM multiply
    int M, N, K;
    M = 1 << numTensor2FreeIndices;
    N = 1 << numTensor1FreeIndices;
    K = 1 << numContractions;
    qreal alpha[2] = {1.0, 0.0};
    qreal beta[2] = {0.0, 0.0};

    printf("MM multiply dimensions: M:%d N:%d K:%d\n", M, N, K);

    printf("Do MM multiply\n");
    // tensor2 x tensor1.
    cblas_zgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans, M, N, K, alpha,
    		tensor2StateVecPermuted, K, tensor1StateVecPermuted, N, beta, outputQureg, N);

    // Copy output into qureg object
    Qureg contractedQureg = createQureg(totalNumQ, env);
    for(long long int index = 0; index < contractedQureg.numAmpsPerChunk; index++) {
      contractedQureg.stateVec.real[index] = outputQureg[2*index];
      contractedQureg.stateVec.imag[index] = outputQureg[2*index+1];
    }

    printf("Free memory\n");
    free(tensor1StateVecPermuted);
    free(tensor2StateVecPermuted);
    free(outputQureg);

    free(tensor1Perm);
    free(tensor2Perm);

    // Free old quregs
    destroyQureg(tensor1.qureg, env);
    destroyQureg(tensor2.qureg, env);

    // Create output tensor
    Tensor outputTensor;
    outputTensor.qureg = contractedQureg;
    //! TODO: We may need to set total number of qubits here
    //outputTensor.numPq = totalNumPq;
    //outputTensor.numVq = totalNumVq;

    outputTensor.qureg = contractedQureg;
    return outputTensor;
}

// ----- operations ------------------------------------------------------------

/** Place target virtual qubit in the zero state
 * @param[in,out] tensor the tensor object
 */
void initVirtualTarget(Tensor tensor, int virtualTargetIndex){
         int vqIndex = virtualTargetIndex + tensor.numPq;
     printf("vqIndex: %d\n", vqIndex);

    Qureg qureg = tensor.qureg;

    long long int stateVecSize;
    long long int index;

    // dimension of the state vector
    // TODO: This won't work in parallel
    stateVecSize = 1LL << vqIndex;

    qreal *stateVecReal = qureg.stateVec.real;
    qreal *stateVecImag = qureg.stateVec.imag;

# ifdef _OPENMP
# pragma omp parallel \
    default  (none) \
    shared   (stateVecSize, stateVecReal, stateVecImag) \
    private  (index)
# endif
    {
# ifdef _OPENMP
# pragma omp for schedule (static)
# endif
        for (index=0; index<stateVecSize; index++) {
            stateVecReal[index+stateVecSize] = 0;
            stateVecImag[index+stateVecSize] = 0;
        }
    }
}

/** Initialize virtual qubit to |0> + |1>
 * NOTE: not (1/sqrt(2))(|0> + |1>)
 * @param[in,out] tensor the tensor object
 * @param[in] virtual qubit to initialize. Index is local to a tensor but includes all physical qubits in the tensor
 */
void initVirtualControl(Tensor tensor, int virtualControlIndex){
    int vqIndex = virtualControlIndex + tensor.numPq;
    Qureg qureg = tensor.qureg;

    long long int stateVecSize;
    long long int index;

    // dimension of the state vector
    // TODO: This won't work in parallel
    stateVecSize = 1LL << vqIndex;

    qreal *stateVecReal = qureg.stateVec.real;
    qreal *stateVecImag = qureg.stateVec.imag;

# ifdef _OPENMP
# pragma omp parallel \
    default  (none) \
    shared   (stateVecSize, stateVecReal, stateVecImag) \
    private  (index)
# endif
    {
# ifdef _OPENMP
# pragma omp for schedule (static)
# endif
        for (index=0; index<stateVecSize; index++) {
            stateVecReal[index+stateVecSize] = stateVecReal[index];
            stateVecImag[index+stateVecSize] = stateVecImag[index];
        }
    }
}

max_stars_repo_name: aniabrown/QuEST-TN

max_issues_repo_name: aniabrown/QuEST-TN

max_forks_repo_name: aniabrown/QuEST-TN


text:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/block/gsl_block.h>

#define BASE_DOUBLE
#include <gsl/templates_on.h>
#include <gsl/block/fwrite_source.c>
#include <gsl/block/fprintf_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:
#pragma once

#include "file_descriptor.h"
#include <filesystem>
#include <gsl/span>
#include <random>

namespace dogbox
{
    inline void create_randomness(gsl::span<std::byte> const into)
    {
        file_descriptor const random = open_file_for_reading("/dev/urandom").value();
        ptrdiff_t written = 0;
        while (written < into.size())
        {
            size_t const reading = static_cast<size_t>(into.size() - written);
            ssize_t const read_result = read(random.handle, into.data() + written, reading);
            if (read_result < 0)
            {
                TO_DO();
            }
            written += reading;
        }
    }
}

max_stars_repo_name: TyRoXx/dogbox

max_issues_repo_name: TyRoXx/dogbox

max_forks_repo_name: TyRoXx/dogbox


text:
#include "jobqueue.h"
#include "degnome.h"
#include "fitfunc.h"
#include "flagparse.c"
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <time.h>
#include <pthread.h>
#include <unistd.h>
#include <limits.h>

typedef struct JobData JobData;
struct JobData {
	Degnome* child;
	Degnome* p1;
	Degnome* p2;
};

void usage(void);
void help_menu(void);
int jobfunc(void* p, void* tdat);

const char* usageMsg =
	"Usage: polygensim [-h] [-c chromosome_length] [-e mutation_effect]\n"
	"\t\t  [-g num_generations] [-m mutation_rate]\n"
	"\t\t  [-o crossover_rate] [-p population_size]\n"
	"\t\t  [-t num_threads] [--seed rngseed]\n"
	"\t\t  [--target hat_height target]\n"
	"\t\t  [--sqrt | --linear | --close | --ceiling | --log]\n";

const char* helpMsg =
	"OPTIONS\n"
	"\t -c chromosome_length\n"
	"\t\t Set chromosome length for the current simulation.\n"
	"\t\t Default chromosome length is 10.\n\n"
	"\t -e mutation_effect\n"
	"\t\t Set how much a mutation will effect a gene on average.\n"
	"\t\t Default mutation effect is 2.\n\n"
	"\t -g num_generations\n"
	"\t\t Set how many generations this simulation will run for.\n"
	"\t\t Default number of generations is 1000.\n\n"
	"\t -h\t Display this help menu.\n\n"
	"\t -m mutation_rate\n"
	"\t\t Set the mutation rate for the current simulation.\n"
	"\t\t Default mutation rate is 1.\n\n"
	"\t -o crossover_rate\n"
	"\t\t Set the crossover rate for the current simulation.\n"
	"\t\t Default crossover rate is 2.\n\n"
	"\t -p population_size\n"
	"\t\t Set the population size for the current simulation.\n"
	"\t\t Default population size is 10.\n"
	"\t -t num_threads\n"
	"\t\t Select the number of threads to be used in the current run.\n"
	"\t\t Default is 0 (which will result in 3/4 of cores being used).\n"
	"\t\t Must be 1 if a seed is used in order to prevent race conditions.\n\n"
	"\t --seed rngseed\n"
	"\t\t Select the seed used by the RNG in the current run.\n"
	"\t\t Default seed is 0 (which will result in a random seed).\n\n"
	"\t --target hat_height target\n"
	"\t\t Sets the ideal hat height for the current simulation\n"
	"\t\t Used for fitness functions that have an \"ideal\" value.\n\n"
	"\t --sqrt\t\t fitness will be sqrt(hat_height)\n\n"
	"\t --linear\t fitness will be hat_height\n\n"
	"\t --close\t fitness will be (target - abs(target - hat_height))\n\n"
	"\t --ceiling\t fitness will quickly level off after passing target\n\n";

pthread_mutex_t seedLock = PTHREAD_MUTEX_INITIALIZER;
unsigned long rngseed = 0;

//	there is no need for the line 'int chrom_size' as it is declared as a global variable in degnome.h
int pop_size;
int num_gens;
int mutation_rate;
int mutation_effect;
int crossover_rate;

int num_threads = 0;
JobQueue* jq;

void *ThreadState_new(void *notused);
void ThreadState_free(void *rng);

void *ThreadState_new(void *notused) {
	// Lock seed, initialize random number generator, increment seed,
	// and unlock.
	gsl_rng *rng = gsl_rng_alloc(gsl_rng_taus);

	pthread_mutex_lock(&seedLock);
	gsl_rng_set(rng, rngseed);
	rngseed = (rngseed == ULONG_MAX ? 0 : rngseed + 1);
	pthread_mutex_unlock(&seedLock);

	return rng;
}

void ThreadState_free(void *rng) {
	gsl_rng_free((gsl_rng *) rng);
}

int jobfunc(void* p, void* tdat) {
	gsl_rng* rng = (gsl_rng*) tdat;
	JobData* data = (JobData*) p;																				//get data out
	Degnome_mate(data->child, data->p1, data->p2, rng, mutation_rate, mutation_effect, crossover_rate);			//mate

	return 0;		//exited without error
}

void usage(void) {
	fputs(usageMsg, stderr);
	exit(EXIT_FAILURE);
}

void help_menu(void) {
	fputs(helpMsg, stderr);
	exit(EXIT_FAILURE);
}

int main(int argc, char **argv) {

	int * flags = NULL;

	if (parse_flags(argc, argv, 1, &flags) == -1) {
		free(flags);
		usage();
	}

	if (flags[2] == 1) {
        free(flags);
		help_menu();
	}

	chrom_size = flags[6];
	mutation_effect = flags[7];
	num_gens = flags[8];
	mutation_rate = flags[9];
	crossover_rate = flags[10];
	pop_size = flags[11];

	num_threads = flags[12];

	if(flags[13] == 0){
		set_function("linear");
	}
	else if(flags[13] == 1){
		set_function("sqrt");
	}
	else if(flags[13] == 2){
		set_function("close");
	}
	else if(flags[13] == 3){
		set_function("ceiling");
	}
	else if(flags[13] == 4){
		set_function("log");
	}
	target_num = flags[14];

	if(flags[15] <= 0) {
		time_t currtime = time(NULL);                  // time
		unsigned long pid = (unsigned long) getpid();  // process id
		rngseed = currtime ^ pid;                      // random seed
	}
	else{
		rngseed = flags[15];
	}
	gsl_rng* rng = gsl_rng_alloc(gsl_rng_taus);    // rand generator
	gsl_rng_set(rng, rngseed);

	free(flags);

	if (num_threads <= 0) {
		if (num_threads < 0) {
			#ifdef DEBUG_MODE
				fprintf(stderr, "Error invalid number of threads: %u\n", num_threads);
			#endif
		}
		num_threads = (3*getNumCores()/4);
	}
	#ifdef DEBUG_MODE
		fprintf(stderr, "Final number of threads: %u\n", num_threads);
	#endif

	Degnome* parents;
	Degnome* children;
	Degnome* temp;

	printf("%u, %u, %u\n", chrom_size, pop_size, num_gens);

	parents = malloc(pop_size*sizeof(Degnome));
	children = malloc(pop_size*sizeof(Degnome));

	for (int i = 0; i < pop_size; i++) {
		parents[i].dna_array = malloc(chrom_size*sizeof(double));
		children[i].dna_array = malloc(chrom_size*sizeof(double));
		parents[i].hat_size = 0;

		for (int j = 0; j < chrom_size; j++) {
			parents[i].dna_array[j] = (i+j);	//children isn't initiilized
			parents[i].hat_size += (i+j);
		}
	}

	printf("Generation 0:\n");
	for (int i = 0; i < pop_size; i++) {
		printf("Degnome %u\n", i);
		for (int j = 0; j < chrom_size; j++) {
			printf("%lf\t", parents[i].dna_array[j]);
		}
		printf("\nTOTAL HAT SIZE: %lg\n\n", parents[i].hat_size);
	}

	jq = JobQueue_new(num_threads, NULL, ThreadState_new, ThreadState_free);

	JobData* dat = malloc(pop_size*sizeof(JobData));

	for (int i = 0; i < num_gens; i++) {

		double fit = get_fitness(parents[0].hat_size);

		double total_hat_size = fit;
		double cum_hat_size[pop_size];
		cum_hat_size[0] = fit;

		for (int j = 1; j < pop_size; j++) {
			fit = get_fitness(parents[j].hat_size);

			total_hat_size += fit;
			cum_hat_size[j] = (cum_hat_size[j-1] + fit);
		}

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

			int m, d;

			double win_m = gsl_rng_uniform(rng);
			win_m *= total_hat_size;
			double win_d = gsl_rng_uniform(rng);
			win_d *= total_hat_size;

			// printf("win_m:%lf, wind:%lf, max: %lf\n", win_m,win_d,total_hat_size);

			for (m = 0; cum_hat_size[m] < win_m; m++) {
				continue;
			}

			for (d = 0; cum_hat_size[d] < win_d; d++) {
				continue;
			}

			// printf("m:%u, d:%u\n", m,d);
			dat[j].child = (children + j);
			dat[j].p1 = (parents + m);
			dat[j].p2 = (parents + d);

			JobQueue_addJob(jq, jobfunc, dat + j);
		}

		JobQueue_waitOnJobs(jq);
		temp = children;
		children = parents;
		parents = temp;
	}

	JobQueue_noMoreJobs(jq);

	printf("Generation %u:\n", num_gens);
	for (int i = 0; i < pop_size; i++) {
		printf("Degnome %u\n", i);
		for (int j = 0; j < chrom_size; j++) {
			printf("%lf\t", parents[i].dna_array[j]);
		}
		printf("\nTOTAL HAT SIZE: %lg\n\n", parents[i].hat_size);
	}

	//free everything

	JobQueue_free(jq);
	free(dat);

	for (int i = 0; i < pop_size; i++) {
		free(parents[i].dna_array);
		free(children[i].dna_array);
		parents[i].hat_size = 0;
	}

	free(parents);
	free(children);

	gsl_rng_free (rng);
}

max_stars_repo_name: masalemi/PopGenSim

max_issues_repo_name: masalemi/PopGenSim

max_forks_repo_name: masalemi/PopGenSim


text:
#pragma once

#include <memory>

#include <gsl/gsl_rng.h>

extern std::unique_ptr<gsl_rng, decltype(&gsl_rng_free)> rng;

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include "arcana/type_traits.h"
#include "arcana/utils/serialization/base_stream.h"
#include "arcana/utils/serialization/serializable.h"

#include <gsl/gsl>
#include <memory>
#include <type_traits>

namespace mira
{
    template<typename LambdaT>
    class binary_iterator
    {
        template<typename T>
        using is_serializable = std::is_base_of<serializable<std::remove_cv_t<T>>, T>;

        template<typename L>
        struct check_custom_iterator
        {
            template<typename T,
                     typename = decltype(binary_iterate(instance_of<const T>(),
                                                        instance_of<binary_iterator<LambdaT>>()))>
            static std::true_type test(T* = nullptr);
            static std::false_type test(...);
        };

        template<typename T>
        using has_custom_iterator = decltype(check_custom_iterator<LambdaT>::test(&instance_of<T>()));

        template<typename T>
        using is_valid_pod =
            std::integral_constant<bool,
                                   !has_custom_iterator<T>::value && std::is_pod<std::decay_t<T>>::value &&
                                       !std::is_pointer<std::decay_t<T>>::value>;

    public:
        template<typename... ArgTs>
        void iterate_all(ArgTs&&... args)
        {
            int unused[] = { (iterate(std::forward<ArgTs>(args)), 0)... };
            (void)unused;
        }

        template<typename T>
        std::enable_if_t<has_custom_iterator<std::decay_t<T>>::value> iterate(const T& value)
        {
            static_assert(!is_valid_pod<T>::value, "we've got some overlap");
            binary_iterate(value, *this);
        }

        template<typename T>
        std::enable_if_t<is_valid_pod<std::decay_t<T>>::value> iterate(const T& ptr)
        {
            static_assert(!has_custom_iterator<T>::value, "we've got some overlap");
            stream{ *this }.write(ptr);
        }

        void iterate(const void* data, size_t size)
        {
            m_iter(data, size);
        }

        template<typename ItrT>
        void iterate(ItrT beg, ItrT end)
        {
            while (beg != end)
            {
                iterate(*beg);
                beg++;
            }
        }

        template<typename A, typename B>
        void iterate(const std::pair<A, B>& pair)
        {
            iterate(pair.first);
            iterate(pair.second);
        }

        template<typename T>
        void iterate(const std::shared_ptr<T>& ptr)
        {
            iterate(ptr.get());
        }

        template<typename T, ptrdiff_t E = -1>
        void iterate(const gsl::span<T, E>& span)
        {
            iterate(span.begin(), span.end());
        }

        template<typename T, typename D>
        void iterate(const std::unique_ptr<T, D>& ptr)
        {
            iterate(ptr.get());
        }

        template<typename T>
        void iterate(const T* ptr)
        {
            if (!ptr)
                iterate((intptr_t)0);
            else
                iterate(*ptr);
        }

        template<typename T>
        void iterate(const serializable<T>& ptr)
        {
            stream s{ *this };
            ptr.serialize(s);
        }

        template<typename ItrT>
        explicit binary_iterator(ItrT&& i)
            : m_iter{ std::forward<ItrT>(i) }
        {}

        LambdaT& iterator()
        {
            return m_iter;
        }

        const LambdaT& iterator() const
        {
            return m_iter;
        }

    private:
        struct stream : public base_stream
        {
            template<typename T>
            void write(const T& obj)
            {
                m_itr.iterate(&obj, sizeof(T));
            }

            template<typename T>
            void write(const T* obj, size_t N)
            {
                m_itr.iterate(obj, N * sizeof(T));
            }

            stream(binary_iterator<LambdaT>& itr)
                : m_itr{ itr }
            {}

        private:
            binary_iterator<LambdaT>& m_itr;
        };

        LambdaT m_iter;
    };

    template<typename LambdaT>
    binary_iterator<std::decay_t<LambdaT>> make_binary_iterator(LambdaT&& lambda)
    {
        return binary_iterator<std::decay_t<LambdaT>>{ std::forward<LambdaT>(lambda) };
    }
}

namespace std
{
    template<typename Iter>
    void binary_iterate(const std::string& str, mira::binary_iterator<Iter>& iter)
    {
        iter.iterate(str.data(), str.size());
    }

    template<typename T, typename Iter>
    void binary_iterate(const std::vector<T>& vec, mira::binary_iterator<Iter>& iter)
    {
        iter.iterate(vec.size());
        iter.iterate(vec.begin(), vec.end());
    }

    template<typename IterT, typename FirstT, typename SecondT>
    void binary_iterate(const std::pair<FirstT, SecondT>& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate_all(data.first, data.second);
    }

    template<typename Iter, typename T>
    void binary_iterate(const std::shared_ptr<T>& ptr, mira::binary_iterator<Iter>& iter)
    {
        iter.iterate(ptr.get());
    }

    template<typename Iter, typename T>
    void binary_iterate(const std::unique_ptr<T>& ptr, mira::binary_iterator<Iter>& iter)
    {
        iter.iterate(ptr.get());
    }
}

namespace gsl
{
    template<typename Iter, typename T, ptrdiff_t E = -1>
    void binary_iterate(const gsl::span<T, E>& span, mira::binary_iterator<Iter>& iter)
    {
        iter.iterate(span.begin(), span.end());
    }
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam