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

#ifndef __GSL_MATRIX_COMPLEX_FLOAT_H__
#define __GSL_MATRIX_COMPLEX_FLOAT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_complex_float.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  float * data;
  gsl_block_complex_float * block;
  int owner;
} gsl_matrix_complex_float ;

typedef struct
{
  gsl_matrix_complex_float matrix;
} _gsl_matrix_complex_float_view;

typedef _gsl_matrix_complex_float_view gsl_matrix_complex_float_view;

typedef struct
{
  gsl_matrix_complex_float matrix;
} _gsl_matrix_complex_float_const_view;

typedef const _gsl_matrix_complex_float_const_view gsl_matrix_complex_float_const_view;

/* Allocation */

gsl_matrix_complex_float *
gsl_matrix_complex_float_alloc (const size_t n1, const size_t n2);

gsl_matrix_complex_float *
gsl_matrix_complex_float_calloc (const size_t n1, const size_t n2);

gsl_matrix_complex_float *
gsl_matrix_complex_float_alloc_from_block (gsl_block_complex_float * b,
                                           const size_t offset,
                                           const size_t n1, const size_t n2, const size_t d2);

gsl_matrix_complex_float *
gsl_matrix_complex_float_alloc_from_matrix (gsl_matrix_complex_float * b,
                                            const size_t k1, const size_t k2,
                                            const size_t n1, const size_t n2);

gsl_vector_complex_float *
gsl_vector_complex_float_alloc_row_from_matrix (gsl_matrix_complex_float * m,
                                                const size_t i);

gsl_vector_complex_float *
gsl_vector_complex_float_alloc_col_from_matrix (gsl_matrix_complex_float * m,
                                                const size_t j);

void gsl_matrix_complex_float_free (gsl_matrix_complex_float * m);

/* Views */

