[mri] add the scripts.

mri/run_mri.py

@@ -2,60 +2,79 @@
 
 from __future__ import division, print_function
 
-import os
-import sys
+
 import numpy as np
-import argparse
+import multiprocessing as mp
 
 # add location where to find the forward model wrapper
+#import sys
+#import os
 #sys.path.append(os.path.abspath("./"))
 import pymsbern
 
-
-# the forward model
-def run_mri(random_params):
-    paramnames = sorted(["MRSignal.T1TissuePreContrast", "InflowCurve.a", "InflowCurve.b", "InflowCurve.tp",
-                         "VesselWall.DiffusionCoefficient", "MRSignal.KappaB", "MRSignal.KappaT"])
+def run_mri(x):
+    """
+    This function runs the forward model with the given parameters x.
+
+    Parameters
+    ----------
+    x : array of shape (n_params)
+        Input parameters.
+
+    Returns
+    -------
+    outputs : 2d-array of shape (2, n_timestep)
+        An array containing the simulation results.
+
+    """
+    x = np.atleast_2d(x)
+    paramnames = ["InflowCurve.a",
+                  "InflowCurve.b",
+                  "InflowCurve.tp",
+                  "MRSignal.KappaB",
+                  "MRSignal.KappaT",
+                  "MRSignal.T1TissuePreContrast",
+                  "VesselWall.DiffusionCoefficient"]
 
     # the parameters passed to the model
     params = {}
 
     # convert to dictionary again
     for i, name in enumerate(paramnames):
-        params[name] = str(random_params[i])
+        params[name] = str(x[0,i])
 
     # create output parameters and dump parameters
-    params["Problem.LesionIntensity"] = "1" # healthy: 0, sick: 1
-    params["Problem.SolveFlow"] = "false" # only transport
-    params["VesselWall.FiltrationCoefficient"] = "9.746174233413226E-011" # Lp
-    params["MRSignal.T2TissuePreContrast"] = "8.000000000000000E-002" # seconds
+    params["Problem.LesionIntensity"] = "1"  # healthy: 0, sick: 1
+    params["Problem.SolveFlow"] = "false"  # only transport
+    params["VesselWall.FiltrationCoefficient"] = "9.746174233413226E-011"
+    params["MRSignal.T2TissuePreContrast"] = "8.000000000000000E-002"
     # convert back from log scale
-    params["VesselWall.DiffusionCoefficient"] = str(10**(-float(params["VesselWall.DiffusionCoefficient"])))
+    diff_coeff = str(10**(-float(params["VesselWall.DiffusionCoefficient"])))
+    params["VesselWall.DiffusionCoefficient"] = diff_coeff
 
     params["Vessel.Grid.File"] = "./data/ratbrain.dgf"
 
-    verbose = False
-    data = np.array(pymsbern.run_forward_model(params, "./data/msbern.input", verbose))
+    data = np.array(pymsbern.run_forward_model(params,
+                                               "./data/msbern.input",
+                                               verbose=False))
+    # Prepare the outputs
+    t = [i*1.4 for i in range(0, 80)]
+    outputs = np.vstack((t, data))
 
-    return data
+    return outputs
 
 
 # main routine
 if __name__ == '__main__':
 
-    # choose an initial parameter set
-    initial = {}
-    initial["MRSignal.T1TissuePreContrast"] = 1.4 # seconds
-    initial["InflowCurve.a"] = 3.812619924254471E+001 # scaling params mol*s/m3
-    initial["InflowCurve.b"] = 1.934938693976044E+000 # baseline mol/m3
-    initial["InflowCurve.tp"] = 5.000000000000000E+000 # time to peak in s
-    initial["VesselWall.DiffusionCoefficient"] = 14.0 # D wall
-    initial["MRSignal.KappaB"] = 1.502480256948131E+001
-    initial["MRSignal.KappaT"] = 0.1
-
-    # convert to plain array for emcee
-    theta = np.array([v for k,v in sorted(initial.items())])
-
-    print(theta)
-    sim = run_model(theta)
-    print(sim)
+    # Define an initial parameter set
+    theta = np.array([
+        3.812619924254471E+001,  # "InflowCurve.a" scaling params mol*s/m3
+        1.934938693976044E+000,  # "InflowCurve.b" baseline mol/m3
+        5.000000000000000E+000,  # "InflowCurve.tp" time to peak in s
+        1.502480256948131E+001,  # "MRSignal.KappaB"
+        0.1,  # "MRSignal.KappaT"
+        1.4,  # "MRSignal.T1TissuePreContrast" in seconds
+        14.0,  # "VesselWall.DiffusionCoefficient" D wall
+        ])
+    sim = run_mri(theta)

