Start adding mcmc diagnotstic plots

.github/worflows/CI.yml

@@ -0,0 +1,31 @@
+name: CI
+
+on:
+  push:
+    branches: [main, master]
+  pull_request:
+    branches: [main, master]
+
+jobs:
+  build:
+    runs-on: ${{ matrix.os }}
+    continue-on-error: ${{ matrix.os == 'windows-latest' }}
+    strategy:
+      matrix:
+        os: [ubuntu-22.04, macos-latest, windows-latest]
+        python-version: [3.6, 3.7, 3.8, 3.9, '3.10']
+
+    steps:
+    - uses: actions/checkout@v3
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v4
+      with:
+        python-version: ${{ matrix.python-version }}
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install -r requirements.txt
+        pip install .
+    - name: Run Tests
+      run: |
+        python -m unittest --verbose

CHANGES.md

@@ -0,0 +1,5 @@
+# impala 0.1.9999
+* refactor to new PEP standards
+
+# impala 0.1
+* initial version of package

README.md

@@ -13,6 +13,10 @@ pip install git+https://github.com/lanl/impala.git
 
 ## References
 
+## Tests
+```console
+> python -m unittest
+```
 
 ************
 

impala/__init__.py

@@ -0,0 +1,3 @@
+from .physical_models_vec import *
+from .impala_noProbit_emu import *
+from .models_withLik import *

impala/superCal/impala_noProbit_emu.py → impala/impala_noProbit_emu.py

@@ -8,12 +8,7 @@
 from scipy.special import erf, erfinv, gammaln
 from scipy.stats import invwishart
 from numpy.linalg import cholesky, slogdet
-#from itertools import repeat
-#import multiprocessing as mp
-#import pandas as pd
-#import impala.superCal.pbar as pbar
 from . import pbar
-#import pbar 
 np.seterr(under='ignore')
 
 # no probit tranform for hierarchical and DP versions

impala/superCal/models_withLik.py → impala/models_withLik.py

@@ -1,8 +1,7 @@
 import numpy as np
 import pyBASS as pb
-#import pyBayesPPR as pbppr
-#import physical_models_vec as pm_vec
-from impala import physics as pm_vec
+import pyBayesPPR as pbppr
+import impala.physical_models_vec as pm_vec
 from itertools import cycle
 from scipy.interpolate import interp1d
 from scipy.linalg import cho_factor, cho_solve, cholesky
@@ -233,7 +232,6 @@ def step(self):
         self.ii = np.random.choice(range(self.nmcmc), 1).item()
         self.emu_vars = self.mod_s2[self.ii]
         return
-    #@profile
     def discrep_sample(self, yobs, pred, cov, itemp):
         #if self.nd>0:
         S = np.linalg.inv(
@@ -244,7 +242,6 @@ def discrep_sample(self, yobs, pred, cov, itemp):
         discrep_vars = chol_sample(S @ m, S/itemp)
         #self.discrep = self.D @ self.discrep_vars
         return discrep_vars
-    #@profile
     def eval(self, parmat, pool = None, nugget=False):
         """
         parmat : ~
@@ -259,32 +256,23 @@ def eval(self, parmat, pool = None, nugget=False):
             return np.concatenate([pred[np.ix_(np.arange(i, nrep*self.nexp, self.nexp), np.where(self.exp_ind==i)[0])] for i in range(self.nexp)], 1)
             # this is evaluating all experiments for all thetas, which is overkill
 
-    #@profile
     def llik(self, yobs, pred, cov):
         vec = yobs - pred 
         out = -.5*(cov['ldet'] + vec.T @ cov['inv'] @ vec)
         return out
-    #@profile
     def lik_cov_inv(self, s2vec):
         vec = self.trunc_error_var + s2vec
-        #mat = np.diag(vec) + self.basis @ np.diag(self.emu_vars) @ self.basis.T
-        #inv = np.linalg.inv(mat)
-        #ldet = np.linalg.slogdet(mat)[1]
-        #out = {'inv' : inv, 'ldet' : ldet}
         Ainv = np.diag(1/vec)
         Aldet = np.log(vec).sum()
         out = self.swm(Ainv, self.basis, np.diag(1/self.emu_vars), self.basis.T, Aldet, np.log(self.emu_vars).sum())
         return out
-    #@profile
     def chol_solve(self, x):
         mat = cho_factor(x)
         ldet = 2 * np.sum(np.log(np.diag(mat[0])))
-        ##la.dpotri(mat, overwrite_c=True) # overwrites mat with original matrix inverse, but not correct
         #inv = cho_solve(mat, np.eye(x.shape[0])) # correct, but slower for small dimension
         inv = np.linalg.inv(x)
         out = {'inv' : inv, 'ldet' : ldet}
         return out
-    #@profile
     def swm(self, Ainv, U, Cinv, V, Aldet, Cldet): # sherman woodbury morrison (A+UCV)^-1 and |A+UCV|
         in_mat = self.chol_solve(Cinv + V @ Ainv @ U)
         inv = Ainv - Ainv @ U @ in_mat['inv'] @ V @ Ainv
@@ -339,18 +327,14 @@ def step(self):
         self.ii = np.random.choice(range(self.nmcmc), 1).item()
         self.emu_vars = self.mod_s2[self.ii]
         return
-    #@profile
     def discrep_sample(self, yobs, pred, cov, itemp):
-        #if self.nd>0:
         S = np.linalg.inv(
             np.eye(self.nd) / self.discrep_tau 
             + self.D.T @ cov['inv'] @ self.D
             )
         m = self.D.T @ cov['inv'] @ (yobs - pred)
         discrep_vars = chol_sample(S @ m, S/itemp)
-        #self.discrep = self.D @ self.discrep_vars
         return discrep_vars
-    #@profile
     def eval(self, parmat, pool = None, nugget=False):
         """
         parmat : ~
@@ -365,32 +349,23 @@ def eval(self, parmat, pool = None, nugget=False):
             return np.concatenate([pred[np.ix_(np.arange(i, nrep*self.nexp, self.nexp), np.where(self.exp_ind==i)[0])] for i in range(self.nexp)], 1)
             # this is evaluating all experiments for all thetas, which is overkill
 
-    #@profile
     def llik(self, yobs, pred, cov):
         vec = yobs - pred 
         out = -.5*(cov['ldet'] + vec.T @ cov['inv'] @ vec)
         return out
-    #@profile
     def lik_cov_inv(self, s2vec):
         vec = self.trunc_error_var + s2vec
-        #mat = np.diag(vec) + self.basis @ np.diag(self.emu_vars) @ self.basis.T
-        #inv = np.linalg.inv(mat)
-        #ldet = np.linalg.slogdet(mat)[1]
-        #out = {'inv' : inv, 'ldet' : ldet}
         Ainv = np.diag(1/vec)
         Aldet = np.log(vec).sum()
         out = self.swm(Ainv, self.basis, np.diag(1/self.emu_vars), self.basis.T, Aldet, np.log(self.emu_vars).sum())
         return out
-    #@profile
     def chol_solve(self, x):
         mat = cho_factor(x)
         ldet = 2 * np.sum(np.log(np.diag(mat[0])))
-        ##la.dpotri(mat, overwrite_c=True) # overwrites mat with original matrix inverse, but not correct
         #inv = cho_solve(mat, np.eye(x.shape[0])) # correct, but slower for small dimension
         inv = np.linalg.inv(x)
         out = {'inv' : inv, 'ldet' : ldet}
         return out
-    #@profile
     def swm(self, Ainv, U, Cinv, V, Aldet, Cldet): # sherman woodbury morrison (A+UCV)^-1 and |A+UCV|
         in_mat = self.chol_solve(Cinv + V @ Ainv @ U)
         inv = Ainv - Ainv @ U @ in_mat['inv'] @ V @ Ainv
@@ -399,9 +374,6 @@ def swm(self, Ainv, U, Cinv, V, Aldet, Cldet): # sherman woodbury morrison (A+UC
         return out
 
 
-
-
-
 class ModelF(AbstractModel):
     def __init__(self, f, input_names, exp_ind=None, s2='gibbs'): # not sure if this is vectorized
         self.mod = f

impala/superCal/plots.py → impala/plots.py

@@ -15,6 +15,10 @@
 from numpy.random import uniform
 
 
+def autocorr(x):
+    result = np.correlate(x, x, mode='full')
+    return result[int(result.size/2):]
+
 # diagnostic plots to include:
 #    - traces of theta, s2, but for hier which ones? All?
 #   - tempering and MCMC counts
@@ -25,6 +29,44 @@ def __init__(self, setup, out):
         self.out   = out
         return
 
+    def plot_log_likehood(save_fn=None):
+        figure = plt.figure()
+        plt.plot(self.out.llik)
+        plt.title('Log Likelihood')
+        if save_fn is not None:
+            plt.savefig(save_fn + 'platinumnp_llike.png', dpi=300, transparent=True)  
+        return
+
+    def plot_chains(save_fn=None, input_names = None):
+        if input_names is None:
+            input_names = self.out.theta_native.keys()
+        nplots = len(input_names)
+        ns1 = int(np.ceil(np.sqrt(nplots)))
+        ns2 = int(np.round(np.sqrt(nplots)))
+        fig, axs = plt.subplots(ns1, ns2)
+        for i, ax in enumerate(fig.axes):
+            if i < nplots:
+                ax.plot(self.out.theta_native[input_names[i]])
+                ax.set(title=input_names[i])
+        plt.tight_layout()
+        if save_fn is not None:
+            plt.savefig(save_fn + 'platinumnp_chains.png', dpi=300, transparent=True)  
+
+        return
+
+    def plot_acf(save_fn=None, input_names = None):
+        if input_names is None:
+            input_names = self.out.theta_native.keys()
+        nplots = len(input_names)
+        fig, axs = plt.subplots(nplots, 1)
+        for i in range(nplots):
+            axs[i].plot(autocorr(self.out.theta_native[input_names[i]]))
+            axs[i].set(title='ACF '+ input_names[i])
+        plt.tight_layout()
+        if save_fn is not None:
+            plt.savefig(save_fn + 'platinumnp_acf.png', dpi=300, transparent=True)  
+        return
+
     def pooled_trace_plots():
         pass
 

pyproject.toml

@@ -0,0 +1,38 @@
+[project]
+name = "impala"
+version = "0.1"
+description = "Bayesian model calibration"
+authors = [
+    {name = "Devin Francom", email = "[email protected]"},
+    {name = "Arthur Lui", email = "[email protected]"},
+    {name = "Peter Trubey", email = "[email protected]"},
+	{name = "J. Derek Tucker", email = "[email protected]"}
+]
+license = {text = "BSD 3-Clause"}
+readme = "README.md"
+requires-python = ">=3.6"
+
+keywords = ["modelcalibration"]
+
+dependencies=[
+        "numpy",
+        "matplotlib",
+        "scipy",
+        "pyBASS",
+]
+
+classifiers = [
+    'License :: OSI Approved :: BSD License',
+    'Operating System :: OS Independent',
+    'Programming Language :: Python',
+    'Topic :: Scientific/Engineering',
+    'Topic :: Scientific/Engineering :: Mathematics',
+    'Programming Language :: Python :: 3',
+    'Programming Language :: Python :: 3.6',
+]
+
+[project.urls]
+repository = "http://www.github.com/lanl/impala"
+
+[tool.setuptools.packages.find]
+exclude = ["data*", "tests*"]

requirements.txt

@@ -0,0 +1,4 @@
+numpy
+pyBASS
+scipy
+matplotlib

tests/img/.gitignore

@@ -0,0 +1,2 @@
+*
+!.gitignore