text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include <nubspline.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>

extern void* xcalloc(int NumPoints, int size);

/*void eval_NUBspline_1d_d(NUBspline_1d_d * restrict spline, double x,
	double* restrict val);

void eval_NUBspline_1d_d_vg(NUBspline_1d_d * restrict spline, double x,
	double* restrict val, double* restrict grad);

void eval_NUBspline_1d_d_vgl(NUBspline_1d_d * restrict spline, double x,
	double* restrict val, double* restrict grad, double* restrict lapl);

void InterpolateWithSplines(int NumPointsToInterpolate, int NumPointsToEvaluate, double* x, double* y, double* NodalX, double* NodalY, double *func, double *NodalFunc, double* NodalDeriv)
{

	double *s = xcalloc(NumPointsToInterpolate, sizeof(double));
	double *NodalS = xcalloc(NumPointsToEvaluate, sizeof(double));
	for (int j = 0; j < NumPointsToInterpolate; j++)
	{
		s[i] = sqrt(pow(x[j]-x[0],2)+pow(y[j]-y[0],2));

	}

	for(int j = 0; j < NumPointsToEvaluate, j++)
	{
		NodalS[i] = sqrt(pow(NodalX[j] -x[0],2) + pow(NodalY[j]-y[0],2));
	}

	NUgrid* myGrid = create_general_grid(s, NumPointsToInterpolate);
	BCtype_d myBC;
	myBC.lCode = NATURAL;
	myBC.rCode = NATURAL;
	NUBspline_1d_d *mySpline = create_NUBspline_1d_d(myGrid, myBC, func);
	for(int i =0; i < NumPointsToEvaluate; i++)
	{
		if (NodalDeriv != NULL)
			eval_NUBspline_1d_d_vg(mySpline, NodalS[i], &NodalFunc[i], &NodalDeriv[i]);
		else
			eval_NUBspline_1d_d(mySpline, NodalS[i], &NodalFunc[i]);
	}

	free(s);
	free(NodalS);

}
*/

void InterpolateWithGSLSplines(int NumPointsToInterpolate, int NumPointsToEvaluate, double* x, double* y, double* NodalX, double* NodalY, double* func, double* NodalFunc, double* NodalDeriv)
{
	double *s = xcalloc(NumPointsToInterpolate, sizeof(double));
	double *NodalS = xcalloc(NumPointsToEvaluate, sizeof(double));
	for (int j = 0; j < NumPointsToInterpolate; j++)
	{
		s[j] = sqrt(pow(x[j]-x[0],2)+pow(y[j]-y[0],2));

	}

	for(int j = 0; j < NumPointsToEvaluate; j++)
	{
		NodalS[j] = sqrt(pow(NodalX[j] -x[0],2) + pow(NodalY[j]-y[0],2));
	}

	gsl_interp_accel *acc = gsl_interp_accel_alloc();
	gsl_spline* spline = gsl_spline_alloc(gsl_interp_cspline, NumPointsToInterpolate);
	gsl_spline_init(spline, s, func, NumPointsToInterpolate);

	for(int j = 0; j < NumPointsToEvaluate; j++)
	{
		NodalFunc[j] = gsl_spline_eval(spline, NodalS[j], acc);
		if (NodalDeriv != NULL)
			NodalDeriv[j] = gsl_spline_eval_deriv(spline, NodalS[j], acc);
	}

	gsl_interp_accel_free(acc);
	gsl_spline_free(spline);
	free(s);
	free(NodalS);

}

ext: c

max_stars_repo_name: evalseth/DG-RAIN

max_issues_repo_name: evalseth/DG-RAIN

max_forks_repo_name: evalseth/DG-RAIN


text:
/* eigen/schur.h
 *
 * Copyright (C) 2006 Patrick Alken
 *
 * 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_SCHUR_H__
#define __GSL_SCHUR_H__

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

/*
 * Prototypes
 */

void gsl_schur_standardize(gsl_matrix *T, size_t row, gsl_complex *eval1,
                           gsl_complex *eval2, int update_t, gsl_matrix *Z);

#endif /* __GSL_SCHUR_H__ */

ext: h

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef _POLIFITGSL_H
#define _POLIFITGSL_H
#include <gsl_multifit.h>
#include <stdbool.h>
#include <math.h>
bool polynomialfit(int obs, int degree,
		   double *dx, double *dy, double *store); /* n, p */
#endif

max_stars_repo_name: ruslankuzmin/julia

max_issues_repo_name: ruslankuzmin/julia

max_forks_repo_name: ruslankuzmin/julia


text:
#ifndef DART_POINTERS_H
#define DART_POINTERS_H

/*----- System Includes -----*/

#include <atomic>
#include <memory>
#include <gsl/gsl>
#include <stdint.h>
#include <functional>
#include <type_traits>

/*----- System Includes -----*/

#include "../refcount_traits.h"

/*----- Type Declarations -----*/

namespace dart {

  namespace detail {

    template <class T>
    struct managed_ptr_deleter {
      void operator ()(T ptr);
    };

    template <class T>
    struct managed_ptr_deleter<T[]> {
      void operator ()(std::remove_extent_t<T>* ptr);
    };

    // Type-erasure mechanism for counted_ptr_base.
    template <class Counter>
    struct managed_ptr {
      managed_ptr(void* ptr, int64_t use_count) : ptr(ptr), use_count(use_count) {}
      virtual ~managed_ptr() = default;
      virtual void destroy() = 0;

      void* ptr;
      Counter use_count;
    };
    template <class Counter, class PtrType, class Deleter>
    struct managed_ptr_eraser : managed_ptr<Counter> {
      managed_ptr_eraser(void const* ptr, int64_t use_count) : managed_ptr<Counter>(const_cast<void*>(ptr), use_count) {}
      managed_ptr_eraser(void const* ptr, int64_t use_count, Deleter const& deleter) :
        managed_ptr<Counter>(const_cast<void*>(ptr), use_count),
        deleter(deleter)
      {}
      managed_ptr_eraser(void const* ptr, int64_t use_count, Deleter&& deleter) :
        managed_ptr<Counter>(const_cast<void*>(ptr), use_count),
        deleter(std::move(deleter))
      {}
      virtual void destroy() override;

      Deleter deleter;
    };

    template <class T, class Counter>
    class counted_ptr_base {

      public:

        /*----- Public Types -----*/

        using element_type = T;
        using counter_type = Counter;
        template <class U>
        using ptr_type = std::conditional_t<std::is_array<U>::value, U, U*>;

        /*----- Lifecycle Functions -----*/

        counted_ptr_base() noexcept : value(nullptr) {}
        counted_ptr_base(std::nullptr_t) noexcept : value(nullptr) {}
        explicit counted_ptr_base(ptr_type<T> ptr) :
          value(new managed_ptr_eraser<counter_type, ptr_type<T>, managed_ptr_deleter<ptr_type<T>>>(ptr, 1))
        {}
        template <class Deleter>
        counted_ptr_base(ptr_type<T> ptr, Deleter&& deleter) :
          value(new managed_ptr_eraser<counter_type, ptr_type<T>, Deleter>(ptr, 1, std::move(deleter)))
        {}
        template <class U, class EnableIf =
          std::enable_if_t<
            std::is_convertible<ptr_type<U>, ptr_type<T>>::value
          >
        >
        counted_ptr_base(counted_ptr_base<U, counter_type> const& other) noexcept;
        template <class U, class EnableIf =
          std::enable_if_t<
            std::is_convertible<ptr_type<U>, ptr_type<T>>::value
          >
        >
        counted_ptr_base(counted_ptr_base<U, counter_type>&& other) noexcept;
        template <class U, class Deleter, class EnableIf =
          std::enable_if_t<
            std::is_same<ptr_type<U>, ptr_type<T>>::value
            ||
            std::is_convertible<ptr_type<U>, ptr_type<T>>::value
          >
        >
        counted_ptr_base(std::unique_ptr<U, Deleter>&& ptr);
        counted_ptr_base(counted_ptr_base const& other) noexcept;
        counted_ptr_base(counted_ptr_base&& other) noexcept;
        ~counted_ptr_base() noexcept;

        /*----- Operators -----*/

        // Assignment.
        template <class U, class EnableIf =
          std::enable_if_t<
            std::is_convertible<ptr_type<U>, ptr_type<T>>::value
          >
        >
        counted_ptr_base& operator =(counted_ptr_base<U, counter_type> const& other) noexcept;
        template <class U, class EnableIf =
          std::enable_if_t<
            std::is_convertible<ptr_type<U>, ptr_type<T>>::value
          >
        >
        counted_ptr_base& operator =(counted_ptr_base<U, counter_type>&& other) noexcept;
        counted_ptr_base& operator =(counted_ptr_base const& other) noexcept;
        counted_ptr_base& operator =(counted_ptr_base&& other) noexcept;
        template <class U, class Deleter, class EnableIf =
          std::enable_if_t<
            std::is_convertible<ptr_type<U>, ptr_type<T>>::value
          >
        >
        counted_ptr_base& operator =(std::unique_ptr<U, Deleter>&& ptr)
          noexcept(std::is_nothrow_move_constructible<Deleter>::value);

        // Comparison.
        bool operator ==(counted_ptr_base const& other) const noexcept;
        bool operator !=(counted_ptr_base const& other) const noexcept;
        bool operator <(counted_ptr_base const& other) const noexcept;
        bool operator <=(counted_ptr_base const& other) const noexcept;
        bool operator >(counted_ptr_base const& other) const noexcept;
        bool operator >=(counted_ptr_base const& other) const noexcept;

        // Conversion.
        explicit operator bool() const noexcept;

        /*----- Public API -----*/

        // Accessors.
        std::remove_extent_t<element_type>* get() const noexcept;
        bool unique() const noexcept;
        int64_t use_count() const noexcept;

        // Mutators.
        void reset() noexcept;
        template <class U, class EnableIf =
          std::enable_if_t<
            std::is_convertible<U, ptr_type<T>>::value
          >
        >
        void reset(U ptr);
        template <class U, class Deleter, class EnableIf =
          std::enable_if_t<
            std::is_convertible<U, ptr_type<T>>::value
          >
        >
        void reset(U ptr, Deleter&& deleter);

      protected:

        /*----- Protected Members -----*/

        gsl::owner<managed_ptr<counter_type>*> value;

        /*----- Friends -----*/

        template <class U, class C>
        friend class counted_ptr_base;

    };

    template <class T, class Counter>
    struct counted_ptr_impl : counted_ptr_base<T, Counter> {

      /*----- Types -----*/

      using typename counted_ptr_base<T, Counter>::element_type;

      /*----- Lifecycle Functions -----*/

      using counted_ptr_base<T, Counter>::counted_ptr_base;

      /*----- Operators -----*/

      // Dereference.
      auto operator *() const noexcept -> element_type&;
      auto operator ->() const noexcept -> element_type*;

      // Comparison.
      using counted_ptr_base<T, Counter>::operator ==;
      using counted_ptr_base<T, Counter>::operator !=;
      using counted_ptr_base<T, Counter>::operator <;
      using counted_ptr_base<T, Counter>::operator <=;
      using counted_ptr_base<T, Counter>::operator >;
      using counted_ptr_base<T, Counter>::operator >=;
      using counted_ptr_base<T, Counter>::operator bool;

    };

    template <class T, class Counter>
    struct counted_array_ptr_impl : counted_ptr_base<T, Counter> {

      /*----- Types -----*/

      using index_type = std::ptrdiff_t;
      using element_type = std::remove_extent_t<typename counted_ptr_base<T, Counter>::element_type>;

      /*----- Lifecycle Functions -----*/

      using counted_ptr_base<T, Counter>::counted_ptr_base;

      /*----- Operators -----*/

      // Derefence.
      element_type& operator [](index_type idx) const noexcept;

      // Comparison.
      using counted_ptr_base<T, Counter>::operator ==;
      using counted_ptr_base<T, Counter>::operator !=;
      using counted_ptr_base<T, Counter>::operator <;
      using counted_ptr_base<T, Counter>::operator <=;
      using counted_ptr_base<T, Counter>::operator >;
      using counted_ptr_base<T, Counter>::operator >=;
      using counted_ptr_base<T, Counter>::operator bool;

    };

  }

