Better structure of data analysis code, including code to generate

all used images and information for my data.

.gitignore

@@ -16,7 +16,7 @@ MATLAB/*.asv
 MATLAB/**/*.asv
 
 # Python generated files
-py/__pycache__
+py/**/__pycache__
 
 # Automatic generated file on Windows
 desktop.ini

py/datafit.py

-"""
-In this file, the four measurements from Julian are analysed and temperature
-dependent models fitted onto these results.
-"""
-
-import numpy as np
-import matplotlib.pyplot as plt
-import scipy.optimize
-from scripts import linreg, chisq_ndf
-from measurementclasses import Experiment
-
-# As their layer was really thin, assume 2D theories
-dimensions = 2;
-
-# Literature values for melting temperature, latent heat and interfacial energy
-class Literature:
-    Tm = 817
-    L = 878.9e6;
-    IE = 5.5e-2 *3
-
-# Some physial constants
-kB = 1.380649e-23
-eV = 1.6021766208e-19
-zeroC = 273.15
-
-# Function for bulk free energy
-dGv = lambda T,L,Tm: L*(Tm-T)/Tm*2*T/(Tm+T)
-
-
-# The temperatures and geometries the data was taken at
-Ts = np.array([115,120,125,130]) + zeroC
-height = 30e-9
-
-# Read all the data
-A = Experiment("../data/experiment/Sb2Te_asd_20W_115C.csv", height)
-B = Experiment("../data/experiment/Sb2Te_asd_20W_120C.csv", height)
-C = Experiment("../data/experiment/Sb2Te_asd_20W_125C.csv", height)
-D = Experiment("../data/experiment/Sb2Te_asd_20W_130C.csv", height)
-
-
-# Get all JMAK fits and growth velocities and calculate nucleation rate from
-# them; get all nucleation rates and calculate nucleation rates from them too
-Js, vs, uvs = np.ones(Ts.shape), np.ones(Ts.shape), np.ones(Ts.shape)
-for i,M in enumerate([A,B,C,D]):
-    v, uv = M.all_rfits(show=False)
-    vs[i] = v
-    uvs[i] = uv
-    [(dNdt,_)] = M.nucleationratefit(show=False)
-    if dimensions == 2:
-        J = dNdt / M.A
-    elif dimensions == 3:
-        J = dNdt / (M.A * height)
-    else:
-        print("Dimensions need to be 2 or 3.")
-    Js[i] = J
-
-
-def velocityfit(show=True):
-    """ Plot velocities over temperature and fit arrhenius behaviour. 
-    Argument show can be 2 to show the arrhenius plot. """
-    fitfunc = lambda x,a,b: a*np.exp(-b/x)
-    m0, _, c0, _, _, _ = linreg(1/Ts, np.log(vs), uvs)
-    popt, pcov = scipy.optimize.curve_fit(fitfunc, Ts, vs, p0=[np.exp(c0),-m0],
-                                          sigma=uvs, absolute_sigma=True)
-    if show==2:
-        plt.errorbar(1/Ts, np.log(vs), uvs/vs, linestyle="none", marker="_",
-                     label="data")
-        m, c = -popt[1], np.log(popt[0])
-        plt.plot(1/Ts, m/Ts+c, label="fit")
-        plt.xlabel("$1/(T \\ [{}^{\circ}C])$")
-        plt.ylabel("$\\log(v \\ [m/s])$")
-    elif show:
-        plt.errorbar(Ts, vs, uvs, linestyle="none", marker="_", label="data")
-        pxs = np.linspace(np.min(Ts)-1,np.max(Ts)+1,100)
-        plt.plot(pxs,fitfunc(pxs,*popt), label="fit")
-        plt.xlabel("$T \\ [{}^{\circ}C]$")
-        plt.ylabel("$v \\ [m/s]$")
-    if show:
-        print("Fitting v = a*exp(-b/T)")
-        print("a:",popt[0],"+-",np.sqrt(pcov[0,0]))
-        print("b:",popt[1],"+-",np.sqrt(pcov[1,1]))
-        Ea, uEa = [x/eV*kB for x in [popt[1], np.sqrt(pcov[1,1])]]
-        print("   = (",Ea,"+-",uEa,")eV/kB")
-        print("chisq/ndf:", chisq_ndf(Ts,vs,uvs,fitfunc,popt))
-        plt.legend()
-        plt.tight_layout()
-        plt.show()
-    return popt, pcov
-
-def nucleationratefit(plotrange=None, show=True):
-    """ 
-    Fit classical nucleation theory formula to J(T).
-    Use either JMAK theory fit  results(second argument False) or counted
-    nuclei (second argument True).
-    A range of temperatures to plot the calculated rate over can be
-    given manually as a 2 tuple.
-    """
-    
-    # Formulas for critical energy and nucleation rate
-    dGc = lambda T,IE,L,Tm: 16*np.pi*IE**3/(3*dGv(T,L,Tm)**2)
-    J = lambda T,b,IE,L,Tm : (b*dGv(T,L,Tm)**2/(8*np.pi) *
-                              1/np.sqrt(IE**3*kB*T) * 
-                              np.exp(-dGc(T,IE,L,Tm)/kB/T))
-    # Fit with a free prefactor and interfacial energy, other should be known
-    fitfunc = lambda T,b,IE : J(T,b,IE,Literature.L,Literature.Tm)
-#    fitfunc = lambda T,b,IE,L,Tm : J(T,b,IE,L,Tm)
-    
-    # Guess for IE is literature value, choose initial prefactor to fit some
-    # arbitary point. Do the fit.
-    b0 = Js[3]/fitfunc(Ts[3],1,Literature.IE)
-    popt, pcov = scipy.optimize.curve_fit(fitfunc, Ts, Js, maxfev=60000,
-                                          p0=[b0,Literature.IE])
-#                              p0=[b0,Literature.IE,Literature.L,Literature.Tm]
-    if show:
-        print("Fitting J = b*dGv**2/(8*pi)/sqrt(IE**3*kB*T) * exp(-dGc/kB/T)")
-        print("        dGc = 16*pi*IE^3/(3*dGv)")
-        print("        dGv = L*(Tm-T)/Tm*2*T/(Tm+T)")
-        print("b: ",popt[0],"+-",np.sqrt(pcov[0,0]))
-        print("IE:",popt[1],"+-",np.sqrt(pcov[1,1]))
-#        print("L: ",popt[2]," +- ",np.sqrt(pcov[2,2]))
-#        print("Tm:",popt[3]," +- ",np.sqrt(pcov[3,3]))
-        plt.plot(Ts, Js, linestyle="none", marker="x", label="mean rates")
-        if plotrange is None:
-            pxs = np.linspace(np.min(Ts)-3,np.max(Ts)+3,50)
-        else:
-            pxs = np.linspace(plotrange[0],plotrange[-1],200)
-        plt.plot(pxs,fitfunc(pxs,*popt), label="fit")
-        plt.xlabel("$T \\ [{}^{\circ}C]$")
-        plt.ylabel("$J \\ [1/m^{}/s]$".format(dimensions))
-        plt.tight_layout()
-        plt.legend()
-        plt.show()
-    return popt, pcov
-
-def mobilityfit(show=True):
-    """ Calculate mobility (as defined in the AIST paper), plot over
-    temperature and fit arrhenius behaviour """
-    
-    M = v/dGv(Ts,Literature.L,Literature.Tm)
-    
-    plt.plot(Ts, M, linestyle="none", marker="x", label="data")
-    expfunc = lambda x,a,b: a*np.exp(b/x)
-    m0, _, c0, _, _, _ = linreg(1/Ts, np.log(M), uvs)
-    popt, pcov = scipy.optimize.curve_fit(expfunc, Ts, M, p0=[np.exp(c0),m0],
-                                          sigma=uvs)
-    pxs = np.linspace(np.min(Ts)-3,np.max(Ts)+3,50)
-    plt.plot(pxs,expfunc(pxs,np.exp(c0),m0),label="fit")
-    plt.xlabel("T [C]")
-    plt.ylabel("M [m/s / (J/m^3)]")
-
-#A.nucleationratefit(splitat=[1.2e3,10e3,19e3,28e3])
-#B.nucleationratefit(splitat=[1000,8888])
-#C.nucleationratefit(splitat=[2000,3500])
-#D.nucleationratefit(splitat=[250,1100,1900])

