Giter Site home page Giter Site logo

vindaar / nim-mpfit Goto Github PK

View Code? Open in Web Editor NEW
17.0 3.0 1.0 198 KB

A wrapper for the cMPFIT library for the Nim programming language, https://vindaar.github.io/nim-mpfit/

License: MIT License

Makefile 0.27% C 74.70% Nim 25.02%
nim nonlinear least-squares fitting levenberg-marquardt non-linear-optimization

nim-mpfit's Introduction

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: media/muon_lifetime_measurement.png

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.

nim-mpfit's People

Contributors

vindaar avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar

Forkers

yaohuic

nim-mpfit's Issues

add ability to fit with uncertainty in x

Add the "effective variance" method to deal with both uncertainties on y as well as x following Jay Orear:
https://aapt.scitation.org/doi/abs/10.1119/1.12972

This is also what's used in gnuplot 5. It may be somewhat worse than a total least squares approach, but it's easy to implement. See gnuplots implementation here:
https://sourceforge.net/p/gnuplot/gnuplot-main/ci/master/tree/src/fit.c#l529

Additional explanation:
https://ned.ipac.caltech.edu/level5/Sept01/Orear/Orear21.html

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.