diff --git a/docs/examples/core/linefit/03-SrFit-demo-constraints-restraints.ipynb b/docs/examples/core/linefit/03-SrFit-demo-constraints-restraints.ipynb new file mode 100644 index 0000000..070cd68 --- /dev/null +++ b/docs/examples/core/linefit/03-SrFit-demo-constraints-restraints.ipynb @@ -0,0 +1,522 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# SrFit example for a simple linear fit to a noisy data." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Simulate linear data with some random Gaussian noise and plot the generated \"observed\" data (xobs, yobs)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from __future__ import print_function\n", + "%matplotlib inline\n", + "from matplotlib.pyplot import plot, title\n", + "\n", + "import numpy as np\n", + "xobs = np.arange(-10, 10.1)\n", + "dyobs = 0.3 * np.ones_like(xobs)\n", + "yobs = 0.5 * xobs + 3 + dyobs * np.random.randn(xobs.size)\n", + "plot(xobs, yobs, 'x')\n", + "title('y = 0.5*x + 3 generated with a normal noise at sigma=0.3');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We are going to define a line fitting regression using SrFit.\n", + "At first we create a SrFit Profile object that holds the observed data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from diffpy.srfit.fitbase import Profile\n", + "linedata = Profile()\n", + "linedata.setObservedProfile(xobs, yobs, dyobs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The second step is to create a FitContribution object, which associates observed profile with a mathematical model for the dependent variable." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from diffpy.srfit.fitbase import FitContribution\n", + "linefit = FitContribution('linefit')\n", + "linefit.setProfile(linedata)\n", + "linefit.setEquation(\"A * x + B\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " SrFit objects can be examined by calling their **show()** function. SrFit\n", + " parses the model equation and finds two parameters A, B at independent\n", + " variable x. The values of parameters A, B are at this stage undefined.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "linefit.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can set A and B to some specific values and calculate model\n", + "observations. The x and y attributes of the FitContribution are \n", + "the observed values, which may be re-sampled or truncated to a shorter \n", + "fitting range." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "linefit.A\n", + "linefit.A = 3\n", + "linefit.B = 5\n", + "print(linefit.A, linefit.A.value)\n", + "print(linefit.B, linefit.B.value)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`linefit.evaluate()` returns the modeled values and `linefit.residual()`,\n", + "the difference between observed and modeled data scaled by estimated\n", + "standard deviations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"linefit.evaluate() =\", linefit.evaluate())\n", + "print(\"linefit.residual() =\", linefit.residual())\n", + "plot(xobs, yobs, 'x', linedata.x, linefit.evaluate(), '-')\n", + "title('Line simulated at A=3, B=5')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " We want to find the optimum model parameters that fit the simulated curve\n", + " to the observations. This is done by associating FitContribution with\n", + " a FitRecipe object. FitRecipe can manage multiple fit contributions and\n", + " optimize all models to fit their respective profiles." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from diffpy.srfit.fitbase import FitRecipe\n", + "rec = FitRecipe()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The `clearFitHooks()` function suppresses printout of iteration numbers. The `addContribution()` function includes the specified FitContribution in the FitRecipe, which acts as a top-level manager of all associated fits. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rec.clearFitHooks()\n", + "rec.addContribution(linefit)\n", + "rec.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " FitContributions may have many parameters. We need to tell the recipe\n", + " which of them should be controlled and potentially optimized in a fit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rec.addVar(rec.linefit.A)\n", + "rec.addVar(rec.linefit.B)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " The call of the addVar function also created two attributes A and B for the rec object,\n", + " which link to the A and B parameters of the linefit contribution.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"rec.A =\", rec.A)\n", + "print(\"rec.A.value =\", rec.A.value)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " The names of the declared variables are stored in the `rec.names` attribute\n", + " and the corresponding values in `rec.values`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"rec.values =\", rec.values)\n", + "print(\"rec.names =\", rec.names)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " Finally the recipe objects provides a residual() function to calculate\n", + " the difference between the observed and simulated values. The residual\n", + " function can accept a list of new variable values in the same order as\n", + " rec.names." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"rec.residual() =\", rec.residual())\n", + "print(\"rec.residual([2, 4]) =\", rec.residual([2, 4]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The FitRecipe.residual function can be directly used with the scipy\n", + "leastsq function for minimizing a sum of squares." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.optimize import leastsq\n", + "leastsq(rec.residual, rec.values)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " Recipe variables and the linked line-function parameters are set to the\n", + " new optimized values." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(rec.names, \"-->\", rec.values)\n", + "linefit.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " The calculated function is available in the `ycalc` attribute of the profile.\n", + " It can be also accessed from the `linefit` contribution attribute of the\n", + " recipe as `rec.linefit.profile.ycalc`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot(linedata.x, linedata.y, 'x', linedata.x, linedata.ycalc, '-')\n", + "title('Line fit using the leastsq least-squares optimizer');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The `FitRecipe.scalarResidual()` function returns the sum of squares and can\n", + "be used with a minimizer that requires a scalar function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.optimize import fmin\n", + "fmin(rec.scalarResidual, [1, 1])\n", + "print(rec.names, \"-->\", rec.values)\n", + "plot(linedata.x, linedata.y, 'x', linedata.x, linedata.ycalc, '-')\n", + "title('Line fit using the fmin scalar optimizer');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For a converged fit recipe, the details of the fit can be extracted\n", + " with the FitResults class." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from diffpy.srfit.fitbase import FitResults\n", + "res = FitResults(rec)\n", + "print(res)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Variables defined in the recipe can be fixed to a constant value." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rec.fix(B=0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The fixed variables can be checked using the \"fixednames\" and\n", + " \"fixedvalues\" attributes of a recipe." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"free:\", rec.names, \"-->\", rec.names)\n", + "print(\"fixed:\", rec.fixednames, \"-->\", rec.fixedvalues)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The fit can be rerun with a constant variable B." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "leastsq(rec.residual, rec.values)\n", + "print(FitResults(rec))\n", + "plot(linedata.x, linedata.y, 'x', linedata.x, linedata.ycalc, '-')\n", + "title('Line fit for variable B fixed to B=0');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Fixed variables may be released with the `free()` function.\n", + " Calling it as `free(\"all\")` releases all fixed variables." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rec.free('all')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Variables may be constrained to a result of an expression." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rec.constrain(rec.A, \"2 * B\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Perform linear fit where the slope value must be two times the offset." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "leastsq(rec.residual, rec.values)\n", + "print(FitResults(rec))\n", + "plot(linedata.x, linedata.y, 'x', linedata.x, linedata.ycalc, '-')\n", + "title('Line fit for variable A constrained to A = 2*B');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Constraint expressions can be removed by calling the unconstrain function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rec.unconstrain(rec.A)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " Variables may be restrained to a specific range. Here \"ub\" is the upper\n", + " boundary and \"sig\" acts as a standard deviation for ((x - ub)/sig)**2\n", + " penalty function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "arst = rec.restrain(rec.A, ub=0.2, sig=0.001)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Perform fit with the line slope restrained to a maximum value of 0.2:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "leastsq(rec.residual, rec.values)\n", + "print(FitResults(rec))\n", + "plot(linedata.x, linedata.y, 'x', linedata.x, linedata.ycalc, '-')\n", + "title('Line fit with A restrained to an upper bound of 0.2');" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "cmi", + "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.13.5" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/docs/examples/core/linefit/SrFit-demo-constraints-restraints.py b/docs/examples/core/linefit/SrFit-demo-constraints-restraints.py new file mode 100644 index 0000000..1a4a298 --- /dev/null +++ b/docs/examples/core/linefit/SrFit-demo-constraints-restraints.py @@ -0,0 +1,142 @@ +from __future__ import print_function + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.pyplot import plot, title +from scipy.optimize import fmin, leastsq + +from diffpy.srfit.fitbase import ( + FitContribution, + FitRecipe, + FitResults, + Profile, +) + + +def main(): + # ---------------------------------------------------------------------- + # Generate synthetic noisy data: y = 0.5 * x + 3 + noise + # ---------------------------------------------------------------------- + xobs = np.arange(-10, 10.1) + dyobs = 0.3 * np.ones_like(xobs) + yobs = 0.5 * xobs + 3 + dyobs * np.random.randn(xobs.size) + + plot(xobs, yobs, "x") + title("y = 0.5*x + 3 with Gaussian noise (σ=0.3)") + plt.show() + # ---------------------------------------------------------------------- + # Create a Profile object to hold the data + # ---------------------------------------------------------------------- + linedata = Profile() + linedata.setObservedProfile(xobs, yobs, dyobs) + + # ---------------------------------------------------------------------- + # Define a FitContribution: linear model A*x + B + # ---------------------------------------------------------------------- + linefit = FitContribution("linefit") + linefit.setProfile(linedata) + linefit.setEquation("A * x + B") + + linefit.show() + + # Assign initial guesses for parameters + linefit.A = 3 + linefit.B = 5 + print("Initial A:", linefit.A, "value:", linefit.A.value) + print("Initial B:", linefit.B, "value:", linefit.B.value) + + # Evaluate model with initial parameters + print("linefit.evaluate() =", linefit.evaluate()) + print("linefit.residual() =", linefit.residual()) + + plot(xobs, yobs, "x", linedata.x, linefit.evaluate(), "-") + title("Line simulated at A=3, B=5") + plt.show() + + # ---------------------------------------------------------------------- + # Create a FitRecipe to manage fitting + # ---------------------------------------------------------------------- + rec = FitRecipe() + rec.clearFitHooks() + rec.addContribution(linefit) + rec.show() + + # Add variables to be refined + rec.addVar(rec.linefit.A) + rec.addVar(rec.linefit.B) + + print("rec.A =", rec.A) + print("rec.A.value =", rec.A.value) + print("rec.values =", rec.values) + print("rec.names =", rec.names) + print("rec.residual() =", rec.residual()) + print("rec.residual([2, 4]) =", rec.residual([2, 4])) + + # ---------------------------------------------------------------------- + # Fit using least squares optimizer + # ---------------------------------------------------------------------- + leastsq(rec.residual, rec.values) + print("After leastsq:", rec.names, "-->", rec.values) + linefit.show() + + plot(linedata.x, linedata.y, "x", linedata.x, linedata.ycalc, "-") + title("Line fit using leastsq optimizer") + plt.show() + + # ---------------------------------------------------------------------- + # Fit using scalar optimizer (fmin) + # ---------------------------------------------------------------------- + fmin(rec.scalarResidual, [1, 1]) + print("After fmin:", rec.names, "-->", rec.values) + + plot(linedata.x, linedata.y, "x", linedata.x, linedata.ycalc, "-") + title("Line fit using fmin optimizer") + plt.show() + + # Display fit results + res = FitResults(rec) + print(res) + + # ---------------------------------------------------------------------- + # Example: Fixing a parameter + # ---------------------------------------------------------------------- + rec.fix(B=0) + print("Free:", rec.names, "-->", rec.values) + print("Fixed:", rec.fixednames, "-->", rec.fixedvalues) + + leastsq(rec.residual, rec.values) + print("Fit with B fixed to 0:", FitResults(rec)) + + plot(linedata.x, linedata.y, "x", linedata.x, linedata.ycalc, "-") + title("Line fit with B fixed at 0") + plt.show() + + rec.free("all") + + # ---------------------------------------------------------------------- + # Example: Adding a constraint (A = 2*B) + # ---------------------------------------------------------------------- + rec.constrain(rec.A, "2 * B") + leastsq(rec.residual, rec.values) + print("Fit with A constrained to 2*B:", FitResults(rec)) + + plot(linedata.x, linedata.y, "x", linedata.x, linedata.ycalc, "-") + title("Line fit with constraint A=2*B") + plt.show() + + rec.unconstrain(rec.A) + + # ---------------------------------------------------------------------- + # Example: Adding a restraint (A close to <= 0.2 with penalty) + # ---------------------------------------------------------------------- + rec.restrain(rec.A, ub=0.2, sig=0.001) + leastsq(rec.residual, rec.values) + print("Fit with A restrained to ub=0.2:", FitResults(rec)) + + plot(linedata.x, linedata.y, "x", linedata.x, linedata.ycalc, "-") + title("Line fit with restraint on A (ub=0.2)") + plt.show() + + +if __name__ == "__main__": + main() diff --git a/news/core-ex.rst b/news/core-ex.rst new file mode 100644 index 0000000..07232e0 --- /dev/null +++ b/news/core-ex.rst @@ -0,0 +1,23 @@ +**Added:** + +* Add line fitting example. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +*