Reflow print_alpha_parameters

hybrid.pl

@@ -210,15 +210,13 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
   sort(Species_list, Sorted_Species), % To eliminate duplication.
   append(Reactant_stoch_list, Product_stoch_list, StochSpecies),
   sort(StochSpecies, SortedStochSpecies), % To eliminate duplication.
-  % build StochReactants(with '_total' or '_stoch') from Sorted_Species.
-  build_reactants(Sorted_Species, SortedODESpecies, StochReactants),
 
   %===================================================================%
   %==== Build Data Lists for Reactions and print events in Stream ====%
   %===================================================================%
   prepare_change_list(Reactions2,ODESpecies,List),
   build_reaction_data(List, Sorted_Species, Reaction_data_list),
-  print_event(Reaction_data_list, SortedStochSpecies, StochReactants,Hybrid_reactants_list, Constraints_list, Stream),
+  print_event(Reaction_data_list, SortedStochSpecies, Hybrid_reactants_list, Constraints_list, Stream),
   
   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",[]),
@@ -459,17 +457,18 @@ convert_into_hybrid(Molecule, _List_species, _Mode, HybridMolecule) :-
   !.
 
 
-%! prepare_change_list
-%
+%! prepare_change_list(+Reaction_list, +Species_list, -Change_list)
 %
+% Extract the reactant, product and kinetics of the reactions
+% TODO: Only MA are converted, need to do others
 
-prepare_change_list(Reactions2,ODESpecies,List) :- % Only MA are converted, need to do others
+prepare_change_list(Reaction_list, Species_list, Change_list) :- 
   findall([Reactants,Products,Kinetics],
          (
-           member(Reaction,Reactions2),
-           convert_MA2(Reaction,ODESpecies,Kinetics,Reactants,Products)
+           member(Reaction, Reaction_list),
+           convert_MA2(Reaction, Species_list, Kinetics, Reactants, Products)
          ),
-         List
+         Change_list
   ).
 
 
@@ -486,7 +485,7 @@ build_stoechiometry(Reactant_list, Product_list, [Sp|TailSp], [Stoech|TailStoech
 	sum_list(CR_list, CRT),
 	sum_list(CP_list, CPT),
 	Stoech_number is CPT-CRT,
-  (
+  ( % Explicitely include the sign of the stoechiometric number
     Stoech_number < 0
   ->
     atom_number(Stoech, Stoech_number)
@@ -512,64 +511,48 @@ build_reaction_data([[RList, PList, Coef]|ReacTail], Species_list, [[Stoech, Coe
 %
 % write the events in the output file
 
-print_event(Reactions_data, Stoch_species_order, _Stoch_reactants_order, Hybrid_reactants_list, Constraints_list, Stream) :- 
+print_event(Reactions_data, Stoch_list, Hybrid_reactants_list, Constraints_list, Stream) :- 
   get_option(stochastic_conversion,Rate), 
   length(Reactions_data, Number_of_alpha),
-  print_alpha_parameters(Number_of_alpha, 1, Stream),
+  length(Stoch_list, Number_of_species),
+  print_alpha_parameters(Number_of_alpha, Stream),
   print_constraints_parameters(Constraints_list, 1, Stream),
-% 1. Build Starting String for the event.
-  write(Stream, 'add_event(Time > tau,'),nl(Stream),
-  print_static_alpha(Number_of_alpha,Reactions_data, 1, Stream),
-  write(Stream,' tau = (if alpha_sum <= 0 then inf else Time+((-1 / alpha_sum) * log(random_float))/'),write(Stream,Rate),write(Stream,'),'),nl(Stream),
-  write(Stream,' ran = random_float,'),nl(Stream),
+  
+  % 1. Build Starting String for the event.
+  format(Stream, "add_event(Time > tau,~n",[]),
+  print_static_alpha(Number_of_alpha, Reactions_data, 1, Stream),
+  format(Stream,
+    "tau = (if alpha_sum <= 0 then inf else Time+((-1 / alpha_sum) * log(random_float))/~w),~n",
+    [Rate]),
+  format(Stream,"ran = random_float,~n",[]),
   print_hybrid_species_constraint(Hybrid_reactants_list, 1, Stream),
 
 % 2. Build alpha list and alpha subsums.
-  length(Reactions_data, Number_of_alpha),
-
-  % Build Alpha_list
   build_alpha_list(1, [], Number_of_alpha, Alpha_list),
-
-  % Build Alpha subsums
   nth1(1,Alpha_list, First_alpha),
   append([], [First_alpha], Init_list),
   alpha_accumulation(Alpha_list, Init_list, Alpha_accumulation_list, 2),
 
-
 % 3. Print reactions in event
-  length(Stoch_species_order, Number_of_species),
   constant_list(Number_of_species,[],Init_list2),
   % Build the print_helper_list.
-  print_normal(Reactions_data, Init_list2, Print_helper_list, Alpha_list, Alpha_accumulation_list, Constraints_list, Stoch_species_order),
-
+  print_normal(Reactions_data, Init_list2, Print_helper_list, Alpha_list, Alpha_accumulation_list, Constraints_list, Stoch_list),
   % Use print_helper_list to really print out the event part.
-  length(Stoch_species_order,NumberOfSpecies),
-  print_change_list(NumberOfSpecies,Print_helper_list, Stoch_species_order, 1, Number_of_species, Stream).
-  /*
-  % 4. Print Macros
-%  nth1(Number_of_alpha, Alpha_accumulation_list, Alpha_constitution),
-  write(Stream, 'macro(alpha_summ'),
-  write(Stream, ', '),
-  print_alpha_sum(Number_of_alpha, 1, Stream),
-%  write(Stream, Alpha_constitution),
-  write(Stream, ').'),
-  nl(Stream),
-  write(Stream, 'macro(delta_t'),
-  write(Stream, ', (-1 / alpha_summ) * log(ran)).'),
-  nl(Stream).*/
+  print_change_list(Number_of_species,Print_helper_list, Stoch_list, 1, Number_of_species, Stream).
+
+
+%! print_alpha_parameters(+N_alpha, +Stream)
+%
+%
+
+print_alpha_parameters(0, Stream) :-
+  format(Stream, "parameter(alpha_sum = 0).~n", []),!.
+
+print_alpha_parameters(N_alpha, Stream) :-
+  format(Stream, "parameter(alpha~d = 0).~n", [N_alpha]),
+  Nm_alpha is N_alpha - 1,
+  print_alpha_parameters(Nm_alpha, Stream).
 
-print_alpha_parameters(Number_of_alpha, Current_count, Stream) :-
-  write(Stream, 'parameter(alpha'),
-  write(Stream, Current_count),
-  write(Stream, ' = 0).\n'),
-  (
-   Current_count < Number_of_alpha
-  ->
-   Next_count is Current_count + 1,
-   print_alpha_parameters(Number_of_alpha, Next_count, Stream)
-  ;
-   write(Stream, 'parameter(alpha_sum = 0).\n')
-  ).
 
 print_constraints_parameters(Constraints_list, Current_counter, Stream) :-
   length(Constraints_list, Total_count),
@@ -1274,7 +1257,7 @@ print_dynamic_event_condition_ifThenElse([],_,_).
 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),
-  print_alpha_parameters(Number_of_alpha, 1, Stream),
+  print_alpha_parameters(Number_of_alpha, Stream),
   print_constraints_parameters(Constraints_list, 1, Stream),
 % 1. Build Starting String for the event.
   write(Stream, 'add_event(Time > tau,'),nl(Stream),