text:
/*

cp ../../../../native_ref/<target binary> libnetlib.so

gcc-mp-4.8 -O3 ddottest.c common.c -o ddottest -L. -lnetlib -I../../../../netlib/CBLAS
./ddottest  > ../../../results/mac_os_x-x86_64-ddot-CBLAS.csv

gcc-mp-4.8 -O3 ddottest.c common.c -o ddottest -I/System/Library/Frameworks/vecLib.framework/Headers -framework veclib
./ddottest  > ../../../results/mac_os_x-x86_64-ddot-veclib.csv

gcc-mp-4.8 -O3 ddottest.c common.c -o ddottest -I/opt/local/include /opt/local/lib/libatlas.a /opt/local/lib/libcblas.a /opt/local/lib/liblapack.a /opt/local/lib/libf77blas.a -lgfortran
./ddottest  > ../../../results/mac_os_x-x86_64-ddot-atlas.csv

gcc-mp-4.8 -O3 ddottest.c common.c -o ddottest -I../../../../netlib/CBLAS -L/opt/intel/composerxe/mkl/lib -lmkl_rt
export DYLD_LIBRARY_PATH=/opt/intel/composerxe/mkl/lib:/opt/intel/composerxe/lib/
./ddottest  > ../../../results/mac_os_x-x86_64-ddot-mkl.csv

gcc-mp-4.8 -O3 ddottest.c cudawrapper.c common.c -o ddottest -I../../../../netlib/CBLAS -I/usr/local/cuda/include/ -L/usr/local/cuda/lib -lcublas
export DYLD_LIBRARY_PATH=/usr/local/cuda/lib
./ddottest  > ../../../results/mac_os_x-x86_64-ddot-cuda.csv

*/

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

long benchmark(int size) {
	long requestStart, requestEnd;

    double* a = random_array(size);
    double* b = random_array(size);

	requestStart = currentTimeNanos();

    cblas_ddot(size, a, 1, b, 1);

	requestEnd = currentTimeNanos();

    free(a);
    free(b);

    return (requestEnd - requestStart);
  }

main()
{
	srand(time(NULL));

    double factor = 6.0 / 100.0;
    int i, j;
    for (i = 0 ; i < 10 ; i++) {
        for (j = 1 ; j <= 100 ; j++) {
            int size = (int) pow(10.0, factor * j);
            if (size < 10) continue;
            long took = benchmark(size);
            printf("\"%d\",\"%lu\"\n", size, took);
        }
    }
}


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

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

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*! \file allvars.h
 *  \brief declares global variables.
 *
 *  This file declares all global variables. Further variables should be added here, and declared as
 *  'extern'. The actual existence of these variables is provided by the file 'allvars.c'. To produce
 *  'allvars.c' from 'allvars.h', do the following:
 *
 *     - Erase all #define's, typedef's, and enum's
 *     - add #include "allvars.h", delete the #ifndef ALLVARS_H conditional
 *     - delete all keywords 'extern'
 *     - delete all struct definitions enclosed in {...}, e.g.
 *        "extern struct global_data_all_processes {....} All;"
 *        becomes "struct global_data_all_processes All;"
 */

#ifndef ALLVARS_H
#define ALLVARS_H

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include "tags.h"

#define  GADGETVERSION   "2.0"   /*!< code version string */

#define  TIMEBASE        (1<<28) /*!< The simulated timespan is mapped onto the integer interval [0,TIMESPAN],
                                  *   where TIMESPAN needs to be a power of 2. Note that (1<<28) corresponds to 2^29
                                  */

#define  MAXTOPNODES     200000   /*!< Maximum number of nodes in the top-level tree used for domain decomposition */

typedef  long long  peanokey;    /*!< defines the variable type used for Peano-Hilbert keys */

#define  BITS_PER_DIMENSION 18	 /*!< Bits per dimension available for Peano-Hilbert order.
				      Note: If peanokey is defined as type int, the allowed maximum is 10.
				      If 64-bit integers are used, the maximum is 21 */

#define  PEANOCELLS (((peanokey)1)<<(3*BITS_PER_DIMENSION))  /*!< The number of different Peano-Hilbert cells */

#define  RNDTABLE         3000   /*!< gives the length of a table with random numbers, refreshed at every timestep.
				      This is used to allow application of random numbers to a specific particle
				      in a way that is independent of the number of processors used. */
#define  MAX_REAL_NUMBER  1e37
#define  MIN_REAL_NUMBER  1e-37

#define  MAXLEN_FILENAME  4096    /*!< Maximum number of characters for filenames (including the full path) */

#ifdef   ISOTHERM_EQS
#define  GAMMA         (1.0)     /*!< index for isothermal gas */
#else
#define  GAMMA         (5.0/3)   /*!< adiabatic index of simulated gas */
#endif

#define  GAMMA_MINUS1  (GAMMA-1)

#define  HYDROGEN_MASSFRAC 0.76  /*!< mass fraction of hydrogen, relevant only for radiative cooling */

/* Some physical constants in cgs units */

#define  GRAVITY           6.672e-8   /*!< Gravitational constant (in cgs units) */
#define  SOLAR_MASS        1.989e33
#define  SOLAR_LUM         3.826e33
#define  RAD_CONST         7.565e-15
#define  AVOGADRO          6.0222e23
#define  BOLTZMANN         1.3806e-16
#define  GAS_CONST         8.31425e7
#define  C                 2.9979e10
#define  PLANCK            6.6262e-27
#define  CM_PER_MPC        3.085678e24
#define  PROTONMASS        1.6726e-24
#define  ELECTRONMASS      9.10953e-28
#define  THOMPSON          6.65245e-25
#define  ELECTRONCHARGE    4.8032e-10
#define  HUBBLE            3.2407789e-18	/* in h/sec */

/* Some conversion factors */

#define  SEC_PER_MEGAYEAR  3.155e13
#define  SEC_PER_YEAR      3.155e7

#ifndef ASMTH
#define ASMTH 1.25  /*!< ASMTH gives the scale of the short-range/long-range force split in units of FFT-mesh cells */
#endif

#ifndef RCUT
#define RCUT  4.5   /*!< RCUT gives the maximum distance (in units of the scale used for the force split) out to
                         which short-range forces are evaluated in the short-range tree walk. */
#endif

#define MAX_NGB             20000  /*!< defines maximum length of neighbour list */

#define MAXLEN_OUTPUTLIST   500	   /*!< maxmimum number of entries in list of snapshot output times */

#define DRIFT_TABLE_LENGTH  1000   /*!< length of the lookup table used to hold the drift and kick factors */

#define MAXITER             150    /*!< maxmimum number of steps for SPH neighbour iteration */

#ifdef DOUBLEPRECISION             /*!< If defined, the variable type FLOAT is set to "double", otherwise to FLOAT */
#define FLOAT double
#else
#define FLOAT float
#endif

#ifndef  TWODIMS
#define  NUMDIMS 3                                      /*!< For 3D-normalized kernel */
#define  KERNEL_COEFF_1  2.546479089470                 /*!< Coefficients for SPH spline kernel and its derivative */
#define  KERNEL_COEFF_2  15.278874536822
#define  KERNEL_COEFF_3  45.836623610466
#define  KERNEL_COEFF_4  30.557749073644
#define  KERNEL_COEFF_5  5.092958178941
#define  KERNEL_COEFF_6  (-15.278874536822)
#define  NORM_COEFF      4.188790204786                 /*!< Coefficient for kernel normalization. Note:  4.0/3 * PI = 4.188790204786 */
#else
#define  NUMDIMS 2                                      /*!< For 2D-normalized kernel */
#define  KERNEL_COEFF_1  (5.0/7*2.546479089470)         /*!< Coefficients for SPH spline kernel and its derivative */
#define  KERNEL_COEFF_2  (5.0/7*15.278874536822)
#define  KERNEL_COEFF_3  (5.0/7*45.836623610466)
#define  KERNEL_COEFF_4  (5.0/7*30.557749073644)
#define  KERNEL_COEFF_5  (5.0/7*5.092958178941)
#define  KERNEL_COEFF_6  (5.0/7*(-15.278874536822))
#define  NORM_COEFF      M_PI                           /*!< Coefficient for kernel normalization. */
#endif

extern int ThisTask;		/*!< the rank of the local processor */
extern int NTask;               /*!< number of processors */
extern int PTask;	        /*!< smallest integer such that NTask <= 2^PTask */

extern int NumPart;		/*!< number of particles on the LOCAL processor */
extern int N_gas;		/*!< number of gas particles on the LOCAL processor  */
extern long long Ntype[6];      /*!< total number of particles of each type */
extern int NtypeLocal[6];       /*!< local number of particles of each type */

extern int NumForceUpdate;      /*!< number of active particles on local processor in current timestep  */
extern int NumSphUpdate;        /*!< number of active SPH particles on local processor in current timestep  */

extern double CPUThisRun;	/*!< Sums the CPU time for the process (current submission only) */

extern int RestartFlag;         /*!< taken from command line used to start code. 0 is normal start-up from
                                     initial conditions, 1 is resuming a run from a set of restart files, while 2
                                     marks a restart from a snapshot file. */

extern char *Exportflag;        /*!< Buffer used for flagging whether a particle needs to be exported to another process */

extern int  *Ngblist;           /*!< Buffer to hold indices of neighbours retrieved by the neighbour search routines */

extern int TreeReconstructFlag; /*!< Signals that a new tree needs to be constructed */

extern int Flag_FullStep;       /*!< This flag signals that the current step involves all particles */

extern int ZeroTimestepEncountered;  /*!< Flag used by AMUSE. When a particle is assigned a timestep of zero, an
                                     exception is raised instead of forcing the application to exit. */

extern gsl_rng *random_generator; /*!< the employed random number generator of the GSL library */

extern double RndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */

extern double DomainCorner[3];    /*!< gives the lower left corner of simulation volume */
extern double DomainCenter[3];    /*!< gives the center of simulation volume */
extern double DomainLen;          /*!< gives the (maximum) side-length of simulation volume */
extern double DomainFac;          /*!< factor used for converting particle coordinates to a Peano-Hilbert mesh covering the simulation volume */
extern int    DomainMyStart;      /*!< first domain mesh cell that resides on the local processor */
extern int    DomainMyLast;       /*!< last domain mesh cell that resides on the local processor */
extern int    *DomainStartList;   /*!< a table that lists the first domain mesh cell for all processors */
extern int    *DomainEndList;     /*!< a table that lists the last domain mesh cell for all processors */
extern double *DomainWork;        /*!< a table that gives the total "work" due to the particles stored by each processor */
extern int    *DomainCount;       /*!< a table that gives the total number of particles held by each processor */
extern int    *DomainCountSph;    /*!< a table that gives the total number of SPH particles held by each processor */

extern int    *DomainTask;        /*!< this table gives for each leaf of the top-level tree the processor it was assigned to */
extern int    *DomainNodeIndex;   /*!< this table gives for each leaf of the top-level tree the corresponding node of the gravitational tree */
extern FLOAT  *DomainTreeNodeLen; /*!< this table gives for each leaf of the top-level tree the side-length of the corresponding node of the gravitational tree */
extern FLOAT  *DomainHmax;        /*!< this table gives for each leaf of the top-level tree the maximum SPH smoothing length among the particles of the corresponding node of the gravitational tree */

extern struct DomainNODE
{
  FLOAT s[3];                     /*!< center-of-mass coordinates */
  FLOAT vs[3];                    /*!< center-of-mass velocities */
  FLOAT mass;                     /*!< mass of node */
#ifdef UNEQUALSOFTENINGS
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
  int   bitflags;                 /*!< this bit-field encodes the particle type with the largest softening among the particles of the nodes, and whether there are particles with different softening in the node */
#else
  FLOAT maxsoft;                  /*!< hold the maximum gravitational softening of particles in the
                                       node if the ADAPTIVE_GRAVSOFT_FORGAS option is selected */
#endif
#endif
}
 *DomainMoment;                   /*!< this table stores for each node of the top-level tree corresponding node data from the gravitational tree */

extern peanokey *DomainKeyBuf;    /*!< this points to a buffer used during the exchange of particle data */

extern peanokey *Key;             /*!< a table used for storing Peano-Hilbert keys for particles */
extern peanokey *KeySorted;       /*!< holds a sorted table of Peano-Hilbert keys for all particles, used to construct top-level tree */

extern int NTopnodes;             /*!< total number of nodes in top-level tree */
extern int NTopleaves;            /*!< number of leaves in top-level tree. Each leaf can be assigned to a different processor */

