# Copyright (C) 2008-today The SG++ project
# This file is part of the SG++ project. For conditions of distribution and
# use, please see the copyright notice provided with SG++ or at
# sgpp.sparsegrids.org

try:
  import pysgpp
  from pysgpp.extensions.datadriven.uq.operations.sparse_grid import getBasis
  import numpy as np
  import random
  import matplotlib.pyplot as plt
  
except ImportError as e:
    print(e.__class__.__name__ + ": " + e.msg)
    print("Skipping example...")
    exit(0)

def test_base():
    n = 100000
    b = pysgpp.SLinearModifiedBase()
    # print(sum([b.eval(0,0,x) for x in np.linspace(0, 1, n)]) / n)
    for l in range(1,3,1):
        print(l)
        for i in range(1, 2**l, 1):
            s = 1
            s = sum([b.eval(l,i,x) for x in np.linspace(0, 1, n)]) / n
            print(("Level:",l , "index:", i))
            epsilon = abs(s - b.getIntegral(l,i))
            if(epsilon > 10**-6):
                print("s: {} res: {}".format(s, b.getIntegral(l,i)))
                print("error: {}".format(epsilon))

def test_LTwoDot(grid, l):
    res = 10000
    b = grid.getBasis();
    grid.getGenerator().regular(l)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    # print(size)
    m = pysgpp.DataMatrix(size, size)
    opMatrix = pysgpp.createOperationLTwoDotExplicit(m, grid)
    # print m
    # m_comp = pysgpp.DataMatrix(size, size)
    for i in range(gridStorage.getSize()):
        for j in range(i, gridStorage.getSize()):
            gpi = gridStorage.getPoint(i)
            gpj = gridStorage.getPoint(j)
            sol = 1
            # print "--------"
            # print "i:{} j:{}".format(i, j)
            for k in range(d):
                lik = gpi.getLevel(k)
                iik = gpi.getIndex(k)
                ljk = gpj.getLevel(k)
                ijk = gpj.getIndex(k)
                # print "i l,i: {},{}   j l,i: {},{}".format(lik, iik, ljk, ijk)
                xs = np.linspace(0, 1, res)
                tmp = sum([b.eval(lik, iik, x)*b.eval(ljk, ijk, x) for x in xs]) / res
                sol *= tmp
                # print("lik:{} iik:{} ljk:{} ijk:{} k:{} tmp: {}".format(lik, iik, ljk, ijk, k,tmp))
            # print(sol)
            error = abs(m.get(i,j) - sol)
            # print error
            if(error >= 10**-4):
                print("i:{} j:{} error: {}".format(i, j, error))
                print("iik:{} lik:{} ijk:{} ljk:{} error: {}".format(iik, lik, ijk, ljk, error))
                print("is:{} should:{}".format(m.get(i,j), sol))

def test_LTwoDotImplicit(grid, l):
    grid.getGenerator().regular(l)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    # print(size)
    m = pysgpp.DataMatrix(size, size)
    opExplicit = pysgpp.createOperationLTwoDotExplicit(m, grid)
    op = pysgpp.createOperationLTwoDotProduct(grid)
    alpha = pysgpp.DataVector(size)
    resultExplicit = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)
    for i in range(size):
        alpha[i] = 1
    opExplicit.mult(alpha, resultExplicit)
    op.mult(alpha, result)
    for i in range(size):
        if result[i] != resultExplicit[i]:
            print("Error result entry {} differs".format(i))

        if abs(result[i] - resultExplicit[i]) > 1e-16:
            # print result[i] - resultExplicit[i]
            print("result:{}".format(result[i]))
            print("resultExplicit:{}".format(resultExplicit[i]))

def test_laplace(grid, lmax):
    resolution = 100000
    grid.getGenerator().regular(lmax)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    b = getBasis(grid)
    op = pysgpp.createOperationLaplace(grid)
    alpha = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)

    for point_i in range(size):
      for point_j in range(size):
        gp_i = gridStorage.getPoint(point_i)
        gp_j = gridStorage.getPoint(point_j)
        print("--------")
        for i in range(0, size):
          alpha[i] = 0
        alpha[point_i] = 1
        op.mult(alpha, result)
        xs = np.linspace(0, 1, resolution)
        approx = sum([b.evalDx(gp_i.getLevel(0), gp_i.getIndex(0), x) * b.evalDx(gp_j.getLevel(0), gp_j.getIndex(0), x) for x in xs]) / resolution
        print("i,j: {},{} result: {} approx:{}".format(point_i, point_j, result[point_j], approx))
        if(abs(result.get(point_j) - approx) > 1e-1):
          print("--------")
          print("points: {},{} ".format(point_i, point_j))
          print("approx:{}".format(approx))
          print("result:{}".format(result.get(point_j)))
          # print result
          print("--------")

