text:
/*
 #  File            : condat_l1ballproj.c
 #
 #  Version History : 1.0, Aug. 15, 2014
 #
 #  Author          : Laurent Condat, PhD, CNRS research fellow in France.
 #
 #  Description     : This file contains an implementation in the C language
 #                    of algorithms described in the research paper:
 #
 #                    L. Condat, "Fast Projection onto the Simplex and the
 #                    l1 Ball", preprint Hal-01056171, 2014.
 #
 #                    This implementation comes with no warranty: due to the
 #                    limited number of tests performed, there may remain
 #                    bugs. In case the functions would not do what they are
 #                    supposed to do, please email the author (contact info
 #                    to be found on the web).
 #
 #                    If you use this code or parts of it for any purpose,
 #                    the author asks you to cite the paper above or, in
 #                    that event, its published version. Please email him if
 #                    the proposed algorithms were useful for one of your
 #                    projects, or for any comment or suggestion.
 #
 #  Usage rights    : Copyright Laurent Condat.
 #                    This file is distributed under the terms of the CeCILL
 #                    licence (compatible with the GNU GPL), which can be
 #                    found at the URL "http://www.cecill.info".
 #
 #  This software is governed by the CeCILL license under French law and
 #  abiding by the rules of distribution of free software. You can  use,
 #  modify and or redistribute the software under the terms of the CeCILL
 #  license as circulated by CEA, CNRS and INRIA at the following URL :
 #  "http://www.cecill.info".
 #
 #  As a counterpart to the access to the source code and rights to copy,
 #  modify and redistribute granted by the license, users are provided only
 #  with a limited warranty  and the software's author,  the holder of the
 #  economic rights,  and the successive licensors  have only  limited
 #  liability.
 #
 #  In this respect, the user's attention is drawn to the risks associated
 #  with loading,  using,  modifying and/or developing or reproducing the
 #  software by the user in light of its specific status of free software,
 #  that may mean  that it is complicated to manipulate,  and  that  also
 #  therefore means  that it is reserved for developers  and  experienced
 #  professionals having in-depth computer knowledge. Users are therefore
 #  encouraged to load and test the software's suitability as regards their
 #  requirements in conditions enabling the security of their systems and/or
 #  data to be ensured and,  more generally, to use and operate it in the
 #  same conditions as regards security.
 #
 #  The fact that you are presently reading this means that you have had
 #  knowledge of the CeCILL license and that you accept its terms.
*/

/* This code was compiled using
gcc -march=native -O2 condat_l1ballproj.c -o main  -I/usr/local/include/
	-lm -lgsl  -L/usr/local/lib/
On my machine, gcc is actually a link to the compiler Apple LLVM version 5.1
(clang-503.0.40) */

/* The following functions are implemented:
l1ballproj_IBIS
l1ballproj_Condat (proposed algorithm)
These functions take the same parameters. They project the vector y onto
the closest vector x of same length (parameter N in the paper) with
sum_{n=0}^{N-1}|x[n]|<=a.
We must have length>=1.
For l1ballproj_Condat only:
we can have x==y (projection done in place). If x!=y, the arrays x and y must
not overlap, as x is used for temporary calculations before y is accessed.
*/

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#define datatype double /* type of the elements in y */

/* Proposed algorithm */
static void l1ballproj_Condat(datatype* y, datatype* x,
const unsigned int length, const double a) {
	if (a<=0.0) {
		if (a==0.0) memset(x,0,length*sizeof(datatype));
		return;
	}
	datatype*	aux = (x==y ? (datatype*)malloc(length*sizeof(datatype)) : x);
	int		auxlength=1;
	int		auxlengthold=-1;
	double	tau=(*aux=(*y>=0.0 ? *y : -*y))-a;
	int 	i=1;
	for (; i<length; i++)
		if (y[i]>0.0) {
			if (y[i]>tau) {
				if ((tau+=((aux[auxlength]=y[i])-tau)/(auxlength-auxlengthold))
				<=y[i]-a) {
					tau=y[i]-a;
					auxlengthold=auxlength-1;
				}
				auxlength++;
			}
		} else if (y[i]!=0.0) {
			if (-y[i]>tau) {
				if ((tau+=((aux[auxlength]=-y[i])-tau)/(auxlength-auxlengthold))
				<=aux[auxlength]-a) {
					tau=aux[auxlength]-a;
					auxlengthold=auxlength-1;
				}
				auxlength++;
			}
		}
	if (tau<=0) {	/* y is in the l1 ball => x=y */
		if (x!=y) memcpy(x,y,length*sizeof(datatype));
		else free(aux);
	} else {
		datatype*  aux0=aux;
		if (auxlengthold>=0) {
			auxlength-=++auxlengthold;
			aux+=auxlengthold;
			while (--auxlengthold>=0)
				if (aux0[auxlengthold]>tau)
					tau+=((*(--aux)=aux0[auxlengthold])-tau)/(++auxlength);
		}
		do {
			auxlengthold=auxlength-1;
			for (i=auxlength=0; i<=auxlengthold; i++)
				if (aux[i]>tau)
					aux[auxlength++]=aux[i];
				else
					tau+=(tau-aux[i])/(auxlengthold-i+auxlength);
		} while (auxlength<=auxlengthold);
		for (i=0; i<length; i++)
			x[i]=(y[i]-tau>0.0 ? y[i]-tau : (y[i]+tau<0.0 ? y[i]+tau : 0.0));
		if (x==y) free(aux0);
	}
}

/* function eplb of the package SLEP, implementing the algorithm IBIS of the paper
J. Liu and J. Ye. Efficient Euclidean Projections in Linear Time, ICML 2009.
Version without the initial guess lambda0 of lambda (tau in our notations) */
static void l1ballproj_IBIS(datatype* y, datatype* x,
const unsigned int length, const double a) {
    #define delta 1e-13
    int i, j, flag=0;
    int rho_1, rho_2, rho, rho_T, rho_S;
    int V_i_b, V_i_e, V_i;
    double lambda_1, lambda_2, lambda_T, lambda_S, lambda;
    double s_1, s_2, s, s_T, s_S, y_max, temp;
    double f_lambda_1, f_lambda_2, f_lambda, f_lambda_T, f_lambda_S;
    int iter_step=0;

    /* find the maximal absolute value in y
     * and copy the (absolute) values from y to x
     */

	/*if (a< 0){
		printf("\n a should be nonnegative!");
		return;
	}*/

    V_i=0;
    if (y[0] !=0){
        rho_1=1;
        s_1=x[V_i]=y_max=fabs(y[0]);
        V_i++;
    }
    else{
        rho_1=0;
        s_1=y_max=0;
    }

    for (i=1;i<length; i++){
        if (y[i]!=0){
            x[V_i]=fabs(y[i]); s_1+= x[V_i]; rho_1++;

            if (x[V_i] > y_max)
                y_max=x[V_i];
            V_i++;
        }
    }

    /* If ||y||_1 <= a, then y is the solution  */
    if (s_1 <= a){
        flag=1;        lambda=0;
        for(i=0;i<length;i++){
            x[i]=y[i];
        }
        /**root=lambda;
        *steps=iter_step;*/
        return;
    }

    lambda_1=0; lambda_2=y_max;
    f_lambda_1=s_1 -a;
    /*f_lambda_1=s_1-rho_1* lambda_1 -a;*/
    rho_2=0; s_2=0; f_lambda_2=-a;
    V_i_b=0; V_i_e=V_i-1;

    /*lambda=lambda0;
    if ( (lambda<lambda_2) && (lambda> lambda_1) ){
         i=V_i_b; j=V_i_e; rho=0; s=0;
        while (i <= j){
            while( (i <= V_i_e) && (x[i] <= lambda) ){
                i++;
            }
            while( (j>=V_i_b) && (x[j] > lambda) ){
                s+=x[j];
                j--;
            }
            if (i<j){
                s+=x[i];

                temp=x[i];  x[i]=x[j];  x[j]=temp;
                i++;  j--;
            }
		}

        rho=V_i_e-j;  rho+=rho_2;  s+=s_2;
		f_lambda=s-rho*lambda-a;

        if ( fabs(f_lambda)< delta ){
            flag=1;
		}

		if (f_lambda <0){
			lambda_2=lambda; s_2=s;	rho_2=rho; f_lambda_2=f_lambda;

			V_i_e=j;  V_i=V_i_e-V_i_b+1;
		}
		else{
			lambda_1=lambda; rho_1=rho;	s_1=s; f_lambda_1=f_lambda;

			V_i_b=i; V_i=V_i_e-V_i_b+1;
		}

		if (V_i==0){
			lambda=(s - a)/ rho;
			flag=1;
		}
    }*/

    while (!flag){
        iter_step++;

        /* compute lambda_T  */
        lambda_T=lambda_1 + f_lambda_1 /rho_1;
        if(rho_2 !=0){
            if (lambda_2 + f_lambda_2 /rho_2 >	lambda_T)
                lambda_T=lambda_2 + f_lambda_2 /rho_2;
        }

        /* compute lambda_S */
        lambda_S=lambda_2 - f_lambda_2 *(lambda_2-lambda_1)/(f_lambda_2-f_lambda_1);

        if (fabs(lambda_T-lambda_S) <= delta){
            lambda=lambda_T; flag=1;
            break;
        }

        /* set lambda as the middle point of lambda_T and lambda_S */
        lambda=(lambda_T+lambda_S)/2;

        s_T=s_S=s=0;
        rho_T=rho_S=rho=0;
        i=V_i_b; j=V_i_e;
        while (i <= j){
            while( (i <= V_i_e) && (x[i] <= lambda) ){
                if (x[i]> lambda_T){
                    s_T+=x[i]; rho_T++;
                }
                i++;
            }
            while( (j>=V_i_b) && (x[j] > lambda) ){
                if (x[j] > lambda_S){
                    s_S+=x[j]; rho_S++;
                }
                else{
                    s+=x[j];  rho++;
                }
                j--;
            }
            if (i<j){
                if (x[i] > lambda_S){
                    s_S+=x[i]; rho_S++;
                }
                else{
                    s+=x[i]; rho++;
                }

                if (x[j]> lambda_T){
                    s_T+=x[j]; rho_T++;
                }

                temp=x[i]; x[i]=x[j];  x[j]=temp;
                i++; j--;
            }
		}

        s_S+=s_2; rho_S+=rho_2;
        s+=s_S; rho+=rho_S;
        s_T+=s; rho_T+=rho;
        f_lambda_S=s_S-rho_S*lambda_S-a;
        f_lambda=s-rho*lambda-a;
        f_lambda_T=s_T-rho_T*lambda_T-a;

        /*printf("\n %d & %d  & %5.6f & %5.6f & %5.6f & %5.6f & %5.6f \\\\ \n \\hline ", iter_step, V_i, lambda_1, lambda_T, lambda, lambda_S, lambda_2);*/

        if ( fabs(f_lambda)< delta ){
            /*printf("\n lambda");*/
            flag=1;
            break;
        }
        if ( fabs(f_lambda_S)< delta ){
           /* printf("\n lambda_S");*/
            lambda=lambda_S; flag=1;
            break;
        }
        if ( fabs(f_lambda_T)< delta ){
           /* printf("\n lambda_T");*/
            lambda=lambda_T; flag=1;
            break;
        }

        /*
        printf("\n\n f_lambda_1=%5.6f, f_lambda_2=%5.6f, f_lambda=%5.6f",f_lambda_1,f_lambda_2, f_lambda);
        printf("\n lambda_1=%5.6f, lambda_2=%5.6f, lambda=%5.6f",lambda_1, lambda_2, lambda);
        printf("\n rho_1=%d, rho_2=%d, rho=%d ",rho_1, rho_2, rho);
         */

        if (f_lambda <0){
            lambda_2=lambda;  s_2=s;  rho_2=rho;
            f_lambda_2=f_lambda;

            lambda_1=lambda_T; s_1=s_T; rho_1=rho_T;
            f_lambda_1=f_lambda_T;

            V_i_e=j;  i=V_i_b;
            while (i <= j){
                while( (i <= V_i_e) && (x[i] <= lambda_T) ){
                    i++;
                }
                while( (j>=V_i_b) && (x[j] > lambda_T) ){
                    j--;
                }
                if (i<j){
                    x[j]=x[i];
                    i++;   j--;
                }
            }
            V_i_b=i; V_i=V_i_e-V_i_b+1;
        }
        else{
            lambda_1=lambda;  s_1=s; rho_1=rho;
            f_lambda_1=f_lambda;

            lambda_2=lambda_S; s_2=s_S; rho_2=rho_S;
            f_lambda_2=f_lambda_S;

            V_i_b=i;  j=V_i_e;
            while (i <= j){
                while( (i <= V_i_e) && (x[i] <= lambda_S) ){
                    i++;
                }
                while( (j>=V_i_b) && (x[j] > lambda_S) ){
                    j--;
                }
                if (i<j){
                    x[i]=x[j];
                    i++;   j--;
                }
            }
            V_i_e=j; V_i=V_i_e-V_i_b+1;
        }

        if (V_i==0){
            lambda=(s - a)/ rho; flag=1;
            /*printf("\n V_i=0, lambda=%5.6f",lambda);*/
            break;
        }
    }/* end of while */

    for(i=0;i<length;i++){
        if (y[i] > lambda)
            x[i]=y[i]-lambda;
        else
            if (y[i]< -lambda)
                x[i]=y[i]+lambda;
            else
                x[i]=0;
    }
    /**root=lambda;
    *steps=iter_step;*/
}

