Commit b3853783 authored by HEMERY Mathieu's avatar HEMERY Mathieu
Browse files

Start reflow of hybrid

parent be9b1b66
:- module(hybrid,
[
hybrid_dynamic_simulation/3,
hybrid_static_simulation/2,
hybrid_static_simulation/5,
hybrid_dynamic_simulation/6
]
[
hybrid_static_simulation/2,
hybrid_static_simulation/5,
hybrid_dynamic_simulation/3,
hybrid_dynamic_simulation/6
]
).
% Only for separate compilation/linting
......@@ -27,14 +27,25 @@ It seems like species names like ABC-DFG-{HIJ} make the simulation not work, typ
and concentration is 1. The simulation method currently used is the
Rosenbrock method.').
%! hybrid_static_simulation(+ODEFilename, +StochFilename)
%
% Fast front-end for the hybrid_static_simulation/5 predicate
hybrid_static_simulation(ODEFilename, StochFilename) :-
biocham_command,
doc('This is a hybrid simulation with a static partitioning of the reactions.
The first input file contains the informations concerning the continous
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).
biocham_command,
doc('This is a hybrid simulation with a static partitioning of the reactions.
The first input file contains the informations concerning the continous
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, +OutFilename, +Volume, +Time)
%
% implement the static separation between continuous and stochastic species
% blend the two files in OutFilename and read it to execute a call to the rosenbrock
% integrator.
hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time) :-
biocham_command,
......@@ -46,21 +57,21 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
load_biocham(ODEFilename),
all_items([kind: reaction],Reactions1),
findall(R1,
(
member(Reaction,Reactions1),
reaction(Reaction, [kinetics: _,reactants: Reactants,inhibitors: _, products: _]),
member(Reactant,Reactants),
Reactant = _*R1
),
Reactants1),
(
member(Reaction,Reactions1),
reaction(Reaction, [kinetics: _,reactants: Reactants,inhibitors: _, products: _]),
member(Reactant,Reactants),
Reactant = _*R1
),
Reactants1),
findall(P1,
(
member(Reaction,Reactions1),
reaction(Reaction, [kinetics: _,reactants: _,inhibitors: _, products: Products]),
member(Product,Products),
Product = _*P1
),
Products1),
(
member(Reaction,Reactions1),
reaction(Reaction, [kinetics: _,reactants: _,inhibitors: _, products: Products]),
member(Product,Products),
Product = _*P1
),
Products1),
append(Reactants1, Products1, ODESpecies),
sort(Reactants1,SortedReactants), % Used for later hybrid-reactants finding
sort(ODESpecies,SortedODESpecies),
......@@ -71,26 +82,26 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
%===================================%
% 1. Copy all the parameters.
% 2. Change all "present" into "parameter", with the same initial value.
% Change the species name by appending "_stoch" to the original name
% Change the species name by appending "_stoch" to the original name
% 3. Add two parameters for stochastic simulation:
% parameter (tau, 0).
% parameter (ran, 0).
% parameter (tau, 0).
% parameter (ran, 0).
% 4. Build Macro for the species appearing in both paradigms.
% 5. User-defined volume.
% 6. Build hybrid-reactants list.
% 6. Build hybrid-reactants list.
% process the parameters and presents
load_biocham(StochFilename),
open(OutFilename, write, Stream),
all_items([kind: parameter],Parameters),
findall(Parameter,
(
member(Parameter,Parameters),
write(Stream, Parameter),
write(Stream,'.'),
nl(Stream)
),
_),
(
member(Parameter,Parameters),
write(Stream, Parameter),
write(Stream,'.'),
nl(Stream)
),
_),
write(Stream, 'parameter(tau = 0).'), nl(Stream),
Seed is 0,
write(Stream, 'parameter(ran = '),
......@@ -100,52 +111,52 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
%write(Stream, 'floor(1664525 * ran + 1013904223 / 4294967296)).'),nl(Stream),
all_items([kind:initial_state],InitState),
findall(Hybrid_original,(
(member(present(Hybrid_original),InitState);member(present(Hybrid_original,_),InitState)),
member(Hybrid_original, SortedODESpecies)
), Hybrid_originals),
(member(present(Hybrid_original),InitState);member(present(Hybrid_original,_),InitState)),
member(Hybrid_original, SortedODESpecies)
), Hybrid_originals),
%format("Hybrid_originals is ~w ~n",[Hybrid_originals]),
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_pre, SortedODESpecies)
->
(
atom_concat(Present_pre, '_stoch', Present),
Init_value = 0,
atom_concat(Present_pre, '_total', 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)
)
;
(
atom_concat(Present_pre, '_stoch', Present),
Init_value = Init_value_pre,
atom_concat(Present_pre, '_total', Hybrid_species),
write(Stream, 'function('),
write(Stream, Hybrid_species),
write(Stream,' = '),
write(Stream, Present),
write(Stream, ').'),
nl(Stream)
)
),
write(Stream, 'parameter('), write(Stream, Present),
write(Stream, ' = '), write(Stream, Init_value),
write(Stream, ').'), nl(Stream)
),_Hybrids),
((member(present(Present_pre),InitState),Init_value_pre = 1);member(present(Present_pre,Init_value_pre),InitState)),
(
member(Present_pre, SortedODESpecies)
->
(
atom_concat(Present_pre, '_stoch', Present),
Init_value = 0,
atom_concat(Present_pre, '_total', 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)
)
;
(
atom_concat(Present_pre, '_stoch', Present),
Init_value = Init_value_pre,
atom_concat(Present_pre, '_total', Hybrid_species),
write(Stream, 'function('),
write(Stream, Hybrid_species),
write(Stream,' = '),
write(Stream, Present),
write(Stream, ').'),
nl(Stream)
)
),
write(Stream, 'parameter('), write(Stream, Present),
write(Stream, ' = '), write(Stream, Init_value),
write(Stream, ').'), nl(Stream)
),_Hybrids),
clear_model,
load_biocham(ODEFilename),
......@@ -159,31 +170,31 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
nb_setval(reactions,Reactions2),
findall(R2,
(
member(Reaction,Reactions2),
reaction(Reaction, [kinetics: _,reactants: Reactants,inhibitors: _, products: _]),
member(Reactant,Reactants),
Reactant = _*R2_pre,
atom_concat(R2_pre, '_stoch', R2)
),
Reactants2),
(
member(Reaction,Reactions2),
reaction(Reaction, [kinetics: _,reactants: Reactants,inhibitors: _, products: _]),
member(Reactant,Reactants),
Reactant = _*R2_pre,
atom_concat(R2_pre, '_stoch', 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: _]),
(
Reactants = []
->
Reactant = none
;
member(Reactant,Reactants)
)
),
Reactants_by_rules),
(
member(Reaction,Reactions2),
reaction(Reaction, [kinetics: _,reactants: Reactants,inhibitors: _, products: _]),
(
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),
......@@ -195,14 +206,14 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
% print(Constraints_list), print('\n'),
findall(P2,
(
member(Reaction,Reactions2),
reaction(Reaction, [kinetics: _,reactants: _,inhibitors: _, products: Products]),
member(Product,Products),
Product = _*P2_pre,
atom_concat(P2_pre, '_stoch', P2)
),
Products2),
(
member(Reaction,Reactions2),
reaction(Reaction, [kinetics: _,reactants: _,inhibitors: _, products: Products]),
member(Product,Products),
Product = _*P2_pre,
atom_concat(P2_pre, '_stoch', P2)
),
Products2),
append(Reactants2, Products2, StochSpecies),
sort(StochSpecies,SortedStochSpecies), % To eliminate duplication.
%format("StochSpecies is ~w ~n",[SortedStochSpecies]),
......@@ -219,31 +230,31 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
prepare_change_list(Reactions2,ODESpecies,List),
%format("List is ~w ~n",[List]),
findall([Species_change, Alpha],
(
member(Element,List),
% To get [RSpecies, PSpecies] from [Reactants, Products]
Element = [Reactants,Products,Alpha],
findall([Species1, Minus_number],
(
member(R, Reactants),
R = Minus_number * Rname,
atom_concat(Rname, '_stoch', Species1)
),
RSpecies),
findall([Species2, Add_number],
(
member(P, Products),
P = Add_number * Pname,
atom_concat(Pname, '_stoch', Species2)
),
PSpecies),
%== Build list (Species_change) to record species change ==%
append_reactants(RSpecies, 1, [], Temp_list, SortedStochSpecies),
%format("Temp_list are ~w ~n",[Temp_list]),
append_products(PSpecies, 1, Temp_list, Changes, SortedStochSpecies),
build_change_list(Changes, 1, Species_number, [], Species_change)
),
Reactions_data
(
member(Element,List),
% To get [RSpecies, PSpecies] from [Reactants, Products]
Element = [Reactants,Products,Alpha],
findall([Species1, Minus_number],
(
member(R, Reactants),
R = Minus_number * Rname,
atom_concat(Rname, '_stoch', Species1)
),
RSpecies),
findall([Species2, Add_number],
(
member(P, Products),
P = Add_number * Pname,
atom_concat(Pname, '_stoch', Species2)
),
PSpecies),
%== Build list (Species_change) to record species change ==%
append_reactants(RSpecies, 1, [], Temp_list, SortedStochSpecies),
%format("Temp_list are ~w ~n",[Temp_list]),
append_products(PSpecies, 1, Temp_list, Changes, SortedStochSpecies),
build_change_list(Changes, 1, Species_number, [], Species_change)
),
Reactions_data
),
%format("Reactions_data are ~w ~n",[Reactions_data]),
......@@ -260,7 +271,7 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
%write(StreamAppend, 'show_macros.'), nl(StreamAppend),
write(StreamAppend,'% Plot part'),nl(StreamAppend),nl(StreamAppend),
format("Event printed ~w ~n",[Hybrid_originals]),
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),
......@@ -293,40 +304,40 @@ build_constraints_by_rules(Reactants_by_rules, ODEReactants, Current_counter, Or
nth1(Current_counter, Reactants_by_rules, This_reactants_list),
%format("This_reactants_list is ~w ~n",[This_reactants_list]),
findall(Constraints,(
(
\+(list(This_reactants_list))
->
This_reactants_list = Least_number * Reactant_name,
member(Reactant_name, ODEReactants)
;
member(Reactants_pair, This_reactants_list),
Reactants_pair = Least_number * Reactant_name,
% print('ode_reactants: '), print(ODEReactants), print('\n'),
member(Reactant_name, ODEReactants)
)
->
(
atom_concat(Reactant_name, '_total >= ', Temp),
Reactants_pair = Least_number * Reactant_name,
atom_number(Least_number_atom,Least_number),
atom_concat(Temp, Least_number_atom, Constraints)
)
;
Constraints = 0
(
\+(list(This_reactants_list))
->
This_reactants_list = Least_number * Reactant_name,
member(Reactant_name, ODEReactants)
;
member(Reactants_pair, This_reactants_list),
Reactants_pair = Least_number * Reactant_name,
% print('ode_reactants: '), print(ODEReactants), print('\n'),
member(Reactant_name, ODEReactants)
)
->
(
atom_concat(Reactant_name, '_total >= ', Temp),
Reactants_pair = Least_number * Reactant_name,
atom_number(Least_number_atom,Least_number),
atom_concat(Temp, Least_number_atom, Constraints)
)
;
Constraints = 0
),This_constraints_list),
% print(This_constraints_list),print('###\n'),
% print(This_constraints_list),print('###\n'),
append(Orig_list, [This_constraints_list], Next_list),
(
Current_counter < Number_of_rules
->
(
Next_counter is Current_counter + 1,
build_constraints_by_rules(Reactants_by_rules, ODEReactants, Next_counter, Next_list, Results_list)
Next_counter is Current_counter + 1,
build_constraints_by_rules(Reactants_by_rules, ODEReactants, Next_counter, Next_list, Results_list)
)
;
(
Results_list = Next_list,
!
Results_list = Next_list,
!
)
).
......@@ -336,13 +347,13 @@ build_constraints_by_rules(Reactants_by_rules, ODEReactants, Current_counter, Or
This_list \== [0]
->
(
atom_number(Reaction_counter_atom,Reaction_counter),
atom_concat('lower_level', Reaction_counter_atom, Corresponding_constraint),
append(Init_list, [Corresponding_constraint], Next_list)
atom_number(Reaction_counter_atom,Reaction_counter),
atom_concat('lower_level', Reaction_counter_atom, Corresponding_constraint),
append(Init_list, [Corresponding_constraint], Next_list)
)
;
(
append(Init_list, ['0'], Next_list)
append(Init_list, ['0'], Next_list)
)
),
(
......@@ -350,8 +361,8 @@ build_constraints_by_rules(Reactants_by_rules, ODEReactants, Current_counter, Or
Reaction_counter < Total_count
->
(
Next_counter is Reaction_counter + 1,
build_constraints_list(Hybrid_reactants_list, Next_counter, Next_list, Constraints_list)
Next_counter is Reaction_counter + 1,
build_constraints_list(Hybrid_reactants_list, Next_counter, Next_list, Constraints_list)
)
;
Constraints_list = Next_list
......@@ -360,20 +371,20 @@ build_constraints_by_rules(Reactants_by_rules, ODEReactants, Current_counter, Or
build_reactants(StochSpecies, ODESpecies, Species_number, Counter, Init, StochReactants) :-
nth1(Counter, StochSpecies, This_species),
(
member(This_species, ODESpecies)
member(This_species, ODESpecies)
->
atom_concat(This_species, '_total', This_reactant)
atom_concat(This_species, '_total', This_reactant)
;
atom_concat(This_species, '_stoch', This_reactant)
atom_concat(This_species, '_stoch', This_reactant)
),
append(Init, [This_reactant], Next),
(
Counter < Species_number
Counter < Species_number
->
Next_counter is Counter + 1,
build_reactants(StochSpecies, ODESpecies, Species_number, Next_counter, Next, StochReactants)
Next_counter is Counter + 1,
build_reactants(StochSpecies, ODESpecies, Species_number, Next_counter, Next, StochReactants)
;
StochReactants = Next
StochReactants = Next
).
write_hybrid_ode(Stream,Hybrid_originals) :-
......@@ -389,43 +400,43 @@ write_hybrid_ode(Stream,Hybrid_originals) :-
list_hybrid_ode(Stream,Hybrid_originals) :- % Only MA are converted, need to do others
all_items([kind: reaction], Reactions),
forall(
member(Reaction,Reactions),
(
Reaction =.. [Op,_,_],
Op = (=>)
->
convert_MA(Reaction,Reactions,Hybrid_originals,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),
NewReaction =.. [for,NewKinetics,SecondPart],
write(Stream,NewReaction),write(Stream,'.'),nl(Stream)
)
member(Reaction,Reactions),
(
Reaction =.. [Op,_,_],
Op = (=>)
->
convert_MA(Reaction,Reactions,Hybrid_originals,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),
NewReaction =.. [for,NewKinetics,SecondPart],
write(Stream,NewReaction),write(Stream,'.'),nl(Stream)
)
).
/* 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 */
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) :-
reaction(Reaction, [kinetics: Kinetics,reactants: Reactants,inhibitors: _, products: _]),
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),!,
(
nb_current(dynamic_parameters,_)
nb_current(dynamic_parameters,_)
->
nth1(ReactionNumber,Reactions,Reaction),
atom_concat('k',ReactionNumber,Kdiff1),
atom_concat(Kdiff1,'_diff',Kdiff),
NewPartialKinetics = K * NewValue,
assertz(reaction_kinetics(Reaction,NewPartialKinetics)),
NewKinetics = Kdiff * NewPartialKinetics
nth1(ReactionNumber,Reactions,Reaction),
atom_concat('k',ReactionNumber,Kdiff1),
atom_concat(Kdiff1,'_diff',Kdiff),
NewPartialKinetics = K * NewValue,
assertz(reaction_kinetics(Reaction,NewPartialKinetics)),
NewKinetics = Kdiff * NewPartialKinetics
;
NewKinetics = K * NewValue
NewKinetics = K * NewValue
).
%format("NewKinetics is ~w ~n",[NewKinetics]).
......@@ -466,25 +477,25 @@ convert_into_hybrid(Molecule,Hybrid_originals,HybridMolecule) :-
prepare_change_list(Reactions2,ODESpecies,List) :- % Only MA are converted, need to do others
findall([Reactants,Products,Kinetics],
(
member(Reaction,Reactions2),
convert_MA2(Reaction,ODESpecies,Kinetics,Reactants,Products)
),
List
(
member(Reaction,Reactions2),
convert_MA2(Reaction,ODESpecies,Kinetics,Reactants,Products)
),
List
).
convert_MA2(Reaction,ODESpecies,NewKinetics,Reactants,Products) :-
reaction(Reaction, [kinetics: Kinetics,reactants: Reactants,inhibitors: _, products: Products]),
reaction(Reaction, [kinetics: Kinetics,reactants: Reactants,inhibitors: _, products: Products]),
Kinetics =.. ['MA',K],
kinetics:eval_kinetics(Reactants, _, product('S'*'M' in [reactants], 'M'^'S'), Value),
convert_into_hybrid2(Value,ODESpecies,NewValue),!,
NewKinetics = K * NewValue,
(
nb_current(dynamic_parameters,_)
nb_current(dynamic_parameters,_)
->
true
true
;
assertz(reaction_kinetics(Reaction,NewKinetics))
assertz(reaction_kinetics(Reaction,NewKinetics))
).
%format("NewKinetics is ~w ~n",[NewKinetics]).
......@@ -528,27 +539,27 @@ build_change_list(Changes, Counter, Length, Init, Result) :-
Counter =< Length
->
(
(
\+ member([Counter, _], Changes)
->
(
\+ member([Counter, _], Changes)
->
(
atom_number(Zero,0),
atom_concat('+', Zero, No_change),
append(Init, [No_change], Temp_list)
)
;
(
findall(N, member([Counter, N], Changes), N_list),
summation(N_list, 1, '', Sum),
append(Init, [Sum], Temp_list)
)
atom_number(Zero,0),
atom_concat('+', Zero, No_change),
append(Init, [No_change], Temp_list)
)
;
(
findall(N, member([Counter, N], Changes), N_list),
summation(N_list, 1, '', Sum),
append(Init, [Sum], Temp_list)
)
),
Next_counter is (Counter + 1),
build_change_list(Changes, Next_counter, Length, Temp_list, Result)
)
;
(
Result = Init
Result = Init
).
summation(Input_list, Count, Init, Final) :-
......@@ -559,12 +570,12 @@ summation(Input_list, Count, Init, Final) :-