def test_poly_evaldx():
    l = 3
    i = 1
    x = 0.12
    eps = 0.0001
    b = pysgpp.SPolyModifiedClenshawCurtisBase(3)
    tang = b.evalDx(l, i, x)
    sec = (b.eval(l, i, x + eps) -  b.eval(l, i, x - eps)) / (2*eps)
    print("evalDx:{}".format(tang))
    print("sekante:{}".format(sec))
    print("evals: {} {}".format( b.eval(l, i, x - eps), b.eval(l, i, x + eps) ))

def plot_evaldx():
    l = 3
    i = 1
    xs = np.linspace(0, 1, 50)
    b = pysgpp.SPolyClenshawCurtisBase(3)
    plt.plot(xs, [b.evalDx(l, i, x) for x in xs])
    plt.show()

def plot_evaldx_prod(grid, lmax, i, j):
    grid.getGenerator().regular(lmax)
    gridStorage = grid.getStorage()
    b = getBasis(grid)
    gp_i = gridStorage.getPoint(i)
    gp_j = gridStorage.getPoint(j)
    lik = gp_i.getLevel(0)
    ljk = gp_j.getLevel(0)
    iik = gp_i.getIndex(0)
    ijk = gp_j.getIndex(0)
    xs = np.linspace(0, 1, 100000)

    plt.plot(xs, [b.evalDx(lik, iik, x) * b.evalDx(ljk, ijk, x) for x in xs])
    plt.show()


def test_laplace2(grid, lmax):
  # grid.getGenerator().regular(lmax)
  gridStorage = grid.getStorage()
  size = gridStorage.getSize()
  m = pysgpp.DataMatrix(size, size)
  op = pysgpp.createOperationLaplaceExplicit(m, grid)
  print(m)



def test_firstMoment(grid, lmax):
  grid.getGenerator().regular(lmax)
  resolution = 100000
  gridStorage = grid.getStorage()
  b = grid.getBasis()
  op = pysgpp.createOperationFirstMoment(grid)
  alpha = pysgpp.DataVector(grid.getSize(), 1.0)
  bounds = pysgpp.DataMatrix(1, 2, 0.0)
  bounds.set(0, 1, 1.0)
  res = 0.0
  for i in range(grid.getSize()):
    lev = gridStorage.getPoint(i).getLevel(0)
    ind = gridStorage.getPoint(i).getIndex(0)
    temp_res = 0.0
    for c in range(resolution):
        x = float(c) / resolution
        temp_res += x * b.eval(lev, ind, x)
    res += alpha.get(i) * temp_res / resolution
  print("--FirstMoment--")
  print(res)
  print(op.doQuadrature(alpha, bounds))
  print (res - op.doQuadrature(alpha, bounds))

def test_secondtMoment(grid, lmax):
  gridStorage = grid.getStorage()
  gridStorage.clear()
  grid.getGenerator().regular(lmax)
  resolution = 1000000
  b = grid.getBasis()
  op = pysgpp.createOperationSecondMoment(grid)
  alpha = pysgpp.DataVector(grid.getSize(), 1.0)
  bounds = pysgpp.DataMatrix(1, 2, 0.0)
  bounds.set(0, 1, 1.0)
  res = 0.0
  for i in range(grid.getSize()):
    lev = gridStorage.getPoint(i).getLevel(0)
    ind = gridStorage.getPoint(i).getIndex(0)
    temp_res = 0.0
    for c in range(resolution):
        x = float(c) / resolution
        temp_res += x * x *b.eval(lev, ind, x)
    res += alpha.get(i) * temp_res / resolution

  print("--SecondMoment--")
  print(res)
  print(op.doQuadrature(alpha, bounds))
  print (res - op.doQuadrature(alpha, bounds))

# test_poly_evaldx()
# plot_evaldx()
# test_base()
d = 1
l = 3
grid = pysgpp.Grid.createModLinearGrid(d)
# plot_evaldx_prod(grid, 4, 1, 4)
# test_laplace(grid, l)
# test_laplace2(grid, l)
# test_LTwoDot(grid, l)
# test_LTwoDotImplicit(grid, l)
test_firstMoment(grid, l)
test_secondtMoment(grid, l)

a67d7c0d2a919c6275851afb38fac9f12c63b0a4
<mask token>


def test_base():
    n = 100000
    b = pysgpp.SLinearModifiedBase()
    for l in range(1, 3, 1):
        print(l)
        for i in range(1, 2 ** l, 1):
            s = 1
            s = sum([b.eval(l, i, x) for x in np.linspace(0, 1, n)]) / n
            print(('Level:', l, 'index:', i))
            epsilon = abs(s - b.getIntegral(l, i))
            if epsilon > 10 ** -6:
                print('s: {} res: {}'.format(s, b.getIntegral(l, i)))
                print('error: {}'.format(epsilon))


<mask token>


