text:
/* randist/gsl-randist.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, 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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_test.h>

void error (const char * s);

int
main (int argc, char *argv[])
{
  size_t i,j;
  size_t n = 0;
  double mu = 0, nu = 0, nu1 = 0, nu2 = 0, sigma = 0, a = 0, b = 0, c = 0;
  double zeta = 0, sigmax = 0, sigmay = 0, rho = 0;
  double p = 0;
  double x = 0, y =0, z=0  ;
  unsigned int N = 0, t = 0, n1 = 0, n2 = 0 ;
  unsigned long int seed = 0 ;
  const char * name ;
  gsl_rng * r ;

  if (argc < 4)
    {
      printf (
"Usage: gsl-randist seed n DIST param1 param2 ...\n"
"Generates n samples from the distribution DIST with parameters param1,\n"
"param2, etc. Valid distributions are,\n\n");

      printf(
"  beta\n"
"  binomial\n"
"  bivariate-gaussian\n"
"  cauchy\n"
"  chisq\n"
"  dir-2d\n"
"  dir-3d\n"
"  dir-nd\n"
"  erlang\n"
"  exponential\n"
"  exppow\n"
"  fdist\n"
"  flat\n"
"  gamma\n"
"  gaussian-tail\n"
"  gaussian\n"
"  geometric\n"
"  gumbel1\n"
"  gumbel2\n"
"  hypergeometric\n"
"  laplace\n"
"  landau\n"
"  levy\n"
"  levy-skew\n"
"  logarithmic\n"
"  logistic\n"
"  lognormal\n"
"  negative-binomial\n"
"  pareto\n"
"  pascal\n"
"  poisson\n"
"  rayleigh-tail\n"
"  rayleigh\n"
"  tdist\n"
"  ugaussian-tail\n"
"  ugaussian\n"
"  weibull\n") ;
      exit (0);
    }

  argv++ ; seed = atol (argv[0]); argc-- ;
  argv++ ; n = atol (argv[0]); argc-- ;
  argv++ ; name = argv[0] ; argc-- ; argc-- ;

  gsl_rng_env_setup() ;

  if (gsl_rng_default_seed != 0) {
    fprintf(stderr,
            "overriding GSL_RNG_SEED with command line value, seed = %ld\n",
            seed) ;
  }

  gsl_rng_default_seed = seed ;

  r = gsl_rng_alloc(gsl_rng_default) ;

#define NAME(x) !strcmp(name,(x))
#define OUTPUT(x) for (i = 0; i < n; i++) { printf("%g\n", (x)) ; }
#define OUTPUT1(a,x) for(i = 0; i < n; i++) { a ; printf("%g\n", x) ; }
#define OUTPUT2(a,x,y) for(i = 0; i < n; i++) { a ; printf("%g %g\n", x, y) ; }
#define OUTPUT3(a,x,y,z) for(i = 0; i < n; i++) { a ; printf("%g %g %g\n", x, y, z) ; }
#define INT_OUTPUT(x) for (i = 0; i < n; i++) { printf("%d\n", (x)) ; }
#define ARGS(x,y) if (argc != x) error(y) ;
#define DBL_ARG(x) if (argc) { x=atof((++argv)[0]);argc--;} else {error( #x);};
#define INT_ARG(x) if (argc) { x=atoi((++argv)[0]);argc--;} else {error( #x);};

  if (NAME("bernoulli"))
    {
      ARGS(1, "p = probability of success");
      DBL_ARG(p)
      INT_OUTPUT(gsl_ran_bernoulli (r, p));
    }
  else if (NAME("beta"))
    {
      ARGS(2, "a,b = shape parameters");
      DBL_ARG(a)
      DBL_ARG(b)
      OUTPUT(gsl_ran_beta (r, a, b));
    }
  else if (NAME("binomial"))
    {
      ARGS(2, "p = probability, N = number of trials");
      DBL_ARG(p)
      INT_ARG(N)
      INT_OUTPUT(gsl_ran_binomial (r, p, N));
    }
  else if (NAME("cauchy"))
    {
      ARGS(1, "a = scale parameter");
      DBL_ARG(a)
      OUTPUT(gsl_ran_cauchy (r, a));
    }
  else if (NAME("chisq"))
    {
      ARGS(1, "nu = degrees of freedom");
      DBL_ARG(nu)
      OUTPUT(gsl_ran_chisq (r, nu));
    }
  else if (NAME("erlang"))
    {
      ARGS(2, "a = scale parameter, b = order");
      DBL_ARG(a)
      DBL_ARG(b)
      OUTPUT(gsl_ran_erlang (r, a, b));
    }
  else if (NAME("exponential"))
    {
      ARGS(1, "mu = mean value");
      DBL_ARG(mu) ;
      OUTPUT(gsl_ran_exponential (r, mu));
    }
  else if (NAME("exppow"))
    {
      ARGS(2, "a = scale parameter, b = power (1=exponential, 2=gaussian)");
      DBL_ARG(a) ;
      DBL_ARG(b) ;
      OUTPUT(gsl_ran_exppow (r, a, b));
    }
  else if (NAME("fdist"))
    {
      ARGS(2, "nu1, nu2 = degrees of freedom parameters");
      DBL_ARG(nu1) ;
      DBL_ARG(nu2) ;
      OUTPUT(gsl_ran_fdist (r, nu1, nu2));
    }
  else if (NAME("flat"))
    {
      ARGS(2, "a = lower limit, b = upper limit");
      DBL_ARG(a) ;
      DBL_ARG(b) ;
      OUTPUT(gsl_ran_flat (r, a, b));
    }
  else if (NAME("gamma"))
    {
      ARGS(2, "a = order, b = scale");
      DBL_ARG(a) ;
      DBL_ARG(b) ;
      OUTPUT(gsl_ran_gamma (r, a, b));
    }
  else if (NAME("gaussian"))
    {
      ARGS(1, "sigma = standard deviation");
      DBL_ARG(sigma) ;
      OUTPUT(gsl_ran_gaussian (r, sigma));
    }
  else if (NAME("gaussian-tail"))
    {
      ARGS(2, "a = lower limit, sigma = standard deviation");
      DBL_ARG(a) ;
      DBL_ARG(sigma) ;
      OUTPUT(gsl_ran_gaussian_tail (r, a, sigma));
    }
  else if (NAME("ugaussian"))
    {
      ARGS(0, "unit gaussian, no parameters required");
      OUTPUT(gsl_ran_ugaussian (r));
    }
  else if (NAME("ugaussian-tail"))
    {
      ARGS(1, "a = lower limit");
      DBL_ARG(a) ;
      OUTPUT(gsl_ran_ugaussian_tail (r, a));
    }
  else if (NAME("bivariate-gaussian"))
    {
      ARGS(3, "sigmax = x std.dev., sigmay = y std.dev., rho = correlation");
      DBL_ARG(sigmax) ;
      DBL_ARG(sigmay) ;
      DBL_ARG(rho) ;
      OUTPUT2(gsl_ran_bivariate_gaussian (r, sigmax, sigmay, rho, &x, &y),
              x, y);
    }
  else if (NAME("dir-2d"))
    {
      OUTPUT2(gsl_ran_dir_2d (r, &x, &y), x, y);
    }
  else if (NAME("dir-3d"))
    {
      OUTPUT3(gsl_ran_dir_3d (r, &x, &y, &z), x, y, z);
    }
  else if (NAME("dir-nd"))
    {
      double *xarr;
      ARGS(1, "n1 = number of dimensions of hypersphere");
      INT_ARG(n1) ;
      xarr = (double *)malloc(n1*sizeof(double));

      for(i = 0; i < n; i++) {
        gsl_ran_dir_nd (r, n1, xarr) ;
        for (j = 0; j < n1; j++) {
          if (j) putchar(' ');
          printf("%g", xarr[j]) ;
        }
        putchar('\n');
      } ;

      free(xarr);
    }
  else if (NAME("geometric"))
    {
      ARGS(1, "p = bernoulli trial probability of success");
      DBL_ARG(p) ;
      INT_OUTPUT(gsl_ran_geometric (r, p));
    }
  else if (NAME("gumbel1"))
    {
      ARGS(2, "a = order, b = scale parameter");
      DBL_ARG(a) ;
      DBL_ARG(b) ;
      OUTPUT(gsl_ran_gumbel1 (r, a, b));
    }
  else if (NAME("gumbel2"))
    {
      ARGS(2, "a = order, b = scale parameter");
      DBL_ARG(a) ;
      DBL_ARG(b) ;
      OUTPUT(gsl_ran_gumbel2 (r, a, b));
    }
  else if (NAME("hypergeometric"))
    {
      ARGS(3, "n1 = tagged population, n2 = untagged population, t = number of trials");
      INT_ARG(n1) ;
      INT_ARG(n2) ;
      INT_ARG(t) ;
      INT_OUTPUT(gsl_ran_hypergeometric (r, n1, n2, t));
    }
  else if (NAME("laplace"))
    {
      ARGS(1, "a = scale parameter");
      DBL_ARG(a) ;
      OUTPUT(gsl_ran_laplace (r, a));
    }
  else if (NAME("landau"))
    {
      ARGS(0, "no arguments required");
      OUTPUT(gsl_ran_landau (r));
    }
  else if (NAME("levy"))
    {
      ARGS(2, "c = scale, a = power (1=cauchy, 2=gaussian)");
      DBL_ARG(c) ;
      DBL_ARG(a) ;
      OUTPUT(gsl_ran_levy (r, c, a));
    }
  else if (NAME("levy-skew"))
    {
      ARGS(3, "c = scale, a = power (1=cauchy, 2=gaussian), b = skew");
      DBL_ARG(c) ;
      DBL_ARG(a) ;
      DBL_ARG(b) ;
      OUTPUT(gsl_ran_levy_skew (r, c, a, b));
    }
  else if (NAME("logarithmic"))
    {
      ARGS(1, "p = probability");
      DBL_ARG(p) ;
      INT_OUTPUT(gsl_ran_logarithmic (r, p));
    }
  else if (NAME("logistic"))
    {
      ARGS(1, "a = scale parameter");
      DBL_ARG(a) ;
      OUTPUT(gsl_ran_logistic (r, a));
    }
  else if (NAME("lognormal"))
    {
      ARGS(2, "zeta = location parameter, sigma = scale parameter");
      DBL_ARG(zeta) ;
      DBL_ARG(sigma) ;
      OUTPUT(gsl_ran_lognormal (r, zeta, sigma));
    }
  else if (NAME("negative-binomial"))
    {
      ARGS(2, "p = probability, a = order");
      DBL_ARG(p) ;
      DBL_ARG(a) ;
      INT_OUTPUT(gsl_ran_negative_binomial (r, p, a));
    }
  else if (NAME("pareto"))
    {
      ARGS(2, "a = power, b = scale parameter");
      DBL_ARG(a) ;
      DBL_ARG(b) ;
      OUTPUT(gsl_ran_pareto (r, a, b));
    }
  else if (NAME("pascal"))
    {
      ARGS(2, "p = probability, n = order (integer)");
      DBL_ARG(p) ;
      INT_ARG(N) ;
      INT_OUTPUT(gsl_ran_pascal (r, p, N));
    }
  else if (NAME("poisson"))
    {
      ARGS(1, "mu = scale parameter");
      DBL_ARG(mu) ;
      INT_OUTPUT(gsl_ran_poisson (r, mu));
    }
  else if (NAME("rayleigh"))
    {
      ARGS(1, "sigma = scale parameter");
      DBL_ARG(sigma) ;
      OUTPUT(gsl_ran_rayleigh (r, sigma));
    }
  else if (NAME("rayleigh-tail"))
    {
      ARGS(2, "a = lower limit, sigma = scale parameter");
      DBL_ARG(a) ;
      DBL_ARG(sigma) ;
      OUTPUT(gsl_ran_rayleigh_tail (r, a, sigma));
    }
  else if (NAME("tdist"))
    {
      ARGS(1, "nu = degrees of freedom");
      DBL_ARG(nu) ;
      OUTPUT(gsl_ran_tdist (r, nu));
    }
  else if (NAME("weibull"))
    {
      ARGS(2, "a = scale parameter, b = exponent");
      DBL_ARG(a) ;
      DBL_ARG(b) ;
      OUTPUT(gsl_ran_weibull (r, a, b));
    }
  else
    {
      fprintf(stderr,"Error: unrecognized distribution: %s\n", name) ;
    }

  return 0 ;
}

void
error (const char * s)
{
  fprintf(stderr, "Error: arguments should be %s\n",s) ;
  exit (EXIT_FAILURE) ;
}


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_poly.h>

/* Find a minimum in x=[0,1] of the interpolating quadratic through
 * (0,f0) (1,f1) with derivative fp0 at x=0.  The interpolating
 * polynomial is q(x) = f0 + fp0 * z + (f1-f0-fp0) * z^2
 */

