module Numeric.Lbfgsb
( minimize
, minimizeV
) where
import Control.Arrow hiding (loop)
import Control.Monad
import Data.Char
import qualified Data.Vector.Generic as V
import qualified Data.Vector.Storable as SV
import Foreign hiding (unsafePerformIO)
import System.IO.Unsafe (unsafePerformIO)
readTask :: Int -> Ptr Word8 -> IO String
readTask n a = do
s <- peekArray n a
return $ map (chr . fromIntegral) s
expandConstraints :: [(Maybe Double, Maybe Double)]
-> ([Double], [Double], [Int])
expandConstraints cs = unzip3 . map conv $ cs
where
conv (Nothing, Nothing) = (47, 47, 0)
conv (Just x, Nothing) = (x, 47, 1)
conv (Just x, Just y) = (x, y, 2)
conv (Nothing, Just y) = (47, y, 3)
vectorize :: ([Double] -> (Double, [Double])) -> SV.Vector Double -> (Double, SV.Vector Double)
vectorize fg = second V.fromList . fg . V.toList
unvectorize :: (SV.Vector Double -> (Double, SV.Vector Double)) -> [Double] -> (Double, [Double])
unvectorize fg = second V.toList . fg . V.fromList
minimizeV :: Int
-> Double
-> Double
-> SV.Vector Double
-> [(Maybe Double, Maybe Double)]
-> (SV.Vector Double -> (Double, SV.Vector Double))
-> SV.Vector Double
minimizeV m factr pgtol x bounds fg = unsafePerformIO $ minimizeIO (1) m factr pgtol x bounds fg
minimize :: Int
-> Double
-> Double
-> [Double]
-> [(Maybe Double, Maybe Double)]
-> ([Double] -> (Double, [Double]))
-> [Double]
minimize m factr pgtol x bounds fg = V.toList $ minimizeV m factr pgtol (V.fromList x) bounds (vectorize fg)
minimizeIO :: Int
-> Int
-> Double
-> Double
-> SV.Vector Double
-> [(Maybe Double, Maybe Double)]
-> (SV.Vector Double -> (Double, SV.Vector Double))
-> IO (SV.Vector Double)
minimizeIO verbosity m factr pgtol x bounds fg =
with n $ \n' ->
with m $ \m' ->
withArray (V.toList x) $ \x' ->
withArray ls $ \l ->
withArray us $ \u ->
withArray nbds $ \nbd ->
alloca $ \f ->
allocaArray n $ \g ->
with factr $ \factr' ->
with pgtol $ \pgtol' ->
allocaArray ((2*m+5)*n + 11*m*m + 8*m) $ \wa ->
allocaArray (3*n + 470) $ \iwa ->
withArray sTART $ \task ->
with verbosity $ \iprint ->
allocaArray 60 $ \csave ->
allocaArray 4 $ \lsave ->
allocaArray 44 $ \isave ->
allocaArray 29 $ \dsave -> do
res <- loop n' m' x' l u nbd f g factr' pgtol' wa iwa task iprint csave lsave isave dsave
return res
where
loop n' m' x' l u nbd f g factr' pgtol' wa iwa task iprint csave lsave isave dsave = loop'
where
loop' = do
setulb_ n' m' x' l u nbd f g factr' pgtol' wa iwa task iprint csave lsave isave dsave 60 60
when (verbosity > 1) $ print =<< readTask 50 task
readTask 5 task >>= \t -> case t of
'F':'G':_ -> updateFg >> loop'
"NEW_X" -> (when (verbosity > 1) $ print =<< readDoubles n x') >>
updateFg >> loop'
_ -> readDoubles n x'
updateFg = do
xs <- readDoubles n x'
let (fnew, gnew) = fg xs
poke f fnew
pokeArray g (V.toList gnew)
n = V.length x
readDoubles :: Int -> Ptr Double -> IO (SV.Vector Double)
readDoubles n' ds = (return . V.fromList) =<< peekArray n' ds
(ls, us, nbds) = expandConstraints . take n $ bounds ++ repeat (Nothing, Nothing)
sTART = map (fromIntegral . ord) . take 60 $ "START" ++ repeat ' '
foreign import ccall setulb_
:: Ptr Int
-> Ptr Int
-> Ptr Double
-> Ptr Double
-> Ptr Double
-> Ptr Int
-> Ptr Double
-> Ptr Double
-> Ptr Double
-> Ptr Double
-> Ptr Double
-> Ptr Int
-> Ptr Word8
-> Ptr Int
-> Ptr Word8
-> Ptr Bool
-> Ptr Int
-> Ptr Double
-> Int
-> Int
-> IO ()