hTensor 0.1.1 → 0.5.0

raw patch · 14 files changed

+1023/−422 lines, 14 filesdep +randomdep ~hmatrixPVP ok

version bump matches the API change (PVP)

Dependencies added: random

Dependency ranges changed: hmatrix

API changes (from Hackage documentation)

Files

• examples/array.hs +0/−59
• examples/exterior.hs +0/−41
• examples/geom.hs +0/−29
• hTensor.cabal +8/−7
• lib/Numeric/LinearAlgebra/Array.hs +4/−3
• lib/Numeric/LinearAlgebra/Array/Decomposition.hs +192/−0
• lib/Numeric/LinearAlgebra/Array/Display.hs +171/−0
• lib/Numeric/LinearAlgebra/Array/Internal.hs +183/−254
• lib/Numeric/LinearAlgebra/Array/Simple.hs +11/−4
• lib/Numeric/LinearAlgebra/Array/Solve.hs +314/−0
• lib/Numeric/LinearAlgebra/Array/Util.hs +111/−6
• lib/Numeric/LinearAlgebra/Exterior.hs +8/−8
• lib/Numeric/LinearAlgebra/Multivector.hs +2/−1
• lib/Numeric/LinearAlgebra/Tensor.hs +19/−10

@@ -1,59 +0,0 @@-import Numeric.LinearAlgebra.Array-import Numeric.LinearAlgebra.Array.Util-import Control.Applicative---- 'listArray' specialized for Array Double-infixl 9 #-(#) :: [Int] -> [Double] -> Array Double-(#) = listArray--(<|) :: Name -> [Array Double] -> Array Double-infixl 8 <|-n <| ls = index n ls--i = ("i" <|)-j = ("j" <|)-k = ("k" <|)--sh x = putStrLn . formatFixed 2 $ x--a = [3,4,2] # [1..]-b = [2,3] # [5,6]-c = 7 :: Array Double-s = [2,2,2,2] # [1..]--t = [3,3,3,3]#[1 ..]!"ijkl"--q = [2,4,3] # (fun <$> r 2 <*> r 4 <*> r 3) !"ijk"-    where r k = [1..k]-          fun = \i j k -> i*2*j-k--mk ds f = ds # map f (sequence $ map (enumFromTo 1 . fromIntegral) $ ds)--m = j [i[2,0,0], i[1,0,1], i[0,3,0]] ~> "ij"--main = do-    putStrLn "8-dimensional array"-    sh $ (replicate 8 2)#[1::Double ..]!(take 8 ['a'..])-    -------------------------    putStrLn "different display formats"-    sh $ a!"ijk"-    printA "%7.3f" a-    putStrLn . formatScaled 2 $ a!"ijk"-    sh . noIdx $ a-    -------------------------    putStrLn "array defined using a function"-    sh q-    sh $ mk [2,4,3] (\[i,j,k] -> i*2*j-k) !"ijk"-    -------------------------    putStrLn "contraction"-    sh t-    sh $ t!"ijkk"-    -------------------------    putStrLn "tensor product"-    sh $ m-    sh $ (t !"pqrs" * m!"kr") ~> "pqks"-    -------------------------    putStrLn "automatic conformability"-    sh $ j[1,2,3] + k[10,20]-    sh $ k [m, 3*m-1, 7, i[1,2,3]]

@@ -1,41 +0,0 @@-import Numeric.LinearAlgebra.Exterior-import Numeric.LinearAlgebra.Array.Util(formatFixed,asMatrix)-import Numeric.LinearAlgebra (det,(><))--printAS = print . asMultivector--sh = putStrLn . formatFixed 2---- 'listTensor' specialized for Tensor Double-infixl 9 #-(#) :: [Int] -> [Double] -> Tensor Double-(#) = listTensor--m = [3,-3]#[ 1,2,5,-             1,2,8,-            -2,0,4]--eps = leviCivita 3--a = vector [1,0,0] /\ vector [0,1,0]-b = vector [2,0,100,0] /\ vector [0,3,0,0] /\ vector [0,0,4,0]--im = eps!"ijb"* m!"pi" * m!"qj" * cov eps!"apq"--main = do-    putStrLn "exterior product"-    print a-    sh a-    printAS a-    printAS b-    putStrLn "\ndeterminant"-    printAS $ cov eps!"pqr" * m!"pi" * m!"qj" * m!"rk"-    print $ det (asMatrix m)-    putStrLn "\ninverse"-    sh $ im-    sh $ im!"ik" * m!"kj"-    putStrLn "\nmeet and join"-    printAS $ (vector [1,0,1] /\ vector [0,1,0]) \/ (vector [1,1,0] /\ vector [0,0,1])-    putStrLn "\nEuclidean inner product of r-vectors"-    printAS $ (vector [1,0,1] /\ vector [0,1,0]) `inner` (vector[1,1,1])-    print $ vector[3,5] `inner` vector[2,1]

@@ -1,29 +0,0 @@-import Numeric.LinearAlgebra.Multivector--o  = e 4-x = vector[1,0,0,1]-y = o + e 2-z = vector[0,0,1] + e 4--p = x /\ y /\ z--p1 = o /\ x /\ y-p2 = z /\ vector[1,0,1,1] /\ vector[0,1,1,1]--l = o /\ vector[1,1,1,1]--rot = rotor 3 (pi/4) (e 3)--l' = rot * l * rever rot--inh v = v / (v -| e 4) - e 4--main = do-     print l-     print $ l \/ p1-     print $ l \/ p2-     print $ l \/ p-     print $ inh $ l \/ p1-     print $ inh $ l \/ p2-     print $ inh $ l' \/ p1-     print $ inh $ l' \/ p2

@@ -1,5 +1,5 @@ Name:               hTensor-Version:            0.1.1+Version:            0.5.0 License:            GPL License-file:       LICENSE Author:             Alberto Ruiz@@ -21,19 +21,17 @@  and operations typically work on whole structures which can be assembled  and decomposed using simple primitives. Arguments are automatically made conformant  by replicating them along extra dimensions appearing in an operation.- There is preliminary support for Geometric Algebra.+ There is preliminary support for Geometric Algebra and for solving multilinear systems.  .  A tutorial can be found in the website of the project.  Category:           Math-tested-with:        GHC ==6.10.3+tested-with:        GHC ==6.12.1  cabal-version:      >=1.2 build-type:         Simple -extra-source-files: examples/array.hs-                    examples/exterior.hs-                    examples/geom.hs+extra-source-files:  flag splitBase     description:    Choose the new smaller, split-up base package.@@ -44,7 +42,7 @@     else       build-depends:    base < 3 -    Build-Depends:      haskell98, hmatrix >= 0.5, containers+    Build-Depends:      haskell98, hmatrix >= 0.8.2, containers, random      hs-source-dirs:     lib     Exposed-modules:    Numeric.LinearAlgebra.Array.Util@@ -52,8 +50,11 @@                         Numeric.LinearAlgebra.Tensor                         Numeric.LinearAlgebra.Exterior                         Numeric.LinearAlgebra.Multivector+                        Numeric.LinearAlgebra.Array.Solve+                        Numeric.LinearAlgebra.Array.Decomposition      other-modules:      Numeric.LinearAlgebra.Array.Internal+                        Numeric.LinearAlgebra.Array.Display                         Numeric.LinearAlgebra.Array.Simple      ghc-prof-options:   -auto-all

@@ -21,14 +21,15 @@     listArray,     scalar,     index,-    (!),(~>),+    (!),(!>),(~>),     (.*),     printA ) where  import Numeric.LinearAlgebra.Array.Simple import Numeric.LinearAlgebra.Array.Util-import Numeric.LinearAlgebra.Array.Internal(printA)+import Numeric.LinearAlgebra.Array.Internal(namesR)+import Numeric.LinearAlgebra.Array.Display(printA) import Data.Packed(Vector)  -- | Create an 'Array' from a list of parts (@index = 'newIndex' 'None'@).@@ -42,7 +43,7 @@ (.*) = zipArray (*)  instance (Coord t, Compat i) => Eq (NArray i t) where-    t1 == t2 = sameStructure t1 t2 && coords t1 == coords (reorder (names t1) t2)+    t1 == t2 = sameStructure t1 t2 && coords t1 == coords (reorder (namesR t1) t2)  instance (Show (NArray i t), Coord t, Compat i) => Num (NArray i t) where     (+) = zipArray (+)

