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

#include "core_allvars_grid.h"
#include "core_proto_grid.h"

int32_t init(void)
{
  int32_t i, status;

  sum_photons = 0.0;

  random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);
  gsl_rng_set(random_generator, 42);	 // start-up seed

  set_units(); // Change units to code units where necessary.
  srand((unsigned) time(NULL));

  read_snap_list(); // Read the snapshots required.

  for(i = 0; i < Snaplistlen; i++)
  {
    ZZ[i] = 1 / AA[i] - 1; // Redshift array.
    Age[i] = time_to_present(ZZ[i]); // Age array.
  }

  if (fescPrescription == 0)
  {
    printf("\n\nUsing a constant escape fraction of %.4f\n", fesc);
  }
  else if (fescPrescription == 1)
  {
    fprintf(stderr, "\n\nDeprecated! Use a new one!\n");
    return EXIT_FAILURE;
  }
  else if (fescPrescription == 2)
  {
    printf("\n\nUsing an fesc prescription that scales with halo mass.\n\n");
    determine_fesc_constants();
  }
  else if (fescPrescription == 3)
  {
    printf("\n\nUsing an fesc prescription that scales with the fraction of ejected mass in the galaxy.\n");
    printf("\nThis takes the form A*fej + B with A = %.4e and B = %.4e\n", alpha, beta);
  }
  else if (fescPrescription == 4)
  {
    printf("\n\nUsing an fesc prescription that depends upon quasar activity.\n");
    printf("\nFor a galaxy that had a quasar event within %.2f dynamical times go the escape fraction will be %.2f.  Otherwise it will have a constant value of %.2f\n", N_dyntime, quasar_boosted, quasar_baseline);
  }
  else if (fescPrescription == 5)
  {
    printf("\n\nUsing Anne's functional form for an escape fraction that decreases for increasing halo mass.\n");
    XASSERT(fesc_low > fesc_high, "Input file contain fesc_low = %.2f and fesc_high = %.2f. For this prescription (fescPrescription == 5), we require fesc_low > fesc_high\n", fesc_low, fesc_high);
  }
  else if (fescPrescription == 6)
  {
    printf("\n\nUsing Anne's functional form for an escape fraction that increases for increasing halo mass.\n");
    XASSERT(fesc_low < fesc_high, "Input file contain fesc_low = %.2f and fesc_high = %.2f. For this prescription (fescPrescription == 6), we require fesc_low < fesc_high\n", fesc_low, fesc_high);
  }
  else if (fescPrescription == 7)
  {
    printf("\n\nUsing an fesc prescription that scales as a power law with the fraction of ejected gas.\n");
    determine_fesc_constants();
  }
  else
  {
    printf("\n\nOnly escape fraction prescriptions 0 to 6 (exlucding 1) are permitted.\n");
    return EXIT_FAILURE;
  }

  Grid = malloc(sizeof(struct GRID_STRUCT));
  if (Grid == NULL)
  {
    fprintf(stderr, "Could not allocate memory for the high level grid structure\n");
    return EXIT_FAILURE;
  }

  status = init_grid(Grid);
  if (status == EXIT_FAILURE)
  {
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;

  QuasarEventsAbovePartCut = 0;
  QuasarEventsBelowPartCut = 0;

}

