text:
#ifndef __REGIONAL_H_
#define __REGIONAL_H_

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

typedef struct stochasticity_grid {
    double *arr;
    int height, width;
    double weight_factor;
    double stddev;
    double mean;
    int depth;
    gsl_rng *rng;
    double *weights;
} t_stochasticity_grid;

t_stochasticity_grid *create_stochasticity_grid(double *arr, int height, int width, double weight_factor, double mean, double stddev, long int seed);
void free_stochasticity_grid(t_stochasticity_grid *g);
void generate_stochasticity(t_stochasticity_grid *g);

#endif /* __REGIONAL_H_ */

max_stars_repo_name: rybicki/corridor-spom

max_issues_repo_name: rybicki/corridor-spom

max_forks_repo_name: rybicki/corridor-spom


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

/***************************************************************************//**
 *
 **/
void QUARK_CORE_cgemv(Quark *quark, Quark_Task_Flags *task_flags,
                      PLASMA_enum trans, int m, int n,
                      PLASMA_Complex32_t alpha, const PLASMA_Complex32_t *A, int lda,
                                                const PLASMA_Complex32_t *x, int incx,
                      PLASMA_Complex32_t beta,        PLASMA_Complex32_t *y, int incy)
{
    DAG_CORE_GEMV;
    QUARK_Insert_Task(quark, CORE_cgemv_quark, task_flags,
        sizeof(PLASMA_enum),             &trans,  VALUE,
        sizeof(int),                     &m,      VALUE,
        sizeof(int),                     &n,      VALUE,
        sizeof(PLASMA_Complex32_t),      &alpha,  VALUE,
        sizeof(PLASMA_Complex32_t)*m*n,  A,               INPUT,
        sizeof(int),                     &lda,    VALUE,
        sizeof(PLASMA_Complex32_t)*n,    x,               INPUT,
        sizeof(int),                     &incx,   VALUE,
        sizeof(PLASMA_Complex32_t),      &beta,   VALUE,
        sizeof(PLASMA_Complex32_t)*m,    y,               INOUT,
        sizeof(int),                     &incy,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cgemv_quark = PCORE_cgemv_quark
#define CORE_cgemv_quark PCORE_cgemv_quark
#endif
void CORE_cgemv_quark(Quark *quark)
{
    PLASMA_enum trans;
    int m, n, lda, incx, incy;
    PLASMA_Complex32_t alpha, beta;
    const PLASMA_Complex32_t *A, *x;
    PLASMA_Complex32_t *y;

    quark_unpack_args_11( quark, trans, m, n, alpha, A, lda, x, incx, beta, y, incy );
    cblas_cgemv(
        CblasColMajor,
        (CBLAS_TRANSPOSE)trans,
        m, n,
        CBLAS_SADDR(alpha), A, lda,
                            x, incx,
        CBLAS_SADDR(beta),  y, incy);
}

ext: c


text:
#ifndef __GSLSTREAM_H__
#define __GSLSTREAM_H__

#include "RngStream.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

namespace rng {
  class GSLStream : public RngStream {
    public:
      GSLStream() : type(gsl_rng_taus2) {}
      virtual ~GSLStream() { free(); }

      inline void alloc(unsigned long seed = time(NULL)) {
        if (! is_alloc) {
          rng = gsl_rng_alloc(type);
          is_alloc = 1;
          gsl_rng_set(rng,seed);
        }
      }

      inline void free() {
        if (is_alloc) {
          gsl_rng_free(rng);
          is_alloc = 0;
        }
      }

      inline void uniform(size_t n, double* r, double a = 0.0, double b = 1.0) {
        for (size_t i = 0; i < n; ++i) {
          r[i] = a + (b-a)*gsl_rng_uniform(rng);
        }
      }

      inline void uniform_int(size_t n, int* r, int a = 0, int b = 10) {
        for (size_t i = 0; i < n; ++i) {
          r[i] = a + gsl_rng_uniform_int(rng,b-a);
        }
      }

      inline void shuffle(int* x, size_t n) { gsl_ran_shuffle(rng,x,n,sizeof(int)); }
      inline void shuffle(double* x, size_t n) { gsl_ran_shuffle(rng,x,n,sizeof(double)); }

      inline void multinomial(size_t k, size_t n, const double* p, unsigned* a) {
        gsl_ran_multinomial(rng,k,n,p,a);
      }

      inline void gaussian(size_t n, double* r, double mu = 0.0, double sigma = 1.0) {
        *r = gsl_ran_gaussian(rng,sigma) + mu;
      };

      inline void poisson(size_t n, int* k, double lambda) {
        *k = gsl_ran_poisson(rng,lambda);
      }

    protected:
      const gsl_rng_type* type;
      gsl_rng* rng;
  };
}

#endif

ext: h

max_stars_repo_name: gaberoo/cdream

max_issues_repo_name: gaberoo/cdream

max_forks_repo_name: gaberoo/cdream


text:
/*
 findcluster_ternary_par.c
 written by Eunseok Lee

 function: find the most representative cluster function

 v1: Feb 2, 2018
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/stat.h>
//#include <cem.h>
#include <mpi.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

void randperm(int*, int);
void load_double_mat(double*, int, int, int, char*);
void load_int_mat(int*, int, int, int, char*);
void load_double_array(double*, int, char*);
void load_int_array(int*, int, char*);
void load_double_value(int,char);
void load_int_value(int,char);
int mat_size_row(int*);
int find_value_in_array(int*,int,int,int);
void mat2d_prod(double*, int, int, double*, int, int, double*, int, int);
double obtainerr2_par(int,int,double*, int, int, int*, int, double*);
void least_square_solver(double*, int, int, double*, double*);

#define MASTER 0
#define FROM_MASTER 1
#define FROM_WORKER 2

int main(int argc, char **argv)
{
    int n1, n2, n3;
    int np, nfu;
    double kv[3][3], pk[3][3], pp[8][3];
    double *rp, *rpL, *rpN, *rpO;
    int *map_to_cluster1, *map_to_cluster2, *map_to_cluster3;
    int *nlist;
    int neighbor_num;
    int ncluster, ncluster2, ncluster3, ncorr_col;
    int c2start, c3start;
    int data;
    double *E, *Ef;
    int ndata;
    int i, j, k, nj, nk, ctn;
    int iter, max_iter;
    double kT=0.0256;   // Default value of kT
    int howmanycluster;
    int dispfreq;
    char buff_line[200], dummy[200], param_name[100];
    int param_name_len;
    double cvs_tol = 0.01;  // Default value of cvs_tol
    int construct_corr_mat = 0; //Default value of construct_corr_mat

    // mpi parameters and initialization
    int numprocs, rank, mtype;
    int row_dist_size;
    int row_ini, row_end, row_offset;
    MPI_Status status;
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);  // Get # processors
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);      // Get my rank (id)

    if (numprocs>1)
        printf("Parallel computing: the input files are read by all nodes. Multiple messages are not due to error!\n",numprocs);

    // Construct conrrelation matrix. It can also be loaded.

    int n_corr_mat_ug_row;
    int n_non_singular_col;
    FILE *fp;
    char dirname[100]="dir_inputs";
    char paramfilename[100];
    sprintf(paramfilename,"%s/param.dat",dirname);
    char corr_mat_ugs_filename[200];
    fp = fopen(paramfilename, "r");
    while(fgets(buff_line,sizeof(buff_line),fp) != NULL) {
        if (buff_line[0] == '#') {
            //            printf("Comment line: %s",buff_line);
            continue;
        }
        strcpy(param_name,"corr_mat_ugs_filename");
        param_name_len = strlen(param_name);
        strncpy(dummy,buff_line,param_name_len);
        dummy[param_name_len] = '\0';
        if (strcmp(dummy,param_name)==0) {
            sscanf(buff_line,"%s %s",dummy,corr_mat_ugs_filename);
            if (rank == MASTER)
                printf("corr_mat_ugs_filename = %s\n",corr_mat_ugs_filename);
        }
        strcpy(param_name,"n_corr_mat_ug_row");
        param_name_len = strlen(param_name);
        strncpy(dummy,buff_line,param_name_len);
        dummy[param_name_len] = '\0';
        if (strcmp(dummy,param_name)==0) {
            sscanf(buff_line,"%s %d",dummy,&n_corr_mat_ug_row);
            if (rank == MASTER)
                printf("n_corr_mat_ug_row = %d\n",n_corr_mat_ug_row);
        }
        strcpy(param_name,"n_non_singular_col");
        param_name_len = strlen(param_name);
        strncpy(dummy,buff_line,param_name_len);
        dummy[param_name_len] = '\0';
        if (strcmp(dummy,param_name)==0) {
            sscanf(buff_line,"%s %d",dummy,&n_non_singular_col);
            if (rank == MASTER)
                printf("n_non_singular_col = %d\n",n_non_singular_col);
        }
        strcpy(param_name,"howmanycluster");
        param_name_len = strlen(param_name);
        strncpy(dummy,buff_line,param_name_len);
        dummy[param_name_len] = '\0';
        if (strcmp(dummy,param_name)==0) {
            sscanf(buff_line,"%s %d",dummy,&howmanycluster);
            if (rank == MASTER)
                printf("howmanycluster = %d\n",howmanycluster);
        }
        strcpy(param_name,"max_iter");
        param_name_len = strlen(param_name);
        strncpy(dummy,buff_line,param_name_len);
        dummy[param_name_len] = '\0';
        if (strcmp(dummy,param_name)==0) {
            sscanf(buff_line,"%s %d",dummy,&max_iter);
            if (rank == MASTER)
                printf("max_iter = %d\n",max_iter);
        }
        strcpy(param_name,"kT");
        param_name_len = strlen(param_name);
        strncpy(dummy,buff_line,param_name_len);
        dummy[param_name_len] = '\0';
        if (strcmp(dummy,param_name)==0) {
            sscanf(buff_line,"%s %lf",dummy,&kT);
            if (rank == MASTER)
                printf("kT = %f\n",kT);
        }
        strcpy(param_name,"cvs_tol");
        param_name_len = strlen(param_name);
        strncpy(dummy,buff_line,param_name_len);
        dummy[param_name_len] = '\0';
        if (strcmp(dummy,param_name)==0) {
            sscanf(buff_line,"%s %lf",dummy,&cvs_tol);
            if (rank == MASTER)
                printf("cvs_tol = %f\n",cvs_tol);
        }
        strcpy(param_name,"dispfreq");
        param_name_len = strlen(param_name);
        strncpy(dummy,buff_line,param_name_len);
        dummy[param_name_len] = '\0';
        if (strcmp(dummy,param_name)==0) {
            sscanf(buff_line,"%s %d",dummy,&dispfreq);
            if (rank == MASTER)
                printf("dispfreq = %d\n",dispfreq);
        }
    }
    fclose(fp);

    char loadfilename[100];
    double *corr_mat_ugs;
    corr_mat_ugs = (double*) malloc(n_corr_mat_ug_row*n_non_singular_col*sizeof(double));
    sprintf(loadfilename,"%s/%s",dirname,corr_mat_ugs_filename);
    load_double_mat(corr_mat_ugs,n_corr_mat_ug_row,n_non_singular_col,1,loadfilename);

    int usefulcorr_col[n_non_singular_col];
    sprintf(loadfilename,"%s/usefulcorr_col.dat",dirname);
    load_int_array(usefulcorr_col,n_non_singular_col,loadfilename);

    double Ef_ug[n_corr_mat_ug_row];
    sprintf(loadfilename,"%s/Ef_ug.dat",dirname);
    load_double_array(Ef_ug,n_corr_mat_ug_row,loadfilename);

/*    int nC_ug[n_corr_mat_ug_row];
    sprintf(loadfilename,"%s/nC_ug.dat",dirname);
    load_int_array(nC_ug,n_corr_mat_ug_row,loadfilename);
*/
    row_offset = (int) ceil(n_corr_mat_ug_row*1.0/numprocs);
    if (numprocs>1 && rank==MASTER)
        printf("Correlation Matrix is calculated over %d processors: row_offset = %d\n",numprocs,row_offset);
    double err2_sum_from_master;
    double err2_sum_from_worker;
    int cluster_set1[howmanycluster];

    /**** Start of MC simulation ****/
    if (rank == MASTER) {
        int cluster_set1_old[howmanycluster], cluster_set1_min[howmanycluster], cluster_set2[n_non_singular_col-howmanycluster];
        double cvs, cvs_old, cvs_min;
        int select1, select2, target, candidate;

        for (i=0;i<howmanycluster;i++)
            cluster_set1[i] = i;
        for (i=0;i<n_non_singular_col-howmanycluster;i++)
            cluster_set2[i] = i+howmanycluster;

        mtype = FROM_MASTER;
        for (i=1;i<numprocs;i++) {
            MPI_Send(&cluster_set1, howmanycluster, MPI_INT, i, mtype, MPI_COMM_WORLD);
        }
        row_ini = rank * row_offset;
        row_end = row_ini +row_offset;
        err2_sum_from_master = obtainerr2_par(row_ini,row_end,corr_mat_ugs,n_corr_mat_ug_row,n_non_singular_col,cluster_set1,howmanycluster,Ef_ug);
//        printf("initial err2 = %f\n",err2_sum_from_master);
        for (i=1;i<numprocs;i++) {
            mtype = FROM_WORKER;
            MPI_Recv(&err2_sum_from_worker, 1, MPI_DOUBLE, i, mtype, MPI_COMM_WORLD, &status);
            err2_sum_from_master += err2_sum_from_worker;
        }
        printf("initial err2 = %f\n",err2_sum_from_master);
        cvs = sqrt(err2_sum_from_master/n_corr_mat_ug_row);
        cvs_min = cvs;
        for (i=0;i<howmanycluster;i++)
            cluster_set1_min[i] = cluster_set1[i];
        printf("initial cvs = %f\n",cvs);
        for (iter=0;iter<max_iter;iter++) {
            select1 = (int) floor(drand48()*howmanycluster);
            select2 = (int) floor(drand48()*(n_non_singular_col-howmanycluster));
            target = cluster_set1[select1];
            candidate = cluster_set2[select2];
            cvs_old = cvs;
            cluster_set1[select1] = candidate;
            mtype = FROM_MASTER;
            for (i=1;i<numprocs;i++) {
                MPI_Send(&cluster_set1, howmanycluster, MPI_INT, i, mtype, MPI_COMM_WORLD);
            }
            row_ini = rank * row_offset;
            row_end = row_ini +row_offset;
            err2_sum_from_master = obtainerr2_par(row_ini,row_end,corr_mat_ugs,n_corr_mat_ug_row,n_non_singular_col,cluster_set1,howmanycluster,Ef_ug);
            for (i=1;i<numprocs;i++) {
                mtype = FROM_WORKER;
                MPI_Recv(&err2_sum_from_worker, 1, MPI_DOUBLE, i, mtype, MPI_COMM_WORLD, &status);
                err2_sum_from_master += err2_sum_from_worker;
            }
            cvs = sqrt(err2_sum_from_master/n_corr_mat_ug_row);
            if (cvs < cvs_old)
                cluster_set2[select2] = target;
            else if (exp(-(cvs-cvs_old)/kT) > drand48())
                cluster_set2[select2] = target;
            else {
                cluster_set1[select1] = target;
                cvs = cvs_old;
            }
            if (cvs < cvs_min) {
                for (i=0;i<howmanycluster;i++)
                    cluster_set1_min[i] = cluster_set1[i];
                cvs_min = cvs;
            }
            if (cvs < cvs_tol)
                break;
            if (iter%dispfreq==0)
                printf("iter=%8d, cvs=%f, cvs_min=%f, cvs_min/fu=%f\n",iter,cvs,cvs_min,cvs_min/72);
        }

        // print the result
        printf("End of MC iteration. The selected clusters are...\n");
        printf("ID in non-singular | ID in entire\n");
        for (i=0;i<howmanycluster;i++)
            printf("    %4d                %4d\n",cluster_set1_min[i],usefulcorr_col[cluster_set1_min[i]]);

        double err_pred, rms;
        double Ef_pred[n_corr_mat_ug_row];
        double *corr_mat_ugs_reduced;
        corr_mat_ugs_reduced = (double*) malloc(n_corr_mat_ug_row*howmanycluster*sizeof(double));
        double eci[howmanycluster];
        for (i=0;i<n_corr_mat_ug_row;i++) {
            for (j=0;j<howmanycluster;j++)
                *(corr_mat_ugs_reduced+howmanycluster*i+j) = *(corr_mat_ugs+n_non_singular_col*i+cluster_set1_min[j]);
        }
        least_square_solver(corr_mat_ugs_reduced,n_corr_mat_ug_row,howmanycluster,Ef_ug,eci);
        mat2d_prod(Ef_pred,n_corr_mat_ug_row,1,corr_mat_ugs_reduced,n_corr_mat_ug_row,howmanycluster,eci,howmanycluster,1);
        err_pred = 0.0;
        for (i=0;i<n_corr_mat_ug_row;i++)
            err_pred += (Ef_ug[i] - Ef_pred[i])*(Ef_ug[i] - Ef_pred[i]);
        rms = sqrt(err_pred/n_corr_mat_ug_row);

        // Write result into files
        FILE *fp2, *fp3;
        char cluster_set_filename[200];
        char eci_filename[100];
        strcpy(dirname,"dir_result");
        struct stat st = {0};
        if (stat(dirname, &st) == -1) {
            mkdir(dirname,0777);
            printf("Created a new directory for storing result.\n");
        }
        sprintf(cluster_set_filename, "%s/cluster_set_N%d.dat",dirname,howmanycluster);
        sprintf(eci_filename,"%s/x.dat",dirname);
        fp2 = fopen(cluster_set_filename,"w");
        fp3 = fopen(eci_filename,"w");
        for (i=0;i<howmanycluster;i++) {
            fprintf(fp2,"%4d\t%4d\n",cluster_set1_min[i],usefulcorr_col[cluster_set1_min[i]]);
            fprintf(fp3,"%.4f\n",eci[i]);
        }
    }
    else {
        for (iter=0;iter<(max_iter+1);iter++) {
            mtype = FROM_MASTER;
            MPI_Recv(&cluster_set1, howmanycluster, MPI_INT, MASTER, mtype, MPI_COMM_WORLD, &status);
            row_ini = rank * row_offset;
            row_end = row_ini + row_offset;
            if (rank == numprocs-1)
                row_end = n_corr_mat_ug_row;  // last processor till the last row, -1 is taken inside obtainerr2_par
            err2_sum_from_worker = obtainerr2_par(row_ini,row_end,corr_mat_ugs,n_corr_mat_ug_row,n_non_singular_col,cluster_set1,howmanycluster,Ef_ug);
            mtype = FROM_WORKER;
            MPI_Send(&err2_sum_from_worker, 1, MPI_DOUBLE, MASTER, mtype, MPI_COMM_WORLD);
        }
    }
    MPI_Finalize();

}

void load_double_mat(double *A, int Arow, int Acol, int Apgs, char *datfilename) {

    FILE *fp;
    int i, j, k;
    double tmp;

    printf("filename: %s\n", datfilename);
    fp = fopen(datfilename, "r");
    for (k=0;k<Apgs;k++)
        for (i=0;i<Arow;i++)
            for (j=0;j<Acol;j++) {
                fscanf(fp, "%lf", &tmp);
                *(A+Arow*Acol*k+Acol*i+j) = tmp;
            }

    fclose(fp);

}

void load_int_mat(int *A, int Arow, int Acol, int Apgs, char *datfilename) {

    FILE *fp;
    int i, j, k, tmp;

    printf("filename: %s\n", datfilename);
    fp = fopen(datfilename, "r");
    for (k=0;k<Apgs;k++)
        for (i=0;i<Arow;i++)
            for (j=0;j<Acol;j++) {
                fscanf(fp, "%d", &tmp);
                *(A+Arow*Acol*k+Acol*i+j) = tmp;
            }

    fclose(fp);

}

void load_int_array(int *array, int arraysize, char *datfilename) {

    FILE *fp;
    int i;

    printf("filename: %s\n", datfilename);
    fp = fopen(datfilename, "r");
    //    printf("step2\n");
    for (i=0;i<arraysize;i++) {
        fscanf(fp, "%d", &array[i]);
        //        printf("%d\n",array[i]);
    }

    fclose(fp);
    //    printf("step3\n");
}

void load_double_array(double *array, int arraysize, char *datfilename) {

    FILE *fp;
    int i;

    printf("filename: %s\n", datfilename);
    fp = fopen(datfilename, "r");
    //    printf("step4\n");
    for (i=0;i<arraysize;i++) {
        fscanf(fp, "%lf", &array[i]);
        //        printf("%f\n",array[i]);
    }

    fclose(fp);
    //    printf("step5\n");
}

void randperm(int *a, int np) {

    int i, id1, id2, id3, tmpr, num_perm;

    num_perm = 100*np;

    for (i=0;i<np;i++)
        *(a+i) = i;

    for (i=0;i<num_perm;i++) {
        id1 = rand() % np;
        id3 = rand() % np;
        id2 = (id1 + id3) % np;
        if (id1 != id2) {
            tmpr = *(a+id2);
            *(a+id2) = *(a+id1);
            *(a+id1) = tmpr;
        }
    }
}

int find_value_in_array(int *A, int target_rowA, int colA, int target_value) {
    int i;

    for (i=0;i<colA;i++)
        if (*(A+colA*target_rowA+i)==target_value)
            return(i);
    if (i==colA)
        return(-1);
}


text:
// Copyright 2019 Victor Hugo Schulz

// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at

//     http://www.apache.org/licenses/LICENSE-2.0

// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef ATTITUDE_DETERMINATION_H
#define ATTITUDE_DETERMINATION_H

#include <iostream> // NEEDED FOR DEBUG ONLY
#include <vector>
#include <cmath>
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
extern "C" {
#include <cblas.h>
#include <clapack.h>
}

#include "structures_st.h"
#include "standard_structures.h"

namespace st{

class AttitudeDetermination
{
public:
    AttitudeDetermination();
    std_str::Sky process(std_str::Sky identified);
private:
    double trace(double *m, int size);
    void eye(double *m, int size, double alpha);
};

}

#endif // ATTITUDE_DETERMINATION_H


text:
/*
 *  emu.c
 *  Edited and renamed to fr_emu.c for cosmo_pmc by Martin Kilbinger 2013.
 *
 *
 *  Created by Earl Lawrence on 9/17/09.
 *  Update 10/9/2012
 *
 *  This program was prepared by Los Alamos National Security, LLC at Los Alamos National Laboratory (LANL)
 *  under contract No. DE-AC52-06NA25396 with the U.S. Department of Energy (DOE). All rights in the program
 *  are reserved by the DOE and Los Alamos National Security, LLC.  Permission is granted to the public to
 *  copy and use this software without charge, provided that this Notice and any statement of authorship are
 *  reproduced on all copies.  Neither the U.S. Government nor LANS makes any warranty, express or implied,
 *  or assumes any liability or responsibility for the use of this software.
 *
 *
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include "coyote.h"
#include "fr_constants.h"

/* Copy parameters to emulator input array */
void fill_xstar6_wo_z(double omega_m, double omega_b, double n_spec, double sigma_8, double w0_de, double h_100,
                      double xstar[])
{
   /* Orders reversed:  (omega_b <-> omega_m), (w0_de <-> sigma_8), compared to v1 */
   xstar[0] = omega_b;
   xstar[1] = omega_m;
   xstar[2] = n_spec;
   xstar[3] = h_100 * 100;
   xstar[4] = w0_de;
   xstar[5] = sigma_8;
}

/* ============================================================ *
 * Returns the Coyote non-linear power spectrum. k is in units  *
 * of h/Mpc, h_100 is now an independent parameter (Coyote v2). *
 * It is also used to transform it to [1/Mpc] as                *
 * needed by the Coyote emulator.                               *
 * Six parameters, for coyote13.				*
 * ============================================================ */
double P_NL_coyote6(double omega_m, double omega_b, double n_spec, double sigma_8, double w0_de,
		    double h_100, double a, double k, double **ystar_allz, error **err)
{
   double xstar[p], ystar[2*fr_nsim], val, k_invMpc;
   int ihi, ilo, i;

   /* Copy cosmological parameters */
   fill_xstar6_wo_z(omega_m, omega_b, n_spec, sigma_8, w0_de, h_100, xstar);

   /* Redshift */
   xstar[p] = 1.0/a - 1.0;

   /* Input k is [h/Mpc], Coyote interpolates on k_invMpc in [1/Mpc] */
   k_invMpc  = k * h_100;

   /* Check k-range */
   testErrorRetVA(k_invMpc < fr_ksim[0] || k_invMpc > fr_ksim[fr_nsim-1], coyote_range,
           "Fourier mode k=%g 1/Mpc out of range [%g;%g]",
           *err, __LINE__, -1.0, k_invMpc, fr_ksim[0], fr_ksim[fr_nsim-1]);

   if (*ystar_allz == NULL) {
      *ystar_allz = malloc_err(sizeof(double) * fr_rs * fr_nsim, err);
      forwardError(*err, __LINE__, -1.0);

      /* The actual emulator, filling the array in (k, z) */
      //testErrorRet(ystar_allz == NULL, io_null, "Coyote13: array for (k,z) 'ystar_allz' not initialised",
      //	   *err, __LINE__,);
      fr_fill_ystar_allz(*ystar_allz, xstar, err);
      forwardError(*err, __LINE__, -1.0);
   }

   /* ystar: first half = k, second half = P(k) */
   fr_emu(xstar, ystar, *ystar_allz, err);
   forwardError(*err, __LINE__, -1.0);

   /* Looking for right k-index with bisection */
   ilo = 0;
   ihi = fr_nsim-1;
   while (ihi - ilo > 1) {
      i = (ihi + ilo) >> 1;
      if (ystar[i] > k_invMpc) ihi = i;
      else ilo = i;
   }
   testErrorRetVA(ihi == ilo, math_wrongValue, "Bisection failed, both indices equal (%d)", *err, __LINE__, -1.0, ihi);

   /* Linear interpolation */
   val = (k_invMpc - ystar[ilo]) / (ystar[ihi] - ystar[ilo]) * (ystar[fr_nsim + ihi]
         - ystar[fr_nsim + ilo]) + ystar[fr_nsim + ilo];

   /* Coyote P(k) [Mpc^3] -> output P(k) [(Mpc/h)^3] */
   val *= h_100 * h_100 * h_100;

   return val;
}

void fr_check_range(const double *xstar, error **err)
{
   int i;

   for(i=0; i<p; i++) {
      testErrorRetVA((xstar[i] < fr_xmin[i]) || (xstar[i] > fr_xmin[i]+fr_xrange[i]), coyote_range,
		     "Coyote FrankenEmu parameter #%d (%g) out of range [%g; %g]",
		     *err, __LINE__,, i, xstar[i], fr_xmin[i], fr_xmin[i] + fr_xrange[i]);
   } // for(i=0; i<p; i++)

   // Check redshift to make sure we're interpolating
   testErrorRetVA((xstar[p] < 0) || (xstar[p] > 4), coyote_range,
		  "Redshift %g out of range [%g;%g]", *err, __LINE__,, xstar[p], 0.0, 4.0);
}

/* ============================================================ *
 * The actual emulation. Fills the array ystar_allz for (k, z). *
 * Cosmological parameters, placeholder for the output, type of *
 * output.
 * xstar[p] = redshift is not used here, only the p=6 cosmo.    *
 * parameters.							*
 * ============================================================ */
void fr_fill_ystar_allz(double *ystar_allz, const double *xstar, error **err)
{
   const double sd = 0.16002;
   int i, j, k;
   double xstarstd[p], wstar[peta], Sigmastar[peta][m], ystaremu[neta];
   double logc;

   // Interpolation stuff for k
   gsl_spline *lininterp_k = gsl_spline_alloc(gsl_interp_linear, neta/rs);
   gsl_interp_accel *accel = gsl_interp_accel_alloc();

   //fprintf(stderr, "MKDEBUG: Filling ystar_allz for (%g, %g, %g, %g, %g, %g)\n",
   //	   xstar[0], xstar[1], xstar[2], xstar[3], xstar[4], xstar[5]);

    // Standardize the inputs
    for(i=0; i<p; i++) {
        xstarstd[i] = (xstar[i] - fr_xmin[i]) / fr_xrange[i];
    }

    // Compute the covariances between the new input and sims for all PCs
    for(i=0; i<peta; i++) {
        for(j=0; j<m; j++) {
            logc = 0.0;
            for(k=0; k<p; k++) {
                logc -= fr_beta[i][k]*dsqr(fr_x[j][k]-xstarstd[k]);
            }
            Sigmastar[i][j] = exp(logc)/fr_lamz[i];
        }
    }

    // Compute wstar, the predicted PC weights for the new input
    for(i=0; i<peta; i++) {
        wstar[i]=0.0;
        for(j=0; j<m; j++) {
            wstar[i] += Sigmastar[i][j] * fr_KrigBasis[i][j];
        }
    }

    // Compute ystar, the new output
    for(i=0; i<neta; i++) {
        ystaremu[i] = 0.0;
        for(j=0; j<peta; j++) {
            ystaremu[i] += fr_K[i][j]*wstar[j];
        }
        ystaremu[i] = ystaremu[i]*sd + fr_mean[i];
    }

    // Interpolate the emulated output onto the original domain.
    for(i=0; i<rs; i++) {
        gsl_spline_init(lininterp_k, fr_kemu, &ystaremu[i*neta/rs], neta/rs);
        for(j=0; j<fr_nsim; j++) {
            ystar_allz[i*fr_nsim+j] = gsl_spline_eval(lininterp_k, fr_ksim[j], accel);
        }
        gsl_interp_accel_reset(accel);
    }

    gsl_spline_free(lininterp_k);
    gsl_interp_accel_free(accel);
}

/* ============================================================ *
 * This used to be the actual emulation, now fr_fill_ystar_allz *
 * is called from here.						*
 * xstar[p+1] contains the six cosmological parameters plus the *
 * redshift.							*
 * ============================================================ */
void fr_emu(const double *xstar, double *ystar, const double *ystar_allz, error **err)
{
    int i, j;
    //double ystar_allz[rs*fr_nsim];
    double zemu[rs], ybyz[rs];

    /* CosmoPMC version of Coyote: no initialisation, no static variables */

    // Interpolation stuff z
    gsl_spline *lininterp_z = gsl_spline_alloc(gsl_interp_linear, rs);
    gsl_interp_accel *accel = gsl_interp_accel_alloc();

    // Check the inputs to make sure we're interpolating.
    fr_check_range(xstar, err);
    forwardError(*err, __LINE__,);

    // Fill in the k values for the final output; first half of ystar
    for(i=0; i<fr_nsim; i++) {
        ystar[i] = fr_ksim[i];
    }

    // Interpolate on to the desired redshift
    // The order needs to be reversed here.
    for(i=0; i<rs; i++) {
        zemu[i] = (1/fr_aemu[rs-i-1]) - 1.0;
    }
    for(i=0; i<fr_nsim; i++) {
        // Build an array with the values of y for a given value of z
        // Reverse the order
        for(j=0; j<rs; j++) {
            ybyz[rs-j-1] = ystar_allz[j*fr_nsim+i];
        }
        gsl_spline_init(lininterp_z, zemu, ybyz, rs);
        ystar[fr_nsim+i] = gsl_spline_eval(lininterp_z, xstar[p], accel);
        gsl_interp_accel_reset(accel);
    }

    //printf("z interped\n");

    // Transform to P(k); second half of ystar
    for(j=0; j<fr_nsim; j++) {
       ystar[fr_nsim+j] = ystar[fr_nsim+j] - 1.5*log10(fr_ksim[j]);
       ystar[fr_nsim+j] = pow(10.0,ystar[fr_nsim+j])*2.0*M_PI*M_PI;
    }

    // Free some stuff.  I always forget this.
    // Thanks to Tim Eifler for discovering it (apparently the hard way).
    gsl_spline_free(lininterp_z);
    gsl_interp_accel_free(accel);
}

/* Undefine constants from fr_constants.h */
#undef m
#undef neta
#undef p
#undef peta
#undef rs

max_stars_repo_name: danielgruen/ccv

max_issues_repo_name: danielgruen/ccv

max_forks_repo_name: danielgruen/ccv


text:
/* multifit/multireg.c
 *
 * Copyright (C) 2015 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.
 */

/*
 * References:
 *
 * [1] P. C. Hansen & D. P. O'Leary, "The use of the L-curve in
 * the regularization of discrete ill-posed problems",  SIAM J. Sci.
 * Comput. 14 (1993), pp. 1487-1503.
 *
 * [2] P. C. Hansen, "Discrete Inverse Problems: Insight and Algorithms,"
 * SIAM Press, 2010.
 */

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

#include "linear_common.c"

int
gsl_multifit_linear_solve (const double lambda,
                           const gsl_matrix * X,
                           const gsl_vector * y,
                           gsl_vector * c,
                           double *rnorm,
                           double *snorm,
                           gsl_multifit_linear_workspace * work)
{
  size_t rank;
  int status;

  status = multifit_linear_solve(X, y, GSL_DBL_EPSILON, lambda, &rank, c,
                                 rnorm, snorm, work);

  return status;
} /* gsl_multifit_linear_solve() */

/*
gsl_multifit_linear_applyW()
  Apply weight matrix to (X,y) LS system

Inputs: X    - least squares matrix n-by-p
        w    - weight vector n-by-1 or NULL for W = I
        y    - right hand side n-by-1
        WX   - (output) sqrt(W) X, n-by-p
        Wy   - (output) sqrt(W) y, n-by-1

Notes:
1) If w = NULL, on output WX = X and Wy = y
2) It is allowed for WX = X and Wy = y for in-place transform
*/

int
gsl_multifit_linear_applyW(const gsl_matrix * X,
                           const gsl_vector * w,
                           const gsl_vector * y,
                           gsl_matrix * WX,
                           gsl_vector * Wy)
{
  const size_t n = X->size1;
  const size_t p = X->size2;

  if (n != y->size)
    {
      GSL_ERROR("y vector does not match X", GSL_EBADLEN);
    }
  else if (w != NULL && n != w->size)
    {
      GSL_ERROR("weight vector does not match X", GSL_EBADLEN);
    }
  else if (n != WX->size1 || p != WX->size2)
    {
      GSL_ERROR("WX matrix dimensions do not match X", GSL_EBADLEN);
    }
  else if (n != Wy->size)
    {
      GSL_ERROR("Wy vector must be length n", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* copy WX = X; Wy = y if distinct pointers */
      if (WX != X)
        gsl_matrix_memcpy(WX, X);
      if (Wy != y)
        gsl_vector_memcpy(Wy, y);

      if (w != NULL)
        {
          /* construct WX = sqrt(W) X and Wy = sqrt(W) y */
          for (i = 0; i < n; ++i)
            {
              double wi = gsl_vector_get(w, i);
              double swi;
              gsl_vector_view row = gsl_matrix_row(WX, i);
              double *yi = gsl_vector_ptr(Wy, i);

              if (wi < 0.0)
                wi = 0.0;

              swi = sqrt(wi);
              gsl_vector_scale(&row.vector, swi);
              *yi *= swi;
            }
        }

      return GSL_SUCCESS;
    }
}

/*
gsl_multifit_linear_wstdform1()
  Using regularization matrix
L = diag(l_1,l_2,...,l_p), transform to Tikhonov standard form:

X~ = sqrt(W) X L^{-1}
y~ = sqrt(W) y
c~ = L c

Inputs: L    - Tikhonov matrix as a vector of diagonal elements p-by-1;
               or NULL for L = I
        X    - least squares matrix n-by-p
        y    - right hand side vector n-by-1
        w    - weight vector n-by-1; or NULL for W = I
        Xs   - least squares matrix in standard form X~ n-by-p
        ys   - right hand side vector in standard form y~ n-by-1
        work - workspace

Return: success/error

Notes:
1) It is allowed for X = Xs and y = ys
*/

int
gsl_multifit_linear_wstdform1 (const gsl_vector * L,
                               const gsl_matrix * X,
                               const gsl_vector * w,
                               const gsl_vector * y,
                               gsl_matrix * Xs,
                               gsl_vector * ys,
                               gsl_multifit_linear_workspace * work)
{
  const size_t n = X->size1;
  const size_t p = X->size2;

  if (n > work->nmax || p > work->pmax)
    {
      GSL_ERROR("observation matrix larger than workspace", GSL_EBADLEN);
    }
  else if (L != NULL && p != L->size)
    {
      GSL_ERROR("L vector does not match X", GSL_EBADLEN);
    }
  else if (n != y->size)
    {
      GSL_ERROR("y vector does not match X", GSL_EBADLEN);
    }
  else if (w != NULL && n != w->size)
    {
      GSL_ERROR("weight vector does not match X", GSL_EBADLEN);
    }
  else if (n != Xs->size1 || p != Xs->size2)
    {
      GSL_ERROR("Xs matrix dimensions do not match X", GSL_EBADLEN);
    }
  else if (n != ys->size)
    {
      GSL_ERROR("ys vector must be length n", GSL_EBADLEN);
    }
  else
    {
      int status = GSL_SUCCESS;

      /* compute Xs = sqrt(W) X and ys = sqrt(W) y */
      status = gsl_multifit_linear_applyW(X, w, y, Xs, ys);
      if (status)
        return status;

      if (L != NULL)
        {
          size_t j;

          /* construct X~ = sqrt(W) X * L^{-1} matrix */
          for (j = 0; j < p; ++j)
            {
              gsl_vector_view Xj = gsl_matrix_column(Xs, j);
              double lj = gsl_vector_get(L, j);

              if (lj == 0.0)
                {
                  GSL_ERROR("L matrix is singular", GSL_EDOM);
                }

              gsl_vector_scale(&Xj.vector, 1.0 / lj);
            }
        }

      return status;
    }
}

/*
gsl_multifit_linear_stdform1()
  Using regularization matrix L = diag(l_1,l_2,...,l_p),
and W = I, transform to Tikhonov standard form:

X~ = X L^{-1}
y~ = y
c~ = L c

Inputs: L    - Tikhonov matrix as a vector of diagonal elements p-by-1
        X    - least squares matrix n-by-p
        y    - right hand side vector n-by-1
        Xs   - least squares matrix in standard form X~ n-by-p
        ys   - right hand side vector in standard form y~ n-by-1
        work - workspace

Return: success/error

Notes:
1) It is allowed for X = Xs
*/

int
gsl_multifit_linear_stdform1 (const gsl_vector * L,
                              const gsl_matrix * X,
                              const gsl_vector * y,
                              gsl_matrix * Xs,
                              gsl_vector * ys,
                              gsl_multifit_linear_workspace * work)
{
  int status;

  status = gsl_multifit_linear_wstdform1(L, X, NULL, y, Xs, ys, work);

  return status;
}

int
gsl_multifit_linear_L_decomp (gsl_matrix * L, gsl_vector * tau)
{
  const size_t m = L->size1;
  const size_t p = L->size2;
  int status;

  if (tau->size != GSL_MIN(m, p))
    {
      GSL_ERROR("tau vector must be min(m,p)", GSL_EBADLEN);
    }
  else if (m >= p)
    {
      /* square or tall L matrix */
      status = gsl_linalg_QR_decomp(L, tau);
      return status;
    }
  else
    {
      /* more columns than rows, compute qr(L^T) */
      gsl_matrix_view LTQR = gsl_matrix_view_array(L->data, p, m);
      gsl_matrix *LT = gsl_matrix_alloc(p, m);

      /* XXX: use temporary storage due to difficulties in transforming
       * a rectangular matrix in-place */
      gsl_matrix_transpose_memcpy(LT, L);
      gsl_matrix_memcpy(&LTQR.matrix, LT);
      gsl_matrix_free(LT);

      status = gsl_linalg_QR_decomp(&LTQR.matrix, tau);

      return status;
    }
}

/*
gsl_multifit_linear_wstdform2()
  Using regularization matrix L which is m-by-p, transform to Tikhonov
standard form. This routine is separated into two cases:

Case 1: m >= p, here we can use the QR decomposition of L = QR, and note
that ||L c|| = ||R c|| where R is p-by-p. Therefore,

X~ = X R^{-1} is n-by-p
y~ = y is n-by-1
c~ is p-by-1
M is not used

Case 2: m < p

X~ is (n - p + m)-by-m
y~ is (n - p + m)-by-1
c~ is m-by-1
M is n-by-p (workspace)

Inputs: LQR  - output from gsl_multifit_linear_L_decomp()
        Ltau - output from gsl_multifit_linear_L_decomp()
        X    - least squares matrix n-by-p
        w    - weight vector n-by-1; or NULL for W = I
        y    - right hand side vector n-by-1
        Xs   - (output) least squares matrix in standard form
               case 1: n-by-p
               case 2: (n - p + m)-by-m
        ys   - (output) right hand side vector in standard form
               case 1: n-by-1
               case 2: (n - p + m)-by-1
        M    - (output) workspace matrix needed to reconstruct solution vector
               case 1: not used
               case 2: n-by-p
        work - workspace

Return: success/error

Notes:
1) If m >= p, on output:
   Xs = X R^{-1}
   ys = y

2) If m < p, on output:
   M(:,1:pm) contains QR decomposition of A * K_o, needed to reconstruct
   solution vector, where pm = p - m; M(:,p) contains Householder scalars
*/

int
gsl_multifit_linear_wstdform2 (const gsl_matrix * LQR,
                               const gsl_vector * Ltau,
                               const gsl_matrix * X,
                               const gsl_vector * w,
                               const gsl_vector * y,
                               gsl_matrix * Xs,
                               gsl_vector * ys,
                               gsl_matrix * M,
                               gsl_multifit_linear_workspace * work)
{
  const size_t m = LQR->size1;
  const size_t n = X->size1;
  const size_t p = X->size2;

  if (n > work->nmax || p > work->pmax)
    {
      GSL_ERROR("observation matrix larger than workspace", GSL_EBADLEN);
    }
  else if (p != LQR->size2)
    {
      GSL_ERROR("LQR and X matrices have different numbers of columns", GSL_EBADLEN);
    }
  else if (n != y->size)
    {
      GSL_ERROR("y vector does not match X", GSL_EBADLEN);
    }
  else if (w != NULL && n != w->size)
    {
      GSL_ERROR("weights vector must be length n", GSL_EBADLEN);
    }
  else if (m >= p) /* square or tall L matrix */
    {
      /* the sizes of Xs and ys depend on whether m >= p or m < p */
      if (n != Xs->size1 || p != Xs->size2)
        {
          GSL_ERROR("Xs matrix must be n-by-p", GSL_EBADLEN);
        }
      else if (n != ys->size)
        {
          GSL_ERROR("ys vector must have length n", GSL_EBADLEN);
        }
      else
        {
          int status;
          size_t i;
          gsl_matrix_const_view R = gsl_matrix_const_submatrix(LQR, 0, 0, p, p);

          /* compute Xs = sqrt(W) X and ys = sqrt(W) y */
          status = gsl_multifit_linear_applyW(X, w, y, Xs, ys);
          if (status)
            return status;

          /* compute X~ = X R^{-1} using QR decomposition of L */
          for (i = 0; i < n; ++i)
            {
              gsl_vector_view v = gsl_matrix_row(Xs, i);

              /* solve: R^T y = X_i */
              gsl_blas_dtrsv(CblasUpper, CblasTrans, CblasNonUnit, &R.matrix, &v.vector);
            }

          return GSL_SUCCESS;
        }
    }
  else /* L matrix with m < p */
    {
      const size_t pm = p - m;
      const size_t npm = n - pm;

      /*
       * This code closely follows section 2.6.1 of Hansen's
       * "Regularization Tools" manual
       */

      if (npm != Xs->size1 || m != Xs->size2)
        {
          GSL_ERROR("Xs matrix must be (n-p+m)-by-m", GSL_EBADLEN);
        }
      else if (npm != ys->size)
        {
          GSL_ERROR("ys vector must be of length (n-p+m)", GSL_EBADLEN);
        }
      else if (n != M->size1 || p != M->size2)
        {
          GSL_ERROR("M matrix must be n-by-p", GSL_EBADLEN);
        }
      else
        {
          int status;
          gsl_matrix_view A = gsl_matrix_submatrix(work->A, 0, 0, n, p);
          gsl_vector_view b = gsl_vector_subvector(work->t, 0, n);

          gsl_matrix_view LTQR = gsl_matrix_view_array(LQR->data, p, m);           /* qr(L^T) */
          gsl_matrix_view Rp = gsl_matrix_view_array(LQR->data, m, m);             /* R factor of L^T */
          gsl_vector_const_view LTtau = gsl_vector_const_subvector(Ltau, 0, m);

          /*
           * M(:,1:p-m) will hold QR decomposition of A K_o; M(:,p) will hold
           * Householder scalars
           */
          gsl_matrix_view MQR = gsl_matrix_submatrix(M, 0, 0, n, pm);
          gsl_vector_view Mtau = gsl_matrix_subcolumn(M, p - 1, 0, GSL_MIN(n, pm));

          gsl_matrix_view AKo, AKp, HqTAKp;
          gsl_vector_view v;
          size_t i;

          /* compute A = sqrt(W) X and b = sqrt(W) y */
          status = gsl_multifit_linear_applyW(X, w, y, &A.matrix, &b.vector);
          if (status)
            return status;

          /* compute: A <- A K = [ A K_p ; A K_o ] */
          gsl_linalg_QR_matQ(&LTQR.matrix, &LTtau.vector, &A.matrix);
          AKp = gsl_matrix_submatrix(&A.matrix, 0, 0, n, m);
          AKo = gsl_matrix_submatrix(&A.matrix, 0, m, n, pm);

          /* compute QR decomposition [H,T] = qr(A * K_o) and store in M */
          gsl_matrix_memcpy(&MQR.matrix, &AKo.matrix);
          gsl_linalg_QR_decomp(&MQR.matrix, &Mtau.vector);

          /* AKp currently contains A K_p; apply H^T from the left to get H^T A K_p */
          gsl_linalg_QR_QTmat(&MQR.matrix, &Mtau.vector, &AKp.matrix);

          /* the last npm rows correspond to H_q^T A K_p */
          HqTAKp = gsl_matrix_submatrix(&AKp.matrix, pm, 0, npm, m);

          /* solve: Xs R_p^T = H_q^T A K_p for Xs */
          gsl_matrix_memcpy(Xs, &HqTAKp.matrix);
          for (i = 0; i < npm; ++i)
            {
              gsl_vector_view x = gsl_matrix_row(Xs, i);
              gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, &Rp.matrix, &x.vector);
            }

          /*
           * compute: ys = H_q^T b; this is equivalent to computing
           * the last q elements of H^T b (q = npm)
           */
          v = gsl_vector_subvector(&b.vector, pm, npm);
          gsl_linalg_QR_QTvec(&MQR.matrix, &Mtau.vector, &b.vector);
          gsl_vector_memcpy(ys, &v.vector);

          return GSL_SUCCESS;
        }
    }
}

int
gsl_multifit_linear_stdform2 (const gsl_matrix * LQR,
                              const gsl_vector * Ltau,
                              const gsl_matrix * X,
                              const gsl_vector * y,
                              gsl_matrix * Xs,
                              gsl_vector * ys,
                              gsl_matrix * M,
                              gsl_multifit_linear_workspace * work)
{
  int status;

  status = gsl_multifit_linear_wstdform2(LQR, Ltau, X, NULL, y, Xs, ys, M, work);

  return status;
}

/*
gsl_multifit_linear_genform1()
  Backtransform regularized solution vector using matrix
L = diag(L)
*/

int
gsl_multifit_linear_genform1 (const gsl_vector * L,
                              const gsl_vector * cs,
                              gsl_vector * c,
                              gsl_multifit_linear_workspace * work)
{
  if (L->size > work->pmax)
    {
      GSL_ERROR("L vector does not match workspace", GSL_EBADLEN);
    }
  else if (L->size != cs->size)
    {
      GSL_ERROR("cs vector does not match L", GSL_EBADLEN);
    }
  else if (L->size != c->size)
    {
      GSL_ERROR("c vector does not match L", GSL_EBADLEN);
    }
  else
    {
      /* compute true solution vector c = L^{-1} c~ */
      gsl_vector_memcpy(c, cs);
      gsl_vector_div(c, L);

      return GSL_SUCCESS;
    }
}

/*
gsl_multifit_linear_wgenform2()
  Backtransform regularized solution vector in standard form to recover
original vector

Inputs: LQR  - output from gsl_multifit_linear_L_decomp()
        Ltau - output from gsl_multifit_linear_L_decomp()
        X    - original least squares matrix n-by-p
        w    - original weight vector n-by-1 or NULL for W = I
        y    - original rhs vector n-by-1
        cs   - standard form solution vector
        c    - (output) original solution vector p-by-1
        M    - matrix computed by gsl_multifit_linear_wstdform2()
        work - workspace
*/

int
gsl_multifit_linear_wgenform2 (const gsl_matrix * LQR,
                               const gsl_vector * Ltau,
                               const gsl_matrix * X,
                               const gsl_vector * w,
                               const gsl_vector * y,
                               const gsl_vector * cs,
                               const gsl_matrix * M,
                               gsl_vector * c,
                               gsl_multifit_linear_workspace * work)
{
  const size_t m = LQR->size1;
  const size_t n = X->size1;
  const size_t p = X->size2;

  if (n > work->nmax || p > work->pmax)
    {
      GSL_ERROR("X matrix does not match workspace", GSL_EBADLEN);
    }
  else if (p != LQR->size2)
    {
      GSL_ERROR("LQR matrix does not match X", GSL_EBADLEN);
    }
  else if (p != c->size)
    {
      GSL_ERROR("c vector does not match X", GSL_EBADLEN);
    }
  else if (w != NULL && n != w->size)
    {
      GSL_ERROR("w vector does not match X", GSL_EBADLEN);
    }
  else if (n != y->size)
    {
      GSL_ERROR("y vector does not match X", GSL_EBADLEN);
    }
  else if (m >= p)                    /* square or tall L matrix */
    {
      if (p != cs->size)
        {
          GSL_ERROR("cs vector must be length p", GSL_EBADLEN);
        }
      else
        {
          int s;
          gsl_matrix_const_view R = gsl_matrix_const_submatrix(LQR, 0, 0, p, p); /* R factor of L */

          /* solve R c = cs for true solution c, using QR decomposition of L */
          gsl_vector_memcpy(c, cs);
          s = gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, &R.matrix, c);

          return s;
        }
    }
  else                                /* rectangular L matrix with m < p */
    {
      if (m != cs->size)
        {
          GSL_ERROR("cs vector must be length m", GSL_EBADLEN);
        }
      else if (n != M->size1 || p != M->size2)
        {
          GSL_ERROR("M matrix must be size n-by-p", GSL_EBADLEN);
        }
      else
        {
          int status;
          const size_t pm = p - m;
          gsl_matrix_view A = gsl_matrix_submatrix(work->A, 0, 0, n, p);
          gsl_vector_view b = gsl_vector_subvector(work->t, 0, n);
          gsl_matrix_view Rp = gsl_matrix_view_array(LQR->data, m, m); /* R_p */
          gsl_matrix_view LTQR = gsl_matrix_view_array(LQR->data, p, m);
          gsl_vector_const_view LTtau = gsl_vector_const_subvector(Ltau, 0, m);
          gsl_matrix_const_view MQR = gsl_matrix_const_submatrix(M, 0, 0, n, pm);
          gsl_vector_const_view Mtau = gsl_matrix_const_subcolumn(M, p - 1, 0, GSL_MIN(n, pm));
          gsl_matrix_const_view To = gsl_matrix_const_submatrix(&MQR.matrix, 0, 0, pm, pm);
          gsl_vector_view workp = gsl_vector_subvector(work->xt, 0, p);
          gsl_vector_view v1, v2;

          /* compute A = sqrt(W) X and b = sqrt(W) y */
          status = gsl_multifit_linear_applyW(X, w, y, &A.matrix, &b.vector);
          if (status)
            return status;

          /* initialize c to zero */
          gsl_vector_set_zero(c);

          /* compute c = L_inv cs = K_p R_p^{-T} cs */

          /* set c(1:m) = R_p^{-T} cs */
          v1 = gsl_vector_subvector(c, 0, m);
          gsl_vector_memcpy(&v1.vector, cs);
          gsl_blas_dtrsv(CblasUpper, CblasTrans, CblasNonUnit, &Rp.matrix, &v1.vector);

          /* c <- K R_p^{-T} cs = [ K_p R_p^{_T} cs ; 0 ] */
          gsl_linalg_QR_Qvec(&LTQR.matrix, &LTtau.vector, c);

          /* compute: b1 = b - A L_inv cs */
          gsl_blas_dgemv(CblasNoTrans, -1.0, &A.matrix, c, 1.0, &b.vector);

          /* compute: b2 = H^T b1 */
          gsl_linalg_QR_QTvec(&MQR.matrix, &Mtau.vector, &b.vector);

          /* compute: b3 = T_o^{-1} b2 */
          v1 = gsl_vector_subvector(&b.vector, 0, pm);
          gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, &To.matrix, &v1.vector);

          /* compute: b4 = K_o b3 */
          gsl_vector_set_zero(&workp.vector);
          v2 = gsl_vector_subvector(&workp.vector, m, pm);
          gsl_vector_memcpy(&v2.vector, &v1.vector);
          gsl_linalg_QR_Qvec(&LTQR.matrix, &LTtau.vector, &workp.vector);

          /* final solution vector */
          gsl_vector_add(c, &workp.vector);

          return GSL_SUCCESS;
        }
    }
}

