Multifit Model Example

// <copyright file="MyMultifitModel.cs" company="rhd instruments GmbH and Co. KG">
// Copyright (c) rhd instruments GmbH and Co. KG. All rights reserved.
// Licensed under the MIT No Attribution (MIT-0) license. See section 'License' in the 'SDK Examples / Tutorials' topic for full license information.
// </copyright>

namespace RelaxIS_SDK_Examples.Plugins
{
    using System;
    using System.Collections.Generic;
    using System.Linq;
    using System.Text;
    using System.Windows.Forms;
    using RelaxIS_SDK.Common;
    using RelaxIS_SDK.libMath;
    using RelaxIS_SDK.Plugins;

    /***
     * MultifitModels are a special type of fitting model used by the Multifit-function of RelaxIS.
     * It is used to fit not just one but multiple spectra simultaneously to one combined model.
     * The model can use spectras metadata as additional values for the calculation.
     * The model can have parameters that apply to all spectra equally and parameters individual for each spectrum.
     *
     * A good example (also implemented here) is to determine one of the parameters via an
     * Arrhenius-function. That way, the model has one set of Arrhenius parameters that is used to
     * calculate this value for each of the spectra depending on the Temperature metadata field.
     * This implicitely forces e.g. a bulk resistance parameter to adhere strictly to Arrhenius
     * behavior, instead of moving the Arrhenius fit into a separate step after the individual spectrum fits.
     */

