Parameters for MEM Inversions

The MEM inversion is an iterative process that alternates between optimizations of the function Q for given values of the Lagrange multipliers and updates to the multiplier values for the next optimization. Iteration continues up to a maximum number of steps (MXFUNC), until C reaches a target value (CHI2AM), or until C stops changing more than a specified amount (CHI2S).

The following paragraphs describe, in order, the line-by-line information needed in an inversion .def file. Remember that running MemExp in simple mode produces a '.mem' file that contains what are likely to be acceptable values for these inversion parameters.

The names of the variables whose values are to be read from the given line are written in bold capital letters. If multiple variables are to be input, do not separate them by commas. A semicolon can be used following the data to precede comments. The information in a few lines specifies a variable number of optional lines that may immediately follow.

The first line specifies IOFORM, the format of the data, which is either RW or Fr (rdwrt or free format). Fr format is almost certainly the format you will use, requiring a file without any header information with two or three numbers per line, specifying t, y(t), and optionally $\sigma(t)$.

The second line specifies the integer ITAKE. If 1, fit to every point in the data files. If 2, every other point, etc.

XSCALE, YSCALE scale the x/y axis of the data (eg., 0.001 would convert from milliseconds to seconds).

MLOG is set to T, true, only for data in transmittance space. For most data, MLOG is F.

IGNORE should be set to 0, 1, 2, or 3. If IGNORE is set to 0, the standard errors $\sigma(t)$ read from the data file are used (not ignored). If IGNORE is set to 1, any uncertainties $\sigma(t)$ (standard errors) in the data file are ignored and replaced as specified by the value of INSIG that is read from the next line. If INSIG $= -2.0$ (recommended, new to version 5.0), then time-dependent standard errors are estimated from the differences between the raw data and the data after smoothing. If INSIG $= -1.0$, then uniform standard errors are estimated from the variance of the last three data points. Otherwise, if INSIG is negative, it specifies the absolute error common to all data, i.e., $\sigma(t) = -INSIG$. If INSIG is positive, it represents the per-cent error: $\sigma(t) = 0.01 \times INSIG \times y(t)$. If INSIG is zero, the data are assumed free of noise. In this case, the data are actually assigned uniform errors (equal to 0.001 times the second data value), and CHI2AM and DCHI2S are reduced so that the calculation proceeds to very small values of CHI2. If IGNORE is set to 2, Poisson statistics are used to treat the noise, and two values are then read from the next line: The Hessian of the Poisson deviance is updated only every HFREQ calls to FUNC until the deviance drops to HDEV, when the Hessian is updated every step. If IGNORE is set to 3, time-independent errors are used.

The integer D0OPT specifies how D0 is to be estimated from the data if the input value of D0 is zero. If D0 = 0 and D0OPT = 1, then D0 is set to the difference between the first and last data values. If D0 = 0 and D0OPT = 2, then D0 is set to the first data value, i.e., for data that would decay to zero after a sufficiently long time. WARNING: These approximations will underestimate D0 for data that have appreciable slope at short times when plotted versus log time. For such data, D0 should be set on the command line to a larger value and LAGMUL should be set to Chi2. Severely underestimating D0 can induce spurious maxima in f [10].

LAGMUL is either Chi2 or Chi2,Sumf specifying that only C is to be constrained or that the image normalization is also to be constrained (equation 3).

CONST specifies the value of the initial uniform value of $f(log \tau)$. If set to zero, CONST is reset automatically.
Recommendation: Set CONST to 0.001.

CHI2AM is the desired (aimed at) value of C. Iteration terminates if this value is reached.

If convergence is not reached, iteration will be terminated after MXFUNC evaluations of Q.

EPSL and EPSU are the logarithms of $\epsilon$ used to limit updates to the Lagrange multipliers $\lambda$ and $\alpha$ (Cornwell and Evans, 1985). They are lower and upper values of the tolerance allowed on the gradient of Q, the function being optimized. Large values speed up the calculation but do not optimize Q as well. EPSL is used at low values of C, and EPSU is used at larger values of C. Thus, by setting EPSL $<$ EPSU, the program proceeds quickly to low values of C before demanding a more precise optimization of Q. If $C \ge$ CHI2U, then $\epsilon = 10^{EPSU}$. Otherwise, $log(\epsilon)$ is ramped from EPSU down to EPSL as C drops to CHI2L. If $C \le$ CHI2L, then $\epsilon = 10^{EPSL}$. Clearly, if EPSL and EPSU are equal, the values of CHI2L and CHI2U are irrelevant.

