text:
#pragma once

#include <type_traits>

#include <cuda/define_specifiers.hpp>
#include <cuda/runtime_api.hpp>

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

#include <thrustshift/fill.h>

namespace thrustshift {

namespace kernel {

template <typename SrcT, typename DstT>
__global__ void copy(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[gtid];
	}
}

template <typename SrcT, typename DstT, typename T, typename I>
__global__ void copy_find(gsl_lite::span<const SrcT> src,
                          gsl_lite::span<DstT> dst,
                          T value,
                          I* pos) {

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

} // namespace kernel

namespace async {

//! thrust uses sometimes a cudaMemcpyAsync instead of a copy kernel
template <class SrcRange, class DstRange>
void copy(cuda::stream_t& stream, SrcRange&& src, DstRange&& dst) {
	gsl_Expects(src.size() == dst.size());

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

	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::copy<src_value_type, dst_value_type>;
	cuda::enqueue_launch(k, stream, c, src, dst);
}

//! Copy and search for an element. If the element occurs more than once, it
//! is undefined which of the valid positions is returned. If the element does
//! not occur `pos` is unchanged.
template <class SrcRange, class DstRange, typename T, typename I>
void copy_find(cuda::stream_t& stream,
               SrcRange&& src,
               DstRange&& dst,
               const T& value,
               I* pos) {
	gsl_Expects(src.size() == dst.size());
	gsl_Expects(pos != nullptr);

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

	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::copy_find<src_value_type, dst_value_type, T, I>;
	cuda::enqueue_launch(k, stream, c, src, dst, value, pos);
}

} // namespace async

namespace detail {

template <int BLOCK_DIM,
          int NUM_ELEMENTS,
          int LD_FIRST = NUM_ELEMENTS,
          int LD_RESULT = NUM_ELEMENTS,
          int NUM_PER_ROW_FIRST = NUM_ELEMENTS,
          int NUM_PER_ROW_RESULT = NUM_ELEMENTS>
struct helper_t {
	template <class iteratorA_t, class iteratorB_t, class F>
	CUDA_FHD static void block_copy_even(iteratorA_t first,
	                                     iteratorB_t result,
	                                     int tid,
	                                     F f) {
		if (BLOCK_DIM < NUM_ELEMENTS) {
// pragma unroll is device code specific
#ifdef __CUDA_ARCH__
#pragma unroll
#endif
			for (int i = 0; i < NUM_ELEMENTS - BLOCK_DIM; i += BLOCK_DIM) {
				f(first,
				  result,
				  i + ((i + tid) / NUM_PER_ROW_FIRST) *
				          (LD_FIRST - NUM_PER_ROW_FIRST) + tid,
				  i + ((i + tid) / NUM_PER_ROW_RESULT) *
				          (LD_RESULT - NUM_PER_ROW_RESULT) + tid);
			}
		}
	}

	template <class iteratorA_t, class iteratorB_t, class F>
	CUDA_FHD static void block_copy_tail(iteratorA_t first,
	                                     iteratorB_t result,
	                                     int tid,
	                                     F f) {
		if (NUM_ELEMENTS % BLOCK_DIM == 0) {
			f(first,
			  result,
			  (NUM_ELEMENTS - BLOCK_DIM) +
			      (((NUM_ELEMENTS - BLOCK_DIM) + tid) / NUM_PER_ROW_FIRST) *
			          (LD_FIRST - NUM_PER_ROW_FIRST) + tid,
			  (NUM_ELEMENTS - BLOCK_DIM) +
			      (((NUM_ELEMENTS - BLOCK_DIM) + tid) / NUM_PER_ROW_RESULT) *
			          (LD_RESULT - NUM_PER_ROW_RESULT) + tid);
		}
		else if (tid + (NUM_ELEMENTS / BLOCK_DIM) * BLOCK_DIM < NUM_ELEMENTS) {
			constexpr int j = (NUM_ELEMENTS / BLOCK_DIM) * BLOCK_DIM;
			f(first,
			  result,
			  j + ((j + tid) / NUM_PER_ROW_FIRST) *
			          (LD_FIRST - NUM_PER_ROW_FIRST) + tid,
			  j + ((j + tid) / NUM_PER_ROW_RESULT) *
			          (LD_RESULT - NUM_PER_ROW_RESULT) + tid);
		}
	}
}; // helper

template <int BLOCK_DIM, int NUM_ELEMENTS>
struct helper_t<BLOCK_DIM,
                NUM_ELEMENTS,
                NUM_ELEMENTS,
                NUM_ELEMENTS,
                NUM_ELEMENTS,
                NUM_ELEMENTS> {
	template <class iteratorA_t, class iteratorB_t, class F>
	CUDA_FHD static void block_copy_even(iteratorA_t first,
	                                     iteratorB_t result,
	                                     int tid,
	                                     F f) {
		if (BLOCK_DIM < NUM_ELEMENTS) {
// pragma unroll is device code specific
#ifdef __CUDA_ARCH__
#pragma unroll
#endif
			for (int i = 0; i < NUM_ELEMENTS - BLOCK_DIM; i += BLOCK_DIM) {
				const int j = i + tid;
				f(first, result, j, j);
			}
		}
	}

	template <class iteratorA_t, class iteratorB_t, class F>
	CUDA_FHD static void block_copy_tail(iteratorA_t first,
	                                     iteratorB_t result,
	                                     int tid,
	                                     F f) {
		if (NUM_ELEMENTS % BLOCK_DIM == 0) {
			const int j = NUM_ELEMENTS - BLOCK_DIM + tid;
			f(first,
			  result,
			  j,
			  j);
		}
		else if (tid + (NUM_ELEMENTS / BLOCK_DIM) * BLOCK_DIM < NUM_ELEMENTS) {
			const int j = (NUM_ELEMENTS / BLOCK_DIM) * BLOCK_DIM + tid;
			f(first, result, j, j);
		}
	}
}; // struct helper_t

//! Specialization for equal block dimension and number of elements
template <int BD_AND_NE>
struct helper_t<BD_AND_NE,
                BD_AND_NE,
                BD_AND_NE,
                BD_AND_NE,
                BD_AND_NE,
                BD_AND_NE> {
	template <class iteratorA_t, class iteratorB_t, class F>
	CUDA_FHD static void block_copy_even(iteratorA_t first,
	                                     iteratorB_t result,
	                                     int tid,
	                                     F f) {
	}

	template <class iteratorA_t, class iteratorB_t, class F>
	CUDA_FHD static void block_copy_tail(iteratorA_t first,
	                                     iteratorB_t result,
	                                     int tid,
	                                     F f) {
		f(first, result, tid, tid);
	}

}; // struct helper_t

}; // namespace detail

/* \brief Copy data blockwise with efficient loop unrolling.
 *
 * This function implements an efficient memory copy. This is in particular
 * useful if you use only a few threads to copy data, since this is limited
 * by register dependencies if trivially implemented. The algorithm is
 * implemented in two steps. First we copy elements subsequently and
 * block wise with the whole thread block:
 *
 *     N = 11, bdim = 4: x x x x | x x x x | o o o
 *     N = 12, bdim = 4: x x x x | x x x x | o o o o
 *     N = 13, bdim = 4: x x x x | x x x x | x x x x | o
 *
 * Thus the last elements are not copied (denoted by `o`).
 * Afterwards we must copy the remaining elements, which is
 * trivial if N % bdim == 0.
 *
 * \param BLOCK_DIM Current thread block dimension
 * \param NUM_ELEMENTS total amount of elements to copy
 * \param LD_RESULT Leading dimension on result. Only used in case of two dimensional
 *        arranged data. This is equal to the number of elements, which
 *        must be skipped to get an element in the next row. This is
 *        useful if you want to skip some of the columns of the 2D data:
 *
 *           |    ld     |
 *            x x x x o o
 *            x x x x o o
 *            x x x x o o
 *            x x x x o o
 *           |  nepr |
 *
 *        `nepr` = num elements per row
 *        `ld` = leading dimension
 *        If you want to copy only the `x` elements.
 *
 * \param NUM_PER_ROW Only used in case of two dimensional data copies
 * \param first iterator to first read element
 * \param result iterator to first result element
 * \param f A functor to write the data, e.g.
 *      ```
 *      auto default_f = [](iteratorA_t first,
 *                          iteratorB_t result,
 *                          int i_first,
 *                          int i_result) {
 *          result[i_result] = first[i_first];
 *      };
 *      ```
 *     `i_first` and `i_result` are the global indices on iterator `first` and `result`. With such
 *     a functor it is possible to modify the read and write process. E.g. you can read
 *     the data in a permuted way.
 *
 */
template <int BLOCK_DIM,
          int NUM_ELEMENTS,
          class iteratorA_t,
          class iteratorB_t,
          int LD_FIRST = NUM_ELEMENTS,
          int LD_RESULT = NUM_ELEMENTS,
          int NUM_PER_ROW_FIRST = NUM_ELEMENTS,
          int NUM_PER_ROW_RESULT = NUM_ELEMENTS,
          class F>
CUDA_FHD void block_copy(iteratorA_t first,
                         iteratorB_t result,
                         int tid,
                         F f) {

	detail::helper_t<BLOCK_DIM,
	                 NUM_ELEMENTS,
	                 LD_FIRST,
	                 LD_RESULT,
	                 NUM_PER_ROW_FIRST,
	                 NUM_PER_ROW_RESULT>::block_copy_even(first,
	                                                      result,
	                                                      tid,
	                                                      f);
	detail::helper_t<BLOCK_DIM,
	                 NUM_ELEMENTS,
	                 LD_FIRST,
	                 LD_RESULT,
	                 NUM_PER_ROW_FIRST,
	                 NUM_PER_ROW_RESULT>::block_copy_tail(first,
	                                                      result,
	                                                      tid,
	                                                      f);
}

template <int BLOCK_DIM,
          int NUM_ELEMENTS,
          int LD_FIRST = NUM_ELEMENTS,
          int LD_RESULT = NUM_ELEMENTS,
          int NUM_PER_ROW_FIRST = NUM_ELEMENTS,
          int NUM_PER_ROW_RESULT = NUM_ELEMENTS,
          class iteratorA_t,
          class iteratorB_t>
CUDA_FHD void block_copy(iteratorA_t first,
                         iteratorB_t result,
                         int tid = threadIdx.x) {

	auto default_f = [](iteratorA_t first,
	                    iteratorB_t result,
	                    int i_first,
	                    int i_result) {
		result[i_result] = first[i_first];
	};

	block_copy<BLOCK_DIM,
	           NUM_ELEMENTS,
	           iteratorA_t,
	           iteratorB_t,
	           LD_FIRST,
	           LD_RESULT,
	           NUM_PER_ROW_FIRST,
	           NUM_PER_ROW_RESULT>(first, result, tid, default_f);
}

//! copy data without loops unrolled
template <class iteratorA_t, class iteratorB_t>
CUDA_FHD void block_copy(iteratorA_t first,
                        int num_elements,
                        iteratorB_t result,
                        int group_dim,
                        int ld_first,
                        int ld_result,
                        int num_per_row_first,
                        int num_per_row_result,
                        int tid) {
	auto f = [](iteratorA_t first,
	            iteratorB_t result,
	            int i_first,
	            int i_result) { result[i_result] = first[i_first]; };
	//////////
	// HEAD //
	//////////
	if (group_dim < num_elements) {
		for (int i = 0; i < num_elements - group_dim; i += group_dim) {
			f(first,
			  result,
			  i + ((i + tid) / num_per_row_first) *
			          (ld_first - num_per_row_first) + tid,
			  i + ((i + tid) / num_per_row_result) *
			          (ld_result - num_per_row_result) + tid);
		}
	}
	//////////
	// TAIL //
	//////////
	if (num_elements % group_dim == 0) {
		f(first,
		  result,
		  (num_elements - group_dim) +
		      (((num_elements - group_dim) + tid) / num_per_row_first) *
		          (ld_first - num_per_row_first) + tid,
		  (num_elements - group_dim) +
		      (((num_elements - group_dim) + tid) / num_per_row_result) *
		          (ld_result - num_per_row_result) + tid);
	}
	else if (tid + (num_elements / group_dim) * group_dim < num_elements) {
		const int j = (num_elements / group_dim) * group_dim;
		f(first,
		  result,
		  j + ((j + tid) / num_per_row_first) * (ld_first - num_per_row_first) + tid,
		  j + ((j + tid) / num_per_row_result) *
		          (ld_result - num_per_row_result) + tid);
	}
}

} // namespace thrustshift

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
#include "cones.h"
#ifdef LAPACK_LIB_FOUND
#include <cblas.h>
#include <lapacke.h>
#endif
void projectsdc(double * X, int n, Work * w);

