{-# LANGUAGE TypeOperators #-}
{- |
Do not import this module. It is only for demonstration purposes.
-}
module Numeric.LAPACK.Example.EconomicAllocation where

import qualified Numeric.LAPACK.Matrix.Square as Square
import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
import qualified Numeric.LAPACK.Matrix as Matrix
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Matrix.Square ((|=|))
import Numeric.LAPACK.Matrix (ShapeInt, shapeInt, (#-#), (#*|), (#\|), (\\#))
import Numeric.LAPACK.Vector ((|-|))
import Numeric.LAPACK.Format ((##))

import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Shape ((::+)((::+)))

import qualified Data.Stream as Stream


{- $setup
>>> import Numeric.LAPACK.Example.EconomicAllocation
>>> import Test.Utility (approxVector)
>>>
>>> import qualified Numeric.LAPACK.Vector as Vector
>>> import Numeric.LAPACK.Vector ((+++))
>>>
>>> import qualified Data.Array.Comfort.Storable as Array
-}


type ZeroInt2 = ShapeInt::+ShapeInt
type Vector sh = Vector.Vector sh Double
type Matrix height width = Matrix.General height width Double
type SquareMatrix size = Square.Square size Double


balances0 :: Vector ZeroInt2
balances0 :: Vector ZeroInt2
balances0 =
   ZeroInt2 -> [Double] -> Vector ZeroInt2
forall sh a. (C sh, Storable a) => sh -> [a] -> Vector sh a
Vector.fromList (Int -> ShapeInt
shapeInt Int
2 ShapeInt -> ShapeInt -> ZeroInt2
forall sh0 sh1. sh0 -> sh1 -> sh0 ::+ sh1
::+ Int -> ShapeInt
shapeInt Int
2)
      [Double
100000, Double
90000, -Double
50000, -Double
120000]

expenses0 :: Matrix ShapeInt ZeroInt2
expenses0 :: Matrix ShapeInt ZeroInt2
expenses0 =
   ShapeInt -> ZeroInt2 -> [Double] -> Matrix ShapeInt ZeroInt2
forall height width a.
(C height, C width, Storable a) =>
height -> width -> [a] -> General height width a
Matrix.fromList (Int -> ShapeInt
shapeInt Int
2) (Int -> ShapeInt
shapeInt Int
2 ShapeInt -> ShapeInt -> ZeroInt2
forall sh0 sh1. sh0 -> sh1 -> sh0 ::+ sh1
::+ Int -> ShapeInt
shapeInt Int
2) ([Double] -> Matrix ShapeInt ZeroInt2)
-> [Double] -> Matrix ShapeInt ZeroInt2
forall a b. (a -> b) -> a -> b
$
   [Double
16000,  Double
4000,  Double
8000, Double
12000,
    Double
10000, Double
30000, Double
40000, Double
20000]

normalize ::
   (Eq height, Shape.C height, Shape.C width) =>
   Matrix height width -> Matrix height width
normalize :: forall height width.
(Eq height, C height, C width) =>
Matrix height width -> Matrix height width
normalize Matrix height width
x = Matrix height width -> Vector height Double
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a -> Vector height a
Matrix.rowSums Matrix height width
x Vector height Double -> Matrix height width -> Matrix height width
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 Floating a) =>
Vector height a
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
\\# Matrix height width
x

normalizeSplit ::
   (Shape.C sh0, Shape.C sh1, Eq sh1) =>
   Matrix sh1 (sh0::+sh1) -> (Matrix sh0 sh1, SquareMatrix sh1)
normalizeSplit :: forall sh0 sh1.
(C sh0, C sh1, Eq sh1) =>
Matrix sh1 (sh0 ::+ sh1) -> (Matrix sh0 sh1, SquareMatrix sh1)
normalizeSplit Matrix sh1 (sh0 ::+ sh1)
expenses =
   let a :: Matrix
  (Array Unpacked Arbitrary)
  ()
  ()
  Filled
  Filled
  Size
  Big
  Big
  (sh0 ::+ sh1)
  sh1
  Double
a = Matrix sh1 (sh0 ::+ sh1)
-> Matrix
     (Array Unpacked Arbitrary)
     ()
     ()
     Filled
     Filled
     Size
     Big
     Big
     (sh0 ::+ sh1)
     sh1
     Double
forall typ xl xu meas vert horiz height width a lower upper.
(Transpose typ, TransposeExtra typ xl, TransposeExtra typ xu,
 Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Matrix typ xl xu lower upper meas vert horiz height width a
-> Matrix typ xu xl upper lower meas horiz vert width height a
forall xl xu meas vert horiz height width a lower upper.
(TransposeExtra (Array Unpacked Arbitrary) xl,
 TransposeExtra (Array Unpacked Arbitrary) xu, Measure meas, C vert,
 C horiz, C height, C width, Floating a) =>
Matrix
  (Array Unpacked Arbitrary)
  xl
  xu
  lower
  upper
  meas
  vert
  horiz
  height
  width
  a
-> Matrix
     (Array Unpacked Arbitrary)
     xu
     xl
     upper
     lower
     meas
     horiz
     vert
     width
     height
     a
Matrix.transpose (Matrix sh1 (sh0 ::+ sh1)
 -> Matrix
      (Array Unpacked Arbitrary)
      ()
      ()
      Filled
      Filled
      Size
      Big
      Big
      (sh0 ::+ sh1)
      sh1
      Double)
-> Matrix sh1 (sh0 ::+ sh1)
-> Matrix
     (Array Unpacked Arbitrary)
     ()
     ()
     Filled
     Filled
     Size
     Big
     Big
     (sh0 ::+ sh1)
     sh1
     Double
forall a b. (a -> b) -> a -> b
$ Matrix sh1 (sh0 ::+ sh1) -> Matrix sh1 (sh0 ::+ sh1)
forall height width.
(Eq height, C height, C width) =>
Matrix height width -> Matrix height width
normalize Matrix sh1 (sh0 ::+ sh1)
expenses
   in (Matrix
  (Array Unpacked Arbitrary)
  ()
  ()
  Filled
  Filled
  Size
  Big
  Big
  (sh0 ::+ sh1)
  sh1
  Double
-> Full Size Big Big sh0 sh1 Double
forall vert height0 height1 width a.
(C vert, C height0, C height1, C width, Floating a) =>
Full Size vert Big (height0 ::+ height1) width a
-> Full Size vert Big height0 width a
Matrix.takeTop Matrix
  (Array Unpacked Arbitrary)
  ()
  ()
  Filled
  Filled
  Size
  Big
  Big
  (sh0 ::+ sh1)
  sh1
  Double
a, Full Size Big Big sh1 sh1 Double -> Square sh1 Double
forall meas vert horiz sh a.
(Measure meas, C vert, C horiz, Eq sh) =>
Full meas vert horiz sh sh a -> Square sh a
Square.fromFull (Full Size Big Big sh1 sh1 Double -> Square sh1 Double)
-> Full Size Big Big sh1 sh1 Double -> Square sh1 Double
forall a b. (a -> b) -> a -> b
$ Matrix
  (Array Unpacked Arbitrary)
  ()
  ()
  Filled
  Filled
  Size
  Big
  Big
  (sh0 ::+ sh1)
  sh1
  Double
-> Full Size Big Big sh1 sh1 Double
forall vert height0 height1 width a.
(C vert, C height0, C height1, C width, Floating a) =>
Full Size vert Big (height0 ::+ height1) width a
-> Full Size vert Big height1 width a
Matrix.takeBottom Matrix
  (Array Unpacked Arbitrary)
  ()
  ()
  Filled
  Filled
  Size
  Big
  Big
  (sh0 ::+ sh1)
  sh1
  Double
a)


completeIdSquare ::
   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>
   Matrix sh1 (sh0::+sh1) -> SquareMatrix (sh0::+sh1)
completeIdSquare :: forall sh0 sh1.
(C sh0, Eq sh0, C sh1, Eq sh1) =>
Matrix sh1 (sh0 ::+ sh1) -> SquareMatrix (sh0 ::+ sh1)
completeIdSquare Matrix sh1 (sh0 ::+ sh1)
x =
   let (Matrix sh0 sh1
p,SquareMatrix sh1
k) = Matrix sh1 (sh0 ::+ sh1) -> (Matrix sh0 sh1, SquareMatrix sh1)
forall sh0 sh1.
(C sh0, C sh1, Eq sh1) =>
Matrix sh1 (sh0 ::+ sh1) -> (Matrix sh0 sh1, SquareMatrix sh1)
normalizeSplit Matrix sh1 (sh0 ::+ sh1)
x
   in (Matrix sh0 sh1 -> Square sh0 Double
forall height width a.
(C height, C width, Floating a) =>
General height width a -> Square height a
Square.identityFromHeight Matrix sh0 sh1
p, Matrix sh0 sh1
p)
      (Square sh0 Double, Matrix sh0 sh1)
-> (Full Size Big Big sh1 sh0 Double, SquareMatrix sh1)
-> Square (sh0 ::+ sh1) Double
forall meas vert horiz sizeA sizeB a.
(Measure meas, C vert, C horiz, C sizeA, Eq sizeA, C sizeB,
 Eq sizeB, Floating a) =>
(Square sizeA a, Full meas vert horiz sizeA sizeB a)
-> (Full meas horiz vert sizeB sizeA a, Square sizeB a)
-> Square (sizeA ::+ sizeB) a
|=|
      (Omni Unpacked Arbitrary Filled Filled Size Big Big sh1 sh0
-> Full Size Big Big sh1 sh0 Double
forall property meas vert horiz height width a pack lower upper.
(Homogeneous property, Measure meas, C vert, C horiz, C height,
 C width, Floating a) =>
Omni pack property lower upper meas vert horiz height width
-> ArrayMatrix
     pack property lower upper meas vert horiz height width a
Matrix.zero (Omni Unpacked Arbitrary Filled Filled Size Big Big sh1 sh0
 -> Full Size Big Big sh1 sh0 Double)
-> Omni Unpacked Arbitrary Filled Filled Size Big Big sh1 sh0
-> Full Size Big Big sh1 sh0 Double
forall a b. (a -> b) -> a -> b
$ Full Size Big Big sh1 sh0 Double
-> Omni Unpacked Arbitrary Filled Filled Size Big Big sh1 sh0
forall pack property lower upper meas vert horiz height width a.
ArrayMatrix
  pack property lower upper meas vert horiz height width a
-> Omni pack property lower upper meas vert horiz height width
ArrMatrix.shape (Full Size Big Big sh1 sh0 Double
 -> Omni Unpacked Arbitrary Filled Filled Size Big Big sh1 sh0)
-> Full Size Big Big sh1 sh0 Double
-> Omni Unpacked Arbitrary Filled Filled Size Big Big sh1 sh0
forall a b. (a -> b) -> a -> b
$ Matrix sh0 sh1 -> Full Size Big Big sh1 sh0 Double
forall typ xl xu meas vert horiz height width a lower upper.
(Transpose typ, TransposeExtra typ xl, TransposeExtra typ xu,
 Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Matrix typ xl xu lower upper meas vert horiz height width a
-> Matrix typ xu xl upper lower meas horiz vert width height a
forall xl xu meas vert horiz height width a lower upper.
(TransposeExtra (Array Unpacked Arbitrary) xl,
 TransposeExtra (Array Unpacked Arbitrary) xu, Measure meas, C vert,
 C horiz, C height, C width, Floating a) =>
Matrix
  (Array Unpacked Arbitrary)
  xl
  xu
  lower
  upper
  meas
  vert
  horiz
  height
  width
  a
-> Matrix
     (Array Unpacked Arbitrary)
     xu
     xl
     upper
     lower
     meas
     horiz
     vert
     width
     height
     a
Matrix.transpose Matrix sh0 sh1
p, SquareMatrix sh1
k)

iterated ::
   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>
   Matrix sh1 (sh0::+sh1) -> Vector (sh0::+sh1) -> Vector (sh0::+sh1)
iterated :: forall sh0 sh1.
(C sh0, Eq sh0, C sh1, Eq sh1) =>
Matrix sh1 (sh0 ::+ sh1)
-> Vector (sh0 ::+ sh1) -> Vector (sh0 ::+ sh1)
iterated Matrix sh1 (sh0 ::+ sh1)
expenses =
   Stream (Vector (sh0 ::+ sh1)) -> Vector (sh0 ::+ sh1)
forall a. Stream a -> a
Stream.head (Stream (Vector (sh0 ::+ sh1)) -> Vector (sh0 ::+ sh1))
-> (Vector (sh0 ::+ sh1) -> Stream (Vector (sh0 ::+ sh1)))
-> Vector (sh0 ::+ sh1)
-> Vector (sh0 ::+ sh1)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (Vector (sh0 ::+ sh1) -> Bool)
-> Stream (Vector (sh0 ::+ sh1)) -> Stream (Vector (sh0 ::+ sh1))
forall a. (a -> Bool) -> Stream a -> Stream a
Stream.dropWhile ((Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>=Double
1e-5) (Double -> Bool)
-> (Vector (sh0 ::+ sh1) -> Double) -> Vector (sh0 ::+ sh1) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector sh1 Double -> Double
Vector sh1 Double -> RealOf Double
forall sh a. (C sh, Floating a) => Vector sh a -> RealOf a
Vector.normInf (Vector sh1 Double -> Double)
-> (Vector (sh0 ::+ sh1) -> Vector sh1 Double)
-> Vector (sh0 ::+ sh1)
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (sh0 ::+ sh1) -> Vector sh1 Double
forall sh0 sh1 a.
(C sh0, C sh1, Storable a) =>
Array (sh0 ::+ sh1) a -> Array sh1 a
Vector.takeRight) (Stream (Vector (sh0 ::+ sh1)) -> Stream (Vector (sh0 ::+ sh1)))
-> (Vector (sh0 ::+ sh1) -> Stream (Vector (sh0 ::+ sh1)))
-> Vector (sh0 ::+ sh1)
-> Stream (Vector (sh0 ::+ sh1))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (Vector (sh0 ::+ sh1) -> Vector (sh0 ::+ sh1))
-> Vector (sh0 ::+ sh1) -> Stream (Vector (sh0 ::+ sh1))
forall a. (a -> a) -> a -> Stream a
Stream.iterate (Matrix sh1 (sh0 ::+ sh1) -> SquareMatrix (sh0 ::+ sh1)
forall sh0 sh1.
(C sh0, Eq sh0, C sh1, Eq sh1) =>
Matrix sh1 (sh0 ::+ sh1) -> SquareMatrix (sh0 ::+ sh1)
completeIdSquare Matrix sh1 (sh0 ::+ sh1)
expenses SquareMatrix (sh0 ::+ sh1)
-> Vector (sh0 ::+ sh1) -> Vector (sh0 ::+ sh1)
forall typ lower upper xl xu meas vert horiz height width a.
(MultiplyVector typ, Strip lower, Strip upper,
 MultiplyVectorExtra typ xl, MultiplyVectorExtra typ xu,
 Measure meas, C vert, C horiz, C height, C width, Eq width,
 Floating a) =>
Matrix typ xl xu lower upper meas vert horiz height width a
-> Vector width a -> Vector height a
#*|)


compensated ::
   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>
   Matrix sh1 (sh0::+sh1) -> Vector (sh0::+sh1) -> Vector sh0
compensated :: forall sh0 sh1.
(C sh0, Eq sh0, C sh1, Eq sh1) =>
Matrix sh1 (sh0 ::+ sh1) -> Vector (sh0 ::+ sh1) -> Vector sh0
compensated Matrix sh1 (sh0 ::+ sh1)
expenses Vector (sh0 ::+ sh1)
balances =
   let (Matrix sh0 sh1
p,SquareMatrix sh1
k) = Matrix sh1 (sh0 ::+ sh1) -> (Matrix sh0 sh1, SquareMatrix sh1)
forall sh0 sh1.
(C sh0, C sh1, Eq sh1) =>
Matrix sh1 (sh0 ::+ sh1) -> (Matrix sh0 sh1, SquareMatrix sh1)
normalizeSplit Matrix sh1 (sh0 ::+ sh1)
expenses
       x :: Array sh0 Double
x = Vector (sh0 ::+ sh1) -> Array sh0 Double
forall sh0 sh1 a.
(C sh0, C sh1, Storable a) =>
Array (sh0 ::+ sh1) a -> Array sh0 a
Vector.takeLeft Vector (sh0 ::+ sh1)
balances
       y :: Array sh1 Double
y = Vector (sh0 ::+ sh1) -> Array sh1 Double
forall sh0 sh1 a.
(C sh0, C sh1, Storable a) =>
Array (sh0 ::+ sh1) a -> Array sh1 a
Vector.takeRight Vector (sh0 ::+ sh1)
balances
   in Array sh0 Double
x Array sh0 Double -> Array sh0 Double -> Array sh0 Double
forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
|-| Matrix sh0 sh1
p Matrix sh0 sh1 -> Array sh1 Double -> Array sh0 Double
forall typ lower upper xl xu meas vert horiz height width a.
(MultiplyVector typ, Strip lower, Strip upper,
 MultiplyVectorExtra typ xl, MultiplyVectorExtra typ xu,
 Measure meas, C vert, C horiz, C height, C width, Eq width,
 Floating a) =>
Matrix typ xl xu lower upper meas vert horiz height width a
-> Vector width a -> Vector height a
#*| (SquareMatrix sh1
k SquareMatrix sh1 -> SquareMatrix sh1 -> SquareMatrix sh1
forall meas vert horiz typ xl xu height width a lower upper.
(Measure meas, C vert, C horiz, Subtractive typ,
 SubtractiveExtra typ xl, SubtractiveExtra typ xu, C height,
 Eq height, C width, Eq width, Floating a) =>
Matrix typ xl xu lower upper meas vert horiz height width a
-> Matrix typ xl xu lower upper meas vert horiz height width a
-> Matrix typ xl xu lower upper meas vert horiz height width a
#-# SquareMatrix sh1 -> SquareMatrix sh1
forall sh a. (C sh, Floating a) => Square sh a -> Square sh a
Square.identityFrom SquareMatrix sh1
k) SquareMatrix sh1 -> Array sh1 Double -> Array sh1 Double
forall typ xl xu lower upper meas height width a.
(Solve typ, ToQuadratic typ, SolveExtra typ xl, SolveExtra typ xu,
 BoxExtra typ xl, BoxExtra typ xu, Strip lower, Strip upper,
 Measure meas, C height, C width, Eq height, Floating a) =>
QuadraticMeas typ xl xu lower upper meas height width a
-> Vector height a -> Vector width a
#\| Array sh1 Double
y


{- |
prop> let result = iterated expenses0 balances0 in approxVector result $ compensated expenses0 balances0 +++ Vector.zero (Array.shape $ Vector.takeRight result)
-}
main :: IO ()
main :: IO ()
main = do
   Matrix ShapeInt ZeroInt2 -> Vector ZeroInt2 -> Vector ZeroInt2
forall sh0 sh1.
(C sh0, Eq sh0, C sh1, Eq sh1) =>
Matrix sh1 (sh0 ::+ sh1)
-> Vector (sh0 ::+ sh1) -> Vector (sh0 ::+ sh1)
iterated Matrix ShapeInt ZeroInt2
expenses0 Vector ZeroInt2
balances0 Vector ZeroInt2 -> String -> IO ()
forall a. Format a => a -> String -> IO ()
## String
"%10.2f"
   Matrix ShapeInt ZeroInt2 -> Vector ZeroInt2 -> Vector ShapeInt
forall sh0 sh1.
(C sh0, Eq sh0, C sh1, Eq sh1) =>
Matrix sh1 (sh0 ::+ sh1) -> Vector (sh0 ::+ sh1) -> Vector sh0
compensated Matrix ShapeInt ZeroInt2
expenses0 Vector ZeroInt2
balances0 Vector ShapeInt -> String -> IO ()
forall a. Format a => a -> String -> IO ()
## String
"%10.2f"