int
gsl_multifit_linear_genform2 (const gsl_matrix * LQR,
                              const gsl_vector * Ltau,
                              const gsl_matrix * X,
                              const gsl_vector * y,
                              const gsl_vector * cs,
                              const gsl_matrix * M,
                              gsl_vector * c,
                              gsl_multifit_linear_workspace * work)
{
  int status;

  status = gsl_multifit_linear_wgenform2(LQR, Ltau, X, NULL, y, cs, M, c, work);

  return status;
}

/*
gsl_multifit_linear_lreg()
  Calculate regularization parameters to use in L-curve
analysis

Inputs: smin      - smallest singular value of LS system
        smax      - largest singular value of LS system > 0
        reg_param - (output) vector of regularization parameters
                    derived from singular values

Return: success/error
*/

int
gsl_multifit_linear_lreg (const double smin, const double smax,
                          gsl_vector * reg_param)
{
  if (smax <= 0.0)
    {
      GSL_ERROR("smax must be positive", GSL_EINVAL);
    }
  else
    {
      const size_t N = reg_param->size;

      /* smallest regularization parameter */
      const double smin_ratio = 16.0 * GSL_DBL_EPSILON;
      const double new_smin = GSL_MAX(smin, smax*smin_ratio);
      double ratio;
      size_t i;

      gsl_vector_set(reg_param, N - 1, new_smin);

      /* ratio so that reg_param(1) = s(1) */
      ratio = pow(smax / new_smin, 1.0 / ((double)N - 1.0));

      /* calculate the regularization parameters */
      for (i = N - 1; i > 0 && i--; )
        {
          double rp1 = gsl_vector_get(reg_param, i + 1);
          gsl_vector_set(reg_param, i, ratio * rp1);
        }

      return GSL_SUCCESS;
    }
}