extern struct topnode_data
{
  int Daughter;                   /*!< index of first daughter cell (out of 8) of top-level node */
  int Pstart;                     /*!< for the present top-level node, this gives the index of the first node in the concatenated list of topnodes collected from all processors */
  int Blocks;                     /*!< for the present top-level node, this gives the number of corresponding nodes in the concatenated list of topnodes collected from all processors */
  int Leaf;                       /*!< if the node is a leaf, this gives its number when all leaves are traversed in Peano-Hilbert order */
  peanokey Size;                  /*!< number of Peano-Hilbert mesh-cells represented by top-level node */
  peanokey StartKey;              /*!< first Peano-Hilbert key in top-level node */
  long long Count;                /*!< counts the number of particles in this top-level node */
}
 *TopNodes;                       /*!< points to the root node of the top-level tree */

extern double TimeOfLastTreeConstruction; /*!< holds what it says, only used in connection with FORCETEST */

/* variables for input/output, usually only used on process 0 */

extern char ParameterFile[MAXLEN_FILENAME];  /*!< file name of parameterfile used for starting the simulation */

extern FILE *FdInfo;       /*!< file handle for info.txt log-file. */
extern FILE *FdEnergy;     /*!< file handle for energy.txt log-file. */
extern FILE *FdTimings;    /*!< file handle for timings.txt log-file. */
extern FILE *FdCPU;        /*!< file handle for cpu.txt log-file. */

#ifdef FORCETEST
extern FILE *FdForceTest;  /*!< file handle for forcetest.txt log-file. */
#endif

extern double DriftTable[DRIFT_TABLE_LENGTH];      /*!< table for the cosmological drift factors */
extern double GravKickTable[DRIFT_TABLE_LENGTH];   /*!< table for the cosmological kick factor for gravitational forces */
extern double HydroKickTable[DRIFT_TABLE_LENGTH];  /*!< table for the cosmological kick factor for hydrodynmical forces */

extern void *CommBuffer;   /*!< points to communication buffer, which is used in the domain decomposition, the
                                parallel tree-force computation, the SPH routines, etc. */

/*! This structure contains data which is the SAME for all tasks (mostly code parameters read from the
 * parameter file).  Holding this data in a structure is convenient for writing/reading the restart file, and
 * it allows the introduction of new global variables in a simple way. The only thing to do is to introduce
 * them into this structure.
 */
extern struct global_data_all_processes
{
  long long TotNumPart;		/*!< total particle numbers (global value) */
  long long TotN_gas;		/*!< total gas particle number (global value) */

  int MaxPart;			/*!< This gives the maxmimum number of particles that can be stored on one processor. */
  int MaxPartSph;		/*!< This gives the maxmimum number of SPH particles that can be stored on one processor. */

  double BoxSize;               /*!< Boxsize in case periodic boundary conditions are used */

  int ICFormat;			/*!< selects different versions of IC file-format */

  int SnapFormat;		/*!< selects different versions of snapshot file-formats */

  int NumFilesPerSnapshot;      /*!< number of files in multi-file snapshot dumps */
  int NumFilesWrittenInParallel;/*!< maximum number of files that may be written simultaneously when
                                     writing/reading restart-files, or when writing snapshot files */

  int BufferSize;		/*!< size of communication buffer in MB */
  int BunchSizeForce;		/*!< number of particles fitting into the buffer in the parallel tree-force algorithm  */
  int BunchSizeDensity;         /*!< number of particles fitting into the communication buffer in the density computation */
  int BunchSizeHydro;           /*!< number of particles fitting into the communication buffer in the SPH hydrodynamical force computation */
  int BunchSizeDomain;          /*!< number of particles fitting into the communication buffer in the domain decomposition */

  double PartAllocFactor;	/*!< in order to maintain work-load balance, the particle load will usually
				     NOT be balanced.  Each processor allocates memory for PartAllocFactor times
				     the average number of particles to allow for that */

  double TreeAllocFactor;	/*!< Each processor allocates a number of nodes which is TreeAllocFactor times
				     the maximum(!) number of particles.  Note: A typical local tree for N
				     particles needs usually about ~0.65*N nodes. */

  /* some SPH parameters */

  double DesNumNgb;             /*!< Desired number of SPH neighbours */
  double MaxNumNgbDeviation;    /*!< Maximum allowed deviation neighbour number */

  double ArtBulkViscConst;      /*!< Sets the parameter \f$\alpha\f$ of the artificial viscosity */
  double ArtBulkViscBeta;      /*!< Sets the parameter \f$\beta\f$ of the artificial viscosity */
  double InitGasTemp;		/*!< may be used to set the temperature in the IC's */
  double MinGasTemp;		/*!< may be used to set a floor for the gas temperature */
  double MinEgySpec;            /*!< the minimum allowed temperature expressed as energy per unit mass */

  /* some force counters  */

  long long TotNumOfForces;	             /*!< counts total number of force computations  */
  long long NumForcesSinceLastDomainDecomp;  /*!< count particle updates since last domain decomposition */

  /* system of units  */

  double G;                        /*!< Gravity-constant in internal units */
  double UnitTime_in_s;   	   /*!< factor to convert internal time unit to seconds/h */
  double UnitMass_in_g;            /*!< factor to convert internal mass unit to grams/h */
  double UnitVelocity_in_cm_per_s; /*!< factor to convert intqernal velocity unit to cm/sec */
  double UnitLength_in_cm;         /*!< factor to convert internal length unit to cm/h */
  double UnitPressure_in_cgs;      /*!< factor to convert internal pressure unit to cgs units (little 'h' still around!) */
  double UnitDensity_in_cgs;       /*!< factor to convert internal length unit to g/cm^3*h^2 */
  double UnitCoolingRate_in_cgs;   /*!< factor to convert internal cooling rate to cgs units */
  double UnitEnergy_in_cgs;        /*!< factor to convert internal energy to cgs units */
  double UnitTime_in_Megayears;    /*!< factor to convert internal time to megayears/h */
  double GravityConstantInternal;  /*!< If set to zero in the parameterfile, the internal value of the
                                        gravitational constant is set to the Newtonian value based on the system of
                                        units specified. Otherwise the value provided is taken as internal gravity constant G. */

  /* Cosmological parameters */

  double Hubble;       /*!< Hubble-constant in internal units */
  double Omega0;       /*!< matter density in units of the critical density (at z=0)*/
  double OmegaLambda;  /*!< vaccum energy density relative to crictical density (at z=0) */
  double OmegaBaryon;  /*!< baryon density in units of the critical density (at z=0)*/
  double HubbleParam;  /*!< little `h', i.e. Hubble constant in units of 100 km/s/Mpc.  Only needed to get absolute physical values for cooling physics */

  /* Code options */

  int ComovingIntegrationOn;	/*!< flags that comoving integration is enabled */
  int PeriodicBoundariesOn;     /*!< flags that periodic boundaries are enabled */
  int ResubmitOn;               /*!< flags that automatic resubmission of job to queue system is enabled */
  int TypeOfOpeningCriterion;   /*!< determines tree cell-opening criterion: 0 for Barnes-Hut, 1 for relative criterion */
  int TypeOfTimestepCriterion;  /*!< gives type of timestep criterion (only 0 supported right now - unlike gadget-1.1) */
  int OutputListOn;             /*!< flags that output times are listed in a specified file */

  /* Parameters determining output frequency */

  int SnapshotFileCount;        /*!< number of snapshot that is written next */
  double TimeBetSnapshot;       /*!< simulation time interval between snapshot files */
  double TimeOfFirstSnapshot;   /*!< simulation time of first snapshot files */
  double CpuTimeBetRestartFile; /*!< cpu-time between regularly generated restart files */
  double TimeLastRestartFile;   /*!< cpu-time when last restart-file was written */
  double TimeBetStatistics;     /*!< simulation time interval between computations of energy statistics */
  double TimeLastStatistics;    /*!< simulation time when the energy statistics was computed the last time */
  int NumCurrentTiStep;         /*!< counts the number of system steps taken up to this point */

  /* Current time of the simulation, global step, and end of simulation */

  double Time;                  /*!< current time of the simulation */
  double TimeBegin;             /*!< time of initial conditions of the simulation */
  double TimeStep;              /*!< difference between current times of previous and current timestep */
  double TimeMax;	        /*!< marks the point of time until the simulation is to be evolved */

  /* variables for organizing discrete timeline */

  double Timebase_interval;     /*!< factor to convert from floating point time interval to integer timeline */
  int Ti_Current;               /*!< current time on integer timeline */
  int Ti_nextoutput;            /*!< next output time on integer timeline */
#ifdef FLEXSTEPS
  int PresentMinStep;           /*!< If FLEXSTEPS is used, particle timesteps are chosen as multiples of the present minimum timestep. */
  int PresentMaxStep;		/*!< If FLEXSTEPS is used, this is the maximum timestep in timeline units, rounded down to the next power 2 division */
#endif
#ifdef PMGRID
  int PM_Ti_endstep;            /*!< begin of present long-range timestep */
  int PM_Ti_begstep;            /*!< end of present long-range timestep */
#endif

  /* Placement of PM grids */

#ifdef PMGRID
  double Asmth[2];              /*!< Gives the scale of the long-range/short-range split (in mesh-cells), both for the coarse and the high-res mesh */
  double Rcut[2];               /*!< Gives the maximum radius for which the short-range force is evaluated with the tree (in mesh-cells), both for the coarse and the high-res mesh */
  double Corner[2][3];          /*!< lower left corner of coarse and high-res PM-mesh */
  double UpperCorner[2][3];     /*!< upper right corner of coarse and high-res PM-mesh */
  double Xmintot[2][3];         /*!< minimum particle coordinates both for coarse and high-res PM-mesh */
  double Xmaxtot[2][3];         /*!< maximum particle coordinates both for coarse and high-res PM-mesh */
  double TotalMeshSize[2];      /*!< total extension of coarse and high-res PM-mesh */
#endif

  /* Variables that keep track of cumulative CPU consumption */

  double TimeLimitCPU;          /*!< CPU time limit as defined in parameterfile */
  double CPU_TreeConstruction;  /*!< time spent for constructing the gravitational tree */
  double CPU_TreeWalk;          /*!< actual time spent for pure tree-walks */
  double CPU_Gravity;           /*!< cumulative time used for gravity computation (tree-algorithm only) */
  double CPU_Potential;         /*!< time used for computing gravitational potentials */
  double CPU_Domain;            /*!< cumulative time spent for domain decomposition */
  double CPU_Snapshot;          /*!< time used for writing snapshot files */
  double CPU_Total;             /*!< cumulative time spent for domain decomposition */
  double CPU_CommSum;           /*!< accumulated time used for communication, and for collecting partial results, in tree-gravity */
  double CPU_Imbalance;         /*!< cumulative time lost accross all processors as work-load imbalance in gravitational tree */
  double CPU_HydCompWalk;       /*!< time used for actual SPH computations, including neighbour search */
  double CPU_HydCommSumm;       /*!< cumulative time used for communication in SPH, and for collecting partial results */
  double CPU_HydImbalance;      /*!< cumulative time lost due to work-load imbalance in SPH */
  double CPU_Hydro;             /*!< cumulative time spent for SPH related computations */
  double CPU_EnsureNgb;         /*!< time needed to iterate on correct neighbour numbers */
  double CPU_Predict;           /*!< cumulative time to drift the system forward in time, including dynamic tree updates */
  double CPU_TimeLine;          /*!< time used for determining new timesteps, and for organizing the timestepping, including kicks of active particles */
  double CPU_PM;                /*!< time used for long-range gravitational force */
  double CPU_Peano;             /*!< time required to establish Peano-Hilbert order */

  /* tree code opening criterion */

  double ErrTolTheta;		/*!< BH tree opening angle */
  double ErrTolForceAcc;	/*!< parameter for relative opening criterion in tree walk */

  /* adjusts accuracy of time-integration */

  double ErrTolIntAccuracy;	/*!< accuracy tolerance parameter \f$ \eta \f$ for timestep criterion. The
                                     timestep is \f$ \Delta t = \sqrt{\frac{2 \eta eps}{a}} \f$ */

  double MinSizeTimestep;       /*!< minimum allowed timestep. Normally, the simulation terminates if the
                                     timestep determined by the timestep criteria falls below this limit. */
  double MaxSizeTimestep;       /*!< maximum allowed timestep */

  double MaxRMSDisplacementFac; /*!< this determines a global timestep criterion for cosmological simulations
                                     in comoving coordinates.  To this end, the code computes the rms velocity
                                     of all particles, and limits the timestep such that the rms displacement
                                     is a fraction of the mean particle separation (determined from the
                                     particle mass and the cosmological parameters). This parameter specifies
                                     this fraction. */

  double CourantFac;		/*!< SPH-Courant factor */

  int BunchSizeTime;		/*!< number of particles fitting into the communication buffer in the timestep communication */

  /* frequency of tree reconstruction/domain decomposition */

  double TreeDomainUpdateFrequency; /*!< controls frequency of domain decompositions  */

  /* Gravitational and hydrodynamical softening lengths (given in terms of an `equivalent' Plummer softening length).
   * Five groups of particles are supported 0="gas", 1="halo", 2="disk", 3="bulge", 4="stars", 5="bndry"
   */

  double MinGasHsmlFractional;  /*!< minimum allowed SPH smoothing length in units of SPH gravitational softening length */
  double MinGasHsml;            /*!< minimum allowed SPH smoothing length */

  double SofteningGas;          /*!< comoving gravitational softening lengths for type 0 */
  double SofteningHalo;         /*!< comoving gravitational softening lengths for type 1 */
  double SofteningDisk;         /*!< comoving gravitational softening lengths for type 2 */
  double SofteningBulge;        /*!< comoving gravitational softening lengths for type 3 */
  double SofteningStars;        /*!< comoving gravitational softening lengths for type 4 */
  double SofteningBndry;        /*!< comoving gravitational softening lengths for type 5 */

  double SofteningGasMaxPhys;   /*!< maximum physical softening length for type 0 */
  double SofteningHaloMaxPhys;  /*!< maximum physical softening length for type 1 */
  double SofteningDiskMaxPhys;  /*!< maximum physical softening length for type 2 */
  double SofteningBulgeMaxPhys; /*!< maximum physical softening length for type 3 */
  double SofteningStarsMaxPhys; /*!< maximum physical softening length for type 4 */
  double SofteningBndryMaxPhys; /*!< maximum physical softening length for type 5 */

  double SofteningTable[6];     /*!< current (comoving) gravitational softening lengths for each particle type */
  double ForceSoftening[6];     /*!< the same, but multiplied by a factor 2.8 - at that scale the force is Newtonian */

  double MassTable[6];          /*!< Table with particle masses for particle types with equal mass.
                                     If particle masses are all equal for one type, the corresponding entry in MassTable
                                     is set to this value, allowing the size of the snapshot files to be reduced. */

  /* some filenames */

  char InitCondFile[MAXLEN_FILENAME];          /*!< filename of initial conditions */
  char OutputDir[MAXLEN_FILENAME];             /*!< output directory of the code */
  char SnapshotFileBase[MAXLEN_FILENAME];      /*!< basename to construct the names of snapshotf files */
  char EnergyFile[MAXLEN_FILENAME];            /*!< name of file with energy statistics */
  char CpuFile[MAXLEN_FILENAME];               /*!< name of file with cpu-time statistics */
  char InfoFile[MAXLEN_FILENAME];              /*!< name of log-file with a list of the timesteps taken */
  char TimingsFile[MAXLEN_FILENAME];           /*!< name of file with performance metrics of gravitational tree algorithm */
  char RestartFile[MAXLEN_FILENAME];           /*!< basename of restart-files */
  char ResubmitCommand[MAXLEN_FILENAME];       /*!< name of script-file that will be executed for automatic restart */
  char OutputListFilename[MAXLEN_FILENAME];    /*!< name of file with list of desired output times */

  double OutputListTimes[MAXLEN_OUTPUTLIST];   /*!< table with desired output times */
  int OutputListLength;                        /*!< number of output times stored in the table of desired output times */

}
 All;                                          /*!< a container variable for global variables that are equal on all processors */