def test_LTwoDotImplicit(grid, l):
    grid.getGenerator().regular(l)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    m = pysgpp.DataMatrix(size, size)
    opExplicit = pysgpp.createOperationLTwoDotExplicit(m, grid)
    op = pysgpp.createOperationLTwoDotProduct(grid)
    alpha = pysgpp.DataVector(size)
    resultExplicit = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)
    for i in range(size):
        alpha[i] = 1
    opExplicit.mult(alpha, resultExplicit)
    op.mult(alpha, result)
    for i in range(size):
        if result[i] != resultExplicit[i]:
            print('Error result entry {} differs'.format(i))
        if abs(result[i] - resultExplicit[i]) > 1e-16:
            print('result:{}'.format(result[i]))
            print('resultExplicit:{}'.format(resultExplicit[i]))


def test_laplace(grid, lmax):
    resolution = 100000
    grid.getGenerator().regular(lmax)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    b = getBasis(grid)
    op = pysgpp.createOperationLaplace(grid)
    alpha = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)
    for point_i in range(size):
        for point_j in range(size):
            gp_i = gridStorage.getPoint(point_i)
            gp_j = gridStorage.getPoint(point_j)
            print('--------')
            for i in range(0, size):
                alpha[i] = 0
            alpha[point_i] = 1
            op.mult(alpha, result)
            xs = np.linspace(0, 1, resolution)
            approx = sum([(b.evalDx(gp_i.getLevel(0), gp_i.getIndex(0), x) *
                b.evalDx(gp_j.getLevel(0), gp_j.getIndex(0), x)) for x in xs]
                ) / resolution
            print('i,j: {},{} result: {} approx:{}'.format(point_i, point_j,
                result[point_j], approx))
            if abs(result.get(point_j) - approx) > 0.1:
                print('--------')
                print('points: {},{} '.format(point_i, point_j))
                print('approx:{}'.format(approx))
                print('result:{}'.format(result.get(point_j)))
                print('--------')


<mask token>


def test_firstMoment(grid, lmax):
    grid.getGenerator().regular(lmax)
    resolution = 100000
    gridStorage = grid.getStorage()
    b = grid.getBasis()
    op = pysgpp.createOperationFirstMoment(grid)
    alpha = pysgpp.DataVector(grid.getSize(), 1.0)
    bounds = pysgpp.DataMatrix(1, 2, 0.0)
    bounds.set(0, 1, 1.0)
    res = 0.0
    for i in range(grid.getSize()):
        lev = gridStorage.getPoint(i).getLevel(0)
        ind = gridStorage.getPoint(i).getIndex(0)
        temp_res = 0.0
        for c in range(resolution):
            x = float(c) / resolution
            temp_res += x * b.eval(lev, ind, x)
        res += alpha.get(i) * temp_res / resolution
    print('--FirstMoment--')
    print(res)
    print(op.doQuadrature(alpha, bounds))
    print(res - op.doQuadrature(alpha, bounds))


<mask token>

<mask token>


def test_base():
    n = 100000
    b = pysgpp.SLinearModifiedBase()
    for l in range(1, 3, 1):
        print(l)
        for i in range(1, 2 ** l, 1):
            s = 1
            s = sum([b.eval(l, i, x) for x in np.linspace(0, 1, n)]) / n
            print(('Level:', l, 'index:', i))
            epsilon = abs(s - b.getIntegral(l, i))
            if epsilon > 10 ** -6:
                print('s: {} res: {}'.format(s, b.getIntegral(l, i)))
                print('error: {}'.format(epsilon))


<mask token>


def test_LTwoDotImplicit(grid, l):
    grid.getGenerator().regular(l)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    m = pysgpp.DataMatrix(size, size)
    opExplicit = pysgpp.createOperationLTwoDotExplicit(m, grid)
    op = pysgpp.createOperationLTwoDotProduct(grid)
    alpha = pysgpp.DataVector(size)
    resultExplicit = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)
    for i in range(size):
        alpha[i] = 1
    opExplicit.mult(alpha, resultExplicit)
    op.mult(alpha, result)
    for i in range(size):
        if result[i] != resultExplicit[i]:
            print('Error result entry {} differs'.format(i))
        if abs(result[i] - resultExplicit[i]) > 1e-16:
            print('result:{}'.format(result[i]))
            print('resultExplicit:{}'.format(resultExplicit[i]))


def test_laplace(grid, lmax):
    resolution = 100000
    grid.getGenerator().regular(lmax)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    b = getBasis(grid)
    op = pysgpp.createOperationLaplace(grid)
    alpha = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)
    for point_i in range(size):
        for point_j in range(size):
            gp_i = gridStorage.getPoint(point_i)
            gp_j = gridStorage.getPoint(point_j)
            print('--------')
            for i in range(0, size):
                alpha[i] = 0
            alpha[point_i] = 1
            op.mult(alpha, result)
            xs = np.linspace(0, 1, resolution)
            approx = sum([(b.evalDx(gp_i.getLevel(0), gp_i.getIndex(0), x) *
                b.evalDx(gp_j.getLevel(0), gp_j.getIndex(0), x)) for x in xs]
                ) / resolution
            print('i,j: {},{} result: {} approx:{}'.format(point_i, point_j,
                result[point_j], approx))
            if abs(result.get(point_j) - approx) > 0.1:
                print('--------')
                print('points: {},{} '.format(point_i, point_j))
                print('approx:{}'.format(approx))
                print('result:{}'.format(result.get(point_j)))
                print('--------')


