diff --git a/.coverage.DESKTOP-ATMEKSV.10780.XNWpwYfx b/.coverage.DESKTOP-ATMEKSV.10780.XNWpwYfx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.12404.XzQmGlCx b/.coverage.DESKTOP-ATMEKSV.12404.XzQmGlCx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.12452.XgsSHyzx b/.coverage.DESKTOP-ATMEKSV.12452.XgsSHyzx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.21288.XCZBWtjx b/.coverage.DESKTOP-ATMEKSV.21288.XCZBWtjx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.21748.XBOljdvx b/.coverage.DESKTOP-ATMEKSV.21748.XBOljdvx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.24064.XFYCqtPx b/.coverage.DESKTOP-ATMEKSV.24064.XFYCqtPx
deleted file mode 100644
index e1f7e6f75f5aa9a7a33884193b67b7cc22082200..0000000000000000000000000000000000000000
Binary files a/.coverage.DESKTOP-ATMEKSV.24064.XFYCqtPx and /dev/null differ
diff --git a/.coverage.DESKTOP-ATMEKSV.26128.XrTwIHQx b/.coverage.DESKTOP-ATMEKSV.26128.XrTwIHQx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.26336.XXhjOUCx b/.coverage.DESKTOP-ATMEKSV.26336.XXhjOUCx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.27368.XKadoIex b/.coverage.DESKTOP-ATMEKSV.27368.XKadoIex
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.27856.XpjNuPSx b/.coverage.DESKTOP-ATMEKSV.27856.XpjNuPSx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.28136.XQytruKx b/.coverage.DESKTOP-ATMEKSV.28136.XQytruKx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.30088.XBTNpNAx b/.coverage.DESKTOP-ATMEKSV.30088.XBTNpNAx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.31644.XFWeLHFx b/.coverage.DESKTOP-ATMEKSV.31644.XFWeLHFx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.32652.XxCzZPCx b/.coverage.DESKTOP-ATMEKSV.32652.XxCzZPCx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.32976.XFBoIPNx b/.coverage.DESKTOP-ATMEKSV.32976.XFBoIPNx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.33768.XqlATslx b/.coverage.DESKTOP-ATMEKSV.33768.XqlATslx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.33992.XFaaFKJx b/.coverage.DESKTOP-ATMEKSV.33992.XFaaFKJx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.35124.XBfpxDZx b/.coverage.DESKTOP-ATMEKSV.35124.XBfpxDZx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.35624.XJVIwqYx b/.coverage.DESKTOP-ATMEKSV.35624.XJVIwqYx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/Outputs_Bayes_None_Calib/emcee_sampler.h5 b/Outputs_Bayes_None_Calib/emcee_sampler.h5
new file mode 100644
index 0000000000000000000000000000000000000000..bd7c2effd1b69997af86c87a8380da8ba2b9ad3d
Binary files /dev/null and b/Outputs_Bayes_None_Calib/emcee_sampler.h5 differ
diff --git a/docs/diagrams/.$Structure_BayesInf.drawio.dtmp b/docs/diagrams/.$Structure_BayesInf.drawio.dtmp
new file mode 100644
index 0000000000000000000000000000000000000000..bb32f9077e4d990a49743e90a0e4ea0837c77a5d
--- /dev/null
+++ b/docs/diagrams/.$Structure_BayesInf.drawio.dtmp
@@ -0,0 +1,449 @@
+<mxfile host="Electron" modified="2024-04-10T14:28:07.605Z" agent="Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) draw.io/22.1.11 Chrome/114.0.5735.289 Electron/25.9.8 Safari/537.36" etag="Yia9ahrH6Vnw8ZCcm-TH" version="22.1.11" type="device" pages="2">
+ <diagram name="Class and function structure" id="efOe0Jku58RX-i1bv-3b">
+ <mxGraphModel dx="1686" dy="1114" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="850" pageHeight="1100" math="0" shadow="0">
+ <root>
+ <mxCell id="0" />
+ <mxCell id="1" parent="0" />
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-1" value="_kernel_rbf" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="1070" y="200" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-2" value="_logpdf" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="860" y="130" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-10" value="<p style="margin:0px;margin-top:4px;text-align:center;"><b>BayesInf</b></p><hr size="1"><div style="height:2px;"></div>" style="verticalAlign=top;align=left;overflow=fill;fontSize=12;fontFamily=Helvetica;html=1;whiteSpace=wrap;" vertex="1" parent="1">
+ <mxGeometry x="-80" y="280" width="1270" height="820" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-24" value="if self.bootstrap <br>or self.bayes_loocv <br>or self.just_analysis" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0;exitY=0.5;exitDx=0;exitDy=0;entryX=1;entryY=0.5;entryDx=0;entryDy=0;labelBackgroundColor=#ffae00;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-13">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-31" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-18">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-42" value="if self.name != 'valid'<br>and self.inference_method != 'rejection'" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=default;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-31">
+ <mxGeometry x="0.5646" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-32" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0;entryY=0.5;entryDx=0;entryDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-22">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-43" value="if self.inference_method == 'mcmc'" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-32">
+ <mxGeometry x="-0.0958" y="-1" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-33" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;entryX=0;entryY=0.5;entryDx=0;entryDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-19">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-52" value="always" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=#C2C2C2;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-33">
+ <mxGeometry x="-0.112" y="1" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-34" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-21">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-47" value="if self.plot_post_pred" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-34">
+ <mxGeometry x="0.2399" y="-1" relative="1" as="geometry">
+ <mxPoint y="1" as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-35" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-20">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-46" value="if self.plot_map_pred" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-35">
+ <mxGeometry x="0.4183" y="-1" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-54" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-55" value="if self.bootstrap" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=#FF9A03;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-54">
+ <mxGeometry x="0.1816" y="3" relative="1" as="geometry">
+ <mxPoint x="1" as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-57" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;entryX=0;entryY=0.5;entryDx=0;entryDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-56">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-58" value="always" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=#FF9A03;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-57">
+ <mxGeometry x="0.7182" y="2" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-60" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.25;exitY=1;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-59">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-61" value="if self.error_model<br>and self.name == 'calib'" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-60">
+ <mxGeometry x="0.3024" y="2" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-9" value="create_inference" style="html=1;whiteSpace=wrap;strokeWidth=2;" vertex="1" parent="1">
+ <mxGeometry x="405" y="539" width="110" height="50" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-25" value="if len(self.perturbed_data) == 0" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-13" target="xary-zVek9Bg-A1b1ZmA-14">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-27" value="if not self.emulator" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-13" target="xary-zVek9Bg-A1b1ZmA-15">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-29" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;entryX=0;entryY=0.5;entryDx=0;entryDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-13" target="xary-zVek9Bg-A1b1ZmA-16">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-44" value="always" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=#cdcbcb;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-29">
+ <mxGeometry x="0.4722" y="1" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-30" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;entryX=0;entryY=0.5;entryDx=0;entryDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-13" target="xary-zVek9Bg-A1b1ZmA-17">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-41" value="if self.emulator" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-30">
+ <mxGeometry x="0.6143" y="-3" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-62" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-13" target="xary-zVek9Bg-A1b1ZmA-59">
+ <mxGeometry relative="1" as="geometry">
+ <mxPoint x="340" y="680" as="targetPoint" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-63" value="if self.error_model<br>and self.name == 'valid'" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=default;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-62">
+ <mxGeometry x="-0.3906" relative="1" as="geometry">
+ <mxPoint y="156" as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-13" value="perform_bootstrap" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
+ <mxGeometry x="40" y="310" width="110" height="50" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-14" value="_perturb_data" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
+ <mxGeometry x="-40" y="560" width="110" height="50" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-15" value="_eval_model" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
+ <mxGeometry x="1050" y="660" width="110" height="50" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-38" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0.5;entryY=1;entryDx=0;entryDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-16" target="xary-zVek9Bg-A1b1ZmA-1">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-49" value="if hasattr bias_inputs&nbsp;<br>and not hasattr error_model" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=#ffae00;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-38">
+ <mxGeometry x="0.3126" y="-3" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-39" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-16" target="xary-zVek9Bg-A1b1ZmA-2">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-16" value="normpdf" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
+ <mxGeometry x="650" y="390" width="110" height="50" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-40" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-17" target="xary-zVek9Bg-A1b1ZmA-2">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-50" value="always" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=#cdcbcb;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-40">
+ <mxGeometry x="-0.6073" y="-5" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-17" value="_corr_factor_BME" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
+ <mxGeometry x="650" y="450" width="110" height="50" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-18" value="_rejection_sampling" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
+ <mxGeometry x="400" y="990" width="120" height="50" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-26" value="if not self.emulator&nbsp;<br>and not self.inference_method == 'rejection'&nbsp;<br>and self.name == 'calib" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;entryX=0;entryY=0.5;entryDx=0;entryDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-19" target="xary-zVek9Bg-A1b1ZmA-15">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-37" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-19" target="xary-zVek9Bg-A1b1ZmA-1">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-48" value="if sigma2_prior is not None<br>and if hasattr bias_inputs<br>and if not hasattr error_model" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=#ffae00;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-37">
+ <mxGeometry x="-0.5544" y="-1" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-19" value="_posterior_predictive" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
+ <mxGeometry x="675" y="590" width="130" height="50" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-28" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-20" target="xary-zVek9Bg-A1b1ZmA-15">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-45" value="always" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=#cdcbcb;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-28">
+ <mxGeometry x="0.0517" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-20" value="_plot_max_a_posteriori" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
+ <mxGeometry x="495" y="790" width="140" height="50" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-21" value="plot_post_predictive" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
+ <mxGeometry x="630" y="720" width="120" height="50" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-22" value="<p style="margin:0px;margin-top:4px;text-align:center;"><b>MCMC</b></p><hr size="1"><div style="height:2px;"></div>" style="verticalAlign=top;align=left;overflow=fill;fontSize=12;fontFamily=Helvetica;html=1;whiteSpace=wrap;" vertex="1" parent="1">
+ <mxGeometry x="1230" y="425" width="760" height="275" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-36" value="Note: Arrows indicate function calls, beginning calls the end" style="text;html=1;strokeColor=none;fillColor=none;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;" vertex="1" parent="1">
+ <mxGeometry x="10" y="10" width="190" height="30" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-51" value="Color meanings:<br><span style="white-space: pre;">	</span>red: wrong, change<br><span style="white-space: pre;">	</span>orange: seems off, look at again" style="text;html=1;strokeColor=none;fillColor=none;align=left;verticalAlign=middle;whiteSpace=wrap;rounded=0;" vertex="1" parent="1">
+ <mxGeometry x="20" y="70" width="220" height="30" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-53" value="plot_log_BME" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="220" y="980" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-56" value="plot_post_params" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="660" y="840" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-59" value="create_error_model" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="35" y="840" width="120" height="60" as="geometry" />
+ </mxCell>
+ </root>
+ </mxGraphModel>
+ </diagram>
+ <diagram id="ME5gyYpVqUByTnAIOcMV" name="Parameter and function interaction">
+ <mxGraphModel dx="2049" dy="1366" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="850" pageHeight="1100" math="0" shadow="0">
+ <root>
+ <mxCell id="0" />
+ <mxCell id="1" parent="0" />
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-33" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0.5;entryY=1;entryDx=0;entryDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-1" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-54" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0.5;entryY=1;entryDx=0;entryDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-1" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-61" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0.5;entryY=1;entryDx=0;entryDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-1" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-1" value="engine" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="160" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-3" value="Discrepancy" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="240" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-71" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-4" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-4" value="emulator" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="320" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-37" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-5" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-57" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-5" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-65" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-5" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-5" value="name" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="400" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-47" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0.5;entryY=1;entryDx=0;entryDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-6" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-6" value="bootstrap" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="480" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-7" value="req_outputs" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="560" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-79" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-8" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-8" value="selected_indices" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="640" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-35" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-9" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-55" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-9" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-67" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-9" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-9" value="prior_samples" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="720" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-36" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-11" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-68" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-11" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-11" value="n_prior_samples" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="800" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-38" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-12" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-80" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-12" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-12" value="measured_data" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="880" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-58" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-13" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-13" value="inference_method" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="960" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-14" value="mcmc_params" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1040" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-63" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-15" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-15" value="perturbed_data" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1120" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-45" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-16" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-77" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-16" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-16" value="bayes_loocv" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1200" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-64" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-17" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-17" value="n_bootstrap_itrs" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1280" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-18" value="bootstrap_noise" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1360" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-46" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-19" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-78" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-19" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-19" value="just_analysis" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1440" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-20" value="valid_metrics" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1520" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-52" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-21" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-21" value="plot_post_pred" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1600" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-51" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-22" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-22" value="plot_map_pred" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1680" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-23" value="max_a_posteriori" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1760" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-24" value="corner_title_fmt" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1840" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-34" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-25" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-25" value="out_dir" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1920" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-50" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-26" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-66" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-26" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-26" value="error_model" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2000" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-56" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-27" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-72" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-27" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-27" value="bias_inputs" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2080" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-41" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-28" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-28" value="measurement_error" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2160" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-44" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-29" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-81" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-29" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-29" value="sigma2s" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2240" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-30" value="log_likes" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2320" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-82" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-31" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-31" value="dtype" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2400" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-32" value="create_inference" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="400" y="20" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-40" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-39" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-39" value="n_tot_measurement" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2480" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-43" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-42" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-42" value="Discrepancy" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2560" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-49" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-48" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-59" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-48" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-48" value="posterior_df" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2640" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-53" value="create_error_model" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="560" y="20" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-60" value="perform_bootstrap" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="720" y="20" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-75" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-69" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-69" value="__mean_pce_prior_pred" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2720" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-76" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-70" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-70" value="_std_pce_prior_pred" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2800" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-74" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-73" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-73" value="__model_prior_pred" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2880" width="120" height="60" as="geometry" />
+ </mxCell>
+ </root>
+ </mxGraphModel>
+ </diagram>
+</mxfile>
diff --git a/docs/diagrams/Structure_BayesInf.drawio b/docs/diagrams/Structure_BayesInf.drawio
index 1651230d37cc9da9120a7f0ae60ff3450d8f225a..ddabd26efb05c7db4ee7a49855f1a3585a883d58 100644
--- a/docs/diagrams/Structure_BayesInf.drawio
+++ b/docs/diagrams/Structure_BayesInf.drawio
@@ -1,6 +1,6 @@
-<mxfile host="Electron" modified="2024-04-10T11:36:50.177Z" agent="Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) draw.io/22.1.11 Chrome/114.0.5735.289 Electron/25.9.8 Safari/537.36" etag="zH0uklPxig_FBrKRFuLP" version="22.1.11" type="device">
- <diagram name="Page-1" id="efOe0Jku58RX-i1bv-3b">
- <mxGraphModel dx="836" dy="1114" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="850" pageHeight="1100" math="0" shadow="0">
+<mxfile host="Electron" modified="2024-04-10T14:28:03.443Z" agent="Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) draw.io/22.1.11 Chrome/114.0.5735.289 Electron/25.9.8 Safari/537.36" etag="3hqTB6OXl10kbEaaeQUl" version="22.1.11" type="device" pages="2">
+ <diagram name="Class and function structure" id="efOe0Jku58RX-i1bv-3b">
+ <mxGraphModel dx="1686" dy="1114" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="850" pageHeight="1100" math="0" shadow="0">
<root>
<mxCell id="0" />
<mxCell id="1" parent="0" />
@@ -11,7 +11,7 @@
<mxGeometry x="860" y="130" width="120" height="60" as="geometry" />
</mxCell>
<mxCell id="xary-zVek9Bg-A1b1ZmA-10" value="<p style="margin:0px;margin-top:4px;text-align:center;"><b>BayesInf</b></p><hr size="1"><div style="height:2px;"></div>" style="verticalAlign=top;align=left;overflow=fill;fontSize=12;fontFamily=Helvetica;html=1;whiteSpace=wrap;" vertex="1" parent="1">
- <mxGeometry x="40" y="280" width="1150" height="620" as="geometry" />
+ <mxGeometry x="-80" y="280" width="1270" height="820" as="geometry" />
</mxCell>
<mxCell id="xary-zVek9Bg-A1b1ZmA-24" value="if self.bootstrap <br>or self.bayes_loocv <br>or self.just_analysis" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0;exitY=0.5;exitDx=0;exitDy=0;entryX=1;entryY=0.5;entryDx=0;entryDy=0;labelBackgroundColor=#ffae00;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-13">
<mxGeometry relative="1" as="geometry" />
@@ -56,6 +56,30 @@
<mxPoint as="offset" />
</mxGeometry>
</mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-54" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-55" value="if self.bootstrap" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=#FF9A03;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-54">
+ <mxGeometry x="0.1816" y="3" relative="1" as="geometry">
+ <mxPoint x="1" as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-57" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;entryX=0;entryY=0.5;entryDx=0;entryDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-56">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-58" value="always" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=#FF9A03;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-57">
+ <mxGeometry x="0.7182" y="2" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-60" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.25;exitY=1;exitDx=0;exitDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-9" target="xary-zVek9Bg-A1b1ZmA-59">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-61" value="if self.error_model<br>and self.name == 'calib'" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-60">
+ <mxGeometry x="0.3024" y="2" relative="1" as="geometry">
+ <mxPoint as="offset" />
+ </mxGeometry>
+ </mxCell>
<mxCell id="xary-zVek9Bg-A1b1ZmA-9" value="create_inference" style="html=1;whiteSpace=wrap;strokeWidth=2;" vertex="1" parent="1">
<mxGeometry x="405" y="539" width="110" height="50" as="geometry" />
</mxCell>
@@ -81,11 +105,21 @@
<mxPoint as="offset" />
</mxGeometry>
</mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-62" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-13" target="xary-zVek9Bg-A1b1ZmA-59">
+ <mxGeometry relative="1" as="geometry">
+ <mxPoint x="340" y="680" as="targetPoint" />
+ </mxGeometry>
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-63" value="if self.error_model<br>and self.name == 'valid'" style="edgeLabel;html=1;align=center;verticalAlign=middle;resizable=0;points=[];labelBackgroundColor=default;" vertex="1" connectable="0" parent="xary-zVek9Bg-A1b1ZmA-62">
+ <mxGeometry x="-0.3906" relative="1" as="geometry">
+ <mxPoint y="156" as="offset" />
+ </mxGeometry>
+ </mxCell>
<mxCell id="xary-zVek9Bg-A1b1ZmA-13" value="perform_bootstrap" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
- <mxGeometry x="150" y="310" width="110" height="50" as="geometry" />
+ <mxGeometry x="40" y="310" width="110" height="50" as="geometry" />
</mxCell>
<mxCell id="xary-zVek9Bg-A1b1ZmA-14" value="_perturb_data" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
- <mxGeometry x="150" y="760" width="110" height="50" as="geometry" />
+ <mxGeometry x="-40" y="560" width="110" height="50" as="geometry" />
</mxCell>
<mxCell id="xary-zVek9Bg-A1b1ZmA-15" value="_eval_model" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
<mxGeometry x="1050" y="660" width="110" height="50" as="geometry" />
@@ -116,7 +150,7 @@
<mxGeometry x="650" y="450" width="110" height="50" as="geometry" />
</mxCell>
<mxCell id="xary-zVek9Bg-A1b1ZmA-18" value="_rejection_sampling" style="html=1;whiteSpace=wrap;" vertex="1" parent="1">
- <mxGeometry x="330" y="790" width="120" height="50" as="geometry" />
+ <mxGeometry x="400" y="990" width="120" height="50" as="geometry" />
</mxCell>
<mxCell id="xary-zVek9Bg-A1b1ZmA-26" value="if not self.emulator&nbsp;<br>and not self.inference_method == 'rejection'&nbsp;<br>and self.name == 'calib" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.5;exitY=1;exitDx=0;exitDy=0;entryX=0;entryY=0.5;entryDx=0;entryDy=0;" edge="1" parent="1" source="xary-zVek9Bg-A1b1ZmA-19" target="xary-zVek9Bg-A1b1ZmA-15">
<mxGeometry relative="1" as="geometry" />
@@ -155,6 +189,260 @@
<mxCell id="xary-zVek9Bg-A1b1ZmA-51" value="Color meanings:<br><span style="white-space: pre;">	</span>red: wrong, change<br><span style="white-space: pre;">	</span>orange: seems off, look at again" style="text;html=1;strokeColor=none;fillColor=none;align=left;verticalAlign=middle;whiteSpace=wrap;rounded=0;" vertex="1" parent="1">
<mxGeometry x="20" y="70" width="220" height="30" as="geometry" />
</mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-53" value="plot_log_BME" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="220" y="980" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-56" value="plot_post_params" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="660" y="840" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="xary-zVek9Bg-A1b1ZmA-59" value="create_error_model" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="35" y="840" width="120" height="60" as="geometry" />
+ </mxCell>
+ </root>
+ </mxGraphModel>
+ </diagram>
+ <diagram id="ME5gyYpVqUByTnAIOcMV" name="Parameter and function interaction">
+ <mxGraphModel dx="2049" dy="1366" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="850" pageHeight="1100" math="0" shadow="0">
+ <root>
+ <mxCell id="0" />
+ <mxCell id="1" parent="0" />
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-33" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0.5;entryY=1;entryDx=0;entryDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-1" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-54" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0.5;entryY=1;entryDx=0;entryDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-1" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-61" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0.5;entryY=1;entryDx=0;entryDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-1" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-1" value="engine" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="160" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-3" value="Discrepancy" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="240" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-71" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-4" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-4" value="emulator" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="320" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-37" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-5" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-57" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-5" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-65" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-5" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-5" value="name" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="400" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-47" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0.5;entryY=1;entryDx=0;entryDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-6" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-6" value="bootstrap" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="480" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-7" value="req_outputs" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="560" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-79" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-8" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-8" value="selected_indices" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="640" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-35" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-9" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-55" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-9" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-67" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-9" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-9" value="prior_samples" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="720" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-36" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-11" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-68" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-11" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-11" value="n_prior_samples" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="800" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-38" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-12" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-80" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-12" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-12" value="measured_data" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="880" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-58" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-13" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-13" value="inference_method" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="960" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-14" value="mcmc_params" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1040" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-63" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-15" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-15" value="perturbed_data" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1120" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-45" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-16" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-77" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-16" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-16" value="bayes_loocv" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1200" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-64" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-17" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-17" value="n_bootstrap_itrs" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1280" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-18" value="bootstrap_noise" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1360" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-46" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-19" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-78" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-19" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-19" value="just_analysis" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1440" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-20" value="valid_metrics" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1520" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-52" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-21" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-21" value="plot_post_pred" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1600" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-51" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-22" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-22" value="plot_map_pred" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1680" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-23" value="max_a_posteriori" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1760" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-24" value="corner_title_fmt" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1840" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-34" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-25" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-25" value="out_dir" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="1920" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-50" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-26" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-66" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-26" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-26" value="error_model" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2000" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-56" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-27" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-72" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-27" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-27" value="bias_inputs" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2080" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-41" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-28" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-28" value="measurement_error" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2160" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-44" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-29" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-81" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-29" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-29" value="sigma2s" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2240" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-30" value="log_likes" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2320" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-82" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-31" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-31" value="dtype" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2400" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-32" value="create_inference" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="400" y="20" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-40" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-39" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-39" value="n_tot_measurement" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2480" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-43" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-42" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-42" value="Discrepancy" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2560" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-49" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-48" target="K5oJ7VEt7dPmeK6pba1f-32">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-59" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-48" target="K5oJ7VEt7dPmeK6pba1f-53">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-48" value="posterior_df" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2640" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-53" value="create_error_model" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="560" y="20" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-60" value="perform_bootstrap" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="720" y="20" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-75" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-69" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-69" value="__mean_pce_prior_pred" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2720" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-76" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-70" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-70" value="_std_pce_prior_pred" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2800" width="120" height="60" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-74" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;" edge="1" parent="1" source="K5oJ7VEt7dPmeK6pba1f-73" target="K5oJ7VEt7dPmeK6pba1f-60">
+ <mxGeometry relative="1" as="geometry" />
+ </mxCell>
+ <mxCell id="K5oJ7VEt7dPmeK6pba1f-73" value="__model_prior_pred" style="rounded=0;whiteSpace=wrap;html=1;" vertex="1" parent="1">
+ <mxGeometry x="40" y="2880" width="120" height="60" as="geometry" />
+ </mxCell>
</root>
</mxGraphModel>
</diagram>
diff --git a/.coverage b/examples/.coverage
similarity index 81%
rename from .coverage
rename to examples/.coverage
index e61eed4b7c43a3f551e07789589a565c8682ff7d..254e10e4371d703eefec0f0437f9c0575be3f5ec 100644
Binary files a/.coverage and b/examples/.coverage differ
diff --git a/examples/analytical-function/example_analytical_function.py b/examples/analytical-function/example_analytical_function.py
index 9ec23f8cdaae003c33bac4cfe75e9a1a4d82bc42..37900127ab39018eef4af9dc5fa3f0326e2f4eeb 100644
--- a/examples/analytical-function/example_analytical_function.py
+++ b/examples/analytical-function/example_analytical_function.py
@@ -227,7 +227,8 @@ if __name__ == "__main__":
MetaModelOpts.ExpDesign = ExpDesign
engine = Engine(MetaModelOpts, Model, ExpDesign)
engine.start_engine()
- engine.train_sequential()
+ #engine.train_sequential()
+ engine.train_normal()
# Load the objects
# with open(f"PCEModel_{Model.name}.pkl", "rb") as input:
@@ -266,9 +267,9 @@ if __name__ == "__main__":
# BayesOpts.selected_indices = [0, 3, 5, 7, 9]
# BME Bootstrap
- # BayesOpts.bootstrap = True
- # BayesOpts.n_bootstrap_itrs = 500
- # BayesOpts.bootstrap_noise = 100
+ BayesOpts.bootstrap = True
+ BayesOpts.n_bootstrap_itrs = 500
+ BayesOpts.bootstrap_noise = 100
# Bayesian cross validation
BayesOpts.bayes_loocv = True # TODO: test what this does
@@ -296,31 +297,34 @@ if __name__ == "__main__":
BayesOpts.Discrepancy = DiscrepancyOpts
# -- (Option C) --
- # DiscOutputOpts = Input()
- # # # OutputName = 'Z'
- # DiscOutputOpts.add_marginals()
- # DiscOutputOpts.Marginals[0].Nnme = '$\sigma^2_{\epsilon}$'
- # DiscOutputOpts.Marginals[0].dist_type = 'uniform'
- # DiscOutputOpts.Marginals[0].parameters = [0, 10]
- # BayesOpts.Discrepancy = {'known': DiscrepancyOpts,
- # 'infer': Discrepancy(DiscOutputOpts)}
-
- # BayesOpts.bias_inputs = {'Z':np.arange(0, 10, 1.).reshape(-1,1) / 9}
- # DiscOutputOpts = Input()
- # # OutputName = 'lambda'
- # DiscOutputOpts.add_marginals()
- # DiscOutputOpts.Marginals[0].name = '$\lambda$'
- # DiscOutputOpts.Marginals[0].dist_type = 'uniform'
- # DiscOutputOpts.Marginals[0].parameters = [0, 1]
-
- # # OutputName = 'sigma_f'
- # DiscOutputOpts.add_marginals()
- # DiscOutputOpts.Marginals[1].Name = '$\sigma_f$'
- # DiscOutputOpts.Marginals[1].dist_type = 'uniform'
- # DiscOutputOpts.Marginals[1].parameters = [0, 1e-4]
- # BayesOpts.Discrepancy = Discrepancy(DiscOutputOpts)
- # BayesOpts.Discrepancy = {'known': DiscrepancyOpts,
- # 'infer': Discrepancy(DiscOutputOpts)}
+ if 0:
+ DiscOutputOpts = Input()
+ # # # OutputName = 'Z'
+ DiscOutputOpts.add_marginals()
+ DiscOutputOpts.Marginals[0].Nnme = '$\sigma^2_{\epsilon}$'
+ DiscOutputOpts.Marginals[0].dist_type = 'uniform'
+ DiscOutputOpts.Marginals[0].parameters = [0, 10]
+ #BayesOpts.Discrepancy = {'known': DiscrepancyOpts,
+ # 'infer': Discrepancy(DiscOutputOpts)}
+
+ BayesOpts.bias_inputs = {'Z':np.arange(0, 10, 1.).reshape(-1,1) / 9}
+
+ DiscOutputOpts = Input()
+ # OutputName = 'lambda'
+ DiscOutputOpts.add_marginals()
+ DiscOutputOpts.Marginals[0].name = '$\lambda$'
+ DiscOutputOpts.Marginals[0].dist_type = 'uniform'
+ DiscOutputOpts.Marginals[0].parameters = [0, 1]
+
+ # # OutputName = 'sigma_f'
+ DiscOutputOpts.add_marginals()
+ DiscOutputOpts.Marginals[1].Name = '$\sigma_f$'
+ DiscOutputOpts.Marginals[1].dist_type = 'uniform'
+ 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()
diff --git a/src/bayesvalidrox/bayes_inference/bayes_inference.py b/src/bayesvalidrox/bayes_inference/bayes_inference.py
index 44e90b77c69b4c2c7bd27639727c5d2bd4fe8720..888a689f82b4840a37459423d89fd5c5182e2293 100644
--- a/src/bayesvalidrox/bayes_inference/bayes_inference.py
+++ b/src/bayesvalidrox/bayes_inference/bayes_inference.py
@@ -221,7 +221,6 @@ class BayesInference:
corner_title_fmt='.2e', out_dir = ''):
self.engine = engine
- self.MetaModel = engine.MetaModel
self.Discrepancy = discrepancy
self.emulator = emulator
self.name = name
@@ -244,32 +243,43 @@ class BayesInference:
self.max_a_posteriori = max_a_posteriori
self.corner_title_fmt = corner_title_fmt
self.out_dir = out_dir
+
+ # Other properties and parameters (found in code, but never set)
+ self.error_model = False # TODO: no example or use case for this!
+ self.bias_inputs = None
+ self.measurement_error = None # TODO: what is this?
+ self.sigma2s = None
+ self.log_likes = None
+ self.n_tot_measurement = None
+ self.Discrepancy = None
+ self.posterior_df = None
+ self.error_MetaModel = None
+ self._mean_pce_prior_pred = None
+ self._std_pce_prior_pred = None
+ self.__model_prior_pred = None
+
+ # System settings
+ if os.name == 'nt':
+ print('')
+ print('WARNING: Performing the inference on windows can lead to reduced accuracy!')
+ print('')
+ self.dtype=np.longdouble
+ else:
+ self.dtype=np.float128
- # -------------------------------------------------------------------------
- def create_inference(self):
+ def setup_inference(self):
"""
- Starts the inference.
-
- Returns
- -------
- BayesInference : obj
- The Bayes inference object.
-
- # TODO: should this function really return the class?
-
+ This function sets up the inference by checking the inputs and getting
+ needed data.
"""
- # Set some variables
- MetaModel = self.MetaModel
Model = self.engine.Model
- n_params = MetaModel.n_params
- output_names = Model.Output.names
- par_names = self.engine.ExpDesign.par_names
# Create output directory
if self.out_dir == '':
- self.out_dir = f'Outputs_Bayes_{Model.name}_{self.name}'
+ self.out_dir = f'Outputs_Bayes_{self.engine.Model.name}_{self.name}'
os.makedirs(self.out_dir, exist_ok=True)
+
# If the prior is set by the user, take it, else generate from ExpDes
if self.prior_samples is None:
self.prior_samples = self.engine.ExpDesign.generate_samples(
@@ -287,7 +297,8 @@ class BayesInference:
self.n_prior_samples = self.prior_samples.shape[0]
# ---------- Preparation of observation data ----------
- # Read observation data and perturb it if requested.
+ # Read observation data
+ # TODO: later use valid #of measurements. but here only get the model observations?
if self.measured_data is None:
self.measured_data = Model.read_observation(case=self.name)
# Convert measured_data to a data frame
@@ -297,11 +308,13 @@ class BayesInference:
# Extract the total number of measurement points
if self.name.lower() == 'calib':
self.n_tot_measurement = Model.n_obs
- else:
+ elif self.name.lower() == 'valid':
self.n_tot_measurement = Model.n_obs_valid
+ else:
+ raise AttributeError('The set inference type is not known! Use either `calib` or `valid`')
# Find measurement error (if not given) for post predictive plot
- if not hasattr(self, 'measurement_error'):
+ if self.measurement_error is None:
if isinstance(self.Discrepancy, dict):
Disc = self.Discrepancy['known']
else:
@@ -314,62 +327,95 @@ class BayesInference:
self.measurement_error = np.sqrt(Disc.parameters)
except TypeError:
pass
+ # TODO: need a transformation for given measurement error?
+
+ # Get Discrepancy type
+ opt_sigma_flag = isinstance(self.Discrepancy, dict)
+ opt_sigma = None
+ # Option A: known error with unknown bias term
+ if opt_sigma_flag and opt_sigma is None:
+ opt_sigma = 'A'
+ # Option B: The sigma2 is known (no bias term)
+ elif self.Discrepancy.parameters is not None:
+ opt_sigma = 'B'
+ # Option C: The sigma2 is unknown (bias term including error)
+ elif not isinstance(self.Discrepancy.InputDisc, str):
+ opt_sigma = 'C'
+ self.Discrepancy.opt_sigma = opt_sigma
+
+ # Set MCMC params if used
+ if self.inference_method.lower() == 'mcmc':
+ if self.mcmc_params is None:
+ self.mcmc_params = {}
+ par_list = ['init_samples', 'n_walkers', 'n_burn', 'n_steps',
+ 'moves', 'multiprocessing', 'verbose']
+ init_val = [None, 100, 200, 100000, None, False, False]
+ for i in range(len(par_list)):
+ if par_list[i] not in list(self.mcmc_params.keys()):
+ self.mcmc_params[par_list[i]] = init_val[i]
+
+
+ # -------------------------------------------------------------------------
+ def create_inference(self):
+ """
+ Starts the inference.
+ Returns
+ -------
+ BayesInference : obj
+ The Bayes inference object.
+
+ # TODO: should this function really return the class?
+
+ """
+ # Do general set up and check some parameters
+ self.setup_inference()
+
# ---------- Preparation of variance for covariance matrix ----------
- # Independent and identically distributed
+ # Independent and identically distributed # TODO: ??
total_sigma2 = dict()
- opt_sigma_flag = isinstance(self.Discrepancy, dict)
- opt_sigma = None
- for key_idx, key in enumerate(output_names):
+ opt_sigma = self.Discrepancy.opt_sigma
+ for key_idx, key in enumerate(self.engine.Model.Output.names):
# Find opt_sigma
- if opt_sigma_flag and opt_sigma is None:
- # Option A: known error with unknown bias term
- opt_sigma = 'A'
- known_discrepancy = self.Discrepancy['known'] # TODO: the syntax here looks different from expected
- self.Discrepancy = self.Discrepancy['infer'] # TODO: the syntax here looks different from expected
+ if opt_sigma == 'A':
+ known_discrepancy = self.Discrepancy['known']
+ self.Discrepancy = self.Discrepancy['infer']
sigma2 = np.array(known_discrepancy.parameters[key])
- elif self.Discrepancy.parameters is not None:
- # Option B: The sigma2 is known (no bias term)
- opt_sigma = 'B'
+ elif opt_sigma == 'B':
sigma2 = np.array(self.Discrepancy.parameters[key])
- elif not isinstance(self.Discrepancy.InputDisc, str):
- # Option C: The sigma2 is unknown (bias term including error)
- opt_sigma = 'C'
+ elif opt_sigma == 'C':
n_measurement = self.measured_data[key].values.shape
sigma2 = np.zeros((n_measurement[0]))
-
total_sigma2[key] = sigma2
- self.Discrepancy.opt_sigma = opt_sigma
self.Discrepancy.total_sigma2 = total_sigma2
# If inferred sigma2s obtained from e.g. calibration are given
try:
self.sigma2s = self.Discrepancy.get_sample(self.n_prior_samples)
except:
- pass #TODO: should an error be raised in this case?
+ pass #TODO: should an error be raised in this case? Should this at least be checked against opt_sigma?
# ---------------- Bootstrap & TOM --------------------
-
if self.bootstrap or self.bayes_loocv or self.just_analysis:
- self.perform_bootstrap(opt_sigma, total_sigma2)
+ self.perform_bootstrap(total_sigma2)
else:
print('No bootstrap for TOM performed!') # TODO: stop the code? Use n_bootstrap = 1?
# ---------------- Parameter Bayesian inference ----------------
+ # Convert to a dataframe if samples are provided after calibration.
if self.name.lower() == 'valid':
- # Convert to a dataframe if samples are provided after calibration.
- self.posterior_df = pd.DataFrame(self.prior_samples, columns=par_names)
+ self.posterior_df = pd.DataFrame(self.prior_samples, columns=self.engine.ExpDesign.par_names)
+ # Instantiate the MCMC object
elif self.inference_method.lower() == 'mcmc':
- # Instantiate the MCMC object
MCMC_Obj = MCMC(self)
self.posterior_df = MCMC_Obj.run_sampler(
self.measured_data, total_sigma2
)
+ # Rejection sampling
elif self.inference_method.lower() == 'rejection':
- # Rejection sampling
self.posterior_df = self._rejection_sampling()
else:
raise AttributeError('The chosen inference method is not available!')
@@ -383,100 +429,93 @@ class BayesInference:
print('-'*50)
# -------- Model Discrepancy -----------
- if hasattr(self, 'error_model') and self.error_model \
- and self.name.lower() == 'calib':
+ if self.error_model and self.name.lower() == 'calib': # TODO: where is this used and what does it actually do there?
+ self.create_error_model(opt_sigma = opt_sigma,
+ type_='posterior', sampler = MCMC_Obj)
+
+ # -------- Posterior predictive -----------
+ self._posterior_predictive()
+
+ # ------------------ Visualization --------------------
+ # Posterior parameters
+ self.plot_post_params(opt_sigma)
+
+ # Plot MAP
+ if self.plot_map_pred:
+ self._plot_max_a_posteriori()
+
+ # Plot log_BME dist
+ if self.bootstrap:
+ self.plot_log_BME()
+
+ # Plot posterior predictive
+ if self.plot_post_pred:
+ self._plot_post_predictive()
+
+ return self
+
+ def create_error_model(self, type_ = 'posterior', opt_sigma = 'B', sampler = None):
+ """
+ Creates an error model in the engine.MetaModel based on input dist
+ samples of the chosen type
+
+ Parameters
+ ----------
+ opt_sigma : string, optional
+ Type of uncertainty description, only used if type_=='posterior'.
+ The default is 'B'
+ type_ : string
+ Type of parameter samples to use, either 'prior' or 'posterior'.
+ The default is 'posterior'.
+ sampler : MCMC, optional
+ Should be an MCMC object if type=='posterior' and MCMC is used in
+ the inference.In al other cases this parameter is not needed.
+
+ Returns
+ -------
+ None.
+
+ """
+ n_params = self.engine.MetaModel.n_params
+
+ # Get MAP estimate from prior samples
+ if type_ == 'prior':
+ # Select prior ? mean as MAP
+ MAP_theta = self.prior_samples.mean(axis=0).reshape((1, n_params))
+
+ # Evaluate the (meta-)model at the MAP
+ y_MAP, y_std_MAP = self.engine.MetaModel.eval_metamodel(samples=MAP_theta)
+
+ # Train a GPR meta-model using MAP
+ self.error_MetaModel = self.engine.MetaModel.create_model_error(
+ self.bias_inputs, y_MAP, self.measured_data, name=self.name
+ )
+
+ # Get MAP estimate from posterior samples
+ if type_ == 'posterior':
if self.inference_method.lower() == 'mcmc':
- self.error_MetaModel = MCMC_Obj.error_MetaModel
+ self.error_MetaModel = sampler.error_MetaModel
else:
# Select posterior mean as MAP
if opt_sigma == "B":
posterior_df = self.posterior_df.values
else:
- posterior_df = self.posterior_df.values[:, :-Model.n_outputs]
+ posterior_df = self.posterior_df.values[:, :-self.engine.Model.n_outputs]
# Select posterior mean as Maximum a posteriori
map_theta = posterior_df.mean(axis=0).reshape((1, n_params))
# map_theta = stats.mode(Posterior_df,axis=0)[0]
# Evaluate the (meta-)model at the MAP
- y_MAP, y_std_MAP = MetaModel.eval_metamodel(samples=map_theta)
+ y_MAP, y_std_MAP = self.engine.MetaModel.eval_metamodel(samples=map_theta)
# Train a GPR meta-model using MAP
- self.error_MetaModel = MetaModel.create_model_error(
- self.bias_inputs, y_MAP, Name=self.name
+ self.error_MetaModel = self.engine.MetaModel.create_model_error(
+ self.bias_inputs, y_MAP, self.measured_data, name=self.name
)
- # -------- Posterior perdictive -----------
- self._posterior_predictive()
-
- # -----------------------------------------------------
- # ------------------ Visualization --------------------
- # -----------------------------------------------------
- # Create Output directory, if it doesn't exist already.
-
- # -------- Posterior parameters --------
- if opt_sigma != "B":
- par_names.extend(
- [self.Discrepancy.InputDisc.Marginals[i].name for i
- in range(len(self.Discrepancy.InputDisc.Marginals))]
- )
- # Pot with corner
- figPosterior = corner.corner(self.posterior_df.to_numpy(),
- labels=par_names,
- quantiles=[0.15, 0.5, 0.85],
- show_titles=True,
- title_fmt=self.corner_title_fmt,
- labelpad=0.2,
- use_math_text=True,
- title_kwargs={"fontsize": 28},
- plot_datapoints=False,
- plot_density=False,
- fill_contours=True,
- smooth=0.5,
- smooth1d=0.5)
-
- # Loop over axes and set x limits
- if opt_sigma == "B":
- axes = np.array(figPosterior.axes).reshape(
- (len(par_names), len(par_names))
- )
- for yi in range(len(par_names)):
- ax = axes[yi, yi]
- ax.set_xlim(self.engine.ExpDesign.bound_tuples[yi])
- for xi in range(yi):
- ax = axes[yi, xi]
- ax.set_xlim(self.engine.ExpDesign.bound_tuples[xi])
- plt.close()
-
- # Turn off gridlines
- for ax in figPosterior.axes:
- ax.grid(False)
-
- if self.emulator:
- plotname = f'/Posterior_Dist_{Model.name}_emulator'
- else:
- plotname = f'/Posterior_Dist_{Model.name}'
-
- figPosterior.set_size_inches((24, 16))
- figPosterior.savefig(f'./{self.out_dir}{plotname}.pdf',
- bbox_inches='tight')
-
- # -------- Plot MAP --------
- if self.plot_map_pred:
- self._plot_max_a_posteriori()
-
- # -------- Plot log_BME dist --------
- if self.bootstrap:
- self.plot_log_BME()
-
- # -------- Posterior perdictives --------
- if self.plot_post_pred:
- # Plot the posterior predictive
- self._plot_post_predictive()
-
- return self
- def perform_bootstrap(self, opt_sigma, total_sigma2):
+ def perform_bootstrap(self, total_sigma2):
"""
Perform bootstrap to get TOM (??)
@@ -491,9 +530,9 @@ class BayesInference:
None.
"""
- MetaModel = self.MetaModel
- n_params = MetaModel.n_params
+ MetaModel = self.engine.MetaModel
output_names = self.engine.Model.Output.names
+ opt_sigma = self.Discrepancy.opt_sigma
# Adding some zero mean noise to the observation function
if len(self.perturbed_data) == 0:
@@ -504,19 +543,8 @@ class BayesInference:
self.n_bootstrap_itrs = len(self.perturbed_data)
# -------- Model Discrepancy -----------
- if hasattr(self, 'error_model') and self.error_model \
- and self.name.lower() != 'calib':# TODO: what should be set so that this is tested?
- # Select prior ? mean as MAP
- MAP_theta = self.prior_samples.mean(axis=0).reshape((1, n_params))
-
- # Evaluate the (meta-)model at the MAP
- y_MAP, y_std_MAP = MetaModel.eval_metamodel(samples=MAP_theta)
-
- # Train a GPR meta-model using MAP
- self.error_MetaModel = MetaModel.create_model_error(
- self.bias_inputs, y_MAP, Name=self.name
- )
-
+ if self.error_model and self.name.lower() == 'valid':# TODO: what should be set so that this is tested?
+ self.create_error_model(type_='prior')
# -----------------------------------------------------
# ----- Loop over the perturbed observation data ------
# -----------------------------------------------------
@@ -536,10 +564,9 @@ class BayesInference:
self._std_pce_prior_pred = y_std
# Correct the predictions with Model discrepancy
- if hasattr(self, 'error_model') and self.error_model: # TODO this does not check for calib?
+ if self.error_model: # TODO this does not check for calib?
y_hat_corr, y_std = self.error_MetaModel.eval_model_error(
- self.bias_inputs, self.__mean_pce_prior_pred
- )
+ self.bias_inputs, self.__mean_pce_prior_pred)
self.__mean_pce_prior_pred = y_hat_corr
self._std_pce_prior_pred = y_std
@@ -548,12 +575,13 @@ class BayesInference:
surrError = MetaModel.rmse
else:
surrError = None
-
+ model_evals = self.__mean_pce_prior_pred
+
+ # Evaluate the model
else:
- # Evaluate the original model
self.__model_prior_pred = self._eval_model(
- samples=self.prior_samples, key='PriorPred'
- )
+ samples=self.prior_samples, key='PriorPred')
+ model_evals = self.__model_prior_pred
surrError = None
# Start the likelihood-BME computations for the perturbed data
@@ -562,17 +590,9 @@ class BayesInference:
total=self.n_bootstrap_itrs,
desc="Bootstrapping the BME calculations", ascii=True
):
- print('')
- #print(itr_idx, data)
- #print(np.nonzero(data))
-
+
# ---------------- Likelihood calculation ----------------
- if self.emulator: # TODO: do this outside of the loop?
- model_evals = self.__mean_pce_prior_pred
- else:
- model_evals = self.__model_prior_pred
-
- # Leave one out # TODO: why is this loo? It just looks at perturbed data?
+ # Leave one out (see form of perturbed data)
if self.bayes_loocv or self.just_analysis:
# Consider only non-zero entries
self.selected_indices = np.nonzero(data)[0]
@@ -587,7 +607,7 @@ class BayesInference:
}
# Unknown sigma2
- if opt_sigma == 'C' or hasattr(self, 'sigma2s'):
+ if opt_sigma == 'C' or self.sigma2s is not None:
logLikelihoods[:, itr_idx] = self.normpdf(
model_evals, data_dict, total_sigma2,
sigma2=self.sigma2s, std=surrError
@@ -602,7 +622,7 @@ class BayesInference:
# BME (log)
log_BME[itr_idx] = np.log(
np.nanmean(np.exp(logLikelihoods[:, itr_idx],
- dtype=np.longdouble))#float128))
+ dtype=self.dtype))
)
# BME correction when using Emulator
@@ -656,18 +676,15 @@ class BayesInference:
self.log_BME = log_BME
# BMECorrFactor (log) (Size: 1,n_bootstrap_itr)
+ # BME = BME + BMECorrFactor
if self.emulator:
- self.log_BME_corr_factor = BME_Corr
- # BME = BME + BMECorrFactor
- self.log_BME += self.log_BME_corr_factor
+ self.log_BME += BME_Corr
if 'kld' in list(map(str.lower, self.valid_metrics)):
self.KLD = KLD
if 'inf_entropy' in list(map(str.lower, self.valid_metrics)):
self.inf_entropy = inf_entropy
-
-
# -------------------------------------------------------------------------
def _perturb_data(self, data, output_names):
@@ -859,8 +876,7 @@ class BayesInference:
covMatrix = np.diag(tot_sigma2s)
# Check the type error term
- if hasattr(self, 'bias_inputs') and \
- not hasattr(self, 'error_model'):
+ if self.bias_inputs != None and self.error_model == None:
# Infer a Bias model usig Gaussian Process Regression
bias_inputs = np.hstack(
(self.bias_inputs[out],
@@ -906,7 +922,7 @@ class BayesInference:
A dictionary with known values of the covariance diagonal entries,
a.k.a sigma^2.
logBME : ??
- ??
+ The log_BME obtained from the estimated likelihoods
Returns
-------
@@ -915,7 +931,7 @@ class BayesInference:
"""
# Extract the requested model outputs for likelihood calulation
- MetaModel = self.MetaModel
+ MetaModel = self.engine.MetaModel
samples = self.engine.ExpDesign.X
model_outputs = self.engine.ExpDesign.Y
n_samples = samples.shape[0]
@@ -997,7 +1013,7 @@ class BayesInference:
if self.prior_samples is None:
raise AttributeError('No prior samples available!')
- if not hasattr(self, 'log_likes'):
+ if self.log_likes is None:
raise AttributeError('No log-likelihoods available!')
# Get sigmas # TODO: is this data uncertainty?
@@ -1018,11 +1034,7 @@ class BayesInference:
index = 0
# Use longdouble on windows, float128 on linux
- if os.name != 'nt':
- likelihoods = np.exp(self.log_likes[:, index], dtype=np.float128)
- else:
- print('WARNING: Performing the inference on windows can lead to reduced accuracy!')
- likelihoods = np.exp(self.log_likes[:, index], dtype=np.longdouble)
+ likelihoods = np.exp(self.log_likes[:, index], dtype=self.dtype)
n_samples = len(likelihoods)
norm_likelihoods = likelihoods / np.max(likelihoods)
@@ -1057,7 +1069,7 @@ class BayesInference:
"""
- MetaModel = self.MetaModel
+ MetaModel = self.engine.MetaModel
Model = self.engine.Model
# Read observation data and perturb it if requested # TODO: where is the perturbation?
@@ -1067,11 +1079,10 @@ class BayesInference:
if not isinstance(self.measured_data, pd.DataFrame):
self.measured_data = pd.DataFrame(self.measured_data)
- # X_values
+ # X_values and prior sigma2
x_values = self.engine.ExpDesign.x_values
-
try:
- sigma2_prior = self.Discrepancy.sigma2_prior
+ sigma2_prior = self.Discrepancy.sigma2_prior # TODO: what is this? Looks to be built for a different Discrepancy structure
except:
sigma2_prior = None
@@ -1079,7 +1090,7 @@ class BayesInference:
posterior_df = self.posterior_df
# Take care of the sigma2
- if sigma2_prior is not None:
+ if sigma2_prior is not None: # TODO: why is this the if for this code?
try:
sigma2s = posterior_df[self.Discrepancy.name].values # TODO: what is Discrepancy.name?
posterior_df = posterior_df.drop(
@@ -1090,10 +1101,10 @@ class BayesInference:
# Posterior predictive
if self.emulator:
- if self.inference_method == 'rejection': # TODO: combine these two?
- prior_pred = self.__mean_pce_prior_pred
- if self.name.lower() != 'calib':
- post_pred = self.__mean_pce_prior_pred
+ if self.inference_method == 'rejection': # TODO: combine these two? Why is there no post_pred_std for rejection sampling?
+ prior_pred = self._mean_pce_prior_pred
+ if self.name.lower() == 'valid':
+ post_pred = self._mean_pce_prior_pred
post_pred_std = self._std_pce_prior_pred
else:
post_pred, post_pred_std = MetaModel.eval_metamodel( # TODO: recheck if this is needed
@@ -1103,7 +1114,7 @@ class BayesInference:
else: # TODO: see emulator version
if self.inference_method == 'rejection':
prior_pred = self.__model_prior_pred
- if self.name.lower() != 'calib':
+ if self.name.lower() == 'valid':
post_pred = self.__mean_pce_prior_pred,
post_pred_std = self._std_pce_prior_pred
else:
@@ -1111,13 +1122,13 @@ class BayesInference:
samples=posterior_df.values, key='PostPred'
)
# Correct the predictions with Model discrepancy
- if hasattr(self, 'error_model') and self.error_model:
+ if self.error_model:
y_hat, y_std = self.error_MetaModel.eval_model_error(
self.bias_inputs, post_pred
)
post_pred, post_pred_std = y_hat, y_std
- # Add discrepancy from likelihood Sample to the current posterior runs
+ # Add discrepancy from likelihood samples to the current posterior runs
total_sigma2 = self.Discrepancy.total_sigma2
post_pred_withnoise = copy.deepcopy(post_pred)
for varIdx, var in enumerate(Model.Output.names):
@@ -1129,11 +1140,10 @@ class BayesInference:
tot_sigma2 = clean_sigma2[:len(pred)]
cov = np.diag(tot_sigma2)
- # Check the type error term
+ # Account for additional error terms
if sigma2_prior is not None:
# Inferred sigma2s
- if hasattr(self, 'bias_inputs') and \
- not hasattr(self, 'error_model'):
+ if self.bias_inputs != None and self.error_model == None:
# TODO: Infer a Bias model usig GPR
bias_inputs = np.hstack((
self.bias_inputs[var], pred.reshape(-1, 1)))
@@ -1149,6 +1159,7 @@ class BayesInference:
# Convert biasSigma2s to a covMatrix
cov += sigma2 * np.eye(len(pred))
+ # Add predictive metamodel error/uncertainty
if self.emulator:
if hasattr(MetaModel, 'rmse') and \
MetaModel.rmse is not None:
@@ -1159,7 +1170,7 @@ class BayesInference:
cov += np.diag(stdPCE**2)
# Sample a multivariate normal distribution with mean of
- # prediction and variance of cov
+ # posterior prediction and variance of cov
post_pred_withnoise[var][i] = np.random.multivariate_normal(
pred, cov, 1
)
@@ -1242,10 +1253,9 @@ class BayesInference:
"""
- MetaModel = self.MetaModel
+ MetaModel = self.engine.MetaModel
Model = self.engine.Model
opt_sigma = self.Discrepancy.opt_sigma
-
# -------- Find MAP and run MetaModel and origModel --------
# Compute the MAP
if self.max_a_posteriori.lower() == 'mean':
@@ -1330,7 +1340,80 @@ class BayesInference:
pdf.close()
+ def plot_post_params(self, opt_sigma):
+ """
+ Plots the multivar. posterior parameter distribution.
+
+
+ Parameters
+ ----------
+ opt_sigma : string
+ Type of uncertainty description available.
+
+ Returns
+ -------
+ None.
+
+ """
+ par_names = self.engine.ExpDesign.par_names
+ if opt_sigma != "B":
+ par_names.extend(
+ [self.Discrepancy.InputDisc.Marginals[i].name for i
+ in range(len(self.Discrepancy.InputDisc.Marginals))]
+ )
+ # Pot with corner
+ figPosterior = corner.corner(self.posterior_df.to_numpy(),
+ labels=par_names,
+ quantiles=[0.15, 0.5, 0.85],
+ show_titles=True,
+ title_fmt=self.corner_title_fmt,
+ labelpad=0.2,
+ use_math_text=True,
+ title_kwargs={"fontsize": 28},
+ plot_datapoints=False,
+ plot_density=False,
+ fill_contours=True,
+ smooth=0.5,
+ smooth1d=0.5)
+
+ # Loop over axes and set x limits
+ if opt_sigma == "B":
+ axes = np.array(figPosterior.axes).reshape(
+ (len(par_names), len(par_names))
+ )
+ for yi in range(len(par_names)):
+ ax = axes[yi, yi]
+ ax.set_xlim(self.engine.ExpDesign.bound_tuples[yi])
+ for xi in range(yi):
+ ax = axes[yi, xi]
+ ax.set_xlim(self.engine.ExpDesign.bound_tuples[xi])
+ plt.close()
+
+ # Turn off gridlines
+ for ax in figPosterior.axes:
+ ax.grid(False)
+
+ if self.emulator:
+ plotname = f'/Posterior_Dist_{self.engine.Model.name}_emulator'
+ else:
+ plotname = f'/Posterior_Dist_{self.engine.Model.name}'
+
+ figPosterior.set_size_inches((24, 16))
+ figPosterior.savefig(f'./{self.out_dir}{plotname}.pdf',
+ bbox_inches='tight')
+
+ plt.clf()
+
+
def plot_log_BME(self):
+ """
+ Plots the log_BME if bootstrap is active.
+
+ Returns
+ -------
+ None.
+
+ """
# Computing the TOM performance
self.log_BME_tom = stats.chi2.rvs(
@@ -1353,13 +1436,14 @@ class BayesInference:
ax.legend(handles=legend_elements)
if self.emulator:
- plotname = f'/BME_hist_{self.Model.name}_emulator'
+ plotname = f'/BME_hist_{self.engine.Model.name}_emulator'
else:
- plotname = f'/BME_hist_{self.Model.name}'
+ plotname = f'/BME_hist_{self.engine.Model.name}'
- plt.savefig(f'./{self.self.out_dir}{plotname}.pdf', bbox_inches='tight')
+ plt.savefig(f'./{self.out_dir}{plotname}.pdf', bbox_inches='tight')
+
plt.show()
- plt.close()
+ plt.clf()
# -------------------------------------------------------------------------
@@ -1382,11 +1466,10 @@ class BayesInference:
# --- Read prior and posterior predictive ---
if self.inference_method == 'rejection' and \
- self.name.lower() != 'valid':
+ self.name.lower() == 'calib':
# --- Prior ---
# Load posterior predictive
- f = h5py.File(
- f'{self.out_dir}/priorPredictive.hdf5', 'r+')
+ f = h5py.File(f'{self.out_dir}/priorPredictive.hdf5', 'r+')
try:
x_coords = np.array(f[f"x_values/{out_name}"])
@@ -1516,3 +1599,5 @@ class BayesInference:
fig.savefig(f'./{self.out_dir}{plotname}_{out_name}.pdf',
bbox_inches='tight')
+
+ plt.clf()
diff --git a/src/bayesvalidrox/bayes_inference/mcmc.py b/src/bayesvalidrox/bayes_inference/mcmc.py
index d78d15b5fd90dc4477da7d0fd58da835acc75310..a672583bae6e1c273bf7d080062502f96bbf9199 100755
--- a/src/bayesvalidrox/bayes_inference/mcmc.py
+++ b/src/bayesvalidrox/bayes_inference/mcmc.py
@@ -15,6 +15,182 @@ import shutil
os.environ["OMP_NUM_THREADS"] = "1"
+# -------------------------------------------------------------------------
+def _check_ranges(theta, ranges): # TODO: this is a replica of exp_designs.check_ranges
+ """
+ This function checks if theta lies in the given ranges.
+
+ Parameters
+ ----------
+ theta : array
+ Proposed parameter set.
+ ranges : nested list
+ List of the praremeter ranges.
+
+ Returns
+ -------
+ c : bool
+ If it lies in the given range, it return True else False.
+
+ """
+ c = True
+ # traverse in the list1
+ for i, bounds in enumerate(ranges):
+ x = theta[i]
+ # condition check
+ if x < bounds[0] or x > bounds[1]:
+ c = False
+ return c
+ return c
+
+# -------------------------------------------------------------------------
+def gelman_rubin(chain, return_var=False):
+ """
+ The potential scale reduction factor (PSRF) defined by the variance
+ within one chain, W, with the variance between chains B.
+ Both variances are combined in a weighted sum to obtain an estimate of
+ the variance of a parameter \\( \\theta \\).The square root of the
+ ratio of this estimates variance to the within chain variance is called
+ the potential scale reduction.
+ For a well converged chain it should approach 1. Values greater than
+ 1.1 typically indicate that the chains have not yet fully converged.
+
+ Source: http://joergdietrich.github.io/emcee-convergence.html
+
+ https://github.com/jwalton3141/jwalton3141.github.io/blob/master/assets/posts/ESS/rwmh.py
+
+ Parameters
+ ----------
+ chain : array (n_walkers, n_steps, n_params)
+ The emcee ensamples.
+
+ Returns
+ -------
+ R_hat : float
+ The Gelman-Robin values.
+
+ """
+ chain = np.array(chain)
+ m_chains, n_iters = chain.shape[:2]
+
+ # Calculate between-chain variance
+ θb = np.mean(chain, axis=1)
+ θbb = np.mean(θb, axis=0)
+ B_over_n = ((θbb - θb)**2).sum(axis=0)
+ B_over_n /= (m_chains - 1)
+
+ # Calculate within-chain variances
+ ssq = np.var(chain, axis=1, ddof=1)
+ W = np.mean(ssq, axis=0)
+
+ # (over) estimate of variance
+ var_θ = W * (n_iters - 1) / n_iters + B_over_n
+
+ if return_var:
+ return var_θ
+ else:
+ # The square root of the ratio of this estimates variance to the
+ # within chain variance
+ R_hat = np.sqrt(var_θ / W)
+ return R_hat
+
+
+# -------------------------------------------------------------------------
+def _iterative_scheme(N1, N2, q11, q12, q21, q22, r0, neff, tol,
+ maxiter, criterion):
+ """
+ Iterative scheme as proposed in Meng and Wong (1996) to estimate the
+ marginal likelihood
+
+ """
+ l1 = q11 - q12
+ l2 = q21 - q22
+ # To increase numerical stability,
+ # subtracting the median of l1 from l1 & l2 later
+ lstar = np.median(l1)
+ s1 = neff/(neff + N2)
+ s2 = N2/(neff + N2)
+ r = r0
+ r_vals = [r]
+ logml = np.log(r) + lstar
+ criterion_val = 1 + tol
+
+ i = 0
+ while (i <= maxiter) & (criterion_val > tol):
+ rold = r
+ logmlold = logml
+ numi = np.exp(l2 - lstar)/(s1 * np.exp(l2 - lstar) + s2 * r)
+ deni = 1/(s1 * np.exp(l1 - lstar) + s2 * r)
+ if np.sum(~np.isfinite(numi))+np.sum(~np.isfinite(deni)) > 0:
+ warnings.warn(
+ """Infinite value in iterative scheme, returning NaN.
+ Try rerunning with more samples.""")
+ r = (N1/N2) * np.sum(numi)/np.sum(deni)
+ r_vals.append(r)
+ logml = np.log(r) + lstar
+ i += 1
+ if criterion == 'r':
+ criterion_val = np.abs((r - rold)/r)
+ elif criterion == 'logml':
+ criterion_val = np.abs((logml - logmlold)/logml)
+
+ if i >= maxiter:
+ return dict(logml=np.NaN, niter=i, r_vals=np.asarray(r_vals))
+ else:
+ return dict(logml=logml, niter=i)
+
+# -------------------------------------------------------------------------
+def _my_ESS(x):
+ """
+ Compute the effective sample size of estimand of interest.
+ Vectorised implementation.
+ https://github.com/jwalton3141/jwalton3141.github.io/blob/master/assets/posts/ESS/rwmh.py
+
+
+ Parameters
+ ----------
+ x : array of shape (n_walkers, n_steps)
+ MCMC Samples.
+
+ Returns
+ -------
+ int
+ Effective sample size.
+
+ """
+ # Format the intput
+ try:
+ x = np.array(x)
+ except(AttributeError):
+ print('The input cannot be formatted!')
+ if x.ndim >2 or x.ndim <2:
+ raise AttributeError('The input should be 2D!')
+ m_chains, n_iters = x.shape
+
+ def variogram(t):
+ variogram = ((x[:, t:] - x[:, :(n_iters - t)])**2).sum()
+ variogram /= (m_chains * (n_iters - t))
+ return variogram
+
+ post_var = gelman_rubin(x, return_var=True)
+
+ t = 1
+ rho = np.ones(n_iters)
+ negative_autocorr = False
+
+ # Iterate until the sum of consecutive estimates of autocorrelation is
+ # negative
+ while not negative_autocorr and (t < n_iters):
+ rho[t] = 1 - variogram(t) / (2 * post_var)
+
+ if not t % 2:
+ negative_autocorr = sum(rho[t-1:t+1]) < 0
+
+ t += 1
+
+ return int(m_chains*n_iters / (1 + 2*rho[1:t].sum()))
+
+
class MCMC:
"""
A class for bayesian inference via a Markov-Chain Monte-Carlo (MCMC)
@@ -41,62 +217,64 @@ class MCMC:
"""
def __init__(self, BayesOpts):
-
+ # Inputs
self.BayesOpts = BayesOpts
+
+ # Param inits
+ self.counter = 0
+ self.observation = None
+ self.total_sigma2 = None
+
+ # Get general params from BayesOpts
+ self.out_dir = self.BayesOpts.out_dir
+
+ # Get MCMC parameters ftom BayesOpts
+ pars = self.BayesOpts.mcmc_params
+ self.initsamples = pars['init_samples']
+ if isinstance(self.initsamples, pd.DataFrame):
+ self.initsamples = self.initsamples.values
+ self.nsteps = int(pars['n_steps'])
+ self.nwalkers = int(pars['n_walkers'])
+ self.nburn = pars['n_burn']
+ self.moves = pars['moves']
+ self.mp = pars['multiprocessing']
+ self.verbose = pars['verbose']
def run_sampler(self, observation, total_sigma2):
+ """
+ Run the MCMC sampler for the given observations and stdevs.
+ Parameters
+ ----------
+ observation : TYPE
+ DESCRIPTION.
+ total_sigma2 : TYPE
+ DESCRIPTION.
+
+ Returns
+ -------
+ Posterior_df : TYPE
+ DESCRIPTION.
+
+ """
+ # Get init values
BayesObj = self.BayesOpts
- MetaModel = BayesObj.engine.MetaModel
- Model = BayesObj.engine.Model
Discrepancy = self.BayesOpts.Discrepancy
- n_cpus = Model.n_cpus
- priorDist = BayesObj.engine.ExpDesign.JDist
- ndim = MetaModel.n_params
- self.counter = 0
- output_dir = f'Outputs_Bayes_{Model.name}_{self.BayesOpts.name}'
- if not os.path.exists(output_dir):
- os.makedirs(output_dir)
+ n_cpus = BayesObj.engine.Model.n_cpus
+ ndim = BayesObj.engine.MetaModel.n_params
+ if not os.path.exists(self.out_dir):
+ os.makedirs(self.out_dir)
+ # Save inputs
self.observation = observation
self.total_sigma2 = total_sigma2
- # Unpack mcmc parameters given to BayesObj.mcmc_params
- self.initsamples = None
- self.nwalkers = 100
- self.nburn = 200
- self.nsteps = 100000
- self.moves = None
- self.mp = False
- self.verbose = False
-
- # Extract initial samples
- if 'init_samples' in BayesObj.mcmc_params:
- self.initsamples = BayesObj.mcmc_params['init_samples']
- if isinstance(self.initsamples, pd.DataFrame):
- self.initsamples = self.initsamples.values
-
- # Extract number of steps per walker
- if 'n_steps' in BayesObj.mcmc_params:
- self.nsteps = int(BayesObj.mcmc_params['n_steps'])
- # Extract number of walkers (chains)
- if 'n_walkers' in BayesObj.mcmc_params:
- self.nwalkers = int(BayesObj.mcmc_params['n_walkers'])
- # Extract moves
- if 'moves' in BayesObj.mcmc_params:
- self.moves = BayesObj.mcmc_params['moves']
- # Extract multiprocessing
- if 'multiprocessing' in BayesObj.mcmc_params:
- self.mp = BayesObj.mcmc_params['multiprocessing']
- # Extract verbose
- if 'verbose' in BayesObj.mcmc_params:
- self.verbose = BayesObj.mcmc_params['verbose']
-
# Set initial samples
np.random.seed(0)
if self.initsamples is None:
try:
- initsamples = priorDist.sample(self.nwalkers).T
+ initsamples = BayesObj.engine.ExpDesign.JDist.sample(self.nwalkers).T
+ initsamples = np.swapaxes(np.array([initsamples]),0,1) # TODO: test if this still works with multiple input dists
except:
# when aPCE selected - gaussian kernel distribution
inputSamples = self.BayesOpts.engine.ExpDesign.raw_data.T
@@ -125,7 +303,7 @@ class MCMC:
initsamples[:, idx_dim] = dist.rvs(size=self.nwalkers)
# Update lower and upper
- MetaModel.ExpDesign.bound_tuples = bound_tuples
+ BayesObj.engine.MetaModel.ExpDesign.bound_tuples = bound_tuples
# Check if sigma^2 needs to be inferred
if Discrepancy.opt_sigma != 'B':
@@ -133,8 +311,6 @@ class MCMC:
# Update initsamples
initsamples = np.hstack((initsamples, sigma2_samples))
-
- # Update ndim
ndim = initsamples.shape[1]
# Discrepancy bound
@@ -146,10 +322,8 @@ class MCMC:
print("\n>>>> Bayesian inference with MCMC for "
f"{self.BayesOpts.name} started. <<<<<<")
- # Set up the backend
- filename = f"{output_dir}/emcee_sampler.h5"
- backend = emcee.backends.HDFBackend(filename)
- # Clear the backend in case the file already exists
+ # Set up the backend and clear it in case the file already exists
+ backend = emcee.backends.HDFBackend(f"{self.out_dir}/emcee_sampler.h5")
backend.reset(self.nwalkers, ndim)
# Define emcee sampler
@@ -176,8 +350,8 @@ class MCMC:
)
# Reset sampler
- sampler.reset()
pos = pos.coords
+ sampler.reset()
else:
pos = initsamples
@@ -252,13 +426,13 @@ class MCMC:
# output current autocorrelation estimate
if self.verbose:
print(f"Mean autocorr. time estimate: {np.nanmean(tau):.3f}")
- list_gr = np.round(self.gelman_rubin(sampler.chain), 3)
+ list_gr = np.round(gelman_rubin(sampler.chain), 3)
print("Gelman-Rubin Test*: ", list_gr)
# check convergence
converged = np.all(tau*autocorreverynsteps < sampler.iteration)
converged &= np.all(np.abs(tauold - tau) / tau < 0.01)
- converged &= np.all(self.gelman_rubin(sampler.chain) < 1.1)
+ converged &= np.all(gelman_rubin(sampler.chain) < 1.1)
if converged:
break
@@ -277,7 +451,7 @@ class MCMC:
thin = int(0.5*np.nanmin(tau)) if int(0.5*np.nanmin(tau)) != 0 else 1
finalsamples = sampler.get_chain(discard=burnin, flat=True, thin=thin)
acc_fr = np.nanmean(sampler.acceptance_fraction)
- list_gr = np.round(self.gelman_rubin(sampler.chain[:, burnin:]), 3)
+ list_gr = np.round(gelman_rubin(sampler.chain[:, burnin:]), 3)
# Print summary
print('\n')
@@ -307,7 +481,7 @@ class MCMC:
# Plot traces
if self.verbose and self.nsteps < 10000:
- pdf = PdfPages(output_dir+'/traceplots.pdf')
+ pdf = PdfPages(self.out_dir+'/traceplots.pdf')
fig = plt.figure()
for parIdx in range(ndim):
# Set up the axes with gridspec
@@ -334,7 +508,6 @@ class MCMC:
# Destroy the current plot
plt.clf()
-
pdf.close()
# plot development of autocorrelation estimate
@@ -348,33 +521,9 @@ class MCMC:
plt.ylim(0, np.nanmax(taus)+0.1*(np.nanmax(taus)-np.nanmin(taus)))
plt.xlabel("number of steps")
plt.ylabel(r"mean $\hat{\tau}$")
- fig1.savefig(f"{output_dir}/autocorrelation_time.pdf",
+ fig1.savefig(f"{self.out_dir}/autocorrelation_time.pdf",
bbox_inches='tight')
- # logml_dict = self.marginal_llk_emcee(sampler, self.nburn, logp=None,
- # maxiter=5000)
- # print('\nThe Bridge Sampling Estimation is "
- # f"{logml_dict['logml']:.5f}.')
-
- # # Posterior-based expectation of posterior probablity
- # postExpPostLikelihoods = np.mean(sampler.get_log_prob(flat=True)
- # [self.nburn*self.nwalkers:])
-
- # # Posterior-based expectation of prior densities
- # postExpPrior = np.mean(self.log_prior(emcee_trace.T))
-
- # # Posterior-based expectation of likelihoods
- # postExpLikelihoods_emcee = postExpPostLikelihoods - postExpPrior
-
- # # Calculate Kullback-Leibler Divergence
- # KLD_emcee = postExpLikelihoods_emcee - logml_dict['logml']
- # print("Kullback-Leibler divergence: %.5f"%KLD_emcee)
-
- # # Information Entropy based on Entropy paper Eq. 38
- # infEntropy_emcee = logml_dict['logml'] - postExpPrior -
- # postExpLikelihoods_emcee
- # print("Information Entropy: %.5f" %infEntropy_emcee)
-
Posterior_df = pd.DataFrame(finalsamples, columns=par_names)
return Posterior_df
@@ -397,8 +546,7 @@ class MCMC:
returned otherwise an array.
"""
-
- MetaModel = self.BayesOpts.MetaModel
+ MetaModel = self.BayesOpts.engine.MetaModel
Discrepancy = self.BayesOpts.Discrepancy
# Find the number of sigma2 parameters
@@ -417,7 +565,7 @@ class MCMC:
for i in range(nsamples):
# Check if the sample is within the parameters' range
- if self._check_ranges(theta[i], params_range):
+ if _check_ranges(theta[i], params_range):
# Check if all dists are uniform, if yes priors are equal.
if all(MetaModel.input_obj.Marginals[i].dist_type == 'uniform'
for i in range(MetaModel.n_params)):
@@ -429,7 +577,7 @@ class MCMC:
# Check if bias term needs to be inferred
if Discrepancy.opt_sigma != 'B':
- if self._check_ranges(theta[i, -n_sigma2:],
+ if _check_ranges(theta[i, -n_sigma2:],
disc_bound_tuples):
if all('unif' in disc_marginals[i].dist_type for i in
range(Discrepancy.ExpDesign.ndim)):
@@ -463,21 +611,20 @@ class MCMC:
"""
BayesOpts = self.BayesOpts
- MetaModel = BayesOpts.MetaModel
+ MetaModel = BayesOpts.engine.MetaModel
Discrepancy = self.BayesOpts.Discrepancy
# Find the number of sigma2 parameters
if Discrepancy.opt_sigma != 'B':
disc_bound_tuples = Discrepancy.ExpDesign.bound_tuples
n_sigma2 = len(disc_bound_tuples)
- else:
- n_sigma2 = -len(theta)
- # Check if bias term needs to be inferred
- if Discrepancy.opt_sigma != 'B':
+ # Check if bias term should be inferred
sigma2 = theta[:, -n_sigma2:]
theta = theta[:, :-n_sigma2]
else:
+ n_sigma2 = -len(theta)
sigma2 = None
+
theta = theta if theta.ndim != 1 else theta.reshape((1, -1))
# Evaluate Model/MetaModel at theta
@@ -561,12 +708,11 @@ class MCMC:
"""
BayesObj = self.BayesOpts
- MetaModel = BayesObj.MetaModel
Model = BayesObj.engine.Model
if BayesObj.emulator:
# Evaluate the MetaModel
- mean_pred, std_pred = MetaModel.eval_metamodel(samples=theta)
+ mean_pred, std_pred = BayesObj.engine.MetaModel.eval_metamodel(samples=theta)
else:
# Evaluate the origModel
mean_pred, std_pred = dict(), dict()
@@ -610,8 +756,7 @@ class MCMC:
A error model.
"""
- BayesObj = self.BayesOpts
- MetaModel = BayesObj.MetaModel
+ MetaModel = self.BayesOpts.engine.MetaModel
# Prepare the poster samples
try:
@@ -636,60 +781,10 @@ class MCMC:
# Train a GPR meta-model using MAP
error_MetaModel = MetaModel.create_model_error(
- BayesObj.BiasInputs, y_map, name='Calib')
+ self.BayesOpts.BiasInputs, y_map, name='Calib')
return error_MetaModel
- # -------------------------------------------------------------------------
- def gelman_rubin(self, chain, return_var=False):
- """
- The potential scale reduction factor (PSRF) defined by the variance
- within one chain, W, with the variance between chains B.
- Both variances are combined in a weighted sum to obtain an estimate of
- the variance of a parameter \\( \\theta \\).The square root of the
- ratio of this estimates variance to the within chain variance is called
- the potential scale reduction.
- For a well converged chain it should approach 1. Values greater than
- 1.1 typically indicate that the chains have not yet fully converged.
-
- Source: http://joergdietrich.github.io/emcee-convergence.html
-
- https://github.com/jwalton3141/jwalton3141.github.io/blob/master/assets/posts/ESS/rwmh.py
-
- Parameters
- ----------
- chain : array (n_walkers, n_steps, n_params)
- The emcee ensamples.
-
- Returns
- -------
- R_hat : float
- The Gelman-Robin values.
-
- """
- m_chains, n_iters = chain.shape[:2]
-
- # Calculate between-chain variance
- θb = np.mean(chain, axis=1)
- θbb = np.mean(θb, axis=0)
- B_over_n = ((θbb - θb)**2).sum(axis=0)
- B_over_n /= (m_chains - 1)
-
- # Calculate within-chain variances
- ssq = np.var(chain, axis=1, ddof=1)
- W = np.mean(ssq, axis=0)
-
- # (over) estimate of variance
- var_θ = W * (n_iters - 1) / n_iters + B_over_n
-
- if return_var:
- return var_θ
- else:
- # The square root of the ratio of this estimates variance to the
- # within chain variance
- R_hat = np.sqrt(var_θ / W)
- return R_hat
-
# -------------------------------------------------------------------------
def marginal_llk_emcee(self, sampler, nburn=None, logp=None, maxiter=1000):
"""
@@ -748,7 +843,7 @@ class MCMC:
samples_4_iter[var, :] = x2.flatten()
# effective sample size of samples_4_iter, scalar
- effective_n[var] = self._my_ESS(x2)
+ effective_n[var] = _my_ESS(x2)
# median effective sample size (scalar)
neff = np.median(effective_n)
@@ -772,7 +867,7 @@ class MCMC:
q21 = logp(gen_samples.T)
# Run iterative scheme:
- tmp = self._iterative_scheme(
+ tmp = _iterative_scheme(
N1, N2, q11, q12, q21, q22, r0, neff, tol1, maxiter, 'r'
)
if ~np.isfinite(tmp['logml']):
@@ -782,7 +877,7 @@ class MCMC:
" usual.")
# use geometric mean as starting value
r0_2 = np.sqrt(tmp['r_vals'][-2]*tmp['r_vals'][-1])
- tmp = self._iterative_scheme(
+ tmp = _iterative_scheme(
q11, q12, q21, q22, r0_2, neff, tol2, maxiter, 'logml'
)
@@ -791,119 +886,3 @@ class MCMC:
q11=q11, q12=q12, q21=q21, q22=q22
)
return marg_llk
-
- # -------------------------------------------------------------------------
- def _iterative_scheme(self, N1, N2, q11, q12, q21, q22, r0, neff, tol,
- maxiter, criterion):
- """
- Iterative scheme as proposed in Meng and Wong (1996) to estimate the
- marginal likelihood
-
- """
- l1 = q11 - q12
- l2 = q21 - q22
- # To increase numerical stability,
- # subtracting the median of l1 from l1 & l2 later
- lstar = np.median(l1)
- s1 = neff/(neff + N2)
- s2 = N2/(neff + N2)
- r = r0
- r_vals = [r]
- logml = np.log(r) + lstar
- criterion_val = 1 + tol
-
- i = 0
- while (i <= maxiter) & (criterion_val > tol):
- rold = r
- logmlold = logml
- numi = np.exp(l2 - lstar)/(s1 * np.exp(l2 - lstar) + s2 * r)
- deni = 1/(s1 * np.exp(l1 - lstar) + s2 * r)
- if np.sum(~np.isfinite(numi))+np.sum(~np.isfinite(deni)) > 0:
- warnings.warn(
- """Infinite value in iterative scheme, returning NaN.
- Try rerunning with more samples.""")
- r = (N1/N2) * np.sum(numi)/np.sum(deni)
- r_vals.append(r)
- logml = np.log(r) + lstar
- i += 1
- if criterion == 'r':
- criterion_val = np.abs((r - rold)/r)
- elif criterion == 'logml':
- criterion_val = np.abs((logml - logmlold)/logml)
-
- if i >= maxiter:
- return dict(logml=np.NaN, niter=i, r_vals=np.asarray(r_vals))
- else:
- return dict(logml=logml, niter=i)
-
- # -------------------------------------------------------------------------
- def _my_ESS(self, x):
- """
- Compute the effective sample size of estimand of interest.
- Vectorised implementation.
- https://github.com/jwalton3141/jwalton3141.github.io/blob/master/assets/posts/ESS/rwmh.py
-
-
- Parameters
- ----------
- x : array of shape (n_walkers, n_steps)
- MCMC Samples.
-
- Returns
- -------
- int
- Effective sample size.
-
- """
- m_chains, n_iters = x.shape
-
- def variogram(t):
- variogram = ((x[:, t:] - x[:, :(n_iters - t)])**2).sum()
- variogram /= (m_chains * (n_iters - t))
- return variogram
-
- post_var = self.gelman_rubin(x, return_var=True)
-
- t = 1
- rho = np.ones(n_iters)
- negative_autocorr = False
-
- # Iterate until the sum of consecutive estimates of autocorrelation is
- # negative
- while not negative_autocorr and (t < n_iters):
- rho[t] = 1 - variogram(t) / (2 * post_var)
-
- if not t % 2:
- negative_autocorr = sum(rho[t-1:t+1]) < 0
-
- t += 1
-
- return int(m_chains*n_iters / (1 + 2*rho[1:t].sum()))
-
- # -------------------------------------------------------------------------
- def _check_ranges(self, theta, ranges):
- """
- This function checks if theta lies in the given ranges.
-
- Parameters
- ----------
- theta : array
- Proposed parameter set.
- ranges : nested list
- List of the praremeter ranges.
-
- Returns
- -------
- c : bool
- If it lies in the given range, it return True else False.
-
- """
- c = True
- # traverse in the list1
- for i, bounds in enumerate(ranges):
- x = theta[i]
- # condition check
- if x < bounds[0] or x > bounds[1]:
- c = False
- return c
- return c
diff --git a/src/bayesvalidrox/pylink/pylink.py b/src/bayesvalidrox/pylink/pylink.py
index 227a51ab38cd834e7e85f6193d83563c7ed3437a..637f42317e6f97815e51ce0b331b61c03f26a85b 100644
--- a/src/bayesvalidrox/pylink/pylink.py
+++ b/src/bayesvalidrox/pylink/pylink.py
@@ -231,7 +231,7 @@ class PyLinkForwardModel(object):
self.observations_valid = self.observations_valid
else:
raise Exception("Please provide the observation data as a "
- "dictionary via observations attribute or pass"
+ "dictionary via observations_valid attribute or pass"
" the csv-file path to MeasurementFile "
"attribute")
# Compute the number of observation
diff --git a/src/bayesvalidrox/surrogate_models/surrogate_models.py b/src/bayesvalidrox/surrogate_models/surrogate_models.py
index a13a96cc53da289b08f05683234c7d7203b6ed8e..5968d3a06bcf0281c39f06ea43f647a19d4cfaa0 100644
--- a/src/bayesvalidrox/surrogate_models/surrogate_models.py
+++ b/src/bayesvalidrox/surrogate_models/surrogate_models.py
@@ -1305,6 +1305,7 @@ class MetaModel():
extracted data.
y : array of shape (n_outputs,)
The model response for the MAP parameter set.
+ MeasuredData :
name : str, optional
Calibration or validation. The default is `'Calib'`.
@@ -1325,6 +1326,8 @@ class MetaModel():
# Fitting GPR based bias model
for out in outputNames:
+ print(out)
+ print(X[out])
nan_idx = ~np.isnan(MeasuredData[out])
# Select data
try:
diff --git a/tests/test_BayesInference.py b/tests/test_BayesInference.py
index 7fbb79195799724a8ae820a2c777c99c0b86d9b0..70bcdbc40c7748f5d58936981b354fcff2b933f8 100644
--- a/tests/test_BayesInference.py
+++ b/tests/test_BayesInference.py
@@ -6,17 +6,21 @@ Tests are available for the following functions
_logpdf - x
_kernel_rbf - x
class BayesInference:
- create_inference
- perform_bootstrap
+ setup_inference - x
+ create_inference - x
+ perform_bootstrap Need working model for tests without emulator
_perturb_data - x
+ create_error_model Error in the MetaModel
_eval_model Need working model to test this
- normpdf
+ normpdf - x
_corr_factor_BME_old - removed
_corr_factor_BME - x
_rejection_sampling - x
- _posterior_predictive
- _plot_max_a_posteriori
- _plot_post_predictive
+ _posterior_predictive - x
+ plot_post_params - x
+ plot_log_BME - x
+ _plot_max_a_posteriori Need working model to test this
+ _plot_post_predictive - x
"""
import sys
sys.path.append("src/")
@@ -360,6 +364,184 @@ def test_normpdf_allsigmas() -> None:
bi = BayesInference(engine)
bi.normpdf(expdes.Y, obs_data, total_sigma2s, sigma2=sigma2, std=total_sigma2s)
+#%% Test setup_inference
+
+def test_setup_inference_noobservation() -> None:
+ """
+ Test the object setup without given observations
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+
+ mod = PL()
+ mod.Output.names = ['Z']
+
+ engine = Engine(mm, mod, expdes)
+
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = DiscrepancyOpts
+ with pytest.raises(Exception) as excinfo:
+ bi.setup_inference()
+ assert str(excinfo.value) == 'Please provide the observation data as a dictionary via observations attribute or pass the csv-file path to MeasurementFile attribute'
+
+def test_setup_inference() -> None:
+ """
+ Test the object setup with observations
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+
+ mod = PL()
+ mod.observations = {'Z':np.array([0.45])}
+ mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])}
+ mod.Output.names = ['Z']
+
+ engine = Engine(mm, mod, expdes)
+
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = DiscrepancyOpts
+ bi.setup_inference()
+
+def test_setup_inference_priorsamples() -> None:
+ """
+ Test the object setup with prior samples set by hand
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+
+ mod = PL()
+ mod.observations = {'Z':np.array([0.45])}
+ mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])}
+ mod.Output.names = ['Z']
+
+ engine = Engine(mm, mod, expdes)
+
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+
+ bi = BayesInference(engine)
+ bi.prior_samples = np.swapaxes(np.array([np.random.normal(0,1,100)]),0,1)
+ bi.Discrepancy = DiscrepancyOpts
+ bi.setup_inference()
+
+def test_setup_inference_valid() -> None:
+ """
+ Test the object setup for valid
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+
+ mod = PL()
+ mod.observations_valid = {'Z':np.array([0.45])}
+ mod.observations_valid = {'Z':np.array([0.45]), 'x_values':np.array([0])}
+ mod.Output.names = ['Z']
+
+ engine = Engine(mm, mod, expdes)
+
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = DiscrepancyOpts
+ bi.name = 'valid'
+ bi.setup_inference()
+
+def test_setup_inference_noname() -> None:
+ """
+ Test the object setup for an invalid inference name
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+
+ mod = PL()
+ mod.observations = {'Z':np.array([0.45])}
+ mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])}
+ mod.Output.names = ['Z']
+
+ engine = Engine(mm, mod, expdes)
+
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = DiscrepancyOpts
+ bi.name = ''
+ with pytest.raises(Exception) as excinfo:
+ bi.setup_inference()
+ assert str(excinfo.value) == 'The set inference type is not known! Use either `calib` or `valid`'
+
+
#%% Test perform_bootstrap
@@ -367,12 +549,51 @@ def test_perform_bootstrap() -> None:
"""
Do bootstrap
"""
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
- opt_sigma = 'B'
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+ expdes.n_max_samples = 4
+ expdes.X = np.array([[0],[1],[0.5]])
+ expdes.Y = {'Z':[[0.4],[0.5],[0.45]]}
+ expdes.x_values = np.array([0]) # Error in plots if this is not available
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mm.fit(expdes.X, expdes.Y)
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+
+ mod = PL()
+ mod.observations = {'Z':np.array([0.45])}
+ mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])} # Error if x_values not given
+ mod.Output.names = ['Z']
+ mod.n_obs = 1
+
+ engine = Engine(mm, mod, expdes)
+
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = DiscrepancyOpts
+ bi.bootstrap = True
+ bi.plot_post_pred = False
total_sigma2s = {'Z':np.array([0.15])}
- bi.perform_bootstrap(opt_sigma, total_sigma2s)
+ bi.setup_inference()
+ bi.perform_bootstrap(total_sigma2s)
-if __name__ == '__main__':
+
+def test_perform_bootstrap_bayesloocv() -> None:
+ """
+ Do bootstrap
+ """
inp = Input()
inp.add_marginals()
inp.Marginals[0].dist_type = 'normal'
@@ -388,35 +609,493 @@ if __name__ == '__main__':
mm = MetaModel(inp)
mm.n_params = 1
mm.fit(expdes.X, expdes.Y)
- expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(mm.pce_deg))
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+
+ mod = PL()
+ mod.observations = {'Z':np.array([0.45])}
+ mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])} # Error if x_values not given
+ mod.Output.names = ['Z']
+ mod.n_obs = 1
+
+ engine = Engine(mm, mod, expdes)
+
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = DiscrepancyOpts
+ bi.bootstrap = True
+ bi.plot_post_pred = False
+ total_sigma2s = {'Z':np.array([0.15])}
+ bi.setup_inference()
+ bi.bayes_loocv = True
+ bi.perform_bootstrap(total_sigma2s)
+
+#%% Test create_error_model
+
+def create_error_model_prior() -> None:
+ """
+ Test creating MetaModel error-model for 'prior'
+ """
+ # TODO: there are issues with the expected formats from the MetaModel
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+ expdes.n_max_samples = 4
+ expdes.X = np.array([[0],[1],[0.5]])
+ expdes.Y = {'Z':[[0.4],[0.5],[0.45]]}
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mm.fit(expdes.X, expdes.Y)
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+
+ mod = PL()
+ mod.observations = {'Z':np.array([0.45])}
+ mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])}
+ mod.Output.names = ['Z']
+ mod.n_obs = 1
+
+ engine = Engine(mm, mod, expdes)
+
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = DiscrepancyOpts
+ bi.bootstrap = True
+ bi.setup_inference()
+ bi.bias_inputs = expdes.X
+ bi.create_error_model(type_ = 'prior', opt_sigma = 'B', sampler = None)
+
+def create_error_model_posterior() -> None:
+ """
+ Test creating MetaModel error-model for 'posterior'
+ """
+ # TODO: there are issues with the expected formats from the MetaModel
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+ expdes.n_max_samples = 4
+ expdes.X = np.array([[0],[1],[0.5]])
+ expdes.Y = {'Z':[[0.4],[0.5],[0.45]]}
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mm.fit(expdes.X, expdes.Y)
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+
+ mod = PL()
+ mod.observations = {'Z':np.array([0.45])}
+ mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])}
+ mod.Output.names = ['Z']
+ mod.n_obs = 1
+
+ engine = Engine(mm, mod, expdes)
+
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+
+ posterior = pd.DataFrame()
+ posterior[None] = [0,1,0.5]
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = DiscrepancyOpts
+ bi.bootstrap = True
+ bi.setup_inference()
+ bi.bias_inputs = expdes.X
+ bi.posterior_df = posterior
+ bi.create_error_model(type_ = 'posterior', opt_sigma = 'B', sampler = None)
+
+#%% Test _posterior_predictive
+
+def test_posterior_predictive() -> None:
+ """
+ Test posterior predictions
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+ prior_samples = np.swapaxes(np.array([np.random.normal(0,1,10)]),0,1)
+
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+ expdes.n_max_samples = 4
+ expdes.X = np.array([[0],[1],[0.5]])
+ expdes.Y = {'Z':[[0.4],[0.5],[0.45]]}
+ expdes.x_values = np.array([0]) # Error in plots if this is not available
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mm.fit(expdes.X, expdes.Y)
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+ y_hat, y_std = mm.eval_metamodel(prior_samples)
+
+ mod = PL()
+ mod.observations = {'Z':np.array([0.45])}
+ mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])} # Error if x_values not given
+ mod.Output.names = ['Z']
+ mod.n_obs = 1
+
+ engine = Engine(mm, mod, expdes)
+
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ total_sigma2s = {'Z':np.array([0.15])}
+ posterior = pd.DataFrame()
+ posterior[None] = [0,1,0.5]
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = DiscrepancyOpts
+ bi.bootstrap = True
+ bi.plot_post_pred = False
+ bi.posterior_df = posterior
+ bi.bias_inputs = expdes.X
+ bi._mean_pce_prior_pred = y_hat
+ bi._std_pce_prior_pred = y_std
+ bi.Discrepancy.total_sigma2 = total_sigma2s
+ bi.setup_inference()
+ bi._posterior_predictive()
+def test_posterior_predictive_rejection() -> None:
+ """
+ Test posterior predictions with rejection inference
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+ prior_samples = np.swapaxes(np.array([np.random.normal(0,1,10)]),0,1)
+
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+ expdes.n_max_samples = 4
+ expdes.X = np.array([[0],[1],[0.5]])
+ expdes.Y = {'Z':[[0.4],[0.5],[0.45]]}
+ expdes.x_values = np.array([0]) # Error in plots if this is not available
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mm.fit(expdes.X, expdes.Y)
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+ y_hat, y_std = mm.eval_metamodel(prior_samples)
mod = PL()
mod.observations = {'Z':np.array([0.45])}
mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])} # Error if x_values not given
mod.Output.names = ['Z']
+ mod.n_obs = 1
engine = Engine(mm, mod, expdes)
sigma2Dict = {'Z':np.array([0.05])}
sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ total_sigma2s = {'Z':np.array([0.15])}
+ posterior = pd.DataFrame()
+ posterior[None] = [0,1,0.5]
obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
DiscrepancyOpts = Discrepancy('')
DiscrepancyOpts.type = 'Gaussian'
DiscrepancyOpts.parameters = (obsData*0.15)**2
+ bi = BayesInference(engine)
+ bi.Discrepancy = DiscrepancyOpts
+ bi.bootstrap = True
+ bi.plot_post_pred = False
+ bi.posterior_df = posterior
+ bi.bias_inputs = expdes.X
+ bi._mean_pce_prior_pred = y_hat
+ bi._std_pce_prior_pred = y_std
+ bi.Discrepancy.total_sigma2 = total_sigma2s
+ bi.inference_method = 'rejection'
+ bi.setup_inference()
+ bi._posterior_predictive()
+
+#%% Test plot_post_params
+
+def test_plot_post_params() -> None:
+ """
+ Plot posterior dist
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.init_param_space(max_deg=1)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mod = PL()
+ engine = Engine(mm, mod, expdes)
+
+ bi = BayesInference(engine)
+ posterior = pd.DataFrame()
+ posterior[None] = [0,1,0.5]
+ bi.posterior_df = posterior
+ bi.plot_post_params('B')
+
+def test_plot_post_params_noemulator() -> None:
+ """
+ Plot posterior dist with emulator = False
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.init_param_space(max_deg=1)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mod = PL()
+ engine = Engine(mm, mod, expdes)
+
+ bi = BayesInference(engine)
+ posterior = pd.DataFrame()
+ posterior[None] = [0,1,0.5]
+ bi.posterior_df = posterior
+ bi.emulator = False
+ bi.plot_post_params('B')
+
+
+#%% Test plot_log_BME
+
+def test_plot_log_BME() -> None:
+ """
+ Show the log_BME from bootstrapping
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.init_param_space(max_deg=1)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mod = PL()
+ engine = Engine(mm, mod, expdes)
+
+ bi = BayesInference(engine)
+ bi.log_BME = np.array([[0,0.2],[0,0.2]])
+ bi.n_tot_measurement = 1
+ bi.plot_log_BME()
+
+def test_plot_log_BME_noemulator() -> None:
+ """
+ Show the log_BME from bootstrapping with emulator = False
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.init_param_space(max_deg=1)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mod = PL()
+ engine = Engine(mm, mod, expdes)
+
+ bi = BayesInference(engine)
+ bi.log_BME = np.array([[0,0.2],[0,0.2]])
+ bi.n_tot_measurement = 1
+ bi.emulator = False
+ bi.plot_log_BME()
+
+#%% Test _plot_max_a_posteriori
+
+def test_plot_max_a_posteriori_rejection() -> None:
+ """
+ Plot MAP estimate for rejection
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.init_param_space(max_deg=1)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mod = PL()
+ engine = Engine(mm, mod, expdes)
+
+ bi = BayesInference(engine)
+ bi.inference_method = 'rejection'
+ bi._plot_post_predictive()
+
+
+def test_plot_max_a_posteriori() -> None:
+ """
+ Plot MAP estimate
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.init_param_space(max_deg=1)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mod = PL()
+ engine = Engine(mm, mod, expdes)
+
+ bi = BayesInference(engine)
+ bi._plot_post_predictive()
+
+
+#%% Test _plot_post_predictive
+
+
+def test_plot_post_predictive_rejection() -> None:
+ """
+ Plot posterior predictions for rejection
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.init_param_space(max_deg=1)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mod = PL()
+ engine = Engine(mm, mod, expdes)
+
+ bi = BayesInference(engine)
+ bi.inference_method = 'rejection'
+ bi._plot_post_predictive()
+
+
+def test_plot_post_predictive() -> None:
+ """
+ Plot posterior predictions
+ """
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.init_param_space(max_deg=1)
+ expdes.n_init_samples = 2
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mod = PL()
+ engine = Engine(mm, mod, expdes)
+
+ bi = BayesInference(engine)
+ bi._plot_post_predictive()
+
+
+#%% Main runs
+if __name__ == '__main__':
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+ #prior_samples = np.swapaxes(np.array([np.random.normal(0,1,10)]),0,1)
+
+ expdes = ExpDesigns(inp)
+ expdes.init_param_space(max_deg=1)
+ expdes.n_init_samples = 2
+ #expdes.n_max_samples = 4
+ expdes.X = np.array([[0],[1],[0.5]])
+ expdes.Y = {'Z':[[0.4],[0.5],[0.45]]}
+ #expdes.x_values = np.array([0]) # Error in plots if this is not
+
+ mm = MetaModel(inp)
+ mm.n_params = 1
+ mm.fit(expdes.X, expdes.Y)
+ #expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(1))
+ #y_hat, y_std = mm.eval_metamodel(prior_samples)
+
+ mod = PL()
+ #mod.observations = {'Z':np.array([0.45])}
+ #mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])} # Error if x_values not given
+ #mod.Output.names = ['Z']
+ #mod.n_obs = 1
+
+ engine = Engine(mm, mod, expdes)
+
+ #sigma2Dict = {'Z':np.array([0.05])}
+ #sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ obsData = pd.DataFrame({'Z':np.array([0.45]), 'x_values':np.array([0])}, columns=mod.Output.names)
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+ DiscrepancyOpts.opt_sigma = 'B'
+
bi = BayesInference(engine)
bi.Discrepancy = DiscrepancyOpts # Error if this not class 'DiscrepancyOpts' or dict(?)
- bi.bootstrap = True # Error if this and bayes_loocv and just_analysis are all False?
- bi.plot_post_pred = False # Remaining issue in the violinplot
+ #bi.bootstrap = True # Error if this and bayes_loocv and just_analysis are all False?
+ #bi.plot_post_pred = False # Remaining issue in the violinplot
#bi.error_model = True
#bi.bayes_loocv = True
- bi.create_inference()
- stop
- opt_sigma = 'B'
- total_sigma2s = {'Z':np.array([0.15])}
- data = pd.DataFrame()
- data['Z'] = [0.45]
- data['x_values'] = [0.3]
- bi.measured_data = data
- bi.perform_bootstrap(opt_sigma, total_sigma2s)
-
\ No newline at end of file
+ #if 0:
+ # bi.create_inference()
+ #opt_sigma = 'B'
+ #total_sigma2s = {'Z':np.array([0.15])}
+ #data = pd.DataFrame()
+ #data['Z'] = [0.45]
+ #data['x_values'] = [0.3]
+ #bi.setup_inference()
+ #bi.perform_bootstrap(total_sigma2s)
+ posterior = pd.DataFrame()
+ posterior[None] = [0,1,0.5]
+ bi.posterior_df = posterior
+ #bi.bias_inputs = expdes.X
+ #bi._mean_pce_prior_pred = y_hat
+ #bi._std_pce_prior_pred = y_std
+ #bi.Discrepancy.total_sigma2 = total_sigma2s
+ #bi.create_error_model(type_ = 'posterior', opt_sigma = 'B', sampler = None)
+ #bi._posterior_predictive()
+ #bi.plot_post_params('B')
+ #bi.log_BME = np.array([[0,0.2],[0,0.2]])
+ #bi.n_tot_measurement = 1
+ #bi.plot_log_BME()
+ bi.inference_method = 'rejection'
+ bi._plot_max_a_posteriori()
\ No newline at end of file
diff --git a/tests/test_MCMC.py b/tests/test_MCMC.py
index 50f6d69b67449ef0647e777f876dd10210566d8f..762ff1ae8feb735bf91a77c554f3437f03fa9712 100644
--- a/tests/test_MCMC.py
+++ b/tests/test_MCMC.py
@@ -2,6 +2,10 @@
"""
Test the MCM class of bayesvalidrox
Tests are available for the following functions
+_check_ranges - x
+gelmain_rubin
+_iterative_scheme
+_my_ESS - x
Class MCMC:
run_sampler
log_prior
@@ -9,20 +13,25 @@ Class MCMC:
log_posterior
eval_model
train_error_model
- gelmain_rubin
marginal_llk_emcee
- _iterative_scheme
- _my_ESS
- _check_ranges
"""
+import emcee
import sys
sys.path.append("src/")
sys.path.append("../src/")
+import pandas as pd
import pytest
import numpy as np
+from bayesvalidrox.surrogate_models.inputs import Input
+from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
+from bayesvalidrox.surrogate_models.surrogate_models import MetaModel
+from bayesvalidrox.pylink.pylink import PyLinkForwardModel as PL
+from bayesvalidrox.surrogate_models.engine import Engine
+from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
from bayesvalidrox.bayes_inference.mcmc import MCMC
from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
+from bayesvalidrox.bayes_inference.mcmc import _check_ranges, gelman_rubin, _my_ESS, _iterative_scheme
#%% Test MCMC init
@@ -30,15 +39,229 @@ def test_MCMC() -> None:
"""
Construct an MCMC object
"""
- MCMC('')
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
-def test_MCMC() -> None:
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+ expdes.n_max_samples = 4
+ expdes.X = np.array([[0],[1],[0.5]])
+ expdes.Y = {'Z':[[0.4],[0.5],[0.45]]}
+ expdes.x_values = np.array([0])
+
+ mm = MetaModel(inp)
+ mm.fit(expdes.X, expdes.Y)
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(mm.pce_deg))
+
+ mod = PL()
+ mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])}
+ mod.Output.names = ['Z']
+ engine = Engine(mm, mod, expdes)
+
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ disc = Discrepancy('')
+ disc.type = 'Gaussian'
+ disc.parameters = (obsData*0.15)**2
+ disc.opt_sigma = 'B'
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = disc
+ bi.inference_method = 'mcmc'
+ bi.setup_inference()
+ MCMC(bi)
+
+#%% Test run_sampler
+
+def test_run_sampler() -> None:
"""
- Construct an MCMC object
+ Run short MCMC
+
+ Returns
+ -------
+ None
+ DESCRIPTION.
+
"""
- MCMC('')
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+ expdes.n_max_samples = 4
+ expdes.X = np.array([[0],[1],[0.5]])
+ expdes.Y = {'Z':[[0.4],[0.5],[0.45]]}
+ expdes.x_values = np.array([0])
+
+ mm = MetaModel(inp)
+ mm.fit(expdes.X, expdes.Y)
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(mm.pce_deg))
+
+ mod = PL()
+ mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])}
+ mod.Output.names = ['Z']
+ engine = Engine(mm, mod, expdes)
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ disc = Discrepancy('')
+ disc.type = 'Gaussian'
+ disc.parameters = (obsData*0.15)**2
+ disc.opt_sigma = 'B'
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = disc
+ bi.inference_method = 'mcmc'
+ bi.setup_inference()
+ total_sigma2s = {'Z':np.array([0.15])}
+ mcmc = MCMC(bi)
+ mcmc.nburn = 10
+ mcmc.nsteps = 50
+ mcmc.run_sampler(mod.observations, total_sigma2s)
+
+
+#%% Test log_prior
+#%% Test log_likelihood
+
+#%% Test log_posterior
+
+#%% Test eval_model
+
+#%% Test train_error_model
+
+#%% Test gelmain_rubin
+
+def test_gelman_rubin() -> None:
+ """
+ Calculate gelman-rubin
+ """
+
+ chain = [[[1],[2]]]
+ assert gelman_rubin(chain) == np.array([None])
+
+def test_gelman_rubin_returnvar() -> None:
+ """
+ Calculate gelman-rubin returning var
+ """
+
+ chain = [[[1],[2]]]
+ assert gelman_rubin(chain, return_var = True) == np.array([None])
+
+#%% Test marginal_llk_emcee
+
+#%% Test _iterative_scheme
+
+def test_iterative_scheme_r() -> None:
+ """
+ Run iterative calc in mode r
+ """
+ _iterative_scheme(1, 1, 0, 1, 2, 3, 4, 1, 1, 10, 'r')
+
+
+def test_iterative_scheme_0iter() -> None:
+ """
+ Run iterative calc in mode r with 0 iterations
+ """
+ _iterative_scheme(1, 1, 0, 1, 2, 3, 4, 1, 1, 0, 'r')
+
+
+def test_iterative_scheme_logml() -> None:
+ """
+ Run iterative calc in mode logml
+ """
+ _iterative_scheme(1, 1, 0, 1, 2, 3, 4, 1, 1, 10, 'logml')
+
+#%% Test _my_ESS
+
+def test_my_ESS() -> None:
+ """
+ Get ESS
+ """
+ chain = [[[1],[2]]][0]
+ _my_ESS(chain)
+
+def test_my_ESS_1d() -> None:
+ """
+ Get ESS with wrong dimensions
+ """
+ chain = [[[1],[2]]][0][0]
+ with pytest.raises(AttributeError) as excinfo:
+ _my_ESS(chain)
+ assert str(excinfo.value) == 'The input should be 2D!'
+
+#%% Test _check_ranges
+
+def test_check_ranges() -> None:
+ """
+ Check to see if theta lies in expected ranges
+ """
+ theta = [0.5,1.2]
+ ranges = [[0,1],[1,2]]
+ assert _check_ranges(theta, ranges) == True
+
+def test_check_ranges_inv() -> None:
+ """
+ Check to see if theta lies not in expected ranges
+ """
+ theta = [1.5,1.2]
+ ranges = [[0,1],[1,2]]
+ assert _check_ranges(theta, ranges) == False
+
+
+#%% Main
+
+if __name__ == '__main__':
+ inp = Input()
+ inp.add_marginals()
+ inp.Marginals[0].dist_type = 'normal'
+ inp.Marginals[0].parameters = [0,1]
+
+ expdes = ExpDesigns(inp)
+ expdes.n_init_samples = 2
+ expdes.n_max_samples = 4
+ expdes.X = np.array([[0],[1],[0.5]])
+ expdes.Y = {'Z':[[0.4],[0.5],[0.45]]}
+ #expdes.x_values = np.array([0]) # Error in plots if this is not available
+
+ mm = MetaModel(inp)
+ mm.fit(expdes.X, expdes.Y)
+ expdes.generate_ED(expdes.n_init_samples, transform=True, max_pce_deg=np.max(mm.pce_deg))
+
+ mod = PL()
+ mod.observations = {'Z':np.array([0.45]), 'x_values':np.array([0])}
+ mod.Output.names = ['Z']
+
+ engine = Engine(mm, mod, expdes)
+
+ sigma2Dict = {'Z':np.array([0.05])}
+ sigma2Dict = pd.DataFrame(sigma2Dict, columns = ['Z'])
+ obsData = pd.DataFrame(mod.observations, columns=mod.Output.names)
+ disc = Discrepancy('')
+ disc.type = 'Gaussian'
+ disc.parameters = (obsData*0.15)**2
+ disc.opt_sigma = 'B'
+
+ bi = BayesInference(engine)
+ bi.Discrepancy = disc
+ bi.inference_method = 'mcmc'
+ bi.setup_inference()
+
+ #chain = [[[1],[2]]]
+ total_sigma2s = {'Z':np.array([0.15])}
+ mcmc = MCMC(bi)
+ mcmc.nsteps = 50
+ mcmc.nburn = 10
+ mcmc.run_sampler(mod.observations, total_sigma2s)
+ #mcmc.gelmain_rubin(chain)
+ chain = [[[1],[2]]]
+ gelman_rubin(chain, return_var=True)
+ _iterative_scheme(1, 1, 0, 1, 2, 3, 4, 1, 1, 10, 'r')
\ No newline at end of file