Coder Social home page Coder Social logo

scinim / numericalnim Goto Github PK

View Code? Open in Web Editor NEW
88.0 9.0 9.0 664 KB

A collection of numerical methods written in Nim

Home Page: https://scinim.github.io/numericalnim/

License: MIT License

Nim 100.00%
nim nim-lang ode ode-solver integration integral numerical-methods numerical-integration numerical-computation ode-integrators

numericalnim's Introduction

NumericalNim

Monthly Test

NumericalNim is a collection of numerical methods written in Nim. Currently it has support for integration, optimization, curve-fitting, interpolation and ODEs. It can operate on floats and custom structures, such as vectors and tensors (if they support a set of operators).

Table of contents

Documentation & Tutorials

  • NumericalNim's documentation
  • SciNim's getting-started site has a bunch of tutorials. Specific tutorial are linked below in their respective sections.

Installation

Install NumericalNim using Nimble:

nimble install numericalnim

ODE

Tutorials

The integrators

These are the implemented ODE integrators:

  • rk21 - Heun's Adaptive 2nd order method.
  • BS32 - Bogacki–Shampine 3rd order adaptive method.
  • DOPRI54 - Dormand & Prince's adaptive 5th order method.
  • Heun2 - Heun's 2nd order fixed timestep method.
  • Ralston2 - Ralston's 2nd order fixed timestep method.
  • Kutta3 - Kutta's 3rd order fixed timestep method.
  • Heun3 - Heuns's 3rd order fixed timestep method.
  • Ralston3 - Ralston's 3rd order fixed timestep method.
  • SSPRK3 - Strong Stability Preserving Runge-Kutta 3rd order fixed timestep method.
  • Ralston4 - Ralston's 4th order fixed timestep method.
  • Kutta4 - Kutta's 4th order fixed timestep method.
  • RK4 - The standard 4th order, fixed timestep method we all know and love.
  • Tsit54 - Tsitouras adaptive 5th order method.
  • Vern65 - Verner's "most efficient" 6th order adaptive timestep method. (https://www.sfu.ca/~jverner/)

Dense Output

All integrators support dense output using a 3rd order Hermite interpolant. Method specific interpolants may be added in the future.

Integration

Tutorials

The methods

  • adaptiveGauss - Uses Gauss-Kronrod Quadrature to adaptivly integrate function. This is the recommended proc to use! Typically this is both the fastest and most accurate of the methods. This is the only one of the methods that supports using Inf and -Inf as integration limits.
  • gaussQuad - Uses Gauss-Legendre Quadrature to integrate functions. Choose between 20 different accuracies by setting how many function evaluations should be made on each subinterval with the nPoints parameter (1 - 20 is valid options).
  • trapz - Uses the trapezoidal rule to integrate both functions and discrete points. 2nd order method.
  • simpson - Uses Simpson's rule to integrate both functions and discrete points. 4th order method.
  • adaptiveSimpson - Uses a adaptive Simpson's rule to subdivide the integration interval in finer pieces where the integrand is changing a lot and wider pieces in intervals where it doesn't change much. This allows it to perform the integral efficiently and still have accuracy. The error is approximated using Richardson extrapolation. So technically the values it outputs are actually not Simpson's, but Booles' method.
  • romberg - Uses Romberg integration to integrate both functions and discrete points. Note: If discrete points are provided they must be equally spaced and the number of points must be of the form 2^k + 1 ie 3, 5, 9, 17, 33, 65, 129 etc.
  • cumtrapz - Uses the trapezoidal rule to integrate both functions and discrete points but it outputs a seq of the integral values at provided x-values.
  • cumsimpson - Uses Simpson's rule to integrate both functions and discrete points but it outputs a seq of the integral values at provided x-values.

Optimization

Tutorials:

One Dimensional optimization methods

  • steepest_descent - Basic method for local minimum finding.
  • conjugate_gradient - iterative implementation of solving Ax = b
  • newtons - Newton-Raphson implementation for 1D functions.
  • secant - The secant method for 1D functions.