    /// <summary>
    /// Defines an example <see cref="RelaxISPlugin_MultifitModel"/> class.
    /// </summary>
    public class MyMultifitModel
        : RelaxISPlugin_MultifitModel
    {
        /// <summary>
        /// Defines the metadata required on spectra for this plugin to be available.
        /// </summary>
        private static readonly string[] Metadata = new[] { MetadataNames.TEMPERATURE };

        /***
         * First, implement the default plugin properties Name and Description that describe the plugin.
         * This is mainly used for display purposes.
         */

        /// <inheritdoc/>
        public override string Name
        {
            get { return "MyMultifitModel"; }
        }

        /// <inheritdoc/>
        public override string Description
        {
            get { return "Implements a MultiFit model that calculates the series resistance of an RP model based on an Arrhenius function"; }
        }

        /***
         * Implement properties that tell RelaxIS which metadata is required and which transfer function should be used.
         */

        /// <inheritdoc/>
        public override List<string> RequiredMetadata
        {
            get
            {
                return new List<string>(Metadata);
            }
        }

        /// <inheritdoc/>
        public override string TransferFunction
        {
            get
            {
                // Refer to the RelaxIS settings -> Plugins window for the internal names of the transfer functions.
                return "Impedance";
            }
        }

        /***
         * This function gets a set of individual parameters for a specific spectrum. In most cases all spectra have the same set of
         * individual parameters, but that is not a requirement. Given that the resistance of the RP model in this example is
         * a shared parameter (see below), the individual parameters are just the CPE parameters.
         * Note, that this function also determines minimal and maxima, as well as initial values for the parameters.
         * This is helpful to make initialization of the parameters easier in the user interface, since these need to be initialized
         * manually for each of the spectra! Having a good automatic initialization algorithm is therefore very helpful.
         */

        /// <inheritdoc/>
        public override List<Fitparameter> GetIndividualParameters(ImpedanceSpectrum Spectrum, List<ImpedanceSpectrum> AllSpectra)
        {
            // Calculate the complex capacitances from the data
            var caps = Spectrum.Data.Select(p => 1.0 / (new Complex(0, 2.0 * Math.PI * p.Impedance.Frequency) * new Complex(p.Impedance.Real, p.Impedance.Imaginary)));

            // Take the minimum and maximum real capacitances and scale them by a factor for the new limits
            var cAvg = caps.Average(c => c.Real);
            var cMin = caps.Min(c => c.Real) / 10000;
            var cMax = caps.Max(c => c.Real) * 10000;

            // Individual parameters are a CPE Q and alpha for each spectrum
            var res = new List<Fitparameter>()
            {
                new Fitparameter("Q", false, cAvg, 0, cMin, cMax),
                new Fitparameter("alpha", false, 0.95, 0, 0.3, 1),
            };

            return res;
        }

        /***
         * This function defines the shared parameters. These will only appear once in the fit model calculation and hence apply to all spectra.
         * In this example these are the Arrhenius parameters that will be used to calculate a series resistance value.
         * This function also initializes the values via a linear fit to the logarithmic average real part of the spectra vs. the inverse temperature.
         */

        /// <inheritdoc/>
        public override List<Fitparameter> GetSharedParameters(List<ImpedanceSpectrum> Spectra)
        {
            // Select inverse temperature
            var temperatures = Spectra.Select(spectrum => 1.0 / spectrum.Metadata.Where(m => m.Name == "Temperature").First().Value).ToList();

            // Select log of average real parts of each spectrum
            var resistances = Spectra.Select(spectrum => spectrum.Data.Select(d => Math.Log(d.Impedance.Real)).Average()).ToList();

            // R = A*Exp(-Ea/(RT)) => ln(R) = ln(A) - Ea/(RT)
            // Linear fit: Ea = -Slope * R, A = Exp(Intercept)
            var reg = GetLinearRegression(temperatures, resistances);
            var ea = 1e-4;
            var ar = 100.0;
            if (!double.IsNaN(reg.Item1) && !double.IsNaN(reg.Item2))
            {
                ea = -reg.Item1 * 8.314;
                ar = Math.Exp(reg.Item2);
            }

            // Arrhenius pre-exponential factor and activation energy is shared across all spectra.
            var res = new List<Fitparameter>()
            {
                new Fitparameter("A_R", false, ar, 0, 1e-15, 1e15),
                new Fitparameter("Ea", false, ea, 0, 1e-15, 1e15),
            };

            return res;
        }

        /***
         * Now the actual fit function needs to be implemented.
         * The X[] parameters contain the frequency in index 0, the spectrum index at index 1, and the metadata as defined in the
         * RequiredMetadata in the subsequent indizes. Here we therefore find the temperature at index 2.
         * The Parameters[] argument contains the full list of parameters: First the shared parameters in the order defined in
         * GetSharedParameters(), and then the individual parameters as returned by GetIndividualParameters() for each spectrum.
         * The order of spectra is the same over all function. If GetIndividualParameters() returned the same number of parameters for
         * each spectrum, the GetIndividualParameterIndex() helper function can be used to determine the parameter indizes of a spectrums
         * parameters in the overall lists.
         * If spectra have different numbers of individual parameters, you can alternatively use the ParameterStartIndizes property that is
         * populated by RelaxIS during the setup phase.
         */

        /// <inheritdoc/>
        public override Complex FitFunction(double[] X, double[] Parameters)
        {
            var f = X[0];
            var iw = new Complex(0, 2 * Math.PI * f);
            var spectrumIdx = (int)X[1];
            var t = X[2];
            var ar = Parameters[0];
            var ea = Parameters[1];

            // Determine the parameters of the individual parameters

            // Given that all spectra have 2 individual parameters we can use the helper function:
            var qIdx = this.GetIndividualParameterIndex(2, 2, 0, spectrumIdx);
            var aIdx = this.GetIndividualParameterIndex(2, 2, 1, spectrumIdx);

            // Alternativly the ParameterStartIndizes dictionary can be used:
            var startIdx = this.ParameterStartIndizes[spectrumIdx];
            qIdx = startIdx + 0;
            aIdx = startIdx + 1;

            var q = Parameters[qIdx];
            var a = Parameters[aIdx];

            // Calculate R using the Arrhenius function
            var r = ar * Math.Exp(-ea / (8.314 * t));

            // Calculate the CPE value normally.
            var cpe = 1.0 / (Complex.Pow(iw, a) * q);
            return r + cpe;
        }

        /***
         * The PerformSetup function is called by RelaxIS before the GetIndividualParameters and GetSharedParameters functions are called and
         * can be used to e.g. perform initialization functions or show a user dialog with the given parent window.
         * Setting Cancel to true will stop the fit execution.
         */

        /// <inheritdoc/>
        public override void PerformSetup(IEnumerable<ImpedanceSpectrum> Spectra, IWin32Window parent, ref bool Cancel)
        {
            var count = Spectra.Count();
            var sb = new StringBuilder();
            sb.AppendLine("You could do some setup steps here that run before the parameter initialization step.");
            sb.AppendFormat("You have selected {0} spectra for fitting", count);
            sb.AppendFormat("This will result in a total of 2 shared and {0} individual fit parameters", count * 2);
            MessageBox.Show(parent, sb.ToString(), "Information", MessageBoxButtons.OK, MessageBoxIcon.Information);
        }

        /***
         * The AfterInitialization function is called by RelaxIS after the GetIndividualParameters and GetSharedParameters functions are called and
         * can be used to e.g. perform post-initialization functions such as saving settings or to show a user dialog with the given parent window.
         * Setting Cancel to true will stop the fit execution.
         */

        /// <inheritdoc/>
        public override void AfterInitialization(IEnumerable<ImpedanceSpectrum> Spectra, IWin32Window parent, List<Fitparameter> SharedParameters, List<Tuple<ImpedanceSpectrum, List<Fitparameter>>> IndividualParameters, ref bool Cancel)
        {
            var sb = new StringBuilder();
            sb.AppendLine("You could do some parameter validation here and possibly cancel the fitting process.");
            MessageBox.Show(parent, sb.ToString(), "Information", MessageBoxButtons.OK, MessageBoxIcon.Information);
        }

        /***
         * This function is called after the fit has completed, allowing the calculation of additional values from the fit results.
         */

        /// <inheritdoc/>
        public override Dictionary<string, double> GetAdditionalEvaluations(ImpedanceSpectrum Spectrum, List<ImpedanceSpectrum> AllSpectra, List<Fitparameter> SharedParameters, List<Fitparameter> IndividualParameters)
        {
            // With SI units the activation energy would be in J/mol
            var ea = SharedParameters[1].Value;
            var ev = ea * 6.2415093433e+18 / 6.022e23;
            var res = new Dictionary<string, double>
            {
                { "Ea [kJ/mol]", ea / 1000.0 },
                { "Ea [eV]", ev },
            };
            return res;
        }

        /***
         * Helper function to perform a linear least squares fit.
         */

        /// <summary>
        /// Calculate the unweighted linear regression over the data.
        /// </summary>
        /// <param name="x">The x (independent) values.</param>
        /// <param name="y">The y (dependent values).</param>
        /// <returns>The regression result, Slope | Intercept.</returns>
        private static Tuple<double, double> GetLinearRegression(IList<double> x, IList<double> y)
        {
            if (x.Count != y.Count)
            {
                throw new ArgumentException("The data lists have different length");
            }

            var xAvg = x.Average();
            var yAvg = y.Average();
            double numerator, denominator;
            numerator = 0;
            denominator = 0;
            var dCount = x.Count;
            for (int i = 0; i < dCount; i++)
            {
                var dx = x[i] - xAvg;
                numerator += dx * (y[i] - yAvg);
                denominator += dx * dx;
            }

            double slope;
            bool nan = false;
            if (numerator != 0 & denominator == 0)
            {
                slope = 0;
            }
            else if (numerator == 0 & denominator == 0)
            {
                nan = true;
                slope = double.PositiveInfinity;
            }
            else
            {
                slope = numerator / denominator;
            }

            if (!nan)
            {
                double b = yAvg - (slope * xAvg);
                return Tuple.Create(slope, b);
            }
            else
            {
                return Tuple.Create(double.NaN, double.NaN);
            }
        }
    }
}