text:
#ifndef VECTORBOOST
#define VECTORBOOST

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/qvm/mat_operations.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <lapacke.h>

#include <random>

#include <iostream>

class matrixBoost{
        private:

                boost::numeric::ublas::matrix<double> value;
		boost::numeric::ublas::range r1;
		boost::numeric::ublas::range r2;
		boost::numeric::ublas::matrix<double>valBack;
		bool subMatrixCopy = false;

		//http://www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl?LU_Matrix_Inversion
		bool InvertMatrix(const boost::numeric::ublas::matrix<double>& input, boost::numeric::ublas::matrix<double>& inverse) {
			using namespace boost::numeric::ublas;
			typedef permutation_matrix<std::size_t> pmatrix;
			// create a working copy of the input
			matrix<double> A(input);
			// create a permutation matrix for the LU-factorization
			pmatrix pm(A.size1());
			// perform LU-factorization
			int res = lu_factorize(A,pm);
			if( res != 0 ) return false;
			// create identity matrix of "inverse"
			inverse.assign(boost::numeric::ublas::identity_matrix<double>(A.size1()));
			// backsubstitute to get the inverse
			lu_substitute(A, pm, inverse);
			return true;
		}

		void setSystemValue(boost::numeric::ublas::matrix<double> v){
			value = v;
		}

        public:
                matrixBoost(){}

		matrixBoost(int rowCount, int colCount){
			value = boost::numeric::ublas::zero_matrix<double>(rowCount,colCount);
		}

                matrixBoost(boost::numeric::ublas::matrix<double> value):value(value){}

                static matrixBoost Identity(int dimension,int dimension2){
                        return matrixBoost(boost::numeric::ublas::identity_matrix<double>(dimension,dimension2));
                }

                matrixBoost transpose(){
                        matrixBoost t(boost::numeric::ublas::trans(value));
                        return t;
                }

                matrixBoost inverse(){
			boost::numeric::ublas::matrix<double> inverse(value.size1(),value.size2());
			InvertMatrix(value,inverse);
                        matrixBoost t(inverse);
                        return t;
                }

                boost::numeric::ublas::matrix<double> getSystemValue() const{
                        return value;
                }

		matrixBoost & block(int i1, int i2, int l1, int l2){
			if(subMatrixCopy){
				std::swap(valBack,value);
			}
			subMatrixCopy = true;
			r1 = boost::numeric::ublas::range(i1,i1+l1);
			r2 = boost::numeric::ublas::range(i2,i2+l2);
			valBack = boost::numeric::ublas::project(value,r1,r2);
			std::swap(valBack,value);
			return *this;
		}

		matrixBoost &operator=( const matrixBoost other){
			if(subMatrixCopy){
				subMatrixCopy = false;
				boost::numeric::ublas::project(valBack,r1,r2) = other.getSystemValue();
				std::swap(valBack,value);
			}else{
				value = other.getSystemValue();
			}
			return *this;//better be no chaining of equal signs
		}

		static void lowerTriangulate(matrixBoost & val){
			boost::numeric::ublas::matrix<double> in= val.getSystemValue();
			int rows = in.size1();
			int cols = in.size2();
			for(int i = 0; i < rows; i++){
				int length = cols - i;
				boost::numeric::ublas::vector<double> v = (row(in,i));
				v = boost::numeric::ublas::subrange(v,i,cols);
				double mew = boost::numeric::ublas::norm_2(v);
				if(mew!=0){
					double beta  = v(0)  + (-2*(v(0) < 0)+1)*mew;
					v = v/beta;
				}
				v(0) = 1;
				boost::numeric::ublas::range t1(i,rows);
				boost::numeric::ublas::range t2(i,cols);
				boost::numeric::ublas::vector<double> w = -2.0/
					(boost::numeric::ublas::norm_2(v)*boost::numeric::ublas::norm_2(v))*
					(prod(project(in,t1,t2),v));

				boost::numeric::ublas::matrix<double> vM(v.size(),1);
				for(int i = 0; i < v.size(); i++){
					vM(i,0) = v(i);
				}
				boost::numeric::ublas::matrix<double> wM(w.size(),1);
				for(int i = 0; i < w.size(); i++){
					wM(i,0) = w(i);
				}
				project(in,t1,t2) = project(in,t1,t2) +  prod(wM,boost::numeric::ublas::trans(vM));
			}
			val.setSystemValue(in);
		}
};

class vectorBoost{
        private:
		boost::numeric::ublas::vector<double> value;
        public:
                vectorBoost(){
                }

                vectorBoost(boost::numeric::ublas::vector<double> value):value(value){}

                boost::numeric::ublas::vector<double> getSystemValue() const{
                        return value;
                }

		static vectorBoost randomVector(int dimension,uint32_t seed){
			std::mt19937 generator(seed);
			std::normal_distribution<double> distribution(0.0,1.0);
			boost::numeric::ublas::vector<double> tmp(dimension);
			for(int i = 0; i < dimension; i++){
				tmp(i) = distribution(generator);
			}
			return vectorBoost(tmp);
		}

		static vectorBoost solve(matrixBoost A, vectorBoost b){
			boost::numeric::ublas::matrix<double> Alocal = A.getSystemValue();
			boost::numeric::ublas::vector<double> blocal = b.getSystemValue();
			boost::numeric::ublas::permutation_matrix<size_t> pm(Alocal.size1());

			boost::numeric::ublas::lu_factorize(Alocal,pm);
			boost::numeric::ublas::lu_substitute(Alocal, pm, blocal);
			return vectorBoost(blocal);
		}
};

vectorBoost operator-(const vectorBoost &lhs,const vectorBoost &rhs){
        vectorBoost res(lhs.getSystemValue() - rhs.getSystemValue());
        return res;
};

vectorBoost operator+(const vectorBoost &lhs,const vectorBoost &rhs){
        vectorBoost res(lhs.getSystemValue() + rhs.getSystemValue());
        return res;
};

std::ostream& operator<<(std::ostream& stream, const vectorBoost & vec) {
	boost::numeric::ublas::vector<double> v = vec.getSystemValue();
	for(int i = 0; i < v.size()-1; i++){
		stream<<double(v(i))<<'\n';
	}
	stream<<v(v.size()-1);
	return stream;
}

matrixBoost operator*(const matrixBoost &lhs, const matrixBoost &rhs){
        matrixBoost res(boost::numeric::ublas::prod(lhs.getSystemValue(),rhs.getSystemValue()));
        return res;
}

matrixBoost operator+(const matrixBoost &lhs,const matrixBoost &rhs){
        matrixBoost res(lhs.getSystemValue() + rhs.getSystemValue());
        return res;
};

matrixBoost operator-(const matrixBoost &lhs,const matrixBoost &rhs){
        matrixBoost res(lhs.getSystemValue() - rhs.getSystemValue());
        return res;
};

std::ostream& operator<<(std::ostream& stream, const matrixBoost & mat) {
	boost::numeric::ublas::matrix<double> m = mat.getSystemValue();
	for(int i = 0; i < m.size1(); i++){
		for(int j = 0; j < m.size2(); j++){
			stream<<m(i,j)<<" ";
		}stream<<'\n';
	}
	return stream;
}

vectorBoost operator*(const matrixBoost &lhs, const vectorBoost &rhs){
        vectorBoost res(boost::numeric::ublas::prod(lhs.getSystemValue(),rhs.getSystemValue()));
        return res;
}

vectorBoost operator*(const vectorBoost &lhs, const matrixBoost &rhs){
        vectorBoost res(boost::numeric::ublas::prod(lhs.getSystemValue(),rhs.getSystemValue()));
        return res;
}
#endif

max_stars_repo_name: mannyray/KalmanFilter

max_issues_repo_name: mannyray/KalmanFilter

max_forks_repo_name: mannyray/KalmanFilter


text:
/* matrix/gsl_matrix_ulong.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_ULONG_H__
#define __GSL_MATRIX_ULONG_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_ulong.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;
  unsigned long * data;
  gsl_block_ulong * block;
  int owner;
} gsl_matrix_ulong;

typedef struct
{
  gsl_matrix_ulong matrix;
} _gsl_matrix_ulong_view;

typedef _gsl_matrix_ulong_view gsl_matrix_ulong_view;

typedef struct
{
  gsl_matrix_ulong matrix;
} _gsl_matrix_ulong_const_view;

typedef const _gsl_matrix_ulong_const_view gsl_matrix_ulong_const_view;

/* Allocation */