PPLT is the number of equally spaced values per unit log $\tau$ used to discretize the $g(log \tau)$ distribution that ranges from LTMIN to LTMAX.
Recommendation: Output image files can be kept small without loss of precision by choosing PPLT so that the spacing in $log \tau$ is an integer multiple of 0.01 or 0.001. The $log \tau$ values are then written in F6.2 or F7.3 format, respectively. Otherwise, they are written as unformatted output, resulting in larger files. Examples: PPLT $= 40$ or PPLT $= 25$. The same holds for PPLT2, LTMIN2, and LTMAX2 (described in the following paragraph).

NDIST, NBASLN, VBASLN, TBASLN: NDIST is an integer, either 1 or 2, specifying the number of distributions to be determined. If NDIST = 1, a single distribution $g(log \tau)$ is used, and the sign of D0 determines whether the data to be fit are rising or decaying (D0 $>$ 0). If NDIST = 2, an additional line is required immediately following the current line. This additional line must specify PPLT2, LTMIN2, and LTMAX2 to be used for the $h(log \tau)$ distribution that describes rising signals. (NOTE: MemExp ensures the spacing in $log \tau$ used to discretize the h distribution is the same as the spacing used for the g distribution.) NBASLN is an integer ranging from 0 (to omit baseline correction) to 4 (to include a cubic baseline). A constant, linear, and parabolic baseline are computed when NBASLN is 1, 2, and 3, respectively. Since the entropy involves a logarithm, only positive variables are possible, and so two variables must be used per polynomial degree of freedeom. The total number of 'image pixels' is therefore LTPTS $+$ 2*NBASLN if NDIST is 1 and LTPTS $+$ LTPTS2 $+$ 2*NBASLN if NDIST is 2. New to version 5.0, NBASLN can be set to -1, followed immediately by the value of OFFSET, which specifies a fixed, constant baseline value. OFFSET is only specified if NBASLN = -1. If VBASLN $\le$ 0, then the baseline parameters are fixed at the MEM-determined values during NLS/ML fits. If TBASLN $<$ 0, then the baseline is determined by the MEM in a simultaneous fit of distributed kinetics plus baseline drift. If TBASLN $>$ 0, the baseline is prefit at all times greater than TBASLN using linear least squares.
Recommendations: VBASLN $> 0$ and TBASLN $< 0$.

If RESPON $=$ 0, then an instrument response function (IRF) will not be deconvolved. If RESPON $<$ 0, the IRF (with values equally spaced in time) is read from the file named IRF-file-name specified on the following line, and EPSIRF, IRSMOTH, and LTSCAT are read from a second additional line. Then TSHIFT is read from a third additional line. If IRSMOTH $>$ 0, the IRF is smoothed by convolution with a bell curve ranging over IRSMOTH points, and the result is plotted in a PostScript file named data-file-name.irf.ps. All values of the (smoothed but not yet renormalized) IRF less than EPSIRF are neglected in the deconvolution from the data. LTSCAT is an integer; if positive, then the light-scattering correction is used (equation 2) by varying the value of the parameter $\xi$. However, if $\xi$ becomes less than 1.0E-30 in the MEM fits, then this correction is not used in the ML fits which can become unstable if $\xi$ adopts unphysical, negative values.

Set TSHIFT to zero to forego use of a zero-time shift ($\delta$ in equation 2) between the data and IRF. Set TSHIFT $> 10$ (recommended when deconvolving an IRF) to specify that the program should use a value of $\delta$ obtained from the recommended discrete fit.

CONSTB and SBIGF are relevant only if NBASLN $>$ 1, but values must always be specified. Together, they determine the value of BIGF to be used for the baseline coefficients other than the constant term. Setting CONSTB to a very small number (e.g., 1.0E-20) is a way to bias the baseline toward flatness; values of CONSTB as large as 0.001 can result in unwarranted time dependence in the baseline. SBIGF is a multiplicative factor less than or equal to one that progressively suppresses the BIGF value for the jth term of the baseline: BIGF(j) = CONSTB*SBIGF**j (j = 1 for linear, 2 for quadratic, 3 for cubic).

