← Back to branch summary

~daniele-bigoni/tensortoolbox/tt-docs

« back to all changes in this revision

Viewing changes to Development/Examples/Ordering/TTdmrgOrdering.py

  • Committer: Daniele Bigoni
  • Date: 2015-01-19 11:10:20 UTC
  • Revision ID: dabi@dtu.dk-20150119111020-p0uckg4ab3xqzf47
merged with research

Show diffs side-by-side

added added

removed removed

Lines of Context:
Development/Examples/Ordering/TTdmrgOrdering.py
 
1
#
 
2
# This file is part of TensorToolbox.
 
3
#
 
4
# TensorToolbox is free software: you can redistribute it and/or modify
 
5
# it under the terms of the LGNU Lesser General Public License as published by
 
6
# the Free Software Foundation, either version 3 of the License, or
 
7
# (at your option) any later version.
 
8
#
 
9
# TensorToolbox is distributed in the hope that it will be useful,
 
10
# but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
# LGNU Lesser General Public License for more details.
 
13
#
 
14
# You should have received a copy of the LGNU Lesser General Public License
 
15
# along with TensorToolbox.  If not, see <http://www.gnu.org/licenses/>.
 
16
#
 
17
# DTU UQ Library
 
18
# Copyright (C) 2014-2015 The Technical University of Denmark
 
19
# Scientific Computing Section
 
20
# Department of Applied Mathematics and Computer Science
 
21
#
 
22
# Author: Daniele Bigoni
 
23
#
 
24
 
 
25
import sys
 
26
 
 
27
import operator
 
28
import time
 
29
import logging
 
30
 
 
31
import numpy as np
 
32
import numpy.linalg as npla
 
33
import numpy.random as npr
 
34
 
 
35
from scipy import stats
 
36
import scipy.cluster.hierarchy as cluster
 
37
 
 
38
import networkx as nx
 
39
import openopt as oo
 
40
 
 
41
import TensorToolbox as DT
 
42
import TensorToolbox.multilinalg as mla
 
43
 
 
44
from SpectralToolbox import Spectral1D as S1D
 
45
from SpectralToolbox import SpectralND as SND
 
46
from SpectralToolbox import Misc
 
47
from UQToolbox import UncertaintyQuantification as UQ
 
48
from UQToolbox import CutANOVA
 
49
 
 
50
import matplotlib.pyplot as plt
 
51
from matplotlib import cm
 
52
from mpl_toolkits.mplot3d import Axes3D
 
53
 
 
54
npr.seed(1)
 
55
 
 
56
##########################################################
 
57
# TEST 0
 
58
#  Corner peak with interleaving dependencies
 
59
#
 
60
# TEST 1
 
61
#  Try function with given covariance (f4 in STT paper)
 
62
#  
 
63
 
 
64
logging.basicConfig(level=logging.INFO)
 
65
 
 
66
TESTS = [0]
 
67
IS_PLOTTING = True
 
68
STORE_FIG = True
 
69
FORMATS = ['pdf','png','eps']
 
70
 
 
71
MCestVarLimit = 1e-1
 
72
MCestMinIter = 100
 
73
MCestMaxIter = 1e4
 
74
MCstep = 10000
 
75
 
 
76
if 0 in TESTS:
 
77
    ##########################################################
 
78
    # TEST 0
 
79
 
 
80
    import itertools
 
81
 
 
82
    maxvoleps = 1e-4
 
83
    delta = 1e-4
 
84
    eps = 1e-10
 
85
    lr_maxit = 50
 
86
 
 
87
    xspan = [0.,1.]
 
88
    d = 20
 
89
    size1D = 7
 
90
    vol = (xspan[1]-xspan[0])**d
 
91
 
 
92
    def f(X,params):
 
93
        groups = params['groups']
 
94
        cs = params['cs']
 
95
 
 
96
        if isinstance(X,np.ndarray):
 
97
            ndim = X.ndim
 
98
            if ndim == 1:
 
99
                X = np.array([X,])
 
100
        else:
 
101
            shape = None
 
102
 
 
103
        out = 1.
 
104
        for g in groups:
 
105
            out *= (1. + np.dot( X[:,g], cs[g] ) ) ** (-len(g)+1) # Check correct indexing of cs
 
106
        
 
107
        if ndim == 1:
 
108
            return out[0]
 
109
        else:
 
110
            return out
 
111
 
 
112
    # Construct groups
 
113
    dim_list = range(d)
 
114
    dim_list_shuffle = dim_list[:]
 
115
    npr.shuffle(dim_list_shuffle)
 
116
    groups = []
 
