Merge branch 'feature/units' into develop

biocham_jupyter/guinbextension/src/config/workflows/commands.py

@@ -148,6 +148,7 @@ commands = [
     "list_stable_states",
     "list_tables",
     "list_tscc_candidates",
+    "list_units",
     "load",
     "load_biocham",
     "load_ginml",
@@ -214,6 +215,7 @@ commands = [
     "set_model_name",
     "set_ode_system_name",
     "set_p_m_rate",
+    "set_units",
     "test_rate_independence",
     "test_rate_independence_invariants",
     "transition",

biocham_jupyter/kernel/biocham_kernel/images.py

@@ -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',

modules/sbml/sbml.pl

@@ -12,6 +12,9 @@
     model_getListOfParameters/2,
     model_getListOfRules/2,
     model_getSpeciesById/3,
+    model_getTimeUnits/2,
+    model_getSubstanceUnits/2,
+    model_getVolumeUnits/2,
     listOf_get/3,
     listOf_size/2,
     reaction_getKineticLaw/2,

modules/sbml/sbml_swiprolog.c

@@ -567,6 +567,63 @@ pl_compartment_getVolume(term_t compartment_term, term_t volume_term) {
 }
 
 
+static foreign_t
+pl_model_getTimeUnits(
+    term_t model_term, term_t units_term
+) {
+    const char *units;
+    Model_t *model;
+    UnitDefinition_t *unitdef;
+    PL_check(PL_get_model(model_term, &model));
+    if (Model_isSetTimeUnits(model)) {
+        PL_check(units = Model_getTimeUnits(model));
+    } else {
+        PL_check(unitdef = Model_getUnitDefinitionById(model, "time"));
+        PL_check(units = UnitDefinition_getName(unitdef))
+    }
+    PL_check(PL_unify_atom_chars(units_term, units));
+    PL_succeed;
+}
+
+
+static foreign_t
+pl_model_getSubstanceUnits(
+    term_t model_term, term_t units_term
+) {
+    const char *units;
+    Model_t *model;
+    UnitDefinition_t *unitdef;
+    PL_check(PL_get_model(model_term, &model));
+    if (Model_isSetSubstanceUnits(model)) {
+        PL_check(units = Model_getSubstanceUnits(model));
+    } else {
+        PL_check(unitdef = Model_getUnitDefinitionById(model, "substance"));
+        PL_check(units = UnitDefinition_getName(unitdef))
+    }
+    PL_check(PL_unify_atom_chars(units_term, units));
+    PL_succeed;
+}
+
+
+static foreign_t
+pl_model_getVolumeUnits(
+    term_t model_term, term_t units_term
+) {
+    const char *units;
+    Model_t *model;
+    UnitDefinition_t *unitdef;
+    PL_check(PL_get_model(model_term, &model));
+    if (Model_isSetVolumeUnits(model)) {
+        PL_check(units = Model_getVolumeUnits(model));
+    } else {
+        PL_check(unitdef = Model_getUnitDefinitionById(model, "volume"));
+        PL_check(units = UnitDefinition_getName(unitdef))
+    }
+    PL_check(PL_unify_atom_chars(units_term, units));
+    PL_succeed;
+}
+
+
 PL_extension sbml_predicates[] = {
     { "readSBML", 2, (pl_function_t) pl_readSBML, 0 },
     { "sbmlDocument_getNumErrors", 2, (pl_function_t) pl_sbmlDocument_getNumErrors, 0 },
@@ -580,6 +637,9 @@ PL_extension sbml_predicates[] = {
     { "model_getListOfRules", 2, (pl_function_t) pl_model_getListOfRules, 0 },
     { "model_getListOfFunctions", 2, (pl_function_t) pl_model_getListOfFunctions, 0 },
     { "model_getSpeciesById", 3, (pl_function_t) pl_model_getSpeciesById, 0 },
+    { "model_getTimeUnits", 2, (pl_function_t) pl_model_getTimeUnits, 0 },
+    { "model_getSubstanceUnits", 2, (pl_function_t) pl_model_getSubstanceUnits, 0 },
+    { "model_getVolumeUnits", 2, (pl_function_t) pl_model_getVolumeUnits, 0 },
     { "listOf_get", 3, (pl_function_t) pl_listOf_get, 0 },
     { "listOf_size", 2, (pl_function_t) pl_listOf_size, 0 },
     { "reaction_getKineticLaw", 2, (pl_function_t) pl_reaction_getKineticLaw, 0 },

numerical_simulation.pl

@@ -133,6 +133,7 @@ filter(only_extrema).
 :- initial(option(filter: no_filter)).
 :- initial(option(stats: no)).
 
+:- dynamic(last_method/1).
 
 :- devdoc('\\section{Commands}').
 
@@ -179,7 +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
@@ -204,7 +207,7 @@ numerical_simulation :-
   ;
     Method = rsbk
   ->
-  assertz(rosenbrock_running),
+    assertz(rosenbrock_running),
     with_current_ode_system(
       with_clean(
         [
@@ -219,7 +222,6 @@ numerical_simulation :-
         solve2
       )
     ),
-    assertz(last_method_rsbk),
     retractall(rosenbrock_running),
     retractall(rosenbrock_running2)
   ;
@@ -234,8 +236,7 @@ numerical_simulation :-
         ],
         solve
       )
-    ),
-    retractall(last_method_rsbk)
+    )
   ),
   statistics(runtime,[_,T]),
   (
@@ -259,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),

