text:
/*
 * author: Achim Gaedke
 * created: January 2002
 * file: pygsl/src/multiminmodule.c
 * $Id: multimin.c,v 1.6 2008/10/25 20:45:04 schnizer Exp $
 */
#include <pygsl/solver.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_multimin.h>
#include "multiminmodule_doc.h"

const char  * filename = __FILE__;
PyObject *module = NULL;

/*
 * I have two different types to minimize and I am too lazy to implement the same
 * twice. Therefore I use unions and add a flag to the minimizer for the type
 */

/* The Callbacks */
double
PyGSL_multimin_function_f(const gsl_vector* x, void* params)
{
     double result;
     int flag;
     int i;

     FUNC_MESS_BEGIN();
     PyGSL_solver *min_o;
     min_o = (PyGSL_solver *) params;
     assert(PyGSL_solver_check(min_o));
     for(i = 0; i<x->size; i++){
	  DEBUG_MESS(2, "Got a x[%d] of %f", i, gsl_vector_get(x, i));
     }
     assert(min_o->mstatic->n_cbs >= 1);
     flag = PyGSL_function_wrap_On_O(x, min_o->cbs[0], min_o->args, &result,
				     NULL, x->size, __FUNCTION__);
     if (flag!= GSL_SUCCESS){
	  result = gsl_nan();
	  if(min_o->isset == 1) longjmp(min_o->buffer,flag);
     }
     DEBUG_MESS(2, "Got a result of %f", result);
     FUNC_MESS_END();
     return result;
}

void
PyGSL_multimin_function_df(const gsl_vector* x, void* params, gsl_vector *df)
{
     int flag, i;
     PyGSL_solver *min_o;
     min_o = (PyGSL_solver *) params;

     FUNC_MESS_BEGIN();
     assert(PyGSL_solver_check(min_o));
     for(i = 0; i<x->size; i++){
	  DEBUG_MESS(2, "Got a x[%d] of %f", i, gsl_vector_get(x, i));
     }
     assert(min_o->mstatic->n_cbs >= 2);
     flag = PyGSL_function_wrap_Op_On(x, df, min_o->cbs[1], min_o->args,
				      x->size, x->size, __FUNCTION__);
     for(i = 0; i<df->size; i++){
	  DEBUG_MESS(2, "Got df x[%d] of %f", i, gsl_vector_get(df, i));
     }
     if(flag!=GSL_SUCCESS){
	  if(min_o->isset == 1) longjmp(min_o->buffer,flag);
     }
     FUNC_MESS_END();
}

void
PyGSL_multimin_function_fdf(const gsl_vector* x, void* params, double *f, gsl_vector *df)
{
     int flag, i;
     PyGSL_solver *min_o;

     FUNC_MESS_BEGIN();
     min_o = (PyGSL_solver *) params;
     assert(PyGSL_solver_check(min_o));
     for(i = 0; i<x->size; i++){
	  DEBUG_MESS(2, "Got a x[%d] of %f", i, gsl_vector_get(x, i));
     }
     assert(min_o->mstatic->n_cbs >= 3);
     flag = PyGSL_function_wrap_On_O(x, min_o->cbs[2], min_o->args, f,
				     df, x->size, __FUNCTION__);
     DEBUG_MESS(2, "Got a result of %f", *f);
     for(i = 0; i<df->size; i++){
	  DEBUG_MESS(2, "Got df x[%d] of %f", i, gsl_vector_get(df, i));
     }
     if (flag!= GSL_SUCCESS){
	  *f = gsl_nan();
	  if(min_o->isset == 1) longjmp(min_o->buffer,flag);
     }
     FUNC_MESS_END();
     return;
}

static PyObject*
PyGSL_multimin_set_f(PyGSL_solver *self, PyObject *pyargs, PyObject *kw)
{

     PyObject* func = NULL, * args = Py_None, * x = NULL, * steps = NULL;
     PyArrayObject * xa = NULL, * stepsa = NULL;
     int n=0, flag=GSL_EFAILED;
     PyGSL_array_index_t stride_recalc;
     gsl_vector_view gsl_x;
     gsl_vector_view gsl_steps;
     gsl_multimin_function * c_sys = NULL;
     static const char *kwlist[] = {"f", "x0", "args", "steps", NULL};

     FUNC_MESS_BEGIN();
     assert(PyGSL_solver_check(self));
     if (self->solver == NULL) {
	  pygsl_error("Got a NULL Pointer of min.f",  filename, __LINE__ - 3, GSL_EFAULT);
	  return NULL;
     }

     assert(pyargs);
     /* arguments PyFunction, Parameters, start Vector, step Vector */
     if (0==PyArg_ParseTupleAndKeywords(pyargs,kw,"OOOO", (char **)kwlist, &func,
					&x,&args,&steps))
	  return NULL;

     if(!PyCallable_Check(func)){
	  pygsl_error("First argument must be callable",  filename, __LINE__ - 3, GSL_EBADFUNC);
	  return NULL;
     }
     n=self->problem_dimensions[0];
     xa  = PyGSL_vector_check(x, n, PyGSL_DARRAY_INPUT(4), &stride_recalc, NULL);
     if (xa == NULL){
	  PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__ - 1);
	  goto fail;
     }

     gsl_x = gsl_vector_view_array_with_stride((double *)(xa->data), stride_recalc, xa->dimensions[0]);
     stepsa =  PyGSL_vector_check(steps, n, PyGSL_DARRAY_INPUT(5), &stride_recalc, NULL);
     if (stepsa == NULL){
	  PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__ - 1);
	  goto fail;
     }
     gsl_steps = gsl_vector_view_array_with_stride((double *)(stepsa->data), stride_recalc, stepsa->dimensions[0]);

     if (self->c_sys != NULL) {
	  /* free the previous function and args */
	  c_sys = self->c_sys;
     } else {
	  /* allocate function space */
	  c_sys=calloc(1, sizeof(gsl_multimin_function));
	  if (c_sys == NULL) {
	       pygsl_error("Could not allocate the object for the minimizer function",
			 filename, __LINE__ - 3, GSL_ENOMEM);
	       goto fail;
	  }
     }
     DEBUG_MESS(3, "Everything allocated args = %p", args);
     /* add new function  and parameters */

     if(PyGSL_solver_func_set(self, args, func, NULL, NULL) != GSL_SUCCESS)
	  goto fail;

     /* initialize the function struct */
     c_sys->n=n;
     c_sys->f=PyGSL_multimin_function_f;
     c_sys->params=(void*)self;

     DEBUG_MESS(3, "Setting jmp buffer isset = % d", self->isset);
     if((flag = setjmp(self->buffer)) == 0){
	  self->isset = 1;
	  flag = gsl_multimin_fminimizer_set(self->solver,c_sys, &gsl_x.vector, &gsl_steps.vector);
	  if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS){
	       goto fail;
	  }
     } else {
	  goto fail;
     }
     DEBUG_MESS(4, "Set evaluated. flag = %d", flag);
     self->c_sys = c_sys;
     Py_DECREF(xa);
     Py_DECREF(stepsa);
     self->set_called = 1;

     Py_INCREF(Py_None);
     self->isset = 0;
     FUNC_MESS_END();
     return Py_None;

 fail:
     FUNC_MESS("Fail");
     PyGSL_ERROR_FLAG(flag);
     self->isset = 0;
     Py_XDECREF(xa);
     Py_XDECREF(stepsa);
     return NULL;

}

