PetabJL integration

.github/workflows/ci.yml

@@ -100,7 +100,7 @@ jobs:
     - name: Install julia
       uses: julia-actions/setup-julia@v1
       with:
-        version: 1.7
+        version: 1.9
 
     - name: Cache
       uses: actions/cache@v3
@@ -113,8 +113,11 @@ jobs:
     - name: Install dependencies
       run: .github/workflows/install_deps.sh
 
+    - name: Install PEtabJL dependencies
+      run: julia -e 'using Pkg; Pkg.add("PEtab"); Pkg.add("OrdinaryDiffEq")'
+
     - name: Run tests
-      timeout-minutes: 20
+      timeout-minutes: 25
       run: tox -e julia
 
     - name: Coverage

doc/example/boehm_JProteomeRes2014/Boehm_validation/Gradient.csv

@@ -0,0 +1,2 @@
+Epo_degradation_BaF3,k_exp_hetero,k_exp_homo,k_imp_hetero,k_imp_homo,k_phos,sd_pSTAT5A_rel,sd_pSTAT5B_rel,sd_rSTAT5A_rel
+-226.01431847863805,0.03263295726678998,0.05074978988074719,-273.3251522968395,-4.9353233689552945e-5,85.33308210815652,-57.82373022138993,6.47910944243548,1.5327970845176904

doc/example/boehm_JProteomeRes2014/Boehm_validation/Hessian.csv

@@ -0,0 +1,10 @@
+Epo_degradation_BaF3,k_exp_hetero,k_exp_homo,k_imp_hetero,k_imp_homo,k_phos,sd_pSTAT5A_rel,sd_pSTAT5B_rel,sd_rSTAT5A_rel
+2348.781845248937,-0.034169780647049425,-103.24202100672939,2772.675796487632,7.194796260472683e-5,-938.3560377189202,996.659728857624,31.143959976782426,13.030712117968559
+-0.034169780647049425,0.07739858784002018,-0.15899501912479183,0.7386693785162849,-4.187803035168086e-8,-0.17134098751539498,-0.2564122376627509,0.03372069395138544,0.07241122210747836
+-103.24202100672939,-0.15899501912479183,132.944794430612,95.49012010074641,-2.508362770407691e-6,-22.250883365424308,11.34374075166351,-3.886786248018129,-7.690665859911962
+2772.675796487632,0.7386693785162849,95.49012010074641,3563.661917245381,1.092950395076315e-5,-1110.2044314822922,1202.406495619125,41.717823566223906,14.58452306207908
+7.194796260472683e-5,-4.187803035168086e-8,-2.508362770407691e-6,1.092950395076315e-5,0.00011364116982094509,-0.00016740979955329535,0.0001089554236616266,0.00011027946613432109,8.045150373694503e-6
+-938.3560377189202,-0.17134098751539498,-22.250883365424308,-1110.2044314822922,-0.00016740979955329535,965.7762449769793,-416.0316410076297,-6.62255869644657,29.6808324053416
+996.659728857624,-0.2564122376627509,11.34374075166351,1202.406495619125,0.0001089554236616266,-416.0316410076297,435.9488576495015,0.0,0.0
+31.143959976782426,0.03372069395138544,-3.886786248018129,41.717823566223906,0.00011027946613432109,-6.62255869644657,0.0,139.82333729022352,0.0
+13.030712117968559,0.07241122210747836,-7.690665859911962,14.58452306207908,8.045150373694503e-6,29.6808324053416,0.0,0.0,162.6019478340014

doc/example/boehm_JProteomeRes2014/Boehm_validation/ObjectiveValue.csv

@@ -0,0 +1,2 @@
+ObjectiveValue
+149.05022585636098

doc/example/boehm_JProteomeRes2014/Boehm_validation/Parameter.csv

@@ -0,0 +1,2 @@
+Epo_degradation_BaF3,k_exp_hetero,k_exp_homo,k_imp_hetero,k_imp_homo,k_phos,sd_pSTAT5A_rel,sd_pSTAT5B_rel,sd_rSTAT5A_rel
+-1.665827601408055,-4.999704893599998,-2.2096987817000167,-1.78600654750001,4.9901140088,4.1977354885,0.5857552705999998,0.8189828191999999,0.49868440400000047