/* Example of code to compare the computation times
Result on my machine: my algorithm is 2.5x faster than IBIS */
int main() {
	datatype*	y;
	datatype*	x;
	double	val=0.0;
	double 	aux;
	double 	epsi;
    int 	i,j;
    int 	Nbrea;
    unsigned int length;
    clock_t	start, end;
    const gsl_rng_type *T;
  	gsl_rng *r;
  	gsl_rng_env_setup();
	gsl_rng_default_seed=32067;
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);
	double timemean=0.0, timevar=0.0, timedelta;
	const double a = 1.0;
	srand((unsigned int)pow(time(NULL)%100,3));
    double*	 timetable;

    length = 1000000;
	Nbrea = 100;
	y=(datatype*)malloc(length*sizeof(datatype));
	x=(datatype*)malloc(length*sizeof(datatype));
	timetable=(double*)malloc((Nbrea+1)*sizeof(double));
	for (j=0; j<=Nbrea; j++) {
		for (i=0; i<length; i++) {
	    	y[i]=gsl_ran_gaussian(r, 0.01);
	    }
	    y[(int)(rand() / (((double)RAND_MAX+1.0)/length))]+=a;
	    start=clock();
	    l1ballproj_Condat(y,x,length,a);
	    /*l1ballproj_IBIS(y,x,length,a); */
	    end=clock();
	    timetable[j]=(double)(end-start)/CLOCKS_PER_SEC;
	}
	/* we discard the first value, because the computation time is higher,
	probably because of cache operations */
	for (j=1; j<=Nbrea; j++) {
		timedelta=timetable[j]-timemean;
		timemean+=timedelta/j;
		timevar+=timedelta*(timetable[j]-timemean);
	}
	timevar=sqrt(timevar/Nbrea);
	printf("av. time: %e std dev: %e\n",timemean,timevar);
	free(x);
    free(y);
    return 0;
}


text:
//  This matrix class is a C++ wrapper for the GNU Scientific Library
//  Copyright (C)  ULP-IPB Strasbourg

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software
//  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

#ifndef __matrix_vector_operators_h
#define __matrix_vector_operators_h

#include "gsl/gsl_blas.h"
#include <gslwrap/matrix_double.h>
#include <gslwrap/matrix_float.h>
#include <gslwrap/vector_double.h>

namespace gsl
{

inline
vector_float operator*(const matrix_float& m, const vector_float& v)
{
	vector_float y(m.get_rows());
	gsl_blas_sgemv(CblasNoTrans, 1.0, m.gslobj(), v.gslobj(), 0.0, y.gslobj());
	return y;
}

inline
vector operator*(const matrix& m, const vector& v)
{
	vector y(m.get_rows());
	gsl_blas_dgemv(CblasNoTrans, 1.0, m.gslobj(), v.gslobj(), 0.0, y.gslobj());
	return y;
}

}

#endif //__matrix_vector_operators_h

ext: h

max_stars_repo_name: entn-at/GlottDNN

max_issues_repo_name: entn-at/GlottDNN

max_forks_repo_name: entn-at/GlottDNN


text:
#include <assert.h>
#include <cblas.h>
#include <complex.h>
#include <errno.h>
#include <fftw3.h>
#include <getopt.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <math.h>
#include <omp.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#define C 299792458
#define PI 3.141592654
#define ANSI_COLOR_RED     "\x1b[31m"
#define ANSI_COLOR_GREEN   "\x1b[32m"
#define ANSI_COLOR_YELLOW  "\x1b[33m"
#define ANSI_COLOR_BLUE    "\x1b[34m"
#define ANSI_COLOR_MAGENTA "\x1b[35m"
#define ANSI_COLOR_CYAN    "\x1b[36m"
#define ANSI_COLOR_RESET   "\x1b[0m"
#define _MODEL_DEBUG 0
#define MAX_TX 1
#define MAX_RX 1
#define MAX_RIBBON_SIZE 20
#define MAX_SURFACES 20
// if max values exceeded, may cause memory bugs

struct simulation {
        struct environment *env;
        struct file_reader *fr;
};

struct spatial_motion_model {
        double velocity[3];
        double position[3];
};

struct transmission_model {
        double power_in_dBm;
};

struct propagation_model {
        double distance;
};

struct general_node {
        struct spatial_motion_model *smm;
        struct transmission_model *tm;
        int id;
};

struct transmitter {
        struct general_node *gn;
        double complex baseband_signal;
};

struct receiver {
        struct general_node *gn;
        double recv_noise_power;
        double complex rx_signal;
        struct receiver_ray_ribbon_ll_node *rlln;
};

struct perfect_reflector {
        double unit_normal[3];
        double unit_length_normal[3];
        double unit_width_normal[3];
        double center_point[3];
        double length, width;
};

struct environment {
        struct receiver **receivers_array;
        struct transmitter **transmitters_array;
        struct perfect_reflector **prarray;
        struct general_node **node_array;
        struct ray_ribbon_array **env_paths;
        struct ray_ribbon_array **tx_paths;

        double recv_unit_normal[3];
        double time;
        double complex *unit_power_gaussian_noise;
        double frequency;
        double wavelength;
        double delta_time;
        double end_time;
        double max_limit;
        double min_limit;
	int refresh_time;
	int awgn_supplied;
        int time_index;

        int num_transmitters;
        int num_receivers;
        int num_reflectors;

        int sz_array_tx;
        int sz_array_rx;
        int sz_array_gn;
        int sz_array_pr;

        // flags
        bool read_in_nodes;
        /* should be true by the time TX or RX is read in */
        bool node_memory_allocated;
        bool tx_paths_updated;
        bool tx_paths_updated_rx_paths_updated;
};

struct file_reader {
        FILE *infile;
        FILE *outfile;

        // filenames longer than 999 will be truncated
        char input_filename[1000];
        char output_filename[1000];
};

int id();
struct simulation *init_simulation();
struct spatial_motion_model *init_spatial_motion_model();
struct transmission_model *init_transmission_model();
struct propagation_model *init_propagation_model();
struct general_node *init_general_node();
struct transmitter *init_transmitter();
struct receiver *init_receiver();
struct perfect_reflector *init_perfect_reflector(const double *normal,
                                         const double *center_point,
                                         const double *length_normal,
                                         double length, double width);
struct perfect_reflector **init_perfect_reflectorarray(int number);
struct environment *init_environment();
int malloc_environment(struct environment *env);

void print_vector(const double *db);
void print_spatial_motion_model(const struct spatial_motion_model *smm);
void print_transmission_model(const struct transmission_model *tm);
void print_propagation_model(const struct propagation_model *pm);
void print_general_node(const struct general_node *gn);
void print_transmitter(const struct transmitter *tx);
void print_receiver(const struct receiver *rx);
void print_perfect_reflectors(const struct perfect_reflector *pr);
void print_environment(const struct environment *env);
void print_env_paths(const struct environment *env);
void print_tx_paths(const struct environment *env);

void destroy_spatial_motion_model(struct spatial_motion_model *smm);
void destroy_simulation(struct simulation *sim);
void destroy_transmission_model(struct transmission_model *tm);
void destroy_propagation_model(struct propagation_model *pm);
void destroy_transmitter(struct transmitter *tx);
void destroy_receiver(struct receiver *tx);
void destroy_environment(struct environment *env);
void destroy_file_reader(struct file_reader *fr);
void destroy_perfect_reflector(struct perfect_reflector *pr);
void destroy_perfect_reflectorarray(struct perfect_reflector **pr_begin);
struct file_reader *init_file_reader(int argc, char *argv[]);
void cross_product(const double *v1, const double *v2, double *v3);
double normalize_unit_vector(double *v1);
void diff(const double *v1, const double *v2, double *v3);
int find_len(void **ptr);

void add_receiver_patch(struct environment *env, int length);
void destroy_last_reflector(struct environment *env);
double distance(const struct general_node *gn1,
                const struct general_node *gn2);

bool update_environment_from_file_sim(struct simulation *sim);
bool handle_request(struct environment *env, FILE *fp, const char *req_type);
bool custom_fscanf(FILE *fp, const char *str, void *ptr);

// Deprecated
bool update_environment_from_file(struct environment *env,
                                FILE *fp);


text:
#pragma once

#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#else
#include "FakeGSL.h"
#endif

#include <vector>

#include <fstream>
#include <iostream>
#include "BasicError.h"
//#define _NOSNOPT 1
#ifndef _NOSNOPT
#include "Snopt.h"
#else
#include "CustomSolver.h"
#endif

using namespace std;
typedef int (* MAX_SOLVER_DF_TYPE)(integer*, doublereal*, integer*, doublereal*, integer*, doublereal*, char*, integer*);

class MaxSolver {
	integer n; integer neF; integer lenG;

	integer *iGfun;
	integer *jGvar;
	doublereal* G;

	doublereal *x;
	doublereal *xlow;
	doublereal *xupp;

	doublereal *F;
	doublereal *Flow;
	doublereal *Fupp;

	gsl_vector* result;
	double objectiveVal;

	MAX_SOLVER_DF_TYPE df;
	char* workspace;

	doublereal obj;
	doublereal* grad;
	doublereal* oldx;

	doublereal stepSize = 0.1;

	int maxSteps = 5000;

public:
	MaxSolver(integer n_, integer neF_): n(n_), neF(neF_){
		lenG = n*neF;
		iGfun = new integer[lenG];
		jGvar = new integer[lenG];
		G = new doublereal[lenG];

		grad = new doublereal[n];

		x = new doublereal[n];
		xlow = new doublereal[n];
		xupp = new doublereal[n];

		F = new doublereal[neF];
		Flow = new doublereal[neF];
		Fupp = new doublereal[neF];

		oldx = new doublereal[n];
		result = gsl_vector_alloc(n);
	}

	~MaxSolver() {
		delete []iGfun; delete []jGvar;
		delete []x; delete []xlow; delete []xupp;
		delete []F; delete []Flow; delete []Fupp;
	}

	void init(char* workspace, integer nef_, MAX_SOLVER_DF_TYPE _df, doublereal *xlow_, doublereal *xupp_, doublereal *Flow_, doublereal *Fupp_);

