text:
/******************************************************************************
 *
 * INTEL CONFIDENTIAL
 *
 * Copyright 2015 Intel Corporation All Rights Reserved.
 *
 * The source code contained or described herein and all documents related
 * to the source code (Material) are owned by Intel Corporation or its
 * suppliers or licensors. Title to the Material remains with
 * Intel Corporation or its suppliers and licensors. The Material contains
 * trade secrets and proprietary and confidential information of Intel or
 * its suppliers and licensors. The Material is protected by worldwide
 * copyright and trade secret laws and treaty provisions. No part of the
 * Material may be used, copied, reproduced, modified, published, uploaded,
 * posted, transmitted, distributed, or disclosed in any way without Intel's
 * prior express written permission.
 *
 * No license under any patent, copyright, trade secret or other intellectual
 * property right is granted to or conferred upon you by disclosure or
 * delivery of the Materials, either expressly, by implication, inducement,
 * estoppel or otherwise. Any license under such intellectual property rights
 * must be express and approved by Intel in writing.
 *
 *
 *   Workfile:   FwUpdateDebug.c
 *
 *   Abstract:   macros to wrap dbglog for simple FwUpdate debugging
 *
 ******************************************************************************/

#ifndef __FW_UPDATE_DEBUG_H__
#define __FW_UPDATE_DEBUG_H__

#include <gsl/span>
#include <iostream>

typedef enum
{
    PRINT_NONE = 0,
    PRINT_CRITICAL = 1,
    PRINT_ERROR,
    PRINT_WARNING,
    PRINT_INFO,
    PRINT_DEBUG,
    PRINT_DEBUG2,
    PRINT_ALL,
} dbg_level;

#define FWCRITICAL(MSG) PRINT(PRINT_CRITICAL, MSG)
#define FWERROR(MSG) PRINT(PRINT_ERROR, MSG)
#define FWWARN(MSG) PRINT(PRINT_WARNING, MSG)
#define FWINFO(MSG) PRINT(PRINT_INFO, MSG)
#define FWDEBUG(MSG) PRINT(PRINT_DEBUG, MSG)
#define FWDEBUG2(MSG) PRINT(PRINT_DEBUG2, MSG)

#define FWDUMP(D, L) DUMP(PRINT_DEBUG, D, L)

#define FW_UPDATE_DEBUG 1

#ifdef FW_UPDATE_DEBUG

extern dbg_level fw_update_get_dbg_level(void);
extern void fw_update_set_dbg_level(dbg_level l);

#define PRINT(LEVEL, MSG)                                                      \
    do                                                                         \
    {                                                                          \
        if ((LEVEL) <= fw_update_get_dbg_level())                              \
        {                                                                      \
            std::stringstream ss;                                              \
            ss << '<' << LEVEL << '>' << __FUNCTION__ << ":" << __LINE__       \
               << ": " << MSG;                                                 \
            std::cerr << ss.str() << std::endl;                                \
        }                                                                      \
    } while (0)

void _dump(dbg_level lvl, const char* fn, int lineno, const char* bname,
           const gsl::span<const uint8_t>& buf);

void _dump(dbg_level lvl, const char* fn, int lineno, const char* bname,
           const void* buf, size_t len);