def test_poly_evaldx():
    l = 3
    i = 1
    x = 0.12
    eps = 0.0001
    b = pysgpp.SPolyModifiedClenshawCurtisBase(3)
    tang = b.evalDx(l, i, x)
    sec = (b.eval(l, i, x + eps) - b.eval(l, i, x - eps)) / (2 * eps)
    print('evalDx:{}'.format(tang))
    print('sekante:{}'.format(sec))
    print('evals: {} {}'.format(b.eval(l, i, x - eps), b.eval(l, i, x + eps)))


def plot_evaldx():
    l = 3
    i = 1
    xs = np.linspace(0, 1, 50)
    b = pysgpp.SPolyClenshawCurtisBase(3)
    plt.plot(xs, [b.evalDx(l, i, x) for x in xs])
    plt.show()


<mask token>


def test_firstMoment(grid, lmax):
    grid.getGenerator().regular(lmax)
    resolution = 100000
    gridStorage = grid.getStorage()
    b = grid.getBasis()
    op = pysgpp.createOperationFirstMoment(grid)
    alpha = pysgpp.DataVector(grid.getSize(), 1.0)
    bounds = pysgpp.DataMatrix(1, 2, 0.0)
    bounds.set(0, 1, 1.0)
    res = 0.0
    for i in range(grid.getSize()):
        lev = gridStorage.getPoint(i).getLevel(0)
        ind = gridStorage.getPoint(i).getIndex(0)
        temp_res = 0.0
        for c in range(resolution):
            x = float(c) / resolution
            temp_res += x * b.eval(lev, ind, x)
        res += alpha.get(i) * temp_res / resolution
    print('--FirstMoment--')
    print(res)
    print(op.doQuadrature(alpha, bounds))
    print(res - op.doQuadrature(alpha, bounds))


<mask token>

<mask token>


def test_base():
    n = 100000
    b = pysgpp.SLinearModifiedBase()
    for l in range(1, 3, 1):
        print(l)
        for i in range(1, 2 ** l, 1):
            s = 1
            s = sum([b.eval(l, i, x) for x in np.linspace(0, 1, n)]) / n
            print(('Level:', l, 'index:', i))
            epsilon = abs(s - b.getIntegral(l, i))
            if epsilon > 10 ** -6:
                print('s: {} res: {}'.format(s, b.getIntegral(l, i)))
                print('error: {}'.format(epsilon))


def test_LTwoDot(grid, l):
    res = 10000
    b = grid.getBasis()
    grid.getGenerator().regular(l)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    m = pysgpp.DataMatrix(size, size)
    opMatrix = pysgpp.createOperationLTwoDotExplicit(m, grid)
    for i in range(gridStorage.getSize()):
        for j in range(i, gridStorage.getSize()):
            gpi = gridStorage.getPoint(i)
            gpj = gridStorage.getPoint(j)
            sol = 1
            for k in range(d):
                lik = gpi.getLevel(k)
                iik = gpi.getIndex(k)
                ljk = gpj.getLevel(k)
                ijk = gpj.getIndex(k)
                xs = np.linspace(0, 1, res)
                tmp = sum([(b.eval(lik, iik, x) * b.eval(ljk, ijk, x)) for
                    x in xs]) / res
                sol *= tmp
            error = abs(m.get(i, j) - sol)
            if error >= 10 ** -4:
                print('i:{} j:{} error: {}'.format(i, j, error))
                print('iik:{} lik:{} ijk:{} ljk:{} error: {}'.format(iik,
                    lik, ijk, ljk, error))
                print('is:{} should:{}'.format(m.get(i, j), sol))


def test_LTwoDotImplicit(grid, l):
    grid.getGenerator().regular(l)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    m = pysgpp.DataMatrix(size, size)
    opExplicit = pysgpp.createOperationLTwoDotExplicit(m, grid)
    op = pysgpp.createOperationLTwoDotProduct(grid)
    alpha = pysgpp.DataVector(size)
    resultExplicit = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)
    for i in range(size):
        alpha[i] = 1
    opExplicit.mult(alpha, resultExplicit)
    op.mult(alpha, result)
    for i in range(size):
        if result[i] != resultExplicit[i]:
            print('Error result entry {} differs'.format(i))
        if abs(result[i] - resultExplicit[i]) > 1e-16:
            print('result:{}'.format(result[i]))
            print('resultExplicit:{}'.format(resultExplicit[i]))


