Commit 338380c4 authored by Mathieu Hemery's avatar Mathieu Hemery Committed by SOLIMAN Sylvain
Browse files

Merge branch 'develop' into feature/xpp

parents d7336399 f1ad8e9e
......@@ -199,7 +199,7 @@ gadagne_restart: gadagne
ssh -t gadagne 'cd /opt/tmpnb/ && sudo ./cleanup.sh && sudo ./tmpnb_biocham.sh'
release:
echo "current version: " && grep '^version' about.pl
@echo "current version: " && grep '^version' about.pl
@read -p "Version number (e.g. 4.1.0): " version && \
git flow release start "$$version" && \
sed -i"" -e "s/^version('.*')/version('$$version')/" about.pl && \
......
......@@ -8,7 +8,7 @@
about/0
]).
version('4.4.13').
version('4.4.14').
copyright(
'Copyright (C) 2003-2020 Inria, EPI Lifeware, Saclay-Île de France, France'
......
......@@ -389,7 +389,7 @@ insert_coef(A, Coef, Result) :-
Coef > 0,
Coef < 1,
CoefInv is 1 / Coef,
0.0 is float_fractional_part(CoefInv)
0.0 =:= float_fractional_part(CoefInv)
->
normalize_number(CoefInv, CoefInvNorm),
Result = A / CoefInvNorm
......
FROM registry.gitlab.inria.fr/lifeware/biocham:v4.4.13
FROM registry.gitlab.inria.fr/lifeware/biocham:v4.4.14
......@@ -37,6 +37,10 @@ reduction_methods(no).
reduction_methods(native).
reduction_methods(fast).
reduction_methods(fastnSAT).
reduction_methods(sat_species).
reduction_methods(sat_reactions).
......@@ -123,7 +127,6 @@ list_reduction([N,PODE,IV], Name_list, [NewN,NewPODE,NewIV], Modified_Names) :-
Order > 2
->
generate_sufficient_variables(PODE, ListAllVariable),
length(ListAllVariable, Len), writeln(Len),
breadth_search(PODE, ListAllVariable, ListVariable),
length(ListVariable,NewN),
% binomial_branch_bound(ListAllVariable, PODE, ListVariable, NewN),
......@@ -230,6 +233,8 @@ highest([],[],[]).
% strictly_smaller(+ExpM,+Exp)
%
% Exp is smaller than ExpM
strictly_smaller(ExpM,Exp) :-
smaller(ExpM,Exp),
......@@ -551,6 +556,21 @@ breadth_search(PODE, ListVariable, Redvar) :-
atom_concat('sat', _, Reduction)
->
maxsat_binomial_reduction(ListVariable, ListDerivative, First_Var, Redvar)
;
member(Reduction, [fast, fastnSAT])
->
fast_search(PODE, First_Var, ListVariable, Redvar1),
(
Reduction = fast
->
Redvar = Redvar1
;
construct_all_derivatives(PODE, Redvar1, ListDerivative1),
with_option(
binomial_reduction:sat_species,
maxsat_binomial_reduction(Redvar1, ListDerivative1, First_Var, Redvar)
)
)
;
Test_list = [[First_Var]],
breadth_search_subroutine(Test_list, ListVariable, ListDerivative, _Next_level, Redvar_t0),
......@@ -746,7 +766,7 @@ convert_to_ode(_Id, [], [], _VL) :- !.
convert_poly_to_ode([Monomial], VarList, Expr) :-
!,
convert_mono_to_ode(Monomial, VarList, Expr).
convert_poly_to_ode([Monomial|Tail], VarList, Expr+Remainder) :-
convert_poly_to_ode([Monomial|Tail], VarList, Remainder+Expr) :-
!,
convert_mono_to_ode(Monomial, VarList, Expr),
convert_poly_to_ode(Tail, VarList, Remainder).
......@@ -1061,3 +1081,138 @@ is_molecule(A) :-
is_numeric(N) :-
item([kind:parameter, key:N]);
number(N).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% FAST BINOMIAL REDUCTION %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fast_search(+PODE, +First_var, +ListAllVariable, -BinomialSubSet)
%
% perform a fast but inoptimal reduction
fast_search(PODE, V0, AllVariables, BinomialSubSet) :-
delete(AllVariables, V0, TailVariables),
% construct_all_derivatives(PODE, TailVariables, ListDerivative),
add_variable_to_exponents([], PODE, V0, [], Exponents),
fast_search_sr([V0], BinomialSubSet, Exponents, [V0|TailVariables], PODE).
fast_search_sr(Acc, RAcc, [], _LAV, _LD) :-
reverse(Acc,RAcc).
fast_search_sr(OldVar, Acc, Exponents, LAV, PODE) :-
heuristic_few_expo(Exponents, OldVar, LAV, NewVar, PODE),
add_variable_to_exponents(Exponents, PODE, NewVar, OldVar, NewExponents),
fast_search_sr([NewVar|OldVar], Acc, NewExponents, LAV, PODE).
%! add_variable_to_exponents(+ExpOld, +PODE, +NewVar, +OldVar, -ExpNew)
%
% compute the new set of unreachable exponents (ExpNew) once we add NewVar to the set of
% variables
add_variable_to_exponents(ExpOld, PODE, NewVar, OldVar, ExpNew) :-
compute_derivative(NewVar, PODE, Derivative),
concatenate_exponents(ExpOld, Derivative, ExpTempo),
sieve_exponents(ExpTempo, [NewVar|OldVar], ExpNew).
%! heuristic_few_expo(+Exponents, +OldVar, +AllVar, -NewVar, +PODE)
%
% determine the NewVar to be added to the set of variable by looking first at all the
% variable that are useful alone (either because they are in the set exponent, or because
% added to an old variable they form a useful exponent) if this set is empty, it try any
% remainding variable. Among this set, it choose the one that will introduce as few new
% exponents as possible.
heuristic_few_expo(Exponents, OldVar, AllVar, NewVar, PODE) :-
maybe_useful(AllVar, Exponents, OldVar, Useful, PODE),
(
Useful = []
->
maybe_useless(AllVar, Exponents, OldVar, Useless, PODE),
keysort(Useless, [_K-NewVar|_UselessSortedTail])
;
keysort(Useful, [_K-NewVar|_UsefulSortedTail])
).
%! maybe_useful(+AllVar, +Exponents, +OldVar, -ChoiceList, +PODE)
%
% Determine the list and cost of all variable that are useful alone
% Choicelist is a list of the form [Cost-Variable|Tail]
maybe_useful([], _Exponents, _OldVar, [], _PODE).
maybe_useful([Head|TailAll], Exponents, OldVar, [Key-Head|TailCut], PODE) :-
(
member(Head, Exponents)
;
member(Other, OldVar),
add_list(Other, Head, Expo),
member(Expo, Exponents)
),!,
estimate_value(Exponents, PODE, Head, OldVar, Key),
maybe_useful(TailAll, Exponents, OldVar, TailCut, PODE).
maybe_useful([_Head|TailAll], Exponents, OldVar, TailCut, PODE) :-
maybe_useful(TailAll, Exponents, OldVar, TailCut, PODE).
%! maybe_useless(+AllVar, +Exponents, +OldVar, -ChoiceList, +PODE)
%
% Determine the list and cost of all variable that are not already used
% Choicelist is a list of the form [Cost-Variable|Tail]
maybe_useless([], _Exponents, _OldVar, [], _PODE).
maybe_useless([Head|TailAll], Exponents, OldVar, [Key-Head|TailCut], PODE) :-
\+(member(Head, OldVar)),
member(Expo, Exponents),
strictly_smaller(Expo, Head),
!,
estimate_value(Exponents, PODE, Head, OldVar, Key),
maybe_useless(TailAll, Exponents, OldVar, TailCut, PODE).
maybe_useless([_Head|TailAll], Exponents, OldVar, TailCut, PODE) :-
maybe_useless(TailAll, Exponents, OldVar, TailCut, PODE).
%! estimate_value(+Exponents, +PODE, +NewVar, +OldVar, -Cost)
%
% Determine the cost (number of new exponents to be computed) of NewVar
estimate_value(ExpOld, PODE, NewVar, OldVar, Value) :-
add_variable_to_exponents(ExpOld, PODE, NewVar, OldVar, ExpNew),
length(ExpNew, Value).
%! concatenate_exponents(+Exponents, +Derivative, -Exponent2)
%
% add all the exponent of the derivative to the list of Exponents
concatenate_exponents(Exponents, Derivative, Exponent2):-
extract_exponent(Derivative, TempoDer),
append(Exponents, TempoDer, Exponent2).
extract_exponent([], []) :- !.
extract_exponent([[_R,E]|TailD], [E|TailE]) :- extract_exponent(TailD, TailE).
%! sieve_exponents(+Exponents, -LV, +NewExponents)
%
% remove from the list of exponents all the ones that are binomially reachable with the
% set of variables LV
sieve_exponents([], _Var_list, []) :- !.
sieve_exponents([Exp|Tail1], Var_list, Tail2) :-
(
constant(Exp)
;
member(Exp,Var_list)
;
member(X, Var_list), member(Y, Var_list),
add_list(X, Y, Exp)
),
!,
sieve_exponents(Tail1, Var_list, Tail2).
sieve_exponents([Exp|Tail1], Var_list, [Exp|Tail2]) :-
sieve_exponents(Tail1, Var_list, Tail2).
......@@ -6,11 +6,11 @@ list_reduction:- use_module(library(plunit)).
test(polynomial_ODE_True) :-
command(a => b),
command(3*a^2 for a+a => c),
reduction:polynomial_ODE.
with_output_to(atom(_), reduction:polynomial_ODE).
test(polynomial_ODE_False) :-
command(a/b for a => b),
\+( reduction:polynomial_ODE ).
with_output_to(atom(_), \+ reduction:polynomial_ODE).
test(generate_names) :-
reduction:generate_names([a,'Out'], [[1,2],[1,0]], [aOut2, a]).
......@@ -84,6 +84,6 @@ test(clean_writing) :-
test(convert_poly_to_ode) :-
reduction:convert_poly_to_ode([[1, [1,1]]], [a,b], a*b),
reduction:convert_poly_to_ode([[1, [1,1]], [-2, [0,2]]], [a,b], a*b + -(2*b^2)).
reduction:convert_poly_to_ode([[1, [1,1]], [-2, [0,2]]], [a,b], -(2*b^2) + a*b).
:- end_tests(binomial_reduction).
{
"name": "gui",
"version": "4.4.13",
"version": "4.4.14",
"description": "biocham gui in jupyter notebook",
"main": "src/index.js",
"scripts": {
......
......@@ -147,6 +147,7 @@ commands = [
"list_stable_states",
"list_tables",
"list_tscc_candidates",
"list_units",
"load",
"load_biocham",
"load_ginml",
......@@ -188,7 +189,9 @@ commands = [
"rate_independence_reduction_inputs_sinks",
"reaction_graph",
"reaction_model",
"read_xpp",
"reduce_model",
"remove_fraction",
"rename_column",
"rename_table",
"reset_options",
......@@ -213,6 +216,7 @@ commands = [
"set_model_name",
"set_ode_system_name",
"set_p_m_rate",
"set_units",
"test_rate_independence",
"test_rate_independence_inputs_sinks",
"test_rate_independence_invariants",
......@@ -223,4 +227,5 @@ commands = [
"variation",
"which_p_m_mode",
"with_timer",
"write_xpp",
]
"""Example magic"""
__version__ = '4.4.13'
__version__ = '4.4.14'
......@@ -69,6 +69,7 @@ def get_data_and_meta(filename):
to_show = metadata[1].split(' ')[1:]
against = metadata[2].split(' ')[1]
axes = metadata[3].split(' ')[1:5]
axes_names = (' '.join(metadata[4].split(' ')[1:]), ' '.join(metadata[5].split(' ')[1:]))
with open(filename, newline='') as csvfile:
lines = list(csv.reader(csvfile))
......@@ -78,7 +79,8 @@ def get_data_and_meta(filename):
'logscale': logscale,
'to_show': to_show,
'against': against,
'axes': axes}
'axes': axes,
'axes_names': axes_names}
def create_bokeh_plot(data, gui=False):
......@@ -100,7 +102,8 @@ def create_bokeh_plot(data, gui=False):
bokeh_loaded = True
fig = figure(
x_axis_label=against,
x_axis_label=data['axes_names'][0],
y_axis_label=data['axes_names'][1],
x_axis_type='log' if 'x' in logscale else 'linear',
y_axis_type='log' if 'y' in logscale else 'linear',
tools='pan,wheel_zoom,box_zoom,reset,save,undo',
......
......@@ -35,7 +35,7 @@ The following files can also be used to export some Biocham objects:
:- doc('\\section{Biocham call options}').
do_arguments :-
current_prolog_flag(argv, [_Exec | Args]),
current_prolog_flag(os_argv, [_Exec | Args]),
do_arguments(Args).
......
......@@ -673,9 +673,9 @@ extract_names((_IC, d(Name)/dt = _Deriv; Remainder), [Name|Tail]) :-
% A monomial is a unique object of the form: [K,[Exp1,Exp2,Exp3]]
% corresponding to the term K v1^Exp1 V2^Exp2 etc.
:- dynamic(input_name/1,
output_name/1,
name_spec/0).
:- dynamic(input_name/1).
:- dynamic(output_name/1).
:- dynamic(name_spec/0).
%! compare_exponants(+ExponentSet1,+ExponentSet2,-Comparator)
%
......
......@@ -75,6 +75,8 @@ find_invar_aux(ForcedMax, Type, OtherType, Operator) :-
Goal2 =.. [Type, P],
Goal2,
\+(member(P, Vars)),
% can be added by normalized_path
\+(base_mol([P])),
assertz(base_mol([P]))
),
[P]
......
......@@ -136,7 +136,7 @@ load(InputFile) :-
type(InputFile, input_file),
doc('
acts as the corresponding \\command{load_biocham/1} / \\command{load_sbml/1}
/ \\command{load_ode/1} / \\command{load_table/1},
/ \\command{import_ode/1} / \\command{load_table/1},
depending on the file extension
(respectively \\texttt{.bc}, \\texttt{.xml}, \\texttt{.ode}, \\texttt{.csv}
-- assuming no extension is \\texttt{.bc}).'),
......
......@@ -86,7 +86,7 @@ normalized([], []) :- !.
normalized(A, Anorm) :-
leading(A, Lc),
mult(float(1 rdiv Lc), A, Anorm).
mult(float(1 / Lc), A, Anorm).
% subtracts two polynomials
subtract(A, B, C) :-
......@@ -123,7 +123,7 @@ long_division(A, B, Q, R) :-
leading(A, LcA),
leading(B, LcB),
D is DegA - DegB,
C is float(LcA rdiv LcB),
C is float(LcA / LcB),
monomial_mult(D, C, B, B1),
subtract(A, B1, A1),
long_division(A1, B, Q1, R), !,
......@@ -163,7 +163,7 @@ simplify_fraction(A, B, Aprime, Bprime) :-
normalized_euclid(A, B, Unorm, Vnorm, Dnorm) :-
euclid(A, B, U, V, D),
leading(D, LcD),
C is float(1 rdiv LcD),
C is float(1 / LcD),
mult(C, U, Unorm),
mult(C, V, Vnorm),
mult(C, D, Dnorm), !.
......@@ -271,7 +271,7 @@ partial_fraction(P, Q, A, QL, ML) :-
% the following three lines ensure Q is unitary
% without altering the fraction P/Q
leading(Q, Lead),
PCoef is float(1 rdiv Lead),
PCoef is float(1 / Lead),
mult(PCoef, P, P_normal),
roots(Q, QLraw, ML),
......
......@@ -133,6 +133,7 @@ filter(only_extrema).
:- initial(option(filter: no_filter)).
:- initial(option(stats: no)).
:- dynamic(last_method/1).
:- devdoc('\\section{Commands}').
......@@ -179,6 +180,9 @@ numerical_simulation :-
'filtering function for the trace'
),
doc('performs a numerical simulation from time 0 up to a given time.'),
retractall(last_method(_)),
debug(numsim, "numerical simulation method: ~w", [Method]),
assertz(last_method(Method)),
statistics(runtime,_),
(
Method = ssa
......@@ -203,7 +207,7 @@ numerical_simulation :-
;
Method = rsbk
->
assertz(rosenbrock_running),
assertz(rosenbrock_running),
with_current_ode_system(
with_clean(
[
......@@ -218,7 +222,6 @@ numerical_simulation :-
solve2
)
),
assertz(last_method_rsbk),
retractall(rosenbrock_running),
retractall(rosenbrock_running2)
;
......@@ -233,8 +236,7 @@ numerical_simulation :-
],
solve
)
),
retractall(last_method_rsbk)
)
),
statistics(runtime,[_,T]),
(
......@@ -258,7 +260,8 @@ continue :-
get_initial_values(Variables, InitialValues),
set_initial_values(Variables, Values),
% FIXME uses global options, not the ones used for previous simulation
( last_method_rsbk
(
last_method(rsbk)
->
nb_setval(rsbk_continue,yes),
with_option([time: Time,method: rsbk], numerical_simulation),
......@@ -305,6 +308,7 @@ prepare_numerical_simulation_options(Options) :-
enumerate_variables,
convert_ode,
gather_equations(Equations),
debug(numsim, "ODEs: ~w", [Equations]),
(
rosenbrock_running
->
......@@ -319,12 +323,16 @@ prepare_numerical_simulation_options(Options) :-
true
),
gather_initial_values(InitialValues),
debug(numsim, "Initial values: ~w", [InitialValues]),
gather_events(RegularEvents, event),
gather_events(TimeEvents, time_event),
append(RegularEvents, TimeEvents, Events),
debug(numsim, "Events: ~w", [Events]),
gather_fields(Events, Fields),
gather_initial_parameter_values(InitialParameterValues),
debug(numsim, "Parameter values: ~w", [InitialParameterValues]),
jacobian(Equations, Jacobian),
debug(numsim, "Jacobian: ~w", [Jacobian]),
Options = [
fields: Fields,
equations: Equations,
......@@ -371,7 +379,9 @@ solve2 :-
solve :-
prepare_numerical_simulation_options(Options),
debug(numsim, "simulation options: ~w", [Options]),
solve(Options, Table),
debug(numsim, "Done", []),
add_table('numerical_simulation', Table).
......@@ -408,6 +418,7 @@ enumerate_nonconstant_parameters(Events, Header: Parameter) :-
),
NCParameters
),
debug(numsim, "NC parameters ~w", [NCParameters]),
member(ParameterIndex, NCParameters),
Parameter = p(ParameterIndex),
convert_identifier(Header, Parameter).
......@@ -417,7 +428,9 @@ enumerate_nonparametric_functions(Header: expression(Expr)) :-
% (item([kind: function, item: function(Header = Body)]) ;
get_ode_function(Header, Body),
atomic(Header),
debug(numsim, "NP function ~w with body ~w", [Header, Body]),
ode:substitute_functions(Body, NoFuncBody),
debug(numsim, "substituted body ~w", [NoFuncBody]),
convert_expression(NoFuncBody, Expr).
......@@ -473,11 +486,13 @@ gather_initial_values(InitialValues) :-
gather_initial_parameter_values(InitialParameterValues) :-
peek_count(parameter_counter, ParameterCount),
ParameterMax is ParameterCount - 1,
debug(numsim, "~d parameters", [ParameterCount]),
findall(
InitialParameterValue,
(
between(0, ParameterMax, ParameterIndex),
parameter_index(Parameter, ParameterIndex),
debug(numsim, "Parameter #~d: ~w", [ParameterIndex, Parameter]),
(%parameter_value(Parameter, InitialParameterValue)
ode_parameter_value(Parameter, InitialParameterValue)
; conditional_parameter_value(ParameterIndex,InitialParameterValue) % virtual event parameters for conditional kinetics
......
......@@ -22,6 +22,7 @@
add_reactions_from_ode_system/0,
add_influences_from_ode_system/0,
(init)/1,
remove_fraction/0,
% load_ode_system/1,
load_reactions_from_ode/1,
% add_ode_system/1,
......@@ -52,7 +53,9 @@
ode/2,
ode/3,
ode_add_expression_to_molecule/2,
ode_predicate/1
ode_predicate/1,
remove_fraction/2,
normalize_ode/1
]).
% Only for separate compilation/linking
......@@ -125,7 +128,9 @@ list_ode :-
(
get_current_ode_system(Id),
list_items([parent: Id, kind: ode]),
list_items([parent: Id, kind: initial_concentration])
list_items([parent: Id, kind: initial_concentration]),
list_items([parent: Id, kind: parameter]),
list_items([parent: Id, kind: function])
)
)
;
......@@ -415,7 +420,31 @@ load_influences_from_ode_system :-
add_influences_from_ode_system.
remove_fraction :-
biocham_command,
doc('Remove the rational fraction in the current ODE system by multiplying the derivative
by the GCD. WARNING: the resulting ODE system is NOT equivalent to the starting one.'),
new_ode_system(NewId),
set_ode_system_name(NewId, without_fraction_ode),
normalize_ode(ListOde),
maplist(remove_fraction, ListOde, NewListOde),
maplist(add_ode(NewId), NewListOde),
with_current_ode_system((
get_current_ode_system(Id),
forall(
item([parent: Id, kind:parameter, item: Item]),
add_item([parent: NewId, kind:parameter, item: Item])
),
forall(
item([parent: Id, kind: function, item: Item]),
add_item([parent: NewId, kind: function, item: Item])
),
forall(
item([parent: Id, kind:initial_concentration, item: Item]),
add_item([parent: NewId, kind:initial_concentration, item: Item])
)
)),
select_ode_system(without_fraction_ode).
:- devdoc('\\section{Public API}').
......@@ -560,6 +589,18 @@ add_ode(Id, ODE) :-
add_item([parent: Id, kind: ode, key: X, item: Item]).
change_ode(Id, ODE) :-
(
parse_ode(ODE, X, Item)
->
true
;
throw(error(illformed_ode(ODE)))
),
delete_item([parent: Id, kind: ode, key: X]),
add_item([parent: Id, kind: ode, key: X, item: Item]).
import_reactions_from_ode_system(Id) :-
check_cleaned(ode:assoc/2),
enumerate_terms_in_odes(Id),
......@@ -663,25 +704,10 @@ non_decomposable_term_in_expression(Coefficient, Term, E) :-
substitute_functions(E, ENoFunc) :-
(
get_current_ode_system(Id)
->
all_items([parent: Id, kind: function], Functions),
maplist(split_func, Functions, Definitions, Bodies),
maplist(
rewrite(substitute(Definitions, Bodies)),
Bodies,
NoFuncBods
),
substitute(Definitions, NoFuncBods, E, ENoFunc)
;
% case of add_function, there is no current system, but no functions
% either
ENoFunc = E
).