Merge branch 'feature/Hybrid' into develop

arithmetic_rules.pl

@@ -9,6 +9,9 @@
     normalize_number/2
   ]).
 
+% Insert here for separate compilation and linting
+:- use_module(doc).
+:- use_module(util).
 
 :- dynamic(canonical/3).
 

events.pl

@@ -3,6 +3,7 @@
   [
     % Public API
     add_event/2,
+    add_time_event/2,
     list_events/0,
     list_model_events/0
   ]
@@ -51,10 +52,22 @@ add_event(Condition, ParameterValues) :-
   add_item([kind: event, item: event(Condition, ParameterValues)]).
 
 
+add_time_event(Condition, ParameterValues) :-
+  biocham_command(*),
+  type(Condition, condition),
+  Condition =.. [_Functor, 'Time', Threshold],
+  type(ParameterValues, '*'(parameter_name = arithmetic_expression)),
+  doc('Introduce a special kind of event when Time is the only variable, they may be used
+    for efficiency reason during numerical integration. The only allowed syntax is of the
+    form: Time > a or Time >= b.'),
+  add_item([kind: time_event, item: event('Time'>Threshold, ParameterValues)]).
+
+
 list_events :-
   biocham_command,
   doc('lists all the declared events.'),
-  list_items([kind: event]).
+  list_items([kind: event]),
+  list_items([kind: time_event]).
 
 
 :- devdoc('\\section{Private predicates}').
@@ -65,6 +78,7 @@ list_model_events :-
     lists all the events in a loadable syntax
     (auxiliary predicate of list_model).
   '),
+  % List the regular events
   \+ (
     item(
       [no_inheritance, kind: event, item: event(Condition, ParameterValues)]
@@ -78,4 +92,19 @@ list_model_events :-
       ),
       write(').\n')
     )
+  ),
+  % List the time events
+  \+ (
+    item(
+      [no_inheritance, kind: time_event, item: event(Condition, ParameterValues)]
+    ),
+    \+ (
+      format('add_time_event(~w, ', [Condition]),
+      write_successes(
+        member(Parameter = Value, ParameterValues),
+        write(', '),
+        format('~w = ~w', [Parameter, Value])
+      ),
+      write(').\n')
+    )
   ).

hybrid.plt

+:- use_module(hybrid).
+
+example_reactions([[[1*gen],[1*tem],c1*'gen_stoch'],[[1*tem],[],c2*('tem_total'/(100*1))],[[1*tem],[1*tem,1*gen],c3*('tem_total'/(100*1))],[[1*gen,1*struc],[1*virus],c4*('gen_stoch'*('struc_total'/(100*1)))]]).
+
+:- begin_tests(hybrid).
+
+test(build_constraints_list) :-
+  hybrid:build_constraints_list([[a],[0],[b,c]], [lower_level1, '0', lower_level3]).
+
+test(build_reaction_data) :-
+	Sorted_species = [gen,struc,tem,virus],
+  example_reactions(List),
+	hybrid:build_reaction_data(List, Sorted_species, Data),
+  Data = [[['-1','+0','+1','+0'],c1*'gen_stoch'],[['+0','+0','-1','+0'],c2*('tem_total'/(100*1))],[['+1','+0','+0','+0'],c3*('tem_total'/(100*1))],[['-1','-1','+0','+1'],c4*('gen_stoch'*('struc_total'/(100*1)))]].
+
+test(compute_maximum_change) :-
+  example_reactions(List),
+  hybrid:compute_maximum_change(List,2).
+
+:- end_tests(hybrid).

library/examples/hybrid/Goutsias_dynamic.bc

+MA(c1) for a => a + b.
+MA(c2) for b => _.
+MA(c3) for c => a + c.
+MA(c4) for a => _.
+MA(c5) for d + e => c.
+MA(c6) for c => d + e.
+MA(c7) for c + e => f.
+MA(c8) for f => c + e.
+MA(c9) for b + b => e.
+MA(c10) for e => b + b.
+
+parameter(c1 = 0.043).
+parameter(c2 = 7e-4).
+parameter(c3 = 71.5).
+parameter(c4 = 3.9e-6).
+parameter(c5 = 0.02).
+parameter(c6 = 0.48).
+parameter(c7 = 4e-4).
+parameter(c8 = 9e-12).
+parameter(c9 = 0.08).
+parameter(c10 = 0.5).
+
+present(a,0).
+present(b,10).
+present(c,0).
+present(d,2).
+present(e,10).
+present(f,0).

