7 files changed, 143 insertions, 48 deletions

diff --git a/packages/base/src/Internal/CG.hs b/packages/base/src/Internal/CG.hs
index fd14212..758d130 100644
--- a/packages/base/src/Internal/CG.hs
+++ b/packages/base/src/Internal/CG.hs
@@ -45,13 +45,13 @@ cg sym at a (CGState p r r2 x _) = CGState p' r' r'2 x' rdx
     ap1 = a p
     ap  | sym       = ap1
         | otherwise = at ap1
-    pap | sym       = p <·> ap1
+    pap | sym       = p <.> ap1
         | otherwise = norm2 ap1 ** 2
     alpha = r2 / pap
     dx = scale alpha p
     x' = x + dx
     r' = r - scale alpha ap
-    r'2 = r' <·> r'
+    r'2 = r' <.> r'
     beta = r'2 / r2
     p' = r' + scale beta p
 
@@ -75,9 +75,9 @@ solveG mat ma meth rawb x0' ϵb ϵx
     b = mat rawb
     x0  = if x0' == 0 then konst 0 (dim b) else x0'
     r0  = b - a x0
-    r20 = r0 <·> r0
+    r20 = r0 <.> r0
     p0  = r0
-    nb2 = b <·> b
+    nb2 = b <.> b
     ok CGState {..}
         =  cgr2 <nb2*ϵb**2
         || cgdx < ϵx
diff --git a/packages/base/src/Internal/Container.hs b/packages/base/src/Internal/Container.hs
index 7fe5758..1c158ff 100644
--- a/packages/base/src/Internal/Container.hs
+++ b/packages/base/src/Internal/Container.hs
@@ -53,28 +53,23 @@ linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n
 
 --------------------------------------------------------------------------------
 
-infixl 7 <.>
--- | An infix synonym for 'dot'
-(<.>) :: Numeric t => Vector t -> Vector t -> t
-(<.>) = dot
+infixr 8 <.>
+{- | An infix synonym for 'dot'
 
+>>> vector [1,2,3,4] <.> vector [-2,0,1,1]
+5.0
 
-infixr 8 <·>, #>
+>>> let 𝑖 = 0:+1 :: C
+>>> fromList [1+𝑖,1] <.> fromList [1,1+𝑖]
+2.0 :+ 0.0
 
-{- | infix synonym for 'dot'
+-}
 
->>> vector [1,2,3,4] <·> vector [-2,0,1,1]
-5.0
+(<.>) :: Numeric t => Vector t -> Vector t -> t
+(<.>) = dot
 
->>> let 𝑖 = 0:+1 :: ℂ
->>> fromList [1+𝑖,1] <·> fromList [1,1+𝑖]
-2.0 :+ 0.0
 
-(the dot symbol "·" is obtained by Alt-Gr .)
 
--}
-(<·>) :: Numeric t => Vector t -> Vector t -> t
-(<·>) = dot
 
 
 {- | infix synonym for 'app'
@@ -91,6 +86,7 @@ infixr 8 <·>, #>
 fromList [140.0,320.0]
 
 -}
