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

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

/*! \file second_order.c
 *  \brief produce actual ICs from special 2nd order lpt ICs
 */

double F1_Omega(double a)
{
  double omega_a;

  omega_a = All.Omega0 / (All.Omega0 + a * (1 - All.Omega0 - All.OmegaLambda) + a * a * a * All.OmegaLambda);

  return pow(omega_a, 5.0 / 9);
}

double F2_Omega(double a)
{
  double omega_a;

  omega_a = All.Omega0 / (All.Omega0 + a * (1 - All.Omega0 - All.OmegaLambda) + a * a * a * All.OmegaLambda);

  return 2.0 * pow(omega_a, 6.0 / 11);
}

void second_order_ics(void)
{
  int i, j;
  double a, hubble, fac1, fac2, fac3;
  double ax, ay, az;

  if(ThisTask == 0)
    {
      printf("\nNow producing ICs based on 2nd-order LPT\n\n");
      fflush(stdout);
    }

  a = All.TimeBegin;

  hubble = hubble_function(a);

  fac1 = 1.0 / (a * a * hubble * F1_Omega(a));	/* this factor converts the Zeldovich peculiar velocity
						   (expressed as  a^2*dX/dt  to comoving displacement */
  fac2 = header.lpt_scalingfactor;

  fac3 = fac2 * a * a * hubble * F2_Omega(a);

  if(ThisTask == 0)
    printf("fac1=%g  fac2=%g  fac3=%g\n", fac1, fac2, fac3);

  for(i = 0; i < NumPart; i++)
    {
      for(j = 0; j < 3; j++)
	P[i].Pos[j] += fac1 * P[i].Vel[j];	/* Zeldovich displacement */

#ifdef PMGRID
      for(j = 0; j < 3; j++)
	P[i].Pos[j] += fac2 * (P[i].GravPM[j] + P[i].g.GravAccel[j]);	/* second order lpt displacement */

      for(j = 0; j < 3; j++)
	P[i].Vel[j] += fac3 * (P[i].GravPM[j] + P[i].g.GravAccel[j]);	/* second order lpt velocity correction */
#else
      for(j = 0; j < 3; j++)
	P[i].Pos[j] += fac2 * (P[i].g.GravAccel[j]);	/* second order lpt displacement */

      for(j = 0; j < 3; j++)
	P[i].Vel[j] += fac3 * (P[i].g.GravAccel[j]);	/* second order lpt velocity correction */
#endif
    }

  /* now set the masses correctly */

  for(i = 0; i < NumPart; i++)
    {
      P[i].Mass = P[i].OldAcc;

      if(All.MassTable[P[i].Type] != 0)
	P[i].Mass = All.MassTable[P[i].Type];
    }

  if(All.ComovingIntegrationOn)
    if(All.PeriodicBoundariesOn == 1)
      check_omega();		/* check again whether we have plausible mass density */

  /* recompute long range force */

  All.DoDynamicUpdate = 1;
  domain_Decomposition();	/* redo domain decomposition because particles may have shifted */

  CPU_Step[CPU_MISC] += measure_time();

#ifdef PMGRID
  long_range_force();
#endif

  CPU_Step[CPU_MESH] += measure_time();

  gravity_tree();

  for(i = 0; i < NumPart; i++)
    {
#ifdef PMGRID
      ax = (P[i].g.GravAccel[0] + P[i].GravPM[0]) / All.G;
      ay = (P[i].g.GravAccel[1] + P[i].GravPM[1]) / All.G;
      az = (P[i].g.GravAccel[2] + P[i].GravPM[2]) / All.G;
#else
      ax = P[i].g.GravAccel[0] / All.G;
      ay = P[i].g.GravAccel[1] / All.G;
      az = P[i].g.GravAccel[2] / All.G;
#endif
      P[i].OldAcc = sqrt(ax * ax + ay * ay + az * az);
    }

  if(All.TypeOfOpeningCriterion == 1)
    {
      All.ErrTolTheta = 0;	/* This will switch to the relative opening criterion for the following force computations */
      gravity_tree();
    }

}

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
#include <stdio.h>
#include <stdlib.h>
#include <lapacke.h>

void inverse(double* A, lapack_int N)
{
    lapack_int *IPIV = (lapack_int *) malloc(N * sizeof(lapack_int));
    int INFO;

    /*
        DGETRF computes an LU factorization of a general M-by-N matrix A
        using partial pivoting with row interchanges.

        The factorization has the form
            A = P * L * U
        where P is a permutation matrix, L is lower triangular with unit
        diagonal elements (lower trapezoidal if m > n), and U is upper
        triangular (upper trapezoidal if m < n).

        This is the right-looking Level 3 BLAS version of the algorithm.

        lapack_int LAPACKE_dgetrf( int matrix_layout, lapack_int m, lapack_int n,
                        double* a, lapack_int lda, lapack_int* ipiv );

        INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
                has been completed, but the factor U is exactly
                singular, and division by zero will occur if it is used
                to solve a system of equations.
    */

    INFO = LAPACKE_dgetrf(LAPACK_ROW_MAJOR,N,N,A,N,IPIV);
    printf("%d ", INFO);

    //DGETRI computes the inverse of a matrix using the LU factorization
    //computed by DGETRF.

    //This method inverts U and then computes inv(A) by solving the system
    //inv(A)*L = inv(U) for inv(A).
    //lapack_int LAPACKE_dgetri( int matrix_layout, lapack_int n, double* a,
    //          lapack_int lda, const lapack_int* ipiv );
    INFO = LAPACKE_dgetri(LAPACK_ROW_MAJOR,N,A,N,IPIV);
    printf("%d \n", INFO);

    free(IPIV);
}

int main(){

    double A [2*2] = {
        1,2,
        3,4
    };

    inverse(A, 2);

    printf("%f %f\n", A[0], A[1]);
    printf("%f %f\n", A[2], A[3]);

    return 0;
}

ext: c

max_stars_repo_name: agnusfeec/lc

max_issues_repo_name: agnusfeec/lc

max_forks_repo_name: agnusfeec/lc


text:
/* poly/zsolve_cubic.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.
 */

/* zsolve_cubic.c - finds the complex roots of x^3 + a x^2 + b x + c = 0 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_poly.h>

#define SWAP(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0)

int
gsl_poly_complex_solve_cubic (double a, double b, double c,
                              gsl_complex *z0, gsl_complex *z1,
                              gsl_complex *z2)
{
  double q = (a * a - 3 * b);
  double r = (2 * a * a * a - 9 * a * b + 27 * c);

  double Q = q / 9;
  double R = r / 54;

  double Q3 = Q * Q * Q;
  double R2 = R * R;

  double CR2 = 729 * r * r;
  double CQ3 = 2916 * q * q * q;

  if (R == 0 && Q == 0)
    {
      GSL_REAL (*z0) = -a / 3;
      GSL_IMAG (*z0) = 0;
      GSL_REAL (*z1) = -a / 3;
      GSL_IMAG (*z1) = 0;
      GSL_REAL (*z2) = -a / 3;
      GSL_IMAG (*z2) = 0;
      return 3;
    }
  else if (CR2 == CQ3)
    {
      /* this test is actually R2 == Q3, written in a form suitable
         for exact computation with integers */

      /* Due to finite precision some double roots may be missed, and
         will be considered to be a pair of complex roots z = x +/-
         epsilon i close to the real axis. */

      double sqrtQ = sqrt (Q);

      if (R > 0)
	{
	  GSL_REAL (*z0) = -2 * sqrtQ - a / 3;
	  GSL_IMAG (*z0) = 0;
	  GSL_REAL (*z1) = sqrtQ - a / 3;
	  GSL_IMAG (*z1) = 0;
	  GSL_REAL (*z2) = sqrtQ - a / 3;
	  GSL_IMAG (*z2) = 0;
	}
      else
	{
	  GSL_REAL (*z0) = -sqrtQ - a / 3;
	  GSL_IMAG (*z0) = 0;
	  GSL_REAL (*z1) = -sqrtQ - a / 3;
	  GSL_IMAG (*z1) = 0;
	  GSL_REAL (*z2) = 2 * sqrtQ - a / 3;
	  GSL_IMAG (*z2) = 0;
	}
      return 3;
    }
  else if (CR2 < CQ3)  /* equivalent to R2 < Q3 */
    {
      double sqrtQ = sqrt (Q);
      double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
      double theta = acos (R / sqrtQ3);
      double norm = -2 * sqrtQ;
      double r0 = norm * cos (theta / 3) - a / 3;
      double r1 = norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3;
      double r2 = norm * cos ((theta - 2.0 * M_PI) / 3) - a / 3;

      /* Sort r0, r1, r2 into increasing order */

      if (r0 > r1)
	SWAP (r0, r1);

      if (r1 > r2)
	{
	  SWAP (r1, r2);

	  if (r0 > r1)
	    SWAP (r0, r1);
	}

      GSL_REAL (*z0) = r0;
      GSL_IMAG (*z0) = 0;

      GSL_REAL (*z1) = r1;
      GSL_IMAG (*z1) = 0;

      GSL_REAL (*z2) = r2;
      GSL_IMAG (*z2) = 0;

      return 3;
    }
  else
    {
      double sgnR = (R >= 0 ? 1 : -1);
      double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0 / 3.0);
      double B = Q / A;

      if (A + B < 0)
	{
	  GSL_REAL (*z0) = A + B - a / 3;
	  GSL_IMAG (*z0) = 0;

	  GSL_REAL (*z1) = -0.5 * (A + B) - a / 3;
	  GSL_IMAG (*z1) = -(sqrt (3.0) / 2.0) * fabs(A - B);

	  GSL_REAL (*z2) = -0.5 * (A + B) - a / 3;
	  GSL_IMAG (*z2) = (sqrt (3.0) / 2.0) * fabs(A - B);
	}
      else
	{
	  GSL_REAL (*z0) = -0.5 * (A + B) - a / 3;
	  GSL_IMAG (*z0) = -(sqrt (3.0) / 2.0) * fabs(A - B);

	  GSL_REAL (*z1) = -0.5 * (A + B) - a / 3;
	  GSL_IMAG (*z1) = (sqrt (3.0) / 2.0) * fabs(A - B);

	  GSL_REAL (*z2) = A + B - a / 3;
	  GSL_IMAG (*z2) = 0;
	}

      return 3;
    }
}


text:
// Copyright 2019, Winfield Chen and Lloyd T. Elliott.

#ifndef __STATISTICS_H_
#  define __STATISTICS_H_

#include <gsl/gsl_sf.h>
#include "matrix.h"

double tcdf1m(double t, double nu);
double log_tcdf1m(double t, double nu);
void regression(t_matrix g,
                t_matrix y,
                t_matrix yt,
                t_matrix obs,
                t_matrix denom,
                t_matrix beta,
                t_matrix se,
                t_matrix tstat,
                t_matrix pval,
                t_matrix b1,
                t_matrix w1,
                t_matrix w2);
#endif

ext: h

max_stars_repo_name: wnfldchen/agent

max_issues_repo_name: wnfldchen/agent

max_forks_repo_name: wnfldchen/agent


text:
#pragma once
/*
 * (C) Copyright 2020-2021 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.
 */
/*! \addtogroup ioda_cxx_variable
 *
 * @{
 * \file FillPolicy.h
 * \brief Default fill values for ioda files.
 */

#include <gsl/gsl-lite.hpp>
#include <memory>
#include <string>

#include "ioda/Variables/Fill.h"
#include "ioda/Exception.h"
#include "ioda/defs.h"