Up to MAXEXP components are to be used in the least-squares fits. The initial values of the discrete $log(\tau)$ are confined to the range EXPMIN and EXPMAX. If NDIST = 2, the initial values of the discrete $log(\tau)$ associated with negative amplitudes (assuming D0 $>$ 0) are confined to the range EXPMN2 and EXPMAX. Note, EXPMN2 is required input even if NDIST = 1.

Data to be included in the fit, those that contribute to C, are within times t1 and t2.

PLOTS is the maximum number of rows (three 2-D axes are plotted per row) in the output PostScript files. FSCALE can be used to restrict the maximum scale of (zoom in on) the lifetime-distribution plot in the right-hand column of the PostScript output. If FSCALE exceeds the maximum value of the distribution to be plotted, then the plot's y axis is automatically scaled to the distribution.

When greater than 1, IPLOT and IPLOT2 reduce the size of the PostScript files created without affecting the calculation. When plotting the data, fit, and baseline (when used) in the second column of the output PostScript file(s), every IPLOT$^{th}$ point is plotted. When plotting the residuals and the autocorrelation function of the residuals, the first three points are plotted and then one out of every IPLOT2 points is plotted. The point plotted for each group of IPLOT2 points is the one with maximum absolute value, ensuring that the worst-fit points are represented. Setting both IPLOT and IPLOT2 to 1 plots all points, producing the most accurate (and largest) graphics output.
Recommendation: Try IPLOT = 10 and IPLOT2 = 1, and adjust the values based on the frequency at which your data are sampled and your need to conserve disk space.

Values of $\lambda$, C, ${\cal S}$, Q, etc. are written to the output file every NPRINT optimization steps, and adjustments to $\lambda$ and $\alpha$ are at least as often as every NPRINT optimization steps.

The MEM image f is written to a file every NWRITE optimization steps, once C has reached CHIWRT. New to version 4.0: If CHIWRT $<$ 0, MemExp attempts to determine the values of NWRITE ( = NPRINT) and CHIWRT that will output results in PostScript format at appropriate intervals to allow the user to assess the results by visual inspection. The objective is to start writing images when the data are still under-fit ($\lambda$ too small) and to continue until the data are over-fit ($\lambda$ too large).
Recommendation: First, try CHIWRT = -1 (see above). Otherise, set CHIWRT = 1.2 or less.

BIGF is name of the array storing the 'prior' distribution that defines maximum entropy ${\cal S}$. When IBIGF = 1 (see below), BIGF is created by convolving F with a gaussian with a full width at half maximum (FWHM) of BFFWHM pixels. If BFFWHM $<$ 0 it is interpreted as the percentage of one log unit in $\tau$. This conveniently scales the blurring to the resolution of the image. Thus, setting BFFWHM to -20.0 will blur with a FWHM of one fifth of a log unit.

ASPIKE, SPIKE, SPIKE2, SPIKE3, DAREA, and PKFWHM are the default parameters that control the differential blurring of F into BIGF when IBIGF = 1 (below). For more convenient manipulation of BIGF, these default parameters can be overruled by setting NSPIKE $>$ 0 (below). But let's assume NSPIKE = 0. Then, if a MEM peak has area greater than ASPIKE, it will be accounted for in one of three ways, depending on the isolation of the peak(s) to be treated differently. Be sure to set SPIKE $>$ SPIKE2 $>$ SPIKE3. 1) Well separated MEM peaks that have a maximum-to-minimum ratio exceeding SPIKE on both sides of the peak are simply subtracted before the remainder of F is blurred uniformly. 2) Partially separated peaks have a maximum-to-minimum ratio exceeding SPIKE on one side and exceeding SPIKE2 on the other side. Moreover, the peak it overlaps is either small in area (less than ASPIKE) or has a maximum-to-minimum ratio exceeding SPIKE3 on its near side. These peaks are "folded over" from the well resolved side before subtraction. 3) Peaks that appear atop a broad feature are identified as having a maximum-to-minimum ratio exceeding SPIKE2 on both sides of the peak, but its neighboring peaks are not sharp (maximum-to-minimum ratios less than SPIKE3). The peak overlapping the broad feature is extracted prior to blurring by assuming a linear baseline that is determined automatically [3]. To compensate for the suppression of shoulders by the MEM, the linear baseline can be shifted up to reduce the area to be sharpened by approximately the multiplicative factor DAREA. Each peak subtracted in one of the above three ways is incorporated into BIGF either unchanged (PKFWHM = 1.0) or as a sharpened Gaussian having a FWHM scaled from the FWHM of the MEM peak by the multiplicative factor PKFWHM ($0 <$ PKFWHM $< 1$). To preclude differential blurring of any MEM peaks, set ASPIKE to a large number (e.g., 9999) and NSPIKE to 0. See references [1,3] for discussions of the pros and cons associated with creating BIGF by differentially blurring F.

