{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nested sampling\n",
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Problem definition\n",
    "from probeye.definition.inverse_problem import InverseProblem\n",
    "from probeye.definition.forward_model import ForwardModelBase\n",
    "from probeye.definition.distribution import Normal, Uniform\n",
    "from probeye.definition.sensor import Sensor\n",
    "from probeye.definition.likelihood_model import GaussianLikelihoodModel\n",
    "from probeye.definition.correlation_model import ExpModel\n",
    "\n",
    "# Inference\n",
    "from probeye.inference.dynesty.solver import DynestySolver"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fixed parameters\n",
    "I = 1e9  # mm^4\n",
    "L = 10_000  # mm\n",
    "\n",
    "# Measurements\n",
    "x_sensors = [2500, 5000]  # mm  (always from lower to higher, bug in inv_cov_vec_1D in tripy)\n",
    "d_sensors = [35, 50]  # mm\n",
    "sigma_model = 2.5  # mm\n",
    "pearson = 0.5\n",
    "l_corr = -np.abs(x_sensors[1] - x_sensors[0]) / np.log(pearson)  # mm (assuming exponential correlation)\n",
    "\n",
    "# Prior\n",
    "E_mean = 60  # GPa\n",
    "E_std = 20  # GPa\n",
    "Q_mean = 60  # kN\n",
    "Q_std = 30  # kN\n",
    "Q_loc_low = 0  # mm\n",
    "Q_loc_high = 10000  # mm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define forward model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def beam_deflection(E, Q, a, x):  # a is load position, x is sensor position\n",
    "    if x < a:\n",
    "        b = L - a\n",
    "        return Q * b * x * (L ** 2 - b ** 2 - x ** 2) / (6 * E * I * L)\n",
    "    \n",
    "    return Q * a * (L - x) * (2 * L * x - x ** 2 - a ** 2) / (6 * E * I * L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BeamModel(ForwardModelBase):\n",
    "    def interface(self):\n",
    "        self.parameters = [\"E\", \"Q\", \"a\"]\n",
    "        self.input_sensors = Sensor(\"x\")\n",
    "        self.output_sensors = Sensor(\"y\", std_model=\"sigma\")\n",
    "\n",
    "    def response(self, inp: dict) -> dict:\n",
    "        E = inp[\"E\"]\n",
    "        Q = inp[\"Q\"]\n",
    "        a = inp[\"a\"]\n",
    "        x = inp[\"x\"]\n",
    "        return {\"y\": [beam_deflection(E, Q, a, float(x[0])), beam_deflection(E, Q, a, float(x[1]))]}  # float() needed, probably a bug in probeye"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define inverse problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate the inverse problem\n",
    "problem = InverseProblem(\"Beam model with two sensors\", print_header=False)\n",
    "\n",
    "# Add latent parameters\n",
    "problem.add_parameter(\n",
    "    \"E\",\n",
    "    tex=\"$E$\",\n",
    "    info=\"Elastic modulus of the beam (GPa)\",\n",
    "    prior=Normal(mean=E_mean, std=E_std),\n",
    ")\n",
    "problem.add_parameter(\n",
    "    \"Q\",\n",
    "    tex=\"$Q$\",\n",
    "    info=\"Load applied to the beam (kN)\",\n",
    "    prior=Normal(mean=Q_mean, std=Q_std),\n",
    ")\n",
    "problem.add_parameter(\n",
    "    \"a\",\n",
    "    tex=\"$a$\",\n",
    "    info=\"Position of the load (mm)\",\n",
    "    prior=Uniform(low=Q_loc_low, high=Q_loc_high)\n",
    ")\n",
    "\n",
    "# Add fixed parameters\n",
    "problem.add_parameter(\n",
    "    \"sigma\",\n",
    "    tex=\"$\\sigma$\",\n",
    "    info=\"Standard deviation of the model error (mm)\",\n",
    "    value=sigma_model,\n",
    ")\n",
    "problem.add_parameter(\n",
    "    \"l_corr\",\n",
    "    tex=\"$l_{corr}$\",\n",
    "    info=\"Correlation length of the model error (mm)\",\n",
    "    value=l_corr,\n",
    ")\n",
    "\n",
    "# Add measurement data\n",
    "problem.add_experiment(\n",
    "    name=\"TestSeries_1\",\n",
    "    sensor_data={\"x\": x_sensors, \"y\": d_sensors}\n",
    ")\n",
    "                 \n",
    "# Add forward model\n",
    "problem.add_forward_model(BeamModel(\"BeamModel\"), experiments=\"TestSeries_1\")\n",
    "\n",
    "# Add likelihood model\n",
    "likelihood_model = GaussianLikelihoodModel(\n",
    "        experiment_name=\"TestSeries_1\", \n",
    "        model_error=\"additive\",\n",
    "        correlation=ExpModel(x=\"l_corr\")\n",
    "        )\n",
    "problem.add_likelihood_model(likelihood_model)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve with Dynesty\n",
    "### Dynamic method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Static method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Static with more live points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}