namespace ioda {

/// \brief This option describes the default fill values that will be used if the user does not
/// manually specify a fill value.
/// \ingroup ioda_cxx_variable
enum class FillValuePolicy {
  HDF5,    ///< Set all fill values to zero or null strings.
  NETCDF4  ///< Use NetCDF4 default fill values. This is the default option for ioda files.
};

/// \brief Holds the different default fill values used in ioda files produced by different
///   backends.
/// \ingroup ioda_cxx_variable
/// \details This matters for netCDF4 vs HDF5-produced files. They have different default
///   fill values.
namespace FillValuePolicies {
template <class T>
T HDF5_default() {
  return 0;
}
template <>
inline std::string HDF5_default<std::string>() {
  return std::string();
}

/// \ingroup ioda_cxx_variable
/// \see netcdf.h, starting around line 62, for these values
///   netcdf uses "ints" and "shorts", but these are all defined as fixed-width types.
template <class T>
T netCDF4_default() {
  return 0;
}
template <>
inline std::string netCDF4_default<std::string>() {
  return std::string();
}
template <>
inline signed char netCDF4_default<signed char>() {
  return static_cast<signed char>(-127);
}
template <>
inline char netCDF4_default<char>() {
  return static_cast<char>(0);
}
template <>
inline int16_t netCDF4_default<int16_t>() {
  return static_cast<int16_t>(-32767);
}
template <>
inline int32_t netCDF4_default<int32_t>() {
  return -2147483647;
}
template <>
inline float netCDF4_default<float>() {
  return 9.9692099683868690e+36f;
}
template <>
inline double netCDF4_default<double>() {
  return 9.9692099683868690e+36;
}
template <>
inline unsigned char netCDF4_default<unsigned char>() {
  return static_cast<unsigned char>(255);
}
template <>
inline uint16_t netCDF4_default<uint16_t>() {
  return static_cast<unsigned short>(65535);
}
template <>
inline uint32_t netCDF4_default<uint32_t>() {
  return 4294967295U;
}
template <>
inline int64_t netCDF4_default<int64_t>() {
  return -9223372036854775806LL;
}
template <>
inline uint64_t netCDF4_default<uint64_t>() {
  return 18446744073709551614ULL;
}

/// \brief Applies the fill value policy. This sets default fill values when fill values are not
///   already provided.
/// \ingroup ioda_cxx_variable
template <class T>
void applyFillValuePolicy(FillValuePolicy pol, detail::FillValueData_t& fvd) {
  if (fvd.set_) return;  // If already set, then do nothing.
  if (pol == FillValuePolicy::HDF5)
    detail::assignFillValue(fvd, HDF5_default<T>());
  else if (pol == FillValuePolicy::NETCDF4)
    detail::assignFillValue(fvd, netCDF4_default<T>());
  else
    throw Exception("Unsupported fill value policy.", ioda_Here());
}

}  // namespace FillValuePolicies
}  // namespace ioda

/// @}

max_stars_repo_name: Gibies/ioda

max_issues_repo_name: Gibies/ioda

max_forks_repo_name: Gibies/ioda


text:
/**
 * Subroutines to correct for the intra-pixel sensitivity
 * variation in e.g. NICMOS pixels. There are routines
 * to corect both, an extracted SPC's and a list
 * of PET pixels.
 *
 */
#include "fitsio.h"
#include <gsl/gsl_interp.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_multifit_nlin.h>
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "aXe_errors.h"
#include "spc_spc.h"
#include "spce_pathlength.h"
#include "fringe_conf.h"
#include "spc_wl_calib.h"
#include "trfit_utils.h"
#include "ipixcorr_utils.h"

#define MAX(x,y) (((x)>(y))?(x):(y))
#define MIN(x,y) (((x)<(y))?(x):(y))

/**
 * Function: get_intpix_corr
 * The function determines and returns a correction factor with
 * values stored in an interpolator at a certain wavelength.
 * To derive the result the dispersion solution is inverted,
 * then the trace position for the wavelength is determined.
 * The the fractional pixel value in y is computed and the
 * corresponding correction factor is returned.
 *
 * Parameters:
 * @param actbeam - beam to correct
 * @param cdisp   - the dispersion coefficients for the beam
 * @param ipcorr  - the correction values stored in an interpolator
 * @param lambda  - the wavelength to correct
 *
 * Returns:
 * @return cfactor - the correction factor
 */
double
get_intpix_corr(beam *actbeam, gsl_vector *cdisp, interpolator *ipcorr,
		float lambda)
{
  double tlength;

  double xvalue;
  double yfract;

  double cfactor=0.0;

  // get the trace length for the dispersion
  // solution and wavelength
  tlength = get_tlength_from_dispersion(cdisp, lambda);

  // compute the x-offset from the reference point
  // for this trace length value
  xvalue = get_xvalue_from_tlength(actbeam, tlength);

  // compute the fractional y-value
  yfract = get_yfract_for_xvalue(actbeam, xvalue);

  // get the correction factor
  cfactor = eval_interp(ipcorr, yfract);

  fprintf(stdout, "xdiff: %e, yfrac: %e, factor: %e ", xvalue, yfract, cfactor);

  //fprintf(stdout, "lambda: %f, tlength: %f, xval: %f, yval: %f, yfract: %f, cfactor: %f\n",
  // 	  lambda, tlength, actbeam->refpoint.x+xvalue, actbeam->refpoint.y+actbeam->spec_trace->func (xvalue, actbeam->spec_trace->data),yfract, cfactor);

  // return the correction factor
  return cfactor;
}

/**
 * Function: get_yfract_for_xvalue
 * The function computes the fractional y-value for a trace
 * position in a beam. The trace position is given as the
 * x-offset position with respect to the reference position
 *
 * Parameters:
 * @param actbeam - beam to correct
 * @param xvalue  - the x-value
 *
 * Returns:
 * @return yfract - the fractional y-value
 */
double
get_yfract_for_xvalue(const beam *actbeam, double xvalue)
{
  double dy;
  double yabs;
  double yfract;

  // compute the y-value relative to the reference point
  dy = actbeam->spec_trace->func (xvalue, actbeam->spec_trace->data);

  // compute the absolute y-value on the chip
  yabs = actbeam->refpoint.y + dy;

  // compute the fractional pixel of this y-value
  yfract = yabs - floor(yabs);

  // return the fractional y-value
  return yfract;
}

/**
 * Function: get_xvalue_from_tlength
 * The function computes the x-offset position from the reference point
 * for a given tracelength in a given beam.
 * The solution is numerically derived and therefore
 * does work for any reasonable tracefunction.
 *
 * Parameters:
 * @param actbeam - the beam to compute the x-offset
 * @param tlength - tracelength
 *
 * Returns:
 * @return xdiff - the x-offset
 */
double
get_xvalue_from_tlength(beam *actbeam, double tlength)
{
  int iter=0;
  int status;

  double xdiff=0.0;
  d_point x_interv;

  // define and initialize the solver
  const gsl_root_fsolver_type *T = gsl_root_fsolver_brent;
  gsl_root_fsolver            *s = gsl_root_fsolver_alloc (T);

  gsl_function F;
  tlength_pars *tpars;

  // derive the x-interval from the beam boundaries
  x_interv = get_xinterv_from_beam(actbeam);
  //  fprintf(stdout, "xpos: %f, ypos: %f\n", x_interv.x, x_interv.y);

  // allocate and fill the parameters
  tpars = (tlength_pars *) malloc(sizeof(tlength_pars));
  tpars->actbeam = actbeam;
  tpars->tlength = tlength;

  // fille the GSL-function
  F.function = &zero_tlength;
  F.params =   tpars;

  // set the boundaries for the solver
  gsl_root_fsolver_set (s, &F, x_interv.x, x_interv.y);

  // iterate to find the zeropoint
  do
    {
      // increment the iteration counter
      iter++;

      // iterate on the solver
      status = gsl_root_fsolver_iterate (s);

      // get a new guess from the solver
      xdiff = gsl_root_fsolver_root (s);

      // derive and set new boundaries
      x_interv.x = gsl_root_fsolver_x_lower (s);
      x_interv.y = gsl_root_fsolver_x_upper (s);

      // check the accuracy
      status = gsl_root_test_interval (x_interv.x, x_interv.y,
				       0, 0.0001);
      //--------------CAVEAT---------------------------------------
      // until March 28th 08 the code, wriongly was:
      //status = gsl_root_test_interval (x_interv.x, x_interv.x,
      //                                  0, 0.0001);
      // somehow this made no difference.....
      //-----------------------------------------------------------
    }
  // check for the break condition
  while (status == GSL_CONTINUE && iter < MAX_ITER);

  // free the memory
  free(tpars);

  // free the memory
  gsl_root_fsolver_free (s);

  // return the result
  return xdiff;
}

/**
 * Function: zero_tlength
 * Helper function for the 1D-root finding routine.
 * Evaluates a functional value at the position given
 * in the first parameter using the values given
 * as second parameter.
 *
 * Parameters:
 * @param x      - the independent value
 * @param params - ingredients to evaluate the function value
 *
 * Returns:
 * @return tzero - the function value
 */
double
zero_tlength (double x, void *params)
{
  double tzero;

  // extract the components from the parameter-structure
  beam  *actbeam = ((tlength_pars *) params)->actbeam;
  double tlength = ((tlength_pars *) params)->tlength;

  // compute the function value
  tzero = actbeam->spec_trace->path_len (x, actbeam->spec_trace->data) - tlength;

  // return the function value
  return tzero;
}

/**
 * Function: get_xinterv_from_beam
 * Determines the maximum and minimum extension of a given
 * beam along the x-axis. For security the interval is extended
 * at both ends by a fixed amount (quantified in the header-file).
 *
 * Parameters:
 * @param actbeam - the beam
 *
 * Returns:
 * @return x_interv - xmin, xmax covered by the beam
 */
d_point
get_xinterv_from_beam(const beam *actbeam)
{
  d_point x_interv;

  // compute the beam extension
  // towards the negative x-axis
  x_interv.x = (double)MIN(actbeam->corners[0].x,
			   MIN(actbeam->corners[1].x,
			       MIN(actbeam->corners[2].x,
				   actbeam->corners[3].x)))-actbeam->refpoint.x - INTERVEXT;

  // compute the beam extension
  // towards the positive x-axis
  x_interv.y = (double)MAX(actbeam->corners[0].x,
			   MAX(actbeam->corners[1].x,
			       MAX(actbeam->corners[2].x,
				   actbeam->corners[3].x)))-actbeam->refpoint.x + INTERVEXT;

  // return the interval
  return x_interv;
}

/**
 * Function: get_tlength_from_dispersion
 * The function computes and returnes the trace length for
 * a given diseprsion solution and wavelength. The tracelength
 * value is returned. Currently only linear and quadratic
 * dispersio solution are considered.
 *
 * Parameters:
 * @param cdisp  - the dispersion coefficient polynomial
 * @param lambda - the wavelength
 *
 * Returns:
 * @return tlength - the trace length
 */
double
get_tlength_from_dispersion(gsl_vector *cdisp, float lambda)
{
  double tlength=0.0;
  //double tlength_a=0.0;
  double det;

  // check whether the polynomial is linear
  if (cdisp->size == 2)
    {
      // compute the tracelength for linear dispersion
      tlength = (lambda-gsl_vector_get(cdisp, 0))/gsl_vector_get(cdisp, 1);
    }
  // check whether the polynomial is quadratic
  else if (cdisp->size == 3)
    {
      // compute the determinante
      det = gsl_vector_get(cdisp, 1) * gsl_vector_get(cdisp, 1)
		 - 4.0 * gsl_vector_get(cdisp, 0) * gsl_vector_get(cdisp, 2);

      // gove error if det < 0.0
      if (det < 0.0)
	aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		     "get_tlength_from_dispersion: Can not determine the tracelength since det < 0.0!\n");

      // compute the tracelength
      tlength   = (-gsl_vector_get(cdisp, 1)+sqrt(det))/(2.0*gsl_vector_get(cdisp, 2));
      //    tlength_a = (-gsl_vector_get(cdisp, 1)-sqrt(det))/(2.0*gsl_vector_get(cdisp, 2));
      //    fprintf(stdout, "plus: %f, minus: f\n", tlength, tlength_a);
    }
  else
    {
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "get_tlength_from_dispersion: The dispersion solution has %i coefficients and can not be inverted!\n", cdisp->size);
    }

  // return the tracelength
  return tlength;
}

/**
 * Function: condense_dispersion
 * The function strips off leading zeroes in the high order terms
 * of the dispersion solution derived from the configuration file.
 * These usually become a problem in the inversion to derive the
 * tracelength.
 *
 * Parameters:
 * @param disp__input - the configuration file dispersion
 *
 * Returns:
 * @return cdisp - the effective dispersion relation
 */
gsl_vector *
condense_dispersion(gsl_vector *disp_input)
{
  int index;
  int i;

  gsl_vector *cdisp;

  // get the size of the incomingg dispersion
  // vector
  index = disp_input->size;

  // determine the true size by subtracting
  // high order zeros
  while (index > -1 && !gsl_vector_get(disp_input, index-1))
    index--;

  // allocate space for the new vector
  cdisp = gsl_vector_alloc(index);

  // tranasfer the relevant values
  // from the old to the new vector
  for (i=0; i < index; i++)
    gsl_vector_set(cdisp, i, gsl_vector_get(disp_input, i));

  // return the new dispersion vector
  return cdisp;
}

/**
 * Function: fitting_ipc_corr
 *
 * Parameters:
 * @param actbeam        - the beam to examine
 * @param conf_file_path - the full pathname too the configuration file
 * @param ipcorr         - the interpolator with the correction factor
 * @param SPC            - the full spectrum to correct
 * @param obs            - the grism image
 * @param data_matrix    - background subtracted grism image
 * @param ipc_file_path  - file name for the phase data
 *
 * Returns:
 * @return icorr - marker whether the correction was applied or not (1/0)
 */
int
fitting_ipc_corr(beam act_beam, char conf_file_path[], interpolator *ipcorr,
		 full_spectr *SPC, observation *obs, gsl_matrix *data_matrix,
		 char ipc_file_path[], beam *beam_ptr)
{
  aperture_conf *conf;

  d_point        fit_qual;

  int            icorr     = 0;

  trace_func    *tracefun  = act_beam.spec_trace;
  double        *tracedata =  (double *)(tracefun->data);

  // load the configuration file
  conf = get_aperture_descriptor(conf_file_path);

  // determine the phase data and fix the phase
  // by fitting a function to it
  fit_qual = kappa_sigma_klipp_ipc(conf, obs, data_matrix, act_beam, ipc_file_path);

  // report on the shift and the quality
  fprintf(stdout, "shift-difference: %e, quality: %e\n", fit_qual.x, fit_qual.y);

  // check whether the quality is good enough
  if (fit_qual.y < MAXQUALITY)
    {
      // apply a correction to the
      // reference point to achive the
      // right corrections
      act_beam.refpoint.y += fit_qual.x;
      beam_ptr->refpoint.y += fit_qual.x;
      beam_ptr->ignore = -100;

      // report the new y-shift
      fprintf(stdout, "new y-shift: %e\n", act_beam.refpoint.y - (double)(int)(act_beam.refpoint.y + tracedata[1]));

      // apply the sensitivity correction to the spectrum
      intpix_corr_beam(act_beam, conf_file_path, ipcorr, SPC);

      // set the correction marker
      icorr=1;
    }

  // free allocated memory
  free_aperture_conf(conf);

  // return the marker
  return icorr;
}

/**
 * Function: intpix_corr_beam
 * The functio corrects a full spectrum for the intrapixel sensitivity
 * variations.
 * For every spectral element its fractional y-value is determined,
 * and then the correction factor for this y-value is computed
 * and applied to the appropriate elements of the spectral bin.
 *
 * Parameters:
 * @param actbeam        - the beam to examine
 * @param conf_file_path - the full pathname too the configuration file
 * @param ipcorr         - the interpolator with the correction factor
 * @param SPC            - the full spectrum to correct
 *
 * Returns:
 * @return -
 */
void
intpix_corr_beam(beam actbeam, char conf_file_path[], interpolator *ipcorr,
		 full_spectr *SPC)
{
  int            index=0;
  int            for_grism=1;

  double         cfactor;

  double         lambda;

  gsl_vector     *cdisp;

  d_point        pixel;

  dispstruct     *beam_disp;

  aperture_conf  *conf;

  // load the configuration file
  conf = get_aperture_descriptor(conf_file_path);

  // check whether it is grism (for_grism=1)
  // or prism (for_grism=0) data
  // give an error if there is a prism solution
  for_grism = check_for_grism (conf_file_path, actbeam.ID);
  if (!for_grism)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "intpix_corr_beam: Only grism dispersion solution can be corrected.\n");

  // get the wavelength dispersion relation at
  // position "refpoint". conf->refx and conf->refy
  // are used at this point to allow for a non (0,0) centered
  // 2D field dependence.
  pixel.x = actbeam.refpoint.x - conf->refx;
  pixel.y = actbeam.refpoint.y - conf->refy;

  // derive the dispersion at the object position
  beam_disp = get_dispstruct_at_pos(conf_file_path, for_grism,
				    actbeam.ID, pixel);

  // skipp high order zeroes in the dispersion solution
  cdisp = condense_dispersion(beam_disp->pol);

  for (index=0; index < SPC->nelems; index++)
    {
      lambda = SPC->fgr_spec->spec[index].lambda_mean;

      if (!gsl_isnan (lambda) && lambda)
	{
	  cfactor = get_intpix_corr(&actbeam, cdisp, ipcorr, lambda);

    	  fprintf(stdout, "lambda: %e, factor: %e\n", lambda, cfactor);
	  /*
	  SPC->fgr_spec->spec[index].count  = SPC->fgr_spec->spec[index].count / cfactor;
	  SPC->fgr_spec->spec[index].error  = SPC->fgr_spec->spec[index].error / cfactor;
	  SPC->fgr_spec->spec[index].flux   = SPC->fgr_spec->spec[index].flux / cfactor;
	  SPC->fgr_spec->spec[index].ferror = SPC->fgr_spec->spec[index].ferror / cfactor;

	  SPC->bck_spec->spec[index].count  = SPC->bck_spec->spec[index].count / cfactor;
	  SPC->bck_spec->spec[index].error  = SPC->bck_spec->spec[index].error / cfactor;
	  SPC->bck_spec->spec[index].flux   = SPC->bck_spec->spec[index].flux / cfactor;
	  SPC->bck_spec->spec[index].ferror = SPC->bck_spec->spec[index].ferror / cfactor;

	  SPC->obj_spec->spec[index].count  = SPC->obj_spec->spec[index].count / cfactor;
	  SPC->obj_spec->spec[index].error  = SPC->obj_spec->spec[index].error / cfactor;
	  SPC->obj_spec->spec[index].flux   = SPC->obj_spec->spec[index].flux / cfactor;
	  SPC->obj_spec->spec[index].ferror = SPC->obj_spec->spec[index].ferror / cfactor;

	  this is the wqrong version!!
	  SPC->fgr_spec->spec[index].count  = SPC->fgr_spec->spec[index].count * cfactor;
	  SPC->fgr_spec->spec[index].error  = SPC->fgr_spec->spec[index].error * cfactor;
	  SPC->fgr_spec->spec[index].flux   = SPC->fgr_spec->spec[index].flux * cfactor;
	  SPC->fgr_spec->spec[index].ferror = SPC->fgr_spec->spec[index].ferror * cfactor;

	  SPC->bck_spec->spec[index].count  = SPC->bck_spec->spec[index].count * cfactor;
	  SPC->bck_spec->spec[index].error  = SPC->bck_spec->spec[index].error * cfactor;
	  SPC->bck_spec->spec[index].flux   = SPC->bck_spec->spec[index].flux * cfactor;
	  SPC->bck_spec->spec[index].ferror = SPC->bck_spec->spec[index].ferror * cfactor;

	  SPC->obj_spec->spec[index].count  = SPC->obj_spec->spec[index].count * cfactor;
	  SPC->obj_spec->spec[index].error  = SPC->obj_spec->spec[index].error * cfactor;
	  SPC->obj_spec->spec[index].flux   = SPC->obj_spec->spec[index].flux * cfactor;
	  SPC->obj_spec->spec[index].ferror = SPC->obj_spec->spec[index].ferror * cfactor;
	  */
	}
    }

  // free the configuration structure
  free_aperture_conf(conf);

  // free the dispersion struct
  free_dispstruct(beam_disp);

  // free the memory for the dispersion
  gsl_vector_free(cdisp);
}

/**
 * Function: get_ipc_lambdas
 * The function computes the wavelengths for a list of x-offsets from the
 * reference point of a certain beam. The wavelength values are
 * returned as a gsl-vector.
 *
 * Parameters:
 * @param actbeam        - the beam
 * @param conf_file_path - the full path to the aXe configuration file
 * @param xvalues        - the list of x-offsets
 *
 * Returns:
 * @return lambdas - the list of wavelength values
 */
gsl_vector *
get_ipc_lambdas(const beam actbeam, char conf_file_path[], gsl_vector *xvalues)
{
  aperture_conf  *conf;

  dispstruct     *disp;
  calib_function *wl_calib;

  gsl_vector *lambdas;

  d_point pixel;
  int for_grism;

  int i;

  // load the configuration file
  conf = get_aperture_descriptor(conf_file_path);

  // check whether it is grism (for_grism=1)
  // or prism (for_grism=0) data
  for_grism = check_for_grism (conf_file_path, actbeam.ID);

  // determine the referencee point position
  pixel.x = actbeam.refpoint.x - conf->refx;
  pixel.y = actbeam.refpoint.y - conf->refy;

  // determine the dispersion at the reference point
  disp = get_dispstruct_at_pos(conf_file_path, for_grism,
			       actbeam.ID,pixel);

  // make a calibration structure from the dispersion
  wl_calib = create_calib_from_gsl_vector(for_grism, disp->pol);

  // convert the x-values to tracelength-values
  abscissa_to_pathlength (actbeam.spec_trace, xvalues);

  // allocate memory for the wavelengths
  lambdas = gsl_vector_alloc(xvalues->size);

  // go over all tracelength values
  for (i=0; i < (int)xvalues->size; i++)
    {
      // determine and store the wavelength for each tracelength
      gsl_vector_set(lambdas, i,
		     wl_calib->func(gsl_vector_get(xvalues, i), wl_calib->order,
				    wl_calib->coeffs));
    }

  // free the configuration structure
  free_aperture_conf(conf);

  // free the memory in the calibration structure
  free_calib(wl_calib);

  // free the dispersion structure
  free_dispstruct(disp);

  // return the wavelengths
  return lambdas;
}

/**
 * Function: get_ipc_cvalues
 * The function computes the intra-pixel correction factors for a certain
 * beam on a set of trace positions specified by their x-offset from
 * the reference point.
 * For every x-offset position the trace positon and its fractional y-pixel
 * (in absolute coordinates) is evaluated. Then the correction factor is
 * determined using the input interpolator.
 *
 * Parameters:
 * @param actbeam - the beam to correct
 * @param ipcorr  - the correction values depending on fractional y-pixel
 * @param xvalues - the list x-offsets from the reference point
 *
 * Returns:
 * @return cvalues - the list of correction values
 */
gsl_vector *
get_ipc_cvalues(const beam actbeam, interpolator *ipcorr, gsl_vector *xvalues)
{
  gsl_vector *cvalues;
  double yfract=0.0;
  int i;

  // allocate memory
  cvalues = gsl_vector_alloc(xvalues->size);

  // go over all x-values
  for (i=0; i < (int)xvalues->size; i++)
    {
      // determine the y-fraction for the x-value
      yfract = get_yfract_for_xvalue(&actbeam, gsl_vector_get(xvalues, i));

      // detyermine and store the correction factor in the array
      gsl_vector_set(cvalues, i, eval_interp(ipcorr, yfract));
    }

  // return the array
  return cvalues;
}

/**
 * Function: get_ipc_xvalues
 * The function computes a list of x-offsets from the reference
 * point of a beam. the regularly space offsets span the range
 * of x-values covered by the pixels of the beam.
 *
 * Parameters:
 * @param actbeam  - the beam to compute the x-offsets for
 *
 * Returns:
 * @return xvalues - the list of x-offsets
 */
gsl_vector *
get_ipc_xvalues(const beam actbeam)
{
  gsl_vector *xvalues;
  d_point xrange;

  int npoints;
  int i;

  // determine the x-range covered by the beam
  xrange =  get_xinterv_from_beam(&actbeam);

  // determine the number of points within the x-range
  npoints = (xrange.y - xrange.x) / XSTEPSIZE + 1;

  // allocate and fill a proper vector
  // with the x-values
  xvalues = gsl_vector_alloc(npoints);
  for (i=0; i < npoints; i++)
    {
      gsl_vector_set(xvalues, i, xrange.x + (double)i * XSTEPSIZE);
    }

  // return the vector with the x-values
  return xvalues;
}

/**
 * Function: get_ipclambda
 * The function creates an interpolator for the intra-pixel correction
 * as a function of wavelength based on this correction as function
 * of fractional y-pixel and the full calibration information
 * on a beam.
 * The interpolator is created after stepping along the beam trace
 * and combining the wavelength values with the correction values
 * at the trace positions.
 *
 * Parameters:
 * @param actbeam        - the beam to find the correction values for
 * @param conf_file_path - the full path to the aXe cofiguration file
 * @param ipcorr         - the correction values as function of y-fraction
 *
 * Returns:
 * @return ipclambda     - the correction values as function of wavelength
 */
