def computeTrilinearFormByList(self,
                                   gpsk, basisk, alphak,
                                   gpsi, basisi,
                                   gpsj, basisj):
        """
        Compute trilinear form for two lists of grid points
        @param gpsk: list of HashGridIndex
        @param basisk: SG++ basis for grid indices gpsk
        @param alphak: coefficients for kth grid
        @param gpsi: list of HashGridIndex
        @param basisi: SG++ basis for grid indices gpsi
        @param gpsj: list of HashGridIndex
        @param basisj: SG++ basis for grid indices gpsj
        @return: DataMatrix
        """
        print "# evals: %i^2 * %i = %i" % (len(gpsi), len(gpsk), len(gpsi) ** 2 * len(gpsk))
        A = DataMatrix(len(gpsi), len(gpsj))
        err = 0.
        # run over all rows
        for i, gpi in enumerate(gpsi):
            # run over all columns
            for j, gpj in enumerate(gpsj):
                # run over all gpks
                b, erri = self.computeTrilinearFormByRow(gpsk, basisk,
                                                         gpi, basisi,
                                                         gpj, basisj)
                # get the overall contribution in the current dimension
                value = alphak.dotProduct(b)
                A.set(i, j, value)

                # error statistics
                err += erri
        return A, err

 def merge(cls, containerList):
     if len(containerList) == 0:
         return None
     
     # determine the total number of entries
     size = 0
     for container in containerList:
         size += len(container.getValues())
         
     dim = container.getPoints().getNcols()
     
     # Copy data to the new DataVector's entry by entry
     allPoints = DataMatrix(size, dim)
     allValues = DataVector(size) 
     tmpVector = DataVector(dim)
     i = 0
     for container in containerList:
         points = container.getPoints()
         values = container.getValues()
         for j in xrange(len(values)):
             points.getRow(j, tmpVector)
             allPoints.setRow(i, tmpVector)
             allValues[i] = values[j]
             i += 1
     
     # return new DataContainer
     return DataContainer(points=allPoints, values=allValues)

    def __getSamples(self, W, T, n):
        if self.samples is None:
            # draw n ans
            ans = W.rvs(n)

            # transform them to the unit hypercube
            ans = DataMatrix(n, W.getDim())
            for i, sample in enumerate(ans):
                p = T.probabilisticToUnit(sample)
                ans.setRow(i, DataVector(p))

            return ans
        else:
            if self.samples.shape[0] == n:
                dataSamples = self.samples
            else:
                ixs = np.random.randint(0, len(self.samples), n)
                dataSamples = self.samples[ixs, :]

            # check if there are just samples for a subset of the random
            # variables. If so, add the missing ones
            if self.__ixs is not None:
                ans = W.rvs(n)

                # transform them to the unit hypercube
                for i, sample in enumerate(dataSamples):
                    ans[i, :] = T.probabilisticToUnit(ans[i, :])
                    ans[i, self.__ixs] = sample
            else:
                ans = dataSamples

            return DataMatrix(ans)

def checkPositivity(grid, alpha):
    # define a full grid of maxlevel of the grid
    gs = grid.getStorage()
    fullGrid = Grid.createLinearGrid(gs.dim())
    fullGrid.createGridGenerator().full(gs.getMaxLevel())
    fullHashGridStorage = fullGrid.getStorage()
    A = DataMatrix(fullHashGridStorage.size(), fullHashGridStorage.dim())
    p = DataVector(gs.dim())
    for i in xrange(fullHashGridStorage.size()):
        fullHashGridStorage.get(i).getCoords(p)
        A.setRow(i, p)

    res = evalSGFunctionMulti(grid, alpha, A)
    ymin, ymax, cnt = 0, -1e10, 0
    for i, yi in enumerate(res.array()):
        if yi < 0. and abs(yi) > 1e-13:
            cnt += 1
            ymin = min(ymin, yi)
            ymax = max(ymax, yi)
            A.getRow(i, p)
            print "  %s = %g" % (p, yi)
    if cnt > 0:
        print "warning: function is not positive"
        print "%i/%i: [%g, %g]" % (cnt, fullHashGridStorage.size(), ymin, ymax)
    return cnt == 0

def computePiecewiseConstantBF(grid, U, admissibleSet):
    # create bilinear form of the grid
    gs = grid.getStorage()
    A = DataMatrix(gs.size(), gs.size())
    createOperationLTwoDotExplicit(A, grid)
    # multiply the entries with the pdf at the center of the support
    p = DataVector(gs.dim())
    q = DataVector(gs.dim())

    B = DataMatrix(admissibleSet.getSize(), gs.size())
    b = DataVector(admissibleSet.getSize())
