```{-# LANGUAGE TupleSections, Rank2Types #-}
module Numeric.MaxEnt.Linear where
import Numeric.Optimization.Algorithms.HagerZhang05
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Storable as S
import GHC.IO                   (unsafePerformIO)
import Data.Traversable
import Data.List (transpose)
import Control.Applicative

probs :: (RealFloat b) => [[b]] -> [b] -> [b]
probs matrix ls = map (pOfK matrix ls) [0..length matrix - 1]

pOfK :: RealFloat a => [[a]] -> [a] -> Int -> a
pOfK matrix ls k = exp (dot (transpose matrix !! k) ls) / partitionFunc matrix ls

-- have the column and rows backwards
partitionFunc :: RealFloat a => [[a]] -> [a] -> a
partitionFunc matrix ws = sum \$ [ exp (dot as ws) | as <- transpose matrix]

-- This is almost the sam as the objectiveFunc
objectiveFunc :: RealFloat a => [[a]] -> [a]-> [a] -> a
objectiveFunc as moments ls = log (partitionFunc as ls) - dot ls moments

data LinearConstraints a = LC {
matrix :: [[a]],
output :: [a]
}
deriving (Show, Eq)

-- | This is for the linear case Ax = b
--   @x@ in this situation is the vector of probablities.
--
--  Consider the 1 dimensional circular convolution using hmatrix.
--
--  >>> import Numeric.LinearAlgebra
--  >>> fromLists [[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]] <> fromLists [[0.2], [0.5], [0.3]]
--  (3><1) [0.276, 0.426, 0.298]
--
--   Now if we were given just the convolution and the output, we can use 'linear' to infer the input.
--
--   >>> linear 3.0e-17 \$ LC [[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]] [0.276, 0.426, 0.298]
--   Right [0.20000000000000004,0.4999999999999999,0.3]
--
--
--
--
linear :: Double -- ^ Tolerance for the numerical solver
-> (forall a. RealFloat a => LinearConstraints a) -- ^ The matrix A and column vector b
-> Either (Result, Statistics) [Double] -- ^ Either the a discription of what wrong or the probability distribution
linear tolerance constraints = result where
obj :: RealFloat a => [a] -> a
obj = objectiveFunc (matrix constraints) (output constraints)

as :: [[Double]]
as = matrix constraints

result = probs as <\$> solve tolerance (length as) obj

--Just for testing
linear1 :: Double -> (forall a. RealFloat a => LinearConstraints a) -- ^ a matrix A and column vector b
-> Either (Result, Statistics) [Double] -- ^ Either the a discription of what wrong or the probability distribution
linear1 precision constraints = result where
obj :: RealFloat a => [a] -> a
obj = objectiveFunc (matrix constraints) (output constraints)

fs :: [[Double]]
fs = matrix constraints

guess :: RealFloat a => [a]
guess = replicate
(length fs) ((1.0) / (fromIntegral (length fs)))

result = Right . probs fs . last \$ conjugateGradientDescent obj guess

{-
I need to clean this up and submit it
submit the sharpen

I need to figure out how to make linear2
what is the
it might help to clean up the notation for the problem
basically I need to symbolically take the derative of the problem
and get rid of the automatic differentation
which I can do hand

-}

```