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

///
/// @file bfgsl.h
///

#ifndef BFGSL_H
#define BFGSL_H

#if defined(__clang__) || defined(__GNUC__)
#pragma GCC system_header
#endif

#include <string.h>

/// @cond

#define concat1(a,b) a ## b
#define concat2(a,b) concat1(a,b)
#define ___ concat2(dont_care, __COUNTER__)

/// @endcond

#ifndef NEED_GSL_LITE

#include <gsl/gsl>

namespace gsl
{

/// @cond

#define expects(cond) Expects(cond)
#define ensures(cond) Ensures(cond)

template <class F>
class final_act_success
{
public:
    explicit final_act_success(F f) noexcept : f_(std::move(f)), invoke_(true) {}

    final_act_success(final_act_success &&other) noexcept : f_(std::move(other.f_)), invoke_(other.invoke_)
    {
        other.invoke_ = false;
    }

    final_act_success(const final_act_success &) = delete;
    final_act_success &operator=(const final_act_success &) = delete;

    ~final_act_success() noexcept
    {
        if (std::uncaught_exception()) {
            return;
        }

        if (invoke_) { f_(); }
    }

private:
    F f_;
    bool invoke_;
};

template <class F>
class final_act_failure
{
public:
    explicit final_act_failure(F f) noexcept : f_(std::move(f)), invoke_(true) {}

    final_act_failure(final_act_failure &&other) noexcept : f_(std::move(other.f_)), invoke_(other.invoke_)
    {
        other.invoke_ = false;
    }

    final_act_failure(const final_act_failure &) = delete;
    final_act_failure &operator=(const final_act_failure &) = delete;

    ~final_act_failure() noexcept
    {
        if (!std::uncaught_exception()) {
            return;
        }

        if (invoke_) { f_(); }
    }

private:
    F f_;
    bool invoke_;
};

template <class F>
inline final_act_success<F> on_success(const F &f) noexcept
{
    return final_act_success<F>(f);
}

template <class F>
inline final_act_success<F> on_success(F &&f) noexcept
{
    return final_act_success<F>(std::forward<F>(f));
}

template <class F>
inline final_act_failure<F> on_failure(const F &f) noexcept
{
    return final_act_failure<F>(f);
}

template <class F>
inline final_act_failure<F> on_failure(F &&f) noexcept
{
    return final_act_failure<F>(std::forward<F>(f));
}

/// @endcond

/// Memset
///
/// Same as std::memset, but for spans
///
/// @param dst The span to memset
/// @param val The value to set the span to
/// @return Returns dst
///
template<class DstElementType, std::ptrdiff_t DstExtent, class T>
auto memset(span<DstElementType, DstExtent> dst, T val)
{
    expects(dst.size() > 0);

    return std::memset(
               dst.data(),
               static_cast<int>(val),
               static_cast<std::size_t>(dst.size())
           );
}

}

#else

#ifdef NEED_STD_LITE
#include <bfstd.h>
#endif

/// @cond

#if defined(__clang__) || defined(__GNUC__)
#define gsl_likely(x) __builtin_expect(!!(x), 1)
#define gsl_unlikely(x) __builtin_expect(!!(x), 0)
#else
#define gsl_likely(x) (x)
#define gsl_unlikely(x) (x)
#endif

#ifndef GSL_ABORT
#define GSL_ABORT abort
#endif

#define expects(cond)                                                                              \
    if (gsl_unlikely(!(cond))) {                                                                   \
        GSL_ABORT();                                                                               \
    }
#define ensures(cond)                                                                              \
    if (gsl_unlikely(!(cond))) {                                                                   \
        GSL_ABORT();                                                                               \
    }

/// @endcond

namespace gsl
{

/// Narrow Cast
///
/// A rename of static_cast to indicate a narrow (e.g. 64bit to 32bit)
///
/// @param u the value to narrow
/// @return static_cast<T>(u)
///
template<class T, class U>
inline constexpr T
narrow_cast(U &&u) noexcept
{
    return static_cast<T>(std::forward<U>(u));
}

/// At
///
/// Returns a reference to an element in an array given an index. Unlike
/// the [] operator, if the indexis out-of-bounds, an exception is thrown,
/// or std::terminate() is called.
///
/// @param arr the array
/// @param index the index of the element to retrieve.
/// @return a reference to
///
template<class T, size_t N, class I>
constexpr T &
at(T(&arr)[N], I index)
{
    expects(index >= 0 && index < narrow_cast<I>(N));
    return arr[static_cast<size_t>(index)];
}

/// At
///
/// Returns a reference to an element in an array given an index. Unlike
/// the [] operator, if the indexis out-of-bounds, an exception is thrown,
/// or std::terminate() is called.
///
/// @param arr the array
/// @param index the index of the element to retrieve.
/// @return a reference to
///
template<class T, size_t N, class I>
const constexpr T &
at(const T(&arr)[N], I index)
{
    expects(index >= 0 && index < narrow_cast<I>(N));
    return arr[static_cast<size_t>(index)];
}

/// At
///
/// Returns a reference to an element in an array given an index. Unlike
/// the [] operator, if the indexis out-of-bounds, an exception is thrown,
/// or std::terminate() is called.
///
/// @param arr the array
/// @param N the size of the array
/// @param index the index of the element to retrieve.
/// @return a reference to
///
template<class T, class I>
constexpr T &
at(T *arr, size_t N, I index)
{
    expects(index >= 0 && index < narrow_cast<I>(N));
    return arr[static_cast<size_t>(index)];
}

/// At
///
/// Returns a reference to an element in an array given an index. Unlike
/// the [] operator, if the indexis out-of-bounds, an exception is thrown,
/// or std::terminate() is called.
///
/// @param arr the array
/// @param N the size of the array
/// @param index the index of the element to retrieve.
/// @return a reference to
///
template<class T, class I>
const constexpr T &
at(const T *arr, size_t N, I index)
{
    expects(index >= 0 && index < narrow_cast<I>(N));
    return arr[static_cast<size_t>(index)];
}

}

#endif

#endif

max_stars_repo_name: chp-io/hypervisor

max_issues_repo_name: chp-io/hypervisor

max_forks_repo_name: chp-io/hypervisor


text:
#pragma once

#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include "core/common.h"
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

void batch_linterp_continuous_axis(float** yy, double* x, size_t size, size_t stride, size_t count);

void linterp_continuous_axis(float** yy, double* x, size_t size, double dt);

void gsl_interp_continuous_axis(float** yy, double* x, size_t size, double dx, gsl_interp_accel **acc, gsl_spline **spline);

void batch_gsl_interp_continuous_axis(float** yy, double* x, size_t size, size_t stride, size_t count);

ext: h

max_stars_repo_name: godsic/semargl-ng

max_issues_repo_name: godsic/semargl-ng

max_forks_repo_name: godsic/semargl-ng


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

#ifndef __GSL_MULTIROOTS_H__
#define __GSL_MULTIROOTS_H__

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

/* Definition of vector-valued functions with parameters based on gsl_vector */

struct gsl_multiroot_function_struct
{
  int (* f) (const gsl_vector * x, void * params, gsl_vector * f);
  size_t n;
  void * params;
};

typedef struct gsl_multiroot_function_struct gsl_multiroot_function ;

#define GSL_MULTIROOT_FN_EVAL(F,x,y) (*((F)->f))(x,(F)->params,(y))

int gsl_multiroot_fdjacobian (gsl_multiroot_function * F,
                              const gsl_vector * x, const gsl_vector * f,
                              double epsrel, gsl_matrix * jacobian);