library/examples/hybrid/Goutsias_input_ode.bc

+MA(c9) for b + b => e.
+MA(c10) for e => b + b.
+
+parameter(c9 = 0.08).
+parameter(c10 = 0.5).
+
+present(b,10).
+present(e,10).

library/examples/hybrid/Goutsias_input_stoch.bc

+MA(c1) for a => a + b.
+MA(c2) for b => _.
+MA(c3) for c => a + c.
+MA(c4) for a => _.
+MA(c5) for d + e => c.
+MA(c6) for c => d + e.
+MA(c7) for c + e => f.
+MA(c8) for f => c + e.
+
+parameter(c1 = 0.043).
+parameter(c2 = 7e-4).
+parameter(c3 = 71.5).
+parameter(c4 = 3.9e-6).
+parameter(c5 = 0.02).
+parameter(c6 = 0.48).
+parameter(c7 = 4e-4).
+parameter(c8 = 9e-12).
+
+present(a,0).
+present(b,10).
+present(c,0).
+present(d,2).
+present(e,10).
+present(f,0).

library/examples/hybrid/T7.bc

+MA(c1) for gen => tem.
+MA(c2) for tem => _.
+MA(c3) for tem => tem + gen.
+MA(c4) for gen + struc => virus.
+MA(c5) for tem => tem + struc.
+MA(c6) for struc => _.
+
+parameter(c1 = 0.025).
+parameter(c2 = 0.25).
+parameter(c3 = 1).
+parameter(c4 = 0.0000075).
+parameter(c5 = 1000).
+parameter(c6 = 1.99).
+
+present(gen,0).
+present(tem).
+present(struc,0).
+present(virus,0).
\ No newline at end of file

library/examples/hybrid/T7_hybrid.bc

+parameter(gen_stoch = 0).
+parameter(virus_stoch = 0).
+parameter(alpha1 = 0).
+parameter(alpha2 = 0).
+parameter(alpha3 = 0).
+parameter(alpha4 = 0).
+parameter(alpha_sum = 0).
+
+% Hybrid species
+present(tem).
+parameter(tem_stoch = 0).
+function(tem_total =tem+tem_stoch).
+parameter(tem_total =tem+tem_stoch).
+
+function(virus = virus_stoch).
+
+parameter(struc_stoch = 0).
+function(struc_total = struc+struc_stoch).
+parameter(struc_total = struc+struc_stoch).
+
+tem_total * 1000 for tem => tem + struc.
+struc_total * 1.99 for struc => _.
+
+% stochastic
+parameter(tau = 0).
+parameter(ran = 0).
+add_event(Time > tau,
+ alpha1 = gen_stoch * 0.025,
+ alpha2 = tem_total * 0.25,
+ alpha3 = tem_total,
+ alpha4 = gen_stoch * struc_total * 0.0000075,
+ alpha_sum =  alpha1 + alpha2 + alpha3 + alpha4,
+ tau = (if alpha_sum <= 0 then inf else Time+(-1 / alpha_sum) * log(random_float)),
+ ran = random_float,
+ gen_stoch = (if alpha_sum*ran>0 and alpha_sum*ran<=alpha1 then gen_stoch-1
+ else if alpha_sum*ran>alpha1+alpha2 and alpha_sum*ran<=alpha1+alpha2+alpha3 then gen_stoch+1
+ else if alpha_sum*ran>alpha1+alpha2+alpha3 and alpha_sum*ran<=alpha1+alpha2+alpha3+alpha4 then gen_stoch-1
+ else gen_stoch),
+ tem_stoch = (if alpha_sum*ran>0 and alpha_sum*ran<=alpha1 then tem_stoch+1
+ else if alpha_sum*ran>alpha1 and alpha_sum*ran<=alpha1+alpha2 then tem_stoch-1
+ else tem_stoch),
+ struc_stoch = (if alpha_sum*ran>alpha1+alpha2+alpha3 and alpha_sum*ran<=alpha1+alpha2+alpha3+alpha4 then struc_stoch-1
+ else struc_stoch),
+ virus_stoch = (if alpha_sum*ran>alpha1+alpha2+alpha3 and alpha_sum*ran<=alpha1+alpha2+alpha3+alpha4 then virus_stoch+1
+ else virus_stoch)).
+
+option(show: {struc_total,virus}).
\ No newline at end of file