+infixr 8 #>
 (#>) :: Numeric t => Matrix t -> Vector t -> Vector t
 (#>) = mXv
 
@@ -264,6 +260,24 @@ instance Numeric Z
 sortVector :: (Ord t, Element t) => Vector t -> Vector t
 sortVector = sortV
 
+{- |
+
+>>> m <- randn 4 10
+>>> disp 2 m
+4x10
+-0.31   0.41   0.43  -0.19  -0.17  -0.23  -0.17  -1.04  -0.07  -1.24
+ 0.26   0.19   0.14   0.83  -1.54  -0.09   0.37  -0.63   0.71  -0.50
+-0.11  -0.10  -1.29  -1.40  -1.04  -0.89  -0.68   0.35  -1.46   1.86
+ 1.04  -0.29   0.19  -0.75  -2.20  -0.01   1.06   0.11  -2.09  -1.58
+
+>>> disp 2 $ m ?? (All, Pos $ sortIndex (m!1))
+4x10
+-0.17  -1.04  -1.24  -0.23   0.43   0.41  -0.31  -0.17  -0.07  -0.19
+-1.54  -0.63  -0.50  -0.09   0.14   0.19   0.26   0.37   0.71   0.83
+-1.04   0.35   1.86  -0.89  -1.29  -0.10  -0.11  -0.68  -1.46  -1.40
+-2.20   0.11  -1.58  -0.01   0.19  -0.29   1.04   1.06  -2.09  -0.75
+
+-}
 sortIndex :: (Ord t, Element t) => Vector t -> Vector I
 sortIndex = sortI
 
@@ -273,11 +287,50 @@ ccompare = ccompare'
 cselect :: (Container c t) => c I -> c t -> c t -> c t -> c t
 cselect = cselect'
 
+{- | Extract elements from positions given in matrices of rows and columns.
+
+>>> r
+(3><3)
+ [ 1, 1, 1
+ , 1, 2, 2
+ , 1, 2, 3 ]
+>>> c
+(3><3)
+ [ 0, 1, 5
+ , 2, 2, 1
+ , 4, 4, 1 ]
+>>> m
+(4><6)
+ [  0,  1,  2,  3,  4,  5
+ ,  6,  7,  8,  9, 10, 11
+ , 12, 13, 14, 15, 16, 17
+ , 18, 19, 20, 21, 22, 23 ]
+
+>>> remap r c m
+(3><3)
+ [  6,  7, 11
+ ,  8, 14, 13
+ , 10, 16, 19 ]
+
+The indexes are autoconformable.
+
+>>> c'
+(3><1)
+ [ 1
+ , 2
+ , 4 ]
+>>> remap r c' m
+(3><3)
+ [  7,  7,  7
+ ,  8, 14, 14
+ , 10, 16, 22 ]
+
+-}
 remap :: Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t
 remap i j m
     | minElement i >= 0 && maxElement i < fromIntegral (rows m) &&
       minElement j >= 0 && maxElement j < fromIntegral (cols m) = remapM i' j' m
-    | otherwise = error $ "out of range index in rmap"
+    | otherwise = error $ "out of range index in remap"
   where
     [i',j'] = conformMs [i,j]
     
diff --git a/packages/base/src/Internal/Element.hs b/packages/base/src/Internal/Element.hs
index 4007491..51d5686 100644
--- a/packages/base/src/Internal/Element.hs
+++ b/packages/base/src/Internal/Element.hs
@@ -80,7 +80,7 @@ breakAt c l = (a++[c],tail b) where
     (a,b) = break (==c) l
 
 --------------------------------------------------------------------------------
-
+-- | Specification of indexes for the operator '??'.
 data Extractor
     = All
     | Range Int Int Int
@@ -102,7 +102,32 @@ ppext (Drop n) = printf "Drop %d" n
 ppext (TakeLast n) = printf "TakeLast %d" n
 ppext (DropLast n) = printf "DropLast %d" n
 
+{- | General matrix slicing.
+
+>>> m
+(4><5)
+ [  0,  1,  2,  3,  4
+ ,  5,  6,  7,  8,  9
+ , 10, 11, 12, 13, 14
+ , 15, 16, 17, 18, 19 ]
+
+>>> m ?? (Take 3, DropLast 2)
+(3><3)
+ [  0,  1,  2
+ ,  5,  6,  7
+ , 10, 11, 12 ]
+
+>>> m ?? (Pos (idxs[2,1]), All)
+(2><5)
+ [ 10, 11, 12, 13, 14
+ ,  5,  6,  7,  8,  9 ]
+
+>>> m ?? (PosCyc (idxs[-7,80]), Range 4 (-2) 0)
+(2><3)
+ [ 9, 7, 5
+ , 4, 2, 0 ]
 
+-}
 infixl 9 ??
 (??)  :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t
 
@@ -328,27 +353,30 @@ r >< c = f where
 
 ----------------------------------------------------------------
 
--- | Creates a matrix with the first n rows of another matrix
 takeRows :: Element t => Int -> Matrix t -> Matrix t
 takeRows n mt = subMatrix (0,0) (n, cols mt) mt
+
 -- | Creates a matrix with the last n rows of another matrix
 takeLastRows :: Element t => Int -> Matrix t -> Matrix t
 takeLastRows n mt = subMatrix (rows mt - n, 0) (n, cols mt) mt
--- | Creates a copy of a matrix without the first n rows
+
 dropRows :: Element t => Int -> Matrix t -> Matrix t
 dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt
+
 -- | Creates a copy of a matrix without the last n rows
 dropLastRows :: Element t => Int -> Matrix t -> Matrix t
 dropLastRows n mt = subMatrix (0,0) (rows mt - n, cols mt) mt
--- |Creates a matrix with the first n columns of another matrix
+
 takeColumns :: Element t => Int -> Matrix t -> Matrix t
 takeColumns n mt = subMatrix (0,0) (rows mt, n) mt
+
 -- |Creates a matrix with the last n columns of another matrix
 takeLastColumns :: Element t => Int -> Matrix t -> Matrix t
 takeLastColumns n mt = subMatrix (0, cols mt - n) (rows mt, n) mt
--- | Creates a copy of a matrix without the first n columns
+
 dropColumns :: Element t => Int -> Matrix t -> Matrix t
 dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt
+
 -- | Creates a copy of a matrix without the last n columns
 dropLastColumns :: Element t => Int -> Matrix t -> Matrix t
 dropLastColumns n mt = subMatrix (0,0) (rows mt, cols mt - n) mt
diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs
index e4b1226..75e92a5 100644
--- a/packages/base/src/Internal/Matrix.hs
+++ b/packages/base/src/Internal/Matrix.hs
@@ -378,7 +378,6 @@ foreign import ccall unsafe "transL" ctransL :: TMM Z
 
 ----------------------------------------------------------------------
 
--- | Extracts a submatrix from a matrix.
 subMatrix :: Element a
           => (Int,Int) -- ^ (r0,c0) starting position
           -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs
index 3b27310..0d906bb 100644
--- a/packages/base/src/Internal/Modular.hs
+++ b/packages/base/src/Internal/Modular.hs
@@ -9,8 +9,8 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE GADTs #-}
-{-# LANGUAGE TypeFamilies #-}
-
+{-# LANGUAGE TypeFamilies  #-}
+{-# LANGUAGE TypeOperators #-}
 
 {- |
 Module      :  Internal.Modular
@@ -23,7 +23,7 @@ Proof of concept of statically checked modular arithmetic.
 -}
 
 module Internal.Modular(
-    Mod
+    Mod, type (./.)
 ) where
 
 import Internal.Vector
@@ -49,6 +49,9 @@ import Data.Complex
 newtype Mod (n :: Nat) t = Mod {unMod:: t}
   deriving (Storable)
 
+infixr 5 ./.
+type (./.) x n = Mod n x
+
 instance (Integral t, Enum t, KnownNat m) => Enum (Mod m t)
   where
     toEnum = l0 (\m x -> fromIntegral $ x `mod` (fromIntegral m))
diff --git a/packages/base/src/Internal/Numeric.hs b/packages/base/src/Internal/Numeric.hs
index ca17c23..efcde2c 100644
--- a/packages/base/src/Internal/Numeric.hs
+++ b/packages/base/src/Internal/Numeric.hs
@@ -481,14 +481,16 @@ step = step'
 -- ,   0.0, 100.0,   7.0,  8.0
 -- ,   0.0,   0.0, 100.0, 12.0 ]
 --
+-- >>> let chop x = cond (abs x) 1E-6 0 0 x
+--
 cond
-    :: (Ord e, Container c e)
+    :: (Ord e, Container c e, Container c x)
     => c e -- ^ a
     -> c e -- ^ b
-    -> c e -- ^ l
-    -> c e -- ^ e
-    -> c e -- ^ g
-    -> c e -- ^ result
+    -> c x -- ^ l
+    -> c x -- ^ e
+    -> c x -- ^ g
+    -> c x -- ^ result
 cond a b l e g = cselect' (ccompare' a b) l e g
 
 
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index 09ba21c..f08f710 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -85,7 +85,7 @@ type ℤ = Int
 type ℂ = Complex Double
 
 -- | imaginary unit
-iC :: ℂ
+iC :: C
 iC = 0:+1
 
 {- | Create a real vector.
@@ -94,7 +94,7 @@ iC = 0:+1
 fromList [1.0,2.0,3.0,4.0,5.0]
 
 -}
-vector :: [ℝ] -> Vector ℝ
+vector :: [R] -> Vector R
 vector = fromList
 
 {- | Create a real matrix.
@@ -108,8 +108,8 @@ vector = fromList
 -}
 matrix
   :: Int -- ^ number of columns
-  -> [ℝ] -- ^ elements in row order
-  -> Matrix ℝ
+  -> [R] -- ^ elements in row order
+  -> Matrix R
 matrix c = reshape c . fromList
 
 
@@ -260,34 +260,34 @@ norm = pnorm PNorm2
 
 class Normed a
   where
-    norm_0   :: a -> ℝ
-    norm_1   :: a -> ℝ
-    norm_2   :: a -> ℝ
-    norm_Inf :: a -> ℝ
+    norm_0   :: a -> R
+    norm_1   :: a -> R
+    norm_2   :: a -> R
+    norm_Inf :: a -> R
 
 
-instance Normed (Vector ℝ)
+instance Normed (Vector R)
   where
     norm_0 v = sumElements (step (abs v - scalar (eps*normInf v)))
     norm_1 = pnorm PNorm1
     norm_2 = pnorm PNorm2
     norm_Inf = pnorm Infinity
 
-instance Normed (Vector ℂ)
+instance Normed (Vector C)
   where
     norm_0 v = sumElements (step (fst (fromComplex (abs v)) - scalar (eps*normInf v)))
     norm_1 = pnorm PNorm1
     norm_2 = pnorm PNorm2
     norm_Inf = pnorm Infinity
 
-instance Normed (Matrix ℝ)
+instance Normed (Matrix R)
   where
     norm_0 = norm_0 . flatten
     norm_1 = pnorm PNorm1
     norm_2 = pnorm PNorm2
     norm_Inf = pnorm Infinity
 
-instance Normed (Matrix ℂ)
+instance Normed (Matrix C)
   where
     norm_0 = norm_0 . flatten
     norm_1 = pnorm PNorm1
@@ -323,12 +323,22 @@ instance Normed (Vector (Complex Float))
     norm_Inf = norm_Inf . double
 
 
-norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> ℝ
+norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R
 norm_Frob = norm_2 . flatten
 
-norm_nuclear :: Field t => Matrix t -> ℝ
+norm_nuclear :: Field t => Matrix t -> R
 norm_nuclear = sumElements . singularValues
 
+{- | Check if the absolute value or complex magnitude is greater than a given threshold
+
+>>> magnit 1E-6 (1E-12 :: R)
+False
+>>> magnit 1E-6 (3+iC :: C)
+True
+>>> magnit 0 (3 :: I ./. 5)
+True
+
+-}
 magnit :: (Element t, Normed (Vector t)) => R -> t -> Bool
 magnit e x = norm_1 (fromList [x]) > e