	bool optimize(const gsl_vector* initState, const gsl_vector* initDir, bool suppressPrint = false);

	gsl_vector* getResults() {
		return result;
	}

	double getObjectiveVal() {
		return objectiveVal;
	}

	bool inLimits();

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
// Wrapper of the GSL interpolation functions that lets you define a
// 1D function by points and then query points on that function.
//
// Author: Mark Desnoyer (mdesnoyer@gmail.com)
// Date: Oct 2012

#ifndef __INTERPOLATION_H__
#define __INTERPOLATION_H__
#include <vector>
#include <gsl/gsl_spline.h>
#include <assert.h>
#include <algorithm>
#include <limits>
#include <exception>

template <typename T>
class SortOrder {
public:
  SortOrder(const std::vector<T>& sortArray) : sortArray_(sortArray) {}

  bool operator()(int lhs, int rhs) const {
    return sortArray_[lhs] < sortArray_[rhs];
  }

private:
  const std::vector<T>& sortArray_;
};

// Classes that inherit this abstract class must intialize spline_ in
// their constructors.
class Interpolator {
public:
  class out_of_bounds : public std::exception {};

  template <typename InputIterator>
  Interpolator(InputIterator xFirst, InputIterator xLast,
               InputIterator yFirst, InputIterator yLast)
    : x_(), y_(), acc_(gsl_interp_accel_alloc()) {
    std::vector<double> x(xFirst, xLast);
    std::vector<double> y(yFirst, yLast);
    std::vector<int> idx;
    for (unsigned int i = 0u; i < x.size(); ++i) {
      idx.push_back(i);
    }

    assert(x.size() == y.size());

    // Get the sorting indicies based on the x matrix
    std::sort(idx.begin(), idx.end(), SortOrder<double>(x));

    // Populate the data in ascending x order skipping equal values
    double curX = std::numeric_limits<double>::quiet_NaN();
    for (unsigned int i = 0u; i < idx.size(); ++i) {
      if (curX != x[idx[i]]) {
        x_.push_back(x[idx[i]]);
        y_.push_back(y[idx[i]]);
      }
      curX = x[idx[i]];
    }

  }

  virtual ~Interpolator();

  // Evaluate the interpolation at xi
  double operator()(double xi) const {
    if (xi < x_[0] || xi > x_.back()) {
      throw out_of_bounds();
    }
    return gsl_spline_eval(spline_, xi, acc_);
  }

  // Evaluate the derivative at xi
  double EvalDeriv(double xi) const {
    if (xi < x_[0] || xi > x_.back()) {
      throw out_of_bounds();
    }
    return gsl_spline_eval_deriv(spline_, xi, acc_);
  }

  double minY() const { return y_.front(); }
  double maxY() const { return y_.back(); }

protected:
  std::vector<double> x_;
  std::vector<double> y_;

  gsl_interp_accel* acc_;
  gsl_spline* spline_;

};

class LinearInterpolator : public Interpolator {
public:
  LinearInterpolator(const std::vector<double>& x, const std::vector<double>& y)
    : Interpolator(x.begin(), x.end(), y.begin(), y.end()) {
    InitSpline();
  }

  template <typename InputIterator>
  LinearInterpolator(InputIterator xFirst, InputIterator xLast,
                     InputIterator yFirst, InputIterator yLast)
    : Interpolator(xFirst, xLast, yFirst, yLast) {
    InitSpline();
  }

  virtual ~LinearInterpolator();

private:
  void InitSpline();

};

class SplineInterpolator : public Interpolator {
public:
  SplineInterpolator(const std::vector<double>& x, const std::vector<double>& y)
    : Interpolator(x.begin(), x.end(), y.begin(), y.end()) {
    InitSpline();
  }

  template <typename InputIterator>
  SplineInterpolator(InputIterator xFirst, InputIterator xLast,
                     InputIterator yFirst, InputIterator yLast)
    : Interpolator(xFirst, xLast, yFirst, yLast) {
    InitSpline();
  }

  virtual ~SplineInterpolator();

private:
  void InitSpline();

};

#endif // __INTERPOLATION_H__


text:
#pragma once

#include <memory>
#include <random>
#include <gsl/gsl_rng.h>

namespace EERAModel {
/**
 * @brief Namespace containing objects related to random number generation
 *
 * This namespace contains a wrapper interace to existing random number
 * generation libraries in order for them to be used within all
 * model functions.
 */
namespace Random {

/**
 * @class RNGInterface
 * @brief Abstract interface to a random number generator
 */
class RNGInterface {
public:
    using Sptr = std::shared_ptr<RNGInterface>;

    virtual ~RNGInterface() = default;

    /**
     * @brief Flat distribution
     *
     * Return a random number from a flat distribution
     *
     * @param a Lower limit of distribution interval
     * @param b Upper limit of distribution interval
     *
     * @return Random number
     */
    virtual double Flat(double a, double b) = 0;

    /**
     * @brief Multinomial distribution
     *
     * Return a random number from a multinomial distribution
     */
    virtual void Multinomial(size_t K, unsigned int N, const double p[], unsigned int n[]) = 0;

    /**
     * @brief Poisson distribution
     *
     * Return a random number from a Poisson distribution
     *
     * @param mu Mean of the distribution
     *
     * @return Random number
     */
    virtual double Poisson(double mu) = 0;

    /**
     * @brief Gamma distribution
     *
     * Return a random number from a Gamma distribution
     *
     * @param a Gamma a parameter
     * @param b Gamma b parameter
     *
     * @return Random number
     */
    virtual double Gamma(double a, double b) = 0;

    /**
     * @brief Beta distribution
     *
     * Return a random number from a Beta distribution
     *
     * @param a Beta a parameter
     * @param b Beta b parameter
     *
     * @return Random number
     */
    virtual double Beta(double a, double b) = 0;

   /**
    * @brief Binomial distribution
    *
    * Return a random number from a Binomial distribution
    *
    * @param p Binomial p parameter: probability
    * @param n Binomial n parameter: the population to sample
    *
    * @return Random number
    */
   virtual unsigned int Binomial(double p, unsigned int n) = 0;

    /**
     * @brief MT19937 generator
     *
     * Returns a reference to a std::mt19937 object
     *
     * @return mt19937 reference
     */
    virtual std::mt19937& MT19937() = 0;
};

/**
 * @class RNG
 * @brief Wrapper for the GNU Scientific Library random number generator
 */
class RNG : public RNGInterface
{
public:
    using Sptr = std::shared_ptr<RNG>;

    /**
     * @brief Contructor
     *
     * Sets up a handle for the underlying GSL and STL random number generator. Seeds the randomiser
     * with the supplied seed
     *
     * @param seed Seed value
     */
    RNG(unsigned long int seed);

    /** @brief Flat distribution override */
    virtual double Flat(double a, double b) override;

    /** @brief Multinomial distribution override */
    virtual void Multinomial(size_t K, unsigned int N, const double p[], unsigned int n[]) override;

    /** @brief Poisson distribution override */
    virtual double Poisson(double mu) override;

    /** @brief Gamma distribution override */
    virtual double Gamma(double a, double b) override;

    /** @brief Beta distribution override */
    virtual double Beta(double a, double b) override;

    /** @brief Binomial distribution override */
    virtual unsigned int Binomial(double p, unsigned int n) override;

    /** @brief MT19937 override */
    virtual std::mt19937& MT19937() override;

private:
    /**
     * @private
     * @brief Handle for the GSL random number generator
     */
    gsl_rng* r_;

    /**
     * @private
     * @brief Mersenne twister generator
     */
    std::mt19937 gen_;
};

} // namespace Random
} // namespace EERAModel


text:
#ifndef __LIB_MATH_OPTIMIZER__
#define __LIB_MATH_OPTIMIZER__

#include <iostream>
#include <limits>
#include <gsl/gsl_multimin.h>
#include "math/gsl_vector.h"
#include "util/flags.h"

using std::cout;
using std::endl;

DEFINE_double(multimin_convergence_threshold,
	      1e-5,
	      "Convergence threshold for conjugate gradient.");
DEFINE_size(max_multimin_iterations,
	    40,
	    "Maximum number of conjugate gradient iterations to perform.");

class Optimizer {
 public:
  Optimizer(size_t size) : size_(size) {
  }

  void Optimize() {
    gsl_multimin_function_fdf my_func;
    my_func.n = size_;
    my_func.f = &MultiminObjectiveWrapper;
    my_func.df = &MultiminGradientWrapper;
    my_func.fdf = &MultiminObjectiveGradientWrapper;
    my_func.params = this;

    gsl_multimin_fdfminimizer* s =
      gsl_multimin_fdfminimizer_alloc(gsl_multimin_fdfminimizer_conjugate_fr, size_);
    GslVector initial_guess(size_);
    MultiminInitialGuess(initial_guess.mutable_ptr());

    // step_size, tol
    //    gsl_multimin_fdfminimizer_set(s, &my_func, initial_guess.ptr(), 0.1, 1.0);
    gsl_multimin_fdfminimizer_set(s, &my_func, initial_guess.ptr(), 0.01, 0.01);

    size_t iter = 0;
    int status;

    double value = std::numeric_limits<double>::infinity();
    double prev_value;
    do {
      prev_value = value;
      iter++;
      status = gsl_multimin_fdfminimizer_iterate(s);
      if (status) {
	cout << "Error: " << gsl_strerror(status) << endl;
	break;
      }
      status = gsl_multimin_test_gradient(s->gradient, 1e-3);
      if (status == GSL_SUCCESS) {
	cout << "Minimum found." << endl;
      }
      value = s->f;
      cout << "Iteration: " << iter << " Value: " <<
	value << " dValue:" << (prev_value - value)/fabs(value) << " " <<
	gsl_strerror(status) << endl;
    } while (status == GSL_CONTINUE &&
	       iter < FLAGS_max_multimin_iterations &&
	     (prev_value - value) / fabs(value) > FLAGS_multimin_convergence_threshold);
    MultiminResult(s->x);
    gsl_multimin_fdfminimizer_free(s);
  }

  virtual void MultiminObjectiveGradient(const gsl_vector* x,
					 double* objective,
					 gsl_vector* gradient) = 0;

  virtual void MultiminInitialGuess(gsl_vector* v) = 0;

  virtual void MultiminResult(gsl_vector* x) = 0;

  virtual ~Optimizer() { }
 protected:
  static double MultiminObjectiveWrapper(const gsl_vector* x, void* params) {
    double objective;
    reinterpret_cast<Optimizer*>(params)->MultiminObjectiveGradient(x, &objective, NULL);
    return objective;
  }

  static void MultiminGradientWrapper(const gsl_vector* x, void* params, gsl_vector* g) {
    reinterpret_cast<Optimizer*>(params)->MultiminObjectiveGradient(x, NULL, g);
  }

  static void MultiminObjectiveGradientWrapper(const gsl_vector* x,
					       void* params,
					       double* f,
					       gsl_vector* g) {
    reinterpret_cast<Optimizer*>(params)->MultiminObjectiveGradient(x, f, g);
  }

 private:
  size_t size_;
};
#endif  // __LIB_MATH_OPTIMIZER__

max_stars_repo_name: boomsbloom/dtm-fmri

max_issues_repo_name: boomsbloom/dtm-fmri

max_forks_repo_name: boomsbloom/dtm-fmri


text:
/**
 *
 * @file testing_dgemm.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Emmanuel Agullo
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:45:18 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_dmain.h"

#undef COMPLEX
#define REAL

static int check_solution(PLASMA_enum transA, PLASMA_enum transB, int M, int N, int K,
                          double alpha, double *A, int LDA,
                          double *B, int LDB,
                          double beta, double *Cref, double *Cplasma, int LDC);

int testing_dgemm(int argc, char **argv)
{
    /* Check for number of arguments*/
    if ( argc != 8) {
        USAGE("GEMM", "alpha beta M N K LDA LDB LDC",
              "   - alpha  : alpha coefficient\n"
              "   - beta   : beta coefficient\n"
              "   - M      : number of rows of matrices A and C\n"
              "   - N      : number of columns of matrices B and C\n"
              "   - K      : number of columns of matrix A / number of rows of matrix B\n"
              "   - LDA    : leading dimension of matrix A\n"
              "   - LDB    : leading dimension of matrix B\n"
              "   - LDC    : leading dimension of matrix C\n");
        return -1;
    }