plot.pl

@@ -192,13 +192,16 @@ gnu_plot_csv :-
   get_option(ymax,Ymax),
   maplist(unquoted_term_to_atom, Show, ShowAtoms),
   atomic_list_concat(ShowAtoms, ' ', ShowAtom),
+  get_axes_names(Against, XAxis, YAxis),
   with_output_to_file(
     MetaDataFilename,
     (
       format('logscale: ~a~n', [Axes]),
       format('show: ~a~n', [ShowAtom]),
       format('against: ~a~n', [Against]),
-      format('axes: ~w ~w ~w ~w~n', [Xmin, Xmax, Ymin, Ymax])
+      format('axes: ~w ~w ~w ~w~n', [Xmin, Xmax, Ymin, Ymax]),
+      format('xaxis_name: ~w~n', [XAxis]),
+      format('yaxis_name: ~w~n', [YAxis])
     )
   ),
   view_image(Filename).
@@ -261,9 +264,11 @@ set datafile separator ","
     format(Stream, 'set logscale ~a~n', [Axes])
   ),
   against(Against),
+  get_axes_names(Against, XAxis, YAxis),
   % heisenbug if choicepoints before process_create
   once(nth1(IndexX, ['Time' | Headers], Against)),
-  format(Stream, 'set xlabel "~a"~n', [Against]),
+  format(Stream, 'set xlabel "~a"~n', [XAxis]),
+  format(Stream, 'set ylabel "~a"~n', [YAxis]),
   nb_setval(plotted, 0),
   forall(
     (

sbml_files.pl

@@ -73,6 +73,7 @@ add_sbml_document(SBML) :-
 
 
 add_sbml_model(Model) :-
+  add_units(Model),
   add_compartments(Model),
   add_global_parameters(Model),
   add_assignment_rules(Model),
@@ -215,6 +216,35 @@ add_reactions(Model) :-
     )
   ).
 
+
+add_units(Model) :-
+  debug(sbml, "getting units…", []),
+  (
+    model_getTimeUnits(Model, Time)
+  ->
+    debug(sbml, "getTimeUnits(~w)", [Time]),
+    set_units(time, Time)
+  ;
+    debug(sbml, "no timeUnits defined", [])
+  ),
+  (
+    model_getSubstanceUnits(Model, Substance)
+  ->
+    debug(sbml, "getSubstanceUnits(~w)", [Substance]),
+    set_units(substance, Substance)
+  ;
+    debug(sbml, "no substanceUnits defined", [])
+  ),
+  (
+    model_getVolumeUnits(Model, Volume)
+  ->
+    debug(sbml, "getVolumeUnits(~w)", [Volume]),
+    set_units(volume, Volume)
+  ;
+    debug(sbml, "no volumeUnits defined", [])
+  ).
+
+
 % Store SBML-definde compartments/locations
 :- dynamic(compartment/1).
 :- dynamic(single_compartment/1).