def test_laplace(grid, lmax):
    resolution = 100000
    grid.getGenerator().regular(lmax)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    b = getBasis(grid)
    op = pysgpp.createOperationLaplace(grid)
    alpha = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)
    for point_i in range(size):
        for point_j in range(size):
            gp_i = gridStorage.getPoint(point_i)
            gp_j = gridStorage.getPoint(point_j)
            print('--------')
            for i in range(0, size):
                alpha[i] = 0
            alpha[point_i] = 1
            op.mult(alpha, result)
            xs = np.linspace(0, 1, resolution)
            approx = sum([(b.evalDx(gp_i.getLevel(0), gp_i.getIndex(0), x) *
                b.evalDx(gp_j.getLevel(0), gp_j.getIndex(0), x)) for x in xs]
                ) / resolution
            print('i,j: {},{} result: {} approx:{}'.format(point_i, point_j,
                result[point_j], approx))
            if abs(result.get(point_j) - approx) > 0.1:
                print('--------')
                print('points: {},{} '.format(point_i, point_j))
                print('approx:{}'.format(approx))
                print('result:{}'.format(result.get(point_j)))
                print('--------')


def test_poly_evaldx():
    l = 3
    i = 1
    x = 0.12
    eps = 0.0001
    b = pysgpp.SPolyModifiedClenshawCurtisBase(3)
    tang = b.evalDx(l, i, x)
    sec = (b.eval(l, i, x + eps) - b.eval(l, i, x - eps)) / (2 * eps)
    print('evalDx:{}'.format(tang))
    print('sekante:{}'.format(sec))
    print('evals: {} {}'.format(b.eval(l, i, x - eps), b.eval(l, i, x + eps)))


def plot_evaldx():
    l = 3
    i = 1
    xs = np.linspace(0, 1, 50)
    b = pysgpp.SPolyClenshawCurtisBase(3)
    plt.plot(xs, [b.evalDx(l, i, x) for x in xs])
    plt.show()


<mask token>


def test_firstMoment(grid, lmax):
    grid.getGenerator().regular(lmax)
    resolution = 100000
    gridStorage = grid.getStorage()
    b = grid.getBasis()
    op = pysgpp.createOperationFirstMoment(grid)
    alpha = pysgpp.DataVector(grid.getSize(), 1.0)
    bounds = pysgpp.DataMatrix(1, 2, 0.0)
    bounds.set(0, 1, 1.0)
    res = 0.0
    for i in range(grid.getSize()):
        lev = gridStorage.getPoint(i).getLevel(0)
        ind = gridStorage.getPoint(i).getIndex(0)
        temp_res = 0.0
        for c in range(resolution):
            x = float(c) / resolution
            temp_res += x * b.eval(lev, ind, x)
        res += alpha.get(i) * temp_res / resolution
    print('--FirstMoment--')
    print(res)
    print(op.doQuadrature(alpha, bounds))
    print(res - op.doQuadrature(alpha, bounds))


<mask token>

try:
    import pysgpp
    from pysgpp.extensions.datadriven.uq.operations.sparse_grid import getBasis
    import numpy as np
    import random
    import matplotlib.pyplot as plt
except ImportError as e:
    print(e.__class__.__name__ + ': ' + e.msg)
    print('Skipping example...')
    exit(0)


def test_base():
    n = 100000
    b = pysgpp.SLinearModifiedBase()
    for l in range(1, 3, 1):
        print(l)
        for i in range(1, 2 ** l, 1):
            s = 1
            s = sum([b.eval(l, i, x) for x in np.linspace(0, 1, n)]) / n
            print(('Level:', l, 'index:', i))
            epsilon = abs(s - b.getIntegral(l, i))
            if epsilon > 10 ** -6:
                print('s: {} res: {}'.format(s, b.getIntegral(l, i)))
                print('error: {}'.format(epsilon))


def test_LTwoDot(grid, l):
    res = 10000
    b = grid.getBasis()
    grid.getGenerator().regular(l)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    m = pysgpp.DataMatrix(size, size)
    opMatrix = pysgpp.createOperationLTwoDotExplicit(m, grid)
    for i in range(gridStorage.getSize()):
        for j in range(i, gridStorage.getSize()):
            gpi = gridStorage.getPoint(i)
            gpj = gridStorage.getPoint(j)
            sol = 1
            for k in range(d):
                lik = gpi.getLevel(k)
                iik = gpi.getIndex(k)
                ljk = gpj.getLevel(k)
                ijk = gpj.getIndex(k)
                xs = np.linspace(0, 1, res)
                tmp = sum([(b.eval(lik, iik, x) * b.eval(ljk, ijk, x)) for
                    x in xs]) / res
                sol *= tmp
            error = abs(m.get(i, j) - sol)
            if error >= 10 ** -4:
                print('i:{} j:{} error: {}'.format(i, j, error))
                print('iik:{} lik:{} ijk:{} ljk:{} error: {}'.format(iik,
                    lik, ijk, ljk, error))
                print('is:{} should:{}'.format(m.get(i, j), sol))