  template <class T>
  class shareable_ptr {

    public:

      /*----- Public Types -----*/

      using value_type = T;
      using element_type = typename refcount_traits<T>::element_type;

    private:

      struct is_nothrow_assignable :
        meta::conjunction<
          typename refcount_traits<T>::is_nothrow_copyable,
          typename refcount_traits<T>::is_nothrow_moveable
        >
      {};

      // XXX: MSVC gets confused if this is used directly in a noexcept specification.
      struct is_nothrow_default_constructible :
        std::integral_constant<bool,
          refcount_traits<T>::template can_nothrow_take<element_type*, std::default_delete<element_type>>::value
        >
      {};

    public:

      /*----- Lifecycle Functions -----*/

      // Converting constructors.
      shareable_ptr(std::nullptr_t)
          noexcept(is_nothrow_default_constructible::value) :
        shareable_ptr()
      {}
      template <class U, class D, class EnableIf =
        std::enable_if_t<
          refcount::can_take<T, U*, D>::value
        >
      >
      explicit shareable_ptr(std::unique_ptr<U, D>&& ptr)
          noexcept(refcount_traits<T>::template can_nothrow_take<U*, D>::value);
      shareable_ptr(T const& other) noexcept(refcount_traits<T>::is_nothrow_copyable::value);
      shareable_ptr(T&& other) noexcept(refcount_traits<T>::is_nothrow_moveable::value);

      // Ownership transfer constructors.
      template <class U, class =
        std::enable_if_t<
          refcount::can_take<T, U*, std::default_delete<U>>::value
        >
      >
      explicit shareable_ptr(U* owner)
          noexcept(refcount_traits<T>::template can_nothrow_take<U*, std::default_delete<U>>::value) :
        shareable_ptr(owner, std::default_delete<U> {})
      {}
      template <class U, class D, class EnableIf =
        std::enable_if_t<
          refcount::can_take<T, U*, D>::value
        >
      >
      explicit shareable_ptr(U* owner, D&& del)
          noexcept(refcount_traits<T>::template can_nothrow_take<U*, D>::value);

      // Special Constructors.
      shareable_ptr() noexcept(is_nothrow_default_constructible::value);
      shareable_ptr(shareable_ptr const& other) noexcept(refcount_traits<T>::is_nothrow_copyable::value);
      shareable_ptr(shareable_ptr&& other) noexcept(refcount_traits<T>::is_nothrow_moveable::value);
      ~shareable_ptr() noexcept;

      /*----- Operators -----*/

      // Assignment.
      template <class U, class D, class EnableIf =
        std::enable_if_t<
          refcount::can_take<T, U*, D>::value
        >
      >
      shareable_ptr& operator =(std::unique_ptr<U, D>&& ptr)
          noexcept(refcount_traits<T>::template can_nothrow_take<U*, D>::value);
      shareable_ptr& operator =(T const& other) noexcept(is_nothrow_assignable::value);
      shareable_ptr& operator =(T&& other) noexcept;
      shareable_ptr& operator =(std::nullptr_t) noexcept(is_nothrow_default_constructible::value);
      shareable_ptr& operator =(shareable_ptr const& other) noexcept(is_nothrow_assignable::value);
      shareable_ptr& operator =(shareable_ptr&& other) noexcept;

      // Dereference.
      auto operator *() const noexcept(refcount_traits<T>::is_nothrow_unwrappable::value) -> element_type&;
      auto operator ->() const noexcept(refcount_traits<T>::is_nothrow_unwrappable::value) -> element_type*;

      // Comparison.
      bool operator ==(shareable_ptr const& other) const noexcept(refcount_traits<T>::is_nothrow_unwrappable::value);
      bool operator !=(shareable_ptr const& other) const noexcept(refcount_traits<T>::is_nothrow_unwrappable::value);
      bool operator <(shareable_ptr const& other) const noexcept(refcount_traits<T>::is_nothrow_unwrappable::value);
      bool operator <=(shareable_ptr const& other) const noexcept(refcount_traits<T>::is_nothrow_unwrappable::value);
      bool operator >(shareable_ptr const& other) const noexcept(refcount_traits<T>::is_nothrow_unwrappable::value);
      bool operator >=(shareable_ptr const& other) const noexcept(refcount_traits<T>::is_nothrow_unwrappable::value);

      // Conversion.
      explicit operator bool() const noexcept(refcount_traits<T>::is_nothrow_unwrappable::value);

      /*----- Public API -----*/

      // Accessors.
      auto get() const noexcept(refcount_traits<T>::is_nothrow_unwrappable::value) -> element_type*;
      bool unique() const noexcept(refcount_traits<T>::has_nothrow_use_count::value);
      int64_t use_count() const noexcept(refcount_traits<T>::has_nothrow_use_count::value);

      // Mutators.
      void reset() noexcept(refcount_traits<T>::has_nothrow_reset::value);

      // Extremely awkward "copy-constructor" that allows us to copy our
      // wrapped value into a raw instance of T.
      // Function has to be written this way since we can't assume T is copyable/moveable
      // The only thing we can definitely do with T is pass it into dart::refcount_traits::copy.
      void share(T& ptr) const noexcept(refcount_traits<T>::is_nothrow_copyable::value);

      // Extremely awkward "move-constructor" that allows us to move our
      // wrapped value into a raw instance of T.
      // Function has to be written this way since we can't assume T is copyable/moveable
      // The only thing we can definitely do with T is pass it into dart::refcount_traits::move.
      void transfer(T& ptr) && noexcept(refcount_traits<T>::is_nothrow_moveable::value);

      auto raw() noexcept -> value_type&;
      auto raw() const noexcept -> value_type const&;

    private:

      /*----- Private Types -----*/

      struct partial_construction_tag {};

      /*----- Private Lifecycle Functions -----*/

      shareable_ptr(partial_construction_tag) {}

      /*----- Private Members -----*/

      union {
        value_type impl;
      };

      /*----- Friends -----*/

      template <class U, class... Args>
      friend shareable_ptr<U> make_shareable(Args&&...);

  };

  template <template <class> class RefCount>
  struct view_ptr_context {

    template <class T>
    class view_ptr {

      public:

        /*----- Public Types -----*/

        using refcount_type = RefCount<T>;
        using is_nonowning = refcount_type;
        using element_type = typename refcount_traits<refcount_type>::element_type;

        // XXX: Necessary to make logic in refcount_traits work.
        template <template <template <class> class> class Binder>
        using refcount_rebind = Binder<RefCount>;

        /*----- Lifecycle Functions -----*/

        // Functions must exist for refcount concept, but will throw if called
        // as view_ptr doesn't have ownership semantics
        explicit view_ptr(T*);
        template <class Del>
        explicit view_ptr(T*, Del&&);

        // I'm not currently disallowing construction from temporaries
        // I can imagine scenarios where it could be useful, and generally,
        // you'll need to know what you're doing to work with the view types
        view_ptr(std::nullptr_t) noexcept : view_ptr() {}
        view_ptr(refcount_type const& owner) noexcept : impl(&owner) {}

        // Lifecycle functions mostly do nothing
        view_ptr() noexcept : impl(nullptr) {}
        view_ptr(view_ptr const&) = default;
        view_ptr(view_ptr&& other) noexcept;
        ~view_ptr() = default;

        /*----- Operators -----*/

        auto operator =(refcount_type const& owner) noexcept -> view_ptr&;

        auto operator =(view_ptr const&) -> view_ptr& = default;
        auto operator =(view_ptr&& other) noexcept -> view_ptr&;

        auto operator *() const noexcept(refcount_traits<refcount_type>::is_nothrow_unwrappable::value) -> element_type&;
        auto operator ->() const noexcept(refcount_traits<refcount_type>::is_nothrow_unwrappable::value) -> element_type*;

        bool operator ==(view_ptr const& other) const noexcept(refcount_traits<refcount_type>::is_nothrow_unwrappable::value);
        bool operator !=(view_ptr const& other) const noexcept(refcount_traits<refcount_type>::is_nothrow_unwrappable::value);
        bool operator <(view_ptr const& other) const noexcept(refcount_traits<refcount_type>::is_nothrow_unwrappable::value);
        bool operator <=(view_ptr const& other) const noexcept(refcount_traits<refcount_type>::is_nothrow_unwrappable::value);
        bool operator >(view_ptr const& other) const noexcept(refcount_traits<refcount_type>::is_nothrow_unwrappable::value);
        bool operator >=(view_ptr const& other) const noexcept(refcount_traits<refcount_type>::is_nothrow_unwrappable::value);

        explicit operator bool() const noexcept(refcount_traits<refcount_type>::is_nothrow_unwrappable::value);

        operator refcount_type const&() const;

        /*----- Public API -----*/

        auto get() const noexcept(refcount_traits<refcount_type>::is_nothrow_unwrappable::value) -> element_type*;
        size_t use_count() const noexcept(refcount_traits<refcount_type>::has_nothrow_use_count::value);

        void reset() noexcept;

        auto raw() const noexcept -> refcount_type const&;

      private:

        /*----- Private Members -----*/

        refcount_type const* impl;

    };

  };

  template <class T>
  using unsafe_ptr = std::conditional_t<
    std::is_same<
      std::remove_extent_t<T>,
      T
    >::value,
    detail::counted_ptr_impl<T, int64_t>,
    detail::counted_array_ptr_impl<T, int64_t>
  >;

  template <class T, class EnableIf =
    std::enable_if_t<
      std::is_array<T>::value
    >
  >
  unsafe_ptr<T> make_unsafe(size_t idx);

  template <class T, class... Args>
  std::enable_if_t<!std::is_array<T>::value, unsafe_ptr<T>> make_unsafe(Args&&... the_args);

  template <class T>
  using skinny_ptr = std::conditional_t<
    std::is_same<
      std::remove_extent_t<T>,
      T
    >::value,
    detail::counted_ptr_impl<T, std::atomic<int64_t>>,
    detail::counted_array_ptr_impl<T, std::atomic<int64_t>>
  >;

  template <class T, class EnableIf =
    std::enable_if_t<
      std::is_array<T>::value
    >
  >
  skinny_ptr<T> make_skinny(size_t idx);

  template <class T, class... Args>
  std::enable_if_t<!std::is_array<T>::value, skinny_ptr<T>> make_skinny(Args&&... the_args);

}

/*----- Template Implementations -----*/

#include "ptrs.tcc"

#endif


text:
#ifndef COSMOLOGY_H
#define COSMOLOGY_H

#include <cmath>
//nbody code
#include <NBody.h>
#include <NBodyMath.h>
//gsl code
#include <gsl/gsl_integration.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

using namespace std;
using namespace Math;

namespace Cosmology
{
    //scale factor function
    Double_t ScaleFactorFunc(const Double_t a, const Double_t om, const Double_t xla);
    //{return 1/sqrt(1.+om*(1./a-1.)+xla*(pow(a,2.0)-1.));};

    //Hubble expansion factor, omegak is diff omega m omegala
    Double_t HubbleFunc(const Double_t a, const Double_t Ho, const Double_t om, const Double_t ola);
    //Generalized Hubble expansion factor with omega_r=0 and omega_k specified
    Double_t HubbleFunc(const Double_t a, const Double_t Ho, const Double_t om, const Double_t ode, const Double_t ok, Double_t alpha=0.0, Double_t wde=-1.0);

	//wrapper for integration
/*	Double_t HubbleIntFunc(Double_t a, void* params);
	double HubbleIntFuncGSL(double a, void* params);
	double HubbleIntFuncGSLMonte(double a, size_t dim, void* params);

	//wrapper for integration of 1/aH for time calculation.
	Double_t aHIntFunc(Double_t a, void* params);
	double aHIntFuncGSL(double a, void* params);
	double* aHIntFuncGSLMonte(double a, size_t dim, void* params);
*/
    //Evaluates the time t (in units of t0) and Hubble parameter (units 1/t0)
    //given the expansion factor a (=1 at t=1)
    //for cosmological parameters omega0 and lambda0 (evaluated at t=1).
    Double_t Timet(const Double_t a, const Double_t om, const Double_t ola);
    Double_t Timeh(const Double_t a, const Double_t om, const Double_t ola);