/*! This structure holds all the information that is
 * stored for each particle of the simulation.
 */
extern struct particle_data
{
  FLOAT Pos[3];			/*!< particle position at its current time */
  FLOAT Mass;			/*!< particle mass */
  FLOAT Vel[3];			/*!< particle velocity at its current time */
  FLOAT GravAccel[3];		/*!< particle acceleration due to gravity */
#ifdef PMGRID
  FLOAT GravPM[3];		/*!< particle acceleration due to long-range PM gravity force*/
#endif
#ifdef FORCETEST
  FLOAT GravAccelDirect[3];	/*!< particle acceleration when computed with direct summation */
#endif
  FLOAT Potential;		/*!< gravitational potential */
  FLOAT OldAcc;			/*!< magnitude of old gravitational force. Used in relative opening criterion */
#ifndef LONGIDS
  unsigned int ID;		/*!< particle identifier */
#else
  unsigned long long ID;        /*!< particle identifier */
#endif

  int Type;		        /*!< flags particle type.  0=gas, 1=halo, 2=disk, 3=bulge, 4=stars, 5=bndry */
  int Ti_endstep;               /*!< marks start of current timestep of particle on integer timeline */
  int Ti_begstep;               /*!< marks end of current timestep of particle on integer timeline */
#ifdef TIMESTEP_LIMITER
  int Ti_sizestep;
#endif
#ifdef FLEXSTEPS
  int FlexStepGrp;		/*!< a random 'offset' on the timeline to create a smooth groouping of particles */
#endif
  float GravCost;		/*!< weight factor used for balancing the work-load */
#ifdef PSEUDOSYMMETRIC
  float AphysOld;               /*!< magnitude of acceleration in last timestep. Used to make a first order
                                     prediction of the change of acceleration expected in the future, thereby
                                     allowing to guess whether a decrease/increase of the timestep should occur
                                     in the timestep that is started. */
#endif
}
 *P,              /*!< holds particle data on local processor */
 *DomainPartBuf;  /*!< buffer for particle data used in domain decomposition */

/* the following struture holds data that is stored for each SPH particle in addition to the collisionless
 * variables.
 */
extern struct sph_particle_data
{
  FLOAT Entropy;                /*!< current value of entropy (actually entropic function) of particle */
  FLOAT Density;		/*!< current baryonic mass density of particle */
  FLOAT Hsml;			/*!< current smoothing length */
  FLOAT Left;                   /*!< lower bound in iterative smoothing length search */
  FLOAT Right;                  /*!< upper bound in iterative smoothing length search */
  FLOAT NumNgb;                 /*!< weighted number of neighbours found */
  FLOAT Pressure;		/*!< current pressure */
  FLOAT DtEntropy;              /*!< rate of change of entropy */
  FLOAT HydroAccel[3];		/*!< acceleration due to hydrodynamical force */
  FLOAT VelPred[3];		/*!< predicted SPH particle velocity at the current time */
  FLOAT DivVel;			/*!< local velocity divergence */
  FLOAT CurlVel;		/*!< local velocity curl */
  FLOAT Rot[3];		        /*!< local velocity curl */
  FLOAT DhsmlDensityFactor;     /*!< correction factor needed in the equation of motion of the conservative entropy formulation of SPH */
  FLOAT MaxSignalVel;           /*!< maximum "signal velocity" occuring for this particle */
#ifdef TIMESTEP_UPDATE
  int   FeedbackFlag;
  FLOAT FeedAccel[3];  /*!< acceleration due to feedback force */
#endif
#ifdef MORRIS97VISC
  FLOAT Alpha;		        /*!< viscosity coefficient */
  FLOAT DAlphaDt;       		/*!< time rate of change of viscosity coefficient */
#endif
}
 *SphP,                        	/*!< holds SPH particle data on local processor */
 *DomainSphBuf;                 /*!< buffer for SPH particle data in domain decomposition */

/*  Variables for Tree
 */

extern int MaxNodes;		/*!< maximum allowed number of internal nodes */
extern int Numnodestree;	/*!< number of (internal) nodes in each tree */

extern struct NODE
{
  FLOAT len;			/*!< sidelength of treenode */
  FLOAT center[3];		/*!< geometrical center of node */
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
  FLOAT maxsoft;                /*!< hold the maximum gravitational softening of particles in the
                                     node if the ADAPTIVE_GRAVSOFT_FORGAS option is selected */
#endif
  union
  {
    int suns[8];		/*!< temporary pointers to daughter nodes */
    struct
    {
      FLOAT s[3];               /*!< center of mass of node */
      FLOAT mass;               /*!< mass of node */
      int bitflags;             /*!< a bit-field with various information on the node */
      int sibling;              /*!< this gives the next node in the walk in case the current node can be used */
      int nextnode;             /*!< this gives the next node in case the current node needs to be opened */
      int father;               /*!< this gives the parent node of each node (or -1 if we have the root node) */
    }
    d;
  }
  u;
}
 *Nodes_base,                   /*!< points to the actual memory allocted for the nodes */
 *Nodes;                        /*!< this is a pointer used to access the nodes which is shifted such that Nodes[All.MaxPart]
 				     gives the first allocated node */

extern int *Nextnode;	        /*!< gives next node in tree walk */
extern int *Father;	        /*!< gives parent node in tree    */

extern struct extNODE           /*!< this structure holds additional tree-node information which is not needed in the actual gravity computation */
{
  FLOAT hmax;			/*!< maximum SPH smoothing length in node. Only used for gas particles */
  FLOAT vs[3];			/*!< center-of-mass velocity */
}
 *Extnodes_base,                /*!< points to the actual memory allocted for the extended node information */
 *Extnodes;                     /*!< provides shifted access to extended node information, parallel to Nodes/Nodes_base */

/*! Header for the standard file format.
 */
extern struct io_header
{
  int npart[6];                        /*!< number of particles of each type in this file */
  double mass[6];                      /*!< mass of particles of each type. If 0, then the masses are explicitly
                                            stored in the mass-block of the snapshot file, otherwise they are omitted */
  double time;                         /*!< time of snapshot file */
  double redshift;                     /*!< redshift of snapshot file */
  int flag_sfr;                        /*!< flags whether the simulation was including star formation */
  int flag_feedback;                   /*!< flags whether feedback was included (obsolete) */
  unsigned int npartTotal[6];          /*!< total number of particles of each type in this snapshot. This can be
                                            different from npart if one is dealing with a multi-file snapshot. */
  int flag_cooling;                    /*!< flags whether cooling was included  */
  int num_files;                       /*!< number of files in multi-file snapshot */
  double BoxSize;                      /*!< box-size of simulation in case periodic boundaries were used */
  double Omega0;                       /*!< matter density in units of critical density */
  double OmegaLambda;                  /*!< cosmological constant parameter */
  double HubbleParam;                  /*!< Hubble parameter in units of 100 km/sec/Mpc */
  int flag_stellarage;                 /*!< flags whether the file contains formation times of star particles */
  int flag_metals;                     /*!< flags whether the file contains metallicity values for gas and star particles */
  unsigned int npartTotalHighWord[6];  /*!< High word of the total number of particles of each type */
  int  flag_entropy_instead_u;         /*!< flags that IC-file contains entropy instead of u */
  char fill[60];	               /*!< fills to 256 Bytes */
}
 header;                               /*!< holds header for snapshot files */

#define IO_NBLOCKS 11   /*!< total number of defined information blocks for snapshot files.
                             Must be equal to the number of entries in "enum iofields" */

enum iofields           /*!< this enumeration lists the defined output blocks in snapshot files. Not all of them need to be present. */
{
  IO_POS,
  IO_VEL,
  IO_ID,
  IO_MASS,
  IO_U,
  IO_RHO,
  IO_HSML,
  IO_POT,
  IO_ACCEL,
  IO_DTENTR,
  IO_TSTP,
};

extern char Tab_IO_Labels[IO_NBLOCKS][4];   /*<! This table holds four-byte character tags used for fileformat 2 */

/* global state of system, used for global statistics
 */
extern struct state_of_system
{
  double Mass;
  double EnergyKin;
  double EnergyPot;
  double EnergyInt;
  double EnergyTot;
  double Momentum[4];
  double AngMomentum[4];
  double CenterOfMass[4];
  double MassComp[6];
  double EnergyKinComp[6];
  double EnergyPotComp[6];
  double EnergyIntComp[6];
  double EnergyTotComp[6];
  double MomentumComp[6][4];
  double AngMomentumComp[6][4];
  double CenterOfMassComp[6][4];
}
 SysState;                       /*<! Structure for storing some global statistics about the simulation. */

/* Various structures for communication
 */
extern struct gravdata_in
{
  union
  {
    FLOAT Pos[3];
    FLOAT Acc[3];
    FLOAT Potential;
  }
  u;
#ifdef UNEQUALSOFTENINGS
  int Type;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
  FLOAT Soft;
#endif
#endif
  union
  {
    FLOAT OldAcc;
    int Ninteractions;
  }
  w;
}
 *GravDataIn,                   /*!< holds particle data to be exported to other processors */
 *GravDataGet,                  /*!< holds particle data imported from other processors */
 *GravDataResult,               /*!< holds the partial results computed for imported particles. Note: We use GravDataResult = GravDataGet, such that the result replaces the imported data */
 *GravDataOut;                  /*!< holds partial results received from other processors. This will overwrite the GravDataIn array */

extern struct gravdata_index
{
  int Task;
  int Index;
  int SortIndex;
}
 *GravDataIndexTable;           /*!< the particles to be exported are grouped by task-number. This table allows the results to be disentangled again and to be assigned to the correct particle */

extern struct densdata_in
{
  FLOAT Pos[3];
  FLOAT Vel[3];
  FLOAT Hsml;
  int Index;
  int Task;
}
 *DensDataIn,                   /*!< holds particle data for SPH density computation to be exported to other processors */
 *DensDataGet;                  /*!< holds imported particle data for SPH density computation */

