numeric-optimization-ad (empty) → 0.1.0.0

raw patch · 9 files changed

+594/−0 lines, 9 filesdep +HUnitdep +addep +basesetup-changed

Dependencies added: HUnit, ad, base, containers, data-default-class, hmatrix, hspec, numeric-optimization, numeric-optimization-ad, primitive, reflection, vector

Files

CHANGELOG.md +11/−0
LICENSE +30/−0
README.md +25/−0
Setup.hs +2/−0
examples/rosenbrock.hs +17/−0
numeric-optimization-ad.cabal +85/−0
src/Numeric/Optimization/AD.hs +270/−0
test/IsClose.hs +134/−0
test/Spec.hs +20/−0

+ CHANGELOG.md view

@@ -0,0 +1,11 @@+# Changelog for `numeric-optimization-ad`++All notable changes to this project will be documented in this file.++The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),+and this project adheres to the+[Haskell Package Versioning Policy](https://pvp.haskell.org/).++## Unreleased++## 0.1.0.0 - YYYY-MM-DD

+ LICENSE view

@@ -0,0 +1,30 @@+Copyright Masahiro Sakai (c) 2023++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Masahiro Sakai nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

+ README.md view

@@ -0,0 +1,25 @@+# numeric-optimization-ad++Wrapper of [numeric-optimization](https://hackage.haskell.org/package/numeric-optimization) package for using with [ad](https://hackage.haskell.org/package/ad) package++## Example Usage++```haskell+{-# LANGUAGE FlexibleContexts #-}+import Numeric.Optimization.AD++main :: IO ()+main = do+  result <- minimize LBFGS def rosenbrock Nothing [] [-3,-4]+  print (resultSuccess result)  -- True+  print (resultSolution result)  -- [0.999999999009131,0.9999999981094296]+  print (resultValue result)  -- 1.8129771632403013e-18++-- https://en.wikipedia.org/wiki/Rosenbrock_function+rosenbrock :: Floating a => [a] -> a+-- rosenbrock :: Reifies s Tape => [Reverse s Double] -> Reverse s Double+rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)++sq :: Floating a => a -> a+sq x = x ** 2+```

+ Setup.hs view

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

+ examples/rosenbrock.hs view

@@ -0,0 +1,17 @@+{-# LANGUAGE FlexibleContexts #-}+import Numeric.Optimization.AD++main :: IO ()+main = do+  result <- minimize LBFGS def rosenbrock Nothing [] [-3,-4]+  print (resultSuccess result)  -- True+  print (resultSolution result)  -- [0.999999999009131,0.9999999981094296]+  print (resultValue result)  -- 1.8129771632403013e-18++-- https://en.wikipedia.org/wiki/Rosenbrock_function+rosenbrock :: Floating a => [a] -> a+-- rosenbrock :: Reifies s Tape => [Reverse s Double] -> Reverse s Double+rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)++sq :: Floating a => a -> a+sq x = x ** 2

+ numeric-optimization-ad.cabal view

@@ -0,0 +1,85 @@+cabal-version: 1.12++-- This file has been generated from package.yaml by hpack version 0.35.1.+--+-- see: https://github.com/sol/hpack++name:           numeric-optimization-ad+version:        0.1.0.0+synopsis:       Wrapper of numeric-optimization package for using with AD package+description:    Please see the README on GitHub at <https://github.com/msakai/nonlinear-optimization-ad/tree/master/numeric-optimization-ad#readme>+category:       Math, Algorithms, Optimisation, Optimization+homepage:       https://github.com/msakai/numeric-optimization-ad#readme+bug-reports:    https://github.com/msakai/numeric-optimization-ad/issues+author:         Masahiro Sakai+maintainer:     masahiro.sakai@gmail.com+copyright:      Masahiro Sakai &lt;masahiro.sakai@gmail.com&gt;+license:        BSD3+license-file:   LICENSE+build-type:     Simple+tested-with:+    GHC == 9.4.5+  , GHC == 9.2.7+  , GHC == 9.0.2+  , GHC == 8.10.7+  , GHC == 8.8.4+  , GHC == 8.6.5+extra-source-files:+    README.md+    CHANGELOG.md++source-repository head+  type: git+  location: https://github.com/msakai/numeric-optimization-ad++library+  exposed-modules:+      Numeric.Optimization.AD+  other-modules:+      Paths_numeric_optimization_ad+  hs-source-dirs:+      src+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints+  build-depends:+      ad >=4.3.6 && <4.6+    , base >=4.12 && <5+    , data-default-class+    , hmatrix >=0.20.0.0+    , numeric-optimization >=0.1.0.0 && <0.2.0.0+    , primitive >=0.6.4.0+    , reflection >=2.1.5+    , vector >=0.12.0.2 && <0.14+  default-language: Haskell2010++executable rosenbrock-ad+  main-is: rosenbrock.hs+  other-modules:+      Paths_numeric_optimization_ad+  hs-source-dirs:+      examples+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N+  build-depends:+      base >=4.12 && <5+    , data-default-class+    , numeric-optimization >=0.1.0.0 && <0.2.0.0+    , numeric-optimization-ad+  default-language: Haskell2010++test-suite numeric-optimization-ad-test+  type: exitcode-stdio-1.0+  main-is: Spec.hs+  other-modules:+      IsClose+      Paths_numeric_optimization_ad+  hs-source-dirs:+      test+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N+  build-depends:+      HUnit >=1.6.0.0 && <1.7+    , base >=4.12 && <5+    , containers >=0.6.0.1 && <0.7+    , data-default-class+    , hspec >=2.7.1 && <3.0+    , numeric-optimization >=0.1.0.0 && <0.2.0.0+    , numeric-optimization-ad+  default-language: Haskell2010

+ src/Numeric/Optimization/AD.hs view

@@ -0,0 +1,270 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.Optimization.AD+-- Copyright   :  (c) Masahiro Sakai 2023+-- License     :  BSD-style+--+-- Maintainer  :  masahiro.sakai@gmail.com+-- Stability   :  provisional+-- Portability :  non-portable+--+-- This module is a wrapper of "Numeric.Optimization" that uses+-- [ad](https://hackage.haskell.org/package/ad)'s automatic differentiation.+--+-----------------------------------------------------------------------------+module Numeric.Optimization.AD+  (+  -- * Main function+    minimize+  , minimizeReverse+  , minimizeSparse++  -- * Problem specification+  , Constraint (..)++  -- * Algorithm selection+  , Method (..)+  , isSupportedMethod+  , Params (..)++  -- * Result+  , Result (..)+  , Statistics (..)+  , OptimizationException (..)++  -- * Utilities and Re-exports+  , Default (..)+  , AD+  , auto+  , Reverse+  , Reifies+  , Tape+  , Sparse+  ) where+++import Control.Monad.Primitive+import Data.Default.Class+import Data.Foldable (foldlM, toList)+import Data.Functor.Contravariant+import Data.Reflection (Reifies)+import Data.Traversable+import qualified Data.Vector as V+import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Generic.Mutable as VGM+import Numeric.AD (AD, auto)+import Numeric.AD.Internal.Reverse (Tape)+import Numeric.AD.Mode.Reverse (Reverse)+import qualified Numeric.AD.Mode.Reverse as Reverse+import Numeric.AD.Mode.Sparse (Sparse)+import qualified Numeric.AD.Mode.Sparse as Sparse+import qualified Numeric.LinearAlgebra as LA+import qualified Numeric.Optimization as Opt+import Numeric.Optimization hiding (minimize, IsProblem (..))++-- ------------------------------------------------------------------------++data ProblemReverse f+  = ProblemReverse+      (forall s. Reifies s Tape => f (Reverse s Double) -> Reverse s Double)+      (Maybe (V.Vector (Double, Double)))+      [Constraint]+      Int+      (f Int)++instance Traversable f => Opt.IsProblem (ProblemReverse f) where+  func (ProblemReverse f _bounds _constraints _size template) x =+    fst $ Reverse.grad' f (fromVector template x)++  bounds (ProblemReverse _f bounds _constraints _size _template) = bounds++  constraints (ProblemReverse _f _bounds constraints _size _template) = constraints++instance Traversable f => Opt.HasGrad (ProblemReverse f) where+  grad (ProblemReverse func _bounds _constraints size template) =+    toVector size . Reverse.grad func . fromVector template++  grad'M (ProblemReverse f _bounds _constraints _size template) x gvec = do+    case Reverse.grad' f (fromVector template x) of+      (y, g) -> do+        writeToMVector g gvec+        return y++instance Traversable f => Opt.Optionally (Opt.HasGrad (ProblemReverse f)) where+  optionalDict = hasOptionalDict++instance Opt.Optionally (Opt.HasHessian (ProblemReverse f)) where+  optionalDict = Nothing++-- | Minimization of scalar function of one or more variables.+--+-- This is a wrapper of 'Opt.minimize' and use "Numeric.AD.Mode.Reverse" to compute gradient.+--+-- It cannot be used with methods that requires hessian (e.g. 'Newton').+--+-- Example:+--+-- > {-# LANGUAGE FlexibleContexts #-}+-- > import Numeric.Optimization.AD+-- > +-- > main :: IO ()+-- > main = do+-- >   (x, result, stat) <- minimizeReverse LBFGS def rosenbrock Nothing [] [-3,-4]+-- >   print (resultSuccess result)  -- True+-- >   print (resultSolution result)  -- [0.999999999009131,0.9999999981094296]+-- >   print (resultValue result)  -- 1.8129771632403013e-18+-- > +-- > -- https://en.wikipedia.org/wiki/Rosenbrock_function+-- > rosenbrock :: Floating a => [a] -> a+-- > -- rosenbrock :: Reifies s Tape => [Reverse s Double] -> Reverse s Double+-- > rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)+-- > +-- > sq :: Floating a => a -> a+-- > sq x = x ** 2+minimizeReverse+  :: forall f. Traversable f+  => Method  -- ^ Numerical optimization algorithm to use+  -> Params (f Double)  -- ^ Parameters for optimization algorithms. Use 'def' as a default.+  -> (forall s. Reifies s Tape => f (Reverse s Double) -> Reverse s Double)  -- ^ Function to be minimized.+  -> Maybe (f (Double, Double))  -- ^ Bounds+  -> [Constraint]  -- ^ Constraints+  -> f Double -- ^ Initial value+  -> IO (Result (f Double))+minimizeReverse method params f bounds constraints x0 = do+  let size :: Int+      template :: f Int+      (size, template) = mapAccumL (\i _ -> i `seq` (i+1, i)) 0 x0++      bounds' :: Maybe (V.Vector (Double, Double))+      bounds' = fmap (toVector size) bounds++      prob = ProblemReverse f bounds' constraints size template+      params' = contramap (fromVector template) params++  result <- Opt.minimize method params' prob (toVector size x0)+  return $ fmap (fromVector template) result++-- ------------------------------------------------------------------------++data ProblemSparse f+  = ProblemSparse+      (forall s. f (AD s (Sparse Double)) -> AD s (Sparse Double))+      (Maybe (V.Vector (Double, Double)))+      [Constraint]+      Int+      (f Int)++instance Traversable f => Opt.IsProblem (ProblemSparse f) where+  func (ProblemSparse f _bounds _constraints _size template) x =+    fst $ Sparse.grad' f (fromVector template x)++  bounds (ProblemSparse _f bounds _constraints _size _template) = bounds++  constraints (ProblemSparse _f _bounds constraints _size _template) = constraints++instance Traversable f => Opt.HasGrad (ProblemSparse f) where+  grad (ProblemSparse func _bounds _constraints size template) =+    toVector size . Sparse.grad func . fromVector template++  grad'M (ProblemSparse f _bounds _constraints _size template) x gvec = do+    case Sparse.grad' f (fromVector template x) of+      (y, g) -> do+        writeToMVector g gvec+        return y++instance Traversable f => Opt.HasHessian (ProblemSparse f) where+  hessian (ProblemSparse func _bounds _constraints size template) =+    toMatrix size . Sparse.hessian func . fromVector template+    where+      toMatrix n xss = (n LA.>< n) $ concat $ map toList $ toList xss++instance Traversable f => Opt.Optionally (Opt.HasGrad (ProblemSparse f)) where+  optionalDict = hasOptionalDict++instance Traversable f => Opt.Optionally (Opt.HasHessian (ProblemSparse f)) where+  optionalDict = hasOptionalDict++-- | Minimization of scalar function of one or more variables.+--+-- This is a wrapper of 'Opt.minimize' and use "Numeric.AD.Mode.Sparse" to compute gradient+-- and hessian.+--+-- Unlike 'minimizeReverse', it can be used with methods that requires hessian (e.g. 'Newton').+--+-- Example:+--+-- > {-# LANGUAGE FlexibleContexts #-}+-- > import Numeric.Optimization.AD+-- >+-- > main :: IO ()+-- > main = do+-- >   (x, result, stat) <- minimizeSparse Newton def rosenbrock Nothing [] [-3,-4]+-- >   print (resultSuccess result)  -- True+-- >   print (resultSolution result)  -- [0.9999999999999999,0.9999999999999998]+-- >   print (resultValue result)  -- 1.232595164407831e-32+-- >+-- > -- https://en.wikipedia.org/wiki/Rosenbrock_function+-- > rosenbrock :: Floating a => [a] -> a+-- > -- rosenbrock :: [AD s (Sparse Double)] -> AD s (Sparse Double)+-- > rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)+-- >+-- > sq :: Floating a => a -> a+-- > sq x = x ** 2+minimizeSparse+  :: forall f. Traversable f+  => Method  -- ^ Numerical optimization algorithm to use+  -> Params (f Double)  -- ^ Parameters for optimization algorithms. Use 'def' as a default.+  -> (forall s. f (AD s (Sparse Double)) -> AD s (Sparse Double))  -- ^ Function to be minimized.+  -> Maybe (f (Double, Double))  -- ^ Bounds+  -> [Constraint]  -- ^ Constraints+  -> f Double -- ^ Initial value+  -> IO (Result (f Double))+minimizeSparse method params f bounds constraints x0 = do+  let size :: Int+      template :: f Int+      (size, template) = mapAccumL (\i _ -> i `seq` (i+1, i)) 0 x0++      bounds' :: Maybe (V.Vector (Double, Double))+      bounds' = fmap (toVector size) bounds++      prob = ProblemSparse f bounds' constraints size template+      params' = contramap (fromVector template) params++  result <- Opt.minimize method params' prob (toVector size x0)+  return $ fmap (fromVector template) result++-- ------------------------------------------------------------------------++-- | Synonym of 'minimizeReverse'+minimize+  :: forall f. Traversable f+  => Method  -- ^ Numerical optimization algorithm to use+  -> Params (f Double)  -- ^ Parameters for optimization algorithms. Use 'def' as a default.+  -> (forall s. Reifies s Tape => f (Reverse s Double) -> Reverse s Double)  -- ^ Function to be minimized.+  -> Maybe (f (Double, Double))  -- ^ Bounds+  -> [Constraint]  -- ^ Constraints+  -> f Double -- ^ Initial value+  -> IO (Result (f Double))+minimize = minimizeReverse++-- ------------------------------------------------------------------------++fromVector :: (Functor f, VG.Vector v a) => f Int -> v a -> f a+fromVector template x = fmap (x VG.!) template++toVector :: (Traversable f, VG.Vector v a) => Int -> f a -> v a+toVector size x = VG.create $ do+  vec <- VGM.new size+  writeToMVector x vec+  return vec++writeToMVector :: (PrimMonad m, VGM.MVector mv a, Traversable f) => f a -> mv (PrimState m) a -> m ()+writeToMVector x vec = do+  _ <- foldlM (\i v -> VGM.write vec i v >> return (i+1)) 0 x+  return ()++-- ------------------------------------------------------------------------

+ test/IsClose.hs view

@@ -0,0 +1,134 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+module IsClose+  (+  -- Tolerance type+    Tol (..)++  -- AllClose class+  , AllClose (..)+  , allCloseRawUnit+  , allCloseRawRealFrac+  , allCloseRawRealFloat++  -- * Re-exports+  , Default (..)++  -- * HUnit+  , assertAllClose+  ) where++import Data.Default.Class+import Data.List.NonEmpty (NonEmpty (..))+import Data.Map (Map)+import qualified Data.Map as Map+import Data.Monoid+import Data.Semigroup+import GHC.Stack (HasCallStack)+import Test.HUnit+import Text.Printf++-- ------------------------------------------------------------------------++-- | Tolerance+--+-- Values @a@ and @b@ are considered /close/ if @abs (a - b) <= atol + rtol * abs b@.+data Tol a+  = Tol+  { rtol :: a -- ^ The relative tolerance parameter (default: @1e-05@)+  , atol :: a -- ^ The absolute tolerance parameter (default: @1e-08@)+  , equalNan :: Bool -- ^ Whether to compare NaN’s as equal (default: @False@)+  } deriving (Show)++instance RealFrac a => Default (Tol a) where+  def = Tol+    { rtol = 1e-05+    , atol = 1e-08+    , equalNan = False+    }++-- ------------------------------------------------------------------------++class Real r => AllClose r a where+  -- | Returns number of mismatches, number of elements, maximal absolute difference, and maximal relative difference.+  -- Returns @'Ap' 'Nothing'@ if given values are incomparable.+  allCloseRaw :: Tol r -> a -> a -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)++  -- | Returns 'True' if the two arrays are equal within the given tolerance; 'False' otherwise.+  allClose :: Tol r -> a -> a -> Bool+  allClose tol x y =+    case getAp (allCloseRaw tol x y) of+      Nothing -> False+      Just (Sum numMismatched, _, _, _) -> numMismatched == 0++allCloseRawRealFrac :: RealFrac r => Tol r -> r -> r -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)+allCloseRawRealFrac t a b = Ap $ Just $+  ( Sum $ if abs (a - b) <= atol t + rtol t * abs b then 0 else 1+  , Sum 1+  , Max (abs (a - b))+  , Max (abs (a - b) / abs b)+  )++allCloseRawRealFloat :: RealFloat r => Tol r -> r -> r -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)+allCloseRawRealFloat t a b+  | isNaN a /= isNaN b = Ap Nothing+  | otherwise = Ap $ Just $+      ( Sum $ if (equalNan t && isNaN a && isNaN b) || a == b || abs (a - b) <= atol t + rtol t * abs b then 0 else 1+      , Sum 1+      , Max (abs (a - b))+      , Max (abs (a - b) / abs b)+      )++allCloseRawUnit :: Num r => Ap Maybe (Sum Int, Sum Int, Max r, Max r)+allCloseRawUnit = Ap (Just (Sum 0, Sum 0, Max 0, Max 0))++instance AllClose Rational Rational where+  allCloseRaw = allCloseRawRealFrac++instance AllClose Double Double where+  allCloseRaw = allCloseRawRealFloat++instance (AllClose r a) => AllClose r (Maybe a) where+  allCloseRaw tol (Just a) (Just b) = allCloseRaw tol a b+  allCloseRaw _ Nothing Nothing = allCloseRawUnit+  allCloseRaw _ _ _ = Ap Nothing++instance (AllClose r v) => AllClose r [v] where+  allCloseRaw tol xs ys+    | length xs == length ys = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip xs ys])+    | otherwise = Ap Nothing++instance (Ord k, AllClose r v) => AllClose r (Map k v) where+  allCloseRaw tol m1 m2+    | Map.keys m1 == Map.keys m2 = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip (Map.elems m1) (Map.elems m2)])+    | otherwise = Ap Nothing++-- ------------------------------------------------------------------------++-- | Assert that two objects are equal up to desired tolerance.+assertAllClose+  :: (HasCallStack, AllClose r a, Show r, Show a)+  => Tol r+  -> a -- ^ actual+  -> a -- ^ desired+  -> Assertion+assertAllClose tol a b =+  case getAp (allCloseRaw tol a b) of+    Nothing ->+      assertString $ unlines $ header ++ ["x and y nan location mismatch:"] ++ footer+    Just (Sum numMismatch, Sum numTotal, Max absDiff, Max relDiff)+      | numMismatch == 0 -> return ()+      | otherwise ->+          assertString $ unlines $+            header +++            [ printf "Mismatched elements: %d / %d (%f%%)" numMismatch numTotal (fromIntegral numMismatch * 100 / fromIntegral numTotal :: Double)+            , " Max absolute difference: " ++ show absDiff+            , " Max relative difference: " ++ show relDiff+            ] ++ footer+   where+     header, footer :: [String]+     header = [printf "Not equal to tolerance rtol=%s, atol=%s" (show (rtol tol)) (show (atol tol)), ""]+     footer = [" x: " ++ show a, " y: " ++ show b]++-- ------------------------------------------------------------------------

+ test/Spec.hs view

@@ -0,0 +1,20 @@+import Test.Hspec++import Numeric.Optimization.AD+import IsClose+++main :: IO ()+main = hspec $ do+  describe "minimize" $ do+    context "when given rosenbrock function" $+      it "returns the global optimum" $ do+        result <- minimize LBFGS def rosenbrock Nothing [] [-3,-4]+        resultSuccess result `shouldBe` True+        assertAllClose (def :: Tol Double) (resultSolution result) [1,1]+++-- https://en.wikipedia.org/wiki/Rosenbrock_function+rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)+  where+    sq x = x ** 2