def test_LTwoDotImplicit(grid, l):
    grid.getGenerator().regular(l)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    m = pysgpp.DataMatrix(size, size)
    opExplicit = pysgpp.createOperationLTwoDotExplicit(m, grid)
    op = pysgpp.createOperationLTwoDotProduct(grid)
    alpha = pysgpp.DataVector(size)
    resultExplicit = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)
    for i in range(size):
        alpha[i] = 1
    opExplicit.mult(alpha, resultExplicit)
    op.mult(alpha, result)
    for i in range(size):
        if result[i] != resultExplicit[i]:
            print('Error result entry {} differs'.format(i))
        if abs(result[i] - resultExplicit[i]) > 1e-16:
            print('result:{}'.format(result[i]))
            print('resultExplicit:{}'.format(resultExplicit[i]))


def test_laplace(grid, lmax):
    resolution = 100000
    grid.getGenerator().regular(lmax)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    b = getBasis(grid)
    op = pysgpp.createOperationLaplace(grid)
    alpha = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)
    for point_i in range(size):
        for point_j in range(size):
            gp_i = gridStorage.getPoint(point_i)
            gp_j = gridStorage.getPoint(point_j)
            print('--------')
            for i in range(0, size):
                alpha[i] = 0
            alpha[point_i] = 1
            op.mult(alpha, result)
            xs = np.linspace(0, 1, resolution)
            approx = sum([(b.evalDx(gp_i.getLevel(0), gp_i.getIndex(0), x) *
                b.evalDx(gp_j.getLevel(0), gp_j.getIndex(0), x)) for x in xs]
                ) / resolution
            print('i,j: {},{} result: {} approx:{}'.format(point_i, point_j,
                result[point_j], approx))
            if abs(result.get(point_j) - approx) > 0.1:
                print('--------')
                print('points: {},{} '.format(point_i, point_j))
                print('approx:{}'.format(approx))
                print('result:{}'.format(result.get(point_j)))
                print('--------')


def test_poly_evaldx():
    l = 3
    i = 1
    x = 0.12
    eps = 0.0001
    b = pysgpp.SPolyModifiedClenshawCurtisBase(3)
    tang = b.evalDx(l, i, x)
    sec = (b.eval(l, i, x + eps) - b.eval(l, i, x - eps)) / (2 * eps)
    print('evalDx:{}'.format(tang))
    print('sekante:{}'.format(sec))
    print('evals: {} {}'.format(b.eval(l, i, x - eps), b.eval(l, i, x + eps)))


def plot_evaldx():
    l = 3
    i = 1
    xs = np.linspace(0, 1, 50)
    b = pysgpp.SPolyClenshawCurtisBase(3)
    plt.plot(xs, [b.evalDx(l, i, x) for x in xs])
    plt.show()


def plot_evaldx_prod(grid, lmax, i, j):
    grid.getGenerator().regular(lmax)
    gridStorage = grid.getStorage()
    b = getBasis(grid)
    gp_i = gridStorage.getPoint(i)
    gp_j = gridStorage.getPoint(j)
    lik = gp_i.getLevel(0)
    ljk = gp_j.getLevel(0)
    iik = gp_i.getIndex(0)
    ijk = gp_j.getIndex(0)
    xs = np.linspace(0, 1, 100000)
    plt.plot(xs, [(b.evalDx(lik, iik, x) * b.evalDx(ljk, ijk, x)) for x in xs])
    plt.show()


def test_laplace2(grid, lmax):
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    m = pysgpp.DataMatrix(size, size)
    op = pysgpp.createOperationLaplaceExplicit(m, grid)
    print(m)


def test_firstMoment(grid, lmax):
    grid.getGenerator().regular(lmax)
    resolution = 100000
    gridStorage = grid.getStorage()
    b = grid.getBasis()
    op = pysgpp.createOperationFirstMoment(grid)
    alpha = pysgpp.DataVector(grid.getSize(), 1.0)
    bounds = pysgpp.DataMatrix(1, 2, 0.0)
    bounds.set(0, 1, 1.0)
    res = 0.0
    for i in range(grid.getSize()):
        lev = gridStorage.getPoint(i).getLevel(0)
        ind = gridStorage.getPoint(i).getIndex(0)
        temp_res = 0.0
        for c in range(resolution):
            x = float(c) / resolution
            temp_res += x * b.eval(lev, ind, x)
        res += alpha.get(i) * temp_res / resolution
    print('--FirstMoment--')
    print(res)
    print(op.doQuadrature(alpha, bounds))
    print(res - op.doQuadrature(alpha, bounds))