117
    groups_shuffle = []
 
118
    d_left = d
 
119
    while d_left > 0:
 
120
        partition = npr.randint(1,min(10,d_left+1))
 
121
        groups.append( dim_list[d-d_left:d-d_left+partition] )
 
122
        groups_shuffle.append( dim_list_shuffle[d-d_left:d-d_left+partition] )
 
123
        d_left -= partition
 
124
 
 
125
    print "Groups:         %s" % str(groups)    
 
126
    print "Groups_shuffle: %s" % str(groups_shuffle)
 
127
 
 
128
    # Generate coefficients
 
129
    unif_dist = stats.uniform()
 
130
    cs_shuffle = unif_dist.rvs(d)
 
131
    cs = cs_shuffle[ list(itertools.chain( *groups_shuffle )) ]
 
132
    
 
133
    # Prepare grid
 
134
    X = [ (S1D.JACOBI,S1D.GAUSS,(0.,0.),[0.,1.]) ] * d
 
135
    
 
136
    # Build function wrapper (shuffled)
 
137
    params_shuffle = {'groups': groups_shuffle,
 
138
                      'cs'    : cs_shuffle }
 
139
 
 
140
    # Build function wrapper (reordered)
 
141
    params = {'groups': groups,
 
142
              'cs'    : cs }
 
143
 
 
144
    ###############################################
 
145
    # Test Surrogate construction on Ordered
 
146
    STTapprox = DT.SQTT(f, X, params, eps=eps, method='ttdmrg', 
 
147
                        surrogateONOFF=True, surrogate_type=DT.PROJECTION, orders=[size1D]*d)
 
148
    STTapprox.build()
 
149
 
 
150
    print "TTdmrg: grid size %d" % STTapprox.TW.get_size()
 
151
    print "TTdmrg: function evaluations %d" % STTapprox.TW.get_fill_level()
 
152
 
 
153
    # Check approximation error using MC
 
154
    var = 1.
 
155
    dist = stats.uniform()
 
156
    DIST = UQ.MultiDimDistribution([dist] * d)
 
157
    intf = []
 
158
    values = []
 
159
    while (len(values) < MCestMinIter or var > mean**2. * MCestVarLimit) and len(values) < MCestMaxIter:
 
160
        # Monte Carlo
 
161
        xx = DIST.rvs(MCstep)
 
162
(??)
 
163
(??)        VsI = None
 
164
(??)        VsI = [P.GradVandermonde1D(xx[:,i]*2.-1.,size1D,0,norm=True) for i in range(d)]
 
165
(??)        TTval = TT_four.interpolate(VsI)
 
166
(??)
 
167
(??)        TTvals = [ TTval[tuple([i]*d)] for i in range(MCstep) ]
 
168
 
 
169
        fval = f(xx,params)
 
170
 
 
171
        intf.extend( list(fval**2.) )
 
172
        values.extend( [(fval[i]-TTvals[i])**2. for i in range(MCstep)])
 
173
        mean = vol * np.mean(values)
 
174
        var = vol**2. * np.var(values) / len(values)
 
175
 
 
176
        sys.stdout.write("L2err estim. iter: %d Var: %e VarLim: %e \r" % (len(values) , var, mean**2. * MCestVarLimit))
 
177
        sys.stdout.flush()
 
178
 
 
179
    sys.stdout.write("\n")
 
180
    sys.stdout.flush()
 
181
    L2err = np.sqrt(mean/np.mean(intf))
 
182
    print "TTdmrg: L2 err TTapprox: %e" % L2err
 
183
 
 
184
    ###############################################
 
185
    # Test Surrogate construction on Shuffled
 
186
    STTapprox_shuffle = DT.SQTT(f, X, params_shuffle,eps=eps,method='ttdmrg',
 
187
                                surrogateONOFF=True, surrogate_type=DT.PROJECTION, orders=[size1D]*d)
 
188
    # STTapprox_shuffle.build()
 
189
 
 
190
    # print "TTdmrg_shuffle: grid size %d" % TW_shuffle.get_size()
 
191
    # print "TTdmrg_shuffle: function evaluations %d" % TW_shuffle.get_fill_level()
 
192
 
 
193
    # # Compute Fourier coefficients TT
 
194
    # TT_four_shuffle = TTapprox_shuffle.project(V,W)
 
195
 
 
196
    # # Check approximation error using MC
 
197
    # var = 1.
 
198
    # dist = stats.uniform()
 
199
    # DIST = UQ.MultiDimDistribution([dist] * d)
 
200
    # intf = []
 
201
    # values = []
 
