text:
#include <pygsl/solver.h>
#include <pygsl/utils.h>
#include <gsl/gsl_odeiv.h>

static const char odeiv_step_type_name   [] = "Odeiv-Step";
static const char odeiv_control_type_name[] = "Odeiv-Control";
static const char odeiv_evolve_type_name [] = "Odeiv-Evolve";

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

struct _mycontrol{
     gsl_odeiv_control *control;
     gsl_odeiv_step    *step;
     PyGSL_solver      *step_ob;
};
typedef struct _mycontrol mycontrol;

struct _myevolve{
     gsl_odeiv_evolve  *evolve;
     gsl_odeiv_control *control;
     gsl_odeiv_step    *step;
     PyGSL_solver      *control_ob;
     PyGSL_solver      *step_ob;
};
typedef struct _myevolve myevolve;

void
_mycontrol_free(mycontrol *c)
{
     FUNC_MESS_BEGIN();
     gsl_odeiv_control_free(c->control);
     if(c->step_ob){
	  DEBUG_MESS(3, "Decreasing step @ %p refcont %d", c->step_ob,
		     c->step_ob->ob_refcnt);
	  Py_DECREF(c->step_ob);
     }else{
	  DEBUG_MESS(3, "Freeing GSL Step @ %p", c->step);
	  gsl_odeiv_step_free(c->step);
     }
     memset(c, 0, sizeof(mycontrol));
     free(c);
     c = NULL;
     FUNC_MESS_END();
}

void
_myevolve_free(myevolve *c)
{
     FUNC_MESS_BEGIN();
     gsl_odeiv_evolve_free(c->evolve);
     if(c->control_ob){
	  DEBUG_MESS(3, "Decreasing control @ %p refcont %d", c->control_ob,
		     c->control_ob->ob_refcnt);
	  Py_DECREF(c->control_ob);
     }else{
	  DEBUG_MESS(3, "Freeing GSL Control @ %p", c->control);
	  gsl_odeiv_control_free(c->control);
     }
     if(c->step_ob){
	  DEBUG_MESS(3, "Decreasing step @ %p refcont %d", c->step_ob,
		     c->step_ob->ob_refcnt);
	  Py_DECREF(c->step_ob);
     }else{
	  DEBUG_MESS(3, "Freeing GSL Step @ %p", c->step);
	  gsl_odeiv_step_free(c->step);
     }
     memset(c, 0, sizeof(myevolve));
     free(c);
     c = NULL;
     FUNC_MESS_END();
}

const char *
PyGSL_mycontrol_getname(void *v)
{
     mycontrol *c = (mycontrol *) v;
     return gsl_odeiv_control_name(c->control);
}

#define _PyGSL_ODEIV_GENERIC_CHECK(ob, type_desc) \
  ((PyGSL_solver_check((ob))) && (((PyGSL_solver *)ob)->mstatic->type_name == (type_desc)))

#define PyGSL_ODEIV_STEP_Check(ob)    _PyGSL_ODEIV_GENERIC_CHECK((ob), odeiv_step_type_name)
#define PyGSL_ODEIV_CONTROL_Check(ob) _PyGSL_ODEIV_GENERIC_CHECK((ob), odeiv_control_type_name)
#define PyGSL_ODEIV_EVOLVE_Check(ob)  _PyGSL_ODEIV_GENERIC_CHECK((ob), odeiv_evolve_type_name)

static const char *this_file = __FILE__;
static char odeiv_step_init_err_msg[] = "odeiv_step.__init__";
static char odeiv_step_apply_doc[] =  "XXX Documentation missing\n";
static char odeiv_step_order_doc[] =  "XXX Documentation missing\n";

static char odeiv_control_hadjust_doc[] =  "XXX Documentation missing\n";
static char odeiv_evolve_apply_doc[] =  "XXX Documentation missing\n";

static char  PyGSL_odeiv_step_doc[] = "XXX Documentation missing\n";
static char  PyGSL_odeiv_control_doc[] = "XXX Documentation missing\n";
static char  PyGSL_odeiv_evolve_doc[] = "XXX Documentation missing\n";

/*---------------------------------------------------------------------------
  Wrapper functions to push call the approbriate Python Objects
  ---------------------------------------------------------------------------*/
static int
PyGSL_odeiv_func(double t, const double y[], double f[], void *params)
{
    int dimension, flag = GSL_FAILURE;

    PyObject *arglist = NULL, *result = NULL;
    PyArrayObject *yo = NULL;
    PyGSL_solver * step;
    gsl_vector_view yv, fv;
    PyGSL_error_info  info;

    FUNC_MESS_BEGIN();

    step = (PyGSL_solver *) params;
    if(!PyGSL_ODEIV_STEP_Check(step)){
	  PyGSL_add_traceback(module, this_file, __FUNCTION__,
			      __LINE__ - 2);
	  pygsl_error("Param not a step type!",
		    this_file, __LINE__ -2, GSL_EFAULT);
	  goto fail;
    }
    dimension =  step->problem_dimensions[0];

    /* Do I need to copy the array ??? */
    yv = gsl_vector_view_array((double *) y, dimension);
    yo = PyGSL_copy_gslvector_to_pyarray(&yv.vector);
    if (yo == NULL) goto fail;

    FUNC_MESS("\t\tBuild args");
    arglist = Py_BuildValue("(dOO)", t, yo, step->args);
    FUNC_MESS("\t\tEnd Build args");

    info.callback = step->cbs[0];
    info.message  = "odeiv_func";
    result  = PyEval_CallObject(step->cbs[0], arglist);

    if((flag = PyGSL_CHECK_PYTHON_RETURN(result, 1, &info)) != GSL_SUCCESS){
	  goto fail;
     }
    info.argnum = 1;
    fv = gsl_vector_view_array(f, dimension);
    if((flag = PyGSL_copy_pyarray_to_gslvector(&fv.vector, result, dimension,
					       &info)) != GSL_SUCCESS){
	  goto fail;
     }

    Py_DECREF(arglist);    arglist = NULL;
    Py_DECREF(yo);         yo = NULL;
    Py_DECREF(result);     result = NULL;
    FUNC_MESS_END();
    return GSL_SUCCESS;

 fail:
    FUNC_MESS("    IN Fail BEGIN");
    Py_XDECREF(yo);
    Py_XDECREF(result);
    Py_XDECREF(arglist);
    assert(flag != GSL_SUCCESS);
    FUNC_MESS("    IN Fail END");
    if(step->isset)
	 longjmp(step->buffer, flag);
    return flag;
}

static int
PyGSL_odeiv_jac(double t, const double y[], double *dfdy, double dfdt[],
		void *params)
{
    int dimension, flag = GSL_FAILURE;
    PyGSL_solver *step;
    PyGSL_error_info  info;

    PyObject *arglist = NULL, *result = NULL, *tmp=NULL;
    PyArrayObject *yo = NULL;

    gsl_vector_view yv, dfdtv;
    gsl_matrix_view dfdyv;

    FUNC_MESS_BEGIN();

    step = (PyGSL_solver *) params;
    if(!PyGSL_ODEIV_STEP_Check(step)){
	  PyGSL_add_traceback(module, this_file, __FUNCTION__,
			      __LINE__ - 2);
	  pygsl_error("Param not a step type!",
		    this_file, __LINE__ -2, GSL_EFAULT);
	  goto fail;
    }
    dimension = step->problem_dimensions[0];

    yv = gsl_vector_view_array((double *) y, dimension);
    yo = PyGSL_copy_gslvector_to_pyarray(&yv.vector);
    if (yo == NULL) goto fail;

    arglist = Py_BuildValue("(dOO)", t, yo, step->args);

    result  = PyEval_CallObject(step->cbs[1], arglist);

    info.callback = step->cbs[1];
    info.message  = "odeiv_jac";
    if((flag = PyGSL_CHECK_PYTHON_RETURN(result, 2, &info)) != GSL_SUCCESS){
	  goto fail;
     }

    info.argnum = 1;
    tmp = PyTuple_GET_ITEM(result, 0);
    dfdyv = gsl_matrix_view_array((double *) dfdy, dimension, dimension);
    if((flag = PyGSL_copy_pyarray_to_gslmatrix(&dfdyv.matrix, tmp, dimension, dimension, &info)) != GSL_SUCCESS){
	  goto fail;
    }

    info.argnum = 2;
    tmp = PyTuple_GET_ITEM(result, 1);
    dfdtv = gsl_vector_view_array((double *) dfdt, dimension);
    if((flag = PyGSL_copy_pyarray_to_gslvector(&dfdtv.vector, tmp, dimension, &info)) != GSL_SUCCESS){
	  goto fail;
    }

    Py_DECREF(arglist);    arglist = NULL;
    Py_DECREF(result);     result = NULL;
    Py_DECREF(yo);         yo = NULL;
    FUNC_MESS_END();
    return GSL_SUCCESS;
 fail:
    FUNC_MESS("IN Fail");
    assert(flag != GSL_SUCCESS);
    longjmp(step->buffer, flag);
    return flag;
}

/* Wrappers for the evaluation of the system */
static PyObject *
PyGSL_odeiv_step_apply(PyGSL_solver *self, PyObject *args)
{
    PyObject *result = NULL;
    PyObject *y0_o = NULL, *dydt_in_o = NULL;
    PyArrayObject *volatile y0 = NULL, * volatile yerr = NULL,
	 *volatile dydt_in = NULL, *volatile dydt_out = NULL,
	 *volatile yout = NULL;

    double t=0, h=0, *volatile dydt_in_d;
    int  r, flag;
    PyGSL_array_index_t dimension;

    FUNC_MESS_BEGIN();
    assert(PyGSL_ODEIV_STEP_Check(self));
    if(! PyArg_ParseTuple(args, "ddOOO", &t, &h, &y0_o,  &dydt_in_o)){
      return NULL;
    }

    dimension = self->problem_dimensions[0];
    y0 = PyGSL_vector_check(y0_o, dimension, PyGSL_DARRAY_CINPUT(1), NULL, NULL);
    if(y0 == NULL) goto fail;

    if (Py_None == dydt_in_o){
	 dydt_in_d = NULL;
    }else{
	 dydt_in = PyGSL_vector_check(dydt_in_o, dimension, PyGSL_DARRAY_CINPUT(2), NULL, NULL);
	 if(dydt_in == NULL) goto fail;
	 dydt_in_d = (double *) dydt_in->data;
    }

    dydt_out =  PyGSL_New_Array(1, &dimension, PyArray_DOUBLE);
    if (dydt_out == NULL) goto fail;

    yerr = PyGSL_New_Array(1, &dimension, PyArray_DOUBLE);
    if(yerr == NULL) goto fail;

    yout = (PyArrayObject *) PyGSL_Copy_Array(y0);
    if(yout == NULL) goto fail;

    self->isset = 0;
    if((flag=setjmp(self->buffer)) == 0){
	  FUNC_MESS("\t\t Setting Jmp Buffer");
	  self->isset = 1;
    } else {
	  FUNC_MESS("\t\t Returning from Jmp Buffer");
	  self->isset = 0;
	  goto fail;
    }

    r = gsl_odeiv_step_apply(self->solver, t, h,
			     (double *) yout->data,
			     (double *) yerr->data,
			     dydt_in_d,
			     (double *) dydt_out->data,
			     ((gsl_odeiv_system *)self->c_sys));
    self->isset = 0;
    if (GSL_SUCCESS != r){
	PyErr_SetString(PyExc_TypeError, "Error While evaluating gsl_odeiv");
      goto fail;
    }

    FUNC_MESS("    Returnlist create ");
    assert(yout != NULL);
    assert(yerr != NULL);
    assert(dydt_out != NULL);

    result = Py_BuildValue("(OOO)", yout, yerr, dydt_out);

    FUNC_MESS("    Memory free ");
    /* Deleting the arrays */
    Py_DECREF(y0);           y0 = NULL;
    Py_DECREF(yout);         yout = NULL;
    Py_DECREF(yerr);         yerr = NULL;
    Py_DECREF(dydt_out);     dydt_out = NULL;
    /* This array does not need to exist ... */
    Py_XDECREF(dydt_in);	 dydt_in=NULL;

    FUNC_MESS_END();
    return result;

    fail:
    FUNC_MESS("IN Fail");
    self->isset = 0;
    Py_XDECREF(y0);
    Py_XDECREF(yout);
    Py_XDECREF(yerr);
    Py_XDECREF(dydt_in);
    Py_XDECREF(dydt_out);
    FUNC_MESS("IN Fail End");
    return NULL;
}

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

  PyObject *result = NULL;
  PyObject *y0_o = NULL, *yerr_o = NULL, *dydt_o = NULL;
  PyArrayObject *y0 = NULL, *yerr = NULL, *dydt = NULL;
  double h = 0;
  int r = 0;
  mycontrol *c;

  PyGSL_array_index_t dimension = 0;

  FUNC_MESS_BEGIN();
  assert(PyGSL_ODEIV_CONTROL_Check(self));
  if(!PyArg_ParseTuple(args, "OOOd",  &y0_o, &yerr_o, &dydt_o, &h)){
    return NULL;
  }

  dimension = self->problem_dimensions[0];

  y0 = PyGSL_vector_check(y0_o, dimension, PyGSL_DARRAY_CINPUT(1), NULL, NULL);
  if(y0 == NULL)   goto fail;
  yerr = PyGSL_vector_check(yerr_o, dimension, PyGSL_DARRAY_CINPUT(2), NULL, NULL);
  if(yerr == NULL) goto fail;
  dydt = PyGSL_vector_check(dydt_o, dimension, PyGSL_DARRAY_CINPUT(3), NULL, NULL);
  if(dydt == NULL) goto fail;

  FUNC_MESS("      Array Pointers End");

  c = (mycontrol *) self->solver;
  r = gsl_odeiv_control_hadjust(c->control, c->step,
				(double *) y0->data,
				(double *) yerr->data,
				(double *) dydt->data, &h);

  FUNC_MESS("      Function End");
  Py_DECREF(y0);       y0 = NULL;
  Py_DECREF(yerr);     yerr = NULL;
  Py_DECREF(dydt);     dydt = NULL;

  result = Py_BuildValue("di",h,r);
  FUNC_MESS_END();
  return result;

 fail:
  FUNC_MESS("IN Fail");
  Py_XDECREF(y0);
  Py_XDECREF(yerr);
  Py_XDECREF(dydt);
  FUNC_MESS("IN Fail END");
  return NULL;
}