/* in place projection (with branches) */
void projCone(double *x, Cone * k, Work * w, int iter)
{
    int i;
    int count;

    /* project onto positive orthant */
    for(i = k->f; i < k->f+k->l; ++i)
    {
        if(x[i] < 0.0) x[i] = 0.0;
        //x[i] = (x[i] < 0.0) ? 0.0 : x[i];
    }
    count = k->l+k->f;
    /* project onto SOC */
    for(i = 0; i < k->qsize; ++i)
    {
        double v1 = x[count];
        double s = calcNorm(&(x[count+1]),k->q[i]-1);
        double alpha = (s + v1)/2.0;

        if(s <= v1) { /* do nothing */ }
        else if (s <= - v1) {
            memset(&(x[count]), 0, k->q[i]*sizeof(double));
        } else {
            x[count] = alpha;
            scaleArray(&(x[count+1]), alpha/s, k->q[i]-1);
            //cblas_dscal(k->q[i]-1, alpha/s, &(x[count+1]),1);
        }
        count += k->q[i];
    }
#ifdef LAPACK_LIB_FOUND
    /* project onto PSD cone */
    for (i=0; i < k->ssize; ++i){
        projectsdc(&(x[count]),k->s[i],w);
        count += (k->s[i])*(k->s[i]);
    }
#else
    if(k->ssize > 0){
        coneOS_printf("WARNING: solving SDP, no lapack library specified in makefile!\n");
        coneOS_printf("ConeOS will return a wrong answer!\n");
    }
#endif

    /*
    * exponential cone is not self dual, if s \in K
    * then y \in K^* and so if K is the primal cone
    * here we project onto K^*, via Moreau
    */
    scaleArray(&(x[count]), -1, 3*k->ep); // x = -x;
    double r,s,t;
    int idx;
#pragma omp parallel for private(r,s,t,idx)
    for (i=0; i < k->ep; ++i) {
        idx = count + 3*i;
        r = x[idx];
        s = x[idx+1];
        t = x[idx+2];

        projExpCone(&(x[idx]));

        x[idx] -= r;
        x[idx+1] -= s;
        x[idx+2] -= t;
    }
    count += 3*k->ep;

    // exponential cone:
#pragma omp parallel for
    for (i=0; i < k->ed; ++i) {
        projExpCone(&(x[count + 3*i]));
    }
    count += 3*k->ed;
    /* project onto OTHER cones */
}

double expNewtonOneD(double rho, double y_hat, double z_hat) {
    double t = fmax( -z_hat , 1e-6);
    double f, fp;
    for (int i=0; i<100; ++i){

        f = t * (t + z_hat) / rho / rho - y_hat/rho + log(t/rho) + 1;
        fp = (2 * t + z_hat) / rho / rho + 1/t;

        t = t - f/fp;

        if (t <= -z_hat) {
            return 0;
        } else if (t <= 0) {
            return z_hat;
        } else if ( fabs(f) < 1e-9 ) {
            break;
        }
    }
    return t + z_hat;
}

void expSolveForXWithRho(double * v, double * x, double rho) {
    x[2] = expNewtonOneD(rho, v[1], v[2]);
    x[1] = (x[2] - v[2]) * x[2] / rho;
    x[0] = v[0] - rho;
}

double expCalcGrad(double * v, double * x, double rho) {
    expSolveForXWithRho(v, x, rho);
    if (x[1] <= 1e-12) {
        return x[0];
    } else {
        return x[0] + x[1] * log( x[1] / x[2] );
    }
}

void expGetRhoUb(double * v, double * x, double * ub, double * lb) {
    *lb = 0;
    *ub = 0.125;
    while(expCalcGrad(v, x, *ub) > 0) {
        *lb = *ub;
        (*ub) *= 2;
    }
}

// project onto the exponential cone, v has dimension *exactly* 3
void projExpCone(double * v) {

    double r = v[0], s = v[1], t = v[2];
    // v in cl(Kexp)
    if( (s*exp(r/s) <= t && s > 0) || (r <= 0 && s == 0 && t >= 0) ) {
        return;
    }

    // -v in Kexp^*
    if ( (-r < 0 && r*exp(s/r) <= -exp(1)*t) || (-r == 0 && -s >= 0 && -t >= 0) ) {
        memset(v, 0, 3*sizeof(double));
        return;
    }

    // special case with analytical solution
    if(r < 0 && s < 0) {
        v[1] = 0.0;
        v[2] = fmax(v[2],0);
        return;
    }

    double ub, lb, rho, g, x[3];
    expGetRhoUb(v, x, &ub, &lb);
    for(int i = 0; i < 100; ++i){
        rho = (ub + lb)/2;
        g = expCalcGrad(v, x, rho);
        if (g > 0) {
            lb = rho;
        } else{
            ub = rho;
        }
        if (ub - lb < 1e-9) {
            break;
        }
    }
    v[0] = x[0];
    v[1] = x[1];
    v[2] = x[2];
}

#ifdef LAPACK_LIB_FOUND
void projectsdc(double *X, int n, Work * w)
{ /* project onto the positive semi-definite cone */
  if (n == 1) {
    if(X[0] < 0.0) X[0] = 0.0;
    return;
  }

  int i, j, m=0;
  double * Xs = w->Xs;
  double * Z = w->Z;
  double * e = w->e;
  memcpy(Xs,X,n*n*sizeof(double));

  // Xs = X + X', save div by 2 for eigen-recomp
  for (i = 0; i < n; ++i){
    cblas_daxpy(n, 1, &(X[i]), n, &(Xs[i*n]), 1);
    //b_daxpy(n, 1, &(X[i]), n, &(Xs[i*n]), 1);
  }

  double EIG_TOL = 1e-8;
  double vupper = calcNorm(Xs,n*n);
  LAPACKE_dsyevr( LAPACK_COL_MAJOR, 'V', 'V', 'U', n, Xs, n, 0.0, vupper, -1, -1, EIG_TOL, &m, e, Z, n , NULL);
  //printf("m is %i, n is %i\n", m ,n);
  //printf("vupper is %f, max eig is %f\n",vupper, e[m>0 ? m-1:0]/2);
  memset(X, 0, n*n*sizeof(double));
  for (i = 0; i < m; ++i) {
    cblas_dsyr(CblasColMajor, CblasLower, n, e[i]/2, &(Z[i*n]), 1, X, n);
    //b_dsyr('L', n, -e[i]/2, &(Z[i*n]), 1, Xs, n);
  }
  // fill in upper half
  for (i = 0; i < n; ++i){
    for (j = i+1; j < n; ++j){
      X[i + j*n] = X[j + i*n];
    }
  }
}
#endif

ext: c

max_stars_repo_name: cvxgrp/coneos

max_issues_repo_name: cvxgrp/coneos

max_forks_repo_name: cvxgrp/coneos


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

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_siman.h>
#include <gsl/gsl_ieee_utils.h>
#include <stdio.h>

/* set up parameters for this simulated annealing run */
#define N_TRIES 200		/* how many points do we try before stepping */
#define ITERS_FIXED_T 1000	/* how many iterations for each T? */
#define STEP_SIZE 1.0		/* max step size in random walk */
#define K 1.0			/* Boltzmann constant */
#define T_INITIAL 0.008		/* initial temperature */
#define MU_T 1.003		/* damping factor for temperature */
#define T_MIN 2.0e-6

gsl_siman_params_t params = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
			     K, T_INITIAL, MU_T, T_MIN};

inline double square (double x) ;
inline double square (double x) { return x * x ; }

double test_E_1D(Element x);
void test_step_1D(const gsl_rng * r, Element *x_p, gsl_siman_params_t p);
void print_pos_1D(Element x);
double distance_1D(Element x, Element y);

double test_E_1D(Element x)
{
  double val = x.D1;
  double u = val - 1 ;
  return exp(-u * u)*sin(8*val);
/* return sin(sin(val*val) - cos(val)) + cos(sin(val) + sin(val)*sin(val));*/
/* return 1.0/(square(x-1.0) + 1.0)*sin(8.0*x); */
/* return 1.0/(square(x-1.0) + 1.0)*sin(8.0*x) + x*x/100.0; */
}

/* takes a step for the test function; max distance: step_size.
 * the new point is put in x_p and returned.
 */
void test_step_1D(const gsl_rng * r, Element *x_p, gsl_siman_params_t p)
{
  double old_x = x_p->D1;
  double new_x;

  new_x = gsl_rng_uniform(r);
  new_x = new_x*2*p.step_size;
  new_x = new_x - p.step_size + old_x;

  x_p->D1 = new_x;
}

/* simple routine to print out a position value */
void print_pos_1D(Element x)
{
  printf("%12g", x.D1);
}

void debug_pos_1D(void * x)
{
  printf("%12g", ((Element *)x)->D1);
}

/* a metric for the 2D space */
double distance_1D(Element x, Element y)
{
  return fabs(y.D1 - x.D1);
}

double test_E_2D(Element x);
void test_step_2D(const gsl_rng * r, Element *x_p, gsl_siman_params_t p);
void print_pos_2D(Element x);
double distance_2D(Element x, Element y);

/* a 2-D function to be minimized */
double test_E_2D(Element x)
{
  double old_x = x.D2[0], old_y = x.D2[1];
  double u = old_x-1 ;
  double v = old_y - 0.8 ;
  return exp(-u * u - v * v)*sin(8*old_x + 8 * old_y);
}

/* takes a step for the test function; max distance: step_size.  the
   new point is put in x_p and returned. */
void test_step_2D(const gsl_rng * r, Element *x_p, gsl_siman_params_t p)
{
  double old_x =  x_p->D2[0], old_y = x_p->D2[1], new_x, new_y;

  new_x = gsl_rng_uniform(r);
  new_x = new_x*2*p.step_size;
  new_x = new_x - p.step_size + old_x;

  new_y = gsl_rng_uniform(r);
  new_y = new_y*2*p.step_size;
  new_y = new_y - p.step_size + old_y;

  x_p->D2[0] = new_x;
  x_p->D2[1] = new_y;
}

/* simple routine to print out a position value */
void print_pos_2D(Element x)
{
  printf("%g::%g", x.D2[0], x.D2[1]);
}

