text:
// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
// reserved. See files LICENSE and NOTICE for details.
//
// This file is part of CEED, a collection of benchmarks, miniapps, software
// libraries and APIs for efficient high-order finite element and spectral
// element discretizations for exascale applications. For more information and
// source code availability see http://github.com/ceed.
//
// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
// a collaborative effort of two U.S. Department of Energy organizations (Office
// of Science and the National Nuclear Security Administration) responsible for
// the planning and preparation of a capable exascale ecosystem, including
// software, applications, hardware, advanced system engineering and early
// testbed platforms, in support of the nation's exascale computing imperative.

#ifndef setup_h
#define setup_h

#include <stdbool.h>
#include <string.h>

#include <petsc.h>
#include <petscdmplex.h>
#include <petscksp.h>
#include <petscfe.h>

#include <ceed.h>

#if PETSC_VERSION_LT(3,14,0)
#  define DMAddBoundary(a,b,c,d,e,f,g,h,i,j,k,l) DMAddBoundary(a,b,c,d,e,f,g,h,j,k,l)
#endif

#ifndef PHYSICS_STRUCT
#define PHYSICS_STRUCT
typedef struct Physics_private *Physics;
struct Physics_private {
  CeedScalar   nu;      // Poisson's ratio
  CeedScalar   E;       // Young's Modulus
};
#endif

// -----------------------------------------------------------------------------
// Command Line Options
// -----------------------------------------------------------------------------
// Problem options
typedef enum {
  ELAS_LIN = 0, ELAS_HYPER_SS = 1, ELAS_HYPER_FS = 2
} problemType;
static const char *const problemTypes[] = {"linElas",
                                           "hyperSS",
                                           "hyperFS",
                                           "problemType","ELAS_",0
                                          };
static const char *const problemTypesForDisp[] = {"Linear elasticity",
                                                  "Hyper elasticity small strain",
                                                  "Hyper elasticity finite strain"
                                                 };

// Forcing function options
typedef enum {
  FORCE_NONE = 0, FORCE_CONST = 1, FORCE_MMS = 2
} forcingType;
static const char *const forcingTypes[] = {"none",
                                           "constant",
                                           "mms",
                                           "forcingType","FORCE_",0
                                          };
static const char *const forcingTypesForDisp[] = {"None",
                                                  "Constant",
                                                  "Manufactured solution"
                                                 };

// Multigrid options
typedef enum {
  MULTIGRID_LOGARITHMIC = 0, MULTIGRID_UNIFORM = 1, MULTIGRID_NONE = 2
} multigridType;
static const char *const multigridTypes [] = {"logarithmic",
                                              "uniform",
                                              "none",
                                              "multigridType","MULTIGRID",0
                                             };
static const char *const multigridTypesForDisp[] = {"P-multigrid, logarithmic coarsening",
                                                    "P-multigrind, uniform coarsening",
                                                    "No multigrid"
                                                   };

typedef PetscErrorCode BCFunc(PetscInt, PetscReal, const PetscReal *, PetscInt,
                              PetscScalar *, void *);
// Note: These variables should be updated if additional boundary conditions
//         are added to boundary.c.
BCFunc BCMMS, BCZero, BCClamp;

// MemType Options
static const char *const memTypes[] = {"host","device","memType",
                                       "CEED_MEM_",0
                                      };

// -----------------------------------------------------------------------------
// Structs
// -----------------------------------------------------------------------------
// Units
typedef struct Units_private *Units;
struct Units_private {
  // Fundamental units
  PetscScalar meter;
  PetscScalar kilogram;
  PetscScalar second;
  // Derived unit
  PetscScalar Pascal;
};

// Application context from user command line options
typedef struct AppCtx_private *AppCtx;
struct AppCtx_private {
  char          ceedResource[PETSC_MAX_PATH_LEN];     // libCEED backend
  char          meshFile[PETSC_MAX_PATH_LEN];         // exodusII mesh file
  PetscBool     testMode;
  PetscBool     viewSoln;
  PetscBool     viewFinalSoln;
  problemType   problemChoice;
  forcingType   forcingChoice;
  multigridType multigridChoice;
  PetscScalar   nuSmoother;
  PetscInt      degree;
  PetscInt      qextra;
  PetscInt      numLevels;
  PetscInt      *levelDegrees;
  PetscInt      numIncrements;                        // Number of steps
  PetscInt      bcClampCount;
  PetscInt      bcClampFaces[16];
  PetscScalar   bcClampMax[16][7];
  PetscInt      bcTractionCount;
  PetscInt      bcTractionFaces[16];
  PetscScalar   bcTractionVector[16][3];
  PetscScalar   forcingVector[3];
  PetscBool     petscHaveCuda, setMemTypeRequest;
  CeedMemType   memTypeRequested;
};

// Problem specific data
// *INDENT-OFF*
typedef struct {
  CeedInt           qdatasize;
  CeedQFunctionUser setupgeo, apply, jacob, energy, diagnostic;
  const char        *setupgeofname, *applyfname, *jacobfname, *energyfname,
                    *diagnosticfname;
  CeedQuadMode      qmode;
} problemData;
// *INDENT-ON*

// Data specific to each problem option
extern problemData problemOptions[3];

// Forcing function data
typedef struct {
  CeedQFunctionUser setupforcing;
  const char        *setupforcingfname;
} forcingData;

extern forcingData forcingOptions[3];

// Data for PETSc Matshell
typedef struct UserMult_private *UserMult;
struct UserMult_private {
  MPI_Comm        comm;
  DM              dm;
  Vec             Xloc, Yloc, NBCs;
  CeedVector      Xceed, Yceed;
  CeedOperator    op;
  CeedQFunction   qf;
  Ceed            ceed;
  PetscScalar     loadIncrement;
  CeedQFunctionContext ctxPhys, ctxPhysSmoother;
  CeedMemType     memType;
  int (*VecGetArray)(Vec, PetscScalar **);
  int (*VecGetArrayRead)(Vec, const PetscScalar **);
  int (*VecRestoreArray)(Vec, PetscScalar **);
  int (*VecRestoreArrayRead)(Vec, const PetscScalar **);
};

// Data for Jacobian setup routine
typedef struct FormJacobCtx_private *FormJacobCtx;
struct FormJacobCtx_private {
  UserMult     *jacobCtx;
  PetscInt     numLevels;
  SNES         snesCoarse;
  Mat          *jacobMat, jacobMatCoarse;
  Vec          Ucoarse;
};

// Data for PETSc Prolongation/Restriction Matshell
typedef struct UserMultProlongRestr_private *UserMultProlongRestr;
struct UserMultProlongRestr_private {
  MPI_Comm     comm;
  DM           dmC, dmF;
  Vec          locVecC, locVecF;
  CeedVector   ceedVecC, ceedVecF;
  CeedOperator opProlong, opRestrict;
  Ceed         ceed;
  CeedMemType   memType;
  int (*VecGetArray)(Vec, PetscScalar **);
  int (*VecGetArrayRead)(Vec, const PetscScalar **);
  int (*VecRestoreArray)(Vec, PetscScalar **);
  int (*VecRestoreArrayRead)(Vec, const PetscScalar **);
};

// libCEED data struct for level
typedef struct CeedData_private *CeedData;
struct CeedData_private {
  Ceed                ceed;
  CeedBasis           basisx, basisu, basisCtoF, basisEnergy, basisDiagnostic;
  CeedElemRestriction Erestrictx, Erestrictu, Erestrictqdi,
                      ErestrictGradui, ErestrictEnergy, ErestrictDiagnostic,
                      ErestrictqdDiagnostici;
  CeedQFunction       qfApply, qfJacob, qfEnergy, qfDiagnostic;
  CeedOperator        opApply, opJacob, opRestrict, opProlong, opEnergy,
                      opDiagnostic;
  CeedVector          qdata, qdataDiagnostic, gradu, xceed, yceed, truesoln;
};

// -----------------------------------------------------------------------------
// Process command line options
// -----------------------------------------------------------------------------
// Process general command line options
PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx appCtx);

// Process physics options
PetscErrorCode ProcessPhysics(MPI_Comm comm, Physics phys, Units units);

// -----------------------------------------------------------------------------
// Setup DM
// -----------------------------------------------------------------------------
PetscErrorCode CreateBCLabel(DM dm, const char name[]);

// Create FE by degree
PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc,
                                     PetscBool isSimplex, const char prefix[],
                                     PetscInt order, PetscFE *fem);

// Read mesh and distribute DM in parallel
PetscErrorCode CreateDistributedDM(MPI_Comm comm, AppCtx appCtx, DM *dm);

// Setup DM with FE space of appropriate degree
PetscErrorCode SetupDMByDegree(DM dm, AppCtx appCtx, PetscInt order,
                               PetscBool boundary, PetscInt ncompu);

// -----------------------------------------------------------------------------
// libCEED Functions
// -----------------------------------------------------------------------------
// Destroy libCEED objects
PetscErrorCode CeedDataDestroy(CeedInt level, CeedData data);

// Utility function - essential BC dofs are encoded in closure indices as -(i+1)
PetscInt Involute(PetscInt i);

// Utility function to create local CEED restriction from DMPlex
PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P,
    CeedInt height, DMLabel domainLabel, CeedInt value,
    CeedElemRestriction *Erestrict);

// Utility function to get Ceed Restriction for each domain
PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height,
                                       DMLabel domainLabel, PetscInt value, CeedInt P, CeedInt Q, CeedInt qdatasize,
                                       CeedElemRestriction *restrictq, CeedElemRestriction *restrictx,
                                       CeedElemRestriction *restrictqdi);

// Set up libCEED for a given degree
PetscErrorCode SetupLibceedFineLevel(DM dm, DM dmEnergy, DM dmDiagnostic,
                                     Ceed ceed, AppCtx appCtx,
                                     CeedQFunctionContext physCtx,
                                     CeedData *data, PetscInt fineLevel,
                                     PetscInt ncompu, PetscInt Ugsz,
                                     PetscInt Ulocsz, CeedVector forceCeed,
                                     CeedVector neumannCeed);

// Set up libCEED multigrid level for a given degree
PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx appCtx,
                                 CeedData *data, PetscInt level,
                                 PetscInt ncompu, PetscInt Ugsz,
                                 PetscInt Ulocsz, CeedVector fineMult);

// Setup context data for Jacobian evaluation
PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx appCtx, DM dm, Vec V,
                                Vec Vloc, CeedData ceedData, Ceed ceed,
                                CeedQFunctionContext ctxPhys,
                                CeedQFunctionContext ctxPhysSmoother,
                                UserMult jacobianCtx);

// Setup context data for prolongation and restriction operators
PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx appCtx, DM dmC,
                                       DM dmF, Vec VF, Vec VlocC, Vec VlocF,
                                       CeedData ceedDataC, CeedData ceedDataF,
                                       Ceed ceed,
                                       UserMultProlongRestr prolongRestrCtx);

