diff --git a/Ipopt.chs b/Ipopt.chs
deleted file mode 100644
--- a/Ipopt.chs
+++ /dev/null
@@ -1,359 +0,0 @@
-{-#  LANGUAGE ForeignFunctionInterface, PatternGuards, RankNTypes, TypeFamilies #-}
-{- |
-
-Copyright: (C) 2013 Adam Vogt
-Maintainer: Adam Vogt <vogt.adam@gmail.com>
-Stability: unstable
-
-Binding to ipopt <http://projects.coin-or.org/Ipopt>. Uses "Numeric.AD" to compute
-derivatives required by ipopt. Current limitations include:
-
-* derivatives are computed and stored without taking advantage of sparsity
-
-* copying is done in converting between "Data.Vector.Storable" and "Data.Vector" might be unnecessary. Currently it is done because AD needs a Traversable structure, but Storable vectors are not traversable.
-
-* probably doesn't work if @coin\/IpStdCInterface.h@ has Number =/= 'CDouble'
-
-* no binding to SetIntermediateCallback
-
-* garbage collection of 'IpProblem' won't free C-side resources
-
-* specifying problems might be made easier by following (extending) the approach taken by glpk-hs
-
-Refer to @Test1.hs@ for an example where the derivatives are computed by hand,
-and @Test2.hs@ for the use of 'createIpoptProblemAD'.
-
--}
-module Ipopt (
-    -- * specifying problem
-    createIpoptProblemAD,
-
-    -- ** solve
-    ipoptSolve,
-    IpOptSolved(..),
-
-    -- ** solver options
-    addIpoptNumOption,
-    addIpoptStrOption,
-    addIpoptIntOption,
-    openIpoptOutputFile,
-
-    -- * types
-    Vec,
-
-    IpNumber(..),
-    IpIndex(..),
-    IpInt(..),
-    IpBool(..),
-
-    IpF(..),
-    IpGradF(..),
-    IpG(..),
-    IpJacG(..),
-    IpH(..),
-
-    IpProblem(..),
-
-    ApplicationReturnStatus(..),
-
-    -- * lower-level parts of the binding
-    createIpoptProblem,
-
-    freeIpoptProblem,
-    setIpoptProblemScaling,
-
-
-    -- ** marshalling functions
-    wrapIpF,
-    wrapIpGradF,
-    wrapIpG,
-    wrapIpJacG,
-    wrapIpH,
-
-    wrapIpF1,
-    wrapIpGradF1,
-    wrapIpG1,
-    wrapIpJacG1,
-    wrapIpH1,
-
-    wrapIpF2,
-    wrapIpGradF2,
-    wrapIpG2,
-    wrapIpJacG2,
-    wrapIpH2,
-
-    ) where
-
-import C2HS
-import Control.Exception
-import Control.Monad
-import Data.IORef
-import Foreign.C
-import Foreign.ForeignPtr
-import Foreign.Ptr
-import Foreign.Storable
-import Numeric.AD
-import qualified Data.Vector as V
-import qualified Data.Vector.Generic as VG
-import qualified Data.Vector.Storable as VS
-import qualified Data.Vector.Storable.Mutable as VM
-
-
-
-#include "coin/IpStdCInterface.h"
-
-type IpNumber = {# type Number #}
-type IpIndex = {# type Index #}
-type IpInt = {# type Int #}
-type IpBool = {# type Bool #}
-
-type IpF = {# type Eval_F_CB #}
-type IpGradF = {# type Eval_Grad_F_CB #}
-type IpG = {# type Eval_G_CB #}
-type IpJacG = {# type Eval_Jac_G_CB #}
-type IpH = {# type Eval_H_CB #}
-
-{#enum ApplicationReturnStatus as ^ {underscoreToCase} deriving (Show) #}
-{#enum AlgorithmMode as ^ {underscoreToCase} deriving (Show) #}
-
-
-newtype IpProblem = IpProblem { unIpProblem :: Ptr ()}
-
-type family UnFunPtr a
-type instance UnFunPtr (FunPtr a) = a
-
-ipTrue = 1 :: IpBool
-ipFalse = 0 :: IpBool
-
--- | likely an unsafe method for getting a "Data.Vector.Storable.Mutable" out of a 'Ptr'
-ptrToVS n p = do
-    fp <- newForeignPtr_ p
-    return (VM.unsafeFromForeignPtr0 fp (fromIntegral n))
-
-foreign import ccall "wrapper" wrapIpF1 :: UnFunPtr IpF -> IO IpF
-foreign import ccall "wrapper" wrapIpG1 :: UnFunPtr IpG -> IO IpG
-foreign import ccall "wrapper" wrapIpGradF1 :: UnFunPtr IpGradF -> IO IpGradF
-foreign import ccall "wrapper" wrapIpJacG1 :: UnFunPtr IpJacG -> IO IpJacG
-foreign import ccall "wrapper" wrapIpH1 :: UnFunPtr IpH -> IO IpH
-
-
-toB x =  either (\ e@SomeException {} -> print e >> return ipFalse)
-                (\ _ -> return ipTrue ) =<< try x
-
-wrapIpF2' fun n xin new_x obj_val _userData = do
-    toB $ poke obj_val =<< fun =<< ptrToVS n xin
-
-wrapIpF2 fun n xin new_x obj_val _userData = do
-    toB $ poke obj_val =<< fun =<< ptrToVS n xin
-
-wrapIpG2 fun n xin new_x m gout _userData = do
-    toB $ join $ liftM2 VM.copy (ptrToVS m gout) (fun =<< ptrToVS n xin)
-
-wrapIpGradF2 fun n x new_x grad_f _userData = do
-    toB $ join $ liftM2 VM.copy (ptrToVS n grad_f) (fun =<< ptrToVS n x)
-
-wrapIpJacG2 fun1 fun2 n x new_x m nj iRow jCol jacs _userData
-    | jacs == nullPtr = do
-            toB $ join $ liftM2 fun1 (ptrToVS nj iRow) (ptrToVS nj jCol)
-    | otherwise = do
-            toB $ join $ liftM2 fun2 (ptrToVS n x) (ptrToVS nj jacs)
-
-wrapIpH2 funSparsity funEval n x new_x obj_factor m lambda new_lambda nHess iRow jCol values _userData
-    | iRow == nullPtr = do
-            toB $ join $ liftM3 (funEval obj_factor)
-                        (ptrToVS m lambda)
-                        (ptrToVS n x)
-                        (ptrToVS nHess values)
-    | otherwise = do
-            toB $ join $ liftM2 funSparsity
-                        (ptrToVS nHess iRow)
-                        (ptrToVS nHess jCol)
-
-wrapIpF f = wrapIpF1 (wrapIpF2 f)
-wrapIpG f = wrapIpG1 (wrapIpG2 f)
-wrapIpGradF f = wrapIpGradF1 (wrapIpGradF2 f)
-wrapIpJacG f1 f2 = wrapIpJacG1 (wrapIpJacG2 f1 f2)
-wrapIpH fSparsity fEval = wrapIpH1 (wrapIpH2 fSparsity fEval)
-
-
-vmUnsafeWith = VM.unsafeWith
-
--- | Vector of numbers
-type Vec = VM.IOVector IpNumber
-
-createIpoptProblem :: Vec -> Vec -> Vec -> Vec
-    -> Int -> Int -> IpF -> IpG -> IpGradF -> IpJacG -> IpH -> IO IpProblem
-createIpoptProblem xL xU gL gU nJac nHess f g gradF jacG hess
-    | lx <- VM.length xL,
-      lx == VM.length xU,
-      lg <- VM.length gL,
-      lg == VM.length gU = createIpoptProblem3 lx xL xU lg gL gU nJac nHess 0 f g gradF jacG hess
-   | otherwise = error "dimensions wrong!"
-
-{#fun CreateIpoptProblem as createIpoptProblem3
-        { `Int', vmUnsafeWith* `Vec', vmUnsafeWith* `Vec',
-          `Int', vmUnsafeWith* `Vec', vmUnsafeWith* `Vec',
-          `Int', `Int', `Int', id `IpF', id `IpG', id `IpGradF',
-          id `IpJacG', id `IpH' } -> `IpProblem' IpProblem #}
-_ = {#fun AddIpoptNumOption as ^
-    { unIpProblem `IpProblem', `String', `Double' } -> `Bool' #}
-
-_ = {#fun AddIpoptStrOption as ^
-    { unIpProblem `IpProblem', `String', `String' } -> `Bool' #}
-
-_ = {#fun AddIpoptIntOption as ^
-    { unIpProblem `IpProblem', `String', `Int' } -> `Bool' #}
-
-_ = {#fun FreeIpoptProblem as ^
-    { unIpProblem `IpProblem' } -> `()' #}
-
-_ = {#fun OpenIpoptOutputFile as ^
-    { unIpProblem `IpProblem', `String', `Int' } -> `Bool' #}
-
-_ = {#fun SetIpoptProblemScaling as ^
-    { unIpProblem `IpProblem',
-    `Double',
-     vmUnsafeWith* `Vec',
-     vmUnsafeWith* `Vec'
-     } -> `Bool' #}
-
-
-data IpOptSolved = IpOptSolved
-  { ipOptSolved_status :: ApplicationReturnStatus,
-    ipOptSolved_objective :: Double,
-    ipOptSolved_x,
-    ipOptSolved_g,
-    ipOptSolved_mult_g,
-    ipOptSolved_mult_x_L,
-    ipOptSolved_mult_x_U :: Vec }
-
-ipoptSolve :: IpProblem
-    -> Vec -- ^ starting point @x@. Note that the value is overwritte with the final @x@.
-    -> IO IpOptSolved
-ipoptSolve problem x = do
-    g <- VM.new (VM.length x)
-    mult_g <- VM.new (VM.length x)
-    mult_x_L <- VM.new (VM.length x)
-    mult_x_U <- VM.new (VM.length x)
-
-    out <- ipoptSolve2
-     problem
-     x
-     g
-     mult_g
-     mult_x_L
-     mult_x_U
-     nullPtr
-    return $ IpOptSolved
-      (fst out)
-      (snd out)
-      x
-      g
-      mult_g
-      mult_x_L
-      mult_x_U
-
-_ = {#fun IpoptSolve as ipoptSolve2
-        { unIpProblem `IpProblem',
-        vmUnsafeWith* `Vec',
-        vmUnsafeWith* `Vec',
-        alloca- `Double' peekFloatConv*,
-        vmUnsafeWith* `Vec',
-        vmUnsafeWith* `Vec',
-        vmUnsafeWith* `Vec',
-        id `Ptr ()' } -> `ApplicationReturnStatus' cToEnum #}
-
-
-{- | Set-up an 'IpProblem' to be solved later. Only objective function (@f@)
-and constraint functions (@g@) need to be specified. Derivatives needed by ipopt
-are computed by "Numeric.AD".
-
-To solve the optimization problem:
-
->              min f(x)
->   such that
->           xL <=  x     <= xU
->           gL <=  g(x)  <= gU
-
-First create an opaque 'IpProblem' object (nlp):
-
-> nlp <- createIpOptProblemAD xL xU gL gU f g
-
-Then pass it off to 'ipoptSolve'.
-
-> ipoptSolve nlp x0
-
-Refer to the example @Test2.hs@ for details of setting up the vectors supplied.
--}
-createIpoptProblemAD
-    :: Vec -- ^ @xL@ 'VM.Vector' of lower bounds for decision variables with length @n@
-    -> Vec -- ^ @xU@ 'VM.Vector' of upper bounds for decision variables
-    -> Vec -- ^ @gL@ 'VM.Vector' of lower bounds for constraint functions @g(x)@ with length @m@
-    -> Vec -- ^ @gU@ 'VM.Vector' of upper bounds for constraint functions @g(x)@
-    -> (forall a. Num a => V.Vector a -> a) -- ^ objective function @f : R^n -> R@
-    -> (forall a. Num a => V.Vector a -> V.Vector a) -- ^ constraint functions @g : R^n -> R^m@
-    -> IO IpProblem
-createIpoptProblemAD xL xU gL gU f g
-    | n <- VM.length xL,
-      n == VM.length xU,
-      m <- VM.length gL,
-      m == VM.length gU = do
-    (eval_f, eval_grad_f, eval_g, eval_jac_g, eval_h) <- mkFs n m f g
-    createIpoptProblem xL xU gL gU (n*m) (((n+1)*n) `div` 2)
-            eval_f eval_g eval_grad_f eval_jac_g eval_h
-
-mkFs :: Int -- ^ @n@ number of variables
-    -> Int -- ^ @m@ number of constraints
-    -> (forall a. Num a => V.Vector a -> a) -- ^ objective function @R^n -> R@
-    -> (forall a. Num a => V.Vector a -> V.Vector a) -- ^ constraint functions @R^n -> R^m@
-    -> IO (IpF, IpGradF, IpG, IpJacG, IpH)
-mkFs n m f g = do
-  ipF <- wrapIpF $ \x -> do
-    x <- VG.convert `fmap` VS.unsafeFreeze x
-    return $ f x
-
-  ipGradF <- wrapIpGradF $ \x -> do
-    x <- VG.convert `fmap` VS.unsafeFreeze x
-    VS.unsafeThaw $ VG.convert (grad f x)
-
-  ipG <- wrapIpG $ \x -> do
-    x <- VG.convert `fmap` VS.unsafeFreeze x
-    VS.unsafeThaw $ VG.convert (g x)
-
-  ipJacG <- wrapIpJacG (denseIJ n m) $ \x y -> do
-    x <- VG.convert `fmap` VS.unsafeFreeze x
-    jac <- VS.unsafeThaw $ VG.convert $ VG.concat $ VG.toList $ jacobian g x
-    VM.copy y jac
-
-  ipH <- wrapIpH (denseIJh n m)
-     ( \ obj_factor lambda x values -> do
-
-        x <- VG.convert `fmap` VS.unsafeFreeze x
-        lambda <- VG.convert `fmap` VS.unsafeFreeze lambda
-
-        let tri = VG.concat . VG.toList . V.imap (\n -> V.take (n+1))
-            obj = V.map (*obj_factor) $ tri $ hessian f x
-            gj = V.zipWith (\l v -> V.map (l*) v) lambda (V.map tri (hessianF g x))
-            lagrangian = V.foldl (V.zipWith (+)) obj gj
-
-        VM.copy values =<< VS.unsafeThaw (VG.convert lagrangian)
-        )
-
-  return (ipF, ipGradF, ipG, ipJacG, ipH)
-
-
--- | indexes the same as http://www.coin-or.org/Ipopt/documentation/node40.html
-denseIJ n m iRow jCol = do
-    VM.copy iRow =<< VS.unsafeThaw (VS.generate (n*m) (\x -> fromIntegral $ x `div` n))
-    VM.copy jCol =<< VS.unsafeThaw (VS.generate (n*m) (\x -> fromIntegral $ x `mod` n))
-
--- | indexes the same as http://www.coin-or.org/Ipopt/documentation/node41.html
-denseIJh n m iRow jCol = do
-    i <- newIORef 0
-    forM_ [0 .. fromIntegral n-1] $ \ row ->
-        forM_ [ 0 .. row ] $ \col -> do
-            ii <- readIORef i
-            VM.write iRow ii row
-            VM.write jCol ii col
-            writeIORef i (ii+1)
-
diff --git a/Ipopt.hs b/Ipopt.hs
new file mode 100644
--- /dev/null
+++ b/Ipopt.hs
@@ -0,0 +1,36 @@
+-- | Description: exports most things you should need 
+--
+-- This module exports most things you should need.
+-- Also take a look at "Ipopt.NLP" and "Ipopt.Raw" and @examples/@
+module Ipopt (
+  -- * high-level
+  NLPT, nlpstate0,
+  module Control.Monad.State,
+  solveNLP',
+  -- ** variables
+  var', var, varFresh,
+  AnyRF(..),
+  -- ** functions
+  addG, addF,
+  -- * low-level bits still needed
+  IpOptSolved(..),
+  ApplicationReturnStatus(..),
+
+  -- ** solver options
+  -- $solverOptRef
+  addIpoptNumOption,
+  addIpoptStrOption, 
+  addIpoptIntOption,
+  setIpoptProblemScaling,
+  openIpoptOutputFile,
+
+  -- *** types
+  Vec, IpNumber,
+) where
+
+import Ipopt.NLP
+import Ipopt.Raw
+import Control.Monad.State
+
+-- $solverOptRef
+-- see <http://www.coin-or.org/Ipopt/documentation/node39.html ipopt's options reference>
diff --git a/Ipopt/NLP.hs b/Ipopt/NLP.hs
new file mode 100644
--- /dev/null
+++ b/Ipopt/NLP.hs
@@ -0,0 +1,315 @@
+{-# LANGUAGE FlexibleInstances, GADTs, RankNTypes, TemplateHaskell #-}
+-- | Description : a EDSL for describing nonlinear programs
+--
+-- see usage in @examples/Test3.hs@
+module Ipopt.NLP where
+
+import Control.Applicative
+import Control.Lens
+import Control.Monad
+import Control.Monad.Identity
+import Control.Monad.State
+import Data.Foldable (toList)
+import Data.IntMap (IntMap)
+import Data.List
+import Data.Monoid
+import Data.Monoid (First)
+import Data.Sequence (Seq)
+import Data.Vector (Vector)
+import Foreign.C.Types (CDouble(..))
+import qualified Data.Foldable as F
+import qualified Data.IntMap as IM
+import qualified Data.Map as M
+import qualified Data.Sequence as Seq
+import qualified Data.Set as S
+import qualified Data.Vector as V
+import qualified Data.Vector.Storable as VS
+
+import Ipopt.Raw
+
+-- * state
+data NLPFun = NLPFun
+    { _funF, _funG :: AnyRF Seq,
+      _boundX, _boundG :: Seq (Double,Double) }
+instance Show NLPFun where
+    show (NLPFun f g a b) = "NLPFun <f> <g> {" ++ show a ++ "}{" ++ show b ++ "}"
+
+data NLPState = NLPState
+    { -- | current maximum index
+      _nMax :: Ix,
+      -- | what namespace (see 'inEnv')
+      _currentEnv :: [String],
+      _variables :: M.Map String Ix,
+      -- | human-readable descriptions for the constraint, objective and
+      -- variables
+      _constraintLabels, _objLabels, _varLabels :: IntMap String,
+      _varEnv :: IntMap (S.Set [String]),
+      _constraintEnv, _objEnv :: IntMap [String],
+      _nlpfun :: NLPFun,
+      _defaultBounds :: (Double,Double),
+      -- | inital state variable for the solver
+      _initX :: Vector Double }
+    deriving (Show)
+
+-- | solver deals with arrays. This type is for indexes into the array
+-- for the current variables that the solver is trying to find.
+newtype Ix = Ix { _varIx :: Int } deriving Show
+
+-- | the initial state to use when you actually have to get to IO
+-- with the solution
+nlpstate0 = NLPState (Ix (-1)) mempty mempty
+    mempty mempty mempty -- labels
+    mempty mempty mempty -- env
+    (mempty :: NLPFun)
+    (-1/0, 1/0)
+    V.empty
+
+type NLPT = StateT NLPState
+type NLP = NLPT IO
+   
+-- ** representing functions
+
+{- | this wrapper holds functions that can be used
+for the objective (`f`) or for constraints (`g`). Many functions
+in the instances provided are partial: this seems to be unavoidable
+because the input variables haven't been decided yet, so you should
+not be allowed to use 'compare' on these. But for now just use the
+standard Prelude classes, and unimplementable functions (which
+would not produce an 'AnyRF') are calls to 'error'
+
+generate these using 'var', or perhaps by directly using the constructor:
+@AnyRF $ Identity . V.sum@, would for example give the sum of all variables.
+-}
+data AnyRF cb = AnyRF (forall a. RealFloat a => Vector a -> cb a)
+
+-- *** helpers for defining instances
+liftOp0 :: (forall a. RealFloat a => a) -> AnyRF Identity
+liftOp0 op = AnyRF $ \x -> Identity op
+
+liftOp1 :: (forall a. RealFloat a => a -> a) -> AnyRF Identity -> AnyRF Identity
+liftOp1 op (AnyRF a) = AnyRF $ \x -> Identity (op (runIdentity (a x)))
+
+liftOp2 :: (forall a. RealFloat a => a -> a -> a) -> AnyRF Identity -> AnyRF Identity -> AnyRF Identity
+liftOp2 op (AnyRF a) (AnyRF b) = AnyRF $ \x -> Identity (runIdentity (a x) `op` runIdentity (b x))
+
+instance Num (AnyRF Identity) where
+    (+) = liftOp2 (+)
+    (*) = liftOp2 (*)
+    (-) = liftOp2 (-)
+    abs = liftOp1 abs
+    signum = liftOp1 signum
+    fromInteger n = liftOp0 (fromInteger n)
+
+instance Fractional (AnyRF Identity) where
+    (/) = liftOp2 (/)
+    recip = liftOp1 recip
+    fromRational n = liftOp0 (fromRational n)
+
+instance Floating (AnyRF Identity) where
+    pi = liftOp0 pi
+    exp = liftOp1 exp
+    sqrt = liftOp1 sqrt
+    log = liftOp1 log
+    sin = liftOp1 sin
+    tan = liftOp1 tan
+    cos = liftOp1 cos
+    asin = liftOp1 asin
+    atan = liftOp1 atan
+    acos = liftOp1 acos
+    sinh = liftOp1 sinh
+    tanh = liftOp1 tanh
+    cosh = liftOp1 cosh
+    asinh = liftOp1 asinh
+    atanh = liftOp1 atanh
+    acosh = liftOp1 acosh
+    (**) = liftOp2 (**)
+    logBase = liftOp2 logBase
+
+instance Real (AnyRF Identity) where
+    toRational _ = error "Real AnyRF Identity"
+
+instance Ord (AnyRF Identity)
+instance Eq (AnyRF Identity)
+instance RealFrac (AnyRF Identity) where
+    properFraction = error "properFraction AnyRF"
+instance RealFloat (AnyRF Identity) where
+    isInfinite = error "isInfinite AnyRF"
+    isNaN = error "isNaN AnyRF"
+    decodeFloat = error "decodeFloat AnyRF"
+    floatRange = error "floatRange AnyRF"
+    isNegativeZero = error "isNegativeZero AnyRF"
+    isIEEE = error "isIEEE AnyRF"
+    isDenormalized = error "isDenormalized AnyRF"
+    floatDigits _ = floatDigits (error "RealFrac AnyRF Identity floatDigits" :: Double)
+    floatRadix _ = floatRadix   (error "RealFrac AnyRF Identity floatRadix" :: Double)
+    atan2 = liftOp2 atan2
+    significand = liftOp1 significand
+    scaleFloat n = liftOp1 (scaleFloat n)
+    encodeFloat a b = liftOp0 (encodeFloat a b)
+
+instance Monoid (AnyRF Seq) where
+    AnyRF f `mappend` AnyRF g = AnyRF (f `mappend` g)
+    mempty = AnyRF mempty
+
+instance Monoid NLPFun where
+    NLPFun f g bx bg `mappend` NLPFun f' g' bx' bg' = NLPFun (f <> f') (g <> g') (bx <> bx') (bg <> bg')
+    mempty = NLPFun mempty mempty mempty mempty
+
+-- ** low level lenses to NLPState
+makeLenses ''NLPFun
+makeLenses ''Ix
+makeLenses ''NLPState
+
+-- | @to@ is one of 'varEnv', 'constraintEnv', 'objEnv'
+copyEnv to n = do
+    ev <- use currentEnv
+    cloneLens to %= IM.insert n ev
+
+-- | @to@ should be 'constraintLabels', 'objLabels', 'varLabels'
+addDesc to Nothing n = return ()
+addDesc to (Just x) n = do
+    to %= IM.insert n x
+
+-- * high-level functions
+
+-- | calls 'createIpoptProblemAD' and 'ipoptSolve'. To be used at the
+-- end of a do-block.
+solveNLP' :: MonadIO m =>
+    (IpProblem -> IO ()) -- ^ set ipopt options (using functions from "Ipopt.Raw")
+    -> NLPT m IpOptSolved
+solveNLP' setOpts = do
+    (xl,xu) <- join $ uses (nlpfun . boundX) seqToVecs
+    (gl,gu) <- join $ uses (nlpfun . boundG) seqToVecs
+    AnyRF fs <- use (nlpfun . funF)
+    AnyRF gs <- use (nlpfun . funG)
+
+    p <- liftIO (createIpoptProblemAD xl xu gl gu (F.sum . fs) (V.fromList . toList . gs))
+    liftIO (setOpts p)
+    x0 <- uses initX (V.convert . V.map CDouble)
+    r <- liftIO (ipoptSolve p =<< VS.thaw x0)
+    liftIO $ freeIpoptProblem p
+    return r
+
+-- | add a constraint
+addG :: Monad m
+    => Maybe String -- ^ optional description
+    -> (Double,Double) -- ^ bounds @(gl,gu)@ for the single inequality @gl_i <= g_i(x) <= gu_i@
+    -> AnyRF Identity -- ^ @g_i(x)@
+    -> NLPT m ()
+addG d b (AnyRF f) = do
+    nlpfun . boundG %= (Seq.|> b)
+    nlpfun . funG %= \(AnyRF fs) -> AnyRF $ \x -> fs x Seq.|> runIdentity (f x)
+    n <- use (nlpfun . boundG . to Seq.length)
+    copyEnv constraintEnv n
+    addDesc constraintLabels d n
+
+{- | add a piece of the objective function, which is added in the form
+`f_1 + f_2 + ...`, to make it easier to understand (at some point)
+which components are responsible for the majority of the cost, and
+which are irrelevant.
+-}
+addF :: Monad m
+    => Maybe String -- ^ description
+    -> AnyRF Identity -- ^ `f_i(x)`
+    -> NLPT m ()
+addF d (AnyRF f) = do
+    nlpfun . funF %= \(AnyRF fs) -> AnyRF $ \x -> fs x Seq.|> runIdentity (f x)
+    n <- use (objEnv . to ((+1) . IM.size))
+    copyEnv objEnv n
+    addDesc objLabels d n
+
+-- | add a variable, or get a reference to the the same variable if it has
+-- already been used
+var' :: (Monad m, Functor m)
+    => Maybe (Double,Double) -- ^ bounds @(xl,xu)@ to request that @xl <= x <= xu@.
+                             -- if Nothing, you get whatever is in 'defaultBounds'
+    -> String -- ^ variable name (namespace from the 'pushEnv' / 'popEnv' can
+              -- make an @"x"@ you request here different from one you
+              -- previously requested
+    -> NLPT m (AnyRF Identity, Ix) -- ^ the value, and index (into the raw
+                                   -- vector of variables that the solver
+                                   -- sees)
+var' bs s = do
+    ev <- use currentEnv
+    m <- use variables
+    let s' = intercalate "." (reverse (s:ev))
+    n <- case M.lookup s' m of
+        Nothing -> do
+            nMax %= over varIx (+1)
+            db <- use defaultBounds
+            nlpfun . boundX %= (Seq.|> db)
+            n' <- use nMax
+            variables %= M.insert s' n'
+            return n'
+        Just n -> return n
+    varEnv %= IM.insert (view varIx n) (S.singleton ev)
+    F.traverse_ (narrowBounds n) bs
+    return (AnyRF $ \x -> Identity $ x V.! view varIx n, n)
+
+-- | 'var'' without the usually unnecessary 'Ix'
+var bs s = fmap fst (var' bs s)
+
+{- | 'var', except this causes the solver to get a new variable,
+so that you can use:
+
+> [a,b,c,d,e] <- replicateM 5 (varFresh (Just (0, 10)) "x")
+
+and the different letters can take different values (between 0 and 10)
+in the optimal solution (depending on what you do with @a@ and similar
+in the objective function).
+-}
+varFresh bs s = do
+    n <- uses variables ((+1) . M.size)
+    var bs (s ++ show n)
+
+-- *** namespace
+{- $namespace
+
+When you build up an optimization problem, it may be composed of pieces.
+Functions in this section help to ease the pain of making unique variables.
+To illustrate:
+
+> m <- inEnv "A" (var b "x")
+> n <- var b "A.x" 
+
+@m@ and @n@ above should refer to the same variable. In some sense this
+is \"better\" that using 'varFresh' all the time, since perhaps you would
+like to add dependencies between components (say the size of a header pipe,
+refridgeration unit, foundation etc. has to satisfy sizes of individual
+components).
+
+-}
+
+-- | combination of 'pushEnv' and 'popEnv'
+inEnv :: Monad m => String -> NLPT m a -> NLPT m a
+inEnv s action = do
+    pushEnv s
+    r <- action
+    popEnv
+    return r
+
+pushEnv :: Monad m => String -> NLPT m ()
+pushEnv s = currentEnv %= (s:)
+
+popEnv :: Monad m => NLPT m String
+popEnv = do
+    n : ns <- use currentEnv
+    currentEnv .= ns
+    return n
+
+-- *** bounds
+
+-- | override bounds. Should be unnecessary given 'var' takes bounds.
+setBounds :: Monad m => Ix -> (Double,Double) -> NLPT m ()
+setBounds (Ix i) bs = nlpfun . boundX %= Seq.update i bs
+
+-- | shrink the interval in which that variable is allowed.
+narrowBounds :: Monad m => Ix -> (Double,Double) -> NLPT m ()
+narrowBounds (Ix i) (a,b) = nlpfun . boundX . ix i %= \(a',b') -> (max a a', min b b')
+
+-- * internal
+seqToVecs :: MonadIO m => Seq (Double,Double) -> m (Vec,Vec)
+seqToVecs x = let (a,b) = unzip (toList x) in liftIO $ do
+    a' <- VS.thaw (VS.map CDouble (VS.fromList a))
+    b' <- VS.thaw (VS.map CDouble (VS.fromList b))
+    return (a',b')
diff --git a/Ipopt/Raw.chs b/Ipopt/Raw.chs
new file mode 100644
--- /dev/null
+++ b/Ipopt/Raw.chs
@@ -0,0 +1,357 @@
+{-#  LANGUAGE ForeignFunctionInterface, PatternGuards, RankNTypes, TypeFamilies #-}
+{- |
+
+Copyright: (C) 2013 Adam Vogt
+Maintainer: Adam Vogt <vogt.adam@gmail.com>
+Stability: unstable
+Description: lowest-level parts of the binding
+
+
+-}
+module Ipopt.Raw (
+    -- * specifying problem
+    createIpoptProblemAD,
+
+    -- ** solve
+    ipoptSolve,
+    IpOptSolved(..),
+
+    -- ** solver options
+    addIpoptNumOption,
+    addIpoptStrOption,
+    addIpoptIntOption,
+    openIpoptOutputFile,
+
+    -- * types
+    Vec,
+
+    IpNumber(..),
+    IpIndex(..),
+    IpInt(..),
+    IpBool(..),
+
+    IpF(..),
+    IpGradF(..),
+    IpG(..),
+    IpJacG(..),
+    IpH(..),
+
+    IpProblem(..),
+
+    ApplicationReturnStatus(..),
+
+    -- * lower-level parts of the binding
+    createIpoptProblem,
+
+    freeIpoptProblem,
+    setIpoptProblemScaling,
+
+
+    -- ** marshalling functions
+    wrapIpF,
+    wrapIpGradF,
+    wrapIpG,
+    wrapIpJacG,
+    wrapIpH,
+
+    wrapIpF1,
+    wrapIpGradF1,
+    wrapIpG1,
+    wrapIpJacG1,
+    wrapIpH1,
+
+    wrapIpF2,
+    wrapIpGradF2,
+    wrapIpG2,
+    wrapIpJacG2,
+    wrapIpH2,
+
+    ) where
+
+import Ipopt.Sparsity
+
+import C2HS
+import Control.Exception
+import Control.Monad
+import Data.IORef
+import Foreign.C
+import Foreign.ForeignPtr
+import Foreign.Ptr
+import Foreign.Storable
+import Numeric.AD
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as VG
+import qualified Data.Vector.Storable as VS
+import qualified Data.Vector.Storable.Mutable as VM
+
+#include "IpStdCInterface.h"
+
+type IpNumber = {# type Number #}
+type IpIndex = {# type Index #}
+type IpInt = {# type Int #}
+type IpBool = {# type Bool #}
+
+type IpF = {# type Eval_F_CB #}
+type IpGradF = {# type Eval_Grad_F_CB #}
+type IpG = {# type Eval_G_CB #}
+type IpJacG = {# type Eval_Jac_G_CB #}
+type IpH = {# type Eval_H_CB #}
+
+{#enum ApplicationReturnStatus as ^ {underscoreToCase} deriving (Show) #}
+{#enum AlgorithmMode as ^ {underscoreToCase} deriving (Show) #}
+
+
+newtype IpProblem = IpProblem { unIpProblem :: Ptr ()}
+
+type family UnFunPtr a
+type instance UnFunPtr (FunPtr a) = a
+
+ipTrue = 1 :: IpBool
+ipFalse = 0 :: IpBool
+
+-- | likely an unsafe method for getting a "Data.Vector.Storable.Mutable" out of a 'Ptr'
+ptrToVS n p = do
+    fp <- newForeignPtr_ p
+    return (VM.unsafeFromForeignPtr0 fp (fromIntegral n))
+
+foreign import ccall "wrapper" wrapIpF1 :: UnFunPtr IpF -> IO IpF
+foreign import ccall "wrapper" wrapIpG1 :: UnFunPtr IpG -> IO IpG
+foreign import ccall "wrapper" wrapIpGradF1 :: UnFunPtr IpGradF -> IO IpGradF
+foreign import ccall "wrapper" wrapIpJacG1 :: UnFunPtr IpJacG -> IO IpJacG
+foreign import ccall "wrapper" wrapIpH1 :: UnFunPtr IpH -> IO IpH
+
+
+toB x =  either (\ e@SomeException {} -> print e >> return ipFalse)
+                (\ _ -> return ipTrue ) =<< try x
+
+wrapIpF2' fun n xin new_x obj_val _userData = do
+    toB $ poke obj_val =<< fun =<< ptrToVS n xin
+
+wrapIpF2 fun n xin new_x obj_val _userData = do
+    toB $ poke obj_val =<< fun =<< ptrToVS n xin
+
+wrapIpG2 fun n xin new_x m gout _userData = do
+    toB $ join $ liftM2 VM.copy (ptrToVS m gout) (fun =<< ptrToVS n xin)
+
+wrapIpGradF2 fun n x new_x grad_f _userData = do
+    toB $ join $ liftM2 VM.copy (ptrToVS n grad_f) (fun =<< ptrToVS n x)
+
+wrapIpJacG2 fun1 fun2 n x new_x m nj iRow jCol jacs _userData
+    | jacs == nullPtr = do
+            toB $ join $ liftM2 fun1 (ptrToVS nj iRow) (ptrToVS nj jCol)
+    | otherwise = do
+            toB $ join $ liftM2 fun2 (ptrToVS n x) (ptrToVS nj jacs)
+
+wrapIpH2 funSparsity funEval n x new_x obj_factor m lambda new_lambda nHess iRow jCol values _userData
+    | iRow == nullPtr = do
+            toB $ join $ liftM3 (funEval obj_factor)
+                        (ptrToVS m lambda)
+                        (ptrToVS n x)
+                        (ptrToVS nHess values)
+    | otherwise = do
+            toB $ join $ liftM2 funSparsity
+                        (ptrToVS nHess iRow)
+                        (ptrToVS nHess jCol)
+
+wrapIpF f = wrapIpF1 (wrapIpF2 f)
+wrapIpG f = wrapIpG1 (wrapIpG2 f)
+wrapIpGradF f = wrapIpGradF1 (wrapIpGradF2 f)
+wrapIpJacG f1 f2 = wrapIpJacG1 (wrapIpJacG2 f1 f2)
+wrapIpH fSparsity fEval = wrapIpH1 (wrapIpH2 fSparsity fEval)
+
+
+vmUnsafeWith = VM.unsafeWith
+
+-- | Vector of numbers
+type Vec = VM.IOVector IpNumber
+
+createIpoptProblem :: Vec -> Vec -> Vec -> Vec
+    -> Int -> Int -> IpF -> IpG -> IpGradF -> IpJacG -> IpH -> IO IpProblem
+createIpoptProblem xL xU gL gU nJac nHess f g gradF jacG hess
+    | lx <- VM.length xL,
+      lx == VM.length xU,
+      lg <- VM.length gL,
+      lg == VM.length gU = createIpoptProblem3 lx xL xU lg gL gU nJac nHess 0 f g gradF jacG hess
+   | otherwise = error "dimensions wrong!"
+
+{#fun CreateIpoptProblem as createIpoptProblem3
+        { `Int', vmUnsafeWith* `Vec', vmUnsafeWith* `Vec',
+          `Int', vmUnsafeWith* `Vec', vmUnsafeWith* `Vec',
+          `Int', `Int', `Int', id `IpF', id `IpG', id `IpGradF',
+          id `IpJacG', id `IpH' } -> `IpProblem' IpProblem #}
+_ = {#fun AddIpoptNumOption as ^
+    { unIpProblem `IpProblem', `String', `Double' } -> `Bool' #}
+
+_ = {#fun AddIpoptStrOption as ^
+    { unIpProblem `IpProblem', `String', `String' } -> `Bool' #}
+
+_ = {#fun AddIpoptIntOption as ^
+    { unIpProblem `IpProblem', `String', `Int' } -> `Bool' #}
+
+_ = {#fun FreeIpoptProblem as ^
+    { unIpProblem `IpProblem' } -> `()' #}
+
+_ = {#fun OpenIpoptOutputFile as ^
+    { unIpProblem `IpProblem', `String', `Int' } -> `Bool' #}
+
+_ = {#fun SetIpoptProblemScaling as ^
+    { unIpProblem `IpProblem',
+    `Double',
+     vmUnsafeWith* `Vec',
+     vmUnsafeWith* `Vec'
+     } -> `Bool' #}
+
+
+data IpOptSolved = IpOptSolved
+  { ipOptSolved_status :: ApplicationReturnStatus,
+    ipOptSolved_objective :: Double,
+    ipOptSolved_x,
+    ipOptSolved_g,
+    ipOptSolved_mult_g,
+    ipOptSolved_mult_x_L,
+    ipOptSolved_mult_x_U :: Vec }
+
+ipoptSolve :: IpProblem
+    -> Vec -- ^ starting point @x@. Note that the value is overwritten with the final @x@.
+    -> IO IpOptSolved
+ipoptSolve problem x = do
+    g <- VM.new (VM.length x)
+    mult_g <- VM.new (VM.length x)
+    mult_x_L <- VM.new (VM.length x)
+    mult_x_U <- VM.new (VM.length x)
+
+    out <- ipoptSolve2
+     problem
+     x
+     g
+     mult_g
+     mult_x_L
+     mult_x_U
+     nullPtr
+    return $ IpOptSolved
+      (fst out)
+      (snd out)
+      x
+      g
+      mult_g
+      mult_x_L
+      mult_x_U
+
+_ = {#fun IpoptSolve as ipoptSolve2
+        { unIpProblem `IpProblem',
+        vmUnsafeWith* `Vec',
+        vmUnsafeWith* `Vec',
+        alloca- `Double' peekFloatConv*,
+        vmUnsafeWith* `Vec',
+        vmUnsafeWith* `Vec',
+        vmUnsafeWith* `Vec',
+        id `Ptr ()' } -> `ApplicationReturnStatus' cToEnum #}
+
+
+{- | Set-up an 'IpProblem' to be solved later. Only objective function (@f@)
+and constraint functions (@g@) need to be specified. Derivatives needed by ipopt
+are computed by "Numeric.AD".
+
+To solve the optimization problem:
+
+>              min f(x)
+>   such that
+>           xL <=  x     <= xU
+>           gL <=  g(x)  <= gU
+
+First create an opaque 'IpProblem' object (nlp):
+
+> nlp <- createIpOptProblemAD xL xU gL gU f g
+
+Then pass it off to 'ipoptSolve'.
+
+> ipoptSolve nlp x0
+
+Refer to the example @Test2.hs@ for details of setting up the vectors supplied.
+-}
+createIpoptProblemAD
+    :: Vec -- ^ @xL@ 'VM.Vector' of lower bounds for decision variables with length @n@
+    -> Vec -- ^ @xU@ 'VM.Vector' of upper bounds for decision variables
+    -> Vec -- ^ @gL@ 'VM.Vector' of lower bounds for constraint functions @g(x)@ with length @m@
+    -> Vec -- ^ @gU@ 'VM.Vector' of upper bounds for constraint functions @g(x)@
+    -> (forall a. RealFloat a => V.Vector a -> a) -- ^ objective function @f : R^n -> R@
+    -> (forall a. RealFloat a => V.Vector a -> V.Vector a) -- ^ constraint functions @g : R^n -> R^m@
+    -> IO IpProblem
+createIpoptProblemAD xL xU gL gU f g
+    | n <- VM.length xL,
+      n == VM.length xU,
+      m <- VM.length gL,
+      m == VM.length gU = do
+    (eval_f, eval_grad_f, eval_g, eval_jac_g, eval_h) <- mkFs (True,True) n m f g
+    createIpoptProblem xL xU gL gU (n*m) (((n+1)*n) `div` 2)
+            eval_f eval_g eval_grad_f eval_jac_g eval_h
+
+mkFs ::
+    (Bool,Bool) -- ^ grad, hessian are sparse?
+    -> Int -- ^ @n@ number of variables
+    -> Int -- ^ @m@ number of constraints
+    -> (forall a. RealFloat a => V.Vector a -> a) -- ^ objective function @R^n -> R@
+    -> (forall a. RealFloat a => V.Vector a -> V.Vector a) -- ^ constraint functions @R^n -> R^m@
+    -> IO (IpF, IpGradF, IpG, IpJacG, IpH)
+mkFs (jacGs,hessGs) n m f g = do
+  ipF <- wrapIpF $ \x -> do
+    x <- VG.convert `fmap` VS.unsafeFreeze x
+    return $ f x
+
+  ipGradF <- wrapIpGradF $ \x -> do
+    x <- VG.convert `fmap` VS.unsafeFreeze x
+    VS.unsafeThaw $ VG.convert (grad f x)
+
+  ipG <- wrapIpG $ \x -> do
+    x <- VG.convert `fmap` VS.unsafeFreeze x
+    VS.unsafeThaw $ VG.convert (g x)
+
+  let ijJ = nanPropagateJ n g
+  ipJacG <- 
+   if jacGs then wrapIpJacG (sparseIJ ijJ) $ \x y -> do
+    x <- VG.convert `fmap` VS.unsafeFreeze x
+    jac <- VS.unsafeThaw $ VG.convert $ jacobianSS ijJ g x
+    VM.copy y jac
+   else wrapIpJacG (denseIJ n m) $ \x y -> do
+    x <- VG.convert `fmap` VS.unsafeFreeze x
+    jac <- VS.unsafeThaw $ VG.convert $ VG.concat $ VG.toList $ jacobian g x
+    VM.copy y jac
+
+  ipH <- wrapIpH (denseIJh n m)
+     ( \ obj_factor lambda x values -> do
+
+        x <- VG.convert `fmap` VS.unsafeFreeze x
+        lambda <- VG.convert `fmap` VS.unsafeFreeze lambda
+
+        let tri = VG.concat . VG.toList . V.imap (\n -> V.take (n+1))
+            obj = V.map (*obj_factor) $ tri $ hessian f x
+            gj = V.zipWith (\l v -> V.map (l*) v) lambda (V.map tri (hessianF g x))
+            lagrangian = V.foldl (V.zipWith (+)) obj gj
+
+        VM.copy values =<< VS.unsafeThaw (VG.convert lagrangian)
+        )
+
+  return (ipF, ipGradF, ipG, ipJacG, ipH)
+
+
+sparseIJ ijJ iRow jCol = do
+    let (iJ,jJ) = V.unzip ijJ
+        conv = VS.unsafeThaw . V.convert . V.map fromIntegral
+    VM.copy iRow =<< conv iJ
+    VM.copy jCol =<< conv jJ
+
+-- | indexes the same as http://www.coin-or.org/Ipopt/documentation/node40.html
+denseIJ n m iRow jCol = do
+    VM.copy iRow =<< VS.unsafeThaw (VS.generate (n*m) (\x -> fromIntegral $ x `div` n))
+    VM.copy jCol =<< VS.unsafeThaw (VS.generate (n*m) (\x -> fromIntegral $ x `mod` n))
+
+-- | indexes the same as http://www.coin-or.org/Ipopt/documentation/node41.html
+denseIJh n m iRow jCol = do
+    i <- newIORef 0
+    forM_ [0 .. fromIntegral n-1] $ \ row ->
+        forM_ [ 0 .. row ] $ \col -> do
+            ii <- readIORef i
+            VM.write iRow ii row
+            VM.write jCol ii col
+            writeIORef i (ii+1)
+
diff --git a/Ipopt/Sparsity.hs b/Ipopt/Sparsity.hs
new file mode 100644
--- /dev/null
+++ b/Ipopt/Sparsity.hs
@@ -0,0 +1,155 @@
+{-# LANGUAGE PatternGuards, RankNTypes, TupleSections #-}
+{- |
+
+Description: find out derivatives are unnecessary
+
+gradients and hessians tend to be sparse. This can be probed by
+calculating these with NaNs as inputs. http://www.gpops2.com/ uses this
+strategy, but perhaps there are relatively common cases (calls to
+blas/lapack for example) that do not propagate NaNs correctly.
+
+All functions provide indices that are affected (and should thus be included)
+
+functions ending in 1 set one variable at a time to NaN, and additionally
+provide a hint as to which variable to change.
+
+functions not ending in 1 set all input variables to NaN
+
+it seems that adding a 1 achieves the same thing as doing one more derivative.
+In other words, nanPropagateG1 tells the same thing as nanPropagateH.
+
+-}
+module Ipopt.Sparsity (
+    -- $setup
+    -- * gradient sparsity
+    nanPropagate1,
+    nanPropagateG,
+
+    -- * jacobian
+    nanPropagateJ,
+    -- * hessian sparsity
+    nanPropagateH,
+    nanPropagateG1,
+
+    nanPropagateHF,
+
+    -- * sparse derivatives (faked for now)
+    jacobianSS,
+    ) where
+
+import Control.Applicative
+import Control.Monad
+import qualified Data.Vector as V
+import Data.Vector ((!))
+import Numeric.AD
+
+import qualified Data.IntMap as M
+import qualified Data.IntSet as S
+import Data.Monoid
+
+{- $setup
+>>> let g x = x!0 + x!1 * x!1 * x!2
+
+-}
+
+{- | a nonzero gradient when inputs are NaN ==> no need to include
+that row/column in the hessian, since it will be zero
+
+>>> nanPropagateG1 4 g
+[(1,fromList [1,2]),(2,fromList [1])]
+
+
+
+-}
+nanPropagateG1 :: Int -- ^ size of input vector
+    -> (forall a. RealFloat a => V.Vector a -> a)
+    -> [(Int, V.Vector Int)] -- ^ @(i,js)@
+        -- shows that a NaN at index @i@ produces NaNs in gradient indexes
+        -- @js@... so the hessians only need to include
+        -- @i@ and @js@
+nanPropagateG1 nx f | x0:_ <- dropWhile (isNaN . f) (trialV nx [0,0.5,1])
+    = [ (i,j) | i <- [0 .. nx-1],
+            let g = grad f (x0 V.// [(i,(0/0))])
+                j = V.findIndices isNaN g,
+            not (V.null j) ]
+
+{- |
+
+>>> nanPropagateG 4 g
+fromList [0,1,2]
+
+-}
+nanPropagateG :: Int -- ^ size of input vector
+    -> (forall a. RealFloat a => V.Vector a -> a)
+    -> V.Vector Int
+nanPropagateG nx f = V.findIndices (\x -> x /= 0)
+        (grad f (V.replicate nx (0/0)))
+
+{- |
+>>> nanPropagate1 4 g
+[0,1,2]
+
+variable 3 isn't even used.
+-}
+nanPropagate1 :: Int -- ^ size of input vector
+    -> (forall a. RealFloat a => V.Vector a -> a)
+    -> [Int] -- ^ inputs that don't become NaN
+nanPropagate1 nx f | x0:_ <- dropWhile (isNaN . f) (trialV nx [0,0.5,1])
+    = [ i | i <- [0 .. nx-1],
+            let g = f (x0 V.// [(i,(0/0))]),
+            isNaN g ]
+
+
+{- |
+>>> nanPropagateH 4 g
+fromList [(1,1),(1,2),(2,1)]
+
+-}
+nanPropagateH :: Int
+    -> (forall a. RealFloat a => V.Vector a -> a)
+    -> V.Vector (Int, Int) -- ^ (i,j) indexes
+nanPropagateH nx f = nonzeroIxs $ hessian f (V.replicate nx (0/0::Double))
+
+nanPropagateJ :: Int
+    -> (forall a. RealFloat a => V.Vector a -> V.Vector a)
+    -> V.Vector (Int,Int)
+nanPropagateJ nx f = nonzeroIxs $ jacobian f (V.replicate nx (0/0::Double))
+
+nanPropagateHF :: Int
+    -> (forall a. RealFloat a => V.Vector a -> V.Vector a)
+    -> V.Vector (V.Vector (Int,Int))
+nanPropagateHF nx f =
+        V.map nonzeroIxs $ hessianF f (V.replicate nx (0/0::Double))
+
+{-
+complement :: (Int,Int) -> V.Vector (Int,Int) -> V.Vector (Int,Int)
+complement (mx,my) = M.fromList . M.toList
+-}
+
+collapse :: V.Vector (V.Vector (Int,Int)) -> V.Vector (Int,Int)
+collapse = V.fromList . concatMap (\(a,bs) -> map (a,) (S.toList bs)) . M.toList .
+            M.mapWithKey (\k s -> fst (S.split (k+1) s)) .
+            M.fromListWith (<>) . map (\(a,b) -> (a, S.singleton b)) . concatMap V.toList . V.toList
+
+nonzeroIxs = V.concatMap id . V.imap (\i -> V.map (i,) . V.findIndices (/= 0))
+
+trialV :: Int -> [Double] -> [V.Vector Double]
+trialV nx cand = V.fromList <$> replicateM nx cand
+
+
+
+jacobianSS :: RealFloat a => V.Vector (Int,Int)
+    -> (forall a. RealFloat a => V.Vector a -> V.Vector a)
+    -> V.Vector a
+    -> V.Vector a
+jacobianSS ijs f x = let v = jacobian f x
+    in V.map (\(i,j) -> v V.! i V.! j) ijs
+
+{-
+hessianFSS :: RealFloat a => V.Vector (Int,Int)
+    -> (forall a. RealFloat a => V.Vector a -> V.Vector a)
+    -> V.Vector a
+    -> V.Vector a
+hessianFSS = undefined
+-}
+
diff --git a/Test1.hs b/Test1.hs
deleted file mode 100644
--- a/Test1.hs
+++ /dev/null
@@ -1,123 +0,0 @@
-{- | 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
-
diff --git a/Test2.hs b/Test2.hs
deleted file mode 100644
--- a/Test2.hs
+++ /dev/null
@@ -1,52 +0,0 @@
-{-# LANGUAGE RankNTypes #-}
-{- | Translation of the "number 71 from the Hock-Schittkowsky test suite"
-<http://www.coin-or.org/Ipopt/documentation/node28.html> testing AD
-
--}
-module Main where
-
-import qualified Data.Vector.Storable as VS
-import Text.Printf ( printf )
-import Data.Vector ((!))
-import qualified Data.Vector as V
-import qualified Data.Vector.Storable as VS
-import Ipopt
-    ( addIpoptNumOption,
-      addIpoptStrOption,
-      ipoptSolve,
-      createIpoptProblemAD )
-
-main = do
-    let n = 4 -- number of variables
-        fromList = VS.unsafeThaw . VS.fromList
-
-    xL <- fromList (replicate n 1)
-    xU <- fromList (replicate n 5)
-
-    gL <- fromList [25, 40]
-    gU <- fromList [2.0e19,40]
-
-    let f x = x!0 * x!3 * (x!0 + x!1 + x!2) + x!2
-        g x = V.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]
-
-    nlp <- createIpoptProblemAD xL xU gL gU f g
-
-    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
-
-
diff --git a/examples/Test1.hs b/examples/Test1.hs
new file mode 100644
--- /dev/null
+++ b/examples/Test1.hs
@@ -0,0 +1,123 @@
+{- | 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.Raw
+    ( 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
+
diff --git a/examples/Test2.hs b/examples/Test2.hs
new file mode 100644
--- /dev/null
+++ b/examples/Test2.hs
@@ -0,0 +1,52 @@
+{-# LANGUAGE RankNTypes #-}
+{- | Translation of the "number 71 from the Hock-Schittkowsky test suite"
+<http://www.coin-or.org/Ipopt/documentation/node28.html> testing AD
+
+-}
+module Main where
+
+import qualified Data.Vector.Storable as VS
+import Text.Printf ( printf )
+import Data.Vector ((!))
+import qualified Data.Vector as V
+import qualified Data.Vector.Storable as VS
+import Ipopt.Raw
+    ( addIpoptNumOption,
+      addIpoptStrOption,
+      ipoptSolve,
+      createIpoptProblemAD )
+
+main = do
+    let n = 4 -- number of variables
+        fromList = VS.unsafeThaw . VS.fromList
+
+    xL <- fromList (replicate n 1)
+    xU <- fromList (replicate n 5)
+
+    gL <- fromList [25, 40]
+    gU <- fromList [2.0e19,40]
+
+    let f x = x!0 * x!3 * (x!0 + x!1 + x!2) + x!2
+        g x = V.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]
+
+    nlp <- createIpoptProblemAD xL xU gL gU f g
+
+    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
+
+
diff --git a/examples/Test3.hs b/examples/Test3.hs
new file mode 100644
--- /dev/null
+++ b/examples/Test3.hs
@@ -0,0 +1,20 @@
+module Main where
+import Ipopt
+import qualified Data.Vector as V
+import Control.Lens
+import Control.Monad
+import Control.Monad.State
+
+-- | the same problem as solved in Test1 and Test2, except
+-- showing how bits from Ipopt.NLP help
+test22 = do
+   [x0, x1, x2, x3] <- replicateM 4 (varFresh (Just (1,5)) "x")
+   addF (Just "test22")      $ x0*x3*(x0 + x1 + x2) + x2
+   addG Nothing (25,2.0e19)  $ x0*x1*x2*x3
+   addG Nothing (40,40)      $ x0*x0 + x1*x1 + x2*x2 + x3*x3
+   initX .= V.fromList [1,5,5,1]
+   solveNLP' (const (return ()))
+ `runStateT` nlpstate0
+
+main = test22
+
diff --git a/ipopt-hs.cabal b/ipopt-hs.cabal
--- a/ipopt-hs.cabal
+++ b/ipopt-hs.cabal
@@ -1,25 +1,75 @@
 name:                ipopt-hs
-version:             0.0.0.0
+version:             0.2.0.0
 synopsis:            haskell binding to ipopt including automatic differentiation
-description:         
+description:  a haskell binding to the nonlinear programming solver ipopt
+              <http://projects.coin-or.org/Ipopt>
+              .
+              [@installation@]
+              needs ipopt installed: For example if you have a
+              `/usr/include/coin/IpStdCInterface.h` from your ipopt installation,
+              use:
+             .
+              > cabal install ipopt-hs --extra-include-dirs=/usr/include/coin
+             .
+              A embedded language, similar to the one provided by glpk-hs, is
+              defined in "Ipopt.NLP". The goal is to define problems at a level
+              similar to other "algebraic modeling languages", but retain some
+              of the safety and flexibility available in haskell.
+             .
+              Refer to @examples/Test1.hs@ for an example where the derivatives
+              are computed by hand, @Test2.hs@ for the use of
+              'createIpoptProblemAD' and @Test3.hs@ for the highest level.
+             .
+              Current limitations include:
+             .
+              * copying in every iteration happens between between
+                "Data.Vector.Storable" and "Data.Vector" might be avoidable
+                somehow.  Currently it is done because AD needs a Traversable
+                structure, but Storable vectors are not traversable.
+             .
+              * sparseness of derivatives isn't used to decide which way to calculate (forward vs. backward mode)
+             .
+              * probably doesn't work if @IpStdCInterface.h@ has Number =/= 'CDouble'
+             .
+              * no binding to SetIntermediateCallback
+             .
+              * garbage collection of 'IpProblem' won't free C-side resources
 license:             BSD3
 license-file:        LICENSE
 author:              Adam Vogt <vogt.adam@gmail.com>
 maintainer:          Adam Vogt <vogt.adam@gmail.com>
-category:            Data
+category:            Optimisation Math Numeric
 build-type:          Simple
 cabal-version:       >=1.8
 
+source-repository head
+    type:   darcs
+    location: http://code.haskell.org/~aavogt/ipopt-hs
+
 library
-  exposed-modules:     Ipopt
+  exposed-modules:     Ipopt, Ipopt.Raw, Ipopt.Sparsity, Ipopt.NLP
   other-modules:       C2HS
-  build-depends:       base ==4.6.*, vector ==0.10.*, ad ==3.4.*
+  build-depends:       base ==4.6.*,
+                       vector ==0.10.*,
+                       ad ==3.4.*,
+                       containers,
+                       mtl == 2.*,
+                       lens
   pkgconfig-depends:   ipopt
+  build-tools: c2hs
 
-executable ipopt-hs_Test1 
+executable ipopt-hs_Test1
   main-is:             Test1.hs
   build-depends:       base ==4.6.*, vector ==0.10.*, ipopt-hs
+  hs-source-dirs:      examples
+  other-modules: Paths_ipopt_hs
 
 executable ipopt-hs_Test2
   main-is:             Test2.hs
   build-depends:       base ==4.6.*, vector ==0.10.*, ipopt-hs
+  hs-source-dirs:      examples
+
+executable ipopt-hs_Test3
+  main-is:             Test3.hs
+  build-depends:       base ==4.6.*, vector ==0.10.*, ipopt-hs, lens, mtl
+  hs-source-dirs:      examples