source: tol/tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.Grid.tol @ b050cc43

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

New class hierarchy

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@1762 be4846ae-f051-483a-831f-b64e7d08527f
  • Property mode set to 100644
File size: 11.4 KB
Line 
1/* -*- mode: C++ -*- */
2//////////////////////////////////////////////////////////////////////////////
3// FILE   : DesignC2.tol
4// PURPOSE: Defines Class @DesignC2
5//////////////////////////////////////////////////////////////////////////////
6
7
8//////////////////////////////////////////////////////////////////////////////
9Text _.autodoc.member.DesignC2.Monomial.Spheric.Grid =
10"";
11Class @DesignC2.Monomial.Spheric.Grid :
12  @DesignC2.Monomial,
13  @DesignC2.Spheric
14//////////////////////////////////////////////////////////////////////////////
15{
16//////////////////////////////////////////////////////////////////////////////
17//Configuration parameters
18//////////////////////////////////////////////////////////////////////////////
19 
20  Real force.reference.near.nodal = True;
21  Real _.radius = ?;
22  Real _.grid.size = ?;
23
24//////////////////////////////////////////////////////////////////////////////
25//Auxiliar members
26//////////////////////////////////////////////////////////////////////////////
27  Real _.K0 = ?;
28  Set _.X.col.info = Copy(Empty);
29  Set _.range = Copy(Empty);
30  VMatrix _.grid = Rand(0,0,0,0);
31  Real _.G = ?;
32  Matrix _.grid.active = Rand(0,0,0,0);
33  Matrix _.grid.inactive = Rand(0,0,0,0);
34  VMatrix _.Y.match = Rand(0,0,0,0);
35  VMatrix _.Y.h = Rand(0,0,0,0);
36  VMatrix _.Y.neighbourhood = Rand(0,0,0,0);
37
38  ////////////////////////////////////////////////////////////////////////////
39  Static @DesignC2.Monomial.Spheric.Grid New(VMatrix X, Real radius.margin)
40  ////////////////////////////////////////////////////////////////////////////
41  {
42    @DesignC2.Monomial.Spheric.Grid new =
43    [[
44      VMatrix _.X = X;
45      Real    _.n = VColumns(X);
46      Real    _.J = VRows(X);
47      Real _.radius.margin = radius.margin
48    ]]
49  };
50
51  //////////////////////////////////////////////////////////////////////////////
52  VMatrix closest.grid.index(VMatrix X)
53  //////////////////////////////////////////////////////////////////////////////
54  {
55    Real m = VRows(X);
56    VMatrix cumulated = Rand(m,1,1,1);
57    VMatrix X.position = Rand(m,1,1,1)+SetSum(For(1,_.n,VMatrix(Real i_)
58    {
59      Real i = _.n-i_+1;
60      Real cRi = Card(_.range[i]);
61      Set info = _.X.col.info[i];
62      VMatrix Xi = SubCol(X,[[i]]);
63      VMatrix Yi.min = Rand(m,1,_.range[i][1],_.range[i][1]);
64      VMatrix Ri = Round((Xi-Yi.min)/_.grid.size);
65      VMatrix minRi = Rand(m,1,0,0);
66      VMatrix maxRi = Rand(m,1,cRi-1,cRi-1);
67      VMatrix aux = Min(maxRi,Max(minRi,Ri));
68      VMatrix pos = aux $* cumulated;
69      VMatrix cumulated := cumulated*cRi;
70      pos
71    }))
72  };
73
74  //////////////////////////////////////////////////////////////////////////////
75  VMatrix closest.reference.index(VMatrix X)
76  //////////////////////////////////////////////////////////////////////////////
77  {
78    Real m = VRows(X);
79    VMatrix X.closest.grid.index =
80    {
81      Matrix aux.1 = VMat2Mat(closest.grid.index(X));
82      Matrix aux.2 = -MatReplace(aux.1,_.grid.active,-VMat2Mat(_.Y.h));
83      Mat2VMat(aux.2)
84    };
85    SubRow(_.Y.neighbourhood, VCol2Set(X.closest.grid.index))
86  };
87
88  //////////////////////////////////////////////////////////////////////////////
89  Real build.reference.points(Real void)
90  //////////////////////////////////////////////////////////////////////////////
91  {
92    Text _MID = "[@DesignC2.Monomial.Spheric.Grid::build.reference.points] ";
93    WriteLn("TRACE "+_MID+"1");
94    Set _.X.col.info := For(1,_.n,Set(Real i)
95    {[[
96      VMatrix col = SubCol(_.X,[[i]]);
97      Real min = VMatMin(col);
98      Real max = VMatMax(col);
99      Real length = max-min
100    ]]});
101   
102    Real If(Or(IsUnknown(_.grid.size),_.grid.size<=0),
103    {
104      Real length.max = SetMax(Traspose(Extract(_.X.col.info,4))[1]);
105      Real density = _.J/SetProd(Traspose(Extract(_.X.col.info,4))[1]);
106      _.grid.size := 0.5*(_.K0/density)^(1/_.n)
107    });
108   
109    Real _.radius := (1+_.radius.margin)*_.grid.size/Sqrt(2);
110   
111    WriteLn("TRACE "+_MID+"2");
112    Set _.range := For(1,_.n,Set(Real i)
113    {
114      Set info = _.X.col.info[i];
115      Real points = 1+Ceil(info::length/_.grid.size);
116      Real exceed = (points-1)*_.grid.size-info::length;
117      Matrix rnk = Tra(DifEq(1/(1-B),Rand(points,1,1,1))-Rand(points,1,1,1));   
118      Set MatSet((rnk*_.grid.size)+info::min-exceed/2)[1]
119    });
120    WriteLn("TRACE "+_MID+"3");
121   
122    VMatrix _.grid := Mat2VMat(SetMat(Group("CartProd",_.range)));
123    Real _.G := VRows(_.grid);
124    WriteLn("TRACE "+_MID+"4");
125    Matrix _.grid.active := If(!force.reference.near.nodal,
126    {
127      DifEq(1/(1-B),Rand(_.G,1,1,1))
128    },
129    {
130      VMatrix X.closest.grid.index = closest.grid.index(_.X);
131      SetCol(Sort(Unique(VCol2Set(X.closest.grid.index)),Compare))
132    });
133    WriteLn("TRACE "+_MID+"5");
134    Matrix _.grid.inactive := If(!force.reference.near.nodal,
135    {
136      Rand(0,1,0,0)
137    },
138    {
139      SetCol(Range(1,_.G,1)-Col2Set(_.grid.active))
140    });
141    WriteLn("TRACE "+_MID+"6");
142    VMatrix _.Y := SubRow(_.grid,Col2Set(_.grid.active));
143    Real _.H := VRows(_.Y);
144    WriteLn("TRACE "+_MID+"7");
145    VMatrix _.Y.match := Convert(Triplet(_.grid.active|Rand(Rows(_.grid.active),2,1,1),_.G,1),
146      "Cholmod.R.Sparse");
147   
148    VMatrix Y.min = Group("ConcatColumns", For(1,_.n,VMatrix(Real i)
149    {
150      Real aux = _.range[i][1];
151      Rand(_.H,1,aux,aux)
152    }));
153   
154    VMatrix Y.max = Group("ConcatColumns", For(1,_.n,VMatrix(Real i)
155    {
156      Real aux = _.range[i][Card(_.range[i])];
157      Rand(_.H,1,aux,aux)
158    }));
159   
160    VMatrix Y.index = Mat2VMat(SetMat(Group("CartProd",For(1,_.n,Set(Real i)
161    {
162      Range(1,Card(_.range[i]),1)
163    }))));
164    WriteLn("TRACE "+_MID+"8");
165   
166    VMatrix Y.closest.grid.index = closest.grid.index(_.Y);
167   
168    VMatrix Y.neighbour.pos = Min(Y.max, _.Y + Rand(_.H,_.n,_.grid.size,_.grid.size));
169    VMatrix Y.neighbour.neg = Max(Y.min, _.Y - Rand(_.H,_.n,_.grid.size,_.grid.size));
170   
171    VMatrix _.Y.h := DifEq(1/(1-B),Rand(_.H,1,1,1));
172   
173    WriteLn("TRACE "+_MID+"9");
174    VMatrix _.Y.neighbourhood := _.Y.h | If(!force.reference.near.nodal,
175    {
176      Set aux.1 = BinGroup("<<", For(1,_.n,Set(Real i)
177      {
178        VMatrix pre  = SubCol(_.Y,Range(1,  i-1,1));
179        VMatrix post = SubCol(_.Y,Range(i+1,_.n,  1));
180        SetOfVMatrix(pre | SubCol(Y.neighbour.neg,[[i]]) | post,
181                     pre | SubCol(Y.neighbour.pos,[[i]]) | post)
182      }));
183      Group("ConcatColumns",EvalSet(aux.1,closest.grid.index))
184    },
185    {
186      Set aux.1 = BinGroup("<<", For(1,_.n,Set(Real i)
187      {
188        VMatrix pre    = SubCol(_.Y,Range(1,  i-1,1));
189        VMatrix post   = SubCol(_.Y,Range(i+1,_.n,  1));
190        VMatrix pre.1  = pre | SubCol(Y.neighbour.neg,[[i]]) | post;
191        VMatrix post.1 = pre | SubCol(Y.neighbour.pos,[[i]]) | post;
192        [[pre.1, post.1]]
193      }));
194      Matrix Y.cls = VMat2Mat(_.Y.h);
195      Matrix Y.inact0 = Rand(Rows(_.grid.inactive),1,0,0);
196      Group("ConcatColumns",EvalSet(aux.1,VMatrix(VMatrix g)
197      {
198        Matrix aux.1 = VMat2Mat(closest.grid.index(g));
199        Matrix aux.2 = MatReplace(aux.1,_.grid.active,-Y.cls);
200        Matrix aux.3 = MatReplace(aux.2,_.grid.inactive,Y.inact0);
201        Matrix aux.4 = IfMat(aux.3,-aux.3,Y.cls);
202        Mat2VMat(aux.4)
203      }))
204    });
205    _.H
206  };
207
208  //////////////////////////////////////////////////////////////////////////////
209  Real build.stress.points(Real void)
210  //////////////////////////////////////////////////////////////////////////////
211  {
212    Text _MID = "[@DesignC2.Monomial.Spheric.Grid::build.stress.points] ";
213    WriteLn("TRACE "+_MID+"1");
214    Matrix Y.neighbour.triplet = Group("ConcatRows", For(1,_.H,Matrix(Real h)
215    {
216      Set neighbours = Unique(VRow2Set(SubRow(_.Y.neighbourhood,[[h]])));
217      Set neighbours.up = Select(neighbours,Real(Real s) { s>h });
218      Matrix s = SetCol(Sort(neighbours.up,Compare));
219      Rand(Rows(s),1,h,h)|s|Rand(Rows(s),1,1,1)
220    }));
221   
222    WriteLn("TRACE "+_MID+"2");
223    VMatrix Y.neighbour.1.map = Triplet(Y.neighbour.triplet,_.H,_.H);
224    VMatrix Y.neighbour.2.map = Y.neighbour.1.map*Y.neighbour.1.map;
225    VMatrix Y.neighbour.12.map = Max(Y.neighbour.1.map,Y.neighbour.2.map);
226    Matrix Y.neighbour.12.triplet = VMat2Triplet(Y.neighbour.12.map);
227    VMatrix Z.h.a = Mat2VMat(SubCol(Y.neighbour.12.triplet,[[1]]));
228    VMatrix Z.h.b = Mat2VMat(SubCol(Y.neighbour.12.triplet,[[2]]));
229    VMatrix Y.neighbour.a = SubRow(_.Y,VCol2Set(Z.h.a));
230    VMatrix Y.neighbour.b = SubRow(_.Y,VCol2Set(Z.h.b));
231    VMatrix Z.all = (Y.neighbour.a + Y.neighbour.b)*0.5;
232   
233    WriteLn("TRACE "+_MID+"3");
234    VMatrix Z.interior.a =
235    {
236      VMatrix Z.norm2 = RPow(Y.neighbour.a-Z.all,2)*Rand(_.n,1,1,1);
237      VMatrix Z.tao = RPow(Z.norm2,0.5)/_.radius;
238      LE(Z.tao,Rand(VRows(Z.tao),1,1,1))
239    };
240    VMatrix Z.interior.b =
241    {
242      VMatrix Z.norm2 = RPow(Y.neighbour.b-Z.all,2)*Rand(_.n,1,1,1);
243      VMatrix Z.tao = RPow(Z.norm2,0.5)/_.radius;
244      LE(Z.tao,Rand(VRows(Z.tao),1,1,1))
245    };
246    VMatrix Z.interior = Z.interior.a $* Z.interior.b;
247   
248    WriteLn("TRACE "+_MID+"4");
249    VMatrix _.Z := SubRow(Z.all,MatQuery::SelectMatch(Z.interior));
250    VMatrix _.Z.h := SubRow(Z.h.a | Z.h.b,MatQuery::SelectMatch(Z.interior));
251   
252    Real _.Tz := VRows(_.Z);
253    Real _.Mz := Case(
254      prior.Z.order< 0, 0,
255      prior.Z.order==0, _.Tz,
256      prior.Z.order==1, _.Tz*(1+_.n),
257      prior.Z.order==2, _.Tz*(1+_.n*3/2+(1/2)*_.n^2));
258    _.Tz
259  };
260
261  //////////////////////////////////////////////////////////////////////////////
262    Real build.model.design(Real void)
263  //////////////////////////////////////////////////////////////////////////////
264  {
265    Text _MID = "[@DesignC2.Monomial.Spheric.Grid::build.model.design] ";
266    FunRn2R::@BaseC2 base.max =
267      FunRn2R::@BaseC2.Monomial::New(_.n,max.deg);
268    Set Append(_.base, [[base.max]]);
269    Real _.base.sel := Card(_.base);
270    Real _.K0 := base.max::_.K;
271    _build.model.design(void)
272  };
273
274  ///////////////////////////////////////////////////////////////////////////
275    VMatrix generate.internal.points(
276      Real m,
277      Real force.reference.near.nodal)
278  ///////////////////////////////////////////////////////////////////////////
279  {
280    If(force.reference.near.nodal,
281    {
282      VMatrix h = Round(Rand(m,1,0.5,_.H+0.499999999));
283      VMatrix y = SubRow(_.Y,VCol2Set(h));
284      VMatrix v = Gaussian(m,_.n,0,1);
285      VMatrix u = v*(VMatDat(MtMSqr(v),1,1)^-0.5);
286      y+u*(Rand(0,1)*_.radius)
287    },
288    {
289      VMatrix h = Round(Rand(m,1,0.5,_.G+0.499999999));
290      VMatrix y = SubRow(_.grid,VCol2Set(h));
291      VMatrix v = Gaussian(m,_.n,0,1);
292      VMatrix u = v*(VMatDat(MtMSqr(v),1,1)^-0.5);
293      y+u*(Rand(0,1)*_.radius)
294    })
295  };
296
297  ////////////////////////////////////////////////////////////////////////////
298    VMatrix match.evaluable(VMatrix X)
299  ////////////////////////////////////////////////////////////////////////////
300  {
301    Real m = VRows(X);
302    Real G.0 = Rows(_.grid.inactive);
303    VMatrix X.closest.grid.index = closest.grid.index(X);
304    If(G.0==0,
305    {
306      VMatrix X.h=SubRow(SubCol(_.Y.neighbourhood,[[1]]),
307                         VCol2Set(X.closest.grid.index));
308      VMatrix X.dif = X - SubRow(_.Y,VCol2Set(X.h));
309      VMatrix X.norm2 = RPow(X.dif,2)*Rand(_.n,1,1,1);
310      VMatrix X.tao = RPow(X.norm2,0.5)/_.radius;
311      LE(X.tao,Rand(m,1,1,1))
312    },
313    {
314      Matrix true_ = Rand(_.H,1,1,1);
315      Matrix false_ = Rand(G.0,1,0,0);
316      Matrix aux.1 = VMat2Mat(X.closest.grid.index);
317      Matrix aux.2 = MatReplace(aux.1,_.grid.active,-true_);
318      Matrix aux.3 = -MatReplace(aux.2,_.grid.inactive,false_);
319      Mat2VMat(aux.3)
320    })
321  }
322
323};
Note: See TracBrowser for help on using the repository browser.