text:
/*
 Copyright (C) 2006-2007 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*/

#include "util.h"
#include "xc.h"

void func(double *x, int n, void *ex)
{
  int i;
  for(i=0; i<n;i++)
    x[i] = cos(x[i]);
}

void test_integration()
{
  double a, b, result;

  for(b=1e-8; b<5; b+=0.001){
    result = xc_integrate(func, NULL, a, b);
    printf("%lf %lf\n", b, result);
  }
}

#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_sf_bessel.h>

void test_expi()
{
  double b, r1, r2, r3, r4;
  int n;

  n = 1;
  for(b=1e-3; b<50; b+=0.01){
    r1 = xc_bessel_K0(b);
    //r2 = gsl_sf_bessel_K0(b);
    //r3 = bessi1(b);
    printf("%5.3lf %12.10lf %12.10lf %12.10lf\n", b, r1, r2, r3);
  }
}

void test_lda()
{
  xc_func_type l1, l2, l3;
  int i;

  //xc_func_init(&l1, XC_LDA_XC_KSDT2, XC_POLARIZED);
  //xc_func_init(&l1, XC_LDA_C_PW, XC_UNPOLARIZED);
  xc_func_init(&l3, XC_LDA_X, XC_UNPOLARIZED);

  //xc_lda_xc_ksdt_set_params(&l1, 1000.0);
  //xc_lda_c_1d_csc_set_params(&l2, 1, 1.0);

  for(i=0; i<=1000; i++){
    double dens, rs, zeta, rho[2];
    double ec1, vc1[2], fxc1[3], kxc1[4];
    double ec2, vc2[2], fxc2[3], kxc2[4];
    double ec3, vc3[2], fxc3[3], kxc3[4];

    //rs   = 0.5 + i/500.0;
    //zeta = -1.0 + 2.0*i/1000000000.0;

    //dens = 1.0/(4.0/3.0*M_PI*pow(rs,3)); /* 3D */
    //dens = 1.0/(2.0*rs); /* 1D */

    //rho[0] = dens*(1.0 + zeta)/2.0;
    //rho[1] = dens*(1.0 - zeta)/2.0;

    rho[0] = 0.05 + i/1000.0;
    rho[1] = 0.001;

    //rho[0] = 1.0/(2.0*rs);
    //rho[1] = 0.0;

    //dens = rho[0] + rho[1];

    //xc_lda_exc_vxc_fxc(&l1, 1, rho, &ec1, vc1, NULL, NULL);
    //xc_lda(&l2, 1, rho, &ec2, vc2, NULL, NULL);
    //xc_lda_fxc_fd(&l2, rho, fxc2);
    //xc_lda_kxc_fd(&l2, rho, kxc2);

    //rho[0] = dens; rho[1] = 0.0;
    //xc_lda(&l3, rho, &ec3, vc3, fxc3, kxc3);

    // printf("%e\t%e\t%e\n", dens, (fxc1[0]+2.0*fxc1[1]+fxc1[2])/4.0, fxc3[0]);
    // printf("%e\t%e\t%e\n", dens, (kxc1[0]+3.0*kxc1[1]+3.0*kxc1[2]+kxc1[3])/8.0, kxc3[0]);

    printf("%e\t%e\t%e\n", rho[0], (rho[0] + rho[1])*ec1, vc1[0]);
  }
}

void test_ak13()
{
  xc_func_type gga;
  double beta = 0.13;
  double x, rho[2] = {0.0, 0.0}, sigma[3] = {0.0, 0.0, 0.0}, zk, vrho[2], vsigma[3];
  double tmp1, tmp2;
  int i;

  xc_func_init(&gga,  XC_GGA_X_AK13,  XC_POLARIZED);

  for(i=0; i<=10000; i++){

    x = 500.0*i/(10000.0);
    rho[0]   = 0.12*exp(-beta * x);
    sigma[0] = 0.12*0.12*beta*beta * rho[0]*rho[0];

    xc_gga_exc_vxc(&gga,  1, rho, sigma, &zk, vrho, vsigma);

    tmp2 = 1.74959015598863046792081721182*beta*x/3.0- 1.62613336586517367779736042170*log(x);
    fprintf(stderr, "%16.10lf\t%16.10lf\t%16.10lf\t%16.10lf\n", x, vrho[0], vsigma[0]*sqrt(sigma[0]),
	    -X_FACTOR_C*X2S*tmp2/2.0);
  }

}

void test_enhance()
{
  double x, f, dfdx, d2fdx2, d3fdx3;

  xc_func_type gga1, gga2, gga3, gga4;

  xc_func_init(&gga1,  XC_GGA_X_B88, XC_POLARIZED);
  xc_func_init(&gga2, XC_GGA_X_GG99, XC_POLARIZED);

  for(x=0.01; x<200; x+=0.01){
    //printf("%le", x);
    //xc_gga_x_b88_enhance(&gga1, 1, x, &f, &dfdx, &d2fdx2, &d3fdx3);
    //printf("\t%le", -f*X_FACTOR_C);
    //xc_gga_x_gg99_enhance(&gga2, 1, x, &f, &dfdx, &d2fdx2, &d3fdx3);
    printf("%le\t%le\t%le\n", x, f, dfdx);
  }
}

void test_gga()
{
  xc_func_type gga, ggap;
  int i, npoints = 1;
  double *rho, *sigma;
  double *zk, zkp, *vrho, vrhop[2], *vsigma, vsigmap[3];
  double *v2rho2, *v2rhosigma, *v2sigma2;
  double *v3rho3, *v3rho2sigma, *v3rhosigma2, *v3sigma3;

  rho         = libxc_malloc( 2*npoints*sizeof(double));
  sigma       = libxc_malloc( 3*npoints*sizeof(double));
  zk          = libxc_malloc( 1*npoints*sizeof(double));
  vrho        = libxc_malloc( 2*npoints*sizeof(double));
  vsigma      = libxc_malloc( 3*npoints*sizeof(double));
  v2rho2      = libxc_malloc( 3*npoints*sizeof(double));
  v2rhosigma  = libxc_malloc( 6*npoints*sizeof(double));
  v2sigma2    = libxc_malloc( 6*npoints*sizeof(double));
  v3rho3      = libxc_malloc( 4*npoints*sizeof(double));
  v3rho2sigma = libxc_malloc( 9*npoints*sizeof(double));
  v3rhosigma2 = libxc_malloc(12*npoints*sizeof(double));
  v3sigma3    = libxc_malloc(10*npoints*sizeof(double));

  xc_func_init(&gga,  XC_GGA_C_PW91,  XC_POLARIZED);

  /*
  for(i=1; i<=10000; i++){
    double x = 30.0*i/(10000.0), f, df, d2f, d3f;
    double c2 = 10.0/81.0;

    xc_gga_x_b88_enhance(&gga, 3, x, &f, &df, &d2f, &d3f);
    printf("%20.14e %20.14e", X2S*X2S*x*x, (f-1.0)/(c2*X2S*X2S*x*x));
    xc_gga_x_ev93_enhance(&gga, 3, x, &f, &df, &d2f, &d3f);
    printf(" %20.14e\n", (f-1.0)/(c2*X2S*X2S*x*x));
  }
  exit(0);
  */

   for(i=0; i<=10000; i++){
     rho[0]   = .01;
     rho[1]   = 0.2;
     sigma[0] = 0.1;
     sigma[1] = 0.00002;
     sigma[2] = 0.5 + i/1000.0;

     //xc_gga(&gga, 1, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2, NULL, v3rho2sigma, v3rhosigma2, v3sigma3);

     fprintf(stderr, "%16.10le\t%16.10le\t%16.10le\n", sigma[2], vsigma[2], v2sigma2[5]);
   }

  /*
  for(i=0; i<npoints; i++){
    rho[2*i + 0]   = 0.048 + i/10000.0;
    rho[2*i + 1]   = 0.025;
    sigma[3*i + 0] = 0.0046;
    sigma[3*i + 1] = 0.0044;
    sigma[3*i + 2] = 0.0041;
  }

  xc_gga(&gga,  npoints, rho, sigma, zk,  vrho,  vsigma,  v2rho2,  v2rhosigma,  v2sigma2, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3);

  for(i=0; i<npoints; i++){
    fprintf(stderr, "%16.10lf\t%16.10lf\t%16.10lf\n", rho[2*i + 0], vrho[2*i + 0], v2rho2[3*i + 0]);
  }
  */

  xc_func_end(&gga);

  libxc_free(rho); libxc_free(sigma);
  libxc_free(zk); libxc_free(vrho); libxc_free(vsigma);
  libxc_free(v2rho2); libxc_free(v2rhosigma); libxc_free(v2sigma2);
  libxc_free(v3rho3); libxc_free(v3rho2sigma); libxc_free(v3rhosigma2); libxc_free(v3sigma3);
}

void test_mgga()
{
  xc_func_type mgga1, mgga2;
  int i;

  xc_func_init(&mgga1, XC_GGA_C_LYP, XC_POLARIZED);
  xc_func_init(&mgga2, XC_MGGA_X_SCAN, XC_POLARIZED);
  //xc_mgga_c_tpss_init(tpss2.mgga);

  for(i=0; i<=1000; i++){
    double rho[2], sigma[3], tau[2], lapl[2];
    double zk,   vrho[2],  vsigma[3],  vtau[2],  vlapl[2];
    double zk2, vrho2[2], vsigma2[3], vtau2[2], vlapl2[2];
    double v2rho2[3], v2sigma2[6], v2lapl2[3], v2tau2[3];
    double v2rhosigma[6], v2rholapl[3], v2rhotau[3];
    double v2sigmalapl[6], v2sigmatau[6], v2lapltau[3];

    rho[0]   = 0.3;
    rho[1]   = 0.4;
    sigma[0] = 0.2032882206468622;
    sigma[1] = 0.11;
    sigma[2] = 0.7 + i/10.0;
    tau[0]   = 1.0;
    tau[1]   = 0.15;
    lapl[0]  = -0.1518421131246519;
    lapl[1]  = 0.12;

    //xc_mgga(&mgga1, 1, rho, sigma, lapl, tau,
    //	     &zk,  vrho, vsigma, vlapl, vtau,
    //	     v2rho2, v2sigma2, v2lapl2, v2tau2, v2rhosigma, v2rholapl, v2rhotau,
    //	     v2sigmalapl, v2sigmatau, v2lapltau);
    //xc_mgga(&mgga2, 1, rho, sigma, lapl, tau,
    //	     &zk2,  vrho2, vsigma2, vlapl2, vtau2,
    // 	     NULL, v2sigma2, v2lapl2, v2tau2, v2rhosigma, v2rholapl, v2rhotau,
    //	     v2sigmalapl, v2sigmatau, v2lapltau);
    //xc_mgga_exc(&mgga2, 1, rho, sigma, lapl, tau,
    //		 &zk2);
    xc_gga_exc_vxc(&mgga1, 1, rho, sigma,
		     &zk,  vrho, vsigma);
    xc_mgga_exc_vxc(&mgga2, 1, rho, sigma, lapl, tau,
		     &zk2,  vrho2, vsigma2, vlapl2, vtau2);

    fprintf(stderr, "%16.10lf\t%16.10lf\t%16.10lf\n", sigma[2], zk2*(rho[0]+rho[1]), vsigma2[2]);
  }

  xc_func_end(&mgga1);
  xc_func_end(&mgga2);
}