#define DUMP(LEVEL, BUF, ...)                                                  \
    do                                                                         \
    {                                                                          \
        if ((LEVEL) <= fw_update_get_dbg_level())                              \
        {                                                                      \
            _dump(LEVEL, __FUNCTION__, __LINE__, #BUF, BUF, ##__VA_ARGS__);    \
        }                                                                      \
    } while (0)

#else /* !FW_UPDATE_DEBUG */

#define PRINT(...)
#define DUMP(...)

#endif /* FW_UPDATE_DEBUG */

#endif /* __FW_UPDATE_DEBUG_H__ */

max_stars_repo_name: Intel-BMC/mtd-util

max_issues_repo_name: Intel-BMC/mtd-util

max_forks_repo_name: Intel-BMC/mtd-util


text:
/* -*- C -*- */
/**
 * Author: Pierre Schnizer
 * Date : January 2003
 */
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <pygsl/general_helpers.h>
#include <pygsl/utils.h>
#include <pygsl/profile.h>

PyGSL_API_EXTERN int
PyGSL_set_error_string_for_callback(PyGSL_error_info * info)
{
     PyObject *name_o = NULL;
     PyObject *callback;
     const char * message = "";
     const char * error_description = "";
     const char * mmesg;
     char * name, msg[1024];
     char *formatstring = "For the callback %s evaluted  for function %s, an error occured : %s";

     FUNC_MESS_BEGIN();
     callback = info->callback;
     if(info->message){
	  message = info->message;
     }
     if(info->error_description){
	  error_description = info->error_description;
     }

     if (message == NULL){
	  mmesg = "Unknown";
     }else{
	  mmesg = message;
     }
     assert(callback != NULL);
     name_o = PyObject_GetAttrString(callback, "__name__");
     if (name_o == NULL){
	  name_o = PyObject_GetAttrString(callback, "func_name");
     }
     if (name_o == NULL){
	 PyErr_SetString(PyExc_AttributeError,
			 "While I was treating an errornous callback object,"
			 " I found that it had no attribute '__name__'!");
	 PyGSL_ERROR("Could not get the name of the callback!", GSL_EBADFUNC);
	 goto fail;
     }
     if(!PyString_Check(name_o)){
	  PyErr_SetString(PyExc_TypeError,
			  " For an errornous callback object,"
			  " the attribute '__name__' was not a Python string!");
	 PyGSL_ERROR("Nameobject of the callback was not a string!", GSL_EBADFUNC);
	 goto fail;
     }
     name = PyString_AsString(name_o);

     /* A non completly standard but safe function. */
     FUNC_MESS("\tmakeing string");
     snprintf(msg, 1024, formatstring, name, mmesg, error_description);
     if(DEBUG > 2){
	  fprintf(stderr, "ERROR: \t%s", msg);
     }
     PyGSL_ERROR(msg, GSL_EBADFUNC);

     FUNC_MESS_END();
     return GSL_EBADFUNC;
     /* Py_XDECREF(name_o);  ??? */
 fail:
     return GSL_EBADFUNC;
     /* Py_XDECREF(name_o);  ??? */
}

PyGSL_API_EXTERN int
PyGSL_pyfloat_to_double(PyObject *object, double *result, PyGSL_error_info *info)
{

     PyObject *object1;
     char *msg="The object returned to the GSL Function could not be converted to float";

     FUNC_MESS_BEGIN();
     object1 = PyNumber_Float(object);
     if(object1 == NULL){
	 *result = gsl_nan();
	 if(info){
	      info->error_description = msg;
	      return PyGSL_set_error_string_for_callback(info);
	 }
	 DEBUG_MESS(2, "Not from call back treatment, normal error. info = %p", info);
	 PyGSL_ERROR(msg, GSL_EBADFUNC);
     }

     *result   = PyFloat_AsDouble(object1);
     DEBUG_MESS(3, "found a double of %f\n", *result);
     Py_DECREF(object1);

     PyGSL_INCREASE_float_transform_counter();
     FUNC_MESS_END();
     return GSL_SUCCESS;
}

PyGSL_API_EXTERN int
PyGSL_pylong_to_uint(PyObject *object, unsigned int *result, PyGSL_error_info *info)
{
     int flag;
     unsigned long int tmp;
     flag =PyGSL_pylong_to_ulong(object, &tmp, info);
     *result = (unsigned int) tmp;
     return flag;
}

PyGSL_API_EXTERN int
PyGSL_pylong_to_ulong(PyObject *object, unsigned long *result, PyGSL_error_info *info)
{

     PyObject *object1;
     char *msg="The object returned to the GSL Function could not be converted to unsigned long";

     object1 = PyNumber_Long(object);
     if(object1 == NULL){
	 *result = 0;
	 if(info){
	      info->error_description = msg;
	      return PyGSL_set_error_string_for_callback(info);
	 }
	 PyGSL_ERROR(msg, GSL_EINVAL);
     }

     *result   = PyLong_AsUnsignedLong(object1);
     if(DEBUG>2){
	  fprintf(stderr, "\t\t%s found a double of %ld\n", __FUNCTION__, *result);
     }
     Py_DECREF(object1);

     PyGSL_INCREASE_float_transform_counter();

     return GSL_SUCCESS;
}

PyGSL_API_EXTERN int
PyGSL_pyint_to_int(PyObject *object, int *result, PyGSL_error_info *info)
{

     PyObject *object1;
     char *msg="The object returned to the GSL Function could not be converted to int";
     long tmp;

     FUNC_MESS_BEGIN();
     object1 = PyNumber_Int(object);
     if(object1 == NULL){
	 *result = INT_MIN;
	 if(info){
	      info->error_description = msg;
	      return PyGSL_set_error_string_for_callback(info);
	 }
	 DEBUG_MESS(2, "Not from call back treatment, normal error. info = %p", info);
	 PyGSL_ERROR(msg, GSL_EINVAL);
     }

     tmp = PyInt_AsLong(object1);
     if(tmp > INT_MAX)
	  PyGSL_ERROR("Number too big for int", GSL_EINVAL);
     else if (tmp < INT_MIN)
	  PyGSL_ERROR("Number too small for int", GSL_EINVAL);

     *result = (int) tmp;
     DEBUG_MESS(3, "found a int of %d\n", *result);
     Py_DECREF(object1);
     FUNC_MESS_END();
     return GSL_SUCCESS;
}

/*
 * Checks following conditions:
 *  For No Arguments: Got Py_None and No Error
 *  For 1  Argument:  Got an Object, No None  and No Error
 *         (Is None a legal return for one object? I think so.) On the other hand its a
 *         callback and Conversions are waiting, so its good not to accept None.
 * For 2  Arguments: Got a tuple of approbriate size
 */
#define PyGSL_CHECK_PYTHON_RETURN(object, nargs, info)                              \
  (                                                                                \
        (  ( (nargs) == 0 ) && ( object ) && ( Py_None == (object) ) && ( !PyErr_Occurred() ) )   \
    ||  (  ( (nargs) == 1 ) && ( object ) && ( Py_None != (object) ) && ( !PyErr_Occurred() ) )   \
    ||  (  ( (nargs) >  1 ) && ( object ) && ( PyTuple_Check((object)) ) &&                       \
                 ( (nargs) == PyTuple_GET_SIZE((object)) ) )                        \
 )                                                                                \
 ?                                                                                \
    GSL_SUCCESS                                                                     \
 :                                                                                \
   PyGSL_check_python_return((object), (nargs), (info))

PyGSL_API_EXTERN int
PyGSL_check_python_return(PyObject *object, int nargs, PyGSL_error_info  *info)
{
     int tuple_size, flag=-1;
     char *msg;

     FUNC_MESS_BEGIN();

     assert(info);

     if(object == NULL && PyErr_Occurred()){
	  /*
	   * Error was apparently raised by the function, so lets just add a
	   * traceback frame ....
	   */
	  info->error_description = "User function raised exception!";
	  PyGSL_add_traceback(NULL, "Unknown file", info->message, __LINE__);
	  return GSL_EBADFUNC;
     }
     if(PyErr_Occurred()){
	  info->error_description = "Function raised an exception.";
	  PyGSL_add_traceback(NULL, "Unknown file", info->message, __LINE__);
	  return GSL_EBADFUNC;
	  /* return PyGSL_set_error_string_for_callback(info); */
     }

     /* Expected No argumets */
     if(nargs == 0){
	  if(object != Py_None){
	       info->error_description = "I expected 0 arguments, but I got an object different from None.";
	       return PyGSL_set_error_string_for_callback(info);
	  } else {
	       return GSL_SUCCESS;
	  }
     }

     if(nargs == 1){
	  if(object == Py_None){
	       info->error_description = "Expected 1 argument, but None was returned. This value is not acceptable for"
		    " the following arithmetic calculations.";
	       return PyGSL_set_error_string_for_callback(info);
	  } else {
	       return GSL_SUCCESS;
	  }
     }

     if(nargs > 1){
	  msg = (char *) malloc(256 * sizeof(char));

	  if(object == Py_None){
	       snprintf(msg, 256, "I expected %d arguments, but the function returned None!", nargs);
	       info->error_description = msg;
	       flag = PyGSL_set_error_string_for_callback(info);

	  } else if(!PyTuple_Check(object)){
	       snprintf(msg, 256, "Expected %d arguments, but I didn't get a tuple! "
			"Did you just return one argument?.", nargs);
	       info->error_description = msg;
	       flag = PyGSL_set_error_string_for_callback(info);

	  } else {
	       tuple_size = PyTuple_GET_SIZE(object);
	       if(tuple_size != nargs){
		    snprintf(msg, 256, "I expected %d arguments, but the function returned %d arguments! ",
			     nargs, tuple_size);
	       info->error_description = msg;
	       flag = PyGSL_set_error_string_for_callback(info);

	       } else {
		    flag = GSL_SUCCESS;
	       }
	  }
	  free(msg);
     }
     FUNC_MESS_END();
     return flag;
}

PyGSL_API_EXTERN void
PyGSL_clear_name(char *name, int size)
{
     int j;
     for(j = 0; j<size; j++){
	  if(name[j] == '-')
	       name[j] = '_';
     }
}

ext: c

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#include <gsl/gsl_errno.h>
#include <gsl/block/gsl_block.h>

#define BASE_DOUBLE
#include <gsl/templates_on.h>
#include <gsl/block/block_source.c>
#include <gsl/templates_off.h>
#undef BASE_DOUBLE

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
/*
 *  Copyright 2010, Weill Medical College of Cornell University
 *  All rights reserved.
 *
 *  This software is distributed WITHOUT ANY WARRANTY
 *  under license "license.txt" included with distribution and
 *  at http://neurodatabase.org/src/license.
 */

#include <gsl/gsl_vector.h>
#include<gsl/gsl_matrix.h>
#include<gsl/gsl_permutation.h>
#include<gsl/gsl_linalg.h>
#include<gsl/gsl_sf.h>
#include "../shared/toolkit_c.h"

#define BUB_MAXWORDS_LIMIT 100000

struct bub_options{
	double lambda_N;
	int use_weighting;
	double mesh_min;
	double mesh_max;
	double mesh_inc;
	int opt_use_weighting;
	double opt_mesh_min;
	double opt_mesh_max;
	double opt_mesh_inc;
	int bag_threshold;
};

int bag1(int N,double m,double *a,double lambda_0,struct bub_options *bubopts);
int bub1(int N,double m,double *a,double lambda_0,int K,int compat,struct bub_options *bubopts);
void make_binom(int N,int meshsize,double *p,double **B_jN);
void make_weight(double *p,int meshsize,double m,int use_weighting,double *f,double *c);
double compute_error(double *a,int N,double m,int compat,struct bub_options *bubopts);

int entropy_bub(struct hist1d *in,struct options_entropy *opts,struct estimate *entropy)
{
	double *a;
	int *h;
	double H;
	struct bub_options *bubopts;
	double eps,epsm;
	int N,i,j,max_bin;
	double m;
	int status;

	N = (*in).P;

	if(opts->possible_words_flag==0)
	{
		if(opts->bub_possible_words_strategy==0)
			m = (double)in->C;
		else if(opts->bub_possible_words_strategy==1)
			m = (double)in->N;
		else if(opts->bub_possible_words_strategy==2)
			m = MIN(max_possible_words(in,1),BUB_MAXWORDS_LIMIT);
		else
			m = (double)in->C;
	}
	else
		if(opts->possible_words>0)
			m = opts->possible_words;
		else
			switch((int)(opts->possible_words))
			{
				case 0: /* Inf */
					entropy->value = NAN; /* this value based on piece-meal derivation when m=INFINITY */
					return EXIT_SUCCESS;
				case -1: /* recommended */
					m = (double)in->C;
					break;
				case -2: /* unique */
					m = (double)in->C;
					break;
				case -3: /* total */
					m = (double)in->P;
					break;
				case -4: /* possible */
					m = max_possible_words(in,0);
					break;
				case -5: /* min_tot_pos */
					m = max_possible_words(in,1);
					break;
				case -6: /* min_lim_tot_pos */
					m = MIN(BUB_MAXWORDS_LIMIT,max_possible_words(in,1));
					break;
			}

	/* First find a vector */
	a = (double *)malloc((N+1)*sizeof(double));

	bubopts = (struct bub_options *)malloc(sizeof(struct bub_options));
	(*bubopts).bag_threshold = 20;
	if(N<(*bubopts).bag_threshold)
	{
		(*bubopts).lambda_N = 0;
		(*bubopts).mesh_min = 0;
		(*bubopts).mesh_max = 1;
		(*bubopts).mesh_inc = 1/(5*(double)N);
		(*bubopts).use_weighting = 0;
		status = bag1(N,m,a,(*opts).bub_lambda_0,bubopts);
	}
	else
	{
		eps = (1/(double)N)*1e-10;
		epsm = (1/m)*1e-10;
		(*bubopts).lambda_N = 1;
		(*bubopts).mesh_min = eps;
		(*bubopts).mesh_max = MIN(1,3/(double)N)-eps;
		(*bubopts).mesh_inc = MIN(1,3/(double)N)/100;
		(*bubopts).lambda_N = 1;
		(*bubopts).use_weighting = 0;
		(*bubopts).opt_mesh_min = epsm;
		(*bubopts).opt_mesh_max = MIN(1,3/m)-epsm;
		(*bubopts).opt_mesh_inc = MIN(1,3/m)/100;
		(*bubopts).opt_use_weighting = 0;
		status = bub1(N,m,a,(*opts).bub_lambda_0,(*opts).bub_K,(*opts).bub_compat,bubopts);
	}

	/* Get counts of counts */
	h = (int *)calloc(N+1,sizeof(int));
	for(i=0; i<(*in).C; i++)
		h[(int)(*in).wordcnt[i]]++;

	/* Then compute entropy */
	H = 0;
	for(j=0; j<=N; j++)
		H += a[j]*h[j];
	entropy->value = NAT2BIT(H);

	free(a);
	free(h);
	free(bubopts);

	return EXIT_SUCCESS; /*	TODO return status; reports errors with staverify.m - need to look into these errors and possibly pass them to message structure */
}

int bag1(int N,double m,double *a,double lambda_0,struct bub_options *bubopts)
{
	int meshsize;
	int i,j;
	double *p,*Y,**X,**B_jN,*A_data,**D,*I0,*IN;
	double *f,c;
	gsl_vector *b,*x,*tau,*norm,*S,*work;
	gsl_matrix *A,*V;
	int *signum;
	int d_status, s_status, status = EXIT_SUCCESS;

	meshsize = floor(((*bubopts).mesh_max-(*bubopts).mesh_min)/(*bubopts).mesh_inc)+1;

	A = gsl_matrix_calloc(meshsize+(N+1)+2,N+1);
	A_data = (*A).data;
	X = VectorToMatrixDouble(&A_data[0],meshsize,N+1);
	D = VectorToMatrixDouble(&A_data[meshsize*(N+1)],N+1,N+1);
	I0 = &A_data[(meshsize+(N+1))*(N+1)];
	IN = &A_data[(meshsize+(N+1)+1)*(N+1)];

	b = gsl_vector_calloc(meshsize+(N+1)+2);
	Y = (*b).data;

	/* Make mesh */
	p = (double *)malloc(meshsize*sizeof(double));
	for(i=0; i<meshsize; i++)
		p[i] = (*bubopts).mesh_min + i*(*bubopts).mesh_inc;

	/* Make weighting */
	f = (double *)malloc(meshsize*sizeof(double));
	make_weight(p,meshsize,m,(*bubopts).use_weighting,f,&c);

	/* Make H */
	for(i=0; i<meshsize; i++)
		Y[i] = f[i]*XLOGX(p[i]);

	/* Make binomial coefficents */
	B_jN = MatrixDouble(meshsize,N+1);
	make_binom(N,meshsize,p,B_jN);

	for(i=0; i<meshsize; i++)
		for(j=0; j<N+1; j++)
			X[i][j] = f[i]*B_jN[i][j];

	/* Make D */
	for(j=0; j<N; j++)
	{
		D[j][j]=-1*sqrt((double)N)/c;
		D[j][j+1]=1*sqrt((double)N)/c;
	}

	/* Make I0 and IN */
	I0[0] = sqrt(((double)N)*lambda_0)/c;
	IN[N] = sqrt(((double)N)*(*bubopts).lambda_N)/c;

	/* Solve the least squares problem */
	V = gsl_matrix_calloc(N+1,N+1);
	S = gsl_vector_calloc(N+1);
	work = gsl_vector_calloc(N+1);
	d_status = gsl_linalg_SV_decomp(A,V,S,work);
	if(d_status)
		status = EXIT_FAILURE;
	x = gsl_vector_calloc(N+1);
	s_status = gsl_linalg_SV_solve(A,V,S,b,x);
	if(s_status)
		status = EXIT_FAILURE;

	/* Read out the answer from x */
	memcpy(a,(*x).data,(N+1)*sizeof(double));

	/* Free memory */
	gsl_matrix_free(V);
	gsl_vector_free(S);
	gsl_vector_free(work);
	gsl_vector_free(x);

	FreeMatrixDouble(B_jN);
	free(X);
	free(D);
	free(f);
	free(p);
	gsl_matrix_free(A);
	gsl_vector_free(b);

	return status;
}

int bub1(int N,double m,double *a,double lambda_0,int K,int compat,struct bub_options *bubopts)
{
	int meshsize;
	int i,j,k;
	double *p,*Y,**X,**B_jN,*A_data,*b_data,**D,*I0,*IN;
	double *f,c;
	gsl_vector *b,*x,*tau,*norm,*S,*work;
	gsl_matrix *A,*V;
	int d_status, s_status, status = EXIT_SUCCESS;
	double best_error,max_error,h_MM;
	double *a_MM,*a_temp;
	int k1;

	/* Compute the a coefficients using the Miller-Madow method */
	a_MM = (double *)calloc(N+1,sizeof(double));
	for(j=0; j<=N; j++)
		a_MM[j] = XLOGX((double)j/(double)N) + (1-((double)j/(double)N))/(2*(double)N);

	meshsize = floor(((*bubopts).mesh_max-(*bubopts).mesh_min)/(*bubopts).mesh_inc)+1;

	/* Make mesh */
	p = (double *)malloc(meshsize*sizeof(double));
	for(i=0; i<meshsize; i++)
		p[i] = (*bubopts).mesh_min + i*(*bubopts).mesh_inc;

	/* Make weighting */
	f = (double *)malloc(meshsize*sizeof(double));
	make_weight(p,meshsize,m,(*bubopts).use_weighting,f,&c);

	B_jN = MatrixDouble(meshsize,N+1);
	make_binom(N,meshsize,p,B_jN);

	best_error = HUGE_VAL;

	for(k=0; k<MIN(K,N); k++)
	{
		k1 = k+1;

		a_temp = (double *)calloc(N+1,sizeof(double));
		memcpy(a_temp,a_MM,(N+1)*sizeof(double));

		/*** MAKE A ***/
		A = gsl_matrix_calloc(meshsize+k1+2,k1);
		A_data = (*A).data;
		X = VectorToMatrixDouble(&A_data[0],meshsize,k1);
		D = VectorToMatrixDouble(&A_data[meshsize*k1],k1,k1);
		I0 = &A_data[(meshsize+k1)*k1];
		IN = &A_data[(meshsize+k1+1)*k1];

		/*** MAKE B ***/
		b = gsl_vector_calloc(meshsize+k1+2);
		b_data = (*b).data;
		Y = &b_data[0];
		if(compat)
			b_data[meshsize+k1+1] = (sqrt((double)N)/c)*a_MM[k];
		else
			b_data[meshsize+k1+1] = (sqrt((double)N)/c)*a_MM[k1];

		for(i=0; i<meshsize; i++)
		{
			h_MM = 0;
			for(j=k1; j<=N; j++)
				h_MM += B_jN[i][j]*a_MM[j];

			Y[i] = f[i]*(XLOGX(p[i]) - h_MM);
		}

		for(i=0; i<meshsize; i++)
			for(j=0; j<k1; j++)
				X[i][j] = f[i]*B_jN[i][j];

		/* Make D */
		for(j=0; j<k1-1; j++)
		{
			D[j][j]=-1*sqrt((double)N)/c;
			D[j][j+1]=1*sqrt((double)N)/c;
		}

		/* Make I0 and IN */
		I0[0] = sqrt(((double)N)*lambda_0)/c;
		IN[k1-1] = sqrt(((double)N)*(*bubopts).lambda_N)/c;

		/* Solve the least squares problem */
		V = gsl_matrix_calloc(k1,k1);
		S = gsl_vector_calloc(k1);
		work = gsl_vector_calloc(k1);
		d_status = gsl_linalg_SV_decomp(A,V,S,work);
		if(d_status)
			status = EXIT_FAILURE;
		x = gsl_vector_calloc(k1);
		s_status = gsl_linalg_SV_solve(A,V,S,b,x);
		if(s_status)
			status = EXIT_FAILURE;

		/* Read out the answer from x */
		memcpy(a_temp,(*x).data,k1*sizeof(double));

#ifdef DEBUG
		printf("k1=%d\n",k1);
		for(j=0; j<=k; j++)
			printf("%f ",a_temp[j]);
		printf("\n");
#endif

		/* Find the error */
		max_error = compute_error(a_temp,N,m,compat,bubopts);

		if(max_error<best_error)
		{
			best_error = max_error;
			memcpy(a,a_temp,(N+1)*sizeof(double));
		}

		free(a_temp);

		/* Free memory */
		gsl_matrix_free(V);
		gsl_vector_free(S);
		gsl_vector_free(work);
		gsl_vector_free(x);
		free(X);
		free(D);
		gsl_matrix_free(A);
		gsl_vector_free(b);
	}

	free(f);
	free(p);
	FreeMatrixDouble(B_jN);
	free(a_MM);
}

void make_binom(int N,int meshsize,double *p,double **B_jN)
{
	int i,j;
	double *NCj;

	NCj = (double *)calloc(N+1,sizeof(double));

	for(j=0; j<N+1; j++)
		NCj[j] = gsl_sf_lngamma(N+1)-gsl_sf_lngamma(j+1)-gsl_sf_lngamma(N-j+1);

	for(i=0; i<meshsize; i++)
		if(p[i]==0)
		{
			B_jN[i][0]=1;
			for(j=1; j<N+1; j++)
				B_jN[i][j]=0;
		}
		else if(p[i]==1)
		{
			for(j=0; j<N; j++)
				B_jN[i][j]=0;
			B_jN[i][N]=1;
		}
		else
		{
			for(j=0; j<N+1; j++)
				B_jN[i][j] = exp(NCj[j] + log(p[i])*j + log(1-p[i])*(N-j));
		}

	free(NCj);
}

void make_weight(double *p,int meshsize,double m,int use_weighting,double *f,double *c)
{
	int i;

	for(i=0; i<meshsize; i++)
		if(use_weighting && p[i]>1/m)
			f[i] = 1./p[i];
		else
			f[i] = m;

	if(use_weighting)
		*c = 2;
	else
		*c = 1;
}

double compute_error(double *a,int N,double m,int compat,struct bub_options *bubopts)
{
	int i,j;
	double max_bias,max_error,max_variance0,max_variance1;
	double *p,*f,**B_jN;
	double temp1,temp2,temp3,temp4,temp5;
	int meshsize;
	double c;

	meshsize = floor(((*bubopts).opt_mesh_max-(*bubopts).opt_mesh_min)/(*bubopts).opt_mesh_inc)+1;

	/* Make mesh */
	p = (double *)malloc(meshsize*sizeof(double));
	for(i=0; i<meshsize; i++)
		p[i] = (*bubopts).opt_mesh_min + i*(*bubopts).opt_mesh_inc;

	/* Make weighting */
	f = (double *)malloc(meshsize*sizeof(double));
	make_weight(p,meshsize,m,(*bubopts).opt_use_weighting,f,&c);

	/* Make B_jN */
	B_jN = MatrixDouble(meshsize,N+1);
	make_binom(N,meshsize,p,B_jN);

	/* Compute the worst-case bias */
	max_bias = 0;
	for(i=0; i<meshsize; i++)
	{
		temp1 = 0;
		for(j=0; j<N+1; j++)
			temp1+=B_jN[i][j]*a[j];
		temp2 = f[i]*fabs(temp1 - XLOGX(p[i]));
		if(temp2>max_bias)
			max_bias = temp2;
	}

	/* Steele loose bound on variance */
	max_variance1 = 0;
	for(i=0; i<meshsize; i++)
	{
		temp3 = 0;
		for(j=1; j<N+1; j++)
			temp3 += B_jN[i][j]*((double)j/(double)N)*pow(a[j]-a[j-1],2);
		temp4 = f[i]*temp3;
		if(temp4>max_variance1)
			max_variance1 = temp4;
	}
	if(compat)
		max_variance1 *= 4*N*c;
	else
		max_variance1 *= 2*N*c;

	/* McDiarmid bound on variance */
	max_variance0 = 0;
	for(j=1; j<N+1; j++)
	{
		temp5 = fabs(a[j]-a[j-1]);
		if(temp5>max_variance0)
			max_variance0=temp5;
	}
	max_variance0 = N*pow(max_variance0,2);

#ifdef DEBUG
	printf("maxbias=%f maxvariance0=%f maxvariance1=%f\n",max_bias,max_variance0,max_variance1);
#endif

	max_error = NAT2BIT(sqrt(pow(max_bias,2) + MIN(max_variance0,max_variance1)));
	free(p);
	free(f);
	FreeMatrixDouble(B_jN);

	return max_error;
}

max_stars_repo_name: joaodornas/random-code

max_issues_repo_name: joaodornas/random-code

max_forks_repo_name: joaodornas/random-code


text:
/**
 * This file is part of yasimSBML (http://www.labri.fr/perso/ghozlane/metaboflux/about/yasimSBML.php)
 * Copyright (C) 2010 Amine Ghozlane from LaBRI and University of Bordeaux 1
 *
 * yasimSBML is free software: you can redistribute it and/or modify
 * it under the terms of the Lesser GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * yasimSBML 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 Lesser GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/**
 * \file yasimSBML.c
 * \brief Main program
 * \author {Amine Ghozlane}
 * \version 1.0
 * \date 27 octobre 2009
 */
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <assert.h>
#include <unistd.h>
#include <gsl/gsl_rng.h>
#include <sbml/SBMLTypes.h>
#include "especes.h"
#include "simulation.h"

/**
 * \fn int main(int argc, char **argv)
 * \author Amine Ghozlane
 * \brief  Enter in the program for simulated annealing
 * \param  argc Number of arguments
 * \param  argv List of arguments
 * \return EXIT_SUCCESS Normal stop of the program
 */
int main(int argc, char **argv) {
  SBMLDocument_t *file=NULL;
  Model_t *mod=NULL;
  unsigned int errors = 0;
  /*unsigned int level;
  unsigned int version;*/

  /* File adress */
  if (argc <= 2) {
      fprintf(stderr, "Usage ./yasimSBML.exe sbml_file.xml number_of_simulation\n");
      return 1;
  }

  /* Load File  */
  file = readSBML(argv[1]);
  errors = SBMLDocument_getNumErrors(file);

  /* Error File */
  if (errors > 0) {
      printf("Encountered the following SBML error(s):\n");
      SBMLDocument_printErrors(file, stdout);
      printf("Simulation skipped.  Please correct the problems above first.\n");
      return errors;
  }

  /*level = SBMLDocument_getLevel(file);
  version = SBMLDocument_getVersion(file);
  printf("SBML Level : %d, version : %d\n", level, version);*/

  /* Testing mod file */
  mod = SBMLDocument_getModel(file);
  if (mod == NULL)
    fprintf(stderr, "Error: file %s doesn't exist!\n", argv[1]);
  /*compute_simulation(mod);*/
  SBML_compute_simulation_mean(mod,atoi(argv[2]));

  /* free memory */
  Model_free(mod);
  /*SBMLDocument_free(file);*/

  return EXIT_SUCCESS;
}

max_stars_repo_name: aghozlane/yasimSBML

max_issues_repo_name: aghozlane/yasimSBML

max_forks_repo_name: aghozlane/yasimSBML


text:
/*
 This file is part of the Princeton FCMA Toolbox
 Copyright (c) 2013 the authors (see AUTHORS file)
 For license terms, please see the LICENSE file.
*/

#ifndef COMMON_H
#define COMMON_H

#include <vector>
#include <map>
#include <list>
#include <algorithm>
#include <sys/types.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <fcntl.h>
#include <unistd.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <climits>
#include <iomanip>
#include <stdlib.h>
// define NDEBUG here to disable assertions
// if/when "release" builds are supported
// #define NDEBUG
#include <cassert>

// nifti1_io.h is from nifticlib http://nifti.nimh.nih.gov
// see README in fcma-toolbox.git/deps/
#include <nifti1_io.h>

// SSE/AVX instrinsics
#if defined(_MSC_VER)
// Microsoft C/C++-compatible compiler */
#include <intrin.h>
#elif defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
// GCC-compatible compiler (gcc,clang,icc) targeting x86/x86-64
#include <x86intrin.h>
#endif

// OpenMP : gcc >= 4.7 or clang-omp http://clang-omp.github.io
#include <omp.h>

// Aligned memory blocks
#ifdef __INTEL_COMPILER
#include <offload.h>
#else
#include <mm_malloc.h>
#endif
#if defined(__GNUC__)
// gcc supports this syntax
#define ALIGNED(x) __attribute__((aligned(x)))
#else
// visual c++, clang, icc
#define ALIGNED(x) __declspec(align(x))
#endif

// BLAS dependency
#ifdef USE_MKL
#include <mkl.h>
#else
typedef int MKL_INT;
#if defined __APPLE__
#include <Accelerate/Accelerate.h>
#else
extern "C" {
#include <cblas.h>
}
#endif
#endif

// Matrix multiplication parameters
#define TINYNUM 1e-4
#define LOGISTICTRHEHOLD 1e-6
#define MAXFILENAMELENGTH 300
#define MAXSUBJS 100
#define MAXTRIALPERSUBJ 100
#define MAXSVMITERATION 20000000

// MPI communication tags
#define COMPUTATIONTAG 1
#define LENGTHTAG 2
#define VOXELCLASSIFIERTAG 3
#define ELAPSETAG 4
#define POSITIONTAG 5
#define SECONDORDERTAG 6

enum Task {
  Corr_Based_SVM = 0,
  Corr_Based_Dis,
  Acti_Based_SVM,
  Corr_Sum,
  Corr_Mask_Classification,
  Corr_Mask_Cross_Validation,
  Acti_Mask_Classification,
  Acti_Mask_Cross_Validation,
  Corr_Visualization,
  Marginal_Screening,
  Error_Type = -1
};

typedef unsigned long long uint64;
typedef unsigned short uint16;

typedef struct raw_matrix_t {
  std::string sname;  // subject name (file name without extension)
  int sid;            // subject id
  int row;
  int col;
  int nx, ny, nz;
  float* matrix;
} RawMatrix;

typedef struct trial_data_t {
  int nTrials;
  int nVoxels;
  int nCols;
  int* trialLengths;  // keep all trials data, subject by subject
  int* scs;
  float* data;  // normalized data for correlation
  trial_data_t(int x, int y) : nTrials(x), nVoxels(y) {}
} TrialData;

typedef struct corr_matrix_t {
  int sid;      // subject id
  int tlabel;   // trial label
  int sr;       // starting row id
  int step;     // row of this matrix
  int nVoxels;  // col of this matrix
  float* matrix;
} CorrMatrix;

// Trial: data structure for the start and end point of a trial
typedef struct trial_t {
  int tid;
  int sid;
  int label;
  int sc, ec;
  // tid_withinsubj: block id within each subject,
  // to be used in averaged matrix of searchlight
  int tid_withinsubj;
} Trial;

typedef struct voxel_t {
  int* vid;  // contains global voxel ids
  float* corr_vecs;
  float* kernel_matrices;  // contains precomputed kernel matrix
  int nTrials;  // row
  int nVoxels;  // col
  // voxel_t(int x, int y, int z) : vid(x), nTrials(y), nVoxels(z) {}
} Voxel;

typedef struct voxel_Score_t {
  int vid;
  float score;
} VoxelScore;

// VoxelXYZ: voxel's 3-d coordinates in mm
typedef struct voxelxyz_t {
  int x, y, z;
} VoxelXYZ;

typedef ALIGNED(64) struct WIDELOCK_T {
  omp_lock_t lock;
} widelock_t;

extern unsigned long long total_count;

#endif

ext: h

max_stars_repo_name: PrincetonUniversity/fcma-toolbox

max_issues_repo_name: PrincetonUniversity/fcma-toolbox

max_forks_repo_name: PrincetonUniversity/fcma-toolbox


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

/*
 * This module implements the sequential TSQR algorithm
 * described in
 *
 * [1] Demmel, J., Grigori, L., Hoemmen, M. F., and Langou, J.
 *     "Communication-optimal parallel and sequential QR and LU factorizations",
 *     UCB Technical Report No. UCB/EECS-2008-89, 2008.
 *
 * The algorithm operates on a tall least squares system:
 *
 * [ A_1 ] x = [ b_1 ]
 * [ A_2 ]     [ b_2 ]
 * [ ... ]     [ ... ]
 * [ A_k ]     [ b_k ]
 *
 * as follows:
 *
 * 1. Initialize
 *    a. [Q_1,R_1] = qr(A_1)
 *    b. z_1 = Q_1^T b_1
 * 2. Loop i = 2:k
 *    a. [Q_i,R_i] = qr( [ R_{i-1} ; A_i ] )
 *    b. z_i = Q_i^T [ z_{i-1} ; b_i ]
 * 3. Output:
 *    a. R = R_k
 *    b. Q^T b = z_k
 *
 * Step 2(a) is optimized to take advantage
 * of the sparse structure of the matrix
 */

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

typedef struct
{
  size_t p;             /* number of columns of LS matrix */
  int init;             /* QR system has been initialized */
  int svd;              /* SVD of R has been computed */
  double normb;         /* || b || for computing residual norm */

  gsl_vector *tau;      /* Householder scalars, p-by-1 */
  gsl_matrix *R;        /* [ R ; A_i ], size p-by-p */
  gsl_vector *QTb;      /* [ Q^T b ; b_i ], size p-by-1 */

  gsl_multifit_linear_workspace *multifit_workspace_p;
} tsqr_state_t;

static void *tsqr_alloc(const size_t p);
static void tsqr_free(void *vstate);
static int tsqr_reset(void *vstate);
static int tsqr_accumulate(gsl_matrix * A, gsl_vector * b,
                           void * vstate);
static int tsqr_solve(const double lambda, gsl_vector * x,
                      double * rnorm, double * snorm,
                      void * vstate);
static int tsqr_rcond(double * rcond, void * vstate);
static int tsqr_lcurve(gsl_vector * reg_param, gsl_vector * rho,
                       gsl_vector * eta, void * vstate);
static int tsqr_svd(tsqr_state_t * state);
static double tsqr_householder_transform (double *v0, gsl_vector * v);
static int tsqr_householder_hv (const double tau, const gsl_vector * v, double *w0,
                                gsl_vector * w);
static int tsqr_householder_hm (const double tau, const gsl_vector * v, gsl_matrix * R,
                                gsl_matrix * A);
static int tsqr_QR_decomp (gsl_matrix * R, gsl_matrix * A, gsl_vector * tau);

/*
tsqr_alloc()
  Allocate workspace for solving large linear least squares
problems using the TSQR approach

Inputs: p    - number of columns of LS matrix

Return: pointer to workspace
*/

static void *
tsqr_alloc(const size_t p)
{
  tsqr_state_t *state;

  if (p == 0)
    {
      GSL_ERROR_NULL("p must be a positive integer",
                     GSL_EINVAL);
    }

  state = calloc(1, sizeof(tsqr_state_t));
  if (!state)
    {
      GSL_ERROR_NULL("failed to allocate tsqr state", GSL_ENOMEM);
    }

  state->p = p;
  state->init = 0;
  state->svd = 0;
  state->normb = 0.0;

  state->R = gsl_matrix_alloc(p, p);
  if (state->R == NULL)
    {
      tsqr_free(state);
      GSL_ERROR_NULL("failed to allocate R matrix", GSL_ENOMEM);
    }

  state->QTb = gsl_vector_alloc(p);
  if (state->QTb == NULL)
    {
      tsqr_free(state);
      GSL_ERROR_NULL("failed to allocate QTb vector", GSL_ENOMEM);
    }

  state->tau = gsl_vector_alloc(p);
  if (state->tau == NULL)
    {
      tsqr_free(state);
      GSL_ERROR_NULL("failed to allocate tau vector", GSL_ENOMEM);
    }

  state->multifit_workspace_p = gsl_multifit_linear_alloc(p, p);
  if (state->multifit_workspace_p == NULL)
    {
      tsqr_free(state);
      GSL_ERROR_NULL("failed to allocate multifit workspace", GSL_ENOMEM);
    }

  return state;
}

static void
tsqr_free(void *vstate)
{
  tsqr_state_t *state = (tsqr_state_t *) vstate;

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

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

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

  if (state->multifit_workspace_p)
    gsl_multifit_linear_free(state->multifit_workspace_p);

  free(state);
}

static int
tsqr_reset(void *vstate)
{
  tsqr_state_t *state = (tsqr_state_t *) vstate;

  gsl_matrix_set_zero(state->R);
  gsl_vector_set_zero(state->QTb);
  state->init = 0;
  state->svd = 0;
  state->normb = 0.0;

  return GSL_SUCCESS;
}

/*
tsqr_accumulate()
  Add a new block of rows to the QR system

Inputs: A      - new block of rows, n-by-p
        b      - new rhs vector n-by-1
        vstate - workspace

Return: success/error

Notes:
1) On output, the upper triangular portion of state->R(1:p,1:p)
contains current R matrix

2) state->QTb(1:p) contains current Q^T b vector

3) A and b are destroyed
*/

static int
tsqr_accumulate(gsl_matrix * A, gsl_vector * b, void * vstate)
{
  tsqr_state_t *state = (tsqr_state_t *) vstate;
  const size_t n = A->size1;
  const size_t p = A->size2;

  if (p != state->p)
    {
      GSL_ERROR("columns of A do not match workspace", GSL_EBADLEN);
    }
  else if (n != b->size)
    {
      GSL_ERROR("A and b have different numbers of rows", GSL_EBADLEN);
    }
  else if (state->init == 0)
    {
      int status;
      const size_t npmin = GSL_MIN(n, p);
      gsl_vector_view tau = gsl_vector_subvector(state->tau, 0, npmin);
      gsl_matrix_view R = gsl_matrix_submatrix(state->R, 0, 0, npmin, p);
      gsl_matrix_view Av = gsl_matrix_submatrix(A, 0, 0, npmin, p);
      gsl_vector_view QTb = gsl_vector_subvector(state->QTb, 0, npmin);
      gsl_vector_view bv = gsl_vector_subvector(b, 0, npmin);

      /* this is the first matrix block A_1, compute its (dense) QR decomposition */

      /* compute QR decomposition of A */
      status = gsl_linalg_QR_decomp(A, &tau.vector);
      if (status)
        return status;

      /* store upper triangular R factor in state->R */
      gsl_matrix_tricpy('U', 1, &R.matrix, &Av.matrix);

      /* compute ||b|| */
      state->normb = gsl_blas_dnrm2(b);

      /* compute Q^T b and keep the first p elements */
      gsl_linalg_QR_QTvec(A, &tau.vector, b);
      gsl_vector_memcpy(&QTb.vector, &bv.vector);

      state->init = 1;

      return GSL_SUCCESS;
    }
  else
    {
      int status;

      /* compute QR decomposition of [ R_{i-1} ; A_i ], accounting for
       * sparse structure */
      status = tsqr_QR_decomp(state->R, A, state->tau);
      if (status)
        return status;

      /* update ||b|| */
      state->normb = gsl_hypot(state->normb, gsl_blas_dnrm2(b));

      /*
       * compute Q^T [ QTb_{i - 1}; b_i ], accounting for the sparse
       * structure of the Householder reflectors
       */
      {
        size_t i;

        for (i = 0; i < p; i++)
          {
            const double ti = gsl_vector_get (state->tau, i);
            gsl_vector_const_view h = gsl_matrix_const_column (A, i);
            double *wi = gsl_vector_ptr(state->QTb, i);
            tsqr_householder_hv (ti, &(h.vector), wi, b);
          }
      }

      return GSL_SUCCESS;
    }
}

/*
tsqr_solve()
  Solve the least squares system:

chi^2 = || QTb - R x ||^2 + lambda^2 || x ||^2

using the SVD of R

Inputs: lambda - regularization parameter
        x      - (output) solution vector p-by-1
        rnorm  - (output) residual norm ||b - A x||
        snorm  - (output) solution norm ||x||
        vstate - workspace

Return: success/error
*/

static int
tsqr_solve(const double lambda, gsl_vector * x,
           double * rnorm, double * snorm,
           void * vstate)
{
  tsqr_state_t *state = (tsqr_state_t *) vstate;
  const size_t p = x->size;

  if (p != state->p)
    {
      GSL_ERROR("solution vector does not match workspace", GSL_EBADLEN);
    }
  else
    {
      int status;

      /* compute SVD of R if not already computed */
      if (state->svd == 0)
        {
          status = tsqr_svd(state);
          if (status)
            return status;
        }

      status = gsl_multifit_linear_solve(lambda, state->R, state->QTb, x, rnorm, snorm,
                                         state->multifit_workspace_p);
      if (status)
        return status;

      /*
       * Since we're solving a reduced square system above, we need
       * to account for the full residual vector:
       *
       * rnorm = || [ Q1^T b - R x ; Q2^T b ] ||
       *
       * where Q1 is the thin Q factor of X, and Q2
       * are the remaining columns of Q. But:
       *
       * || Q2^T b ||^2 = ||b||^2 - ||Q1^T b||^2
       *
       * so add this into the rnorm calculation
       */
      {
        double norm_Q1Tb = gsl_blas_dnrm2(state->QTb);
        double ratio = norm_Q1Tb / state->normb;
        double diff = 1.0 - ratio*ratio;

        if (diff > GSL_DBL_EPSILON)
          {
            double norm_Q2Tb = state->normb * sqrt(diff);
            *rnorm = gsl_hypot(*rnorm, norm_Q2Tb);
          }
      }

      return GSL_SUCCESS;
    }
}

/*
tsqr_lcurve()
  Compute L-curve of least squares system

Inputs: reg_param - (output) vector of regularization parameters
        rho       - (output) vector of residual norms
        eta       - (output) vector of solution norms
        vstate    - workspace

Return: success/error
*/

static int
tsqr_lcurve(gsl_vector * reg_param, gsl_vector * rho,
            gsl_vector * eta, void * vstate)
{
  tsqr_state_t *state = (tsqr_state_t *) vstate;
  int status;

  /* compute SVD of R if not already computed */
  if (state->svd == 0)
    {
      status = tsqr_svd(state);
      if (status)
        return status;
    }

  status = gsl_multifit_linear_lcurve(state->QTb, reg_param, rho, eta,
                                      state->multifit_workspace_p);

  /* now add contribution to rnorm from Q2 factor */
  {
    double norm_Q1Tb = gsl_blas_dnrm2(state->QTb);
    double ratio = norm_Q1Tb / state->normb;
    double diff = 1.0 - ratio*ratio;
    size_t i;

    if (diff > GSL_DBL_EPSILON)
      {
        double norm_Q2Tb = state->normb * sqrt(diff);

        for (i = 0; i < rho->size; ++i)
          {
            double *rhoi = gsl_vector_ptr(rho, i);
            *rhoi = gsl_hypot(*rhoi, norm_Q2Tb);
          }
      }
  }

  return status;
}

static int
tsqr_rcond(double * rcond, void * vstate)
{
  tsqr_state_t *state = (tsqr_state_t *) vstate;

  /* compute SVD of R if not already computed */
  if (state->svd == 0)
    {
      int status = tsqr_svd(state);
      if (status)
        return status;
    }

  *rcond = gsl_multifit_linear_rcond(state->multifit_workspace_p);

  return GSL_SUCCESS;
}

/*
tsqr_svd()
  Compute the SVD of the upper triangular
R factor. This allows us to compute the upper/lower
bounds on the regularization parameter and compute
the matrix reciprocal condition number.

Inputs: state - workspace

Return: success/error
*/

static int
tsqr_svd(tsqr_state_t * state)
{
  int status;

  status = gsl_multifit_linear_svd(state->R, state->multifit_workspace_p);
  if (status)
    {
      GSL_ERROR("error computing SVD of R", status);
    }

  state->svd = 1;

  return GSL_SUCCESS;
}

/*
tsqr_householder_transform()
  This routine is an optimized version of
gsl_linalg_householder_transform(), designed for the QR
decomposition of M-by-N matrices of the form:

T = [ R ]
    [ A ]

where R is N-by-N upper triangular, and A is (M-N)-by-N dense.
This routine computes a householder transformation (tau,v) of a
x so that P x = [ I - tau*v*v' ] x annihilates x(1:n-1). x will
be a subcolumn of the matrix T, and so its structure will be:

x = [ x0 ] <- 1 nonzero value for the diagonal element of R
    [ 0  ] <- N - j - 1 zeros, where j is column of matrix in [0,N-1]
    [ x  ] <- M-N nonzero values for the dense part A

Inputs: v0 - pointer to diagonal element of R
             on input, v0 = x0;
        v  - on input, x vector
             on output, householder vector v
*/

static double
tsqr_householder_transform (double *v0, gsl_vector * v)
{
  /* replace v[0:M-1] with a householder vector (v[0:M-1]) and
     coefficient tau that annihilate v[1:M-1] */

  double alpha, beta, tau ;

  /* compute xnorm = || [ 0 ; v ] ||, ignoring zero part of vector */
  double xnorm = gsl_blas_dnrm2(v);

  if (xnorm == 0)
    {
      return 0.0; /* tau = 0 */
    }

  alpha = *v0;
  beta = - (alpha >= 0.0 ? +1.0 : -1.0) * hypot(alpha, xnorm) ;
  tau = (beta - alpha) / beta ;

  {
    double s = (alpha - beta);

    if (fabs(s) > GSL_DBL_MIN)
      {
        gsl_blas_dscal (1.0 / s, v);
        *v0 = beta;
      }
    else
      {
        gsl_blas_dscal (GSL_DBL_EPSILON / s, v);
        gsl_blas_dscal (1.0 / GSL_DBL_EPSILON, v);
        *v0 = beta;
      }
  }

  return tau;
}

/*
tsqr_householder_hv()
  Apply Householder reflector to a vector. The Householder
reflectors are for the QR decomposition of the matrix

  [ R ]
  [ A ]

where R is p-by-p upper triangular and A is n-by-p dense.
Therefore all relevant components of the Householder
vector are stored in the columns of A, while the components
in R are 0, except for diag(R) which are 1.

The vector w to be transformed is partitioned as

  [ w1 ]
  [ w2 ]

where w1 is p-by-1 and w2 is n-by-1. The w2 portion
of w is transformed by v, but most of w1 remains unchanged
except for the first element, w0

Inputs: tau - Householder scalar
        v   - Householder vector, n-by-1
        w0  - (input/output)
              on input, w1(0);
              on output, transformed w1(0)
        w   - (input/output) n-by-1
              on input, vector w2;
              on output, P*w2
*/

static int
tsqr_householder_hv (const double tau, const gsl_vector * v, double *w0, gsl_vector * w)
{
  /* applies a householder transformation v to vector w */

  if (tau == 0)
    return GSL_SUCCESS ;

  {
    double d1, d;

    /* compute d1 = v(2:n)' w(2:n) */
    gsl_blas_ddot (v, w, &d1);

    /* compute d = v'w = w(1) + d1 since v(1) = 1 */
    d = *w0 + d1;

    /* compute w = w - tau (v) (v'w) */
    *w0 -= tau * d;
    gsl_blas_daxpy (-tau * d, v, w);
  }

  return GSL_SUCCESS;
}

/*
tsqr_householder_hm()
  Apply Householder reflector to a submatrix of

  [ R ]
  [ A ]

where R is p-by-p upper triangular and A is n-by-p dense.
The diagonal terms of R are already transformed by
tsqr_householder_transform(), so we just need to operate
on the submatrix A(:,i:p) as well as the superdiagonal
elements of R

Inputs: tau - Householder scalar
        v   - Householder vector
        R   - upper triangular submatrix of R, (p-i)-by-(p-i-1)
        A   - dense submatrix of A, n-by-(p-i)
*/

static int
tsqr_householder_hm (const double tau, const gsl_vector * v, gsl_matrix * R,
                     gsl_matrix * A)
{
  /* applies a householder transformation v,tau to matrix [ R ; A ] */

  if (tau == 0.0)
    {
      return GSL_SUCCESS;
    }
  else
    {
      size_t j;

      for (j = 0; j < A->size2; j++)
        {
          double R0j = gsl_matrix_get (R, 0, j);
          double wj;
          gsl_vector_view A1j = gsl_matrix_column(A, j);

          gsl_blas_ddot (&A1j.vector, v, &wj);
          wj += R0j;

          gsl_matrix_set (R, 0, j, R0j - tau *  wj);

          gsl_blas_daxpy (-tau * wj, v, &A1j.vector);
        }

      return GSL_SUCCESS;
    }
}

/*
tsqr_QR_decomp()
  Compute the QR decomposition of the matrix

  [ R ]
  [ A ]

where R is p-by-p upper triangular and A is n-by-p dense.

Inputs: R   - upper triangular p-by-p matrix
        A   - dense n-by-p matrix
        tau - Householder scalars
*/

static int
tsqr_QR_decomp (gsl_matrix * R, gsl_matrix * A, gsl_vector * tau)
{
  const size_t n = A->size1;
  const size_t p = R->size2;

  if (R->size2 != A->size2)
    {
      GSL_ERROR ("R and A have different number of columns", GSL_EBADLEN);
    }
  else if (tau->size != p)
    {
      GSL_ERROR ("size of tau must be p", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      for (i = 0; i < p; i++)
        {
          /* Compute the Householder transformation to reduce the j-th
             column of the matrix [ R ; A ] to a multiple of the j-th unit vector,
             taking into account the sparse structure of R */

          gsl_vector_view c = gsl_matrix_column(A, i);
          double *Rii = gsl_matrix_ptr(R, i, i);
          double tau_i = tsqr_householder_transform(Rii, &c.vector);

          gsl_vector_set (tau, i, tau_i);

          /* Apply the transformation to the remaining columns and
             update the norms */

          if (i + 1 < p)
            {
              gsl_matrix_view Rv = gsl_matrix_submatrix(R, i, i + 1, p - i, p - (i + 1));
              gsl_matrix_view Av = gsl_matrix_submatrix(A, 0, i + 1, n, p - (i + 1));
              tsqr_householder_hm (tau_i, &(c.vector), &(Rv.matrix), &(Av.matrix));
            }
        }

      return GSL_SUCCESS;
    }
}

static const gsl_multilarge_linear_type tsqr_type =
{
  "tsqr",
  tsqr_alloc,
  tsqr_reset,
  tsqr_accumulate,
  tsqr_solve,
  tsqr_rcond,
  tsqr_lcurve,
  tsqr_free
};

const gsl_multilarge_linear_type * gsl_multilarge_linear_tsqr =
  &tsqr_type;

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#define NOMINMAX

#include <iostream>
#include <vector>
#include <string>
#include <string_view>
#include <sstream>
#include <variant>
#include <codecvt>
#include <locale>
#include <filesystem>
#include <fstream>
#include <chrono>
#include <functional>
#include <numeric>

#include <wrl/client.h>
#include <wil/result.h>
#include <Windows.h>
#include <d3d12.h>
#include <dxcore.h>

#include <gsl/gsl>

#define DML_TARGET_VERSION_USE_LATEST
#include <DirectML.h>
#include <directx/d3dx12.h>
#include "DirectMLX.h"

#include <rapidjson/document.h>
#include <rapidjson/istreamwrapper.h>

#include <fmt/format.h>

max_stars_repo_name: jstoecker/dxdispatch

max_issues_repo_name: jstoecker/dxdispatch

max_forks_repo_name: jstoecker/dxdispatch


text:
/* Copyright 2017 Peter Williams and collaborators
 * Licensed under the MIT License.
 */

#include <Python.h>

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

#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_poly.h>

#define COUNT(arr) (sizeof(arr) / sizeof(arr[0]))

typedef enum result_t {
    RESULT_OK = 0,
    RESULT_LAMBDA_FAILED = 1,
    RESULT_X_FAILED = 2,
    RESULT_GSL_ERROR = 3,
    RESULT_HIT_CRITICAL_POINT = 4,
} result_t;

typedef enum handedness_t {
    R_MODE = 1,
    L_MODE = -1,
} handedness_t;

typedef enum wave_filtering_t {
    WF_ALL_WAVES = 0,
    WF_NO_FORWARD_MODES = 1 << 0,
    WF_NO_BACKWARD_MODES = 1 << 1,
} wave_filtering_t;

typedef enum integration_type_t {
    ITYPE_BOUNCE_AVERAGED,
    ITYPE_LOCAL,
} integration_type_t;

typedef struct parameters_t {
    double E; /* dimensionless kinetic energy */
    double sin_alpha; /* sine of the pitch angle */
    double Omega_e; /* local electron gyrofrequency in rad/s */
    handedness_t handedness; /* handedness of wave; `s` in Summers/Shprits notation */
    double alpha_star; /* (cyclotron freq / plasma freq)**2 */
    double R; /* magnetic wave energy perturbation: (delta B / B)**2 */
    double x_m; /* center of wave frequency spectrum in units of cyclotron freq */
    double delta_x; /* width of wave frequency spectrum in units of cyclotron freq */
    double max_wave_latitude; /* maximum latitude at which waves are found, in radians */
    wave_filtering_t wave_filtering; /* which kinds of waves to filter out */
} parameters_t;

typedef struct state_t {
    parameters_t p;
    integration_type_t itype;
    result_t last_error;

    double gamma; /* Lorentz factor */
    double a; /* See Summers 2005, just after Eqn 24 */
    double b; /* See Summers 2005, just after Eqn 25 */
    double beta; /* velocity over speed of light */

    double latitude; /* current magnetic latitude, radians */
    double sin_alpha; /* current sin(alpha); varies since mu is conserved */
    double mu; /* current cos(alpha) */
    double Omega_e; /* current Omega_e; varies with latitude since it scales as B */

    int n_xy; /* Number of resonances encountered */
    double x[3]; /* x values of resonances */
    double y[3]; /* y values of resonances */
    double Q[3]; /* (1 - x cos alpha / y beta) */
    double R[3]; /* complex expression in summand */
} state_t;

typedef struct coefficients_t {
    double dimensionless_p; /* momentum over (m_0 * c) */
    double Daa;
    double err_Daa;
    double Dap_on_p;
    double err_Dpp_on_p2;
    double Dpp_on_p2;
    double err_Dap_on_p;
} coefficients_t;

static state_t *global_context = NULL;
static char global_err_msg[1024] = "";

static void
s05_error_handler(const char *reason, const char *file, int line, int gsl_errno)
{
    /* Go ahead and segfault if global_context is NULL. */

    int n;

    n = snprintf(global_err_msg, COUNT(global_err_msg),
                 "%s (%s:%d; %d = %s)", reason, file, line, gsl_errno, gsl_strerror(gsl_errno));
    if ((size_t) n > COUNT(global_err_msg))
        n = COUNT(global_err_msg) - 1;
    global_err_msg[n] = '\0';

    global_context->last_error = RESULT_GSL_ERROR;
}

void
summers2005_debug_hook(void)
{
    printf("Debug hook called.\n");
}

static gsl_integration_workspace *integ_workspace = NULL; /* We're lame and leak these */
static gsl_poly_complex_workspace *poly5_workspace = NULL;
static gsl_poly_complex_workspace *poly7_workspace = NULL;

const size_t INTEG_WS_SIZE = 512;
const double integ_abs_tol = 0.0;
const double integ_rel_tol = 1e-6;

const int lambda = -1; /* signifies that we're looking at electrons */
const double epsilon = 0.0005446; /* electron-to-proton mass ratio */
const double epsm1 =  -0.9994554; /* epsilon - 1; C is lame and can't do the math with consts! */
const double sigma = 2.; /* range of `x` that we consider */
const double pi_on_2nu = 0.89039; /* pi/2nu, nu = sqrt(pi)erf(sigma); see after Shprits 06 Eqn 5 */

static inline result_t
apply_latitude(double latitude, state_t *state)
{
    int i;

    if (state->last_error != RESULT_OK)
        return state->last_error;

    state->latitude = latitude;

    double scl = cos(latitude); /* sin of colatitude = cos of latitude */
    double r = sqrt(4 - 3 * scl * scl) / pow(scl, 6); /* B(lam) / B_eq */

    state->Omega_e = state->p.Omega_e * r;
    state->sin_alpha = state->p.sin_alpha * sqrt(r);
    state->mu = sqrt(1 - state->sin_alpha * state->sin_alpha); /* cos(alpha) */

    if (state->mu == 0.) /* happens at the very top of our bounce trajectory */
        state->mu = 1e-20;

    /* Find the critical `x` and `y` values. We do this following Appendix A
     * and Equation 24 of Summers 2005. */

    state->n_xy = 0;

    double c[5] = { 0 };
    double z[8] = { 0 };

    const int s = state->p.handedness;
    const double a = state->a;
    const double a2 = a * a;
    const double bm2 = pow(state->beta * state->mu, 2);

    c[0] = -a2 * epsilon;
    c[1] = a2 * s * epsm1 - 2 * a * epsilon;
    c[2] = a2 + 2 * a * s * epsm1 - epsilon + bm2 * (state->b + epsilon);
    c[3] = 2 * a + s * epsm1 - bm2 * s * epsm1;
    c[4] = 1 - bm2;

    gsl_poly_complex_solve(c, 5, poly5_workspace, z);

    if (state->last_error != RESULT_OK) {
        strncat(global_err_msg, " (while finding x/y roots)", COUNT(global_err_msg) - 1);
        return state->last_error;
    }

    for (i = 0; i < 4; i++) {
        if (z[2*i + 1] != 0.) /* Imaginary root? */
            continue;

        double x = z[2*i];

        if (x < 0) /* Non-physical root? */
            continue;

        if (x < state->p.x_m - sigma * state->p.delta_x || x > state->p.x_m + sigma * state->p.delta_x)
            continue; /* Out of waveband? */

        double y = (x + state->a) / (state->beta * state->mu);

        if ((state->p.wave_filtering & WF_NO_FORWARD_MODES) && y > 0)
            continue;

        if ((state->p.wave_filtering & WF_NO_BACKWARD_MODES) && y < 0)
            continue;

        state->x[state->n_xy] = x;
        state->y[state->n_xy] = y;
        state->n_xy++;
    }

    if (state->n_xy > 3) {
        snprintf(global_err_msg, COUNT(global_err_msg), "expect 0 to 3 X/Y solutions; got %d", state->n_xy);
        return RESULT_X_FAILED;
    }

    /* Now we calculate F(x,y) = dx/dy, following Appendix C of Summers 2005,
     * and related values that contribute to the sums in the integrands. `Q`
     * is `1 - x cosa / y beta` and `R` is the factor that is in common
     * between all three calculations. */

    for (i = 0; i < state->n_xy; i++) {
        double x = state->x[i];
        double y = state->y[i];

        c[0] = epsilon * (state->b + epsilon);
        c[1] = -0.5 * s * epsm1 * (state->b + 4 * epsilon);
        c[2] = 1 - 4 * epsilon + epsilon * epsilon;
        c[3] = 2 * s * epsm1;
        c[4] = 1.;

        double g = gsl_poly_eval(c, 5, x);
        double F = y * pow((x - s) * (x + s * epsilon), 2) / (x * g);

        state->Q[i] = 1 - x * state->mu / (y * state->beta);

        double k0 = pow((x - state->p.x_m) / state->p.delta_x, 2);
        double k1 = fabs(state->beta * state->mu - F);

        if (k1 < 1e-5) {
            snprintf(global_err_msg, COUNT(global_err_msg), "hit a critical point (%f, %f, %f)",
                     state->sin_alpha, x, y);
            return RESULT_HIT_CRITICAL_POINT;
        }

        state->R[i] = state->p.R * fabs(F) * exp(-k0) / (state->p.delta_x * k1);
    }

    return RESULT_OK;
}

typedef double (*gsl_signature)(double, void *);

static double
daa_integrand(double latitude, state_t *state)
{
    result_t r;

    if (state->last_error != RESULT_OK)
        return 0;

    if ((r = apply_latitude(latitude, state)) != RESULT_OK) {
        state->last_error = r;
        return 0;
    }

    int i;
    double v = 0;

    for (i = 0; i < state->n_xy; i++)
        v += state->R[i] * pow(state->Q[i], 2);

    if (state->itype == ITYPE_BOUNCE_AVERAGED)
        v *= state->mu * pow(cos(state->latitude), 7);

    return v * state->Omega_e;
}

static double
dap_on_p_integrand(double latitude, state_t *state)
{
    result_t r;

    if (state->last_error != RESULT_OK)
        return 0;

    if ((r = apply_latitude(latitude, state)) != RESULT_OK) {
        state->last_error = r;
        return 0;
    }

    int i;
    double v = 0;

    for (i = 0; i < state->n_xy; i++)
        v += state->R[i] * state->Q[i] * state->x[i] / state->y[i];

    /* Here we transform the sqrt(1 + 3 sin^2 lambda) term to use the cos
     * instead */

    double cos_lam = cos(state->latitude);
    double sa = state->sin_alpha;

    if (state->itype == ITYPE_BOUNCE_AVERAGED)
        v *= state->mu * cos_lam * sqrt(4 - 3 * cos_lam * cos_lam) / sa;

    return v * state->Omega_e * sa;
}

static double
dpp_on_p2_integrand(double latitude, state_t *state)
{
    result_t r;

    if (state->last_error != RESULT_OK)
        return 0;

    if ((r = apply_latitude(latitude, state)) != RESULT_OK) {
        state->last_error = r;
        return 0;
    }

    int i;
    double v = 0;

    for (i = 0; i < state->n_xy; i++)
        v += state->R[i] * pow(state->x[i] / state->y[i], 2);

    /* Same transform as in dap_on_p. */

    double cos_lam = cos(state->latitude);
    double sa = state->sin_alpha;

    if (state->itype == ITYPE_BOUNCE_AVERAGED)
        v *= cos_lam * sqrt(4 - 3 * cos_lam * cos_lam) / state->mu;

    return v * state->Omega_e * sa * sa;
}

static result_t
calc_coefficients(parameters_t *params, coefficients_t *coeffs)
{
    gsl_error_handler_t *prev_handler;
    gsl_function integrand;
    state_t state;
    double c[7] = { 0 };
    double z[12] = { 0 };
    double lambda_m = -1;
    double result, abserr, k;
    int i;

    if (params->sin_alpha > 0.99985) /* ~= 89 degrees */
        params->sin_alpha = 0.99985;

    global_context = &state;
    prev_handler = gsl_set_error_handler(s05_error_handler);

    if (integ_workspace == NULL) {
        integ_workspace = gsl_integration_workspace_alloc(INTEG_WS_SIZE);
        poly5_workspace = gsl_poly_complex_workspace_alloc(5);
        poly7_workspace = gsl_poly_complex_workspace_alloc(7);
    }

    state.p = *params;
    state.itype = ITYPE_BOUNCE_AVERAGED;
    state.last_error = RESULT_OK;

    /* Figure out the mirroring latitude; Shprits (2006) equation 10. */

    double f0 = pow(state.p.sin_alpha, 4);
    c[0] = -4 * f0;
    c[1] = 3 * f0;
    c[6] = 1.;

    gsl_poly_complex_solve(c, 7, poly7_workspace, z);

    if (state.last_error != RESULT_OK) {
        strncat(global_err_msg, " (while finding lambda_m)", COUNT(global_err_msg) - 1);
        return state.last_error;
    }

    for (i = 0; i < 6; i++) {
        if (z[2*i + 1] == 0. && z[2*i] > 0) {
            lambda_m = acos(sqrt(z[2*i]));
            break;
        }
    }

    if (!(lambda_m >= 0 && lambda_m <= M_PI_2)) {
        gsl_set_error_handler(prev_handler);
        global_context = NULL;
        snprintf(global_err_msg, COUNT(global_err_msg), "failed to compute lambda_m (%.16lf)", lambda_m);
        return RESULT_LAMBDA_FAILED;
    }

    /* It can happen that lambda_m is just a tiiiiiny bit too large such that
     * when we compute the sine of the pitch angle at lambda_m, we get
     * something bigger than one. Witness our beautiful workaround.
     */

    double cos_lat = cos(lambda_m);

    if (sqrt(4 - 3 * cos_lat * cos_lat) / (pow(cos_lat, 6) * state.p.sin_alpha * state.p.sin_alpha) > 1.)
        lambda_m -= 1e-7;

    /* If a latitude limit on waves is imposed, truncate lambda_m further. */

    if (lambda_m > state.p.max_wave_latitude)
        lambda_m = state.p.max_wave_latitude;

    /* Pre-compute quantities that do not depend on lambda. When lambda
     * changes, B and alpha change. Consequently Omega_e changes too. In our
     * model the other parameters stay fixed.
     */

    state.gamma = state.p.E + 1;
    state.a = ((int) state.p.handedness) * lambda / state.gamma;
    state.b = (1 + epsilon) / state.p.alpha_star;
    state.beta = sqrt(state.p.E * (state.p.E + 2)) / (state.p.E + 1);
    coeffs->dimensionless_p = state.gamma * state.beta;

    /* Shprits 2006 eqn 9, credited to Lencheck+ 1971 and Shultz & Lanzerotti 1974: */

    double s0 = 1.38 - 0.32 * (params->sin_alpha + sqrt(params->sin_alpha));
    double f1 = pi_on_2nu / (pow(state.p.E + 1, 2) * s0);

    /* Let's integrate! First, D_aa. */

    integrand.function = (gsl_signature) daa_integrand;
    integrand.params = &state;

    gsl_integration_qag(
        &integrand,
        0., /* lower limit */
        lambda_m, /* upper limit */
        integ_abs_tol,
        integ_rel_tol,
        INTEG_WS_SIZE, /* allocated workspace size */
        GSL_INTEG_GAUSS51, /* integration rule */
        integ_workspace, /* workspace */
        &result, /* output: value of the integral */
        &abserr /* estimated absolute error */
    );

    if (state.last_error != RESULT_OK) {
        strncat(global_err_msg, " (while computing Daa)", COUNT(global_err_msg) - 1);
        gsl_set_error_handler(prev_handler);
        global_context = NULL;
        return state.last_error;
    }

    k = f1 / (1 - state.p.sin_alpha * state.p.sin_alpha);
    coeffs->Daa = result * k;
    coeffs->err_Daa = abserr * k;

    /* D_ap/p. */

    integrand.function = (gsl_signature) dap_on_p_integrand;

    gsl_integration_qag(
        &integrand,
        0., /* lower limit */
        lambda_m, /* upper limit */
        integ_abs_tol,
        integ_rel_tol,
        INTEG_WS_SIZE, /* allocated workspace size */
        GSL_INTEG_GAUSS51, /* integration rule */
        integ_workspace, /* workspace */
        &result, /* output: value of the integral */
        &abserr /* estimated absolute error */
    );

    if (state.last_error != RESULT_OK) {
        strncat(global_err_msg, " (while computing Dap)", COUNT(global_err_msg) - 1);
        gsl_set_error_handler(prev_handler);
        global_context = NULL;
        return state.last_error;
    }

    k = -f1 * state.p.sin_alpha / (state.beta * sqrt(1 - state.p.sin_alpha * state.p.sin_alpha));
    coeffs->Dap_on_p = result * k;
    coeffs->err_Dap_on_p = abserr * fabs(k);

    /* D_pp/p^2. */

    integrand.function = (gsl_signature) dpp_on_p2_integrand;

    gsl_integration_qag(
        &integrand,
        0., /* lower limit */
        lambda_m, /* upper limit */
        integ_abs_tol,
        integ_rel_tol,
        INTEG_WS_SIZE, /* allocated workspace size */
        GSL_INTEG_GAUSS51, /* integration rule */
        integ_workspace, /* workspace */
        &result, /* output: value of the integral */
        &abserr /* estimated absolute error */
    );

    if (state.last_error != RESULT_OK) {
        strncat(global_err_msg, " (while computing Dpp)", COUNT(global_err_msg) - 1);
        gsl_set_error_handler(prev_handler);
        global_context = NULL;
        return state.last_error;
    }

    k = f1 * pow(state.beta, -2);
    coeffs->Dpp_on_p2 = result * k;
    coeffs->err_Dpp_on_p2 = abserr * k;

    /* Huzzah, all done. */

    gsl_set_error_handler(prev_handler);
    global_context = NULL;
    return RESULT_OK;
}

static result_t
calc_unaveraged_coefficients(parameters_t *params, coefficients_t *coeffs)
{
    gsl_error_handler_t *prev_handler;
    state_t state;

    global_context = &state;
    prev_handler = gsl_set_error_handler(s05_error_handler);

    if (integ_workspace == NULL) {
        integ_workspace = gsl_integration_workspace_alloc(INTEG_WS_SIZE);
        poly5_workspace = gsl_poly_complex_workspace_alloc(5);
        poly7_workspace = gsl_poly_complex_workspace_alloc(7);
    }

    state.p = *params;
    state.itype = ITYPE_LOCAL;
    state.last_error = RESULT_OK;

    /* The structure here is paralleling calc_coefficients for maintainability. */

    state.gamma = state.p.E + 1;
    state.a = ((int) state.p.handedness) * lambda / state.gamma;
    state.b = (1 + epsilon) / state.p.alpha_star;
    state.beta = sqrt(state.p.E * (state.p.E + 2)) / (state.p.E + 1);
    coeffs->dimensionless_p = state.gamma * state.beta;

    double f1 = pi_on_2nu / pow(state.p.E + 1, 2);

    /* Compute! It all works if we pretend latitude = 0 and set `itype` to LOCAL. */

    coeffs->Daa = daa_integrand(0., &state);

    if (state.last_error != RESULT_OK) {
        strncat(global_err_msg, " (while computing Daa)", COUNT(global_err_msg) - 1);
        gsl_set_error_handler(prev_handler);
        global_context = NULL;
        return state.last_error;
    }

    coeffs->Daa *= f1;
    coeffs->err_Daa = 0;

    /* D_ap/p. */

    coeffs->Dap_on_p = dap_on_p_integrand(0., &state);

    if (state.last_error != RESULT_OK) {
        strncat(global_err_msg, " (while computing Dap)", COUNT(global_err_msg) - 1);
        gsl_set_error_handler(prev_handler);
        global_context = NULL;
        return state.last_error;
    }

    coeffs->Dap_on_p *= f1 / state.beta;
    coeffs->err_Dap_on_p = 0;

    /* D_pp/p^2. */

    coeffs->Dpp_on_p2 = dpp_on_p2_integrand(0., &state);

    if (state.last_error != RESULT_OK) {
        strncat(global_err_msg, " (while computing Dpp)", COUNT(global_err_msg) - 1);
        gsl_set_error_handler(prev_handler);
        global_context = NULL;
        return state.last_error;
    }

    coeffs->Dpp_on_p2 *= f1 * pow(state.beta, -2);
    coeffs->err_Dpp_on_p2 = 0.;

    /* Huzzah, all done. */

    gsl_set_error_handler(prev_handler);
    global_context = NULL;
    return RESULT_OK;
}

static PyObject*
get_coeffs(PyObject *self, PyObject* args)
{
    int modespec, handspec, wave_filtering;
    parameters_t params;
    coefficients_t coeffs = { 0 };
    result_t r;

    if (!PyArg_ParseTuple(args, "iiidddddddd", &modespec, &handspec, &wave_filtering,
                          &params.E,
                          &params.sin_alpha,
                          &params.Omega_e,
                          &params.alpha_star,
                          &params.R,
                          &params.x_m,
                          &params.delta_x,
                          &params.max_wave_latitude))
        return NULL;

    if (handspec == 0)
        params.handedness = R_MODE;
    else if (handspec == 1)
        params.handedness = L_MODE;
    else {
        PyErr_SetString(PyExc_RuntimeError, "unexpected handedness magic constant");
        return NULL;
    }

    params.wave_filtering = (wave_filtering_t) wave_filtering;

    if (modespec == 0)
        r = calc_coefficients(&params, &coeffs);
    else if (modespec == 1)
        r = calc_unaveraged_coefficients(&params, &coeffs);
    else {
        PyErr_SetString(PyExc_RuntimeError, "unexpected mode magic constant");
        return NULL;
    }

    if (r != RESULT_OK) {
        PyErr_SetString(PyExc_RuntimeError, global_err_msg);
        return NULL;
    }

    return Py_BuildValue("ddddddd", coeffs.dimensionless_p,
                         coeffs.Daa, coeffs.err_Daa,
                         coeffs.Dap_on_p, coeffs.err_Dap_on_p,
                         coeffs.Dpp_on_p2, coeffs.err_Dpp_on_p2);
}

static PyMethodDef methods[] = {
    { "get_coeffs", get_coeffs, METH_VARARGS|METH_KEYWORDS,
      "(mode, handedness, wfilt, E, sa, Oe, a*, R, xm, dx, mwl) -> "
      "(dp, daa, udaa, dap/p, udap/p, dpp/p2, udpp/p2)" },
    { NULL, NULL, METH_NOARGS, NULL },
};

#if PY_MAJOR_VERSION >= 3
static struct PyModuleDef module_def = {
    PyModuleDef_HEAD_INIT,
    "_summers2005",
    NULL,
    0,
    methods,
    NULL,
    NULL,
    NULL,
    NULL
};

# define INIT_RET_TYPE PyObject *
#else
# define INIT_RET_TYPE void
#endif

INIT_RET_TYPE
PyInit__summers2005(void)
{
#if PY_MAJOR_VERSION >= 3
    return PyModule_Create(&module_def);
#else
    PyImport_AddModule("_summers2005");
    Py_InitModule("_summers2005", methods);
#endif
}

max_stars_repo_name: pkgw/vernon

max_issues_repo_name: pkgw/vernon

max_forks_repo_name: pkgw/vernon


text:
// for license information, see the accompanying LICENSE file

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <assert.h>

// CBLAS, LAPACKE are f2c translated Netlib F77 reference routines.
// 'cblas.h' is widely supported, and is the f2c translated Netlib F77 refernce version.
// There are known issues with name mangling schemes regarding wrapped code and inter-language operability between Fortran and C.
//
// Netlib CBLAS / LAPACKE target
// #define LISE_LA_REF
// #include <lapacke.h>
// #include <cblas.h>
//
// IBM ESSL target
// #define LISE_LA_ESSL
// #include <essl.h>
//
// Intel MKL target
// #define LISE_LA_MKL
// #include <mkl.h>

#ifdef LISE_LA_MKL
#include <mkl.h>
double ddot(const int *, const double *, const int *, const double *, const int *);
void daxpy(const int *, const double *, const double *, const int *, double *, const int *);
void dgemv(const char *, const int *, const int *, const double *, const double *, const int *, const double *, const int *, const double *, double *, const int *);
void dcopy(const int *, const double *, const int *, double *, const int *);
void dgemm(const char *, const char *, const int *, const int *, const int *, const double *, const double *, const int *, const double *, const int *, const double *, double *, const int *);
void dgesdd(const char *, const int *, const int *, double *, const int *, double *, double *, const int *, double *, const int *, double *, const int *, int *, int *);
void dger(const int *, const int *, const double *, const double *, const int *, const double *, const int *, double *, const int *);
#endif

#ifdef LISE_LA_ESSL
#include <essl.h>
double ddot(int, const double *, int, const double *,int);
void daxpy(int,  double,  double *, int,  double *,  int);
void dgemv(const char *, int, int,  double, const void *, int, const  double *, int,  double,  double *, int);
void dcopy(int,  double *, int,  double *, int);
void dgemm(const char *, const char *, int, int, int,  double, const void *, int, const void *, int,  double, void *, int);
void dgesdd(const char *, int, int, void *, int, double *, void *, int, void *, int, double *, int, int *, int *);
void dger(int, int,  double, const  double *, int, const  double *, int, void *, int);
// void dger1(int, int,  double, const  double *, int, const  double *, int, void *, int);
#endif

#ifdef LISE_LA_REF
#include <lapacke.h>
#include <cblas.h>
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY);
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY);
void cblas_dgemv(CBLAS_LAYOUT layout, CBLAS_TRANSPOSE TransA, const int M, const int N, const double alpha, const double *A, const int lda, const double *X, const int incX, const double beta, double *Y, const int incY);
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY);
void cblas_dgemm(CBLAS_LAYOUT layout, CBLAS_TRANSPOSE TransA, CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const double alpha, const double *A, const int lda, const double *B, const int ldb, const double beta, double *C, const int ldc);
void cblas_dger(CBLAS_LAYOUT layout, const int M, const int N, const double alpha, const double *X, const int incX, const double *Y, const int incY, double *A, const int lda);
lapack_int LAPACKE_dgesdd( int matrix_layout, char jobz, lapack_int m, lapack_int n, double* a, lapack_int lda, double* s, double* u, lapack_int ldu, double* vt, lapack_int ldvt );
#endif