/* a metric for the 2D space */
double distance_2D(Element x, Element y)
{
  double u = y.D2[0]-x.D2[0] ;
  double v = y.D2[1]-x.D2[1] ;
  return sqrt(u*u + v*v);
}

/**********************************************/
/************ 3-dimensional search ************/
/**********************************************/

double test_E_3D(Element x);
void test_step_3D(const gsl_rng * r, Element *x_p, gsl_siman_params_t p);
void print_pos_3D(Element x);
double distance_3D(Element x, Element y);

/* a 3-D function to be minimized */
double test_E_3D(Element x)
{
  return exp(-square(x.D3[0]-1) - square(x.D3[1] - 0.8)
    - square(x.D3[2] - 0.8)) * sin(8*x.D3[0] + 8*x.D3[1] + 8*x.D3[2])
      + (square(x.D3[0]) + square(x.D3[1]) + square(x.D3[2]))/10000.0;
}

/* takes a step for the test function; max distance: step_size.
 * the new point is put in x_p and returned.
 */
void test_step_3D(const gsl_rng * r, Element *x_p, gsl_siman_params_t p)
{
  double old_x = x_p->D3[0], old_y = x_p->D3[1], old_z = x_p->D3[2];
  double new_x, new_y, new_z;

  new_x = gsl_rng_uniform(r);
  new_x = new_x*2*p.step_size;
  new_x = new_x - p.step_size + old_x;

  new_y = gsl_rng_uniform(r);
  new_y = new_y*2*p.step_size;
  new_y = new_y - p.step_size + old_y;

  new_z = gsl_rng_uniform(r);
  new_z = new_z*2*p.step_size;
  new_z = new_z - p.step_size + old_z;

  x_p->D3[0] = new_x;
  x_p->D3[1] = new_y;
  x_p->D3[2] = new_z;
}

/* simple routine to print out a position value */
void print_pos_3D(Element x)
{
  printf("%g::%g::%g", x.D3[0], x.D3[1], x.D3[2]);
}

/* a metric for the 2D space */
double distance_3D(Element x, Element y)
{
  return sqrt(square(y.D3[0]-x.D3[0]) + square(y.D3[1]-x.D3[1])
	      + square(y.D3[2]-x.D3[2]));
}

double E1(void *xp);
double M1(void *xp, void *yp);
void S1(const gsl_rng * r, void *xp, double step_size);
void P1(void *xp);

/* now some functions to test in one dimension */
double E1(void *xp)
{
  double x = * ((double *) xp);

/*   return exp(-square(x-1))*sin(8*x); */
  return exp(-square(x-1))*sin(8*x) - exp(-square(x-1000))*0.89;
}

double M1(void *xp, void *yp)
{
  double x = *((double *) xp);
  double y = *((double *) yp);

  return fabs(x - y);
}

void S1(const gsl_rng * r, void *xp, double step_size)
{
  double old_x = *((double *) xp);
  double new_x;

  new_x = gsl_rng_uniform(r)*2*step_size - step_size + old_x;
/*   new_x = new_x*2*step_size; */
/*   new_x = new_x - step_size + old_x; */

  /* printf("test step from old_x = %g to new_x = %g\n", old_x, new_x); */

  memcpy(xp, &new_x, sizeof(new_x));
}

void P1(void *xp)
{
  printf(" %12g ", *((double *) xp));
}

int main(void)
{
  double x_min = 1.36312999455315182 ;
  double x ;

  gsl_rng * r = gsl_rng_alloc (gsl_rng_env_setup()) ;

  gsl_ieee_env_setup ();

  /* The function tested here has multiple mimima.
     The global minimum is at    x = 1.36312999, (f = -0.87287)
     There is a local minimum at x = 0.60146196, (f = -0.84893) */

  gsl_ieee_env_setup ();

  x = -10.0 ;
  gsl_siman_solve(r, &x, E1, S1, M1, NULL, NULL, NULL, NULL,
		  sizeof(double), params);
  gsl_test_rel(x, x_min, 1e-3, "f(x)= exp(-(x-1)^2) sin(8x), x0=-10") ;

  x = +10.0 ;
  gsl_siman_solve(r, &x, E1, S1, M1, NULL, NULL, NULL, NULL,
		  sizeof(double), params);
  gsl_test_rel(x, x_min, 1e-3, "f(x)= exp(-(x-1)^2) sin(8x), x0=10") ;

  /* Start at the false minimum */

  x = +0.6 ;
  gsl_siman_solve(r, &x, E1, S1, M1, NULL, NULL, NULL, NULL,
		  sizeof(double), params);
  gsl_test_rel(x, x_min, 1e-3, "f(x)= exp(-(x-1)^2) sin(8x), x0=0.6") ;

  x = +0.5 ;
  gsl_siman_solve(r, &x, E1, S1, M1, NULL, NULL, NULL, NULL,
		  sizeof(double), params);
  gsl_test_rel(x, x_min, 1e-3, "f(x)= exp(-(x-1)^2) sin(8x), x0=0.5") ;

  x = +0.4 ;
  gsl_siman_solve(r, &x, E1, S1, M1, NULL, NULL, NULL, NULL,
		  sizeof(double), params);
  gsl_test_rel(x, x_min, 1e-3, "f(x)= exp(-(x-1)^2) sin(8x), x0=0.4") ;

  exit (gsl_test_summary ());

#ifdef JUNK
  x0.D1 = 12.0;
  printf("#one dimensional problem, x0 = %f\n", x0.D1);
  gsl_siman_Usolve(r, &x0, test_E_1D, test_step_1D, distance_1D,
		   print_pos_1D, params);

  x0.D2[0] = 12.0;
  x0.D2[1] = 5.5;
  printf("#two dimensional problem, (x0,y0) = (%f,%f)\n",
	 x0.D2[0], x0.D2[1]);
  gsl_siman_Usolve(r, &x0, test_E_2D, test_step_2D, distance_2D,
		   print_pos_2D, params);

  x0.D3[0] = 12.2;
  x0.D3[1] = 5.5;
  x0.D3[2] = -15.5;
  printf("#three dimensional problem, (x0,y0,z0) = (%f,%f,%f)\n",
	 x0.D3[0], x0.D3[1], x0.D3[2]);
  gsl_siman_Usolve(r, &x0, test_E_3D, test_step_3D, distance_3D,
		   print_pos_3D, params);

  x0.D2[0] = 12.2;
  x0.D2[1] = 5.5;

  gsl_siman_solve(r, &x0, test_E_2D, test_step_2D, distance_2D, print_pos_2D, params);

  x0.D3[0] = 12.2;
  x0.D3[1] = 5.5;
  x0.D3[2] = -15.5;

  gsl_siman_solve(r, &x0, test_E_3D, test_step_3D, distance_3D, print_pos_3D, params);

  return 0;
#endif
}


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

#ifndef UTILS_H_
#define UTILS_H_

#include <gsl/gsl_permutation.h>
#include <gsl/gsl_randist.h>

#include <vector>

using namespace std;

namespace hlda {

// This class provides functionality for summing values,
// for reading data from files and also provides
// an interface to gsl specific methods.
// This class is not thread-safe (see private static gsl_rng* RANDNUMGEN).
class Utils {
 public:
  // Sum up the values in a vector.
  static double Sum(const vector<double>& v);

  // An approximation of the logarithm of a sum.
  // Given log(a) denoted with log_a and log(b) denoted with log_b
  // it approximates log(a + b) as:
  // if log_a < log_b then log(a + b) = log_b + log(1 + exp(log_a - log_b))
  // otherwise log(a + b) = log_a + log(1 + exp(log_b - log_a)).
  static double LogSum(double log_a, double log_b);

  // Sample log probabilities.
  // The log_pr vector keeps for each level
  // the current log probability of the word + the current
  // log probability of the level in the tree.
  // These values are used to sample the new level.
  // The vector should contain at least one element.
  static int SampleFromLogPr(const vector<double>& log_pr);

  // Shuffle the values in a gsl_permutation.
  static void Shuffle(gsl_permutation* permutation, int size);

  // Initialize the random number generator.
  // rng_seed is the random number generator seed.
  static void InitRandomNumberGen(long rng_seed);

  // Return a gsl Gaussian random variate with mean and stdev as parameters
  static double RandGauss(double mean, double stdev);

  // Return a random number using the gsl random number generator.
  static double RandNo();

 private:
  static gsl_rng* RANDNUMGEN;
};

}  // namespace hlda

#endif  // UTILS_H_

ext: h


text:
/**
 * @file constants.h
 * @brief Public header file defining constants
 */

#ifndef CONSTANTS_H
#define CONSTANTS_H

#include <petsc.h>
#define epsilon 1E-8 /**< epsilon value for finite differences */
#define MAXLINE 1000 /**< Max. number of characters in a line */
#define ISOLATED_BUS 4 /**< Isolated bus */
#define REF_BUS 3 /**< Reference bus (swing bus) */
#define PV_BUS 2 /**< PV (voltage-controlled) bus */
#define PQ_BUS 1 /**< PQ bus */
#define NGEN_AT_BUS_MAX 15 /**< Maximum number of generators allowed at a bus */
#define NLOAD_AT_BUS_MAX 10 /**< Maximum number of loads allowed at a bus */
#define NPHASE 1 /**< Per-phase analysis */

#endif

max_stars_repo_name: pflow-team/PFLOW

max_issues_repo_name: pflow-team/PFLOW

max_forks_repo_name: pflow-team/PFLOW


text:
#ifndef OPENMC_TALLIES_FILTER_CELL_INSTANCE_H
#define OPENMC_TALLIES_FILTER_CELL_INSTANCE_H

#include <cstdint>
#include <unordered_map>

#include <gsl/gsl>

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

namespace openmc {

//==============================================================================
//! Specifies cell instances that tally events reside in.
//==============================================================================

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

  CellInstanceFilter() = default;
  CellInstanceFilter(gsl::span<CellInstance> instances);
  ~CellInstanceFilter() = default;

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

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

  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

  const vector<CellInstance>& cell_instances() const { return cell_instances_; }

  const std::unordered_set<int32_t>& cells() const { return cells_; }

  void set_cell_instances(gsl::span<CellInstance> instances);

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

  //! The indices of the cells binned by this filter.
  vector<CellInstance> cell_instances_;

  //! The set of cells used in this filter
  std::unordered_set<int32_t> cells_;

  //! A map from cell/instance indices to filter bin indices.
  std::unordered_map<CellInstance, gsl::index, CellInstanceHash> map_;

  //! Indicates if filter uses only material-filled cells
  bool material_cells_only_;
};

} // namespace openmc

#endif // OPENMC_TALLIES_FILTER_CELL_INSTANCE_H


text:
#include <stdio.h>
#include <math.h>
#include <gbpLib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_min.h>
#include <gbpInterpolate.h>

double interpolate_maximum_function(double x, void *params);
double interpolate_maximum_function(double x, void *params) {
    return (-interpolate((interp_info *)params, x));
}

void interpolate_maximum(interp_info *interp,
                         double       x_lo_in,
                         double       x_guess_in,
                         double       x_hi_in,
                         double       threshold,
                         double *     x_maxima,
                         double *     y_maxima) {
    double                         x_lo;
    double                         x_hi;
    double                         x_guess;
    double                         r_val;
    const gsl_min_fminimizer_type *T;
    gsl_min_fminimizer *           s;
    gsl_function                   F;
    int                            status;
    int                            iter = 0;
    int                            max_iter;

    x_lo       = x_lo_in;
    x_guess    = x_guess_in;
    x_hi       = x_hi_in;
    max_iter   = GBP_MAX(100, interp->n);
    F.function = interpolate_maximum_function;
    F.params   = (void *)interp;
    T          = gsl_min_fminimizer_brent;
    s          = gsl_min_fminimizer_alloc(T);
    gsl_min_fminimizer_set(s, &F, x_guess, x_lo, x_hi);

    do {
        iter++;
        status  = gsl_min_fminimizer_iterate(s);
        x_lo    = gsl_min_fminimizer_x_lower(s);
        x_guess = gsl_min_fminimizer_x_minimum(s);
        x_hi    = gsl_min_fminimizer_x_upper(s);
        status  = gsl_min_test_interval(x_lo, x_hi, 0., threshold);
    } while(status == GSL_CONTINUE && iter < max_iter);
    gsl_min_fminimizer_free(s);
    (*x_maxima) = x_guess;
    (*y_maxima) = interpolate(interp, x_guess);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#define ORIGIN_DATE       "2012.01.05"
#define ORIGIN_VERSION    "2.1"

#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <pthread.h>
#include <assert.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "common.h"
#include "ivector.h"
#include "species.h"
#include "specieslist.h"
#include "graph.h"
#include "utils.h"

#define MODEL_BDM_NEUTRAL      0
#define MODEL_BDM_SELECTION    1

// Parameters of the simulations.
typedef struct
{
    int m;             // Type of model.
    int communities;   // Number of communities.
    int j_per_c;       // Number of individuals per community.
    int k_gen;         // Number of generations (in thousands).
    int init_species;  // Initial number of species.
    double mu;         // Mutation rate.
    double omega;      // Weight of the links between communities.
    double s;          // Selection coefficient.
    double r;          // Radius (for random geometric graphs).
    double w;          // Width (for rectangle random geometric graphs).
    char *ofilename;   // Name of the output files.
    char *shape;       // Shape of the metacommunity.
}
Params;

// Prototype for the function used by the threads
void *sim(void *parameters);
// The function used to setup the cumulative jagged array from the graph.
double **setup_cumulative_list(const graph *g, double omega);

/////////////////////////////////////////////////////////////
// Main                                                    //
/////////////////////////////////////////////////////////////
int main(int argc, const char *argv[])
{
    // Set default values;
    Params p;
    p.m = MODEL_BDM_SELECTION;
    p.communities = 10;
    p.j_per_c = 10000;
    p.init_species = 20;
    p.k_gen = 100;
    p.mu = 1e-4;
    p.omega = 5e-4;
    p.s = 0.15;
    p.r = 0.25;
    p.w = 0.25;
    p.ofilename = (char*)malloc(50);
    p.shape = (char*)malloc(20);

    // Number of simulations;
    int n_threads = 1;

    // Options;
    if (argc > 1 && argv[1][0] == '-' && argv[1][1] == '-')
    {
        // --help
        if (argv[1][2] == 'h')
        {
            printf("Usage: ./origin [options]\n");
            printf("Example: ./origin -x=8 -shape=circle -model=1 -s=0.10\n");
            printf("General options:\n");
            printf("  --help          How you got here...\n");
            printf("  --version       Display version.\n");
            printf("  --ref           Display reference.\n");
            printf("Simulation parameters:\n");
            printf("  -x\n");
            printf("    description:  Number of simulations to run. Each simulation\n");
            printf("                  will use a distinct POSIX thread. \n");
            printf("    values:       Any unsigned integer.\n");
            printf("    default:      1\n");
            printf("  -model\n");
            printf("    description:  Set the model used.\n");
            printf("    values:       0, 1.\n");
            printf("    details:      0 = Neutral BDM speciation.\n");
            printf("                  1 = BDM speciation with selection.\n");
            printf("    default:      1\n");
            printf("  -shape\n");
            printf("    description:  Set the shape of the metacommunity.\n");
            printf("    values:       circle, complete, random, rectangle, star.\n");
            printf("    default:      random\n");
            printf("  -r\n");
            printf("    description:  Threshold radius for random geometric\n");
            printf("                  graphs.\n");
            printf("    values:       Any double greater than 0.\n");
            printf("    default:      0.25\n");
            printf("  -w\n");
            printf("    description:  The width of a rectangle random geometric\n");
            printf("                  graph.\n");
            printf("    values:       0 < w < 1\n");
            printf("    default:      0.25\n");
            printf("  -g\n");
            printf("    description:  Number of generations in thousands.\n");
            printf("    values:       Any positive integer.\n");
            printf("    default:      100 (i.e.: 100 000 generations)\n");
            printf("  -c\n");
            printf("    description:  Number of local communities.\n");
            printf("    values:       Any integer greater than 0.\n");
            printf("    default:      10\n");
            printf("  -jpc\n");
            printf("    description:  Number of individuals per communities.\n");
            printf("    values:       Any integer greater than 0.\n");
            printf("    default:      10000\n");
            printf("  -sp\n");
            printf("    description:  Initial number of species. Must be a factor\n");
            printf("                  of the number of individuals/community.\n");
            printf("    values:       Any positive integer greater than 0.\n");
            printf("    default:      20\n");
            printf("  -s\n");
            printf("    description:  Set the selection coefficient (model 1 only).\n");
            printf("    values:       Any double in the (-1.0, 1.0) interval.\n");
            printf("    default:      0.15\n");
            printf("  -mu\n");
            printf("    description:  Set the mutation rate.\n");
            printf("    values:       Any positive double.\n");
            printf("    default:      1e-4\n");
            printf("  -omega\n");
            printf("    description:  The weight of the proper edges.\n");
            printf("    values:       Any nonnegative double.\n");
            printf("    default:      5e-4\n");
            printf("  -o\n");
            printf("    description:  Name of the output files.\n");
            printf("    values:       Any string.\n");
            return EXIT_SUCCESS;
        }
        else if (argv[1][2] == 'v') // --version
        {
            printf("origin v.%s (%s)\n", ORIGIN_VERSION, ORIGIN_DATE);
            printf("Copyright (c) 2010-2011 Philippe Desjardins-Proulx <philippe.d.proulx@gmail.com>\n");
            printf("GPLv2 license  (see LICENSE)\n");
            return EXIT_SUCCESS;
        }
        else if (argv[1][2] == 'r') // --ref
        {
            printf("P. Desjardins-Proul and D. Gravel. How likely is speciation in\n  neutral ecology? The American Naturalist 179(1).\n");
            return EXIT_SUCCESS;
        }
    } // end '--' options

    // Read options
    read_opt_i("g", argv, argc, &p.k_gen);
    read_opt_i("jpc", argv, argc, &p.j_per_c);
    read_opt_i("model", argv, argc, &p.m);
    read_opt_i("c", argv, argc, &p.communities);
    read_opt_i("sp", argv, argc, &p.init_species);
    read_opt_i("x", argv, argc, &n_threads);
    read_opt_d("mu", argv, argc, &p.mu);
    read_opt_d("omega", argv, argc, &p.omega);
    read_opt_d("r", argv, argc, &p.r);
    read_opt_d("w", argv, argc, &p.w);
    read_opt_d("s", argv, argc, &p.s);
    if (read_opt_s("o", argv, argc, p.ofilename) == false)
    {
        sprintf(p.ofilename, "");
    }
    if (read_opt_s("shape", argv, argc, p.shape) == false)
    {
        sprintf(p.shape, "random");
    }

    printf("<?xml version=\"1.0\"?>\n");
    printf("<origin_ssne>\n");
    printf("  <model>");
    if (p.m == MODEL_BDM_NEUTRAL)
    {
        printf("Neutral BDM speciation</model>\n");
    }
    else if (p.m == MODEL_BDM_SELECTION)
    {
        printf("BDM speciation with selection</model>\n");
    }
    printf("  <n>%d</n>\n", n_threads);
    printf("  <shape_metacom>%s</shape_metacom>\n", p.shape);
    printf("  <metacom_size>%d</metacom_size>\n", p.j_per_c * p.communities);
    printf("  <k_gen>%d</k_gen>\n", p.k_gen);
    printf("  <num_comm>%d</num_comm>\n", p.communities);
    printf("  <individuals_per_comm>%d</individuals_per_comm>\n", p.j_per_c);
    printf("  <initial_num_species>%d</initial_num_species>\n", p.init_species);
    printf("  <mutation_rate>%.2e</mutation_rate>\n", p.mu);
    printf("  <omega>%.2e</omega>\n", p.omega);
    if (p.m == MODEL_BDM_SELECTION)
    {
        printf("  <selection>%.2e</selection>\n", p.s);
    }
    if (p.shape[0] == 'r')
    {
        printf("  <radius>%.4f</radius>\n", p.r);
    }
    if (p.shape[0] == 'r' && p.shape[1] == 'e')
    {
        printf("  <width>%.4f</width>\n", p.w);
    }
    printf("  <filename>%s</filename>\n", p.ofilename);

    // The threads:
    pthread_t threads[n_threads];

    const time_t start = time(NULL);

    // Create the threads:
    for (int i = 0; i < n_threads; ++i)
    {
        pthread_create(&threads[i], NULL, sim, (void*)&p);
    }

    // Wait for the threads to end:
    for (int i = 0; i < n_threads; ++i)
    {
        pthread_join(threads[i], NULL);
    }

    const time_t end_t = time(NULL);
    printf("  <seconds>%lu</seconds>\n", (unsigned long)(end_t - start));
    printf("  <time>%s</time>\n", sec_to_string(end_t - start));
    printf("</origin_ssne>\n");
    return EXIT_SUCCESS; // yeppie !
}

///////////////////////////////////////////
// The simulation                        //
///////////////////////////////////////////
void *sim(void *parameters)
{
    const Params P = *((Params*)parameters);

    char *shape = P.shape;
    const int k_gen = P.k_gen;
    const int communities = P.communities;
    const int j_per_c = P.j_per_c;
    const int init_species = P.init_species;
    const int init_pop_size = j_per_c / init_species;
    const double omega = P.omega;
    const double mu = P.mu;
    const double s = P.s;
    const double radius = P.r;
    const double width = P.w;

    // GSL's Taus generator:
    gsl_rng *rng = gsl_rng_alloc(gsl_rng_taus2);
    // Initialize the GSL generator with /dev/urandom:
    const unsigned int seed = devurandom_get_uint();
    gsl_rng_set(rng, seed); // Seed with time
    printf("  <seed>%u</seed>\n", seed);
    // Used to name the output file:
    char *buffer = (char*)malloc(100);
    // Nme of the file:
    sprintf(buffer, "%s%u.xml", P.ofilename, seed);
    // Open the output file:
    FILE *restrict out = fopen(buffer, "w");
    // Store the total num. of species/1000 generations:
    int *restrict total_species = (int*)malloc(k_gen * sizeof(int));
    // Number of speciation events/1000 generations:
    int *restrict speciation_events = (int*)malloc(k_gen * sizeof(int));
    // Number of extinctions/1000 generations:
    int *restrict extinction_events = (int*)malloc(k_gen * sizeof(int));
    // Number of speciation events/vertex:
    int *restrict speciation_per_c = (int*)malloc(communities * sizeof(int));
    // Number of local extinction events/vertex:
    int *restrict extinction_per_c = (int*)malloc(communities * sizeof(int));
    // Store the lifespan of the extinct species:
    ivector lifespan;
    ivector_init0(&lifespan);
    // Store the population size at speciation event:
    ivector pop_size;
    ivector_init0(&pop_size);
    // (x, y) coordinates for the spatial graph:
    double *restrict x = (double*)malloc(communities * sizeof(double));
    double *restrict y = (double*)malloc(communities * sizeof(double));
    for (int i = 0; i < communities; ++i)
    {
        speciation_per_c[i] = 0;
        extinction_per_c[i] = 0;
    }
    // Initialize an empty list of species:
    species_list *restrict list = species_list_init();
    // Initialize the metacommunity and fill them with the initial species evenly:
    for (int i = 0; i < init_species; ++i)
    {
        // Intialize the species and add it to the list:
        species_list_add(list, species_init(communities, 0, 3));
    }
    // To iterate the list;
    slnode *it = list->head;
    // Fill the communities:
    while (it != NULL)
    {
        for (int i = 0; i < communities; ++i)
        {
            it->sp->n[i] = init_pop_size;
            it->sp->genotypes[0][i] = init_pop_size;
        }
        it = it->next;
    }

    // To iterate the list;
    const int remainder = j_per_c - (init_species * init_pop_size);
    for (int i = 0; i < communities; ++i)
    {
        it = list->head;
        for (int j = 0; j < remainder; ++j, it = it->next)
        {
            ++(it->sp->n[i]);
            ++(it->sp->genotypes[0][i]);
        }
    }

    // Test (will be removed for v2.0):
    int sum = 0;
    it = list->head;
    while (it != NULL)
    {
        sum += species_total(it->sp);
        it = it->next;
    }
    assert(sum == j_per_c * communities);

    // Create the metacommunity;
    graph g;
    switch(shape[0])
    {
    case 's':
        shape = "star";
        graph_get_star(&g, communities);
        break;
    case 'c':
        if (shape[1] == 'o')
        {
            shape = "complete";
            graph_get_complete(&g, communities);
            break;
        }
        else
        {
            shape = "circle";
            graph_get_circle(&g, communities);
            break;
        }
    case 'r':
        if (shape[1] == 'e')
        {
            shape = "rectangle";
            graph_get_rec_crgg(&g, communities, width, radius, x, y, rng);
            break;
        }
        else
        {
            shape = "random";
            graph_get_crgg(&g, communities, radius, x, y, rng);
            break;
        }
    default:
        shape = "random";
        graph_get_crgg(&g, communities, radius, x, y, rng);
    }
    // Setup the cumulative jagged array for migration:
    double **cumul = setup_cumulative_list(&g, omega);

    fprintf(out, "<?xml version=\"1.0\"?>\n");
    fprintf(out, "<simulation>\n");
    fprintf(out, "  <model>");
    if (P.m == MODEL_BDM_NEUTRAL)
    {
        fprintf(out, "Neutral BDM speciation</model>\n");
    }
    else if (P.m == MODEL_BDM_SELECTION)
    {
        fprintf(out, "BDM speciation with selection</model>\n");
    }
    fprintf(out, "  <seed>%u</seed>\n", seed);
    fprintf(out, "  <shape_metacom>%s</shape_metacom>\n", shape);
    fprintf(out, "  <metacom_size>%d</metacom_size>\n", j_per_c * communities);
    fprintf(out, "  <k_gen>%d</k_gen>\n", k_gen);
    fprintf(out, "  <num_comm>%d</num_comm>\n", communities);
    fprintf(out, "  <individuals_per_comm>%d</individuals_per_comm>\n", j_per_c);
    fprintf(out, "  <initial_num_species>%d</initial_num_species>\n", init_species);
    fprintf(out, "  <mutation_rate>%.2e</mutation_rate>\n", mu);
    fprintf(out, "  <omega>%.2e</omega>\n", omega);
    if (P.m == MODEL_BDM_SELECTION)
    {
        fprintf(out, "  <selection>%.2e</selection>\n", s);
    }
    if (shape[0] == 'r')
    {
        fprintf(out, "  <radius>%.4f</radius>\n", radius);
    }
    if (shape[0] == 'r' && shape[1] == 'e')
    {
        fprintf(out, "  <width>%.4f</width>\n", width);
    }
    // To select the species and genotypes to pick and replace:
    slnode *s0 = list->head; // species0
    slnode *s1 = list->head; // species1
    int g0 = 0;
    int g1 = 0;
    int v1 = 0; // Vertex of the individual 1

    /////////////////////////////////////////////
    // Groups of 1 000 generations             //
    /////////////////////////////////////////////
    for (int k = 0; k < k_gen; ++k)
    {
        extinction_events[k] = 0;
        speciation_events[k] = 0;

        /////////////////////////////////////////////
        // 1 000 generations                       //
        /////////////////////////////////////////////
        for (int gen = 0; gen < 1000; ++gen)
        {
            const int current_date = (k * 1000) + gen;
            /////////////////////////////////////////////
            // A single generation                     //
            /////////////////////////////////////////////
            for (int t = 0; t < j_per_c; ++t)
            {
                /////////////////////////////////////////////
                // A single time step (for each community) //
                /////////////////////////////////////////////
                for (int c = 0; c < communities; ++c)
                {
                    // Select the species and genotype of the individual to be replaced
                    int position = (int)(gsl_rng_uniform(rng) * j_per_c);
                    s0 = list->head;
                    int index = s0->sp->n[c];
                    while (index <= position)
                    {
                        s0 = s0->next;
                        index += s0->sp->n[c];
                    }
                    position = (int)(gsl_rng_uniform(rng) * s0->sp->n[c]);
                    if (position < s0->sp->genotypes[0][c])
                    {
                        g0 = 0;
                    }
                    else if (position < (s0->sp->genotypes[0][c] + s0->sp->genotypes[1][c]))
                    {
                        g0 = 1;
                    }
                    else
                    {
                        g0 = 2;
                    }
                    // Choose the vertex for the individual
                    const double r_v1 = gsl_rng_uniform(rng);
                    v1 = 0;
                    while (r_v1 > cumul[c][v1])
                    {
                        ++v1;
                    }
                    v1 = g.adj_list[c][v1];
                    // species of the new individual
                    position = (int)(gsl_rng_uniform(rng) * j_per_c);
                    s1 = list->head;
                    index = s1->sp->n[v1];
                    while (index <= position)
                    {
                        s1 = s1->next;
                        index += s1->sp->n[v1];
                    }
                    if (v1 == c) // local remplacement
                    {
                        const double r = gsl_rng_uniform(rng);
                        const int aa = s1->sp->genotypes[0][v1];
                        const int Ab = s1->sp->genotypes[1][v1];
                        const int AB = s1->sp->genotypes[2][v1];

                        // The total fitness of the population 'W':
                        const double w = aa + Ab * (1.0 + s) + AB * (1.0 + s) * (1.0 + s);

                        if (r < aa / w)
                        {
                            g1 = gsl_rng_uniform(rng) < mu ? 1 : 0;
                        }
                        else
                        {
                            if (AB == 0 || r < (aa + Ab * (1.0 + s)) / w)
                            {
                                g1 = gsl_rng_uniform(rng) < mu ? 2 : 1;
                            }
                            else
                            {
                                g1 = 2;
                            }
                        }
                    }
                    else
                    { // Migration event
                        g1 = 0;
                    }
                    // Apply the changes
                    s0->sp->n[c]--;
                    s0->sp->genotypes[g0][c]--;
                    s1->sp->n[c]++;
                    s1->sp->genotypes[g1][c]++;

                    ////////////////////////////////////////////
                    // Check for local extinction             //
                    ////////////////////////////////////////////
                    if (s0->sp->n[c] == 0)
                    {
                        extinction_per_c[c]++;
                    }
                    ////////////////////////////////////////////
                    // Check for speciation                   //
                    ////////////////////////////////////////////
                    else if (s0->sp->genotypes[2][c] > 0 && s0->sp->genotypes[0][c] == 0 && s0->sp->genotypes[1][c] == 0)
                    {
                        species_list_add(list, species_init(communities, current_date, 3)); // Add the new species

                        const int pop = s0->sp->n[c];
                        list->tail->sp->n[c] = pop;
                        list->tail->sp->genotypes[0][c] = pop;
                        s0->sp->n[c] = 0;
                        s0->sp->genotypes[2][c] = 0;

                        // To keep info on patterns of speciation...
                        ivector_add(&pop_size, pop);
                        ++speciation_events[k];
                        ++speciation_per_c[c];
                    }

                } // End 'c'

            } // End 't'

            // Remove extinct species from the list and store the number of extinctions.
            extinction_events[k] += species_list_rmv_extinct2(list, &lifespan, current_date);

        } // End 'g'

        total_species[k] = list->size;

    } // End 'k'

    //////////////////////////////////////////////////
    // PRINT THE FINAL RESULTS                      //
    //////////////////////////////////////////////////
    fprintf(out, "  <global>\n");
    fprintf(out, "    <proper_edges>%d</proper_edges>\n", graph_edges(&g));
    fprintf(out, "    <links_per_c>%.4f</links_per_c>\n", (double)graph_edges(&g) / communities);

    fprintf(out, "    <avr_lifespan>%.4f</avr_lifespan>\n", imean(lifespan.array, lifespan.size));
    fprintf(out, "    <median_lifespan>%.4f</median_lifespan>\n", imedian(lifespan.array, lifespan.size));

    fprintf(out, "    <avr_pop_size_speciation>%.4f</avr_pop_size_speciation>\n", imean(pop_size.array, pop_size.size));
    fprintf(out, "    <median_pop_size_speciation>%.4f</median_pop_size_speciation>\n", imedian(pop_size.array, pop_size.size));

    fprintf(out, "    <speciation_per_k_gen>");
    int i = 0;
    for (; i < k_gen - 1; ++i)
    {
        fprintf(out, "%d ", speciation_events[i]);
    }
    fprintf(out, "%d</speciation_per_k_gen>\n", speciation_events[i]);

    fprintf(out, "    <extinctions_per_k_gen>");
    for (i = 0; i < k_gen - 1; ++i)
    {
        fprintf(out, "%d ", extinction_events[i]);
    }
    fprintf(out, "%d</extinctions_per_k_gen>\n", extinction_events[i]);

    fprintf(out, "    <extant_species_per_k_gen>");
    for (i = 0; i < k_gen - 1; ++i)
    {
        fprintf(out, "%d ", total_species[i]);
    }
    fprintf(out, "%d</extant_species_per_k_gen>\n", total_species[i]);

    // Print global distribution
    fprintf(out, "    <species_distribution>");
    ivector species_distribution;
    ivector_init1(&species_distribution, 128);
    it = list->head;
    while (it != NULL)
    {
        ivector_add(&species_distribution, species_total(it->sp));
        it = it->next;
    }
    ivector_sort_asc(&species_distribution);
    ivector_print(&species_distribution, out);
    fprintf(out, "</species_distribution>\n");

    double *octaves;
    int oct_num = biodiversity_octaves(species_distribution.array, species_distribution.size, &octaves);
    fprintf(out, "    <octaves>");
    for (i = 0; i < oct_num; ++i)
    {
        fprintf(out, "%.2f ", octaves[i]);
    }
    fprintf(out, "%.2f</octaves>\n", octaves[i]);
    fprintf(out, "  </global>\n");

    // Print info on all vertices
    double *restrict ric_per_c = (double*)malloc(communities * sizeof(double));
    for (int c = 0; c < communities; ++c)
    {
        fprintf(out, "  <vertex>\n");
        fprintf(out, "    <id>%d</id>\n", c);
        if (shape[0] == 'r')
        {
            fprintf(out, "    <xcoor>%.4f</xcoor>\n", x[c]);
            fprintf(out, "    <ycoor>%.4f</ycoor>\n", y[c]);
        }
        fprintf(out, "    <degree>%d</degree>\n", g.num_e[c] + 1);
        fprintf(out, "    <speciation_events>%d</speciation_events>\n", speciation_per_c[c]);
        fprintf(out, "    <extinction_events>%d</extinction_events>\n", extinction_per_c[c]);

        int vertex_richess = 0;
        ivector_rmvall(&species_distribution);
        it = list->head;
        while (it != NULL)
        {
            ivector_add(&species_distribution, it->sp->n[c]);
            it = it->next;
        }
        // Sort the species distribution and remove the 0s
        ivector_sort_asc(&species_distribution);
        ivector_trim_small(&species_distribution, 1);

        ric_per_c[c] = (double)species_distribution.size;
        fprintf(out, "    <species_richess>%d</species_richess>\n", species_distribution.size);
        fprintf(out, "    <species_distribution>");
        ivector_print(&species_distribution, out);
        fprintf(out, "</species_distribution>\n");

        // Print octaves
        free(octaves);
        oct_num = biodiversity_octaves(species_distribution.array, species_distribution.size, &octaves);
        fprintf(out, "    <octaves>");
        for (i = 0; i < oct_num - 1; ++i)
        {
            fprintf(out, "%.2f ", octaves[i]);
        }
        fprintf(out, "%.2f</octaves>\n", octaves[i]);
        fprintf(out, "  </vertex>\n");
    }

    fprintf(out, "</simulation>\n");

    // GraphML output:
    sprintf(buffer, "%s%u.graphml", P.ofilename, seed);
    FILE *outgml = fopen(buffer, "w");
    graph_graphml(&g, outgml, seed);

    // Print to SVG files.
    sprintf(buffer, "%s%u.svg", P.ofilename, seed);
    FILE *outsvg = fopen(buffer, "w");
    graph_svg(&g, x, y, 400.0, 20.0, outsvg);

    sprintf(buffer, "%s%u-speciation.svg", P.ofilename, seed);
    FILE *outsvgspe = fopen(buffer, "w");
    double *spe_per_c = (double*)malloc(communities * sizeof(double));
    for (int c = 0; c < communities; ++c)
    {
        spe_per_c[c] = (double)speciation_per_c[c];
    }
    scale_0_1(spe_per_c, communities);
    graph_svg_abun(&g, x, y, 400.0, 20.0, spe_per_c, 2, outsvgspe);

    sprintf(buffer, "%s%u-richness.svg", P.ofilename, seed);
    FILE *outsvgric = fopen(buffer, "w");
    scale_0_1(ric_per_c, communities);
    graph_svg_abun(&g, x, y, 400.0, 20.0, ric_per_c, 1, outsvgric);

    //////////////////////////////////////////////////
    // EPILOGUE...                                  //
    //////////////////////////////////////////////////
    // Close files;
    fclose(out);
    fclose(outgml);
    fclose(outsvg);
    fclose(outsvgspe);
    fclose(outsvgric);
    // Free arrays;
    free(x);
    free(y);
    free(ric_per_c);
    free(spe_per_c);
    free(buffer);
    free(total_species);
    free(octaves);
    free(speciation_per_c);
    free(extinction_per_c);
    free(speciation_events);
    free(extinction_events);
    // Free structs;
    species_list_free(list);
    ivector_free(&species_distribution);
    ivector_free(&lifespan);
    ivector_free(&pop_size);
    graph_free(&g);
    gsl_rng_free(rng);

    return NULL;
}

double **setup_cumulative_list(const graph *g, double omega)
{
    const int num_v = g->num_v;
    double **cumul = (double**)malloc(num_v * sizeof(double*));
    for (int i = 0; i < num_v; ++i)
    {
        cumul[i] = (double*)malloc(g->num_e[i] * sizeof(double));
    }
    for (int i = 0; i < num_v; ++i)
    {
        for (int j = 0; j < g->num_e[i]; ++j)
        {
            cumul[i][j] = g->w_list[i][j];
        }
    }
    for (int i = 0; i < num_v; ++i)
    {
        for (int j = 0; j < g->num_e[i]; ++j)
        {
            if (i != g->adj_list[i][j])
            {
                cumul[i][j] = omega;
            }
            else
            {
                cumul[i][j] = 1.0;
            }
        }
    }
    for (int i = 0; i < num_v; ++i)
    {
        double sum = 0.0;
        const int num_e = g->num_e[i];
        for (int j = 0; j < num_e; ++j)
        {
            sum += cumul[i][j];
        }
        for (int j = 0; j < num_e; ++j)
        {
            cumul[i][j] /= sum;
        }
    }
    for (int i = 0; i < num_v; ++i)
    {
        const int num_e = g->num_e[i];
        for (int j = 1; j < num_e - 1; ++j)
        {
            cumul[i][j] += cumul[i][j - 1];
        }
        cumul[i][num_e - 1] = 1.0;
    }
    /*
    graph_print(g, stdout);
    for (int i = 0; i < num_v; ++i)
    {
        const int num_e = g->num_e[i];
        for (int j = 0; j < num_e; ++j)
        {
            printf("%.4f   ", cumul[i][j]);
        }
        printf("\n");
    }
    */
    return cumul;
}

max_stars_repo_name: PhDP/Origin

max_issues_repo_name: PhDP/Origin

max_forks_repo_name: PhDP/Origin


text:
#ifndef MOVE_H
#define MOVE_H
#include <cmath>
#include "param.h"
#include "mpiutils.h"
#include <gsl/gsl_rng.h>

namespace weakarg
{
extern gsl_rng * rng;

/**
    @brief This virtual class represents a move of the MCMC
*/
class Move
{
protected:
    std::string name;///<The moves name
    std::string description;///<The moves description
    Param*param;///<Parameter set on which the move is performed
    double alpha;///<Weighting for this move (probability is alpha/sum_i(alpha_i) )
    int intDist; //FMA_CHANGES: Number of intermediate distns for AISRJ
    int numcalls;///<Number of attempts at the move
    int numaccept;///<Number of acceptances of the move
public:
    Move(Param*p,double a=1.0);///<Constructor specifying on which Parameter the move is to be performed
    virtual Move * clone()=0;
    virtual ~Move()=0;
    virtual int move()=0;///<Performs the move once
    virtual inline int getTopoCounts(){return(-1);};
    virtual inline int getTopoAcc(){return(-1);};
    virtual inline int move(vector<int> *samplespace){cout<<"Problem in Move!"<<endl;return(move());};///<Performs the move once
    inline void setParam(Param * param)
    {
        this->param=param;
    }///<Change the Parameter on which the move is to be performed
    inline std::string desc()
    {
        return description;
    }///<Returns a description of the move
    inline std::string getName()
    {
        return name;
    }///<Returns a description of the move
    inline void setName(std::string newname)
    {
        name=newname;
    }
    inline void setAlpha(double a)
    {
        alpha=a;
    }
    ;///<Sets alpha for the move
    inline double getAlpha()
    {
        return(alpha);
    }
    inline int getCounts()
    {
        return(numcalls);
    }
    inline int getAcc()
    {
        return(numaccept);
    }

    ;///<Returns alpha for the move
};

} // end namespace weakarg
#endif


text:
#pragma once

#include <gsl/gsl>

#include <string>
#include <string_view>
#include <memory>

#include "aas/aas_math.h"

namespace aas {

	/**
	 * Animations can either be advanced by time (automatic), distance moved, or object rotation.
	 */
	enum class SkelAnimDriver : uint8_t {
		Time = 0,
		Distance,
		Rotation
	};

	struct SkelAnimStream {
		uint16_t frames;
		int16_t variationId; // This might actually be a uint8_t
		float frameRate;
		float dps;
		// Offset to the start of the key-frame stream in relation to SkaAnimation start
		uint32_t dataOffset;
	};

	struct SkelAnimStreamData;

	struct SkelAnimEvent {
		int16_t frame;
		char type[48];
		char action[128];
	};

	struct SkelAnim {
		char name[64];
		SkelAnimDriver driveType;
		uint8_t loopable;
		uint16_t eventCount; // seems to be uint8_t in reality...
		uint32_t eventOffset;
		uint16_t streamCount;
		uint16_t unk;
		SkelAnimStream streams[10];

		const gsl::span<SkelAnimEvent> GetEventData() const {
			if (eventCount == 0) {
				return {};
			} else {
				auto data = (SkelAnimEvent*)(((char*)this) + eventOffset);
				return gsl::span(data, eventCount);
			}
		}

		const SkelAnimStreamData* GetStreamKeyframes(int streamIdx = 0) const {
			return (SkelAnimStreamData*)(((char*)this) + streams[streamIdx].dataOffset);
		}
	};

	struct SkelBoneState {
		DX::XMFLOAT4 scale;
		DX::XMFLOAT4 rotation;
		DX::XMFLOAT4 translation;

		static void Lerp(gsl::span<SkelBoneState> out,
			const gsl::span<SkelBoneState> from,
			const gsl::span<SkelBoneState> to,
			int count,
			float fraction);
	};

	struct SkelBone {
		int16_t flags;
		int16_t parentId;
		char name[40];
		int field_2C;
		float field_30;
		SkelBoneState initialState;
	};

	class Skeleton {
	public:
		Skeleton(std::string filename, std::vector<uint8_t> data);

		gsl::span<SkelBone> GetBones() const {
			return bones_;
		}

		gsl::span<SkelAnim> GetAnimations() const {
			return animations_;
		}

		/**
		 * Searches for a bone in this skeleton by name (case-insensitive).
		 * Returns nullptr if no bone is found.
		 */
		const SkelBone* FindBoneByName(std::string_view name) const;

		/**
		* Searches for the index of a bone in this skeleton by name (case-insensitive).
		* Returns -1 if no bone is found.
		*/
		int FindBoneIdxByName(std::string_view name) const;

		/**
		 * Searches for an animation in this skeleton by it's name (case insensitive)
		 * and returns the index of the animation in the skeleton.
		 * If no animation is found, returns -1.
		 */
		int FindAnimIdxByName(std::string_view name) const;

		const std::string &GetFilename() const {
			return filename_;
		}

	private:
		std::string filename_;
		std::vector<uint8_t> data_;

		// These are views into the data buffer
		gsl::span<SkelBone> bones_;
		gsl::span<SkelAnim> animations_;
	};

	std::unique_ptr<Skeleton> LoadSkeletonFile(std::string_view filename);

}

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
/* ============================================================ *
 * decomp_eb.h							*
 * Martin Kilbinger, Liping Fu 2008, 2009			*
 * ============================================================ */

#ifndef __DECOMP_EB_H
#define __DECOMP_EB_H

#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_combination.h>

#include "maths.h"
#include "io.h"
#include "errorlist.h"

/* Error codes */
#define mr_base		-1500
#define mr_range	-1 + mr_base
#define mr_poly		-2 + mr_base
#define mr_func		-3 + mr_base
#define mr_type		-4 + mr_base
#define mr_incompatible -5 + mr_base
#define mr_null         -6 + mr_base
#define mr_dimension    -7 + mr_base
#define mr_file         -8 + mr_base

/* Maximum COSEBI mode. The code is accurate up *
 * to Nmax_cosebi = 13			        */
#define NMAX_COSEBI     20

typedef enum {cheby, cheby2, legen, cosebi} poly_t;
#define spoly_t(i) ( \
  i==cheby  ? "cheby" : \
  i==cheby2 ? "cheby2" : \
  i==legen  ? "legen" : \
  i==cosebi  ? "cosebi" : \
  "")
