text:
--
--
-- <BASE HREF="http://www.numeric-quest.com/haskell/Orthogonals.html">

-- <title>
--         Indexless linear algebra algorithms
-- </title>
--
--
--

--
--

--                         ***
--

--

--         Indexless linear algebra algorithms
--

--
--

--         Orthogonalization, linear equations, eigenvalues and eigenvectors
--

--         Literate Haskell module Orthogonals.lhs
--
--

--         Jan Skibinski,
--         Numeric Quest Inc., Huntsville, Ontario, Canada
--

--         1998.09.19, last modified 1998.12.28
--
--

--

--         It has been argued that the functional paradigm offers more
--         support for scientific computing than the traditional imperative
--         programming, such as greater similarity of functional implementation
--         to mathematical specification of a problem. However, efficiency
--         of scientific algorithms implemented in Haskell is very low compared
--         to efficiencies of C or Fortran implementations - notwithstanding
--         the exceptional descriptive power of Haskell.
--

--         It has been also argued that tradition and inertia are partially
--         responsible for this sore state and that many functional algorithms
--         are direct translations of their imperative counterparts.
--

--         Arrays - with their indexing schemes and destructive updating
--         are basic tools of imperative programming. But pure functional
--         languages, which prohibit variable reassignments, cannot compete
--         with imperative languages by using the same tools and following
--         similar reasoning and patterns - unless the functional arrays
--         themselves are designed with performance in mind. This is
--         a case with Clean, where efficiency of one kind of their arrays
--         -- strict unboxed array, approaches efficiency of C.
--

--         But this has not been done for Haskell arrays yet. They are
--         lazy, boxed and use auxiliary association lists (index, value)
--         for initialization -- the latter being mostly responsible for
--         low efficiency of those algorithms that create many interim
--         arrays.
--

--         It appears, that -- as long as indexing scheme is not used
--         for lookups and updates -- Haskell lists are more efficient
--         than arrays -- at least at the currents state of Haskell.
--

--         With this in mind, we are attempting to demonstrate here
--         that the indexing traps can be successfully avoided.
--         This module implements afresh several typical problems from linear
--         algebra. Standard Haskell lists are employed instead of arrays
--         and not a single algorithm ever uses indices for lookups
--         or updates.
--

--         We do not claim high efficiency of these algorithms; consider
--         them exploratory. However, we do claim that the clarity of
--         these algorithms is significantly better than of those functionally
--         similar algorithms that employ indexing schemes.
--

--         Two major algorithms have been invented and implemented in Haskell:
--         one for solving systems of linear equations and one for finding
--         eigenvalues and eigenvectors of almost any type of a square matrix.
--         This includes symmetric, hermitian, general complex or nonsymmetric
--         matrices with real eigenvalues.
--

--         Amazingly, both methods are based on the same factorization, akin
--         to QR method, but not exactly the same as the standard QR one.
--         A simple trick allows to extend this method to nonsymmetric real
--         matrices with complex eigenvalues and thus one method applies to
--         all types of matrices.
--         It follows that the eigenvalue/eigenvector problem can be consistently
--         treated all across the board. In addition, no administrative
--         (housekeeping) boring trivia is required here and that helps to
--         clearly explain the mechanisms employed.

--
--

--

--

--
--         Contents
--
--

--

--

--         Notation
--

--         Scalar products and vector normalization
--

--

--             bra_ket, scalar product
--

--             sum_product, a cousin of bra_ket
--

--             norm, vector norm
--

--             normalized, vector normalized to one
--

--

--         Transposition and adjoining of matrices
--

--

--             transposed, transposed matrix
--

--             adjoint, transposed and conjugated matrix
--

--

--         Products involving matrices
--

--

--             matrix_matrix, product of two matrices as list of rows
--

--             matrix_matrix', product of two matrices as list of columns
--

--             triangle_matrix', upper triangular matrix times square matrix
--

--             matrix_ket, matrix times ket vector
--

--             bra_matrix, bra vector times matrix
--

--             bra_matrix_ket, matrix multiplied on both sides by vectors
--

--             scalar_matrix, scalar times matrix
--

--

--         Orthogonalization process
--

--

--             orthogonals, set of orthogonal vectors
--

--             gram_schmidt, vector perpendicular to a hyperplane
--

--

--         Solutions of linear equations by orthogonalization
--

--

--             one_ket_triangle, triangularization of one vector equation
--

--             one_ket_solution, solution for one unknown vector
--

--             many_kets_triangle, triangularization of several vector equations
--

--             many_kets_solution, solution for several unknown vectors
--

--

--         Matrix inversion
--

--

--             inverse, inverse of a matrix
--

--

--         QR factorization of matrices provided by "many_kets_triangle"
--

--

--             factors_QR, QR alike factorization of matrices
--

--             determinant, computation of the determinant based on the QR factorization
--

--

--         Similarity transformations and eigenvalues
--

--

--             similar_to, matrix obtained by similarity transformation
--

--             iterated_eigenvalues, list of approximations of eigenvalues
--

--             eigenvalues, final approximation of eigenvalues
--

--

--         Preconditioning of real nonsymmetric matrices
--

--

--             add_to_diagonal, simple preconditioning method
--

--

--         Examples of iterated eigenvalues
--

--

--             Symmetric real matrix
--

--             Hermitian complex matrix
--

--             General complex matrix
--

--             Nonsymmetric real matrix with real eigenvalues
--

--             Nonsymmetric real matrix with complex eigenvalues
--

--

--         Eigenvectors for distinct eigenvalues
--

--

--                 eigenkets, eigenvectors for distinct eigenvalues
--

--

--         Eigenvectors for degenerated eigenvalues
--

--

--                 eigenket', eigenvector based on a trial vector
--

--

--         Auxiliary functions
--

--

--         unit_matrix, a unit matrix with 1's on a diagonal
--

--         unit_vector, a vector with one non-zero component
--

--         diagonals, vector made of a matrix diagonal
--

--

--

--

--

--
--         Notation
--
--

--         What follows is written in Dirac's notation, as used
--         in Quantum Mechanics. Matrices are represented by capital
--         letters, while vectors come in two varieties:
--

--

--         Bra vector x, written < x |, is represented by one-row matrix
--

 Ket vector y, written | y >, is represented by one-column matrix
--

--

--         Bra vectors can be obtained from ket vectors by transposition
--         and conjugation of their components. Conjugation is only
--         important for complex vectors.
--

--         Scalar product of two vectors | x > and | y > is written
--         as
--

--         < x | y >
--

--         which looks like a bracket and is sometimes called a "bra_ket".
--         This justifies "bra" and "ket" names introduced by Dirac. There
--         is a good reason for conjugating the components of "bra-vector":
--         the scalar product of
--

--         < x | x >
--

--         should be a square of the norm of the vector "x", and that
--         means that it should be represented by a real number, or complex
--         number but with its imaginary part equal to zero.
--

--

--

--

 module Orthogonals where
 import Data.Complex
 import Data.Ratio
 import qualified Data.List as List

--

--
--         Scalar product and vector normalization
--
--

--         The scalar product "bra_ket" is a basis of many algorithms
--         presented here.

--

 bra_ket :: (Scalar a, Num a) => [a] -> [a] -> a
 bra_ket u v =
       --
       -- Scalar product of two vectors u and v,
       -- or  in Dirac's notation.
       -- This is equally valid for both: real and complex vectors.
       --
       sum_product u (map coupled v)

--

--         Notice the call to function "coupled" in the above implementation
--         of scalar product. This function conjugates its argument
--         if it is complex, otherwise does not change it. It is defined
--         in the class Scalar - specifically designed for this purpose
--         mainly.
--

--         This class also defines a norm of a vector that might be used
--         by some algorithms. So far we have been able to avoid this.
--

 class Eq a => Scalar a where
     coupled    :: a->a
     norm       :: [a] -> a
     almostZero :: a -> Bool
     scaled     :: [a] -> [a]

 instance Scalar Double where
     coupled x    = x
     norm u       = sqrt (bra_ket u u)
     almostZero x = (abs x) < 1.0e-8
     scaled       = scaled'

 instance Scalar Float where
     coupled x    = x
     norm u       = sqrt (bra_ket u u)
     almostZero x = (abs x) < 1.0e-8
     scaled       = scaled'

 instance (Integral a) => Scalar (Ratio a) where
     coupled x    = x
     -- norm u    = fromDouble ((sqrt (bra_ket u u))::Double)
     -- Intended hack to silently convert to and from Double.
     -- But I do not know how to declare it properly.
     --
     -- Our type Fraction, when used instead of Ratio a, has its own
     -- definition of sqrt. No hack would be needed here.
     almostZero x = abs x < 1e-8
     scaled       = scaled'

 instance (RealFloat a) => Scalar (Complex a) where
     coupled (x:+y) = x:+(-y)
     norm u         = sqrt (realPart (bra_ket u u)) :+ 0
     almostZero z   = (realPart (abs z)) < 1.0e-8
     scaled u       = [(x/m):+(y/m) | x:+y <- u]
        where m = maximum [max (abs x) (abs y) | x:+y <- u]

 norm1 :: (Num a) => [a] -> a
 norm1 = sum . map abs

 norminf :: (Num a, Ord a) => [a] -> a
 norminf = maximum . map abs

 matnorm1 :: (Num a, Ord a) => [[a]] -> a
 matnorm1 = matnorminf . transposed

 matnorminf :: (Num a, Ord a) => [[a]] -> a
 matnorminf = maximum . map norm1

--

--         But we also need a slightly different definition of
--         scalar product that will appear in multiplication of matrices
--         by vectors (or vice versa): a straightforward accumulated product
--         of two lists, where no complex conjugation takes place.
--         We will call it a 'sum_product".
--

 sum_product :: Num a => [a] -> [a] -> a
 sum_product u v =
       --
       -- Similar to scalar product but without
       -- conjugations of | u > components
       -- Used in matrix-vector or vector-matrix products
       --
       sum (zipWith (*) u v)

--

--         Some algorithms might need vectors normalized to one, although
--         we'll try to avoid the normalizations due to its high cost
--         or its inapplicability to rational numbers. Instead, we wiil
--         scale vectors by their maximal components.
--

 normalized :: (Scalar a, Fractional a) => [a] -> [a]
 normalized u =
       map (/norm u) u

 scaled' :: (Fractional t, Ord t) => [t] -> [t]
 scaled' u =
       map (/norminf u) u

--

--

--

--
--         Transposition and adjoining of matrices
--
--

--         Matrices are represented here by lists of lists.
--         Function "transposed" converts from row-wise to column-wise
--         representation, or vice versa.
--

--         When transposition is combined with complex conjugation
--         the resulting matrix is called "adjoint".
--

--         A square matrix is called symmetric if it is equal to its transpose
--

--         A = AT
--

--         It is called Hermitian, or self-adjoint, if it equals to
--         its adjoint
--

--         A = A+

 transposed :: [[a]] -> [[a]]
 transposed a
     | null (head a) = []
     | otherwise = ([head mi| mi <- a])
                   :transposed ([tail mi| mi <- a])

 adjoint :: Scalar a => [[a]] -> [[a]]
 adjoint a
     | null (head a) = []
     | otherwise = ([coupled (head mi)| mi <- a])
                   :adjoint ([tail mi| mi <- a])

--

--

--

--

--
--         Linear combination and sum of two matrices
--
--

--         One can form a linear combination of two matrices, such
--         as
--

--         C = alpha A + beta B
--         where
--             alpha and beta are scalars
--

--         The most generic form of any combination, not necessary
--         linear, of components of two matrices is given by "matrix_zipWith"
--         function below, which accepts a function "f" describing such
--         combination. For the linear combination with two scalars
--         the function "f" could be defined as:
--

--         f alpha beta a b = alpha*a + beta*b
--