void test_neg_rho()
{
  xc_func_type func;
  double rho[5][2] = {
    {9.03897273e-06, -1.00463992e-06},
    {8.48383564e-06, -3.51231267e-07},
    {1.45740621e-08, -2.94546705e-09},
    {2.62778445e-07, -1.00191745e-07},
    {2.55745103e-06, -1.54789964e-06}
  };
  double sigma[5][3] = {
    {1.20122271e-08, 4.83240746e-09, 6.24774836e-09},
    {1.54146602e-07, 1.41584609e-07, 1.36663204e-07},
    {2.75312438e-08, 2.75224049e-08, 2.75135719e-08},
    {1.90251649e-07, 1.91241798e-07, 1.92240989e-07},
    {9.29562712e-09, 7.83940082e-09, 8.05714636e-09}
  };
  double vsigma[5][3];
  double zk[5], vk[5][2];
  int i, func_id;

  for(func_id=1; func_id<1000; func_id++){
    if(xc_func_init(&func, func_id, XC_POLARIZED) != 0) continue;
    if(func_id == XC_LDA_C_2D_PRM || func_id == XC_GGA_X_LB) goto end;

    printf("\n%s:\n", func.info->name);

    switch(func.info->family){
    case XC_FAMILY_LDA:
      xc_lda_exc_vxc(&func, 5, &rho[0][0], zk, &vk[0][0]);
      break;
    case XC_FAMILY_GGA:
      xc_gga_exc_vxc(&func, 5, &rho[0][0], &sigma[0][0], zk, &vk[0][0], &vsigma[0][0]);
      break;
    }

    switch(func.info->family){
    case XC_FAMILY_LDA:
      for(i=0; i<5; i+=1)
        printf("%.8e %.8e %.8e %.8e %.8e\n",
               rho[i][0], rho[i][1], zk[i], vk[i][0], vk[i][1]);
      break;
    case XC_FAMILY_GGA:
      for(i=0; i<5; i+=1)
        printf("%.8e %.8e %.8e %.8e %.8e %.8e %.8e %.8e %.8e %.8e %.8e\n",
               rho[i][0], rho[i][1], sigma[i][0], sigma[i][1], sigma[i][2],
               zk[i], vk[i][0], vk[i][1], vsigma[i][0], vsigma[i][1], vsigma[i][2]);
      break;
    }

  end:
    xc_func_end(&func);
  }
}

double xc_mgga_x_mbrxc_get_x(double Q);
double xc_mgga_x_br89_get_x(double Q);
void test_mbrxc()
{
  double Q, x, rhs, Q2;

  for(Q=-10.1e4; Q<10.1e4; Q+=.01e4){
  printf("Q = %lf ", Q);
  fflush(stdout);

  //x = xc_mgga_x_mbrxc_get_x(Q);
  //rhs = pow(1.0 + x, 5.0/3.0)*exp(-2.0*x/3.0)/(x - 3.0);
  //Q2 = pow(32.0*M_PI, 2.0/3.0)/(6.0*rhs);

  x = xc_mgga_x_br89_get_x(Q);
  rhs = x*exp(-2.0*x/3.0)/(x - 2.0);
  Q2 = 2.0/3.0*pow(M_PI, 2.0/3.0)/rhs;
  printf("Q2 = %lf x = %lf\n", Q2, x);
  }
}

int main()
{
  //test_expi();
  //test_integration();
  //test_neg_rho();

  //test_lda();
  //test_enhance();
  //test_gga();
  //test_ak13();
  //test_mgga();
  test_mbrxc();

  //printf("number = '%d'; key = '%s'", 25, xc_functional_get_name(25));

  return 0;
}


text:
/**
 * Set of routines to handle absolute flux calibration of a spectrum
 * structure.
 * The response curves can be in ASCII or (preferably) in FITS
 * binary format (first extension). ASCI format is simply three
 * space separated colums (WAVELENGTH, THROUGHPUT, ERROR). The FITS
 * binary table must contain three columns with names WAVELENGTH,
 * SENSITIVITY, ERROR. The throughput should be DN/s per Erg/cm^2/s/A.
 */

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

#include "aXe_grism.h"
#include "disp_conf.h"
#include "fringe_conf.h"
#include "spc_wl_calib.h"
#include "spc_resp.h"

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

/**
 * Function: get_response_function_from_FITS
 * This function reads a FITS table, reads the WAVELENGTH,
 * SENSITIVITY, and ERROR columns and returns a spectrum
 * structure which has been populated with this throughput.
 *
 * Parameters:
 * @param filename - a pointer to a char array containing the
 *                   name of an existing FITS file
 * @param hdunum   - the extension number to read (2 for STSDAS tables)
 *
 * Returns:
 * @return res     - a pointer to a newly allocated spectrum structure
 */
spectrum *get_response_function_from_FITS(char filename[], int hdunum)
{
  fitsfile *input;
  int f_status = 0;
  int hdutype, anynul;
  double *lambda, *resp, *error;
  double tmp;
  long numrows;
  int colnum, i;
  spectrum *res;

  //  Open the file for reading
  fits_open_file (&input, filename, READONLY, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "get_response_function_fromFITS: "
                   "Could not open" " file: %s",
                   filename);
    }
  /* Move to first hdu */
  fits_movabs_hdu (input, hdunum, &hdutype, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "get_response_function_fromFITS: "
                   "Could not read extention %d from file: %s",
                   hdunum, filename);
    }

  /* Get number of rows */
  fits_get_num_rows (input, &numrows, &f_status);
  if (f_status) {
    ffrprt (stderr, f_status);
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                 "get_response_function_fromFITS: "
                 "Could not determine the number of rows in"
                 " table %s",filename);
  }

  /* Allocate temporary memory space */
  lambda = (double *) malloc(numrows*sizeof(double));
  if (!lambda) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
                 "Memory allocation failed");
  }
  resp = (double *) malloc(numrows*sizeof(double));
  if (!resp) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
                 "Memory allocation failed");
  }
  error = (double *) malloc(numrows*sizeof(double));
  if (!error) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
                 "Memory allocation failed");
  }

  /**************************/
  /* Read the WAVELENGTH column */
  /**************************/
  /* Get column number */
  fits_get_colnum (input, CASEINSEN, "WAVELENGTH", &colnum, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "get_response_function_fromFITS: "
                   "Could not determine WAVELENGTH column number in "
                   " table %s",filename);
    }

  /* Read the data */
  fits_read_col (input, TDOUBLE, colnum, 1, 1, numrows, NULL, lambda,
                    &anynul, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "get_response_function_fromFITS: "
                   "Could not read content of WAVELENGTH column "
                   " from BINARY table %s",filename);
    }

  /**************************/
  /* Read the SENSITIVITY column */
  /**************************/
  /* Get column number */
  fits_get_colnum (input, CASEINSEN, "SENSITIVITY", &colnum, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "get_response_function_fromFITS: "
                   "Could not determine SENSITIVITY column number in "
                   " table %s",filename);
    }
  /* Read the data */
  fits_read_col (input, TDOUBLE, colnum, 1, 1, numrows, NULL, resp,
                    &anynul, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "get_response_function_fromFITS: "
                   "Could not read content of SENSITIVITY column "
                   " from BINARY table %s",filename);
    }
  /**************************/
  /* Read the ERROR column */
  /**************************/
  /* Get column number */
  fits_get_colnum (input, CASEINSEN, "ERROR", &colnum, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "get_response_function_fromFITS: "
                   "Could not determine ERROR column number in "
                   " table %s",filename);
    }
  /* Read the data */
  fits_read_col (input, TDOUBLE, colnum, 1, 1, numrows, NULL, error,
                    &anynul, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "get_response_function_fromFITS: "
                   "Could not read content of ERROR column "
                   " from BINARY table %s",filename);
    }

  /* Allocate memory */
  res = (spectrum *) malloc(sizeof(spectrum));
  if(!res) {
      aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
                   "Memory allocation failed");
  }
  res->spec = (spc_entry *) malloc(numrows*sizeof(spc_entry));
  if(!(res->spec)) {
      aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
                   "Memory allocation failed");
  }

  for (i=0;i<numrows;i++) {
    res->spec[i].lambda_mean = lambda[i];
    res->spec[i].flux = resp[i];
    res->spec[i].ferror = error[i];
  }

  res->spec_len=i;

  // define the minimum and maximum wavelength
  res->lambdamin=res->spec[0].lambda_mean;
  res->lambdamax=res->spec[numrows-1].lambda_mean;

  // check whether the wavelength is raising.
  // switch minimum and maximum for falling spectrum
  if (res->lambdamin > res->lambdamax)
    {
      tmp = res->lambdamin;
      res->lambdamin = res->lambdamax;
      res->lambdamax = tmp;
    }

  fits_close_file(input,&f_status);
  if (f_status) {
      aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
                   "Could not close %s",filename);
  }

  free(lambda);
  free(resp);
  free(error);

  return res;
}

/**
 * Function: apply_response_function
 * This function applies a thoughput curve to a given spectrum
 * by dividing the count attribute of the spectrum. The flux
 * and ferror attributes of the spectrum are populated.
 * gsl spline interpolation is used to compute the throughput
 * curve at the needed spectral wavelengths.
 *
 * Parameters:
 * @param spec    - a pointer to an existing spectrum array
 * @param resp    - a pointer to a spectrum structure containing
 *                  a throughput curve previously loaded
 * @param exptime - the exposure time
 * @param gain    - the gain
 *
 * Returns
 * @return -
 */
void
apply_response_function(spectrum *spec, spectrum *resp, const int quant_cont)
{
  gsl_interp_accel *acc1 = gsl_interp_accel_alloc();
  gsl_spline *spline1 = gsl_spline_alloc (gsl_interp_cspline, resp->spec_len);

  gsl_interp_accel *acc2 = gsl_interp_accel_alloc();
  gsl_spline *spline2 = gsl_spline_alloc (gsl_interp_cspline, resp->spec_len);

  double *x1,*y1, *x2, *y2;
  int i,j;
  double r1, r2, fr, fr1, fr2;
  double min_l, max_l;
  // test
  // int nguess;
  // double tval=0.0;

  if ((int) resp->spec_len <= 0){
    aXe_message(aXe_M_ERROR, __FILE__, __LINE__,
                "apply_response_function: response spec length invalid: %i", resp->spec_len);
      gsl_spline_free(spline1);
      gsl_interp_accel_free(acc1);
      gsl_spline_free(spline2);
      gsl_interp_accel_free(acc2);
      return;
  }
  if ( spec==NULL)
    {
      aXe_message (aXe_M_WARN4, __FILE__, __LINE__,
                   "apply_response_function: spectra empty.");
      gsl_spline_free(spline1);
      gsl_interp_accel_free(acc1);
      gsl_spline_free(spline2);
      gsl_interp_accel_free(acc2);
      return;
    }

  x1 = malloc((double)(resp->spec_len) *sizeof(double));
  y1 = malloc((double)(resp->spec_len) *sizeof(double));
  x2 = malloc((double)(resp->spec_len) *sizeof(double));
  y2 = malloc((double)(resp->spec_len) *sizeof(double));

  min_l = 1e32;
  max_l = -1e32;

  for (i=0;i<resp->spec_len; i++) {
    if (resp->spec[i].lambda_mean > max_l) {
     max_l = resp->spec[i].lambda_mean;
    }
    if (resp->spec[i].lambda_mean < min_l) {
     min_l = resp->spec[i].lambda_mean;
    }

    x1[i] = resp->spec[i].lambda_mean;
    y1[i] = resp->spec[i].flux - resp->spec[i].ferror;
    x2[i] = resp->spec[i].lambda_mean;
    y2[i] = resp->spec[i].flux + resp->spec[i].ferror;
  }

  gsl_spline_init (spline1, x1, y1, resp->spec_len);
  gsl_spline_init (spline2, x2, y2, resp->spec_len);
  // test
  // nguess = 0;

  for (j=0; j<spec->spec_len; j++) {
    fprintf(stderr,"%f %f %f\n",spec->spec[j].lambda_mean, min_l, max_l);
    if ((spec->spec[j].lambda_mean>=min_l) && (spec->spec[j].lambda_mean <= max_l)) {
    //     fprintf(stderr,"lambda mean between min and max\n\t%f : %f %f\n", spec->spec[j].lambda_mean, min_l, max_l);
    r1 = gsl_spline_eval(spline1, spec->spec[j].lambda_mean, acc1);
    r2 = gsl_spline_eval(spline2, spec->spec[j].lambda_mean, acc2);
   }
   else {
     r1 = 0;
     r2 = 0;
   }

    /* Need to divide by the pixel width of the
       wavelength calibration at the given lambda */

    if ((r2+r1)!=0) {

      // test
      // spec->spec[j].flux = spec->spec[j].count / tval;
      spec->spec[j].flux = spec->spec[j].count / ((r2+r1)/2.);
      spec->spec[j].flux = spec->spec[j].flux/spec->spec[j].dlambda;
      // fprintf(stdout, "# %f %f %f\n", spec->spec[j].lambda_mean, spec->spec[j].dlambda, spec->spec[j].weight);

      if (quant_cont && (int)spec->spec[j].contam != -1)
        {
          // test
          //      spec->spec[j].contam = spec->spec[j].contam / tval;
          spec->spec[j].contam = spec->spec[j].contam / ((r2+r1)/2.);
          spec->spec[j].contam = spec->spec[j].contam/spec->spec[j].dlambda;
        }
    } else {
      spec->spec[j].flux = GSL_NAN;
    }

    fr1 = spec->spec[j].error / spec->spec[j].count; /* frac. error */
    fr2 = fabs((r2-r1)/((r1+r2)/2.0)); /* frac. error in resp. */
    fr = sqrt(fr1*fr1+fr2*fr2);
    spec->spec[j].ferror = fabs(spec->spec[j].flux) * fr;
  }

  gsl_spline_free(spline1);
  gsl_interp_accel_free(acc1);
  gsl_spline_free(spline2);
  gsl_interp_accel_free(acc2);
  free(x1);
  free(x2);
  free(y1);
  free(y2);

}