@@ -0,0 +1,192 @@+{-# LANGUAGE FlexibleContexts #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Packed.Array.Decomposition+-- Copyright   :  (c) Alberto Ruiz 2009+-- License     :  GPL+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+--+-- Common multidimensional array decompositions. See the paper by Kolda & Balder.+--+-----------------------------------------------------------------------------++module Numeric.LinearAlgebra.Array.Decomposition (+    -- * HOSVD+    hosvd, hosvd', truncateFactors,+    -- * CP+    cpAuto, cpRun, cpInitRandom, cpInitSvd,+    -- * Utilities+    ALSParam(..), defaultParameters+) where++import Numeric.LinearAlgebra.Array+import Numeric.LinearAlgebra.Array.Internal(seqIdx,namesR,sizesR,renameRaw)+import Numeric.LinearAlgebra.Array.Util+import Numeric.LinearAlgebra.Array.Solve+import Numeric.LinearAlgebra hiding ((.*), scalar)+import Data.List+import System.Random+--import Control.Parallel.Strategies++{- | Full version of 'hosvd'.++    The first element in the result pair is a list with the core (head) and rotations so that+    t == product (fst (hsvd' t)).++    The second element is a list of rank and singular values along each mode,+    to give some idea about core structure.+-}+hosvd' :: Array Double -> ([Array Double],[(Int,Vector Double)])+hosvd' t = (factors,ss) where+    (rs,ss) = unzip $ map usOfSVD $ flats t+    n = length rs+    dummies = take n $ seqIdx (2*n) "" \\ (namesR t)+    axs = zipWith (\a b->[a,b]) dummies (namesR t)+    factors = renameRaw core dummies : zipWith renameRaw (map (fromMatrix None None . trans) rs) axs+    core = product $ renameRaw t dummies : zipWith renameRaw (map (fromMatrix None None) rs) axs++{- | Multilinear Singular Value Decomposition (or Tucker's method, see Lathauwer et al.).++    The result is a list with the core (head) and rotations so that+    t == product (hsvd t).++    The core and the rotations are truncated to the rank of each mode.++    Use 'hosvd'' to get full transformations and rank information about each mode.++-}+hosvd :: Array Double -> [Array Double]+hosvd a = truncateFactors rs h where+    (h,info) = hosvd' a+    rs = map fst info+++-- get the matrices of the flattened tensor for all dimensions+flats t = map (flip fibers t) (namesR t)+++--check trans/ctrans+usOfSVD m = if rows m < cols m+        then let (s2,u) = eigSH' $ m <> ctrans m+                 s = sqrt (abs s2)+              in (u,r s)+        else let (s2,v) = eigSH' $ ctrans m <> m+                 s = sqrt (abs s2)+                 u = m <> v <> pinv (diag s)+              in (u,r s)+    where r s = (ranksv (sqrt eps) (max (rows m) (cols m)) (toList s), s)+                -- (rank m, sv m) where sv m = s where (_,s,_) = svd m+++ttake ns t = (foldl1' (.) $ zipWith (onIndex.take) ns (namesR t)) t++-- | Truncate a 'hosvd' decomposition from the desired number of principal components in each dimension.+truncateFactors :: [Int] -> [Array Double] -> [Array Double]+truncateFactors _ [] = []+truncateFactors ns (c:rs) = ttake ns c : zipWith f rs ns+    where f r n = onIndex (take n) (head (namesR r)) r++------------------------------------------------------------------------++frobT = pnorm PNorm2 . coords++------------------------------------------------------------------------++unitRows [] = error "unitRows []"+unitRows (c:as) = foldl1' (.*) (c:xs) : as' where+    (xs,as') = unzip (map g as)+    g a = (x,a')+        where n = head (namesR a) -- hmmm+              rs = parts a n+              scs = map frobT rs+              x = diagT scs (order c) `renameRaw` (namesR c)+              a' = (zipWith (.*) (map (scalar.recip) scs)) `onIndex` n $ a+++{- | Basic CP optimization for a given rank. The result includes the obtained sequence of errors.++For example, a rank 3 approximation can be obtained as follows, where initialization+is based on the hosvd:++@+(y,errs) = cpRank 3 t+     where cpRank r t = cpRun (cpInitSvd (fst $ hosvd' t) r) defaultParameters t+@++-}+cpRun :: [Array Double] -- ^ starting point+      -> ALSParam None Double     -- ^ optimization parameters+      -> Array Double -- ^ input array+      -> ([Array Double], [Double]) -- ^ factors and error history+cpRun s0 params t = (unitRows $ head s0 : sol, errs) where+    (sol,errs) = mlSolve params [head s0] (tail s0) t++++{- | Experimental implementation of the CP decomposition, based on alternating+     least squares. We try approximations of increasing rank, until the relative reconstruction error is below a desired percent of Frobenius norm (epsilon).++     The approximation of rank k is abandoned if the error does not decrease at least delta% in an iteration.++    Practical usage can be based on something like this:++@+cp finit d e t = cpAuto (finit t) defaultParameters {delta = d, epsilon = e} t++cpS = cp (InitSvd . fst . hosvd')+cpR s = cp (cpInitRandom s)+@++     So we can write++@+ \-\- initialization based on hosvd+y = cpS 0.01 1E-6 t++ \-\- (pseudo)random initialization+z = cpR seed 0.1 0.1 t+@++-}+cpAuto :: (Int -> [Array Double]) -- ^ Initialization function for each rank+       -> ALSParam None Double    -- ^ optimization parameters+       -> Array Double -- ^ input array+       -> [Array Double] -- ^ factors+cpAuto finit params t = fst . head . filter ((<epsilon params). head . snd)+                      . map (\r->cpRun (finit r) params t) $ [1 ..]++----------------------++-- | cp initialization based on the hosvd+cpInitSvd :: [NArray None Double] -- ^ hosvd decomposition of the target array+          -> Int                  -- ^ rank+          -> [NArray None Double] -- ^ starting point+cpInitSvd (hos) k = d:as+    where c:rs = hos+          as = trunc (replicate (order c) k) rs+          d = diagT (replicate k 1) (order c) `renameO` (namesR c)+          trunc ns xs = zipWith f xs ns+              where f r n = onIndex (take n . cycle) (head (namesR r)) r++cpInitSeq rs t k = ones:as where+    n = order t+    auxIndx = take n $ seqIdx (2*n) "" \\ namesR t+              --take (order t) $ map return ['a'..] \\ namesR t+    ones = diagT (replicate k 1) (order t) `renameO` auxIndx+    ts = takes (map (*k) (sizesR t)) rs+    as = zipWith4 f ts auxIndx (namesR t) (sizesR t)+    f c n1 n2 p = (listArray [k,p] c) `renameO` [n1,n2]++takes [] _ = []+takes (n:ns) xs = take n xs : takes ns (drop n xs)++-- | pseudorandom cp initialization from a given seed+cpInitRandom :: Int        -- ^ seed+             -> NArray i t -- ^ target array to decompose+             -> Int        -- ^ rank+             -> [NArray None Double] -- ^ random starting point+cpInitRandom seed = cpInitSeq (randomRs (-1,1) (mkStdGen seed))++----------------------------------------------------------------------

@@ -0,0 +1,171 @@+-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Packed.Array.Display+-- Copyright   :  (c) Alberto Ruiz 2009+-- License     :  GPL+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Formatting utilities++-----------------------------------------------------------------------------++module Numeric.LinearAlgebra.Array.Display (+    formatArray, formatFixed, formatScaled, printA, dummyAt, noIdx, showBases,+) where++import Numeric.LinearAlgebra.Array.Internal+import Data.Packed+import Data.List+import Text.Printf++showBases x = f $ concatMap (shbld) x+    where   shbld (c,[]) = shsign c ++ showc c+            shbld (c,l) = shsign c ++ g (showc c) ++ "{"++ concatMap show l++"}"+            shsign c = if c < 0 then " - " else " + "+            showc c+                | abs (fromIntegral (round c :: Int) - c) <1E-10  = show (round $ abs c::Int)+                | otherwise = printf "%.3f" (abs c)+            f (' ':'+':' ':rs) = rs+            f (' ':'-':' ':rs) = '-':rs+            f a = a+            g "1" = ""+            g a = a++---------------------------------------------------------++data Rect = Rect { li :: Int, co :: Int, els :: [String] }++rect s = pad r c (Rect r 0 ss)+    where ss = lines s+          r  = length ss+          c  = maximum (map length ss)++pad nr nc (Rect r c ss) = Rect (r+r') (c+c') ss'' where+    r' = max 0 (nr-r)+    c' = max 0 (nc-c)+    ss' = map (padH nc) ss+    ss'' = replicate r' (replicate nc '-') ++ ss'+    padH l s = take (l-length s) (" | "++repeat ' ') ++ s++dispH :: Int -> [Rect] -> Rect+dispH k rs = Rect nr nc nss where+    nr = maximum (map li rs)+    nss' = mapTail (\x-> pad nr (co x + k) x) rs+    nss = foldl1' (zipWith (++)) (map els nss')+    nc = length (head nss)++dispV :: Int -> [Rect] -> Rect+dispV k rs = Rect nr nc nss where+    nc = maximum (map co rs)+    nss' = mapTail (\x-> pad (li x + k) nc x) rs+    nss = concatMap els nss'+    nr = length nss++mapTail f (a:b) = a : map f b+mapTail _ x     = x+++formatAux f x = unlines . addds . els . fmt ms $ x where+    fmt [] _ = undefined -- cannot happen+    fmt (g:gs) t+        | order t == 0 = rect (f (coords t @> 0))+        | order t == 1 =  rect $ unwords $ map f (toList $ coords t)+        | order t == 2 =  decor t $ rect $ w1 $ format " " f (reshape (iDim $ last $ dims t) (coords t))+        | otherwise    = decor t (g ps)+      where ps = map (fmt gs ) (partsRaw t (head (namesR t)))+    ds = showNice (filter ((/='*').head.iName) $ dims x)+    addds = if null ds then (showRawDims (dims x) :) else (ds:)+    w1 = unlines . map (' ':) . lines+    ms = cycle [dispV 1, dispH 2]+    decor t | odd (order t) = id+            | otherwise = decorLeft  (namesR t!!0) . decorUp (namesR t!!1)+++showNice x = unwords . intersperse "x" . map show $ x+showRawDims = showNice . map iDim . filter ((/="*").iName)++------------------------------------------------------++-- | Show a multidimensional array as a nested 2D table.+formatArray :: (Coord t, Compat i)+      => (t -> String) -- ^ format function (eg. printf \"5.2f\")+      -> NArray i t+      -> String+formatArray f t | odd (order t) = formatAux f (dummyAt 0 t)+            | otherwise    = formatAux f t+++decorUp s rec+    | head s == '*' = rec+    | otherwise     = dispV 0 [rs,rec]+  where+    c = co rec+    c1 = (c - length s) `div` 2+    c2 = c - length s - c1+    rs = rect $ replicate c1 ' ' ++ s ++ replicate c2 ' '++decorLeft s rec+    | head s == '*' = rec+    | otherwise     = dispH 0 [rs,rec]+  where+    c = li rec+    r1 = (c - length s+1) `div` 2+    r2 = c - length s - r1+    rs = rect $ unlines $ replicate r1 spc ++ s : replicate (r2) spc+    spc = replicate (length s) ' '++------------------------------------------------------++-- | Print the array as a nested table with the desired format (e.g. %7.2f) (see also 'formatArray', and 'formatScaled').+printA :: (Coord t, Compat i, PrintfArg t) => String -> NArray i t -> IO ()+printA f t = putStrLn (formatArray (printf f) t)+++-- | Show the array as a nested table with autoscaled entries.+formatScaled :: (Compat i)+      => Int -- ^ number of of decimal places+      -> NArray i Double+      -> String+formatScaled dec t = unlines (('(':d++")  E"++show o) : m)+    where ss = formatArray (printf fmt. g) t+          d:m = lines ss+          g x | o >= 0    = x/10^(o::Int)+              | otherwise = x*10^(-o)+          o = floor $ maximum $ map (logBase 10 . abs) $ toList $ coords t+          fmt = '%':show (dec+3) ++ '.':show dec ++"f"++-- | Show the array as a nested table with a \"\%.nf\" format. If all entries+-- are approximate integers the array is shown without the .00.. digits.+formatFixed :: (Compat i)+      => Int -- ^ number of of decimal places+      -> NArray i Double+      -> String+formatFixed dec t+    | isInt t   = formatArray (printf ('%': show (width t) ++".0f")) t+    | otherwise = formatArray (printf ('%': show (width t+dec+1) ++"."++show dec ++"f")) t++isInt = all lookslikeInt . toList . coords+lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx+    where shx = show x+-- needsSign t = vectorMin (coords t) < 0+-- width :: Compat i => NArray i Double -> Int+width = maximum . map (length . (printf "%.0f"::Double->String)) . toList . coords+-- width t = k + floor (logBase 10 (max 1 $ vectorMax (abs $ coords t))) :: Int+--      where k | needsSign t = 2+--              | otherwise = 1++------------------------------------------------------++-- | Insert a dummy index of dimension 1 at a given level (for formatting purposes).+dummyAt :: Int -> NArray i t -> NArray i t+dummyAt k t = mkNArray d' (coords t) where+    (d1,d2) = splitAt k (dims t)+    d' = d1 ++ d : d2+    d = Idx (iType (head (dims t))) 1 "*"++-- | Rename indices so that they are not shown in formatted output.+noIdx :: Compat i => NArray i t -> NArray i t+noIdx t = renameSuperRaw t (map ('*':) (namesR t))

@@ -24,63 +24,67 @@ module Numeric.LinearAlgebra.Array.Internal (     -- * Data structures     NArray, Idx(..), Name,-    rank, names, size, typeOf , dims, coords,+    order, namesR, names, size, sizesR, sizes, typeOf , dims, coords,     Compat(..),     -- * Array creation     scalar,     mkNArray,     fromVector, fromMatrix,     -- * Array manipulation-    rename,(!),-    parts,+    renameRaw,+    parts, partsRaw,     (|*|),+    analyzeProduct,+    smartProduct,     zipArray,     mapArray,     extract,     onIndex,     -- * Utilities-    reorder, (~>),+    seqIdx,+    reorder,     sameStructure,     conformable,     makeConformant,-    mapTypes,-    renameRaw,-    formatArray, formatFixed, formatScaled, printA,-    showBases,+    mapTypes, mapNames,+    renameSuperRaw, renameExplicit,     newIndex,-    dummyAt, noIdx,     basisOf,     common,+    selDims, mapDims,+    takeDiagT, atT,+    firstIdx, fibers, matrixator, matrixatorFree,     Coord,-    asMatrix, asVector, asScalar+    asMatrix, asVector, asScalar,+    debug ) where  import Data.Packed import Data.List-import Numeric.LinearAlgebra(outer,multiply,Field)+import Numeric.LinearAlgebra((<>),Field) import Control.Applicative import Data.Function(on)-import Text.Printf+-- import Control.Parallel.Strategies+import Debug.Trace +debug m f x = trace (m ++ show (f x)) x+ -- | Types that can be elements of the multidimensional arrays. class (Num (Vector t), Field t) => Coord t instance Coord Double instance Coord (Complex Double) --- import Debug.Trace--- --- debug s f x = trace (s ++ ": " ++ show (f x)) x- -- | indices are denoted by strings, (frequently single-letter) type Name = String  -- | Dimension descriptor.-data Idx i = Idx { iDim  :: Int+data Idx i = Idx { iType :: i+                 , iDim  :: Int                  , iName :: Name-                 , iType :: i                  } deriving (Eq) -+instance Eq i => Ord (Idx i) where+    compare = compare `on` iName  -- | A multidimensional array with index type i and elements t. data NArray i t = A { dims   :: [Idx i]   -- ^ Get detailed dimension information about the array.@@ -88,7 +92,6 @@                                           -- flattened structure (in the order specified by 'dims').                     } - -- | development function not intended for the end user mkNArray :: [Idx i] -> Vector a -> NArray i a mkNArray [] _ = error "array with empty dimensions, use scalar"@@ -104,56 +107,92 @@ scalar :: Coord t => t -> NArray i t scalar x = A [] (fromList [x]) ---- | 'rename' the indices with single-letter names. Equal indices of compatible type are contracted out.-infixl 8 !-(!) :: (Coord t, Compat i)-       => NArray i t-       -> String   -- ^ new indices-       -> NArray i t-t ! ns = rename t (map return ns)---- | Rename indices. Equal indices are contracted out.-rename :: (Coord t, Compat i)+-- | Rename indices (in the internal order). Equal indices are contracted out.+renameRaw :: (Coord t, Compat i)        => NArray i t        -> [Name]     -- ^ new names        -> NArray i t-rename t ns = reorder orig (contract t')-    where t' = renameRaw t ns-          orig = nub (names t') \\ common1 t'-+renameRaw t ns = contract (renameSuperRaw t ns) -renameRaw (A d v) l | length l == length d = A d' v-                    | otherwise = error $ "rename " ++ show d ++ " with " ++ show l-    where d' = zipWith f d l-          f i n = i {iName=n}+renameSuperRaw (A d v) l+    | length l == length d = A d' v+    | otherwise = error $ "renameRaw " ++ show d ++ " with " ++ show l+   where d' = zipWith f d l+         f i n = i {iName=n}  mapDims f (A d v) = A (map f d) v -mapTypes :: (i1 -> i) -> NArray i1 t -> NArray i t+mapTypes :: (i1 -> i2) -> NArray i1 t -> NArray i2 t mapTypes f = mapDims (\i -> i {iType = f (iType i)}) --- mapNames f = mapDims (\i -> i {iName = f (iName i)})+mapNames :: (Name -> Name) -> NArray i t -> NArray i t+mapNames f = mapDims (\i -> i {iName = f (iName i)}) +-- | Rename indices using an association list.+renameExplicit :: (Compat i, Coord t) => [(Name,Name)] -> NArray i t -> NArray i t+renameExplicit al = g . mapNames f where+    f n = maybe n id (lookup n al)+    g t = reorder orig (contract t) where orig = nub (namesR t) \\ common1 t++-- | Index names (in internal order).+namesR :: NArray i t -> [Name]+namesR = map iName . dims+ -- | Index names. names :: NArray i t -> [Name]-names = map iName . dims+names = sort . namesR  -- | Dimension of given index. size :: Name -> NArray i t -> Int size n t = (iDim . head) (filter ((n==).iName) (dims t)) +sizesR :: NArray i t -> [Int]+sizesR = map iDim . dims++sizes :: NArray i t -> [Int]+sizes t = map (flip size t) (names t)+ -- | Type of given index. typeOf :: Compat i => Name -> NArray i t -> i typeOf n t = (iType . head) (filter ((n==).iName) (dims t)) - -- | The number of dimensions of a multidimensional array.-rank :: NArray i t -> Int-rank = length . dims+order :: NArray i t -> Int+order = length . dims +selDims ds = map f where+    f n = head $ filter ((n==).iName) ds+ ---------------------------------------------------------- +common2 t1 t2 = [ n1 | n1 <- namesR t1, n2 <- namesR t2, n1==n2]++analyzeProduct :: (Coord t, Compat i) => NArray i t -> NArray i t -> Maybe (NArray i t, Int)+analyzeProduct a b = r where+    nx  = common2 a b+    dx1 = selDims (dims a) nx+    dx2 = selDims (dims b) nx+    ok  = and $ zipWith compat dx1 dx2+    (tma,na) = matrixatorFree a nx+    (mb,nb)  = matrixatorFree b nx+    mc  = trans tma <> mb+    da  = selDims (dims a) na+    db  = selDims (dims b) nb+    dc  = da ++ db+    c   = A dc (flatten mc)+    sz  = product (map iDim dc)+    r | ok = Just (c, sz)+      | otherwise = Nothing++infixl 5 |*|+-- | Tensor product with automatic contraction of repeated indices, following Einstein summation convention.+(|*|) :: (Coord t, Compat i) => NArray i t -> NArray i t -> NArray i t+t1 |*| t2 = case analyzeProduct t1 t2 of+    Nothing -> error $ "wrong contraction2: "++(show $ dims t1)++" and "++(show $ dims t2)+    Just (r,_) -> r++----------------------------------------------------------+ lastIdx name t = ((d1,d2),m) where     (d1,d2) = span (\d -> iName d /= name) (dims t)     c = product (map iDim d2)@@ -165,15 +204,37 @@           nd = d2++d1           c = dim (coords t) `div` (iDim $ head d2) +-- | Obtain a matrix whose columns are the fibers of the array in the given dimension. The column order depends on the selected index (see 'matrixator').+fibers :: Coord t => Name -> NArray i t -> Matrix t+fibers n = snd . firstIdx n +-- | Reshapes an array as a matrix with the desired dimensions as flattened rows and flattened columns.+matrixator :: (Coord t) => NArray i t -- ^ input array+                        -> [Name]    -- ^ row dimensions+                        -> [Name]    -- ^ column dimensions+                        -> Matrix t   -- ^ result+matrixator t nr nc = reshape s (coords q) where+    q = reorder (nr++nc) t+    s = product (map (flip size t) nc)++-- | Reshapes an array as a matrix with the desired dimensions as flattened rows and flattened columns. We do not force the order of the columns.+matrixatorFree :: (Coord t)+               => NArray i t          -- ^ input array+               -> [Name]              -- ^ row dimensions+               -> (Matrix t, [Name])  -- ^ (result, column dimensions)+matrixatorFree t nr = (reshape s (coords q), nc) where+    q = tridx nr t+    nc = drop (length nr) (map iName (dims q))+    s = product (map (flip size t) nc)+ -- | Create a list of the substructures at the given level. parts :: (Coord t)        => NArray i t       -> Name        -- ^ index to expand       -> [NArray i t]-parts a name | name `elem` (names a) = map (reorder orig) (partsRaw a name)-             | otherwise = error $ "parts: " ++ show name ++ " is not a dimension of "++(show $ names a)-    where orig = names a \\ [name]+parts a name | name `elem` (namesR a) = map (reorder orig) (partsRaw a name)+             | otherwise = error $ "parts: " ++ show name ++ " is not a dimension of "++(show $ namesR a)+    where orig = namesR a \\ [name]  partsRaw a name = map f (toRows m)     where (_:ds,m) = firstIdx name a@@ -191,22 +252,10 @@ -- | Change the internal layout of coordinates. -- The array, considered as an abstract object, does not change. reorder :: (Coord t) => [Name] -> NArray i t -> NArray i t-reorder ns b | ns == names b = b-             | sort ns == sort (names b) = tridx ns b+reorder ns b | ns == namesR b = b+             | sort ns == sort (namesR b) = tridx ns b              | otherwise = error $ "wrong index sequence " ++ show ns-                                    ++ " to reorder "++(show $ names b)----- | 'reorder' (transpose) the dimensions of the array (with single letter names).------ Operations are defined by named indices, so the transposed array is operationally equivalent to the original one.-infixl 8 ~>-(~>) :: (Coord t) => NArray i t -> String -> NArray i t-t ~> ns = reorder (map return ns) t-----------------------------------------------------------------rawProduct (A d1 v1) (A d2 v2) = A (d1++d2) (flatten (outer v1 v2))+                                    ++ " to reorder "++(show $ namesR b)  ---------------------------------------------------------------------- @@ -222,6 +271,7 @@ -- | Class of compatible indices for contractions. class (Eq a, Show (Idx a)) => Compat a where     compat :: Idx a -> Idx a -> Bool+    opos   :: Idx a -> Idx a   @@ -241,38 +291,15 @@  contract1c t n = contract1 renamed n n'     where n' = " "++n++" " -- forbid spaces in names...-          renamed = renameRaw (t) auxnames+          renamed = renameSuperRaw (t) auxnames           auxnames = h ++ (n':r)-          (h,_:r) = break (==n) (names t)+          (h,_:r) = break (==n) (namesR t)  common1 t = [ n1 | (a,n1) <- x , (b,n2) <- x, a>b, n1==n2]-    where x = zip [0 ::Int ..] (names t)+    where x = zip [0 ::Int ..] (namesR t)  contract t = foldl' contract1c t (common1 t) -------------------------------------------------------------------------contract2 t1 t2 n | ok = A (tail ds1 ++ tail ds2) (flatten m)-                  | otherwise = error $ "wrong contraction2: "++ n ++ " of "++-                                      (show $ dims t1)++" and "++ (show $ dims t2)-  where ok = (compat <$> getName t1 n <*> getName t2 n) == Just True-        (ds1,m1) = firstIdx n t1-        (ds2,m2) = firstIdx n t2-        m = (trans m1) `multiply` m2--common2 t1 t2 = [ n1 | n1 <- names t1, n2 <- names t2, n1==n2]--infixl 5 |*|--- | Tensor product with automatic contraction of repeated indices, following Einstein summation convention.-(|*|) :: (Coord t, Compat i)-      => NArray i t -> NArray i t -> NArray i t-t1 |*| t2 = r where-    cs = common2 t1 t2-    r = case cs of-        [] -> rawProduct t1 t2-        n:_ -> reorder orig $ contract (contract2 t1 t2 n)-    orig = nub (names t1 ++ names t2) \\ cs- -------------------------------------------------------------  -- | Check if two arrays have the same structure.@@ -292,145 +319,6 @@  ------------------------------------------------------- -showBases x = f $ concatMap (shbld) x-    where   shbld (c,[]) = shsign c ++ showc c-            shbld (c,l) = shsign c ++ g (showc c) ++ "{"++ concatMap show l++"}"-            shsign c = if c < 0 then " - " else " + "-            showc c-                | abs (fromIntegral (round c :: Int) - c) <1E-10  = show (round $ abs c::Int)-                | otherwise = printf "%.3f" (abs c)-            f (' ':'+':' ':rs) = rs-            f (' ':'-':' ':rs) = '-':rs-            f a = a-            g "1" = ""-            g a = a-------------------------------------------------------------data Rect = Rect { li :: Int, co :: Int, els :: [String] }--rect s = pad r c (Rect r 0 ss)-    where ss = lines s-          r  = length ss-          c  = maximum (map length ss)--pad nr nc (Rect r c ss) = Rect (r+r') (c+c') ss'' where-    r' = max 0 (nr-r)-    c' = max 0 (nc-c)-    ss' = map (padH nc) ss-    ss'' = replicate r' (replicate nc '-') ++ ss'-    padH l s = take (l-length s) (" | "++repeat ' ') ++ s--dispH :: Int -> [Rect] -> Rect-dispH k rs = Rect nr nc nss where-    nr = maximum (map li rs)-    nss' = mapTail (\x-> pad nr (co x + k) x) rs-    nss = foldl1' (zipWith (++)) (map els nss')-    nc = length (head nss)--dispV :: Int -> [Rect] -> Rect-dispV k rs = Rect nr nc nss where-    nc = maximum (map co rs)-    nss' = mapTail (\x-> pad (li x + k) nc x) rs-    nss = concatMap els nss'-    nr = length nss--mapTail f (a:b) = a : map f b-mapTail _ x     = x-----formatAux f x = unlines . addds . els . fmt ms $ x where-    fmt [] _ = undefined -- cannot happen-    fmt (g:gs) t-        | rank t == 0 = rect (f (coords t @> 0))-        | rank t == 1 =  rect $ unwords $ map f (toList $ coords t)-        | rank t == 2 =  decor t $ rect $ w1 $ format " " f (reshape (iDim $ last $ dims t) (coords t))-        | otherwise    = decor t (g ps)-      where ps = map (fmt gs ) (partsRaw t (head (names t)))-    ds = showNice (filter ((/='*').head.iName) $ dims x)-    addds = if null ds then (showRawDims (dims x) :) else (ds:)-    w1 = unlines . map (' ':) . lines-    ms = cycle [dispV 1, dispH 2]-    decor t | odd (rank t) = id-            | otherwise = decorLeft  (names t!!0) . decorUp (names t!!1)---showNice x = unwords . intersperse "x" . map show $ x-showRawDims = showNice . map iDim . filter ((/="*").iName)------------------------------------------------------------ | Show a multidimensional array as a nested 2D table.-formatArray :: (Coord t, Compat i)-      => (t -> String) -- ^ format function (eg. printf \"5.2f\")-      -> NArray i t-      -> String-formatArray f t | odd (rank t) = formatAux f (dummyAt 0 t)-            | otherwise    = formatAux f t---decorUp s rec-    | head s == '*' = rec-    | otherwise     = dispV 0 [rs,rec]-  where-    c = co rec-    c1 = (c - length s) `div` 2-    c2 = c - length s - c1-    rs = rect $ replicate c1 ' ' ++ s ++ replicate c2 ' '--decorLeft s rec-    | head s == '*' = rec-    | otherwise     = dispH 0 [rs,rec]-  where-    c = li rec-    r1 = (c - length s+1) `div` 2-    r2 = c - length s - r1-    rs = rect $ unlines $ replicate r1 spc ++ s : replicate (r2) spc-    spc = replicate (length s) ' '------------------------------------------------------------ | Print the array as a nested table with the desired format (e.g. %7.2f) (see also 'formatArray', and 'formatScaled').-printA :: (Coord t, Compat i, PrintfArg t) => String -> NArray i t -> IO ()-printA f t = putStrLn (formatArray (printf f) t)----- | Show the array as a nested table with autoscaled entries.-formatScaled :: (Compat i)-      => Int -- ^ number of of decimal places-      -> NArray i Double-      -> String-formatScaled dec t = unlines (('(':d++")  E"++show o) : m)-    where ss = formatArray (printf fmt. g) t-          d:m = lines ss-          g x = x/10^(o::Int)-          o = floor $ maximum $ map (logBase 10 . abs) $ toList $ coords t-          fmt = '%':show (dec+3) ++ '.':show dec ++"f"---- | Show the array as a nested table with a \"\%.nf\" format. If all entries--- are approximate integers the array is shown without the .00.. digits.-formatFixed :: (Compat i)-      => Int -- ^ number of of decimal places-      -> NArray i Double-      -> String-formatFixed dec t-    | isInt t   = formatArray (printf ('%': show (width t) ++".0f")) t-    | otherwise = formatArray (printf ('%': show (width t+dec+1) ++"."++show dec ++"f")) t--isInt = all lookslikeInt . toList . coords-lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx-    where shx = show x--- needsSign t = vectorMin (coords t) < 0--- width :: Compat i => NArray i Double -> Int-width = maximum . map (length . (printf "%.0f"::Double->String)) . toList . coords--- width t = k + floor (logBase 10 (max 1 $ vectorMax (abs $ coords t))) :: Int---      where k | needsSign t = 2---              | otherwise = 1--------------------------------------------------------- -- | Create an array from a list of subarrays. (The inverse of 'parts'.) newIndex:: (Coord t, Compat i) =>      i  -- ^ index type@@ -438,21 +326,11 @@      -> [NArray i t]      -> NArray i t newIndex i name ts = r where-    ds = Idx (length ts) name i : (dims (head cts))+    ds = Idx i (length ts) name : (dims (head cts))     cts = makeConformant ts     r = mkNArray ds (join $ map coords cts) ---- | Insert a dummy index of dimension 1 at a given level (for formatting purposes).-dummyAt :: Int -> NArray i t -> NArray i t-dummyAt k t = mkNArray d' (coords t) where-    (d1,d2) = splitAt k (dims t)-    d' = d1 ++ d : d2-    d = Idx 1 "*" undefined---- | Rename indices so that they are not shown in formatted output.-noIdx :: Compat i => NArray i t -> NArray i t-noIdx t = renameRaw t (map ('*':) (names t))+-------------------------------------------------------  -- | Obtain a canonical base for the array. basisOf :: Coord t => NArray i t -> [NArray i t]@@ -469,6 +347,10 @@     real = mapArray real     complex = mapArray complex ++-- instance (NFData t, Element t) => NFData (NArray i t) where+--     rnf = rnf . coords+ ----------------------------------------------------------------------  -- | obtains the common value of a property of a list@@ -484,30 +366,31 @@ -- | Extract the 'Matrix' corresponding to a two-dimensional array, -- in the rows,cols order. asMatrix :: (Coord t) => NArray i t -> Matrix t-asMatrix a | rank a == 2 = reshape c (coords a)-           | otherwise = error $ "asMatrix requires a rank 2 array."-    where c = size (last (names a)) a+asMatrix a | order a == 2 = reshape c (coords a')+           | otherwise = error $ "asMatrix requires a 2nd order array."+    where c = size (last (namesR a')) a'+          a' = reorder (sort (namesR a)) a  -- | Extract the 'Vector' corresponding to a one-dimensional array. asVector :: (Coord t) => NArray i t -> Vector t-asVector a | rank a == 1 = coords a-           | otherwise = error $ "asVector requires a rank 1 array."+asVector a | order a == 1 = coords a+           | otherwise = error $ "asVector requires a 1st order array."  -- | Extract the scalar element corresponding to a 0-dimensional array. asScalar :: (Coord t) => NArray i t -> t-asScalar a | rank a == 0 = coords a @>0-           | otherwise = error $ "asScalar requires a rank 0 array."+asScalar a | order a == 0 = coords a @>0+           | otherwise = error $ "asScalar requires a 0th order array."  ------------------------------------------------------------------------ --- | Create a rank-1 array from an hmatrix 'Vector'.+-- | Create a 1st order array from a 'Vector'. fromVector :: Compat i => i -> Vector t -> NArray i t-fromVector i v = mkNArray [Idx (dim v) "1" i ] v+fromVector i v = mkNArray [Idx i (dim v) "1"] v --- | Create a rank-2 array from an hmatrix 'Matrix'.+-- | Create a 2nd order array from a 'Matrix'. fromMatrix :: (Compat i, Coord t) => i -> i -> Matrix t -> NArray i t-fromMatrix ir ic m = mkNArray [Idx (rows m) "1" ir,-                               Idx (cols m) "2" ic] (flatten m)+fromMatrix ir ic m = mkNArray [Idx ir (rows m) "1",+                               Idx ic (cols m) "2"] (flatten m)  ------------------------------------------------------------------------ @@ -517,7 +400,7 @@         -> Name         -> NArray i t         -> NArray i t-extract f name arr = reorder (names arr)+extract f name arr = reorder (namesR arr)                    . newIndex (typeOf name arr) name                    . map snd . filter (uncurry f)                    $ zip [1..] (parts arr name)@@ -528,7 +411,11 @@      -> Name      -> NArray i a      -> NArray i b-onIndex f name t = reorder (names t) $ newIndex (typeOf name t) name (f (parts t name))+onIndex f name t = r where+     r = if sort (namesR x) == sort (namesR t)+            then reorder (namesR t) x+            else x+     x = newIndex (typeOf name t) name (f (parts t name))  ------------------------------------------------------------------------ @@ -560,3 +447,45 @@         Just alldims -> (extend alldims t1, extend alldims t2)         Nothing -> error $ "makeConformantT with inconsistent dimensions "                          ++ show (dims t1, dims t2)++---------------------------------------------++takeDiagT :: (Compat i, Coord t) => NArray i t -> [t]+takeDiagT t = map (asScalar . atT t) cds where+    n = minimum (sizesR t)+    o = order t+    cds = map (replicate o) [0..n-1]++atT :: (Compat i, Coord t) => NArray i t -> [Int] -> NArray i t+atT t c = atT' c t where+    atT' cs = foldl1' (.) (map fpart cs)+    fpart k q = parts q (head (namesR q)) !! k++----------------------------------------------++-- not very smart...++-- | This is equivalent to the regular 'product', but in the order that minimizes the size of the+-- intermediate factors.+smartProduct :: (Coord t, Compat i, Num (NArray i t)) => [NArray i t] -> NArray i t+smartProduct [] = 1+smartProduct [a] = a+smartProduct [a,b] = a*b+smartProduct ts = r where+    n = length ts+    ks = [0 .. n-1]+    xs = zip ks ts+    g a b = case analyzeProduct a b of+              Nothing -> error $ "inconsistent dimensions in smartProduct: "++(show $ dims a)++" and "++(show $ dims b)+              Just (_,c) -> c+    pairs = [ ((i,j), g a b) | (i,a) <- init xs, (j,b) <- drop (i+1) xs ]+    (p,q) = fst $ minimumBy (compare `on` snd) pairs+    r = smartProduct (ts!!p * ts!!q : (dropElemPos p . dropElemPos q) ts)++dropElemPos k xs = take k xs ++ drop (k+1) xs++----------------------------------------------++-- | sequence of n indices with given prefix+seqIdx :: Int -> String -> [Name]+seqIdx n prefix = [prefix ++ show k | k <- [1 .. n] ]

@@ -23,18 +23,20 @@  import Numeric.LinearAlgebra.Array.Internal import Data.Packed+import Data.List(intersperse)   instance Show (Idx None) where-    show (Idx n s _t) = show n ++ ":" ++ s+    show (Idx _t n s) = s ++ ":" ++ show n  -- | Unespecified coordinate type. Contractions only -- require equal dimension.-data None = None deriving Eq+data None = None deriving (Eq,Show)   instance Compat None where     compat d1 d2 = iDim d1 == iDim d2+    opos = id   -- | Multidimensional array with unespecified coordinate type.@@ -42,14 +44,19 @@  instance (Coord t) => Show (Array t) where     show t | null (dims t) = "scalar "++ show (coords t @>0)-           | otherwise = "listArray "++ show (dims t) ++ " "++ show (toList $ coords t)+           | order t == 1 = "index " ++ show n ++" " ++ (show . toList . coords $ t)+           | otherwise = "index "++ show n ++ " [" ++ ps ++ "]"+      where n = head (namesR t)+            ps = concat $ intersperse ", " $ map show (parts t n) +-- ++ " "++ show (toList $ coords t)+ -- | Construction of an 'Array' from a list of dimensions and a list of elements in left to right order. listArray :: (Coord t)     => [Int] -- ^ dimensions     -> [t]   -- ^ elements     -> Array t listArray ds cs = mkNArray dms (product ds |> (cs ++ repeat 0))-    where dms = zipWith3 Idx ds (map show [1::Int ..]) (repeat None)+    where dms = zipWith3 Idx (repeat None) ds (map show [1::Int ..])  

@@ -0,0 +1,314 @@+{-# LANGUAGE FlexibleContexts #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Packed.Array.Solve+-- Copyright   :  (c) Alberto Ruiz 2009+-- License     :  GPL+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+--+-- Solution of general multidimensional linear and multilinear systems.+--+-----------------------------------------------------------------------------++module Numeric.LinearAlgebra.Array.Solve (+-- * Linear systems+    solve,+    solveHomog, solveHomog1, solveH,+    solveP,+-- *  Multilinear systems+-- ** General+    ALSParam(..), defaultParameters,+    mlSolve, mlSolveH, mlSolveP,+-- ** Factorized+    solveFactors, solveFactorsH,+-- * Utilities+    eps, eqnorm, infoRank,+    solve', solveHomog', solveHomog1', solveP'+) where++import Numeric.LinearAlgebra.Array.Util+import Numeric.LinearAlgebra.Exterior+import Numeric.LinearAlgebra.Array.Internal(mkNArray, selDims, debug, namesR)+import Numeric.LinearAlgebra hiding ((.*), scalar)+import Data.List+import System.Random+++-- | Solution of the linear system a x = b, where a and b are+-- general multidimensional arrays. The structure and dimension names+-- of the result are inferred from the arguments.+solve :: (Compat i, Coord t)+        => NArray i t -- ^ coefficients (a)+        -> NArray i t -- ^ target       (b)+        -> NArray i t -- ^ result       (x)+solve = solve' id++solve' g a b = x where+    nx = namesR a \\ namesR b+    na = namesR a \\ nx+    nb = namesR b \\ namesR a+    aM = g $ matrixator a na nx+    bM = matrixator b na nb+    xM = linearSolveSVD aM bM+    dx = map opos (selDims (dims a) nx) ++ selDims (dims b) nb+    x = mkNArray dx (flatten xM)+++-- | Solution of the homogeneous linear system a x = 0, where a is a+-- general multidimensional array.+--+-- If the system is overconstrained we may provide the theoretical rank to get a MSE solution.+solveHomog :: (Compat i, Coord t)+           =>  NArray i t    -- ^ coefficients (a)+           -> [Name]         -- ^ desired dimensions for the result+                             --   (a subset selected from the target).+           -> Either Double Int -- ^ Left \"numeric zero\" (e.g. eps), Right \"theoretical\" rank+           -> [NArray i t] -- ^ basis for the solutions (x)+solveHomog = solveHomog' id++solveHomog' g a nx' hint = xs where+    nx = filter (`elem` (namesR a)) nx'+    na = namesR a \\ nx+    aM = g $ matrixator a na nx+    vs = nullspaceSVD hint aM (rightSV aM)+    dx = map opos (selDims (dims a) nx)+    xs = map (mkNArray dx) vs++-- | A simpler way to use 'solveHomog', which returns just one solution.+-- If the system is overconstrained it returns the MSE solution.+solveHomog1 :: (Compat i, Coord t)+            => NArray i t+            -> [Name]+            -> NArray i t+solveHomog1 = solveHomog1' id++solveHomog1' g m ns = head $ solveHomog' g m ns (Right (k-1))+    where k = product $ map iDim $ selDims (dims m) ns++-- | 'solveHomog1' for single letter index names.+solveH :: (Compat i, Coord t) => NArray i t -> [Char] -> NArray i t+solveH m ns = solveHomog1 m (map return ns)+++-- | Solution of the linear system a x = b, where a and b are+-- general multidimensional arrays, with homogeneous equality along a given index.+solveP :: Tensor Double   -- ^ coefficients (a)+       -> Tensor Double   -- ^ desired result (b)+       -> Name            -- ^ the homogeneous dimension+       -> Tensor Double   -- ^ result (x)+solveP = solveP' id++solveP' g a b h = mapTat (solveP1 g h a) (namesR b \\ (h:namesR a)) b++-- solveP for a single right hand side+solveP1 g nh a b = solveHomog1' g ou ns where+    k = size nh b+    epsi = t $ leviCivita k `renameO` (nh : (take (k-1) $ (map (('e':).(:[])) ['2'..])))+    ou = a .* b' * epsi+    ns = (namesR a \\ namesR b) ++ x+    b' = renameExplicit [(nh,"e2")] b+    x = if nh `elem` (namesR a) then [] else [nh]+    t = if typeOf nh b == Co then contrav else cov+        -- mapTypes (const (opos $ typeOf nh b))++-----------------------------------------------------------------------++-- | optimization parameters for alternating least squares+data ALSParam i t = ALSParam+    { nMax  ::   Int     -- ^ maximum number of iterations+    , delta ::   Double  -- ^ minimum relative improvement in the optimization (percent, e.g. 0.1)+    , epsilon :: Double  -- ^ maximum relative error. For nonhomogeneous problems it is+                         --   the reconstruction error in percent (e.g.+                         --   1E-3), and for homogeneous problems is the frobenius norm of the+                         --  expected zero structure in the right hand side.+    , post :: [NArray i t] -> [NArray i t] -- ^ post-processing function after each full iteration (e.g. 'id')+    , postk :: Int -> NArray i t -> NArray i t-- ^ post-processing function for the k-th argument (e.g. 'const' 'id')+    , presys :: Matrix t -> Matrix t -- ^ preprocessing function for the linear systems (eg. 'id', or 'infoRank')+    }+++optimize :: (x -> x)      -- ^ method+         -> (x -> Double) -- ^ error function+         -> x             -- ^ starting point+         -> ALSParam i t     -- ^ optimization parameters+         -> (x, [Double]) -- ^ solution and error history+optimize method errfun s0 p = (sol,e) where+    sols = take (max 1 (nMax p)) $ iterate method s0+    errs = map errfun sols+    (sol,e) = convergence (zip sols errs) []+    convergence [] _  = error "impossible"+    convergence [(s,err)] prev = (s, err:prev)+    convergence ((s1,e1):(s2,e2):ses) prev+        | e1 < epsilon p = (s1, e1:prev)+        | abs (100*(e1 - e2)/e1) < delta p = (s2, e2:prev)+        | otherwise = convergence ((s2,e2):ses) (e1:prev)++percent t s = 100 * frobT (t - smartProduct s) / frobT t++percentP h t s = 100 * frobT (t' - s') / frobT t' where+    t' = f t+    s' = f (smartProduct s)+    f = mapTat g (namesR t \\ [h])+    g v = v / atT v [n]+    n = size h t - 1++frobT t = pnorm PNorm2 . coords $ t+--unitT t = t / scalar (frobT t)++dropElemPos k xs = take k xs ++ drop (k+1) xs+replaceElemPos k v xs = take k xs ++ v : drop (k+1) xs++takes [] _ = []+takes (n:ns) xs = take n xs : takes ns (drop n xs)++----------------------------------------------------------------------++alsStep f params a x = (foldl1' (.) (map (f params a) [n,n-1 .. 0])) x+    where n = length x - 1++-----------------------------------------------------------------------++-- | Solution of a multilinear system a x y z ... = b based on alternating least squares.+mlSolve+  :: (Compat i, Coord t, Num (NArray i t), Normed (Vector t))+     => ALSParam i t     -- ^ optimization parameters+     -> [NArray i t]  -- ^ coefficients (a), given as a list of factors.+     -> [NArray i t]  -- ^ initial solution [x,y,z...]+     -> NArray i t    -- ^ target (b)+     -> ([NArray i t], [Double]) -- ^ Solution and error history+mlSolve params a x0 b+    = optimize (post params . alsStep (alsArg b) params a) (percent b . (a++)) x0 params++alsArg _ _ _ _ [] = error "alsArg _ _ []"+alsArg b params a k xs = sol where+    p = smartProduct (a ++ dropElemPos k xs)+    x = solve' (presys params) p b+    x' = postk params k x+    sol = replaceElemPos k x' xs++----------------------------------------------------------++-- | Solution of the homogeneous multilinear system a x y z ... = 0 based on alternating least squares.+mlSolveH+  :: (Compat i, Coord t, Num (NArray i t), Normed (Vector t))+     => ALSParam  i t    -- ^ optimization parameters+     -> [NArray i t]  -- ^ coefficients (a), given as a list of factors.+     -> [NArray i t]  -- ^ initial solution [x,y,z...]+     -> ([NArray i t], [Double]) -- ^ Solution and error history+mlSolveH params a x0+    = optimize (post params . alsStep alsArgH params a) (frobT . smartProduct . (a++)) x0 params++alsArgH _ _ _ [] = error "alsArgH _ _ []"+alsArgH params a k xs = sol where+    p = smartProduct (a ++ dropElemPos k xs)+    x = solveHomog1' (presys params) p (namesR (xs!!k))+    x' = postk params k x+    sol = replaceElemPos k x' xs++----------------------------------------------------------++-- | Solution of a multilinear system a x y z ... = b, with a homogeneous index, based on alternating least squares.+mlSolveP+     :: ALSParam Variant Double     -- ^ optimization parameters+     -> [Tensor Double]  -- ^ coefficients (a), given as a list of factors.+     -> [Tensor Double]  -- ^ initial solution [x,y,z...]+     -> Tensor Double    -- ^ target (b)+     -> Name             -- ^ homogeneous index+     -> ([Tensor Double], [Double]) -- ^ Solution and error history+mlSolveP params a x0 b h+    = optimize (post params . alsStep (alsArgP b h) params a) (percentP h b . (a++)) x0 params++alsArgP _ _ _ _ _ [] = error "alsArgP _ _ []"+alsArgP b h params a k xs = sol where+    p = smartProduct (a ++ dropElemPos k xs)+    x = solveP' (presys params)  p b h+    x' = postk params k x+    sol = replaceElemPos k x' xs++-------------------------------------------------------------++{- | Given two arrays a (source) and  b (target), we try to compute linear transformations x,y,z,... for each dimension, such that product [a,x,y,z,...] == b.+(We can use 'eqnorm' for 'post' processing, or 'id'.)+-}+solveFactors :: (Coord t, Random t, Compat i, Num (NArray i t), Normed (Vector t))+             => Int          -- ^ seed for random initialization+             -> ALSParam i t     -- ^ optimization parameters+             -> [NArray i t] -- ^ source (also factorized)+             -> String       -- ^ index pairs for the factors separated by spaces+             -> NArray i t   -- ^ target+             -> ([NArray i t],[Double]) -- ^ solution and error history+solveFactors seed params a pairs b =+    mlSolve params a (initFactorsRandom seed (smartProduct a) pairs b) b++initFactorsSeq rs a pairs b | ok = as+                            | otherwise = error "solveFactors index pairs"+  where+    (ia,ib) = unzip (map (\[x,y]->([x],[y])) (words pairs))+    ic = intersect (namesR a) (namesR b)+    ok = sort (namesR b\\ic) == sort ib && sort (namesR a\\ic) == sort ia+    db = selDims (dims b) ib+    da = selDims (dims a) ia+    nb = map iDim db+    na = map iDim da+    ts = takes (zipWith (*) nb na) rs+    as = zipWith5 f ts ib ia db da+    f c i1 i2 d1 d2 = (mkNArray [d1,opos d2] (fromList c)) `renameO` [i1,i2]++initFactorsRandom seed a b = initFactorsSeq (randomRs (-1,1) (mkStdGen seed)) a b+++-- | Homogeneous factorized system. Given an array a,+-- given as a list of factors as, and a list of pairs of indices+-- [\"pi\",\"qj\", \"rk\", etc.], we try to compute linear transformations+-- x!\"pi\", y!\"pi\", z!\"rk\", etc. such that product [a,x,y,z,...] == 0.+solveFactorsH+  :: (Coord t, Random t, Compat i, Num (NArray i t), Normed (Vector t))+     => Int -- ^ seed for random initialization+     -> ALSParam  i t    -- ^ optimization parameters+     -> [NArray i t] -- ^ coefficient array (a), (also factorized)+     -> String       -- ^ index pairs for the factors separated by spaces+     -> ([NArray i t], [Double]) -- ^ solution and error history+solveFactorsH seed params a pairs =+    mlSolveH params a (initFactorsHRandom seed (smartProduct a) pairs)++initFactorsHSeq rs a pairs = as where+    (ir,it) = unzip (map (\[x,y]->([x],[y])) (words pairs))+    nr = map (flip size a) ir+    nt = map (flip size a) it+    ts = takes (zipWith (*) nr nt) rs+    as = zipWith5 f ts ir it (selDims (dims a) ir) (selDims (dims a) it)+    f c i1 i2 d1 d2 = (mkNArray (map opos [d1,d2]) (fromList c)) `renameO` [i1,i2]++initFactorsHRandom seed a pairs = initFactorsHSeq (randomRs (-1,1) (mkStdGen seed)) a pairs++----------------------------------++-- | post processing function that modifies a list of tensors so that they+-- have equal frobenius norm+eqnorm :: (Coord t, Coord (Complex t), Compat i, Num (NArray i t), Normed (Vector t) )+       => [NArray i t] -> [NArray i t]++eqnorm [] = error "eqnorm []"+eqnorm as = as' where+    n = length as+    fs = map (frobT) as+    s = product fs ** (1/fromIntegral n)+    as' = zipWith g as fs where g a f = a * real (scalar (s/f))++-- | nMax = 20, epsilon = 1E-3, delta = 1, post = id, postk = const id, presys = id+defaultParameters :: ALSParam i t+defaultParameters = ALSParam {+    nMax = 20,+    epsilon = 1E-3,+    delta = 1,+    post = id,+    postk = const id,+    presys = id+  }++-- | debugging function (e.g. for 'presys'), which shows rows, columns and rank of the+-- coefficient matrix of a linear system.+infoRank :: Field t => Matrix t -> Matrix t+infoRank a = debug "" (const (rows a, cols a, rank a)) a

@@ -17,16 +17,16 @@     Coord, Compat(..),     NArray, Idx(..), Name,     scalar,-    rank, names, size, typeOf, dims, coords,+    order, names, size, sizes, typeOf, dims, coords, -    rename, (!),+    renameExplicit, (!>), renameO, (!),      parts,     newIndex, -    mapArray, zipArray, (|*|),+    mapArray, zipArray, (|*|), smartProduct, outers, -    extract, onIndex,+    extract, onIndex, mapTat,      reorder, (~>),     formatArray, formatFixed, formatScaled,@@ -35,10 +35,115 @@     sameStructure,     makeConformant,     basisOf,-    asScalar, asVector, asMatrix,+    atT, takeDiagT, diagT,+    mkFun, mkAssoc, setType,+    renameParts,+    asScalar, asVector, asMatrix, applyAsMatrix,+    fibers, matrixator, matrixatorFree, analyzeProduct,     fromVector, fromMatrix,     Container(..), ) where  import Numeric.LinearAlgebra.Array.Internal-import Data.Packed(Container(..))+import Numeric.LinearAlgebra.Array.Display+import Data.Packed(Container(..),Matrix)+import Numeric.LinearAlgebra.Array.Simple+import Data.List(intersperse,sort,foldl1')++-- infixl 9 #+-- (#) :: [Int] -> [Double] -> Array Double+-- (#) = listArray++-- | Multidimensional diagonal of given order.+diagT :: [Double] -> Int -> Array Double+diagT v n = replicate n k `listArray` concat (intersperse z (map return v))+    where k = length v+          tot = k^n+          nzeros = (tot - k) `div` (k-1)+          z = replicate nzeros 0+++-- | Explicit renaming of single letter index names.+--+-- For instance, @t >\@> \"pi qj\"@ changes index \"p\" to \"i\" and \"q\" to \"j\".+(!>) :: (Compat i, Coord t) => NArray i t -> [Char] -> NArray i t+infixl 9 !>+t !> s = renameExplicit (map (\[a,b]->([a],[b])) (words s)) t+++-- | Rename indices in alphabetical order. Equal indices of compatible type are contracted out.+renameO :: (Coord t, Compat i)+    => NArray i t+    -> [Name]+    -> NArray i t+renameO t ns = renameExplicit (zip od ns) t+    where od = map iName (sort (dims t))+++-- | Rename indices in alphabetical order ('renameO') using single letter names.+(!) :: (Compat i, Coord t) => NArray i t -> [Char] -> NArray i t+infixl 9 !+t ! s = renameExplicit (zip od (map return s)) t+    where od = map iName (sort (dims t))+++-- -- | 'renameRaw' the indices (in the internal order) with single-letter names. Equal indices of compatible type are contracted out.+-- infixl 8 !!!+-- (!!!) :: (Coord t, Compat i)+--        => NArray i t+--        -> String   -- ^ new indices+--        -> NArray i t+-- t !!! ns = renameRaw t (map return ns)+++-- | 'reorder' (transpose) dimensions of the array (with single letter names).+--+-- Operations are defined by named indices, so the transposed array is operationally equivalent to the original one.+infixl 8 ~>+(~>) :: (Coord t) => NArray i t -> String -> NArray i t+t ~> ns = reorder (map return ns) t+++-- | Map a function at the internal level selected by a set of indices+mapTat :: (Coord a, Coord b, Compat i)+         => (NArray i a -> NArray i b)+         -> [Name]+         -> NArray i a+         -> NArray i b+mapTat f [] = f+mapTat f (a:as) = onIndex (map $ mapTat f as) a++-- | Outer product of a list of arrays along the common indices.+outers :: (Coord a, Compat i) => [NArray i a] -> NArray i a+outers = foldl1' (zipArray (*))++-- | Define an array using a function.+mkFun :: [Int] -> ([Int] -> Double) -> Array Double+mkFun ds f = listArray ds $ map f (sequence $ map (enumFromTo 0 . subtract 1. fromIntegral) $ ds)++-- | Define an array using an association list.+mkAssoc :: [Int] -> [([Int], Double)] -> Array Double+mkAssoc ds ps = mkFun ds f where+   f = maybe 0 id . flip lookup ps++-- | Change type of index.+setType :: (Compat i, Coord t) => Name -> i -> NArray i t -> NArray i t+setType n t a = mapDims f a where+    f i | iName i == n = i {iType = t}+        | otherwise    = i++-- | Extract the 'parts' of an array, and renameRaw one of the remaining indices+-- with succesive integers.+renameParts :: (Compat i, Coord t)+            => Name         -- ^ index of the parts to extract+            -> NArray i t   -- ^ input array+            -> Name         -- ^ index to renameRaw+            -> String       -- ^ prefix for the new names+            -> [NArray i t] -- ^ list or results+renameParts p t x pre = zipWith renameExplicit [[(x,pre ++ show k)] | k<-[1::Int ..] ] (parts t p)+++applyAsMatrix :: (Coord t, Compat i) => (Matrix t -> Matrix t) -> (NArray i t -> NArray i t)+applyAsMatrix f t = flip renameRaw nms . fromMatrix r c . f . asMatrix $ t+    where [r,c] = map (flip typeOf t) nms+          nms = sort . namesR $ t

@@ -43,7 +43,7 @@             | otherwise                 = -1  gsym f t = mkNArray (dims t) (coords $ sum ts) where-    ns = map show [1 .. rank t]+    ns = map show [1 .. order t]     t' = cov $ renameRaw t ns     per = permutations ns     ts  = map (flip renameRaw ns . f . flip reorder t') per@@ -51,11 +51,11 @@ -- symmetrize t = gsym id t  antisymmetrize t = gsym scsig t-    where scsig x = scalar (signature (names x)) * x+    where scsig x = scalar (signature (namesR x)) * x  fact n = product [1..n] -wedge a b = antisymmetrize (a*b) * (recip . fromIntegral) (fact (rank a) * fact (rank b))+wedge a b = antisymmetrize (a*b) * (recip . fromIntegral) (fact (order a) * fact (order b))  infixl 5 /\ -- | The exterior (wedge) product of two tensors. Obtains the union of subspaces.@@ -75,7 +75,7 @@  levi n = listTensor (replicate n n) $ map signature $ sequence (replicate n [1..n]) --- | The full antisymmetric tensor of rank n (contravariant version).+-- | The full antisymmetric tensor of order n (contravariant version). leviCivita :: Int -> Tensor Double leviCivita = (map levi [0..] !!) @@ -103,12 +103,12 @@       => Tensor t       -> Tensor t       -> Tensor t-inner a b | rank a < rank b = switch (renseq a) * renseq b * k+inner a b | order a < order b = switch (renseq a) * renseq b * k           | otherwise       = renseq a * switch (renseq b) * k-    where k = recip . fromIntegral $ fact $ min (rank a) (rank b)+    where k = recip . fromIntegral $ fact $ min (order a) (order b) -renseq t = renameRaw t (map show [1..rank t])-renseq' t = renameRaw t (map ((' ':).show) [1..rank t])+renseq t = renameRaw t (map show [1..order t])+renseq' t = renameRaw t (map ((' ':).show) [1..order t])  isScalar = null . dims 

@@ -26,6 +26,7 @@  import Numeric.LinearAlgebra(toList,reshape,(<\>),(@>)) import Numeric.LinearAlgebra.Array.Internal hiding (scalar,coords)+import Numeric.LinearAlgebra.Array.Display (showBases) import Numeric.LinearAlgebra.Tensor hiding (scalar,vector) import qualified Numeric.LinearAlgebra.Array.Internal as Array import Data.List@@ -228,7 +229,7 @@  mat rowidx t = reshape c $ Array.coords t'     where c = iDim $ last (dims t')-          t' = reorder (rowidx: (names t\\[rowidx])) t+          t' = reorder (rowidx: (namesR t\\[rowidx])) t  -- on the right pmat k b = mat "k" $ g!"ijk" * tb!"j"

@@ -27,24 +27,33 @@ ) where  import Numeric.LinearAlgebra.Array.Internal-import Numeric.LinearAlgebra hiding (rank)+import Numeric.LinearAlgebra import Numeric.LinearAlgebra.Array+import Data.List(intersperse)  type Tensor t = NArray Variant t -data Variant = Co | Contra deriving (Eq)+data Variant = Contra | Co deriving (Eq,Show)  instance Compat Variant where     compat d1 d2 = iDim d1 == iDim d2 && iType d1 /= iType d2+    opos (Idx x n s) = Idx (flipV x) n s  instance Show (Idx Variant) where-    show (Idx n s Co)     = show n ++ "_" ++ s-    show (Idx n s Contra) = show n ++ "^" ++ s+    show (Idx Co n s)     = s ++ "_" ++ show n+    show (Idx Contra n s) = s ++ "^" ++ show n  instance (Coord t) => Show (Tensor t) where-    show t | null (dims t) = show (coords t @>0)-           | otherwise = "listTensor " ++ show (dims t) ++ " "++ show (toList (coords t))+    show t | null (dims t) = "scalar "++ show (coords t @>0)+           | order t == 1 = ixn ++ show n ++" " ++ (show . toList . coords $ t)+           | otherwise = ixn ++ show n ++ " [" ++ ps ++ "]"+      where n = head (namesR t)+            ps = concat $ intersperse ", " $ map show (parts t n)+            ixn = idxn (typeOf n t)+            idxn Co     = "subindex "+            idxn Contra = "superindex " + flipV Co = Contra flipV Contra = Co @@ -56,7 +65,7 @@            -> [t]   -- ^ coordinates            -> Tensor t listTensor ds cs = mkNArray dms (product ds' |> (cs ++ repeat 0))-    where dms = zipWith3 Idx ds' (map show [1::Int ..]) (map f ds)+    where dms = zipWith3 Idx (map f ds) ds' (map show [1::Int ..])           ds' = map abs ds           f n | n>0       = Contra               | otherwise = Co@@ -89,14 +98,14 @@  -------------------------------------------------------------- --- | Create a contravariant rank-1 tensor from a list of coordinates.+-- | Create a contravariant 1st order tensor from a list of coordinates. vector :: [Double] -> Tensor Double vector   = fromVector Contra . fromList --- | Create a covariant rank-1 tensor from a list of coordinates.+-- | Create a covariant 1st order tensor from a list of coordinates. covector :: [Double] -> Tensor Double covector = fromVector Co . fromList --- | Create a 1-contravariant, 1-covariant rank-2 tensor from list of lists of coordinates.+-- | Create a 1-contravariant, 1-covariant 2nd order from list of lists of coordinates. transf :: [[Double]] -> Tensor Double transf   = fromMatrix Contra Co . fromLists