import numpy as nm

try:
    import sympy as sm
except ImportError:
    sm = None

from sfepy.base.testing import TestCommon
from sfepy.base.base import ordered_iteritems

def symarray(prefix, shape):
    """
    Copied from SymPy so that the tests pass for its different versions.
    """
    arr = nm.empty(shape, dtype=object)
    for index in nm.ndindex(shape):
        arr[index] = sm.Symbol('%s_%s' % (prefix, '_'.join(map(str, index))))
    return arr

def get_poly(order, dim, is_simplex=False):
    """
    Construct a polynomial of given `order` in space dimension `dim`,
    and integrate it symbolically over a rectangular or simplex domain
    for coordinates in [0, 1].
    """
    xs = symarray('x', dim)

    opd = max(1, int((order + 1) / dim))

    poly = 1.0
    oo = 0
    for x in xs:
        if (oo + opd) > order:
            opd = order - oo
            
        poly *= (x**opd + 1)
        oo += opd

    limits = [[xs[ii], 0, 1] for ii in range(dim)]
    if is_simplex:
        for ii in range(1, dim):
            for ip in range(0, ii):
                limits[ii][2] -= xs[ip]

    integral = sm.integrate(poly, *reversed(limits))

    return xs, poly, limits, integral

class Test(TestCommon):

    @staticmethod
    def from_conf( conf, options ):

        test = Test(conf=conf, options=options)
        return test

    def test_weight_consistency(self):
        """
        Test if integral of 1 (= sum of weights) gives the domain volume.
        """
        from sfepy.fem.quadratures import quadrature_tables

        ok = True
        for geometry, qps in ordered_iteritems(quadrature_tables):
            self.report('geometry:', geometry)
            for order, qp in ordered_iteritems(qps):
                _ok = nm.allclose(qp.weights.sum(), qp.volume,
                                  rtol=0.0, atol=1e-15)
                self.report('order %d: %s' % (order, _ok))

                ok = ok and _ok

        return ok

    def test_quadratures(self):
        """
        Test if the quadratures have orders they claim to have, using
        symbolic integration by sympy.
        """
        from sfepy.fem.quadratures import quadrature_tables, tp_geometries
        from sfepy.fem.integrals import Integral

        if sm is None:
            self.report('cannot import sympy, skipping')
            return True

        quad = Integral('test_integral')

        ok = True
        failed = []
        for geometry, qps in ordered_iteritems(quadrature_tables):
            self.report('geometry:', geometry)

            if geometry in tp_geometries:
                iter_qp = xrange(1, 11)

            else:
                iter_qp = sorted(qps.keys())

            for order in iter_qp:
                self.report('order:', order)

                aux = self._test_quadratures(quad, geometry, order)
                _ok, integral, val = aux

                if not _ok:
                    failed.append((geometry, order, float(integral), val))

                ok = ok and _ok

        if not ok:
            self.report('failed:')
            for aux in failed:
                self.report(aux)

        return ok

    def _test_quadratures(self, quad, geometry, order):
        dim = int(geometry[0])
        n_v = int(geometry[2])
        is_simplex = n_v == (dim + 1)

        porder = order if is_simplex else dim * order

        xs, poly, limits, integral = get_poly(porder, dim,
                                              is_simplex=is_simplex)

        self.report('  polynomial:', poly)
        self.report('  limits:', limits)
        self.report('  integral:', integral)

        def fun(coors):
            vals = nm.empty(coors.shape[0], dtype=nm.float64)

            subs = {}
            for ir, cc in enumerate(coors):
                for ic, x in enumerate(xs):
                    subs[x] = coors[ir,ic]
                vals[ir] = float(poly.subs(subs))

            return vals

        val = quad.integrate(fun, order=order, geometry=geometry)
        ok = nm.allclose(val, float(integral), rtol=0.0, atol=1e-14)

        self.report('  sym. == num.: %f == %f -> %s' %
                    (float(integral), val, ok))

        return ok, integral, val

@@ -145,16 +145,25 @@
 e import
+ assert_,
  ordered
@@ -802,16 +802,31 @@
  for
+ ii,
  x in 
+enumerate(
 xs
+)
 :%0A  
@@ -834,16 +834,17 @@
      if 
+(
 (oo + op
@@ -853,16 +853,42 @@
  %3E order
+) or (ii == (len(xs) - 1))
 :%0A      
@@ -899,16 +899,20 @@
   opd = 
+max(
 order - 
@@ -913,29 +913,21 @@
 der - oo
-%0A            
+, 0)%0A
 %0A       
@@ -967,16 +967,42 @@
 += opd%0A%0A
+    assert_(oo == order)%0A%0A
     limi

fix polynomial order in get_poly()