<commit_before>{-# LANGUAGE FlexibleInstances #-}

-- | Exploring optimisation as presented in "An introduciton to the
-- conjugate gradient method without the agonizing pain"

module CG1 where

import Prelude hiding (id)

-- | S for Scalar
type S = Double

-- | for simplicity only 2 by 2 matrices for now
type M2 = (S,S
          ,S,S)

-- | vectors of dimension 2
type V2 = (S,S)

instance Num a => Num (a,a) where
  fromInteger x = (fromInteger x,fromInteger x)
  (a,b) + (x,y) = (a + x,b + y)
  (a,b) * (x,y) = (a * x,b * y)
  abs (a,b) = (abs a,abs b)
  signum = undefined
  negate (a,b) = (-a,-b)

instance Fractional a => Fractional (a,a,a,a) where
  fromRational a = (fromRational a,fromRational a,fromRational a,fromRational a)
  recip (a,b,c,d) = (recip a,recip b,recip c,recip d)

instance Num a => Num (a,a,a,a)  where
  fromInteger x = (fromInteger x,fromInteger x
                  ,fromInteger x,fromInteger x)
  (a,b,c,d) + (x,y,z,w) = (a+x,b+y
                          ,c+z,d+w)
  (a,b,c,d) * (x,y,z,w) = (a*b+b*z,a*y+b*w
                          ,c*x+d*z,c*y+d*w) -- matrix multiplication
  abs (a,b,c,d) = (abs a,abs b,abs c,abs d)
  signum = undefined
  negate (a,b,c,d) = (negate a,negate b,negate c,negate d)

instance Fractional a => Fractional (a,a) where
  fromRational a = (fromRational a,fromRational a)
  recip (a,b) = (recip a,recip b)

id :: M2
id = (1,0
     ,0,1)

innerMul :: V2 -> V2 -> S
innerMul (a,b) (c,d) = a*c + b*d

preMul :: V2 -> M2 -> V2
preMul (a,b) (x,y
             ,z,w) = (a * x + b * z
                     ,a * y + b * w)

postMul :: M2 -> V2 -> V2
postMul (x,y
        ,z,w) (a,b) = (x * a + y * b
                      ,z * a + w * b)

-- | v^T A v = ppm v A
--
prePostMul :: V2 -> M2 -> S
prePostMul v a = innerMul (preMul v a) v

-- matMul :: M2 -> M2 -> M2
-- matMul (a,b
--        ,c,d) (x,y
--              ,z,w) = (a*b+b*z,a*y+b*w
--                      ,c*x+d*z,c*y+d*w)

transpose :: M2 -> M2
transpose (x,y
          ,z,w) = (x,z
                  ,y,w)

scaleM :: S -> M2 -> M2
scaleM x (a,b
         ,c,d) = (x * a,x * b
                 ,x * c,x * d)

scaleV :: S -> V2 -> V2
scaleV x (a,b) = (x*a,x*b)

det :: M2 -> S
det (a,b
    ,c,d) = a*d - b*c

-- a matrix m is postive definite if for all x . x^Tmx > 0
prop_positive_definite :: M2 -> V2 -> Bool
prop_positive_definite m = \v -> prePostMul v m > 0

-- | Example vect and matrix
a :: M2
a = (3,2
    ,2,6)

b :: V2
b = (2,8)

c :: S
c = 0

-- quadratic form
f x = 0.5 * prePostMul x a - innerMul b x + c

-- gradient of quadradic form
f' x = 0.5 * postMul (transpose a) x + 0.5 * postMul a x - b

-- error, how far from goal are we
-- e_i = x_i - x
-- residual, how far from the correct value of b
-- r_i = b - Ax_i = -Ae_i
-- "the error transformed by A into the same space as b"


-- steepest descent takes a step
-- x_1 = x_0 + alpha*r_0
-- how big should alpha be? p. 6
sdStep x = let r = b - postMul a x -- r = b - Ax
               alpha = innerMul r r / prePostMul r a -- a = (r^T r)/(r^T A r)
               x' = x + scaleV alpha r -- x' = x + a r
           in x'


-- Jacobi iterations, p 11

-- a matrix is diagonal if all but the diagonal is 0
prop_diagonal :: M2 -> Bool
prop_diagonal (_,b
              ,c,_) = b == 0 && c == 0 -- possibly with some epsilon

-- invert a diagonal matrix
invertDiag :: M2 -> M2
invertDiag (a,_
           ,_,d) = (recip a,0
                   ,0      ,recip d)

-- split matrix into two parts: D whose diagonal elemts are equal to
-- A's and off diagonal elements are 0 and E whose diagonal
-- elements are 0 and off diagonal elements equal to A's
jacobiSplit (a,b
            ,c,d) = ((a,0
                     ,0,d)
                    ,(0,b
                     ,c,0))

jacobiStep x =
  let (d,e) = jacobiSplit a
      f = (negate (invertDiag d)) * e -- can precompute this
      z = postMul (invertDiag d) b
      x' = postMul f x + z
  in x'<commit_msg>prop: A orthogonal vectors<commit_after>{-# LANGUAGE FlexibleInstances #-}

-- | Exploring optimisation as presented in "An introduciton to the
-- conjugate gradient method without the agonizing pain"

module CG1 where

import Prelude hiding (id)

-- | S for Scalar
type S = Double

-- | for simplicity only 2 by 2 matrices for now
type M2 = (S,S
          ,S,S)

-- | vectors of dimension 2
type V2 = (S,S)

instance Num a => Num (a,a) where
  fromInteger x = (fromInteger x,fromInteger x)
  (a,b) + (x,y) = (a + x,b + y)
  (a,b) * (x,y) = (a * x,b * y)
  abs (a,b) = (abs a,abs b)
  signum = undefined
  negate (a,b) = (-a,-b)

instance Fractional a => Fractional (a,a,a,a) where
  fromRational a = (fromRational a,fromRational a,fromRational a,fromRational a)
  recip (a,b,c,d) = (recip a,recip b,recip c,recip d)

instance Num a => Num (a,a,a,a)  where
  fromInteger x = (fromInteger x,fromInteger x
                  ,fromInteger x,fromInteger x)
  (a,b,c,d) + (x,y,z,w) = (a+x,b+y
                          ,c+z,d+w)
  (a,b,c,d) * (x,y,z,w) = (a*b+b*z,a*y+b*w
                          ,c*x+d*z,c*y+d*w) -- matrix multiplication
  abs (a,b,c,d) = (abs a,abs b,abs c,abs d)
  signum = undefined
  negate (a,b,c,d) = (negate a,negate b,negate c,negate d)

instance Fractional a => Fractional (a,a) where
  fromRational a = (fromRational a,fromRational a)
  recip (a,b) = (recip a,recip b)

id :: M2
id = (1,0
     ,0,1)

innerMul :: V2 -> V2 -> S
innerMul (a,b) (c,d) = a*c + b*d

preMul :: V2 -> M2 -> V2
preMul (a,b) (x,y
             ,z,w) = (a * x + b * z
                     ,a * y + b * w)

postMul :: M2 -> V2 -> V2
postMul (x,y
        ,z,w) (a,b) = (x * a + y * b
                      ,z * a + w * b)

-- | v^T A v = ppm v A
--
prePostMul :: V2 -> M2 -> S
prePostMul v a = innerMul (preMul v a) v

-- matMul :: M2 -> M2 -> M2
-- matMul (a,b
--        ,c,d) (x,y
--              ,z,w) = (a*b+b*z,a*y+b*w
--                      ,c*x+d*z,c*y+d*w)

transpose :: M2 -> M2
transpose (x,y
          ,z,w) = (x,z
                  ,y,w)

scaleM :: S -> M2 -> M2
scaleM x (a,b
         ,c,d) = (x * a,x * b
                 ,x * c,x * d)

scaleV :: S -> V2 -> V2
scaleV x (a,b) = (x*a,x*b)

det :: M2 -> S
det (a,b
    ,c,d) = a*d - b*c

-- a matrix m is postive definite if for all x . x^Tmx > 0
prop_positive_definite :: M2 -> V2 -> Bool
prop_positive_definite m = \v -> prePostMul v m > 0

-- | Example vect and matrix
a :: M2
a = (3,2
    ,2,6)

b :: V2
b = (2,8)

c :: S
c = 0

-- quadratic form
f x = 0.5 * prePostMul x a - innerMul b x + c

-- gradient of quadradic form
f' x = 0.5 * postMul (transpose a) x + 0.5 * postMul a x - b

-- error, how far from goal are we
-- e_i = x_i - x
-- residual, how far from the correct value of b
-- r_i = b - Ax_i = -Ae_i
-- "the error transformed by A into the same space as b"


-- steepest descent takes a step
-- x_1 = x_0 + alpha*r_0
-- how big should alpha be? p. 6
sdStep x = let r = b - postMul a x -- r = b - Ax
               alpha = innerMul r r / prePostMul r a -- a = (r^T r)/(r^T A r)
               x' = x + scaleV alpha r -- x' = x + a r
           in x'


-- Jacobi iterations, p 11

-- a matrix is diagonal if all but the diagonal is 0
prop_diagonal :: M2 -> Bool
prop_diagonal (_,b
              ,c,_) = b == 0 && c == 0 -- possibly with some epsilon

-- invert a diagonal matrix
invertDiag :: M2 -> M2
invertDiag (a,_
           ,_,d) = (recip a,0
                   ,0      ,recip d)

-- split matrix into two parts: D whose diagonal elemts are equal to
-- A's and off diagonal elements are 0 and E whose diagonal
-- elements are 0 and off diagonal elements equal to A's
jacobiSplit (a,b
            ,c,d) = ((a,0
                     ,0,d)
                    ,(0,b
                     ,c,0))

jacobiStep x =
  let (d,e) = jacobiSplit a
      f = (negate (invertDiag d)) * e -- can precompute this
      z = postMul (invertDiag d) b
      x' = postMul f x + z
  in x'


-----

-- | v1^T A v2 = 0 <=> v1, v2 are A orthogonal
prop_A_orthogonal :: M2 -> V2 -> V2 -> Bool
prop_A_orthogonal a v1 v2 = innerMul (preMul v1 a) v2 == 0
{'fraction_non_alphanumeric': 0.1569115815691158, 'fraction_numerical': 0.019053549190535493, 'mean_word_length': 1.9417582417582417, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7204104', 'n_tokens_mistral': 3451, 'n_tokens_neox': 3316, 'n_words': 1092}