library/examples/hybrid/T7_input_ode.bc

+MA(c5) for tem => tem + struc.
+MA(c6) for struc => _.
+
+parameter(c5 = 1000). 
+parameter(c6 = 1.99). 
+
+present(tem, 1). 
+present(struc, 0). 
\ No newline at end of file

library/examples/hybrid/T7_input_stoch.bc

+MA(c1) for gen => tem.
+MA(c2) for tem => _.
+MA(c3) for tem => tem + gen.
+MA(c4) for gen + struc => virus.
+
+parameter(c1 = 0.025).
+parameter(c2 = 0.25).
+parameter(c3 = 1).
+parameter(c4 = 0.0000075).
+
+present(gen, 0).
+present(tem, 1).
+present(struc, 0).
+present(virus, 0).
\ No newline at end of file

library/examples/hybrid/Test_dynamic.bc

+MA(k1) for A + 2*B => C.
+MA(k2) for C => 2*A.
+
+parameter(k1 = 1).
+parameter(k2 = 1).
+
+present(A,1).
+present(B,1).
+present(C,0).
\ No newline at end of file

library/examples/hybrid/Test_dynamic_out.bc

+parameter(tau = 0).
+parameter(ran = 0).
+parameter(delta_t = 0).
+
+parameter(k1_diff = 0).
+parameter(k2_diff = 0).
+
+function(kk1 = k1_diff).
+function(kk2 = k2_diff).
+
+
+function(A_total = floor(100*1*A + A_stoch)).
+parameter(A_stoch = 0).
+function(AA = A_stoch).
+function(B_total = floor(100*1*B + B_stoch)).
+parameter(B_stoch = 0).
+function(BB = B_stoch).
+function(C_total = floor(100*1*C + C_stoch)).
+parameter(C_stoch = 0).
+function(CC = C_stoch).
+
+k1_diff*(k1*((A_total/(100*1))*(B_total/(100*1))^2))for A+2*B=>C.
+k2_diff*(k2*(C_total/(100*1)))for C=>2*A.
+present(A,1).
+present(B,1).
+present(C,0).
+parameter(
+  k1 = 1,
+  k2 = 1
+).
+parameter(alpha1 = 0).
+parameter(alpha2 = 0).
+parameter(alpha_sum = 0).
+parameter(lower_level1 = 0).
+parameter(lower_level2 = 0).
+parameter(lower_level3 = 0).
+add_event(Time > tau,
+ k1_diff = (if k1*(A_total*B_total^2) >= 1/(step_size * 20) and A_total >= 10 and B_total >= 10 then 1 else 0),
+ k2_diff = (if k2*C_total >= 1/(step_size * 20) and C_total >= 20 then 1 else 0),
+ alpha1 = (1-k1_diff)*(k1/2*A_total/100*B_total/100*(B_total+ -1)/100),
+ alpha2 = (1-k2_diff)*(k2*C_total/100),
+ alpha_sum = alpha1 + alpha2,
+ delta_t = (if alpha_sum <= 0 then step_size else (-1 / alpha_sum) * log(random_float)),
+ tau = Time+delta_t/100,
+ ran = random_float,
+ lower_level1 = (if A_total >= 1 then 1 else 0),
+ lower_level2 = (if B_total >= 2 then 1 else 0),
+ lower_level3 = (if C_total >= 1 then 1 else 0),
+ A_stoch = (if alpha_sum*ran>0 and alpha_sum*ran<=alpha1 and lower_level1>0 then A_stoch-1
+ else if alpha_sum*ran>alpha1 and alpha_sum*ran<=alpha1 + alpha2 and lower_level2>0 then A_stoch+2
+ else A_stoch),
+ B_stoch = (if alpha_sum*ran>0 and alpha_sum*ran<=alpha1 and lower_level1>0 then B_stoch-2
+ else B_stoch),
+ C_stoch = (if alpha_sum*ran>0 and alpha_sum*ran<=alpha1 and lower_level1>0 then C_stoch+1
+ else if alpha_sum*ran>alpha1 and alpha_sum*ran<=alpha1 + alpha2 and lower_level2>0 then C_stoch-1
+ else C_stoch)
+).
+
+% Plot part
+
+option(show: {A_total, B_total, C_total}). % Choose here the species to plot
+numerical_simulation(method:rsbk, time: 20).
+plot.
+option(show: {}).
\ No newline at end of file

