text:
#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define MAX(a, b) ((a) > (b) ? (a) : (b))

#define DefaultWorkSize 1e9 // 1GB
#define NumSplitDefaultMax 20

#define OUTPUT stdout // stderr

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <typeinfo>

#define __STDC_FORMAT_MACROS
#include <inttypes.h>
#include <sys/time.h>
#include <float.h>

#if defined (MKL)
#include <mkl_cblas.h>
#include <mkl_trans.h>
#include <mkl_spblas.h>
#else
#include <cblas.h>
#endif

#include <omp.h>

#if defined (FLOAT128)
#if defined (ARM)
#define __float128 long double
#define FLT128_MAX LDBL_MAX
#define FLT128_MIN LDBL_MIN
#else
#include <quadmath.h>
#endif
//#ifdef __INTEL_COMPILER
//__float128 std::abs (const __float128 x);
//#endif
#include <mplapack/mpblas__Float128.h>
#endif

#if defined (CUBLAS)
#include <cuda.h>
#include <cuda_runtime_api.h>
#include <cublas_v2.h>
#endif

#include "../../include/ozblas.h"
#include "ozblas_internal.h"

// -------------------------------------
// constexpr functions
// -------------------------------------

// -------------------------------------
// getEmin
// -------------------------------------
template <typename TYPE>
constexpr int32_t getEmin () {
	fprintf (OUTPUT, "OzBLAS error: TYPE is not specified in getEmin.\n");
	exit (1);
	return 0;
}
template <>
constexpr int32_t getEmin <float> () {
	return -126;
}
template <>
constexpr int32_t getEmin <double> () {
	return -1022;
}
#if defined (FLOAT128)
template <>
constexpr int32_t getEmin <__float128> () {
	return -16382;
}
#endif

// -------------------------------------
// getEpse
// -------------------------------------
template <typename TYPE>
constexpr int32_t getEpse () {
	fprintf (OUTPUT, "OzBLAS error: TYPE is not specified in getEpse.\n");
	exit (1);
	return 0;
}
template <>
constexpr int32_t getEpse <float> () {
	return 24;
}
template <>
constexpr int32_t getEpse <double> () {
	return 53;
}
#if defined (FLOAT128)
template <>
constexpr int32_t getEpse <__float128> () {
	return 113;
}
#endif

// -------------------------------------
// getTypeMax
// -------------------------------------
template <typename TYPE>
constexpr TYPE getTypeMax () {
	fprintf (OUTPUT, "OzBLAS error: TYPE is not specified in getTypeMax.\n");
	exit (1);
	return 0;
}
template <>
constexpr float getTypeMax <float> () {
	return FLT_MAX;
}
template <>
constexpr double getTypeMax <double> () {
	return DBL_MAX;
}
#if defined (FLOAT128)
template <>
constexpr __float128 getTypeMax <__float128> () {
	return FLT128_MAX;
}
#endif

// -------------------------------------
// getTypeMin
// -------------------------------------
template <typename TYPE>
constexpr TYPE getTypeMin () {
	fprintf (OUTPUT, "OzBLAS error: TYPE is not specified in getTypeMin.\n");
	exit (1);
	return 0;
}
template <>
constexpr float getTypeMin <float> () {
	return FLT_MIN;
}
template <>
constexpr double getTypeMin <double> () {
	return DBL_MIN;
}
#if defined (FLOAT128)
template <>
constexpr __float128 getTypeMin <__float128> () {
	return FLT128_MIN;
}
#endif

/*
#define SCALBN scalbnl
#define LOG2 log2l
#define LOG logl
#define CEIL ceill
#define FABS fabsl
*/

max_stars_repo_name: wsmoses/Riken-blas

max_issues_repo_name: wsmoses/Riken-blas

max_forks_repo_name: wsmoses/Riken-blas


text:
////////////////////////////////////////////////////////////
//
// Copyright (c) 2018 Jan Filipowicz, Filip Turobos
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//
////////////////////////////////////////////////////////////

#ifndef GENETIC_ALGORITHM_LIBRARY_GENETIC_ALGORITHM_H
#define GENETIC_ALGORITHM_LIBRARY_GENETIC_ALGORITHM_H

#include <algorithm>
#include <functional>
#include <iterator>
#include <utility>
#include <vector>
#include <gsl/gsl_assert>
#include "evaluated_specimen.h"
#include "repeat.h"

template<class Specimen, class Rating>
class genetic_algorithm {
public:
	using specimen_type = Specimen;
	using rating_type = Rating;
	using evaluated_specimen_type = evaluated_specimen<Specimen, Rating>;
	struct context_type;
	enum struct stage_type { generated, selected, bred };
	explicit genetic_algorithm(const context_type& context);
	explicit genetic_algorithm(context_type&& context);
	template<class... Functions>
	evaluated_specimen_type operator()(Functions&&... observers) const;
private:
	template<class... Functions>
	void communicate_stage(stage_type stage, const std::vector<evaluated_specimen_type>& specimens, Functions&&... observers) const;
	void evaluate(std::vector<evaluated_specimen_type>& specimens) const;
	std::vector<evaluated_specimen_type> breed(const std::vector<evaluated_specimen_type>& specimens) const;
	context_type context;
};

template<class Specimen, class Rating>
struct genetic_algorithm<Specimen, Rating>::context_type {
	std::size_t initial_population_size = 1;
	std::size_t breeding_population_size = 1;
	std::size_t max_iterations = 0;
	std::function<specimen_type()> generator;
	std::function<rating_type(const specimen_type&)> evaluator;
	std::function<void(std::vector<evaluated_specimen_type>&, std::size_t)> selector;
	std::function<specimen_type(const specimen_type&, const specimen_type&)> breeder;
	std::function<bool(rating_type, rating_type)> comparator;
};

template<class Context>
constexpr bool valid(const Context& context) noexcept {
	return context.initial_population_size > 0
		&& context.breeding_population_size > 0
		&& context.generator != nullptr
		&& context.evaluator != nullptr
		&& context.selector != nullptr
		&& context.breeder != nullptr
		&& context.comparator != nullptr
	;
}

template<class Specimen, class Rating>
inline genetic_algorithm<Specimen, Rating>::genetic_algorithm(const context_type& context)
	: context(context) {}

template<class Specimen, class Rating>
inline genetic_algorithm<Specimen, Rating>::genetic_algorithm(context_type&& context)
	: context(std::move(context)) {}

template<class Specimen, class Rating>
template<class... Functions>
inline auto genetic_algorithm<Specimen, Rating>::operator()(Functions&&... observers) const -> evaluated_specimen_type {
	Expects(valid(context));
	std::vector<evaluated_specimen_type> specimens(context.initial_population_size);
	for (auto&& specimen : specimens) {
		specimen.value() = context.generator();
	}
	evaluate(specimens);
	communicate_stage(stage_type::generated, specimens, std::forward<Functions>(observers)...);
	repeat(context.max_iterations, [&] {
		context.selector(specimens, context.breeding_population_size);
		communicate_stage(stage_type::selected, specimens, std::forward<Functions>(observers)...);
		specimens = breed(specimens);
		evaluate(specimens);
		communicate_stage(stage_type::bred, specimens, std::forward<Functions>(observers)...);
	});
	return std::move(*std::max_element(specimens.begin(), specimens.end(), [this](const evaluated_specimen_type& lhs, const evaluated_specimen_type& rhs) {
		return context.comparator(lhs.rating(), rhs.rating());
	}));
}

template<class Specimen, class Rating>
template<class... Functions>
inline void genetic_algorithm<Specimen, Rating>::communicate_stage(stage_type stage, const std::vector<evaluated_specimen_type>& specimens, Functions&&... observers) const {
	(void)stage;
	(void)specimens;
	((void)observers(*this, stage, specimens), ...);
}

template<class Specimen, class Rating>
inline void genetic_algorithm<Specimen, Rating>::evaluate(std::vector<evaluated_specimen_type>& specimens) const {
	for (auto&& specimen : specimens) {
		specimen.evaluate(context.evaluator);
	}
}

template<class Specimen, class Rating>
inline auto genetic_algorithm<Specimen, Rating>::breed(const std::vector<evaluated_specimen_type>& specimens) const -> std::vector<evaluated_specimen_type> {
	const std::size_t specimen_count = specimens.size();
	std::vector<evaluated_specimen_type> result;
	const std::size_t reserve = specimen_count * (specimen_count - 1) / 2;
	result.reserve(reserve);
	for (auto it = specimens.begin(); it != specimens.end(); ++it) {
		const specimen_type& father = it->value();
		std::for_each(std::next(it), specimens.end(), [&](const evaluated_specimen_type& mother) {
			result.emplace_back(evaluated_specimen_type {context.breeder(father, mother.value())});
		});
	}
	Ensures(result.size() == reserve);
	return result;
}

#endif

ext: h

max_stars_repo_name: SirEmentaler/Eugenics-Wars

max_issues_repo_name: SirEmentaler/Eugenics-Wars

max_forks_repo_name: SirEmentaler/Eugenics-Wars


text:
#ifndef PyGSL_BLOCK_HELPERS_H
#define PyGSL_BLOCK_HELPERS_H 1
/*
 * Author:  Pierre Schnizer
 *
 *
 */

/*****************************************************************************
 * The functions in this file are implemented in src/init/block_helpers*     *
 * These functions try to simplify the conversion between the PyArrayObjects *
 * and the gsl_vector and gsl_matrix types. Further they try to provide a    *
 * compatibility layer for the different implementations of the              *
 * PyArrayObjects (Numeric, numarray, numpy). All code using the functions   *
 * must:                                                                     *
 *       Use these functions here :-)					     *
 *       Declare all varibales referencing dimensions and indexing arrays    *
 *       "PyGSL_array_index_t"						     *
 *       Expect these variables to be either of type int or long. This most  *
 *       probable is only noticable on 64bit architectures.		     *
 *****************************************************************************/
/*
#define PY_ARRAY_UNIQUE_SYMBOL PyGSL_PY_ARRAY_API
*/

#ifndef _PyGSL_API_MODULE
/*
 * Not possible to just import the declaration of the PyArrayObject struct. The
 * pointer comes along as well ...
 *
 */
#ifndef PyGSL_IMPORT_ARRAY
#define NO_IMPORT_ARRAY
#endif
#endif  /* _PyGSL_API_MODULE */
#include <pygsl/intern.h>
#include <pygsl/utils.h>
#include <pygsl/general_helpers.h>
#include <limits.h>
/*
 * Build array info out of various parts.
 *
 * array_flag
 */

#if ULONG_MAX < 0xffffff
#error "unsigned long not big enough for array info"
#else
typedef unsigned long PyGSL_array_info_t;
#endif

/* typedef Py_ssize_t PyGSL_array_index_t; */
/*
#if PY_VERSION_HEX < 0x02050000
typedef int Py_ssize_t;
#endif
*/
#include <pygsl/arrayobject.h>
#ifdef PyGSL_NUMPY
#include <numpy/arrayobject.h>
#endif
#ifdef PyGSL_NUMERIC
#include <Numeric/arrayobject.h>
#endif
#ifdef PyGSL_NUMARRAY
#include <numarray/libnumarray.h>
#endif
#if (!defined PyGSL_NUMPY) && (!defined  PyGSL_NUMERIC) && (!defined PyGSL_NUMARRAY)
#error "Neither numpy nor numarray nor Numeric is defined!"
#endif
/*
 *  required for 64 bit machines; backward compability. I rather perfer
 *  to declare it as a long as its just my functions.
 *
 *  For now I keep it as the same type as the underlaying array package
 *  uses for its conversions.
 */
#ifdef PyGSL_NUMERIC
typedef int PyGSL_array_index_t;
#endif
#ifdef PyGSL_NUMARRAY
typedef maybelong PyGSL_array_index_t;
#endif
#ifdef PyGSL_NUMPY
/* typedef intp PyGSL_array_index_t; */
typedef npy_intp PyGSL_array_index_t;
#endif

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

enum PyGSL_Array_Flags {
     PyGSL_NON_CONTIGUOUS = 0,
     PyGSL_CONTIGUOUS = 1,
     /* Additional flags needed for numarray and numpy */
     PyGSL_INPUT_ARRAY = 2,
     PyGSL_OUTPUT_ARRAY = 4,
     PyGSL_IO_ARRAY = 8,
};

/*
 * Convienience functions for handling the flags
 */
#define PyGSL_GET_ARRAYFLAG(flag) (( (flag) >>  0) & 0x000000ff)
#define PyGSL_GET_ARRAYTYPE(flag) (( (flag) >>  8) & 0x000000ff)
#define PyGSL_GET_TYPESIZE(flag)  (( (flag) >> 16) & 0x000000ff)
#define PyGSL_GET_ARGNUM(flag)    (( (flag) >> 24) & 0x000000ff)

#define PyGSL_BUILD_ARRAY_INFO(array_flag, array_type, type_size, argnum) \
(\
 (  ((array_flag) & 0x000000ff) <<  0) + \
 (  ((array_type) & 0x000000ff) <<  8) + \
 (  ((type_size)  & 0x000000ff) << 16) + \
 (  ((argnum)     & 0x000000ff) << 24)   \
)

#define PyGSL_DARRAY_INFO(array_flag, argnum)  \
   PyGSL_BUILD_ARRAY_INFO(array_flag, PyArray_DOUBLE, sizeof(double), argnum)
#define PyGSL_DARRAY_CINPUT(argnum)  PyGSL_DARRAY_INFO(PyGSL_CONTIGUOUS     | PyGSL_INPUT_ARRAY, argnum)
#define PyGSL_DARRAY_INPUT(argnum)   PyGSL_DARRAY_INFO(PyGSL_NON_CONTIGUOUS | PyGSL_INPUT_ARRAY, argnum)

/*
 * PyGSL_New_Array:
 *                Generate an new array with the specified dimensions.
 *
 *                The numpy backend expects an array of int, whereas
 *                the numarray backend expects  an array of long
 */
PyGSL_API_EXTERN PyArrayObject *
PyGSL_New_Array(int nd,  PyGSL_array_index_t *dimensions, int type);

/*
 * PyGSL_Copy_Array:
 *
 *                 Copy an array. The output array will have the same
 *                 size as the input array.
 */
PyGSL_API_EXTERN PyArrayObject *
PyGSL_Copy_Array(PyArrayObject *, int type);

/*
 * PyGSL_STRIDE_RECALC:
 *                                Recalc a stride and check if it is okay
 *
 * Numpy calculates strides in bytes, gsl as multiple of the basis type size
 *
 * Return value:
 *         -1: Conversion Failed
 *         pos : recalculated stride
 */
PyGSL_API_EXTERN int
PyGSL_stride_recalc(PyGSL_array_index_t strides, int basis_type_size,
		    PyGSL_array_index_t * stride_recalc);

/*
 * PyGSL_PyArray_generate_gsl_vector_view :
 *                                 Generate an new array of given dimensions .
 *
 *  This function will try to convert the object, to a Python integer and
 *  generate an apropriate one dimensional numpy array.
 *
 * Input :
 *         object              : a general python object
 *         array_type          : the required C type for the array
 *         argument number     : The argument number. Used for error reporting
 *
 * Output:
 *                             : a pointer to a PyArrayObject or NULL in case
 *                              of error. This object must be dereferenced.
 */
PyGSL_API_EXTERN PyArrayObject *
PyGSL_PyArray_generate_gsl_vector_view(PyObject *src,
				       int array_type,
				       int argnum);

/*
 * PyGSL_PyArray_generate_gsl_matrix_view :
 *                                 Generate an new array of given dimensions .
 *
 *  This function will try to convert the object, to a sequence of two Python
 *  integer and generate an apropriate two dimensional numpy array.
 *
 * Input :
 *         object              : a general python object
 *         array_type          : the required C type for the array
 *         argument number     : The argument number. Used for error reporting
 *
 * Output:
 *                             : a pointer to a PyArrayObject or NULL in case
 *                              of error. This object must be dereferenced.
 */
PyGSL_API_EXTERN PyArrayObject *
PyGSL_PyArray_generate_gsl_matrix_view(PyObject *src,
				       int array_type,
				       int argnum);

/*
 * Check if a python object can be used as a vector. If so return an
 *  approbriate python array object.
 *
 *
 * Input:
 *      src        ... the python object
 *      size       ... the required size. Pass -1 if it should not be checked.
 *
 *      array_type ... the required type.
 *      flag       ... shall be an input or output or inout array?
 *                 ... shall it be contiguous?
 *      argnum     ... the positional argument number. Used when reporting an
 *                     error
 *
 *                     Information if the stride of the vector should be
 *                     recalulated to the native C type.
 *      type_size  ... the size of the native C type (typically sizeof(type)).
 *                     If stride is NULL this value is without meaning.
 *
 *      *stride    ... the recalculated stride. If NULL stride will not be
 *                      recalulated, and the following items can be of any value
 *
 *      info       ... info structure passed by callbacks for reporting an error
 *
 *     Convienience macros:
 *                PyGSL_DVECTOR_CHECK(src, size, flag, argnum, stride, info)
 */
PyGSL_API_EXTERN PyArrayObject *
PyGSL_vector_check(PyObject *src, PyGSL_array_index_t size, PyGSL_array_info_t ainfo,
		   PyGSL_array_index_t *stride, PyGSL_error_info * info);

/*
 * Check if a python object can be used as a vector. If so return an
 * approbriate python array object.
 *
 *
 * Input:
 *      src        ... the python object
 *
 *      size1      ... the required size of the first dimension.  Pass -1 if it
 *                     should not be checked.
 *      size2     ...  the required size of the second dimension. Pass -1 if it
 *                     should not be checked.
 *
 *      array_type ... the required type.
 *      flag       ... shall be an input or output or inout array?
 *                     shall it be contiguous?
 *      argnum     ... the positional argument number. Used when reporting an
 *                     error
 *
 *                     Information if the stride of the vector should be
 *                     recalulated to the native C type.
 *      type_size  ... the size of the native C type (typically sizeof(type)).
 *                     If stride is NULL this value is without meaning.
 *
 *      *stride1   ... the recalculated stride. If NULL stride will not be
 *                      recalulated, and the following items can be of any value
 *      *stride2   ... the recalculated stride. If NULL stride will not be
 *                      recalulated, and the following items can be of any value
 *
 *      info       ... info structure passed by callbacks for reporting an error
 *
 *     Convienience macros:
 *       PyGSL_DMATRIX_CHECK(src, size1, size2, flag, argnum, stride1, stride2, info)
 */
PyGSL_API_EXTERN PyArrayObject *
PyGSL_matrix_check(PyObject *src, PyGSL_array_index_t size1, PyGSL_array_index_t size2,
		   PyGSL_array_info_t ainfo,
		   PyGSL_array_index_t *stride1, PyGSL_array_index_t *stride2,
		   PyGSL_error_info * info);

/*
 * PyGSL_copy_pyarray_to_gslvector
 *                           Copy the contents of a pyarray to a gsl_vector.
 *
 *
 * Input :
 *         f                   : a pointer to the target vector
 *         object              : a general python object
 *         n                   : number of elements
 *         info                : a PyGSL_Error_info struct. Used for error
 *                               reporting during callback evaluation. Pass
 *                               NULL if not needed.
 *
 * Output:
 *                             : GSL_SUCCESS | GSL_FAILURE
 *
 */
PyGSL_API_EXTERN int
PyGSL_copy_pyarray_to_gslvector(gsl_vector *f, PyObject *object, PyGSL_array_index_t n,
				PyGSL_error_info * info);

/*
 * PyGSL_copy_pyarray_to_gslmatrix
 *                           Copy the contents of a pyarray to a gsl_matrix.
 *
 *
 * Input :
 *         f                   : a pointer to the target gsl_vector
 *         object              : a general python object referring to a numpy
 *                               array
 *         n                   : number of elements in the first dimension
 *         p                   : number of elements in the second dimension
 *         info                : a PyGSL_error_info struct. Used for error
 *                               reporting during callback evaluation. Pass
 *                               NULL if not needed.
 *
 * Output:
 *                             : GSL_SUCCESS | GSL_FAILURE
 *
 */
PyGSL_API_EXTERN int
PyGSL_copy_pyarray_to_gslmatrix(gsl_matrix *f, PyObject *object,  PyGSL_array_index_t n,
				PyGSL_array_index_t p, PyGSL_error_info * info);

PyGSL_API_EXTERN PyArrayObject *
PyGSL_vector_or_double(PyObject *src, PyGSL_array_info_t ainfo, PyGSL_array_index_t size,
		       PyGSL_error_info * info);

/*
 * PyGSL_copy_gslvector_to_pyarrary :
 *                Generate a new numpy array of approbriate size and copy the
 *                data of the vector to it.
 *
 *
 * Input :
 *              x              : a gsl_vector
 * Output:
 *                             : a pointer to a PyArrayObject or NULL in case
 *                              of error. This object must be dereferenced.
 */
PyGSL_API_EXTERN PyArrayObject *
PyGSL_copy_gslvector_to_pyarray(const gsl_vector *x);