202
    # while (len(values) < MCestMinIter or var > mean**2. * MCestVarLimit) and len(values) < MCestMaxIter:
 
203
    #     # Monte Carlo
 
204
    #     xx = DIST.rvs(MCstep)
 
205
 
 
206
    #     VsI = None
 
207
    #     VsI = [P.GradVandermonde1D(xx[:,i]*2.-1.,size1D,0,norm=True) for i in range(d)]
 
208
    #     TTval = TT_four_shuffle.interpolate(VsI)
 
209
 
 
210
    #     TTvals = [ TTval[tuple([i]*d)] for i in range(MCstep) ]
 
211
 
 
212
    #     fval = f(xx,params_shuffle)
 
213
 
 
214
    #     intf.extend( list(fval**2.) )
 
215
    #     values.extend( [(fval[i]-TTvals[i])**2. for i in range(MCstep)])
 
216
    #     mean = vol * np.mean(values)
 
217
    #     var = vol**2. * np.var(values) / len(values)
 
218
 
 
219
    #     sys.stdout.write("L2err estim. iter: %d Var: %e VarLim: %e \r" % (len(values) , var, mean**2. * MCestVarLimit))
 
220
    #     sys.stdout.flush()
 
221
 
 
222
    # sys.stdout.write("\n")
 
223
    # sys.stdout.flush()
 
224
    # L2err_shuffle = np.sqrt(mean/np.mean(intf))
 
225
    # print "TTdmrg_shuffle: L2 err TTapprox: %e" % L2err_shuffle
 
226
 
 
227
    
 
228
    # #####################################################
 
229
    # # Construct Cut-ANOVA-HDMR expansion ORD 2
 
230
    # CUT_ORDER = 2
 
231
    # Ns = [ 10 ] * d
 
232
    # polys = [ S1D.Poly1D(S1D.JACOBI,[0.,0.]) ] * d
 
233
    # tol = 2. * Misc.machineEpsilon()
 
234
    # X_cut = np.zeros((1,d))
 
235
    # cut_HDMR = CutANOVA.CutHDMR(polys,Ns,CUT_ORDER,X_cut,tol)
 
236
    # def transformFunc(X):
 
237
    #     # from [-1,1] to [0,1]
 
238
    #     return (X+1.)/2.
 
239
    # def f_hdmr(X):
 
240
    #     return f(X,params_shuffle)
 
241
    # cut_HDMR.evaluateFun(f_hdmr,transformFunc)
 
242
    # cut_HDMR.computeCutHDMR()
 
243
    # cut_HDMR.computeANOVA_HDMR()
 
244
    # # Fill correlation matrix
 
245
    # Sigma = np.zeros((d,d))
 
246
    # k = 0
 
247
    # minval = np.inf
 
248
    # maxval = 0.
 
249
    # for i in range(d-1):
 
250
    #     for j in range(i+1,d):
 
251
    #         Sigma[i,j] = np.dot( cut_HDMR.grids[2][k].ANOVA_HDMR_vals**2., cut_HDMR.grids[2][k].WF )
 
252
    #         if minval > Sigma[i,j]: minval = Sigma[i,j]
 
253
    #         if maxval < Sigma[i,j]: maxval = Sigma[i,j]
 
254
    #         k += 1
 
255
    # Sigma += Sigma.T
 
256
    # # Plot covariance matrix
 
257
    # plt.figure()
 
258
    # plt.imshow(Sigma, interpolation='none', origin='lower')
 
259
    # plt.colorbar()
 
260
    # plt.title('Estimated Covariance')
 
261
    # plt.show(block=False)
 
262
    # # Plot correlation network
 
263
    # import networkx as nx
 
264
    # plt.figure()
 
265
    # dt = [('weight', float)]
 
266
    # A = Sigma.view(dt)
 
267
    # G = nx.from_numpy_matrix(A)
 
268
    # pos=nx.circular_layout(G)
 
269
    # weights = [ G.get_edge_data(G.edges()[i][0],G.edges()[i][1])['weight'] for i in range(len(G.edges())) ]
 
270
    # nx.draw(G,pos,edge_color=weights,edge_cmap=cm.gray,edge_vmin=minval, edge_vmax=maxval, width=4)
 
271
    # plt.show(block=False) # display
 
272
 
 
273
    ##########################################################
 
274
    # Construct vicinity matrix
 
275
    TW_shuffle = STTapprox_shuffle.TW
 
276
    CutIdx = (size1D+1)//2
 
277
    
 
278
    Sigma = np.zeros((d,d))
 
279
    minval = np.inf
 
280
    maxval = 0.
 