#define Npoly_t 4

/*
typedef enum {cosebi} filter_t;
#define sfilter_t(i) ( \
  i==cosebu  ? "cosebi" : \
  "")
#define Nfilter_t 1
*/

typedef enum {comb_none, all_sc} mring_combinef_t;
#define smring_combinef_t(i) ( \
   i==comb_none    ? "comb_none" : \
   i==all_sc       ? "all_sc" : \
   "")
#define Nmring_combinef_t 2

typedef enum {inner, outer, fromZ} cov_mode_t;

/* === CFHTLS Wide 3rd data release, Fu&Kilbinger (2010) === */
#define N_FK10 6
#define eta_FK10_SN 0.02
#define eta_FK10_FoM_eta10 0.1
#define eta_FK10_FoM_eta50 0.02

double Cheby(double x, int n, error **err);
double Cheby2(double x, int n, error **err);
double Legen(double x, int n);
double C(double x, int n, poly_t poly, error **err);
double Tp(double x, const double *a, int N, poly_t poly, error **err);

double Fn0(double x, int n, poly_t poly, error **err);
void Fnnu(double x, int n,  poly_t poly, double Fn[], error **err);
double alpha_explicit(double x, int n, double R, poly_t poly, error **err);
double Tm(double x, const double *a, int N, poly_t poly, double R, error **err);