/*
gsl_multifit_linear_lcurve()
  Calculate L-curve using regularization parameters estimated
from singular values of least squares matrix

Inputs: y         - right hand side vector
        reg_param - (output) vector of regularization parameters
                    derived from singular values
        rho       - (output) vector of residual norms ||y - X c||
        eta       - (output) vector of solution norms ||lambda c||
        work      - workspace

Return: success/error

Notes:
1) SVD of X must be computed first by calling multifit_linear_svd();
   work->n and work->p are initialized by this function
*/

int
gsl_multifit_linear_lcurve (const gsl_vector * y,
                            gsl_vector * reg_param,
                            gsl_vector * rho, gsl_vector * eta,
                            gsl_multifit_linear_workspace * work)
{
  const size_t n = y->size;
  const size_t N = rho->size; /* number of points on L-curve */

  if (n != work->n)
    {
      GSL_ERROR("y vector does not match workspace", GSL_EBADLEN);
    }
  else if (N < 3)
    {
      GSL_ERROR ("at least 3 points are needed for L-curve analysis",
                 GSL_EBADLEN);
    }
  else if (N != eta->size)
    {
      GSL_ERROR ("size of rho and eta vectors do not match",
                 GSL_EBADLEN);
    }
  else if (reg_param->size != eta->size)
    {
      GSL_ERROR ("size of reg_param and eta vectors do not match",
                 GSL_EBADLEN);
    }
  else
    {
      int status = GSL_SUCCESS;
      const size_t p = work->p;

      size_t i, j;

      gsl_matrix_view A = gsl_matrix_submatrix(work->A, 0, 0, n, p);
      gsl_vector_view S = gsl_vector_subvector(work->S, 0, p);
      gsl_vector_view xt = gsl_vector_subvector(work->xt, 0, p);
      gsl_vector_view workp = gsl_matrix_subcolumn(work->QSI, 0, 0, p);
      gsl_vector_view workp2 = gsl_vector_subvector(work->D, 0, p); /* D isn't used for regularized problems */

      const double smax = gsl_vector_get(&S.vector, 0);
      const double smin = gsl_vector_get(&S.vector, p - 1);

      double dr; /* residual error from projection */
      double normy = gsl_blas_dnrm2(y);
      double normUTy;

      /* compute projection xt = U^T y */
      gsl_blas_dgemv (CblasTrans, 1.0, &A.matrix, y, 0.0, &xt.vector);

      normUTy = gsl_blas_dnrm2(&xt.vector);
      dr = normy*normy - normUTy*normUTy;

      /* calculate regularization parameters */
      gsl_multifit_linear_lreg(smin, smax, reg_param);

      for (i = 0; i < N; ++i)
        {
          double lambda = gsl_vector_get(reg_param, i);
          double lambda_sq = lambda * lambda;

          for (j = 0; j < p; ++j)
            {
              double sj = gsl_vector_get(&S.vector, j);
              double xtj = gsl_vector_get(&xt.vector, j);
              double f = sj / (sj*sj + lambda_sq);

              gsl_vector_set(&workp.vector, j, f * xtj);
              gsl_vector_set(&workp2.vector, j, (1.0 - sj*f) * xtj);
            }

          gsl_vector_set(eta, i, gsl_blas_dnrm2(&workp.vector));
          gsl_vector_set(rho, i, gsl_blas_dnrm2(&workp2.vector));
        }

      if (n > p && dr > 0.0)
        {
          /* add correction to residual norm (see eqs 6-7 of [1]) */
          for (i = 0; i < N; ++i)
            {
              double rhoi = gsl_vector_get(rho, i);
              double *ptr = gsl_vector_ptr(rho, i);

              *ptr = sqrt(rhoi*rhoi + dr);
            }
        }

      /* restore D to identity matrix */
      gsl_vector_set_all(work->D, 1.0);

      return status;
    }
} /* gsl_multifit_linear_lcurve() */