    PLASMA_Set(PLASMA_TILE_SIZE, 128);
    double alpha = (double) atol(argv[0]);
    double beta = (double) atol(argv[1]);
    int M     = atoi(argv[2]);
    int N     = atoi(argv[3]);
    int K     = atoi(argv[4]);
    int LDA   = atoi(argv[5]);
    int LDB   = atoi(argv[6]);
    int LDC   = atoi(argv[7]);

    double eps;
    int info_solution;
    int i, j, ta, tb;
    int LDAxK = LDA*max(M,K);
    int LDBxN = LDB*max(K,N);
    int LDCxN = LDC*N;

    double *A      = (double *)malloc(LDAxK*sizeof(double));
	#pragma omp register( [LDAxK]A )
    double *B      = (double *)malloc(LDBxN*sizeof(double));
	#pragma omp register( [LDBxN]B )
    double *C      = (double *)malloc(LDCxN*sizeof(double));
	#pragma omp register( [LDCxN]C )
    double *Cinit  = (double *)malloc(LDCxN*sizeof(double));
	#pragma omp register( [LDCxN]Cinit )
    double *Cfinal = (double *)malloc(LDCxN*sizeof(double));
	#pragma omp register( [LDCxN]Cfinal )

    /* Check if unable to allocate memory */
    if ((!A)||(!B)||(!Cinit)||(!Cfinal)){
        printf("Out of Memory \n ");
        return -2;
    }

    eps = LAPACKE_dlamch_work('e');
    printf("\n");
    printf("------ TESTS FOR PLASMA DGEMM ROUTINE -------  \n");
    printf("            Size of the Matrix %d by %d\n", M, N);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n",eps);
    printf(" Computational tests pass if scaled residuals are less than 10.\n");

    /*----------------------------------------------------------
     *  TESTING DGEMM
     */

    /* Initialize A, B, C */
    LAPACKE_dlarnv_work(IONE, ISEED, LDAxK, A);
    LAPACKE_dlarnv_work(IONE, ISEED, LDBxN, B);
    LAPACKE_dlarnv_work(IONE, ISEED, LDCxN, C);

#ifdef COMPLEX
    for (ta=0; ta<3; ta++) {
        for (tb=0; tb<3; tb++) {
#else
    for (ta=0; ta<2; ta++) {
        for (tb=0; tb<2; tb++) {
#endif
            for ( i = 0; i < M; i++)
                for (  j = 0; j < N; j++)
                    Cinit[LDC*j+i] = C[LDC*j+i];
            for ( i = 0; i < M; i++)
                for (  j = 0; j < N; j++)
                    Cfinal[LDC*j+i] = C[LDC*j+i];

            /* PLASMA DGEMM */
            PLASMA_dgemm(trans[ta], trans[tb], M, N, K, alpha, A, LDA, B, LDB, beta, Cfinal, LDC);

            /* Check the solution */
            info_solution = check_solution(trans[ta], trans[tb], M, N, K,
                                           alpha, A, LDA, B, LDB, beta, Cinit, Cfinal, LDC);

            if (info_solution == 0) {
                printf("***************************************************\n");
                printf(" ---- TESTING DGEMM (%s, %s) ............... PASSED !\n", transstr[ta], transstr[tb]);
                printf("***************************************************\n");
            }
            else {
                printf("************************************************\n");
                printf(" - TESTING DGEMM (%s, %s) ... FAILED !\n", transstr[ta], transstr[tb]);
                printf("************************************************\n");
            }
        }
    }
#ifdef _UNUSED_
    }}
#endif
    free(A); free(B); free(C);
    free(Cinit); free(Cfinal);

    return 0;
}

/*--------------------------------------------------------------
 * Check the solution
 */

static int check_solution(PLASMA_enum transA, PLASMA_enum transB, int M, int N, int K,
                          double alpha, double *A, int LDA,
                          double *B, int LDB,
                          double beta, double *Cref, double *Cplasma, int LDC)
{
    int info_solution;
    double Anorm, Bnorm, Cinitnorm, Cplasmanorm, Clapacknorm, Rnorm, result;
    double eps;
    double beta_const;

    double *work = (double *)malloc(max(K,max(M, N))* sizeof(double));
    int Am, An, Bm, Bn;

    beta_const  = -1.0;

    if (transA == PlasmaNoTrans) {
        Am = M; An = K;
    } else {
        Am = K; An = M;
    }
    if (transB == PlasmaNoTrans) {
        Bm = K; Bn = N;
    } else {
        Bm = N; Bn = K;
    }

    Anorm       = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), Am, An, A,       LDA, work);
    Bnorm       = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), Bm, Bn, B,       LDB, work);
    Cinitnorm   = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M,  N,  Cref,    LDC, work);
    Cplasmanorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M,  N,  Cplasma, LDC, work);

    cblas_dgemm(CblasColMajor, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, M, N, K,
                (alpha), A, LDA, B, LDB, (beta), Cref, LDC);

    Clapacknorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Cref, LDC, work);

    cblas_daxpy(LDC * N, (beta_const), Cplasma, 1, Cref, 1);

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Cref, LDC, work);

    eps = LAPACKE_dlamch_work('e');

    printf("Rnorm %e, Anorm %e, Bnorm %e, Cinitnorm %e, Cplasmanorm %e, Clapacknorm %e\n",
           Rnorm, Anorm, Bnorm, Cinitnorm, Cplasmanorm, Clapacknorm);

    result = Rnorm / ((Anorm + Bnorm + Cinitnorm) * N * eps);
    printf("============\n");
    printf("Checking the norm of the difference against reference DGEMM \n");
    printf("-- ||Cplasma - Clapack||_oo/((||A||_oo+||B||_oo+||C||_oo).N.eps) = %e \n",
           result);

    if (  isnan(Rnorm) || isinf(Rnorm) || isnan(result) || isinf(result) || (result > 10.0) ) {
         printf("-- The solution is suspicious ! \n");
         info_solution = 1;
    }
    else {
         printf("-- The solution is CORRECT ! \n");
         info_solution= 0 ;
    }

    free(work);

    return info_solution;
}