double ddotII( int * , double * , int * , double * , int * ) ;
void daxpyII( int * , double * , double * , int * , double * , int * ) ;
int broyden_minB( double * v_in , double * v_in_old , double * f_old , double * b , double * v_out , const int n , int * it_out , const double alpha ) ;
int broyden_min( double * v_in , double * v_in_old , double * f_old , double * b_bra , double * b_ket , double * d , double * v_out , const int n , int * it_out , const double alpha ) ;
int broydenMod_min( double * v_in , double * v_in_old , double * f_old , double * b_bra , double * b_ket , double * d , double * v_out , const int n , int * it_out , const double alpha ) ;

void broyden_min_( double * v_in , double * v_in_old , double * f_old , double * b_bra , double * b_ket , double * d , double * v_out , int * n , int * it_out , double * alpha )

{
  broydenMod_min( v_in , v_in_old , f_old , b_bra , b_ket , d , v_out , * n , it_out , * alpha ) ;
}

void broyden_minb_( double * v_in , double * v_in_old , double * f_old , double * b_bra , double * v_out , int * n , int * it_out , double * alpha )

{

  broyden_minB( v_in , v_in_old , f_old , b_bra , v_out , * n , it_out , * alpha ) ;

}

int broyden_min( double * v_in , double * v_in_old , double * f_old , double * b_bra , double * b_ket , double * d , double * v_out , const int n , int * it_out , const double alpha )

{

  double * f , * dv , * df ;

  int i , ii , it1 ;

  int ione = 1 ;

  double norm , norm1 , norm2 ;

  int n1 , ishift ;

  int it ;

  it = * it_out ;

  n1 = n ;

  if( it == 0 )

    {

      * it_out = 1 ;

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

	v_in_old[ i ] = v_in[ i ] ;

	f_old[ i ] = alpha * ( v_in[ i ] - v_out[ i ] ) ;

	v_in[ i ]  -= f_old[ i ] ;

	v_out[ i ] = v_in[ i ] ;

      }

      return( 0 ) ;

    }

  assert( f = malloc( n * sizeof( double ) ) ) ;

  assert( df = malloc( n * sizeof( double ) ) ) ;

  assert( dv = malloc( n * sizeof( double ) ) ) ;

  it1 = it - 1 ;

  ishift = it1 * n ;

  norm = 0. ; /* || F || */

  norm1 = 0. ; /* || F_old || */

  for( i = 0 ; i < n ; i++ )

    {

      f[ i ] =  alpha * ( v_in[ i ] - v_out[ i ] ) ;

      norm += f[ i ] * f[ i ] ;

      norm1 += f_old[ i ] * f_old[ i ] ;

      dv[ i ] = v_in[ i ] - v_in_old[ i ] ;

      df[ i ] = f[ i ] - f_old[ i ] ;

      b_ket[ ishift + i ] = df[ i ] ;

      b_bra[ ishift + i ] = dv[ i ] ;

    }

  for( ii = 0 ; ii < it1 ; ++ii )

    {

#ifdef LISE_LA_MKL
      norm =  d[ ii ] * ddot( &n1 , df , &ione , b_bra + ii * n , &ione ) ;
      daxpy( &n1 , &norm , b_ket + ii * n , &ione , b_ket + ishift , &ione ) ;
      norm = d[ ii ] * ddot( &n1 , dv , &ione , b_ket + ii * n , &ione ) ;
      daxpy( &n1 , &norm , b_bra + ii * n , &ione , b_bra + ishift ,&ione ) ;
#endif
#ifdef LISE_LA_ESSL
      norm =  d[ ii ] * ddot( n1 , df , ione , b_bra + ii * n , ione ) ;
      daxpy( n1 , norm , b_ket + ii * n , ione , b_ket + ishift , ione ) ;
      norm = d[ ii ] * ddot( n1 , dv , ione , b_ket + ii * n , ione ) ;
      daxpy( n1 , norm , b_bra + ii * n , ione , b_bra + ishift ,ione ) ;
#endif
#ifdef LISE_LA_REF
//double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY);
//void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY);
      norm =  d[ ii ] * cblas_ddot( n1 , df , ione , b_bra + ii * n , ione ) ;
      cblas_daxpy( n1 , norm , b_ket + ii * n , ione , b_ket + ishift , ione ) ;
      norm = d[ ii ] * cblas_ddot( n1 , dv , ione , b_ket + ii * n , ione ) ;
      cblas_daxpy( n1 , norm , b_bra + ii * n , ione , b_bra + ishift ,ione ) ;
#endif

    }