doc/example/conversion_reaction/PEtabJl_module.jl

@@ -0,0 +1,19 @@
+module MyPEtabJlModule
+
+using OrdinaryDiffEq
+using PEtab
+
+pathYaml = "doc/example/conversion_reaction/conversion_reaction.yaml"
+petabModel = readPEtabModel(pathYaml, verbose=true)
+
+# A full list of options for createPEtabODEProblem can be found at https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.setupPEtabODEProblem
+petabProblem = createPEtabODEProblem(
+	petabModel,
+	odeSolverOptions=ODESolverOptions(Rodas5P(), abstol=1e-08, reltol=1e-08, maxiters=Int64(1e4)),
+	gradientMethod=:ForwardDiff,
+	hessianMethod=:ForwardDiff,
+	sparseJacobian=false,
+	verbose=true
+)
+
+end

pypesto/objective/julia/__init__.py

@@ -3,3 +3,4 @@
 ===============
 """
 from .base import JuliaObjective, display_source_ipython
+from .petabJl import PEtabJlObjective

pypesto/objective/julia/petabJl.py

@@ -0,0 +1,214 @@
+"""Interface to PEtab.jl."""
+import logging
+import os
+
+import numpy as np
+
+from .base import JuliaObjective, _read_source
+
+logger = logging.getLogger(__name__)
+
+PEtabProblemJl = ["PEtab.jl::PEtabODEProblem"]
+
+
+class PEtabJlObjective(JuliaObjective):
+    """
+    Wrapper around an objective defined in PEtab.jl.
+
+    Parameters
+    ----------
+    module:
+        Name of the julia module containing the objective.
+    source_file:
+        Julia source file. Defaults to "{module_name}.jl".
+    petab_problem_name:
+        Name of the petab problem variable in the julia module.
+    """
+
+    def __init__(
+        self,
+        module: str,
+        source_file: str = None,
+        petab_problem_name: str = "petabProblem",
+        precompile: bool = True,
+        force_compile: bool = False,
+    ):
+        """Initialize objective."""
+        # lazy imports
+        try:
+            from julia import Main, Pkg  # noqa: F401
+
+            Pkg.activate(".")
+        except ImportError:
+            raise ImportError(
+                "Install PyJulia, e.g. via `pip install pypesto[julia]`, "
+                "and see the class documentation",
+            )
+
+        self.module = module
+        self.source_file = source_file
+        self._petab_problem_name = petab_problem_name
+        if precompile:
+            self.precompile_model(force_compile=force_compile)
+
+        if self.source_file is None:
+            self.source_file = f"{module}.jl"
+
+        # Include module if not already included
+        _read_source(module, source_file)
+
+        petab_jl_problem = self.get(petab_problem_name)
+        self.petab_jl_problem = petab_jl_problem
+
+        # get functions
+        fun = self.petab_jl_problem.computeCost
+        grad = self.petab_jl_problem.computeGradient
+        hess = self.petab_jl_problem.computeHessian
+        x_names = np.asarray(self.petab_jl_problem.θ_estNames)
+
+        # call the super super constructor
+        super(JuliaObjective, self).__init__(
+            fun=fun, grad=grad, hess=hess, x_names=x_names
+        )
+
+    def __getstate__(self):
+        """Get state for pickling."""
+        # if not dumped, dump it via JLD2
+        return {
+            'module': self.module,
+            'source_file': self.source_file,
+            '_petab_problem_name': self._petab_problem_name,
+        }
+
+    def __setstate__(self, state):
+        """Set state from pickling."""
+        for key, value in state.items():
+            setattr(self, key, value)
+        # lazy imports
+        try:
+            from julia import Main  # noqa: F401
+            from julia import Pkg
+
+            Pkg.activate(".")
+        except ImportError:
+            raise ImportError(
+                "Install PyJulia, e.g. via `pip install pypesto[julia]`, "
+                "and see the class documentation",
+            )
+        # Include module if not already included
+        _read_source(self.module, self.source_file)
+
+        petab_jl_problem = self.get(self._petab_problem_name)
+        self.petab_jl_problem = petab_jl_problem
+
+        # get functions
+        fun = self.petab_jl_problem.computeCost
+        grad = self.petab_jl_problem.computeGradient
+        hess = self.petab_jl_problem.computeHessian
+        x_names = np.asarray(self.petab_jl_problem.θ_estNames)
+
+        # call the super super constructor
+        super(JuliaObjective, self).__init__(fun, grad, hess, x_names)
+
+    def __deepcopy__(self, memodict=None):
+        """Deepcopy."""
+        return PEtabJlObjective(
+            module=self.module,
+            source_file=self.source_file,
+            petab_problem_name=self._petab_problem_name,
+            precompile=False,
+        )
+
+    def precompile_model(self, force_compile: bool = False):
+        """
+        Use Julias PrecompilationTools to precompile the relevant code.
+
+        Only needs to be done once, and speeds up Julia loading drastically.
+        """
+        directory = os.path.dirname(self.source_file)
+        # check whether precompilation is necessary, if the directory exists
+        if (
+            os.path.exists(f"{directory}/{self.module}_pre")
+            and not force_compile
+        ):
+            logger.info("Precompilation module already exists.")
+            return None
+        # lazy imports
+        try:
+            from julia import Main  # noqa: F401
+        except ImportError:
+            raise ImportError(
+                "Install PyJulia, e.g. via `pip install pypesto[julia]`, "
+                "and see the class documentation",
+            )
+        # setting up a local project, where the precompilation will be done in
+        from julia import Pkg
+
+        Pkg.activate(".")
+        # create a Project f"{self.module}_pre".
+        try:
+            Pkg.generate(f"{directory}/{self.module}_pre")
+        except Exception:
+            logger.info("Module is already generated. Skipping generate...")
+        # Adjust the precompilation file
+        write_precompilation_module(
+            module=self.module,
+            source_file_orig=self.source_file,
+        )
+        # add a new line at the top of the original module to use the
+        # precompiled module
+        with open(self.source_file, "r") as read_f:
+            if read_f.readline().endswith("_pre\n"):
+                with open("dummy_temp_file.jl", "w+") as write_f:
+                    write_f.write(f"using {self.module}_pre\n\n")
+                    write_f.write(read_f.read())
+                os.remove(self.source_file)
+                os.rename("dummy_temp_file.jl", self.source_file)
+
+        try:
+            Pkg.develop(path=f"{directory}/{self.module}_pre")
+        except Exception:
+            logger.info("Module is already developed. Skipping develop...")
+        Pkg.activate(f"{directory}/{self.module}_pre/")
+        # add dependencies
+        Pkg.add("PrecompileTools")
+        Pkg.add("OrdinaryDiffEq")
+        Pkg.add("PEtab")
+        Pkg.precompile()
+
+
+def write_precompilation_module(module, source_file_orig):
+    """Write the precompilation module for the PEtabJl module."""
+    # read the original source file
+    with open(source_file_orig) as f:
+        lines = np.array(f.readlines())
+    # path to the yaml file
+    yaml_path = "\t".join(lines[["yaml" in line for line in lines]])
+    # packages
+    packages = "\t\t".join(
+        lines[[line.startswith("using ") for line in lines]]
+    )
+    # get everything in between the packages and the end line
+    start = int(np.argwhere([line.startswith("using ") for line in lines])[-1])
+    end = int(np.argwhere([line.startswith("end") for line in lines])[0])
+    petab_loading = "\t\t".join(lines[start:end])
+
+    content = (
+        f"module {module}_pre\n\n"
+        f"using PrecompileTools\n\n"
+        f"# Reduce time for reading a PEtabModel and for "
+        f"building a PEtabODEProblem\n"
+        f"@setup_workload begin\n"
+        f"\t{yaml_path}"
+        f"\t@compile_workload begin\n"
+        f"\t\t{packages}"
+        f"\t\t{petab_loading}"
+        f"\tend\n"
+        f"end\n\n"
+        f"end\n"
+    )
+    # get the directory of the source file
+    directory = os.path.dirname(source_file_orig)
+    # write file
+    with open(f"{directory}/{module}_pre/src/{module}_pre.jl", "w") as f:
+        f.write(content)