static PyObject*
PyGSL_multimin_set_fdf(PyGSL_solver *self, PyObject *pyargs, PyObject *kw)
{

     PyObject * f = NULL, * df = NULL, * fdf = NULL, * args = Py_None,
	      * x = NULL;
     PyArrayObject * xa = NULL;
     int n=0, flag=GSL_EFAILED;
     PyGSL_array_index_t stride_recalc=-1;
     double step=0.01, tol=1e-4;
     gsl_vector_view gsl_x;
     gsl_multimin_function_fdf * c_sys;
     static const char *kwlist[] = {"f", "df", "fdf",  "x0", "args", "step", "tol", NULL};

     FUNC_MESS_BEGIN();
     assert(PyGSL_solver_check(self));
     if (self->solver == NULL) {
	  pygsl_error("Got a NULL Pointer of min.fdf", filename, __LINE__ - 3, GSL_EFAULT);
	       return NULL;
	  }

     /* arguments PyFunction, Parameters, start Vector, step Vector */
     if (0==PyArg_ParseTupleAndKeywords(pyargs, kw, "OOOO|Odd", (char **)kwlist,
					&f, &df, &fdf, &x, &args, &step, &tol))
	  return NULL;

     n=self->problem_dimensions[0];
     xa  = PyGSL_vector_check(x, n, PyGSL_DARRAY_INPUT(4), &stride_recalc, NULL);
     if (xa == NULL){
	  PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__ - 1);
	  goto fail;
     }
     gsl_x = gsl_vector_view_array_with_stride((double *)(xa->data), stride_recalc, xa->dimensions[0]);

     if (self->c_sys != NULL) {
	  /* free the previous function and args */
	  c_sys = self->c_sys;
     } else {
	  /* allocate function space */
	  c_sys=malloc(sizeof(gsl_multimin_function_fdf));
	  if (c_sys==NULL) {
	       pygsl_error("Could not allocate the object for the minimizer function",
			 filename, __LINE__ - 3, GSL_ENOMEM);
	       goto fail;
	  }
     }
     if(PyGSL_solver_func_set(self, args, f, df, fdf) != GSL_SUCCESS)
	  goto fail;

     /* initialize the function struct */
     c_sys->n=n;
     c_sys->f  =PyGSL_multimin_function_f;
     c_sys->df =PyGSL_multimin_function_df;
     c_sys->fdf=PyGSL_multimin_function_fdf;
     c_sys->params=(void*)self;

     if((flag = setjmp(self->buffer)) == 0){
	  self->isset = 1;
	  flag = gsl_multimin_fdfminimizer_set(self->solver, c_sys, &gsl_x.vector, step, tol);
	  if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS){
	       goto fail;
	  }
     }else{
	  goto fail;
     }
     self->c_sys = c_sys;
     self->isset = 0;
     Py_DECREF(xa);
     self->set_called = 1;

     Py_INCREF(Py_None);
     FUNC_MESS_END();
     return Py_None;

 fail:
     PyGSL_ERROR_FLAG(flag);
     self->isset = 0;
     Py_XDECREF(xa);
     return NULL;

}

static PyObject*
PyGSL_multimin_dx(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_vec(self, args, (ret_vec)gsl_multimin_fdfminimizer_dx);
}

static PyObject*
PyGSL_multimin_gradient(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_vec(self, args, (ret_vec)gsl_multimin_fdfminimizer_gradient);
}

static PyObject*
PyGSL_multimin_f_x(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_vec(self, args, (ret_vec)gsl_multimin_fminimizer_x);
}

static PyObject*
PyGSL_multimin_fdf_x(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_vec(self, args, (ret_vec)gsl_multimin_fdfminimizer_x);
}

static PyObject*
PyGSL_multimin_f_minimum(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_multimin_fminimizer_minimum);
}

static PyObject*
PyGSL_multimin_fdf_minimum(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_multimin_fdfminimizer_minimum);
}

static PyObject*
PyGSL_multimin_f_size(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_multimin_fminimizer_size);
}

static PyObject*
PyGSL_multimin_test_size_method(PyGSL_solver *self, PyObject *args)
{
     int flag=GSL_EFAILED;
     double epsabs;

     FUNC_MESS_BEGIN();

     if (0==PyArg_ParseTuple(args,"d", &epsabs))
	  return NULL;
     flag = gsl_multimin_test_size(gsl_multimin_fminimizer_size(self->solver), epsabs);
     FUNC_MESS_END();
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);
}

static PyObject*
PyGSL_multimin_test_gradient_method(PyGSL_solver * self, PyObject *args)
{

     double epsabs;
     int flag;
     FUNC_MESS_BEGIN();

     assert(PyGSL_solver_check(self));
     if (0==PyArg_ParseTuple(args,"d", &epsabs))
	  return NULL;

     flag = gsl_multimin_test_gradient(gsl_multimin_fdfminimizer_gradient(self->solver), epsabs);
     FUNC_MESS_END();
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);

}

static PyMethodDef PyGSL_multimin_fmethods[] = {
     {"x",           (PyCFunction)PyGSL_multimin_f_x,      METH_NOARGS,(char *)multimin_x_doc,      },
     {"minimum",     (PyCFunction)PyGSL_multimin_f_minimum,METH_NOARGS,(char *)multimin_minimum_doc,},
     {"set",      (PyCFunction)PyGSL_multimin_set_f,           METH_VARARGS|METH_KEYWORDS, (char *)multimin_set_f_doc    },
     {"size",     (PyCFunction)PyGSL_multimin_f_size,            METH_NOARGS,  (char *)multimin_size_doc     },
     {"test_size",(PyCFunction)PyGSL_multimin_test_size_method,METH_VARARGS, (char *)multimin_test_size_doc},
     {NULL, NULL, 0, NULL}           /* sentinel */
};

static PyMethodDef PyGSL_multimin_fdfmethods[] = {
     {"x",           (PyCFunction)PyGSL_multimin_fdf_x,      METH_NOARGS,(char *)multimin_x_doc,      },
     {"minimum",     (PyCFunction)PyGSL_multimin_fdf_minimum,METH_NOARGS,(char *)multimin_minimum_doc,},
     {"set",      (PyCFunction)PyGSL_multimin_set_fdf,                 METH_VARARGS|METH_KEYWORDS, (char *)multimin_set_fdf_doc      },
     {"dx",       (PyCFunction)PyGSL_multimin_dx,                      METH_NOARGS,  (char *)multimin_dx_doc           },
     {"gradient", (PyCFunction)PyGSL_multimin_gradient,                METH_NOARGS,  (char *)multimin_gradient_doc     },
     {"test_gradient",(PyCFunction)PyGSL_multimin_test_gradient_method,METH_VARARGS, (char *)multimin_test_gradient_doc},
     {NULL, NULL, 0, NULL}           /* sentinel */
};

const struct _GSLMethods
multimin_f   = { (void_m_t) gsl_multimin_fminimizer_free,
		/* gsl_multimin_fminimizer_restart */  (void_m_t) NULL,
		  (name_m_t) gsl_multimin_fminimizer_name,
		 (int_m_t) gsl_multimin_fminimizer_iterate};

const struct _GSLMethods
multimin_fdf = {(void_m_t) gsl_multimin_fdfminimizer_free,
		(void_m_t) gsl_multimin_fdfminimizer_restart,
		(name_m_t) gsl_multimin_fdfminimizer_name,
		(int_m_t)  gsl_multimin_fdfminimizer_iterate};

static const char multimin_f_type_name[] = "F-MultiMinimizer";
static const char multimin_fdf_type_name[] = "FdF-MultiMinimizer";

const struct _SolverStatic
multimin_solver_f   = {{(void_m_t) gsl_multimin_fminimizer_free,
			(void_m_t) NULL,
			(name_m_t) gsl_multimin_fminimizer_name,
			(int_m_t) gsl_multimin_fminimizer_iterate},
		       1, PyGSL_multimin_fmethods, multimin_f_type_name},
multimin_solver_fdf = {{(void_m_t) gsl_multimin_fdfminimizer_free,
			(void_m_t) gsl_multimin_fdfminimizer_restart,
			(name_m_t) gsl_multimin_fdfminimizer_name,
			(int_m_t)  gsl_multimin_fdfminimizer_iterate},
		       3, PyGSL_multimin_fdfmethods, multimin_fdf_type_name};

static PyObject*
PyGSL_multimin_f_init(PyObject *self, PyObject *args,
		      const gsl_multimin_fminimizer_type * type)
{

     PyObject *tmp=NULL;
     solver_alloc_struct s = {type, (void_an_t) gsl_multimin_fminimizer_alloc,
			      &multimin_solver_f};
     FUNC_MESS_BEGIN();
     tmp = PyGSL_solver_dn_init(self, args, &s, 1);
     FUNC_MESS_END();
     return tmp;
}

static PyObject*
PyGSL_multimin_fdf_init(PyObject *self, PyObject *args,
		      const gsl_multimin_fdfminimizer_type * type)
{

     PyObject *tmp=NULL;
     solver_alloc_struct s = {type, (void_an_t) gsl_multimin_fdfminimizer_alloc,
			      &multimin_solver_fdf};
     FUNC_MESS_BEGIN();
     tmp = PyGSL_solver_dn_init(self, args, &s, 1);
     FUNC_MESS_END();
     return tmp;
}

#define AMINIMIZER_F(name)                                                  \
static PyObject* PyGSL_multimin_init_ ## name (PyObject *self, PyObject *args)\
{                                                                             \
     PyObject *tmp = NULL;                                                    \
     FUNC_MESS_BEGIN();                                                       \
     tmp = PyGSL_multimin_f_init(self, args,  gsl_multimin_fminimizer_ ## name); \
     if (tmp == NULL){                                                        \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__); \
     }                                                                        \
     FUNC_MESS_END();                                                         \
     return tmp;                                                              \
}

#define AMINIMIZER_FDF(name)                                                  \
static PyObject* PyGSL_multimin_init_ ## name (PyObject *self, PyObject *args)\
{                                                                             \
     PyObject *tmp = NULL;                                                    \
     FUNC_MESS_BEGIN();                                                       \
     tmp = PyGSL_multimin_fdf_init(self, args,  gsl_multimin_fdfminimizer_ ## name); \
     if (tmp == NULL){                                                        \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__); \
     }                                                                        \
     FUNC_MESS_END();                                                         \
     return tmp;                                                              \
}

AMINIMIZER_F(nmsimplex)
AMINIMIZER_FDF(steepest_descent)
AMINIMIZER_FDF(vector_bfgs)
AMINIMIZER_FDF(conjugate_pr)
AMINIMIZER_FDF(conjugate_fr)

static PyObject*
PyGSL_multimin_test_size(PyObject * self, PyObject *args)
{
     double size, epsabs;
     int flag = GSL_EFAILED;
     FUNC_MESS_BEGIN();
     if (0==PyArg_ParseTuple(args,"dd", &size, &epsabs))
	  return NULL;
     flag = gsl_multimin_test_size(size, epsabs);
     FUNC_MESS_END();
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);

}

static PyObject*
PyGSL_multimin_test_gradient(PyObject * self, PyObject *args)
{
     return PyGSL_solver_vd_i (self, args, gsl_multimin_test_gradient);
}

static PyMethodDef mMethods[] = {
     /* solvers */
     {"nmsimplex",        PyGSL_multimin_init_nmsimplex,        METH_VARARGS, (char *)nmsimplex_doc       },
     {"steepest_descent", PyGSL_multimin_init_steepest_descent, METH_VARARGS, (char *)steepest_descent_doc},
     {"vector_bfgs",      PyGSL_multimin_init_vector_bfgs,      METH_VARARGS, (char *)vector_bfgs_doc     },
     {"conjugate_pr",     PyGSL_multimin_init_conjugate_pr,     METH_VARARGS, (char *)conjugate_pr_doc    },
     {"conjugate_fr",     PyGSL_multimin_init_conjugate_fr,     METH_VARARGS, (char *)conjugate_fr_doc    },
     /* multimin funcs */
     {"test_size",        PyGSL_multimin_test_size,             METH_VARARGS, (char *)test_size_doc       },
     {"test_gradient",    PyGSL_multimin_test_gradient,         METH_VARARGS, (char *)test_gradient_doc   },
     {NULL, NULL, 0, NULL}
};

void
initmultimin(void)
{
     PyObject* m, *dict, *item;
     FUNC_MESS_BEGIN();

     m=Py_InitModule("multimin", mMethods);
     module = m;
     assert(m);
     dict = PyModule_GetDict(m);
     if(!dict)
	  goto fail;

     init_pygsl()
     import_pygsl_solver();
     assert(PyGSL_API);

     if (!(item = PyString_FromString((char*)PyGSL_multimin_module_doc))){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not generate module doc string!");
	  goto fail;
     }
     if (PyDict_SetItemString(dict, "__doc__", item) != 0){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not init doc string!");
	  goto fail;
     }

     FUNC_MESS_END();

 fail:
     FUNC_MESS("FAIL");
     return;
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */
/*
 *
 *
 *       Filename:  write.c
 *
 *    Description:  Write the data into file .dat
 *
 *        Version:  1.0
 *        Created:  07/05/2014 22:46:50
 *       Revision:  none
 *        License:  BSD
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 *
 */

#include "funcs.h"
#include <stdio.h>
#include <gsl/gsl_matrix.h>

/*
 *      FUNCTION
 *         Name:  save_integrals
 *  Description:
 *
 */
int save_integrals ( void* params )
{
	/* perform the integrals */
  	double integrals[12] ;
	void* p = (void*) params ;
	int status = integration ( p, integrals ) ;

	/* write them into a file */
	FILE* f_integ = fopen ( "INTEGRALS.dat", "w+" ) ;
	int i ;
	for ( i = 0 ; i < 12 ; i++ )
		fprintf( f_integ, "integrals[%d]: %.9f\n", i, integrals[i] ) ;
	fclose( f_integ ) ;

	return status;
}		/* -----  end of function save_integrals  ----- */

/*
 *      FUNCTION
 *         Name:  save_matrices
 *  Description:
 *
 */
int save_matrices ( void* params )
{
	void* p = (void *) params ;

	/* create the Redfield matrix and save it into a file */
	gsl_matrix* red_matrix = gsl_matrix_calloc ( 4, 4 ) ;
	int status1 = red_mat ( red_matrix, p ) ;
	int status2 = mat_write ( red_matrix, "REDFIELD_MATRIX" ) ;

	/* create the CP matrix and save it into a file */
	gsl_matrix* cp_matrix = gsl_matrix_calloc ( 4, 4 ) ;
	int status3 = cp_mat ( cp_matrix, p ) ;
	int status4 = mat_write ( cp_matrix, "CP_MATRIX" ) ;

	gsl_matrix_free (red_matrix) ;
	gsl_matrix_free (cp_matrix) ;

	return status1+status2+status3+status4 ;
}		/* -----  end of function save_matrices  ----- */

ext: c


text:
#include <petsc.h>

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

    PetscErrorCode ierr;
    PetscMPIInt    rank;
    PetscInt       i;
    PetscReal      localval, globalsum;
    PetscReal      x = 1; // default value

    ierr = PetscInitialize(&argc, &argv, NULL,
                           "Compute e in parallel with PETSc.\n\n");
    if (ierr) return ierr;

    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);

    ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL,
                             "Maclaurin Series for exp(x)",
                             NULL); CHKERRQ(ierr);
      ierr = PetscOptionsReal("-x", "input to exp(x) function",
                              NULL, x, &x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    // compute  x^n / n!  where n = (rank of process) + 1
    localval = 1.;
    for (i = 1; i < rank+1; i++) localval *= x/i;

    // sum the contributions over all processes
    ierr = MPI_Allreduce(&localval, &globalsum, 1, MPIU_REAL, MPIU_SUM,
                         PETSC_COMM_WORLD); CHKERRQ(ierr);

    // output estimate of e and report on work from each process
    ierr = PetscPrintf(PETSC_COMM_WORLD,
                       "exp(x) is about %17.15f\n",
                       globalsum); CHKERRQ(ierr);

    ierr = PetscPrintf(PETSC_COMM_SELF,
                       "rank %d did %d flops\n",
                       rank, (rank > 0) ? rank-1 : 0); CHKERRQ(ierr);

    return PetscFinalize();
}