    //Growth factor
    Double_t GrowthFunc(const Double_t a, const Double_t om, const Double_t ola);
    Double_t GrowthFunc(const Double_t a, const Double_t Ho, const Double_t om, const Double_t ode, const Double_t ok, Double_t alpha=0.0, Double_t wde=-1.0);

	//window function
	Double_t WKR2(const Double_t k, const Double_t R);

	//Power Spectra
	Double_t PK(const Double_t k, const Double_t ns, const Double_t Amp);
	//----Based on Bardeen et al 1986 BBKS
	//qBBKS
	Double_t qBBSK(const Double_t k, const Double_t ns, const Double_t Omegam, const Double_t h, const Double_t Omegab);
	Double_t TBBKS(const Double_t k, const Double_t ns, const Double_t Omegam, const Double_t h, const Double_t Omegab);
	//PBBSK here k is in Mpc^-1
	Double_t PBBKS(const Double_t k, const Double_t ns, const Double_t Amp, const Double_t Omegam, const Double_t h, const Double_t Omegab=0.);

	//----Based on Eisenstien and Hu 1998
	//q98
	Double_t qEH98(const Double_t k, const Double_t Omegam, const Double_t h, const Double_t Theta27);
	//L98
	Double_t LEH98(const Double_t k, const Double_t Omegam, const Double_t h, const Double_t Theta27);
	//C98
	Double_t CEH98(const Double_t k, const Double_t Omegam, const Double_t h, const Double_t Theta27);
	//T98
	Double_t TEH98(const Double_t k, const Double_t Omegam, const Double_t h, const Double_t Theta27);
	//P98 here k is in Mpc^-1.
	Double_t PEH98(const Double_t k, const Double_t ns, const Double_t deltaH, const Double_t clight, const Double_t Ho, const Double_t Omegam, const Double_t h, const Double_t Theta27);
	//P98 here k is in Mpc^-1, but just pass amplitude in
	Double_t PEH98(const Double_t k, const Double_t ns, const Double_t Amp, const Double_t Omegam, const Double_t h, const Double_t Theta27);
	//----P98 with Warm Dark Matter dampenning like term, here k is in Mpc^-1.
	Double_t PWDMEH98(const Double_t k, const Double_t ns, const Double_t AMP, const Double_t Omegam, const Double_t h, const Double_t Theta27, const Double_t Rd);

	//effective index in k units of Mpc^-1. neff =dlnP/dlnk
	//diff L98 by k
	Double_t LEH98diff(const Double_t k, const Double_t Omegam, const Double_t h, const Double_t Theta27);
	//diff C98 by k
	Double_t CEH98diff(const Double_t k, const Double_t Omegam, const Double_t h, const Double_t Theta27);
	Double_t neffEH98(const Double_t k, const Double_t ns, const Double_t Omegam, const Double_t h, const Double_t Theta27);
	//----same as above but with Warm Dark Matter dampenning like term, here k is in Mpc^-1.
	Double_t neffWDMEH98(const Double_t k, const Double_t ns, const Double_t Omegam, const Double_t h, const Double_t Theta27, const Double_t Rd);

	//----Based on preprint of Green et al 04 for zeq>>z>zb and k>kb, zb=150 and kb~1e3 Mpc^-1.
	Double_t PGreen04(const Double_t k, const Double_t Amp, const Double_t wm);

	//----Based on data from a file not certain how to construct it yet (could use void * to store k and pk arrays.
	Double_t PFromFile(const Double_t k, void* params);

	// For integrals, void params has form of params[0]=ns;params[1]=Amp;params[2]=Omegam;params[3]=h;
	// Must eventually alter so that Omega_Lambda is not assumed to be 1-Omegam

	//integrate P(k)k^2dk*4pi/(2pi)^3
	Double_t IntegralkkPower(int powerspectype, void *params, int itype=0, bool logint=true);

	//integrate W^2(kR)P(k)k^2dk*4pi/(2pi)^3
	Double_t IntegralSigma(int powerspectype, Double_t R, void *params, int itype=0, bool logint=true);

    //-------------------------------------------------------------------
    //  Does open-ended romberg integration of the scale factor integral from a1 to a2
    //  Based on a Numerical Recipes routine
    Double_t aintegral(const Double_t a1,const Double_t a2, const Double_t omega, const Double_t xlambda);

}

#endif

max_stars_repo_name: ICRAR/NBodylib

max_issues_repo_name: ICRAR/NBodylib

max_forks_repo_name: ICRAR/NBodylib


text:
#include <config.h>
#include <stddef.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>

#include <gsl/gsl_dft_complex.h>
#include <gsl/gsl_dft_complex_float.h>

#include "complex_internal.h"

#define BASE_DOUBLE
#include "templates_on.h"
#include "dft_source.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "dft_source.c"
#include "templates_off.h"
#undef  BASE_FLOAT

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <errno.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include <signal.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>
#include <setjmp.h>

#include <glib.h>
#if GLIB_CHECK_VERSION (2, 6, 0)
#  include <glib/gstdio.h>
#endif
#include <gsl/gsl_spline.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_math.h>

#include "asf.h"
#include "float_image.h"
#include "asf_endian.h"
#include "asf_jpeg.h"

#include "asf_view.h"

#ifndef linux
#ifndef darwin
#ifndef win32
static double
round (double arg)
{
  return floor (arg + 0.5);
}
#endif // #ifndef win32
#endif // #ifndef darwin
#endif // #ifndef linux

#include "asf_glib.h"

// at 64MB tiles, this is ~1.5GB
static const int MAX_TILES = 24;

// quit blathering?
int quiet = FALSE;

static int data_size(CachedImage *self)
{
    switch (self->data_type) {
        case GREYSCALE_FLOAT:
            return 4;
        case RGB_BYTE:
            return 3;
        case GREYSCALE_BYTE:
            return 1;
        case RGB_FLOAT:
            return 12;
        default:
            assert(FALSE);
            return 0;
    }
}

static void print_cache_size(CachedImage *self)
{
    int i;
    int size=0;
    int ds = data_size(self);
    for (i=0; i<self->n_tiles; ++i)
        size += ds*self->rows_per_tile*self->ns;

    asfPrintStatus("Cache size is %.1f megabytes.\n",
        (float)size/1024./1024.);
}

static unsigned char *get_pixel(CachedImage *self, int line, int samp)
{
    // check if outside the image
    static unsigned char zero = 0;
    if (line<0 || samp<0 || line >= self->nl || samp >= self->ns)
        return &zero;

    // size of each pixel
    int ds = data_size(self);

    int i;
    for (i=0; i<self->n_tiles; ++i) {
        int rs = self->rowstarts[i];
        if (rs >= 0) {
            if (line >= rs && line < rs+self->rows_per_tile) {
                // found the right cache
                assert(self->cache[i]);

                // this probably won't ever happen, but here we go anyway
                if (self->n_access > 1024*1024*1024) {
                    asfPrintStatus("Resetting n_access.\n");
                    for (i=0; i<self->n_tiles; ++i)
                        self->access_counts[i] = 0;
                    self->n_access = 1;
                }

                // mark this as the most recently accessed
                self->access_counts[i] = self->n_access++;

                // return pointer to the cached value
                return &self->cache[i][((line-rs)*self->ns + samp)*ds];
            }
        }
    }

    int spot = 0;
    if (!self->reached_max_tiles) {
        assert(self->cache[self->n_tiles] == NULL);
        unsigned char *data = malloc(ds*self->ns*self->rows_per_tile);
        if (!data) {
            // if this is the first tile -- abort, we are out of memory
            if (self->n_tiles == 0)
                asfPrintError("Failed to allocate cache of %ld bytes.\n"
                              "Out of memory.\n",
                              ds*self->ns*self->rows_per_tile);
            // couldn't allocate the next tile -- must dump existing
            if (!quiet)
                asfPrintStatus("reached max # of tiles: %d\n", self->n_tiles);
            print_cache_size(self);
            self->reached_max_tiles = TRUE;
        } else {
            spot = self->n_tiles;
            self->cache[spot] = data;
            ++self->n_tiles;
        }
    }

    if (self->reached_max_tiles) {
        // dump an existing cached tile
        // not found in the cache -- find least used spot
        int least_access_count = self->access_counts[0];
        for (i=0; i<self->n_tiles; ++i) {
            if (self->access_counts[i] < least_access_count) {
                least_access_count = self->access_counts[i];
                spot = i;
            }
        }
    }

    if (!self->reached_max_tiles && self->n_tiles == MAX_TILES) {
        if (!quiet)
            asfPrintStatus("Fully loaded with %d tiles.\n", self->n_tiles);
        print_cache_size(self);
        self->reached_max_tiles = TRUE;
    }

    // load info from file
    assert(spot >= 0 && spot < self->n_tiles);
    assert(self->cache[spot] != NULL);

    // clear out the cache -- we may not fill up the tile, if
    // we are near the end of the file, and we don't want old data
    // to appear
    int ns = self->ns;
    memset(self->cache[spot], 0, ds*ns*self->rows_per_tile);

    // update where this cache entry starts
    int rs = (line / self->rows_per_tile) * self->rows_per_tile;
    self->rowstarts[spot] = rs;

    // mark this tile as the most recently accessed
    self->access_counts[spot] = self->n_access++;

    //printf("Updated entry #%d:\n"
    //    "     row start: %d\n"
    //    "     row end: %d\n"
    //    "     access count: %d\n", spot,
    //    rs, rs+self->rows_per_tile-1, self->access_counts[spot]);

    // ensure we don't read past the end of the file
    int rows_to_get = self->rows_per_tile;
    if (rs + self->rows_per_tile > self->nl)
        rows_to_get = self->nl - rs;

    if (!quiet) {
        asfPrintStatus("Cache: loading into spot #%d: rows %d-%d\n",
            spot, rs, rs+rows_to_get);
        //print_cache_size(self);
    }

    self->client->read_fn(rs, rows_to_get, (void*)(self->cache[spot]),
        self->client->read_client_info, self->meta, self->client->data_type);

    assert((line-rs)*self->ns + samp <= self->ns*self->rows_per_tile);
    return &self->cache[spot][((line-rs)*self->ns + samp)*ds];
}

void load_thumbnail_data(CachedImage *self, int thumb_size_x, int thumb_size_y,
                         void *dest_void)
{
    if (self->entire_image_fits || !self->client->thumb_fn) {
        // Either we don't have thumbnailing support from the client,
        // or the image will fit entirely in memory.  In both cases, we
        // can just call get_pixel() on the subset necessary to show
        // the image, which will load the entire image into cache.
        // (... well, unless it is so huge that our tiles fit between
        //      the thumbnail grid ... but that's crazy talk)
        int ds = data_size(self);
        unsigned char *dest = (unsigned char*)dest_void;

        // this will fill the cache with the image data
        int sf = self->meta->general->line_count / thumb_size_y;
        //assert(sf == self->meta->general->sample_count / thumb_size_x);

        // supress the "populating cache" msgs when loading the whole thing
        quiet=TRUE;

        int i,j;
        for (i=0; i<thumb_size_y; ++i) {
            for (j=0; j<thumb_size_x; ++j) {
                // make this independent of the data type
                unsigned char *p = get_pixel(self,i*sf,j*sf);
                memcpy(dest+(i*thumb_size_x+j)*ds, p, ds);
            }
            asfPercentMeter((float)i/(thumb_size_y-1));
        }

        quiet=FALSE;
    } else {
        self->client->thumb_fn(thumb_size_x, thumb_size_y,
            self->meta, self->client->read_client_info, dest_void,
            self->client->data_type);
    }
}