static double
interp_quad (double f0, double fp0, double f1, double zl, double zh)
{
  double fl = f0 + zl*(fp0 + zl*(f1 - f0 -fp0));
  double fh = f0 + zh*(fp0 + zh*(f1 - f0 -fp0));
  double c = 2 * (f1 - f0 - fp0);       /* curvature */

  double zmin = zl, fmin = fl;

  if (fh < fmin) { zmin = zh; fmin = fh; }

  if (c > 0)  /* positive curvature required for a minimum */
    {
      double z = -fp0 / c;      /* location of minimum */
      if (z > zl && z < zh) {
        double f = f0 + z*(fp0 + z*(f1 - f0 -fp0));
        if (f < fmin) { zmin = z; fmin = f; };
      }
    }

  return zmin;
}

/* Find a minimum in x=[0,1] of the interpolating cubic through
 * (0,f0) (1,f1) with derivatives fp0 at x=0 and fp1 at x=1.
 *
 * The interpolating polynomial is:
 *
 * c(x) = f0 + fp0 * z + eta * z^2 + xi * z^3
 *
 * where eta=3*(f1-f0)-2*fp0-fp1, xi=fp0+fp1-2*(f1-f0).
 */

static double
cubic (double c0, double c1, double c2, double c3, double z)
{
  return c0 + z * (c1 + z * (c2 + z * c3));
}

static void
check_extremum (double c0, double c1, double c2, double c3, double z,
                double *zmin, double *fmin)
{
  /* could make an early return by testing curvature >0 for minimum */

  double y = cubic (c0, c1, c2, c3, z);

  if (y < *fmin)
    {
      *zmin = z;  /* accepted new point*/
      *fmin = y;
    }
}

static double
interp_cubic (double f0, double fp0, double f1, double fp1, double zl, double zh)
{
  double eta = 3 * (f1 - f0) - 2 * fp0 - fp1;
  double xi = fp0 + fp1 - 2 * (f1 - f0);
  double c0 = f0, c1 = fp0, c2 = eta, c3 = xi;
  double zmin, fmin;
  double z0, z1;

  zmin = zl; fmin = cubic(c0, c1, c2, c3, zl);
  check_extremum (c0, c1, c2, c3, zh, &zmin, &fmin);

  {
    int n = gsl_poly_solve_quadratic (3 * c3, 2 * c2, c1, &z0, &z1);

    if (n == 2)  /* found 2 roots */
      {
        if (z0 > zl && z0 < zh)
          check_extremum (c0, c1, c2, c3, z0, &zmin, &fmin);
        if (z1 > zl && z1 < zh)
          check_extremum (c0, c1, c2, c3, z1, &zmin, &fmin);
      }
    else if (n == 1)  /* found 1 root */
      {
        if (z0 > zl && z0 < zh)
          check_extremum (c0, c1, c2, c3, z0, &zmin, &fmin);
      }
  }

  return zmin;
}

static double
interpolate (double a, double fa, double fpa,
             double b, double fb, double fpb, double xmin, double xmax,
             int order)
{
  /* Map [a,b] to [0,1] */
  double z, alpha, zmin, zmax;

  zmin = (xmin - a) / (b - a);
  zmax = (xmax - a) / (b - a);

  if (zmin > zmax)
    {
      double tmp = zmin;
      zmin = zmax;
      zmax = tmp;
    };

  if (order > 2 && GSL_IS_REAL(fpb)) {
    z = interp_cubic (fa, fpa * (b - a), fb, fpb * (b - a), zmin, zmax);
  } else {
    z = interp_quad (fa, fpa * (b - a), fb, zmin, zmax);
  }

  alpha = a + z * (b - a);

  return alpha;
}

/* recommended values from Fletcher are
   rho = 0.01, sigma = 0.1, tau1 = 9, tau2 = 0.05, tau3 = 0.5 */

static int
minimize (gsl_function_fdf * fn, double rho, double sigma,
          double tau1, double tau2, double tau3,
          int order, double alpha1, double *alpha_new)
{
  double f0, fp0, falpha, falpha_prev, fpalpha, fpalpha_prev, delta,
    alpha_next;
  double alpha = alpha1, alpha_prev = 0.0;
  double a, b, fa, fb, fpa, fpb;
  const size_t bracket_iters = 100, section_iters = 100;
  size_t i = 0;

  GSL_FN_FDF_EVAL_F_DF (fn, 0.0, &f0, &fp0);
  falpha_prev = f0;
  fpalpha_prev = fp0;

  /* Avoid uninitialized variables morning */
  a = 0.0; b = alpha;
  fa = f0; fb = 0.0;
  fpa = fp0; fpb = 0.0;

  /* Begin bracketing */

  while (i++ < bracket_iters)
    {
      falpha = GSL_FN_FDF_EVAL_F (fn, alpha);

      /* Fletcher's rho test */

      if (falpha > f0 + alpha * rho * fp0 || falpha >= falpha_prev)
        {
          a = alpha_prev; fa = falpha_prev; fpa = fpalpha_prev;
          b = alpha; fb = falpha; fpb = GSL_NAN;
          break;                /* goto sectioning */
        }

      fpalpha = GSL_FN_FDF_EVAL_DF (fn, alpha);

      /* Fletcher's sigma test */

      if (fabs (fpalpha) <= -sigma * fp0)
        {
          *alpha_new = alpha;
          return GSL_SUCCESS;
        }

      if (fpalpha >= 0)
        {
          a = alpha; fa = falpha; fpa = fpalpha;
          b = alpha_prev; fb = falpha_prev; fpb = fpalpha_prev;
          break;                /* goto sectioning */
        }

      delta = alpha - alpha_prev;

      {
        double lower = alpha + delta;
        double upper = alpha + tau1 * delta;

        alpha_next = interpolate (alpha_prev, falpha_prev, fpalpha_prev,
                             alpha, falpha, fpalpha, lower, upper, order);

      }

      alpha_prev = alpha;
      falpha_prev = falpha;
      fpalpha_prev = fpalpha;
      alpha = alpha_next;
    }

  /*  Sectioning of bracket [a,b] */

  while (i++ < section_iters)
    {
      delta = b - a;

      {
        double lower = a + tau2 * delta;
        double upper = b - tau3 * delta;

        alpha = interpolate (a, fa, fpa, b, fb, fpb, lower, upper, order);
      }

      falpha = GSL_FN_FDF_EVAL_F (fn, alpha);

      if ((a-alpha)*fpa <= GSL_DBL_EPSILON) {
        /* roundoff prevents progress */
        return GSL_ENOPROG;
      };

      if (falpha > f0 + rho * alpha * fp0 || falpha >= fa)
        {
          /*  a_next = a; */
          b = alpha; fb = falpha; fpb = GSL_NAN;
        }
      else
        {
          fpalpha = GSL_FN_FDF_EVAL_DF (fn, alpha);

          if (fabs(fpalpha) <= -sigma * fp0)
            {
              *alpha_new = alpha;
              return GSL_SUCCESS;  /* terminate */
            }

          if ( ((b-a) >= 0 && fpalpha >= 0) || ((b-a) <=0 && fpalpha <= 0))
            {
              b = a; fb = fa; fpb = fpa;
              a = alpha; fa = falpha; fpa = fpalpha;
            }
          else
            {
              a = alpha; fa = falpha; fpa = fpalpha;
            }
        }
    }

  return GSL_SUCCESS;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <gsl/gsl>
#include <halley/text/halleystring.h>
#include "lua_reference.h"

#include "halley/data_structures/hash_map.h"

struct lua_State;

namespace Halley {
	class Resources;
	class LuaState;

	class LuaState {
	public:
		LuaState(Resources& resources);
		~LuaState();

		const LuaReference* tryGetModule(const String& moduleName) const;
		const LuaReference& getModule(const String& moduleName) const;
		const LuaReference& getOrLoadModule(const String& moduleName);
		const LuaReference& loadModule(const String& moduleName, gsl::span<const gsl::byte> data);
		void unloadModule(const String& moduleName);

		void call(int nArgs, int nRets);

		lua_State* getRawState();

		void pushCallback(LuaCallback&& callback);

		void pushErrorHandler();
		void popErrorHandler();
		void pushLuaState(lua_State* lua);
		void popLuaState();
		String errorHandler(String message);

	private:
		lua_State* lua;
		Vector<lua_State*> pushedStates;
		Resources* resources;

		HashMap<String, LuaReference> modules;
		Vector<std::unique_ptr<LuaCallback>> closures;
		std::unique_ptr<LuaReference> errorHandlerRef;
		Vector<int> errorHandlerStackPos;

		LuaReference loadScript(const String& chunkName, gsl::span<const gsl::byte> data);

		void print(String string);
		const LuaReference& packageLoader(String moduleName);
		String printVariableAtTop(int maxDepth = 2, bool quote = true);
	};

	class LuaStateOverrider {
	public:
		LuaStateOverrider(LuaState& state, lua_State* rawState);
		~LuaStateOverrider();

	private:
		LuaState& state;
	};
}

ext: h

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


text:
#ifndef _SPATIAL_CONVOLUTION_
#define _SPATIAL_CONVOLUTION_

#include "Mem_Alloc.h"
#include "Time.h"
#include "Global_Setting.h"
//#include <cblas.h>

extern int spatial_convolution_float32(float* input_ptr, float* kernel_ptr, float* output_ptr);
extern int spatial_convolution_float32_mm(float* input_ptr, float* kernel_ptr, float* output_ptr);
extern int spatial_convolution_float32_mm_blas(float* input_ptr, float* kernel_ptr, float* output_ptr);

#endif	//	_SPATIAL_CONVOLUTION_


text:
/* Function in common for mlfit1, mlfitN, etc ...
 */

#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <stdint.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>

#include "mlfit.h"
#include "blit3.h"

#ifndef verbose
#define verbose 0
#endif

double estimateNphot(double *V, size_t Vm)
  // V is a Vm x Vm matrix
  // The number of photons locally as sum(V-min(V))
{
  double s = 0;
  double minV = INFINITY;

  for(size_t kk=0; kk<Vm*Vm; kk++)
    s+=V[kk];

  for(size_t kk=0; kk<Vm*Vm; kk++)
    minV = GSL_MIN(V[kk], minV);

  s = s-minV*Vm*Vm;

  //printf("estimateNphot: %f\n", s);
  return s;
}

int double_cmp(const void * a, const void * b)
{
  if (*(const double*)a > *(const double*)b) {
    return 1;
  } else {
    if (*(const double*)a < *(const double*)b) {
      return -1;
    } else {
    return 0;
  }
}
}

double estimateBGV(double *V, size_t Vm, size_t Vn, size_t Vp, double * D)
  // estimate the background level in V [VmxVnxVp] around the point D
  // (x,y,z)
{

  int32_t radius = 7;
  double * ePixels = malloc(sizeof(double)*6*(2*radius+1)*(2*radius+1)); // a few more than we need

  int32_t x = nearbyint(D[0]);
  int32_t y = nearbyint(D[1]);
  int32_t z = nearbyint(D[2]);

  // Loop over the box, and copy the edge pixels to ePixels
  size_t nPixels = 0;
  for(int32_t zz = z-radius; zz<=z+radius; zz++)
    for(int32_t yy = y-radius; yy<=y+radius; yy++)
     for(int32_t xx = x-radius; xx<=x+radius; xx++)
       if((abs(zz-z)+abs(yy-y)+abs(xx-x)) == radius)
         if(xx>=0 && yy>=0 && zz>=0 && xx< (int32_t) Vm && yy<(int32_t) Vn && zz<(int32_t) Vp)
           ePixels[nPixels++] = V[xx + yy*Vm + zz*Vm*Vn];

#if verbose > 0
  printf("%lu ePixels\n", nPixels);
#endif

  double median = 0;
  if(nPixels>0)
  {
  qsort(ePixels, nPixels, sizeof(double),  double_cmp);
  gsl_sort(ePixels, 1, nPixels);
  median = gsl_stats_median_from_sorted_data(ePixels, 1, nPixels);
  //median = (double) nPixels;
#if verbose > 0
  printf("Median: %f\n", median);
#endif
  }
  else
  { median = 0;}

  free(ePixels);
  return(median);
  //return((double) q);
}

double estimateBG(double *V, size_t Vm)
  // Estimate the background level as the median of the edge pixels
{

  size_t nEdge = 4*Vm-4;
  double * ePixels = malloc(nEdge*sizeof(double));

  size_t pos = 0;
  for(uint32_t kk = 0; kk<Vm; kk++, pos++) // Left
    ePixels[pos] = V[kk];
  for(uint32_t kk = 0; kk<Vm; kk++, pos++) // Right
    ePixels[pos] = V[kk+(Vm-1)*Vm];
  for(uint32_t kk = 1; kk<Vm-1; kk++, pos++) // Top
    ePixels[pos] = V[Vm*kk];
  for(uint32_t kk = 1; kk<Vm-1; kk++, pos++) // Bottom
    ePixels[pos] = (double) V[Vm-1+Vm*kk];

  assert(pos == 4*Vm-4);

  if(0) {
    for(uint32_t kk = 0; kk<nEdge; kk++)
      printf("%2d :%f\n", kk, ePixels[kk]);
    printf("\n");
  }

//  qsort(ePixels, nEdge, sizeof(double),  double_cmp);

  if(0) {
    for(uint32_t kk = 0; kk<nEdge; kk++)
      printf("%2d %f\n", kk, ePixels[kk]);
    printf("\n");
  }

  gsl_sort(ePixels, 1, nEdge);
  double median = gsl_stats_median_from_sorted_data(ePixels, 1, nEdge);
  // double mean = gsl_stats_mean(ePixels, 1, nEdge);
#if verbose > 0
  printf("Median: %f\n", median);
#endif
  free(ePixels);
  return(median);
}

int getZLine(double *W, size_t Ws,
    double * V, size_t Vm, size_t Vn, size_t Vp,
    double *D)
  /* Copy a line from V into W.
   * The line will be D[x, y, z-hWs:z+hWs]
   * where hWs = (Ws-1)/2 and D = [x,y,z]
   *
   * returns 1 on failure (i.e. line is out of bounds)
   * returns 0 if ok.
   *
   * See also: getRegion
   */
{

  int64_t x = nearbyint(D[0]);
  int64_t y = nearbyint(D[1]);
  int64_t z = nearbyint(D[2]);

  size_t Ws2 = (Ws-1)/2;

  if(z+Ws2 >= Vp)
    return 1;
  if(z<(int64_t) Ws2)
    return 1;

  size_t pos = 0;
  for(size_t zz = z-Ws2; zz<=z+Ws2; zz++)
  {
    W[pos++] = V[x + y*Vm + zz*Vm*Vn];
  }
  return 0;
}

int getRegion(double * W, size_t Ws,
    double * V, size_t Vm, size_t Vn, size_t Vp,
    double * D)
  // Get a 2D region for constant z or returns 0 if out of bounds
  // Vm, Vn, Vp is the centre of the region.
{
  size_t Ws2 = (Ws-1)/2;
  size_t pos = 0;

  int64_t x = nearbyint(D[0]);
  int64_t y = nearbyint(D[1]);
  int64_t z = nearbyint(D[2]);

#if verbose > 0
  printf("getRegion round(D): (%lu %lu %lu)\n", x, y, z);
  printf("getRegion size(V)   (%lu %lu %lu)\n", Vm, Vn, Vp);
#endif

  if(x-(int64_t) Ws2 < 0)
    return 1;
  if(y-(int64_t) Ws2 < 0)
    return 1;
  if(x+(uint64_t) Ws2  >= Vm)
    return 1;
  if(y+(uint64_t) Ws2  >= Vn)
    return 1;
  if(z<0)
    return 1;
  if((size_t) z>=Vp)
    return 1;

#if verbose > 1
  printf("Region valid\n");
  printf("%lu %lu %lu\n", Vm, Vn, Vp);
#endif

  uint32_t zz = z;
  for(uint32_t yy = y-Ws2; yy<=y+Ws2; yy++) {
    for(uint32_t xx = x-Ws2; xx<=x+Ws2; xx++) {
      assert(xx<Vm); assert(yy<Vn); assert(zz<Vp);
    // printf("x %lu y %lu z %lu\n", x, y, z);

      W[pos++] = V[xx + yy*Vm + zz*Vm*Vn];
    }
  }
  //  showRegion(W, Ws);
 return 0;
}

void showRegion(double * W, size_t Ws)
  // Dump some quadratic region to the terminal
{
printf("Show region ...\n");
  for(size_t xx = 0 ; xx<Ws; xx++) {
    for(size_t yy = 0; yy<Ws; yy++) {
      printf("%f ", W[yy*Ws+xx]);
    }
    printf("\n");
  }
}

max_stars_repo_name: elgw/dotter

max_issues_repo_name: elgw/dotter

max_forks_repo_name: elgw/dotter


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

#include "tests.h"

void
test_axpy (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int N = 1;
   float alpha = 0.0f;
   float X[] = { 0.018f };
   int incX = 1;
   float Y[] = { -0.417f };
   int incY = -1;
   float expected[] = { -0.417f };
   cblas_saxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected[i], flteps, "saxpy(case 64)");
     }
   };
  };

  {
   int N = 1;
   double alpha = 0;
   double X[] = { 0.071 };
   int incX = 1;
   double Y[] = { -0.888 };
   int incY = -1;
   double expected[] = { -0.888 };
   cblas_daxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected[i], dbleps, "daxpy(case 65)");
     }
   };
  };

  {
   int N = 1;
   float alpha[2] = {1.0f, 0.0f};
   float X[] = { -0.542f, -0.362f };
   int incX = 1;
   float Y[] = { -0.459f, -0.433f };
   int incY = -1;
   float expected[] = { -1.001f, -0.795f };
   cblas_caxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected[2*i], flteps, "caxpy(case 66) real");
       gsl_test_rel(Y[2*i+1], expected[2*i+1], flteps, "caxpy(case 66) imag");
     };
   };
  };

  {
   int N = 1;
   double alpha[2] = {-1, 0};
   double X[] = { 0.003, -0.514 };
   int incX = 1;
   double Y[] = { -0.529, 0.743 };
   int incY = -1;
   double expected[] = { -0.532, 1.257 };
   cblas_zaxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected[2*i], dbleps, "zaxpy(case 67) real");
       gsl_test_rel(Y[2*i+1], expected[2*i+1], dbleps, "zaxpy(case 67) imag");
     };
   };
  };

  {
   int N = 1;
   float alpha = 0.1f;
   float X[] = { 0.771f };
   int incX = -1;
   float Y[] = { 0.507f };
   int incY = 1;
   float expected[] = { 0.5841f };
   cblas_saxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected[i], flteps, "saxpy(case 68)");
     }
   };
  };

  {
   int N = 1;
   double alpha = -0.3;
   double X[] = { 0.029 };
   int incX = -1;
   double Y[] = { -0.992 };
   int incY = 1;
   double expected[] = { -1.0007 };
   cblas_daxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected[i], dbleps, "daxpy(case 69)");
     }
   };
  };

  {
   int N = 1;
   float alpha[2] = {-0.3f, 0.1f};
   float X[] = { 0.194f, -0.959f };
   int incX = -1;
   float Y[] = { 0.096f, 0.032f };
   int incY = 1;
   float expected[] = { 0.1337f, 0.3391f };
   cblas_caxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected[2*i], flteps, "caxpy(case 70) real");
       gsl_test_rel(Y[2*i+1], expected[2*i+1], flteps, "caxpy(case 70) imag");
     };
   };
  };

  {
   int N = 1;
   double alpha[2] = {0, 1};
   double X[] = { 0.776, -0.671 };
   int incX = -1;
   double Y[] = { 0.39, 0.404 };
   int incY = 1;
   double expected[] = { 1.061, 1.18 };
   cblas_zaxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected[2*i], dbleps, "zaxpy(case 71) real");
       gsl_test_rel(Y[2*i+1], expected[2*i+1], dbleps, "zaxpy(case 71) imag");
     };
   };
  };

  {
   int N = 1;
   float alpha = 1.0f;
   float X[] = { 0.647f };
   int incX = -1;
   float Y[] = { 0.016f };
   int incY = -1;
   float expected[] = { 0.663f };
   cblas_saxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected[i], flteps, "saxpy(case 72)");
     }
   };
  };

  {
   int N = 1;
   double alpha = -1;
   double X[] = { -0.558 };
   int incX = -1;
   double Y[] = { 0.308 };
   int incY = -1;
   double expected[] = { 0.866 };
   cblas_daxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected[i], dbleps, "daxpy(case 73)");
     }
   };
  };

  {
   int N = 1;
   float alpha[2] = {-0.3f, 0.1f};
   float X[] = { 0.899f, -0.624f };
   int incX = -1;
   float Y[] = { 0.155f, -0.33f };
   int incY = -1;
   float expected[] = { -0.0523f, -0.0529f };
   cblas_caxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected[2*i], flteps, "caxpy(case 74) real");
       gsl_test_rel(Y[2*i+1], expected[2*i+1], flteps, "caxpy(case 74) imag");
     };
   };
  };

  {
   int N = 1;
   double alpha[2] = {0, 1};
   double X[] = { -0.451, 0.768 };
   int incX = -1;
   double Y[] = { 0.007, 0.732 };
   int incY = -1;
   double expected[] = { -0.761, 0.281 };
   cblas_zaxpy(N, alpha, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected[2*i], dbleps, "zaxpy(case 75) real");
       gsl_test_rel(Y[2*i+1], expected[2*i+1], dbleps, "zaxpy(case 75) imag");
     };
   };
  };

}