/*
 * PyGSL_copy_gslmatrix_to_pyarrary :
 *                Generate a new numpy array of approbriate size and copy the
 *                data of the matrix to it.
 *
 *
 * Input :
 *              x              : a gsl_matrix
 * Output:
 *                             : a pointer to a PyArrayObject or NULL in case
 *                              of error. This object must be dereferenced.
 */
PyGSL_API_EXTERN PyArrayObject *
PyGSL_copy_gslmatrix_to_pyarray(const gsl_matrix *x);
#ifndef _PyGSL_API_MODULE
#define PyGSL_stride_recalc \
(*(int (*)(PyGSL_array_index_t, int, PyGSL_array_index_t * ))                                        PyGSL_API[PyGSL_stride_recalc_NUM])
#define  PyGSL_New_Array  \
(*(PyArrayObject * (*)(int, PyGSL_array_index_t *, int))                              PyGSL_API[PyGSL_PyArray_new_NUM])
#define  PyGSL_Copy_Array  \
(*(PyArrayObject * (*)(PyArrayObject *))                               PyGSL_API[PyGSL_PyArray_copy_NUM])
#ifdef NO_PyGSL_LEGACY
#define  PyGSL_PyArray_prepare_gsl_vector_view  \
(*(PyArrayObject * (*)(PyObject *, int, int,  PyGSL_array_index_t, int, PyGSL_error_info *)) \
                                                                     PyGSL_API[PyGSL_PyArray_prepare_gsl_vector_view_NUM])
#define  PyGSL_PyArray_prepare_gsl_matrix_view  \
(*(PyArrayObject * (*)(PyObject *, int, int,  PyGSL_array_index_t, PyGSL_array_index_t, int, PyGSL_error_info *)) \
                                                                     PyGSL_API[PyGSL_PyArray_prepare_gsl_matrix_view_NUM])
#endif /* NO_PyGSL_LEGACY */
#define PyGSL_PyArray_generate_gsl_vector_view \
(*(PyArrayObject *(*)(PyObject *, int, int))          PyGSL_API[PyGSL_PyArray_generate_gsl_vector_view_NUM])

#define PyGSL_PyArray_generate_gsl_matrix_view \
(*(PyArrayObject *(*)(PyObject *, int, int))          PyGSL_API[PyGSL_PyArray_generate_gsl_matrix_view_NUM])

#define PyGSL_copy_pyarray_to_gslvector \
(*(int (*) (gsl_vector *, PyObject *, PyGSL_array_index_t, PyGSL_error_info *))      PyGSL_API[PyGSL_copy_pyarray_to_gslvector_NUM])
#define PyGSL_copy_pyarray_to_gslmatrix \
(*(int (*) (gsl_matrix *, PyObject *, PyGSL_array_index_t, PyGSL_array_index_t, PyGSL_error_info *)) PyGSL_API[PyGSL_copy_pyarray_to_gslmatrix_NUM])

#define PyGSL_copy_gslvector_to_pyarray \
 (*(PyArrayObject * (*)(const gsl_vector *))                         PyGSL_API[ PyGSL_copy_gslvector_to_pyarray_NUM])

#define PyGSL_copy_gslmatrix_to_pyarray \
 (*(PyArrayObject * (*)(const gsl_matrix *))                         PyGSL_API[ PyGSL_copy_gslmatrix_to_pyarray_NUM])

#define  PyGSL_vector_or_double  \
(*(PyArrayObject * (*)(PyObject *, PyGSL_array_info_t , PyGSL_array_index_t, PyGSL_error_info *)) \
                                                                     PyGSL_API[PyGSL_vector_or_double_NUM])

#define PyGSL_vector_check  \
(*(PyArrayObject * (*) (PyObject *, PyGSL_array_index_t, \
 PyGSL_array_info_t, PyGSL_array_index_t *, PyGSL_error_info *)) \
                                                                      PyGSL_API[PyGSL_vector_check_NUM])

#define PyGSL_matrix_check  \
(*(PyArrayObject * (*) (PyObject *, PyGSL_array_index_t, PyGSL_array_index_t, \
 PyGSL_array_info_t, PyGSL_array_index_t *, PyGSL_array_index_t *, PyGSL_error_info *)) \
                                                                      PyGSL_API[PyGSL_matrix_check_NUM])

#define PyGSL_array_check (*(int (*)(PyObject *))  PyGSL_API[PyGSL_array_check_NUM])
#define PyGSL_array_return(ob) ((PyObject *) ob)
#endif /* _PyGSL_API_MODULE */

#define PyGSL_STRIDE_RECALC(strides, basis_type_size, stride_recalc) \
      ( \
           (((strides) % (basis_type_size)) == 0) \
         ? \
           ( ((*(stride_recalc)) = (strides) / (basis_type_size)), GSL_SUCCESS ) \
         : \
           PyGSL_stride_recalc(strides, basis_type_size, stride_recalc) \
     )
/*
#define PyGSL_DVECTOR_CHECK(src, size, info, stride, info) \
PyGSL_vector_check(src, size, PyArray_DOUBLE, sizeof(double), flag, argnum, stride, info)

#define PyGSL_DMATRIX_CHECK(src, size1, size2, flag, argnum, stride1, stride2, info) \
PyGSL_vector_check(src, size1, size2, PyArray_DOUBLE, sizeof(double), flag, argnum, stride1, stride2, info)
*/

#ifdef PyGSL_NUMPY
#include <pygsl/block_helpers_numpy.h>
#endif

#ifdef PyGSL_NUMERIC
#include <pygsl/block_helpers_numeric.h>
#endif

#ifdef PyGSL_NUMARRAY
#include <pygsl/block_helpers_numarray.h>
#endif

#endif /* PyGSL_BLOCK_HELPERS_H */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#ifndef SVD_H
#define SVD_H

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

/**
 * Computes, stores, and maintains a singular value decomposition of a gsl_matrix.
 *
 * A matrix-pointer is passed to the SVD on construction. This pointer is kept through
 * the life-time of the SVD object (but ownership is not assumed). If the contents
 * of the matrix is changed, the SVD can be update to reflect the change using the
 * abstract UpdateFromMatrix function. This is not automatically done on construction.
 */
class SVD{
 protected:
  const int m, n; ///< Dimensions of matrix

  /** Decomposes M into U*S*V^t.  */
  SVD(gsl_matrix* M);
 public:

  gsl_matrix * const matrix;  //TODO: Make private
  gsl_matrix * const U;       //TODO: Make private
  gsl_matrix * const V;       //TODO: Make private
  gsl_vector * const S;       //TODO: Make private

  /** Constructs an SVD object. If MKL is available it will be an MKLSVD and if GSL
   * is available it will be GSLSVD
   */
  static SVD* createSVD(gsl_matrix* M);

  virtual ~SVD();

  /** Update U, V_t and S to reflect the decomposition of matrix. */
  virtual void UpdateFromMatrix() = 0;

  /** Compute the pseudo-inverse of matrix */
  gsl_matrix* PseudoInverse() const;

  /** Compute the pseudo-inverse of matrix using a damping factor of lambda. */
  gsl_matrix* PseudoInverse(double lambda) const;

  /** Print the SVD to standard out */
  void print() const;

};

#endif


text:
/*
 * agentlib.h
 *
 *  Created on: Jul. 5, 2019
 *      Author: takis
 */

#ifndef MY_INCLUDES_AGENTLIB_H_
#define MY_INCLUDES_AGENTLIB_H_

/********************************************************************************************************
 * 											INCLUDES
 * ******************************************************************************************************
 */
#include <vector>
#include <string>
#include <fstream>

// my own includes
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>
#include <math.h>

/********************************************************************************************************
 * 											MACROS
 * ******************************************************************************************************
 */
#define		Nx				100		// number of columns in physical environment
#define		Ny				100 	// number of rows in physical environment
#define		MaxMapEntries	48		// maximum number of contexts
#define		MaxIterations	9999	// maximum number of times main mapping loop can execute
#define		SenseRange		75		// units of distance that the agent can "see"
#define		MinVectorNorm	0.01	// a landmark must be further away than 1 cm from the agent
#define		MATSIZE			8		// matrix dimension (MATSIZE x MATSIZE) for analysis

/********************************************************************************************************
 * 											TYPEDEFs
 * ******************************************************************************************************
 */
typedef	std::vector<std::vector<std::string> >	vvs_t;	// 2D vector type for the
														// physical environment

typedef	vvs_t::size_type	vvs_sz_t;	// size_type for 2D vector type

typedef	enum 	e_region_t 	// enumerate the possible regions
	{
		R0,
		R1,
		R2,
		R3,
		R4,
		R5
	}	region_t;

typedef	enum	e_orient_t	// enumerate the possible directions
	{
		North,
		NorthEast,
		East,
		SouthEast,
		South,
		SouthWest,
		West,
		NorthWest
	}	orient_t;

typedef	enum	e_scan_direction_t	// scan goes from top left to bottom right,
									// or from bottom right to top left
	{
		SE,
		NW
	} 	scan_direction_t;

struct	loc_t	// type for a location in the physical environment
{
	vvs_sz_t		X;
	vvs_sz_t		Y;
};

struct	read_data_t		// type for storing data read from file
{
	std::string		context;
	region_t		region;
};

struct dprec_t // define a type to store dot products and corresponding directions
{
	orient_t dir;
	double dp;
};

/********************************************************************************************************
 * 											FUNCTION PROTOTYPES
 * ******************************************************************************************************
 */

/*
 * 		function: 	sgn
 * 		purpose:	determines the sign of a double variable
 * 		parameters:	a real number (double)
 * 		returns:	+1 or -1 (short int)
 *
 */
short int	sgn(double);

/*
 * 		function: 	load_known_contexts_from_file
 * 		purpose:	reads pairs of tokens from a file of the form <string, unsigned_integer> <newline>
 * 		parameters:	string consisting of the file name
 * 		returns:	a vector of read_data_t elements
 *
 */
std::vector<read_data_t> loadKnownContextsfromFile(const std::string);

/*
 * 		function: 	load_global_map_with_file_data
 * 		purpose:	fills contextual/associative map data structure with information in vector of read_data_t elements
 * 		parameters:	the associate map (passed by reference), and the vector of read_data_t elements, passed by constant reference
 * 		returns:	nothing (void function)space
 *
 */
void loadGlobalMapwithFileData(std::map<std::string, region_t>&, const std::vector<read_data_t>&);

/*
 * 		function: 	convertContextStringtoMatrix
 * 		purpose:	takes a key string associated with the spatial associative map, and converts it into a matrix
 * 		parameters:	the context string
 * 		returns:	a matrix (of type gsl_matrix_float*)
 */
