{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Effect of measurement noise on accuracy of fit\n\nExperimental noise often corrupts concentration profiles. When the diffusion\nmatrix is determined from noisy concentration profiles, the accuracy of the\nestimation decreases when the intensity of the noise increases. Furthermore,\nminor eigenvectors and eigenvalues are more affected by the noise than the\ndominant eigenvector.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nfrom multidiff import compute_diffusion_matrix, create_diffusion_profiles"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Our diffusion system has three components (e.g. oxides). The diffusion\nmatrix has two eigenvalues with quite different magnitudes. Exchange vectors\ncorrespond to the exchange of the three possible pairs of components\n(meaning that each endmember is enriched in one component and poorer in\nanother, compared to the other endmember).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "diags = np.array([0.5, 5.])\nP = np.matrix([[1, 1], [-1, 0]])\nxpoints_exp1 = np.linspace(-10, 10, 100)\nx_points = [xpoints_exp1] * 3\n\n\nexchange_vectors = np.array([[0, 1, 1],\n                             [1, -1, 0],\n                             [-1, 0, -1]])\n\n\n# Initialization for the diffusion matrix\ndiags_init = np.array([1, 2])\nP_init = np.matrix([[1.5, 0.9], [-1, -0.7]])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now iterate over increasing values of measurement Gaussian noise. For\neach value of the noise, synthetic diffusion profiles are created, and the\ndiffusion matrix is computed from the profiles. We plot the relative error\non the eigenvalues, defined as\nabs((lambda_computed - lambda_real)/lambda_real).\nThe relative error is consistently larger for the smaller eigenvalue, since\nerrors on the larger (dominant) eigenvalue and eigenvectors result in errors\nof comparable magnitude on the smaller eigenvalue(s) and eigenvector(s).\n\nNote that for the dominant eigenvalue, the relative error is smaller than the\nintensity of the noise. Indeed, we have enough experiments to overconstrain\nthe system, plus a large number of measurement points.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "noises = np.arange(0, 0.2, 0.02)\nn_real = 50\nerrors = np.empty((len(noises), n_real, 2))\n\n\nfor i, noise in enumerate(noises):\n    for j in range(n_real):\n        concentration_profiles = create_diffusion_profiles((diags, P), \n                                    x_points,\n                                    exchange_vectors, noise_level=noise,\n                                    seed=i * n_real + j)\n        diags_res, eigvecs, _, _, _ = (\n            compute_diffusion_matrix((diags_init, P_init), x_points,\n                                        concentration_profiles, plot=False))\n        errors[i, j, :] = ((np.sort(diags_res) - diags) / diags)**2.\n\n\nerrors = np.sqrt(errors.mean(axis=1))\n\nplt.figure()\nplt.plot(noises, errors[:, 0], 'co--',\n        label='relative error on minor eigenvalue')\nplt.plot(noises, errors[:, 1], 'rs-', ms=10,\n         label='relative error on major eigenvalue')\nplt.xlabel('intensity of Gaussian noise', fontsize=16)\nplt.title('Relative error on fitted eigenvalues')\nplt.legend(loc='best')\n\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.8.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}