text:
/*
** gsl routines for double data type
**
** G.Lohmann, July 2004
*/

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

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

#define dvset gsl_vector_set
#define dvget gsl_vector_get
#define dmset gsl_matrix_set
#define dmget gsl_matrix_get

#define fvset gsl_vector_float_set
#define fvget gsl_vector_float_get
#define fmset gsl_matrix_float_set
#define fmget gsl_matrix_float_get

#define ABS(x) ((x) > 0 ? (x) : -(x))

/*
**  y = Ax
*/
gsl_vector *
dmat_x_vector(gsl_matrix *A, gsl_vector *x, gsl_vector *y)
{
  int i,j,nrows,ncols;
  double *ptr1,*ptr2,*ptr3,sum;

  nrows = A->size1;
  ncols = A->size2;

  if (y == NULL) {
     y = gsl_vector_alloc (nrows);
  }

  if (x->size != ncols || y->size != nrows) {
    fprintf(stderr," dmat_x_vector: incongruent dimensions\n");
    exit(0);
  }

  ptr1 = A->data;
  ptr3 = y->data;

  for (i=0; i<nrows; i++) {
    sum = 0;
    ptr2 = x->data;
    for (j=0; j<ncols; j++) {
      sum += (*ptr1++) * (*ptr2++);
    }
    *ptr3++ = sum;
  }
  return y;
}

/*
**  y = x^T A
*/
gsl_vector *
dvector_x_mat(gsl_vector *x, gsl_matrix *A, gsl_vector *y)
{
  int i,j,nrows,ncols;
  double *ptr1,*ptr2,*ptr3,sum;

  nrows = A->size1;
  ncols = A->size2;

  if (y == NULL) {
     y = gsl_vector_alloc (nrows);
  }

  if (x->size != ncols || y->size != nrows) {
    fprintf(stderr," dmat_x_vector: incongruent dimensions\n");
    exit(0);
  }

  ptr1 = A->data;
  ptr3 = y->data;

  for (i=0; i<nrows; i++) {
    sum = 0;
    ptr2 = x->data;
    for (j=0; j<ncols; j++) {
      sum += (*ptr1++) * (*ptr2++);
    }
    *ptr3++ = sum;
  }
  return y;
}

/*
**  z = x^T y
*/
double
dskalarproduct(gsl_vector *x,gsl_vector *y)
{
  int i,n;
  double *ptr1,*ptr2,sum;

  n = x->size;
  if (y->size != n) {
    fprintf(stderr," dskalarproduct: incongruent vector sizes: %d %d",n,y->size);
    exit(0);
  }

  ptr1 = x->data;
  ptr2 = y->data;
  sum = 0;
  for (i=0; i<n; i++) {
    sum += (*ptr1) * (*ptr2);
    ptr1++;
    ptr2++;
  }
  return sum;
}

/*
**    B = A^-1 = A^+
*/
gsl_matrix *
dmat_PseudoInv(gsl_matrix *A,gsl_matrix *B)
{
  int j,k,l,n,m;
  double u,x;
  static gsl_matrix *U=NULL,*V=NULL;
  static gsl_vector *w=NULL;

  m = A->size1;
  n = A->size2;

  if (B == NULL) {
    B = gsl_matrix_alloc (n,m);
  }
  else if (B->size1 != n || B->size2 != m) {
    gsl_matrix_free(B);
    B = gsl_matrix_alloc (n, m);
  }

  if (U == NULL) {
    U = gsl_matrix_alloc (m, n);
    V = gsl_matrix_alloc (n, n);
    w = gsl_vector_alloc (n);
  }
  else if (U->size1 != m || w->size != n) {
    gsl_matrix_free(U);
    gsl_matrix_free(V);
    gsl_vector_free(w);
    U = gsl_matrix_alloc (m, n);
    V = gsl_matrix_alloc (n, n);
    w = gsl_vector_alloc (n);
  }

  /* singular value decomposition */
  gsl_matrix_memcpy(U,A);
  gsl_linalg_SV_decomp_jacobi(U,V,w);

  gsl_matrix_set_zero(B);

  for (k=0; k<n; k++) {
    for (l=0; l<m; l++) {
      for (j=0; j<n; j++) {
	u = dvget(w,j);
        if (ABS(u) > 1.0e-6) {
          x = dmget(B,k,l);
	  x += dmget(V,k,j)*dmget(U,l,j)/u;
          dmset(B,k,l,x);
        }
      }
    }
  }
  return B;
}

/*
**    C = A x B^T
*/
gsl_matrix *
dmat_x_matT(gsl_matrix *A,gsl_matrix *B,gsl_matrix *C)
{
  int i,j,k;
  int n,r,m;
  double *ptr1,*ptr2,*ptr3,sum;

  n = A->size1;
  r = A->size2;
  m = B->size1;

  if (B->size2 != r) {
    fprintf(stderr,"fmat_x_matT: incongruent matrix dimensions (A,B).\n");
    exit(0);
  }

  if (C == NULL) {
    C = gsl_matrix_alloc (n,m);
  }
  else {
    if (C->size1 != n || C->size2 != m) {
      fprintf(stderr,"fmat_x_matT: incongruent matrix dimensions(C).\n");
      exit(0);
    }
  }
  ptr1 = C->data;
  for (i=0; i<n; i++) {
    for (j=0; j<m; j++) {

      ptr2 = gsl_matrix_ptr (A,i,0);
      ptr3 = gsl_matrix_ptr (B,j,0);

      sum = 0;
      for (k=0; k<r; k++) {
	sum  += (*ptr2) * (*ptr3);
	ptr2++;
	ptr3++;
      }
      *ptr1++ = sum;
    }
  }
  return C;
}

/*
**    C = A^T x B
*/
gsl_matrix *
dmatT_x_mat(gsl_matrix *A,gsl_matrix *B,gsl_matrix *C)
{
  int i,j,k;
  int n,m,r;
  double *ptr1,*ptr2,*ptr3,sum;

  n = A->size1;
  r = A->size2;
  m = B->size2;

  if (B->size1 != n) {
    fprintf(stderr,"dmatT_x_mat: incongruent matrix dimensions (A,B).\n");
    exit(0);
  }

  if (C == NULL) {
    C = gsl_matrix_alloc (r,m);
  }
  else {
    if (C->size1 != r || C->size2 != m) {
      fprintf(stderr,"dmatT_x_mat: incongruent matrix dimensions (C, %d %d, %d %d).\n",
	      C->size1,C->size2,r,m);
      exit(0);
    }
  }

  ptr1 = C->data;
  for (i=0; i<r; i++) {
    for (j=0; j<m; j++) {

      ptr2 = gsl_matrix_ptr (A,0,i);
      ptr3 = gsl_matrix_ptr (B,0,j);

      sum = 0;
      for (k=0; k<n; k++) {
	sum  += (*ptr2) * (*ptr3);
	ptr2 += A->tda;
	ptr3 += B->tda;
      }
      *ptr1++ = sum;
    }
  }

  return C;
}