extern struct densdata_out
{
  FLOAT Rho;
  FLOAT Div, Rot[3];
  FLOAT DhsmlDensity;
  FLOAT Ngb;
}
 *DensDataResult,               /*!< stores the locally computed SPH density results for imported particles */
 *DensDataPartialResult;        /*!< imported partial SPH density results from other processors */

extern struct hydrodata_in
{
  FLOAT Pos[3];
  FLOAT Vel[3];
  FLOAT Hsml;
  FLOAT Mass;
  FLOAT Density;
  FLOAT Pressure;
  FLOAT F1;
  FLOAT DhsmlDensityFactor;
  int   Timestep;
  int   Task;
  int   Index;
#ifdef MORRIS97VISC
  FLOAT Alpha;
#endif
}
 *HydroDataIn,                  /*!< holds particle data for SPH hydro-force computation to be exported to other processors */
 *HydroDataGet;                 /*!< holds imported particle data for SPH hydro-force computation */

extern struct hydrodata_out
{
  FLOAT Acc[3];
  FLOAT DtEntropy;
  FLOAT MaxSignalVel;
}
 *HydroDataResult,              /*!< stores the locally computed SPH hydro results for imported particles */
 *HydroDataPartialResult;       /*!< imported partial SPH hydro-force results from other processors */

#ifdef TIMESTEP_LIMITER
  extern struct timedata_in
  {
    FLOAT  Pos[3];
    FLOAT  Hsml;
    int    Size;
    int    Begin;
    int    Index;
    int    Task;
  }
  *TimeDataIn,
  *TimeDataGet;
#endif

#ifndef NOMPI
#ifdef MPI_VERSION
    extern MPI_Comm GADGET_WORLD;
#endif
#endif

#endif

ext: h

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
/*
 * Copyright (c) 2016 Kengo Tateishi (@tkengo)
 * Licensed under MIT license.
 *   http://www.opensource.org/licenses/mit-license.php
 *
 * MxRuby is really inspired by numpy.
 *
 * This file is implementation of extended array functions for ruby. The MX class is defined on ruby
 * interpreter by compiling it.
 */
#include <ruby.h>
/* #include <cblas.h> */
#include "mx.h"

/*
 * MX class object in C. This is initialized in Init_mxruby function.
 */
VALUE rb_cMx;
VALUE rb_eDataTypeError;

MX *mxx_initialize(VALUE shape, DTYPE dtype)
{
    MX *mx = MX_ALLOC(MX);
    mxx_setup(mx, dtype);
    mxx_initialize_shape(mx, shape);

    if (mx->size > 0) {
        mx->elptr = MX_ALLOC_N(char, DTYPE_SIZES[mx->dtype] * mx->size);
    }

    return mx;
}

void mxx_setup(MX *mx, DTYPE dtype)
{
    mx->dim = 0;
    mx->size = 0;
    mx->dtype = dtype;
    mx->shape = NULL;
    mx->elptr = NULL;
}

void mxx_initialize_shape(MX *mx, VALUE shape)
{
    char shape_is_array = TYPE(shape) == T_ARRAY;

    mx->size = 1;
    mx->dim = shape_is_array ? RARRAY_LEN(shape) : 1;
    mx->shape = MX_ALLOC_N(size_t, mx->dim);

    for (int i = 0; i < mx->dim; i++) {
        mx->shape[i] = FIX2INT(shape_is_array ? RARRAY_AREF(shape, i) : shape);
        mx->size *= mx->shape[i];
    }
}

void mxx_free(MX *mx)
{
    if (mx->shape != NULL) {
        MX_FREE(mx->shape);
    }
    if (mx->elptr != NULL) {
        MX_FREE(mx->elptr);
    }

    MX_FREE(mx);
}

static VALUE mx_alloc(VALUE klass)
{
    MX *mx = MX_ALLOC(MX);
    mxx_setup(mx, DTYPE_FLOAT64);

    return Data_Wrap_Struct(klass, 0, mxx_free, mx);
}

static VALUE mx_initialize(VALUE self, VALUE shape, VALUE initial_value, VALUE opt)
{
    MX *mx;
    Data_Get_Struct(self, MX, mx);

    if (!NIL_P(opt)) {
        VALUE rb_dtype = rb_hash_aref(opt, ID2SYM(rb_intern("dtype")));
        if (!NIL_P(rb_dtype)) {
            mx->dtype = mxx_dtype_from_symbol(rb_dtype);
        }
    }

    mxx_initialize_shape(mx, shape);

    if (mx->size > 0) {
        mx->elptr = MX_ALLOC_N(char, DTYPE_SIZES[mx->dtype] * mx->size);

        if (TYPE(initial_value) == T_ARRAY) {
            for (int i = 0; i < mx->size; i++) {
                VALUE v = RARRAY_AREF(initial_value, i);
                mxx_rb_to_c(mx->elptr + i * DTYPE_SIZES[mx->dtype], mx->dtype, v);
            }
        } else if (!NIL_P(initial_value)) {
            for (int i = 0; i < mx->size; i++) {
                mxx_rb_to_c(mx->elptr + i * DTYPE_SIZES[mx->dtype], mx->dtype, initial_value);
            }
        }
    }

    return self;
}

static VALUE mx_shape(VALUE self)
{
    MX *mx = MX_DATA_PTR(self);
    VALUE *shape = MX_ALLOC_N(VALUE, mx->dim);
    for (int i = 0; i < mx->dim; i++) {
        shape[i] = INT2FIX(mx->shape[i]);
    }

    return rb_ary_new4(mx->dim, shape);
}

static VALUE mx_dim(VALUE self)
{
    return INT2FIX(MX_DATA_DIM(self));
}

static VALUE mx_size(VALUE self)
{
    return INT2FIX(MX_DATA_SIZE(self));
}

static VALUE mx_dtype(VALUE self)
{
    return ID2SYM(rb_intern(DTYPE_NAMES[MX_DATA_PTR(self)->dtype]));
}

static VALUE mx_dot(VALUE self, VALUE target)
{
    int N = 2;
    MX *mx = MX_DATA_PTR(self);
    void *c = MX_ALLOC_N(char, mx->shape[0] * mx->shape[1] * DTYPE_SIZES[mx->dtype]);

    /* cblas_sgemm( */
    /*         CblasRowMajor, CblasNoTrans, CblasNoTrans, */
    /*         N, N, N, */
    /*         1, mx->elptr, N, */
    /*         MX_DATA_ELPTR(target), N, */
    /*         0, c, N */
    /*         ); */

    MX_FREE(mx->elptr);
    mx->elptr = c;
    return self;
}

static VALUE mx_astype(int argc, VALUE *argv, VALUE self)
{
    VALUE rb_new_dtype, rb_opt;
    rb_scan_args(argc, argv, "1:", &rb_new_dtype, &rb_opt);

    DTYPE new_dtype = mxx_dtype_from_symbol(rb_new_dtype);
    MX *dest = MX_ALLOC(MX);
    mxx_copy_cast(MX_DATA_PTR(self), dest, new_dtype);

    return Data_Wrap_Struct(CLASS_OF(self), 0, mxx_free, dest);
}

static void mxs_to_a(VALUE *ary, size_t shape_index, size_t array_index, MX *mx)
{
    if (shape_index == mx->dim - 1) {
        for (int i = 0; i < mx->shape[shape_index]; i++) {
            void *p = (char *)mx->elptr + (i + array_index) * DTYPE_SIZES[mx->dtype];
            ary[i] = mxx_c_to_rb(p, mx->dtype);
        }
    } else {
        size_t skip = 1;
        for (int i = shape_index + 1; i < mx->dim; i++) {
            skip *= mx->shape[i];
        }

        for (int i = 0; i < mx->shape[shape_index]; i++) {
            VALUE *v = MX_ALLOC_N(VALUE, mx->shape[shape_index + 1]);
            mxs_to_a(v, shape_index + 1, skip * i + array_index, mx);
            ary[i] = rb_ary_new4(mx->shape[shape_index + 1], v);
        }
    }
}

static VALUE mx_to_a(VALUE self)
{
    MX *mx = MX_DATA_PTR(self);
    VALUE *ary = MX_ALLOC_N(VALUE, mx->shape[0]);
    mxs_to_a(ary, 0, 0, mx);
    return rb_ary_new4(mx->shape[0], ary);
}

static VALUE mx_flatten(VALUE self)
{
    MX *mx = MX_DATA_PTR(self);
    VALUE *el = MX_ALLOC_N(VALUE, mx->size);
    size_t dsize = DTYPE_SIZES[mx->dtype];
    for (int i = 0; i < mx->size; i++) {
        el[i] = mxx_c_to_rb((char *)mx->elptr + i * dsize, mx->dtype);
    }
    return rb_ary_new4(mx->size, el);
}

static VALUE mx_first(VALUE self)
{
    MX *mx = MX_DATA_PTR(self);
    if (mx->size > 0) {
        return mxx_c_to_rb(mx->elptr, mx->dtype);
    }
    return Qnil;
}

static VALUE mx_second(VALUE self)
{
    MX *mx = MX_DATA_PTR(self);
    if (mx->size > 1) {
        return mxx_c_to_rb(mx->elptr + DTYPE_SIZES[mx->dtype], mx->dtype);
    }
    return Qnil;
}

static VALUE mx_last(VALUE self)
{
    MX *mx = MX_DATA_PTR(self);
    if (mx->size > 0) {
        return mxx_c_to_rb(mx->elptr + (mx->size - 1) * DTYPE_SIZES[mx->dtype], mx->dtype);
    }
    return Qnil;
}

static VALUE mx_reshape(VALUE self, VALUE shape)
{
    MX *mx = MX_DATA_PTR(self);

    MX *reshaped_mx = MX_ALLOC(MX);
    mxx_initialize_shape(reshaped_mx, shape);

    if (mx->size != reshaped_mx->size) {
        rb_raise(rb_eArgError, "Cannot do reshape operation between shape size of [ %ld ] and [ %ld ]. Total size of new array must be unchanged.", mx->size, reshaped_mx->size);
    }

    reshaped_mx->dtype = mx->dtype;
    reshaped_mx->elptr = MX_ALLOC_N(char, reshaped_mx->size * DTYPE_SIZES[reshaped_mx->dtype]);

    mxx_copy_elptr(mx, reshaped_mx);
    return Data_Wrap_Struct(CLASS_OF(self), 0, mxx_free, reshaped_mx);
}

static VALUE mx_sum(int argc, VALUE *argv, VALUE self)
{
    MX *mx = MX_DATA_PTR(self);

    if (DTYPE_IS_INT(mx->dtype)) {
        long sum = (long)mxx_sum(mx);
        return INT2FIX(sum);
    } else {
        return rb_float_new(mxx_sum(mx));
    }
}

static VALUE mx_is_empty(VALUE self)
{
    MX *mx = MX_DATA_PTR(self);
    return mx->size == 0 ? Qtrue : Qfalse;
}

static VALUE mx_set(int argc, VALUE *argv, VALUE self)
{
    return self;
}

/**
 * Slice the `src` matrix, then the sliced matrix will be stored into the dest matrix.
 * @params src   Source matrix to be sliced
 * @params dest  Destination matrix
 * @params src_pos   Copy start position in the source matrix
 * @params count Copy element count
 * @params depth ?
 * @params idx
 * @params idx_len
 */
static void _mxs_slice(MX *src, MX *dest, VALUE idx, size_t idx_len, int depth, size_t size,
        size_t src_pos, size_t *dest_pos)
{
    if (depth == idx_len) {
        mxx_copy_by_size(
            MX_DATA_POS(src,  src_pos),
            MX_DATA_POS(dest, *dest_pos),
            size,
            dest->dtype
        );

        *dest_pos += size;
    } else {
        size_t skip = 1;
        for (int i = depth + 1; i < src->dim; i++) {
            skip *= src->shape[i];
        }

        long beg, len;
        VALUE v = RARRAY_AREF(idx, depth);
        if (FIXNUM_P(v)) {
            _mxs_slice(src, dest, idx, idx_len, depth + 1, skip, src_pos + FIX2INT(v) * skip, dest_pos);
        } else if (rb_range_beg_len(v, &beg, &len, src->shape[depth], 0) == Qtrue) {
            for (long i = 0; i < len; i++) {
                _mxs_slice(src, dest, idx, idx_len, depth + 1, skip, src_pos + (beg + i) * skip, dest_pos);
            }
        }
    }
}

