Commit 583f9519 authored by Mathieu Hemery's avatar Mathieu Hemery
Browse files

Cleaning the code of verbous juttering

parent e0a90cda
......@@ -83,7 +83,6 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
append(Reactants1, Products1, ODESpecies),
sort(Reactants1,SortedReactants), % Used for later hybrid-reactants finding
sort(ODESpecies,SortedODESpecies),
%format("SortedODESpecies is ~w ~n",[SortedODESpecies]),
clear_model,
%===================================%
%== Preprocessing Stochastic file ==%
......@@ -123,7 +122,6 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
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)),
......@@ -206,12 +204,7 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
format("Reactants_by_rules is ~w ~n",[Reactants_by_rules]),
build_constraints_by_rules(Reactants_by_rules, SortedReactants, 1, [], Hybrid_reactants_list),
% print('Hybrid_reactants_list: '),
% print(Hybrid_reactants_list), print('\n'),
build_constraints_list(Hybrid_reactants_list, 1, [], Constraints_list),
%format("Constraints_list is ~w ~n",[Constraints_list]),
% print('out from build_constraints_list\n'),
% print(Constraints_list), print('\n'),
findall(P2,
(
......@@ -224,19 +217,15 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
Products2),
append(Reactants2, Products2, StochSpecies),
sort(StochSpecies,SortedStochSpecies), % To eliminate duplication.
%format("StochSpecies is ~w ~n",[SortedStochSpecies]),
% build StochReactants(with '_total') from SortedStochSpecies.
length(SortedStochSpecies, Species_number),
build_reactants(SortedStochSpecies, SortedODESpecies, Species_number, 1, [], StochReactants),
%format("StochReactants is ~w ~n",[StochReactants]),
%========================================%
%==== Build Data Lists for Reactions ====%
%========================================%
prepare_change_list(Reactions2,ODESpecies,List),
%format("List is ~w ~n",[List]),
findall([Species_change, Alpha],
(
member(Element,List),
......@@ -258,18 +247,13 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
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]),
% print('Change list: '), print(Reactions_data), print('\n'),
%== Build Events and Calaculate Alphas==%
%== Build Events and Calaculate Alphas==%
close(Stream),
open(OutFilename, append, StreamAppend),
format("Hybrid_reactants_list is ~w ~n",[Hybrid_reactants_list]),
......@@ -277,7 +261,6 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
nl(StreamAppend),
%print_common_alpha_macros(Reactions_data, 1, StreamAppend),
%write(StreamAppend, 'show_macros.'), nl(StreamAppend),
write(StreamAppend,'% Plot part'),nl(StreamAppend),nl(StreamAppend),
format("Event printed ~w ~n",[Hybrid_originals]),
%print_show_option(Hybrid_originals,StreamAppend),
......@@ -325,7 +308,6 @@ species_to_total(Species, Species_total) :-
build_constraints_by_rules(Reactants_by_rules, ODEReactants, Current_counter, Orig_list, Results_list) :-
length(Reactants_by_rules, Number_of_rules),
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))
......@@ -335,7 +317,6 @@ build_constraints_by_rules(Reactants_by_rules, ODEReactants, Current_counter, Or
;
member(Reactants_pair, This_reactants_list),
Reactants_pair = Least_number * Reactant_name,
% print('ode_reactants: '), print(ODEReactants), print('\n'),
member(Reactant_name, ODEReactants)
)
->
......@@ -349,7 +330,6 @@ build_constraints_by_rules(Reactants_by_rules, ODEReactants, Current_counter, Or
;
Constraints = 0
),This_constraints_list),
% print(This_constraints_list),print('###\n'),
append(Orig_list, [This_constraints_list], Next_list),
(
Current_counter < Number_of_rules
......@@ -413,13 +393,10 @@ build_reactants(StochSpecies, ODESpecies, Species_number, Counter, Init, StochRe
write_hybrid_ode(Stream,Hybrid_originals) :-
list_hybrid_ode(Stream,Hybrid_originals),
%list_model_influences,
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_model_macros(ltl_pattern),
%list_model_options.
list_hybrid_ode(Stream,Hybrid_originals) :- % Only MA are converted, need to do others
all_items([kind: reaction], Reactions),
......@@ -462,7 +439,6 @@ convert_MA(Reaction,Reactions,Hybrid_originals,NewKinetics) :-
;
NewKinetics = K * NewValue
).
%format("NewKinetics is ~w ~n",[NewKinetics]).
convert_into_hybrid(Value * 1, ODESpecies,NewValue) :-
convert_into_hybrid(Value,ODESpecies,NewValue).
......@@ -522,7 +498,6 @@ convert_MA2(Reaction,ODESpecies,NewKinetics,Reactants,Products) :-
;
assertz(reaction_kinetics(Reaction,NewKinetics))
).
%format("NewKinetics is ~w ~n",[NewKinetics]).
convert_into_hybrid2(Value * 1, ODESpecies,NewValue) :-
convert_into_hybrid2(Value,ODESpecies,NewValue).
......@@ -612,23 +587,11 @@ print_event(Reactions_data, Stoch_species_order, _Stoch_reactants_order, Hybrid_
print_constraints_parameters(Constraints_list, 1, Stream),
% 1. Build Starting String for the event.
write(Stream, 'add_event(Time > tau,'),nl(Stream),
%%% print('Checkpoint 1\n'),
print_static_alpha(Number_of_alpha,Reactions_data, 1, Stream),
format("Checkpoint ~n",[]),
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),
%print_event_constraints(Constraints_list, 1, Stream),
%concat_species(Stoch_species_order, 1, '', Species_for_event),
%write(Stream, Species_for_event),
%write(Stream, '], [rand, Time + delta_t, rand'),
%%% print('Checkpoint 2\n'),
%print_updating_alpha(Number_of_alpha, 1, Stream),
%%% print('Checkpoint 3\n'),
print_hybrid_species_constraint(Hybrid_reactants_list, 1, Stream),
%%% print('Checkpoint 4\n'),
% 2. Build alpha list and alpha subsums.
length(Reactions_data, Number_of_alpha),
......@@ -639,7 +602,6 @@ print_event(Reactions_data, Stoch_species_order, _Stoch_reactants_order, Hybrid_
nth1(1,Alpha_list, First_alpha),
append([], [First_alpha], Init_list),
alpha_accumulation(Alpha_list, Init_list, Alpha_accumulation_list, 2),
%format("Alpha_accumulation_list is ~w ~n",[Alpha_accumulation_list]),
% 3. Print reactions in event
......@@ -647,14 +609,10 @@ print_event(Reactions_data, Stoch_species_order, _Stoch_reactants_order, Hybrid_
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),
format("Checkpoint ~n",[]),
%format("Print_helper_list is ~w ~n",[Print_helper_list]),
%%% print('Checkpoint 5\n'),
% 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).
%%% print('Checkpoint 6\n'),
/*
% 4. Print Macros
% nth1(Number_of_alpha, Alpha_accumulation_list, Alpha_constitution),
......@@ -711,8 +669,6 @@ print_static_alpha(Number_of_alpha,Reactions_data, Current_count, Stream) :-
nb_getval(reactions,Reactions),
nth1(Current_count,Reactions,Reaction),
reaction_kinetics(Reaction,NewKinetics),
format("Reaction ~w is ~w Kinetics is ~w ",[Current_count,Reaction,NewKinetics]),
(
nb_current(dynamic_parameters,_)
->
......@@ -857,7 +813,6 @@ print_normal(Reactions_data, Init_list, N_finish_list, Alpha_list, Alpha_accumul
% Recursive call for all Normal reactions.
%format("Constraints_list is ~w ~n",[Constraints_list]),
print_normal_middle(Reactions_data, Later_list, N_finish_list,Alpha_list, Alpha_accumulation_list, Constraints_list,Stoch_species_order, 2).
check_change1(First_alpha, Constraint, Stoch_species_order, Species_counter, First_species_change, Orig_list, Later_list) :-
......@@ -928,9 +883,7 @@ print_normal_middle(Reactions_data, Orig_list, Result_list, Alpha_list, Alpha_ac
nth1(Previous_reaction, Alpha_accumulation_list, First_alpha),
nth1(Counter_for_reaction, Alpha_accumulation_list, Second_alpha),
check_change(First_alpha, Second_alpha, Constraint, Stoch_species_order,1, Species_change, Orig_list, Middle_list),
format("Checkpoint ~w for ~w steps ~n",[Counter_for_reaction,Last_reaction]),
%%% print('middle list: '), print(Middle_list), print('\n'),
% Recursive condition.
(
length(Reactions_data, Last_reaction),
......@@ -940,9 +893,7 @@ print_normal_middle(Reactions_data, Orig_list, Result_list, Alpha_list, Alpha_ac
Next_reaction is (Counter_for_reaction + 1),
print_normal_middle(Reactions_data, Middle_list, Result_list,Alpha_list, Alpha_accumulation_list, Constraints_list,Stoch_species_order, Next_reaction)
;
format("FinalStep ~n",[]),
Result_list = Middle_list
%%% print(Result_list), print('\n')
).
check_change(First_alpha, Second_alpha, Constraint, Stoch_species_order,Species_counter, Species_change, Orig_list, Result_list) :-
......@@ -1009,7 +960,7 @@ print_change_list(NumberOfSpecies,Print_helper_list, Species_order, Current_coun
write(Stream, Constraint),
write(Stream, '>0')
;
debug(debug, "No constraint for this\n", [])
debug(hybrid, "No constraint for this\n", [])
),
write(Stream, ' then '),
......@@ -1106,7 +1057,6 @@ print_hybrid_species_constraint(Hybrid_reactants_list, Current_counter, Stream)
length(Hybrid_reactants_list, Total_count),
nth1(Current_counter, Hybrid_reactants_list, This_reaction_constraints),
length(This_reaction_constraints, Number_of_constraints),
%format("This_reaction_constraints is ~w ~n",[This_reaction_constraints]),
(
This_reaction_constraints \== [0]
->
......@@ -1232,7 +1182,6 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
length(Reactions,NumberOfReactions),
print_dynamic_parameters(NumberOfReactions,Stream),
nl(Stream),nl(Stream),
format("Checkpoint ~w ~n",[SortedSpecies]),
all_items([kind:initial_state],InitState),
get_option(stochastic_conversion, Rate),
......@@ -1291,18 +1240,14 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
species_to_stoch(Species,StochSpecie)
),
SortedStochSpecies),
%format("SortedStochSpecies is ~w ~n",[SortedStochSpecies]),
build_constraints_by_rules(Reactants_by_rules, SortedReactants, 1, [], Hybrid_reactants_list),
% print('Hybrid_reactants_list: '),
% print(Hybrid_reactants_list), print('\n'),
build_constraints_list(Hybrid_reactants_list, 1, [], Constraints_list),
length(SortedSpecies,Species_number),
build_reactants(SortedSpecies, SortedSpecies, Species_number, 1, [], StochReactants),
prepare_change_list(Reactions,SortedSpecies,List),
%format("List is ~w ~n",[List]),
% Because inf/0 is not defined in old versions of SWI-Prolog, we use
% something else
current_prolog_flag(max_tagged_integer, Inf),
......@@ -1315,7 +1260,6 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
member(Element,List),
% To get [RSpecies, PSpecies] from [Reactants, Products]
Element = [Reactants,Products,Alpha],
%format("Reactants are ~w Products are ~w ~n",[Reactants,Products]),
findall([Rname, Minus_number],
(
member(R, Reactants),
......@@ -1354,8 +1298,6 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
compute_maximum_change(MaxParticleChange),
ParticuleTreshold is MaxParticleChange*PartTresh,
format("Hybrid_dynamic_reactants is ~w ~n",[HybridDynamicReactants]),
%format("Reactions_data is ~w ~n",[Reactions_data]),
close(Stream),
open(OutFilename,append,StreamAppend),
print_dynamic_event(SortedSpecies,Reactions_data, HybridDynamicReactants,ParticuleTreshold,PropTresh, SortedStochSpecies, StochReactants,Hybrid_reactants_list, Constraints_list, StreamAppend),
......@@ -1429,7 +1371,6 @@ compute_maximum_change(MaxParticleChange) :-
particle_sum_list([],0).
particle_sum_list([Element|List],Sum) :-
%format("Element is ~w ~n",[Element]),
nth1(2,Element,TrueElement),
particle_sum_list(List,PartialSum),
Sum = PartialSum + TrueElement.
......@@ -1452,7 +1393,6 @@ print_dynamic_event_condition(SortedSpecies,HybridDynamicReactants,ReactionNumbe
nl(Stream),
(
length(HybridDynamicReactants,NumberOfReactions),
format("ReactionNumber is ~w NumberOfReactions is ~w ~n",[ReactionNumber,NumberOfReactions]),
ReactionNumber < NumberOfReactions
->
NextReactionNumber is ReactionNumber + 1,
......@@ -1490,26 +1430,14 @@ print_dynamic_event(SortedSpecies,Reactions_data, HybridDynamicReactants,MaxPart
print_constraints_parameters(Constraints_list, 1, Stream),
% 1. Build Starting String for the event.
write(Stream, 'add_event(Time > tau,'),nl(Stream),
%%% print('Checkpoint 1\n'),
format("Checkpoint ~n",[]),
print_dynamic_event_condition(SortedSpecies,HybridDynamicReactants,1,PropTresh,MaxParticleChange,Stream),
print_static_alpha(Number_of_alpha,Reactions_data, 1, 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/'),write(Stream,Rate),write(Stream,','),nl(Stream),
write(Stream,' ran = random_float,'),nl(Stream),
%print_event_constraints(Constraints_list, 1, Stream),
%concat_species(Stoch_species_order, 1, '', Species_for_event),
%write(Stream, Species_for_event),
%write(Stream, '], [rand, Time + delta_t, rand'),
%%% print('Checkpoint 2\n'),
%print_updating_alpha(Number_of_alpha, 1, Stream),
%%% print('Checkpoint 3\n'),
print_hybrid_species_constraint(Hybrid_reactants_list, 1, Stream),
%%% print('Checkpoint 4\n'),
% 2. Build alpha list and alpha subsums.
% 2. Build alpha list and alpha subsums.
length(Reactions_data, Number_of_alpha),
% Build Alpha_list
......@@ -1519,7 +1447,6 @@ print_dynamic_event(SortedSpecies,Reactions_data, HybridDynamicReactants,MaxPart
nth1(1,Alpha_list, First_alpha),
append([], [First_alpha], Init_list),
alpha_accumulation(Alpha_list, Init_list, Alpha_accumulation_list, 2),
%format("Alpha_accumulation_list is ~w ~n",[Alpha_accumulation_list]),
% 3. Print reactions in event
......@@ -1527,10 +1454,7 @@ print_dynamic_event(SortedSpecies,Reactions_data, HybridDynamicReactants,MaxPart
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),
format("Checkpoint ~n",[]),
%format("Print_helper_list is ~w ~n",[Print_helper_list]),
%%% print('Checkpoint 5\n'),
% 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).
......@@ -932,7 +932,7 @@ apply_changes(State,Time, Parameter = Expression) :-
%! euler(+Equations, +Time, +CurrentState, +CurrentTable)
euler(Equations,Time,CurrentState,CurrentTable) :-
format("One event occured ~n",[]),
% format("One event occured ~n",[]),
retractall(euler_last_state(_,_)),
assertz(euler_last_state(Time,CurrentState)),
nb_getval(euler_step_size,EulerStepSize),
......
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