source: tol/tol_pkg/MulVarFunAprox/DesignC2.tol @ 174f4b6

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

Counting number of nodal, stress and border points for each reference point

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@1792 be4846ae-f051-483a-831f-b64e7d08527f
  • Property mode set to 100644
File size: 10.9 KB
Line 
1/* -*- mode: C++ -*- */
2//////////////////////////////////////////////////////////////////////////////
3// FILE   : DesignC2.tol
4// PURPOSE: Defines Class @DesignC2
5//////////////////////////////////////////////////////////////////////////////
6
7
8//////////////////////////////////////////////////////////////////////////////
9Text _.autodoc.member.DesignC2 =
10"";
11Class @DesignC2
12//////////////////////////////////////////////////////////////////////////////
13{
14//Nodal points
15  VMatrix _.X;
16  Real _.n;
17  Real _.J;
18
19//////////////////////////////////////////////////////////////////////////////
20//Configuration parameters
21//////////////////////////////////////////////////////////////////////////////
22  Real prior.Z.order = 2;
23  Real min.contrast.surface = 3;
24
25//////////////////////////////////////////////////////////////////////////////
26//Auxiliar members
27//////////////////////////////////////////////////////////////////////////////
28
29  Real _.X.maxNeighbour = ?;
30  Real _.Z.maxNeighbour = ?;
31  VMatrix _.Y = Rand(0,2,0,0);
32  VMatrix _.Yj = Rand(0,1,0,0);
33  VMatrix _.Yt = Rand(0,1,0,0);
34  VMatrix _.Yb = Rand(0,1,0,0);
35  VMatrix _.Ym = Rand(0,1,0,0);
36  VMatrix _.Yn = Rand(0,1,0,0);
37  VMatrix _.Z = Rand(0,2,0,0);
38  VMatrix _.Z.h = Rand(0,1,0,0);
39  VMatrix _.Z.Y.neighbour= Rand(0,1,0,0);
40  VMatrix _.B = Rand(0,2,0,0);
41  VMatrix _.B.h = Rand(0,0,0,0);
42  Real _.K = 0;
43  Real _.H = 0;
44  Real _.Tz = 0;
45  Real _.Mz = 0;
46  Real _.Tb = 0;
47  Real _.Mb = 0;
48  Real _.M = 0;
49  Real _.N = 0;
50  Set _.base = Copy(Empty);
51  Real _.base.sel = ?;
52  VMatrix _.QX = Rand(0,0,0,0);
53  VMatrix _.QB = Rand(0,0,0,0);
54  VMatrix _.PZ = Rand(0,0,0,0);
55  VMatrix _.Q = Rand(0,0,0,0);
56  VMatrix _.LQ = Rand(0,0,0,0);
57  VMatrix _.Q.nonNullCellsByColumns = Rand(0,0,0,0);
58
59  //////////////////////////////////////////////////////////////////////////////
60  NameBlock get.base(Real sel)
61  //////////////////////////////////////////////////////////////////////////////
62  {
63    Real sel_ = If(Or(IsUnknown(sel),sel<=0,sel>Card(_.base)),_.base.sel);
64    _.base[sel_]
65  };
66
67  //////////////////////////////////////////////////////////////////////////////
68  //Virtual methods
69  //////////////////////////////////////////////////////////////////////////////
70  Real build.base(Real void);
71  Real build.reference.points(Real void);
72  Real build.stress.points(Real void) { False };
73  Real build.border.points(Real void) { False };
74  VMatrix closest.reference.index(VMatrix X);
75  Set get.X.weight(Set X.dif);
76  VMatrix match.evaluable(VMatrix X);
77
78  //////////////////////////////////////////////////////////////////////////////
79    Real _build.model.design(Real void)
80  //////////////////////////////////////////////////////////////////////////////
81  {
82    Real _.X.maxNeighbour := 1+2*_.n;
83    Real _.Z.maxNeighbour := If(prior.Z.order>=0,2,0);
84    Real build.reference.points(0);
85    Real build.stress.points(0);
86    Real build.border.points(0);
87    VMatrix _.Ym := _.Yj+((_.Yt+_.Yb)*Case(
88      prior.Z.order< 0, 0,
89      prior.Z.order==0, 1,
90      prior.Z.order==1, 1*(1+_.n),
91      prior.Z.order==2, 1*(1+_.n*3/2+(1/2)*_.n^2)));
92    Real _.M := _.J+_.Mb+_.Mz;
93    Real build.base(0);
94    Real _.N := _.H*_.K;
95    _.N<_.M
96  };
97
98  //////////////////////////////////////////////////////////////////////////////
99    Real build.model.design(Real void)
100  //////////////////////////////////////////////////////////////////////////////
101  {
102    _build.model.design(void)
103  };
104
105  //////////////////////////////////////////////////////////////////////////////
106    VMatrix build.obs.eval(VMatrix X, VMatrix X.h, FunRn2R::@BaseC2.Eval eval)
107  //////////////////////////////////////////////////////////////////////////////
108  {
109    Text _MID = "[@DesignC2::build.obs] ";
110    Real m = VRows(X);
111    WriteLn("TRACE "+_MID+"1 m="<<m);
112    If(!m, Rand(0,_.N,0,0), {
113    WriteLn("TRACE "+_MID+"2");
114    WriteLn("TRACE "+_MID+"3");
115    VMatrix X.j=DifEq(1/(1-B),Rand(m,1,1,1));
116    WriteLn("TRACE "+_MID+"4");
117    Set X.dif = For(1,_.X.maxNeighbour,VMatrix(Real i)
118    {
119      Set X.h.i = VCol2Set(SubCol(X.h,[[i]]));
120      X - SubRow(_.Y,X.h.i)
121    });
122    WriteLn("TRACE "+_MID+"5");
123    Set X.weight = get.X.weight(X.dif);
124    WriteLn("TRACE "+_MID+"6");
125    Matrix QX.triplet = Group("ConcatRows", For(1,_.X.maxNeighbour,Matrix(Real i)
126    {
127      VMatrix w.i = X.weight[i];
128      Set sel.rows = MatQuery::SelectMatch(NE(Drop(Abs(w.i),1.E-15),w.i*0));
129      VMatrix wi = SubRow(w.i,sel.rows);
130      VMatrix X.hiK = SubRow((SubCol(X.h,[[i]])-1)*_.K,sel.rows);
131      VMatrix Xdif.i = SubRow(X.dif[i],sel.rows);
132      VMatrix X.j.i = SubRow(X.j,sel.rows);
133      Real eval::upgrade.cache(Xdif.i);
134      Group("ConcatRows", For(1,_.K,Matrix(Real k)
135      {
136        VMatrix pk = eval::evaluate.v(k,Xdif.i);
137        VMatrix qk = wi$*pk;
138        VMat2Mat(X.j.i | (X.hiK+k) | qk)
139      }))
140    }));
141    WriteLn("TRACE "+_MID+"7");
142    Pack(Convert(Triplet(QX.triplet, m, _.N),"Cholmod.R.Sparse"))
143  })};
144
145  //////////////////////////////////////////////////////////////////////////////
146    VMatrix build.obs(VMatrix X, VMatrix X.h)
147  //////////////////////////////////////////////////////////////////////////////
148  {
149    FunRn2R::@BaseC2.Eval.Function eval =
150     FunRn2R::@BaseC2.Eval.Function::New(get.base(0));
151    build.obs.eval(X,X.h,eval)
152  };
153  //////////////////////////////////////////////////////////////////////////////
154    VMatrix build.obs.partial.derivative.1(VMatrix X, VMatrix X.h, Real s)
155  //////////////////////////////////////////////////////////////////////////////
156  {
157    FunRn2R::@BaseC2.Eval.Partial.Derivative.1 eval =
158     FunRn2R::@BaseC2.Eval.Partial.Derivative.1::New(get.base(0), s);
159    build.obs.eval(X,X.h,eval)
160  };
161  //////////////////////////////////////////////////////////////////////////////
162    VMatrix build.obs.partial.derivative.2(VMatrix X, VMatrix X.h,
163                                           Real s1, Real s2)
164  //////////////////////////////////////////////////////////////////////////////
165  {
166    FunRn2R::@BaseC2.Eval.Partial.Derivative.2 eval =
167     FunRn2R::@BaseC2.Eval.Partial.Derivative.2::New(get.base(0), s1, s2);
168    build.obs.eval(X,X.h,eval)
169  };
170  //////////////////////////////////////////////////////////////////////////////
171    VMatrix build.obs.partial.integrate.1(VMatrix X, VMatrix X.h, Real s)
172  //////////////////////////////////////////////////////////////////////////////
173  {
174    FunRn2R::@BaseC2.Eval.Partial.Integrate.1 eval =
175     FunRn2R::@BaseC2.Eval.Partial.Integrate.1::New(get.base(0), s);
176    build.obs.eval(X,X.h,eval)
177  };
178  //////////////////////////////////////////////////////////////////////////////
179    VMatrix build.border.prior(Real void)
180  //////////////////////////////////////////////////////////////////////////////
181  {
182    If(!_.Tb, Rand(0,_.N,0,0), {
183    VMatrix B.t=DifEq(1/(1-B),Rand(_.Tb,1,1,1));
184    VMatrix Y = SubRow(_.Y, VCol2Set(_.B.h));
185    VMatrix B.dif = _.B - Y;
186    FunRn2R::@BaseC2 b = get.base(0);   
187    Matrix PB.triplet =
188    {
189      VMatrix B.hK = (_.B.h-1)*_.K;
190      Real b::upgrade.cache(B.dif,prior.Z.order);
191      Real lastRow = 0;
192      Group("ConcatRows", For(1,_.K,Matrix(Real k)
193      {
194        Real cumRows = lastRow;
195        VMatrix pk = b::function.v(k,B.dif);
196        VMatrix col.idx = (B.hK+k);
197        {
198          VMatrix row.idx = B.t+cumRows;
199          Real cumRows := cumRows + _.Tb;
200          VMat2Mat(row.idx | col.idx | pk)
201        }  << If(prior.Z.order<1,Rand(0,3,0,0),
202        Group("ConcatRows",For(1,_.n,Matrix(Real d)
203        {
204          VMatrix pkd = b::partial.derivative.1.v(k,B.dif,d);
205          VMatrix row.idx = B.t+cumRows;
206          Real cumRows := cumRows + _.Tb;
207          VMat2Mat(row.idx | col.idx | pkd)
208        })))<< If(prior.Z.order<2,Rand(0,3,0,0),
209        Group("ConcatRows",For(1,_.n,Matrix(Real d1)
210        {
211          Group("ConcatRows",For(1,d1,Matrix(Real d2)
212          {
213            VMatrix pkd2 = b::partial.derivative.2.v(k,B.dif,d1,d2);
214            VMatrix row.idx = B.t+cumRows;
215            Real cumRows := cumRows + _.Tb;
216            VMat2Mat(row.idx | col.idx | pkd2)
217          }))
218        })))
219      }))
220    };
221    Convert(Triplet(PB.triplet, _.Mb, _.N),"Cholmod.R.Sparse")
222  })};
223  //////////////////////////////////////////////////////////////////////////////
224    VMatrix build.stress.prior(Real void)
225  //////////////////////////////////////////////////////////////////////////////
226  {
227    If(!_.Tz, Rand(0,_.N,0,0), {
228    VMatrix Z.t=DifEq(1/(1-B),Rand(_.Tz,1,1,1));
229   
230    Set Z.dif = For(1,_.Z.maxNeighbour,VMatrix(Real i)
231    {
232      _.Z - SubRow(_.Y,VCol2Set(SubCol(_.Z.h,[[i]])))
233    });
234    FunRn2R::@BaseC2 b = get.base(0);   
235    Matrix PZ.triplet = Group("ConcatRows", For(1,_.Z.maxNeighbour,Matrix(Real i)
236    {
237      VMatrix Z.hiK = (SubCol(_.Z.h,[[i]])-1)*_.K;
238      VMatrix Zdif.i = Z.dif[i];
239      VMatrix Z.t.i = Z.t;
240      Real b::upgrade.cache(Zdif.i,prior.Z.order);
241      Real lastRow = 0;
242      Group("ConcatRows", For(1,_.K,Matrix(Real k)
243      {
244        Real cumRows = lastRow;
245        VMatrix pk = b::function.v(k,Zdif.i)*If(i==1,+1,-1);
246        VMatrix col.idx = (Z.hiK+k);
247        {
248          VMatrix row.idx = Z.t.i+cumRows;
249          Real cumRows := cumRows + _.Tz;
250          VMat2Mat(row.idx | col.idx | pk)
251        }  << If(prior.Z.order<1,Rand(0,3,0,0),
252        Group("ConcatRows",For(1,_.n,Matrix(Real d)
253        {
254          Real sign = If(i==1,+1/2,-1/2);
255          VMatrix pkd = b::partial.derivative.1.v(k,Zdif.i,d)*sign;
256          VMatrix row.idx = Z.t.i+cumRows;
257          Real cumRows := cumRows + _.Tz;
258          VMat2Mat(row.idx | col.idx | pkd)
259        })))<< If(prior.Z.order<2,Rand(0,3,0,0),
260        Group("ConcatRows",For(1,_.n,Matrix(Real d1)
261        {
262          Group("ConcatRows",For(1,d1,Matrix(Real d2)
263          {
264            Real sign = If(i==1,+1/4,-1/4);
265            VMatrix pkd2 = b::partial.derivative.2.v(k,Zdif.i,d1,d2)*sign;
266            VMatrix row.idx = Z.t.i+cumRows;
267            Real cumRows := cumRows + _.Tz;
268            VMat2Mat(row.idx | col.idx | pkd2)
269          }))
270        })))
271      }))
272    }));
273    Convert(Triplet(PZ.triplet, _.Mz, _.N),"Cholmod.R.Sparse")
274  })};
275  //////////////////////////////////////////////////////////////////////////////
276    Real build.decomposition(Real void)
277  //////////////////////////////////////////////////////////////////////////////
278  {
279    Text _MID = "[@DesignC2::build.decomposition] ";
280    WriteLn("TRACE "+_MID+"1");
281    VMatrix X.h=closest.reference.index(_.X);
282    VMatrix _.QX := build.obs(_.X, X.h);
283    WriteLn("TRACE "+_MID+"2");
284    VMatrix _.QB := build.border.prior(0);
285    WriteLn("TRACE "+_MID+"3");
286    VMatrix _.PZ := build.stress.prior(0);
287    WriteLn("TRACE "+_MID+"4");
288    VMatrix _.Q := _.QX << _.QB << _.PZ;
289    WriteLn("TRACE "+_MID+"5");
290    VMatrix _.Q.nonNullCellsByColumns :=
291    {
292      VMatrix aux = Tra(Rand(1,_.M,1,1)*Or(_.Q,Zeros(_.M,_.N)));
293      Group("ConcatColumns",For(1,_.K,VMatrix(Real k)
294      {
295        SubRow(aux,Range(k,_.N,_.K))
296      }))
297    };
298    WriteLn("TRACE "+_MID+"6");
299    VMatrix _.LQ := CholeskiFactor(_.Q,"XtX");
300    WriteLn("TRACE "+_MID+"7");
301    True
302  }
303};
304
Note: See TracBrowser for help on using the repository browser.