static MX *mxs_slice(MX *src, VALUE idx, size_t idx_len)
{
    /**
     * mx = [ [ [1,1], [2,2], [3,3], [4,4] ],
     *        [ [0,0], [0,0], [0,0], [0,0] ],
     *        [ [5,5], [6,6], [7,7], [8,8] ] ], shape = [3, 4, 2]
     *
     * mx[0]    =   [ [1,1], [2,2], [3,3], [4,4] ],   shape = [4, 2]
     * mx[0..1] = [ [ [1,1], [2,2], [3,3], [4,4] ],
     *              [ [0,0], [0,0], [0,0], [0,0] ] ], shape = [2, 4, 2]
     *
     * mx[0, 0..1] = [ [1,1], [2,2] ],   shape = [2, 2]
     *
     * mx[0..1, 0]    = [ [1,1],
     *                    [0,0] ], shape = [2, 2]
     * mx[0..1, 0..1] = [ [ [1,1], [2,2] ],
     *                    [ [0,0], [0,0] ] ], shape = [2, 2, 2]
     *
     * mx[0..1, 0..1, 0] = [ [ 1, 2 ],
     *                       [ 0, 0 ] ], shape = [2, 2]
     *
     * mx = [ [1, 2, 3],
     *        [4, 5, 6],
     *        [7, 8, 9] ], shape = [3, 3]
     *
     * mx[0]         =   [1, 2, 3],   shape = [3]
     * mx[0..1]      = [ [1, 2, 3]
     *                   [4, 5, 6] ], shape = [2, 3]
     * mx[0..1, 0]   = [ 1
     *                   4 ],         shape = [2, 1]
     * mx[0..1, [0]] = [ [1]
       *                 [4] ],       shape = [2, 1]
     *
     * mx[0, 0]    =     1,  shape = nil
     * mx[0, 0..1] = [1, 2], shape = [2]
     *
     * mx[[0]]    = [ [1, 2, 3] ], shape = [1, 3]
     * mx[0..0]   = [ [1, 2, 3]
     *                [4, 5, 6] ], shape = [2, 3]
     * mx[[0, 1]] = [ [1, 2, 3]
     *                [4, 5, 6] ], shape = [2, 3]
     *
     * 引数がFIXNUMだったらdimが1個減る
     * 引数がrangeまたはarrayだったら次元は減らずにshapeがその長さになる
     */
    long beg, len;
    size_t new_dim = 0;
    size_t *tmp = MX_ALLOC_N(size_t, src->dim);
    for (int i = 0; i < src->dim; i++) {
        if (i < idx_len) {
            VALUE v = RARRAY_AREF(idx, i);
            if (rb_range_beg_len(v, &beg, &len, src->shape[i], 0) == Qtrue) {
                tmp[new_dim++] = len;
            }
        } else {
            tmp[new_dim++] = src->shape[i];
        }
    }

    VALUE *new_shape = MX_ALLOC_N(VALUE, new_dim);
    for (int i = 0; i < new_dim; i++) {
        new_shape[i] = INT2FIX(tmp[i]);
    }

    MX *new_mx = MX_ALLOC(MX);
    mxx_initialize_shape(new_mx, rb_ary_new4(new_dim, new_shape));
    new_mx->dtype = src->dtype;
    new_mx->elptr = MX_ALLOC_N(char, new_mx->size * DTYPE_SIZES[new_mx->dtype]);
    MX_FREE(tmp);

    size_t dest_pos = 0;
    _mxs_slice(src, new_mx, idx, idx_len, 0, 0, 0, &dest_pos);
}

static VALUE mx_get(int argc, VALUE *argv, VALUE self)
{
    VALUE idx;
    rb_scan_args(argc, argv, "0*", &idx);
    MX *mx = MX_DATA_PTR(self);

    size_t idx_len = RARRAY_LEN(idx);
    if (idx_len > mx->dim) {
        return Qnil;
    }

    MX *new_mx = mxs_slice(mx, idx, idx_len);
    return Data_Wrap_Struct(CLASS_OF(self), 0, mxx_free, new_mx);

    return self;
}

static VALUE mx_ewadd(VALUE self, VALUE other)
{
    MX *mx = MX_DATA_PTR(self);
    MX *new_mx = MX_ALLOC(MX);

    if (IS_MX(other)) {
        MX *other_mx = MX_DATA_PTR(other);
        if (!mxx_is_same_shape(mx, other_mx)) {
            rb_raise(rb_eDataTypeError, "Cannot do addition operation between different shapes");
        }
        mxx_ewadd_array(mx, other_mx, new_mx);
    } else {
        if (TYPE(other) == T_FLOAT) {
            mxx_ewadd_dblscalar(mx, NUM2DBL(other), new_mx);
        } else {
            mxx_ewadd_intscalar(mx, FIX2INT(other), new_mx);
        }
    }

    return Data_Wrap_Struct(CLASS_OF(self), 0, mxx_free, new_mx);
}

static VALUE mx_ewsub(VALUE self, VALUE other)
{
    MX *mx = MX_DATA_PTR(self);
    MX *new_mx = MX_ALLOC(MX);

    if (IS_MX(other)) {
        MX *other_mx = MX_DATA_PTR(other);
        if (!mxx_is_same_shape(mx, other_mx)) {
            rb_raise(rb_eDataTypeError, "Cannot do subtraction operation between different shapes");
        }
        mxx_ewsub_array(mx, other_mx, new_mx);
    } else {
        if (TYPE(other) == T_FLOAT) {
            mxx_ewsub_dblscalar(mx, NUM2DBL(other), new_mx);
        } else {
            mxx_ewsub_intscalar(mx, FIX2INT(other), new_mx);
        }
    }

    return Data_Wrap_Struct(CLASS_OF(self), 0, mxx_free, new_mx);
}

static VALUE mx_ewmul(VALUE self, VALUE other)
{
    MX *mx = MX_DATA_PTR(self);
    MX *new_mx = MX_ALLOC(MX);

    if (IS_MX(other)) {
        MX *other_mx = MX_DATA_PTR(other);
        if (!mxx_is_same_shape(mx, other_mx)) {
            rb_raise(rb_eDataTypeError, "Cannot do multiplication operation between different shapes");
        }
        mxx_ewmul_array(mx, other_mx, new_mx);
    } else {
        if (TYPE(other) == T_FLOAT) {
            mxx_ewmul_dblscalar(mx, NUM2DBL(other), new_mx);
        } else {
            mxx_ewmul_intscalar(mx, FIX2INT(other), new_mx);
        }
    }

    return Data_Wrap_Struct(CLASS_OF(self), 0, mxx_free, new_mx);
}

static VALUE mx_ewpow(VALUE self, VALUE other)
{
    MX *mx = MX_DATA_PTR(self);
    MX *new_mx = MX_ALLOC(MX);

    if (IS_MX(other)) {
        MX *other_mx = MX_DATA_PTR(other);
        if (!mxx_is_same_shape(mx, other_mx)) {
            rb_raise(rb_eDataTypeError, "Cannot do power operation between different shapes");
        }

        if (DTYPE_IS_INT(other_mx->dtype)) {
            mxx_ewintpow_array(mx, other_mx, new_mx);
        } else {
            mxx_ewpow_array(mx, other_mx, new_mx);
        }
    } else {
        if (TYPE(other) == T_FIXNUM) {
            mxx_ewpow_intscalar(mx, FIX2INT(other), new_mx);
        } else {
            mxx_ewpow_dblscalar(mx, NUM2DBL(other), new_mx);
        }
    }

    return Data_Wrap_Struct(CLASS_OF(self), 0, mxx_free, new_mx);
}

static VALUE mx_sing_eye(int argc, VALUE *argv, VALUE klass)
{
    int row, col, dim = 2;
    VALUE rb_row, rb_col, rb_opt;
    rb_scan_args(argc, argv, "11:", &rb_row, &rb_col, &rb_opt);
    DTYPE dtype = mxx_dtype_from_opt(rb_opt, DTYPE_FLOAT64);

    VALUE *shape = MX_ALLOC_N(VALUE, dim);
    shape[0] = rb_row;
    if (NIL_P(rb_col)) {
        shape[1] = rb_row;
    } else {
        shape[1] = rb_col;
    }
    int count = MX_MIN(FIX2INT(shape[0]), FIX2INT(shape[1]));
    if (count < 0) {
        MX_FREE(shape);
        rb_raise(rb_eArgError, "Negative dimensions are not allowed");
    }

    MX *mx = mxx_initialize(rb_ary_new4(dim, shape), dtype);

    memset(mx->elptr, 0x00, mx->size * DTYPE_SIZES[mx->dtype]);
    mxx_eye(mx, count);
    MX_FREE(shape);
    return Data_Wrap_Struct(klass, 0, mxx_free, mx);;
}

static VALUE mx_sing_arange(int argc, VALUE *argv, VALUE klass)
{
    double start, stop, step;
    int start_type, stop_type, step_type;
    VALUE rb_start, rb_stop, rb_step, rb_opt;
    DTYPE dtype;

    rb_scan_args(argc, argv, "12:", &rb_start, &rb_stop, &rb_step, &rb_opt);

    if (NIL_P(rb_stop)) {
        start      = 0;
        stop       = NUM2DBL(rb_start);
        start_type = T_FIXNUM;
        stop_type  = TYPE(rb_start);
    } else {
        start      = NUM2DBL(rb_start);
        stop       = NUM2DBL(rb_stop);
        start_type = TYPE(rb_start);
        stop_type  = TYPE(rb_stop);
    }

    step      = NIL_P(rb_step) ? 1 : NUM2DBL(rb_step);
    step_type = NIL_P(rb_step) ? T_FIXNUM : TYPE(rb_step);

    if (step == 0) {
        rb_raise(rb_eArgError, "Cannot specify 0 to the step argument");
    } else if (step > 0 && start > stop) {
        rb_raise(rb_eArgError, "Confuse arguments order. start > stop.");
    } else if (step < 0 && start < stop) {
        rb_raise(rb_eArgError, "Confuse arguments order. start < stop.");
    } else if (start == stop) {
        MX *mx = mxx_initialize(INT2FIX(0), DTYPE_INT64);
        return Data_Wrap_Struct(klass, 0, mxx_free, mx);
    }

    // T_FIXNUM => FIXNUM_P
    if (step_type == T_FIXNUM && start_type == T_FIXNUM && stop_type == T_FIXNUM) {
        dtype = DTYPE_INT64;
    } else {
        dtype = DTYPE_FLOAT64;
    }

    size_t shape = (size_t)fabs(ceil((stop - start) / step));
    MX *mx = mxx_initialize(INT2FIX(shape), mxx_dtype_from_opt(rb_opt, dtype));
    mxx_linspace(mx, start, step);

    return Data_Wrap_Struct(klass, 0, mxx_free, mx);
}

static VALUE mx_sing_linspace(int argc, VALUE *argv, VALUE klass)
{
    VALUE rb_start, rb_stop, rb_num, rb_opt;
    rb_scan_args(argc, argv, "21:", &rb_start, &rb_stop, &rb_num, &rb_opt);

    if (NIL_P(rb_num)) {
        rb_num = INT2FIX(100);
    }

    if (FIX2INT(rb_num) < 0) {
        rb_raise(rb_eArgError, "Cannot specify negative number to the shape argument");
    }

    DTYPE dtype = mxx_dtype_from_opt(rb_opt, DTYPE_FLOAT64);
    MX *mx = mxx_initialize(rb_num, dtype);
    double start = NUM2DBL(rb_start);
    double stop  = NUM2DBL(rb_stop);
    double step  = (stop - start) * (1.0 / (double)(mx->size - 1));
    mxx_linspace(mx, start, step);

    return Data_Wrap_Struct(klass, 0, mxx_free, mx);
}

void Init_mxruby()
{
    rb_cMx = rb_define_class("MX", rb_cObject);
    rb_eDataTypeError = rb_define_class("DataTypeError", rb_eStandardError);

    rb_define_alloc_func(rb_cMx, mx_alloc);

    rb_define_protected_method(rb_cMx, "__set__", mx_initialize, 3);

    rb_define_method(rb_cMx, "shape", mx_shape, 0);
    rb_define_method(rb_cMx, "dim", mx_dim, 0);
    rb_define_method(rb_cMx, "size", mx_size, 0);
    rb_define_method(rb_cMx, "dtype", mx_dtype, 0);
    rb_define_method(rb_cMx, "dot", mx_dot, 1);
    rb_define_method(rb_cMx, "astype", mx_astype, -1);
    rb_define_method(rb_cMx, "to_a", mx_to_a, 0);
    rb_define_method(rb_cMx, "flatten", mx_flatten, 0);
    rb_define_method(rb_cMx, "first", mx_first, 0);
    rb_define_method(rb_cMx, "second", mx_second, 0);
    rb_define_method(rb_cMx, "last", mx_last, 0);
    rb_define_method(rb_cMx, "reshape", mx_reshape, 1);
    rb_define_method(rb_cMx, "sum", mx_sum, -1);

    rb_define_method(rb_cMx, "empty?", mx_is_empty, 0);

    rb_define_method(rb_cMx, "[]=", mx_set, -1);
    rb_define_method(rb_cMx, "[]", mx_get, -1);
    rb_define_method(rb_cMx, "+", mx_ewadd, 1);
    rb_define_method(rb_cMx, "-", mx_ewsub, 1);
    rb_define_method(rb_cMx, "*", mx_ewmul, 1);
    rb_define_method(rb_cMx, "**", mx_ewpow, 1);

    rb_define_singleton_method(rb_cMx, "eye", mx_sing_eye, -1);
    rb_define_singleton_method(rb_cMx, "arange", mx_sing_arange, -1);
    rb_define_singleton_method(rb_cMx, "linspace", mx_sing_linspace, -1);

    Init_random();
}