// -----------------------------------------------------------------------------
// Jacobian setup
// -----------------------------------------------------------------------------
PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat Jpre, void *ctx);

// -----------------------------------------------------------------------------
// Solution output
// -----------------------------------------------------------------------------
PetscErrorCode ViewSolution(MPI_Comm comm, Vec U, PetscInt increment,
                            PetscScalar loadIncrement);

PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU,
                                        UserMult user, Vec U,
                                        CeedElemRestriction ErestrictDiagnostic);

// -----------------------------------------------------------------------------
// libCEED Operators for MatShell
// -----------------------------------------------------------------------------
// This function uses libCEED to compute the local action of an operator
PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user);

// This function uses libCEED to compute the non-linear residual
PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx);

// This function uses libCEED to apply the Jacobian for assembly via a SNES
PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx);

// This function uses libCEED to compute the action of the Jacobian
PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y);

// This function uses libCEED to compute the action of the prolongation operator
PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y);

// This function uses libCEED to compute the action of the restriction operator
PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y);

// This function returns the computed diagonal of the operator
PetscErrorCode GetDiag_Ceed(Mat A, Vec D);

// This function calculates the strain energy in the final solution
PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user,
                                   CeedOperator opEnergy, Vec X,
                                   PetscReal *energy);

// -----------------------------------------------------------------------------
// Boundary Functions
// -----------------------------------------------------------------------------
// Note: If additional boundary conditions are added, an update is needed in
//         elasticity.h for the boundaryOptions variable.

// BCMMS - boundary function
// Values on all points of the mesh is set based on given solution below
// for u[0], u[1], u[2]
PetscErrorCode BCMMS(PetscInt dim, PetscReal loadIncrement,
                     const PetscReal coords[], PetscInt ncompu,
                     PetscScalar *u, void *ctx);

// BCClamp - fix boundary values with affine transformation at fraction of load
//   increment
PetscErrorCode BCClamp(PetscInt dim, PetscReal loadIncrement,
                       const PetscReal coords[], PetscInt ncompu,
                       PetscScalar *u, void *ctx);

#endif //setup_h


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_histogram2d.h>

#include "find.c"

int
gsl_histogram2d_pdf_sample (const gsl_histogram2d_pdf * p,
                            double r1, double r2,
                            double *x, double *y)
{
  size_t k;
  int status;

/* Wrap the exclusive top of the bin down to the inclusive bottom of
   the bin. Since this is a single point it should not affect the
   distribution. */

  if (r2 == 1.0)
    {
      r2 = 0.0;
    }
  if (r1 == 1.0)
    {
      r1 = 0.0;
    }

  status = find (p->nx * p->ny, p->sum, r1, &k);

  if (status)
    {
      GSL_ERROR ("cannot find r1 in cumulative pdf", GSL_EDOM);
    }
  else
    {
      size_t i = k / p->ny;
      size_t j = k - (i * p->ny);
      double delta = (r1 - p->sum[k]) / (p->sum[k + 1] - p->sum[k]);
      *x = p->xrange[i] + delta * (p->xrange[i + 1] - p->xrange[i]);
      *y = p->yrange[j] + r2 * (p->yrange[j + 1] - p->yrange[j]);
      return GSL_SUCCESS;
    }
}