/*
**    C = A x B
*/
gsl_matrix *
dmat_x_mat(gsl_matrix *A,gsl_matrix *B,gsl_matrix *C)
{
  int i,j,k,m,n,r;
  int nrowsA,ncolsA;
  int nrowsB,ncolsB;
  int nrowsC,ncolsC;
  double *ptr1,*ptr2,*ptr3,sum;

  nrowsA = A->size1;
  ncolsA = A->size2;

  nrowsB = B->size1;
  ncolsB = B->size2;

  if (C == NULL)
    C = gsl_matrix_alloc (nrowsA,ncolsB);

  nrowsC = C->size1;
  ncolsC = C->size2;

  if (ncolsA != nrowsB || nrowsA != nrowsC || ncolsB != ncolsC) {
    fprintf(stderr,"dmat_x_mat: incongruent matrix dimensions.\n");
    exit(0);
  }

  m = A->size1;
  n = A->size2;
  r = B->size2;

  ptr1 = C->data;

  for (i=0; i<m; i++) {
    for (j=0; j<r; j++) {

      ptr2 = gsl_matrix_ptr (A,i,0);
      ptr3 = gsl_matrix_ptr (B,0,j);

      sum = 0;
      for (k=0; k<n; k++) {
	sum += (*ptr2++) * (*ptr3);
	ptr3 += B->tda;
      }
      *ptr1++ = sum;
    }
  }
  return C;
}

gsl_matrix *
dtranspose(gsl_matrix *A,gsl_matrix *B)
{
  int i,j,n,m;

  n = A->size1;
  m = A->size2;

  if (B == NULL) {
    B = gsl_matrix_alloc(m,n);
  }
  else if (B->size1 != m || B->size2 != n) {
    gsl_matrix_free(B);
    B = gsl_matrix_alloc(m,n);
  }

  for (i=0; i<n; i++) {
    for (j=0; j<m; j++) {
      dmset(B,j,i,dmget(A,i,j));
    }
  }
  return B;
}

/*

int main (void)
{
  int i,j;
  int m,n,r;
  gsl_matrix *a=NULL,*b=NULL,*c=NULL;
  double x;

  n = 648;
  a = gsl_matrix_alloc(n,n);

  for (i=0; i<n; i++) {
    for (j=0; j<n; j++) {
      x = i+j+1;
      dmset(a,i,j,x);
    }
  }

  b = gsl_matrix_alloc(n,n);
  for (i=0; i<n; i++) {
    for (j=0; j<n; j++) {
      x = i-j+1;
      dmset(b,i,j,x);
    }
  }

  c = dmat_x_mat(a,b,NULL);
  exit(0);
}
*/

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/*
    hyper.c
*/
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_sf_psi.h>

static double update_hyper(double param, int *vec, int **mat, int nrow, int ncol){
    double numer = 0.0;
    double denom = 0.0;
    int i, j;

    for(i = 0;i < nrow;i++)
        for(j = 0;j < ncol;j++)
            numer += gsl_sf_psi((double)mat[i][j] + param);
    numer -= (double)nrow * (double)ncol * gsl_sf_psi(param);

    for(i = 0;i < nrow;i++)
        denom += (double)ncol * gsl_sf_psi((double)vec[i] + param * (double)ncol);
    denom -= (double)nrow * (double)ncol * gsl_sf_psi(param * (double)ncol);

    return param * numer / denom;
}

void update_alpha(double *alpha, int *n_m, int **n_mz, int ndocs, int nclass){
    double sum_alpha = 0.0;
    double numer = 0.0;
    double denom = 0.0;
    int d, k;
    for(k = 0;k < nclass;k++)
        sum_alpha += alpha[k];
    for(d = 0;d < ndocs;d++)
        denom += gsl_sf_psi((double)n_m[d] + sum_alpha);
    denom -= (double)ndocs * gsl_sf_psi(sum_alpha);

    for(k = 0;k < nclass;k++){
        numer = 0.0;
        for(d = 0;d < ndocs;d++)
            numer += gsl_sf_psi((double)n_mz[d][k] + alpha[k]);
        numer -= (double)ndocs * gsl_sf_psi(alpha[k]);
        alpha[k] = alpha[k] * numer / denom;
    }
}

double update_beta(double beta, int *n_z, int **n_zw, int nclass, int nlex){
    double new_beta = update_hyper(beta, n_z, n_zw, nclass, nlex);
    return new_beta;
}


text:
// build commmand:
// gcc -c test_fresnel.c -O2 -std=c99 -I../LISAsim -I/opt/local/include
// gcc -o test_fresnel test_fresnel.o struct.o ../LISAsim/LISAnoise.o -lgsl -lgslcblas

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

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "splinecoeffs.h"

/* Implementation of the Thomas algorithm to solve a tridiagonal system */
/* Note: assumes numerical stability (e.g. diagonal-dominated matrix) */
static void SolveTridiagThomas(
  gsl_vector* vectx,        /* Output: solution vector, length n - already allocated */
  gsl_vector* vecta,        /* Diagonal of the matrix, length n */
  gsl_vector* vectb,        /* Lower diagonal, length n-1 */
  gsl_vector* vectc,        /* Upper diagonal, length n-1 */
  gsl_vector* vecty,        /* Right-hand-side vector, length n */
  int n)                    /* Length of vectx, vecty */
{
  /* Check lengths */
  if(!(vectx->size==n && vecty->size==n && vecta->size==n && vectb->size==n-1 && vectc->size==n-1)) {
    printf("Error: incompatible lengths in SolveTridiagThomas.\n");
    exit(1);
  }

  /* Note: we modify the vectors c, y in place */
  double* chat = vectc->data;
  double* yhat = vecty->data;
  double* x = vectx->data;
  double* a = vecta->data;
  double* b = vectb->data;

  /* Sweep forward, computing the chat and yhat values */
  chat[0] = chat[0] / a[0];
  yhat[0] = yhat[0] / a[0];
  double factor;
  for(int i=1; i<=n-2; i++) {
    factor = 1./(a[i] - b[i-1]*chat[i-1]);
    chat[i] = chat[i] * factor;
    yhat[i] = (yhat[i] - b[i-1]*yhat[i-1]) * factor;
  }
  factor = 1./(a[n-1] - b[n-2]*chat[n-2]);
  yhat[n-1] = (yhat[n-1] - b[n-2]*yhat[n-2]) * factor;

  /* Solve for x going backward */
  x[n-1] = yhat[n-1];
  for(int i=n-2; i>=0; i--) {
    x[i] = yhat[i] - chat[i] * x[i+1];
  }
}

