Update to HDF5 and PEtab v2

Project.toml

@@ -6,6 +6,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"

src/julia/create_testdata.jl

@@ -1,4 +1,4 @@
-using Lux, CSV, DataFrames, StableRNGs, YAML
+using Lux, CSV, DataFrames, StableRNGs, YAML, HDF5
 
 include(joinpath(@__DIR__, "helper.jl"))
 

src/python/create_testdata.sh

@@ -0,0 +1,17 @@
+# To activate correct conda environment. TODO: Refactor
+eval "$(conda shell.bash hook)"
+conda activate petab_sciml
+export PYTHONPATH="$pwd:$PYTHONPATH"
+
+# Run all PyTorch scripts
+for ((i = 1 ; i < 52 ; i++)); do
+    echo "Test case $i"
+    if [ $i -lt 10 ]; then
+        path="./test_cases/net_00$i/create_testdata/net.py"
+    else
+        path="./test_cases/net_0$i/create_testdata/net.py"
+    fi
+    python $path
+done
+
+exit 0

src/python/helper.py

@@ -1,7 +1,7 @@
 import torch
 import sys
-import pandas as pd
 import os
+import h5py
 
 sys.path.insert(1, os.path.join(os.getcwd(), 'mkstd', "examples", "petab_sciml"))
 from petab_sciml_standard import Input, MLModel, PetabScimlStandard
@@ -19,19 +19,19 @@ def test_nn(net, dir_save, layer_names, dropout=False, atol=1e-3):
         if not layer_names is None:
             for layer_name in layer_names:
                 layer = getattr(net, layer_name)
-                df = pd.read_csv(os.path.join(dir_save, "net_ps_" + str(i) + ".tsv"), delimiter='\t')
-                ps_weight = get_ps_layer(df, layer_name, "weight")
+                ps_h5 = h5py.File(os.path.join(dir_save, "net_ps_" + str(i) + ".h5"), "r")
+                ps_weight = ps_h5[layer_name]["weight"][:]
                 with torch.no_grad():
-                    layer.weight[:] = ps_weight
+                    layer.weight[:] = torch.from_numpy(ps_weight)
                 if hasattr(layer, "bias") and (not layer.bias is None):
-                    ps_bias = get_ps_layer(df, layer_name, "bias")
+                    ps_bias = ps_h5[layer_name]["bias"][:]
                     with torch.no_grad():
-                        layer.bias[:] = ps_bias
+                        layer.bias[:] = torch.from_numpy(ps_bias)
 
-        df_input = pd.read_csv(os.path.join(dir_save, "net_input_" + str(i) + ".tsv"), delimiter='\t')
-        df_output = pd.read_csv(os.path.join(dir_save, "net_output_" + str(i) + ".tsv"), delimiter='\t')
-        input = read_array(df_input)
-        output_ref = read_array(df_output)
+        input_h5 = h5py.File(os.path.join(dir_save, "net_input_" + str(i) + ".h5"))
+        output_h5 = h5py.File(os.path.join(dir_save, "net_output_" + str(i) + ".h5"))
+        input = torch.from_numpy(input_h5["input"][:])
+        output_ref = torch.from_numpy(output_h5["output"][:])
         if dropout == False:    
             output = net.forward(input)
         else:
@@ -40,36 +40,3 @@ def test_nn(net, dir_save, layer_names, dropout=False, atol=1e-3):
                 output += net.forward(input)
             output /= 50000
         torch.testing.assert_close(output_ref, output, atol=atol, rtol=0.0)
-
-
-def extract_numbers(series):
-    out = [[int(part) for part in s.split('_') if part] for s in series]
-    return out
-
-def get_dim(ix):
-    max_values = [max(column) + 1 for column in zip(*ix)]
-    return max_values
-
-def get_ps_layer(df, layer_name, ps_name):
-    i_layer = df.loc[:, "parameterId"].str.startswith("net_" + layer_name + "_" + ps_name)
-    df_layer = df.loc[i_layer, :]
-    df_layer.reset_index(drop=True, inplace=True)
-    ix = df_layer.loc[:, "parameterId"]
-    ix = extract_numbers(ix.str.replace("net_" + layer_name + "_" + ps_name, ""))
-    dims = get_dim(ix)
-    out = torch.ones(*dims)
-    for i in range(df_layer.shape[0]):
-        _ix = ix[i]
-        out[*_ix] = df_layer.loc[i, "value"]
-
-    return out
-
-def read_array(df):
-    ix = df.loc[:, "ix"].astype("string")
-    ix = ix.apply(lambda x: [int(i) for i in x.split(';')]).tolist()
-    dims = get_dim(ix)
-    out = torch.ones(*dims)
-    for i in range(df.shape[0]):
-        _ix = ix[i]
-        out[*_ix] = df.loc[i, "value"]
-    return out

