« back to all changes in this revision
Viewing changes to FIAT/polynomial.py
- browse files at revision 2
- compare with another revision
- download diff
- download tarball
- view history from revision 2
- Committer: kirby
- Date: 2005-05-31 15:03:57 UTC
- Revision ID: devnull@localhost-20050531150357-n7ka5987yo6ylpaz
[project @ 2005-05-31 15:03:57 by kirby]
Initial revision
Initial revision
- files added:
- FIAT
- doc
added
removed
FIAT/polynomial.py
1
# Written by Robert C. Kirby
2
# Copyright 2005 by The University of Chicago
3
# Distributed under the LGPL license
4
# This work is partially supported by the US Department of Energy
5
# under award number DE-FG02-04ER25650
6
7
# last modified 2 May 2005
8
9
import expansions, shapes, Numeric, curry, points, LinearAlgebra
10
from LinearAlgebra import singular_value_decomposition as svd
11
class PolynomialBase( object ):
12
"""Defines an object that can tabulate the orthonormal polynomials
13
of a particular degree on a particular shape """
14
def __init__( self , shape , n ):
15
"""shape is a shape code defined in shapes.py
16
n is the polynomial degree"""
17
self.bs = expansions.make_expansion( shape , n )
18
self.shape,self.n = shape,n
19
d = shapes.dimension( shape )
20
if n == 0:
21
self.dmats = [ Numeric.array( [ [ 0.0 ] ] , "d" ) ] * d
22
else:
23
dtildes = [ Numeric.zeros( (len(self.bs),len(self.bs)) , "d" ) \
24
for i in range(0,d) ]
25
pts = points.make_points( shape , d , 0 , n + d + 1 )
26
v = Numeric.transpose( expansions.tabulators[shape](n,pts) )
27
vinv = LinearAlgebra.inverse( v )
28
dtildes = [Numeric.transpose(a) \
29
for a in expansions.deriv_tabulators[shape](n,pts)]
30
self.dmats = [ Numeric.dot( vinv , dtilde ) for dtilde in dtildes ]
31
return
32
33
def eval_all( self , x ):
34
"""Returns the array A[i] = phi_i(x)."""
35
return self.tabulate( Numeric.array([x]))[:,0]
36
37
def tabulate( self , xs ):
38
"""xs is an iterable object of points.
39
returns an array A[i,j] where i runs over the basis functions
40
and j runs over the elements of xs."""
41
return expansions.tabulators[self.shape](self.n,Numeric.array(xs))
42
43
def degree( self ):
44
return self.n
45
46
def spatial_dimension( self ):
47
return shapes.dimension( self.shape )
48
49
def domain_shape( self ):
50
return self.shape
51
52
def __len__( self ):
53
"""Returns the number of items in the set."""
54
return len( self.bs )
55
56
def __getitem__( self , i ):
57
return self.bs[i]
58
59
60
class AbstractPolynomialSet( object ):
61
"""A PolynomialSet is a collection of polynomials defined as
62
linear combinations of a particular base set of polynomials.
63
In this code, these are the orthonormal polynomials defined
64
on simplices in one, two, or three spatial dimensions.
65
PolynomialSets may contain array-valued elements."""
66
def __init__( self , base , coeffs ):
67
"""Base is an instance of PolynomialBase.
68
coeffs[i][j] is the coefficient (possibly array-valued)
69
of the j:th orthonormal function in the i:th member of
70
the polynomial set. If rank("""
71
self.base, self.coeffs = base, coeffs
72
return
73
def __getitem__( self , i ):
74
"""If i is an integer, returns the i:th member of the set
75
as the appropriate subtype of AbstractPolynomial. If
76
i is a slice, returns the the implied subset of polynomials
77
as the appropriate subtype of AbstractPolynomialSet."""
78
if type(i) == type(1): # single item
79
return Polynomial( self.base , self.coeffs[i] )
80
elif type(i) == type(slice(1)):
81
return PolynomialSet( self.base , self.coeffs[i] )
82
def eval_all( self , x ):
83
"""Returns the array A[i] of the i:th member of the set
84
at point x."""
85
pass
86
def tabulate( self , xs ):
87
"""Returns the array A[i][j] with the i:th member of the
88
set evaluated at the j:th member of xs."""
89
pass
90
def degree( self ):
91
"""Returns the polynomial degree of the space. If the polynomial
92
set lies between two degrees, such as the Raviart-Thomas space, then
93
the degree of the smallest space of complete degree containing the
94
set is returned. For example, the degree of RT0 is 1 since some
95
linear functions are required to represent the basis."""
96
return self.base.degree()
97
def rank( self ):
98
"""Returns the tensor rank of the range of the functions (0
99
for scalar, two for vector, etc"""
100
return Numeric.rank( self.coeffs ) - 2
101
def tensor_dim( self ):
102
if self.rank() == 0:
103
return tuple()
104
else:
105
return self.coeffs.shape[1:-1]
106
def tensor_shape( self ):
107
pass
108
def spatial_dimension( self ):
109
return shapes.dimension( self.base.shape )
110
def domain_shape( self ):
111
"""Returns the code for the element shape of the domain. This
112
is the code given in module shapes.py"""
113
return self.base.shape
114
def __len__( self ):
115
"""Returns the number of items in the set."""
116
return len( self.coeffs )
117
def take( self , items ):
118
"""Extracts the subset of polynomials given by indices stored
119
in iterable items."""
120
return PolynomialSet( self.base , Numeric.take( self.coeffs , items ) )
121
122
# ScalarPolynomialSet can be made with either
123
# -- ScalarPolynomialSet OR
124
# -- PolynomialBase
125
# along with a rank 2 array, where each row is the
126
# coefficients of the members of the base for
127
# a member of the new set
128
class ScalarPolynomialSet( AbstractPolynomialSet ):
129
def __init__( self , base , coeffs ):
130
if Numeric.rank( coeffs ) != 2:
131
raise RuntimeError, \
132
"Illegal coeff matrix: ScalarPolynomialSet"
133
if isinstance( base, PolynomialBase ):
134
AbstractPolynomialSet.__init__( self , base , coeffs )
135
elif isinstance( base , ScalarPolynomialSet ):
136
new_base = base.base
137
new_coeffs = Numeric.dot( coeffs , base.coeffs )
138
AbstractPolynomialSet.__init__( self , new_base , new_coeffs )
139
return
140
def eval_all( self , x ):
141
"""Returns the array A[i] = psi_i(x)."""
142
bvals = self.base.eval_all( x )
143
return Numeric.dot( self.coeffs , bvals )
144
145
def tabulate( self , xs ):
146
"""Returns the array A[i][j] with i running members of the set
147
and j running over the members of xs."""
148
bvals = self.base.tabulate( xs )
149
return Numeric.dot( self.coeffs , bvals )
150
151
def deriv_all( self , i ):
152
"""Returns the PolynomialSet containing the partial derivative
153
in the i:th direction of each component."""
154
D = self.base.dmats[i]
155
new_coeffs = Numeric.dot( self.coeffs , Numeric.transpose(D) )
156
return ScalarPolynomialSet( self.base , new_coeffs )
157
158
def multi_deriv_all( self , alpha ):
159
"""Returns the PolynomialSet containing the alpha partial
160
derivative of everything. alpha is a multiindex of the same
161
size as the spatial dimension."""
162
U = self
163
for c in range(len(alpha)):
164
for i in range(alpha[c]):
165
U = U.deriv_all( c )
166
return U
167
168
def tabulate_jet( self , order , xs ):
169
"""Computes all partial derivatives of the members of the set
170
up to order. Returns an array of dictionaries
171
a[i][mi] where i is the order of partial differentiation
172
and mi is a multiindex with |mi| = i. The value of
173
a[i][mi] is an array A[i][j] containing the appropriate derivative
174
of the i:th member of the set at the j:th member of xs."""
175
alphas = [ mis( shapes.dimension( self.base.shape ) , i ) \
176
for i in range(order+1) ]
177
a = [ None ] * len(alphas)
178
for i in range(len(alphas)):
179
a[i] = {}
180
for alpha in alphas[i]:
181
a[i][alpha] = self.multi_deriv_all( alpha ).tabulate( xs )
182
return a
183
184
def tensor_shape( self ):
185
return (1,)
186
187
# A VectorPolynomialSet may be made with either
188
# -- PolynomialBase and rank 3 array of coefficients OR
189
# -- VectorPolynomialSet and rank 2 array
190
191
# coeffs for a VectorPolynomialSet is an array
192
# C[i,j,k] where i runs over the members of the set,
193
# j runs over the components of the vectors, and
194
# k runs over the members of the base set.
195
196
class VectorPolynomialSet( AbstractPolynomialSet ):
197
"""class modeling sets of vector-valued polynomials."""
198
def __init__( self , base , coeffs ):
199
if isinstance( base , PolynomialBase ):
200
if Numeric.rank( coeffs ) != 3:
201
raise RuntimeError, \
202
"Illegal coeff matrix: VectorPolynomialSet"
203
AbstractPolynomialSet.__init__( self , base , coeffs )
204
pass
205
elif isinstance( base , VectorPolynomialSet ):
206
if Numeric.rank( coeffs ) != 2:
207
raise RuntimeError, \
208
"Illegal coeff matrix: VectorPolynomialSet"
209
new_base_shape = (base.coeffs.shape[0], \
210
reduce(lambda a,b:a*b, \
211
base.coeffs.shape[1:] ) )
212
base_coeffs_reshaped = Numeric.reshape( base.coeffs , \
213
new_base_shape )
214
new_coeffs_shape = tuple( [coeffs.shape[0]] \
215
+ list(base.coeffs.shape[1:]) )
216
new_coeffs_flat = Numeric.dot( coeffs , base_coeffs_reshaped )
217
new_coeffs = Numeric.reshape( new_coeffs_flat , new_coeffs_shape )
218
AbstractPolynomialSet.__init__( self , base.base , new_coeffs )
219
return
220
else:
221
raise RuntimeError, "Illegal base set: VectorPolynomialSet"
222
223
def eval_all( self , x ):
224
"""Returns the array A[i,j] where i runs over the members of the
225
set and j runs over the components of each member."""
226
bvals = self.base.eval_all( x )
227
old_shape = self.coeffs.shape
228
flat_coeffs = Numeric.reshape( self.coeffs , \
229
(old_shape[0]*old_shape[1] , \
230
old_shape[2] ) )
231
flat_dot = Numeric.dot( flat_coeffs , bvals )
232
return Numeric.reshape( flat_dot , old_shape[:2] )
233
234
def tabulate( self , xs ):
235
"""xs is an iterable object of points.
236
returns an array A[i,j,k] where i runs over the members of the
237
set, j runs over the components of the vectors, and k runs
238
over the points."""
239
bvals = self.base.tabulate( xs )
240
old_shape = self.coeffs.shape
241
flat_coeffs = Numeric.reshape( self.coeffs , \
242
( old_shape[0]*old_shape[1] , \
243
old_shape[2] ) )
244
flat_dot = Numeric.dot( flat_coeffs , bvals )
245
unflat_dot = Numeric.reshape( flat_dot , \
246
( old_shape[0] , old_shape[1] , len(xs) ) )
247
return unflat_dot
248
249
def select_vector_component( self , i ):
250
"""Returns the ScalarPolynomialSetconsisting of the i:th
251
component of allthe vectors. It keeps zeros around. This
252
could change in later versions."""
253
return ScalarPolynomialSet( self.base , self.coeffs[:,i,:] )
254
255
def tensor_shape( self ):
256
return self.coeffs.shape[1:2]
257
258
def tabulate_jet( self , order , xs ):
259
return [ self.select_vector_component( i ).tabulate_jet( order , xs ) \
260
for i in range(self.tensor_shape()[0]) ]
261
262
263
# The code for TensorPolynomialSet will look just like VectorPolynomialSet
264
# except that we will flatten all of the coefficients to make them
265
# look like vectors instead of tensors.
266
267
class TensorPolynomialSet( AbstractPolynomialSet ):
268
pass
269
270
def PolynomialSet( base , coeffs ):
271
"""Factory function that takes some PolynomialBase or
272
AbstractPolynomialSet and a collection of coefficients and
273
either returns the appropriate kind of subclass of
274
AbstractPolynomialSet or else raises an exception."""
275
if len( coeffs.shape ) == 2:
276
if isinstance( base , PolynomialBase ):
277
return ScalarPolynomialSet( base , coeffs )
278
elif isinstance( base , VectorPolynomialSet ):
279
return VectorPolynomialSet( base , coeffs )
280
elif isinstance( base , ScalarPolynomialSet ):
281
return ScalarPolynomialSet( base , coeffs )
282
else:
283
raise RuntimeError, "???: PolynomialSet"
284
elif len( coeffs.shape ) == 3:
285
return VectorPolynomialSet( base , coeffs )
286
elif len( coeffs.shape ) > 3:
287
return TensorPolynomialSet( base , coeffs )
288
else:
289
raise RuntimeError, "Unknown error, PolynomialSet"
290
291
292
class AbstractPolynomial( object ):
293
"""Base class from which scalar- vector- and tensor- valued
294
polynomials are implemented. At least, all types of
295
polynomials must support:
296
-- evaluation via __call__().
297
All polynomials are represented as a linear combination of
298
some base set of polynomials."""
299
def __init__( self , base , dof ):
300
"""Constructor for AbstractPolynomial -- must be provided
301
with PolynomialBase instance and a Numeric.array of
302
coefficients."""
303
if not isinstance( base , PolynomialBase ):
304
raise RuntimeError, "Illegal base: AbstractPolynomial"
305
self.base , self.dof = base , dof
306
return
307
def __call__( self , x ): pass
308
def __getitem__( self , x ): pass
309
def degree( self ): return self.base.degree()
310
311
class ScalarPolynomial( AbstractPolynomial ):
312
"""class of scalar-valued polynomials supporting
313
evaluation and differentiation."""
314
def __init__( self , base , dof ):
315
if Numeric.rank( dof ) != 1:
316
raise RuntimeError, \
317
"This isn't a scalar polynomial."
318
AbstractPolynomial.__init__( self , base , dof )
319
return
320
def __call__( self , x ):
321
"""Evaluates the polynomial at point x"""
322
bvals = self.base.eval_all( x )
323
return Numeric.dot( self.dof , bvals )
324
def deriv( self , i ):
325
"""Computes the partial derivative in the i:th direction,
326
represented as a polynomial over the same base."""
327
b = self.base
328
D = b.dmats[i]
329
deriv_dof = Numeric.dot( D , self.dof )
330
return ScalarPolynomial( b , deriv_dof )
331
def __getitem__( self , i ):
332
if i != 0: raise RuntimeError, "Illegal indexing into ScalarPolynomial"
333
return self
334
335
class VectorPolynomial( AbstractPolynomial ):
336
"""class of vector-valued polynomials supporting evaluation
337
and component selection."""
338
def __init__( self , base , dof ):
339
if Numeric.rank( dof ) != 2 :
340
raise RuntimeError, "This isn't a vector polynomial."
341
AbstractPolynomial.__init__( self , base , dof )
342
return
343
def __call__( self , x ):
344
bvals = self.base.eval_all( x )
345
return Numeric.dot( self.dof , bvals )
346
def __getitem__( self , i ):
347
if type( i ) != type( 1 ):
348
raise RuntimeError, "Illegal input type."
349
return ScalarPolynomial( self.base , self.dof[i] )
350
351
class TensorPolynomial( AbstractPolynomial ):
352
pass
353
354
# factory function that determines whether to instantiate
355
# a scalar, vector, or general tensor polynomial
356
def Polynomial( base , dof ):
357
"""Returns a instance of the appropriate subclass of
358
AbstractPolynomial."""
359
if not isinstance( base , PolynomialBase ) and \
360
not isinstance( base , PolynomialSet ):
361
raise RuntimeError, "Illegal types, Polynomial"
362
if Numeric.rank( dof ) == 1:
363
return ScalarPolynomial( base , dof )
364
elif Numeric.rank( dof ) == 2:
365
return VectorPolynomial( base , dof )
366
elif Numeric.rank( dof ) > 2:
367
return TensorPolynomial( base , dof )
368
else:
369
raise RuntimeError, "Illegal shape dimensions"
370
371
def OrthogonalPolynomialSet( element_shape , degree ):
372
"""Returns a ScalarPolynomialSet that models the
373
orthormal basis functions on element_shape. This allows
374
us to arbitrarily differentiate the orthogonal polynomials
375
if we need to."""
376
b = PolynomialBase( element_shape , degree )
377
coeffs = Numeric.identity( shapes.poly_dims[element_shape](degree) , \
378
"d" )
379
return ScalarPolynomialSet( b , coeffs )
380
381
def OrthogonalPolynomialArraySet( element_shape , degree , nc = None ):
382
"""Returns a VectorPolynomialSet that models the orthnormal basis
383
for vector-valued functions with d components, where d is the spatial
384
dimension of element_shape"""
385
b = PolynomialBase( element_shape , degree )
386
space_dim = shapes.dims[ element_shape ]
387
if nc == None:
388
nc = space_dim
389
M = shapes.poly_dims[ element_shape ]( degree )
390
coeffs = Numeric.zeros( (nc * M , nc , M ) , "d" )
391
ident = Numeric.identity( M , "d" )
392
for i in range(nc):
393
coeffs[(i*M):(i+1)*M,i,:] = ident[:,:]
394
395
return PolynomialSet( b , coeffs )
396
397
class FiniteElement( object ):
398
"""class modeling the Ciarlet abstraction of a finite element"""
399
def __init__( self , Udual , U ):
400
self.Udual = Udual
401
v = outer_product( Udual.get_functional_set() , U )
402
vinv = LinearAlgebra.inverse( v )
403
self.U = PolynomialSet( U , Numeric.transpose(vinv) )
404
return
405
def domain_shape( self ): return self.U.domain_shape()
406
def function_space( self ): return self.U
407
def dual_basis( self ): return self.Udual
408
409
def ConstrainedPolynomialSet( fset ):
410
# takes a FunctionalSet object and constructs the PolynomialSet
411
# consisting of the intersection of the null spaces of its member
412
# functionals acting on its function set
413
tol = 1.e-12
414
L = outer_product( fset , fset.function_space() )
415
(U_L , Sig_L , Vt_L) = svd( L , 1 ) # full svd
416
417
# some of the constraint functionals may be redundant. I
418
# must check to see how many singular values there are, as this
419
# is what determines the rank and nullity of L
420
Sig_L_nonzero = Numeric.array( [ a for a in Sig_L if abs(a) > tol ] )
421
num_nonzero_svs = len( Sig_L_nonzero )
422
vcal = Vt_L[ num_nonzero_svs: ]
423
return PolynomialSet( fset.function_space() , vcal )
424
425
def outer_product(flist,plist):
426
# if the functions are vector-valued, I need to reshape the coefficient
427
# arrays so they are two-dimensional to do the dot product
428
shp = flist.mat.shape
429
if len(shp) > 2:
430
num_cols = reduce( lambda a,b:a*b , shp[1:] )
431
new_shp = (shp[0],num_cols)
432
A = Numeric.reshape( flist.mat,(flist.mat.shape[0],num_cols) )
433
B = Numeric.reshape( plist.coeffs,(plist.coeffs.shape[0],num_cols) )
434
else:
435
A = flist.mat
436
B = plist.coeffs
437
438
return Numeric.dot( A , Numeric.transpose( B ) )
439
440
def gradient( u ):
441
if not isinstance( u , ScalarPolynomial ):
442
raise RuntimeError, "Illegal input to gradient"
443
new_dof = Numeric.zeros( (u.base.spatial_dimension() , len(u.dof) ) , "d" )
444
for i in range(u.base.spatial_dimension()):
445
new_dof[i,:] = Numeric.dot( u.base.dmats[i] , \
446
u.dof )
447
return VectorPolynomial( u.base , new_dof )
448
449
def gradients( U ):
450
"""For U a ScalarPolynomialSet, computes the gradient of each member
451
of U, returning a VectorPolynomialSet."""
452
if not isinstance( U , ScalarPolynomialSet ):
453
raise RuntimeError, "Illegal input to gradients"
454
# new matrix is len(U) by spatial_dimension by length of poly base
455
new_dofs = Numeric.zeros( (len(U),U.spatial_dimension(),len(U.base)) , "d" )
456
for i in range(U.spatial_dimension()):
457
new_dofs[:,i,:] = Numeric.dot( U.coeffs , U.base.dmats[i] )
458
return VectorPolynomialSet( U.base , new_dofs )
459
460
def divergence( u ):
461
if not isinstance( u , VectorPolynomial ):
462
raise RuntimeError, "Illegal input to divergence"
463
new_dof = Numeric.zeros( len(u.base) , "d" )
464
for i in range(u.base.spatial_dimension()):
465
new_dof[:] += Numeric.dot( u.base.dmats[i] , \
466
u.dof[i] )
467
return ScalarPolynomial( u.base , new_dof )
468
469
def curl( u ):
470
if not isinstance( u , VectorPolynomial ):
471
raise RuntimeError, "Illegal input to curl"
472
if u.base.domain_shape() != shapes.TETRAHEDRON:
473
raise RuntimeError, "Illegal shape to curl"
474
new_dof = Numeric.zeros( (3,len(u.base)) , "d" )
475
new_dof[0] = u[2].deriv(1).dof - u[1].deriv(2).dof
476
new_dof[1] = u[0].deriv(2).dof - u[2].deriv(0).dof
477
new_dof[2] = u[1].deriv(0).dof - u[0].deriv(1).dof
478
return VectorPolynomial( u.base, new_dof )
479
480
# | i j k |
481
# | d0 d1 d2 |
482
# | v0 v1 v2 |
483
# -->
484
# (d1*v2-d2*v1,d2*v0-d0*v2,d0*v1-d1*v0)
485
# OR
486
# | 0 -d2 d1 | | v0 |
487
# | d2 0 -d0 | | v1 |
488
# | -d1 d0 0 | | v2 |
489
# SO...
490
# curl^t:
491
# | 0 d2t -d1t |
492
# | -d2t 0 d0t |
493
# | d1t -d0t 0 |
494
495
def curl_transpose( u ):
496
if not isinstance( u , VectorPolynomial ):
497
raise RuntimeError, "Illegal input to curl"
498
if u.base.domain_shape() != shapes.TETRAHEDRON:
499
raise RuntimeError, "Illegal shape to curl"
500
new_dof = Numeric.zeros( (3,len(u.base)) , "d" )
501
new_dof[0] = Numeric.dot( Numeric.transpose( u.base.dmats[2] ) , \
502
u[1].dof ) \
503
- Numeric.dot( Numeric.transpose( u.base.dmats[1] ) , \
504
u[2].dof )
505
new_dof[1] = Numeric.dot( Numeric.transpose( u.base.dmats[0] ) , \
506
u[2].dof ) \
507
- Numeric.dot( Numeric.transpose( u.base.dmats[2]) , \
508
u[0].dof )
509
new_dof[2] = Numeric.dot( Numeric.transpose( u.base.dmats[1] ) , \
510
u[0].dof ) \
511
- Numeric.dot( Numeric.transpose( u.base.dmats[0] ) , \
512
u[1].dof )
513
return VectorPolynomial( u.base , new_dof )
514
515
516
# this is used in tabulating jets.
517
def mis(m,n):
518
"""returns all m-tuples of nonnegative integers that sum up to n."""
519
if m==1:
520
return [(n,)]
521
elif n==0:
522
return [ tuple([0]*m) ]
523
else:
524
return [ tuple([n-i]+list(foo)) \
525
for i in range(n+1) \
526
for foo in mis(m-1,i) ]
527
528
# projection of some function onto a scalar polynomial set.
529
530
def projection( U , f , Q ):
531
f_at_qps = [ f(x) for x in Q.get_points() ]
532
phis_at_qps = U.tabulate( Q.get_points() )
533
return ScalarPolynomial( U.base , \
534
Numeric.array( [ sum( Q.get_weights() \
535
* f_at_qps \
536
* phi ) \
537
for phi in phis_at_qps ] ) )
538
539
# C --> u,sig,vt. first r columns of u span column range,
540
# so first r columns of v span row range
541
# so take first r rows of vt
542
543
def poly_set_union( U , V ):
544
"""Takes the union of two polynomial sets by appending their
545
coefficient tensors and computing the range of the resulting set
546
by the SVD."""
547
new_coeffs = Numeric.array( list( U.coeffs ) + list( V.coeffs ) )
548
func_shape = new_coeffs.shape[1:]
549
if len( func_shape ) == 1:
550
(u,sig,vt) = svd( new_coeffs )
551
num_sv = len( [ s for s in sig if abs( s ) > 1.e-10 ] )
552
return PolynomialSet( U.base , vt[:num_sv] )
553
else:
554
new_shape0 = new_coeffs.shape[0]
555
new_shape1 = reduce(lambda a,b:a*b,func_shape)
556
nc = Numeric.reshape( new_coeffs , (new_shape0,new_shape1) )
557
(u,sig,vt) = svd( nc , 1 )
558
num_sv = len( [ s for s in sig if abs( s ) > 1.e-10 ] )
559
coeffs = vt[:num_sv]
560
return PolynomialSet( U.base , \
561
Numeric.reshape( coeffs , \
562
tuple( [len(coeffs)] \
563
+ list( func_shape ) ) ) )
Loggerhead is a web-based interface for Breezy
Version: 2.0.1