{-# LANGUAGE NoImplicitPrelude,RebindableSyntax #-}
{-# LANGUAGE ScopedTypeVariables,FlexibleContexts #-}
module Algebra.Linear.Subspace
  (
    Subspace
  , fromGenerators
  , span
  , line
  , empty
  , fromRelations
  
  , inside
  , basis
  , union
  , intersection
  , pullback
  , image
  , kernel
  , dimension
  ) where


import           Algebra.Linear
  (
    Vector
  , Matrix
  
  , Relation(Relation)
  , getRelation
  
  , satisfies
  , solve
  , equations
  , matrixProduct
  )

import qualified Algebra.Field
import           NumericPrelude hiding (span)

import           Data.List       (transpose,genericLength)
import           Data.Proxy      (Proxy(Proxy))
import           Data.Reflection (Reifies,reflect)


data Subspace k
  = Generators Integer [Vector k]  -- These are linearly independent.
  | Relations Integer [Relation k] -- These might be redundant.

degree :: Subspace k -> Integer
degree (Generators d _) = d
degree (Relations  d _) = d

instance (Show k) => Show (Subspace k) where
  show (Generators d gs) = "Subspace (of space of dimension " ++ show d ++ ") generated by " ++ show gs
  show (Relations d rs)  = "Subspace (of space of dimension " ++ show d ++ ") defined by relations " ++ show (map getRelation rs)

fromGenerators :: (Algebra.Field.C k,Eq k) => Integer -> [Vector k] -> Subspace k
fromGenerators d gs = Relations d (equations d gs)

span :: (Algebra.Field.C k,Eq k) => [Vector k] -> Subspace k
span gs = fromGenerators d gs where
  d = case gs of
    g : _ -> genericLength g
    []    -> error "Algebra.Linear.Subspace.fromGenerators: no generators"

line :: (Algebra.Field.C k,Eq k) => Vector k -> Subspace k
line = span . (: [])

empty :: Integer -> Subspace k
empty d = Generators d []

fromRelations :: Integer -> [Relation k] -> Subspace k
fromRelations = Relations


inside :: (Algebra.Field.C k,Eq k) => Subspace k -> Vector k -> Bool
inside s v = all (v `satisfies`) (toRelations s)

basis :: (Algebra.Field.C k,Eq k) => Subspace k -> [Vector k]
basis (Generators _ gs) = gs
basis (Relations d rs) = solve d rs

toRelations :: (Algebra.Field.C k,Eq k) => Subspace k -> [Relation k]
toRelations (Generators d gs) = equations d gs
toRelations (Relations _ rs) = rs

union :: (Algebra.Field.C k,Eq k) => Subspace k -> Subspace k -> Subspace k
union a b = fromGenerators (sameDegree a b) $ basis a ++ basis b

intersection :: (Algebra.Field.C k,Eq k) => Subspace k -> Subspace k -> Subspace k
intersection a b = Relations (sameDegree a b) $ toRelations a ++ toRelations b where

sameDegree :: Subspace k -> Subspace k -> Integer
sameDegree a b
  | degree a == degree b = degree a
  | otherwise            = error "Algebra.Linear.Subspace.sameDegree: subspaces of different degree"

pullback :: (Algebra.Field.C k,Eq k) => Matrix k -> Subspace k -> Subspace k
pullback [] = error "Algebra.Linear.Subspace.pullback: empty matrix"
pullback m  =
    Relations d
  . map Relation
  . transpose
  . matrixProduct (transpose m)
  . transpose
  . map getRelation
  . toRelations
  . intersection (image m)
 where
  d = genericLength (head m)

image :: (Algebra.Field.C k,Eq k) => Matrix k -> Subspace k
image rows = fromGenerators d . transpose $ rows where
  d = genericLength rows

kernel :: Matrix k -> Subspace k
kernel []   = error "Algebra.Linear.Subspace.kernel: empty matrix"
kernel rows = Relations d . map Relation $ rows where
  d = genericLength (head rows)

dimension :: (Algebra.Field.C k,Eq k) => Subspace k -> Integer
dimension = genericLength . basis