#ifdef LISE_LA_MKL
  norm = ddot(  &n1 , dv , &ione , b_ket + ishift , &ione ) ;
#endif
#ifdef LISE_LA_ESSL
  norm = ddot(  n1 , dv , ione , b_ket + ishift , ione ) ;
#endif
#ifdef LISE_LA_REF
  norm = cblas_ddot(  n1 , dv , ione , b_ket + ishift , ione ) ;
#endif

  printf( "<dv|B|df> =%12.6e\n", norm ) ;

  if( fabs( norm ) < 1.e-15 ) return( 0 ) ;

  for( i = 0 ; i < n ; ++i )

    b_ket[ ishift + i ] = ( dv[ i ] - b_ket[ ishift + i ] ) / norm ;

#ifdef LISE_LA_MKL
  norm1 = sqrt( ddot( &n1 , b_ket + ishift , &ione , b_ket + ishift , &ione ) ) ;
  norm2 = sqrt( ddot( &n1 , b_bra + ishift , &ione , b_bra + ishift , &ione ) ) ;
#endif
#ifdef LISE_LA_ESSL
  norm1 = sqrt( ddot( n1 , b_ket + ishift , ione , b_ket + ishift , ione ) ) ;
  norm2 = sqrt( ddot( n1 , b_bra + ishift , ione , b_bra + ishift , ione ) ) ;
#endif
#ifdef LISE_LA_REF
  norm1 = sqrt( cblas_ddot( n1 , b_ket + ishift , ione , b_ket + ishift , ione ) ) ;
  norm2 = sqrt( cblas_ddot( n1 , b_bra + ishift , ione , b_bra + ishift , ione ) ) ;
#endif

  for( i = 0 ; i < n ; i++ )

    {

      b_ket[ ishift + i ] = b_ket[ ishift + i ] / norm1 ;

      b_bra[ ishift + i ] = b_bra[ ishift + i ] / norm2 ;

    }

  d[ it1 ] = norm1 * norm2 ;

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

      v_in_old[ i ] = v_in[ i ] ;

      v_in[ i ] -= f[ i ] ;

    }

  for( ii = 0 ; ii < it ; ii++ )

    {

#ifdef LISE_LA_MKL
      norm = - d[ ii ] * ddot( &n1 , b_bra + ii * n , &ione , f , &ione ) ;
      daxpy( &n1 , &norm , b_ket + ii * n , &ione , v_in , &ione ) ;
#endif
#ifdef LISE_LA_ESSL
      norm = - d[ ii ] * ddot( n1 , b_bra + ii * n , ione , f , ione ) ;
      daxpy( n1 , norm , b_ket + ii * n , ione , v_in , ione ) ;
#endif
#ifdef LISE_LA_REF
      norm = - d[ ii ] * cblas_ddot( n1 , b_bra + ii * n , ione , f , ione ) ;
      cblas_daxpy( n1 , norm , b_ket + ii * n , ione , v_in , ione ) ;
#endif

    }

  for( i = 0 ; i < n ; i++ )

    {

      f_old[ i ] = f[ i ] ;

      v_out[ i ] = v_in[ i ] ;

    }

  ( * it_out ) ++ ;

  free( f ) ; free( dv ) ; free( df ) ;

  return( 0 ) ;

}

int broyden_minB( double * v_in , double * v_in_old , double * f_old , double * b , double * v_out , const int n , int * it_out , const double alpha )

{

  double * f , * dv , * df , * bdf , * bdv ;

  int i , ii , it1 ;

  int ione = 1 ;

  double norm , norm1 , norm2 ;

  int n1 , ishift ;

  int it , j ;

  double * u , * vt , * s , * work , wk_[ 2 ] ;

  int lwork = -1 , info , * iwork ;

  it = * it_out ;

  printf( "Starting Broyden %d \n\n" , it ) ;

  n1 = n ;

  if( it == 0 )

    {

      * it_out = 1 ;

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

	v_in_old[ i ] = v_in[ i ] ;

	f_old[ i ] = alpha * ( v_in[ i ] - v_out[ i ] ) ;

	v_in[ i ]  -= f_old[ i ] ;

	v_out[ i ] = v_in[ i ] ;

	for( j = 0 ; j < n ; j++ )

	  b[ i + n * j ] = 0. ;

	b[ i + n * i ] = 1. ;

      }

      return( 0 ) ;

    }

  assert( f = malloc( n * sizeof( double ) ) ) ;

  assert( df = malloc( n * sizeof( double ) ) ) ;

  assert( dv = malloc( n * sizeof( double ) ) ) ;

  assert( bdv = malloc( n * sizeof( double ) ) ) ;

  assert( bdf = malloc( n * sizeof( double ) ) ) ;

  it1 = it - 1 ;

  ishift = it1 * n ;

  for( i = 0 ; i < n ; i++ )

    {

      f[ i ] =  alpha * ( v_in[ i ] - v_out[ i ] ) ;

      dv[ i ] = v_in[ i ] - v_in_old[ i ] ;

      df[ i ] = f[ i ] - f_old[ i ] ;

    }

  norm = 0. ;

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

    bdv[ i ] = 0. ;

    bdf[ i ] = dv[ i ] ;

    for( j = 0 ; j < n ; j++ ){

      bdv[ i ] += b[ j + n * i ] * dv[ j ] ;

      bdf[ i ] -= b[ i + n * j ] * df[ j ] ;

    }

    norm += df[ i ] * bdv[ i ] ;

  }

  printf( "<dv|B|df> =%12.6e\n", norm ) ;

  norm = 1./norm ;

#ifdef LISE_LA_REF
  cblas_dger( CBLASColMajor, n1 , n1 , norm , bdf , ione , bdv , ione , b , n1 ) ;
#endif
#ifdef LISE_LA_ESSL
  dger( n1 , n1 , norm , bdf , ione , bdv , ione , b , n1 ) ;
#endif
#ifdef LISE_LA_MKL
// void dger(const int *, const int *, const double *, const double *, const int *, const double *, const int *, void *, const int *);
  dger( &n1 , &n1 , &norm , bdf , &ione , bdv , &ione , b , &n1 ) ;
#endif

  if( it > 5 ){

  /* add a SVD decomposition of matrix B */

    assert( s = malloc( n * sizeof( double ) ) ) ;

    assert( u = malloc( n * n  * sizeof( double ) ) ) ;

    assert( vt = malloc( n * n * sizeof( double ) ) ) ;

    assert( iwork = malloc( 8 * n * sizeof( double ) ) ) ;

#ifdef LISE_LA_MKL
    dgesdd( "A", &n1 , &n1 , b , &n1 , s , u , &n1 , vt , &n1 , wk_ , &lwork , iwork , &info ) ;
#endif
#ifdef LISE_LA_ESSL
    dgesdd( "A", n1 , n1 , b , n1 , s , u , n1 , vt , n1 , wk_ , lwork , iwork , &info ) ;
#endif
#ifdef LISE_LA_REF
    info = LAPACKE_dgesdd(LAPACK_COL_MAJOR, 'A', n1 , n1 , b , n1 , s , u , n1 , vt , n1 ) ;
#endif

    lwork = ( int ) wk_[ 0 ] ;

    assert( work = malloc( lwork * sizeof( double ) ) ) ;

#ifdef LISE_LA_ESSL
    dgesdd( "A" , n1 , n1 , b , n1 , s , u , n1 , vt , n1 , work , lwork , iwork , &info ) ;
#endif
#ifdef LISE_LA_REF
    info = LAPACKE_dgesdd(LAPACK_COL_MAJOR, 'A', n1 , n1 , b , n1 , s , u , n1 , vt , n1 ) ;
#endif
#ifdef LISE_LA_MKL
    dgesdd( "A", &n1 , &n1 , b , &n1 , s , u , &n1 , vt , &n1 , work , &lwork , iwork , &info ) ;
#endif

    free( work ) ; free( iwork ) ;

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

      if( s[ i ] < 1.e-6 ) {

	printf( " s[ %d ] = %g \n" , i , s[ i ] ) ;

	s[ i ] = 0. ;

      }

      for( j = 0 ; j < n ; j++ )

	vt[ i + j * n ] *= s[ i ] ;

    }

    free( s ) ;

    norm = 0. ;

    norm1 = 1. ;

#ifdef LISE_LA_ESSL
    dgemm( "N" , "N" , n1 , n1 , n1 , norm1 , u , n1 , vt , n1 , norm , b , n1 ) ;
#endif
#ifdef LISE_LA_MKL
// void dgemm(const char *, const char *, const int *, const int *, const int *, const double *, const double *, const int *, const double *, const int *, const double *, double *, const int *);
    dgemm( "N" , "N" , &n1 , &n1 , &n1 , &norm1 , u , &n1 , vt , &n1 , &norm , b , &n1 ) ;
