vector norms moved to Vectors

4 files changed, 67 insertions, 126 deletions

diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 0f3ff0e..28063a9 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -61,9 +61,10 @@ module Numeric.LinearAlgebra.Algorithms (
     nullspacePrec,
     nullVector,
     nullspaceSVD,
+-- * Norms
+    Normed(..), NormType(..),
 -- * Misc
     eps, peps, i,
-    Normed(..), NormType(..),
 -- * Util
     haussholder,
     unpackQR, unpackHess,
@@ -172,9 +173,6 @@ thinSVD = {-# SCC "thinSVD" #-} thinSVD'
 singularValues :: Field t => Matrix t -> Vector Double
 singularValues = {-# SCC "singularValues" #-} sv'
 
-instance (Field t, RealOf t ~ Double) => Norm2 Matrix t where
-    norm2 m = singularValues m @> 0
-
 -- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
 --
 -- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
@@ -546,7 +544,7 @@ epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]]
 geps delta = head [ k | (k,g) <- epslist, g<delta]
 
 expGolub m = iterate msq f !! j
-    where j = max 0 $ floor $ log2 $ normInf m
+    where j = max 0 $ floor $ log2 $ pnorm Infinity m
           log2 x = log x / log 2
           a = m */ fromIntegral ((2::Int)^j)
           q = geps eps -- 7 steps
@@ -569,7 +567,7 @@ expGolub m = iterate msq f !! j
 {- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,
      based on a scaled Pade approximation.
 -}
-expm :: (Norm Vector t, Field t) => Matrix t -> Matrix t
+expm :: (Normed (Matrix t), Field t) => Matrix t -> Matrix t
 expm = expGolub
 
 --------------------------------------------------------------
@@ -585,11 +583,11 @@ It only works with invertible matrices that have a real solution. For diagonaliz
  [ 2.0, 2.25
  , 0.0,  2.0 ]@
 -}
-sqrtm ::  (Norm Vector t, Field t) => Matrix t -> Matrix t
+sqrtm ::  (Normed (Matrix t), Field t) => Matrix t -> Matrix t
 sqrtm = sqrtmInv
 
 sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x))
-    where fixedPoint (a:b:rest) | norm1 (fst a |-| fst b) < peps   = a
+    where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps   = a
                                 | otherwise = fixedPoint (b:rest)
           fixedPoint _ = error "fixedpoint with impossible inputs"
           f (y,z) = (0.5 .* (y |+| inv z),
@@ -632,33 +630,56 @@ luFact (l_u,perm) | r <= c    = (l ,u ,p, s)
 
 ---------------------------------------------------------------------------
 
-data NormType = Infinity | PNorm1 | PNorm2
+data NormType = Infinity | PNorm1 | PNorm2 | Frobenius
 
--- | Old class
 class Normed t where
     pnorm :: NormType -> t -> Double
 
-instance Norm Vector t => Normed (Vector t) where
-    pnorm PNorm1   = realToFrac . norm1
-    pnorm PNorm2   = realToFrac . normFrob
-    pnorm Infinity = realToFrac . normInf
+instance Normed (Vector Double) where
+    pnorm PNorm1    = norm1
+    pnorm PNorm2    = norm2
+    pnorm Infinity  = normInf
+    pnorm Frobenius = normInf
+
+instance Normed (Vector (Complex Double)) where
+    pnorm PNorm1    = realPart . norm1
+    pnorm PNorm2    = realPart . norm2
+    pnorm Infinity  = realPart . normInf
+    pnorm Frobenius = realPart . normInf
+
+instance Normed (Vector Float) where
+    pnorm PNorm1    = realToFrac . norm1
+    pnorm PNorm2    = realToFrac . norm2
+    pnorm Infinity  = realToFrac . normInf
+    pnorm Frobenius = realToFrac . normInf
+
+instance Normed (Vector (Complex Float)) where
+    pnorm PNorm1    = realToFrac . realPart . norm1
+    pnorm PNorm2    = realToFrac . realPart . norm2
+    pnorm Infinity  = realToFrac . realPart . normInf
+    pnorm Frobenius = realToFrac . realPart . normInf
+
 
 instance Normed (Matrix Double) where
-    pnorm PNorm1   = norm1
-    pnorm PNorm2   = norm2
-    pnorm Infinity = normInf
+    pnorm PNorm1    = maximum . map norm1 . toColumns
+    pnorm PNorm2    = (@>0) . singularValues
+    pnorm Infinity  = pnorm PNorm1 . trans
+    pnorm Frobenius = norm2 . flatten
 
 instance Normed (Matrix (Complex Double)) where
-    pnorm PNorm1   = norm1
-    pnorm PNorm2   = norm2
-    pnorm Infinity = normInf
+    pnorm PNorm1    = maximum . map (realPart.norm1) . toColumns
+    pnorm PNorm2    = (@>0) . singularValues
+    pnorm Infinity  = pnorm PNorm1 . trans
+    pnorm Frobenius = realPart . norm2 . flatten
 
 instance Normed (Matrix Float) where
-    pnorm PNorm1   = realToFrac . norm1
-    pnorm PNorm2   = norm2 . double -- not yet optimized for Float
-    pnorm Infinity = realToFrac . normInf
+    pnorm PNorm1    = realToFrac . maximum . map norm1 . toColumns
+    pnorm PNorm2    = realToFrac . (@>0) . singularValues . double
+    pnorm Infinity  = realToFrac . pnorm PNorm1 . trans
+    pnorm Frobenius = realToFrac . norm2 . flatten
 
 instance Normed (Matrix (Complex Float)) where
-    pnorm PNorm1   = realToFrac . norm1
-    pnorm PNorm2   = norm2 . double -- not yet optimized for Float
-    pnorm Infinity = realToFrac . normInf
+    pnorm PNorm1    = realToFrac . maximum . map (realPart.norm1) . toColumns
+    pnorm PNorm2    = realToFrac . (@>0) . singularValues . double
+    pnorm Infinity  = realToFrac . pnorm PNorm1 . trans
+    pnorm Frobenius = realToFrac . realPart . norm2 . flatten
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs
index ec08694..13d175b 100644
--- a/lib/Numeric/LinearAlgebra/Interface.hs
+++ b/lib/Numeric/LinearAlgebra/Interface.hs
@@ -19,7 +19,7 @@ In the context of the standard numeric operators, one-component vectors and matr
 -----------------------------------------------------------------------------
 
 module Numeric.LinearAlgebra.Interface(
-    (<>),(<.>),mulG, Adapt, adaptElements,
+    (<>),(<.>),
     (<\>),
     (.*),(*/),
     (<|>),(<->),
@@ -29,28 +29,20 @@ import Data.Packed.Vector
 import Data.Packed.Matrix
 import Numeric.LinearAlgebra.Algorithms
 import Numeric.LinearAlgebra.Linear
-import Data.Complex
-import Control.Arrow((***))
-
---import Numeric.GSL.Vector
 
 class Mul a b c | a b -> c where
  infixl 7 <>
  -- | Matrix-matrix, matrix-vector, and vector-matrix products.
  (<>)  :: Product t => a t -> b t -> c t
- mulG :: (Element r, Element s, Adapt r s t t, Product t) => a r -> b s -> c t
 
 instance Mul Matrix Matrix Matrix where
     (<>) = mXm
-    mulG a b = uncurry mXm (curry adapt a b)
 
 instance Mul Matrix Vector Vector where
     (<>) m v = flatten $ m <> (asColumn v)
-    mulG m v = flatten $ m `mulG` (asColumn v)
 
 instance Mul Vector Matrix Vector where
     (<>) v m = flatten $ (asRow v) <> m
-    mulG v m = flatten $ (asRow v) `mulG` m
 
 ---------------------------------------------------
 
@@ -124,48 +116,3 @@ a <|> b = joinH a b
 -- -- | Vertical concatenation of matrices and vectors.
 -- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
 a <-> b = joinV a b
-
-----------------------------------------------------
-
-class Adapt a b c d | a b -> c, a b -> d where
-    adapt :: Container k => (k a, k b) -> (k c, k d) 
-
---instance Adapt a a a a where
---    adapt = id *** id
-
-instance Adapt Float Float Float Float where
-    adapt = id *** id
-
-instance Adapt Double Double Double Double where
-    adapt = id *** id
-
-instance Adapt (Complex Float) (Complex Float) (Complex Float) (Complex Float) where
-    adapt = id *** id
-
-instance Adapt Float Double Double Double where
-    adapt = double *** id
-
-instance Adapt Double Float Double Double where
-    adapt = id *** double
-
-instance Adapt Float (Complex Float) (Complex Float) (Complex Float) where
-    adapt = complex *** id
-
-instance Adapt (Complex Float) Float (Complex Float) (Complex Float) where
-    adapt = id *** complex
-
-instance (Convert a, Convert (DoubleOf a), ComplexOf (DoubleOf a) ~ Complex Double) => Adapt a (Complex Double) (Complex Double) (Complex Double) where
-    adapt = complex.double *** id
-
-instance (Convert a, Convert (DoubleOf a), ComplexOf (DoubleOf a) ~ Complex Double) => Adapt (Complex Double) a (Complex Double) (Complex Double) where
-    adapt = id *** complex.double
-
-instance Adapt Double (Complex Float) (Complex Double) (Complex Double) where
-    adapt = complex *** double
-
-instance Adapt (Complex Float) Double (Complex Double) (Complex Double) where
-    adapt = double *** complex
-
-adaptElements:: (Adapt a b c d, Container k) => (k a, k b) -> (k c, k d)
-adaptElements p = adapt p
-
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
index 6c21a16..4d7f608 100644
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -24,8 +24,6 @@ module Numeric.LinearAlgebra.Linear (
     Product(..),
     mXm,mXv,vXm,
     outer, kronecker,
-    -- * Norms
-    Norm(..), Norm2(..),
     -- * Creation of numeric vectors
     constant, linspace
 ) where
@@ -42,53 +40,64 @@ import Control.Monad(ap)
 class Num e => Vectors a e where
     -- the C functions sumX are twice as fast as using foldVector
     vectorSum :: a e -> e
-    euclidean :: a e -> e
     absSum    :: a e -> e
     vectorMin :: a e -> e
     vectorMax :: a e -> e
     minIdx    :: a e -> Int
     maxIdx    :: a e -> Int
     dot       :: a e -> a e -> e
+    norm1     :: a e -> e
+    norm2     :: a e -> e
+    normInf   :: a e -> e
+
 
 instance Vectors Vector Float where
     vectorSum = sumF
-    euclidean = toScalarF Norm2
+    norm2     = toScalarF Norm2
     absSum    = toScalarF AbsSum
     vectorMin = toScalarF Min
     vectorMax = toScalarF Max
     minIdx    = round . toScalarF MinIdx
     maxIdx    = round . toScalarF MaxIdx
     dot       = dotF
+    norm1     = toScalarF AbsSum
+    normInf   = vectorMax . vectorMapF Abs
 
 instance Vectors Vector Double where
     vectorSum = sumR
-    euclidean = toScalarR Norm2
+    norm2     = toScalarR Norm2
     absSum    = toScalarR AbsSum
     vectorMin = toScalarR Min
     vectorMax = toScalarR Max
     minIdx    = round . toScalarR MinIdx
     maxIdx    = round . toScalarR MaxIdx
     dot       = dotR
+    norm1     = toScalarR AbsSum
+    normInf   = vectorMax . vectorMapR Abs
 
 instance Vectors Vector (Complex Float) where
     vectorSum = sumQ
-    euclidean = (:+ 0) . toScalarQ Norm2
+    norm2     = (:+ 0) . toScalarQ Norm2
     absSum    = (:+ 0) . toScalarQ AbsSum
     vectorMin = ap (@>) minIdx
     vectorMax = ap (@>) maxIdx
     minIdx    = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
     maxIdx    = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
     dot       = dotQ
+    norm1     = (:+ 0) . vectorSum . fst . fromComplex . vectorMapQ Abs
+    normInf   = (:+ 0) . vectorMax . fst . fromComplex . vectorMapQ Abs
 
 instance Vectors Vector (Complex Double) where
-    vectorSum  = sumC
-    euclidean = (:+ 0) . toScalarC Norm2
+    vectorSum = sumC
+    norm2     = (:+ 0) . toScalarC Norm2
     absSum    = (:+ 0) . toScalarC AbsSum
     vectorMin = ap (@>) minIdx
     vectorMax = ap (@>) maxIdx
     minIdx    = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
     maxIdx    = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
     dot       = dotC
+    norm1     = (:+ 0) . vectorSum . fst . fromComplex . vectorMapC Abs
+    normInf   = (:+ 0) . vectorMax . fst . fromComplex . vectorMapC Abs
 
 ----------------------------------------------------
 
@@ -268,39 +277,3 @@ kronecker a b = fromBlocks
               $ flatten a `outer` flatten b
 
 --------------------------------------------------
-
--- | simple norms
-class (Element t, RealFloat (RealOf t)) => Norm c t where
-    norm1   :: c t -> RealOf t
-    normInf  :: c t -> RealOf t
-    normFrob :: c t -> RealOf t
-
-instance Norm Vector Double where
-    normFrob = toScalarR Norm2
-    norm1 = toScalarR AbsSum
-    normInf = vectorMax . vectorMapR Abs
-
-instance Norm Vector Float where
-    normFrob = toScalarF Norm2
-    norm1 = toScalarF AbsSum
-    normInf = vectorMax . vectorMapF Abs
-
-instance (Norm Vector t, Vectors Vector t, RealElement t
-         , RealOf t ~ t, RealOf (Complex t) ~ t
-         ) => Norm Vector (Complex t) where
-    normFrob = normFrob . asReal
-    norm1 = norm1 . mapVector magnitude
-    normInf = vectorMax . mapVector magnitude
-
-instance Norm Vector t => Norm Matrix t where
-    normFrob = normFrob . flatten
-    norm1 = maximum . map norm1 . toColumns
-    normInf = norm1 . trans
-
-class Norm2 c t where
-    norm2 :: c t -> RealOf t
-
-instance Norm Vector t => Norm2 Vector t where
-    norm2 = normFrob
-
--- (the instance Norm2 Matrix t requires singular values and is defined later)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index d6bb338..d312e52 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -54,9 +54,9 @@ complex x = complex'' x
 debug x = trace (show x) x
 
 -- relative error
---dist :: (Normed t, Num t) => t -> t -> Double
+dist :: (Normed t, Num t) => t -> t -> Double
 dist a b = r
-    where norm = normInf
+    where norm = pnorm Infinity
           na = norm a
           nb = norm b
           nab = norm (a-b)