CachedImage * cached_image_new_from_file(
    const char *file, meta_parameters *meta, ClientInterface *client,
    ImageStats *stats, ImageStatsRGB *stats_r, ImageStatsRGB *stats_g,
    ImageStatsRGB *stats_b)
{
    CachedImage *self = MALLOC(sizeof(CachedImage));

    asfPrintStatus("Opening cache: %s\n", file);

    self->data_type = client->data_type;
    assert(self->data_type != UNDEFINED);

    self->client = client;     // take ownership of this
    self->meta = meta;         // do NOT take ownership of this

    self->stats = stats;       // do NOT take ownership of this

    self->stats_r = stats_r;   // do NOT take ownership of this
    self->stats_g = stats_g;   // do NOT take ownership of this
    self->stats_b = stats_b;   // do NOT take ownership of this

    // line line_count may have been fudges, if we are multilooking
    self->nl = meta->general->line_count;
    self->ns = meta->general->sample_count;
    asfPrintStatus("Image is %dx%d LxS\n", self->nl, self->ns);

    if (client->require_full_load) {
        // Use only 1 tile -- load entire image into it
        // (client tells us we should do it this way)
        // will it fit?  Who knows.  Assume it will, MALLOC will fail
        // if it actually does not.
        self->rows_per_tile = self->nl;
    } else {
        // how many rows per tile?
        // We will use ~64 Meg tiles
        self->rows_per_tile = 64*1024*1024 / (self->ns*data_size(self));

        // test line -- uncomment this for very small tiles
        //self->rows_per_tile = 2*1024*1024 / (self->ns*data_size(self));
    }

    asfPrintStatus("Using %d rows per tile.\n", self->rows_per_tile);

    int n_tiles_required = (int)ceil((double)self->nl / self->rows_per_tile);
    self->entire_image_fits = n_tiles_required <= MAX_TILES;
    // self->entire_image_fits = FALSE; // uncomment to test thumb_fn

    // at the beginning, we have no tiles
    self->n_tiles = 0;
    self->reached_max_tiles = FALSE;

    int i;
    self->rowstarts = MALLOC(sizeof(int)*MAX_TILES);
    self->cache = MALLOC(sizeof(float*)*MAX_TILES);
    self->access_counts = MALLOC(sizeof(int)*MAX_TILES);
    for (i=0; i<MAX_TILES; ++i) {
        self->rowstarts[i] = -1;
        self->cache[i] = NULL;
        self->access_counts[i] = 0;
    }

    self->n_access = 0;

    asfPrintStatus("Number of tiles required for the entire image: %d\n",
        n_tiles_required);
    asfPrintStatus("Fits in memory: %s\n",
        self->entire_image_fits ? "Yes" : "No");

    return self;
}

float cached_image_get_pixel (CachedImage *self, int line, int samp)
{
    if (self->data_type == GREYSCALE_FLOAT) {
        return *((float*)get_pixel(self, line, samp));
    }
    else if (self->data_type == GREYSCALE_BYTE) {
        return (float) *(get_pixel(self, line, samp));
    }
    else if (self->data_type == RGB_BYTE || self->data_type == RGB_FLOAT) {
        unsigned char r, g, b;
        cached_image_get_rgb(self, line, samp, &r, &g, &b);
        return ((float)r + (float)g + (float)b)/3.;
    }

    // not reached
    assert(0);
    return 0;
}

void cached_image_get_rgb(CachedImage *self, int line, int samp,
                          unsigned char *r, unsigned char *g,
                          unsigned char *b)
{
    if (self->data_type == GREYSCALE_FLOAT) {
        float f = cached_image_get_pixel(self, line, samp);
        if (have_lut()) {
            // do not scale in the case of a lut
            apply_lut((int)f, r, g, b);
        } else {
            *r = *g = *b =
              (unsigned char)calc_scaled_pixel_value(self->stats, f);
        }
    }
    else if (self->data_type == GREYSCALE_BYTE) {
        if (have_lut()) {
            apply_lut((int)(*(get_pixel(self, line, samp))), r, g, b);
        }
        else {
            float f = (float) *(get_pixel(self, line, samp));
            *r = *g = *b =
                (unsigned char)calc_scaled_pixel_value(self->stats, f);
        }
    }
    else if (self->data_type == RGB_BYTE) {
        unsigned char *uc = get_pixel(self, line, samp);

        *r = (unsigned char)calc_rgb_scaled_pixel_value(self->stats_r,
                                                        (float)uc[0]);
        *g = (unsigned char)calc_rgb_scaled_pixel_value(self->stats_g,
                                                        (float)uc[1]);
        *b = (unsigned char)calc_rgb_scaled_pixel_value(self->stats_b,
                                                        (float)uc[2]);
    }
    else if (self->data_type == RGB_FLOAT) {
        float *f = (float*)get_pixel(self, line, samp);

        *r = (unsigned char)calc_rgb_scaled_pixel_value(self->stats_r,f[0]);
        *g = (unsigned char)calc_rgb_scaled_pixel_value(self->stats_g,f[1]);
        *b = (unsigned char)calc_rgb_scaled_pixel_value(self->stats_b,f[2]);
    }
    else {
        // impossible!
        assert(0);
        *r = *g = *b = 0;
    }
}

void cached_image_get_rgb_float(CachedImage *self, int line, int samp,
                                float *r, float *g, float *b)
{
    if (self->data_type == GREYSCALE_FLOAT) {
        float f = cached_image_get_pixel(self, line, samp);
        *r = *g = *b = calc_scaled_pixel_value(self->stats, f);
    }
    else if (self->data_type == GREYSCALE_BYTE) {
        *r = *g = *b = (float)(*(get_pixel(self, line, samp)));
    }
    else if (self->data_type == RGB_BYTE) {
        unsigned char *uc = get_pixel(self, line, samp);

        *r = (float)(uc[0]);
        *g = (float)(uc[1]);
        *b = (float)(uc[2]);
    }
    else if (self->data_type == RGB_FLOAT) {
        float *f = (float*)get_pixel(self, line, samp);

        *r = f[0];
        *g = f[1];
        *b = f[2];
    }
    else {
        // impossible!
        assert(0);
        *r = *g = *b = 0.;
    }
}

void cached_image_free (CachedImage *self)
{
    int i;
    for (i=0; i<self->n_tiles; ++i) {
        if (self->cache[i])
            free(self->cache[i]);
    }

    if (self->client->free_fn)
      self->client->free_fn(self->client->read_client_info);

    free(self->rowstarts);
    free(self->access_counts);
    free(self->cache);
    free(self->client);

    // we do not own the metadata -- don't free it!

    free(self);
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
/* addmatrix_gsl.c */
/*
  Program to add two variable dimension matrices input by user
*/
// https://www.gnu.org/software/gsl/manual/html_node/Example-programs-for-matrices.html

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_matrix.h>

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

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

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

  gsl_matrix_free(m);

  return 0;
}

max_stars_repo_name: ernestyalumni/CompPhys

max_issues_repo_name: ernestyalumni/CompPhys

max_forks_repo_name: ernestyalumni/CompPhys


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2018-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */
#pragma once

#include <gsl/gsl-lite.hpp>
#include <platform/sized_buffer.h>
#include <array>
#include <optional>
#include <type_traits>