/*
gsl_multifit_linear_lcorner()
  Determine point on L-curve of maximum curvature. For each
set of 3 points on the L-curve, the circle which passes through
the 3 points is computed. The radius of the circle is then used
as an estimate of the curvature at the middle point. The point
with maximum curvature is then selected.

Inputs: rho - vector of residual norms ||A x - b||
        eta - vector of solution norms ||L x||
        idx - (output) index i such that
              (log(rho(i)),log(eta(i))) is the point of
              maximum curvature

Return: success/error
*/

int
gsl_multifit_linear_lcorner(const gsl_vector *rho,
                            const gsl_vector *eta,
                            size_t *idx)
{
  const size_t n = rho->size;

  if (n < 3)
    {
      GSL_ERROR ("at least 3 points are needed for L-curve analysis",
                 GSL_EBADLEN);
    }
  else if (n != eta->size)
    {
      GSL_ERROR ("size of rho and eta vectors do not match",
                 GSL_EBADLEN);
    }
  else
    {
      int s = GSL_SUCCESS;
      size_t i;
      double x1, y1;      /* first point of triangle on L-curve */
      double x2, y2;      /* second point of triangle on L-curve */
      double rmin = -1.0; /* minimum radius of curvature */

      /* initial values */
      x1 = log(gsl_vector_get(rho, 0));
      y1 = log(gsl_vector_get(eta, 0));

      x2 = log(gsl_vector_get(rho, 1));
      y2 = log(gsl_vector_get(eta, 1));

      for (i = 1; i < n - 1; ++i)
        {
          /*
           * The points (x1,y1), (x2,y2), (x3,y3) are the previous,
           * current, and next point on the L-curve. We will find
           * the circle which fits these 3 points and take its radius
           * as an estimate of the curvature at this point.
           */
          double x3 = log(gsl_vector_get(rho, i + 1));
          double y3 = log(gsl_vector_get(eta, i + 1));

          double x21 = x2 - x1;
          double y21 = y2 - y1;
          double x31 = x3 - x1;
          double y31 = y3 - y1;
          double h21 = x21*x21 + y21*y21;
          double h31 = x31*x31 + y31*y31;
          double d = fabs(2.0 * (x21*y31 - x31*y21));
          double r = sqrt(h21*h31*((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2))) / d;

          /* if d =~ 0 then there are nearly colinear points */
          if (gsl_finite(r))
            {
              /* check for smallest radius of curvature */
              if (r < rmin || rmin < 0.0)
                {
                  rmin = r;
                  *idx = i;
                }
            }

          /* update previous/current L-curve values */
          x1 = x2;
          y1 = y2;
          x2 = x3;
          y2 = y3;
        }

      /* check if a minimum radius was found */
      if (rmin < 0.0)
        {
          /* possibly co-linear points */
          GSL_ERROR("failed to find minimum radius", GSL_EINVAL);
        }

      return s;
    }
} /* gsl_multifit_linear_lcorner() */