281
    for i in range(d-1):
 
282
        for j in range(i+1,d):
 
283
            idxs = [ ( CutIdx if ((k != i) and (k != j)) else slice(None,None,None) ) \
 
284
                     for k in range(d) ]
 
285
            CutMat = TW_shuffle[ tuple(idxs) ]
 
286
 
 
287
            # Sigma[i,j] = npla.norm(CutMat,'fro')
 
288
            
 
289
            (u,s,v) = npla.svd(CutMat)
 
290
            eps_svd = 1e-10
 
291
            rank = (i for i,e in enumerate(s) if e < eps_svd * s[0]).next()
 
292
            Sigma[i,j] = float(rank)
 
293
 
 
294
            if minval > Sigma[i,j]: minval = Sigma[i,j]
 
295
            if maxval < Sigma[i,j]: maxval = Sigma[i,j]
 
296
    Sigma += Sigma.T
 
297
    # Plot matrix
 
298
    plt.figure(figsize=(6,5))
 
299
    plt.imshow(Sigma, interpolation='none', origin='lower')
 
300
    plt.colorbar()
 
301
    # plt.title('Frobenius norm')
 
302
    plt.tight_layout()
 
303
    plt.show(block=False)
 
304
    if STORE_FIG:
 
305
        for ff in FORMATS:
 
306
            plt.savefig('figs/VicinityMatrix.'+ff,format=ff)
 
307
    # Plot correlation network
 
308
    plt.figure(figsize=(6,5))
 
309
    dt = [('weight', float)]
 
310
    sigtmp = Sigma.copy()
 
311
    np.fill_diagonal(sigtmp,np.inf)
 
312
    A = ( np.ones(Sigma.shape) / (sigtmp/maxval) ).view(dt)
 
313
    G = nx.from_numpy_matrix(A)
 
314
    pos=nx.circular_layout(G)
 
315
    weights = [ G.get_edge_data(G.edges()[i][0],G.edges()[i][1])['weight'] for i in range(len(G.edges())) ]
 
316
    nx.draw(G,pos,edge_color=weights,edge_cmap=cm.gray, width=2)
 
317
    plt.show(block=False) # display
 
318
    if STORE_FIG:
 
319
        for ff in FORMATS:
 
320
            plt.savefig('figs/ConnectivityNetwork.'+ff,format=ff)
 
321
 
 
322
    # # Solve Traveling Sales Man problem on the network 
 
323
    # # (Need to install openopt, FuncDesigner, cvxopt, [glpk, python-glpk])
 
324
    # problem = oo.TSP(G, objective='weight')
 
325
    # tsp_sol = problem.solve('glpk')
 
326
 
 
327
    # tspG = nx.Graph()
 
328
    # tspG.add_edges_from(tsp_sol.Edges)
 
329
    # tspW = [ tsp_sol.Edges[i][2]['weight'] for i in range(len(tsp_sol.Edges)) ]
 
330
    # plt.figure()
 
331
    # nx.draw(tspG,pos,edge_color=tspW,edge_cmap=cm.gray, width=2)
 
332
    # plt.show(block=False) # display
 
333
    # # Determine the reordering: Leave out the weakest connection (lowest frobenious norm correspond to max edge)
 
334
    # tspW = [ tsp_sol.Edges[i][2]['weight'] for i in range(len(tsp_sol.Edges)) ]
 
335
    # nodes = tsp_sol.nodes
 
336
    # max_edge = np.argmax(tspW)
 
337
    # shuffle = nodes[max_edge+1:] + nodes[1:max_edge+1]
 
338
 
 
339
    # Solve clustering problem
 
340
    tmp = np.ones(Sigma.shape) / (sigtmp/maxval)
 
341
    cls = cluster.linkage(tmp, method='complete')
 
342
    plt.figure(figsize=(6,5))
 
343
    cluster.dendrogram(cls)
 
344
    plt.tight_layout()
 
345
    plt.show(block=False)
 
346
    if STORE_FIG:
 
347
        for ff in FORMATS:
 
348
            plt.savefig('figs/Clustering.'+ff,format=ff)
 
349
    cls_groups = cluster.fcluster(cls,0.5*np.max(tmp),'distance')
 
350
    cls_idx_groups = [[] for i in range(np.max(cls_groups)) ]
 
351
    for i,c in enumerate(cls_groups): cls_idx_groups[c-1].append(i)
 
352
    shuffle = list(itertools.chain( *cls_idx_groups ))
 
353
 
 
354
    # Print shuffling
 
355
    print "groups_shuffle:   %s" % (str(groups_shuffle))
 
356
    print "dim_list_shuffle: %s" % (str(dim_list_shuffle))
 
357
    print "cluster_groups:   %s" % (str(cls_idx_groups))
 
358
    print "shuffle:          %s" % (str(shuffle))
 
359
    
 
360
    # Define the new Tensor Wrapper for the unshuffled function
 
361
    # unsfl_idxs = np.argsort(dim_list_shuffle) # To be changed with the shuffle
 
362
    unsfl_idxs = np.argsort(shuffle)
 
363
    X_usfl = [ np.arange(size1D+1,dtype=int) for i in range(d) ]
 
364
    def f_usfl(X,params):
 
365
        usfl = params['unshuffle']
 
366
        return f( X[usfl], params_shuffle )
 
367
    params_usfl = {'unshuffle': unsfl_idxs}
 
368
    TW_usfl = DT.TensorWrapper(f_usfl, X_usfl, params_usfl, twtype='view')
 
369
    
 
370
    STTapprox_unshuffle = DT.SQTT(f_usfl, X, params_usfl, eps=eps,method='ttdmrg',
 
371
                                  surrogateONOFF=True, surrogate_type=DT.PROJECTION, orders=[size1D]*d)
 
372
    STTapprox_unshuffle.build()
 
373
    
 
374
    print "TTdmrg_unshuffle: grid size %d" % TW_shuffle.get_size()
 
375
    print "TTdmrg_unshuffle: function evaluations %d" % TW_shuffle.get_fill_level()
 
376
 
 
377
    # Check approximation error using MC
 
378
    var = 1.
 
379
    dist = stats.uniform()
 
380
    DIST = UQ.MultiDimDistribution([dist] * d)
 
381
    intf = []
 
382
    values = []
 
383
    while (len(values) < MCestMinIter or var > mean**2. * MCestVarLimit) and len(values) < MCestMaxIter:
 
384
        # Monte Carlo
 
385
        xx = DIST.rvs(MCstep)
 
386
 
 
387
        TTvals = STTapprox_unshuffle(xx[:,shuffle])
 
388
        
 
389
        fval = f(xx,params_shuffle)
 
390
 
 
391
        intf.extend( list(fval**2.) )
 
392
        values.extend( [(fval[i]-TTvals[i])**2. for i in range(MCstep)])
 
393
        mean = vol * np.mean(values)
 
394
        var = vol**2. * np.var(values) / len(values)
 
395
 
 
396
        sys.stdout.write("L2err estim. iter: %d Var: %e VarLim: %e \r" % (len(values) , var, mean**2. * MCestVarLimit))
 
397
        sys.stdout.flush()
 
398
 
 
399
    sys.stdout.write("\n")
 
400
    sys.stdout.flush()
 
401
    L2err_shuffle = np.sqrt(mean/np.mean(intf))
 
402
    print "TTdmrg_unshuffle: L2 err TTapprox: %e" % L2err_shuffle
 
403
 
 
404
if 1 in TESTS:
 
405
    ##########################################################
 
406
    # TEST 1
 
407
 
 
408
    import itertools
 
409
 
 
410
    maxvoleps = 1e-10
 
411
    delta = 1e-10
 
412
    eps = 1e-13
 
413
 
 
414
    xspan = [-1.,1.]
 
415
    d = 6
 
416
    size1D = 10
 
417
    maxordF = 10
 
418
    vol = (xspan[1]-xspan[0])**d
 
419
 
 
420
    # Construct Sigma
 
421
 
 
422
    # Uniform
 
423
    # sigma_dist = stats.uniform()
 
424
    # Sigma = sigma_dist.rvs((d,d)) - 0.5
 
425
 
 
426
    # Beta 0.1,0.1 to push extreme vals
 
427
    # sigma_dist = stats.beta(0.5,0.5)
 
428
    # sigma_dist = stats.uniform(loc=-1,scale=2)
 
429
    # Sigma = np.zeros((d,d))
 
430
    # for i in range(d-1): Sigma[i,i+1:] = sigma_dist.rvs( d-1-i )
 
431
    # Sigma += Sigma.T
 
432
    # Sigma += np.eye(d) # * float(d)
 
433
 
 
434
    # sigma_dist = stats.beta(0.5,0.5)
 
435
    # Sigma1 = (sigma_dist.rvs((d,d))) * 2. - 1.
 
436
    # Sigma2 = (sigma_dist.rvs((d,d))) * 2. - 1.
 
437
    # Sigma = (Sigma1 + Sigma2)/2.
 
438
    # Sigma += np.eye(d) * float(d)
 