void BuildNotAKnotSpline(
  gsl_matrix* splinecoeffs,   /* Output: matrix containing all the spline coeffs (already allocated) */
  gsl_vector* vectx,          /* Input: vector x*/
  gsl_vector* vecty,          /* Input: vector y */
  int n)                      /* Size of x, y, and of output matrix */
{
  /* Check lengths */
  if(!(vectx->size==n && vecty->size==n && splinecoeffs->size1==n && splinecoeffs->size2==5)) {
    printf("Error: incompatible lengths in NotAKnotSpline.\n");
    exit(1);
  }
  double* x = vectx->data;
  double* y = vecty->data;

  /* Computing vecth and vectDeltay */
  gsl_vector* vecth = gsl_vector_alloc(n-1);
  gsl_vector* vectDeltay = gsl_vector_alloc(n-1);
  gsl_vector* vectDeltayoverh = gsl_vector_alloc(n-1);
  double* h = vecth->data;
  double* Deltay = vectDeltay->data;
  double* Deltayoverh = vectDeltayoverh->data;
  for(int i=0; i<n-1; i++) {
    h[i] = x[i+1] - x[i];
    Deltay[i] = y[i+1] - y[i];
    Deltayoverh[i] = Deltay[i] / h[i];
  }

  /* Structures for the tridiagonal system */
  gsl_vector* vectY = gsl_vector_alloc(n-2);
  gsl_vector* vecta = gsl_vector_alloc(n-2);
  gsl_vector* vectb = gsl_vector_alloc(n-3);
  gsl_vector* vectc = gsl_vector_alloc(n-3);
  double* Y = vectY->data;
  double* a = vecta->data;
  double* b = vectb->data;
  double* c = vectc->data;
  for(int i=0; i<=n-3; i++) {
    Y[i] = 3.*(Deltayoverh[i+1] - Deltayoverh[i]);
    a[i] = 2.*(h[i+1] + h[i]);
  }
  for(int i=0; i<=n-4; i++) {
    b[i] = h[i+1];
    c[i] = h[i+1];
  }
  /* Adjusting for the not-a-knot condition */
  a[0] += h[0] + h[0]*h[0]/h[1];
  c[0] += -h[0]*h[0]/h[1];
  a[n-3] += h[n-2] + h[n-2]*h[n-2]/h[n-3];
  b[n-4] += -h[n-2]*h[n-2]/h[n-3];

  /* Solving the tridiagonal system */
  gsl_vector* vectp2 = gsl_vector_alloc(n);
  gsl_vector_view viewp2trunc = gsl_vector_subvector(vectp2, 1, n-2);
  SolveTridiagThomas(&viewp2trunc.vector, vecta, vectb, vectc, vectY, n-2);
  double* p2 = vectp2->data;
  p2[0] = p2[1] - h[0]/h[1] * (p2[2] - p2[1]);
  p2[n-1] = p2[n-2] + h[n-2]/h[n-3] * (p2[n-2] - p2[n-3]);

  /* Deducing the p1's and the p3's */
  gsl_vector* vectp1 = gsl_vector_alloc(n);
  gsl_vector* vectp3 = gsl_vector_alloc(n);
  double* p1 = vectp1->data;
  double* p3 = vectp3->data;
  for(int i=0; i<=n-2; i++) {
    p1[i] = Deltayoverh[i] - h[i]/3. * (p2[i+1] + 2.*p2[i]);
    p3[i] = (p2[i+1] - p2[i]) / (3*h[i]);
  }
  /* Note: p1[n-1], p2[n-1], p3[n-1] are set to values coherent with the derivatives of the spline at the last point, but they are not stricly speaking coefficients of the spline. */
  p1[n-1] = p1[n-2] + 2.*p2[n-2]*h[n-2] + 3.*p3[n-2]*h[n-2]*h[n-2];
  p3[n-1] = p3[n-2];

  /* Copying the results in the output matrix */
  gsl_vector_view viewx = gsl_matrix_column(splinecoeffs, 0);
  gsl_vector_view viewp0 = gsl_matrix_column(splinecoeffs, 1);
  gsl_vector_view viewp1 = gsl_matrix_column(splinecoeffs, 2);
  gsl_vector_view viewp2 = gsl_matrix_column(splinecoeffs, 3);
  gsl_vector_view viewp3 = gsl_matrix_column(splinecoeffs, 4);
  gsl_vector_memcpy(&viewx.vector, vectx);
  gsl_vector_memcpy(&viewp0.vector, vecty);
  gsl_vector_memcpy(&viewp1.vector, vectp1);
  gsl_vector_memcpy(&viewp2.vector, vectp2);
  gsl_vector_memcpy(&viewp3.vector, vectp3);

  /* Cleanup*/
  gsl_vector_free(vecth);
  gsl_vector_free(vectDeltay);
  gsl_vector_free(vectDeltayoverh);
  gsl_vector_free(vectY);
  gsl_vector_free(vecta);
  gsl_vector_free(vectb);
  gsl_vector_free(vectc);
  gsl_vector_free(vectp1);
  gsl_vector_free(vectp2);
  gsl_vector_free(vectp3);
}

void BuildQuadSpline(
  gsl_matrix* splinecoeffs,   /* Output: matrix containing all the spline coeffs (already allocated) */
  gsl_vector* vectx,          /* Input: vector x*/
  gsl_vector* vecty,          /* Input: vector y */
  int n)                      /* Size of x, y, and of output matrix */
{
  /* Check lengths */
  if(!(vectx->size==n && vecty->size==n && splinecoeffs->size1==n && splinecoeffs->size2==4)) {
    printf("Error: incompatible lengths in NotAKnotSpline.\n");
    exit(1);
  }
  double* x = vectx->data;
  double* y = vecty->data;

  /* Computing vecth and vectDeltay */
  gsl_vector* vecth = gsl_vector_alloc(n-1);
  gsl_vector* vectDeltay = gsl_vector_alloc(n-1);
  gsl_vector* vectDeltayoverh = gsl_vector_alloc(n-1);
  double* h = vecth->data;
  double* Deltay = vectDeltay->data;
  double* Deltayoverh = vectDeltayoverh->data;
  for(int i=0; i<n-1; i++) {
    h[i] = x[i+1] - x[i];
    Deltay[i] = y[i+1] - y[i];
    Deltayoverh[i] = Deltay[i] / h[i];
  }

  /* Solving for p1 */
  gsl_vector* vectp1 = gsl_vector_alloc(n);
  double* p1 = vectp1->data;
  double ratio = h[n-2] / h[n-3];
  p1[n-3] = ((2. + ratio)*Deltayoverh[n-3] - Deltayoverh[n-2]) / (1. + ratio);
  p1[n-2] = -p1[n-3] + 2.*Deltayoverh[n-3];
  for(int i=n-4; i>=0; i--) {
    p1[i] = -p1[i+1] + 2.*Deltayoverh[i];
  }
  p1[n-1] = (1. + ratio)*p1[n-2] - ratio*p1[n-3];

  /* Deducing the p2's */
  gsl_vector* vectp2 = gsl_vector_alloc(n);
  double* p2 = vectp2->data;
  for(int i=0; i<=n-2; i++) {
    p2[i] = (p1[i+1] - p1[i]) / (2.*h[i]);
  }
  /* Note: p2[n-1] is set to values coherent with the derivatives of the spline at the last point, but not stricly speaking a coefficient of the spline. */
  p2[n-1] = p2[n-2];

  /* Copying the results in the output matrix */
  gsl_vector_view viewx = gsl_matrix_column(splinecoeffs, 0);
  gsl_vector_view viewp0 = gsl_matrix_column(splinecoeffs, 1);
  gsl_vector_view viewp1 = gsl_matrix_column(splinecoeffs, 2);
  gsl_vector_view viewp2 = gsl_matrix_column(splinecoeffs, 3);
  gsl_vector_memcpy(&viewx.vector, vectx);
  gsl_vector_memcpy(&viewp0.vector, vecty);
  gsl_vector_memcpy(&viewp1.vector, vectp1);
  gsl_vector_memcpy(&viewp2.vector, vectp2);

  /* Cleanup*/
  gsl_vector_free(vecth);
  gsl_vector_free(vectDeltay);
  gsl_vector_free(vectDeltayoverh);
  gsl_vector_free(vectp1);
  gsl_vector_free(vectp2);
}

void BuildSplineCoeffs(
  CAmpPhaseSpline** splines,                  /*  */
  CAmpPhaseFrequencySeries* freqseries)       /*  */
{
  /* Initialize output structure */
  int n = (int) freqseries->freq->size;
  CAmpPhaseSpline_Init(splines, n);

  /* Build the splines */
  BuildNotAKnotSpline((*splines)->spline_amp_real, freqseries->freq, freqseries->amp_real, n);
  BuildNotAKnotSpline((*splines)->spline_amp_imag, freqseries->freq, freqseries->amp_imag, n);
  BuildQuadSpline((*splines)->quadspline_phase, freqseries->freq, freqseries->phase, n);
}

void BuildListmodesCAmpPhaseSpline(
  ListmodesCAmpPhaseSpline** listspline,              /* Output: list of modes of splines in matrix form */
  ListmodesCAmpPhaseFrequencySeries* listh)           /* Input: list of modes in amplitude/phase form */
{
    if(*listspline){ /* We don't allow for the case where listspline already points to something */
      printf("Error: Tried to add a mode to an already existing ListmodesCAmpPhaseSpline ");
      exit(1);
    }
    else {
      ListmodesCAmpPhaseFrequencySeries* listelementh = listh;
      while(listelementh) {
	int l = listelementh->l;
	int m = listelementh->m;
	CAmpPhaseSpline* splines = NULL;
	BuildSplineCoeffs(&splines, listelementh->freqseries);
	*listspline = ListmodesCAmpPhaseSpline_AddModeNoCopy(*listspline, splines, l, m);
	listelementh = listelementh->next;
      }
    }
}