double RR_data(const double *xip, const double *xim, const double *th, const int Nxi,
	       double THETA_MIN, double THETA_MAX, const double *a, int N,
	       poly_t poly, int pm, error **err);

double cov_RR(const double *THETA_MIN, const double *THETA_MAX, const double *a, int N, poly_t poly,
	      const double *theta, const double *cov_xi, int Ntheta,
	      cov_mode_t cov_mode, int n, int m, double fac, error **err);
double cov_RR_diag_xi(const double *THETA_MIN, const double *THETA_MAX, const double *a, int N, poly_t poly,
		      const double *theta, const double *var_xi, int Ntheta, int islog,
		      error **err);

double chi2_RB_null(const double *RB, const double *covRB, int NRB);

double *read_zeros_norm_cosebi_auto_check(double Psimin, double Psimax, const char *path, error **err);
double *read_zeros_norm_cosebi(const char *rname, double *psimin, double *psimax, error **err);
double Tplog(double x, const double *r, int n, error **err);
double Tplog_c(double z, const double *c, int n, error **err);
double Tmlog(double x, const double *c, int n, error **err);
double an2(int n,  const double *c);
double an4(int n, const double *c);
double dnm(int n, const int m, const double *c);
double sum_combinations(int j, int n, const double *r, error **err);
void xipmEB(double theta, double THETA_MIN, double THETA_MAX, const double *c,
	    const double *E, const double *B, int N, double xi_pm_EB[4], error **err);