Multidimensional optimization methods

  • lbfgs: Limited-memory BFGS, a lighter version of BFGS. This is the recommended method for larger problems.
  • bfgs: Broyden–Fletcher–Goldfarb–Shanno algorithm for optimization. Quasi-Newton method.
  • newton: The classic Newton method. Fast for small problems but struggles with bigger ones.
  • steepestDescent: The basic gradient descent method. Only use this if you have specific reasons to, as it is slower than the others for almost all problems.

Curve fitting

Tutorials

Levenberg-Marquardt

Levenberg-Marquardt (levmarq) is an algorithm for solving non-linear least squares problems, like curve fitting. You provide the data points (x- & y-values), and optionally the errors in your data, along with a function describing the curve you want to fit and a starting guess of the parameters. levmarq will then find a (locally) optimal solution that minimizes the squared error (or the χ² if you provided errors). The uncertainties of the parameters can be approximated using paramUncertainties.

Interpolation

Tutorials

Natural Cubic Splines

Cubic splines are piecewise polynomials of degree 3 ie. it is defined differently on different intervals. It passes through all the supplied points and has a continuos derivative. To find which interval a certain x-value is in we use a binary search algorithm to find it instead of looping over all intervals one after another.

Cubic Hermite Splines

Cubic Hermite Splines are piecewise polynomials of degree 3, the same as Natural Cubic Splines. The difference is that we can pass the the derivative at each point as well as the function value. If the derivatives are not passed, a three-point finite difference will be used instead but this will not give as accurate results compared to with derivatives. It may be better to use Natural Cubic Splines then instead. The advantage Hermite Splines have over Natural Cubic Splines in NumericalNim is that it can handle other types of y-values than floats. For example if you want to interpolate data (dependent on one variable) in Arraymancer Tensors you can do it by passing those as a seq[Tensor]. Hermite Splines' main mission in NumericalNim is to interpolate data points you get from solving ODEs as both the function value and the derivative is known at every point in time.

Radial basis function interpolation

Radial basis function interpolation conceptually works by placing a Gaussian at every data point in the dataset and scale them all such that the linear combination of them interpolates the values in each point. This means that this method don't have to construct a mesh. It also means the interpolant will be smooth as it simply is a sum of smooth Gaussians. Internally a partition of unity method is employed though to improve performance and numerical stability. Simply put, the domain is divided into patches and only points in neighboring patches affect each other.

Utils

I have included a few handy tools in numericalnim/utils.

Documentation

linspace

linspace creates a seq of N evenly spaced points between two numbers x1 and x2.

echo linspace(0.0, 10.0, 11)
@[0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]

arange

arange creates a seq of float between x1 and x2 separated by dx. You can choose to include or exclude the start- and endpoint (unless the steps goes exactly to x2).

echo arange(0.0, 5.0, 0.5)
echo arange(0.0, 4.9, 0.5)
@[0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0]
@[0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5]

meshgrid

Meshgrid allows you to create a Tensor containing the points in an N-dimensional grid.

chi2

Calculates the fit of a curve to data points.

Tips and tricks

Suggested compilation flags

There are many optimizations the C/C++ compiler can do but doesn't do by default. These can be activated by passing flags to the underlying compiler with -t:<flag>. These can be both of the SIMD variety where multiple operations are done at the same time, or fast-math which breaks the standards of floating point arithmetics to exchange a small bit of accuracy for performance.

To enable these, nim has to pass some flags to the C compiler. Below is a small table, listing the flags to add when compiling with nim c for some widely used compilers:

Compiler Flags
clang -t:-march=native -t:-ffast-math
gcc -t:-march=native -t:-ffast-math
icc -t:-march=core-avx2 -t:-fast
msvc -t:arch:AVX2 -t:fp:fast

Using enums as keys for NumContext

If you want to avoid KeyErrors regarding mistyped keys, you can use enums for the keys instead. The enum value will be converted to a string internally so there is no constraint that all keys must be from the same enum. Here is one example of how to use it:

type
  MyKeys = enum
    key1
    key2
    key3
var ctx = newNumContext[float]()
ctx[key1] = 3.14
ctx[key2] = 6.28

Status

This is a hobby project of mine, so please keep that in mind. I have tried to cover most of the use-cases in the tests but there's always the risk that I have missed something. If you spot a bug (or something worse) please open an Issue. If you want to help improve this project I would very much appreciate it.

numericalnim's People

Contributors

barroff avatar clonkk avatar hugogranstrom avatar joe-shyft avatar narimiran avatar ringabout avatar vindaar avatar

Stargazers

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

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

numericalnim's Issues

interpolate.nim Error: cannot instantiate: 'T'

Many thanks for this very interesting package.
I've tried the examples in your documentation.
The following CubicSpline example does not compile on my Nim (devel)

import numericalnim
var X = @[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0] # seconds
var Y = @[0.0, 4.0, 6.0, 6.5, 6.4, 6.2, 0.0] # meters/second
let spline = newCubicSpline(X, Y)
echo spline.eval(5.0)
#[
NumNim_BUG.nim(6, 6) template/generic instantiation from here
$HOME/.nimble/pkgs/numericalnim-0.6.0.1/numericalnim/interpolate.nim(182, 27) Error: cannot instantiate: 'T'
]#

Nim Compiler Version 1.5.1 [Linux: amd64]
Compiled at 2021-01-07
Copyright (c) 2006-2020 by Andreas Rumpf

git hash: 89a21e4ec71e705833d2aacd069e291cf41a19c6
active boot switches: -d:release

Barycentric interpolator produce a codegen error with --gc:arc

var ptsCoord = newSeq[array[2, float]]()
# corners
ptsCoord.add([-200.0,-150.0])
ptsCoord.add([-200.0,150.0])
ptsCoord.add([200.0,-150.0])
ptsCoord.add([200.0,150.0])
var ptsValue = @[0.0,1.0,4.0,5.0]
# extra pts
ptsCoord.add([0.0,-150.0])
ptsCoord.add([0.0,150.0])
ptsCoord.add([0.0,0.0])
ptsValue.add(8.0)
ptsValue.add(1.0)
ptsValue.add(1.0)
let bary = newBarycentric2D(ptsCoord.toTensor(), ptsValue.toTensor())

Compiled with nim cpp --gc:arc with Nim 1.6.0

It seems related to Vector2 in cdt

Improve documentation

Now that we finally have a documentation you can really see what procs are documented and which are not.

It could also be nice with a paragraph at the top of each module as well to describe what it provides and recommended methods.

  • Differentiate
  • integrate
  • interpolate
  • ode
  • optimize
  • rbf
  • utils

Check for duplicate x-values when sorting dataset

If there are duplicates in x-value in a dataset there are two scenarios:

  1. If the y-values are the same, remove all but one of the duplicates
  2. If the y-values aren't the same, raise exception.

While we are at it we should also revamp sortDataset to accept a variable number of input arrays. Should we use array[N: static int, seq[T]] as both input and output then? Nah I'll just do varargs[seq[T]] and sort along the first sequence only and keep track of the indices for the rest. That should be efficient enough. T prepare for VectorLike I should probably not use varargs though, will go with openarray.

Improve error handling and error messages

Stage 1: Improve the error messages (specifically trimAndSortDataset)

Stage 2: improve error handling. Some universal @SciNim way of reporting back non-fatal error to the user. For example that the maximum number of iterations was reached.

Update README

  • Remove recent major changes, they aren't very recent anymore.
  • The Gitter channels aren't really relevant either. Link to Science channel on Matrix and Discord instead.
  • Add link to scinim/getting-started

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.