/* Note: for the spines in matrix form, the first column contains the x values, so the coeffs start at 1 */
double EvalCubic(
  gsl_vector* coeffs,  /**/
  double eps,          /**/
  double eps2,         /**/
  double eps3)         /**/
{
  double p0 = gsl_vector_get(coeffs, 1);
  double p1 = gsl_vector_get(coeffs, 2);
  double p2 = gsl_vector_get(coeffs, 3);
  double p3 = gsl_vector_get(coeffs, 4);
  return p0 + p1*eps + p2*eps2 + p3*eps3;
}

double EvalQuad(
  gsl_vector* coeffs,  /**/
  double eps,          /**/
  double eps2)         /**/
{
  double p0 = gsl_vector_get(coeffs, 1);
  double p1 = gsl_vector_get(coeffs, 2);
  double p2 = gsl_vector_get(coeffs, 3);
  return p0 + p1*eps + p2*eps2;
}

void EvalCAmpPhaseSpline(
  CAmpPhaseSpline* splines,                    //input
  CAmpPhaseFrequencySeries* freqseries)  //in/out defines CAmpPhase from defined freqs
{
  int ispline=0;
  //printf("Enter: n=%i\n",freqseries->freq->size);
  for(int i=0;i<freqseries->freq->size;i++){
    double f=gsl_vector_get(freqseries->freq,i);
    //printf(" f:%g < %g < %g\n",gsl_matrix_get(splines->quadspline_phase, 0, 0),f,gsl_matrix_get(splines->quadspline_phase, splines->quadspline_phase->size1-1, 0));

    /* Adjust the index in the spline if necessary and compute */
    while(gsl_matrix_get(splines->quadspline_phase, ispline+1, 0)<f)ispline++;

    double eps = f - gsl_matrix_get(splines->quadspline_phase, ispline, 0);
    double eps2 = eps*eps;
    double eps3 = eps2*eps;
    gsl_vector_view coeffsampreal = gsl_matrix_row(splines->spline_amp_real, ispline);
    gsl_vector_view coeffsampimag = gsl_matrix_row(splines->spline_amp_imag, ispline);
    gsl_vector_view coeffsphase = gsl_matrix_row(splines->quadspline_phase, ispline);
    double Ar = EvalCubic(&coeffsampreal.vector, eps, eps2, eps3);
    double Ai = EvalCubic(&coeffsampimag.vector, eps, eps2, eps3);
    double Ph = EvalQuad(&coeffsphase.vector, eps, eps2);

    gsl_vector_set(freqseries->amp_real,i,Ar);
    gsl_vector_set(freqseries->amp_imag,i,Ai);
    gsl_vector_set(freqseries->phase,i,Ph);
  }
};

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
/**
 * \file Component.h
 */

#ifndef ATK_MODELLING_COMPONENT_H
#define ATK_MODELLING_COMPONENT_H

#include <ATK/Modelling/config.h>
#include <ATK/Modelling/Types.h>

#include <gsl/gsl>

#include <Eigen/Eigen>

#include <tuple>
#include <vector>

namespace ATK
{
  template<typename DataType_>
  class DynamicModellerFilter;

  enum class Linear {
    IsLinear,
    IsNotLinear
  };

  /// Base class for all components
  template<typename DataType_>
  class ATK_MODELLING_EXPORT Component
  {
  public:
    using DataType = DataType_;

  protected:
    /// Local pins to which this component is connected to
    std::vector<std::tuple<PinType, gsl::index>> pins;

    /// dt value, not useful for every component
    DataType dt;

    /// The current modeller where the component is located
    DynamicModellerFilter<DataType>* modeller;

  public:
    /// Virtual destructor
    virtual ~Component();

    Component() = default;
    Component(const Component&) = delete;
    Component& operator=(const Component&) = delete;

    virtual Linear isLinear() const;

    /**
     * sets the pins for the component
     * @params pins is the set of pins for this component
     */
    void set_pins(std::vector<std::tuple<PinType, gsl::index>> pins);

    /// Returns the pins of this component
    const std::vector<std::tuple<PinType, gsl::index>>& get_pins() const
    {
      return pins;
    }

    /**
     * Used to indicate if the modeller needs to update its set of equations with those provided by this component
     * @param modeller the modeller to update
     */
    virtual void update_model(DynamicModellerFilter<DataType>* modeller);

    /**
     * Update the component for its steady state condition
     * @param dt is the delat that will be used in following updates
     */
    virtual void update_steady_state(DataType dt);

    /**
     * Update the component for its current state condition
     */
    virtual void update_state();

    /**
     * Precompute internal value before asking current and gradients
     * @param steady_state is a flag to indcate steady state computation (used for some components)
     */
    virtual void precompute(bool steady_state);

    /**
     * Get current for the given pin based on the state
     * @param pin_index is the pin from which to compute the current
     * @param steady_state is a flag to indcate steady state computation (used for some components)
     */
    virtual DataType get_current(gsl::index pin_index, bool steady_state) const = 0;

    /**
     * Get current gradient for the given pins based on the state
     * @param pin_index_ref is the pin of the current from which the gradient is computed
     * @param pin_index is the pin from which to compute the gradient of the pin_index current
     * @param steady_state is a flag to indcate steady state computation (used for some components)
     */
    virtual DataType get_gradient(gsl::index pin_index_ref, gsl::index pin_index, bool steady_state) const = 0;

    /**
     * Add a new equation to the modeller
     * @param eq_index is equation row in the state vector and the jacobian
     * @param eq_number is the equation number to add
     * @param eqs is the state vector to update
     * @param jacobian is the jacobian to update
     * @param steady_state is a flag to indcate steady state computation (used for some components)
     */
    virtual void add_equation(gsl::index eq_index, gsl::index eq_number, Eigen::Matrix<DataType, Eigen::Dynamic, 1>& eqs, Eigen::Matrix<DataType, Eigen::Dynamic, Eigen::Dynamic>& jacobian, bool steady_state) const;

    virtual gsl::index get_number_parameters() const;

    /// Get the name of a parameter
    virtual std::string get_parameter_name(gsl::index identifier) const;

    /// Get the value of a parameter
    virtual DataType_ get_parameter(gsl::index identifier) const;

    /// Set the value of a parameter
    virtual void set_parameter(gsl::index identifier, DataType_ value);
  };
}

#endif

ext: h

max_stars_repo_name: AudioTK/ATK-Modelling-free

max_issues_repo_name: AudioTK/ATK-Modelling-free

max_forks_repo_name: AudioTK/ATK-Modelling-free


text:
#include <iostream>
#include <vector>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <time.h>
#include <math.h>
//#include <C:/Program Files (x86)/Intel/Composer XE/compiler/include/cilk/cilk.h>
#ifndef _HEART_H
#define _HEART_H

class heart
{
private:
	int L;  //sytem size length for now square

	std::vector<std::vector<int>> state;  //value indicates whether in excited state refactory or resting
	std::vector<std::vector<int>> new_state;  //using new_state and state to distinguish newly excited cells from previously excited
	std::vector<std::vector<int>> functional; //1 if functional 0 if disfunctional
	std::vector<std::vector<int>> bond; //1 if bond to site above current cell, 0 if no bond above site
	std::vector<std::vector<int>> last_excited;
	std::vector<std::vector<int>> last_excited_counter; // have a dummy counter so as to not interfere with APD dynamics
	std::vector<std::vector<int>> refractory;
	std::vector<std::vector<int>> total_rest; //have counter for total time cells spend in rest or refractory

