diff --git a/Demo.hs b/Demo.hs --- a/Demo.hs +++ b/Demo.hs @@ -1,994 +1,1451 @@ -- 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 +module Main where + +import LevMar ( levmar + + , Model + , Jacobian + + , Options(..), defaultOpts + + , LinearConstraints, noLinearConstraints + + , LevMarError + + , Info(..), CovarMatrix + + , S, Z + , SizedList(..) + ) + + +import qualified LevMar.AD as AD +import qualified LevMar.Fitting as Fitting +import qualified LevMar.Fitting.AD as Fitting.AD + +import qualified SizedList as SL (replicate) + + +-------------------------------------------------------------------------------- + +type Result n = Either LevMarError + ( SizedList n Double + , Info Double + , CovarMatrix n Double + ) + +printInteresting :: Result n -> IO () +printInteresting (Left err) = putStrLn ("Error: " ++ show err) +printInteresting (Right (ps, inf, covar)) = + do putStrLn ("infStopReason = " ++ show (infStopReason inf)) + putStrLn ("infNorm2E = " ++ show (infNorm2E inf)) + putStrLn ("infNumIter = " ++ show (infNumIter inf)) + putStrLn ("ps = " ++ show ps) + +sqr :: Num a => a -> a +sqr x = x*x + +-------------------------------------------------------------------------------- +-- 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 +type N6 = S N5 + +-------------------------------------------------------------------------------- +-- Default options: + +opts :: Options Double +opts = defaultOpts { optStopNormInfJacTe = 1e-15 + , optStopNorm2Dp = 1e-15 + , optStopNorm2E = 1e-20 + } + +-------------------------------------------------------------------------------- +-- Rosenbrock function, +-- global minimum at (1, 1) + +ros :: Floating r => Model N2 N2 r +ros p0 p1 = SL.replicate (sqr (1.0 - p0) + ros_d*sqr m) + where + m = p1 - sqr p0 + +ros_jac :: Floating r => Jacobian N2 N2 r +ros_jac p0 p1 = SL.replicate ( -2 + 2*p0 - 4*ros_d*m*p0 + ::: 2*ros_d*m + ::: Nil + ) + where + m = p1 - sqr p0 + +ros_d :: Floating r => r +ros_d = 105.0 + +ros_params :: Floating r => SizedList N2 r +ros_params = -1.2 ::: 1.0 ::: Nil + +ros_samples :: Floating r => SizedList N2 r +ros_samples = SL.replicate 0.0 + +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +-- !! TODO: Find out why these return with: infStopReason = MaxIterations !! +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +run_ros :: IO () +run_ros = printInteresting $ + levmar ros + Nothing + ros_params + ros_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +run_ros_jac :: IO () +run_ros_jac = printInteresting $ + levmar ros + (Just ros_jac) + ros_params + ros_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +run_ros_autojac :: IO () +run_ros_autojac = printInteresting $ + AD.levmar ros + ros_params + ros_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Modified Rosenbrock problem, +-- global minimum at (1, 1) + +modros :: Floating r => Model N2 N3 r +modros p0 p1 = 10*(p1 - sqr p0) + ::: 1.0 - p0 + ::: modros_lam + ::: Nil + +modros_jac :: Floating r => Jacobian N2 N3 r +modros_jac p0 _ = (-20*p0 ::: 10.0 ::: Nil) + ::: (-1.0 ::: 0.0 ::: Nil) + ::: (0.0 ::: 0.0 ::: Nil) + ::: Nil + +modros_lam :: Floating r => r +modros_lam = 1e02 + +modros_params :: Floating r => SizedList N2 r +modros_params = -1.2 ::: 1.0 ::: Nil + +modros_samples :: Floating r => SizedList N3 r +modros_samples = SL.replicate 0.0 + +run_modros :: IO () +run_modros = printInteresting $ + levmar modros + Nothing + modros_params + modros_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +run_modros_jac :: IO () +run_modros_jac = printInteresting $ + levmar modros + (Just modros_jac) + modros_params + modros_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +run_modros_autojac :: IO () +run_modros_autojac = printInteresting $ + AD.levmar modros + modros_params + modros_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Powell's function, +-- minimum at (0, 0) + +powell :: Floating r => Model N2 N2 r +powell p0 p1 = p0 + ::: 10.0*p0 / m + 2*sqr p1 + ::: Nil + where + m = p0 + 0.1 + +powell_jac :: Floating r => Jacobian N2 N2 r +powell_jac p0 p1 = (1.0 ::: 0.0 ::: Nil) + ::: (1.0 / sqr m ::: 4.0*p1 ::: Nil) + ::: Nil + where + m = p0 + 0.1 + +powell_params :: Floating r => SizedList N2 r +powell_params = -1.2 ::: 1.0 ::: Nil + +powell_samples :: Floating r => SizedList N2 r +powell_samples = SL.replicate 0.0 + +run_powell :: IO () +run_powell = printInteresting $ + levmar powell + Nothing + powell_params + powell_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +run_powell_jac :: IO () +run_powell_jac = printInteresting $ + levmar powell + (Just powell_jac) + powell_params + powell_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +-- !! TODO: Here the automatic jacobian does not seem right because !! +-- !! infNorm2E is very high compared to the manual jacobian! !! +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +run_powell_autojac :: IO () +run_powell_autojac = printInteresting $ + AD.levmar powell + powell_params + powell_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- Wood's function, +-- minimum at (1, 1, 1, 1) + +wood :: Floating r => Model N4 N6 r +wood p0 p1 p2 p3 = 10.0*(p1 - sqr p0) + ::: 1.0 - p0 + ::: sqrt 90.0*(p3 - sqr p2) + ::: 1.0 - p2 + ::: sqrt 10.0*(p1 + p3 - 2.0) + ::: (p1 - p3) / sqrt 10.0 + ::: Nil + +wood_params :: Floating r => SizedList N4 r +wood_params = -3.0 ::: -1.0 ::: -3.0 ::: -1.0 ::: Nil + +wood_samples :: Floating r => SizedList N6 r +wood_samples = SL.replicate 0.0 + +run_wood :: IO () +run_wood = printInteresting $ + levmar wood + Nothing + wood_params + wood_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +run_wood_autojac :: IO () +run_wood_autojac = printInteresting $ + AD.levmar wood + 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 :: Floating r => Fitting.SimpleModel N3 r +meyer p0 p1 p2 x = p0*exp (10.0*p1 / (ui + p2) - 13.0) + where + ui = 0.45 + 0.05*x + +meyer_jac :: Floating r => Fitting.SimpleJacobian N3 r +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_params :: Floating r => SizedList N3 r +meyer_params = 8.85 ::: 4.0 ::: 2.5 ::: Nil + +-- TODO: Unfortunately 'zip [0..] ...' won't work because (:~>) +-- doesn't have an Enum instance: +meyer_samples :: (Num a, Floating r) => [(a, r)] +meyer_samples = [ ( 0, 34.780) + , ( 1, 28.610) + , ( 2, 23.650) + , ( 3, 19.630) + , ( 4, 16.370) + , ( 5, 13.720) + , ( 6, 11.540) + , ( 7, 9.744) + , ( 8, 8.261) + , ( 9, 7.030) + , (10, 6.005) + , (11, 5.147) + , (12, 4.427) + , (13, 3.820) + , (14, 3.307) + , (15, 2.872) + ] + +run_meyer :: IO () +run_meyer = printInteresting $ + Fitting.levmar meyer + Nothing + meyer_params + meyer_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +run_meyer_jac :: IO () +run_meyer_jac = printInteresting $ + Fitting.levmar meyer + (Just meyer_jac) + meyer_params + meyer_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +-- !! TODO: Here the automatic jacobian does not seem right because !! +-- !! infNorm2E is very high compared to the manual jacobian! !! +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +run_meyer_autojac :: IO () +run_meyer_autojac = printInteresting $ + Fitting.AD.levmar meyer + meyer_params + meyer_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-------------------------------------------------------------------------------- +-- helical valley function, +-- minimum at (1.0, 0.0, 0.0) + +helval :: (Ord r, Floating r) => Model N3 N3 r +helval p0 p1 p2 = 10.0*(p2 - 10.0*theta) + ::: 10.0*sqrt tmp - 1.0 + ::: p2 + ::: Nil + where + m = atan (p1 / p0) / (2.0*pi) + + tmp = sqr p0 + sqr p1 + + theta | p0 < 0.0 = m + 0.5 + | 0.0 < p0 = m + | p1 >= 0 = 0.25 + | otherwise = -0.25 + +heval_jac :: Floating r => Jacobian N3 N3 r +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) + ::: Nil + where + tmp = sqr p0 + sqr p1 + +helval_params :: Floating r => SizedList N3 r +helval_params = -1.0 ::: 0.0 ::: 0.0 ::: Nil + +helval_samples :: Floating r => SizedList N3 r +helval_samples = SL.replicate 0.0 + +run_helval :: IO () +run_helval = printInteresting $ + levmar helval + Nothing + helval_params + helval_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +run_helval_jac :: IO () +run_helval_jac = printInteresting $ + levmar helval + (Just heval_jac) + helval_params + helval_samples + 1000 + opts + Nothing + Nothing + noLinearConstraints + Nothing + +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +-- !! TODO: This function exits with the following error: !! +-- !! : (==): No overloading for function !! +-- !! : interrupted !! +-- !! : warning: too many hs_exit()s !! +-- !! !! +-- !! Process haskell exited abnormally with code 252 !! +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +run_helval_autojac :: IO () +run_helval_autojac = printInteresting $ + AD.levmar helval + 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 :: Floating r => Model N5 N5 r +bt3 p0 p1 p2 p3 p4 = SL.replicate ( sqr t1 + + sqr t2 + + sqr t3 + + sqr t4 + ) + where + t1 = p0 - p1 + t2 = p1 + p2 - 2.0 + t3 = p3 - 1.0 + t4 = p4 - 1.0 + +bt3_jac :: Floating r => Jacobian N5 N5 r +bt3_jac p0 p1 p2 p3 p4 = SL.replicate ( 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_params :: Floating r => SizedList N5 r +bt3_params = 2.0 ::: 2.0 ::: 2.0 :::2.0 ::: 2.0 ::: Nil + +bt3_samples :: Floating r => SizedList N5 r +bt3_samples = SL.replicate 0.0 + +bt3_linear_constraints :: Floating r => LinearConstraints N3 N5 r +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 :: IO () +run_bt3 = printInteresting $ + levmar bt3 + Nothing + bt3_params + bt3_samples + 1000 + opts + Nothing + Nothing + (Just bt3_linear_constraints) + Nothing + +run_bt3_jac :: IO () +run_bt3_jac = printInteresting $ + levmar bt3 + (Just bt3_jac) + bt3_params + bt3_samples + 1000 + opts + Nothing + Nothing + (Just bt3_linear_constraints) + Nothing + +run_bt3_autojac :: IO () +run_bt3_autojac = printInteresting $ + AD.levmar bt3 + 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 :: Floating r => Model N3 N3 r +hs28 p0 p1 p2 = SL.replicate ( sqr t1 + + sqr t2 + ) + where + t1 = p0 + p1 + t2 = p1 + p2 + +hs28_jac :: Floating r => Jacobian N3 N3 r +hs28_jac p0 p1 p2 = SL.replicate ( 2.0*t1 + ::: 2.0*(t1 + t2) + ::: 2.0*t2 + ::: Nil + ) + where + t1 = p0 + p1 + t2 = p1 + p2 + +hs28_params :: Floating r => SizedList N3 r +hs28_params = -4.0 ::: 1.0 ::: 1.0 ::: Nil + +hs28_samples :: Floating r => SizedList N3 r +hs28_samples = SL.replicate 0.0 + +hs28_linear_constraints :: Floating r => LinearConstraints N1 N3 r +hs28_linear_constraints = ( ((1.0 ::: 2.0 ::: 3.0 ::: Nil) ::: Nil) + , 1.0 ::: Nil + ) + +run_hs28 :: IO () +run_hs28 = printInteresting $ + levmar hs28 + Nothing + hs28_params + hs28_samples + 1000 + opts + Nothing + Nothing + (Just hs28_linear_constraints) + Nothing + +run_hs28_jac :: IO () +run_hs28_jac = printInteresting $ + levmar hs28 + (Just hs28_jac) + hs28_params + hs28_samples + 1000 + opts + Nothing + Nothing + (Just hs28_linear_constraints) + Nothing + +run_hs28_autojac :: IO () +run_hs28_autojac = printInteresting $ + AD.levmar hs28 + 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 :: Floating r => Model N5 N5 r +hs48 p0 p1 p2 p3 p4 = SL.replicate ( sqr t1 + + sqr t2 + + sqr t3 + ) + where + t1 = p0 - 1.0 + t2 = p1 - p2 + t3 = p3 - p4 + +hs48_jac :: Floating r => Jacobian N5 N5 r +hs48_jac p0 p1 p2 p3 p4 = SL.replicate ( 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_params :: Floating r => SizedList N5 r +hs48_params = 3.0 ::: 5.0 ::: -3.0 ::: 2.0 ::: -2.0 ::: Nil + +hs48_samples :: Floating r => SizedList N5 r +hs48_samples = SL.replicate 0.0 + +hs48_linear_constraints :: Floating r => LinearConstraints N2 N5 r +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 :: IO () +run_hs48 = printInteresting $ + levmar hs48 + Nothing + hs48_params + hs48_samples + 1000 + opts + Nothing + Nothing + (Just hs48_linear_constraints) + Nothing + +run_hs48_jac :: IO () +run_hs48_jac = printInteresting $ + levmar hs48 + (Just hs48_jac) + hs48_params + hs48_samples + 1000 + opts + Nothing + Nothing + (Just hs48_linear_constraints) + Nothing + +run_hs48_autojac :: IO () +run_hs48_autojac = printInteresting $ + AD.levmar hs48 + 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 :: Floating r => Model N5 N5 r +hs51 p0 p1 p2 p3 p4 = SL.replicate ( sqr t1 + + sqr t2 + + sqr t3 + + sqr t4 + ) + where + t1 = p0 - p1 + t2 = p1 + p2 - 2.0 + t3 = p3 - 1.0 + t4 = p4 - 1.0 + +hs51_jac :: Floating r => Jacobian N5 N5 r +hs51_jac p0 p1 p2 p3 p4 = SL.replicate ( 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_params :: Floating r => SizedList N5 r +hs51_params = 2.5 ::: 0.5 ::: 2.0 ::: -1.0 ::: 0.5 ::: Nil + +hs51_samples :: Floating r => SizedList N5 r +hs51_samples = SL.replicate 0.0 + +hs51_linear_constraints :: Floating r => LinearConstraints N3 N5 r +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 :: IO () +run_hs51 = printInteresting $ + levmar hs51 + Nothing + hs51_params + hs51_samples + 1000 + opts + Nothing + Nothing + (Just hs51_linear_constraints) + Nothing + +run_hs51_jac :: IO () +run_hs51_jac = printInteresting $ + levmar hs51 + (Just hs51_jac) + hs51_params + hs51_samples + 1000 + opts + Nothing + Nothing + (Just hs51_linear_constraints) + Nothing + +run_hs51_autojac :: IO () +run_hs51_autojac = printInteresting $ + AD.levmar hs51 + 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 :: Floating r => Model N2 N2 r +hs01 p0 p1 = 10.0*(p1 - sqr p0) + ::: 1.0 - p0 + ::: Nil + +hs01_jac :: Floating r => Jacobian N2 N2 r +hs01_jac p0 _ = (-20.0*p0 ::: 10.0 ::: Nil) + ::: (-1.0 ::: 0.0 ::: Nil) + ::: Nil + +hs01_params :: Floating r => SizedList N2 r +hs01_params = -2.0 ::: 1.0 ::: Nil + +hs01_samples :: Floating r => SizedList N2 r +hs01_samples = SL.replicate 0.0 + +hs01_lb, hs01_ub :: Floating r => SizedList N2 r +hs01_lb = -_DBL_MAX ::: -1.5 ::: Nil +hs01_ub = _DBL_MAX ::: _DBL_MAX ::: Nil + +_DBL_MAX :: Floating r => r +_DBL_MAX = 1e+37 -- TODO: Get this directly from . + +run_hs01 :: IO () +run_hs01 = printInteresting $ + levmar hs01 + Nothing + hs01_params + hs01_samples + 1000 + opts + (Just hs01_lb) + (Just hs01_ub) + noLinearConstraints + Nothing + +run_hs01_jac :: IO () +run_hs01_jac = printInteresting $ + levmar hs01 + (Just hs01_jac) + hs01_params + hs01_samples + 1000 + opts + (Just hs01_lb) + (Just hs01_ub) + noLinearConstraints + Nothing + +run_hs01_autojac :: IO () +run_hs01_autojac = printInteresting $ + AD.levmar hs01 + 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 :: Floating r => Model N2 N2 r +hs21 p0 p1 = p0 / 10.0 + ::: p1 + ::: Nil + +hs21_jac :: Floating r => Jacobian N2 N2 r +hs21_jac _ _ = (0.1 ::: 0.0 ::: Nil) + ::: (0.0 ::: 1.0 ::: Nil) + ::: Nil + +hs21_params :: Floating r => SizedList N2 r +hs21_params = -1.0 ::: -1.0 ::: Nil + +hs21_samples :: Floating r => SizedList N2 r +hs21_samples = SL.replicate 0.0 + +hs21_lb, hs21_ub :: Floating r => SizedList N2 r +hs21_lb = 2.0 ::: -50.0 ::: Nil +hs21_ub = 50.0 ::: 50.0 ::: Nil + +run_hs21 :: IO () +run_hs21 = printInteresting $ + levmar hs21 + Nothing + hs21_params + hs21_samples + 1000 + opts + (Just hs21_lb) + (Just hs21_ub) + noLinearConstraints + Nothing + +run_hs21_jac :: IO () +run_hs21_jac = printInteresting $ + levmar hs21 + (Just hs21_jac) + hs21_params + hs21_samples + 1000 + opts + (Just hs21_lb) + (Just hs21_ub) + noLinearConstraints + Nothing + +run_hs21_autojac :: IO () +run_hs21_autojac = printInteresting $ + AD.levmar hs21 + 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 :: Floating r => Model N4 N4 r +hatfldb p0 p1 p2 p3 = p0 - 1.0 + ::: p0 - sqrt p1 + ::: p1 - sqrt p2 + ::: p2 - sqrt p3 + ::: Nil + +hatfldb_jac :: Floating r => Jacobian N4 N4 r +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) + ::: Nil + +hatfldb_params :: Floating r => SizedList N4 r +hatfldb_params = 0.1 ::: 0.1 ::: 0.1 ::: 0.1 ::: Nil + +hatfldb_samples :: Floating r => SizedList N4 r +hatfldb_samples = SL.replicate 0.0 + +hatfldb_lb, hatfldb_ub :: Floating r => SizedList N4 r +hatfldb_lb = 0.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil +hatfldb_ub = _DBL_MAX ::: 0.8 ::: _DBL_MAX ::: _DBL_MAX ::: Nil + +run_hatfldb :: IO () +run_hatfldb = printInteresting $ + levmar hatfldb + Nothing + hatfldb_params + hatfldb_samples + 1000 + opts + (Just hatfldb_lb) + (Just hatfldb_ub) + noLinearConstraints + Nothing + +run_hatfldb_jac :: IO () +run_hatfldb_jac = printInteresting $ + levmar hatfldb + (Just hatfldb_jac) + hatfldb_params + hatfldb_samples + 1000 + opts + (Just hatfldb_lb) + (Just hatfldb_ub) + noLinearConstraints + Nothing + +run_hatfldb_autojac :: IO () +run_hatfldb_autojac = printInteresting $ + AD.levmar hatfldb + 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 :: Floating r => Model N4 N4 r +hatfldc p0 p1 p2 p3 = p0 - 1.0 + ::: p0 - sqrt p1 + ::: p1 - sqrt p2 + ::: p3 - 1.0 + ::: Nil + +hatfldc_jac :: Floating r => Jacobian N4 N4 r +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) + ::: Nil + +hatfldc_params :: Floating r => SizedList N4 r +hatfldc_params = 0.9 ::: 0.9 ::: 0.9 ::: 0.9 ::: Nil + +hatfldc_samples :: Floating r => SizedList N4 r +hatfldc_samples = SL.replicate 0.0 + +hatfldc_lb, hatfldc_ub :: Floating r => SizedList N4 r +hatfldc_lb = 0.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil +hatfldc_ub = 10.0 ::: 10.0 ::: 10.0 ::: 10.0 ::: Nil + +run_hatfldc :: IO () +run_hatfldc = printInteresting $ + levmar hatfldc + Nothing + hatfldc_params + hatfldc_samples + 1000 + opts + (Just hatfldc_lb) + (Just hatfldc_ub) + noLinearConstraints + Nothing + +run_hatfldc_jac :: IO () +run_hatfldc_jac = printInteresting $ + levmar hatfldc + (Just hatfldc_jac) + hatfldc_params + hatfldc_samples + 1000 + opts + (Just hatfldc_lb) + (Just hatfldc_ub) + noLinearConstraints + Nothing + +run_hatfldc_autojac :: IO () +run_hatfldc_autojac = printInteresting $ + AD.levmar hatfldc + 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 :: Floating r => Model N5 N4 r +modhs52 p0 p1 p2 p3 p4 = 4.0*p0 - p1 + ::: p1 + p2 - 2.0 + ::: p3 - 1.0 + ::: p4 - 1.0 + ::: Nil + +modhs52_jac :: Floating r => Jacobian N5 N4 r +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) + ::: Nil + +modhs52_params :: Floating r => SizedList N5 r +modhs52_params = 2.0 ::: 2.0 ::: 2.0 ::: 2.0 ::: 2.0 ::: Nil + +modhs52_samples :: Floating r => SizedList N4 r +modhs52_samples = SL.replicate 0.0 + +modhs52_linear_constraints :: Floating r => LinearConstraints N3 N5 r +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 :: Floating r => SizedList N5 r +modhs52_weights = 2000.0 ::: 2000.0 ::: 2000.0 ::: 2000.0 ::: 2000.0 ::: Nil + +modhs52_lb, modhs52_ub :: Floating r => SizedList N5 r +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 :: IO () +run_modhs52 = printInteresting $ + levmar modhs52 + Nothing + modhs52_params + modhs52_samples + 1000 + opts + (Just modhs52_lb) + (Just modhs52_ub) + (Just modhs52_linear_constraints) + (Just modhs52_weights) + +run_modhs52_jac :: IO () +run_modhs52_jac = printInteresting $ + levmar modhs52 + (Just modhs52_jac) + modhs52_params + modhs52_samples + 1000 + opts + (Just modhs52_lb) + (Just modhs52_ub) + (Just modhs52_linear_constraints) + (Just modhs52_weights) + +run_modhs52_autojac :: IO () +run_modhs52_autojac = printInteresting $ + AD.levmar modhs52 + 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 :: Floating r => Model N3 N2 r +mods235 p0 p1 _ = 0.1*(p0 - 1.0) + ::: p1 - sqr p0 + ::: Nil + +mods235_jac :: Floating r => Jacobian N3 N2 r +mods235_jac p0 _ _ = (0.1 ::: 0.0 ::: 0.0 ::: Nil) + ::: (-2.0*p0 ::: 1.0 ::: 0.0 ::: Nil) + ::: Nil + +mods235_params :: Floating r => SizedList N3 r +mods235_params = -2.0 ::: 3.0 ::: 1.0 ::: Nil + +mods235_samples :: Floating r => SizedList N2 r +mods235_samples = SL.replicate 0.0 + +mods235_linear_constraints :: Floating r => LinearConstraints N2 N3 r +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 :: Floating r => SizedList N3 r +mods235_lb = -_DBL_MAX ::: 0.1 ::: 0.7 ::: Nil +mods235_ub = _DBL_MAX ::: 2.9 ::: _DBL_MAX ::: Nil + +run_mods235 :: IO () +run_mods235 = printInteresting $ + levmar mods235 + Nothing + mods235_params + mods235_samples + 1000 + opts + (Just mods235_lb) + (Just mods235_ub) + (Just mods235_linear_constraints) + Nothing + +run_mods235_jac :: IO () +run_mods235_jac = printInteresting $ + levmar mods235 + (Just mods235_jac) + mods235_params + mods235_samples + 1000 + opts + (Just mods235_lb) + (Just mods235_ub) + (Just mods235_linear_constraints) + Nothing + + +run_mods235_autojac :: IO () +run_mods235_autojac = printInteresting $ + AD.levmar mods235 + 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 :: Floating r => Model N5 N5 r +modbt7 p0 p1 _ _ _ = SL.replicate (100.0*sqr m + sqr n) + where + m = p1 - sqr p0 + n = p0 - 1.0 + +modbt7_jac :: Floating r => Jacobian N5 N5 r +modbt7_jac p0 p1 _ _ _ = SL.replicate + ( -400.0*m*p0 + 2.0*p0 - 2.0 + ::: 200.0*m + ::: 0.0 + ::: 0.0 + ::: 0.0 + ::: Nil + ) + where + m = p1 - sqr p0 + +modbt7_params :: Floating r => SizedList N5 r +modbt7_params = -2.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: Nil + +modbt7_samples :: Floating r => SizedList N5 r +modbt7_samples = SL.replicate 0.0 + +modbt7_linear_constraints :: Floating r => LinearConstraints N3 N5 r +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 :: Floating r => SizedList N5 r +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 + +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +-- !! TODO: Find out why these return with: infStopReason = MaxIterations !! +-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +run_modbt7 :: IO () +run_modbt7 = printInteresting $ + levmar modbt7 + Nothing + modbt7_params + modbt7_samples + 1000 + opts + (Just modbt7_lb) + (Just modbt7_ub) + (Just modbt7_linear_constraints) + Nothing + +run_modbt7_jac :: IO () +run_modbt7_jac = printInteresting $ + levmar modbt7 + (Just modbt7_jac) + modbt7_params + modbt7_samples + 1000 + opts + (Just modbt7_lb) + (Just modbt7_ub) + (Just modbt7_linear_constraints) + Nothing + +run_modbt7_autojac :: IO () +run_modbt7_autojac = printInteresting $ + AD.levmar modbt7 + modbt7_params + modbt7_samples + 1000 + opts + (Just modbt7_lb) + (Just modbt7_ub) + (Just modbt7_linear_constraints) + Nothing + +-------------------------------------------------------------------------------- +-- 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 :: Floating r => Model N5 N5 r +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 + ::: Nil + +r, r5, r6, r7, r8, r9, r10 :: Floating r => r +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 :: Floating r => Jacobian N5 N5 r +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 + ) + ::: Nil + +combust_params :: Floating r => SizedList N5 r +combust_params = 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: Nil + +combust_samples :: Floating r => SizedList N5 r +combust_samples = SL.replicate 0.0 + +combust_lb, combust_ub :: Floating r => SizedList N5 r +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 :: IO () +run_combust = printInteresting $ + levmar combust + Nothing + combust_params + combust_samples + 1000 + opts + (Just combust_lb) + (Just combust_ub) + noLinearConstraints + Nothing + +run_combust_jac :: IO () +run_combust_jac = printInteresting $ + levmar combust + (Just combust_jac) + combust_params + combust_samples + 1000 + opts + (Just combust_lb) + (Just combust_ub) + noLinearConstraints + Nothing + +run_combust_autojac :: IO () +run_combust_autojac = printInteresting $ + AD.levmar combust + combust_params + combust_samples + 1000 + opts + (Just combust_lb) + (Just combust_ub) + noLinearConstraints + Nothing + -- The End --------------------------------------------------------------------- diff --git a/LevMar.hs b/LevMar.hs --- a/LevMar.hs +++ b/LevMar.hs @@ -52,10 +52,18 @@ import qualified LevMar.Intermediate as LMA_I -import TypeLevelNat (Z, S, Nat) -import SizedList (SizedList(..), toList, unsafeFromList) -import NFunction (NFunction, ($*)) +import LevMar.Utils ( LinearConstraints + , noLinearConstraints + , Matrix + , CovarMatrix + , convertLinearConstraints + , convertResult + ) +import TypeLevelNat ( Z, S, Nat ) +import SizedList ( SizedList(..), toList, unsafeFromList ) +import NFunction ( NFunction, ($*) ) + import Data.Either @@ -63,25 +71,28 @@ -- Model & Jacobian. -------------------------------------------------------------------------------- -{- | A function from @n@ parameters of type @r@ to a list of @r@. +{- | A functional relation describing measurements represented as a function +from @m@ parameters to @n@ expected measurements. 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 - ] +hatfldc :: Model N4 N4 Double +hatfldc p0 p1 p2 p3 = p0 - 1.0 + ::: p0 - sqrt p1 + ::: p1 - sqrt p2 + ::: p3 - 1.0 + ::: Nil @ -} -type Model n r = NFunction n r [r] +type Model m n r = NFunction m r (SizedList n 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@ +{- | The jacobian of the 'Model' function. Expressed as a function +from @m@ parameters to a @n@/x/@m@ matrix which for each of the @n@ +expected measurement describes the @m@ partial derivatives of the +parameters. See: @@ -90,16 +101,16 @@ @ 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 - ] +hatfldc_jac :: Jacobian N4 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) + ::: Nil @ -} -type Jacobian n r = NFunction n r [SizedList n r] +type Jacobian m n r = NFunction m r (Matrix n m r) -------------------------------------------------------------------------------- @@ -107,58 +118,33 @@ -------------------------------------------------------------------------------- -- | 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 +levmar :: forall m n k r. (Nat m, Nat n, Nat k, LMA_I.LevMarable r) + => (Model m n r) -- ^ Model + -> Maybe (Jacobian m n r) -- ^ Optional jacobian + -> SizedList m r -- ^ Initial parameters + -> SizedList n 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) + -> Maybe (SizedList m r) -- ^ Optional lower bounds + -> Maybe (SizedList m r) -- ^ Optional upper bounds + -> Maybe (LinearConstraints k m r) -- ^ Optional linear constraints + -> Maybe (SizedList m r) -- ^ Optional weights + -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m 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 + (toList ys) itMax opts (fmap toList mLowBs) (fmap toList mUpBs) - (fmap convertLinC mLinC) + (fmap convertLinearConstraints 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 + convertModel f = \ps -> toList (f $* (unsafeFromList ps :: SizedList m r) :: SizedList n r) + convertJacob f = \ps -> toList (fmap toList (f $* (unsafeFromList ps :: SizedList m r) :: Matrix n m r)) -- The End --------------------------------------------------------------------- diff --git a/LevMar/AD.hs b/LevMar/AD.hs new file mode 100644 --- /dev/null +++ b/LevMar/AD.hs @@ -0,0 +1,142 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} + +-------------------------------------------------------------------------------- +-- | +-- Module : LevMar.AD +-- 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 +-- +-- A levmar variant that uses Automatic Differentiation to +-- automatically compute the Jacobian. +-- +-- For additional documentation see the documentation of the levmar C +-- library which this library is based on: +-- +-- +-------------------------------------------------------------------------------- + + +module LevMar.AD + ( -- * Model + LMA.Model + + -- * 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 as LMA +import qualified LevMar.Intermediate as LMA_I + +import LevMar.Utils ( LinearConstraints + , noLinearConstraints + , Matrix + , CovarMatrix + , convertLinearConstraints + , convertResult + ) + +import TypeLevelNat ( Z, S, Nat ) +import SizedList ( SizedList(..), toList, unsafeFromList ) +import NFunction ( NFunction, ($*) ) + +import LevMar.Utils.AD ( value, firstDeriv, constant, idDAt ) + +-- From vector-space: +import Data.Derivative ( (:~>) ) +import Data.VectorSpace ( VectorSpace, Scalar ) +import Data.Basis ( HasBasis, Basis ) + +import Data.List ( transpose ) + + +-------------------------------------------------------------------------------- +-- Levenberg-Marquardt algorithm. +-------------------------------------------------------------------------------- + +-- | The Levenberg-Marquardt algorithm that automatically computes the +-- 'Jacobian' using automatic differentiation of the model function. +-- +-- /Warning/: Don't apply 'levmar' to 'LMA.Model's that apply methods of +-- the 'Eq' and 'Ord' classes to the parameters. These methods are +-- undefined for ':~>'!!! +levmar :: forall m n k r. + ( Nat m + , Nat n + , Nat k + , HasBasis r + , Basis r ~ () + , VectorSpace (Scalar r) + , LMA_I.LevMarable r + ) + => (LMA.Model m n (r :~> r)) -- ^ Model. Note that ':~>' + -- is overloaded for all the + -- numeric classes. + -> SizedList m r -- ^ Initial parameters + -> SizedList n r -- ^ Samples + -> Integer -- ^ Maximum number of iterations + -> LMA_I.Options r -- ^ Minimization options + -> Maybe (SizedList m r) -- ^ Optional lower bounds + -> Maybe (SizedList m r) -- ^ Optional upper bounds + -> Maybe (LinearConstraints k m r) -- ^ Optional linear constraints + -> Maybe (SizedList m r) -- ^ Optional weights + -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m r) + +levmar model params ys itMax opts mLowBs mUpBs mLinC mWghts = + fmap convertResult $ LMA_I.levmar (convertModel model) + (Just $ jacobianOf model) + (toList params) + (toList ys) + itMax + opts + (fmap toList mLowBs) + (fmap toList mUpBs) + (fmap convertLinearConstraints mLinC) + (fmap toList mWghts) + where + convertModel :: LMA.Model m n (r :~> r) -> LMA_I.Model r + (convertModel mdl) ps = fmap value $ toList + (mdl $* pDs :: SizedList n (r :~> r)) + where + pDs :: SizedList m (r :~> r) + pDs = unsafeFromList $ fmap constant ps + + jacobianOf :: LMA.Model m n (r :~> r) -> LMA_I.Jacobian r + (jacobianOf mdl) ps = fmap (\fs -> zipWith (firstDeriv .) fs ps) + . transpose + . fmap (\pD -> toList (mdl $* (pD :: SizedList m (r :~> r)) :: SizedList n (r :~> r))) + $ pDs + where + pDs :: [SizedList m (r :~> r)] + pDs = [unsafeFromList $ idDAt n ps | n <- [0 .. length ps - 1]] + + +-- The End --------------------------------------------------------------------- diff --git a/LevMar/Fitting.hs b/LevMar/Fitting.hs --- a/LevMar/Fitting.hs +++ b/LevMar/Fitting.hs @@ -21,48 +21,56 @@ module LevMar.Fitting ( -- * Model & Jacobian. Model + , SimpleModel , Jacobian + , SimpleJacobian -- * Levenberg-Marquardt algorithm. - , LMA.LevMarable + , LMA_I.LevMarable , levmar - , LMA.LinearConstraints - , LMA.noLinearConstraints - , LMA.Matrix + , LinearConstraints + , noLinearConstraints + , Matrix -- * Minimization options. - , LMA.Options(..) - , LMA.defaultOpts + , LMA_I.Options(..) + , LMA_I.defaultOpts -- * Output - , LMA.Info(..) - , LMA.StopReason(..) - , LMA.CovarMatrix + , LMA_I.Info(..) + , LMA_I.StopReason(..) + , CovarMatrix - , LMA.LevMarError(..) + , LMA_I.LevMarError(..) -- *Type-level machinery , Z, S, Nat , SizedList(..) , NFunction - , ComposeN ) where -import qualified LevMar as LMA +import qualified LevMar.Intermediate.Fitting as LMA_I +import LevMar.Utils ( LinearConstraints + , noLinearConstraints + , convertLinearConstraints + , Matrix + , CovarMatrix + , convertResult + ) -import TypeLevelNat (Z, S, Nat, witnessNat) -import SizedList (SizedList) -import NFunction (NFunction, ComposeN, compose) +import TypeLevelNat ( Z, S, Nat ) +import SizedList ( SizedList(..), toList, unsafeFromList ) +import NFunction ( NFunction, ($*) ) -------------------------------------------------------------------------------- -- Model & Jacobian. -------------------------------------------------------------------------------- -{- | A function from @n@ parameters of type @r@ and an x-value of type -@a@ to a value of type @r@. +{- | A functional relation describing measurements represented as a function +from @m@ parameters and an x-value to an expected measurement. For example, the quadratic function @f(x) = a*x^2 + b*x + c@ can be written as: @@ -74,12 +82,15 @@ quad a b c x = a*x^2 + b*x + c @ -} -type Model n r a = NFunction n r (a -> r) +type Model m r a = NFunction m 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@. +-- | This type synonym expresses that usually the @a@ in @'Model' m r a@ +-- equals the type of the parameters. +type SimpleModel m r = Model m r r +{- | The jacobian of the 'Model' function. Expressed as a function from @n@ +parameters and an x-value to the @m@ partial derivatives of the parameters. + See: For example, the jacobian of the quadratic function @f(x) = a*x^2 + @@ -97,40 +108,44 @@ Notice you don't have to differentiate for @x@. -} -type Jacobian n r a = NFunction n r (a -> SizedList n r) +type Jacobian m r a = NFunction m r (a -> SizedList m r) +-- | This type synonym expresses that usually the @a@ in @'Jacobian' m r a@ +-- equals the type of the parameters. +type SimpleJacobian m r = Jacobian m r 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 +levmar :: forall m k r a. (Nat m, Nat k, LMA_I.LevMarable r) + => (Model m r a) -- ^ Model + -> Maybe (Jacobian m r a) -- ^ Optional jacobian + -> SizedList m 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 + -> LMA_I.Options r -- ^ Minimization options + -> Maybe (SizedList m r) -- ^ Optional lower bounds + -> Maybe (SizedList m r) -- ^ Optional upper bounds + -> Maybe (LinearConstraints k m r) -- ^ Optional linear constraints + -> Maybe (SizedList m r) -- ^ Optional weights + -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m 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 convertLinearConstraints mLinC) + (fmap toList mWghts) 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) + convertModel mdl = \ps -> mdl $* (unsafeFromList ps :: SizedList m r) + convertJacob jac = \ps x -> toList ((jac $* (unsafeFromList ps :: SizedList m r)) x :: SizedList m r) -- The End --------------------------------------------------------------------- diff --git a/LevMar/Fitting/AD.hs b/LevMar/Fitting/AD.hs new file mode 100644 --- /dev/null +++ b/LevMar/Fitting/AD.hs @@ -0,0 +1,139 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} + +-------------------------------------------------------------------------------- +-- | +-- Module : LevMar.Fitting.AD +-- 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 that uses Automatic Differentiation to +-- automatically compute the Jacobian. +-- +-- For additional documentation see the documentation of the levmar C +-- library which this library is based on: +-- +-- +-------------------------------------------------------------------------------- + +module LevMar.Fitting.AD + ( -- * Model. + LMA.Model + , LMA.SimpleModel + + -- * 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.Fitting as LMA +import qualified LevMar.Intermediate.Fitting as LMA_I + +import LevMar.Utils ( LinearConstraints + , noLinearConstraints + , convertLinearConstraints + , Matrix + , CovarMatrix + , convertResult + ) + +import TypeLevelNat ( Z, S, Nat ) +import SizedList ( SizedList(..), toList, unsafeFromList ) +import NFunction ( NFunction, ($*) ) + +import LevMar.Utils.AD ( value, firstDeriv, constant, idDAt ) + +-- From vector-space: +import Data.Derivative ( (:~>) ) +import Data.VectorSpace ( VectorSpace, Scalar ) +import Data.Basis ( HasBasis, Basis ) + + +-------------------------------------------------------------------------------- +-- Levenberg-Marquardt algorithm. +-------------------------------------------------------------------------------- + +-- | The Levenberg-Marquardt algorithm specialised for curve-fitting +-- that automatically computes the 'Jacobian' using automatic +-- differentiation of the model function. +-- +-- /Warning/: Don't apply 'levmar' to 'LMA_I.Model's that apply methods of +-- the 'Eq' and 'Ord' classes to the parameters. These methods are +-- undefined for ':~>'!!! +levmar :: forall m k r a. + ( Nat m + , Nat k + , HasBasis r + , Basis r ~ () + , VectorSpace (Scalar r) + , LMA_I.LevMarable r + ) + => LMA.Model m (r :~> r) a -- ^ Model. Note that + -- ':~>' is overloaded + -- for all the numeric + -- classes. + -> SizedList m r -- ^ Initial parameters + -> [(a, r)] -- ^ Samples + -> Integer -- ^ Maximum number of iterations + -> LMA_I.Options r -- ^ Minimization options + -> Maybe (SizedList m r) -- ^ Optional lower bounds + -> Maybe (SizedList m r) -- ^ Optional upper bounds + -> Maybe (LinearConstraints k m r) -- ^ Optional linear constraints + -> Maybe (SizedList m r) -- ^ Optional weights + -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m r) + +levmar model params ys itMax opts mLowBs mUpBs mLinC mWghts = + fmap convertResult $ LMA_I.levmar (convertModel model) + (Just $ jacobianOf model) + (toList params) + ys + itMax + opts + (fmap toList mLowBs) + (fmap toList mUpBs) + (fmap convertLinearConstraints mLinC) + (fmap toList mWghts) + where + convertModel :: LMA.Model m (r :~> r) a -> LMA_I.Model r a + (convertModel f) ps x = value $ (f $* pDs :: a -> r :~> r) x + where + pDs :: SizedList m (r :~> r) + pDs = unsafeFromList $ fmap constant ps + + jacobianOf :: LMA.Model m (r :~> r) a -> LMA_I.Jacobian r a + (jacobianOf f) ps x = fmap combine $ zip [0..] ps + where + combine (ix, p) = firstDeriv $ (f $* pDs :: a -> r :~> r) x p + where + pDs :: SizedList m (r :~> r) + pDs = unsafeFromList $ idDAt ix ps + + +-- The End --------------------------------------------------------------------- diff --git a/LevMar/Intermediate.hs b/LevMar/Intermediate.hs --- a/LevMar/Intermediate.hs +++ b/LevMar/Intermediate.hs @@ -42,12 +42,12 @@ ) 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 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 @@ -57,9 +57,51 @@ -- Model & Jacobian. -------------------------------------------------------------------------------- +{- | A functional relation describing measurements represented as a function +from a list of parameters to a list of expected measurements. + + * Ensure that the length of the parameters list equals the length of the + initial parameters list in 'levmar'. + + * Ensure that the length of the ouput list equals the length of the samples + list in 'levmar'. + +For example: + +@ +hatfldc :: Model Double +hatfldc [p0, p1, p2, p3] = [ p0 - 1.0 + , p0 - sqrt p1 + , p1 - sqrt p2 + , p3 - 1.0 + ] +@ +-} type Model r = [r] -> [r] --- | See: +{- | The jacobian of the 'Model' function. Expressed as a function from a list +of parameters to a list of lists which for each expected measurement describes +the partial derivatives of the parameters. + +See: + + * Ensure that the length of the parameter list equals the length of the initial + parameter list in 'levmar'. + + * Ensure that the output matrix has the dimension @n@/x/@m@ where @n@ is the + number of samples and @m@ is the number of parameters. + +For example the jacobian of the above @hatfldc@ model is: + +@ +hatfldc_jac :: Jacobian Double +hatfldc_jac _ p1 p2 _ = [ [1.0, 0.0, 0.0, 0.0] + , [1.0, -0.5 / sqrt p1, 0.0, 0.0] + , [0.0, 1.0, -0.5 / sqrt p2, 0.0] + , [0.0, 0.0, 0.0, 1.0] + ] +@ +-} type Jacobian r = [r] -> [[r]] @@ -144,7 +186,7 @@ f_blec_der f_blec_dif model mJac ps ys itMax opts mLowBs mUpBs mLinC mWeights - = unsafePerformIO $ + = unsafePerformIO . withArray (map realToFrac ps) $ \psPtr -> withArray (map realToFrac ys) $ \ysPtr -> withArray (map realToFrac $ optsToList opts) $ \optsPtr -> @@ -188,7 +230,7 @@ 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 + then return . Left $ convertLevMarError r else do result <- peekArray lenPs psPtr info <- peekArray LMA_C._LM_INFO_SZ infoPtr @@ -223,22 +265,21 @@ withLinConstraints f g = withArray (map realToFrac $ concat cMat) $ \cMatPtr -> withArray (map realToFrac rhcVec) $ \rhcVecPtr -> - f $ g cMatPtr rhcVecPtr $ fromIntegral $ length cMat + f . g cMatPtr rhcVecPtr . fromIntegral $ length cMat - withWeights f g = maybeWithArray ((fmap . fmap) realToFrac mWeights) $ \weightsPtr -> - f $ g weightsPtr + withWeights f g = maybeWithArray ((fmap . fmap) realToFrac mWeights) $ f . g 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 + 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 + 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 @@ -333,7 +374,8 @@ -------------------------------------------------------------------------------- data LevMarError - = LapackError -- ^ A call to a lapack subroutine failed in the underlying C levmar library. + = LevMarError -- ^ Generic error (not one of the others) + | 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. @@ -346,7 +388,8 @@ levmarCErrorToLevMarError :: [(CInt, LevMarError)] levmarCErrorToLevMarError = - [ (LMA_C._LM_ERROR_LAPACK_ERROR, LapackError) + [ (LMA_C._LM_ERROR, LevMarError) + , (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) diff --git a/LevMar/Intermediate/AD.hs b/LevMar/Intermediate/AD.hs new file mode 100644 --- /dev/null +++ b/LevMar/Intermediate/AD.hs @@ -0,0 +1,101 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} + +-------------------------------------------------------------------------------- +-- | +-- Module : LevMar.Intermediate.AD +-- 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 +-- +-- A levmar variant that uses Automatic Differentiation to +-- automatically compute the Jacobian. +-- +-- For additional documentation see the documentation of the levmar C +-- library which this library is based on: +-- +-- +-------------------------------------------------------------------------------- + +module LevMar.Intermediate.AD + ( -- * Model. + LMA_I.Model + + -- * 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 + +import LevMar.Utils.AD ( value, firstDeriv, constant, idDAt ) + +-- From vector-space: +import Data.Derivative ( (:~>) ) +import Data.VectorSpace ( VectorSpace, Scalar ) +import Data.Basis ( HasBasis, Basis ) + +import Data.List ( transpose ) + + +-------------------------------------------------------------------------------- +-- Levenberg-Marquardt algorithm. +-------------------------------------------------------------------------------- + +-- | The Levenberg-Marquardt algorithm that automatically computes the +-- 'Jacobian' using automatic differentiation of the model function. +-- +-- /Warning/: Don't apply 'levmar' to 'LMA_I.Model's that apply methods of +-- the 'Eq' and 'Ord' classes to the parameters. These methods are +-- undefined for ':~>'!!! +levmar :: forall r. + ( HasBasis r + , Basis r ~ () + , VectorSpace (Scalar r) + , LMA_I.LevMarable r + ) + => LMA_I.Model (r :~> r) -- ^ Model. Note that + -- ':~>' is overloaded + -- for all the numeric + -- classes. + -> [r] -- ^ Initial parameters + -> [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 = LMA_I.levmar (convertModel model) . Just $ jacobianOf model + where + convertModel :: LMA_I.Model (r :~> r) -> LMA_I.Model r + convertModel mdl = map value . mdl . map constant + + jacobianOf :: LMA_I.Model (r :~> r) -> LMA_I.Jacobian r + (jacobianOf mdl) ps = map (\fs -> zipWith (firstDeriv .) fs ps) + . transpose $ map mdl pDs + where + pDs = [idDAt n ps | n <- [0 .. length ps - 1]] + + +-- The End --------------------------------------------------------------------- diff --git a/LevMar/Intermediate/Fitting.hs b/LevMar/Intermediate/Fitting.hs --- a/LevMar/Intermediate/Fitting.hs +++ b/LevMar/Intermediate/Fitting.hs @@ -19,7 +19,9 @@ module LevMar.Intermediate.Fitting ( -- * Model & Jacobian. Model + , SimpleModel , Jacobian + , SimpleJacobian -- * Levenberg-Marquardt algorithm. , LMA_I.LevMarable @@ -47,12 +49,56 @@ -- Model & Jacobian. -------------------------------------------------------------------------------- +{- | A functional relation describing measurements represented as a function +from a list of parameters and an x-value to an expected measurement. + + * Ensure that the length of the parameters list equals the lenght of the initial + parameters list in 'levmar'. + +For example, the quadratic function @f(x) = a*x^2 + b*x + c@ can be +written as: + +@ +quad :: 'Num' r => 'Model' r r +quad [a, b, c] x = a*x^2 + b*x + c +@ +-} type Model r a = [r] -> a -> r --- | See: +-- | This type synonym expresses that usually the @a@ in @'Model' r a@ +-- equals the type of the parameters. +type SimpleModel r = Model r r + +{- | The jacobian of the 'Model' function. Expressed as a function from a list +of parameters and an x-value to the partial derivatives of the parameters. + +See: + + * Ensure that the length of the parameters list equals the lenght of the initial + parameters list in 'levmar'. + + * Ensure that the length of the output parameter derivatives list equals the + length of the input parameters list. + +For example, the jacobian of the above @quad@ model can be written as: + +@ +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 + ] +@ + +Notice you don't have to differentiate for @x@. +-} type Jacobian r a = [r] -> a -> [r] +-- | This type synonym expresses that usually the @a@ in @'Jacobian' r a@ +-- equals the type of the parameters. +type SimpleJacobian r = Jacobian r r + -------------------------------------------------------------------------------- -- Levenberg-Marquardt algorithm. -------------------------------------------------------------------------------- @@ -70,13 +116,16 @@ -> 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 +levmar model mJac params samples = + LMA_I.levmar (convertModel model) + (fmap convertJacob mJac) + params ys where (xs, ys) = unzip samples + + convertModel mdl = \ps -> map (mdl ps) xs + convertJacob jac = \ps -> map (jac ps) xs -- The End --------------------------------------------------------------------- diff --git a/LevMar/Intermediate/Fitting/AD.hs b/LevMar/Intermediate/Fitting/AD.hs new file mode 100644 --- /dev/null +++ b/LevMar/Intermediate/Fitting/AD.hs @@ -0,0 +1,100 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} + +-------------------------------------------------------------------------------- +-- | +-- Module : LevMar.Intermediate.Fitting.AD +-- 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 +-- +-- A levmar variant specialised for curve-fitting that uses Automatic +-- Differentiation to automatically compute the Jacobian. +-- +-- For additional documentation see the documentation of the levmar C +-- library which this library is based on: +-- +-- +-------------------------------------------------------------------------------- + +module LevMar.Intermediate.Fitting.AD + ( -- * Model. + LMA_I.Model + , LMA_I.SimpleModel + + -- * 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.Fitting as LMA_I + +import LevMar.Utils.AD ( value, firstDeriv, constant, idDAt ) + +-- From vector-space: +import Data.Derivative ( (:~>) ) +import Data.VectorSpace ( VectorSpace, Scalar ) +import Data.Basis ( HasBasis, Basis ) + + +-------------------------------------------------------------------------------- +-- Levenberg-Marquardt algorithm. +-------------------------------------------------------------------------------- + +-- | The Levenberg-Marquardt algorithm specialised for curve-fitting +-- that automatically computes the 'Jacobian' using automatic +-- differentiation of the model function. +-- +-- /Warning/: Don't apply 'levmar' to 'LMA_I.Model's that apply methods of +-- the 'Eq' and 'Ord' classes to the parameters. These methods are +-- undefined for ':~>'!!! +levmar :: forall r a. + ( HasBasis r + , Basis r ~ () + , VectorSpace (Scalar r) + , LMA_I.LevMarable r + ) + => LMA_I.Model (r :~> r) a -- ^ Model. Note that + -- ':~>' is overloaded + -- for all the numeric + -- classes. + -> [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 = LMA_I.levmar (convertModel model) . Just $ jacobianOf model + where + convertModel :: LMA_I.Model (r :~> r) a -> LMA_I.Model r a + convertModel mdl = \ps -> value . mdl (map constant ps) + + jacobianOf :: LMA_I.Model (r :~> r) a -> LMA_I.Jacobian r a + jacobianOf mdl = + \ps x -> map (\(ix, p) -> firstDeriv $ mdl (idDAt ix ps) x p) + $ zip [0..] ps + + +-- The End --------------------------------------------------------------------- diff --git a/LevMar/Utils.hs b/LevMar/Utils.hs new file mode 100644 --- /dev/null +++ b/LevMar/Utils.hs @@ -0,0 +1,44 @@ +module LevMar.Utils + ( LinearConstraints + , noLinearConstraints + , Matrix + , CovarMatrix + , convertLinearConstraints + , convertResult + ) where + +import qualified LevMar.Intermediate as LMA_I + +import TypeLevelNat ( Nat, Z ) +import SizedList ( SizedList, toList, unsafeFromList ) + +-- | 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 + +convertLinearConstraints :: (Nat k, Nat n) => LinearConstraints k n r -> LMA_I.LinearConstraints r +convertLinearConstraints (cMat, rhcVec) = ( map toList $ toList cMat + , toList rhcVec + ) + +convertResult :: (Nat n) + => ([r], LMA_I.Info r, LMA_I.CovarMatrix r) + -> (SizedList n r, LMA_I.Info r, CovarMatrix n r) +convertResult (psResult, info, covar) = ( unsafeFromList psResult + , info + , unsafeFromList $ map unsafeFromList covar + ) diff --git a/LevMar/Utils/AD.hs b/LevMar/Utils/AD.hs new file mode 100644 --- /dev/null +++ b/LevMar/Utils/AD.hs @@ -0,0 +1,42 @@ +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} + +module LevMar.Utils.AD where + +import Data.Derivative ( (:~>), (:>), powVal, idD, pureD, derivAtBasis ) +import Data.VectorSpace ( VectorSpace, Scalar, AdditiveGroup ) +import Data.Basis ( HasBasis, Basis ) +import Data.MemoTrie ( HasTrie ) + + +value :: a :~> b -> b +value m = powVal $ m undefined + +-- | @firstDeriv f@ returns the first derivative of @f@. +firstDeriv :: (HasBasis a, Basis a ~ (), AdditiveGroup b) + => (a :> b) -> b +firstDeriv f = powVal $ derivAtBasis f () + +-- | A constant infinitely differentiable function. +constant :: (AdditiveGroup b, HasBasis a, HasTrie (Basis a)) + => b -> a:~>b +constant = const . pureD + +-- | @idDAt n ps@ maps each parameter in @ps@ to a /constant/ +-- infinitely differentiable function (@const . pureD@), except the @n@th +-- parameter is replaced with the differentiable /identity/ function +-- (@idD@). +idDAt :: (HasBasis r, HasTrie (Basis r), VectorSpace (Scalar r)) + => Int -> [r] -> [r :~> r] +idDAt n = replace n idD . map constant + +-- | @replace i r xs@ replaces the @i@th element in @xs@ with @r@. +replace :: Int -> a -> [a] -> [a] +replace i r xs + | i < 0 = xs + | otherwise = rep i xs + where rep _ [] = [] + rep j (y:ys) + | j > 0 = y : rep (j - 1) ys + | otherwise = r : ys diff --git a/NFunction.hs b/NFunction.hs --- a/NFunction.hs +++ b/NFunction.hs @@ -8,8 +8,8 @@ , compose ) where -import TypeLevelNat (Z(..), S(..), Nat) -import SizedList (SizedList(..)) +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@. diff --git a/SizedList.hs b/SizedList.hs --- a/SizedList.hs +++ b/SizedList.hs @@ -1,27 +1,41 @@ {-# LANGUAGE GADTs #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE Rank2Types #-} module SizedList ( SizedList(..) + , foldr + , foldrN , toList + , length , fromList , unsafeFromList - , length , replicate ) where -import Prelude hiding (replicate, length) -import Data.Maybe (fromMaybe) -import TypeLevelNat (Z(..), S(..), Nat, induction, witnessNat, N(..)) +import Prelude hiding ( foldr, 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 +instance Functor (SizedList n) where + fmap _ Nil = Nil + fmap f (x ::: xs) = f x ::: fmap f xs + infixr 5 ::: -- Same precedence and associativity as (:) + +-------------------------------------------------------------------------------- + consPrecedence :: Int consPrecedence = 5 @@ -32,31 +46,46 @@ . 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) +-------------------------------------------------------------------------------- + +-- | Fold a binary operator over a @SizedList@. +foldr :: forall a b n. (a -> b -> b) -> b -> SizedList n a -> b +foldr f z = foldr_f_z where - tl0 :: SizedList Z a -> [a] - tl0 Nil = [] + foldr_f_z :: forall k. SizedList k a -> b + foldr_f_z Nil = z + foldr_f_z (x ::: xs) = f x $ foldr_f_z xs - tlS :: forall x. Nat x => (SizedList x a -> [a]) -> SizedList (S x) a -> [a] - tlS f (x ::: xs) = x : f xs +-- | Fold a binary operator yielding a value with a natural number +-- indexed type over a @SizedList@. +foldrN :: forall a b n. (forall m. a -> b m -> b (S m)) -> b Z -> SizedList n a -> b n +foldrN f z = foldrN_f_z + where + foldrN_f_z :: forall k. SizedList k a -> b k + foldrN_f_z Nil = z + foldrN_f_z (x ::: xs) = f x $ foldrN_f_z xs -newtype FromList a n = FromList { unFromList :: [a] -> Maybe (SizedList n a) } +-- | Convert a @SizedList@ to a normal list. +toList :: SizedList n a -> [a] +toList = foldr (:) [] --- | Convert a normal list to a @SizeList@. If the length of the given +-- | Returns the length of the @SizedList@. +length :: SizedList n a -> N n +length = foldrN (const Succ) Zero + + +-------------------------------------------------------------------------------- + +newtype FromList a n = FL { unFL :: [a] -> Maybe (SizedList n a) } + +-- | Convert a normal list to a @SizedList@. 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) +fromList = unFL $ induction (witnessNat :: n) (FL flZ) (FL . flS . unFL) where - fl0 [] = Just Nil - fl0 _ = Nothing + flZ [] = Just Nil + flZ _ = Nothing flS _ [] = Nothing flS k (x:xs) = fmap (x :::) $ k xs @@ -67,10 +96,14 @@ 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 +-------------------------------------------------------------------------------- + +newtype Replicate a n = R { unR :: SizedList n a} + +-- | @replicate x :: SizedList n a@ returns a @SizedList@ of @n@ @x@s. +replicate :: forall a n. Nat n => a -> SizedList n a +replicate x = unR $ induction (witnessNat :: n) (R Nil) (R . (x :::) . unR) + + +-- The End --------------------------------------------------------------------- diff --git a/TypeLevelNat.hs b/TypeLevelNat.hs --- a/TypeLevelNat.hs +++ b/TypeLevelNat.hs @@ -16,6 +16,7 @@ , witnessNat , N(..) + , nat ) where @@ -76,6 +77,9 @@ data N n where Zero :: N Z Succ :: N n -> N (S n) + +nat :: forall n. Nat n => n -> N n +nat n = induction n Zero Succ {- Template Haskell code to construct a type synonym for an arbitrary diff --git a/levmar.cabal b/levmar.cabal --- a/levmar.cabal +++ b/levmar.cabal @@ -1,14 +1,15 @@ name: levmar -version: 0.1 +version: 0.2 cabal-version: >= 1.6 build-type: Simple stability: experimental +tested-with: GHC ==6.10.4 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 +category: Numerical, Math 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 @@ -40,12 +41,26 @@ * LevMar: A high-level layer that uses type-level programming to add extra type safety. . - Each layer also has special data-fitting variants: + Each layer also has special curve-fitting variants: . * LevMar.Intermediate.Fitting . * LevMar.Fitting . + Each layer also has special variants that automatically compute + the jacobian using automatic differentiation using Conal + Elliott's vector-space library: + . + * LevMar.Intermediate.AD + . + * LevMar.Intermediate.Fitting.AD + . + * LevMar.AD + . + * LevMar.Fitting.AD + . + Note however that this feature is still very experimental! + . All modules are self-contained; i.e. each module re-exports all the things you need to work with it. . @@ -68,12 +83,20 @@ library build-depends: base >= 3 && < 4.2 - , bindings-levmar < 0.2 + , bindings-levmar == 0.1.* + , vector-space >= 0.5.7 && < 0.6 + , MemoTrie >= 0.4.5 && < 0.5 exposed-modules: LevMar + , LevMar.AD , LevMar.Fitting + , LevMar.Fitting.AD , LevMar.Intermediate + , LevMar.Intermediate.AD , LevMar.Intermediate.Fitting + , LevMar.Intermediate.Fitting.AD , TypeLevelNat , SizedList , NFunction + other-modules: LevMar.Utils + , LevMar.Utils.AD ghc-options: -Wall -O2