def test_secondtMoment(grid, lmax):
    gridStorage = grid.getStorage()
    gridStorage.clear()
    grid.getGenerator().regular(lmax)
    resolution = 1000000
    b = grid.getBasis()
    op = pysgpp.createOperationSecondMoment(grid)
    alpha = pysgpp.DataVector(grid.getSize(), 1.0)
    bounds = pysgpp.DataMatrix(1, 2, 0.0)
    bounds.set(0, 1, 1.0)
    res = 0.0
    for i in range(grid.getSize()):
        lev = gridStorage.getPoint(i).getLevel(0)
        ind = gridStorage.getPoint(i).getIndex(0)
        temp_res = 0.0
        for c in range(resolution):
            x = float(c) / resolution
            temp_res += x * x * b.eval(lev, ind, x)
        res += alpha.get(i) * temp_res / resolution
    print('--SecondMoment--')
    print(res)
    print(op.doQuadrature(alpha, bounds))
    print(res - op.doQuadrature(alpha, bounds))


d = 1
l = 3
grid = pysgpp.Grid.createModLinearGrid(d)
test_firstMoment(grid, l)
test_secondtMoment(grid, l)

# Copyright (C) 2008-today The SG++ project
# This file is part of the SG++ project. For conditions of distribution and
# use, please see the copyright notice provided with SG++ or at
# sgpp.sparsegrids.org

try:
  import pysgpp
  from pysgpp.extensions.datadriven.uq.operations.sparse_grid import getBasis
  import numpy as np
  import random
  import matplotlib.pyplot as plt
  
except ImportError as e:
    print(e.__class__.__name__ + ": " + e.msg)
    print("Skipping example...")
    exit(0)

def test_base():
    n = 100000
    b = pysgpp.SLinearModifiedBase()
    # print(sum([b.eval(0,0,x) for x in np.linspace(0, 1, n)]) / n)
    for l in range(1,3,1):
        print(l)
        for i in range(1, 2**l, 1):
            s = 1
            s = sum([b.eval(l,i,x) for x in np.linspace(0, 1, n)]) / n
            print(("Level:",l , "index:", i))
            epsilon = abs(s - b.getIntegral(l,i))
            if(epsilon > 10**-6):
                print("s: {} res: {}".format(s, b.getIntegral(l,i)))
                print("error: {}".format(epsilon))

def test_LTwoDot(grid, l):
    res = 10000
    b = grid.getBasis();
    grid.getGenerator().regular(l)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    # print(size)
    m = pysgpp.DataMatrix(size, size)
    opMatrix = pysgpp.createOperationLTwoDotExplicit(m, grid)
    # print m
    # m_comp = pysgpp.DataMatrix(size, size)
    for i in range(gridStorage.getSize()):
        for j in range(i, gridStorage.getSize()):
            gpi = gridStorage.getPoint(i)
            gpj = gridStorage.getPoint(j)
            sol = 1
            # print "--------"
            # print "i:{} j:{}".format(i, j)
            for k in range(d):
                lik = gpi.getLevel(k)
                iik = gpi.getIndex(k)
                ljk = gpj.getLevel(k)
                ijk = gpj.getIndex(k)
                # print "i l,i: {},{}   j l,i: {},{}".format(lik, iik, ljk, ijk)
                xs = np.linspace(0, 1, res)
                tmp = sum([b.eval(lik, iik, x)*b.eval(ljk, ijk, x) for x in xs]) / res
                sol *= tmp
                # print("lik:{} iik:{} ljk:{} ijk:{} k:{} tmp: {}".format(lik, iik, ljk, ijk, k,tmp))
            # print(sol)
            error = abs(m.get(i,j) - sol)
            # print error
            if(error >= 10**-4):
                print("i:{} j:{} error: {}".format(i, j, error))
                print("iik:{} lik:{} ijk:{} ljk:{} error: {}".format(iik, lik, ijk, ljk, error))
                print("is:{} should:{}".format(m.get(i,j), sol))

def test_LTwoDotImplicit(grid, l):
    grid.getGenerator().regular(l)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    # print(size)
    m = pysgpp.DataMatrix(size, size)
    opExplicit = pysgpp.createOperationLTwoDotExplicit(m, grid)
    op = pysgpp.createOperationLTwoDotProduct(grid)
    alpha = pysgpp.DataVector(size)
    resultExplicit = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)
    for i in range(size):
        alpha[i] = 1
    opExplicit.mult(alpha, resultExplicit)
    op.mult(alpha, result)
    for i in range(size):
        if result[i] != resultExplicit[i]:
            print("Error result entry {} differs".format(i))

        if abs(result[i] - resultExplicit[i]) > 1e-16:
            # print result[i] - resultExplicit[i]
            print("result:{}".format(result[i]))
            print("resultExplicit:{}".format(resultExplicit[i]))