test_cases/001/create_testdata/testdata.jl

@@ -9,8 +9,8 @@ Random.seed!(123)
 
 function compute_nllh(x, oprob::ODEProblem, solver, measurements::DataFrame; abstol = 1e-9,
                       reltol = 1e-9)::Real
-    mprey = measurements[measurements[!, :observableId] .== "prey", :measurement]
-    mpredator = measurements[measurements[!, :observableId] .== "predator", :measurement]
+    mprey = measurements[measurements[!, :observableId] .== "prey_o", :measurement]
+    mpredator = measurements[measurements[!, :observableId] .== "predator_o", :measurement]
     tsave = unique(measurements.time)
 
     _oprob = remake(oprob, p = x)
@@ -32,8 +32,7 @@ measurements = CSV.read(joinpath(@__DIR__, "..", "petab", "measurements.tsv"), D
 # Objective function
 x = deepcopy(p_ode)
 # Read neural net parameters, and assign to x
-df_ps_nn = CSV.read(joinpath(@__DIR__, "..", "petab", "parameters_nn.tsv"), DataFrame)
-set_ps_net!(x.p_net1, df_ps_nn, :net1, nn_model)
+set_ps_net!(x.p_net1, joinpath(@__DIR__, "..", "petab", "net1_ps.hf5"), nn_model)
 
 ## Compute model values
 _f = (x) -> compute_nllh(x, oprob_nn, Vern9(), measurements; abstol = 1e-12, reltol = 1e-12)
@@ -46,13 +45,15 @@ sol = solve(oprob_nn, Vern9(), abstol = 1e-9, reltol = 1e-9,
             saveat = unique(measurements.time))
 simulated_values = vcat(sol[1, :], sol[2, :])
 
-## Write values for saving to file
+## Write test values
 # YAML problem file
 solutions = Dict(:llh => llh,
                  :tol_llh => 1e-3,
                  :tol_grad_llh => 1e-1,
                  :tol_simulations => 1e-3,
-                 :grad_llh_files => ["grad_llh.tsv"],
+                 :grad_llh_files => Dict(
+                    :mech => "grad_mech.tsv",
+                    :net1 => "grad_net1.hf5"),
                  :simulation_files => ["simulations.tsv"])
 YAML.write_file(joinpath(@__DIR__, "..", "solutions.yaml"), solutions)
 # Simulated values
@@ -61,27 +62,33 @@ rename!(simulations_df, "measurement" => "simulation")
 simulations_df.simulation .= simulated_values
 CSV.write(joinpath(@__DIR__, "..", "simulations.tsv"), simulations_df, delim = '\t')
 # Gradient values
-df_net = deepcopy(df_ps_nn)
-df_net.value .= llh_grad.p_net1
 df_mech = DataFrame(parameterId = ["alpha", "delta", "beta"],
                     value = llh_grad[1:3])
-df_grad = vcat(df_mech, df_net)
-CSV.write(joinpath(@__DIR__, "..", "grad_llh.tsv"), df_grad, delim = '\t')
-# Write problem yaml
+CSV.write(joinpath(@__DIR__, "..", "grad_mech.tsv"), df_mech, delim = '\t')
+nn_ps_to_h5(nn_model, llh_grad.p_net1, joinpath(@__DIR__, "..", "grad_net1.hf5"))
+
+## Write PEtabProblem files
+mapping_table = DataFrame(
+    Dict("petab.MODEL_ENTITY_ID" => ["net1.input1", "net1.input2", "net1.output1"],
+         "petab.PETAB_ENTITY_ID" => ["prey", "predator", "gamma"]))
+CSV.write(joinpath(@__DIR__, "..", "petab", "mapping_table.tsv"), mapping_table; delim = '\t')
 problem_yaml = Dict(
-    :format_version => 1,
+    :format_version => 2,
     :parameter_file => "parameters_ude.tsv",
     :problems => [Dict(
         :condition_files => ["conditions.tsv"],
         :measurement_files => ["measurements.tsv"],
         :observable_files => ["observables.tsv"],
-        :sbml_files => ["lv.xml"],
-        :mapping_tables => "mapping_table.tsv")],
+        :model_files => Dict(
+            :language => "sbml",
+            :location => "lv.xml"),
+        :mapping_files => ["mapping_table.tsv"])],
     :extensions => Dict(
         :petab_sciml => Dict(
-            :net_files => ["net1.yaml"],
-            :hybridization => Dict(
-                :net1 => Dict(
+            :net1 => Dict(
+                :file => "net1.yaml",
+                :parameters => "net1_ps.h5",
+                :hybridization => Dict(
                     :input => "ode",
                     :output => "ode")))))
 YAML.write_file(joinpath(@__DIR__, "..", "petab", "problem_ude.yaml"), problem_yaml)

test_cases/001/create_testdata/train_net.jl

@@ -54,6 +54,4 @@ sol1 = solve(prob, OptimizationOptimisers.Adam(0.001), maxiters = 10000)
 x0 .= sol1.u
 sol2 = solve(prob, Optimization.LBFGS(), maxiters = 4000)
 
-# Write neural-net parameters to file
-ps_df = nn_ps_to_tidy(nn_model, sol2.u, :net1)
-CSV.write(joinpath(@__DIR__, "..", "petab", "parameters_nn.tsv"), ps_df, delim = '\t')
+nn_ps_to_h5(nn_model, sol2.u, joinpath(@__DIR__, "..", "petab", "net1_ps.hf5"))

test_cases/001/grad_mech.tsv

@@ -0,0 +1,4 @@
+parameterId	value
+alpha	1427.7245805979464
+delta	-1089.8948747895438
+beta	-3711.6409538321677

test_cases/001/petab/mapping_table.tsv

@@ -1,4 +1,4 @@
-netId	ioId	ioValue
-net1	input1	prey
-net1	input2	predator
-net1	output1	gamma
+petab.MODEL_ENTITY_ID	petab.PETAB_ENTITY_ID
+net1.input1	prey
+net1.input2	predator
+net1.output1	gamma

test_cases/001/petab/measurements.tsv

@@ -1,21 +1,21 @@
 observableId	simulationConditionId	measurement	time
-prey	cond1	0.2400576308360028	1.0
-prey	cond1	0.42814983671389123	2.0
-prey	cond1	1.551119400580878	3.0
-prey	cond1	5.473290024101726	4.0
-prey	cond1	3.0925928860569543	5.0
-prey	cond1	0.2066446962341261	6.0
-prey	cond1	0.2754091191501447	7.0
-prey	cond1	0.8363612728837162	8.0
-prey	cond1	3.0547210667300426	9.0
-prey	cond1	8.886774826412893	10.0
-predator	cond1	0.9810266619495472	1.0
-predator	cond1	0.17933792963889872	2.0
-predator	cond1	0.19548793632652334	3.0
-predator	cond1	0.09809325439896757	4.0
-predator	cond1	6.661319440631508	5.0
-predator	cond1	1.98105039531585	6.0
-predator	cond1	0.3984814127945589	7.0
-predator	cond1	0.10102862759201464	8.0
-predator	cond1	0.12027022434787174	9.0
-predator	cond1	1.1961696974794869	10.0
+prey_o	cond1	0.2400576308360028	1.0
+prey_o	cond1	0.42814983671389123	2.0
+prey_o	cond1	1.551119400580878	3.0
+prey_o	cond1	5.473290024101726	4.0
+prey_o	cond1	3.0925928860569543	5.0
+prey_o	cond1	0.2066446962341261	6.0
+prey_o	cond1	0.2754091191501447	7.0
+prey_o	cond1	0.8363612728837162	8.0
+prey_o	cond1	3.0547210667300426	9.0
+prey_o	cond1	8.886774826412893	10.0
+predator_o	cond1	0.9810266619495472	1.0
+predator_o	cond1	0.17933792963889872	2.0
+predator_o	cond1	0.19548793632652334	3.0
+predator_o	cond1	0.09809325439896757	4.0
+predator_o	cond1	6.661319440631508	5.0
+predator_o	cond1	1.98105039531585	6.0
+predator_o	cond1	0.3984814127945589	7.0
+predator_o	cond1	0.10102862759201464	8.0
+predator_o	cond1	0.12027022434787174	9.0
+predator_o	cond1	1.1961696974794869	10.0

test_cases/001/petab/observables.tsv

@@ -1,3 +1,3 @@
 observableId	observableName	observableFormula	noiseFormula	observableTransformation	noiseDistribution
-prey		prey	0.05	lin	normal
-predator		predator	0.05	lin	normal
+prey_o		prey	0.05	lin	normal
+predator_o		predator	0.05	lin	normal

test_cases/001/petab/parameters_ude.tsv

@@ -2,4 +2,4 @@ parameterId	parameterName	parameterScale	lowerBound	upperBound	nominalValue	esti
 alpha	alpha	lin	0	15	1.3	1
 beta	beta	lin	0	15	0.9	1
 delta	delta	lin	0	15	1.8	1
-p_net1	net1	lin	-inf	inf	0.0	1
+net1	net1	lin	-inf	inf	0.0	1