diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -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
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -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
+```
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/examples/rosenbrock.hs b/examples/rosenbrock.hs
new file mode 100644
--- /dev/null
+++ b/examples/rosenbrock.hs
@@ -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
diff --git a/numeric-optimization-ad.cabal b/numeric-optimization-ad.cabal
new file mode 100644
--- /dev/null
+++ b/numeric-optimization-ad.cabal
@@ -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
diff --git a/src/Numeric/Optimization/AD.hs b/src/Numeric/Optimization/AD.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Optimization/AD.hs
@@ -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 ()
+
+-- ------------------------------------------------------------------------
diff --git a/test/IsClose.hs b/test/IsClose.hs
new file mode 100644
--- /dev/null
+++ b/test/IsClose.hs
@@ -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]
+
+-- ------------------------------------------------------------------------
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -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