text:
#pragma once

#include "Library.h"

#include <gsl/span>

#include <cstdint>
#include <memory>
#include <optional>
#include <utility>
#include <variant>
#include <vector>

namespace CesiumGeometry {

namespace AvailabilityUtilities {
uint8_t countOnesInByte(uint8_t _byte);
uint32_t countOnesInBuffer(gsl::span<const std::byte> buffer);
} // namespace AvailabilityUtilities

struct CESIUMGEOMETRY_API ConstantAvailability {
  bool constant;
};

struct CESIUMGEOMETRY_API SubtreeBufferView {
  uint32_t byteOffset;
  uint32_t byteLength;
  uint8_t buffer;
};

typedef std::variant<ConstantAvailability, SubtreeBufferView> AvailabilityView;

struct CESIUMGEOMETRY_API AvailabilitySubtree {
  AvailabilityView tileAvailability;
  AvailabilityView contentAvailability;
  AvailabilityView subtreeAvailability;
  std::vector<std::vector<std::byte>> buffers;
};

/**
 * @brief Availability nodes wrap subtree objects and link them together to
 * form a downwardly traversable availability tree.
 */
struct CESIUMGEOMETRY_API AvailabilityNode {
  /**
   * @brief The subtree data for this node.
   *
   * If a node exists but its subtree does not exist, it indicates that the
   * subtree is known to be available and is actively in the process of loading.
   */
  std::optional<AvailabilitySubtree> subtree;

  /**
   * @brief The child nodes for this subtree node.
   */
  std::vector<std::unique_ptr<AvailabilityNode>> childNodes;

  /**
   * @brief Creates an empty instance;
   */
  AvailabilityNode() noexcept;

  /**
   * @brief Sets the loaded subtree for this availability node.
   *
   * @param subtree_ The loaded subtree to set for this node.
   * @param maxChildrenSubtrees The maximum number of children this subtree
   * could possible have if all of them happen to be available.
   */
  void setLoadedSubtree(
      AvailabilitySubtree&& subtree_,
      uint32_t maxChildrenSubtrees) noexcept;
};

struct CESIUMGEOMETRY_API AvailabilityTree {
  std::unique_ptr<AvailabilityNode> pRoot;
};

class CESIUMGEOMETRY_API AvailabilityAccessor {
public:
  AvailabilityAccessor(
      const AvailabilityView& view,
      const AvailabilitySubtree& subtree) noexcept;

  bool isBufferView() const noexcept {
    return pBufferView != nullptr && bufferAccessor;
  }

  bool isConstant() const noexcept { return pConstant != nullptr; }

  /**
   * @brief Unsafe if isConstant is false.
   */
  bool getConstant() const { return pConstant->constant; }

  /**
   * @brief Unsafe is isBufferView is false.
   */
  const gsl::span<const std::byte>& getBufferAccessor() const {
    return *bufferAccessor;
  }

  /**
   * @brief Unsafe if isBufferView is false.
   */
  const std::byte& operator[](size_t i) const {
    return bufferAccessor.value()[i];
  }

  /**
   * @brief Unsafe if isBufferView is false;
   */
  size_t size() const { return pBufferView->byteLength; }

private:
  const SubtreeBufferView* pBufferView;
  const ConstantAvailability* pConstant;
  std::optional<gsl::span<const std::byte>> bufferAccessor;
};
} // namespace CesiumGeometry

max_stars_repo_name: yieryi/cesium-native

max_issues_repo_name: yieryi/cesium-native

max_forks_repo_name: yieryi/cesium-native


text:
/*
** LISA single subject statistical inference using 1st level GLM (general linear modeling).
** The prewhitening approach is adapted from Worsley et al (2002), Neuroimage 15(1).
**
** G.Lohmann, MPI-KYB, 2018
*/
#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>

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

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

extern gsl_matrix *PseudoInv(gsl_matrix *A,gsl_matrix *B);
extern void printmat(gsl_matrix *R,char *str);
extern void printvec(gsl_vector *x,char *str);

extern void prewhite(gsl_matrix_float** pinvX,gsl_matrix_float** invM,gsl_matrix_float* X,int);

extern void dfs(VFloat** Dfs,gsl_matrix_float* pinvX,gsl_matrix_float* X,gsl_matrix_float* con,int);

extern void whitecov(VImage rho_vol,gsl_matrix_float* Y,gsl_matrix_float* invM,
		     gsl_matrix_float* pinvX,gsl_matrix_float* X,int,int);

extern void whitecov2(VImage effect_image, VImage rho_vol,
		      gsl_matrix_float* Y, gsl_matrix_float* X, gsl_matrix_float* con,
		      VFloat* Dfs, int, int slice);

void CopySlice(gsl_matrix *Data,gsl_matrix_float *Y,VImage map,int slice)
{
  size_t nvox;
  size_t i=0,j=0;
  double u=0;
  int ncols = VPixel(map,0,3,2,VShort);

  gsl_matrix_float_set_zero(Y);

  for (nvox=0; nvox < Data->size1; nvox++) {
    int b = VPixel(map,0,0,nvox,VShort);
    if (b != slice) continue;
    int r = VPixel(map,0,1,nvox,VShort);
    int c = VPixel(map,0,2,nvox,VShort);
    i = c + ncols*r;
    if (i >= Y->size2) continue;

    for (j=0; j<Data->size2; j++) {
      u = gsl_matrix_get(Data,nvox,j);
      gsl_matrix_float_set(Y,j,i,(float)u);
    }
  }
}

gsl_matrix_float *GslFloatCpy(gsl_matrix *X)
{
  int i,j;
  double u;
  gsl_matrix_float *Z = gsl_matrix_float_calloc(X->size2,X->size1);
  for (i=0; i<X->size1; i++) {
    for (j=0; j<X->size2; j++) {
      u = gsl_matrix_get(X,i,j);
      gsl_matrix_float_set(Z,j,i,(float)u);
    }
  }
  return Z;
}

/*
** general linear model (GLM) using whitening
*/
void VWhiteGLM(gsl_matrix *Data,VImage map,gsl_matrix *X0,gsl_vector *con,int numlags,VImage zmap)
{
  int i,slice;
  int numcon=1;
  int nslices = VPixel(map,0,3,0,VShort);
  int nrows = VPixel(map,0,3,1,VShort);
  int ncols = VPixel(map,0,3,2,VShort);
  int npix = (nrows*ncols);
  int ntimesteps = Data->size2;

  gsl_set_error_handler_off();
  VFillImage(zmap,VAllBands,0);

  /* contrast vector */
  gsl_matrix_float *contrast = gsl_matrix_float_alloc(1,con->size);
  for (i=0; i<con->size; i++) gsl_matrix_float_set(contrast,0,i,(float)con->data[i]);

  /* copy design matrix */
  gsl_matrix_float *X = GslFloatCpy(X0);

  /* initialize rho_vol */
  VImage rho_vol = VCreateImage(npix,nslices,numlags,VFloatRepn);
  VFillImage(rho_vol,VAllBands,0);

  /* prewhitening */
  gsl_matrix_float *pinvX=NULL,*invM=NULL;
  prewhite(&pinvX,&invM,X,numlags);

  /* degrees of freedom */
  VFloat *dfs_gsl = (VFloat *)VMalloc(sizeof(VFloat)* (numcon+1));
  dfs(&dfs_gsl,pinvX,X,contrast,numlags);

  /* first pass */
  gsl_matrix_float *Y = gsl_matrix_float_calloc(ntimesteps,npix);
  for (slice=0; slice<nslices; slice++) {
    CopySlice(Data,Y,map,(int)slice);
    whitecov(rho_vol,Y,invM,pinvX,X,numlags,slice);
  }

  /* second pass */
  VImage dest = VCreateImage(nslices,nrows,ncols,VFloatRepn);
  VFillImage(dest,VAllBands,0);
  for (slice=0; slice<nslices; slice++) {
    CopySlice(Data,Y,map,(int)slice);
    whitecov2(dest,rho_vol,Y,X,contrast,dfs_gsl,numlags,slice);
  }

  /* output */
  zmap = VCopyImagePixels(dest,zmap,VAllBands);

  /* cleanup */
  gsl_matrix_float_free(contrast);
  gsl_matrix_float_free(X);
  gsl_matrix_float_free(Y);
  gsl_matrix_float_free(pinvX);
  gsl_matrix_float_free(invM);
  VFree(dfs_gsl);
  VDestroyImage(rho_vol);
  VDestroyImage(dest);
}

ext: c

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


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

// #include <stdlib.h>
// #include <time.h>

double boxcar_overlap(double rA, double radiiA, double rB, double radiiB)
{

    double F;
    F = pow(rB - rA, 2) / pow(radiiA + radiiB, 2);

    return F;

}

size_t gen_pts_rsa_1d(double *x,
    size_t npoints, double radius, int step_limit,
    unsigned long randSeed)
{

    // Setup GSL random number generator
    const gsl_rng_type * T;
    gsl_rng * r;
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    // Set the seed
    // srand ( time(NULL) );
    // unsigned long randSeed = rand();
    gsl_rng_set(r, randSeed);

    // Set the initial position
    double xn = gsl_rng_uniform (r) * (1 - 2 * radius) + radius;
    x[0] = xn;

    size_t valid_pts;
    double F;
    int k, flag, step;

    step = 0;
    valid_pts = 1;

    while ((valid_pts < npoints) & (step < step_limit))
    {

        xn = gsl_rng_uniform (r) * (1 - 2 * radius) + radius;

        flag = 1;
        for (k = 0; k < valid_pts; k++)
        {

            F = boxcar_overlap(xn, radius, x[k], radius);
            if (F < 1.0)
            {

                flag = 0;
                break;

            }
        }
        if (flag == 1)
        {

           x[valid_pts] = xn;
           valid_pts += 1;

        }

        step += 1;

    }

    gsl_rng_free (r);

    return valid_pts;

}

unsigned int metro_md_1d(double *x,
    double radius, size_t npoints, unsigned int step_limit,
    unsigned long randSeed)
{

    /* Setup GSL random number generator */
    const gsl_rng_type * T;
    gsl_rng * r;
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    /* Set the seed */
    // srand ( time(NULL) );
    // unsigned long randSeed = rand();
    gsl_rng_set(r, randSeed);

    double diameter = 2 * radius;

    double dx, xn, F;
    unsigned int step, i, k, flag, success_steps;

    step = 0;
    success_steps = 0;
    while (step < step_limit)
    {

        i = step % npoints;

        /* Generate new position */
        while (1)
        {

            dx = diameter * (gsl_rng_uniform (r) - 0.5);
            xn = x[i] + dx;

            if (((xn > radius) & (xn < 1 - radius)))
            {
               break;
            }

        }

        /* Determine if new position overlaps with other positions */
        flag = 1;
        for (k = 0; k < npoints; k++)
        {

            F = boxcar_overlap(xn, radius, x[k], radius);
            if ((F < 1.0) & (i != k))
            {
                flag = 0;
                break;
            }

        }

        if (flag == 1)
        {

            x[i] = xn;

            success_steps = success_steps + 1;

        }

        step = step + 1;

    }

    gsl_rng_free (r);

    return success_steps;

}