def test_laplace(grid, lmax):
    resolution = 100000
    grid.getGenerator().regular(lmax)
    gridStorage = grid.getStorage()
    size = gridStorage.getSize()
    b = getBasis(grid)
    op = pysgpp.createOperationLaplace(grid)
    alpha = pysgpp.DataVector(size)
    result = pysgpp.DataVector(size)

    for point_i in range(size):
      for point_j in range(size):
        gp_i = gridStorage.getPoint(point_i)
        gp_j = gridStorage.getPoint(point_j)
        print("--------")
        for i in range(0, size):
          alpha[i] = 0
        alpha[point_i] = 1
        op.mult(alpha, result)
        xs = np.linspace(0, 1, resolution)
        approx = sum([b.evalDx(gp_i.getLevel(0), gp_i.getIndex(0), x) * b.evalDx(gp_j.getLevel(0), gp_j.getIndex(0), x) for x in xs]) / resolution
        print("i,j: {},{} result: {} approx:{}".format(point_i, point_j, result[point_j], approx))
        if(abs(result.get(point_j) - approx) > 1e-1):
          print("--------")
          print("points: {},{} ".format(point_i, point_j))
          print("approx:{}".format(approx))
          print("result:{}".format(result.get(point_j)))
          # print result
          print("--------")

def test_poly_evaldx():
    l = 3
    i = 1
    x = 0.12
    eps = 0.0001
    b = pysgpp.SPolyModifiedClenshawCurtisBase(3)
    tang = b.evalDx(l, i, x)
    sec = (b.eval(l, i, x + eps) -  b.eval(l, i, x - eps)) / (2*eps)
    print("evalDx:{}".format(tang))
    print("sekante:{}".format(sec))
    print("evals: {} {}".format( b.eval(l, i, x - eps), b.eval(l, i, x + eps) ))

def plot_evaldx():
    l = 3
    i = 1
    xs = np.linspace(0, 1, 50)
    b = pysgpp.SPolyClenshawCurtisBase(3)
    plt.plot(xs, [b.evalDx(l, i, x) for x in xs])
    plt.show()

def plot_evaldx_prod(grid, lmax, i, j):
    grid.getGenerator().regular(lmax)
    gridStorage = grid.getStorage()
    b = getBasis(grid)
    gp_i = gridStorage.getPoint(i)
    gp_j = gridStorage.getPoint(j)
    lik = gp_i.getLevel(0)
    ljk = gp_j.getLevel(0)
    iik = gp_i.getIndex(0)
    ijk = gp_j.getIndex(0)
    xs = np.linspace(0, 1, 100000)

    plt.plot(xs, [b.evalDx(lik, iik, x) * b.evalDx(ljk, ijk, x) for x in xs])
    plt.show()


def test_laplace2(grid, lmax):
  # grid.getGenerator().regular(lmax)
  gridStorage = grid.getStorage()
  size = gridStorage.getSize()
  m = pysgpp.DataMatrix(size, size)
  op = pysgpp.createOperationLaplaceExplicit(m, grid)
  print(m)



def test_firstMoment(grid, lmax):
  grid.getGenerator().regular(lmax)
  resolution = 100000
  gridStorage = grid.getStorage()
  b = grid.getBasis()
  op = pysgpp.createOperationFirstMoment(grid)
  alpha = pysgpp.DataVector(grid.getSize(), 1.0)
  bounds = pysgpp.DataMatrix(1, 2, 0.0)
  bounds.set(0, 1, 1.0)
  res = 0.0
  for i in range(grid.getSize()):
    lev = gridStorage.getPoint(i).getLevel(0)
    ind = gridStorage.getPoint(i).getIndex(0)
    temp_res = 0.0
    for c in range(resolution):
        x = float(c) / resolution
        temp_res += x * b.eval(lev, ind, x)
    res += alpha.get(i) * temp_res / resolution
  print("--FirstMoment--")
  print(res)
  print(op.doQuadrature(alpha, bounds))
  print (res - op.doQuadrature(alpha, bounds))

def test_secondtMoment(grid, lmax):
  gridStorage = grid.getStorage()
  gridStorage.clear()
  grid.getGenerator().regular(lmax)
  resolution = 1000000
  b = grid.getBasis()
  op = pysgpp.createOperationSecondMoment(grid)
  alpha = pysgpp.DataVector(grid.getSize(), 1.0)
  bounds = pysgpp.DataMatrix(1, 2, 0.0)
  bounds.set(0, 1, 1.0)
  res = 0.0
  for i in range(grid.getSize()):
    lev = gridStorage.getPoint(i).getLevel(0)
    ind = gridStorage.getPoint(i).getIndex(0)
    temp_res = 0.0
    for c in range(resolution):
        x = float(c) / resolution
        temp_res += x * x *b.eval(lev, ind, x)
    res += alpha.get(i) * temp_res / resolution

  print("--SecondMoment--")
  print(res)
  print(op.doQuadrature(alpha, bounds))
  print (res - op.doQuadrature(alpha, bounds))

# test_poly_evaldx()
# plot_evaldx()
# test_base()
d = 1
l = 3
grid = pysgpp.Grid.createModLinearGrid(d)
# plot_evaldx_prod(grid, 4, 1, 4)
# test_laplace(grid, l)
# test_laplace2(grid, l)
# test_LTwoDot(grid, l)
# test_LTwoDotImplicit(grid, l)
test_firstMoment(grid, l)
test_secondtMoment(grid, l)