/**
 * Function: apply_response_functionII
 * This function applies a thoughput curve to a given spectrum
 * by dividing the count attribute of the spectrum. The flux
 * and ferror attributes of the spectrum are populated.
 *
 * Parameters:
 * @param spec       - a pointer to an existing spectrum array
 * @param resp_func  - a pointer to a spectrum structure containing
 *                     a throughput curve previously loaded
 * @param quant_cont - the exposure time
 *
 * Returns
 * @return -
 */
void
apply_response_functionII(spectrum *spec, response_function *resp_func, const int quant_cont)
{
  long j;
  double fr, fr1, fr2;

  double *resp_vals = malloc( 2 * sizeof(double));

  // immediately return
  // an empty spectrum
  if ( spec==NULL)
    {
      aXe_message (aXe_M_WARN4, __FILE__, __LINE__,
                   "apply_response_function: spectra empty.");
      return;
    }

  // go over all elements
  for (j=0;j<spec->spec_len;j++)
    {
      // get the sensitivity value plus error
      get_response_values(resp_func, spec->spec[j].lambda_mean, resp_vals);

      // make sure you can compute the flux value
      if (resp_vals[0] != 0.0)
        {
          // compute the flux value
          spec->spec[j].flux = spec->spec[j].count / resp_vals[0];
          spec->spec[j].flux = spec->spec[j].flux/spec->spec[j].dlambda;

          // compute the flux error
          fr1 = spec->spec[j].error / spec->spec[j].count; /* frac. error */
          fr2 = fabs(resp_vals[1] / resp_vals[0]); /* frac. error in resp. */
          fr = sqrt(fr1*fr1+fr2*fr2);
          spec->spec[j].ferror = fabs(spec->spec[j].flux) * fr;

          // check whether the contamination value is necessary
          if (quant_cont && (int)spec->spec[j].contam != -1)
            {
              // compute the contamination value
              spec->spec[j].contam = spec->spec[j].contam / resp_vals[0];
              spec->spec[j].contam = spec->spec[j].contam/spec->spec[j].dlambda;
            }
        }
      else
        {
          // set the flux and the error to NaN
          spec->spec[j].flux = GSL_NAN;
          spec->spec[j].ferror = GSL_NAN;
        }
    }

  // free the memory
  free(resp_vals);
}

/**
 * Function: get_smooth_pars_for_beam
 * The function determines the smoothing parameters for an object beam
 * and a given indicator. If smoothed calibration is requested, the source
 * size (in dispersion direction) is compared to the poin-like object size.
 * For non poin-like objects, the adjustment factor is extracted and
 * the smoothing size from the object is derived.
 *
 * Parameters:
 * @param conf        - the configuarion file structure
 * @param smooth_conv - boolean for smooth conversion
 * @param actbeam     - the beam to be calibrated
 *
 * Returns
 * @return smooth_pars - the paramters used in the smoothing
 */
d_point
get_smooth_pars_for_beam(const aperture_conf *conf, const int smooth_conv, beam actbeam)
  {
    d_point smooth_pars;

    // initialize the return
    smooth_pars.x = -1.0;
    smooth_pars.y = -1.0;

    // if there is nothing to
    // do, return the default
    if (!smooth_conv)
      return smooth_pars;

    // check whether the relevant information
    // does exist
    if (smooth_conv && actbeam.slitgeom[2] < 0.0)
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
          "Smoothed flux conversion is impossible"
          " since the OAF\ndoes not contain the slit width!\n");

    //fprintf(stdout, "smoothing factor: %f, point-like size: %f\n", conf->smfactor, conf->pobjsize);
    // check whether the objects is smaller than a psf-object
    // in dispersion direction
    if (conf->pobjsize < actbeam.slitgeom[2])
      {
        // compute the smoothing size from the objects;
        // transfer the adjustment factor
        smooth_pars.x = pow((actbeam.slitgeom[2]* actbeam.slitgeom[2]) - (conf->pobjsize*conf->pobjsize), 0.5);
        smooth_pars.y = conf->smfactor;
      }

    // return the
    // smoothing parameters
    return smooth_pars;
  }

/**
 * Function: check_conf_for_smoothing
 * The functions checks whether a smoothed flux conversion
 * is possible or not. In case that keywords in the configuration
 * files are missing, it is NOT possible, and 0 is returned.
 *
 * Parameters:
 * @param conf        - the configuarion file structure
 * @param smooth_conv - boolean for smooth conversion
 *
 * Returns
 * @return is_possible - the paramters used in the smoothing
 */
int
check_conf_for_smoothing(const aperture_conf *conf, const int smooth_conv)
  {
    int is_possible=1;

    // make sure that, if smoothed covnersion is requested,
    // a point-like objects size and an adjustment factor is defined
    if (smooth_conv && (conf->pobjsize < 0.0 || conf->smfactor < 0.0))
      // change the switch
      is_possible = 0;

    // return the pointer
    return is_possible;
  }

/**
 * Function: apply_smoothed_response
 * Convert a spectrum from e/s to flux units in f_lambda. The conversion is done
 * using smoothed sensitivity values that take the object width into account
 * by smoothing the sensitivity function. Gaussian smoothing is applied, taking
 * into account that the original calibration function was derived for a point
 * like object with a finite size.
 *
 * Parameters:
 * @param spec           - a pointer to an existing spectrum array
 * @param resp_func      - a pointer to a spectrum structure containing
 *                         a throughput curve
 * @param quant_cont     - indicator for quantitative contamination
 * @param bwidth         - the B_IMAGE value of the beam
 * @param psf_ext        - the width of a point-like object
 * @param wl_calibration - the wavelength calibration structure
 *
 * Returns
 * @return -
 */
void
apply_smoothed_response(const calib_function *wl_calibration, const int for_grism,
                        const int quant_cont, response_function *resp_func,
                        const d_point smooth_pars, spectrum *spec)
{
  long j;
  double fr, fr1, fr2;
  //double sigma_wav=10.0;

  double *resp_vals = malloc( 2 * sizeof(double));
  //double *tmp_vals  = malloc( 2 * sizeof(double));

  gsl_vector *weights;

  // allocate memory for the vectors
  weights = gsl_vector_alloc(2 * RESP_SMOOTH_LENGTH + 1);

  // immediately return
  // if an empty spectrum
  if ( (spec==NULL) || (spec->spec_len == 0))
    {
      aXe_message (aXe_M_WARN4, __FILE__, __LINE__,
                   "apply_response_function: spectra empty.");
      return;
    }

  // fill the weight vector
  // for smoothing
  fill_weight(weights);

  // output to screen
  fprintf(stdout, "aXe_PET2SPC: Gaussian smoothing in tracelength: %.2f * %.2f = %.2f pix\n", smooth_pars.x, smooth_pars.y, smooth_pars.x * smooth_pars.y);

  // go over all elements
  for (j=0;j<spec->spec_len;j++)
    {
      // get the smoothed sensitivity value plus error
      get_smoothed_response(spec->spec[j].lambda_mean, smooth_pars, wl_calibration,
                            for_grism, weights, resp_func, resp_vals);
      // make sure you can compute the flux value
      if (resp_vals[0] != 0.0)
        {
          // compute the flux value
          spec->spec[j].flux = spec->spec[j].count / resp_vals[0];
          spec->spec[j].flux = spec->spec[j].flux/spec->spec[j].dlambda;

          // compute the flux error
          fr1 = spec->spec[j].error / spec->spec[j].count; /* frac. error */
          fr2 = fabs(resp_vals[1] / resp_vals[0]); /* frac. error in resp. */
          fr = sqrt(fr1*fr1+fr2*fr2);
          spec->spec[j].ferror = fabs(spec->spec[j].flux) * fr;

          // check whether the contamination value is necessary
          if (quant_cont && (int)spec->spec[j].contam != -1)
            {
              // compute the contamination value
              spec->spec[j].contam = spec->spec[j].contam / resp_vals[0];
              spec->spec[j].contam = spec->spec[j].contam/spec->spec[j].dlambda;
            }
        }
      else
        {
          // set the flux and the error to NaN
          spec->spec[j].flux = GSL_NAN;
          spec->spec[j].ferror = GSL_NAN;
        }
    }

  // free the memory
  free(resp_vals);
  gsl_vector_free(weights);
}

/**
 * Function: find_wavelength
 * For a given trace-length value, the function returns the difference
 * between the actual and the targeted wavelength.
 * The function is designed t be called by a gsl root-solving procedure.
 *
 * Parameters:
 * @param x        - guess value for the tracelength
 * @param params   - structure with input parameters
 *
 * Returns
 * @return wav_zero - difference between actual and targeted wavelength
 */
double
find_wavelength(double x, void *params)
  {
    double wav_zero;

    // extract the components from the parameter-structure
    const calib_function *wl_calib = ((trlength_search *) params)->wl_calibration;
    double wavelength        = ((trlength_search *) params)->wavelength;

    // compute the function value
    wav_zero = wl_calib->func (x, wl_calib->order, wl_calib->coeffs) - wavelength;

    // return the function value
    return wav_zero;
  }

/**
 * Function: get_tlength_prismwav
 * The function determines the tracelength value for a given wavelength
 * value and dispersion function. Should be called only for a prism dispersion.
 * The gsl root solver is used.
 *
 * Parameters:
 * @param wavelength     - the wavelength value
 * @param wl_calibration - the dispersion solution
  *
 * Returns
 * @return wav_zero - difference between actual and targeted wavelength
 */
double
get_tlength_prismwav(const double wavelength, const calib_function *wl_calibration)
{
  int iter=0;
  int status;

  d_point x_interv;

  double tr_length=0.0;

  // 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;
  trlength_search *tr_pars;

  // allocate and fill the parameters
  tr_pars = (trlength_search *) malloc(sizeof(trlength_search));
  tr_pars->wl_calibration = wl_calibration;
  tr_pars->wavelength     = wavelength;

  // fille the GSL-function
  F.function = &find_wavelength;
  F.params =   tr_pars;

  // initialize the intervall variable
  x_interv.x = gsl_vector_get(wl_calibration->pr_range, 0) + wl_calibration->coeffs[0];
  x_interv.y = gsl_vector_get(wl_calibration->pr_range, 1) + wl_calibration->coeffs[0];

  // 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
      tr_length = 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);
    }
  // check for the break condition
  while (status == GSL_CONTINUE && iter < MAX_ITER_TL);

  // free the memory
  free(tr_pars);

  // free the memory
  gsl_root_fsolver_free (s);

  // return the result
  return tr_length;
}

/**
 * Function: get_tlength_prismwav
 * The function determines the tracelength value for a given wavelength
 * value and dispersion function. Should be called only for a prism dispersion.
 * The gsl root solver is used.
 *
 * Parameters:
 * @param wavelength     - the wavelength value
 * @param wl_calibration - the dispersion solution
 * @param for_grism      - boolean to indicate grism/prism dispersion solution
 *
 * Returns
 * @return wav_zero - difference between actual and targeted wavelength
 */
double
get_central_tracelength(const double wavelength, const calib_function *wl_calibration,
                        const int for_grism)
  {
    double tl_central=0.0;
    double a, b, c;

    // checc whether it's a grism
    if (for_grism)
      {
        // solve a linear dispersion
        if (wl_calibration->order == 1)
          {
            tl_central = (wavelength - wl_calibration->coeffs[0]) / wl_calibration->coeffs[1];
          }
        // solve a quadratic dispersion
        else if (wl_calibration->order == 2)
          {
            // make the quantities
            a = wl_calibration->coeffs[2];
            b = wl_calibration->coeffs[1];
            c = wl_calibration->coeffs[0] - wavelength;

            // compute the central tracelength
            //tl_central = (-b + sqrt(b*b - 4.0*a*c)) / (2.0*c);
            tl_central = (-b + sqrt(b*b - 4.0*a*c)) / (2.0*a);
            //fprintf(stdout, "I should be here... %g %g %g %g\n", tl_central, a, b, c);
          }
        else
          {
            // give error message
            aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                "Smoothed flux conversion is impossible for grism "
                "data with dispersion solution of order: %i", wl_calibration->order);
          }
      }
    else
      {
        // use the root finder to get the tace-length value
        // for the given wavelength
        if (isnan(wavelength))
          tl_central = GSL_NAN;
        else
          tl_central = get_tlength_prismwav(wavelength, wl_calibration);
      }

    // return the central
    // tracelength value
    return tl_central;
  }

/**
 * Function: get_smoothed_response
 * The function determine a smoothe sensitivity value and its error.
 * Gaussian smoothing is applied, but the Gaussian weights are delivered
 * as a function parameter (for perfomrance reasons). The function detemines
 * the response values over the weight window and then combines these values
 * with the weights to get one weighted sensitivity value at the desired
 * wavelength.
 *
 * Parameters:
 * @param wavelength     - the central wavelength
 * @param smooth_pars    - the smoothing parameters
 * @param calib_function - the wavelength calibration function
 * @param for_grism      - boolean to indicate grism/prism solution
 * @param weights        - the weight vector
 * @param resp_func      - the response function
 * @param resp_vals      - vector for the response values
 *
 * Returns
 * @return -
 */
void
get_smoothed_response(const double wavelength, const d_point smooth_pars,
                           const calib_function *wl_calibration, const int for_grism,
                           const gsl_vector *weights, response_function *resp_func, double *resp_vals)
{
  int j;

  double tl_incr;
  double tl_act;
  double tl_central;
  double lambda_act;

  gsl_vector *pixvalues;
  gsl_vector *errvalues;
  gsl_vector *pmask;

  double *tmp_vals = malloc( 2 * sizeof(double));

  // allocate memory for the vectors
  pixvalues = gsl_vector_alloc(2 * RESP_SMOOTH_LENGTH + 1);
  errvalues = gsl_vector_alloc(2 * RESP_SMOOTH_LENGTH + 1);
  pmask     = gsl_vector_alloc(2 * RESP_SMOOTH_LENGTH + 1);

  // set all mask points to zero
  gsl_vector_set_all(pmask, 0.0);

  // determine the tracelength value
  tl_central = get_central_tracelength(wavelength, wl_calibration, for_grism);

  // determine the wavelength increments
  //wav_incr = (double)RESP_SMOOTH_NSIG * sigma_wav / (double)RESP_SMOOTH_LENGTH;
  tl_incr = smooth_pars.x * smooth_pars.y * (double)RESP_SMOOTH_NSIG / (double)RESP_SMOOTH_LENGTH;

  // set the initial wavelength value
  //value = wavelength - sigma_wav * (double)RESP_SMOOTH_NSIG;
  tl_act = tl_central - 1.0 * smooth_pars.x * smooth_pars.y * (double)RESP_SMOOTH_NSIG;

  // print the starting and central wavelength
  //fprintf(stdout, " %.5g < %.5g",
  //    wl_calibration->func(tl_act, wl_calibration->order, wl_calibration->coeffs),
  //    wl_calibration->func(tl_central, wl_calibration->order, wl_calibration->coeffs));

  // go over the array
  for (j=0; j < (int)pixvalues->size; j++)
    {
      // determine the wavelength value
      lambda_act = wl_calibration->func(tl_act, wl_calibration->order, wl_calibration->coeffs);

      // check whether the actual wavelength is inside
      // the defined sensitivity interval
      if  (lambda_act < resp_func->xmin || lambda_act > resp_func->xmax)
        {
          // for outside values:
          // set everything to zero
          gsl_vector_set(pixvalues, j, 0.0);
          gsl_vector_set(errvalues, j, 0.0);
          gsl_vector_set(pmask, j, 0.0);
        }
      else
        {
          // for inside values:
          // compute and set the response
          get_response_values(resp_func, lambda_act, tmp_vals);
          gsl_vector_set(pixvalues, j, tmp_vals[0]);
          gsl_vector_set(errvalues, j, tmp_vals[1]);
          gsl_vector_set(pmask, j, 1.0);
        }

      // increment trace length
      tl_act += tl_incr;
    }
  // print the ending wavelength
  //fprintf(stdout, " > %.5g;", wl_calibration->func(tl_act, wl_calibration->order, wl_calibration->coeffs));

  // get the weighted response values
  get_weighted_sensitivity(pixvalues, errvalues, weights, pmask, resp_vals);

  // free the memory
  free(tmp_vals);
  gsl_vector_free(pixvalues);
  gsl_vector_free(errvalues);
  gsl_vector_free(pmask);
}

/**
 * Function: get_weighted_sensitivity
 * Combines the weigh values, the sensitivity and error at the weight point
 * to a single sensitivity value plus error. A mask array defines at which
 * points values do exist.
 *
 * Parameters:
 * @param pixvalues - the central wavelength
 * @param errvalues - the Gaussian sigma used in smoothing
 * @param weights   - the weight vector
 * @param pmask     - the response function
 * @param resp_vals - vector for the response values
 *
 * Returns
 * @return
 */
void
get_weighted_sensitivity(const gsl_vector *pixvalues, const gsl_vector *errvalues, const gsl_vector *weights,
                         const gsl_vector *pmask, double *resp_vals)
{
  int index;

  // initialize the total
  // sum and weight
  double val_sum = 0.0;
  double err_sum = 0.0;
  double www=0.0;

  // go over all indices of an array
  for (index=0; index < (int)pixvalues->size; index ++)
    {
      // check whether the index is masked
      if (gsl_vector_get(pmask, index) != 0.0)
        {
          // enhance the total sum
          val_sum += gsl_vector_get(pixvalues, index) * gsl_vector_get(weights, index);
          err_sum += gsl_vector_get(errvalues, index) * gsl_vector_get(errvalues, index)
            * gsl_vector_get(weights, index) * gsl_vector_get(weights, index);

          // enhance the total weight
          www += gsl_vector_get(weights, index);
        }
    }

  // check whether there
  // were weights at all
  if (www)
    {
      // return the total sum,
      // divided by the total weight
      resp_vals[0] =  val_sum / www;
      resp_vals[1] =  sqrt(err_sum) / www;
    }
  else
    {
      // set zero sensitivity
      resp_vals[0] = 0.0;
      resp_vals[1] = 0.0;
    }
}

/**
 * Function: fill_weight
 * The function fills an array with Gaussian values.
 * The number of points and the extend are all defined in
 * global variables in the h-file.
 *
 * Parameters:
 * @param weights   - the weight vector
 *
 * Returns
 * @return -
 */
void
fill_weight(gsl_vector *weights)
{
  int j;

  double d_incr;
  double value;

  // determine the increments
  d_incr = (double)RESP_SMOOTH_NSIG / (double)RESP_SMOOTH_LENGTH;

  // compute the initial exponent
  // go over the weights vector
  value = -(double)RESP_SMOOTH_NSIG;
  for (j=0; j < (int)weights->size; j++)
    {
      // compute and set the weight
      gsl_vector_set(weights, j, exp(-1.0*value * value));

      // increment the exponent
      value += d_incr;
    }
}

/**
 * Function: get_troughput_table_name
 * Parses a configuration file and get the name of the SENSITIVITY
 * curve defined for a wanted beamID
 *
 * Parameters:
 * @param filename   - a pointer to a string containing the name of the
 *                     configuration file
 * @param beamID     - the beam ID number (see BEAM())
 * @param table_name - a pointer pointing to a string to contain the name
 *                     of the file.
 *
 * Returns:
 * @return -
 */
void
get_troughput_table_name(char *filename, int beamID, char *table_name)
{
  char beam[MAXCHAR];
  char name[MAXCHAR];
  int i=0;

  struct CfgStrings ThrConfig[] = {
    {NULL, NULL},
    {NULL,NULL}
  };
  ThrConfig[0].name = beam;
  sprintf (beam, "SENSITIVITY_%c", BEAM(beamID));

  CfgRead (filename, ThrConfig);
  if ((ThrConfig[0].name == beam) && (ThrConfig[0].data != NULL))
    {
      strcpy(name,ThrConfig[0].data);
    } else {
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "get_throughput_table_name: %s was not found in "
                   "%s",beam,filename);
    }

  // release memory
  i=0;
  while(ThrConfig[i].name!=NULL)
    free(ThrConfig[i++].data);

  strcpy(table_name,name);
}

/**
 * Function: get_response_value_plus
 * The function determines the sensitivity at a given wavelength
 * in a sensitivity curve. A guess value helps to start searching
 * around the estimated position of the wavelength in the vector
 * of the sensitifity curve. This should considerably speed up the
 * code. the guess value is updated in case that a sensitivity
 * value could be determined. The return value is interpolated
 * between the adjacent tabulated values.
 *
 * Parameters:
 * @param  resp       - the sensitivity spectrum
 * @param  wavelength - the wavelength to determine the sensitivity for
 * @param  nguess     - the guess value for the array to start searching from
 *
 * Returns:
 * @return ret        - the sensitivity value at the input wavelength
 */
double
get_response_value_plus(const spectrum *resp, const double wavelength,
                        int *nguess)
{
  double ret;
  double factor;

  int nact;

  // check whether the wavelength is within the range of
  //the sensitivity curve return GSL_NAN if not
  if (wavelength < resp->lambdamin || wavelength > resp->lambdamax)
    return GSL_NAN;

  // check whether you have to search upwards or downwards
  if (wavelength >= resp->spec[*nguess].lambda_mean)
    {

      // in case that you search upwards, go up
      // the spectrum until you find the right interval
      nact = *nguess + 1;
      while(wavelength > resp->spec[nact].lambda_mean)
        nact++;
    }
  else
    {

      // in case that you search upwards, go up
      // the spectrum until you find the right interval
      nact = *nguess;
      while(wavelength < resp->spec[nact-1].lambda_mean)
        nact--;
    }

  // interpolate within the interval to calculate the
  // sensitivity
  factor = (wavelength - resp->spec[nact-1].lambda_mean) /
    (resp->spec[nact].lambda_mean - resp->spec[nact-1].lambda_mean);
  ret = resp->spec[nact-1].flux + factor * (resp->spec[nact].flux - resp->spec[nact-1].flux);

  // update the guess value
  *nguess = nact;

  // return the sensitivity value
  return ret;
}

/**
 * Function: get_response_values
 * The function evaluates the response function at the desired
 * independent value. The sensitivity and it error are returned
 * in an array parameter.
 *
 * Parameters:
 * @param resp_func  - the response function
 * @param wavelength - the independent value
 * @param resp_vals  - the dependent value and its error
 *
 * Returns:
 * @return -
 */
void
get_response_values(response_function *resp_func, double wavelength, double* resp_vals)
{
  // check whether the independent value is outside the deifned range
  if (wavelength < resp_func->xmin || wavelength > resp_func->xmax)
    {
      // give the default
      // outside values
      resp_vals[0] = 0.0;
      resp_vals[1] = 0.0;
    }
  else
    {
      // get the response value and its error
      resp_vals[0] = eval_interp(resp_func->resp_values, wavelength);
      resp_vals[1] = eval_interp(resp_func->resp_errors, wavelength);
    }
}

/**
 * Function: create_response_function
 * Load the sensitivity file and creates a new response
 * function. Relies heavily on other subroutines
 *
 * Parameters:
 * @param filename - sensitivity file name
 *
 * Returns:
 * @return resp_func - the response function created
 */
response_function *
create_response_function(char *filename)
{
  response_function *resp_func;

  // allocate space for the return structure;
  // complain if this fails
  resp_func = (response_function *)malloc (sizeof (response_function));
  if (resp_func == NULL)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                 "Could not allocate memory for response function");

  // load the two intgerpolators
  resp_func->resp_values = create_interp_ftable(filename, 2, "WAVELENGTH", "SENSITIVITY", RESP_FUNC_INTERP_TYPE);
  resp_func->resp_errors = create_interp_ftable(filename, 2, "WAVELENGTH", "ERROR", RESP_FUNC_INTERP_TYPE);

  // transfer the min and max values from the interpolator
  resp_func->xmin = resp_func->resp_values->xmin;
  resp_func->xmax = resp_func->resp_values->xmax;

  // return the respons function
  return resp_func;
}

/**
 * Function: free_response_function
 * Releases the memory allocated in a response function
 *
 * Parameters:
 * @param resp_func - the response function
 *
 * Returns:
 * @return -
  */
void
free_response_function(response_function *resp_func)
{
  // free the two interpolators
  free_interp(resp_func->resp_values);
  free_interp(resp_func->resp_errors);

  // free the rest
  free(resp_func);

  // set it to NULL
  resp_func = NULL;

}

/**
 * Function: get_calfunc_for_beam
 * The subroutine determines and returns the calibration function
 * for a certain beam and configuration file. The coefficients of
 * the wavelength calibration are evaluated at the position of the
 * beam object.
 *
 * Parameters:
 * @param actbeam   - the current beam
 * @param for_grism - indicator for grism or prism dispersion
 * @param CONF_file - full pathname to configuration file
 * @param conf      - the configuration structure
 *
 * Returns:
 * @return wl_calibration - the calibration function
  */
calib_function *
get_calfunc_for_beam(const beam actbeam, const int for_grism, char CONF_file[],
                     const aperture_conf * conf)
{
  calib_function *wl_calibration;
  dispstruct     *disp;
  d_point         pixel;

  // get the reference point right
  pixel.x = actbeam.refpoint.x - conf->refx;
  pixel.y = actbeam.refpoint.y - conf->refy;

  // get the dispersion structure
  disp = get_dispstruct_at_pos(CONF_file, for_grism, actbeam.ID, pixel);

  // transform the dispersion structure into the
  // wavelength calibration
  wl_calibration = create_calib_from_gsl_vector(for_grism, disp->pol);
  if (!for_grism)
    wl_calibration->pr_range = get_prange (CONF_file, actbeam.ID);

  // free the whole structure
  free_dispstruct(disp);

  // return the wavelength calibration
  return wl_calibration;
}

ext: c

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/* linalg/lu.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2009 Gerard Jungman, Brian Gough
 * Copyright (C) 2019 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 <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

#include "recurse.h"

static int LU_decomp_L2 (gsl_matrix * A, gsl_vector_uint * ipiv);
static int LU_decomp_L3 (gsl_matrix * A, gsl_vector_uint * ipiv);
static int singular (const gsl_matrix * LU);
static int apply_pivots(gsl_matrix * A, const gsl_vector_uint * ipiv);

/* Factorise a general N x N matrix A into,
 *
 *   P A = L U
 *
 * where P is a permutation matrix, L is unit lower triangular and U
 * is upper triangular.
 *
 * L is stored in the strict lower triangular part of the input
 * matrix. The diagonal elements of L are unity and are not stored.
 *
 * U is stored in the diagonal and upper triangular part of the
 * input matrix.
 *
 * P is stored in the permutation p. Column j of P is column k of the
 * identity matrix, where k = permutation->data[j]
 *
 * signum gives the sign of the permutation, (-1)^n, where n is the
 * number of interchanges in the permutation.
 *
 * See Golub & Van Loan, Matrix Computations, Algorithm 3.4.1 (Gauss
 * Elimination with Partial Pivoting).
 */

int
gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * p, int *signum)
{
  const size_t M = A->size1;

  if (p->size != M)
    {
      GSL_ERROR ("permutation length must match matrix size1", GSL_EBADLEN);
    }
  else
    {
      int status;
      const size_t N = A->size2;
      const size_t minMN = GSL_MIN(M, N);
      gsl_vector_uint * ipiv = gsl_vector_uint_alloc(minMN);
      gsl_matrix_view AL = gsl_matrix_submatrix(A, 0, 0, M, minMN);
      size_t i;

      status = LU_decomp_L3 (&AL.matrix, ipiv);

      /* process remaining right matrix */
      if (M < N)
        {
          gsl_matrix_view AR = gsl_matrix_submatrix(A, 0, M, M, N - M);

          /* apply pivots to AR */
          apply_pivots(&AR.matrix, ipiv);

          /* AR = AL^{-1} AR */
          gsl_blas_dtrsm(CblasLeft, CblasLower, CblasNoTrans, CblasUnit, 1.0, &AL.matrix, &AR.matrix);
        }

      /* convert ipiv array to permutation */

      gsl_permutation_init(p);
      *signum = 1;

      for (i = 0; i < minMN; ++i)
        {
          unsigned int pivi = gsl_vector_uint_get(ipiv, i);

          if (p->data[pivi] != p->data[i])
            {
              size_t tmp = p->data[pivi];
              p->data[pivi] = p->data[i];
              p->data[i] = tmp;
              *signum = -(*signum);
            }
        }

      gsl_vector_uint_free(ipiv);

      return status;
    }
}

/*
LU_decomp_L2
  LU decomposition with partial pivoting using Level 2 BLAS

Inputs: A    - on input, matrix to be factored; on output, L and U factors
        ipiv - (output) array containing row swaps

Notes:
1) Based on LAPACK DGETF2
*/

static int
LU_decomp_L2 (gsl_matrix * A, gsl_vector_uint * ipiv)
{
  const size_t M = A->size1;
  const size_t N = A->size2;
  const size_t minMN = GSL_MIN(M, N);

  if (ipiv->size != minMN)
    {
      GSL_ERROR ("ipiv length must equal MIN(M,N)", GSL_EBADLEN);
    }
  else
    {
      size_t i, j;

      for (j = 0; j < minMN; ++j)
        {
          /* find maximum in the j-th column */
          gsl_vector_view v = gsl_matrix_subcolumn(A, j, j, M - j);
          size_t j_pivot = j + gsl_blas_idamax(&v.vector);
          gsl_vector_view v1, v2;

          gsl_vector_uint_set(ipiv, j, j_pivot);

          if (j_pivot != j)
            {
              /* swap rows j and j_pivot */
              v1 = gsl_matrix_row(A, j);
              v2 = gsl_matrix_row(A, j_pivot);
              gsl_blas_dswap(&v1.vector, &v2.vector);
            }

          if (j < M - 1)
            {
              double Ajj = gsl_matrix_get(A, j, j);

              if (fabs(Ajj) >= GSL_DBL_MIN)
                {
                  v1 = gsl_matrix_subcolumn(A, j, j + 1, M - j - 1);
                  gsl_blas_dscal(1.0 / Ajj, &v1.vector);
                }
              else
                {
                  for (i = 1; i < M - j; ++i)
                    {
                      double * ptr = gsl_matrix_ptr(A, j + i, j);
                      *ptr /= Ajj;
                    }
                }
            }

          if (j < minMN - 1)
            {
              gsl_matrix_view A22 = gsl_matrix_submatrix(A, j + 1, j + 1, M - j - 1, N - j - 1);
              v1 = gsl_matrix_subcolumn(A, j, j + 1, M - j - 1);
              v2 = gsl_matrix_subrow(A, j, j + 1, N - j - 1);

              gsl_blas_dger(-1.0, &v1.vector, &v2.vector, &A22.matrix);
            }
        }

      return GSL_SUCCESS;
    }
}

/*
LU_decomp_L3
  LU decomposition with partial pivoting using Level 3 BLAS

Inputs: A    - on input, matrix to be factored; on output, L and U factors
        ipiv - (output) array containing row swaps

Notes:
1) Based on ReLAPACK DGETRF
*/

