1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
|
/*1:*/
#line 6 "./normal_moments.cweb"
#include "normal_moments.h"
#include "permutation.h"
#include "kron_prod.h"
#include "tl_static.h"
/*2:*/
#line 18 "./normal_moments.cweb"
UNormalMoments::UNormalMoments(int maxdim,const TwoDMatrix&v)
:TensorContainer<URSingleTensor> (1)
{
if(maxdim>=2)
generateMoments(maxdim,v);
}
/*:2*/
#line 12 "./normal_moments.cweb"
;
/*3:*/
#line 35 "./normal_moments.cweb"
void UNormalMoments::generateMoments(int maxdim,const TwoDMatrix&v)
{
TL_RAISE_IF(v.nrows()!=v.ncols(),
"Variance-covariance matrix is not square in UNormalMoments constructor");
int nv= v.nrows();
URSingleTensor*mom2= new URSingleTensor(nv,2);
mom2->getData()= v.getData();
insert(mom2);
URSingleTensor*kronv= new URSingleTensor(nv,2);
kronv->getData()= v.getData();
for(int d= 4;d<=maxdim;d+= 2){
URSingleTensor*newkronv= new URSingleTensor(nv,d);
KronProd::kronMult(ConstVector(v.getData()),
ConstVector(kronv->getData()),
newkronv->getData());
delete kronv;
kronv= newkronv;
URSingleTensor*mom= new URSingleTensor(nv,d);
/*4:*/
#line 70 "./normal_moments.cweb"
mom->zeros();
const EquivalenceSet eset= ebundle.get(d);
for(EquivalenceSet::const_iterator cit= eset.begin();
cit!=eset.end();cit++){
if(selectEquiv(*cit)){
Permutation per(*cit);
per.inverse();
for(Tensor::index it= kronv->begin();it!=kronv->end();++it){
IntSequence ind(kronv->dimen());
per.apply(it.getCoor(),ind);
Tensor::index it2(mom,ind);
mom->get(*it2,0)+= kronv->get(*it,0);
}
}
}
/*:4*/
#line 55 "./normal_moments.cweb"
;
insert(mom);
}
delete kronv;
}
/*:3*/
#line 13 "./normal_moments.cweb"
;
/*5:*/
#line 88 "./normal_moments.cweb"
bool UNormalMoments::selectEquiv(const Equivalence&e)
{
if(2*e.numClasses()!=e.getN())
return false;
for(Equivalence::const_seqit si= e.begin();
si!=e.end();++si){
if((*si).length()!=2)
return false;
}
return true;
}
/*:5*/
#line 14 "./normal_moments.cweb"
;
/*6:*/
#line 103 "./normal_moments.cweb"
FNormalMoments::FNormalMoments(const UNormalMoments&moms)
:TensorContainer<FRSingleTensor> (1)
{
for(UNormalMoments::const_iterator it= moms.begin();
it!=moms.end();++it){
FRSingleTensor*fm= new FRSingleTensor(*((*it).second));
insert(fm);
}
}
/*:6*/
#line 15 "./normal_moments.cweb"
;
/*:1*/
|