/*
gsl_multifit_linear_lcorner2()
  Determine point on L-curve (lambda^2, ||c||^2) of maximum curvature.
For each set of 3 points on the L-curve, the circle which passes through
the 3 points is computed. The radius of the circle is then used
as an estimate of the curvature at the middle point. The point
with maximum curvature is then selected.

This routine is based on the paper

M. Rezghi and S. M. Hosseini, "A new variant of L-curve for Tikhonov
regularization", J. Comp. App. Math., 231 (2009).

Inputs: reg_param - vector of regularization parameters
        eta       - vector of solution norms ||L x||
        idx       - (output) index i such that
                    (lambda(i)^2,eta(i)^2) is the point of
                    maximum curvature

Return: success/error
*/

int
gsl_multifit_linear_lcorner2(const gsl_vector *reg_param,
                             const gsl_vector *eta,
                             size_t *idx)
{
  const size_t n = reg_param->size;

  if (n < 3)
    {
      GSL_ERROR ("at least 3 points are needed for L-curve analysis",
                 GSL_EBADLEN);
    }
  else if (n != eta->size)
    {
      GSL_ERROR ("size of reg_param and eta vectors do not match",
                 GSL_EBADLEN);
    }
  else
    {
      int s = GSL_SUCCESS;
      size_t i;
      double x1, y1;      /* first point of triangle on L-curve */
      double x2, y2;      /* second point of triangle on L-curve */
      double rmin = -1.0; /* minimum radius of curvature */

      /* initial values */
      x1 = gsl_vector_get(reg_param, 0);
      x1 *= x1;
      y1 = gsl_vector_get(eta, 0);
      y1 *= y1;

      x2 = gsl_vector_get(reg_param, 1);
      x2 *= x2;
      y2 = gsl_vector_get(eta, 1);
      y2 *= y2;

      for (i = 1; i < n - 1; ++i)
        {
          /*
           * The points (x1,y1), (x2,y2), (x3,y3) are the previous,
           * current, and next point on the L-curve. We will find
           * the circle which fits these 3 points and take its radius
           * as an estimate of the curvature at this point.
           */
          double lamip1 = gsl_vector_get(reg_param, i + 1);
          double etaip1 = gsl_vector_get(eta, i + 1);
          double x3 = lamip1 * lamip1;
          double y3 = etaip1 * etaip1;

          double x21 = x2 - x1;
          double y21 = y2 - y1;
          double x31 = x3 - x1;
          double y31 = y3 - y1;
          double h21 = x21*x21 + y21*y21;
          double h31 = x31*x31 + y31*y31;
          double d = fabs(2.0 * (x21*y31 - x31*y21));
          double r = sqrt(h21*h31*((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2))) / d;

          /* if d =~ 0 then there are nearly colinear points */
          if (gsl_finite(r))
            {
              /* check for smallest radius of curvature */
              if (r < rmin || rmin < 0.0)
                {
                  rmin = r;
                  *idx = i;
                }
            }

          /* update previous/current L-curve values */
          x1 = x2;
          y1 = y2;
          x2 = x3;
          y2 = y3;
        }

      /* check if a minimum radius was found */
      if (rmin < 0.0)
        {
          /* possibly co-linear points */
          GSL_ERROR("failed to find minimum radius", GSL_EINVAL);
        }

      return s;
    }
} /* gsl_multifit_linear_lcorner2() */

#define GSL_MULTIFIT_MAXK      100

/*
gsl_multifit_linear_L()
  Compute discrete approximation to derivative operator of order
k on a regular grid of p points, ie: L is (p-k)-by-p
*/

int
gsl_multifit_linear_Lk(const size_t p, const size_t k, gsl_matrix *L)
{
  if (p <= k)
    {
      GSL_ERROR("p must be larger than derivative order", GSL_EBADLEN);
    }
  else if (k >= GSL_MULTIFIT_MAXK - 1)
    {
      GSL_ERROR("derivative order k too large", GSL_EBADLEN);
    }
  else if (p - k != L->size1 || p != L->size2)
    {
      GSL_ERROR("L matrix must be (p-k)-by-p", GSL_EBADLEN);
    }
  else
    {
      double c_data[GSL_MULTIFIT_MAXK];
      gsl_vector_view cv = gsl_vector_view_array(c_data, k + 1);
      size_t i, j;

      /* zeroth derivative */
      if (k == 0)
        {
          gsl_matrix_set_identity(L);
          return GSL_SUCCESS;
        }

      gsl_matrix_set_zero(L);

      gsl_vector_set_zero(&cv.vector);
      gsl_vector_set(&cv.vector, 0, -1.0);
      gsl_vector_set(&cv.vector, 1, 1.0);

      for (i = 1; i < k; ++i)
        {
          double cjm1 = 0.0;

          for (j = 0; j < k + 1; ++j)
            {
              double cj = gsl_vector_get(&cv.vector, j);

              gsl_vector_set(&cv.vector, j, cjm1 - cj);
              cjm1 = cj;
            }
        }

      /* build L, the c_i are the entries on the diagonals */
      for (i = 0; i < k + 1; ++i)
        {
          gsl_vector_view v = gsl_matrix_superdiagonal(L, i);
          double ci = gsl_vector_get(&cv.vector, i);

          gsl_vector_set_all(&v.vector, ci);
        }

      return GSL_SUCCESS;
    }
} /* gsl_multifit_linear_Lk() */