text:
/**
 * @file batchg_zdotu.c
 *
 * Part of API test for Batched BLAS routines.
 *
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date
 *
 * @precisions normal z -> c d s
 *
 **/

#include <cblas.h>
#include "bblas.h"

#define COMPLEX

void batchg_zdotu_sub(
		  const int *n,
		  BBLAS_Complex64_t const * const * x,
		  const int *incx,
		  BBLAS_Complex64_t const * const * y,
		  const int *incy,
		  BBLAS_Complex64_t *dotu,
		  const int group_count, const int *group_size,
		  int* info)
{
	/* Local variables */
	char func_name[15] = "batchg_zdotu";
	int group_iter = 0;
	int offset = 0; // How many subproblems solved so far
	int end_of_group = 0; // End of the current group

	/* Check group_count */
	if (group_count < 0)
	  {
	    xerbla_batch(func_name, BBLAS_ERR_GROUP_COUNT, -1);
	    return;
	  }

	/* Check group_size and call fixed batch computation */
	for (group_iter = 0; group_iter < group_count; group_iter++)
	{
		if (group_size[group_iter] < 0)
		{
			end_of_group = end_of_group + group_size[group_iter];
			xerbla_batch(func_name, BBLAS_ERR_GROUP_SIZE, group_iter);
			info[group_iter] = BBLAS_ERR_GROUP_SIZE;
			continue;
		}

		/* Call fixed batch computation on the group */
		batchf_zdotu_sub(
				 n[group_iter],
				 x+offset,
				 incx[group_iter],
				 y+offset,
				 incy[group_iter],
				 dotu+offset,
				 group_size[group_iter],
				 info+group_iter);
		/* Update offset */
		offset += group_size[group_iter];
	} // End of group loop
}
#undef COMPLEX

max_stars_repo_name: mawussi/BBLAS-group


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

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

// http://stackoverflow.com/questions/5167269
#ifdef __MACH__
#include <mach/clock.h>
#include <mach/mach.h>
#endif

// http://stackoverflow.com/questions/7034930
double drand() {
  return (rand()+1.0)/(RAND_MAX+1.0);
}

double random_normal()  {
  return sqrt(-2*log(drand())) * cos(2*M_PI*drand());
}

double* random_array(int size) {
	double *a = malloc(sizeof(double) * size);
	int i;

	for (i = 0 ; i < size; i++) {
		a[i] = random_normal();
	}

	return a;
}

long currentTimeNanos() {
struct timespec ts;

#ifdef __MACH__
clock_serv_t cclock;
mach_timespec_t mts;
host_get_clock_service(mach_host_self(), CALENDAR_CLOCK, &cclock);
clock_get_time(cclock, &mts);
mach_port_deallocate(mach_task_self(), cclock);
ts.tv_sec = mts.tv_sec;
ts.tv_nsec = mts.tv_nsec;

#else
clock_gettime(CLOCK_REALTIME, &ts);
#endif

return (ts.tv_sec * 1000000000) + ts.tv_nsec;
}