namespace cb::mcbp {

/**
 * Helper code for encode and decode of LEB128 values.
 * - mcbp encodes collection-ID as an unsigned LEB128
 * - see https://en.wikipedia.org/wiki/LEB128
 */

// Empty, non specialised version of the decoder class
template <class T, class Enable = void>
class unsigned_leb128 {};

/**
 * For encoding a unsigned T leb128, class constructs from a T value and
 * provides a const_byte_buffer for access to the encoded
 */
template <class T>
class unsigned_leb128<
        T,
        typename std::enable_if<std::is_unsigned<T>::value>::type> {
public:
    explicit unsigned_leb128(T in) {
        while (in > 0) {
            auto byte = gsl::narrow_cast<uint8_t>(in & 0x7full);
            in >>= 7;

            // In has more data?
            if (in > 0) {
                byte |= 0x80;
                encodedData[encodedSize - 1] = byte;
                // Increase the size
                encodedSize++;
            } else {
                encodedData[encodedSize - 1] = byte;
            }
        }
    }

    cb::const_byte_buffer get() const {
        return {encodedData.data(), encodedSize};
    }

    const uint8_t* begin() const {
        return encodedData.data();
    }

    const uint8_t* end() const {
        return encodedData.data() + encodedSize;
    }

    const uint8_t* data() const {
        return encodedData.data();
    }

    size_t size() const {
        return encodedSize;
    }

    constexpr static size_t getMaxSize() {
        return maxSize;
    }

    /**
     * decode returns the decoded T and a const_byte_buffer initialised with the
     * data following the leb128 data.
     *
     * @param buf buffer containing a leb128 encoded value (of size T). This can
     *            be a prefix on some other data, the decode will only process
     *            up to the maximum number of bytes permitted for the type T.
     *            E.g. uint32_t use 5 bytes maximum. buf.size must be >= 1
     *
     * @returns A std::pair where first is the decoded value and second is a
     *          buffer initialised with the data following the leb128 data. Note
     *          if the input buf was 100% only a leb128, the returned buffer
     *          will point outside of the input buf, but size will be 0.
     *
     * @throws std::invalid_argument if the input is not a valid leb128, this
     *         means decode processed 'getMaxSize' bytes without a stop byte.
     */
    static std::pair<T, cb::const_byte_buffer> decode(
            cb::const_byte_buffer buf);

    /**
     * decodeCanonical returns the decoded value of type T and a
     * const_byte_buffer initialised with the data following the leb128 data.
     *
     * This version does not throw an exception, but returns failure for two
     * reasons.
     *  - no-stop byte found
     *  - non-canonical encoding was used, e.g. 0x81.00 instead of 0x01
     *
     * The caller will have to inspect the input to determine the error.
     *
     * @param buf buffer containing a leb128 encoded value (of size T). This can
     *            be a prefix on some other data, the decode will only process
     *            up to the maximum number of bytes permitted for the type T.
     *            E.g. uint32_t use 5 bytes maximum. buf.size must be >= 1
     *
     * @returns On error a std::pair where the second buffer has a nullptr and
     *          zero size, first is set to 0. On success a std::pair where first
     *          is the decoded value and second is a buffer initialised with the
     *          data following the leb128 data. Note if the input buf was 100%
     *          only a leb128, the returned buffer will point outside of the
     *          input buf, but size will be 0.
     */
    static std::pair<T, cb::const_byte_buffer> decodeCanonical(
            cb::const_byte_buffer buf);

    /**
     * decodeNoThrow returns the decoded value of type T and a const_byte_buffer
     * initialised with the data following the leb128 data.
     *
     * This version does not throw an exception, but returns failure if no-stop
     * is byte found.
     *
     * @param buf buffer containing a leb128 encoded value (of size T). This can
     *            be a prefix on some other data, the decode will only process
     *            up to the maximum number of bytes permitted for the type T.
     *            E.g. uint32_t use 5 bytes maximum. buf.size must be >= 1
     *
     * @returns On error a std::pair where the second buffer has a nullptr and
     *          zero size, first is set to 0. On success a std::pair where first
     *          is the decoded value and second is a buffer initialised with the
     *          data following the leb128 data. Note if the input buf was 100%
     *          only a leb128, the returned buffer will point outside of the
     *          input buf, but size will be 0.
     */
    static std::pair<T, cb::const_byte_buffer> decodeNoThrow(
            cb::const_byte_buffer buf);

protected:
    struct NoThrow {};
    /**
     * decode returns the decoded value of type T and a const_byte_buffer
     * initialised with the data following the leb128 data.
     *
     * This is the protected inner method used by the public decode, does not
     * throw for bad input (allows public methods to decide error fate)
     */
    static std::pair<T, cb::const_byte_buffer> decode(cb::const_byte_buffer buf,
                                                      NoThrow);

    template <typename, typename>
    friend class unsigned_leb128;

    /**
     * Test that a decoded value was encoded in the canonical format.
     *
     * The test works by examining the length and comparing against a constant.
     * The constant is the maximum value that can be encoded as leb128 in
     * 'encodedLength - 1' bytes.
     *
     * For example if the encodedLength was 2 and the value was less than or
     * equal to 127, a non-canonical encoding was used, 127 and less can and
     * must be encoded in only 1 byte.
     *
     * So the test when encoded length is 2 is that the value is greater than
     * 127. If the encoded length is 3 the value must be greater than 16383 and
     * so on.
     *
     * @param value The integer that was decoded
     * @param encodedLength How many bytes the leb128 encoding used
     */
    static inline bool is_canonical(uint64_t value, size_t encodedLength);

private:
    // Larger T may need a larger array
    static_assert(sizeof(T) <= 8, "Class is only valid for uint 8/16/64");

    // value is large enough to store ~0 as leb128
    static constexpr size_t maxSize = sizeof(T) + (((sizeof(T) + 1) / 8) + 1);
    std::array<uint8_t, maxSize> encodedData{};
    uint8_t encodedSize{1};
};

// Generate the maximum value that can be encoded in nbytes
#define MAX_LEB128(nbytes) \
    ((0x7full << ((nbytes - 1) * 7)) | ((1ull << ((nbytes - 1) * 7)) - 1ull))

template <>
inline bool unsigned_leb128<uint8_t>::is_canonical(uint64_t value,
                                                   size_t encodedLength) {
    return (encodedLength == 2 && value > MAX_LEB128(1)) || encodedLength == 1;
}

template <>
inline bool unsigned_leb128<uint16_t>::is_canonical(uint64_t value,
                                                    size_t encodedLength) {
    if (unsigned_leb128<uint8_t>::is_canonical(value, encodedLength)) {
        return true;
    }
    return encodedLength == 3 && value > MAX_LEB128(2);
}

template <>
inline bool unsigned_leb128<uint32_t>::is_canonical(uint64_t value,
                                                    size_t encodedLength) {
    if (unsigned_leb128<uint16_t>::is_canonical(value, encodedLength)) {
        return true;
    }

    switch (encodedLength) {
    case 4:
        return value > MAX_LEB128(3);
    case 5:
        return value > MAX_LEB128(4);
    }

    return false;
}

template <>
inline bool unsigned_leb128<uint64_t>::is_canonical(uint64_t value,
                                                    size_t encodedLength) {
    // We first have to ask if this is non-canonical for the lower size, e.g.
    // u32. Each size asks the lower size first.
    if (unsigned_leb128<uint32_t>::is_canonical(value, encodedLength)) {
        return true;
    }

    switch (encodedLength) {
    case 6:
        return value > MAX_LEB128(5);
    case 7:
        return value > MAX_LEB128(6);
    case 8:
        return value > MAX_LEB128(7);
    case 9:
        return value > MAX_LEB128(8);
    case 10:
        return value > MAX_LEB128(9);
    }

    return false;
}

template <class T>
std::pair<T, cb::const_byte_buffer>
unsigned_leb128<T, typename std::enable_if<std::is_unsigned<T>::value>::type>::
        decode(cb::const_byte_buffer buf, NoThrow) {
    T rv = buf[0] & 0x7full;
    size_t end = 0;
    // Process up to the end of buf, or the max size for T, this ensures that
    // bad input, e.g. no stop-byte avoids invalid shifts (where shift just
    // keeps getting better). Primarily this gives us much better control over
    // invalid input, e.g. 20 bytes of 0x80 with a stop byte, would of
    // previously decoded to 0, but is really not valid input.
    size_t size =
            std::min(buf.size(), cb::mcbp::unsigned_leb128<T>::getMaxSize());
    if ((buf[0] & 0x80) == 0x80ull) {
        T shift = 7;
        // shift in the remaining data
        for (end = 1; end < size; end++) {
            rv |= (buf[end] & 0x7full) << shift;
            if ((buf[end] & 0x80ull) == 0) {
                break; // no more
            }
            shift += 7;
        }
        // We should of stopped for a stop byte, not the end of the buffer or
        // max encoding
        if (end == size) {
            return {0, cb::const_byte_buffer{}};
        }
    }
    // Return the decoded value and a buffer for any remaining data
    return {rv,
            cb::const_byte_buffer{buf.data() + end + 1,
                                  buf.size() - (end + 1)}};
}

template <class T>
std::pair<T, cb::const_byte_buffer>
unsigned_leb128<T, typename std::enable_if<std::is_unsigned<T>::value>::type>::
        decode(cb::const_byte_buffer buf) {
    if (buf.size() > 0) {
        auto rv = unsigned_leb128<T>::decode(buf, NoThrow{});
        if (rv.second.data()) {
            return rv;
        }
    }
    throw std::invalid_argument(
            "`unsigned_leb128::decode invalid leb128 of size:" +
            std::to_string(buf.size()));
}

template <class T>
std::pair<T, cb::const_byte_buffer>
unsigned_leb128<T, typename std::enable_if<std::is_unsigned<T>::value>::type>::
        decodeCanonical(cb::const_byte_buffer buf) {
    auto rv = unsigned_leb128<T>::decode(buf, NoThrow{});

    if (rv.second.data() &&
        !is_canonical(rv.first, size_t(rv.second.data() - buf.data()))) {
        return {0, cb::const_byte_buffer{}};
    }

    return rv;
}

template <class T>
std::pair<T, cb::const_byte_buffer>
unsigned_leb128<T, typename std::enable_if<std::is_unsigned<T>::value>::type>::
        decodeNoThrow(cb::const_byte_buffer buf) {
    return unsigned_leb128<T>::decode(buf, NoThrow{});
}

/**
 * @return a buffer to the data after the leb128 prefix
 */
template <class T>
typename std::enable_if<std::is_unsigned<T>::value, cb::const_byte_buffer>::type
skip_unsigned_leb128(cb::const_byte_buffer buf) {
    return unsigned_leb128<T>::decode(buf).second;
}

/// @return the index of the stop byte within buf
static inline std::optional<size_t> unsigned_leb128_get_stop_byte_index(
        cb::const_byte_buffer buf) {
    // If buf does not contain a stop-byte, invalid
    size_t stopByte = 0;
    for (auto c : buf) {
        if ((c & 0x80ull) == 0) {
            return stopByte;
        }
        stopByte++;
    }
    return {};
}

} // namespace cb::mcbp


text:
/* ode-initval/gsl_odeiv.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Author:  G. Jungman
 */
#ifndef __GSL_ODEIV_H__
#define __GSL_ODEIV_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_types.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

/* Description of a system of ODEs.
 *
 * y' = f(t,y) = dydt(t, y)
 *
 * The system is specified by giving the right-hand-side
 * of the equation and possibly a jacobian function.
 *
 * Some methods require the jacobian function, which calculates
 * the matrix dfdy and the vector dfdt. The matrix dfdy conforms
 * to the GSL standard, being a continuous range of floating point
 * values, in row-order.
 *
 * As with GSL function objects, user-supplied parameter
 * data is also present.
 */

typedef struct
{
  int (* function) (double t, const double y[], double dydt[], void * params);
  int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);
  size_t dimension;
  void * params;
}
gsl_odeiv_system;

#define GSL_ODEIV_FN_EVAL(S,t,y,f)  (*((S)->function))(t,y,f,(S)->params)
#define GSL_ODEIV_JA_EVAL(S,t,y,dfdy,dfdt)  (*((S)->jacobian))(t,y,dfdy,dfdt,(S)->params)

/* General stepper object.
 *
 * Opaque object for stepping an ODE system from t to t+h.
 * In general the object has some state which facilitates
 * iterating the stepping operation.
 */

typedef struct
{
  const char * name;
  int can_use_dydt_in;
  int gives_exact_dydt_out;
  void * (*alloc) (size_t dim);
  int  (*apply)  (void * state, size_t dim, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * dydt);
  int  (*reset) (void * state, size_t dim);
  unsigned int  (*order) (void * state);
  void (*free)  (void * state);
}
gsl_odeiv_step_type;

typedef struct {
  const gsl_odeiv_step_type * type;
  size_t dimension;
  void * state;
}
gsl_odeiv_step;

/* Available stepper types.
 *
 * rk2    : embedded 2nd(3rd) Runge-Kutta
 * rk4    : 4th order (classical) Runge-Kutta
 * rkck   : embedded 4th(5th) Runge-Kutta, Cash-Karp
 * rk8pd  : embedded 8th(9th) Runge-Kutta, Prince-Dormand
 * rk2imp : implicit 2nd order Runge-Kutta at Gaussian points
 * rk4imp : implicit 4th order Runge-Kutta at Gaussian points
 * gear1  : M=1 implicit Gear method
 * gear2  : M=2 implicit Gear method
 */

GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2;
GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk4;
GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rkf45;
GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rkck;
GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk8pd;
GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2imp;
GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2simp;
GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk4imp;
GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_bsimp;
GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_gear1;
GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_gear2;

/* Constructor for specialized stepper objects.
 */
GSL_FUN gsl_odeiv_step * gsl_odeiv_step_alloc(const gsl_odeiv_step_type * T, size_t dim);
GSL_FUN int  gsl_odeiv_step_reset(gsl_odeiv_step * s);
GSL_FUN void gsl_odeiv_step_free(gsl_odeiv_step * s);

/* General stepper object methods.
 */
GSL_FUN const char * gsl_odeiv_step_name(const gsl_odeiv_step * s);
GSL_FUN unsigned int gsl_odeiv_step_order(const gsl_odeiv_step * s);

GSL_FUN int  gsl_odeiv_step_apply(gsl_odeiv_step * s, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * dydt);

/* General step size control object.
 *
 * The hadjust() method controls the adjustment of
 * step size given the result of a step and the error.
 * Valid hadjust() methods must return one of the codes below.
 *
 * The general data can be used by specializations
 * to store state and control their heuristics.
 */

typedef struct
{
  const char * name;
  void * (*alloc) (void);
  int  (*init) (void * state, double eps_abs, double eps_rel, double a_y, double a_dydt);
  int  (*hadjust) (void * state, size_t dim, unsigned int ord, const double y[], const double yerr[], const double yp[], double * h);
  void (*free) (void * state);
}
gsl_odeiv_control_type;

typedef struct
{
  const gsl_odeiv_control_type * type;
  void * state;
}
gsl_odeiv_control;

/* Possible return values for an hadjust() evolution method.
 */
#define GSL_ODEIV_HADJ_INC   1  /* step was increased */
#define GSL_ODEIV_HADJ_NIL   0  /* step unchanged     */
#define GSL_ODEIV_HADJ_DEC (-1) /* step decreased     */

GSL_FUN gsl_odeiv_control * gsl_odeiv_control_alloc(const gsl_odeiv_control_type * T);
GSL_FUN int gsl_odeiv_control_init(gsl_odeiv_control * c, double eps_abs, double eps_rel, double a_y, double a_dydt);
GSL_FUN void gsl_odeiv_control_free(gsl_odeiv_control * c);
GSL_FUN int gsl_odeiv_control_hadjust (gsl_odeiv_control * c, gsl_odeiv_step * s, const double y[], const double yerr[], const double dydt[], double * h);
GSL_FUN const char * gsl_odeiv_control_name(const gsl_odeiv_control * c);

/* Available control object constructors.
 *
 * The standard control object is a four parameter heuristic
 * defined as follows:
 *    D0 = eps_abs + eps_rel * (a_y |y| + a_dydt h |y'|)
 *    D1 = |yerr|
 *    q  = consistency order of method (q=4 for 4(5) embedded RK)
 *    S  = safety factor (0.9 say)
 *
 *                      /  (D0/D1)^(1/(q+1))  D0 >= D1
 *    h_NEW = S h_OLD * |
 *                      \  (D0/D1)^(1/q)      D0 < D1
 *
 * This encompasses all the standard error scaling methods.
 *
 * The y method is the standard method with a_y=1, a_dydt=0.
 * The yp method is the standard method with a_y=0, a_dydt=1.
 */

GSL_FUN gsl_odeiv_control * gsl_odeiv_control_standard_new(double eps_abs, double eps_rel, double a_y, double a_dydt);
GSL_FUN gsl_odeiv_control * gsl_odeiv_control_y_new(double eps_abs, double eps_rel);
GSL_FUN gsl_odeiv_control * gsl_odeiv_control_yp_new(double eps_abs, double eps_rel);

/* This controller computes errors using different absolute errors for
 * each component
 *
 *    D0 = eps_abs * scale_abs[i] + eps_rel * (a_y |y| + a_dydt h |y'|)
 */
GSL_FUN gsl_odeiv_control * gsl_odeiv_control_scaled_new(double eps_abs, double eps_rel, double a_y, double a_dydt, const double scale_abs[], size_t dim);

/* General evolution object.
 */
typedef struct {
  size_t dimension;
  double * y0;
  double * yerr;
  double * dydt_in;
  double * dydt_out;
  double last_step;
  unsigned long int count;
  unsigned long int failed_steps;
}
gsl_odeiv_evolve;

/* Evolution object methods.
 */
GSL_FUN gsl_odeiv_evolve * gsl_odeiv_evolve_alloc(size_t dim);
GSL_FUN int gsl_odeiv_evolve_apply(gsl_odeiv_evolve * e, gsl_odeiv_control * con, gsl_odeiv_step * step, const gsl_odeiv_system * dydt, double * t, double t1, double * h, double y[]);
GSL_FUN int gsl_odeiv_evolve_reset(gsl_odeiv_evolve * e);
GSL_FUN void gsl_odeiv_evolve_free(gsl_odeiv_evolve * e);

__END_DECLS

#endif /* __GSL_ODEIV_H__ */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/******************************************************************************
CosmoLike Configuration Space Covariances for Projected Galaxy 2-Point Statistics
https://github.com/CosmoLike/CosmoCov
by CosmoLike developers
******************************************************************************/

#include <math.h>
#include <stdlib.h>
#if !defined(__APPLE__)
#include <malloc.h>
#endif
#include <stdio.h>
#include <assert.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

#include <fftw3.h>

#include "../cosmolike_core/2dfftlog/twobessel.h"
#include "../cosmolike_core/2dfftlog/utils.h"

#include "../cosmolike_core/theory/basics.c"
#include "../cosmolike_core/theory/structs.c"
#include "../cosmolike_core/theory/recompute.c"
#include "../cosmolike_core/theory/cosmo3D.c"
#include "../cosmolike_core/theory/redshift_spline.c"
#include "../cosmolike_core/theory/halo.c"
#include "../cosmolike_core/theory/cosmo2D_fourier.c"

#include "../cosmolike_core/theory/covariances_3D.c"
#include "../cosmolike_core/theory/covariances_fourier.c"
#include "../cosmolike_core/theory/covariances_real_bin_fft.c"
#include "../cosmolike_core/theory/run_covariances_real_bin_fft.c"
#include "init.c"

#include "../cosmolike_core/2dfftlog/utils_complex.h"

#include <time.h>

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

  if (argc != 3){
    fprintf(stderr, "Syntax: %s  block_number  config_file\n", argv[0]);
    exit(1);
  }

  int hit=atoi(argv[1]);
  char *inifile;
  inifile = argv[2];

  FILE *F1,*F2;
  int i,l,m,n,o,s,p,output;
  double ktmp;
  char OUTFILE[400],filename[400];

  // set this to one to output details about inputs for diagnostics
  output = 0;
  FILE *F;

  set_cosmological_parameters(inifile,1);
  set_survey_parameters(inifile,1);
  set_cov_parameters(inifile,1);

  init_source_sample(redshift.shear_REDSHIFT_FILE,tomo.shear_Nbin);
  init_lens_sample(redshift.clustering_REDSHIFT_FILE,tomo.clustering_Nbin);

  if (covparams.lin_bins){
    printf("covariances_real_flat_fft does not support linear angular binning\nEXIT\n");exit(1);
  }
  if (covparams.full_tomo){
    printf("covariances_real_flat_fft does not support full cross-zbins clustering yet\nEXIT\n");exit(1);
  }

  int NG, cNG;
  if (covparams.ng==1){
    NG = 1;
    if (covparams.cng==1){
      cNG = 1;
    }
    else{
      cNG = 0;
    }
  }
  else {
    NG = 0;
    cNG = 0;
  }

  printf("running multi_covariance_real with NG = %d, cNG = %d\n",NG, cNG);

  like.Ntheta=covparams.ntheta;
  like.vtmin = covparams.tmin;
  like.vtmax = covparams.tmax;
  int Ntheta = like.Ntheta;
  double *thetamin,*dtheta,*theta;
  thetamin=create_double_vector(0,Ntheta);
  dtheta=create_double_vector(0,Ntheta-1);
  theta=create_double_vector(0,Ntheta-1);
  set_angular_binning(thetamin,dtheta);
  for (i =0; i< Ntheta; i++){
    theta[i] = 2./3.*(pow(thetamin[i+1],3.)-pow(thetamin[i],3.))/(pow(thetamin[i+1],2.)-pow(thetamin[i],2.));
  }

  int k=1;
  if (strcmp(covparams.ss,"true")==0)
  {
    sprintf(OUTFILE,"%s_ssss_++_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);

    for (l=0;l<tomo.shear_Npowerspectra; l++){
      for (m=l;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_shear_shear_real_binned(OUTFILE,covparams.outdir,theta,thetamin,dtheta,Ntheta,l,m,1,1,k);
        }
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_ssss_--_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.shear_Npowerspectra; l++){
      for (m=l;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_shear_shear_real_binned(OUTFILE,covparams.outdir,theta,thetamin,dtheta,Ntheta,l,m,0,0,k);
        }
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_ssss_+-_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.shear_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_shear_shear_real_binned(OUTFILE,covparams.outdir,theta,thetamin,dtheta,Ntheta,l,m,1,0,k);
        }
        k=k+1;
      }
    }
  }
  if (strcmp(covparams.ll,"true")==0)
  {
    sprintf(OUTFILE,"%s_llll_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=l;m<tomo.clustering_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_clustering_real_binned(OUTFILE,covparams.outdir,theta,thetamin,dtheta,Ntheta,l,m,k);
        }
        k=k+1;
      }
    }
  }
  if (strcmp(covparams.ls,"true")==0)
  {
    sprintf(OUTFILE,"%s_lsls_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=l;m<tomo.ggl_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_ggl_real_binned(OUTFILE,covparams.outdir,theta,thetamin,dtheta,Ntheta,l,m,k);
        }
        k=k+1;
      }
    }
  }
  if (strcmp(covparams.ls,"true")==0 && strcmp(covparams.ss,"true")==0)
  {
    sprintf(OUTFILE,"%s_lsss_+_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_ggl_shear_real_binned(OUTFILE,covparams.outdir,theta,thetamin,dtheta,Ntheta,l,m,1,k);
        }
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_lsss_-_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_ggl_shear_real_binned(OUTFILE,covparams.outdir,theta,thetamin,dtheta,Ntheta,l,m,0,k);
        }
        k=k+1;
      }
    }
  }
  if (strcmp(covparams.ll,"true")==0 && strcmp(covparams.ss,"true")==0)
  {
    sprintf(OUTFILE,"%s_llss_+_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_clustering_shear_real_binned(OUTFILE,covparams.outdir,theta,thetamin,dtheta,Ntheta,l,m,1,k);
        }
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_llss_-_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_clustering_shear_real_binned(OUTFILE,covparams.outdir,theta,thetamin,dtheta,Ntheta,l,m,0,k);
        }
        k=k+1;
      }
    }
  }
  if (strcmp(covparams.ll,"true")==0 && strcmp(covparams.ls,"true")==0)
  {
    sprintf(OUTFILE,"%s_llls_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.ggl_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_clustering_ggl_real_binned(OUTFILE,covparams.outdir,theta,thetamin,dtheta,Ntheta,l,m,k);
        }
        k=k+1;
      }
    }
  }

  if (hit==0)
  {
    sprintf(OUTFILE,"%s%s",covparams.outdir,covparams.filename);
    // write_gglensing_zbins(OUTFILE);

    sprintf(OUTFILE,"%s%s.blocks",covparams.outdir,covparams.filename);
    F1 = fopen(OUTFILE,"w");
    fprintf(F1,"%d\n",k-1);
    fclose(F1);
  }

  printf("number of cov blocks for parallelization: %d\n",k-1);
  printf("-----------------\n");
  printf("PROGRAM EXECUTED\n");
  printf("-----------------\n");
  return 0;
}

max_stars_repo_name: joezuntz/CosmoCov

max_issues_repo_name: joezuntz/CosmoCov

max_forks_repo_name: joezuntz/CosmoCov


text:
#ifndef VERBOSITY_INCLUDE
#define VERBOSITY_INCLUDE

#include <iostream>

#include <gsl/string_span>

#include "config/fleetingOptionsInterface.h"
#include "config/variablesMap.h"
#include "core/task.h"

namespace execHelper {
namespace plugins {
using Verbosity = bool;

const gsl::czstring<> VERBOSITY_KEY = "verbose";

/**
 * \brief Extends the functionality to include the _verbose_ config parameter and processes this parameter, using the --verbose flag
 */
struct VerbosityLong {
    /*! @copydoc JobsLong::getVariables(config::VariablesMap&, const config::FleetingOptionsInterface&)
     */
    static void
    getVariables(config::VariablesMap& variables,
                 const config::FleetingOptionsInterface& options) noexcept;

    /*! @copydoc JobsLong::apply(core::Task&, const config::VariablesMap&)
     */
    inline static void apply(core::Task& task,
                             const config::VariablesMap& variables) noexcept {
        if(*(variables.get<Verbosity>(VERBOSITY_KEY))) {
            task.append("--verbose");
        }
    }
};

/**
 * \brief Extends the functionality to include the _verbose_ config parameter and processes this parameter, using the --debug flag
 */
struct VerbosityDebug {
    /*! @copydoc JobsLong::getVariables(config::VariablesMap&, const config::FleetingOptionsInterface&)
     */
    inline static void
    getVariables(config::VariablesMap& variables,
                 const config::FleetingOptionsInterface& options) noexcept {
        VerbosityLong::getVariables(variables, options);
    }

    /*! @copydoc JobsLong::apply(core::Task&, const config::VariablesMap&)
     */
    inline static void apply(core::Task& task,
                             const config::VariablesMap& variables) noexcept {
        if(*(variables.get<Verbosity>(VERBOSITY_KEY))) {
            task.append("--debug");
        }
    }
};

/**
 * \brief Extends the functionality to include the _verbose_ config parameter and processes this parameter, using the -v flag
 */
struct VerbosityShort {
    /*! @copydoc JobsLong::getVariables(config::VariablesMap&, const config::FleetingOptionsInterface&)
     */
    inline static void
    getVariables(config::VariablesMap& variables,
                 const config::FleetingOptionsInterface& options) noexcept {
        VerbosityLong::getVariables(variables, options);
    }

    /*! @copydoc JobsLong::apply(core::Task&, const config::VariablesMap&)
     */
    inline static void apply(core::Task& task,
                             const config::VariablesMap& variables) noexcept {
        if(*(variables.get<Verbosity>(VERBOSITY_KEY))) {
            task.append("-v");
        }
    }
};
} // namespace plugins
} // namespace execHelper

#endif /* VERBOSITY_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


text:
#pragma once

#include <Python.h>  // has to be the first included header

#ifdef HAVE_LAPACK
#include <lapacke.h>
#endif

#include <stdlib.h>
#include <stdio.h>
#include <vector>
#include <string>

using namespace std;

class SvdUtility
{

 public:

  static void getSVD(double input[], int rows, int cols, double leftSingVec[]);

  static void getSVD(double input[], int rows, int cols, double leftSingVec[], double sigma[], double rightSingVecT[]);

  static void getSVD(double _Complex input[], int rows, int cols, double _Complex leftSingVec[], double sigma[], double _Complex rightSingVecT[]);

  static void reconstructSnapshots(int rows, int cols, double leftSingVec[], double sigma[], double rightSingVecT[], double output[]);

  static void printMatrix(std::string name, double input[], int rows, int cols);

  static vector<double> readCSV(string filename);

  static vector<double> readCSV(string filename, int rows);

