Reflow of hybrid_dynamic_simulation

hybrid.pl

@@ -4,7 +4,7 @@
 * continuous).
 * The main purpose of his module is thus to write a .bc file describing the simulation.
 *
-* Coders: P. Remondeau
+* Coders: P. Remondeau, M. Hemery
 */
 :- module(hybrid,
   [
@@ -1090,7 +1090,6 @@ hybrid_dynamic_simulation(InputFile,PropTresh,PartTresh) :-
 
 hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh) :-
   biocham_command,
-  %get_option(stochastic_conversion, Rate)
   nb_setval(volume,Volume),
   load_biocham(InputFile),
   open(OutFilename,write,Stream),
@@ -1115,21 +1114,19 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
   append(Reactants1, Products1, Species1),
   sort(Reactants1,SortedReactants), % Used for later hybrid-reactants finding and constraints lists 
   sort(Species1,SortedSpecies),
+  maplist(species_to_stoch, SortedSpecies, SortedStochSpecies),
 
   format(Stream, "% Hybrid version of ~w~n", InputFile),
-  write(Stream, 'parameter(tau = 0).'), nl(Stream), 
+  format(Stream, "parameter(tau = 0).~n", []),
   Seed is 0,
-  format(Stream, "seed(~d).~n", [Seed]),
-  write(Stream, 'parameter(ran = '),
-  write(Stream, Seed),
-  write(Stream, ').'), nl(Stream),
-  write(Stream,'parameter(delta_t = 0).'),nl(Stream),nl(Stream),
-
+  format(Stream, "seed(~d).~nparameter(ran = ~d).~n", [Seed, Seed]),
+  format(Stream,"parameter(delta_t = 0).~n~n",[]),
 
   set_counter(dynamic_parameters,1), % Will count the number of ki_diff parameters
   length(Reactions,NumberOfReactions),
-  print_dynamic_parameters(NumberOfReactions,Stream),
-  nl(Stream),nl(Stream),
+  numlist(1, NumberOfReactions, Numbered_reac),
+  maplist(format(Stream, "parameter(k~g_diff = 0).~n"), Numbered_reac),
+  format(Stream, "~n~n", []),
 
   all_items([kind:initial_state],InitState),
   get_option(stochastic_conversion, Rate),
@@ -1176,57 +1173,31 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
              )
            ),
           Reactants_by_rules),
-
-  findall(StochSpecie,
-           (
-             member(Species,SortedSpecies),
-             species_to_stoch(Species,StochSpecie)
-           ),
-          SortedStochSpecies),
   
   build_constraints_by_rules(Reactants_by_rules, SortedReactants, 1, [], Hybrid_reactants_list),  
   build_constraints_list(Hybrid_reactants_list, Constraints_list),
-  
   build_reactants(SortedSpecies, SortedSpecies, StochReactants),
 
   prepare_change_list(Reactions, SortedSpecies, Reaction_list),
   build_reaction_data(Reaction_list, SortedSpecies, Reactions_data),
+  build_hybrid_reactants(Reaction_list, HybridDynamicReactants),
   compute_maximum_change(Reaction_list, MaxParticleChange),
   ParticuleTreshold is MaxParticleChange*PartTresh,
+
+  print_dynamic_event(SortedSpecies,Reactions_data, HybridDynamicReactants,ParticuleTreshold,PropTresh, SortedStochSpecies, StochReactants,Hybrid_reactants_list, Constraints_list, Stream),
+  retractall(reaction_kinetics(_,_)),
   
-  findall([RSpecies,Alpha],
-        (
-           member(Element,Reaction_list),
-           Element = [Reactants,_,Alpha],
-           findall(Rname1,
-                  (
-                     member(R, Reactants),
-                     R = _ * Rname,
-                     species_to_total(Rname, Rname1)
-                   ),
-                   RSpecies)
-         ),
-         HybridDynamicReactants),
-           
+  format(Stream, "~n% Plot part~n~n", []),
+  print_show_option(SortedSpecies,Stream),
+  format(Stream, "numerical_simulation(method:rsbk, time: ~g).~n", [Time]),
+  format(Stream, "plot.~noption(show: {}).~n",[]),
   close(Stream),
-  open(OutFilename,append,StreamAppend),
-  print_dynamic_event(SortedSpecies,Reactions_data, HybridDynamicReactants,ParticuleTreshold,PropTresh, SortedStochSpecies, StochReactants,Hybrid_reactants_list, Constraints_list, StreamAppend),
-  nb_delete(dynamic_parameters),
-  nb_delete(reactions),
-  retractall(reaction_kinetics(_,_)),
-  nl(StreamAppend),
   
-  write(StreamAppend,'% Plot part'),nl(StreamAppend),nl(StreamAppend),
-  print_show_option(SortedSpecies,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),
   nb_delete(volume),
+  nb_delete(dynamic_parameters),
+  nb_delete(reactions),
+
+  % Load the new file and execute the simulation
   load_biocham(OutFilename),
   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",[]).
 
@@ -1237,6 +1208,8 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
 
 
 %! compute_maximum_change(+Reaction_list, -MaxParticleChange) :-
+%
+% Determine the maximum change that a reaction may provoke in the number of particles
 
 compute_maximum_change(Reaction_list, MaxParticleChange) :-
   extract_stoechiometry_list(Reaction_list, Reactant_list, Product_list),
@@ -1250,6 +1223,7 @@ compute_maximum_change(Reaction_list, MaxParticleChange) :-
 %! extract_stoechiometry_list(+Reaction_list, -Reactant_list, -Product_list)
 %
 % Extract a list of the total stoechiometry of each reaction 
+% to be used in compute_maximum_change/2
 
 extract_stoechiometry_list([],[],[]) :- !.
 
@@ -1260,13 +1234,20 @@ extract_stoechiometry_list([[Reactants, Products, _A]|TReac], [SR|TReactant], [S
   sum_list(CoeffP_list, SP),
   extract_stoechiometry_list(TReac, TReactant, TProduct).
 
+% build_hybrid_reactants(+Reaction_list, -HybridDynamicReactants)
+
+build_hybrid_reactants([], []) :- !.
 
-particle_sum_list([],0).
+build_hybrid_reactants([[Reactants, _Products, Alpha]|ReacTail], [[ReacName, Alpha]|OutTail]) :-
+  build_hybrid_subroutine(Reactants, ReacName),
+  build_hybrid_reactants(ReacTail, OutTail).
+
+build_hybrid_subroutine([], []) :- !.
+
+build_hybrid_subroutine([_N*X|InputTail], [Name|OutputTail]) :-
+  species_to_total(X, Name),
+  build_hybrid_subroutine(InputTail, OutputTail).
 
-particle_sum_list([Element|List],Sum) :-
-  nth1(2,Element,TrueElement),
-  particle_sum_list(List,PartialSum),
-  Sum = PartialSum + TrueElement.
 
 print_dynamic_event_condition(SortedSpecies,HybridDynamicReactants,ReactionNumber,PropTresh,MaxParticleChange,Stream) :-
   nth1(ReactionNumber,HybridDynamicReactants,ReactionInfo),
@@ -1301,21 +1282,9 @@ print_dynamic_event_condition_ifThenElse([Reactant|Reactants],MaxParticleChange,
   write(Stream,MaxParticleChange),
   print_dynamic_event_condition_ifThenElse(Reactants,MaxParticleChange,Stream).
 
-
 print_dynamic_event_condition_ifThenElse([],_,_).
 
 
-print_dynamic_parameters(NumberOfReactions,Stream) :-
-  repeat,
-  count(dynamic_parameters,Counter),
-  write(Stream,'parameter(k'),
-  write(Stream,Counter),
-  write(Stream,'_diff = 0).'),
-  nl(Stream),
-  NextCounter is Counter + 1,
-  NextCounter > NumberOfReactions,
-  !.
-
 print_dynamic_event(SortedSpecies,Reactions_data, HybridDynamicReactants,MaxParticleChange,PropTresh,Stoch_species_order, _Stoch_reactants_order, Hybrid_reactants_list, Constraints_list, Stream) :-  
   get_option(stochastic_conversion,Rate),
   length(Reactions_data, Number_of_alpha),