/*
gsl_multifit_linear_Lsobolev()
  Construct Sobolev smoothing norm operator

L = [ a_0 I; a_1 L_1; a_2 L_2; ...; a_k L_k ]

by computing the Cholesky factor of L^T L

Inputs: p     - number of columns of L
        kmax  - maximum derivative order (< p)
        alpha - vector of weights; alpha_k multiplies L_k, size kmax + 1
        L     - (output) upper triangular Sobolev matrix p-by-p,
                stored in upper triangle
        work  - workspace

Notes:
1) work->Q is used to store intermediate L_k matrices
*/

int
gsl_multifit_linear_Lsobolev(const size_t p, const size_t kmax,
                             const gsl_vector *alpha, gsl_matrix *L,
                             gsl_multifit_linear_workspace *work)
{
  if (p > work->pmax)
    {
      GSL_ERROR("p is larger than workspace", GSL_EBADLEN);
    }
  else if (p <= kmax)
    {
      GSL_ERROR("p must be larger than derivative order", GSL_EBADLEN);
    }
  else if (kmax + 1 != alpha->size)
    {
      GSL_ERROR("alpha must be size kmax + 1", GSL_EBADLEN);
    }
  else if (p != L->size1)
    {
      GSL_ERROR("L matrix is wrong size", GSL_EBADLEN);
    }
  else if (L->size1 != L->size2)
    {
      GSL_ERROR("L matrix is not square", GSL_ENOTSQR);
    }
  else
    {
      int s;
      size_t j, k;
      gsl_vector_view d = gsl_matrix_diagonal(L);
      const double alpha0 = gsl_vector_get(alpha, 0);

      /* initialize L to alpha0^2 I */
      gsl_matrix_set_zero(L);
      gsl_vector_add_constant(&d.vector, alpha0 * alpha0);

      for (k = 1; k <= kmax; ++k)
        {
          gsl_matrix_view Lk = gsl_matrix_submatrix(work->Q, 0, 0, p - k, p);
          double ak = gsl_vector_get(alpha, k);

          /* compute a_k L_k */
          s = gsl_multifit_linear_Lk(p, k, &Lk.matrix);
          if (s)
            return s;

          gsl_matrix_scale(&Lk.matrix, ak);

          /* LTL += L_k^T L_k */
          gsl_blas_dsyrk(CblasLower, CblasTrans, 1.0, &Lk.matrix, 1.0, L);
        }

      s = gsl_linalg_cholesky_decomp(L);
      if (s)
        return s;

      /* copy Cholesky factor to upper triangle and zero out bottom */
      gsl_matrix_transpose_tricpy('L', 1, L, L);

      for (j = 0; j < p; ++j)
        {
          for (k = 0; k < j; ++k)
            gsl_matrix_set(L, j, k, 0.0);
        }

      return GSL_SUCCESS;
    }
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <cblas.h>

#include <aml/defs.h>
#include <aml/handle.h>

namespace aml {
namespace impl {
namespace cpu {

inline CBLAS_TRANSPOSE convert_op(char op) {
  if (op == 'n') {
    return CblasNoTrans;
  } else {
    return CblasTrans;
  }
}

/** BLAS LEVEL 1 **************************************************************/

inline float nrm2(aml::Handle h,
                  Index n,
                  const float *x) {
  AML_ASSERT_INT(n);
  auto tic = h.tic("cpu_snrm2_" + std::to_string(n));

  float res = cblas_snrm2(n, x, 1);

  tic.stop();
  return res;
}

inline double nrm2(aml::Handle h,
                   Index n,
                   const double *x) {
  AML_ASSERT_INT(n);
  auto tic = h.tic("cpu_dnrm2_" + std::to_string(n));

  double res = cblas_dnrm2(n, x, 1);

  tic.stop();
  return res;
}

/** BLAS LEVEL 2 **************************************************************/

inline void gemv(aml::Handle h,
                 char op,
                 Index m,
                 Index n,
                 float alpha,
                 const float *a,
                 Index lda,
                 const float *x,
                 float beta,
                 float *y) {
  AML_ASSERT_INT(m, n, lda);
  auto tic = h.tic("cpu_sgemv_" + std::to_string(m) + "_" + std::to_string(n));

  cblas_sgemv(CblasColMajor,
      convert_op(op), m, n, alpha, a, lda, x, 1, beta, y, 1);

  tic.stop();
}

inline void gemv(aml::Handle h,
                 char op,
                 Index m,
                 Index n,
                 double alpha,
                 const double *a,
                 Index lda,
                 const double *x,
                 double beta,
                 double *y) {
  AML_ASSERT_INT(m, n, lda);
  auto tic = h.tic("cpu_dgemv_" + std::to_string(m) + "_" + std::to_string(n));

  cblas_dgemv(CblasColMajor,
      convert_op(op), m, n, alpha, a, lda, x, 1, beta, y, 1);

  tic.stop();
}

inline void trsv(aml::Handle h,
                 char op,
                 Index m,
                 const float *a,
                 Index lda,
                 float *x) {
  AML_ASSERT_INT(m, lda);
  auto tic = h.tic("cpu_strsv_" + std::to_string(m));

  cblas_strsv(CblasColMajor, CblasLower,
      convert_op(op), CblasNonUnit, m, a, lda, x, 1);

  tic.stop();
}

inline void trsv(aml::Handle h,
                 char op,
                 Index m,
                 const double *a,
                 Index lda,
                 double *x) {
  AML_ASSERT_INT(m, lda);

  auto tic = h.tic("cpu_dtrsv_" + std::to_string(m));

  cblas_dtrsv(CblasColMajor, CblasLower,
      convert_op(op), CblasNonUnit, m, a, lda, x, 1);

  tic.stop();
}

/** BLAS LEVEL 3 **************************************************************/

inline void gemm(aml::Handle h,
                 char op_a,
                 char op_b,
                 Index m,
                 Index n,
                 Index k,
                 float alpha,
                 const float *a,
                 Index lda,
                 const float *b,
                 Index ldb,
                 float beta,
                 float *c,
                 Index ldc) {
  AML_ASSERT_INT(m, n, k, lda, ldb, ldc);
  auto tic = h.tic("cpu_sgemm_"
      + std::to_string(m) + "_" + std::to_string(n) + "_" + std::to_string(k));

  cblas_sgemm(CblasColMajor, convert_op(op_a), convert_op(op_b),
      m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);

  tic.stop();
}

inline void gemm(aml::Handle h,
                 char op_a,
                 char op_b,
                 Index m,
                 Index n,
                 Index k,
                 double alpha,
                 const double *a,
                 Index lda,
                 const double *b,
                 Index ldb,
                 double beta,
                 double *c,
                 Index ldc) {
  AML_ASSERT_INT(m, n, k, lda, ldb, ldc);
  auto tic = h.tic("cpu_dgemm_"
      + std::to_string(m) + "_" + std::to_string(n) + "_" + std::to_string(k));

  cblas_dgemm(CblasColMajor, convert_op(op_a), convert_op(op_b),
      m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);

  tic.stop();
}

inline void syrk(aml::Handle h,
                 char op,
                 Index n,
                 Index k,
                 float alpha,
                 const float *a,
                 Index lda,
                 float beta,
                 float *c,
                 Index ldc) {
  AML_ASSERT_INT(n, k, lda, ldc);
  auto tic = h.tic("cpu_ssyrk_" + std::to_string(n) + "_" + std::to_string(k));

  cblas_ssyrk(CblasColMajor, CblasLower,
      convert_op(op), n, k, alpha, a, lda, beta, c, ldc);

  tic.stop();
}

inline void syrk(aml::Handle h,
                 char op,
                 Index n,
                 Index k,
                 double alpha,
                 const double *a,
                 Index lda,
                 double beta,
                 double *c,
                 Index ldc) {
  AML_ASSERT_INT(n, k, lda, ldc);
  auto tic = h.tic("cpu_dsyrk_" + std::to_string(n) + "_" + std::to_string(k));

  cblas_dsyrk(CblasColMajor, CblasLower,
      convert_op(op), n, k, alpha, a, lda, beta, c, ldc);

  tic.stop();
}

inline void trsm(aml::Handle h,
                 char op,
                 Index m,
                 Index n,
                 float alpha,
                 const float *a,
                 Index lda,
                 float *b,
                 Index ldb) {
  AML_ASSERT_INT(m, n, lda, ldb);
  auto tic = h.tic("cpu_strsm_" + std::to_string(m) + "_" + std::to_string(n));

  cblas_strsm(
      CblasColMajor, CblasRight, CblasLower, convert_op(op), CblasNonUnit,
      m, n, alpha, a, lda, b, ldb);

  tic.stop();
}

inline void trsm(aml::Handle h,
                 char op,
                 Index m,
                 Index n,
                 double alpha,
                 const double *a,
                 Index lda,
                 double *b,
                 Index ldb) {
  AML_ASSERT_INT(m, n, lda, ldb);
  auto tic = h.tic("cpu_dtrsm_" + std::to_string(m) + "_" + std::to_string(n));

  cblas_dtrsm(
      CblasColMajor, CblasRight, CblasLower, convert_op(op), CblasNonUnit,
      m, n, alpha, a, lda, b, ldb);

  tic.stop();
}

}  // namespace cpu
}  // namespace impl
}  // namespace aml

max_stars_repo_name: foges/aml

max_issues_repo_name: foges/aml

max_forks_repo_name: foges/aml


text:
/*! \file deleter.h
    \brief gsl_interp_accelとgsl_splineのデリータを宣言・定義したヘッダファイル

    Copyright © 2015 @dc1394 All Rights Reserved.
    This software is released under the BSD 2-Clause License.
*/

#ifndef _DELETER_H_
#define _DELETER_H_

#pragma once

#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

namespace getdata {
    //! A function.
    /*!
        gsl_interp_accelへのポインタを解放するラムダ式
        \param acc gsl_interp_accelへのポインタ
    */
    static auto const gsl_interp_accel_deleter = [](gsl_interp_accel * acc) {
        gsl_interp_accel_free(acc);
    };

    //! A function.
    /*!
        gsl_splineへのポインタを解放するラムダ式
        \param spline gsl_splineへのポインタ
    */
    static auto const gsl_spline_deleter = [](gsl_spline * spline) {
        gsl_spline_free(spline);
    };
}

#endif  // _DELETER_H_


text:
 // MIT License

 // Copyright (c) [2017] [Vinay Yuvashankar]

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

 // The above copyright notice and this permission notice shall be included in all
 // copies or substantial portions of the Software.

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

/** \file ERSP.cc

	\brief This file contains all of the function required to generate the Event Related Spectral Pertubation of EEG signals.
*/
#include "wavelet.h"
#include <gsl/gsl_statistics.h>

#define NORMALIZATION_FACTOR 0.375402913609157562

int ERSP (double * raw_data, double* scales, const int sampling_frequency, const int n,
	const int J, int const trials, const int padding_type,
	double * output)
{
	int i, j, x;

	//Calculate the necessary constants for the Continuous Wavelet Transform.
    const int    PADDED_SIZE = CalculatePaddingSize(n, padding_type);
    const int    m           = PRE_EVENT_TIME * sampling_frequency;
    const double dw          = (2 * M_PI * sampling_frequency)/(PADDED_SIZE); //NOT IN RAD/SEC in Hz

    fftw_init_threads();
    #pragma omp parallel private(i, j, x) shared(raw_data, output, scales) default(none)
    {
    	//Array Inits
    	double * pre_stimulus, *wavelet_out, *baseline_out;
    	fftw_complex *data_in, *fft_data, *filter_convolution, *fftw_result;
    	fftw_plan plan_forward, plan_backward;

	    //Memory Allocations
	    wavelet_out  = (double*) malloc( n * J * sizeof(double) );
	    baseline_out = (double*) malloc( n * J * sizeof(double) );
	    pre_stimulus = (double*) malloc( m     * sizeof(double) );

	    //FFTW Memory Allocations
	    data_in            = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * PADDED_SIZE );
		fft_data           = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * PADDED_SIZE );
		filter_convolution = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * PADDED_SIZE );
		fftw_result        = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * PADDED_SIZE );

		#pragma omp critical (make_plan)
		{
			fftw_plan_with_nthreads(1);
			if (fftw_import_wisdom_from_filename("FFTW_Plan.wise") == 0)
			{
				printf("No FFTW Plan, using an unoptimized method\n");
				plan_forward  = fftw_plan_dft_1d(PADDED_SIZE, data_in,            fft_data,    FFTW_FORWARD,  FFTW_ESTIMATE);
				plan_backward = fftw_plan_dft_1d(PADDED_SIZE, filter_convolution, fftw_result, FFTW_BACKWARD, FFTW_ESTIMATE);
			}
			else
			{
				plan_forward  = fftw_plan_dft_1d(PADDED_SIZE, data_in,            fft_data,    FFTW_FORWARD,  FFTW_EXHAUSTIVE);
				plan_backward = fftw_plan_dft_1d(PADDED_SIZE, filter_convolution, fftw_result, FFTW_BACKWARD, FFTW_EXHAUSTIVE);
			}
		}

		/*Begin ERSP*/
		#pragma omp for
		for ( x = 0; x < trials; ++x)
		{
			memset(data_in,            0.0, sizeof( fftw_complex ) * PADDED_SIZE);
			memset(fft_data,           0.0, sizeof( fftw_complex ) * PADDED_SIZE);
			memset(filter_convolution, 0.0, sizeof( fftw_complex ) * PADDED_SIZE);
			memset(fftw_result,        0.0, sizeof( fftw_complex ) * PADDED_SIZE);

			/*Begin Wavelet Analysis*/
			PopulateDataArray(raw_data, n, x,
							  PADDED_SIZE, padding_type, data_in);
			fftw_execute(plan_forward);

			for (i = 0; i < J; ++i)
			{
				FrequencyMultiply(fft_data, PADDED_SIZE, scales[i], dw,
								  filter_convolution);

				//Take the inverse FFT and store it in fftw_result
				fftw_execute(plan_backward);

				//Calculate the power and store it in result this may need to be changed to accomodate for phase
				for (j = 0; j < n; ++j)
				{
					wavelet_out[i * n + j] = MAGNITUDE(fftw_result[j][0], fftw_result[j][1]) / (NORMALIZATION_FACTOR * sqrt(scales[i]) );
				}
			}
			/*End Wavelet Analysis*/

			//Remove the baseline
			RemoveBaseline(pre_stimulus, wavelet_out,
						   n, J, m,
						   baseline_out);

			for ( i = 0; i < n * J; ++i)
			{
				output[i] += fabs(baseline_out[i]);
			}
		}

		#pragma omp for simd
		for ( i = 0; i < n * J; ++i)
		{
			output[i] /= trials;
		}
		/*End ERSP*/

		//Sanitation Engineering
		fftw_destroy_plan(plan_forward); fftw_destroy_plan(plan_backward);
		fftw_free(data_in); fftw_free(fft_data); fftw_free(filter_convolution); fftw_free(fftw_result);
		free(pre_stimulus); free(baseline_out); free(wavelet_out);

    }/*End of OpenMP*/
	return(0);
}

int RemoveBaseline(double* pre_stimulus, double* pre_baseline_array,
	const int n, const int J, const int m,
	double* output)
{
	double value;
	double mean, sDeviation;
	const int stride = 1;
	int i, j;

	for ( i = 0; i < J; ++i)
	{
		//Copy the pre trial results from each frequency block into pre_stimulus.
		for ( j = 0; j < m; ++j)
		{
			pre_stimulus[j] = pre_baseline_array[i * n + j];
		}

		//Calculate mean and standard deviation
		mean       = gsl_stats_mean(pre_stimulus, stride, m);
    	sDeviation = gsl_stats_sd_m(pre_stimulus, stride, m, mean);

    	//Remove the Baseline
	    for ( j = 0; j < n; ++j)
	    {
	    	value = pre_baseline_array[i * n + j] * pre_baseline_array[i * n + j];
	        output[i * n + j] = (fabs(value) - mean) / sDeviation;
	        // output[i * n + j] = ( pre_baseline_array[i * n + j] - mean ) / sDeviation;
	    }
	}

	return(0);
}

int Generate_FFTW_Wisdom(int padded_size)
{
	int success_flag = 1;

	//Array Inits
	fftw_complex *data_in, *fft_data, *filter_convolution, *fftw_result;
	fftw_plan plan_forward, plan_backward;

	fftw_init_threads();

	//FFTW Memory Allocations
    data_in            = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * padded_size );
	fft_data           = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * padded_size );
	filter_convolution = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * padded_size );
	fftw_result        = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * padded_size );

	//FFTW Planning
	fftw_plan_with_nthreads(1);
	printf("Generating an Exhaustive FFTW Plan\n");
	plan_forward  = fftw_plan_dft_1d(padded_size, data_in,            fft_data,    FFTW_FORWARD,  FFTW_EXHAUSTIVE);
	plan_backward = fftw_plan_dft_1d(padded_size, filter_convolution, fftw_result, FFTW_BACKWARD, FFTW_EXHAUSTIVE);

	printf("Writing FFTW plan to FFTW_Plan.wise\n");
	if (fftw_export_wisdom_to_filename("FFTW_Plan.wise") != 0)
	{
		success_flag = 0;
	}

	fftw_destroy_plan(plan_forward); fftw_destroy_plan(plan_backward);
	fftw_free(data_in); fftw_free(fft_data); fftw_free(filter_convolution); fftw_free(fftw_result);

	return(success_flag);
}

int FrequencyMultiply(fftw_complex* fft_data,
	const int data_size, const double scale, const double dw,
	fftw_complex* filter_convolution)
{
	int j;
	double value;
	//Compute the Fourier Morlet at 0 and N/2
	double norm = sqrt(scale);
	value = CompleteFourierMorlet(0.0, scale, norm);

	filter_convolution[0][0] = (fft_data[0][0]/data_size) * value;
	filter_convolution[0][1] = (fft_data[0][1]/data_size) * value;

	filter_convolution[data_size/2][0] = 0.0;
	filter_convolution[data_size/2][1] = 0.0;

	//Compute the Fourier Morlet Convolution in between
	for (j = 1; j < data_size/2 - 1; ++j)
	{
		value = CompleteFourierMorlet( j * dw , scale, norm);
		filter_convolution[j][0] = (fft_data[j][0]/data_size) * value;
		filter_convolution[j][1] = (fft_data[j][1]/data_size) * value;

		filter_convolution[data_size- j][0] = 0.0;
		filter_convolution[data_size- j][1] = 0.0;
	}

	return(0);
}

int PopulateDataArray(double* input_data, const int data_size, const int trial_number,
	const int padded_size, const int padding_type,
	fftw_complex* output_data)
{
	const double ramp = 2.0/data_size; // = 1.0/ n / 2
	double gain;
	int i;

	int output_counter = 0;
	int input_counter = 0;

	switch(padding_type)
	{
		case 0: //No Padding what so ever
			//populate the FFTW data vector.
			for (i = 0; i < data_size; ++i)
		    {
		    	output_data[i][0] = input_data[trial_number * data_size + i];
		    	output_data[i][1] = 0.0;
		    }
		    break;

		case 1: //Zero - Padding
			//populate the FFTW data vector.
			for (i = 0; i < data_size; ++i)
		    {
		    	output_data[i][0] = input_data[trial_number * data_size + i];
		    	output_data[i][1] = 0.0;
		    }

		    //Horse the rest of the data vector to zero just in case
		    for (i = data_size; i < padded_size; ++i)
		    {
		    	output_data[i][0] = 0.0;
		    	output_data[i][1] = 0.0;
		    }
		    break;

		case 2: //Duplicate array and ramp up and ramp down output
			for (i = 0; i < data_size; ++i)
		    {
		    	output_data[i][0] = input_data[trial_number * data_size + i];
		    	output_data[i][1] = 0.0;
		    }

		    for (i = 0; i < (int) (0.5* data_size); ++i)
		    {
		    	output_counter = data_size + i;
		    	input_counter = (int) (trial_number * data_size + (0.5*data_size) + i);
		    	gain = i * ramp;

		    	output_data[output_counter][0] = (1.0 - gain) * input_data[input_counter];
		    	output_data[output_counter][1] = 0.0;

		    	output_data[output_counter + (int) (0.5 * data_size)][0] = gain * input_data[trial_number * data_size + i];
		    	output_data[output_counter + (int) (0.5 * data_size)][1] = 0.0;
		    }
			break;

		default: //Return just the array no padding.
			//populate the FFTW data vector.
			for (i = 0; i < data_size; ++i)
		    {
		    	output_data[i][0] = input_data[trial_number * data_size + i];
		    	output_data[i][1] = 0.0;
		    }
		    break;
	}
    return(padding_type);
}

max_stars_repo_name: yuvashankar/Research

max_issues_repo_name: yuvashankar/Research

max_forks_repo_name: yuvashankar/Research


text:
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

typedef struct ornsteinuhlenbeckprocess_{
  int steps;
  double dt;
  double sigma;
  double *x;
  gsl_rng *rng;
} OrnsteinUhlenbeckProcess;

void OrnsteinUhlenbeckProcess_run(OrnsteinUhlenbeckProcess *self) {
  int i;
  double sgm = self->sigma * sqrt(self->dt);
  for (i = 1; i < self->steps; ++i){
    self->x[i] = (1 - self->dt) * self->x[i-1] \
      + gsl_ran_gaussian_ziggurat(self->rng, sgm);
  }
}

max_stars_repo_name: tkf/comparatist

max_issues_repo_name: tkf/comparatist

max_forks_repo_name: tkf/comparatist


text:
#pragma once
#include <gsl/gsl>
#include <d3d11.h>
#undef min
#undef max

namespace Halley {
	class DX11Video;

	class DX11Buffer {
    public:
		enum class Type
		{
			Vertex,
			Index,
			Constant
		};

        DX11Buffer(DX11Video& video, Type type, size_t initialSize = 0);
		DX11Buffer(DX11Buffer&& other) noexcept;
		~DX11Buffer();

		DX11Buffer(const DX11Buffer& other) = delete;

		DX11Buffer& operator=(const DX11Buffer& other) = delete;
		DX11Buffer& operator=(DX11Buffer&& other) = delete;

		void setData(gsl::span<const gsl::byte> data);
		ID3D11Buffer*& getBuffer();
		UINT getOffset() const;
		UINT getLastSize() const;
		bool canFit(size_t size) const;

		void reset();
		void clear();

	private:
		DX11Video& video;
		Type type;
		ID3D11Buffer* buffer = nullptr;
		size_t curSize = 0;
		size_t curPos = 0;
		size_t lastSize = 0;
		size_t lastPos = 0;
		bool waitingReset = false;