  static void writeCSV(string filename, vector<double> values, int m, int n);

  static void writeCSV(string filename, double values[], int m, int n, bool columnWise = false);

  static int getCSVRowCount(string filename);

  static int getCSVColumnCount(string filename);
};

max_stars_repo_name: maierbn/opendihu

max_issues_repo_name: maierbn/opendihu

max_forks_repo_name: maierbn/opendihu


text:
#ifndef GSL_WRAPPERS_H
#define GSL_WRAPPERS_H

// #include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <math.h>
#include <assert.h>
#include <time.h>
#include <sys/stat.h>
#include <sys/types.h>

#define outlog(format, args...) \
    fprintf(stderr, format, args); \
    fprintf(stderr, "\n");

double safe_log(double);
double log_sum(double, double);

static inline double vget(const gsl_vector* v, int i)
{ return(gsl_vector_get(v, i)); };

static inline void vset(gsl_vector* v, int i, double x)
{ gsl_vector_set(v, i, x); };

// Increment a vector element by a double.
void vinc(gsl_vector*, int, double);

static inline double mget(const gsl_matrix* m, int i, int j)
{ return(gsl_matrix_get(m, i, j)); };

static inline void mset(gsl_matrix* m, int i, int j, double x)
{ gsl_matrix_set(m, i, j, x); };

void msetcol(gsl_matrix* m, int r, const gsl_vector* val);

// Increment a matrix element by a double.
void minc(gsl_matrix*, int, int, double);
void msetrow(gsl_matrix*, int, const gsl_vector*);

void col_sum(gsl_matrix*, gsl_vector*);

void vct_printf(const gsl_vector* v);
void mtx_printf(const gsl_matrix* m);
void vct_fscanf(const char*, gsl_vector* v);
void mtx_fscanf(const char*, gsl_matrix* m);
void vct_fprintf(const char* filename, gsl_vector* v);
void mtx_fprintf(const char* filename, const gsl_matrix* m);

double log_det(gsl_matrix*);

void matrix_inverse(gsl_matrix*, gsl_matrix*);

void sym_eigen(gsl_matrix*, gsl_vector*, gsl_matrix*);

double sum(const gsl_vector* v);

double norm(gsl_vector * v);

void vct_log(gsl_vector* v);
void vct_exp(gsl_vector* x);

void choose_k_from_n(int k, int n, int* result);

void log_normalize(gsl_vector* x);
void normalize(gsl_vector* x);

void optimize(int dim,
              gsl_vector* x,
              void* params,
              void (*fdf)(const gsl_vector*, void*, double*, gsl_vector*),
              void (*df)(const gsl_vector*, void*, gsl_vector*),
              double (*f)(const gsl_vector*, void*));

void optimize_fdf(int dim,
                  gsl_vector* x,
                  void* params,
                  void (*fdf)(const gsl_vector*, void*, double*, gsl_vector*),
                  void (*df)(const gsl_vector*, void*, gsl_vector*),
                  double (*f)(const gsl_vector*, void*),
                  double* f_val,
                  double* conv_val,
                  int* niter);

void log_write(FILE* f, char* string);
int directory_exist(const char *dname);
void make_directory(char* name);

gsl_rng* new_random_number_generator();

#endif

max_stars_repo_name: iwangjian/dtm-lab

max_issues_repo_name: iwangjian/topic-extractor

max_forks_repo_name: iwangjian/topic-extractor


text:
// Based on code found in https://github.com/deeplearning4j/libnd4j/blob/master/blas/cpu/NativeBlas.cpp

#include <cblas.h>

#ifdef _WIN32
#include <Windows.h>
#else
#include <dlfcn.h>
#endif

static int maxThreads = -1;
static int vendor = 0;

static void blas_set_num_threads(int num) {
    typedef void* (*void_int)(int);
    typedef int* (*int_int)(int);
    typedef int* (*int_int_int)(int, int);

    maxThreads = num;
#ifdef __MKL
    // if we're linked against mkl - just go for it
    MKL_Set_Num_Threads(num);
    MKL_Domain_Set_Num_Threads(num, 0); // MKL_DOMAIN_ALL
    MKL_Domain_Set_Num_Threads(num, 1); // MKL_DOMAIN_BLAS
    MKL_Set_Num_Threads_Local(num);
#elif __OPENBLAS
#ifdef _WIN32
    // for win32 we just check for mkl_rt.dll
    HMODULE handle = LoadLibrary("mkl_rt.dll");
    if (handle != NULL) {
        void_int mkl_global = (void_int) GetProcAddress(handle, "MKL_Set_Num_Threads");
        if (mkl_global != NULL) {
            mkl_global(num);

            vendor = 3;

            int_int_int mkl_domain = (int_int_int) GetProcAddress(handle, "MKL_Domain_Set_Num_Threads");
            if (mkl_domain != NULL) {
                mkl_domain(num, 0); // DOMAIN_ALL
                mkl_domain(num, 1); // DOMAIN_BLAS
            }

            int_int mkl_local = (int_int) GetProcAddress(handle, "MKL_Set_Num_Threads_Local");
            if (mkl_local != NULL) {
                mkl_local(num);
            }
        } else {
            printf("Unable to tune runtime. Please set OMP_NUM_THREADS manually.\n");
        }
        //FreeLibrary(handle);
    } else {
      // OpenBLAS path
      handle = LoadLibrary("libopenblas.dll");
      if (handle != NULL) {
        void_int oblas = (void_int) GetProcAddress(handle, "openblas_set_num_threads");
        if (oblas != NULL) {
            vendor = 2;
            oblas(num);
        } else {
            printf("Unable to tune runtime. Please set OMP_NUM_THREADS manually.\n");
        }
        //FreeLibrary(handle);
      } else {
        printf("Unable to guess runtime. Please set OMP_NUM_THREADS manually.\n");
      }
    }
#else
    // it's possible to have MKL being loaded at runtime
    void *handle = dlopen("libmkl_rt.so", RTLD_NOW|RTLD_GLOBAL);
    if (handle == NULL) {
        handle = dlopen("libmkl_rt.dylib", RTLD_NOW|RTLD_GLOBAL);
    }
    if (handle != NULL) {

        // we call for openblas only if libmkl isn't loaded, and openblas_set_num_threads exists
        void_int mkl_global = (void_int) dlsym(handle, "MKL_Set_Num_Threads");
        if (mkl_global != NULL) {
            // we're running against mkl
            mkl_global((int) num);

            vendor = 3;

            int_int_int mkl_domain = (int_int_int) dlsym(handle, "MKL_Domain_Set_Num_Threads");
            if (mkl_domain != NULL) {
                mkl_domain(num, 0); // DOMAIN_ALL
                mkl_domain(num, 1); // DOMAIN_BLAS
            }

            int_int mkl_local = (int_int) dlsym(handle, "MKL_Set_Num_Threads_Local");
            if (mkl_local != NULL) {
                mkl_local(num);
            }
        } else {
            printf("Unable to tune runtime. Please set OMP_NUM_THREADS manually.\n");
        }
        dlclose(handle);
    } else {
        // we're falling back to bundled OpenBLAS opening libopenblas.so.0
        handle = dlopen("libopenblas.so.0", RTLD_NOW|RTLD_GLOBAL);
        if (handle == NULL) {
            handle = dlopen("libopenblas.so", RTLD_NOW|RTLD_GLOBAL);
        }
        if (handle == NULL) {
            handle = dlopen("libopenblas.dylib", RTLD_NOW|RTLD_GLOBAL);
        }

        if (handle != NULL) {
            void_int oblas = (void_int) dlsym(handle, "openblas_set_num_threads");
            if (oblas != NULL) {
                vendor = 2;
                // we're running against openblas
                oblas((int) num);
            } else {
                printf("Unable to tune runtime. Please set OMP_NUM_THREADS manually.\n");
            }

            dlclose(handle);
        } else printf("Unable to guess runtime. Please set OMP_NUM_THREADS manually.\n");
    }
#endif

#else
    printf("Unable to guess runtime. Please set OMP_NUM_THREADS or equivalent manually.\n");
#endif
    fflush(stdout);
}

static int blas_get_num_threads() {
    return maxThreads;
}

/**
 *  0 - Unknown
 *  1 - cuBLAS
 *  2 - OpenBLAS
 *  3 - MKL
 */
static int blas_get_vendor() {
    return vendor;
}

max_stars_repo_name: MaxKelsen/javacpp-presets

max_issues_repo_name: MaxKelsen/javacpp-presets

max_forks_repo_name: MaxKelsen/javacpp-presets


text:
/*
Project Name : OpenMEEG

© INRIA and ENPC (contributors: Geoffray ADDE, Maureen CLERC, Alexandre
GRAMFORT, Renaud KERIVEN, Jan KYBIC, Perrine LANDREAU, Théodore PAPADOPOULO,
Emmanuel OLIVI
Maureen.Clerc.AT.sophia.inria.fr, keriven.AT.certis.enpc.fr,
kybic.AT.fel.cvut.cz, papadop.AT.sophia.inria.fr)

The OpenMEEG software is a C++ package for solving the forward/inverse
problems of electroencephalography and magnetoencephalography.

This software is governed by the CeCILL-B license under French law and
abiding by the rules of distribution of free software.  You can  use,
modify and/ or redistribute the software under the terms of the CeCILL-B
license as circulated by CEA, CNRS and INRIA at the following URL
"http://www.cecill.info".

As a counterpart to the access to the source code and  rights to copy,
modify and redistribute granted by the license, users are provided only
with a limited warranty  and the software's authors,  the holders of the
economic rights,  and the successive licensors  have only  limited
liability.

In this respect, the user's attention is drawn to the risks associated
with loading,  using,  modifying and/or developing or reproducing the
software by the user in light of its specific status of free software,
that may mean  that it is complicated to manipulate,  and  that  also
therefore means  that it is reserved for developers  and  experienced
professionals having in-depth computer knowledge. Users are therefore
encouraged to load and test the software's suitability as regards their
requirements in conditions enabling the security of their systems and/or
data to be ensured and,  more generally, to use and operate it in the
same conditions as regards security.

The fact that you are presently reading this means that you have had
knowledge of the CeCILL-B license and that you accept its terms.
*/

#pragma once

//  cmake configuration.
#include <OpenMEEGConfigure.h>
#include <DLLDefinesOpenMEEGMaths.h>

#ifdef USE_MATIO
#include "matio.h"
#endif

#if !defined(USE_MKL)
    #if defined (WIN32) || !defined(USE_ATLAS)
        #define FC_GLOBAL(x,X) x ## _
    #else
        #include <FC.h>
    #endif
#endif

//#define inline __forceinline
//#define inline __attribute__((always_inline))
//#define inline __attribute__((weak)) inline

#if WIN32
    #pragma inline_recursion (on)
    #pragma inline_depth (255)
    #pragma warning( disable : 4530)    //MSVC standard library can't be inlined
    #pragma warning( disable : 4996)    //MSVC warning C4996: declared deprecated
#endif

#ifdef USE_ATLAS
#ifdef __APPLE__
    extern "C" {
        #include <cblas.h>
        #include <clapack.h>
    }
    #define BLAS(x,X) cblas_ ## x
    #define LAPACK(x,X) x ## _
    // #define LAPACK(x,X) FC_GLOBAL(x,X)
#else
    extern "C" {
        #include <cblas.h>
        #include <clapack.h>
    }
    #define BLAS(x,X) cblas_ ## x
    #define LAPACK(x,X) clapack_ ## x
#endif
#endif

#ifdef USE_MKL
// Hack to avoid the MKL declarations of Lapack Functions which do not use the power of C++ references
    #define _MKL_LAPACK_H_
    #include <mkl.h>
    #define BLAS(x,X) cblas_ ## x
    #define LAPACK(x,X) x
    #define FC_GLOBAL(x,X) x
#endif

#ifdef USE_ACML
    #include <acml.h>
    //  Those macros are not used yet
    #define BLAS(x,X) x
    #define LAPACK(x,X) FC_GLOBAL(x,X)
    extern "C" void vrda_sin (int n, double *t, double *p);
    extern "C" void vrda_cos (int n, double *t, double *p);
    extern "C" void vrda_exp (int n, double *t, double *p);
    extern "C" void vrda_log (int n, double *t, double *p);