static PyObject *
PyGSL_odeiv_evolve_apply(PyGSL_solver *self, PyObject *args)
{
    PyObject *result = NULL;
    PyObject *y0_o = NULL,  *myargs = NULL;
    PyArrayObject *volatile y0 = NULL, *volatile yout = NULL;
    myevolve *e = NULL;

    double t=0, h=0, t1 = 0, flag;

    int dimension = self->problem_dimensions[0], r;

    assert(PyGSL_ODEIV_EVOLVE_Check(self));
    FUNC_MESS_BEGIN();

    if(!PyArg_ParseTuple(args, "dddOO", &t, &t1, &h, &y0_o, &myargs)){
      return NULL;
    }

    DEBUG_MESS(3, "y0_o @ %p", y0_o);

    y0 = PyGSL_vector_check(y0_o, dimension, PyGSL_DARRAY_CINPUT(1), NULL, NULL);
    if(y0 == NULL) goto fail;

    yout = (PyArrayObject *)  PyGSL_Copy_Array(y0);
    if(yout == NULL) goto fail;

    e = (myevolve *) self->solver;

    if((flag=setjmp(e->step_ob->buffer)) == 0){
	 e->step_ob->isset = 1;
	  FUNC_MESS("\t\t Setting Jmp Buffer");
     } else {
	  FUNC_MESS("\t\t Returning from Jmp Buffer");
	  e->step_ob->isset = 0;
	  goto fail;
     }
    DEBUG_MESS(3, "evolve @ %p\t control @ %p\t step @ %p", e, e->control, e->step);

    r = gsl_odeiv_evolve_apply(e->evolve,
			       e->control,
			       e->step,
			       e->step_ob->c_sys, &t, t1, &h,
			       (double * )yout->data);
   e->step_ob->isset = 0;
    if (GSL_SUCCESS != r){
	 goto fail;
    }

    assert(yout != NULL);

    result = Py_BuildValue("(ddO)", t, h, yout);

    /* Deleting the arrays */
    Py_DECREF(yout);     yout = NULL;
    Py_DECREF(y0);       y0=NULL;
    FUNC_MESS_END();
    return result;

 fail:
    FUNC_MESS("IN Fail");
    e->step_ob->isset = 0;
    PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__);
    Py_XDECREF(y0);
    Py_XDECREF(yout);
    FUNC_MESS("IN Fail End");
    return NULL;
}

GETINT(odeiv_step_order)

static struct PyMethodDef PyGSL_odeiv_step_methods[] = {
     {"apply", (PyCFunction) PyGSL_odeiv_step_apply, METH_VARARGS, odeiv_step_apply_doc},
     {"order", (PyCFunction) PyGSL_odeiv_step_order, METH_VARARGS, odeiv_step_order_doc},
     {NULL, NULL}
};
static struct PyMethodDef PyGSL_odeiv_control_methods[] = {
     {"hadjust", (PyCFunction) PyGSL_odeiv_control_hadjust, METH_VARARGS, odeiv_control_hadjust_doc},
     {NULL, NULL}
};

static struct PyMethodDef PyGSL_odeiv_evolve_methods[] = {
     {"apply", (PyCFunction) PyGSL_odeiv_evolve_apply, METH_VARARGS, odeiv_evolve_apply_doc},
     {NULL, NULL}
};

static const struct _SolverStatic
_StepMethods     = {{(void_m_t) gsl_odeiv_step_free,
		   (void_m_t) gsl_odeiv_step_reset,
		   (name_m_t) gsl_odeiv_step_name,
		   (int_m_t) NULL},
		    3, PyGSL_odeiv_step_methods, odeiv_step_type_name},
_ControlMethods  = {{(void_m_t) _mycontrol_free,
		      (void_m_t) NULL,
		      (name_m_t) PyGSL_mycontrol_getname,
		      (int_m_t) NULL},
		    3, PyGSL_odeiv_control_methods, odeiv_control_type_name},
_EvolveMethods   = {{(void_m_t) _myevolve_free,
		     (void_m_t) gsl_odeiv_evolve_reset,
		     (name_m_t) NULL,
		     (int_m_t) NULL,},
		    3, PyGSL_odeiv_evolve_methods, odeiv_evolve_type_name};

static PyObject *
PyGSL_odeiv_step_init(PyObject *self, PyObject *args, PyObject *kwdict, const gsl_odeiv_step_type * odeiv_type)
{

     PyObject *func=NULL, *jac=NULL, *o_params=NULL;
     PyGSL_solver *solver = NULL;

     static char * kwlist[] = {"dimension", "func", "jac", "args", NULL};
     int dim, has_jacobian = 0;
     gsl_odeiv_system * c_sys;

     solver_alloc_struct s = {odeiv_type, (void_an_t) gsl_odeiv_step_alloc,
			      &_StepMethods};

     FUNC_MESS_BEGIN();

     assert(args);
     if (0 == PyArg_ParseTupleAndKeywords(args, kwdict, "iOOO:odeiv_step.__init__", kwlist,
					  &dim, &func, &jac, &o_params)){
	  PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__ - 2);
	  return NULL;
     }
     if (dim <= 0){
	  PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__ - 1);
	  pygsl_error("The dimension of the problem must be at least 1",
		    this_file, __LINE__ -2, GSL_EDOM);
	  return NULL;
     }
     if(!PyCallable_Check(func)){
	  PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__ - 1);
	  pygsl_error("The function object is not callable!",
		    this_file, __LINE__ -2, GSL_EBADFUNC);
	  goto fail;
     }

     if(jac == Py_None){
	  if(odeiv_type == gsl_odeiv_step_bsimp){
	       PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__ - 1);
	       pygsl_error("The bsimp method needs a jacobian! You supplied None.",
			 this_file, __LINE__ -2, GSL_EBADFUNC);
	       goto fail;
	  }
     }else{
	  if(!PyCallable_Check(jac)){
	       PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__ - 1);
	       pygsl_error("The jacobian object must be None or callable!",
			 this_file, __LINE__ -2, GSL_EBADFUNC);
	       goto fail;
	  }
	  has_jacobian = 1;

     }

     solver = (PyGSL_solver *) PyGSL_solver_dn_init(self, args, &s, 3);

     if(solver == NULL){
	  goto fail;
     }
     DEBUG_MESS(3, "solver @ %p", solver);

     solver->solver =  gsl_odeiv_step_alloc(odeiv_type, dim);
     if(solver->solver == NULL){
	  goto fail;
     }
     DEBUG_MESS(3, "step @ %p", solver->solver);
     c_sys = (gsl_odeiv_system *)calloc(1, sizeof(gsl_odeiv_system));
     if(c_sys == NULL){
	  PyErr_NoMemory();
	  goto fail;
     }

     /* Need for cleanup in fail : */
     solver->c_sys = c_sys;
     DEBUG_MESS(3, "c_sys @ %p", solver->c_sys);
     solver->problem_dimensions[0] = dim;
     if(has_jacobian){
	  c_sys->jacobian = PyGSL_odeiv_jac;
	  if(!PyCallable_Check(jac))
	       goto fail;
	  solver->cbs[1] = jac;
     }else{
	  c_sys->jacobian = NULL;
	  solver->cbs[1] = NULL;
     }
     c_sys->function = PyGSL_odeiv_func;
     if(!PyCallable_Check(func))
	       goto fail;

     solver->cbs[0] = func;
     c_sys->params = (void *) solver;
     DEBUG_MESS(3, "params @ %p", c_sys->params);
     Py_INCREF(solver->cbs[0]);
     Py_XINCREF(solver->cbs[1]);
     Py_XINCREF(solver->args);

     solver->args  = o_params;
     Py_INCREF(solver->args);

     FUNC_MESS_END();
     return (PyObject *) solver;

 fail:
     FUNC_MESS("FAIL");
     Py_XDECREF(solver);
     return NULL;
}
#define ADD_ODESTEPPER(mytype)                                                   \
static PyObject *                                                                \
PyGSL_odeiv_step_init_ ## mytype (PyObject * self, PyObject * args, PyObject *kwdic)   \
{                                                                                \
     return PyGSL_odeiv_step_init(self, args, kwdic, gsl_odeiv_step_ ## mytype); \
}
ADD_ODESTEPPER(rk2)
ADD_ODESTEPPER(rk4)
ADD_ODESTEPPER(rkf45)
ADD_ODESTEPPER(rkck)
ADD_ODESTEPPER(rk8pd)
ADD_ODESTEPPER(rk2imp)
ADD_ODESTEPPER(rk4imp)
ADD_ODESTEPPER(bsimp)
ADD_ODESTEPPER(gear1)
ADD_ODESTEPPER(gear2)

static PyObject *
PyGSL_odeiv_control_init(PyObject *self, PyObject *args, void * type)
{
     int nargs = -1, tmp=0;
     double eps_abs, eps_rel, a_y, a_dydt;
     PyGSL_solver *step=NULL, *solver=NULL;
     mycontrol * c;

     gsl_odeiv_control *(*evaluator_5)(double , double , double , double ) = NULL;
     gsl_odeiv_control *(*evaluator_3)(double , double ) = NULL;

     solver_alloc_struct s = {type, (void_an_t) gsl_odeiv_control_alloc,
			      &_ControlMethods};

     FUNC_MESS_BEGIN();
     /* The arguments depend on the type of control */
     if(type == (void *) gsl_odeiv_control_standard_new){
	  /* step, eps_abs, eps_rel, a_y, a_dydt */
	  nargs = 5;
     }else if(type == (void *) gsl_odeiv_control_y_new ||
	      type == (void *) gsl_odeiv_control_yp_new){
	  nargs = 3;
     }else{
	  PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg,
			      __LINE__ - 2);
	  pygsl_error("Unknown control type",
		    this_file, __LINE__ -2, GSL_EFAULT);
	  goto fail;
     }
     assert(nargs > -1);

     switch(nargs){
     case 5:
	  tmp == PyArg_ParseTuple(args, "Odddd:odeiv_control.__init__",
				   &step, &eps_abs, &eps_rel, &a_y, &a_dydt);
	  break;
     case 3:
	  tmp == PyArg_ParseTuple(args, "Odd:odeiv_control.__init__",
				  &step, &eps_abs, &eps_rel);
	  break;
     default:
	  fprintf(stderr, "nargs = %d\n", nargs);
	  pygsl_error("Unknown number of arguments",
		    this_file, __LINE__ -2, GSL_EFAULT);
	  goto fail; break;
     }
     if(!PyGSL_ODEIV_STEP_Check(step)){
	  int flag;
	  flag = PyGSL_solver_check(step);
	  DEBUG_MESS(3, "is solver?  %d, %p %p ", flag,  PyGSL_API[PyGSL_solver_type_NUM], step->ob_type);
	  if(flag){
	       DEBUG_MESS(3, "solver = %s, %p !=  %p", step->mstatic->type_name, step->mstatic->type_name,
			  odeiv_step_type_name);
	       pygsl_error("First argument must be a step solver!", __FILE__, __LINE__, GSL_EINVAL);
	  }
	  goto fail;
     }

     if(tmp){
	  PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg,
			      __LINE__ - 2);
	  return NULL;
     }

     solver =  (PyGSL_solver *) PyGSL_solver_dn_init(self, args, &s, 3);
     if (NULL == solver){
	  PyErr_NoMemory();
	  goto fail;
     }
     c = calloc(1, sizeof(mycontrol));
     if(c == NULL){
	  PyErr_NoMemory();
	  goto fail;
     }
     solver->solver = c;
     switch(nargs){
     case 5:
	  evaluator_5 = type;
	  c->control = evaluator_5(eps_abs, eps_rel, a_y, a_dydt);
	  break;
     case 3:
	  evaluator_3 = type;
	  c->control = evaluator_3(eps_abs, eps_rel);
	  break;
     default:
	  goto fail;
     }
     if (NULL == c->control){
	  PyErr_NoMemory();
	  goto fail;
     }
     DEBUG_MESS(3, "c->control @ %p", c->control);
     c->step =  step->solver;
     c->step_ob = step;
     Py_INCREF(step);
     FUNC_MESS_END();
     return (PyObject *) solver;

 fail:
     FUNC_MESS("FAIL");
     Py_XDECREF(solver);
     return NULL;

}

