ipopt-hs-0.0.0.0: Test1.hs

{- | Translation of the "number 71 from the Hock-Schittkowsky test suite"
<http://www.coin-or.org/Ipopt/documentation/node28.html> to check that the binding works.


-}
module Main where


import Ipopt
    ( wrapIpF,
      wrapIpG,
      wrapIpGradF,
      wrapIpJacG,
      wrapIpH,
      createIpoptProblem,
      addIpoptNumOption,
      addIpoptStrOption,
      ipoptSolve )
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VM
import Data.Vector.Storable ( (!) )
import Control.Monad ( forM_, forM )
import Data.IORef ( writeIORef, readIORef, newIORef )
import Text.Printf ( printf )

main = do
    let n = 4
        m = 2

        fromList = VS.unsafeThaw . VS.fromList
    xL <- fromList (replicate n 1)
    xU <- fromList (replicate n 5)

    gL <- fromList [25, 40]
    gU <- fromList [2.0e19,40]

    eval_f <- wrapIpF $ \y -> do
                x <- VS.unsafeFreeze y
                return (x!0 * x!3 * (x!0 + x!1 + x!2) + x!2)

    eval_grad_f <- wrapIpGradF $ \y -> do
        x <- VS.unsafeFreeze y
        fromList
            [ x!0 * x!3 + x!3 * (x!0 + x!1 + x!2),
              x!0 * x!3,
              x!0 * x!3 + 1,
              x!0 * (x!0 + x!1 + x!2)]

    eval_g <- wrapIpG $ \y -> do
          x <- VS.unsafeFreeze y
          fromList [x!0 * x!1 * x!2 * x!3,
                    x!0*x!0 + x!1*x!1 + x!2*x!2 + x!3*x!3]

    eval_jac_g <- wrapIpJacG
      (\ iRow jCol -> do
            print (VM.length iRow)
            VM.set (VM.slice 0 4 iRow) 0
            VM.set (VM.slice 4 4 iRow) 1
            forM_ (zip [0 .. 7] $ cycle [0 .. 3]) $ \(a,b) ->
                VM.write jCol a b
        )
       (\ x values -> do
            x <- VS.unsafeFreeze x
            values_new <- fromList [
                x!1*x!2*x!3,
                x!0*x!2*x!3,
                x!0*x!1*x!3,
                x!0*x!1*x!2,

                2*x!0,
                2*x!1,
                2*x!2,
                2*x!3 ]
            VM.copy values values_new
            )


    eval_h <- wrapIpH
     ( \ iRow jCol -> do
        i <- newIORef 0
        forM [0 .. 3] $ \ row ->
            forM [ 0 .. row ] $ \col -> do
                ii <- readIORef i
                VM.write iRow ii row
                VM.write jCol ii col
                writeIORef i (ii+1)
       )

     ( \ obj_factor lambda x values -> do

        lambda <- VS.unsafeFreeze lambda
        x <- VS.unsafeFreeze x

        VM.copy values =<< fromList
            [ obj_factor * (2*x!3) + lambda!1 * 2,
              obj_factor * (x!3) + lambda!0 * (x!2 * x!3),
              lambda!1 * 2,
              obj_factor * (x!3) + lambda!0 * (x!1 * x!3),
              lambda!0 * (x!0 * x!3),
              lambda!1 * 2,
              obj_factor * (2*x!0 + x!1 + x!2) + lambda!0 * (x!1 * x!2),
              obj_factor * (x!0) + lambda!0 * (x!0 * x!2),
              obj_factor * (x!0) + lambda!0 * (x!0 * x!1),
              lambda!1 * 2 ]
        )

    nlp <- createIpoptProblem xL xU gL gU 8 10 eval_f eval_g eval_grad_f eval_jac_g eval_h

    addIpoptNumOption nlp "tol" 1.0e-9
    addIpoptStrOption nlp "mu_strategy" "adaptive"

    x <- fromList [1,5, 5,1]

    ans <- ipoptSolve nlp x

    let x_official = VS.fromList [1, 4.74299964, 3.82114998, 1.37940829]
    x <- VS.freeze x

    let sse :: Double
        sse = realToFrac $ VS.sum $ VS.zipWith (\a b -> (a-b)^2) x x_official

    printf "\n\n||x - x_official||_2 = %f\n" sse