Commit b313399f authored by Mathieu Hemery's avatar Mathieu Hemery
Browse files

Reformat the hybrid_static_simulation

parent a1d1a828
......@@ -20,9 +20,6 @@
:- use_module(biocham).
:- use_module(numerical_simulation).
/* I don't know what to do with the stochastic conversion rate
It seems like species names like ABC-DFG-{HIJ} make the simulation not work, typical example MAPK*/
:- doc('This section describes a hybrid numerical simulation
method. Some reactions which occurs more often and which have
a lot of particles reactants will be seen as continuous,
......@@ -46,7 +43,7 @@ hybrid_static_simulation(ODEFilename, StochFilename) :-
reactions. The second input file contains the informations concerning the
stochastic reactions. Please list out all the species initial value
(by present(_,_)) in input files. '),
hybrid_static_simulation(ODEFilename, StochFilename, 'out_static.bc', '1', 20).
hybrid_static_simulation(ODEFilename, StochFilename, 'out_static.bc', 1, 20).
%! hybrid_static_simulation(+ODEFilename, +StochFilename, +OutFilename, +Volume, +Time)
......@@ -84,6 +81,7 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
sort(Reactants1,SortedReactants), % Used for later hybrid-reactants finding
sort(ODESpecies,SortedODESpecies),
clear_model,
%===================================%
%== Preprocessing Stochastic file ==%
%===================================%
......@@ -100,124 +98,104 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
% process the parameters and presents
load_biocham(StochFilename),
open(OutFilename, write, Stream),
format(Stream, "% Hybrid version of ~w and ~w ~n", [ODEFilename, StochFilename]),
format(Stream, "parameter(tau = 0).~n", []),
format(Stream, "seed(0).~nparameter(ran = 0).~n", []),
all_items([kind: parameter],Parameters),
findall(Parameter,
(
member(Parameter,Parameters),
write(Stream, Parameter),
write(Stream,'.'),
nl(Stream)
),
_),
format(Stream, "% Hybrid version of ~w and ~w ~n", ODEFilename, StochFilename),
write(Stream, 'parameter(tau = 0).'), nl(Stream),
Seed is 0,
format(Stream, "seed(~d).~n", [Seed]),
write(Stream, 'parameter(ran = '),
write(Stream, Seed),
write(Stream, ').'), nl(Stream),
%write(Stream, 'macro(rand, 1664525 * ran + 1013904223 / 4294967296 - '),
%write(Stream, 'floor(1664525 * ran + 1013904223 / 4294967296)).'),nl(Stream),
maplist(format(Stream, "~w.~n"), Parameters),
all_items([kind:initial_state],InitState),
findall(Hybrid_original,(
(member(present(Hybrid_original),InitState);member(present(Hybrid_original,_),InitState)),
(
member(present(Hybrid_original),InitState)
;
member(present(Hybrid_original,_),InitState)
),
member(Hybrid_original, SortedODESpecies)
), Hybrid_originals),
),
Hybrid_original_list),
get_option(stochastic_conversion, Rate),
findall(Hybrid_species,(
((member(present(Present_pre),InitState),Init_value_pre = 1);member(present(Present_pre,Init_value_pre),InitState)),
(
(member(present(Present_pre),InitState),Init_value_pre = 1)
;
member(present(Present_pre,Init_value_pre),InitState)
),
(
member(Present_pre, SortedODESpecies)
->
->
(
species_to_stoch(Present_pre, Present),
Init_value = 0,
species_to_total(Present_pre, Hybrid_species),
write(Stream, 'function('),
write(Stream, Hybrid_species),
write(Stream,' = floor('),
write(Stream,Rate),
write(Stream,'*'),
write(Stream, Volume),
write(Stream, '*'),
write(Stream, Present_pre),
write(Stream, ' + '),
write(Stream, Present),
write(Stream, ')).'),
nl(Stream)
format(Stream,"function(~w = floor((~g*~g*~w) + ~w)).~n",
[Hybrid_species, Rate, Volume, Present_pre, Present])
)
;
;
(
species_to_stoch(Present_pre, Present),
Init_value = Init_value_pre,
species_to_total(Present_pre, Hybrid_species),
write(Stream, 'function('),
write(Stream, Hybrid_species),
write(Stream,' = '),
write(Stream, Present),
write(Stream, ').'),
nl(Stream)
format(Stream,"function(~w = ~w).~n",
[Hybrid_species, Present])
)
),
write(Stream, 'parameter('), write(Stream, Present),
write(Stream, ' = '), write(Stream, Init_value),
write(Stream, ').'), nl(Stream)
),_Hybrids),
format(Stream, "parameter(~w = ~g).~n",
[Present, Init_value])
),_Hybrid_list),
clear_model,
load_biocham(ODEFilename),
nl(Stream),write(Stream,'% ODE part'),nl(Stream),nl(Stream),
write_hybrid_ode(Stream,Hybrid_originals),
format(Stream, "~n% ODE part~n~n", []),
write_hybrid_ode(Stream,Hybrid_original_list),
clear_model,
nl(Stream),write(Stream,'% Hybrid event part'),nl(Stream),nl(Stream),
format(Stream, "~n% Hybrid event part~n~n", []),
load_biocham(StochFilename),
all_items([kind: reaction],Reactions2),
nb_setval(reactions,Reactions2),
findall(R2,
(
member(Reaction,Reactions2),
reaction(Reaction, [kinetics: _,reactants: Reactants,inhibitors: _, products: _]),
member(Reactant,Reactants),
Reactant = _*R2_pre,
species_to_stoch(R2_pre, R2)
),
Reactants2),
% Add constraints for hybrid species! (28/9/12)
% For hybrid species, the determinisic part can make the particle number
% lower than one discrete step needed for stochastic simulation.
findall(Reactant,
(
member(Reaction,Reactions2),
reaction(Reaction, [kinetics: _,reactants: Reactants,inhibitors: _, products: _]),
(
member(Reaction,Reactions2),
reaction(Reaction, [kinetics: _,reactants: Reactants,inhibitors: _, products: _]),
(
Reactants = []
->
Reactant = none
;
member(Reactant,Reactants)
)
),
Reactants_by_rules),
Reactants = []
->
Reactant = none
;
member(Reactant,Reactants)
)
),
Reactants_by_rules),
format("Reactants_by_rules is ~w ~n",[Reactants_by_rules]),
build_constraints_by_rules(Reactants_by_rules, SortedReactants, 1, [], Hybrid_reactants_list),
build_constraints_list(Hybrid_reactants_list, 1, [], Constraints_list),
findall(P2,
findall(Reactant_stoch,
(
member(Reaction,Reactions2),
reaction(Reaction, [kinetics: _,reactants: Reactants,inhibitors: _, products: _]),
member(Reactant,Reactants),
Reactant = _*Reactant_pre,
species_to_stoch(Reactant_pre, Reactant_stoch)
),
Reactant_stoch_list),
findall(Product_stoch,
(
member(Reaction,Reactions2),
reaction(Reaction, [kinetics: _,reactants: _,inhibitors: _, products: Products]),
member(Product,Products),
Product = _*P2_pre,
species_to_stoch(P2_pre, P2)
Product = _*Product_pre,
species_to_stoch(Product_pre, Product_stoch)
),
Products2),
append(Reactants2, Products2, StochSpecies),
Product_stoch_list),
append(Reactant_stoch_list, Product_stoch_list, StochSpecies),
sort(StochSpecies,SortedStochSpecies), % To eliminate duplication.
% build StochReactants(with '_total') from SortedStochSpecies.
......@@ -229,7 +207,7 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
%========================================%
prepare_change_list(Reactions2,ODESpecies,List),
findall([Species_change, Alpha],
(
(
member(Element,List),
% To get [RSpecies, PSpecies] from [Reactants, Products]
Element = [Reactants,Products,Alpha],
......@@ -247,35 +225,22 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
species_to_stoch(Pname, Species2)
),
PSpecies),
%== Build list (Species_change) to record species change ==%
%== Build list (Species_change) to record species change ==%
append_reactants(RSpecies, 1, [], Temp_list, SortedStochSpecies),
append_products(PSpecies, 1, Temp_list, Changes, SortedStochSpecies),
build_change_list(Changes, 1, Species_number, [], Species_change)
),
Reactions_data
),
),
Reactions_data),
%== Build Events and Calaculate Alphas==%
close(Stream),
open(OutFilename, append, StreamAppend),
format("Hybrid_reactants_list is ~w ~n",[Hybrid_reactants_list]),
print_event(Reactions_data, SortedStochSpecies, StochReactants,Hybrid_reactants_list, Constraints_list, StreamAppend),
nl(StreamAppend),
%print_common_alpha_macros(Reactions_data, 1, StreamAppend),
print_event(Reactions_data, SortedStochSpecies, StochReactants,Hybrid_reactants_list, Constraints_list, Stream),
write(StreamAppend,'% Plot part'),nl(StreamAppend),nl(StreamAppend),
format("Event printed ~w ~n",[Hybrid_originals]),
%print_show_option(Hybrid_originals,StreamAppend),
write(StreamAppend,'option(show: {}). % You can choose here which sepcies you want to plot : option(show: {a_total,b_total,c_total,d_total}) to plot species a,b,c and d.'),
nl(StreamAppend),
write(StreamAppend, 'numerical_simulation(method:rsbk, time: '),
write(StreamAppend, Time),
write(StreamAppend, ').'),
nl(StreamAppend),
write(StreamAppend, 'plot.'),
nl(StreamAppend),
write(StreamAppend,'option(show: {}).'),
close(StreamAppend),
format(Stream, "~n% Plot part~n~n", []),
format(Stream,"option(show: {}). % You can choose here which species you want to plot : option(show: {a_total,b_total,c_total,d_total}) to plot species a,b,c and d.~n",[]),
format(Stream, "numerical_simulation(method:rsbk, time: ~g).~n", Time),
format(Stream, "plot.~noption(show: {}).~n", []),
close(Stream),
nb_delete(volume),
nb_delete(reactions),
retractall(reaction_kinetics(_,_)),
......@@ -393,14 +358,14 @@ build_reactants(StochSpecies, ODESpecies, Species_number, Counter, Init, StochRe
StochReactants = Next
).
write_hybrid_ode(Stream,Hybrid_originals) :-
list_hybrid_ode(Stream,Hybrid_originals),
write_hybrid_ode(Stream,Hybrid_original_list) :-
list_hybrid_ode(Stream,Hybrid_original_list),
with_output_to(Stream, models:list_model_initial_state),
with_output_to(Stream, models:list_model_parameters),
with_output_to(Stream, models:list_model_events),
with_output_to(Stream, models:list_model_macros(function)).
list_hybrid_ode(Stream,Hybrid_originals) :- % Only MA are converted, need to do others
list_hybrid_ode(Stream,Hybrid_original_list) :- % Only MA are converted, need to do others
all_items([kind: reaction], Reactions),
forall(
member(Reaction,Reactions),
......@@ -408,13 +373,13 @@ list_hybrid_ode(Stream,Hybrid_originals) :- % Only MA are converted, need to do
Reaction =.. [Op,_,_],
Op = (=>)
->
convert_MA(Reaction,Reactions,Hybrid_originals,NewKinetics),
convert_MA(Reaction,Reactions,Hybrid_original_list,NewKinetics),
NewReaction =.. [for,NewKinetics,Reaction],
write(Stream,NewReaction),write(Stream,'.'),nl(Stream)
;
Reaction =.. [Op,_,_], % Otherwise says Singleton variable in branch: Op
Reaction =.. [Op,_,SecondPart],
convert_MA(Reaction,Reactions,Hybrid_originals,NewKinetics),
convert_MA(Reaction,Reactions,Hybrid_original_list,NewKinetics),
NewReaction =.. [for,NewKinetics,SecondPart],
write(Stream,NewReaction),write(Stream,'.'),nl(Stream)
)
......@@ -424,11 +389,11 @@ list_hybrid_ode(Stream,Hybrid_originals) :- % Only MA are converted, need to do
/* Can try to get the kinetics function in BIOCHAM 4 directly and translate, but i didn't know how
Also can direclty adapt the kinetics with the hybrid species list, but need to translate the kinetics first if needed */
convert_MA(Reaction,Reactions,Hybrid_originals,NewKinetics) :-
convert_MA(Reaction,Reactions,Hybrid_original_list,NewKinetics) :-
reaction(Reaction, [kinetics: Kinetics,reactants: Reactants,inhibitors: _, products: _]),
Kinetics =.. ['MA',K],
kinetics:eval_kinetics(Reactants, _, product('S'*'M' in [reactants], 'M'^'S'), Value),
convert_into_hybrid(Value,Hybrid_originals,NewValue),!,
convert_into_hybrid(Value,Hybrid_original_list,NewValue),!,
(
nb_current(dynamic_parameters,_)
->
......@@ -449,26 +414,26 @@ convert_into_hybrid(Value ^ 1, ODESpecies,NewValue) :-
convert_into_hybrid(Value,ODESpecies,NewValue).
convert_into_hybrid(Value1 * Value2, Hybrid_originals,NewValue1 * NewValue2) :-
convert_into_hybrid(Value1,Hybrid_originals,NewValue1),
convert_into_hybrid(Value2,Hybrid_originals,NewValue2).
convert_into_hybrid(Value1 * Value2, Hybrid_original_list,NewValue1 * NewValue2) :-
convert_into_hybrid(Value1,Hybrid_original_list,NewValue1),
convert_into_hybrid(Value2,Hybrid_original_list,NewValue2).
convert_into_hybrid(Value1 ^ Value2, Hybrid_originals,NewValue1 ^ NewValue2) :-
convert_into_hybrid(Value1,Hybrid_originals,NewValue1),
convert_into_hybrid(Value2,Hybrid_originals,NewValue2).
convert_into_hybrid(Value1 ^ Value2, Hybrid_original_list,NewValue1 ^ NewValue2) :-
convert_into_hybrid(Value1,Hybrid_original_list,NewValue1),
convert_into_hybrid(Value2,Hybrid_original_list,NewValue2).
convert_into_hybrid(Number,_,Number) :-
number(Number),
!.
convert_into_hybrid(Molecule,Hybrid_originals,Molecule) :-
\+(member(Molecule,Hybrid_originals)),
convert_into_hybrid(Molecule,Hybrid_original_list,Molecule) :-
\+(member(Molecule,Hybrid_original_list)),
!.
convert_into_hybrid(Molecule,Hybrid_originals,HybridMolecule) :-
convert_into_hybrid(Molecule,Hybrid_original_list,HybridMolecule) :-
get_option(stochastic_conversion, Rate),
nb_getval(volume,Volume),
member(Molecule,Hybrid_originals),
member(Molecule,Hybrid_original_list),
species_to_total(Molecule, Molecule_total),
atom_concat('(',Molecule_total,Molecule1),
atom_concat(Molecule1,'/(',Molecule2),
......@@ -524,10 +489,10 @@ convert_into_hybrid2(Molecule,ODESpecies,StochMolecule) :-
species_to_stoch(Molecule,StochMolecule),
!.
convert_into_hybrid2(Molecule,Hybrid_originals,HybridMolecule) :-
convert_into_hybrid2(Molecule,Hybrid_original_list,HybridMolecule) :-
get_option(stochastic_conversion, Rate),
nb_getval(volume,Volume),
member(Molecule,Hybrid_originals),
member(Molecule,Hybrid_original_list),
species_to_total(Molecule, Molecule_total),
atom_concat('(',Molecule_total,Molecule1),
atom_concat(Molecule1,'/(',Molecule2),
......@@ -1102,13 +1067,13 @@ concat_species2(Species, Init, Results) :-
atom_concat(Init, First1, Newsubstring),
concat_species(Species, 2, Newsubstring, Results).
print_show_option(Hybrid_originals,Stream) :-
print_show_option(Hybrid_original_list,Stream) :-
(
Hybrid_originals = []
Hybrid_original_list = []
->
write(Stream,'%option(show: {}).'),nl(Stream)
;
concat_species2(Hybrid_originals,'',SpeciesList),
concat_species2(Hybrid_original_list,'',SpeciesList),
write(Stream,'option(show: {'),
write(Stream,SpeciesList),
write(Stream,'}). % Choose here the species to plot'),
......
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