int32_t init_grid(struct GRID_STRUCT *grid)
{

#define ALLOCATE_GRID_MEMORY(name, length) \
{                                          \
  name = calloc(length, sizeof(*(name)));  \
  if (name == NULL)                        \
  {                                        \
    fprintf(stderr, "Out of memory allocating %ld bytes, could not allocate"#name".\n", sizeof(*(name)* length)); \
    return EXIT_FAILURE;                   \
  }                                        \
}

  int32_t i;
  uint64_t cell_idx;

  if (grid == NULL)
  {
    fprintf(stderr, "init_grid was called with a GRID_STRUCT pointer that has not been initialized\n");
    return EXIT_FAILURE;
  }

  grid->GridSize = GridSize;
  grid->NumCellsTotal = CUBE(GridSize);

  grid->NumGrids = NGrid;

  grid->GridProperties = malloc(sizeof(struct GRID_PROPERTIES_STRUCT) * grid->NumGrids);
  if (grid->GridProperties == NULL)
  {
    fprintf(stderr, "Could not allocate memory for the grid properties struct\n");
    return EXIT_FAILURE;
  }

  for (i = 0; i < grid->NumGrids; ++i)
  {
#ifdef MPI
    printf("Allocating grid %d on Task %d.\n", i, ThisTask);
#endif

    /* For now the self-consistent model only using the number of ionizing photons. */
    /* Also use the number of galaxies for sanity check. */

    ALLOCATE_GRID_MEMORY(grid->GridProperties[i].Nion_HI, grid->NumCellsTotal);
    ALLOCATE_GRID_MEMORY(grid->GridProperties[i].GalCount, grid->NumCellsTotal);

    if (self_consistent == 0)
    {
      ALLOCATE_GRID_MEMORY(grid->GridProperties[i].SFR, grid->NumCellsTotal);
      ALLOCATE_GRID_MEMORY(grid->GridProperties[i].StellarMass, grid->NumCellsTotal);
      ALLOCATE_GRID_MEMORY(grid->GridProperties[i].Nion_HeI, grid->NumCellsTotal);
      ALLOCATE_GRID_MEMORY(grid->GridProperties[i].Nion_HeII, grid->NumCellsTotal);

    }

    // All the memory has been allocated for the inner grid.
    // Now let's initialize the values.

    for (cell_idx = 0; cell_idx < Grid->NumCellsTotal; ++cell_idx)
    {
      grid->GridProperties[i].Nion_HI[cell_idx] = 0.0;
      grid->GridProperties[i].GalCount[cell_idx] = 0;

      if (self_consistent == 0)
      {
        grid->GridProperties[i].SFR[cell_idx] = 0.0;
        grid->GridProperties[i].StellarMass[cell_idx] = 0.0;
        grid->GridProperties[i].Nion_HeI[cell_idx] = 0.0;
        grid->GridProperties[i].Nion_HeII[cell_idx] = 0.0;

      }
    }

  }

  printf("All grids have been initialized\n");
  return EXIT_SUCCESS;
}

int32_t free_grid(void)
{

  int32_t i;

  for (i = 0; i < Grid->NumGrids; ++i)
  {
    free(Grid->GridProperties[i].Nion_HI);
    free(Grid->GridProperties[i].GalCount);

    if (self_consistent == 0)
    {
      free(Grid->GridProperties[i].Nion_HeII);
      free(Grid->GridProperties[i].Nion_HeI);
      free(Grid->GridProperties[i].StellarMass);
      free(Grid->GridProperties[i].SFR);
    }
  }

  free(Grid->GridProperties);
  free(Grid);

  return EXIT_SUCCESS;

}

void set_units(void)
{

  UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s;
  UnitTime_in_Megayears = UnitTime_in_s / SEC_PER_MEGAYEAR;
  G = GRAVITY / pow(UnitLength_in_cm, 3) * UnitMass_in_g * pow(UnitTime_in_s, 2);
  UnitDensity_in_cgs = UnitMass_in_g / pow(UnitLength_in_cm, 3);
  UnitPressure_in_cgs = UnitMass_in_g / UnitLength_in_cm / pow(UnitTime_in_s, 2);
  UnitCoolingRate_in_cgs = UnitPressure_in_cgs / UnitTime_in_s;
  UnitEnergy_in_cgs = UnitMass_in_g * pow(UnitLength_in_cm, 2) / pow(UnitTime_in_s, 2);

  // convert some physical input parameters to internal units
  Hubble = HUBBLE * UnitTime_in_s;

  // compute a few quantitites
  RhoCrit = 3 * Hubble * Hubble / (8 * M_PI * G);

}

void read_snap_list(void)
{
  FILE *fd;
  char fname[MAXLEN];

  snprintf(fname, MAXLEN, "%s", FileWithSnapList);

  if(!(fd = fopen(fname, "r")))
  {
    printf("can't read output list in file '%s'\n", fname);
    ABORT(0);
  }

  Snaplistlen = 0;
  do
  {
    if(fscanf(fd, " %lg ", &AA[Snaplistlen]) == 1)
      Snaplistlen++;
    else
      break;
  }
  while(Snaplistlen < MAXSNAPS);

  fclose(fd);

#ifdef MPI
  if(ThisTask == 0)
#endif
  printf("found %d defined times in snaplist\n", Snaplistlen);
}

double time_to_present(double z)
{
#define WORKSIZE 1000
  gsl_function F;
  gsl_integration_workspace *workspace;
  double time, result, abserr;

  workspace = gsl_integration_workspace_alloc(WORKSIZE);
  F.function = &integrand_time_to_present;

  gsl_integration_qag(&F, 1.0 / (z + 1), 1.0, 1.0 / Hubble,
    1.0e-8, WORKSIZE, GSL_INTEG_GAUSS21, workspace, &result, &abserr);

  time = 1 / Hubble * result;

  gsl_integration_workspace_free(workspace);

  // return time to present as a function of redshift
  return time;
}

double integrand_time_to_present(double a, void *param)
{
  return 1 / sqrt(Omega / a + (1 - Omega - OmegaLambda) + OmegaLambda * a * a);
}

void determine_fesc_constants(void)
{

  double A, B, log_A;

  log_A = (log10(fesc_high) - (log10(fesc_low)*log10(MH_high)/log10(MH_low))) * pow(1 - (log10(MH_high) / log10(MH_low)), -1);
  B = (log10(fesc_low) - log_A) / log10(MH_low);
  A = pow(10, log_A);

  alpha = A;
  beta = B;

  printf("Fixing the points (%.4e, %.2f) and (%.4e, %.2f)\n", MH_low, fesc_low, MH_high, fesc_high);
  printf("This gives a power law with constants A = %.4e, B = %.4e\n", alpha, beta);

}
/*
void estimate_grid_memory(void)
{

  if (Verbose == 1)
  {
    printf("The size of the Grid Struct (that contains all the pointers) is %lu bytes\nEach grid cell has a struct, in addition to data that the pointers point to.\n", sizeof(struct GRID));
    printf("The total memory taken up by the Grid will be the sizeof a grid cell (%lu bytes) times the number of grid cells (%d) times the number of output snapshots (%d) times 2 because we need memory for both the pointers for the grid and the actual numbers themselves.\n", sizeof(struct GRID), CUBE(GridSize), NGrid);
  }

  printf("Approximate total memory for the grid ===== %lu mb\n", sizeof(struct GRID)*CUBE(GridSize)*NGrid*2/1024/1024);

}

void estimate_gal_memory(int NtotGals)
{
  if (Verbose == 1)
  {
    printf("The size of the Galaxy struct is %lu bytes and the size of the GalaxyGrid struct (that contains all the pointers for the SnapShot history) is %lu\n", sizeof(struct GALAXY_INPUT), sizeof(struct GALAXY_GRID));
    printf("The total memory taken up by this file of galaxies will be the sizeof the galaxy input (%lu bytes) times the number of galaxies (%d) plus the sizeof the galaxygrid (%lu bytes) times the number of snapshots in the simulation (%d) times the number of galaxies (%d)\n", sizeof(struct GALAXY_INPUT), NtotGals, sizeof(struct GALAXY_GRID), LastSnapShotNr, NtotGals);
  }

  printf("Approximate memory occupied for this file of galaxies ===== %lu mb\n", (sizeof(struct GALAXY_INPUT)*NtotGals + sizeof(struct GALAXY_GRID)*NtotGals*LastSnapShotNr)/1024/1024);

}
*/

max_stars_repo_name: jacobseiler/rsage

max_issues_repo_name: jacobseiler/rsage

max_forks_repo_name: jacobseiler/rsage


text:
#pragma once
#include "halley/maths/vector3.h"
#include <gsl/span>

namespace Halley
{
	class AudioListenerData;

	class AudioChannelData
	{
	public:
		float pan; // TODO, do this right
	};

	class AudioPosition
	{
	public:
		static AudioPosition makeUI(float pan = 0.0f); // -1.0f = left, 1.0f = right
		static AudioPosition makePositional(Vector2f pos, float referenceDistance = 200.0f, float maxDistance = 400.0f);
		static AudioPosition makePositional(Vector3f pos, float referenceDistance = 200.0f, float maxDistance = 400.0f);
		static AudioPosition makeFixed();

		void setMix(size_t srcChannels, gsl::span<const AudioChannelData> dstChannels, gsl::span<float, 16> dst, float gain, const AudioListenerData& listener) const;

		AudioPosition();

	private:
		Vector3f pos;
		float referenceDistance;
		float maxDistance;
		bool isUI;
		bool isPannable;
	};
}

ext: h

max_stars_repo_name: Healthire/halley

max_issues_repo_name: Healthire/halley

max_forks_repo_name: Healthire/halley


text:
#include <math.h>
#include <gsl/gsl_linalg.h>
#include "getNextIterate.h"

void getNextIterate(const grid *grd, const bool eos_func,
		    const wksp *w
#if AA_M > 0
		    , const unsigned long itCount, const long nHist,
		    double *residMax
#endif
		    ) {

  long i;
#if AA_M > 0
  long j, mk, colptr;
  gsl_matrix *resid;
  gsl_vector *tau, *llsresid, *aa_wgts, *constraint;
  gsl_vector_view constraint_view;
#endif

#if AA_M == 0
/* No Anderson acceleration, so just move temporaries into next
   guess */
  for (i=0; i<grd->nr; i++) {
    w->colNew[i] = w->colTmp[i];
    w->presNew_g[i+1] = gsl_vector_get(w->presTmp_g, i+1);
    if (eos_func==1)
      w->eIntNew[i] = w->eIntTmp[i];
  }
#else

  /* Anderson acceleration code */

  /* Set up pointers */
  mk = itCount < AA_M+1 ? itCount : AA_M+1;
  colptr = (itCount - 1) % (AA_M+1);

  /* Store the current guesses in the history array */
  for (i=0; i<grd->nr; i++) {
    w->colHist[i + colptr*grd->nr] = w->colTmp[i];
    w->presHist[i + colptr*grd->nr] = gsl_vector_get(w->presTmp_g, i+1);
    if (eos_func)
      w->eIntHist[i + colptr*grd->nr] = w->eIntTmp[i];
  }

  /* Compute residuals and store them in the residual matrix */
  for (i=0; i<grd->nr; i++) {
    w->colResid[i + colptr*grd->nr] =
      (w->colHist[i+colptr*grd->nr]-w->colNew[i]) /
      w->colHist[i+colptr*grd->nr];
    if (fabs(w->colResid[i + colptr*grd->nr]) > *residMax) {
      *residMax = fabs(w->colResid[i + colptr*grd->nr]);
    }
    w->presResid[i + colptr*grd->nr] =
      (w->presHist[i+colptr*grd->nr]-w->presNew_g[i+1]) /
      w->presHist[i+colptr*grd->nr];
    if (fabs(w->presResid[i + colptr*grd->nr]) > *residMax) {
      *residMax = fabs(w->presResid[i + colptr*grd->nr]);
    }
    if (eos_func) {
      w->eIntResid[i + colptr*grd->nr] =
	(w->eIntHist[i+colptr*grd->nr]-w->eIntNew[i]) /
	w->eIntHist[i+colptr*grd->nr];
      if (fabs(w->eIntResid[i + colptr*grd->nr]) > *residMax)
	*residMax = fabs(w->eIntResid[i + colptr*grd->nr]);
    }
  }

  /* If we're in the first iteration, set move temporaries to next
     guess */
  if (itCount == 1) {
    for (i=0; i<grd->nr; i++) {
      w->colNew[i] = w->colHist[i+colptr*grd->nr];
      w->presNew_g[i+1] = w->presHist[i+colptr*grd->nr];
      if (eos_func)
	w->eIntNew[i] = w->eIntHist[i+colptr*grd->nr];
    }
  } else {
    /* On subsequent iterations, solve the linear least squares
       optimization problem to get a weighted sum of previous
       guesses for the new guess */

    /* Allocate workspace for this iteration */
    resid = gsl_matrix_alloc(nHist+1, mk);
    llsresid = gsl_vector_alloc(nHist+1);
    aa_wgts = gsl_vector_alloc(mk);
    if (nHist+1 < mk) {
      tau = gsl_vector_alloc(nHist+1);
    } else {
      tau = gsl_vector_alloc(mk);
    }

    /* Build matrix of residuals */
    for (i=0; i<grd->nr; i++) {
      for (j=0; j<mk; j++) {
	gsl_matrix_set(resid, i, j, w->colResid[i+j*grd->nr]);
	gsl_matrix_set(resid, i+grd->nr, j, w->presResid[i+j*grd->nr]);
	if (eos_func)
	  gsl_matrix_set(resid, i+2*grd->nr, j, w->eIntResid[i+j*grd->nr]);
      }
    }

    /* Set up constraints on coefficients */
    for (j=0; j<mk; j++) gsl_matrix_set(resid, nHist, j, 1.0e3*(*residMax));
    constraint_view = gsl_vector_subvector(w->constraint, 0, nHist+1);
    constraint = &(constraint_view.vector);
    gsl_vector_set(constraint, nHist, 1e3*(*residMax));

    /* Solve linear least squares system by QR decomposition */
    gsl_linalg_QR_decomp(resid, tau);
    gsl_linalg_QR_lssolve(resid, tau, constraint, aa_wgts, llsresid);

    /* Build the new guess from the computed weights */
    for (i=0; i<grd->nr; i++) {
      w->colNew[i] =
	gsl_vector_get(aa_wgts, 0)*w->colHist[i];
      w->presNew_g[i+1] =
	gsl_vector_get(aa_wgts, 0)*w->presHist[i];
      if (eos_func)
	w->eIntNew[i] =
	  gsl_vector_get(aa_wgts, 0)*w->eIntHist[i];
    }
    for (j=1; j<mk; j++) {
      for (i=0; i<grd->nr; i++) {
	w->colNew[i] += gsl_vector_get(aa_wgts, j) *
	  w->colHist[i+j*grd->nr];
	w->presNew_g[i+1] += gsl_vector_get(aa_wgts, j) *
	  w->presHist[i+j*grd->nr];
	if (eos_func)
	  w->eIntNew[i] += gsl_vector_get(aa_wgts, j) *
	    w->eIntHist[i+j*grd->nr];
      }
    }

    /* Free memory for this iteration */
    gsl_vector_free(llsresid);
    gsl_vector_free(tau);
    gsl_vector_free(aa_wgts);
    gsl_matrix_free(resid);
  }

#endif
  /* End Anderson acceleration code */
}

ext: c

max_stars_repo_name: franciscaconcha/amuse-vader

max_issues_repo_name: franciscaconcha/amuse-vader

max_forks_repo_name: franciscaconcha/amuse-vader


text:
#include <petsc.h>
#include "poissonfunctions.h"

PetscErrorCode Poisson1DFunctionLocal(DMDALocalInfo *info, double *au,
                                   double *aF, PoissonCtx *user) {
    PetscErrorCode ierr;
    int          i;
    double       xmax[1], xmin[1], h, x, ue, uw;
    ierr = DMDAGetBoundingBox(info->da,xmin,xmax); CHKERRQ(ierr);
    h = (xmax[0] - xmin[0]) / (info->mx - 1);
    for (i = info->xs; i < info->xs + info->xm; i++) {
        x = xmin[0] + i * h;
        if (i==0 || i==info->mx-1) {
            aF[i] = au[i] - user->g_bdry(x,0.0,0.0,user);
            aF[i] *= user->cx * (2.0 / h);
        } else {
            ue = (i+1 == info->mx-1) ? user->g_bdry(x+h,0.0,0.0,user)
                                     : au[i+1];
            uw = (i-1 == 0)          ? user->g_bdry(x-h,0.0,0.0,user)
                                     : au[i-1];
            aF[i] = user->cx * (2.0 * au[i] - uw - ue) / h
                    - h * user->f_rhs(x,0.0,0.0,user);
        }
    }
    ierr = PetscLogFlops(9.0*info->xm);CHKERRQ(ierr);
    return 0;
}

//STARTFORM2DFUNCTION
PetscErrorCode Poisson2DFunctionLocal(DMDALocalInfo *info, double **au,
                                   double **aF, PoissonCtx *user) {
    PetscErrorCode ierr;
    int     i, j;
    double  xymin[2], xymax[2], hx, hy, darea, scx, scy, scdiag, x, y,
            ue, uw, un, us;
    ierr = DMDAGetBoundingBox(info->da,xymin,xymax); CHKERRQ(ierr);
    hx = (xymax[0] - xymin[0]) / (info->mx - 1);
    hy = (xymax[1] - xymin[1]) / (info->my - 1);
    darea = hx * hy;
    scx = user->cx * hy / hx;
    scy = user->cy * hx / hy;
    scdiag = 2.0 * (scx + scy);    // diagonal scaling
    for (j = info->ys; j < info->ys + info->ym; j++) {
        y = xymin[1] + j * hy;
        for (i = info->xs; i < info->xs + info->xm; i++) {
            x = xymin[0] + i * hx;
            if (i==0 || i==info->mx-1 || j==0 || j==info->my-1) {
                aF[j][i] = au[j][i] - user->g_bdry(x,y,0.0,user);
                aF[j][i] *= scdiag;
            } else {
                ue = (i+1 == info->mx-1) ? user->g_bdry(x+hx,y,0.0,user)
                                         : au[j][i+1];
                uw = (i-1 == 0)          ? user->g_bdry(x-hx,y,0.0,user)
                                         : au[j][i-1];
                un = (j+1 == info->my-1) ? user->g_bdry(x,y+hy,0.0,user)
                                         : au[j+1][i];
                us = (j-1 == 0)          ? user->g_bdry(x,y-hy,0.0,user)
                                         : au[j-1][i];
                aF[j][i] = scdiag * au[j][i]
                           - scx * (uw + ue) - scy * (us + un)
                           - darea * user->f_rhs(x,y,0.0,user);
            }
        }
    }
    ierr = PetscLogFlops(11.0*info->xm*info->ym);CHKERRQ(ierr);
    return 0;
}
//ENDFORM2DFUNCTION

PetscErrorCode Poisson3DFunctionLocal(DMDALocalInfo *info, double ***au,
                                   double ***aF, PoissonCtx *user) {
    PetscErrorCode ierr;
    int    i, j, k;
    double xyzmin[3], xyzmax[3], hx, hy, hz, dvol, scx, scy, scz, scdiag,
           x, y, z, ue, uw, un, us, uu, ud;
    ierr = DMDAGetBoundingBox(info->da,xyzmin,xyzmax); CHKERRQ(ierr);
    hx = (xyzmax[0] - xyzmin[0]) / (info->mx - 1);
    hy = (xyzmax[1] - xyzmin[1]) / (info->my - 1);
    hz = (xyzmax[2] - xyzmin[2]) / (info->mz - 1);
    dvol = hx * hy * hz;
    scx = user->cx * dvol / (hx*hx);
    scy = user->cy * dvol / (hy*hy);
    scz = user->cz * dvol / (hz*hz);
    scdiag = 2.0 * (scx + scy + scz);
    for (k = info->zs; k < info->zs + info->zm; k++) {
        z = xyzmin[2] + k * hz;
        for (j = info->ys; j < info->ys + info->ym; j++) {
            y = xyzmin[1] + j * hy;
            for (i = info->xs; i < info->xs + info->xm; i++) {
                x = xyzmin[0] + i * hx;
                if (   i==0 || i==info->mx-1
                    || j==0 || j==info->my-1
                    || k==0 || k==info->mz-1) {
                    aF[k][j][i] = au[k][j][i] - user->g_bdry(x,y,z,user);
                    aF[k][j][i] *= scdiag;
                } else {
                    ue = (i+1 == info->mx-1) ? user->g_bdry(x+hx,y,z,user)
                                             : au[k][j][i+1];
                    uw = (i-1 == 0)          ? user->g_bdry(x-hx,y,z,user)
                                             : au[k][j][i-1];
                    un = (j+1 == info->my-1) ? user->g_bdry(x,y+hy,z,user)
                                             : au[k][j+1][i];
                    us = (j-1 == 0)          ? user->g_bdry(x,y-hy,z,user)
                                             : au[k][j-1][i];
                    uu = (k+1 == info->mz-1) ? user->g_bdry(x,y,z+hz,user)
                                             : au[k+1][j][i];
                    ud = (k-1 == 0)          ? user->g_bdry(x,y,z-hz,user)
                                             : au[k-1][j][i];
                    aF[k][j][i] = scdiag * au[k][j][i]
                        - scx * (uw + ue) - scy * (us + un) - scz * (uu + ud)
                        - dvol * user->f_rhs(x,y,z,user);
                }
            }
        }
    }
    ierr = PetscLogFlops(14.0*info->xm*info->ym*info->zm);CHKERRQ(ierr);
    return 0;
}

PetscErrorCode Poisson1DJacobianLocal(DMDALocalInfo *info, PetscScalar *au,
                                   Mat J, Mat Jpre, PoissonCtx *user) {
    PetscErrorCode  ierr;
    int          i,ncols;
    double       xmin[1], xmax[1], h, v[3];
    MatStencil   col[3],row;

    ierr = DMDAGetBoundingBox(info->da,xmin,xmax); CHKERRQ(ierr);
    h = (xmax[0] - xmin[0]) / (info->mx - 1);
    for (i = info->xs; i < info->xs+info->xm; i++) {
        row.i = i;
        col[0].i = i;
        ncols = 1;
        if (i==0 || i==info->mx-1) {
            v[0] = user->cx * 2.0 / h;
        } else {
            v[0] = user->cx * 2.0 / h;
            if (i-1 > 0) {
                col[ncols].i = i-1;  v[ncols++] = - user->cx / h;
            }
            if (i+1 < info->mx-1) {
                col[ncols].i = i+1;  v[ncols++] = - user->cx / h;
            }
        }
        ierr = MatSetValuesStencil(Jpre,1,&row,ncols,col,v,INSERT_VALUES); CHKERRQ(ierr);
    }

    ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != Jpre) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}