text:
/*	$Id$ */
/*
 * Copyright (c) 2015 Kristaps Dzonsons <kristaps@kcons.eu>
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#include <assert.h>
#include <math.h>
#include <stdint.h>
#include <stdlib.h>

#include <gtk/gtk.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_histogram.h>
#include <kplot.h>

#include "extern.h"

struct simbuf *
simbuf_alloc_warm(struct kdata *hot, size_t bufsz)
{
	struct simbuf	*buf;

	g_assert(NULL != hot);
	buf = g_malloc0(sizeof(struct simbuf));
	buf->cold = kdata_buffer_alloc(bufsz);
	g_assert(NULL != buf->cold);
	return(buf);
}

struct simbuf *
simbuf_alloc(struct kdata *hot, size_t bufsz)
{
	struct simbuf	*buf;

	g_assert(NULL != hot);
	buf = g_malloc0(sizeof(struct simbuf));
	buf->hot = hot;
	buf->hotlsb = kdata_buffer_alloc(bufsz);
	g_assert(NULL != buf->hotlsb);
	buf->warm = kdata_buffer_alloc(bufsz);
	g_assert(NULL != buf->warm);
	buf->cold = kdata_buffer_alloc(bufsz);
	g_assert(NULL != buf->cold);
	return(buf);
}

void
simbuf_free(struct simbuf *buf)
{

	kdata_destroy(buf->hot);
	kdata_destroy(buf->hotlsb);
	kdata_destroy(buf->warm);
	kdata_destroy(buf->cold);
	free(buf);
}

void
simbuf_copy_hotlsb(struct simbuf *buf)
{
	int	rc;

	rc = kdata_buffer_copy(buf->hotlsb, buf->hot);
	g_assert(0 != rc);
}

void
simbuf_copy_warm(struct simbuf *buf)
{
	int	rc;

	rc = kdata_buffer_copy(buf->warm, buf->hotlsb);
	g_assert(0 != rc);
}

void
simbuf_copy_cold(struct simbuf *buf)
{
	int	rc;

	rc = kdata_buffer_copy(buf->cold, buf->warm);
	g_assert(0 != rc);
}

max_stars_repo_name: kristapsdz/bmigrate

max_issues_repo_name: kristapsdz/bmigrate

max_forks_repo_name: kristapsdz/bmigrate


text:
#pragma once

#include <type_traits>

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

namespace thrustshift {

namespace kernel {

template <typename SrcT, typename MapT, typename DstT>
__global__ void scatter(gsl_lite::span<const SrcT> src,
                        gsl_lite::span<const MapT> map,
                        gsl_lite::span<DstT> dst) {

	const auto gtid = threadIdx.x + blockIdx.x * blockDim.x;
	if (gtid < src.size()) {
		dst[map[gtid]] = src[gtid];
	}
}

} // namespace kernel

namespace async {

template <class SrcRange, class MapRange, class DstRange>
void scatter(cuda::stream_t& stream,
             SrcRange&& src,
             MapRange&& map,
             DstRange&& dst) {

	gsl_Expects(src.size() == dst.size());
	gsl_Expects(src.size() == map.size());
	gsl_Expects(src.data() != dst.data());

	if (src.empty()) {
		return;
	}

	using src_value_type =
	    typename std::remove_reference<SrcRange>::type::value_type;
	using map_index_type =
	    typename std::remove_reference<MapRange>::type::value_type;
	using dst_value_type =
	    typename std::remove_reference<DstRange>::type::value_type;

	constexpr cuda::grid::block_dimension_t block_dim = 128;
	const cuda::grid::dimension_t grid_dim =
	    (src.size() + block_dim - 1) / block_dim;
	auto c = cuda::make_launch_config(grid_dim, block_dim);
	auto k = kernel::scatter<src_value_type, map_index_type, dst_value_type>;
	cuda::enqueue_launch(k, stream, c, src, map, dst);
}

} // namespace async

} // namespace thrustshift

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
#pragma once

#include "CesiumGltf/IExtensionJsonHandler.h"
#include "CesiumGltf/Model.h"
#include "CesiumGltf/ReaderLibrary.h"
#include <functional>
#include <gsl/span>
#include <map>
#include <memory>
#include <optional>
#include <string>
#include <unordered_map>
#include <vector>

namespace CesiumGltf {

struct ReaderContext;

/**
 * @brief The result of reading a glTF model with
 * {@link GltfReader::readModel}.
 */
struct CESIUMGLTFREADER_API ModelReaderResult {
  /**
   * @brief The read model, or std::nullopt if the model could not be read.
   */
  std::optional<Model> model;

  /**
   * @brief Errors, if any, that occurred during the load process.
   */
  std::vector<std::string> errors;

  /**
   * @brief Warnings, if any, that occurred during the load process.
   */
  std::vector<std::string> warnings;
};

/**
 * @brief The result of reading an image with
 * {@link GltfReader::readImage}.
 */
struct CESIUMGLTFREADER_API ImageReaderResult {

  /**
   * @brief The {@link ImageCesium} that was read.
   *
   * This will be `std::nullopt` if the image could not be read.
   */
  std::optional<ImageCesium> image;

  /**
   * @brief Error messages that occurred while trying to read the image.
   */
  std::vector<std::string> errors;

  /**
   * @brief Warning messages that occurred while reading the image.
   */
  std::vector<std::string> warnings;
};

/**
 * @brief The state of a glTF extension.
 */
enum class ExtensionState {
  /**
   * @brief The extension is enabled.
   *
   * If a statically-typed class is available for the extension, it will be
   * used. Otherwise the extension will be represented as a
   * {@link CesiumUtility::JsonValue}.
   */
  Enabled,

  /**
   * @brief The extension is enabled but will always be deserialized as a
   * {@link CesiumUtility::JsonValue}.
   *
   * Even if a statically-typed class is available for the extension, it will
   * not be used.
   */
  JsonOnly,

  /**
   * @brief The extension is disabled.
   *
   * It will not be represented in the loaded model at all.
   */
  Disabled
};

/**
 * @brief Options for how to read a glTF.
 */
struct CESIUMGLTFREADER_API ReadModelOptions {
  /**
   * @brief Whether data URLs in buffers and images should be automatically
   * decoded as part of the load process.
   */
  bool decodeDataUrls = true;

  /**
   * @brief Whether data URLs should be cleared after they are successfully
   * decoded.
   *
   * This reduces the memory usage of the model.
   */
  bool clearDecodedDataUrls = true;

  /**
   * @brief Whether embedded images in {@link Model::buffers} should be
   * automatically decoded as part of the load process.
   *
   * The {@link ImageSpec::mimeType} property is ignored, and instead the
   * [stb_image](https://github.com/nothings/stb) library is used to decode
   * images in `JPG`, `PNG`, `TGA`, `BMP`, `PSD`, `GIF`, `HDR`, or `PIC` format.
   */
  bool decodeEmbeddedImages = true;

  /**
   * @brief Whether geometry compressed using the `KHR_draco_mesh_compression`
   * extension should be automatically decoded as part of the load process.
   */
  bool decodeDraco = true;
};

/**
 * @brief Reads glTF models and images.
 */
class CESIUMGLTFREADER_API GltfReader {
public:
  /**
   * @brief Constructs a new instance.
   */
  GltfReader();

  /**
   * @brief Registers an extension for a glTF object.
   *
   * @tparam TExtended The glTF object to extend.
   * @tparam TExtensionHandler The extension's
   * {@link CesiumJsonReader::JsonHandler}.
   * @param extensionName The name of the extension.
   */
  template <typename TExtended, typename TExtensionHandler>
  void registerExtension(const std::string& extensionName) {
    auto it =
        this->_extensions.emplace(extensionName, ObjectTypeToHandler()).first;
    it->second.insert_or_assign(
        TExtended::TypeName,
        ExtensionReaderFactory([](const ReaderContext& context) {
          return std::make_unique<TExtensionHandler>(context);
        }));
  }

  /**
   * @brief Registers an extension for a glTF object.
   *
   * The extension name is obtained from `TExtensionHandler::ExtensionName`.
   *
   * @tparam TExtended The glTF object to extend.
   * @tparam TExtensionHandler The extension's
   * {@link CesiumJsonReader::JsonHandler}.
   */
  template <typename TExtended, typename TExtensionHandler>
  void registerExtension() {
    auto it =
        this->_extensions
            .emplace(TExtensionHandler::ExtensionName, ObjectTypeToHandler())
            .first;
    it->second.insert_or_assign(
        TExtended::TypeName,
        ExtensionHandlerFactory([](const ReaderContext& context) {
          return std::make_unique<TExtensionHandler>(context);
        }));
  }

  /**
   * @brief Enables or disables a glTF extension.
   *
   * By default, all extensions are enabled. When an enabled extension is
   * encountered in the source glTF, it is read into a statically-typed
   * extension class, if one is registered, or into a
   * {@link CesiumUtility::JsonValue} if not.
   *
   * When a disabled extension is encountered in the source glTF, it is ignored
   * completely.
   *
   * An extension may also be set to `ExtensionState::JsonOnly`, in which case
   * it will be read into a {@link CesiumUtility::JsonValue} even if a
   * statically-typed extension class is registered.
   *
   * @param extensionName The name of the extension to be enabled or disabled.
   * @param newState The new state for the extension.
   */
  void
  setExtensionState(const std::string& extensionName, ExtensionState newState);

  /**
   * @brief Reads a glTF or binary glTF (GLB) from a buffer.
   *
   * @param data The buffer from which to read the glTF.
   * @param options Options for how to read the glTF.
   * @return The result of reading the glTF.
   */
  ModelReaderResult readModel(
      const gsl::span<const std::byte>& data,
      const ReadModelOptions& options = ReadModelOptions()) const;

  /**
   * @brief Reads an image from a buffer.
   *
   * The [stb_image](https://github.com/nothings/stb) library is used to decode
   * images in `JPG`, `PNG`, `TGA`, `BMP`, `PSD`, `GIF`, `HDR`, or `PIC` format.
   *
   * @param data The buffer from which to read the image.
   * @return The result of reading the image.
   */
  ImageReaderResult readImage(const gsl::span<const std::byte>& data) const;

  std::unique_ptr<IExtensionJsonHandler> createExtensionHandler(
      const ReaderContext& context,
      const std::string_view& extensionName,
      const std::string& extendedObjectType) const;

private:
  using ExtensionHandlerFactory =
      std::function<std::unique_ptr<IExtensionJsonHandler>(
          const ReaderContext&)>;
  using ObjectTypeToHandler = std::map<std::string, ExtensionHandlerFactory>;
  using ExtensionNameMap = std::map<std::string, ObjectTypeToHandler>;

  ExtensionNameMap _extensions;
  std::unordered_map<std::string, ExtensionState> _extensionStates;
};

} // namespace CesiumGltf

max_stars_repo_name: zrkcode/cesium-native

max_issues_repo_name: zrkcode/cesium-native

max_forks_repo_name: zrkcode/cesium-native


text:
/*
 *
 * ML fitting of dots using a 2D Gaussian model provided by gaussianInt2.
 * I.e., only x and y are fitted.
 *
 * Compilation:
 * For debugging:
 *  gcc  3.c -Wall -std=c99 -I/usr/local/include -L/usr/local/lib -lgsl -lgslcblas -lm -g -o 3
 *  valgrind ./3
 * For speed:
 *  gcc  -Wall -std=c99 3.c -I/usr/local/include -L/usr/local/lib -lgsl -lgslcblas -lm -O3 -o 3
 *  ./3
 * Memory accesses can also be checked with Valgrind or the gcc -fmudflap memory protection option.
 *
 * MATLAB interface in df_mlfit1.c
 * 2017.04.12. Valgrind ok.
 *
 * TODO:
 *  - What are good step sizes?
 *  - Set the initial positions with higher than integer precision if possible.
 *
 */

/*
   From the GSL documentation:

   If necessary you can turn off range checking completely without
   modifying any source files by recompiling your program with the
   preprocessor definition GSL_RANGE_CHECK_OFF. Provided your compiler
   supports inline functions the effect of turning off range checking is
   to replace calls to gsl_vector_get(v,i) by v->data[i*v->stride] and
   calls to gsl_vector_set(v,i,x) by v->data[i*v->stride]=x. Thus there
   should be no performance penalty for using the range checking
   functions when range checking is turned off.

*/

#include <assert.h>
#include <inttypes.h>
#include <stdint.h>
#include <stdio.h>

#include <gsl/gsl_multimin.h>
#include <gsl/gsl_statistics.h>

#include "gaussianInt2.h"
#include "mlfit.h"
#include "blit3.h"
#include "mlfit1.h"

// use the compile flag -D verbose=0 etc
// 0 = no information
// 1 = per dot information
// 2 = per iteration information
#ifndef verbose
#define verbose 1
#endif

// Globals
const size_t maxIterations = 10000; // for z-fitting
const double convCriteria = 1e-7; // for z-fitting

// Possibly, add the window size, denoted Ws in most places

// Headers
double lxy (const gsl_vector *, void *);
double lz (const gsl_vector *, void *);
int localizeDotXY(double *, size_t, double *,  double *,
    double); // sigma
int localizeDotZ(double *, size_t, double *,  double *,
    double); // sigma
int localize(double *, size_t, size_t, size_t, double *, size_t, double *,
    double, double); // sigma_xy, sigma_z -- size of dot
int unit_tests(void);

void lxy_df (const gsl_vector *v, void *params, gsl_vector * df);
void lxy_fdf (const gsl_vector *v, void *params, double * f, gsl_vector *df);

// Optimization constants
typedef struct {
  double * R; // local region
  double * G; // A already allocated temporary space for the Gaussian kernel
  size_t Rw; // R and G has size Rw*Rw
  double sigma; // For constant sigma fitting
  double bg;
} optParams;

double lz (const gsl_vector *v, void *params)
{
  double x;

  optParams *p = (optParams *) params;

  size_t Rw = p->Rw;
  double * R = p->R;
  double * GI = p->G;
  double sigma = p->sigma;
  double bg = p->bg;

#if verbose > 2
  printf("Rw: %lu\n", Rw);
  printf("Nphot: %lu\n", Nphot);
  printf("sigma: %f\n", sigma);
  printf("bg   : %f\n", bg);
#endif

  // Get the other parameters ...
  x = gsl_vector_get(v, 0);
  double Nphot = gsl_vector_get(v,1);

  //  printf("bg: %f, Nphot: %f, sigma: %f z: %f\n", bg, Nphot, sigma, x);

  // Create Gaussian ...
  double mu[] = {x};
  gaussianInt1(GI, mu, &sigma, Rw);
  for(size_t kk = 0; kk<Rw; kk++)
    GI[kk] = bg+Nphot*GI[kk];
#if verbose>2
  printf("GI:\n");
  showRegion(GI, Rw);
#endif

  /*
   * from LL2PG.m
   *   model = x(3)+x(4)*gaussianInt2([x(1), x(2)], x(5), (size(patch, 1)-1)/2);
   *   mask = disk2d((size(patch,1)-1)/2);
   *   %L = -sum(sum(-(patch-model).^2./model - .5*log(model)));
   *   L = -sum(sum(mask.*(-(patch-model).^2./model - .5*log(model))));
   */

  double E = 0;
  for (size_t kk=0; kk<Rw; kk++)
  {
    E+= (GI[kk]-R[kk])*(GI[kk]-R[kk])/GI[kk] - .5*log(GI[kk]);
  }
  //  E = -E;

  if(0){
    printf("R=[");
    for (size_t kk=0; kk<Rw; kk++)
    {
      printf("%f ", R[kk]);
    }
    printf("]\n");
    printf("G = [");
    for (size_t kk=0; kk<Rw; kk++)
    {
      printf("%f ", GI[kk]);
    }
    printf("]\n");
    printf("E: %f, mu: %f\n", E, mu[0]);
  }

  /* Quadratic
   * for (size_t kk=0; kk<Rw*Rw; kk++)
   * E+= (GI[kk]-R[kk])*(GI[kk]-R[kk]);
   */

  if(!isfinite(E))
    E = 10e99;

  return E;
}

double lxy (const gsl_vector *v, void *params)
  // likelihood in xy, function to optimize
{

  double x, y;

  optParams *p = (optParams *) params;

  size_t Rw = p->Rw;
  double * R = p->R;
  double * GI = p->G;
  double sigma = p->sigma;
  double bg = p->bg;

#if verbose > 2
  printf("Rw: %lu\n", Rw);
  printf("Nphot: %lu\n", Nphot);
  printf("sigma: %f\n", sigma);
  printf("bg   : %f\n", bg);
#endif

  // Get the other parameters ...
  x = gsl_vector_get(v, 0);
  y = gsl_vector_get(v, 1);
  double Nphot = gsl_vector_get(v,2);

  // Create Gaussian ...
  double mu[] = {x,y};
  gaussianInt2(GI, mu, &sigma, Rw);
  for(size_t kk = 0; kk<Rw*Rw; kk++)
    GI[kk] = bg+Nphot*GI[kk];
#if verbose>2
  printf("GI:\n");
  showRegion(GI, Rw);
#endif

  /*
   * from LL2PG.m
   *   model = x(3)+x(4)*gaussianInt2([x(1), x(2)], x(5), (size(patch, 1)-1)/2);
   *   mask = disk2d((size(patch,1)-1)/2);
   *   %L = -sum(sum(-(patch-model).^2./model - .5*log(model)));
   *   L = -sum(sum(mask.*(-(patch-model).^2./model - .5*log(model))));
   */

  double E = 0;
  for (size_t kk=0; kk<Rw*Rw; kk++)
    E+= (GI[kk]-R[kk])*(GI[kk]-R[kk])/GI[kk] - .5*log(GI[kk]);
  //  E = -E;

  /* Quadratic
   * for (size_t kk=0; kk<Rw*Rw; kk++)
   * E+= (GI[kk]-R[kk])*(GI[kk]-R[kk]);
   */
//  printf("E = %f\n", E);
  return E;
}

void lxy_df (const gsl_vector *v, void *params, gsl_vector * df)
{
  double h = 1.49e-8; // step size
  gsl_vector * dv = gsl_vector_alloc(3);
  gsl_vector_memcpy(dv, v);

  for(int kk = 0; kk<3; kk++)
  {
    gsl_vector_set(dv, kk, gsl_vector_get(v, kk) + h); // Set
    gsl_vector_set(df, kk, (lxy(dv, params)-lxy(v, params))/h);
   // printf("%f ", gsl_vector_get(df, kk));
    gsl_vector_set(dv, kk, gsl_vector_get(v, kk)); // Reset
  }
//  printf("\n");

  gsl_vector_free(dv);

}

void lxy_fdf (const gsl_vector *v, void *params, double * f, gsl_vector *df)
{
  *f = lxy(v, params);
  lxy_df(v, params, df);
}

int localizeDotXY(double * V, size_t Vm,
    double * D,  double * F, double sigma)
  // Localization for a dot roughly centered in V of size Vm x Vm
  // D[0], D[1], D[3] are the global coordinates of the dot
  // F are the fitted coordinates
{

  /* Parameters that will not be optimized */
  optParams par;
  par.R = V;
  par.Rw = Vm;
  par.sigma = sigma;
  par.bg = estimateBG(V, Vm);
  par.G = (double *) malloc(Vm*Vm*sizeof(double));

  const gsl_multimin_fdfminimizer_type *T =
    gsl_multimin_fdfminimizer_vector_bfgs;
  //const gsl_multimin_fminimizer_type *T =
    //  gsl_multimin_fminimizer_nmsimplex2;

  gsl_multimin_fdfminimizer *s = NULL;
  gsl_vector *ss, *x;
  gsl_multimin_function_fdf minex_func;

  size_t iter = 0;
  int status;
  double size;

  /* Starting point for parameters to be optimized */
  x = gsl_vector_alloc (3);
  gsl_vector_set(x, 0, 0); // x position
  gsl_vector_set(x, 1, 0); // y position
  gsl_vector_set(x, 2, estimateNphot(V, Vm)); // Nphot

  /* Grid search to find a reasonable starting point
   * for (x,y) to avoid local minima */

  if(1) {
  double xmin = 0;
  double ymin = 0;
  double emin = 10e99;
  int ng = 5; // ng*ng grid points
  for(double dx = -.4; dx<.4; dx=dx+0.8/ng)
    for(double dy = -.4; dy<.4; dy=dy+0.8/ng)
    {
      gsl_vector_set(x, 0, dx); // x position
      gsl_vector_set(x, 1, dy); // y position

      double err = lxy(x, (void *) &par);
//    printf("%f %f :  %f\n", dx, dy, err);
      if( err < emin)
      {
        emin = err;
        xmin = dx;
        ymin = dy;
      }
    }
   gsl_vector_set(x, 0, xmin); // x position
   gsl_vector_set(x, 1, ymin); // y position
 //  printf("Initial position: %f %f\n", xmin, ymin);
  }

  /* Set initial step sizes */
  ss = gsl_vector_alloc(3);
  gsl_vector_set_all(ss, 0.1);

  /* Initialize method and iterate */
  minex_func.n = 3;
  minex_func.f = lxy; // Error
  minex_func.df = lxy_df; // Derivative
  minex_func.fdf = lxy_fdf; // Error and Derivative
  minex_func.params = ∥

  s = gsl_multimin_fdfminimizer_alloc (T, 3);
  gsl_multimin_fdfminimizer_set (s, // Optimizer
      &minex_func, // Function to optimize
      x, // Start point
      0.1, // Default step size of the line search
      0.1); // tol. 0 = exact line search 0.1 = recommended from GSL manual

  /* tol=0.1 was about twice as fast as tol=0 for this problem with no noticeable
   * regression in accuracy */

size_t maxIterationsQN = 500; // abort if more iterations

  do
  {
    iter++;
    status = gsl_multimin_fdfminimizer_iterate(s);

    if (status)
      break;

    size = 0; //gsl_multimin_fdfminimizer_size(s);
    status = gsl_multimin_test_gradient (s->gradient, 1e-4);

    if (status == GSL_SUCCESS)
    {
#if verbose > 0
      printf ("converged to minimum at\n");
      printf ("%5lu x:%10.3e y:%10.3e NP:%6.1f f() = %7.3f size = %10.3e\n",
          iter,
          gsl_vector_get (s->x, 0),
          gsl_vector_get (s->x, 1),
          gsl_vector_get (s->x, 2),
          0.0, size);
#endif
    }
  }
  while (status == GSL_CONTINUE && iter < maxIterationsQN);

  // F is always written even if the convergence failed
  F[0] = gsl_vector_get (s->x, 0)+nearbyint(D[0]);
  F[1] = gsl_vector_get (s->x, 1)+nearbyint(D[1]);
  F[2] = D[2]; // Not optimized

  gsl_vector_free(x);
  gsl_vector_free(ss);
  gsl_multimin_fdfminimizer_free (s);

  free(par.G);

  return status;
}

int localizeDotZ(double * V, size_t Vm,
    double * D,  double * F, double sigma)
  // Localization for a dot roughly centered in V of size Vm x Vm
  // D[0], D[1], D[3] are the global coordinates of the dot
  // F are the fitted coordinates
{

  // Non-optimized parameters
  optParams par;
  par.R = V;
  par.Rw = Vm;
  par.sigma = sigma;
  par.bg = estimateBG(V, Vm);
  par.G = (double *) malloc(Vm*sizeof(double));

  const gsl_multimin_fminimizer_type *T =
    gsl_multimin_fminimizer_nmsimplex2;
  gsl_multimin_fminimizer *s = NULL;
  gsl_vector *ss, *x;
  gsl_multimin_function minex_func;

  size_t iter = 0;
  int status;
  double size;

  /* Starting point */
  x = gsl_vector_alloc (2);
  gsl_vector_set(x, 0, 0); // z position
  gsl_vector_set(x, 1, estimateNphot(V, Vm)); // Nphot

  /* Set initial step sizes */
  ss = gsl_vector_alloc(2);
  gsl_vector_set_all(ss, 0.1);
  gsl_vector_set(ss, 1, 10); // Number of photons

  /* Initialize method and iterate */
  minex_func.n = 2;
  minex_func.f = lz;
  minex_func.params = ∥

  //printf(".\n"); fflush(stdout);

  s = gsl_multimin_fminimizer_alloc (T, 2);

  //  printf("..\n"); fflush(stdout);
  int error = gsl_multimin_fminimizer_set(s,
      &minex_func,
      x, // starting point
      ss); // Step sizes
  if(error>0)
    printf("error: %d\n", error);
  // printf("...\n"); fflush(stdout);

  do
  {
    iter++;
    status = gsl_multimin_fminimizer_iterate(s);
    //  printf("....\n"); fflush(stdout);

    if (status)
      break;

    size = gsl_multimin_fminimizer_size(s);
    status = gsl_multimin_test_size(size, convCriteria);
    //  printf(".....\n"); fflush(stdout);

    if (status == GSL_SUCCESS)
    {
#if verbose > 0
      printf ("converged to minimum at\n");
      printf ("%5lu z:%10.3e NP:%6.1f f() = %7.3f size = %10.3e\n",
          iter,
          gsl_vector_get (s->x, 0),
          gsl_vector_get (s->x, 1),
          s->fval, size);
#endif
    }
  }
  while (status == GSL_CONTINUE && iter < maxIterations);

  // Always Update Z position -- even if convergence failed
  F[2] = gsl_vector_get (s->x, 0)+nearbyint(D[2]);

  gsl_vector_free(x);
  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free (s);

  free(par.G);

  return status;
}
int localize(double * V, size_t Vm, size_t Vn, size_t Vp,
    double * D, size_t Dm, double * F, double sigma_xy, double sigma_z)
  // run the localization routine for a list of dots
{

  /* The window size is a very important parameter. Small windows can not capture the full signals and that results in
   * bad fittings.
   * Ws = 7 is good until sigma = 1.1
   * Ws = 9 is good at least until sigma = 1.6
   */
  int Ws = 9; // Window size
  double * W = malloc(Ws*Ws*sizeof(double));

  for(size_t kk =0; kk<Dm; kk++)
  {
    // Copy the neighbourhood around each D into W
    if(getRegion(W, Ws,
          V, Vm, Vn, Vp, D+kk*3) == 0)
    {
#if verbose > 1
      printf("W:\n");
      showRegion(W, Ws);
#endif
      // Local fitting in W
#if verbose > 0
      int status = localizeDotXY(W,  Ws, D+kk*3, F+kk*3, sigma_xy);
      printf("Status: %d\n", status);
#else
      localizeDotXY(W,  Ws, D+kk*3, F+kk*3, sigma_xy);
#endif

      if(getZLine(W,Ws,
            V, Vm, Vn, Vp, D+kk*3) == 0)
      {
        localizeDotZ(W, Ws, D+kk*3, F+kk*3, sigma_z);
      }

    }
    else
    {
      // TODO: add also Nphot and status to the output
      F[kk*3] = -1;
      F[kk*3+1] = -1;
      F[kk*3+2] = -1;
    }
  }
  free(W);
  return 0;
}