max_stars_repo_name: tkengo/mxruby

max_issues_repo_name: tkengo/mxruby

max_forks_repo_name: tkengo/mxruby


text:
#ifndef __GSL_SORT_VECTOR_H__
#define __GSL_SORT_VECTOR_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_sort_vector_long_double.h>
#include <gsl/gsl_sort_vector_double.h>
#include <gsl/gsl_sort_vector_float.h>

#include <gsl/gsl_sort_vector_ulong.h>
#include <gsl/gsl_sort_vector_long.h>

#include <gsl/gsl_sort_vector_uint.h>
#include <gsl/gsl_sort_vector_int.h>

#include <gsl/gsl_sort_vector_ushort.h>
#include <gsl/gsl_sort_vector_short.h>

#include <gsl/gsl_sort_vector_uchar.h>
#include <gsl/gsl_sort_vector_char.h>

#endif /* __GSL_SORT_VECTOR_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
#include <sstream>
#include <iostream>
#include <string>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_rng.h>
#include <iomanip>
#include <complex>
#include <cmath>
#include <random>

using namespace std::complex_literals;

// pretty print a gsl vector.
std::string ppVector(const gsl_vector *v) {
    std::ostringstream str;
    int r;

    str << "(";
    for (r = 0; r < (int) v->size; r++) {
        str << gsl_vector_get(v, r);
        if (r < v->size - 1) {
            str << ",";
        }
    }
    str << ")";
    return str.str();
}

void
vectorOpsCheck() {
    fprintf(stdout, "\n Vector Operations ..............................................\n");
    int N = 3;

    gsl_vector *v1 = gsl_vector_alloc(N);
    gsl_vector *v2 = gsl_vector_alloc(N);
    gsl_vector *v3 = gsl_vector_alloc(N);
    gsl_vector *v4 = gsl_vector_alloc(N);
    gsl_vector *v5 = gsl_vector_alloc(N);

    gsl_vector_set_all(v1, 1.0);
    gsl_vector_set_zero(v2);

    gsl_vector_set(v2, 1, 1.0);

    gsl_vector_set(v3, 0, 0.1);
    gsl_vector_set(v3, 1, 0.6);
    gsl_vector_set(v3, 2, 0.3);

    gsl_vector_set(v4, 0, 1.3);
    gsl_vector_set(v4, 1, 0.0);
    gsl_vector_set(v4, 2, 0.9);

    gsl_vector_set(v5, 0, 0.8);
    gsl_vector_set(v5, 1, 1.2);
    gsl_vector_set(v5, 2, 1.6);

    fprintf(stdout, "v1 = %s\n", ppVector(v1).c_str());
    fprintf(stdout, "v2 = %s\n", ppVector(v2).c_str());
    gsl_vector_add(v1, v2);
    fprintf(stdout, "v1 + v2 = %s\n", ppVector(v1).c_str());
    fprintf(stdout, "\n");

    fprintf(stdout, "v1 = %s\n", ppVector(v3).c_str());
    fprintf(stdout, "v2 = %s\n", ppVector(v5).c_str());
    gsl_vector_sub(v3, v5);
    fprintf(stdout, "v1 - v2 = %s\n", ppVector(v3).c_str());
    fprintf(stdout, "\n");

    fprintf(stdout, "v1 = %s\n", ppVector(v2).c_str());
    fprintf(stdout, "v2 = %s\n", ppVector(v4).c_str());
    gsl_vector_mul(v2, v4);
    fprintf(stdout, "v1 dot v2 = %s\n", ppVector(v2).c_str());
    fprintf(stdout, "\n");

    fprintf(stdout, "v1 = %s\n",  ppVector(v5).c_str());
    fprintf(stdout, "norm(v1) = %f\n", gsl_blas_dnrm2(v5));
    fprintf(stdout, "\n");

    gsl_vector_free(v1);
    gsl_vector_free(v2);
    gsl_vector_free(v3);
    gsl_vector_free(v4);
    gsl_vector_free(v5);

    fprintf(stdout, "\n");
}

void
complexNumOpsCheck() {
    fprintf(stdout, "\n Complex Number Operations ..............................................\n");

    std::complex<double> c1 = 1.1 + 2.0i;
    std::complex<double> c2 = 2.3 - 0.4i;
    fprintf(stdout, "let c1 = Complex(real: %1.1f, imag: %1.1f)\n", real(c1), imag(c1));
    fprintf(stdout, "let c2 = Complex(real: %1.1f, imag: %1.1f)\n", real(c2), imag(c2));
    fprintf(stdout, "XCTAssertEqual(c1+c2, c2+c1)\n");
    fprintf(stdout, "XCTAssertEqual(c2*c1, c2*c1)\n");

    fprintf(stdout, "XCTAssertEqual(c1+c2, Complex(real:%1.1f, imag: %1.1f))\n", real(c1+c2), imag(c1+c2));
    fprintf(stdout, "XCTAssertEqual(c1-c2, Complex(real:%1.1f, imag: %1.1f))\n", real(c1-c2), imag(c1-c2));
    fprintf(stdout, "XCTAssertEqual(c1*c2, Complex(real:%1.1f, imag: %1.1f))\n", real(c1*c2), imag(c1*c2));
    fprintf(stdout, "XCTAssertEqual(c1/c2, Complex(real:%1.6f, imag: %1.6f))\n", real(c1/c2), imag(c1/c2));
    fprintf(stdout, "XCTAssertEqual(c2/c1, Complex(real:%1.6f, imag: %1.5f))\n", real(c2/c1), imag(c2/c1));

    fprintf(stdout, "XCTAssertEqual(c1.conjugate(), Complex(real:%1.1f, imag: %1.1f))\n", real(conj(c1)), imag(conj(c1)));
    fprintf(stdout, "XCTAssertEqual(c1.modulus(), %1.14f\n", norm(c1));
    fprintf(stdout, "XCTAssertEqual(c1.norm(), %1.14f\n", norm(c1));
    fprintf(stdout, "XCTAssertEqualWithAccuracy(c1.arg(), %1.15f, accuracy:0.0001)\n", arg(c1));

    fprintf(stdout, "XCTAssertEqual(c2.conjugate(), Complex(real:%1.1f, imag: %1.1f))\n", real(conj(c2)), imag(conj(c2)));
    fprintf(stdout, "XCTAssertEqual(c2.modulus(), %1.14f\n", norm(c2));
    fprintf(stdout, "XCTAssertEqual(c2.norm(), %1.14f\n", norm(c2));
    fprintf(stdout, "XCTAssertEqualWithAccuracy(c2.arg(), %1.15f, accuracy:0.0001)\n", arg(c2));

    fprintf(stdout, "\n");
}

void
printComplexMathChecks(std::complex<double> c1) {

     fprintf(stdout, "var c1 = Complex(real: %1.1f, imag: %1.1f)\n", real(c1), imag(c1));
    fprintf(stdout, "let cosc1 = cos(value: c1)\n");
    fprintf(stdout, "let sinc1 = sin(value: c1)\n");
    fprintf(stdout, "let tanc1 = tan(value: c1)\n");
    fprintf(stdout, "XCTAssertEqual(cosc1.real(), %1.15f)\n", real(cos(c1)));
    fprintf(stdout, "XCTAssertEqual(cosc1.imag(), %1.15f)\n", imag(cos(c1)));
    fprintf(stdout, "XCTAssertEqual(sinc1.real(), %1.15f)\n", real(sin(c1)));
    fprintf(stdout, "XCTAssertEqual(sinc1.imag(), %1.15f)\n", imag(sin(c1)));
    fprintf(stdout, "XCTAssertEqual(tanc1.real(), %1.15f)\n", real(tan(c1)));
    fprintf(stdout, "XCTAssertEqual(tanc1.imag(), %1.15f)\n", imag(tan(c1)));

    fprintf(stdout, "let coshc1 = cosh(value: c1)\n");
    fprintf(stdout, "let sinhc1 = sinh(value: c1)\n");
    fprintf(stdout, "let tanhc1 = tanh(value: c1)\n");
    fprintf(stdout, "XCTAssertEqual(coshc1.real(), %1.15f)\n", real(cosh(c1)));
    fprintf(stdout, "XCTAssertEqual(coshc1.imag(), %1.15f)\n", imag(cosh(c1)));
    fprintf(stdout, "XCTAssertEqual(sinhc1.real(), %1.15f)\n", real(sinh(c1)));
    fprintf(stdout, "XCTAssertEqual(sinhc1.imag(), %1.15f)\n", imag(sinh(c1)));
    fprintf(stdout, "XCTAssertEqual(tanhc1.real(), %1.15f)\n", real(tanh(c1)));
    fprintf(stdout, "XCTAssertEqual(tanhc1.imag(), %1.15f)\n", imag(tanh(c1)));

    fprintf(stdout, "let acosc1 = acos(value: c1)\n");
    fprintf(stdout, "let asinc1 = asin(value: c1)\n");
    fprintf(stdout, "let atanc1 = atan(value: c1)\n");
    fprintf(stdout, "XCTAssertEqual(acosc1.real(), %1.15f)\n", real(acos(c1)));
    fprintf(stdout, "XCTAssertEqual(acosc1.imag(), %1.15f)\n", imag(acos(c1)));
    fprintf(stdout, "XCTAssertEqual(asinc1.real(), %1.15f)\n", real(asin(c1)));
    fprintf(stdout, "XCTAssertEqual(asinc1.imag(), %1.15f)\n", imag(asin(c1)));
    fprintf(stdout, "XCTAssertEqual(atanc1.real(), %1.15f)\n", real(atan(c1)));
    fprintf(stdout, "XCTAssertEqual(atanc1.imag(), %1.15f)\n", imag(atan(c1)));

    fprintf(stdout, "let acoshc1 = acosh(value: c1)\n");
    fprintf(stdout, "let asinhc1 = asinh(value: c1)\n");
    fprintf(stdout, "let atanhc1 = atanh(value: c1)\n");
    fprintf(stdout, "XCTAssertEqual(acoshc1.real(), %1.15f)\n", real(acosh(c1)));
    fprintf(stdout, "XCTAssertEqual(acoshc1.imag(), %1.15f)\n", imag(acosh(c1)));
    fprintf(stdout, "XCTAssertEqual(asinhc1.real(), %1.15f)\n", real(asinh(c1)));
    fprintf(stdout, "XCTAssertEqual(asinhc1.imag(), %1.15f)\n", imag(asinh(c1)));
    fprintf(stdout, "XCTAssertEqual(atanhc1.real(), %1.15f)\n", real(atanh(c1)));
    fprintf(stdout, "XCTAssertEqual(atanhc1.imag(), %1.15f)\n", imag(atanh(c1)));

    fprintf(stdout, "let expc1 = exp(value: c1)\n");
    fprintf(stdout, "let logc1 = log(value: c1)\n");
    fprintf(stdout, "let log10c1 = log10(value: c1)\n");
    fprintf(stdout, "let sqrtc1 = sqrt(value: c1)\n");
    fprintf(stdout, "XCTAssertEqual(expc1.real(), %1.15f)\n", real(exp(c1)));
    fprintf(stdout, "XCTAssertEqual(expc1.imag(), %1.15f)\n", imag(exp(c1)));
    fprintf(stdout, "XCTAssertEqual(logc1.real(), %1.15f)\n", real(log(c1)));
    fprintf(stdout, "XCTAssertEqual(logc1.imag(), %1.15f)\n", imag(log(c1)));
    fprintf(stdout, "XCTAssertEqual(log10c1.real(), %1.15f)\n", real(log10(c1)));
    fprintf(stdout, "XCTAssertEqual(log10c1.imag(), %1.15f)\n", imag(log10(c1)));
    fprintf(stdout, "XCTAssertEqual(sqrtc1.real(), %1.15f)\n", real(sqrt(c1)));
    fprintf(stdout, "XCTAssertEqual(sqrtc1.imag(), %1.15f)\n", imag(sqrt(c1)));
}

void
complexMathCheck() {

    fprintf(stdout, "\n Complex Number Trig Operations ..............................................\n");

    std::complex<double> c1 = 1.1 + 2.0i;
    printComplexMathChecks(c1);

    c1 = 2.3 - 0.4i;
    printComplexMathChecks(c1);

    fprintf(stdout, "\n");
}

