-- This module is a Haskell translation of lmdemo.c from the C levmar library. 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 ---------------------------------------------------------------------