library/examples/hybrid/biomodel9_input_ode.bc

+MA(r2a__a2) for P_KKK+E2=>E2_P_KKK.
+MA(r2a__d2) for E2_P_KKK=>P_KKK+E2.
+MA(r2b__k2) for E2_P_KKK=>E2+KKK.
+MA(r3a__a3) for KK+P_KKK=>P_KKK_KK.
+MA(r3a__d3) for P_KKK_KK=>KK+P_KKK.
+MA(r3b__k3) for P_KKK_KK=>P_KK+P_KKK.
+MA(r4a__a4) for P_KK+KKPase=>KKPase_P_KK.
+MA(r4a__d4) for KKPase_P_KK=>P_KK+KKPase.
+MA(r4b__k4) for KKPase_P_KK=>KK+KKPase.
+MA(r5a__a5) for P_KK+P_KKK=>P_KKK_P_KK.
+MA(r5a__d5) for P_KKK_P_KK=>P_KK+P_KKK.
+MA(r5b__k5) for P_KKK_P_KK=>PP_KK+P_KKK.
+MA(r6a__a6) for PP_KK+KKPase=>KKPase_PP_KK.
+MA(r6a__d6) for KKPase_PP_KK=>PP_KK+KKPase.
+MA(r6b__k6) for KKPase_PP_KK=>P_KK+KKPase.
+MA(r7a__a7) for K+PP_KK=>PP_KK_K.
+MA(r7a__d7) for PP_KK_K=>K+PP_KK.
+MA(r7b__k7) for PP_KK_K=>P_K+PP_KK.
+MA(r8a__a8) for P_K+KPase=>KPase_P_K.
+MA(r8a__d8) for KPase_P_K=>P_K+KPase.
+MA(r8b__k8) for KPase_P_K=>K+KPase.
+MA(r9a__a9) for P_K+PP_KK=>PP_KK_P_K.
+MA(r9a__d9) for PP_KK_P_K=>P_K+PP_KK.
+MA(r9b__k9) for PP_KK_P_K=>PP_KK+PP_K.
+MA(r10a__a10) for PP_K+KPase=>KPase_PP_K.
+MA(r10a__d10) for KPase_PP_K=>PP_K+KPase.
+MA(r10b__k10) for KPase_PP_K=>P_K+KPase.
+present(E2,0.0003).
+present(KKK,0.003).
+present(P_KKK,0.0).
+present(KK,1.2).
+present(P_KK,0.0).
+present(PP_KK,0.0).
+present(K,1.2).
+present(P_K,0.0).
+present(PP_K,0.0).
+present(KPase,0.12).
+present(KKPase,0.0003).
+
+present(E2_P_KKK,0.0).
+present(P_KKK_KK,0.0).
+present(P_KKK_P_KK,0.0).
+present(PP_KK_K,0.0).
+present(PP_KK_P_K,0.0).
+present(KKPase_PP_KK,0.0).
+present(KKPase_P_KK,0.0).
+present(KPase_PP_K,0.0).
+present(KPase_P_K,0.0).
+parameter(
+  compartment = 4.0e-12,
+  K_PP_norm_max = 0.900049,
+
+  r2a__a2 = 1000.0,
+  r2a__d2 = 150.0,
+  r2b__k2 = 150.0,
+  r3a__a3 = 1000.0,
+  r3a__d3 = 150.0,
+  r3b__k3 = 150.0,
+  r4a__a4 = 1000.0,
+  r4a__d4 = 150.0,
+  r4b__k4 = 150.0,
+  r5a__a5 = 1000.0,
+  r5a__d5 = 150.0,
+  r5b__k5 = 150.0,
+  r6a__a6 = 1000.0,
+  r6a__d6 = 150.0,
+  r6b__k6 = 150.0,
+  r7a__a7 = 1000.0,
+  r7a__d7 = 150.0,
+  r7b__k7 = 150.0,
+  r8a__a8 = 1000.0,
+  r8a__d8 = 150.0,
+  r8b__k8 = 150.0,
+  r9a__a9 = 1000.0,
+  r9a__d9 = 150.0,
+  r9b__k9 = 150.0,
+  r10a__a10 = 1000.0,
+  r10a__d10 = 150.0,
+  r10b__k10 = 150.0
+  ).
+
+function(
+  K_PP_norm = (PP_K+KPase_PP_K)/(PP_K+P_K+K+PP_KK_K+PP_KK_P_K+KPase_PP_K+KPase_P_K),
+  rel_K_PP_max = K_PP_norm/K_PP_norm_max,
+  KK_PP_norm = (PP_KK+PP_KK_K+PP_KK_P_K+KKPase_PP_KK)/(PP_KK+P_KK+KK+PP_KK_K+PP_KK_P_K+P_KKK_KK+P_KKK_P_KK+KKPase_PP_KK+KKPase_P_KK),
+  KKK_P_norm = (P_KKK+P_KKK_KK+P_KKK_P_KK)/(KKK+P_KKK+P_KKK_KK+P_KKK_P_KK)
+).
\ No newline at end of file