gsl_matrix_float *convertContextStringtoMatrix(const std::string&);

/*
 * 		function:	mat2Norm
 * 		purpose:	compute the 2-norm of a matrix from the definition on page 291 of
 * 					Matrix Analysis, Horn & Johnson
 * 		parameters:	a matrix of type gsl_matrix_float*
 * 		returns:	the matrix 2-norm value
 */
float mat2Norm(const gsl_matrix_float);

/*
 * 		function: 	matTr
 * 		purpose: 	yields the transformation matrix we need, T_k, where k=0,+/-1,+/-2,...
 * 		parameters:	an integer (k)
 * 		returns:	a matrix of type gsl_matrix_float*
 *
 */
gsl_matrix_float *matTr(int);

/*
 * 		function:	compareContexts
 * 		purpose:	compare context matrices, return the closeness value
 * 		parameters:	two matrices to be compared, of type gsl_matrix_float*
 * 		returns:	a closeness value, the greater, the closer the matrices are
 */
float	compareContexts(const gsl_matrix_float*, const gsl_matrix_float*);

/*
 * 		function: 	alphaK
 * 		purpose:	implements a continuous monotonically decreasing function
 * 		parameters: a real number, x (float)
 * 		returns:	(float) exp(-x)
 */
float alphaK(const float);

/*
 * 		function:	agentScan360
 * 		purpose:	at a fixed position (loc_t, first parameter) in the physical environment, deduce the context string,
 * 					which consists of relative position data, e.g., lm_A to East, lm_B to SW, to be identified
 * 					as the key AEBSW for the associative map;
 * 		parameters: the agent's location, the associative map, the physical map, and the landmark locations
 * 		returns:	nothing (changes are made to the associative map)
 */
void agentScan360(loc_t, std::map<std::string, region_t>&, const vvs_t&, const std::vector<loc_t>&, const std::vector<std::string>&);

/*
 * 		function:	scan
 * 		purpose:	survey a rectangular subset of 2D space, either in a SE or a NW sense, using recursion
 * 		parameters:	rA (top left corner of scan), rB (bottom right corner of scan), r0 (agent current location)
 * 					spatial_mem (associative map of space), physmap (physical map of space)
 * 		returns: 	nothing (void type, but changes are made to the spatial memory, a.k.a., the associate map)
 */
void scan(bool& done, const scan_direction_t, 		loc_t, 		loc_t, 		loc_t, 		std::map<std::string, region_t>&, 				const vvs_t&, 			const std::map<char, loc_t>&, 			const std::vector<std::string>&);
//   scan(bool& done, const scan_direction_t dir, 	loc_t rA, 	loc_t rB, 	loc_t r0, 	std::map<std::string, region_t>& spatial_mem, 	const vvs_t& physmap, 	const std::map<char, loc_t>& landmarks, const std::vector<std::string>& orient_str)

#endif /* MY_INCLUDES_AGENTLIB_H_ */

max_stars_repo_name: takiszourntos/accelerated-cpp

max_issues_repo_name: takiszourntos/accelerated-cpp

max_forks_repo_name: takiszourntos/accelerated-cpp


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>
#include <time.h>
#include <math.h>
#include <stdlib.h>

#ifndef RANFUN_H
#define RANFUN_H

#define min(x,y) ((x) < (y) ? (x) : (y))
#define max(x,y) ((x) > (y) ? (x) : (y))

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

//interface for gsl_compare
// x==y, 0, x<y, -1, x>y, 1.
int fcmp (double x, double y, double epsilon = 0.00000000001);

/*Interface with gsl*/

/*Available random number generators from gsl
 gsl_rng_taus

       see http://sourceware.cygnus.com/gsl/ref/gsl-ref_5.html*/

#define GENERATOR gsl_rng_taus

/*to use inline functions in gsl (in c++)*/
/*define  DHAVE_INLINE at compile time*/
#ifdef DHAVE_INLINE
#define  INLINE inline
#endif

#ifndef  DHAVE_INLINE
#define  INLINE /*nothing*/
#endif

void Seed(unsigned long int s);

unsigned long int GetSeed();

double Un01();  /*Un01() */

double Unab(double a, double b);  /*U(a,b]*/

unsigned int Un0n(int n); /*Uniform integer in { 0, 1, ..., n-1 }*/

int Sign(); /*random -1 or 1*/

double NorSim(double m, double sigma); /*Normal*/

double ExpSim(double mu); /*p(x) dx = {1 \over \mu} \exp(-x/\mu) dx*/

double GammaSim(double a, double b); /*Gamma: p(x) dx = K x^{a-1} e^{-x/b} dx*/
double GammaDens(double x, double a, double b); /*density for beta, as above */

double BetaSim(double a, double b); /*Beta: p(x) dx = K x^{a-1} (1-x)^{b-1} dx */
double BetaDens(double x, double a, double b); /*density for beta, as above */

unsigned int BiSim( double p, int n); /*Binomial: p^k (1-p)^{n-k} */

/*Log of Special function Gamma*/

double LogGamma(double x);

#endif

max_stars_repo_name: SFlantua/Workshops

max_issues_repo_name: SFlantua/Workshops

max_forks_repo_name: SFlantua/Workshops


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

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

#include <gsl/gsl_ieee_utils.h>

#include "test_funcs.h"
int test_fdf (const char * desc, gsl_multiroot_function_fdf * function, initpt_function initpt, double factor, const gsl_multiroot_fdfsolver_type * T);
int test_f (const char * desc, gsl_multiroot_function_fdf * fdf, initpt_function initpt, double factor, const gsl_multiroot_fsolver_type * T);

int
main (void)
{
  const gsl_multiroot_fsolver_type * fsolvers[5] ;
  const gsl_multiroot_fsolver_type ** T1 ;

  const gsl_multiroot_fdfsolver_type * fdfsolvers[5] ;
  const gsl_multiroot_fdfsolver_type ** T2 ;

  double f;

  fsolvers[0] = gsl_multiroot_fsolver_dnewton;
  fsolvers[1] = gsl_multiroot_fsolver_broyden;
  fsolvers[2] = gsl_multiroot_fsolver_hybrid;
  fsolvers[3] = gsl_multiroot_fsolver_hybrids;
  fsolvers[4] = 0;

  fdfsolvers[0] = gsl_multiroot_fdfsolver_newton;
  fdfsolvers[1] = gsl_multiroot_fdfsolver_gnewton;
  fdfsolvers[2] = gsl_multiroot_fdfsolver_hybridj;
  fdfsolvers[3] = gsl_multiroot_fdfsolver_hybridsj;
  fdfsolvers[4] = 0;

  gsl_ieee_env_setup();

  f = 1.0 ;

  T1 = fsolvers ;

  while (*T1 != 0)
    {
      test_f ("Rosenbrock", &rosenbrock, rosenbrock_initpt, f, *T1);
      test_f ("Roth", &roth, roth_initpt, f, *T1);
      test_f ("Powell badly scaled", &powellscal, powellscal_initpt, f, *T1);
      test_f ("Brown badly scaled", &brownscal, brownscal_initpt, f, *T1);
      test_f ("Powell singular", &powellsing, powellsing_initpt, f, *T1);
      test_f ("Wood", &wood, wood_initpt, f, *T1);
      test_f ("Helical", &helical, helical_initpt, f, *T1);
      test_f ("Discrete BVP", &dbv, dbv_initpt, f, *T1);
      test_f ("Trig", &trig, trig_initpt, f, *T1);
      T1++;
    }

  T2 = fdfsolvers ;

  while (*T2 != 0)
    {
      test_fdf ("Rosenbrock", &rosenbrock, rosenbrock_initpt, f, *T2);
      test_fdf ("Roth", &roth, roth_initpt, f, *T2);
      test_fdf ("Powell badly scaled", &powellscal, powellscal_initpt, f, *T2);
      test_fdf ("Brown badly scaled", &brownscal, brownscal_initpt, f, *T2);
      test_fdf ("Powell singular", &powellsing, powellsing_initpt, f, *T2);
      test_fdf ("Wood", &wood, wood_initpt, f, *T2);
      test_fdf ("Helical", &helical, helical_initpt, f, *T2);
      test_fdf ("Discrete BVP", &dbv, dbv_initpt, f, *T2);
      test_fdf ("Trig", &trig, trig_initpt, f, *T2);
      T2++;
    }

  exit (gsl_test_summary ());
}

void scale (gsl_vector * x, double factor);

void
scale (gsl_vector * x, double factor)
{
  size_t i, n = x->size;

  if (gsl_vector_isnull(x))
    {
      for (i = 0; i < n; i++)
        {
          gsl_vector_set (x, i, factor);
        }
    }
  else
    {
      for (i = 0; i < n; i++)
        {
          double xi = gsl_vector_get(x, i);
          gsl_vector_set(x, i, factor * xi);
        }
    }
}

int
test_fdf (const char * desc, gsl_multiroot_function_fdf * function,
          initpt_function initpt, double factor,
          const gsl_multiroot_fdfsolver_type * T)
{
  int status;
  double residual = 0;
  size_t i, n = function->n, iter = 0;

  gsl_vector *x = gsl_vector_alloc (n);
  gsl_matrix *J = gsl_matrix_alloc (n, n);

  gsl_multiroot_fdfsolver *s;

  (*initpt) (x);

  if (factor != 1.0) scale(x, factor);

  s = gsl_multiroot_fdfsolver_alloc (T, n);
  gsl_multiroot_fdfsolver_set (s, function, x);

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

      if (status)
        break ;

      status = gsl_multiroot_test_residual (s->f, 0.0000001);
    }
  while (status == GSL_CONTINUE && iter < 1000);

#ifdef DEBUG
  printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n");
  printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n");
#endif

#ifdef TEST_JACOBIAN
 {
    double r,sum; size_t j;

    gsl_multiroot_function f1 ;
    f1.f = function->f ;
    f1.n = function->n ;
    f1.params = function->params ;

    gsl_multiroot_fdjacobian (&f1, s->x, s->f, GSL_SQRT_DBL_EPSILON, J);

    /* compare J and s->J */

    r=0;sum=0;
    for (i = 0; i < n; i++)
      for (j = 0; j< n ; j++)
        {
          double u = gsl_matrix_get(J,i,j);
          double su = gsl_matrix_get(s->J, i, j);
          r = fabs(u - su)/(1e-6 + 1e-6 * fabs(u)); sum+=r;
          if (fabs(u - su) > 1e-6 + 1e-6 * fabs(u))
            printf("broken jacobian %g\n", r);
        }
    printf("avg r = %g\n", sum/(n*n));
  }