py/experimentalfits.py

+import numpy as np
+import matplotlib.pyplot as plt
+import fits
+
+# Set font size and stuff for images in the tex document
+plt.rcParams.update({'font.size': 16,
+                     'lines.linewidth': 2,
+                     'lines.markersize': 10,
+                     'lines.markeredgewidth': 2})
+
+# Read the experiments
+height = 30e-9
+Ex115 = fits.Experiment("../data/experiment/Sb2Te_asd_20W_115C.csv", height)
+Ex120 = fits.Experiment("../data/experiment/Sb2Te_asd_20W_120C.csv", height)
+Ex125 = fits.Experiment("../data/experiment/Sb2Te_asd_20W_125C.csv", height)
+Ex130 = fits.Experiment("../data/experiment/Sb2Te_asd_20W_130C.csv", height)
+
+# Methods performing fits on this data, giving images and information ready to
+# use in tex
+
+def show_incubationtime():
+    axes = Ex120.plot_N()
+    fig = axes.figure
+    line = axes.lines[0]
+    line.set_linestyle(":")
+    axes.set_xlim(-120, 3500)
+    axes.set_ylim(-5, 120)
+    axes.set_title("Sb${}_2$Te, 120°C")
+    fig.canvas.draw()
+    fig.tight_layout()
+
+def fit_v(Ex=Ex130, i=2):
+    vfit = Ex.velocityFit()
+    vfit.fitall()
+    vfit.info()
+    axes = vfit.plot(i)
+    Ex.plot_r(i, axes)
+    [fitline, expline] = axes.lines
+    expline.set_color("red")
+    fitline.set_color("black")
+    axes.legend()
+    fig = axes.figure
+    fig.tight_layout()
+
+def fit_J_115():
+    jfit = Ex115.nucleationRateFit([1.7e3,9.7e3,20e3]).fit()
+    jfit.info()
+    axes = jfit.plotJ()
+    axes.set_xlim(-500,31800)
+    axes.set_ylim(-0.001, 0.048)
+    fig = axes.figure
+    fig.tight_layout()
+
+def fit_Nc_115():
+    jfit = Ex115.nucleationRateFit([1.7e3,9.7e3,20e3]).fit()
+    jfit.info()
+    axes = jfit.plotN()
+    axes.figure.tight_layout()
+
+def comp_120():
+    axes = Ex120.plot_fc()
+    axes.lines[0].set_markersize(6)
+    fit = Ex120.jmakFit().fit(2)
+    rec = Ex120.jmakBuild(2)
+    fit.plot(axes)
+    rec.plot(axes)
+    axes.legend()
+    axes.figure.tight_layout()
+
+def comp_125():
+    axes = Ex125.plot_fc()
+    axes.lines[0].set_markersize(6)
+    fit = Ex125.jmakFit().fit(2)
+    rec = Ex125.jmakBuild(2)
+    fit.plot(axes)
+    rec.plot(axes)
+    axes.legend()
+    axes.figure.tight_layout()
+
+# Combine the measurements and perform fits to their relations
+
+S = fits.SeriesOfMeasurements(np.array([115,120,125,130]) + 273.15,
+                              [Ex115,Ex120,Ex125,Ex130], 2)
+S.load()
+V = fits.ArrheniusFit(S).fit()
+C = fits.SimpleCntFit(S, 296.4e6, 837, 9e-2).fit()
+
+def fit_cnt():
+    C.info()
+    fig, axes = plt.subplots(1,1)
+    axes.set_xlim(384.7,405.3)
+    axes.set_ylim(-1e12,4e13)
+    C.plot(axes)
+    axes.legend()
+    fig.tight_layout()
+
+def fit_arrhenius():
+    V.info()
+    axes = V.plot()
+    axes.legend()
+    axes.figure.tight_layout()
+
+def fit_arrhenius_lin():
+    V.info()
+    axes = V.plot_lin()
+    axes.set_xticks(np.linspace(2.48e-3, 2.58e-3, 3))
+    axes.legend()
+    axes.figure.tight_layout()
+
+def fit_mobility(L=296.4e6, Tm=837):
+    mf = S.mobilityFit(lambda T: L*(Tm-T)/(2*Tm)*2*T/(Tm+T) ).fit()
+    mf.info()
+    axes = mf.plot()
+    axes.legend()
+    axes.figure.tight_layout()

py/fits/__init__.py

+"""
+Analysis of experiments on phase change materials.
+"""
+
+# Import classes for single measurements and fits
+from .measurements import *
+from .measurements import __all__
+
+# Import classes for combining several measurements
+from .series import *
+from .series import __all__ as __all2__
+
+# Correct list of available classes
+__all__.extend(__all2__)
+del __all2__

py/scripts.py → py/fits/scripts.py

-# Some math stuff
+"""
+This file provides some general mathematical methods:
+    linreg - Linear regression
+    weightedmean - Weighted mean value of values with uncertainties
+    chisq_ndf - Chisquared per degree of freedom for given data, function and
+        parameters
+"""
 
 import numpy as np
 
@@ -41,4 +47,13 @@ def chisq_ndf(x,y,uncert_y,f,popt):
     Returns the chisquared per degree of freedom for data x, y with
     uncertainty uncert_y and a function f(x,*popt) with parameters popt.
     """
-    return np.sum(((y-f(x,*popt))/uncert_y)**2) / (len(x) - len(popt))
\ No newline at end of file
+    return np.sum(((y-f(x,*popt))/uncert_y)**2) / (len(x) - len(popt))
+
+def enlargedLinspace(left,right,num=50,factor=1):
+    """
+    Get evenly spaced numbers over the interval from left to right, enlarged by
+    a factor.
+    Default factor is one, default amount of points is 50.
+    """
+    shift = (right-left)*factor
+    return np.linspace(left-shift, right+shift, num)
\ No newline at end of file

py/fits/series.py

+"""
+This file provides classes to fit temperature dependent models to data obtained
+from several measurements at different temperatures.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.optimize
+from .scripts import linreg, chisq_ndf
+
+# Some physial constants
+kB = 1.380649e-23
+eV = 1.6021766208e-19
+
+class SeriesOfMeasurements:
+    """
+    A class for a series of several measurements at different temperatures.
+    """
+    
+    def __init__(self, temperatures, measurements, dimensions):
+        """
+        Initialise a series of measurements (given as a list) at different
+        temperatures (given as a list).
+        Analyse assuming a given number of dimensions.
+        """
+        self.Ts = np.array(temperatures)
+        self.Ms = measurements
+        self.dimensions = dimensions
+    
+    def load(self):
+        """
+        Do all the fits for all the measurements.
+        """
+        self.vs = np.zeros(len(self.Ms))
+        self.uvs  = np.zeros(len(self.Ms))
+        self.Js  = np.zeros(len(self.Ms))
+        self.uJs  = np.zeros(len(self.Ms))
+        for i,M in enumerate(self.Ms):
+            vFit = M.velocityFit().fitall()
+            self.vs[i], self.uvs[i] = vFit.v, vFit.uv
+            volume = M.A * M.h
+            jFit = M.nucleationRateFit().fit()
+            [(J,uJ)] = jFit.Js
+            self.Js[i], self.uJs[i] = J/volume, uJ/volume
+        return self
+    
+    def arrheniusFit(self):
+        """
+        Return an arrhenius behaviour fit for the growth velocities of this
+        series of measurements.
+        """
+        return ArrheniusFit(self)
+    
+    def simpleCntFit(self, L, Tm, IE):
+        """
+        Return an object to fit classical nucleation theory to the rates of
+        this series, using given values for the latent heat L, the melting
+        temperature Tm and an initial guess for the interfacial energy IE.
+        """
+        return SimpleCntFit(self, L, Tm, IE)
+    
+    def mobilityFit(self, dGv):
+        """
+        Return an arrhenius behaviour fit for the mobility v/dGv of this
+        series of measurements, using a given temperature dependent function
+        for the free energy.
+        """
+        return MobilityFit(self, dGv)
+
+
+class ArrheniusFit:
+    """ Class to fit an Arrhenius behaviour to the growth velocity. """
+    
+    def __init__(self, series):
+        """ Build an object of a given series of measurements. """
+        self.Ts = series.Ts
+        self.vs = series.vs
+        self.uvs = series.uvs
+    
+    def linfit(self):
+        """ Fit in the linearised problem. """
+        m, um, c, uc, _, _ = linreg(1/self.Ts, np.log(self.vs), self.uvs)
+        self.vinf, self.uvinf = np.exp(c), np.exp(c)*uc
+        self.Ea, self.uEa = -m, um
+        return self
+    
+    def fit(self):
+        """ Perform a fit. """
+        self.linfit()
+        fitfunc = lambda x,a,b: a*np.exp(-b/x)
+        Ts, vs, uvs = self.Ts, self.vs, self.uvs
+        p0 = [self.vinf, self.Ea]
+        popt, pcov = scipy.optimize.curve_fit(fitfunc, Ts, vs, p0=p0,
+                                              sigma=uvs, absolute_sigma=True)
+        self.vinf, self.Ea = popt[0], popt[1]
+        self.uvinf, self.uEa = np.sqrt(pcov[0,0]), np.sqrt(pcov[1,1])
+        return self
+    
+    def v(self,T):
+        """
+        Velocity as a function of temperature according to the fit result.
+        """
+        return self.vinf * np.exp(-self.Ea/T)
+    
+    def info(self):
+        """ Print some information. """
+        print("Fitting v = vinf*exp(-Ea/T)")
+        print("vinf = {} +- {}".format(self.vinf, self.uvinf))
+        print("Ea   = {} +- {}".format(self.Ea, self.uEa))
+        Ea, uEa = [x/eV*kB for x in [self.Ea, self.uEa]]
+        print("   = (",Ea,"+-",uEa,")eV/kB")
+        print("chisq/ndf:", chisq_ndf(self.Ts, self.vs, self.uvs,
+                                      (lambda x,a,b: a*np.exp(-b/x)),
+                                      [self.vinf, self.Ea]))
+        return self
+    
+    def plot_lin(self, axes=None):
+        """ Plot the linearised problem. """
+        if axes is None:
+            fig, axes = plt.subplots(1,1)
+        invTs = 1/self.Ts
+        m, c = -self.Ea, np.log(self.vinf)
+        axes.errorbar(invTs, np.log(self.vs), self.uvs/self.vs,
+                      linestyle="none", marker="_", label="data")
+        axes.plot(invTs, m*invTs+c, label="fit")
+        axes.set_xlabel("$1/(T \\ [K])$")
+        axes.set_ylabel("$\\log(v \\ [m/s])$")
+        return axes
+    
+    def plot(self, axes=None):
+        """ Show the data and fit. """
+        if axes is None:
+            fig, axes = plt.subplots(1,1)
+        axes.errorbar(self.Ts, self.vs, self.uvs, linestyle="none", marker="_",
+                     label="data")
+        pxs = np.linspace(np.min(self.Ts)-1,np.max(self.Ts)+1,100)
+        axes.plot(pxs,self.v(pxs), label="fit")
+        axes.set_xlabel("$T \\ [K]$")
+        axes.set_ylabel("$v \\ [m/s]$")
+        return axes
+
+
+# Function for bulk free energy and critical energy
+def bulkFreeEnergy(T,L,Tm):
+    return L*(Tm-T)/Tm*2*T/(Tm+T)
+
+def criticalFreeEnergy(T,IE,L,Tm):
+    return 16*np.pi*IE**3/(3*bulkFreeEnergy(T,L,Tm)**2)
+
+class SimpleCntFit:
+    """
+    A simple fit of classical nucleation theory onto the nucleation rates,
+    using given values for latent heat L and melting temperature Tm and leaving
+    the prefactor b and interfacial energy IE as fit parameters.
+    """
+    
+    def __init__(self, series, L, Tm, IE):
+        """
+        Build an object of a given series of measurements, using given values
+        for the latent heat L, melting temperature Tm and an initial guess for
+        the interfacial energy IE.
+        """
+        self.Ts = series.Ts
+        self.Js = series.Js
+        self.uJs = series.uJs
+        self.L, self.Tm, self.IE = L, Tm, IE
+    
+    def fit(self):
+        """ Perform the fit. """
+        Ts, Js, uJs = self.Ts, self.Js, self.uJs
+        # Functions for nucleation rate from CNT
+        J = lambda T,b,IE,L,Tm: (b*bulkFreeEnergy(T,L,Tm)**2/(8*np.pi) *
+                                1/np.sqrt(IE**3*kB*T) * 
+                                np.exp(-criticalFreeEnergy(T,IE,L,Tm)/kB/T))
+        # Fitfunction with L and Tm from literature
+        fitfunc = lambda T,b,IE: J(T,b,IE,self.L,self.Tm)
+        # Guess for IE is literature value, choose initial prefactor to fit some
+        # arbitary point. Do the fit.
+        b0 = Js[len(Js)//2]/fitfunc(Ts[len(Js)//2],1,self.IE)
+        popt, pcov = scipy.optimize.curve_fit(fitfunc, Ts, Js, maxfev=60000,
+                                              p0=[b0,self.IE])
+        self.b, self.IE = popt
+        self.ub, self.uIE = np.sqrt(pcov[0,0]), np.sqrt(pcov[1,1])
+#        self.b, self.IE = b0, self.IE
+#        self.ub, self.uIE = 1, 1
+        return self
+    
+    def J(self,T):
+        """
+        Nucleation rate as a temperature dependent function from this fit.
+        """
+        b, IE, L, Tm = self.b, self.IE, self.L, self.Tm
+        return (b*bulkFreeEnergy(T,L,Tm)**2/(8*np.pi) *
+                1/np.sqrt(IE**3*kB*T) * 
+                np.exp(-criticalFreeEnergy(T,IE,L,Tm)/kB/T))
+    
+    def info(self):
+        """ Print some information. """
+        print("Fitting J = b*dGv**2/(8*pi)/sqrt(IE**3*kB*T) * exp(-dGc/kB/T)")
+        print("        dGc = 16*pi*IE^3/(3*dGv)")
+        print("        dGv = L*(Tm-T)/Tm*2*T/(Tm+T)")
+        print("b: ",self.b,"+-",self.ub)
+        print("IE:",self.IE,"+-",self.uIE)
+        return self
+    
+    def plot(self, axes=None):
+        """ Plot the problem and fit. """
+        Ts, Js = self.Ts, self.Js
+        if axes is None:
+            fig, axes = plt.subplots(1,1)
+            plotrange = np.min(Ts), np.max(Ts)
+        else:
+            plotrange = axes.get_xlim()
+        axes.plot(Ts, Js, linestyle="none", marker=".", label="mean rates")
+        pxs = np.linspace(plotrange[0], plotrange[1], 200)
+        axes.plot(pxs, self.J(pxs), label="fit")
+        axes.set_xlabel("$T \\ [K]$")
+        axes.set_ylabel("$J \\ [1/m^3/s]$")
+        return axes
+
+class MobilityFit:
+    """ Fitting an Arrhenius behaviour to the mobility M = v / dGv """
+    
+    def __init__(self, series, dGv):
+        """
+        Build an object of a given series of measurements and a given model for
+        the free energy dGv as a function fo temperature.
+        """
+        self.Ts = series.Ts
+        self.vs = series.vs
+        self.uvs = series.uvs
+        self.Ms = self.vs / dGv(self.Ts)
+        self.uMs = self.uvs / dGv(self.Ts)
+    
+    def linfit(self):
+        """ Fit the linear problem. """
+        m, um, c, uc, _, _ = linreg(1/self.Ts, np.log(self.Ms), self.uMs)
+        self.Minf, self.uMinf = np.exp(c), np.exp(c)*uc
+        self.Ea, self.uEa = -m, um
+        return self
+    
+    def fit(self):
+        """ Perform a fit. """
+        self.linfit()
+        fitfunc = lambda x,a,b: a*np.exp(-b/x)
+        Ts, Ms, uMs = self.Ts, self.Ms, self.uMs
+        p0 = [self.Minf, self.Ea]
+        popt, pcov = scipy.optimize.curve_fit(fitfunc, Ts, Ms, p0=p0,
+                                              sigma=uMs, absolute_sigma=True)
+        self.Minf, self.Ea = popt[0], popt[1]
+        self.uMinf, self.uEa = np.sqrt(pcov[0,0]), np.sqrt(pcov[1,1])
+        return self
+    
+    def M(self,T):
+        """
+        Mobility as temperature dependent function according to this fit results.
+        """
+        return self.Minf * np.exp(-self.Ea/T)
+    
+    def info(self):
+        """ Print some information. """
+        print("Fitting M = Minf*exp(-Ea/T)")
+        print("Minf = {} +- {}".format(self.Minf, self.uMinf))
+        print("Ea   = {} +- {}".format(self.Ea, self.uEa))
+        Ea, uEa = [x/eV*kB for x in [self.Ea, self.uEa]]
+        print("   = (",Ea,"+-",uEa,")eV/kB")
+        print("chisq/ndf:", chisq_ndf(self.Ts, self.Ms, self.uMs,
+                                      (lambda x,a,b: a*np.exp(-b/x)),
+                                      [self.Minf, self.Ea]))
+        return self
+    
+    def plot_lin(self, axes=None):
+        """ Plot the linearised problem. """
+        if axes is None:
+            fig, axes = plt.subplots(1,1)
+        invTs = 1/self.Ts
+        m, c = -self.Ea, np.log(self.Minf)
+        axes.errorbar(invTs, np.log(self.Ms), self.uMs/self.Ms,
+                      linestyle="none", marker="_", label="data")
+        axes.plot(invTs, m*invTs+c, label="fit")
+        axes.set_xlabel("$1/(T \\ [K])$")
+        axes.set_ylabel("$\\log(M \\ \\left[\\frac{m^3}{J s}\\right])$")
+        return axes
+    
+    def plot(self, axes=None):
+        """ Show a plot. """
+        if axes is None:
+            fig, axes = plt.subplots(1,1)
+        axes.errorbar(self.Ts, self.Ms, self.uMs, linestyle="none", marker="_",
+                     label="data")
+        pxs = np.linspace(np.min(self.Ts)-1,np.max(self.Ts)+1,100)
+        axes.plot(pxs,self.M(pxs), label="fit")
+        axes.set_xlabel("$T \\ [K]$")
+        axes.set_ylabel("$M \\ \\left[\\frac{m^3}{J s}\\right]$")
+        return axes        
+
+__all__ = ["SeriesOfMeasurements", "ArrheniusFit", "SimpleCntFit",
+           "MobilityFit"]

py/measurementclasses.py

-"""
-This file provides classes to analyse experimental or simulated data on the
-crystallisation of a phase change material.
-"""
-
-import csv
-import numpy as np
-import matplotlib.pyplot as plt
-import scipy.optimize
-from scripts import linreg, weightedmean, chisq_ndf
-
-
-class Measurement:
-    """
-    A class for a measurement for one temperature as I got from Julian or as
-    returned from one simulation run.
-    Properties:
-      times - times of measurements [s]
-      nuclei - number of seen nuclei
-      radii - 2D array of radii [m]
-      afract - amorphous fraction
-      cfract - crystallized fraction
-    Functions:
-      plot_fc
-      plot_N
-      plot_r
-      JMAK_linfit
-      JMAK_fit
-      rfit
-      all_rfits
-      nucleationratefit
-    """
-
-    def plot_fc(self):
-        """ Plot the crystallised fraction over time. """
-        plt.plot(self.times, self.cfract, linestyle="none", marker="x")
-        plt.xlabel("$t \\ [s]$")
-        plt.ylabel("$f_c$")
-        plt.tight_layout()
-        
-    def plot_N(self):
-        """ Plot the counted nuclei over time. """
-        plt.plot(self.times, self.nuclei, linestyle="none", marker="x")
-        plt.xlabel("$t \\ [s]$")
-        plt.ylabel("$N$")
-        plt.tight_layout()
-    
-    def plot_r(self, i):
-        """ Plot the size of the i-th nucleus over time. """
-        plt.plot(self.times, self.radii[i]*1e6, linestyle="none", marker="x")
-        plt.xlabel("$t \\ [s]$")
-        plt.ylabel("$r \\ [\\mu m]$")
-        plt.tight_layout()
-    
-    def JMAK_linfit(self, show=True):
-        """
-        Linearise the fit to a JMAK model and use linear regression.
-        Uncertainties are given so that chisquared per degree of freedom is 1.
-        Returns (n, k) for f_a = exp(-k*t^n)
-        """
-        
-        # Extract points where log() is well defined
-        cond = np.logical_and( np.logical_and( self.afract > 0,
-                                               self.afract < 1),
-                               self.times > 0)
-        xs = np.log( np.compress(cond, self.times))
-        fs = np.compress(cond, self.afract)
-        ys = np.log( - np.log( fs ))
-        # Assume equal uncertainty on area values, transform to y
-        sigmas = np.abs(np.ones(xs.shape)/np.log(fs)/fs)
-        n, un, c, uc, corr, chisq = linreg(xs, ys, sigmas)
-        chisqndf = chisq / (len(xs)-2)
-        # correct uncertainties, calculate k
-        un, uc = un*np.sqrt(chisqndf), uc*np.sqrt(chisqndf)
-        k, uk = np.exp(c), np.exp(c)*uc
-        if show:
-            print("f_a = exp(-k*t^n), linearise:")
-            print("ln(-ln(f_a)) = n * ln(t) + ln(k)")
-            print("n:", n, "+-", un)
-            print("k:", k, "+-", uk)
-            plt.plot(xs, ys, linestyle="none", marker="x", label="data")
-            plt.plot(xs, n*xs+c, label="fit")
-            plt.xlabel('$\ln(t \\ [s])$')
-            plt.ylabel('$\ln( - \ln( 1 - f_c ))$')
-            plt.legend()
-            plt.tight_layout()
-        return n, k
-        
-    
-    def JMAK_fit(self, dimensions=None, show=True):
-        """
-        Fit to a JMAK model.
-        Either with a fixed number of dimensions or using this as another free
-        parameter.
-        Uncertainties are given so that chisquared per degree of freedom is 1.
-        Returns ([n,k],pcov) where pcov is the covariance matix.
-        """
-        
-        # Initial values from linear problem, function to fit onto
-        n0, k0 = self.JMAK_linfit(show=False)
-        fitfunc = lambda t,n,k: np.exp(-k*t**n)
-        if dimensions is None:
-            popt, pcov = scipy.optimize.curve_fit(fitfunc, self.times,
-                                self.afract, p0=[n0,k0])
-        else:
-            # If dimension is given, calculate exponent and put into fitfunction
-            n = dimensions + 1
-            fitfunc2 = lambda t,k: fitfunc(t,n,k)
-            popt, pcov = scipy.optimize.curve_fit(fitfunc2, self.times, 
-                                                  self.afract, p0=[k0])
-            # Fill popt, pcov to the same size as from above
-            popt = np.insert(popt, 0, n)
-            pcov = np.array([[0, 0], [0, pcov[0,0]]])
-        # Get chisquared and correct uncertainties
-        chisq = chisq_ndf(self.times, self.afract, np.ones(self.afract.shape),
-                          fitfunc, popt)
-        pcov = pcov*chisq
-        if show:
-            print("f_a = exp(-k*t^n)")
-            print("n:", popt[0], "+-", np.sqrt(pcov[0,0]))
-            print("k:", popt[1], "+-", np.sqrt(pcov[1,1]))
-            print("covariance: ", pcov[0,1])
-            plt.plot(self.times, self.afract, linestyle="none", marker="x",
-                     label="data")
-            plottimes = np.linspace(np.min(self.times), np.max(self.times), 100)
-            plt.plot(plottimes, fitfunc(plottimes, *popt), label="fit")
-            plt.xlabel("$t \\ [s]$")
-            plt.ylabel("$1-f_c$")
-            plt.legend()
-            plt.tight_layout()
-        return popt, pcov
-    
-    def rfit(self, i, show=True):
-        """
-        Fit a linear growth to the i-th grain.
-        Uncertainties are given so that chisquared per degree of freedom is 1.
-        Returns (v, sigma_v).
-        """
-        
-        cond = np.logical_not( np.isnan(self.radii[i]) )
-        ts = np.compress(cond, self.times)
-        rs = np.compress(cond, self.radii[i])
-        v, uv, c, uc, corr, chiq = linreg(ts, rs, np.ones(rs.shape))
-        # Scale to make chisquared/ndf equals 1
-        chiq = chiq/(len(ts)-2)
-        uv, uc = [x * np.sqrt(chiq) for x in [uv, uc]]
-        if show:
-            print("v:", v, "+-", uv)
-            plt.plot(ts, rs*1e6, linestyle="none", marker="x", label="data")
-            plt.plot(ts, (v*ts+c)*1e6, label="fit", color="black", alpha=0.7)
-            plt.xlabel("$t \\ [s]$")
-            plt.ylabel("$r \\ [\\mu m]$")
-            if show != 2:
-                plt.legend()
-            plt.tight_layout()
-            plt.show()
-        return v, uv
-    
-    def all_rfits(self, show=True):
-        """ Fit all grains and calculate mean growht velocity """
-        
-        # Initialise results array and do first fit
-        vs = np.zeros(len(self.radii))
-        vs[0] = self.rfit(0, show=show)[0]
-        # For all other fits, do not show a legend anymore
-        for i in range(1,len(self.radii)):
-            vs[i] = self.rfit(i, show=show*2)[0]
-        v, uv = np.mean(vs), np.std(vs)
-        return v, uv
-    
-    def nucleationratefit(self, splitat=None, show=True):
-        """
-        Fit a nucleation rate to the number of grains. Either to the whole
-        domain, or give points to split at. Will return a list of 2-tuples
-        each containing a rate and an uncertainty.
-        """
-        
-        # Calculate the new nuclei per time step
-        dN = np.diff(self.nuclei)
-        dt = np.diff(self.times)
-        meantimes = (self.times[1:]+self.times[:-1])/2
-        # Scale with fraction of material where nucleation can occour
-        fa = (self.afract[1:] + self.afract[:-1]) / 2
-        Js = dN/dt * 1/fa
-        # Integrate again
-        N2 = np.cumsum(Js * dt)
-        N2 = np.insert(N2,0,0)
-        # Assume poisson uncertainty
-        uN2 = np.maximum(np.sqrt(N2), 1)
-        
-        # If no points to split at are given, create an empty list
-        if splitat is None:
-            splitat = []
-        # Calculate domain boundarys: Insert first time at front, cut out every
-        # time thats out of available times and insert last time at end
-        splitat = np.insert(splitat, 0, self.times[0])
-        splitat = np.compress(splitat<self.times[-1], splitat)
-        splitat = np.append(splitat,self.times[-1])
-        
-        # Show the data
-        if show:
-            plt.errorbar(self.times, N2, uN2, linestyle=":", label="data")
-            plt.xlabel("$t \\ [s]$")
-            plt.ylabel("$N_{corrected}$")
-            plt.tight_layout()
-        
-        # For each of the domains
-        res = []
-        y0, y1 = 0, 0
-        for tmin, tmax in zip(splitat[0:-1],splitat[1:]):
-            # extract data in this domain
-            cond = np.logical_and(meantimes >= tmin,
-                                  meantimes <= tmax)
-            Js_i = np.compress(cond, Js)
-            dN_i = np.compress(cond, dN)
-            dt_i = np.compress(cond, dt)
-            fa_i = np.compress(cond,fa)
-            
-            # Assume poisson uncertainty with the mean rate on all points
-            uJ_i = np.sqrt(np.maximum(np.mean(dN_i),1))/dt_i * 1/fa_i
-            
-            # Mean rate
-            J, uJ = weightedmean(Js_i, uJ_i)
-            
-            if show:
-                print("from",tmin,"to",tmax,":")
-                print("   rate:",J,"+-",uJ)
-                y0, y1 = y1, y1+J*(tmax-tmin)
-                plt.plot([tmin,tmax], [y0, y1], label="integrated mean rate")
-                
-            res.append((J,uJ))
-            
-        return res
-
-
-# Calculate a length from the values in pixels, according to the mail I got.
-# First method is the correction I got in the second mail
-correctedpx = lambda px: 300/72*px
-px2len = lambda p: (0.3876 * correctedpx(p) - 0.3802) * 1e-6
-
-class Experiment(Measurement):
-    """
-    Class for experimental data as I got from Julian.
-    Additional properties:
-      columntitles - column titles from the csv
-      aarea - amorphous area [m^2]
-      A - overall observed area [m^2]
-      h - thickness of the material [m]
-    Additional methods:
-      compare_results
-    """
-    
-    def __init__(self, filename, height):
-        """ Read all data from the file and build variables of interest. """
-         
-        # Open file and build csv reader
-        file = open(filename, newline="")
-        reader = csv.reader(file, delimiter=",")
-        # First line has titles to the columns
-        self.columntitles = list(map(lambda x: x.strip(), next(reader)))
-        # Rest is data structured in columns
-        columns = [ [] for _ in self.columntitles]
-        for row in reader:
-            for column, value in zip(columns, row):
-                # Replace empty cells with nan to parse them to float
-                if value == "":
-                    value = "nan"
-                column.append(value.replace(",","."))
-        # First column are times in minutes, convert to seconds
-        self.times = np.array(columns[0], dtype=np.float) * 60
-        # Second column is amporphous area. The square root is a length and can
-        # be converted according to the mail, then square the result again
-        self.aarea = (px2len(np.sqrt(np.array(columns[1], dtype=np.float))))**2
-        # Third column has number of grains, here no value means probably zero
-        self.nuclei = np.array(columns[2], dtype=np.float)
-        self.nuclei[np.isnan(self.nuclei)] = 0
-        # Rest of the file are radii of grains
-        self.radii = px2len( np.array(columns[3:], dtype=np.float) )
-        # Overall observed area
-        self.A = max(self.aarea)
-        self.h = height
-        # Fractions of amorphous and crystaline area
-        self.afract = self.aarea/self.A
-        self.cfract = 1-self.afract
-    
-    def compare_results(self, dimensions):
-        """
-        Compare the results from the JMAK fit with the ones obtained from
-        fitting on grain count and sizes.
-        """
-        v, uv = self.all_rfits(show=False)
-        [(dNdt, udNdt)] = self.nucleationratefit(show=False)
-        J  =  dNdt / (self.A*self.h)
-        uJ = udNdt / (self.A*self.h)
-        [n,k1], pcov = self.JMAK_fit(dimensions,show=False)
-        
-        if dimensions == 1:
-            func = lambda t: np.exp(-J*self.h**2*v*t**2)
-            k2 = J*self.h**2*v
-            uk2 = k2 * np.sqrt(uJ**2/J**2 + uv**2/v**2)
-        elif dimensions == 2:
-            func = lambda t: np.exp(-np.pi/3*J*self.h*v**2*t**3)
-            k2 = np.pi/3*J*self.h*v**2
-            uk2 = k2 * np.sqrt(uJ**2/J**2 + 4*uv**2/v**2)
-        elif dimensions == 3:
-            func = lambda t: np.exp(-np.pi/3*J*v**3*t**4)
-            k2 = np.pi/3*J*v**3
-            uk2 = k2 * np.sqrt(uJ**2/J**2 + 9*uv**2/v**2)
-        
-        print("JMAK fit:   k = {} +- {}".format(k1,np.sqrt(pcov[0,0])))
-        print("v & J fits: k = {} +- {}".format(k2,uk2))
-        
-        plottimes = np.linspace(np.min(self.times), np.max(self.times), 100)
-        plt.plot(self.times, self.afract, linestyle="none", marker="x", 
-                 label="data")
-        plt.plot(plottimes, np.exp(-k1*plottimes**n),
-                 label="fit on $f_c$")
-        plt.plot(plottimes, func(plottimes),
-                 label="fits on $\\frac{dN}{dt}$ and $v$")
-        plt.legend()
-
-
-class Simulation(Measurement):
-    """
-    Class for simulation results as received from my code.
-    Additional properties:
-        columntitles - column titles from the csv
-    """
-    
-    def __init__(self, filename):
-        """ Read all data from the file and build variables of interest. """
-        
-        # Open file and build csv reader
-        file = open(filename, newline="")
-        reader = csv.reader(file, delimiter=",")
-        # First line has titles to the columns
-        self.columntitles = list(map(lambda x: x.strip(), next(reader)))
-        # Rest is data structured in columns
-        lines = []
-        for row in reader:
-            lines.append(row)
-        data = np.array(lines, dtype=np.float)
-        # First column are times
-        self.times = data[:,0]
-        # Second column has number of grains
-        self.nuclei = data[:,1]
-        # Third column is crystalline fraction
-        self.cfract = data[:,2]
-        self.afract = 1 - self.cfract
-        # Rest of the file are sizes of grains
-        self.radii = np.sqrt(data[:,3:]/np.pi).transpose()