diff --git a/mri/run_mri.py b/mri/run_mri.py
index 086fbe2c8453de6dd5a318b9fb5f72109abaf96d..c9e7251b38c0f52c056dfa18e02388c75973e7d3 100755
--- a/mri/run_mri.py
+++ b/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)
diff --git a/mri/uq_mri.py b/mri/uq_mri.py
new file mode 100755
index 0000000000000000000000000000000000000000..58cdfa5445e238e6906a3ba1435abda4412fd864
--- /dev/null
+++ b/mri/uq_mri.py
@@ -0,0 +1,313 @@
+#!/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)