lbfgs (empty) → 0.0.1
raw patch · 6 files changed
+2071/−0 lines, 6 filesdep +arraydep +basesetup-changed
Dependencies added: array, base
Files
- LICENSE +234/−0
- Numeric/LBFGS.hs +259/−0
- Numeric/LBFGS/Raw.hsc +186/−0
- Setup.hs +2/−0
- cbits/lbfgs.c +1371/−0
- lbfgs.cabal +19/−0
+ LICENSE view
@@ -0,0 +1,234 @@+liblbfgs (MIT license):++Copyright (c) 1990 Jorge Nocedal+Copyright (c) 2007-2010 Naoaki Okazaki++Haskell module (Apache License version 2.0):++Copyright (c) 2010 Daniël de Kok++---++The MIT License++Permission is hereby granted, free of charge, to any person obtaining a+copy of this software and associated documentation files (the "Software"),+to deal in the Software without restriction, including without limitation+the rights to use, copy, modify, merge, publish, distribute, sublicense,+and/or sell copies of the Software, and to permit persons to whom the+Software is furnished to do so, subject to the following conditions:++The above copyright notice and this permission notice shall be included in+all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN+THE SOFTWARE.++---++ Apache License+ Version 2.0, January 2004+ http://www.apache.org/licenses/++ TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION++ 1. Definitions.++ "License" shall mean the terms and conditions for use, reproduction,+ and distribution as defined by Sections 1 through 9 of this document.++ "Licensor" shall mean the copyright owner or entity authorized by+ the copyright owner that is granting the License.++ "Legal Entity" shall mean the union of the acting entity and all+ other entities that control, are controlled by, or are under common+ control with that entity. For the purposes of this definition,+ "control" means (i) the power, direct or indirect, to cause the+ direction or management of such entity, whether by contract or+ otherwise, or (ii) ownership of fifty percent (50%) or more of the+ outstanding shares, or (iii) beneficial ownership of such entity.++ "You" (or "Your") shall mean an individual or Legal Entity+ exercising permissions granted by this License.++ "Source" form shall mean the preferred form for making modifications,+ including but not limited to software source code, documentation+ source, and configuration files.++ "Object" form shall mean any form resulting from mechanical+ transformation or translation of a Source form, including but+ not limited to compiled object code, generated documentation,+ and conversions to other media types.++ "Work" shall mean the work of authorship, whether in Source or+ Object form, made available under the License, as indicated by a+ copyright notice that is included in or attached to the work+ (an example is provided in the Appendix below).++ "Derivative Works" shall mean any work, whether in Source or Object+ form, that is based on (or derived from) the Work and for which the+ editorial revisions, annotations, elaborations, or other modifications+ represent, as a whole, an original work of authorship. For the purposes+ of this License, Derivative Works shall not include works that remain+ separable from, or merely link (or bind by name) to the interfaces of,+ the Work and Derivative Works thereof.++ "Contribution" shall mean any work of authorship, including+ the original version of the Work and any modifications or additions+ to that Work or Derivative Works thereof, that is intentionally+ submitted to Licensor for inclusion in the Work by the copyright owner+ or by an individual or Legal Entity authorized to submit on behalf of+ the copyright owner. For the purposes of this definition, "submitted"+ means any form of electronic, verbal, or written communication sent+ to the Licensor or its representatives, including but not limited to+ communication on electronic mailing lists, source code control systems,+ and issue tracking systems that are managed by, or on behalf of, the+ Licensor for the purpose of discussing and improving the Work, but+ excluding communication that is conspicuously marked or otherwise+ designated in writing by the copyright owner as "Not a Contribution."++ "Contributor" shall mean Licensor and any individual or Legal Entity+ on behalf of whom a Contribution has been received by Licensor and+ subsequently incorporated within the Work.++ 2. Grant of Copyright License. Subject to the terms and conditions of+ this License, each Contributor hereby grants to You a perpetual,+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable+ copyright license to reproduce, prepare Derivative Works of,+ publicly display, publicly perform, sublicense, and distribute the+ Work and such Derivative Works in Source or Object form.++ 3. Grant of Patent License. Subject to the terms and conditions of+ this License, each Contributor hereby grants to You a perpetual,+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable+ (except as stated in this section) patent license to make, have made,+ use, offer to sell, sell, import, and otherwise transfer the Work,+ where such license applies only to those patent claims licensable+ by such Contributor that are necessarily infringed by their+ Contribution(s) alone or by combination of their Contribution(s)+ with the Work to which such Contribution(s) was submitted. If You+ institute patent litigation against any entity (including a+ cross-claim or counterclaim in a lawsuit) alleging that the Work+ or a Contribution incorporated within the Work constitutes direct+ or contributory patent infringement, then any patent licenses+ granted to You under this License for that Work shall terminate+ as of the date such litigation is filed.++ 4. Redistribution. You may reproduce and distribute copies of the+ Work or Derivative Works thereof in any medium, with or without+ modifications, and in Source or Object form, provided that You+ meet the following conditions:++ (a) You must give any other recipients of the Work or+ Derivative Works a copy of this License; and++ (b) You must cause any modified files to carry prominent notices+ stating that You changed the files; and++ (c) You must retain, in the Source form of any Derivative Works+ that You distribute, all copyright, patent, trademark, and+ attribution notices from the Source form of the Work,+ excluding those notices that do not pertain to any part of+ the Derivative Works; and++ (d) If the Work includes a "NOTICE" text file as part of its+ distribution, then any Derivative Works that You distribute must+ include a readable copy of the attribution notices contained+ within such NOTICE file, excluding those notices that do not+ pertain to any part of the Derivative Works, in at least one+ of the following places: within a NOTICE text file distributed+ as part of the Derivative Works; within the Source form or+ documentation, if provided along with the Derivative Works; or,+ within a display generated by the Derivative Works, if and+ wherever such third-party notices normally appear. The contents+ of the NOTICE file are for informational purposes only and+ do not modify the License. You may add Your own attribution+ notices within Derivative Works that You distribute, alongside+ or as an addendum to the NOTICE text from the Work, provided+ that such additional attribution notices cannot be construed+ as modifying the License.++ You may add Your own copyright statement to Your modifications and+ may provide additional or different license terms and conditions+ for use, reproduction, or distribution of Your modifications, or+ for any such Derivative Works as a whole, provided Your use,+ reproduction, and distribution of the Work otherwise complies with+ the conditions stated in this License.++ 5. Submission of Contributions. Unless You explicitly state otherwise,+ any Contribution intentionally submitted for inclusion in the Work+ by You to the Licensor shall be under the terms and conditions of+ this License, without any additional terms or conditions.+ Notwithstanding the above, nothing herein shall supersede or modify+ the terms of any separate license agreement you may have executed+ with Licensor regarding such Contributions.++ 6. Trademarks. This License does not grant permission to use the trade+ names, trademarks, service marks, or product names of the Licensor,+ except as required for reasonable and customary use in describing the+ origin of the Work and reproducing the content of the NOTICE file.++ 7. Disclaimer of Warranty. Unless required by applicable law or+ agreed to in writing, Licensor provides the Work (and each+ Contributor provides its Contributions) on an "AS IS" BASIS,+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or+ implied, including, without limitation, any warranties or conditions+ of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A+ PARTICULAR PURPOSE. You are solely responsible for determining the+ appropriateness of using or redistributing the Work and assume any+ risks associated with Your exercise of permissions under this License.++ 8. Limitation of Liability. In no event and under no legal theory,+ whether in tort (including negligence), contract, or otherwise,+ unless required by applicable law (such as deliberate and grossly+ negligent acts) or agreed to in writing, shall any Contributor be+ liable to You for damages, including any direct, indirect, special,+ incidental, or consequential damages of any character arising as a+ result of this License or out of the use or inability to use the+ Work (including but not limited to damages for loss of goodwill,+ work stoppage, computer failure or malfunction, or any and all+ other commercial damages or losses), even if such Contributor+ has been advised of the possibility of such damages.++ 9. Accepting Warranty or Additional Liability. While redistributing+ the Work or Derivative Works thereof, You may choose to offer,+ and charge a fee for, acceptance of support, warranty, indemnity,+ or other liability obligations and/or rights consistent with this+ License. However, in accepting such obligations, You may act only+ on Your own behalf and on Your sole responsibility, not on behalf+ of any other Contributor, and only if You agree to indemnify,+ defend, and hold each Contributor harmless for any liability+ incurred by, or claims asserted against, such Contributor by reason+ of your accepting any such warranty or additional liability.++ END OF TERMS AND CONDITIONS++ APPENDIX: How to apply the Apache License to your work.++ To apply the Apache License to your work, attach the following+ boilerplate notice, with the fields enclosed by brackets "[]"+ replaced with your own identifying information. (Don't include+ the brackets!) The text should be enclosed in the appropriate+ comment syntax for the file format. We also recommend that a+ file or class name and description of purpose be included on the+ same "printed page" as the copyright notice for easier+ identification within third-party archives.++ Copyright [yyyy] [name of copyright owner]++ Licensed under the Apache License, Version 2.0 (the "License");+ you may not use this file except in compliance with the License.+ You may obtain a copy of the License at++ http://www.apache.org/licenses/LICENSE-2.0++ Unless required by applicable law or agreed to in writing, software+ distributed under the License is distributed on an "AS IS" BASIS,+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.+ See the License for the specific language governing permissions and+ limitations under the License.
+ Numeric/LBFGS.hs view
@@ -0,0 +1,259 @@+-- |+-- Module : Numeric.LBFGS+-- Copyright : (c) 2010 Daniël de Kok+-- License : Apache 2+--+--+-- Maintainer : Daniël de Kok <me@danieldk.eu>+-- Stability : experimental+--+-- Binding for the liblbfgs library, much implements the Limited-memory+-- Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method for solving+-- unconstrained minimization problems. The original C library is+-- available from:+--+-- <http://www.chokkan.org/software/liblbfgs/>++module Numeric.LBFGS (LineSearchAlgorithm(..), EvaluateFun,+ ProgressFun, LBFGSResult, lbfgs) where++import Data.Array.Storable (StorableArray,+ unsafeForeignPtrToStorableArray)+import Foreign.C.Types (CDouble, CInt)+import Foreign.ForeignPtr (newForeignPtr_)+import Foreign.Marshal.Alloc (malloc, free)+import Foreign.Ptr (Ptr, freeHaskellFunPtr, nullPtr, plusPtr)+import Foreign.Storable (Storable(..), peek, poke, sizeOf)++import qualified Numeric.LBFGS.Raw as R+import Numeric.LBFGS.Raw (CEvaluateFun, CProgressFun, CLBFGSParameter(..),+ defaultCParam, CLBFGSResult(..),+ c_lbfgs_malloc, c_lbfgs_free,+ c_lbfgs_evaluate_t_wrap, c_lbfgs_progress_t_wrap,+ c_lbfgs+ )++-- |+-- Various line search algorithms. Wolfe backtracking algorithms require+-- a coefficient.+data LineSearchAlgorithm = DefaultLineSearch+ | MoreThuente+ | BacktrackingArmijo+ | Backtracking+ | BacktrackingWolfe {coeff :: Double }+ | BacktrackingStrongWolfe {coeff :: Double }++mergeLineSearchAlgorithm :: CLBFGSParameter -> LineSearchAlgorithm ->+ CLBFGSParameter+mergeLineSearchAlgorithm p DefaultLineSearch =+ p {R.linesearch = R.defaultLineSearch}+mergeLineSearchAlgorithm p MoreThuente =+ p { R.linesearch = R.moreThuente }+mergeLineSearchAlgorithm p BacktrackingArmijo =+ p { R.linesearch = R.backtrackingArmijo }+mergeLineSearchAlgorithm p Backtracking =+ p { R.linesearch = R.backtracking }+mergeLineSearchAlgorithm p (BacktrackingWolfe c) =+ p { R.linesearch = R.backtrackingWolfe,+ R.wolfe = realToFrac c }+mergeLineSearchAlgorithm p (BacktrackingStrongWolfe c) =+ p { R.linesearch = R.backtrackingStrongWolfe,+ R.wolfe = realToFrac c }++withParam :: LineSearchAlgorithm -> CLBFGSParameter+withParam lineSearch =+ mergeLineSearchAlgorithm defaultCParam lineSearch+++data LBFGSResult+ = Success+ | Stop+ | AlreadyMinimized+ | UnknownError+ | LogicError+ | OutOfMemory+ | Canceled+ | InvalidN+ | InvalidNSSE+ | InvalidXSSE+ | InvalidEpsilon+ | InvalidTestPeriod+ | InvalidDelta+ | InvalidLineSearch+ | InvalidMinStep+ | InvalidMaxStep+ | InvalidFtol+ | InvalidWolfe+ | InvalidGtol+ | InvalidXtol+ | InvalidMaxLineSearch+ | InvalidOrthantwise+ | InvalidOrthantwiseStart+ | InvalidOrthantwiseEnd+ | OutOfInterval+ | IncorrectTMinMax+ | RoundingError+ | MinimumStep+ | MaximumStep+ | MaximumLineSearch+ | MaximumIteration+ | WidthTooSmall+ | InvalidParameters+ | IncreaseGradient+ deriving (Eq, Show)++deriveResult :: CLBFGSResult -> LBFGSResult+deriveResult r+ | r == R.lbfgsSuccess = Success+ | r == R.lbfgsStop = Stop+ | r == R.lbfgsAlreadyMinimized = AlreadyMinimized+ | r == R.lbfgserrUnknownerror = UnknownError+ | r == R.lbfgserrLogicerror = LogicError+ | r == R.lbfgserrOutofmemory = OutOfMemory+ | r == R.lbfgserrCanceled = Canceled+ | r == R.lbfgserrInvalidN = InvalidN+ | r == R.lbfgserrInvalidNSse = InvalidNSSE+ | r == R.lbfgserrInvalidXSse = InvalidXSSE+ | r == R.lbfgserrInvalidEpsilon = InvalidEpsilon+ | r == R.lbfgserrInvalidTestperiod = InvalidTestPeriod+ | r == R.lbfgserrInvalidDelta = InvalidDelta+ | r == R.lbfgserrInvalidLinesearch = InvalidLineSearch+ | r == R.lbfgserrInvalidMinstep = InvalidMinStep+ | r == R.lbfgserrInvalidMaxstep = InvalidMaxStep+ | r == R.lbfgserrInvalidFtol = InvalidFtol+ | r == R.lbfgserrInvalidWolfe = InvalidWolfe+ | r == R.lbfgserrInvalidGtol = InvalidGtol+ | r == R.lbfgserrInvalidXtol = InvalidXtol+ | r == R.lbfgserrInvalidMaxlinesearch = InvalidMaxLineSearch+ | r == R.lbfgserrInvalidOrthantwise = InvalidOrthantwise+ | r == R.lbfgserrInvalidOrthantwiseStart = InvalidOrthantwiseStart+ | r == R.lbfgserrInvalidOrthantwiseEnd = InvalidOrthantwiseEnd+ | r == R.lbfgserrOutofinterval = OutOfInterval+ | r == R.lbfgserrIncorrectTminmax = IncorrectTMinMax+ | r == R.lbfgserrRoundingError = RoundingError+ | r == R.lbfgserrMinimumstep = MinimumStep+ | r == R.lbfgserrMaximumstep = MaximumStep+ | r == R.lbfgserrMaximumlinesearch = MaximumLineSearch+ | r == R.lbfgserrMaximumiteration = MaximumIteration+ | r == R.lbfgserrWidthtoosmall = WidthTooSmall+ | r == R.lbfgserrInvalidparameters = InvalidParameters+ | r == R.lbfgserrIncreasegradient = IncreaseGradient++cDoublePlusPtr :: Ptr CDouble -> Int -> Ptr CDouble+cDoublePlusPtr ptr n = plusPtr ptr (n * sizeOf (undefined :: CDouble))++listToVector :: [Double] -> IO (CInt, Ptr CDouble)+listToVector l = do+ v <- c_lbfgs_malloc n+ copyList l v+ return (n, v)+ where n = fromIntegral . length $ l++copyList :: [Double] -> Ptr CDouble -> IO ()+copyList [] _ = return ()+copyList l p = do+ poke p $ realToFrac $ head l+ copyList (tail l) (cDoublePlusPtr p 1)+++freeVector :: Ptr CDouble -> IO ()+freeVector = c_lbfgs_free++vectorToList :: CInt -> Ptr CDouble -> IO ([Double])+vectorToList cn p = vectorToList_ p (cDoublePlusPtr p n) []+ where n = fromIntegral cn++vectorToList_ :: Ptr CDouble -> Ptr CDouble -> [Double] -> IO ([Double])+vectorToList_ pStart pCur l+ | pCur >= pStart = do+ cval <- peek pCur+ let val = realToFrac cval+ vectorToList_ pStart (cDoublePlusPtr pCur (-1)) (val:l)+ | otherwise = return l+++-- |+-- Type signature for the objective function and gradient evaluations.+type EvaluateFun a =+ a -- ^ Instance data+ -> StorableArray Int CDouble -- ^ Current variables (should not be+ -- modified by the function)+ -> StorableArray Int CDouble -- ^ Gradients+ -> CInt -- ^ Number of variables+ -> CDouble -- ^ Step of the line search algorithm+ -> IO (CDouble) -- ^ Value of the objective function++wrapEvaluateFun :: (Storable a) => EvaluateFun a -> Ptr a -> Ptr CDouble ->+ Ptr CDouble -> CInt -> CDouble -> IO (CDouble)+wrapEvaluateFun fun inst x g n step = do+ let nInt = fromIntegral n+ instV <- peek inst+ xFp <- newForeignPtr_ x+ xArr <- unsafeForeignPtrToStorableArray xFp (0, nInt - 1)+ gFp <- newForeignPtr_ g+ gArr <- unsafeForeignPtrToStorableArray gFp (0, nInt - 1)+ fun instV xArr gArr n step++-- |+-- Type signature for a function reporting on the progress of the+-- optimization.+type ProgressFun a =+ a -- ^ Instance data+ -> StorableArray Int CDouble -- ^ Variables (should not be modified+ -- by the function)+ -> StorableArray Int CDouble -- ^ Gradients (should not be modified+ -- by the function)+ -> CDouble -- ^ Value of the objective function+ -> CDouble -- ^ Euclidean norm of the variables+ -> CDouble -- ^ Eucledian norm of the gradients+ -> CDouble -- ^ Step of the line search algorithm+ -> CInt -- ^ Number of variables+ -> CInt -- ^ Iteration count+ -> CInt -- ^ Number of evaluations for this iteration+ -> IO (CInt) -- ^ Return zero to continue the evaluation,+ -- non-zero otherwise++wrapProgressFun :: (Storable a) => ProgressFun a -> Ptr a -> Ptr CDouble ->+ Ptr CDouble-> CDouble -> CDouble -> CDouble -> CDouble ->+ CInt -> CInt -> CInt -> IO (CInt)+wrapProgressFun fun inst x g fx xn gn step n k ls = do+ let nInt = fromIntegral n+ instV <- peek inst+ xFp <- newForeignPtr_ x+ xArr <- unsafeForeignPtrToStorableArray xFp (0, nInt - 1)+ gFp <- newForeignPtr_ g+ gArr <- unsafeForeignPtrToStorableArray gFp (0, nInt - 1)+ fun instV xArr gArr fx xn gn step n k ls++-- |+-- Start a L-BFGS optimization. The initial variables should be+-- provided as a list of doubles.+lbfgs :: (Storable a) =>+ LineSearchAlgorithm -- ^ The line search algorithm+ -> EvaluateFun a -- ^ Objective function+ -> ProgressFun a -- ^ Progress report function+ -> a -- ^ Instance data+ -> [Double] -- ^ Initial variable values+ -> IO(LBFGSResult, [Double]) -- ^ Result and variable values+lbfgs ls evalFun progressFun inst p = lbfgs_ ls (wrapEvaluateFun evalFun)+ (wrapProgressFun progressFun) inst p++lbfgs_ :: (Storable a) => LineSearchAlgorithm -> CEvaluateFun a ->+ CProgressFun a -> a -> [Double] -> IO(LBFGSResult, [Double])+lbfgs_ ls evalFun progressFun inst p = do+ (n, pVec) <- listToVector p+ let param = withParam ls+ instP <- malloc+ poke instP inst+ paramP <- malloc+ poke paramP param+ evalW <- c_lbfgs_evaluate_t_wrap evalFun+ progressW <- c_lbfgs_progress_t_wrap progressFun+ r <- c_lbfgs n pVec nullPtr evalW progressW instP paramP+ freeHaskellFunPtr progressW+ freeHaskellFunPtr evalW+ free paramP+ free instP+ freeVector pVec+ rl <- vectorToList n pVec+ return (deriveResult $ CLBFGSResult r, rl)
+ Numeric/LBFGS/Raw.hsc view
@@ -0,0 +1,186 @@+{-# LANGUAGE ForeignFunctionInterface, GeneralizedNewtypeDeriving #-}++#include "lbfgs.h"+#let alignment t = "%lu", (unsigned long)offsetof(struct {char x__; t (y__); }, y__)++module Numeric.LBFGS.Raw (CLineSearchAlgorithm, CLBFGSParameter(..),+ CEvaluateFun, CProgressFun,+ defaultCParam, c_lbfgs, c_lbfgs_malloc,+ c_lbfgs_free, c_lbfgs_evaluate_t_wrap,+ c_lbfgs_progress_t_wrap,++ defaultLineSearch, moreThuente, backtrackingArmijo,+ backtracking, backtrackingWolfe,+ backtrackingStrongWolfe,++ CLBFGSResult(..),+ lbfgsSuccess,+ lbfgsConvergence,+ lbfgsStop,+ lbfgsAlreadyMinimized,+ lbfgserrUnknownerror,+ lbfgserrLogicerror,+ lbfgserrOutofmemory,+ lbfgserrCanceled,+ lbfgserrInvalidN,+ lbfgserrInvalidNSse,+ lbfgserrInvalidXSse,+ lbfgserrInvalidEpsilon,+ lbfgserrInvalidTestperiod,+ lbfgserrInvalidDelta,+ lbfgserrInvalidLinesearch,+ lbfgserrInvalidMinstep,+ lbfgserrInvalidMaxstep,+ lbfgserrInvalidFtol,+ lbfgserrInvalidWolfe,+ lbfgserrInvalidGtol,+ lbfgserrInvalidXtol,+ lbfgserrInvalidMaxlinesearch,+ lbfgserrInvalidOrthantwise,+ lbfgserrInvalidOrthantwiseStart,+ lbfgserrInvalidOrthantwiseEnd,+ lbfgserrOutofinterval,+ lbfgserrIncorrectTminmax,+ lbfgserrRoundingError,+ lbfgserrMinimumstep,+ lbfgserrMaximumstep,+ lbfgserrMaximumlinesearch,+ lbfgserrMaximumiteration,+ lbfgserrWidthtoosmall,+ lbfgserrInvalidparameters,+ lbfgserrIncreasegradient++) where++import Foreign.Storable (Storable(..))+import Foreign.C.Types (CDouble, CInt)+import Foreign.Ptr (FunPtr, Ptr)++newtype CLineSearchAlgorithm =+ CLineSearchAlgorithm { unCLineSearchAlgorithm :: CInt }+ deriving (Storable, Show)++#{enum CLineSearchAlgorithm, CLineSearchAlgorithm,+ defaultLineSearch = LBFGS_LINESEARCH_DEFAULT,+ moreThuente = LBFGS_LINESEARCH_MORETHUENTE,+ backtrackingArmijo = LBFGS_LINESEARCH_BACKTRACKING_ARMIJO,+ backtracking = LBFGS_LINESEARCH_BACKTRACKING,+ backtrackingWolfe = LBFGS_LINESEARCH_BACKTRACKING_WOLFE,+ backtrackingStrongWolfe = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE+}++newtype CLBFGSResult =+ CLBFGSResult { unCLBFGSResult :: CInt }+ deriving (Eq, Show)++#{enum CLBFGSResult, CLBFGSResult,+ LBFGS_SUCCESS, LBFGS_CONVERGENCE, LBFGS_STOP,+ LBFGS_ALREADY_MINIMIZED, LBFGSERR_UNKNOWNERROR, LBFGSERR_LOGICERROR,+ LBFGSERR_OUTOFMEMORY, LBFGSERR_CANCELED, LBFGSERR_INVALID_N,+ LBFGSERR_INVALID_N_SSE, LBFGSERR_INVALID_X_SSE,+ LBFGSERR_INVALID_EPSILON, LBFGSERR_INVALID_TESTPERIOD,+ LBFGSERR_INVALID_DELTA, LBFGSERR_INVALID_LINESEARCH,+ LBFGSERR_INVALID_MINSTEP, LBFGSERR_INVALID_MAXSTEP,+ LBFGSERR_INVALID_FTOL, LBFGSERR_INVALID_WOLFE,+ LBFGSERR_INVALID_GTOL, LBFGSERR_INVALID_XTOL,+ LBFGSERR_INVALID_MAXLINESEARCH, LBFGSERR_INVALID_ORTHANTWISE,+ LBFGSERR_INVALID_ORTHANTWISE_START,+ LBFGSERR_INVALID_ORTHANTWISE_END, LBFGSERR_OUTOFINTERVAL,+ LBFGSERR_INCORRECT_TMINMAX, LBFGSERR_ROUNDING_ERROR,+ LBFGSERR_MINIMUMSTEP, LBFGSERR_MAXIMUMSTEP,+ LBFGSERR_MAXIMUMLINESEARCH, LBFGSERR_MAXIMUMITERATION,+ LBFGSERR_WIDTHTOOSMALL, LBFGSERR_INVALIDPARAMETERS,+ LBFGSERR_INCREASEGRADIENT }++data CLBFGSParameter = CLBFGSParameter {+ m :: CInt,+ epsilon :: CDouble,+ past :: CInt,+ delta :: CDouble,+ max_iterations :: CInt,+ linesearch :: CLineSearchAlgorithm,+ max_linesearch :: CInt,+ min_step :: CDouble,+ max_step :: CDouble,+ ftol :: CDouble,+ wolfe :: CDouble,+ gtol :: CDouble,+ xtol :: CDouble,+ orthantwise_c :: CDouble,+ orthantwise_start :: CDouble,+ orthantwise_end :: CDouble+} deriving Show++defaultCParam :: CLBFGSParameter+defaultCParam = CLBFGSParameter 6 1e-5 0 1e-5 0 defaultLineSearch 40 1e-20+ 1e20 1e-4 0.9 0.9 1.0e-16 0.0 0.0 (-1.0)++instance Storable CLBFGSParameter where+ sizeOf _ = #{size lbfgs_parameter_t}+ alignment _ = #{alignment lbfgs_parameter_t}+ peek ptr = do+ m <- (#peek lbfgs_parameter_t, m) ptr+ epsilon <- (#peek lbfgs_parameter_t, epsilon) ptr+ past <- (#peek lbfgs_parameter_t, past) ptr+ delta <- (#peek lbfgs_parameter_t, delta) ptr+ max_iterations <- (#peek lbfgs_parameter_t, max_iterations) ptr+ linesearch <- (#peek lbfgs_parameter_t, linesearch) ptr+ max_linesearch <- (#peek lbfgs_parameter_t, max_linesearch) ptr+ min_step <- (#peek lbfgs_parameter_t, min_step) ptr+ max_step <- (#peek lbfgs_parameter_t, max_step) ptr+ ftol <- (#peek lbfgs_parameter_t, ftol) ptr+ wolfe <- (#peek lbfgs_parameter_t, wolfe) ptr+ gtol <- (#peek lbfgs_parameter_t, gtol) ptr+ xtol <- (#peek lbfgs_parameter_t, xtol) ptr+ orthantwise_c <- (#peek lbfgs_parameter_t, orthantwise_c) ptr+ orthantwise_start <- (#peek lbfgs_parameter_t, orthantwise_start) ptr+ orthantwise_end <- (#peek lbfgs_parameter_t, orthantwise_end) ptr+ return $ CLBFGSParameter m epsilon past delta max_iterations+ linesearch max_linesearch min_step max_step+ ftol wolfe gtol xtol orthantwise_c+ orthantwise_start orthantwise_end+ poke ptr (CLBFGSParameter m epsilon past delta max_iterations+ linesearch max_linesearch min_step max_step+ ftol wolfe gtol xtol orthantwise_c+ orthantwise_start orthantwise_end+ ) = do+ (#poke lbfgs_parameter_t, m) ptr m+ (#poke lbfgs_parameter_t, epsilon) ptr epsilon+ (#poke lbfgs_parameter_t, past) ptr past+ (#poke lbfgs_parameter_t, delta) ptr delta+ (#poke lbfgs_parameter_t, max_iterations) ptr max_iterations+ (#poke lbfgs_parameter_t, linesearch) ptr linesearch+ (#poke lbfgs_parameter_t, max_linesearch) ptr max_linesearch+ (#poke lbfgs_parameter_t, min_step) ptr min_step+ (#poke lbfgs_parameter_t, max_step) ptr max_step+ (#poke lbfgs_parameter_t, ftol) ptr ftol+ (#poke lbfgs_parameter_t, wolfe) ptr wolfe+ (#poke lbfgs_parameter_t, gtol) ptr gtol+ (#poke lbfgs_parameter_t, xtol) ptr xtol+ (#poke lbfgs_parameter_t, orthantwise_c) ptr orthantwise_c+ (#poke lbfgs_parameter_t, orthantwise_start) ptr orthantwise_start+ (#poke lbfgs_parameter_t, orthantwise_end) ptr orthantwise_end++type CEvaluateFun a = (Ptr a -> Ptr CDouble -> Ptr CDouble -> CInt ->+ CDouble -> IO (CDouble))++type CProgressFun a = (Ptr a -> Ptr CDouble -> Ptr CDouble -> CDouble ->+ CDouble -> CDouble -> CDouble -> CInt -> CInt ->+ CInt -> IO (CInt))++foreign import ccall "wrapper"+ c_lbfgs_evaluate_t_wrap :: CEvaluateFun a -> IO (FunPtr (CEvaluateFun a))++foreign import ccall "wrapper"+ c_lbfgs_progress_t_wrap :: CProgressFun a -> IO (FunPtr (CProgressFun a))++foreign import ccall safe "lbfgs.h lbfgs" c_lbfgs ::+ CInt -> Ptr CDouble -> Ptr CDouble -> FunPtr (CEvaluateFun a) ->+ FunPtr (CProgressFun a) -> Ptr a -> Ptr (CLBFGSParameter) -> IO (CInt)++foreign import ccall unsafe "lbfgs.h lbfgs_malloc" c_lbfgs_malloc ::+ CInt -> IO (Ptr CDouble)++foreign import ccall unsafe "lbfgs.h lbfgs_free" c_lbfgs_free ::+ Ptr CDouble -> IO ()+
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ cbits/lbfgs.c view
@@ -0,0 +1,1371 @@+/*+ * Limited memory BFGS (L-BFGS).+ *+ * Copyright (c) 1990, Jorge Nocedal+ * Copyright (c) 2007-2010 Naoaki Okazaki+ * All rights reserved.+ *+ * Permission is hereby granted, free of charge, to any person obtaining a copy+ * of this software and associated documentation files (the "Software"), to deal+ * in the Software without restriction, including without limitation the rights+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell+ * copies of the Software, and to permit persons to whom the Software is+ * furnished to do so, subject to the following conditions:+ *+ * The above copyright notice and this permission notice shall be included in+ * all copies or substantial portions of the Software.+ *+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN+ * THE SOFTWARE.+ */++/* $Id: lbfgs.c 65 2010-01-29 12:19:16Z naoaki $ */++/*+This library is a C port of the FORTRAN implementation of Limited-memory+Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method written by Jorge Nocedal.+The original FORTRAN source code is available at:+http://www.ece.northwestern.edu/~nocedal/lbfgs.html++The L-BFGS algorithm is described in:+ - Jorge Nocedal.+ Updating Quasi-Newton Matrices with Limited Storage.+ <i>Mathematics of Computation</i>, Vol. 35, No. 151, pp. 773--782, 1980.+ - Dong C. Liu and Jorge Nocedal.+ On the limited memory BFGS method for large scale optimization.+ <i>Mathematical Programming</i> B, Vol. 45, No. 3, pp. 503-528, 1989.++The line search algorithms used in this implementation are described in:+ - John E. Dennis and Robert B. Schnabel.+ <i>Numerical Methods for Unconstrained Optimization and Nonlinear+ Equations</i>, Englewood Cliffs, 1983.+ - Jorge J. More and David J. Thuente.+ Line search algorithm with guaranteed sufficient decrease.+ <i>ACM Transactions on Mathematical Software (TOMS)</i>, Vol. 20, No. 3,+ pp. 286-307, 1994.++This library also implements Orthant-Wise Limited-memory Quasi-Newton (OWL-QN)+method presented in:+ - Galen Andrew and Jianfeng Gao.+ Scalable training of L1-regularized log-linear models.+ In <i>Proceedings of the 24th International Conference on Machine+ Learning (ICML 2007)</i>, pp. 33-40, 2007.++I would like to thank the original author, Jorge Nocedal, who has been+distributing the effieicnt and explanatory implementation in an open source+licence.+*/++#ifdef HAVE_CONFIG_H+#include <config.h>+#endif/*HAVE_CONFIG_H*/++#include <stdio.h>+#include <stdlib.h>+#include <math.h>++#include <lbfgs.h>++#ifdef _MSC_VER+#define inline __inline+typedef unsigned int uint32_t;+#endif/*_MSC_VER*/++#if defined(USE_SSE) && defined(__SSE2__) && LBFGS_FLOAT == 64+/* Use SSE2 optimization for 64bit double precision. */+#include "arithmetic_sse_double.h"++#elif defined(USE_SSE) && defined(__SSE__) && LBFGS_FLOAT == 32+/* Use SSE optimization for 32bit float precision. */+#include "arithmetic_sse_float.h"++#else+/* No CPU specific optimization. */+#include "arithmetic_ansi.h"++#endif++#define min2(a, b) ((a) <= (b) ? (a) : (b))+#define max2(a, b) ((a) >= (b) ? (a) : (b))+#define max3(a, b, c) max2(max2((a), (b)), (c));++struct tag_callback_data {+ int n;+ void *instance;+ lbfgs_evaluate_t proc_evaluate;+ lbfgs_progress_t proc_progress;+};+typedef struct tag_callback_data callback_data_t;++struct tag_iteration_data {+ lbfgsfloatval_t alpha;+ lbfgsfloatval_t *s; /* [n] */+ lbfgsfloatval_t *y; /* [n] */+ lbfgsfloatval_t ys; /* vecdot(y, s) */+};+typedef struct tag_iteration_data iteration_data_t;++static const lbfgs_parameter_t _defparam = {+ 6, 1e-5, 0, 1e-5,+ 0, LBFGS_LINESEARCH_DEFAULT, 40,+ 1e-20, 1e20, 1e-4, 0.9, 0.9, 1.0e-16,+ 0.0, 0, -1,+};++/* Forward function declarations. */++typedef int (*line_search_proc)(+ int n,+ lbfgsfloatval_t *x,+ lbfgsfloatval_t *f,+ lbfgsfloatval_t *g,+ lbfgsfloatval_t *s,+ lbfgsfloatval_t *stp,+ const lbfgsfloatval_t* xp,+ const lbfgsfloatval_t* gp,+ lbfgsfloatval_t *wa,+ callback_data_t *cd,+ const lbfgs_parameter_t *param+ );+ +static int line_search_backtracking(+ int n,+ lbfgsfloatval_t *x,+ lbfgsfloatval_t *f,+ lbfgsfloatval_t *g,+ lbfgsfloatval_t *s,+ lbfgsfloatval_t *stp,+ const lbfgsfloatval_t* xp,+ const lbfgsfloatval_t* gp,+ lbfgsfloatval_t *wa,+ callback_data_t *cd,+ const lbfgs_parameter_t *param+ );++static int line_search_backtracking_owlqn(+ int n,+ lbfgsfloatval_t *x,+ lbfgsfloatval_t *f,+ lbfgsfloatval_t *g,+ lbfgsfloatval_t *s,+ lbfgsfloatval_t *stp,+ const lbfgsfloatval_t* xp,+ const lbfgsfloatval_t* gp,+ lbfgsfloatval_t *wp,+ callback_data_t *cd,+ const lbfgs_parameter_t *param+ );++static int line_search_morethuente(+ int n,+ lbfgsfloatval_t *x,+ lbfgsfloatval_t *f,+ lbfgsfloatval_t *g,+ lbfgsfloatval_t *s,+ lbfgsfloatval_t *stp,+ const lbfgsfloatval_t* xp,+ const lbfgsfloatval_t* gp,+ lbfgsfloatval_t *wa,+ callback_data_t *cd,+ const lbfgs_parameter_t *param+ );++static int update_trial_interval(+ lbfgsfloatval_t *x,+ lbfgsfloatval_t *fx,+ lbfgsfloatval_t *dx,+ lbfgsfloatval_t *y,+ lbfgsfloatval_t *fy,+ lbfgsfloatval_t *dy,+ lbfgsfloatval_t *t,+ lbfgsfloatval_t *ft,+ lbfgsfloatval_t *dt,+ const lbfgsfloatval_t tmin,+ const lbfgsfloatval_t tmax,+ int *brackt+ );++static lbfgsfloatval_t owlqn_x1norm(+ const lbfgsfloatval_t* x,+ const int start,+ const int n+ );++static void owlqn_pseudo_gradient(+ lbfgsfloatval_t* pg,+ const lbfgsfloatval_t* x,+ const lbfgsfloatval_t* g,+ const int n,+ const lbfgsfloatval_t c,+ const int start,+ const int end+ );++static void owlqn_project(+ lbfgsfloatval_t* d,+ const lbfgsfloatval_t* sign,+ const int start,+ const int end+ );+++#if defined(USE_SSE) && (defined(__SSE__) || defined(__SSE2__))+static int round_out_variables(int n)+{+ n += 7;+ n /= 8;+ n *= 8;+ return n;+}+#endif/*defined(USE_SSE)*/++lbfgsfloatval_t* lbfgs_malloc(int n)+{+#if defined(USE_SSE) && (defined(__SSE__) || defined(__SSE2__))+ n = round_out_variables(n);+#endif/*defined(USE_SSE)*/+ return (lbfgsfloatval_t*)vecalloc(sizeof(lbfgsfloatval_t) * n);+}++void lbfgs_free(lbfgsfloatval_t *x)+{+ vecfree(x);+}++void lbfgs_parameter_init(lbfgs_parameter_t *param)+{+ memcpy(param, &_defparam, sizeof(*param));+}++int lbfgs(+ int n,+ lbfgsfloatval_t *x,+ lbfgsfloatval_t *ptr_fx,+ lbfgs_evaluate_t proc_evaluate,+ lbfgs_progress_t proc_progress,+ void *instance,+ lbfgs_parameter_t *_param+ )+{+ int ret;+ int i, j, k, ls, end, bound;+ lbfgsfloatval_t step;++ /* Constant parameters and their default values. */+ lbfgs_parameter_t param = (_param != NULL) ? (*_param) : _defparam;+ const int m = param.m;++ lbfgsfloatval_t *xp = NULL;+ lbfgsfloatval_t *g = NULL, *gp = NULL, *pg = NULL;+ lbfgsfloatval_t *d = NULL, *w = NULL, *pf = NULL;+ iteration_data_t *lm = NULL, *it = NULL;+ lbfgsfloatval_t ys, yy;+ lbfgsfloatval_t xnorm, gnorm, beta;+ lbfgsfloatval_t fx = 0.;+ lbfgsfloatval_t rate = 0.;+ line_search_proc linesearch = line_search_morethuente;++ /* Construct a callback data. */+ callback_data_t cd;+ cd.n = n;+ cd.instance = instance;+ cd.proc_evaluate = proc_evaluate;+ cd.proc_progress = proc_progress;++#if defined(USE_SSE) && (defined(__SSE__) || defined(__SSE2__))+ /* Round out the number of variables. */+ n = round_out_variables(n);+#endif/*defined(USE_SSE)*/++ /* Check the input parameters for errors. */+ if (n <= 0) {+ return LBFGSERR_INVALID_N;+ }+#if defined(USE_SSE) && (defined(__SSE__) || defined(__SSE2__))+ if (n % 8 != 0) {+ return LBFGSERR_INVALID_N_SSE;+ }+ if (((unsigned short)x & 0x000F) != 0) {+ return LBFGSERR_INVALID_X_SSE;+ }+#endif/*defined(USE_SSE)*/+ if (param.epsilon < 0.) {+ return LBFGSERR_INVALID_EPSILON;+ }+ if (param.past < 0) {+ return LBFGSERR_INVALID_TESTPERIOD;+ }+ if (param.delta < 0.) {+ return LBFGSERR_INVALID_DELTA;+ }+ if (param.min_step < 0.) {+ return LBFGSERR_INVALID_MINSTEP;+ }+ if (param.max_step < param.min_step) {+ return LBFGSERR_INVALID_MAXSTEP;+ }+ if (param.ftol < 0.) {+ return LBFGSERR_INVALID_FTOL;+ }+ if (param.linesearch == LBFGS_LINESEARCH_BACKTRACKING_WOLFE ||+ param.linesearch == LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE) {+ if (param.wolfe <= param.ftol || 1. <= param.wolfe) {+ return LBFGSERR_INVALID_WOLFE;+ }+ }+ if (param.gtol < 0.) {+ return LBFGSERR_INVALID_GTOL;+ }+ if (param.xtol < 0.) {+ return LBFGSERR_INVALID_XTOL;+ }+ if (param.max_linesearch <= 0) {+ return LBFGSERR_INVALID_MAXLINESEARCH;+ }+ if (param.orthantwise_c < 0.) {+ return LBFGSERR_INVALID_ORTHANTWISE;+ }+ if (param.orthantwise_start < 0 || n < param.orthantwise_start) {+ return LBFGSERR_INVALID_ORTHANTWISE_START;+ }+ if (param.orthantwise_end < 0) {+ param.orthantwise_end = n;+ }+ if (n < param.orthantwise_end) {+ return LBFGSERR_INVALID_ORTHANTWISE_END;+ }+ if (param.orthantwise_c != 0.) {+ switch (param.linesearch) {+ case LBFGS_LINESEARCH_BACKTRACKING:+ linesearch = line_search_backtracking_owlqn;+ break;+ default:+ /* Only the backtracking method is available. */+ return LBFGSERR_INVALID_LINESEARCH;+ }+ } else {+ switch (param.linesearch) {+ case LBFGS_LINESEARCH_MORETHUENTE:+ linesearch = line_search_morethuente;+ break;+ case LBFGS_LINESEARCH_BACKTRACKING_ARMIJO:+ case LBFGS_LINESEARCH_BACKTRACKING_WOLFE:+ case LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE:+ linesearch = line_search_backtracking;+ break;+ default:+ return LBFGSERR_INVALID_LINESEARCH;+ }+ }++ /* Allocate working space. */+ xp = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));+ g = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));+ gp = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));+ d = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));+ w = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));+ if (xp == NULL || g == NULL || gp == NULL || d == NULL || w == NULL) {+ ret = LBFGSERR_OUTOFMEMORY;+ goto lbfgs_exit;+ }++ if (param.orthantwise_c != 0.) {+ /* Allocate working space for OW-LQN. */+ pg = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));+ if (pg == NULL) {+ ret = LBFGSERR_OUTOFMEMORY;+ goto lbfgs_exit;+ }+ }++ /* Allocate limited memory storage. */+ lm = (iteration_data_t*)vecalloc(m * sizeof(iteration_data_t));+ if (lm == NULL) {+ ret = LBFGSERR_OUTOFMEMORY;+ goto lbfgs_exit;+ }++ /* Initialize the limited memory. */+ for (i = 0;i < m;++i) {+ it = &lm[i];+ it->alpha = 0;+ it->ys = 0;+ it->s = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));+ it->y = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t));+ if (it->s == NULL || it->y == NULL) {+ ret = LBFGSERR_OUTOFMEMORY;+ goto lbfgs_exit;+ }+ }++ /* Allocate an array for storing previous values of the objective function. */+ if (0 < param.past) {+ pf = (lbfgsfloatval_t*)vecalloc(param.past * sizeof(lbfgsfloatval_t));+ }++ /* Evaluate the function value and its gradient. */+ fx = cd.proc_evaluate(cd.instance, x, g, cd.n, 0);+ if (0. != param.orthantwise_c) {+ /* Compute the L1 norm of the variable and add it to the object value. */+ xnorm = owlqn_x1norm(x, param.orthantwise_start, param.orthantwise_end);+ fx += xnorm * param.orthantwise_c;+ owlqn_pseudo_gradient(+ pg, x, g, n,+ param.orthantwise_c, param.orthantwise_start, param.orthantwise_end+ );+ }++ /* Store the initial value of the objective function. */+ if (pf != NULL) {+ pf[0] = fx;+ }++ /*+ Compute the direction;+ we assume the initial hessian matrix H_0 as the identity matrix.+ */+ if (param.orthantwise_c == 0.) {+ vecncpy(d, g, n);+ } else {+ vecncpy(d, pg, n);+ }++ /*+ Make sure that the initial variables are not a minimizer.+ */+ vec2norm(&xnorm, x, n);+ if (param.orthantwise_c == 0.) {+ vec2norm(&gnorm, g, n);+ } else {+ vec2norm(&gnorm, pg, n);+ }+ if (xnorm < 1.0) xnorm = 1.0;+ if (gnorm / xnorm <= param.epsilon) {+ ret = LBFGS_ALREADY_MINIMIZED;+ goto lbfgs_exit;+ }++ /* Compute the initial step:+ step = 1.0 / sqrt(vecdot(d, d, n))+ */+ vec2norminv(&step, d, n);++ k = 1;+ end = 0;+ for (;;) {+ /* Store the current position and gradient vectors. */+ veccpy(xp, x, n);+ veccpy(gp, g, n);++ /* Search for an optimal step. */+ if (param.orthantwise_c == 0.) {+ ls = linesearch(n, x, &fx, g, d, &step, xp, gp, w, &cd, ¶m);+ } else {+ ls = linesearch(n, x, &fx, g, d, &step, xp, pg, w, &cd, ¶m);+ owlqn_pseudo_gradient(+ pg, x, g, n,+ param.orthantwise_c, param.orthantwise_start, param.orthantwise_end+ );+ }+ if (ls < 0) {+ /* Revert to the previous point. */+ veccpy(x, xp, n);+ veccpy(g, gp, n);+ ret = ls;+ goto lbfgs_exit;+ }++ /* Compute x and g norms. */+ vec2norm(&xnorm, x, n);+ if (param.orthantwise_c == 0.) {+ vec2norm(&gnorm, g, n);+ } else {+ vec2norm(&gnorm, pg, n);+ }++ /* Report the progress. */+ if (cd.proc_progress) {+ if (ret = cd.proc_progress(cd.instance, x, g, fx, xnorm, gnorm, step, cd.n, k, ls)) {+ goto lbfgs_exit;+ }+ }++ /*+ Convergence test.+ The criterion is given by the following formula:+ |g(x)| / \max(1, |x|) < \epsilon+ */+ if (xnorm < 1.0) xnorm = 1.0;+ if (gnorm / xnorm <= param.epsilon) {+ /* Convergence. */+ ret = LBFGS_SUCCESS;+ break;+ }++ /*+ Test for stopping criterion.+ The criterion is given by the following formula:+ (f(past_x) - f(x)) / f(x) < \delta+ */+ if (pf != NULL) {+ /* We don't test the stopping criterion while k < past. */+ if (param.past <= k) {+ /* Compute the relative improvement from the past. */+ rate = (pf[k % param.past] - fx) / fx;++ /* The stopping criterion. */+ if (rate < param.delta) {+ ret = LBFGS_STOP;+ break;+ }+ }++ /* Store the current value of the objective function. */+ pf[k % param.past] = fx;+ }++ if (param.max_iterations != 0 && param.max_iterations < k+1) {+ /* Maximum number of iterations. */+ ret = LBFGSERR_MAXIMUMITERATION;+ break;+ }++ /*+ Update vectors s and y:+ s_{k+1} = x_{k+1} - x_{k} = \step * d_{k}.+ y_{k+1} = g_{k+1} - g_{k}.+ */+ it = &lm[end];+ vecdiff(it->s, x, xp, n);+ vecdiff(it->y, g, gp, n);++ /*+ Compute scalars ys and yy:+ ys = y^t \cdot s = 1 / \rho.+ yy = y^t \cdot y.+ Notice that yy is used for scaling the hessian matrix H_0 (Cholesky factor).+ */+ vecdot(&ys, it->y, it->s, n);+ vecdot(&yy, it->y, it->y, n);+ it->ys = ys;++ /*+ Recursive formula to compute dir = -(H \cdot g).+ This is described in page 779 of:+ Jorge Nocedal.+ Updating Quasi-Newton Matrices with Limited Storage.+ Mathematics of Computation, Vol. 35, No. 151,+ pp. 773--782, 1980.+ */+ bound = (m <= k) ? m : k;+ ++k;+ end = (end + 1) % m;++ /* Compute the steepest direction. */+ if (param.orthantwise_c == 0.) {+ /* Compute the negative of gradients. */+ vecncpy(d, g, n);+ } else {+ vecncpy(d, pg, n);+ }++ j = end;+ for (i = 0;i < bound;++i) {+ j = (j + m - 1) % m; /* if (--j == -1) j = m-1; */+ it = &lm[j];+ /* \alpha_{j} = \rho_{j} s^{t}_{j} \cdot q_{k+1}. */+ vecdot(&it->alpha, it->s, d, n);+ it->alpha /= it->ys;+ /* q_{i} = q_{i+1} - \alpha_{i} y_{i}. */+ vecadd(d, it->y, -it->alpha, n);+ }++ vecscale(d, ys / yy, n);++ for (i = 0;i < bound;++i) {+ it = &lm[j];+ /* \beta_{j} = \rho_{j} y^t_{j} \cdot \gamma_{i}. */+ vecdot(&beta, it->y, d, n);+ beta /= it->ys;+ /* \gamma_{i+1} = \gamma_{i} + (\alpha_{j} - \beta_{j}) s_{j}. */+ vecadd(d, it->s, it->alpha - beta, n);+ j = (j + 1) % m; /* if (++j == m) j = 0; */+ }++ /*+ Constrain the search direction for orthant-wise updates.+ */+ if (param.orthantwise_c != 0.) {+ for (i = param.orthantwise_start;i < param.orthantwise_end;++i) {+ if (d[i] * pg[i] >= 0) {+ d[i] = 0;+ }+ }+ }++ /*+ Now the search direction d is ready. We try step = 1 first.+ */+ step = 1.0;+ }++lbfgs_exit:+ /* Return the final value of the objective function. */+ if (ptr_fx != NULL) {+ *ptr_fx = fx;+ }++ vecfree(pf);++ /* Free memory blocks used by this function. */+ if (lm != NULL) {+ for (i = 0;i < m;++i) {+ vecfree(lm[i].s);+ vecfree(lm[i].y);+ }+ vecfree(lm);+ }+ vecfree(pg);+ vecfree(w);+ vecfree(d);+ vecfree(gp);+ vecfree(g);+ vecfree(xp);++ return ret;+}++++static int line_search_backtracking(+ int n,+ lbfgsfloatval_t *x,+ lbfgsfloatval_t *f,+ lbfgsfloatval_t *g,+ lbfgsfloatval_t *s,+ lbfgsfloatval_t *stp,+ const lbfgsfloatval_t* xp,+ const lbfgsfloatval_t* gp,+ lbfgsfloatval_t *wp,+ callback_data_t *cd,+ const lbfgs_parameter_t *param+ )+{+ int ret = 0, count = 0;+ lbfgsfloatval_t width, dg, norm = 0.;+ lbfgsfloatval_t finit, dginit = 0., dgtest;+ const lbfgsfloatval_t dec = 0.5, inc = 2.1;++ /* Check the input parameters for errors. */+ if (*stp <= 0.) {+ return LBFGSERR_INVALIDPARAMETERS;+ }++ /* Compute the initial gradient in the search direction. */+ vecdot(&dginit, g, s, n);++ /* Make sure that s points to a descent direction. */+ if (0 < dginit) {+ return LBFGSERR_INCREASEGRADIENT;+ }++ /* The initial value of the objective function. */+ finit = *f;+ dgtest = param->ftol * dginit;++ for (;;) {+ veccpy(x, xp, n);+ vecadd(x, s, *stp, n);++ /* Evaluate the function and gradient values. */+ *f = cd->proc_evaluate(cd->instance, x, g, cd->n, *stp);++ ++count;++ if (*f > finit + *stp * dgtest) {+ width = dec;+ } else {+ /* The sufficient decrease condition (Armijo condition). */+ if (param->linesearch == LBFGS_LINESEARCH_BACKTRACKING_ARMIJO) {+ /* Exit with the Armijo condition. */+ return count;+ }++ /* Check the Wolfe condition. */+ vecdot(&dg, g, s, n);+ if (dg < param->wolfe * dginit) {+ width = inc;+ } else {+ if(param->linesearch == LBFGS_LINESEARCH_BACKTRACKING_WOLFE) {+ /* Exit with the regular Wolfe condition. */+ return count;+ }++ /* Check the strong Wolfe condition. */+ if(dg > -param->wolfe * dginit) {+ width = dec;+ } else {+ /* Exit with the strong Wolfe condition. */+ return count;+ }+ }+ }++ if (*stp < param->min_step) {+ /* The step is the minimum value. */+ return LBFGSERR_MINIMUMSTEP;+ }+ if (*stp > param->max_step) {+ /* The step is the maximum value. */+ return LBFGSERR_MAXIMUMSTEP;+ }+ if (param->max_linesearch <= count) {+ /* Maximum number of iteration. */+ return LBFGSERR_MAXIMUMLINESEARCH;+ }++ (*stp) *= width;+ }+}++++static int line_search_backtracking_owlqn(+ int n,+ lbfgsfloatval_t *x,+ lbfgsfloatval_t *f,+ lbfgsfloatval_t *g,+ lbfgsfloatval_t *s,+ lbfgsfloatval_t *stp,+ const lbfgsfloatval_t* xp,+ const lbfgsfloatval_t* gp,+ lbfgsfloatval_t *wp,+ callback_data_t *cd,+ const lbfgs_parameter_t *param+ )+{+ int i, ret = 0, count = 0;+ lbfgsfloatval_t width = 0.5, norm = 0.;+ lbfgsfloatval_t finit = *f, dgtest;++ /* Check the input parameters for errors. */+ if (*stp <= 0.) {+ return LBFGSERR_INVALIDPARAMETERS;+ }++ /* Choose the orthant for the new point. */+ for (i = 0;i < n;++i) {+ wp[i] = (xp[i] == 0.) ? -gp[i] : xp[i];+ }++ for (;;) {+ /* Update the current point. */+ veccpy(x, xp, n);+ vecadd(x, s, *stp, n);++ /* The current point is projected onto the orthant. */+ owlqn_project(x, wp, param->orthantwise_start, param->orthantwise_end);++ /* Evaluate the function and gradient values. */+ *f = cd->proc_evaluate(cd->instance, x, g, cd->n, *stp);++ /* Compute the L1 norm of the variables and add it to the object value. */+ norm = owlqn_x1norm(x, param->orthantwise_start, param->orthantwise_end);+ *f += norm * param->orthantwise_c;++ ++count;++ dgtest = 0.;+ for (i = 0;i < n;++i) {+ dgtest += (x[i] - xp[i]) * gp[i];+ }++ if (*f <= finit + param->ftol * dgtest) {+ /* The sufficient decrease condition. */+ return count;+ }++ if (*stp < param->min_step) {+ /* The step is the minimum value. */+ return LBFGSERR_MINIMUMSTEP;+ }+ if (*stp > param->max_step) {+ /* The step is the maximum value. */+ return LBFGSERR_MAXIMUMSTEP;+ }+ if (param->max_linesearch <= count) {+ /* Maximum number of iteration. */+ return LBFGSERR_MAXIMUMLINESEARCH;+ }++ (*stp) *= width;+ }+}++++static int line_search_morethuente(+ int n,+ lbfgsfloatval_t *x,+ lbfgsfloatval_t *f,+ lbfgsfloatval_t *g,+ lbfgsfloatval_t *s,+ lbfgsfloatval_t *stp,+ const lbfgsfloatval_t* xp,+ const lbfgsfloatval_t* gp,+ lbfgsfloatval_t *wa,+ callback_data_t *cd,+ const lbfgs_parameter_t *param+ )+{+ int count = 0;+ int brackt, stage1, uinfo = 0;+ lbfgsfloatval_t dg;+ lbfgsfloatval_t stx, fx, dgx;+ lbfgsfloatval_t sty, fy, dgy;+ lbfgsfloatval_t fxm, dgxm, fym, dgym, fm, dgm;+ lbfgsfloatval_t finit, ftest1, dginit, dgtest;+ lbfgsfloatval_t width, prev_width;+ lbfgsfloatval_t stmin, stmax;++ /* Check the input parameters for errors. */+ if (*stp <= 0.) {+ return LBFGSERR_INVALIDPARAMETERS;+ }++ /* Compute the initial gradient in the search direction. */+ vecdot(&dginit, g, s, n);++ /* Make sure that s points to a descent direction. */+ if (0 < dginit) {+ return LBFGSERR_INCREASEGRADIENT;+ }++ /* Initialize local variables. */+ brackt = 0;+ stage1 = 1;+ finit = *f;+ dgtest = param->ftol * dginit;+ width = param->max_step - param->min_step;+ prev_width = 2.0 * width;++ /*+ The variables stx, fx, dgx contain the values of the step,+ function, and directional derivative at the best step.+ The variables sty, fy, dgy contain the value of the step,+ function, and derivative at the other endpoint of+ the interval of uncertainty.+ The variables stp, f, dg contain the values of the step,+ function, and derivative at the current step.+ */+ stx = sty = 0.;+ fx = fy = finit;+ dgx = dgy = dginit;++ for (;;) {+ /*+ Set the minimum and maximum steps to correspond to the+ present interval of uncertainty.+ */+ if (brackt) {+ stmin = min2(stx, sty);+ stmax = max2(stx, sty);+ } else {+ stmin = stx;+ stmax = *stp + 4.0 * (*stp - stx);+ }++ /* Clip the step in the range of [stpmin, stpmax]. */+ if (*stp < param->min_step) *stp = param->min_step;+ if (param->max_step < *stp) *stp = param->max_step;++ /*+ If an unusual termination is to occur then let+ stp be the lowest point obtained so far.+ */+ if ((brackt && ((*stp <= stmin || stmax <= *stp) || param->max_linesearch <= count + 1 || uinfo != 0)) || (brackt && (stmax - stmin <= param->xtol * stmax))) {+ *stp = stx;+ }++ /*+ Compute the current value of x:+ x <- x + (*stp) * s.+ */+ veccpy(x, xp, n);+ vecadd(x, s, *stp, n);++ /* Evaluate the function and gradient values. */+ *f = cd->proc_evaluate(cd->instance, x, g, cd->n, *stp);+ vecdot(&dg, g, s, n);++ ftest1 = finit + *stp * dgtest;+ ++count;++ /* Test for errors and convergence. */+ if (brackt && ((*stp <= stmin || stmax <= *stp) || uinfo != 0)) {+ /* Rounding errors prevent further progress. */+ return LBFGSERR_ROUNDING_ERROR;+ }+ if (*stp == param->max_step && *f <= ftest1 && dg <= dgtest) {+ /* The step is the maximum value. */+ return LBFGSERR_MAXIMUMSTEP;+ }+ if (*stp == param->min_step && (ftest1 < *f || dgtest <= dg)) {+ /* The step is the minimum value. */+ return LBFGSERR_MINIMUMSTEP;+ }+ if (brackt && (stmax - stmin) <= param->xtol * stmax) {+ /* Relative width of the interval of uncertainty is at most xtol. */+ return LBFGSERR_WIDTHTOOSMALL;+ }+ if (param->max_linesearch <= count) {+ /* Maximum number of iteration. */+ return LBFGSERR_MAXIMUMLINESEARCH;+ }+ if (*f <= ftest1 && fabs(dg) <= param->gtol * (-dginit)) {+ /* The sufficient decrease condition and the directional derivative condition hold. */+ return count;+ }++ /*+ In the first stage we seek a step for which the modified+ function has a nonpositive value and nonnegative derivative.+ */+ if (stage1 && *f <= ftest1 && min2(param->ftol, param->gtol) * dginit <= dg) {+ stage1 = 0;+ }++ /*+ A modified function is used to predict the step only if+ we have not obtained a step for which the modified+ function has a nonpositive function value and nonnegative+ derivative, and if a lower function value has been+ obtained but the decrease is not sufficient.+ */+ if (stage1 && ftest1 < *f && *f <= fx) {+ /* Define the modified function and derivative values. */+ fm = *f - *stp * dgtest;+ fxm = fx - stx * dgtest;+ fym = fy - sty * dgtest;+ dgm = dg - dgtest;+ dgxm = dgx - dgtest;+ dgym = dgy - dgtest;++ /*+ Call update_trial_interval() to update the interval of+ uncertainty and to compute the new step.+ */+ uinfo = update_trial_interval(+ &stx, &fxm, &dgxm,+ &sty, &fym, &dgym,+ stp, &fm, &dgm,+ stmin, stmax, &brackt+ );++ /* Reset the function and gradient values for f. */+ fx = fxm + stx * dgtest;+ fy = fym + sty * dgtest;+ dgx = dgxm + dgtest;+ dgy = dgym + dgtest;+ } else {+ /*+ Call update_trial_interval() to update the interval of+ uncertainty and to compute the new step.+ */+ uinfo = update_trial_interval(+ &stx, &fx, &dgx,+ &sty, &fy, &dgy,+ stp, f, &dg,+ stmin, stmax, &brackt+ );+ }++ /*+ Force a sufficient decrease in the interval of uncertainty.+ */+ if (brackt) {+ if (0.66 * prev_width <= fabs(sty - stx)) {+ *stp = stx + 0.5 * (sty - stx);+ }+ prev_width = width;+ width = fabs(sty - stx);+ }+ }++ return LBFGSERR_LOGICERROR;+}++++/**+ * Define the local variables for computing minimizers.+ */+#define USES_MINIMIZER \+ lbfgsfloatval_t a, d, gamma, theta, p, q, r, s;++/**+ * Find a minimizer of an interpolated cubic function.+ * @param cm The minimizer of the interpolated cubic.+ * @param u The value of one point, u.+ * @param fu The value of f(u).+ * @param du The value of f'(u).+ * @param v The value of another point, v.+ * @param fv The value of f(v).+ * @param du The value of f'(v).+ */+#define CUBIC_MINIMIZER(cm, u, fu, du, v, fv, dv) \+ d = (v) - (u); \+ theta = ((fu) - (fv)) * 3 / d + (du) + (dv); \+ p = fabs(theta); \+ q = fabs(du); \+ r = fabs(dv); \+ s = max3(p, q, r); \+ /* gamma = s*sqrt((theta/s)**2 - (du/s) * (dv/s)) */ \+ a = theta / s; \+ gamma = s * sqrt(a * a - ((du) / s) * ((dv) / s)); \+ if ((v) < (u)) gamma = -gamma; \+ p = gamma - (du) + theta; \+ q = gamma - (du) + gamma + (dv); \+ r = p / q; \+ (cm) = (u) + r * d;++/**+ * Find a minimizer of an interpolated cubic function.+ * @param cm The minimizer of the interpolated cubic.+ * @param u The value of one point, u.+ * @param fu The value of f(u).+ * @param du The value of f'(u).+ * @param v The value of another point, v.+ * @param fv The value of f(v).+ * @param du The value of f'(v).+ * @param xmin The maximum value.+ * @param xmin The minimum value.+ */+#define CUBIC_MINIMIZER2(cm, u, fu, du, v, fv, dv, xmin, xmax) \+ d = (v) - (u); \+ theta = ((fu) - (fv)) * 3 / d + (du) + (dv); \+ p = fabs(theta); \+ q = fabs(du); \+ r = fabs(dv); \+ s = max3(p, q, r); \+ /* gamma = s*sqrt((theta/s)**2 - (du/s) * (dv/s)) */ \+ a = theta / s; \+ gamma = s * sqrt(max2(0, a * a - ((du) / s) * ((dv) / s))); \+ if ((u) < (v)) gamma = -gamma; \+ p = gamma - (dv) + theta; \+ q = gamma - (dv) + gamma + (du); \+ r = p / q; \+ if (r < 0. && gamma != 0.) { \+ (cm) = (v) - r * d; \+ } else if (a < 0) { \+ (cm) = (xmax); \+ } else { \+ (cm) = (xmin); \+ }++/**+ * Find a minimizer of an interpolated quadratic function.+ * @param qm The minimizer of the interpolated quadratic.+ * @param u The value of one point, u.+ * @param fu The value of f(u).+ * @param du The value of f'(u).+ * @param v The value of another point, v.+ * @param fv The value of f(v).+ */+#define QUARD_MINIMIZER(qm, u, fu, du, v, fv) \+ a = (v) - (u); \+ (qm) = (u) + (du) / (((fu) - (fv)) / a + (du)) / 2 * a;++/**+ * Find a minimizer of an interpolated quadratic function.+ * @param qm The minimizer of the interpolated quadratic.+ * @param u The value of one point, u.+ * @param du The value of f'(u).+ * @param v The value of another point, v.+ * @param dv The value of f'(v).+ */+#define QUARD_MINIMIZER2(qm, u, du, v, dv) \+ a = (u) - (v); \+ (qm) = (v) + (dv) / ((dv) - (du)) * a;++/**+ * Update a safeguarded trial value and interval for line search.+ *+ * The parameter x represents the step with the least function value.+ * The parameter t represents the current step. This function assumes+ * that the derivative at the point of x in the direction of the step.+ * If the bracket is set to true, the minimizer has been bracketed in+ * an interval of uncertainty with endpoints between x and y.+ *+ * @param x The pointer to the value of one endpoint.+ * @param fx The pointer to the value of f(x).+ * @param dx The pointer to the value of f'(x).+ * @param y The pointer to the value of another endpoint.+ * @param fy The pointer to the value of f(y).+ * @param dy The pointer to the value of f'(y).+ * @param t The pointer to the value of the trial value, t.+ * @param ft The pointer to the value of f(t).+ * @param dt The pointer to the value of f'(t).+ * @param tmin The minimum value for the trial value, t.+ * @param tmax The maximum value for the trial value, t.+ * @param brackt The pointer to the predicate if the trial value is+ * bracketed.+ * @retval int Status value. Zero indicates a normal termination.+ * + * @see+ * Jorge J. More and David J. Thuente. Line search algorithm with+ * guaranteed sufficient decrease. ACM Transactions on Mathematical+ * Software (TOMS), Vol 20, No 3, pp. 286-307, 1994.+ */+static int update_trial_interval(+ lbfgsfloatval_t *x,+ lbfgsfloatval_t *fx,+ lbfgsfloatval_t *dx,+ lbfgsfloatval_t *y,+ lbfgsfloatval_t *fy,+ lbfgsfloatval_t *dy,+ lbfgsfloatval_t *t,+ lbfgsfloatval_t *ft,+ lbfgsfloatval_t *dt,+ const lbfgsfloatval_t tmin,+ const lbfgsfloatval_t tmax,+ int *brackt+ )+{+ int bound;+ int dsign = fsigndiff(dt, dx);+ lbfgsfloatval_t mc; /* minimizer of an interpolated cubic. */+ lbfgsfloatval_t mq; /* minimizer of an interpolated quadratic. */+ lbfgsfloatval_t newt; /* new trial value. */+ USES_MINIMIZER; /* for CUBIC_MINIMIZER and QUARD_MINIMIZER. */++ /* Check the input parameters for errors. */+ if (*brackt) {+ if (*t <= min2(*x, *y) || max2(*x, *y) <= *t) {+ /* The trival value t is out of the interval. */+ return LBFGSERR_OUTOFINTERVAL;+ }+ if (0. <= *dx * (*t - *x)) {+ /* The function must decrease from x. */+ return LBFGSERR_INCREASEGRADIENT;+ }+ if (tmax < tmin) {+ /* Incorrect tmin and tmax specified. */+ return LBFGSERR_INCORRECT_TMINMAX;+ }+ }++ /*+ Trial value selection.+ */+ if (*fx < *ft) {+ /*+ Case 1: a higher function value.+ The minimum is brackt. If the cubic minimizer is closer+ to x than the quadratic one, the cubic one is taken, else+ the average of the minimizers is taken.+ */+ *brackt = 1;+ bound = 1;+ CUBIC_MINIMIZER(mc, *x, *fx, *dx, *t, *ft, *dt);+ QUARD_MINIMIZER(mq, *x, *fx, *dx, *t, *ft);+ if (fabs(mc - *x) < fabs(mq - *x)) {+ newt = mc;+ } else {+ newt = mc + 0.5 * (mq - mc);+ }+ } else if (dsign) {+ /*+ Case 2: a lower function value and derivatives of+ opposite sign. The minimum is brackt. If the cubic+ minimizer is closer to x than the quadratic (secant) one,+ the cubic one is taken, else the quadratic one is taken.+ */+ *brackt = 1;+ bound = 0;+ CUBIC_MINIMIZER(mc, *x, *fx, *dx, *t, *ft, *dt);+ QUARD_MINIMIZER2(mq, *x, *dx, *t, *dt);+ if (fabs(mc - *t) > fabs(mq - *t)) {+ newt = mc;+ } else {+ newt = mq;+ }+ } else if (fabs(*dt) < fabs(*dx)) {+ /*+ Case 3: a lower function value, derivatives of the+ same sign, and the magnitude of the derivative decreases.+ The cubic minimizer is only used if the cubic tends to+ infinity in the direction of the minimizer or if the minimum+ of the cubic is beyond t. Otherwise the cubic minimizer is+ defined to be either tmin or tmax. The quadratic (secant)+ minimizer is also computed and if the minimum is brackt+ then the the minimizer closest to x is taken, else the one+ farthest away is taken.+ */+ bound = 1;+ CUBIC_MINIMIZER2(mc, *x, *fx, *dx, *t, *ft, *dt, tmin, tmax);+ QUARD_MINIMIZER2(mq, *x, *dx, *t, *dt);+ if (*brackt) {+ if (fabs(*t - mc) < fabs(*t - mq)) {+ newt = mc;+ } else {+ newt = mq;+ }+ } else {+ if (fabs(*t - mc) > fabs(*t - mq)) {+ newt = mc;+ } else {+ newt = mq;+ }+ }+ } else {+ /*+ Case 4: a lower function value, derivatives of the+ same sign, and the magnitude of the derivative does+ not decrease. If the minimum is not brackt, the step+ is either tmin or tmax, else the cubic minimizer is taken.+ */+ bound = 0;+ if (*brackt) {+ CUBIC_MINIMIZER(newt, *t, *ft, *dt, *y, *fy, *dy);+ } else if (*x < *t) {+ newt = tmax;+ } else {+ newt = tmin;+ }+ }++ /*+ Update the interval of uncertainty. This update does not+ depend on the new step or the case analysis above.++ - Case a: if f(x) < f(t),+ x <- x, y <- t.+ - Case b: if f(t) <= f(x) && f'(t)*f'(x) > 0,+ x <- t, y <- y.+ - Case c: if f(t) <= f(x) && f'(t)*f'(x) < 0, + x <- t, y <- x.+ */+ if (*fx < *ft) {+ /* Case a */+ *y = *t;+ *fy = *ft;+ *dy = *dt;+ } else {+ /* Case c */+ if (dsign) {+ *y = *x;+ *fy = *fx;+ *dy = *dx;+ }+ /* Cases b and c */+ *x = *t;+ *fx = *ft;+ *dx = *dt;+ }++ /* Clip the new trial value in [tmin, tmax]. */+ if (tmax < newt) newt = tmax;+ if (newt < tmin) newt = tmin;++ /*+ Redefine the new trial value if it is close to the upper bound+ of the interval.+ */+ if (*brackt && bound) {+ mq = *x + 0.66 * (*y - *x);+ if (*x < *y) {+ if (mq < newt) newt = mq;+ } else {+ if (newt < mq) newt = mq;+ }+ }++ /* Return the new trial value. */+ *t = newt;+ return 0;+}++++++static lbfgsfloatval_t owlqn_x1norm(+ const lbfgsfloatval_t* x,+ const int start,+ const int n+ )+{+ int i;+ lbfgsfloatval_t norm = 0.;++ for (i = start;i < n;++i) {+ norm += fabs(x[i]);+ }++ return norm;+}++static void owlqn_pseudo_gradient(+ lbfgsfloatval_t* pg,+ const lbfgsfloatval_t* x,+ const lbfgsfloatval_t* g,+ const int n,+ const lbfgsfloatval_t c,+ const int start,+ const int end+ )+{+ int i;++ /* Compute the negative of gradients. */+ for (i = 0;i < start;++i) {+ pg[i] = g[i];+ }++ /* Compute the psuedo-gradients. */+ for (i = start;i < end;++i) {+ if (x[i] < 0.) {+ /* Differentiable. */+ pg[i] = g[i] - c;+ } else if (0. < x[i]) {+ /* Differentiable. */+ pg[i] = g[i] + c;+ } else {+ if (g[i] < -c) {+ /* Take the right partial derivative. */+ pg[i] = g[i] + c;+ } else if (c < g[i]) {+ /* Take the left partial derivative. */+ pg[i] = g[i] - c;+ } else {+ pg[i] = 0.;+ }+ }+ }++ for (i = end;i < n;++i) {+ pg[i] = g[i];+ }+}++static void owlqn_project(+ lbfgsfloatval_t* d,+ const lbfgsfloatval_t* sign,+ const int start,+ const int end+ )+{+ int i;++ for (i = start;i < end;++i) {+ if (d[i] * sign[i] <= 0) {+ d[i] = 0;+ }+ }+}
+ lbfgs.cabal view
@@ -0,0 +1,19 @@+Name: lbfgs+Version: 0.0.1+License: OtherLicense+License-File: LICENSE+Copyright: Daniël de Kok+Maintainer: Daniël de Kok <me@danieldk.eu>+Author: Daniël de Kok <me@danieldk.eu>+Category: Numeric+Synopsis: L-BFGS optimization+Description: Limited memory BFGS solver for non-linear optimization+ problems.+Build-Type: Simple+Cabal-Version: >= 1.4++Library+ Build-Depends: base >= 4 && < 5, array >= 0.3.0.0+ Exposed-modules: Numeric.LBFGS.Raw, Numeric.LBFGS+ Include-Dirs: cbits+ C-Sources: cbits/lbfgs.c