int timing_dgemm(int argc, char **argv)
{

	int transa; int transb;
	int m, n, k;
	double alpha, beta;
	int lda, ldb, ldc;
	double *A, *B, *C;
	double *C0;
	int bs, rep;
	int i;

	double start, end;
	double elapsed;
	FILE *log;
	int num_threads;

	//gemm transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc
	if (argc != 12) {
		fprintf(stderr, "GEMMs transa transb m n k alpha lda ldb beta ldc bs rep\n");
		return 1;;
	}

	sscanf(argv[0], "%d", &transa);
	sscanf(argv[1], "%d", &transb);

	if ( !transa )
		transa = PlasmaNoTrans;
	else
		transa = PlasmaTrans;
	if ( !transb )
		transb = PlasmaNoTrans;
	else
		transb = PlasmaTrans;

	sscanf(argv[2], "%d", &m);
	sscanf(argv[3], "%d", &n);
	sscanf(argv[4], "%d", &k);
	sscanf(argv[5], "%lf", &alpha);
	sscanf(argv[6], "%d", &lda);
	sscanf(argv[7], "%d", &ldb);
	sscanf(argv[8], "%lf", &beta);
	sscanf(argv[9], "%d", &ldc);
	sscanf(argv[10], "%d", &bs);
	sscanf(argv[11], "%d", &rep);

	int dimA = max(m,k) * lda;
	int dimB = max(n,k) * ldb;
	int dimC = max(m,n) * ldc;
	A = malloc( dimA * sizeof(double));
	#pragma omp register ([dimA]A)
	B = malloc( dimB * sizeof(double));
	#pragma omp register ([dimB]B)
	C = malloc( dimC * sizeof(double));
	#pragma omp register ([dimC]C)

	LAPACKE_dlarnv(IONE, ISEED, dimA, A);
	LAPACKE_dlarnv(IONE, ISEED, dimB, B);
	LAPACKE_dlarnv(IONE, ISEED, dimC, C);

	PLASMA_Set( PLASMA_TILE_SIZE, bs);

	elapsed = 0.0;
	for ( i = 0; i < rep; i++ ) {
		start = gtime();
		PLASMA_dgemm(transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
		end = gtime();
		elapsed += end - start;
	}

	num_threads = omp_get_max_threads();
	dump_info("plasma_dgemm.log", num_threads, elapsed, rep);

	free(A); free(B); free(C);
    return 0;
}

ext: c


text:
#ifndef S3DCV_DEMO_DEMO_UTILS_H
#define S3DCV_DEMO_DEMO_UTILS_H

#include <s3d/cv/disparity/disparity_analyzer_stan.h>

#include <gsl/gsl>

class BadNumberOfArgumentsException {};
class FileNotFoundException {};

/**
 * Creates dataset loader from program input arguments.
 *
 * First input argument must be the dataset folder, e.g. 'Arch'
 * Second input argument must be the dataset path e.g. '/home/you/dataset'
 *
 */
inline DatasetLoaderVSG createDatasetLoader(int argc, char* argv[]) {
  auto input_args = gsl::span<char*>(argv, argc);
  if (input_args.size() != 3) {
    throw BadNumberOfArgumentsException{};
  }
  const std::string datasetName = input_args[1];
  const std::string datasetPath = input_args[2];

  return DatasetLoaderVSG{datasetName, datasetPath};
}

inline s3d::robust::Lmeds<s3d::StanFundamentalMatrixSolver, s3d::SampsonDistanceFunction> createRansac(
        s3d::Size imageSize) {
  int h = imageSize.getHeight();
  int w = imageSize.getWidth();

  s3d::robust::Parameters params;
  params.nbTrials = 500;
  params.distanceThreshold = 0.01 * sqrt(h * h + w * w);

  return s3d::DisparityAnalyzerSTAN::RansacAlgorithmSTAN(params);
}

#endif //S3DCV_DEMO_DEMO_UTILS_H


text:
#pragma once
#include <stdio.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>

typedef  unsigned int  uint;

typedef struct{
  double mu;
  double sigma;
} Gauss_params;

void Gauss_params_print( Gauss_params params );

void GSLfun_setup();

double GSLfun_ran_beta( double a, double b );
double GSLfun_ran_beta_Jeffreys();
uint   GSLfun_ran_binomial( double p, uint n );
double GSLfun_ran_gamma( double a, double theta );
double GSLfun_ran_gaussian( Gauss_params params );
double GSLfun_ran_gaussian_pdf( double x, Gauss_params params );

double gsl_ran_flat01();


text:
/* psl.c
////
Provide a library of gsl functions.

Features:
- request calculation via message
- function lookup
- inlets-on-demand

Author: shakfu
Repo: https://github.com/shakfu/pd-psl.git

*/
#include <math.h>
#include <string.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf_airy.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_clausen.h>
#include <gsl/gsl_sf_dawson.h>
#include <gsl/gsl_sf_debye.h>

#include "m_pd.h"
#include "tinyexpr.h"

// macros and defines
//  ---------------------------------------------------------------------------

#define MAX_ARGS 6
#define STR_BUF_SIZE 1000

// function lookup infratructure
// ---------------------------------------------------------------------------

unsigned long hash(const char *str) {
    unsigned int h = 0;
    int c;

    while ((c = *str++))
        h += (h << 1) + c;

    return h;
}

enum FUNC {
    ADD = 1273,
    LOG1P = 12931,
    EXPM1 = 12805,
    HYPOT = 13148,
    HYPOT3 = 39495,
    ACOSH = 11978,
    ASINH = 12341,
    ATANH = 12296,
    LDEXP = 12829,
    POW_INT = 122216,
    POW_2 = 13475,
    POW_3 = 13476,
    POW_4 = 13477,
    POW_5 = 13478,
    POW_6 = 13479,
    POW_7 = 13480,
    POW_8 = 13481,
    POW_9 = 13482,
    RANDO = 13254,
    FCMP = 4084,
    AIRY_AI = 109980,
    AIRY_BI = 109983,
    BESSEL_J0 = 987963,
    BESSEL_J1 = 987964,
    BESSEL_JN = 988025,
    BESSEL_Y0 = 988008,
    BESSEL_Y1 = 988009,
    BESSEL_YN = 988070,
    BESSEL_I0 = 987960,
    BESSEL_I1 = 987961,
    BESSEL_IN = 988022,
    CLAUSEN = 110879,
    DAWSON = 36848,
    DEBYE_1 = 109891,
    DEBYE_2 = 109892,
    DEBYE_3 = 109893,
    DEBYE_4 = 109894,
};

// forward declarations / prototypes
// ---------------------------------------------------------------------------

typedef struct _psl t_psl;

void select_default_function(t_psl *x, t_symbol *s);

// psl class objects
// ---------------------------------------------------------------------------

static t_class *psl_class;

static t_class *psl_inlet_class;

// psl class struct (data-space)
// ---------------------------------------------------------------------------

typedef void (*unary_func)(t_psl *, t_floatarg);
typedef void (*binary_func)(t_psl *, t_floatarg, t_floatarg);
typedef void (*tri_func)(t_psl *, t_floatarg, t_floatarg, t_floatarg);

typedef struct _psl_inlet
{
    t_class *x_pd;  // minimal pd object.
    t_psl   *owner; // the owning object to forward inlet messages to.
    int     id;     // the number of this inlet.
} t_psl_inlet;

typedef struct _psl {
    t_object x_obj;

    // assigned function
    t_symbol *func_name;
    int nargs;

    // function slots
    unary_func ufunc;
    binary_func bfunc;
    tri_func tfunc;

    // param_array
    t_float *arg_array;

    // for expression
    char expr_buffer[MAXPDSTRING];

    // inlets
    int inlets;          // # of extra inlets in addition to default
    t_psl_inlet *ins;    // the inlets themselves

    // outlets
    t_outlet *out_f;
} t_psl;

// psl class methods (operation-space)
// ---------------------------------------------------------------------------

// typed-methods

void psl_bang(t_psl *x) {
    if (x->nargs == 1 && x->inlets == 0) {
        x->ufunc(x, x->arg_array[0]);
    }

    if (x->nargs == 2 && x->inlets == 1) {
        x->bfunc(x, x->arg_array[0], x->arg_array[1]);
    }

    if (x->nargs == 3 && x->inlets == 2) {
        x->tfunc(x, x->arg_array[0], x->arg_array[1], x->arg_array[2]);

    }

}

void psl_float(t_psl *x, t_floatarg f) {
    post("psl_float: %f", f);
    if (x->nargs > 0) {
        x->arg_array[0] = f;
        psl_bang(x);
    } else {
        post("nothing to do: no function selected.");
        outlet_float(x->out_f, f);
    }
}

void psl_list(t_psl *x, t_symbol *s, int argc, t_atom *argv) {

    // atom_post("psl_list: ", argc, argv);

    if (s == gensym("list")) {
        post("s: list");

        if (argc == 0) {
            return;
        }

        if (argc == 1) {
            if (argv->a_type == A_FLOAT) {
                float f = atom_getfloat(argv);
                post("got float: %f", f);
                x->ufunc(x, f);
                return;
            }
        }

        if (argc == 2) {
            char buf[STR_BUF_SIZE];
            for (int i = 0; i < argc; i++) {
                atom_string((argv+i), buf, STR_BUF_SIZE);
                post("arg+%i: %s", i, buf);
            }
            if (argv->a_type == A_FLOAT && (argv + 1)->a_type == A_FLOAT) {
                float f1 = atom_getfloat(argv+0);
                float f2 = atom_getfloat(argv+1);
                post("f(%.2f, %.2f)", f1, f2);
                x->bfunc(x, f1, f2);
                return;
            }
        }

        if (argc == 3) {
            char buf[STR_BUF_SIZE];
            for (int i = 0; i < argc; i++) {
                atom_string((argv+i), buf, STR_BUF_SIZE);
                post("arg+%i: %s", i, buf);
            }
            if (argv->a_type == A_FLOAT && (argv+1)->a_type == A_FLOAT && (argv+2)->a_type == A_FLOAT) {
                float f1 = atom_getfloat(argv+0);
                float f2 = atom_getfloat(argv+1);
                float f3 = atom_getfloat(argv+2);
                post("f(%.2f, %.2f, %.2f)", f1, f2, f3);
                x->tfunc(x, f1, f2, f3);
                return;
            }
        }
    }

    post("list body");
    return;

    // error:
    //     pd_error(x, "psl_list error: incorrect arg type");
}

void psl_symbol(t_psl *x, t_symbol *s) {
    post("s: %s", s->s_name);

    // local buffer
    int length = strlen(s->s_name);
    char *buf = (char *)malloc(length * sizeof(char));
    strcpy(buf, s->s_name);
    post("buf: %s", buf);

    // clear expr_buffer
    memset(x->expr_buffer, 0, MAXPDSTRING);

    int j = 0;
    for (int i = 0; i < length; i++) {
        // remove escape `\` required for commas
        if (buf[i] != '\\') {
            x->expr_buffer[j++] = buf[i];
        } else if (x->expr_buffer[j - 1] == ' ') {
            j--;
        }
    }
    x->expr_buffer[length] = '\0';
    free(buf);

    post("x->expr_buffer: %s", x->expr_buffer);

    te_variable vars[] = {
        {"hypot", gsl_hypot, TE_FUNCTION2, NULL} /* TE_FUNCTION2 used because my_sum takes two arguments. */
    };

    te_expr *expr = te_compile(x->expr_buffer, vars, 2, 0);
    const double res = te_eval(expr);
    te_free(expr);
    outlet_float(x->out_f, res);
}

// message-methods

void psl_add(t_psl *x, t_floatarg f1, t_floatarg f2) {
    outlet_float(x->out_f, f1+f2);
}

void psl_rando(t_psl *x, t_floatarg n, t_floatarg seed) {
    post("rando: n:%.2f seed:%.2f", n, seed);
    gsl_rng_env_setup();

    int argc = (int)n;
    int _seed = (int)seed;

    gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(r, _seed);

    t_atom av[argc];
    for (int i = 0; i < argc; i++) {
        SETFLOAT(av + i, gsl_rng_uniform(r));
    }
    outlet_list(x->out_f, gensym("list"), argc, av);

    gsl_rng_free(r);
}

void psl_airy_ai(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_airy_Ai(f, GSL_PREC_APPROX));
}

void psl_airy_bi(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_airy_Bi(f, GSL_PREC_APPROX));
}

void psl_log1p(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_log1p(f));
}

void psl_expm1(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_expm1(f));
}

void psl_hypot(t_psl *x, t_floatarg f1, t_floatarg f2) {
    outlet_float(x->out_f, gsl_hypot(f1, f2));
}

void psl_hypot3(t_psl *x, t_floatarg f1, t_floatarg f2, t_floatarg f3) {
    outlet_float(x->out_f, gsl_hypot3(f1, f2, f3));
}

void psl_acosh(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_acosh(f));
}

void psl_asinh(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_asinh(f));
}

void psl_atanh(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_atanh(f));
}

void psl_ldexp(t_psl *x, t_floatarg f1, t_floatarg f2) {
    outlet_float(x->out_f, gsl_ldexp(f1, f2));
}

void psl_pow_int(t_psl *x, t_floatarg f1, t_floatarg f2) {
    outlet_float(x->out_f, gsl_pow_int(f1, f2));
}

void psl_pow_2(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_pow_2(f));
}

void psl_pow_3(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_pow_3(f));
}

void psl_pow_4(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_pow_4(f));
}

void psl_pow_5(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_pow_5(f));
}

void psl_pow_6(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_pow_6(f));
}

void psl_pow_7(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_pow_7(f));
}

void psl_pow_8(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_pow_8(f));
}

void psl_pow_9(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_pow_9(f));
}

void psl_fcmp(t_psl *x, t_floatarg f1, t_floatarg f2, t_floatarg f3) {
    outlet_float(x->out_f, gsl_fcmp(f1, f2, f3));
}

void psl_bessel_j0(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_bessel_J0(f));
}

void psl_bessel_j1(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_bessel_J1(f));
}

void psl_bessel_jn(t_psl *x, t_floatarg f1, t_floatarg f2) {
    outlet_float(x->out_f, gsl_sf_bessel_Jn(f1, f2));
}

void psl_bessel_y0(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_bessel_Y0(f));
}

void psl_bessel_y1(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_bessel_Y1(f));
}

void psl_bessel_yn(t_psl *x, t_floatarg f1, t_floatarg f2) {
    outlet_float(x->out_f, gsl_sf_bessel_Yn(f1, f2));
}

void psl_bessel_i0(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_bessel_I0(f));
}

void psl_bessel_i1(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_bessel_I1(f));
}

void psl_bessel_in(t_psl *x, t_floatarg f1, t_floatarg f2) {
    outlet_float(x->out_f, gsl_sf_bessel_In(f1, f2));
}

void psl_clausen(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_clausen(f));
}

void psl_dawson(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_dawson(f));
}

void psl_debye_1(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_debye_1(f));
}

void psl_debye_2(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_debye_2(f));
}

void psl_debye_3(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_debye_3(f));
}

void psl_debye_4(t_psl *x, t_floatarg f) {
    outlet_float(x->out_f, gsl_sf_debye_4(f));
}

// function selection
//---------------------------------------------------------------------------

