Skip to content

Commit 73aea27

Browse files
ZocoLiniZocoLini
committed
printing the original times in the output file
1 parent d3e36d2 commit 73aea27

File tree

2 files changed

+7
-124
lines changed

2 files changed

+7
-124
lines changed

‎modules/cuqdyn-c/src/cuqdyn.c‎

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -74,12 +74,14 @@ CuqdynResult *cuqdyn_algo(const char *data_file, const char *sacess_conf_file, c
7474
exit(1);
7575
}
7676

77+
N_Vector scaled_times = N_VNew_Serial(NV_LENGTH_S(times), get_sundials_ctx());
78+
7779
for (int i = 0; i < NV_LENGTH_S(times); ++i)
7880
{
79-
NV_Ith_S(times, i) = NV_Ith_S(times, i) * config->time_scaling;
81+
NV_Ith_S(scaled_times, i) = NV_Ith_S(times, i) * config->time_scaling;
8082
}
8183

82-
const sunrealtype t0 = NV_Ith_S(times, 0);
84+
const sunrealtype t0 = NV_Ith_S(scaled_times, 0);
8385
N_Vector initial_condition = NULL;
8486

8587
if (config->y0.len == 0)
@@ -107,7 +109,7 @@ CuqdynResult *cuqdyn_algo(const char *data_file, const char *sacess_conf_file, c
107109
resid_loo = NewDenseMatrix(m, n);
108110
predicted_params_matrix = NewDenseMatrix(m, config->ode_expr.p_count);
109111

110-
N_Vector texp = copy_vector_remove_indices(times, create_array((long[]){}, 0));
112+
N_Vector texp = copy_vector_remove_indices(scaled_times, create_array((long[]){}, 0));
111113
SUNMatrix yexp = copy_matrix_remove_rows(observed_data, create_array((long[]){}, 0));
112114

113115
N_Vector tmp_initial_condition = copy_vector_remove_indices(initial_condition, create_array((long[]){}, 0));
@@ -164,7 +166,7 @@ CuqdynResult *cuqdyn_algo(const char *data_file, const char *sacess_conf_file, c
164166
{
165167
LongArray indices_to_remove = create_array((long[]){i + 1}, 1);
166168

167-
N_Vector texp = copy_vector_remove_indices(times, indices_to_remove);
169+
N_Vector texp = copy_vector_remove_indices(scaled_times, indices_to_remove);
168170
SUNMatrix yexp = copy_matrix_remove_columns(observed_data, indices_to_remove);
169171

170172
N_Vector tmp_initial_condition = copy_vector_remove_indices(initial_condition, create_array((long[]){}, 0));
@@ -173,7 +175,7 @@ CuqdynResult *cuqdyn_algo(const char *data_file, const char *sacess_conf_file, c
173175
execute_ess_solver(sacess_conf_file, output_file, texp, yexp, tmp_initial_condition, initial_params);
174176

175177
// Saving the ode solution data obtained with the predicted params
176-
TransposedStates ode_solution_states = solve_ode(predicted_params, initial_condition, t0, times);
178+
TransposedStates ode_solution_states = solve_ode(predicted_params, initial_condition, t0, scaled_times);
177179
ObservablesTransposedStates predicted_obs_states = transform_states(ode_solution_states);
178180

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

‎modules/cuqdyn-c/src/ode_solver.c‎

Lines changed: 0 additions & 119 deletions
Original file line numberDiff line numberDiff line change
@@ -12,121 +12,8 @@
1212
#include "cvodes/cvodes_ls.h"
1313
#include "sundials/sundials_types.h"
1414

15-
#ifndef JE
16-
#define JE(J,i,j) SM_ELEMENT_D((J),(i),(j))
17-
#endif
18-
1915
static int check_retval(void *, const char *, int);
2016

21-
int nfk_jac(sunrealtype t, N_Vector y, N_Vector fy,
22-
SUNMatrix J, void *user_data,
23-
N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
24-
{
25-
(void)t; (void)fy; (void)tmp1; (void)tmp2; (void)tmp3;
26-
27-
const sunrealtype *Y = NV_DATA_S(y);
28-
const sunrealtype y1=Y[0], y2=Y[1], y3=Y[2], y4=Y[3], y5=Y[4],
29-
y6=Y[5], y7=Y[6], y8=Y[7], y9=Y[8], y10=Y[9],
30-
y11=Y[10], y12=Y[11], y13=Y[12], y14=Y[13], y15=Y[14];
31-
32-
const sunrealtype *p = (const sunrealtype*) user_data;
33-
const sunrealtype p1 = p[0], p2 = p[1], p3 = p[2], p4 = p[3], p5 = p[4],
34-
p6 = p[5], p7 = p[6], p8 = p[7], p9 = p[8], p10= p[9],
35-
p11= p[10], p12= p[11], p13= p[12], p14= p[13], p15= p[14],
36-
p16= p[15], p17= p[16], p18= p[17], p19= p[18], p20= p[19],
37-
p21= p[20], p22= p[21], p23= p[22], p24= p[23], p25= p[24],
38-
p26= p[25], p27= p[26], p28= p[27], p29= p[28];
39-
40-
// Limpia la matriz (si tu versión de SUNDIALS no garantiza cero)
41-
for (int j = 0; j < SM_COLUMNS_D(J); ++j)
42-
{
43-
for (int i = 0; i < SM_ROWS_D(J); ++i)
44-
{
45-
SM_ELEMENT_D(J,i,j) = 0.0;
46-
}
47-
}
48-
49-
// f1
50-
JE(J, 0, 0) = -p21 -p17;
51-
52-
// f2
53-
JE(J, 1, 0) = p17;
54-
JE(J, 1, 1) = -p19 - p18*y8 - p21 - p2*y10 - p4*y13;
55-
JE(J, 1, 3) = p3;
56-
JE(J, 1, 4) = p5;
57-
JE(J, 1, 7) = -p18*y2;
58-
JE(J, 1, 9) = -p2*y2;
59-
JE(J, 1, 12) = -p4*y2;
60-
61-
// f3
62-
JE(J, 2, 1) = p19 + p18*y8;
63-
JE(J, 2, 2) = -p21;
64-
JE(J, 2, 7) = p18*y2;
65-
66-
// f4
67-
JE(J, 3, 1) = p2*y10;
68-
JE(J, 3, 3) = -p3;
69-
JE(J, 3, 9) = p2*y2;
70-
71-
// f5
72-
JE(J, 4, 1) = p4*y13;
73-
JE(J, 4, 4) = -p5;
74-
JE(J, 4, 12) = p4*y2;
75-
76-
// f6
77-
JE(J, 5, 4) = p5;
78-
JE(J, 5, 5) = -p1*y10 -p23;
79-
JE(J, 5, 9) = -p1*y6;
80-
JE(J, 5, 12) = p11;
81-
82-
// f7
83-
JE(J, 6, 5) = p23*p22;
84-
JE(J, 6, 6) = -p1*y11;
85-
JE(J, 6, 10) = -p1*y7;
86-
87-
// f8
88-
JE(J, 7, 7) = -p16;
89-
JE(J, 7, 8) = p15;
90-
91-
// f9
92-
JE(J, 8, 6) = p12;
93-
JE(J, 8, 8) = -p14;
94-
95-
// f10
96-
JE(J, 9, 1) = -p2*y10;
97-
JE(J, 9, 5) = -p1*y10;
98-
JE(J, 9, 9) = -p2*y2 - p1*y6 - p10 - p25;
99-
JE(J, 9, 10) = p26;
100-
JE(J, 9, 11) = p9;
101-
102-
// f11
103-
JE(J, 10, 6) = -p1*y11;
104-
JE(J, 10, 9) = p25*p22;
105-
JE(J, 10, 10) = -p1*y7 - p26*p22;
106-
107-
// f12
108-
JE(J, 11, 6) = p6;
109-
JE(J, 11, 11) = -p8;
110-
111-
// f13
112-
JE(J, 12, 1) = -p4*y13;
113-
JE(J, 12, 5) = p1*y10;
114-
JE(J, 12, 9) = p1*y6;
115-
JE(J, 12, 12) = -p11 -p4*y2;
116-
JE(J, 12, 13) = p24;
117-
118-
// f14
119-
JE(J, 13, 6) = p1*y11;
120-
JE(J, 13, 10) = p1*y7;
121-
JE(J, 13, 13) = -p24*p22;
122-
123-
// f15
124-
JE(J, 14, 6) = p27;
125-
JE(J, 14, 14) = -p29;
126-
127-
return 0; // éxito
128-
}
129-
13017

13118
TransposedStates solve_ode(N_Vector parameters, N_Vector initial_values, sunrealtype t0, N_Vector times)
13219
{
@@ -169,12 +56,6 @@ TransposedStates solve_ode(N_Vector parameters, N_Vector initial_values, sunreal
16956
{
17057
return NULL;
17158
}
172-
173-
retval = CVodeSetJacFn(cvode_mem, nfk_jac);
174-
if (check_retval(&retval, "CVodeSetJacFn", 1))
175-
{
176-
return NULL;
177-
}
17859

17960
char *min_step_s = getenv("SACESS_CVODES_MIN_STEP");
18061
if (min_step_s != NULL)

0 commit comments

Comments
 (0)