void
rngCheck() {

    fprintf(stdout, "\n Random Number Generator Tests ..............................................\n");
    fprintf(stdout, "\n     Mersenne Twister 64 bit ..............................................\n");

    int n = 10;

    std::mt19937_64 generator (12345678);
    std::uniform_real_distribution<double> dis(0.0, 1.0);
    fprintf(stdout, "[%f",  dis(generator));
    for (int i = 1; i < n; i++) {
        fprintf(stdout, ", %f",  dis(generator));
    }
    fprintf(stdout, "]");

    fprintf(stdout, "\n");
}

int
main(void) {
    complexNumOpsCheck();
    vectorOpsCheck();
    rngCheck();
    complexMathCheck();

    return 0;
}

max_stars_repo_name: anthonyohare/SwiftScientificLibrary

max_issues_repo_name: anthonyohare/SwiftScientificLibrary

max_forks_repo_name: anthonyohare/SwiftScientificLibrary


text:
/*
/
/ adkGSL.h
/
/ homebrewed addons to gsl */

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include <math.h>
#include <assert.h>

int gsl_vector_includes(gsl_vector *aVec, double aValue);
double gsl_vector_sum(gsl_vector *aVec, int aVecSize);
double gsl_matrix_row_sum(gsl_matrix *mat, int row, int rowSize);
double log_sum(gsl_vector *vec);
double  *gsl_matrix_2_lapack(gsl_matrix *m);
double gsl_vector_dot_product(gsl_vector *vec1, gsl_vector *vec2);
void gsl_vector_outer_product(gsl_vector *vec1, gsl_vector *vec2, gsl_matrix *result);

gsl_matrix *gsl_matrix_power(gsl_matrix *A, int k);
gsl_matrix *gsl_matrix_power_logs(gsl_matrix *A, int k);
void gsl_matrix_lower_tri_copy(gsl_matrix *src, gsl_matrix *dest);
void gsl_matrix_upper_tri_copy(gsl_matrix *src, gsl_matrix *dest);
void gsl_matrix_lower_tri(gsl_matrix *x);
void gsl_matrix_upper_tri(gsl_matrix *x);
double gsl_matrix_trace(gsl_matrix *x);
void gsl_matrix_covariance(gsl_matrix *data, gsl_matrix *cov);

void fillMatrixFromArray(double *numbers, gsl_matrix *dest, int nrow, int ncol);
void fillMatrixFromVector(gsl_vector *numbers, gsl_matrix *dest, int nrow, int ncol);
void fillArrayFromMatrix(gsl_matrix *src, gsl_vector *dest);

void fillCholArrayFromMatrix(gsl_matrix *src, gsl_vector *dest);
void fillMatrixFromCholArray(double *numbers, gsl_matrix *dest, int nrow, int ncol);

//lapack wrappers
int dgeev(char jobvl, char jobvr, int n,  double *a, int lda, double *wr, double *wi, double *vl,
	  int ldvl,double *vr, int ldvr, double *work, int lwork);
int dhseqr(char job, char compz, int n, int ilo, int ihi,double *h, int ldh, double *wr, double *wi, double *z, int ldz, double *work, int lwork);

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


text:
// Copyright 2021 atframework
// Created by owent on

#pragma once

#include <stdint.h>
#include <cstddef>
#include <cstring>
#include <memory>
#include <string>
#include <type_traits>

#include <config/compiler/protobuf_prefix.h>

#include <rapidjson/document.h>

#include <config/compiler/protobuf_suffix.h>

#include <gsl/select-gsl.h>

#include <log/log_wrapper.h>

#include <libatbus.h>

#include "atframe/atapp_config.h"

namespace google {
namespace protobuf {
class Message;
class Timestamp;
class Duration;
}  // namespace protobuf
}  // namespace google

namespace atapp {

struct LIBATAPP_MACRO_API_HEAD_ONLY rapidsjon_loader_string_mode {
  enum type {
    RAW = 0,
    URI,
    URI_COMPONENT,
  };
};

struct LIBATAPP_MACRO_API_HEAD_ONLY rapidsjon_loader_load_options {
  bool reserve_empty;
  bool convert_large_number_to_string;  // it's friendly to JSON.parse(...) in javascript
  rapidsjon_loader_string_mode::type string_mode;

  inline rapidsjon_loader_load_options()
      : reserve_empty(false), convert_large_number_to_string(false), string_mode(rapidsjon_loader_string_mode::RAW) {}
};

struct LIBATAPP_MACRO_API_HEAD_ONLY rapidsjon_loader_dump_options {
  rapidsjon_loader_string_mode::type string_mode;
  bool convert_number_from_string;  // it's friendly to JSON.parse(...) in javascript
  inline rapidsjon_loader_dump_options()
      : string_mode(rapidsjon_loader_string_mode::RAW), convert_number_from_string(true) {}
};

LIBATAPP_MACRO_API std::string rapidsjon_loader_stringify(const rapidjson::Document &doc, size_t more_reserve_size = 0);
LIBATAPP_MACRO_API bool rapidsjon_loader_unstringify(rapidjson::Document &doc, const std::string &json);
LIBATAPP_MACRO_API const char *rapidsjon_loader_get_type_name(rapidjson::Type t);

LIBATAPP_MACRO_API std::string rapidsjon_loader_stringify(
    const ATBUS_MACRO_PROTOBUF_NAMESPACE_ID::Message &src,
    const rapidsjon_loader_load_options &options = rapidsjon_loader_load_options());
LIBATAPP_MACRO_API bool rapidsjon_loader_parse(
    ATBUS_MACRO_PROTOBUF_NAMESPACE_ID::Message &dst, const std::string &src,
    const rapidsjon_loader_dump_options &options = rapidsjon_loader_dump_options());

LIBATAPP_MACRO_API void rapidsjon_loader_mutable_set_member(rapidjson::Value &parent, gsl::string_view key,
                                                            rapidjson::Value &&val, rapidjson::Document &doc);
LIBATAPP_MACRO_API void rapidsjon_loader_mutable_set_member(rapidjson::Value &parent, gsl::string_view key,
                                                            const rapidjson::Value &val, rapidjson::Document &doc);
LIBATAPP_MACRO_API void rapidsjon_loader_mutable_set_member(rapidjson::Value &parent, gsl::string_view key,
                                                            gsl::string_view val, rapidjson::Document &doc);

LIBATAPP_MACRO_API void rapidsjon_loader_append_to_list(rapidjson::Value &list_parent, gsl::string_view val,
                                                        rapidjson::Document &doc);

LIBATAPP_MACRO_API void rapidsjon_loader_dump_to(
    const rapidjson::Document &src, ATBUS_MACRO_PROTOBUF_NAMESPACE_ID::Message &dst,
    const rapidsjon_loader_dump_options &options = rapidsjon_loader_dump_options());

LIBATAPP_MACRO_API void rapidsjon_loader_load_from(
    rapidjson::Document &dst, const ATBUS_MACRO_PROTOBUF_NAMESPACE_ID::Message &src,
    const rapidsjon_loader_load_options &options = rapidsjon_loader_load_options());

LIBATAPP_MACRO_API void rapidsjon_loader_dump_to(
    const rapidjson::Value &src, ATBUS_MACRO_PROTOBUF_NAMESPACE_ID::Message &dst,
    const rapidsjon_loader_dump_options &options = rapidsjon_loader_dump_options());

LIBATAPP_MACRO_API void rapidsjon_loader_load_from(
    rapidjson::Value &dst, rapidjson::Document &doc, const ATBUS_MACRO_PROTOBUF_NAMESPACE_ID::Message &src,
    const rapidsjon_loader_load_options &options = rapidsjon_loader_load_options());

// ============ template implement ============

template <class TVAL, bool>
struct rapidsjon_loader_mutable_member_helper;

template <class TVAL>
struct rapidsjon_loader_mutable_member_helper<TVAL, true> {
  static LIBATAPP_MACRO_API_HEAD_ONLY inline void set(rapidjson::Value &parent, gsl::string_view key, TVAL &&val,
                                                      rapidjson::Document &doc) {
    rapidsjon_loader_mutable_set_member(parent, key, gsl::string_view(std::forward<TVAL>(val)), doc);
  }

  static LIBATAPP_MACRO_API_HEAD_ONLY inline void append(rapidjson::Value &parent, TVAL &&val,
                                                         rapidjson::Document &doc) {
    rapidsjon_loader_append_to_list(parent, gsl::string_view(std::forward<TVAL>(val)), doc);
  }
};

template <class TVAL>
struct rapidsjon_loader_mutable_member_helper<TVAL, false> {
  static LIBATAPP_MACRO_API_HEAD_ONLY inline void set(rapidjson::Value &parent, gsl::string_view key, TVAL &&val,
                                                      rapidjson::Document &doc) {
    if (parent.IsNull()) {
      parent.SetObject();
    }

    if (!parent.IsObject()) {
      FWLOGERROR("parent should be a object, but we got {}.", rapidsjon_loader_get_type_name(parent.GetType()));
      return;
    }

    rapidjson::Value jkey;
    jkey.SetString(rapidjson::StringRef(key.data(), static_cast<rapidjson::SizeType>(key.size())));
    rapidjson::Value::MemberIterator iter = parent.FindMember(jkey);
    if (iter != parent.MemberEnd()) {
      iter->value.Set(std::forward<TVAL>(val), doc.GetAllocator());
    } else {
      rapidjson::Value k;
      k.SetString(key.data(), static_cast<rapidjson::SizeType>(key.size()), doc.GetAllocator());
      parent.AddMember(k, std::forward<TVAL>(val), doc.GetAllocator());
    }
  }

  static LIBATAPP_MACRO_API_HEAD_ONLY inline void append(rapidjson::Value &parent, TVAL &&val,
                                                         rapidjson::Document &doc) {
    if (parent.IsNull()) {
      parent.SetArray();
    }

    if (!parent.IsArray()) {
      FWLOGERROR("parent should be a array, but we got {}.", rapidsjon_loader_get_type_name(parent.GetType()));
      return;
    }

    parent.PushBack(std::forward<TVAL>(val), doc.GetAllocator());
  }
};

template <class TVAL, class = typename std::enable_if<
                          !std::is_same<typename std::decay<TVAL>::type, gsl::string_view>::value>::type>
LIBATAPP_MACRO_API_HEAD_ONLY void rapidsjon_loader_mutable_set_member(rapidjson::Value &parent, gsl::string_view key,
                                                                      TVAL &&val, rapidjson::Document &doc) {
  rapidsjon_loader_mutable_member_helper<TVAL, std::is_convertible<TVAL, gsl::string_view>::value>::set(
      parent, key, std::forward<TVAL>(val), doc);
}

template <class TVAL, class = typename std::enable_if<
                          !std::is_same<typename std::decay<TVAL>::type, gsl::string_view>::value>::type>
LIBATAPP_MACRO_API_HEAD_ONLY void rapidsjon_loader_append_to_list(rapidjson::Value &parent, TVAL &&val,
                                                                  rapidjson::Document &doc) {
  rapidsjon_loader_mutable_member_helper<TVAL, std::is_convertible<TVAL, gsl::string_view>::value>::append(
      parent, std::forward<TVAL>(val), doc);
}
}  // namespace atapp

max_stars_repo_name: atframework/libatapp

max_issues_repo_name: atframework/libatapp

max_forks_repo_name: atframework/libatapp


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include "ccl.h"
#include "ccl_params.h"
#include "ccl_redshifts.h"

// ---- LSST redshift distributions & current specs -----
// ---- Consider spline for input dN/dz - pending

/*------ ROUTINE: ccl_create_photoz_info ------
INPUT: void * params, (double *) pz_func (double, double, void *)
TASK: create a structure amalgamating the user-input information on the photo-z model.
The structure holds a pointer to the function which returns the probability of getting a certain z_ph (first double)
given a z_spec (second double), and a pointer to the parameters which get passed to that function (other than z_ph and z_sp); */
pz_info* ccl_create_photoz_info(void * params,
					   double (*pz_func)(double, double,void*, int*))
{
  pz_info * this_info = malloc(sizeof(pz_info));
  this_info ->your_pz_params = params;
  this_info -> your_pz_func = pz_func;

  return this_info;
}

/*------ ROUTINE: ccl_photoz -----
INPUT: double z_ph, void *params
TASK:  Returns the value of p(z_photo, z). Change this function to
       change the way true-z and photo-z's are related.
       This has to be in a form that gsl can integrate.
*/
// struct of parameters to pass to photo_z
struct pz_params{
  double z_true; // Gives the true redshift at which to evaluate
  pz_info * pz_information; //Calls the photo-z scatter model
  int *status;
};

static double ccl_photoz(double z_ph, void * params)
{
  struct pz_params * p = (struct pz_params *) params;
  pz_info * user_stuff = (pz_info*) p->pz_information;
  double z_s = p->z_true;

  return (user_stuff->your_pz_func)(z_ph, z_s, user_stuff->your_pz_params,p->status);
}

