source: tol/tol_pkg/MulVarFunAprox/DesignC2.Spheric.tol @ 77c938e

1941-docker-image-to-testdevel
Last change on this file since 77c938e was 77c938e, checked in by Víctor de Buen <vdebuen@…>, 17 years ago

Using efficient grouping functions

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@2048 be4846ae-f051-483a-831f-b64e7d08527f
  • Property mode set to 100644
File size: 8.1 KB
Line 
1/* -*- mode: C++ -*- */
2//////////////////////////////////////////////////////////////////////////////
3// FILE   : DesignC2.tol
4// PURPOSE: Defines Class @DesignC2
5//////////////////////////////////////////////////////////////////////////////
6
7
8//////////////////////////////////////////////////////////////////////////////
9Text _.autodoc.member.DesignC2.Spheric =
10"";
11Class @DesignC2.Spheric : @DesignC2
12//////////////////////////////////////////////////////////////////////////////
13{
14
15//////////////////////////////////////////////////////////////////////////////
16//Configuration parameters
17//////////////////////////////////////////////////////////////////////////////
18  Real _.radius.margin;
19
20//////////////////////////////////////////////////////////////////////////////
21//Auxiliar members
22//////////////////////////////////////////////////////////////////////////////
23  VMatrix _.R =Rand(0,0,0,0);
24  MatQuery::@ANN.KDTree _.ann =
25    MatQuery::@ANN.KDTree::Default(0);
26  MatQuery::@KNN.Neighbourhood _.ngbh =
27    MatQuery::@KNN.Neighbourhood::Default(0);
28
29  //////////////////////////////////////////////////////////////////////////////
30  //Methods
31  //////////////////////////////////////////////////////////////////////////////
32
33  //////////////////////////////////////////////////////////////////////////////
34  Real set.reference.points(VMatrix Y)
35  //////////////////////////////////////////////////////////////////////////////
36  {
37    VMatrix _.Y := Y;
38    Real _.H := VRows(_.Y);
39    Real _.ann::set.data(VMat2Mat(_.Y));
40    Real _.ann::update.search(_.ann::_.X,_.X.maxNeighbour,_.ngbh);
41    VMatrix _.R :=
42      Mat2VMat(Sqrt(SubCol(_.ngbh::_.distance2,[[_.X.maxNeighbour]])))*
43      (1+Min(0,_.radius.margin));
44    True
45  };
46
47  ////////////////////////////////////////////////////////////////////////////
48    VMatrix match.evaluable(VMatrix X)
49  ////////////////////////////////////////////////////////////////////////////
50  {
51    MatQuery::@KNN.Neighbourhood ngbh =
52      _.ann::search(VMat2Mat(X),1);   
53    VMatrix dist = Mat2VMat(Sqrt(ngbh::_.distance2));
54    VMatrix X.h = Mat2VMat(ngbh::_.neighbour);
55    LE(dist,get.radius(X.h))$*
56    LE(dist,Constant(VRows(dist),1,stress.maxDistance))
57  };
58
59  //////////////////////////////////////////////////////////////////////////////
60  VMatrix closest.reference.index(VMatrix X)
61  //////////////////////////////////////////////////////////////////////////////
62  {
63    MatQuery::@KNN.Neighbourhood ngbh =
64      _.ann::search(VMat2Mat(X),_.X.maxNeighbour);
65    Mat2VMat(ngbh::_.neighbour)
66  };
67
68  //////////////////////////////////////////////////////////////////////////////
69  VMatrix get.radius(VMatrix Y.idx)
70  //////////////////////////////////////////////////////////////////////////////
71  {
72    SubRow(_.R,VCol2Set(Y.idx))
73  };
74
75  //////////////////////////////////////////////////////////////////////////////
76    VMatrix build.obs(VMatrix X)
77  //////////////////////////////////////////////////////////////////////////////
78  {
79    FunRn2R::@BaseC2.Eval.Function eval =
80       FunRn2R::@BaseC2.Eval.Function::New(get.base(0));
81    MatQuery::@KNN.Neighbourhood ngbh =
82      _.ann::search(VMat2Mat(X),_.X.maxNeighbour);
83    Real m = VRows(X);
84    VMatrix dist = Mat2VMat(Sqrt(ngbh::_.distance2));
85    VMatrix X.h = Mat2VMat(ngbh::_.neighbour);
86    VMatrix X.idx = DifEq(1/(1-B),Constant(m,1,1));
87    VMatrix triplet.all = Group("ConcatRows", For(1,_.X.maxNeighbour,VMatrix(Real k)
88    {
89      VMatrix X.h.k = SubCol(X.h,[[k]]);
90      VMatrix dist.k = SubCol(dist,[[k]]);
91      VMatrix R = get.radius(X.h.k);
92      VMatrix tao = dist.k $/ R;
93      VMatrix lambda = (-tao^2+1)^4;
94      VMatrix match = LT(dist.k,R);
95      SubRow(X.idx | X.h.k | dist.k | tao| lambda, MatQuery::SelectMatch(match))
96    }));
97    VMatrix x.idx = SubCol(triplet.all,[[1]]);
98    VMatrix y.idx = SubCol(triplet.all,[[2]]);
99    VMatrix xy.idx = x.idx | y.idx;
100    VMatrix lambda.col = SubCol(triplet.all,[[5]]);
101    VMatrix triplet.lambda = xy.idx | lambda.col;
102    VMatrix lambda = Triplet(VMat2Mat(triplet.lambda),m,_.H);
103    VMatrix sum.lambda = lambda*Constant(_.H,1,1);
104    VMatrix sum.lambda.inv = Convert(sum.lambda,"Blas.R.Dense")^-1;
105    VMatrix weight = Eye(m,m,0,sum.lambda.inv)*lambda;
106    VMatrix dif =
107      SubRow(  X,VCol2Set(SubCol(triplet.all,[[1]]))) -
108      SubRow(_.Y,VCol2Set(SubCol(triplet.all,[[2]]))) ;
109    VMatrix P.row = Group("ConcatRows", For(1,_.K,VMatrix(Real k)
110    {
111      x.idx
112    }));
113    VMatrix P.y = Group("ConcatRows", For(1,_.K,VMatrix(Real k)
114    {
115      y.idx
116    }));
117    VMatrix YnCum = DifEq(1/(1-B),_.Yn);
118    VMatrix YnCumB = Sub(SetCol([[0]])<<YnCum,1,1,_.H,1);
119    VMatrix P.YnCum  = SubRow(YnCum, VCol2Set(P.y));
120    VMatrix P.YnCumB = SubRow(YnCumB,VCol2Set(P.y));
121    VMatrix P.k = Group("ConcatRows", For(1,_.K,VMatrix(Real k)
122    {
123      y.idx*0+k
124    }));
125    VMatrix P.col = P.YnCumB + P.k;
126    VMatrix P.match = LE(P.col,P.YnCum);
127    Real eval::upgrade.cache(dif);
128    VMatrix P.ev = Group("ConcatRows", For(1,_.K,VMatrix(Real k)
129    {
130      lambda.col $* eval::evaluate.v(k,dif)
131    }));
132    Matrix P.triplet = VMat2Mat(
133      SubRow(P.row | P.col | P.ev, MatQuery::SelectMatch(P.match)));
134    VMatrix Px = Pack(Convert(Triplet(P.triplet, m, _.N),"Cholmod.R.Sparse"));
135    VMatrix Qx = Eye(m,m,0,sum.lambda.inv)*Px
136  };
137
138
139  ///////////////////////////////////////////////////////////////////////////
140    VMatrix generate.internal.points(Real m)
141  ///////////////////////////////////////////////////////////////////////////
142  {
143    Text _MID = "[@DesignC2.Monomial.Spheric.MeshFree::generate.internal.points] ";
144    WriteLn("TRACE "+_MID+"1 m="<<m);
145    VMatrix h = Round(Rand(m,1,0.5,_.H+0.499999999));
146    WriteLn("TRACE "+_MID+"2 h:\n"<<h);
147    VMatrix y = SubRow(_.Y,VCol2Set(h));
148    WriteLn("TRACE "+_MID+"3 y:\n"<<y);
149    VMatrix v = Gaussian(m,_.n,0,1);
150    WriteLn("TRACE "+_MID+"4 v:\n"<<v);
151    VMatrix u = v*(VMatDat(MtMSqr(v),1,1)^-0.5);
152    WriteLn("TRACE "+_MID+"5 u:\n"<<u);
153    VMatrix rnd = get.radius(h)$*Rand(m,1,0,1);
154    WriteLn("TRACE "+_MID+"5 rnd:\n"<<rnd);
155    VMatrix r = Eye(m,m,0,rnd);
156    WriteLn("TRACE "+_MID+"6 r:\n"<<r);
157    y + r * u
158  };
159
160  ///////////////////////////////////////////////////////////////////////////
161  Set get.X.weight(Set X.dif, Set X.h)
162  ///////////////////////////////////////////////////////////////////////////
163  {
164    Real m = VRows(X.dif[1]);
165    Real numNeighbour = Card(X.h);
166    Set X.tao = For(1,numNeighbour,VMatrix(Real i)
167    {
168      VMatrix X.norm2 = RPow(X.dif[i],2)*Rand(_.n,1,1,1);
169      RPow(X.norm2,0.5)$/get.radius(X.h[i])
170    });
171    Set X.tao2 = For(1,numNeighbour,VMatrix(Real i)
172    {
173      RPow(X.tao[i],2)
174    });
175    Set X.match = For(1,numNeighbour,VMatrix(Real i)
176    {
177      LE(X.tao[i],Rand(m,1,1,1))
178    });
179    Set X.lambda = For(1,numNeighbour,VMatrix(Real i)
180    {
181      RPow(X.match[i] $* (Rand(m,1,1,1)-X.tao2[i]),4)
182    });
183    VMatrix X.lambda.sum = SetSum(X.lambda);
184   
185    For(1,numNeighbour,VMatrix(Real i)
186    {
187      X.lambda[i] $/ X.lambda.sum
188    })
189  };
190  ///////////////////////////////////////////////////////////////////////////
191  Set get.X.weight.partial.derivative.1(Set X.dif, Set X.h, Real s)
192  ///////////////////////////////////////////////////////////////////////////
193  {
194    Real m = VRows(X.dif[1]);
195    Real numNeighbour = Card(X.h);
196    Set X.tao = For(1,numNeighbour,VMatrix(Real i)
197    {
198      VMatrix X.norm2 = RPow(X.dif[i],2)*Rand(_.n,1,1,1);
199      RPow(X.norm2,0.5)$/get.radius(X.h[i])
200    });
201    Set X.tao2 = For(1,numNeighbour,VMatrix(Real i)
202    {
203      RPow(X.tao[i],2)
204    });
205    Set X.match = For(1,numNeighbour,VMatrix(Real i)
206    {
207      LE(X.tao[i],Rand(m,1,1,1))
208    });
209    Set X.lambda = For(1,numNeighbour,VMatrix(Real i)
210    {
211      RPow(X.match[i] $* (Rand(m,1,1,1)-X.tao2[i]),4)
212    });
213    Set X.lambda.der = For(1,numNeighbour,VMatrix(Real i)
214    {
215      RPow(X.match[i] $*
216           (Rand(m,1,1,1)-X.tao2[i]),3) $*
217           (SubCol(X.dif[i],[[s]])*8) $/ get.radius(X.h[i])
218    });
219    VMatrix X.lambda.sum = SetSum(X.lambda);
220    VMatrix X.lambda.sum.2 = RPow(X.lambda.sum,2);
221    VMatrix X.lambda.der.sum = SetSum(X.lambda.der);
222    For(1,numNeighbour,VMatrix(Real i)
223    {
224      (X.lambda.der[i]$*X.lambda.sum-X.lambda[i]$*X.lambda.der.sum)
225      $/ X.lambda.sum.2
226    })
227  }
228
229};
Note: See TracBrowser for help on using the repository browser.