#endif

  for (i = 0; i < n ; i++)
    {
      residual += fabs(gsl_vector_get(s->f, i));
    }

  gsl_multiroot_fdfsolver_free (s);
  gsl_matrix_free(J);
  gsl_vector_free(x);

  gsl_test(status, "%s on %s (%g), %u iterations, residual = %.2g", T->name, desc, factor, iter, residual);

  return status;
}

int
test_f (const char * desc, gsl_multiroot_function_fdf * fdf,
        initpt_function initpt, double factor,
        const gsl_multiroot_fsolver_type * T)
{
  int status;
  size_t i, n = fdf->n, iter = 0;
  double residual = 0;

  gsl_vector *x;

  gsl_multiroot_fsolver *s;
  gsl_multiroot_function function;

  function.f = fdf->f;
  function.params = fdf->params;
  function.n = n ;

  x = gsl_vector_alloc (n);

  (*initpt) (x);

  if (factor != 1.0) scale(x, factor);

  s = gsl_multiroot_fsolver_alloc (T, n);
  gsl_multiroot_fsolver_set (s, &function, x);

/*   printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n"); */
/*   printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n"); */

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

      if (status)
        break ;

      status = gsl_multiroot_test_residual (s->f, 0.0000001);
    }
  while (status == GSL_CONTINUE && iter < 1000);

#ifdef DEBUG
  printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n");
  printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n");
#endif

  for (i = 0; i < n ; i++)
    {
      residual += fabs(gsl_vector_get(s->f, i));
    }

  gsl_multiroot_fsolver_free (s);
  gsl_vector_free(x);

  gsl_test(status, "%s on %s (%g), %u iterations, residual = %.2g", T->name, desc, factor, iter, residual);

  return status;
}

ext: c


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <sys/stat.h>
#include <sys/types.h>

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

#ifdef SUBFIND

#include "fof.h"
#include "subfind.h"

/*! Structure for communication during the density computation. Holds data that is sent to other processors.
 */
static struct densdata_in
{
  MyDouble Pos[3];
  MyFloat Hsml;
  int NodeList[NODELISTLENGTH];
}
 *DensDataIn, *DensDataGet;

static struct densdata_out
{
  MyFloat Rho;
  MyFloat VelDisp, Vx, Vy, Vz;
  int Ngb;
}
 *DensDataResult, *DensDataOut;

static MyFloat *DM_Vx, *DM_Vy, *DM_Vz;
static long long Ntotal;

