• MPFIT for Nim - Non-linear least squares fitting Wrapper of the [[https://www.physics.wisc.edu/~craigm/idl/cmpfit.html][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 [[https://github.com/Vindaar/ggplotnim][ggplotnim]].

#+BEGIN_SRC nim :tangle examples/fit_half_life.nim 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 ) #+END_SRC

which outputs the following: #+BEGIN_SRC sh χ² = 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 #+END_SRC

and creates this plot: [[file:media/muon_lifetime_measurement.png]]

** Dependencies & Installation

The library depends on the cMPFIT library as a shared object. Either get the source code from [[https://www.physics.wisc.edu/~craigm/idl/cmpfit.html][here]] or use the code in [[file:c_src/]]. Compile the C library as follows: #+BEGIN_SRC sh gcc -c -Wall -Werror -fpic mpfit.c mpfit.h gcc -shared -o libmpfit.so mpfit.o #+END_SRC 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): #+BEGIN_SRC cp libmpfit.so /usr/lib/x86_64-linux-gnu #+END_SRC

Once the shared library is available, there shouldn't be anything else to do to use the library. Note: the example given in [[file: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: #+BEGIN_SRC nim proc fit*[T](f: FuncProto[T], pS: openArray[T], x, y, ey: openArray[T]): (seq[T], mp_result) = #+END_SRC 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 [[file:src/wrapper/mpfit_wrapper.nim]].

The =FuncProto[T]= type is the following: #+BEGIN_SRC nim proc [T](p: seq[T], x: T): T #+END_SRC defined in [[file: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 [[file:c_src/DISCLAIMER]]. The Nim code is published under the MIT license.