library/examples/hybrid/biomodel9_input_stoch.bc

+MA(r1a__a1) for KKK+E1=>E1_KKK.
+MA(r1a__d1) for E1_KKK=>KKK+E1.
+MA(r1b__k2) for E1_KKK=>E1+P_KKK.
+
+present(E1,3.0e-5).
+present(KKK,0.003).
+present(E1_KKK,0.0).
+present(P_KKK,0.0).
+
+parameter(  
+  r1a__a1 = 1000.0,
+  r1a__d1 = 150.0,
+  r1b__k2 = 150.0).
\ No newline at end of file

library/examples/hybrid/cell_input_ode.bc

+
+
+
+
+
+MA(k3)   for a => b.
+
+
+MA(k4)   for b + 2* c => 3*c.
+
+
+MA(k8)   for h => d.
+MA(k9)   for d => h.
+
+present(h = 1).
+
+%function(YT = g + f + [b] + [c]).
+% function(CT = b + c + d + h).
+
+
+parameter(k3  = 200).   % k3*CT = 200
+
+parameter(k4  = 180).   % 10 < k4/CT^2 < 1000
+
+parameter(k8  = 10000). % k8 >> k9
+parameter(k9  = 100).   % k9 >> k6
\ No newline at end of file

library/examples/hybrid/cell_input_stoch.bc

+MA(k1)   for _=>g.
+MA(k2)   for g=>_.
+MA(k4p)  for b => c.
+MA(k5)   for c => b.
+
+MA(k6)   for c => e.
+MA(k7)   for f =>_.
+
+
+parameter(k1  = 0.1). % k1*CT = 0.015
+parameter(k2  = 0).
+parameter(k5  = 0).
+parameter(k6  = 1).     % 0.1 < k6 < 10
+parameter(k7  = 0.6).
+parameter(k4p = 0.018).

library/initial.bc