gsl_histogram2d_pdf *
gsl_histogram2d_pdf_alloc (const size_t nx, const size_t ny)
{
  const size_t n = nx * ny;

  gsl_histogram2d_pdf *p;

  if (n == 0)
    {
      GSL_ERROR_VAL ("histogram2d pdf length n must be positive integer",
                        GSL_EDOM, 0);
    }

  p = (gsl_histogram2d_pdf *) malloc (sizeof (gsl_histogram2d_pdf));

  if (p == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf struct",
                        GSL_ENOMEM, 0);
    }

  p->xrange = (double *) malloc ((nx + 1) * sizeof (double));

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

      GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf xranges",
                        GSL_ENOMEM, 0);
    }

  p->yrange = (double *) malloc ((ny + 1) * sizeof (double));

  if (p->yrange == 0)
    {
      free (p->xrange);
      free (p);         /* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf yranges",
                        GSL_ENOMEM, 0);
    }

  p->sum = (double *) malloc ((n + 1) * sizeof (double));

  if (p->sum == 0)
    {
      free (p->yrange);
      free (p->xrange);
      free (p);         /* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf sums",
                        GSL_ENOMEM, 0);
    }

  p->nx = nx;
  p->ny = ny;

  return p;
}

int
gsl_histogram2d_pdf_init (gsl_histogram2d_pdf * p, const gsl_histogram2d * h)
{
  size_t i;
  const size_t nx = p->nx;
  const size_t ny = p->ny;
  const size_t n = nx * ny;

  if (nx != h->nx || ny != h->ny)
    {
      GSL_ERROR ("histogram2d size must match pdf size", GSL_EDOM);
    }

  for (i = 0; i < n; i++)
    {
      if (h->bin[i] < 0)
        {
          GSL_ERROR ("histogram bins must be non-negative to compute"
                     "a probability distribution", GSL_EDOM);
        }
    }

  for (i = 0; i < nx + 1; i++)
    {
      p->xrange[i] = h->xrange[i];
    }

  for (i = 0; i < ny + 1; i++)
    {
      p->yrange[i] = h->yrange[i];
    }

  {
    double mean = 0, sum = 0;

    for (i = 0; i < n; i++)
      {
        mean += (h->bin[i] - mean) / ((double) (i + 1));
      }

    p->sum[0] = 0;

    for (i = 0; i < n; i++)
      {
        sum += (h->bin[i] / mean) / n;
        p->sum[i + 1] = sum;
      }
  }

  return GSL_SUCCESS;
}

void
gsl_histogram2d_pdf_free (gsl_histogram2d_pdf * p)
{
  RETURN_IF_NULL (p);
  free (p->xrange);
  free (p->yrange);
  free (p->sum);
  free (p);
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

int
main (void)
{
  double c  = GSL_CONST_MKSA_SPEED_OF_LIGHT;
  double au = GSL_CONST_MKSA_ASTRONOMICAL_UNIT;
  double minutes = GSL_CONST_MKSA_MINUTE;

  /* distance stored in meters */
  double r_earth = 1.00 * au;
  double r_mars  = 1.52 * au;

  double t_min, t_max;

  t_min = (r_mars - r_earth) / c;
  t_max = (r_mars + r_earth) / c;

  printf ("light travel time from Earth to Mars:\n");
  printf ("minimum = %.1f minutes\n", t_min / minutes);
  printf ("maximum = %.1f minutes\n", t_max / minutes);

  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified


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

#include "JoinComp.h"
#include "Lambda.h"
#include "LDADocWordTopicAssignment.h"
#include "LDA/LDATopicWordProb.h"
#include "LambdaCreationFunctions.h"
#include "LDADocument.h"
#include "IntDoubleVectorPair.h"
#include "PDBVector.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <random>
#include <gsl/gsl_vector.h>
#include <algorithm>
#include <iostream>
#include <math.h>

/* This class implements the join between documents, doc-topic probability and word-topic probability */
using namespace pdb;

class LDADocWordTopicJoin : public JoinComp<LDADocWordTopicAssignment,
                                            LDADocument,
                                            IntDoubleVectorPair,
                                            LDATopicWordProb> {

private:
    Handle<Vector<char>> myMem;
    unsigned numWords;

public:
    ENABLE_DEEP_COPY

    LDADocWordTopicJoin() {}

    LDADocWordTopicJoin(unsigned numWords) : numWords(numWords) {

	/* Set up the random number generator */
        /* Set up the gsl_rng *src */
        gsl_rng* src = gsl_rng_alloc(gsl_rng_mt19937);
        std::random_device rd;
        std::mt19937 gen(rd());
        gsl_rng_set(src, gen());

        /* Allocate space needed for myRand */
        int spaceNeeded = sizeof(gsl_rng) + src->type->size;
        myMem = makeObject<Vector<char>>(spaceNeeded, spaceNeeded);

        /* Copy src over */
        memcpy(myMem->c_ptr(), src, sizeof(gsl_rng));
        memcpy(myMem->c_ptr() + sizeof(gsl_rng), src->state, src->type->size);

        gsl_rng_free(src);
    }

    /* Join condition */
    Lambda<bool> getSelection(Handle<LDADocument> doc,
                              Handle<IntDoubleVectorPair> DocTopicProb,
                              Handle<LDATopicWordProb> WordTopicProb) override {
        return (makeLambdaFromMethod(doc, getDoc) ==
                makeLambdaFromMethod(DocTopicProb, getUnsigned)) &&
            (makeLambdaFromMethod(doc, getWord) == makeLambdaFromMethod(WordTopicProb, getKey));
    }

    Lambda<Handle<LDADocWordTopicAssignment>> getProjection(
        Handle<LDADocument> doc,
        Handle<IntDoubleVectorPair> DocTopicProb,
        Handle<LDATopicWordProb> WordTopicProb) override {

        return makeLambda(
            doc,
            DocTopicProb,
            WordTopicProb,
            [&](Handle<LDADocument>& doc,
                Handle<IntDoubleVectorPair>& DocTopicProb,
                Handle<LDATopicWordProb>& WordTopicProb) {
                int size = (DocTopicProb->getVector()).size();

                /* Compute the posterior probailities of the word coming from each topic */
                double* myProb = new double[size];
                double* topicProbs = DocTopicProb->getVector().c_ptr();
                double* wordProbs = WordTopicProb->getVector().c_ptr();
                for (int i = 0; i < size; ++i) {
                    myProb[i] = topicProbs[i] * wordProbs[i];
                }

                unsigned* topics = new unsigned[size]{0};

                /* Sample the topics (multinomial sampling) */
                gsl_rng* rng = getRng();

                unsigned counts = doc->getCount();
                double* random_values = new double[counts];

                double sum = 0.0;
                for (int i = 0; i < size; ++i) {
                    sum += myProb[i];
                }

                for (int i = 0; i < counts; ++i) {
                    random_values[i] = gsl_rng_uniform(rng) * sum;
                }

                std::sort(random_values, random_values + counts);

                int j = 0;
                double accumuProb = 0.0;
                for (int i = 0; i < size; i++) {
                    accumuProb += myProb[i];
                    while ((j < counts) && (random_values[j] < accumuProb)) {
                        topics[i]++;
                        j++;
                    }
                }

                /* Get the container for the return value */
                Handle<LDADocWordTopicAssignment> retVal = makeObject<LDADocWordTopicAssignment>();
                retVal->setup();
                LDADocWordTopicAssignment& myGuy = *retVal;

                unsigned myDoc = doc->getDoc();
                unsigned myWord = doc->getWord();

                /* Create the doc assignment and topic assignment from the sampled topics, and put them in the return value */
                for (int i = 0; i < size; ++i) {
                    if (topics[i] != 0) {
                        Handle<DocAssignment> whichDoc =
                            makeObject<DocAssignment>(size, myDoc, i, topics[i]);
                        Handle<TopicAssignment> whichTopic =
                            makeObject<TopicAssignment>(numWords, i, myWord, topics[i]);
                        myGuy.push_back(whichDoc);
                        myGuy.push_back(whichTopic);
                    }
                }

                delete[] myProb;
                delete[] topics;
                delete[] random_values;

                return retVal;
            });
    }

    /* Get the GSL RNG from myMem */
    gsl_rng* getRng() {
        gsl_rng* dst = (gsl_rng*)myMem->c_ptr();
        dst->state = (void*)(myMem->c_ptr() + sizeof(gsl_rng));
        dst->type = gsl_rng_mt19937;
        return dst;
    }
};

#endif

ext: h

max_stars_repo_name: dimitrijejankov/pdb

max_issues_repo_name: dimitrijejankov/pdb

max_forks_repo_name: dimitrijejankov/pdb


text:
#pragma once

#include <cuda/runtime_api.hpp>

#include <gsl-lite/gsl-lite.hpp>

#include <cub/cub.cuh>

#include <thrustshift/not-a-vector.h>

namespace thrustshift {

namespace async {

template <class ValuesInRange,
          class ValuesOutRange,
          class ScanOp,
          class MemoryResource>
void inclusive_scan(cuda::stream_t& stream,
                    ValuesInRange&& values_in,
                    ValuesOutRange&& values_out,
                    ScanOp scan_op,
                    MemoryResource& delayed_memory_resource) {

	const std::size_t N = values_in.size();

	gsl_Expects(values_out.size() == N);

	size_t tmp_bytes_size = 0;
	void* tmp_ptr = nullptr;

	auto exec = [&] {
		cuda::throw_if_error(cub::DeviceScan::InclusiveScan(tmp_ptr,
		                                                    tmp_bytes_size,
		                                                    values_in.data(),
		                                                    values_out.data(),
		                                                    scan_op,
		                                                    N,
		                                                    stream.handle()));
	};
	exec();
	auto tmp =
	    make_not_a_vector<uint8_t>(tmp_bytes_size, delayed_memory_resource);
	tmp_ptr = tmp.to_span().data();
	exec();
}

} // namespace async

} // namespace thrustshift

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
/*
** Read input data into array
**
** G.Lohmann, MPI-KYB, 2018
*/
#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

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

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

/*
** get global mean regressor
*/
void GlobalMean(gsl_matrix *Data,gsl_matrix *covariates,int column)
{
  int i=0,j=0;
  int nt = Data->size2;

  double sum=0,nx=(double)Data->size1;
  for (i=0; i<nt; i++) {
    for (j=0; j<Data->size1; j++) {
      sum += gsl_matrix_get(Data,j,i);
    }
    gsl_matrix_set(covariates,i,column,sum/nx);
  }

  /* normalize */
  double u=0,s1=0,s2=0,mx=(double)Data->size2;
  for (i=0; i<nt; i++) {
    u = gsl_matrix_get(covariates,i,column);
    s1 += u;
    s2 += u*u;
  }
  double mean = s1/mx;
  double var = (s2 - mx * mean * mean) / (mx - 1.0);
  double sd = sqrt(var);

  for (i=0; i<nt; i++) {
    u = gsl_matrix_get(covariates,i,column);
    gsl_matrix_set(covariates,i,column,(u-mean/sd));
  }
}

/* normalize time courses */
void VRowNormalize(gsl_matrix *Data)
{
  size_t i,j;
  double u=0,s1=0,s2=0,mean=0,sd=0;
  double nx=(double)Data->size2;

  for (i=0; i<Data->size1; i++) {
    s1 = s2 = 0;
    for (j=0; j<Data->size2; j++) {
      u = gsl_matrix_get(Data,i,j);
      s1 += u;
      s2 += u*u;
    }
    mean = s1/nx;
    sd = sqrt((s2 - nx * mean * mean) / (nx - 1.0));

    for (j=0; j<Data->size2; j++) {
      u = gsl_matrix_get(Data,i,j);
      gsl_matrix_set(Data,i,j,(u-mean)/sd);
    }
  }
}

/* get timing information */
void VGetTimeInfos(VAttrList *list,int nlists,double *mtr,float *run_duration)
{
  VImage src=NULL;
  VAttrListPosn posn;
  int k,ntimesteps=0;
  VLong itr=0;
  double tr=-1.0;

  int ntt=0;
  for (k=0; k<nlists; k++) {
    ntt = 0;
    for (VFirstAttr (list[k], & posn); VAttrExists (& posn); VNextAttr (& posn)) {
      if (VGetAttrRepn (& posn) != VImageRepn) continue;
      VGetAttrValue (& posn, NULL,VImageRepn, & src);
      if (tr < 0) {
	if (VGetAttr(VImageAttrList(src),"repetition_time",NULL,VLongRepn,&itr) != VAttrFound)
	  VError(" TR info missing in header");
	tr = (double) itr / 1000.0;
      }
      ntt = VImageNBands(src);
      run_duration[k] = tr*(float)ntt;
    }
    ntimesteps += ntt;
  }
  if (tr < 0) VError(" TR info missing");
  *mtr = tr;
}

/* read data block */
gsl_matrix *VReadImageData(VAttrList *list,int nlists)
{
  VAttrListPosn posn;
  int i=0,k,ntt=0,slice,row,col,nrows=0,ncols=0,ntimesteps=0;
  double u=0,s1=0,s2=0,nx=0,sd=0,mean=0,tiny=1.0e-8;

  int nslices = VAttrListNumImages(list[0]);

  /* alloc src image */
  VImage **src = (VImage **) VCalloc(nlists,sizeof(VImage *));
  int *nt = (int *) VCalloc(nlists,sizeof(int));
  for (k=0; k<nlists; k++) {
    src[k] = (VImage *) VCalloc(nslices,sizeof(VImage));
    i = 0;
    for (VFirstAttr (list[k], & posn); VAttrExists (& posn); VNextAttr (& posn)) {
      if (VGetAttrRepn (& posn) != VImageRepn) continue;
      if (i >= nslices) VError(" inconsistent number of slices, %d %d %d",k,i,nslices);
      VGetAttrValue (& posn, NULL,VImageRepn, & src[k][i]);
      i++;
    }
  }

  /* get image dimensions */
  nrows = VImageNRows(src[0][0]);
  ncols = VImageNColumns(src[0][0]);
  ntimesteps = 0;
  for (k=0; k<nlists; k++) {
    ntt = VImageNBands(src[k][0]);
    ntimesteps += ntt;
    nt[k] = ntt;
    if (nrows != VImageNRows(src[k][0])) VError(" inconsistent number of rows ");
    if (ncols != VImageNColumns(src[k][0])) VError(" inconsistent number of columns ");
  }

  /* number of nonzero voxels */
  size_t nvox=0;
  for (slice=0; slice<nslices; slice++) {
    for (row=0; row<nrows; row++) {
      for (col=0; col<ncols; col++) {
	u = VGetPixel(src[0][slice],0,row,col);
	if (gsl_isnan(u) || gsl_isinf(u)) continue;
	if (fabs(u) < TINY) continue;
	nvox++;
      }
    }
  }

  /* alloc new data struct */
  gsl_matrix *Data = gsl_matrix_calloc(nvox,ntimesteps);
  if (!Data) VError(" error allocating Data");

  /* fill map and data */
  for (k=0; k<nlists; k++) {

    ntt = 0;
    for (i=0; i<k; i++) ntt += nt[i];

    nvox=0;
    for (slice=0; slice<nslices; slice++) {
      for (row=0; row<nrows; row++) {
	for (col=0; col<ncols; col++) {
	  u = VGetPixel(src[0][slice],0,row,col);
	  if (gsl_isnan(u) || gsl_isinf(u)) continue;
	  if (fabs(u) < TINY) continue;

	  s1=s2=nx=0;
	  for (i=0; i<nt[k]; i++) {
	    u = VGetPixel(src[k][slice],i,row,col);
	    s1 += u;
	    s2 += u*u;
	    nx++;
	  }
	  mean = s1/nx;
	  sd = sqrt((s2 - nx * mean * mean) / (nx - 1.0));
	  if (sd < tiny) {
	    u = 0.0;
	    for (i=0; i<nt[k]; i++) {
	      gsl_matrix_set(Data,nvox,ntt+i,u);
	    }
	  }
	  else {
	    for (i=0; i<nt[k]; i++) {
	      u = VGetPixel(src[k][slice],i,row,col);
	      u = (u-mean)/sd;
	      gsl_matrix_set(Data,nvox,ntt+i,u);
	    }
	  }
	  nvox++;
	}
      }
    }
  }
  VFree(nt);
  return Data;
}

/* get map of voxel addresses */
VImage VoxelMap(VAttrList list)
{
  VAttrListPosn posn;
  int i=0,slice,row,col,nrows=0,ncols=0;

  int nslices = VAttrListNumImages(list);

  /* alloc src image */
  i = nrows = ncols = 0;
  VImage *src = (VImage *) VCalloc(nslices,sizeof(VImage));
  i = 0;
  for (VFirstAttr (list, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
    if (VGetAttrRepn (& posn) != VImageRepn) continue;
    if (i >= nslices) VError(" inconsistent number of slices, %d %d",i,nslices);
    VGetAttrValue (& posn, NULL,VImageRepn, & src[i]);
    if (VImageNRows(src[i])  > nrows)  nrows = VImageNRows(src[i]);
    if (VImageNColumns(src[i]) > ncols) ncols = VImageNColumns(src[i]);
    i++;
  }

  /* number of nonzero voxels */
  double u=0;
  size_t nvox=0;
  for (slice=0; slice<nslices; slice++) {
    for (row=0; row<nrows; row++) {
      for (col=0; col<ncols; col++) {
	u = VGetPixel(src[slice],0,row,col);
	if (fabs(u) < TINY) continue;
	nvox++;
      }
    }
  }
  if (nvox < 1) VError(" No non-zero voxels in input image");

  /* voxel addresses */
  VImage map = VCreateImage(1,4,nvox,VShortRepn);
  if (map == NULL) VError(" error allocating addr map");
  VFillImage(map,VAllBands,0);
  VCopyImageAttrs (src[0],map);
  VSetAttr(VImageAttrList(map),"nvoxels",NULL,VLongRepn,(VLong)nvox);
  VSetAttr(VImageAttrList(map),"nslices",NULL,VLongRepn,(VLong)nslices);
  VSetAttr(VImageAttrList(map),"nrows",NULL,VLongRepn,(VLong)nrows);
  VSetAttr(VImageAttrList(map),"ncols",NULL,VLongRepn,(VLong)ncols);
  VPixel(map,0,3,0,VShort) = nslices;
  VPixel(map,0,3,1,VShort) = nrows;
  VPixel(map,0,3,2,VShort) = ncols;

  /* fill map */
  nvox=0;
  for (slice=0; slice<nslices; slice++) {
    for (row=0; row<nrows; row++) {
      for (col=0; col<ncols; col++) {
	u = VGetPixel(src[slice],0,row,col);
	if (fabs(u) < TINY) continue;
	VPixel(map,0,0,nvox,VShort) = slice;
	VPixel(map,0,1,nvox,VShort) = row;
	VPixel(map,0,2,nvox,VShort) = col;
	nvox++;
      }
    }
  }
  return map;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#ifndef VKST_PLAT_FS_NOTIFY_LINUX_H
#define VKST_PLAT_FS_NOTIFY_LINUX_H

#ifndef VKST_PLAT_FS_NOTIFY_H
#error "Include fs_notify.h only"
#endif

#include "fs_notify.h"
#include <gsl.h>
#include <vector>

namespace plat {

class fs_notify final : public impl::fs_notify<fs_notify> {
public:

private:
  //std::vector<gsl::unique_ptr<watch>> _watches;

  watch_id do_add(plat::filesystem::path path,
                  impl::fs_notify<fs_notify>::notify_delegate delegate,
                  bool recursive, std::error_code& ec) noexcept;

  void do_remove(watch_id id) noexcept;

  void do_tick() noexcept;

  friend class impl::fs_notify<fs_notify>;
}; // class fs_notify

} // namespace plat

#endif // VKST_PLAT_FS_NOTIFY_LINUX_H

max_stars_repo_name: wesleygriffin/vkst

max_issues_repo_name: wesleygriffin/vkst

max_forks_repo_name: wesleygriffin/vkst


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

/*
Adapted from C code from Darren Wilkinson's blog
https://darrenjw.wordpress.com/2011/07/16/gibbs-sampler-in-various-languages-revisited/
To build:
gcc -O3 gibbs_gsl.c -I/usr/local/include -L/usr/local/lib -lgsl -lgslcblas -lm -o gibbs_gsl
To run:
time ./gibbs_gsl
*/

/* Arrays xa and ya are passed by reference */
void gibbs(int N, int thin, double **xa, double **ya)
{
  int i,j;
  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
  double x=0;
  double y=0;
  for (i=0;i<N;i++) {
    for (j=0;j<thin;j++) {
      x=gsl_ran_gamma(r,3.0,1.0/(y*y+4));
      y=1.0/(x+1)+gsl_ran_gaussian(r,1.0/sqrt(2*x+2));
    }
    (*xa)[i]=x; /* Note how the arrays are dereferenced */
    (*ya)[i]=y;
  }
  free(r);
}

int main()
{
  int N=50000;
  int thin=1000;
  double *xa,*ya;
  xa = (double *)malloc(sizeof(double)*N);
  ya = (double *)malloc(sizeof(double)*N);
  gibbs(N,thin,&xa,&ya);
  return(0);
}

max_stars_repo_name: ruqianl/distributions

max_issues_repo_name: ruqianl/distributions

max_forks_repo_name: ruqianl/distributions


text:
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

// y' = exp(x)
// https://www.wolframalpha.com/input/?i=y%27+%3D+2y+%2B+x%2C+y%280%29+%3D+1%2C+y%281%29
// https://www.wolframalpha.com/input/?i=1%2F4+*+%285+*+e%5E2+-+3%29
int odefunc (double x, const double y[], double f[], void *params)
{
    // f[0] = x+2*y[0];
	// pow(): 멱급수(a, b) = a^b (a의 b승)
	// math.h에 있으며 pow를 활용하기 위해서는 -lm 옵션이 추가로 필요합니다.
	// GSL 라이브러리를 활용하기 위해선 -lgsl 옵션이 필요합니다.
	f[0] = pow(M_E, x);

    return GSL_SUCCESS;
}

int * jac;

int main(void)
{
    int dim = 1;
	// 미방을 풀기 위한 데이터 타입
    gsl_odeiv2_system sys = {odefunc, NULL, dim, NULL};

	// 드라이버는 미방을 쉽게 풀 수 있도록 시간에 따른 변화, 제어, 스텝에 대한 래퍼
	// gsl_odeiv2_step_rkf45는 룬게쿠타 4-5차를 의미함
	// 다음 인자는 미분방정식의 시작지점
	// 절대적인 에러 바운드와 상대적인 에러 바운드
    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rkf45, 1e-6, 1e-6, 0.0);
    int i;
    double x0 = 0.0,  xf = 10.0; /* start and end of integration interval */
    double x = x0;
    double y[1] = { 1 };  /* initial value */
	double tmp = 0.0;

    for (i = 0; tmp <= xf; tmp += 0.1)
    {
		// for문을 포함하여 아래 한 줄을 통해 상황에 따라
		// 간소화된 형태로 결과를 출력할 수도 있고
		// 시간이 좀 더 오래걸리더라도 더 정확한 데이터를 산출할 수도 있음
        double xi = x0 + tmp * (xf-x0) / xf;
        int status = gsl_odeiv2_driver_apply (d, &x, xi, y);

        if (status != GSL_SUCCESS)
        {
            printf ("error, return value=%d\n", status);
            break;
        }

        printf ("%.8e %.8e\n", x, y[0]);
    }

    gsl_odeiv2_driver_free (d);
    return 0;
}


text:
/**
* Subroutines for fitting functions to data.
 * Most of the functions are used for the
 * background determination.
 */

#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_multifit.h>

#include "spce_fitting.h"

/**
 * Function: det_vector_polyN
 * A function using GSL fitting functions to fit a set x,y,w of n elements
 * such that an Nth degree polynomial is fitted to the x and y vectors
 * using the weights w. This function avoids NaN values in ys AND elements
 * with an associated weight that is zero.
 *
 * Parameters:
 * @param m   - order of the fit
 * @param xs  - double vector containing the x values
 * @param ys  - double vector containing the y values
 * @param ws  - double vector containing the weights associated with ys
 * @param n   - number of points in xs,ys, and ws (must be greater than m!)
 * @param c   - the vector with the fitted coefficients
 * @param cov - the covariance matrix
 *
 * Returns:
 * @return interp - vector with additional fitting information
 */
gsl_vector *
det_vector_polyN (int m, const double *const xs, double *const ys,
		   double *const ws, const int n, gsl_vector *c,
		   gsl_matrix *cov)
{
  int i, j, nn, ii;
  double chisq;
  gsl_matrix *X;
  gsl_vector *y, *w;

  double *tmp, *xmp;
  double median, xean;

  gsl_vector *interp;

  // allocate the return vector
  interp = gsl_vector_alloc(5);

  // allocate temporary vectors
  tmp = (double *) malloc (n * sizeof (double));
  xmp = (double *) malloc (n * sizeof (double));

  // fill weights and independent values
  // for the background pixels into the arrays
  nn = 0;
  for (i = 0; i < n; i++)
    {
      if (ws[i] > 0.0)
	{
	  //	  tmp[nn] = ws[i];
	  tmp[nn] = 1.0/(ws[i]*ws[i]);
	  xmp[nn] = xs[i];
	  nn++;
	}
    }

  // compute the median weight
  // to be used for all pixels
  gsl_sort( tmp, 1, m);
  median = gsl_stats_median_from_sorted_data(tmp, 1, nn);

  // determine the anchor point in the
  // independent variable
  xean = gsl_stats_mean (xmp, 1, nn);

  // check whether the desired
  // intepolation degree is feasible
  if (nn < m)
    {
      aXe_message (aXe_M_WARN4, __FILE__, __LINE__,
		   "Not enough points (%d)  found to perform the %d order fit. Doing %d order instead.\n",
		   nn, m, nn);
      m = nn;
    }

  // allocate more temporary space
  X = gsl_matrix_alloc (nn, m);
  y = gsl_vector_alloc (nn);
  w = gsl_vector_alloc (nn);
  //  c = gsl_vector_alloc (m);
  //  cov = gsl_matrix_alloc (m, m);

  // transfer independent/dependent and weight
  // values from the background pixels to the
  // temporary arrays
  ii = 0;
  for (i = 0; i < n; i++)
    {
      if (ws[i] > 0.0)
	{
	  // shift the x-values
	  // and store them in the matrix
	  for (j = 0; j < m; j++)
	    {
	      gsl_matrix_set (X, ii, j, pow (xs[i] - xean, j));
	    }
	  gsl_vector_set (y, ii, ys[i]);
	  gsl_vector_set (w, ii, median);
	  ii++;
	}

    }

  // allocate space and do the fit; release the space
  gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (nn, m);
  gsl_multifit_wlinear (X, w, y, c, cov, &chisq, work);
  gsl_multifit_linear_free (work);

  // release memory
  gsl_matrix_free (X);
  gsl_vector_free (y);
  gsl_vector_free (w);

  // free the tmp-arrays
  free(tmp);
  free(xmp);

  // fill the mean x-value
  // and the order into the return vector
  gsl_vector_set(interp, 0, xean);
  gsl_vector_set(interp, 1, (double)m);

  return interp;
}

/**
 * Function: fill_polyN_interp
 * Evaluates a polynomial of any order and fills the
 * result plus the error, compute using the covariance
 * matrix, in vectors.
 *
 * Parameters:
 * @param xs     - double vector containing the x values
 * @param ys     - double vector containing the y values
 * @param yi     - double vector containing the intermediate y-values
 * @param ws     - double vector containing the weights associated with ys
 * @param n      - number of points in xs,ys, and ws (must be greater than m!)
 * @param c      - the vector with the fitted coefficients
 * @param cov    - the covariance matrix
 * @param interp - additional fitting ifnormation
 * @param final  - flagg to indicate final fit and
 *                 therefore a replacement
 */
void
fill_polyN_interp(const double *const xs, double *const ys,
		  double *const ws, double *yi, const int n,
		  gsl_vector *coeffs, gsl_matrix *cov,
		  gsl_vector *interp, const int final)
{
  int i, j, m;
  double xf, yf, sq_yf_err, yf_err;
  double xean;

  // get the mean x-value and
  // polynomial order from the vector
  xean = gsl_vector_get(interp, 0);
  m    = (int)gsl_vector_get(interp, 1);

  // put the interpolated values
  // for object pixels back into the
  // inut arrays
  for (i = 0; i < n; i++)
    {

      // compute and store interpolated
      // values an corresponding error
      xf = xs[i]-xean;
      yf = 0.0;
      yf_err = 0.0;
      sq_yf_err = 0.0;

      // compute the polynimials
      for (j = 0; j < m; j++)
	{
	  // first compute the y-value
	  yf += gsl_vector_get(coeffs,j)*pow (xf, j);

	  // compute the associated error
	  sq_yf_err += gsl_matrix_get(cov,j,j) * pow (xf, j) * pow (xf, j);
	}

      // put the intepolated value in the intermediate vector
      yi[i] = yf;

      // if requested, fill the interpolated values
      // into the origina; arrays
      if (final && ws[i] == 0.0)
	{
	  ys[i] = yf;
	  ws[i] = sqrt(sq_yf_err);
	}
    }
}

/**
 * Function: comp_vector_average
 * The function provides the computation of the average of
 * the y-values and store that average in temporary
 * or the original y-vector.
 *
 * @param xs    - absissa
 * @param ys    - value at xs[]
 * @param ws    - weight in ys[]
 * @param yi    - temprary y-values
 * @param n     - number of points in xs, ys, ws
 * @param final - flagg to store to temp (0) or final (1) vectors
 */
void
comp_vector_average (const double *xs, double *ys,
		     double *ws, double *yi, const int n, const int final)
{
  double ave, std;
  //int i, m;

  det_vector_average (xs, ys, ws, n, &ave, &std);
  fill_const_value(ys, ws, yi, n, ave, std, final);
}

/**
 * Function: comp_vector_median
 * The function provides the computation of the median of
 * the y-values and store that median in temporary
 * or the original y-vector
 *
 * @param xs - absissa
 * @param ys - value at xs[]
 * @param ws - weight in ys[]
 * @param yi    - temprary y-values
 * @param n     - number of points in xs, ys, ws
 * @param final - flagg to store to temp (0) or final (1) vectors
 */
void
comp_vector_median (const double *xs, double *ys,
		    double *ws, double *yi, const int n, const int final)
{
  double med, std;

  det_vector_median (xs, ys, ws, n, &med, &std);
  fill_const_value(ys, ws, yi, n, med, std, final);

}

/**
 * Function: comp_vector_linear
 * The function fits a linear function to the y-values
 * and fills the linear inteprolated values in a temporary
 * vector and, if requested, in the original vector as well.
 *
 * @param xs    - absissa
 * @param ys    - value at xs[]
 * @param ws    - weight in ys[]
 * @param yi    - temprary y-values
 * @param n     - number of points in xs, ys, ws
 * @param final - flagg to store to temp (0) or final (1) vectors
 */
void
comp_vector_linear (const double *xs, double *ys,
		    double *ws, double *yi, const int n, const int final)
{
  gsl_vector *interp;

  // make a weighted linear fit
  interp = det_vector_linear(xs, ys, ws, n, 1);

  // make the interpolations,
  // using the linear fit
  fill_linear_interp(xs, ys, ws, yi, n, interp, final);

  gsl_vector_free(interp);
}

/**
 * Function: comp_vector_polyN
 * The function fits a polynomial function to the y-values
 * and fills the evaluated polynomial values in a temporary
 * vector and, if requested, in the original vector as well.
 *
 * @param xs    - absissa
 * @param ys    - value at xs[]
 * @param ws    - weight in ys[]
 * @param yi    - temprary y-values
 * @param n     - number of points in xs, ys, ws
 * @param final - flagg to store to temp (0) or final (1) vectors
 */
void
comp_vector_polyN (const int m, const double *xs, double *ys,
		    double *ws, double *yi, const int n, const int final)
{
  gsl_vector *interp, *coeffs;
  gsl_matrix *cov;

  coeffs = gsl_vector_alloc(m);
  cov = gsl_matrix_alloc(m,m);

  interp = det_vector_polyN (m, xs, ys, ws, n, coeffs, cov);
  fill_polyN_interp(xs, ys, ws, yi, n, coeffs, cov, interp, final);

  gsl_vector_free(interp);
  gsl_vector_free(coeffs);
  gsl_matrix_free(cov);
}

/**
 * Function: det_vector_linear
 * Performs a linear fit of the xs, ys, and ws arrays.
 * Ignores NaN values. Fitted values are returned in
 * ys and errors in ws.
 *
 * Parameters:
 * @param xs     - absissa
 * @param ys     - value at xs[]
 * @param ws     - error in ys[]
 * @param n      - number of points in xs, ys, ws
 * @param weight - make weighted fit (=1) or not (=0)
 *
 * Returns:
 * @return ret - the covariance values
 */
gsl_vector *
det_vector_linear(const double *xs, double *ys, double *ws,
		  const int n, const int weight)
{
  int i, m;
  double c0, c1, cov00, cov01, cov11, chisq;
  double *xx, *yy, *ww, *tmp, *xmp;
  double median, xean;
  //double xf, yf, yf_err;

  gsl_vector *ret;

  ret = gsl_vector_alloc(10);

  // allocate space for temporary vectors
  tmp = (double *) malloc (n * sizeof (double));
  xmp = (double *) malloc (n * sizeof (double));

  // fill the temporary vectors
  // with background values and weights
  m = 0;
  for (i = 0; i < n; i++)
    {
      if (ws[i] > 0.0)
	{
	  //	  tmp[m] = ws[i];
	  tmp[m] = 1.0/(ws[i]*ws[i]);
	  xmp[m] = xs[i];
	  m++;
	}
    }

  // determine the median weight which
  // is used for all data
  gsl_sort( tmp, 1, m);
  median = gsl_stats_median_from_sorted_data(tmp, 1, m);

  // determine the mean in the independent variable
  xean = gsl_stats_mean (xmp, 1, m);

  // allocate more temporary vectors
  xx = (double *) malloc (m * sizeof (double));
  yy = (double *) malloc (m * sizeof (double));
  ww = (double *) malloc (m * sizeof (double));

  // fill the temporary vectors with
  // background independent/dependent variables
  // plus the weight
  m = 0;
  for (i = 0; i < n; i++)
    {
      if (ws[i] > 0.0)
	{
	 xx[m] = xs[i]-xean;
	 yy[m] = ys[i];
	 ww[m] = median;
	 m++;
	}
    }

  // check for weighted fit
  if (weight)
    // doe the linear fit
    gsl_fit_wlinear (xx, 1, ww, 1, yy, 1, m, &c0, &c1, &cov00,
		     &cov01, &cov11, &chisq);
  else
    // does the linear fit
    gsl_fit_linear (xx, 1, yy, 1, m, &c0, &c1, &cov00,
		    &cov01, &cov11, &chisq);

  gsl_vector_set(ret, 0, xean);
  gsl_vector_set(ret, 1, c0);
  gsl_vector_set(ret, 2, c1);
  gsl_vector_set(ret, 3, cov00);
  gsl_vector_set(ret, 4, cov01);
  gsl_vector_set(ret, 5, cov11);
  gsl_vector_set(ret, 6, chisq);

  // free the arrays
  free(tmp);
  free(xmp);
  free(xx);
  free(yy);
  free(ww);

  return ret;
}

/**
 * Function: fill_linear_interp
 * The function computes and stores linear interpolated values
 * in one or several output vectors.
 *
 * Parameters:
 * @param xs     - double vector containing the x values
 * @param ys     - double vector containing the y values
 * @param ws     - double vector containing the weights associated with ys
 * @param yi     - double vector containing the intermediate y-values
 * @param n      - number of points in xs,ys, and ws (must be greater than m!)
 * @param interp - vector with the fitted values
 * @param final  - flagg to indicate final fit and
 *                 therefore a replacement
 */
void
fill_linear_interp(const double *const xs, double *const ys, double *const ws,
		   double *yi, const int n, gsl_vector *interp,
		   const int final)
{
  double c0, c1, cov00, cov01, cov11, xean;
  double xf, yf, yf_err;
  int i;

  xean  = gsl_vector_get(interp, 0);
  c0    = gsl_vector_get(interp, 1);
  c1    = gsl_vector_get(interp, 2);
  cov00 = gsl_vector_get(interp, 3);
  cov01 = gsl_vector_get(interp, 4);
  cov11 = gsl_vector_get(interp, 5);

  // determine the interpolated values
  // and store them in the input arrays
  for (i = 0; i < n; i++)
    {

      xf = xs[i]-xean;
      gsl_fit_linear_est (xf, c0, c1, cov00, cov01, cov11, &yf, &yf_err);
      yi[i] = yf;

      // for object pixels,
      // fill in the interpolated values
      if (final && ws[i] == 0.0)
	{
	  ys[i] = yf;
	  ws[i] = yf_err;
	}
    }
}

/**
 * Function: det_vector_average
 * The function computes the average and the standard
 * deviation from the values in a vector
 *
 * @param xs  - absissa
 * @param ys  - value at xs[]
 * @param ws  - weight in ys[]
 * @param n   - number of points in xs, ys, ws
 * @param avg - the average
 * @param std - the standard deviation
 */
void
det_vector_average (const double *xs, double *ys,
		    double *ws, const int n, double *avg, double *std)
{
  double *tmp, sum = 0.0;
  int nn = 0;
  int j;

  // count the number of
  // interpolation points
  nn = 0;
  for (j = 0; j < n; j++)
    {
      if (ws[j] != 0.0)
	nn++;
    }

  // allocate a tmp-array
  tmp = malloc (nn * sizeof (double));

  // sum up all values in the
  // interpolation points,
  // and transfer the values to the
  // tmp-array
  nn = 0;
  for (j = 0; j < n; j++)
    {
      if (ws[j] != 0.0)
	{
	  sum += ys[j];
	  tmp[nn] = ys[j];
	  nn++;
	}
    }

  // compute the average
  if (nn > 0)
    *avg = sum / (double)nn;

  // compute the standard deviation
  *std = gsl_stats_sd (tmp, 1, nn);
}

/**
 * Function: det_vector_median
 * The function computes the median and the standard
 * deviation from the values in a vector
 *
 * @param xs  - absissa
 * @param ys  - value at xs[]
 * @param ws  - weight in ys[]
 * @param n   - number of points in xs, ys, ws
 * @param med - the median
 * @param std - the standard deviation
 */
void
det_vector_median (const double *xs, double *ys,
		   double *ws, const int n, double *med, double *std)
{
  double *tmp;
  int nn = 0;
  int j;

  // count the number of
  // interpolation points
  for (j = 0; j < n; j++)
    {
      if (ws[j] != 0.0)
	  nn++;
    }

  // allocate a tmp-array
  tmp = malloc (nn * sizeof (double));

  // transfer the values of the interpolation
  // points to the tmp-array
  nn = 0;
  for (j = 0; j < n; j++)
    {
      if (ws[j] != 0.0)
	{
	  tmp[nn] = ys[j];
	  nn++;
	}
    }

  // sort the tmp array
  gsl_sort (tmp, 1, nn);

  // compute the median
  *med = gsl_stats_median_from_sorted_data (tmp, 1, nn);

  // compute the standard deviation
  *std = gsl_stats_sd (tmp, 1, nn);
}

/**
 * Function:fill_const_value
 * The function fills two constant values into
 * up to three vectors.
 *
 * @param ys    - value at xs[]
 * @param ws    - weight in ys[]
 * @param yi    - value at xs[]
 * @param n     - number of points in xs, ys, ws
 * @param cval  - number of points in xs, ys, ws
 * @param stdev - number of points in xs, ys, ws
 * @param final - flagg which indicates which vectors to fill
 */
void
fill_const_value(double *ys, double *ws, double *yi, const int n,
		 double cval, double stdev, const int final)
{

  int j;
  int m=0;

  // fill the constant value and the stdev
  // in the whole value and error
  // arrays, respectively
  for (j = 0; j < n; j++)
    {
      yi[j] = cval;

      if (final && ws[j] == 0.0)
	{
	  ys[j] = cval;
	  ws[j] = stdev;
	  m++;
	}
    }
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_elementary.h"

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

int
gsl_sf_multiply_e(const double x, const double y, gsl_sf_result * result)
{
  const double ax = fabs(x);
  const double ay = fabs(y);

  if(x == 0.0 || y == 0.0) {
    /* It is necessary to eliminate this immediately.
     */
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if((ax <= 1.0 && ay >= 1.0) || (ay <= 1.0 && ax >= 1.0)) {
    /* Straddling 1.0 is always safe.
     */
    result->val = x*y;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double f = 1.0 - 2.0 * GSL_DBL_EPSILON;
    const double min = GSL_MIN_DBL(fabs(x), fabs(y));
    const double max = GSL_MAX_DBL(fabs(x), fabs(y));
    if(max < 0.9 * GSL_SQRT_DBL_MAX || min < (f * DBL_MAX)/max) {
      result->val = GSL_COERCE_DBL(x*y);
      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      CHECK_UNDERFLOW(result);
      return GSL_SUCCESS;
    }
    else {
      OVERFLOW_ERROR(result);
    }
  }
}

int
gsl_sf_multiply_err_e(const double x, const double dx,
                         const double y, const double dy,
                         gsl_sf_result * result)
{
  int status = gsl_sf_multiply_e(x, y, result);
  result->err += fabs(dx*y) + fabs(dy*x);
  return status;
}

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

#include "eval.h"

double gsl_sf_multiply(const double x, const double y)
{
  EVAL_RESULT(gsl_sf_multiply_e(x, y, &result));
}


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for functions windowing and computing FFT/IFFT of time/frequency series.
 *
 */

#ifndef __GENERATETDIFD_H__
#define __GENERATETDIFD_H__ 1

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#define _XOPEN_SOURCE 500

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "waveform.h"
#include "fft.h"
#include "LISAgeometry.h"
#include "LISAFDresponse.h"
#include "LISAnoise.h"

/********************************** Structures ******************************************/

typedef enum GenTDIFDtag {
  TDIhlm,
  TDIFD,
  TDIFFT,
  h22FFT,
  yslrFFT
} GenTDIFDtag;

/* Parameters for the generation of a LISA waveform (in the form of a list of modes) */
typedef struct tagGenTDIFDparams {
  double tRef;               /* reference time (s) - GPS time at the frequency representing coalescence */
  double torb;               /* reference orbital time entering orbital response (s) */
  double phiRef;             /* reference phase (rad) - phase at the frequency representing coalescence (or at fRef if specified) */
  double fRef;               /* reference frequency at which phiRef is set (Hz, default 0 which is interpreted as Mf=0.14) */
  double m1;                 /* mass of companion 1 (solar masses, default 2e6) */
  double m2;                 /* mass of companion 2 (solar masses, default 1e6) */
  double distance;           /* distance of source (Mpc, default 1e3) */
  double inclination;        /* inclination of source (rad, default pi/3) */
  double lambda;             /* first angle for the position in the sky (rad, default 0) */
  double beta;               /* second angle for the position in the sky (rad, default 0) */
  double polarization;       /* polarization angle (rad, default 0) */

  int nbmode;                /* number of modes to generate (starting with 22) - defaults to 5 (all modes) */
  double minf;               /* Minimal frequency (Hz, default=0) - when set to 0, use the lowest frequency where the detector noise model is trusted __LISASimFD_Noise_fLow (set somewhat arbitrarily)*/
  double maxf;               /* Maximal frequency (Hz, default=0) - when set to 0, use the highest frequency where the detector noise model is trusted __LISASimFD_Noise_fHigh (set somewhat arbitrarily)*/
  double deltatobs;          /* Observation duration (years, default=2) */
  int tagextpn;              /* Tag to allow PN extension of the waveform at low frequencies (default=1) */
  double Mfmatch;            /* When PN extension allowed, geometric matching frequency: will use ROM above this value. If <=0, use ROM down to the lowest covered frequency (default=0.) */
  int setphiRefatfRef;       /* Flag for adjusting the FD phase at phiRef at the given fRef, which depends also on tRef - if false, treat phiRef simply as an orbital phase shift (minus an observer phase shift) (default=1) */
  double deltaf;             /* When generating frequency series from the mode contributions, deltaf for the output (0 to set automatically at 1/2*1/(2T)) */
  double twindowbeg;         /* When generating frequency series from file by FFT, twindowbeg (0 to set automatically at 0.05*duration) */
  double twindowend;         /* When generating frequency series from file by FFT, twindowend (0 to set automatically at 0.01*duration) */
  int tagh22fromfile;        /* Tag choosing wether to load h22 FD downsampled Amp/Phase from file (default 0) */
  int tagtdi;                /* Tag selecting the desired output format */
  int frozenLISA;            /* Freeze the orbital configuration to the time of peak of the injection (default 0) */
  ResponseApproxtag responseapprox;    /* Approximation in the GAB and orb response - choices are full (full response, default), lowfL (keep orbital delay frequency-dependence but simplify constellation response) and lowf (simplify constellation and orbital response) - WARNING : at the moment noises are not consistent, and TDI combinations from the GAB are unchanged */
  int taggenwave;            /* Tag selecting the desired output format */
  int restorescaledfactor;   /* If 1, restore the factors that were scaled out of TDI observables */
  int FFTfromtdfile;         /* Option for loading time series and FFTing (default: false) */
  int nsamplesinfile;        /* Number of lines of input file */
  int binaryin;              /* Tag for loading the data in gsl binary form instead of text (default false) */
  int binaryout;             /* Tag for outputting the data in gsl binary form instead of text (default false) */
  char indir[256];           /* Path for the input directory */
  char infile[256];          /* Path for the input file */
  char outdir[256];          /* Path for the output directory */
  char outfileprefix[256];   /* Path for the output file */
} GenTDIFDparams;

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _GENERATETDIFD_H */

max_stars_repo_name: JohnGBaker/flare

max_issues_repo_name: JohnGBaker/flare

max_forks_repo_name: JohnGBaker/flare


text:
#pragma once

#include <Babylon/JsRuntime.h>
#include <napi/env.h>
#include <gsl/gsl>
#include <cassert>

namespace Babylon
{
    class NativeDataStream final : public Napi::ObjectWrap<NativeDataStream>
    {
        static constexpr auto JS_CLASS_NAME = "_NativeDataStream";
        static constexpr auto JS_ENGINE_CONSTRUCTOR_NAME = "NativeDataStream";

        static constexpr auto VALIDATION_ENABLED = false;

        enum class ValidationType : uint32_t
        {
            Uint32,
            Int32,
            Float32,
            Uint32Array,
            Int32Array,
            Float32Array,
            NativeData,
            Boolean,
        };

        template<ValidationType T, typename ReaderT>
        static inline void Validate(ReaderT& reader)
        {
            if constexpr (VALIDATION_ENABLED)
            {
                uint32_t value{reader.template Read<uint32_t>()};
                if (value != static_cast<uint32_t>(T))
                {
                    throw std::runtime_error{"Data stream validation error: data type mismatch."};
                }
            }
        }

    public:
        class Reader final
        {
        public:
            Reader(const Reader&) = delete;
            Reader(Reader&&) = delete;
            Reader operator=(const Reader&) = delete;
            Reader operator=(const Reader&&) = delete;

            bool CanRead() const
            {
                assert(m_position <= static_cast<size_t>(m_buffer.size()));
                return m_position < static_cast<size_t>(m_buffer.size());
            }

            uint32_t ReadUint32()
            {
                Validate<ValidationType::Uint32>(*this);
                return Read<uint32_t>();
            }

            int32_t ReadInt32()
            {
                Validate<ValidationType::Int32>(*this);
                return Read<int32_t>();
            }

            float ReadFloat32()
            {
                Validate<ValidationType::Float32>(*this);
                return Read<float>();
            }

            gsl::span<uint32_t> ReadUint32Array()
            {
                Validate<ValidationType::Uint32Array>(*this);
                return ReadArray<uint32_t>();
            }

            gsl::span<int32_t> ReadInt32Array()
            {
                Validate<ValidationType::Int32Array>(*this);
                return ReadArray<int32_t>();
            }

            gsl::span<float> ReadFloat32Array()
            {
                Validate<ValidationType::Float32Array>(*this);
                return ReadArray<float>();
            }

            template<typename T>
            T ReadNativeData()
            {
                Validate<ValidationType::NativeData>(*this);
                static_assert(sizeof(T) % 4 == 0);
                auto span = gsl::make_span(reinterpret_cast<uint32_t*>(m_buffer.data() + m_position), sizeof(T) / 4);
                m_position += sizeof(T) / 4;
                return *reinterpret_cast<T*>(span.data());
            }

            template<typename T, class = typename std::enable_if<!std::is_pointer<T>::value>::type>
            auto ReadPointer()
            {
                return ReadNativeData<typename std::conditional<std::is_member_pointer<T>::value, T, T*>::type>();
            }

        private:
            gsl::span<uint32_t> m_buffer{};
            size_t m_position{0};
            const gsl::final_action<std::function<void()>> m_scopeGuard;

            friend class NativeDataStream;

            template<typename CallableT>
            Reader(gsl::span<uint32_t> buffer, CallableT&& callable)
                : m_buffer{buffer}
                , m_scopeGuard{std::forward<CallableT>(callable)}
            {
            }

            template<typename T>
            T Read()
            {
                static_assert(sizeof(T) % 4 == 0);
                T t{*reinterpret_cast<T*>(m_buffer.data() + m_position)};
                m_position += sizeof(T) / 4;
                return t;
            }

            template<typename T>
            gsl::span<T> ReadArray()
            {
                static_assert(sizeof(T) % 4 == 0);

                // The first 32 bit number is a length
                uint32_t length = Read<uint32_t>();

                auto span = gsl::make_span<T>(reinterpret_cast<T*>(m_buffer.data() + m_position), length * (sizeof(T) / 4));
                m_position += length;
                return span;
            }
        };

        static void Initialize(Napi::Env env)
        {
            Napi::HandleScope scope{env};
            if constexpr (VALIDATION_ENABLED)
            {
                Napi::Function func = DefineClass(
                    env,
                    JS_CLASS_NAME,
                    {
                        InstanceMethod("writeBuffer", &NativeDataStream::WriteBuffer),

                        StaticValue("VALIDATION_ENABLED", Napi::Boolean::From(env, VALIDATION_ENABLED)),
                        StaticValue("VALIDATION_UINT_32", Napi::Number::From(env, static_cast<uint32_t>(ValidationType::Uint32))),
                        StaticValue("VALIDATION_INT_32", Napi::Number::From(env, static_cast<uint32_t>(ValidationType::Int32))),
                        StaticValue("VALIDATION_FLOAT_32", Napi::Number::From(env, static_cast<uint32_t>(ValidationType::Float32))),
                        StaticValue("VALIDATION_UINT_32_ARRAY", Napi::Number::From(env, static_cast<uint32_t>(ValidationType::Uint32Array))),
                        StaticValue("VALIDATION_INT_32_ARRAY", Napi::Number::From(env, static_cast<uint32_t>(ValidationType::Int32Array))),
                        StaticValue("VALIDATION_FLOAT_32_ARRAY", Napi::Number::From(env, static_cast<uint32_t>(ValidationType::Float32Array))),
                        StaticValue("VALIDATION_NATIVE_DATA", Napi::Number::From(env, static_cast<uint32_t>(ValidationType::NativeData))),
                        StaticValue("VALIDATION_BOOLEAN", Napi::Number::From(env, static_cast<uint32_t>(ValidationType::Boolean))),
                    });
                JsRuntime::NativeObject::GetFromJavaScript(env).Set(JS_ENGINE_CONSTRUCTOR_NAME, func);
            }
            else
            {
                Napi::Function func = DefineClass(
                    env,
                    JS_CLASS_NAME,
                    {
                        InstanceMethod("writeBuffer", &NativeDataStream::WriteBuffer)
                    });
                JsRuntime::NativeObject::GetFromJavaScript(env).Set(JS_ENGINE_CONSTRUCTOR_NAME, func);
            }
        }

        NativeDataStream(const Napi::CallbackInfo& info)
            : Napi::ObjectWrap<NativeDataStream>(info)
            , m_requestFlushCallback{Napi::Persistent(info[0].As<Napi::Function>())}
        {
        }

        void WriteBuffer(const Napi::CallbackInfo& info)
        {
            assert(!m_locked); // Cannot write bytes while the stream is locked for reading.

            const auto& buffer = info[0].As<Napi::ArrayBuffer>();
            const auto& length = info[1].ToNumber().Uint32Value();

            auto span = gsl::make_span(reinterpret_cast<uint32_t*>(buffer.Data()), static_cast<ptrdiff_t>(length));
            m_buffer.insert(m_buffer.end(), span.begin(), span.end());
        }

        Reader GetReader()
        {
            assert(!m_locked);
            m_requestFlushCallback.Call({});
            m_locked = true;
            return {m_buffer, [this]() {
                m_buffer.clear();
                m_locked = false;
            }};
        }

    private:
        std::vector<uint32_t> m_buffer{};
        Napi::FunctionReference m_requestFlushCallback{};
        bool m_locked{false};
    };
}

max_stars_repo_name: SergioRZMasson/BabylonNative

max_issues_repo_name: SergioRZMasson/BabylonNative

max_forks_repo_name: SergioRZMasson/BabylonNative


text:
/* this program is in the public domain */

#include "bench-user.h"
#include <math.h>

#ifdef HAVE_GSL_GSL_VERSION_H
#include <gsl/gsl_version.h>
#endif

static const char *mkvers(void)
{
#ifdef HAVE_GSL_GSL_VERSION_H
     return gsl_version;
#else
     return "unknown"
#endif
}

BEGIN_BENCH_DOC
BENCH_DOC("name", "gsl-mixed-radix")
BENCH_DOC("package", "GNU Scientific Library (GSL)")
BENCH_DOC("author", "Brian Gough")
BENCH_DOC("year", "2002")
BENCH_DOC("url", "http://sources.redhat.com/gsl")
BENCH_DOC("url-was-valid-on", "Thu Jul 26 07:48:45 EDT 2001")
BENCH_DOC("copyright",
  "Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough\n"
  "\n"
  "This program is free software; you can redistribute it and/or modify\n"
  "it under the terms of the GNU General Public License as published by\n"
  "the Free Software Foundation; either version 2 of the License, or (at\n"
  "your option) any later version.\n")
BENCH_DOC("language", "C")
BENCH_DOCF("version", mkvers)
END_BENCH_DOC

#ifdef BENCHFFT_SINGLE
#include <gsl/gsl_fft_complex_float.h>
#include <gsl/gsl_fft_real_float.h>
#include <gsl/gsl_fft_halfcomplex_float.h>
#define C(x) gsl_fft_complex_float_##x
#define CWv(x) gsl_fft_complex_wavetable_float_##x
#define CWk(x) gsl_fft_complex_workspace_float_##x
#define R(x) gsl_fft_real_float_##x
#define RWv(x) gsl_fft_real_wavetable_float_##x
#define RWk(x) gsl_fft_real_workspace_float_##x
#define H(x) gsl_fft_halfcomplex_float_##x
#define HWv(x) gsl_fft_halfcomplex_wavetable_float_##x
#else
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>
#define C(x) gsl_fft_complex_##x
#define CWv(x) gsl_fft_complex_wavetable_##x
#define CWk(x) gsl_fft_complex_workspace_##x
#define R(x) gsl_fft_real_##x
#define RWv(x) gsl_fft_real_wavetable_##x
#define RWk(x) gsl_fft_real_workspace_##x
#define H(x) gsl_fft_halfcomplex_##x
#define HWv(x) gsl_fft_halfcomplex_wavetable_##x
#endif

int can_do(struct problem *p)
{
     return (p->rank == 1 && problem_in_place(p));
}

/* fftpack-style real transforms */
void copy_h2c(struct problem *p, bench_complex *out)
{
     copy_h2c_1d_fftpack(p, out, -1.0);
}

void copy_c2h(struct problem *p, bench_complex *in)
{
     copy_c2h_1d_fftpack(p, in, -1.0);
}

static void *Wv = 0, *Wk = 0;

void setup(struct problem *p)
{
     int n;

     BENCH_ASSERT(can_do(p));
     n = p->n[0];

     if (p->kind == PROBLEM_COMPLEX) {
	  Wv = CWv(alloc)(n);
	  Wk = CWk(alloc)(n);
     } else {
	  if (p->sign == -1)
	       Wv = RWv(alloc)(n);
	  else
	       Wv = HWv(alloc)(n);
	  Wk = RWk(alloc)(n);
     }
}

void doit(int iter, struct problem *p)
{
     int i;
     int n = p->n[0];
     void *in = p->in;
     void *wv = Wv, *wk = Wk;

     if (p->kind == PROBLEM_COMPLEX) {
	  if (p->sign == -1) {
	    for (i = 0; i < iter; ++i)
		 C(forward)(in, 1, n, wv, wk);
	  } else {
	    for (i = 0; i < iter; ++i)
		 C(backward)(in, 1, n, wv, wk);
	  }
     } else {
	  if (p->sign == -1) {
	    for (i = 0; i < iter; ++i)
		 R(transform)(in, 1, n, wv, wk);
	  } else {
	    for (i = 0; i < iter; ++i)
		 H(transform)(in, 1, n, wv, wk);
	  }
     }
}

void done(struct problem *p)
{
     UNUSED(p);
     if (p->kind == PROBLEM_COMPLEX) {
	  CWv(free)(Wv);
     } else {
	  if (p->sign == -1)
	       RWv(free)(Wv);
	  else
	       HWv(free)(Wv);
	  RWk(free)(Wk);
     }
}

max_stars_repo_name: mreineck/benchfft

max_issues_repo_name: mreineck/benchfft

max_forks_repo_name: mreineck/benchfft


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

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>

#include <gsl/gsl_ieee_utils.h>

int
main (void)
{
  double y, y_expected;

  gsl_ieee_env_setup ();

  /* Test for expm1 */

  y = gsl_expm1 (0.0); y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(0.0)");

  y = gsl_expm1 (1e-10); y_expected = 1.000000000050000000002e-10;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(1e-10)");

  y = gsl_expm1 (-1e-10); y_expected = -9.999999999500000000017e-11;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(-1e-10)");

  y = gsl_expm1 (0.1); y_expected = 0.1051709180756476248117078264902;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(0.1)");

  y = gsl_expm1 (-0.1); y_expected = -0.09516258196404042683575094055356;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(-0.1)");

  y = gsl_expm1 (10.0); y_expected = 22025.465794806716516957900645284;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(10.0)");

  y = gsl_expm1 (-10.0); y_expected = -0.99995460007023751514846440848444;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(-10.0)");

  /* Test for log1p */

  y = gsl_log1p (0.0); y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_log1p(0.0)");

  y = gsl_log1p (1e-10); y_expected = 9.9999999995000000000333333333308e-11;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_log1p(1e-10)");

  y = gsl_log1p (0.1); y_expected = 0.095310179804324860043952123280765;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_log1p(0.1)");

  y = gsl_log1p (10.0); y_expected = 2.3978952727983705440619435779651;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_log1p(10.0)");

  /* Test for gsl_hypot */

  y = gsl_hypot (0.0, 0.0) ; y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(0.0, 0.0)");

  y = gsl_hypot (1e-10, 1e-10) ; y_expected = 1.414213562373095048801688e-10;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e-10, 1e-10)");

  y = gsl_hypot (1e-38, 1e-38) ; y_expected = 1.414213562373095048801688e-38;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e-38, 1e-38)");

  y = gsl_hypot (1e-10, -1.0) ; y_expected = 1.000000000000000000005;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e-10, -1)");

  y = gsl_hypot (-1.0, 1e-10) ; y_expected = 1.000000000000000000005;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(-1, 1e-10)");

  y = gsl_hypot (1e307, 1e301) ; y_expected = 1.000000000000499999999999e307;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e307, 1e301)");

  y = gsl_hypot (1e301, 1e307) ; y_expected = 1.000000000000499999999999e307;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e301, 1e307)");

  y = gsl_hypot (1e307, 1e307) ; y_expected = 1.414213562373095048801688e307;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e307, 1e307)");

  /* Test for acosh */

  y = gsl_acosh (1.0); y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_acosh(1.0)");

  y = gsl_acosh (1.1); y_expected = 4.435682543851151891329110663525e-1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_acosh(1.1)");

  y = gsl_acosh (10.0); y_expected = 2.9932228461263808979126677137742e0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_acosh(10.0)");

  y = gsl_acosh (1e10); y_expected = 2.3718998110500402149594646668302e1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_acosh(1e10)");

  /* Test for asinh */

  y = gsl_asinh (0.0); y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(0.0)");

  y = gsl_asinh (1e-10); y_expected = 9.9999999999999999999833333333346e-11;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(1e-10)");

  y = gsl_asinh (-1e-10); y_expected = -9.9999999999999999999833333333346e-11;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(1e-10)");

  y = gsl_asinh (0.1); y_expected = 9.983407889920756332730312470477e-2;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(0.1)");

  y = gsl_asinh (-0.1); y_expected = -9.983407889920756332730312470477e-2;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(-0.1)");

  y = gsl_asinh (1.0); y_expected = 8.8137358701954302523260932497979e-1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(1.0)");

  y = gsl_asinh (-1.0); y_expected = -8.8137358701954302523260932497979e-1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(-1.0)");

  y = gsl_asinh (10.0); y_expected = 2.9982229502979697388465955375965e0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(10)");

  y = gsl_asinh (-10.0); y_expected = -2.9982229502979697388465955375965e0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(-10)");

  y = gsl_asinh (1e10); y_expected = 2.3718998110500402149599646668302e1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(1e10)");

  y = gsl_asinh (-1e10); y_expected = -2.3718998110500402149599646668302e1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(-1e10)");

  /* Test for atanh */

  y = gsl_atanh (0.0); y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(0.0)");

  y = gsl_atanh (1e-20); y_expected = 1e-20;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(1e-20)");

  y = gsl_atanh (-1e-20); y_expected = -1e-20;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(-1e-20)");

  y = gsl_atanh (0.1); y_expected = 1.0033534773107558063572655206004e-1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(0.1)");

  y = gsl_atanh (-0.1); y_expected = -1.0033534773107558063572655206004e-1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(-0.1)");

  y = gsl_atanh (0.9); y_expected = 1.4722194895832202300045137159439e0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(0.9)");

  y = gsl_atanh (-0.9); y_expected = -1.4722194895832202300045137159439e0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(0.9)");

  /* Test for pow_int */

  y = gsl_pow_2 (-3.14); y_expected = pow(-3.14, 2.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_2(-3.14)");

  y = gsl_pow_3 (-3.14); y_expected = pow(-3.14, 3.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_3(-3.14)");

  y = gsl_pow_4 (-3.14); y_expected = pow(-3.14, 4.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_4(-3.14)");

  y = gsl_pow_5 (-3.14); y_expected = pow(-3.14, 5.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_5(-3.14)");

  y = gsl_pow_6 (-3.14); y_expected = pow(-3.14, 6.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_6(-3.14)");

  y = gsl_pow_7 (-3.14); y_expected = pow(-3.14, 7.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_7(-3.14)");

  y = gsl_pow_8 (-3.14); y_expected = pow(-3.14, 8.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_8(-3.14)");

  y = gsl_pow_9 (-3.14); y_expected = pow(-3.14, 9.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_9(-3.14)");

  {
    int n;
    for (n = -9; n < 10; n++) {
      y = gsl_pow_int (-3.14, n); y_expected = pow(-3.14, n);
      gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_n(-3.14,%d)", n);
    }
  }

  /* Test for isinf, isnan, finite*/

  {
    double zero, one, inf, nan;
    int s;

    zero = 0.0;
    one = 1.0;
    inf = exp(1.0e10);
    nan = inf / inf;

    s = gsl_isinf(zero);
    gsl_test_int (s, 0, "gsl_isinf(0)");

    s = gsl_isinf(one);
    gsl_test_int (s, 0, "gsl_isinf(1)");

    s = gsl_isinf(inf);
    gsl_test_int (s, 1, "gsl_isinf(inf)");

    s = gsl_isinf(-inf);
    gsl_test_int (s, -1, "gsl_isinf(-inf)");

    s = gsl_isinf(nan);
    gsl_test_int (s, 0, "gsl_isinf(nan)");

    s = gsl_isnan(zero);
    gsl_test_int (s, 0, "gsl_isnan(0)");

    s = gsl_isnan(one);
    gsl_test_int (s, 0, "gsl_isnan(1)");

    s = gsl_isnan(inf);
    gsl_test_int (s, 0, "gsl_isnan(inf)");

    s = gsl_isnan(nan);
    gsl_test_int (s, 1, "gsl_isnan(nan)");

    s = gsl_finite(zero);
    gsl_test_int (s, 1, "gsl_finite(0)");

    s = gsl_finite(one);
    gsl_test_int (s, 1, "gsl_finite(1)");

    s = gsl_finite(inf);
    gsl_test_int (s, 0, "gsl_finite(inf)");

    s = gsl_finite(nan);
    gsl_test_int (s, 0, "gsl_finite(nan)");
  }

  {
    double x = gsl_fdiv (2.0, 3.0);
    gsl_test_rel (x, 2.0/3.0, 4*GSL_DBL_EPSILON, "gsl_fdiv(2,3)");
  }

  exit (gsl_test_summary ());
}


text:
#ifndef _OPTIMIZE_STATE_RECONSTRUCTOR_GSL_H_
#define _OPTIMIZE_STATE_RECONSTRUCTOR_GSL_H_

using namespace std;

#include <armadillo>
using namespace arma;

#include "rate_model.h"
#include "state_reconstructor.h"

#include <gsl/gsl_vector.h>

class OptimizeStateReconstructor{
private:
    RateModel * rm;
    StateReconstructor * sr;
    mat * free_variables;
    int nfree_variables;
    int maxiterations;
    double stoppingprecision;

    double GetLikelihoodWithOptimized(const gsl_vector * variables);
    static double GetLikelihoodWithOptimized_gsl(const gsl_vector * variables, void *obj);

public:
    OptimizeStateReconstructor(RateModel *,StateReconstructor *, mat *, int);
    mat optimize();
};

#endif /* _OPTIMIZE_STATE_RECONSTRUCTOR_GSL_H_ */


text:
/*****************************************************************************
*
* Rokko: Integrated Interface for libraries of eigenvalue decomposition
*
* Copyright (C) 2015 by Rokko Developers https://github.com/t-sakashita/rokko
*
* Distributed under the Boost Software License, Version 1.0. (See accompanying
* file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*
*****************************************************************************/

#ifndef ROKKO_LAPACK_H
#define ROKKO_LAPACK_H

#include <rokko/mangling.h>
#include <lapacke.h>

#ifdef __cplusplus

#include <boost/numeric/ublas/traits.hpp>

#define LAPACKE_GEQRF_DECL(NAME, TYPE) \
lapack_int LAPACKE_geqrf(int matrix_order, lapack_int m, lapack_int n, TYPE * a, lapack_int lda, TYPE * tau);
LAPACKE_GEQRF_DECL(sgeqrf, float);
LAPACKE_GEQRF_DECL(dgeqrf, double);
LAPACKE_GEQRF_DECL(cgeqrf, lapack_complex_float);
LAPACKE_GEQRF_DECL(zgeqrf, lapack_complex_double);
#undef LAPACKE_GEQRF_DECL

#define LAPACKE_GESVD_DECL(NAME, TYPE) \
lapack_int LAPACKE_gesvd(int matrix_order, char jobu, char jobvt, lapack_int m, lapack_int n, TYPE * a, lapack_int lda, boost::numeric::ublas::type_traits<TYPE>::real_type * s, TYPE * u, lapack_int ldu, TYPE * vt, lapack_int ldvt, boost::numeric::ublas::type_traits<TYPE>::real_type * superb);
LAPACKE_GESVD_DECL(sgesvd, float);
LAPACKE_GESVD_DECL(dgesvd, double);
LAPACKE_GESVD_DECL(cgesvd, lapack_complex_float);
LAPACKE_GESVD_DECL(zgesvd, lapack_complex_double);
#undef LAPACKE_GESVD_DECL

#define LAPACKE_HEEV_DECL(NAME, TYPE) \
lapack_int LAPACKE_heev(int matrix_order, char jobz, char uplo, lapack_int n, TYPE * a, lapack_int lda, boost::numeric::ublas::type_traits<TYPE>::real_type * w);
LAPACKE_HEEV_DECL(cheev, lapack_complex_float);
LAPACKE_HEEV_DECL(zheev, lapack_complex_double);
#undef LAPACKE_HEEV_DECL

#define LAPACKE_ORGQR_DECL(NAME, TYPE) \
lapack_int LAPACKE_orgqr(int matrix_order, lapack_int m, lapack_int n, lapack_int k, TYPE * a, lapack_int lda, const TYPE * tau); \
lapack_int LAPACKE_ungqr(int matrix_order, lapack_int m, lapack_int n, lapack_int k, TYPE * a, lapack_int lda, const TYPE * tau);
LAPACKE_ORGQR_DECL(sorgqr, float);
LAPACKE_ORGQR_DECL(dorgqr, double);
#undef LAPACKE_ORGQR_DECL

#define LAPACKE_SYEV_DECL(NAME, TYPE) \
lapack_int LAPACKE_syev(int matrix_order, char jobz, char uplo, lapack_int n, TYPE * a, lapack_int lda, TYPE * w ); \
lapack_int LAPACKE_heev(int matrix_order, char jobz, char uplo, lapack_int n, TYPE * a, lapack_int lda, TYPE * w );
LAPACKE_SYEV_DECL(ssyev, float);
LAPACKE_SYEV_DECL(dsyev, double);
#undef LAPACKE_SYEV_DECL

#define LAPACKE_UNGQR_DECL(NAME, TYPE) \
lapack_int LAPACKE_ungqr(int matrix_order, lapack_int m, lapack_int n, lapack_int k, TYPE * a, lapack_int lda, const TYPE * tau);
LAPACKE_UNGQR_DECL(cungqr, lapack_complex_float);
LAPACKE_UNGQR_DECL(zungqr, lapack_complex_double);
#undef LAPACKE_UNGQR_DECL

#define LAPACKE_SYGV_DECL(NAME, TYPE) \
lapack_int LAPACKE_sygv(int matrix_order, lapack_int itype, char jobz, char uplo, lapack_int n, TYPE * a, lapack_int lda, TYPE * b, lapack_int ldb, TYPE * w ); \
lapack_int LAPACKE_heev(int matrix_order, lapack_int itype, char jobz, char uplo, lapack_int n, TYPE * a, lapack_int lda, TYPE * b, lapack_int ldb, TYPE * w );
LAPACKE_SYGV_DECL(ssygv, float);
LAPACKE_SYGV_DECL(dsygv, double);
#undef LAPACKE_SYGV_DECL

#endif // __cplusplus

#endif // ROKKO_LAPACK_H

max_stars_repo_name: wistaria/rokko

max_issues_repo_name: wistaria/rokko

max_forks_repo_name: wistaria/rokko