#endif

max_stars_repo_name: danielgruen/ccv

max_issues_repo_name: danielgruen/ccv

max_forks_repo_name: danielgruen/ccv


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

int
main (void)
{
  const gsl_rng_type * T;
  gsl_rng * r;

  int i, n = 10;
  double mu = 3.0;

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

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  /* print n random variates chosen from
     the poisson distribution with mean
     parameter mu */

  for (i = 0; i < n; i++)
    {
      unsigned int k = gsl_ran_poisson (r, mu);
      printf (" %u", k);
    }

  printf ("\n");
  gsl_rng_free (r);
  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* Utility functions for working with GSL data types. */

#include <math.h>
#include <stdio.h>
#include <unistd.h>

#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_statistics.h>
#include <hdf5.h>
#include <hdf5_hl.h>
#include <osqp/osqp.h>

#include "qdm.h"

gsl_vector *
qdm_vector_seq(double from, double to, double by)
{
  size_t size = fabs(to - from) / by;
  gsl_vector *s = gsl_vector_alloc(size + 1);

  double value = from;
  for (size_t i = 0; i < s->size; i++) {
    gsl_vector_set(s, i, value);
    value += by;
  }

  return s;
}

void
qdm_vector_set_seq(gsl_vector *v, double from, double to)
{
  double by = (to - from) / (double)(v->size - 1);

  double value = from;
  for (size_t i = 0; i < v->size; i++) {
    gsl_vector_set(v, i, value);
    value += by;
  }
}

size_t
qdm_vector_search(const gsl_vector *v, double needle)
{
  double x = 0;
  size_t i = v->size - 1;

  for (; i > 0; i--) {
    x = gsl_vector_get(v, i);

    if (needle >= x) {
      break;
    }
  }

  return i;
}

void
qdm_vector_csv_fwrite(FILE *f, const gsl_vector *v)
{
  for (size_t i = 0; i < v->size; i++) {
    fprintf(f, "%.17g", gsl_vector_get(v, i));
    if (i < v->size - 1) {
      fprintf(f, ",");
    }
  }
  fprintf(f, "\n");
}

void
qdm_matrix_csv_fwrite(FILE *f, const gsl_matrix *m)
{
  if (m != NULL) {
    for (size_t i = 0; i < m->size1; i++) {
      gsl_vector_const_view row = gsl_matrix_const_row(m, i);
      qdm_vector_csv_fwrite(f, &row.vector);
    }
  }
}

/* Compute M^T * M */
int
qdm_matrix_tmm(gsl_matrix *m, gsl_matrix *result)
{
  return gsl_blas_dgemm(
      CblasTrans , CblasNoTrans , 1.0 ,
      m          , m            ,
      0.0        , result
  );
}

/* Compute det(M^T * M) */
int
qdm_matrix_det_tmm(gsl_matrix *m, double *det)
{
  int status = 0;
  gsl_matrix *c = gsl_matrix_alloc(m->size2, m->size2);
  gsl_permutation *p = gsl_permutation_alloc(c->size1);

  status = gsl_blas_dgemm(
      CblasTrans , CblasNoTrans , 1.0 ,
      m          , m            ,
      0.0        , c
  );
  if (status != 0) {
    goto cleanup;
  }

  int signum = 0;
  status = gsl_linalg_LU_decomp(c, p, &signum);
  if (status != 0) {
    goto cleanup;
  }

  *det = gsl_linalg_LU_det(c, signum);

cleanup:
  gsl_permutation_free(p);
  gsl_matrix_free(c);

  return status;
}

/* Create a sorted copy of the vector v. */
gsl_vector *
qdm_vector_sorted(const gsl_vector *v)
{
  gsl_vector *s = gsl_vector_alloc(v->size);

  gsl_vector_memcpy(s, v);
  gsl_sort_vector(s);

  return s;
}

gsl_vector *
qdm_vector_quantile(gsl_vector *data, gsl_vector *probs)
{
  gsl_vector *sorted = qdm_vector_sorted(data);
  gsl_vector *quantiles = gsl_vector_alloc(probs->size);

  for (size_t i = 0; i < probs->size; i++) {
    double q = gsl_stats_quantile_from_sorted_data(
        sorted->data,
        sorted->stride,
        sorted->size,
        gsl_vector_get(probs, i)
    );

    gsl_vector_set(quantiles, i, q);
  }

  free(sorted);

  return quantiles;
}

/* Compute the residual sum of squares:
 *
 * sum((y[i] - f(x[i])) ^ 2)
 */
double
qdm_vector_rss(const gsl_vector *y, const gsl_vector *fx)
{
  int status = 0;
  double rss = 0.0;

  gsl_vector *se = gsl_vector_alloc(y->size);
  status = gsl_vector_memcpy(se, y);
  if (status != 0) {
    goto cleanup;
  }

  status = gsl_vector_sub(se, fx);
  if (status != 0) {
    goto cleanup;
  }

  status = gsl_vector_mul(se, se);
  if (status != 0) {
    goto cleanup;
  }

  rss = qdm_vector_sum(se);

cleanup:
  gsl_vector_free(se);

  return rss;
}

/* Compute the summation of the vector.
 *
 * This uses the "iterative Kahan-Babuska algorithm" (aka Klein summation) as
 * outlined here:
 *
 * https://en.wikipedia.org/wiki/Kahan_summation_algorithm
 */
double
qdm_vector_sum(gsl_vector *v)
{
  double s = 0.0;
  double cs = 0.0;
  double ccs = 0.0;

  for (size_t i = 0; i < v->size; i++) {
    double t = s + gsl_vector_get(v, i);

    double c;
    if (fabs(s) >= fabs(gsl_vector_get(v, i))) {
      c = (s - t) + gsl_vector_get(v, i);
    } else {
      c = (gsl_vector_get(v, i) - t) + s;
    }

    s = t;
    t = cs + c;

    double cc;
    if (fabs(cs) >= fabs(c)) {
      cc = (cs - t) + c;
    } else {
      cc = (c - t) + cs;
    }

    cs = t;
    ccs = ccs + cc;
  }

  return s + cs + ccs;
}

/* Return the first index with a value greater than the value. */
size_t
qdm_vector_greater_than(const gsl_vector *v, double value)
{
  for (size_t i = 0; i < v->size; i++) {
    if (gsl_vector_get(v, i) > value) {
      return i;
    }
  }

  return 0;
}

/* Select the elements in the upper triangle of the matrix. All other elements
 * will be set to zero.
 */
void
qdm_matrix_select_upper_triangle(gsl_matrix *m)
{
  for (size_t i = 0; i < m->size1; i++) {
    for (size_t j = 0; j < m->size2; j++) {
      if (j < i) {
        gsl_matrix_set(m, i, j, 0);
      }
    }
  }
}

int
qdm_matrix_to_csc_matrix(csc **result, gsl_matrix *input)
{
  int status = 0;

  gsl_spmatrix *sm = gsl_spmatrix_alloc(input->size1, input->size2);

  status = gsl_spmatrix_d2sp(sm, input);
  if (status != 0) {
    goto cleanup_sm;
  }

  gsl_spmatrix *sm_csc = gsl_spmatrix_compress(sm, GSL_SPMATRIX_CSC);

  int m = sm_csc->size1;
  int n = sm_csc->size2;
  int nzmax = gsl_spmatrix_nnz(sm_csc);

  size_t x_size = sizeof(double) * nzmax;
  double *x = malloc(x_size);
  memcpy(x, sm_csc->data, x_size);

  size_t i_size = sizeof(int) * nzmax;
  int *i = malloc(i_size);
  memcpy(i, sm_csc->i, i_size);

  size_t p_size = sizeof(int) * (n + 1);
  int *p = malloc(p_size);
  memcpy(p, sm_csc->p, p_size);

  *result = csc_matrix(
      m,     // m     First dimension (rows)
      n,     // n     Second dimension (columns)
      nzmax, // nzmax Maximum number of nonzero elements
      x,     // x     Vector of data (size nzmax)
      i,     // i     Vector of row indices (size nzmax)
      p      // p     Vector of column pointers (size n+1)
  );

  gsl_spmatrix_free(sm_csc);

cleanup_sm:
  gsl_spmatrix_free(sm);

  return status;
}

gsl_vector *
qdm_vector_copy(const gsl_vector *src)
{
  gsl_vector *dst = gsl_vector_alloc(src->size);
  gsl_vector_memcpy(dst, src);

  return dst;
}

gsl_matrix *
qdm_matrix_copy(const gsl_matrix *src)
{
  gsl_matrix *dst = gsl_matrix_alloc(src->size1, src->size2);
  gsl_matrix_memcpy(dst, src);

  return dst;
}

int create_hd5(
    const char *file_path,
    const char *group_path
)
{
  int status = 0;

  hid_t file  = -1;
  hid_t group = -1;
  hid_t gcpl  = -1;

  /* Save old error handler */
  herr_t (*old_func)(hid_t, void*) = NULL;
  void *old_client_data = NULL;
  H5Eget_auto(H5E_DEFAULT, &old_func, &old_client_data);

  /* Turn off error handling */
  H5Eset_auto(H5E_DEFAULT, NULL, NULL);

  file = H5Fopen(file_path, H5F_ACC_RDWR, H5P_DEFAULT);
  if (file < 0) {
    file = H5Fcreate(file_path, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
  }
  if (file < 0) {
    status = file;
    goto cleanup;
  }

  group = H5Gopen(file, group_path, H5P_DEFAULT);
  if (group < 0) {
    gcpl = H5Pcreate(H5P_LINK_CREATE);
    if (gcpl < 0) {
      status = gcpl;
      goto cleanup;
    }

    status = H5Pset_create_intermediate_group(gcpl, 1);
    if (status < 0) {
      goto cleanup;
    }

    group = H5Gcreate(file, group_path, gcpl, H5P_DEFAULT, H5P_DEFAULT);

    if (group < 0) {
      status = group;
      goto cleanup;
    }
  }

cleanup:
  if (gcpl >= 0) {
    H5Pclose(gcpl);
  }

  if (group >= 0) {
    H5Gclose(group);
  }

  if (file >= 0) {
    H5Fclose(file);
  }

  /* Restore previous error handler */
  H5Eset_auto(H5E_DEFAULT, old_func, old_client_data);

  if (status < 0) {
    H5Eprint(H5E_DEFAULT, stderr);
  }

  return status;
}

int
qdm_vector_hd5_read(
    hid_t id,
    const char *name,
    gsl_vector **v
)
{
  int status = 0;

  int rank = 0;

  status = H5LTget_dataset_ndims(id, name, &rank);
  if (status < 0) {
    return status;
  }

  hsize_t dims[rank];

  status = H5LTget_dataset_info(id, name, dims, NULL, NULL);
  if (status < 0) {
    return status;
  }

  size_t size = 1;
  for (int i = 0; i < rank; i++) {
    size *= dims[i];
  }

  gsl_vector *tmp = gsl_vector_alloc(size);

  status = H5LTread_dataset_double(id, name, tmp->data);
  if (status < 0) {
    gsl_vector_free(tmp);

    return status;
  }

  *v = tmp;

  return status;
}

int
qdm_vector_hd5_write(
    hid_t id,
    const char *name,
    const gsl_vector *v
)
{
  int status = 0;

  hid_t datatype  = -1;
  hid_t dataspace = -1;
  hid_t dataset   = -1;
  hid_t dcpl      = -1;

  if (v == NULL) {
    goto cleanup;
  }

  datatype = H5Tcopy(H5T_NATIVE_DOUBLE);

  status = H5Tset_order(datatype, H5T_ORDER_LE);
  if (status != 0) {
    goto cleanup;
  }

  hsize_t dims[1] = {
    v->size,
  };
  dataspace = H5Screate_simple(1, dims, NULL);
  if (dataspace < 0) {
    status = dataspace;
    goto cleanup;
  }

  dcpl = H5Pcreate(H5P_DATASET_CREATE);
  if (dcpl < 0) {
    status = dcpl;

    goto cleanup;
  }

  /* Only enable compression on larger vectors. */
  /*
  if (v->size > 1024) {
    status = H5Pset_deflate(dcpl, 9);
    if (status != 0) {
      goto cleanup;
    }

    hsize_t chunk_dims[1] = {
      v->size,
    };
    status = H5Pset_chunk(dcpl, 1, chunk_dims);
    if (status != 0) {
      goto cleanup;
    }
  }
  */

  dataset = H5Dcreate(id, name, datatype, dataspace, H5P_DEFAULT, dcpl, H5P_DEFAULT);
  if (dataset < 0) {
    status = dataset;
    goto cleanup;
  }

  status = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, v->data);
  if (status != 0) {
    goto cleanup;
  }

cleanup:
  if (dcpl >= 0) {
    H5Pclose(dcpl);
  }

  if (dataset >= 0) {
    H5Dclose(dataset);
  }

  if (dataspace >= 0) {
    H5Sclose(dataspace);
  }

  if (datatype >= 0) {
    H5Tclose(datatype);
  }

  H5Oflush(id);

  return status;
}

int
qdm_matrix_hd5_read(
    hid_t id,
    const char *name,
    gsl_matrix **m
)
{
  int status = 0;

  int rank = 0;

  status = H5LTget_dataset_ndims(id, name, &rank);
  if (status < 0) {
    return status;
  }

  if (rank < 2) {
    return -1;
  }

  hsize_t dims[rank];

  status = H5LTget_dataset_info(id, name, dims, NULL, NULL);
  if (status < 0) {
    return status;
  }

  size_t size1 = dims[0];
  size_t size2 = 1;
  for (int i = 1; i < rank; i++) {
    size2 *= dims[i];
  }

  gsl_matrix *tmp = gsl_matrix_alloc(size1, size2);

  status = H5LTread_dataset_double(id, name, tmp->data);
  if (status < 0) {
    gsl_matrix_free(tmp);

    return status;
  }

  *m = tmp;

  return status;
}

int
qdm_matrix_hd5_write(
    hid_t id,
    const char *name,
    const gsl_matrix *m
)
{
  int status = 0;

  hid_t datatype  = -1;
  hid_t dataspace = -1;
  hid_t dataset   = -1;
  hid_t dcpl      = -1;

  if (m == NULL) {
    goto cleanup;
  }

  datatype = H5Tcopy(H5T_NATIVE_DOUBLE);

  status = H5Tset_order(datatype, H5T_ORDER_LE);
  if (status != 0) {
    goto cleanup;
  }

  hsize_t dims[2] = {
    m->size1,
    m->size2,
  };
  dataspace = H5Screate_simple(2, dims, NULL);
  if (dataspace < 0) {
    status = dataspace;

    goto cleanup;
  }

  dcpl = H5Pcreate(H5P_DATASET_CREATE);
  if (dcpl < 0) {
    status = dcpl;

    goto cleanup;
  }

  /* Only enable compression on larger matrices. */
  /*
  if (m->size1 * m->size2 > 1024) {
    status = H5Pset_deflate(dcpl, 9);
    if (status != 0) {
      goto cleanup;
    }

    hsize_t chunk_dims[2] = {
      m->size1,
      m->size2
    };
    status = H5Pset_chunk(dcpl, 2, chunk_dims);
    if (status != 0) {
      goto cleanup;
    }
  }
  */

  dataset = H5Dcreate(id, name, datatype, dataspace, H5P_DEFAULT, dcpl, H5P_DEFAULT);
  if (dataset < 0) {
    status = dataset;

    goto cleanup;
  }

  status = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, m->data);
  if (status != 0) {
    goto cleanup;
  }

cleanup:
  if (dcpl >= 0) {
    H5Pclose(dcpl);
  }

  if (dataset >= 0) {
    H5Dclose(dataset);
  }

  if (dataspace >= 0) {
    H5Sclose(dataspace);
  }

  if (datatype >= 0) {
    H5Tclose(datatype);
  }

  H5Oflush(id);

  return status;
}