typedef struct
  {
    const char *name;
    size_t size;
    int (*alloc) (void *state, size_t n);
    int (*set) (void *state, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
    int (*iterate) (void *state, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
    void (*free) (void *state);
  }
gsl_multiroot_fsolver_type;

typedef struct
  {
    const gsl_multiroot_fsolver_type * type;
    gsl_multiroot_function * function ;
    gsl_vector * x ;
    gsl_vector * f ;
    gsl_vector * dx ;
    void *state;
  }
gsl_multiroot_fsolver;

gsl_multiroot_fsolver *
gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * T,
                                     size_t n);

void gsl_multiroot_fsolver_free (gsl_multiroot_fsolver * s);

int gsl_multiroot_fsolver_set (gsl_multiroot_fsolver * s,
                               gsl_multiroot_function * f, gsl_vector * x);

int gsl_multiroot_fsolver_iterate (gsl_multiroot_fsolver * s);

const char * gsl_multiroot_fsolver_name (const gsl_multiroot_fsolver * s);
gsl_vector * gsl_multiroot_fsolver_root (const gsl_multiroot_fsolver * s);
gsl_vector * gsl_multiroot_fsolver_dx (const gsl_multiroot_fsolver * s);
gsl_vector * gsl_multiroot_fsolver_f (const gsl_multiroot_fsolver * s);

/* Definition of vector-valued functions and gradient with parameters
   based on gsl_vector */

struct gsl_multiroot_function_fdf_struct
{
  int (* f) (const gsl_vector * x, void * params, gsl_vector * f);
  int (* df) (const gsl_vector * x, void * params, gsl_matrix * df);
  int (* fdf) (const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix *df);
  size_t n;
  void * params;
};

typedef struct gsl_multiroot_function_fdf_struct gsl_multiroot_function_fdf ;

#define GSL_MULTIROOT_FN_EVAL_F(F,x,y) ((*((F)->f))(x,(F)->params,(y)))
#define GSL_MULTIROOT_FN_EVAL_DF(F,x,dy) ((*((F)->df))(x,(F)->params,(dy)))
#define GSL_MULTIROOT_FN_EVAL_F_DF(F,x,y,dy) ((*((F)->fdf))(x,(F)->params,(y),(dy)))

typedef struct
  {
    const char *name;
    size_t size;
    int (*alloc) (void *state, size_t n);
    int (*set) (void *state, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
    int (*iterate) (void *state, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
    void (*free) (void *state);
  }
gsl_multiroot_fdfsolver_type;

typedef struct
  {
    const gsl_multiroot_fdfsolver_type * type;
    gsl_multiroot_function_fdf * fdf ;
    gsl_vector * x;
    gsl_vector * f;
    gsl_matrix * J;
    gsl_vector * dx;
    void *state;
  }
gsl_multiroot_fdfsolver;

gsl_multiroot_fdfsolver *
gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * T,
                                      size_t n);

int
gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver * s,
                             gsl_multiroot_function_fdf * fdf,
                             gsl_vector * x);

int
gsl_multiroot_fdfsolver_iterate (gsl_multiroot_fdfsolver * s);

void
gsl_multiroot_fdfsolver_free (gsl_multiroot_fdfsolver * s);

const char * gsl_multiroot_fdfsolver_name (const gsl_multiroot_fdfsolver * s);
gsl_vector * gsl_multiroot_fdfsolver_root (const gsl_multiroot_fdfsolver * s);
gsl_vector * gsl_multiroot_fdfsolver_dx (const gsl_multiroot_fdfsolver * s);
gsl_vector * gsl_multiroot_fdfsolver_f (const gsl_multiroot_fdfsolver * s);

int gsl_multiroot_test_delta (const gsl_vector * dx, const gsl_vector * x,
                              double epsabs, double epsrel);

int gsl_multiroot_test_residual (const gsl_vector * f, double epsabs);

GSL_VAR const gsl_multiroot_fsolver_type * gsl_multiroot_fsolver_dnewton;
GSL_VAR const gsl_multiroot_fsolver_type * gsl_multiroot_fsolver_broyden;
GSL_VAR const gsl_multiroot_fsolver_type * gsl_multiroot_fsolver_hybrid;
GSL_VAR const gsl_multiroot_fsolver_type * gsl_multiroot_fsolver_hybrids;

GSL_VAR const gsl_multiroot_fdfsolver_type * gsl_multiroot_fdfsolver_newton;
GSL_VAR const gsl_multiroot_fdfsolver_type * gsl_multiroot_fdfsolver_gnewton;
GSL_VAR const gsl_multiroot_fdfsolver_type * gsl_multiroot_fdfsolver_hybridj;
GSL_VAR const gsl_multiroot_fdfsolver_type * gsl_multiroot_fdfsolver_hybridsj;

__END_DECLS

#endif /* __GSL_MULTIROOTS_H__ */

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


text:
#pragma once
/*
 * (C) Copyright 2020 UCAR
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 */
/// \file Copying.h
/// \brief Generic copying facility
/// \note Feature is under development. This is a placeholder header file.

#include <algorithm>
#include <gsl/gsl-lite.hpp>
#include <memory>
#include <set>
#include <utility>
#include <vector>

#include "ioda/Group.h"
#include "ioda/defs.h"

namespace ioda {
class Group;
class Has_Variables;

class ObjectSelection;
struct ScaleMapping;

/// \todo This function is not yet complete.
///   It exists in its present form to provide the guts for
///   a timing test. Do not use.
/// \brief Generic data copying function
/// \param from contains the objects that we are copying
/// \param to contains the destination(s) of the copy
/// \param scale_map contains settings regarding how dimension scales
///   are propagated
IODA_DL void copy(const ObjectSelection& from, ObjectSelection& to, const ScaleMapping& scale_map);

/// \brief Allows you to select objects for a copy operation
class IODA_DL ObjectSelection {
  // friend IODA_DL void copy(const ObjectSelection&, ObjectSelection&, const ScaleMapping&);
public:
  Group g_;
  bool recurse_ = false;

public:
  ~ObjectSelection();
  ObjectSelection();
  ObjectSelection(const Group& g, bool recurse = true);
  /*
  ObjectSelection(const Variable& v) : ObjectSelection() { insert(v); }
  ObjectSelection(const std::vector<Variable>& v) : ObjectSelection() { insert(v); }
  ObjectSelection(const Has_Variables& v) : ObjectSelection() { insert(v); }
  ObjectSelection(const Group& g, const std::vector<std::string>& v) : ObjectSelection() {
  insert(g,v); }

  void insert(const ObjectSelection&);
  void insert(const Variable&);
  void insert(const std::vector<Variable>&);
  void insert(const Has_Variables&);
  void insert(const Group&, const std::vector<std::string>&);
  void insert(const Group&, bool recurse = true);

  ObjectSelection operator+(const ObjectSelection&) const;
  ObjectSelection operator+(const Variable&) const;
  ObjectSelection operator+(const std::vector<Variable>&) const;
  ObjectSelection operator+(const Has_Variables&) const;
  // Group insertions need to be wrapped by an ObjectSelection.

  ObjectSelection& operator+=(const ObjectSelection&);
  ObjectSelection& operator+=(const Variable&);
  ObjectSelection& operator+=(const std::vector<Variable>&);
  ObjectSelection& operator+=(const Has_Variables&);
  */
};

/// \brief Settings for how to remap dimension scales
struct IODA_DL ScaleMapping {
public:
  ~ScaleMapping();

  std::vector<std::pair<Variable, Variable>> map_from_to;
  std::vector<Variable> map_new;
  bool autocreate = false;
};

}  // namespace ioda

max_stars_repo_name: NOAA-EMC/ioda

max_issues_repo_name: NOAA-EMC/ioda

max_forks_repo_name: NOAA-EMC/ioda


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

#ifndef __GSL_CHECK_RANGE_H__
#define __GSL_CHECK_RANGE_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 <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

GSL_VAR int gsl_check_range;

/* Turn range checking on by default, unless the user defines
   GSL_RANGE_CHECK_OFF, or defines GSL_RANGE_CHECK to 0 explicitly */

#ifdef GSL_RANGE_CHECK_OFF
# ifndef GSL_RANGE_CHECK
#  define GSL_RANGE_CHECK 0
# else
#  error "cannot set both GSL_RANGE_CHECK and GSL_RANGE_CHECK_OFF"
# endif
#else
# ifndef GSL_RANGE_CHECK
#  define GSL_RANGE_CHECK 1
# endif
#endif

__END_DECLS

#endif /* __GSL_CHECK_RANGE_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
/*! \file allvars.c
 *  \brief provides instances of all global variables.
 */

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include "tags.h"
#include "allvars.h"

int ThisTask;		/*!< the rank of the local processor */
int NTask;               /*!< number of processors */
int PTask;	        /*!< smallest integer such that NTask <= 2^PTask */

int NumPart;		/*!< number of particles on the LOCAL processor */
int N_gas;		/*!< number of gas particles on the LOCAL processor  */
long long Ntype[6];      /*!< total number of particles of each type */
int NtypeLocal[6];       /*!< local number of particles of each type */

int NumForceUpdate;      /*!< number of active particles on local processor in current timestep  */
int NumSphUpdate;        /*!< number of active SPH particles on local processor in current timestep  */

double CPUThisRun;	/*!< Sums the CPU time for the process (current submission only) */

int RestartFlag;         /*!< taken from command line used to start code. 0 is normal start-up from
                                     initial conditions, 1 is resuming a run from a set of restart files, while 2
                                     marks a restart from a snapshot file. */

char *Exportflag;        /*!< Buffer used for flagging whether a particle needs to be exported to another process */

int  *Ngblist;           /*!< Buffer to hold indices of neighbours retrieved by the neighbour search routines */

int TreeReconstructFlag; /*!< Signals that a new tree needs to be constructed */

int Flag_FullStep;       /*!< This flag signals that the current step involves all particles */

int ZeroTimestepEncountered;  /*!< Flag used by AMUSE. When a particle is assigned a timestep of zero, an
                                     exception is raised instead of forcing the application to exit. */

gsl_rng *random_generator; /*!< the employed random number generator of the GSL library */

double RndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */

double DomainCorner[3];    /*!< gives the lower left corner of simulation volume */
double DomainCenter[3];    /*!< gives the center of simulation volume */
double DomainLen;          /*!< gives the (maximum) side-length of simulation volume */
double DomainFac;          /*!< factor used for converting particle coordinates to a Peano-Hilbert mesh covering the simulation volume */
int    DomainMyStart;      /*!< first domain mesh cell that resides on the local processor */
int    DomainMyLast;       /*!< last domain mesh cell that resides on the local processor */
int    *DomainStartList;   /*!< a table that lists the first domain mesh cell for all processors */
int    *DomainEndList;     /*!< a table that lists the last domain mesh cell for all processors */
double *DomainWork;        /*!< a table that gives the total "work" due to the particles stored by each processor */
int    *DomainCount;       /*!< a table that gives the total number of particles held by each processor */
int    *DomainCountSph;    /*!< a table that gives the total number of SPH particles held by each processor */

int    *DomainTask;        /*!< this table gives for each leaf of the top-level tree the processor it was assigned to */
int    *DomainNodeIndex;   /*!< this table gives for each leaf of the top-level tree the corresponding node of the gravitational tree */
FLOAT  *DomainTreeNodeLen; /*!< this table gives for each leaf of the top-level tree the side-length of the corresponding node of the gravitational tree */
FLOAT  *DomainHmax;        /*!< this table gives for each leaf of the top-level tree the maximum SPH smoothing length among the particles of the corresponding node of the gravitational tree */

struct DomainNODE
 *DomainMoment;                    /*!< this table stores for each node of the top-level tree corresponding node data from the gravitational tree */

peanokey *DomainKeyBuf;    /*!< this points to a buffer used during the exchange of particle data */

peanokey *Key;             /*!< a table used for storing Peano-Hilbert keys for particles */
peanokey *KeySorted;       /*!< holds a sorted table of Peano-Hilbert keys for all particles, used to construct top-level tree */

int NTopnodes;             /*!< total number of nodes in top-level tree */
int NTopleaves;            /*!< number of leaves in top-level tree. Each leaf can be assigned to a different processor */

struct topnode_data
 *TopNodes;                      /*!< points to the root node of the top-level tree */

double TimeOfLastTreeConstruction; /*!< holds what it says, only used in connection with FORCETEST */

/* variables for input/output, usually only used on process 0 */

char ParameterFile[MAXLEN_FILENAME];  /*!< file name of parameterfile used for starting the simulation */

FILE *FdInfo;       /*!< file handle for info.txt log-file. */
FILE *FdEnergy;     /*!< file handle for energy.txt log-file. */
FILE *FdTimings;    /*!< file handle for timings.txt log-file. */
FILE *FdCPU;        /*!< file handle for cpu.txt log-file. */

#ifdef FORCETEST
FILE *FdForceTest;  /*!< file handle for forcetest.txt log-file. */
#endif

double DriftTable[DRIFT_TABLE_LENGTH];      /*!< table for the cosmological drift factors */
double GravKickTable[DRIFT_TABLE_LENGTH];   /*!< table for the cosmological kick factor for gravitational forces */
double HydroKickTable[DRIFT_TABLE_LENGTH];  /*!< table for the cosmological kick factor for hydrodynmical forces */

void *CommBuffer;   /*!< points to communication buffer, which is used in the domain decomposition, the
                                parallel tree-force computation, the SPH routines, etc. */

/*! This structure contains data which is the SAME for all tasks (mostly code parameters read from the
 * parameter file).  Holding this data in a structure is convenient for writing/reading the restart file, and
 * it allows the introduction of new global variables in a simple way. The only thing to do is to introduce
 * them into this structure.
 */
struct global_data_all_processes
 All;

/*! This structure holds all the information that is
 * stored for each particle of the simulation.
 */
struct particle_data
 *P,              /*!< holds particle data on local processor */
 *DomainPartBuf;  /*!< buffer for particle data used in domain decomposition */

/* the following struture holds data that is stored for each SPH particle in addition to the collisionless
 * variables.
 */
struct sph_particle_data
 *SphP,                        	/*!< holds SPH particle data on local processor */
 *DomainSphBuf;                 /*!< buffer for SPH particle data in domain decomposition */

/*  Variables for Tree
 */

int MaxNodes;		/*!< maximum allowed number of internal nodes */
int Numnodestree;	/*!< number of (internal) nodes in each tree */

struct NODE
 *Nodes_base,                   /*!< points to the actual memory allocted for the nodes */
 *Nodes;                        /*!< this is a pointer used to access the nodes which is shifted such that Nodes[All.MaxPart]
 				     gives the first allocated node */

int *Nextnode;	        /*!< gives next node in tree walk */
int *Father;	        /*!< gives parent node in tree    */

struct extNODE           /*!< this structure holds additional tree-node information which is not needed in the actual gravity computation */
 *Extnodes_base,                /*!< points to the actual memory allocted for the extended node information */
 *Extnodes;                     /*!< provides shifted access to extended node information, parallel to Nodes/Nodes_base */

/*! Header for the standard file format.
 */
struct io_header
 header;  /*!< holds header for snapshot files */

char Tab_IO_Labels[IO_NBLOCKS][4];   /*<! This table holds four-byte character tags used for fileformat 2 */

/* global state of system, used for global statistics
 */
struct state_of_system
 SysState;      /*<! Structure for storing some global statistics about the simulation. */

/* Various structures for communication
 */
struct gravdata_in
 *GravDataIn,                   /*!< holds particle data to be exported to other processors */
 *GravDataGet,                  /*!< holds particle data imported from other processors */
 *GravDataResult,               /*!< holds the partial results computed for imported particles. Note: We use GravDataResult = GravDataGet, such that the result replaces the imported data */
 *GravDataOut;                  /*!< holds partial results received from other processors. This will overwrite the GravDataIn array */

struct gravdata_index
 *GravDataIndexTable;           /*!< the particles to be exported are grouped by task-number. This table allows the results to be disentangled again and to be assigned to the correct particle */

struct densdata_in
 *DensDataIn,                   /*!< holds particle data for SPH density computation to be exported to other processors */
 *DensDataGet;                  /*!< holds imported particle data for SPH density computation */

struct densdata_out
 *DensDataResult,               /*!< stores the locally computed SPH density results for imported particles */
 *DensDataPartialResult;        /*!< imported partial SPH density results from other processors */

struct hydrodata_in
 *HydroDataIn,                  /*!< holds particle data for SPH hydro-force computation to be exported to other processors */
 *HydroDataGet;                 /*!< holds imported particle data for SPH hydro-force computation */

struct hydrodata_out
 *HydroDataResult,              /*!< stores the locally computed SPH hydro results for imported particles */
 *HydroDataPartialResult;       /*!< imported partial SPH hydro-force results from other processors */

#ifdef TIMESTEP_LIMITER
struct timedata_in *TimeDataIn, *TimeDataGet;
#endif

#ifndef NOMPI
#include <mpi.h>
MPI_Comm GADGET_WORLD;
#endif

ext: c

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
/*
 * Copyright (c) 2009-2010 Simon van Heeringen <s.vanheeringen@ncmls.ru.nl>
 *
 * This module is free software. You can redistribute it and/or modify it under
 * the terms of the MIT License, see the file COPYING included with this
 * distribution.
 *
 * This module contains all the code to compare motifs
 *
 */

#include <Python.h>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <gsl/gsl_cdf.h>

void fill_matrix(double matrix[][4], PyObject *matrix_o) {
	// Parse a PyObject matrix into a C array
	// Structure of the matrix is as follows: [[freqA, freqC, freqG, freqT], ...]
	int matrix_len = PyList_Size(matrix_o);
	int i,j;
	PyObject *item;
	for (i = 0; i < matrix_len; i++) {

		PyObject *row = PyList_GetItem(matrix_o, i);
		if (!PyList_Check(row)) {
			PyErr_SetString( PyExc_TypeError, "wrong matrix structure");
			return;
		}
		if (PyList_Size(row) != 4) {
			PyErr_SetString( PyExc_TypeError, "4 nucleotides!");
			return;
		}
		for (j = 0; j < 4; j++) {
			item = PyList_GetItem(row, j);
			matrix[i][j] = PyFloat_AsDouble(item);
		}
	}
}

double mean(double array[], double length) {
	double sum = 0;
	int i;
	for (i = 0; i < length; i++) {
		sum += array[i];
	}
	return sum / length;
}

double sum(double array[], double length) {
	double sum = 0;
	int i;
	for (i = 0; i < length; i++) {
		sum += array[i];
	}
	return sum;
}

double wic(double col1[], double col2[]) {
		int n;
		double x,y;
		double log2 = log(2.0);
		double factor = 2.5;
		double score = 0;;
		double score_a = 0.0;
		double score_b = 0.0;
		double pseudo = 0.0000001;

		for (n = 0; n < 4; n++) {
			x = col1[n] * log((col1[n] + pseudo) / 0.25) / log2;
			y = col2[n] * log((col2[n] + pseudo) / 0.25) / log2;
			score += fabs(x - y);
			score_a += x;
			score_b += y;
		}
		return sqrt(score_a * score_b) - factor * score;
}

double matrix_wic_mean(double matrix1[][4], double matrix2[][4], int length) {
	// Return the mean pcc of two matrices
	int i;
	double result[length];

	for (i = 0; i < length; i++ ) {
		result[i] = wic(matrix1[i], matrix2[i]);
	}
	return mean(result, length);
}

double matrix_wic_sum(double matrix1[][4], double matrix2[][4], int length) {
	// Return the mean pcc of two matrices
	int i;
	double result[length];

	for (i = 0; i < length; i++ ) {
		result[i] = wic(matrix1[i], matrix2[i]);
	}
	return sum(result, length);
}

double pcc(double col1[], double col2[]) {
	// Return the Pearson correlation coefficient of two motif columns
	int n;
	double sum1 = 0;
	double sum2 = 0;
	double mean1, mean2;

	for (n = 0; n < 4; n++) {
		sum1 += col1[n];
		sum2 += col2[n];
	}
	mean1 = sum1 / 4;
	mean2 = sum2 / 4;

	float a = 0;
	float x = 0;
	float y = 0;
	for (n = 0; n < 4; n++) {
		if (((col1[n] - mean1) != 0) && ((col2[n] - mean2) != 0)) {
			a += (col1[n] - mean1) * (col2[n] - mean2);
			x += pow((col1[n] - mean1), 2);
			y += pow((col2[n] - mean2), 2);
		}
		else {
			return 0;
		}
	}

	return  a / sqrt(x * y);
}

double matrix_pcc_mean(double matrix1[][4], double matrix2[][4], int length) {
	// Return the mean pcc of two matrices
	int i;
	double result[length];

	for (i = 0; i < length; i++ ) {
		result[i] = pcc(matrix1[i], matrix2[i]);
	}
	return mean(result, length);
}

double ed(double col1[], double col2[]) {
	// Return the euclidian distance of two motif columns
	int n;
	double score = 0;

	for (n = 0; n < 4; n++) {
		score = score + pow((col1[n] - col2[n]), 2);
	}

	return  -sqrt(score);
}

double matrix_ed_mean(double matrix1[][4], double matrix2[][4], int length) {
	// Return the mean euclidian distance of two matrices
	int i;
	double result[length];

	for (i = 0; i < length; i++ ) {
		result[i] = ed(matrix1[i], matrix2[i]);
	}
	return mean(result, length);
}

double distance(double col1[], double col2[]) {
	// Return the distance between two motifs (Harbison et al.)
	int n;
	double d = 0;

	for (n = 0; n < 4; n++) {
		d = d + pow(col1[n] - col2[n], 2);
	}
	d = d * 1 / sqrt(2);

	return  -d;
}

double matrix_distance_mean(double matrix1[][4], double matrix2[][4], int length) {
	// Return the mean distance (Harbison et al.) of two matrices
	int i;
	double result[length];

	for (i = 0; i < length; i++ ) {
		result[i] = distance(matrix1[i], matrix2[i]);
	}
	return mean(result, length);
}

double chisq(double col1[], double col2[]) {
	// Return the distance between two motifs (Harbison et al.)
	int n;
	double sum1 = 0;
	double sum2 = 0;
	double e1, e2;
	double c = 0;

	for (n = 0; n < 4; n++) {
		sum1 += col1[n];
		sum2 += col2[n];
	}

	for (n = 0; n < 4; n++) {
		e1 = (sum1 * (col1[n] + col2[n])) / (sum1 + sum2);
		e2 = (sum2 * (col1[n] + col2[n])) / (sum1 + sum2);
		c += pow((col1[n] - e1), 2) / e1;
		c += pow((col2[n] - e2), 2) / e2;

	}

	return 1 - gsl_cdf_chisq_P(c, 3.0);
}

double matrix_chisq_mean(double matrix1[][4], double matrix2[][4], int length) {
	// Return the mean distance (Harbison et al.) of two matrices
	int i;
	double result[length];

	for (i = 0; i < length; i++ ) {
		result[i] = chisq(matrix1[i], matrix2[i]);
	}
	return mean(result, length);
}

static PyObject * c_metrics_score(PyObject *self, PyObject * args)
{
	PyObject *matrix1_o;
	PyObject *matrix2_o;
	const char *metric;
	const char *combine;
	int matrix1_len;
	int matrix2_len;

	if (!PyArg_ParseTuple(args, "OOss", &matrix1_o, &matrix2_o, &metric, &combine)) {
		PyErr_SetString( PyExc_TypeError, "Error parsing arguments");
		return NULL;
	}

	// Retrieve frequency matrix
	if (!PyList_Check(matrix1_o)) {
		PyErr_SetString( PyExc_TypeError, "Error: missing required argument matrix1");
		return NULL;
	}

	if (!PyList_Check(matrix2_o)) {
		PyErr_SetString( PyExc_TypeError, "Error: missing required argument matrix2");
		return NULL;
	}
	// Fill matrix1
	matrix1_len = PyList_Size(matrix1_o);
	double matrix1[matrix1_len][4];
	fill_matrix(matrix1, matrix1_o);

	// Fill matrix2
	matrix2_len = PyList_Size(matrix2_o);
	double matrix2[matrix2_len][4];
	fill_matrix(matrix2, matrix2_o);

	if (matrix1_len != matrix2_len) {
		PyErr_SetString( PyExc_TypeError, "Error: matrices have different sizes");
		return NULL;
	}
	int i;
	double result[matrix1_len];
	if (!strcmp(metric, "chisq")) {
		for (i = 0; i < matrix1_len; i++ ) {
			result[i] = chisq(matrix1[i], matrix2[i]);
		}
	}
	else if (!strcmp(metric, "wic")) {
		for (i = 0; i < matrix1_len; i++ ) {
			result[i] = wic(matrix1[i], matrix2[i]);
		}
	}
	else if (!strcmp(metric, "ed")) {
		for (i = 0; i < matrix1_len; i++ ) {
			result[i] = ed(matrix1[i], matrix2[i]);
		}
	}
	else if (!strcmp(metric, "pcc")) {
		for (i = 0; i < matrix1_len; i++ ) {
			result[i] = pcc(matrix1[i], matrix2[i]);
		}
	}
	else if (!strcmp(metric, "distance")) {
		for (i = 0; i < matrix1_len; i++ ) {
			result[i] = distance(matrix1[i], matrix2[i]);
		}
	}
	else {
			PyErr_SetString( PyExc_TypeError, "Unknown metric");
			return NULL;
	}

	if (!strcmp(combine, "mean")) {
		return Py_BuildValue("f", mean(result, matrix1_len));
	}
	else if (!strcmp(combine, "sum")) {
		return Py_BuildValue("f", sum(result, matrix1_len));
	}
	else {
			PyErr_SetString( PyExc_TypeError, "Unknown combine");
			return NULL;
	}

}

int get_truncate_len(len1, len2, pos) {
	//

	if (pos < 0) {
		len2 += pos;
	}
	else if (pos > 0) {
		len1 -= pos;
	}

	if (len1 > len2) {
		return len2;
	}
	else {
		return len1;
	}

}

void fill_tmp_matrices(double matrix1[][4], double matrix2[][4], int pos, int l, double tmp_matrix1[][4], double tmp_matrix2[][4]) {
	int start1 = 0;
	int start2 = 0;
	int i, n;

	if (pos < 0) {
		start2 = -pos;
	}
	else {
		start1 = pos;
	}

	for (i = 0; i < l; i++) {
		for (n = 0; n < 4; n++) {
			tmp_matrix1[i][n] = matrix1[start1 + i][n];
			tmp_matrix2[i][n] = matrix2[start2 + i][n];
		}
	}

}

int index_at_max(double scores[], int len) {
	double max = scores[0];
	int i_at_max = 0;
	int i;

	if (len == 1) {
		return 0;
	}

	for (i = 1; i < len; i++) {
		if (scores[i] > max) {
			max = scores[i];
			i_at_max = i;
		}
	}
	return i_at_max;
}

void fill_rc_matrix(double matrix[][4], int len, double rc_matrix[][4]) {
	int i = 0;

	for (i = 0; i < len; i++) {
		rc_matrix[i][0] = matrix[len - i - 1][3];
		rc_matrix[i][1] = matrix[len - i - 1][2];
		rc_matrix[i][2] = matrix[len - i - 1][1];
		rc_matrix[i][3] = matrix[len - i - 1][0];
	}

}

static PyObject * c_metrics_max_subtotal(PyObject *self, PyObject * args)
{
	PyObject *matrix1_o;
	PyObject *matrix2_o;
	const char *metric;
	const char *combine;
	int matrix1_len;
	int matrix2_len;
	if (!PyArg_ParseTuple(args, "OOss", &matrix1_o, &matrix2_o, &metric, &combine)) {
		PyErr_SetString( PyExc_TypeError, "Error parsing arguments");
		return NULL;
	}

	// Retrieve frequency matrix
	if (!PyList_Check(matrix1_o)) {
		PyErr_SetString( PyExc_TypeError, "Error: missing required argument matrix1");
		return NULL;
	}

	if (!PyList_Check(matrix2_o)) {
		PyErr_SetString( PyExc_TypeError, "Error: missing required argument matrix2");
		return NULL;
	}
	// Fill matrix1
	matrix1_len = PyList_Size(matrix1_o);
	double matrix1[matrix1_len][4];
	fill_matrix(matrix1, matrix1_o);

	// Fill matrix2
	matrix2_len = PyList_Size(matrix2_o);
	double matrix2[matrix2_len][4];
	fill_matrix(matrix2, matrix2_o);

	// Fill rc matrix2
	double rc_matrix2[matrix2_len][4];
	fill_rc_matrix(matrix2, matrix2_len, rc_matrix2);

	// Assign correct function

	int pos = -2;
	int l;
	float max_score;
	int min_overlap = 6;
	int nr_matches;
	int start_pos = -(matrix2_len - min_overlap);
	nr_matches = matrix1_len + matrix2_len - 2 * min_overlap + 1;
	int i;
	int positions[nr_matches * 2];
	double scores[nr_matches * 2];
	int orients[nr_matches * 2];

	double (*ptr_metric_function)(double[][4], double[][4], int) = NULL;

	if ((!strcmp(metric, "chisq")) && (!strcmp(combine, "mean"))) {
		ptr_metric_function = &matrix_chisq_mean;
	}
	else if ((!strcmp(metric, "wic")) && (!strcmp(combine, "mean"))) {
		ptr_metric_function = &matrix_wic_mean;
	}
	else if ((!strcmp(metric, "wic")) && (!strcmp(combine, "sum"))) {
		ptr_metric_function = &matrix_wic_sum;
	}
	else if ((!strcmp(metric, "pcc")) && (!strcmp(combine, "mean"))) {
		ptr_metric_function = &matrix_pcc_mean;
	}
	else if ((!strcmp(metric, "ed")) && (!strcmp(combine, "mean"))) {
		ptr_metric_function = &matrix_ed_mean;
	}
	else if ((!strcmp(metric, "distance")) && (!strcmp(combine, "mean"))) {
		ptr_metric_function = &matrix_distance_mean;
	}
	else {
		PyErr_SetString( PyExc_TypeError, "Unknown metric or combination");
		return NULL;
	}

	for (i = 0; i < nr_matches; i++) {
		pos = i + start_pos;
		l = get_truncate_len(matrix1_len, matrix2_len, pos);

		// Normal matrix
		double tmp_matrix1[l][4];
		double tmp_matrix2[l][4];

		fill_tmp_matrices(matrix1, matrix2, pos, l, tmp_matrix1, tmp_matrix2);
		max_score = (*ptr_metric_function) (tmp_matrix1, tmp_matrix2, l);
		positions[i] = pos;
		scores[i] = max_score;
		orients[i] = 1;
		//printf("CC 1 Len %i score %f\n", l, max_score);

		// Reverse complement matrix
		fill_tmp_matrices(matrix1, rc_matrix2, pos, l, tmp_matrix1, tmp_matrix2);
		max_score = (*ptr_metric_function) (tmp_matrix1, tmp_matrix2, l);
		positions[i + nr_matches] = pos;
		scores[i + nr_matches] = max_score;
		orients[i+ nr_matches] = -1;

		//printf("CC -1 Len %i score %f\n", l, max_score);
	}

	i = index_at_max(scores, nr_matches * 2);
	return Py_BuildValue("fii", scores[i], positions[i], orients[i]);

}
static PyObject * c_metrics_pwmscan(PyObject *self, PyObject * args)
{
	PyObject *seq_o;
	PyObject *pwm_o;
	PyObject *cutoff_o;
	char *seq;
	int seq_len;
	int n_report;
	int pwm_len;
	int i, j;
	double zeta = 0.01;
	int scan_rc;

	if (!PyArg_ParseTuple(args, "OOOii", &seq_o, &pwm_o, &cutoff_o, &n_report, &scan_rc))
		return NULL;

	// Sequence and length
	if (!PyString_Check(seq_o))
		return NULL;
	seq = PyString_AsString(seq_o);
	seq_len = PyString_Size(seq_o);
	//Py_DECREF(seq_o);

	// Retrieve frequency matrix
	if (!PyList_Check(pwm_o))
		return NULL;

	// Weight matrices
	pwm_len = PyList_Size(pwm_o);
	double pwm[pwm_len][4];
	fill_matrix(pwm, pwm_o);
	//Py_DECREF(pwm_o);

	// Cutoff for every spacer length
	double cutoff;
	cutoff = PyFloat_AsDouble(cutoff_o);

	// Scan sequence
	int j_max = seq_len - pwm_len + 1;
	double score_matrix[j_max];
	double rc_score_matrix[j_max];
	double score, rc_score;

	if (j_max < 0) { j_max = 0;}
	int m;
	double g = 0.25;
	double z = 0.01;
	for (j = 0; j < j_max; j++) {
		score = 0;
		rc_score = 0;
		for (m = 0; m < pwm_len; m++) {
			switch(seq[j + m]) {
				case 'A':
					score += log(pwm[m][0] / g + z);
					rc_score += log(pwm[pwm_len - m - 1][3] / g + z);
					break;
				case 'C':
					score += log(pwm[m][1] / g + z);
					rc_score += log(pwm[pwm_len - m - 1][2] / g + z);
					break;
				case 'G':
					score += log(pwm[m][2] / g + z);
					rc_score += log(pwm[pwm_len - m - 1][1] / g + z);
					break;
				case 'T':
					score += log(pwm[m][3] / g + z);
					rc_score += log(pwm[pwm_len - m - 1][0] / g + z);
					break;
			}

		}
		score_matrix[j] = score;
		rc_score_matrix[j] = rc_score;
	}

	// Initialize matrices of n_report highest scores and corresponding positions + strands
	double maxScores[n_report];
	double maxPos[n_report];
	int maxStrand[n_report];
	for (j = 0; j < n_report; j++) {
		maxScores[j] = -100;
		maxPos[j] = -1;
		maxStrand[j] = 1;
	}
	PyObject*  return_list = PyList_New(0);

	int p,q;
	PyObject *x;
	for (j = 0; j < j_max; j++) {
		score = score_matrix[j];
		if (n_report > 0) {
			if (score >= cutoff) {
				p = n_report - 1;
				while ((p >= 0) && (score > maxScores[p])) {
					p--;
				}
				if (p < (n_report-1)) {
					for (q = n_report - 1; q > (p + 1); q--) {
						maxScores[q] = maxScores[q - 1];
						maxPos[q] = maxPos[q - 1];
						maxStrand[q] = maxStrand[q - 1];
					}
					maxScores[p + 1] = score;
					maxPos[p + 1] = j;
					maxStrand[p + 1] = 1;
				}
			}
		}
		else {
			if (score >= cutoff) {
				PyObject* row = PyList_New(0);
				x = PyFloat_FromDouble(score);
				PyList_Append(row, x); Py_DECREF(x);
				x = PyInt_FromLong((long) j);
				PyList_Append(row, x); Py_DECREF(x);
				x = PyInt_FromLong((long) 1);
				PyList_Append(row, x); Py_DECREF(x);
				PyList_Append(return_list, row);
				Py_DECREF(row);
			}
		}
	}

	if (scan_rc) {
		for (j = 0; j < j_max; j++) {
			score = rc_score_matrix[j];
			if (n_report > 0) {
				if (score >= cutoff) {
					p = n_report - 1;
					while ((p >= 0) && (score > maxScores[p])) {
						p--;
					}
					if (p < (n_report-1)) {
						for (q = n_report - 1; q > (p + 1); q--) {
							maxScores[q] = maxScores[q - 1];
							maxPos[q] = maxPos[q - 1];
							maxStrand[q] = maxStrand[q - 1];
						}
						maxScores[p + 1] = score;
						maxPos[p + 1] = j;
						maxStrand[p + 1] = -1;
					}
				}
			}
			else {
				if (score >= cutoff) {
					PyObject* row = PyList_New(0);
					x = PyFloat_FromDouble(score);
					PyList_Append(row, x); Py_DECREF(x);
					x = PyInt_FromLong((long) j);
					PyList_Append(row, x); Py_DECREF(x);
					x =  PyInt_FromLong((long) 1);
					PyList_Append(row, x); Py_DECREF(x);
					PyList_Append(return_list, row);
					Py_DECREF(row);
				}
			}
		}
	}

	for (i = 0; i < n_report; i++) {
		if (maxPos[i] > - 1) {
			PyObject* row = PyList_New(0);
			x = PyFloat_FromDouble(maxScores[i]);
			PyList_Append(row, x);  Py_DECREF(x);
			x = PyInt_FromLong((long)maxPos[i]);
			PyList_Append(row, x); Py_DECREF(x);
			x = PyInt_FromLong((long)maxStrand[i]);
			PyList_Append(row, x); Py_DECREF(x);
			PyList_Append(return_list, row);
			Py_DECREF(row);
		}
	}

	return return_list;
}

static PyMethodDef CoreMethods[] = {
	{"score", c_metrics_score, METH_VARARGS,"Test"},
	{"c_max_subtotal", c_metrics_max_subtotal, METH_VARARGS,"Test"},
	{"pwmscan", c_metrics_pwmscan, METH_VARARGS,"Test"},
	{NULL, NULL, NULL, 0, NULL}
};

PyMODINIT_FUNC
initc_metrics(void)
{
	(void) Py_InitModule("c_metrics", CoreMethods);
};

max_stars_repo_name: kipkurui/gimmemotifs

max_issues_repo_name: kipkurui/gimmemotifs

max_forks_repo_name: kipkurui/gimmemotifs


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

#include <sys/ieeefp.h>
#include <floatingpoint.h>
#include <signal.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>

int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
  char * out ;

  switch (precision)
    {
    case GSL_IEEE_SINGLE_PRECISION:
      ieee_flags ("set", "precision", "single", out) ;
      break ;
    case GSL_IEEE_DOUBLE_PRECISION:
      ieee_flags ("set", "precision", "double", out) ;
      break ;
    case GSL_IEEE_EXTENDED_PRECISION:
      ieee_flags ("set", "precision", "extended", out) ;
      break ;
    default:
      ieee_flags ("set", "precision", "extended", out) ;
    }

  switch (rounding)
    {
    case GSL_IEEE_ROUND_TO_NEAREST:
      ieee_flags ("set", "direction", "nearest", out) ;
      break ;
    case GSL_IEEE_ROUND_DOWN:
      ieee_flags ("set", "direction", "negative", out) ;
      break ;
    case GSL_IEEE_ROUND_UP:
      ieee_flags ("set", "direction", "positive", out) ;
      break ;
    case GSL_IEEE_ROUND_TO_ZERO:
      ieee_flags ("set", "direction", "tozero", out) ;
      break ;
    default:
      ieee_flags ("set", "direction", "nearest", out) ;
    }

  if (exception_mask & GSL_IEEE_MASK_INVALID)
    {
      ieee_handler ("set", "invalid", SIGFPE_IGNORE) ;
    }
  else
    {
      ieee_handler ("set", "invalid", SIGFPE_ABORT) ;
    }

  if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
    {
      ieee_handler ("set", "denormalized", SIGFPE_IGNORE) ;
    }
  else
    {
      GSL_ERROR ("sunos4 does not support the denormalized operand exception. "
                 "Use 'mask-denormalized' to work around this.",
                 GSL_EUNSUP) ;
    }

  if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
    {
      ieee_handler ("set", "division", SIGFPE_IGNORE) ;
    }
  else
    {
      ieee_handler ("set", "division", SIGFPE_ABORT) ;
    }

  if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
    {
      ieee_handler ("set", "overflow", SIGFPE_IGNORE) ;
    }
  else
    {
      ieee_handler ("set", "overflow", SIGFPE_ABORT) ;
    }

  if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
    {
      ieee_handler ("set", "underflow", SIGFPE_IGNORE) ;
    }
  else
    {
      ieee_handler ("set", "underflow", SIGFPE_ABORT) ;
    }

  if (exception_mask & GSL_IEEE_TRAP_INEXACT)
    {
      ieee_handler ("set", "inexact", SIGFPE_ABORT) ;
    }
  else
    {
      ieee_handler ("set", "inexact", SIGFPE_IGNORE) ;
    }

  return GSL_SUCCESS ;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


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

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_sf.h>

#include "test_sf.h"

#define WKB_TOL (1.0e+04 * TEST_SQRT_TOL0)

/*
 * Test the identities:
 *
 * Sum_{k=0}^n (n choose k) (2 y)^{n - k} H_k(x) = H_n(x + y)
 *
 * and
 *
 * Sum_{k=0}^n (n choose k) y^{n-k} He_k(x) = He_n(x + y)
 *
 * see: http://mathworld.wolfram.com/HermitePolynomial.html (Eq. 55)
 */
void
test_hermite_id1(const int n, const double x, const double y)
{
  double *a = malloc((n + 1) * sizeof(double));
  double *b = malloc((n + 1) * sizeof(double));
  double lhs, rhs;
  int k;

  a[0] = gsl_pow_int(2.0 * y, n);
  b[0] = gsl_pow_int(y, n);
  for (k = 1; k <= n; ++k)
    {
      double fac = (n - k + 1.0) / (k * y);
      a[k] = 0.5 * fac * a[k - 1];
      b[k] = fac * b[k - 1];
    }

  lhs = gsl_sf_hermite_phys_series(n, x, a);
  rhs = gsl_sf_hermite_phys(n, x + y);
  gsl_test_rel(lhs, rhs, TEST_TOL4, "identity1 phys n=%d x=%g y=%g", n, x, y);

  lhs = gsl_sf_hermite_prob_series(n, x, b);
  rhs = gsl_sf_hermite_prob(n, x + y);
  gsl_test_rel(lhs, rhs, TEST_TOL3, "identity1 prob n=%d x=%g y=%g", n, x, y);

  free(a);
  free(b);
}

int
test_hermite(void)
{
  gsl_sf_result r;

  int s = 0;
  int m, n, sa;
  double res[256];
  double x;
  const double aizero1 = -2.3381074104597670384891972524467; /* first zero of the Airy function Ai */

  /* test some known identities */
  test_hermite_id1(10, 0.75, 0.33);
  test_hermite_id1(8, -0.75, 1.20);
  test_hermite_id1(7, 2.88, -3.2);

  x = 0.75;

  TEST_SF(s, gsl_sf_hermite_prob_e, (0, 0.75, &r),  1.,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_prob_e, (1, 0.75, &r),  x,  TEST_TOL0, GSL_SUCCESS);

  n = 25;
  TEST_SF(s, gsl_sf_hermite_prob_e, (n, 0., &r),  0.,  TEST_TOL0, GSL_SUCCESS);
  n = 28;
  TEST_SF(s, gsl_sf_hermite_prob_e, (n, 0., &r),  213458046676875,  TEST_TOL0, GSL_SUCCESS);
  n = 25;
  TEST_SF(s, gsl_sf_hermite_prob_e, (n, 0.75, &r),  -1.08128685847680748265939328423e12,  TEST_TOL0, GSL_SUCCESS);
  n = 28;
  TEST_SF(s, gsl_sf_hermite_prob_e, (n, 0.75, &r),  -1.60620252094658918105511125135e14,  TEST_TOL0, GSL_SUCCESS);

#if 0
  n = 10025;
  x = ((sqrt(2*n+1.)+aizero1/pow(8.*n,1/6.))/2)*M_SQRT2;
  TEST_SF(s, gsl_sf_hermite_prob_e, (n, x, &r),  -3.3090527852387782540121569578e18961,  TEST_TOL0, GSL_SUCCESS);

  n = 10028;
  x = ((sqrt(2*n+1.)+aizero1/pow(8.*n,1/6.))/2)*M_SQRT2;
  TEST_SF(s, gsl_sf_hermite_prob_e, (n, x, &r),  -7.515478445930242044360704363e18967,  TEST_TOL0, GSL_SUCCESS);

  n = 10025;
  x = (sqrt(2*n+1.)-(aizero1/pow(8.*n,1/6.))/2)*M_SQRT2;
  TEST_SF(s, gsl_sf_hermite_prob_e, (n, x, &r),  4.1369269649092456235914193753e22243,  TEST_TOL0, GSL_SUCCESS);

  n = 10028;
  x = (sqrt(2*n+1.)-(aizero1/pow(8.*n,1/6.))/2)*M_SQRT2;
  TEST_SF(s, gsl_sf_hermite_prob_e, (n, x, &r),  8.363694992558646923734666303e22250,  TEST_TOL0, GSL_SUCCESS);

  n = 10025;
  x = (sqrt(2*n+1.)-2*(aizero1/pow(8.*n,1/6.)))*M_SQRT2;
  TEST_SF(s, gsl_sf_hermite_prob_e, (n, x, &r),  7.398863979737363164340057757e22273,  TEST_TOL0, GSL_SUCCESS);

  n = 10028;
  x = (sqrt(2*n+1.)-2*(aizero1/pow(8.*n,1/6.)))*M_SQRT2;
  TEST_SF(s, gsl_sf_hermite_prob_e, (n, x, &r),  1.507131397474022356488976968e22281,  TEST_TOL0, GSL_SUCCESS);
#endif

  x = 0.75;

  n = 128;
  m = 225;
  TEST_SF(s, gsl_sf_hermite_prob_der_e, (m, n, x, &r),  0.,  TEST_TOL0, GSL_SUCCESS);

  n = 128;
  m = 5;
  TEST_SF(s, gsl_sf_hermite_prob_der_e, (m, n, x, &r),  -3.0288278964712702882066404e112,  TEST_TOL1, GSL_SUCCESS);

  sa = 0;
  gsl_sf_hermite_prob_array(0, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_array(0, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array(1, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[1],   +0.0,   x,   TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_array(1, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array(100, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   823.810509681701660156250, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,  1.03749254986255755872498e78, TEST_TOL1);
  gsl_test(sa, "gsl_sf_hermite_prob_array(100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array_der(0, 100, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   823.810509681701660156250, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,  1.03749254986255755872498e78, TEST_TOL1);
  gsl_test(sa, "gsl_sf_hermite_prob_array_der(0, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array_der(1000, 100, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   0.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_array_der(1000, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array_der(23, 100, x, res);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[37], +0.0,  2.3592417210568968566591219172e37, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0, 2.6503570965896336273549197e100, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_array(23, 37, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_der_array(100, 50, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   -3.88338863813139372375411561e31, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   7.9614368698398116765703194e38, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   0.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_der_array(100, 50, 0.75)");
  s += sa;

  n = 128;
  res[0] = 1.;
  for(m=1; m<=n; m++){
    res[m] = res[m-1]/2.;
  }
  TEST_SF(s, gsl_sf_hermite_prob_series_e, (n, x, res, &r),  -4.0451066556993485405907339548e68,  TEST_TOL0, GSL_SUCCESS);

  /* phys */

  x = 0.75;

  TEST_SF(s, gsl_sf_hermite_phys_e, (0, 0.75, &r),  1.,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_phys_e, (1, 0.75, &r),  2.*x,  TEST_TOL0, GSL_SUCCESS);

  n = 25;
  TEST_SF(s, gsl_sf_hermite_phys_e, (n, 0., &r),  0.,  TEST_TOL0, GSL_SUCCESS);
  n = 28;
  TEST_SF(s, gsl_sf_hermite_phys_e, (n, 0., &r),  3497296636753920000,  TEST_TOL0, GSL_SUCCESS);
  n = 25;
  TEST_SF(s, gsl_sf_hermite_phys_e, (n, 0.75, &r),  -9.7029819451106077507781088352e15,  TEST_TOL0, GSL_SUCCESS);
  n = 28;
  TEST_SF(s, gsl_sf_hermite_phys_e, (n, 0.75, &r),  3.7538457078067672096408339776e18,  TEST_TOL0, GSL_SUCCESS);

#if 0
  n = 10025;
  x = ((sqrt(2*n+1.)+aizero1/pow(8.*n,1/6.))/2);
  TEST_SF(s, gsl_sf_hermite_phys_e, (n, x, &r),  -2.7074282783315424535693575770e20470,  TEST_TOL0, GSL_SUCCESS);

  n = 10028;
  x = ((sqrt(2*n+1.)+aizero1/pow(8.*n,1/6.))/2);
  TEST_SF(s, gsl_sf_hermite_phys_e, (n, x, &r),  -1.7392214893577690864150561850e20477,  TEST_TOL0, GSL_SUCCESS);

  n = 10025;
  x = (sqrt(2*n+1.)-(aizero1/pow(8.*n,1/6.))/2);
  TEST_SF(s, gsl_sf_hermite_phys_e, (n, x, &r),  3.3847852473526979744366379542e23752,  TEST_TOL0, GSL_SUCCESS);

  n = 10028;
  x = (sqrt(2*n+1.)-(aizero1/pow(8.*n,1/6.))/2);
  TEST_SF(s, gsl_sf_hermite_phys_e, (n, x, &r),  1.9355145738418079256435712027e23760,  TEST_TOL0, GSL_SUCCESS);

  n = 10025;
  x = (sqrt(2*n+1.)-2*(aizero1/pow(8.*n,1/6.)));
  TEST_SF(s, gsl_sf_hermite_phys_e, (n, x, &r),  6.053663953512337293393128307e23782,  TEST_TOL0, GSL_SUCCESS);

  n = 10028;
  x = (sqrt(2*n+1.)-2*(aizero1/pow(8.*n,1/6.)));
  TEST_SF(s, gsl_sf_hermite_phys_e, (n, x, &r),  3.487782358276961096026268141e23790,  TEST_TOL0, GSL_SUCCESS);
#endif

  x = 0.75;

  n = 128;
  m = 225;
  TEST_SF(s, gsl_sf_hermite_phys_der_e, (m, n, x, &r),  0.,  TEST_TOL0, GSL_SUCCESS);

  n = 128;
  m = 5;
  TEST_SF(s, gsl_sf_hermite_phys_der_e, (m, n, x, &r),  2.89461215568095657569833e132,  TEST_TOL0, GSL_SUCCESS);

  sa = 0;
  gsl_sf_hermite_phys_array(0, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_phys_array(0, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_phys_array(1, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[1],   +0.0,   2.*x,   TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_phys_array(1, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_phys_array(100, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   38740.4384765625, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,  -1.4611185395125104593177790757e93, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_phys_array(100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_phys_array_der(0, 100, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   38740.4384765625, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,  -1.4611185395125104593177790757e93, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_phys_array_der(0, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_phys_array_der(1000, 100, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   0.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_phys_array_der(1000, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_phys_array_der(23, 100, x, res);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[37],  +0.0,   1.930387357696033719818118732e46, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   2.957966000491202678467161e118, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_phys_array(23, 37, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_phys_der_array(100, 50, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   -8.2663221830586310072686183e38, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   1.52281030265187793605875604e49, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   0.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_phys_der_array(100, 50, 0.75)");
  s += sa;

  n = 128;
  /* arbitrary weights */
  res[0] = 1.;
  for(m=1; m<=n; m++){
    res[m] = res[m-1]/2.;
  }
  TEST_SF(s, gsl_sf_hermite_phys_series_e, (n, x, res, &r),  1.07772223811696567390619566842e88,  TEST_TOL0, GSL_SUCCESS);

  x = 0.75;
  n = 28;
  TEST_SF(s, gsl_sf_hermite_func_e, (n, 0, &r),  0.290371943657199641200016132937,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  0.23526280808621240649319140441,  TEST_TOL0, GSL_SUCCESS);
  n = 100028;
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  -0.02903467369856961147236598086,  TEST_TOL4, GSL_SUCCESS);

  n = 10025;
  x = ((sqrt(2*n+1.)+aizero1/pow(8.*n,1/6.))/2.5);
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  -0.05301278920004176459797680403,  TEST_TOL4, GSL_SUCCESS);

  n = 10028;
  x = ((sqrt(2*n+1.)+aizero1/pow(8.*n,1/6.))/2.5);
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  0.06992968509693993526829596970,  TEST_TOL3, GSL_SUCCESS);

  n = 10025;
  x = (sqrt(2*n+1.)-(aizero1/pow(8.*n,1/6.))/2.5);
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  0.08049000991742150521366671021,  TEST_TOL4, GSL_SUCCESS);

  n = 10028;
  x = (sqrt(2*n+1.)-(aizero1/pow(8.*n,1/6.))/2.5);
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  0.08048800667512084252723933250,  TEST_TOL4, GSL_SUCCESS);

  n = 10025;
  x = (sqrt(2*n+1.)-2.5*(aizero1/pow(8.*n,1/6.)));
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  7.97206830806663013555068100e-6,  TEST_TOL4, GSL_SUCCESS);

  n = 10028;
  x = (sqrt(2*n+1.)-2.5*(aizero1/pow(8.*n,1/6.)));
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  7.97188517397786729928465829e-6,  TEST_TOL4, GSL_SUCCESS);

  x = 0.75;

  sa = 0;
  gsl_sf_hermite_func_array(100, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   0.566979307027693616978839335983, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.360329854170806945032958735574, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   -0.07616422890563462489003733382, TEST_TOL1);
  gsl_test(sa, "gsl_sf_hermite_func_array(100, 0.75)");
  s += sa;

  n = 128;
  /* arbitrary weights */
  res[0] = 1.;
  for(m=1; m<=n; m++){
    res[m] = res[m-1]/2.;
  }
  TEST_SF(s, gsl_sf_hermite_func_series_e, (n, x, res, &r),  0.81717103529960997134154552556,  TEST_TOL0, GSL_SUCCESS);

  m = 5;
  n = 28;
  TEST_SF(s, gsl_sf_hermite_func_der_e, (0, n, x, &r),  0.235262808086212406493191404,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_der_e, (m, n, x, &r),  4035.32788513029308540826835,  TEST_TOL1, GSL_SUCCESS);

  {
    /* positive zeros of the probabilists' Hermite polynomial of order 17 */
    double He17z[8] = { 0.751842600703896170737870774614, 1.50988330779674075905491513417, 2.28101944025298889535537879396,
                        3.07379717532819355851658337833, 3.90006571719800990903311840097, 4.77853158962998382710540812497,
                        5.74446007865940618125547815768,6.88912243989533223256205432938 };

    n = 17;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_prob_zero_e, (n, m, &r),  He17z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the probabilists' Hermite polynomial of order 18 */
    double He18z[9] = { 0.365245755507697595916901619097, 1.09839551809150122773848360538, 1.83977992150864548966395498992,
                        2.59583368891124032910545091458, 3.37473653577809099529779309480, 4.18802023162940370448450911428,
                        5.05407268544273984538327527397, 6.00774591135959752029303858752, 7.13946484914647887560975631213 };

    n = 18;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_prob_zero_e, (n, m, &r),  He18z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the probabilists' Hermite polynomial of order 23 */
    double He23z[11] = { 0.648471153534495816722576841197, 1.29987646830397886997876116860, 1.95732755293342410739100839243,
                         2.62432363405918177067330340783, 3.30504002175296456723204112903, 4.00477532173330406712238633738,
                         4.73072419745147329426707133987, 5.49347398647179412855289367747, 6.31034985444839982842886078177,
                         7.21465943505186138595859194492, 8.29338602741735258945596770157 };

    n = 23;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_prob_zero_e, (n, m, &r),  He23z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the probabilists' Hermite polynomial of order 24 */
    double He24z[12] = { 0.317370096629452319318170455994, 0.953421922932109084904629632351, 1.59348042981642010695074168129,
                         2.24046785169175236246653790858, 2.89772864322331368932008199475, 3.56930676407356024709649151613,
                         4.26038360501990548884317727406, 4.97804137463912033462166468006, 5.73274717525120114834341822330,
                         6.54167500509863444148277523364, 7.43789066602166310850331715501, 8.50780351919525720508386233432 };

    n = 24;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_prob_zero_e, (n, m, &r),  He24z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the physicists' Hermite polynomial of order 17 */
    double H17z[8] = { 0.531633001342654731349086553718, 1.06764872574345055363045773799, 1.61292431422123133311288254454,
                       2.17350282666662081927537907149, 2.75776291570388873092640349574, 3.37893209114149408338327069289,
                       4.06194667587547430689245559698, 4.87134519367440308834927655662 };

    n = 17;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_phys_zero_e, (n, m, &r),  H17z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the physicists' Hermite polynomial of order 18 */
    double H18z[9] = { 0.258267750519096759258116098711, 0.776682919267411661316659462284, 1.30092085838961736566626555439,
                       1.83553160426162889225383944409, 2.38629908916668600026459301424, 2.96137750553160684477863254906,
                       3.57376906848626607950067599377, 4.24811787356812646302342016090, 5.04836400887446676837203757885 };

    n = 18;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_phys_zero_e, (n, m, &r),  H18z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the physicists' Hermite polynomial of order 23 */
    double H23z[11] = { 0.458538350068104797757887329284, 0.919151465442563765431719239593, 1.38403958568249523732634717118,
                        1.85567703767137106251504753718, 2.33701621147445578644623502174, 2.83180378712615690144806140734,
                        3.34512715994122457247439814585, 3.88447270810610186607248760288, 4.46209117374000667673186157071,
                        5.10153461047667712968749766165, 5.86430949898457256538748413474 };

    n = 23;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_phys_zero_e, (n, m, &r),  H23z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the physicists' Hermite polynomial of order 24 */
    double H24z[12] = { 0.224414547472515585151136715527, 0.674171107037212236000245923730, 1.12676081761124507213306126773,
                        1.58425001096169414850563336202, 2.04900357366169891178708399532, 2.52388101701142697419907602333,
                        3.01254613756556482565453858421, 3.52000681303452471128987227609, 4.05366440244814950394766297923,
                        4.62566275642378726504864923776, 5.25938292766804436743072304398, 6.01592556142573971734857350899 };

    n = 24;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_phys_zero_e, (n, m, &r),  H24z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  n = 2121;
  x = -1.*n;
  res[0] = (double) n;
  sa = 0;
  for (m=1; m<=n/2; m++) {
    gsl_sf_hermite_prob_zero_e(n, m, &r);
    if (x>=r.val) {
      sa += TEST_SF_INCONS;
      printf("sanity check failed! (gsl_sf_hermite_prob_zero)\n");
    }
    res[0] = GSL_MIN(res[0],fabs(x-r.val));
    x = r.val;
  }
  gsl_test(sa, "gsl_sf_hermite_prob_zero(n, m, r)");

  n = 2121;
  x = -1.*n;
  res[0] = (double) n;
  sa = 0;
  for (m=1; m<=n/2; m++) {
    gsl_sf_hermite_phys_zero_e(n, m, &r);
    if (x>=r.val) {
      sa += TEST_SF_INCONS;
      printf("sanity check failed! (gsl_sf_hermite_phys_zero)\n");
    }
    res[0] = GSL_MIN(res[0],fabs(x-r.val));
    x = r.val;
  }
  gsl_test(sa, "gsl_sf_hermite_phys_zero(n, m, r)");

  return s;
}


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

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multiroots.h>

int
gsl_multiroot_test_delta (const gsl_vector * dx, const gsl_vector * x,
                     double epsabs, double epsrel)
{
  size_t i;
  int ok = 1;
  const size_t n = x->size ;

  if (epsrel < 0.0)
    {
      GSL_ERROR ("relative tolerance is negative", GSL_EBADTOL);
    }

  for (i = 0 ; i < n ; i++)
    {
      double xi = gsl_vector_get(x,i);
      double dxi = gsl_vector_get(dx,i);
      double tolerance = epsabs + epsrel * fabs(xi)  ;

      if (fabs(dxi) < tolerance)
        {
          ok = 1;
        }
      else
        {
          ok = 0;
          break;
        }
    }

  if (ok)
    return GSL_SUCCESS ;

  return GSL_CONTINUE;
}

int
gsl_multiroot_test_residual (const gsl_vector * f, double epsabs)
{
  size_t i;

  double residual = 0;

  const size_t n = f->size;

  if (epsabs < 0.0)
    {
      GSL_ERROR ("absolute tolerance is negative", GSL_EBADTOL);
    }

  for (i = 0 ; i < n ; i++)
    {
      double fi = gsl_vector_get(f, i);

      residual += fabs(fi);
    }

  if (residual < epsabs)
    {
      return GSL_SUCCESS;
    }

  return GSL_CONTINUE ;
}


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
/*

Given a block system
A x = b,
or
( K G )(u) = (f)
( D C )(p)   (h)

We define a symmetrically scaled system,
L A R R^-1 x = L b,		A'x' = b'
where L = diag(L1, 1), R = diag(R1,1) are block diagonal
and L_i and R_i are both diagonal matrirces.

The scaling produces
( L1   )( K G )( R1   )( R1   )^-1(u) = ( L1   )(f)
Du + Cp = (h)
and we solve
( L1   )( K G )( R1   )(v) = ( L1   )(f)
 Du + Cp  = (h)
or
( L1*K*R1, L1*G )(v) = (L1*f)
( D*R1  ,  C    )(p) = ( h  )

The solution u,p is recovered via
(u) = (R1*v)
(p) = (p)

*/

#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>

#include <StGermain/libStGermain/src/StGermain.h>
#include <StgDomain/libStgDomain/src/StgDomain.h>

#include "common-driver-utils.h"
#include "stokes_Kblock_scaling.h"

/* private prototypes */
PetscErrorCode BSSCR_MatStokesKBlock_ApplyScaling( MatStokesBlockScaling BA, Mat A, Vec b, Vec x, Mat S, PetscTruth sym );
PetscErrorCode BSSCR_mat_kblock_invert_scalings( MatStokesBlockScaling BA );

/* A x = b -> A'x' = b' */
/*
Note this routine actually modifies the matrix and rhs b.
*/

// updated
PetscErrorCode BSSCR_MatStokesKBlock_ApplyScaling( MatStokesBlockScaling BA, Mat A, Vec b, Vec x, Mat S, PetscTruth sym )
{
	Mat K,G,D;
	Vec L1,R1;
	Vec f,u;

	/* Get the scalings out the block mat data */
	VecNestGetSubVec( BA->Lz, 0, &L1 );
	VecNestGetSubVec( BA->Rz, 0, &R1 );

	/* get the subblock solution and rhs */
	if( x != PETSC_NULL ) {
		VecNestGetSubVec( x, 0, &u );
		VecPointwiseDivide( u, u,R1); /* x <- x * 1/R1 */
	}
	if( b != PETSC_NULL ) {
		VecNestGetSubVec( b, 0, &f );
		VecPointwiseMult( f, f,L1); /* f <- f * L1 */
	}

	/* Scale matrices */
	MatNestGetSubMat( A, 0,0, &K );
	MatNestGetSubMat( A, 0,1, &G );
	MatNestGetSubMat( A, 1,0, &D );

	if( K != PETSC_NULL ) {		MatDiagonalScale( K, L1,R1 );		}
	if( G != PETSC_NULL ) {		MatDiagonalScale( G, L1,PETSC_NULL );		}
	if( D != PETSC_NULL && !sym ) {	MatDiagonalScale( D, PETSC_NULL,R1 );		}

	PetscFunctionReturn(0);
}

// updated
PetscErrorCode BSSCR_MatStokesKBlockScalingCreate( MatStokesBlockScaling *_BA )
{
	MatStokesBlockScaling BA;
	PetscErrorCode ierr;

	ierr = PetscMalloc( sizeof(struct _p_MatStokesBlockScaling), &BA );CHKERRQ(ierr);

	BA->Lz = PETSC_NULL;
	BA->Rz = PETSC_NULL;

	BA->scaling_exists              = PETSC_FALSE;
	BA->scalings_have_been_inverted = PETSC_FALSE;
	BA->system_has_been_scaled      = PETSC_FALSE;

	*_BA = BA;

	PetscFunctionReturn(0);
}

// updated
PetscErrorCode BSSCR_MatStokesKBlockScalingDestroy( MatStokesBlockScaling BA )
{

	if( BA->scaling_exists == PETSC_FALSE ) PetscFunctionReturn(0);

	if( BA->Lz != PETSC_NULL ) {
		Stg_VecDestroy(&BA->Lz );
		BA->Lz = PETSC_NULL;
	}

	if( BA->Rz != PETSC_NULL ) {
		Stg_VecDestroy(&BA->Rz );
		BA->Rz = PETSC_NULL;
	}

	PetscFree( BA );

	PetscFunctionReturn(0);
}

/*
A is 2x2 block matrix
b and x are 2x1 block vectors
*/
// updated
PetscErrorCode BSSCR_MatKBlock_ConstructScaling( MatStokesBlockScaling BA, Mat A, Vec b, Vec x, PetscTruth sym )
{

	if( BA->scaling_exists == PETSC_FALSE ) {
      //PetscInt M,N;
		//PetscTruth is_block;

		/* check A is 2x2 block matrix */
		/* Stg_PetscTypeCompare( (PetscObject)A, "block", &is_block ); */
		/* if (is_block==PETSC_FALSE) { */
		/* 	Stg_SETERRQ( PETSC_ERR_SUP, "Only valid for MatType = block" ); */
		/* } */
		/* MatGetSize( A, &M, &N ); */
		/* if ( (M!=2) || (N!=2) ) { */
		/* 	Stg_SETERRQ2( PETSC_ERR_SUP, "Only valid for 2x2 block. Yours has dimension %Dx%D", M,N ); */
		/* } */

		VecDuplicate( x, &BA->Lz );
		VecDuplicate( x, &BA->Rz );

		BA->scaling_exists = PETSC_TRUE;
	}

//	if( BA->user_build_scaling != PETSC_NULL ) {
//		BA->user_build_scaling( A,b,x,BA->scaling_ctx);
//	}
//	else {
	BSSCR_MatStokesKBlockDefaultBuildScaling(BA,A,b,x, sym);
//	}

	BA->scalings_have_been_inverted = PETSC_FALSE;

	PetscFunctionReturn(0);
}

// updated
PetscErrorCode BSSCR_mat_kblock_invert_scalings( MatStokesBlockScaling BA )
{
	Vec L1,R1;

	VecNestGetSubVec( BA->Lz, 0, &L1 );
	VecNestGetSubVec( BA->Rz, 0, &R1 );

	VecReciprocal(L1);
	VecReciprocal(R1);

	/* toggle inversion flag */
	if( BA->scalings_have_been_inverted == PETSC_TRUE ) {   BA->scalings_have_been_inverted = PETSC_FALSE;  }
	if( BA->scalings_have_been_inverted == PETSC_FALSE ) {  BA->scalings_have_been_inverted = PETSC_TRUE;   }

	PetscFunctionReturn(0);
}

/* updated */
PetscErrorCode BSSCR_MatStokesKBlockScaleSystem( MatStokesBlockScaling BA, Mat A, Vec b, Vec x, Mat S, PetscTruth sym )
{
	if( BA->scaling_exists == PETSC_FALSE ) {
	    BSSCR_MatKBlock_ConstructScaling( BA,A,b,x, sym );
	}

	if( BA->scalings_have_been_inverted == PETSC_TRUE ) {
		BSSCR_mat_kblock_invert_scalings(BA); /* to undo inversion */
	}

	BSSCR_MatStokesKBlock_ApplyScaling(BA,A,b,x,S,sym);
	BA->system_has_been_scaled = PETSC_TRUE;

	/*
	PetscPrintf( PETSC_COMM_WORLD, "Post Scaling \n");
	MatBlock_ReportOperatorScales(A);
	*/

	PetscFunctionReturn(0);
}

// updated
PetscErrorCode BSSCR_MatStokesKBlockUnScaleSystem( MatStokesBlockScaling BA, Mat A, Vec b, Vec x, Mat S, PetscTruth sym )
{

	if( BA->system_has_been_scaled == PETSC_FALSE ) {
		printf("Warning: MatBlock has not been scaled !! \n");
		PetscFunctionReturn(0);
	}

	if( BA->scalings_have_been_inverted == PETSC_FALSE ) {
		BSSCR_mat_kblock_invert_scalings(BA);
	}

	BSSCR_MatStokesKBlock_ApplyScaling(BA,A,b,x,S,sym);
	BA->system_has_been_scaled = PETSC_FALSE;

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_MatStokesKBlockReportOperatorScales( Mat A, PetscTruth sym )
{
	Vec rA, rG;
	PetscInt loc;
	PetscReal min, max;
	Mat K,G,D,C;
	//PetscTruth is_block;

	/* check A is 2x2 block matrix */
	/* Stg_PetscTypeCompare( (PetscObject)A, "block", &is_block ); */
	/* if (is_block==PETSC_FALSE) { */
	/* 	Stg_SETERRQ( PETSC_ERR_SUP, "Only valid for MatType = block" ); */
	/* } */
	/* MatGetSize( A, &M, &N ); */
	/* if ( (M!=2) || (N!=2) ) { */
	/* 	Stg_SETERRQ2( PETSC_ERR_SUP, "Only valid for 2x2 block. Yours has dimension %Dx%D", M,N ); */
	/* } */

	MatNestGetSubMat( A, 0,0, &K );
	MatNestGetSubMat( A, 0,1, &G );
	MatNestGetSubMat( A, 1,0, &D );
	MatNestGetSubMat( A, 1,1, &C );

	MatGetVecs( K, PETSC_NULL, &rA );
	VecDuplicate( rA, &rG );

	/* Report the row max and mins */
	if (K!=PETSC_NULL) {
		MatGetRowMaxAbs( K, rA, PETSC_NULL );
		VecMax( rA, &loc, &max );
		PetscPrintf( PETSC_COMM_WORLD, "Sup_max(K) = %g \n", max );

		MatGetRowMinAbs( K, rA, PETSC_NULL );
		VecMin( rA, &loc, &min );
		PetscPrintf( PETSC_COMM_WORLD, "Sup_min(K) = %g \n\n", min );
	}

	if( G != PETSC_NULL ) {
		MatGetRowMaxAbs( G, rG, PETSC_NULL );
		VecMax( rG, &loc, &max );
		PetscPrintf( PETSC_COMM_WORLD, "Sup_max(G) = %g \n", max );

		MatGetRowMinAbs( G, rG, PETSC_NULL );
		VecMin( rG, &loc, &min );
		PetscPrintf( PETSC_COMM_WORLD, "Sup_min(G) = %g \n", min );
	}

	if( D != PETSC_NULL && !sym ) {
                Vec rD;

                MatGetVecs( D, PETSC_NULL, &rD );
		MatGetRowMaxAbs( D, rD, PETSC_NULL );
		VecMax( rD, &loc, &max );
		PetscPrintf( PETSC_COMM_WORLD, "Sup_max(D) = %g \n", max );

		MatGetRowMinAbs( D, rD, PETSC_NULL );
		VecMin( rD, &loc, &min );
		PetscPrintf( PETSC_COMM_WORLD, "Sup_min(D) = %g \n", min );

                Stg_VecDestroy(&rD );
	}

	if( C != PETSC_NULL ) {
		Vec cG;

		MatGetVecs( G, &cG, PETSC_NULL );
		MatGetRowMaxAbs( C, cG, PETSC_NULL );
		VecMax( cG, &loc, &max );
		PetscPrintf( PETSC_COMM_WORLD, "Sup_max(C) = %g \n", max );

		MatGetRowMinAbs( C, cG, PETSC_NULL );
		VecMin( cG, &loc, &min );
		PetscPrintf( PETSC_COMM_WORLD, "Sup_min(C) = %g \n\n", min );

		Stg_VecDestroy(&cG);
	}

	Stg_VecDestroy(&rA );
	Stg_VecDestroy(&rG );

	PetscFunctionReturn(0);
}

// updated
PetscErrorCode BSSCR_MatStokesKBlockDefaultBuildScaling( MatStokesBlockScaling BA, Mat A, Vec b, Vec x, PetscTruth sym )
{
	Mat K;
	Vec rA;
	Vec L1, R1;

	VecNestGetSubVec( BA->Lz, 0, &L1 );
	VecNestGetSubVec( BA->Rz, 0, &R1 );

	rA = L1;

	MatNestGetSubMat( A, 0,0, &K );

	/* Get diag of K */
	MatGetDiagonal( K, rA);
	VecSqrt( rA );
	VecReciprocal( rA );

	/* leave Shat alone as the Schur Compliment is unchanged by any consistent scaling of K */
	//MatGetDiagonal( Shat, rC);
	//VecSqrt( rC );
	//VecReciprocal( rC );

	/* no scaling for pressure */
	//VecSet( rC, 1.0 );

	VecCopy( L1, R1 );

	PetscFunctionReturn(0);
}


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

double interpolate_integral(interp_info *interp, double x_lo, double x_hi) {
    double x_lo_interp;
    double x_hi_interp;
    double x_lo_tmp;
    double x_hi_tmp;
    double x_min;
    double x_max;
    double y_x_min;
    double y_x_max;
    double r_val;
    double alpha;
    int    flag_error = GBP_FALSE;

    x_min   = ((interp_info *)interp)->x[0];
    x_max   = ((interp_info *)interp)->x[((interp_info *)interp)->n - 1];
    y_x_min = ((interp_info *)interp)->y[0];
    y_x_max = ((interp_info *)interp)->y[((interp_info *)interp)->n - 1];

    r_val       = 0.;
    x_lo_interp = x_lo;
    x_hi_interp = x_hi;

    if(x_hi > x_lo) {
        // Extrapolate integral to low-x assuming dI prop. to x^alpha
        if(x_lo < x_min) {
            x_lo_tmp = x_lo;
            x_hi_tmp = GBP_MIN(x_hi, x_min);
            alpha    = interpolate_derivative(interp, ((interp_info *)interp)->x[0]);
            if(alpha < -1. && x_lo_tmp * x_hi_tmp == 0.) {
                fprintf(stderr, "ERROR: integration extrapolation to low-x is divergent! (alpha=%lf)\n", alpha);
                flag_error = GBP_TRUE;
            } else
                r_val += (y_x_min / (alpha + 1.)) * (pow(x_hi_tmp / x_min, alpha + 1.) - pow(x_lo_tmp / x_min, alpha + 1.));
            x_lo_interp = x_min;
        }

        // Extrapolate integral to high-x assuming dI prop. to x^alpha
        if(x_hi > x_max) {
            x_lo_tmp = GBP_MAX(x_max, x_lo);
            x_hi_tmp = x_hi;
            alpha    = interpolate_derivative(interp, ((interp_info *)interp)->x[((interp_info *)interp)->n - 1]);
            if(alpha < 0. && x_lo_tmp * x_hi_tmp == 0.) {
                fprintf(stderr, "ERROR: integration extrapolation to high-x is divergent! (alpha=%lf)\n", alpha);
                flag_error = GBP_TRUE;
            } else
                r_val += (y_x_min / (alpha + 1.)) * (pow(x_hi_tmp / x_min, alpha + 1.) - pow(x_lo_tmp / x_min, alpha + 1.));
            x_hi_interp = x_max;
        }

        r_val += gsl_interp_eval_integ(((interp_info *)interp)->interp,
                                       ((interp_info *)interp)->x,
                                       ((interp_info *)interp)->y,
                                       x_lo_interp,
                                       x_hi_interp,
                                       ((interp_info *)interp)->accel);
    }
    return (r_val);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode