text:
//Converts LSFs (line spectral frequencies) to polynomial representation along cols or rows of X.
//Works by getting p_sym and p_asym and then roots and then lsf2poly.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cblas.h>

#ifdef __cplusplus
namespace ov {
extern "C" {
#endif

int lsf2poly_s (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim);
int lsf2poly_d (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim);

int lsf2poly_s (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim)
{
    const float z = 0.0f, o = 1.0f;
    const int F = (dim==0) ? R : C;  //num LSFs
    const int nr = 2*F + 2;          //length roots
    const int P = F + 2;             //length psym and pasym
    const int Po = P - 1;            //num output poly coeffs
    int r, c, p, f;
    float *roots, *psym, *pasym;

    //Checks
    if (R<1) { fprintf(stderr,"error in lsf2poly_s: nrows X must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in lsf2poly_s: ncols X must be positive\n"); return 1; }

    //Initialize
    if (!(psym=(float *)malloc((size_t)(2*P)*sizeof(float)))) { fprintf(stderr,"error in lsf2poly_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(pasym=(float *)malloc((size_t)(2*P)*sizeof(float)))) { fprintf(stderr,"error in lsf2poly_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(roots=(float *)malloc((size_t)(2*nr)*sizeof(float)))) { fprintf(stderr,"error in lsf2poly_s: problem with malloc. "); perror("malloc"); return 1; }
    roots[0] = o; roots[nr] = -o; roots[1] = roots[nr+1] = z;

    if (dim==0)
    {
        if (iscolmajor)
        {
            //cblas_scopy(C,&o,0,Y,Po);
            for (c=0; c<C; c++)
            {
                for (f=0; f<F; f++)
                {
                    roots[2*f+2] = roots[4*F-2*f+2] = cosf(X[f+c*R]);
                    roots[2*f+3] = sinf(X[f+c*R]); roots[4*F-2*f+3] = -roots[2*f+3];
                }
                cblas_scopy(2*P,&z,0,psym,1); cblas_scopy(2*P,&z,0,pasym,1);
                psym[0] = pasym[0] = o;
                for (r=0; r<P-1; r++)
                {
                    for (p=r+1; p>0; p--)
                    {
                        pasym[2*p] = fmaf(-roots[4*r],pasym[2*p-2],fmaf(roots[4*r+1],pasym[2*p-1],pasym[2*p]));
                        pasym[2*p+1] = fmaf(-roots[4*r],pasym[2*p-1],fmaf(-roots[4*r+1],pasym[2*p-2],pasym[2*p+1]));
                        psym[2*p] = fmaf(-roots[4*r+2],psym[2*p-2],fmaf(roots[4*r+3],psym[2*p-1],psym[2*p]));
                        psym[2*p+1] = fmaf(-roots[4*r+2],psym[2*p-1],fmaf(-roots[4*r+3],psym[2*p-2],psym[2*p+1]));
                    }
                }
                for (r=0; r<Po; r++) { Y[r+c*Po] = 0.5f*(psym[2*(Po-r)]-pasym[2*(Po-r)]); }
            }
        }
        else
        {
            //cblas_scopy(C,&o,0,Y,1);
            for (c=0; c<C; c++)
            {
                for (f=0; f<F; f++)
                {
                    roots[2*f+2] = roots[4*F-2*f+2] = cosf(X[c+f*C]);
                    roots[2*f+3] = sinf(X[c+f*C]); roots[4*F-2*f+3] = -roots[2*f+3];
                }
                cblas_scopy(2*P,&z,0,psym,1); cblas_scopy(2*P,&z,0,pasym,1);
                psym[0] = pasym[0] = o;
                for (r=0; r<P-1; r++)
                {
                    for (p=r+1; p>0; p--)
                    {
                        pasym[2*p] = fmaf(-roots[4*r],pasym[2*p-2],fmaf(roots[4*r+1],pasym[2*p-1],pasym[2*p]));
                        pasym[2*p+1] = fmaf(-roots[4*r],pasym[2*p-1],fmaf(-roots[4*r+1],pasym[2*p-2],pasym[2*p+1]));
                        psym[2*p] = fmaf(-roots[4*r+2],psym[2*p-2],fmaf(roots[4*r+3],psym[2*p-1],psym[2*p]));
                        psym[2*p+1] = fmaf(-roots[4*r+2],psym[2*p-1],fmaf(-roots[4*r+3],psym[2*p-2],psym[2*p+1]));
                    }
                }
                for (r=0; r<Po; r++) { Y[c+r*C] = 0.5f*(psym[2*(Po-r)]-pasym[2*(Po-r)]); }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            //cblas_scopy(R,&o,0,Y,1);
            for (r=0; r<R; r++)
            {
                for (f=0; f<F; f++)
                {
                    roots[2*f+2] = roots[4*F-2*f+2] = cosf(X[r+f*R]);
                    roots[2*f+3] = sinf(X[r+f*R]); roots[4*F-2*f+3] = -roots[2*f+3];
                }
                cblas_scopy(2*P,&z,0,psym,1); cblas_scopy(2*P,&z,0,pasym,1);
                psym[0] = pasym[0] = o;
                for (c=0; c<P-1; c++)
                {
                    for (p=c+1; p>0; p--)
                    {
                        pasym[2*p] = fmaf(-roots[4*c],pasym[2*p-2],fmaf(roots[4*c+1],pasym[2*p-1],pasym[2*p]));
                        pasym[2*p+1] = fmaf(-roots[4*c],pasym[2*p-1],fmaf(-roots[4*c+1],pasym[2*p-2],pasym[2*p+1]));
                        psym[2*p] = fmaf(-roots[4*c+2],psym[2*p-2],fmaf(roots[4*c+3],psym[2*p-1],psym[2*p]));
                        psym[2*p+1] = fmaf(-roots[4*c+2],psym[2*p-1],fmaf(-roots[4*c+3],psym[2*p-2],psym[2*p+1]));
                    }
                }
                for (c=0; c<Po; c++) { Y[r+c*R] = 0.5f*(psym[2*(Po-c)]-pasym[2*(Po-c)]); }
            }
        }
        else
        {
            //cblas_scopy(R,&o,0,Y,Po);
            for (r=0; r<R; r++)
            {
                for (f=0; f<F; f++)
                {
                    roots[2*f+2] = roots[4*F-2*f+2] = cosf(X[f+r*C]);
                    roots[2*f+3] = sinf(X[f+r*C]); roots[4*F-2*f+3] = -roots[2*f+3];
                }
                cblas_scopy(2*P,&z,0,psym,1); cblas_scopy(2*P,&z,0,pasym,1);
                psym[0] = pasym[0] = o;
                for (c=0; c<P-1; c++)
                {
                    for (p=c+1; p>0; p--)
                    {
                        pasym[2*p] = fmaf(-roots[4*c],pasym[2*p-2],fmaf(roots[4*c+1],pasym[2*p-1],pasym[2*p]));
                        pasym[2*p+1] = fmaf(-roots[4*c],pasym[2*p-1],fmaf(-roots[4*c+1],pasym[2*p-2],pasym[2*p+1]));
                        psym[2*p] = fmaf(-roots[4*c+2],psym[2*p-2],fmaf(roots[4*c+3],psym[2*p-1],psym[2*p]));
                        psym[2*p+1] = fmaf(-roots[4*c+2],psym[2*p-1],fmaf(-roots[4*c+3],psym[2*p-2],psym[2*p+1]));
                    }
                }
                for (c=0; c<Po; c++) { Y[c+r*Po] = 0.5f*(psym[2*(Po-c)]-pasym[2*(Po-c)]); }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in lsf2poly_s: dim must be 0 or 1.\n"); return 1;
    }

    //Exit
    free(roots); free(psym); free(pasym);
    return 0;
}

int lsf2poly_d (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim)
{
    const double z = 0.0, o = 1.0;
    const int F = (dim==0) ? R : C;  //num LSFs
    const int nr = 2*F + 2;          //length roots
    const int P = F + 2;             //length psym and pasym
    const int Po = P - 1;            //num output poly coeffs
    int r, c, p, f;
    double *roots, *psym, *pasym;

    //Checks
    if (R<1) { fprintf(stderr,"error in lsf2poly_d: nrows X must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in lsf2poly_d: ncols X must be positive\n"); return 1; }

    //Initialize
    if (!(psym=(double *)malloc((size_t)(2*P)*sizeof(double)))) { fprintf(stderr,"error in lsf2poly_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(pasym=(double *)malloc((size_t)(2*P)*sizeof(double)))) { fprintf(stderr,"error in lsf2poly_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(roots=(double *)malloc((size_t)(2*nr)*sizeof(double)))) { fprintf(stderr,"error in lsf2poly_d: problem with malloc. "); perror("malloc"); return 1; }
    roots[0] = o; roots[nr] = -o; roots[1] = roots[nr+1] = z;

    if (dim==0)
    {
        if (iscolmajor)
        {
            for (c=0; c<C; c++)
            {
                for (f=0; f<F; f++)
                {
                    roots[2*f+2] = roots[4*F-2*f+2] = cos(X[f+c*R]);
                    roots[2*f+3] = sin(X[f+c*R]); roots[4*F-2*f+3] = -roots[2*f+3];
                }
                cblas_dcopy(2*P,&z,0,psym,1); cblas_dcopy(2*P,&z,0,pasym,1);
                psym[0] = pasym[0] = o;
                for (r=0; r<P-1; r++)
                {
                    for (p=r+1; p>0; p--)
                    {
                        pasym[2*p] = fma(-roots[4*r],pasym[2*p-2],fma(roots[4*r+1],pasym[2*p-1],pasym[2*p]));
                        pasym[2*p+1] = fma(-roots[4*r],pasym[2*p-1],fma(-roots[4*r+1],pasym[2*p-2],pasym[2*p+1]));
                        psym[2*p] = fma(-roots[4*r+2],psym[2*p-2],fma(roots[4*r+3],psym[2*p-1],psym[2*p]));
                        psym[2*p+1] = fma(-roots[4*r+2],psym[2*p-1],fma(-roots[4*r+3],psym[2*p-2],psym[2*p+1]));
                    }
                }
                for (r=0; r<Po; r++) { Y[r+c*Po] = 0.5*(psym[2*(Po-r)]-pasym[2*(Po-r)]); }
            }
        }
        else
        {
            for (c=0; c<C; c++)
            {
                for (f=0; f<F; f++)
                {
                    roots[2*f+2] = roots[4*F-2*f+2] = cos(X[c+f*C]);
                    roots[2*f+3] = sin(X[c+f*C]); roots[4*F-2*f+3] = -roots[2*f+3];
                }
                cblas_dcopy(2*P,&z,0,psym,1); cblas_dcopy(2*P,&z,0,pasym,1);
                psym[0] = pasym[0] = o;
                for (r=0; r<P-1; r++)
                {
                    for (p=r+1; p>0; p--)
                    {
                        pasym[2*p] = fma(-roots[4*r],pasym[2*p-2],fma(roots[4*r+1],pasym[2*p-1],pasym[2*p]));
                        pasym[2*p+1] = fma(-roots[4*r],pasym[2*p-1],fma(-roots[4*r+1],pasym[2*p-2],pasym[2*p+1]));
                        psym[2*p] = fma(-roots[4*r+2],psym[2*p-2],fma(roots[4*r+3],psym[2*p-1],psym[2*p]));
                        psym[2*p+1] = fma(-roots[4*r+2],psym[2*p-1],fma(-roots[4*r+3],psym[2*p-2],psym[2*p+1]));
                    }
                }
                for (r=0; r<Po; r++) { Y[c+r*C] = 0.5*(psym[2*(Po-r)]-pasym[2*(Po-r)]); }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (r=0; r<R; r++)
            {
                for (f=0; f<F; f++)
                {
                    roots[2*f+2] = roots[4*F-2*f+2] = cos(X[r+f*R]);
                    roots[2*f+3] = sin(X[r+f*R]); roots[4*F-2*f+3] = -roots[2*f+3];
                }
                cblas_dcopy(2*P,&z,0,psym,1); cblas_dcopy(2*P,&z,0,pasym,1);
                psym[0] = pasym[0] = o;
                for (c=0; c<P-1; c++)
                {
                    for (p=c+1; p>0; p--)
                    {
                        pasym[2*p] = fma(-roots[4*c],pasym[2*p-2],fma(roots[4*c+1],pasym[2*p-1],pasym[2*p]));
                        pasym[2*p+1] = fma(-roots[4*c],pasym[2*p-1],fma(-roots[4*c+1],pasym[2*p-2],pasym[2*p+1]));
                        psym[2*p] = fma(-roots[4*c+2],psym[2*p-2],fma(roots[4*c+3],psym[2*p-1],psym[2*p]));
                        psym[2*p+1] = fma(-roots[4*c+2],psym[2*p-1],fma(-roots[4*c+3],psym[2*p-2],psym[2*p+1]));
                    }
                }
                for (c=0; c<Po; c++) { Y[r+c*R] = 0.5*(psym[2*(Po-c)]-pasym[2*(Po-c)]); }
            }
        }
        else
        {
            for (r=0; r<R; r++)
            {
                for (f=0; f<F; f++)
                {
                    roots[2*f+2] = roots[4*F-2*f+2] = cos(X[f+r*C]);
                    roots[2*f+3] = sin(X[f+r*C]); roots[4*F-2*f+3] = -roots[2*f+3];
                }
                cblas_dcopy(2*P,&z,0,psym,1); cblas_dcopy(2*P,&z,0,pasym,1);
                psym[0] = pasym[0] = o;
                for (c=0; c<P-1; c++)
                {
                    for (p=c+1; p>0; p--)
                    {
                        pasym[2*p] = fma(-roots[4*c],pasym[2*p-2],fma(roots[4*c+1],pasym[2*p-1],pasym[2*p]));
                        pasym[2*p+1] = fma(-roots[4*c],pasym[2*p-1],fma(-roots[4*c+1],pasym[2*p-2],pasym[2*p+1]));
                        psym[2*p] = fma(-roots[4*c+2],psym[2*p-2],fma(roots[4*c+3],psym[2*p-1],psym[2*p]));
                        psym[2*p+1] = fma(-roots[4*c+2],psym[2*p-1],fma(-roots[4*c+3],psym[2*p-2],psym[2*p+1]));
                    }
                }
                for (c=0; c<Po; c++) { Y[c+r*Po] = 0.5*(psym[2*(Po-c)]-pasym[2*(Po-c)]); }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in lsf2poly_d: dim must be 0 or 1.\n"); return 1;
    }

    //Exit
    free(roots); free(psym); free(pasym);
    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
/*

Excited States software: KGS
Contributors: See CONTRIBUTORS.txt
Contact: kgs-contact@simtk.org

Copyright (C) 2009-2017 Stanford University

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
of the Software, and to permit persons to whom the Software is furnished to do
so, subject to the following conditions:

This entire text, including the above copyright notice and this permission notice
shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
IN THE SOFTWARE.

*/

#ifndef JACOBIANRELATED_H
#define JACOBIANRELATED_H

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multimin.h>

///These values have to be chosen according to the numerical analysis
static const double SINGVAL_TOL = 0.000000000001; // only generic 10^-12
//static const double SINGVAL_TOL = 0.000001; // include non-generic 10^-10
static const double RIGID_TOL = 0.0000000001; //depends on molecule, but 10^-10 seems a good fit!
//static const double RIGID_TOL = 0.000001; //depends on molecule ToDo: Adapt to obtain on-the-fly parameter selection

class NullSpaceRet
{
public:
//	static gsl_matrix* Ut;//The Ut from the SVD
	static gsl_matrix* V;//The V from the SVD
	static gsl_vector* singularValues;
	int nullspaceSize;
	int m, n;

//	gsl_matrix* P;//Null-space projection matrix
	gsl_matrix* m_nullspaceBasis; //basis of the nullspace
//	gsl_matrix* Jinv;//Jacobian pseudo-inverse

    gsl_vector* m_rigidAngles;
    gsl_vector* m_rigidHBonds;

    int m_numCoordinated;
    int m_numRigid;
    int m_numRigidHBonds;

	NullSpaceRet (int input_m, int input_n);
	~NullSpaceRet ();
	void print();

	// A function to make the projection operation more efficient;
	void NullSpacePostCompute();

	/// A function that identifies rigidified dihedral angles
	void RigidityAnalysis(gsl_matrix* HBondJacobian);

	//A function to project a vector on the nullspace
	void ProjectOnNullSpace (gsl_vector *to_project, gsl_vector *after_project);

	// To manipulate owner configations of NullSpaceRet instances
	void numOwners(unsigned int numOwners);
	unsigned int numOwners() const;
private:
	unsigned int numOwners_;
};

void nr_svd(gsl_matrix* a_g, int m, int n, double w[], gsl_matrix* v_out);
void ComputeNullSpace(gsl_matrix* Jac, NullSpaceRet* Ret);
//void ProjectOnNullSpace (NullSpaceRet* nullspace, gsl_vector *e, gsl_vector *to_project, gsl_vector *after_project);
void TorsionUpdate(NullSpaceRet* nullspace, gsl_vector *e, double lambda, gsl_vector *dTheta);

#endif

ext: h


text:
#ifndef PRECISION_RECALL_H_G2FJDYVV
#define PRECISION_RECALL_H_G2FJDYVV

#include <gsl/gsl>
#include <opencv2/core/base.hpp>
#include <opencv2/core/types.hpp>
#include <opencv2/imgproc.hpp>
#include <optional>
#include <string_view>
#include <util/common_structures.h>

namespace sens_loc::apps {

/// Defines the color and line strength used for backprojection.
struct backproject_style {
    backproject_style() = default;
    backproject_style(int r, int g, int b, int strength)
        : color{CV_RGB(r, g, b)}
        , strength{strength} {
        Expects(strength > 0);
        Expects(r >= 0);
        Expects(g >= 0);
        Expects(b >= 0);
        Expects(r < 256);
        Expects(g < 256);
        Expects(b < 256);
    }

    cv::Scalar color    = CV_RGB(0, 0, 0);
    int        strength = 1;
};

/// Composition of styles on how the backprojections are plotted.
/// \sa backproject_style
struct backproject_config {
    backproject_style true_positive;
    backproject_style false_negative;
    backproject_style false_positive;
};

struct recognition_analysis_input {
    std::string_view                depth_image_pattern;
    std::string_view                pose_file_pattern;
    std::string_view                intrinsic_file;
    std::optional<std::string_view> mask_file;
    cv::NormTypes                   matching_norm;
    float                           keypoint_distance_threshold;

    /// Unit-Conversion for the ICP of kinect images. No other depth images
    /// are treated with ICP, so this is dirty set to a constant.
    const float unit_factor = 0.001F;
};

struct recognition_analysis_output_options {
    std::optional<std::string> backproject_pattern;
    std::optional<std::string> original_files;
    std::optional<std::string> stat_file;
    std::optional<std::string> backprojection_selected_histo;
    std::optional<std::string> relevant_histo;
    std::optional<std::string> true_positive_histo;
    std::optional<std::string> false_positive_histo;
    std::optional<std::string> true_positive_distance_histo;
    std::optional<std::string> false_positive_distance_histo;
};

int analyze_recognition_performance(
    util::processing_input                     in,
    const recognition_analysis_input&          required_data,
    const recognition_analysis_output_options& output_options,
    const backproject_config&                  backproject_config);

}  // namespace sens_loc::apps

#endif /* end of include guard: PRECISION_RECALL_H_G2FJDYVV */

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_pow_int.h>
#include "bessel.h"
#include "bessel_amp_phase.h"
#include "bessel_olver.h"
#include <gsl/gsl_sf_bessel.h>

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_bessel_Jn_e(int n, double x, gsl_sf_result * result)
{
  int sign = 1;

  if(n < 0) {
    /* reduce to case n >= 0 */
    n = -n;
    if(GSL_IS_ODD(n)) sign = -sign;
  }

  if(x < 0.0) {
    /* reduce to case x >= 0. */
    x = -x;
    if(GSL_IS_ODD(n)) sign = -sign;
  }

  /* CHECK_POINTER(result) */

  if(n == 0) {
    gsl_sf_result b0;
    int stat_J0 = gsl_sf_bessel_J0_e(x, &b0);
    result->val = sign * b0.val;
    result->err = b0.err;
    return stat_J0;
  }
  else if(n == 1) {
    gsl_sf_result b1;
    int stat_J1 = gsl_sf_bessel_J1_e(x, &b1);
    result->val = sign * b1.val;
    result->err = b1.err;
    return stat_J1;
  }
  else {
    if(x == 0.0) {
      result->val = 0.0;
      result->err = 0.0;
      return GSL_SUCCESS;
    }
    else if(x*x < 10.0*(n+1.0)*GSL_ROOT5_DBL_EPSILON) {
      gsl_sf_result b;
      int status = gsl_sf_bessel_IJ_taylor_e((double)n, x, -1, 50, GSL_DBL_EPSILON, &b);
      result->val  = sign * b.val;
      result->err  = b.err;
      result->err += GSL_DBL_EPSILON * fabs(result->val);
      return status;
    }
    else if(GSL_ROOT4_DBL_EPSILON * x > (n*n+1.0)) {
      int status = gsl_sf_bessel_Jnu_asympx_e((double)n, x, result);
      result->val *= sign;
      return status;
    }
    else if(n > 50) {
      int status = gsl_sf_bessel_Jnu_asymp_Olver_e((double)n, x, result);
      result->val *= sign;
      return status;
    }
    else if(x > 1000.0)
    {
      /* We need this to avoid feeding large x to CF1; note that
       * due to the above check, we know that n <= 50.
       */
      int status = gsl_sf_bessel_Jnu_asympx_e((double)n, x, result);
      result->val *= sign;
      return status;
    }
    else {
      double ans;
      double err;
      double ratio;
      double sgn;
      int stat_b;
      int stat_CF1 = gsl_sf_bessel_J_CF1((double)n, x, &ratio, &sgn);

      /* backward recurrence */
      double Jkp1 = GSL_SQRT_DBL_MIN * ratio;
      double Jk   = GSL_SQRT_DBL_MIN;
      double Jkm1;
      int k;

      for(k=n; k>0; k--) {
        Jkm1 = 2.0*k/x * Jk - Jkp1;
        Jkp1 = Jk;
        Jk   = Jkm1;
      }

      if(fabs(Jkp1) > fabs(Jk)) {
        gsl_sf_result b1;
        stat_b = gsl_sf_bessel_J1_e(x, &b1);
        ans = b1.val/Jkp1 * GSL_SQRT_DBL_MIN;
        err = b1.err/Jkp1 * GSL_SQRT_DBL_MIN;
      }
      else {
        gsl_sf_result b0;
        stat_b = gsl_sf_bessel_J0_e(x, &b0);
        ans = b0.val/Jk * GSL_SQRT_DBL_MIN;
        err = b0.err/Jk * GSL_SQRT_DBL_MIN;
      }

      result->val = sign * ans;
      result->err = fabs(err);
      return GSL_ERROR_SELECT_2(stat_CF1, stat_b);
    }
  }
}

int
gsl_sf_bessel_Jn_array(int nmin, int nmax, double x, double * result_array)
{
  /* CHECK_POINTER(result_array) */

  if(nmin < 0 || nmax < nmin) {
    int n;
    for(n=nmax; n>=nmin; n--) {
      result_array[n-nmin] = 0.0;
    }
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else if(x == 0.0) {
    int n;
    for(n=nmax; n>=nmin; n--) {
      result_array[n-nmin] = 0.0;
    }
    if(nmin == 0) result_array[0] = 1.0;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result r_Jnp1;
    gsl_sf_result r_Jn;
    int stat_np1 = gsl_sf_bessel_Jn_e(nmax+1, x, &r_Jnp1);
    int stat_n   = gsl_sf_bessel_Jn_e(nmax,   x, &r_Jn);
    int stat = GSL_ERROR_SELECT_2(stat_np1, stat_n);

    double Jnp1 = r_Jnp1.val;
    double Jn   = r_Jn.val;
    double Jnm1;
    int n;

    if(stat == GSL_SUCCESS) {
      for(n=nmax; n>=nmin; n--) {
        result_array[n-nmin] = Jn;
        Jnm1 = -Jnp1 + 2.0*n/x * Jn;
        Jnp1 = Jn;
        Jn   = Jnm1;
      }
    }
    else {
      for(n=nmax; n>=nmin; n--) {
        result_array[n-nmin] = 0.0;
      }
    }

    return stat;
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_Jn(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_bessel_Jn_e(n, x, &result));
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include "../EngineTest.h"

#include <Catch2>

#include <gsl/gsl_math.h>

TEST_CASE("kurtosis Function Evaluation Tests", "[kurtosis]") {

    // SECTION("Empty Test"){
    //     requireIsEqual("kurtosis()", "Insufficient Number of Arguments for Function: kurtosis");
    // }

    SECTION("`kurtosis` Test 1"){
        requireIsEqual("kurtosis(-9.75, -6.89, 8.71, 9.56)", -2.4054875147);
    }

    SECTION("`kurtosis` Test 2"){
        requireIsEqual("kurtosis(-0.03, -0.66, 8.08, 7.46, -9.47, 5.48)", -1.3867096466);
    }

    SECTION("`kurt` Test 3"){
        requireIsEqual("kurt(6.88, 2.73, 0.91, 7.03)", -2.3097016758);
    }

    SECTION("`kurt` Test 4"){
        requireIsEqual("kurt(8.51, 4.35)", -2.75);
    }

    SECTION("`kurtosis` Test 5"){
        requireIsEqual("kurtosis(9.24, -5.27, 8.46, 8.69, 3.49, 7.69)", -0.7621422371);
    }

    SECTION("`kurt` Test 6"){
        requireIsEqual("kurt(0.94, 2.67, -0.5, -4.74)", -1.8915463404);
    }

    SECTION("`kurt` Test 7"){
        requireIsEqual("kurt(-3.29, -3.04, 2.81, 6.5, -2.0)", -1.8864498289);
    }

    SECTION("`kurtosis` Test 8"){
        requireIsEqual("kurtosis(-9.95, -8.32, 3.62, -9.95)", -1.7095361309);
    }

    SECTION("`kurtosis` Test 9"){
        requireIsEqual("kurtosis(-2.47, 7.98, 7.08, 5.36, -0.13, 0.57, -5.41, -6.09, -4.31)", -1.7566598859);
    }

    SECTION("`kurtosis` Test 10"){
        requireIsEqual("kurtosis(-3.21, -8.75, 6.86, 9.35, -3.09, 4.99)", -1.932477903);
    }

}

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
/**
 *
 * @file qwrapper_sgemm.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:56 2014
 *
 **/
#include <cblas.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_sgemm(Quark *quark, Quark_Task_Flags *task_flags,
                      PLASMA_enum transA, int transB,
                      int m, int n, int k, int nb,
                      float alpha, const float *A, int lda,
                                                const float *B, int ldb,
                      float beta, float *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_sgemm_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(float),         &alpha,     VALUE,
        sizeof(float)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(float)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(float),         &beta,      VALUE,
        sizeof(float)*nb*nb,    C,                 INOUT,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_sgemm2( Quark *quark, Quark_Task_Flags *task_flags,
                        PLASMA_enum transA, int transB,
                        int m, int n, int k, int nb,
                        float alpha, const float *A, int lda,
                        const float *B, int ldb,
                        float beta, float *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_sgemm_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(float),         &alpha,     VALUE,
        sizeof(float)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(float)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(float),         &beta,      VALUE,
        sizeof(float)*nb*nb,    C,                 INOUT | LOCALITY | GATHERV,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_sgemm_quark = PCORE_sgemm_quark
#define CORE_sgemm_quark PCORE_sgemm_quark
#endif
void CORE_sgemm_quark(Quark *quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int m;
    int n;
    int k;
    float alpha;
    float *A;
    int lda;
    float *B;
    int ldb;
    float beta;
    float *C;
    int ldc;

    quark_unpack_args_13(quark, transA, transB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
    cblas_sgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        m, n, k,
        (alpha), A, lda,
        B, ldb,
        (beta), C, ldc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_sgemm_f2(Quark *quark, Quark_Task_Flags *task_flags,
                         PLASMA_enum transA, int transB,
                         int m, int n, int k, int nb,
                         float alpha, const float *A, int lda,
                                                   const float *B, int ldb,
                         float beta, float *C, int ldc,
                         float *fake1, int szefake1, int flag1,
                         float *fake2, int szefake2, int flag2)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_sgemm_f2_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(float),         &alpha,     VALUE,
        sizeof(float)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(float)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(float),         &beta,      VALUE,
        sizeof(float)*nb*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        sizeof(float)*szefake1, fake1,             flag1,
        sizeof(float)*szefake2, fake2,             flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_sgemm_f2_quark = PCORE_sgemm_f2_quark
#define CORE_sgemm_f2_quark PCORE_sgemm_f2_quark
#endif
void CORE_sgemm_f2_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    float alpha;
    float *A;
    int LDA;
    float *B;
    int LDB;
    float beta;
    float *C;
    int LDC;
    void *fake1, *fake2;

    quark_unpack_args_15(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC, fake1, fake2);
    cblas_sgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        (alpha), A, LDA,
        B, LDB,
        (beta), C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_sgemm_p2(Quark *quark, Quark_Task_Flags *task_flags,
                         PLASMA_enum transA, int transB,
                         int m, int n, int k, int nb,
                         float alpha, const float *A, int lda,
                         const float **B, int ldb,
                         float beta, float *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_sgemm_p2_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(float),         &alpha,     VALUE,
        sizeof(float)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(float*),         B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(float),         &beta,      VALUE,
        sizeof(float)*ldc*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_sgemm_p2_quark = PCORE_sgemm_p2_quark
#define CORE_sgemm_p2_quark PCORE_sgemm_p2_quark
#endif
void CORE_sgemm_p2_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    float alpha;
    float *A;
    int LDA;
    float **B;
    int LDB;
    float beta;
    float *C;
    int LDC;

    quark_unpack_args_13(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC);
    cblas_sgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        (alpha), A, LDA,
        *B, LDB,
        (beta), C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_sgemm_p3(Quark *quark, Quark_Task_Flags *task_flags,
                           PLASMA_enum transA, int transB,
                           int m, int n, int k, int nb,
                           float alpha, const float *A, int lda,
                           const float *B, int ldb,
                           float beta, float **C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_sgemm_p3_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(float),         &alpha,     VALUE,
        sizeof(float)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(float)*ldb*nb,   B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(float),         &beta,      VALUE,
        sizeof(float*),         C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_sgemm_p3_quark = PCORE_sgemm_p3_quark
#define CORE_sgemm_p3_quark PCORE_sgemm_p3_quark
#endif
void CORE_sgemm_p3_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    float alpha;
    float *A;
    int LDA;
    float *B;
    int LDB;
    float beta;
    float **C;
    int LDC;

    quark_unpack_args_13(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC);
    cblas_sgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        (alpha), A, LDA,
        B, LDB,
        (beta), *C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_sgemm_p2f1(Quark *quark, Quark_Task_Flags *task_flags,
                           PLASMA_enum transA, int transB,
                           int m, int n, int k, int nb,
                           float alpha, const float *A, int lda,
                           const float **B, int ldb,
                           float beta, float *C, int ldc,
                           float *fake1, int szefake1, int flag1)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_sgemm_p2f1_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(float),         &alpha,     VALUE,
        sizeof(float)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(float*),         B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(float),         &beta,      VALUE,
        sizeof(float)*ldc*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        sizeof(float)*szefake1, fake1,             flag1,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_sgemm_p2f1_quark = PCORE_sgemm_p2f1_quark
#define CORE_sgemm_p2f1_quark PCORE_sgemm_p2f1_quark
#endif
void CORE_sgemm_p2f1_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    float alpha;
    float *A;
    int LDA;
    float **B;
    int LDB;
    float beta;
    float *C;
    int LDC;
    void *fake1;

    quark_unpack_args_14(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC, fake1);
    cblas_sgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        (alpha), A, LDA,
        *B, LDB,
        (beta), C, LDC);
}


text:
/*
Copyright (c) 2015, Patrick Weltevrede
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics_double.h>
#include "psrsalsa.h"
#define KSTEST 1
#define KSFLAT 2
#define KSSIN 3
#define KSFLAT_WITH_MINMAX 4
#define CORREL 20
#define PEARSON 21
#define MOMENTS 30
#define MEDIAN 31
#define CHI2TEST_HIST 50
#define CHI2TEST_CDF 51
int main(int argc, char **argv)
{
  psrsalsaApplication application;
  long i, j;
  int file1_column1, file1_column2, file1_column3, file2_column1, file2_column2, file2_column3, typetest, read_log, output_idx, skiplines;
  double threshold1, threshold2, threshold3, ksflat_min, ksflat_max;
  initApplication(&application, "pstat", "[options] inputfile(s)");
  application.switch_libversions = 1;
  application.switch_verbose = 1;
  application.switch_debug = 1;
  file1_column1 = 0;
  file1_column2 = 0;
  file1_column3 = 0;
  file2_column1 = 0;
  file2_column2 = 0;
  file2_column3 = 0;
  typetest = 0;
  read_log = 0;
  threshold1 = 0;
  threshold2 = 0;
  threshold3 = 0;
  output_idx = 0;
  skiplines = 0;
  if(argc < 2) {
    printf("Program to perform various statistical tests on input data. One or two input\n");
    printf("files are required depending on the statistical test to be performed. The input\n");
    printf("files should be ascii files with each line having an equal number of columns.\n");
    printf("Lines starting with a # will be ignored. Usage:\n\n");
    printApplicationHelp(&application);
    printf("Input/Output options:\n");
    printf("-col \"c1 c2 c3\"  Specify one, two or three column numbers (counting from 1) to\n");
    printf("                 be read in from the first input file. Some statistical tests\n");
    printf("                 require only one input column, others two or more.\n");
    printf("-col1            Equivalent to -col\n");
    printf("-col2            Like -col1, but for second input file.\n");
    printf("-log             The base 10 log of the input values is used.\n");
    printf("-output          Specify output filename to use rather than the stdout.\n");
    printf("-skiplines nr    Specify number of lines that should be skipped from start file(s).\n");
    printf("\nStatistical tests:\n");
    printf("-chi2hist \"t1 t2 t3\" Input should be two files, each being a histogram and\n");
    printf("                     having two (or three, see below) columns: bin location and\n");
    printf("                 the height. The histograms should have equal bin widths,\n");
    printf("                 overlap at least partially and have aligned bins. Only bins\n");
    printf("                 in the 1st input file with a height >= t1 are considered (so if\n");
    printf("                 negative, everything should be included. t2 is the theshold for\n");
    printf("                 the 2nd file, while t3 is that for the sum of heights of the\n");
    printf("                 two histograms. Specifying t2 and t3 is optional. The optional\n");
    printf("                 3rd column is the std. dev. on each bin, which is set to zero\n");
    printf("                 for bins outside the covered range of the histograms, and the\n");
    printf("                 height is set to zero as well. With only two columns, the\n");
    printf("                 chi-square test is unweighted, i.e. each bin has an equal\n");
    printf("                 weight (set to 1) regardless of the height of the bin(s).\n");
    printf("                 Examples: \n");
    printf("                 pstat -chi2hist \"-1 -1 0.1\" -col1 \"1 2\" -col2 \"1 2\" file1 file2\n");
    printf("                 pstat -chi2hist -1 -col1 \"1 2 3\" -col2 \"1 2 3\" file1 file2\n");
    printf("-chi2cdf         This test is similar to -chi2hist, except that the inputs are\n");
    printf("                 not histograms, but the distribution of points itself (i.e. one\n");
    printf("                 column). This test therefore does not depend on binning. The\n");
    printf("                 vertical difference in their CDF is quantified with an non-\n");
    printf("                 weighted chi-square test (weights are set to 1). This is done\n");
    printf("                 for each point in the CDF of the shortest distribution and\n");
    printf("                 using a linear interpolation of the other. This test can be\n");
    printf("                 useful to decide which model distribution fits an observed\n");
    printf("                 distribution best, but the chi square itself is not immediately\n");
    printf("                 meaningful. Example: pstat -chi2cdf file1 file2\n");
    printf("-correl          Produces the cross correlation function (as function of lag)\n");
    printf("                 calculated in the Fourier domain assuming equal sampling.\n");
    printf("                 Examples: pstat -correl -col1 1 -col2 1 file1 file2\n");
    printf("                 or:       pstat -correl -col \"1 2\" file1\n");
    printf("-ks              Kolmogorov-Smirnov test: Assess difference between two\n");
    printf("                 distributions of points. The made approximations are similar to\n");
    printf("                 the implementation in Numerical Recipes in C, 2nd edition and.\n");
    printf("                 does not depend on any binning.\n");
    printf("                 Examples: pstat -ks -col1 1 -col2 1 file1 file2\n");
    printf("                 or:       pstat -ks -col \"1 2\" file1\n");
    printf("-ksflat          Like -ks, but compare single input distribution with a flat\n");
    printf("                 distribution covering the range of input values.\n");
    printf("                 Example:  pstat -ksflat -col 1 file1\n");
    printf("-ksflat2 \"min max\"  Like -ksflat, but set the range to the specified values.\n");
    printf("                 Example:  pstat -ksflat2 \"0 1\" -col 1 file1\n");
    printf("-kssin           Like -ks, but compare single input distribution with a\n");
    printf("                 sinusoidal distribution between 0 and 90 deg.\n");
    printf("                 Example:  pstat -kssin -col 1 file1\n");
    printf("-median          Compute the median of the distribution.\n");
    printf("                 Example:  pstat -median -col 1 file1\n");
    printf("-moments         Compute different moments of the distribution (mean, variance\n");
    printf("                 etc.)\n");
    printf("                 Example:  pstat -moments -col 1 file1\n");
    printf("-pearson         Pearson product-moment correlation coefficient:\n");
    printf("                 Measures the linear correlation between two variables.\n");
    printf("                 Examples: pstat -pearson -col1 1 -col2 1 file1 file2\n");
    printf("                 or:       pstat -pearson -col \"1 2\" file1\n");
    printf("\n");
    printCitationInfo();
    terminateApplication(&application);
    return 0;
  }else {
    for(i = 1; i < argc; i++) {
      int index;
      index = i;
      if(processCommandLine(&application, argc, argv, &index)) {
 i = index;
      }else if(strcmp(argv[i], "-col") == 0 || strcmp(argv[i], "-col1") == 0) {
 int ret;
 ret = parse_command_string(application.verbose_state, argc, argv, i+1, 0, 1, "%d %d %d", &file1_column1, &file1_column2, &file1_column3, NULL);
 if(ret == 1) {
   file1_column2 = 0;
   file1_column3 = 0;
 }else if(ret == 2) {
   file1_column3 = 0;
 }else if(ret != 3) {
   printerror(application.verbose_state.debug, "ERROR pstat: Cannot parse '%s' option, 1, 2 or 3 values.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-col2") == 0) {
 int ret;
 ret = parse_command_string(application.verbose_state, argc, argv, i+1, 0, 1, "%d %d %d", &file2_column1, &file2_column2, &file2_column3, NULL);
 if(ret == 1) {
   file2_column2 = 0;
   file2_column3 = 0;
 }else if(ret == 2) {
   file2_column3 = 0;
 }else if(ret != 3) {
   printerror(application.verbose_state.debug, "ERROR pstat: Cannot parse '%s' option, 1, 2 or 3 values.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-skiplines") == 0) {
 int ret;
 ret = parse_command_string(application.verbose_state, argc, argv, i+1, 0, 1, "%d", &skiplines, NULL);
 if(ret != 1) {
   printerror(application.verbose_state.debug, "ERROR pstat: Cannot parse '%s' option, 1 integer value was expected.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-output") == 0) {
 output_idx = i+1;
 i++;
      }else if(strcmp(argv[i], "-correl") == 0) {
 if(typetest != 0) {
   printerror(application.verbose_state.debug, "pstat: Cannot specify more than one type of statistical test at the time");
   return 0;
 }
 typetest = CORREL;
      }else if(strcasecmp(argv[i], "-pearson") == 0) {
 if(typetest != 0) {
   printerror(application.verbose_state.debug, "pstat: Cannot specify more than one type of statistical test at the time");
   return 0;
 }
 typetest = PEARSON;
      }else if(strcmp(argv[i], "-ks") == 0) {
 if(typetest != 0) {
   printerror(application.verbose_state.debug, "pstat: Cannot specify more than one type of statistical test at the time");
   return 0;
 }
 typetest = KSTEST;
      }else if(strcmp(argv[i], "-ksflat") == 0) {
 if(typetest != 0) {
   printerror(application.verbose_state.debug, "pstat: Cannot specify more than one type of statistical test at the time");
   return 0;
 }
 typetest = KSFLAT;
      }else if(strcmp(argv[i], "-ksflat2") == 0) {
 if(typetest != 0) {
   printerror(application.verbose_state.debug, "pstat: Cannot specify more than one type of statistical test at the time");
   return 0;
 }
 int ret;
 ret = parse_command_string(application.verbose_state, argc, argv, i+1, 0, 1, "%lf %lf", &ksflat_min, &ksflat_max, NULL);
 if(ret != 2) {
   printerror(application.verbose_state.debug, "ERROR pstat: Cannot parse '%s' option, 2 values were expected.", argv[i]);
   return 0;
 }
 typetest = KSFLAT_WITH_MINMAX;
 i++;
      }else if(strcmp(argv[i], "-moments") == 0) {
 if(typetest != 0) {
   printerror(application.verbose_state.debug, "pstat: Cannot specify more than one type of statistical test at the time");
   return 0;
 }
 typetest = MOMENTS;
      }else if(strcmp(argv[i], "-median") == 0) {
 if(typetest != 0) {
   printerror(application.verbose_state.debug, "pstat: Cannot specify more than one type of statistical test at the time");
   return 0;
 }
 typetest = MEDIAN;
      }else if(strcmp(argv[i], "-kssin") == 0) {
 if(typetest != 0) {
   printerror(application.verbose_state.debug, "pstat: Cannot specify more than one type of statistical test at the time");
   return 0;
 }
 typetest = KSSIN;
      }else if(strcmp(argv[i], "-chi2hist") == 0) {
 int ret;
 ret = parse_command_string(application.verbose_state, argc, argv, i+1, 0, 1, "%lf %lf %lf", &threshold1, &threshold2, &threshold3, NULL);
 if(ret == 1) {
   threshold2 = threshold3 = -1;
 }else if(ret == 2) {
   threshold3 = -1;
 }else if(ret != 3) {
   printerror(application.verbose_state.debug, "ERROR pstat: Cannot parse '%s' option, 1, 2 or 3 values.", argv[i]);
   return 0;
 }
 if(typetest != 0) {
   printerror(application.verbose_state.debug, "pstat: Cannot specify more than one type of statistical test at the time");
   return 0;
 }
 typetest = CHI2TEST_HIST;
 i++;
      }else if(strcmp(argv[i], "-chi2cdf") == 0) {
 if(typetest != 0) {
   printerror(application.verbose_state.debug, "pstat: Cannot specify more than one type of statistical test at the time");
   return 0;
 }
 typetest = CHI2TEST_CDF;
      }else if(strcmp(argv[i], "-log") == 0) {
 read_log = 1;
      }else {
 if(argv[i][0] == '-') {
   printerror(application.verbose_state.debug, "pstat: Unknown option: %s\n\nRun pstat without command line arguments to show help", argv[i]);
   terminateApplication(&application);
   return 0;
 }else {
   if(applicationAddFilename(i, application.verbose_state) == 0)
     return 0;
 }
      }
    }
  }
  if(applicationFilenameList_checkConsecutive(argv, application.verbose_state) == 0) {
    return 0;
  }
  {
    int nrInputFiles;
    int nrInputColumns;
    nrInputFiles = numberInApplicationFilenameList(&application, argv, application.verbose_state);
    if(nrInputFiles < 1) {
      printerror(application.verbose_state.debug, "ERROR pstat: No files specified");
      return 0;
    }
    if(nrInputFiles > 2) {
      printerror(application.verbose_state.debug, "ERROR pstat: Cannot specify more than two input files. %d input files are currently specified.", nrInputFiles);
      if(application.verbose_state.debug) {
 printerror(0, "  These are:");
 char *filename_ptr;
 for(i = 0; i < nrInputFiles; i++) {
   filename_ptr = getNextFilenameFromList(&application, argv, application.verbose_state);
   if(filename_ptr != NULL) {
     printerror(0, "    %s", filename_ptr);
   }
 }
      }
      return 0;
    }
    if(file1_column1 == 0) {
      {
 printwarning(application.verbose_state.debug, "WARNING pstat: No -col1 option. Using the default (-col1 1).");
 file1_column1 = 1;
 file1_column2 = 0;
 file1_column3 = 0;
      }
    }
    if(nrInputFiles == 2 && file2_column1 == 0) {
      {
 printwarning(application.verbose_state.debug, "WARNING pstat: No -col2 option, while there are two input files. Default is -col2 1.");
 file2_column1 = 1;
 file2_column2 = 0;
 file2_column3 = 0;
      }
    }
    if(file2_column1 && nrInputFiles != 2) {
      printerror(application.verbose_state.debug, "ERROR pstat: -col2 option suggest that there should be two input files, but there is only %d defined on command line.", nrInputFiles);
      return 0;
    }
    nrInputColumns = 0;
    if(file1_column1)
      nrInputColumns++;
    if(file1_column2)
      nrInputColumns++;
    if(file1_column3)
      nrInputColumns++;
    if(file2_column1)
      nrInputColumns++;
    if(file2_column2)
      nrInputColumns++;
    if(file2_column3)
      nrInputColumns++;
    if(typetest == KSTEST) {
      if(nrInputColumns != 2) {
 printerror(application.verbose_state.debug, "ERROR pstat: KS-test requires two columns of data to be read in. Example: pstat -ks -col1 1 -col2 1 file1 file2 or pstat -ks -col \"1 2\" file1");
 return 0;
      }
    }else if(typetest == CORREL) {
      if(nrInputColumns != 2) {
 printerror(application.verbose_state.debug, "ERROR pstat: Correlation requires two columns of data to be read in. Example: pstat -correl -col1 1 -col2 1 file1 file2 or pstat -correl -col \"1 2\" file1");
 return 0;
      }
    }else if(typetest == PEARSON) {
      if(nrInputColumns != 2) {
 printerror(application.verbose_state.debug, "ERROR pstat: Correlation requires two columns of data to be read in. Example: pstat -pearson -col1 1 -col2 1 file1 file2 or pstat -pearson -col \"1 2\" file1");
 return 0;
      }
    }else if(typetest == MOMENTS) {
      if(nrInputColumns != 1) {
 printerror(application.verbose_state.debug, "ERROR pstat: Computation of the moments of a distribution requires one columns of data to be read in. Example: pstat -moments -col 1 file1");
 return 0;
      }
    }else if(typetest == MEDIAN) {
      if(nrInputColumns != 1) {
 printerror(application.verbose_state.debug, "ERROR pstat: Computation of the median of a distribution requires one columns of data to be read in. Example: pstat -median -col 1 file1");
 return 0;
      }
    }else if(typetest == KSFLAT) {
      if(nrInputColumns != 1) {
 printerror(application.verbose_state.debug, "ERROR pstat: KS-test comparison with a flat distribution requires one columns of data to be read in. Example: pstat -ksflat -col 1 file1");
 return 0;
      }
    }else if(typetest == KSFLAT_WITH_MINMAX) {
      if(nrInputColumns != 1) {
 printerror(application.verbose_state.debug, "ERROR pstat: KS-test comparison with a flat distribution requires one columns of data to be read in. Example: pstat -ksflat2 \"0 1\" -col 1 file1");
 return 0;
      }
    }else if(typetest == KSSIN) {
      if(nrInputColumns != 1) {
 printerror(application.verbose_state.debug, "ERROR pstat: KS-test comparison with a sinusoidal distribution requires one columns of data to be read in. Example: pstat -kssin -col 1 file1");
 return 0;
      }
    }else if(typetest == CHI2TEST_HIST) {
      if(nrInputColumns != 4 && nrInputColumns != 6) {
 printerror(application.verbose_state.debug, "ERROR pstat: The chi-square histogram test requires four or six columns of data to be read in. Examples: pstat -chi2hist -1 -col1 \"1 2\" -col2 \"1 2\" file1 file2 or pstat -chi2hist -1 -col1 \"1 2 3\" -col2 \"1 2 3\" file1 file2");
 return 0;
      }
    }else if(typetest == CHI2TEST_CDF) {
      if(nrInputColumns != 2) {
 printerror(application.verbose_state.debug, "ERROR pstat: The chi-square CDF test requires two columns of data to be read in. Examples: pstat -chi2cdf -col1 1 -col2 1 file1 file2 or pstat -chi2cdf -col1 \"1 2\" file1");
 return 0;
      }
    }else {
      printerror(application.verbose_state.debug, "ERROR pstat: No statistical test has been specified, nothing to do.");
      return 0;
    }
  }
  char *filename_ptr;
  double *input_array[6];
  long number_values[6];
  int number_input_arrays;
  number_input_arrays = 0;
  number_values[0] = 0;
  filename_ptr = getNextFilenameFromList(&application, argv, application.verbose_state);
  if(filename_ptr == NULL) {
    printerror(application.verbose_state.debug, "ERROR pstat: Bug!");
    return 0;
  }
  if(file1_column1) {
    double min_x, max_x, avrg;
    if(read_ascii_column_double(filename_ptr, skiplines, '#', -1, 1, &number_values[number_input_arrays], file1_column1, 1.0, read_log, &input_array[number_input_arrays], &min_x, &max_x, &avrg, application.verbose_state, 0) == 0) {
      printerror(application.verbose_state.debug, "ERROR pstat: cannot load file.\n");
      return 0;
    }
    number_input_arrays++;
  }
  if(file1_column2) {
    double min_x, max_x, avrg;
    if(read_ascii_column_double(filename_ptr, skiplines, '#', -1, 1, &number_values[number_input_arrays], file1_column2, 1.0, read_log, &input_array[number_input_arrays], &min_x, &max_x, &avrg, application.verbose_state, 0) == 0) {
      printerror(application.verbose_state.debug, "ERROR pstat: cannot load file.\n");
      return 0;
    }
    number_input_arrays++;
  }
  if(file1_column3) {
    double min_x, max_x, avrg;
    if(read_ascii_column_double(filename_ptr, skiplines, '#', -1, 1, &number_values[number_input_arrays], file1_column3, 1.0, read_log, &input_array[number_input_arrays], &min_x, &max_x, &avrg, application.verbose_state, 0) == 0) {
      printerror(application.verbose_state.debug, "ERROR pstat: cannot load file.\n");
      return 0;
    }
    number_input_arrays++;
  }
  if(file2_column1 || file2_column2 || file2_column3) {
    filename_ptr = getNextFilenameFromList(&application, argv, application.verbose_state);
    if(filename_ptr == NULL) {
      printerror(application.verbose_state.debug, "ERROR pstat: Bug!");
      return 0;
    }
  }
  if(file2_column1) {
    double min_x, max_x, avrg;
    if(read_ascii_column_double(filename_ptr, skiplines, '#', -1, 1, &number_values[number_input_arrays], file2_column1, 1.0, read_log, &input_array[number_input_arrays], &min_x, &max_x, &avrg, application.verbose_state, 0) == 0) {
      printerror(application.verbose_state.debug, "ERROR pstat: cannot load file.\n");
      return 0;
    }
    number_input_arrays++;
  }
  if(file2_column2) {
    double min_x, max_x, avrg;
    if(read_ascii_column_double(filename_ptr, skiplines, '#', -1, 1, &number_values[number_input_arrays], file2_column2, 1.0, read_log, &input_array[number_input_arrays], &min_x, &max_x, &avrg, application.verbose_state, 0) == 0) {
      printerror(application.verbose_state.debug, "ERROR pstat: cannot load file.\n");
      return 0;
    }
    number_input_arrays++;
  }
  if(file2_column3) {
    double min_x, max_x, avrg;
    if(read_ascii_column_double(filename_ptr, skiplines, '#', -1, 1, &number_values[number_input_arrays], file2_column3, 1.0, read_log, &input_array[number_input_arrays], &min_x, &max_x, &avrg, application.verbose_state, 0) == 0) {
      printerror(application.verbose_state.debug, "ERROR pstat: cannot load file.\n");
      return 0;
    }
    number_input_arrays++;
  }
  FILE *fout;
  if(output_idx) {
    fout = fopen(argv[output_idx], "w");
    if(fout == NULL) {
      printerror(application.verbose_state.debug, "ERROR pstat: Cannot open %s", argv[output_idx]);
      return 0;
    }
  }else {
    fout = stdout;
  }
  if(typetest == KSTEST) {
    double max_diff, prob;
    if(number_input_arrays != 2) {
      printerror(application.verbose_state.debug, "ERROR pstat: KS-test requires two input arrays to be specified");
      return 0;
    }
    kstest(input_array[0], number_values[0], input_array[1], number_values[1], 0, 0, 0, NULL, &max_diff, &prob, application.verbose_state);
    if(application.verbose_state.verbose == 0) {
      fprintf(fout, "probability = %e\n", prob);
    }
  }else if(typetest == KSFLAT) {
    double max_diff, prob;
    if(number_input_arrays != 1) {
      printerror(application.verbose_state.debug, "ERROR pstat: KS-test comparison with a flat distribution requires one input array to be specified");
      return 0;
    }
    kstest(input_array[0], number_values[0], NULL, 0, 1, 0, 0, NULL, &max_diff, &prob, application.verbose_state);
    if(application.verbose_state.verbose == 0) {
      fprintf(fout, "probability = %e\n", prob);
    }
  }else if(typetest == KSFLAT_WITH_MINMAX) {
    double max_diff, prob;
    if(number_input_arrays != 1) {
      printerror(application.verbose_state.debug, "ERROR pstat: KS-test comparison with a flat distribution requires one input array to be specified");
      return 0;
    }
    kstest(input_array[0], number_values[0], NULL, 0, 3, ksflat_min, ksflat_max, NULL, &max_diff, &prob, application.verbose_state);
    if(application.verbose_state.verbose == 0) {
      fprintf(fout, "probability = %e\n", prob);
    }
  }else if(typetest == KSSIN) {
    double max_diff, prob;
    if(number_input_arrays != 1) {
      printerror(application.verbose_state.debug, "ERROR pstat: KS-test comparison with a sinusoidal distribution requires one input array to be specified");
      return 0;
    }
    kstest(input_array[0], number_values[0], NULL, 0, 2, 0, 0, NULL, &max_diff, &prob, application.verbose_state);
    if(application.verbose_state.verbose == 0) {
      fprintf(fout, "probability = %e\n", prob);
    }
  }else if(typetest == MOMENTS) {
    double mean, variance, skew, kurt;
    if(number_input_arrays != 1) {
      printerror(application.verbose_state.debug, "ERROR pstat: Computation of the moments of a distribution requires one columns of data to be read in.");
      return 0;
    }
    mean = gsl_stats_mean(input_array[0], 1, number_values[0]);
    fprintf(fout, "Mean               = %e\n", mean);
    variance = gsl_stats_variance_m(input_array[0], 1, number_values[0], mean);
    fprintf(fout, "Variance           = %e (unbiased, normalised by 1/(N-1))\n", variance);
    fprintf(fout, "Standard deviation = %e (unbiased, normalised by 1/(N-1))\n", sqrt(variance));
    variance *= (number_values[0]-1.0)/(double)number_values[0];
    fprintf(fout, "Variance           = %e (normalised by 1/N)\n", variance);
    fprintf(fout, "Standard deviation = %e (normalised by 1/N)\n", sqrt(variance));
    skew = gsl_stats_skew(input_array[0], 1, number_values[0]);
    fprintf(fout, "Skewness           = %e\n", skew);
    kurt = gsl_stats_kurtosis(input_array[0], 1, number_values[0]);
    fprintf(fout, "Kurtosis           = %e\n", kurt);
  }else if(typetest == MEDIAN) {
    double median;
    if(number_input_arrays != 1) {
      printerror(application.verbose_state.debug, "ERROR pstat: Computation of the moments of a distribution requires one columns of data to be read in.");
      return 0;
    }
    gsl_sort(input_array[0], 1, number_values[0]);
    median = gsl_stats_median_from_sorted_data(input_array[0], 1, number_values[0]);
    fprintf(fout, "Median = %e\n", median);
  }else if(typetest == PEARSON) {
    double cc;
    if(number_input_arrays != 2) {
      printerror(application.verbose_state.debug, "ERROR pstat: Correlation requires two input arrays to be specified");
      return 0;
    }
    if(number_values[0] != number_values[1]) {
      printerror(application.verbose_state.debug, "ERROR pstat: Correlation requires two input arrays of equal length");
      return 0;
    }
#if GSL_VERSION_NUMBER >= 110
    cc = gsl_stats_correlation(input_array[0], 1, input_array[1], 1, number_values[0]);
    fprintf(fout, "Pearson correlation coefficient = %e\n", cc);
#else
    printerror(application.verbose_state.debug, "ERROR pstat: Pearson correlation coefficient cannot be calculated if GSL < 1.10");
    return 0;
#endif
  }else if(typetest == CORREL) {
    if(number_input_arrays != 2) {
      printerror(application.verbose_state.debug, "ERROR pstat: Correlation requires two input arrays to be specified");
      return 0;
    }
    if(number_values[0] != number_values[1]) {
      printerror(application.verbose_state.debug, "ERROR pstat: Correlation requires two input arrays of equal length");
      return 0;
    }
    if(number_values[0] > 2147483646) {
      printerror(application.verbose_state.debug, "ERROR pstat: Input array too long.");
      return 0;
    }
    float *x1, *y1;
    x1 = malloc(number_values[0]*sizeof(float));
    y1 = malloc(number_values[1]*sizeof(float));
    if(x1 == NULL || y1 == NULL) {
      printerror(application.verbose_state.debug, "ERROR pstat: Memory allocation error");
      return 0;
    }
    for(i = 0; i < number_values[0]; i++)
      x1[i] = (input_array[0])[i];
    for(i = 0; i < number_values[1]; i++)
      y1[i] = (input_array[1])[i];
    int extrazeropad, cc_length;
    float *ans;
    extrazeropad = 0;
    if(crosscorrelation_fft_padding(x1, y1, number_values[0], extrazeropad, &ans, &cc_length, application.verbose_state) == 0) {
      printerror(application.verbose_state.debug, "ERROR pstat: Cross correlation failed.");
      return 0;
    }
    free(x1);
    free(y1);
    int lagoutput;
    lagoutput = 1;
    if(lagoutput == 0) {
      for(i = 0; i < cc_length; i++)
 fprintf(fout, "%ld %e\n", i, ans[i]);
    }else {
      for(i = cc_length/2; i <= cc_length-1; i++)
 fprintf(fout, "%ld %e\n", i-cc_length, ans[i]);
      for(i = 0; i < cc_length/2; i++)
 fprintf(fout, "%ld %e\n", i, ans[i]);
    }
    free(ans);
  }else if(typetest == CHI2TEST_HIST) {
    double binwidth, ratio, offset, chi2;
    long offset_binnr, file2_x_col, dof, i2, nr_overlapping_bins;
    if(number_input_arrays != 4 && number_input_arrays != 6) {
      printerror(application.verbose_state.debug, "ERROR pstat: The chi-square histogram test requires four or six columns of data to be specified.");
      return 0;
    }
    if(number_input_arrays == 4) {
      file2_x_col = 2;
      printwarning(application.verbose_state.debug, "WARNING pstat: Since no column numbers with error-bars are provided, uniform weighting of the different bins is assumed with sigma=1. This is unlikely to be correct.");
    }else {
      file2_x_col = 3;
    }
    binwidth = (input_array[0])[1] - (input_array[0])[0];
    ratio = binwidth/((input_array[file2_x_col])[1] - (input_array[file2_x_col])[0]);
    if(application.verbose_state.verbose) {
      printf("Ratio bin widths of two histograms is %lf (should be very close to 1)\n", binwidth);
    }
    if(ratio < 0.999 || ratio > 1.001) {
      printerror(application.verbose_state.debug, "ERROR pstat: The binwidths of the two histograms appear to be different (%e != %e).", (input_array[0])[1] - (input_array[0])[0], (input_array[file2_x_col])[1] - (input_array[file2_x_col])[0]);
      return 0;
    }
    offset = (input_array[file2_x_col])[0];
    offset -= (input_array[0])[0];
    offset /= binwidth;
    if(application.verbose_state.verbose) {
      printf("Offset between two histograms is %lf bins (should be very close to an integer value)\n", offset);
    }
    offset_binnr = round(offset);
    offset = fabs(offset - offset_binnr);
    if(offset > 0.001) {
      printerror(application.verbose_state.debug, "ERROR pstat: The bins of the two histograms do not appear to be aligned, but have an offset of %lf.", offset);
      return 0;
    }
    double height1, height2;
    height1 = 0;
    height2 = 0;
    for(i = 0; i < number_values[0]; i++)
      height1 += (input_array[1])[i];
    for(i = 0; i < number_values[file2_x_col]; i++)
      height2 += (input_array[file2_x_col+1])[i];
    if(application.verbose_state.verbose) {
      printf("Ratio of integrals of two histograms is %lf (should be very close to 1)\n", height1/height2);
    }
    if(height1/height2 > 1.001 || height2/height1 > 1.001) {
      printwarning(application.verbose_state.debug, "WARNING pstat: The two histograms appear to be normalised differently, so the derived numbers are unlikely to give useful results.");
    }
    chi2 = 0;
    dof = 0;
    nr_overlapping_bins = 0;
    for(i = -labs(offset_binnr)-10; i < number_values[0]+labs(offset_binnr)+10; i++) {
      int hist1_exist, hist2_exist;
      hist1_exist = hist2_exist = 0;
      if(i >= 0 && i < number_values[0]) {
 height1 = (input_array[1])[i];
 hist1_exist = 1;
      }else {
 height1 = 0;
      }
      i2 = i - offset_binnr;
      if(i2 >= 0 && i2 < number_values[file2_x_col]) {
 height2 = (input_array[file2_x_col+1])[i2];
 hist2_exist = 1;
 if(hist1_exist) {
   offset = (input_array[file2_x_col])[i2];
   offset -= (input_array[0])[i];
   offset /= binwidth;
   if(fabs(offset) > 0.001) {
     printerror(application.verbose_state.debug, "ERROR pstat: The bins of the two histograms do not appear to be aligned, but bins %ld and %ld have an offset of %lf.", i+1, i2+1, offset);
     return 0;
   }
   nr_overlapping_bins++;
 }
      }else {
 height2 = 0;
      }
      if(hist1_exist || hist2_exist) {
 if(height1 >= threshold1 && height2 >= threshold2 && (height1+height2) >= threshold3) {
   double delta_y, var;
   delta_y = height2 - height1;
   if(number_input_arrays == 4) {
     chi2 += delta_y*delta_y;
     dof++;
   }else {
     var = 0;
     if(hist1_exist)
       var += (input_array[2])[i] * (input_array[2])[i];
     if(hist2_exist)
       var += (input_array[file2_x_col+2])[i2] * (input_array[file2_x_col+2])[i2];
     chi2 += delta_y*delta_y/var;
     dof++;
   }
 }
      }
    }
    if(nr_overlapping_bins == 0) {
      printerror(application.verbose_state.debug, "ERROR pstat: There appears to be no overlap between the two histograms.");
      return 0;
    }
    if(application.verbose_state.verbose) {
      printf("Number of overlapping bins between two distributions:     %ld\n", nr_overlapping_bins);
      printf("Total number of bins considered with specified threshold: %ld\n", dof);
    }
    if(number_input_arrays != 6)
      fprintf(fout, "Total non-weighted chi square:   %f = %e\n", chi2, chi2);
    else if(application.verbose_state.verbose)
      fprintf(fout, "Total chi square:   %f = %e\n", chi2, chi2);
    if(number_input_arrays == 6) {
      fprintf(fout, "Reduced chi square: %f = %e\n", chi2/(double)dof, chi2/(double)dof);
    }
  }else if(typetest == CHI2TEST_CDF) {
    double x1, cdf1, cdf2, chi2;
    long start2, dof;
    if(number_input_arrays != 2) {
      printerror(application.verbose_state.debug, "ERROR pstat: The chi-square CDF test requires two columns of data to be specified.");
      return 0;
    }
    if(number_values[0] > number_values[1]) {
      i = number_values[0];
      number_values[0] = number_values[1];
      number_values[1] = i;
      input_array[2] = input_array[0];
      input_array[0] = input_array[1];
      input_array[1] = input_array[2];
    }
    gsl_sort(input_array[0], 1, number_values[0]);
    gsl_sort(input_array[1], 1, number_values[1]);
    start2 = 0;
    chi2 = 0;
    dof = 0;
    for(i = 0; i < number_values[0] - 1; i++) {
      cdf1 = (i+1)/(double)number_values[0];
      x1 = (input_array[0])[i] + 0.5*((input_array[0])[i+1]-(input_array[0])[i]);
      if(application.verbose_state.debug)
 printf("Going to find chi2 of observation distribution cdf point: (%e, %e)\n", x1, cdf1);
      if(x1 <= (input_array[1])[0]) {
 if(application.verbose_state.debug)
   printf("  model cdf starts after point of interest\n");
 cdf2 = 0;
 if(x1 == (input_array[1])[0]) {
   cdf2 = (j+0.5)/(double)number_values[1];
 }
      }else if(x1 >= (input_array[1])[number_values[1]-1]) {
 if(application.verbose_state.debug)
   printf("  model cdf ends before point of interest\n");
 cdf2 = 1;
 if(x1 == (input_array[1])[number_values[1]-1]) {
   cdf2 = (number_values[1]-1 +0.5)/(double)number_values[1];
 }
      }else {
 for(j = start2; j < number_values[1]; j++) {
   if((input_array[1])[j] >= x1) {
     if(application.verbose_state.debug)
       printf("  First model cdf point after point of interest: (%e %e)\n", (input_array[1])[j], (j+1)/(double)number_values[1]);
     long index1, index2;
     index1 = j-1;
     index2 = j;
     if(j == 0) {
       printerror(application.verbose_state.debug, "ERROR pstat: Bug!\n");
       return 0;
     }
     while((input_array[1])[index2] == (input_array[1])[index1]) {
       if(index2 < number_values[1] - 1) {
  index2++;
       }else if(index1 > 1) {
  index1--;
       }else {
  printerror(application.verbose_state.debug, "ERROR pstat: Something is wrong with the second input distribution, as all input values appear to be identical.\n");
  return 0;
       }
     }
     if(application.verbose_state.debug)
       printf("  Going to interpolate following model points: (%e %e) and (%e %e)\n", (input_array[1])[index1], (index1+1-0.5)/(double)number_values[1], (input_array[1])[index2], (index2+1-0.5)/(double)number_values[1]);
     if(index1 < 0 || index2 >= number_values[1]) {
       printerror(application.verbose_state.debug, "ERROR pstat: Interpolation failed - trying to access second array out of limits (index1=%ld, index2=%ld).\n", index1, index2);
       return 0;
     }
     cdf2 = ((index2+1-0.5)/(double)number_values[1] - (index1+1-0.5)/(double)number_values[1]) * (x1 - (input_array[1])[index1]) / ((input_array[1])[index2] - (input_array[1])[index1]) + (index1+1-0.5)/(double)number_values[1];
     if(!isfinite(cdf2)) {
       printerror(application.verbose_state.debug, "ERROR pstat: Interpolation failed: i=%ld, j=%ld (array2=[%e, %e, %e]), \n", i, j, (input_array[1])[j-1], (input_array[1])[j], (input_array[1])[j+1]);
       return 0;
     }
     start2 = j - 2;
     if(start2 < 0)
       start2 = 0;
     break;
   }
 }
      }
      chi2 += (cdf2-cdf1)*(cdf2-cdf1);
      dof++;
      if(application.verbose_state.debug) {
 printf("  cdf2=%e, cdf1=%e, diff=%e\n", cdf2, cdf1, (cdf2-cdf1));
 printf("  new chi2 = %e\n", chi2);
      }
    }
    if(application.verbose_state.verbose) {
      printf("\nTotal number of bins considered: %ld\n", dof);
    }
    fprintf(fout, "Non-weighted total chi square = %f = %e\n", chi2, chi2);
  }else {
    printerror(application.verbose_state.debug, "ERROR pstat: Bug!");
    return 0;
  }
  if(output_idx) {
    fclose(fout);
  }
  for(i = 0; i < number_input_arrays; i++) {
    free(input_array[i]);
  }
  terminateApplication(&application);
  return 0;
}

max_stars_repo_name: David-McKenna/psrsalsa

max_issues_repo_name: David-McKenna/psrsalsa

max_forks_repo_name: David-McKenna/psrsalsa


text:
/*************************************************************************

description:
This file contains the class that can do the actual radiative transport,
including routines for grid calculation, physics and the actual transport

Copyright Jan-Pieter Paardekooper and Chael Kruip October 2011

This file is part of SimpleX.

SimpleX 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.

SimpleX 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 SimpleX.  If not, see <http://www.gnu.org/licenses/>.

**************************************************************************/

#ifndef SIMPLEX_H
#define SIMPLEX_H

#ifndef NOMPI
#include "mpi.h"
#endif

#include "rates.h"
#include "Common.h"
#include "Structs.h"
#include <gsl/gsl_rng.h>  //random number generator
#include "configfile.h"   //keyvalue input file

#ifdef HEAL_PIX
  #include "healpix_base.h" //healpix header
#endif

#include "tree_structures.h" //octree
#include "hilbert.h"         //hilbert curve

#include "h5w_serial.h"      //hdf5 header

#include <algorithm>

#if defined(__cplusplus)
extern "C"
{
#endif
#include <stdio.h>
#include <stdlib.h>
#include <libqhull.h>
#include <mem.h>
#include <qset.h>
#include <geom.h>
#include <merge.h>
#include <poly.h>
#include <io.h>
#include <stat.h>
#if defined(__cplusplus)
}
#endif

using namespace std;

#define VERTEX_ITERATOR vector< Vertex >::iterator
#define SITE_ITERATOR vector< Site >::iterator
#define SIMPL_ITERATOR vector< Simpl >::iterator

//! Class in which all functions that work on the SimpleX grid are stored

//! This is the main class of the simulation!
class SimpleX{

  public:
    //! constructor
    SimpleX(const string& output_path = ".", const string& data_path = ".");
    //! destructor
    ~SimpleX();

  //================ Grid Initilisation Functions ==================================//

    //! Create triangulation

    //! Initialize the simulation by creating/reading the vertex list and compute the triangulation
    //! Result is a list of Sites and a list of Simplices
    void init_triangulation(char* inputName);

    //TEST
    void reinit_triangulation();
    //END TEST

    //! Read in parameter file
    void read_parameters( char* inputName );

    //! Create orientations and mappings from correct header file
    void set_direction_bins();

    //! Create a simple homogeneous point distribution

    //! Points are placed using the gsl randoom number generator
    void poisson_square();

    //! Read in hdf5 file with site information
    void read_vertex_list();

    //! Compute the largest index of all the vertices
    unsigned long long int computeMaxId();

    //! Create boundary around computational domain

    //! Width of the boundary is variable borderBox.
    //! Number of points in the boundary is variable borderSites.
    void create_boundary();

    //!Create octree containing all vertices

    //!Size of the tree depends on the number of subboxes used, which
    //!is determined by the hilbert order m
    void create_vertex_tree();

    //! Decompose the domain

    //! Count the number of sites contained in one subbox and add those until
    //! number of points is approximately equal per processor
    void decompose_domain();

    //! Assign the correct process to vertices
    void assign_process();

    //! Check if site is in unit domain
    bool inDomain(const float& x, const float& y, const float& z, const unsigned int& rank);

    //! Compute the triangulation

    //! Call to QHull to perform the triangulation of vertices in one subbox
    //! The result is a vector of simplices
    //! Also, a check is done whether the subbox and domain boundaries are
    //! sufficiently large. If not, subbox boundaries are extended. If the
    //! boundary around the unity domain is too small, the code exits,
    //! to avoid memory issues when this number would be increased
    void compute_triangulation();

    //! Create the sites array on which radiative transfer is performed

    //! From the list of simplices the vertices relevant for this proc are selected and
    //! put in the sites vector
    void create_sites();

    //! Give each site an id corresponding to place in local sites vector
    void assign_site_ids();

    //! Determine which sites use ballistic transport and which sites
    //! should use direction conserving transport at the start of the
    //! simulation.
    void initiate_ballistic_sites();

    //! Shuffle list of sites

    //! This is necessary since the sites list is created from the simplices returned by QHull,
    //! so they are ordered.
    void shuffle_sites();

  //================ Site Properties Functions ==================================//

    //! Compute properties of the sites from the triangulation

    //! Important properties like the neighbours, the volumes and
    //! the most straight paths are computed
    void compute_site_properties();

    //! Compute the neighbours of vertex

    //! Generate the neighbour vectors for each vertex from the simplex vector
    //! Care should be taken so that no vertex is put twice in the neighbour list.
    //! Also get rid of vertices in the sub box boundary that are not connected to a site
    //! in the computational domain
    void compute_neighbours();

    //! Compute the volume of the Voronoi cell

    //! Use the volumes of the simplices that the vertex is part of to compute
    //! the volume of the Voronoi cell around the vertex
    void compute_volumes();

    //! Calculate the most straight paths for every incoming direction

    //! Needed for fast radiative transport
    void calculate_straight();

    //! Identify the place of every ballistic site in neighbour array
    //! of neighbouring site
    void match_neighbours();

    //! Compute the solid angles into which the intensity is distributed
    void compute_solid_angles( const bool& rotate );

    //! Check for cyclic connections in the grid and remove them if possible
    void check_cyclic_connections();

    //! Calculate mean distance to neighbours
    void calculate_line_lengths();

    //! Create list without border simplicesc
    void remove_border_simplices();

    //================ Grid Update Routines ============================//

    //! Rotate the solid angles of the unit sphere
    //! to avoid preferential directions
    void rotate_solid_angles();

    //! calculate normal distributed random numbers using the Box-Muller transformation
    void randGauss( float& y1, float& y2 );

    //! Store the intensities for use in later run, directions are kept.
    void store_intensities();

    //! Get the ballistic sites that have intensity to be stored
    const vector< unsigned long int > get_ballistic_sites_to_store();

    //! Store ballistic intensities in unit sphere tesselation
    void store_ballistic_intensities( const vector< unsigned long int >& sites_to_store );

    //! Return the intensities from the previous run
    void return_intensities();

    //! Return the intensities to ballistic sites
    void return_ballistic_intensities();

    //! Check which sites are ballistic and which not
    //! according to user specified switch

    //! Only used is case of combined transport
    void check_ballistic_sites();

    //! Create new sites array without the removed sites
    void store_site_properties();

    //! Create updated vertex list from site_properties vector
    void create_new_vertex_list();

    //================ Send Routines ===================================//

    //! Send vertices to processors
    void send_vertices();

    //! Send the domain decomposition to all procs
    void send_dom_dec();

    //! Fill the list of indices of sites that contain information to be send
    //! to other procs
    void fill_send_list();

    //! Give every vertex its correct process (needed for vertices in boundary between procs)
    void send_site_properties();

    //! Give every ballistic site its correct neighbour information

    //! This is only needed for ballistic sites in boundary between procs, and is necessary because
    //! the order of the neighbours in the neighbour vector is not the same among procs.
    void send_neighbour_properties();

    //! Send teh updated ballistic sites to all procs
    void send_site_ballistics( const vector< unsigned long int >& sites_to_send );

     //! Send the intensities among procs
    void send_intensities();

    //! Send site physics to all procs

    //! In case sites change from proc send physics
    void send_site_physics();

    //! In case sites change from proc send intensities
    void send_site_intensities();

    //! Send stored intensities among procs
    void send_stored_intensities();

    //! Send the updated vertex list to all procs
    void send_new_vertex_list();

    //================ Physics Functions ==================================//

    //! Compute all the physics needed to do radiative transport
    void compute_physics( const unsigned int& run );

    //! Initialize physical parameters
    void initialise_physics_params();

    //! Read metal line cooling data and initialize vectors
    void read_metals();

    //! set homogeneous number density.

    //! Use only in combination with Poisson_Square()
    void set_homogeneous_number_density( const float& nH );

    //! Give the masses and fluxes that were read in to the sites

    //! Use only in combination with read_vertex_list()
    void assign_read_properties();

    //! Give the updates sites their correct physical quantities from stored values
    void return_physics();

    //! Calculate the bounds of the frequency bins
    vector<double> calc_freq_bounds( const double& lowerBound, const double& upperBound );

    //! calculate effective cross section and spectrum of black body source
    void black_body_source( const double& tempSource );

    //! Compute the total number of atoms on this proc, used for output
    void output_number_of_atoms();

    //! Return mean optical depth of total grid
    double output_optical_depth();

    //! output relevant physical parameters to screen
    void parameter_output( const unsigned int& run );

    //! Add more neighbours to vertices that have a flux, by adding the neighbours
    //! of neighbours.
    void make_source_isotropic();

    //! calculate recombination
    double recombine();

    //! Rate of energy loss in cell
    double cooling_rate( Site& site );

    //! Rate of energy gain in cell
    double heating_rate( const vector<double>& N_ion, const double& t_end );

    //! compute mean molecular weight
    double compute_mu( Site& );

    //!Convert internal energy to temperature
    double u_to_T( const double& u, const double& mu );

    //!Convert temperature to internal energy
    double T_to_u( const double& T, const double& mu );

    //! update the temperature by including heating and cooling processes
    double update_temperature( Site& site, const vector<double>& N_ion, const double& t_end );

 //================= Radiative Transfer Functions ======================//

    //! Solve the rate equation to determine ionisations and recombinations
    vector<double> solve_rate_equation( Site& site );

    //! Send source photons
    void source_transport( Site& site );

    //! Transport photons diffusely
    void diffuse_transport( Site& site, vector<double>& N_out_total );

    //! Redistribute intensity: actual radiative transport
    void non_diffuse_transport( Site& site, vector<double>& N_out_total );

    //! Call to do radiative transfer
    void radiation_transport( const unsigned int& run );

  //================== Output Functions ================================//

    //! Call the output routines
    void generate_output( const unsigned int& run );

    //! Write hdf5 output file
    void write_hdf5_output(char* name, const unsigned int& run);

    //! calculate position of the I-Front for source in centre
    double calc_IFront( const unsigned int& run );

    //! calculate escape fraction
    double calc_escape_fraction( const unsigned int& run, const unsigned int& times );

    //================ Generic Functions ==============================//

    //! clear all temporary structures for the next run
    void clear_temporary();

    //! Function to give diagnostic output about total number of photons in simulation
    double count_photons();

    //! return number of runs
    const unsigned int& get_numRuns() const{ return numRuns; }

    //! Return number of outputs
    const unsigned int& get_numOutputs() const{ return numOutputs; }

    // ================ Public lists and variables ===================//

    //! List of vertices that will be used to perform triangulation
    vector< Vertex > vertices;

    //! List of sites for doing radiative transfer
    vector< Site > sites;

    //! Vector holding the domain decomposition
    vector< unsigned int > dom_dec;

    //! List of sites from the previous run, whose properties are still needed
    vector< Site_Update > site_properties;

    //! List of entries in site_intensities array of every site
    vector<unsigned long int> intens_ids;
    //! List of site intensities from previous run, that are still needed
    vector<float> site_intensities;

    //! List of indices of sites needed for transport along procs
    vector< unsigned long int > send_list;

    //! List of site indices that might be removed
    vector< unsigned long int > sites_marked_for_removal;
    //! List of sites and properties to be updated
    vector< Site_Remove > total_sites_marked_for_removal;

    //! List of simplices that results from the QHull call
    vector< Simpl > simplices;
    //! Temporary list of neutral number densities that are read in
    vector< float > temp_n_HI_list;
    //! Temporary list of ionised number densities that are read in
    vector< float > temp_n_HII_list;
    //! Temporary list of fluxes that are read in
    vector< float > temp_flux_list;
    //! Temporary list of internalEnergies that are read in
    vector< float > temp_u_list;
    //! Temporary list of dudt that are read in
    vector< float > temp_dudt_list;
    //! Temporary list of clumping factors that are read in
    vector< float > temp_clumping_list;
    //! Temporary list of metallicity
    vector<float> temp_metallicity_list;

    //! Cross section of hydrogen
    vector<double> cross_H;
    //! The photon excess energy used to heat the gas
    vector<double> photon_excess_energy;
    //! The elemental abundances of metals (and Helium)
    vector<double> abundances;
    //! Vector of cooling curves
    vector< cooling_curve > curves;

    //! Variable used by random number generator
    gsl_rng * ran;
    //! Output to logfile
    ofstream simpleXlog;

  protected:

    unsigned long int numSites;            //!< Total number of sites in the simulation
    unsigned long int origNumSites;        //!< Original number of sites
    unsigned int borderSites;      //!< Number of sites in the boundary around computational domain
    unsigned int hilbert_order;    //!< Hilbert order that determines the number of subboxes
    unsigned int numSweeps;        //!< Number of grid sweeps during simulation
    unsigned int numRuns;          //!< Number of runs of the simulation
    unsigned int numOutputs;       //!< Number of outputs of the simulation
    unsigned int numPixels;        //!< Number of pixels in HEALPix calculation
    unsigned int numStraight;      //!< Number of straight directions
    int randomSeed;                //!< Seed for random number generator
    short dimension;               //!< Dimension of the simulation (default = 3)
    string inputFileName;          //!< File which holds densities and fluxes to be read in
    string dataPath;               //!< Path where the metal line cooling tables are
    bool blackBody;                //!< Use blackbody for source spectrum?
    bool recombination;            //!< Include recombination or not
    bool rec_rad;                  //!< Include recombination radiation or not
    bool coll_ion;                 //!< Include collisional ionisations?
    bool heat_cool;                //!< Include heating and cooling?
    bool metal_cooling;            //!< Include metal line cooling?
    short int numFreq;             //!< Number of frequencies
    short int freq_spacing;        //!< Spacing of the frequency bins
    short RTmethod;                //!< Method for RT
    bool diffuseTransport;         //!< Set this to 1 to do only ballistic transport
    bool ballisticTransport;       //!< Set this to 1 to do only ballistic transport
    bool dirConsTransport;         //!< Set this to 1 to do only direction conserving transport
    bool combinedTransport;        //!< Set this to 1 to do combined transport
    bool photon_conservation;      //!< Use temporal photon conservation scheme or not?
    double subcycle_frac;          //!< Fraction of the characteristic time scale at which subcycling is done
    double sizeBox;                //!< Size of the box in parsec
    float borderBox;               //!< Size of boundary around computational domain
    float padding_subbox;          //!< Subbox boundary in which extra points are triangulated
    double simTime;                //!< Total simulation time in Myr
    double time_conversion;        //!< Time conversion for output
    double sourceTeff;             //!< Effective temperature of the sources in K
    double gasIonTemp;             //!< Temperature of the ionised gas in K
    double gasNeutralTemp;         //!< Temperature of the neutral gas in K
    double switchTau;              //!< Optical depth below which ballistic transport is no longer trusted
    double totalVolume;            //!< Total volume of all Voronoi cells added up
    double totalAtoms;             //!< Total number of atoms in the simulation domain
    int localMaxRes;               //!< Maximum resolution on the local processor
    unsigned int COMM_RANK;        //!< Rank of this proc
    unsigned int COMM_SIZE;        //!< Total number of procs
    unsigned int chunk_size;       //!< Maximum size of chunks written in one piece to hdf5 file
    unsigned int max_msg_to_send;  //!< Maximum number of messages to send in MPI call
    unsigned long int vertex_id_max;    //!< Maximum vertex index
    short int orientation_index;      //!< Orientation of the unit sphere tesselation
    short int orientation_index_old;  //!< Orientation of the unit sphere tesselation in previous run
    float euler_phi;               //!< Euler angle phi for random rotation of coordinate system
    float euler_theta;             //!< Euler angle theta for random rotation of coordinate system
    float euler_psi;               //!< Euler angle psi for random rotation of coordinate system
    bool heal_pix;                 //!< Use heal_pix for associating directions and neighbours or not?
    int num_ref_pix_HP;            //!< Number of pixels to use for associating neighbours in compute_solid_angles()
    float straightAngle;           //!< Maximum angle allowed between straight direction and Delaunay direction.
    bool straight_from_tess;       //!< In DCT, calculate straightest neighbours wrt Delaunay edge or direction bin?
    double UNIT_L;                 //!< Unit length in box
    double UNIT_M;                 //!< Unit mass in box, converts number density into number of atoms
    double UNIT_I;                 //!< Unit number of ionising photons in box
    double UNIT_T;                 //!< Time of 1 sweep in seconds
    double UNIT_T_MYR;             //!< Time of 1 sweep in Myr
    double UNIT_D;                 //!< Unit number density
    double UNIT_V;                 //!< Unit volume
    float straight_correction_factor; //!< Correction factor to correct for the fact that no straight lines exist on the grid
    vertex_tree vert_tree;         //!< Octree that will contain the vertices
    bool give_IFront;                   //!< Calculate IFront position for a source in centre?
    ofstream IFront_output;        //!< Output of IFront data to file
    float*** orient;               //!< Array to hold the direction bins
    unsigned int*** maps;          //!< Array to hold the mappings between direction bins
    unsigned int number_of_directions; //!< Number of discretizations of the unit sphere for DCT
    unsigned int number_of_orientations;//!< Number of random rotations of the DCT discretization
};

#endif

max_stars_repo_name: joshuawall/amuse

max_issues_repo_name: joshuawall/amuse

max_forks_repo_name: joshuawall/amuse


text:
#include <stdio.h>
#include <gsl/gsl_blas.h>

int main()
{
  double a[] = { 0.11, 0.12, 0.13,
                0.21, 0.22, 0.23 };

  double b[] = { 1011, 1012,
                 1021, 1022,
                 1031, 1032 };

  double c[] = { 0.00, 0.00,
                 0.00, 0.00 };

  gsl_matrix_view A = gsl_matrix_view_array(a, 2, 3);
  gsl_matrix_view B = gsl_matrix_view_array(b, 3, 2);
  gsl_matrix_view C = gsl_matrix_view_array(c, 2, 2);

  gsl_matrix *D = gsl_matrix_alloc(2, 3);
  int i, j;
  for (i=1; i<3; i++)
    for (j=1; j<4; j++)
      gsl_matrix_set (D, i-1, j-1, ((float) i)/10 + ((float) j)/100);

  for (i=0; i<2; i++)
    for (j=0; j<3; j++)
      printf("D(%d,%d) = %f\n", i, j, gsl_matrix_get(D, i, j));

  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
                 1.0, D, &B.matrix,
                 0.0, &C.matrix);

  printf ("[ %g, %g\n", c[0], c[1]);
  printf ("  %g, %g ]\n", c[2], c[3]);

  return 0;
}

max_stars_repo_name: tcrundall/chronostar

max_issues_repo_name: tcrundall/chronostar

max_forks_repo_name: tcrundall/chronostar