PetscErrorCode Poisson2DJacobianLocal(DMDALocalInfo *info, PetscScalar **au,
                                   Mat J, Mat Jpre, PoissonCtx *user) {
    PetscErrorCode  ierr;
    double       xymin[2], xymax[2], hx, hy, scx, scy, scdiag, v[5];
    int          i,j,ncols;
    MatStencil   col[5],row;

    ierr = DMDAGetBoundingBox(info->da,xymin,xymax); CHKERRQ(ierr);
    hx = (xymax[0] - xymin[0]) / (info->mx - 1);
    hy = (xymax[1] - xymin[1]) / (info->my - 1);
    scx = user->cx * hy / hx;
    scy = user->cy * hx / hy;
    scdiag = 2.0 * (scx + scy);
    for (j = info->ys; j < info->ys+info->ym; j++) {
        row.j = j;
        col[0].j = j;
        for (i = info->xs; i < info->xs+info->xm; i++) {
            row.i = i;
            col[0].i = i;
            ncols = 1;
            v[0] = scdiag;
            if (i>0 && i<info->mx-1 && j>0 && j<info->my-1) {
                if (i-1 > 0) {
                    col[ncols].j = j;    col[ncols].i = i-1;  v[ncols++] = - scx;  }
                if (i+1 < info->mx-1) {
                    col[ncols].j = j;    col[ncols].i = i+1;  v[ncols++] = - scx;  }
                if (j-1 > 0) {
                    col[ncols].j = j-1;  col[ncols].i = i;    v[ncols++] = - scy;  }
                if (j+1 < info->my-1) {
                    col[ncols].j = j+1;  col[ncols].i = i;    v[ncols++] = - scy;  }
            }
            ierr = MatSetValuesStencil(Jpre,1,&row,ncols,col,v,INSERT_VALUES); CHKERRQ(ierr);
        }
    }

    ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != Jpre) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}

PetscErrorCode Poisson3DJacobianLocal(DMDALocalInfo *info, PetscScalar ***au,
                                   Mat J, Mat Jpre, PoissonCtx *user) {
    PetscErrorCode  ierr;
    double       xyzmin[3], xyzmax[3], hx, hy, hz, dvol, scx, scy, scz, scdiag, v[7];
    int          i,j,k,ncols;
    MatStencil   col[7],row;

    ierr = DMDAGetBoundingBox(info->da,xyzmin,xyzmax); CHKERRQ(ierr);
    hx = (xyzmax[0] - xyzmin[0]) / (info->mx - 1);
    hy = (xyzmax[1] - xyzmin[1]) / (info->my - 1);
    hz = (xyzmax[2] - xyzmin[2]) / (info->mz - 1);
    dvol = hx * hy * hz;
    scx = user->cx * dvol / (hx*hx);
    scy = user->cy * dvol / (hy*hy);
    scz = user->cz * dvol / (hz*hz);
    scdiag = 2.0 * (scx + scy + scz);
    for (k = info->zs; k < info->zs+info->zm; k++) {
        row.k = k;
        col[0].k = k;
        for (j = info->ys; j < info->ys+info->ym; j++) {
            row.j = j;
            col[0].j = j;
            for (i = info->xs; i < info->xs+info->xm; i++) {
                row.i = i;
                col[0].i = i;
                ncols = 1;
                v[0] = scdiag;
                if (i>0 && i<info->mx-1 && j>0 && j<info->my-1 && k>0 && k<info->mz-1) {
                    if (i-1 > 0) {
                        col[ncols].k = k;    col[ncols].j = j;    col[ncols].i = i-1;
                        v[ncols++] = - scx;
                    }
                    if (i+1 < info->mx-1) {
                        col[ncols].k = k;    col[ncols].j = j;    col[ncols].i = i+1;
                        v[ncols++] = - scx;
                    }
                    if (j-1 > 0) {
                        col[ncols].k = k;    col[ncols].j = j-1;  col[ncols].i = i;
                        v[ncols++] = - scy;
                    }
                    if (j+1 < info->my-1) {
                        col[ncols].k = k;    col[ncols].j = j+1;  col[ncols].i = i;
                        v[ncols++] = - scy;
                    }
                    if (k-1 > 0) {
                        col[ncols].k = k-1;  col[ncols].j = j;    col[ncols].i = i;
                        v[ncols++] = - scz;
                    }
                    if (k+1 < info->mz-1) {
                        col[ncols].k = k+1;  col[ncols].j = j;    col[ncols].i = i;
                        v[ncols++] = - scz;
                    }
                }
                ierr = MatSetValuesStencil(Jpre,1,&row,ncols,col,v,INSERT_VALUES); CHKERRQ(ierr);
            }
        }
    }
    ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != Jpre) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}

PetscErrorCode InitialState(DM da, InitialType it, PetscBool gbdry,
                            Vec u, PoissonCtx *user) {
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    PetscRandom    rctx;
    switch (it) {
        case ZEROS:
            ierr = VecSet(u,0.0); CHKERRQ(ierr);
            break;
        case RANDOM:
            ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx); CHKERRQ(ierr);
            ierr = VecSetRandom(u,rctx); CHKERRQ(ierr);
            ierr = PetscRandomDestroy(&rctx); CHKERRQ(ierr);
            break;
        default:
            SETERRQ(PETSC_COMM_WORLD,4,"invalid InitialType ... how did I get here?\n");
    }
    if (!gbdry) {
        return 0;
    }
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    switch (info.dim) {
        case 1:
        {
            int    i;
            double xmax[1], xmin[1], h, x, *au;
            ierr = DMDAVecGetArray(da, u, &au); CHKERRQ(ierr);
            ierr = DMDAGetBoundingBox(da,xmin,xmax); CHKERRQ(ierr);
            h = (xmax[0] - xmin[0]) / (info.mx - 1);
            for (i = info.xs; i < info.xs + info.xm; i++) {
                if (i==0 || i==info.mx-1) {
                    x = xmin[0] + i * h;
                    au[i] = user->g_bdry(x,0.0,0.0,user);
                }
            }
            ierr = DMDAVecRestoreArray(da, u, &au); CHKERRQ(ierr);
            break;
        }
        case 2:
        {
            int     i, j;
            double  xymin[2], xymax[2], hx, hy, x, y, **au;
            ierr = DMDAVecGetArray(da, u, &au); CHKERRQ(ierr);
            ierr = DMDAGetBoundingBox(da,xymin,xymax); CHKERRQ(ierr);
            hx = (xymax[0] - xymin[0]) / (info.mx - 1);
            hy = (xymax[1] - xymin[1]) / (info.my - 1);
            for (j = info.ys; j < info.ys + info.ym; j++) {
                y = xymin[1] + j * hy;
                for (i = info.xs; i < info.xs + info.xm; i++) {
                    if (i==0 || i==info.mx-1 || j==0 || j==info.my-1) {
                        x = xymin[0] + i * hx;
                        au[j][i] = user->g_bdry(x,y,0.0,user);
                    }
                }
            }
            ierr = DMDAVecRestoreArray(da, u, &au); CHKERRQ(ierr);
            break;
        }
        case 3:
        {
            int     i, j, k;
            double  xyzmin[3], xyzmax[3], hx, hy, hz, x, y, z, ***au;
            ierr = DMDAVecGetArray(da, u, &au); CHKERRQ(ierr);
            ierr = DMDAGetBoundingBox(da,xyzmin,xyzmax); CHKERRQ(ierr);
            hx = (xyzmax[0] - xyzmin[0]) / (info.mx - 1);
            hy = (xyzmax[1] - xyzmin[1]) / (info.my - 1);
            hz = (xyzmax[2] - xyzmin[2]) / (info.mz - 1);
            for (k = info.zs; k < info.zs+info.zm; k++) {
                z = xyzmin[2] + k * hz;
                for (j = info.ys; j < info.ys + info.ym; j++) {
                    y = xyzmin[1] + j * hy;
                    for (i = info.xs; i < info.xs + info.xm; i++) {
                        if (i==0 || i==info.mx-1 || j==0 || j==info.my-1
                                 || k==0 || k==info.mz-1) {
                            x = xyzmin[0] + i * hx;
                            au[k][j][i] = user->g_bdry(x,y,z,user);
                        }
                    }
                }
            }
            ierr = DMDAVecRestoreArray(da, u, &au); CHKERRQ(ierr);
            break;
        }
        default:
            SETERRQ(PETSC_COMM_WORLD,5,"invalid dim from DMDALocalInfo\n");
    }
    return 0;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <inttypes.h>
#include <math.h>
#include <assert.h>
#include <nlopt.h>
// The "procedural" interface to Reduce/CSL.
#include "C_call_CSL.h"

// activate by setting to 1 or 2
int DEBUG = 0;

// The optimization "object":
static nlopt_opt opt       = NULL;
// The result computed by NLOPT_optimize:
static double *result;
// The local optimization sub-algorithm/object, for algorithms that use it:
static nlopt_opt local_opt = NULL;

// Nlopt return codes, see nlopt_result in nlopt.h.  Used only in few
// places. nlopt_main.red does error checking at the Reduce level.
char *nlopt_result_mess(nlopt_result r) {
  char *s;
  switch (r) {
  case -1: s = "nlopt failure!";  break;
  case -2: s = "invalid arguments or unsupported constraints (perhaps for this algorithm)!";  break;
  case -3: s = "out of memory!";  break;
  case -4: s = "progress limited by roundoff";  break;
  case -5: s = "forced stop";  break;
  case 1:  s = "nlopt success";  break;
  case 2:  s = "stop value reached";  break;
  case 3:  s=  "f-tolerance reached";  break;
  case 4:  s = "x-tolerance reached";  break;
  case 5:  s = "max evaluations reached";  break;
  case 6:  s = "max time reached";  break;
  };
  return s;
}

#define ERRORCHECK(who,r) if (r != NLOPT_SUCCESS) \
    {fprintf(stderr, "*%s: %s\n", who, nlopt_result_mess(r));}

// =======================================================================
// Basics
// =======================================================================

void NLOPT_create(int32_t anum, int32_t n) {
  opt = nlopt_create(anum, n);
  if (opt == NULL)
    fprintf(stderr, "Couldn't create '%s'!\n", nlopt_algorithm_name(anum));
    // nlopt_algorithm_name is built-in to the NLopt library
}

void NLOPT_destroy() {
  nlopt_destroy(opt);
}

int32_t NLOPT_get_algorithm() {
  // The only reason for this 'NLOPT_' is to convert from nlopt_algorithm to 'int32_t'!
  // Isn't there a better way?
  int32_t a = nlopt_get_algorithm(opt);
  return a;
}

int32_t NLOPT_get_dimension() {
  // The only reason for this 'NLOPT_' is to convert from "unsigned" to 'int32_t'!
  int32_t d = nlopt_get_dimension(opt);
  return d;
}

// =======================================================================
// Setting upper and lower bounds (hyper-rectangle)
// =======================================================================

int32_t NLOPT_set_lower_bounds(const double *lb) {
  int32_t r = nlopt_set_lower_bounds(opt, lb);
  return r;
}
int32_t NLOPT_set_upper_bounds(const double *ub) {
  int32_t r = nlopt_set_upper_bounds(opt, ub);
  return r;
}

// acts like a Icon generator
double NLOPT_get_lower_bound(int32_t i, int32_t n) {
  static double *slb = NULL;
  if (i == 0) {
    // Contortion 1:
    slb = (double *) malloc(n*sizeof(double));
    int32_t r = nlopt_get_lower_bounds(opt, slb);
    ERRORCHECK("NLOPT_get_lower_bound",r);
  }
  // Contortion 2:
  double lbi = slb[i];
  if (i == n-1) free(slb);
  return lbi;
}

// acts like a Icon generator
double NLOPT_get_upper_bound(int32_t i, int32_t n) {
  static double *sub = NULL;
  if (i == 0) {
    // Contortion 1:
    sub = (double *) malloc(n*sizeof(double));
    int32_t r = nlopt_get_upper_bounds(opt, sub);
    ERRORCHECK("NLOPT_get_upper_bound",r);
  }
  // Contortion 2:
  double ubi = sub[i];
  if (i == n-1) free(sub);
  return ubi;
}

// =======================================================================
// Stopping criteria
// =======================================================================

int32_t NLOPT_set_xtol_rel(double t) {
  return nlopt_set_xtol_rel(opt, t);
}
int32_t NLOPT_set_ftol_rel(double t) {
  return nlopt_set_ftol_rel(opt, t);
}
int32_t NLOPT_set_ftol_abs(double t) {
  return nlopt_set_ftol_abs(opt, t);
}
int32_t NLOPT_set_maxeval(int32_t m) {
  return nlopt_set_maxeval(opt, m);
}
int32_t NLOPT_set_maxtime(double m) {
  return nlopt_set_maxtime(opt, m);
}
int32_t NLOPT_set_stopval(double v) {
  return nlopt_set_stopval(opt, v);
}

double NLOPT_get_xtol_rel() {
  return nlopt_get_xtol_rel(opt);
}
double NLOPT_get_ftol_rel() {
  return nlopt_get_ftol_rel(opt);
}
double NLOPT_get_ftol_abs() {
  return nlopt_get_ftol_abs(opt);
}
int32_t NLOPT_get_maxeval() {
  return nlopt_get_maxeval(opt);
}
double NLOPT_get_maxtime() {
  return nlopt_get_maxtime(opt);
}

// ===========================================================================
// Functions *internal* to the interface for dealing with arrays of doubles
// To understand what's going on see the comments in nlopt_set_lower_bounds()
// in nlopt_main.red.
// ===========================================================================
double *newDoubleArray(int32_t len) {
  double *a;
  a = (double *) malloc(len*sizeof(double));
  return a;
}
void setDoubleArray(double a[], int32_t i, double x) {
  a[i] = x;
}
void freeDoubleArray(double a[]) {
  free(a);
}
// This internal function turns a C array of doubles into a Lisp list and
// leaves it at the top of the stack.  The "stack" is used by the "procedural"
// interface to Reduce, see file "C_call_CSL.h".
void C_darray_to_Lisp_list(int32_t n, const double x[]) {
  // Push the x[i] on the stack.
  for (int i = 0; i < n; i++) PROC_push_floating(x[i]);
  // Make a list out of the x[i]; it is left at the top of the stack.
  PROC_make_function_call("list", n);
  // Now having put 'list on the front, get rid of it.
  // (This perhaps shows that that ACN should really extend the procedural
  // interface to make this slightly nicer. But with luck this will suffice
  // for now.)
  PROC_lisp_eval();
}

// ================================================================
// Setting the optimization objective
// ================================================================
// This is the most complex and difficult part of the interface so far.

// This function is internal to the interface (static), and is of type nlopt_func.
// It is modelled after f77_func_wrap in api/f77api.c.  It is used as a wrapper
// for scalar functions, which can represent either the optimization objective,
// or a general equality or inequality constraint.
// The PROC_ code is borrowed from cuba_integrand in redcuba.c of the 'cuba' package.
static double red_func_wrap(unsigned n, const double *x, double *grad, void *f_data) {
  double z;
  char *fn = (char *) f_data;
  if (DEBUG > 1) {
    fprintf(stderr, "red_func_wrap: %s", fn);
    if (grad) fprintf(stderr, " [G]");
    fprintf(stderr, "\n");
  }
  // Convert the function name (string) to a symbol, push it on the stack, and
  // wrap it in a quote:
  PROC_push_symbol(fn);
  PROC_make_function_call("quote", 1);
  // Turn the array 'x' into a Lisp list and leave it at the top of the "stack".
  C_darray_to_Lisp_list(n,x);
  // Wrap it in a quote:
  PROC_make_function_call("quote", 1);
  // Evaluate a call to the Reduce algebraic function 'fn' that takes an n-list
  // of reals 'x', and delivers either a single real, f(x), or, if 'grad' is not
  // NULL, the flattened (n+1)-list of reals {f(x),grad(x)}.
  int32_t r = PROC_make_function_call("call_alg_fn_list2list", 2);
  PROC_lisp_eval();
  if (DEBUG > 1) {
    fprintf(stderr, "  %s(", fn);
    for (int i=0; i<n; i++) {
      fprintf(stderr, "%.4g", x[i]);
      if (i < n-1) fprintf(stderr, ",");
    }
  };
  if (grad) {
    // If 'grad' is not NULL, the Reduce function 'fn' is expected to return a
    // flattened list {f(x),grad(x)}, where grad(x) is an n-list which should
    // be assigned to 'grad'. (Also see the function red_mfunc_wrap.)
    PROC_handle y = PROC_get_value();
    z = PROC_floating_value(PROC_first(y));
    if (DEBUG > 1) {
      fprintf(stderr, ") = %.5g\n", z);
      fprintf(stderr, "  grad = {");
    };
    y = PROC_rest(y);
    for (int i = 0; i < n; i++) {
      grad[i] = PROC_floating_value(PROC_first(y));
      y = PROC_rest(y);
      if (DEBUG > 1) {
	fprintf(stderr, "%.4g", grad[i]);
	if (i < n-1) fprintf(stderr, ",");
      }
    }
    if (DEBUG > 1) fprintf(stderr, "}\n");
  }
  else {  // grad is NULL, result is just f(x)
    z = PROC_floating_value(PROC_get_value());
    if (DEBUG > 1) fprintf(stderr, ") = %.5g\n", z);
  }
  return z;
}

// Modelled after F77_nlo_set_min_objective in api/f77funcs_.h.
int32_t NLOPT_set_min_objective(char *fname) {
  // In the call to nlopt_set_min_objective:
  //   1. The wrapper red_func_wrap is the nlopt_func argument,
  //   2. The name of the Reduce function defining the objective is the *f_data argument,
  //   3. We can't just pass fname to to this argument, we need the string rfun on the heap!
  char *rfun = malloc(strlen(fname)*sizeof(char)+1);
  strcpy(rfun,fname);
  return nlopt_set_min_objective(opt, red_func_wrap, rfun);
}

int32_t NLOPT_set_max_objective(char *fname) {
  char *rfun = malloc(strlen(fname)*sizeof(char)+1);
  strcpy(rfun,fname);
  return nlopt_set_max_objective(opt, red_func_wrap, rfun);
}

int32_t NLOPT_optimize(double *x0, int32_t n) {
  double *x, f;
  x = (double *) malloc(n*sizeof(double));
  for (int j=0; j<n; j++) x[j] = x0[j];
  if (DEBUG) fprintf(stderr, "NLOPT_optimize:\n");
  int32_t r = nlopt_optimize(opt, x, &f);
  // assign results to the static 'result' array:
  result = (double *) malloc((n+1)*sizeof(double));
  result[0] = f;
  for (int j=0; j<n; j++) result[j+1] = x[j];
  free(x);
  if (DEBUG) {
    fprintf(stderr, "  f = %f\n", result[0]);
    fprintf(stderr, "  result =");
    for (int j=0; j<n; j++)
      fprintf(stderr, " %.5g", result[j+1]);
    fprintf(stderr, "\n");
  }
  return r;
}

double get_result(int32_t i) {
  return result[i];
}

void free_result() {
  free(result);
}

// ================================================================
// Setting the local optimization sub-algorithm
// ================================================================

int32_t NLOPT_set_local_optimizer(int32_t anum) {
  local_opt = nlopt_create(anum, nlopt_get_dimension(opt));
  return nlopt_set_local_optimizer(opt, local_opt);
}

// =============================================================================
// Defining scalar constraints
// =============================================================================

// Note: in the next two functions
// 1. The wrapper function red_func_wrap is the nlopt_func argument,
// 2. The name of the Reduce function defining the constraint is the *f_data argument,
// 3. We can't just pass c_fname to to this argument, we need the string cfun on the heap!

int32_t NLOPT_add_eq_constraint(char *c_fname, double tol) {
  char *cfun = malloc(strlen(c_fname)*sizeof(char)+1);
  assert(cfun != NULL);
  strcpy(cfun,c_fname);
  if (DEBUG) fprintf(stderr, "added '%s' equality constraint, %.3g\n", cfun, tol);
  return nlopt_add_equality_constraint(opt, red_func_wrap, cfun, tol);
}
int32_t NLOPT_add_ineq_constraint(char *c_fname, double tol) {
  char *cfun = malloc(strlen(c_fname)*sizeof(char)+1);
  assert(cfun != NULL);
  strcpy(cfun,c_fname);
  if (DEBUG) fprintf(stderr, "added '%s' inequality constraint, %.3g\n", cfun, tol);
  return nlopt_add_inequality_constraint(opt, red_func_wrap, cfun, tol);
}

// =============================================================================
// Defining vector-valued constraints
// =============================================================================

// This function is internal to the interface (static), and is of type nlopt_mfunc.
// It is modelled after f77_mfunc_wrap in api/f77api.c.  It is used as a wrapper for
// functions that, given an n-vector 'x', return an m-vector in 'result'.
static void red_mfunc_wrap(unsigned m, double *result, unsigned n, const double *x,
			   double *grad, void *f_data) {
  char *fn = (char *) f_data;
  if (DEBUG > 1) {
    fprintf(stderr, "red_mfunc_wrap: %s", fn);
    if (grad) fprintf(stderr, " [G]");
    fprintf(stderr, "\n");
  }
  // Convert the function name (string) to a symbol, push it on the stack, and
  // wrap it in a quote:
  PROC_push_symbol(fn);
  PROC_make_function_call("quote", 1);
  // Turn the array x into a Lisp list and leave it at the top of the "stack".
  C_darray_to_Lisp_list(n,x);
  // Wrap it in a quote:
  PROC_make_function_call("quote", 1);
  // Evaluate a call to the Reduce m-dimensional algebraic function fn: R^n -> R^m
  int32_t r = PROC_make_function_call("call_alg_fn_list2list", 2);
  PROC_lisp_eval();
  if (DEBUG > 1) {
    fprintf(stderr, "  %s(", fn);
    for (int i=0; i<n; i++) {
      fprintf(stderr, "%.4g", x[i]);
      if (i < n-1) fprintf(stderr, ",");
    }
    fprintf(stderr, ") = {");
  }
  // If 'grad' is NULL, the result is a an m-list, which we assign to the C
  // array 'result' element-by-element.  If 'grad' is non-NULL, then it points
  // to an array of length m*n, and the result of 'fn' also contains a list of
  // length m*n, the gradients of the constraint functions with respect to x,
  // and 'grad' should be set to this list.  It is assumed that ∂c_i/∂x_j is
  // stored in grad[i*n + j].
  // Irrespective of 'grad', we first get the first m elements of the list and
  // assign them to 'result'.
  PROC_handle y = PROC_get_value();
  for (int i = 0; i < m; i++) {
    result[i] = PROC_floating_value(PROC_first(y));
    y = PROC_rest(y);
    if (DEBUG > 1) {
      fprintf(stderr, "%.4g", result[i]);
      if (i < m-1)  fprintf(stderr, ",");
    }
  }
  if (DEBUG > 1) fprintf(stderr, "}\n");
  if (grad) {
    // Now handle the last m*n elements of what 'fn' has returned
    if (DEBUG > 1) fprintf(stderr, "  grad = {");
    for (int i = 0; i < m*n; i++) {
      grad[i] = PROC_floating_value(PROC_first(y));
      y = PROC_rest(y);
      if (DEBUG > 1) {
	fprintf(stderr, "%.4g", grad[i]);
	if (i < m*n-1) fprintf(stderr, ",");
      }
    }
    if (DEBUG > 1) fprintf(stderr, "}\n");
  }
}

// Note: look at the notes for the functions for scalar constraints
int32_t NLOPT_add_eq_mconstraint(char *c_fname, int32_t m, double *tol) {
  char *cfun = malloc(strlen(c_fname)*sizeof(char)+1);
  assert(cfun != NULL);
  strcpy(cfun,c_fname);
  if (DEBUG) fprintf(stderr, "added '%s' equality m-constraint\n", cfun);
  return nlopt_add_equality_mconstraint(opt, m, red_mfunc_wrap, cfun, tol);
}
int32_t NLOPT_add_ineq_mconstraint(char *c_fname, int32_t m, double *tol) {
  char *cfun = malloc(strlen(c_fname)*sizeof(char)+1);
  assert(cfun != NULL);
  strcpy(cfun,c_fname);
  if (DEBUG) fprintf(stderr, "added '%s' inequality m-constraint\n", cfun);
  return nlopt_add_inequality_mconstraint(opt, m, red_mfunc_wrap, cfun, tol);
}

// =============================================================================
// Removing constraints
// =============================================================================

int32_t NLOPT_remove_eq_constraints() {
  return nlopt_remove_equality_constraints(opt);
}
int32_t NLOPT_remove_ineq_constraints() {
  return nlopt_remove_inequality_constraints(opt);
}

// =============================================================================
// Step size for algorithms without derivatives
// =============================================================================

int32_t NLOPT_set_initial_step(const double *dx) {
  return nlopt_set_initial_step(opt, dx);
}
int32_t NLOPT_set_initial_step1(double dx) {
  return nlopt_set_initial_step1(opt, dx);
}

// ============================================================================
// Stochastic search algorithms
// ============================================================================

int32_t NLOPT_set_population(int p) {
  return nlopt_set_population(opt, (unsigned) p);
}

void NLOPT_srand(int64_t seed) {
  nlopt_srand((unsigned long) seed);
}

max_stars_repo_name: arthurcnorman/general

max_issues_repo_name: arthurcnorman/general

max_forks_repo_name: arthurcnorman/general


text:
#include "jfftw_complex_nd_Plan.h"
#include <fftw.h>

/*
 * Class:     jfftw_complex_nd_Plan
 * Method:    createPlan
 * Signature: ([III)V
 */
JNIEXPORT void JNICALL Java_jfftw_complex_nd_Plan_createPlan( JNIEnv *env, jobject obj, jintArray dim, jint dir, jint flags )
{
	jclass clazz;
	jfieldID id;
	jbyteArray arr;
	unsigned char* carr;
	int rank;
	int *cdim;

	if( sizeof( jdouble ) != sizeof( fftw_real ) )
	{
		(*env)->ThrowNew( env, (*env)->FindClass( env, "java/lang/RuntimeException" ), "jdouble and fftw_real are incompatible" );
		return;
	}

	clazz = (*env)->GetObjectClass( env, obj );
	id    = (*env)->GetFieldID( env, clazz, "plan", "[B" );
	arr   = (*env)->NewByteArray( env, sizeof( fftwnd_plan ) );
	carr  = (*env)->GetByteArrayElements( env, arr, 0 );
	rank  = (*env)->GetArrayLength( env, dim );
	cdim  = (*env)->GetIntArrayElements( env, dim, 0 );

	(*env)->MonitorEnter( env, (*env)->FindClass( env, "jfftw/Plan" ) );

	*(fftwnd_plan*)carr = fftwnd_create_plan( rank, cdim, dir, flags );

	(*env)->MonitorExit( env, (*env)->FindClass( env, "jfftw/Plan" ) );

	(*env)->ReleaseIntArrayElements( env, dim, cdim, 0 );
	(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
	(*env)->SetObjectField( env, obj, id, arr );
}
/*
 * Class:     jfftw_complex_nd_Plan
 * Method:    createPlanSpecific
 * Signature: ([III[DI[DI)V
 */
JNIEXPORT void JNICALL Java_jfftw_complex_nd_Plan_createPlanSpecific( JNIEnv *env, jobject obj, jintArray dim, jint dir, jint flags, jdoubleArray in, jint idist, jdoubleArray out, jint odist )
{
	jclass clazz;
	jfieldID id;
	jbyteArray arr;
	unsigned char* carr;
	int rank;
	int *cdim;
	double *cin, *cout;

	if( sizeof( jdouble ) != sizeof( fftw_real ) )
	{
		(*env)->ThrowNew( env, (*env)->FindClass( env, "java/lang/RuntimeException" ), "jdouble and fftw_real are incompatible" );
		return;
	}

	clazz = (*env)->GetObjectClass( env, obj );
	id    = (*env)->GetFieldID( env, clazz, "plan", "[B" );
	arr   = (*env)->NewByteArray( env, sizeof( fftwnd_plan ) );
	carr  = (*env)->GetByteArrayElements( env, arr, 0 );
	rank  = (*env)->GetArrayLength( env, dim );
	cdim  = (*env)->GetIntArrayElements( env, dim, 0 );
	cin   = (*env)->GetDoubleArrayElements( env, in, 0 );
	cout  = (*env)->GetDoubleArrayElements( env, out, 0 );

	(*env)->MonitorEnter( env, (*env)->FindClass( env, "jfftw/Plan" ) );

	*(fftwnd_plan*)carr = fftwnd_create_plan_specific( rank, cdim, dir, flags, (fftw_complex*)cin, idist, (fftw_complex*)cout, odist );

	(*env)->MonitorExit( env, (*env)->FindClass( env, "jfftw/Plan" ) );

	(*env)->ReleaseDoubleArrayElements( env, in, cin, 0 );
	(*env)->ReleaseDoubleArrayElements( env, out, cout, 0 );
	(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
	(*env)->SetObjectField( env, obj, id, arr );
}
/*
 * Class:     jfftw_complex_nd_Plan
 * Method:    destroyPlan
 * Signature: ()V
 */
JNIEXPORT void JNICALL Java_jfftw_complex_nd_Plan_destroyPlan( JNIEnv* env, jobject obj )
{
	jclass clazz = (*env)->GetObjectClass( env, obj );
	jfieldID id = (*env)->GetFieldID( env, clazz, "plan", "[B" );
	jbyteArray arr = (jbyteArray)(*env)->GetObjectField( env, obj, id );
	unsigned char* carr = (*env)->GetByteArrayElements( env, arr, 0 );

	fftwnd_destroy_plan( *(fftwnd_plan*)carr );

	(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
	(*env)->SetObjectField( env, obj, id, NULL );
}
/*
 * Class:     jfftw_complex_nd_Plan
 * Method:    transform
 * Signature: ([D)[D
 */
JNIEXPORT jdoubleArray JNICALL Java_jfftw_complex_nd_Plan_transform___3D( JNIEnv* env, jobject obj, jdoubleArray in )
{
	jdouble *cin, *cout;
	jdoubleArray out;
	int i;

	jclass clazz = (*env)->GetObjectClass( env, obj );
	jfieldID id = (*env)->GetFieldID( env, clazz, "plan", "[B" );
	jbyteArray arr = (jbyteArray)(*env)->GetObjectField( env, obj, id );
	unsigned char* carr = (*env)->GetByteArrayElements( env, arr, 0 );
	fftwnd_plan plan = *(fftwnd_plan*)carr;
	int length = 1;
	for( i = 0; i < plan->rank; ++i ) length *= plan->n[i];
	if( length * 2 != (*env)->GetArrayLength( env, in ) )
	{
		(*env)->ThrowNew( env, (*env)->FindClass( env, "java/lang/IndexOutOfBoundsException" ), "the Plan was created for a different length" );
		(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
		return NULL;
	}

	cin = (*env)->GetDoubleArrayElements( env, in, 0 );

	if( plan->rank > 0 && ! plan->plans[0]->flags & FFTW_THREADSAFE )
	{
		// synchronization
		(*env)->MonitorEnter( env, obj );
	}

	if( plan->is_in_place )
	{
		out = in;

		fftwnd_one( plan, (fftw_complex*)cin, NULL );
	}
	else
	{
		out = (*env)->NewDoubleArray( env, length * 2 );
		cout = (*env)->GetDoubleArrayElements( env, out, 0 );

		fftwnd_one( plan, (fftw_complex*)cin, (fftw_complex*)cout );

		(*env)->ReleaseDoubleArrayElements( env, out, cout, 0 );
	}

	if( plan->rank > 0 && ! plan->plans[0]->flags & FFTW_THREADSAFE )
	{
		// synchronization
		(*env)->MonitorExit( env, obj );
	}

	(*env)->ReleaseDoubleArrayElements( env, in, cin, 0 );
	(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
	return out;
}
/*
 * Class:     jfftw_complex_nd_Plan
 * Method:    transform
 * Signature: (I[DII[DII)V
 */
JNIEXPORT void JNICALL Java_jfftw_complex_nd_Plan_transform__I_3DII_3DII( JNIEnv *env, jobject obj, jint howmany, jdoubleArray in, jint istride, jint idist, jdoubleArray out, jint ostride, jint odist )
{
	jdouble *cin, *cout;
	int i;

	jclass clazz = (*env)->GetObjectClass( env, obj );
	jfieldID id = (*env)->GetFieldID( env, clazz, "plan", "[B" );
	jbyteArray arr = (jbyteArray)(*env)->GetObjectField( env, obj, id );
	unsigned char* carr = (*env)->GetByteArrayElements( env, arr, 0 );
	fftwnd_plan plan = *(fftwnd_plan*)carr;
	int length = 1;
	for( i = 0; i < plan->rank; ++i ) length *= plan->plans[i]->n;
	if( (howmany - 1) * idist * 2 + length * istride * 2 != (*env)->GetArrayLength( env, in ) )
	{
		(*env)->ThrowNew( env, (*env)->FindClass( env, "java/lang/IndexOutOfBoundsException" ), "the Plan was created for a different length (in)" );
		(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
		return;
	}
	if( (howmany - 1) * odist * 2 + length * ostride * 2 != (*env)->GetArrayLength( env, out ) )
	{
		(*env)->ThrowNew( env, (*env)->FindClass( env, "java/lang/IndexOutOfBoundsException" ), "the Plan was created for a different length (out)" );
		(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
		return;
	}

	cin = (*env)->GetDoubleArrayElements( env, in, 0 );
	cout = (*env)->GetDoubleArrayElements( env, out, 0 );

	if( plan->rank > 0 && ! plan->plans[0]->flags & FFTW_THREADSAFE )
	{
		// synchronization
		(*env)->MonitorEnter( env, obj );
	}

	fftwnd( plan, howmany, (fftw_complex*)cin, istride, idist, (fftw_complex*)cout, ostride, odist );

	if( plan->rank > 0 && ! plan->plans[0]->flags & FFTW_THREADSAFE )
	{
		// synchronization
		(*env)->MonitorExit( env, obj );
	}

	(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
	(*env)->ReleaseDoubleArrayElements( env, in, cin, 0 );
	(*env)->ReleaseDoubleArrayElements( env, out, cout, 0 );
}

max_stars_repo_name: CSCSI/Triana

max_issues_repo_name: CSCSI/Triana

max_forks_repo_name: CSCSI/Triana


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

#include <gsl/gsl_vector.h>

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

    gsl_vector * v1 = gsl_vector_calloc(10);
    gsl_vector * v2 = gsl_vector_calloc(10);

    gsl_vector_set_all(v1, 2.0);
    gsl_vector_set_all(v2, 6.0);

    gsl_vector_mul(v1, v2);

    gsl_vector_fprintf(stdout, v1, "%.2f");

    gsl_vector_free(v1);
    gsl_vector_free(v2);

    return 0;
}

max_stars_repo_name: jeremiedecock/snippets

max_issues_repo_name: jeremiedecock/snippets

max_forks_repo_name: jeremiedecock/snippets


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

#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])

int
main (void)
{
  int i;
  const int n = 630;
  double data[2*n];

  gsl_fft_complex_wavetable * wavetable;
  gsl_fft_complex_workspace * workspace;

  for (i = 0; i < n; i++)
    {
      REAL(data,i) = 0.0;
      IMAG(data,i) = 0.0;
    }

  data[0] = 1.0;

  for (i = 1; i <= 10; i++)
    {
      REAL(data,i) = REAL(data,n-i) = 1.0;
    }

  for (i = 0; i < n; i++)
    {
      printf ("%d: %e %e\n", i, REAL(data,i),
                                IMAG(data,i));
    }
  printf ("\n");

  wavetable = gsl_fft_complex_wavetable_alloc (n);
  workspace = gsl_fft_complex_workspace_alloc (n);

  for (i = 0; i < (int) wavetable->nf; i++)
    {
       printf ("# factor %d: %zu\n", i,
               wavetable->factor[i]);
    }

  gsl_fft_complex_forward (data, 1, n,
                           wavetable, workspace);

  for (i = 0; i < n; i++)
    {
      printf ("%d: %e %e\n", i, REAL(data,i),
                                IMAG(data,i));
    }

  gsl_fft_complex_wavetable_free (wavetable);
  gsl_fft_complex_workspace_free (workspace);
  return 0;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

/* Author:  G. Jungman */

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

#include <gsl/gsl_linalg.h>

#include "givens.c"
#include "apply_givens.c"

/* Factorise a general M x N matrix A into
 *
 *   A = Q R
 *
 * where Q is orthogonal (M x M) and R is upper triangular (M x N).
 *
 * Q is stored as a packed set of Householder transformations in the
 * strict lower triangular part of the input matrix.
 *
 * R is stored in the diagonal and upper triangle of the input matrix.
 *
 * The full matrix for Q can be obtained as the product
 *
 *       Q = Q_k .. Q_2 Q_1
 *
 * where k = MIN(M,N) and
 *
 *       Q_i = (I - tau_i * v_i * v_i')
 *
 * and where v_i is a Householder vector
 *
 *       v_i = [1, m(i+1,i), m(i+2,i), ... , m(M,i)]
 *
 * This storage scheme is the same as in LAPACK.  */

int
gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (tau->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      for (i = 0; i < GSL_MIN (M, N); i++)
        {
          /* Compute the Householder transformation to reduce the j-th
             column of the matrix to a multiple of the j-th unit vector */

          gsl_vector_view c_full = gsl_matrix_column (A, i);
          gsl_vector_view c = gsl_vector_subvector (&(c_full.vector), i, M-i);

          double tau_i = gsl_linalg_householder_transform (&(c.vector));

          gsl_vector_set (tau, i, tau_i);

          /* Apply the transformation to the remaining columns and
             update the norms */

          if (i + 1 < N)
            {
              gsl_matrix_view m = gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i + 1));
              gsl_linalg_householder_hm (tau_i, &(c.vector), &(m.matrix));
            }
        }

      return GSL_SUCCESS;
    }
}

/* Solves the system A x = b using the QR factorisation,

 *  R x = Q^T b
 *
 * to obtain x. Based on SLATEC code.
 */

int
gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x)
{
  if (QR->size1 != QR->size2)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (QR->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (QR->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Copy x <- b */

      gsl_vector_memcpy (x, b);

      /* Solve for x */

      gsl_linalg_QR_svx (QR, tau, x);

      return GSL_SUCCESS;
    }
}

/* Solves the system A x = b in place using the QR factorisation,

 *  R x = Q^T b
 *
 * to obtain x. Based on SLATEC code.
 */

int
gsl_linalg_QR_svx (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * x)
{

  if (QR->size1 != QR->size2)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (QR->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match x/rhs size", GSL_EBADLEN);
    }
  else
    {
      /* compute rhs = Q^T b */

      gsl_linalg_QR_QTvec (QR, tau, x);

      /* Solve R x = rhs, storing x in-place */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);

      return GSL_SUCCESS;
    }
}

/* Find the least squares solution to the overdetermined system
 *
 *   A x = b
 *
 * for M >= N using the QR factorization A = Q R.
 */

int
gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (M < N)
    {
      GSL_ERROR ("QR matrix must have M>=N", GSL_EBADLEN);
    }
  else if (M != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (N != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else if (M != residual->size)
    {
      GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_const_view R = gsl_matrix_const_submatrix (QR, 0, 0, N, N);
      gsl_vector_view c = gsl_vector_subvector(residual, 0, N);

      gsl_vector_memcpy(residual, b);

      /* compute rhs = Q^T b */

      gsl_linalg_QR_QTvec (QR, tau, residual);

      /* Solve R x = rhs */

      gsl_vector_memcpy(x, &(c.vector));

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, &(R.matrix), x);

      /* Compute residual = b - A x = Q (Q^T b - R x) */

      gsl_vector_set_zero(&(c.vector));

      gsl_linalg_QR_Qvec(QR, tau, residual);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const gsl_vector * b, gsl_vector * x)
{
  if (QR->size1 != QR->size2)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (QR->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (QR->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
    }
  else
    {
      /* Copy x <- b */

      gsl_vector_memcpy (x, b);

      /* Solve R x = b, storing x in-place */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector * x)
{
  if (QR->size1 != QR->size2)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (QR->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match rhs size", GSL_EBADLEN);
    }
  else
    {
      /* Solve R x = b, storing x in-place */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_R_solve (const gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
{
  if (R->size1 != R->size2)
    {
      GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
    }
  else if (R->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (R->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Copy x <- b */

      gsl_vector_memcpy (x, b);

      /* Solve R x = b, storing x inplace in b */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector * x)
{
  if (R->size1 != R->size2)
    {
      GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
    }
  else if (R->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Solve R x = b, storing x inplace in b */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);

      return GSL_SUCCESS;
    }
}

/* Form the product Q^T v  from a QR factorized matrix
 */

int
gsl_linalg_QR_QTvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (tau->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (v->size != M)
    {
      GSL_ERROR ("vector size must be N", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* compute Q^T v */

      for (i = 0; i < GSL_MIN (M, N); i++)
        {
          gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
          gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector), i, M - i);
          gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
          double ti = gsl_vector_get (tau, i);
          gsl_linalg_householder_hv (ti, &(h.vector), &(w.vector));
        }
      return GSL_SUCCESS;
    }
}

int
gsl_linalg_QR_Qvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (tau->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (v->size != M)
    {
      GSL_ERROR ("vector size must be N", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* compute Q^T v */

      for (i = GSL_MIN (M, N); i-- > 0;)
        {
          gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
          gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector),
                                                                i, M - i);
          gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
          double ti = gsl_vector_get (tau, i);
          gsl_linalg_householder_hv (ti, &h.vector, &w.vector);
        }
      return GSL_SUCCESS;
    }
}

/* Form the product Q^T A  from a QR factorized matrix */

int
gsl_linalg_QR_QTmat (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (tau->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (A->size1 != M)
    {
      GSL_ERROR ("matrix must have M rows", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* compute Q^T A */

      for (i = 0; i < GSL_MIN (M, N); i++)
        {
          gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
          gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector), i, M - i);
          gsl_matrix_view m = gsl_matrix_submatrix(A, i, 0, M - i, A->size2);
          double ti = gsl_vector_get (tau, i);
          gsl_linalg_householder_hm (ti, &(h.vector), &(m.matrix));
        }
      return GSL_SUCCESS;
    }
}

/*  Form the orthogonal matrix Q from the packed QR matrix */

int
gsl_linalg_QR_unpack (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (Q->size1 != M || Q->size2 != M)
    {
      GSL_ERROR ("Q matrix must be M x M", GSL_ENOTSQR);
    }
  else if (R->size1 != M || R->size2 != N)
    {
      GSL_ERROR ("R matrix must be M x N", GSL_ENOTSQR);
    }
  else if (tau->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else
    {
      size_t i, j;

      /* Initialize Q to the identity */

      gsl_matrix_set_identity (Q);

      for (i = GSL_MIN (M, N); i-- > 0;)
        {
          gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
          gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector,
                                                                i, M - i);
          gsl_matrix_view m = gsl_matrix_submatrix (Q, i, i, M - i, M - i);
          double ti = gsl_vector_get (tau, i);
          gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
        }

      /*  Form the right triangular matrix R from a packed QR matrix */

      for (i = 0; i < M; i++)
        {
          for (j = 0; j < i && j < N; j++)
            gsl_matrix_set (R, i, j, 0.0);

          for (j = i; j < N; j++)
            gsl_matrix_set (R, i, j, gsl_matrix_get (QR, i, j));
        }

      return GSL_SUCCESS;
    }
}

/* Update a QR factorisation for A= Q R ,  A' = A + u v^T,

 * Q' R' = QR + u v^T
 *       = Q (R + Q^T u v^T)
 *       = Q (R + w v^T)
 *
 * where w = Q^T u.
 *
 * Algorithm from Golub and Van Loan, "Matrix Computations", Section
 * 12.5 (Updating Matrix Factorizations, Rank-One Changes)
 */

int
gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R,
                      gsl_vector * w, const gsl_vector * v)
{
  const size_t M = R->size1;
  const size_t N = R->size2;

  if (Q->size1 != M || Q->size2 != M)
    {
      GSL_ERROR ("Q matrix must be M x M if R is M x N", GSL_ENOTSQR);
    }
  else if (w->size != M)
    {
      GSL_ERROR ("w must be length M if R is M x N", GSL_EBADLEN);
    }
  else if (v->size != N)
    {
      GSL_ERROR ("v must be length N if R is M x N", GSL_EBADLEN);
    }
  else
    {
      size_t j, k;
      double w0;

      /* Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0)

         J_1^T .... J_(n-1)^T w = +/- |w| e_1

         simultaneously applied to R,  H = J_1^T ... J^T_(n-1) R
         so that H is upper Hessenberg.  (12.5.2) */

      for (k = M - 1; k > 0; k--)  /* loop from k = M-1 to 1 */
        {
          double c, s;
          double wk = gsl_vector_get (w, k);
          double wkm1 = gsl_vector_get (w, k - 1);

          create_givens (wkm1, wk, &c, &s);
          apply_givens_vec (w, k - 1, k, c, s);
          apply_givens_qr (M, N, Q, R, k - 1, k, c, s);
        }

      w0 = gsl_vector_get (w, 0);

      /* Add in w v^T  (Equation 12.5.3) */

      for (j = 0; j < N; j++)
        {
          double r0j = gsl_matrix_get (R, 0, j);
          double vj = gsl_vector_get (v, j);
          gsl_matrix_set (R, 0, j, r0j + w0 * vj);
        }

      /* Apply Givens transformations R' = G_(n-1)^T ... G_1^T H
         Equation 12.5.4 */

      for (k = 1; k < GSL_MIN(M,N+1); k++)
        {
          double c, s;
          double diag = gsl_matrix_get (R, k - 1, k - 1);
          double offdiag = gsl_matrix_get (R, k, k - 1);

          create_givens (diag, offdiag, &c, &s);
          apply_givens_qr (M, N, Q, R, k - 1, k, c, s);

          gsl_matrix_set (R, k, k - 1, 0.0);    /* exact zero of G^T */
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
{
  const size_t M = R->size1;
  const size_t N = R->size2;

  if (M != N)
    {
      return GSL_ENOTSQR;
    }
  else if (Q->size1 != M || b->size != M || x->size != M)
    {
      return GSL_EBADLEN;
    }
  else
    {
      /* compute sol = Q^T b */

      gsl_blas_dgemv (CblasTrans, 1.0, Q, b, 0.0, x);

      /* Solve R x = sol, storing x in-place */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);

      return GSL_SUCCESS;
    }
}


text:
#define _XOPEN_SOURCE 600

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <math.h>
#include <time.h>

void
setupInitialSeed(
		 int id,
		 gsl_rng *rgen)
{

  struct timespec spec;

  clock_gettime( CLOCK_REALTIME, &spec );

  gsl_rng_set( rgen, (unsigned long int) ( round( spec.tv_nsec / 1.0e6 ) + 20 * id ) );

}

void createInitialSample(
			 double *x,
			 long N,
			 gsl_rng *rgen,
			 int state_dim )
{

  for ( long i = 0; i < N; i++)
    for ( int d = 0; d < state_dim; d++ )
      x[ state_dim * i + d ] = gsl_ran_gaussian_ziggurat( rgen, 1.0 );

}

max_stars_repo_name: heinekmp/AIRPF

max_issues_repo_name: heinekmp/AIRPF

max_forks_repo_name: heinekmp/AIRPF


text:
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"
#include "cblas.h"

void
cblas_drotmg (double *d1, double *d2, double *b1, const double b2, double *P)
{
#define BASE double
#include "source_rotmg.h"
#undef BASE
}


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef SIM_PHYSICS_WIND_H_
#define SIM_PHYSICS_WIND_H_

#include <gsl/gsl_vector.h>

#include <string>

#include "common/c_math/mat3.h"
#include "common/c_math/vec3.h"
#include "common/macros.h"

// Functions related to wind shear and turbulence.
//
// Many of these functions define a "mean wind frame" denoted by
// the "mw" subscript.  This should not be confused with the
// wind frame defined in sim/physics/aero_frame.h, which related to
// the apparent wind at a body.
//
// The origin of the mean wind frame lies at ground level directly
// below the ground station origin.  The axes are defined by:
//   x: Points along the predominate wind vector.
//   y: Completes a right-handed frame.
//   z: Points up.
// These are chosen to be consistent with NREL's TurbSim software
// and various IEC standards.

namespace sim {

namespace physics {

namespace wind {

double CalcWindShear(double height_agl, double ref_height_agl,
                     double ref_wind_speed, double wind_shear_exponent);

class FrozenTurbulenceWindDatabase {
  friend class FrozenTurbulenceWindDatabaseTest;

 public:
  FrozenTurbulenceWindDatabase(double t0, double y0,
                               const std::string &filename);
  ~FrozenTurbulenceWindDatabase() {
    gsl_vector_free(u_);
    gsl_vector_free(v_);
    gsl_vector_free(w_);
  }
  void CalcWind(double t, const Vec3 &pos_mw, Vec3 *wind_mw) const;

 private:
  // Time offset [s] to apply.  At t = 0, positions whose x coordinate
  // (in the mw frame) is zero will fetch the database entry given at
  // time t0_.
  double t0_;
  // Y offset [m] to apply from the center of the wind database.
  double y0_;
  // Mean wind speed [m/s] that sets the rate at which the database
  // is "advected" along the mean wind coordinate x axis.
  double mean_wind_speed_;
  // Number [#] of t, y, and z grid points.
  int32_t num_t_, num_y_, num_z_;
  // Duration [s], width [m], and height [m] of the grid defining
  // the database.
  double duration_, width_, height_;
  // Vectors storing the components of the wind velocity [m/s].  Each
  // vector has num_t * num_y * num_z elements.  The u, v, and w
  // components point along the x, y, and z axes of the mean wind
  // coordinate system.
  gsl_vector *u_;
  gsl_vector *v_;
  gsl_vector *w_;

  DISALLOW_COPY_AND_ASSIGN(FrozenTurbulenceWindDatabase);
};

enum class IecTurbulenceCategory { kA, kB, kC };

enum class IecWindTurbineClass { kI, kII, kIII };

class IecWindModel {
 public:
  IecWindModel(double hub_height_agl__, double rotor_diameter,
               double hub_wind_speed, IecTurbulenceCategory turbulence_category,
               IecWindTurbineClass wind_turbine_class);

  double hub_height_agl() const { return hub_height_agl_; }
  double ntm_sigma_1() const { return ntm_sigma_1_; }
  double ewm_sigma_1() const { return ewm_sigma_1_; }
  double etm_sigma_1() const { return etm_sigma_1_; }

  // IEC 61400-1 Design situations.

  // IEC 61400-1 Sec. 6.3.1.2 -- Normal Wind Profile Model (NWP).
  double CalcNormalWindProfile(double height_agl) const;

  // IEC 61400-1 Sec. 6.3.2.1 -- Extreme wind speed model (EWM).
  double CalcExtremeWindSpeed1YearRecurrence(double height_agl,
                                             bool is_turbulent) const;
  double CalcExtremeWindSpeed50YearRecurrence(double height_agl,
                                              bool is_turbulent) const;

  // IEC 61400-1 Sec. 6.3.2.2 -- Extreme Operating Gust (EOG).
  double CalcExtremeOperatingGust(double t, double height_agl) const;

  // IEC 61400-1 Sec. 6.3.2.4 -- Extreme direction change (EDC).
  double CalcExtremeDirectionChange(double t) const;

  // IEC 61400-1 Sec. 6.3.2.5 -- Extreme coherent gust with direction
  // change (ECD).
  double CalcExtremeCoherentGustWithDirectionChange(
      double t, double height_agl, double *wind_direction) const;

  // IEC 61400-1 Sec. 6.3.2.6 -- Extreme wind shear (EWS).
  double CalcExtremeWindShearVertical(double t, double height_agl) const;
  double CalcExtremeWindShearHorizontal(double t, double height_agl,
                                        double y) const;

 private:
  // Hub height above ground level [m].
  const double hub_height_agl_;

  // Rotor diameter [m].
  const double rotor_diameter_;

  // Nominal wind speed [m/s] at hub height.
  const double hub_wind_speed_;

  // Reference wind speed [m/s], see IEC 61400-1 Sec. 6.2.
  const double v_ref_;

  // Expected value of the turbulence intensity [#] at 15 m/s.
  const double I_ref_;

  // Longitudinal turbulence scale parameter [#].
  const double lambda_1_;

  // Normal turbulence model (NTM) standard deviation [m/s].
  const double ntm_sigma_1_;

  // Extreme wind model (EWM) standard deviation [m/s].
  const double ewm_sigma_1_;

  // Extreme turbulence model (ETM) standard deviation [m/s].
  const double etm_sigma_1_;
};

}  // namespace wind

}  // namespace physics

}  // namespace sim

#endif  // SIM_PHYSICS_WIND_H_


text:
#pragma once
#include <array>
#include <variant>
#include <gsl/gsl-lite.hpp>

#include <clrie/method_info.h>
#include <clrie/instruction_factory.h>

#include <corhdr.h>

#include "cil.h"
#include "signature.h"

namespace appmap {
    template <typename T>
    constexpr COR_SIGNATURE type_signature() {
        if constexpr (std::is_same_v<T, void>) {
            return ELEMENT_TYPE_VOID;
        } else if constexpr (std::is_same_v<T, int32_t>) {
            return ELEMENT_TYPE_I4;
        } else if constexpr (std::is_same_v<T, int64_t>) {
            return ELEMENT_TYPE_I8;
        } else if constexpr (std::is_same_v<T, uint32_t>) {
            return ELEMENT_TYPE_U4;
        } else if constexpr (std::is_same_v<T, uint64_t>) {
            return ELEMENT_TYPE_U8;
        } else if constexpr (std::is_same_v<T, char *> || std::is_same_v<T, const char *>) {
            return ELEMENT_TYPE_STRING;
        } else if constexpr (std::is_pointer_v<T>) {
            static_assert(sizeof(void *) == 8);
            return ELEMENT_TYPE_I8;
        } else if constexpr (std::is_same_v<T, bool>) {
            return ELEMENT_TYPE_BOOLEAN;
        } else {
            static_assert(std::is_same_v<T, void>, "unhandled native type");
        }
    }

    template <typename F>
    struct func_traits;

    template <typename Ret, typename... Args>
    struct func_traits<Ret(*)(Args...)>
    {
        static constexpr std::size_t arity = sizeof...(Args);
        static constexpr std::array<COR_SIGNATURE, 3 + arity> signature = {
            IMAGE_CEE_UNMANAGED_CALLCONV_STDCALL,
            arity,
            type_signature<Ret>(),
            type_signature<Args>()...
        };
    };

    template <typename C, typename Ret, typename... Args>
    struct func_traits<Ret(C::*)(Args...)>
    {
        using instance_type = C;
        static constexpr std::size_t arity = sizeof...(Args);
        static constexpr std::array<COR_SIGNATURE, 4 + arity> signature = {
            IMAGE_CEE_UNMANAGED_CALLCONV_STDCALL,
            arity + 1,
            type_signature<Ret>,
            type_signature<C*>,
            type_signature<Args>...
        };
    };

    struct instrumentation : public clrie::instruction_factory {
        instrumentation(clrie::method_info method_info) :
            clrie::instruction_factory(method_info.instruction_factory()),
            method(method_info),
            module(method.module_info()),
            module_id(module.module_id()),
            metadata(module.meta_data_emit()),
            type_factory(module.create_type_factory()),
            locals(method.get(&IMethodInfo::GetLocalVariables))
        {}

        clrie::method_info method;
        clrie::module_info module;
        const ModuleID module_id;
        com::ptr<IMetaDataEmit> metadata;
        com::ptr<ITypeCreator> type_factory;
        com::ptr<ILocalVariableCollection> locals;

        static com::ptr<ISignatureBuilder> signature_builder;

        template <class F>
        instruction_sequence make_call(F f) const {
            return make_call_sig(reinterpret_cast<void *>(f), gsl::make_span(func_traits<F>::signature));
        }

        instruction_sequence create_call_to_string(const clrie::type &type) const noexcept;

        // Note capture_value takes a reference; in case of a composite type, it dereferences
        // it to a primitive that can be then passed onto the correct native function.
        // The argument is updated to reflect the resulting simple type.
        instruction_sequence capture_value(clrie::type &type) const noexcept;

        template <typename T>
        uint64_t add_local()
        {
            auto t = type_factory.get(&ITypeCreator::FromCorElement, static_cast<CorElementType>(type_signature<T>()));
            return locals.get(&ILocalVariableCollection::AddLocal, t);
        }

        // emit metadata and return reference tokens

        mdAssemblyRef assembly_reference(const char16_t *assembly);

        mdTypeRef type_reference(mdAssemblyRef assembly, const char16_t *type);
        mdTypeRef type_reference(const char16_t *assembly, const char16_t *type) {
            return type_reference(assembly_reference(assembly), type);
        }

        mdMemberRef member_reference(mdTypeRef type, const char16_t *member,
            gsl::span<const COR_SIGNATURE> signature);

        mdMemberRef member_reference(mdAssemblyRef assembly, const char16_t *type,
            const char16_t *member, gsl::span<const COR_SIGNATURE> signature) {
            return member_reference(type_reference(assembly, type), member, signature);
        }

        mdMemberRef member_reference(const char16_t *assembly, const char16_t *type,
            const char16_t *member, gsl::span<const COR_SIGNATURE> signature) {
            return member_reference(type_reference(assembly, type), member, signature);
        }

        mdTypeDef define_type(const char16_t *name);
        mdFieldDef define_field(mdTypeDef type, const char16_t *name, gsl::span<const COR_SIGNATURE> signature);

        mdMethodDef define_method(
            mdTypeDef type,
            const char16_t *name,
            gsl::span<const COR_SIGNATURE> signature,
            std::initializer_list<appmap::signature::type> locals,
            std::vector<cil::instruction> code
        );

        mdMethodDef define_method(mdTypeDef type, const char16_t *name, gsl::span<const COR_SIGNATURE> signature, std::vector<cil::instruction> code) {
            return define_method(type, name, signature, {}, std::move(code));
        }

        template <class Ret, class... Args>
        mdTypeSpec native_type(Ret (*)(Args...)) {
            constexpr auto sig = func_traits<Ret(*)(Args...)>::signature;
            return metadata.get(&IMetaDataEmit::GetTokenFromSig, sig.data(), sig.size());
        }

        mdTypeSpec type_token(gsl::span<const COR_SIGNATURE> sig) {
            return metadata.get(&IMetaDataEmit::GetTokenFromTypeSpec, sig.data(), sig.size());
        }

    protected:
        instruction_sequence make_call_sig(void *fn, gsl::span<const COR_SIGNATURE> signature) const;
    };
}

max_stars_repo_name: applandinc/appmap-dotnet

max_issues_repo_name: applandinc/appmap-dotnet

max_forks_repo_name: applandinc/appmap-dotnet


text:
/* BalSelSims.c

Simulation of balancing selection simulation with sampling (finite population effects)
Selection; reproduction; mutation based on deterministic recursions, then sampling
Repeats for 'Length' generations

Simulation uses routines found with the GNU Scientific Library (GSL)
(http://www.gnu.org/software/gsl/)
Since GSL is distributed under the GNU General Public License
(http://www.gnu.org/copyleft/gpl.html), you must download it
separately from this file.

This program can be compiled in e.g. GCC using a command like:
gcc BalSelSims -lm -lgsl -lgslcblas -I/usr/local/include -L/usr/local/lib BalSelSims.c

Then run by executing:
./BalSelSims N s rec sex self gc reps
Where:
- N is the population size
- s is the fitness disadvantage of homozygotes
- rec is recombination rate
- sex is rate of sex (a value between 0 = obligate asex, and 1 = obligate sex)
- self is selfing rate
- gc is gene conversion
- reps is how many times to introduce linked neutral allele

Note that haplotypes are defined as:
x1 = ab
x2 = Ab
x3 = aB
x4 = AB

Genotypes defined as:
g11 = g1 = ab/ab
g12 = g2 = Ab/ab
g13 = g3 = aB/ab
g14 = g4 = AB/ab
g22 = g5 = Ab/Ab
g23 = g6 = Ab/aB
g24 = g7 = Ab/AB
g33 = g8 = aB/aB
g34 = g9 = aB/AB
g44 = g10 = AB/AB

*/

/* Preprocessor statements */
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <stddef.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* Function prototypes */
void geninit(double *geninit);
void selection(double *geninit);
void reproduction(double *geninit);
void gconv(double *geninit);
void neutinit(double *geninit,const gsl_rng *r);
double ncheck(double *geninit);
double pcheck(double *geninit);

/* Global variable declaration */
unsigned int N = 0;		/* Pop size */
double s = 0;			/* Fitness disadvantage of homozyogtes */
double rec = 0;			/* Recombination rate */
double sex = 0;			/* Rate of sexual reproduction */
double self = 0;		/* Rate of self-fertilisation */
double gc = 0;			/* Rate of gene conversion */

/* Main program */
int main(int argc, char *argv[]){
	unsigned int g, i; 				/* Counters. Reps counter, geno counter */
	unsigned int reps;				/* Length of simulation (no. of introductions of neutral site) */
	double Bcheck = 0;				/* Frequency of B after each reproduction */
	double Acheck = 0;				/* Frequency of polymorphism */
	double Hsum = 0;				/* Summed heterozygosity over transit time of neutral allele */

	/* GSL random number definitions */
	const gsl_rng_type * T;
	gsl_rng * r;

	/* This reads in data from command line. */
	if(argc != 8){
		fprintf(stderr,"Invalid number of input values.\n");
		exit(1);
	}
	N = strtod(argv[1],NULL);
	s = strtod(argv[2],NULL);
	rec = strtod(argv[3],NULL);
	sex = strtod(argv[4],NULL);
	self = strtod(argv[5],NULL);
	gc = strtod(argv[6],NULL);
	reps = strtod(argv[7],NULL);

	/* Arrays definition and memory assignment */
	double *genotype = calloc(10,sizeof(double));				/* Genotype frequencies */
	unsigned int *gensamp = calloc(10,sizeof(unsigned int));	/* New population samples */

	/* create a generator chosen by the
    environment variable GSL_RNG_TYPE */

	gsl_rng_env_setup();
	if (!getenv("GSL_RNG_SEED")) gsl_rng_default_seed = time(0);
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);

	/* Initialising genotypes */
	geninit(genotype);

    /* Run simulation for 2000 generations to create a burn in */
    for(g = 0; g < 2000; g++){

    	/* Selection routine */
    	selection(genotype);

    	/* Reproduction routine */
    	reproduction(genotype);

       	/* Gene conversion routine */
       	gconv(genotype);

       	/* Sampling based on new frequencies */
       	gsl_ran_multinomial(r,10,N,genotype,gensamp);
       	for(i = 0; i < 10; i++){
			*(genotype + i) = (*(gensamp + i))/(1.0*N);
       	}

       	/* Printing out results (for testing) */
		/*
	   	for(i = 0; i < 10; i++){
			printf("%.10lf ", *(genotype + i));
		}
		printf("\n");
		*/
    }

	/* Reintroducing neutral genotype, resetting hap sum */
	neutinit(genotype,r);
    Bcheck = ncheck(genotype);
    Hsum = Bcheck*(1-Bcheck);
    /* printf("%.10lf %.10lf\n",Bcheck,Hsum); */

    /* Introduce and track neutral mutations 'reps' times */
    g = 0;
    while(g < reps){

    	/* Selection routine */
    	selection(genotype);

    	/* Reproduction routine */
    	reproduction(genotype);

    	/* Gene conversion routine */
       	gconv(genotype);

    	/* Sampling based on new frequencies */
       	gsl_ran_multinomial(r,10,N,genotype,gensamp);
       	for(i = 0; i < 10; i++){
       		*(genotype + i) = (*(gensamp + i))/(1.0*N);
       	}

       	/* Checking state of haplotypes: if B fixed reset so can start fresh next time */
		Bcheck = ncheck(genotype);
		Hsum += Bcheck*(1-Bcheck);
		/* printf("%.10lf %.10lf\n",Bcheck,Hsum); */

		/* If polymorphism fixed then abandon simulation */
       	Acheck = pcheck(genotype);
       	if(Acheck == 0){
       		g = reps;
       	}

       	if(Bcheck == 0 || Bcheck == 1){
       		printf("%.10lf\n",Hsum);
       		g++;
       		/* printf("Rep Number %d\n",g); */

       		if(Bcheck == 1){
       			/* Reset genotypes so B becomes ancestral allele */
       			*(genotype + 0) = *(genotype + 7);
       			*(genotype + 1) = *(genotype + 8);
       			*(genotype + 4) = *(genotype + 9);
       			*(genotype + 7) = 0;
       			*(genotype + 8) = 0;
       			*(genotype + 9) = 0;
       		}

     		/* Reintroducing neutral genotype, resetting hap sum */
	    	neutinit(genotype,r);
			Bcheck = ncheck(genotype);
       		Hsum = Bcheck*(1-Bcheck);
       	}

	}	/* End of simulation */

	/* Freeing memory and wrapping up */
 	gsl_rng_free(r);
 	free(gensamp);
	free(genotype);
	/* printf("The End!\n"); */
	return 0;
}

/* Initialising genotypes */
void geninit(double *geninit){

	/* Basic idea: before neutral mutation introduced, g0 = 1/4; g1 = 1/2; g4 = 1/4. Deterministic frequencies. */

	*(geninit + 0) = 0.25;
	*(geninit + 1) = 0.50;
	*(geninit + 2) = 0;
	*(geninit + 3) = 0;
	*(geninit + 4) = 0.25;
	*(geninit + 5) = 0;
	*(geninit + 6) = 0;
	*(geninit + 7) = 0;
	*(geninit + 8) = 0;
	*(geninit + 9) = 0;

}	/* End of gen initiation routine */

/* Initialising NEUTRAL allele */
void neutinit(double *geninit, const gsl_rng *r){

	double *probin = calloc(4,sizeof(double));					/* Probability inputs, determine location of neutral allele */
	unsigned int *probout = calloc(4,sizeof(unsigned int));		/* Output from multinomial sampling */

	/* Basic idea: g11 at freq p^2; g12 at freq 2pq; g22 at freq q^2.
	So weighted sampling to determine which genotype the neutral allele arises on
	*/

	/* Prob definitions */
	*(probin + 0) = *(geninit + 0);
	*(probin + 1) = (*(geninit + 1))/2.0;
	*(probin + 2) = (*(geninit + 1))/2.0;
	*(probin + 3) = *(geninit + 4);

	gsl_ran_multinomial(r,4,1,probin,probout);

	/* Redefining genotypes depending on outcome */
	if(*(probout + 0) == 1){
		*(geninit + 0) = (*(geninit + 0)) - 1/(1.0*N);
		*(geninit + 2) = (*(geninit + 2)) + 1/(1.0*N);
	}
	else if(*(probout + 1) == 1){
		*(geninit + 1) = (*(geninit + 1)) - 1/(1.0*N);
		*(geninit + 3) = (*(geninit + 3)) + 1/(1.0*N);
	}
	else if(*(probout + 2) == 1){
		*(geninit + 1) = (*(geninit + 1)) - 1/(1.0*N);
		*(geninit + 5) = (*(geninit + 5)) + 1/(1.0*N);
	}
	else if(*(probout + 3) == 1){
		*(geninit + 4) = (*(geninit + 4)) - 1/(1.0*N);
		*(geninit + 6) = (*(geninit + 6)) + 1/(1.0*N);
	}

 	free(probout);
	free(probin);

}	/* End of gen initiation routine */

/* Selection routine */
void selection(double *geninit){
	double Waa, WAa, WAA;		/* Fitness of locus A (bal sel locus) */
	double Wmean;				/* Mean fitness */

	Waa = 1-s;
	WAa = 1;
	WAA = 1-s;

	/* Mean fitness calculation */
	Wmean = ((*(geninit + 0))*Waa) + ((*(geninit + 1))*WAa) + ((*(geninit + 2))*Waa) + ((*(geninit + 3))*WAa) + ((*(geninit + 4))*WAA) + ((*(geninit + 5))*WAa) + ((*(geninit + 6))*WAA) + ((*(geninit + 7))*Waa) + ((*(geninit + 8))*WAa) + ((*(geninit + 9))*WAA);

	/* Changing frequencies by selection */
	*(geninit + 0) = ((*(geninit + 0))*Waa)/Wmean;
	*(geninit + 1) = ((*(geninit + 1))*WAa)/Wmean;
	*(geninit + 2) = ((*(geninit + 2))*Waa)/Wmean;
	*(geninit + 3) = ((*(geninit + 3))*WAa)/Wmean;
	*(geninit + 4) = ((*(geninit + 4))*WAA)/Wmean;
	*(geninit + 5) = ((*(geninit + 5))*WAa)/Wmean;
	*(geninit + 6) = ((*(geninit + 6))*WAA)/Wmean;
	*(geninit + 7) = ((*(geninit + 7))*Waa)/Wmean;
	*(geninit + 8) = ((*(geninit + 8))*WAa)/Wmean;
	*(geninit + 9) = ((*(geninit + 9))*WAA)/Wmean;

}	/* End of selection routine */

/* Reproduction routine */
void reproduction(double *geninit){
	/* Fed-in genotype frequencies (for ease of programming) */
	double g11s, g12s, g13s, g14s, g22s, g23s, g24s, g33s, g34s, g44s;
	/* Genotype frequencies after sex (outcross and selfing) */
	double g11SX, g12SX, g13SX, g14SX, g22SX, g23SX, g24SX, g33SX, g34SX, g44SX;
	/* Genotype frequencies after ASEX */
	double g11AS, g12AS, g13AS, g14AS, g22AS, g23AS, g24AS, g33AS, g34AS, g44AS;
	/* Haplotypes */
	double x1, x2, x3, x4;

	/* Initial definition of genotypes */
	g11s = *(geninit + 0);
	g12s = *(geninit + 1);
	g13s = *(geninit + 2);
	g14s = *(geninit + 3);
	g22s = *(geninit + 4);
	g23s = *(geninit + 5);
	g24s = *(geninit + 6);
	g33s = *(geninit + 7);
	g34s = *(geninit + 8);
	g44s = *(geninit + 9);

	/* Baseline change in haplotype frequencies */
	x1 = g11s + (g12s + g13s + g14s)/2.0 - ((g14s - g23s)*rec)/2.0;
	x2 = g22s + (g12s + g23s + g24s)/2.0 + ((g14s - g23s)*rec)/2.0;
	x3 = g33s + (g13s + g23s + g34s)/2.0 + ((g14s - g23s)*rec)/2.0;
	x4 = g44s + (g14s + g24s + g34s)/2.0 - ((g14s - g23s)*rec)/2.0;

	/* Change in SEXUAL frequencies (both outcrossing and selfing) */
	g11SX = (g11s + (g12s + g13s + g14s*pow((1 - rec),2) + g23s*pow(rec,2))/4.0)*self*sex + (1 - self)*pow(x1,2)*sex;
	g22SX = (g22s + (g12s + g24s + g23s*pow((1 - rec),2) + g14s*pow(rec,2))/4.0)*self*sex + (1 - self)*pow(x2,2)*sex;
	g33SX = (g33s + (g13s + g34s + g23s*pow((1 - rec),2) + g14s*pow(rec,2))/4.0)*self*sex + (1 - self)*pow(x3,2)*sex;
	g44SX = (g44s + (g24s + g34s + g14s*pow((1 - rec),2) + g23s*pow(rec,2))/4.0)*self*sex + (1 - self)*pow(x4,2)*sex;
	g12SX = ((g12s + (g14s + g23s)*(1 - rec)*rec)*self*sex)/2.0 + 2.0*(1 - self)*x1*x2*sex;
	g13SX = ((g13s + (g14s + g23s)*(1 - rec)*rec)*self*sex)/2.0 + 2.0*(1 - self)*x1*x3*sex;
	g14SX = ((g14s*pow((1 - rec),2) + g23s*pow(rec,2))*self*sex)/2.0 + 2.0*(1 - self)*x1*x4*sex;
	g23SX = ((g23s*pow((1 - rec),2) + g14s*pow(rec,2))*self*sex)/2.0 + 2.0*(1 - self)*x2*x3*sex;
	g24SX = ((g24s + (g14s + g23s)*(1 - rec)*rec)*self*sex)/2.0 + 2.0*(1 - self)*x2*x4*sex;
	g34SX = ((g34s + (g14s + g23s)*(1 - rec)*rec)*self*sex)/2.0 + 2.0*(1 - self)*x3*x4*sex;

	/* Change in ASEXUAL frequencies */
	g11AS = g11s*(1 - sex);
	g12AS = g12s*(1 - sex);
	g13AS = g13s*(1 - sex);
	g14AS = g14s*(1 - sex);
	g22AS = g22s*(1 - sex);
	g23AS = g23s*(1 - sex);
	g24AS = g24s*(1 - sex);
	g33AS = g33s*(1 - sex);
	g34AS = g34s*(1 - sex);
	g44AS = g44s*(1 - sex);

	/* Combining to give overall frequency change following reproduction */
	*(geninit + 0) = g11AS + g11SX;
	*(geninit + 1) = g12AS + g12SX;
	*(geninit + 2) = g13AS + g13SX;
	*(geninit + 3) = g14AS + g14SX;
	*(geninit + 4) = g22AS + g22SX;
	*(geninit + 5) = g23AS + g23SX;
	*(geninit + 6) = g24AS + g24SX;
	*(geninit + 7) = g33AS + g33SX;
	*(geninit + 8) = g34AS + g34SX;
	*(geninit + 9) = g44AS + g44SX;

}	/* End of reproduction routine */

/* Gene conversion routine */
void gconv(double *geninit){

	/* Fed-in genotype frequencies (for ease of programming) */
	double g11r, g12r, g13r, g14r, g22r, g23r, g24r, g33r, g34r, g44r;
	/* Frequencies after gene conversion */
	double g11gc, g12gc, g13gc, g14gc, g22gc, g23gc, g24gc, g33gc, g34gc, g44gc;

	/* Initial definition of genotypes */
	g11r = *(geninit + 0);
	g12r = *(geninit + 1);
	g13r = *(geninit + 2);
	g14r = *(geninit + 3);
	g22r = *(geninit + 4);
	g23r = *(geninit + 5);
	g24r = *(geninit + 6);
	g33r = *(geninit + 7);
	g34r = *(geninit + 8);
	g44r = *(geninit + 9);

	/* Gene conversion equations */
	g11gc = g11r + (gc*g12r)/4.0 + (gc*g13r)/4.0;
	g12gc = g12r*(1 - gc/2.0) + ((g14r + g23r)*gc)/4.0;
	g13gc = g13r*(1 - gc/2.0) + ((g14r + g23r)*gc)/4.0;
	g14gc = (1-gc)*g14r;
	g22gc = g22r + (g12r*gc)/4.0 + (g24r*gc)/4.0;
	g23gc = (1-gc)*g23r;
	g24gc = g24r*(1 - gc/2.0) + ((g14r + g23r)*gc)/4.0;
	g33gc = g33r + (g13r*gc)/4.0 + (g34r*gc)/4.0;
	g34gc = g34r*(1 - gc/2.0) + ((g14r + g23r)*gc)/4.0;
	g44gc = g44r + (g24r*gc)/4.0 + (g34r*gc)/4.0;

	/* Output */
	*(geninit + 0) = g11gc;
	*(geninit + 1) = g12gc;
	*(geninit + 2) = g13gc;
	*(geninit + 3) = g14gc;
	*(geninit + 4) = g22gc;
	*(geninit + 5) = g23gc;
	*(geninit + 6) = g24gc;
	*(geninit + 7) = g33gc;
	*(geninit + 8) = g34gc;
	*(geninit + 9) = g44gc;
}		/* End of gene conversion routine */

/* Has neutral allele fixed or not?  Measuring freq of B */
double ncheck(double *geninit){
	/* Fed-in genotype frequencies (for ease of programming) */
	double g11s, g12s, g13s, g14s, g22s, g23s, g24s, g33s, g34s, g44s;
	/* Haplotypes ONLY CONTAINING B */
	double x3, x4;
	double Btot = 0;        /* Total frequency of B */

	/* Initial definition of genotypes */
	g11s = *(geninit + 0);
	g12s = *(geninit + 1);
	g13s = *(geninit + 2);
	g14s = *(geninit + 3);
	g22s = *(geninit + 4);
	g23s = *(geninit + 5);
	g24s = *(geninit + 6);
	g33s = *(geninit + 7);
	g34s = *(geninit + 8);
	g44s = *(geninit + 9);

	/* Calculation of haplotypes containing B */

	x3 = g33s + (g13s + g23s + g34s)/2.0;
	x4 = g44s + (g14s + g24s + g34s)/2.0;

	/* Checking */
	Btot = x3 + x4;
	return Btot;

}	/* End of B check routine */

/* Checking if balancing polymorphism lost or not */
double pcheck(double *geninit){
	/* Fed-in genotype frequencies (for ease of programming) */
	double g11s, g12s, g13s, g14s, g22s, g23s, g24s, g33s, g34s, g44s;
	double x2, x4;
	double Atot = 0;        /* Total frequency of A */

	/* Initial definition of genotypes */
	g11s = *(geninit + 0);
	g12s = *(geninit + 1);
	g13s = *(geninit + 2);
	g14s = *(geninit + 3);
	g22s = *(geninit + 4);
	g23s = *(geninit + 5);
	g24s = *(geninit + 6);
	g33s = *(geninit + 7);
	g34s = *(geninit + 8);
	g44s = *(geninit + 9);

	/* Calculation of haplotypes containing B */

	x2 = g22s + (g12s + g23s + g24s)/2.0;
	x4 = g44s + (g14s + g24s + g34s)/2.0;

	/* Checking */
	Atot = x2 + x4;
	return Atot;

}	/* End of A check routine */

/* End of program */

max_stars_repo_name: MattHartfield/BalSelSims

max_issues_repo_name: MattHartfield/BalSelSims

max_forks_repo_name: MattHartfield/BalSelSims


text:
/*
 *  Copyright (C) 2021 FISCO BCOS.
 *  SPDX-License-Identifier: Apache-2.0
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 * @file LogEntry.h
 * @author: yujiechen
 * @date: 2021-03-18
 */
#pragma once
#include "../interfaces/crypto/CryptoSuite.h"
#include "../libutilities/Common.h"
#include "../libutilities/FixedBytes.h"
#include <gsl/span>

namespace bcos
{
namespace protocol
{
class LogEntry
{
public:
    using Ptr = std::shared_ptr<LogEntry>;
    LogEntry() = default;
    LogEntry(bytes const& _address, h256s _topics, bytes _data)
      : m_address(_address), m_topics(std::move(_topics)), m_data(std::move(_data))
    {}

    ~LogEntry() {}

    std::string_view address() const
    {
        return std::string_view((char*)m_address.data(), m_address.size());
    }
    gsl::span<const h256> topics() const { return gsl::span(m_topics.data(), m_topics.size()); }
    bytesConstRef data() const { return ref(m_data); }
    // Define the scale decode method, which cannot be modified at will
    template <class Stream, typename = std::enable_if_t<Stream::is_decoder_stream>>
    friend Stream& operator>>(Stream& _stream, LogEntry& _logEntry)
    {
        return _stream >> _logEntry.m_address >> _logEntry.m_topics >> _logEntry.m_data;
    }

    // Define the scale encode method, which cannot be modified at will
    template <class Stream, typename = std::enable_if_t<Stream::is_encoder_stream>>
    friend Stream& operator<<(Stream& _stream, LogEntry const& _logEntry)
    {
        return _stream << _logEntry.m_address << _logEntry.m_topics << _logEntry.m_data;
    }

private:
    bcos::bytes m_address;
    bcos::h256s m_topics;
    bytes m_data;
};

using LogEntries = std::vector<LogEntry>;
using LogEntriesPtr = std::shared_ptr<std::vector<LogEntry>>;
}  // namespace protocol
}  // namespace bcos