#define ADD_ODECONTROL(name)                                                  \
static PyObject *                                                             \
PyGSL_odeiv_control_init_ ## name (PyObject * self, PyObject * args)          \
{                                                                             \
     return PyGSL_odeiv_control_init(self, args, (void *) gsl_odeiv_control_ ## name);    \
}
ADD_ODECONTROL(standard_new)
ADD_ODECONTROL(y_new)
ADD_ODECONTROL(yp_new)

static PyObject *
PyGSL_odeiv_evolve_init(PyObject *self, PyObject *args)
{
     PyGSL_solver *step, *control, *a_ev = NULL;
     myevolve *e;
     solver_alloc_struct s = {NULL, (void_an_t) gsl_odeiv_evolve_alloc,
			      &_EvolveMethods};

     /* step, control */
     FUNC_MESS_BEGIN();
     if(0== PyArg_ParseTuple(args, "OO:odeiv_evolve.__init__",
			     &step, &control)){
	  return NULL;
     }
     if(!PyGSL_ODEIV_STEP_Check(step)){
	  pygsl_error("First argument must be a step solver!", __FILE__, __LINE__, GSL_EINVAL);
	  goto fail;
     }

     if(!PyGSL_ODEIV_CONTROL_Check(control)){
	  pygsl_error("Second argument must be a control solver!", __FILE__, __LINE__, GSL_EINVAL);
	  goto fail;
     }

     a_ev =  (PyGSL_solver *) PyGSL_solver_dn_init(self, args, &s, 3);
     if(NULL == a_ev){
	  PyErr_NoMemory();
	  return NULL;
     }
     a_ev->problem_dimensions[0] = step->problem_dimensions[0];

     e = (myevolve *) calloc(1, sizeof(myevolve));
     if(e == NULL){
	  PyErr_NoMemory();
	  goto fail;
     }
     a_ev->solver = e;
     e->step_ob =  step;
     e->control_ob = control;
     Py_INCREF(step);
     Py_INCREF(control);
     e->step = step->solver;
     e->control = ((mycontrol *)control->solver)->control;
     e->evolve = gsl_odeiv_evolve_alloc(step->problem_dimensions[0]);
     if(NULL == e->evolve){
	  PyErr_NoMemory();
	  goto fail;
     }
     FUNC_MESS_END();
     return (PyObject *) a_ev;
 fail:
     FUNC_MESS("FAIL");
     Py_XDECREF(a_ev);
     return NULL;
}

static const char PyGSL_odeiv_module_doc [] = "XXX Missing ";
static PyMethodDef mMethods[] = {
     {"step_rk2",    (PyCFunction)PyGSL_odeiv_step_init_rk2,    METH_VARARGS|METH_KEYWORDS, PyGSL_odeiv_step_doc},
     {"step_rk4",    (PyCFunction)PyGSL_odeiv_step_init_rk4,    METH_VARARGS|METH_KEYWORDS, PyGSL_odeiv_step_doc},
     {"step_rkf45",  (PyCFunction)PyGSL_odeiv_step_init_rkf45,  METH_VARARGS|METH_KEYWORDS, PyGSL_odeiv_step_doc},
     {"step_rkck",   (PyCFunction)PyGSL_odeiv_step_init_rkck,   METH_VARARGS|METH_KEYWORDS, PyGSL_odeiv_step_doc},
     {"step_rk8pd",  (PyCFunction)PyGSL_odeiv_step_init_rk8pd,  METH_VARARGS|METH_KEYWORDS, PyGSL_odeiv_step_doc},
     {"step_rk2imp", (PyCFunction)PyGSL_odeiv_step_init_rk2imp, METH_VARARGS|METH_KEYWORDS, PyGSL_odeiv_step_doc},
     {"step_rk4imp", (PyCFunction)PyGSL_odeiv_step_init_rk4imp, METH_VARARGS|METH_KEYWORDS, PyGSL_odeiv_step_doc},
     {"step_bsimp",  (PyCFunction)PyGSL_odeiv_step_init_bsimp,  METH_VARARGS|METH_KEYWORDS, PyGSL_odeiv_step_doc},
     {"step_gear1",  (PyCFunction)PyGSL_odeiv_step_init_gear1,  METH_VARARGS|METH_KEYWORDS, PyGSL_odeiv_step_doc},
     {"step_gear2",  (PyCFunction)PyGSL_odeiv_step_init_gear2,  METH_VARARGS|METH_KEYWORDS, PyGSL_odeiv_step_doc},
     {"control_standard_new", PyGSL_odeiv_control_init_standard_new, METH_VARARGS, PyGSL_odeiv_control_doc},
     {"control_y_new",        PyGSL_odeiv_control_init_y_new,        METH_VARARGS, PyGSL_odeiv_control_doc},
     {"control_yp_new",       PyGSL_odeiv_control_init_yp_new,       METH_VARARGS, PyGSL_odeiv_control_doc},
     {"evolve", PyGSL_odeiv_evolve_init, METH_VARARGS, PyGSL_odeiv_evolve_doc},
     {NULL, NULL, 0, NULL}
};

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

     m=Py_InitModule("odeiv", 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_odeiv_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:
#include <stdio.h>
#include <gsl_rng.h>
#include <gsl_randist.h>

int main(int argv, char* argc[])
{
  // Regression test for https://github.com/conda-forge/ipopt-feedstock/issues/57
  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);

  if (r == NULL)
  {
    return 1;
  }

  gsl_rng_free(r);

  return 0;
}

ext: c

max_stars_repo_name: hmaarrfk/gsl-feedstock

max_issues_repo_name: hmaarrfk/gsl-feedstock

max_forks_repo_name: hmaarrfk/gsl-feedstock


text:
#pragma once
#include <memory>
#include <gsl/span>
#include "D3DManager.h"

namespace vrutil
{

class OverlayItem
{
    class OverlayItemImpl *m_impl;
    OverlayItem(class OverlayItemImpl *impl);

public:
    ~OverlayItem();
    static std::shared_ptr<OverlayItem> Create(uint64_t handle);
    void ProcessEvents();
    void SetTexture(const ComPtr<ID3D11Texture2D> &texture);
    void SetSharedHandle(void *handle);
    void SetOverlayWidthInMeters(float meter);
    void Show();
    void SetMatrix34(const float *matrix34);
};
using OverlayItemPtr = std::shared_ptr<OverlayItem>;

} // namespace vrutil

ext: h

max_stars_repo_name: ousttrue/VROverlaySample

max_issues_repo_name: ousttrue/VROverlaySample

max_forks_repo_name: ousttrue/VROverlaySample


text:
/* cdf/cdf_beta.c
 *
 * Copyright (C) 2003 Brian Gough.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_math.h>

#include "beta_inc.c"

double
gsl_cdf_beta_P (const double x, const double a, const double b)
{
  double P;

  if (x <= 0.0 )
    {
      return 0.0;
    }

  if ( x >= 1.0 )
    {
      return 1.0;
    }

  P = beta_inc_AXPY (1.0, 0.0, a, b, x);

  return P;
}

double
gsl_cdf_beta_Q (const double x, const double a, const double b)
{
  double Q;

  if ( x >= 1.0)
    {
      return 0.0;
    }

  if ( x <= 0.0 )
    {
      return 1.0;
    }

  Q = beta_inc_AXPY (-1.0, 1.0, a, b, x);

  return Q;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <unistd.h>
#include <stdio.h>
#include <gsl/gsl_math.h>

int main (void) {
    for ( ; ; ) {
    double x = 5.0;
    double y = gsl_expm1(100.0);
    printf ("J0(%g) = %.18e\n", x, y);
    sleep(2);
    }
    return 0;
}


text:
#ifndef SCALING_H_SAL4XICS
#define SCALING_H_SAL4XICS

#include <algorithm>
#include <gsl/gsl>
#include <type_traits>
#include <utility>

namespace sens_loc::math {

/// This class is a small wrapper to clarify the mathematical notation
/// of ranges. It does NOT distinguish between \f$[min, max)\f$,
/// \f$[min, max]\f$, \f$(min, max)\f$.
/// \tparam Real underlying type of the range, integer type for example
/// \warning this struct does not ensure that \p min is always smaller \p max.
template <typename Real>
struct numeric_range {
    static_assert(std::is_arithmetic_v<Real>);

    Real min;  ///< lower bound of the range
    Real max;  ///< uppoer bound of the range
};

/// This function scales \p value from the range \p source_range to the
/// new range \p target_range.
/// \note \p value is clamped to \p source_range
/// \pre the ranges are well formed (\p min is smaller \p max)
/// \returns the scaled value that is clamped to \p target_range.
template <typename Real>
inline constexpr Real scale(const numeric_range<Real>& source_range,
                            const numeric_range<Real>& target_range,
                            Real                       value) noexcept {
    static_assert(std::is_arithmetic_v<Real>);

    Expects(source_range.min < source_range.max);
    Expects(target_range.min < target_range.max);

    const Real clamped = std::clamp(value, source_range.min, source_range.max);

    Ensures(clamped <= source_range.max);
    Ensures(clamped >= source_range.min);

    const Real scaled =
        ((target_range.max - target_range.min) * (value - source_range.min)) /
            (source_range.max - source_range.min) +
        target_range.min;

    const Real target_clamped =
        std::clamp(scaled, target_range.min, target_range.max);

    Ensures(target_clamped <= target_range.max);
    Ensures(target_clamped >= target_range.min);

    return target_clamped;
}
}  // namespace sens_loc::math

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

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
#pragma once

#include <halley/maths/vector2.h>
#include <halley/maths/rect.h>
#include <memory>
#include <halley/resources/resource.h>
#include <halley/text/halleystring.h>
#include <halley/data_structures/hash_map.h>
#include <gsl/span>
#include "halley/maths/vector4.h"

namespace Halley
{
	class Resources;
	class Serializer;
	class Deserializer;
	class ResourceDataStatic;
	class Texture;
	class ResourceLoader;

	class SpriteSheetEntry
	{
	public:
		Vector2f pivot;
		Vector2i origPivot;
		Vector2f size;
		Rect4f coords;
		Vector4s trimBorder;
		Vector4s slices;
		int duration = 0;
		bool rotated = false;
		bool sliced = false;

		void serialize(Serializer& s) const;
		void deserialize(Deserializer& s);
	};

	class SpriteSheetFrameTag
	{
	public:
		String name;
		int from = 0;
		int to = 0;

		void serialize(Serializer& s) const;
		void deserialize(Deserializer& s);
	};

	class SpriteSheet : public Resource
	{
	public:
		const std::shared_ptr<const Texture>& getTexture() const;
		const SpriteSheetEntry& getSprite(const String& name) const;
		const SpriteSheetEntry& getSprite(size_t idx) const;

		const std::vector<SpriteSheetFrameTag>& getFrameTags() const;
		std::vector<String> getSpriteNames() const;

		size_t getSpriteCount() const;
		size_t getIndex(const String& name) const;
		bool hasSprite(const String& name) const;

		void loadJson(gsl::span<const gsl::byte> data);

		void addSprite(String name, const SpriteSheetEntry& sprite);
		void setTextureName(String name);

		static std::unique_ptr<SpriteSheet> loadResource(ResourceLoader& loader);
		constexpr static AssetType getAssetType() { return AssetType::SpriteSheet; }
		void reload(Resource&& resource) override;

		void serialize(Serializer& s) const;
		void deserialize(Deserializer& s);

	private:
		Resources* resources = nullptr;

		mutable std::shared_ptr<const Texture> texture;
		std::vector<SpriteSheetEntry> sprites;
		HashMap<String, uint32_t> spriteIdx;
		std::vector<SpriteSheetFrameTag> frameTags;
		String textureName;

		void loadTexture(Resources& resources) const;
	};

	class SpriteResource : public Resource
	{
	public:
		SpriteResource(std::shared_ptr<const SpriteSheet> spriteSheet, size_t idx);

		const SpriteSheetEntry& getSprite() const;
		size_t getIdx() const;
		std::shared_ptr<const SpriteSheet> getSpriteSheet() const;

		constexpr static AssetType getAssetType() { return AssetType::Sprite; }
		static std::unique_ptr<SpriteResource> loadResource(ResourceLoader& loader);
		void reload(Resource&& resource) override;

	private:
		std::weak_ptr<const SpriteSheet> spriteSheet;
		size_t idx = -1;
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


text:
/**
 *  \file matrix.h
 *  \brief a template class to describe the Matrix
 *  \author wbd
 *
 *  This class implements the template matrix. It contains the array "data" to
 *  hold the matrix elements.
 *  The matrix can be initialize with the matrix_initializer_list
 *  (a.k, recursive initializer_list).
 *  For example:
 *    Matrix<float, 2> amat {{1,2}, {2, 4}};
 *  The matrix can also initialize with the MatrixShape and the initialize value.
 *  For example:
 *    Matrix<float 2> bmat({2,2}, 0);
 *
 */
#ifndef MATRIX_H_
#define MATRIX_H_

#include <iostream>
#include <fstream>
#include <cblas.h>
#include <cstring>
#include <iomanip>
#include <memory>
#include <initializer_list>
#include <type_traits>
#include <array>
#include "expr.h"
#include "matrix_shape.h"
#include "sub_matrix.h"

//#define DEBUG
namespace matrix {

// ============================================================================
// recursive initializer_list for Matrix Class
template<typename DataType, size_t N>
struct MatrixInit {
  using type = std::initializer_list<typename MatrixInit<DataType, N-1>::type>;
};

template<typename DataType>
struct MatrixInit<DataType, 1> {
  using type = std::initializer_list<DataType>;
};

template<typename DataType>
struct MatrixInit<DataType, 0> ;

template<typename DataType, size_t N>
using matrix_initializer_list = typename MatrixInit<DataType, N>::type;

template<typename DataType, size_t N>
class Matrix : public ExprBase<Matrix<DataType, N>, DataType> {
 public:

  /**
   * default constructor
   */
  Matrix()
      : data(nullptr),
        stride(size_t(0)),
        capicity(size_t(0)),
        row(size_t(0)),
        column(size_t(0)) {
  }
  ;

  /**
   * constructor with matrix shape and initial value
   * @param s is shape of matrix
   * @param n is initial value for the elements in the matrix
   */
  inline Matrix(const MatrixShape<N> &s, const DataType n);

  /**
   * constructor with matrix shape, stride and  initial value
   * @param s is shape of the matrix
   * @param st is stride for each row of the multi-array(Note that: change the shape[N-1] memeory )
   * @param n is initial value for the elements in the matrix
   */
  inline Matrix(const MatrixShape<N> &s, const size_t st, const DataType n);

  /**
   * copy constructor
   * @param m is another matrix to be cloned to this matrix
   */
  inline Matrix(const Matrix<DataType, N> &m);

  /**
   * copy assignment
   * @param m is another matrix to be cloned to this matrix
   * @return the matrix that has copied data from another matrix
   */
  inline Matrix<DataType, N> & operator =(const Matrix<DataType, N> &m);

  /**
   * copy constructor from SubMatrix
   * @param m is the submatrix to be cloned to this matrix
   */
  inline Matrix(const SubMatrix<DataType, N> &m);

  /**
   * copy assignment from SubMatrix
   * @param m is another submatrix to be cloned to this matrix
   * @return the matrix that has copied data from another matrix
   */
  inline Matrix<DataType, N> & operator =(const SubMatrix<DataType, N> &m);

  /**
   * move constructor
   * @param m is another matrix to be moved to this matrix
   */
  inline Matrix(Matrix<DataType, N> &&m);

  /**
   * move assignment operator
   * @param m is matrix to be moved to this matrix
   * @return the matrix that has gotten data from another matrix
   */
  inline Matrix<DataType, N> & operator =(Matrix<DataType, N> &&m);

  /**
   *  constructor with matrix_initializer_list
   * @param t is initializer_list for the matrix
   */
  inline Matrix(matrix_initializer_list<DataType, N> t);

  /**
   *  assignment operator with matrix_initializer_list
   * @param t is initializer_list for the matrix
   * @return the matrix that has filled the data with the initializer_list
   */
  inline Matrix & operator =(matrix_initializer_list<DataType, N> t);
  template<typename T>
  inline Matrix(std::initializer_list<T> t) = delete;
  template<typename T>
  inline Matrix & operator =(std::initializer_list<T> t) = delete;

  /**
   * de-constructor: delete the array
   */
  ~Matrix() {
#ifdef DEBUG
    std::cout<<"decons!!!!!!!!!!!!!!!!!!!!" << this << std::endl;
#endif
    delete[] data;
  }

  /**
   *  index operation for the matrix
   * @param i is the index
   * @return a SubMatrix Object which is refer to some elements in the Matrix Object
   */
  inline SubMatrix<DataType, N - 1> operator[](size_t i) const;
  inline SubMatrix<DataType, N - 1> operator[](size_t i);
  /**
   * slice function for the matrix
   * @param i is the start index
   * @param j is the end index
   * @return a SubMatrix Object which is refer to some elements in the Matrix Object
   */
  inline SubMatrix<DataType, N> slice(size_t i, size_t j) const;
  inline SubMatrix<DataType, N> slice(size_t i, size_t j);

  /**
   * Scalar add Operator
   * @param n is the scalar added to matrix
   * @return the matrix in which the elements are all added the scalar value
   */
  inline Matrix<DataType, N> & operator +=(const DataType & n);

  /**
   * Scalar sub Operator
   * @param n is the scalar subscribed to matrix
   * @return the matrix in which the elements are all subscribed the scalar value
   */
  inline Matrix<DataType, N> & operator -=(const DataType & n);

  /**
   * Scalar multiply Operator
   * @param n is the scalar multiplied to matrix
   * @return the matrix in which the elements are all multiplied the scalar value
   */
  inline Matrix<DataType, N> & operator *=(const DataType & n);

  /**
   * Scalar devision Operator
   * @param n is the scalar divided to matrix
   * @return the matrix in which the elements are all devided the scalar value
   */
  inline Matrix<DataType, N> & operator /=(const DataType & n);

  /**
   * add Operator
   * @param t is the matrix added to this matrix
   * @return the matrix in which the elements: m1(i,j) += m2(i,j)
   */
  inline Matrix<DataType, N> & operator +=(const Matrix<DataType, N> & t);

  /**
   * sub Operator
   * @param t is the matrix subscribed to this matrix
   * @return the matrix in which the elements: m1(i,j) -= m2(i,j)
   */
  inline Matrix<DataType, N> & operator -=(const Matrix<DataType, N> & t);

  /**
   * multiply Operator
   * @param t is the matrix multiplied to this matrix
   * @return the matrix in which the elements: m1(i,j) *= m2(i,j)
   */
  inline Matrix<DataType, N> & operator *=(const Matrix<DataType, N> & t);

  /**
   * division opertor
   * @param t is the matrix divided to this matrix
   * @return the matrix in which the elements: m1(i,j) /= m2(i,j)
   */
  inline Matrix<DataType, N>& operator /=(const Matrix<DataType, N> & t);

  /**
   * Assign Operator
   * @param e is the right matrix object or a scalar
   * @return the result matrix
   */
  template<typename SubType>
  inline Matrix<DataType, N>& operator=(const ExprBase<SubType, DataType> &e);

  /**
   * get the basic element in index (i,j)
   * @param i is the row index
   * @param j is the column index
   * @return the element in the index(i,j)
   */
  inline DataType eval(size_t i, size_t j) const;

  inline void set_row_ele(int i, const Matrix<DataType, N> & s);

  /**
   * get the column size
   * @return the column size
   */
  inline size_t getColumn() const {
    return column;
  }

  /**
   * get the row size
   * @return the row size
   */
  inline size_t getRow() const {
    return row;
  }

  /**
   * get the capicity of the matrix
   * @return the capicity of the matrix
   */
  inline size_t getCapicity() const {
    return capicity;
  }

  /**
   * get the memory length for the matrix
   * @return the capicity of the matrix
   */
  inline int get_capicity() const;
  /**
   * get the number of the elemments in the matrix
   * @return the matrix elements number
   */
  inline int get_size() const;
  /**
   * get the length of the shape
   * @param s is the shape object
   * @param stride is the row length
   * @return the length of the shape with the stride
   */
  inline size_t get_length(const MatrixShape<N> &s, size_t stride) const;
  /**
   * @return the pointer to the data of the matrix
   */
  inline DataType * get_data() {
    return data;
  }
  /**
   * @return the pointer to the data of the matrix
   */
  inline DataType * get_data() const {
    return data;
  }
  //private:
  MatrixShape<N> shape;
  size_t stride;
  size_t capicity;
  DataType * data;
  size_t row;
  size_t column;
};

//=============================================================================
/**
 * One-dimension Matrix Class
 * For example:
 *  Matrix<float, 1> m {1,2,3,4}
 *  support  operator [], such as m[0].
 */
template<typename DataType>
class Matrix<DataType, 1> : public ExprBase<Matrix<DataType, 1>, DataType> {
 public:
  Matrix()
      : data(nullptr),
        stride(size_t(0)) {
  }
  ;
  inline Matrix(const MatrixShape<1> &s, const DataType n);
  inline Matrix(const MatrixShape<1> &s, const size_t st, const DataType n);
  inline Matrix(const Matrix<DataType, 1> &m);
  inline Matrix<DataType, 1> & operator =(const Matrix<DataType, 1> &m);
  inline Matrix(const SubMatrix<DataType, 1> &m);
  inline Matrix<DataType, 1> & operator =(const SubMatrix<DataType, 1> &m);
  inline Matrix(Matrix<DataType, 1> &&m);
  inline Matrix<DataType, 1> & operator =(Matrix<DataType, 1> &&m);
  inline Matrix(matrix_initializer_list<DataType, 1> t);
  inline Matrix & operator =(matrix_initializer_list<DataType, 1> t);
  ~Matrix() {
    delete[] data;
  }
  inline int get_capicity() const;
  inline int get_size() const;
  inline size_t get_length(const MatrixShape<1> &s, size_t stride) const;
  inline DataType * get_data() {
    return data;
  }
  inline DataType * get_data() const {
    return data;
  }
  inline const DataType & operator[](size_t i) const;
  inline DataType & operator[](size_t i);
  inline SubMatrix<DataType, 1> slice(size_t i, size_t j) const;
  inline SubMatrix<DataType, 1> slice(size_t i, size_t j);
  inline Matrix<DataType, 1> & operator +=(const DataType & n);
  inline Matrix<DataType, 1> & operator -=(const DataType & n);
  inline Matrix<DataType, 1> & operator *=(const DataType & n);
  inline Matrix<DataType, 1> & operator /=(const DataType & n);
  inline Matrix<DataType, 1> & operator +=(const Matrix<DataType, 1> & t);
  inline Matrix<DataType, 1> & operator -=(const Matrix<DataType, 1> & t);
  inline Matrix<DataType, 1> & operator *=(const Matrix<DataType, 1> & t);
  inline Matrix<DataType, 1>& operator /=(const Matrix<DataType, 1> & t);
  template<typename SubType>
  inline Matrix<DataType, 1>& operator=(const ExprBase<SubType, DataType> &e);
  inline DataType eval(size_t i, size_t j) const;
  inline void set_row_ele(int i, const Matrix<DataType, 1> & s);
  inline size_t getColumn() const {
    return column;
  }
  inline size_t getRow() const {
    return row;
  }
  inline size_t getCapicity() const {
    return capicity;
  }
  //private:
  MatrixShape<1> shape;
  size_t stride;
  size_t capicity;
  DataType * data;
  size_t row;
  size_t column;
};

//=============================================================================
//matrix operation for two dimensional Matrix and SubMatrix
// matrix multiplication operator

// overload equal Operator
// Ruturn true if all elements in the matrix is equal and stride is equal
template<typename DataType, size_t N>
inline bool operator==(const Matrix<DataType, N>& m1,
                       const Matrix<DataType, N>& m2) {
  DataType * d1 = m1.data;
  DataType * d2 = m2.data;
  if (m1.get_size() != m2.get_size())
    return false;
  for (int i = 0; i < m1.get_size(); ++i) {
    if (d1[i] != d2[i])
      return false;
  }
  return true;
}

template<typename DataType, size_t N>
inline bool operator==(const SubMatrix<DataType, N>& m1,
                       const Matrix<DataType, N>& m2) {
  DataType * d1 = m1.data;
  DataType * d2 = m2.data;
  if (m1.get_size() != m2.get_size())
    return false;
  for (int i = 0; i < m1.get_size(); ++i) {
    if (d1[i] != d2[i])
      return false;
  }
  return true;
}

template<typename DataType, size_t N>
inline bool operator==(const Matrix<DataType, N>& m1,
                       const SubMatrix<DataType, N>& m2) {
  DataType * d1 = m1.data;
  DataType * d2 = m2.data;
  if (m1.get_size() != m2.get_size())
    return false;
  for (int i = 0; i < m1.get_size(); ++i) {
    if (d1[i] != d2[i])
      return false;
  }
  return true;
}

template<typename DataType, size_t N>
inline bool operator==(const SubMatrix<DataType, N>& m1,
                       const SubMatrix<DataType, N>& m2) {
  DataType * d1 = m1.data;
  DataType * d2 = m2.data;
  if (m1.get_size() != m2.get_size())
    return false;
  for (int i = 0; i < m1.get_size(); ++i) {
    if (d1[i] != d2[i])
      return false;
  }
  return true;
}

}  //namespace matrix
#ifdef MATRIX_SCALAR_TYPE_
#error "MATRIX_SCALAR_TYPE_ Should not be defined!"
#endif
#define MATRIX_SCALAR_TYPE_ float
#include "expr-inl.h"
#undef MATRIX_SCALAR_TYPE_
#define MATRIX_SCALAR_TYPE_ double
#include "expr-inl.h"
#undef MATRIX_SCALAR_TYPE_
#define MATRIX_SCALAR_TYPE_ int
#include "expr-inl.h"
#undef MATRIX_SCALAR_TYPE_
#include "matrix-inl.h"
#include "sub_matrix-inl.h"
#include "matrix_math.h"
#include "random.h"
#endif // MATRIX_MATRIX_H

max_stars_repo_name: seaslee/snoopy

max_issues_repo_name: seaslee/snoopy

max_forks_repo_name: seaslee/snoopy


text:
#include <stdio.h>
#include "sweeny_uf.h"
#include <ctype.h>
#include <strings.h>
#include <stdio.h>
#include "../src/uf.h"
#include "../src/queue.h"
#include <math.h>
#include <gsl/gsl_rng.h> //Verwendung von Mersenne-Twister Pseudozufallszahlengenerator
static __s8 dN;

#define MIN(a,b) a <= b ? a: b
#define MAX(a,b) a >= b? a: b

static char setup=0;
static char verbose=0;
static double rcWeights[4]; // array with precalculated mc weights
static double p_min_del,p_max_del, p_min_ins,p_max_ins;
static struct queue collect1,collect2;
static struct queue_node *todo_pool, *collect_pool;
static __u32 offset_1=1,offset_2=2;
static __u64 *sec_cs_moment,*four_cs_moment;
static __u32 *size_giant;
static __u32 *num_bonds;
static __u32 *num_cluster;
static __u32 DX;
static __u32 N;
static __u32 seed;
static double q;
static double coupling;
static double beta;
static double v;
static double K;
static __u32 steps;
static __s8 *bonds;
static gsl_rng * r;
static __u32 edge[2];
static __s32 adj1[4];
static __s32 adj2[4];
static __s32 *uf1;
static __u32 cs1,cs2;
static __u32 activeEdges=0;
static __u32 cutoff;
/******************************************************************************
 *****************************************************************************/
static char init(void)
{
  K=coupling*beta;
  v = exp(K) - 1.0;
  N = DX*DX;
  rcWeights[0] = MIN(pow(v,-1),1); //db==-1,dN==0
  rcWeights[1] = MIN(pow(v,-1)*q,1); //db==-1,dN=1
  rcWeights[2] = MIN(v,1); //db==1,dN==0
  rcWeights[3] = MIN(v*pow(q,-1),1); //db==1,dN==-1
  p_min_del = MIN(rcWeights[0],rcWeights[1]);
  p_max_del = MAX(rcWeights[0],rcWeights[1]);
  p_min_ins = MIN(rcWeights[2],rcWeights[3]);
  p_max_ins = MAX(rcWeights[2],rcWeights[3]);
  if(verbose) {
	  printf("p_min_del = %f\n",p_min_del);
	  printf("p_max_del = %f\n",p_max_del);
	  printf("p_min_ins = %f\n",p_min_ins);
	  printf("p_max_ins = %f\n",p_max_ins);

  }
  r = gsl_rng_alloc(gsl_rng_mt19937);
  if(!r)
      return 0;
  gsl_rng_set(r,seed);
  __u32 i;
  uf1 = initUF(N);
  todo_pool = malloc(N*sizeof(struct queue_node));
  if(!todo_pool) {
      return 0;
  }
  collect_pool = malloc(N*sizeof(struct queue_node));
  if(!collect_pool) {
      return 0;
  }
  for(i=0;i<N;++i) {

      collect_pool[i].next = todo_pool[i].next = NULL;
      collect_pool[i].data = todo_pool[i].data= i;
      collect_pool[i].visited = todo_pool[i].visited = 0;
    }

  bonds = malloc(sizeof(*bonds)*2*N);
  if(!bonds) {
    return 0;
  }
  for(i=0;i<2*N;i++) bonds[i] = -1;
  return 1;
}
/******************************************************************************
 *****************************************************************************/
static void destroy(void)
{
  destroyUF(uf1);
  free(bonds);
  gsl_rng_free(r);
  free(collect_pool);
  free(todo_pool);
}
/******************************************************************************
 *****************************************************************************/
static inline __u32 ltcXnext(__u32 idx)
{
  if((idx+1) % DX)
    return idx+1;
  else
    return idx +1- DX ;
}
/******************************************************************************
 *****************************************************************************/
static inline __u32 ltcXprev(__u32 idx)
{
  if((idx%DX))
    return idx-1;
  else
    return idx + DX  -1;
}
/******************************************************************************
 *****************************************************************************/
static inline __u32 ltcYnext(__u32 idx)
{

  if(idx <= N -1 && idx >= N - DX)
    return idx +DX- N ;
  else
    return idx+DX;
}
/******************************************************************************
 *****************************************************************************/
static inline __u32 ltcYprev(__u32 idx)
{
  if(idx <= DX - 1)
    return idx + N - DX;
  else
    return idx - DX;

}
/******************************************************************************
 *****************************************************************************/
static void Adjacent(__u32 bidx)
{

  if(bidx%2) {
    bidx = (bidx -1)/2;
    edge[0] = bidx;
    edge[1] = ltcXnext(bidx);
  }
  else {
    bidx = bidx/2;
    edge[0] = bidx;
    edge[1] = ltcYprev(bidx);
  }

}
/******************************************************************************
 *****************************************************************************/
static void Adjacent2(__u32 idx, __u8 a) {
  if(a == 1)
  {
    if(bonds[2*idx] == 1)  adj1[0] = ltcYprev(idx);
    else  adj1[0] = -1;
    if(bonds[(2*idx)+1] == 1)  adj1[1] = ltcXnext(idx);
    else adj1[1] = -1;
    if(bonds[2*ltcYnext(idx)] == 1)	adj1[2] = ltcYnext(idx);
    else adj1[2] = -1;
    if(bonds[(2*ltcXprev(idx))+1] ==1)adj1[3] = ltcXprev(idx);
    else adj1[3] = -1;
  }
  else {
    if(bonds[2*idx] == 1)adj2[0] = ltcYprev(idx);
    else adj2[0] = -1;
    if(bonds[(2*idx)+1] == 1)adj2[1] = ltcXnext(idx);
    else adj2[1] = -1;
    if(bonds[2*ltcYnext(idx)] == 1)	adj2[2] = ltcYnext(idx);
    else adj2[2] = -1;
    if(bonds[(2*ltcXprev(idx))+1] == 1)adj2[3] = ltcXprev(idx);
    else adj2[3] = -1;
  }

}
/******************************************************************************
 *****************************************************************************/
static __u8 breadthFirstSearch(__u32 start, __u32 goal,__u8 accept_split)
{

    static struct queue todo1,todo2;
  __u32 i=0,activeP1=0,activeP2=0;
  init_queue(&todo1,todo_pool);
  init_queue(&todo2,todo_pool);
  init_queue(&collect1,collect_pool);
  init_queue(&collect2,collect_pool);
  cs1=0;
  cs2=0;
  enqueue(&todo1,start); // Put starting point onto queue 1
  enqueue(&collect1,start);
  cs1++; // Increase cluster size of bfs 1
  collect_pool[start].visited = offset_1;
  enqueue(&todo2,goal); // Put goal point onto queue 2
  enqueue(&collect2,goal);
  collect_pool[goal].visited = offset_2;
  cs2++; // increase cluster size of bfs 2
    while(!queue_empty_p(&todo1) && !queue_empty_p(&todo2)) {
        dequeue(&todo2,&activeP2); // get next vertex of bfs 2
        Adjacent2(activeP2,2);	 // get all adjacent vertices of current vertex
        for(i=0;i<4;i++) {
            if(adj2[i] != -1) { // if accessable (edge exists)
                if(collect_pool[adj2[i]].visited == offset_2) continue; // already visited
                if((__u32)adj2[i] == start || collect_pool[adj2[i]].visited == offset_1) { // reconnected
                    while(!queue_empty_p(&todo1)) dequeue(&todo1,&activeP1); // empty queue 1
                    while(!queue_empty_p(&todo2)) dequeue(&todo2,&activeP2); // empty queue 2
                    while(!queue_empty_p(&collect1)) dequeue(&collect1,&activeP1);
                    while(!queue_empty_p(&collect2)) dequeue(&collect2,&activeP2);
                    offset_1+=2;
                    offset_2+=2;
                    return 1; // 1 indicates success
                }
                enqueue(&todo2,adj2[i]);
                enqueue(&collect2,adj2[i]);
                collect_pool[adj2[i]].visited = offset_2; // mark as visited
                cs2++; // increase cluster size
            }
        }
        dequeue(&todo1,&activeP1);
        Adjacent2(activeP1,1);
        for(i=0;i<4;i++) {
            if(adj1[i] != -1) {
                if(collect_pool[adj1[i]].visited == offset_1) continue;
                if((__u32)adj1[i] == goal || collect_pool[adj1[i]].visited == offset_2) {
                    while(!queue_empty_p(&todo1)) dequeue(&todo1,&activeP1);
                    while(!queue_empty_p(&todo2)) dequeue(&todo2,&activeP2);
                    while(!queue_empty_p(&collect1)) dequeue(&collect1,&activeP1);
                    while(!queue_empty_p(&collect2)) dequeue(&collect2,&activeP2);
                    offset_1+=2;
                    offset_2+=2;
                    return 1;
                }
                enqueue(&todo1,adj1[i]);
                enqueue(&collect1,adj1[i]);
                collect_pool[adj1[i]].visited = offset_1;
                cs1++;
            }
        }
    }
    /* Only if a deletion of a cutting edge has to be accepted,
        * completely traverse the unfinished cluster to
        * be able to update the UF datastructure afterwards based on
        * the collect1 and collect2 queues.
    */
    if(accept_split) {
        while(!queue_empty_p(&todo1)) {
            dequeue(&todo1,&activeP1);
            Adjacent2(activeP1,1);
            for(i=0;i<4;i++) {
                if(adj1[i] != -1) {
                    if(collect_pool[adj1[i]].visited == offset_1) continue;
                    enqueue(&todo1,adj1[i]);
                    enqueue(&collect1,adj1[i]);
                    collect_pool[adj1[i]].visited = offset_1;
                    cs1++;
                }
            }

        }
        while(!queue_empty_p(&todo2)) {
            dequeue(&todo2,&activeP2);
            Adjacent2(activeP2,2);
            for(i=0;i<4;i++) {
                if(adj2[i] != -1) {
                    if(collect_pool[adj2[i]].visited == offset_2) continue;
                    enqueue(&todo2,adj2[i]);
                    enqueue(&collect2,adj2[i]);
                    collect_pool[adj2[i]].visited = offset_2;
                    cs2++;
                }
            }
        }
    }
    else {
        while(!queue_empty_p(&todo1)) dequeue(&todo1,&activeP1);
        while(!queue_empty_p(&todo2)) dequeue(&todo2,&activeP2);
        while(!queue_empty_p(&collect1)) dequeue(&collect1,&activeP1);
        while(!queue_empty_p(&collect2)) dequeue(&collect2,&activeP2);
    }
    offset_1+=2;
    offset_2+=2;

    return 0;
}
/******************************************************************************
 *****************************************************************************/
static __u8 removeBond(__u32 a, __u32 b,__u8 accept_split,__u8 accept_non_split) {
    __u32 vertex;
    if(breadthFirstSearch(a,b,accept_split)) //replacement edge found
        return accept_non_split;
    else {  // no replacement edge available and hence cluster will be splitted if move accepted
        if(accept_split) {  // accept cluster splitting case? if yes rebuild uf
            while(!queue_empty_p(&collect1)) {
                dequeue(&collect1,&vertex);
                uf1[vertex] = a;
            }
            uf1[a] = -1*cs1;
            while(!queue_empty_p(&collect2)) {
                dequeue(&collect2,&vertex);
                uf1[vertex] = b;
            }
            uf1[b] = -1*cs2;
            return 1;
        } // not accepted
        else {
            return 0;
        }
    }

}
/******************************************************************************
 *****************************************************************************/
static inline __u8 acceptChange(__s8 d_numCl, __s8 d_numAB, double randV)
{
  return d_numAB == -1 ? (randV < rcWeights[d_numCl]) : (randV < rcWeights[2 - d_numCl]);
}
/******************************************************************************
 *****************************************************************************/
static inline void   mcStep(void)
{

  static __u32 bidx;
  static double rnd_num;
  bidx = gsl_rng_uniform_int(r,2*N);
  rnd_num  = gsl_rng_uniform(r);
  Adjacent(bidx);
  if(bonds[bidx] == 1) { //edge is active hence delete it
    if(rnd_num < p_min_del) {
   	  bonds[bidx] = -1;
      if(removeBond(edge[0],edge[1],1,1))
        activeEdges--;
      else {
          if(verbose)fprintf(stderr,"ERROR!\n");
          exit(-1);
      }
    }
    else {
	if(rnd_num < p_max_del) {
	  bonds[bidx] = -1;
	  if(removeBond(edge[0],edge[1], acceptChange(1,-1,rnd_num),acceptChange(0,-1,rnd_num)))
          activeEdges--;
      else
          bonds[bidx] = 1;
    }

    }

  }
  else {
    if(rnd_num < p_min_ins) {
        bonds[bidx]  = 1;
        activeEdges++;
        unite(edge[0],edge[1],uf1);
      }
    else {
    if(rnd_num < p_max_ins) {
      dN = connected(edge[0],edge[1],uf1) ? 0: -1;
      if(acceptChange(dN,1,rnd_num)) {
        bonds[bidx]=1;
        activeEdges++;
        if(dN == -1)
            unite(edge[0],edge[1],uf1);
      }
    }

  }
  }
}
/******************************************************************************
 *****************************************************************************/
static void extract_observables(__u32 cnt)
{
  __u32 i=0,tn=0,nclust=0;
  size_giant[cnt] = 0;
  sec_cs_moment[cnt] = 0;
  four_cs_moment[cnt] = 0;
  for(i=0;i<N && tn <= N;i++) {
    if(uf1[i] < 0) {
      tn -= uf1[i];
      nclust+=1;
      sec_cs_moment[cnt]+= (__u64)pow(uf1[i],2);
      four_cs_moment[cnt]+= (__u64)pow(uf1[i],4);
      if((__u32)(-uf1[i]) > size_giant[cnt]) size_giant[cnt] = (__u32)(-uf1[i]);
    }
  }
  num_bonds[cnt] = activeEdges;
  num_cluster[cnt] = nclust;
}
/******************************************************************************
 *****************************************************************************/
static void generateTimeSeries(void)
{
  __u32 i=0,j=0;
  for(i=0;i<cutoff;i++){
      for(j=0;j<2*N;j++)mcStep();
      // reset visited array
      // otherwise it will overflow at some point
      // could also be done less frequent
      for(j=0;j<N;j++)collect_pool[j].visited=0;
  }
  for(i=0;i<steps;i++) {
    for(j=0;j<2*N;j++)
      mcStep();
      extract_observables(i);
      for(j=0;j<N;j++)collect_pool[j].visited=0;
    }
}
/******************************************************************************
 *****************************************************************************/
char init_sweeny_uf(double _q,unsigned int _l,double _beta,double _coupl,
        unsigned int _cutoff,unsigned _tslength,unsigned int rng_seed,
        void *ts_0,void *ts_1,void * ts_2,void *ts_3, void *ts_4) {
    q = _q;
    DX = _l;
    beta = _beta;
    coupling = _coupl;
    cutoff = _cutoff;
    steps = _tslength;
    seed = rng_seed;
    num_bonds = (__u32 *)ts_0;
    num_cluster = (__u32 *)ts_1;
    size_giant = (__u32 *)ts_2;
    sec_cs_moment = (__u64 *)ts_3;
    four_cs_moment = (__u64 *)ts_4;
    return setup=init();

}
/******************************************************************************
 *****************************************************************************/
void destroy_sweeny_uf(void) {
    if(setup)
        destroy();
    setup=0;
}
/******************************************************************************
 *****************************************************************************/
char simulate_sweeny_uf(void) {
    if(!setup)return 0;
    generateTimeSeries();
    return 1;
}
/******************************************************************************
 *****************************************************************************/

max_issues_repo_name: ernmeel/sweeny


text:
#pragma once

#include <array>

#include <gsl/span>

#include "Quiver/Animation/Rect.h"

namespace qvr {

struct ViewBuffer {
	int viewCount = 0;
	std::array<Animation::Rect, 8> views;

	// TODO: Consider making initial value of viewCount == 1.
	// TODO: Consider making viewCount < 1 illegal.
};

// TODO: Consider making these public member functions of ViewBuffer and
// making the data members private.

inline void SetView(
	ViewBuffer& vb,
	const Animation::Rect& singleView)
{
	vb.viewCount = 1;
	vb.views[0] = singleView;
}

inline void SetViews(
	ViewBuffer& target,
	const gsl::span<const Animation::Rect> newViews)
{
	target.viewCount =
		std::min(
		(int)newViews.length(),
			(int)target.views.max_size());

	std::copy(
		std::begin(newViews),
		std::begin(newViews) + target.viewCount,
		std::begin(target.views));
}

inline auto GetViews(const ViewBuffer& vb) -> gsl::span<const Animation::Rect> {
	return gsl::make_span(&vb.views[0], vb.viewCount);
}

inline const Animation::Rect& CalculateView(
	const ViewBuffer& vb,
	const float objectAngle,
	const float viewAngle)
{
	assert(vb.viewCount > 0);

	const float Tau = 6.28318530718f;

	const int viewIndex =
		(int)((fmodf((objectAngle - viewAngle + Tau), Tau) / Tau) * (vb.viewCount));

	// Make sure my maths is right.
	assert(viewIndex >= 0);
	assert(viewIndex < vb.viewCount);

	return vb.views[viewIndex];
}

}

max_stars_repo_name: rachelnertia/Quarrel

max_issues_repo_name: rachelnertia/Quarrel

max_forks_repo_name: rachelnertia/Quiver


text:
/* ///////////////////////////////////////////////////////////////////// */
/*!
  \file
  \brief Contains basic functions for problem initialization.

  The init.c file collects most of the user-supplied functions useful
  for problem configuration.
  It is automatically searched for by the makefile.

  \author A. Mignone (mignone@ph.unito.it)
  \date   March 5, 2017
*/
/* ///////////////////////////////////////////////////////////////////// */
#include "pluto.h"
#include "local_pluto.h"
#include <math.h>
#include "coolit.h"
#include "init_interp.h"
#include <stdbool.h>
//#include <gsl.h>

double mu, mue, mui, Rcl;
double rho0, P0, cs0, v0, s0, cool_0, tcool0, r_crit;
double rhobnd[2], pbnd[2], vbnd[2]; //boundary conditions at inner and outer radius

int q = 1, counter = 0;
double T0 , Mach0, n0, chi, singular;

/* ********************************************************************* */
void Init (double *v, double x1, double x2, double x3)
/*!
 * The Init() function can be used to assign initial conditions as
 * as a function of spatial position.
 *
 * \param [out] v   a pointer to a vector of primitive variables
 * \param [in] x1   coordinate point in the 1st dimension
 * \param [in] x2   coordinate point in the 2nd dimension
 * \param [in] x3   coordinate point in the 3rdt dimension
 *
 * The meaning of x1, x2 and x3 depends on the geometry:
 * \f[ \begin{array}{cccl}
 *    x_1  & x_2    & x_3  & \mathrm{Geometry}    \\ \noalign{\medskip}
 *     \hline
 *    x    &   y    &  z   & \mathrm{Cartesian}   \\ \noalign{\medskip}
 *    R    &   z    &  -   & \mathrm{cylindrical} \\ \noalign{\medskip}
 *    R    & \phi   &  z   & \mathrm{polar}       \\ \noalign{\medskip}
 *    r    & \theta & \phi & \mathrm{spherical}
 *    \end{array}
 *  \f]
 *
 * Variable names are accessed by means of an index v[nv], where
 * nv = RHO is density, nv = PRS is pressure, nv = (iVR, VX2, VX3) are
 * the three components of velocity, and so forth.
 *
 *********************************************************************** */
{

  #if PHYSICS == MHD || PHYSICS == RMHD
  v[BX1] = 0.0;
  v[BX2] = 0.0;
  v[BX3] = 0.0;

  v[AX1] = 0.0;
  v[AX2] = 0.0;
  v[AX3] = 0.0;
  #endif
}

/* ********************************************************************* */
void InitDomain (Data *d, Grid *grid)
/*!
 * Assign initial condition by looping over the computational domain.
 * Called after the usual Init() function to assign initial conditions
 * on primitive variables.
 * Value assigned here will overwrite those prescribed during Init().
 *
 *
 *********************************************************************** */
{
	int i, j, k;
	int id;
	double *x1 = grid->x[IDIR];
	double *x2 = grid->x[JDIR];
	double *x3 = grid->x[KDIR];
	double ne0, ni0;
	double* muemui = MeanElIonWeight((double*)d->Vc);
    mue = muemui[0];
	mui = muemui[1];
	mu = MeanMolecularWeight((double*)d->Vc);

	T0 = g_inputParam[T_CRIT]; //K
	Mach0 = g_inputParam[MACH_CRIT];
	n0 = g_inputParam[N_CRIT]; //cm^-3
	chi = g_inputParam[CHI];
	singular = g_inputParam[SINGULAR];

	rho0 = n0*mu*CONST_mp; //g cm^-3
	ne0 = n0*mu/mue;
	ni0 = n0*mu/mui;
	P0 = n0*CONST_kB*T0; //dyne cm^-2
	cs0 = sqrt(g_gamma*P0/rho0); //cm s^-1
	v0 = -Mach0*cs0; //cm s^-1
	s0 = P0/pow(rho0,g_gamma); //CGS
	cool_0 = lambda(T0); //CGS
	tcool0 = (1./(g_gamma-1))*n0*CONST_kB*T0/(ne0*ni0*cool_0); //s

	r_crit = -q*v0*g_gamma*tcool0; //cm
	Rcl = singular*r_crit*pow(UNIT_LENGTH,-1); // code units (pc usually)

	//print("Test: %.3e\t%.3e\t%.3e\n",r_crit,UNIT_LENGTH,r_crit*pow(UNIT_LENGTH,-1));
	print("T0=%.1e\tM0=%.1f\tn0=%.2e\t\tchi=%d\tsingular=%.2f\n",T0,Mach0, n0, (int)chi, singular);
	print("rho0=%.1e\t\tP0=%.1e\t\tv0=%.2e\t\tr_crit=%.2e kpc\tRcl=%.2e kpc\t\ttcool0=%.2e Myr\n",rho0/UNIT_DENSITY,P0/(UNIT_DENSITY*pow(UNIT_VELOCITY,2)), v0/UNIT_VELOCITY, r_crit*pow(UNIT_LENGTH,-1), Rcl, tcool0*pow(CONST_PI*1e13,-1.));
	print("Starting from %.2f kpc (r/ro=%.2e) to %.2f kpc (r/r0=%.2e)\n",g_domBeg[IDIR],g_domBeg[IDIR]/r_crit*UNIT_LENGTH,g_domEnd[IDIR],g_domEnd[IDIR]/r_crit*UNIT_LENGTH);
	print("Using Solar metallicity\n");
	print("Useful info: mu=%.2f mu_e=%.2f mu_i=%.2f\n",mu,mue,mui);
	print("Mass flux: %.2e Msun/yr\n",-4*CONST_PI*r_crit*r_crit*rho0*v0/(CONST_Msun/(365*24*60*60)));

	#if SUBSONIC_INI != NO
	readICfile("subsonic_rdpv.txt");
	#else
	readICfile("transonic_rdpv.txt");
	#endif

	TOT_LOOP(k,j,i){
		//if (prank==0 && (x1[i]<g_domBeg[IDIR] || x1[i]>g_domEnd[IDIR])) printf("%.2f pc\n",x1[i]);
		double* dpv = initdpv(x1[i]*UNIT_LENGTH/r_crit);
		if (i==0){
			rhobnd[0] = dpv[0]; pbnd[0] = dpv[1]; vbnd[0] = dpv[2];
		}
		else if(i==NX1-10){
			rhobnd[1] = dpv[0]; pbnd[1] = dpv[1]; vbnd[1] = dpv[2];
		}
		//print("%.2e\t\t%.2e\t\t%.2e\t\t%.2e\n",x1[i]*UNIT_LENGTH/r_crit,dpv[0],dpv[1],dpv[2]);
		d->Vc[RHO][k][j][i] = dpv[0]*rho0/UNIT_DENSITY; //(x1[i]/Rcl<=1.)?chi*n0*mu*CONST_mp/UNIT_DENSITY:n0*mu*CONST_mp/UNIT_DENSITY;
		d->Vc[PRS][k][j][i] = dpv[1]*P0/(UNIT_DENSITY*pow(UNIT_VELOCITY,2));
		EXPAND(d->Vc[iVR][k][j][i] = dpv[2]*abs(v0)/UNIT_VELOCITY; ,d->Vc[iVZ][k][j][i] = 0.;, d->Vc[VX3][k][j][i] = 0.;)
		d->Vc[TRC][k][j][i] = 1.0;
	}
	print("Done reading IC from file!\n");
	print("Left fixed boundary: rho/rho0=%.2e\t\tp/p0=%.2e\t\tv/v0=%.2e\n",rhobnd[0],pbnd[0],vbnd[0]);
	print("Right fixed boundary: rho/rho0=%.2e\t\tp/p0=%.2e\t\tv/v0=%.2e\n",rhobnd[1],pbnd[1],vbnd[1]);
	#if VX1RHO_OUTFLOW ==YES
	print("iVR and Density are outflowing.\n");
	#endif
}

/* ********************************************************************* */
void Analysis (const Data *d, Grid *grid)
/*!
 *  Perform runtime data analysis.
 *
 * \param [in] d the PLUTO Data structure
 * \param [in] grid   pointer to array of Grid structures
 *
 *********************************************************************** */
{
	/*
	int i, j, k;
	double *dx, *dy, *dz;*/
	/* ---- Parallel data reduction ---- */
	/*
	#ifdef PARALLEL
	int transfer_size = 3;
	int transfer = 0;
	double sendArray [transfer_size], recvArray[transfer_size];*/

	/* ---- Set pointer shortcuts ---- */
	/*
	dx = grid->dx[IDIR];
	dy = grid->dx[JDIR];
	dz = grid->dx[KDIR];*/

}
#if PHYSICS == MHD
/* ********************************************************************* */
void BackgroundField (double x1, double x2, double x3, double *B0)
/*!
 * Define the component of a static, curl-free background
 * magnetic field.
 *
 * \param [in] x1  position in the 1st coordinate direction \f$x_1\f$
 * \param [in] x2  position in the 2nd coordinate direction \f$x_2\f$
 * \param [in] x3  position in the 3rd coordinate direction \f$x_3\f$
 * \param [out] B0 array containing the vector componens of the background
 *                 magnetic field
 *********************************************************************** */
{
   B0[0] = 0.0;
   B0[1] = 0.0;
   B0[2] = 0.0;
}
#endif

/* ********************************************************************* */
void UserDefBoundary (const Data *d, RBox *box, int side, Grid *grid)
/*!
 *  Assign user-defined boundary conditions.
 *
 * \param [in,out] d  pointer to the PLUTO data structure containing
 *                    cell-centered primitive quantities (d->Vc) and
 *                    staggered magnetic fields (d->Vs, when used) to
 *                    be filled.
 * \param [in] box    pointer to a RBox structure containing the lower
 *                    and upper indices of the ghost zone-centers/nodes
 *                    or edges at which data values should be assigned.
 * \param [in] side   specifies the boundary side where ghost zones need
 *                    to be filled. It can assume the following
 *                    pre-definite values: X1_BEG, X1_END,
 *                                         X2_BEG, X2_END,
 *                                         X3_BEG, X3_END.
 *                    The special value side == 0 is used to control
 *                    a region inside the computational domain.
 * \param [in] grid  pointer to an array of Grid structures.
 *
 *********************************************************************** */
{
  int   i, j, k, nv;
  double  *x1, *x2, *x3;

  x1 = grid->x[IDIR];
  x2 = grid->x[JDIR];
  x3 = grid->x[KDIR];

  if (side == X1_BEG) {

    BOX_LOOP(box,k,j,i) {
		d->Vc[PRS][k][j][i] = pbnd[0]*P0/(UNIT_DENSITY*pow(UNIT_VELOCITY,2)); //only fix the pressure
		#if VX1RHO_OUTFLOW ==YES
		d->Vc[RHO][k][j][i] = d->Vc[RHO][k][j][IBEG];   //outflow density
		d->Vc[iVR][k][j][i] = d->Vc[iVR][k][j][IBEG];   //outflow velocity
		#else
		d->Vc[RHO][k][j][i] = rhobnd[0]*rho0/UNIT_DENSITY;
		d->Vc[iVR][k][j][i] = vbnd[0]*abs(v0)/UNIT_VELOCITY;
		#endif
	}

  }else if (side == X1_END){

    BOX_LOOP(box,k,j,i) {
		d->Vc[PRS][k][j][i] = pbnd[1]*P0/(UNIT_DENSITY*pow(UNIT_VELOCITY,2));
		#if VX1RHO_OUTFLOW ==YES
		d->Vc[RHO][k][j][i] = d->Vc[RHO][k][j][IEND];   //outflow density
		d->Vc[iVR][k][j][i] = d->Vc[iVR][k][j][IEND]; //outflow velocity
		#else
		d->Vc[RHO][k][j][i] = rhobnd[1]*rho0/UNIT_DENSITY;
		d->Vc[iVR][k][j][i] = vbnd[1]*abs(v0)/UNIT_VELOCITY;
		#endif
	}

  }

  /*
  if (side == 0) {  */  /* -- check solution inside domain -- */
  /* TOT_LOOP(k,j,i){

    if ((x1[i]*UNIT_LENGTH/r_crit>=(1.-1e-1)) && (x1[i]*UNIT_LENGTH/r_crit<=(1.+1e-1))) {
		d->Vc[RHO][k][j][i] = rho0/UNIT_DENSITY;
		d->Vc[iVR][k][j][i] = v0/UNIT_VELOCITY;
		d->Vc[PRS][k][j][i] = P0/(UNIT_DENSITY*pow(UNIT_VELOCITY,2));
	}/*
	if (x1[i]/Rcl<1.)
		d->Vc[iVR][k][j][i] = 0.; *//*
   }
  }*/
}

#if BODY_FORCE != NO
/* ********************************************************************* */
void BodyForceVector(double *v, double *g, double x1, double x2, double x3)
/*!
 * Prescribe the acceleration vector as a function of the coordinates
 * and the vector of primitive variables *v.
 *
 * \param [in] v  pointer to a cell-centered vector of primitive
 *                variables
 * \param [out] g acceleration vector
 * \param [in] x1  position in the 1st coordinate direction \f$x_1\f$
 * \param [in] x2  position in the 2nd coordinate direction \f$x_2\f$
 * \param [in] x3  position in the 3rd coordinate direction \f$x_3\f$
 *
 *********************************************************************** */
{
  g[IDIR] = 0.0;
  g[JDIR] = 0.0;
  g[KDIR] = 0.0;
}
/* ********************************************************************* */
double BodyForcePotential(double x1, double x2, double x3)
/*!
 * Return the gravitational potential as function of the coordinates.
 *
 * \param [in] x1  position in the 1st coordinate direction \f$x_1\f$
 * \param [in] x2  position in the 2nd coordinate direction \f$x_2\f$
 * \param [in] x3  position in the 3rd coordinate direction \f$x_3\f$
 *
 * \return The body force potential \f$ \Phi(x_1,x_2,x_3) \f$.
 *
 *********************************************************************** */
{
  return 0.0;
}
#endif

max_stars_repo_name: dutta-alankar/cooling-flow-model

max_issues_repo_name: dutta-alankar/cooling-flow-model

max_forks_repo_name: dutta-alankar/cooling-flow-model


text:
//! \file neutronics_driver.h
//! Base class for single-physics neutronics solver
#ifndef NEUTRONICS_DRIVER_H
#define NEUTRONICS_DRIVER_H

#include "enrico/driver.h"
#include "enrico/geom.h"
#include "enrico/mpi_types.h"

#include <gsl/gsl>
#include <xtensor/xtensor.hpp>

#include <vector>

namespace enrico {

using CellHandle = gsl::index;

//! Base class for driver that controls a neutronics solve
class NeutronicsDriver : public Driver {
public:
  explicit NeutronicsDriver(MPI_Comm comm)
    : Driver(comm)
  {}

  virtual ~NeutronicsDriver() = default;

  //! Get energy deposition in each material normalized to a given power
  //! \param power User-specified power in [W]
  //! \return Heat source in each material as [W/cm3]
  virtual xt::xtensor<double, 1> heat_source(double power) const = 0;

  //! Find cells corresponding to a vector of positions
  //! \param positions (x,y,z) coordinates to search for
  //! \return Handles to cells
  virtual std::vector<CellHandle> find(const std::vector<Position>& positions) = 0;

  //! Set the density of the material in a cell
  //! \param cell Handle to a cell
  //! \param rho Density in [g/cm^3]
  virtual void set_density(CellHandle cell, double rho) const = 0;

  //! Set the temperature of a cell
  //! \param cell Handle to a cell
  //! \param T Temperature in [K]
  virtual void set_temperature(CellHandle cell, double T) const = 0;

  //! Get the density of a cell
  //! \param cell Handle to a cell
  //! \return Cell density in [g/cm^3]
  virtual double get_density(CellHandle cell) const = 0;

  //! Get the temperature of a cell
  //! \param cell Handle to a cell
  //! \return Temperature in [K]
  virtual double get_temperature(CellHandle cell) const = 0;

  //! Get the volume of a cell
  //! \param cell Handle to a cell
  //! \return Volume in [cm^3]
  virtual double get_volume(CellHandle cell) const = 0;

  //! Detemrine whether a cell contains fissionable nuclides
  //! \param cell Handle to a cell
  //! \return Whether the cell contains fissionable nuclides
  virtual bool is_fissionable(CellHandle cell) const = 0;

  //! Determine number of cells participating in coupling
  //! \return Number of cells
  virtual std::size_t n_cells() const = 0;

  //! Create energy production tallies
  virtual void create_tallies() = 0;

  //! Get a label for a cell
  //! \param cell Handle to a clel
  //! \return Label for the cell
  virtual std::string cell_label(CellHandle cell) const = 0;
};

} // namespace enrico

#endif // NEUTRONICS_DRIVER_H

max_stars_repo_name: pshriwise/enrico

max_issues_repo_name: pshriwise/enrico

max_forks_repo_name: lebuller/enrico


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_randist.h>
#include <errno.h>
#include <string.h>
#include <time.h>

const int SAMPLENUM = 1000;
int obsnum;
gsl_rng * r;

typedef struct {
	int prey;
	int predator;
} lvstate;

typedef struct {
	double time;
	double rawdata;
} obsdata;

void simPrior(lvstate *simData);
void rowsample(int *rows, double *w);
void stepLV(lvstate *state, double *t0p, double *dtp, double *lvParam);
void peturb(double *lvParam);

double obsLik(lvstate *mystate, double obs) {
	const double SIGMA = 10.0;
	return log(gsl_ran_gaussian_pdf((obs - mystate->prey), SIGMA));
}

void callStepLV(lvstate *simData, double *t0p, double *dtp, double *lvParam) {
	int j;
	for (j=0;j<SAMPLENUM;j++) {
		stepLV(simData + j, t0p, dtp, lvParam);
	}

}

//particle filter: output new estimated lvstate and its likehood
void pfPropPar(lvstate ppstate[16], obsdata *myobsdata, double *lvParam,
		double *ll) {
	//printf("Starting pfPropPar...\n");
	typedef struct {
		lvstate *vpptuple;
	} stateVec;

	int i, j, k, l;

	double timeUnit = 2.0;
	int deltaTimeUnit = 1;
	double delTime = timeUnit * deltaTimeUnit;
	double curTime = 0.0;

	double lw[SAMPLENUM];
	double w[SAMPLENUM];
	double wSum;
	double maxLw;
	double likehood = 0.0;
	unsigned int rows[SAMPLENUM];

	//printf("Starting pfPropPar 1 ...\n");
	stateVec *myStateVec = malloc(obsnum * sizeof(stateVec));
	if (myStateVec == NULL) {
		printf("Fail to initiate myStateVec...Exit");
		goto freeMem;
	}
	memset(myStateVec, 0, obsnum * sizeof(stateVec));

	stateVec *tempStateVec = malloc(obsnum * sizeof(stateVec));
	if (tempStateVec == NULL) {
		printf("Fail to initiate tempStateVec...Exit");
		goto freeMem;
	}
	memset(tempStateVec, 0, obsnum * sizeof(stateVec));

	lvstate *simData = malloc(SAMPLENUM * sizeof(lvstate));
	if (simData == NULL) {
		printf("Fail to initiate simData...EXIT!");
		goto freeMem;
	}
	memset(simData, 0, SAMPLENUM * sizeof(lvstate));
	//initiate state: get simulated prey-predator tuple
	simPrior(simData);

	*ll = 0.0;

	if ((tempStateVec + 0)->vpptuple == NULL) {
		(tempStateVec + 0)->vpptuple = (lvstate*) malloc(
				SAMPLENUM * sizeof(lvstate));
		if ((tempStateVec + 0)->vpptuple == NULL) {
			printf("Fail to initiate myStateVec->vpptuple...Exit");
			goto freeMem;
			//free(tempStateVec);
			//goto error;
		}
	}

	memset(tempStateVec[0].vpptuple, 0, SAMPLENUM * sizeof(lvstate));
	for (j = 0; j < SAMPLENUM; j++) {
		((tempStateVec + 0)->vpptuple + j)->predator = (simData + j)->predator;
		((tempStateVec + 0)->vpptuple + j)->prey = (simData + j)->prey;
	}

	for (i = 1; i < obsnum; i++) {
		curTime = myobsdata[i - 1].time;

		wSum = 0.0;

		if ((tempStateVec + i)->vpptuple == NULL) {
			//printf("init temp \n");
			(tempStateVec + i)->vpptuple = (lvstate*) malloc(
					SAMPLENUM * sizeof(lvstate));
			if ((tempStateVec + i)->vpptuple == NULL) {
				printf("Fail to initiate myStateVec->vpptuple...Exit");
				goto freeMem;
			}
			memset((tempStateVec + i)->vpptuple, 0,
					SAMPLENUM * sizeof(lvstate));
		}

		//printf("tempvec 1 %i\n", ((tempStateVec+0)->vpptuple+1)->prey);
		callStepLV(simData, &curTime, &delTime, lvParam);

		for (j = 0; j < SAMPLENUM; j++) {
			//stepLV(simData + j, &curTime, &delTime, lvParam);
			((tempStateVec + i)->vpptuple + j)->predator =
					(simData + j)->predator;
			((tempStateVec + i)->vpptuple + j)->prey = (simData + j)->prey;
			//printf("tempState %i vpptuple %i , %p\n", i,j,(tempStateVec+i)->vpptuple+j);
			lw[j] = obsLik(((tempStateVec + i)->vpptuple + j),
					myobsdata[i].rawdata);

			if (j == 0) {
				maxLw = lw[0];
			} else {
				if (lw[j] > maxLw) {
					maxLw = lw[j];
				}
			}
		}

		for (l = 0; l < SAMPLENUM; l++) {
			w[l] = exp(lw[l] - maxLw);
			wSum = wSum + w[l];
		}
		rowsample(rows, w);

		//reorganize the vector with only keeping the rows in 'rows'
		for (k = 0; k <= i; k++) {
			if ((myStateVec + k)->vpptuple == NULL) {
				//printf("init mystatevec  \n");
				(myStateVec + k)->vpptuple = (lvstate*) malloc(
						SAMPLENUM * sizeof(lvstate));
				if ((myStateVec + k)->vpptuple == NULL) {
					printf("Fail to initiate myStateVec->vpptuple...Exit");
					goto freeMem;
				}
			}
//			memset(myStateVec[k].vpptuple,0,SAMPLENUM * sizeof(lvstate));

			for (j = 0; j < SAMPLENUM; j++) {
				((myStateVec + k)->vpptuple + j)->predator =
						((tempStateVec + k)->vpptuple + rows[j])->predator;
				((myStateVec + k)->vpptuple + j)->prey =
						((tempStateVec + k)->vpptuple + rows[j])->prey;
				//printf("row %i %i,%i,%i,%i,%i\n", i,k,j,rows[j],((myStateVec+k)->vpptuple+j)->prey,((tempStateVec+k)->vpptuple+rows[j])->prey);
			}
			for (j = 0; j < SAMPLENUM; j++) {
				((tempStateVec + k)->vpptuple + j)->predator =
						((myStateVec + k)->vpptuple + j)->predator;
				((tempStateVec + k)->vpptuple + j)->prey =
						((myStateVec + k)->vpptuple + j)->prey;
			}
		}
		for (j = 0; j < SAMPLENUM; j++) {
			simData[j].predator = ((myStateVec + i)->vpptuple + j)->predator;
			simData[j].prey = ((myStateVec + i)->vpptuple + j)->prey;
		}
		likehood = likehood + maxLw + log(wSum / SAMPLENUM);
		*ll = likehood;
		rows[SAMPLENUM] = NULL;
	}

	for (i = 0; i < obsnum; i++) {
		ppstate[i].prey = ((myStateVec + i)->vpptuple + 0)->prey;
		ppstate[i].predator = ((myStateVec + i)->vpptuple + 0)->predator;
	}

	freeMem:
	//printf("free mem...");
	if (myStateVec != NULL) {
		//printf("free mem 2...");
		for (i = 0; ((myStateVec + i)->vpptuple != NULL) && i < obsnum; i++) {
			//printf("free myStateVec...i %i\n", i);
			free((myStateVec + i)->vpptuple);
			(myStateVec + i)->vpptuple = NULL;
		}
		free(myStateVec);
		myStateVec = NULL;
	}
	if (tempStateVec != NULL) {
		for (i = 0; NULL != (tempStateVec + i)->vpptuple && i < obsnum; i++) {
			//printf("free tempStateVec...i %i\n", i);
			free((tempStateVec + i)->vpptuple);
			(tempStateVec + i)->vpptuple = NULL;
		}
		free(tempStateVec);
		tempStateVec = NULL;
	}
	if (simData != NULL) {
		free(simData);
		simData = NULL;
	}

}

void stepLV(lvstate *state, double *t0p, double *dtp, double *lvParam) {
	//printf("Starting stepLV...\n");
	register double t0 = *t0p, dt = *dtp, t;
	register double h0, h1, h2, h3, u, l0, l1, l2;
	register int prey, predator;
	prey = state->prey;
	predator = state->predator;
	l0 = lvParam[0];
	l1 = lvParam[1];
	l2 = lvParam[2];

	while (dt > 0) {
		h1 = l0 * prey;
		h2 = l1 * prey * predator;
		h3 = l2 * predator;
		h0 = h1 + h2 + h3;

		if ((h0 < (1e-10)) || (prey >= 1000000))
			t = 1e99;
		else {
			t = gsl_ran_exponential(r, 1/h0);
		}
		if (t > dt) {
			state->prey = prey;
			state->predator = predator;
			return;
		} else {
			u = gsl_rng_uniform(r);
			if (u < (h1 / h0)) {
				prey = prey + 1;
				t0 = t0 + t;
				dt = dt - t;
			} else if (u < ((h1 + h2) / h0)) {
				prey = prey - 1;
				predator = predator + 1;
				t0 = t0 + t;
				dt = dt - t;
			} else {
				predator = predator - 1;
				t0 = t0 + t;
				dt = dt - t;
			}
		}
	}

	return;
}

//mcmc process
void runPmmhPath(int its, double *lvParam, double *obslik, lvstate ppstate[16],
		obsdata *myobsdata) {
	//printf("Starting runPmmhPath...\n");
	int i, j;
	double propParam[3];
	double curParam[3];
	double ll;
	double propMll, curMll = -1e99;
	lvstate curPath[obsnum];
	lvstate propPath[obsnum];

	FILE *fp = fopen("mcmc-out.txt", "wab+");
	if (fp == NULL) {
		printf("file is null.\n");
		return;
	}
	fprintf(fp, "th1,th2,th3,");
	for (i = 0; i <= 30; i++) {
		if ((i % 2) == 0) {
			if (i == 30) {
				fprintf(fp, "x%i,y%i", i, i);
			} else {
				fprintf(fp, "x%i,y%i,", i, i);
			}
		}
	}
	fprintf(fp, "\n");

	memcpy(curPath, ppstate, sizeof(lvstate) * obsnum);
	memcpy(curParam, lvParam, sizeof(double) * 3);
	for (i = its; i > 0; i--) {
		//calculate measurements and decide moving to next state or not
		memcpy(propParam, curParam, sizeof(double) * 3);
		peturb(propParam);

		pfPropPar(propPath, myobsdata, propParam, &ll);
		propMll = ll;
		if (log(gsl_ran_flat(r, 0.0, 1.0)) < (propMll - curMll)) {
			curMll = propMll;
			memcpy(curParam, propParam, sizeof(double) * 3);
			memcpy(curPath, propPath, sizeof(lvstate) * obsnum);
		}
		//write current parameters and state sequence into output
		fprintf(fp, "%f,%f,%f,", curParam[0], curParam[1], curParam[2]);

		for (j = 0; j < obsnum; j++) {
			if (j == (obsnum - 1)) {
				fprintf(fp, "%i,%i", (int)curPath[j].prey, (int)curPath[j].predator);
			} else {
				fprintf(fp, "%i,%i,", (int)curPath[j].prey, (int)curPath[j].predator);
			}
		}
		fprintf(fp, "\n");
	}
	fclose(fp);
}

void rowsample(int *rows, double *w) {
	//printf("Starting rowsample...\n");
	gsl_ran_discrete_t * grdp;
	int row;
	int i;
	grdp = gsl_ran_discrete_preproc(SAMPLENUM, w);
	for (i = 0; i < SAMPLENUM; i++) {

		row = (int) gsl_ran_discrete(r, grdp);
		rows[i] = row;

	}
	gsl_ran_discrete_free(grdp);
	return;
}

void peturb(double *lvParam) {
	//printf("Starting peturb...\n");

	const double SIGMA = 0.01;
	//printf("lvParam_before %f, %f, %f\n", lvParam[0], lvParam[1], lvParam[2]);
	lvParam[0] = lvParam[0] * exp(gsl_ran_gaussian(r, SIGMA));
	lvParam[1] = lvParam[1] * exp(gsl_ran_gaussian(r, SIGMA));
	lvParam[2] = lvParam[2] * exp(gsl_ran_gaussian(r, SIGMA));
	//printf("lvParam_after  %f, %f, %f\n", lvParam[0], lvParam[1], lvParam[2]);
}

//to get the prior simulate value of prey and predator
void simPrior(lvstate *simData) {
	//printf("Starting simPrior...\n");
	int i;

	const double PREY_MEAN = 50.0;
	const double PRED_MEAN = 100.0;

	for (i = 0; i < SAMPLENUM; i++) {
		simData[i].prey = gsl_ran_poisson(r, PREY_MEAN);
		simData[i].predator = gsl_ran_poisson(r, PRED_MEAN);
		//printf("simPrior %i\n", simData[i].prey);
	}
	//printf("End simPrior...\n");

}

void runModel(int its) {
	obsnum = 16;
	obsdata myobsdata[obsnum];
	lvstate mylvstate[obsnum];
	double ll = 0.0;

	double lvParam[3] = { 1.0, 0.005, 0.6 };

	//produce time list(with even number in)
	int i;
	int j = 0;

	for (i = 0; i <= 30; i++) {
		if ((i % 2) == 0) {
			myobsdata[j].time = (double) i;
			j++;
		}
	}

	//read external txt file
	FILE *file = fopen("LVpreyNoise10.txt", "r");
	char line[1024];
	if (file == NULL) {
		printf("file is null.\n");
		exit(EXIT_FAILURE);
	}
	if (file != NULL) {
		i = 0;
		while (fgets(line, 1024, file)) {
			myobsdata[i].rawdata = atof(line);
			i++;
		}
		fclose(file);
	}
	runPmmhPath(its, lvParam, &ll, mylvstate, myobsdata);
	return;
}

int main(int argc, char *argv[]) {
	const gsl_rng_type * T;
	long seed;
	gsl_rng_env_setup();
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);
	seed = time(NULL) * getpid();    // set seed
	gsl_rng_set(r, seed);
	clock_t begin, end;
	double time_spent;

	begin = clock();
	int its;
	printf("Starting main...\n");

	if (argc == 1) {
		its = 5;
	} else {
		its = atoi(argv[1]);
	}
	runModel(its);

	gsl_rng_free(r);

	end = clock();
	time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
	printf("Time consuming in total is %f\n", time_spent);

	printf("Running for %i", its);
	printf(" iterations, Done.");
	return 0;
}

max_stars_repo_name: scania/scania-pmcmc

max_issues_repo_name: scania/scania-pmcmc

max_forks_repo_name: scania/scania-pmcmc


text:
/* bst.c
 *
 * Copyright (C) 2018 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.
 */

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_bst.h>
#include <gsl/gsl_errno.h>

static void * bst_malloc(size_t size, void * params);
static void bst_free(void * block, void * params);

static const gsl_bst_allocator bst_default_allocator =
{
  bst_malloc,
  bst_free
};

/*
gsl_bst_alloc()
  Allocate binary search tree

Inputs: T         - tree type
        allocator - memory allocator
        compare   - comparison function
        params    - parameters to pass to allocator and compare
*/

gsl_bst_workspace *
gsl_bst_alloc(const gsl_bst_type * T, const gsl_bst_allocator * allocator,
              gsl_bst_cmp_function * compare, void * params)
{
  int status;
  gsl_bst_workspace *w;

  w = calloc(1, sizeof(gsl_bst_workspace));
  if (w == NULL)
    {
      GSL_ERROR_NULL("failed to allocate bst workspace", GSL_ENOMEM);
    }

  w->type = T;

  status = (w->type->init)(allocator != NULL ? allocator : &bst_default_allocator, compare, params, (void *) &w->table);
  if (status)
    {
      gsl_bst_free(w);
      GSL_ERROR_NULL("failed to initialize bst", GSL_EFAILED);
    }

  return w;
}

void
gsl_bst_free(gsl_bst_workspace * w)
{
  /* free tree nodes */
  gsl_bst_empty(w);

  free(w);
}

/* delete all nodes from tree */
int
gsl_bst_empty(gsl_bst_workspace * w)
{
  return (w->type->empty)((void *) &w->table);
}

/*
gsl_bst_insert()
  Inserts |item| into tree

AVL: if duplicate found, returns pointer to item without inserting
AVLmult: if duplicate found, increase multiplicity for that node

If no duplicate found, insert item and return pointer to item.
Returns NULL if a memory allocation error occurred.
*/

void *
gsl_bst_insert(void * item, gsl_bst_workspace * w)
{
  return (w->type->insert)(item, (void *) &w->table);
}

void *
gsl_bst_find(const void * item, const gsl_bst_workspace * w)
{
  return (w->type->find)(item, (const void *) &w->table);
}

void *
gsl_bst_remove(const void * item, gsl_bst_workspace * w)
{
  return (w->type->remove)(item, (void *) &w->table);
}

/* return number of nodes in tree */
size_t
gsl_bst_nodes(const gsl_bst_workspace * w)
{
  return (w->type->nodes)((const void *) &w->table);
}

/* return size (in bytes) of each node in tree */
size_t
gsl_bst_node_size(const gsl_bst_workspace * w)
{
  return w->type->node_size;
}

const char *
gsl_bst_name(const gsl_bst_workspace * w)
{
  return w->type->name;
}

/**********************************************
 * INTERNAL ROUTINES                          *
 **********************************************/

static void *
bst_malloc(size_t size, void * params)
{
  (void) params; /* avoid unused parameter warning */
  return malloc(size);
}

static void
bst_free(void * block, void * params)
{
  (void) params; /* avoid unused parameter warning */
  free(block);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense