MPFIT for Nim - Non-linear least squares fitting

Wrapper of the cMPFIT library for Nim.

Documentation

For the (mostly autogenerated) documentation, see:

https://vindaar.github.io/nim-mpfit/

Example

The following example shows the basic usage of the library. Note that the actual code related to mpfit-nim is only the definition of proc expH(...) and proc fitHalfLife(...).

In addition an echoResult proc is defined, which can be used to pretty print the result of the fit. Alternatively, you can call pretty taking the final parameters and the mp_result (i.e. the tuple returned by fit as separate arguments) to get a string representation of it.

Note: we use the mpfit/plotting submodule here to generate a plot in one line. It depends on ggplotnim.

import std / [strutils, sequtils, strformat]
import pkg / [zero_functional, seqmath]
import mpfit

const
  filename = "data/half_life_muon.txt"

func expH(p: seq[float], x: float): float =
  ## the function we'd like to fit. Any user defined function needs to
  ## be of the signature
  ## `proc[T](p: seq[T], x: T): T`
  ## i.e. conform to the `FuncProto[T]` type defined in mpfit_nim
  result = p[0] * exp(-p[1] * x)
  
proc parseHalfLifeData(filename: string): (seq[float], seq[float]) =
  ## Parse the input file. First create seq of tuple of floats
  ## then convert that to tuple of seq[float]
  let s = readFile(filename).splitLines --> filter('#' notin it and it.len > 0).
                                             map(it.splitWhitespace).
                                             map((it[0].parseFloat,
                                                  it[1].parseFloat))
  result[0] = s --> map(it[0])
  result[1] = s --> map(it[1])  
                                            
proc fitHalfLife(bins, counts, countsErr: seq[float]): (seq[float], mp_result) =
  ## the actual code which performs the fitting. Call the `fit` proc
  ## with the user defined function to be fitted as the first argument,
  ## the initial parameter guess as the second and finally x, y and y_err
  # start parameters
  let p = [1400.0, 1.0]
  # now just call fit
  let (pRes, res) = fit(expH, p, bins, counts, countsErr)
  echoResult(pRes, res = res)
  result = (pRes, res)
  echo &"The lifetime of the muon is ~ {1.0 / pRes[1]:.2f} µs"

when isMainModule:
  # first parse the data from the file
  let (bins, counts) = parseHalfLifeData(filename)
  # calculates errors: poisson errors on the counts
  let countsErr = counts.mapIt(sqrt(it))
  # perform the fit and echo results
  let (pRes, res) = fitHalfLife(bins, counts, countsErr)

  # plot the data and the fit
  import mpfit / plotting # import plotting convenience function
  plot(
    expH, pRes, # the function we fit and resulting fit params
    bins, counts, countsErr, # the input data & errors
    res, # the `mp_result` returned from the `fit` call
    xMin = 0.0, xMax = 10.0, # customize range of plot
    xlabel = "time / μs", ylabel = "# counts", title = "Muon half life measurement", # and labels
    outfile = "../media/muon_lifetime_measurement.png", # save as png
    verbose = false # we set `verbose` to false, as we already `echoResult` manually
  )

which outputs the following:

  χ²      = 74.0799    (8 DOF)
  χ²/dof  = 9.25999
  NPAR    = 2
  NFREE   = 2
  NPEGGED = 0
  NITER   = 11
  NFEV    = 33
  P[0] = 1937.1   +/- 45.6327
  P[1] = 0.515508 +/- 0.00876106
The lifetime of the muon is ~ 1.94 µs

and creates this plot:

Dependencies & Installation

The library depends on the cMPFIT library as a shared object. Either get the source code from here or use the code in c_src/. Compile the C library as follows:

gcc -c -Wall -Werror -fpic mpfit.c mpfit.h
gcc -shared -o libmpfit.so mpfit.o

which should create a libmpfit.so file in the same directory. The Nim library will link against it. Either copy the shared library to the location of your Nim code, in which you use mpfit-nim, or install it system wide, depending on your system it may look like the following (Ubuntu x64):

cp libmpfit.so /usr/lib/x86_64-linux-gnu

Once the shared library is available, there shouldn’t be anything else to do to use the library. Note: the example given in examples/fit_half_life.nim requires a few Nim libraries, which are not dependencies, since they are only used in the example, notably:

• seqmath (for linspace)
• plotly (+ chroma) (to plot the data and fit)
• zero_functional (to parse the data)

Usage

The library consists of a single exported fit procedure, which has the following signature:

proc fit*[T](f: FuncProto[T], pS: openArray[T], x, y, ey: openArray[T]): (seq[T], mp_result) =

the first argument is a user defined function (see below), the following arguments are:

• pS: the first guess for the parameters
• x: data for x
• y: data for y
• ey: errors for y

Note: currently the ey may not be an empty sequence, nor 0, since we use it as a weight. (TODO: change that!)

The mp_result object contains the chi^2 values for the fit, the errors on the parameters and additional information about the internal fitting process (e.g. number of times the user defined function was called). The type is defined in src/wrapper/mpfit_wrapper.nim.

The FuncProto[T] type is the following:

proc [T](p: seq[T], x: T): T

defined in src/mpfit_nim.nim. The user defined function needs to conform to that (see the example above).

License

The C code is governed by the licence as shown in c_src/DISCLAIMER. The Nim code is published under the MIT license.