interpolator *
get_ipclambda(beam actbeam, char conf_file_path[],
	      interpolator *ipcorr)
{
  gsl_vector *xvalues;
  gsl_vector *cvalues;
  gsl_vector *lambdas;

  double *cv;
  double *lv;

  interpolator *ipclambda;

  int i=0;
  int j=0;

  // get the x values around the reference point
  xvalues = get_ipc_xvalues(actbeam);

  // get the correction factors for these x-values
  cvalues = get_ipc_cvalues(actbeam, ipcorr, xvalues);

  // get the wavelength at the correction factors
  lambdas = get_ipc_lambdas(actbeam, conf_file_path, xvalues);

  // allocate memory for the dependent values
  cv = (double *) malloc(cvalues->size * sizeof(double));
  if (!cv) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		 "Memory allocation failed");
  }

  // allocate memory for the independent values
  lv = (double *) malloc(cvalues->size * sizeof(double));
  if (!lv) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		 "Memory allocation failed");
  }

  if (gsl_vector_get(lambdas, lambdas->size-1) > gsl_vector_get(lambdas, 0))
    {
      // transfer the values from the
      // gsl vectors to the c-vectors
      for (i = 0; i < (int)cvalues->size; i++)
	{
	  cv[i] = gsl_vector_get(cvalues, i);
	  lv[i] = gsl_vector_get(lambdas, i);
	}
    }
  else
    {
      // transfer the values from the
      // gsl vectors to the c-vectors
      // invert the order from the
      // gsl-vectors
      j = cvalues->size - 1;
      for (i = 0; i < (int)cvalues->size; i++)
	{
	  cv[j] = gsl_vector_get(cvalues, i);
	  lv[j] = gsl_vector_get(lambdas, i);
	  j--;
	}
    }

  // create the interpolator
  ipclambda = create_interp(cvalues->size, FILTER_INTERP_TYPE, lv, cv);

  // free the memory allocated to
  // the vectors
  gsl_vector_free(xvalues);
  gsl_vector_free(cvalues);
  gsl_vector_free(lambdas);

  // return the interpolator
  return ipclambda;
}

/**
 * Function: apply_corr_pet
 * The function applies an intra-pixel sensitivity correction to
 * a list of PET pixels. The correction values are given depending
 * on the wavelength of the pixel, and are applied to the PET pixels
 * in place.
 *
 * Parameters:
 * @param ipclambda - the ipc correction factor as function of wavelength
 * @param PET       - the list of PET pixels to correct
 *
 * Returns:
 * -
 */
void
apply_corr_pet(interpolator *ipclambda, ap_pixel *PET)
{
  int j = 0;

  double cvalue=0.0;

  // go along the PET pixels until you meet
  // the last
  while (PET[j].p_x != -1)
    {
      // get the correction factor
      cvalue = eval_interp(ipclambda, PET[j].lambda);

      // apply the corection factor
      // to the counts and the error
      PET[j].count = PET[j].count / cvalue;
      PET[j].error = PET[j].error / cvalue;

      //      fprintf(stdout, "lambda: %f, correction: %f\n", PET[j].lambda, cvalue);
      j++;
    }

}

/**
 * Function: intpix_corr_pet
 * The function applies the intra-pixel sensitivity correction to
 * a list of PET pixels. The values in the pixels are corrected in
 * place using the correction function, gemoetrical parameters and
 * dispersion relation in the various parameters.
 *
 * Parameters:
 * @param actbeam        - the beam to correct
 * @param conf_file_path - full path to the axe configuration file
 * @param ipcorr         - the ipc correction function
 * @param PET            - the list of PET pixels to correct
 *
 * Returns:
 * -
 */
void
intpix_corr_pet(beam actbeam, char conf_file_path[],
		interpolator *ipcorr, ap_pixel *PET)
{
  interpolator *ipclambda;
  aperture_conf  *conf;

  int for_grism=0;

  // load the configuration file
  conf = get_aperture_descriptor(conf_file_path);

  // check whether it is grism (for_grism=1)
  // or prism (for_grism=0) data
  // give an error if there is a prism solution
  for_grism = check_for_grism (conf_file_path, actbeam.ID);
  if (!for_grism)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "intpix_corr_beam: Only grism dispersion solution can be corrected.\n");

  // determine the correction factor versus wavelength
  ipclambda = get_ipclambda(actbeam, conf_file_path, ipcorr);

  print_interp(ipclambda);

  // apply the correction to the PET
  apply_corr_pet(ipclambda, PET);

  // free the configuration structure
  free_aperture_conf(conf);

  // free the memory in the interpolator
  free_interp(ipclambda);
}

/**
 * Function: is_pointlike
 * The function evaluates all criteria to decide whether an object
 * is pointlike or not. In the current implementation an object
 * is considered pointlike if a special OAF is given and the
 * the beam is NOT excluded OR if the spatial extension of the
 * object is smaller than the maximal extension.
 *
 * Parameters:
 * @param actbeam   - the beam to examine
 * @param spec_OAF  - indicates a special OAF file
 * @param max_ext   - the maximal allowed extension
 *
 * Returns:
 * @return point_like - pointer to the opened fits file
 */
int
is_pointlike(beam actbeam, int spec_OAF, double max_ext)
{
  int point_like=0;

  int OAF_crit=0;
  int EXT_crit=0;

  // check whether a special OAF is given
  // and the beam should NOT be ignored
  if (spec_OAF && !actbeam.ignore)
    OAF_crit=1;

  // check whether the maximal extension
  // is given and the object extension
  // is smaller;
  // also check whether non-default values
  // for the object width are set
  if (actbeam.awidth > -1.0 && actbeam.bwidth > -1.0 && max_ext && actbeam.awidth < max_ext && actbeam.bwidth < max_ext)
    EXT_crit=1;

  // set to pointlike if at least
  // one of the criteria is set
  if (OAF_crit || EXT_crit)
    point_like=1;

  // return the result
  return point_like;
}

/**
 * Function: get_SPC_opened
 * The function opens an existing SPC file and returns the pointer
 * to it.
 * As of now, the mode of opening it is automatically READWRITE.
 * Later on a differentiation using the free parameter "mode"
 * might be added  to generalize the function.
 *
 * Parameters:
 * @param SPCname - name of the SPC file
 * @param mode    - mode to open it (not yet used)
 *
 * Returns:
 * @return SPC_ptr - pointer to the opened fits file
 *
 */
fitsfile *
get_SPC_opened(char SPCname[], int mode)
{
    fitsfile *SPC_ptr;
    int f_status=0;

  // Open the OPET file for reading/writing
  fits_open_file (&SPC_ptr, SPCname, READWRITE, &f_status);
  if (f_status)
    {
      ffrprt (stdout, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "aXe_INTPIXCORR: Could not open file: %s\n",
		   SPCname);
    }

  // return the pointer to the fits file
  return SPC_ptr;
}

/**
 * Function: get_ALL_from_next_in_SPC
 * The function creates and fills a full spectrum structure with
 * the content of a SPC table extension. This is only done when,
 * according to the beam ID, this beam should be corrected.
 * For extensions with higher order beams which are not corrected
 * an emply structure is returned.
 *
 * Parameters:
 * @param SPCn_ptr - pointer to the opened SPC file
 * @param aperID   - pointer to aperture identification number
 * @param beamID   - pointer to beam identification number
 *
 * Returns:
 * @return SPC - the full spectrum structure
 */
full_spectr *
get_ALL_from_next_in_SPC(fitsfile *SPC_ptr, int *aperID, int *beamID)
{
  int f_status=0, hdutype;

  long tmp;
  //long nrows=0;
  char comment[FLEN_COMMENT];

  full_spectr *SPC;

  fits_movrel_hdu (SPC_ptr, 1, &hdutype, &f_status);

  if (f_status)
    {
      *aperID = -1;
      *beamID = -1;
      SPC = NULL;
      return SPC;
    }

  // read the beam ID number
  fits_read_key_lng (SPC_ptr, "BEAMID", &tmp, comment, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "get_ALL_from_next_in_SPC: Error getting index keyword OBJECTID");
    }
  *beamID = (int)tmp;

  // check whether this beam should be correct
  if (*beamID > CORRMAX)
    {
      // set it to NULL and return
      SPC = NULL;
      //      fprintf (stdout, "aXe_PETFF: Skipping beam: %c.\n", BEAM(*beamID));
      return SPC;
    }

  // the beam shall be corrected and
  // first must be read in
  SPC = (full_spectr *) malloc (sizeof (full_spectr ));

  // transfer the beam ID
  SPC->beamID = (int)tmp;

  // read the aperture number
  fits_read_key_lng (SPC_ptr, "OBJECTID", &tmp, comment, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "get_ALL_from_next_in_SPC: Error getting index keyword OBJECTID");
    }
  // transfer the aperture ID
  *aperID = (int)tmp;
  SPC->aperID = (int)tmp;

  // Get the number of rows
  fits_get_num_rows (SPC_ptr, &tmp, &f_status);
  if (f_status) {
    ffrprt (stderr, f_status);
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "get_ALL_from_next_in_SPC: "
		 "Could not determine the number of rows in"
		 " correction function table!");
  }
  SPC->nelems = (int)tmp;

  // load the background subtracted object spectrum
  SPC->obj_spec = get_spectrum_from_SPC(SPC_ptr, "COUNT", "ERROR", SPC->nelems);

  // load the total object spectrum
  SPC->fgr_spec = get_spectrum_from_SPC(SPC_ptr, "TCOUNT", "TERROR", SPC->nelems);

  // load the background spectrum
  SPC->bck_spec = get_spectrum_from_SPC(SPC_ptr, "BCOUNT", "BERROR", SPC->nelems);

  // return the filled structure
  return SPC;
}

/**
 * Function: free_full_spectr
 * The function frees the memory allocated in a full
 * spectrum structure.
 *
 * Parameters:
 * @param SPC - the full spectrum structure
 *
 * Returns:
 * @return -
 */
void
free_full_spectr(full_spectr *SPC)
{
  // free the three spectra in the
  // full spectrum structure
  free_spectrum(SPC->obj_spec);
  free_spectrum(SPC->fgr_spec);
  free_spectrum(SPC->bck_spec);

  // free the full spectrum
  free(SPC);

  // set the structure to NULL
  SPC = NULL;
}

/**
 * Function: get_spectrum_from_SPC
 * The function fills a spectrum structure with data in an SPC extension.
 * There is more data in an SPC extension than can be stored in a
 * spectrum structure. Two parameters in this function specify partly which
 * data should be loaded.
 *
 * Parameters:
 * @param SPC_ptr   - pointer to the SPC extension
 * @param count_col - column to load for 'count'
 * @param error_col - column to load for 'error'
 * @param nelems    - the number ol elements in the spectrum
 *
 * Returns:
 * @return act_spec - the filled spectrum structure
 */
spectrum *
get_spectrum_from_SPC(fitsfile *SPC_ptr, char count_col[],
		      char error_col[], int nelems)
{
  int index=0;

  spectrum *act_spec;

  double *count;
  double *lambda;
  double *error;
  double *flux;
  double *ferror;
  double *weight;
  double *contam;
  long   *dq;

  // allocate the spectrum
  act_spec = allocate_spectrum (nelems);

  // transfer the column entries into a vector
  lambda = get_dcolumn_from_SPC_opened(SPC_ptr, "LAMBDA", nelems);
  count  = get_dcolumn_from_SPC_opened(SPC_ptr, count_col, nelems);
  error  = get_dcolumn_from_SPC_opened(SPC_ptr, error_col, nelems);
  flux   = get_dcolumn_from_SPC_opened(SPC_ptr, "FLUX", nelems);
  ferror = get_dcolumn_from_SPC_opened(SPC_ptr, "FERROR", nelems);
  weight = get_dcolumn_from_SPC_opened(SPC_ptr, "WEIGHT", nelems);
  contam = get_dcolumn_from_SPC_opened(SPC_ptr, "CONTAM", nelems);
  dq     = get_lcolumn_from_SPC_opened(SPC_ptr, "DQ", nelems);

  // transfer the data from the vector into the
  // spectrum structure
  for (index = 0; index < nelems; index++)
    {
      act_spec->spec[index].lambda_mean = lambda[index];
      act_spec->spec[index].count       = count[index];
      act_spec->spec[index].error       = error[index];
      act_spec->spec[index].flux        = flux[index];
      act_spec->spec[index].ferror      = ferror[index];
      act_spec->spec[index].weight      = weight[index];
      act_spec->spec[index].contam      = contam[index];
      act_spec->spec[index].dq          = dq[index];
    }

  // free the arrays
  free(lambda);
  free(count);
  free(error);
  free(flux);
  free(ferror);
  free(weight);
  free(contam);
  free(dq);

  // return the spectrum
  return act_spec;
}