static int
LU_decomp_L3 (gsl_matrix * A, gsl_vector_uint * ipiv)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M < N)
    {
      GSL_ERROR ("matrix must have M >= N", GSL_EBADLEN);
    }
  else if (ipiv->size != GSL_MIN(M, N))
    {
      GSL_ERROR ("ipiv length must equal MIN(M,N)", GSL_EBADLEN);
    }
  else if (N <= CROSSOVER_LU)
    {
      /* use Level 2 algorithm */
      return LU_decomp_L2(A, ipiv);
    }
  else
    {
      /*
       * partition matrix:
       *
       *       N1  N2
       * N1  [ A11 A12 ]
       * M2  [ A21 A22 ]
       *
       * and
       *      N1  N2
       * M  [ AL  AR  ]
       */
      int status;
      const size_t N1 = GSL_LINALG_SPLIT(N);
      const size_t N2 = N - N1;
      const size_t M2 = M - N1;
      gsl_matrix_view A11 = gsl_matrix_submatrix(A, 0, 0, N1, N1);
      gsl_matrix_view A12 = gsl_matrix_submatrix(A, 0, N1, N1, N2);
      gsl_matrix_view A21 = gsl_matrix_submatrix(A, N1, 0, M2, N1);
      gsl_matrix_view A22 = gsl_matrix_submatrix(A, N1, N1, M2, N2);

      gsl_matrix_view AL = gsl_matrix_submatrix(A, 0, 0, M, N1);
      gsl_matrix_view AR = gsl_matrix_submatrix(A, 0, N1, M, N2);

      /*
       * partition ipiv = [ ipiv1 ] N1
       *                  [ ipiv2 ] N2
       */
      gsl_vector_uint_view ipiv1 = gsl_vector_uint_subvector(ipiv, 0, N1);
      gsl_vector_uint_view ipiv2 = gsl_vector_uint_subvector(ipiv, N1, N2);

      size_t i;

      /* recursion on (AL, ipiv1) */
      status = LU_decomp_L3(&AL.matrix, &ipiv1.vector);
      if (status)
        return status;

      /* apply ipiv1 to AR */
      apply_pivots(&AR.matrix, &ipiv1.vector);

      /* A12 = A11^{-1} A12 */
      gsl_blas_dtrsm(CblasLeft, CblasLower, CblasNoTrans, CblasUnit, 1.0, &A11.matrix, &A12.matrix);

      /* A22 = A22 - A21 * A12 */
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0, &A21.matrix, &A12.matrix, 1.0, &A22.matrix);

      /* recursion on (A22, ipiv2) */
      status = LU_decomp_L3(&A22.matrix, &ipiv2.vector);
      if (status)
        return status;

      /* apply pivots to A21 */
      apply_pivots(&A21.matrix, &ipiv2.vector);

      /* shift pivots */
      for (i = 0; i < N2; ++i)
        {
          unsigned int * ptr = gsl_vector_uint_ptr(&ipiv2.vector, i);
          *ptr += N1;
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_LU_solve (const gsl_matrix * LU, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else if (singular (LU))
    {
      GSL_ERROR ("matrix is singular", GSL_EDOM);
    }
  else
    {
      int status;

      /* copy x <- b */
      gsl_vector_memcpy (x, b);

      /* solve for x */
      status = gsl_linalg_LU_svx (LU, p, x);

      return status;
    }
}

int
gsl_linalg_LU_svx (const gsl_matrix * LU, const gsl_permutation * p, gsl_vector * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution/rhs size", GSL_EBADLEN);
    }
  else if (singular (LU))
    {
      GSL_ERROR ("matrix is singular", GSL_EDOM);
    }
  else
    {
      /* apply permutation to RHS */
      gsl_permute_vector (p, x);

      /* solve for c using forward-substitution, L c = P b */
      gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasUnit, LU, x);

      /* perform back-substitution, U x = c */
      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LU, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_LU_refine (const gsl_matrix * A, const gsl_matrix * LU, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x, gsl_vector * work)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("matrix a must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (A->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be decomposition of a", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else if (LU->size1 != work->size)
    {
      GSL_ERROR ("matrix size must match workspace size", GSL_EBADLEN);
    }
  else if (singular (LU))
    {
      GSL_ERROR ("matrix is singular", GSL_EDOM);
    }
  else
    {
      int status;

      /* compute residual = (A * x  - b) */
      gsl_vector_memcpy (work, b);
      gsl_blas_dgemv (CblasNoTrans, 1.0, A, x, -1.0, work);

      /* find correction, delta = - (A^-1) * residual, and apply it */
      status = gsl_linalg_LU_svx (LU, p, work);
      gsl_blas_daxpy (-1.0, work, x);

      return status;
    }
}

int
gsl_linalg_LU_invert (const gsl_matrix * LU, const gsl_permutation * p, gsl_matrix * inverse)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (inverse->size1 != LU->size1 || inverse->size2 != LU->size2)
    {
      GSL_ERROR ("inverse matrix must match LU matrix dimensions", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_memcpy(inverse, LU);
      return gsl_linalg_LU_invx (inverse, p);
    }
}

int
gsl_linalg_LU_invx (gsl_matrix * LU, const gsl_permutation * p)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (singular (LU))
    {
      GSL_ERROR ("matrix is singular", GSL_EDOM);
    }
  else
    {
      int status;
      const size_t N = LU->size1;
      size_t i;

      /* compute U^{-1} */
      status = gsl_linalg_tri_invert(CblasUpper, CblasNonUnit, LU);
      if (status)
        return status;

      /* compute L^{-1} */
      status = gsl_linalg_tri_invert(CblasLower, CblasUnit, LU);
      if (status)
        return status;

      /* compute U^{-1} L^{-1} */
      status = gsl_linalg_tri_UL(LU);
      if (status)
        return status;

      /* apply permutation to columns of A^{-1} */
      for (i = 0; i < N; ++i)
        {
          gsl_vector_view v = gsl_matrix_row(LU, i);
          gsl_permute_vector_inverse(p, &v.vector);
        }

      return GSL_SUCCESS;
    }
}

double
gsl_linalg_LU_det (gsl_matrix * LU, int signum)
{
  size_t i, n = LU->size1;

  double det = (double) signum;

  for (i = 0; i < n; i++)
    {
      det *= gsl_matrix_get (LU, i, i);
    }

  return det;
}

double
gsl_linalg_LU_lndet (gsl_matrix * LU)
{
  size_t i, n = LU->size1;

  double lndet = 0.0;

  for (i = 0; i < n; i++)
    {
      lndet += log (fabs (gsl_matrix_get (LU, i, i)));
    }

  return lndet;
}

int
gsl_linalg_LU_sgndet (gsl_matrix * LU, int signum)
{
  size_t i, n = LU->size1;

  int s = signum;

  for (i = 0; i < n; i++)
    {
      double u = gsl_matrix_get (LU, i, i);

      if (u < 0)
        {
          s *= -1;
        }
      else if (u == 0)
        {
          s = 0;
          break;
        }
    }

  return s;
}

static int
singular (const gsl_matrix * LU)
{
  size_t i, n = LU->size1;

  for (i = 0; i < n; i++)
    {
      double u = gsl_matrix_get (LU, i, i);
      if (u == 0) return 1;
    }

  return 0;
}

static int
apply_pivots(gsl_matrix * A, const gsl_vector_uint * ipiv)
{
  if (A->size1 < ipiv->size)
    {
      GSL_ERROR("matrix does not match pivot vector", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      for (i = 0; i < ipiv->size; ++i)
        {
          size_t pi = gsl_vector_uint_get(ipiv, i);

          if (i != pi)
            {
              /* swap rows i and pi */
              gsl_vector_view v1 = gsl_matrix_row(A, i);
              gsl_vector_view v2 = gsl_matrix_row(A, pi);
              gsl_blas_dswap(&v1.vector, &v2.vector);
            }
        }

      return GSL_SUCCESS;
    }
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include "scf.h"

/* Computing Phi_nl(r), which is the radial part of the basis functions
 *
 * The Gegenbauer polynomials for the host halo and subhalos can either be computed
 * using my own hard coded implementation (see gegenbauer.c) or GSL's implementation.
 *
 */

double int_power(float x, int n)
{
    return gsl_pow_int(x, n);
}

#ifdef USE_GSL_GEGENBAUER_HOST
float Phi_nl(double Cna, int l, float r)
{
    return -1.0f * int_power(r,l)/int_power(1.0f+r, 2.0*l+1) * (float)Cna * ROOT_FOUR_PI;
}

float dPhi_nl_dr(double Cna, float phinl, int n, int l, float r)
{
    int two_l = 2*l;
	float r1 = r+1.f;

	float val = phinl * ( l/r - (two_l+1.0f)/r1 );

	if(n>0)
	{
		val -= ROOT_FOUR_PI * int_power(r,l)*4.0f*(two_l+1.5f) * (float)Cna / int_power(r1,two_l+3);
	}

	return val;
}

#else
float Phi_nl(int n, int l, float r)
{
    float xi = (r-1.0f)/(r+1.0f);

    return -1.0f * int_power(r,l)/int_power(1.0f+r, 2.0*l+1) * Cna(n,2.0f*l+1.5f,xi) * ROOT_FOUR_PI;
}

float dPhi_nl_dr(float phinl, int n, int l, float r)
{
    int two_l = 2*l;
	float r1 = r+1.f;

	float val = phinl * ( l/r - (two_l+1.0f)/r1 );

	if(n>0)
	{
		float xi = (r-1.f)/r1;

		val -= ROOT_FOUR_PI * int_power(r,l)*4.0f*(two_l+1.5f) * Cna(n-1,two_l+2.5f,xi) / int_power(r1,two_l+3);
	}

	return val;
}
#endif

#ifdef USE_GSL_GEGENBAUER_SUBHALO
float Phi_n0(double Cna, float r)
{
    return -1.0f / (1.0f+r) * (float)Cna * ROOT_FOUR_PI;
}

float dPhi_n0_dr(double Cna, float phin0, int n, float r)
{
    float r1 = r+1.f;

    float val = phin0 * -1.f/r1;

    if(n>0)
    {
        val -= ROOT_FOUR_PI * 6.f/r1/r1/r1 * (float)Cna;
    }

    return val;
}

#else
float Phi_n0(int n, float r)
{
    float xi = (r-1.0f)/(r+1.0f);

    return -1.0f / (1.0f+r) * Cna(n,1.5f,xi) * ROOT_FOUR_PI;
}

float dPhi_n0_dr(float phin0, int n, float r)
{
    float r1 = r+1.f;

    float val = phin0 * -1.f/r1;

    if(n>0)
    {
		float xi = (r-1.f)/r1;
        val -= ROOT_FOUR_PI * 6.f/r1/r1/r1 * Cna(n-1,2.5f,xi);
    }

    return val;
}

#endif


text:
/*
 * C version of Diffusive Nested Sampling (DNest4) by Brendon J. Brewer
 *
 * Yan-Rong Li, liyanrong@mail.ihep.ac.cn
 * Jun 30, 2016
 *
 */

#ifndef _MODEL2_H
#define _MODEL2_H

#include <stdbool.h>
#include <gsl/gsl_rng.h>

/*  data type */
typedef struct
{
  double x;
  double y;
}DataType;

/* number of model parameters */
extern int num_params;

/* data storage */
extern int num_data_points;
extern DataType *data;

extern int which_particle_update; // which particule to be updated
extern int which_level_update;
extern double *limits;
extern int thistask, totaltask;

extern DNestFptrSet *fptrset_thismodel2;

/* functions */
void from_prior_thismodel2(void *model);
void data_load_thismodel2();
void print_particle_thismodel2(FILE *fp, const void *model);
double log_likelihoods_cal_thismodel2(const void *model);
double perturb_thismodel2(void *model);
void restart_action_model2(int iflag);

#endif

ext: h

max_stars_repo_name: LiyrAstroph/CDNest


text:
/* linalg/test_cholesky.c
 *
 * Copyright (C) 2016 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_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>

#include "test_common.c"

static int test_cholesky_decomp_eps(const int scale, const gsl_matrix * m,
                                    const double expected_rcond, const double eps,
                                    const char * desc);
static int test_cholesky_decomp(gsl_rng * r);
int test_cholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc);
int test_cholesky_invert(gsl_rng * r);
static int test_pcholesky_decomp_eps(const int scale, const gsl_matrix * m,
                                     const double expected_rcond, const double eps,
                                     const char * desc);
static int test_pcholesky_decomp(gsl_rng * r);
int test_pcholesky_solve_eps(const int scale, const gsl_matrix * m, const gsl_vector * rhs,
                             const gsl_vector * sol, const double eps,
                             const char * desc);
static int test_pcholesky_solve(gsl_rng * r);
int test_pcholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc);
int test_pcholesky_invert(gsl_rng * r);

static int test_mcholesky_decomp_eps(const int posdef, const int scale, const gsl_matrix * m,
                                     const double expected_rcond, const double eps, const char * desc);

/* Hilbert matrix condition numbers, as calculated by LAPACK DPOSVX */
double hilb_rcond[] = { 1.000000000000e+00, 3.703703703704e-02, 1.336898395722e-03,
                        3.524229074890e-05, 1.059708198754e-06, 3.439939465186e-08,
                        1.015027593823e-09, 2.952221630602e-11, 9.093751565191e-13,
                        2.828277420229e-14, 8.110242564869e-16, 2.409320075800e-17 };

static int
test_cholesky_decomp_eps(const int scale, const gsl_matrix * m,
                         const double expected_rcond, const double eps,
                         const char * desc)
{
  int s = 0;
  size_t i, j, N = m->size2;

  gsl_matrix * V  = gsl_matrix_alloc(N, N);
  gsl_matrix * A  = gsl_matrix_alloc(N, N);
  gsl_matrix * L  = gsl_matrix_calloc(N, N);
  gsl_matrix * LT = gsl_matrix_calloc(N, N);
  gsl_vector * S = gsl_vector_alloc(N);

  gsl_matrix_memcpy(V, m);

  if (scale)
    s += gsl_linalg_cholesky_decomp2(V, S);
  else
    s += gsl_linalg_cholesky_decomp1(V);

  /* compute L and LT */
  gsl_matrix_tricpy('L', 1, L, V);
  gsl_matrix_transpose_tricpy('L', 1, LT, L);

  if (scale)
    {
      /* L <- S^{-1} L, LT <- LT S^{-1} */
      for (i = 0; i < N; ++i)
        {
          double Si = gsl_vector_get(S, i);
          gsl_vector_view v = gsl_matrix_row(L, i);
          gsl_vector_view w = gsl_matrix_column(LT, i);

          gsl_vector_scale(&v.vector, 1.0 / Si);
          gsl_vector_scale(&w.vector, 1.0 / Si);
        }
    }

  /* compute A = L LT */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, L, LT, 0.0, A);

  for (i = 0; i < N; i++)
    {
      for (j = 0; j < N; j++)
        {
          double Aij = gsl_matrix_get(A, i, j);
          double mij = gsl_matrix_get(m, i, j);

          gsl_test_rel(Aij, mij, eps,
                       "%s: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, Aij, mij);
        }
    }

  if (expected_rcond > 0 && !scale)
    {
      gsl_vector *work = gsl_vector_alloc(3 * N);
      double rcond;

      gsl_linalg_cholesky_rcond(V, &rcond, work);

      gsl_test_rel(rcond, expected_rcond, 1.0e-6,
                   "%s rcond: (%3lu,%3lu): %22.18g   %22.18g\n",
                   desc, N, N, rcond, expected_rcond);

      gsl_vector_free(work);
    }

  gsl_matrix_free(V);
  gsl_matrix_free(A);
  gsl_matrix_free(L);
  gsl_matrix_free(LT);
  gsl_vector_free(S);

  return s;
}

static int
test_cholesky_decomp(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 50;
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);
      test_cholesky_decomp_eps(0, m, -1.0, 1.0e2 * N * GSL_DBL_EPSILON, "cholesky_decomp unscaled random");
      test_cholesky_decomp_eps(1, m, -1.0, 1.0e2 * N * GSL_DBL_EPSILON, "cholesky_decomp scaled random");

      if (N <= 12)
        {
          double expected_rcond = -1.0;

          if (hilb_rcond[N - 1] > 1.0e-12)
            expected_rcond = hilb_rcond[N - 1];

          create_hilbert_matrix2(m);

          test_cholesky_decomp_eps(0, m, expected_rcond, N * GSL_DBL_EPSILON, "cholesky_decomp unscaled hilbert");
          test_cholesky_decomp_eps(1, m, expected_rcond, N * GSL_DBL_EPSILON, "cholesky_decomp scaled hilbert");
        }

      gsl_matrix_free(m);
    }

  return s;
}