mri/uq_mri.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+This test shows a surrogate-assisted Bayesian calibration of a time dependent
+    mri model.
+
+Author: Farid Mohammadi, M.Sc.
+E-Mail: farid.mohammadi@iws.uni-stuttgart.de
+Department of Hydromechanics and Modelling of Hydrosystems (LH2)
+Institute for Modelling Hydraulic and Environmental Systems (IWS), University
+of Stuttgart, www.iws.uni-stuttgart.de/lh2/
+Pfaffenwaldring 61
+70569 Stuttgart
+
+Created on Mon Feb 14 2022
+
+"""
+
+import numpy as np
+import pandas as pd
+import scipy.stats as stats
+import sys
+import joblib
+# import bayesvalidrox
+sys.path.append("../src/bayesvalidrox/")
+from pylink.pylink import PyLinkForwardModel
+from surrogate_models.inputs import Input
+from surrogate_models.surrogate_models import MetaModel
+from post_processing.post_processing import PostProcessing
+from bayes_inference.bayes_inference import BayesInference, Discrepancy
+
+import matplotlib
+matplotlib.use('agg')
+
+if __name__ == "__main__":
+
+    # =====================================================
+    # =============   COMPUTATIONAL MODEL  ================
+    # =====================================================
+    Model = PyLinkForwardModel()
+
+    Model.link_type = 'Function'
+    Model.py_file = 'run_mri'
+    Model.name = 'mri'
+
+    Model.Output.names = ['NMR signal']  # (normalized to baseline)
+
+    # Pass data as dict for Bayesian inversion
+    Model.observations = {}
+    Model.observations['Time [s]'] = [i*1.4 for i in range(0, 80)]
+    Model.observations['NMR signal'] = np.loadtxt('data/msbern-sick.obs')
+
+    # =====================================================
+    # =========   PROBABILISTIC INPUT MODEL  ==============
+    # =====================================================
+    # Define the uncertain parameters with their mean and
+    # standard deviation
+    Inputs = Input()
+
+    # "InflowCurve.a" scaling params mol*s/m3 [10, 100]
+    Inputs.addMarginals()
+    Inputs.Marginals[0].name = '$a$'
+    Inputs.Marginals[0].dist_type = 'uniform'
+    Inputs.Marginals[0].parameters = [10, 100]
+
+    # "InflowCurve.b" baseline mol/m3 [0, 2]
+    Inputs.addMarginals()
+    Inputs.Marginals[1].name = '$b$'
+    Inputs.Marginals[1].dist_type = 'uniform'
+    Inputs.Marginals[1].parameters = [0, 2]
+
+    # "InflowCurve.tp" time to peak in s [1, 8]
+    Inputs.addMarginals()
+    Inputs.Marginals[2].name = '$t_p$'
+    Inputs.Marginals[2].dist_type = 'uniform'
+    Inputs.Marginals[2].parameters = [1, 8]
+
+    # "MRSignal.KappaB" [0, 100]
+    Inputs.addMarginals()
+    Inputs.Marginals[3].name = '$\\kappa_B$'
+    Inputs.Marginals[3].dist_type = 'uniform'
+    Inputs.Marginals[3].parameters = [0, 100]
+
+    # "MRSignal.KappaT" [0, 50]
+    Inputs.addMarginals()
+    Inputs.Marginals[4].name = '$\\kappa_T$'
+    Inputs.Marginals[4].dist_type = 'uniform'
+    Inputs.Marginals[4].parameters = [0, 50]
+
+    # "MRSignal.T1TissuePreContrast" in seconds [0.8, 3.0]
+    Inputs.addMarginals()
+    Inputs.Marginals[5].name = '$T_{1, \\textrm{pre}}$'
+    Inputs.Marginals[5].dist_type = 'uniform'
+    Inputs.Marginals[5].parameters = [0.8, 3.0]
+
+    # "VesselWall.DiffusionCoefficient" D wall [6, 10]
+    Inputs.addMarginals()
+    Inputs.Marginals[6].name = '$Log_{10}(D_{\\omega})$'
+    Inputs.Marginals[6].dist_type = 'uniform'
+    Inputs.Marginals[6].parameters = [6, 10]
+
+    # =====================================================
+    # ==========  DEFINITION OF THE METAMODEL  ============
+    # =====================================================
+    MetaModelOpts = MetaModel(Inputs)
+
+    # Select if you want to preserve the spatial/temporal depencencies
+    MetaModelOpts.dim_red_method = 'PCA'
+    MetaModelOpts.var_pca_threshold = 99.95
+    #MetaModelOpts.n_pca_components = 12
+
+    # Select your metamodel method
+    # 1) PCE (Polynomial Chaos Expansion) 2) aPCE (arbitrary PCE)
+    # 3) GPE (Gaussian Process Emulator)
+    MetaModelOpts.meta_model_type = 'aPCE'
+
+    # ------------------------------------------------
+    # ------------- PCE Specification ----------------
+    # ------------------------------------------------
+    # Select the sparse least-square minimization method for
+    # the PCE coefficients calculation:
+    # 1)OLS: Ordinary Least Square  2)BRR: Bayesian Ridge Regression
+    # 3)LARS: Least angle regression  4)ARD: Bayesian ARD Regression
+    # 5)FastARD: Fast Bayesian ARD Regression
+    # 6)VBL: Variational Bayesian Learning
+    # 7)EBL: Emperical Bayesian Learning
+    MetaModelOpts.pce_reg_method = 'FastARD'
+
+    # Specify the max degree to be compared by the adaptive algorithm:
+    # The degree with the lowest Leave-One-Out cross-validation (LOO)
+    # error (or the highest score=1-LOO) estimator is chosen.
+    # pce_deg accepts degree as a scalar or a range.
+    MetaModelOpts.pce_deg = 8  # np.arange(7)
+
+    # q-quasi-norm 0<q<1 (default=1)
+    MetaModelOpts.pce_q_norm = 0.85
+
+    # Print summary of the regression results
+    # MetaModelOpts.verbose = True
+
+    # ------------------------------------------------
+    # ------ Experimental Design Configuration -------
+    # ------------------------------------------------
+    MetaModelOpts.add_ExpDesign()
+
+    # One-shot (normal) or Sequential Adaptive (sequential) Design
+    MetaModelOpts.ExpDesign.Method = 'normal'
+    MetaModelOpts.ExpDesign.n_init_samples = 400
+
+    # Sampling methods
+    # 1) random 2) latin_hypercube 3) sobol 4) halton 5) hammersley 6) korobov
+    # 7) chebyshev(FT) 8) grid(FT) 9) nested_grid(FT) 10)user
+    MetaModelOpts.ExpDesign.sampling_method = 'user'
+
+    # Provide the experimental design object with a hdf5 file
+    MetaModelOpts.ExpDesign.hdf5_file = 'ExpDesign_mri.hdf5'
+
+    # Sequential experimental design (needed only for sequential ExpDesign)
+    MetaModelOpts.ExpDesign.n_new_samples = 1
+    MetaModelOpts.ExpDesign.n_max_samples = 15  # 150
+    MetaModelOpts.ExpDesign.mod_LOO_threshold = 1e-16
+
+    # Defining the measurement error, if it's known a priori
+    obsData = pd.DataFrame(Model.observations, columns=Model.Output.names)
+    DiscrepancyOpts = Discrepancy('')
+    DiscrepancyOpts.Type = 'Gaussian'
+    DiscrepancyOpts.Parameters = obsData**2
+    MetaModelOpts.Discrepancy = DiscrepancyOpts
+
+    # ------------------------------------------------
+    # ------- Sequential Design configuration --------
+    # ------------------------------------------------
+    MetaModelOpts.adapt_verbose = True
+    # 1) None 2) 'equal' 3)'epsilon-decreasing' 4) 'adaptive'
+    MetaModelOpts.ExpDesign.tradeoff_scheme = None
+    # MetaModelOpts.ExpDesign.nReprications = 20
+    # -------- Exploration ------
+    # 1)'Voronoi' 2)'random' 3)'latin_hypercube' 4)'LOOCV' 5)'dual annealing'
+    MetaModelOpts.ExpDesign.explore_method = 'latin_hypercube'
+
+    # Use when 'dual annealing' chosen
+    MetaModelOpts.ExpDesign.MaxFunItr = 200
+
+    # Use when 'Voronoi' or 'random' or 'latin_hypercube' chosen
+    MetaModelOpts.ExpDesign.n_canddidate = 5000
+    MetaModelOpts.ExpDesign.n_cand_groups = 4  # 8
+
+    # -------- Exploitation ------
+    # 1)'BayesOptDesign' 2)'BayesActDesign' 3)'VarOptDesign' 4)'alphabetic'
+    # 5)'Space-filling'
+    MetaModelOpts.ExpDesign.exploit_method = 'VarOptDesign'
+
+    # BayesOptDesign -> when data is available
+    # 1)DKL (Kullback-Leibler Divergence) 2)DPP (D-Posterior-percision)
+    # 3)APP (A-Posterior-percision)  # ['DKL', 'BME', 'infEntropy']
+    # MetaModelOpts.ExpDesign.util_func = 'DKL'
+
+    # VarBasedOptDesign -> when data is not available
+    # Only with Vornoi >>> 1)Entropy 2)EIGF, 3)LOOCV
+    # or a combination as a list
+    MetaModelOpts.ExpDesign.util_func = 'Entropy'
+
+    # alphabetic
+    # 1)D-Opt (D-Optimality) 2)A-Opt (A-Optimality)
+    # 3)K-Opt (K-Optimality) or a combination as a list
+    # MetaModelOpts.ExpDesign.util_func = 'D-Opt'
+
+    # For calculation of validation error for SeqDesign
+    # valid_set = np.load(f"data/Prior_{ndim}.npy")
+    # valid_set_outputs = np.load(f"data/origModelOutput_{ndim}.npy")
+    # MetaModelOpts.valid_samples = valid_set
+    # MetaModelOpts.valid_model_runs = {'NMR signal': valid_set_outputs}
+
+    # >>>>>>>>>>>>>>>>>>>>>> Build Surrogate <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+    # Train the meta model
+    PCEModel = MetaModelOpts.create_metamodel(Model)
+
+    # Save PCE models
+    with open(f'PCEModel_{Model.name}.pkl', 'wb') as output:
+        joblib.dump(PCEModel, output, 2)
+
+    # Load PCEModel
+    # with open(f'PCEModel_{Model.name}.pkl', 'rb') as input:
+    #     PCEModel = joblib.load(input)
+    print(PCEModel.score_dict)
+    # =====================================================
+    # =========  POST PROCESSING OF METAMODELS  ===========
+    # =====================================================
+    PostPCE = PostProcessing(PCEModel)
+
+    # Compute the moments and compare with the Monte-Carlo reference
+    PostPCE.plotMoments()
+
+    # Plot the evolution of the KLD,BME, and Modified LOOCV error
+    if MetaModelOpts.ExpDesign.Method == 'sequential':
+        PostPCE.seqDesignDiagnosticPlots()
+
+    # Plot the sobol indices
+    sobol_cell, total_sobol = PostPCE.sobolIndicesPCE()
+
+    # Plot to check validation visually.
+    PostPCE.validMetamodel(nValidSamples=3)
+
+    # Compute and print RMSE error
+    PostPCE.accuracyCheckMetaModel(nSamples=10)
+    sys.exit()
+    # =====================================================
+    # ========  Bayesian inference with Emulator ==========
+    # =====================================================
+    BayesOpts = BayesInference(PCEModel)
+    BayesOpts.emulator = True
+    BayesOpts.PlotPostPred = True
+
+    # BayesOpts.selectedIndices = [0, 3, 5,  7, 9]
+    # BME Bootstrap
+    # BayesOpts.Bootstrap = True
+    # BayesOpts.BootstrapItrNr = 500
+    # BayesOpts.BootstrapNoise = 100
+
+    # Bayesian cross validation
+    # BayesOpts.BayesLOOCV = True
+
+    # Select the inference method
+    BayesOpts.SamplingMethod = 'MCMC'
+    BayesOpts.MCMCnwalkers = 30
+    # BayesOpts.MCMCinitSamples = np.zeros((1,10))+1e-3*np.random.randn(200, 10)
+    # BayesOpts.MCMCnSteps = 1000
+    import emcee
+    BayesOpts.MCMCmoves = emcee.moves.KDEMove()
+    # BayesOpts.MultiProcessMCMC = True
+
+    # ----- Define the discrepancy model -------
+    BayesOpts.MeasurementError = obsData
+
+    # # -- (Option B) --
+    DiscrepancyOpts = bayesvalidrox.Discrepancy('')
+    DiscrepancyOpts.Type = 'Gaussian'
+    DiscrepancyOpts.Parameters = obsData**2
+    # DiscrepancyOpts.Parameters = {'Z':np.ones(10)}
+    BayesOpts.Discrepancy = DiscrepancyOpts
+
+    # -- (Option C) --
+    # DiscOutputOpts = Input()
+    # # # OutputName = 'Z'
+    # DiscOutputOpts.addMarginals()
+    # DiscOutputOpts.Marginals[0].Name = '$\sigma^2_{\epsilon}$'
+    # DiscOutputOpts.Marginals[0].DistType = 'unif'
+    # DiscOutputOpts.Marginals[0].Parameters =  [0, 10]
+    # BayesOpts.Discrepancy = {'known': DiscrepancyOpts,
+    #                           'infer': Discrepancy(DiscOutputOpts)}
+
+    # BayesOpts.BiasInputs = {'Z':np.arange(0, 10, 1.).reshape(-1,1) / 9}
+    # DiscOutputOpts = Input()
+    # # OutputName = 'lambda'
+    # DiscOutputOpts.addMarginals()
+    # DiscOutputOpts.Marginals[0].Name = '$\lambda$'
+    # DiscOutputOpts.Marginals[0].DistType = 'unif'
+    # DiscOutputOpts.Marginals[0].Parameters = [0, 1]
+
+    # # OutputName = 'sigma_f'
+    # DiscOutputOpts.addMarginals()
+    # DiscOutputOpts.Marginals[1].Name = '$\sigma_f$'
+    # DiscOutputOpts.Marginals[1].DistType = 'unif'
+    # DiscOutputOpts.Marginals[1].Parameters = [0, 1e-4]
+    # BayesOpts.Discrepancy = Discrepancy(DiscOutputOpts)
+    # BayesOpts.Discrepancy = {'known': DiscrepancyOpts,
+    #                           'infer': Discrepancy(DiscOutputOpts)}
+    # Start the calibration/inference
+    Bayes_PCE = BayesOpts.create_Inference()
+
+    # Save class objects
+    with open(f'Bayes_{Model.name}.pkl', 'wb') as output:
+        joblib.dump(Bayes_PCE, output, 2)