_gsl_matrix_complex_float_view
gsl_matrix_complex_float_submatrix (gsl_matrix_complex_float * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

_gsl_vector_complex_float_view
gsl_matrix_complex_float_row (gsl_matrix_complex_float * m, const size_t i);

_gsl_vector_complex_float_view
gsl_matrix_complex_float_column (gsl_matrix_complex_float * m, const size_t j);

_gsl_vector_complex_float_view
gsl_matrix_complex_float_diagonal (gsl_matrix_complex_float * m);

_gsl_vector_complex_float_view
gsl_matrix_complex_float_subdiagonal (gsl_matrix_complex_float * m, const size_t k);

_gsl_vector_complex_float_view
gsl_matrix_complex_float_superdiagonal (gsl_matrix_complex_float * m, const size_t k);

_gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_array (float * base,
                             const size_t n1,
                             const size_t n2);

_gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_array_with_tda (float * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

_gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_vector (gsl_vector_complex_float * v,
                              const size_t n1,
                              const size_t n2);

_gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_vector_with_tda (gsl_vector_complex_float * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

_gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_submatrix (const gsl_matrix_complex_float * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

_gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_row (const gsl_matrix_complex_float * m,
                            const size_t i);

_gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_column (const gsl_matrix_complex_float * m,
                               const size_t j);

_gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_diagonal (const gsl_matrix_complex_float * m);

_gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_subdiagonal (const gsl_matrix_complex_float * m,
                                    const size_t k);

_gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_superdiagonal (const gsl_matrix_complex_float * m,
                                      const size_t k);

_gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_array (const float * base,
                                   const size_t n1,
                                   const size_t n2);

_gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_array_with_tda (const float * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

_gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_vector (const gsl_vector_complex_float * v,
                                    const size_t n1,
                                    const size_t n2);

_gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_vector_with_tda (const gsl_vector_complex_float * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

gsl_complex_float gsl_matrix_complex_float_get(const gsl_matrix_complex_float * m, const size_t i, const size_t j);
void gsl_matrix_complex_float_set(gsl_matrix_complex_float * m, const size_t i, const size_t j, const gsl_complex_float x);

gsl_complex_float * gsl_matrix_complex_float_ptr(gsl_matrix_complex_float * m, const size_t i, const size_t j);
const gsl_complex_float * gsl_matrix_complex_float_const_ptr(const gsl_matrix_complex_float * m, const size_t i, const size_t j);

void gsl_matrix_complex_float_set_zero (gsl_matrix_complex_float * m);
void gsl_matrix_complex_float_set_identity (gsl_matrix_complex_float * m);
void gsl_matrix_complex_float_set_all (gsl_matrix_complex_float * m, gsl_complex_float x);

int gsl_matrix_complex_float_fread (FILE * stream, gsl_matrix_complex_float * m) ;
int gsl_matrix_complex_float_fwrite (FILE * stream, const gsl_matrix_complex_float * m) ;
int gsl_matrix_complex_float_fscanf (FILE * stream, gsl_matrix_complex_float * m);
int gsl_matrix_complex_float_fprintf (FILE * stream, const gsl_matrix_complex_float * m, const char * format);

int gsl_matrix_complex_float_memcpy(gsl_matrix_complex_float * dest, const gsl_matrix_complex_float * src);
int gsl_matrix_complex_float_swap(gsl_matrix_complex_float * m1, gsl_matrix_complex_float * m2);

int gsl_matrix_complex_float_swap_rows(gsl_matrix_complex_float * m, const size_t i, const size_t j);
int gsl_matrix_complex_float_swap_columns(gsl_matrix_complex_float * m, const size_t i, const size_t j);
int gsl_matrix_complex_float_swap_rowcol(gsl_matrix_complex_float * m, const size_t i, const size_t j);

int gsl_matrix_complex_float_transpose (gsl_matrix_complex_float * m);
int gsl_matrix_complex_float_transpose_memcpy (gsl_matrix_complex_float * dest, const gsl_matrix_complex_float * src);

int gsl_matrix_complex_float_isnull (const gsl_matrix_complex_float * m);

int gsl_matrix_complex_float_add (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
int gsl_matrix_complex_float_sub (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
int gsl_matrix_complex_float_mul_elements (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
int gsl_matrix_complex_float_div_elements (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
int gsl_matrix_complex_float_scale (gsl_matrix_complex_float * a, const gsl_complex_float x);
int gsl_matrix_complex_float_add_constant (gsl_matrix_complex_float * a, const gsl_complex_float x);
int gsl_matrix_complex_float_add_diagonal (gsl_matrix_complex_float * a, const gsl_complex_float x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
int gsl_matrix_complex_float_get_row(gsl_vector_complex_float * v, const gsl_matrix_complex_float * m, const size_t i);
int gsl_matrix_complex_float_get_col(gsl_vector_complex_float * v, const gsl_matrix_complex_float * m, const size_t j);
int gsl_matrix_complex_float_set_row(gsl_matrix_complex_float * m, const size_t i, const gsl_vector_complex_float * v);
int gsl_matrix_complex_float_set_col(gsl_matrix_complex_float * m, const size_t j, const gsl_vector_complex_float * v);

#ifdef HAVE_INLINE

extern inline
gsl_complex_float
gsl_matrix_complex_float_get(const gsl_matrix_complex_float * m,
                     const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  gsl_complex_float zero = {{0,0}};

  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, zero) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, zero) ;
    }
#endif
  return *(gsl_complex_float *)(m->data + 2*(i * m->tda + j)) ;
}

extern inline
void
gsl_matrix_complex_float_set(gsl_matrix_complex_float * m,
                     const size_t i, const size_t j, const gsl_complex_float x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  *(gsl_complex_float *)(m->data + 2*(i * m->tda + j)) = x ;
}

extern inline
gsl_complex_float *
gsl_matrix_complex_float_ptr(gsl_matrix_complex_float * m,
                             const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (gsl_complex_float *)(m->data + 2*(i * m->tda + j)) ;
}

extern inline
const gsl_complex_float *
gsl_matrix_complex_float_const_ptr(const gsl_matrix_complex_float * m,
                                   const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const gsl_complex_float *)(m->data + 2*(i * m->tda + j)) ;
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_MATRIX_COMPLEX_FLOAT_H__ */

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


text:
 /**
  * File:    utils.h
  *
  * Author:  Vishal R (vishalr@pesu.pes.edu or vishalramesh01@gmail.com)
  *
  * Summary of File:
  *   Matrix operation library required for cdnn. Contains all the header files for matrix operations.
  */

#ifndef UTILS_H
#define UTILS_H

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <omp.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include <errno.h>
#include <cblas.h>
#include <signal.h>
#include <sys/time.h>
#include <sys/stat.h>
#include <termcap.h>

#ifdef __cplusplus
extern "C"{
#endif

  typedef struct array{
    float * matrix;
    int shape[2];
  }dARRAY;

  dARRAY * zeros(int * dims);
  dARRAY * ones(int * dims);
  dARRAY * eye(int * dims);
  dARRAY * randn(int * dims);
  dARRAY * add(dARRAY * MatrixA, dARRAY * MatrixB);
  dARRAY * addScalar(dARRAY * matrix, float scalar);
  dARRAY * subtract(dARRAY * MatrixA, dARRAY * MatrixB);
  dARRAY * subScalar(dARRAY * matrix, float scalar);
  dARRAY * sum(dARRAY * matrix, int axis);
  dARRAY * transpose(dARRAY *  Matrix);
  dARRAY * dot(dARRAY *  MatrixA, dARRAY *  MatrixB);
  dARRAY * multiply(dARRAY * MatrixA, dARRAY * MatrixB);
  dARRAY * mulScalar(dARRAY * matrix, float scalar);
  dARRAY * divison(dARRAY * MatrixA, dARRAY * MatrixB);
  dARRAY * divScalar(dARRAY * matrix, float scalar);
  dARRAY * power(dARRAY * matrix, float power);
  dARRAY * squareroot(dARRAY * matrix);
  dARRAY * exponentional(dARRAY * matrix);
  dARRAY * b_cast(dARRAY * MatrixA, dARRAY * MatrixB);
  dARRAY * reshape(dARRAY * matrix, int * dims);

  int * permutation(int length);
  float mean(dARRAY * matrix);
  float var(dARRAY * matrix, char * type);
  float std(dARRAY * matrix, char * type);
  float gaussGenerator(float * cache, int * return_cache);
  float gaussRandom();
  float rand_norm(float mu, float std);

  float frobenius_norm(dARRAY * matrix);
  float Manhattan_distance(dARRAY * matrix);

  int size(dARRAY * A);
  void shape(dARRAY * A);
  void free2d(dARRAY * matrix);
  void sleep_my(int milliseconds);
  void cleanSTDIN();
  void get_safe_nn_threads();
#ifdef __cplusplus
}
#endif
#endif //UTILS_H

ext: h


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * 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 SIM_MATH_ODE_SOLVER_GSL_H_
#define SIM_MATH_ODE_SOLVER_GSL_H_

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#include <stdint.h>

#include <vector>

#include "common/macros.h"
#include "sim/math/ode_solver.h"
#include "sim/sim_types.h"

namespace sim {

// Wrapper for the GSL ODE library.
class GslOdeSolver : public OdeSolver {
 public:
  explicit GslOdeSolver(const OdeSystem &ode_system,
                        const SimOdeSolverParams &params);
  ~GslOdeSolver() {
    if (ode_driver_ != nullptr) gsl_odeiv2_driver_free(ode_driver_);
  }

  OdeSolverStatus Integrate(double t0, double tf, const std::vector<double> &x0,
                            double *t_int, std::vector<double> *x) override;

 private:
  // Static callback function for the GSL ODE solver.
  //
  // Args:
  //   t: Time at which derivative will be evaluated.
  //   x: Array of length num_states() containing the state at which
  //       the derivative will be evaluated.
  //   dx: Array of length num_states() into which the derivative is stored.
  //   context: Pointer to the OdeSolver class containing the OdeSystem.
  //
  // Returns:
  //   GSL_SUCCESS if the derivative was calculated successfully,
  //   GSL_FAILURE if the time step was too large, or
  //   GSL_EBADFUNC if the integration should be aborted immediately.
  static int32_t GslCallback(double t, const double x[], double dx[],
                             void *context);

  // Solver parameters.
  const SimOdeSolverParams &params_;

  // ODE to be solved.
  const OdeSystem &ode_system_;

  // Parameters for GSL.
  gsl_odeiv2_system sys_;
  gsl_odeiv2_driver *ode_driver_;

  DISALLOW_COPY_AND_ASSIGN(GslOdeSolver);
};

}  // namespace sim

#endif  // SIM_MATH_ODE_SOLVER_GSL_H_


text:
/*
*/

#ifndef MINIMIZER_H
#define MINIMIZER_H

#include <iostream>
#include "Vector.h"
#include <stdexcept>

#ifdef USE_GSL
extern "C" {
#include <gsl/gsl_multimin.h>
}

namespace PsimagLite {

template<typename FunctionType>
typename FunctionType::FieldType myFunction(const gsl_vector *v, void *params)
{
	FunctionType* ft = (FunctionType *)params;
	return ft->operator()(v->data,v->size);
}

template<typename FunctionType>
void myDfunction(const gsl_vector *v,
                 void *params,
                 gsl_vector* df)
{
	FunctionType* ft = (FunctionType *)params;
	ft->df(v->data,v->size,df->data,df->size);
}

template<typename FunctionType>
void myFdFunction(const gsl_vector *v,
                  void *params,
                  double *f,
                  gsl_vector *df)
{
	*f = myFunction<FunctionType>(v, params);
	myDfunction<FunctionType>(v, params, df);
}

template<typename RealType,typename FunctionType>
class Minimizer {

	typedef typename FunctionType::FieldType FieldType;
	typedef typename Vector<FieldType>::Type VectorType;
	typedef Minimizer<RealType,FunctionType> ThisType;

public:

	enum {GSL_SUCCESS=::GSL_SUCCESS, GSL_CONTINUE=::GSL_CONTINUE};

	Minimizer(FunctionType& function,SizeType maxIter, bool verbose = false)
	    : function_(function),
	      maxIter_(maxIter),
	      verbose_(verbose),
	      status_(100),
	      gslT_(gsl_multimin_fminimizer_nmsimplex2),
	      gslS_(gsl_multimin_fminimizer_alloc(gslT_,function_.size())),
	      gslDt_(gsl_multimin_fdfminimizer_conjugate_fr),
	      gslDs_(gsl_multimin_fdfminimizer_alloc (gslDt_, function_.size()))
	{}

	~Minimizer()
	{
		gsl_multimin_fminimizer_free(gslS_);
		gsl_multimin_fdfminimizer_free(gslDs_);
	}

	int simplex(VectorType& minVector,RealType delta=1e-3,RealType tolerance=1e-3)
	{
		gsl_vector *x;
		/* Starting point,  */
		x = gsl_vector_alloc (function_.size());
		for (SizeType i=0;i<minVector.size();i++)
			gsl_vector_set (x, i, minVector[i]);

		gsl_vector *xs;
		xs = gsl_vector_alloc (function_.size());
		for (SizeType i=0;i<minVector.size();i++)
			gsl_vector_set (xs, i, delta);

		gsl_multimin_function func;
		func.f= myFunction<FunctionType>;
		func.n = function_.size();
		func.params = &function_;
		gsl_multimin_fminimizer_set (gslS_, &func, x, xs);

		SizeType iter = 0;
		RealType prevValue = 0;

		for (;iter<maxIter_;iter++) {
			status_ = gsl_multimin_fminimizer_iterate (gslS_);

			if (status_) {
				String gslError(gsl_strerror(status_));
				String msg("Minimizer::simplex(...): GSL Error: ");
				msg += gslError + "\n";
				throw RuntimeError(msg);
			}

			RealType size = gsl_multimin_fminimizer_size(gslS_);
			status_ = gsl_multimin_test_size(size, tolerance);

			if (verbose_) {
				RealType thisValue = function_(gslS_->x->data,func.n);
				RealType diff = fabs(thisValue - prevValue);
				std::cerr<<"simplex: "<<iter<<" "<<thisValue<<" diff= "<<diff;
				std::cerr<<" status= "<<status_<<" size="<<size<<"\n";
				prevValue = thisValue;
			}

			if (status_ == GSL_SUCCESS) break;
		}

		found(minVector,gslS_->x,iter);
		gsl_vector_free (x);
		gsl_vector_free (xs);
		return iter;
	}

	int conjugateGradient(VectorType& minVector,
	                      RealType delta=1e-3,
	                      RealType delta2=1e-3,
	                      RealType tolerance=1e-3,
	                      SizeType saveEvery = 0)
	{
		gsl_vector *x;
		/* Starting point,  */
		x = gsl_vector_alloc (function_.size());
		for (SizeType i=0;i<minVector.size();i++)
			gsl_vector_set (x, i, minVector[i]);

		gsl_multimin_function_fdf func;
		func.f= myFunction<FunctionType>;
		func.df = myDfunction<FunctionType>;
		func.fdf = myFdFunction<FunctionType>;
		func.n = function_.size();
		func.params = &function_;

		gsl_multimin_fdfminimizer_set(gslDs_, &func, x, delta, delta2);

		RealType prevValue = 0;
		SizeType iter = 0;
		for (;iter<maxIter_;iter++) {
			status_ = gsl_multimin_fdfminimizer_iterate(gslDs_);

			if (status_) {
				String gslError(gsl_strerror(status_));
				String msg("Minimizer::conjugateGradient(...): GSL Error: ");
				msg += gslError + "\n";
				throw RuntimeError(msg);
			}

			status_ = gsl_multimin_test_gradient (gslDs_->gradient, tolerance);

			if (verbose_) {
				RealType thisValue = function_(gslDs_->x->data,func.n);
				RealType diff = fabs(thisValue - prevValue);
				std::cerr<<"conjugateGradient: "<<iter<<" "<<thisValue;
				std::cerr<<" diff= "<<diff;
				std::cerr<<" gradientNorm= "<<gradientNorm(gslDs_->x);
				std::cerr<<" status= "<<status_<<"\n";
				prevValue = thisValue;
			}

			if (status_ == GSL_SUCCESS) break;

			if (saveEvery > 0 && iter%saveEvery==0)
				printIntermediate(gslDs_->x, iter);
		}

		found(minVector,gslDs_->x,iter);
		gsl_vector_free (x);
		return iter;
	}

	int status() const { return status_; }

	String statusString() const
	{
		switch (status_) {
		case GSL_SUCCESS:
			return "GSL_SUCCESS";
			break;
		case GSL_CONTINUE:
			return "GSL_CONTINUE";
			break;
		default:
			return "UNKNOWN";
			break;
		}
	}

private:

	void found(VectorType& minVector,gsl_vector* x,SizeType iter)
	{
		for (SizeType i=0;i<minVector.size();i++)
			minVector[i] = gsl_vector_get(x,i);
	}

	void printIntermediate(gsl_vector* x, SizeType iter) const
	{
		std::cerr<<"INTERMEDIATE "<<iter<<"\n";
		std::cerr<<x->size<<"\n";
		for (SizeType i=0;i<x->size;i++)
			std::cerr<<gsl_vector_get(x,i)<<"\n";
	}

	RealType gradientNorm(const gsl_vector *v) const
	{
		gsl_vector* df = gsl_vector_alloc (function_.size());
		myDfunction<FunctionType>(v,&function_,df);
		RealType sum = 0;
		for (int i = 0; i < df->size; ++i) {
			RealType tmp = gsl_vector_get(df,i);
			sum += PsimagLite::conj(tmp) * tmp;
		}

		gsl_vector_free (df);
		return sqrt(sum);
	}

	FunctionType& function_;
	SizeType maxIter_;
	bool verbose_;
	int status_;
	const gsl_multimin_fminimizer_type *gslT_;
	gsl_multimin_fminimizer *gslS_;
	const gsl_multimin_fdfminimizer_type *gslDt_;
	gsl_multimin_fdfminimizer *gslDs_;

}; // class Minimizer
}

#else

namespace PsimagLite {

template<typename RealType,typename FunctionType>
class Minimizer {

	typedef typename FunctionType::FieldType FieldType;
	typedef typename Vector<FieldType>::Type VectorType;

public:

	enum {GSL_SUCCESS=0, GSL_CONTINUE=1};

	Minimizer(FunctionType&,SizeType, bool = false)
	{
		String str("Minimizer needs the gsl\n");
		throw RuntimeError(str);
	}

	int simplex(VectorType&,RealType=1e-3,RealType=1e-3)
	{
		String str("Minimizer needs the gsl\n");
		throw RuntimeError(str);
	}

	int status() const { return 1; }

	String statusString() const
	{
		return "Minimizer needs the gsl";
	}
};

} // namespace PsimagLite
#endif

#endif // MINIMIZER_H


text:
#pragma once

#include <memory>
#include <gsl/span>

#include "halley/time/halleytime.h"
#include "../session/network_session.h"
#include "entity_network_remote_peer.h"
#include "halley/bytes/serialization_dictionary.h"
#include "halley/entity/system.h"
#include "halley/entity/world.h"

namespace Halley {
	class EntityFactory;
	class Resources;
	class World;
	class NetworkSession;

	class EntitySessionSharedData : public SharedData {
	public:

	};

	class EntityClientSharedData : public SharedData {
	public:
		std::optional<Rect4i> viewRect;

		void serialize(Serializer& s) const override;
		void deserialize(Deserializer& s) override;
	};

	class EntityNetworkSession : NetworkSession::IListener, NetworkSession::ISharedDataHandler, public IWorldNetworkInterface {
    public:
		class IEntityNetworkSessionListener {
		public:
			virtual ~IEntityNetworkSessionListener() = default;
			virtual void onStartSession(NetworkSession::PeerId myPeerId) = 0;
			virtual void onRemoteEntityCreated(EntityRef entity, NetworkSession::PeerId peerId) {}
			virtual void setupInterpolators(DataInterpolatorSet& interpolatorSet, EntityRef entity, bool remote) = 0;
			virtual bool isEntityInView(EntityRef entity, const EntityClientSharedData& clientData) = 0;
		};

		EntityNetworkSession(std::shared_ptr<NetworkSession> session, Resources& resources, std::set<String> ignoreComponents, IEntityNetworkSessionListener* listener);
		~EntityNetworkSession() override;

		void setWorld(World& world, SystemMessageBridge bridge);

		void sendUpdates(Time t, Rect4i viewRect, gsl::span<const std::pair<EntityId, uint8_t>> entityIds); // Takes pairs of entity id and owner peer id
		void receiveUpdates();

		World& getWorld() const;
		EntityFactory& getFactory() const;
		NetworkSession& getSession() const;
		bool hasWorld() const;

		const EntityFactory::SerializationOptions& getEntitySerializationOptions() const;
		const EntityDataDelta::Options& getEntityDeltaOptions() const;
		const SerializerOptions& getByteSerializationOptions() const;
		SerializationDictionary& getSerializationDictionary();

		Time getMinSendInterval() const;

		void onRemoteEntityCreated(EntityRef entity, NetworkSession::PeerId peerId);
		void requestSetupInterpolators(DataInterpolatorSet& interpolatorSet, EntityRef entity, bool remote);
		void setupOutboundInterpolators(EntityRef entity);

		bool isReadyToStart() const;
		bool isEntityInView(EntityRef entity, const EntityClientSharedData& clientData) const;

		Vector<Rect4i> getRemoteViewPorts() const;

		bool isHost() override;
		bool isRemote(ConstEntityRef entity) const override;
		void sendEntityMessage(EntityRef entity, int messageType, Bytes messageData) override;
		void sendSystemMessage(String targetSystem, int messageType, Bytes messageData, SystemMessageDestination destination, SystemMessageCallback callback) override;

		void sendToAll(EntityNetworkMessage msg);
		void sendToPeer(EntityNetworkMessage msg, NetworkSession::PeerId peerId);

	protected:
		void onStartSession(NetworkSession::PeerId myPeerId) override;
		void onPeerConnected(NetworkSession::PeerId peerId) override;
		void onPeerDisconnected(NetworkSession::PeerId peerId) override;
		std::unique_ptr<SharedData> makeSessionSharedData() override;
		std::unique_ptr<SharedData> makePeerSharedData() override;

	private:
		struct QueuedMessage {
			NetworkSession::PeerId fromPeerId;
			EntityNetworkMessage message;
		};

		struct PendingSysMsgResponse {
			SystemMessageCallback callback;
		};

		Resources& resources;
		std::shared_ptr<EntityFactory> factory;
		IEntityNetworkSessionListener* listener = nullptr;
		SystemMessageBridge messageBridge;
		uint32_t systemMessageId = 0;
		HashMap<uint32_t, PendingSysMsgResponse>  pendingSysMsgResponses;

		EntityFactory::SerializationOptions entitySerializationOptions;
		EntityDataDelta::Options deltaOptions;
		SerializerOptions byteSerializationOptions;
		SerializationDictionary serializationDictionary;

		std::shared_ptr<NetworkSession> session;
		Vector<EntityNetworkRemotePeer> peers;

		Vector<QueuedMessage> queuedPackets;

		HashMap<int, Vector<EntityNetworkMessage>> outbox;

		bool readyToStart = false;

		bool canProcessMessage(const EntityNetworkMessage& msg) const;
		void processMessage(NetworkSession::PeerId fromPeerId, EntityNetworkMessage msg);
		void onReceiveEntityUpdate(NetworkSession::PeerId fromPeerId, EntityNetworkMessage msg);
		void onReceiveReady(NetworkSession::PeerId fromPeerId, const EntityNetworkMessageReadyToStart& msg);
		void onReceiveMessageToEntity(NetworkSession::PeerId fromPeerId, const EntityNetworkMessageEntityMsg& msg);
		void onReceiveSystemMessage(NetworkSession::PeerId fromPeerId, const EntityNetworkMessageSystemMsg& msg);
		void onReceiveSystemMessageResponse(NetworkSession::PeerId fromPeerId, const EntityNetworkMessageSystemMsgResponse& msg);

		void sendMessages();

		void setupDictionary();
	};
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


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/constant.h"
#include "chainerx/error.h"
#include "chainerx/optional_container_arg.h"
#include "chainerx/stack_vector.h"

namespace chainerx {

class Axes : public StackVector<int8_t, kMaxNdim> {
    using BaseVector = StackVector<int8_t, kMaxNdim>;

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.

    Axes() = default;

    ~Axes() = default;

    // by iterators
    template <typename InputIt>
    Axes(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 Axes(absl::Span<const int8_t> axes) : Axes{axes.begin(), axes.end()} {}

    // by initializer list
    Axes(std::initializer_list<int8_t> axes) : Axes{axes.begin(), axes.end()} {}

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

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

    std::string ToString() const;

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

    int8_t& operator[](int8_t index) {
        if (!(0 <= index && static_cast<size_t>(index) < size())) {
            throw DimensionError{"index out of bounds"};
        }
        return this->StackVector::operator[](index);
    }

    const int8_t& operator[](int8_t index) const {
        if (!(0 <= index && static_cast<size_t>(index) < size())) {
            throw DimensionError{"index out of bounds"};
        }
        return this->StackVector::operator[](index);
    }

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

std::ostream& operator<<(std::ostream& os, const Axes& axes);

using OptionalAxes = OptionalContainerArg<Axes>;

namespace internal {

bool IsAxesPermutation(const Axes& axes, int8_t ndim);

// Normalizes possibly-negative axis to non-negative axis in [0, ndim).
// If `axis` does not fit in [-ndim, ndim), DimensionError is thrown.
int8_t NormalizeAxis(int8_t axis, int8_t ndim);

// Resolves the axis argument of many operations.
// Negative axes are converted to non-negative ones (by wrapping at ndim).
Axes GetNormalizedAxes(const Axes& axis, int8_t ndim);

// Resolves the axis argument of many operations.
// Negative axes are converted to non-negative ones (by wrapping at ndim).
// Axes are then sorted.
Axes GetSortedAxes(const Axes& axis, int8_t ndim);

// Resolves the axis argument of many operations.
// Negative axes are converted to non-negative ones (by wrapping at ndim).
// Axes are then sorted.
// nullopt is converted to a vector of all axes.
Axes GetSortedAxesOrAll(const OptionalAxes& axis, int8_t ndim);

}  // namespace internal
}  // namespace chainerx

max_stars_repo_name: zaltoprofen/chainer

max_issues_repo_name: zaltoprofen/chainer

max_forks_repo_name: zaltoprofen/chainer


text:
// Alessandro Casalino
//
// Compile with (Mac with default homebrew gsl 2.6) gcc-9 -O2 coupled_quintessence_evolve.c -o coupled_quintessence_evolve.exe -L/usr/local/Cellar/gsl/2.6/lib -I/usr/local/Cellar/gsl/2.6/include -lgsl
// Run with ./coupled_quintessence_evolve.exe

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

// PHYSICAL PARAMETERS VALUES

// INITIAL CONDITIONS
// Initial value of the scale factor
double a_init = 1e-14;
// Final value of the scale factor
double a_end = 1.1;

// Values of fractional density for the cosmological matter today
double Omega_rad_0 = 8e-5;
double Omega_b_0 = 0.0486;//0.0490439;
double Omega_Lambda_0 = 0.;//0.6911;
double Omega_cdm_0 = 0.2589;//0.264201;

// Physical constants for output conversions
double _G_ = 6.67428e-11;                 /**< Newton constant in m^3/Kg/s^2 */
double _MPc_over_m_ = 3.085677581282e22;  // Conversion factor from Mpc to m
double _Gyr_over_Mpc_ = 3.06601394e2;     // Conversion factor from Gyr to Mpc
double _c_ = 2.99792458e8;                // Speed of light in m/s
double _H0_ = 67.74;                      // H0 of LCDM in Km/s/Mpc (h:0.6731)

double fourpiG = 0.5;

double T_CONF = 1.;

// COMPUTATION PARAMETERS VALUES

// Number of points used in the computation
int points = (int) 1e6;
// Max particle horizon value
double max_tau = 40.;
// Min particle horizon value
double min_tau = 1e-9;

// Raise this value to make the csv file smaller, but decreasing resolution
// The value inserted is the ratio between the number of values written in a full resolution file / decreased resolution file
int csv_resolution = 10;

// BISECTION PARAMTERS
// Bisection precision
double delta_bisection = 1e-2;

// Bisection target
//  1: mg_field
//  2: mg_field_p
//  3: potential_constant
int bisection_type = 3;

double bisection_min = 0.;
double bisection_max = 20.;

// Initial conditions (considered accordingly to the bisection_type choice)
double mg_field_bk_0 = 1e-20;
double mg_field_p_bk_0 = 1.;

// TEST MODE
// Provides some informations on the terminal and the .csv file during the computation (1: on, others: off)
int TEST_MODE = 1;

// Definition of the POTENTIAL
// For a list of the models see above
double mg_pot_const = 1.;
double mg_pot_exp = -0.5;

double mg_pot(const double mg_field_bk) {

  return mg_pot_const * pow(mg_field_bk, mg_pot_exp);

}

double mg_pot_p(const double mg_field_bk) {

  return mg_pot_exp * mg_pot_const * pow(mg_field_bk, mg_pot_exp - 1.);

}

// Definition of the COUPLING FUNCION
// For a list of the models see above
double mg_coupl_const = 0.1;

double mg_coupl(const double mg_field_bk) {

  // 2 factor to be consistent with hi_class units (8 pi G = 1)
  return mg_coupl_const * mg_field_bk * mg_field_bk;

}

double mg_coupl_p(const double mg_field_bk) {

  return 2.*mg_coupl_const * mg_field_bk;

}

double mg_coupl_pp(const double mg_field_bk) {

  return 2.*mg_coupl_const;

}

// Definitions of DIFFERENTIAL EQUATION system
// Need to have first order equations to use Runge Kutta
//
double a_p_rk4(const double a, const double a_p, const double mg_field_bk, const double mg_field_p_bk)
{
	return a_p;
}
double a_pp_rk4(const double a, const double a_p, const double mg_field_bk, const double mg_field_p_bk)
{
  double a3 = a * a * a;
  double rhoa3 = (Omega_cdm_0 + Omega_b_0) + (Omega_Lambda_0 * a3) + (Omega_rad_0 / a);
  double Pa3 = - (Omega_Lambda_0 * a3) + 1./3. * (Omega_rad_0 / a );

  return fourpiG / 3. * ( (rhoa3 - 3. * Pa3) - a * mg_field_p_bk * mg_field_p_bk * (1. + 3./2./fourpiG * mg_coupl_pp(mg_field_bk)) + 4. * a3 * mg_pot(mg_field_bk) + 3./2./fourpiG * a3 * mg_pot_p(mg_field_bk) * mg_coupl_p(mg_field_bk))/(1. + mg_coupl(mg_field_bk) + 3./4./fourpiG * mg_coupl_p(mg_field_bk) * mg_coupl_p(mg_field_bk));
}
double mg_field_bk_p_rk4(const double a, const double a_p, const double mg_field_bk, const double mg_field_p_bk)
{
	return mg_field_p_bk;
}
double mg_field_bk_pp_rk4(const double a, const double a_p, const double mg_field_bk, const double mg_field_p_bk)
{
  double a2 = a * a;
  double rhoa2 = (Omega_cdm_0 + Omega_b_0) / a + (Omega_Lambda_0 * a2) + (Omega_rad_0 / a / a);
  double Pa2 = - (Omega_Lambda_0 * a2) + 1./3. * (Omega_rad_0 / a / a );

	return - 2. * a_p/a * mg_field_p_bk + ( - mg_pot_p(mg_field_bk) * a2 * ( 1. + mg_coupl(mg_field_bk) ) + ( rhoa2 - 3. * Pa2 + 4. * a2 * mg_pot(mg_field_bk) - mg_field_p_bk * mg_field_p_bk * ( 1. + 3./2./fourpiG * mg_coupl_pp(mg_field_bk) ) ) /2. * mg_coupl_p(mg_field_bk) ) / (1. + mg_coupl(mg_field_bk) + 3./4./fourpiG * mg_coupl_p(mg_field_bk) * mg_coupl_p(mg_field_bk) );
}
double Hconf(const double a, const double mg_field_bk, const double mg_field_p_bk)
{
  double a2 = a * a;
  double rhoa2 = (Omega_cdm_0 + Omega_b_0) / a + (Omega_Lambda_0 * a2) + (Omega_rad_0 / a / a);

  return - mg_field_p_bk * mg_coupl_p(mg_field_bk)/(1. + mg_coupl(mg_field_bk))/2. + sqrt((2. * fourpiG / 3.) * (1. + mg_coupl(mg_field_bk)) * ( rhoa2 + (mg_field_p_bk * mg_field_p_bk / 2.) + (a2 * mg_pot(mg_field_bk)) ) + mg_field_p_bk * mg_field_p_bk/4. * mg_coupl_p(mg_field_bk) * mg_coupl_p(mg_field_bk))/(1. + mg_coupl(mg_field_bk));
}

// Integrand for the particle horizon integral
double particleHorizonIntegrand(double a, double mg_field_bk, double mg_field_p_bk)
{
	return 1. / ( a * Hconf(a, mg_field_bk, mg_field_p_bk) );
}
// Particle horizon integral step
double particleHorizon(const int i, double * a, double * mg_field_bk, double * mg_field_p_bk) {
  double h = a[i]-a[i-1];
  double fa = particleHorizonIntegrand(a[i-1],mg_field_bk[i-1],mg_field_p_bk[i-1]);
  double fb = particleHorizonIntegrand(a[i],mg_field_bk[i],mg_field_p_bk[i]);

  return h*(fa+fb)/2.;
}

int scan_for_a0 (double * a, double ae) {

  int i=0;

  for(i=0; a[i] <= ae && i < points; i++);

  return i;

}

// Function used to print results stored in vectors as a csv file
void csv(double * t, double * a, double * a_p, double * mg_field_bk, double * mg_field_p_bk, double * particleHorizonVec, char * filename) {

    FILE *fp;
    fp = fopen (filename, "w+");
    fprintf(fp, "%s, %s, %s, %s, %s, %s, %s, %s", "t", "a(t)", "H(t)/H0", "H_prime(t)/H0^2", "Omega_df", "Omega_r", "Omega_b", "Omega_cdm");
    fprintf(fp, ", %s", "omega_df");
    fprintf(fp, ", %s", "PH"); // Particle horizon
    fprintf(fp, ", %s", "mg_field");
    fprintf(fp, ", %s", "mg_field_p");
    if(TEST_MODE==1) fprintf(fp, ", %s", "H_check(t)/H0");
    fprintf(fp, "\n");

    // Time conversion factor
    // double tcf = 1./_H0_/(60.*60.*24.*365.*1e9)*_MPc_over_m_/1000.;

    int i = 0;
    double particleHorizonRes = 0.;

    for(i=1;i<points;i++){
      particleHorizonRes += particleHorizon(i,a,mg_field_bk,mg_field_p_bk);
      particleHorizonVec[i] = particleHorizonRes;
    }

    i = 0;

    int i_a0 = scan_for_a0(a,1.);
    double H0 = Hconf(a[i_a0],mg_field_bk[i_a0],mg_field_p_bk[i_a0]);

    int w1=0, w2=0, w3=0;

    while( a[i] <= a_end ){

      double H = Hconf(a[i],mg_field_bk[i],mg_field_p_bk[i]); // this is the Hubble constant with conformal time
      double H_prime = a_pp_rk4(a[i], a_p[i], mg_field_bk[i], mg_field_p_bk[i]) / a[i] - H * H; // this is Hubble prime with conformal time
      double phi_pp = mg_field_bk_pp_rk4(a[i], a_p[i], mg_field_bk[i], mg_field_p_bk[i]);

      double P_df   = - mg_pot(mg_field_bk[i]) + mg_field_p_bk[i] * mg_field_p_bk[i]/2./a[i]/a[i] + 1./(2.*fourpiG)*( mg_coupl(mg_field_bk[i])*(H*H+2.*H_prime) + (H*mg_field_p_bk[i]+phi_pp) * mg_coupl_p(mg_field_bk[i]) + mg_field_p_bk[i] * mg_field_p_bk[i] * mg_coupl_pp(mg_field_bk[i]))/a[i]/a[i];
      double rho_df = mg_pot(mg_field_bk[i]) + mg_field_p_bk[i] * mg_field_p_bk[i]/2./a[i]/a[i] - 3./2./fourpiG * ( H * mg_field_p_bk[i] * mg_coupl_p(mg_field_bk[i]) + H * H * mg_coupl(mg_field_bk[i]))/a[i]/a[i];

      double Omega_df  = 2. * fourpiG / 3. * a[i] * a[i] * rho_df /H /H;
      double Omega_rad = 2. * fourpiG / 3. * Omega_rad_0 /a[i] /a[i] /H /H;
      double Omega_b = 2. * fourpiG / 3. * Omega_b_0 /a[i] /H /H;
      double Omega_cdm = 2. * fourpiG / 3. * Omega_cdm_0 /a[i] /H /H;

      if(rho_df<0 && w1==0 && i>0){
        printf("WARNING: the energy density is negative at a = %e\n", a[i]);
        w1++;
      }
      if( (Omega_df+Omega_rad+Omega_b+Omega_cdm > 1.01 || Omega_df+Omega_rad+Omega_b+Omega_cdm < 0.99) && w2==0){
        printf("WARNING: the fractional density sum is %e at a = %e\n", Omega_df+Omega_rad+Omega_b+Omega_cdm,a[i]);
        w2++;
      }
      if( (Omega_df> 1. || Omega_rad > 1. || Omega_b > 1. || Omega_cdm > 1.) && w3==0){
        printf("WARNING: one fractional density is > 1 at a = %e\n", a[i]);
        w3++;
      }

      double omega_df = P_df/rho_df;

      if(T_CONF==0) H_prime = ( H_prime - H * H )/a[i]/a[i]; // this is Hubble prime with cosmological time
      if(T_CONF==0) H = H / a[i]; // this it the Hubble constant with cosmological time

      fprintf(fp, "%e, %e, %e, %e, %e, %e, %e, %e", t[i], a[i], H / H0, H_prime /H0 /H0, Omega_df, Omega_rad, Omega_b, Omega_cdm );
      fprintf(fp, ", %e", omega_df);

      fprintf(fp, ", %e", particleHorizonVec[i]);

      fprintf(fp, ", %e", mg_field_bk[i]);
      fprintf(fp, ", %e", mg_field_p_bk[i]);

      if(TEST_MODE==1){

        double H_test;

        H_test = a_p[i]/a[i];
        if(T_CONF==0) H_test = H_test / a[i];

        fprintf(fp, ", %e", H_test / sqrt(2. * fourpiG / 3. / (1. + mg_coupl(mg_field_bk[i_a0])) ));

      }

      fprintf(fp, "\n");

      if(a[i+csv_resolution]<= a_end){
        i=i+csv_resolution;
      }
      else{
        i++;
      }

    }

    fclose(fp);

}

// This divides an interval in a logarithmic scale of basis 10, and store results in input vector v
void logscale10 (double * v, double A, double B, int points){

    int i;

    double a = log10(A);
    double b = log10(B);

    double h = (b - a) / (points-1.0);

    for(i=0;i<points;i++){
        v[i] = pow(10.0, a + i * h);
    }

}

// This function is a temporal step of the Runge-Kutta 4 method:
// evolves the system for a time equal to dtau (h in the program)
void mg_rungekutta4bg(double * f, const double dtau)
{

  double a = f[1];
  double a_p = f[2];
  double mg_field_bk = f[3];
  double mg_field_p_bk = f[4];

	double k1a, k2a, k3a, k4a;
	double k1ap, k2ap, k3ap, k4ap;
	double k1f, k2f, k3f, k4f;
	double k1fp, k2fp, k3fp, k4fp;

	k1a 	= a_p_rk4             (a, a_p, mg_field_bk, mg_field_p_bk);
	k1ap 	= a_pp_rk4            (a, a_p, mg_field_bk, mg_field_p_bk);
	k1f 	= mg_field_bk_p_rk4   (a, a_p, mg_field_bk, mg_field_p_bk);
	k1fp 	= mg_field_bk_pp_rk4  (a, a_p, mg_field_bk, mg_field_p_bk);

	k2a 	= a_p_rk4             (a + k1a * dtau / 2., a_p + k1ap * dtau / 2., mg_field_bk + k1f * dtau / 2., mg_field_p_bk + k1fp * dtau / 2.);
	k2ap 	= a_pp_rk4            (a + k1a * dtau / 2., a_p + k1ap * dtau / 2., mg_field_bk + k1f * dtau / 2., mg_field_p_bk + k1fp * dtau / 2.);
	k2f 	= mg_field_bk_p_rk4   (a + k1a * dtau / 2., a_p + k1ap * dtau / 2., mg_field_bk + k1f * dtau / 2., mg_field_p_bk + k1fp * dtau / 2.);
	k2fp 	= mg_field_bk_pp_rk4  (a + k1a * dtau / 2., a_p + k1ap * dtau / 2., mg_field_bk + k1f * dtau / 2., mg_field_p_bk + k1fp * dtau / 2.);

	k3a 	= a_p_rk4             (a + k2a * dtau / 2., a_p + k2ap * dtau / 2., mg_field_bk + k2f * dtau / 2., mg_field_p_bk + k2fp * dtau / 2.);
	k3ap 	= a_pp_rk4            (a + k2a * dtau / 2., a_p + k2ap * dtau / 2., mg_field_bk + k2f * dtau / 2., mg_field_p_bk + k2fp * dtau / 2.);
	k3f 	= mg_field_bk_p_rk4   (a + k2a * dtau / 2., a_p + k2ap * dtau / 2., mg_field_bk + k2f * dtau / 2., mg_field_p_bk + k2fp * dtau / 2.);
	k3fp 	= mg_field_bk_pp_rk4  (a + k2a * dtau / 2., a_p + k2ap * dtau / 2., mg_field_bk + k2f * dtau / 2., mg_field_p_bk + k2fp * dtau / 2.);

	k4a 	= a_p_rk4             (a + k3a * dtau, a_p + k3ap * dtau, mg_field_bk + k3f * dtau, mg_field_p_bk + k3fp * dtau);
	k4ap 	= a_pp_rk4            (a + k3a * dtau, a_p + k3ap * dtau, mg_field_bk + k3f * dtau, mg_field_p_bk + k3fp * dtau);
	k4f 	= mg_field_bk_p_rk4   (a + k3a * dtau, a_p + k3ap * dtau, mg_field_bk + k3f * dtau, mg_field_p_bk + k3fp * dtau);
	k4fp 	= mg_field_bk_pp_rk4  (a + k3a * dtau, a_p + k3ap * dtau, mg_field_bk + k3f * dtau, mg_field_p_bk + k3fp * dtau);

	f[1] += dtau * (k1a + 2. * k2a + 2. * k3a + k4a) / 6.;
	f[2] += dtau * (k1ap + 2. * k2ap + 2. * k3ap + k4ap) / 6.;
	f[3] += dtau * (k1f + 2. * k2f + 2. * k3f + k4f) / 6.;
	f[4] += dtau * (k1fp + 2. * k2fp + 2. * k3fp + k4fp) / 6.;

}

// This function evolves the system with the Runge-Kutta 4 method until t_stop
double rk4(double * t, double * a, double * a_p, double * mg_field_bk, double * mg_field_p_bk, double Omega_f_0) {

    double f[5];
    int i = 0;

    // Call the function for a logarithmic scale of time
    // We don't start from t = 0 to avoid problems with quintessence potentials
    logscale10(t,min_tau,max_tau,points);

    while( i < points - 1 ){

      f[0] = t[i];
      f[1] = a[i];
      f[2] = a_p[i];
      f[3] = mg_field_bk[i];
      f[4] = mg_field_p_bk[i];

      mg_rungekutta4bg(f, t[i+1]-t[i]);

      a[i+1] = f[1];
      a_p[i+1] = f[2];
      mg_field_bk[i+1] = f[3];
      mg_field_p_bk[i+1] = f[4];
      i++;

    }

    i = scan_for_a0(a,1.);

    double H         = Hconf(a[i], mg_field_bk[i], mg_field_p_bk[i]);

    double rho_df = mg_pot(mg_field_bk[i]) + mg_field_p_bk[i] * mg_field_p_bk[i]/2./a[i]/a[i] - 3./2./fourpiG * ( H * mg_field_p_bk[i] * mg_coupl_p(mg_field_bk[i]) + H * H * mg_coupl(mg_field_bk[i]))/a[i]/a[i];
    double Omega_df  = 2. * fourpiG / 3. * a[i] * a[i] * rho_df /H /H;
    return Omega_df  - Omega_f_0;

}

double bisection (double min, double max, double * t, double * a, double * a_p, double * mg_field_bk, double * mg_field_p_bk, const double Omega_f_0, int type) {

  double result = 0.;

  double C = (max+min)/2.;

  if(type==1){

    printf("\t-> Bisection target: mg_field_bk.\n", result);
    if(TEST_MODE==1) printf("\n");

    while(fabs((max-min)/min)>delta_bisection){

      a_p[0] = a_init * Hconf(a[0], min, mg_field_p_bk[0]);
      mg_field_bk[0] = min;
      double rk4_min = rk4(t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0);
      a_p[0] = a_init * Hconf(a[0], C, mg_field_p_bk[0]);
      mg_field_bk[0] = C;
      double rk4_C = rk4(t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0);

      if(rk4_min*rk4_C>=0) {
        min=C;
      }
      else {
        max=C;
      }

      C=(max+min)/2.;

      if (TEST_MODE == 1) printf("TEST_MODE ON - min: %e , max: %e, C: %e, rk4_min: %e , rk4_C: %e \n", min, max, C, rk4_min, rk4_C);

    }

    result = (max+min)/2.;

    if(TEST_MODE==1) printf("\n");
    printf("\t-> Result of bisection method is mg_field_bk: %e (internal units).\n", result);

  }

  else if(type==2){

    printf("\t-> Bisection target: mg_field_p_bk.\n", result);
    if(TEST_MODE==1) printf("\n");

    while(fabs((max-min)/min)>delta_bisection){

      a_p[0] = a_init * Hconf(a[0], mg_field_bk[0], min);
      mg_field_p_bk[0] = min;
      double rk4_min = rk4(t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0);
      a_p[0] = a_init * Hconf(a[0], mg_field_bk[0], C);
      mg_field_p_bk[0] = C;
      double rk4_C = rk4(t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0);

      if(rk4_min*rk4_C>=0) {
        min=C;
      }
      else {
        max=C;
      }

      C=(max+min)/2.;

      if (TEST_MODE == 1) printf("TEST_MODE ON - min: %e , max: %e, C: %e, rk4_min: %e , rk4_C: %e \n", min, max, C, rk4_min, rk4_C);

    }

    result = (max+min)/2.;

    if(TEST_MODE==1) printf("\n");
    printf("\t-> Result of bisection method is mg_field_p_bk: %e (internal units).\n", result);

  }

  else if(type==3){

    printf("\t-> Bisection target: mg_pot_const.\n", result);
    if(TEST_MODE==1) printf("\n");

    while(fabs((max-min)/min)>delta_bisection){

      mg_pot_const = min;
      a_p[0] = a_init * Hconf(a[0], mg_field_bk[0], mg_field_p_bk[0]);
      double rk4_min = rk4(t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0);
      mg_pot_const = C;
      a_p[0] = a_init * Hconf(a[0], mg_field_bk[0], mg_field_p_bk[0]);
      double rk4_C = rk4(t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0);

      if(rk4_min*rk4_C>=0) {
        min=C;
      }
      else {
        max=C;
      }

      C=(max+min)/2.;

      if (TEST_MODE == 1) printf("TEST_MODE ON - min: %e , max: %e, C: %e, rk4_min: %e , rk4_C: %e \n", min, max, C, rk4_min, rk4_C);

    }

    result = (max+min)/2.;

    if(TEST_MODE==1) printf("\n");
    printf("\t-> Result of bisection method is mg_pot_const: %e (internal units).\n", result);

  }

  else{
    result = 0.;
  }

  return result;

}

int main() {

    double Omega_f_0 = 1. - Omega_Lambda_0 - Omega_rad_0 - Omega_b_0 - Omega_cdm_0;

    printf("\n\t\t----------------------------------------\n\n");

    // Definition of the vector needed for the evolution functions
    double * t; double * a; double * a_p; double * mg_field_bk; double * mg_field_p_bk; double * particleHorizonVec;
    t                   = (double *) malloc(sizeof(double) * points);
    a                   = (double *) malloc(sizeof(double) * points);
    a_p                 = (double *) malloc(sizeof(double) * points);
    mg_field_bk         = (double *) malloc(sizeof(double) * points);
    mg_field_p_bk       = (double *) malloc(sizeof(double) * points);
    particleHorizonVec  = (double *) malloc(sizeof(double) * points);

    if(!t||!a||!a_p||!mg_field_bk||!mg_field_p_bk||!particleHorizonVec){
      printf("Error! The memory cannot be allocated. The program will be terminated.\n");
      exit(1);
    }

    // Initial conditions (tau=0)
    a[0] = a_init;

    mg_field_bk[0] = mg_field_bk_0;
    mg_field_p_bk[0] = mg_field_p_bk_0;

    printf("\n Extended quintessence evolution\n");

    printf("\n Parameters used: \n");
    printf("\t- Potential constant: %e\n", mg_pot_const);
    printf("\t- Potential exponent: %e\n", mg_pot_exp);
    printf("\t- Function constant: %e\n", mg_coupl_const);

    printf(" Searching for best initial value for the dark fluid ... \n \n");

    double bisection_result = bisection(bisection_min, bisection_max, t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0,bisection_type);

    if (bisection_type==1) {
      mg_field_bk[0] = bisection_result;
      if(mg_field_bk[0] > bisection_max*(1.-0.05)) printf("\n\n WARNING: it seems that the initial value is too near the bounds! \n\n");
    }
    else if (bisection_type==2) {
      mg_field_p_bk[0] = bisection_result;
      if(mg_field_p_bk[0] > bisection_max*(1.-0.05)) printf("\n\n WARNING: it seems that the initial value is too near the bounds! \n\n");
    }
    else if (bisection_type==3) {
      mg_pot_const = bisection_result;
      if(mg_pot_const > bisection_max*(1.-0.05)) printf("\n\n WARNING: it seems that the initial value is too near the bounds! \n\n");
    }
    else {
      printf("No bisection type chosen. The evolution will be computed with provided initial conditions.");
    }

    a_p[0] = a_init * Hconf(a_init, mg_field_bk[0], mg_field_p_bk[0]);

    printf("\n\n Evolving the system ...\n");

    rk4(t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0);

    int last_int = scan_for_a0(a,a_end);
    int last_int_print = scan_for_a0(a,1.);

    if(last_int >= points-1) printf("\n\n WARNING: the a = a_end value seems not to be reached in the result vectors! \n\n");

    double t_age = 0.;
    int i = 0;
    for(i=2;i<last_int_print;i++){
      t_age += (a[i]+a[i-1])*(t[i]-t[i-1])/2.;
    }

    printf("\n Results:\n");
    printf("\t- H0: %f \n", a_p[last_int_print]/a[last_int_print]/sqrt(2. * fourpiG / 3.) );
    printf("\t- Age of the Universe: %f Gyr\n", sqrt(2. * fourpiG / 3.) * t_age/_H0_/(60.*60.*24.*365.*1e9)*_MPc_over_m_/1000.);//t[last_int_print]/((4*M_PI)/3.)
    if(TEST_MODE==1) printf("\t- Number of points: %d \n", last_int);

    char filename[50];
    sprintf (filename, "coupl_mg_bk.csv");

    csv(t, a, a_p, mg_field_bk, mg_field_p_bk, particleHorizonVec, filename);

    printf("\n The results are saved in '.csv' files.\n");

    if(TEST_MODE==1) printf("\n TEST_MODE ON: check the values of H in .csv file. They must be equal!");

    printf("\n\t\t----------------------------------------\n");

    free(a);free(a_p);free(mg_field_bk);free(mg_field_p_bk);free(particleHorizonVec);

    printf("\n");

    exit(0);

}

ext: c

max_stars_repo_name: alessandrocasalino/mgevolution

max_issues_repo_name: alessandrocasalino/mgevolution

max_forks_repo_name: alessandrocasalino/mgevolution


text:
#include <gsl/gsl_sf_gamma.h>

//
// TIMING
//

// Clock functions
static vector<double> ClockTimes;
static double C0;

double Clock()  { return double(clock()) / CLOCKS_PER_SEC; }

void StartClocking()
{
    ClockTimes.clear();
    C0 = Clock();
}

void ClockCheckpoint(unsigned int cp)
{
    double C1 = Clock();
    if (cp >= ClockTimes.size()) ClockTimes.resize(cp + 1, 0.0);
    ClockTimes[cp] += C1 - C0;
    C0 = C1;
}

void ShowClockInfo()
{
    double total = 0;
    for (unsigned int i = 0; i < ClockTimes.size(); ++i)
    {
        cout << "Checkpoint " << i << ": " << ClockTimes[i] << "\n";
        total += ClockTimes[i];
    }
    cout << "Total time: " << total << " seconds.\n";
}

//
// DISCRETE DISTRIBUTION
//

struct Parameters;

class MNApprox
{
public:
    static const unsigned int NVariants = 32;
    static const unsigned int VariantMask = 31;

    void Set(Parameters& P, Randomizer& Rand, vector<double>& p);

    void operator()(unsigned int N, vector<unsigned int>& out)
    {
        out.assign(out.size(), 0);
        for (unsigned int shift = 0; shift < 32; ++shift)
        {
            if ((N >> shift) & 1)
            {
                for (unsigned int i = 0; i < out.size(); ++i)
                {
                    out[i] += x[shift][cycle[shift]][i];
                }
                cycle[shift] = (cycle[shift] + 1) & VariantMask;
            }
        }
    }

private:
    vector<vector<vector<unsigned int>>> x;
    vector<unsigned int> cycle;
};

// An arbitrary discrete distribution, used for maturation times in compartments
class Discrete
{
public:
    // Create distribution from set of unnormalised weights
    void operator=(std::vector<double> uw)
    {
        weights = uw;

        double total = accumulate(weights.begin(), weights.end(), 0.0);
        for (auto& w : weights)
            w /= total;

        storage.assign(weights.size(), 0);
    }

    // Normalised weights
    std::vector<double> weights;

    // Storage for draws of integers
    std::vector<unsigned int> storage;

    // Fast approximate multinomial draws
    MNApprox mn_approx;
};

//
// MATRIX
//

struct Matrix
{
    Matrix()
     : nc(0)
    { }

    Matrix(double X, unsigned int nrow, unsigned int ncol)
     : x(nrow * ncol, X), nc(ncol)
    { }

    Matrix(vector<double> X, unsigned int nrow, unsigned int ncol)
     : x(X), nc(ncol)
    {
        if (x.size() != nrow * ncol)
            throw runtime_error("Improperly sized matrix.");
    }

    double& operator()(unsigned int i, unsigned int j)
    {
        return x[i * nc + j];
    }

    unsigned int NCol() { return nc; }
    unsigned int NRow() { return x.size() / nc; }

    vector<double> x;
    unsigned int nc;
};

void MNApprox::Set(Parameters& P, Randomizer& Rand, vector<double>& p)
{
    x.assign(32, vector<vector<unsigned int>>(NVariants, vector<unsigned int>(p.size(), 0)));
    cycle.assign(NVariants, 0);
    for (unsigned int shift = 0; shift < 32; ++shift)
        for (unsigned int variant = 0; variant < NVariants; ++variant)
            gsl_ran_multinomial(Rand.GSL_RNG(), p.size(), 1 << shift, &p[0], &x[shift][variant][0]);
}

max_stars_repo_name: nicholasdavies/covid-tiers


text:
/* randist/chisq.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 <math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* The chisq distribution has the form

   p(x) dx = (1/(2*Gamma(nu/2))) (x/2)^(nu/2 - 1) exp(-x/2) dx

   for x = 0 ... +infty */

double
gsl_ran_chisq (const gsl_rng * r, const double nu)
{
  double chisq = 2 * gsl_ran_gamma (r, nu / 2, 1.0);
  return chisq;
}

double
gsl_ran_chisq_pdf (const double x, const double nu)
{
  if (x < 0)
    {
      return 0 ;
    }
  else
    {
      if(nu == 2.0)
        {
          return exp(-x/2.0) / 2.0;
        }
      else
        {
          double p;
          double lngamma = gsl_sf_lngamma (nu / 2);

          p = exp ((nu / 2 - 1) * log (x/2) - x/2 - lngamma) / 2;
          return p;
        }
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include "rb_sm.h"
#include <stdio.h>
#include <gsl/gsl_nan.h>
#include <options/options.h>

struct sm_params rb_sm_params;
struct sm_result rb_sm_result;

void rb_sm_init_journal(const char*journal_file) {
	jj_set_stream(open_file_for_writing(journal_file));
/*	sm_journal_open(journal_file);*/
}

void rb_sm_close_journal() {
	FILE * s = jj_get_stream();
	if(s) fclose(s);

	jj_set_stream(0);
}

void rb_sm_odometry(double x, double y, double theta){
	rb_sm_params.first_guess[0]=x;
	rb_sm_params.first_guess[1]=y;
	rb_sm_params.first_guess[2]=theta;
}

void rb_sm_odometry_cov(double cov_x, double cov_y, double cov_theta){

}

struct option* ops = 0;

int rb_sm_set_configuration(const char*name, const char*value) {
	if(!ops) {
		ops = options_allocate(30);
		sm_options(&rb_sm_params, ops);
	}

	if(!options_try_pair(ops, name, value)) {

		return 0;
	} else
	return 1;
}

const char *rb_result_to_json() {
	static char buf[5000];
	JO jo = result_to_json(&rb_sm_params, &rb_sm_result);
	strcpy(buf, jo_to_string(jo));
	jo_free(jo);
	return buf;
}

int rb_sm_icp() {
	sm_icp(&rb_sm_params, &rb_sm_result);
	return rb_sm_result.valid;
}

int rb_sm_gpm() {
	sm_gpm(&rb_sm_params, &rb_sm_result);
	return rb_sm_result.valid;
}

void rb_set_laser_ref(const char*s) {
	rb_sm_params.laser_ref = string_to_ld(s);
/*	fprintf(stderr, "Set laser_ref to %p\n ", rb_sm_params.laser_ref );*/
}

void rb_set_laser_sens(const char*s) {
	rb_sm_params.laser_sens = string_to_ld(s);
/*	fprintf(stderr, "Set laser_sens to %p\n ", rb_sm_params.laser_sens );*/
}

void rb_sm_cleanup() {
	if(rb_sm_params.laser_ref)
	ld_free(rb_sm_params.laser_ref);
	if(rb_sm_params.laser_sens)
	ld_free(rb_sm_params.laser_sens);
}

LDP string_to_ld(const char*s) {
	JO jo = json_parse(s);
	if(!jo) {
		fprintf(stderr, "String passed from Ruby is invalid JSON: \n\n%s\n", s);
		return 0;
	}
	LDP ld = json_to_ld(jo);
	if(!ld) {
		fprintf(stderr, "String passed from Ruby is valid JSON, "
			"but can't load laser_data. \n\n%s\n", s);
		return 0;
	}
	jo_free(jo);
	return ld;
}


text:
#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>

/* Compile all the inline matrix functions */

#define COMPILE_INLINE_STATIC
#include "build.h"
#include <gsl/gsl_matrix.h>

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

	1. Redistributions of source code must retain the above copyright notice,
		this list of conditions and the following disclaimer.

	2. 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.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``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 Javier Burguete Tolosa OR CONTRIBUTORS 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.
*/

/**
 * \file rk_6_2.c
 * \brief Source file to optimize Runge-Kutta 6 steps 2nd order methods.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <string.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#include "config.h"
#include "utils.h"
#include "optimize.h"
#include "rk.h"
#include "rk_6_2.h"

#define DEBUG_RK_6_2 0          ///< macro to debug.

/**
 * Function to obtain the coefficients of a 6 steps 2nd order Runge-Kutta
 * method.
 */
int
rk_tb_6_2 (Optimize * optimize) ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_6_2
  fprintf (stderr, "rk_tb_6_2: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t6 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b31 (tb) = r[4];
  b32 (tb) = r[5];
  t4 (tb) = r[6];
  b41 (tb) = r[7];
  b42 (tb) = r[8];
  b43 (tb) = r[9];
  t5 (tb) = r[10];
  b51 (tb) = r[11];
  b52 (tb) = r[12];
  b53 (tb) = r[13];
  b54 (tb) = r[14];
  b62 (tb) = r[15];
  b63 (tb) = r[16];
  b64 (tb) = r[17];
  b65 (tb) = r[18];
  b61 (tb) = (0.5L - b62 (tb) * t2 (tb) - b63 (tb) * t3 (tb)
              - b64 (tb) * t4 (tb) - b65 (tb) * t5 (tb)) / t1 (tb);
  if (isnan (b61 (tb)))
    return 0;
  rk_b_6 (tb);
#if DEBUG_RK_6_2
  rk_print_tb (optimize, "rk_tb_6_2", stderr);
  fprintf (stderr, "rk_tb_6_2: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 6 steps 2nd order, 3rd order in
 * equations depending only in time, Runge-Kutta method.
 */
int
rk_tb_6_2t (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_6_2
  fprintf (stderr, "rk_tb_6_2t: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t6 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b31 (tb) = r[4];
  b32 (tb) = r[5];
  t4 (tb) = r[6];
  b41 (tb) = r[7];
  b42 (tb) = r[8];
  b43 (tb) = r[9];
  t5 (tb) = r[10];
  b51 (tb) = r[11];
  b52 (tb) = r[12];
  b53 (tb) = r[13];
  b54 (tb) = r[14];
  b61 (tb) = r[15];
  b62 (tb) = r[16];
  b63 (tb) = r[17];
  b64 (tb) = (1.L / 3.L - 0.5L * t5 (tb)
              - b61 (tb) * t1 (tb) * (t1 (tb) - t5 (tb))
              - b62 (tb) * t2 (tb) * (t2 (tb) - t5 (tb))
              - b63 (tb) * t3 (tb) * (t3 (tb) - t5 (tb)))
    / (t4 (tb) * (t4 (tb) - t5 (tb)));
  if (isnan (b64 (tb)))
    return 0;
  b65 (tb) = (0.5L - b61 (tb) * t1 (tb) - b62 (tb) * t2 (tb)
              - b63 (tb) * t3 (tb) - b64 (tb) * t4 (tb)) / t5 (tb);
  if (isnan (b65 (tb)))
    return 0;
  rk_b_6 (tb);
#if DEBUG_RK_6_2
  rk_print_tb (optimize, "rk_tb_6_2t", stderr);
  fprintf (stderr, "rk_tb_6_2t: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 6 steps 1st-2nd order Runge-Kutta
 * pair.
 */
int
rk_tb_6_2p (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_6_2
  fprintf (stderr, "rk_tb_6_2p: start\n");
#endif
  if (!rk_tb_6_2 (optimize))
    return 0;
  tb = optimize->coefficient;
  e61 (tb) = e62 (tb) = e63 (tb) = e64 (tb) = 0.L;
  rk_e_6 (tb);
#if DEBUG_RK_6_2
  fprintf (stderr, "rk_tb_6_2p: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 6 steps 1st-2nd order, 1st-3rd order
 * in equations depending only in time, Runge-Kutta pair.
 */
int
rk_tb_6_2tp (Optimize * optimize)       ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_6_2
  fprintf (stderr, "rk_tb_6_2tp: start\n");
#endif
  if (!rk_tb_6_2t (optimize))
    return 0;
  tb = optimize->coefficient;
  e61 (tb) = e62 (tb) = e63 (tb) = e64 (tb) = 0.L;
  rk_e_6 (tb);
#if DEBUG_RK_6_2
  fprintf (stderr, "rk_tb_6_2tp: end\n");
#endif
  return 1;
}

/**
 * Function to calculate the objective function of a 6 steps 2nd order
 * Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_6_2 (RK * rk)   ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_6_2
  fprintf (stderr, "rk_objective_tb_6_2: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b60 (tb) < 0.L)
    o += b60 (tb);
  if (b61 (tb) < 0.L)
    o += b61 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L,
             fmaxl (t1 (tb),
                    fmaxl (t2 (tb),
                           fmaxl (t3 (tb), fmaxl (t4 (tb), t5 (tb))))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_6_2
  fprintf (stderr, "rk_objective_tb_6_2: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_6_2: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 6 steps 2nd order, third
 * order in equations depending only on time, Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_6_2t (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_6_2
  fprintf (stderr, "rk_objective_tb_6_2t: start\n");
#endif
  tb = rk->tb->coefficient;
#if DEBUG_RK_6_2
  rk_print_tb (optimize, "rk_objective_tb_6_2t", stderr);
#endif
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b60 (tb) < 0.L)
    o += b60 (tb);
  if (b64 (tb) < 0.L)
    o += b64 (tb);
  if (b65 (tb) < 0.L)
    o += b65 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L,
             fmaxl (t1 (tb),
                    fmaxl (t2 (tb),
                           fmaxl (t3 (tb), fmaxl (t4 (tb), t5 (tb))))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_6_2
  fprintf (stderr, "rk_objective_tb_6_2t: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_6_2t: end\n");
#endif
  return o;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
/*
 * optimize_tnc.h
 *
 *  Created on: Feb 9, 2010
 *      Author: smitty
 */

#ifndef _OPTIMIZE_CONT_MODELS_NLOPT_H_
#define _OPTIMIZE_CONT_MODELS_NLOPT_H_

#include <nlopt.h>

#include "cont_models.h"

#include <armadillo>
using namespace arma;

double nlopt_bm_sr(unsigned n, const double *x, double *grad, void *data);
double nlopt_bm_sr_log(unsigned n, const double *x, double *grad, void *data);
double nlopt_ou_sr_log(unsigned n, const double *x, double *grad, void *data);
double nlopt_bm_bl(unsigned n, const double *x, double *grad, void *data);
vector<double> optimize_single_rate_bm_nlopt(rowvec & _x, mat & _vcv,bool log);
vector<double> optimize_single_rate_bm_ou_nlopt(rowvec & _x, mat & _vcv);
vector<double> optimize_single_rate_bm_bl(Tree * tr);

#endif /* _OPTIMIZE_CONT_MODELS_NLOPT_H_ */


text:
#define USE_GSL_INTEGRATION 0

#if USE_GSL_INTEGRATION
	#include <gsl/gsl_integration.h>
#endif

struct integrand_params {
	u32 n[2];
	u32 component_count;
	u32 coeff_count;
	f64* coeff;

	nlse_operator_func* op;
	void* op_userdata;

	struct basis basis;
};

static void sbmf_log_integration_result(struct quadgk_result* res) {
	sbmf_log_info("integral: %.10e", res->integral);
	sbmf_log_info("error: %.10e", res->error);
	//sbmf_log_info("performed evals: %d", res->performed_evals);
	sbmf_log_info("converged: %s", (res->converged) ? "yes" : "no");
}

/* functions to handle spatial guesses */

struct guess_integrand_params {
    u32 n;
	nlse_spatial_guess_func* func;
	struct basis b;
};

void guess_integrand(f64* out, f64* in, u32 len, void* data) {
    struct guess_integrand_params* params = data;

    f64 eig[len];
	params->b.eigenfunc(params->n, len, eig, in);

	f64 sample[len];
	params->func(sample, in, len, NULL);

    for (u32 i = 0; i < len; ++i) {
        out[i] = eig[i] * sample[i];
    }
}

/* function to sample linear and non-linear matrix elements */
static void linear_me_integrand(f64* out, f64* in, u32 len, void* data) {
	struct integrand_params* params = data;

	f64 eig1[len];
	f64 eig2[len];
	params->basis.eigenfunc(params->n[0], len, eig1, in);
	params->basis.eigenfunc(params->n[1], len, eig2, in);

	f64 op[len];
	params->op(len, op, in, 0, NULL, params->op_userdata);

	for (u32 i = 0; i < len; ++i) {
		out[i] = eig1[i]*eig2[i]*op[i];
	}
}

static void nonlinear_me_integrand(f64* out, f64* in, u32 len, void* data) {
	struct integrand_params* params = data;

	f64 eig1[len];
	f64 eig2[len];
	params->basis.eigenfunc(params->n[0], len, eig1, in);
	params->basis.eigenfunc(params->n[1], len, eig2, in);

	f64 sample[len*params->component_count];
	for (u32 i = 0; i < params->component_count; ++i) {
		params->basis.sample(params->coeff_count, &params->coeff[i*params->coeff_count], len, &sample[i*len], in);
	}

	f64 op[len];
	params->op(len, op, in, params->component_count, sample, params->op_userdata);

	for (u32 i = 0; i < len; ++i) {
		out[i] = eig1[i]*eig2[i]*op[i];
	}
}

#if USE_GSL_INTEGRATION
static f64 nonlinear_me_integrand_gsl(f64 in, void* data) {
	struct integrand_params* params = data;

	f64 eig1;
	f64 eig2;
	params->basis.eigenfunc(params->n[0], 1, &eig1, &in);
	params->basis.eigenfunc(params->n[1], 1, &eig2, &in);

	f64 sample[params->component_count];
	for (u32 i = 0; i < params->component_count; ++i) {
		params->basis.sample(params->coeff_count, &params->coeff[i*params->coeff_count], 1, &sample[i], &in);
	}

	f64 op;
	params->op(1, &op, &in, params->component_count, sample, params->op_userdata);

	return eig1*eig2*op;
}
#endif

struct nlse_result nlse_solver(struct nlse_settings settings, const u32 component_count, struct nlse_component* component) {
	/* Lazy */
	const u32 N = settings.num_basis_funcs;

	/* Setup results struct */
	struct nlse_result res = {
		.iterations = 0,
		.component_count = component_count,
		.coeff_count = N,
		.coeff 		 = sbmf_stack_push(component_count*(N*sizeof(f64))),
		.abs_error 	 = sbmf_stack_push(component_count*sizeof(f64)),
		.energy 	 = sbmf_stack_push(component_count*sizeof(f64)),
		.residual 	 = sbmf_stack_push(component_count*sizeof(f64)),
		.hamiltonian = sbmf_stack_push(component_count * sizeof(struct symmetric_bandmat)),
		.converged = false,
	};
	for (u32 i = 0; i < component_count; ++i) {
		res.hamiltonian[i] = symmetric_bandmat_new(N,N);
	}

	/* Place to store coeffs from previous iterations */
	f64* old_coeff = sbmf_stack_push(component_count*(N*sizeof(f64)));
	f64* old_energy = sbmf_stack_push(component_count * sizeof(f64));

	struct quadgk_settings int_settings = {
		.gk = (settings.gk.gauss_size > 0) ? settings.gk : gk15,
		.abs_error_tol = 1e-15,
		.rel_error_tol = 1e-8,
		.max_iters = settings.max_quadgk_iters,
	};

	/* Setup intial guess values for coeffs */
	for (u32 i = 0; i < component_count; ++i) {
		switch (component[i].guess.type) {
			case DEFAULT_GUESS: {
				for (u32 j = 0; j < res.coeff_count; ++j) {
					res.coeff[i*res.coeff_count + j] = 0;
				}
				/* Initialize the i:th component to the i:th eigenfunction */
				res.coeff[i*res.coeff_count + i] = 1;
				break;
			}
			case SPATIAL_GUESS: {
				struct guess_integrand_params p = {
					.func = component[i].guess.data.spatial_guess,
					.b = settings.basis,
				};
				int_settings.userdata = &p;
				f64* out = &res.coeff[i*res.coeff_count];
				for (u32 j = 0; j < res.coeff_count; ++j) {
					p.n = j;
					struct quadgk_result r;

					u8 quadgk_memory[quadgk_required_memory_size(&int_settings)];
					quadgk_infinite_interval(guess_integrand, &int_settings, quadgk_memory, &r);
					out[j] = r.integral;
				}
				f64_normalize(out, out, res.coeff_count);

				break;
			}
			case COEFF_GUESS: {
				component[i].guess.data.coeff_guess(&res.coeff[i*res.coeff_count], res.coeff_count, i);
				break;
			}
			case RANDOM_GUESS: {
				struct symmetric_bandmat bm = symmetric_bandmat_new_zero(N,N);
				SYMMETRIC_BANDMAT_FOREACH(bm, r,c) {
					u32 index = symmetric_bandmat_index(bm, r,c);
					bm.data[index] = 2.0 * ((f64)rand()/(f64)RAND_MAX) - 1.0;
				}
				struct eigen_result_real eigres = find_eigenpairs_sparse_real(bm, 1, EV_SMALLEST);
				for (u32 j = 0; j < res.coeff_count; ++j) {
					res.coeff[i*res.coeff_count + j] = eigres.eigenvectors[j];
				}
				break;
			}
			default:
				sbmf_log_error("nlse_solver: component %u has invalid guess!", i);
				return res;
		};
	}

	struct integrand_params params = {
		.component_count = res.component_count,
		.coeff_count = res.coeff_count,
		.coeff = res.coeff,
		.op = NULL,
		.basis = settings.basis,
	};

	/*
	 * Unique number of me's that need to be calculated in
	 * a symmetric matrix
	 */
	const u32 matrix_element_count = symmetric_bandmat_element_count(N);

	/* Precompute indices for matrix elements
	 * 	This way we get a single array looping
	 * 	over the matrix which is easier to
	 * 	parallelize well.
	 * 		Picking the first hamiltonian as a
	 * 	dummny matrix just to get the correct
	 * 	iteration.
	 * */
	u32 matrix_element_rows[matrix_element_count];
	u32 matrix_element_cols[matrix_element_count];
	{
		u32 matrix_element_index = 0;
		SYMMETRIC_BANDMAT_FOREACH(res.hamiltonian[0], r,c) {
			matrix_element_rows[matrix_element_index] = r;
			matrix_element_cols[matrix_element_index] = c;
			matrix_element_index++;
		}
	}

	/* Construct linear hamiltonian */
	struct symmetric_bandmat linear_hamiltonian = symmetric_bandmat_new_zero(N,N);
	{
		if (settings.spatial_pot_perturbation) {
			params.op = settings.spatial_pot_perturbation;
#pragma omp parallel for firstprivate(params, int_settings) shared(res)
			for (u32 i = 0; i < matrix_element_count; ++i) {
				u32 r = matrix_element_rows[i];
				u32 c = matrix_element_cols[i];

				int_settings.userdata = &params;
				params.n[0] = r;
				params.n[1] = c;

				u8 quadgk_memory[quadgk_required_memory_size(&int_settings)];
				struct quadgk_result int_res;
				quadgk_infinite_interval(linear_me_integrand, &int_settings, quadgk_memory, &int_res);

				if (fabs(int_res.integral) <= settings.zero_threshold)
					int_res.integral = 0.0;

				if (!int_res.converged) {
					sbmf_log_error("In construction of linear hamiltonian:");
					sbmf_log_error("\tIntegration failed for %d,%d", r,c);
					sbmf_log_integration_result(&int_res);
				}
				assert(int_res.converged);

				u32 me_index = symmetric_bandmat_index(linear_hamiltonian, r,c);
				linear_hamiltonian.data[me_index] = int_res.integral;
			}
		}

		for (u32 i = 0; i < res.coeff_count; ++i) {
			u32 me_index = symmetric_bandmat_index(linear_hamiltonian, i,i);
			linear_hamiltonian.data[me_index] += settings.basis.eigenval(i);
		}

		assert(symmetric_bandmat_is_valid(linear_hamiltonian));
	}

#if USE_GSL_INTEGRATION
	const u32 num_threads = omp_get_max_threads();
	gsl_integration_workspace* ws[num_threads];
	/* Create gsl workspaces if necessary */
	for (u32 i = 0; i < num_threads; ++i) {
		ws[i] = gsl_integration_workspace_alloc(settings.max_iterations);
	}
#endif

	struct symmetric_bandmat old_Hs[res.component_count];
	struct symmetric_bandmat premix_Hs[res.component_count];
	for (u32 i = 0; i < res.component_count; ++i) {
		old_Hs[i] = symmetric_bandmat_new(N,N);
		premix_Hs[i] = symmetric_bandmat_new(N,N);
	}

	/* Do the actual iterations */
	for (; res.iterations < settings.max_iterations; ++res.iterations) {
		sbmf_log_info("nlse starting iteration %u", res.iterations);
		memcpy(old_energy, res.energy, res.component_count * sizeof(f64));

		/* Call debug callback if requested by user */
		if (settings.measure_every > 0 &&
				settings.debug_callback &&
				res.iterations % settings.measure_every == 0) {
			settings.debug_callback(settings, res);
		}

		/* Should we still apply mixing? */
		if (settings.orbital_mixing > 0 && settings.mix_until_iteration > 0 && res.iterations == settings.mix_until_iteration) {
			settings.orbital_mixing = 0;
		}
		if (settings.hamiltonian_mixing > 0 && settings.mix_until_iteration > 0 && res.iterations == settings.mix_until_iteration) {
			settings.hamiltonian_mixing = 0;
		}

		if (settings.orbital_choice != NLSE_ORBITAL_MAXIMUM_OVERLAP && settings.mom_enable_at_iteration > 0 && res.iterations == settings.mom_enable_at_iteration) {
			settings.orbital_choice = NLSE_ORBITAL_MAXIMUM_OVERLAP;
		}

		/* Apply orbital mixing */
		for (u32 i = 0; i < component_count; ++i) {
			if (settings.orbital_mixing > 0.0) {
				//res.energy[i] = (1.0 - settings.orbital_mixing) * res.energy[i] + settings.orbital_mixing*old_energy[i];

				for (u32 j = 0; j < res.coeff_count; ++j) {
					res.coeff[i*res.coeff_count + j] =
						(1.0 - settings.orbital_mixing) * res.coeff[i*res.coeff_count + j] + settings.orbital_mixing * old_coeff[i*res.coeff_count + j];
				}
				f64_normalize(&res.coeff[i*res.coeff_count], &res.coeff[i*res.coeff_count], res.coeff_count);
			}
		}
		memcpy(old_coeff, res.coeff, res.component_count * N * sizeof(f64));

		/*
		 * Construct all the Hamiltonians
		 */
		for (u32 i = 0; i < component_count; ++i) {
			params.op = component[i].op;
			params.op_userdata = component[i].userdata;

			memcpy(old_Hs[i].data, res.hamiltonian[i].data, N*N*sizeof(f64));

			/* Construct the i:th component's hamiltonian */
#pragma omp parallel for firstprivate(params, int_settings) shared(res, linear_hamiltonian)
			for (u32 j = 0; j < matrix_element_count; ++j) {
				u32 r = matrix_element_rows[j];
				u32 c = matrix_element_cols[j];
				int_settings.userdata = &params;
				params.n[0] = r;
				params.n[1] = c;

#if USE_GSL_INTEGRATION
				integration_result int_res;
				gsl_function F;
				F.function = nonlinear_me_integrand_gsl;
				F.params = &params;
				gsl_integration_qagi(&F, 1e-10, 1e-7, settings.max_iterations, ws[omp_get_thread_num()], &int_res.integral, &int_res.error);
				int_res.converged = true;
#else
				u8 quadgk_memory[quadgk_required_memory_size(&int_settings)];

				struct quadgk_result int_res;
				quadgk_infinite_interval(nonlinear_me_integrand, &int_settings, quadgk_memory, &int_res);
#endif

				/* Check if the resultant integral is less than what we can resolve,
				 * if so, zero it.
				 */
				if (fabs(int_res.integral) <= settings.zero_threshold)
					int_res.integral = 0.0;

				if (!int_res.converged) {
					sbmf_log_error("In construction of component %u's hamiltonian:", i);
					sbmf_log_error("\tIntegration failed for %d,%d", r,c);
					sbmf_log_integration_result(&int_res);
				}
				assert(int_res.converged);

				u32 me_index = symmetric_bandmat_index(res.hamiltonian[i], r,c);
				res.hamiltonian[i].data[me_index] = linear_hamiltonian.data[me_index] + int_res.integral;
				premix_Hs[i].data[me_index] = res.hamiltonian[i].data[me_index];

				if (settings.hamiltonian_mixing > 0) {
					res.hamiltonian[i].data[me_index] = (1.0 - settings.hamiltonian_mixing)*res.hamiltonian[i].data[me_index] + settings.hamiltonian_mixing*old_Hs[i].data[me_index];
				}
			}

			assert(symmetric_bandmat_is_valid(res.hamiltonian[i]));
		}

		/*
		 * Solve and normalize all the Hamiltonian
		 * eigenvalue problems
		 */
		for (u32 i = 0; i < component_count; ++i) {
			struct eigen_result_real eigres;
			u32 new_orbital_index = 0;
			switch (settings.orbital_choice) {
				case NLSE_ORBITAL_LOWEST_ENERGY: {
					eigres = find_eigenpairs_sparse_real(res.hamiltonian[i], 1, EV_SMALLEST);
					f64_normalize(&eigres.eigenvectors[0], &eigres.eigenvectors[0], res.coeff_count);
					break;
				}
				case NLSE_ORBITAL_MAXIMUM_OVERLAP: {
					eigres = find_eigenpairs_sparse_real(res.hamiltonian[i], settings.mom_orbitals_to_consider, EV_SMALLEST);
					f64 maximum_overlap = -INFINITY;
					for (u32 j = 0; j < settings.mom_orbitals_to_consider; ++j) {
						f64_normalize(&eigres.eigenvectors[j*res.coeff_count], &eigres.eigenvectors[j*res.coeff_count], res.coeff_count);
						f64 overlap = 0.0;
						for (u32 k = 0; k < res.coeff_count; ++k) {
							f64 a = eigres.eigenvectors[j*res.coeff_count + k];
							f64 b = old_coeff[i*res.coeff_count + k];
							overlap += a*b;
						}
						overlap = fabs(overlap);

						if (overlap > maximum_overlap) {
							new_orbital_index = j;
							maximum_overlap = overlap;
						}
					}
					break;
				}
				default: {
					sbmf_log_panic("nlse_solver(...): Unspecified energy selection method!");
					assert(0);
				}
			};

			/* Copy energies and coeffs */
			res.energy[i] = eigres.eigenvalues[new_orbital_index];
			for (u32 j = 0; j < res.coeff_count; ++j) {
				res.coeff[i*res.coeff_count + j] = eigres.eigenvectors[res.coeff_count*new_orbital_index + j];
			}
			f64_normalize(&res.coeff[i*res.coeff_count], &res.coeff[i*res.coeff_count], res.coeff_count);

			///* Apply orbital mixing */
			//if (settings.orbital_mixing > 0.0) {
			//	res.energy[i] = (1.0 - settings.orbital_mixing) * res.energy[i] + settings.orbital_mixing*old_energy[i];

			//	for (u32 j = 0; j < res.coeff_count; ++j) {
			//		res.coeff[i*res.coeff_count + j] =
			//			(1.0 - settings.orbital_mixing) * res.coeff[i*res.coeff_count + j] + settings.orbital_mixing * old_coeff[i*res.coeff_count + j];
			//	}
			//	f64_normalize(&res.coeff[i*res.coeff_count], &res.coeff[i*res.coeff_count], res.coeff_count);
			//}

		}

		/* Calculate error */
		for (u32 i = 0; i < component_count; ++i) {
			f64 sum_diff = 0.0;
			f64 sum_prev = 0.0;
			for (u32 j = 0; j < res.coeff_count; ++j) {
				f64 diff = fabs(res.coeff[i*res.coeff_count + j]) - fabs(old_coeff[i*res.coeff_count + j]);
				sum_diff += diff*diff;
				sum_prev += old_coeff[i*res.coeff_count + j]*old_coeff[i*res.coeff_count + j];
			}
			res.abs_error[i] = sqrt(sum_diff);
		}

		/* Break condition */
		bool should_exit = true;
		for (u32 i = 0; i < component_count; ++i) {
			if (res.abs_error[i] > settings.abs_error_tol)
				should_exit = false;
			sbmf_log_info("\t[%u] -- abs error: %.10e energy: %.10e", i, res.abs_error[i], res.energy[i]);
		}

		if (should_exit) {
			res.converged = true;
			break;
		}
	}

#if USE_GSL_INTEGRATION
	/* free gsl workspaces if necessary */
	for (u32 i = 0; i < num_threads; ++i) {
		gsl_integration_workspace_free(ws[i]);
	}
#endif

	/* Compute residuals */
	for (u32 i = 0; i < component_count; ++i) {
		f64 ans1[N], ans2[N];
		symmetric_bandmat_mulv(ans1, premix_Hs[i], &res.coeff[i*N]);

		for (u32 j = 0; j < N; ++j) {
			ans2[j] = res.energy[i] * res.coeff[i*N + j];
		}

		f64 sum = 0.0;
		for (u32 j = 0; j < N; ++j) {
			sum += fabs(ans1[j] - ans2[j]);
		}
		sum = sqrt(sum);

		res.residual[i] = sum;
		sbmf_log_info("\t[%u] residual %e", i, sum);
	}

	return res;
}

/*
 * basic serialization
 */

void nlse_write_to_binary_file(const char* file, struct nlse_result res) {
	FILE* fd = fopen(file, "w");
	fwrite(&res.iterations, 	 sizeof(u32), 1, fd);
	fwrite(&res.component_count, sizeof(u32), 1, fd);
	fwrite(&res.coeff_count, 	 sizeof(u32), 1, fd);

	fwrite(res.coeff,  sizeof(f64), res.coeff_count*res.component_count, fd);
	fwrite(res.abs_error,  sizeof(f64), res.component_count, fd);
	fwrite(res.energy, sizeof(f64), res.component_count, fd);

	for (u32 i = 0; i < res.component_count; ++i) {
		fwrite(&res.hamiltonian[i].size, 	  sizeof(u32), 1, fd);
		fwrite(&res.hamiltonian[i].bandcount, sizeof(u32), 1, fd);
		u32 elements_written = fwrite(res.hamiltonian[i].data, sizeof(f64),
				res.hamiltonian[i].bandcount*res.hamiltonian[i].size,
				fd);
		assert(elements_written == res.hamiltonian[i].bandcount*res.hamiltonian[i].size);
	}
	fwrite(&res.converged, sizeof(bool), 1, fd);
	fclose(fd);
}

struct nlse_result nlse_read_from_binary_file(const char* file) {
	struct nlse_result res;

	FILE* fd = fopen(file, "r");
	fread(&res.iterations, 	 	sizeof(u32), 1, fd);
	fread(&res.component_count, sizeof(u32), 1, fd);
	fread(&res.coeff_count,  	sizeof(u32), 1, fd);

	res.coeff = sbmf_stack_push(sizeof(f64)*res.coeff_count*res.component_count);
	fread(res.coeff,  sizeof(f64), res.coeff_count*res.component_count, fd);

	res.abs_error = sbmf_stack_push(sizeof(f64)*res.component_count);
	fread(res.abs_error,  sizeof(f64), res.component_count, fd);

	res.energy = sbmf_stack_push(sizeof(f64)*res.component_count);
	fread(res.energy, sizeof(f64), res.component_count, fd);

	res.hamiltonian = sbmf_stack_push(res.component_count * sizeof(struct symmetric_bandmat));
	for (u32 i = 0; i < res.component_count; ++i) {
		u32 size, bandcount;
		fread(&size, 	  sizeof(u32), 1, fd);
		fread(&bandcount, sizeof(u32), 1, fd);
		res.hamiltonian[i] = symmetric_bandmat_new(bandcount, size);
		u32 bytes_written = fread(res.hamiltonian[i].data, sizeof(f64),
				bandcount*size,
				fd);
		assert(bytes_written == bandcount*size);
	}
	fread(&res.converged, sizeof(bool), 1, fd);
	fclose(fd);

	sbmf_log_info("loaded:");
	sbmf_log_info("    iterations: %u", res.iterations);
	sbmf_log_info("    component_count: %u", res.component_count);
	sbmf_log_info("    coeff_count: %u", res.coeff_count);
	for (u32 i = 0; i < res.component_count; ++i) {
		sbmf_log_info("    [%u]: size: %u", i, res.hamiltonian[i].size);
		sbmf_log_info("    [%u]: bandcount: %u", i, res.hamiltonian[i].bandcount);
	}

	return res;
}

max_stars_repo_name: AntonJohansson/sbmf

max_issues_repo_name: AntonJohansson/sbmf

max_forks_repo_name: AntonJohansson/sbmf


text:
/*
Copyright (C) 2019-2020 JingWeiZhangHuai <jingweizhanghuai@163.com>
Licensed under the Apache License, Version 2.0; 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.
*/
#include "morn_util.h"

#include <cblas.h>
#ifdef MORN_USE_CL
#include <CL/cl.h>
#include <clBLAS.h>
cl_context mDeviceContext(int device);
cl_command_queue mDeviceQueue(int device);
#endif

#define MORN_NO_TRANS 0
#define MORN_TRANS    1

struct HandleCLBlas
{
    int clblas_setup;
};
void endCLBlas(struct HandleCLBlas *handle)
{
    NULL;
    #ifdef MORN_USE_CL
    if(handle->clblas_setup) clblasTeardown();
    #endif
}
#define HASH_CLBlas 0x3b4b7c3a
static struct HandleCLBlas *morn_clblas_handle=NULL;
void SetupCLBlas()
{
    MHandle *hdl=mHandle(mMornObject(NULL,DFLT),CLBlas);
    morn_clblas_handle=(struct HandleCLBlas *)(hdl->handle);
    if(hdl->valid==0)
    {
        #ifdef MORN_USE_CL
        mException(clblasSetup()!=CL_SUCCESS,EXIT,"CLBlas not work");
        morn_clblas_handle->clblas_setup=1;
        #endif
        hdl->valid=1;
    }
}

void mSgemm(int device,int a_trans,int b_trans,int m,int n,int k,float alpha,MMemoryBlock *a,int sa,MMemoryBlock *b,int sb,float beta,MMemoryBlock *c,int sc)
{
    #ifdef MORN_USE_CL
    if(device!=MORN_HOST)
    {
        a_trans=(a_trans==MORN_TRANS)?clblasTrans:clblasNoTrans;
        b_trans=(b_trans==MORN_TRANS)?clblasTrans:clblasNoTrans;

        if(morn_clblas_handle==NULL) SetupCLBlas();

        device = c->device;
        mMemoryBlockCopy(a,device);
        mMemoryBlockCopy(b,device);
        cl_event a_event=(cl_event)(a->cl_evt);
        cl_event b_event=(cl_event)(b->cl_evt);
        cl_event c_event=(cl_event)(c->cl_evt);

        cl_event event_list[2] = {a_event,b_event};
        cl_command_queue queue=mDeviceQueue(device);
        int ret=clblasSgemm(clblasRowMajor,a_trans,b_trans,m,n,k,alpha,a->cl_data,0,sa,b->cl_data,0,sb,beta,c->cl_data,0,sc,1,&queue,2,event_list,&c_event);
        mException(ret!=CL_SUCCESS,EXIT,"clblas error");

        c->flag = MORN_DEVICE;
        return;
    }
    #endif
    a_trans=(a_trans==MORN_TRANS)?CblasTrans:CblasNoTrans;
    b_trans=(b_trans==MORN_TRANS)?CblasTrans:CblasNoTrans;
    cblas_sgemm(CblasRowMajor,a_trans,b_trans,m,n,k,alpha,a->data,sa,b->data,sb,beta,c->data,sc);
}

char *saxpby=mString(
__kernel void saxpby(__global const float* a,__global float* b,const float alpha,const float beta,const int sa,const int sb)
{
    const int idx = get_global_id(0);
    const int ia = idx*sa;
    const int ib = idx*sb;
    b[ib]=a[ia]*alpha+b[ib]*beta;
});

void mSaxpby(int device,int n,float alpha,MMemoryBlock *a,int sa,float beta,MMemoryBlock *b,int sb)
{
    #ifdef MORN_USE_CL
    if(device!=MORN_HOST)
    {
        mCLFunction(saxpby,CLSIZE(n),CLIN(a),CLINOUT(b),CLPARA(&alpha,sizeof(float)),CLPARA(&beta,sizeof(float)),CLPARA(&sa,sizeof(int)),CLPARA(&sb,sizeof(int)));
        return;
    }
    #endif
    cblas_saxpby(n,alpha,a->data,sa,beta,b->data,sb);
}

max_stars_repo_name: ishine/Morn

max_issues_repo_name: ishine/Morn

max_forks_repo_name: ishine/Morn


text:
#pragma once

#include "continuation_scheduler.h"

#include <arcana/threading/cancellation.h>
#include <arcana/threading/task.h>

#include <gsl/gsl>

#include <condition_variable>
#include <mutex>
#include <optional>

namespace Babylon::Graphics
{
    class SafeTimespanGuarantor
    {
    public:
        SafeTimespanGuarantor(std::optional<arcana::cancellation_source>&);

        continuation_scheduler<>& OpenScheduler()
        {
            return m_openDispatcher.scheduler();
        }

        continuation_scheduler<>& CloseScheduler()
        {
            return m_closeDispatcher.scheduler();
        }

        using SafetyGuarantee = gsl::final_action<std::function<void()>>;
        SafetyGuarantee GetSafetyGuarantee();

        void Open();
        void RequestClose();
        void Lock();
        void Unlock();

    private:
        enum class State
        {
            Open,
            Closing,
            Closed,
            Locked
        };

        std::optional<arcana::cancellation_source>& m_cancellation;
        State m_state{State::Locked};
        uint32_t m_count{};
        std::mutex m_mutex{};
        std::condition_variable m_condition_variable{};
        continuation_dispatcher<> m_openDispatcher{};
        continuation_dispatcher<> m_closeDispatcher{};
    };
}

max_stars_repo_name: SergioRZMasson/BabylonNative

max_issues_repo_name: SergioRZMasson/BabylonNative

max_forks_repo_name: SergioRZMasson/BabylonNative


text:
#include <stdio.h>
#include <errno.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_permutation.h>
#include "invwishpdf.h"

int main(int argc, char *argv[]){

	int args;
	int m, n;
	double pdf, dof;

	args = 1;

	m = atoi(argv[args++]);
	n = atoi(argv[args++]);

	// allocate space for matrix X, test_copy, and Scale.
	gsl_matrix *X = gsl_matrix_alloc(m, n);
	gsl_matrix *inv = gsl_matrix_alloc(m,n);
	gsl_matrix *Scale = gsl_matrix_alloc(m, n);

	gsl_matrix_set_identity(Scale);
	fill_matrix(X);

	dof = m + 1;
	pdf = iwishpdf(X, Scale, inv, dof);

	// Print the data
	printf("\n%s", "Inverse Wishart matrix");
	print_matrix(X);
	printf("\n%s", "Scale matrix");
	print_matrix(Scale);
	printf("\n%s", "Wishart matrix");
	print_matrix(inv);
	printf("\n%s", "Density");
	printf("\n%f\n", pdf);

	return(0);
}

gsl_matrix fill_matrix(gsl_matrix *X)
{
	int i,j;
	int step = 0;
	int n = X->size1;
	int m = X->size2;

	// Create a test matrix
	double a_vect[] = {0.17955079, -0.09332222, -0.01273225, -0.09332222,  0.28282278,  0.01724885, -0.01273225,  0.01724885,  0.14424177};
	for (i = 0; i < m; ++i){
		for (j = 0; j < n; ++j){
			gsl_matrix_set(X, i, j, a_vect[step++]);
		}
	}
	return(*X);
}

double iwishpdf(gsl_matrix *X, gsl_matrix *Scale, gsl_matrix *inv, double dof)
{
	double X_det, scale_det, denom, pdf, trace, numerator;
	int m = X->size1;
	int n = X->size2;

	// Allocate matrix for inv(X)
	gsl_matrix *for_mult = gsl_matrix_alloc(m, n);

	// Get determinant of X and Scale matrix
	X_det = matrix_determ(X);
	scale_det = matrix_determ(Scale);

	// Invert X
	inv_matrix(X,inv);

	// Multiple Scale * inv(X)
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Scale, inv, 1.0, for_mult);

	// Get trace of above.
	trace = matrix_trace(for_mult);

	numerator =  pow(scale_det, dof / 2.0) * pow(X_det, (-dof-m-1)/ 2.0) * exp(-0.5 * trace);
	denom = pow(2,dof * m / 2) * mv_gamma(dof/2, m);

	pdf = (numerator/denom);

	return(pdf);
}

gsl_matrix inv_matrix(gsl_matrix *X, gsl_matrix *inv)
{
	int s;
	int n = X->size1;
	int m = X->size2;
	gsl_matrix *a_copy = gsl_matrix_alloc(m, n);
	gsl_matrix_memcpy( a_copy, X );
	gsl_permutation *P = gsl_permutation_alloc(n);
	gsl_linalg_LU_decomp(a_copy, P, &s);
	gsl_linalg_LU_invert(a_copy, P, inv);
	return(*inv);
}

double matrix_determ(gsl_matrix *X)
{
	int s;
	int n = X->size1;
	int m = X->size2;
	gsl_matrix *a_copy = gsl_matrix_alloc(m, n);
	gsl_matrix_memcpy(a_copy, X );
	gsl_permutation *P = gsl_permutation_alloc(n);
	gsl_linalg_LU_decomp(a_copy, P, &s);
	double my_det = gsl_linalg_LU_det (a_copy, s);
	return(my_det);
}

double matrix_trace(gsl_matrix *X){
	int i, m;
	m = X->size1;
	double trace = 0.0;
	for(i=0;i<m;i++){
		trace += gsl_matrix_get(X,i,i);
	}
	return(trace);
}

double mv_gamma(double a, double d){
	double val = 1.0;
	int i;
	for(i = 1; i <= d; i++){
		val *= gsl_sf_gamma(a - (0.5 * (i - 1)));
	}
	val *=  pow(M_PI, (d * (d - 1) / 4.0));
	return(val);
}

void print_matrix(gsl_matrix *X)
{
	int i, j;
	int n = X->size1;
	int m = X->size2;
	printf("\n");
	for(i=0; i<m; i++){
		for(j=0; j<n; j++){
			printf("%05.2f ", gsl_matrix_get(X, i, j));
		}
		printf("\n");
	}
}


text:
#ifndef ARRUS_COMMON_FORMAT_H
#define ARRUS_COMMON_FORMAT_H

// String formatting and parsing utilities.
// Currently wraps fmt library calls.
#include <fmt/format.h>
#include <stdexcept>
#include <cctype>
#include <vector>
#include <set>
#include <unordered_set>
#include <optional>
#include <sstream>
#include <gsl/span>

#include <boost/algorithm/string/join.hpp>

#include "arrus/core/api/common/Tuple.h"
#include "arrus/core/api/common/Interval.h"

namespace arrus {

template<typename... Args>
auto format(Args &&... args) {
    return fmt::format(std::forward<Args>(args)...);
}

/**
 * Returns true if the given string contains numeric characters only.
 *
 * @param num string to verify
 * @return true if the given string contains numeric characters only,
 *         false otherwise.
 */
inline bool isDigitsOnly(const std::string &num) {
    return std::all_of(num.begin(), num.end(), isdigit);
}

/**
 * General purpose 'toString' function basing on ostream << operator.
 */
template<typename T>
std::string toString(const T &t)  {
    std::ostringstream ss;
    ss << t;
    return ss.str();
}

template<typename T>
inline std::string toString(const std::vector<T> &values) {
    std::vector<std::string> vStr(values.size());
    std::transform(std::begin(values), std::end(values), std::begin(vStr),
                   [](auto v) { return std::to_string(v); });
    return boost::algorithm::join(vStr, ", ");
}

template<typename T>
inline std::string toString(
        const gsl::span<T> &values) {
    std::vector<std::string> vStr(values.size());
    std::transform(std::begin(values), std::end(values), std::begin(vStr),
                   [](auto v) { return std::to_string(v); });
    return boost::algorithm::join(vStr, ", ");
}

template<typename T>
inline std::string toStringTransform(
        const std::vector<T> &values,
        const std::function<std::string(const T&)> &func) {
    std::vector<std::string> vStr(values.size());
    std::transform(std::begin(values), std::end(values), std::begin(vStr),
                   [&func](T v) { return func(v); });
    return boost::algorithm::join(vStr, ", ");
}

template<typename T>
inline std::string toString(const ::std::set<T> &values) {
    std::vector<std::string> vStr(values.size());
    std::transform(std::begin(values), std::end(values), std::begin(vStr),
                   [](auto v) { return std::to_string(v); });
    return boost::algorithm::join(vStr, ", ");
}

template<typename T>
inline std::string toString(const ::std::unordered_set<T> &values) {
    std::vector<std::string> vStr(values.size());
    std::transform(std::begin(values), std::end(values), std::begin(vStr),
                   [](auto v) { return std::to_string(v); });
    return boost::algorithm::join(vStr, ", ");
}

template<typename T>
inline std::string toString(const std::optional<T> value) {
    if(value.has_value()) {
        return std::to_string(value.value());
    } else return "(no value)";
}

template<typename T>
inline std::string toString(const Tuple<T> tuple) {
    return ::arrus::format("Tuple({})", toString(tuple.getValues()));
}

template<typename T>
inline std::string toString(const Interval<T> i) {
    return ::arrus::format("Interval: start: {}, right: {}",
                           i.start(), i.end());
}

}

#endif //ARRUS_COMMON_FORMAT_H


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

#include <gsl/gsl_math.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

int
main()
{
  const gsl_interp2d_type *T = gsl_interp2d_bilinear;
  const size_t N = 100;             /* number of points to interpolate */
  const double xa[] = { 0.0, 1.0 }; /* define unit square */
  const double ya[] = { 0.0, 1.0 };
  const size_t nx = sizeof(xa) / sizeof(double); /* x grid points */
  const size_t ny = sizeof(ya) / sizeof(double); /* y grid points */
  double *za = malloc(nx * ny * sizeof(double));
  gsl_spline2d *spline = gsl_spline2d_alloc(T, nx, ny);
  gsl_interp_accel *xacc = gsl_interp_accel_alloc();
  gsl_interp_accel *yacc = gsl_interp_accel_alloc();
  size_t i, j;

  /* set z grid values */
  gsl_spline2d_set(spline, za, 0, 0, 0.0);
  gsl_spline2d_set(spline, za, 0, 1, 1.0);
  gsl_spline2d_set(spline, za, 1, 1, 0.5);
  gsl_spline2d_set(spline, za, 1, 0, 1.0);

  /* initialize interpolation */
  gsl_spline2d_init(spline, xa, ya, za, nx, ny);

  /* interpolate N values in x and y and print out grid for plotting */
  for (i = 0; i < N; ++i)
    {
      double xi = i / (N - 1.0);

      for (j = 0; j < N; ++j)
        {
          double yj = j / (N - 1.0);
          double zij = gsl_spline2d_eval(spline, xi, yj, xacc, yacc);

          printf("%f %f %f\n", xi, yj, zij);
        }
      printf("\n");
    }

  gsl_spline2d_free(spline);
  gsl_interp_accel_free(xacc);
  gsl_interp_accel_free(yacc);
  free(za);

  return 0;
}