@@ -7,6 +7,7 @@ option(
   error_epsilon_relative: 1e-6,
   initial_step_size: 1e-6,
   maximum_step_size: 1e-2,
+  minimum_step_size: 1e-5,
   precision: 6, %5
   nusmv_initial_states: all,
   nusmv_counter_example: no,

markov.pl

@@ -26,8 +26,12 @@ stochastic_simulation(Until, Propensities, Time, Boolean) :-
   retractall(last_time(_)),
   assertz(last_time(0.0)),
   get_option(stochastic_thresholding, Threshold),
+  %statistics(runtime,_),
   markov_run(InitState, Conditions, Actions, Until, Propensities, Time,
     Threshold),
+  %statistics(runtime,[_,T]),
+  %TT is T/1000,
+  %format("Simulation time: ~ps~n",[TT]),
   make_state_table,
   restore_parameters.
 

matrix.pl

+:- module(matrix,
+  [
+    constant_list/3,
+    constant_matrix/4,
+    replace/4,
+    get_element/3,
+    replace_element/3,
+    init_matrix/1,
+    write_matrix/1,
+    ludcmp_matrix/2,
+    ludcmp_matrix/0,
+    matrix_sum/4,
+    matrix_prod/3,
+    scalar_prod/3,
+    matrix_apply/3,
+    matrix_sum/3,
+    create_diag_matrix/3,
+    sub_matrix/1,
+    transpose_matrix/2,
+    lusubst/2,
+    solve_linear_problem/3
+  ] 
+).
+
+:- use_module(util).
+
+/* Module Matrix
+*
+* Several operations are done on a "global variable" of SWI-prolog (using nb_set/getval)
+* but this is due to historical version of an efficient implementation in GNU-prolog and
+* may not be necessary nor useful now.
+*/
+
+
+%! constant_list(+Length, +Element, -List)
+%
+% create a list with every element set up to Element
+
+constant_list(0,_,[]) :-
+  !. 
+
+constant_list(N1,X,[X|L]) :- 
+  N1 > 0, N is N1 - 1,    
+  constant_list(N,X,L).
+
+
+%! constant_matrix(+N_row, +N_column, +Element, -List)
+%
+% create a matrix with every element set up to Element
+
+constant_matrix(RowSize,ColumnSize,Value,Matrix) :-
+  constant_list(ColumnSize,Value,Row), % The row size is the number of column
+  constant_list(RowSize,Row,Matrix). 
+
+
+%! replace(+List, +Index, +NewElement, -NewList)
+%
+% Allows to replace the i-th element of the list, where the 1st element is index by 1
+
+replace([_|Tail], 1, NewElement, [NewElement|Tail]) :- 
+  !.
+   
+replace([Element|Tail], I, NewElement, [Element|Tail2]):-
+  I > 1, I1 is I-1,
+  replace(Tail, I1, NewElement, Tail2).   
+
+
+%! init_matrix(+Matrix)
+%
+% Creates the global variable which will stores the matrix on which we are doing operations 
+
+init_matrix([Head|Tail]):-
+  length([Head|Tail],RowSize),
+  length(Head,ColumnSize),
+  nb_setval(stored_matrix,[Head|Tail]),
+  nb_setval(formal_m,RowSize),
+  nb_setval(formal_n,ColumnSize).
+
+
+%! get_element(+RowIndex, +ColumnIndex, -Value)
+%
+% get the (i,j)-th element of the stored_matrix
+
+get_element(RowIndex,ColumnIndex,Value) :- 
+  nb_getval(stored_matrix,Matrix),
+  nth1(RowIndex,Matrix,Row),
+  nth1(ColumnIndex,Row,Value).
+
+
+%! replace_element(+RowIndex, +ColumnIndex, +Value) :- 
+%
+% set the (i,j)-th element of the stored matrix to Value
+
+replace_element(RowIndex,ColumnIndex,Value) :- 
+  nb_getval(stored_matrix,Matrix),
+  nth1(RowIndex,Matrix,Row),
+  replace(Row,ColumnIndex,Value,NewRow),
+  replace(Matrix,RowIndex,NewRow,NewMatrix),
+  nb_setval(stored_matrix,NewMatrix).
+
+
+/* Some tools to write matrix in a pretty way, for debugging for examples*/
+%! write_matrix(+Matrix)
+
+write_matrix([]).
+
+write_matrix([H|T]):-
+  write_tab(H),
+  write_matrix(T).
+
+
+%! write_tab(+List)
+
+write_tab([]) :- 
+  nl.
+
+write_tab([H|T]):-
+  write(H),write('\t'),
+  write_tab(T).
+
+%! ludcmp_matrix(+Matrix, -TransformedMatrix)
+%
+% Perform an LU decomposition, that is find two matrices L and U such that Matrix = L.U,
+% L and U begin respectively lower and upper diagonal.
+% By convention, L has only 1 on its diagonal so that we gather them in a single matrix M
+% such that M[i,j] = L[i,j] if i<j and M[i,j] = U[i,j] if i>=j.  M is the transformed
+% matrix and is stored in the global variable at the end of the call.
+