/**
 * Function: get_dcolumn_from_SPC_opened
 * The function reads the values from a double column specified
 * by its name into a array of doubles. The array is returned.
 *
 * Parameters:
 * @param SPC_ptr   - pointer to the SPC extension
 * @param count_col - name of the column to load
 * @param nelems    - number of elements in the column
 *
 * Returns:
 * @return values - pointer to a filled array
 */
double *
get_dcolumn_from_SPC_opened(fitsfile *SPC_ptr, char colname[], int nelems)
{
  int colnum=0;
  int anynul;
  int f_status=0;
  double *values;

  // allocate the return array;
  // give an error if allocation fails
  values = (double *) malloc(nelems * sizeof(double));
  if (!values) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		 "Memory allocation failed");
  }

  // get the desired column number;
  // give an error if the column name
  // can not be read
  fits_get_colnum (SPC_ptr, CASEINSEN, colname, &colnum, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "get_dcolumn_from_SPC_opened: "
		   "Could not determine column %s in "
		   " table!\n", colname);
    }

  // read all data in the column
  fits_read_col (SPC_ptr, TDOUBLE, colnum, 1, 1, nelems, NULL, values,
                    &anynul, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "get_dcolumn_from_SPC_opened: "
		   "Could not read column %s"
		   " from BINARY table!",colname);
    }

  // return the filled vector
  return values;
}

/**
 * Function: get_lcolumn_from_SPC_opened
 * The function reads the values from a column with long specified
 * by its name into a array of type long. The array is returned.
 *
 * Parameters:
 * @param SPC_ptr   - pointer to the SPC extension
 * @param count_col - name of the column to load
 * @param nelems    - number of elements in the column
 *
 * Returns:
 * @return values - pointer to a filled array
 */
long *
get_lcolumn_from_SPC_opened(fitsfile *SPC_ptr, char colname[], int nelems)
{
  int colnum=0;
  int anynul;
  int f_status=0;
  long *values;

  // allocate the return array;
  // give an error if allocation fails
  values = (long *) malloc(nelems * sizeof(long));
  if (!values) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		 "Memory allocation failed");
  }

  // get the desired column number;
  // give an error if the column name
  // can not be read
  fits_get_colnum (SPC_ptr, CASEINSEN, colname, &colnum, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "get_dcolumn_from_SPC_opened: "
		   "Could not determine column %s in "
		   " table!\n", colname);
    }

  // read all data in the column
  fits_read_col (SPC_ptr, TLONG, colnum, 1, 1, nelems, NULL, values,
                    &anynul, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "get_dcolumn_from_SPC_opened: "
		   "Could not read column %s"
		   " from BINARY table!",colname);
    }

  // return the filled array
  return values;
}

/**
 * Function: create_nlincor
 * The function creates an interpolator for the nonlinearity
 * correction applied to NICMOS data.
 *
 * Parameters:
 *
 * Returns:
 * @return nlincorr - the interpolator created
 */
interpolator *
create_nlincor()
{
  interpolator *nlincorr;
  /*
  double x[14] = {8250.0, 8750.0, 9250.0, 9750.0, 11000.0, 12000.0, 13000.0, 14000.0, 15000.0, 16000.0, 17000.0, 18000.0, 19000.0, 20000.0};
  double y[14] = {  .069,   .057,   .052,   .050,    .049,    .048,      .041,    .023,    .013,    .008,    .004,    .0,      .0,      .0};
  */
  double *xx;
  double *yy;

  xx = (double *) malloc(14 * sizeof(double));
  yy = (double *) malloc(14 * sizeof(double));

  xx[0]  = 8250.0;
  xx[1]  = 8750.0;
  xx[2]  = 9250.0;
  xx[3]  = 9750.0;
  xx[4]  = 11000.0;
  xx[5]  = 12000.0;
  xx[6]  = 13000.0;
  xx[7]  = 14000.0;
  xx[8]  = 15000.0;
  xx[9]  = 16000.0;
  xx[10] = 17000.0;
  xx[11] = 18000.0;
  xx[12] = 19000.0;
  xx[13] = 20000.0;

  yy[0]  = .069;
  yy[1]  = .057;
  yy[2]  = .052;
  yy[3]  = .050;
  yy[4]  = .049;
  yy[5]  = .048;
  yy[6]  = .041;
  yy[7]  = .023;
  yy[8]  = .013;
  yy[9]  = .008;
  yy[10] = .004;
  yy[11] = .0;
  yy[12] = .0;
  yy[13] = .0;

  // create the interpolator
  nlincorr = create_interp(14, NLINCORR_INTERP_TYPE, xx, yy);

  // return the interpolator
  return nlincorr;
}

/**
 * Function: nlin_corr_beam
 *
 * Parameters:
 * @param nlincorr       - the interpolator with the correction factor
 * @param SPC            - the full spectrum to correct
 *
 * Returns:
 * @return -
 */
void
nlin_corr_beam(interpolator *nlincorr, double adcgain, full_spectr *SPC)
{
  int            index=0;
  //  int            for_grism=1;

  double         cfactor;
  double         cps;

  double         lambda;

  for (index=0; index < SPC->nelems; index++)
    {
      // get the independent spectral values,
      // the wavelenth and the cps value
      lambda = SPC->obj_spec->spec[index].lambda_mean;
      cps    = SPC->obj_spec->spec[index].count;

      // check whether the spectral element
      // is not corrupt
      if (!gsl_isnan (lambda) && lambda)
	{
	  if (cps > 0.0)
	    // compute the correction factor
	    cfactor = get_nlin_corr(nlincorr, lambda, cps/adcgain);
	  else
	    // make a dummy factor
	    cfactor = 1.0;
	  // correct what should be corrected
	  SPC->fgr_spec->spec[index].count  = SPC->fgr_spec->spec[index].count / cfactor;
	  SPC->fgr_spec->spec[index].error  = SPC->fgr_spec->spec[index].error / cfactor;
	  SPC->fgr_spec->spec[index].flux   = SPC->fgr_spec->spec[index].flux / cfactor;
	  SPC->fgr_spec->spec[index].ferror = SPC->fgr_spec->spec[index].ferror / cfactor;

	  SPC->bck_spec->spec[index].count  = SPC->bck_spec->spec[index].count / cfactor;
	  SPC->bck_spec->spec[index].error  = SPC->bck_spec->spec[index].error / cfactor;
	  SPC->bck_spec->spec[index].flux   = SPC->bck_spec->spec[index].flux / cfactor;
	  SPC->bck_spec->spec[index].ferror = SPC->bck_spec->spec[index].ferror / cfactor;

	  SPC->obj_spec->spec[index].count  = SPC->obj_spec->spec[index].count / cfactor;
	  SPC->obj_spec->spec[index].error  = SPC->obj_spec->spec[index].error / cfactor;
	  SPC->obj_spec->spec[index].flux   = SPC->obj_spec->spec[index].flux / cfactor;
	  SPC->obj_spec->spec[index].ferror = SPC->obj_spec->spec[index].ferror / cfactor;
	}
    }

}

/*
 * Function: get_nlin_corr
 *
 * Parameters:
 * @param nlincorr - the interpolator with the parameter
 * @param lamda    - the full spectrum to correct
 * @param cps      - the count rate
 *
 * Returns:
 * @return cfactor - the correction factor
 */