int unit_tests(){
  double * V; // image
  int Vm = 1024; int Vn = 1024; int Vp = 60;
  double * D; // list of dots
  size_t Dm = 5; // number of dots
  double * F; // fitted dots

  printf("Image size: %dx%dx%d\n", Vm, Vn, Vp);
  printf("Localizing %lu dots\n", Dm);

  V = (double *) malloc(Vm*Vn*Vp*sizeof(double));
  D = (double *) malloc(3*Dm*sizeof(double));
  F = (double *) malloc(3*Dm*sizeof(double));

  // Initialize the data
  for(int kk=0; kk<Vm*Vn*Vp; kk++)
    V[kk] = 1000+rand_range(0,0);

  for(uint32_t kk=0; kk<Dm; kk++)
  {
    size_t pos = kk*3;
    D[pos] =   rand_range(6, Vm-7);
    D[pos+1] = rand_range(6, Vn-7);
    D[pos+2] = rand_range(6, Vp-7);
    if(D[pos] < 6)
      D[pos] = 6;
    if(D[pos+1] < 6)
      D[pos+1] = 6;
    if(D[pos+2] < 6)
      D[pos+2] = 6;
#if verbose > 0
    printf("D %03d %f %f %f\n", kk, D[pos], D[pos+1], D[pos+2]);
#endif
  }

  for(uint32_t kk=0; kk<Dm; kk++)
  {
    size_t pos = D[kk*3] + D[kk*3+1]*Vm + D[kk*3+2]*Vm*Vn;
    V[pos] = V[pos]+5;
  }

  D[0] = 100; D[1] = 100; D[2] = 30;
  V[(int) D[0]+(int) D[1]*Vm+(int) D[2]*Vm*Vn] = 7;
  V[(int) D[0]+1+(int) D[1]*Vm+(int) D[2]*Vm*Vn] = 6;

  // Run the optimization
  localize(V, Vm, Vn, Vp, D, Dm, F, 1.1, 1.1);

  if(1)
  {
    for(double sigma = 0.1; sigma<2; sigma=sigma+0.1)
      localize(V, Vm, Vn, Vp, D, Dm, F, sigma, sigma);
  }

  // In next version, also supply clustering information

#if verbose >0
  for(size_t kk = 0; kk<Dm; kk++)
  {
    size_t pos = kk*3;
    printf("%6lu (%f, %f, %f) -> (%f, %f, %f)\n", kk,
        D[pos], D[pos+1], D[pos+2],
        F[pos], F[pos+1], F[pos+2]);
  }
#endif
  free(F);
  free(D);
  free(V);
  return 0;
}

int main(int argc, char ** argv)
  // For testing, not used any more, see the MATLAB interface in
  // df_mlfit.c
{

  printf("%s\n", argv[0]);

  if(argc == 1)
    return unit_tests();

  return 0;
}


text:
#pragma once

#include <algorithm>
#include <array>
#include <cstdint>
#include <initializer_list>
#include <iterator>
#include <sstream>
#include <string>

#include <absl/types/span.h>
#include <gsl/gsl>

#include "chainerx/axes.h"
#include "chainerx/constant.h"
#include "chainerx/dims.h"
#include "chainerx/error.h"

namespace chainerx {

class Strides;

class Shape : public Dims {
    using BaseVector = Dims;

public:
    using const_iterator = BaseVector::const_iterator;
    using const_reverse_iterator = BaseVector::const_reverse_iterator;
    // TODO(niboshi): Declare other types required for this class to be a container.

    Shape() = default;

    ~Shape() = default;

    // by iterators
    template <typename InputIt>
    Shape(InputIt first, InputIt last) {
        if (std::distance(first, last) > kMaxNdim) {
            throw DimensionError{"too many dimensions: ", std::distance(first, last)};
        }
        insert(begin(), first, last);
    }

    // by span
    explicit Shape(absl::Span<const int64_t> dims) : Shape{dims.begin(), dims.end()} {}

    // by initializer list
    Shape(std::initializer_list<int64_t> dims) : Shape{dims.begin(), dims.end()} {}

    // copy
    Shape(const Shape&) = default;
    Shape& operator=(const Shape&) = default;

    // move
    Shape(Shape&&) = default;
    Shape& operator=(Shape&&) = default;

    int64_t GetTotalSize() const;

    std::string ToString() const;

    int8_t ndim() const noexcept { return gsl::narrow_cast<int8_t>(size()); }

    const int64_t& operator[](int8_t index) const {
        if (!(0 <= index && static_cast<size_t>(index) < size())) {
            throw IndexError{"Shape index ", index, " out of bounds for shape with ", size(), " size."};
        }
        return this->Dims::operator[](index);
    }

    int64_t& operator[](int8_t index) {
        if (!(0 <= index && static_cast<size_t>(index) < size())) {
            throw IndexError{"Shape index ", index, " out of bounds for shape with ", size(), " size."};
        }
        return this->Dims::operator[](index);
    }

    // span
    absl::Span<const int64_t> span() const { return {*this}; }
};

namespace internal {

bool IsContiguous(const Shape& shape, const Strides& strides, int64_t item_size);

// Returns true if a reduction can take place under the given conditions, only considering the number of dimensions.
// Otherwise, returns false.
//
// TODO(hvy): Check the dimension lengths too and reconsider the interface. E.g. return void and assert inside the function if only used for
// assertions.
bool IsValidReductionShape(const Shape& in_shape, const Axes& axes, const Shape& out_shape, bool allow_keepdims);

int64_t CountItemsAlongAxes(const Shape& shape, const Axes& axes);

Shape BroadcastShapes(const Shape& shape0, const Shape& shape1);

// Returns a shape where axes are reduced.
Shape ReduceShape(const Shape& shape, const Axes& axes, bool keepdims);

// Returns a shape with additional axes, with length 1.
Shape ExpandShape(const Shape& shape, const Axes& axes);

Shape TransposeShape(const Shape& shape, const Axes& axes);

}  // namespace internal

std::ostream& operator<<(std::ostream& os, const Shape& shape);

void CheckEqual(const Shape& lhs, const Shape& rhs);

}  // namespace chainerx

max_stars_repo_name: zaltoprofen/chainer

max_issues_repo_name: zaltoprofen/chainer

max_forks_repo_name: zaltoprofen/chainer


text:
// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Lars Nilse $
// $Authors: Steffen Sass, Holger Plattfaut, Bastian Blank $
// --------------------------------------------------------------------------

#ifndef OPENMS_FILTERING_DATAREDUCTION_SILACFILTER_H
#define OPENMS_FILTERING_DATAREDUCTION_SILACFILTER_H

#include <OpenMS/KERNEL/StandardTypes.h>
#include <OpenMS/FILTERING/DATAREDUCTION/SILACFiltering.h>
#include <OpenMS/FILTERING/DATAREDUCTION/SILACPattern.h>
#include <OpenMS/FILTERING/DATAREDUCTION/IsotopeDistributionCache.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <queue>
#include <list>

namespace OpenMS
{
  /**
   * @brief Filter to use for SILACFiltering
   *
   * A SILACFilter searches for SILAC patterns, which correspond to the defined mass shifts and charge.
   * Only peaks are taken into account, which were not blacklisted by other filters before e.g. are not part of
   * a SILAC pair yet.
   *
   * @see SILACFiltering
   */
  class OPENMS_DLLAPI SILACFilter
  {
private:
    friend class SILACFiltering;

    typedef IsotopeDistributionCache::TheoreticalIsotopePattern TheoreticalIsotopePattern;

    /**
     * @brief mass shift(s) in [Da] to search for
     */
    std::vector<DoubleReal> mass_separations_;

    /**
     * @brief charge of the ions to search for
     */
    Int charge_;

    /**
     * @brief maximal value of which a predicted SILAC feature may deviate from the averagine model
     */
    DoubleReal model_deviation_;

    /**
     * @brief number of peaks per peptide to search for
     */
    Size isotopes_per_peptide_;

    /**
     * @brief minimal intensity of SILAC features
     */
    DoubleReal intensity_cutoff_;

    /**
     * @brief minimal intensity correlation between regions of different peaks
     */
    DoubleReal intensity_correlation_;

    /**
     * @brief flag for missing peaks
     */
    bool allow_missing_peaks_;

    /**
     * Isotope distributions
     */
    static IsotopeDistributionCache * isotope_distribution_;

    /**
     * @brief number of peptides [i.e. number of labelled peptides +1, e.g. for SILAC triplet =3]
     */
    Size number_of_peptides_;

    /**
     * @brief peak positions of SILAC pattern
     */
    std::vector<DoubleReal> peak_positions_;

    /**
     * @brief m/z separation between individual peptides [e.g. {0 Th, 4 Th, 5 Th}]
     */
    std::vector<DoubleReal> mz_peptide_separations_;

    /**
     * @brief m/z shifts relative to mono-isotopic peak of unlabelled peptide
     */
    std::vector<DoubleReal> expected_mz_shifts_;

    /**
     * @brief distance between isotopic peaks of a peptide in [Th]
     */
    DoubleReal isotope_distance_;

    /**
     * @brief holds the recognized features
     */
    std::vector<SILACPattern> elements_;

    /**
     * @brief m/z at which the filter is currently applied to
     */
    DoubleReal current_mz_;

    /**
     * @brief exact m/z shift of isotopic peaks in a SILAC pattern relative to the mono-isotopic peak of the light peptide, peptides (row) x isotope (column)
     */
    std::vector<std::vector<DoubleReal> > exact_shifts_;

    /**
     * @brief m/z positions mz + exact_shifts in a SILAC pattern, where mz is the m/z of the mono-isotopic peak of light peptide
     */
    std::vector<std::vector<DoubleReal> > exact_mz_positions_;

    /**
     * @brief intensities at mz + exact_shifts in a SILAC pattern, where mz is the m/z of the mono-isotopic peak of light peptide
     */
    std::vector<std::vector<DoubleReal> > exact_intensities_;

    /**
     * @brief expected m/z shift of isotopic peaks in a SILAC pattern relative to the mono-isotopic peak of the light peptide, peptides (row) x isotope (column)
     */
    std::vector<std::vector<DoubleReal> > expected_shifts_;

    /**
     * @brief Checks if there exists a SILAC feature at the given position in the raw (interpolated) data, which corresponds to the filter's properties
     * @param rt RT value of the position
     * @param mz m/z value of the position
     */
    bool isSILACPattern_(const MSSpectrum<Peak1D> &, const SILACFiltering::SpectrumInterpolation &, DoubleReal mz, DoubleReal picked_mz, const SILACFiltering &, MSSpectrum<Peak1D> & debug, SILACPattern & pattern);

    /**
     * @brief Checks if there exists a SILAC feature at the given position in the picked data
     */
    bool isSILACPatternPicked_(const MSSpectrum<Peak1D> &, DoubleReal mz, const SILACFiltering &, MSSpectrum<Peak1D> & debug);

    /**
     * @brief Extracts mass shifts and intensities from the raw (interpolated) data
     */
    bool extractMzShiftsAndIntensities_(const MSSpectrum<Peak1D> &, const SILACFiltering::SpectrumInterpolation &, DoubleReal mz, DoubleReal picked_mz, const SILACFiltering &);

    /**
     * @brief Extracts mass shifts and intensities from the picked data
     */
    bool extractMzShiftsAndIntensitiesPicked_(const MSSpectrum<Peak1D> &, DoubleReal mz, const SILACFiltering &);

    /**
     * @brief Extracts mass shifts and intensities from the picked data and returns pattern information
     */
    bool extractMzShiftsAndIntensitiesPickedToPattern_(const MSSpectrum<Peak1D> &, DoubleReal mz, const SILACFiltering &, SILACPattern & pattern);

    /**
     * @brief Checks all peaks against intensity cutoff
     */
    bool intensityFilter_();

    /**
     * @brief Checks peak form correlation between peaks of one isotope
     */
    bool correlationFilter1_(const SILACFiltering::SpectrumInterpolation &, DoubleReal mz, const SILACFiltering &);

    /**
     * @brief Checks peak form correlation between peaks of different isotopes
     */
    bool correlationFilter2_(const SILACFiltering::SpectrumInterpolation &, DoubleReal mz, const SILACFiltering &);

    /**
     * @brief Checks peak intensities against the averagine model
     */
    bool averageneFilter_(DoubleReal mz);

public:
    /**
     * @brief detailed constructor for SILAC pair filtering
     * @param mass_separations all mass shifts of the filter
     * @param charge charge of the ions to search for
     * @param model_deviation maximum deviation from the averagine model
     * @param isotopes_per_peptide number of peaks per peptide to search for
     * @param intensity_cutoff ...
     * @param intensity_correlation minimal intensity correlation between regions of different peaks
     * @param allow_missing_peaks flag for missing peaks
     */
    SILACFilter(std::vector<DoubleReal> mass_separations, Int charge, DoubleReal model_deviation, Int isotopes_per_peptide,
                DoubleReal intensity_cutoff, DoubleReal intensity_correlation, bool allow_missing_peaks);

    /**
     * @brief gets the m/z values of all peaks , which belong the last identified feature
     */
    std::vector<DoubleReal> getPeakPositions();

    /**
     * @brief gets the m/z shifts relative to mono-isotopic peak of unlabelled peptide
     */
    const std::vector<DoubleReal> & getExpectedMzShifts();

    /**
     * @brief returns all identified elements
     */
    std::vector<SILACPattern> & getElements();

    /**
     * @brief returns the charge of the filter
     */
    Int getCharge();

    /**
     * @brief returns the mass shifts of the filter in [Da]
     */
    std::vector<DoubleReal> & getMassSeparations();
  };
}

#endif /* SILACFILTER_H_ */

max_stars_repo_name: open-ms/all-svn-branches

max_issues_repo_name: kreinert/OpenMS

max_forks_repo_name: kreinert/OpenMS


text:
/*
 *  hubble.c
 *
 *  written by christian wagner
 *
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>

#define VC 299792.458
#define TWOPI 6.283185307179586
#define PI    3.141592653589793
#define Tcmb  2.725
#define pow2(x) ((x)*(x))
#define pow3(x) ((x)*(x)*(x))
#define pow4(x) ((x)*(x)*(x)*(x))

#include "hubble.h"

double hubble(struct cosmo mycosmo, double a){
  double H0=mycosmo.hub*100.;
  double Omega_m = mycosmo.wm/pow2(mycosmo.hub);
  double Omega_r  = 4.15e-5/pow2(mycosmo.hub);
  double Omega_X  = 1.-Omega_m;
  double Omega_k = 0.;

  return  H0*sqrt(Omega_r/pow4(a)+Omega_m/pow3(a)+Omega_k/pow2(a)+Omega_X*pow(a,-3*(1+mycosmo.w0+mycosmo.wa))*exp(-3.*mycosmo.wa*(1.-a)));

}

double co_distance_int(double a,void * params){
    double t;
    struct cosmo mycosmo = *(struct cosmo *) params;

    return  1./(a*a*hubble(mycosmo,a));

}

double co_distance(struct cosmo mycosmo, double a){

  gsl_function  F;
   double result;
   double epsabs=0.001;
   double epsrel=0.001;
   double abserr;
   size_t neval;

   F.function = &co_distance_int;
   F.params = &mycosmo;

   gsl_integration_qng (&F, a, 1., epsabs, epsrel, &result, &abserr, &neval);

   return result*VC*mycosmo.hub;

}

double ang_dist(struct cosmo mycosmo, double a){
  return co_distance(mycosmo,a); // assuming flatness
}

int
ode_growth (double a, const double y[], double f[], void *params)
{

  struct cosmo mycosmo = * (struct cosmo *) params;

  double w_a=mycosmo.w0+(1-a)*mycosmo.wa;
  double Omega_m=mycosmo.wm/pow2(mycosmo.hub);
  double Omega_X=1-Omega_m;
  double Xvar=(1-Omega_m)*pow(a,-3.*(1+mycosmo.w0+mycosmo.wa))*exp(-3.*mycosmo.wa*(1.-a))/(Omega_m/(a*a*a)+(1-Omega_m)*pow(a,-3.*(1+mycosmo.w0+mycosmo.wa))*exp(-3.*mycosmo.wa*(1.-a)));

    f[0] = y[1];
    f[1] = -1.5*(1.-w_a)*Xvar*y[0]/(a*a) - (3.5-1.5*w_a*Xvar)*y[1]/a ;

    return GSL_SUCCESS;
}

int
jac_growth (double a, const double y[], double *dfdy,
     double dfdt[], void *params)
{
    return GSL_SUCCESS;
}

double growth(struct cosmo mycosmo, double a){

    const gsl_odeiv_step_type * T
	= gsl_odeiv_step_rkf45;

    gsl_odeiv_step * s
	= gsl_odeiv_step_alloc (T, 2);
    gsl_odeiv_control * c
	= gsl_odeiv_control_y_new (1e-12, 0.0);
    gsl_odeiv_evolve * e
	= gsl_odeiv_evolve_alloc (2);

    double mu = 10;
    gsl_odeiv_system sys = {ode_growth, jac_growth, 2, &mycosmo};

    double  t1=a;
    double  t = 1.E-12;

    double d1;

    double h = 1e-12;
    double y[2] = { 1.0, 0.0 };

    int status;

    while (t < t1){
      status = gsl_odeiv_evolve_apply (e, c, s,
					   &sys,
					   &t, t1,
					   &h, y);
      if  (h>0.001)
	h=0.001;

      if (status != GSL_SUCCESS)
	break;

    }

    d1=y[0]*(t);

    gsl_odeiv_evolve_free (e);
    gsl_odeiv_control_free (c);
    gsl_odeiv_step_free (s);

    return d1;
}

double z_lastscattering(struct cosmo mycosmo){
  /*
  z_lastscattering(cosmo):
  Returns z_LS from Hu & White, DampingTail paper.
  */
  double wm = mycosmo.wm;
  double wb = mycosmo.wb;
  double b1 = 0.0783*pow(wb,(-0.238))/(1+39.5*pow(wb,0.763));
  double b2 = 0.560/(1+21.1*pow(wb,1.81));
  double zls= 1048.*(1+0.00124*pow(wb,(-0.738)))*(1+b1*pow(wm,b2));
  return zls;
}

double soundhorizon(struct cosmo mycosmo){
  /*
  soundhorizon(cosmo):
  A fit to the sound horizon, in Mpc, from Hu & Sugiyama (1995), Eq. B6
  */
  double wm = mycosmo.wm;
  double wb = mycosmo.wb;
  double wg = 2.4888e-5*pow4(Tcmb/2.73);
  double r0 = 0.75*wb/wg;
  double zeq= 5464.0*(wm/0.135)/pow4(Tcmb/2.725)/(1+0.6851)-1.0;
  double req= r0/(1.+zeq);
  double zLS= z_lastscattering(mycosmo);
  double rls= r0/(1+zLS);
  double tmp= (sqrt(1.+rls)+sqrt(rls+req))/(1.+sqrt(req));
  tmp= 3997.0*sqrt(wg/wm/wb)*log(tmp);
  return(tmp);
}

double distls(struct cosmo mycosmo){
  /*
  distls(cosmo):
  Returns the distance to LS, in Mpc.
  */
  double zLS = z_lastscattering(mycosmo);
  double dLS = ang_dist(mycosmo,1./(1.+zLS))/mycosmo.hub;
  return(dLS);
}

double solvehh(double dLS,struct cosmo mycosmo){
  /*
  solvehh(dLS,cosmo):
  Solves for h given the other cosmological parameters and the distance
  to last scattering (in Mpc).
  */
  struct cosmo hmin = mycosmo;
  hmin.hub = 0.3;
  // printf("min_cosmo: %g %g %g %g %g\n",hmin.hub,hmin.w0,hmin.wa,hmin.wm,hmin.wb);
  double ymin = distls(hmin)-dLS;
  //  printf("ymin: %g\n",ymin);
  struct cosmo hmax = mycosmo;
  hmax.hub = 1.0;
  double ymax = distls(hmax)-dLS;
  //  printf("ymax: %g\n",ymax);
  double ymid;
  struct cosmo hmid = mycosmo;
  while (fabs(hmax.hub-hmin.hub)>1e-3){
    hmid.hub = (hmax.hub+hmin.hub)/2.0;
    ymid = distls(hmid)-dLS;
    if (ymin*ymid<0){
	hmax = hmid;
	ymax = ymid;
      } else {
      hmin = hmid;
      ymin = ymid;
    }
  }
  //  printf("hubble from CMB: %g\n",hmid.hub);
  return hmid.hub;
}

//void getH0fromCMB(double *xstar, double myh, double *rs, double *z_lss, double *d_lss, double *h_CMB, int *writeout, char *fname){
void getH0fromCMB(double *xstar, double *stuff){

  FILE *fp;

  struct cosmo mycosmo={ 0.72, xstar[3], 0, xstar[1], xstar[0]};

  double rs=soundhorizon(mycosmo);
  stuff[0] = rs;

  double z_lss=z_lastscattering(mycosmo);
  stuff[1] = z_lss;

  double d_lss=302.4*rs/PI;
  stuff[2] = d_lss;

  double hubble_cmb= solvehh(d_lss,mycosmo);
  stuff[3] = hubble_cmb;

}

// Linker function for use with Fortran

void geth0fromcmb_(double *xstar, double *stuff){
  void getH0fromCMB();
  getH0fromCMB(xstar,stuff);
}


text:
#ifndef __PAR_EQUILIBRIUM_IC_H__
#define __PAR_EQUILIBRIUM_IC_H__

#include "GAMER.h"
#include "TestProb.h"
#include "vector"
#include <iostream>
#include <fstream>
#include <sstream>
using namespace std;

//gsl library
#ifdef SUPPORT_GSL
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#endif

typedef struct Filename_Parameter{
   int Cloud_Num;
   vector<string> Params_Filenames;
}FP;
typedef struct Physical_Parameter{
   int      Cloud_Num;

   string Params_Filenames;

   char     Cloud_Type[MAX_STRING];
   char     Density_Table_Name[MAX_STRING];
   int      AddExtPot;
   char     ExtPot_Table_Name[MAX_STRING];
   int      Cloud_RSeed;
   double   Cloud_Rho0;
   double   Cloud_R0;
   double   Cloud_MaxR;
   double*  Cloud_Center;
   double*  Cloud_BulkVel;
   double   Cloud_Einasto_Power_Factor;
   double   Cloud_Par_Num;

   int      Cloud_MassProfNBin;
}PhysP;

class Par_EquilibriumIC
{
   public:
      Par_EquilibriumIC();
      virtual ~Par_EquilibriumIC();
      void Read_Filenames( const char *filename_para);
      void Load_Physical_Params(const FP filenames,const int cloud_idx, const long NPar_AllRank);
      void Init();
      void Par_SetEquilibriumIC(real *Mass_AllRank, real *Pos_AllRank[3], real *Vel_AllRank[3],const long Par_Idx);

      PhysP params;
      FP   filenames;
   protected:

   private:
      // Derive physical attributes for particles
      double Set_Mass( double x );
      double Set_Density( double x );
      double Set_Velocity(const double x);

      // Initialize physical parameter tables
      void Init_Mass();
      void Init_Pot();
      void Init_Prob_Dens();

      //  Initialization through Table
      void Init_Mass_Table();
      void Init_Pot_Table();

      //  Add External Potential
      void Add_Ext_Pot();

      // Auxiliary functions
      int Aux_CountRow( const char *filename );
      int Aux_Countcolumn( const char *filename );
      int GetParams( const char *filename,const char *keyword,const int para_num,const char *para_type,vector <string> &container);
      void Check_InputFileName();
      void RanVec_FixRadius( const double r, double RanVec[] );

      // Solve Eddington's equation
      double potential(const double x);
      double inverse_psi_to_index (double psi);
      double integration_eng_base(double eng);

      double delta;
      double eng_min;
      double *prob_dens;
      double *int_prob_dens;
      double *psi;

      // statistics
      double slope(double* a,double* b,int start,int fin);
      void smooth_all(double* x,int start,int fin);
      double ave(double* a,int start,int fin);
      double var_n(double* a,int start,int fin);
      double cor(double* x,double* y,int start,int fin);
      void mask(double* x,int start,int fin);
      void add_num(double* x,int start,int fin);

      // Tables of particles' attributes
      double *Table_r;
      double *Table_Enclosed_Mass;
      double *Table_Density;
      double *Table_dRho_dr;
      double *Table_dRho_dx;
      double *Table_Gravity_Field;
      double *Table_Gravity_Potential;

      // Random number generator
      RandomNumber_t *Random_Num_Gen ;
};

#endif //__PAR_EQUILIBRIUM_IC_H__

max_stars_repo_name: CrazyAncestor/gamer

max_issues_repo_name: CrazyAncestor/gamer

max_forks_repo_name: CrazyAncestor/gamer


text:
#pragma once
#include "typetraits.h"

#if defined(__has_include) && __has_include(<version>)
#include <version>
#if defined(__cpp_lib_span) && __cpp_lib_span >= 202002L
#include <span>
#else
#include <gsl/span>
#endif    // __cpp_lib_span >= 202002L
#else
#include <gsl/span>
#endif    //__has_include(<version>)

#include "uuid.h"

#include <array>
#include <chrono>
#include <compare>

template <typename T> struct UuidObjectT;

template <typename T> struct UuidBasedId
{
    static UuidBasedId<T> Create()
    {
        UuidBasedId<T>     uuid;
        std::random_device rd;
        auto               seed_data = std::array<unsigned, std::mt19937::state_size>{};
        std::generate(std::begin(seed_data), std::end(seed_data), std::ref(rd));
        std::seed_seq                seq(std::begin(seed_data), std::end(seed_data));
        std::mt19937                 generator(seq);
        uuids::uuid_random_generator gen{generator};
        uuid.uuid = gen();
        return uuid;
    }

    template <typename TStr> static constexpr UuidBasedId<T> FromString(TStr const& str)
    {
        UuidBasedId<T> uuid;
        uuid.guid = uuids::uuid(str);
    }

    public:
    UuidBasedId() = default;

    uuids::uuid uuid;

    constexpr bool                  Empty() const { return uuid == Invalid().uuid; }
    static constexpr UuidBasedId<T> Invalid() { return UuidBasedId<T>(); }
    friend UuidObjectT<T>;

    friend std::hash<UuidBasedId<T>>;
};

