source: tol/tol_pkg/MulVarFunAprox/test/test_0002/test.tol @ 4a13114

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

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@1896 be4846ae-f051-483a-831f-b64e7d08527f
  • Property mode set to 100644
File size: 8.8 KB
Line 
1/* -*- mode: C++ -*- */
2//////////////////////////////////////////////////////////////////////////////
3// FILE   : test.tol
4// PURPOSE: Package MulVarFunAprox test
5//////////////////////////////////////////////////////////////////////////////
6
7Text email_ = "vdebuen@tol-project.org";
8Text link_ = "";
9Text summary_ = "Test for mutivariant function aproximation";
10
11//Write here the initial test TOL code if needed
12Real numErr0 = Copy(NError);
13Real numWar0 = Copy(NWarning);
14
15//Random seed setting
16Real PutRandomSeed(0);
17Real rndSeed = GetRandomSeed(0);
18WriteLn("Current Random Seed = "<<rndSeed);
19
20#Require MulVarFunAprox;
21//#Embed "../../MulVarFunAprox.tol";
22
23
24//////////////////////////////////////////////////////////////////////////////
25VMatrix generate.nodal.points(Real J, Real n)
26//////////////////////////////////////////////////////////////////////////////
27{
28  Group("ConcatColumns",For(1,n,VMatrix(Real i)
29  {
30    Real min = Rand(2,4);
31    Real len = Rand(5,15)*((i/n)^(1/n));
32    Rand(J,1,min,min+len)
33  }))
34};
35VMatrix X = generate.nodal.points(500,2);
36/* * /
37MulVarFunAprox::@DesignC2 design = {
38  MulVarFunAprox::@DesignC2.Monomial.Spheric.Grid::New(
39    X, radius.margin=0.1) };
40/* */
41WriteLn("TRACE 1");
42MulVarFunAprox::@DesignC2 design = {
43  MulVarFunAprox::@DesignC2.Monomial.Spheric.MeshFree::New(
44    VMatrix X,
45    Real maxNeighbour = 4*VColumns(X),
46    Real stress.maxDistance_ = 5) };
47/* */
48WriteLn("TRACE 2");
49Real design::max.deg := 8;
50Real design::prior.Z.order := 2;
51Real design::use.compact.base := False;
52Real design::min.contrast.surface := 1.0;
53
54Real t0 = Copy(Time);
55WriteLn("Building model design for "<<design::_.J+" points "<<
56  "of dimension "<<design::_.n+" ...");
57Real deg = design::build.model.design(0);
58Real t.building = Copy(Time)-t0;
59WriteLn("Built model design in "<<t.building+" seconds");
60
61Real t0 := Copy(Time);
62WriteLn("Building design decomposition for design matrix ("<<
63  design::_.M+"x"<<design::_.N+") ...");
64Real design::build.decomposition(0);
65Real t.decomposition = Copy(Time)-t0;
66WriteLn("Input decomposed in "<<t.decomposition+" seconds");
67
68FunRn2R::@BaseC2 b = design::get.base(0);
69
70FunRn2R::@FunC2 base.random.fun=FunRn2R::@FunC2.InBaseC2::Random(b, 1);
71
72
73//////////////////////////////////////////////////////////////////////////////
74VMatrix fun.eval.intra.base(VMatrix X)
75//////////////////////////////////////////////////////////////////////////////
76{
77  VMatrix X1 =SubCol(X,[[1]]);
78  VMatrix X2 =SubCol(X,[[2]]);
79  VMatrix one = X1*0+1;
80  one*7 -X1 +X2*2 -X1$*X1*3 +X1$*X2 +X1$*X1$*X1*0.4 -X1$*X2$*X2*0.1
81};
82//////////////////////////////////////////////////////////////////////////////
83VMatrix der.1.eval.intra.base(VMatrix X)
84//////////////////////////////////////////////////////////////////////////////
85{
86  VMatrix X1 =SubCol(X,[[1]]);
87  VMatrix X2 =SubCol(X,[[2]]);
88  VMatrix one = X1*0+1;
89    -one    -X1*6  +X2 +X1$*X1*1.2 - X2$*X2*0.1
90};
91//////////////////////////////////////////////////////////////////////////////
92VMatrix der.2.eval.intra.base(VMatrix X)
93//////////////////////////////////////////////////////////////////////////////
94{
95  VMatrix X1 =SubCol(X,[[1]]);
96  VMatrix X2 =SubCol(X,[[2]]);
97  VMatrix one = X1*0+1;
98   one*2 +X2 -X1$*X2*0.2
99};
100
101//////////////////////////////////////////////////////////////////////////////
102VMatrix fun.eval.extra.base(VMatrix X)
103//////////////////////////////////////////////////////////////////////////////
104{
105  VMatrix X1 =SubCol(X,[[1]]);
106  VMatrix X2 =SubCol(X,[[2]]);
107  Sin(X1)$*Cos(X2)
108};
109//////////////////////////////////////////////////////////////////////////////
110VMatrix der.1.eval.extra.base(VMatrix X)
111//////////////////////////////////////////////////////////////////////////////
112{
113  VMatrix X1 =SubCol(X,[[1]]);
114  VMatrix X2 =SubCol(X,[[2]]);
115  Cos(X1)$*Cos(X2)
116};
117//////////////////////////////////////////////////////////////////////////////
118VMatrix der.2.eval.extra.base(VMatrix X)
119//////////////////////////////////////////////////////////////////////////////
120{
121  VMatrix X1 =SubCol(X,[[1]]);
122  VMatrix X2 =SubCol(X,[[2]]);
123  -Sin(X1)$*Sin(X2)
124};
125
126Real use_extra = True;
127
128Code fun.eval   = If(use_extra,  fun.eval.extra.base,  fun.eval.intra.base);
129Code der.1.eval = If(use_extra,der.1.eval.extra.base,der.1.eval.intra.base);
130Code der.2.eval = If(use_extra,der.2.eval.extra.base,der.2.eval.intra.base);
131
132//Real b::upgrade.cache(design::_.X,0,0);
133//Code fun.eval = base.random.fun::eval.function.v;
134
135VMatrix F_X= fun.eval(design::_.X);
136
137Real box.cox.exp = 1;
138MulVarFunAprox::@ModelC2 model = MulVarFunAprox::@ModelC2::New(
139   design,
140   RPow(F_X,box.cox.exp));
141WriteLn("Running regression model for design matrix "<<design::_.Q);
142Real t0 := Copy(Time);
143Real model::estimate(0);
144Real t.estimating = Copy(Time)-t0;
145WriteLn("Estimated model in "<<t.estimating+" seconds");
146
147/*
148VMatrix dF_X1 = der.1.eval(design::_.X);
149VMatrix dF_X2 = der.2.eval(design::_.X);
150
151VMatrix dF_X1_ = model::evaluate.partial.derivative.1(design::_.X, 1);
152VMatrix dF_X2_ = model::evaluate.partial.derivative.1(design::_.X, 2);
153
154VMatrix dF_X1_cmp = dF_X1 | dF_X1_;
155VMatrix dF_X2_cmp = dF_X2 | dF_X2_;
156
157VMatrix dF_X1_err = dF_X1 - dF_X1_;
158VMatrix dF_X2_err = dF_X2 - dF_X2_;
159
160Real dF_X1_err.avr = VMatAvr (dF_X1_err);
161Real dF_X1_err.ste = VMatStDs(dF_X1_err);
162Real dF_X2_err.avr = VMatAvr (dF_X2_err);
163Real dF_X2_err.ste = VMatStDs(dF_X2_err);
164
165*/
166VMatrix X2 = { design::generate.internal.points(Real m = design::_.J/2)};
167
168Real b::upgrade.cache(X2,0,0);
169VMatrix Fx2 = fun.eval(X2);
170Real t0 := Copy(Time);
171WriteLn("Evaluating forecast for "<<VRows(X2)+" points...");
172VMatrix F2.frc = model::evaluate(X2,0);
173Real t.forecasting = Copy(Time)-t0;
174WriteLn("Evaluated forecast in "<<t.forecasting+" seconds");
175WriteLn("Forecast time by point "<<(t.forecasting/VRows(X2))+" seconds");
176VMatrix F2.cmp = Fx2 | F2.frc;
177VMatrix E2 = Fx2 - F2.frc;
178Real E2.avr = VMatAvr(E2);
179Real E2.ste = VMatStDs(E2);
180
181Real q.1 = 1-model::_.E.ste;
182Real q.2 = 1-E2.ste;
183Set View([[
184  Real design::_.n, design::_.K, design::_.J, design::_.H,
185  design::_.X.maxNeighbour, design::_.Z.maxNeighbour,
186  design::_.Tz, design::_.N, design::_.Mz, design::_.M,
187  model::_.E.avr, model::_.E.ste,
188  "\n"]],"");
189
190Set View([[
191  E2.avr, E2.ste,
192  "\n"]],"");
193
194VMatrix show.X.Y =
195  ((design::_.X*0) <<  design::_.Y   ) |
196  ( design::_.X    << (design::_.Y*0));
197
198VMatrix show.X.Y.Z =
199  ( design::_.X    << (design::_.Y*0) << (design::_.Z*0)) |
200  ((design::_.X*0) <<  design::_.Y    << (design::_.Z*0)) |
201  ((design::_.X*0) << (design::_.Y*0) <<  design::_.Z   );
202
203VMatrix show.Y.Z =
204  ( design::_.Y    << (design::_.Z*0)) |
205  ((design::_.Y*0) <<  design::_.Z   );
206
207
208VMatrix Y.n.m = model::_.design::_.Yn | model::_.design::_.Ym;
209
210VMatrix xyz = X|F_X;
211
212Matrix MatWriteFile("real.xyz.bbm",VMat2Mat(xyz));
213
214VMatrix _.X = X;
215Real _.n = VColumns(_.X);
216    Set _.X.col.info = For(1,_.n,Set(Real i)
217    {[[
218      VMatrix col = SubCol(_.X,[[i]]);
219      Real min = VMatMin(col);
220      Real max = VMatMax(col);
221      Real length = max-min;
222      Set Range(min,max,length/99)
223    ]]});
224VMatrix grid.X = Mat2VMat(SetMat(Group("CartProd",Traspose(Extract(_.X.col.info,5))[1])));
225VMatrix grid.F = RPow(model::evaluate(grid.X,0),1/box.cox.exp);
226
227//Rango de latitudes
228Matrix x = SetCol(_.X.col.info[1][5]);
229//Rango de latitudes
230Matrix y = SetCol(_.X.col.info[2][5]);
231Matrix z =
232{
233  Matrix aux= Group("ConcatColumns",For(1,Rows(x),Matrix(Real j)
234  {
235    VMat2Mat(SubRow(grid.F, Range(1+(j-1)*Rows(y),(j)*Rows(y),1)))
236  }));
237  Tra(aux)
238};
239
240Matrix MatWriteFile("aprox.x.bbm",x);
241Matrix MatWriteFile("aprox.y.bbm",y);
242Matrix MatWriteFile("aprox.z.bbm",z);
243
244Set BMTFile([[z]],"aprox.grid.csv");
245
246Text currentDir = GetAbsolutePath(GetFilePath(GetSourcePath(z)));
247Text WriteFile("show.R",
248"
249library(rgl)
250
251################################################################################
252read.bbm <-
253################################################################################
254function(fbbm) {
255  fd = file(fbbm, open='rb')
256  dim <- readBin(fd,integer(), n=2)
257  data<-readBin(fd,double(),n=dim[1]*dim[2])
258  close(fd)
259  matrix(data, nrow=dim[1], ncol=dim[2], byrow=TRUE)
260}
261
262x<-t(read.bbm(\""+currentDir+"aprox.x.bbm\"))
263 y<-t(read.bbm(\""+currentDir+"aprox.y.bbm\"))
264 z<-read.bbm(\""+currentDir+"aprox.z.bbm\")
265xyz.real <- read.bbm(\""+currentDir+"real.xyz.bbm\")
266open3d()
267bg3d(\"white\")
268
269persp3d(x[,],y[,],z, xlim = range(x, na.rm = TRUE), ylim = range(y, na.rm = TRUE), col=\"green\")
270points3d(xyz.real)
271");
272
273/*
274Real numErr1 = Copy(NError);
275Real numWar1 = Copy(NWarning);
276
277Set partialResults_ = [[numErr0, numErr1,
278                        q.1, q.2]];
279
280//This is a messure of the success of the test
281Real quality_ = And(numErr1 == numErr0,
282                    numWar1 == numWar0)*Min(q.1,q.2);
283
284//Return the results
285Set resultStr_ = @strTestStatus(summary_, link_, quality_,
286                  "Partial results = "<<partialResults_,
287                  "NO DBServerType", "NO DBServerHost", "NO DBAlias",
288                  email_);
289WriteLn(""<<resultStr_);
290resultStr_;
291
292/* */
293
Note: See TracBrowser for help on using the repository browser.