#     s = np.ndarray(gs.dim(), dtype='float')
    for k, gpi in enumerate(admissibleSet.values()):
        i = gs.seq(gpi)
        gpi.getCoords(p)
        for j in xrange(gs.size()):
            gs.get(j).getCoords(q)
#             for d in xrange(gs.dim()):
#                 # get level index
#                 xlow = max(p[0], q[0])
#                 xhigh = min(p[1], q[1])
#                 s[d] = U[d].cdf(xhigh) - U[d].cdf(xlow)

            y = float(A.get(i, j) * U.pdf(p))
            B.set(k, j, y)
            if i == j:
                b[k] = y
    return B, b

 def testOperationB(self):
     from pysgpp import Grid, DataVector, DataMatrix
     factory = Grid.createLinearBoundaryGrid(1)
     gen = factory.createGridGenerator()
     gen.regular(2)
     
     alpha = DataVector(factory.getStorage().size())
     p = DataMatrix(1,1)
     beta = DataVector(1)
     
     
     alpha.setAll(0.0)
     p.set(0,0,0.25)
     beta[0] = 1.0
     
     opb = factory.createOperationB()
     opb.mult(beta, p, alpha)
     
     self.failUnlessAlmostEqual(alpha[0], 0.75)
     self.failUnlessAlmostEqual(alpha[1], 0.25)
     self.failUnlessAlmostEqual(alpha[2], 0.5)
     self.failUnlessAlmostEqual(alpha[3], 1.0)
     self.failUnlessAlmostEqual(alpha[4], 0.0)
     
     alpha.setAll(0.0)
     alpha[2] = 1.0
     
     p.set(0,0, 0.25)
     
     beta[0] = 0.0
     
     opb.multTranspose(alpha, p, beta)
     self.failUnlessAlmostEqual(beta[0], 0.5)

    def writeSurplusesLevelWise(self, filename):
        # locate all knowledge types available
        dtypes = self.__learner.getKnowledgeTypes()
        names = ['level']
        for dtype in dtypes:
            names.append("surplus_%s" % KnowledgeTypes.toString(dtype))

        ts = self.__knowledge.getAvailableTimeSteps()
        for t in ts:
            # collect all the surpluses classifying them by level sum
            data = {}
            n = 0
            for dtype in dtypes:
                data[dtype] = self.computeSurplusesLevelWise(t, dtype)
                n = sum([len(values) for values in data[dtype].values()])

            A = DataMatrix(n, len(names))
            # add them to a matrix structure
            for i, dtype in enumerate(dtypes):
                k = 0
                for level, surpluses in data[dtype].items():
                    for j, surplus in enumerate(surpluses):
                        A.set(k + j, i + 1, surplus)
                        A.set(k + j, 0, level)
                    k += len(surpluses)
            writeDataARFF({'filename': "%s.t%s.surpluses.arff" % (filename, t),
                           'data': A,
                           'names': names})

    def __prepareDataContainer(self, data, name):
        """
        Prepare data for learning
        @param data: dictionary loaded from UQSetting
        @return dictionary {dtype: {t: <DataContainer>}}
        """
        ans = {}
        U = self.getParameters()\
                .activeParams()\
                .getIndependentJointDistribution()
        for dtype in self.getKnowledgeTypes():
            ans[dtype] = {}
            dim = self.grid.getStorage().dim()

            # prepare data container depending on the given knowledge type
            tmp = KnowledgeTypes.transformData(data, U, dtype)

            # load data for all time steps
            for t, values in tmp.items():
                size = len(values)
                mydata = DataMatrix(size, dim)
                sol = DataVector(size)
                for i, (sample, res) in enumerate(values.items()):
                    p = DataVector(sample.getActiveUnit())
                    mydata.setRow(i, p)
                    sol[i] = float(res)
                ans[dtype][t] = DataContainer(points=mydata, values=sol, name=name)
        return ans

    def computeMoments(self, ts=None):
        names = ['time',
                 'iteration',
                 'grid_size',
                 'mean',
                 'meanDiscretizationError',
                 'var',
                 'varDiscretizationError']
        # parameters
        ts = self.__samples.keys()
        nrows = len(ts)
        ncols = len(names)
        data = DataMatrix(nrows, ncols)
        v = DataVector(ncols)

        row = 0
        for t in ts:
            v.setAll(0.0)
            v[0] = t
            v[1] = 0
            v[2] = len(self.__samples[t].values())
            v[3], v[4] = self.mean(ts=[t])
            v[5], v[6] = self.var(ts=[t])

            # write results to matrix
            data.setRow(row, v)
            row += 1

        return {'data': data,
                'names': names}

    def testOperationTest_test(self):
        from pysgpp import Grid, DataVector, DataMatrix

        factory = Grid.createLinearBoundaryGrid(1)
        gen = factory.createGridGenerator()
        gen.regular(1)
        
        alpha = DataVector(factory.getStorage().size())        
        
        data = DataMatrix(1,1)
        data.setAll(0.25)
        classes = DataVector(1)
        classes.setAll(1.0)

        testOP = factory.createOperationTest()

        alpha[0] = 0.0
        alpha[1] = 0.0
        alpha[2] = 1.0
        
        c = testOP.test(alpha, data, classes)
        self.failUnless(c > 0.0)
        
        alpha[0] = 0.0
        alpha[1] = 0.0
        alpha[2] = -1.0
        c = testOP.test(alpha, data, classes)
        self.failUnless(c == 0.0)

 def plotGrid(self, learner, suffix):
     from mpl_toolkits.mplot3d.axes3d import Axes3D
     import matplotlib.pyplot as plt
     xs = np.linspace(0, 1, 30)
     ys = np.linspace(0, 1, 30)
     X, Y = np.meshgrid(xs, ys)
     Z = zeros(np.shape(X))
     input = DataMatrix(np.shape(Z)[0]*np.shape(Z)[1], 2)
     r = 0
     for i in xrange(np.shape(Z)[0]):
         for j in xrange(np.shape(Z)[1]):
             input.set(r, 0, X[i,j])
             input.set(r, 1, Y[i,j])
             r += 1
     result = learner.applyData(input)
     r = 0
     for i in xrange(np.shape(Z)[0]):
         for j in xrange(np.shape(Z)[1]):
             Z[i,j] = result[r]
             r += 1
     fig = plt.figure()
     ax = Axes3D(fig)
     ax.plot_wireframe(X,Y,Z) 
     #plt.draw()
     plt.savefig("grid3d_%s_%i.png" % (suffix, learner.iteration))
     fig.clf()
     plt.close(plt.gcf())

 def setUp(self):
     self.size = 9
     self.level = 4
     points = DataMatrix(self.size, 1)
     values = DataVector(self.size)
     for i in xrange(self.size):
         points.set(i, 0, i)
         values[i] = -1 if i < self.size/2 else 1
     self.dataContainer = DataContainer(points=points, values=values)
     self.policy = StratifiedFoldingPolicy(self.dataContainer, self.level)

