source: tol/tol_pkg/MulVarFunAprox/DesignC2.tol @ 2fb8329

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

Preparing for derivation and integration methods

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