units.pl

 :- module(units,
   [
+    unit_type/1,
+    set_units/2,
+    list_units/0,
+    get_axes_names/3,
     list_dimensions/0,
     set_dimension/3,
     clear_dimensions/0,
@@ -10,9 +14,90 @@
 :- use_module(doc).
 :- use_module(library(clpfd)).
 
+:- grammar(unit_type).
 
-% store declared units
-:- dynamic(k_unit/3).
+unit_type(time).
+
+unit_type(volume).
+
+unit_type(substance).
+
+
+%! set_units(+Type, +Unit) is det.
+%
+% changes the default Type unit to Unit.
+set_units(Type, Unit) :-
+  biocham_command,
+  type(Type, unit_type),
+  type(Unit, name),
+  doc('Set the default units for the given type (time, volume or substance) in the current model.'),
+  change_item([], unit, Type, Unit).
+
+
+%! default_unit(+Type, +Unit) is det.
+%
+% defines the default units for each possible type
+default_unit(time, 's').
+
+default_unit(volume, 'l').
+
+default_unit(substance, 'mol').
+
+
+%! get_units(+Type, -Unit) is det.
+%
+% unifies Unit with the current default unit of type Type.
+get_units(Type, Unit) :-
+  (
+    item([kind: unit, key: Type, item: Unit])
+  ->
+    true
+  ;
+    default_unit(Type, Unit)
+  ).
+
+
+%! list_units is det.
+%
+% Lists the values of all current default units for all types.
+list_units :-
+  biocham_command,
+  doc('Lists the default units of the current model.'),
+  forall(
+    unit_type(Type),
+    (
+      get_units(Type, Unit),
+      format("~w: ~w~n", [Type, Unit])
+    )
+  ).
+
+
+%! get_axes(+Against, -XAxis, -YAxis) is det.
+%
+% unifies XAxis and YAxis with atoms corresponding to the correct axes names
+get_axes_names(Against, XAxis, YAxis) :-
+  units:get_units(substance, SubstUnit),
+  numerical_simulation:last_method(Method),
+  (
+    memberchk(Method, [ssa, spn, pn, sbn, bn])
+  ->
+    YAxis = SubstUnit
+  ;
+    units:get_units(volume, VolUnit),
+    format(atom(YAxis), "~w/~w", [SubstUnit, VolUnit])
+  ),
+  (
+    Against == 'Time'
+  ->
+    units:get_units(time, TimeUnit),
+    format(atom(XAxis), "Time (~w)", [TimeUnit])
+  ;
+    format(atom(XAxis), "~w (~w)", [Against, YAxis])
+  ).
+
+
+% store declared dimensions
+:- dynamic(k_dimension/3).
 
 :- biocham_silent(clear_model).
 
@@ -21,7 +106,7 @@
 %
 % biocham command
 %
-% Infers dimensions (time and volume, without units) of all parameters of the model.
+% Infers dimensions (time and volume, dimensionless) of all parameters of the model.
 %
 % Molecular concentration have dimension volume^-1
 % Reaction rates (and influence forces) have dimension time^-1
@@ -33,7 +118,7 @@
 
 :- doc('
    Dimensional analysis infers dimensions for model parameters and checks their consistency.
-   In Biocham, only time and volume dimensions are considered, without units.
+   In Biocham, only time and volume dimensions are considered.
    The dimension of a molecular concentration is volume$^{-1}$.
    The dimension of a kinetic expression (i.e. reaction rate or influence force) is time$^{-1}$.
    Warning: numbers have no dimension so kinetic parameters should be used in kinetic expressions instead of directly numbers.
@@ -75,28 +160,28 @@ set_dimension(P, U, V) :-
   type(P, parameter_name),
   type(U, number),
   type(V, number),
-  retractall(k_unit(P, _, _)),
-  assertz(k_unit(P, U, V)).
+  retractall(k_dimension(P, _, _)),
+  assertz(k_dimension(P, U, V)).
 
 
 clear_dimensions :-
   biocham_command,
   doc('Clear all declared dimensions'),
-  retractall(k_unit(_ , _, _)).
+  retractall(k_dimension(_ , _, _)).
 
 
 clear_dimension(P) :-
   biocham_command,
   doc('Clear declared dimension for parameter \\argument{P}'),
   type(P, parameter_name),
-  retractall(k_unit(P, _, _)).
+  retractall(k_dimension(P, _, _)).
 
 
 check_set_dim([]).
 
 check_set_dim([((U, V) - P) | L]) :-
   (
-    k_unit(P, U, V),
+    k_dimension(P, U, V),
     !
   ;
     true
@@ -112,7 +197,7 @@ find_parameter_dim_rec([0 | EL], L) :-
   find_parameter_dim_rec(EL, L).
 
 find_parameter_dim_rec([E | EL], L) :-
-  debug(units, "Inferring for ~w with ~w", [E, L]),
+  debug(dims, "Inferring for ~w with ~w", [E, L]),
   % all kinetics are Time^(-1) Volume^0
   nb_setval(current_dim_expr, E),
   nb_setval(current_dim_term, E),
@@ -121,7 +206,7 @@ find_parameter_dim_rec([E | EL], L) :-
 
 
 find_parameter_dim(A, Dim, L) :-
-  debug(units, "checking ~w for ~w with ~w", [A, Dim, L]),
+  debug(dims, "checking ~w for ~w with ~w", [A, Dim, L]),
   fail.
 
 find_parameter_dim(X + Y, Dim, L) :-
@@ -297,7 +382,7 @@ write_dimensions([((DT, DV) - P) | L]) :-
   ->
     format("~w has dimension time^(~w).volume^(~w)~n", [P, DT, DV])
   ;
-    debug(units, "~w has unknown dimension: ~w", [P, (DT, DV)]),
+    debug(dims, "~w has unknown dimension: ~w", [P, (DT, DV)]),
     format("~w has unknown dimension~n", [P])
   ),
   write_dimensions(L).