double
get_nlin_corr(interpolator *nlincorr, const double lambda, const double cps)
{
  double cfactor;
  double bbb;

  // evaluate the parameter
  bbb = eval_interp(nlincorr, lambda);

  // compute the correction factor
  cfactor = 1.0 - 2.0 * bbb + bbb * log10(cps);

  // return the correction factor
  return cfactor;
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/**
 *
 * RV-GOMEA
 *
 * If you use this software for any purpose, please cite the most recent publication:
 * A. Bouter, C. Witteveen, T. Alderliesten, P.A.N. Bosman. 2017.
 * Exploiting Linkage Information in Real-Valued Optimization with the Real-Valued
 * Gene-pool Optimal Mixing Evolutionary Algorithm. In Proceedings of the Genetic
 * and Evolutionary Computation Conference (GECCO 2017).
 * DOI: 10.1145/3071178.3071272
 *
 * Copyright (c) 1998-2017 Peter A.N. Bosman
 *
 * The software in this file is the proprietary information of
 * Peter A.N. Bosman.
 *
 * IN NO EVENT WILL THE AUTHOR OF THIS SOFTWARE BE LIABLE TO YOU FOR ANY
 * DAMAGES, INCLUDING BUT NOT LIMITED TO LOST PROFITS, LOST SAVINGS, OR OTHER
 * INCIDENTIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR THE INABILITY
 * TO USE SUCH PROGRAM, EVEN IF THE AUTHOR HAS BEEN ADVISED OF THE POSSIBILITY
 * OF SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. THE AUTHOR MAKES NO
 * REPRESENTATIONS OR WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE
 * AUTHOR SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED BY ANYONE AS A RESULT OF
 * USING, MODIFYING OR DISTRIBUTING THIS SOFTWARE OR ITS DERIVATIVES.
 *
 * The software in this file is the result of (ongoing) scientific research.
 * The following people have been actively involved in this research over
 * the years:
 * - Peter A.N. Bosman
 * - Dirk Thierens
 * - Jörn Grahl
 * - Anton Bouter
 *
 */
#pragma once

#include <stdio.h>
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <eigen3/Eigen/Dense>
#include <cblas.h>
#include <lapack.h>

using Eigen::MatrixXd;

#define PI 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798
#define FALSE 0
#define TRUE 1

namespace GOMEA
{
void *
Malloc(long size);
void
shrunkCovariance(MatrixXd & emp_cov, const double alpha);
void
shrunkCovarianceOAS(MatrixXd & emp_cov, const int pop_size);
float
choleskyDecomposition(MatrixXd & result, MatrixXd & matrix, int n);
void
lowerTriangularMatrixInverse(MatrixXd & A, const int n);
int *
mergeSort(double * array, int array_size);
void
mergeSortWithinBounds(double * array, int * sorted, int * tosort, int p, int q);
void
mergeSortWithinBoundsInt(int * array, int * sorted, int * tosort, int p, int q);

void
mergeSortMerge(double * array, int * sorted, int * tosort, int p, int r, int q);
int *
mergeSortInt(int * array, int array_size);
void
mergeSortMergeInt(int * array, int * sorted, int * tosort, int p, int r, int q);

int *
getRanks(double * array, int array_size);
int *
getRanksFromSorted(int * sorted, int array_size);

double
randomRealUniform01(void);
int
randomInt(int maximum);
double
random1DNormalUnit(void);
double
random1DNormalParameterized(double mean, double variance);
void
initializeRandomNumberGenerator(void);
int *
randomPermutation(int n);
int **
allPermutations(int length, int * numberOfPermutations);
int **
allPermutationsSubroutine(int from, int length, int * numberOfPermutations);

double
max(double x, double y);
double
min(double x, double y);
double
distanceEuclidean(double * solution_a, double * solution_b, int n);
double
distanceEuclidean2D(double x1, double y1, double x2, double y2);

double *
matrixVectorPartialMultiplication(double ** matrix,
                                  double *  vector,
                                  int       n0,
                                  int       number_of_elements,
                                  int *     element_indices);

extern int64_t random_seed, /* The seed used for the random-number generator. */
  random_seed_changing;     /* Internally used variable for randomly setting a random seed. */

extern double haveNextNextGaussian, /* Internally used variable for sampling the normal distribution. */
  nextNextGaussian;
}

max_stars_repo_name: joasiee/elastix

max_issues_repo_name: joasiee/elastix

max_forks_repo_name: joasiee/elastix


text:
/*
Written by John MacCallum, The Center for New Music and Audio Technologies,
University of California, Berkeley.  Copyright (c) 2006-09, The Regents of
the University of California (Regents).
Permission to use, copy, modify, distribute, and distribute modified versions
of this software and its documentation without fee and without a signed
licensing agreement, is hereby granted, provided that the above copyright
notice, this paragraph and the following two paragraphs appear in all copies,
modifications, and distributions.

IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING
OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF REGENTS HAS
BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE. THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF ANY, PROVIDED
HEREUNDER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION TO PROVIDE
MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
NAME: randdist
DESCRIPTION: Random number generator with over 30 statistical distributions.
AUTHORS: John MacCallum
COPYRIGHT_YEARS: 2006-09
DRUPAL_NODE: /patch/4023
SVN_REVISION: $LastChangedRevision: 587 $
VERSION 1.1: Changed the way the random seed it made
VERSION 1.2: Universal binary
VERSION 1.2.1: Changed the license to be GPL compatible
VERSION 1.3: The object now takes arguments to specify the distribution
VERSION 1.3.1: Use of the buffering system is now optional and off by default
VERSION 1.3.2: Fixed a bug that would cause a crash if randdist was instantiated with a number instead of a symbol as its first arg.
VERSION 1.3.3: Added Gaussdist faker to helpfile, bump version to re-release. -mzed
VERSION 1.4: Default state, changed dump to distlist, user-defined pmfs are now done with the nonparametric message, bugfixes.
VERSION 2.0: Rewritten with common code between randdist and randdist~ put in libranddist.h/c
VERSION 2.1: Re-re-factored.  Nonparametric now uses gsl_ran_discrete().
VERSION 2.1.1: Fixed a bug that occurred when distlist was used with nonparametric
VERSION 2.1.2: Fixed a bug with the multinomial distribution
VERSION 2.1.3: Added multivariate hypergeometric distribution
VERSION 2.1.4: fixed a bug with the way that multinomial arguments were being handled
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
*/

#include "ext.h"
#include "ext_obex.h"
#include "ext_obex_util.h"
#include "libranddist.h"
#include "version.h"
#include "version.c"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#define RDIST_DEFAULT_BUF_SIZE 1024

typedef struct _rdist{
        t_object r_ob;
        void *r_out0;
	void *r_out1;
        gsl_rng *r_rng; // random number generator
	gsl_ran_discrete_t *r_g; // lookup table for nonparametric pdf sampling
        t_atom r_vars[R_MAX_N_VARS];
        t_symbol *r_dist;
        t_atom *r_arIn;
        short r_pmfLength;
        int r_numVars;
	int r_stride;
	void (*r_function)(gsl_rng*, int, void*, int, float*);
	t_atom *r_output_buffer;
} t_rdist;

void *rdist_class;

void rdist_anything(t_rdist *x, t_symbol *msg, short argc, t_atom *argv);
void rdist_bang(t_rdist *x);
void rdist_distlist(t_rdist *x, long n);
void rdist_assist(t_rdist *x, void *b, long m, long a, char *s);
void *rdist_new(t_symbol *msg, short argc, t_atom *argv);
void rdist_anything(t_rdist *x, t_symbol *msg, short argc, t_atom *argv);
void rdist_list(t_rdist *x, t_symbol *msg, short argc, t_atom *argv);
void rdist_nonparametric(t_rdist *x, t_symbol *msg, short argc, t_atom *argv);
void rdist_int(t_rdist *x, long n);
void rdist_float(t_rdist *x, double n);
void rdist_seed(t_rdist *x, long s);
void rdist_free(t_rdist *x);
void rdist_tellmeeverything(t_rdist *x);
void rdist_errorHandler(const char * reason, const char * file, int line, int gsl_errno);

int main(void){
	setup((t_messlist **)&rdist_class, (method)rdist_new, (method)rdist_free, (short)sizeof(t_rdist), 0L, A_GIMME, 0);

	addmess((method)version, "version", 0);
	addmess((method)rdist_assist, "assist", A_CANT, 0);
	addbang((method)rdist_bang);
	addmess((method)rdist_distlist, "distlist", A_DEFLONG, 0);
        addint((method)rdist_int);
        addfloat((method)rdist_float);
        addmess((method)rdist_list, "list", A_GIMME, 0);
	addmess((method)rdist_anything, "anything", A_GIMME, 0);
        addmess((method)rdist_nonparametric, "nonparametric", A_GIMME, 0);
        addmess((method)rdist_seed, "seed", A_LONG, 0);
        addmess((method)rdist_tellmeeverything, "tellmeeverything", 0);

	version(0);

	return 0;
}

void *rdist_new(t_symbol *msg, short argc, t_atom *argv){
	t_rdist *x;
	int i;
	t_atom ar[2];

	x = (t_rdist *)newobject(rdist_class); // create a new instance of this object

	x->r_out0 = outlet_new(x, 0);

	x->r_numVars = 0;

	// set up the random number generator
	gsl_rng_env_setup();

	// waterman14 was the fastest according to my tests
	x->r_rng = gsl_rng_alloc((const gsl_rng_type *)gsl_rng_waterman14);

	// seed it by reading from /dev/random on mac os x and
	// something similar on windows
	gsl_rng_set(x->r_rng, makeseed());

	// this is really fucking important.  if there's an error and the gsl's
	// default handler gets called, it aborts the program!
	gsl_set_error_handler(rdist_errorHandler);

	// setup a workspace
	x->r_output_buffer = (t_atom *)malloc(RDIST_DEFAULT_BUF_SIZE * sizeof(t_atom));

	// init the lib.  just gensyms all the distribution names
	librdist_init();

	// handle the args.  this should be done with attributes.
	if(argc){
		if(argv[0].a_type == A_SYM){
			rdist_anything(x, argv[0].a_w.w_sym, argc - 1, argv + 1);
		}
	} else {
		SETFLOAT(&(ar[0]), 0.);
		SETFLOAT(&(ar[1]), 1.);
		rdist_anything(x, gensym("uniform"), 2, ar);
	}

	return x;
}

void rdist_bang(t_rdist *x){
	int stride = x->r_stride;
	float out[stride];
	int i;
	t_symbol *msg;

	//msg = (stride > 1) ? gensym("list") : gensym("float");

	//post("buffer: %d, pos: %d, val: %f", x->r_whichBuffer, x->r_bufferPos, out[0].a_w.w_float);
	if(x->r_dist == ps_nonparametric){
		x->r_function(x->r_rng, 0, x->r_g, stride, out);
	}else{
		x->r_function(x->r_rng, x->r_numVars, x->r_vars, stride, out);
	}

	if(stride > 1){
		for(i = 0; i < stride; i++){
			SETFLOAT(x->r_output_buffer + i, out[i]);
		}
		outlet_list(x->r_out0, NULL, stride, x->r_output_buffer);
	}else{
		outlet_float(x->r_out0, out[0]);
	}
	//outlet_anything(x->r_out0, msg, stride, x->r_output_buffer);
}

void rdist_anything(t_rdist *x, t_symbol *msg, short argc, t_atom *argv){
	if(argc > R_MAX_N_VARS){
		error("randdist: too many variables");
		return;
	}

	// need to check that this is a valid dist.
	x->r_dist = msg;
	int i;
	if(x->r_dist == gensym("multinomial")){
		double sum = 0;
		for(i = 1; i < argc; i++){
			sum += atom_getfloat(argv + i);
		}
		x->r_vars[0] = argv[0];
		for(i = 1; i < argc; i++){
			atom_setfloat(x->r_vars + i, atom_getfloat(argv + i) / sum);
		}
	}else{
		memcpy(x->r_vars, argv, argc * sizeof(t_atom));
	}

	x->r_numVars = (int)argc;

	if(x->r_dist == ps_bivariate_gaussian){
		x->r_stride = 2;
	}else if(x->r_dist == ps_dirichlet){
		x->r_stride = x->r_numVars;
	}else if(x->r_dist == ps_multinomial || x->r_dist == ps_multivariate_hypergeometric){
		x->r_stride = x->r_numVars - 1;
	}else{
		x->r_stride = 1;
	}

	x->r_function = librdist_get_function(x->r_dist);
}

void rdist_nonparametric(t_rdist *x, t_symbol *msg, short argc, t_atom *argv){
	double f[argc];
	int i;
	x->r_dist = msg;
	for(i = 0; i < argc; i++){
		f[i] = librdist_atom_getfloat(argv + i);
		//post("%d, %f", i, f[i]);
	}
	if(x->r_g){
		gsl_ran_discrete_free(x->r_g);
	}
	x->r_g = gsl_ran_discrete_preproc(argc, f);
	x->r_function = librdist_nonparametric;
}

void rdist_list(t_rdist *x, t_symbol *msg, short argc, t_atom *argv){
	rdist_anything(x, x->r_dist, argc, argv);
}

void rdist_int(t_rdist *x, long n){
	t_atom argv[1];
	SETFLOAT(argv, (double)n);
	rdist_anything(x, x->r_dist, 1, argv);
}

void rdist_float(t_rdist *x, double n){
	t_atom argv[1];
	SETFLOAT(argv, n);
	rdist_anything(x, x->r_dist, 1, argv);
}

void rdist_distlist(t_rdist *x, long n){
	int i, j;
	int stride = x->r_stride;
	float out[n * stride];
	if(stride * n > RDIST_DEFAULT_BUF_SIZE){
		error("randdist: output list is too long (%d > %d)",
		      n * stride,
		      RDIST_DEFAULT_BUF_SIZE);
		return;
	}

	if(n < 1){
		error("randdist: distlist argument must be >= 1.");
		return;
	}

	for(i = 0; i < n; i++){
		if(x->r_dist == gensym("nonparametric")){
			x->r_function(x->r_rng, x->r_numVars, x->r_g, stride, out);
		}else{
			x->r_function(x->r_rng, x->r_numVars, x->r_vars, stride, out);
		}
		for(j = 0; j < stride; j++){
			SETFLOAT(x->r_output_buffer + ((i * stride) + j), out[j]);
		}
	}

	outlet_anything(x->r_out0, gensym("list"), n * stride, x->r_output_buffer);
}

void rdist_assist(t_rdist *x, void *b, long m, long a, char *s){
	if (m == ASSIST_INLET)
		sprintf(s,"Probability distribution and arguments");
	else {
		switch (a) {
		case 0:
			sprintf(s,"Random variate");
			break;
		}
	}
}

void rdist_tellmeeverything(t_rdist *x){
	int i;
	float vars[x->r_numVars];
	post("Distribution: %s", x->r_dist->s_name);
	post("Args:");
	for(i = 0; i < x->r_numVars; i++){
		post("\t%f", x->r_vars[i]);
	}
}

void rdist_free(t_rdist *x){
	if(x->r_rng){
		gsl_rng_free(x->r_rng);
	}
	if(x->r_g){
		gsl_ran_discrete_free(x->r_g);
	}
	if(x->r_arIn){
		free(x->r_arIn);
	}
	if(x->r_output_buffer){
		free(x->r_output_buffer);
	}
}

void rdist_errorHandler(const char * reason, const char * file, int line, int gsl_errno){
	error("randdist: a(n) %s has occured in file %s at line %d (error %d)", reason, file, line, gsl_errno);
}


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.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <mpi.h>
#include <petsc.h>
#include <petscvec.h>
#include <petscmat.h>
#include <petscksp.h>
#include <petscpc.h>
#include <petscsnes.h>

#include <petscversion.h>
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >=3) )
  #if (PETSC_VERSION_MINOR >=6)
     #include "petsc/private/kspimpl.h"
  #else
     #include "petsc-private/kspimpl.h"   /*I "petscksp.h" I*/
  #endif
#else
  #include "private/kspimpl.h"   /*I "petscksp.h" I*/
#endif

#include <StGermain/StGermain.h>
#include <StgDomain/StgDomain.h>
#include <StgFEM/StgFEM.h>
#include <PICellerator/PICellerator.h>
#include <Underworld/Underworld.h>
#include <Solvers/Solvers.h>

#include "petsccompat.h"

#include "TestKSP.h"

#define KSPTEST         "test"

EXTERN_C_BEGIN
EXTERN PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_TEST(KSP);
EXTERN_C_END

