/idmclib/src/traj.c
C | 159 lines | 113 code | 22 blank | 24 comment | 24 complexity | fb4e248e876254e0716edb81ba15aac4 MD5 | raw file
Possible License(s): GPL-2.0
- /*
- iDMC C library
- Adapted from iDMC, Copyright (C) 2004-2006 Marji Lines and Alfredo Medio
- Copyright (C) 2006,2007 Marji Lines and Alfredo Medio.
- 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 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.
- Last modified: $Date: 2007-01-06 14:59:13 +0100 (Sat, 06 Jan 2007) $
- */
- #include <stdlib.h>
- #include <gsl/gsl_errno.h>
- #include <gsl/gsl_odeiv.h>
- #include "model.h"
- #include "traj.h"
- #define SAFE_FREE(p) if(p!=NULL) free(p)
- int idmc_traj_trajectory_alloc(idmc_model *model, double *parValues, double *varValues, idmc_traj_trajectory **out_ans) {
- idmc_traj_trajectory *ans;
- *out_ans = (idmc_traj_trajectory *) calloc(1, sizeof(idmc_traj_trajectory));
- if((*out_ans)==NULL)
- return IDMC_EMEM;
- ans = *out_ans;
- ans->model = idmc_model_clone(model);
- if(!(ans->model)) {
- idmc_traj_trajectory_free(ans);
- return IDMC_EMEM;
- }
- ans->par = (double*) malloc(ans->model->par_len * sizeof(double));
- if(!(ans->par)) {
- idmc_traj_trajectory_free(ans);
- return IDMC_EMEM;
- }
- memcpy(ans->par, parValues, model->par_len * sizeof(double));
- ans->var = (double*) malloc(ans->model->var_len * sizeof(double));
- if(!ans->var) {
- idmc_traj_trajectory_free(ans);
- return IDMC_EMEM;
- }
- memcpy(ans->var, varValues, model->var_len * sizeof(double));
- ans->step = 0;
- *out_ans = ans;
- return IDMC_OK;
- }
- void idmc_traj_trajectory_free(idmc_traj_trajectory *traj){
- if(traj->model != NULL)
- idmc_model_free(traj->model);
- SAFE_FREE(traj->par);
- SAFE_FREE(traj->var);
- SAFE_FREE(traj);
- }
- int idmc_traj_trajectory_step(idmc_traj_trajectory *t) {
- int result = idmc_model_f(t->model, t->par, t->var /*from point*/, t->var /*to point*/);
- (t->step)++;
- return result;
- }
- /*Continuous systems*/
- static int odeiv_function(double t, const double yy[], double dydt[], void * params);
- int idmc_traj_ctrajectory_alloc(idmc_model *model, double *parValues, double *varValues,
- double step_size, gsl_odeiv_step_type *step_function_code, idmc_traj_ctrajectory **ans) {
- idmc_traj_ctrajectory *trajectory;
-
- if(model->type[0] != 'C')
- return IDMC_EMATH;
- trajectory = (idmc_traj_ctrajectory *) malloc(sizeof (idmc_traj_ctrajectory));
- if (trajectory == NULL) {
- return IDMC_EMEM;
- }
- /* model */
- trajectory->model = idmc_model_clone(model);
- if(!trajectory->model) {
- idmc_traj_ctrajectory_free(trajectory);
- return IDMC_EMEM;
- }
- /* parameters */
- trajectory->par = malloc(model->par_len * sizeof(double));
- if (trajectory->par == NULL) {
- idmc_traj_ctrajectory_free(trajectory);
- return IDMC_EMEM;
- }
- memcpy(trajectory->par, parValues, model->par_len * sizeof(double));
- /* initial values */
- trajectory->var = (double*) malloc(model->var_len * sizeof(double));
- if (trajectory->var == NULL) {
- idmc_traj_ctrajectory_free(trajectory);
- return IDMC_EMEM;
- }
- memcpy(trajectory->var, varValues, model->var_len * sizeof(double));
- /* initial values */
- trajectory->error = (double*) malloc(model->var_len * sizeof(double));
- if (trajectory->error == NULL) {
- idmc_traj_ctrajectory_free(trajectory);
- return IDMC_EMEM;
- }
-
- trajectory->step_function_code = step_function_code;
- trajectory->step_function = gsl_odeiv_step_alloc(step_function_code, model->var_len);
- if (trajectory->step_function == NULL) {
- idmc_traj_ctrajectory_free(trajectory);
- return IDMC_EMEM;
- }
-
- trajectory->step_size = step_size;
-
- /* gsl ode system */
- trajectory->system.function = odeiv_function;
- trajectory->system.jacobian = NULL; // FIXME
- trajectory->system.dimension = model->var_len;
- trajectory->system.params = trajectory;
-
- *ans = trajectory;
- return IDMC_OK;
- }
- void idmc_traj_ctrajectory_free(idmc_traj_ctrajectory *trajectory) {
- if(trajectory->model!=NULL)
- idmc_model_free(trajectory->model);
- if (trajectory->step_function!=NULL)
- gsl_odeiv_step_free(trajectory->step_function);
- SAFE_FREE(trajectory->error);
- SAFE_FREE(trajectory->var);
- SAFE_FREE(trajectory->par);
- SAFE_FREE(trajectory);
- }
- int idmc_traj_ctrajectory_step(idmc_traj_ctrajectory *trajectory) {
- return gsl_odeiv_step_apply(trajectory->step_function, 0, trajectory->step_size,
- trajectory->var, trajectory->error, NULL, NULL, &trajectory->system);
- }
- static int odeiv_function(double t, const double yy[], double dydt[], void * params) {
- idmc_traj_ctrajectory *trajectory = params;
- int err;
- err = idmc_model_f(trajectory->model, trajectory->par, yy, dydt);
- if (err != IDMC_OK)
- return GSL_EFAILED;
- return GSL_SUCCESS;
- }