Modeler-Copy1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy.interpolate import interp1d\n",
    "import scipy.optimize as so\n",
    "import tmm_core as tmm\n",
    "import lmfit\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Material:\n",
    "    def __init__(self, path, name, wavelengths):\n",
    "        f = pd.read_csv(path)\n",
    "        wv_raw = np.array(f.wv)\n",
    "        nc = np.array(f.n+f.k*1j)  # complex refractive index\n",
    "        self.name = name\n",
    "        self.f = interp1d(wv_raw, nc, kind='cubic')\n",
    "        self.min_wv = min(wv_raw)\n",
    "        self.max_wv = max(wv_raw)\n",
    "        self.df = pd.Series(data=nc, index=wv_raw, name=name)\n",
    "        # self.df = self.df.reindex(wavelengths)\n",
    "        # self.df = self.df.interpolate('spline', order=3)\n",
    "        self.interp(wavelengths)\n",
    "\n",
    "    def interp(self, wavelengths):\n",
    "        wavelengths = wavelengths[np.nonzero(np.logical_and(wavelengths > self.min_wv, wavelengths < self.max_wv))]\n",
    "        nc = self.f(wavelengths)\n",
    "        self.df = pd.DataFrame(nc,wavelengths, [self.name])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Model:\n",
    "    def __init__(self):\n",
    "        self.wavelength = np.arange(300, 1700, 1)\n",
    "        self.increment = .2\n",
    "        self.index_array = np.array([])\n",
    "        self.materials = []\n",
    "        self.mat_df = pd.DataFrame([], self.wavelength)\n",
    "        self.data = {}\n",
    "        self.thicknesses = [sp.inf, 0, sp.inf]\n",
    "        self.layers = []\n",
    "\n",
    "        self.add_material(\"TCO\", './Materials/Semiconductor/TCO.csv')\n",
    "        self.add_material(\"CdSe\", './Materials/Semiconductor/CdSe.csv')\n",
    "        self.add_material(\"CdTe\", './Materials/Semiconductor/CdTe.csv')\n",
    "        self.add_material(\"Air\", './Materials/Semiconductor/Air.csv')\n",
    "        self.add_material(\"BK7\", './Materials/Dielectric/BK7.csv')\n",
    "        self.add_material(\"SLG\", './Materials/Dielectric/SLG.csv')\n",
    "        self.add_material(\"SS TCO\", './Materials/Dielectric/Sunnyside TCO.csv')\n",
    "        self.add_material(\"SiO2\", './Materials/Dielectric/SiO2-Jaw.csv')\n",
    "        \n",
    "    def add_material(self, film, path):\n",
    "        if film not in self.mat_df:\n",
    "            mat = Material(path, film, self.wavelength)\n",
    "            self.materials.append(mat)\n",
    "            self.mat_df = self.mat_df.join(mat.df)\n",
    "\n",
    "    def set_wavelength(self, low, high, interval):\n",
    "        self.wavelength = np.arange(low, high, interval)\n",
    "        df = self.mat_df.reindex(self.wavelength)\n",
    "        df = df.interpolate('spline', order=3)\n",
    "        self.mat_df = df\n",
    "\n",
    "    def better_bruggeman(self, n1, n2, percent_included):\n",
    "        p = n1/n2\n",
    "        b = .25*((3*percent_included-1)*(1/p-p)+p)\n",
    "        z = b + (b**2 + .5)**0.5\n",
    "        e = z*n1*n2\n",
    "        return {\"e\": e, \"n\": e**0.5, 'conc': percent_included, \"n1\": n1, 'n2': n2}\n",
    "\n",
    "    def brug_transform(self, df, layer, incl, percent):\n",
    "        p = df[layer]/incl\n",
    "        b = .25*((3*percent-1)*(1/p-p)+p)\n",
    "        z = b + (b**2 + .5)**0.5\n",
    "        e = z*df[layer]*incl\n",
    "        n = e**.5\n",
    "        df[layer] = n\n",
    "    \n",
    "    def run(self, wavelengths, void_percent):\n",
    "        mat = self.mat_df.ix[wavelengths]\n",
    "        mat = mat[self.layers]\n",
    "        self.brug_transform(mat, self.layers[1], mat['Air'], void_percent)\n",
    "        self.index_array = np.array(mat)\n",
    "        theta0 = 45*sp.pi/180\n",
    "        self.data = tmm.unpolarized_RT(self.index_array, self.thicknesses, theta0, wavelengths)\n",
    "        \n",
    "    def normalized(a, axis=-1, order=2):\n",
    "        l2 = np.atleast_1d(np.linalg.norm(a, order, axis))\n",
    "        l2[l2==0] = 1\n",
    "        return a / np.expand_dims(l2, axis)\n",
    "    \n",
    "    def get_R(self, wavelengths, thickness, void_percent):\n",
    "        self.thicknesses[1] = thickness\n",
    "        self.run(wavelengths, void_percent) \n",
    "        R = self.data['R']\n",
    "        R -= min(R)\n",
    "        R /= max(R)\n",
    "        return R\n",
    "        \n",
    "    def RMSE(self, thickness, data):\n",
    "        df = pd.DataFrame(data)\n",
    "        self.get_R(thickness)\n",
    "        model = pd.DataFrame(self.data['R'], index=self.wavelength)\n",
    "        df = df.join(model, how='inner')\n",
    "        n = len(df.index)\n",
    "        return (sum((data-model)**2)/n)**0.5\n",
    "    \n",
    "    def norm(self, wavelength):\n",
    "        df = self.df.ix[wavelength]\n",
    "        df = df - df.min()\n",
    "        df = df / df.max()\n",
    "        return df\n",
    "    \n",
    "    def fit(self, wv, data):\n",
    "        mod = lmfit.Model(self.get_R, ['wavelengths'], ['thickness','void_percent'])\n",
    "        mod.set_param_hint('thickness', value = 120, min=50, max=300)\n",
    "        mod.set_param_hint('void_percent', value = .15, min=0, max=1)\n",
    "        \n",
    "        weight = np.array(data.weight.ix[wv])\n",
    "        R = data.series.ix[wv]\n",
    "        result = mod.fit(R, wavelengths=wv, weights=weight)\n",
    "        \n",
    "        RMSE = (sp.sum(result.residual**2)/(result.residual.size-2))**0.5\n",
    "        bf_values = result.best_values\n",
    "        bf_str = 'thk: ' + str(round(bf_values['thickness'])) +\", Void %: \" + str(round(bf_values['void_percent']*100, 2))\n",
    "        txt_spot = wv.min()-200 + (wv.max()-wv.min()) / 2\n",
    "\n",
    "        fig = plt.figure()\n",
    "        ax = fig.add_subplot(111)\n",
    "        ax.text(txt_spot, .9, \"RMSE: \"+str(round(RMSE, 3)))\n",
    "        ax.text(txt_spot, .85, bf_str)\n",
    "        result.plot_fit(yerr=np.zeros(len(data.series.index)), data_kws ={'marker':'+'})\n",
    "\n",
    "        plt.show()\n",
    "        print(result.fit_report())\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Data:\n",
    "    def __init__(self, path, name, wavelengths, c_type):\n",
    "        f = pd.read_csv(path)\n",
    "        wv_raw = np.array(f.wv)\n",
    "        R = np.array(f.R)\n",
    "        self.name = name\n",
    "        self.f = interp1d(wv_raw, R, kind='cubic')\n",
    "        self.min_wv = min(wv_raw)\n",
    "        self.max_wv = max(wv_raw)\n",
    "        self.raw_series = pd.Series(data=R, index=wv_raw, name=name)\n",
    "        self.interp(wavelengths)\n",
    "        self.series = self.norm(wavelengths)\n",
    "        self.weight_amp = 100\n",
    "        self.weight_wid = 30\n",
    "        self.weight_cen = 597\n",
    "        self.weight = pd.Series(1, self.series.index, name='Weighting')\n",
    "        \n",
    "    def interp(self, wavelengths):\n",
    "        wavelengths = wavelengths[np.nonzero(np.logical_and(wavelengths > self.min_wv, wavelengths < self.max_wv))]\n",
    "        data = self.f(wavelengths)\n",
    "        self.raw_series = pd.Series(data=data, index=wavelengths, name=self.name)\n",
    "    \n",
    "    def norm(self, wavelength, data_filter=\"Savitzky-Golay\"):\n",
    "        s = self.raw_series.ix[wavelength]\n",
    "        s = s.dropna()\n",
    "        s = s - s.min()\n",
    "        s = s / s.max()\n",
    "        if data_filter == \"Savitzky-Golay\":\n",
    "            arr = self.savitzky_golay(y=np.array(s), window_size=11, order=3, deriv=0, rate=1)\n",
    "            s = pd.Series(arr, s.index, name=self.raw_series.name+\" processed\")\n",
    "        return s\n",
    "    \n",
    "    def savitzky_golay(self, y, window_size, order, deriv=0, rate=1):\n",
    "        r\"\"\"Smooth (and optionally differentiate) data with a Savitzky-Golay filter.\n",
    "        The Savitzky-Golay filter removes high frequency noise from data.\n",
    "        It has the advantage of preserving the original shape and\n",
    "        features of the signal better than other types of filtering\n",
    "        approaches, such as moving averages techniques.\n",
    "        Parameters\n",
    "        ----------\n",
    "        y : array_like, shape (N,)\n",
    "            the values of the time history of the signal.\n",
    "        window_size : int\n",
    "            the length of the window. Must be an odd integer number.\n",
    "        order : int\n",
    "            the order of the polynomial used in the filtering.\n",
    "            Must be less then `window_size` - 1.\n",
    "        deriv: int\n",
    "            the order of the derivative to compute (default = 0 means only smoothing)\n",
    "        rate: int\n",
    "            don't know\n",
    "        Returns\n",
    "        -------\n",
    "        ys : ndarray, shape (N)\n",
    "            the smoothed signal (or it's n-th derivative).\n",
    "        Notes\n",
    "        -----\n",
    "        The Savitzky-Golay is a type of low-pass filter, particularly\n",
    "        suited for smoothing noisy data. The main idea behind this\n",
    "        approach is to make for each point a least-square fit with a\n",
    "        polynomial of high order over a odd-sized window centered at\n",
    "        the point.\n",
    "        Examples\n",
    "        --------\n",
    "        t = np.linspace(-4, 4, 500)\n",
    "        y = np.exp( -t**2 ) + np.random.normal(0, 0.05, t.shape)\n",
    "        ysg = savitzky_golay(y, window_size=31, order=4)\n",
    "        import matplotlib.pyplot as plt\n",
    "        plt.plot(t, y, label='Noisy signal')\n",
    "        plt.plot(t, np.exp(-t**2), 'k', lw=1.5, label='Original signal')\n",
    "        plt.plot(t, ysg, 'r', label='Filtered signal')\n",
    "        plt.legend()\n",
    "        plt.show()\n",
    "        References\n",
    "        ----------\n",
    "        .. [1] A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of\n",
    "           Data by Simplified Least Squares Procedures. Analytical\n",
    "           Chemistry, 1964, 36 (8), pp 1627-1639.\n",
    "        .. [2] Numerical Recipes 3rd Edition: The Art of Scientific Computing\n",
    "           W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery\n",
    "           Cambridge University Press ISBN-13: 9780521880688\n",
    "        \"\"\"\n",
    "        import numpy as np\n",
    "        from math import factorial\n",
    "\n",
    "        try:\n",
    "            window_size = np.abs(np.int(window_size))\n",
    "            order = np.abs(np.int(order))\n",
    "        except (ValueError, msg):\n",
    "            raise ValueError(\"window_size and order have to be of type int\")\n",
    "        if window_size % 2 != 1 or window_size < 1:\n",
    "            raise TypeError(\"window_size size must be a positive odd number\")\n",
    "        if window_size < order + 2:\n",
    "            raise TypeError(\"window_size is too small for the polynomials order\")\n",
    "        order_range = range(order+1)\n",
    "        half_window = (window_size -1) // 2\n",
    "        # precompute coefficients\n",
    "        b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])\n",
    "        m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)\n",
    "        # pad the signal at the extremes with\n",
    "        # values taken from the signal itself\n",
    "        firstvals = y[0] - np.abs(y[1:half_window+1][::-1] - y[0])\n",
    "        lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])\n",
    "        y = np.concatenate((firstvals, y, lastvals))\n",
    "        return np.convolve(m[::-1], y, mode='valid')\n",
    "    \n",
    "    def get_weighting(self, cen, amp, wid):\n",
    "        x = np.array(self.series.index)\n",
    "        weight = self.weight\n",
    "        for i in cen:\n",
    "            g = amp * np.exp(-(x-i)**2 /wid)\n",
    "            weight += g\n",
    "        return pd.Series(weight, x, name=\"Weighting\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = Model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = Data('C:/WORK!!!!/VF088 - ARC Thickness/Golden Coupon Data/0996 G4 XLS.csv', 'Data', np.arange(300,1700), 'R')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data.weight = pd.Series(1, data.series.index, name='Weighting')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data.weight_cen = [404]\n",
    "data.weight_wid = 30\n",
    "data.weight_amp = 100\n",
    "data.weight = data.get_weighting(data.weight_cen, data.weight_amp, data.weight_wid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data.weight[602]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " - Adding parameter for hint \"thickness\"\n",
      " - Adding parameter for hint \"void_percent\"\n",
      " - Adding parameter \"thickness\"\n",
      " - Adding parameter \"void_percent\"\n",
      "[[Model]]\n",
      "    Model(get_R)\n",
      "[[Fit Statistics]]\n",
      "    # function evals   = 112\n",
      "    # data points      = 350\n",
      "    # variables        = 2\n",
      "    chi-square         = 1.458\n",
      "    reduced chi-square = 0.004\n",
      "[[Variables]]\n",
      "    thickness:      152.000744 +/- 0.884159 (0.58%) (init= 120)\n",
      "    void_percent:   2.3138e-10 +/- 0.013720 (5929638297.50%) (init= 0.15)\n",
      "[[Correlations]] (unreported correlations are <  0.100)\n",
      "    C(thickness, void_percent)   =  0.997 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "model.layers = ['Air', 'SiO2', 'SS TCO']\n",
    "wv = np.arange(350, 700)\n",
    "data.series = data.norm(wv)\n",
    "model.fit(wv, data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model.layers = ['Air', 'SiO2', 'SS TCO']\n",
    "mat = model.mat_df[model.layers]\n",
    "model.thicknesses = [sp.inf, 120, sp.inf]\n",
    "model.index_array = np.array(mat)\n",
    "wv = np.array(model.wavelength)\n",
    "theta0 = 55*sp.pi/180\n",
    "mod_data = tmm.unpolarized_RT(model.index_array, model.thicknesses, theta0, wv)\n",
    "s = tmm.coh_tmm('s', model.index_array, model.thicknesses, theta0, wv)\n",
    "p = tmm.coh_tmm('p', model.index_array, model.thicknesses, theta0, wv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pd.DataFrame(np.array([mod_data['R'], s['R'], p['R']]).T, index=wv, columns=['unpolarized', 's', 'p'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pd.DataFrame(np.array([mat[\"Air\"].real, mat['Air'].imag, mat[\"SiO2\"].real, mat['SiO2'].imag, mat[\"SS TCO\"].real, mat['SS TCO'].imag]).T, mat.index, columns=['n-Air', 'k-Air', 'n-SiO2', 'k-SiO2', 'n-TCO', 'k-TCO']).to_csv(\"Test.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}