#undef __FUNCT__
#define __FUNCT__ "KSPRegisterTEST"
PetscErrorCode PETSCKSP_DLLEXPORT KSPRegisterTEST(const char path[])
{
    PetscErrorCode ierr;

    PetscFunctionBegin;
    ierr = Stg_KSPRegister(KSPTEST, path, "KSPCreate_TEST", KSPCreate_TEST );CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "KSPSolve_TEST"
PetscErrorCode  KSPSolve_TEST(KSP ksp)
{
    PetscPrintf( PETSC_COMM_WORLD, "----- In KSP Solver %s\n",__func__);
    PetscFunctionReturn(0);

}

#undef __FUNCT__
#define __FUNCT__ "KSPDestroy_TEST"
PetscErrorCode KSPDestroy_TEST(KSP ksp)
{
    PetscErrorCode ierr;
    PetscFunctionBegin;
    ierr = KSPDefaultDestroy(ksp);CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "KSPSetFromOptions_TEST"
PetscErrorCode KSPSetFromOptions_TEST(KSP ksp)
{

    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "KSPView_TEST"
PetscErrorCode KSPView_TEST(KSP ksp,PetscViewer viewer)
{
    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "KSPSetUp_TEST"
PetscErrorCode KSPSetUp_TEST(KSP ksp)
{
    PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "KSPCreate_TEST"
PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_TEST(KSP ksp)
{
    PetscErrorCode ierr;
    PetscFunctionBegin;

    #if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >= 2 ) )
    ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr);
    #endif
    /*
       Sets the functions that are associated with this data structure
       (in C++ this is the same as defining virtual functions)
    */
    ksp->ops->setup                = KSPSetUp_TEST;
    ksp->ops->solve                = KSPSolve_TEST;
    ksp->ops->destroy              = KSPDestroy_TEST;
    ksp->ops->view                 = KSPView_TEST;
    ksp->ops->setfromoptions       = KSPSetFromOptions_TEST;
    ksp->ops->buildsolution        = KSPDefaultBuildSolution;
    ksp->ops->buildresidual        = KSPDefaultBuildResidual;

    PetscFunctionReturn(0);
}
EXTERN_C_END


text:
/* siman/siman.c
 *
 * Copyright (C) 2007 Brian Gough
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Mark Galassi
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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 <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

#include <gsl/gsl_machine.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_siman.h>

static inline double
boltzmann(double E, double new_E, double T, gsl_siman_params_t *params)
{
  double x = -(new_E - E) / (params->k * T);
  /* avoid underflow errors for large uphill steps */
  return (x < GSL_LOG_DBL_MIN) ? 0.0 : exp(x);
}

static inline void
copy_state(void *src, void *dst, size_t size, gsl_siman_copy_t copyfunc)
{
  if (copyfunc) {
    copyfunc(src, dst);
  } else {
    memcpy(dst, src, size);
  }
}

/* implementation of a basic simulated annealing algorithm */

void
gsl_siman_solve (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef,
                 gsl_siman_step_t take_step,
                 gsl_siman_metric_t distance,
                 gsl_siman_print_t print_position,
                 gsl_siman_copy_t copyfunc,
                 gsl_siman_copy_construct_t copy_constructor,
                 gsl_siman_destroy_t destructor,
                 size_t element_size,
                 gsl_siman_params_t params)
{
  void *x, *new_x, *best_x;
  double E, new_E, best_E;
  int i;
  double T, T_factor;
  int n_evals = 1, n_iter = 0, n_accepts, n_rejects, n_eless;

  /* this function requires that either the dynamic functions (copy,
     copy_constructor and destrcutor) are passed, or that an element
     size is given */
  assert((copyfunc != NULL && copy_constructor != NULL && destructor != NULL)
         || (element_size != 0));

  distance = 0 ; /* This parameter is not currently used */
  E = Ef(x0_p);

  if (copyfunc) {
    x = copy_constructor(x0_p);
    new_x = copy_constructor(x0_p);
    best_x = copy_constructor(x0_p);
  } else {
    x = (void *) malloc (element_size);
    memcpy (x, x0_p, element_size);
    new_x = (void *) malloc (element_size);
    best_x =  (void *) malloc (element_size);
    memcpy (best_x, x0_p, element_size);
  }

  best_E = E;

  T = params.t_initial;
  T_factor = 1.0 / params.mu_t;

  if (print_position) {
    printf ("#-iter  #-evals   temperature     position   energy\n");
  }

  while (1) {

    n_accepts = 0;
    n_rejects = 0;
    n_eless = 0;

    for (i = 0; i < params.iters_fixed_T; ++i) {

      copy_state(x, new_x, element_size, copyfunc);

      take_step (r, new_x, params.step_size);
      new_E = Ef (new_x);

      if(new_E <= best_E){
        if (copyfunc) {
          copyfunc(new_x,best_x);
        } else {
          memcpy (best_x, new_x, element_size);
        }
        best_E=new_E;
      }

      ++n_evals;                /* keep track of Ef() evaluations */
      /* now take the crucial step: see if the new point is accepted
         or not, as determined by the boltzmann probability */
      if (new_E < E) {

	if (new_E < best_E) {
	  copy_state(new_x, best_x, element_size, copyfunc);
	  best_E = new_E;
	}

        /* yay! take a step */
	copy_state(new_x, x, element_size, copyfunc);
        E = new_E;
        ++n_eless;

      } else if (gsl_rng_uniform(r) < boltzmann(E, new_E, T, &params)) {
        /* yay! take a step */
	copy_state(new_x, x, element_size, copyfunc);
        E = new_E;
        ++n_accepts;

      } else {
        ++n_rejects;
      }
    }

    if (print_position) {
      /* see if we need to print stuff as we go */
      /*       printf("%5d %12g %5d %3d %3d %3d", n_iter, T, n_evals, */
      /*           100*n_eless/n_steps, 100*n_accepts/n_steps, */
      /*           100*n_rejects/n_steps); */
      printf ("%5d   %7d  %12g", n_iter, n_evals, T);
      print_position (x);
      printf ("  %12g  %12g\n", E, best_E);
    }

    /* apply the cooling schedule to the temperature */
    /* FIXME: I should also introduce a cooling schedule for the iters */
    T *= T_factor;
    ++n_iter;
    if (T < params.t_min) {
      break;
    }
  }

  /* at the end, copy the result onto the initial point, so we pass it
     back to the caller */
  copy_state(best_x, x0_p, element_size, copyfunc);

  if (copyfunc) {
    destructor(x);
    destructor(new_x);
    destructor(best_x);
  } else {
    free (x);
    free (new_x);
    free (best_x);
  }
}

/* implementation of a simulated annealing algorithm with many tries */

void
gsl_siman_solve_many (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef,
                      gsl_siman_step_t take_step,
                      gsl_siman_metric_t distance,
                      gsl_siman_print_t print_position,
                      size_t element_size,
                      gsl_siman_params_t params)
{
  /* the new set of trial points, and their energies and probabilities */
  void *x, *new_x;
  double *energies, *probs, *sum_probs;
  double Ex;                    /* energy of the chosen point */
  double T, T_factor;           /* the temperature and a step multiplier */
  int i;
  double u;                     /* throw the die to choose a new "x" */
  int n_iter;

  if (print_position) {
    printf ("#-iter    temperature       position");
    printf ("         delta_pos        energy\n");
  }

  x = (void *) malloc (params.n_tries * element_size);
  new_x = (void *) malloc (params.n_tries * element_size);
  energies = (double *) malloc (params.n_tries * sizeof (double));
  probs = (double *) malloc (params.n_tries * sizeof (double));
  sum_probs = (double *) malloc (params.n_tries * sizeof (double));

  T = params.t_initial;
  T_factor = 1.0 / params.mu_t;

  memcpy (x, x0_p, element_size);

  n_iter = 0;
  while (1)
    {
      Ex = Ef (x);
      for (i = 0; i < params.n_tries - 1; ++i)
        {                       /* only go to N_TRIES-2 */
          /* center the new_x[] around x, then pass it to take_step() */
          sum_probs[i] = 0;
          memcpy ((char *)new_x + i * element_size, x, element_size);
          take_step (r, (char *)new_x + i * element_size, params.step_size);
          energies[i] = Ef ((char *)new_x + i * element_size);
          probs[i] = boltzmann(Ex, energies[i], T, &params);
        }
      /* now add in the old value of "x", so it is a contendor */
      memcpy ((char *)new_x + (params.n_tries - 1) * element_size, x, element_size);
      energies[params.n_tries - 1] = Ex;
      probs[params.n_tries - 1] = boltzmann(Ex, energies[i], T, &params);

      /* now throw biased die to see which new_x[i] we choose */
      sum_probs[0] = probs[0];
      for (i = 1; i < params.n_tries; ++i)
        {
          sum_probs[i] = sum_probs[i - 1] + probs[i];
        }
      u = gsl_rng_uniform (r) * sum_probs[params.n_tries - 1];
      for (i = 0; i < params.n_tries; ++i)
        {
          if (u < sum_probs[i])
            {
              memcpy (x, (char *) new_x + i * element_size, element_size);
              break;
            }
        }
      if (print_position)
        {
          printf ("%5d\t%12g\t", n_iter, T);
          print_position (x);
          printf ("\t%12g\t%12g\n", distance (x, x0_p), Ex);
        }
      T *= T_factor;
      ++n_iter;
      if (T < params.t_min)
	{
	  break;
        }
    }

  /* now return the value via x0_p */
  memcpy (x0_p, x, element_size);

  /*  printf("the result is: %g (E=%g)\n", x, Ex); */

  free (x);
  free (new_x);
  free (energies);
  free (probs);
  free (sum_probs);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef DOCH
#define DOCH

#include <gsl/gsl_vector.h>
#include <assert.h>
#include "utils.h"
#include "topic.h"
#include "typedefs.h"

#define OFFSET 0
#define MH_GEM_STDEV 0.05
#define MH_GEM_MEAN_STDEV 0.05
#define MH_GEM_STDEV 0.05

/*
 * resample the levels of a document
 *
 */

void doc_sample_levels(doc* d, short do_permute, short do_remove);

/*
 * update a level count
 *
 */

void doc_update_level(doc* d, int l, double update);

/*
 * read corpus from data
 *
 */

void read_corpus(char* filename, corpus* c, int depth);

/*
 * allocate a new corpus
 *
 */

corpus* corpus_new(double mean, double scale);

/*
 * score the corpus
 *
 */

double gem_score(corpus* corp);

/*
 * GEM MH updates
 *
 */

void corpus_mh_update_gem(corpus* corp);
void corpus_mh_update_gem_mean(corpus* corp);
void corpus_mh_update_gem_scale(corpus* corp);

void compute_log_p_level(doc* d, double gem_mean, double gem_scale);

/*
 * write the document clustering to a file
 *
 */

void write_corpus_assignment(corpus* corp, FILE* file);
void write_corpus_levels(corpus* corp, FILE* file);

// free a corpus

void free_corpus(corpus* corp);
void free_doc(doc* d);

#endif

max_stars_repo_name: thatindiandude/PGM-Project

max_issues_repo_name: thatindiandude/PGM-Project

max_forks_repo_name: thatindiandude/PGM-Project


text:
#ifndef __GSL_VERSION_H__
#define __GSL_VERSION_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

#define GSL_VERSION "2.2.1"
#define GSL_MAJOR_VERSION 2
#define GSL_MINOR_VERSION 2

GSL_VAR const char * gsl_version;

__END_DECLS

#endif /* __GSL_VERSION_H__ */

max_stars_repo_name: Tjoppen/fmigo

max_issues_repo_name: Tjoppen/fmigo

max_forks_repo_name: Tjoppen/fmigo


text:
/* Copyright (c) Kyrylo Polezhaiev and contributors. All rights reserved.
   Released under the MIT license. See LICENSE file in the project root for full license information. */

#include "leopenblas.h"
#include <assert.h>
#include <le/le.h>
#include <le/tensors/letensor-imp.h>
#include <cblas.h>

LeTensor *
le_openblas_matrix_new_product(const LeTensor *a, bool transpose_a, const LeTensor *b, bool transpose_b)
{
    /// @todo: Take stride into account
    assert(a->element_type == LE_TYPE_FLOAT32);
    assert(b->element_type == LE_TYPE_FLOAT32);
    assert(a->shape->num_dimensions == 2);
    assert(b->shape->num_dimensions == 2);
    assert(le_tensor_contiguous(a));
    assert(le_tensor_contiguous(b));

    unsigned size_a = transpose_a ? a->shape->sizes[0] : a->shape->sizes[1];
    unsigned size_b = transpose_b ? b->shape->sizes[1] : b->shape->sizes[0];
    assert(size_a == size_b);

    unsigned c_height = transpose_a ? a->shape->sizes[1] : a->shape->sizes[0];
    unsigned c_width = transpose_b ? b->shape->sizes[0] : b->shape->sizes[1];

    LeTensor *c = le_matrix_new_uninitialized(LE_TYPE_FLOAT32, c_height, c_width);

    cblas_sgemm(CblasRowMajor,
                transpose_a ? CblasTrans : CblasNoTrans,
                transpose_b ? CblasTrans : CblasNoTrans,
                c_height, c_width, size_a,
                1.0f,
                a->data, a->shape->sizes[1],
                b->data, b->shape->sizes[1],
                0.0f,
                c->data, c->shape->sizes[1]);

    return c;
}

float
le_openblas_dot_product(const LeTensor *a, const LeTensor *b)
{
    assert(a->element_type == LE_TYPE_FLOAT32);
    assert(b->element_type == LE_TYPE_FLOAT32);
    assert(a->shape->num_dimensions == 2);
    assert(b->shape->num_dimensions == 2);
    /** @todo: Test results against transposed a multiplied by b */
    assert(a->shape->sizes[0] == b->shape->sizes[0]);
    assert(a->shape->sizes[1] == 1);
    assert(b->shape->sizes[1] == 1);

    return cblas_sdot(a->shape->sizes[0], a->data, a->stride, b->data, b->stride);
}

max_stars_repo_name: neonkingfr/le

max_issues_repo_name: neonkingfr/le

max_forks_repo_name: neonkingfr/le


text:
//          Copyright Jean Pierre Cimalando 2019.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE.md or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

#pragma once
#include <gsl/gsl>
#include <string>

struct fmidi_smf;

struct SMF_Encoding_Detector {
public:
    void scan(const fmidi_smf &smf);

    std::string general_encoding() const;
    std::string encoding_for_text(gsl::cstring_span input) const;

    std::string decode_to_utf8(gsl::cstring_span input) const;

private:
    static gsl::cstring_span encoding_from_marker(gsl::cstring_span input);
    static gsl::cstring_span strip_encoding_marker(gsl::cstring_span text, gsl::cstring_span enc);

private:
    std::string encoding_;
};

max_stars_repo_name: jpcima/smf-dsp

max_issues_repo_name: jpcima/smf-dsp

max_forks_repo_name: jpcima/smf-dsp


text:
/* specfunc/legendre_P.c
 *
 * Copyright (C) 2009-2013 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_legendre.h>

/*
 * The routines in this module compute associated Legendre functions
 * (ALFs) up to order and degree 2700, using the method described
 * in
 *
 * [1] S. A. Holmes and W. E. Featherstone, A unified approach
 *     to the Clenshaw summation and the recursive computation of very
 *     high degree and order normalised associated Legendre functions,
 *     Journal of Geodesy, 76, pg. 279-299, 2002.
 *
 * Further information on ALFs can be found in
 *
 * [2] Abramowitz and Stegun, Handbook of Mathematical Functions,
 *     Chapter 8, 1972.
 */

static void legendre_sqrts(const size_t lmax, double *array);

#define LEGENDRE
#include "legendre_source.c"
#undef LEGENDRE

#define LEGENDRE_DERIV
#include "legendre_source.c"
#undef LEGENDRE_DERIV

#define LEGENDRE_DERIV_ALT
#include "legendre_source.c"
#undef LEGENDRE_DERIV_ALT

#define LEGENDRE_DERIV2
#include "legendre_source.c"
#undef LEGENDRE_DERIV2

#define LEGENDRE_DERIV2_ALT
#include "legendre_source.c"
#undef LEGENDRE_DERIV2_ALT

/* number of P_{lm} functions for a given lmax */
size_t
gsl_sf_legendre_nlm(const size_t lmax)
{
  return ((lmax + 1) * (lmax + 2) / 2);
}

/*
gsl_sf_legendre_array_n()
  This routine returns the minimum result_array[] size needed
for a given lmax
*/

size_t
gsl_sf_legendre_array_n(const size_t lmax)
{
  size_t nlm = gsl_sf_legendre_nlm(lmax);
  size_t nsqrt = 2 * lmax + 2; /* extra room to precompute sqrt factors */

  return (nlm + nsqrt);
} /* gsl_sf_legendre_array_n() */

/*
gsl_sf_legendre_array_index()
This routine computes the index into a result_array[] corresponding
to a given (l,m)
*/

size_t
gsl_sf_legendre_array_index(const size_t l, const size_t m)
{
  return (l * (l + 1) / 2 + m);
} /* gsl_sf_legendre_array_index() */

/*********************************************************
 *                 INTERNAL ROUTINES                     *
 *********************************************************/

/*
legendre_sqrts()
  Precompute square root factors needed for Legendre recurrence.
On output, array[i] = sqrt(i)
*/

static void
legendre_sqrts(const size_t lmax, double *array)
{
  size_t l;
  for (l = 0; l <= 2 * lmax + 1; ++l)
    array[l] = sqrt((double) l);
}


text:
#pragma once

#ifndef _NOSNOPT
#include <fstream>
#include <iostream>
#include "BasicError.h"
#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#else
#include "FakeGSL.h"
#endif
#include "OptSolver.h"
#include "snopt.hh"
#include "snoptProblem.hh"

using namespace std;

class SnoptSolver : public OptSolver {
	snoptProblem snoptProb;

	integer n; integer neF; integer lenA; integer lenG;

	integer *iAfun;
	integer *jAvar;
	doublereal *A;

	integer *iGfun;
	integer *jGvar;

	doublereal *x;
	doublereal *xlow;
	doublereal *xupp;
	doublereal *xmul;
	integer *xstate;

	doublereal *F;
	doublereal *Flow;
	doublereal *Fupp;
	doublereal *Fmul;
	integer *Fstate;

	integer nxnames;
	integer nFnames = 1;
	char *xnames;
	char *Fnames;

	gsl_vector* result;
	double objectiveVal;

	integer Cold = 0, Basis = 1, Warm = 2;
public:
	SnoptSolver(integer n_, integer neF_, integer lenA_): n(n_), neF(neF_), lenA(lenA_) {
		iAfun = new integer[lenA];
		jAvar = new integer[lenA];
		A = new doublereal[lenA];

		lenG = n*neF;
		iGfun = new integer[lenG];
		jGvar = new integer[lenG];

		x = new doublereal[n];
		xlow = new doublereal[n];
		xupp = new doublereal[n];
		xmul = new doublereal[n];
		xstate = new integer[n];

		F = new doublereal[neF];
		Flow = new doublereal[neF];
		Fupp = new doublereal[neF];
		Fmul = new doublereal[neF];
		Fstate = new integer[neF];

		nxnames = 1;
		nFnames = 1;
		xnames = new char[nxnames*8];
		Fnames = new char[nFnames*8];

		snoptProb.setXNames(xnames, nxnames);
		snoptProb.setFNames(Fnames, nFnames);

		result = gsl_vector_alloc(n);
	}

	~SnoptSolver() {
		delete []iAfun; delete []jAvar; delete []A;
		delete []iGfun; delete []jGvar;
		delete []x; delete []xlow; delete []xupp;
		delete []xmul; delete []xstate;
		delete []F; delete []Flow; delete []Fupp;
		delete []Fmul; delete []Fstate;
		delete []xnames; delete []Fnames;
	}

	virtual void init(char* workspace, integer nef_, DFT df, integer ObjRow, doublereal ObjAdd, doublereal *xlow_, doublereal *xupp_, doublereal *Flow_, doublereal *Fupp_);
	bool optimize(gsl_vector* initState, bool suppressPrint = false);

	virtual gsl_vector* getResults() {
		return result;
	}

	virtual double getObjectiveVal() {
		return objectiveVal;
	}

};

#endif

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


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

#ifndef __GSL_VECTOR_ULONG_H__
#define __GSL_VECTOR_ULONG_H__

#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_block_ulong.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  unsigned long *data;
  gsl_block_ulong *block;
  int owner;
}
gsl_vector_ulong;

typedef struct
{
  gsl_vector_ulong vector;
} _gsl_vector_ulong_view;

typedef _gsl_vector_ulong_view gsl_vector_ulong_view;

typedef struct
{
  gsl_vector_ulong vector;
} _gsl_vector_ulong_const_view;

typedef const _gsl_vector_ulong_const_view gsl_vector_ulong_const_view;

/* Allocation */

gsl_vector_ulong *gsl_vector_ulong_alloc (const size_t n);
gsl_vector_ulong *gsl_vector_ulong_calloc (const size_t n);

gsl_vector_ulong *gsl_vector_ulong_alloc_from_block (gsl_block_ulong * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

gsl_vector_ulong *gsl_vector_ulong_alloc_from_vector (gsl_vector_ulong * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

void gsl_vector_ulong_free (gsl_vector_ulong * v);

/* Views */

_gsl_vector_ulong_view
gsl_vector_ulong_view_array (unsigned long *v, size_t n);

_gsl_vector_ulong_view
gsl_vector_ulong_view_array_with_stride (unsigned long *base,
                                         size_t stride,
                                         size_t n);

_gsl_vector_ulong_const_view
gsl_vector_ulong_const_view_array (const unsigned long *v, size_t n);

_gsl_vector_ulong_const_view
gsl_vector_ulong_const_view_array_with_stride (const unsigned long *base,
                                               size_t stride,
                                               size_t n);

_gsl_vector_ulong_view
gsl_vector_ulong_subvector (gsl_vector_ulong *v,
                            size_t i,
                            size_t n);

_gsl_vector_ulong_view
gsl_vector_ulong_subvector_with_stride (gsl_vector_ulong *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

_gsl_vector_ulong_const_view
gsl_vector_ulong_const_subvector (const gsl_vector_ulong *v,
                                  size_t i,
                                  size_t n);

_gsl_vector_ulong_const_view
gsl_vector_ulong_const_subvector_with_stride (const gsl_vector_ulong *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

unsigned long gsl_vector_ulong_get (const gsl_vector_ulong * v, const size_t i);
void gsl_vector_ulong_set (gsl_vector_ulong * v, const size_t i, unsigned long x);

unsigned long *gsl_vector_ulong_ptr (gsl_vector_ulong * v, const size_t i);
const unsigned long *gsl_vector_ulong_const_ptr (const gsl_vector_ulong * v, const size_t i);

void gsl_vector_ulong_set_zero (gsl_vector_ulong * v);
void gsl_vector_ulong_set_all (gsl_vector_ulong * v, unsigned long x);
int gsl_vector_ulong_set_basis (gsl_vector_ulong * v, size_t i);

int gsl_vector_ulong_fread (FILE * stream, gsl_vector_ulong * v);
int gsl_vector_ulong_fwrite (FILE * stream, const gsl_vector_ulong * v);
int gsl_vector_ulong_fscanf (FILE * stream, gsl_vector_ulong * v);
int gsl_vector_ulong_fprintf (FILE * stream, const gsl_vector_ulong * v,
			      const char *format);

int gsl_vector_ulong_memcpy (gsl_vector_ulong * dest, const gsl_vector_ulong * src);

int gsl_vector_ulong_reverse (gsl_vector_ulong * v);

int gsl_vector_ulong_swap (gsl_vector_ulong * v, gsl_vector_ulong * w);
int gsl_vector_ulong_swap_elements (gsl_vector_ulong * v, const size_t i, const size_t j);

unsigned long gsl_vector_ulong_max (const gsl_vector_ulong * v);
unsigned long gsl_vector_ulong_min (const gsl_vector_ulong * v);
void gsl_vector_ulong_minmax (const gsl_vector_ulong * v, unsigned long * min_out, unsigned long * max_out);

size_t gsl_vector_ulong_max_index (const gsl_vector_ulong * v);
size_t gsl_vector_ulong_min_index (const gsl_vector_ulong * v);
void gsl_vector_ulong_minmax_index (const gsl_vector_ulong * v, size_t * imin, size_t * imax);

int gsl_vector_ulong_add (gsl_vector_ulong * a, const gsl_vector_ulong * b);
int gsl_vector_ulong_sub (gsl_vector_ulong * a, const gsl_vector_ulong * b);
int gsl_vector_ulong_mul (gsl_vector_ulong * a, const gsl_vector_ulong * b);
int gsl_vector_ulong_div (gsl_vector_ulong * a, const gsl_vector_ulong * b);
int gsl_vector_ulong_scale (gsl_vector_ulong * a, const double x);
int gsl_vector_ulong_add_constant (gsl_vector_ulong * a, const double x);

int gsl_vector_ulong_isnull (const gsl_vector_ulong * v);

extern int gsl_check_range;

#ifdef HAVE_INLINE

extern inline
unsigned long
gsl_vector_ulong_get (const gsl_vector_ulong * v, const size_t i)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= v->size)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

extern inline
void
gsl_vector_ulong_set (gsl_vector_ulong * v, const size_t i, unsigned long x)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= v->size)
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

extern inline
unsigned long *
gsl_vector_ulong_ptr (gsl_vector_ulong * v, const size_t i)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (unsigned long *) (v->data + i * v->stride);
}

extern inline
const unsigned long *
gsl_vector_ulong_const_ptr (const gsl_vector_ulong * v, const size_t i)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const unsigned long *) (v->data + i * v->stride);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_ULONG_H__ */