gsl_vector *
qdm_matrix_filter(
    const gsl_matrix *m,
    size_t needle_column,
    double needle_value,
    size_t select_column
)
{
  size_t size = 0;

  for (size_t i = 0; i < m->size1; i++) {
    if (gsl_matrix_get(m, i, needle_column) == needle_value) {
      size++;
    }
  }

  gsl_vector *v = gsl_vector_alloc(size);

  size_t j = 0;
  for (size_t i = 0; i < m->size1; i++) {
    double select_value = gsl_matrix_get(m, i, select_column);

    if (gsl_matrix_get(m, i, needle_column) == needle_value) {
      gsl_vector_set(v, j, select_value);
      j++;
    }
  }

  return v;
}

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
// Copyright 2018-2019 TAP, Inc. All Rights Reserved.

#pragma once

#include "targetver.h"

#define WIN32_LEAN_AND_MEAN
#define NOMINMAX
#include <Windows.h>

#include <cstdint>
#include <cassert>
#include <cstdlib>
#include <cstdio>
#include <corecrt_math_defines.h>

#include <array>
#include <vector>
#include <iostream>
#include <sstream>
#include <fstream>
#include <map>
#include <chrono>
#include <filesystem>
#include <thread>
#include <string>
using namespace std::literals::string_literals;

#pragma warning(disable : 4127)
#pragma warning(disable : 4201)
#define GLM_FORCE_RADIANS
#define GLM_FORCE_DEPTH_ZERO_TO_ONE
#include <glm/glm.hpp>
#include <glm/gtc/quaternion.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtc/type_ptr.hpp>
#include <glm/gtx/string_cast.hpp>
#pragma warning(default : 4201)
#pragma warning(default : 4127)

#pragma warning(disable : 4100)
#pragma warning(disable : 4458)
#include <gli/gli.hpp>
#pragma warning(default : 4458)
#pragma warning(default : 4100)

#define VK_USE_PLATFORM_WIN32_KHR
#include <vulkan/vulkan.h>

#include <vulkan/vk_mem_alloc.h>

#include <entt/entt.hpp>

#include <gsl/gsl>


text:
#pragma once

#include "../configuration.h"
#include <chrono>
#include <gsl/gsl-lite.hpp>
#include <string_view>
#include <utility>
#include <yaml-cpp/yaml.h>

namespace angonoka {
/**
    Load Configuration from a YAML string.

    @param text Null-terminated string

    @return An instance of Configuration
*/
Configuration load_text(gsl::czstring text);

/**
    Load Configuration from a YAML file.

    @param path YAML configuration location

    @return An instance of Configuration
*/
Configuration load_file(std::string_view path);
} // namespace angonoka

namespace angonoka::detail {
/**
    Finds or inserts a group into Configuration.groups.

    @param groups   An array of Groups
    @param group    Group name

    @return The index of the group in Configuration.groups
    and whether the insertion took place.
*/
std::pair<GroupIndex, bool>
find_or_insert_group(Groups& groups, std::string_view group);

/**
    Parses human-readable durations.

    Parses durations with resolution ranging from seconds
    up to and including months. Examples:

    1h 15min
    3 weeks and 5 months
    30 seconds

    @param text A string containing a duration.
    @return The duration in seconds.
*/
std::chrono::seconds parse_duration(std::string_view text);

/**
    Parses agents blocks.

    Parses blocks such as these:

    agents:
      agent 1:
      agent 2:

    @param node     "agents" node
    @param config   An instance of Configuration
*/
void parse_agents(const YAML::Node& node, Configuration& config);

/**
    Parses tasks blocks.

    Parses blocks such as these:

    tasks:
      - name: task 1
        duration:
          min: 1 h
          max: 3 h
      - name: task 2
        duration: 2h

    @param node     "tasks" node
    @param config   An instance of Configuration
*/
void parse_tasks(const YAML::Node& node, Configuration& config);
} // namespace angonoka::detail

max_stars_repo_name: coffee-lord/angonoka

max_issues_repo_name: coffee-lord/angonoka

max_forks_repo_name: coffee-lord/angonoka