int
test_cholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc)
{
  int s = 0;
  size_t i, j, N = m->size1;

  gsl_matrix * v  = gsl_matrix_alloc(N, N);
  gsl_matrix * c  = gsl_matrix_alloc(N, N);

  gsl_matrix_memcpy(v, m);

  s += gsl_linalg_cholesky_decomp1(v);
  s += gsl_linalg_cholesky_invert(v);

  /* c = m m^{-1} */
  gsl_blas_dsymm(CblasLeft, CblasUpper, 1.0, m, v, 0.0, c);

  /* c should be the identity matrix */

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          double cij = gsl_matrix_get(c, i, j);
          double expected = (i == j) ? 1.0 : 0.0;

          gsl_test_rel(cij, expected, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, cij, expected);
        }
    }

  gsl_matrix_free(v);
  gsl_matrix_free(c);

  return s;
}

int
test_cholesky_invert(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 50;
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);

      test_cholesky_invert_eps(m, N * GSL_DBL_EPSILON, "cholesky_invert unscaled random");

      if (N <= 4)
        {
          create_hilbert_matrix2(m);
          test_cholesky_invert_eps(m, N * 256.0 * GSL_DBL_EPSILON, "cholesky_invert unscaled hilbert");
        }

      gsl_matrix_free(m);
    }

  return s;
}

static int
test_mcholesky_decomp_eps(const int posdef, const int scale, const gsl_matrix * m,
                          const double expected_rcond, const double eps, const char * desc)
{
  int s = 0;
  size_t i, j, N = m->size2;

  gsl_matrix * LDLT = gsl_matrix_alloc(N, N);
  gsl_matrix * V  = gsl_matrix_alloc(N, N);
  gsl_matrix * A  = gsl_matrix_alloc(N, N);
  gsl_matrix * L  = gsl_matrix_alloc(N, N);
  gsl_matrix * LT = gsl_matrix_alloc(N, N);
  gsl_vector * S = gsl_vector_alloc(N);
  gsl_vector * E = gsl_vector_alloc(N);
  gsl_permutation * perm = gsl_permutation_alloc(N);

  gsl_vector_view D = gsl_matrix_diagonal(LDLT);

  gsl_matrix_memcpy(LDLT, m);
  s += gsl_linalg_mcholesky_decomp(LDLT, perm, E);

  /* check that the upper triangle of LDLT equals original matrix */
  for (i = 0; i < N; ++i)
    {
      for (j = i + 1; j < N; ++j)
        {
          double mij = gsl_matrix_get(m, i, j);
          double aij = gsl_matrix_get(LDLT, i, j);

          gsl_test_rel(aij, mij, 1.0e-12,
                       "%s upper triangle: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, aij, mij);
        }
    }

  if (posdef)
    {
      /* ||E|| should be 0 */
      double norm = gsl_blas_dnrm2(E);
      s = norm != 0.0;
      gsl_test(s, "%s: (%zu,%zu): ||E|| = %.12e",
               desc, N, N, norm);

      /* check that D is decreasing */
      s = 0;
      for (i = 1; i < N; ++i)
        {
          double dprev = gsl_vector_get(&D.vector, i - 1);
          double di = gsl_vector_get(&D.vector, i);

          if (di > dprev)
            s = 1;
        }

      gsl_test(s, "%s: (%zu,%zu): D is not decreasing",
               desc, N, N);
    }

  /* compute L and LT */
  gsl_matrix_set_identity(L);
  gsl_matrix_set_identity(LT);

  gsl_matrix_tricpy('L', 0, L, LDLT);
  gsl_matrix_transpose_tricpy('L', 0, LT, L);

  /* compute (L sqrt(D)) and (sqrt(D) LT) */
  for (i = 0; i < N; ++i)
    {
      gsl_vector_view v = gsl_matrix_column(L, i);
      gsl_vector_view w = gsl_matrix_row(LT, i);
      double di = gsl_vector_get(&D.vector, i);

      gsl_vector_scale(&v.vector, sqrt(di));
      gsl_vector_scale(&w.vector, sqrt(di));
    }

  /* compute A = L D LT */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, L, LT, 0.0, A);

  /* compute V = P (S M S + E) P^T */

  gsl_matrix_memcpy(V, m);
  D = gsl_matrix_diagonal(V);

  /* compute S M S */
  if (scale)
    {
      gsl_linalg_cholesky_scale_apply(V, S);
      gsl_matrix_transpose_tricpy('L', 0, V, V);
    }

  /* compute S M S + E */
  gsl_vector_add(&D.vector, E);

  /* compute M P^T */
  for (i = 0; i < N; ++i)
    {
      gsl_vector_view v = gsl_matrix_row(V, i);
      gsl_permute_vector(perm, &v.vector);
    }

  /* compute P M P^T */
  for (i = 0; i < N; ++i)
    {
      gsl_vector_view v = gsl_matrix_column(V, i);
      gsl_permute_vector(perm, &v.vector);
    }

  for (i = 0; i < N; i++)
    {
      double Ei = gsl_vector_get(E, i);

      for (j = 0; j < N; j++)
        {
          double Aij = gsl_matrix_get(A, i, j); /* L D L^T */
          double Bij = gsl_matrix_get(V, i, j); /* P M P^T */
          double Cij;                           /* P M P^T + E */

          if (i == j)
            Cij = Bij + Ei*0;
          else
            Cij = Bij;

          gsl_test_rel(Aij, Cij, eps,
                       "%s: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, Aij, Cij);
        }
    }

  if (expected_rcond > 0 && !scale)
    {
      gsl_vector *work = gsl_vector_alloc(3 * N);
      double rcond;

      gsl_linalg_mcholesky_rcond(LDLT, perm, &rcond, work);

      gsl_test_rel(rcond, expected_rcond, 1.0e-6,
                   "%s rcond: (%3lu,%3lu): %22.18g   %22.18g\n",
                   desc, N, N, rcond, expected_rcond);

      gsl_vector_free(work);
    }

  gsl_matrix_free(LDLT);
  gsl_matrix_free(V);
  gsl_matrix_free(A);
  gsl_matrix_free(L);
  gsl_matrix_free(LT);
  gsl_vector_free(S);
  gsl_vector_free(E);
  gsl_permutation_free(perm);

  return s;
}

static int
test_mcholesky_decomp(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 50;
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);
      test_mcholesky_decomp_eps(1, 0, m, -1.0, 128.0 * N * GSL_DBL_EPSILON, "mcholesky_decomp unscaled random posdef");

      create_symm_matrix(m, r);
      test_mcholesky_decomp_eps(0, 0, m, -1.0, 8192.0 * N * GSL_DBL_EPSILON, "mcholesky_decomp unscaled random symm");

      if (N <= 8)
        {
          double expected_rcond = -1.0;

          if (hilb_rcond[N - 1] > 1.0e-12)
            expected_rcond = hilb_rcond[N - 1];

          create_hilbert_matrix2(m);

          test_mcholesky_decomp_eps(1, 0, m, expected_rcond, 128.0 * N * GSL_DBL_EPSILON, "mcholesky_decomp unscaled hilbert");
        }

      gsl_matrix_free(m);
    }

  return s;
}

int
test_mcholesky_solve_eps(const gsl_matrix * m, const gsl_vector * rhs,
                         const gsl_vector * sol, const double eps,
                         const char * desc)
{
  int s = 0;
  size_t i, N = m->size1;
  gsl_matrix * u  = gsl_matrix_alloc(N, N);
  gsl_vector * x = gsl_vector_calloc(N);
  gsl_vector * S = gsl_vector_alloc(N);
  gsl_permutation * perm = gsl_permutation_alloc(N);

  gsl_matrix_memcpy(u, m);

  s += gsl_linalg_mcholesky_decomp(u, perm, NULL);
  s += gsl_linalg_mcholesky_solve(u, perm, rhs, x);

  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(sol, i);

      gsl_test_rel(xi, yi, eps,
                   "%s: %3lu[%lu]: %22.18g   %22.18g\n",
                   desc, N, i, xi, yi);
    }

  gsl_vector_free(x);
  gsl_vector_free(S);
  gsl_matrix_free(u);
  gsl_permutation_free(perm);

  return s;
}