NSPIKE is an integer specifying the number of relatively large, sharp features ("spikes") to be differentially blurred, not according to the defaults specified in the preceeding line, but as specified by the user in the following lines. This allows the user more control over the interpretation inherent in differential blurring F into BIGF. See reference [3].

If NSPIKE $> 0$, then NSPIKE lines are read next, each of which specifies the values of MNSPK,HOWSPK, DARSPK, and FWSPK for one "spike" to be differentially blurred in a customized fashion (not according to the defaults: ASPIKE, etc). MNSPK is the approximate mean (in log $\tau$) of the spike. For the MEM peak located near MNSPK, the parameters on this line override the defaults described two paragraphs above. For MEM peaks not specified by the NSPIKE value(s) of MNSPK, differential blurring will be performed according to the defaults above, when the default conditions are satisfied. HOWSPK = 1 specifies that this peak is well separated and should be replaced directly by a sharpened Gassian. HOWSPK = 2 specifies that this peak is to be "folded over" before being sharpened. HOWSPK = 3 specifies that a linear baseline is to be used to determine the fraction of the peak to be sharpened. DARSPK is the fraction of the area above the baseline initially determined to be sharpened (only relevant if HOWSPK = 3). FWSPK is the multiplicative factor ($0 <$ FWSPK $< 1$) used to narrow the FWHM of this peak. See reference [3] for a detailed discussion of these options.

When IBIGF = 1 (below), then BIGF is to be derived from F one or more times. If DCHI2 = 0, BIGF is revised only once when CHI2 reaches CHI2UP. If DCHI2 is greater than zero, BIGF is revised every time the fractional drop in C exceeds DCHI2, once CHI2 drops below CHI2UP.

DCHI2C is the fractional drop in CHI2 and in TAUC upon which the recommended lifetime distribution is chosen. DCHI2S is the fractional drop in CHI2 upon which the calculation is stopped.

FIXSUM is set to T (true) to fix the sum of discrete amplitudes during NLS fitting. FIXSUM = T is only supported for NDIST = 1 (monotonic kinetics) and for NLS fitting (Gaussian noise).

If WRNORM is set to T (true), the images written to files are normalized first. Note that the exponential amplitudes plotted in the PostScript summary of the discrete fits reflect this normalization when WRNORM is T. If WRNORM = T (true) in allin1 mode, the data's normalization, D0, is iteratively refined along with the error estimates: D0 is adjusted for use in the next MEM inversion according to the normalization of the distribution recommended by MemExp in the current inversion. When the normalization D0 is known (not to be changed by the program), set WRNORM to F (false) to avoid this adjustment of D0 in allin1 mode.

IBIGF should be set to 0, 1, 2, 3, or 4. IBIGF = 0 calls for a constant value of BIGF to be used throughout (a uniform 'prior' model). If IBIGF = 1, BIGF is to be revised by blurring F, either uniformly or differentially. If IBIGF = 2 the BIGF distribution is to be read from a file named on the following (final) line. If IBIGF = 3 or 4, the cutoff LTCUT is read from the following (final) line. When IBIGF = 3, peaks resolved in F at log $\tau <$ LTCUT are not represented in BIGF, and the remaining (slower) processes resolved in F are incorporated into BIGF as Gaussians with each of their widths multiplied by BFFWHM. When IBIG = 4, BIGF is taken to be a single Gaussian of unit area with mean and variance equal to those of the F distribution calculated only at values of log $\tau >=$ LTCUT. When deconvolving an IRF with a light-scattering correction, the value of $\xi$ in the 'prior' model is increased by an amount proportional to the relevant area of F at log $\tau <$ LTCUT (when IBIGF = 3 or 4). In this way, the calculation is resumed with density in F at short lifetimes reinterpreted as scattered excitation light.