#endif

#if defined(HAVE_BLAS) && !defined(USE_ATLAS) && !defined(USE_ACML)

#ifndef USE_MKL
    #define CblasColMajor
    #define CblasTrans 'T'
    #define CblasNoTrans 'N'
    #define CblasRight 'R'
    #define CblasLeft 'L'
    #define CblasUpper 'U'
    #define BLAS(x,X) FC_GLOBAL(x,X)
    #define LAPACK(x,X) FC_GLOBAL(x,X)

    extern "C" {
        void BLAS(dcopy,DCOPY)(const int&,const double*,const int&,double*,const int&);
        void BLAS(daxpy,DAXPY)(const int&,const double&,const double*,const int&,double*,const int&);
        double BLAS(ddot,DDOT)(const int&,const double*,const int&,const double*,const int&);
        double BLAS(dnrm2,DNRM2)(const int&,const double*,const int&);
        void BLAS(dscal,DSCAL)(const int&,const double&,double*,const int&);
        void BLAS(dger,DGER)(const int&,const int&,const double&,const double*,const int&,const double*,const int&,double*,const int&);
        void BLAS(dspmv,DSPMV)(const char&,const int&,const double&,const double*,const double*,const int&,const double&,double*,const int&);
        void BLAS(dtpmv,DTPMV)(const char&,const char&,const char&,const int&,const double*,double*,const int&);
        void BLAS(dsymm,DSYMM)(const char&,const char&,const int&,const int&,const double&,const double*,const int&,const double*,const int&, const double&,double*,const int&);
        void BLAS(dgemm,DGEMM)(const char&,const char&,const int&,const int&,const int&,const double&,const double*,const int&,const double*,const int&,const double&,double*,const int&);
        void BLAS(dtrmm,DTRMM)(const char&,const char&,const char&,const char&,const int&,const int&,const double&,const double*,const int&,const double*,const int&);
        void BLAS(dgemv,DGEMV)(const char&,const int&,const int&,const double&,const double*,const int&,const double*,const int&,const double&,double*,const int&);
    }
#endif

    extern "C" {
        void LAPACK(dgetrf,DGETRF)(const int&,const int&,double*,const int&,int*,int&);
        void LAPACK(dgetri,DGETRI)(const int&,double*,const int&,int*,double*,const int&,int&);
    }

#endif

#if defined(HAVE_LAPACK)
    extern "C" {
        void FC_GLOBAL(dgesdd,DGESDD)(const char&,const int&,const int&,double*,const int&,double*,double*,const int&,double*,const int&,double*,const int&,int*,int&);
        void FC_GLOBAL(dpotf2,DPOTF2)(const char&,const int&,double*,const int&,int&);
        double FC_GLOBAL(dlange,DLANGE)(const char&,const int&,const int&,const double*,const int&,double*);
        void FC_GLOBAL(dsptrf,DSPTRF)(const char&,const int&,double*,int*,int&);
        void FC_GLOBAL(dtptri,DTPTRI)(const char&,const char&,const int&,double*,int&,int&);
        void FC_GLOBAL(dsptri,DSPTRI)(const char&,const int&,double*,int*,double*,int&);
        void FC_GLOBAL(dpptrf,DPPTRF)(const char&,const int&,double*,int&);
        void FC_GLOBAL(dpptri,DPPTRI)(const char&,const int&,double*,int&);
        void FC_GLOBAL(dspevd,DSPEVD)(const char&,const char&,const int&,double*,double*,double*,const int&,double*,const int&,int*,const int&,int&);
        void FC_GLOBAL(dsptrs,DSPTRS)(const char&,const int&,const int&,double*,int*,double*,const int&,int&);
    }
#endif

#define DGESDD FC_GLOBAL(dgesdd,DGESDD)
#define DPOTF2 FC_GLOBAL(dpotf2,DPOTF2)
#define DLANGE FC_GLOBAL(dlange,DLANGE)

#define DSPTRF FC_GLOBAL(dsptrf,DSPTRF)
#define DTPTRI FC_GLOBAL(dtptri,DTPTRI)
#define DPPTRF FC_GLOBAL(dpptrf,DPPTRF)
#define DPPTRI FC_GLOBAL(dpptri,DPPTRI)
#define DSPEVD FC_GLOBAL(dspevd,DSPEVD)
#define DSPTRS FC_GLOBAL(dsptrs,DSPTRS)

#if defined(USE_ATLAS) || defined(USE_MKL)
    #define DGER(X1,X2,X3,X4,X5,X6,X7,X8,X9) BLAS(dger,DGER)(CblasColMajor,X1,X2,X3,X4,X5,X6,X7,X8,X9)
    #define DSPMV(X1,X2,X3,X4,X5,X6,X7,X8,X9) BLAS(dspmv,DSPMV)(CblasColMajor,X1,X2,X3,X4,X5,X6,X7,X8,X9)
    #define DTPMV(X1,X2,X3,X4,X5,X6,X7) BLAS(dtpmv,DTPMV)(CblasColMajor,X1,X2,X3,X4,X5,X6,X7)
    #define DSYMM(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12) BLAS(dsymm,DSYMM)(CblasColMajor,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12)
    #define DGEMV(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11) BLAS(dgemv,DGEMV)(CblasColMajor,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11)
    #define DGEMM(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13) BLAS(dgemm,DGEMM)(CblasColMajor,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13)
    #define DTRMM(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11) BLAS(dtrmm,DTRMM)(CblasColMajor,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11)
    #if defined(USE_ATLAS)
        #ifdef __APPLE__
            #define DGETRF(X1,X2,X3,X4,X5,X6) LAPACK(dgetrf,DGETRF)(&X1,&X2,X3,&X4,X5,&X6)
            // #define DGETRF(X1,X2,X3,X4,X5,X6) LAPACK(dgetrf,DGETRF)(CblasColMajor,X1,X2,X3,X4,X5)
            #define DGETRI(X1,X2,X3,X4,X5,X6,X7) LAPACK(dgetri,DGETRI)(&X1,X2,&X3,X4,X5,&X6,&X7)
            // #define DGETRI(X1,X2,X3,X4,X5,X6,X7) LAPACK(dgetri,DGETRI)(CblasColMajor,X1,X2,X3,X4)
        #else
            #define DGETRF(X1,X2,X3,X4,X5) LAPACK(dgetrf,DGETRF)(CblasColMajor,X1,X2,X3,X4,X5)
            #define DGETRI(X1,X2,X3,X4)    LAPACK(dgetri,DGETRI)(CblasColMajor,X1,X2,X3,X4)
        #endif
    #else
        #define DGETRF(X1,X2,X3,X4,X5,X6) LAPACK(dgetrf,DGETRF)(X1,X2,X3,X4,X5,X6)
        #define DGETRI(X1,X2,X3,X4,X5,X6,X7) LAPACK(dgetri,DGETRI)(X1,X2,X3,X4,X5,X6,X7)
    #endif
    #define DSPTRI(X1,X2,X3,X4,X5,X6) FC_GLOBAL(dsptri,DSPTRI)(X1,X2,X3,X4,X5,X6)
#else
    #define DGER(X1,X2,X3,X4,X5,X6,X7,X8,X9) BLAS(dger,DGER)(X1,X2,X3,X4,X5,X6,X7,X8,X9)
    #define DSPMV(X1,X2,X3,X4,X5,X6,X7,X8,X9) BLAS(dspmv,DSPMV)(X1,X2,X3,X4,X5,X6,X7,X8,X9)
    #define DTPMV(X1,X2,X3,X4,X5,X6,X7) BLAS(dtpmv,DTPMV)(X1,X2,X3,X4,X5,X6,X7)
    #define DSYMM(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12) BLAS(dsymm,DSYMM)(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12)
    #define DGEMV(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11) BLAS(dgemv,DGEMV)(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11)
    #define DGEMM(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13) BLAS(dgemm,DGEMM)(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13)
    #define DTRMM(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11) BLAS(dtrmm,DTRMM)(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11)
    #define DGETRF LAPACK(dgetrf,DGETRF)
    #if defined(USE_ACML)
        #define DGETRI(X1,X2,X3,X4,X5,X6,X7) LAPACK(dgetri,DGETRI)(X1,X2,X3,X4,X7)
        #define DSPTRI(X1,X2,X3,X4,X5,X6) LAPACK(dsptri,DSPTRI)(X1,X2,X3,X4,X6)
    #else
        #define DGETRI(X1,X2,X3,X4,X5,X6,X7) LAPACK(dgetri,DGETRI)(X1,X2,X3,X4,X5,X6,X7)
        #define DSPTRI(X1,X2,X3,X4,X5,X6) LAPACK(dsptri,DSPTRI)(X1,X2,X3,X4,X5,X6)
    #endif
#endif

max_stars_repo_name: neurodebian/openmeeg

max_stars_repo_licenses: Intel CECILL-B

max_issues_repo_name: neurodebian/openmeeg

max_issues_repo_licenses: Intel CECILL-B

max_forks_repo_name: neurodebian/openmeeg

max_forks_repo_licenses: Intel CECILL-B


text:
#ifndef __TEST_MAIN_H__
#define __TEST_MAIN_H__

#define PY_SSIZE_T_CLEAN
#include <Python.h>
#include <stdlib.h>
#include <math.h>
#include <gsl_vector.h>
#include "common.h"
#include "svm_main.h"
#include "kernels.h"

// returns predictions and true values
PyObject* __test_t(PyObject* training, PyObject* testing, PyObject* alphas, double b, int rbf);

#endif  /* __TEST_MAIN_H__ */


text:
#ifndef INCLUDED_BLAS_TRAITS_H
#define INCLUDED_BLAS_TRAITS_H

#include <cblas.h>
#include <type_traits>

namespace blas{

    /*
    Specialize to define the blas traits of an object. Should
    be constructed from an instance of T
    */
    template< typename T >
    struct Traits {

    };

    template< typename T >
    constexpr auto blas_cast( T&& x ){
        return Traits<T>( std::forward<T>( x ) );
    }
    template< typename T >
    constexpr CBLAS_ORDER layout( T&& traits ){
        return (CBLAS_ORDER) traits.layout(); }
    template< typename T >
    constexpr CBLAS_TRANSPOSE trans( T&& traits ){
        return (CBLAS_TRANSPOSE) traits.trans(); }
    template< typename T >
    constexpr CBLAS_UPLO uplo( T&& traits ){
        return (CBLAS_UPLO) traits.uplo(); }
    template< typename T >
    constexpr auto nrows( T&& traits ){ return traits.nrows(); }
    template< typename T >
    constexpr auto ncols( T&& traits ){ return traits.ncols(); }
    template< typename T >
    constexpr auto array( T&& traits ){ return traits.array(); }
    template< typename T >
    constexpr auto ld( T&& traits ){ return traits.ld(); }

    template< typename T >
    struct MatrixTraits {
        /*
        using scalar_t = ...
        static CBLAS_ORDER layout( const T& );
        static CBLAS_TRANSPOSE trans( const T& );
        static CBLAS_UPLO uplo( const T& ); --if symmetric
        static std::size_t nrows( const T& );
        static std::size_t ncols( const T& );
        static scalar_t* array( T& );
        static const scalar_t* array( const T& );
        static std::size_t ld( const T& );
        */
    };

    template< typename T >
    struct VectorTraits{
        /*
        using scalar_t = ...
        static std::size_t length( const T& );
        static std::size_t inc( const T& );
        static scalar_t* array( T& );
        static const scalar_t* array( const T& );
        */
    };

    template< typename T >
    struct BlasCast : public std::false_type{
        /*
        When implemented - extend std::true_type (define ::value = true)
        static _ blas_cast( T& );
        static _ blas_cast( const T& );
        */
    };

} //namespace blas

#endif

max_stars_repo_name: tbepler/telepath

max_issues_repo_name: tbepler/telepath

max_forks_repo_name: tbepler/telepath