Commit e43c4067 authored by REMONDEAU Paul's avatar REMONDEAU Paul
Browse files

small changes to hybrid simulation, need to fix stochastic_conversion rate

parent d0dc0682
......@@ -30,9 +30,10 @@ hybrid_static_simulation(ODEFilename, StochFilename) :-
The first input file contains the informations concerning the continous
reactions. The second input file contains the informations concerning the
stochastic reactions.'),
hybrid_static_simulation(ODEFilename, StochFilename, 'out.bc', 1, 200).
hybrid_static_simulation(ODEFilename, StochFilename, 'out_static.bc', 1, 200).
hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time) :-
biocham_command,
%============================%
%== Preprocessing ODE file ==%
%============================%
......@@ -260,13 +261,15 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
write(StreamAppend, ').'),
nl(StreamAppend),
write(StreamAppend, 'plot.'),
nl(StreamAppend),
write(StreamAppend,'option(show: {}).'),
close(StreamAppend),
load_biocham(OutFilename),
format(" ~n If you want to choose which species are plot, please change the corresponging ~n line in the output file (by default out.bc in the biocham directory). ~n ~n",[]).
format(" ~n If you want to choose which species are plot, please change the corresponding ~n line in the output file (by default out_static.bc in the biocham directory). ~n ~n",[]).
:- doc('\\begin{example}T7 model using hybrid static simulation:').
:- biocham_silent(clear_model).
:- biocham('hybrid_static_simulation(library:examples/hybid/T7_input_ode,library:examples/hybid/T7_input_stoch').
:- biocham(hybrid_static_simulation('library:examples/hybrid/T7_input_ode.bc','library:examples/hybrid/T7_input_stoch.bc')).
:- doc('\\end{example}').
%===========================================================%
......@@ -403,7 +406,9 @@ convert_MA(Reaction,Reactions,Hybrid_originals,NewKinetics) :-
nth1(ReactionNumber,Reactions,Reaction),
atom_concat('k',ReactionNumber,Kdiff1),
atom_concat(Kdiff1,'_diff',Kdiff),
NewKinetics = Kdiff * K * NewValue
NewPartialKinetics = K * NewValue,
assertz(reaction(Reaction,NewPartialKinetics)),
NewKinetics = Kdiff * NewPartialKinetics
;
NewKinetics = K * NewValue
).
......@@ -1100,6 +1105,8 @@ hybrid_dynamic_simulation(InputFile,PropTresh,PartTresh) :-
Then, this new fresh file is loaded and a numerical simulation with the Rosenbrock method takes place */
hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh) :-
biocham_command,
%get_option(stochastic_conversion, Rate)
load_biocham(InputFile),
open(OutFilename,write,Stream),
all_items([kind: reaction],Reactions),
......@@ -1204,7 +1211,10 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
prepare_change_list(Reactions,SortedSpecies,List),
%format("List is ~w ~n",[List]),
assertz(maximum_change(0)),
assertz(minimum_reactant_change(inf)),
assertz(maximum_reactant_change(0)),
assertz(minimum_product_change(inf)),
assertz(maximum_product_change(0)),
findall([Species_change, Alpha],
(
member(Element,List),
......@@ -1226,16 +1236,7 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
%== Build list (Species_change) to record species change ==%
particle_sum_list(RSpecies,RSum),
particle_sum_list(PSpecies,PSum),
NewCoeff is abs(PSum-RSum),
retract(maximum_change(OldMax)),
%format("OldMax is ~w NewCoeff is ~w ~n",[OldMax,NewCoeff]),
(
NewCoeff > OldMax
->
assertz(maximum_change(NewCoeff))
;
assertz(maximum_change(OldMax))
),
compute_partial_maximum_change(RSum,PSum),
append_reactants(RSpecies, 1, [], Temp_list, SortedSpecies),
append_products(PSpecies, 1, Temp_list, Changes, SortedSpecies),
build_change_list(Changes, 1, Species_number, [], Species_change)
......@@ -1256,7 +1257,7 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
),
HybridDynamicReactants),
retract(maximum_change(MaxParticleChange)),
compute_maximum_change(MaxParticleChange),
ParticuleTreshold is MaxParticleChange*PartTresh,
%format("Hybrid_dynamic_reactants is ~w ~n",[HybridDynamicReactants]),
%format("Reactions_data is ~w ~n",[Reactions_data]),
......@@ -1266,19 +1267,7 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
nb_delete(dynamic_parameters),
nb_delete(reactions),
nl(StreamAppend),
%print_common_alpha_macros(Reactions_data, 1, StreamAppend),
/*
% Plotting Part; depends on what the user want to show on the plot.
write(StreamAppend, 'show_macros.'), nl(StreamAppend),
write(StreamAppend, 'show_parameters.'), nl(StreamAppend),
write(StreamAppend, 'show_macros({'),
concat_species2(Hybrids, '', Hybrid_Results),
write(StreamAppend, Hybrid_Results),
write(StreamAppend, '}).'), nl(StreamAppend),
*/
%write(StreamAppend, 'show_macros.'), nl(StreamAppend),
write(StreamAppend,'% Plot part'),nl(StreamAppend),nl(StreamAppend),
print_show_option(SortedSpecies,StreamAppend),
write(StreamAppend, 'numerical_simulation(method:rsbk, time: '),
......@@ -1286,9 +1275,54 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
write(StreamAppend, ').'),
nl(StreamAppend),
write(StreamAppend, 'plot.'),
nl(StreamAppend),
write(StreamAppend,'option(show: {}).'),
close(StreamAppend),
load_biocham(OutFilename),
format(" ~n If you want to choose which species are plot, please change the corresponging ~n line in the output file (by default out_dynamic.bc in the biocham directory). ~n ~n",[]).
format(" ~n If you want to choose which species are plot, please change the corresponding ~n line in the output file (by default out_dynamic.bc in the biocham directory). ~n ~n",[]).
compute_partial_maximum_change(RSum,PSum) :-
retract(minimum_reactant_change(OldRMin)),
retract(maximum_reactant_change(OldRMax)),
retract(minimum_product_change(OldPMin)),
retract(maximum_product_change(OldPMax)),
(
RSum < OldRMin
->
assertz(minimum_reactant_change(RSum))
;
assertz(minimum_reactant_change(OldRMin))
),
(
RSum > OldRMax
->
assertz(maximum_reactant_change(RSum))
;
assertz(maximum_reactant_change(OldRMax))
),
(
PSum < OldPMin
->
assertz(minimum_product_change(PSum))
;
assertz(minimum_product_change(OldPMin))
),
(
PSum > OldPMax
->
assertz(maximum_product_change(PSum))
;
assertz(maximum_product_change(OldPMax))
).
compute_maximum_change(MaxParticleChange) :-
retract(minimum_reactant_change(RMin)),
retract(maximum_reactant_change(RMax)),
retract(minimum_product_change(PMin)),
retract(maximum_product_change(PMax)),
Max1 is abs(PMax-RMin),
Max2 is abs(RMax-PMin),
MaxParticleChange is max(Max1,Max2).
particle_sum_list([],0).
......@@ -1306,10 +1340,8 @@ print_dynamic_event_condition(HybridDynamicReactants,ReactionNumber,PropTresh,Ma
write(Stream,'_diff = (if '),
nb_getval(reactions,Reactions),
nth1(ReactionNumber,Reactions,Reaction),
reaction(Reaction, [kinetics: Kinetics,reactants: _,inhibitors: _,products: _]),
Kinetics =.. [_,K],
write(Stream,K),
retract(reaction(Reaction,NewKinetics)),
write(Stream,NewKinetics),
write(Stream, ' >= 1/(step_size * '),
write(Stream,PropTresh),
write(Stream,')'),
......@@ -1359,7 +1391,7 @@ print_dynamic_event(Reactions_data, HybridDynamicReactants,MaxParticleChange,Pro
format("Checkpoint ~n",[]),
print_dynamic_event_condition(HybridDynamicReactants,1,PropTresh,MaxParticleChange,Stream),
print_static_alpha(Number_of_alpha,Reactions_data, 1, Stream),
write(Stream,' delta_t = (if alpha_sum <= step_size then inf else (-1 / alpha_sum) * log(random_float)),'),nl(Stream),
write(Stream,' delta_t = (if alpha_sum <= 0 then step_size else (-1 / alpha_sum) * log(random_float)),'),nl(Stream),
write(Stream,' tau = Time+delta_t,'),nl(Stream),
write(Stream,' ran = random_float,'),nl(Stream),
%print_event_constraints(Constraints_list, 1, Stream),
......
MA(c1) for a => a + b.
MA(c2) for b => _.
MA(c3) for c => a + c.
MA(c4) for a => _.
MA(c5) for d + e => c.
MA(c6) for c => d + e.
MA(c7) for c + e => f.
MA(c8) for f => c + e.
MA(c9) for b + b => e.
MA(c10) for e => b + b.
parameter(c1 = 0.043).
parameter(c2 = 7e-4).
parameter(c3 = 71.5).
parameter(c4 = 3.9e-6).
parameter(c5 = 0.02).
parameter(c6 = 0.48).
parameter(c7 = 4e-4).
parameter(c8 = 9e-12).
parameter(c9 = 0.08).
parameter(c10 = 0.5).
present(a,0).
present(b,10).
present(c,0).
present(d,2).
present(e,10).
present(f,0).
MA(c9) for b + b => e.
MA(c10) for e => b + b.
parameter(c9 = 0.08).
parameter(c10 = 0.5).
present(b,10).
present(e,10).
MA(c1) for a => a + b.
MA(c2) for b => _.
MA(c3) for c => a + c.
MA(c4) for a => _.
MA(c5) for d + e => c.
MA(c6) for c => d + e.
MA(c7) for c + e => f.
MA(c8) for f => c + e.
parameter(c1 = 0.043).
parameter(c2 = 7e-4).
parameter(c3 = 71.5).
parameter(c4 = 3.9e-6).
parameter(c5 = 0.02).
parameter(c6 = 0.48).
parameter(c7 = 4e-4).
parameter(c8 = 9e-12).
present(a,0).
present(b,10).
present(c,0).
present(d,2).
present(e,10).
present(f,0).
......@@ -198,7 +198,8 @@ update_functions_state(CurrentState,Time) :-
FunctionValue,
(member(function(_,Expression),Functions),
eval_coeff(Expression,CurrentState,Parameters,Time,FunctionValue1),
simplify(FunctionValue1,FunctionValue)
simplify(FunctionValue1,FunctionValue2),
FunctionValue is FunctionValue2
),
FunctionsValues
),
......@@ -431,6 +432,7 @@ rosenbrock_step(Equations,Jacobian,CurrentState,Parameters,NewCurrentState) :-
eval_matrix(Jacobian,TrueCurrentState,Parameters,JacobianCurrentState), % Create a choicepoint which I can't find a way to prevent
nb_getval(h,StepSize),
nb_setval(hdid,StepSize),
format("Ah ok ~n",[]),
length(Equations,Length),
% evaluate Jacobian at current time point
......
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