439
 
 
440
    sigma_dist = stats.uniform(loc=-1,scale=2)
 
441
    Sigma = sigma_dist.rvs((d,d))
 
442
    Sigma += Sigma.T
 
443
    Sigma += np.eye(d) * float(d)
 
444
 
 
445
    # Normalize maximum eigenvalue (direction of maximum decay)
 
446
    (val,vec) = npla.eig(Sigma)
 
447
    # w = w * stats.beta(0.5,0.5).rvs(len(w))
 
448
    val = (val / np.max(val)) # * stats.beta(0.01,0.01).rvs(len(w))
 
449
    Sigma = np.dot(vec,np.dot(np.diag(val),vec.T))
 
450
 
 
451
    # Sigma /= 3
 
452
 
 
453
    print Sigma
 
454
    
 
455
    # Plot Covariance Matrix
 
456
    sigtmp = Sigma.copy()
 
457
    np.fill_diagonal(sigtmp,0.)
 
458
    plt.figure()
 
459
    plt.imshow(np.abs(sigtmp), interpolation='none')
 
460
    plt.colorbar()
 
461
    plt.title('Exact Covariance')
 
462
    plt.show(block=False)
 
463
 
 
464
    # Plot network graph
 
465
    import networkx as nx
 
466
    import string
 
467
    plt.figure()
 
468
    dt = [('weight', float)]
 
469
    A = np.abs(sigtmp) / np.max(np.abs(sigtmp))
 
470
    A = A.view(dt)
 
471
    G = nx.from_numpy_matrix(A)
 
472
    pos=nx.circular_layout(G)
 
473
    weights = [ G.get_edge_data(G.edges()[i][0],G.edges()[i][1])['weight'] for i in range(len(G.edges())) ]
 
474
    nx.draw(G,pos,edge_color=weights,edge_cmap=cm.gray,edge_vmin=0., edge_vmax=1., width=4)
 
475
    plt.show(block=False) # display
 
476
 
 
477
    # Using all idxs \leq maxordF
 
478
    ls_tmp = [ range(maxordF+1) for i in range(d) ]
 
479
    idxs = list(itertools.product( *ls_tmp ))
 
480
 
 
481
    coeffs = np.array([ np.exp(- np.dot( np.array(idx), np.dot(Sigma,np.array(idx))) ) for idx in idxs ])
 
482
 
 
483
    names = ['Exp f4']
 
484
 
 
485
    P = S1D.Poly1D(S1D.JACOBI,[0.,0.])
 
486
    (x,w) = P.Quadrature(size1D)
 
487
    VV = [ P.GradVandermonde1D(x,size1D,0,norm=True)] * d
 
488
    Gammas = [ np.diag(np.dot(VV[i],VV[i].T)) for i in range(d) ]
 
489
 
 
490
    def f(X,params):        
 
491
        V = np.array([1.])
 
492
        for i in range(d):
 
493
            (idx,) = np.where(x == X[i])
 
494
            V = np.kron( V, params['VV'][i][idx,:] )
 
495
        normFact = np.sum( [np.dot( w, params['VV'][i][:,0]) for i in range(d) ] )
 
496
        out = np.dot( params['coeffs'], V.flatten() ) / normFact
 
497
        return out
 
498
 
 
499
    # Used to compute MC
 
500
    def fMC(X,params): 
 
501
        if isinstance(X,np.ndarray):
 
502
            out = np.zeros(X.shape[0])
 
503
        else:
 
504
            out = 0.
 
505
        for idxs,coeff in zip(params['idxs'],params['coeffs']):
 
506
            if coeff > 1e-13:
 
507
                tmp = 1.
 
508
                for i,ii in enumerate(idxs):
 
509
                    tmp *= params['VV'][i][:,ii] 
 
510
 
 
511
                out += coeff * tmp / ( d * np.dot( w, VV[0][:,0]))
 
512
        return out
 
513
 
 
514
    X = [x for i in range(d)]
 
515
    W = [w for i in range(d)]
 
516
    V_func = [ P.GradVandermonde1D(x,maxordF,0,norm=True)] * d
 
517
    params = { 'idxs': list(idxs),
 
518
               'coeffs': coeffs,
 
519
               'VV': V_func,
 
520
               'Gammas': Gammas}
 
521
    TW = DT.TensorWrapper( f, X, params )
 
522
    
 
523
    TTapprox = DT.QTTvec(TW)
 
524
    TTapprox.build(eps=eps,method='ttdmrg')
 
525
 
 
526
    print "TTdmrg: grid size %d" % TW.get_size()
 
527
    print "TTdmrg: function evaluations %d" % TW.get_fill_level()
 
528
 
 
529
    # Compute Fourier coefficients TT
 
530
    Vs = [P.GradVandermonde1D(X[i],size1D,0,norm=True)] * d
 
531
    TT_four = TTapprox.project(Vs,W)
 
532
 
 
533
    # Check approximation error using MC
 
534
    MCestVarLimit = 1e-1
 
535
    MCestMinIter = 100
 
536
    MCestMaxIter = 1e4
 
537
    MCstep = 10000
 
538
    var = 1.
 
539
    dist = stats.uniform()
 
540
    DIST = UQ.MultiDimDistribution([dist] * d)
 
541
    intf = []
 
542
    values = []
 
543
    while (len(values) < MCestMinIter or var > mean**2. * MCestVarLimit) and len(values) < MCestMaxIter:
 
544
        # Monte Carlo
 
545
        xx = DIST.rvs(MCstep)
 
546
 
 
547
        VsI = None
 
548
        VsI = [P.GradVandermonde1D(xx[:,i]*2.-1.,size1D,0,norm=True) for i in range(d)]
 
549
        TTval = TT_four.interpolate(VsI)
 
550
 
 
551
        TTvals = [ TTval[tuple([i]*d)] for i in range(MCstep) ]
 
552
 
 
553
        VsI_func = [P.GradVandermonde1D(xx[:,i]*2.-1.,maxordF,0,norm=True) for i in range(d)]
 
554
        paramsMC = { 'idxs': list(idxs),
 
555
                     'coeffs': coeffs,
 
556
                     'VV': VsI_func,
 
557
                     'Gammas': Gammas}
 
558
        fval = fMC(xx,paramsMC)
 
559
 
 
560
        intf.extend( list(fval**2.) )
 
561
        values.extend( [(fval[i]-TTvals[i])**2. for i in range(MCstep)])
 
562
        mean = vol * np.mean(values)
 
563
        var = vol**2. * np.var(values) / len(values)
 
564
 
 
565
        sys.stdout.write("L2err estim. iter: %d Var: %e VarLim: %e \r" % (len(values) , var, mean**2. * MCestVarLimit))
 
566
        sys.stdout.flush()
 
567
 
 
568
    sys.stdout.write("\n")
 
569
    sys.stdout.flush()
 
570
    L2err = np.sqrt(mean/np.mean(intf))
 
571
    print "TTdmrg: L2 err TTapprox: %e" % L2err
 
572
 
 
573
    # Error of the TT Fourier coeff., with respect to their exact value
 
574
    idx_slice = []
 
575
    for i in range(d):
 
576
        if i in nzdim:
 
577
            idx_slice.append( slice(None,None,None) )
 
578
        else:
 
579
            idx_slice.append( 0 )
 
580
    idx_slice = tuple(idx_slice)
 
581
    four_approx = TT_four[ idx_slice ]
 
582
    print "TTdmrg: L2 norm error TT_four: %e" % npla.norm( four_approx.flatten() - np.asarray(coeffs) ,2)
 
583
    print "TTdmrg: Inf norm error TT_four: %e" % npla.norm( four_approx.flatten() - np.asarray(coeffs) , np.inf)
 
584
 
 
585
    if size1D**d < 1e4:
 
586
        # Construct full tensor
 
587
        full_tens = TW.copy()[ tuple([slice(None,None,None) for i in range(d)]) ]
 
588
 
 
589
        print "TTdmrg: Frobenius err.: %e" % (npla.norm( (TTapprox.to_tensor()-full_tens).flatten()) / npla.norm(full_tens.flatten()))
 
590
        
 
591
        TTapproxSVD = DT.TTvec(full_tens)
 
592
 
 
593
        # Compute Fourier coefficients TT
 
594
        Vs = [P.GradVandermonde1D(X[i],size1D,0,norm=True)] * d
 
595
        TT_fourSVD = TTapproxSVD.project(Vs,W)
 
596
 
 
597
        four_tensSVD = np.zeros( tuple( [size1D+1 for i in range(len(nzdim))] ) )
 
598
        import itertools
 
599
        ls_tmp = [ range(size1D+1) for i in range(len(nzdim)) ]
 
600
        idx_tmp = itertools.product( *ls_tmp )
 
601
        for ii  in idx_tmp:
 
602
            ii_tt = [0]*d
 
603
            for jj, tti in enumerate(nzdim): ii_tt[tti] = ii[jj]
 
604
            ii_tt = tuple(ii_tt)
 
605
            four_tensSVD[ii] = TT_fourSVD[ii_tt]
 
606
        
 
607
        print "TT-SVD: ranks: %s" % str( TTapproxSVD.ranks() )
 
608
        print "TT-SVD: Frobenius norm TT_four:  %e" % mla.norm(TT_fourSVD,'fro')
 
609
        print "TT-SVD: Frobenius norm sub_tens: %e" % npla.norm(four_tensSVD.flatten())
 
610
 
 
611
    VMAX = 0.
 
612
    VMIN = np.inf
 
613
    if size1D**d < 1e4 and IS_PLOTTING and len(nzdim) == 2:
 
614
        # Plot 2D Fourier Coeff
 
615
        TTsvd_fourier_abs = np.maximum(np.abs(four_tensSVD),1e-20*np.ones(four_tens.shape))
 
616
        VMAX = np.max(np.log10(TTsvd_fourier_abs))
 
617
        VMIN = np.min(np.log10(TTsvd_fourier_abs))
 
618
 
 
619
    if IS_PLOTTING and len(nzdim) == 2:
 
620
        # Plot 2D Fourier Coeff
 
621
        TTdmrg_fourier_abs = np.maximum(np.abs(four_tens),1e-20*np.ones(four_tens.shape))
 
622
        VMAX = max(VMAX, np.max(np.log10(TTdmrg_fourier_abs)))
 
623
        VMIN = min(VMIN, np.min(np.log10(TTdmrg_fourier_abs)))
 
624
 
 
625
    if IS_PLOTTING and d == 2:
 
626
        # Plot 2D Fourier Coeff from TensorToolbox.product
 
627
        V2D = np.kron(Vs[0],Vs[1])
 
628
        WW = np.kron(W[0],W[1])
 
629
        fhat = np.dot(V2D.T,WW*TW.copy()[:,:].flatten()).reshape([size1D+1 for s in range(d)])
 
630
        fhat_abs = np.maximum(np.abs(fhat),1e-20*np.ones(fhat.shape))
 
631
        VMAX = max(VMAX, np.max(np.log10(fhat_abs)))
 
632
        VMIN = min(VMIN, np.min(np.log10(fhat_abs)))
 
633
 
 
634
    if size1D**d < 1e4 and IS_PLOTTING and len(nzdim) == 2:
 
635
        # Plot 2D Fourier Coeff
 
636
        fig = plt.figure(figsize=fsize)
 
637
        plt.imshow(np.log10(TTsvd_fourier_abs), interpolation='none', vmin = VMIN, vmax = VMAX, origin='lower')
 
638
        plt.colorbar()
 
639
        plt.title('TT-SVD Fourier - %s' % names[COEFF])
 
640
 
 
641
    if IS_PLOTTING and len(nzdim) == 2:
 
642
        # Plot 2D Fourier Coeff
 
643
        fig = plt.figure(figsize=fsize)
 
644
        plt.imshow(np.log10(TTdmrg_fourier_abs), interpolation='none', vmin = VMIN, vmax = VMAX, origin='lower')
 
645
        plt.colorbar()
 
646
        plt.title('TT-dmrg Fourier - %s' % names[COEFF])
 
647
 
 
648
    if IS_PLOTTING and d == 2:
 
649
        # Plot 2D Fourier Coeff from TensorToolbox.product
 
650
        fig = plt.figure(figsize=fsize)
 
651
        plt.imshow(np.log10(fhat_abs), interpolation='none', vmin = VMIN, vmax = VMAX, origin='lower')
 
652
        plt.colorbar()
 
653
        plt.title('Full Fourier - %s' % names[COEFF])
 
654
 
 
655
    if size1D**d < 1e4 and IS_PLOTTING and len(nzdim) == 2:
 
656
        # TT-svd - TT-dmrg
 
657
        fig = plt.figure(figsize=fsize)
 
658
        plt.imshow(np.log10(np.abs(TTsvd_fourier_abs - TTdmrg_fourier_abs)), interpolation='none', origin='lower')
 
659
        plt.colorbar()
 
660
        plt.title('(TT-svd - TT-dmrg) Fourier - %s' % names[COEFF])
 
661
 
 
662
    if IS_PLOTTING and d == 2:    
 
663
        # Full - TT-dmrg
 
664
        fig = plt.figure(figsize=fsize)
 
665
        plt.imshow(np.log10(np.abs(fhat_abs - TTdmrg_fourier_abs)), interpolation='none', origin='lower')
 
666
        plt.colorbar()
 
667
        plt.title('(Full - TT-dmrg) Fourier - %s' % names[COEFF])
 
668
        
 
669
plt.show(block=False)