static int
test_mcholesky_solve(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 50;
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix * m = gsl_matrix_alloc(N, N);
      gsl_vector * rhs = gsl_vector_alloc(N);
      gsl_vector * sol = gsl_vector_alloc(N);

      create_posdef_matrix(m, r);
      create_random_vector(sol, r);
      gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);

      test_mcholesky_solve_eps(m, rhs, sol, 64.0 * N * GSL_DBL_EPSILON, "mcholesky_solve random");

      if (N <= 3)
        {
          create_hilbert_matrix2(m);
          gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);
          test_mcholesky_solve_eps(m, rhs, sol, 1.0e3 * N * GSL_DBL_EPSILON, "mcholesky_solve hilbert");
        }

      gsl_matrix_free(m);
      gsl_vector_free(rhs);
      gsl_vector_free(sol);
    }

  return s;
}

int
test_mcholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc)
{
  int s = 0;
  size_t i, j, N = m->size1;

  gsl_matrix * v  = gsl_matrix_alloc(N, N);
  gsl_matrix * c  = gsl_matrix_alloc(N, N);
  gsl_matrix * minv = gsl_matrix_alloc(N, N);
  gsl_vector * E = gsl_vector_alloc(N);
  gsl_permutation * p = gsl_permutation_alloc(N);

  gsl_matrix_memcpy(v, m);

  s += gsl_linalg_mcholesky_decomp(v, p, E);
  s += gsl_linalg_mcholesky_invert(v, p, minv);

  /* c = m m^{-1} */
  gsl_blas_dsymm(CblasLeft, CblasUpper, 1.0, m, minv, 0.0, c);

  /* c should be the identity matrix */

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          double cij = gsl_matrix_get(c, i, j);
          double expected = (i == j) ? 1.0 : 0.0;

          gsl_test_rel(cij, expected, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, cij, expected);
        }
    }

  gsl_matrix_free(v);
  gsl_matrix_free(c);
  gsl_matrix_free(minv);
  gsl_vector_free(E);
  gsl_permutation_free(p);

  return s;
}

int
test_mcholesky_invert(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 30;
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);
      test_mcholesky_invert_eps(m, N * GSL_DBL_EPSILON, "mcholesky_invert unscaled random");

      if (N <= 4)
        {
          create_hilbert_matrix2(m);
          test_mcholesky_invert_eps(m, 256.0 * N * GSL_DBL_EPSILON, "mcholesky_invert unscaled hilbert");
        }

      gsl_matrix_free(m);
    }

  return s;
}

static int
test_pcholesky_decomp_eps(const int scale, const gsl_matrix * m,
                          const double expected_rcond, const double eps,
                          const char * desc)
{
  int s = 0;
  size_t i, j, N = m->size2;

  gsl_matrix * LDLT = gsl_matrix_alloc(N, N);
  gsl_matrix * V  = gsl_matrix_alloc(N, N);
  gsl_matrix * A  = gsl_matrix_alloc(N, N);
  gsl_matrix * L  = gsl_matrix_alloc(N, N);
  gsl_matrix * LT = gsl_matrix_alloc(N, N);
  gsl_vector * S = gsl_vector_alloc(N);
  gsl_permutation * perm = gsl_permutation_alloc(N);

  gsl_vector_view D = gsl_matrix_diagonal(LDLT);

  gsl_matrix_memcpy(LDLT, m);

  if (scale)
    s += gsl_linalg_pcholesky_decomp2(LDLT, perm, S);
  else
    s += gsl_linalg_pcholesky_decomp(LDLT, perm);

  /* check that the upper triangle of LDLT equals original matrix */
  for (i = 0; i < N; ++i)
    {
      for (j = i + 1; j < N; ++j)
        {
          double mij = gsl_matrix_get(m, i, j);
          double aij = gsl_matrix_get(LDLT, i, j);

          gsl_test_rel(aij, mij, 1.0e-12,
                       "%s upper triangle: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, aij, mij);
        }
    }

  /* check that D is decreasing */
  s = 0;
  for (i = 1; i < N; ++i)
    {
      double dprev = gsl_vector_get(&D.vector, i - 1);
      double di = gsl_vector_get(&D.vector, i);

      if (di > dprev)
        s = 1;
    }

  gsl_test(s, "%s: (%zu,%zu): D is not decreasing",
           desc, N, N);

  /* compute L and LT */
  gsl_matrix_set_identity(L);
  gsl_matrix_set_identity(LT);

  gsl_matrix_tricpy('L', 0, L, LDLT);
  gsl_matrix_transpose_tricpy('L', 0, LT, L);

  /* compute (L sqrt(D)) and (sqrt(D) LT) */
  for (i = 0; i < N; ++i)
    {
      gsl_vector_view v = gsl_matrix_column(L, i);
      gsl_vector_view w = gsl_matrix_row(LT, i);
      double di = gsl_vector_get(&D.vector, i);

      gsl_vector_scale(&v.vector, sqrt(di));
      gsl_vector_scale(&w.vector, sqrt(di));
    }

  /* compute A = L D LT */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, L, LT, 0.0, A);

  /* compute V = P S M S P^T */

  gsl_matrix_memcpy(V, m);

  /* compute S M S */
  if (scale)
    {
      gsl_linalg_cholesky_scale_apply(V, S);
      gsl_matrix_transpose_tricpy('L', 0, V, V);
    }

  /* compute M P^T */
  for (i = 0; i < N; ++i)
    {
      gsl_vector_view v = gsl_matrix_row(V, i);
      gsl_permute_vector(perm, &v.vector);
    }

  /* compute P M P^T */
  for (i = 0; i < N; ++i)
    {
      gsl_vector_view v = gsl_matrix_column(V, i);
      gsl_permute_vector(perm, &v.vector);
    }

  for (i = 0; i < N; i++)
    {
      for (j = 0; j < N; j++)
        {
          double Aij = gsl_matrix_get(A, i, j); /* L D L^T */
          double Bij = gsl_matrix_get(V, i, j); /* P M P^T */

          gsl_test_rel(Aij, Bij, eps,
                       "%s: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, Aij, Bij);
        }
    }

  if (expected_rcond > 0 && !scale)
    {
      gsl_vector *work = gsl_vector_alloc(3 * N);
      double rcond;

      gsl_linalg_pcholesky_rcond(LDLT, perm, &rcond, work);

      gsl_test_rel(rcond, expected_rcond, 1.0e-6,
                   "%s rcond: (%3lu,%3lu): %22.18g   %22.18g\n",
                   desc, N, N, rcond, expected_rcond);

      gsl_vector_free(work);
    }

  gsl_matrix_free(LDLT);
  gsl_matrix_free(V);
  gsl_matrix_free(A);
  gsl_matrix_free(L);
  gsl_matrix_free(LT);
  gsl_vector_free(S);
  gsl_permutation_free(perm);

  return s;
}

static int
test_pcholesky_decomp(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 50;
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);
      test_pcholesky_decomp_eps(0, m, -1.0, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp unscaled random");
      test_pcholesky_decomp_eps(1, m, -1.0, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp scaled random");

      if (N <= 12)
        {
          double expected_rcond = -1.0;

          if (hilb_rcond[N - 1] > 1.0e-12)
            expected_rcond = hilb_rcond[N - 1];

          create_hilbert_matrix2(m);

          test_pcholesky_decomp_eps(0, m, expected_rcond, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp unscaled hilbert");
          test_pcholesky_decomp_eps(1, m, expected_rcond, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp scaled hilbert");
        }

      gsl_matrix_free(m);
    }

  return s;
}

int
test_pcholesky_solve_eps(const int scale, const gsl_matrix * m, const gsl_vector * rhs,
                         const gsl_vector * sol, const double eps,
                         const char * desc)
{
  int s = 0;
  size_t i, N = m->size1;
  gsl_matrix * u  = gsl_matrix_alloc(N, N);
  gsl_vector * x = gsl_vector_calloc(N);
  gsl_vector * S = gsl_vector_alloc(N);
  gsl_permutation * perm = gsl_permutation_alloc(N);

  gsl_matrix_memcpy(u, m);

  if (scale)
    {
      s += gsl_linalg_pcholesky_decomp2(u, perm, S);
      s += gsl_linalg_pcholesky_solve2(u, perm, S, rhs, x);
    }
  else
    {
      s += gsl_linalg_pcholesky_decomp(u, perm);
      s += gsl_linalg_pcholesky_solve(u, perm, rhs, x);
    }

  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(sol, i);

      gsl_test_rel(xi, yi, eps,
                   "%s: %3lu[%lu]: %22.18g   %22.18g\n",
                   desc, N, i, xi, yi);
    }

  gsl_vector_free(x);
  gsl_vector_free(S);
  gsl_matrix_free(u);
  gsl_permutation_free(perm);

  return s;
}

static int
test_pcholesky_solve(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 50;
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix * m = gsl_matrix_alloc(N, N);
      gsl_vector * rhs = gsl_vector_alloc(N);
      gsl_vector * sol = gsl_vector_alloc(N);

      create_posdef_matrix(m, r);
      create_random_vector(sol, r);
      gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);

      test_pcholesky_solve_eps(0, m, rhs, sol, 64.0 * N * GSL_DBL_EPSILON, "pcholesky_solve unscaled random");
      test_pcholesky_solve_eps(1, m, rhs, sol, 64.0 * N * GSL_DBL_EPSILON, "pcholesky_solve scaled random");

      if (N <= 3)
        {
          create_hilbert_matrix2(m);
          gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);
          test_pcholesky_solve_eps(0, m, rhs, sol, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_solve unscaled hilbert");
          test_pcholesky_solve_eps(1, m, rhs, sol, 2048.0 * N * GSL_DBL_EPSILON, "pcholesky_solve scaled hilbert");
        }

      gsl_matrix_free(m);
      gsl_vector_free(rhs);
      gsl_vector_free(sol);
    }

  return s;
}

int
test_pcholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc)
{
  int s = 0;
  size_t i, j, N = m->size1;

  gsl_matrix * v  = gsl_matrix_alloc(N, N);
  gsl_matrix * c  = gsl_matrix_alloc(N, N);
  gsl_matrix * minv = gsl_matrix_alloc(N, N);
  gsl_permutation * p = gsl_permutation_alloc(N);

  gsl_matrix_memcpy(v, m);

  s += gsl_linalg_pcholesky_decomp(v, p);
  s += gsl_linalg_pcholesky_invert(v, p, minv);

  /* c = m m^{-1} */
  gsl_blas_dsymm(CblasLeft, CblasUpper, 1.0, m, minv, 0.0, c);

  /* c should be the identity matrix */

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          double cij = gsl_matrix_get(c, i, j);
          double expected = (i == j) ? 1.0 : 0.0;

          gsl_test_rel(cij, expected, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, cij, expected);
        }
    }

  gsl_matrix_free(v);
  gsl_matrix_free(c);
  gsl_matrix_free(minv);
  gsl_permutation_free(p);

  return s;
}

int
test_pcholesky_invert(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 30;
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);
      test_pcholesky_invert_eps(m, N * GSL_DBL_EPSILON, "pcholesky_invert unscaled random");

      if (N <= 4)
        {
          create_hilbert_matrix2(m);
          test_pcholesky_invert_eps(m, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_invert unscaled hilbert");
        }

      gsl_matrix_free(m);
    }

  return s;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE