Commit dd130a98 authored by Thierry Martinez's avatar Thierry Martinez

Plotting

parent bbeade81
...@@ -41,6 +41,7 @@ write_gsl(Options) :- ...@@ -41,6 +41,7 @@ write_gsl(Options) :-
write_jacobian(Options), write_jacobian(Options),
write_initial_values(Options), write_initial_values(Options),
write_initial_parameter_values(Options), write_initial_parameter_values(Options),
write_print_headers(Options),
write_print(Options). write_print(Options).
write_headers(Options) :- write_headers(Options) :-
...@@ -228,6 +229,35 @@ initial_parameter_values(double p[]) ...@@ -228,6 +229,35 @@ initial_parameter_values(double p[])
'). ').
write_print_headers(Options) :-
(
memberchk(headers: Headers, Options)
->
write('
void
print_headers(FILE *csv)
{
fprintf(csv, "#Time'),
\+ (
member(Header, Headers),
\+ (
format(',\\"~a\\"', [Header])
)
),
write('\\n"'),
write(');
}
')
;
write('
void
print_headers(FILE *csv)
{
}
')
).
write_print(Options) :- write_print(Options) :-
memberchk(equations: Equations, Options), memberchk(equations: Equations, Options),
length(Equations, VariableCount), length(Equations, VariableCount),
......
...@@ -22,6 +22,7 @@ main(void) ...@@ -22,6 +22,7 @@ main(void)
double x[VARIABLE_COUNT]; double x[VARIABLE_COUNT];
initial_values(x); initial_values(x);
initial_parameter_values(parameters); initial_parameter_values(parameters);
print_headers(csv);
while (t < TIME_FINAL) { while (t < TIME_FINAL) {
int status = gsl_odeiv2_evolve_apply( int status = gsl_odeiv2_evolve_apply(
e, c, s, &sys, &t, TIME_FINAL, &h, x); e, c, s, &sys, &t, TIME_FINAL, &h, x);
......
...@@ -28,10 +28,12 @@ make_ode_system :- ...@@ -28,10 +28,12 @@ make_ode_system :-
solve :- solve :-
gather_headers(Headers),
gather_equations(Equations), gather_equations(Equations),
gather_initial_values(InitialValues), gather_initial_values(InitialValues),
gather_initial_parameter_values(InitialParameterValues), gather_initial_parameter_values(InitialParameterValues),
Options = [ Options = [
headers: Headers,
equations: Equations, equations: Equations,
initial_values: InitialValues, initial_values: InitialValues,
initial_parameter_values: InitialParameterValues, initial_parameter_values: InitialParameterValues,
...@@ -46,6 +48,20 @@ solve :- ...@@ -46,6 +48,20 @@ solve :-
add_trace('numerical_simulation', Table). add_trace('numerical_simulation', Table).
gather_headers(Headers) :-
peek_count(variable_counter, VariableCount),
VariableMax is VariableCount - 1,
findall(
Header,
(
between(0, VariableMax, VariableIndex),
variable(Molecule, VariableIndex),
format(atom(Header), '[~a]', [Molecule])
),
Headers
).
gather_equations(Equations) :- gather_equations(Equations) :-
peek_count(variable_counter, VariableCount), peek_count(variable_counter, VariableCount),
VariableMax is VariableCount - 1, VariableMax is VariableCount - 1,
...@@ -96,6 +112,7 @@ enumerate_variables :- ...@@ -96,6 +112,7 @@ enumerate_variables :-
\+ ( \+ (
ode(Molecule, _), ode(Molecule, _),
\+ ( \+ (
atom(Molecule),
count(variable_counter, VariableIndex), count(variable_counter, VariableIndex),
assertz(variable(Molecule, VariableIndex)) assertz(variable(Molecule, VariableIndex))
) )
......
...@@ -20,8 +20,8 @@ list_ODE :- ...@@ -20,8 +20,8 @@ list_ODE :-
compute_ode :- compute_ode :-
retractall(ode(_, _)), retractall(ode(_, _)),
\+ ( \+ (
item([model: current_model, kind: rule, item: Item]), item([model: current_model, kind: reaction, item: Item]),
rule(Item, Kinetics, Reactants, Products), reaction(Item, Kinetics, Reactants, Products),
\+ ( \+ (
kinetics(Kinetics, Reactants, KineticsExpression), kinetics(Kinetics, Reactants, KineticsExpression),
add_molecules(Reactants, -KineticsExpression), add_molecules(Reactants, -KineticsExpression),
......
:- module(
plot,
[
plot/0,
export_plot/1
]
).
plot :-
biocham_command,
doc('plots the current trace.'),
export_plot(plot),
process_create(path(gnuplot), ['-persist', 'plot.plot'], []).
export_plot(FileTemplate) :-
biocham_command,
doc('
saves the numerical trace of the last simulation into two files:
\\argument{FileTemplate}\\texttt{.csv} and \\texttt{.plot}.
'),
format(atom(PlotFile), '~a.plot', [FileTemplate]),
format(atom(CsvFile), '~a.csv', [FileTemplate]),
export_trace(CsvFile),
setup_call_cleanup(
open(PlotFile, write, Stream),
export_plot_stream(Stream, CsvFile),
close(Stream)
).
export_plot_stream(Stream, CsvFile) :-
format(Stream, '\c
set style data lines
set format y "%g"
set xrange [*:*]
set yrange [*:*]
set ytics nomirror autofreq
set xtics nomirror
set mxtics default
set key outside Left reverse
set datafile separator ","
plot "~a"
', [CsvFile]).
...@@ -55,6 +55,8 @@ ...@@ -55,6 +55,8 @@
- numerical_simulation.pl - numerical_simulation.pl
*** Traces *** Traces
- traces.pl - traces.pl
*** Plotting the result of simulations
- plot.pl
* Index * Index
+ index + index
* Bibliography * Bibliography
......
...@@ -139,7 +139,7 @@ command(Command) :- ...@@ -139,7 +139,7 @@ command(Command) :-
Command = (_ for _) Command = (_ for _)
) )
-> ->
add_reaction(Command) command(add_reaction(Command))
; ;
throw(error(unknown_command(Functor/Arity))) throw(error(unknown_command(Functor/Arity)))
). ).
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment