source: tol/tol_pkg/MulVarFunAprox/DesignC2.tol @ 4bd88af

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