// set default function from symbol
void select_default_function(t_psl *x, t_symbol *s) {
    x->func_name = s;
    post("func %s selected", s->s_name);

    switch (hash(s->s_name)) {
        case ADD:
            x->nargs = 2;
            x->bfunc = &psl_add;
            break;
        case LOG1P:
            x->nargs = 1;
            x->ufunc = &psl_log1p;
            break;
        case EXPM1:
            x->nargs = 1;
            x->ufunc = &psl_expm1;
            break;
        case HYPOT:
            x->nargs = 2;
            x->bfunc = &psl_hypot;
            break;
        case HYPOT3:
            x->nargs = 3;
            x->tfunc = &psl_hypot3;
            break;
        case ACOSH:
            x->nargs = 1;
            x->ufunc = &psl_acosh;
            break;
        case ASINH:
            x->nargs = 1;
            x->ufunc = &psl_asinh;
            break;
        case ATANH:
            x->nargs = 1;
            x->ufunc = &psl_atanh;
            break;
        case LDEXP:
            x->nargs = 2;
            x->bfunc = &psl_ldexp;
            break;
        case POW_INT:
            x->nargs = 2;
            x->bfunc = &psl_pow_int;
            break;
        case POW_2:
            x->nargs = 1;
            x->ufunc = &psl_pow_2;
            break;
        case POW_3:
            x->nargs = 1;
            x->ufunc = &psl_pow_3;
            break;
        case POW_4:
            x->nargs = 1;
            x->ufunc = &psl_pow_4;
            break;
        case POW_5:
            x->nargs = 1;
            x->ufunc = &psl_pow_5;
            break;
        case POW_6:
            x->nargs = 1;
            x->ufunc = &psl_pow_6;
            break;
        case POW_7:
            x->nargs = 1;
            x->ufunc = &psl_pow_7;
            break;
        case POW_8:
            x->nargs = 1;
            x->ufunc = &psl_pow_8;
            break;
        case POW_9:
            x->nargs = 1;
            x->ufunc = &psl_pow_9;
            break;
        case RANDO:
            x->nargs = 2;
            x->bfunc = &psl_rando;
            break;
        case FCMP:
            x->nargs = 3;
            x->tfunc = &psl_fcmp;
            break;
        case AIRY_AI:
            x->nargs = 1;
            x->ufunc = &psl_airy_ai;
            break;
        case AIRY_BI:
            x->nargs = 1;
            x->ufunc = &psl_airy_bi;
            break;
        case BESSEL_J0:
            x->nargs = 1;
            x->ufunc = &psl_bessel_j0;
            break;
        case BESSEL_J1:
            x->nargs = 1;
            x->ufunc = &psl_bessel_j1;
            break;
        case BESSEL_JN:
            x->nargs = 2;
            x->bfunc = &psl_bessel_jn;
            break;
        case BESSEL_Y0:
            x->nargs = 1;
            x->ufunc = &psl_bessel_y0;
            break;
        case BESSEL_Y1:
            x->nargs = 1;
            x->ufunc = &psl_bessel_y1;
            break;
        case BESSEL_YN:
            x->nargs = 2;
            x->bfunc = &psl_bessel_yn;
            break;
        case BESSEL_I0:
            x->nargs = 1;
            x->ufunc = &psl_bessel_i0;
            break;
        case BESSEL_I1:
            x->nargs = 1;
            x->ufunc = &psl_bessel_i1;
            break;
        case BESSEL_IN:
            x->nargs = 2;
            x->bfunc = &psl_bessel_in;
            break;
        case CLAUSEN:
            x->nargs = 1;
            x->ufunc = &psl_clausen;
            break;
        case DAWSON:
            x->nargs = 1;
            x->ufunc = &psl_dawson;
            break;
        case DEBYE_1:
            x->nargs = 1;
            x->ufunc = &psl_debye_1;
            break;
        case DEBYE_2:
            x->nargs = 1;
            x->ufunc = &psl_debye_2;
            break;
        case DEBYE_3:
            x->nargs = 1;
            x->ufunc = &psl_debye_3;
            break;
        case DEBYE_4:
            x->nargs = 1;
            x->ufunc = &psl_debye_4;
            break;
        default:
            post("func selection failed, reverting to defaults");
            break;
   }
}

// psl-inlet funcs
// ---------------------------------------------------------------------------

static void psl_inlet_float(t_psl_inlet *x, float f)
{
    x->owner->arg_array[x->id+1] = f;
    // outlet_float(x->owner->out_f, x->id + f);
    post("x->owner->arg_array[x->id]: %.02f", x->owner->arg_array[x->id+1]);
    psl_bang(x->owner);
}

// psl class constructor
// ---------------------------------------------------------------------------

void *psl_new(t_symbol *s) {
    t_psl *x = (t_psl *)pd_new(psl_class);

    // initialize variables
    x->nargs = 0;
    x->inlets = 0;
    x->ufunc = NULL;
    x->bfunc = NULL;
    x->tfunc = NULL;

    select_default_function(x, s);
    // sets x->nargs to correct number

    // create inlets
    x->inlets = x->nargs - 1;
    x->ins = (t_psl_inlet *)getbytes(x->inlets * sizeof(*x->ins));
    x->arg_array = malloc(x->inlets * sizeof(float));

    for (int i=0; i < x->inlets; i++) {
        x->ins[i].x_pd = psl_inlet_class;
        x->ins[i].owner = x;
        x->ins[i].id = i;
        x->arg_array[i] = 0.0;
        inlet_new((t_object *)x, &(x->ins[i].x_pd), 0, 0);
    }

    // initialize outlets
    x->out_f = outlet_new(&x->x_obj, &s_float);

    return (void *)x;
}

// psl class destructor
// ---------------------------------------------------------------------------

// TODO: not sure if this is correct!
void psl_free(t_psl *x) {
    free(x->arg_array);
    freebytes(x->ins, x->inlets * sizeof(*x->ins));
    post("DONE");
}

// psl class setup
// ---------------------------------------------------------------------------