template <typename T> struct UuidObjectT
{
    using Id = UuidBasedId<T>;
    UuidObjectT() {}
    Id id = Id::Create();
};

template <> struct Stencil::TypeTraits<uuids::uuid>
{
    using Categories = std::tuple<Stencil::Category::Primitive>;
};

template <typename T> struct Stencil::TypeTraits<UuidBasedId<T>>
{
    using Categories = std::tuple<Stencil::Category::Primitive>;
};

max_stars_repo_name: ankurvdev/stencil

max_issues_repo_name: ankurvdev/stencil

max_forks_repo_name: ankurvdev/stencil


text:
/**
 *
 * @file qwrapper_dlacpy.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Julien Langou
 * @author Henricus Bouwmeester
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:56 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlacpy(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int m, int n, int nb,
                       const double *A, int lda,
                       double *B, int ldb)
{
    DAG_CORE_LACPY;
    QUARK_Insert_Task(quark, CORE_dlacpy_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &m,     VALUE,
        sizeof(int),                        &n,     VALUE,
        sizeof(double)*nb*nb,    A,             INPUT,
        sizeof(int),                        &lda,   VALUE,
        sizeof(double)*nb*nb,    B,             OUTPUT,
        sizeof(int),                        &ldb,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlacpy_quark = PCORE_dlacpy_quark
#define CORE_dlacpy_quark PCORE_dlacpy_quark
#endif
void CORE_dlacpy_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int M;
    int N;
    const double *A;
    int LDA;
    double *B;
    int LDB;

    quark_unpack_args_7(quark, uplo, M, N, A, LDA, B, LDB);
    LAPACKE_dlacpy_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        M, N, A, LDA, B, LDB);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlacpy_f1(Quark *quark, Quark_Task_Flags *task_flags,
                          PLASMA_enum uplo, int m, int n, int nb,
                          const double *A, int lda,
                          double *B, int ldb,
                          double *fake1, int szefake1, int flag1)
{
    DAG_CORE_LACPY;
    if ( fake1 == B ) {
        QUARK_Insert_Task(quark, CORE_dlacpy_quark, task_flags,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(int),                        &m,     VALUE,
            sizeof(int),                        &n,     VALUE,
            sizeof(double)*nb*nb,    A,             INPUT,
            sizeof(int),                        &lda,   VALUE,
            sizeof(double)*nb*nb,    B,             OUTPUT | flag1,
            sizeof(int),                        &ldb,   VALUE,
            0);
    }
    else {
        QUARK_Insert_Task(quark, CORE_dlacpy_f1_quark, task_flags,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(int),                        &m,     VALUE,
            sizeof(int),                        &n,     VALUE,
            sizeof(double)*nb*nb,    A,             INPUT,
            sizeof(int),                        &lda,   VALUE,
            sizeof(double)*nb*nb,    B,             OUTPUT,
            sizeof(int),                        &ldb,   VALUE,
            sizeof(double)*szefake1, fake1,         flag1,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlacpy_f1_quark = PCORE_dlacpy_f1_quark
#define CORE_dlacpy_f1_quark PCORE_dlacpy_f1_quark
#endif
void CORE_dlacpy_f1_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int M;
    int N;
    const double *A;
    int LDA;
    double *B;
    int LDB;
    void *fake1;

    quark_unpack_args_8(quark, uplo, M, N, A, LDA, B, LDB, fake1);
    LAPACKE_dlacpy_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        M, N, A, LDA, B, LDB);
}


text:
#pragma once

#include <type_traits>

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

namespace thrustshift {

namespace kernel {

template <typename MapT, typename SrcT, typename DstT>
__global__ void gather(gsl_lite::span<const MapT> map, gsl_lite::span<const SrcT> src, gsl_lite::span<DstT> dst) {

	const auto gtid = threadIdx.x + blockIdx.x * blockDim.x;
	if (gtid < src.size()) {
		dst[gtid] = src[map[gtid]];
	}
}

} // namespace kernel

namespace async {

template <class MapRange, class SrcRange, class DstRange>
void gather(cuda::stream_t& stream, MapRange&& map, SrcRange&& src, DstRange&& dst) {
	gsl_Expects(src.size() == dst.size());
	gsl_Expects(src.size() == map.size());
	gsl_Expects(src.data() != dst.data());

	if (src.empty()) {
		return;
	}

	using map_index_type =
	    typename std::remove_reference<MapRange>::type::value_type;
	using src_value_type =
	    typename std::remove_reference<SrcRange>::type::value_type;
	using dst_value_type =
	    typename std::remove_reference<DstRange>::type::value_type;

	constexpr cuda::grid::block_dimension_t block_dim = 128;
	const cuda::grid::dimension_t grid_dim =
	    (src.size() + block_dim - 1) / block_dim;
	auto c = cuda::make_launch_config(grid_dim, block_dim);
	auto k = kernel::gather<map_index_type, src_value_type, dst_value_type>;
	cuda::enqueue_launch(k, stream, c, map, src, dst);
}

} // namespace async

} // namespace thrustshift


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

int
main (void)
{
  int i, j;
  gsl_matrix * m = gsl_matrix_alloc (10, 3);

  for (i = 0; i < 10; i++)
    for (j = 0; j < 3; j++)
      gsl_matrix_set (m, i, j, 0.23 + 100*i + j);

  for (i = 0; i < 100; i++)  /* OUT OF RANGE ERROR */
    for (j = 0; j < 3; j++)
      printf ("m(%d,%d) = %g\n", i, j,
              gsl_matrix_get (m, i, j));

  gsl_matrix_free (m);

  return 0;
}


text:
#include <bindings.cmacros.h>
#include <gsl/gsl_multifit_nlin.h>

BC_INLINE3VOID(GSL_MULTIFIT_FN_EVAL,gsl_multifit_function*,gsl_vector*,gsl_vector*)
BC_INLINE3(GSL_MULTIFIT_FN_EVAL_F,gsl_multifit_function_fdf*,gsl_vector*,gsl_vector*,int)
BC_INLINE3(GSL_MULTIFIT_FN_EVAL_DF,gsl_multifit_function_fdf*,gsl_vector*,gsl_matrix*,int)
BC_INLINE4(GSL_MULTIFIT_FN_EVAL_F_DF,gsl_multifit_function_fdf*,gsl_vector*,gsl_vector*,gsl_matrix*,int)


text:
#ifndef QDM_LOGL_H
#define QDM_LOGL_H 1

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

#include <qdm/tau.h>

int
qdm_find_tau(
    double *result,

    double v,
    size_t spline_df,
    const gsl_vector *knots,
    const gsl_vector *mmm
);

int
qdm_logl(
    double *log_likelihood,
    double *tau,

    double v,
    size_t spline_df,
    const gsl_vector *knots,
    const gsl_vector *mmm,

    double tau_low,
    double tau_high,
    double xi_low,
    double xi_high
);

int
qdm_logl_2(
    double *log_likelihood,
    double *tau,

    double x,
    double y,

    double tau_low,
    double tau_high,

    double xi_low,
    double xi_high,

    size_t spline_df,

    const gsl_matrix *theta,
    const gsl_vector *knots
);

void
qdm_logl_3(
    double *log_likelihood,
    double *tau,

    double x,
    double y,

    const qdm_tau *t,
    const gsl_vector *xi,

    const gsl_matrix *theta
);

#endif /* QDM_LOGL_H */

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#ifdef _WIN32
#else
#include <sys/resource.h>
#endif
#include <unistd.h>
#include <signal.h>
#include <gsl/gsl_rng.h>
#ifndef NOMPI
#include <mpi.h>
#endif

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

/*! \file system.c
 *  \brief contains miscellaneous routines, e.g. elapsed time measurements
 */

/*! This routine returns a random number taken from a table of random numbers,
 *  which is refilled every timestep.  This method is used to allow random
 *  number application to particles independent of the number of processors
 *  used, and independent of the particular order the particles have. In order
 *  to work properly, the particle IDs should be set properly to unique
 *  integer values.
 */
double get_random_number(int id)
{
  return RndTable[(id % RNDTABLE)];
}

/*! This routine fills the random number table.
 */
void set_random_numbers(void)
{
  int i;

  for(i = 0; i < RNDTABLE; i++)
    RndTable[i] = gsl_rng_uniform(random_generator);
}

/*! returns the number of cpu-ticks in seconds that have elapsed, or the
 *  wall-clock time obtained with MPI_Wtime().
 */
double second(void)
{
#ifdef WALLCLOCK
#ifndef NOMPI
  return MPI_Wtime();
#else
  return ((double) clock()) / CLOCKS_PER_SEC;
#endif
#else
  return ((double) clock()) / CLOCKS_PER_SEC;
#endif

  /* note: on AIX and presumably many other 32bit systems,
   * clock() has only a resolution of 10ms=0.01sec
   */
}

/*! returns the time difference between two measurements obtained with
 *  second(). The routine takes care of the possible overflow of the tick
 *  counter on 32bit systems, but depending on the system, this may not always
 *  work properly. Similarly, in some MPI implementations, the MPI_Wtime()
 *  function may also overflow, in which case a negative time difference would
 *  be returned. The routine returns instead a time difference equal to 0.
 */
double timediff(double t0, double t1)
{
  double dt;

  dt = t1 - t0;

  if(dt < 0)	/* overflow has occured (for systems with 32bit tick counter) */
    {
#ifdef WALLCLOCK
      dt = 0;
#else
      dt = t1 + pow(2, 32) / CLOCKS_PER_SEC - t0;
#endif
    }

  return dt;
}

/*! returns the maximum of two double
 */
double dmax(double x, double y)
{
  if(x > y)
    return x;
  else
    return y;
}

/*! returns the minimum of two double
 */
double dmin(double x, double y)
{
  if(x < y)
    return x;
  else
    return y;
}

/*! returns the maximum of two integers
 */
int imax(int x, int y)
{
  if(x > y)
    return x;
  else
    return y;
}

/*! returns the minimum of two integers
 */
int imin(int x, int y)
{
  if(x < y)
    return x;
  else
    return y;
}

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
#include <stdio.h>
#include <cblas.h>

#define M 2
#define N 2
#define K 3
#define X 4
#define Y 5

void print_mat(const char *name, int r, int c, float *m)
{
    printf("%s =\n", name);
    for (int i = 0; i < r; i++) {
        for (int j = 0; j < c; j++) {
            printf("%.2lf ", m[i * c + j]);
        }
        printf("\n");
    }
}

int main()
{
    float A[M * K] = {
        1, 2, 3,
        4, 5, 6,
    };
    float B[K * N] = {
        0.1, 0.2,
        0.3, 0.4,
        0.5, 0.6,
    };
    float C[X * Y] = { 0 };

    // M x K, K x N -> M x N (single precision, 's'gemm)
    // Save A x B in C starting from (1, 1)
    cblas_sgemm(
        CblasRowMajor, CblasNoTrans, CblasNoTrans,
        M, N, K,
        1,
        A, K,
        B, N,
        0,
        C + Y + 1, Y
    );

    print_mat("A", M, K, A);
    print_mat("B", K, N, B);
    print_mat("C", X, Y, C);

    return 0;
}


text:
#ifndef OPENMC_TALLIES_FILTER_MATERIAL_H
#define OPENMC_TALLIES_FILTER_MATERIAL_H

#include <cstdint>
#include <unordered_map>

#include <gsl/gsl>

#include "openmc/tallies/filter.h"
#include "openmc/vector.h"

namespace openmc {

//==============================================================================
//! Specifies which material tally events reside in.
//==============================================================================

class MaterialFilter : public Filter
{
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors

  ~MaterialFilter() = default;

  //----------------------------------------------------------------------------
  // Methods

  std::string type() const override {return "material";}

  void from_xml(pugi::xml_node node) override;

  void get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match)
  const override;

  void to_statepoint(hid_t filter_group) const override;

  std::string text_label(int bin) const override;

  //----------------------------------------------------------------------------
  // Accessors

  vector<int32_t>& materials() { return materials_; }

  const vector<int32_t>& materials() const { return materials_; }

  void set_materials(gsl::span<const int32_t> materials);

private:
  //----------------------------------------------------------------------------
  // Data members

  //! The indices of the materials binned by this filter.
  vector<int32_t> materials_;

  //! A map from material indices to filter bin indices.
  std::unordered_map<int32_t, int> map_;
};

} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_MATERIAL_H

max_stars_repo_name: cjwyett/openmc

max_issues_repo_name: cjwyett/openmc

max_forks_repo_name: cjwyett/openmc