// Gaussian photo-z function
double gaussian_pz(double z_ph, double z_s, void* params, int *status){
    double sigma_z0 = *((double*) params);

    double sigma_z = sigma_z0 * (1. + z_s);
    return exp(- (z_ph - z_s)*(z_ph - z_s) / (2.*sigma_z*sigma_z)) \
         / (sqrt(2.*M_PI) *sigma_z);
}

/*------ ROUTINE: ccl_specs_create_gaussian_photoz_info ------
INPUT: void * user_pz_params, (double *) user_pz_func (double, double, void *)
TASK: Convenience function for creating a Gaussian photo-z model with error
sigma(z) = sigma_z0 (1 + z). */

pz_info* ccl_create_gaussian_photoz_info(double sigma_z0){

    // Allocate memory so that this value persists
    double* sigma_z0_copy = malloc(sizeof(double));
    *sigma_z0_copy = sigma_z0;

    // Construct pz_info struct
    pz_info * this_info = malloc(sizeof(pz_info));
    this_info->your_pz_params = sigma_z0_copy;
    this_info->your_pz_func = &gaussian_pz;
    return this_info;
}

/* ------ ROUTINE: ccl_free_photoz_info -------
INPUT: pz_info my_photoz_info
TASK: free memory holding the structure containing user-input photoz information */

void ccl_free_photoz_info(pz_info *my_photoz_info)
{
  free(my_photoz_info);
}

/*------ ROUTINE: ccl_create_dNdz_info ------
INPUT: void * params, (double *) dNdz_func (double, void *, int*)
TASK: create a structure amalgamating the information on an analytic true dNdz model.
The structure holds a pointer to the function which returns the analytic dNdz
* and a pointer to the parameters which get passed to that function (other than z); */
dNdz_info* ccl_create_dNdz_info(void * params, double(*dNdz_func)(double,void*,int*))
{
  dNdz_info * this_info = malloc(sizeof(dNdz_info));
  this_info ->your_dN_params = params;
  this_info -> your_dN_func = dNdz_func;

  return this_info;
}

/*------ ROUTINE: dNdz_smail ----------
 * INPUT: z, params (containing: alpha, beta, z0), status
 * OUTPUT: Analytic Smail et al. type dNdz (NOT normalized) */

double dNdz_smail(double z, void* params, int *status){
    double alpha = ((smail_params*) params)->alpha;
    double beta = ((smail_params*) params)->beta;
    double z0 = ((smail_params*) params)->z0;

    return pow(z, alpha) * exp(- pow(z/z0, beta) );
}

/*------ ROUTINE: ccl_create_smail_dNdz_info ------
INPUT: alpha, beta, z0
TASK: Convenience function for creating an analytic dNdz wrt true z
* of the 'smail' form: dNdz ~ z^alpha exp(- (z/z0)^beta) */

dNdz_info* ccl_create_Smail_dNdz_info(double alpha, double beta, double z0){

    // Allocate a smail type parmaeter structure
    smail_params * smail_par = malloc(sizeof(smail_params));
    smail_par->alpha = alpha;
    smail_par->beta = beta;
    smail_par->z0 = z0;

    // Construct dNdz_info struct
    dNdz_info * this_info = malloc(sizeof(dNdz_info));
    this_info->your_dN_params = smail_par;
    this_info->your_dN_func = &dNdz_smail;
    return this_info;
}

/*------ ROUTINE: ccl_dNdz -----
INPUT: double z, void *params
TASK:  Returns the value of dNdz(z). Change this function to
       change the way true-z and photo-z's are related.
       This has to be in a form that gsl can integrate.
*/

static double ccl_dNdz(double z, dNdz_info*  params, int* status)
{

  return (params->your_dN_func)(z, params->your_dN_params, status);
}

/* ------ ROUTINE: ccl_free_dNdz_info -------
INPUT: dNdz_info my_dNdz_info
TASK: free memory holding the structure containing dNdz information */

void ccl_free_dNdz_info(dNdz_info *my_dNdz_info)
{
  free(my_dNdz_info);
}

/*------ ROUTINE: ccl_specs_norm_integrand -----
INPUT: double z_ph, void *params
TASK:  Returns the integrand which is integrated to get the normalization of
       dNdz in a given photometric redshift bin (the denominator from dNdz_sources_tomog).
       This has to be an separate function that gsl can integrate.
*/

// struct of parameters to pass to norm_integrand
struct norm_params {
  double bin_zmin_;
  double bin_zmax_;
  pz_info * pz_information;
  dNdz_info * dN_information;
  int *status;
};

static double ccl_norm_integrand(double z, void* params)
{

  // This is a struct that contains a true redshift and a pointer to the information about the photo_z model
  struct pz_params pz_val_p; // parameters for the photoz pdf wrt true-z

  double pz_int=0; // pointer to the value of the integral over the photoz model
  struct norm_params *p = (struct norm_params *) params; // parameters of the current function (because of form required for gsl integration)

  double z_min = p->bin_zmin_;
  double z_max = p->bin_zmax_;

  // Check whether ccl_splines and ccl_gsl exist; exit gracefully if they
  // can't be loaded
  if(ccl_splines==NULL || ccl_gsl==NULL) ccl_cosmology_read_config();
  if(ccl_splines==NULL || ccl_gsl==NULL) {
    ccl_raise_exception(CCL_ERROR_MISSING_CONFIG_FILE,
                        "ccl_redshift.c: Failed to read config file.");
    return NAN;
  }

  // Set up parameters for the pz part of the intermediary integral.
  pz_val_p.z_true = z;
  pz_val_p.status = p->status;
  pz_val_p.pz_information = p-> pz_information;

  // Do the intermediary integral over the model relating  photo-z to true-z
  gsl_integration_cquad_workspace * workspace = gsl_integration_cquad_workspace_alloc(ccl_gsl->N_ITERATION);
  gsl_function F;
  F.function = ccl_photoz;
  F.params = &pz_val_p;
  int gslstatus = gsl_integration_cquad(&F, z_min, z_max, 0.0,ccl_gsl->INTEGRATION_DNDZ_EPSREL,workspace,&pz_int, NULL, NULL);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_redshifts.c: ccl_specs_norm_integrand():");
    *p->status|= gslstatus;
  }
  gsl_integration_cquad_workspace_free(workspace);

  return ccl_dNdz(z, p->dN_information, p->status) * pz_int ;
}

/*------ ROUTINE: ccl_specs_dNdz_tomog -----
INPUT: double z, , double bin_zmin, double bin_zmax, dNdz function pointer, sigma_z function pointer
       tomographic boundaries are [bin_zmin,bin_zmax]
TASK:  dNdz in a particular tomographic bin,
       convolved with a photo-z model (defined by the user), and normalized.
       returns a status integer 0 if called with an allowable type of dNdz, non-zero otherwise
       (this is different from the regular status handling procedure because we don't pass a cosmology to this function)
*/
void ccl_dNdz_tomog(double z, double bin_zmin, double bin_zmax,
              pz_info * photo_info,  dNdz_info * dN_info, double *tomoout, int *status)
{
  // This uses equation 33 of Joachimi & Schneider 2009, arxiv:0905.0393
  double numerator_integrand=0, denom_integrand=0, dNdz_t;
  // This struct contains a spec redshift and a pointer to a photoz information struct.
  struct pz_params pz_p_val; //parameters for the integral over the photoz's
  struct norm_params norm_p_val;

  // Check whether ccl_splines and ccl_gsl exist; exit gracefully if they
  // can't be loaded
  if(ccl_splines==NULL || ccl_gsl==NULL) ccl_cosmology_read_config();
  if(ccl_splines==NULL || ccl_gsl==NULL) {
    ccl_raise_exception(CCL_ERROR_MISSING_CONFIG_FILE,
                        "ccl_redshifts.c: Failed to read config file.");
    *status = CCL_ERROR_MISSING_CONFIG_FILE;
    return;
  }

  // Set up the parameters to pass to the normalising integral (of type struct norm_params
  norm_p_val.bin_zmin_=bin_zmin;
  norm_p_val.bin_zmax_=bin_zmax;
  norm_p_val.pz_information = photo_info;
  norm_p_val.status = status;
  norm_p_val.dN_information = dN_info;

  dNdz_t = ccl_dNdz(z, dN_info, status);

  // Set up the parameters for the integral over the photo z function in the numerator (of type struct pz_params)
  pz_p_val.z_true = z;
  pz_p_val.status = status;
  pz_p_val.pz_information = photo_info; // pointer to user information

  // Integrate over the assumed pdf of photo-z wrt true-z in this bin (this goes in the numerator of the result):
  gsl_integration_cquad_workspace * workspace = gsl_integration_cquad_workspace_alloc(ccl_gsl->N_ITERATION);
  gsl_function F;
  F.function = ccl_photoz;
  F.params = &pz_p_val;
  int gslstatus = gsl_integration_cquad(&F, bin_zmin, bin_zmax, 0.0,ccl_gsl->INTEGRATION_DNDZ_EPSREL,workspace,&numerator_integrand, NULL, NULL);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_redshifts.c: ccl_norm_integrand():");
    *status |= gslstatus;
  }
  gsl_integration_cquad_workspace_free(workspace);

  // Now get the denominator, which normalizes dNdz over the photometric bin
  workspace = gsl_integration_cquad_workspace_alloc(ccl_gsl->N_ITERATION);
  F.function = ccl_norm_integrand;
  F.params = &norm_p_val;
  gslstatus = gsl_integration_cquad(&F, Z_MIN_SOURCES, Z_MAX_SOURCES, 0.0,ccl_gsl->INTEGRATION_DNDZ_EPSREL,workspace,&denom_integrand, NULL, NULL);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_redshifts.c: ccl_norm_integrand():");
    *status |= gslstatus;
  }

  gsl_integration_cquad_workspace_free(workspace);
  if (*status) {
    *status = CCL_ERROR_INTEG;
    return;
  }
  *tomoout = dNdz_t * numerator_integrand / denom_integrand;

}

max_stars_repo_name: vrastil/CCL

max_issues_repo_name: vrastil/CCL

max_forks_repo_name: vrastil/CCL


text:
// NIT_inspiral - code to rapidly compute extreme mass-ratio inspirals using self-force results
// Copyright (C) 2017  Niels Warburton
//
// 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, see <https://www.gnu.org/licenses/>.

#ifndef __INTERPOLANT_H__
#define __INTERPOLANT_H__

#include <vector>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

using namespace std;

typedef vector<double> Vector;

class Interpolant{
	public:
		// 1D interpolation
		Interpolant(Vector x, Vector f);
		double eval(double x);

		// 2D interpolation
		Interpolant(Vector x, Vector y, Vector f);
		double eval(double x, double y);

		// Destructor
		~Interpolant();

	private:
		int interp_type;	// Set to 1 for 1D interpolation and 2 for 2D interpolation

		gsl_spline *spline;
		gsl_spline2d *spline2d;
		gsl_interp_accel *xacc;
		gsl_interp_accel *yacc;
};
#endif // __INTERPOLANT_H__


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:44:57 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -twiddleinv 2 */

/*
 * This function contains 6 FP additions, 4 FP multiplications,
 * (or, 4 additions, 2 multiplications, 2 fused multiply/add),
 * 10 stack variables, and 8 memory accesses
 */

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftwi_twiddle_2(fftw_complex *A, const fftw_complex *W, int iostride, int m, int dist)
{
     int i;
     fftw_complex *inout;
     inout = A;
     for (i = m; i > 0; i = i - 1, inout = inout + dist, W = W + 1) {
	  fftw_real tmp1;
	  fftw_real tmp8;
	  fftw_real tmp6;
	  fftw_real tmp7;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp1 = c_re(inout[0]);
	  tmp8 = c_im(inout[0]);
	  {
	       fftw_real tmp3;
	       fftw_real tmp5;
	       fftw_real tmp2;
	       fftw_real tmp4;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp3 = c_re(inout[iostride]);
	       tmp5 = c_im(inout[iostride]);
	       tmp2 = c_re(W[0]);
	       tmp4 = c_im(W[0]);
	       tmp6 = (tmp2 * tmp3) + (tmp4 * tmp5);
	       tmp7 = (tmp2 * tmp5) - (tmp4 * tmp3);
	  }
	  c_re(inout[iostride]) = tmp1 - tmp6;
	  c_re(inout[0]) = tmp1 + tmp6;
	  c_im(inout[0]) = tmp7 + tmp8;
	  c_im(inout[iostride]) = tmp8 - tmp7;
     }
}

static const int twiddle_order[] =
{1};
fftw_codelet_desc fftwi_twiddle_2_desc =
{
     "fftwi_twiddle_2",
     (void (*)()) fftwi_twiddle_2,
     2,
     FFTW_BACKWARD,
     FFTW_TWIDDLE,
     55,
     1,
     twiddle_order,
};

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt