packages feed

np-linear-0.1.1.1: src/Algebra/Module/Free.hs

{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE FlexibleInstances,MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE DeriveGeneric #-}
module Algebra.Module.Free
  (
    Free
  , first
  , second
    
  , Map
  , Domain
  , Codomain
  , apply
  
  , matrix
  , vector
  , fromVector
  
  , Enumerable
  , enumerate
  , coefficient
  , basisVector
  , fromList
  , FreeBasis(FreeBasis)
  ) where


import           Algebra.Linear (Matrix,Vector)

import qualified Algebra.Additive
import qualified Algebra.Module
import qualified Algebra.ModuleBasis
import qualified Algebra.Ring
import NumericPrelude

-- import           Data.Bifunctor (Bifunctor,first,second)
import           Data.Binary  (Binary)
import           Data.List    (intercalate)
import qualified Data.Map.Strict as Map
import           Data.Proxy   (Proxy(Proxy))
import           GHC.Generics (Generic)


newtype Free k t
  = Free (Map.Map t k)
  deriving (Generic)

-- The natural Bifunctor Free instance is not possible, because of the constraints.
first :: (k₁ -> k₂) -> Free k₁ t -> Free k₂ t
first f (Free m) = Free (fmap f m)
second :: (Algebra.Additive.C k,Ord t₂) => (t₁ -> t₂) -> Free k t₁ -> Free k t₂
second g (Free m) = Free (Map.mapKeysWith (+) g m)

instance (Algebra.Additive.C k,Ord t) => Algebra.Additive.C (Free k t) where
  zero = Free Map.empty
  negate (Free m) = Free (fmap negate m)
  Free m₁ + Free m₂ = Free $ Map.unionWith (+) m₁ m₂

instance (Algebra.Ring.C k,Ord t) => Algebra.Module.C k (Free k t) where
  c *> (Free m) = Free (fmap (c *) m)

instance (Show t,Ord t,Show k) => Show (Free k t) where
  show (Free m) = intercalate " + " . map (\ (x,c) -> show c ++ " · " ++ show x) $ Map.toList m

instance (Binary t,Binary k) => Binary (Free k t)

-- 'b' is a basis for the module 'm'
class ModuleBasis b m where
  type Scalar b m
  type BasisElement b m
  coefficient :: b -> BasisElement b m -> m -> Scalar b m
  basisVector :: b -> BasisElement b m -> m

fromList :: (ModuleBasis b m,Algebra.Module.C (Scalar b m) m) => b -> [(BasisElement b m,Scalar b m)] -> m
fromList b = sum . map (\ (x,c) -> c *> basisVector b x)

data FreeBasis
  = FreeBasis

instance (Algebra.Ring.C k,Ord t) => ModuleBasis FreeBasis (Free k t) where
  type Scalar FreeBasis (Free k t) = k
  type BasisElement FreeBasis (Free k t) = t
  coefficient FreeBasis x (Free m) = Map.findWithDefault zero x m
  basisVector FreeBasis x = Free (Map.singleton x one)

class (Ord t) => Enumerable p t where
  enumerate :: p -> [t]

-- instance forall k t. (Enumerable t,Algebra.Ring.C k) => Algebra.ModuleBasis.C k (Free k t) where
--   dimension _ _ = length $ enumerate (Proxy :: Proxy t)
--   flatten (Free m) = map (flip coefficient m) $ enumerate (Proxy :: Proxy t)
--   basis _ = map basisVector $ enumerate (Proxy :: Proxy t)


class Map k m where
  type Domain m
  type Codomain m
  apply :: m -> Domain m -> Codomain m

instance (Algebra.Ring.C k,Ord t₁,Ord t₂) => Map k (Free k (t₁,t₂)) where
  type Domain (Free k (t₁,t₂)) = Free k t₁
  type Codomain (Free k (t₁,t₂)) = Free k t₂
  apply (Free m) (Free v₁) = Free $ Map.foldrWithKey
    (\ (x₁,x₂) c -> case Map.lookup x₁ v₁ of
      Nothing -> id
      Just d  -> Map.insertWith (+) x₂ (c * d)
    )
    Map.empty
    m

matrix :: forall k t₁ t₂ p₁ p₂. (Enumerable p₁ t₁,Enumerable p₂ t₂,Algebra.Ring.C k) => p₁ -> p₂ -> Free k (t₁,t₂) -> Matrix k
matrix p₁ p₂ m = [[coefficient FreeBasis (x₁,x₂) m | x₁ <- enumerate p₁] | x₂ <- enumerate p₂]

vector :: forall k t p. (Enumerable p t,Algebra.Ring.C k) => p -> Free k t -> Vector k
vector p v = [coefficient FreeBasis x v | x <- enumerate p]

fromVector :: forall k t p. (Enumerable p t,Algebra.Ring.C k) => p -> Vector k -> Free k t
fromVector p = Free . Map.fromList . zip (enumerate p)

-- data X2 = S1 | S2 deriving (Eq,Ord)
-- instance Enumerable X2 where
--   enumerate _ = [S1,S2]
-- data X3 = T1 | T2 | T3 deriving (Eq,Ord)
-- instance Enumerable X3 where
--   enumerate _ = [T1,T2,T3]
-- 
-- m :: Free Integer (X2,X3)
-- m = (2 :: Integer) *> basisVector FreeBasis (S2,T2) - (3 :: Integer) *> basisVector FreeBasis (S1,T3) :: Free Integer (X2,X3)