Commit 1aecead3 authored by Sylvain Soliman's avatar Sylvain Soliman

Merge branch 'release/4.4.10'

parents f9a5eba5 a801e491
# make jupyter will make biocham without checking unitary tests (avoiding potential problems of access to BioModels in sbml.plt)
ADDITIONAL_MODULES= \
modules/sbml/sbml_utils.pl \
modules/partialfrac/partialfraction.pl
......
......@@ -8,7 +8,7 @@
about/0
]).
version('4.4.9').
version('4.4.10').
copyright(
'Copyright (C) 2003-2020 Inria, EPI Lifeware, Saclay-Île de France, France'
......
FROM registry.gitlab.inria.fr/lifeware/biocham:v4.4.9
FROM registry.gitlab.inria.fr/lifeware/biocham:v4.4.10
{
"name": "gui",
"version": "4.4.9",
"version": "4.4.10",
"description": "biocham gui in jupyter notebook",
"main": "src/index.js",
"scripts": {
......
......@@ -188,6 +188,7 @@ commands = [
"select_model",
"select_ode_system",
"select_table",
"sensitivity",
"set_attribute",
"set_dimension",
"set_graph_name",
......@@ -195,6 +196,7 @@ commands = [
"set_model_name",
"set_ode_system_name",
"set_p_m_rate",
"test_rate_independence",
"transition",
"tropicalize",
"undefined",
......
"""Example magic"""
__version__ = '4.4.9'
__version__ = '4.4.10'
......@@ -198,8 +198,10 @@ def handle_bokeh_plot(fig, gui_fig, against, cds=None, gui_cds=None, old_handle=
gui_layout = gui_fig
if against != 'Time':
button = Button(label="Animate", callback=animation_callback(cds))
gui_button = Button(label="Animate", callback=animation_callback(gui_cds))
button = Button(label="Animate")
button.js_on_click(animation_callback(cds))
gui_button = Button(label="Animate")
gui_button.js_on_click(animation_callback(gui_cds))
layout = column(fig, button)
gui_layout = column(gui_fig, gui_button, sizing_mode="scale_width")
......
......@@ -32,8 +32,9 @@ compile_cpp_program(SourceFiles, Options, ExecutableFilename) :-
;
CXX = 'g++'
),
setenv('OMPI_MPICXX', CXX),
compile_program(
CXX, SourceFiles,
'mpicxx', SourceFiles,
%FF no warning of deprecated PPL code
% ['--std=c++11', '-Dtypeof(x)=decltype(x)' | Options],
['-w', '--std=c++11', '-Dtypeof(x)=decltype(x)' | Options],
......@@ -41,7 +42,7 @@ compile_cpp_program(SourceFiles, Options, ExecutableFilename) :-
).
%% TODO
%% TODO
%% - robust parsing (w.r.t. spaces for instance)
%% - error handling (unknown library)
pkg_compile_options(PkgName, Options) :-
......
......@@ -132,7 +132,7 @@ ctl(E) :-
:- doc('\\emphright{steady(f) is a shorthand for EG(f)}').
:- doc('\\emphright{stable(f) is a shorthand for AG(f)}').
:- doc('\\emphright{checkpoint(f, g) is a shorthand for not EU(not f, g)). Note that such a factual property does not imply any causality relationship from f to g.}').
:- doc('\\emphright{checkpoint2(f, g) is a shorthand for not(g)/\\\\ EF(g)/\\\\checkpoint(f,g)}').
:- doc('\\emphright{checkpoint2(f, g) is a shorthand for not(g)/\\\\ EF(g)/\\\\checkpoint(f,g)}').
:- doc('\\emphright{oscil2(f) is a shorthand for EU(not(f), f /\\\\ EU(f, not(f) /\\\\ EU(not(f), f))) which is true if f has at least two peaks on one path}').
:- doc('\\emphright{oscil3(f) is a shorthand for EU(not(f), f /\\\\ EU(f, not(f) /\\\\ EU(not(f), f /\\\\ EU(f, not(f) /\\\\ EU(not(f), f))))) which is true if f has at least three peaks on one path}').
:- doc('\\emphright{oscil(f) is a shorthand for oscil3(f) /\\\\ EG(EF(f) /\\\\ EF(not f)) which is a necessary (not sufficient) condition for infinite oscillations in CTL}').
......
......@@ -289,7 +289,9 @@ reduction_methods(no).
reduction_methods(native).
reduction_methods(sat).
reduction_methods(sat_species).
reduction_methods(sat_reactions).
%! compile_from_pivp(+PIVP_string, +Output)
......@@ -381,16 +383,20 @@ compile_from_pivp(PIVP, Input, Output) :-
:- doc('
\\begin{example} Compilation of the Hill function of order 2 as a function of an input $h(x)=x^2/(1+x^2)$ \n
\\begin{example} Compilation of the Hill function of order 2 as a function of an input $h(x)=x^2/(1+x^2)$ .
This time \\texttt{sat_species} is used (minimizing the number of variables as is the case for \\texttt{native}) and then
\\texttt{sat_reactions} that minimizes the number of monomes.
').
%:- biocham(option(lazy_negatives:yes)).
:- biocham(option(binomial_reduction:native)).
:- biocham(option(binomial_reduction:sat_species)).
:- biocham(compile_from_pivp((0.0,d(h)/dt= 2*n^2*t;1.0,d(n)/dt= -2*n^2*t;0.0,d(t)/dt=1),x,h)).
:- biocham(parameter(input=2)).
:- biocham(list_model).
:- biocham(numerical_simulation(time:10)). % method:msbdf
:- biocham(plot(show:{h})).
:- biocham(plot(show:{h}, logscale:x)).
:- biocham(compile_from_pivp((0.0,d(h)/dt= 2*n^2*t;1.0,d(n)/dt= -2*n^2*t;0.0,d(t)/dt=1),x,h, binomial_reduction: sat_reactions)).
:- biocham(list_model).
:- doc('
\\end{example}
').
......@@ -440,6 +446,7 @@ main_compiler(PIVP_input, Name_list_raw, Input):-
PIVP_unsigned = PIVP_non_bin,
Name_list_unsigned = Name_list_non_bin
),
debug(pivp, "PIVP ~w reduced to ~w", [PIVP_non_bin, PIVP_unsigned]),
get_option(lazy_negatives, Lazyness),
(
Lazyness = yes
......@@ -448,6 +455,7 @@ main_compiler(PIVP_input, Name_list_raw, Input):-
;
lng:rewrite_PIVP_all_negated(PIVP_unsigned, PIVP, VarNeg)
),
debug(pivp, "PIVP lazy to ~w", [PIVP]),
sort(VarNeg, VarNegSort),
negate_name(Name_list_unsigned, VarNegSort, Species_names),
(
......@@ -457,7 +465,10 @@ main_compiler(PIVP_input, Name_list_raw, Input):-
;
true
),
compile_to_biocham_model(PIVP, Species_names).
debug(pivp, "Negated names ~w", [Species_names]),
compile_to_biocham_model(PIVP, Species_names),
% FIXME several results above???
!.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......@@ -532,21 +543,25 @@ parse_pivp((_,(d(_)/dt=P)), L, [Pol]):-
% parse_monomial/3.
parse_polynomial(+ X, L, P):-
parse_polynomial(X, L, P).
!,
parse_polynomial(X, L, P).
parse_polynomial(- X, L, [[D, R]]):-
parse_monomial(X, L, [C, R]),
D is -C.
!,
parse_monomial(X, L, [C, R]),
D is -C.
parse_polynomial(X + Y, L, R):-
parse_polynomial(X, L, R1),
parse_polynomial(Y, L, R2),
append(R1, R2, R).
!,
parse_polynomial(X, L, R1),
parse_polynomial(Y, L, R2),
append(R1, R2, R).
parse_polynomial(X - Y, L, [[D, M] | R1]):-
parse_polynomial(X, L, R1),
parse_monomial(Y, L, [C, M]),
D is -C.
!,
parse_polynomial(X, L, R1),
parse_monomial(Y, L, [C, M]),
D is -C.
parse_polynomial(X, L, [M]):-
parse_monomial(X, L, M).
......@@ -1300,7 +1315,7 @@ breadth_search(PODE, ListVariable, Redvar) :-
displace_exponent([1],0,Lm,First_Var),
get_option(binomial_reduction, Reduction),
(
Reduction == sat
atom_concat('sat', _, Reduction)
->
maxsat_binomial_reduction(ListVariable, ListDerivative, First_Var, Redvar)
;
......@@ -1808,51 +1823,112 @@ stutter(N, In, Out) :-
% Calls a MAX-Sat solver to find the smallest binomial reduction of PODE
maxsat_binomial_reduction(AllVariables, AllDerivatives, OutVar, Solution) :-
debug(pivp, "all vars: ~w~nall derivs: ~w", [AllVariables, AllDerivatives]),
nb_setval(count_clauses, 0),
tmp_file_stream(text, FileRoot, StreamFileInInit),
close(StreamFileInInit),
% we create a boolean var for each possible ODE var
% except the output one they will all be associated to a soft clause
% hence their number is a good Top value for the hard clauses
length(AllVariables, N),
nth1(OutIndex, AllVariables, OutVar),
derivatives_to_variables(AllDerivatives, AllVariables, Coverings),
% we overestimate the number of clauses, but rc2.py is ok with that
maplist(length, Coverings, NCov),
sum_list(NCov, NClauses),
tmp_file_stream(text, FileRoot, StreamFileInInit),
close(StreamFileInInit),
string_concat(FileRoot, ".wcnf", FileIn),
debug(pivp, "Will be writing clauses to: ~w", [FileIn]),
string_concat(FileRoot, ".wcnf", FileIn),
get_option(binomial_reduction, Reduction),
(
Reduction == sat_species
->
% we have N-1 soft clauses for each variable
Top = N,
NVar = N,
NMono = 0,
M = 0
;
derivatives_to_allmonomials(AllDerivatives, AllMonomials),
debug(pivp, "all mono: ~w", [AllMonomials]),
append(AllMonomials, MonoClauses),
length(MonoClauses, NMono),
sort(MonoClauses, SMonoClauses),
length(SMonoClauses, M),
% we have M soft clauses for each unique monomial
Top is M + 1,
NVar is N + M
),
setup_call_cleanup(
open(FileIn, write, TmpStream),
(
% fake header that is overwrite safe
format(TmpStream, "fake hdr cccccccccccccccccc~n", []),
derivatives_to_variables(AllDerivatives, AllVariables, TmpStream, Top)
),
close(TmpStream)
),
% we overestimate the number of clauses, but rc2.py is ok with that
nb_getval(count_clauses, NCov),
NClauses is NMono + Top + NCov,
setup_call_cleanup(
open(FileIn, update, UpStream),
% header
format(UpStream, "p wcnf ~d ~d ~d~n", [NVar, NClauses, Top]),
close(UpStream)
),
setup_call_cleanup(
open(FileIn, write, Stream),
open(FileIn, append, Stream),
(
% header
format(Stream, "p wcnf ~d ~d ~d~n", [N, NClauses, N]),
write_variable_clauses(Stream, N, OutIndex),
write_derivative_covering_clauses(Stream, N, Coverings)
write_variable_clauses(Stream, N, M, OutIndex, Top, Reduction),
(
Reduction == sat_reactions
->
write_monomial_clauses(Stream, N, Top, SMonoClauses, AllMonomials)
;
true
)
),
close(Stream)
),
debug(pivp, "~d clauses written to: ~w", [NClauses, FileIn]),
debug(pivp, "Now calling MaxSAT", []),
call_cleanup(
(
% there is always a solution, so we get solution and weight
run_sat(FileIn, yes, [ResVariables, _Weight]),
variable_string_to_solution(ResVariables, AllVariables, Solution)
variable_string_to_solution(ResVariables, AllVariables, Solution),
debug(pivp, "Solution: ~w", [Solution])
),
delete_file(FileIn)
(
have_to_delete_temporary_files
->
delete_file(FileIn)
;
true
)
).
%! write_variable_clauses(+Stream, +N, +OutIndex) is det.
%! write_variable_clauses(+Stream, +N, +M, +OutIndex, +Top, +Reduction) is det.
%
% Write to stream one clause per possible variable to include.
% it is a hard clause (weight is N) if it is the output variable
% it is a hard clause (weight is Top) if it is the output variable
% otherwise the negation is a soft clause
write_variable_clauses(Stream, N, OutIndex) :-
write_variable_clauses(Stream, N, M, OutIndex, Top, Reduction) :-
forall(
between(1, N, I),
(
Reduction == sat_species
->
between(1, N, I)
;
(
% if we minimize monomials, no soft clauses for variables
I = OutIndex
;
Start is N + 1,
End is N + M,
between(Start, End, I)
)
),
(
(
I == OutIndex
->
Weight = N,
Weight = Top,
Value = I
;
Weight = 1,
......@@ -1863,40 +1939,35 @@ write_variable_clauses(Stream, N, OutIndex) :-
).
%! write_derivative_covering_clauses(Stream, N, Coverings)
%! write_monomial_clauses(Stream, N, SMonoClauses, AllMonomials) is det.
%
% Write clause that impose that if a variable is present, its derivative is covered by
% present variables
write_derivative_covering_clauses(Stream, N, Coverings) :-
% Write clauses imposing that if a variable is present, all monomials of its derivative are present too
write_monomial_clauses(Stream, N, Top, SMonoClauses, AllMonomials) :-
forall(
between(1, N, I),
(
nth1(I, Coverings, Cov),
II is -I,
nth1(I, AllMonomials, Monomials),
maplist(get_mono(SMonoClauses, N), Monomials, Numbers),
forall(
(
member(Co, Cov),
% Optimisation on trivially satisfied clauses
\+ member(I, Co)
),
(
% II is not at the right position in clause, but rc2.py is ok
% with that too
atomic_list_concat([II | Co], ' ', Clause),
format(Stream, "~d ~w 0~n", [N, Clause])
)
member(Number, Numbers),
format(Stream, "~d -~w ~w 0~n", [Top, I, Number])
)
)
).
get_mono(List, Offset, Element, Position) :-
nth1(Index, List, Element),
Position is Index + Offset.
%! variable_string_to_solution(+String: string, AllVariables: list, -Solution: list) is det.
%
% Convert the output of the MaxSAT solver (-1 2 3 -4 ...) to a list of the
% variables that are present (i.e. positive)
variable_string_to_solution(String, AllVariables, Solution) :-
split_string(String, " ", " ", VarStrings),
maplist(to_int, VarStrings, VarNumbers),
maplist(atom_to_int, VarStrings, VarNumbers),
% keep only the present variables
exclude(>(0), VarNumbers, SolutionIndexes),
findall(
......@@ -1909,32 +1980,52 @@ variable_string_to_solution(String, AllVariables, Solution) :-
).
%! derivatives_to_variables(+AllDerivatives, +AllVariables, -Coverings) is det.
%! derivatives_to_variables(+AllDerivatives, +AllVariables, +TmpStream, +Top) is det.
%
% Return a list corresponding to the clauses for covering each derivative
derivatives_to_variables(AllDerivatives, AllVariables, Coverings) :-
maplist(derivative_to_variables(AllVariables), AllDerivatives, Coverings).
% Write to TmpStream the clauses for covering each derivative
derivatives_to_variables(AllDerivatives, AllVariables, TmpStream, Top) :-
nb_setval(varcount, 1),
maplist(derivative_to_variables(AllVariables, TmpStream, Top), AllDerivatives).
%! derivative_to_variables(+AllVariables, +Derivative, -Covering) is det.
%! derivative_to_variables(+AllVariables, +TmpStream, -Top, +Derivative) is det.
%
% Covering is a list of clauses covering all monomials of Derivative
derivative_to_variables(AllVariables, Derivative, Covering) :-
maplist(monomial_to_variables(AllVariables), Derivative, Coverings),
% we flatten our conjunction of conjunctions
append(Coverings, Covering).
% Write to TmpStream clauses covering all monomials of Derivative
derivative_to_variables(AllVariables, TmpStream, Top, Derivative) :-
maplist(monomial_remove_coeff, Derivative, Monomials),
maplist(monomial_to_variables(AllVariables, TmpStream, Top), Monomials),
nb_getval(varcount, Count),
debug(pivp, "var ~d done", [Count]),
CCount is Count + 1,
nb_setval(varcount, CCount).
%! monomial_to_variables(+AllVariables, +Monomial, -Covering) is det.
%! monomial_remove_coeff([+Coeff, +Monomial], -Monomial) is det.
%
% Remove the coefficient of a monomial
monomial_remove_coeff([_, M], M).
derivatives_to_allmonomials(AllDerivatives, AllMonomials) :-
maplist(derivative_to_allmonomials, AllDerivatives, Monomials),
append(Monomials, SomeMonomials),
sort(SomeMonomials, AllMonomials).
derivative_to_allmonomials(Derivative, Monomials) :-
maplist(monomial_remove_coeff, Derivative, Monomials).
%! monomial_to_variables(+AllVariables, +TmpStream, +Top, +Monomial) is det.
%
% make a list of clauses for the possible coverings of Monomial (ignoring its constant
% part) by variables of AllVariables
monomial_to_variables(_Variables, [_, Monomial], []) :-
monomial_to_variables(_Variables, _Stream, _Top, Monomial) :-
% The derivative is constant, we return the empty conjunction
max_list(Monomial, 0),
!.
monomial_to_variables(AllVariables, [_, Monomial], Covering) :-
monomial_to_variables(AllVariables, Stream, Top, Monomial) :-
% all sums of two different variables
findall(
[I1, I2],
......@@ -1963,27 +2054,27 @@ monomial_to_variables(AllVariables, [_, Monomial], Covering) :-
;
DNFCovering = SQCoverings
),
dnf_cnf(DNFCovering, Covering),
debug(pivp, "Monomial: ~w~nDNF: ~w~nCNF: ~w~n~n", [Monomial, DNFCovering, Covering]).
nb_getval(varcount, Count),
debug(pivp, "Var: ~d~nMonomial: ~w~nDNF: ~w~n~n", [Count, Monomial, DNFCovering]),
dnf_cnf(DNFCovering, Count, Top, Stream).
%! dnf_cnf(+DNF:list, -CNF:list) is det.
%! dnf_cnf(+DNF:list, +Count:integer, +Top:integer, +Stream) is det.
%
% simple DNF to CNF conversion
dnf_cnf(DNF, CNF) :-
findall(
L,
% writing to Stream the clause for variable Count
dnf_cnf(DNF, Count, Top, Stream) :-
forall(
(
maplist(member, L0, DNF),
sort(L0, L)
\+ member(Count, L0)
),
CNF
(
sort(L0, L),
atomic_list_concat([Count | L], ' ', Atom),
format(Stream, "~w -~a 0~n", [Top, Atom]),
nb_getval(count_clauses, C),
CC is C+1,
nb_setval(count_clauses, CC)
)
).
%! to_int(+String:atom, -Int:integer) is det.
%
% Read an Int from a String
to_int(String, Int) :-
read_term_from_atom(String, Int, []),
integer(Int).
......@@ -26,7 +26,7 @@ then
}
brew tap brewsci/science
brewlings='swi-prolog gsl gnuplot graphviz pkg-config ppl libsbml nusmv node'
brewlings='swi-prolog gsl gnuplot graphviz pkg-config ppl libsbml nusmv node open-mpi'
for brewling in $brewlings ; do
install_or_upgrade "$brewling"
done
......@@ -39,7 +39,7 @@ then
### Ubuntu or Debian
sudo apt-get -qy update
sudo apt-get -qy install lsb-release curl gcc g++ make pkg-config libgsl0-dev gnuplot libppl-dev graphviz-dev git nodejs npm
sudo apt-get -qy install lsb-release curl gcc g++ make pkg-config libgsl0-dev gnuplot libppl-dev graphviz-dev git nodejs npm openmpi-bin libopenmpi-dev
release=$(lsb_release -sc)
......
......@@ -38,7 +38,9 @@ a10*Kpp*KPase for Kpp+KPase=>KPase_Kpp.
b10*KPase_Kpp for KPase_Kpp=>Kpp+KPase.
c10*KPase_Kpp for KPase_Kpp=>Kp+KPase.
present(E1,3.0e-5).
present(E1,e1).
parameter(e1=3.0e-5).
present(KKK,0.003).
present(E2,0.0003).
......
This diff is collapsed.
% Cell cycle adapted from Qu et al, Biophysical journal
% kwpcn=5, Kweem=0.2 instead of 2 and 1 because of Cry-deficient mice
% ksmpf=0.35 (instead of 0.3) for shorter cell cycle in autonomous cell cycle model.
% Cell cycle molecules
present(preMPF,0.279939269).
present(MPF,0.747404087).
present(C25,0.629659946).
present(C25P,0.370340054).
present(Wee1,0.00223979517).
present(Wee1p,0.132059534).
present(APC,0.0303504407).
present(IE,0.41856651).
present(mWee1,0.141190132).
present(APCi,0.973282491).
% Cell cycle parameters
parameter(kd25=1).
parameter(kdmpfp=0.25).
%%% kampf values
% 5 gives free period ~20h
% 2.2 ~24h
% 1.35 ~28h -> period doubling
parameter(kampf=2.2).
parameter(kaapc=1).
parameter(ksweemp=0.001).
parameter(kdweem=1).
parameter(Kweem=0.2).
parameter(kwpcn=5).
parameter(Vac=1).
parameter(Kac25=0.01).
parameter(Kic25=0.01).
parameter(kdmpf=2).
parameter(Vaw=1).
parameter(Vapc=0.1).
parameter(Kaw=0.005).
parameter(Kiw=0.005).
parameter(Vapw=0.1).
parameter(Kapc=0.01).
parameter(ks25=1).
%parameter(kswee=0.5).
parameter(kswee=1.5).
%parameter(kimpf=1.5).
parameter(kimpf=2.5).
parameter(Vic=0.2).
%%% period when kapmf is 1.35
% 1.42 -> 19.0h
% 1.4 -> 21.6h
% 1.2 -> 24.9h
% 0.5 -> 28.0h
% 0.3 -> 30.5h
parameter(kiapc=0.5).
parameter(kdwee=1).
parameter(Viw=0.2).
parameter(ksmpf=0.35).
parameter(kdie=0.1).
parameter(ksiep=0).
parameter(ksie=0.3).
parameter(kaapcp=0).
% Reaction rules
% MPF - preMPF
ksmpf for _=>preMPF.
kdmpfp*preMPF for preMPF=>_.
kdmpfp*MPF for MPF=>_.
kdmpf*APC*MPF for MPF=[APC]=>_.
kdmpf*APC*preMPF for preMPF=[APC]=>_.
kampf*C25P*preMPF for preMPF=[C25P]=>MPF.
kimpf*Wee1*MPF for MPF=[Wee1]=>preMPF.
% Wee1
Viw*Wee1p/(Kiw+Wee1p) for Wee1p=>Wee1.
%(Vapw+Vaw*MPF)*Wee1/(Kaw+Wee1) for Wee1=[MPF]=>Wee1p.
Vapw*Wee1/(Kaw+Wee1) for Wee1=[MPF]=>Wee1p.
Vaw*MPF*Wee1/(Kaw+Wee1) for Wee1=[MPF]=>Wee1p.
ksweemp for _=> mWee1.
kdweem*mWee1 for mWee1=>_.
kswee*mWee1 for _=[mWee1]=>Wee1.
kdwee*Wee1 for Wee1=>_.
kdwee*Wee1p for Wee1p=>_.
% Cdc25
Vic*C25P/(Kic25+C25P) for C25P=>C25.
%(Vapc+Vac*MPF)*C25/(Kac25+C25) for C25=[MPF]=>C25P.
Vapc*C25/(Kac25+C25) for C25=[MPF]=>C25P.
Vac*MPF*C25/(Kac25+C25) for C25=[MPF]=>C25P.