diff --git a/Demo.hs b/Demo.hs new file mode 100644 --- /dev/null +++ b/Demo.hs @@ -0,0 +1,994 @@ +-- This module is a Haskell translation of lmdemo.c from the C levmar library. + +module Demo where + +import LevMar ( levmar + + , Model + , Jacobian + + , Options(..), defaultOpts + + , LinearConstraints, noLinearConstraints + + , LevMarError + + , Info, CovarMatrix + + , S, Z + , SizedList(..) + ) + +import qualified LevMar.Fitting as Fitting + +type Result n = Either LevMarError + ( SizedList n Double + , Info Double + , CovarMatrix n Double + ) + +-------------------------------------------------------------------------------- +-- Handy type synonyms for type-level naturals: + +type N0 = Z +type N1 = S N0 +type N2 = S N1 +type N3 = S N2 +type N4 = S N3 +type N5 = S N4 + +-------------------------------------------------------------------------------- +-- Default options: + +opts :: Options Double +opts = defaultOpts { optStopNormInfJacTe = 1e-15 + , optStopNorm2Dp = 1e-15 + , optStopNorm2E = 1e-20 + } + +-------------------------------------------------------------------------------- +-- Rosenbrock function, +-- global minimum at (1, 1) + +ros :: Model N2 Double +ros p0 p1 = replicate ros_n ((1.0 - p0)**2 + ros_d*m**2) + where + m = p1 - p0**2 + +ros_jac :: Jacobian N2 Double +ros_jac p0 p1 = replicate ros_n (p0d ::: p1d ::: Nil) + where + p0d = -2 + 2*p0 - 4*ros_d*m*p0 + p1d = 2*ros_d*m + m = p1 - p0**2 + +ros_d :: Double +ros_d = 105.0 + +ros_n :: Int +ros_n = 2 + +ros_params :: SizedList N2 Double +ros_params = -1.2 ::: 1.0 ::: Nil + +ros_samples :: [Double] +ros_samples = replicate ros_n 0.0 + +run_ros :: Result N2 +run_ros = levmar ros + (Just ros_jac) + ros_params + ros_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Modified Rosenbrock problem, +-- global minimum at (1, 1) + +modros :: Model N2 Double +modros p0 p1 = [ 10*(p1 - p0**2) + , 1.0 - p0 + , modros_lam + ] + +modros_jac :: Jacobian N2 Double +modros_jac p0 _ = [ -20*p0 ::: 10.0 ::: Nil + , -1.0 ::: 0.0 ::: Nil + , 0.0 ::: 0.0 ::: Nil + ] + +modros_lam :: Double +modros_lam = 1e02 + +modros_n :: Int +modros_n = 3 + +modros_params :: SizedList N2 Double +modros_params = -1.2 ::: 1.0 ::: Nil + +modros_samples :: [Double] +modros_samples = replicate modros_n 0.0 + +run_modros :: Result N2 +run_modros = levmar modros + (Just modros_jac) + modros_params + modros_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Powell's function, +-- minimum at (0, 0) + +powell :: Model N2 Double +powell p0 p1 = [ p0 + , 10.0*p0 / m + 2*p1**2 + ] + where + m = p0 + 0.1 + +powell_jac :: Jacobian N2 Double +powell_jac p0 p1 = [ 1.0 ::: 0.0 ::: Nil + , 1.0 / m**2 ::: 4.0*p1 ::: Nil + ] + where + m = p0 + 0.1 + +powell_n :: Int +powell_n = 2 + +powell_params :: SizedList N2 Double +powell_params = -1.2 ::: 1.0 ::: Nil + +powell_samples :: [Double] +powell_samples = replicate powell_n 0.0 + +run_powell :: Result N2 +run_powell = levmar powell + (Just powell_jac) + powell_params + powell_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Wood's function, +-- minimum at (1, 1, 1, 1) + +wood :: Model N4 Double +wood p0 p1 p2 p3 = [ 10.0*(p1 - p0**2) + , 1.0 - p0 + , sqrt 90.0*(p3 - p2**2) + , 1.0 - p2 + , sqrt 10.0*(p1 + p3 - 2.0) + , (p1 - p3) / sqrt 10.0 + ] + +wood_n :: Int +wood_n = 6 + +wood_params :: SizedList N4 Double +wood_params = -3.0 ::: -1.0 ::: -3.0 ::: -1.0 ::: Nil + +wood_samples :: [Double] +wood_samples = replicate wood_n 0.0 + +run_wood :: Result N4 +run_wood = levmar wood + Nothing + wood_params + wood_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Meyer's (reformulated) data fitting problem, +-- minimum at (2.48, 6.18, 3.45) + +meyer :: Fitting.Model N3 Double Double +meyer p0 p1 p2 x = p0*exp (10.0*p1 / (ui + p2) - 13.0) + where + ui = 0.45 + 0.05*x + +meyer_jac :: Fitting.Jacobian N3 Double Double +meyer_jac p0 p1 p2 x = tmp + ::: 10.0*p0*tmp / (ui + p2) + ::: -10.0*p0*p1*tmp / ((ui + p2)*(ui + p2)) + ::: Nil + where + tmp = exp (10.0*p1 / (ui + p2) - 13.0) + ui = 0.45 + 0.05*x + +meyer_n :: Int +meyer_n = 16 + +meyer_params :: SizedList N3 Double +meyer_params = 8.85 ::: 4.0 ::: 2.5 ::: Nil + +meyer_samples :: [(Double, Double)] +meyer_samples = zip [0..] [ 34.780 + , 28.610 + , 23.650 + , 19.630 + , 16.370 + , 13.720 + , 11.540 + , 9.744 + , 8.261 + , 7.030 + , 6.005 + , 5.147 + , 4.427 + , 3.820 + , 3.307 + , 2.872 + ] + +run_meyer_jac :: Result N3 +run_meyer_jac = Fitting.levmar meyer + (Just meyer_jac) + meyer_params + meyer_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +run_meyer :: Result N3 +run_meyer = Fitting.levmar meyer + Nothing + meyer_params + meyer_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- helical valley function, +-- minimum at (1.0, 0.0, 0.0) + +helval :: Model N3 Double +helval p0 p1 p2 = [ 10.0*(p2 - 10.0*theta) + , 10.0*sqrt tmp - 1.0 + , p2 + ] + where + m = atan (p1 / p0) / (2.0*pi) + + tmp = p0**2 + p1**2 + + theta | p0 < 0.0 = m + 0.5 + | 0.0 < p0 = m + | p1 >= 0 = 0.25 + | otherwise = -0.25 + +heval_jac :: Jacobian N3 Double +heval_jac p0 p1 _ = [ 50.0*p1 / (pi*tmp) ::: -50.0*p0 / (pi*tmp) ::: 10.0 ::: Nil + , 10.0*p0 / sqrt tmp ::: 10.0*p1 / sqrt tmp ::: 0.0 ::: Nil + , 0.0 ::: 0.0 ::: 1.0 ::: Nil + ] + where + tmp = p0**2 + p1**2 + +helval_n :: Int +helval_n = 3 + +helval_params :: SizedList N3 Double +helval_params = -1.0 ::: 0.0 ::: 0.0 ::: Nil + +helval_samples :: [Double] +helval_samples = replicate helval_n 0.0 + +run_helval :: Result N3 +run_helval = levmar helval + (Just heval_jac) + helval_params + helval_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Boggs - Tolle problem 3 (linearly constrained), +-- minimum at (-0.76744, 0.25581, 0.62791, -0.11628, 0.25581) +-- +-- constr1: p0 + 3*p1 = 0 +-- constr2: p2 + p3 - 2*p4 = 0 +-- constr3: p1 - p4 = 0 + +bt3 :: Model N5 Double +bt3 p0 p1 p2 p3 p4 = replicate bt3_n (t1**2 + t2**2 + t3**2 + t4**2) + where + t1 = p0 - p1 + t2 = p1 + p2 - 2.0 + t3 = p3 - 1.0 + t4 = p4 - 1.0 + +bt3_jac :: Jacobian N5 Double +bt3_jac p0 p1 p2 p3 p4 = replicate bt3_n ( 2.0*t1 + ::: 2.0*(t2 - t1) + ::: 2.0*t2 + ::: 2.0*t3 + ::: 2.0*t4 + ::: Nil + ) + where + t1 = p0 - p1 + t2 = p1 + p2 - 2.0 + t3 = p3 - 1.0 + t4 = p4 - 1.0 + +bt3_n :: Int +bt3_n = 5 + +bt3_params :: SizedList N5 Double +bt3_params = 2.0 ::: 2.0 ::: 2.0 :::2.0 ::: 2.0 ::: Nil + +bt3_samples :: [Double] +bt3_samples = replicate bt3_n 0.0 + +bt3_linear_constraints :: LinearConstraints N3 N5 Double +bt3_linear_constraints = ( (1.0 ::: 3.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil) + ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil) + ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil) + ::: Nil + , 0.0 ::: 0.0 ::: 0.0 ::: Nil + ) + +run_bt3 :: Result N5 +run_bt3 = levmar bt3 + (Just bt3_jac) + bt3_params + bt3_samples + 1000 + opts + Nothing + Nothing + (Just bt3_linear_constraints) + Nothing + +-------------------------------------------------------------------------------- +-- Hock - Schittkowski problem 28 (linearly constrained), +-- minimum at (0.5, -0.5, 0.5) +-- +-- constr1: p0 + 2*p1 + 3*p2 = 1 + +hs28 :: Model N3 Double +hs28 p0 p1 p2 = replicate hs28_n (t1**2 + t2**2) + where + t1 = p0 + p1 + t2 = p1 + p2 + +hs28_jac :: Jacobian N3 Double +hs28_jac p0 p1 p2 = replicate hs28_n ( 2.0*t1 + ::: 2.0*(t1 + t2) + ::: 2.0*t2 + ::: Nil + ) + where + t1 = p0 + p1 + t2 = p1 + p2 + +hs28_n :: Int +hs28_n = 3 + +hs28_params :: SizedList N3 Double +hs28_params = -4.0 ::: 1.0 ::: 1.0 ::: Nil + +hs28_samples :: [Double] +hs28_samples = replicate hs28_n 0.0 + +hs28_linear_constraints :: LinearConstraints N1 N3 Double +hs28_linear_constraints = ( ((1.0 ::: 2.0 ::: 3.0 ::: Nil) ::: Nil) + , 1.0 ::: Nil + ) + +run_hs28 :: Result N3 +run_hs28 = levmar hs28 + (Just hs28_jac) + hs28_params + hs28_samples + 1000 + opts + Nothing + Nothing + (Just hs28_linear_constraints) + Nothing + +-------------------------------------------------------------------------------- +-- Hock - Schittkowski problem 48 (linearly constrained), +-- minimum at (1.0, 1.0, 1.0, 1.0, 1.0) +-- +-- constr1: sum [p0, p1, p2, p3, p4] = 5 +-- constr2: p2 - 2*(p3 + p4) = -3 + +hs48 :: Model N5 Double +hs48 p0 p1 p2 p3 p4 = replicate hs48_n (t1**2 + t2**2 + t3**2) + where + t1 = p0 - 1.0 + t2 = p1 - p2 + t3 = p3 - p4 + +hs48_jac :: Jacobian N5 Double +hs48_jac p0 p1 p2 p3 p4 = replicate hs48_n ( 2.0*t1 + ::: 2.0*t2 + ::: 2.0*t2 + ::: 2.0*t3 + ::: 2.0*t3 + ::: Nil + ) + where + t1 = p0 - 1.0 + t2 = p1 - p2 + t3 = p3 - p4 + +hs48_n :: Int +hs48_n = 3 + +hs48_params :: SizedList N5 Double +hs48_params = 3.0 ::: 5.0 ::: -3.0 ::: 2.0 ::: -2.0 ::: Nil + +hs48_samples :: [Double] +hs48_samples = replicate hs48_n 0.0 + +hs48_linear_constraints :: LinearConstraints N2 N5 Double +hs48_linear_constraints = ( (1.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: Nil) + ::: (0.0 ::: 0.0 ::: 1.0 ::: -2.0 ::: -2.0 ::: Nil) + ::: Nil + , 5.0 ::: -3.0 ::: Nil + ) + +run_hs48 :: Result N5 +run_hs48 = levmar hs48 + (Just hs48_jac) + hs48_params + hs48_samples + 1000 + opts + Nothing + Nothing + (Just hs48_linear_constraints) + Nothing + +-------------------------------------------------------------------------------- +-- Hock - Schittkowski problem 51 (linearly constrained), +-- minimum at (1.0, 1.0, 1.0, 1.0, 1.0) +-- +-- constr1: p0 + 3*p1 = 4 +-- constr2: p2 + p3 - 2*p4 = 0 +-- constr3: p1 - p4 = 0 + +hs51 :: Model N5 Double +hs51 p0 p1 p2 p3 p4 = replicate hs51_n (t1**2 + t2**2 + t3**2 + t4**2) + where + t1 = p0 - p1 + t2 = p1 + p2 - 2.0 + t3 = p3 - 1.0 + t4 = p4 - 1.0 + +hs51_jac :: Jacobian N5 Double +hs51_jac p0 p1 p2 p3 p4 = replicate hs51_n ( 2.0*t1 + ::: 2.0*(t2 - t1) + ::: 2.0*t2 + ::: 2.0*t3 + ::: 2.0*t4 + ::: Nil + ) + where + t1 = p0 - p1 + t2 = p1 + p2 - 2.0 + t3 = p3 - 1.0 + t4 = p4 - 1.0 + +hs51_n :: Int +hs51_n = 5 + +hs51_params :: SizedList N5 Double +hs51_params = 2.5 ::: 0.5 ::: 2.0 ::: -1.0 ::: 0.5 ::: Nil + +hs51_samples :: [Double] +hs51_samples = replicate hs51_n 0.0 + +hs51_linear_constraints :: LinearConstraints N3 N5 Double +hs51_linear_constraints = ( (1.0 ::: 3.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil) + ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil) + ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil) + ::: Nil + , 4.0 ::: 0.0 ::: 0.0 ::: Nil + ) + +run_hs51 :: Result N5 +run_hs51 = levmar hs51 + (Just hs51_jac) + hs51_params + hs51_samples + 1000 + opts + Nothing + Nothing + (Just hs51_linear_constraints) + Nothing + +-------------------------------------------------------------------------------- +-- Hock - Schittkowski problem 01 (box constrained), +-- minimum at (1.0, 1.0) +-- +-- constr1: p1 >= -1.5 + +hs01 :: Model N2 Double +hs01 p0 p1 = [ 10.0*(p1 - p0**2) + , 1.0 - p0 + ] + +hs01_jac :: Jacobian N2 Double +hs01_jac p0 _ = [ -20.0*p0 ::: 10.0 ::: Nil + , -1.0 ::: 0.0 ::: Nil + ] + +hs01_n :: Int +hs01_n = 2 + +hs01_params :: SizedList N2 Double +hs01_params = -2.0 ::: 1.0 ::: Nil + +hs01_samples :: [Double] +hs01_samples = replicate hs01_n 0.0 + +hs01_lb, hs01_ub :: SizedList N2 Double +hs01_lb = -_DBL_MAX ::: -1.5 ::: Nil +hs01_ub = _DBL_MAX ::: _DBL_MAX ::: Nil + +_DBL_MAX :: Double +_DBL_MAX = 1e+37 -- TODO: Get this directly from . + +run_hs01 :: Result N2 +run_hs01 = levmar hs01 + (Just hs01_jac) + hs01_params + hs01_samples + 1000 + opts + (Just hs01_lb) + (Just hs01_ub) + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Hock - Schittkowski MODIFIED problem 21 (box constrained), +-- minimum at (2.0, 0.0) +-- +-- constr1: 2 <= p0 <=50 +-- constr2: -50 <= p1 <=50 +-- +-- Original HS21 has the additional constraint 10*p0 - p1 >= 10 +-- which is inactive at the solution, so it is dropped here. + +hs21 :: Model N2 Double +hs21 p0 p1 = [p0 / 10.0, p1] + +hs21_jac :: Jacobian N2 Double +hs21_jac _ _ = [ 0.1 ::: 0.0 ::: Nil + , 0.0 ::: 1.0 ::: Nil + ] + +hs21_n :: Int +hs21_n = 2 + +hs21_params :: SizedList N2 Double +hs21_params = -1.0 ::: -1.0 ::: Nil + +hs21_samples :: [Double] +hs21_samples = replicate hs21_n 0.0 + +hs21_lb, hs21_ub :: SizedList N2 Double +hs21_lb = 2.0 ::: -50.0 ::: Nil +hs21_ub = 50.0 ::: 50.0 ::: Nil + +run_hs21 :: Result N2 +run_hs21 = levmar hs21 + (Just hs21_jac) + hs21_params + hs21_samples + 1000 + opts + (Just hs21_lb) + (Just hs21_ub) + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Problem hatfldb (box constrained), +-- minimum at (0.947214, 0.8, 0.64, 0.4096) +-- +-- constri: pi >= 0.0 (i=1..4) +-- constr5: p1 <= 0.8 + +hatfldb :: Model N4 Double +hatfldb p0 p1 p2 p3 = [ p0 - 1.0 + , p0 - sqrt p1 + , p1 - sqrt p2 + , p2 - sqrt p3 + ] + +hatfldb_jac :: Jacobian N4 Double +hatfldb_jac _ p1 p2 p3 = [ 1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil + , 1.0 ::: -0.5 / sqrt p1 ::: 0.0 ::: 0.0 ::: Nil + , 0.0 ::: 1.0 ::: -0.5 / sqrt p2 ::: 0.0 ::: Nil + , 0.0 ::: 0.0 ::: 1.0 ::: -0.5 / sqrt p3 ::: Nil + ] + +hatfldb_n :: Int +hatfldb_n = 4 + +hatfldb_params :: SizedList N4 Double +hatfldb_params = 0.1 ::: 0.1 ::: 0.1 ::: 0.1 ::: Nil + +hatfldb_samples :: [Double] +hatfldb_samples = replicate hatfldb_n 0.0 + +hatfldb_lb, hatfldb_ub :: SizedList N4 Double +hatfldb_lb = 0.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil +hatfldb_ub = _DBL_MAX ::: 0.8 ::: _DBL_MAX ::: _DBL_MAX ::: Nil + +run_hatfldb :: Result N4 +run_hatfldb = levmar hatfldb + (Just hatfldb_jac) + hatfldb_params + hatfldb_samples + 1000 + opts + (Just hatfldb_lb) + (Just hatfldb_ub) + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Problem hatfldc (box constrained), +-- minimum at (1.0, 1.0, 1.0, 1.0) +-- +-- constri: pi >= 0.0 (i=1..4) +-- constri+4: pi <= 10.0 (i=1..4) + +hatfldc :: Model N4 Double +hatfldc p0 p1 p2 p3 = [ p0 - 1.0 + , p0 - sqrt p1 + , p1 - sqrt p2 + , p3 - 1.0 + ] + +hatfldc_jac :: Jacobian N4 Double +hatfldc_jac _ p1 p2 _ = [ 1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil + , 1.0 ::: -0.5 / sqrt p1 ::: 0.0 ::: 0.0 ::: Nil + , 0.0 ::: 1.0 ::: -0.5 / sqrt p2 ::: 0.0 ::: Nil + , 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil + ] + +hatfldc_n :: Int +hatfldc_n = 4 + +hatfldc_params :: SizedList N4 Double +hatfldc_params = 0.9 ::: 0.9 ::: 0.9 ::: 0.9 ::: Nil + +hatfldc_samples :: [Double] +hatfldc_samples = replicate hatfldc_n 0.0 + +hatfldc_lb, hatfldc_ub :: SizedList N4 Double +hatfldc_lb = 0.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil +hatfldc_ub = 10.0 ::: 10.0 ::: 10.0 ::: 10.0 ::: Nil + +run_hatfldc :: Result N4 +run_hatfldc = levmar hatfldc + (Just hatfldc_jac) + hatfldc_params + hatfldc_samples + 1000 + opts + (Just hatfldc_lb) + (Just hatfldc_ub) + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Hock - Schittkowski (modified) problem 52 (box/linearly constrained), +-- minimum at (-0.09, 0.03, 0.25, -0.19, 0.03) +-- +-- constr1: p0 + 3*p1 = 0 +-- constr2: p2 + p3 - 2*p4 = 0 +-- constr3: p1 - p4 = 0 +-- +-- To the above 3 constraints, we add the following 5: +-- constr4: -0.09 <= p0 +-- constr5: 0.0 <= p1 <= 0.3 +-- constr6: p2 <= 0.25 +-- constr7: -0.2 <= p3 <= 0.3 +-- constr8: 0.0 <= p4 <= 0.3 + +modhs52 :: Model N5 Double +modhs52 p0 p1 p2 p3 p4 = [ 4.0*p0 - p1 + , p1 + p2 - 2.0 + , p3 - 1.0 + , p4 - 1.0 + ] + +modhs52_jac :: Jacobian N5 Double +modhs52_jac _ _ _ _ _ = [ 4.0 ::: -1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil + , 0.0 ::: 1.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: Nil + , 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: 0.0 ::: Nil + , 0.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil + ] + +modhs52_n :: Int +modhs52_n = 4 + +modhs52_params :: SizedList N5 Double +modhs52_params = 2.0 ::: 2.0 ::: 2.0 ::: 2.0 ::: 2.0 ::: Nil + +modhs52_samples :: [Double] +modhs52_samples = replicate modhs52_n 0.0 + +modhs52_linear_constraints :: LinearConstraints N3 N5 Double +modhs52_linear_constraints = ( (1.0 ::: 3.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil) + ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil) + ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil) + ::: Nil + , 0.0 ::: 0.0 ::: 0.0 ::: Nil + ) + +modhs52_weights :: SizedList N5 Double +modhs52_weights = 2000.0 ::: 2000.0 ::: 2000.0 ::: 2000.0 ::: 2000.0 ::: Nil + +modhs52_lb, modhs52_ub :: SizedList N5 Double +modhs52_lb = -0.09 ::: 0.0 ::: -_DBL_MAX ::: -0.2 ::: 0.0 ::: Nil +modhs52_ub = _DBL_MAX ::: 0.3 ::: 0.25 ::: 0.3 ::: 0.3 ::: Nil + +run_modhs52 :: Result N5 +run_modhs52 = levmar modhs52 + (Just modhs52_jac) + modhs52_params + modhs52_samples + 1000 + opts + (Just modhs52_lb) + (Just modhs52_ub) + (Just modhs52_linear_constraints) + (Just modhs52_weights) + +-------------------------------------------------------------------------------- +-- Schittkowski (modified) problem 235 (box/linearly constrained), +-- minimum at (-1.725, 2.9, 0.725) +-- +-- constr1: p0 + p2 = -1.0; +-- +-- To the above constraint, we add the following 2: +-- constr2: p1 - 4*p2 = 0 +-- constr3: 0.1 <= p1 <= 2.9 +-- constr4: 0.7 <= p2 + +mods235 :: Model N3 Double +mods235 p0 p1 _ = [ 0.1*(p0 - 1.0) + , p1 - p0**2 + ] + +mods235_jac :: Jacobian N3 Double +mods235_jac p0 _ _ = [ 0.1 ::: 0.0 ::: 0.0 ::: Nil + , -2.0*p0 ::: 1.0 ::: 0.0 ::: Nil + ] + +mods235_n :: Int +mods235_n = 2 + +mods235_params :: SizedList N3 Double +mods235_params = -2.0 ::: 3.0 ::: 1.0 ::: Nil + +mods235_samples :: [Double] +mods235_samples = replicate mods235_n 0.0 + +mods235_linear_constraints :: LinearConstraints N2 N3 Double +mods235_linear_constraints = ( (1.0 ::: 0.0 ::: 1.0 ::: Nil) + ::: (0.0 ::: 1.0 ::: -4.0 ::: Nil) + ::: Nil + , -1.0 ::: 0.0 ::: Nil + ) + +mods235_lb, mods235_ub :: SizedList N3 Double +mods235_lb = -_DBL_MAX ::: 0.1 ::: 0.7 ::: Nil +mods235_ub = _DBL_MAX ::: 2.9 ::: _DBL_MAX ::: Nil + +run_mods235 :: Result N3 +run_mods235 = levmar mods235 + (Just mods235_jac) + mods235_params + mods235_samples + 1000 + opts + (Just mods235_lb) + (Just mods235_ub) + (Just mods235_linear_constraints) + Nothing + +-------------------------------------------------------------------------------- +-- Boggs and Tolle modified problem 7 (box/linearly constrained), +-- minimum at (0.7, 0.49, 0.19, 1.19, -0.2) +-- +-- We keep the original objective function & starting point and use the +-- following constraints: +-- +-- subject to cons1: +-- x[1]+x[2] - x[3] = 1.0; +-- subject to cons2: +-- x[2] - x[4] + x[1] = 0.0; +-- subject to cons3: +-- x[5] + x[1] = 0.5; +-- subject to cons4: +-- x[5]>=-0.3; +-- subject to cons5: +-- x[1]<=0.7; + +modbt7 :: Model N5 Double +modbt7 p0 p1 _ _ _ = replicate modbt7_n (100.0*m**2 + n**2) + where + m = p1 - p0**2 + n = p0 - 1.0 + +modbt7_jac :: Jacobian N5 Double +modbt7_jac p0 p1 _ _ _ = replicate modbt7_n + ( -400.0*m*p0 + 2.0*p0 - 2.0 + ::: 200.0*m + ::: 0.0 + ::: 0.0 + ::: 0.0 + ::: Nil + ) + where + m = p1 - p0**2 + +modbt7_n :: Int +modbt7_n = 5 + +modbt7_params :: SizedList N5 Double +modbt7_params = -2.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: Nil + +modbt7_samples :: [Double] +modbt7_samples = replicate modbt7_n 0.0 + +modbt7_linear_constraints :: LinearConstraints N3 N5 Double +modbt7_linear_constraints = ( (1.0 ::: 1.0 ::: -1.0 ::: 0.0 ::: 0.0 ::: Nil) + ::: (1.0 ::: 1.0 ::: 0.0 ::: -1.0 ::: 0.0 ::: Nil) + ::: (1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil) + ::: Nil + , 1.0 ::: 0.0 ::: 0.5 ::: Nil + ) + +modbt7_lb, modbt7_ub :: SizedList N5 Double +modbt7_lb = -_DBL_MAX ::: -_DBL_MAX ::: -_DBL_MAX ::: -_DBL_MAX ::: -0.3 ::: Nil +modbt7_ub = 0.7 ::: _DBL_MAX ::: _DBL_MAX ::: _DBL_MAX ::: _DBL_MAX ::: Nil + +run_modbt7 :: Result N5 +run_modbt7 = levmar modbt7 + (Just modbt7_jac) + modbt7_params + modbt7_samples + 1000 + opts + (Just modbt7_lb) + (Just modbt7_ub) + (Just modbt7_linear_constraints) + Nothing + +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +-- !! TODO: This returns with: infStopReason = MaxIterations !! +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +-------------------------------------------------------------------------------- +-- Equilibrium combustion problem, constrained nonlinear equation from the book +-- by Floudas et al. +-- +-- Minimum at (0.0034, 31.3265, 0.0684, 0.8595, 0.0370) +-- +-- constri: pi>=0.0001 (i=1..5) +-- constri+5: pi<=100.0 (i=1..5) + +combust :: Model N5 Double +combust p0 p1 p2 p3 p4 = + [ p0*p1 + p0 - 3*p4 + , 2*p0*p1 + p0 + 3*r10*p1*p1 + p1*p2*p2 + r7*p1*p2 + r9*p1*p3 + r8*p1 - r*p4 + , 2*p1*p2*p2 + r7*p1*p2 + 2*r5*p2*p2 + r6*p2-8*p4 + , r9*p1*p3 + 2*p3*p3 - 4*r*p4 + , p0*p1 + p0 + r10*p1*p1 + p1*p2*p2 + r7*p1*p2 + r9*p1*p3 + r8*p1 + r5*p2*p2 + r6*p2 + p3*p3 - 1.0 + ] + +r, r5, r6, r7, r8, r9, r10 :: Double +r = 10 +r5 = 0.193 +r6 = 4.10622*1e-4 +r7 = 5.45177*1e-4 +r8 = 4.4975 *1e-7 +r9 = 3.40735*1e-5 +r10 = 9.615 *1e-7 + +combust_jac :: Jacobian N5 Double +combust_jac p0 p1 p2 p3 _ = + [ p1 + 1 + ::: p0 + ::: 0.0 + ::: 0.0 + ::: -3 + ::: Nil + , 2*p1 + 1 + ::: 2*p0 + 6*r10*p1 + p2*p2 + r7*p2 + r9*p3 + r8 + ::: 2*p1*p2 + r7*p1 + ::: r9*p1 + ::: -r + ::: Nil + , 0.0 + ::: 2*p2*p2 + r7*p2 + ::: 4*p1*p2 + r7*p1 + 4*r5*p2 + r6 + ::: 0.0 + ::: -8 + ::: Nil + , 0.0 + ::: r9*p3 + ::: 0.0 + ::: r9*p1 + 4*p3 + ::: -4*r + ::: Nil + , p1 + 1 + ::: p0 + 2*r10*p1 + p2*p2 + r7*p2 + r9*p3 + r8 + ::: 2*p1*p2 + r7*p1 + 2*r5*p2 + r6 + ::: r9*p1 + 2*p3 + ::: 0.0 + ::: Nil + ] + +combust_n :: Int +combust_n = 5 + +combust_params :: SizedList N5 Double +combust_params = 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: Nil + +combust_samples :: [Double] +combust_samples = replicate combust_n 0.0 + +combust_lb, combust_ub :: SizedList N5 Double +combust_lb = 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: Nil +combust_ub = 100.0 ::: 100.0 ::: 100.0 ::: 100.0 ::: 100.0 ::: Nil + +run_combust :: Result N5 +run_combust = levmar combust + (Just combust_jac) + combust_params + combust_samples + 1000 + opts + (Just combust_lb) + (Just combust_ub) + noLinearConstraints + Nothing + +-- The End --------------------------------------------------------------------- diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,32 @@ +Copyright (c) 2009 Roel van Dijk, Bas van Dijk + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * The name of Roel van Dijk and Bas van Dijk and the names of + contributors may NOT be used to endorse or promote products + derived from this software without specific prior written + permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/LevMar.hs b/LevMar.hs new file mode 100644 --- /dev/null +++ b/LevMar.hs @@ -0,0 +1,164 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE ScopedTypeVariables #-} + +-------------------------------------------------------------------------------- +-- | +-- Module : LevMar +-- Copyright : (c) 2009 Roel van Dijk & Bas van Dijk +-- License : BSD-style (see the file LICENSE) +-- +-- Maintainer : vandijk.roel@gmail.com, v.dijk.bas@gmail.com +-- Stability : Experimental +-- +-- +-- +-- For additional documentation see the documentation of the levmar C +-- library which this library is based on: +-- +-- +-------------------------------------------------------------------------------- + +module LevMar + ( -- * Model & Jacobian. + Model + , Jacobian + + -- * Levenberg-Marquardt algorithm. + , LMA_I.LevMarable + , levmar + + , LinearConstraints + , noLinearConstraints + , Matrix + + -- * Minimization options. + , LMA_I.Options(..) + , LMA_I.defaultOpts + + -- * Output + , LMA_I.Info(..) + , LMA_I.StopReason(..) + , CovarMatrix + + , LMA_I.LevMarError(..) + + -- *Type-level machinery + , Z, S, Nat + , SizedList(..) + , NFunction + ) + where + + +import qualified LevMar.Intermediate as LMA_I + +import TypeLevelNat (Z, S, Nat) +import SizedList (SizedList(..), toList, unsafeFromList) +import NFunction (NFunction, ($*)) + +import Data.Either + + +-------------------------------------------------------------------------------- +-- Model & Jacobian. +-------------------------------------------------------------------------------- + +{- | A function from @n@ parameters of type @r@ to a list of @r@. + +An example from /Demo.hs/: + +@ +type N4 = 'S' ('S' ('S' ('S' 'Z'))) + +hatfldc :: Model N4 Double +hatfldc p0 p1 p2 p3 = [ p0 - 1.0 + , p0 - sqrt p1 + , p1 - sqrt p2 + , p3 - 1.0 + ] +@ +-} +type Model n r = NFunction n r [r] + +{- | The jacobian of the 'Model' function. Expressed as a function from +@n@ parameters of type @r@ to a list of @n@-sized lists of @r@ + +See: + +For example the jacobian of the above @hatfldc@ model is: + +@ +type N4 = 'S' ('S' ('S' ('S' 'Z'))) + +hatfldc_jac :: Jacobian N4 Double +hatfldc_jac _ p1 p2 _ = [ 1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil + , 1.0 ::: -0.5 / sqrt p1 ::: 0.0 ::: 0.0 ::: Nil + , 0.0 ::: 1.0 ::: -0.5 / sqrt p2 ::: 0.0 ::: Nil + , 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil + ] +@ +-} + +type Jacobian n r = NFunction n r [SizedList n r] + + +-------------------------------------------------------------------------------- +-- Levenberg-Marquardt algorithm. +-------------------------------------------------------------------------------- + +-- | The Levenberg-Marquardt algorithm. +levmar :: forall n k r. (Nat n, Nat k, LMA_I.LevMarable r) + => (Model n r) -- ^ Model + -> Maybe (Jacobian n r) -- ^ Optional jacobian + -> SizedList n r -- ^ Initial parameters + -> [r] -- ^ Samples + -> Integer -- ^ Maximum number of iterations + -> LMA_I.Options r -- ^ Minimization options + -> Maybe (SizedList n r) -- ^ Optional lower bounds + -> Maybe (SizedList n r) -- ^ Optional upper bounds + -> Maybe (LinearConstraints k n r) -- ^ Optional linear constraints + -> Maybe (SizedList n r) -- ^ Optional weights + -> Either LMA_I.LevMarError (SizedList n r, LMA_I.Info r, CovarMatrix n r) + +levmar model mJac params ys itMax opts mLowBs mUpBs mLinC mWghts = + fmap convertResult $ LMA_I.levmar (convertModel model) + (fmap convertJacob mJac) + (toList params) + ys + itMax + opts + (fmap toList mLowBs) + (fmap toList mUpBs) + (fmap convertLinC mLinC) + (fmap toList mWghts) + where + convertModel f = \ps -> f $* (unsafeFromList ps :: SizedList n r) + convertJacob f = \ps -> map toList ((f $* (unsafeFromList ps :: SizedList n r)) :: [SizedList n r]) + convertLinC (cMat, rhcVec) = ( map toList $ toList cMat + , toList rhcVec + ) + convertResult (psResult, info, covar) = ( unsafeFromList psResult + , info + , unsafeFromList $ map unsafeFromList covar + ) + +-- | Linear constraints consisting of a constraints matrix, /kxn/ and +-- a right hand constraints vector, /kx1/ where /n/ is the number of +-- parameters and /k/ is the number of constraints. +type LinearConstraints k n r = (Matrix k n r, SizedList k r) + +-- |Value to denote the absense of any linear constraints over the +-- parameters of the model function. Use this instead of 'Nothing' +-- because the type parameter which contains the number of constraints +-- can't be inferred. +noLinearConstraints :: Nat n => Maybe (LinearConstraints Z n r) +noLinearConstraints = Nothing + +-- | A /nxm/ matrix is a sized list of /n/ sized lists of length /m/. +type Matrix n m r = SizedList n (SizedList m r) + +-- | Covariance matrix corresponding to LS solution. +type CovarMatrix n r = Matrix n n r + + +-- The End --------------------------------------------------------------------- diff --git a/LevMar/Fitting.hs b/LevMar/Fitting.hs new file mode 100644 --- /dev/null +++ b/LevMar/Fitting.hs @@ -0,0 +1,136 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +-------------------------------------------------------------------------------- +-- | +-- Module : LevMar.Fitting +-- Copyright : (c) 2009 Roel van Dijk & Bas van Dijk +-- License : BSD-style (see the file LICENSE) +-- +-- Maintainer : vandijk.roel@gmail.com, v.dijk.bas@gmail.com +-- Stability : Experimental +-- +-- This module provides the Levenberg-Marquardt algorithm specialised +-- for curve-fitting. +-- +-- For additional documentation see the documentation of the levmar C +-- library which this library is based on: +-- +-- +-------------------------------------------------------------------------------- + +module LevMar.Fitting + ( -- * Model & Jacobian. + Model + , Jacobian + + -- * Levenberg-Marquardt algorithm. + , LMA.LevMarable + , levmar + + , LMA.LinearConstraints + , LMA.noLinearConstraints + , LMA.Matrix + + -- * Minimization options. + , LMA.Options(..) + , LMA.defaultOpts + + -- * Output + , LMA.Info(..) + , LMA.StopReason(..) + , LMA.CovarMatrix + + , LMA.LevMarError(..) + + -- *Type-level machinery + , Z, S, Nat + , SizedList(..) + , NFunction + , ComposeN + ) where + + +import qualified LevMar as LMA + +import TypeLevelNat (Z, S, Nat, witnessNat) +import SizedList (SizedList) +import NFunction (NFunction, ComposeN, compose) + + +-------------------------------------------------------------------------------- +-- Model & Jacobian. +-------------------------------------------------------------------------------- + +{- | A function from @n@ parameters of type @r@ and an x-value of type +@a@ to a value of type @r@. + +For example, the quadratic function @f(x) = a*x^2 + b*x + c@ can be +written as: + +@ +type N3 = 'S' ('S' ('S' 'Z')) + +quad :: 'Num' r => 'Model' N3 r r +quad a b c x = a*x^2 + b*x + c +@ +-} +type Model n r a = NFunction n r (a -> r) + +{- | The jacobian of the 'Model' function. Expressed as a function from +@n@ parameters of type @r@ and an x-value of type @a@ to a vector +of @n@ values of type @r@. + +See: + +For example, the jacobian of the quadratic function @f(x) = a*x^2 + +b*x + c@ can be written as: + +@ +type N3 = 'S' ('S' ('S' 'Z')) + +quadJacob :: 'Num' r => 'Jacobian' N3 r r +quadJacob _ _ _ x = x^2 -- with respect to a + ::: x -- with respect to b + ::: 1 -- with respect to c + ::: 'Nil' +@ + +Notice you don't have to differentiate for @x@. +-} +type Jacobian n r a = NFunction n r (a -> SizedList n r) + + +-------------------------------------------------------------------------------- +-- Levenberg-Marquardt algorithm. +-------------------------------------------------------------------------------- + +-- | The Levenberg-Marquardt algorithm specialised for curve-fitting. +levmar :: forall n k r a. (Nat n, ComposeN n, Nat k, LMA.LevMarable r) + => (Model n r a) -- ^ Model + -> Maybe (Jacobian n r a) -- ^ Optional jacobian + -> SizedList n r -- ^ Initial parameters + -> [(a, r)] -- ^ Samples + -> Integer -- ^ Maximum number of iterations + -> LMA.Options r -- ^ Options + -> Maybe (SizedList n r) -- ^ Optional lower bounds + -> Maybe (SizedList n r) -- ^ Optional upper bounds + -> Maybe (LMA.LinearConstraints k n r) -- ^ Optional linear constraints + -> Maybe (SizedList n r) -- ^ Optional weights + -> Either LMA.LevMarError (SizedList n r, LMA.Info r, LMA.CovarMatrix n r) +levmar model mJac params samples = LMA.levmar (convertModel model) + (fmap convertJacob mJac) + params + ys + where + (xs, ys) = unzip samples + + convertModel :: Model n r a -> LMA.Model n r + convertModel = compose (witnessNat :: n) (undefined :: r) + (\(f :: a -> r) -> map f xs) + + convertJacob :: Jacobian n r a -> LMA.Jacobian n r + convertJacob = compose (witnessNat :: n) (undefined :: r) + (\(f :: a -> SizedList n r) -> map f xs) + + +-- The End --------------------------------------------------------------------- diff --git a/LevMar/Intermediate.hs b/LevMar/Intermediate.hs new file mode 100644 --- /dev/null +++ b/LevMar/Intermediate.hs @@ -0,0 +1,366 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE FlexibleInstances #-} + +-------------------------------------------------------------------------------- +-- | +-- Module : LevMar.Intermediate +-- Copyright : (c) 2009 Roel van Dijk & Bas van Dijk +-- License : BSD-style (see the file LICENSE) +-- +-- Maintainer : vandijk.roel@gmail.com, v.dijk.bas@gmail.com +-- Stability : Experimental +-- +-- +-- +-- For additional documentation see the documentation of the levmar C +-- library which this library is based on: +-- +-- +-------------------------------------------------------------------------------- + +module LevMar.Intermediate + ( -- * Model & Jacobian. + Model + , Jacobian + + -- * Levenberg-Marquardt algorithm. + , LevMarable + , levmar + + , LinearConstraints + + -- * Minimization options. + , Options(..) + , defaultOpts + + -- * Output + , Info(..) + , StopReason(..) + , CovarMatrix + + , LevMarError(..) + ) where + + +import Foreign.Marshal.Array (allocaArray, peekArray, pokeArray, withArray) +import Foreign.Ptr (Ptr, nullPtr, plusPtr) +import Foreign.Storable (Storable) +import Foreign.C.Types (CInt) +import System.IO.Unsafe (unsafePerformIO) +import Data.Maybe (fromJust, fromMaybe, isJust) +import Control.Monad.Instances -- for 'instance Functor (Either a)' + +import qualified Bindings.LevMar.CurryFriendly as LMA_C + + +-------------------------------------------------------------------------------- +-- Model & Jacobian. +-------------------------------------------------------------------------------- + +type Model r = [r] -> [r] + +-- | See: +type Jacobian r = [r] -> [[r]] + + +-------------------------------------------------------------------------------- +-- Levenberg-Marquardt algorithm. +-------------------------------------------------------------------------------- + +-- | The Levenberg-Marquardt algorithm is overloaded to work on 'Double' and 'Float'. +class LevMarable r where + + -- | The Levenberg-Marquardt algorithm. + levmar :: Model r -- ^ Model + -> Maybe (Jacobian r) -- ^ Optional jacobian + -> [r] -- ^ Initial parameters + -> [r] -- ^ Samples + -> Integer -- ^ Maximum iterations + -> Options r -- ^ Minimization options + -> Maybe [r] -- ^ Optional lower bounds + -> Maybe [r] -- ^ Optional upper bounds + -> Maybe (LinearConstraints r) -- ^ Optional linear constraints + -> Maybe [r] -- ^ Optional weights + -> Either LevMarError ([r], Info r, CovarMatrix r) + +instance LevMarable Float where + levmar = gen_levmar LMA_C.slevmar_der + LMA_C.slevmar_dif + LMA_C.slevmar_bc_der + LMA_C.slevmar_bc_dif + LMA_C.slevmar_lec_der + LMA_C.slevmar_lec_dif + LMA_C.slevmar_blec_der + LMA_C.slevmar_blec_dif + +instance LevMarable Double where + levmar = gen_levmar LMA_C.dlevmar_der + LMA_C.dlevmar_dif + LMA_C.dlevmar_bc_der + LMA_C.dlevmar_bc_dif + LMA_C.dlevmar_lec_der + LMA_C.dlevmar_lec_dif + LMA_C.dlevmar_blec_der + LMA_C.dlevmar_blec_dif + +{- | @gen_levmar@ takes the low-level C functions as arguments and +executes one of them depending on the optional jacobian and constraints. + +Preconditions: + length ys >= length ps + + isJust mLowBs && length (fromJust mLowBs) == length ps + && isJust mUpBs && length (fromJust mUpBs) == length ps + + boxConstrained && (all $ zipWith (<=) (fromJust mLowBs) (fromJust mUpBs)) +-} +gen_levmar :: forall cr r. (Storable cr, RealFrac cr, Real r, Fractional r) + => LMA_C.LevMarDer cr + -> LMA_C.LevMarDif cr + -> LMA_C.LevMarBCDer cr + -> LMA_C.LevMarBCDif cr + -> LMA_C.LevMarLecDer cr + -> LMA_C.LevMarLecDif cr + -> LMA_C.LevMarBLecDer cr + -> LMA_C.LevMarBLecDif cr + + -> Model r -- ^ Model + -> Maybe (Jacobian r) -- ^ Optional jacobian + -> [r] -- ^ Initial parameters + -> [r] -- ^ Samples + -> Integer -- ^ Maximum iterations + -> Options r -- ^ Options + -> Maybe [r] -- ^ Optional lower bounds + -> Maybe [r] -- ^ Optional upper bounds + -> Maybe (LinearConstraints r) -- ^ Optional linear constraints + -> Maybe [r] -- ^ Optional weights + -> Either LevMarError ([r], Info r, CovarMatrix r) +gen_levmar f_der + f_dif + f_bc_der + f_bc_dif + f_lec_der + f_lec_dif + f_blec_der + f_blec_dif + model mJac ps ys itMax opts mLowBs mUpBs mLinC mWeights + = unsafePerformIO $ + withArray (map realToFrac ps) $ \psPtr -> + withArray (map realToFrac ys) $ \ysPtr -> + withArray (map realToFrac $ optsToList opts) $ \optsPtr -> + allocaArray LMA_C._LM_INFO_SZ $ \infoPtr -> + allocaArray covarLen $ \covarPtr -> + LMA_C.withModel (convertModel model) $ \modelPtr -> do + + let runDif :: LMA_C.LevMarDif cr -> IO CInt + runDif f = f modelPtr + psPtr + ysPtr + (fromIntegral lenPs) + (fromIntegral lenYs) + (fromIntegral itMax) + optsPtr + infoPtr + nullPtr + covarPtr + nullPtr + + r <- case mJac of + Just jac -> LMA_C.withJacobian (convertJacobian jac) $ \jacobPtr -> + let runDer :: LMA_C.LevMarDer cr -> IO CInt + runDer f = runDif $ f jacobPtr + in if boxConstrained + then if linConstrained + then withBoxConstraints (withLinConstraints $ withWeights runDer) f_blec_der + else withBoxConstraints runDer f_bc_der + else if linConstrained + then withLinConstraints runDer f_lec_der + else runDer f_der + + Nothing -> if boxConstrained + then if linConstrained + then withBoxConstraints (withLinConstraints $ withWeights runDif) f_blec_dif + else withBoxConstraints runDif f_bc_dif + else if linConstrained + then withLinConstraints runDif f_lec_dif + else runDif f_dif + + if r < 0 + && r /= LMA_C._LM_ERROR_SINGULAR_MATRIX -- we don't treat these two as an error + && r /= LMA_C._LM_ERROR_SUM_OF_SQUARES_NOT_FINITE + then return $ Left $ convertLevMarError r + else do result <- peekArray lenPs psPtr + info <- peekArray LMA_C._LM_INFO_SZ infoPtr + + let covarPtrEnd = plusPtr covarPtr covarLen + let convertCovarMatrix ptr + | ptr == covarPtrEnd = return [] + | otherwise = do row <- peekArray lenPs ptr + rows <- convertCovarMatrix $ plusPtr ptr lenPs + return $ row : rows + + covar <- convertCovarMatrix covarPtr + + return $ Right ( map realToFrac result + , listToInfo info + , map (map realToFrac) covar + ) + where + lenPs = length ps + lenYs = length ys + covarLen = lenPs * lenPs + (cMat, rhcVec) = fromJust mLinC + + -- Whether the parameters are constrained by a linear equation. + linConstrained = isJust mLinC + + -- Whether the parameters are constrained by a bounding box. + boxConstrained = isJust mLowBs || isJust mUpBs + + withBoxConstraints f g = maybeWithArray ((fmap . fmap) realToFrac mLowBs) $ \lBsPtr -> + maybeWithArray ((fmap . fmap) realToFrac mUpBs) $ \uBsPtr -> + f $ g lBsPtr uBsPtr + + withLinConstraints f g = withArray (map realToFrac $ concat cMat) $ \cMatPtr -> + withArray (map realToFrac rhcVec) $ \rhcVecPtr -> + f $ g cMatPtr rhcVecPtr $ fromIntegral $ length cMat + + withWeights f g = maybeWithArray ((fmap . fmap) realToFrac mWeights) $ \weightsPtr -> + f $ g weightsPtr + +convertModel :: (Real r, Fractional r, Storable c, Real c, Fractional c) + => Model r -> LMA_C.Model c +convertModel model = \parPtr hxPtr numPar _ _ -> do + params <- peekArray (fromIntegral numPar) parPtr + pokeArray hxPtr $ map realToFrac $ model $ map realToFrac params + +convertJacobian :: (Real r, Fractional r, Storable c, Real c, Fractional c) + => Jacobian r -> LMA_C.Jacobian c +convertJacobian jac = \parPtr jPtr numPar _ _ -> do + params <- peekArray (fromIntegral numPar) parPtr + pokeArray jPtr $ concatMap (map realToFrac) $ jac $ map realToFrac params + +maybeWithArray :: Storable a => Maybe [a] -> (Ptr a -> IO b) -> IO b +maybeWithArray Nothing f = f nullPtr +maybeWithArray (Just xs) f = withArray xs f + + +-- | Linear constraints consisting of a constraints matrix, /kxm/ and +-- a right hand constraints vector, /kx1/ where /m/ is the number of +-- parameters and /k/ is the number of constraints. +type LinearConstraints r = ([[r]], [r]) + + +-------------------------------------------------------------------------------- +-- Minimization options. +-------------------------------------------------------------------------------- + +-- | Minimization options +data Options r = + Opts { optScaleInitMu :: r -- ^ Scale factor for initial /mu/. + , optStopNormInfJacTe :: r -- ^ Stopping thresholds for @||J^T e||_inf@. + , optStopNorm2Dp :: r -- ^ Stopping thresholds for @||Dp||_2@. + , optStopNorm2E :: r -- ^ Stopping thresholds for @||e||_2@. + , optDelta :: r -- ^ Step used in the difference approximation to the Jacobian. + -- If @optDelta<0@, the Jacobian is approximated + -- with central differences which are more accurate + -- (but slower!) compared to the forward differences + -- employed by default. + } deriving Show + +-- | Default minimization options +defaultOpts :: Fractional r => Options r +defaultOpts = Opts { optScaleInitMu = LMA_C._LM_INIT_MU + , optStopNormInfJacTe = LMA_C._LM_STOP_THRESH + , optStopNorm2Dp = LMA_C._LM_STOP_THRESH + , optStopNorm2E = LMA_C._LM_STOP_THRESH + , optDelta = LMA_C._LM_DIFF_DELTA + } + +optsToList :: Options r -> [r] +optsToList (Opts mu eps1 eps2 eps3 delta) = + [mu, eps1, eps2, eps3, delta] + + +-------------------------------------------------------------------------------- +-- Output +-------------------------------------------------------------------------------- + +-- | Information regarding the minimization. +data Info r = Info { infNorm2initE :: r -- ^ @||e||_2@ at initial parameters. + , infNorm2E :: r -- ^ @||e||_2@ at estimated parameters. + , infNormInfJacTe :: r -- ^ @||J^T e||_inf@ at estimated parameters. + , infNorm2Dp :: r -- ^ @||Dp||_2@ at estimated parameters. + , infMuDivMax :: r -- ^ @\mu/max[J^T J]_ii ]@ at estimated parameters. + , infNumIter :: Integer -- ^ Number of iterations. + , infStopReason :: StopReason -- ^ Reason for terminating. + , infNumFuncEvals :: Integer -- ^ Number of function evaluations. + , infNumJacobEvals :: Integer -- ^ Number of jacobian evaluations. + , infNumLinSysSolved :: Integer -- ^ Number of linear systems solved, i.e. attempts for reducing error. + } deriving Show + +listToInfo :: (RealFrac cr, Fractional r) => [cr] -> Info r +listToInfo [a,b,c,d,e,f,g,h,i,j] = + Info { infNorm2initE = realToFrac a + , infNorm2E = realToFrac b + , infNormInfJacTe = realToFrac c + , infNorm2Dp = realToFrac d + , infMuDivMax = realToFrac e + , infNumIter = floor f + , infStopReason = toEnum $ floor g - 1 + , infNumFuncEvals = floor h + , infNumJacobEvals = floor i + , infNumLinSysSolved = floor j + } +listToInfo _ = error "liftToInfo: wrong list length" + +-- | Reason for terminating. +data StopReason = SmallGradient -- ^ Stopped because of small gradient @J^T e@. + | SmallDp -- ^ Stopped because of small Dp. + | MaxIterations -- ^ Stopped because maximum iterations was reached. + | SingularMatrix -- ^ Stopped because of singular matrix. Restart from current estimated parameters with increased 'optScaleInitMu'. + | SmallestError -- ^ Stopped because no further error reduction is possible. Restart with increased 'optScaleInitMu'. + | SmallNorm2E -- ^ Stopped because of small @||e||_2@. + | InvalidValues -- ^ Stopped because model function returned invalid values (i.e. NaN or Inf). This is a user error. + deriving (Show, Enum) + +-- | Covariance matrix corresponding to LS solution. +type CovarMatrix r = [[r]] + + +-------------------------------------------------------------------------------- +-- Error +-------------------------------------------------------------------------------- + +data LevMarError + = LapackError -- ^ A call to a lapack subroutine failed in the underlying C levmar library. + | FailedBoxCheck -- ^ At least one lower bound exceeds the upper one. + | MemoryAllocationFailure -- ^ A call to @malloc@ failed in the underlying C levmar library. + | ConstraintMatrixRowsGtCols -- ^ The matrix of constraints cannot have more rows than columns. + | ConstraintMatrixNotFullRowRank -- ^ Constraints matrix is not of full row rank. + | TooFewMeasurements -- ^ Cannot solve a problem with fewer measurements than unknowns. + -- In case linear constraints are provided, this error is also returned + -- when the number of measurements is smaller than the number of unknowns + -- minus the number of equality constraints. + deriving Show + +levmarCErrorToLevMarError :: [(CInt, LevMarError)] +levmarCErrorToLevMarError = + [ (LMA_C._LM_ERROR_LAPACK_ERROR, LapackError) + --, (LMA_C._LM_ERROR_NO_JACOBIAN, can never happen) + --, (LMA_C._LM_ERROR_NO_BOX_CONSTRAINTS, can never happen) + , (LMA_C._LM_ERROR_FAILED_BOX_CHECK, FailedBoxCheck) + , (LMA_C._LM_ERROR_MEMORY_ALLOCATION_FAILURE, MemoryAllocationFailure) + , (LMA_C._LM_ERROR_CONSTRAINT_MATRIX_ROWS_GT_COLS, ConstraintMatrixRowsGtCols) + , (LMA_C._LM_ERROR_CONSTRAINT_MATRIX_NOT_FULL_ROW_RANK, ConstraintMatrixNotFullRowRank) + , (LMA_C._LM_ERROR_TOO_FEW_MEASUREMENTS, TooFewMeasurements) + --, (LMA_C._LM_ERROR_SINGULAR_MATRIX, we don't treat this as an error) + --, (LMA_C._LM_ERROR_SUM_OF_SQUARES_NOT_FINITE, we don't treat this as an error) + ] + +convertLevMarError :: CInt -> LevMarError +convertLevMarError err = fromMaybe (error "Unknown levmar error") $ + lookup err levmarCErrorToLevMarError + + +-- The End --------------------------------------------------------------------- diff --git a/LevMar/Intermediate/Fitting.hs b/LevMar/Intermediate/Fitting.hs new file mode 100644 --- /dev/null +++ b/LevMar/Intermediate/Fitting.hs @@ -0,0 +1,82 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : LevMar.Intermediate.Fitting +-- Copyright : (c) 2009 Roel van Dijk & Bas van Dijk +-- License : BSD-style (see the file LICENSE) +-- +-- Maintainer : vandijk.roel@gmail.com, v.dijk.bas@gmail.com +-- Stability : Experimental +-- +-- This module provides the Levenberg-Marquardt algorithm specialised +-- for curve-fitting. +-- +-- For additional documentation see the documentation of the levmar C +-- library which this library is based on: +-- +-- +-------------------------------------------------------------------------------- + +module LevMar.Intermediate.Fitting + ( -- * Model & Jacobian. + Model + , Jacobian + + -- * Levenberg-Marquardt algorithm. + , LMA_I.LevMarable + , levmar + + , LMA_I.LinearConstraints + + -- * Minimization options. + , LMA_I.Options(..) + , LMA_I.defaultOpts + + -- * Output + , LMA_I.Info(..) + , LMA_I.StopReason(..) + , LMA_I.CovarMatrix + + , LMA_I.LevMarError(..) + ) where + + +import qualified LevMar.Intermediate as LMA_I + + +-------------------------------------------------------------------------------- +-- Model & Jacobian. +-------------------------------------------------------------------------------- + +type Model r a = [r] -> a -> r + +-- | See: +type Jacobian r a = [r] -> a -> [r] + + +-------------------------------------------------------------------------------- +-- Levenberg-Marquardt algorithm. +-------------------------------------------------------------------------------- + +-- | The Levenberg-Marquardt algorithm specialised for curve-fitting. +levmar :: LMA_I.LevMarable r + => Model r a -- ^ Model + -> Maybe (Jacobian r a) -- ^ Optional jacobian + -> [r] -- ^ Initial parameters + -> [(a, r)] -- ^ Samples + -> Integer -- ^ Maximum iterations + -> LMA_I.Options r -- ^ Minimization options + -> Maybe [r] -- ^ Optional lower bounds + -> Maybe [r] -- ^ Optional upper bounds + -> Maybe (LMA_I.LinearConstraints r) -- ^ Optional linear constraints + -> Maybe [r] -- ^ Optional weights + -> Either LMA_I.LevMarError ([r], LMA_I.Info r, LMA_I.CovarMatrix r) +levmar model mJac ps samples = + LMA_I.levmar (\ps' -> map (model ps') xs) + (fmap (\jac -> \ps' -> map (jac ps') xs) mJac) + ps + ys + where + (xs, ys) = unzip samples + + +-- The End --------------------------------------------------------------------- diff --git a/NFunction.hs b/NFunction.hs new file mode 100644 --- /dev/null +++ b/NFunction.hs @@ -0,0 +1,54 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module NFunction + ( NFunction + , ($*) + , ComposeN + , compose + ) where + +import TypeLevelNat (Z(..), S(..), Nat) +import SizedList (SizedList(..)) + +-- | A @NFunction n a b@ is a function which takes @n@ arguments of +-- type @a@ and returns a @b@. +-- For example: @NFunction (S (S (S Z))) a b ~ (a -> a -> a -> b)@ +type family NFunction n a b :: * + +type instance NFunction Z a b = b +type instance NFunction (S n) a b = a -> NFunction n a b + +-- | @f $* xs@ applies the /n/-arity function @f@ to each of the arguments in +-- the /n/-sized list @xs@. +($*) :: NFunction n a b -> SizedList n a -> b +f $* Nil = f +f $* (x ::: xs) = f x $* xs + +infixr 0 $* -- same as $ + +class Nat n => ComposeN n where + -- | Composition of NFunctions. + -- + -- Note that the @n@ and @a@ arguments are used by the type + -- checker to select the right @ComposeN@ instance. They are + -- usally given as @(witnessNat :: n)@ and @(undefined :: a)@. + compose :: forall a b c. n -> a + -> (b -> c) -> NFunction n a b -> NFunction n a c + +instance ComposeN Z where + compose Z _ = ($) + +instance ComposeN n => ComposeN (S n) where + compose (S n) (_ :: a) f g = compose n (undefined :: a) f . g + +{- +TODO: The following does not work as expected. +See: http://www.haskell.org/pipermail/haskell-cafe/2009-August/065850.html + +-- | @f .* g@ composes @f@ with the /n/-arity function @g@. +(.*) :: forall n a b c. (ComposeN n) => (b -> c) -> NFunction n a b -> NFunction n a c +(.*) = compose (witnessNat :: n) (undefined :: a) + +infixr 9 .* -- same as . +-} diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,3 @@ +import Distribution.Simple + +main = defaultMain diff --git a/SizedList.hs b/SizedList.hs new file mode 100644 --- /dev/null +++ b/SizedList.hs @@ -0,0 +1,76 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module SizedList + ( SizedList(..) + , toList + , fromList + , unsafeFromList + , length + , replicate + ) where + +import Prelude hiding (replicate, length) +import Data.Maybe (fromMaybe) +import TypeLevelNat (Z(..), S(..), Nat, induction, witnessNat, N(..)) + +-- | A list which is indexed with a type-level natural that denotes the size of +-- the list. +data SizedList n a where + Nil :: SizedList Z a + (:::) :: a -> SizedList n a -> SizedList (S n) a + +infixr 5 ::: -- Same precedence and associativity as (:) + +consPrecedence :: Int +consPrecedence = 5 + +instance Show a => Show (SizedList n a) where + showsPrec _ Nil = showString "Nil" + showsPrec p (x ::: xs) = showParen (p > consPrecedence) + $ showsPrec (consPrecedence + 1) x + . showString " ::: " + . showsPrec consPrecedence xs + +newtype ToList a n = ToList { unToList :: SizedList n a -> [a] } + +-- | Convert a @SizedList@ to a normal list. +toList :: forall a n. Nat n => SizedList n a -> [a] +toList = unToList $ induction (witnessNat :: n) + (ToList tl0) + (ToList . tlS . unToList) + where + tl0 :: SizedList Z a -> [a] + tl0 Nil = [] + + tlS :: forall x. Nat x => (SizedList x a -> [a]) -> SizedList (S x) a -> [a] + tlS f (x ::: xs) = x : f xs + +newtype FromList a n = FromList { unFromList :: [a] -> Maybe (SizedList n a) } + +-- | Convert a normal list to a @SizeList@. If the length of the given +-- list does not equal @n@, @Nothing@ is returned. +fromList :: forall a n. Nat n => [a] -> Maybe (SizedList n a) +fromList = unFromList $ induction (witnessNat :: n) + (FromList fl0) + (FromList . flS . unFromList) + where + fl0 [] = Just Nil + fl0 _ = Nothing + + flS _ [] = Nothing + flS k (x:xs) = fmap (x :::) $ k xs + +-- | Convert a normal list to a @SizeList@. If the length of the given +-- list does not equal @n@, an error is thrown. +unsafeFromList :: forall a n. Nat n => [a] -> SizedList n a +unsafeFromList = fromMaybe (error "unsafeFromList xs: xs does not have the right length ") . + fromList + +replicate :: N n -> a -> SizedList n a +replicate Zero _ = Nil +replicate (Succ n) x = x ::: replicate n x + +length :: SizedList n a -> N n +length Nil = Zero +length (_ ::: xs) = Succ $ length xs diff --git a/TypeLevelNat.hs b/TypeLevelNat.hs new file mode 100644 --- /dev/null +++ b/TypeLevelNat.hs @@ -0,0 +1,99 @@ +-- Thanks to Ryan Ingram who wrote most of this module. +-- See: http://www.haskell.org/pipermail/haskell-cafe/2009-August/065674.html + +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} + +module TypeLevelNat + ( Z(..) + , S(..) + , Nat + , caseNat + , induction + , witnessNat + + , N(..) + ) where + + +-- | Type-level natural denoting zero +data Z = Z deriving Show + +-- | Type-level natural denoting the /S/uccessor of another type-level natural. +newtype S n = S n deriving Show + +-- | Class of all type-level naturals. +class Nat n where + -- | Case analysis on natural numbers. + caseNat :: forall r. + n -- ^ The natural number to case analyse. + -> (n ~ Z => r) -- ^ The result @r@ when @n@ equals zero. + -> (forall p. (n ~ S p, Nat p) => p -> r) -- ^ Function to apply to the predecessor + -- of @n@ to yield the result @r@. + -> r + +instance Nat Z where + caseNat _ z _ = z + +instance Nat p => Nat (S p) where + caseNat (S p) _ s = s p + +-- | The axiom of induction on natural numbers. +-- See: +induction :: forall p n. Nat n + => n + -> p Z + -> (forall m. Nat m => p m -> p (S m)) + -> p n +induction n z s = caseNat n isZ isS + where + isZ :: n ~ Z => p n + isZ = z + + isS :: forall m. (n ~ S m, Nat m) => m -> p n + isS m = s (induction m z s) + +newtype Witness x = Witness { unWitness :: x } + +-- | The value of @witnessNat :: n@ is the natural number of type @n@. +-- For example: +-- +-- @ +-- *TypeLevelNat> witnessNat :: S (S (S Z)) +-- S (S (S Z)) +-- @ +witnessNat :: forall n. Nat n => n +witnessNat = theWitness + where + theWitness = unWitness $ induction (undefined `asTypeOf` theWitness) + (Witness Z) + (Witness . S . unWitness) + +-- | A value-level natural indexed with an equivalent type-level natural. +data N n where + Zero :: N Z + Succ :: N n -> N (S n) + +{- +Template Haskell code to construct a type synonym for an arbitrary +type level natural number. + +Instead of + +> type N6 = S (S (S (S (S (S Z))))) + +you can write + +> $(mkNat "N6" 6) +-} + +-- import Language.Haskell.TH.Syntax + +-- mkNat :: String -> Int -> Q [Dec] +-- mkNat syn = runQ . return . (:[]) . TySynD (mkName syn) [] . go +-- where go 0 = ConT $ mkName "Z" +-- go n = AppT (ConT $ mkName "S") $ go (n - 1) + diff --git a/levmar.cabal b/levmar.cabal new file mode 100644 --- /dev/null +++ b/levmar.cabal @@ -0,0 +1,79 @@ +name: levmar +version: 0.1 +cabal-version: >= 1.6 +build-type: Simple +stability: experimental +author: Roel van Dijk & Bas van Dijk +maintainer: vandijk.roel@gmail.com, v.dijk.bas@gmail.com +copyright: (c) 2009 Roel van Dijk & Bas van Dijk +license: BSD3 +license-file: LICENSE +category: numerical +synopsis: An implementation of the Levenberg-Marquardt algorithm +description: The Levenberg-Marquardt algorithm is an iterative + technique that finds a local minimum of a function that + is expressed as the sum of squares of nonlinear + functions. It has become a standard technique for + nonlinear least-squares problems and can be thought of + as a combination of steepest descent and the + Gauss-Newton method. When the current solution is far + from the correct one, the algorithm behaves like a + steepest descent method: slow, but guaranteed to + converge. When the current solution is close to the + correct solution, it becomes a Gauss-Newton method. + . + Optional box- and linear constraints can be given. Both + single and double precision floating point types are + supported. + . + The actual algorithm is implemented in a C library + which is bundled with bindings-levmar which this + package depends on. See: + . + . + This library consists of two layers: + . + * LevMar.Intermediate: A medium-level layer that wraps + the low-level functions from bindings-levmar to + provide a more Haskell friendly interface. + . + * LevMar: A high-level layer that uses type-level + programming to add extra type safety. + . + Each layer also has special data-fitting variants: + . + * LevMar.Intermediate.Fitting + . + * LevMar.Fitting + . + All modules are self-contained; i.e. each module + re-exports all the things you need to work with it. + . + For an example how to use this library see Demo.hs + which is included in this package. Demo.hs is a Haskell + translation of lmdemo.c from the C levmar library. + . + A note regarding the license: + . + This library depends on bindings-levmar which is + bundled together with a C library which falls under the + GPL. Please be aware of this when distributing programs + linked with this library. For details see the + description and license of bindings-levmar. +extra-source-files: Demo.hs + +source-repository head + Type: darcs + Location: http://code.haskell.org/levmar + +library + build-depends: base >= 3 && < 4.2 + , bindings-levmar < 0.2 + exposed-modules: LevMar + , LevMar.Fitting + , LevMar.Intermediate + , LevMar.Intermediate.Fitting + , TypeLevelNat + , SizedList + , NFunction + ghc-options: -Wall -O2