#endif
#ifdef LISE_LA_REF
    cblas_dgemm(CBLAS_COL_MAJOR,"N" , "N" , n1 , n1 , n1 , norm1 , u , n1 , vt , n1 , norm , b , n1);
#endif

    free( u ) ; free( vt ) ;

  }

#ifdef LISE_LA_ESSL
  dcopy( n1 , v_in , ione , v_in_old , ione ) ;
#endif
#ifdef LISE_LA_MKL
//void dcopy(const int *, const double *, const int *, double *, const int *);
  dcopy( &n1 , v_in , &ione , v_in_old , &ione ) ;
#endif
#ifdef LISE_LA_REF
  cblas_dcopy( n1 , v_in , ione , v_in_old , ione ) ;
#endif

  norm1 = -1. ;

  norm = 1. ;

#ifdef LISE_LA_MKL
//void dgemv(const char *, const int *, const int *, const double *, const double *, const int *, const double *, const int *, const double *, double *, const int *);
  dgemv( "N" , &n1 , &n1 , &norm1 , b , &n1 , f , &ione , &norm , v_in , &ione ) ;
#endif
#ifdef LISE_LA_ESSL
  dgemv( "N" , n1 , n1 , norm1 , b , n1 , f , ione , norm , v_in , ione ) ;
#endif
#ifdef LISE_LA_REF
  cblas_dgemv(CBLAS_COL_MAJOR, "N" , n1 , n1 , norm1 , b , n1 , f , ione , norm , v_in , ione ) ;
#endif

  for( i = 0 ; i < n ; i++ )

    {

      f_old[ i ] = f[ i ] ;

      v_out[ i ] = v_in[ i ] ;

    }

  ( * it_out ) ++ ;

  free( f ) ; free( dv ) ; free( df ) ; free( bdf ) ; free( bdv ) ;

  printf( "Done with Broyden \n\n" ) ;

  return( 0 ) ;

}

int broydenMod_min( double * v_in , double * v_in_old , double * f_old , double * b_bra , double * b_ket , double * d , double * v_out , const int n , int * it_out , const double alpha )
/*

Change in the formula. Use Eq. (10) in J.Chem.Phys. 134 (2011) 134109

 */

{

  double * f , * dv , * df ;

  int i , ii , it1 ;

  int ione = 1 ;

  double norm , norm1 , norm2 ;

  int n1 , ishift ;

  int it ;

  it = * it_out ;

  printf( "Starting Broyden %d \n\n" , it ) ;

  n1 = n ;

  if( it == 0 )

    {

      * it_out = 1 ;

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

	v_in_old[ i ] = v_in[ i ] ;

	f_old[ i ] = alpha * ( v_in[ i ] - v_out[ i ] ) ;

	v_in[ i ]  -= f_old[ i ] ;

	v_out[ i ] = v_in[ i ] ;

      }

      return( 0 ) ;

    }

  assert( f = malloc( n * sizeof( double ) ) ) ;

  assert( df = malloc( n * sizeof( double ) ) ) ;

  assert( dv = malloc( n * sizeof( double ) ) ) ;

  it1 = it - 1 ;

  ishift = it1 * n ;

  norm = 0. ; /* || F || */

  norm1 = 0. ; /* || F_old || */

  for( i = 0 ; i < n ; i++ )

    {

      f[ i ] =  alpha * ( v_in[ i ] - v_out[ i ] ) ;

      norm += f[ i ] * f[ i ] ;

      norm1 += f_old[ i ] * f_old[ i ] ;

      dv[ i ] = v_in[ i ] - v_in_old[ i ] ;

      df[ i ] = f[ i ] - f_old[ i ] ;

      b_ket[ ishift + i ] = df[ i ] ;

      b_bra[ ishift + i ] = df[ i ] ;

    }

  for( ii = 0 ; ii < it1 ; ++ii )

    {

#ifdef LISE_LA_MKL
      norm =  d[ ii ] * ddot( &n1 , df , &ione , b_bra + ii * n , &ione ) ;
      daxpy( &n1 , &norm , b_ket + ii * n , &ione , b_ket + ishift , &ione ) ;
#endif
#ifdef LISE_LA_ESSL
      norm =  d[ ii ] * ddot( n1 , df , ione , b_bra + ii * n , ione ) ;
      daxpy( n1 , norm , b_ket + ii * n , ione , b_ket + ishift , ione ) ;
#endif
#ifdef LISE_LA_REF
      norm =  d[ ii ] * cblas_ddot( n1 , df , ione , b_bra + ii * n , ione ) ;
      cblas_daxpy( n1 , norm , b_ket + ii * n , ione , b_ket + ishift , ione ) ;
#endif

    }

#ifdef LISE_LA_MKL
  norm = ddot( &n1 , df , &ione , df , &ione ) ;
#endif
#ifdef LISE_LA_ESSL
  norm = ddot( n1 , df , ione , df , ione ) ;
#endif
#ifdef LISE_LA_REF
  norm = cblas_ddot( n1 , df , ione , df , ione ) ;
#endif

  printf( "<df|df> =%12.6e\n", norm ) ;

  for( i = 0 ; i < n ; ++i )

    b_ket[ ishift + i ] = ( dv[ i ] - b_ket[ ishift + i ] ) / norm ;

#ifdef LISE_LA_MKL
  norm1 = sqrt( ddot( &n1 , b_ket + ishift , &ione , b_ket + ishift , &ione ) ) ;
  norm2 = sqrt( ddot( &n1 , b_bra + ishift , &ione , b_bra + ishift ,&ione ) ) ;
#endif
#ifdef LISE_LA_ESSL
  norm1 = sqrt( ddot( n1 , b_ket + ishift , ione , b_ket + ishift , ione ) ) ;
  norm2 = sqrt( ddot( n1 , b_bra + ishift , ione , b_bra + ishift ,ione ) ) ;
#endif
#ifdef LISE_LA_REF
  norm1 = sqrt( cblas_ddot( n1 , b_ket + ishift , ione , b_ket + ishift , ione ) ) ;
  norm2 = sqrt( cblas_ddot( n1 , b_bra + ishift , ione , b_bra + ishift ,ione ) ) ;
#endif

  for( i = 0 ; i < n ; i++ )

    {

      b_ket[ ishift + i ] = b_ket[ ishift + i ] / norm1 ;

      b_bra[ ishift + i ] = b_bra[ ishift + i ] / norm2 ;

    }

  d[ it1 ] = norm1 * norm2 ;

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

      v_in_old[ i ] = v_in[ i ] ;

      v_in[ i ] -= f[ i ] ;

    }

  for( ii = 0 ; ii < it ; ii++ )

    {

#ifdef LISE_LA_MKL
      norm = - d[ ii ] * ddot( &n1 , b_bra + ii * n , &ione , f , &ione ) ;
      daxpy( &n1 , &norm , b_ket + ii * n , &ione , v_in , &ione ) ;
#endif
#ifdef LISE_LA_ESSL
      norm = - d[ ii ] * ddot( n1 , b_bra + ii * n , ione , f , ione ) ;
      daxpy( n1 , norm , b_ket + ii * n , ione , v_in , ione ) ;
#endif
#ifdef LISE_LA_REF
      norm = - d[ ii ] * cblas_ddot( n1 , b_bra + ii * n , ione , f , ione ) ;
      cblas_daxpy( n1 , norm , b_ket + ii * n , ione , v_in , ione ) ;
#endif

    }

  for( i = 0 ; i < n ; i++ )

    {

      f_old[ i ] = f[ i ] ;

      v_out[ i ] = v_in[ i ] ;

    }

  ( * it_out ) ++ ;

  free( f ) ; free( dv ) ; free( df ) ;

  printf( "Done with Broyden \n\n" ) ;

  return( 0 ) ;

}

double ddotII( int * n , double *v1 , int * i1 , double * v2 , int * i2 )

{

  double sum = 0. ;

  int i ;

  for( i = 0 ; i < * n ; i++ )

    sum += v1[ i ] * v2[ i ] ;

  return( sum ) ;

}

void daxpyII( int * n , double * alpha , double * x , int * i1 , double * y , int * i2 )

{

  int i ;

  for( i = 0 ; i < * n ; i++ )

    y[ i ] += *alpha * x[ i ] ;

}

max_stars_repo_name: dmageeLANL/LISE

max_issues_repo_name: dmageeLANL/LISE

max_forks_repo_name: dmageeLANL/LISE


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

double* default_vector;
double* default_matrix;
double** results;

int max(int a, int b){
  if(a > b) return a;
  return b;
}

void save_res_to_csv(double** results, size_t results_c, int stride, int trials, char* file_name){
  FILE* csv = fopen(file_name, "w");
  if(csv == NULL) perror("Error creating file!");

  // header of the .csv
  fprintf(csv, "size,trial,result\n" );
  // measurement results
  for(int i = 0; i < results_c; i++){
    for (int j = 0; j < trials; j++) {
      fprintf(csv, "%d,%d,%lf\n", i * stride, j, results[i][j]);
    }
  }
  fclose(csv);

}

void init_vector(int max_size, int stride){
  int actual_size = max_size * stride;
  default_vector = (double*) malloc(actual_size * sizeof(double));
  for (int i = 0; i < actual_size; i++) {
    default_vector[i] = i;
  }
}

void init_matrix(int max_size, int stride){
  int actual_size = max_size * stride;
  default_matrix = (double*) malloc(actual_size* actual_size * sizeof(double));
  for (int i = 0; i < actual_size; i++) {
    for (int j = 0; j < actual_size; j++) {
      printf("%d,%d\n",i,j );
      default_matrix[actual_size * i + j] = actual_size * i + j;
    }
  }
}

void init_results(int max_size, int trials){
  results = (double**) malloc(max_size * sizeof(double*));
  for (int i = 0; i < max_size; i++) {
    results[i] = (double*) malloc(trials * sizeof(double));
  }
}

void cleanup_vector(int max_size, int stride){
  free(default_vector);
}

void cleanup_matrix(int max_size, int stride){

  int actual_size = max_size * stride;
  free(default_matrix);
}

void cleanup_results(int max_size){
  for (int i = 0; i < max_size; i++) {
    free(results[i]);
  }
  free(results);
}

double vector_multi_test(int size){
  gsl_vector_view tmp = gsl_vector_view_array(default_vector, size);
  double res;
  double start = get_real_time();
  gsl_blas_ddot(&tmp.vector, &tmp.vector, &res);
  double end = get_real_time();
  return end - start;
}

double matrix_multi_test(int size){
  gsl_vector_view tmp_v = gsl_vector_view_array(default_vector, size);
  gsl_matrix_view tmp_m = gsl_matrix_view_array(default_matrix, size, size);
  double res;
  double start = get_real_time();
  gsl_blas_dgemv(CblasNoTrans, 1, &tmp_m.matrix, &tmp_v.vector, 0, &tmp_v.vector);

  double end = get_real_time();
  return end - start;

}

void vector_analysis(int max_size, int stride, int trials){
  init_vector(max_size, stride);

  for(int i = 1; i < max_size; i++){
    for (int j = 0; j < trials; j++) {
      double res = vector_multi_test(i * stride);
      results[i][j] = res;
    }
  }

  save_res_to_csv(results, max_size, stride, trials, "vector_multi.csv");
  cleanup_vector(max_size, stride);
}

void matrix_analysis(int max_size, int stride, int trials){
  printf("v\n");
  init_vector(max_size, stride);
  printf("m\n");
  init_matrix(max_size, stride);
  printf("a\n" );
  for(int i = 1; i < max_size; i++){
    for (int j = 0; j < trials; j++) {
      printf("%d,%d\n",i,j );
      double res = matrix_multi_test(i * stride);
      results[i][j] = res;
    }
  }
  printf("b\n" );
  save_res_to_csv(results, max_size, stride, trials, "matrix_multi.csv");
  printf("cleanup_matrix\n");
  cleanup_matrix(max_size, stride);
  printf("cleanup_vector\n");

  cleanup_vector(max_size, stride);
  printf("c\n");

}

int
main (void)
{
  int vector_max_size = 1000;
  int vector_stride = 1000;
  int matrix_max_size = 100; //50
  int matrix_stride = 100;
  int trials = 10;
  int results_size = max(vector_max_size, matrix_max_size);
  init_results(results_size, trials);
  printf("Analyzing vector*vector multiplication.\n" );
  vector_analysis(vector_max_size, vector_stride, trials);
  printf("Analyzing matrix*vector multiplication.\n" );
  // matrix_analysis(matrix_max_size, matrix_stride, trials);
  cleanup_results(results_size);

  return 0;
}


text:
#pragma once

#include "icode_generator.h"
#include <halley/data_structures/hash_map.h>
#include <gsl/gsl>
#include "halley/file/path.h"

namespace YAML
{
	class Node;
}

namespace Halley
{
	class ECSData;
	class ComponentSchema;
	class SystemSchema;
	class MessageSchema;
	class CustomTypeSchema;

	class Codegen
	{
		struct Stats
		{
			int written = 0;
			int skipped = 0;
			std::vector<Path> files;
		};

	public:
		using ProgressReporter = std::function<bool(float, String)>;

		static void run(Path inDir, Path outDir);
		static std::vector<Path> generateCode(const ECSData& data, Path directory);

	private:
		static bool writeFile(Path path, const char* data, size_t dataSize, bool stub);
		static void writeFiles(Path directory, const CodeGenResult& files, Stats& stats);
	};
}

max_stars_repo_name: code-disaster/halley

max_issues_repo_name: code-disaster/halley

max_forks_repo_name: code-disaster/halley