def dehierarchizeOnNewGrid(gridResult, grid, alpha):
    # dehierarchization
    gsResult = gridResult.getStorage()
    ps = DataMatrix(gsResult.size(), gsResult.dim())
    p = DataVector(gsResult.dim())
    for i in xrange(gsResult.size()):
        gsResult.get(i).getCoords(p)
        ps.setRow(i, p)
    nodalValues = evalSGFunctionMulti(grid, alpha, ps)
    return nodalValues

 def setUp(self):
     self.size = 11
     self.level = 10
     points = DataMatrix(self.size, 1)
     values = DataVector(self.size)
     for i in xrange(self.size):
         points.set(i, 0, i)
         values[i] = i
     self.dataContainer = DataContainer(points=points, values=values)
     self.policy = SequentialFoldingPolicy(self.dataContainer, self.level)

def buildTrainingVector(data):
    from pysgpp import DataMatrix
    dim = len(data["data"])
    training = DataMatrix(len(data["data"][0]), dim)

    # i iterates over the data points, d over the dimension of one data point
    for i in xrange(len(data["data"][0])):
        for d in xrange(dim):
            training.set(i, d, data["data"][d][i])

    return training

    def general_test(self, d, l, bb, xs):

        test_desc = "dim=%d, level=%d, len(x)=%s" % (d, l, len(xs))

        print test_desc

        self.grid = Grid.createLinearGrid(d)
        self.grid_gen = self.grid.createGridGenerator()
        self.grid_gen.regular(l)

        alpha = DataVector([self.get_random_alpha() for i in xrange(self.grid.getSize())])

        bb_ = BoundingBox(d)

        for d_k in xrange(d):
            dimbb = DimensionBoundary()
            dimbb.leftBoundary = bb[d_k][0]
            dimbb.rightBoundary = bb[d_k][1]
            bb_.setBoundary(d_k, dimbb)

        # Calculate the expected value without the bounding box

        expected_normal = [self.calc_exp_value_normal(x, d, bb, alpha) for x in xs]
        #expected_transposed = [self.calc_exp_value_transposed(x, d, bb, alpha) for x in xs]

        # Now set the bounding box

        self.grid.getStorage().setBoundingBox(bb_)

        dm = DataMatrix(len(xs), d)
        for k, x in enumerate(xs):
            dv = DataVector(x)
            dm.setRow(k, dv)

        multEval = createOperationMultipleEval(self.grid, dm)

        actual_normal = DataVector(len(xs))
        #actual_transposed = DataVector(len(xs))

        multEval.mult(alpha, actual_normal)
        #multEval.mult(alpha, actual_transposed)

        actual_normal_list = []
        for k in xrange(len(xs)):
            actual_normal_list.append(actual_normal.__getitem__(k))

        #actual_transposed_list = []
        #for k in xrange(len(xs)):
        #    actual_transposed_list.append(actual_transposed.__getitem__(k))

        self.assertAlmostEqual(actual_normal_list, expected_normal)
        #self.assertAlmostEqual(actual_tranposed_list, expected_tranposed)

        del self.grid