	int heart_time; //internal clock
	//model parameters
	int tau; //refactory period
	int T; // pacemaker rate
	double nu, epsilom, delta; //probability vertical connection exists, probability of misfire given faulty and probability cell faulty

	//random number generator
	const gsl_rng_type*Type;
	gsl_rng*r;

	int total_bonds;
	int total_functional;
	int total_excited;
	double refractory_period;

public:

	heart()
	{

	}
	heart(int _L, int _T, int _tau, double _nu, double _epsilom, double _delta)
	{

		refractory_period =0;
		L = _L;
		nu = _nu;
		tau = _tau;
		T = _T;
		epsilom = _epsilom;
		delta = _delta;

		heart_time = T;

		state.resize(L);
		new_state.resize(L);
		functional.resize(L);
		bond.resize(L);
		last_excited.resize(L);
		last_excited_counter.resize(L);
		refractory.resize(L);
		total_rest.resize(L);
		//pacemaker action
		for (unsigned i = 0; i<L; i++)
		{
			state[i].resize(L, 0);  // all cells eligible for excitation
			new_state[i].resize(L, 0);
			functional[i].resize(L, 1);// none disfunctional at the moment
			bond[i].resize(L, 1); // all transverse bonds present at the moment
			last_excited[i].resize(L, T);
			last_excited_counter[i].resize(L, 0);
			refractory[i].resize(L, tau);
			total_rest[i].resize(L,0);
		}
		//	total_functional = L*L;
		//	total_bonds = L*L;
		total_excited = 0;//could be useful avoiding calling size of some vector

		//HAVE NOT USED gsl_rng in class before need to make sure is working correctly

		//  environment variable GSL_RNG_TYPE
		gsl_rng_env_setup();

		Type = gsl_rng_mt19937;//random number generator using the Mersenne Twister algorithm
		r = gsl_rng_alloc(Type);

		//gsl_rng_set(r, time(NULL));//random number generator seeded with time
		gsl_rng_set(r, 10);//random number generator seeded with time

		//assigning cells as having a bond above them or not and functional or not functional
		int healthy = 1;

		for (int j = 0; j<L; j++)
		{
			for (int i = 0; i<L; i++)
			{

				bond[i][j] = (gsl_rng_uniform(r)<nu);
				functional[i][j] = (gsl_rng_uniform(r)>delta);// inequality reversed
				//double x =  gsl_rng_uniform(r);
				//	std::cout <<x<< '\t' << (x<nu) << std::endl;

			}
		}

		gsl_rng_set(r, 20);//random number generator seeded with time
	}
	void operate()//one time step
	{
		total_excited = 0;
		refractory_period = 0;
		//Relax excited cells
		for (int i = 0; i < L; i++)
		{
			for (int j = 0; j<L; j++)
			{
				if (state[i][j] == 1)//searching 4 neighbours of excited cell trying to excite if resting and if don't misfire and if there is a bond
				{
					new_state[((i + 1)<L)*(i + 1)][j] += (new_state[((i + 1)<L)*(i + 1)][j]<1)*(state[((i + 1)<L)*(i + 1)][j]<1)*bond[i][j] * (1 - (1 - functional[((i + 1)<L)*(i + 1)][j])*(gsl_rng_uniform(r)<epsilom));

					new_state[(i - 1)*(i>0) + (L - 1)*(i<1)][j] += (new_state[(i - 1)*(i>0) + (L - 1)*(i<1)][j]<1)*(state[(i - 1)*(i>0) + (L - 1)*(i<1)][j]<1)*bond[(i - 1)*(i>0) + (L - 1)*(i<1)][j] * (1 - (1 - functional[(i - 1)*(i>0) + (L - 1)*(i<1)][j])*(gsl_rng_uniform(r)<epsilom));

					new_state[i][(j + 1)*(j<L - 1)] += (new_state[i][(j + 1)*(j<L - 1)]<1)*(j<L - 1)*(state[i][(j + 1)*(j<L - 1)]<1)*(1 - (1 - functional[i][(j + 1)*(j<L - 1)])*(gsl_rng_uniform(r)<epsilom));

					new_state[i][(j - 1)*(j>0)] += (new_state[i][(j - 1)*(j>0)]<1)*(j>0)*(state[i][(j - 1)*(j>0)]<1)*(1 - (1 - functional[i][(j - 1)*(j>0)])*(gsl_rng_uniform(r)<epsilom));

					refractory[i][j] = restitution(i, j);

					total_rest[i][j] += last_excited_counter[i][j];
				}
			}
		}
		//end relaxing
		//std::vector<std::pair<int,int>> next_excited_cells;

		//evolving cells that were not excited previously or now (some subtlty to those possibly excited now)
		for (int i = 0; i<L; i++)
		{
			for (int j = 0; j<L; j++)
			{
				new_state[i][j] += (state[i][j]>0);
				//if state was excited or refactory before +1 states that were resting may have been newly excited or not but this ignores them
				new_state[i][j] *= (state[i][j]<refractory[i][j]);//promote to resting after end of refactory which is now variable
				total_excited += (new_state[i][j]<2)*new_state[i][j];

				last_excited[i][j] += (state[i][j] < 1) + (state[i][j]>1);//counts if not excited
				last_excited[i][j] *= (state[i][j] < 1) + (state[i][j]>1); //cycles back to 0 if excited;  old states so don't put lastexcited to 0 before calculating new refractory
				last_excited_counter[i][j] += (state[i][j] < 1) + (state[i][j]>1);//dummy counter
				last_excited_counter[i][j] *= (state[i][j] < 1) + (state[i][j]>1); //dummy counter
				refractory_period += refractory[i][j];

			}
		}

		//SAN
		for (int i = 0; i<L; i++)
		{
			//new_state[i][0] += (state[i][0]<1)*(heart_time>T - 1)*(1 - (1 - functional[i][0])*(gsl_rng_uniform(r)<epsilom));
			new_state[i][0] += (state[i][0]<1)*(heart_time>T - 1);

		}

		heart_time++;
		heart_time *= (heart_time<T + 1);
		//end SAN

		state = new_state;

	}
	int get_state(int i, int j)//perhaps have alternative for 2d array
	{

		return state[i][j];

	}
	void set_state(int i, int j, int val)
	{

		state[i][j] = val;
	}
	int get_function(int i, int j)
	{
		return functional[i][j];
	}
	void set_function(int i, int j, int val)
	{
		functional[i][j] = val;
	}
	int get_bond(int i, int j)
	{
		return bond[i][j];
	}
	void set_bond(int i, int j, int val)
	{
		bond[i][j] = val;
	}
	void update_total_rest(int i, int j)
	{
		total_rest[i][j] += last_excited_counter[i][j];
	}
	double get_total_rest(int i,int j)
	{
		return total_rest[i][j];
	}
	void reset_total_rest(int i, int j)
	{
		total_rest[i][j] = 0;
		last_excited_counter[i][j] = 0;
	}

	int get_tau()
	{
		return tau;
	}
	int get_T()
	{
		return T;
	}

	int get_excited()
	{
		return total_excited;
	}
	int get_refractory(int i, int j)
	{
		return refractory[i][j];
	}
	int get_last_excited(int i, int j)
	{
		return last_excited[i][j];
	}
	double colour(int i, int j)
	{
		return (1 - state[i][j] / (1.0*tau))*(state[i][j]>1);
		//return state[i][j]/(1.0*T);

	}
	int restitution(int i, int j)
	{
		//some function of tau T and last_excited[i][j]
		if (last_excited[i][j] < T)
		{
			//return floor((1 - exp(-4.0*last_excited[i][j] / (1.0*T)))*tau);
			return tau;
		}
		else
		{
			return tau;
		}
		//if T since excited should be approx tau
		//if less shoud be less than tau

	}
	double get_avg_refractory()
	{

		return refractory_period/(L*L);
	}
};
#endif

max_stars_repo_path: standard model + apd  100% working/heart_james.h