GSL_EXPORT
gsl_matrix_ulong *
gsl_matrix_ulong_alloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_ulong *
gsl_matrix_ulong_calloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_ulong *
gsl_matrix_ulong_alloc_from_block (gsl_block_ulong * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_EXPORT
gsl_matrix_ulong *
gsl_matrix_ulong_alloc_from_matrix (gsl_matrix_ulong * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
gsl_vector_ulong *
gsl_vector_ulong_alloc_row_from_matrix (gsl_matrix_ulong * m,
                                        const size_t i);

GSL_EXPORT
gsl_vector_ulong *
gsl_vector_ulong_alloc_col_from_matrix (gsl_matrix_ulong * m,
                                        const size_t j);

GSL_EXPORT void gsl_matrix_ulong_free (gsl_matrix_ulong * m);

/* Views */

GSL_EXPORT
_gsl_matrix_ulong_view
gsl_matrix_ulong_submatrix (gsl_matrix_ulong * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_ulong_view
gsl_matrix_ulong_row (gsl_matrix_ulong * m, const size_t i);

GSL_EXPORT
_gsl_vector_ulong_view
gsl_matrix_ulong_column (gsl_matrix_ulong * m, const size_t j);

GSL_EXPORT
_gsl_vector_ulong_view
gsl_matrix_ulong_diagonal (gsl_matrix_ulong * m);

GSL_EXPORT
_gsl_vector_ulong_view
gsl_matrix_ulong_subdiagonal (gsl_matrix_ulong * m, const size_t k);

GSL_EXPORT
_gsl_vector_ulong_view
gsl_matrix_ulong_superdiagonal (gsl_matrix_ulong * m, const size_t k);

GSL_EXPORT
_gsl_matrix_ulong_view
gsl_matrix_ulong_view_array (unsigned long * base,
                             const size_t n1,
                             const size_t n2);

GSL_EXPORT
_gsl_matrix_ulong_view
gsl_matrix_ulong_view_array_with_tda (unsigned long * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_EXPORT
_gsl_matrix_ulong_view
gsl_matrix_ulong_view_vector (gsl_vector_ulong * v,
                              const size_t n1,
                              const size_t n2);

GSL_EXPORT
_gsl_matrix_ulong_view
gsl_matrix_ulong_view_vector_with_tda (gsl_vector_ulong * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_EXPORT
_gsl_matrix_ulong_const_view
gsl_matrix_ulong_const_submatrix (const gsl_matrix_ulong * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_ulong_const_view
gsl_matrix_ulong_const_row (const gsl_matrix_ulong * m,
                            const size_t i);

GSL_EXPORT
_gsl_vector_ulong_const_view
gsl_matrix_ulong_const_column (const gsl_matrix_ulong * m,
                               const size_t j);

GSL_EXPORT
_gsl_vector_ulong_const_view
gsl_matrix_ulong_const_diagonal (const gsl_matrix_ulong * m);

GSL_EXPORT
_gsl_vector_ulong_const_view
gsl_matrix_ulong_const_subdiagonal (const gsl_matrix_ulong * m,
                                    const size_t k);

GSL_EXPORT
_gsl_vector_ulong_const_view
gsl_matrix_ulong_const_superdiagonal (const gsl_matrix_ulong * m,
                                      const size_t k);

GSL_EXPORT
_gsl_matrix_ulong_const_view
gsl_matrix_ulong_const_view_array (const unsigned long * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_EXPORT
_gsl_matrix_ulong_const_view
gsl_matrix_ulong_const_view_array_with_tda (const unsigned long * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_EXPORT
_gsl_matrix_ulong_const_view
gsl_matrix_ulong_const_view_vector (const gsl_vector_ulong * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
_gsl_matrix_ulong_const_view
gsl_matrix_ulong_const_view_vector_with_tda (const gsl_vector_ulong * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_EXPORT unsigned long   gsl_matrix_ulong_get(const gsl_matrix_ulong * m, const size_t i, const size_t j);
GSL_EXPORT void    gsl_matrix_ulong_set(gsl_matrix_ulong * m, const size_t i, const size_t j, const unsigned long x);

GSL_EXPORT unsigned long * gsl_matrix_ulong_ptr(gsl_matrix_ulong * m, const size_t i, const size_t j);
GSL_EXPORT const unsigned long * gsl_matrix_ulong_const_ptr(const gsl_matrix_ulong * m, const size_t i, const size_t j);

GSL_EXPORT void gsl_matrix_ulong_set_zero (gsl_matrix_ulong * m);
GSL_EXPORT void gsl_matrix_ulong_set_identity (gsl_matrix_ulong * m);
GSL_EXPORT void gsl_matrix_ulong_set_all (gsl_matrix_ulong * m, unsigned long x);

GSL_EXPORT int gsl_matrix_ulong_fread (FILE * stream, gsl_matrix_ulong * m) ;
GSL_EXPORT int gsl_matrix_ulong_fwrite (FILE * stream, const gsl_matrix_ulong * m) ;
GSL_EXPORT int gsl_matrix_ulong_fscanf (FILE * stream, gsl_matrix_ulong * m);
GSL_EXPORT int gsl_matrix_ulong_fprintf (FILE * stream, const gsl_matrix_ulong * m, const char * format);

GSL_EXPORT int gsl_matrix_ulong_memcpy(gsl_matrix_ulong * dest, const gsl_matrix_ulong * src);
GSL_EXPORT int gsl_matrix_ulong_swap(gsl_matrix_ulong * m1, gsl_matrix_ulong * m2);

GSL_EXPORT int gsl_matrix_ulong_swap_rows(gsl_matrix_ulong * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_ulong_swap_columns(gsl_matrix_ulong * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_ulong_swap_rowcol(gsl_matrix_ulong * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_ulong_transpose (gsl_matrix_ulong * m);
GSL_EXPORT int gsl_matrix_ulong_transpose_memcpy (gsl_matrix_ulong * dest, const gsl_matrix_ulong * src);

GSL_EXPORT unsigned long gsl_matrix_ulong_max (const gsl_matrix_ulong * m);
GSL_EXPORT unsigned long gsl_matrix_ulong_min (const gsl_matrix_ulong * m);
GSL_EXPORT void gsl_matrix_ulong_minmax (const gsl_matrix_ulong * m, unsigned long * min_out, unsigned long * max_out);

GSL_EXPORT void gsl_matrix_ulong_max_index (const gsl_matrix_ulong * m, size_t * imax, size_t *jmax);
GSL_EXPORT void gsl_matrix_ulong_min_index (const gsl_matrix_ulong * m, size_t * imin, size_t *jmin);
GSL_EXPORT void gsl_matrix_ulong_minmax_index (const gsl_matrix_ulong * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_EXPORT int gsl_matrix_ulong_isnull (const gsl_matrix_ulong * m);

GSL_EXPORT int gsl_matrix_ulong_add (gsl_matrix_ulong * a, const gsl_matrix_ulong * b);
GSL_EXPORT int gsl_matrix_ulong_sub (gsl_matrix_ulong * a, const gsl_matrix_ulong * b);
GSL_EXPORT int gsl_matrix_ulong_mul_elements (gsl_matrix_ulong * a, const gsl_matrix_ulong * b);
GSL_EXPORT int gsl_matrix_ulong_div_elements (gsl_matrix_ulong * a, const gsl_matrix_ulong * b);
GSL_EXPORT int gsl_matrix_ulong_scale (gsl_matrix_ulong * a, const double x);
GSL_EXPORT int gsl_matrix_ulong_add_constant (gsl_matrix_ulong * a, const double x);
GSL_EXPORT int gsl_matrix_ulong_add_diagonal (gsl_matrix_ulong * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_EXPORT int gsl_matrix_ulong_get_row(gsl_vector_ulong * v, const gsl_matrix_ulong * m, const size_t i);
GSL_EXPORT int gsl_matrix_ulong_get_col(gsl_vector_ulong * v, const gsl_matrix_ulong * m, const size_t j);
GSL_EXPORT int gsl_matrix_ulong_set_row(gsl_matrix_ulong * m, const size_t i, const gsl_vector_ulong * v);
GSL_EXPORT int gsl_matrix_ulong_set_col(gsl_matrix_ulong * m, const size_t j, const gsl_vector_ulong * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
unsigned long
gsl_matrix_ulong_get(const gsl_matrix_ulong * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_ulong_set(gsl_matrix_ulong * m, const size_t i, const size_t j, const unsigned long 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
  m->data[i * m->tda + j] = x ;
}

extern inline
unsigned long *
gsl_matrix_ulong_ptr(gsl_matrix_ulong * 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 (unsigned long *) (m->data + (i * m->tda + j)) ;
}

extern inline
const unsigned long *
gsl_matrix_ulong_const_ptr(const gsl_matrix_ulong * 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 unsigned long *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_ULONG_H__ */

ext: h

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/** @file */
#ifndef __CCL_CORE_H_INCLUDED__
#define __CCL_CORE_H_INCLUDED__

#include <stdbool.h>
#include <stdio.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_const_mksa.h>

#include "ccl_utils.h"
#include "ccl_f1d.h"
#include "ccl_f2d.h"

CCL_BEGIN_DECLS

/**
 * Struct to hold physical constants.
 */
typedef struct ccl_physical_constants {
  /**
   * Lightspeed / H0 in units of Mpc/h (from CODATA 2014)
   */
  double CLIGHT_HMPC;

  /**
   * Newton's gravitational constant in units of m^3/Kg/s^2
   */
  double GNEWT;

  /**
   * Solar mass in units of kg (from GSL)
   */
  double SOLAR_MASS;

  /**
   * Mpc to meters (from PDG 2013)
   */
  double MPC_TO_METER;

  /**
   * pc to meters (from PDG 2013)
   */
  double PC_TO_METER;

  /**
   * Rho critical in units of M_sun/h / (Mpc/h)^3
   */
  double RHO_CRITICAL;

  /**
   * Boltzmann constant in units of J/K
  */
  double KBOLTZ;

  /**
   * Stefan-Boltzmann constant in units of kg/s^3 / K^4
   */
  double STBOLTZ;

  /**
   * Planck's constant in units kg m^2 / s
   */
  double HPLANCK;

  /**
   * The speed of light in m/s
   */
  double CLIGHT;

  /**
   * Electron volt to Joules convestion
   */
  double EV_IN_J;

  /**
   * Temperature of the CMB in K
   */
  double T_CMB;

  /**
   * T_ncdm, as taken from CLASS, explanatory.ini
   */
  double TNCDM;

  /**
   * neutrino mass splitting differences
   * See Lesgourgues and Pastor, 2012 for these values.
   * Adv. High Energy Phys. 2012 (2012) 608515,
   * arXiv:1212.6154, page 13
  */
  double DELTAM12_sq;
  double DELTAM13_sq_pos;
  double DELTAM13_sq_neg;
} ccl_physical_constants;

extern ccl_physical_constants ccl_constants;

/**
 * Struct that contains all the parameters needed to create certain splines.
 * This includes splines for the scale factor, masses, and power spectra.
 */
typedef struct ccl_spline_params {
  // scale factor splines
  int A_SPLINE_NA;
  double A_SPLINE_MIN;
  double A_SPLINE_MINLOG_PK;
  double A_SPLINE_MIN_PK;
  double A_SPLINE_MINLOG_SM;
  double A_SPLINE_MIN_SM;
  double A_SPLINE_MAX;
  double A_SPLINE_MINLOG;
  int A_SPLINE_NLOG;

  //Mass splines
  double LOGM_SPLINE_DELTA;
  int LOGM_SPLINE_NM;
  double LOGM_SPLINE_MIN;
  double LOGM_SPLINE_MAX;

  //PS a and k spline
  int A_SPLINE_NA_SM;
  int A_SPLINE_NLOG_SM;
  int A_SPLINE_NA_PK;
  int A_SPLINE_NLOG_PK;

  //k-splines and integrals
  double K_MAX_SPLINE;
  double K_MAX;
  double K_MIN;
  double DLOGK_INTEGRATION;
  int N_K;
  int N_K_3DCOR;

  //Correlation function parameters
  double ELL_MIN_CORR;
  double ELL_MAX_CORR;
  int N_ELL_CORR;

  // interpolation types
  const gsl_interp_type* A_SPLINE_TYPE;
  const gsl_interp_type* K_SPLINE_TYPE;
  const gsl_interp_type* M_SPLINE_TYPE;
  const gsl_interp_type* D_SPLINE_TYPE;
  const gsl_interp2d_type* PNL_SPLINE_TYPE;
  const gsl_interp2d_type* PLIN_SPLINE_TYPE;
  const gsl_interp_type* CORR_SPLINE_TYPE;
} ccl_spline_params;

extern const ccl_spline_params default_spline_params;

/**
 * Struct that contains parameters that control the accuracy of various GSL
 * routines.
 */
typedef struct ccl_gsl_params {
  // General parameters
  size_t N_ITERATION;

  // Integration
  int INTEGRATION_GAUSS_KRONROD_POINTS;
  double INTEGRATION_EPSREL;
  // Limber integration
  int INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS;
  double INTEGRATION_LIMBER_EPSREL;
  // Distance integrals
  double INTEGRATION_DISTANCE_EPSREL;
  // sigma_R integral
  double INTEGRATION_SIGMAR_EPSREL;
  // k_NL integral
  double INTEGRATION_KNL_EPSREL;

  // Root finding
  double ROOT_EPSREL;
  int ROOT_N_ITERATION;

  // ODE
  double ODE_GROWTH_EPSREL;

  // growth
  double EPS_SCALEFAC_GROWTH;

  // halo model
  double HM_MMIN;
  double HM_MMAX;
  double HM_EPSABS;
  double HM_EPSREL;
  size_t HM_LIMIT;
  int HM_INT_METHOD;

} ccl_gsl_params;

extern const ccl_gsl_params default_gsl_params;

/**
 * Struct containing the parameters defining a cosmology
 */
typedef struct ccl_parameters {

  // Densities: CDM, baryons, total matter, neutrinos, curvature
  double Omega_c; /**< Density of CDM relative to the critical density*/
  double Omega_b; /**< Density of baryons relative to the critical density*/
  double Omega_m; /**< Density of all matter relative to the critical density*/
  double Omega_k; /**< Density of curvature relative to the critical density*/
  double sqrtk; /**< Square root of the magnitude of curvature, k */ //TODO check
  int k_sign; /**<Sign of the curvature k */

  // Dark Energy
  double w0;
  double wa;

  // Hubble parameters
  double H0;
  double h;

  // Neutrino properties

  double Neff; // Effective number of relativistic neutrino species in the early universe.
  int N_nu_mass; // Number of species of neutrinos which are nonrelativistic today
  double N_nu_rel;  // Number of species of neutrinos which are relativistic  today
  double *m_nu;  // total mass of massive neutrinos (This is a pointer so that it can hold multiple masses.)
  double sum_nu_masses; // sum of the neutrino masses.
  double Omega_nu_mass; // Omega_nu for MASSIVE neutrinos
  double Omega_nu_rel; // Omega_nu for MASSLESS neutrinos

  // Primordial power spectra
  double A_s;
  double n_s;

  // Radiation parameters
  double Omega_g;
  double T_CMB;

  // BCM baryonic model parameters
  double bcm_log10Mc;
  double bcm_etab;
  double bcm_ks;

  // mu / Sigma quasistatica parameterisation of modified gravity params
  double mu_0;
  double sigma_0;

  // Derived parameters
  double sigma8;
  double Omega_l;
  double z_star;

  //Modified growth rate
  bool has_mgrowth;
  int nz_mgrowth;
  double *z_mgrowth;
  double *df_mgrowth;
} ccl_parameters;

/**
 * Struct containing references to gsl splines for distance and acceleration calculations
 */
typedef struct ccl_data {
  // These are all functions of the scale factor a.

  // Distances are defined in Mpc
  double growth0;
  gsl_spline * chi;
  gsl_spline * growth;
  gsl_spline * fgrowth;
  gsl_spline * E;
  gsl_spline * achi;

  // Function of Halo mass M
  gsl_spline2d * logsigma;

  // real-space splines for RSD
  ccl_f1d_t* rsd_splines[3];
  double rsd_splines_scalefactor;
} ccl_data;

/**
 * Sturct containing references to instances of the above structs, and boolean flags of precomputed values.
 */
typedef struct ccl_cosmology {
  ccl_parameters    params;
  ccl_configuration config;
  ccl_data          data;
  ccl_spline_params spline_params;
  ccl_gsl_params    gsl_params;

  bool computed_distances;
  bool computed_growth;
  bool computed_sigma;

  int status;
  //this is optional - less tedious than tracking all numerical values for status in error handler function
  char status_message[500];

  // other flags?
} ccl_cosmology;

// Initialization and life cycle of objects
ccl_cosmology * ccl_cosmology_create(ccl_parameters params, ccl_configuration config);

/* Internal function to set the status message safely. */
void ccl_cosmology_set_status_message(ccl_cosmology * cosmo, const char * status_message, ...);

// User-facing creation routines
/**
 * Create a cosmology
 * @param Omega_c Omega_c
 * @param Omega_b Omega_b
 * @param Omega_k Omega_k
 * @param Neff Number of relativistic neutrino species in the early universe
 * @param mnu neutrino mass, either sum or list of length 3
 * @param mnu_type determines neutrino mass convention (ccl_mnu_list, ccl_mnu_sum, ccl_mnu_sum_inverted, ccl_mnu_sum_equal)
 * @param w0 Dark energy EoS parameter
 * @param wa Dark energy EoS parameter
 * @param h Hubble constant in units of 100 km/s/Mpc
 * @param norm_pk the normalization of the power spectrum, either A_s or sigma8
 * @param n_s the power-law index of the power spectrum
 * @param bcm_log10Mc log10 cluster mass, one of the parameters of the BCM model
 * @param bcm_etab ejection radius parameter, one of the parameters of the BCM model
 * @param bcm_ks wavenumber for the stellar profile, one of the parameters of the BCM model
 * @param nz_mgrowth the number of redshifts where the modified growth is provided
 * @param zarr_mgrowth the array of redshifts where the modified growth is provided
 * @param dfarr_mgrowth the modified growth function vector provided
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 * For specific cases see documentation for ccl_error.c
 * @return void
 */
ccl_parameters ccl_parameters_create(double Omega_c, double Omega_b, double Omega_k,
				     double Neff, double* mnu, int n_mnu,
				     double w0, double wa, double h, double norm_pk,
				     double n_s, double bcm_log10Mc, double bcm_etab, double bcm_ks,
				     double mu_0, double sigma_0, int nz_mgrowth, double *zarr_mgrowth,
				     double *dfarr_mgrowth, int *status);

/**
 * Free a parameters struct
 * @param params ccl_parameters struct
 * @return void
 */
void ccl_parameters_free(ccl_parameters * params);

/**
 * Write a cosmology parameters object to a file in yaml format, .
 * @param params Cosmological parameters
 * @param filename Name of file to create and write
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 * @return void
 */
void ccl_parameters_write_yaml(ccl_parameters * params, const char * filename, int * status);

/**
 * Read a cosmology parameters object from a file in yaml format, .
 * @param filename Name of existing file to read from
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 * @return cosmo Cosmological parameters
 */
ccl_parameters ccl_parameters_read_yaml(const char * filename, int *status);

/**
 * Free a cosmology struct
 * @param cosmo Cosmological parameters
 * @return void
 */
void ccl_cosmology_free(ccl_cosmology * cosmo);

int ccl_get_pk_spline_na(ccl_cosmology *cosmo);
int ccl_get_pk_spline_nk(ccl_cosmology *cosmo);
void ccl_get_pk_spline_a_array(ccl_cosmology *cosmo,int ndout,double* doutput,int *status);
void ccl_get_pk_spline_lk_array(ccl_cosmology *cosmo,int ndout,double* doutput,int *status);

CCL_END_DECLS

#endif

max_stars_repo_name: borisbolliet/CCL

max_issues_repo_name: borisbolliet/CCL

max_forks_repo_name: borisbolliet/CCL


text:
#pragma once

#include <gsl/gsl>
#include <Metal/Metal.h>

namespace Halley {
	class MetalVideo;

	class MetalBuffer {
	public:
		enum class Type
		{
			Vertex,
			Index,
			Constant
		};

		MetalBuffer(MetalVideo& video, Type type, size_t initialSize = 0);
		MetalBuffer(MetalBuffer&& other) noexcept;
		~MetalBuffer();

		MetalBuffer(const MetalBuffer& other) = delete;

		MetalBuffer& operator=(const MetalBuffer& other) = delete;
		MetalBuffer& operator=(MetalBuffer&& other) = delete;

		void setData(gsl::span<const gsl::byte> data);
		void bindVertex(id<MTLRenderCommandEncoder> encoder, int bindPoint);
		void bindFragment(id<MTLRenderCommandEncoder> encoder, int bindPoint);
		id<MTLBuffer> getBuffer();

	private:
		MetalVideo& video;
		Type type;
		id<MTLBuffer> buffer = nil;
	};
}

max_stars_repo_name: code-disaster/halley

max_issues_repo_name: code-disaster/halley

max_forks_repo_name: code-disaster/halley


text:
/* randist/exppow.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006 James Theiler, Brian Gough
 * Copyright (C) 2006 Giulio Bottazzi
 *
 * 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.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* The exponential power probability distribution is

   p(x) dx = (1/(2 a Gamma(1+1/b))) * exp(-|x/a|^b) dx

   for -infty < x < infty. For b = 1 it reduces to the Laplace
   distribution.

   The exponential power distribution is related to the gamma
   distribution by E = a * pow(G(1/b),1/b), where E is an exponential
   power variate and G is a gamma variate.

   We use this relation for b < 1. For b >=1 we use rejection methods
   based on the laplace and gaussian distributions which should be
   faster.  For b>4 we revert to the gamma method.

   See P. R. Tadikamalla, "Random Sampling from the Exponential Power
   Distribution", Journal of the American Statistical Association,
   September 1980, Volume 75, Number 371, pages 683-686.

*/

double
gsl_ran_exppow (const gsl_rng * r, const double a, const double b)
{
  if (b < 1 || b > 4)
    {
      double u = gsl_rng_uniform (r);
      double v = gsl_ran_gamma (r, 1 / b, 1.0);
      double z = a * pow (v, 1 / b);

      if (u > 0.5)
        {
          return z;
        }
      else
        {
          return -z;
        }
    }
  else if (b == 1)
    {
      /* Laplace distribution */
      return gsl_ran_laplace (r, a);
    }
  else if (b < 2)
    {
      /* Use laplace distribution for rejection method, from Tadikamalla */

      double x, h, u;

      double B = pow (1 / b, 1 / b);

      do
        {
          x = gsl_ran_laplace (r, B);
          u = gsl_rng_uniform_pos (r);
          h = -pow (fabs (x), b) + fabs (x) / B - 1 + (1 / b);
        }
      while (log (u) > h);

      return a * x;
    }
  else if (b == 2)
    {
      /* Gaussian distribution */
      return gsl_ran_gaussian (r, a / sqrt (2.0));
    }
  else
    {
      /* Use gaussian for rejection method, from Tadikamalla */

      double x, h, u;

      double B = pow (1 / b, 1 / b);

      do
        {
          x = gsl_ran_gaussian (r, B);
          u = gsl_rng_uniform_pos (r);
          h = -pow (fabs (x), b) + (x * x) / (2 * B * B) + (1 / b) - 0.5;
        }
      while (log (u) > h);

      return a * x;
    }
}

double
gsl_ran_exppow_pdf (const double x, const double a, const double b)
{
  double p;
  double lngamma = gsl_sf_lngamma (1 + 1 / b);
  p = (1 / (2 * a)) * exp (-pow (fabs (x / a), b) - lngamma);
  return p;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include "Data/Data.h"
#include <Eigen/Dense>
#include <opencv2/core/eigen.hpp>
#include <opencv2/calib3d.hpp>
#include <opencv2/imgproc.hpp>

#include <assert.h>
#include <array>
#include <gsl/gsl>

#include "eigen.h"

namespace mage
{
    using Quaternion = Eigen::Quaternionf;

    inline Quaternion ToQuat(const cv::Matx33f& rotation)
    {
        Eigen::Matrix3f eigenMat;
        cv::cv2eigen(rotation, eigenMat);

        Quaternion eigenQuat(eigenMat);
        // TODO validate if this normalize is even necessary, eigen might already be doing this
        eigenQuat.normalize();

        return eigenQuat;
    }

    template<int N>
    inline Eigen::Map<Eigen::Matrix<float, N, 1>> ToMap(cv::Vec<float, N>& vec)
    {
        return Eigen::Map<Eigen::Vector3f>{ vec.val };
    }

    template<int N>
    inline Eigen::Map<const Eigen::Matrix<float, N, 1>> ToCMap(const cv::Vec<float, N>& vec)
    {
        return Eigen::Map<const Eigen::Matrix<float, N, 1>>{ vec.val };
    }

    inline Eigen::Map<Eigen::Vector3f> ToMap(cv::Matx31f& vec)
    {
        return Eigen::Map<Eigen::Vector3f>{ vec.val };
    }

    inline Eigen::Map<const Eigen::Vector3f> ToCMap(const cv::Matx31f& vec)
    {
        return Eigen::Map<const Eigen::Vector3f>{ vec.val };
    }

    inline Eigen::Map<Eigen::Vector3f> ToMap(cv::Point3f& vec)
    {
        return Eigen::Map<Eigen::Vector3f>{ &vec.x };
    }

    inline Eigen::Map<const Eigen::Vector3f> ToCMap(const cv::Point3f& vec)
    {
        return Eigen::Map<const Eigen::Vector3f>{ &vec.x };
    }

    inline Eigen::Map<Eigen::Vector2f> ToMap(cv::Point2f& vec)
    {
        return Eigen::Map<Eigen::Vector2f>{ &vec.x };
    }

    inline Eigen::Map<const Eigen::Vector2f> ToCMap(const cv::Point2f& vec)
    {
        return Eigen::Map<const Eigen::Vector2f>{ &vec.x };
    }

    template<size_t N>
    inline cv::Vec<float, N> ToVec(gsl::span<const float, N> values)
    {
        return cv::Vec<float, N>(values.data());
    }

    inline Quaternion QuatFromTwoVectors(const cv::Vec3f& from, const cv::Vec3f& to)
    {
        Quaternion quat;
        quat.setFromTwoVectors(Eigen::Vector3f{ from(0), from(1), from(2) }, Eigen::Vector3f{ to(0), to(1), to(2) });
        return quat;
    }

    inline cv::Matx33f ToMat(const Quaternion& quat)
    {
        cv::Matx33f mat;
        cv::eigen2cv(quat.toRotationMatrix(), mat);
        return mat;
    }

    inline cv::Matx33f ToMat(float pitchRads, float rollRads, float yawRads)
    {
        // Formula for (Pitch * Roll * Yaw) from http://www.songho.ca/opengl/gl_anglestoaxes.html
        float a = pitchRads;
        float b = yawRads;
        float c = rollRads;

        float sa = sinf(a);
        float sb = sinf(b);
        float sc = sinf(c);

        float ca = cosf(a);
        float cb = cosf(b);
        float cc = cosf(c);

        return cv::Matx33f{
            cc * cb,                  -sc,         cc * sb,
            ca * sc * cb + sa * sb,   ca * cc,     ca * sc * sb - sa * cb,
            sa * sc * cb - ca * sb,   sa * cc,     sa * sc * sb + ca * cb
        };
    }

    inline cv::Matx44f To4x4(const cv::Matx34f& mat)
    {
        return{
            mat(0, 0), mat(0, 1), mat(0, 2), mat(0, 3),
            mat(1, 0), mat(1, 1), mat(1, 2), mat(1, 3),
            mat(2, 0), mat(2, 1), mat(2, 2), mat(2, 3),
            0,         0,         0,         1
        };
    }

    inline cv::Matx34f To3x4(const cv::Matx44f& mat)
    {
        return{
            mat(0, 0), mat(0, 1), mat(0, 2), mat(0, 3),
            mat(1, 0), mat(1, 1), mat(1, 2), mat(1, 3),
            mat(2, 0), mat(2, 1), mat(2, 2), mat(2, 3)
        };
    }

    inline cv::Matx44f To4x4(const cv::Matx33f& rot, const cv::Vec3f& trans)
    {
        return{
            rot(0, 0), rot(0, 1), rot(0, 2), trans(0),
            rot(1, 0), rot(1, 1), rot(1, 2), trans(1),
            rot(2, 0), rot(2, 1), rot(2, 2), trans(2),
            0,         0,         0,         1
        };
    }

    inline cv::Matx44f FromQuatAndTrans(const Quaternion& quat, const cv::Vec3f& trans)
    {
        return To4x4(ToMat(quat), trans);
    }

    inline cv::Matx33f Rotation(const cv::Matx44f& mat)
    {
        return{
            mat(0, 0), mat(0, 1), mat(0, 2),
            mat(1, 0), mat(1, 1), mat(1, 2),
            mat(2, 0), mat(2, 1), mat(2, 2)
        };
    }

    inline cv::Matx33f Rotation(const cv::Matx34f& mat)
    {
        return{
            mat(0, 0), mat(0, 1), mat(0, 2),
            mat(1, 0), mat(1, 1), mat(1, 2),
            mat(2, 0), mat(2, 1), mat(2, 2)
        };
    }

    inline cv::Vec3f Translation(const cv::Matx44f& mat)
    {
        return{
            mat(0, 3),
            mat(1, 3),
            mat(2, 3)
        };
    }

    inline cv::Vec3f Translation(const cv::Matx34f& mat)
    {
        return{
            mat(0, 3),
            mat(1, 3),
            mat(2, 3)
        };
    }
	//Swap the returned parameters to suppress x64 compile warning c4324
	//c:\program files(x86)\microsoft visual studio\2017\enterprise\vc\tools\msvc\14.11.25503\include\utility(271) : warning C4324 : 'std::pair<mage::Quaternion,cv::Vec3f>' : structure was padded due to alignment specifier
    inline std::pair<cv::Vec3f, Quaternion> Decompose(const cv::Matx44f& matrix)
    {
        return std::make_pair(Translation(matrix), ToQuat(matrix.get_minor<3, 3>(0, 0)));
    }

    template<int ROWS, int COLS>
    inline bool MatxEqual(const cv::Matx<float, ROWS, COLS>& mat0, const cv::Matx<float, ROWS, COLS>& mat1)
    {
        for (int i = 0; i < ROWS; i++)
        {
            for (int j = 0; j < COLS; j++)
            {
                if (mat0(i,j) != mat1(i,j))
                    return false;
            }
        }

        return true;
    }

    inline bool MatEqual(const cv::Mat& mat0, const cv::Mat& mat1)
    {
        if (mat0.size != mat1.size)
            return false;

        for (int i = 0; i < mat0.rows; i++)
        {
            for (int j = 0; j < mat0.cols; j++)
            {
                if (mat0.at<float>(i, j) != mat1.at<float>(i, j))
                    return false;
            }
        }

        return true;
    }

    template<int ROWS, int COLS>
    inline cv::Matx44f Invert(const cv::Matx<float, ROWS, COLS>& transform)
    {
        static_assert((ROWS == 3 || ROWS == 4) && COLS == 4, "Invert only supports 3x4 or 4x4 matrices");

        //invert rotation by transpose on copy of upper 3x3 for inverse rotation
        cv::Matx44f invRotation{
            transform(0, 0), transform(1, 0), transform(2, 0), 0,
            transform(0, 1), transform(1, 1), transform(2, 1), 0,
            transform(0, 2), transform(1, 2), transform(2, 2), 0,
                          0,               0,               0, 1
        };

#ifndef NDEBUG
        // this approach would only work if the upper 3x3 was a pure rotation matrix (no scale)
        // that would mean all rows and columns were unit length and orthogonal. this tests for the former.
        for(int idx=0; idx < 3; idx++)
        {
            assert(abs(1.0 - cv::norm(invRotation.col(idx))) <= 0.001 && "Expecting unit length (no scale) for rotation matrix");
            assert(abs(1.0 - cv::norm(invRotation.row(idx))) <= 0.001 && "Expecting unit length (no scale) for rotation matrix");
        }

        //a valid rotation matrix has determinant 1 (a mirror would be -1)
        assert(abs(1.0 - cv::determinant(invRotation)) < 0.001  && "Expecting determinant of 1 for rotation matrix");
#endif

        //invert translation by negation on copy
        cv::Matx44f invTranslation{
            1, 0, 0, -transform(0, 3),
            0, 1, 0, -transform(1, 3),
            0, 0, 1, -transform(2, 3),
            0, 0, 0,               1
        };

        //inverse matrix
        return invRotation * invTranslation;
    }

    inline cv::Matx44f ComputeFrameTransform(const cv::Matx44f& from, const cv::Matx44f& to)
    {
        return to * Invert(from);
    }

    inline cv::Point3f UnProject(const cv::Matx33f& invCameraMatrix, const cv::Matx44f& viewMatrix, const cv::Point2i& point, float depth)
    {
        cv::Matx31f pixelSpace{ static_cast<float>(point.x), static_cast<float>(point.y), 1.f };
        cv::Matx31f cameraSpace = invCameraMatrix * pixelSpace;
        cameraSpace *= depth;

        cv::Matx41f world = Invert(viewMatrix) * cv::Matx41f{ cameraSpace(0), cameraSpace(1), cameraSpace(2), 1 };
        return{ world(0), world(1), world(2) };
    }

    const float NRM_EPSILON = 0.00001f;

    inline float Length(const cv::Vec3f& vec)
    {
        return std::sqrtf(vec.dot(vec));
    }

    inline cv::Vec3f Normalize(const cv::Vec3f& vec)
    {
        float distance = Length(vec);
        if (distance == 0)
            return vec;

        cv::Vec3f normVec = (vec / distance);
        assert(abs(sqrt(normVec.dot(normVec)) - 1.0f) < NRM_EPSILON);
        return normVec;
    }

    //rightHanded, column vec convention, intended for rotation of vectors (active)
    //https://en.wikipedia.org/wiki/Rotation_matrix
    //vs http://mathworld.wolfram.com/RotationMatrix.html intended for rotation of coordinate frames
    inline cv::Matx44f RotationXForVectors(float radians)
    {
        float cosR = std::cosf(radians);
        float sinR = std::sinf(radians);

        return{ 1.0f, 0.0f, 0.0f, 0.0f,
                 0.0f, cosR, -sinR, 0.0f,
                 0.0f, sinR, cosR, 0.0f,
                 0.0f, 0.0f, 0.0f, 1.0f };
    }

    //rightHanded, column vec convention, intended for rotation of vectors (active)
    //https://en.wikipedia.org/wiki/Rotation_matrix
    //vs http://mathworld.wolfram.com/RotationMatrix.html intended for rotation of coordinate frames
    inline cv::Matx44f RotationYForVectors(float radians)
    {
        float cosR = std::cosf(radians);
        float sinR = std::sinf(radians);

        return{ cosR, 0.0f, sinR, 0.0f,
                0.0f, 1.0f, 0.0f, 0.0f,
                -sinR, 0.0f, cosR, 0.0f,
                0.0f, 0.0f, 0.0f, 1.0f };
    }

    //rightHanded, column vec convention, intended for rotation of vectors (active)
    //https://en.wikipedia.org/wiki/Rotation_matrix
    //vs http://mathworld.wolfram.com/RotationMatrix.html intended for rotation of coordinate frames
    inline cv::Matx44f RotationZForVectors(float radians)
    {
        float cosR = std::cosf(radians);
        float sinR = std::sinf(radians);

        return{ cosR, -sinR, 0.0f, 0.0f,
                sinR, cosR, 0.0f, 0.0f,
                0.0f, 0.0f, 1.0f, 0.0f,
                0.0f, 0.0f, 0.0f, 1.0f };

    }

    //rightHanded, column vec convention, intended for rotation of coordinate frames (passive)
    //https://en.wikipedia.org/wiki/Rotation_matrix
    //ala http://mathworld.wolfram.com/RotationMatrix.html
    inline cv::Matx44f RotationXForCoordinateFrames(float radians)
    {
        float cosR = std::cosf(radians);
        float sinR = std::sinf(radians);

        return{ 1.0f, 0.0f, 0.0f, 0.0f,
            0.0f, cosR, sinR, 0.0f,
            0.0f, -sinR, cosR, 0.0f,
            0.0f, 0.0f, 0.0f, 1.0f };
    }

    //rightHanded, column vec convention,intended for rotation of coordinate frames (passive)
    //https://en.wikipedia.org/wiki/Rotation_matrix
    //ala http://mathworld.wolfram.com/RotationMatrix.html
    inline cv::Matx44f RotationYForCoordinateFrames(float radians)
    {
        float cosR = std::cosf(radians);
        float sinR = std::sinf(radians);

        return{ cosR, 0.0f, -sinR, 0.0f,
            0.0f, 1.0f, 0.0f, 0.0f,
            sinR, 0.0f, cosR, 0.0f,
            0.0f, 0.0f, 0.0f, 1.0f };
    }

    //rightHanded, column vec convention, intended for rotation of coordinate frames (passive)
    //https://en.wikipedia.org/wiki/Rotation_matrix
    //ala http://mathworld.wolfram.com/RotationMatrix.html
    inline cv::Matx44f RotationZForCoordinateFrames(float radians)
    {
        float cosR = std::cosf(radians);
        float sinR = std::sinf(radians);

        return{ cosR, sinR, 0.0f, 0.0f,
            -sinR, cosR, 0.0f, 0.0f,
            0.0f, 0.0f, 1.0f, 0.0f,
            0.0f, 0.0f, 0.0f, 1.0f };
    }

    // encoding of axis/angle in a vec3f
    // http://docs.opencv.org/2.4/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html#rodrigues
    inline cv::Matx33f RotationFromRodrigues(const cv::Matx31f& rotation)
    {
        cv::Matx33f rotationMat;
        cv::Rodrigues(rotation, rotationMat);
        return rotationMat;
    }

    template<typename T, size_t M, size_t N>
    inline std::array<T, M*N> ArrayFromMat(const cv::Matx<T, M, N>& mat)
    {
        std::array<T, M * N> out;
        std::copy(std::begin(mat.val), std::end(mat.val), out.begin());
        return out;
    }

    inline cv::Matx44f ToCVMat4x4(const mage::Matrix& matrix)
    {
        return cv::Matx44f(&matrix.M11);
    }

    inline bool IsEntirelyOffscreen(const cv::Rect& rect, const cv::Size& screenResPixels)
    {
        assert(rect.width > 0 && "invalid rect");
        assert(rect.height > 0 && "invalid rect");
        assert(screenResPixels.width > 0 && "invalid screenResPixels");
        assert(screenResPixels.height > 0 && "invalid screenResPixels");

        int maxHorizPixel = rect.width + rect.x - 1;
        int maxVerticalPixel = rect.height + rect.y - 1;

        bool offscreenX = (maxHorizPixel < 0) || (rect.x >(int)(screenResPixels.width - 1));
        bool offscreenY = (maxVerticalPixel < 0) || (rect.y >(int)(screenResPixels.height - 1));

        return (offscreenX || offscreenY);
    }

    inline bool IsWithinArea(const cv::Point2f& pt, const cv::Rect& crop)
    {
        assert(crop.width > 0 && "width must be positive nonzero");
        assert(crop.height > 0 && "height must be positive nonzero");

        return (pt.x >= crop.x && pt.x < (crop.x+crop.width)) && (pt.y >= crop.y && pt.y <= (crop.y + crop.height));
    }

    inline cv::Rect ToCVRect(const mage::Rect& rect)
    {
        return { (int)rect.X, (int)rect.Y, (int)rect.Width, (int)rect.Height };
    }

    cv::Mat CreateGrayCVMat(const cv::Size& resolution, const PixelFormat& format, const gsl::span<const uint8_t> imageBytes);
    cv::Mat CreateBGRCVMat(const cv::Mat& source, const PixelFormat& format);
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#pragma once
#include <map>
#include <vector>
#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#else
#include "FakeGSL.h"
#endif
#include "BooleanNodes.h"
#include "CommandLineArgs.h"

class Util {
public:
	static map<int, int> getNodeToValMap(map<int, int>& inputMap, vector<int>& inputs) {
		map<int, int> res;
		for(auto i = 0; i < inputs.size(); i++) {
			if (inputs[i] == 0 || inputs[i] == 1) {
				res[inputMap[i]] = inputs[i];
			}
		}
		return res;
	}

    static bool isAbsolute(bool_node* n) {
        if (n->type != bool_node::LT) {
            return false;
        }
        bool_node* v;
        if (n->mother()->type == bool_node::CONST && ((CONST_node*) (n->mother()))->getFval() == 0) {
            v = n->father();
        } else if (n->father()->type == bool_node::CONST && ((CONST_node*) (n->father()))->getFval() == 0) {
            v = n->mother();
        } else {
            return false;
        }
        FastSet<bool_node>& children = n->children;
        if (children.size() > 1) return false;
        for(child_iter it = children.begin(); it != children.end(); ++it) {
            if ((*it)->type == bool_node::ARRACC) {
                ARRACC_node* ac = (ARRACC_node*)(*it);
                if (ac->mother() != n) {
                    return false;
                }
                bool_node* m = ac->arguments(0);
                bool_node* f = ac->arguments(1);
                if (m != v && f != v) {
                    return false;
                }
                if (m == v) {
                    if (f->type != bool_node::NEG || f->mother() != v) {
                        return false;
                    }
                }
                if (f == v) {
                    if (m->type != bool_node::NEG || m->mother() != v) {
                        return false;
                    }
                }
            } else {
                return false;
            }
        }
        cout << "Absolute node: " << n->lprint() << endl;
        return true;
    }

    static bool isSqrt(bool_node* n) {
        if (n->type != bool_node::UFUN) {
            return false;
        }
        UFUN_node* un = (UFUN_node*) n;
        if (un->get_ufname() == "sqrt_math") {
            return true;
        }
        return false;
    }

    static bool hasArraccChild(bool_node* n) {
        FastSet<bool_node>& children = n->children;
        for(child_iter it = children.begin(); it != children.end(); ++it) {
            if ((*it)->type == bool_node::ARRACC && (*it)->mother() == n) {
                return true;
            }
        }
        return false;
    }

    static bool hasAssertChild(bool_node& n) {
        FastSet<bool_node>& children = n.children;
        if (children.size() > 1) return false;
        for(child_iter it = children.begin(); it != children.end(); ++it) {
            if ((*it)->type == bool_node::ASSERT) {
                return true;
            }
        }
        return false;
    }

    static bool hasNotAssertChild(bool_node& n) {
        FastSet<bool_node>& children = n.children;
        if (children.size() > 1) return false;
        for(child_iter it = children.begin(); it != children.end(); ++it) {
            if ((*it)->type == bool_node::NOT) {
                FastSet<bool_node>& grandChildren = (*it)->children;
                if (grandChildren.size() > 1) return false;
                for (child_iter it1 = grandChildren.begin(); it1 != grandChildren.end(); ++it1) {
                    if ((*it1)->type == bool_node::ASSERT) {
                        return true;
                    }
                }
            }
        }
        return false;
    }

    static bool hasAssertChild(bool_node* n) {
        return hasAssertChild(*n);
    }

    static bool hasNotAssertChild(bool_node* n) {
        return hasNotAssertChild(*n);
    }

    static set<int> getRelevantNodes(bool_node* n) {
        set<int> ids;
        set<int> visitedIds;
        vector<bool_node*> toVisit;
        toVisit.push_back(n);

        while(toVisit.size() > 0) {
            bool_node* node = toVisit.back();
            toVisit.pop_back();
            if (visitedIds.find(node->id) == visitedIds.end()) {
                visitedIds.insert(node->id);
                ids.insert(node->id);

				for (auto it = node->p_begin(); it != node->p_end(); ++it) {
                    toVisit.push_back(*it);
                }
            }
        }
        return ids;
    }

    static string print(const gsl_vector* v) {
        stringstream s;
        for (int i = 0; i < v->size; i++) {
            s << gsl_vector_get(v, i) << ";";
        }
        return s.str();
    }
    static string print(const gsl_vector* v, string delimiter) {
        stringstream s;
        for (int i = 0; i < v->size; i++) {
            s << gsl_vector_get(v, i) << delimiter;
        }
        return s.str();
    }

    static string print(double* arr, int len) {
        stringstream s;
        for (int i = 0; i < len; i++) {
            s << arr[i] << ";";
        }
        return s.str();
    }

    static string print(const vector<int>& v) {
        stringstream s;
        for (auto it = v.begin(); it != v.end(); it++) {
            s << *it << ";";
        }
        return s.str();
    }

    static string print(const vector<double>& v) {
        stringstream s;
        for (auto it = v.begin(); it != v.end(); it++) {
            s << *it << ";";
        }
        return s.str();
    }

    static string print(const set<int>& v) {
        stringstream s;
        for (auto it = v.begin(); it != v.end(); it++) {
            s << *it << ";";
        }
        return s.str();
    }

    static double getMin(const vector<double>& vals) {
        double minval = 1e30;
        for (int i = 0; i < vals.size(); i++) {
            if (vals[i] < minval) {
                minval = vals[i];
            }
        }
        return minval;
    }

    static double norm(const gsl_vector* v) {
        return gsl_blas_dnrm2(v);
    }

    static bool sameDir(const gsl_vector* v1, const gsl_vector* v2) {
        double dp = 0.0;
        gsl_blas_ddot(v1, v2, &dp);
        dp = dp/(norm(v1) * norm(v2));
        //cout << "dot product: " << dp << endl;
        if (dp > 0.9) return true;
        return false;
    }

    static string benchName() {
        string s = PARAMS->inputFname;
        int x1 = s.rfind(".sk");
        int x2 = s.rfind("/tmp/");
        return s.substr(x2+5, x1-x2 - 5);

    }

    static vector<string> split(string &text, string sep) {
        vector<string> tokens;
        size_t start = 0, end = 0;
        while ((end = text.find(sep, start)) != string::npos) {
            tokens.push_back(text.substr(start, end - start));
            start = end + 1;
        }
        tokens.push_back(text.substr(start));
        return tokens;
    }

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

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

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

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

/*
 * This function contains 62 FP additions, 38 FP multiplications,
 * (or, 60 additions, 36 multiplications, 2 fused multiply/add),
 * 24 stack variables, and 28 memory accesses
 */
static const fftw_real K445041867 = FFTW_KONST(+0.445041867912628808577805128993589518932711138);
static const fftw_real K1_801937735 = FFTW_KONST(+1.801937735804838252472204639014890102331838324);
static const fftw_real K1_246979603 = FFTW_KONST(+1.246979603717467061050009768008479621264549462);
static const fftw_real K867767478 = FFTW_KONST(+0.867767478235116240951536665696717509219981456);
static const fftw_real K1_949855824 = FFTW_KONST(+1.949855824363647214036263365987862434465571601);
static const fftw_real K1_563662964 = FFTW_KONST(+1.563662964936059617416889053348115500464669037);
static const fftw_real K2_000000000 = FFTW_KONST(+2.000000000000000000000000000000000000000000000);

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

void fftw_hc2real_14(const fftw_real *real_input, const fftw_real *imag_input, fftw_real *output, int real_istride, int imag_istride, int ostride)
{
     fftw_real tmp3;
     fftw_real tmp13;
     fftw_real tmp6;
     fftw_real tmp14;
     fftw_real tmp26;
     fftw_real tmp35;
     fftw_real tmp23;
     fftw_real tmp34;
     fftw_real tmp12;
     fftw_real tmp16;
     fftw_real tmp20;
     fftw_real tmp33;
     fftw_real tmp9;
     fftw_real tmp15;
     fftw_real tmp1;
     fftw_real tmp2;
     ASSERT_ALIGNED_DOUBLE;
     tmp1 = real_input[0];
     tmp2 = real_input[7 * real_istride];
     tmp3 = tmp1 - tmp2;
     tmp13 = tmp1 + tmp2;
     {
	  fftw_real tmp4;
	  fftw_real tmp5;
	  fftw_real tmp24;
	  fftw_real tmp25;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp4 = real_input[2 * real_istride];
	  tmp5 = real_input[5 * real_istride];
	  tmp6 = tmp4 - tmp5;
	  tmp14 = tmp4 + tmp5;
	  tmp24 = imag_input[2 * imag_istride];
	  tmp25 = imag_input[5 * imag_istride];
	  tmp26 = tmp24 - tmp25;
	  tmp35 = tmp24 + tmp25;
     }
     {
	  fftw_real tmp21;
	  fftw_real tmp22;
	  fftw_real tmp10;
	  fftw_real tmp11;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp21 = imag_input[6 * imag_istride];
	  tmp22 = imag_input[imag_istride];
	  tmp23 = tmp21 - tmp22;
	  tmp34 = tmp21 + tmp22;
	  tmp10 = real_input[6 * real_istride];
	  tmp11 = real_input[real_istride];
	  tmp12 = tmp10 - tmp11;
	  tmp16 = tmp10 + tmp11;
     }
     {
	  fftw_real tmp18;
	  fftw_real tmp19;
	  fftw_real tmp7;
	  fftw_real tmp8;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp18 = imag_input[4 * imag_istride];
	  tmp19 = imag_input[3 * imag_istride];
	  tmp20 = tmp18 - tmp19;
	  tmp33 = tmp18 + tmp19;
	  tmp7 = real_input[4 * real_istride];
	  tmp8 = real_input[3 * real_istride];
	  tmp9 = tmp7 - tmp8;
	  tmp15 = tmp7 + tmp8;
     }
     {
	  fftw_real tmp38;
	  fftw_real tmp37;
	  fftw_real tmp29;
	  fftw_real tmp28;
	  ASSERT_ALIGNED_DOUBLE;
	  output[7 * ostride] = tmp3 + (K2_000000000 * (tmp6 + tmp9 + tmp12));
	  tmp38 = (K1_563662964 * tmp35) + (K1_949855824 * tmp33) + (K867767478 * tmp34);
	  tmp37 = tmp3 + (K1_246979603 * tmp6) - (K1_801937735 * tmp12) - (K445041867 * tmp9);
	  output[ostride] = tmp37 - tmp38;
	  output[13 * ostride] = tmp37 + tmp38;
	  {
	       fftw_real tmp40;
	       fftw_real tmp39;
	       fftw_real tmp36;
	       fftw_real tmp32;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp40 = (K867767478 * tmp33) + (K1_563662964 * tmp34) - (K1_949855824 * tmp35);
	       tmp39 = tmp3 + (K1_246979603 * tmp12) - (K1_801937735 * tmp9) - (K445041867 * tmp6);
	       output[5 * ostride] = tmp39 - tmp40;
	       output[9 * ostride] = tmp39 + tmp40;
	       tmp36 = (K1_563662964 * tmp33) - (K1_949855824 * tmp34) - (K867767478 * tmp35);
	       tmp32 = tmp3 + (K1_246979603 * tmp9) - (K445041867 * tmp12) - (K1_801937735 * tmp6);
	       output[11 * ostride] = tmp32 - tmp36;
	       output[3 * ostride] = tmp32 + tmp36;
	  }
	  output[0] = tmp13 + (K2_000000000 * (tmp14 + tmp15 + tmp16));
	  tmp29 = (K867767478 * tmp20) + (K1_563662964 * tmp23) - (K1_949855824 * tmp26);
	  tmp28 = tmp13 + (K1_246979603 * tmp16) - (K1_801937735 * tmp15) - (K445041867 * tmp14);
	  output[12 * ostride] = tmp28 - tmp29;
	  output[2 * ostride] = tmp28 + tmp29;
	  {
	       fftw_real tmp31;
	       fftw_real tmp30;
	       fftw_real tmp27;
	       fftw_real tmp17;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp31 = (K1_563662964 * tmp26) + (K1_949855824 * tmp20) + (K867767478 * tmp23);
	       tmp30 = tmp13 + (K1_246979603 * tmp14) - (K1_801937735 * tmp16) - (K445041867 * tmp15);
	       output[8 * ostride] = tmp30 - tmp31;
	       output[6 * ostride] = tmp30 + tmp31;
	       tmp27 = (K1_563662964 * tmp20) - (K1_949855824 * tmp23) - (K867767478 * tmp26);
	       tmp17 = tmp13 + (K1_246979603 * tmp15) - (K445041867 * tmp16) - (K1_801937735 * tmp14);
	       output[4 * ostride] = tmp17 - tmp27;
	       output[10 * ostride] = tmp17 + tmp27;
	  }
     }
}

fftw_codelet_desc fftw_hc2real_14_desc =
{
     "fftw_hc2real_14",
     (void (*)()) fftw_hc2real_14,
     14,
     FFTW_BACKWARD,
     FFTW_HC2REAL,
     323,
     0,
     (const int *) 0,
};

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


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

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

/* set up parameters for this simulated annealing run */

#define N_TRIES 200		/* how many points do we try before stepping */
#define ITERS_FIXED_T 2000	/* how many iterations for each T? */
#define STEP_SIZE 1.0		/* max step size in random walk */
#define K 1.0			/* Boltzmann constant */
#define T_INITIAL 5000.0	/* initial temperature */
#define MU_T 1.002		/* damping factor for temperature */
#define T_MIN 5.0e-1

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

struct s_tsp_city {
  const char * name;
  double lat, longitude;	/* coordinates */
};
typedef struct s_tsp_city Stsp_city;

void prepare_distance_matrix(void);
void exhaustive_search(void);
void print_distance_matrix(void);
double city_distance(Stsp_city c1, Stsp_city c2);
double Etsp(void *xp);
double Mtsp(void *xp, void *yp);
void Stsp(const gsl_rng * r, void *xp, double step_size);
void Ptsp(void *xp);

/* in this table, latitude and longitude are obtained from the US
   Census Bureau, at http://www.census.gov/cgi-bin/gazetteer */

Stsp_city cities[] = {{"Santa Fe",    35.68,   105.95},
		      {"Phoenix",     33.54,   112.07},
		      {"Albuquerque", 35.12,   106.62},
		      {"Clovis",      34.41,   103.20},
		      {"Durango",     37.29,   107.87},
		      {"Dallas",      32.79,    96.77},
		      {"Tesuque",     35.77,   105.92},
		      {"Grants",      35.15,   107.84},
		      {"Los Alamos",  35.89,   106.28},
		      {"Las Cruces",  32.34,   106.76},
		      {"Cortez",      37.35,   108.58},
		      {"Gallup",      35.52,   108.74}};

#define N_CITIES (sizeof(cities)/sizeof(Stsp_city))

double distance_matrix[N_CITIES][N_CITIES];

/* distance between two cities */
double city_distance(Stsp_city c1, Stsp_city c2)
{
  const double earth_radius = 6375.000;	/* 6000KM approximately */
  /* sin and cos of lat and long; must convert to radians */
  double sla1 = sin(c1.lat*M_PI/180), cla1 = cos(c1.lat*M_PI/180),
    slo1 = sin(c1.longitude*M_PI/180), clo1 = cos(c1.longitude*M_PI/180);
  double sla2 = sin(c2.lat*M_PI/180), cla2 = cos(c2.lat*M_PI/180),
    slo2 = sin(c2.longitude*M_PI/180), clo2 = cos(c2.longitude*M_PI/180);

  double x1 = cla1*clo1;
  double x2 = cla2*clo2;

  double y1 = cla1*slo1;
  double y2 = cla2*slo2;

  double z1 = sla1;
  double z2 = sla2;

  double dot_product = x1*x2 + y1*y2 + z1*z2;

  double angle = acos(dot_product);

  /* distance is the angle (in radians) times the earth radius */
  return angle*earth_radius;
}

/* energy for the travelling salesman problem */
double Etsp(void *xp)
{
  /* an array of N_CITIES integers describing the order */
  int *route = (int *) xp;
  double E = 0;
  unsigned int i;

  for (i = 0; i < N_CITIES; ++i) {
    /* use the distance_matrix to optimize this calculation; it had
       better be allocated!! */
    E += distance_matrix[route[i]][route[(i + 1) % N_CITIES]];
  }

  return E;
}

double Mtsp(void *xp, void *yp)
{
  int *route1 = (int *) xp, *route2 = (int *) yp;
  double distance = 0;
  unsigned int i;

  for (i = 0; i < N_CITIES; ++i) {
    distance += ((route1[i] == route2[i]) ? 0 : 1);
  }

  return distance;
}

/* take a step through the TSP space */
void Stsp(const gsl_rng * r, void *xp, double step_size)
{
  int x1, x2, dummy;
  int *route = (int *) xp;

  step_size = 0 ; /* prevent warnings about unused parameter */

  /* pick the two cities to swap in the matrix; we leave the first
     city fixed */
  x1 = (gsl_rng_get (r) % (N_CITIES-1)) + 1;
  do {
    x2 = (gsl_rng_get (r) % (N_CITIES-1)) + 1;
  } while (x2 == x1);

  dummy = route[x1];
  route[x1] = route[x2];
  route[x2] = dummy;
}

void Ptsp(void *xp)
{
  unsigned int i;
  int *route = (int *) xp;
  printf("  [");
  for (i = 0; i < N_CITIES; ++i) {
    printf(" %d ", route[i]);
  }
  printf("]  ");
}

int main(void)
{
  int x_initial[N_CITIES];
  unsigned int i;

  const gsl_rng * r = gsl_rng_alloc (gsl_rng_env_setup()) ;

  gsl_ieee_env_setup ();

  prepare_distance_matrix();

  /* set up a trivial initial route */
  printf("# initial order of cities:\n");
  for (i = 0; i < N_CITIES; ++i) {
    printf("# \"%s\"\n", cities[i].name);
    x_initial[i] = i;
  }

  printf("# distance matrix is:\n");
  print_distance_matrix();

  printf("# initial coordinates of cities (longitude and latitude)\n");
  /* this can be plotted with */
  /* ./siman_tsp > hhh ; grep city_coord hhh | awk '{print $2 "   " $3}' | xyplot -ps -d "xy" > c.eps */
  for (i = 0; i < N_CITIES+1; ++i) {
    printf("###initial_city_coord: %g %g \"%s\"\n",
	   -cities[x_initial[i % N_CITIES]].longitude,
	   cities[x_initial[i % N_CITIES]].lat,
	   cities[x_initial[i % N_CITIES]].name);
  }

/*   exhaustive_search(); */

  gsl_siman_solve(r, x_initial, Etsp, Stsp, Mtsp, Ptsp, NULL, NULL, NULL,
		  N_CITIES*sizeof(int), params);

  printf("# final order of cities:\n");
  for (i = 0; i < N_CITIES; ++i) {
    printf("# \"%s\"\n", cities[x_initial[i]].name);
  }

  printf("# final coordinates of cities (longitude and latitude)\n");
  /* this can be plotted with */
  /* ./siman_tsp > hhh ; grep city_coord hhh | awk '{print $2 "   " $3}' | xyplot -ps -d "xy" > c.eps */
  for (i = 0; i < N_CITIES+1; ++i) {
    printf("###final_city_coord: %g %g %s\n",
	   -cities[x_initial[i % N_CITIES]].longitude,
	   cities[x_initial[i % N_CITIES]].lat,
	   cities[x_initial[i % N_CITIES]].name);
  }

  printf("# ");
  fflush(stdout);
#if 0
  system("date");
#endif /* 0 */
  fflush(stdout);

  return 0;
}

void prepare_distance_matrix()
{
  unsigned int i, j;
  double dist;

  for (i = 0; i < N_CITIES; ++i) {
    for (j = 0; j < N_CITIES; ++j) {
      if (i == j) {
	dist = 0;
      } else {
	dist = city_distance(cities[i], cities[j]);
      }
      distance_matrix[i][j] = dist;
    }
  }
}

void print_distance_matrix()
{
  unsigned int i, j;

  for (i = 0; i < N_CITIES; ++i) {
    printf("# ");
    for (j = 0; j < N_CITIES; ++j) {
      printf("%15.8f   ", distance_matrix[i][j]);
    }
    printf("\n");
  }
}

/* [only works for 12] search the entire space for solutions */
static double best_E = 1.0e100, second_E = 1.0e100, third_E = 1.0e100;
static int best_route[N_CITIES];
static int second_route[N_CITIES];
static int third_route[N_CITIES];
static void do_all_perms(int *route, int n);

void exhaustive_search()
{
  static int initial_route[N_CITIES] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11};
  printf("\n# ");
  fflush(stdout);
#if 0
  system("date");
#endif
  fflush(stdout);
  do_all_perms(initial_route, 1);
  printf("\n# ");
  fflush(stdout);
#if 0
  system("date");
#endif /* 0 */
  fflush(stdout);

  printf("# exhaustive best route: ");
  Ptsp(best_route);
  printf("\n# its energy is: %g\n", best_E);

  printf("# exhaustive second_best route: ");
  Ptsp(second_route);
  printf("\n# its energy is: %g\n", second_E);

  printf("# exhaustive third_best route: ");
  Ptsp(third_route);
  printf("\n# its energy is: %g\n", third_E);
}

/* James Theiler's recursive algorithm for generating all routes */
static void do_all_perms(int *route, int n)
{
  if (n == (N_CITIES-1)) {
    /* do it! calculate the energy/cost for that route */
    double E;
    E = Etsp(route);		/* TSP energy function */
    /* now save the best 3 energies and routes */
    if (E < best_E) {
      third_E = second_E;
      memcpy(third_route, second_route, N_CITIES*sizeof(*route));
      second_E = best_E;
      memcpy(second_route, best_route, N_CITIES*sizeof(*route));
      best_E = E;
      memcpy(best_route, route, N_CITIES*sizeof(*route));
    } else if (E < second_E) {
      third_E = second_E;
      memcpy(third_route, second_route, N_CITIES*sizeof(*route));
      second_E = E;
      memcpy(second_route, route, N_CITIES*sizeof(*route));
    } else if (E < third_E) {
      third_E = E;
      memcpy(route, third_route, N_CITIES*sizeof(*route));
    }
  } else {
    int new_route[N_CITIES];
    unsigned int j;
    int swap_tmp;
    memcpy(new_route, route, N_CITIES*sizeof(*route));
    for (j = n; j < N_CITIES; ++j) {
      swap_tmp = new_route[j];
      new_route[j] = new_route[n];
      new_route[n] = swap_tmp;
      do_all_perms(new_route, n+1);
    }
  }
}