diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2013, Jonathan Fischoff
+
+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 Jonathan Fischoff 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/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/lagrangian.cabal b/lagrangian.cabal
new file mode 100644
--- /dev/null
+++ b/lagrangian.cabal
@@ -0,0 +1,85 @@
+-- Initial lagrangian.cabal generated by cabal init.  For further 
+-- documentation, see http://haskell.org/cabal/users-guide/
+
+-- The name of the package.
+name:                lagrangian
+
+-- The package version.  See the Haskell package versioning policy (PVP) 
+-- for standards guiding when and how versions should be incremented.
+-- http://www.haskell.org/haskellwiki/Package_versioning_policy
+-- PVP summary:      +-+------- breaking API changes
+--                   | | +----- non-breaking API additions
+--                   | | | +--- code changes with no API change
+version:             0.1.0.0
+
+-- A short (one-line) description of the package.
+synopsis:            Solve lagrangian multiplier problems
+
+-- A longer description of the package.
+-- description:         
+
+-- URL for the project homepage or repository.
+homepage:            http://github.com/jfischoff/lagrangian
+
+-- The license under which the package is released.
+license:             BSD3
+
+-- The file containing the license text.
+license-file:        LICENSE
+
+-- The package author(s).
+author:              Jonathan Fischoff
+
+-- An email address to which users can send suggestions, bug reports, and 
+-- patches.
+maintainer:          jonathangfischoff@gmail.com
+
+-- A copyright notice.
+-- copyright:           
+
+category:            Math
+
+build-type:          Simple
+
+-- Constraint on the version of Cabal needed to build this package.
+cabal-version:       >=1.8
+
+
+library
+  -- Modules exported by the library.
+  exposed-modules: Numeric.AD.Lagrangian
+  
+  -- Modules included in this library but not exported.
+  other-modules: Numeric.AD.Lagrangian.Internal      
+  
+  -- Other library packages from which modules are imported.
+  build-depends:    base ==4.6.*, 
+                    nonlinear-optimization ==0.3.*, 
+                    vector ==0.10.*, 
+                    ad ==3.3.*,
+                    hmatrix == 0.14.*
+  
+  -- Directories containing source files.
+  hs-source-dirs:      src
+
+Test-Suite tests
+  Hs-Source-Dirs: src, tests
+  type:       exitcode-stdio-1.0
+  main-is:    Main.hs
+  build-depends: base ==4.6.*,
+                 nonlinear-optimization ==0.3.*, 
+                 vector ==0.10.*, 
+                 ad ==3.3.*,
+                 hmatrix == 0.14.*, 
+                 test-framework ==0.6.*, 
+                 test-framework-hunit ==0.2.*, 
+                 test-framework-quickcheck2 ==0.2.*,
+                 HUnit == 1.2.*
+
+
+
+
+
+
+
+ 
diff --git a/src/Numeric/AD/Lagrangian.hs b/src/Numeric/AD/Lagrangian.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Lagrangian.hs
@@ -0,0 +1,27 @@
+-- |Numerically solve convex lagrange multiplier problems with conjugate gradient descent. 
+-- 
+--  Convexity is key, otherwise the descent algorithm can return the wrong answer.
+--  
+--  Convexity can be tested by assuring that the hessian of the lagrangian is positive
+--  definite over region the function is defined in. 
+--  
+--  I have provided test that the hessian is positive definite at a point, which is something,
+--  but not enough to ensure that the whole function is convex.
+--  
+--  Be that as it may, if you know what the your lagrangian is convex you can use 'solve' to 
+--  find the minimum.
+--  
+--  For example, find the maximum entropy with the constraint that the probabilities add
+--  up to one. 
+--  
+--  @ 
+--     solve (negate . sum . map (\x -> x * log x), [(sum, 1)]) 3
+--  @
+--  
+--  Gives the answer ([0.33, 0.33, 0.33], [-0.09])
+--  
+--  The first elements of the result pair are the arguments for the objective function at the minimum. 
+--  The second elements are the lagrange multipliers.
+module Numeric.AD.Lagrangian (
+    solve) where
+import Numeric.AD.Lagrangian.Internal (solve, feasible)
diff --git a/src/Numeric/AD/Lagrangian/Internal.hs b/src/Numeric/AD/Lagrangian/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Lagrangian/Internal.hs
@@ -0,0 +1,105 @@
+{-# LANGUAGE Rank2Types #-}
+module Numeric.AD.Lagrangian.Internal where
+import Numeric.Optimization.Algorithms.HagerZhang05
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Storable as S
+import Numeric.AD
+import GHC.IO                   (unsafePerformIO)
+import Numeric.AD.Types
+import Numeric.AD.Internal.Classes
+import Numeric.LinearAlgebra.Algorithms
+import qualified Data.Packed.Vector as V
+import qualified Data.Packed.Matrix as M
+
+-- In general I am fighting against the lack of type inference rank two types.
+-- Hopefully some of the explicit type signatures can be removed.
+
+
+-- The type for the contraints.
+-- Given a constraint g(x, y, ...) = c, we would represent it as (g, c).
+type Constraint a = ([a] -> a, a)
+
+-- | This is not a true feasibility test for the function. I am not sure exactly how to 
+--   implement that. This just checks the feasiblility at point. If this ever returns 
+--   false, 'solve' can fail.
+feasible :: (forall a. Floating a => ([a] -> a, [Constraint a], [a]))
+      -> Bool
+feasible params = result where
+    obj :: Floating a => [a] -> a
+    obj argsAndLams = squaredGrad lang argsAndLams
+
+    lang :: Floating a => (forall s. Mode s => [AD s a] -> AD s a)
+    lang = lagrangian fAndGs (length point)
+    
+    fAndGs :: (forall a. Floating a => ([a] -> a, [Constraint a]))
+    fAndGs = (\(x, y, _) -> (x, y)) params
+    
+    point :: Floating a => [a]
+    point = (\(_, _, x) -> x) params
+    
+    h :: [[Double]]
+    h = hessian obj point
+    -- I want the hessian as a matrix
+    hessianMatrix = M.fromLists h
+
+    -- make sure all of the eigenvalues are positive
+    result = all (>0) . V.toList . eigenvaluesSH $ hessianMatrix 
+
+-- | This is the lagrangrain multiplier solver. It is assumed that the 
+--   objective function and all of the constraints take in the 
+--   same about of arguments.
+solve :: (forall a. Floating a => ([a] -> a, [Constraint a])) -- ^ A pair of the function to minimize and the constraints
+      -> Int -- ^ The arity of the objective function and the constraints.
+      -> Either (Result, Statistics) ([Double], [Double]) -- ^ Either an explaination of why the gradient descent failed or a pair of the arguments at the minimum and the lagrange multipliers
+solve params argCount = result where
+    obj :: Floating a => [a] -> a
+    obj argsAndLams = squaredGrad lang argsAndLams
+
+    lang :: Floating a => (forall s. Mode s => [AD s a] -> AD s a)
+    lang = lagrangian params argCount
+    
+    constraintCount = length (snd params)
+    
+    guess = U.fromList $ replicate (argCount + constraintCount) (1.0 :: Double) 
+
+    result = case unsafePerformIO (optimize (defaultParameters { printFinal = False }) 
+                    0.00001 guess (toFunction obj) (toGradient obj)
+                       Nothing) of
+        
+       (vs, ToleranceStatisfied, _) -> Right (take argCount . S.toList $ vs, 
+                                              drop argCount . S.toList $ vs) 
+       (_, x, y) -> Left (x, y)
+
+-- Convert a objective function and a list of constraints to a lagrangian
+lagrangian :: Floating a
+             => ([a] -> a, [Constraint a]) 
+             -> Int
+             -> [a] 
+             -> a
+lagrangian (f, constraints) argsLength argsAndLams = result where
+    -- L(x, y, ..., lam0, lam1, ...) = f(x, y, ...) + 
+    result = f args + (sum $ zipWith (*) lams appliedConstraints)
+    
+    -- Apply the arguments to the constraint function
+    -- and subtract to set equal to zero
+    -- (g, c) <=> g(x, y, ...) = c <=> g(x, y, ...) - c = 0
+    appliedConstraints = map (\(f, c) -> f args - c) constraints
+
+    -- Split the input by args and lambdas.
+    -- It is assumed that the args for f and g's come before the
+    -- lambdas for the constraints
+    args = take argsLength argsAndLams
+    lams = drop argsLength argsAndLams
+
+sumMap f = sum . map f 
+
+squaredGrad :: Num a 
+            => (forall s. Mode s => [AD s a] -> AD s a) -> [a] -> a
+squaredGrad f vs = sumMap (\x -> x*x) (grad f vs)
+
+toFunction :: (forall a. Floating a => [a] -> a) -> Function Simple
+toFunction f = VFunction (f . U.toList)
+
+toGradient :: (forall a. Floating a => [a] -> a) -> Gradient Simple
+toGradient f = VGradient (U.fromList . grad f . U.toList)
+
diff --git a/tests/Main.hs b/tests/Main.hs
new file mode 100644
--- /dev/null
+++ b/tests/Main.hs
@@ -0,0 +1,35 @@
+module Main where
+import Test.Framework (defaultMain, testGroup, defaultMainWithArgs)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.Framework.Providers.QuickCheck2 (testProperty)
+import Numeric.AD.Lagrangian.Internal
+import Control.Applicative
+
+main = defaultMain [
+        testGroup "trival test" [
+            testCase "noConstraints" noConstraints,
+            testCase "entropyTest" entropyTest
+        ]
+    ] 
+    
+    
+noConstraints = (fst <$> actual) @?= Right expected where
+    actual    = solve (f, []) 1
+    expected  = [1]
+    f [x] = -(x - 1) ^2
+    
+--class Approximate a where
+--    x =~= y :: a -> a -> Bool
+
+
+entropyTest = (sum . map abs $ zipWith (-) actual expected) < 0.02 @?= True  where
+    Right actual = fst <$> solve (f, [(\xs -> sum xs, 1.0)]) 3
+    expected  = [0.33, 0.33, 0.33]
+    f :: Floating a => [a] -> a
+    f = negate . sum . map (\x -> x * log x)
+    
+    
+
+    
+