		void resize(size_t size);
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


text:
#include <math.h>
#include <gsl/gsl_nan.h>

#include "csm_all.h"

#include <egsl/egsl_macros.h>

void find_neighbours(LDP ld, int i, int max_num, int*indexes, size_t*num_found);
void filter_orientation(double theta0, double rho0, size_t n,
 	const double*thetas, const double*rhos, double *alpha, double*cov0_alpha );

/** Requires the "cluster" field to be set */
void ld_compute_orientation(LDP ld, int size_neighbourhood, double sigma) {
	int i;
	for(i=0;i<ld->nrays;i++){
		if(!ld_valid_ray(ld,i) || (ld->cluster[i] == -1)) {
			ld->alpha[i] = GSL_NAN;
			ld->cov_alpha[i] = GSL_NAN;
			ld->alpha_valid[i] = 0;
			continue;
		}

		int neighbours[size_neighbourhood*2];
		size_t num_neighbours;
		find_neighbours(ld, i, size_neighbourhood, neighbours, &num_neighbours);

		if(0==num_neighbours) {
			ld->alpha[i] = GSL_NAN;
			ld->cov_alpha[i] = GSL_NAN;
			ld->alpha_valid[i] = 0;
			continue;
		}

/*		printf("orientation for i=%d:\n",i); */
		double thetas[num_neighbours];
		double readings[num_neighbours];
		size_t a=0;
		for(a=0;a<num_neighbours;a++) {
			thetas[a] = ld->theta[neighbours[a]];
			readings[a] = ld->readings[neighbours[a]];
			/* printf(" j = %d theta = %f rho = %f\n", neighbours[a], thetas[a],readings[a]); */
		}

		double alpha=42, cov0_alpha=32;
		filter_orientation(ld->theta[i],ld->readings[i],num_neighbours,
			thetas,readings,&alpha,&cov0_alpha);

		if(gsl_isnan(alpha)) {
			ld->alpha[i] = GSL_NAN;
			ld->cov_alpha[i] = GSL_NAN;
			ld->alpha_valid[i] = 0;
		} else {
			ld->alpha[i] = alpha;
			ld->cov_alpha[i] = cov0_alpha * square(sigma);
			ld->alpha_valid[i] = 1;
		}
		/* printf("---------- i = %d alpha = %f sigma=%f cov_alpha = %f\n", i, alpha, ld->cov_alpha[i]);*/
	}
}

/** A very cool algorithm for finding the orientation */

void filter_orientation(double theta0, double rho0, size_t n,
 	const double*thetas, const double*rhos, double *alpha, double*cov0_alpha ) {

	egsl_push();
	/* Y = L x + R epsilon */
	val Y = zeros(n,1);
	val L = ones(n,1);
	val R = zeros(n,n+1);

 	size_t i; for(i=0;i<n;i++) {
		*egsl_atmp(Y, i, 0)   = (rhos[i]-rho0)/(thetas[i]-theta0);
		*egsl_atmp(R, i, 0)   =             -1/(thetas[i]-theta0);
		*egsl_atmp(R, i, i+1) =             +1/(thetas[i]-theta0);
	}

	val eRinv = inv(m(R, tr(R)));
	val vcov_f1 = inv(m3(tr(L),eRinv, L));
	val vf1 =   m4(vcov_f1, tr(L), eRinv, Y);

	double cov_f1 = *egsl_atmp(vcov_f1,0,0);
	double f1 = *egsl_atmp(vf1,0,0);

	*alpha = theta0 - atan(f1/rho0);

	if(cos(*alpha)*cos(theta0)+sin(*alpha)*sin(theta0)>0)
		*alpha = *alpha + M_PI;

	double dalpha_df1  = rho0 / (square(rho0) + square(f1));
	double dalpha_drho = -f1 /  (square(rho0) + square(f1));

	*cov0_alpha	= square(dalpha_df1) * cov_f1 + square(dalpha_drho);

	if(gsl_isnan(*alpha)) {
		egsl_print("Y",Y);
		egsl_print("L",L);
		egsl_print("R",R);
		egsl_print("eRinv",eRinv);
		egsl_print("vcov_f1",vcov_f1);

		printf("   f1 = %f cov =%f \n", f1,cov_f1);
		printf("   f1/rho = %f \n", f1/rho0);
		printf("   atan = %f \n", atan(f1/rho0));
		printf("   theta0= %f \n", theta0);
	}

	egsl_pop();
/*
//	printf("dalpha_df1 = %f dalpha_drho = %f\n",dalpha_df1,dalpha_drho);
//	printf("f1 = %f covf1 = %f alpha = %f cov_alpha = %f\n ",f1,cov_f1,*alpha,*cov0_alpha);
//	printf("sotto = %f\n ",(square(rho0) + square(f1)));

//	printf("   alpha = %f sigma= %f°\n", *alpha, rad2deg(0.01*sqrt(*cov0_alpha)));

	printf("l= ");
	gsl_matrix_fprintf(stdout, l, "%f");
	printf("\ny= ");
	gsl_matrix_fprintf(stdout, y, "%f");
	printf("\nr= ");
	gsl_matrix_fprintf(stdout, r, "%f");
	printf("\ninv(r*r)= ");
	gsl_matrix_fprintf(stdout, Rinv, "%f");
	printf("\nf1 = %lf ",f1);
	printf("\ncov_f1 = %lf ",cov_f1);
*/
}

/* indexes: an array of size "max_num*2" */
void find_neighbours(LDP ld, int i, int max_num, int*indexes, size_t*num_found) {
	*num_found = 0;

	int up = i;
	while ((up+1 <= i+max_num) && (up+1<ld->nrays) && ld_valid_ray(ld,up+1)
			&& (ld->cluster[up+1] == ld->cluster[i])) {
		up+=1;
		indexes[(*num_found)++] = up;
	}
	int down = i;
	while ((down >= i-max_num) && (down-1>=0) && ld_valid_ray(ld,down-1) &&
			(ld->cluster[down-1] == ld->cluster[i])) {
		down-=1;
		indexes[(*num_found)++] = down;
	}
}


text:
#ifndef SOLVER_H
#define SOLVER_H

#if PY_MAJOR_VERSION >= 3
#define PY3K
#endif

// include user-defined headers
#include "common.h"
#include "utils.h"

#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>

typedef struct _TOV_Tuple {
    double a;
    double b;
    double c;
} Pair;

#endif


text:
/* examples/C/ssmfe/precond_ssmfe.c */
/* Laplacian on a square grid (using SPRAL_SSMFE routines) */
#include "spral.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <cblas.h>

/* Header that implements Laplacian and preconditioners */
#include "laplace2d.h"

int main(void) {
   const int m   = 20;     /* grid points along each side */
   const int n   = m*m;    /* problem size */
   const int nep = 5;      /* eigenpairs wanted */

   double lambda[2*nep];                  /* eigenvalues */
   double X[2*nep][n];                    /* eigenvectors */
   struct spral_ssmfe_rcid rci;           /* reverse communication data */
   struct spral_ssmfe_options options;    /* options */
   void *keep;                            /* private data */
   struct spral_ssmfe_inform inform;      /* information */

   /* Initialize options to default values */
   spral_ssmfe_default_options(&options);
   /* gap between the last converged eigenvalue and the rest of the spectrum
    * must be at least 0.1 times average gap between computed eigenvalues */
   options.left_gap = -0.1;

   rci.job = 0; keep = NULL;
   while(true) { /* reverse communication loop */
      spral_ssmfe_standard_double(&rci, nep, 2*nep, lambda, n, &X[0][0], n,
         &keep, &options, &inform);
      switch ( rci.job ) {
      case 1:
         apply_laplacian(m, m, rci.nx, rci.x, rci.y);
         break;
      case 2:
         apply_gauss_seidel_step(m, m, rci.nx, rci.x, rci.y);
         break;
      default:
         goto finished;
      }
   }
finished:
   printf("%d eigenpairs converged in %d iterations\n", inform.left, inform.iteration);
   for(int i=0; i<inform.left; i++)
      printf(" lambda[%1d] = %13.7e\n", i, lambda[i]);
   spral_ssmfe_free_double(&keep, &inform);

   /* Success */
   return 0;
}

max_stars_repo_name: mjacobse/spral

max_issues_repo_name: mjacobse/spral

max_forks_repo_name: mjacobse/spral


text:
/* This module contains the routine rotCurveSpline, which takes a
   tabulated rotaiton curve and calculates a B-spline fit to it and
   its derivatives. */

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_multifit.h>
#include "rotCurveSpline.h"

struct dpsidr_params {
  gsl_bspline_workspace *wksp;
#ifdef GSL1
  gsl_bspline_deriv_workspace *dwksp;
#endif
  gsl_vector *B, *coef;
  gsl_matrix *dB, *cov;
};

double dpsidr(double r, void *params) {
  double vphi, err;
  struct dpsidr_params *p = (struct dpsidr_params *) params;

#ifdef GSL1
  gsl_bspline_deriv_eval(log(r), 2, p->dB, p->wksp, p->dwksp);
#else
  gsl_bspline_deriv_eval(log(r), 2, p->dB, p->wksp);
#endif
  gsl_matrix_get_col(p->B, p->dB, 0);
  gsl_multifit_linear_est(p->B, p->coef, p->cov, &vphi, &err);
  return(-vphi*vphi/r);
}

void
rotCurveSpline(const double *rTab, const double *vTab,
	       const unsigned long nTab,
	       const unsigned long bspline_degree,
	       const unsigned long bspline_breakpoints,
	       const double *r, const unsigned long nr,
	       double *vphi, double *psi,
	       double *beta) {

  gsl_vector *xInput, *vInput, *coef, *B, *bkpt_vec;
  gsl_matrix *dB, *cov, *pred;
  gsl_bspline_workspace *wksp;
#ifdef GSL1
  gsl_bspline_deriv_workspace *dwksp;
#endif
  gsl_multifit_linear_workspace *mw;
  gsl_integration_workspace *intwksp;
  gsl_function dpsidr_func;
  struct dpsidr_params params;
  double x0, x1, v0, bkpt_sum_target, sum, x, Bij, chisq;
  double dvphidlogr, err, psiInt;
  unsigned long bspline_ncoef, i, j, ptr;

  /* Safety check: make sure we have enough data points for the choice
     of degree and breakpoints */
  bspline_ncoef = bspline_breakpoints + bspline_degree - 2;
  if (bspline_ncoef > nTab) {
    fprintf(stderr, "Error: %lu data points for %lu coefficients; ndata must be >= ncoef = nbreakpoints + degree - 2\n", nTab, bspline_ncoef);
    exit(1);
  }

  /* Initialize the GSL vectors to hold the data table */
  xInput = gsl_vector_alloc(nTab);
  vInput = gsl_vector_alloc(nTab);
  for (i=0; i<nTab; i++) {
    gsl_vector_set(xInput, i, log(rTab[i]));
    gsl_vector_set(vInput, i, vTab[i]);
  }

  /* Allocate workspace for and associated stuff for B-spline solver */
  wksp = gsl_bspline_alloc(bspline_degree, bspline_breakpoints);
#ifdef GSL1
  dwksp = gsl_bspline_deriv_alloc(bspline_degree);
#endif
  B = gsl_vector_alloc(bspline_ncoef);
  dB = gsl_matrix_alloc(bspline_ncoef, 3);

  /* Compute breakpoint locations using the method of Gans & Gill,
     Applied Spectroscopy, 38, 297, 1984. The formula for
     uniformly-spaced data is that two successive breakpoints t_i and t_i+1
     should be separated by an amount such that
     sum_{x=t_i}^{x=t_i+1} f(x)^0.5 dx = S_T
     where
     S_T = (1/(m+1)) sum f(x_i)^0.5 dx.
     Here the first sum runs over all data points between t_i and
     t_i+1, and the second runs all data points; m is the number of
     knots, and f(x) is the function to be approximated. This formula
     assumes uniform spacing dx, but here it is modified for
     non-uniform data spacing such that dx is taken to be the mean
     size of the data interval, averaging over the left and right
     sides.
  */
  v0 = gsl_vector_get(vInput, 0);
  x0 = gsl_vector_get(xInput, 0);
  x1 = gsl_vector_get(xInput, 1);
  bkpt_sum_target = sqrt(v0) * (x1-x0);
  for (i=1; i<nTab-1; i++) {
    v0 = gsl_vector_get(vInput, i);
    x0 = gsl_vector_get(xInput, i-1);
    x1 = gsl_vector_get(xInput, i+1);
    bkpt_sum_target += sqrt(v0) * 0.5 * (x1-x0);
  }
  v0 = gsl_vector_get(vInput, nTab-1);
  x0 = gsl_vector_get(xInput, nTab-2);
  x1 = gsl_vector_get(xInput, nTab-1);
  bkpt_sum_target += sqrt(v0) * (x1-x0);
  bkpt_sum_target /= (bspline_breakpoints + 1);
  bkpt_vec = gsl_vector_alloc(bspline_breakpoints);
  gsl_vector_set(bkpt_vec, 0, gsl_vector_get(xInput, 0));
  gsl_vector_set(bkpt_vec, bspline_breakpoints-1,
		 gsl_vector_get(xInput, nTab-1));
  ptr=0;
  for (i=1; i<bspline_breakpoints-2; i++) {
    sum = 0;
    while (sum < bkpt_sum_target) {
      if (ptr == 0) {
	v0 = gsl_vector_get(vInput, 0);
	x0 = gsl_vector_get(xInput, 0);
	x1 = gsl_vector_get(xInput, 1);
	sum += sqrt(v0) * (x1-x0);
      } else if (ptr == nTab-1) {
	v0 = gsl_vector_get(vInput, nTab-1);
	x0 = gsl_vector_get(xInput, nTab-2);
	x1 = gsl_vector_get(xInput, nTab-1);
	sum += sqrt(v0) * (x1-x0);
      } else {
	v0 = gsl_vector_get(vInput, ptr);
	x0 = gsl_vector_get(xInput, ptr-1);
	x1 = gsl_vector_get(xInput, ptr+1);
	sum += sqrt(v0) * 0.5 * (x1-x0);
      }
      ptr++;
      if (ptr == nTab) break;
    }
    if (ptr == nTab) break;
    gsl_vector_set(bkpt_vec, i, gsl_vector_get(xInput, ptr-1));
  }
  gsl_vector_set(bkpt_vec, i, gsl_vector_get(xInput, nTab-1));

  /* Construct knots vector */
  gsl_bspline_knots(bkpt_vec, wksp);

  /* Allocate workspace for the least squares fit to find the B spline
     coefficients; pred is the predictor variable matrix, cov is the
     covariance matrix, coef is the vector of best-fitting
     coefficients, and mw is the workspace used by the gsl fitter */
  pred = gsl_matrix_alloc(nTab, bspline_ncoef);
  mw = gsl_multifit_linear_alloc(nTab, bspline_ncoef);
  cov = gsl_matrix_alloc(bspline_ncoef, bspline_ncoef);
  coef = gsl_vector_alloc(bspline_ncoef);

  /* Load the predictor matrix with the B spline functions evaluated
     at the input positions */
  for (i=0; i<nTab; i++) {
    x = gsl_vector_get(xInput, i); /* Position of data point */
    gsl_bspline_eval(x, B, wksp); /* Evaluate B splines at this
				     position */
    /* Load predictor matrix */
    for (j=0; j<bspline_ncoef; j++) {
      Bij = gsl_vector_get(B, j);
      gsl_matrix_set(pred, i, j, Bij);
    }
  }

  /* Perform least squares fit to find coefficients */
  gsl_multifit_linear(pred, vInput, coef, cov, &chisq, mw);

  /* Compute B-spline fits at input positions */
  for (i=0; i<nr; i++) {
#ifdef GSL1
    gsl_bspline_deriv_eval(log(r[i]), 2, dB, wksp, dwksp);
#else
    gsl_bspline_deriv_eval(log(r[i]), 2, dB, wksp);
#endif
    gsl_matrix_get_col(B, dB, 0);
    gsl_multifit_linear_est(B, coef, cov, vphi+i, &err);
    gsl_matrix_get_col(B, dB, 1);
    gsl_multifit_linear_est(B, coef, cov, &dvphidlogr, &err);
    beta[i] = dvphidlogr / vphi[i];
  }

  /* Compute potential at cell centers and edges by numerically
     integrating from grid edge to center */
  intwksp = gsl_integration_workspace_alloc(4096);
  dpsidr_func.function = &dpsidr;
  dpsidr_func.params = &params;
  params.wksp = wksp;
#ifdef GSL1
  params.dwksp = dwksp;
#endif
  params.B = B;
  params.coef = coef;
  params.dB = dB;
  params.cov = cov;
  psi[nr-1] = 0.0;
  for (i=nr-1; i>0; i--) {
    gsl_integration_qag(&dpsidr_func, r[i], r[i-1],
			1e-6*vphi[i]*vphi[i]*
			(r[i+1]-r[i]),
			1e-6, 4096, 6, intwksp, &psiInt, &err);
    psi[i-1] = psi[i] - psiInt;
  }

  /* Free memory */
  gsl_matrix_free(pred);
  gsl_vector_free(bkpt_vec);
  gsl_vector_free(xInput);
  gsl_vector_free(vInput);
  gsl_multifit_linear_free(mw);
  gsl_vector_free(coef);
  gsl_vector_free(B);
  gsl_matrix_free(dB);
  gsl_matrix_free(cov);
  gsl_bspline_free(wksp);
#ifdef GSL1
  gsl_bspline_deriv_free(dwksp);
#endif
  gsl_integration_workspace_free(intwksp);
}

max_stars_repo_name: franciscaconcha/amuse-vader

max_issues_repo_name: franciscaconcha/amuse-vader

max_forks_repo_name: franciscaconcha/amuse-vader