def dehierarchize(grid, alpha):
    # dehierarchization
    gs = grid.getStorage()
    p = DataVector(gs.dim())
    nodalValues = DataVector(gs.size())
    A = DataMatrix(gs.size(), gs.dim())
    for i in xrange(gs.size()):
        gs.get(i).getCoords(p)
        A.setRow(i, p)
    opEval = createOperationMultipleEval(grid, A)
    opEval.mult(alpha, nodalValues)
    return nodalValues

def estimateDiscreteL2Error(grid, alpha, f, n=1000):
    gs = grid.getStorage()
    # create control samples
    samples = DataMatrix(np.random.rand(n, gs.dim()))

    nodalValues = evalSGFunctionMulti(grid, alpha, samples)
    fvalues = DataVector(samples.getNrows())
    for i, sample in enumerate(samples.array()):
        fvalues[i] = f(sample)

    # compute the difference
    nodalValues.sub(fvalues)
    return nodalValues.l2Norm()

    def sampleGrids(self, filename):
        ts = self.__learner.getTimeStepsOfInterest()

        names = self.__params.getNames()
        names.append('f_\\mathcal{I}(x)')

        for t in ts:
            grid, surplus = self.__knowledge.getSparseGridFunction(self._qoi, t)

            # init
            gs = grid.getStorage()
            dim = gs.dim()

            # -----------------------------------------
            # do full grid sampling of sparse grid function
            # -----------------------------------------
            data = eval_fullGrid(4, dim)
            res = evalSGFunctionMulti(grid, surplus, data)

            data.transpose()
            data.appendRow()
            data.setRow(data.getNrows() - 1, res)
            data.transpose()

            # write results
            writeDataARFF({'filename': "%s.t%f.samples.arff" % (filename, t),
                           'data': data,
                           'names': names})

            # -----------------------------------------
            # write sparse grid points to file
            # -----------------------------------------
            data = DataMatrix(gs.size(), dim)
            data.setAll(0.0)

            for i in xrange(gs.size()):
                gp = gs.get(i)
                v = np.array([gp.getCoord(j) for j in xrange(dim)])
                data.setRow(i, DataVector(v))

            # write results
            writeDataARFF({'filename': "%s.t%f.gridpoints.arff" % (filename, t),
                           'data': data,
                           'names': names})

            # -----------------------------------------
            # write alpha
            # -----------------------------------------
            writeAlphaARFF("%s.t%f.alpha.arff" % (filename, t),
                           surplus)

def computeBilinearFormQuad(grid, U):
    gs = grid.getStorage()
    basis = getBasis(grid)
    A = DataMatrix(gs.size(), gs.size())

    level = DataMatrix(gs.size(), gs.dim())
    index = DataMatrix(gs.size(), gs.dim())
    gs.getLevelIndexArraysForEval(level, index)

    s = np.ndarray(gs.dim(), dtype='float')
    # run over all rows
    for i in xrange(gs.size()):
        gpi = gs.get(i)
        # run over all columns
        for j in xrange(i, gs.size()):
            # print "%i/%i" % (i * gs.size() + j + 1, gs.size() ** 2)
            gpj = gs.get(j)
            for d in xrange(gs.dim()):
                # get level index
                lid, iid = level.get(i, d), index.get(i, d)
                ljd, ijd = level.get(j, d), index.get(j, d)

                # compute left and right boundary of the support of both
                # basis functions
                lb = max([(iid - 1) / lid, (ijd - 1) / ljd])
                ub = min([(iid + 1) / lid, (ijd + 1) / ljd])

                # same level, different index
                if lid == ljd and iid != ijd:
                    s[d] = 0.
                # the support does not overlap
                elif lid != ljd and lb >= ub:
                    s[d] = 0.
                else:
                    lid, iid = gpi.getLevel(d), int(iid)
                    ljd, ijd = gpj.getLevel(d), int(ijd)
                    # ----------------------------------------------------
                    # use scipy for integration

                    def f(x):
                        return basis.eval(lid, iid, x) * \
                            basis.eval(ljd, ijd, x) * \
                            U[d].pdf(x)

                    s[d], _ = quad(f, lb, ub, epsabs=1e-8)
                    # ----------------------------------------------------
            A.set(i, j, float(np.prod(s)))
            A.set(j, i, A.get(i, j))

    return A