void psl_setup(void) {

    psl_inlet_class = class_new(gensym("psl-inlet"),
                                0, 0,
                                sizeof(t_psl_inlet),
                                CLASS_PD,
                                0);

    if (psl_inlet_class) {
        class_addfloat(psl_inlet_class, (t_method)psl_inlet_float);
    }

    psl_class = class_new(gensym("psl"),
                        (t_newmethod)psl_new,
                        (t_method)psl_free,  // destructor
                        sizeof(t_psl),
                        CLASS_DEFAULT,
                        A_DEFSYMBOL,
                        0);

    // typed methods
    class_addbang(psl_class, psl_bang);
    class_addfloat(psl_class, psl_float);
    class_addlist(psl_class, psl_list);
    class_addsymbol(psl_class, psl_symbol);

    // message methods

    // class-addmethods
    class_addmethod(psl_class, (t_method)psl_add,  gensym("add"), A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_log1p,  gensym("log1p"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_expm1,  gensym("expm1"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_hypot,  gensym("hypot"), A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_hypot3,  gensym("hypot3"), A_DEFFLOAT, A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_acosh,  gensym("acosh"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_asinh,  gensym("asinh"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_atanh,  gensym("atanh"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_ldexp,  gensym("ldexp"), A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_pow_int,  gensym("pow_int"), A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_pow_2,  gensym("pow_2"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_pow_3,  gensym("pow_3"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_pow_4,  gensym("pow_4"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_pow_5,  gensym("pow_5"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_pow_6,  gensym("pow_6"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_pow_7,  gensym("pow_7"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_pow_8,  gensym("pow_8"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_pow_9,  gensym("pow_9"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_rando,  gensym("rando"), A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_fcmp,  gensym("fcmp"), A_DEFFLOAT, A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_airy_ai,  gensym("airy_ai"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_airy_bi,  gensym("airy_bi"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_bessel_j0,  gensym("bessel_j0"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_bessel_j1,  gensym("bessel_j1"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_bessel_jn,  gensym("bessel_jn"), A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_bessel_y0,  gensym("bessel_y0"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_bessel_y1,  gensym("bessel_y1"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_bessel_yn,  gensym("bessel_yn"), A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_bessel_i0,  gensym("bessel_i0"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_bessel_i1,  gensym("bessel_i1"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_bessel_in,  gensym("bessel_in"), A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_clausen,  gensym("clausen"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_dawson,  gensym("dawson"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_debye_1,  gensym("debye_1"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_debye_2,  gensym("debye_2"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_debye_3,  gensym("debye_3"), A_DEFFLOAT, 0);
    class_addmethod(psl_class, (t_method)psl_debye_4,  gensym("debye_4"), A_DEFFLOAT, 0);

    // create alias
    class_addcreator((t_newmethod)psl_new, gensym("gsl"), A_DEFSYMBOL, 0);

    // set name of default help file
    class_sethelpsymbol(psl_class, gensym("help-psl"));
}

max_stars_repo_name: shakfu/pd-psl

max_issues_repo_name: shakfu/pd-psl

max_forks_repo_name: shakfu/pd-psl


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

/*
 * This module calculates the solution of the normal equations least squares
 * system:
 *
 * [ J~^T J~ + mu D~^T D~ ] p~ = -J~^T f
 *
 * using the modified Cholesky decomposition. Quantities are scaled
 * according to:
 *
 * J~ = J S
 * D~ = D S
 * p~ = S^{-1} p
 *
 * where S is a diagonal matrix and S_jj = || J_j || and J_j is column
 * j of the Jacobian. This balancing transformation seems to be more
 * numerically stable for some Jacobians.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>

#include "common.c"

typedef struct
{
  gsl_matrix *JTJ;           /* J^T J */
  gsl_matrix *work_JTJ;      /* copy of J^T J */
  gsl_vector *rhs;           /* -J^T f, size p */
  gsl_permutation *perm;     /* permutation matrix for modified Cholesky */
  gsl_vector *work3p;        /* workspace, size 3*p */
  double mu;                 /* current regularization parameter */
} mcholesky_state_t;

static void *mcholesky_alloc (const size_t n, const size_t p);
static int mcholesky_init(const void * vtrust_state, void * vstate);
static int mcholesky_presolve(const double mu, const void * vtrust_state, void * vstate);
static int mcholesky_solve(const gsl_vector * f, gsl_vector *x,
                           const  void * vtrust_state, void *vstate);
static int mcholesky_solve_rhs(const gsl_vector * b, gsl_vector *x, mcholesky_state_t *state);
static int mcholesky_regularize(const double mu, const gsl_vector * diag, gsl_matrix * A,
                                mcholesky_state_t * state);

static void *
mcholesky_alloc (const size_t n, const size_t p)
{
  mcholesky_state_t *state;

  (void)n;

  state = calloc(1, sizeof(mcholesky_state_t));
  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate mcholesky state", GSL_ENOMEM);
    }

  state->JTJ = gsl_matrix_alloc(p, p);
  if (state->JTJ == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for JTJ", GSL_ENOMEM);
    }

  state->work_JTJ = gsl_matrix_alloc(p, p);
  if (state->work_JTJ == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for JTJ workspace",
                      GSL_ENOMEM);
    }

  state->rhs = gsl_vector_alloc(p);
  if (state->rhs == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for rhs", GSL_ENOMEM);
    }

  state->perm = gsl_permutation_alloc(p);
  if (state->perm == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for perm", GSL_ENOMEM);
    }

  state->work3p = gsl_vector_alloc(3 * p);
  if (state->work3p == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for work3p", GSL_ENOMEM);
    }

  state->mu = -1.0;

  return state;
}

static void
mcholesky_free(void *vstate)
{
  mcholesky_state_t *state = (mcholesky_state_t *) vstate;

  if (state->JTJ)
    gsl_matrix_free(state->JTJ);

  if (state->work_JTJ)
    gsl_matrix_free(state->work_JTJ);

  if (state->rhs)
    gsl_vector_free(state->rhs);

  if (state->perm)
    gsl_permutation_free(state->perm);

  if (state->work3p)
    gsl_vector_free(state->work3p);

  free(state);
}

static int
mcholesky_init(const void * vtrust_state, void * vstate)
{
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  mcholesky_state_t *state = (mcholesky_state_t *) vstate;

  /* compute J^T J */
  gsl_blas_dsyrk(CblasLower, CblasTrans, 1.0, trust_state->J, 0.0, state->JTJ);

  return GSL_SUCCESS;
}

/*
mcholesky_presolve()
  Compute the modified Cholesky decomposition of J^T J + mu D^T D.
Modified Cholesky is used in case mu = 0 and there are rounding
errors in forming J^T J which could lead to an indefinite matrix.

Inputs: mu     - LM parameter
        vstate - workspace

Notes:
1) On output, state->work_JTJ contains the Cholesky decomposition of
J^T J + mu D^T D
*/

static int
mcholesky_presolve(const double mu, const void * vtrust_state, void * vstate)
{
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  mcholesky_state_t *state = (mcholesky_state_t *) vstate;
  gsl_matrix *JTJ = state->work_JTJ;
  const gsl_vector *diag = trust_state->diag;
  int status;

  /* copy lower triangle of A to workspace */
  gsl_matrix_tricpy(CblasLower, CblasNonUnit, JTJ, state->JTJ);

  /* augment normal equations: A -> A + mu D^T D */
  status = mcholesky_regularize(mu, diag, JTJ, state);
  if (status)
    return status;

  /* compute modified Cholesky decomposition */
  status = gsl_linalg_mcholesky_decomp(JTJ, state->perm, NULL);
  if (status)
    return status;

  state->mu = mu;

  return GSL_SUCCESS;
}

/*
mcholesky_solve()
  Compute (J^T J + mu D^T D) x = -J^T f

Inputs: f      - right hand side vector f
        x      - (output) solution vector
        vstate - mcholesky workspace
*/

static int
mcholesky_solve(const gsl_vector * f, gsl_vector *x,
               const  void * vtrust_state, void *vstate)
{
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  mcholesky_state_t *state = (mcholesky_state_t *) vstate;
  int status;

  /* compute rhs = -J^T f */
  gsl_blas_dgemv(CblasTrans, -1.0, trust_state->J, f, 0.0, state->rhs);

  status = mcholesky_solve_rhs(state->rhs, x, state);
  if (status)
    return status;

  return GSL_SUCCESS;
}

static int
mcholesky_rcond(double * rcond, void * vstate)
{
  int status;
  mcholesky_state_t *state = (mcholesky_state_t *) vstate;
  double rcond_JTJ;

  if (state->mu != 0)
    {
      /*
       * Cholesky decomposition hasn't been computed yet, or was computed
       * with mu > 0 - recompute Cholesky decomposition of J^T J
       */

      /* copy lower triangle of JTJ to workspace */
      gsl_matrix_tricpy(CblasLower, CblasNonUnit, state->work_JTJ, state->JTJ);

      /* compute modified Cholesky decomposition */
      status = gsl_linalg_mcholesky_decomp(state->work_JTJ, state->perm, NULL);
      if (status)
        return status;
    }

  status = gsl_linalg_mcholesky_rcond(state->work_JTJ, state->perm, &rcond_JTJ, state->work3p);
  if (status == GSL_SUCCESS)
    *rcond = sqrt(rcond_JTJ);

  return status;
}

/* solve: (J^T J + mu D^T D) x = b */
static int
mcholesky_solve_rhs(const gsl_vector * b, gsl_vector *x, mcholesky_state_t *state)
{
  int status;
  gsl_matrix *JTJ = state->work_JTJ;

  status = gsl_linalg_mcholesky_solve(JTJ, state->perm, b, x);
  if (status)
    return status;

  return GSL_SUCCESS;
}

/* A <- A + mu D^T D */
static int
mcholesky_regularize(const double mu, const gsl_vector * diag, gsl_matrix * A,
                     mcholesky_state_t * state)
{
  (void) state;

  if (mu != 0.0)
    {
      size_t i;

      for (i = 0; i < diag->size; ++i)
        {
          double di = gsl_vector_get(diag, i);
          double *Aii = gsl_matrix_ptr(A, i, i);
          *Aii += mu * di * di;
        }
    }

  return GSL_SUCCESS;
}

static const gsl_multifit_nlinear_solver mcholesky_type =
{
  "mcholesky",
  mcholesky_alloc,
  mcholesky_init,
  mcholesky_presolve,
  mcholesky_solve,
  mcholesky_rcond,
  mcholesky_free
};

const gsl_multifit_nlinear_solver *gsl_multifit_nlinear_solver_mcholesky = &mcholesky_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#pragma once
#include <gsl/span>
#include "halley/core/api/audio_api.h"
#include "audio_buffer.h"

#if defined(_M_X64) || defined(__x86_64__)
#define HAS_SSE
#if !defined(__linux__)
#define HAS_AVX
#endif
#endif

#if defined(_M_IX86) || defined(__i386)
// Might not be available, but do we really care about such old processors?
#define HAS_SSE
#endif

namespace Halley
{
	class AudioMixer
	{
	public:
		virtual ~AudioMixer() {}

		virtual void mixAudio(gsl::span<const AudioSamplePack> src, gsl::span<AudioSamplePack> dst, float gainStart, float gainEnd);
		virtual void interleaveChannels(gsl::span<AudioSamplePack> dst, gsl::span<AudioBuffer*> src);
		virtual void compressRange(gsl::span<AudioSamplePack> buffer);
		static std::unique_ptr<AudioMixer> makeMixer();
	};
}

max_stars_repo_name: lye/halley

max_issues_repo_name: lye/halley

max_forks_repo_name: lye/halley


text:
#ifndef __asl234n0a2__LIBRARIES_H
  #define __asl234n0a2__LIBRARIES_H

  /* Standard libraries : */
  #include <stdio.h>
  #include <stdlib.h>
  #include <signal.h>
  #include <time.h>
  #include <math.h>
    /* wtf ? */
    #include <fcntl.h>
    #include <unistd.h>
    #include <limits.h>

  /* Sleeping-related */
  #include <time.h>
  #include <errno.h>

  /* GNU Scientific library */
  #include <gsl/gsl_rng.h>

  /* Simple Direct Media layer (SDL) graphics: */
  #include <SDL/SDL.h>

  /* Locally defined */
    /* Headers : */
    #include "dynamics.h"
    #include "simulation.h"
    #include "display.h"
    #include "parsing.h"
    #include "aux_math.h"
    #include "dev_random.h"
    #include "mersenne_twister.h"
    #include "singly_linked_list.h"
    #include "datastructures.h"
    #include "SDL_datastructure.h"
    #include "visualization.h"

    /* Constants and macros : */
    #include "constants.h"
    #include "probabilities.h"

#endif

max_stars_repo_name: gmagannaDevelop/TousAntiCovid

max_issues_repo_name: gmagannaDevelop/TousAntiCovid

max_forks_repo_name: gmagannaDevelop/TousAntiCovid


text:
/**
 * Provides function to fit a polynomial curve to a set of data.
 * Source: https://rosettacode.org/wiki/Polynomial_regression#C
 */

#ifndef _POLIFITGSL_H
#define _POLIFITGSL_H
#include <gsl/gsl_multifit.h>
#include <stdbool.h>
#include <math.h>
bool polynomialfit(int obs, int degree,
        const double *dx, const double *dy, double *store); /* n, p */
#endif


text:
/* monte/gsl_monte_plain.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
 *
 * 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.
 */

/* Plain Monte-Carlo. */

/* Author: MJB */

#ifndef __GSL_MONTE_PLAIN_H__
#define __GSL_MONTE_PLAIN_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdio.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_rng.h>

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

__BEGIN_DECLS

typedef struct {
  size_t dim;
  double *x;
} gsl_monte_plain_state;

GSL_FUN int
gsl_monte_plain_integrate (const gsl_monte_function * f,
                           const double xl[], const double xu[],
                           const size_t dim,
                           const size_t calls,
                           gsl_rng * r,
                           gsl_monte_plain_state * state,
                           double *result, double *abserr);

GSL_FUN gsl_monte_plain_state* gsl_monte_plain_alloc(size_t dim);

GSL_FUN int gsl_monte_plain_init(gsl_monte_plain_state* state);

GSL_FUN void gsl_monte_plain_free (gsl_monte_plain_state* state);

__END_DECLS

#endif /* __GSL_MONTE_PLAIN_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
/*
** Test for irreducibility of a symmetric matrix R = X^T X,
** i.e. there is at least one positive off-diagonal entry in each row and column.
** Since R is symmetric, it suffices to show that each row has one non-negative entry.
** i.e. (X^TX v - v)_i > 0 for all i=0,...,n-1
** exclude main diagonal
**
** G.Lohmann, MPI-KYB, Nov 2018
*/

#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>

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

extern void SignFlip(gsl_matrix_float *X);

/*
** if irreducible return 1, else return 0
*/
int Irreducible(gsl_matrix_float *X,int type)
{
  size_t i;
  size_t nvox = X->size1;
  size_t nt = X->size2;
  float u=0,tiny=1.0e-6;

  /* ini */
  gsl_vector_float *x = gsl_vector_float_calloc(nvox);
  gsl_vector_float *y = gsl_vector_float_calloc(nvox);
  gsl_vector_float *z = gsl_vector_float_calloc(nt);
  for (i=0; i<nvox; i++) x->data[i] = y->data[i] = 1.0;

  /* R = X^t X */
  gsl_blas_sgemv (CblasTrans,1.0,X,x,0.0,z);
  if (type == 7) SignFlip(X);
  gsl_blas_sgemv (CblasNoTrans,1.0,X,z,0.0,y);
  if (type == 7) SignFlip(X);

  /* count off-diagonal entries with value '0' */
  size_t n=0;
  for (i=0; i<nvox; i++) {
    u = y->data[i] - x->data[i];  /* subtract main diagonal */
    if (fabs(u) < tiny) n++;
  }
  gsl_vector_float_free(x);
  gsl_vector_float_free(y);
  gsl_vector_float_free(z);

  /* return code */
  int rtcode = 0;
  if (n > 0) rtcode = 0;
  else rtcode = 1;
  if (rtcode == 0) {
    VWarning(" Matrix is not irreducible,   n=%lu,  nvox= %lu\n",n,nvox);
  }
  return rtcode;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#ifndef CONSTANTS_H
#define CONSTANTS_H

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <limits.h>
#include <math.h>
#include <time.h>
#include <float.h>
#include <utils/version.h>
#include <utils/help.h>
#include <utils/error.h>
#include <gsl/gsl_statistics_double.h>

/*
 * MIN macro
 */
#define MIN(a,b) (((a)<(b))?(a):(b))

/*
 * MAX macro
 */
#define MAX(a,b) (((a)>(b))?(a):(b))

/*
 * STR macro
 */
#define STR(a) ((a > 0)?("-"):("+"))

/*
 * Default value for read_minlen parameter
 */
#define MIN_READ_LEN 0

/*
 * Default value for minlen parameter
 */
#define MIN_LEN 15

/*
 * Default value for maxlen parameter
 */
#define MAX_LEN 30

/*
 * Default value for spacing parameter
 */
#define SPACING 20

/*
 * Default value for minheight parameter
 */
#define MIN_READS 10.0f

/*
 * Default value for trimming threshold parameter
 */
#define TRIM_THRESHOLD 0.05

/*
 * Default value for minimum trimming parameter
 */
#define TRIM_MIN 2

/*
 * Default value for maximum trimming parameter
 */
#define TRIM_MAX 20

/*
 * Default value for IDR cutoff
 */
#define CUTOFF 2.0f

/*
 * Default value for Replicate number
 */
#define REPLICATE_NUMBER 1

/*
 * Replicate treatment options
 */
#define REPLICATE_POOL_STR "pool" // default value
#define REPLICATE_MEAN_STR "mean"
#define REPLICATE_REPLICATE_STR "replicate"
#define REPLICATE_POOL 0
#define REPLICATE_MEAN 1
#define REPLICATE_REPLICATE 2

/*
 * IDR method options
 */
#define IDR_COMMON_STR "common" // default value
#define IDR_NONE_STR "none"
#define IDR_SERE_STR "sere"
#define IDR_IDR_STR "idr"
#define IDR_COMMON 0
#define IDR_NONE 1
#define IDR_SERE 2
#define IDR_IDR 3

/*
 * Maximum length of a contig
 */
#define MAX_CONTIG_LENGTH 200

/*
 * Maximum number of contigs
 */
#define MAX_CONTIGS 10000000

/*
 * Maximum size, in chars, of an error message
 */
#define MAX_ERR_MSG 100

/*
 * Maximum size, in chars, of a given path
 */
#define MAX_PATH 500

/*
 * Maximum size, in chars, of the 4th field (name) in the query BED file
 */
#define MAX_FEATURE 50

/*
 * Maximum number of replicates
 */
#define MAX_REPLICATES 10

/*
 * Maximum number of alignments per heap
 */
#define MAX_ALIGN_HEAP 50000000

/*
 * Maximum read length
 */
#define MAX_READ_LENGTH 200

/*
 * Maximum profile length
 */
#define MAX_PROFILE_LENGTH 500

/*
 * Maximum number of block base pairs
 */
#define MAX_BLOCK 3000

/*
 * Alignment strand
 */
#define FWD_STRAND 0 // forward/watson
#define REV_STRAND 1 // reverse/crick

/*
 * Alignment validity
 */
#define VALID_ALIGNMENT 1
#define INVALID_ALIGNMENT 0

/*
 * Constants for npIDR method
 */
#define ABSOLUTE 0
#define CONDITIONAL 1

/*
 * Path separator
 */
#ifdef __unix__
  #define PATH_SEPARATOR "/"
#else
  #define PATH_SEPARATOR "\\"
#endif

/*
 * Output file suffixes
 */
#define PROFILES_SUFFIX "profiles.dat"
#define CONTIGS_SUFFIX "contigs.dat"
#define CROSSCOR_SUFFIX "crosscor.dat"
#define CLUSTERS_SUFFIX "clusters.neWick"
#define ANNOTATION_O_SUFFIX "annotation.bed"
#define TMPROFILES_SUFFIX "tmprofiles.dat"

/*
 * Maximum N limit for gaussian white noise generation
 */
#define MAX_GNOISE_N 20

/*
 * Cluster cutoff default value
 */
#define CLUSTER_CUTOFF -1.0f

/*
 * Condition for existence of annotation file
 */
#define ANNOTATION_CONDITION 0

/*
 * Condition for existence of additional profiles file
 */
#define ADDITIONAL_P_CONDITION 0

/*
 * Default value for the feature to profile overlap percentage
 */
#define OVERLAP_FTOP 0.9

/*
 * Default value for the profile to feature overlap percentage
 */
#define OVERLAP_PTOF 0.5

/*
 * Maximum number of annotation files
 */
#define MAX_ANNOTATIONS 10

/*
 * Condition for existence of correlations file
 */
#define CORRELATIONS_CONDITION 0

/*
 * Constants for profile category
 */
#define NOVEL 0
#define KNOWN 1

/*
 * Differential processing default p-value
 */
#define P_VALUE 0.01

/*
 * Differential processing default overlap
 */
#define DP_FOLD_CHANGE 7

/*
 * Suffix for differentially processed profiles
 */
#define DIFFPROC_PROFILE_O_SUFFIX "diffprofiles.dat"
/*
 * Suffix for differentially processed clusters
 */
#define DIFFPROC_CLUSTER_O_SUFFIX "diffclusters.dat"
#endif

max_stars_repo_name: comprna/SeRPeNT

max_issues_repo_name: comprna/SeRPeNT

max_forks_repo_name: comprna/SeRPeNT


text:
/*
 *  aplowly.h
 *  lowly_project
 *
 *  Created by Karim Benabed on 11/03/09.
 *  Copyright 2009 Institut d'Astrophysique de Paris. All rights reserved.
 *
 */

#include "lowly_common.h"

#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_int.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_gamma.h>

#ifndef __APLOWLY__
#define __APLOWLY__

typedef struct {
  double *data;
  double *alpha,*beta,*wl,*C,*gl;
  double Nl,logdets2;
  int is_zero;
  size_t ndim;
  gsl_matrix_view  C_view_container;
  gsl_matrix  *C_view;
  gsl_vector_view  gl_view_container;
  gsl_vector *gl_view;

} approx_lowly;

approx_lowly* approx_lowly_init(size_t ndim,int* ell, int lmax,double *clhat, double *fsky, double *wl, double *V,double noisevar, long nside, error **err);
void free_approx_lowly(approx_lowly **self);
double approx_lowly_log_pdf(void *self, double *Cl, error **err);
double approx_zero_lowly_log_pdf(void* pelf,double* pars, error **err);

typedef struct {
  double *data;
  double *al,*bl,*Cm,*gl,*nl;
  int *ell;
  double logdets2;
  int isdiag;
  int nell;
} tease;

typedef struct {
  double *data;
  double *cl,*vl,*nl,*Cm,*gl;
  int *ell;
  double logdets2;
  int isdiag;
  int nell;
} gausslkl;

tease* tease_init(size_t nell, int* ell,double *al, double *bl, double *nl, int nlcst, double *Cm, int Cmisdiag, error **err);
void tease_free(void **ping);
double tease_log_pdf(void* ing, double* pars, error **err);

gausslkl* gausslkl_init(size_t nell, int* ell,double *al, double *bl, double *nl, int nlcst, double *Cm, int Cmisdiag, error **err);
void gausslkl_free(void **ping);
double gausslkl_log_pdf(void* ing, double* pars, error **err);

typedef struct {
  double *data;
  double *ldlh,*sig2,*wl;
  double Nl;
  size_t ndim;
} log_normal;

log_normal* log_normal_init(size_t ndim,int* ell, double *clhat, double *fsky, double *wl, double noisevar, long nside, error **err);
void free_log_normal(log_normal **self);
double log_normal_log_pdf(void *self, double *Cl, error **err);

typedef struct {
  int np;
  int *nb;
  int *obe;
  int *obv;
  double *be;
  double *bv;
} binlkl;

binlkl* init_binlkl(int np, int* nb, double* be, double* bv, error **err);
double binlkl_lkl(void *vbon, double *pars, error **err);
void free_binlkl(void** pvbon);

#define ZEROLKL -1e30

#define lowly_no_correlation -20 + clowly_base

#endif


text:
#include <stdio.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_pow_int.h>
#include <math.h>

double sin_term(int l,int m,double phi,double theta){
 double out = gsl_pow_int(-1.0,m)*gsl_sf_legendre_sphPlm(l,m,gsl_sf_cos(theta))*gsl_sf_sin(m*phi);
 return out;
}

double cos_term(int l,int m, double phi, double theta){
  double out= gsl_pow_int(-1.0,m)*gsl_sf_legendre_sphPlm(l,m,gsl_sf_cos(theta))*gsl_sf_cos(m*phi);
  return out;
}

float acos_fast(float x) {
  float negate = (float)x < 0;
  x =  x < 0 ? -x : x;
  float ret = -0.0187293;
  ret = ret * x;
  ret = ret + 0.0742610;
  ret = ret * x;
  ret = ret - 0.2121144;
  ret = ret * x;
  ret = ret + 1.5707288;
  ret = ret * sqrt(1.0-x);
  ret = ret - 2 * negate * ret;
  return negate * 3.14159265358979 + ret;
}

float fast_grf(int L, float x, float y,float z,float rands[]){
  float result = 0.0;
  float phi =atan2(y,x);
  float theta = acos_fast(z);
  float temp ;
  register int l = 0;
  register int m = 1;
  for (l=0; l<= L;l++){
    result += gsl_sf_legendre_sphPlm(l,0,gsl_sf_cos(theta))*rands[2*L*l];
    temp = 0.0;
    for(m=1;m<=l;m++){
      temp += (1 - ((m & 1) << 1))*gsl_sf_legendre_sphPlm(l,m,gsl_sf_cos(theta))*(rands[2*m+2*L*l+l]*gsl_sf_sin(m*phi)+rands[1+2*m+2*L*l]*gsl_sf_cos(m*phi));
    }
    result += 1.41421*temp;
  }
  return result;
}


text:
#pragma once

#include "arcana/finally_scope.h"
#include "arcana/functional/inplace_function.h"
#include "arcana/sentry.h"
#include "arcana/threading/affinity.h"

#include <algorithm>
#include <tuple>
#include <vector>

#include <gsl/gsl>

namespace arcana
{
    using ticket_seed = int64_t;

    using ticket = gsl::final_action<
        stdext::inplace_function<void(), sizeof(std::aligned_storage_t<sizeof(ticket_seed) + sizeof(void*)>)>>;

    using ticket_scope = finally_scope<ticket>;

    /*
        An event routing class used to dispatch events to multiple listeners.
        Each router class can only handle a certain fixed set of event types defined by EventTs.
    */
    template<typename... EventTs>
    class router
    {
    public:
        static constexpr size_t LISTENER_SIZE = 4 * sizeof(int64_t);

        template<typename EventT>
        using listener_function = stdext::inplace_function<void(const EventT&), LISTENER_SIZE>;

        /*
            Sends an event synchronously to all listeners.
        */
        template<typename EventT>
        void fire(const EventT& evt)
        {
            GSL_CONTRACT_CHECK("thread affinity", m_affinity.check());

            using event = std::decay_t<EventT>;

            auto& listeners = std::get<listener_group<event>>(m_listeners);

            {
                auto guard = std::get<sentry<event>>(m_sentries).take();

                for (listener<event>& listener : listeners)
                {
                    if (listener.valid)
                    {
                        listener.callback(evt);
                    }
                }
            }

            // if we're no longer iterating the listeners list in this stack
            // remove all the unregistered listeners and add the pending ones
            if (!std::get<sentry<event>>(m_sentries).is_active())
            {
                listeners.erase(std::remove_if(listeners.begin(),
                                               listeners.end(),
                                               [](const listener<event>& l) { return !l.valid; }),
                                listeners.end());

                // move the pending listeners to the real list
                // and clear the pending list
                auto& pending = std::get<listener_group<event>>(m_pending);
                std::move(pending.begin(), pending.end(), std::back_inserter(listeners));
                pending.clear();
            }
        }

        /*
            Adds an event listener.
        */
        template<typename EventT, typename T>
        ticket add_listener(T&& listener)
        {
            auto id = internal_add_listener<EventT>(std::forward<T>(listener));

            return ticket{ [id, this] { internal_remove_listener<EventT>(id); } };
        }

        /*
            Sets the routers thread affinity. Once this is set the methods
            on this instance will need to be called by that thread.
        */
        void set_affinity(const affinity& aff)
        {
            m_affinity = aff;
        }

    private:
        /*
        Adds an event listener.
        */
        template<typename EventT, typename T>
        ticket_seed internal_add_listener(T&& listener)
        {
            GSL_CONTRACT_CHECK("thread affinity", m_affinity.check());

            using event = std::decay_t<EventT>;

            // if we're currently firing an event in that group we need to wait until we're done before
            // adding the listener to the list
            auto& listeners = std::get<sentry<event>>(m_sentries).is_active()
                                  ? std::get<listener_group<event>>(m_pending)
                                  : std::get<listener_group<event>>(m_listeners);

            auto id = m_nextId++;
            listeners.emplace_back(std::forward<T>(listener), id);
            return id;
        }

        /*
        Removes an event listener by id.
        */
        template<typename EventT>
        void internal_remove_listener(const ticket_seed& id)
        {
            GSL_CONTRACT_CHECK("thread affinity", m_affinity.check());

            using event = std::decay_t<EventT>;

            auto& listeners = std::get<listener_group<event>>(m_listeners);
            auto found = std::find_if(listeners.begin(), listeners.end(), [id](const listener<event>& listener) {
                return listener.id == id;
            });

            if (found == listeners.end())
            {
                assert(false && "removing item that isn't there");
                return;
            }

            // don't modify the collection while iterating, just disable the listener
            if (std::get<sentry<event>>(m_sentries).is_active())
            {
                found->valid = false;
            }
            else
            {
                listeners.erase(found);
            }
        }

        template<typename EventT>
        struct listener
        {
            using callback_t = listener_function<EventT>;

            callback_t callback;
            ticket_seed id;
            bool valid;

            listener(callback_t&& callback, const ticket_seed& id)
                : callback{ std::move(callback) }
                , id{ id }
                , valid{ true }
            {}

            listener(const callback_t& callback, const ticket_seed& id)
                : callback{ callback }
                , id{ id }
                , valid{ true }
            {}
        };

        template<typename EventT>
        using listener_group = std::vector<listener<EventT>>;

        std::tuple<listener_group<EventTs>...> m_listeners;
        std::tuple<sentry<EventTs>...> m_sentries;
        std::tuple<listener_group<EventTs>...> m_pending;

        affinity m_affinity;
        ticket_seed m_nextId;
    };
}