unsigned int metro_pd_1d(double *x,
    double *radius, size_t npoints, int step_limit,
    unsigned long randSeed)
{

    /* Setup GSL random number generator */
    const gsl_rng_type * T;
    gsl_rng * r;
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    /* Set the seed */
    // srand ( time(NULL) );
    // unsigned long randSeed = rand();
    gsl_rng_set(r, randSeed);

    double dx, xn, diameter, F;
    unsigned int step, i, k, flag, success_steps;

    step = 0;
    success_steps = 0;
    while (step < step_limit)
    {

        i = step % npoints;

        /* Generate new position */
        while (1)
        {

            diameter = 2 * radius[i];

            dx = diameter * (gsl_rng_uniform (r) - 0.5);
            xn = x[i] + dx;

            if (((xn > radius[i]) & (xn < 1 - radius[i])) )
            {
               break;
            }

        }

        /* Determine if new position overlaps with other positions */
        flag = 1;
        for (k = 0; k < npoints; k++)
        {

            F = boxcar_overlap(xn, radius[i], x[k], radius[k]);
            if ((F < 1.0) & (i != k))
            {
                flag = 0;
                break;
            }

        }

        if (flag == 1)
        {

            x[i] = xn;

            success_steps += 1;

        }

        step = step + 1;

    }

    gsl_rng_free (r);

    return success_steps;

}


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

double f(double x, void *p) {
    (void) (p); /* avoid unused parameter warning */
    return pow(0.5, pow(x, 2) + 2 * x);
}

int main(void) {
    int i, j, n = 10000;
    gsl_function F;
    F.function = f;
    F.params = 0;
    gsl_cheb_series *cs = gsl_cheb_alloc(40);
    FILE *file;
    double *x = malloc(2 * n * sizeof(double));
    double *w = malloc(2 * n * sizeof(double));
    double *y = malloc(2 * n * sizeof(double));
    for (i = -n, j = 0; i < n; i++, j++) {
        double xf = i / (double) n;
        x[j] = xf;
        w[j] = 1;
        y[j] = GSL_FN_EVAL (&F, xf);
    }
    double c0, c1, cov00, cov01, cov11, chisq;
    gsl_fit_wlinear(x, 1, w, 1, y, 1, n,
                    &c0, &c1, &cov00, &cov01, &cov11,
                    &chisq);
    gsl_cheb_init(cs, &F, -1.0, 1.0);
    file = fopen("dalsze.dat", "w");
    for (i = -n; i < n; i++) {
        double xf = i / (double) n;
        double yf, yf_err;
        gsl_fit_linear_est(xf,
                           c0, c1,
                           cov00, cov01, cov11,
                           &yf, &yf_err);
        double r10 = gsl_cheb_eval_n(cs, 10, xf);
        double r40 = gsl_cheb_eval(cs, xf);
        fprintf(file, "%g %g %g %g %g\n",
                xf, GSL_FN_EVAL (&F, xf), yf, r10, r40);
    }
    free(x);
    free(y);
    free(w);
    gsl_cheb_free(cs);
    fclose(file);
    return EXIT_SUCCESS;
}

max_stars_repo_name: mistyfiky/agh-mownit

max_issues_repo_name: mistyfiky/agh-mownit

max_forks_repo_name: mistyfiky/agh-mownit


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>

#include "allvars.h"
#include "proto.h"

#ifdef SUBFIND
#include "subfind.h"

/*! Structure for communication during the density computation. Holds data that is sent to other processors.
 */
static struct nearestdata_in
{
  MyDouble Pos[3];
  MyIDType ID;
  MyFloat Hsml;
  MyFloat Density;
  MyFloat Dist[2];
  int Count;
  long long Index[2];
  int NodeList[NODELISTLENGTH];
}
 *NearestDataIn, *NearestDataGet;

static struct nearestdata_out
{
  MyFloat Dist[2];
  long long Index[2];
  int Count;
}
 *NearestDataResult, *NearestDataOut;

void subfind_find_nearesttwo(void)
{
  int i, j, k, l, ndone, ndone_flag, dummy;
  int ngrp, sendTask, recvTask, place, nexport, nimport, maxexport;

  if(ThisTask == 0)
    printf("Start finding nearest two (%d particles on task=%d)\n", NumPartGroup, ThisTask);

  /* allocate buffers to arrange communication */

  Ngblist = (int *) mymalloc(NumPartGroup * sizeof(int));
  Dist2list = (double *) mymalloc(NumPartGroup * sizeof(double));

  All.BunchSize =
    (int) ((All.BufferSize * 1024 * 1024) / (sizeof(struct data_index) + sizeof(struct data_nodelist) +
					     sizeof(struct nearestdata_in) + sizeof(struct nearestdata_out) +
					     sizemax(sizeof(struct nearestdata_in),
						     sizeof(struct nearestdata_out))));
  DataIndexTable = (struct data_index *) mymalloc(All.BunchSize * sizeof(struct data_index));
  DataNodeList = (struct data_nodelist *) mymalloc(All.BunchSize * sizeof(struct data_nodelist));
  maxexport = All.BunchSize - 2 * NTopleaves / NODELISTLENGTH;

  for(i = 0; i < NumPartGroup; i++)
    NgbLoc[i].count = 0;

  /* we will repeat the whole thing for those particles where we didn't find enough neighbours */

  i = 0;			/* begin with this index */

  do
    {
      for(j = 0; j < NTask; j++)
	{
	  Send_count[j] = 0;
	  Exportflag[j] = -1;
	}

      /* do local particles and prepare export list */

      for(nexport = 0; i < NumPartGroup; i++)
	{
	  if(subfind_nearesttwo_evaluate(i, 0, &nexport, Send_count) < 0)
	    break;
	}

      qsort(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);

      MPI_Allgather(Send_count, NTask, MPI_INT, Sendcount_matrix, NTask, MPI_INT, MPI_COMM_WORLD);

      for(j = 0, nimport = 0, Recv_offset[0] = 0, Send_offset[0] = 0; j < NTask; j++)
	{
	  Recv_count[j] = Sendcount_matrix[j * NTask + ThisTask];
	  nimport += Recv_count[j];

	  if(j > 0)
	    {
	      Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
	      Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
	    }
	}

      NearestDataGet = (struct nearestdata_in *) mymalloc(nimport * sizeof(struct nearestdata_in));
      NearestDataIn = (struct nearestdata_in *) mymalloc(nexport * sizeof(struct nearestdata_in));

      /* prepare particle data for export */
      for(j = 0; j < nexport; j++)
	{
	  place = DataIndexTable[j].Index;

	  NearestDataIn[j].Pos[0] = P[place].Pos[0];
	  NearestDataIn[j].Pos[1] = P[place].Pos[1];
	  NearestDataIn[j].Pos[2] = P[place].Pos[2];
	  NearestDataIn[j].Hsml = P[place].DM_Hsml;
	  NearestDataIn[j].ID = P[place].ID;
	  NearestDataIn[j].Density = P[place].u.DM_Density;
	  NearestDataIn[j].Count = NgbLoc[place].count;
	  for(k = 0; k < NgbLoc[place].count; k++)
	    {
	      NearestDataIn[j].Dist[k] = R2Loc[place].dist[k];
	      NearestDataIn[j].Index[k] = NgbLoc[place].index[k];
	    }

	  memcpy(NearestDataIn[j].NodeList,
		 DataNodeList[DataIndexTable[j].IndexGet].NodeList, NODELISTLENGTH * sizeof(int));
	}

      /* exchange particle data */
      for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	{
	  sendTask = ThisTask;
	  recvTask = ThisTask ^ ngrp;

	  if(recvTask < NTask)
	    {
	      if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		{
		  /* get the particles */
		  MPI_Sendrecv(&NearestDataIn[Send_offset[recvTask]],
			       Send_count[recvTask] * sizeof(struct nearestdata_in), MPI_BYTE,
			       recvTask, TAG_DENS_A,
			       &NearestDataGet[Recv_offset[recvTask]],
			       Recv_count[recvTask] * sizeof(struct nearestdata_in), MPI_BYTE,
			       recvTask, TAG_DENS_A, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		}
	    }
	}

      myfree(NearestDataIn);
      NearestDataResult = (struct nearestdata_out *) mymalloc(nimport * sizeof(struct nearestdata_out));
      NearestDataOut = (struct nearestdata_out *) mymalloc(nexport * sizeof(struct nearestdata_out));

      /* now do the particles that were sent to us */
      for(j = 0; j < nimport; j++)
	subfind_nearesttwo_evaluate(j, 1, &dummy, &dummy);

      if(i >= NumPartGroup)
	ndone_flag = 1;
      else
	ndone_flag = 0;

      MPI_Allreduce(&ndone_flag, &ndone, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

      /* get the result */
      for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	{
	  sendTask = ThisTask;
	  recvTask = ThisTask ^ ngrp;
	  if(recvTask < NTask)
	    {
	      if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		{
		  /* send the results */
		  MPI_Sendrecv(&NearestDataResult[Recv_offset[recvTask]],
			       Recv_count[recvTask] * sizeof(struct nearestdata_out),
			       MPI_BYTE, recvTask, TAG_DENS_B,
			       &NearestDataOut[Send_offset[recvTask]],
			       Send_count[recvTask] * sizeof(struct nearestdata_out),
			       MPI_BYTE, recvTask, TAG_DENS_B, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		}
	    }
	}

      /* add the result to the local particles */
      for(j = 0; j < nexport; j++)
	{
	  place = DataIndexTable[j].Index;

	  for(k = 0; k < NearestDataOut[j].Count; k++)
	    {
	      if(NgbLoc[place].count >= 1)
		if(NgbLoc[place].index[0] == NearestDataOut[j].Index[k])
		  continue;

	      if(NgbLoc[place].count == 2)
		if(NgbLoc[place].index[1] == NearestDataOut[j].Index[k])
		  continue;

	      if(NgbLoc[place].count < 2)
		{
		  l = NgbLoc[place].count;
		  NgbLoc[place].count++;
		}
	      else
		{
		  if(R2Loc[place].dist[0] > R2Loc[place].dist[1])
		    l = 0;
		  else
		    l = 1;

		  if(NearestDataOut[j].Dist[k] >= R2Loc[place].dist[l])
		    continue;
		}

	      R2Loc[place].dist[l] = NearestDataOut[j].Dist[k];
	      NgbLoc[place].index[l] = NearestDataOut[j].Index[k];

	      if(NgbLoc[place].count == 2)
		if(NgbLoc[place].index[0] == NgbLoc[place].index[1])
		  {
		    /*
		       printf("taaa=%d i=%d  task_0=%d index_0=%d  task_1=%d index_1=%d\n",
		       ThisTask, place,
		       NgbLoc[place].task[0], NgbLoc[place].index[0],
		       NgbLoc[place].task[1], NgbLoc[place].index[1]);

		       printf
		       ("NearestDataOut[j].Count=%d  l=%d k=%d task_0=%d index_0=%d  task_1=%d index_1=%d\n",
		       NearestDataOut[j].Count, l, k, NearestDataOut[j].Task[0], NearestDataOut[j].Index[0],
		       NearestDataOut[j].Task[1], NearestDataOut[j].Index[1]);
		     */

		    endrun(112);
		  }
	    }
	}

      myfree(NearestDataOut);
      myfree(NearestDataResult);
      myfree(NearestDataGet);
    }
  while(ndone < NTask);

  myfree(DataNodeList);
  myfree(DataIndexTable);

  myfree(Dist2list);
  myfree(Ngblist);
}

/*! This function represents the core of the SPH density computation. The
 *  target particle may either be local, or reside in the communication
 *  buffer.
 */
int subfind_nearesttwo_evaluate(int target, int mode, int *nexport, int *nsend_local)
{
  int j, k, n, count;
  MyIDType ID;
  long long index[2];
  double dist[2];
  int startnode, numngb, numngb_inbox, listindex = 0;
  double hmax;
  double h, h2, hinv, hinv3;
  double r2, density;
  MyDouble *pos;

  numngb = 0;

  if(mode == 0)
    {
      ID = P[target].ID;
      density = P[target].u.DM_Density;
      pos = P[target].Pos;
      h = P[target].DM_Hsml;
      count = NgbLoc[target].count;
      for(k = 0; k < count; k++)
	{
	  dist[k] = R2Loc[target].dist[k];
	  index[k] = NgbLoc[target].index[k];
	}
    }
  else
    {
      ID = NearestDataGet[target].ID;
      density = NearestDataGet[target].Density;
      pos = NearestDataGet[target].Pos;
      h = NearestDataGet[target].Hsml;
      count = NearestDataGet[target].Count;
      for(k = 0; k < count; k++)
	{
	  dist[k] = NearestDataGet[target].Dist[k];
	  index[k] = NearestDataGet[target].Index[k];
	}
    }

  h2 = h * h;
  hinv = 1.0 / h;
  hinv3 = hinv * hinv * hinv;

  if(mode == 0)
    {
      startnode = All.MaxPart;	/* root node */
    }
  else
    {
      startnode = NearestDataGet[target].NodeList[0];
      startnode = Nodes[startnode].u.d.nextnode;	/* open it */
    }

  numngb = 0;

  if(count == 2)
    if(index[0] == index[1])
      {
	/*      printf("T-A-S-K=%d target=%d mode=%d  task_0=%d index_0=%d  task_1=%d index_1=%d\n",
	   ThisTask, target, mode, task[0], index[0], task[1], index[1]);
	 */
	endrun(1232);
      }

  while(startnode >= 0)
    {
      while(startnode >= 0)
	{
	  numngb_inbox =
	    subfind_ngb_treefind_nearesttwo(pos, h, target, &startnode, mode, &hmax, nexport, nsend_local);

	  if(numngb_inbox < 0)
	    return -1;

	  for(n = 0; n < numngb_inbox; n++)
	    {
	      j = Ngblist[n];
	      r2 = Dist2list[n];

	      if(P[j].ID != ID)	/* exclude the self-particle */
		{
		  if(P[j].u.DM_Density > density)	/* we only look at neighbours that are denser */
		    {
		      if(count < 2)
			{
			  dist[count] = r2;
			  index[count] = (((long long) ThisTask) << 32) + j;
			  count++;
			}
		      else
			{
			  if(dist[0] > dist[1])
			    k = 0;
			  else
			    k = 1;

			  if(r2 < dist[k])
			    {
			      dist[k] = r2;
			      index[k] = (((long long) ThisTask) << 32) + j;
			    }
			}
		    }

		  numngb++;
		}
	    }
	}

      if(mode == 1)
	{
	  listindex++;
	  if(listindex < NODELISTLENGTH)
	    {
	      startnode = NearestDataGet[target].NodeList[listindex];
	      if(startnode >= 0)
		startnode = Nodes[startnode].u.d.nextnode;	/* open it */
	    }
	}
    }

  if(mode == 0)
    {
      NgbLoc[target].count = count;
      for(k = 0; k < count; k++)
	{
	  R2Loc[target].dist[k] = dist[k];
	  NgbLoc[target].index[k] = index[k];
	}

      if(count == 2)
	if(NgbLoc[target].index[0] == NgbLoc[target].index[1])
	  {
	    /*
	       printf("TASK=%d target=%d  task_0=%d index_0=%d  task_1=%d index_1=%d\n",
	       ThisTask, target,
	       NgbLoc[target].task[0], NgbLoc[target].index[0],
	       NgbLoc[target].task[1], NgbLoc[target].index[1]);
	     */
	    endrun(2);
	  }
    }
  else
    {
      NearestDataResult[target].Count = count;
      for(k = 0; k < count; k++)
	{
	  NearestDataResult[target].Dist[k] = dist[k];
	  NearestDataResult[target].Index[k] = index[k];
	}

      if(count == 2)
	if(NearestDataResult[target].Index[0] == NearestDataResult[target].Index[1])
	  {
	    /*
	       printf("TASK!!=%d target=%d  task_0=%d index_0=%d  task_1=%d index_1=%d\n",
	       ThisTask, target,
	       NearestDataResult[target].Task[0], NearestDataResult[target].Index[0],
	       NearestDataResult[target].Task[1], NearestDataResult[target].Index[1]);
	     */
	    endrun(22);
	  }
    }

  return 0;
}

int subfind_ngb_treefind_nearesttwo(MyDouble searchcenter[3], double hsml, int target, int *startnode,
				    int mode, double *hmax, int *nexport, int *nsend_local)
{
  int numngb, no, p, task, nexport_save, exported = 0;
  struct NODE *current;
  double dx, dy, dz, dist, r2;

#ifdef PERIODIC
  double xtmp;
#endif

  nexport_save = *nexport;

  *hmax = 0;
  numngb = 0;
  no = *startnode;

  while(no >= 0)
    {
      if(no < All.MaxPart)	/* single particle */
	{
	  p = no;
	  no = Nextnode[no];

#ifdef DENSITY_SPLIT_BY_TYPE
	  if(!((1 << P[p].Type) & (DENSITY_SPLIT_BY_TYPE)))
	    continue;
#else
	  if(!((1 << P[p].Type) & (FOF_PRIMARY_LINK_TYPES)))
	    continue;
#endif

	  dist = hsml;
	  dx = NGB_PERIODIC_LONG_X(P[p].Pos[0] - searchcenter[0]);
	  if(dx > dist)
	    continue;
	  dy = NGB_PERIODIC_LONG_Y(P[p].Pos[1] - searchcenter[1]);
	  if(dy > dist)
	    continue;
	  dz = NGB_PERIODIC_LONG_Z(P[p].Pos[2] - searchcenter[2]);
	  if(dz > dist)
	    continue;
	  if((r2 = (dx * dx + dy * dy + dz * dz)) > dist * dist)
	    continue;

	  Dist2list[numngb] = r2;
	  Ngblist[numngb++] = p;
	}
      else
	{
	  if(no >= All.MaxPart + MaxNodes)	/* pseudo particle */
	    {
	      if(mode == 1)
		endrun(12312);

	      if(target >= 0)	/* if no target is given, export will not occur */
		{
		  exported = 1;
		  if(Exportflag[task = DomainTask[no - (All.MaxPart + MaxNodes)]] != target)
		    {
		      Exportflag[task] = target;
		      Exportnodecount[task] = NODELISTLENGTH;
		    }

		  if(Exportnodecount[task] == NODELISTLENGTH)
		    {
		      if(*nexport >= All.BunchSize)
			{
			  *nexport = nexport_save;
			  if(nexport_save == 0)
			    endrun(13004);	/* in this case, the buffer is too small to process even a single particle */
			  for(task = 0; task < NTask; task++)
			    nsend_local[task] = 0;
			  for(no = 0; no < nexport_save; no++)
			    nsend_local[DataIndexTable[no].Task]++;
			  return -1;
			}
		      Exportnodecount[task] = 0;
		      Exportindex[task] = *nexport;
		      DataIndexTable[*nexport].Task = task;
		      DataIndexTable[*nexport].Index = target;
		      DataIndexTable[*nexport].IndexGet = *nexport;
		      *nexport = *nexport + 1;
		      nsend_local[task]++;
		    }

		  DataNodeList[Exportindex[task]].NodeList[Exportnodecount[task]++] =
		    DomainNodeIndex[no - (All.MaxPart + MaxNodes)];

		  if(Exportnodecount[task] < NODELISTLENGTH)
		    DataNodeList[Exportindex[task]].NodeList[Exportnodecount[task]] = -1;
		}

	      no = Nextnode[no - MaxNodes];
	      continue;
	    }

	  current = &Nodes[no];

	  if(mode == 1)
	    {
	      if(current->u.d.bitflags & (1 << BITFLAG_TOPLEVEL))	/* we reached a top-level node again, which means that we are done with the branch */
		{
		  *startnode = -1;
		  return numngb;
		}
	    }

	  no = current->u.d.sibling;	/* in case the node can be discarded */

	  dist = hsml + 0.5 * current->len;;
	  dx = NGB_PERIODIC_LONG_X(current->center[0] - searchcenter[0]);
	  if(dx > dist)
	    continue;
	  dy = NGB_PERIODIC_LONG_Y(current->center[1] - searchcenter[1]);
	  if(dy > dist)
	    continue;
	  dz = NGB_PERIODIC_LONG_Z(current->center[2] - searchcenter[2]);
	  if(dz > dist)
	    continue;
	  /* now test against the minimal sphere enclosing everything */
	  dist += FACT1 * current->len;
	  if(dx * dx + dy * dy + dz * dz > dist * dist)
	    continue;

	  no = current->u.d.nextnode;	/* ok, we need to open the node */
	}
    }

  *startnode = -1;
  return numngb;
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
/* specfunc/beta_inc.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 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 <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_gamma.h>

#include "error.h"
#include "check.h"

static
int
beta_cont_frac(
  const double a,
  const double b,
  const double x,
  gsl_sf_result * result
  )
{
  const unsigned int max_iter = 512;        /* control iterations      */
  const double cutoff = 2.0 * GSL_DBL_MIN;  /* control the zero cutoff */
  unsigned int iter_count = 0;
  double cf;

  /* standard initialization for continued fraction */
  double num_term = 1.0;
  double den_term = 1.0 - (a+b)*x/(a+1.0);
  if (fabs(den_term) < cutoff) den_term = cutoff;
  den_term = 1.0/den_term;
  cf = den_term;

  while(iter_count < max_iter) {
    const int k  = iter_count + 1;
    double coeff = k*(b-k)*x/(((a-1.0)+2*k)*(a+2*k));
    double delta_frac;

    /* first step */
    den_term = 1.0 + coeff*den_term;
    num_term = 1.0 + coeff/num_term;
    if(fabs(den_term) < cutoff) den_term = cutoff;
    if(fabs(num_term) < cutoff) num_term = cutoff;
    den_term  = 1.0/den_term;

    delta_frac = den_term * num_term;
    cf *= delta_frac;

    coeff = -(a+k)*(a+b+k)*x/((a+2*k)*(a+2*k+1.0));

    /* second step */
    den_term = 1.0 + coeff*den_term;
    num_term = 1.0 + coeff/num_term;
    if(fabs(den_term) < cutoff) den_term = cutoff;
    if(fabs(num_term) < cutoff) num_term = cutoff;
    den_term = 1.0/den_term;

    delta_frac = den_term*num_term;
    cf *= delta_frac;

    if(fabs(delta_frac-1.0) < 2.0*GSL_DBL_EPSILON) break;

    ++iter_count;
  }

  result->val = cf;
  result->err = iter_count * 4.0 * GSL_DBL_EPSILON * fabs(cf);

  if(iter_count >= max_iter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_beta_inc_e(
  const double a,
  const double b,
  const double x,
  gsl_sf_result * result
  )
{
  if(a <= 0.0 || b <= 0.0 || x < 0.0 || x > 1.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x == 1.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result ln_beta;
    gsl_sf_result ln_x;
    gsl_sf_result ln_1mx;
    gsl_sf_result prefactor;
    const int stat_ln_beta = gsl_sf_lnbeta_e(a, b, &ln_beta);
    const int stat_ln_1mx = gsl_sf_log_1plusx_e(-x, &ln_1mx);
    const int stat_ln_x = gsl_sf_log_e(x, &ln_x);
    const int stat_ln = GSL_ERROR_SELECT_3(stat_ln_beta, stat_ln_1mx, stat_ln_x);

    const double ln_pre_val = -ln_beta.val + a * ln_x.val + b * ln_1mx.val;
    const double ln_pre_err =  ln_beta.err + fabs(a*ln_x.err) + fabs(b*ln_1mx.err);
    const int stat_exp = gsl_sf_exp_err_e(ln_pre_val, ln_pre_err, &prefactor);

    if(stat_ln != GSL_SUCCESS) {
      result->val = 0.0;
      result->err = 0.0;
      GSL_ERROR ("error", GSL_ESANITY);
    }

    if(x < (a + 1.0)/(a+b+2.0)) {
      /* Apply continued fraction directly. */
      gsl_sf_result cf;
      const int stat_cf = beta_cont_frac(a, b, x, &cf);
      int stat;
      result->val = prefactor.val * cf.val / a;
      result->err = (fabs(prefactor.err * cf.val) + fabs(prefactor.val * cf.err))/a;

      stat = GSL_ERROR_SELECT_2(stat_exp, stat_cf);
      if(stat == GSL_SUCCESS) {
        CHECK_UNDERFLOW(result);
      }
      return stat;
    }
    else {
      /* Apply continued fraction after hypergeometric transformation. */
      gsl_sf_result cf;
      const int stat_cf = beta_cont_frac(b, a, 1.0-x, &cf);
      int stat;
      const double term = prefactor.val * cf.val / b;
      result->val  = 1.0 - term;
      result->err  = fabs(prefactor.err * cf.val)/b;
      result->err += fabs(prefactor.val * cf.err)/b;
      result->err += 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(term));
      stat = GSL_ERROR_SELECT_2(stat_exp, stat_cf);
      if(stat == GSL_SUCCESS) {
        CHECK_UNDERFLOW(result);
      }
      return stat;
    }
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_beta_inc(const double a, const double b, const double x)
{
  EVAL_RESULT(gsl_sf_beta_inc_e(a, b, x, &result));
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
-----------natrix mul---
/*
 * standard includes
 */
	#include <stdio.h>
	#include <stdlib.h>
/*
 *  includes for GSL components
 *  	- use double precision
 */
	#include <gsl/gsl_vector_double.h>
	#include <gsl/gsl_matrix_double.h>
	#include <gsl/gsl_rng.h>
	#include <gsl/gsl_randist.h>

/*
 * FUNCTIONS
 */

/*
 * simple Fibonacci sequence generator function, using recursion
 */
size_t fib(size_t k)
{
	if (k==0)
	{
		return 0;
	}
	else if (k==1)
	{
		return 1;
	}
	else /* k >= 2 */
	{
		return fib(k-1) + fib(k-2);
	}
}
gsl_matrix *embt_mm(const gsl_matrix *U,const gsl_matrix *V, size_t N)
{
	gsl_matrix *W=gsl_matrix_calloc(N,N);
	double sum=0;
	for (size_t i=0;i!=N;i++)
	{
		for(size_t j=0;j!=N;j++)
		{
			for (size_t k=0;k!=N;k++)
			{
				sum=sum+gsl_matrix_get(U,i,k)*gsl_matrix_get(V,k,j);
			}
			gsl_matrix_set(W,i,j,(double)sum);
			sum=0;
		}
	}
	return W;
}

gsl_matrix *embt_apb(gsl_vector *u, gsl_vector *v,size_t N)
{
	gsl_vector *w=gsl_vector_calloc(N);
	gsl_vector_add(w,u);
	gsl_vector_add(w,v);
	return w;
}

gsl_matrix *embt_m_print(const gsl_matrix *U, size_t N)
{
	for (size_t i = 0; i != N; ++i)
	{
		for (size_t j=0;j!=N;j++)
		{
			printf("%0.2f\t",gsl_matrix_get(U, i, j));
		}
		printf("\n");
	}
}
gsl_matrix *embt_v_print(const gsl_vector *u, size_t N)
{
	for (int i = 0; i != N; ++i)
	{
		printf("%0.2f\n",gsl_vector_get(u, i));
	}
}
int main()
{
	/*
	 * INITIALIZE PARAMETERS
	 */
		/* vectors parameters */
		size_t		N=3; /* index type, vector sizes */
		gsl_vector *a = gsl_vector_alloc(N); /* allocate vector from heap of size N */
		gsl_vector *b = gsl_vector_alloc(N); /* allocate vector from heap of size N */
		gsl_vector *c = gsl_vector_calloc(N); /* allocate vector of size N but initialize entries to zero */

		/* random number generator parameters */
		const gsl_rng_type *T;
		gsl_rng *r; /* handle for our random number generator */

		/* matrix parameters */
		gsl_matrix *A = gsl_matrix_alloc(N,N);
		gsl_matrix *B = gsl_matrix_alloc(N,N);
		gsl_matrix *C = gsl_matrix_calloc(N,N);

	/*
	 * SET UP RANDOM NUMBER GENERATION
	 */
		gsl_rng_env_setup();
		T = gsl_rng_default;
		r = gsl_rng_alloc(T);

	/*
	 *  VECTOR OPERATIONS
	 */
		/* set the vector elements */
		for (size_t i = 0; i != N; ++i)
		{
			gsl_vector_set(a, i, fib(i)); /* set element i of vector a to Fibonacci number i */
			gsl_vector_set(b, i, gsl_ran_flat(r,-1.0,+1.0)); /* set element of vector b to random number */
		}

		/* c = a + b */
		c=embt_apb(a,b,N); /* c += a */
		//gsl_vector_add(c, b); /* c += b */

		/* print results */

		int sum=0;
	/*	MATRIX OPERATIONS */
		for (size_t i=0;i!=N;i++)
			for (size_t j=0;j!=N;j++)
			{
				//size_t x=j+(i*N);
				gsl_matrix_set(A,i,j,(double)fib(j+i*N));
				gsl_matrix_set(B,i,j,gsl_ran_flat(r,-100,+100));
			}
		//gsl_matrix_mul_elements(C, A); /* c += a */
		//gsl_matrix_mul_elements(C, B); /* c += b */
		C=embt_mm(A,B,N);

		printf("a= \n");embt_v_print(a,N);
		printf("b= \n");embt_v_print(b,N);
		printf("c= \n");embt_v_print(c,N);
		printf("A= \n");embt_m_print(A,N);
		printf("B= \n");embt_m_print(B,N);
		printf("C= \n");embt_m_print(C,N);

	/* de-allocate memory */
	gsl_vector_free(a);
	gsl_vector_free(b);
	gsl_vector_free(c);
	gsl_matrix_free(A);
	gsl_matrix_free(B);
	gsl_matrix_free(C);
	return EXIT_SUCCESS;
}

max_stars_repo_path: C program/matric_mul_gsl.c

max_issues_repo_path: C program/matric_mul_gsl.c

max_forks_repo_path: C program/matric_mul_gsl.c


text:
/* Calling DGELS using row-major order */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cblas.h>
#include <lapacke.h>
#include <string.h>
#include "magneticTrap.h"
//#include <sys/stat.h>
#define SIZE 3

//Function to calculate beta. The speed is then v=beta*c
double calculateBeta(double T) {
  double eMass = 9.10938356e-31;
  double c = 299792458.0;
  double factor = pow(T,2.0)/(pow(eMass,2.0)*pow(c,4.0)) + (2.0*T)/(eMass*pow(c,2.0)) + 1;
  double beta = sqrt(1.0 - 1.0/factor);
  return beta;
}

int main (int argc, const char * argv[]) {
    int i, j;
    //FILE *inputFile = fopen("simParms.dat","r");
    //for(int kk=0; kk<10; kk++){
    //}
    //
    double lambda,lambda0,Bx,By,Bz,q,dt,posX,posY,posZ;
    double kinE, mass,gamma,c,Ux,Uy,Uz,magV,magVComputed;
    double dirCosX, dirCosY, dirCosZ, rotationR, solenoidWidth;
    double coilRadius, current, mu, coilPosZMin1, coilPosZMax1, coilPosZMin2, coilPosZMax2;
    double Z_coil1_, Z_coil2_ ;
    int    nTsteps, nCoils;

    //double vecRHS[SIZE]  = {0.0, 0.0, 0.0};
    //double Bfield[SIZE]  = {0.0, 0.0, 0.0 };
    //double Bcgrnd[SIZE]  = {0.0, 0.0, 1.0 };
    double *vecRHS, *Bfield, *background;
    vecRHS = (double*)malloc(SIZE*sizeof(double));
    Bfield = (double*)malloc(SIZE*sizeof(double));
    background = (double*)malloc(SIZE*sizeof(double));
    background[0] = 0.0; background[1] = 0.0; background[2] = 1.0;
    lapack_int info,nRowsA,nColsB,nColsARowsB,ldA,ldB,ipvt[SIZE*SIZE];

    dt= 5.0e-12;
    nTsteps = 2.5e6;
    q = 1.602176634e-19;
    c = 299792458.0;
    mass = 9.10938356e-31;
    kinE = 18575.0*q;
    rotationR = 4.6375e-4;
    solenoidWidth = 0.2;
    coilRadius = 0.05;
    current = 20.0;
    nCoils = 11;
    mu = 1.25663706e-6;
    Z_coil1_ = -0.5;
    Z_coil2_ =  0.5;
    coilPosZMin1 = Z_coil1_ - solenoidWidth/2.0;
    coilPosZMax1 = Z_coil1_ + solenoidWidth/2.0;
    coilPosZMin2 = Z_coil2_ - solenoidWidth/2.0;
    coilPosZMax2 = Z_coil2_ + solenoidWidth/2.0;

    posX=rotationR;posY=0.0;posZ=0.0;
    bathTubFieldSolenoids(coilRadius, current, coilPosZMin1, coilPosZMax1, coilPosZMin2, coilPosZMax2, nCoils,
                           mu, posX, posY, posZ, background, Bfield);
    Bx = Bfield[0]; By = Bfield[1]; Bz = Bfield[2];
    //-----------
    magV = calculateBeta(kinE)*c;
    double pitchAngle = 85.0 * M_PI/180.0;
    dirCosX = 0.0; dirCosZ=cos(pitchAngle); dirCosY=1.0-dirCosZ;
    Ux = magV*dirCosX; Uy=magV*dirCosY; Uz=magV*dirCosZ;
    //double ufuture[SIZE] = {0.0, 0.0, 0.0};
	//double upast[SIZE]   = {Ux, Uy, Uz};
    double *ufuture, *upast;
    ufuture = (double*)malloc(SIZE*sizeof(double));
    upast = (double*)malloc(SIZE*sizeof(double));
    upast[0] = Ux; upast[1] = Uy; upast[2] = Uz;
    printf("# Starting simulation with U(x), U(y), U(z) = %8.5e , %8.5e , %8.5e \n", Ux, Uy, Uz);
    gamma = 1.0/(sqrt(1.0-(pow(magV,2.0))/(pow(c,2.0))));
    lambda0= (0.5*q*dt)/(2.0*mass*gamma); // Careful: this definition is only for the first time step
    lambda = (1.0*q*dt)/(2.0*mass*gamma); // Careful: this definition is only for the first time step

    //writing a and b as vectors, but they are in concept the matrices of the system of linear equations
    //double a[SIZE*SIZE]  = { 1.0, lambda0*Bz,-lambda0*By , -lambda0*Bz,1.0, lambda0*Bx ,  lambda0*By,-lambda0*Bx,1.0 };
    //double b[SIZE*SIZE]  = { 1.0,-lambda0*Bz, lambda0*By ,  lambda0*Bz,1.0,-lambda0*Bx , -lambda0*By, lambda0*Bx,1.0 };
    double *a, *b;
    a = (double*)malloc(SIZE*SIZE*sizeof(double));
    b = (double*)malloc(SIZE*SIZE*sizeof(double));

    a[0] = 1.0;         a[1] = lambda0*Bz;  a[2] = -lambda0*By;
    a[3] = -lambda0*Bz; a[4] = 1.0;         a[5] = lambda0*Bx;
    a[6] = lambda0*By;  a[7] = -lambda0*Bx; a[8] = 1.0;

    b[0] = 1.0;         b[1] = -lambda0*Bz; b[2] = lambda0*By;
    b[3] = lambda0*Bz;  b[4] = 1.0;         b[5] = -lambda0*Bx;
    b[6] = -lambda0*By; b[7] = lambda0*Bx;  b[8] = 1.0;

    nRowsA = SIZE;
    nColsB = 1;
    nColsARowsB = SIZE;
    ldA    = nRowsA;
    ldB    = nColsARowsB;

    // Printing LHS matrix
    /*printf("#LHS Matrix (Traspose -> Column major order -----------\t | \n");
    for(i=0;i<SIZE;i++)
    {
        printf("#| \t");
        for(j=0;j<SIZE;j++)
        {
           printf("%.5e \t",a[i*SIZE + j]);
        }
        printf("\t | \n");
    }
    printf("#------------------------------------------------------\n");
    printf("#RHS Matrix (Traspose -> Column major order -----------\t | \n");
    for(i=0;i<SIZE;i++)
    {
        printf("#| \t");
        for(j=0;j<SIZE;j++)
        {
           printf("%.5e \t",b[i*SIZE + j]);
        }
        printf("\t | \n");
    }
    printf("#----------------------------------------------------------\n");
    printf("# Time       | Ux \t        | Uy \t        | Uz \t        | magV \t        | PosX \t        | PosY \t        | PosZ \t        |\n");
    */
    double timeVal = 0.0;
    for (int tt = 0; tt < nTsteps; tt++) {
        timeVal = (double)tt * dt;
        //using header variables from LAPACK and openblas
        cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,
            nRowsA,nColsB,nColsARowsB,1.0,b,ldA,upast,ldB,0.0,vecRHS,nRowsA);
        for (i=0; i<SIZE; i++){ufuture[i] = vecRHS[i];}; // This is because dgesv rewrites ufuture
        info = LAPACKE_dgesv(LAPACK_COL_MAJOR, SIZE, 1, a, ldA, ipvt, ufuture, ldA);
        for (i=0; i<SIZE; i++){upast[i] = ufuture[i];};
        Ux = ufuture[0];
        Uy = ufuture[1];
        Uz = ufuture[2];
        magVComputed = sqrt(pow(Ux,2.0)+pow(Uy,2.0)+pow(Uz,2.0));
        magV = magVComputed;
        gamma = 1.0/(sqrt(1.0-(pow(magV,2.0))/(pow(c,2.0))));
        lambda = (q*dt)/(2.0*mass*gamma);

        bathTubFieldSolenoids(coilRadius, current, coilPosZMin1, coilPosZMax1,
                coilPosZMin2, coilPosZMax2, nCoils, mu, posX, posY, posZ, background, Bfield);
        Bx = Bfield[0]; By = Bfield[1]; Bz = Bfield[2];

        a[0] = 1.0;        a[1] = lambda*Bz;    a[2] = -lambda*By;
        a[3] = -lambda*Bz; a[4] = 1.0;          a[5] = lambda*Bx;
        a[6] = lambda*By;  a[7] = -lambda*Bx;   a[8] = 1.0;

        b[0] = 1.0;        b[1] = -lambda*Bz;   b[2] = lambda*By;
        b[3] = lambda*Bz;  b[4] = 1.0;          b[5] = -lambda*Bx;
        b[6] = -lambda*By; b[7] = lambda*Bx;    b[8] = 1.0;
        // -----------
        posX = posX + dt*Ux;
        posY = posY + dt*Uy;
        posZ = posZ + dt*Uz;
        if ( tt%1 == 0 ) {    //This alters the number of points written to the text file
            printf(" %7.4e \t %8.5e \t %8.5e \t %8.5e \t %8.5e \t %8.5e \t %8.5e \n ",
                timeVal, posX, posY, posZ, ufuture[0], ufuture[1], ufuture[2]);
        }
    }

    //free(a); free(b);
    free(ufuture); free(upast);
    free(vecRHS); free(Bfield); free(background);
    info = 0;
    return(info);
}


text:
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

#include "allvars.h"
#include "proto.h"

#define WORKSIZE 100000

static double R8;
static double r_tophat;

static double AA, BB, CC;
static double nu;
static double Norm;

static int NPowerTable;

static struct pow_table
{
  double logk, logD;
}
 *PowerTable;

double PowerSpec(double k)
{
  double power = 0;

  switch (WhichSpectrum)
    {
    case 1:
      power = PowerSpec_EH(k);
      break;

    case 2:
      power = PowerSpec_Tabulated(k);
      break;

    default:
      power = PowerSpec_Efstathiou(k);
      break;
    }

  power *= pow(k, PrimordialIndex - 1.0);

  return power;
}

void read_power_table(void)
{
  FILE *fd;
  char buf[500];
  double k, p;

  sprintf(buf, FileWithInputSpectrum);

  if(!(fd = fopen(buf, "r")))
    {
      printf("can't read input spectrum in file '%s' on task %d\n", buf, ThisTask);
      FatalError(17);
    }

  NPowerTable = 0;
  do
    {
      if(fscanf(fd, " %lg %lg ", &k, &p) == 2)

	NPowerTable++;
      else
	break;
    }
  while(1);

  fclose(fd);

  if(ThisTask == 0)
    {
      printf("found %d rows in input spectrum table\n", NPowerTable);
      fflush(stdout);
    }

  PowerTable = malloc(NPowerTable * sizeof(struct pow_table));

  sprintf(buf, FileWithInputSpectrum);

  if(!(fd = fopen(buf, "r")))
    {
      printf("can't read input spectrum in file '%s' on task %d\n", buf, ThisTask);
      FatalError(18);
    }

  NPowerTable = 0;
  do
    {
      double p;

      if(fscanf(fd, " %lg %lg ", &k, &p) == 2)
	{
	  PowerTable[NPowerTable].logk = k;
	  PowerTable[NPowerTable].logD = p;
	  NPowerTable++;
	}
      else
	break;
    }
  while(1);

  fclose(fd);

  qsort(PowerTable, NPowerTable, sizeof(struct pow_table), compare_logk);
}

int compare_logk(const void *a, const void *b)
{
  if(((struct pow_table *) a)->logk < (((struct pow_table *) b)->logk))
    return -1;

  if(((struct pow_table *) a)->logk > (((struct pow_table *) b)->logk))
    return +1;

  return 0;
}

void initialize_powerspectrum(void)
{
  double res;

  InitTime = 1 / (1 + Redshift);

  AA = 6.4 / ShapeGamma * (3.085678e24 / UnitLength_in_cm);
  BB = 3.0 / ShapeGamma * (3.085678e24 / UnitLength_in_cm);
  CC = 1.7 / ShapeGamma * (3.085678e24 / UnitLength_in_cm);
  nu = 1.13;

  R8 = 8 * (3.085678e24 / UnitLength_in_cm);	/* 8 Mpc/h */

  if(WhichSpectrum == 2)
    read_power_table();

  if(ReNormalizeInputSpectrum == 0 && WhichSpectrum == 2)
    {
      Norm = 1.0;
      /* tabulated file is already at the initial redshift */
      Dplus = 1.0;
    }
  else
    {
      Norm = 1.0;
      res = TopHatSigma2(R8);
      if(ThisTask == 0 && WhichSpectrum == 2)
	printf("\nNormalization of spectrum in file:  Sigma8 = %g\n", sqrt(res));

      Norm = Sigma8 * Sigma8 / res;

      if(ThisTask == 0 && WhichSpectrum == 2)
	printf("Normalization adjusted to  Sigma8=%g   (Normfac=%g)\n\n", Sigma8, Norm);

      Dplus = GrowthFactor(InitTime, 1.0);
    }
}

double PowerSpec_Tabulated(double k)
{
  double logk, logD, P, kold, u, dlogk, Delta2;
  int binlow, binhigh, binmid;

  kold = k;

  k *= (InputSpectrum_UnitLength_in_cm / UnitLength_in_cm);	/* convert to h/Mpc */

  logk = log10(k);

  if(logk < PowerTable[0].logk || logk > PowerTable[NPowerTable - 1].logk)
    return 0;

  binlow = 0;
  binhigh = NPowerTable - 1;

  while(binhigh - binlow > 1)
    {
      binmid = (binhigh + binlow) / 2;
      if(logk < PowerTable[binmid].logk)
	binhigh = binmid;
      else
	binlow = binmid;
    }

  dlogk = PowerTable[binhigh].logk - PowerTable[binlow].logk;

  if(dlogk == 0)
    FatalError(777);

  u = (logk - PowerTable[binlow].logk) / dlogk;

  logD = (1 - u) * PowerTable[binlow].logD + u * PowerTable[binhigh].logD;

  Delta2 = pow(10.0, logD);

  P = Norm * Delta2 / (4 * M_PI * kold * kold * kold);

  return P;
}

double PowerSpec_Efstathiou(double k)
{
  return Norm * k / pow(1 + pow(AA * k + pow(BB * k, 1.5) + CC * CC * k * k, nu), 2 / nu);
}

double PowerSpec_EH(double k)	/* Eisenstein & Hu */
{
  return Norm * k * pow(tk_eh(k), 2);
}

double tk_eh(double k)		/* from Martin White */
{
  double q, theta, ommh2, a, s, gamma, L0, C0;
  double tmp;
  double omegam, ombh2, hubble;

  /* other input parameters */
  hubble = HubbleParam;

  omegam = Omega;
  ombh2 = OmegaBaryon * HubbleParam * HubbleParam;

  if(OmegaBaryon == 0)
    ombh2 = 0.04 * HubbleParam * HubbleParam;

  k *= (3.085678e24 / UnitLength_in_cm);	/* convert to h/Mpc */

  theta = 2.728 / 2.7;
  ommh2 = omegam * hubble * hubble;
  s = 44.5 * log(9.83 / ommh2) / sqrt(1. + 10. * exp(0.75 * log(ombh2))) * hubble;
  a = 1. - 0.328 * log(431. * ommh2) * ombh2 / ommh2
    + 0.380 * log(22.3 * ommh2) * (ombh2 / ommh2) * (ombh2 / ommh2);
  gamma = a + (1. - a) / (1. + exp(4 * log(0.43 * k * s)));
  gamma *= omegam * hubble;
  q = k * theta * theta / gamma;
  L0 = log(2. * exp(1.) + 1.8 * q);
  C0 = 14.2 + 731. / (1. + 62.5 * q);
  tmp = L0 / (L0 + C0 * q * q);
  return (tmp);
}

double TopHatSigma2(double R)
{
  double result, abserr;
  gsl_integration_workspace *workspace;
  gsl_function F;

  workspace = gsl_integration_workspace_alloc(WORKSIZE);

  F.function = &sigma2_int;

  r_tophat = R;

  gsl_integration_qag(&F, 0, 500.0 * 1 / R,
		      0, 1.0e-4, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);

  gsl_integration_workspace_free(workspace);

  return result;

  /* note: 500/R is here chosen as (effectively) infinity integration boundary */
}

double sigma2_int(double k, void *param)
{
  double kr, kr3, kr2, w, x;

  kr = r_tophat * k;
  kr2 = kr * kr;
  kr3 = kr2 * kr;

  if(kr < 1e-8)
    return 0;

  w = 3 * (sin(kr) / kr3 - cos(kr) / kr2);
  x = 4 * PI * k * k * w * w * PowerSpec(k);

  return x;
}

double GrowthFactor(double astart, double aend)
{
  return growth(aend) / growth(astart);
}

double growth(double a)
{
  double hubble_a;

  hubble_a = sqrt(Omega / (a * a * a) + (1 - Omega - OmegaLambda) / (a * a) + OmegaLambda);

  double result, abserr;
  gsl_integration_workspace *workspace;
  gsl_function F;

  workspace = gsl_integration_workspace_alloc(WORKSIZE);

  F.function = &growth_int;

  gsl_integration_qag(&F, 0, a, 0, 1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);

  gsl_integration_workspace_free(workspace);

  return hubble_a * result;
}

double growth_int(double a, void *param)
{
  return pow(a / (Omega + (1 - Omega - OmegaLambda) * a + OmegaLambda * a * a * a), 1.5);
}

double F_Omega(double a)
{
  double omega_a;

  omega_a = Omega / (Omega + a * (1 - Omega - OmegaLambda) + a * a * a * OmegaLambda);

  return pow(omega_a, 0.6);
}


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

double f (double x, void * params)
{
  return pow (x, 1.5);
}

int
main (void)
{
  gsl_function F;
  double result, abserr;

  F.function = &f;
  F.params = 0;

  printf ("f(x) = x^(3/2)\n");

  gsl_deriv_central (&F, 2.0, 1e-8, &result, &abserr);
  printf ("x = 2.0\n");
  printf ("f'(x) = %.10f +/- %.10f\n", result, abserr);
  printf ("exact = %.10f\n\n", 1.5 * sqrt(2.0));

  gsl_deriv_forward (&F, 0.0, 1e-8, &result, &abserr);
  printf ("x = 0.0\n");
  printf ("f'(x) = %.10f +/- %.10f\n", result, abserr);
  printf ("exact = %.10f\n", 0.0);

  return 0;
}


text:
/*
 * Test uni_qdblflt() in reprand.cc has a uniform distribution. From:
 * http://rosettacode.org/wiki/Verify_distribution_uniformity/Chi-squared_test
 *
 * Compile with :
 *  gcc -std=c99 -Wall -o reprand_test reprand_test.c -lgsl
 */

#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <gsl/gsl_sf_gamma.h>

#define REPRAND_TEST

#include "reprand.cc"

double chi2UniformDistance( const double const *ds, int dslen)
{
    double expected = 0.0;
    double sum = 0.0;
    int k;

    for (k=0; k<dslen; k++)
      {
        expected += ds[k];
      }
    expected /= k;

    for (k=0; k<dslen; k++)
      {
        double x = ds[k] - expected;
        sum += x*x;
      }
    return sum/expected;
}

double chi2Probability(double dof, double distance)
{
  return gsl_sf_gamma_inc_Q( 0.5*dof, 0.5*distance);
}

int chiIsUniform
( const double * const dset, const int dslen, const double significance)
{
    int dof = dslen -1;
    double dist = chi2UniformDistance( dset, dslen);
    return chi2Probability( dof, dist ) > significance;
}

int main(int argc, char **argv)
{
  /* Example sets */
  const double const dset1[] = { 199809., 200665., 199607., 200270., 199649. };
  const double const dset2[] = { 522573., 244456., 139979.,  71531.,  21461. };

  /* Data set for getRan() */
  const uint32_t N=0x10000;
  double dset_getRan[N];
  /* Initialise with getRan() */
  init_JKISS();
  for(uint32_t k=0;k<N;k++)
    {
      dset_getRan[k]=uni_qdblflt();
    }

  /* Run chi^2 tests. Start with the example sets */
  const double * const dsets[] = { dset1, dset2, dset_getRan };
  const uint32_t const dslens[] = { 5, 5, N };
  for (uint32_t k=0; k<sizeof(dslens)/sizeof(uint32_t); k++)
    {
      double dist = chi2UniformDistance(dsets[k], dslens[k]);
      uint32_t dof = dslens[k]-1;
      double prob = chi2Probability( dof, dist );
      const double significance = 0.001;
      char * const uniform =
        chiIsUniform(dsets[k], dslens[k], significance) ? "Yes":"No";

      printf("dof: %d, dist.: %.4f, prob.: %.6f, uniform? %s, sig. %.3f\n",
             dof, dist, prob, uniform, significance);
    }

    return 0;
}

max_stars_repo_name: robertgj/DesignOfIIRFilters

max_issues_repo_name: robertgj/DesignOfIIRFilters

max_forks_repo_name: robertgj/DesignOfIIRFilters


text:
/*
 * Copyright (c) 2018, Aleksi Kurkela, Aleksas Mazeliauskas, Jean-Francois
 * Paquet, Soeren Schlichting and Derek Teaney
 * All rights reserved.
 *
 * KoMPoST is distributed under MIT license;
 * see the LICENSE file that should be present in the root
 * of the source distribution, or alternately available at:
 * https://github.com/KMPST/KoMPoST/
 */

#include <gsl/gsl_math.h>
#include <iostream>
#ifndef M_HBARC
#define M_HBARC 0.197326979
#endif

namespace GreensFunctions {

class LowKArguments;

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

namespace EnergyPerturbations {

namespace FreeStreaming {
namespace CoordinateSpace {
double Gs(double dX, double dT);
double Gv(double dX, double dT);
double Gd(double dX, double dT);
double Gr(double dX, double dT);
}
} // FreeStreaming

namespace KineticTheory {
namespace CoordinateSpace {
double Gs(double dX, double dT, double ScalingVarOut);
double Gv(double dX, double dT, double ScalingVarOut);
double Gd(double dX, double dT, double ScalingVarOut);
double Gr(double dX, double dT, double ScalingVarOut);
}
} // KineticTheory

namespace LowKLimit {
namespace CoordinateSpace {
double Gs(double dX, double dT, const LowKArguments &in);
double Gv(double dX, double dT, const LowKArguments &in);
double Gd(double dX, double dT, const LowKArguments &in);
double Gr(double dX, double dT, const LowKArguments &in);
}
} // LowKLimit

} // EnergyPerturbations

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

namespace MomentumPerturbations {

namespace FreeStreaming {
namespace CoordinateSpace {
double Hv(double dX, double dT);
double Hd(double dX, double dT);
double Hr(double dX, double dT);
double Htd(double dX, double dT);
double Htm(double dX, double dT);
double Htr(double dX, double dT);
}
} // FreeStreaming

namespace KineticTheory {
namespace CoordinateSpace {
double Hv(double dX, double dT, double ScalingVarOut);
double Hd(double dX, double dT, double ScalingVarOut);
double Hr(double dX, double dT, double ScalingVarOut);
double Htd(double dX, double dT, double ScalingVarOut);
double Htm(double dX, double dT, double ScalingVarOut);
double Htr(double dX, double dT, double ScalingVarOut);
}
} // Kinetic Theory

} // MomentumPerturbations

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

void Setup(double Sigma, int NumberOfPoints, int ENERGY_PERTURBATIONS,
           int MOMENTUM_PERTURBATIONS);

void Output(int ENERGY_PERTURBATIONS, int MOMENTUM_PERTURBATIONS);

class LowKArguments {

  // private:
public:
  // Number of DOF
  const int NuG = 16;

  double tauF;  // Final time in fm
  double tauIn; // Initial time in fm

  double E0;
  double EF;

  double TXX0;
  double TXXF;

  double EtaByS;
  double SigmaBG; // Sigma in fm

public:
  LowKArguments(double tf, double tin, double e0, double ef, double txx0,
                double txxf, double etabys, double sigmabg)
      : tauF(tf), tauIn(tin), E0(e0), EF(ef), TXX0(txx0), TXXF(txxf),
        EtaByS(etabys), SigmaBG(sigmabg) {}

  double S(double dX, double dT) const {
    const double &s = SigmaBG;
    double rs = dX / s;
    return 1. / (M_PI * s * s) * exp(-rs * rs);
  }

  double DS(double dX, double dT) const {
    const double &s = SigmaBG;
    double rs = dX / s;
    return dT / (M_PI * s * s * s) * exp(-rs * rs) * (-2. * rs);
  }

  double DSOverR(double dX, double dT) const {
    const double &s = SigmaBG;
    double rs = dX / s;
    return dT * dT / (M_PI * s * s * s * s) * exp(-rs * rs) * (-2.);
  }

  double DDS(double dX, double dT) const {
    const double &s = SigmaBG;
    double rs = dX / s;
    return dT * dT / (M_PI * s * s * s * s) * exp(-rs * rs) *
           (4 * rs * rs - 2.);
  }

  double Ass() const {
    double dT = tauF - tauIn;
    double s2 = SigmaBG * SigmaBG / (dT * dT);
    return (-0.5 + s2 / 2.) * E0 / EF * (EF + TXXF) / (E0 + TXX0);
  }

  double Asv() const { return 0.5 * E0 / EF * (EF + TXXF) / (E0 + TXX0); }

  double ScalingVariableInverse() const {
    return 3. / 4. * EtaByS * Entropy(EF) * M_HBARC / (EF * tauF);
  }

  double Entropy(double e) const {
    double T = pow(30. / (M_PI * M_PI * NuG) * e, 0.25);
    return 4. / 3. * e / T;
  }
};

} // GreensFunctions

ext: h

max_stars_repo_name: j-f-paquet/kompost

max_issues_repo_name: j-f-paquet/kompost

max_forks_repo_name: j-f-paquet/kompost


text:
#pragma once

#include "concepts/MemoryIntrospection.h"
#include "math/Vector3.h"
#include "pointcloud/PointAttributes.h"

#include <gsl/gsl>
#include <optional>
#include <vector>

// TECH_DEBT Make attributes more dynamic (map<AttributeType, GenericAttribute*>
// or something like that)

/// <summary>
/// Buffer structure that stores point attributes (position, color etc.) for
/// multiple points at once in a structure-of-array fashion. Compared to storing
/// all points as Point structures, this has better performance for data access
/// </summary>
struct PointBuffer
{
  struct PointConstReference;

  /// <summary>
  /// Mutable indirect reference to a single point inside the PointBuffer
  /// </summary>
  struct PointReference
  {
    friend struct PointBuffer;
    friend struct PointConstReference;

    PointReference();
    PointReference(const PointReference&) = default;
    PointReference& operator=(const PointReference&) = default;

    Vector3<double>& position() const;
    Vector3<uint8_t>* rgbColor() const;
    Vector3<float>* normal() const;
    uint16_t* intensity() const;
    uint8_t* classification() const;
    uint8_t* edge_of_flight_line() const;
    double* gps_time() const;
    uint8_t* number_of_returns() const;
    uint8_t* return_number() const;
    uint16_t* point_source_id() const;
    uint8_t* scan_direction_flag() const;
    int8_t* scan_angle_rank() const;
    uint8_t* user_data() const;

  private:
    PointReference(PointBuffer* pointBuffer, size_t index);

    PointBuffer* _pointBuffer;
    size_t _index;
  };

  /// <summary>
  /// Constant indirect reference to a single point inside the PointBuffer
  /// </summary>
  struct PointConstReference
  {
    friend struct PointBuffer;

    PointConstReference();
    PointConstReference(const PointReference& point_reference);
    PointConstReference(const PointConstReference&) = default;
    PointConstReference& operator=(const PointConstReference&) = default;

    const Vector3<double>& position() const;
    const Vector3<uint8_t>* rgbColor() const;
    const Vector3<float>* normal() const;
    const uint16_t* intensity() const;
    const uint8_t* classification() const;
    const uint8_t* edge_of_flight_line() const;
    const double* gps_time() const;
    const uint8_t* number_of_returns() const;
    const uint8_t* return_number() const;
    const uint16_t* point_source_id() const;
    const uint8_t* scan_direction_flag() const;
    const int8_t* scan_angle_rank() const;
    const uint8_t* user_data() const;

  private:
    PointConstReference(PointBuffer const* pointBuffer, size_t index);

    PointBuffer const* _pointBuffer;
    size_t _index;
  };

  /// <summary>
  /// Creates an empty point buffer
  /// </summary>
  PointBuffer();

  /// <summary>
  /// Creates a PointBuffer from the given range of PointReferences
  /// </summary>
  explicit PointBuffer(gsl::span<PointReference> points);
  /// <summary>
  /// Creates a PointBuffer from the given range of PointConstReferences
  /// </summary>
  explicit PointBuffer(gsl::span<PointConstReference> points);

  /// <summary>
  /// Creates a new PointBuffer storing count points. For all passed attribute
  /// vectors, their count has to be equal to the specified count, otherwise a
  /// invalid_argument error is thrown
  /// </summary>
  PointBuffer(size_t count,
              std::vector<Vector3<double>> positions,
              std::vector<Vector3<uint8_t>> rgbColors = {},
              std::vector<Vector3<float>> normals = {},
              std::vector<uint16_t> intensities = {},
              std::vector<uint8_t> classifications = {},
              std::vector<uint8_t> edge_of_flight_lines = {},
              std::vector<double> gps_times = {},
              std::vector<uint8_t> number_of_returns = {},
              std::vector<uint8_t> return_numbers = {},
              std::vector<uint16_t> point_source_ids = {},
              std::vector<uint8_t> scan_direction_flags = {},
              std::vector<int8_t> scan_angle_ranks = {},
              std::vector<uint8_t> user_data = {});

  /**
   * Creates a new PointBuffer storing 'count' default-constructed points with the
   * attributes given by 'attributes'
   */
  PointBuffer(size_t count, const PointAttributes& attributes);

  PointBuffer(const PointBuffer&) = default;
  PointBuffer(PointBuffer&&) = default;

  PointBuffer& operator=(const PointBuffer&) = default;
  PointBuffer& operator=(PointBuffer&&) = default;

  /// <summary>
  /// Push a single point into this PointBuffer. If the point has attributes
  /// that are not defined in this PointBuffer they are ignored. Attributes from
  /// this PointBuffer that are not defined in the given point are filled with
  /// default values
  /// </summary>
  void push_point(PointConstReference point);

  /// <summary>
  /// Returns a constant indirect reference to the point with the given index
  /// </summary>
  PointConstReference get_point(size_t point_index) const;

  /// <summary>
  /// Returns a mutable indirect reference to the point with the given index
  /// </summary>
  PointReference get_point(size_t point_index);

  /// <summary>
  /// Appends the contents of the given PointBuffer to this PointBuffer. This
  /// will copy all attributes that exist from the given buffer into this
  /// buffer. Attributes that exist in this buffer but not in the other buffer
  /// are filled with default values
  /// </summary>
  void append_buffer(const PointBuffer& other);

  /**
   * Applies a new attribute schema to this PointBuffer. Clears all attributes that are not in the
   * new schema and fills all attributes that are in the new schema but weren't in the PointBuffers
   * old schema (i.e. the state of the PointBuffer prior to calling 'apply_schema') with default
   * values
   */
  void apply_schema(const PointAttributes& schema);

  size_t count() const { return _count; }
  bool empty() const { return _count == 0; }
  void clear();
  void shrink_to_fit();
  void resize(size_t new_size);

  std::vector<Vector3<double>>& positions() { return _positions; }
  std::vector<Vector3<uint8_t>>& rgbColors() { return _rgbColors; }
  std::vector<Vector3<float>>& normals() { return _normals; }
  std::vector<uint16_t>& intensities() { return _intensities; }
  std::vector<uint8_t>& classifications() { return _classifications; }

  auto& edge_of_flight_lines() { return _edge_of_flight_lines; }
  auto& gps_times() { return _gps_times; }
  auto& number_of_returns() { return _number_of_returns; }
  auto& return_numbers() { return _return_numbers; }
  auto& point_source_ids() { return _point_source_ids; }
  auto& scan_direction_flags() { return _scan_direction_flags; }
  auto& scan_angle_ranks() { return _scan_angle_ranks; }
  auto& user_data() { return _user_data; }

  const std::vector<Vector3<double>>& positions() const { return _positions; }
  const std::vector<Vector3<uint8_t>>& rgbColors() const { return _rgbColors; }
  const std::vector<Vector3<float>>& normals() const { return _normals; }
  const std::vector<uint16_t>& intensities() const { return _intensities; }
  const std::vector<uint8_t>& classifications() const { return _classifications; }

  const auto& edge_of_flight_lines() const { return _edge_of_flight_lines; }
  const auto& gps_times() const { return _gps_times; }
  const auto& number_of_returns() const { return _number_of_returns; }
  const auto& return_numbers() const { return _return_numbers; }
  const auto& point_source_ids() const { return _point_source_ids; }
  const auto& scan_direction_flags() const { return _scan_direction_flags; }
  const auto& scan_angle_ranks() const { return _scan_angle_ranks; }
  const auto& user_data() const { return _user_data; }

  bool hasColors() const;
  bool hasNormals() const;
  bool hasIntensities() const;
  bool hasClassifications() const;
  bool has_edge_of_flight_lines() const;
  bool has_gps_times() const;
  bool has_number_of_returns() const;
  bool has_return_numbers() const;
  bool has_point_source_ids() const;
  bool has_scan_direction_flags() const;
  bool has_scan_angle_ranks() const;
  bool has_user_data() const;

  void verify() const;

  /// <summary>
  /// Returns the raw size in bytes of the contents (positions, normals etc.) of
  /// this PointBuffer. This does NOT include the in-memory size of a
  /// PointBuffer structure itself, but rather the allocated memory of all the
  /// vectors of the PointBuffer
  /// </summary>
  size_t content_byte_size() const;

  struct PointIterator
  {
    PointIterator(PointBuffer& pointBuffer, size_t idx);

    PointReference operator*() const;
    PointIterator& operator++();
    PointIterator operator++(int);
    PointIterator& operator--();
    PointIterator operator--(int);
    PointIterator operator+(std::ptrdiff_t count) const;
    PointIterator operator-(std::ptrdiff_t count) const;
    PointIterator& operator+=(std::ptrdiff_t count);
    PointIterator& operator-=(std::ptrdiff_t count);
    friend std::ptrdiff_t operator-(const PointIterator& l, const PointIterator& r);
    PointReference operator[](std::ptrdiff_t idx) const;

    bool operator==(const PointIterator& other) const;
    bool operator!=(const PointIterator& other) const;
    friend bool operator<(const PointIterator& l, const PointIterator& r);
    friend bool operator<=(const PointIterator& l, const PointIterator& r);
    friend bool operator>(const PointIterator& l, const PointIterator& r);
    friend bool operator>=(const PointIterator& l, const PointIterator& r);

  private:
    PointBuffer* _pointBuffer;
    size_t _index;
  };

  struct PointConstIterator
  {
    PointConstIterator(PointBuffer const& pointBuffer, size_t idx);

    PointConstReference operator*() const;
    PointConstIterator& operator++();
    PointConstIterator operator++(int);
    PointConstIterator& operator--();
    PointConstIterator operator--(int);
    PointConstIterator operator+(std::ptrdiff_t count) const;
    PointConstIterator operator-(std::ptrdiff_t count) const;
    PointConstIterator& operator+=(std::ptrdiff_t count);
    PointConstIterator& operator-=(std::ptrdiff_t count);
    friend std::ptrdiff_t operator-(const PointConstIterator& l, const PointConstIterator& r);
    PointConstReference operator[](std::ptrdiff_t idx) const;

    bool operator==(const PointConstIterator& other) const;
    bool operator!=(const PointConstIterator& other) const;
    friend bool operator<(const PointConstIterator& l, const PointConstIterator& r);
    friend bool operator<=(const PointConstIterator& l, const PointConstIterator& r);
    friend bool operator>(const PointConstIterator& l, const PointConstIterator& r);
    friend bool operator>=(const PointConstIterator& l, const PointConstIterator& r);

  private:
    PointBuffer const* _pointBuffer;
    size_t _index;
  };

  PointIterator begin();
  PointIterator end();

  PointConstIterator begin() const;
  PointConstIterator end() const;

private:
  size_t _count;
  std::vector<Vector3<double>> _positions;
  std::vector<Vector3<uint8_t>> _rgbColors;
  std::vector<Vector3<float>> _normals;
  std::vector<uint16_t> _intensities;
  std::vector<uint8_t> _classifications;
  std::vector<uint8_t> _edge_of_flight_lines;
  std::vector<double> _gps_times;
  std::vector<uint8_t> _number_of_returns;
  std::vector<uint8_t> _return_numbers;
  std::vector<uint16_t> _point_source_ids;
  std::vector<uint8_t> _scan_direction_flags;
  std::vector<int8_t> _scan_angle_ranks;
  std::vector<uint8_t> _user_data;
};

namespace std {
template<>
struct iterator_traits<::PointBuffer::PointIterator>
{
  using difference_type = std::ptrdiff_t;
  using value_type = ::PointBuffer::PointReference;
  using iterator_category = std::random_access_iterator_tag;
};

template<>
struct iterator_traits<::PointBuffer::PointConstIterator>
{
  using difference_type = std::ptrdiff_t;
  using value_type = ::PointBuffer::PointConstReference;
  using iterator_category = std::random_access_iterator_tag;
};
} // namespace std

namespace concepts {

namespace detail {
template<>
constexpr bool has_constant_size_impl<::PointBuffer> = false;
}

template<>
inline unit::byte
size_in_memory(PointBuffer const& point_buffer)
{
  return (sizeof(size_t) * boost::units::information::byte) +
         size_in_memory(point_buffer.positions()) + size_in_memory(point_buffer.rgbColors()) +
         size_in_memory(point_buffer.normals()) + size_in_memory(point_buffer.intensities()) +
         size_in_memory(point_buffer.classifications());
}
} // namespace concepts

max_stars_repo_name: igd-geo/schwarzwald

max_issues_repo_name: igd-geo/schwarzwald

max_forks_repo_name: igd-geo/schwarzwald


text:
/* ode-initval/odeiv.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include "gsl_odeiv.h"

gsl_odeiv_step *
gsl_odeiv_step_alloc(const gsl_odeiv_step_type * T, size_t dim)
{
  gsl_odeiv_step *s = (gsl_odeiv_step *) malloc (sizeof (gsl_odeiv_step));

  if (s == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for ode struct", GSL_ENOMEM);
    };

  s->type = T;
  s->dimension = dim;

  s->state = s->type->alloc(dim);

  if (s->state == 0)
    {
      free (s);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_NULL ("failed to allocate space for ode state",	GSL_ENOMEM);
    };

  return s;
}

const char *
gsl_odeiv_step_name(const gsl_odeiv_step * s)
{
  return s->type->name;
}

unsigned int
gsl_odeiv_step_order(const gsl_odeiv_step * s)
{
  return s->type->order(s->state);
}

int
gsl_odeiv_step_apply(
  gsl_odeiv_step * s,
  double t,
  double h,
  double y[],
  double yerr[],
  const double dydt_in[],
  double dydt_out[],
  const gsl_odeiv_system * dydt)
{
  return s->type->apply(s->state, s->dimension, t, h, y, yerr, dydt_in, dydt_out, dydt);
}

int
gsl_odeiv_step_reset(gsl_odeiv_step * s)
{
  return s->type->reset(s->state, s->dimension);
}

void
gsl_odeiv_step_free(gsl_odeiv_step * s)
{
  s->type->free(s->state);
  free(s);
}


text:
/*
 * C version of Diffusive Nested Sampling (DNest4) by Brendon J. Brewer
 *
 * Yan-Rong Li, liyanrong@mail.ihep.ac.cn
 * Jun 30, 2016
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include <gsl/gsl_rng.h>

#include "dnestvars.h"

/* output files */
FILE *fsample, *fsample_info;

/* random number generator */
const gsl_rng_type * dnest_gsl_T;
gsl_rng * dnest_gsl_r;

Options options;
char options_file[STR_MAX_LENGTH];

// sampler
bool save_to_disk;
double compression;
unsigned int regularisation;

void *particles;
int dnest_size_of_modeltype;
int particle_offset_size, particle_offset_double;

LikelihoodType *log_likelihoods;
unsigned int *level_assignments;

// number account of unaccepted times
unsigned int *account_unaccepts;

int size_levels;
Level *levels;
unsigned int count_saves, num_saves, num_saves_restart;
int dnest_which_particle_update; // which particle to be updated
int dnest_which_level_update;    // which level to be updated;
unsigned long long int count_mcmc_steps;
LikelihoodType *above;
unsigned int size_above;

double post_logz;
int dnest_num_params;
char dnest_sample_postfix[STR_MAX_LENGTH], dnest_sample_tag[STR_MAX_LENGTH], dnest_sample_dir[STR_MAX_LENGTH];

double *limits, *copies_of_limits;

int *dnest_perturb_accept;
int dnest_root;

int dnest_flag_restart=0, dnest_flag_postprc=0, dnest_flag_sample_info=0, dnest_flag_limits=0;
double dnest_post_temp=1.0;
char file_restart[STR_MAX_LENGTH], file_save_restart[STR_MAX_LENGTH];

void *dnest_arg;

//***********************************************
/*                  functions                  */
double mod(double y, double x);
void dnest_wrap(double *x, double min, double max);
void wrap_limit(double *x, double min, double max);
int mod_int(int y, int x);
int dnest_cmp(const void *pa, const void *pb);

void options_load(int max_num_saves, double ptol);
void setup(int argc, char** argv, DNestFptrSet *fptrset, int num_params, char *sample_dir, int max_num_saves, double ptol);
void finalise();

double dnest(int argc, char **argv, DNestFptrSet *fptrset,  int num_params, char *sample_dir,
             int max_num_saves, double pdff, const void *arg);
void dnest_run();
void dnest_mcmc_run();
void update_particle(unsigned int which);
void update_level_assignment(unsigned int which);
double log_push(unsigned int which_level);
bool enough_levels(Level *l, int size_l);
void do_bookkeeping();
void save_levels();
void save_particle();
void save_limits();
void kill_lagging_particles();
void renormalise_visits();
void recalculate_log_X();
double dnest_randh();
double dnest_rand();
double dnest_randn();
int dnest_rand_int(int size);
void dnest_postprocess(double temperature, int max_num_saves, double ptol);
void postprocess(double temperature);
void initialize_output_file();
void close_output_file();
void dnest_save_restart();
void dnest_restart();
void dnest_restart_action(int iflag);
void dnest_accept_action();
void dnest_kill_action(int i, int i_copy);
void dnest_print_particle(FILE *fp, const void *model, const void *arg);
void dnest_read_particle(FILE *fp, void *model);
int dnest_get_size_levels();
int dnest_get_which_level_update();
int dnest_get_which_particle_update();
void dnest_get_posterior_sample_file(char *fname);
int dnest_check_version(char *verion_str);
unsigned int dnest_get_which_num_saves();
unsigned int dnest_get_count_saves();
unsigned long long int dnest_get_count_mcmc_steps();
void dnest_check_fptrset(DNestFptrSet *fptrset);
DNestFptrSet * dnest_malloc_fptrset();
void dnest_free_fptrset(DNestFptrSet * fptrset);
/*=====================================================*/
// users responsible for following functions
void (*print_particle)(FILE *fp, const void *model, const void *arg);
void (*read_particle)(FILE *fp, void *model);
void (*from_prior)(void *model, const void *arg);
double (*log_likelihoods_cal)(const void *model, const void *arg);
double (*log_likelihoods_cal_initial)(const void *model, const void *arg);
double (*log_likelihoods_cal_restart)(const void *model, const void *arg);
double (*perturb)(void *model, const void *arg);
void (*restart_action)(int iflag);
void (*accept_action)();
void (*kill_action)(int i, int i_copy);
/*=====================================================*/

max_stars_repo_name: LiyrAstroph/PyCALI

max_issues_repo_name: LiyrAstroph/pyCALI

max_forks_repo_name: LiyrAstroph/pyCALI