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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>

gsl_histogram *
gsl_histogram_alloc (size_t n)
{
  gsl_histogram *h;

  if (n == 0)
    {
      GSL_ERROR_VAL ("histogram length n must be positive integer",
                        GSL_EDOM, 0);
    }

  h = (gsl_histogram *) malloc (sizeof (gsl_histogram));

  if (h == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for histogram struct",
                        GSL_ENOMEM, 0);
    }

  h->range = (double *) malloc ((n + 1) * sizeof (double));

  if (h->range == 0)
    {
      free (h);         /* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for histogram ranges",
                        GSL_ENOMEM, 0);
    }

  h->bin = (double *) malloc (n * sizeof (double));

  if (h->bin == 0)
    {
      free (h->range);
      free (h);         /* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for histogram bins",
                        GSL_ENOMEM, 0);
    }

  h->n = n;

  return h;
}

static void
make_uniform (double range[], size_t n, double xmin, double xmax)
{
  size_t i;

  for (i = 0; i <= n; i++)
    {
      double f1 = ((double) (n-i) / (double) n);
      double f2 = ((double) i / (double) n);
      range[i] = f1 * xmin +  f2 * xmax;
    }
}

gsl_histogram *
gsl_histogram_calloc_uniform (const size_t n, const double xmin,
                              const double xmax)
{
  gsl_histogram *h;

  if (xmin >= xmax)
    {
      GSL_ERROR_VAL ("xmin must be less than xmax", GSL_EINVAL, 0);
    }

  h = gsl_histogram_calloc (n);

  if (h == 0)
    {
      return h;
    }

  make_uniform (h->range, n, xmin, xmax);

  return h;
}

gsl_histogram *
gsl_histogram_calloc (size_t n)
{
  gsl_histogram * h = gsl_histogram_alloc (n);

  if (h == 0)
    {
      return h;
    }

  {
    size_t i;

    for (i = 0; i < n + 1; i++)
      {
        h->range[i] = i;
      }

    for (i = 0; i < n; i++)
      {
        h->bin[i] = 0;
      }
  }

  h->n = n;

  return h;
}

void
gsl_histogram_free (gsl_histogram * h)
{
  RETURN_IF_NULL (h);
  free (h->range);
  free (h->bin);
  free (h);
}

/* These initialization functions suggested by Achim Gaedke */

int
gsl_histogram_set_ranges_uniform (gsl_histogram * h, double xmin, double xmax)
{
  size_t i;
  const size_t n = h->n;

  if (xmin >= xmax)
    {
      GSL_ERROR ("xmin must be less than xmax", GSL_EINVAL);
    }

  /* initialize ranges */

  make_uniform (h->range, n, xmin, xmax);

  /* clear contents */

  for (i = 0; i < n; i++)
    {
      h->bin[i] = 0;
    }

  return GSL_SUCCESS;
}

int
gsl_histogram_set_ranges (gsl_histogram * h, const double range[], size_t size)
{
  size_t i;
  const size_t n = h->n;

  if (size != (n+1))
    {
      GSL_ERROR ("size of range must match size of histogram", GSL_EINVAL);
    }

  /* initialize ranges */

  for (i = 0; i <= n; i++)
    {
      h->range[i] = range[i];
    }

  /* clear contents */

  for (i = 0; i < n; i++)
    {
      h->bin[i] = 0;
    }

  return GSL_SUCCESS;
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef _amici_TPL_MODELNAME_h
#define _amici_TPL_MODELNAME_h
#include <cmath>
#include <memory>
#include <gsl/gsl-lite.hpp>

#include "amici/model_ode.h"
#include "amici/solver_cvodes.h"

#include "sundials/sundials_types.h"

namespace amici {

class Solver;

namespace model_TPL_MODELNAME {

extern std::array<const char*, TPL_NP> parameterNames;
extern std::array<const char*, TPL_NK> fixedParameterNames;
extern std::array<const char*, TPL_NX_RDATA> stateNames;
extern std::array<const char*, TPL_NY> observableNames;
extern std::array<const char*, TPL_NW> expressionNames;
extern std::array<const char*, TPL_NP> parameterIds;
extern std::array<const char*, TPL_NK> fixedParameterIds;
extern std::array<const char*, TPL_NX_RDATA> stateIds;
extern std::array<const char*, TPL_NY> observableIds;
extern std::array<const char*, TPL_NW> expressionIds;

TPL_JY_DEF
TPL_DJYDSIGMA_DEF
TPL_DJYDY_DEF
TPL_DJYDY_COLPTRS_DEF
TPL_DJYDY_ROWVALS_DEF
TPL_ROOT_DEF
TPL_DWDP_DEF
TPL_DWDP_COLPTRS_DEF
TPL_DWDP_ROWVALS_DEF
TPL_DWDX_DEF
TPL_DWDX_COLPTRS_DEF
TPL_DWDX_ROWVALS_DEF
TPL_DWDW_DEF
TPL_DWDW_COLPTRS_DEF
TPL_DWDW_ROWVALS_DEF
TPL_DXDOTDW_DEF
TPL_DXDOTDW_COLPTRS_DEF
TPL_DXDOTDW_ROWVALS_DEF
TPL_DXDOTDP_EXPLICIT_DEF
TPL_DXDOTDP_EXPLICIT_COLPTRS_DEF
TPL_DXDOTDP_EXPLICIT_ROWVALS_DEF
TPL_DXDOTDX_EXPLICIT_DEF
TPL_DXDOTDX_EXPLICIT_COLPTRS_DEF
TPL_DXDOTDX_EXPLICIT_ROWVALS_DEF
TPL_DYDX_DEF
TPL_DYDP_DEF
TPL_SIGMAY_DEF
TPL_DSIGMAYDP_DEF
TPL_W_DEF
TPL_X0_DEF
TPL_X0_FIXEDPARAMETERS_DEF
TPL_SX0_DEF
TPL_SX0_FIXEDPARAMETERS_DEF
TPL_XDOT_DEF
TPL_Y_DEF
TPL_STAU_DEF
TPL_DELTAX_DEF
TPL_DELTASX_DEF
TPL_X_RDATA_DEF
TPL_X_SOLVER_DEF
TPL_TOTAL_CL_DEF

/**
 * @brief AMICI-generated model subclass.
 */
class Model_TPL_MODELNAME : public amici::Model_ODE {
  public:
    /**
     * @brief Default constructor.
     */
    Model_TPL_MODELNAME()
        : amici::Model_ODE(
              amici::ModelDimensions(
                  TPL_NX_RDATA,                            // nx_rdata
                  TPL_NXTRUE_RDATA,                        // nxtrue_rdata
                  TPL_NX_SOLVER,                           // nx_solver
                  TPL_NXTRUE_SOLVER,                       // nxtrue_solver
                  TPL_NX_SOLVER_REINIT,                    // nx_solver_reinit
                  TPL_NP,                                  // np
                  TPL_NK,                                  // nk
                  TPL_NY,                                  // ny
                  TPL_NYTRUE,                              // nytrue
                  TPL_NZ,                                  // nz
                  TPL_NZTRUE,                              // nztrue
                  TPL_NEVENT,                              // nevent
                  TPL_NOBJECTIVE,                          // nobjective
                  TPL_NW,                                  // nw
                  TPL_NDWDX,                               // ndwdx
                  TPL_NDWDP,                               // ndwdp
                  TPL_NDWDW,                               // ndwdw
                  TPL_NDXDOTDW,                            // ndxdotdw
                  TPL_NDJYDY,                              // ndjydy
                  0,                                       // nnz
                  TPL_UBW,                                 // ubw
                  TPL_LBW                                  // lbw
              ),
              amici::SimulationParameters(
                  std::vector<realtype>{TPL_FIXED_PARAMETERS}, // fixedParameters
                  std::vector<realtype>{TPL_PARAMETERS}        // dynamic parameters
              ),
              TPL_O2MODE,                                  // o2mode
              std::vector<realtype>(TPL_NX_SOLVER, 0.0),   // idlist
              std::vector<int>{},                          // z2event
              true,                                        // pythonGenerated
              TPL_NDXDOTDP_EXPLICIT,                       // ndxdotdp_explicit
              TPL_NDXDOTDX_EXPLICIT,                       // ndxdotdx_explicit
              TPL_W_RECURSION_DEPTH                        // w_recursion_depth
          ) {}

    /**
     * @brief Clone this model instance.
     * @return A deep copy of this instance.
     */
    virtual amici::Model *clone() const override {
        return new Model_TPL_MODELNAME(*this);
    }

    /** model specific implementation of fJrz
     * @param nllh regularization for event measurements z
     * @param iz event output index
     * @param p parameter vector
     * @param k constant vector
     * @param z model event output at timepoint
     * @param sigmaz event measurement standard deviation at timepoint
     **/
    virtual void fJrz(realtype *nllh, const int iz, const realtype *p,
                      const realtype *k, const realtype *rz,
                      const realtype *sigmaz) override {}

    TPL_JY_IMPL

    /** model specific implementation of fJz
     * @param nllh negative log-likelihood for event measurements z
     * @param iz event output index
     * @param p parameter vector
     * @param k constant vector
     * @param z model event output at timepoint
     * @param sigmaz event measurement standard deviation at timepoint
     * @param mz event measurements at timepoint
     **/
    virtual void fJz(realtype *nllh, const int iz, const realtype *p,
                     const realtype *k, const realtype *z,
                     const realtype *sigmaz, const realtype *mz) override {}

    /** model specific implementation of fdJrzdsigma
     * @param dJrzdsigma Sensitivity of event penalization Jrz w.r.t.
     * standard deviation sigmaz
     * @param iz event output index
     * @param p parameter vector
     * @param k constant vector
     * @param rz model root output at timepoint
     * @param sigmaz event measurement standard deviation at timepoint
     **/
    virtual void fdJrzdsigma(realtype *dJrzdsigma, const int iz,
                             const realtype *p, const realtype *k,
                             const realtype *rz,
                             const realtype *sigmaz) override {}

    /** model specific implementation of fdJrzdz
     * @param dJrzdz partial derivative of event penalization Jrz
     * @param iz event output index
     * @param p parameter vector
     * @param k constant vector
     * @param rz model root output at timepoint
     * @param sigmaz event measurement standard deviation at timepoint
     **/
    virtual void fdJrzdz(realtype *dJrzdz, const int iz, const realtype *p,
                         const realtype *k, const realtype *rz,
                         const realtype *sigmaz) override {}

    TPL_DJYDSIGMA_IMPL

    /** model specific implementation of fdJzdsigma
     * @param dJzdsigma Sensitivity of event measurement
     * negative log-likelihood Jz w.r.t. standard deviation sigmaz
     * @param iz event output index
     * @param p parameter vector
     * @param k constant vector
     * @param z model event output at timepoint
     * @param sigmaz event measurement standard deviation at timepoint
     * @param mz event measurement at timepoint
     **/
    virtual void fdJzdsigma(realtype *dJzdsigma, const int iz,
                            const realtype *p, const realtype *k,
                            const realtype *z, const realtype *sigmaz,
                            const realtype *mz) override {}

    /** model specific implementation of fdJzdz
     * @param dJzdz partial derivative of event measurement negative
     *log-likelihood Jz
     * @param iz event output index
     * @param p parameter vector
     * @param k constant vector
     * @param z model event output at timepoint
     * @param sigmaz event measurement standard deviation at timepoint
     * @param mz event measurement at timepoint
     **/
    virtual void fdJzdz(realtype *dJzdz, const int iz, const realtype *p,
                        const realtype *k, const realtype *z,
                        const realtype *sigmaz, const realtype *mz) override {}

    /** model specific implementation of fdeltasx
     * @param deltaqB sensitivity update
     * @param t current time
     * @param x current state
     * @param p parameter vector
     * @param k constant vector
     * @param h heaviside vector
     * @param ip sensitivity index
     * @param ie event index
     * @param xdot new model right hand side
     * @param xdot_old previous model right hand side
     * @param xB adjoint state
     **/
    virtual void fdeltaqB(realtype *deltaqB, const realtype t,
                          const realtype *x, const realtype *p,
                          const realtype *k, const realtype *h, const int ip,
                          const int ie, const realtype *xdot,
                          const realtype *xdot_old,
                          const realtype *xB) override {}

    TPL_DELTASX_IMPL

    TPL_DELTAX_IMPL

    /** model specific implementation of fdeltaxB
     * @param deltaxB adjoint state update
     * @param t current time
     * @param x current state
     * @param p parameter vector
     * @param k constant vector
     * @param h heaviside vector
     * @param ie event index
     * @param xdot new model right hand side
     * @param xdot_old previous model right hand side
     * @param xB current adjoint state
     **/
    virtual void fdeltaxB(realtype *deltaxB, const realtype t,
                          const realtype *x, const realtype *p,
                          const realtype *k, const realtype *h, const int ie,
                          const realtype *xdot, const realtype *xdot_old,
                          const realtype *xB) override {}

    /** model specific implementation of fdrzdp
     * @param drzdp partial derivative of root output rz w.r.t. model parameters
     *p
     * @param ie event index
     * @param t current time
     * @param x current state
     * @param p parameter vector
     * @param k constant vector
     * @param h heaviside vector
     * @param ip parameter index w.r.t. which the derivative is requested
     **/
    virtual void fdrzdp(realtype *drzdp, const int ie, const realtype t,
                        const realtype *x, const realtype *p, const realtype *k,
                        const realtype *h, const int ip) override {}

    /** model specific implementation of fdrzdx
     * @param drzdx partial derivative of root output rz w.r.t. model states x
     * @param ie event index
     * @param t current time
     * @param x current state
     * @param p parameter vector
     * @param k constant vector
     * @param h heaviside vector
     **/
    virtual void fdrzdx(realtype *drzdx, const int ie, const realtype t,
                        const realtype *x, const realtype *p, const realtype *k,
                        const realtype *h) override {}

    TPL_DSIGMAYDP_IMPL

    /** model specific implementation of fsigmaz
     * @param dsigmazdp partial derivative of standard deviation of event
     *measurements
     * @param t current time
     * @param p parameter vector
     * @param k constant vector
     * @param ip sensitivity index
     **/
    virtual void fdsigmazdp(realtype *dsigmazdp, const realtype t,
                            const realtype *p, const realtype *k,
                            const int ip) override {}

    TPL_DJYDY_IMPL
    TPL_DJYDY_COLPTRS_IMPL
    TPL_DJYDY_ROWVALS_IMPL

    TPL_DWDP_IMPL
    TPL_DWDP_COLPTRS_IMPL
    TPL_DWDP_ROWVALS_IMPL

    TPL_DWDX_IMPL
    TPL_DWDX_COLPTRS_IMPL
    TPL_DWDX_ROWVALS_IMPL

    TPL_DWDW_IMPL
    TPL_DWDW_COLPTRS_IMPL
    TPL_DWDW_ROWVALS_IMPL

    TPL_DXDOTDW_IMPL
    TPL_DXDOTDW_COLPTRS_IMPL
    TPL_DXDOTDW_ROWVALS_IMPL

    TPL_DXDOTDP_EXPLICIT_IMPL
    TPL_DXDOTDP_EXPLICIT_COLPTRS_IMPL
    TPL_DXDOTDP_EXPLICIT_ROWVALS_IMPL

    TPL_DXDOTDX_EXPLICIT_IMPL
    TPL_DXDOTDX_EXPLICIT_COLPTRS_IMPL
    TPL_DXDOTDX_EXPLICIT_ROWVALS_IMPL

    TPL_DYDX_IMPL

    TPL_DYDP_IMPL

    /** model specific implementation of fdzdp
     * @param dzdp partial derivative of event-resolved output z w.r.t. model
     *parameters p
     * @param ie event index
     * @param t current time
     * @param x current state
     * @param p parameter vector
     * @param k constant vector
     * @param h heaviside vector
     * @param ip parameter index w.r.t. which the derivative is requested
     **/
    virtual void fdzdp(realtype *dzdp, const int ie, const realtype t,
                       const realtype *x, const realtype *p, const realtype *k,
                       const realtype *h, const int ip) override {}

    /** model specific implementation of fdzdx
     * @param dzdx partial derivative of event-resolved output z w.r.t. model
     *states x
     * @param ie event index
     * @param t current time
     * @param x current state
     * @param p parameter vector
     * @param k constant vector
     * @param h heaviside vector
     **/
    virtual void fdzdx(realtype *dzdx, const int ie, const realtype t,
                       const realtype *x, const realtype *p, const realtype *k,
                       const realtype *h) override {}

    TPL_ROOT_IMPL

    /** model specific implementation of frz
     * @param rz value of root function at current timepoint (non-output events
     *not included)
     * @param ie event index
     * @param t current time
     * @param x current state
     * @param p parameter vector
     * @param k constant vector
     * @param h heaviside vector
     **/
    virtual void frz(realtype *rz, const int ie, const realtype t,
                     const realtype *x, const realtype *p, const realtype *k,
                     const realtype *h) override {}

    TPL_SIGMAY_IMPL

    /** model specific implementation of fsigmaz
     * @param sigmaz standard deviation of event measurements
     * @param t current time
     * @param p parameter vector
     * @param k constant vector
     **/
    virtual void fsigmaz(realtype *sigmaz, const realtype t, const realtype *p,
                         const realtype *k) override {}

    /** model specific implementation of fsrz
     * @param srz Sensitivity of rz, total derivative
     * @param ie event index
     * @param t current time
     * @param x current state
     * @param p parameter vector
     * @param k constant vector
     * @param sx current state sensitivity
     * @param h heaviside vector
     * @param ip sensitivity index
     **/
    virtual void fsrz(realtype *srz, const int ie, const realtype t,
                      const realtype *x, const realtype *p, const realtype *k,
                      const realtype *h, const realtype *sx,
                      const int ip) override {}

    TPL_STAU_IMPL
    TPL_SX0_IMPL
    TPL_SX0_FIXEDPARAMETERS_IMPL

    /** model specific implementation of fsz
     * @param sz Sensitivity of rz, total derivative
     * @param ie event index
     * @param t current time
     * @param x current state
     * @param p parameter vector
     * @param k constant vector
     * @param h heaviside vector
     * @param sx current state sensitivity
     * @param ip sensitivity index
     **/
    virtual void fsz(realtype *sz, const int ie, const realtype t,
                     const realtype *x, const realtype *p, const realtype *k,
                     const realtype *h, const realtype *sx,
                     const int ip) override {}

    TPL_W_IMPL

    TPL_X0_IMPL

    TPL_X0_FIXEDPARAMETERS_IMPL

    TPL_XDOT_IMPL

    TPL_Y_IMPL

    /** model specific implementation of fz
     * @param z value of event output
     * @param ie event index
     * @param t current time
     * @param x current state
     * @param p parameter vector
     * @param k constant vector
     * @param h heaviside vector
     **/
    virtual void fz(realtype *z, const int ie, const realtype t,
                    const realtype *x, const realtype *p, const realtype *k,
                    const realtype *h) override {}

    TPL_X_RDATA_IMPL

    TPL_X_SOLVER_IMPL

    TPL_TOTAL_CL_IMPL

    std::string getName() const override {
        return "TPL_MODELNAME";
    }

    /**
     * @brief Get names of the model parameters
     * @return the names
     */
    virtual std::vector<std::string> getParameterNames() const override {
        return std::vector<std::string>(parameterNames.begin(),
                                        parameterNames.end());
    }

    /**
     * @brief Get names of the model states
     * @return the names
     */
    virtual std::vector<std::string> getStateNames() const override {
        return std::vector<std::string>(stateNames.begin(), stateNames.end());
    }

    /**
     * @brief Get names of the fixed model parameters
     * @return the names
     */
    virtual std::vector<std::string> getFixedParameterNames() const override {
        return std::vector<std::string>(fixedParameterNames.begin(),
                                        fixedParameterNames.end());
    }

    /**
     * @brief Get names of the observables
     * @return the names
     */
    virtual std::vector<std::string> getObservableNames() const override {
        return std::vector<std::string>(observableNames.begin(),
                                        observableNames.end());
    }

    /**
     * @brief Get names of model expressions
     * @return Expression names
     */
    virtual std::vector<std::string> getExpressionNames() const override {
        return std::vector<std::string>(expressionNames.begin(),
                                        expressionNames.end());
    }

    /**
     * @brief Get ids of the model parameters
     * @return the ids
     */
    virtual std::vector<std::string> getParameterIds() const override {
        return std::vector<std::string>(parameterIds.begin(),
                                        parameterIds.end());
    }

    /**
     * @brief Get ids of the model states
     * @return the ids
     */
    virtual std::vector<std::string> getStateIds() const override {
        return std::vector<std::string>(stateIds.begin(), stateIds.end());
    }

    /**
     * @brief Get ids of the fixed model parameters
     * @return the ids
     */
    virtual std::vector<std::string> getFixedParameterIds() const override {
        return std::vector<std::string>(fixedParameterIds.begin(),
                                        fixedParameterIds.end());
    }

    /**
     * @brief Get ids of the observables
     * @return the ids
     */
    virtual std::vector<std::string> getObservableIds() const override {
        return std::vector<std::string>(observableIds.begin(),
                                        observableIds.end());
    }

    /**
     * @brief Get IDs of model expressions
     * @return Expression IDs
     */
    virtual std::vector<std::string> getExpressionIds() const override {
        return std::vector<std::string>(expressionIds.begin(),
                                        expressionIds.end());
    }

    /**
     * @brief function indicating whether reinitialization of states depending on
     fixed parameters is permissible
     * @return flag indicating whether reinitialization of states depending on
     fixed parameters is permissible
     */
    virtual bool isFixedParameterStateReinitializationAllowed() const override {
        return TPL_REINIT_FIXPAR_INITCOND;
    }

    /**
     * @brief returns the AMICI version that was used to generate the model
     * @return AMICI version string
     */
    virtual std::string getAmiciVersion() const override {
        return "TPL_AMICI_VERSION_STRING";
    }

    /**
     * @brief returns the amici version that was used to generate the model
     * @return AMICI git commit hash
     */
    virtual std::string getAmiciCommit() const override {
        return "TPL_AMICI_COMMIT_STRING";
    }

    virtual bool hasQuadraticLLH() const override {
        return TPL_QUADRATIC_LLH;
    }
};

} // namespace model_TPL_MODELNAME

} // namespace amici

#endif /* _amici_TPL_MODELNAME_h */

ext: h

max_stars_repo_name: jvanhoefer/AMICI

max_issues_repo_name: jvanhoefer/AMICI

max_forks_repo_name: jvanhoefer/AMICI


text:
#include <viaio/VImage.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include "gsl_utils.h"

#define ABS(x) ((x) >= 0) ? (x) : -(x)
#define MIN(a,b) ((a) < (b)) ? (a) : (b)
#define MAX(a,b) ((a) > (b)) ? (a) : (b)
#define ROUND(x) (int)(((x) >= 0) ? ((x) + 0.5) : ((x) - 0.5))

void whitecov2(VImage effect_image,VImage rho_vol,
	       gsl_matrix_float* Y, gsl_matrix_float* X, gsl_matrix_float* con,
	       VFloat* Dfs, int numlags, int slice)
{
  int numcon = 1;
  int n = Y->size1;
  VFloat Df = Dfs[numcon];

  int i=0,j=0,k=0;

  /* some checks and warnings */
  if(X->size2 != n)
    VError("whitecov2: Warning! sizes of design and data matrix do not match, n=%d, %ld",n,X->size2);

  if(X->size1 != con->size2)
    VError("whitecov2: Warning! contrast and design do not match.");

  gsl_vector_float* irho = gsl_vector_float_alloc(VImageNBands(rho_vol));
  if(numlags == 1) {
    float* p = irho->data;
    const float drho = 0.01;
    for(i=0;i<irho->size;i++) {
      float tmp = ROUND((VPixel(rho_vol, i, slice, 0, VFloat) / drho)) * drho;
      tmp = MIN(tmp, (1-drho));
      tmp = MAX(tmp, (-1+drho));
      *p++ = tmp;
    }
  }
  else {
    float* p = irho->data;
    for(i=1;i<=irho->size;i++)
      *p++ = i;
  }

  /* : X=X';  */
  gsl_matrix_float* transX = gsl_matrix_float_alloc(X->size2, X->size1);
  gsl_matrix_float_transpose_memcpy(transX, X);

  /* : Xstar=X; */
  gsl_matrix_float* XStar = gsl_matrix_float_alloc(transX->size1, transX->size2);
  gsl_matrix_float_memcpy(XStar, transX);

  /* END allocate memory buffer */

  /****************************** order == 1 *****************************/
  if (numlags == 1) {

    /* allocate memory buffers */
    gsl_vector_float* pixBig = gsl_vector_float_alloc(irho->size);
    gsl_matrix_float* pinvXStar = gsl_matrix_float_alloc(XStar->size2, XStar->size1);
    gsl_matrix_float* V = gsl_matrix_float_alloc(pinvXStar->size1, pinvXStar->size1);
    gsl_matrix_float* buffer_01 = gsl_matrix_float_alloc(con->size1,V->size2);
    gsl_matrix_float* buffer_02 = gsl_matrix_float_alloc(con->size1,con->size1);
    gsl_matrix_float* cVcinv = gsl_matrix_float_alloc(con->size1, con->size1);

    gsl_vector_float* u = funique(irho);
    float* rho = u->data;
    for(i=0;i<u->size;i++) {

      /* :  pix=int16(find(irho==rho)); */
      float* p = irho->data;
      float* ppB = pixBig->data;
      int counter = 0;
      for(j=0;j<irho->size;j++) {
	if(*p++ == *rho) {
	  *ppB++ = (float)j;
	  counter++;
	}
      }
      gsl_vector_float* pix = gsl_vector_float_alloc(counter);
      ppB = pixBig->data;
      float* pp = pix->data;
      for(j=0;j<counter;j++) *pp++ = *ppB++;

      /* : Ystar=Y(:,pix); */
      gsl_matrix_float* YStar = fmat_subcols(Y,pix);

      /* : factor=1./sqrt(1-rho^2); */
      float factor = pow(1 - *rho * *rho, -0.5);

      /* : Ystar(k1,:)=(Y(k1,pix)-rho*Y(k1-1,pix))*factor; */
      /* point to last element */
      float* pStar = YStar->data+(YStar->size1*YStar->size2)-1;
      /* traverse matrix in reverse order */
      for (j = 0; j < (YStar->size1-1)*YStar->size2; ++j) {
	*pStar = (*pStar - *rho * *(pStar - YStar->size2)) * factor;
	pStar--;
      }

      /* : Xstar(k1,:)=(X(k1,:)-rho*X(k1-1,:))*factor; */
      pStar = XStar->data+XStar->size2;
      float* pX = transX->data+transX->size2;
      for (j = 0; j < (XStar->size1-1)*XStar->size2; ++j) {
	*pStar++ = (*pX - *rho * *(pX - transX->size2)) * factor;
	pX++;
      }

      /* : pinvXstar=pinv(Xstar); */
      fmat_PseudoInv(XStar,pinvXStar);

      /* : betahat=pinvXstar*Ystar; */
      gsl_matrix_float* betahat = fmat_x_mat(pinvXStar, YStar,NULL);

      /* : resid=Ystar-Xstar*betahat; */
      gsl_matrix_float* buffer = fmat_x_mat(XStar, betahat, NULL);
      gsl_matrix_float_sub(YStar, buffer);
      gsl_matrix_float* resid =  YStar;
      gsl_matrix_float_free(buffer);

      /* : SSE=sum(resid.^2,1); */
      buffer = gsl_matrix_float_alloc(resid->size1, resid->size2);
      gsl_matrix_float_memcpy (buffer, resid);
      gsl_matrix_float_mul_elements (buffer, buffer);
      gsl_vector_float* sse = fsum(buffer, 1,NULL);
      gsl_matrix_float_free(buffer);

      /* : sd=sqrt(SSE/Df); */
      gsl_vector_float* sd = gsl_vector_float_alloc(sse->size);
      float* pSSE = sse->data;
      float* pSD = sd->data;
      for (j = 0; j < sse->size; ++j) {
	*pSD++ = (float)sqrt(*pSSE++/Df);
      }

      /* :   V=pinvXstar*pinvXstar'; */
      fmat_x_matT(pinvXStar, pinvXStar, V);

      /* : mag_ef=contrast*betahat; */
      gsl_matrix_float* mag_ef = fmat_x_mat(con, betahat, NULL);

      /* : mag_sd=sqrt(diag(contrast*V*contrast'))*sd; */
      fmat_x_mat(con,V,buffer_01);
      fmat_x_matT(buffer_01, con, buffer_02);

      gsl_vector_float_view diag = gsl_matrix_float_diagonal(buffer_02);
      gsl_matrix_float* mag_sd = gsl_matrix_float_alloc(diag.vector.size, sd->size);
      float* pMag = mag_sd->data;
      for (j = 0; j < diag.vector.size; ++j) {
	float s = sqrt(gsl_vector_float_get(&diag.vector,j));
	pSD = sd->data;
	for (k = 0; k < sd->size; ++k) {
	  *pMag++ = s * (*pSD++);
	}
      }

      /* effect_slice(pix,1)= (mag_ef./(mag_sd+(mag_sd<=0)).*(mag_sd>0))'; */
      VFloat* peff = VPixelPtr(effect_image,slice,0,0);
      float* pmeff = mag_ef->data;
      float* pmsd = mag_sd->data;
      float* ppix = pix->data;
      for (j=0; j<pix->size; ++j) {
	*(peff+(int)*ppix++) = *pmeff / (*pmsd+(*pmsd<=0)) * (*pmsd>0);
	pmeff++; pmsd++;
      }

      /* free some memory */
      gsl_matrix_float_free(betahat);
      gsl_matrix_float_free(YStar);
      gsl_matrix_float_free(mag_ef);
      gsl_matrix_float_free(mag_sd);
      gsl_vector_float_free(pix);
      gsl_vector_float_free(sse);
      gsl_vector_float_free(sd);

      rho++;
    } /* END for(i=0;i<u->size;i++)  */

    /* free some memory */
    gsl_vector_float_free(pixBig);
    gsl_matrix_float_free(pinvXStar);
    gsl_matrix_float_free(V);
    gsl_matrix_float_free(buffer_01);
    gsl_matrix_float_free(buffer_02);
    gsl_matrix_float_free(cVcinv);

  } /******************* END order == 1 **********************************/

/*************************** order > 1 *********************************/
    else{

      /* allocate memory buffer */
      gsl_vector_float* Coradj_pix = gsl_vector_float_alloc(VImageNColumns(rho_vol)+1);
      gsl_matrix_float* Ainv = gsl_matrix_float_alloc(Coradj_pix->size, Coradj_pix->size);
      gsl_matrix* dbuff = gsl_matrix_alloc(Ainv->size1,Ainv->size2);
      gsl_matrix_float* A = gsl_matrix_float_alloc(Ainv->size1,Ainv->size2);
      int nl = Ainv->size2;
      gsl_matrix_float* buffer_01 = gsl_matrix_float_alloc(n-nl,1);
      gsl_matrix_float* buffer_02 = gsl_matrix_float_alloc(A->size1,1);
      gsl_matrix_float* buffer_03 = gsl_matrix_float_alloc(A->size1,transX->size2);
      gsl_matrix_float* B = gsl_matrix_float_alloc(buffer_01->size1, A->size2);
      gsl_matrix_float* Vmhalf = gsl_matrix_float_alloc(n-nl,n);
      gsl_matrix_float* buffer_04 = gsl_matrix_float_alloc(Vmhalf->size1, transX->size2);
      gsl_matrix_float* YStar = gsl_matrix_float_alloc(n,1);
      gsl_matrix_float* pinvXStar = gsl_matrix_float_alloc(XStar->size2, XStar->size1);
      gsl_matrix_float* betahat = gsl_matrix_float_alloc(pinvXStar->size1, YStar->size2);
      gsl_matrix_float* buffer_05 = gsl_matrix_float_alloc(XStar->size1, betahat->size2);
      gsl_vector_float* sse = gsl_vector_float_alloc(1);
      gsl_vector_float* sd = gsl_vector_float_alloc(YStar->size2);
      gsl_matrix_float* V = gsl_matrix_float_alloc(pinvXStar->size1, pinvXStar->size1);
      gsl_matrix_float* mag_ef = gsl_matrix_float_alloc(con->size1, betahat->size2);
      gsl_matrix_float* buffer_06 = gsl_matrix_float_alloc(con->size1, V->size2);
      gsl_matrix_float* buffer_07 = gsl_matrix_float_alloc(con->size1, con->size1);
      gsl_matrix_float* mag_sd = gsl_matrix_float_alloc(con->size1, sd->size);
      gsl_matrix_float* buffer_08 = gsl_matrix_float_alloc(mag_ef->size1, mag_ef->size2);
      gsl_vector_float* sst = gsl_vector_float_alloc(1);
      /* END allocate memory buffer */

     /* we know that irho only contains a list of indices for every voxel so we can skip
        the unique call and iterate over irho itself */
      float* rho = irho->data;
      for(i=0;i<irho->size;i++) {

        /* : Coradj_pix=squeeze(rho_vol(pix,slice,:)); */
        float* pCoradj = Coradj_pix->data+1;
        for(k=0;k<Coradj_pix->size-1;k++) {
          *(pCoradj++) = VPixel(rho_vol,
                                i,
                                slice,
                                k,
                                VFloat);
        }

        /* : [Ainvt posdef]=chol(toeplitz([1 Coradj_pix'])); */

        /* at first the toeplitz matrix */
        Coradj_pix->data[0] = 1;
        fmat_toeplitz(Coradj_pix, Ainv);
        /* the cholesky decomposition */
        /* double buffer */
        for(j=0;j<Ainv->size1;j++) {
          for(k=0;k<Ainv->size2;k++) {
            gsl_matrix_set(dbuff,j,k,(double)gsl_matrix_float_get(Ainv,j,k));
          }
        }

        if(gsl_linalg_cholesky_decomp(dbuff) == GSL_EDOM)
          VError("Calculation error, cholesky decomposition failed for pix=%d",i);

        /* note: the gsl cholesky factorization returns the inverse
         * of Ainvt in the lower triangular part of it's output matrix.
         * Since we will use the inverse matrix of Ainvt in the next steps we
         * will work with the lower triangular part in contrast to the matlab algorithm */

        /* convert back to float and remove upper triangular matrix */
        for(j=0;j<Ainv->size1;j++) {
          for(k=0;k<Ainv->size2;k++) {
            if(j<k) {
              gsl_matrix_float_set(Ainv,j,k,0);
            }
            else {
              gsl_matrix_float_set(Ainv,j,k,(float)gsl_matrix_get(dbuff,j,k));
            }
          }
        }

        /* :  A=inv(Ainvt'); */
        fInv(Ainv, A);

        /* : B=ones(n-nl,1)*A(nl,:); */
	gsl_matrix_float_set_all(buffer_01,1);
        gsl_matrix_float_view subm = gsl_matrix_float_submatrix(A,nl-1,0,1,A->size2);
        fmat_x_mat(buffer_01,&subm.matrix,B);

        /* : Vmhalf=spdiags(double(B),double(1:nl),double(n-nl),double(n)); */
        gsl_matrix_float_set_zero(Vmhalf);
        for(j=0;j<B->size1;j++) {
          for(k=0;k<B->size2;k++) {
            gsl_matrix_float_set(Vmhalf,j,k+1+j,gsl_matrix_float_get(B,j,k));
          }
        }

        /* : Ystar=single(zeros(n,1)); */
        gsl_matrix_float_set_zero(YStar);

        /* : Ystar(1:nl)=A*Y(1:nl,pix); */
        subm = gsl_matrix_float_submatrix(Y,0,i,nl,1);
        fmat_x_mat(A,&subm.matrix,buffer_02);
        for(j=0;j<nl;j++) {
          gsl_matrix_float_set(YStar,j,0,gsl_matrix_float_get(buffer_02,j,0));
        }

	/* : Ystar((nl+1):n)=single(Vmhalf*double(Y(:,pix))); */
	subm = gsl_matrix_float_submatrix(Y,0,i,Y->size1,1);
	fmat_x_mat(Vmhalf,&subm.matrix,buffer_01);
	for(j=nl;j<n;j++) {
	  gsl_matrix_float_set(YStar,j,0,gsl_matrix_float_get(buffer_01,j-nl,0));
	}

        /* : Xstar(1:nl,:)=A*X(1:nl,:); */
        subm = gsl_matrix_float_submatrix(transX,0,0,nl,transX->size2);
        fmat_x_mat(A, &subm.matrix, buffer_03);
        for(j=0;j<buffer_03->size1;j++) {

          for(k=0;k<buffer_03->size2;k++) {
            gsl_matrix_float_set(XStar,j,k,gsl_matrix_float_get(buffer_03,j,k));
          }
        }

        /* : Xstar((nl+1):n,:)=single(Vmhalf*double(X)); */
        fmat_x_mat(Vmhalf,transX,buffer_04);
        for(j=nl;j<n;j++) {
          for(k=0;k<transX->size2;k++) {
            gsl_matrix_float_set(XStar,j,k,gsl_matrix_float_get(buffer_04,j-nl,k));
          }
        }

        /* : pinvXstar=pinv(Xstar); */
        fmat_PseudoInv(XStar,pinvXStar);

        /* : betahat=pinvXstar*Ystar; */
        fmat_x_mat(pinvXStar, YStar,betahat);

        /* : resid=Ystar-Xstar*betahat; */
        fmat_x_mat(XStar, betahat, buffer_05);
        gsl_matrix_float_sub(YStar, buffer_05);
        gsl_matrix_float* resid =  YStar;

        /* : SSE=sum(resid.^2,1); */
        gsl_matrix_float_mul_elements (resid, resid);
        fsum(resid, 1,sse);

        /* : sd=sqrt(SSE/Df); */
        float* pSSE = sse->data;
        float* pSD = sd->data;
        for (j = 0; j < sse->size; ++j) {
          *pSD++ = (float)sqrt(*pSSE++/Df);
        }

        /* :   V=pinvXstar*pinvXstar'; */
        fmat_x_matT(pinvXStar, pinvXStar, V);

        /* : mag_ef=contrast*betahat; */
        fmat_x_mat(con, betahat, mag_ef);

        /* : mag_sd=sqrt(diag(contrast*V*contrast'))*sd; */
        fmat_x_mat(con,V,buffer_06);
        fmat_x_matT(buffer_06, con, buffer_07);

        gsl_vector_float_view diag = gsl_matrix_float_diagonal(buffer_07);
        float* pMag = mag_sd->data;
        for (j = 0; j < diag.vector.size; ++j) {
          float s = sqrt(gsl_vector_float_get(&diag.vector,j));
          pSD = sd->data;
          for (k = 0; k < sd->size; ++k) {
            *pMag++ = s * *pSD++;
          }
        }

	/* effect_slice(pix,1)= (mag_ef./(mag_sd+(mag_sd<=0)).*(mag_sd>0))'; */
	VFloat* peff = VPixelPtr(effect_image, slice, 0, 0);
	float* pmeff = mag_ef->data;
	float* pmsd = mag_sd->data;
	*(peff+i) = *pmeff / (*pmsd+(*pmsd<=0)) * (*pmsd>0);

        /* next value */
        rho++;
      } /* END for(i=0;i<irho->size;i++) */

      /* free some memory */
      gsl_vector_float_free(Coradj_pix);
      gsl_matrix_float_free(Ainv);
      gsl_matrix_free(dbuff);
      gsl_matrix_float_free(A);
      gsl_matrix_float_free(buffer_01);
      gsl_matrix_float_free(buffer_02);
      gsl_matrix_float_free(buffer_03);
      gsl_matrix_float_free(buffer_04);
      gsl_matrix_float_free(B);
      gsl_matrix_float_free(Vmhalf);
      gsl_matrix_float_free(YStar);
      gsl_matrix_float_free(pinvXStar);
      gsl_matrix_float_free(buffer_05);
      gsl_vector_float_free(sd);
      gsl_matrix_float_free(V);
      gsl_matrix_float_free(mag_ef);
      gsl_matrix_float_free(buffer_06);
      gsl_matrix_float_free(buffer_07);
      gsl_matrix_float_free(mag_sd);
      gsl_matrix_float_free(buffer_08);
      gsl_vector_float_free(sse);
      gsl_vector_float_free(sst);

    }/************************************** END order > 1 ***************************************/

  /* :
   * minus = ones(size(effect_slice),'single');
   * minus(find(effect_slice<0))=single(-1.0);
   * effect_slice=single(sqrt(Df .* log(1 + effect_slice.^2 ./ Df) .* (1 - 0.5 ./ Df)));
   * effect_slice=effect_slice .* minus;
   */
  double u=0;
  VFloat* peff = VPixelPtr(effect_image, slice,0,0);
  for (i = 0; i < VImageNColumns(effect_image) * VImageNRows(effect_image); ++i) {
    VFloat oldVal = *peff;
    Df = Dfs[0];
    u = (double)(*peff);
    u = sqrt(Df * log(1.0 + pow(u,2.0) / Df) * (1 - 0.5 / Df));
    u = (oldVal < 0) ? u * (-1.0) : u;
    if (gsl_isnan(u) || gsl_isinf(u)) u = 0;
    (*peff) = (float)u;
    peff++;
  }

  /* free allocated memory */
  gsl_matrix_float_free(XStar);
  gsl_matrix_float_free(transX);
  /* END free allocated memory */
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <math.h>
#include <stdio.h>
#include <mpfr.h>
#include <petsc.h>
#include <stdlib.h>
#include <string.h>

#include "../ellipsoid/ellipsoid.h"

#undef __FUNCT__
#define __FUNCT__ "TestNormalizationMPFR"
PetscErrorCode TestNormalizationMPFR()
{
  PetscErrorCode ierr;
  EllipsoidalSystem e;
  const PetscInt prec = 128;
  PetscReal xA = 3.0;
  PetscReal yB = 2.0;
  PetscReal zC = 1.0;
  mpfr_t a, b, c;
  mpfr_t a2, b2, c2;
  mpfr_t a3, b3, c3;
  mpfr_t temp1, temp2, temp3;
  mpfr_t pi;

  PetscFunctionBegin;
  mpfr_set_default_prec(4*prec);

  mpfr_inits(a, b, c, NULL);
  mpfr_inits(a2, b2, c2, NULL);
  mpfr_inits(a3, b3, c3, NULL);
  mpfr_inits(temp1, temp2, temp3, pi, NULL);

  mpfr_set_d(a, xA, MPFR_RNDN);
  mpfr_set_d(b, yB, MPFR_RNDN);
  mpfr_set_d(c, zC, MPFR_RNDN);

  //init ellipsoidal system
  printf("initEllipsoidalSystem...\n");
  ierr = initEllipsoidalSystem(&e, xA, yB, zC, prec);CHKERRQ(ierr);
  printf("done.\n");

  mpfr_mul(a2, a, a, MPFR_RNDN);
  mpfr_mul(a3, a2, a, MPFR_RNDN);
  mpfr_mul(b2, b, b, MPFR_RNDN);
  mpfr_mul(b3, b2, b, MPFR_RNDN);
  mpfr_mul(c2, c, c, MPFR_RNDN);
  mpfr_mul(c3, c2, c, MPFR_RNDN);

  mpfr_t firstTerm, secondTerm, LambdaD, LambdaDprime;
  mpfr_inits(firstTerm, secondTerm, LambdaD, LambdaDprime, NULL);

  //Dassios (B14)
  //firstTerm = (a*a + b*b + c*c)/3.0;
  mpfr_add(firstTerm, a2, b2, MPFR_RNDN);
  mpfr_add(firstTerm, firstTerm, c2, MPFR_RNDN);
  mpfr_div_d(firstTerm, firstTerm, 3.0, MPFR_RNDN);

  //double secondTerm = sqrt((a*a*a*a - b*b*c*c) + (b*b*b*b - a*a*c*c) + (c*c*c*c - a*a*b*b))/3.0;
  mpfr_mul(temp1, a3, a, MPFR_RNDN);
  mpfr_mul(temp2, b2, c2, MPFR_RNDN);
  mpfr_sub(temp1, temp1, temp2, MPFR_RNDN);
  mpfr_mul(temp2, b3, b, MPFR_RNDN);
  mpfr_mul(temp3, a2, c2, MPFR_RNDN);
  mpfr_sub(temp2, temp2, temp3, MPFR_RNDN);
  mpfr_add(temp1, temp1, temp2, MPFR_RNDN);
  mpfr_mul(temp2, c3, c, MPFR_RNDN);
  mpfr_mul(temp3, a2, b2, MPFR_RNDN);
  mpfr_sub(temp2, temp2, temp3, MPFR_RNDN);
  mpfr_add(temp1, temp1, temp2, MPFR_RNDN);
  mpfr_sqrt(secondTerm, temp1, MPFR_RNDN);

  //double LambdaD = firstTerm + secondTerm;
  mpfr_add(LambdaD, firstTerm, secondTerm, MPFR_RNDN);

  //double LambdaDprime = firstTerm - secondTerm;
  mpfr_sub(LambdaDprime, firstTerm, secondTerm, MPFR_RNDN);

  mpfr_t hx, hy, hz;
  mpfr_inits(hx, hy, hz, NULL);

  //Dassios (B16)-(B20)
  //double hx = sqrt(b*b - c*c);
  mpfr_sub(temp1, b2, c2, MPFR_RNDN);
  mpfr_sqrt(hx, temp1, MPFR_RNDN);

  //double hy = e.k;
  mpfr_set(hy, e.hp_k, MPFR_RNDN);
  //double hz = e.h;
  mpfr_set(hz, e.hp_h, MPFR_RNDN);

  //set pi
  mpfr_const_pi(pi, MPFR_RNDN);

  mpfr_t analytic[9];
  mpfr_t approx[9], doubleApprox[9];
  for(PetscInt k=0; k < 9; ++k)
    mpfr_inits(analytic[k], approx[k], doubleApprox[k], NULL);

  //analytic[0] = 4*PETSC_PI
  mpfr_mul_d(analytic[0], pi, 4.0, MPFR_RNDN);
  printf("analytic[0] = %4.4e\n", mpfr_get_d(analytic[0], MPFR_RNDN));

  //analytic[1] = 4*PETSC_PI/3 * hy*hy*hz*hz
  mpfr_div_d(temp1, pi, 3.0, MPFR_RNDN);
  mpfr_mul_d(temp1, temp1, 4.0, MPFR_RNDN);
  mpfr_mul(temp2, hy, hy, MPFR_RNDN);
  mpfr_mul(temp2, temp2, hz, MPFR_RNDN);
  mpfr_mul(temp2, temp2, hz, MPFR_RNDN);
  mpfr_mul(analytic[1], temp1, temp2, MPFR_RNDN);

  //analytic[2] = 4*PETSC_PI/3 * hx*hx*hz*hz
  mpfr_div_d(temp1, pi, 3.0, MPFR_RNDN);
  mpfr_mul_d(temp1, temp1, 4.0, MPFR_RNDN);
  mpfr_mul(temp2, hx, hx, MPFR_RNDN);
  mpfr_mul(temp2, temp2, hz, MPFR_RNDN);
  mpfr_mul(temp2, temp2, hz, MPFR_RNDN);
  mpfr_mul(analytic[2], temp1, temp2, MPFR_RNDN);

  //analytic[3] = 4*PETSC_PI/3 * hx*hx*hy*hy
  mpfr_div_d(temp1, pi, 3.0, MPFR_RNDN);
  mpfr_mul_d(temp1, temp1, 4.0, MPFR_RNDN);
  mpfr_mul(temp2, hx, hx, MPFR_RNDN);
  mpfr_mul(temp2, temp2, hy, MPFR_RNDN);
  mpfr_mul(temp2, temp2, hy, MPFR_RNDN);
  mpfr_mul(analytic[3], temp1, temp2, MPFR_RNDN);

  //analytic[4] = -8*PETSC_PI/5 * (LambdaD - LambdaDprime)*(LambdaD - a*a)*(LambdaD - b*b)*(LambdaD - c*c)

  //analytic[5] = 8*PETSC_PI/5 * (LambdaD - LambdaDprime)*(LambdaDprime - a*a)*(LambdaDprime - b*b)*(LambdaDprime - c*c),

  ierr = calcNormalizationMPFR(&e, 0, 0, approx+0);CHKERRQ(ierr);
  ierr = calcNormalization    (&e, 0, 0, doubleApprox+0);CHKERRQ(ierr);
  ierr = calcNormalizationMPFR(&e, 1, 0, approx+1);CHKERRQ(ierr);
  ierr = calcNormalizationMPFR(&e, 1, 1, approx+2);CHKERRQ(ierr);
  mpfr_t err[9];
  for(PetscInt k=0; k < 4; ++k) {
    mpfr_init(err[k]);
    mpfr_sub(err[k], approx[k], analytic[k], MPFR_RNDN);
    mpfr_div(err[k], err[k], analytic[k], MPFR_RNDN);
    mpfr_abs(err[k], err[k], MPFR_RNDN);
    mpfr_log10(err[k], err[k], MPFR_RNDN);
    printf("the error is %2.2f\n", mpfr_get_d(err[k], MPFR_RNDN));
  }

  for(PetscInt k=0; k < 9; ++k)
    mpfr_clears(analytic[k], approx[k], err[k], doubleApprox[k], NULL);
  mpfr_clears(hx, hy, hz, NULL);
  mpfr_clears(temp1, temp2, temp3, pi, NULL);
  mpfr_clears(a3, b3, c3, NULL);
  mpfr_clears(a2, b2, c2, NULL);
  mpfr_clears(a, b, c, NULL);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT __ "main"
PetscErrorCode main(int argc, char **argv)
{
  PetscErrorCode ierr;

  PetscFunctionBeginUser;

  ierr = PetscInitialize(&argc, &argv, NULL, NULL);CHKERRQ(ierr);

  ierr = TestNormalizationMPFR();

  ierr = PetscFinalize();

}

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation


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

#include <config.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

/* tests a single evaluator function from the low-level interface */
static int
test_single_low_level(
  double (*evaluator)(const gsl_interp2d *, const double[], const double[],
                      const double[], const double, const double,
                      gsl_interp_accel *, gsl_interp_accel *),
  int (*evaluator_e)(const gsl_interp2d *, const double[], const double[],
                     const double[], const double, const double,
                     gsl_interp_accel *, gsl_interp_accel *, double *),
  const gsl_interp2d * interp, const double xarr[], const double yarr[],
  const double zarr[], const double x, const double y,
  gsl_interp_accel * xa, gsl_interp_accel * ya,
  const double expected_results[], size_t i)
{
  if (expected_results != NULL)
    {
      int status;
      double result = evaluator(interp, xarr, yarr, zarr, x, y, xa, ya);
      gsl_test_rel(result, expected_results[i], 1e-10,
                   "low level %s %d", gsl_interp2d_name(interp), i);

      status = evaluator_e(interp, xarr, yarr, zarr, x, y, xa, ya, &result);
      if (status == GSL_SUCCESS)
        {
          gsl_test_rel(result, expected_results[i], 1e-10,
                       "low level _e %s %d", gsl_interp2d_name(interp), i);
        }
    }

  return 0;
}

/* tests a single evaluator function from the high-level interface */
static int
test_single_high_level(
    double (*evaluator)(const gsl_spline2d *, const double, const double,
                        gsl_interp_accel *, gsl_interp_accel *),
    int (*evaluator_e)(const gsl_spline2d *, const double, const double,
                       gsl_interp_accel *, gsl_interp_accel *, double *),
    const gsl_spline2d * interp, const double x, const double y,
    gsl_interp_accel * xa, gsl_interp_accel * ya,
    const double expected_results[], size_t i)
{
  if (expected_results != NULL)
    {
      int status;
      double result = evaluator(interp, x, y, xa, ya);

      gsl_test_rel(result, expected_results[i], 1e-10,
                   "high level %s %d", gsl_spline2d_name(interp), i);

      status = evaluator_e(interp, x, y, xa, ya, &result);
      if (status == GSL_SUCCESS)
        {
          gsl_test_rel(result, expected_results[i], 1e-10,
                       "high level _e %s %d", gsl_spline2d_name(interp), i);
        }
    }

  return 0;
}

/*
 * Tests that a given interpolation type reproduces the data points
 * it is given, and then tests that it correctly reproduces additional
 * values.
 */
static int
test_interp2d(
const double xarr[], const double yarr[], const double zarr[], /* data */
              size_t xsize, size_t ysize,                /* array sizes */
              const double xval[], const double yval[],  /* test points */
              const double zval[],                       /* expected results */
              const double zxval[], const double zyval[],
              const double zxxval[], const double zyyval[], const double zxyval[],
              size_t test_size,                          /* number of test points */
              const gsl_interp2d_type * T)
{
  gsl_interp_accel *xa = gsl_interp_accel_alloc();
  gsl_interp_accel *ya = gsl_interp_accel_alloc();
  int status = 0;
  size_t xi, yi, zi, i;
  gsl_interp2d * interp = gsl_interp2d_alloc(T, xsize, ysize);
  gsl_spline2d * interp_s = gsl_spline2d_alloc(T, xsize, ysize);
  unsigned int min_size = gsl_interp2d_type_min_size(T);

  gsl_test_int(min_size, T->min_size, "gsl_interp2d_type_min_size on %s", gsl_interp2d_name(interp));

  gsl_interp2d_init(interp, xarr, yarr, zarr, xsize, ysize);
  gsl_spline2d_init(interp_s, xarr, yarr, zarr, xsize, ysize);

  /* First check that the interpolation reproduces the given points */
  for (xi = 0; xi < xsize; xi++)
    {
      double x = xarr[xi];
      for (yi = 0; yi < ysize; yi++)
        {
          double y = yarr[yi];

          zi = gsl_interp2d_idx(interp, xi, yi);
          test_single_low_level(&gsl_interp2d_eval, &gsl_interp2d_eval_e,
                                interp, xarr, yarr, zarr, x, y,
                                xa, ya, zarr, zi);
          test_single_low_level(&gsl_interp2d_eval_extrap,
                                &gsl_interp2d_eval_extrap_e, interp,
                                xarr, yarr, zarr, x, y, xa, ya, zarr, zi);
          test_single_high_level(&gsl_spline2d_eval, &gsl_spline2d_eval_e,
                                 interp_s, x, y, xa, ya, zarr, zi);
          test_single_high_level(&gsl_spline2d_eval_extrap, &gsl_spline2d_eval_extrap_e,
                                 interp_s, x, y, xa, ya, zarr, zi);
        }
    }

  /* Then check additional points provided */
  for (i = 0; i < test_size; i++)
    {
      double x = xval[i];
      double y = yval[i];

      test_single_low_level(&gsl_interp2d_eval,         &gsl_interp2d_eval_e,          interp, xarr, yarr, zarr, x, y, xa, ya, zval, i);
      test_single_low_level(&gsl_interp2d_eval_deriv_x, &gsl_interp2d_eval_deriv_x_e,  interp, xarr, yarr, zarr, x, y, xa, ya, zxval, i);
      test_single_low_level(&gsl_interp2d_eval_deriv_y, &gsl_interp2d_eval_deriv_y_e,  interp, xarr, yarr, zarr, x, y, xa, ya, zyval, i);
      test_single_low_level(&gsl_interp2d_eval_deriv_xx,&gsl_interp2d_eval_deriv_xx_e, interp, xarr, yarr, zarr, x, y, xa, ya, zxxval, i);
      test_single_low_level(&gsl_interp2d_eval_deriv_yy,&gsl_interp2d_eval_deriv_yy_e, interp, xarr, yarr, zarr, x, y, xa, ya, zyyval, i);
      test_single_low_level(&gsl_interp2d_eval_deriv_xy,&gsl_interp2d_eval_deriv_xy_e, interp, xarr, yarr, zarr, x, y, xa, ya, zxyval, i);

      test_single_high_level(&gsl_spline2d_eval,         &gsl_spline2d_eval_e,          interp_s, x, y, xa, ya, zval, i);
      test_single_high_level(&gsl_spline2d_eval_extrap,  &gsl_spline2d_eval_extrap_e,          interp_s, x, y, xa, ya, zval, i);
      test_single_high_level(&gsl_spline2d_eval_deriv_x, &gsl_spline2d_eval_deriv_x_e,  interp_s, x, y, xa, ya, zxval, i);
      test_single_high_level(&gsl_spline2d_eval_deriv_y, &gsl_spline2d_eval_deriv_y_e,  interp_s, x, y, xa, ya, zyval, i);
      test_single_high_level(&gsl_spline2d_eval_deriv_xx,&gsl_spline2d_eval_deriv_xx_e, interp_s, x, y, xa, ya, zxxval, i);
      test_single_high_level(&gsl_spline2d_eval_deriv_yy,&gsl_spline2d_eval_deriv_yy_e, interp_s, x, y, xa, ya, zyyval, i);
      test_single_high_level(&gsl_spline2d_eval_deriv_xy,&gsl_spline2d_eval_deriv_xy_e, interp_s, x, y, xa, ya, zxyval, i);

      test_single_low_level(&gsl_interp2d_eval_extrap, &gsl_interp2d_eval_extrap_e, interp, xarr, yarr, zarr, x, y, xa, ya, zval, i);
    }

  gsl_interp_accel_free(xa);
  gsl_interp_accel_free(ya);
  gsl_interp2d_free(interp);
  gsl_spline2d_free(interp_s);

  return status;
}

/*
 * Tests bilinear interpolation using a symmetric function, f(x,y)==f(y,x),
 * and diagonal interpolation points (x,y) where x==y. If these tests don't
 * pass, something is seriously broken.
 */
static int
test_bilinear_symmetric()
{
  int status;
  double xarr[] = {0.0, 1.0, 2.0, 3.0};
  double yarr[] = {0.0, 1.0, 2.0, 3.0};
  double zarr[] = {1.0, 1.1, 1.2, 1.3,
                   1.1, 1.2, 1.3, 1.4,
                   1.2, 1.3, 1.4, 1.5,
                   1.3, 1.4, 1.5, 1.6};
  double xval[] = {0.0, 0.5, 1.0, 1.5, 2.5, 3.0};
  double yval[] = {0.0, 0.5, 1.0, 1.5, 2.5, 3.0};
  double zval[] = {1.0, 1.1, 1.2, 1.3, 1.5, 1.6};
  size_t xsize = sizeof(xarr) / sizeof(xarr[0]);
  size_t ysize = sizeof(yarr) / sizeof(yarr[0]);
  size_t test_size = sizeof(xval) / sizeof(xval[0]);

  status = test_interp2d(xarr, yarr, zarr, xsize, ysize, xval, yval, zval,
                         NULL, NULL, NULL, NULL, NULL, test_size,
                         gsl_interp2d_bilinear);
  gsl_test(status, "bilinear interpolation with symmetric values");

  return status;
}

/*
 * Tests bilinear interpolation with an asymmetric function, f(x,y)!=f(y,x),
 * and off-diagonal interpolation points (x,y) where x and y may or may not
 * be equal.
 */
static int
test_bilinear_asymmetric_z()
{
  int status;
  double xarr[] = {0.0, 1.0, 2.0, 3.0};
  double yarr[] = {0.0, 1.0, 2.0, 3.0};
  double zarr[] = {1.0, 1.1, 1.2, 1.4,
                   1.3, 1.4, 1.5, 1.7,
                   1.5, 1.6, 1.7, 1.9,
                   1.6, 1.9, 2.2, 2.3};
  double xval[] = { 0.0, 0.5, 1.0, 1.5,  2.5, 3.0,
                    1.3954, 1.6476, 0.824957,
                    2.41108,  2.98619, 1.36485 };
  double yval[] = {0.0, 0.5, 1.0, 1.5,  2.5, 3.0,
                   0.265371, 2.13849, 1.62114,
                   1.22198, 0.724681, 0.0596087 };

  /* results computed using Mathematica 9.0.1.0 */
  double zval[] = {1.0, 1.2, 1.4, 1.55, 2.025, 2.3,
                   1.2191513, 1.7242442248, 1.5067237,
                   1.626612, 1.6146423, 1.15436761};
  size_t xsize = sizeof(xarr) / sizeof(xarr[0]);
  size_t ysize = sizeof(yarr) / sizeof(yarr[0]);
  size_t test_size = sizeof(xval) / sizeof(xval[0]);

  status = test_interp2d(xarr, yarr, zarr, xsize, ysize, xval, yval, zval,
                         NULL, NULL, NULL, NULL, NULL, test_size,
                         gsl_interp2d_bilinear);
  gsl_test(status, "bilinear interpolation with asymmetric z values");

  return status;
}

static int
test_bicubic()
{
  int status;
  double xarr[] = {0.0, 1.0, 2.0, 3.0};
  double yarr[] = {0.0, 1.0, 2.0, 3.0};
  double zarr[] = {1.0, 1.1, 1.2, 1.3,
                   1.1, 1.2, 1.3, 1.4,
                   1.2, 1.3, 1.4, 1.5,
                   1.3, 1.4, 1.5, 1.6};
  double xval[] = {1.0, 1.5, 2.0};
  double yval[] = {1.0, 1.5, 2.0};
  double zval[] = {1.2, 1.3, 1.4};
  size_t xsize = sizeof(xarr) / sizeof(xarr[0]);
  size_t ysize = sizeof(yarr) / sizeof(yarr[0]);
  size_t test_size = sizeof(xval) / sizeof(xval[0]);

  status = test_interp2d(xarr, yarr, zarr, xsize, ysize, xval, yval, zval,
                         NULL, NULL, NULL, NULL, NULL, test_size,
                         gsl_interp2d_bicubic);
  gsl_test(status, "bicubic interpolation on linear function");

  return status;
}

static int
test_bicubic_nonlinear()
{
  int status;
  double xarr[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
  double yarr[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
  /* least common multiple of x and y */
  double zarr[] = { 1,  2,  3,  4,  5,  6,  7,  8,
                    2,  2,  6,  4, 10,  6, 14,  8,
                    3,  6,  3, 12, 15,  6, 21, 24,
                    4,  4, 12,  4, 20, 12, 28,  8,
                    5, 10, 15, 20,  5, 30, 35, 40,
                    6,  6,  6, 12, 30,  6, 42, 24,
                    7, 14, 21, 28, 35, 42,  7, 56,
                    8,  8, 24,  8, 40, 24, 56,  8};
  double xval[] = {1.4, 2.3, 4.7, 3.3, 7.5, 6.6, 5.1};
  double yval[] = {1.0, 1.8, 1.9, 2.5, 2.7, 4.1, 3.3};

  /* results computed using GSL 1D cubic interpolation twice */
  double zval[] = { 1.4, 3.11183531264736, 8.27114315792559, 5.03218982537718,
                    22.13230634702637, 23.63206834997871, 17.28553080971182 };
  size_t xsize = sizeof(xarr) / sizeof(xarr[0]);
  size_t ysize = sizeof(yarr) / sizeof(yarr[0]);
  size_t test_size = sizeof(xval) / sizeof(xval[0]);

  status = test_interp2d(xarr, yarr, zarr, xsize, ysize, xval, yval, zval,
                         NULL, NULL, NULL, NULL, NULL, test_size,
                         gsl_interp2d_bicubic);
  gsl_test(status, "bicubic interpolation on nonlinear symmetric function");

  return status;
}

/* This function contributed by Andrew W. Steiner <awsteiner@gmail.com> */
static int
test_bicubic_nonlinear_nonsq()
{
  int status;
  double xarr[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
  double yarr[] = {1.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0};
  double zarr[] = { 1,  2,  3,  4,  5,  6,  7,  8, 9, 10,
                    2,  2,  6,  4, 10,  6, 14,  8, 11, 12,
                    3,  6,  3, 12, 15,  6, 21, 24, 13, 14,
                    4,  4, 12,  4, 20, 12, 28,  8, 15, 16,
                    5, 10, 15, 20,  5, 30, 35, 40, 17, 18,
                    6,  6,  6, 12, 30,  6, 42, 24, 19, 20,
                    7, 14, 21, 28, 35, 42,  7, 56, 21, 22,
                    8,  8, 24,  8, 40, 24, 56,  8, 23, 24};
  double xval[] = {1.4, 2.3, 9.7, 3.3, 9.5, 6.6, 5.1};
  double yval[] = {1.0, 1.8, 1.9, 2.5, 2.7, 4.1, 3.3};

  /* results computed using GSL 1D cubic interpolation twice */
  double zval[] = { 1.4, 2.46782030941187003, 10.7717721621846465,
                    4.80725067958096375, 11.6747032398627297,
                    11.2619968682970111, 9.00168877916872567};
  size_t xsize = sizeof(xarr) / sizeof(xarr[0]);
  size_t ysize = sizeof(yarr) / sizeof(yarr[0]);
  size_t test_size = sizeof(xval) / sizeof(xval[0]);

  status = test_interp2d(xarr, yarr, zarr, xsize, ysize, xval, yval, zval,
                         NULL, NULL, NULL, NULL, NULL, test_size,
                         gsl_interp2d_bicubic);
  gsl_test(status, "bicubic interpolation on nonlinear symmetric function");

  return status;
}

/* runs all the tests */
int
test_interp2d_main(void)
{
  int status = 0;

  status += test_bilinear_symmetric();
  status += test_bilinear_asymmetric_z();
  status += test_bicubic();
  status += test_bicubic_nonlinear();
  status += test_bicubic_nonlinear_nonsq();

  return status;
}

max_stars_repo_name: igormcoelho/optstats

max_issues_repo_name: igormcoelho/optstats

max_forks_repo_name: igormcoelho/optstats


text:
/***
 * Copyright 2020 The Katla Authors
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#ifndef KATLA_SQLITE_DATABASE_H
#define KATLA_SQLITE_DATABASE_H

#include "katla/core/core.h"
#include "katla/core/error.h"

#include "sqlite3.h"

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

namespace katla {

struct SqliteTableData {
    int nrOfColumns {};
    std::vector<std::string> columnNames;
    std::vector<std::string> data;
};

struct SqliteQueryResult {
    std::optional<SqliteTableData> queryResult;
};

class SqliteDatabase {
public:
    SqliteDatabase();
    virtual ~SqliteDatabase();

    outcome::result<void, Error> init();

    outcome::result<void, Error> open();
    outcome::result<void, Error> close();

    outcome::result<void, Error> create(std::string path);

    outcome::result<SqliteQueryResult, Error> exec(std::string sql);

    outcome::result<SqliteQueryResult, Error> insert(std::string table,  gsl::span<std::pair<std::string, std::string>> values);

    void setPath(std::string path) {
        m_path = path;
    }

    bool isOpen () {
        return m_handle != nullptr;
    }

private:
    sqlite3* m_handle {};
    std::string m_path {};
};

}

#endif

max_stars_repo_name: plok/katla

max_issues_repo_name: plok/katla

max_forks_repo_name: plok/katla


text:
#pragma once
#include "common/to_do.h"
#include "open_file.h"
#include "trees/regular_file.h"
#include <gsl/span>

namespace dogbox::tree
{
    enum class read_caching
    {
        none,
        one_piece
    };

    inline std::ostream &operator<<(std::ostream &out, read_caching const printed)
    {
        switch (printed)
        {
        case read_caching::none:
            return out << "none";
        case read_caching::one_piece:
            return out << "one_piece";
        }
        TO_DO();
    }

    constexpr dogbox::tree::read_caching all_read_caching_modes[] = {
        dogbox::tree::read_caching::none, dogbox::tree::read_caching::one_piece};

    size_t read_file(open_file &file, sqlite3 &database, regular_file::length_type const offset,
                     gsl::span<std::byte> const into, read_caching const caching);
}

max_stars_repo_name: TyRoXx/dogbox

max_issues_repo_name: TyRoXx/dogbox

max_forks_repo_name: TyRoXx/dogbox


text:
/* cdf/tdistinv.c
 *
 * Copyright (C) 2007, 2010 Brian Gough
 * Copyright (C) 2002 Jason H. Stover.
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 */

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

#include <stdio.h>

static double
inv_cornish_fisher (double z, double nu)
{
  double a = 1 / (nu - 0.5);
  double b = 48.0 / (a * a);

  double cf1 = z * (3 + z * z);
  double cf2 = z * (945 + z * z * (360 + z * z * (63 + z * z * 4)));

  double y = z - cf1 / b + cf2 / (10 * b * b);

  double t = GSL_SIGN (z) * sqrt (nu * expm1 (a * y * y));

  return t;
}

double
gsl_cdf_tdist_Pinv (const double P, const double nu)
{
  double x, ptail;

  if (P == 1.0)
    {
      return GSL_POSINF;
    }
  else if (P == 0.0)
    {
      return GSL_NEGINF;
    }

  if (nu == 1.0)
    {
      x = tan (M_PI * (P - 0.5));
      return x;
    }
  else if (nu == 2.0)
    {
      x = (2 * P - 1) / sqrt (2 * P * (1 - P));
      return x;
    }

  ptail = (P < 0.5) ? P : 1 - P;

  if (sqrt (M_PI * nu / 2) * ptail > pow (0.05, nu / 2))
    {
      double xg = gsl_cdf_ugaussian_Pinv (P);
      x = inv_cornish_fisher (xg, nu);
    }
  else
    {
      /* Use an asymptotic expansion of the tail of integral */

      double beta = gsl_sf_beta (0.5, nu / 2);

      if (P < 0.5)
        {
          x = -sqrt (nu) * pow (beta * nu * P, -1.0 / nu);
        }
      else
        {
          x = sqrt (nu) * pow (beta * nu * (1 - P), -1.0 / nu);
        }

      /* Correct nu -> nu/(1+nu/x^2) in the leading term to account
         for higher order terms. This avoids overestimating x, which
         makes the iteration unstable due to the rapidly decreasing
         tails of the distribution. */

      x /= sqrt (1 + nu / (x * x));
    }

  {
    double dP, phi;
    unsigned int n = 0;

  start:
    dP = P - gsl_cdf_tdist_P (x, nu);
    phi = gsl_ran_tdist_pdf (x, nu);

    if (dP == 0.0 || n++ > 32)
      goto end;

    {
      double lambda = dP / phi;
      double step0 = lambda;
      double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0);

      double step = step0;

      if (fabs (step1) < fabs (step0))
        {
          step += step1;
        }

      if (P > 0.5 && x + step < 0)
        x /= 2;
      else if (P < 0.5 && x + step > 0)
        x /= 2;
      else
        x += step;

      if (fabs (step) > 1e-10 * fabs (x))
        goto start;
    }

  end:
    if (fabs(dP) > GSL_SQRT_DBL_EPSILON * P)
      {
        GSL_ERROR_VAL("inverse failed to converge", GSL_EFAILED, GSL_NAN);
      }

    return x;
  }
}

double
gsl_cdf_tdist_Qinv (const double Q, const double nu)
{
  double x, qtail;

  if (Q == 0.0)
    {
      return GSL_POSINF;
    }
  else if (Q == 1.0)
    {
      return GSL_NEGINF;
    }

  if (nu == 1.0)
    {
      x = tan (M_PI * (0.5 - Q));
      return x;
    }
  else if (nu == 2.0)
    {
      x = (1 - 2 * Q) / sqrt (2 * Q * (1 - Q));
      return x;
    }

  qtail = (Q < 0.5) ? Q : 1 - Q;

  if (sqrt (M_PI * nu / 2) * qtail > pow (0.05, nu / 2))
    {
      double xg = gsl_cdf_ugaussian_Qinv (Q);
      x = inv_cornish_fisher (xg, nu);
    }
  else
    {
      /* Use an asymptotic expansion of the tail of integral */

      double beta = gsl_sf_beta (0.5, nu / 2);

      if (Q < 0.5)
        {
          x = sqrt (nu) * pow (beta * nu * Q, -1.0 / nu);
        }
      else
        {
          x = -sqrt (nu) * pow (beta * nu * (1 - Q), -1.0 / nu);
        }

      /* Correct nu -> nu/(1+nu/x^2) in the leading term to account
         for higher order terms. This avoids overestimating x, which
         makes the iteration unstable due to the rapidly decreasing
         tails of the distribution. */

      x /= sqrt (1 + nu / (x * x));
    }

  {
    double dQ, phi;
    unsigned int n = 0;

  start:
    dQ = Q - gsl_cdf_tdist_Q (x, nu);
    phi = gsl_ran_tdist_pdf (x, nu);

    if (dQ == 0.0 || n++ > 32)
      goto end;

    {
      double lambda = - dQ / phi;
      double step0 = lambda;
      double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0);

      double step = step0;

      if (fabs (step1) < fabs (step0))
        {
          step += step1;
        }

      if (Q < 0.5 && x + step < 0)
        x /= 2;
      else if (Q > 0.5 && x + step > 0)
        x /= 2;
      else
        x += step;

      if (fabs (step) > 1e-10 * fabs (x))
        goto start;
    }
  }

end:

  return x;
}


text:
#ifndef __CNN_CONV_H__
#define __CNN_CONV_H__

#include <cblas.h>
#include <string.h>

#include "cnn_macro.h"
#include "cnn_types.h"

#ifdef CNN_WITH_CUDA
#include <cublas_v2.h>
#include <cuda_runtime.h>

#include "cnn_init.h"

void cnn_map_gpu(float* dst, float* src, int* map, int len);
void cnn_map_inv_gpu(float* dst, float* src, int* map, int len);
#endif

static inline void cnn_conv_unroll_2d_valid(int* indexMap, int dstHeight,
                                            int dstWidth, int kSize,
                                            int srcHeight, int srcWidth,
                                            int srcCh)
{
    int __kMemSize = kSize * kSize;
    int __srcImSize = srcHeight * srcWidth;
    int __indexMapCols = __kMemSize * srcCh;

    for (int __h = 0; __h < dstHeight; __h++)
    {
        int __dstRowShift = __h * dstWidth;

        for (int __w = 0; __w < dstWidth; __w++)
        {
            int __indexMapRow = __dstRowShift + __w;
            int __indexMemBase = __indexMapRow * __indexMapCols;

            for (int __ch = 0; __ch < srcCh; __ch++)
            {
                int __indexMemShiftBase = __indexMemBase + __kMemSize * __ch;
                int __srcChShift = __ch * __srcImSize;

                for (int __convH = 0; __convH < kSize; __convH++)
                {
                    int __indexMemShift = __indexMemShiftBase + __convH * kSize;
                    int __srcShift = (__h + __convH) * srcWidth + __srcChShift;

                    for (int __convW = 0; __convW < kSize; __convW++)
                    {
                        indexMap[__indexMemShift + __convW] =
                            __srcShift + (__w + __convW);
                    }
                }
            }
        }
    }
}

static inline void cnn_conv_unroll_2d_same(int* indexMap, int dstHeight,
                                           int dstWidth, int kSize,
                                           int srcHeight, int srcWidth,
                                           int srcCh)
{
    int __kMemSize = kSize * kSize;
    int __srcImSize = srcHeight * srcWidth;
    int __indexMapCols = __kMemSize * srcCh;

    int __convHBase = -kSize / 2;
    int __convWBase = -kSize / 2;

    for (int __h = 0; __h < dstHeight; __h++)
    {
        int __dstRowShift = __h * dstWidth;

        for (int __w = 0; __w < dstWidth; __w++)
        {
            int __indexMapRow = __dstRowShift + __w;
            int __indexMemBase = __indexMapRow * __indexMapCols;

            for (int __ch = 0; __ch < srcCh; __ch++)
            {
                int __indexMemShiftBase = __indexMemBase + __kMemSize * __ch;
                int __srcChShift = __ch * __srcImSize;

                for (int __convH = 0; __convH < kSize; __convH++)
                {
                    int __indexMemShift = __indexMemShiftBase + __convH * kSize;
                    int __convHIndex = __h + __convH + __convHBase;

                    if (__convHIndex >= 0 && __convHIndex < srcHeight)
                    {
                        int __srcShift = __convHIndex * srcWidth + __srcChShift;

                        for (int __convW = 0; __convW < kSize; __convW++)
                        {
                            int __convWIndex = __w + __convW + __convWBase;
                            if (__convWIndex >= 0 && __convWIndex < srcWidth)
                            {
                                int __tmpIndex = __srcShift + __convWIndex;

                                indexMap[__indexMemShift + __convW] =
                                    __tmpIndex;
                            }
                        }
                    }
                }
            }
        }
    }
}

static inline void cnn_conv_2d(float* dst, int dstHeight, int dstWidth,
                               float* kernel, int kSize, int chIn, int chOut,
                               float* src, int srcHeight, int srcWidth)
{
    int __kMemSize = kSize * kSize;
    int __filterSize = chIn * __kMemSize;
    int __dstImSize = dstHeight * dstWidth;
    int __srcImSize = srcHeight * srcWidth;

    for (int __chOut = 0; __chOut < chOut; __chOut++)
    {
        int __filterShift = __chOut * __filterSize;
        int __dstChShift = __chOut * __dstImSize;

        for (int __chIn = 0; __chIn < chIn; __chIn++)
        {
            int __kShiftBase = __chIn * __kMemSize + __filterShift;
            int __srcChShift = __chIn * __srcImSize;

            for (int __h = 0; __h < dstHeight; __h++)
            {
                int __dstShift = __h * dstWidth + __dstChShift;
                for (int __w = 0; __w < dstWidth; __w++)
                {
                    float __conv = 0;
                    for (int __convH = 0; __convH < kSize; __convH++)
                    {
                        int __kShift = __convH * kSize + __kShiftBase;
                        int __srcShift =
                            (__h + __convH) * srcWidth + __srcChShift;

                        for (int __convW = 0; __convW < kSize; __convW++)
                        {
                            __conv += kernel[__kShift + __convW] *
                                      src[__srcShift + (__w + __convW)];
                        }
                    }

                    dst[__dstShift + __w] += __conv;
                }
            }
        }
    }
}

static inline void cnn_conv_2d_grad(float* srcGrad, int srcHeight, int srcWidth,
                                    float* kernel, int kSize, int srcCh,
                                    int lCh, float* lGrad, int lHeight,
                                    int lWidth)
{
    int __kMemSize = kSize * kSize;
    int __filterSize = srcCh * __kMemSize;
    int __lImSize = lHeight * lWidth;
    int __srcImSize = srcHeight * srcWidth;

    for (int __lCh = 0; __lCh < lCh; __lCh++)
    {
        int __filterShift = __lCh * __filterSize;
        int __lChShift = __lCh * __lImSize;

        for (int __srcCh = 0; __srcCh < srcCh; __srcCh++)
        {
            int __srcChShift = __srcCh * __srcImSize;
            int __kShiftBase = __srcCh * __kMemSize + __filterShift;

            for (int __h = 0; __h < lHeight; __h++)
            {
                int __lShift = __h * lHeight + __lChShift;
                for (int __w = 0; __w < lWidth; __w++)
                {
                    for (int __convH = 0; __convH < kSize; __convH++)
                    {
                        int __kShift = __convH * kSize + __kShiftBase;
                        int __srcShift =
                            (__h + __convH) * srcWidth + __srcChShift;

                        for (int __convW = 0; __convW < kSize; __convW++)
                        {
                            srcGrad[__srcShift + (__w + __convW)] +=
                                lGrad[__lShift + __w] *
                                kernel[__kShift + __convW];
                        }
                    }
                }
            }
        }
    }
}

static inline void cnn_conv_2d_kernel_grad(float* lGrad, int lHeight,
                                           int lWidth, float* kGrad, int kSize,
                                           int lCh, int srcCh, float* src,
                                           int srcHeight, int srcWidth)
{
    int __kMemSize = kSize * kSize;
    int __filterSize = srcCh * __kMemSize;
    int __lImSize = lHeight * lWidth;
    int __srcImSize = srcHeight * srcWidth;

    for (int __lCh = 0; __lCh < lCh; __lCh++)
    {
        int __filterShift = __lCh * __filterSize;
        int __lChShift = __lCh * __lImSize;

        for (int __srcCh = 0; __srcCh < srcCh; __srcCh++)
        {
            int __srcChShift = __srcCh * __srcImSize;
            int __kShiftBase = __srcCh * __kMemSize + __filterShift;

            for (int __h = 0; __h < lHeight; __h++)
            {
                int __lShift = __h * lWidth + __lChShift;
                for (int __w = 0; __w < lWidth; __w++)
                {
                    for (int __convH = 0; __convH < kSize; __convH++)
                    {
                        int __kShift = __convH * kSize + __kShiftBase;
                        int __srcShift =
                            (__h + __convH) * srcWidth + __srcChShift;

                        for (int __convW = 0; __convW < kSize; __convW++)
                        {
                            kGrad[__kShift + __convW] +=
                                lGrad[__lShift + __w] *
                                src[__srcShift + (__w + __convW)];
                        }
                    }
                }
            }
        }
    }
}

static inline void cnn_forward_conv(union CNN_LAYER* layerRef,
                                    struct CNN_CONFIG* cfgRef, int layerIndex)
{
#ifdef CNN_WITH_CUDA
    float alpha = 1.0;
    float beta = 0.0;

    struct CNN_LAYER_CONV* layerPtr = &layerRef[layerIndex].conv;

    struct CNN_MAT* outData = &layerPtr->outMat.data;
    struct CNN_MAT* preOutData = &layerRef[layerIndex - 1].outMat.data;

    cnn_assert_cudnn(cudnnConvolutionForward(
        cnnInit.cudnnHandle,                        //
        &alpha,                                     //
        layerPtr->srcTen, preOutData->mat,          //
        layerPtr->kernelTen, layerPtr->kernel.mat,  //
        layerPtr->convDesc, layerPtr->convAlgoFW, cnnInit.wsData,
        cnnInit.wsSize,  //
        &beta,           //
        layerPtr->dstTen, outData->mat));

#if defined(CNN_CONV_BIAS_FILTER)
    beta = 1.0;
    cnn_assert_cudnn(cudnnAddTensor(cnnInit.cudnnHandle,                    //
                                    &alpha,                                 //
                                    layerPtr->biasTen, layerPtr->bias.mat,  //
                                    &beta,                                  //
                                    layerPtr->dstTen, outData->mat));
#endif

#else
    // Cache
    int mapRows =
        layerRef[layerIndex].outMat.width * layerRef[layerIndex].outMat.height;
    int mapCols = layerRef[layerIndex - 1].outMat.channel *
                  cfgRef->layerCfg[layerIndex].conv.size *
                  cfgRef->layerCfg[layerIndex].conv.size;
    int mapSize = mapRows * mapCols;
    int* indexMap = layerRef[layerIndex].conv.indexMap;

    int chOut = layerRef[layerIndex].outMat.channel;
    float* kernel = layerRef[layerIndex].conv.kernel.mat;

    // Clear outputs
    // memset(layerRef[layerIndex].outMat.data.mat, 0,
    //       sizeof(float) * layerRef[layerIndex].outMat.data.rows *
    //           layerRef[layerIndex].outMat.data.cols);

    for (int j = 0; j < cfgRef->batch; j++)
    {
        int srcShift = j * layerRef[layerIndex - 1].outMat.data.cols;
        int dstShift = j * layerRef[layerIndex].outMat.data.cols;
        int mapShift = j * mapSize;

        float* srcPtr = layerRef[layerIndex - 1].outMat.data.mat + srcShift;
        float* dstPtr = layerRef[layerIndex].outMat.data.mat + dstShift;
        float* mapPtr = layerRef[layerIndex].conv.unroll.mat + mapShift;

        for (int k = 0; k < mapSize; k++)
        {
            int tmpIndex = indexMap[k];
            if (tmpIndex >= 0)
            {
                mapPtr[k] = srcPtr[tmpIndex];
            }
        }

        cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, chOut, mapRows,
                    mapCols, 1.0, kernel, mapCols, mapPtr, mapCols, 0.0, dstPtr,
                    mapRows);

        // Add bias
#if defined(CNN_CONV_BIAS_FILTER)
#ifdef DEBUG
#pragma message("cnn_forward_conv(): Enable convolution filter bias")
#endif
        for (int ch = 0; ch < chOut; ch++)
        {
            cblas_saxpy(
                mapRows, 1.0, &layerRef[layerIndex].conv.bias.mat[ch], 0,
                &layerRef[layerIndex].outMat.data.mat[dstShift + ch * mapRows],
                1);
        }
#elif defined(CNN_CONV_BIAS_LAYER)
#ifdef DEBUG
#pragma message("cnn_forward_conv(): Enable convolution layer bias")
#endif
        cblas_saxpy(layerRef[layerIndex].conv.bias.cols, 1.0,
                    layerRef[layerIndex].conv.bias.mat, 1,
                    &layerRef[layerIndex].outMat.data.mat[dstShift], 1);
#endif
    }
#endif
}

static inline void cnn_backward_conv(union CNN_LAYER* layerRef,
                                     struct CNN_CONFIG* cfgRef, int layerIndex)
{
#ifdef CNN_WITH_CUDA
    float alpha = 1.0;
    float beta = 1.0;

    struct CNN_LAYER_CONV* layerPtr = &layerRef[layerIndex].conv;

    struct CNN_MAT* outData = &layerPtr->outMat.data;
    struct CNN_MAT* preOutData = &layerRef[layerIndex - 1].outMat.data;

    cnn_assert_cudnn(cudnnConvolutionBackwardFilter(
        cnnInit.cudnnHandle,                //
        &alpha,                             //
        layerPtr->srcTen, preOutData->mat,  //
        layerPtr->dstTen, outData->grad,    //
        layerPtr->convDesc, layerPtr->convAlgoBWFilter, cnnInit.wsData,
        cnnInit.wsSize,  //
        &beta,           //
        layerPtr->kernelTen, layerPtr->kernel.grad));

    cnn_assert_cudnn(
        cudnnConvolutionBackwardBias(cnnInit.cudnnHandle,              //
                                     &alpha,                           //
                                     layerPtr->dstTen, outData->grad,  //
                                     &beta,                            //
                                     layerPtr->biasTen, layerPtr->bias.grad));

#else
    // Cache
    int mapRows =
        layerRef[layerIndex].outMat.width * layerRef[layerIndex].outMat.height;
    int mapCols = layerRef[layerIndex - 1].outMat.channel *
                  cfgRef->layerCfg[layerIndex].conv.size *
                  cfgRef->layerCfg[layerIndex].conv.size;
    int mapSize = mapRows * mapCols;

    int* indexMap = layerRef[layerIndex].conv.indexMap;

    int chOut = layerRef[layerIndex].outMat.channel;
    float* kernel = layerRef[layerIndex].conv.kernel.mat;
    float* kGrad = layerRef[layerIndex].conv.kernel.grad;

    // Sum gradient
    for (int j = 0; j < cfgRef->batch; j++)
    {
        int gradShift = j * layerRef[layerIndex].outMat.data.cols;
        int mapShift = j * mapSize;

        float* gradPtr = &layerRef[layerIndex].outMat.data.grad[gradShift];
        float* mapPtr = &layerRef[layerIndex].conv.unroll.mat[mapShift];

        cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, chOut, mapCols,
                    mapRows, 1.0, gradPtr, mapRows, mapPtr, mapCols, 1.0, kGrad,
                    mapCols);

        // Sum bias gradient matrix
#if defined(CNN_CONV_BIAS_FILTER)
#ifdef DEBUG
#pragma message("cnn_forward_conv(): Enable convolution filter bias")
#endif
        for (int ch = 0; ch < chOut; ch++)
        {
            cblas_saxpy(mapRows, 1.0, &gradPtr[ch * mapRows], 1,
                        &layerRef[layerIndex].conv.bias.grad[ch], 0);
        }
#elif defined(CNN_CONV_BIAS_LAYER)
#ifdef DEBUG
#pragma message("cnn_forward_conv(): Enable convolution layer bias")
#endif
        cblas_saxpy(layerRef[layerIndex].conv.bias.cols, 1.0, gradPtr, 1,
                    layerRef[layerIndex].conv.bias.grad, 1);
#endif
    }
#endif

    // Find layer gradient
    if (layerIndex > 1)
    {
#ifdef CNN_WITH_CUDA
        cudaMemset(preOutData->grad, 0,
                   sizeof(float) * preOutData->rows * preOutData->cols);
        beta = 0.0;
        cnn_assert_cudnn(cudnnConvolutionBackwardData(
            cnnInit.cudnnHandle,                        //
            &alpha,                                     //
            layerPtr->kernelTen, layerPtr->kernel.mat,  //
            layerPtr->dstTen, outData->grad,            //
            layerPtr->convDesc, layerPtr->convAlgoBWGrad, cnnInit.wsData,
            cnnInit.wsSize,  //
            &beta,           //
            layerPtr->srcTen, preOutData->grad));
#else
        memset(layerRef[layerIndex - 1].outMat.data.grad, 0,
               sizeof(float) * layerRef[layerIndex - 1].outMat.data.rows *
                   layerRef[layerIndex - 1].outMat.data.cols);
        for (int j = 0; j < cfgRef->batch; j++)
        {
            int gradShift = j * layerRef[layerIndex].outMat.data.cols;
            int preGradShift = j * layerRef[layerIndex - 1].outMat.data.cols;
            int mapShift = j * mapSize;

            float* gradPtr = layerRef[layerIndex].outMat.data.grad + gradShift;
            float* preGradPtr =
                layerRef[layerIndex - 1].outMat.data.grad + preGradShift;
            float* mapPtr = layerRef[layerIndex].conv.unroll.grad + mapShift;

            cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, mapRows,
                        mapCols, chOut, 1.0, gradPtr, mapRows, kernel, mapCols,
                        0.0, mapPtr, mapCols);

            for (int i = 0; i < mapSize; i++)
            {
                int tmpIndex = indexMap[i];
                if (tmpIndex >= 0)
                {
                    preGradPtr[tmpIndex] += mapPtr[i];
                }
            }
        }
#endif
    }
}

#endif

max_stars_repo_name: jamesljlster/cnn

max_issues_repo_name: jamesljlster/cnn

max_forks_repo_name: jamesljlster/cnn