--         For a straightforward addition of two matrices this auxiliary
--         function is simply "(+)".
--

 matrix_zipWith :: (a -> b -> c) -> [[a]] -> [[b]] -> [[c]]
 matrix_zipWith f a b =
     --
     -- Matrix made of a combination
     -- of matrices a and b - as specified by f
     --
     [zipWith f ak bk | (ak,bk) <- zip a b]

 add_matrices :: (Num a) => t -> t1 -> [[a]] -> [[a]] -> [[a]]
 add_matrices _ _ = matrix_zipWith (+)

--

--

--

--

--
--         Products involving matrices
--
--

--         Variety of products involving matrices can be defined.
--         Our Haskell implementation is based on lists of lists
--         and therefore is open to interpretation: sublists
--         can either represent the rows or the columns of a matrix.
--

--         The following definitions are somehow arbitrary, since
--         one can choose alternative interpretations of lists
--         representing matrices.
--

--
--         C = A B
--
--

--         Inner product of two matrices A B can be expressed quite simply,
--         providing that matrix A is represented by a list of rows
--         and B - by a list of columns. Function "matrix_matrix"
--         answers list of rows, while "matrix_matrix'" - list
--         of columns.
--

--         Major algorithms of this module make use of "triangle_matrix'",
--         which calculates a product of upper triangular matrix
--         with square matrix and returns a rectangular list of columns.

--

 matrix_matrix :: Num a => [[a]] -> [[a]] -> [[a]]
 matrix_matrix a b
 --
 -- A matrix being an inner product
 -- of matrices A and B, where
 --     A is represented by a list of rows a
 --     B is represented by a list of columns b
 --     result is represented by list of rows
 -- Require: length of a is equal of length of b
 -- Require: all sublists are of equal length

       | null a = []
       | otherwise = ([sum_product (head a) bi | bi <- b])
                  : matrix_matrix (tail a) b

 matrix_matrix' :: (Num a) => [[a]] -> [[a]] -> [[a]]
 matrix_matrix' a b
       --
       -- Similar to "matrix_matrix"
       -- but the result is represented by
       -- a list of columns
       --
       | null b = []
       | otherwise = ([sum_product ai (head b) | ai <- a])
                    : matrix_matrix' a (tail b)

 triangle_matrix' :: Num a => [[a]] -> [[a]] -> [[a]]
 triangle_matrix' r q =
       --
       -- List of columns of of a product of
       -- upper triangular matrix R and square
       -- matrix Q
       -- where
       --     r is a list of rows of R
       --     q is a list of columns of A
       --
       [f r qk | qk <- q]
       where
           f t u
               | null t = []
               | otherwise = (sum_product (head t) u)
                             : (f (tail t) (tail u))

--

--
--         | u > = A | v >
--
--

--         Product of a matrix and a ket-vector is another
--         ket vector. The following implementation assumes
--         that list "a" represents rows of matrix A.
--

 matrix_ket :: Num a => [[a]] -> [a] -> [a]
 matrix_ket a v = [sum_product ai v| ai <- a]

--

--
--         < u | = < v | A
--
--

--         Bra-vector multiplied by a matrix produces
--         another bra-vector. The implementation below
--         assumes that list "a" represents columns
--         of matrix A. It is also assumed that vector
--         "v" is given in its standard "ket" representation,
--         therefore the definition below uses "bra_ket"
--         instead of "sum_product".
--

 bra_matrix :: (Scalar a, Num a) => [a] -> [[a]] -> [a]
 bra_matrix v a = [bra_ket v ai | ai <- a]

--

--
--         alpha = < u | A | v >
--
--

--         This kind of product results in a scalar and is often
--         used to define elements of a new matrix, such as
--

--         B[i,j] = < ei | A | ej >
--

--         The implementation below assumes that list "a" represents
--         rows of matrix A.
--

 bra_matrix_ket :: (Scalar a, Num a) => [a] -> [[a]] -> [a] -> a
 bra_matrix_ket u a v =
     bra_ket u (matrix_ket a v)

--

--
--         B = alpha A
--
--

--         Below is a function which multiplies matrix by a scalar:
--

 scalar_matrix :: Num a => a -> [[a]] -> [[a]]
 scalar_matrix alpha a =
       [[alpha*aij| aij <- ai] | ai<-a]

--

--

--

--

--
--         Orthogonalization process
--
--

--         Gram-Schmidt orthogonalization procedure is used here
--         for calculation of sets of mutually orthogonal vectors.
--

--         Function "orthogonals" computes a set of mutually orthogonal
--         vectors - all orthogonal to a given vector. Such set plus
--         the input vector form a basis of the vector space. Another
--         words, they are the base vectors, although we cannot call them
--         unit vectors since we do not normalize them for two reasons:
--

--

--         None of the algorithms presented here needs this -- quite
--         costly -- normalization.
--

--         Some algorithms can be used either with doubles or with
--         rationals. The neat output of the latter is sometimes desirable
--         for pedagogical or accuracy reasons. But normalization requires "sqrt"
--         function, which is not defined for rational numbers. We could
--         use our module Fraction instead, where "sqrt" is defined,
--         but we'll leave it for a future revision of this module.
--

--

--         Function "gram_schmidt" computes one vector - orthogonal
--         to an incomplete set of orthogonal vectors, which form a hyperplane
--         in the vector space. Another words, "gram_schmidt" vector is
--         perpendicular to such a hyperlane.

--

 orthogonals :: (Scalar a, Fractional a) => [a] -> [[a]]
 orthogonals x =
       --
       -- List of (n-1) linearly independent vectors,
       -- (mutually orthogonal) and orthogonal to the
       -- vector x, but not normalized,
       -- where
       --     n is a length of x.
       --
       orth [x] size (next (-1))
       where
           orth a n m
               | n == 1        = drop 1 (reverse a)
               | otherwise     = orth ((gram_schmidt a u ):a) (n-1) (next m)
               where
                   u = unit_vector m size
           size = length x
           next i = if (i+1) == k then (i+2) else (i+1)
           k = length (takeWhile (== 0) x) -- first non-zero component of x

 gram_schmidt :: (Scalar a, Fractional a) => [[a]] -> [a] -> [a]
 gram_schmidt a u =
       --
       -- Projection of vector | u > on some direction
       -- orthogonal to the hyperplane spanned by the list 'a'
       -- of mutually orthogonal (linearly independent)
       -- vectors.
       --
       gram_schmidt' a u u
       where
           gram_schmidt' [] _ w = w
           gram_schmidt' (b:bs) v w
               | all (== 0) b = gram_schmidt' bs v w
               | otherwise    = gram_schmidt' bs v w'
               where
                   w' = vectorCombination w (-(bra_ket b v)/(bra_ket b b)) b
           vectorCombination x c y
               | null x = []
               | null y = []
               | otherwise = (head x + c * (head y))
                             : (vectorCombination (tail x) c (tail y))

--

--

--

--

--
--         Solutions of linear equations by orthogonalization
--
--

--         A matrix equation for unknown vector | x >
--

--         A | x > = | b >
--

--         can be rewritten as
--

--         x1 | 1 > + x2 | 2 > + x3 | 3 > + ... + xn | n > = | b >     (7.1)
--         where
--                 | 1 >, | 2 >... represent columns of the matrix A
--

--         For any n-dimensional vector, such as "1", there exist
--         n-1 linearly independent vectors "ck" that are orthogonal to "1";
--         that is, each satisfies the relation:
--

--         < ck | 1 > = 0, for k = 1...m, where m = n - 1
--

--         If we could find all such vectors, then we could multiply
--         the equation (7.1) by each of them, and end up with m = n-1
--         following equations
--

--         < c1 | 2 > x2 + < c1 | 3 > x3 + ... < c1 | n > xn = < c1 | b >
--         < c2 | 2 > x2 + < c2 | 3 > x3 + ... < c2 | n > xn = < c2 | b >
--         .......
--         < cm | 2 > x2 + < cm | 3 > x3 + ... < cm | n > xn = < cm | b >
--

--         But the above is nothing more than a new matrix equation
--

--         A' | x' > = | b' >
--         or

--         x2 | 2'> + x3 | 3'> .... + xn | n'> = | b'>
--         where
--             primed vectors | 2' >, etc. are the columns of the new
--             matrix A'.
--

--         with the problem dimension reduced by one.

--

--         Taking as an example a four-dimensional problem and writing
--         down the successive transformations of the original equation
--         we will end up with the following triangular pattern made of
--         four vector equations:

--

--         x1 | 1 > + x2 | 2 > + x3 | 3 >  + x4 | 4 >   = | b >
--                    x2 | 2'> + x3 | 3'>  + x4 | 4'>   = | b'>
--                               x3 | 3''> + x4 | 4''>  = | b''>
--                                           x4 | 4'''> = | b'''>
--

--         But if we premultiply each vector equation by a non-zero vector
--         of our choice, say < 1 | , < 2' |, < 3'' |, and < 4''' | - chosen
--         correspondingly for equations 1, 2, 3 and 4, then the above
--         system of vector equations will be converted to much simpler
--         system of scalar equations. The result is
--         shown below in matrix representation:

--

--         | p11  p12   p13   p14 | | x1 | = | q1 |
--         | 0    p22   p23   p24 | | x2 | = | q2 |
--         | 0    0     p33   p34 | | x3 | = | q3 |
--         | 0    0     0     p44 | | x4 | = | q4 |
--

--         In effect, we have triangularized our original matrix A.
--         Below is a function that does that for any problem size:
--

 one_ket_triangle :: (Scalar a, Fractional a) => [[a]] -> [a] -> [([a],a)]
 one_ket_triangle a b
     --
     -- List of pairs: (p, q) representing
     -- rows of triangular matrix P and of vector | q >
     -- in the equation P | x > = | q >, which
     -- has been obtained by linear transformation
     -- of the original equation A | x > = | b >
     --
     | null a = []
     | otherwise = (p,q):(one_ket_triangle a' b')
     where
         p    = [bra_ket u ak | ak <- a]
         q    = bra_ket u b
         a'   = [[bra_ket ck ai | ck <- orth] | ai <- v]
         b'   = [ bra_ket ck b  | ck <- orth]
         orth = orthogonals u
         u    = head a
         v    = tail a

--

--         The triangular system of equations can be easily solved by
--         successive substitutions - starting with the last equation.

--

 one_ket_solution :: (Scalar a, Fractional a) => [[a]] -> [a] -> [a]
 one_ket_solution a b =
     --
     -- List representing vector |x>, which is
     -- a solution of the matrix equation
     --     A |x> = |b>
     -- where
     --     a is a list of columns of matrix A
     --     b is a list representing vector |b>
     --
     solve' (unzip (reverse (one_ket_triangle a b))) []
     where
         solve' (d, c) xs
             | null d  = xs
             | otherwise = solve' ((tail d), (tail c)) (x:xs)
             where
                 x = (head c - (sum_product (tail u) xs))/(head u)
                 u = head d

--

--         The triangularization procedure can be easily extended
--         to a list of several ket-vectors | b > on the right hand
--         side of the original equation A | x > = | b > -- instead
--         of just one:
--

 many_kets_triangle :: (Scalar a, Fractional a) => [[a]] -> [[a]] -> [([a],[a])]
 many_kets_triangle a b
     --
     -- List of pairs: (p, q) representing
     -- rows of triangular matrix P and of rectangular matrix Q
     -- in the equation P X = Q, which
     -- has been obtained by linear transformation
     -- of the original equation A X = B
     -- where
     --     a is a list of columns of matrix A
     --     b is a list of columns of matrix B
     --
     | null a = []
     | otherwise = (p,q):(many_kets_triangle a' b')
     where
         p    = [bra_ket u ak   | ak <- a]
         q    = [bra_ket u bk   | bk <- b]
         a'   = [[bra_ket ck ai | ck <- orth] | ai <- v]
         b'   = [[bra_ket ck bi | ck <- orth] | bi <- b]
         orth = orthogonals u
         u    = head a
         v    = tail a

--

--         Similarly, function 'one_ket_solution' can be generalized
--         to function 'many_kets_solution' that handles cases with
--         several ket-vectors on the right hand side.
--

 many_kets_solution :: (Scalar a, Fractional a) => [[a]] -> [[a]] -> [[a]]
 many_kets_solution a b =
     --
     -- List of columns of matrix X, which is
     -- a solution of the matrix equation
     --     A X = B
     -- where
     --     a is a list of columns of matrix A
     --     b is a list of columns of matrix B
     --
     solve' p q emptyLists
     where
         (p, q) = unzip (reverse (many_kets_triangle a b))
         emptyLists = [[] | _ <- [1..(length (head q))]]
         solve' a' b' x
             | null a'  = x
             | otherwise = solve' (tail a') (tail b')
                                 [(f vk xk):xk  | (xk, vk) <- (zip x v)]
             where
                 f vk xk = (vk - (sum_product (tail u) xk))/(head u)
                 u = head a'
                 v = head b'

--

--

--

--

--
--         Matrix inversion
--
--

--         Function 'many_kets_solution' can be used to compute
--         inverse of matrix A by specializing matrix B to a unit
--         matrix I:
--

--         A X = I
--

--         It follows that matrix X is an inverse of A; that is X = A-1.
--

 inverse :: (Scalar a, Fractional a) => [[a]] -> [[a]]
 inverse a = many_kets_solution a (unit_matrix (length a))
       --
       -- List of columns of inverse of matrix A
       -- where
       --     a is list of columns of A

--

--

--

--

--
--         QR factorization of matrices
--
--

--         The process described above and implemented by
--         'many_kets_triangle' function transforms the equation
--

--         A X = B
--

--         into another equation for the same matrix X
--

--         R X = S
--

--         where R is an upper triangular matrix. All operations
--         performed on matrices A and B during this process are linear,
--         and therefore we should be able to find a square matrix Q
--         that describes the entire process in one step. Indeed, assuming
--         that matrix A can be decomposed as a product of unknown matrix Q
--         and triangular matrix R and that Q-1 is an inverse of matrix Q
--         we can reach the last equation by following these steps:
--

--         A X       = B
--         (Q R) X   = B
--         Q-1 Q R X = Q-1 B
--         R X       = S
--

--         It follows that during this process a given matrix B
--         transforms to matrix S, as delivered by 'many_kets_triangle':
--

--         S = Q-1 B
--

--         from which the inverse of Q can be found:
--

--         Q-1 = S B-1
--

--         Having a freedom of choice of the right hand side matrix B
--         we can choose the unit matrix I in place of B, and therefore
--         simplify the definition of Q-1:
--

--         Q-1 = S,  if B is unit matrix
--

--         It follows that any non-singular matrix A can be decomposed
--         as a product of a matrix Q and a triangular matrix R

--

--         A = Q R
--

--         where matrices Q-1 and R are delivered by "many_kets_triangle"
--         as a result of triangularization process of equation:
--

--         A X = I
--

--         The function below extracts a pair of matrices Q and R
--         from the answer provided by "many_kets_triangle".
--         During this process it inverts matrix Q-1 to Q.
--         This factorization will be used by a sequence of similarity
--         transformations to be defined in the next section.

--

 factors_QR :: (Scalar a, Fractional a) => [[a]] -> ([[a]],[[a]])
 factors_QR a =
       --
       -- A pair of matrices (Q, R), such that
       -- A = Q R
       -- where
       --     R is upper triangular matrix in row representation
       --     (without redundant zeros)
       --     Q is a transformation matrix in column representation
       --     A is square matrix given as columns
       --
       (inverse (transposed q1),r)
       where
           (r, q1) = unzip (many_kets_triangle a (unit_matrix (length a)))

--

--

--

--

--
--         Computation of the determinant
--

--

--

 determinant :: (Scalar a, Fractional a) => [[a]] -> a
 determinant a =
    let (q,r) = factors_QR a
    -- matrix Q is not normed so we have to respect the norms of its rows
    in  product (map norm q) * product (map head r)

--

-- Naive division-free computation of the determinant by expanding the first column.
-- It consumes n! multiplications.

--

 determinantNaive :: (Num a) => [[a]] -> a
 determinantNaive [] = 1
 determinantNaive m  =
    sum (alternate
       (zipWith (*) (map head m)
           (map determinantNaive (removeEach (map tail m)))))

--

-- Compute the determinant with about n^4 multiplications
-- without division according to the clow decomposition algorithm
-- of Mahajan and Vinay, and Berkowitz
-- as presented by Günter Rote:
--
-- Division-Free Algorithms for the Determinant and the Pfaffian:
-- Algebraic and Combinatorial Approaches.

--

 determinantClow :: (Num a) => [[a]] -> a
 determinantClow [] = 1
 determinantClow m =
    let lm = length m
    in  parityFlip lm (last (newClow m
           (nest (lm-1) (longerClow m)
               (take lm (iterate (0:) [1])))))

--

-- Compute the weights of all clow sequences
-- where the last clow is closed and a new one is started.

--

 newClow :: (Num a) => [[a]] -> [[a]] -> [a]
 newClow a c =
    scanl (-) 0
          (sumVec (zipWith (zipWith (*)) (List.transpose a) c))

--

-- Compute the weights of all clow sequences
-- where the last (open) clow is extended by a new arc.

--

 extendClow :: (Num a) => [[a]] -> [[a]] -> [[a]]
 extendClow a c =
    map (\ai -> sumVec (zipWith scaleVec ai c)) a

--

-- Given the matrix of all weights of clows of length l
-- compute the weight matrix for all clows of length (l+1).
-- Take the result of 'newClow' as diagonal
-- and the result of 'extendClow' as lower triangle
-- of the weight matrix.

--

 longerClow :: (Num a) => [[a]] -> [[a]] -> [[a]]
 longerClow a c =
    let diagonal = newClow a c
        triangle = extendClow a c
    in  zipWith3 (\i t d -> take i t ++ [d]) [0 ..] triangle diagonal

--

-- Auxiliary functions for the clow determinant.

--

 {- | Compositional power of a function,
      i.e. apply the function n times to a value. -}
 nest :: Int -> (a -> a) -> a -> a
 nest 0 _ x = x
 nest n f x = f (nest (n-1) f x)

 {- successively select elements from xs and remove one in each result list -}
 removeEach :: [a] -> [[a]]
 removeEach xs =
    zipWith (++) (List.inits xs) (tail (List.tails xs))

 alternate :: (Num a) => [a] -> [a]
 alternate = zipWith id (cycle [id, negate])

 parityFlip :: Num a => Int -> a -> a
 parityFlip n x = if even n then x else -x

 {-| Weight a list of numbers by a scalar. -}
 scaleVec :: (Num a) => a -> [a] -> [a]
 scaleVec k = map (k*)

 {-| Add corresponding numbers of two lists. -}
 {- don't use zipWith because it clips to the shorter list -}
 addVec :: (Num a) => [a] -> [a] -> [a]
 addVec x [] = x
 addVec [] y = y
 addVec (x:xs) (y:ys) = x+y : addVec xs ys

 {-| Add some lists. -}
 sumVec :: (Num a) => [[a]] -> [a]
 sumVec = foldl addVec []

--

--

--

--

--
--         Similarity transformations and eigenvalues
--
--

--         Two n-square matrices A and B are called similar if there
--         exists a non-singular matrix S such that:
--

--         B = S-1 A S
--

--         It can be proven that:
--

--

--         Any two similar matrices have the same eigenvalues
--

--         Every n-square matrix A is similar to a triangular matrix
--         whose diagonal elements are the eigenvalues of A.
--

--

--         If matrix A can be transformed to a triangular or a diagonal
--         matrix Ak by a sequence of similarity transformations then
--         the eigenvalues of matrix A are the diagonal elements of Ak.

--

--         Let's construct the sequence of matrices similar to A
--

--         A, A1, A2, A3...
--

--         by the following iterations - each of which factorizes a matrix
--         by applying the function 'factors_QR' and then forms a product
--         of the factors taken in the reverse order:
--

--         A = Q R          = Q (R Q) Q-1    = Q A1 Q-1        =
--           = Q (Q1 R1) Q-1 = Q Q1 (R1 Q1) Q1-1 Q-1 = Q Q1 A2 Q1-1 Q-1 =
--           = Q Q1 (Q2 R2) Q1-1 Q-1 = ...

--

--         We are hoping that after some number of iterations some matrix
--         Ak would become triangular and therefore its diagonal
--         elements could serve as eigenvalues of matrix A. As long as
--         a matrix has real eigenvalues only, this method should work well.
--         This applies to symmetric and hermitian matrices. It appears
--         that general complex matrices -- hermitian or not -- can also
--         be handled this way. Even more, this method also works for some
--         nonsymmetric real matrices, which have real eigenvalues only.
--

--         The only type of matrices that cannot be treated by this algorithm
--         are real nonsymmetric matrices, whose some eigenvalues are complex.
--         There is no operation in the process that converts real elements
--         to complex ones, which could find their way into diagonal
--         positions of a triangular matrix. But a simple preconditioning
--         of a matrix -- described in the next section -- replaces
--         a real matrix by a complex one, whose eigenvalues are related
--         to the eigenvalues of the matrix being replaced. And this allows
--         us to apply the same method all across the board.
--

--         It is worth noting that a process known in literature as QR
--         factorization is not uniquely defined and different algorithms
--         are employed for this. The algorithms using QR factorization
--         apply only to symmetric or hermitian matrices, and Q matrix
--         must be either orthogonal or unitary.
--

--         But our transformation matrix Q is not orthogonal nor unitary,
--         although its first row is orthogonal to all other rows. In fact,
--         this factorization is only similar to QR factorization. We just
--         keep the same name to help identify a category of the methods
--         to which it belongs.
--

--         The same factorization is used for tackling two major problems:
--         solving the systems of linear equations and finding the eigenvalues
--         of matrices.
--

--         Below is the function 'similar_to', which makes a new matrix that is
--         similar to a given matrix by applying our similarity transformation.
--

--         Function 'iterated_eigenvalues' applies this transformation n
--         times - storing diagonals of each new matrix as approximations of
--         eigenvalues.
--

--         Function 'eigenvalues' follows the same process but reports the last
--         approximation only.
--

 similar_to :: (Scalar a, Fractional a) => [[a]] -> [[a]]
 similar_to a =
       --
       -- List of columns of matrix A1 similar to A
       -- obtained by factoring A as Q R and then
       -- forming the product A1 = R Q = (inverse Q) A Q
       -- where
       --     a is list of columns of A
       --
       triangle_matrix' r q
       where
           (q,r) = factors_QR a

 iterated_eigenvalues :: (Scalar a1, Fractional a1, Eq a, Num a) => [[a1]] -> a -> [[a1]]
 iterated_eigenvalues a n
       --
       -- List of vectors representing
       -- successive approximations of
       -- eigenvalues of matrix A
       -- where
       --     a is a list of columns of A
       --     n is a number of requested iterations
       --
       | n == 0 = []
       | otherwise = (diagonals a)
                     : iterated_eigenvalues (similar_to a) (n-1)

 eigenvalues :: (Scalar a1, Fractional a1, Eq a, Num a) => [[a1]] -> a -> [a1]
 eigenvalues a n
       --
       -- Eigenvalues of matrix A
       -- obtained by n similarity iterations
       -- where
       --     a are the columns of A
       --
       | n == 0    = diagonals a
       | otherwise = eigenvalues (similar_to a) (n-1)

--

--

--

--

--
--         Preconditioning of real nonsymmetric matrices
--
--

--         As mentioned above, our QR-like factorization method works
--         well with almost all kind of matrices, but with the exception
--         of a class of real nonsymmetric matrices that have
--         complex eigenvalues.
--

--         There is no mechanism in that method that would be able to
--         produce complex eigenvalues out of the real components of
--         this type of nonsymmetric matrices. Simple trivial replacement
--         of real components of a matrix by its complex counterparts
--         does not work because zero-valued imaginary components do
--         not contribute in any way to production of nontrivial
--         imaginary components during the factorization process.
--

--         What we need is a trick that replaces real nonsymmetric matrix
--         by a nontrivial complex matrix in such a way that the results
--         of such replacements could be undone when the series of
--         similarity transformations finally produced the expected
--         effect in a form of a triangular matrix.
--

--         The practical solution is surprisingly simple:
--         it's suffice to add any complex number, such as "i", to the
--         main diagonal of a matrix, and when triangularization is done
--         -- subtract it back from computed eigenvalues.
--         The explanation follows.
--

--         Consider the eigenproblem for real and nonsymmetric matrix A.
--

--         A | x > = a | x >
--

--         Let us now define a new complex matrix B, such that:
--

--         B = A + alpha I
--         where
--             I is a unit matrix and alpha is a complex scalar
--

--         It is obvious that matrices A and B commute; that is:
--

--         A B = B A
--

--         It can be proven that if two matrices commute then they
--         have the same eigenvectors. Therefore we can use vector
--         | x > of matrix A as an eigenvector of B:
--

--         B | x > = b | x >
--         B | x > = A | x > + alpha I | x >
--                 = a | x > + alpha | x >
--                 = (a + alpha) | x >
--

--         It follows that eigenvalues of B are related to the eigenvalues
--         of A by:
--

--         b = a + alpha
--

--         After eigenvalues of complex matrix B have been succesfully
--         computed, all what remains is to subtract "alpha" from them
--         all to obtain eigenvalues of A. And nothing has to be done
--         to eigenvectors of B - they are the same for A as well.
--         Simple and elegant!
--

--         Below is an auxiliary function that adds a scalar to the
--         diagonal of a matrix:

--

 add_to_diagonal :: Num a => a -> [[a]] -> [[a]]
 add_to_diagonal alpha a =
       --
       -- Add constant alpha to diagonal of matrix A
       --
       [f ai ni | (ai,ni) <- zip a [0..(length a -1)]]
       where
           f b k = p++[head q + alpha]++(tail q)
               where
                   (p,q) = splitAt k b

--

--

--

--

--
--         Examples of iterated eigenvalues
--
--

--         Here is an example of a symmetric real matrix with results
--         of application of function 'iterated_eigenvalues'.
--

--         | 7  -2  1 |
--         |-2  10 -2 |
--         | 1  -2  7 |

--          [[7.0,     10.0,    7.0],
--           [8.66667, 9.05752, 6.27582],
--           [10.7928, 7.11006, 6.09718],
--           [11.5513, 6.40499, 6.04367],
--           [11.7889, 6.18968, 6.02142],
--           [11.8943, 6.09506, 6.01068],
--           [11.9468, 6.04788, 6.00534],
--           [11.9733, 6.02405, 6.00267],
--           [11.9866, 6.01206, 6.00134],
--           [11.9933, 6.00604, 6.00067],
--           [11.9966, 6.00302, 6.00034],
--           [11.9983, 6.00151, 6.00017],
--           [11.9992, 6.00076, 6.00008],
--           [11.9996, 6.00038, 6.00004],
--           [11.9998, 6.00019, 6.00002],
--           [11.9999, 6.00010, 6.00001],
--           [11.9999, 6.00005, 6.00001]]

--           The true eigenvalues are:
--           12, 6, 6

--

--         Here is an example of a hermitian matrix. (Eigenvalues of hermitian
--         matrices are real.) The algorithm works well and converges fast.
--

--         | 2   0     i|
--         [ 0   1   0  |
--         [ -i  0   2  |

--         [[2.8     :+ 0.0, 1.0 :+ 0.0, 1.2     :+ 0.0],
--          [2.93979 :+ 0.0, 1.0 :+ 0.0, 1.06021 :+ 0.0],
--          [2.97972 :+ 0.0, 1.0 :+ 0.0, 1.02028 :+ 0.0],
--          [2.9932  :+ 0.0, 1.0 :+ 0.0, 1.0068  :+ 0.0],
--          [2.99773 :+ 0.0, 1.0 :+ 0.0, 1.00227 :+ 0.0],
--          [2.99924 :+ 0.0, 1.0 :+ 0.0, 1.00076 :+ 0.0],
--          [2.99975 :+ 0.0, 1.0 :+ 0.0, 1.00025 :+ 0.0],
--          [2.99992 :+ 0.0, 1.0 :+ 0.0, 1.00008 :+ 0.0],
--          [2.99997 :+ 0.0, 1.0 :+ 0.0, 1.00003 :+ 0.0],
--          [2.99999 :+ 0.0, 1.0 :+ 0.0, 1.00001 :+ 0.0],
--          [3.0     :+ 0.0, 1.0 :+ 0.0, 1.0     :+ 0.0],
--          [3.0     :+ 0.0, 1.0 :+ 0.0, 1.0     :+ 0.0],
--          [3.0     :+ 0.0, 1.0 :+ 0.0, 1.0     :+ 0.0]]

--

--         Here is another example: this is a complex matrix and it is not
--         even hermitian. Yet, the algorithm still works, although its
--         fluctuates around true values.
--

--         | 2-i   0      i |
--         | 0     1+i  0   |
--         |   i   0    2-i |

--         [[2.0     :+ (-1.33333), 1.0 :+ 1.0, 2.0     :+ (-0.666667)],
--          [1.89245 :+ (-1.57849), 1.0 :+ 1.0, 2.10755 :+ (-0.421509)],
--          [1.81892 :+ (-1.80271), 1.0 :+ 1.0, 2.18108 :+ (-0.197289)],
--          [1.84565 :+ (-1.99036), 1.0 :+ 1.0, 2.15435 :+ (-0.00964242)],
--          [1.93958 :+ (-2.07773), 1.0 :+ 1.0, 2.06042 :+ 0.0777281],
--          [2.0173  :+ (-2.06818), 1.0 :+ 1.0, 1.9827  :+ 0.0681793],
--          [2.04357 :+ (-2.02437), 1.0 :+ 1.0, 1.95643 :+ 0.0243654],
--          [2.03375 :+ (-1.99072), 1.0 :+ 1.0, 1.96625 :+ (-0.00928429)],
--          [2.01245 :+ (-1.97875), 1.0 :+ 1.0, 1.98755 :+ (-0.0212528)],
--          [1.99575 :+ (-1.98307), 1.0 :+ 1.0, 2.00425 :+ (-0.0169263)],
--          [1.98938 :+ (-1.99359), 1.0 :+ 1.0, 2.01062 :+ (-0.00640583)],
--          [1.99145 :+ (-2.00213), 1.0 :+ 1.0, 2.00855 :+ 0.00212504],
--          [1.9968  :+ (-2.00535), 1.0 :+ 1.0, 2.0032  :+ 0.00535265],
--          [2.00108 :+ (-2.00427), 1.0 :+ 1.0, 1.99892 :+ 0.0042723],
--          [2.00268 :+ (-2.00159), 1.0 :+ 1.0, 1.99732 :+ 0.00158978],
--          [2.00213 :+ (-1.99946), 1.0 :+ 1.0, 1.99787 :+ (-0.000541867)],
--          [2.00079 :+ (-1.99866), 1.0 :+ 1.0, 1.9992  :+ (-0.00133514)],
--          [1.99973 :+ (-1.99893), 1.0 :+ 1.0, 2.00027 :+ (-0.00106525)],
--          [1.99933 :+ (-1.9996) , 1.0 :+ 1.0, 2.00067 :+ (-0.000397997)],
--          [1.99947 :+ (-2.00013), 1.0 :+ 1.0, 2.00053 :+ 0.000134972]]

--          The true eigenvalues are
--          2 - 2i, 1 + i, 2
--

--         Some nonsymmetric real matrices have all real eigenvalues and
--         our algorithm still works for such cases. Here is one
--         such an example, which traditionally would have to be treated
--         by one of the Lanczos-like algorithms, specifically designed
--         for nonsymmetric real matrices. Evaluation of
--

--
--         iterated_eigenvalues [[2,1,1],[-2,1,3],[3,1,-1::Double]] 20
--
--

--         gives the following results
--

--         [[3.0,     -0.70818,-0.291815],
--          [3.06743, -3.41538, 2.34795],
--          [3.02238, -1.60013, 0.577753],
--          [3.00746, -2.25793, 1.25047],
--          [3.00248, -1.88764, 0.885154],
--          [3.00083, -2.06025, 1.05943],
--          [3.00028, -1.97098, 0.970702],
--          [3.00009, -2.0148,  1.01471],
--          [3.00003, -1.99268, 0.992648],
--          [3.00001, -2.00368, 1.00367],
--          [3.0,     -1.99817, 0.998161],
--          [3.0,     -2.00092, 1.00092],
--          [3.0,     -1.99954, 0.99954],
--          [3.0,     -2.00023, 1.00023],
--          [3.0,     -1.99989, 0.999885],
--          [3.0,     -2.00006, 1.00006],
--          [3.0,     -1.99997, 0.999971],
--          [3.0,     -2.00001, 1.00001],
--          [3.0,     -1.99999, 0.999993],
--          [3.0,     -2.0,     1.0]]

--          The true eigenvalues are:
--          3, -2, 1
--

--         Finally, here is a case of a nonsymmetric real matrix with
--         complex eigenvalues:
--

--         | 2 -3 |
--         | 1  0 |
--

--         The direct application of "iterated_eigenvalues" would
--         fail to produce expected eigenvalues:
--

--         1 + i sqrt(2) and 1 - i sqrt (2)
--

--         But if we first precondition the matrix by adding "i" to its diagonal:
--

--         | 2+i  -3|
--         | 1     i|
--

--         and then compute its iterated eigenvalues:
--

--
--         iterated_eigenvalues [[2:+1,1],[-3,0:+1]] 20
--
--

--         then the method will succeed. Here are the results:
--

--         [[1.0     :+ 1.66667, 1.0     :+   0.333333 ],
--         [0.600936 :+ 2.34977, 1.39906 :+ (-0.349766)],
--         [0.998528 :+ 2.59355, 1.00147 :+ (-0.593555)],
--         [1.06991  :+ 2.413,   0.93009 :+ (-0.412998)],
--         [1.00021  :+ 2.38554, 0.99979 :+ (-0.385543)],
--         [0.988004 :+ 2.41407, 1.012   :+ (-0.414074)],
--         [0.999963 :+ 2.41919, 1.00004 :+ (-0.419191)],
--         [1.00206  :+ 2.41423, 0.99794 :+ (-0.414227)],
--         [1.00001  :+ 2.41336, 0.99999 :+ (-0.413361)],
--         [0.999647 :+ 2.41421, 1.00035 :+ (-0.414211)],
--         [0.999999 :+ 2.41436, 1.0     :+ (-0.41436) ],
--         [1.00006  :+ 2.41421, 0.99993 :+ (-0.414214)],
--         [1.0      :+ 2.41419, 1.0     :+ (-0.414188)],
--         [0.99999  :+ 2.41421, 1.00001 :+ (-0.414213)],
--         [1.0      :+ 2.41422, 1.0     :+ (-0.414218)],
--         [1.0      :+ 2.41421, 0.99999 :+ (-0.414213)],
--         [1.0      :+ 2.41421, 1.0     :+ (-0.414212)],
--         [1.0      :+ 2.41421, 1.0     :+ (-0.414213)],
--         [1.0      :+ 2.41421, 1.0     :+ (-0.414213)],
--         [1.0      :+ 2.41421, 1.0     :+ (-0.414213)]]
--

--         After subtracting "i" from the last result, we will get
--         what is expected.

--

--

--

--
--         Eigenvectors for distinct eigenvalues
--
--

--         Assuming that eigenvalues of matrix A are already found
--         we may now attempt to find the corresponding aigenvectors
--         by solving the following homogeneous equation
--

--         (A - a I) | x > = 0
--

--         for each eigenvalue "a". The matrix
--

--         B = A - a I
--

--         is by definition singular, but in most cases it can be
--         triangularized by the familiar "factors_QR" procedure.
--

--         B | x > = Q R | x > = 0
--

--         It follows that the unknown eigenvector | x > is one of
--         the solutions of the homogeneous equation:

--

--         R | x > = 0
--

--         where R is a singular, upper triangular matrix with at least one
--         zero on its diagonal.
--

--         If | x > is a solution we seek, so is its scaled version
--         alpha | x >. Therefore we have some freedom of scaling choice.
--         Since this is a homogeneous equation, one of the components
--         of | x > can be freely chosen, while the remaining components
--         will depend on that choice.
--

--         To solve the above, we will be working from the bottom up of
--         the matrix equation, as illustrated in the example below:
--

--         | 0     1     1     3     | | x1 |
--         | 0     1     1     2     | | x2 |      /\
--         | 0     0     2     4     | | x3 | = 0  ||
--         | 0     0     0     0     | | x4 |      ||
--

--         Recall that the diagonal elements of any triangular matrix
--         are its eigenvalues.
--         Our example matrix has three distinct eigenvalues:
--         0, 1, 2. The eigenvalue 0 has degree of degeneration two.
--         Presence of degenerated eigenvalues complicates
--         the solution process. The complication arises when we have to
--         make our decision about how to solve the trivial scalar equations
--         with zero coefficients, such as
--

--         0 * x4 = 0
--

--         resulting from multiplication of the bottom row by vector | x >.
--         Here we have two choices: "x4" could be set to 0, or to any
--         nonzero number 1, say. By always choosing the "0" option
--         we might end up with the all-zero trivial vector --  which is
--         obviously not what we want. Persistent choice of the "1" option,
--         might lead to a conflict between some of the equations, such as
--         the equations one and four in our example.
--

--         So the strategy is as follows.
--

--         If there is at least one zero on the diagonal, find the topmost
--         row with zero on the diagonal and choose for it the solution "1".
--         Diagonal zeros in other rows would force the solution "0".
--         If the diagonal element is not zero than simply solve
--         an arithmetic equation that arises from the substitutions of
--         previously computed components of the eigenvector. Since certain
--         inaccuracies acumulate during QR factorization, set to zero all
--         very small elements of matrix R.
--

--         By applying this strategy to our example we'll end up with the
--         eigenvector
--

--         < x | = [1, 0, 0, 0]
--

--

--         If the degree of degeneration of an eigenvalue of A is 1 then the
--         corresponding eigenvector is unique -- subject to scaling.
--         Otherwise an eigenvector found by this method is one of many
--         possible solutions, and any linear combination of such solutions
--         is also an eigenvector. This method is not able to find more than one
--         solution for degenerated eigenvalues. An alternative method, which
--         handles degenerated cases, will be described in the next section.
--

--         The function below calculates eigenvectors corresponding to
--         distinct selected eigenvalues of any square matrix A, provided
--         that the singular matrix B = A - a I can still be factorized as Q R,
--         where R is an upper triangular matrix.

--

 eigenkets :: (Scalar a, Fractional a) => [[a]] -> [a] -> [[a]]
 eigenkets a u
       --
       -- List of eigenkets of a square matrix A
       -- where
       --     a is a list of columns of A
       --     u is a list of eigenvalues of A
       --     (This list does not need to be complete)
       --
       | null u        = []
       | not (null x') = x':(eigenkets a (tail u))
       | otherwise     = (eigenket_UT (reverse b) d []):(eigenkets a (tail u))
       where
           a'  = add_to_diagonal (-(head u)) a
           x'  = unit_ket a' 0 (length a')
           b   = snd (factors_QR a')
           d   = discriminant [head bk | bk <- b] 1
           discriminant v n
               | null v    = []
               | otherwise = x : (discriminant (tail v) m)
               where
                   (x, m)
                       | (head u) == 0     = (n, 0)
                       | otherwise         = (n, n)
           eigenket_UT c e xs
               | null c    = xs
               | otherwise = eigenket_UT (tail c) (tail e) (x:xs)
               where
                   x = solve_row (head c) (head e) xs

           solve_row (v:vs) n x
               | almostZero v = n
               | otherwise    = q/v
               where
                   q
                       | null x = 0
                       | otherwise = -(sum_product vs x)

           unit_ket b' m n
               | null b'              = []
               | all (== 0) (head b') = unit_vector m n
               | otherwise            = unit_ket (tail b') (m+1) n

--

--

--

--

--
--         Eigenvectors for degenerated eigenvalues
--
--

--         Few facts:
--

--

--         Eigenvectors of a general matrix A, which does not have any
--         special symmetry, are not generally orthogonal. However, they
--         are orthogonal, or can be made orthogonal, to another set of
--         vectors that are eigenvectors of adjoint matrix A+;
--         that is the matrix obtained by complex conjugation and transposition
--         of matrix A.
--

--         Eigenvectors corresponding to nondegenerated eigenvalues of
--         hermitian or symmetric matrix are orthogonal.
--

--         Eigenvectors corresponding to degenerated eigenvalues are - in
--         general - neither orthogonal among themselves, nor orthogonal
--         to the remaining eigenvectors corresponding to other
--         eigenvalues. But since any linear combination of such degenerated
--         eigenvectors is also an eigenvector, we can orthogonalize
--         them by Gram-Schmidt orthogonalization procedure.
--

--         Many practical applications deal solely with hermitian
--         or symmetric matrices, and for such cases the orthogonalization
--         is not only possible, but also desired for variety of reasons.
--

--         But the method presented in the previous section is not able
--         to find more than one eigenvector corresponding to a degenerated
--         eigenvalue. For example, the symmetric matrix
--

--             |  7  -2   1 |
--         A = | -2  10  -2 |
--             |  1  -2   7 |
--

--         has two distinct eigenvalues: 12 and 6 -- the latter
--         being degenerated with degree of two. Two corresponding
--         eigenvectors are:
--

--         < x1 | = [1, -2, 1] -- for 12
--         < x2 | = [1,  1, 1] -- for 6
--

--         It happens that those vectors are orthogonal, but this is
--         just an accidental result. However, we are missing a third
--         distinct eigenvector. To find it we need another method.
--         One possibility is presented below and the explanation
--         follows.
--

 eigenket' :: (Scalar a, Fractional a) => [[a]] -> a -> a -> [a] -> [a]
 eigenket' a alpha eps x' =
       --
       -- Eigenket of matrix A corresponding to eigenvalue alpha
       -- where
       --     a is a list of columns of matrix A
       --     eps is a trial inaccuracy factor
       --         artificially introduced to cope
       --         with singularities of A - alpha I.
       --         One might try eps = 0, 0.00001, 0.001, etc.
       --     x' is a trial eigenvector
       --
       scaled [xk' - dk | (xk', dk) <- zip x' d]
       where
           b = add_to_diagonal (-alpha*(1+eps)) a
           d = one_ket_solution b y
           y = matrix_ket (transposed b) x'

--

--         Let us assume a trial vector | x' >, such that
--

--         | x' > = | x > + | d >
--         where
--             | x > is an eigenvector we seek
--             | d > is an error of our estimation of | x >
--

--         We first form a matrix B, such that:
--

--         B = A - alpha I
--

--         and multiply it by the trial vector | x' >, which
--         results in a vector | y >
--

--         B | x' > = |y >
--

--         On another hand:
--

--         B | x' > = B | x > + B | d > = B | d >
--         because
--             B | x > = A | x > - alpha | x > = 0
--

--         Comparing both equations we end up with:
--

--         B | d > = | y >
--

--         that is: with the system of linear equations for unknown error | d >.
--         Finally, we subtract error | d > from our trial vector | x' >
--         to obtain the true eigenvector | x >.
--

--         But there is some problem with this approach: matrix B is
--         by definition singular, and as such, it might be difficult
--         to handle. One of the two processes might fail, and their failures
--         relate to division by zero that might happen during either the
--         QR factorization, or the solution of the triangular system of equations.
--

--         But if we do not insist that matrix B should be exactly singular,
--         but almost singular:
--

--         B = A - alpha (1 + eps) I
--

--         then this method might succeed. However, the resulting eigenvector
--         will be the approximation only, and we would have to experiment
--         a bit with different values of "eps" to extrapolate the true
--         eigenvector.
--

--         The trial vector | x' > can be chosen randomly, although some
--         choices would still lead to singularity problems. Aside from
--         this, this method is quite versatile, because:
--

--

--         Any random vector | x' > leads to the same eigenvector
--         for nondegenerated eigenvalues,
--

--         Different random vectors | x' >, chosen for degenerated
--         eigenvalues, produce -- in most cases -- distinct eigenvectors.
--         And this is what we want. If we need it, we can the always
--         orthogonalize those eigenvectors either internally (always
--         possible) or externally as well (possible only for hermitian
--         or symmetric matrices).
--

--         It might be instructive to compute the eigenvectors for
--         the examples used in demonstration of computation of eigenvalues.
--         We'll leave to the reader, since this module is already too obese.
--

--

--

--
--         Auxiliary functions
--
--

--         The functions below are used in the main algorithms of
--         this module. But they can be also used for testing. For example,
--         the easiest way to test the usage of resources is to use easily
--         definable unit matrices and unit vectors, as in:

--

--         one_ket_solution (unit_matrix n::[[Double]])
--                          (unit_vector 0 n::[Double])
--         where n = 20, etc.

 unit_matrix :: Num a => Int -> [[a]]
 unit_matrix m =
       --
       -- Unit square matrix of with dimensions m x m
       --
       [ [ if j==k then 1 else 0 | j <- [0 .. m-1] ] | k <- [0 .. m-1]]

 unit_vector :: Num a => Int -> Int -> [a]
 unit_vector i m =
       --
       -- Unit vector of length m
       -- with 1 at position i, zero otherwise
       map (\k -> if k==i then 1 else 0) [0 .. m-1]

 diagonals :: [[a]] -> [a]
 diagonals a =
       --
       -- Vector made of diagonal components
       -- of square matrix a
       --
       diagonals' a 0
       where
           diagonals' b n
               | null b = []
               | otherwise =
                   (head $ drop n $ head b) : (diagonals' (tail b) (n+1))

--

--

-- -----------------------------------------------------------------------------
-- --
-- -- Copyright:
-- --
-- --      (C) 1998 Numeric Quest Inc., All rights reserved
-- --
-- -- Email:
-- --
-- --      jans@numeric-quest.com
-- --
-- -- License:
-- --
-- --      GNU General Public License, GPL
-- --
-- -----------------------------------------------------------------------------
--

--

--

--
--

max_stars_repo_name: rzil/wLPAs

max_issues_repo_name: rzil/wLPAs

max_forks_repo_name: rzil/wLPAs


text:
-- |
-- Module      :  Statistics.Distribution.Dirichlet
-- Description :  Multivariate Dirichlet distribution
-- Copyright   :  (c) Dominik Schrempf, 2021
-- License     :  GPL-3.0-or-later
--
-- Maintainer  :  dominik.schrempf@gmail.com
-- Stability   :  unstable
-- Portability :  portable
--
-- Creation date: Tue Oct 20 10:10:39 2020.
module Statistics.Distribution.Dirichlet
  ( -- * Dirichlet distribution
    DirichletDistribution (ddGetParameters),
    dirichletDistribution,
    dirichletDensity,
    dirichletSample,

    -- * Symmetric Dirichlet distribution
    DirichletDistributionSymmetric (ddSymGetParameter),
    dirichletDistributionSymmetric,
    dirichletDensitySymmetric,
    dirichletSampleSymmetric,
  )
where

import Control.Monad.Primitive
import qualified Data.Vector.Unboxed as V
import Numeric.Log
import Numeric.SpecFunctions
import System.Random.MWC
import System.Random.MWC.Distributions

-- | The [Dirichlet distribution](https://en.wikipedia.org/wiki/Dirichlet_distribution).
data DirichletDistribution = DirichletDistribution
  { ddGetParameters :: V.Vector Double,
    _getDimension :: Int,
    _getNormConst :: Log Double
  }
  deriving (Eq, Show)

-- Check if vector is strictly positive.
isNegativeOrZero :: V.Vector Double -> Bool
isNegativeOrZero = V.any (<= 0)

-- Inverse multivariate beta function. Does not check if parameters are valid!
invBeta :: V.Vector Double -> Log Double
invBeta v = Exp $ logDenominator - logNominator
  where
    logNominator = V.sum $ V.map logGamma v
    logDenominator = logGamma (V.sum v)

-- | Create a Dirichlet distribution from the given parameter vector.
--
-- Return Left if:
--
-- - The parameter vector has less then two elements.
--
-- - One or more parameters are negative or zero.
dirichletDistribution :: V.Vector Double -> Either String DirichletDistribution
dirichletDistribution v
  | V.length v < 2 =
    Left "dirichletDistribution: Parameter vector is too short."
  | isNegativeOrZero v =
    Left "dirichletDistribution: One or more parameters are negative or zero."
  | otherwise = Right $ DirichletDistribution v (V.length v) (invBeta v)

-- Tolerance.
eps :: Double
eps = 1e-14

-- Check if vector is normalized with tolerance 'eps'.
isNormalized :: V.Vector Double -> Bool
isNormalized v
  | abs (V.sum v - 1.0) > eps = False
  | otherwise = True

-- | Density of the Dirichlet distribution evaluated at a given value vector.
--
-- Return 0 if:
--
-- - The value vector has a different length than the parameter vector.
--
-- - The value vector has elements being negative or zero.
--
-- - The value vector does not sum to 1.0 (with tolerance @eps = 1e-14@).
dirichletDensity :: DirichletDistribution -> V.Vector Double -> Log Double
dirichletDensity (DirichletDistribution as k c) xs
  | k /= V.length xs = 0
  | isNegativeOrZero xs = 0
  | not (isNormalized xs) = 0
  | otherwise = c * Exp logXsPow
  where
    logXsPow = V.sum $ V.zipWith (\a x -> log $ x ** (a - 1.0)) as xs

-- | Sample a value vector from the Dirichlet distribution.
dirichletSample :: PrimMonad m => DirichletDistribution -> Gen (PrimState m) -> m (V.Vector Double)
dirichletSample (DirichletDistribution as _ _) g = do
  ys <- V.mapM (\a -> gamma a 1.0 g) as
  let s = V.sum ys
  return $ V.map (/ s) ys

-- | See 'DirichletDistribution' but with parameter vector @replicate DIM VAL@.
data DirichletDistributionSymmetric = DirichletDistributionSymmetric
  { ddSymGetParameter :: Double,
    _symGetDimension :: Int,
    _symGetNormConst :: Log Double
  }
  deriving (Eq, Show)

-- Inverse multivariate beta function. Does not check if parameters are valid!
invBetaSym :: Int -> Double -> Log Double
invBetaSym k a = Exp $ logDenominator - logNominator
  where
    logNominator = fromIntegral k * logGamma a
    logDenominator = logGamma (fromIntegral k * a)

-- | Create a symmetric Dirichlet distribution of given dimension and parameter.
--
-- Return Left if:
--
-- - The given dimension is smaller than two.
--
-- - The parameter is negative or zero.
dirichletDistributionSymmetric :: Int -> Double -> Either String DirichletDistributionSymmetric
dirichletDistributionSymmetric k a
  | k < 2 =
    Left "dirichletDistributionSymmetric: The dimension is smaller than two."
  | a <= 0 =
    Left "dirichletDistributionSymmetric: The parameter is negative or zero."
  | otherwise = Right $ DirichletDistributionSymmetric a k (invBetaSym k a)

-- | Density of the symmetric Dirichlet distribution evaluated at a given value
-- vector.
--
-- Return 0 if:
--
-- - The value vector has a different dimension.
--
-- - The value vector has elements being negative or zero.
--
-- - The value vector does not sum to 1.0 (with tolerance @eps = 1e-14@).
dirichletDensitySymmetric :: DirichletDistributionSymmetric -> V.Vector Double -> Log Double
dirichletDensitySymmetric (DirichletDistributionSymmetric a k c) xs
  | k /= V.length xs = 0
  | isNegativeOrZero xs = 0
  | not (isNormalized xs) = 0
  | otherwise = c * Exp logXsPow
  where
    logXsPow = V.sum $ V.map (\x -> log $ x ** (a - 1.0)) xs

-- | Sample a value vector from the symmetric Dirichlet distribution.
dirichletSampleSymmetric ::
  PrimMonad m =>
  DirichletDistributionSymmetric ->
  Gen (PrimState m) ->
  m (V.Vector Double)
dirichletSampleSymmetric (DirichletDistributionSymmetric a k _) g = do
  ys <- V.replicateM k (gamma a 1.0 g)
  let s = V.sum ys
  return $ V.map (/ s) ys

ext: hs

max_stars_repo_name: dschrempf/dirichlet

max_issues_repo_name: dschrempf/dirichlet

max_forks_repo_name: dschrempf/dirichlet


text:
{-# LANGUAGE Rank2Types #-}
module CreateGalaxy
where

import Data.List (sort)
import Control.Monad.State
import qualified Data.Map as M
import qualified Data.Edison.Assoc.StandardMap as E

import Galaxy
import DataFunction
import Statistics
import Math
import Utils

gasGiantByMass :: Flt -> PlanetType
gasGiantByMass mass | mass < 20.0  = SmallGasGiant
                    | mass < 100.0 = MediumGasGiant
                    | mass < 400.0 = LargeGasGiant
                    | otherwise    = VeryLargeGasGiant

createPlanetType :: Flt -> Temperature -> Flt -> Rnd PlanetType
createPlanetType mass startemp orbitradius = do
  if mass < 0.001 then return Planetoid
   else if mass < 0.01 then return NoAtmosphere
   else if mass > 15.0 then return (gasGiantByMass mass)
   else createRockyPlanetAtmosphere mass startemp orbitradius

createRockyPlanetAtmosphere :: Flt -> Temperature -> Flt -> Rnd PlanetType
createRockyPlanetAtmosphere mass startemp orbitradius = do
  weather <- randomRM (1, 100 :: Int)
  let atm1 = if weather < 10
               then RockyPlanet WaterWeatherSystem
               else if weather < 30 then RockyPlanet MethaneWeatherSystem
               else if weather < 50 then RockyPlanet SulphurDioxide
               else if weather < 90 then RockyPlanet CarbonDioxide
               else RockyPlanet Nitrogen
  let ptemp = planetTemperature'' startemp orbitradius atm1
  return $! case atm1 of
    RockyPlanet WaterWeatherSystem   -> if ptemp < waterMaxTemperature && ptemp > waterMinTemperature then atm1 else RockyPlanet Nitrogen
    RockyPlanet MethaneWeatherSystem -> if ptemp < 120 && ptemp > 70  then atm1 else RockyPlanet Nitrogen
    _                                -> atm1

createSatellite :: Flt -> Flt -> (Planet () -> Rnd a) -> Temperature -> String -> Flt -> Rnd (Planet a)
createSatellite minmass maxmass genfunc startemp name orbitradius = do
  orbit <- createOrbit orbitradius
  mass <- randomRM (minmass, maxmass)
  atmosphere <- createPlanetType mass startemp orbitradius
  let ptemp = planetTemperature'' startemp orbitradius atmosphere
  let emptyplanet = Planet name orbit (BodyPhysics mass) atmosphere ptemp M.empty ()
  cont <- genfunc emptyplanet
  return $! Planet name orbit (BodyPhysics mass) atmosphere ptemp M.empty cont

createPlanet :: (Planet () -> Rnd a) -> Temperature -> String -> Flt -> Rnd (Planet a)
createPlanet genfunc startemp name orbitradius = do
  orbit <- createOrbit orbitradius
  gentype <- randomRM (1, 4 :: Int)
  mass <- case gentype of
            1 -> randomRM (50, 400)
            2 -> randomRM (20, 50)
            3 -> randomRM (0.0001, 0.1)
            _ -> randomRM (0.1, 20)
  atmosphere <- createPlanetType mass startemp orbitradius
  numsatellites <- if mass < 1.0 then return 0 else randomRM (0 :: Int, min 16 (floor (sqrt mass)))
  satellites <- zipWithM (createSatellite (0.00001 * mass) (0.01 * mass) genfunc startemp) (bodyNames name) (replicate numsatellites orbitradius)
  let temp = planetTemperature'' startemp orbitradius atmosphere
  cont <- genfunc (Planet name orbit (BodyPhysics mass) atmosphere temp M.empty ())
  return $! Planet name orbit (BodyPhysics mass) atmosphere temp (namedsToMap satellites) cont

bodyNames :: String -> [String]
bodyNames = namesFromBasenameNum

starprobs, starbinaryprobs :: [(Flt, SpectralType)]
starprobs       = [(0.001, SpectralTypeB), (0.004, SpectralTypeA), (0.02, SpectralTypeF), (0.06, SpectralTypeG), (0.12, SpectralTypeK), (1, SpectralTypeM)]
starbinaryprobs = [(0.001, SpectralTypeB), (0.010, SpectralTypeA), (0.10, SpectralTypeF), (0.25, SpectralTypeG), (0.40, SpectralTypeK), (1, SpectralTypeM)]

getProbTableValue :: (Ord a) => a -> [(a, b)] -> b
getProbTableValue n [(_, v)] = v
getProbTableValue n ((k, v):xs) = if n <= k then v else getProbTableValue n xs
getProbTableValue _ []       = error "Given probability >1?"

randomStarTemperature :: SpectralType -> Rnd Temperature
randomStarTemperature s = randomRM (specTempRange s)

createStar :: (Planet () -> Rnd a) -> String -> Flt -> Orbit -> Rnd (Star a)
createStar genfunc name maxplanetorbitradius orbit = do
  r <- randomRM (0, 1)
  let s = getProbTableValue r starprobs
  createStar' genfunc name maxplanetorbitradius orbit s

addPrefixToPlanetName :: String -> Planet a -> Planet a
addPrefixToPlanetName n p =
  let oname = planetname p
      nname = n ++ oname
  in p{planetname = nname, satellites = E.map (addPrefixToPlanetName n) (satellites p)}

createStar' :: (Planet () -> Rnd a) -> String -> Flt -> Orbit -> SpectralType -> Rnd (Star a)
createStar' genfunc name maxplanetorbitradius orbit spectraltype = do
  t <- randomStarTemperature spectraltype
  numplanets <- randomRM (64, 128)
  let planetnames = bodyNames name
  planetorbitradiuses <- separate `fmap` sort `fmap` replicateM numplanets (randomRM (0.0001, min (fromIntegral t * 2 / 120) maxplanetorbitradius))
  planets <- zipWithM (createPlanet genfunc t) (repeat "") planetorbitradiuses
  let planets' = filter (\p -> planetTemperature' t p > 30 && planetTemperature' t p < t `div` 4) planets
  let planets'' = zipWith addPrefixToPlanetName planetnames planets'
  return $! Star name t orbit (namedsToMap planets'')

namesFromBasenameCap :: String -> [String]
namesFromBasenameCap n = zipWith (++) (repeat (n ++ " ")) (map (:[]) ['A' .. 'Z'])

namesFromBasenameMin :: String -> [String]
namesFromBasenameMin n = zipWith (++) (repeat (n ++ " ")) (map (:[]) ['a' .. 'z'])

namesFromBasenameNum :: String -> [String]
namesFromBasenameNum n = zipWith (++) (repeat (n ++ " ")) (map show [(1 :: Int) ..])

createOrbit :: Flt -> Rnd Orbit
createOrbit oradius = return $! Orbit oradius 0 0 1 0

noOrbit :: Orbit
noOrbit = Orbit 0 0 0 1 0

createStars :: (Planet () -> Rnd a) -> String -> Int -> Rnd [Star a]
createStars genfunc basename 1 = do
  s <- createStar genfunc basename 3000 noOrbit
  return [s]
createStars genfunc basename numstars = do
  let names = namesFromBasenameCap basename
  orbitradiuses <- sort `fmap` replicateM numstars (randomRM (0.1, 12000))
  orbits <- mapM createOrbit orbitradiuses
  let dists = map (/4) (distances orbitradiuses)
  spectraltypes <- map (flip getProbTableValue starbinaryprobs) `fmap` replicateM numstars (randomRM (0, 1))
  stars <- zipWith4M (createStar' genfunc) names dists orbits (sort spectraltypes)
  return $! stars

createStarSystem :: (Planet () -> Rnd a) -> String -> Vector3 -> Rnd (StarSystem a)
createStarSystem genfunc ssname sspos = do
  -- http://www.cfa.harvard.edu/news/2006/pr200611.html
  -- "Most milky way stars are single"
  singular <- chance 2 3
  n <- if singular then return 1 else randomRM (2, 6 :: Int)
  stars <- (createStars genfunc) ssname n
  return $! StarSystem ssname sspos (namedsToMap stars)

create2DPoint :: (Flt, Flt) -> Rnd Vector3
create2DPoint (minc, maxc) = do
  x <- randomRM (minc, maxc)
  y <- randomRM (minc, maxc)
  return (x, y, 0)

create3DPoint :: (Flt, Flt) -> Rnd Vector3
create3DPoint (minc, maxc) = do
  x <- randomRM (minc, maxc)
  y <- randomRM (minc, maxc)
  z <- randomRM (minc, maxc)
  return (x, y, z)

ssSpacingCoefficient :: Float
ssSpacingCoefficient = 5

createGalaxy :: (Planet () -> Rnd a) -> String -> [String] -> Rnd (Galaxy a)
createGalaxy genfunc galname ssnames = do
  let numss = length ssnames
  let dim = sqrt (fromIntegral numss) * ssSpacingCoefficient -- TODO: when 3d galaxy, use cbrt
  points <- replicateM numss (create3DPoint (-dim, dim))
  sss <- zipWithM (createStarSystem genfunc) ssnames points
  return $! Galaxy galname (namedsToMap sss)

max_stars_repo_name: anttisalonen/starrover

max_issues_repo_name: anttisalonen/starrover

max_forks_repo_name: anttisalonen/starrover


text:
-- | Module providing definitions of basic quantum gates, including parametrized
-- rotation and controlled operations.
module Hoqus.Gates where

import Numeric.LinearAlgebra.Data

-- | Constant for pure imaginary unit
ii :: C
ii = 0:+1

-- | Pauli gate X, logical negation.
sx :: Matrix C
sx = (2><2) [0, 1, 1, 0] :: Matrix C

-- | Pauli gate Y.
-- sy = -ii ( sx <> sz )
sy :: Matrix C
sy = (2><2) [0, -ii, ii, 0] :: Matrix C

-- | Pauli gate Z.
sz :: Matrix C
sz = (2><2) [1, 0, 0, -1] :: Matrix C

-- | Hadamard gate.
h :: Matrix C
h = 1/sqrt(2) * (2><2) [1, 1, 1, -1] :: Matrix C

max_stars_repo_name: jmiszczak/hoqus

max_issues_repo_name: jmiszczak/hoqus

max_forks_repo_name: jmiszczak/hoqus


text:
{-# LANGUAGE
  DeriveDataTypeable ,
  DeriveFoldable ,
  ScopedTypeVariables,
  RankNTypes,
  TypeOperators,
  DeriveTraversable,
  GeneralizedNewtypeDeriving,
  ImpredicativeTypes,
  TypeFamilies,
  PackageImports,
  StandaloneDeriving,
  FlexibleContexts,
  DeriveFunctor,
  TypeSynonymInstances,
  FlexibleInstances,
  GeneralizedNewtypeDeriving,
  NoMonomorphismRestriction,
  TemplateHaskell #-}

-- Generate a dataset
-- Using a the bicicle model

import System.Directory (doesFileExist)
import qualified MultiLinear.Class as M
import qualified System.IO as IO
import Debug.Trace as D
import Linear.V1
import Data.Distributive
import Linear.V2
import Solver.RungeKutta
import Linear.V3
import Linear.Metric
import Control.Monad (when)
import Sensor.Razor9DOF
import qualified Linear.Vector as V
import Control.Proxy.Prelude.Base hiding(Product)

import Display.Cube
import Data.List(transpose)
import Local
import Exponential.SO3
import Exponential.Class
import Rotation.SO3
import Space.SO3
import Space.Class
import Vectorization
import Kalman
import Prelude hiding(sequence,readFile,writeFile)
import As

import Display.Main

import Data.Data
import Linear.Matrix
import Data.Maybe
import qualified Data.Map as M
import Data.ByteString.Lazy (readFile,writeFile)

import Graphics.Gnuplot.Simple

import System.Environment
import Data.Foldable(toList,Foldable)
import Data.Vector.Unboxed.Base
import qualified Data.Vector.Unboxed as U
import Control.Lens
import Control.Lens.Getter
import Control.Lens.Setter
import Control.Lens.TH
import Control.Monad (liftM,replicateM,(>=>))
import Data.Word
import Data.Traversable
import Control.Applicative
import Data.Distributive
import Data.Functor.Product
import Data.Functor.Compose
import qualified Numeric.LinearAlgebra as L

import System.IO
import Data.Time
import Data.Time.Calendar

import Control.Proxy hiding (Product)
import qualified Control.Proxy.Trans.State as State

import Space.SE3

import Position
import Attitude

s0 :: Position Double
s0 = (0 `Pair` expM (pure 0)) `Pair` 0 `Pair` 0 `Pair` (expM (V3 0 (-0.35) 0 )) `Pair` accLin `Pair` magLin `Pair` gyroLin `Pair` V1 9.81 `Pair` V1 1

init0 =(0,(s0,100 *^^ gyroCov))

main = do
    f:xs <- getArgs
    let flog = "result/" ++ f
    hLog <- openFile ("calibrate/" ++ f ) ReadMode
    fLog <- openFile flog WriteMode
    t <- getCurrentTime
    configureDisplay
    stop
    hClose fLog
    plot   flog

prettyShow x = do
     D.trace (show $ fst $ snd x) (return ())
     --putStrLn "Covariance : "
     --print (snd $ snd x)

imuPipe file x = do
        let stateFile = "state/" ++  file
        existFile <- lift$ doesFileExist stateFile
        if existFile
            then do
                lift $ putStrLn "File Exist"
                sfile <- lift$ IO.readFile stateFile
                let st0 = read sfile
                State.evalStateP (_2.~ st0 $ init0 ) imuKalman
            else
                State.evalStateP init0 imuKalman

state l = _2 ._1 . l
covariance l = _2 . _2 . l
acc  = orientation |.| accCalib |.| gravity
mag  = orientation |.| orientationOffset |.| magCalib |.| magField
gyr  = orientation |.| gyroCalib

linscaleV3 = liftA2 linscale
--imuKalman ::Proxy p => State.StateP (Double, (Attitude Double,Covariance Attitude Double )) p (Attitude Double) (Double, IMU Double) () (Double,[V3 Double] ) IO ()
imuKalman = init where
    init = do
        slast <- State.get
        (t0,mi@(IMU accR gyrR magR)) <- request (snd slast)
        let
            accUn = linscaleV3 (s0 ^. accCalib.deCompose ) accR
            acc = fmap (*9.81) $ fmap (/norm accUn) accUn
            gyr = linscaleV3 (s0 ^. gyroCalib.deCompose) gyrR
            magUn = linscaleV3 (s0 ^. magCalib.deCompose) magR
            mag = fmap (/norm magUn) magUn
            angles = anglesMeasure (IMU acc gyr mag)
        State.modify (_1 .~ t0)
        State.modify (_2._1.orientation.~(rotation angles))
        loop
    loop = do
        slast <- State.get
        (t,mi@(IMU accR gyrR magR)) <- request (snd slast)
        State.modify (predictionE (positionTransition accR gyrR) t gyroCov )
        State.modify (_2  %~ measure magneticField magneticCov magR)
        snew <- State.get
        let
            accUn = linscaleV3 (s0 ^. accCalib.deCompose) accR
            acc = accUn
        when (abs (norm acc ) - 9.81 <  2.0 ) $ State.modify (_2  %~ measure gravityField accelerometerCov accR)
        () <-respond((snd $ snew , mi),[angles (snew ^. _2._1.orientation), anglesMeasure mi,angles (snew ^. _2._1.orientationOffset )])
        loop

plot file = do
    buff<- System.IO.readFile file
    let convert = map read . lines
        input = convert buff :: [(((Attitude Double ,Covariance Attitude Double),IMU Double ),[V3 Double])]
        gyroLin'  = map ( unCompose .  (^. gyroCalib ) . fst . fst . fst)input
        accLin'  = map ( unCompose .  (^. accCalib ) . fst . fst . fst)input
        gyroRawVal =  map ( gyroRaw  . snd . fst)input
        accRawVal =  map ( accRaw  . snd . fst)input
        scaleparam f = map(toList . (\(Compose x )-> fmap _scale x).f . fst . fst . fst) input
        offsetparam f = map(toList . (\(Compose x )-> fmap _offset x).f . fst . fst . fst) input
        gyroVal  = map (toList . fmap (*0.055) ) $ zipWith (liftA2 linscale) gyroLin' gyroRawVal
        accVal  = map (toList . fmap (*0.055) ) $ zipWith (liftA2 linscale) accLin' accRawVal
        accScale = scaleparam (^. accCalib )
        gyroScale = scaleparam (^. gyroCalib )
        magScale = scaleparam (^. magCalib )
        magOffset= offsetparam (^. magCalib )
        gyroOffset= offsetparam (^. gyroCalib )
        accOffset= offsetparam (^. accCalib )
        grav= map(toList .(^. gravity ) . fst .  fst . fst) input
        magneticField = map(toList .(^. magField ) . fst . fst . fst) input
        angles = map (concat . map toList . snd ) input
    plotListsStyle [Title "Angles"] $ zipWith  (\x y ->(PlotStyle Lines . CustomStyle .(\i -> [LineTitle i]) $ x,y)) ["true-roll","true-pitch","true-yaw","raw-roll","raw-pitch","raw-yaw","offset-roll","offset-pitch","offset-yaw"] $ transpose angles
    plotLists [Title "Magnetometer Scale"] $ transpose $ magScale
    plotLists [Title "Accelerometer Scale"] $ transpose $ accScale
    plotLists [Title "Magnetometer Offset"] $ transpose $ magOffset
    plotLists [Title "Accelerometer Offset"] $ transpose $ accOffset
    plotLists [Title "Gyroscope Offset"] $ transpose $ gyroOffset
    plotLists [Title "Gyroscope Scale"] $ transpose $ gyroScale
    plotLists [Title "Gravity"] $ transpose $ grav
    plotLists [Title "Magnetic"] $ transpose $ magneticField
    plotLists [Title "Gyroscope Integrated"] (transpose $ scanl1 (zipWith(+)) gyroVal )
    plotLists [Title "Accelerometer Integrated"] (transpose $ scanl1 (zipWith(+)) accVal )

degrees x = x*180/pi

max_stars_repo_name: massudaw/mtk

max_issues_repo_name: massudaw/mtk

max_forks_repo_name: massudaw/mtk


text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE TupleSections         #-}
{-# LANGUAGE TypeFamilies          #-}

import           Codec.Picture.Bitmap
import           Codec.Picture.Types
import           Control.Monad.Random
import           Data.Serialize
import           Data.Word (Word8)
import           GHC.TypeLits
import           Grenade
import           Options.Applicative
import qualified Data.ByteString as B
import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra.Static as SA

type ZDim      = 5
type TotalDim  = 2+5
type OutDim    = 3
type Hidden    = 100

-- type BatchSize = 500 * 500

type CppnNet = Network '[ FullyConnected TotalDim Hidden , Tanh
                        , FullyConnected Hidden   Hidden , Softmax
                        , FullyConnected Hidden   Hidden , Tanh
                        , FullyConnected Hidden   OutDim , Logit
                        ]
                       '[ 'D1 TotalDim
                        , 'D1 Hidden   , 'D1 Hidden
                        , 'D1 Hidden   , 'D1 Hidden
                        , 'D1 Hidden   , 'D1 Hidden
                        , 'D1 OutDim   , 'D1 OutDim
                        ]

randomNet :: MonadRandom m => m CppnNet
randomNet = randomNetwork

buildInputs :: Int
            -> Int
            -> S ('D1 ZDim)
            -> [ S ('D1 TotalDim) ]
buildInputs h w (S1D z') = vects
  where
    z      = V.toList (SA.unwrap z')
    start  = -1
    end    = 1
    range  = end - start
    -- TODO: Presently only works for a square.
    step   = range / (fromIntegral w - 1)
    pts    = [ start, start + step .. end ]
    points = [ [x, y] ++ z | x <- pts, y <- pts ]

    vects  = map (S1D . SA.vector) points

netForward :: CppnNet
           -> Int
           -> Int
           -> IO (CppnNet, Image PixelRGB8)
netForward net width height = do

  z <- do
    seed <- getRandom
    return $ S1D (SA.randomVector seed SA.Uniform :: SA.R ZDim)

  let inputs       = buildInputs width height z
      outputs      = map (runNet net) inputs
      extractPixel = toPixel . extract

      pixels :: V.Vector (PixelBaseComponent PixelRGB8)
      pixels = V.fromList (concatMap extractPixel outputs)

      img :: Image PixelRGB8
      img = Image width height pixels

  return $ (net, img)

toPixel :: (Double, Double, Double)
        -> [Word8]
toPixel (r,g,b) = [ fromIntegral $ round (r * 255)
                  , fromIntegral $ round (g * 255)
                  , fromIntegral $ round (b * 255) ]

extract :: S ('D1 OutDim) -> (Double, Double, Double)
extract (S1D a) = let [r,g,b] = V.toList (SA.unwrap a)
                    in (r,g,b)

netLoad :: FilePath -> IO CppnNet
netLoad modelPath = do
  modelData <- B.readFile modelPath
  either fail return $ runGet (get :: Get CppnNet) modelData

data CppnOpts = CppnOpts Int Int (Maybe FilePath) (Maybe FilePath)

opts :: Parser CppnOpts
opts =
  CppnOpts <$> option auto (long "width"  <> short 'w' <> value 500)
           <*> option auto (long "height" <> short 'h' <> value 500)
           <*> optional (strOption (long "load"))
           <*> optional (strOption (long "save"))

main :: IO ()
main = do
  CppnOpts width height load save <- execParser (info (opts <**> helper) idm)

  net0 <- case load of
    Just loadFile -> netLoad loadFile
    Nothing -> randomNet

  (net, image) <- netForward net0 width height

  writeBitmap "test.png" image

  case save of
    Just saveFile -> B.writeFile saveFile $ runPut (put net)
    Nothing -> return ()

  putStrLn $ "Done!"

max_stars_repo_name: silky/fashion

max_issues_repo_name: silky/fashion

max_forks_repo_name: silky/fashion


text:
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeFamilies #-}

module Interpreter where

import Data.Array as A
import Data.Complex as DC
import qualified Data.IntMap.Strict as IM
import Data.Map (Map, fromList)
import qualified Data.Map as Map
import Expression
import Utils

type family Output d rc where
    Output Scalar R = Double
    Output Scalar C = DC.Complex Double
    Output One R = Array Int Double
    Output One C = Array Int (DC.Complex Double)

-- |
--
data ValMaps =
    ValMaps
        { vm0 :: Map String Double
        , vm1 :: Map String (Array Int Double)
        }
    deriving (Eq, Show, Ord)

subs :: [(String, Double)] -> [(String, Array Int Double)] -> ValMaps
subs vm0 vm1 = ValMaps (fromList vm0) (fromList vm1)

-- |
--
class Evaluable d rc where
    eval :: ValMaps -> Expression d rc -> Output d rc

-- |
--
instance Evaluable Scalar R where
    eval :: ValMaps -> Expression Scalar R -> Double
    eval valMap e@(Expression n mp) =
        case IM.lookup n mp of
            Just ([], Var name) ->
                case Map.lookup name $ vm0 valMap of
                    Just val -> val
                    _ -> error "no value associated with the variable"
            Just ([], Sum Real [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression Scalar R
                    subExp2 = Expression node2 mp :: Expression Scalar R
                 in eval valMap subExp1 + eval valMap subExp2
            Just ([], Mul Real [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression Scalar R
                    subExp2 = Expression node2 mp :: Expression Scalar R
                 in eval valMap subExp1 * eval valMap subExp2
            Just ([], Scale Real [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression Scalar R
                    subExp2 = Expression node2 mp :: Expression Scalar R
                 in eval valMap subExp1 * eval valMap subExp2
            Just ([], InnerProd Real [node1, node2]) ->
                case IM.lookup node1 mp of
                    Just ([], _) ->
                        let subExp1 = Expression node1 mp :: Expression Scalar R -- shape is [], so must be Scalar R
                            subExp2 = Expression node2 mp :: Expression Scalar R -- shape is [], so must be Scalar R
                         in eval valMap subExp1 * eval valMap subExp2
                    Just ([size], _) ->
                        let subExp1 = Expression node1 mp :: Expression One R -- shape is [size], so must be One R
                            subExp2 = Expression node2 mp :: Expression One R -- shape is [size], so must be One R
                            lst1 = A.elems $ eval valMap subExp1
                            lst2 = A.elems $ eval valMap subExp2
                         in sum $ zipWith (*) lst1 lst2
            _ -> error "expression structure Scalar R is wrong"

-- |
--
instance Evaluable Scalar C where
    eval :: ValMaps -> Expression Scalar C -> DC.Complex Double
    eval valMap e@(Expression n mp) =
        case IM.lookup n mp of
            Just ([], Sum Complex [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression Scalar C
                    subExp2 = Expression node2 mp :: Expression Scalar C
                 in eval valMap subExp1 + eval valMap subExp2
            Just ([], Mul Complex [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression Scalar C
                    subExp2 = Expression node2 mp :: Expression Scalar C
                 in eval valMap subExp1 * eval valMap subExp2
            Just ([], Scale Complex [node1, node2]) ->
                let subExp2 = Expression node2 mp :: Expression Scalar C
                    scale =
                        case nodeNumType . retrieveNode mp $ node1 of
                            Real ->
                                fromReal . eval valMap $
                                (Expression node1 mp :: Expression Scalar R)
                            Complex ->
                                eval
                                    valMap
                                    (Expression node1 mp :: Expression Scalar C)
                 in scale * eval valMap subExp2
            Just ([], RealImg [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression Scalar R
                    subExp2 = Expression node2 mp :: Expression Scalar R
                 in eval valMap subExp1 :+ eval valMap subExp2
            Just ([], InnerProd Complex [node1, node2]) ->
                case IM.lookup node1 mp of
                    Just ([], _) ->
                        let subExp1 = Expression node1 mp :: Expression Scalar C -- shape is [], so must be Scalar C
                            subExp2 = Expression node2 mp :: Expression Scalar C -- shape is [], so must be Scalar C
                         in eval valMap subExp1 * eval valMap subExp2

                    Just ([size], _) ->
                        let subExp1 = Expression node1 mp :: Expression One C -- shape is [size], so must be One C
                            subExp2 = Expression node2 mp :: Expression One C -- shape is [size], so must be One C
                            lst1 = A.elems $ eval valMap subExp1
                            lst2 = A.elems $ eval valMap subExp2
                         in sum $ zipWith (*) lst1 lst2
            _ -> error "expression structure Scalar C is wrong"

-- |
--
instance Evaluable One R where
    eval :: ValMaps -> Expression One R -> Array Int Double
    eval valMap e@(Expression n mp) =
        case IM.lookup n mp of
            Just ([size], Var name) ->
                case Map.lookup name $ vm1 valMap of
                    Just val -> val
                    _ -> error "no value associated with the variable"
            Just ([size], Sum Real [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression One R
                    subExp2 = Expression node2 mp :: Expression One R
                    lst1 = A.elems $ eval valMap subExp1
                    lst2 = A.elems $ eval valMap subExp2
                    lstRes = zipWith (+) lst1 lst2
                 in A.listArray (0, size - 1) lstRes
            Just ([size], Mul Real [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression One R
                    subExp2 = Expression node2 mp :: Expression One R
                    lst1 = A.elems $ eval valMap subExp1
                    lst2 = A.elems $ eval valMap subExp2
                    lstRes = zipWith (*) lst1 lst2
                 in A.listArray (0, size - 1) lstRes
            Just ([size], Scale Real [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression Scalar R
                    subExp2 = Expression node2 mp :: Expression One R
                    scale = eval valMap subExp1
                 in fmap (* scale) $ eval valMap subExp2
            _ -> error "expression structure One R is wrong"

-- |
--
instance Evaluable One C where
    eval :: ValMaps -> Expression One C -> Array Int (DC.Complex Double)
    eval valMap e@(Expression n mp) =
        case IM.lookup n mp of
            Just ([size], Sum Complex [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression One C
                    subExp2 = Expression node2 mp :: Expression One C
                    lst1 = A.elems $ eval valMap subExp1
                    lst2 = A.elems $ eval valMap subExp2
                    lstRes = zipWith (+) lst1 lst2
                 in A.listArray (0, size - 1) lstRes
            Just ([size], Mul Complex [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression One C
                    subExp2 = Expression node2 mp :: Expression One C
                    lst1 = A.elems $ eval valMap subExp1
                    lst2 = A.elems $ eval valMap subExp2
                    lstRes = zipWith (*) lst1 lst2
                 in A.listArray (0, size - 1) lstRes
            Just ([size], Scale Complex [node1, node2]) ->
                let subExp2 = Expression node2 mp :: Expression One C
                    lst = A.elems $ eval valMap subExp2
                    scale =
                        case nodeNumType . retrieveNode mp $ node1 of
                            Real ->
                                fromReal . eval valMap $
                                (Expression node1 mp :: Expression Scalar R)
                            Complex ->
                                eval
                                    valMap
                                    (Expression node1 mp :: Expression Scalar C)
                 in A.listArray (0, size - 1) $ map (* scale) lst
            Just ([size], RealImg [node1, node2]) ->
                let subExp1 = Expression node1 mp :: Expression One R
                    subExp2 = Expression node2 mp :: Expression One R
                    lst1 = A.elems $ eval valMap subExp1
                    lst2 = A.elems $ eval valMap subExp2
                    lstRes = zipWith (:+) lst1 lst2
                 in A.listArray (0, size - 1) lstRes
            _ -> error "expression structure One C is wrong"

max_stars_repo_name: dandoh/NeatExpression

max_issues_repo_name: dandoh/NeatExpression

max_forks_repo_name: dandoh/NeatExpression