#ifdef DENSITY_SPLIT_BY_TYPE
void subfind_density(int j_in)
#else
void subfind_density(void)
#endif
{
  long long ntot;
  int i, j, ndone, ndone_flag, npleft, dummy, iter = 0;
  MyFloat *Left, *Right;
  char *Todo;
  int ngrp, sendTask, recvTask, place, nexport, nimport;
  double vel_to_phys, dmax1, dmax2, t0, t1;

#ifdef DENSITY_SPLIT_BY_TYPE
  int j_target;

  if(j_in < 0)
    j_target = -j_in - 1;
  else
    j_target = j_in;
#endif

  if(ThisTask == 0)
    {
#ifdef DENSITY_SPLIT_BY_TYPE
      printf("finding densities, mode %d \n", j_in);
#else
      printf("finding densities for all particles\n");
#endif
      fflush(stdout);
    }

  /* allocate buffers to arrange communication */

  Ngblist = (int *) mymalloc(NumPart * sizeof(int));
  Dist2list = (double *) mymalloc(NumPart * sizeof(double));

  All.BunchSize =
    (int) ((All.BufferSize * 1024 * 1024) / (sizeof(struct data_index) + sizeof(struct data_nodelist) +
					     sizeof(struct densdata_in) + sizeof(struct densdata_out) +
					     sizemax(sizeof(struct densdata_in),
						     sizeof(struct densdata_out))));
  DataIndexTable = (struct data_index *) mymalloc(All.BunchSize * sizeof(struct data_index));
  DataNodeList = (struct data_nodelist *) mymalloc(All.BunchSize * sizeof(struct data_nodelist));

  Left = mymalloc(sizeof(MyFloat) * NumPart);
  Right = mymalloc(sizeof(MyFloat) * NumPart);
  Todo = mymalloc(sizeof(char) * NumPart);

  DM_Vx = mymalloc(sizeof(MyFloat) * NumPart);
  DM_Vy = mymalloc(sizeof(MyFloat) * NumPart);
  DM_Vz = mymalloc(sizeof(MyFloat) * NumPart);

  for(i = 0; i < NumPart; i++)
    {
      Left[i] = Right[i] = 0;
      P[i].DM_NumNgb = 0;
      Todo[i] = 1;
    }

  /* we will repeat the whole thing for those particles where we didn't find enough neighbours */
  do
    {
      t0 = second();

      i = 0;			/* begin with this index */

      do
	{
	  for(j = 0; j < NTask; j++)
	    {
	      Send_count[j] = 0;
	      Exportflag[j] = -1;
	    }

	  /* do local particles and prepare export list */

	  for(nexport = 0; i < NumPart; i++)
	    {
	      if(Todo[i])
#ifdef DENSITY_SPLIT_BY_TYPE
		if(P[i].Type == j_target)
		  {
		    if(subfind_density_evaluate(i, 0, &nexport, Send_count, j_in) < 0)
		      break;
		  }
#else
		if(((1 << P[i].Type) & (FOF_PRIMARY_LINK_TYPES)))
		  {
		    if(subfind_density_evaluate(i, 0, &nexport, Send_count) < 0)
		      break;
		  }
#endif
	    }

	  qsort(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);

	  MPI_Allgather(Send_count, NTask, MPI_INT, Sendcount_matrix, NTask, MPI_INT, MPI_COMM_WORLD);

	  for(j = 0, nimport = 0, Recv_offset[0] = 0, Send_offset[0] = 0; j < NTask; j++)
	    {
	      Recv_count[j] = Sendcount_matrix[j * NTask + ThisTask];
	      nimport += Recv_count[j];

	      if(j > 0)
		{
		  Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
		  Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
		}
	    }

	  DensDataGet = (struct densdata_in *) mymalloc(nimport * sizeof(struct densdata_in));
	  DensDataIn = (struct densdata_in *) mymalloc(nexport * sizeof(struct densdata_in));

	  /* prepare particle data for export */
	  for(j = 0; j < nexport; j++)
	    {
	      place = DataIndexTable[j].Index;

	      DensDataIn[j].Pos[0] = P[place].Pos[0];
	      DensDataIn[j].Pos[1] = P[place].Pos[1];
	      DensDataIn[j].Pos[2] = P[place].Pos[2];
	      DensDataIn[j].Hsml = P[place].DM_Hsml;

	      memcpy(DensDataIn[j].NodeList,
		     DataNodeList[DataIndexTable[j].IndexGet].NodeList, NODELISTLENGTH * sizeof(int));
	    }

	  /* exchange particle data */
	  for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	    {
	      sendTask = ThisTask;
	      recvTask = ThisTask ^ ngrp;

	      if(recvTask < NTask)
		{
		  if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		    {
		      /* get the particles */
		      MPI_Sendrecv(&DensDataIn[Send_offset[recvTask]],
				   Send_count[recvTask] * sizeof(struct densdata_in), MPI_BYTE,
				   recvTask, TAG_DENS_A,
				   &DensDataGet[Recv_offset[recvTask]],
				   Recv_count[recvTask] * sizeof(struct densdata_in), MPI_BYTE,
				   recvTask, TAG_DENS_A, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		    }
		}
	    }

	  myfree(DensDataIn);
	  DensDataResult = (struct densdata_out *) mymalloc(nimport * sizeof(struct densdata_out));
	  DensDataOut = (struct densdata_out *) mymalloc(nexport * sizeof(struct densdata_out));

	  /* now do the particles that were sent to us */
	  for(j = 0; j < nimport; j++)
#ifdef DENSITY_SPLIT_BY_TYPE
	    subfind_density_evaluate(j, 1, &dummy, &dummy, j_in);
#else
	    subfind_density_evaluate(j, 1, &dummy, &dummy);
#endif
	  if(i >= NumPart)
	    ndone_flag = 1;
	  else
	    ndone_flag = 0;

	  MPI_Allreduce(&ndone_flag, &ndone, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

	  /* get the result */
	  for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	    {
	      sendTask = ThisTask;
	      recvTask = ThisTask ^ ngrp;
	      if(recvTask < NTask)
		{
		  if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		    {
		      /* send the results */
		      MPI_Sendrecv(&DensDataResult[Recv_offset[recvTask]],
				   Recv_count[recvTask] * sizeof(struct densdata_out),
				   MPI_BYTE, recvTask, TAG_DENS_B,
				   &DensDataOut[Send_offset[recvTask]],
				   Send_count[recvTask] * sizeof(struct densdata_out),
				   MPI_BYTE, recvTask, TAG_DENS_B, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		    }
		}
	    }

	  /* add the result to the local particles */
	  for(j = 0; j < nexport; j++)
	    {
	      place = DataIndexTable[j].Index;

	      P[place].DM_NumNgb += DensDataOut[j].Ngb;
	      P[place].u.DM_Density += DensDataOut[j].Rho;
#ifdef DENSITY_SPLIT_BY_TYPE
	      if(j_in >= 0)	/* sum up VelDisp only for own species */
#endif
		P[place].v.DM_VelDisp += DensDataOut[j].VelDisp;
	      DM_Vx[place] += DensDataOut[j].Vx;
	      DM_Vy[place] += DensDataOut[j].Vy;
	      DM_Vz[place] += DensDataOut[j].Vz;
	    }

	  myfree(DensDataOut);
	  myfree(DensDataResult);
	  myfree(DensDataGet);
	}
      while(ndone < NTask);

#ifdef DENSITY_SPLIT_BY_TYPE
      if(ThisTask == 0)
	printf("checking smoothing lenghs ... \n");
#endif

      /* do final operations on results */
      for(i = 0, npleft = 0; i < NumPart; i++)
	{
	  /* now check whether we had enough neighbours */
#ifdef DENSITY_SPLIT_BY_TYPE
	  if(j_in < 0)
	    Todo[i] = 0;
	  if((P[i].Type != j_target) || (j_in < 0))
	    continue;
#else
	  if(!((1 << P[i].Type) & (FOF_PRIMARY_LINK_TYPES)))
	    continue;
#endif

	  if(Todo[i])
	    {
	      if(P[i].DM_NumNgb != All.DesNumNgb &&
		 ((Right[i] - Left[i]) > 1.0e-4 * Left[i] || Left[i] == 0 || Right[i] == 0))
		{
		  /* need to redo this particle */
		  npleft++;

		  if(P[i].DM_NumNgb < All.DesNumNgb)
		    Left[i] = DMAX(P[i].DM_Hsml, Left[i]);
		  else
		    {
		      if(Right[i] != 0)
			{
			  if(P[i].DM_Hsml < Right[i])
			    Right[i] = P[i].DM_Hsml;
			}
		      else
			Right[i] = P[i].DM_Hsml;
		    }

		  if(iter >= MAXITER - 10)
		    {
		      printf
			("i=%d task=%d ID=%d Hsml=%g Left=%g Right=%g Ngbs=%g Right-Left=%g\n   pos=(%g|%g|%g)\n",
			 i, ThisTask, (int) P[i].ID, P[i].DM_Hsml, Left[i], Right[i],
			 (double) P[i].DM_NumNgb, Right[i] - Left[i], P[i].Pos[0], P[i].Pos[1], P[i].Pos[2]);
		      fflush(stdout);
		    }

		  if(Right[i] > 0 && Left[i] > 0)
		    P[i].DM_Hsml = pow(0.5 * (pow(Left[i], 3) + pow(Right[i], 3)), 1.0 / 3);
		  else
		    {
		      if(Right[i] == 0 && Left[i] == 0)
			endrun(8187);	/* can't occur */

		      if(Right[i] == 0 && Left[i] > 0)
			P[i].DM_Hsml *= 1.26;

		      if(Right[i] > 0 && Left[i] == 0)
			P[i].DM_Hsml /= 1.26;
		    }
		}
	      else
		Todo[i] = 0;
	    }
	}

      sumup_large_ints(1, &npleft, &ntot);

      t1 = second();

      if(ntot > 0)
	{
	  iter++;

	  if(iter > 0 && ThisTask == 0)
	    {
	      printf("ngb iteration %d: need to repeat for %d%09d particles. (took %g sec)\n", iter,
		     (int) (ntot / 1000000000), (int) (ntot % 1000000000), timediff(t0, t1));
	      fflush(stdout);
	    }

	  if(iter > MAXITER)
	    {
	      printf("failed to converge in neighbour iteration in density()\n");
	      fflush(stdout);
	      endrun(1155);
	    }
	}
    }
  while(ntot > 0);

  vel_to_phys = 1.0 / All.Time;

#ifdef DENSITY_SPLIT_BY_TYPE
  if(ThisTask == 0)
    printf("final operations ... \n");
#endif

  for(i = 0; i < NumPart; i++)
#ifdef DENSITY_SPLIT_BY_TYPE
    if(P[i].Type == j_target && j_in >= 0)
#else
    if(((1 << P[i].Type) & (FOF_PRIMARY_LINK_TYPES)))
#endif
      {
#ifdef DENSITY_SPLIT_BY_TYPE
	if(P[i].DM_NumNgb == 0)
	  endrun(5621);
#endif
	DM_Vx[i] /= P[i].DM_NumNgb;
	DM_Vy[i] /= P[i].DM_NumNgb;
	DM_Vz[i] /= P[i].DM_NumNgb;
	P[i].v.DM_VelDisp /= P[i].DM_NumNgb;

	P[i].v.DM_VelDisp = vel_to_phys * sqrt(P[i].v.DM_VelDisp -
					       DM_Vx[i] * DM_Vx[i] -
					       DM_Vy[i] * DM_Vy[i] - DM_Vz[i] * DM_Vz[i]);
      }

  myfree(DM_Vz);
  myfree(DM_Vy);
  myfree(DM_Vx);

  myfree(Todo);
  myfree(Right);
  myfree(Left);

  myfree(DataNodeList);
  myfree(DataIndexTable);

  myfree(Dist2list);
  myfree(Ngblist);
}

/*! This function represents the core of the SPH density computation. The
 *  target particle may either be local, or reside in the communication
 *  buffer.
 */
#ifdef DENSITY_SPLIT_BY_TYPE
int subfind_density_evaluate(int target, int mode, int *nexport, int *nsend_local, int tp)
#else
int subfind_density_evaluate(int target, int mode, int *nexport, int *nsend_local)
#endif
{
  int j, n;
  int startnode, numngb, ngb, listindex = 0;
  double hmax;
  double h, h2, hinv, hinv3;
  double rho, wk;
  double r, r2, u, mass_j, v2, vx, vy, vz;
  MyDouble *pos;

  rho = 0;
  numngb = 0;
  v2 = vx = vy = vz = 0;

  if(mode == 0)
    {
      pos = P[target].Pos;
      h = P[target].DM_Hsml;
    }
  else
    {
      pos = DensDataGet[target].Pos;
      h = DensDataGet[target].Hsml;
    }

#ifdef DENSITY_SPLIT_BY_TYPE
  if(h == 0)
    endrun(5622);
#endif

  h2 = h * h;
  hinv = 1.0 / h;
  hinv3 = hinv * hinv * hinv;

  if(mode == 0)
    {
      startnode = All.MaxPart;	/* root node */
    }
  else
    {
      startnode = DensDataGet[target].NodeList[0];
      startnode = Nodes[startnode].u.d.nextnode;	/* open it */
    }

  numngb = 0;

  while(startnode >= 0)
    {
      while(startnode >= 0)
	{
#ifdef DENSITY_SPLIT_BY_TYPE
	  if(tp < 0)
	    ngb =
	      subfind_ngb_treefind_linkpairs(pos, h, target, &startnode, mode, &hmax, nexport, nsend_local);
	  else
	    ngb = subfind_ngb_treefind_linkngb(pos, h, target, &startnode, mode, &hmax, nexport, nsend_local);
#else
	  ngb = subfind_ngb_treefind_linkngb(pos, h, target, &startnode, mode, &hmax, nexport, nsend_local);
#endif
	  if(ngb < 0)
	    return -1;

#ifdef DENSITY_SPLIT_BY_TYPE
	  if(tp >= 0)
#endif
	    if(mode == 0 && hmax > 0)
	      {
		P[target].DM_Hsml = hmax;
		h = hmax;
		h2 = h * h;
		hinv = 1.0 / h;
		hinv3 = hinv * hinv * hinv;

		if(ngb != All.DesNumNgb)
		  endrun(121);
	      }

	  numngb += ngb;

	  for(n = 0; n < ngb; n++)
	    {
	      j = Ngblist[n];

	      r2 = Dist2list[n];

#ifdef DENSITY_SPLIT_BY_TYPE
	      if(tp < 0)
		{
		  h = P[j].DM_Hsml;
		  if(h == 0)
		    endrun(5622);
		  h2 = h * h;
		  hinv = 1.0 / h;
		  hinv3 = hinv * hinv * hinv;
		}
#endif

	      if(r2 < h2)
		{
		  r = sqrt(r2);

		  u = r * hinv;

		  if(u < 0.5)
		    wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
		  else
		    wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);

		  mass_j = P[j].Mass;

		  rho += (mass_j * wk);
		}

	      vx += P[j].Vel[0];
	      vy += P[j].Vel[1];
	      vz += P[j].Vel[2];

	      v2 += P[j].Vel[0] * P[j].Vel[0] + P[j].Vel[1] * P[j].Vel[1] + P[j].Vel[2] * P[j].Vel[2];
	    }
	}

      if(mode == 1)
	{
	  listindex++;
	  if(listindex < NODELISTLENGTH)
	    {
	      startnode = DensDataGet[target].NodeList[listindex];
	      if(startnode >= 0)
		startnode = Nodes[startnode].u.d.nextnode;	/* open it */
	    }
	}
    }

  if(mode == 0)
    {
      P[target].DM_NumNgb = numngb;
      P[target].u.DM_Density = rho;
#ifdef DENSITY_SPLIT_BY_TYPE
      if(tp >= 0)		/* sum up VelDisp only for own species */
#endif
	P[target].v.DM_VelDisp = v2;
      DM_Vx[target] = vx;
      DM_Vy[target] = vy;
      DM_Vz[target] = vz;
    }
  else
    {
      DensDataResult[target].Ngb = numngb;
      DensDataResult[target].Rho = rho;
      DensDataResult[target].VelDisp = v2;
      DensDataResult[target].Vx = vx;
      DensDataResult[target].Vy = vy;
      DensDataResult[target].Vz = vz;
    }

  return 0;
}

#ifdef DENSITY_SPLIT_BY_TYPE
void subfind_setup_smoothinglengths(int j)
#else
void subfind_setup_smoothinglengths(void)
#endif
{
  int i, no, p;

  for(i = 0; i < NumPart; i++)
    {
#ifdef DENSITY_SPLIT_BY_TYPE
      if(P[i].Type == j)
#else
      if(((1 << P[i].Type) & (FOF_PRIMARY_LINK_TYPES)))
#endif
	{
	  no = Father[i];

	  /* Not a good guess for gas/stars component, need more thought ! */
	  while(10 * All.DesNumNgb * P[i].Mass > Nodes[no].u.d.mass)
	    {
	      p = Nodes[no].u.d.father;

	      if(p < 0)
		break;

	      no = p;
	    }

#ifdef DENSITY_SPLIT_BY_TYPE
	  if(P[i].Type == 0)
	    P[i].DM_Hsml = PPP[i].Hsml;
	  else
#ifdef KEEP_DM_HSML_AS_GUESS
	  if(P[i].DM_Hsml < 0)
#endif
	    P[i].DM_Hsml =
	      pow(3.0 / (4 * M_PI) * All.DesNumNgb * P[i].Mass / Nodes[no].u.d.mass, 1.0 / 3) * Nodes[no].len;
#else
	  P[i].DM_Hsml =
	    pow(3.0 / (4 * M_PI) * All.DesNumNgb * P[i].Mass / Nodes[no].u.d.mass, 1.0 / 3) * Nodes[no].len;
#endif
	}
    }
}

static int Nhsml;

static struct hsml_data
{
  float Hsml;
  float Density;
  float VelDisp;
  MyIDType ID;
}
 *Hsml_list;

int subfind_compare_hsml_data(const void *a, const void *b)
{
  if(((struct hsml_data *) a)->ID < ((struct hsml_data *) b)->ID)
    return -1;

  if(((struct hsml_data *) a)->ID > ((struct hsml_data *) b)->ID)
    return +1;

  return 0;
}

void subfind_save_densities(int num)
{
  int i, nprocgroup, masterTask, groupTask;
  char buf[1000];
  double t0, t1;

  if(ThisTask == 0)
    {
      printf("start saving smoothing lengths and densities\n");
      fflush(stdout);
    }

  for(i = 0, Nhsml = 0; i < NumPart; i++)
#ifdef DENSITY_SPLIT_BY_TYPE
    if(((1 << P[i].Type) & (DENSITY_SPLIT_BY_TYPE)))
#else
    if(((1 << P[i].Type) & (FOF_PRIMARY_LINK_TYPES)))
#endif
      Nhsml++;

  MPI_Allgather(&Nhsml, 1, MPI_INT, Send_count, 1, MPI_INT, MPI_COMM_WORLD);
  for(i = 1, Send_offset[0] = 0; i < NTask; i++)
    Send_offset[i] = Send_offset[i - 1] + Send_count[i - 1];

  sumup_large_ints(1, &Nhsml, &Ntotal);

  Hsml_list = mymalloc(Nhsml * sizeof(struct hsml_data));

  for(i = 0, Nhsml = 0; i < NumPart; i++)
#ifdef DENSITY_SPLIT_BY_TYPE
    if(((1 << P[i].Type) & (DENSITY_SPLIT_BY_TYPE)))
#else
    if(((1 << P[i].Type) & (FOF_PRIMARY_LINK_TYPES)))
#endif
      {
	Hsml_list[Nhsml].Hsml = P[i].DM_Hsml;
	Hsml_list[Nhsml].Density = P[i].u.DM_Density;
	Hsml_list[Nhsml].VelDisp = P[i].v.DM_VelDisp;
#ifdef SAVE_HSML_IN_IC_ORDER
	Hsml_list[Nhsml].ID = P[i].ID_ic_order;
#else
	Hsml_list[Nhsml].ID = P[i].ID;
#endif
	Nhsml++;
      }

  t0 = second();
  parallel_sort(Hsml_list, Nhsml, sizeof(struct hsml_data), subfind_compare_hsml_data);
  t1 = second();

  if(ThisTask == 0)
    {
      printf("Sorting of densities in ID sequence took = %g sec\n", timediff(t0, t1));
      fflush(stdout);
    }

  if(ThisTask == 0)
    {
      sprintf(buf, "%s/hsmldir_%03d", All.OutputDir, num);
      mkdir(buf, 02755);
    }
  MPI_Barrier(MPI_COMM_WORLD);

  if(NTask < All.NumFilesWrittenInParallel)
    {
      printf
	("Fatal error.\nNumber of processors must be a smaller or equal than `NumFilesWrittenInParallel'.\n");
      endrun(241931);
    }

  nprocgroup = NTask / All.NumFilesWrittenInParallel;
  if((NTask % All.NumFilesWrittenInParallel))
    nprocgroup++;
  masterTask = (ThisTask / nprocgroup) * nprocgroup;
  for(groupTask = 0; groupTask < nprocgroup; groupTask++)
    {
      if(ThisTask == (masterTask + groupTask))	/* ok, it's this processor's turn */
	subfind_save_local_densities(num);
      MPI_Barrier(MPI_COMM_WORLD);	/* wait inside the group */
    }

  myfree(Hsml_list);

}

void subfind_save_local_densities(int num)
{
  char fname[1000];
  int i;
  float *tmp;
  FILE *fd;

  sprintf(fname, "%s/hsmldir_%03d/%s_%03d.%d", All.OutputDir, num, "hsml", num, ThisTask);
  if(!(fd = fopen(fname, "w")))
    {
      printf("can't open file `%s`\n", fname);
      endrun(1183);
    }

  my_fwrite(&Nhsml, sizeof(int), 1, fd);
  my_fwrite(&Send_offset[ThisTask], sizeof(int), 1, fd);	/* this is the number of IDs in previous files */
  my_fwrite(&Ntotal, sizeof(long long), 1, fd);
  my_fwrite(&NTask, sizeof(int), 1, fd);

  tmp = mymalloc(Nhsml * sizeof(float));

  for(i = 0; i < Nhsml; i++)
    tmp[i] = Hsml_list[i].Hsml;
  my_fwrite(tmp, sizeof(float), Nhsml, fd);

  for(i = 0; i < Nhsml; i++)
    tmp[i] = Hsml_list[i].Density;
  my_fwrite(tmp, sizeof(float), Nhsml, fd);

  for(i = 0; i < Nhsml; i++)
    tmp[i] = Hsml_list[i].VelDisp;
  my_fwrite(tmp, sizeof(float), Nhsml, fd);

  myfree(tmp);

  fclose(fd);
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <Python.h>
#include <numpy/arrayobject.h>

#include <math.h>

#ifdef __INTEL_COMPILER
#include <mkl_cblas.h>
#else
#include <cblas.h>
#define kmp_set_blocktime(k)
#endif

static inline float minf (float a, float b) { return a < b ? a : b; }
static inline float maxf (float a, float b) { return a < b ? b : a; }

PyObject * compute_overlap (float bbx1, float bby1, float bbx2, float bby2,
                            int fdimy, int fdimx, int dimy, int dimx,
                            float scale, int pady, int padx, int h, int w) {
    int im_area = h*w;
    int bbox_area = (bbx2-bbx1)*(bby2-bby1);
    int im_clip = ((double)bbox_area / (double)im_area) < 0.7;
    npy_intp dims[2] = { dimy, dimx };
    int x, y;

    PyArrayObject * pyoverlap = (PyArrayObject*)PyArray_SimpleNew ((npy_intp)2, dims, NPY_FLOAT);

    for (x = 0; x < dimx; ++x) {
        for (y = 0; y < dimy; ++y) {
            float x1 = (x - padx) * scale;
            float y1 = (y - pady) * scale;
            float x2 = x1 + fdimx*scale - 1;
            float y2 = y1 + fdimy*scale - 1;

            if (im_clip) {
                x1 = minf(maxf(x1, 0.0f), w-1);
                y1 = minf(maxf(y1, 0.0f), h-1);
                x2 = minf(maxf(x2, 0.0f), w-1);
                y2 = minf(maxf(y2, 0.0f), h-1);
            }

            float xx1 = maxf(x1, bbx1);
            float yy1 = maxf(y1, bby1);
            float xx2 = minf(x2, bbx2);
            float yy2 = minf(y2, bby2);

            float intw = xx2 - xx1 + 1;
            float inth = yy2 - yy1 + 1;

            if (intw > 0 && inth > 0) {
                float filterw = x2 - x1 + 1;
                float filterh = y2 - y1 + 1;
                float filter_area = filterw*filterh;
                float int_area = intw*inth;
                float union_area = filter_area + bbox_area - int_area;

                *(float*)PyArray_GETPTR2(pyoverlap, y, x) = int_area / union_area;
            } else {
                *(float*)PyArray_GETPTR2(pyoverlap, y, x) = 0;
            }
        }
    }

    return Py_BuildValue("N", pyoverlap);
}

#if 0
PyObject * objective_function (PyObject * pyexamples) {
    int i;
    int j;
    int k;

    for (i = 0; i < PyList_Size(pyexamples); ++i) {
        PyObject * pyexample = PyList_GetItem (pyexamples, i);
        if (!PyList_Check(pyexample)) {
            PyErr_SetString(PyExc_TypeError, "example list elements must be list.");
            return NULL;
        }

        for (j = 0; j < PyList_Size(pyexample); ++j) {

            PyObject * pyentry = PyList_GetItem(pyexample, j);
            if (!PyTuple_Check(pyentry)) {
                PyErr_SetString(PyExc_TypeError, "entries must be tuples.");
                return NULL;
            }

            if (PyTuple_Size(pyentry) != 2) {
                PyErr_SetString(PyExc_TypeError, "entry tuples must be of size 2");
                return NULL;
            }

            PyObject * pyfeatures = PyTuple_GetItem(pyentry, 0);
            if (!PyDict_Check(pyfeatures)) {
                PyErr_SetString(PyExc_TypeError, "features must be a dict");
                return NULL;
            }

            PyObject * pyfeatureslist = PyDict_Items(pyfeatures);

            float score = 0;
            for (k = 0; k < PyDict_Size(pyfeatureslist); ++k) {
                PyObject * pyfeature = PyList_GetItem(pyfeatureslist, k);
                PyObject * pyblock = PyTuple_GetItem(pyfeature, 0);

                PyObject * pyw = PyTuple_GetItem(pyfeature, 1);
                if (!PyArray_Check(pyw)) {
                    PyErr_SetString(PyExc_TypeError, "w must be a numpy array");
                    return NULL;
                }

                PyArrayObject * pywarray = (PyArrayObject*)pyw;

            }

            PyObject * pydetection = PyTuple_GetItem(pyentry, 1);
        }
    }
    return Py_BuildValue("i", 0);
}

PyObject * gradient (PyObject * examples, PyObject * gradients) {
    Py_RETURN_NONE;
}
#endif

static PyObject * ComputeOverlap(PyObject * self, PyObject * args)
{
    float x1, x2, y1, y2;
    int fdimy, fdimx, dimy, dimx;
    float scale;
    int pady, padx;
    int w, h;
    if (!PyArg_ParseTuple(args, "ffffiiiifiiii", &x1, &y1, &x2, &y2, &fdimy, &fdimx, &dimy, &dimx, &scale, &pady, &padx, &h, &w))
        return NULL;
    return compute_overlap(x1, y1, x2, y2, fdimy, fdimx, dimy, dimx, scale, pady, padx, h, w);
}

#if 0
static PyObject * ObjectiveFunction (PyObject * self, PyObject * args) {
    PyObject * pyexamples;
    if (!PyArg_ParseTuple (args, "O!", &PyList_Type, &pyexamples)) {
        return NULL;
    }

    return objective_function (pyexamples);
}

static PyObject * Gradient (PyObject * self, PyObject * args) {
    PyObject * pyexamples;
    PyObject * pygradients;
    if (!PyArg_ParseTuple (args, "O!O!", &PyList_Type, &pyexamples, &PyDict_Type, &pygradients)) {
        return NULL;
    }

    return gradient (pyexamples, pygradients);
}
#endif

#if PY_MAJOR_VERSION >= 3
static struct PyModuleDef moduledef = {
    PyModuleDef_HEAD_INIT,
    "_train",
    "Routines for training the DPM.",
    -1,
    NULL,
    NULL,
    NULL,
    NULL,
    NULL
};
#endif

#if PY_MAJOR_VERSION < 3
static PyMethodDef _train_methods[] = {
    {"compute_overlap", ComputeOverlap, METH_VARARGS, "Compute detection overlaps with bbox."},
#if 0
    {"ObjectiveFunction", ObjectiveFunction, METH_VARARGS, "WL-SSVM objective function."},
    {"Gradient", Gradient, METH_VARARGS, "WL-SSVM gradient."},
#endif
    {NULL}
};
#endif

#if PY_MAJOR_VERSION >= 3
PyMODINIT_FUNC
PyInit__train(void)
#else
PyMODINIT_FUNC
init_train(void)
#endif
{
    import_array();

#if PY_MAJOR_VERSION >= 3
    PyObject *m = PyModule_Create(&moduledef);
#else
    Py_InitModule3("_train", _train_methods, "Compute detection overlaps with bbox.");
#endif

#if PY_MAJOR_VERSION >= 3
    return m;
#endif
}

max_stars_repo_name: caomw/pydro

max_issues_repo_name: caomw/pydro

max_forks_repo_name: kmatzen/pydro


text:
#include "stdio.h"
#include "gsl/gsl_integration.h"
#include "gsl/gsl_errno.h"
//include "limber.h"
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_sf_bessel.h>

// The integrated C_ell are in general allowed to be zero or negative if
// they describe cross-correlations. We use these statuses to describe
// errors where a log was requested also.
// LIMBER_STATUS_NEGATIVE is probably always an error
// but LIMBER_STATUS_ZERO is not necessarily
#define LIMBER_STATUS_OK 0
#define LIMBER_STATUS_ZERO 1
#define LIMBER_STATUS_NEGATIVE 2
#define LIMBER_STATUS_ERROR 3

// These are the options you can set for
// the Limber integrator.
typedef struct cl_config{
	int n_ell;  // Number of ell values you want in the spline
	int * ell;  // The chosen ell values you want
	double prefactor; //Scaling prefactor
    int status; // did everything go okay?
    double delta_range_factor;
    double log_nu_range;
    double absolute_tolerance;
    double relative_tolerance;
} cl_config;

int cl_integral(cl_config * config, gsl_spline * WX_red,
			    gsl_spline * WY_red, gsl_spline2d * P, gsl_spline * D_chi,
			    double * cl_out, double * cl_err_out);


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include "stats.h"

double ttest_equalvar(const double *data1, int n1, const double *data2, int n2)
{
    double mean1, mean2, var1, var2, num, sp, denom, t, dof, p;
    mean1 = gsl_stats_mean(data1, 1, n1);
    mean2 = gsl_stats_mean(data2, 1, n2);
    var1 = gsl_stats_variance(data1, 1, n1);
    var2 = gsl_stats_variance(data2, 1, n2);
    num = mean1 - mean2;
    sp = sqrt((((n1 - 1) * var1) + ((n2 - 1) * var2)) / (n1 + n2 - 2));
    denom = sp * sqrt((1 / (double) n1) + (1 / (double) n2));
    t = fabs(num / denom);
    dof = (double) n1 + n2 - 2;
    p = (1 - gsl_cdf_tdist_P(t, dof)) * 2;
    return p;
}

double ttest_unequalvar(const double *data1, int n1, const double *data2, int n2)
{
    double mean1, mean2, var1, var2, r1, r2, t, dof, p;
    mean1 = gsl_stats_mean(data1, 1, n1);
    mean2 = gsl_stats_mean(data2, 1, n2);
    var1 = gsl_stats_variance(data1, 1, n1);
    var2 = gsl_stats_variance(data2, 1, n2);
    r1 = var1 / n1;
    r2 = var2 / n2;
    // Compute t-statistic
    t = fabs((mean1 - mean2) / sqrt(r1 + r2));
    dof = pow(r1 + r2, 2) / ((pow(r1, 2) / (n1 - 1)) + (pow(r2, 2) / (n2 - 1)));
    p = (1 - gsl_cdf_tdist_P(t, dof)) * 2;
    return p;
}

double ttest_paired(const double *data1, const double *data2, int n)
{
    double data[n];
    double mean, sd, t, dof, p;
    for (int i = 0; i < n; i++)
    {
        data[i] = data1[i] - data2[i];
    }
    p = ttest(data, n);
    return p;
}

double ttest(const double *data, int n)
{
    double mean, sd, t, dof, p;
    mean = gsl_stats_mean(data, 1, n);
    sd = gsl_stats_sd(data, 1, n);
    t = fabs((mean - 0) / (sd / sqrt((double) n)));
    dof = (double) n - 1;
    p = (1 - gsl_cdf_tdist_P(t, dof)) * 2;
    return p;
}

max_stars_repo_name: kgoettler/python-c-ext

max_issues_repo_name: kgoettler/python-c-ext

max_forks_repo_name: kgoettler/python-c-ext


text:
/* randist/wishart.c
 *
 * Copyright (C) 2017 Timothée Flutre
 *
 * 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_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_sf_gamma.h>

/* Generate a random matrix from a Wishart distribution using the Bartlett
 * decomposition, following Smith and Hocking, Journal of the Royal Statistical
 * Society. Series C (Applied Statistics), Vol. 21, No. 3 (1972), pp. 341-345.
 *
 * df      degrees of freedom
 * L       matrix resulting from the Cholesky decomposition of
 *         the scale matrix V = L L^T (dimension d x d)
 * result  output matrix (dimension d x d)
 * work    matrix used for intermediate computations (dimension d x d)
 */
int
gsl_ran_wishart (const gsl_rng * r,
                 const double df,
                 const gsl_matrix * L,
                 gsl_matrix * result,
                 gsl_matrix * work)
{
  if (L->size1 != L->size2)
    {
      GSL_ERROR("L should be a square matrix", GSL_ENOTSQR);
    }
  else if (result->size1 != result->size2)
    {
      GSL_ERROR("result should be a square matrix", GSL_ENOTSQR);
    }
  else if (work->size1 != work->size2)
    {
      GSL_ERROR("work should be a square matrix", GSL_ENOTSQR);
    }
  else if (result->size1 != L->size1)
    {
      GSL_ERROR("incompatible dimensions of result matrix", GSL_EBADLEN);
    }
  else if (work->size1 != L->size1)
    {
      GSL_ERROR("incompatible dimensions of work matrix", GSL_EBADLEN);
    }
  else if (df <= L->size1 - 1)
    {
      GSL_ERROR("incompatible degrees of freedom", GSL_EDOM);
    }
  else
    {
      /* result: X = L A A^T L^T */

      size_t d = L->size1, i, j;

      /* insure the upper part of A is zero before filling its lower part */
      gsl_matrix_set_zero(work);
      for (i = 0; i < d; ++i)
        {
          gsl_matrix_set(work, i, i, sqrt(gsl_ran_chisq(r, df - i)));

          for (j = 0; j < i; ++j)
            {
              gsl_matrix_set(work, i, j, gsl_ran_ugaussian(r));
            }
        }

      /* compute L * A */
      gsl_blas_dtrmm(CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, 1.0,
                     L, work);

      /* compute (L * A) * (L * A)^T */
      gsl_blas_dsyrk(CblasUpper, CblasNoTrans, 1.0, work, 0.0, result);
      for (i = 0; i < d; ++i)
        {
          for (j = 0; j < i; ++j)
            {
              gsl_matrix_set(result, i, j, gsl_matrix_get(result, j, i));
            }
        }

      return GSL_SUCCESS;
    }
}

/* Compute the log of the probability density function at a given quantile
 * matrix for a Wishart distribution using the Cholesky decomposition of the
 * scale matrix.
 *
 * X       quantile matrix at which to evaluate the PDF (dimension d x d)
 * L_X     matrix resulting from the Cholesky decomposition of
 *         of the quantile matrix at which to evaluate the PDF
 *         X = L_X L_X^T (dimension d x d)
 * df      degrees of freedom
 * L       matrix resulting from the Cholesky decomposition of
 *         the scale matrix V = L L^T (dimension d x d)
 * result  output of the density (dimension 1)
 * work    matrix used for intermediate computations (dimension d x d)
 */
int
gsl_ran_wishart_log_pdf (const gsl_matrix * X,
                         const gsl_matrix * L_X,
                         const double df,
                         const gsl_matrix * L,
                         double * result,
                         gsl_matrix * work)
{
  if (L->size1 != L->size2)
    {
      GSL_ERROR("L should be a square matrix", GSL_ENOTSQR);
    }
  else if (X->size1 != X->size2)
    {
      GSL_ERROR("X should be a square matrix", GSL_ENOTSQR);
    }
  else if (L_X->size1 != L_X->size2)
    {
      GSL_ERROR("L_X should be a square matrix", GSL_ENOTSQR);
    }
  else if (X->size1 != L->size1)
    {
      GSL_ERROR("incompatible dimensions of X matrix", GSL_EBADLEN);
    }
  else if (L_X->size1 != L->size1)
    {
      GSL_ERROR("incompatible dimensions of L_X matrix", GSL_EBADLEN);
    }
  else if (df <= L->size1 - 1)
    {
      GSL_ERROR("incompatible degrees of freedom", GSL_EDOM);
    }
  else
    {
      size_t d = L->size1, i;
      int status;
      double log_mv_Ga, log_det_V, log_det_X, tr_Vinv_X;

      /* compute the log of the multivariate Gamma */
      log_mv_Ga = d * (d-1) * 0.25 * log(M_PI);
      for (i = 0; i < d; ++i)
        {
          log_mv_Ga += gsl_sf_lngamma((df - i + 1) * 0.5);
        }

      /* compute the log of the determinant of the scale matrix */
      log_det_V = log(gsl_matrix_get(L, 0, 0));
      for (i = 1; i < d; ++i)
        {
          log_det_V += log(gsl_matrix_get(L, i, i));
        }
      log_det_V = 2 * log_det_V;

      /* compute the log of the determinant of the quantile matrix */
      log_det_X = log(gsl_matrix_get(L_X, 0, 0));
      for (i = 1; i < d; ++i)
        {
          log_det_X += log(gsl_matrix_get(L_X, i, i));
        }
      log_det_X = 2 * log_det_X;

      /* compute the trace of V^(-1) X */
      status = gsl_linalg_cholesky_solve_mat(L, X, work);
      if (status)
        return status;
      tr_Vinv_X = gsl_matrix_get(work, 0, 0);
      for (i = 1; i < d; ++i)
        {
          tr_Vinv_X += gsl_matrix_get(work, i, i);
        }

      /* add all to get the log of the PDF */
      *result = - (0.5 * df * d) * log(2.0)
                - (0.5 * df) * log_det_V
                - log_mv_Ga
                + 0.5 * (df - d - 1) * log_det_X
                - 0.5 * tr_Vinv_X;

      return GSL_SUCCESS;
    }
}

int
gsl_ran_wishart_pdf (const gsl_matrix * X,
                     const gsl_matrix * L_X,
                     const double df,
                     const gsl_matrix * L,
                     double * result,
                     gsl_matrix * work)
{
  double logpdf;
  int status = gsl_ran_wishart_log_pdf(X, L_X, df, L, &logpdf, work);

  if (status == GSL_SUCCESS)
    *result = exp(logpdf);

  return status;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense