Source code for irfpy.util.graminan

''' Module for `Gramian matrix <http://en.wikipedia.org/wiki/Gramian_matrix>`_


'''
import logging
_logger = logging.getLogger(__name__)
import numpy as np


class Graminan:
    def __init__(self, a_m):
        r''' Instance granian.

        :param a_m: M arrays with dimension of N.
            Real space only considered.

        Consider 3-D real space (N=3). Take 2 vector (M=2).
        Then, Graminan can be obtained via 

        >>> a1 = np.array([1., 0, -2])
        >>> a2 = np.array([0, -2., 1])
        >>> G2 = Graminan([a1, a2])
        >>> print(G2.m)
        2
        >>> print(G2.n)
        3
        >>> print(G2.matrix)
        [[ 5. -2.]
         [-2.  5.]]

        Determinant of G highly relate to the area and volume.
        For this case, only two vector povides only area.
        
        .. math::

            S = \sqrt{det(G_2)}


        For m=3 case, you can get volume as

        .. math::

            V = \sqrt{det(G_3)}

        Probably for m=4 case, you can get super volume as

        .. math::

            V_4 = \sqrt{det(G_4)}

        These functionalities are implemented as :meth:`get_volume` function.
        '''

        self.m = len(a_m)
        self.n = len(a_m[0])

        if self.m > self.n:
            _logger.error('The Graminan should have m<=n, while m=%d n=%d' % (self.m, self.n))
            raise ValueError('m=%d n=%d, should be m<=n' % (self.m, self.n))

        m = self.m
 
        self.matrix = np.zeros([m, m])
        for ix in range(m):
            for iy in range(m):
                self.matrix[ix, iy] = np.dot(a_m[ix], a_m[iy])
    

    def get_volume(self):
        ''' Return the volume (generalized).

        For m=2 case, area (of trapezoid) is returned.
        For m=3 case, volume (hexahedron) is returend.
        Probably for m=4 case, super volume is returend.
        '''
        return np.sqrt(np.linalg.det(self.matrix))


import unittest
import doctest


def doctests():
    return unittest.TestSuite((
        doctest.DocTestSuite(),
        ))
if __name__ == '__main__':
    unittest.main(defaultTest='doctests')