Changeset b050cc43 in tol for tol_pkg/MulVarFunAprox
- Timestamp:
- Dec 4, 2009, 4:44:16 PM (17 years ago)
- Branches:
- 1941-docker-image-to-test, devel, master
- Children:
- 4bd88af
- Parents:
- fabac4a
- Location:
- tol_pkg/MulVarFunAprox
- Files:
-
- 4 added
- 5 edited
-
DesignC2.Monomial.Spheric.Disconnected.tol (added)
-
DesignC2.Monomial.Spheric.Grid.tol (modified) (11 diffs)
-
DesignC2.Monomial.tol (added)
-
DesignC2.Spheric.tol (added)
-
DesignC2.tol (modified) (10 diffs)
-
IndexedGrid.tol (added)
-
ModelC2.tol (modified) (3 diffs)
-
MulVarFunAprox.tol (modified) (1 diff)
-
test/test_0001/test.tol (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.Grid.tol
rfabac4a rb050cc43 9 9 Text _.autodoc.member.DesignC2.Monomial.Spheric.Grid = 10 10 ""; 11 Class @DesignC2.Monomial.Spheric.Grid : @DesignC2 11 Class @DesignC2.Monomial.Spheric.Grid : 12 @DesignC2.Monomial, 13 @DesignC2.Spheric 12 14 ////////////////////////////////////////////////////////////////////////////// 13 15 { … … 17 19 18 20 Real force.reference.near.nodal = True; 19 Real radius.margin = 0.05; 20 Real max.deg = 6; 21 Real grid.size = ?; 21 Real _.radius = ?; 22 Real _.grid.size = ?; 22 23 23 24 ////////////////////////////////////////////////////////////////////////////// … … 27 28 Set _.X.col.info = Copy(Empty); 28 29 Set _.range = Copy(Empty); 29 Real _.radius = ?;30 30 VMatrix _.grid = Rand(0,0,0,0); 31 31 Real _.G = ?; … … 37 37 38 38 //////////////////////////////////////////////////////////////////////////// 39 Static @DesignC2.Monomial.Spheric.Grid New(VMatrix X )39 Static @DesignC2.Monomial.Spheric.Grid New(VMatrix X, Real radius.margin) 40 40 //////////////////////////////////////////////////////////////////////////// 41 41 { … … 44 44 VMatrix _.X = X; 45 45 Real _.n = VColumns(X); 46 Real _.J = VRows(X) 46 Real _.J = VRows(X); 47 Real _.radius.margin = radius.margin 47 48 ]] 48 49 }; … … 61 62 VMatrix Xi = SubCol(X,[[i]]); 62 63 VMatrix Yi.min = Rand(m,1,_.range[i][1],_.range[i][1]); 63 VMatrix Ri = Round((Xi-Yi.min)/ grid.size);64 VMatrix Ri = Round((Xi-Yi.min)/_.grid.size); 64 65 VMatrix minRi = Rand(m,1,0,0); 65 66 VMatrix maxRi = Rand(m,1,cRi-1,cRi-1); … … 99 100 ]]}); 100 101 101 Real If(Or(IsUnknown( grid.size),grid.size<=0),102 Real If(Or(IsUnknown(_.grid.size),_.grid.size<=0), 102 103 { 103 104 Real length.max = SetMax(Traspose(Extract(_.X.col.info,4))[1]); 104 105 Real density = _.J/SetProd(Traspose(Extract(_.X.col.info,4))[1]); 105 grid.size := 0.5*(_.K0/density)^(1/_.n)106 }); 107 108 Real _.radius := (1+ radius.margin)*grid.size/Sqrt(2);106 _.grid.size := 0.5*(_.K0/density)^(1/_.n) 107 }); 108 109 Real _.radius := (1+_.radius.margin)*_.grid.size/Sqrt(2); 109 110 110 111 WriteLn("TRACE "+_MID+"2"); … … 112 113 { 113 114 Set info = _.X.col.info[i]; 114 Real points = 1+Ceil(info::length/ grid.size);115 Real exceed = (points-1)* grid.size-info::length;115 Real points = 1+Ceil(info::length/_.grid.size); 116 Real exceed = (points-1)*_.grid.size-info::length; 116 117 Matrix rnk = Tra(DifEq(1/(1-B),Rand(points,1,1,1))-Rand(points,1,1,1)); 117 Set MatSet((rnk* grid.size)+info::min-exceed/2)[1]118 Set MatSet((rnk*_.grid.size)+info::min-exceed/2)[1] 118 119 }); 119 120 WriteLn("TRACE "+_MID+"3"); … … 165 166 VMatrix Y.closest.grid.index = closest.grid.index(_.Y); 166 167 167 VMatrix Y.neighbour.pos = Min(Y.max, _.Y + Rand(_.H,_.n, grid.size,grid.size));168 VMatrix Y.neighbour.neg = Max(Y.min, _.Y - Rand(_.H,_.n, grid.size,grid.size));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)); 169 170 170 171 VMatrix _.Y.h := DifEq(1/(1-B),Rand(_.H,1,1,1)); … … 249 250 VMatrix _.Z.h := SubRow(Z.h.a | Z.h.b,MatQuery::SelectMatch(Z.interior)); 250 251 251 Real _.T := VRows(_.Z);252 Real _.Tz := VRows(_.Z); 252 253 Real _.Mz := Case( 253 254 prior.Z.order< 0, 0, 254 prior.Z.order==0, _.T, 255 prior.Z.order==1, _.T*(1+_.n), 256 prior.Z.order==2, _.T*(1+_.n*3/2+(1/2)*_.n^2)); 257 _.T 258 }; 259 260 ////////////////////////////////////////////////////////////////////////////// 261 Real build.base(Real void) 262 ////////////////////////////////////////////////////////////////////////////// 263 { 264 Text _MID = "[@DesignC2.Monomial.Spheric.Grid::build.base] "; 265 Real max.K = Min(_.K0, Floor(_.M/(_.H*min.contrast.surface))); 266 267 WriteLn("TRACE "+_MID+"1.1 K0="<<_.K0); 268 WriteLn("TRACE "+_MID+"1.2 M="<<_.M); 269 WriteLn("TRACE "+_MID+"1.3 H="<<_.H); 270 WriteLn("TRACE "+_MID+"1.4 SC="<<min.contrast.surface); 271 WriteLn("TRACE "+_MID+"1.5 max.K="<<max.K); 272 Real deg = 273 { 274 Real deg = 0; 275 Real K = 1; 276 Set pair.deg.K = [[ [[Copy(deg),Copy(K),1]] ]]; 277 While(And(deg<=max.deg, K<=max.K), 278 { 279 Real deg := deg+1; 280 Real k = Comb(_.n+deg-1,deg); 281 Real K := K+k; 282 WriteLn("TRACE "+_MID+"2."<<deg+" K="<<K); 283 Set Append(pair.deg.K, [[ [[Copy(deg),Copy(K),Copy(k)]] ]]) 284 }); 285 If(pair.deg.K[deg+1][2]<=max.K,deg,deg-1) 286 }; 287 288 WriteLn("TRACE "+_MID+"3"); 289 FunRn2R::@BaseC2 base.feasible = FunRn2R::@BaseC2.Monomial::New(_.n,deg); 290 Real _.K := base.feasible::_.K; 291 Set Append(_.base, [[base.feasible]]); 292 Real _.base.sel := Card(_.base); 293 _.K 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 294 259 }; 295 260 … … 305 270 Real _.K0 := base.max::_.K; 306 271 _build.model.design(void) 307 };308 309 ///////////////////////////////////////////////////////////////////////////310 Set get.X.weight(Set X.dif)311 ///////////////////////////////////////////////////////////////////////////312 {313 Real m = VRows(X.dif[1]);314 Set X.tao = For(1,_.X.maxNeighbour,VMatrix(Real i)315 {316 VMatrix X.norm2 = RPow(X.dif[i],2)*Rand(_.n,1,1,1);317 RPow(X.norm2,0.5)/_.radius318 });319 Set X.tao2 = For(1,_.X.maxNeighbour,VMatrix(Real i)320 {321 RPow(X.tao[i],2)322 });323 Set X.match = For(1,_.X.maxNeighbour,VMatrix(Real i)324 {325 LE(X.tao[i],Rand(m,1,1,1))326 });327 Set X.lambda = For(1,_.X.maxNeighbour,VMatrix(Real i)328 {329 RPow(X.match[i] $* (Rand(m,1,1,1)-X.tao2[i]),4)330 });331 VMatrix X.lambda.sum = SetSum(X.lambda);332 333 Set X.weight = For(1,_.X.maxNeighbour,VMatrix(Real i)334 {335 X.lambda[i] $/ X.lambda.sum336 })337 272 }; 338 273 -
tol_pkg/MulVarFunAprox/DesignC2.tol
rfabac4a rb050cc43 28 28 Real _.X.maxNeighbour = ?; 29 29 Real _.Z.maxNeighbour = ?; 30 VMatrix _.Y = Rand(0,0,0,0); 31 Real _.H = ?; 32 VMatrix _.Z = Rand(0,0,0,0); 33 VMatrix _.Z.h = Rand(0,0,0,0); 34 Real _.T = ?; 35 Real _.Mz = ?; 36 Real _.M = ?; 37 Real _.K = ?; 38 Real _.N = ?; 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; 39 43 Set _.base = Copy(Empty); 40 44 Real _.base.sel = ?; 41 45 VMatrix _.QX = Rand(0,0,0,0); 46 VMatrix _.QB = Rand(0,0,0,0); 42 47 VMatrix _.PZ = Rand(0,0,0,0); 43 48 VMatrix _.Q = Rand(0,0,0,0); 44 49 VMatrix _.LQ = Rand(0,0,0,0); 45 50 VMatrix _.Q.nonNullCellsByColumns = Rand(0,0,0,0); 46 51 ////////////////////////////////////////////////////////////////////////////// 47 52 NameBlock get.base(Real sel) … … 55 60 //Virtual methods 56 61 ////////////////////////////////////////////////////////////////////////////// 62 Real build.base(Real void); 57 63 Real build.reference.points(Real void); 58 Real build.stress.points(Real void); 59 Real build.base(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); 60 67 Set get.X.weight(Set X.dif); 61 68 VMatrix match.evaluable(VMatrix X); 62 VMatrix closest.reference.index(VMatrix X);63 69 64 70 ////////////////////////////////////////////////////////////////////////////// … … 70 76 Real build.reference.points(0); 71 77 Real build.stress.points(0); 72 Real _.M := _.J+_.Mz; 78 Real build.border.points(0); 79 Real _.M := _.J+_.Mb+_.Mz; 73 80 Real build.base(0); 74 81 Real _.N := _.H*_.K; … … 84 91 85 92 ////////////////////////////////////////////////////////////////////////////// 86 VMatrix build.obs(VMatrix X) 87 ////////////////////////////////////////////////////////////////////////////// 88 { 93 VMatrix build.obs(VMatrix X, VMatrix X.h) 94 ////////////////////////////////////////////////////////////////////////////// 95 { 96 Text _MID = "[@DesignC2::build.obs] "; 89 97 Real m = VRows(X); 90 VMatrix X.h=closest.reference.index(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"); 91 102 VMatrix X.j=DifEq(1/(1-B),Rand(m,1,1,1)); 103 WriteLn("TRACE "+_MID+"2"); 92 104 Set X.dif = For(1,_.X.maxNeighbour,VMatrix(Real i) 93 105 { … … 95 107 X - SubRow(_.Y,X.h.i) 96 108 }); 109 WriteLn("TRACE "+_MID+"2"); 97 110 Set X.weight = get.X.weight(X.dif); 111 WriteLn("TRACE "+_MID+"2"); 98 112 FunRn2R::@BaseC2 b = get.base(0); 113 WriteLn("TRACE "+_MID+"2"); 99 114 Matrix QX.triplet = Group("ConcatRows", For(1,_.X.maxNeighbour,Matrix(Real i) 100 115 { … … 114 129 })); 115 130 Pack(Convert(Triplet(QX.triplet, m, _.N),"Cholmod.R.Sparse")) 116 }; 117 ////////////////////////////////////////////////////////////////////////////// 118 VMatrix build.prior(Real void) 119 ////////////////////////////////////////////////////////////////////////////// 120 { 121 VMatrix Z.t=DifEq(1/(1-B),Rand(_.T,1,1,1)); 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)); 122 183 123 184 Set Z.dif = For(1,_.Z.maxNeighbour,VMatrix(Real i) … … 140 201 { 141 202 VMatrix row.idx = Z.t.i+cumRows; 142 Real cumRows := cumRows + _.T ;203 Real cumRows := cumRows + _.Tz; 143 204 VMat2Mat(row.idx | col.idx | pk) 144 205 } << If(prior.Z.order<1,Rand(0,3,0,0), … … 148 209 VMatrix pkd = b::partial.derivative.1.v(k,Zdif.i,d)*sign; 149 210 VMatrix row.idx = Z.t.i+cumRows; 150 Real cumRows := cumRows + _.T ;211 Real cumRows := cumRows + _.Tz; 151 212 VMat2Mat(row.idx | col.idx | pkd) 152 213 })))<< If(prior.Z.order<2,Rand(0,3,0,0), … … 158 219 VMatrix pkd2 = b::partial.derivative.2.v(k,Zdif.i,d1,d2)*sign; 159 220 VMatrix row.idx = Z.t.i+cumRows; 160 Real cumRows := cumRows + _.T ;221 Real cumRows := cumRows + _.Tz; 161 222 VMat2Mat(row.idx | col.idx | pkd2) 162 223 })) … … 165 226 })); 166 227 Convert(Triplet(PZ.triplet, _.Mz, _.N),"Cholmod.R.Sparse") 167 } ;228 })}; 168 229 ////////////////////////////////////////////////////////////////////////////// 169 230 Real build.decomposition(Real void) 170 231 ////////////////////////////////////////////////////////////////////////////// 171 232 { 172 VMatrix _.QX := build.obs(_.X); 173 VMatrix _.PZ := build.prior(0); 174 VMatrix Q.all = _.QX << _.PZ; 175 Set Q.active = MatQuery::SelectMatch( 176 NE(Tra(Rand(1,_.M,1,1)*Abs(Q.all)),Rand(_.N,1,0,0))); 177 VMatrix _.Q := SubCol(Q.all,Q.active); 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"); 178 253 VMatrix _.LQ := CholeskiFactor(_.Q,"XtX"); 254 WriteLn("TRACE "+_MID+"9"); 179 255 True 180 256 } -
tol_pkg/MulVarFunAprox/ModelC2.tol
rfabac4a rb050cc43 17 17 //Auxiliar members 18 18 ////////////////////////////////////////////////////////////////////////////// 19 VMatrix _.Fz = Rand(0,0,0,0); 20 VMatrix _.F = Rand(0,0,0,0); 21 VMatrix _.alpha = Rand(0,0,0,0); 22 VMatrix _.F.frc = Rand(0,0,0,0); 23 VMatrix _.F.cmp = Rand(0,0,0,0); 24 VMatrix _.E = Rand(0,0,0,0); 19 VMatrix _.F = Rand(0,1,0,0); 20 VMatrix _.Fz = Rand(0,1,0,0); 21 VMatrix _.Fb = Rand(0,1,0,0); 22 VMatrix _.alpha = Rand(0,1,0,0); 23 VMatrix _.F.frc = Rand(0,1,0,0); 24 VMatrix _.F.cmp = Rand(0,2,0,0); 25 VMatrix _.E = Rand(0,1,0,0); 25 26 Real _.E.avr = ?; 26 27 Real _.E.ste = ?; … … 44 45 { 45 46 VMatrix _.Fz := Rand(_.design::_.Mz,1,0,0); 46 VMatrix _.F := _.Fx<<_.Fz; 47 VMatrix _.Fb := Rand(_.design::_.Mb,1,0,0); 48 VMatrix _.F := _.Fx << _.Fb << _.Fz; 47 49 VMatrix _.alpha := CholeskiSolve(_.design::_.LQ,Tra(Tra(_.F)*_.design::_.Q),"PtLLtP"); 48 50 VMatrix _.F.frc := SubRow(_.design::_.Q*_.alpha,Range(1,_.design::_.J,1)); … … 63 65 VMatrix X.1 = SubRow(X,X.select); 64 66 Real m.1 = VRows(X.1); 65 VMatrix QX.1 = _.design::build.obs(X.1); 67 VMatrix X.1.h = _.design::closest.reference.index(X.1); 68 VMatrix QX.1 = _.design::build.obs(X.1, X.1.h); 66 69 VMatrix F.1 = QX.1*_.alpha; 67 70 Matrix triplet.1 = SetCol(X.select)|Rand(m.1,1,1,1)|VMat2Mat(F.1); -
tol_pkg/MulVarFunAprox/MulVarFunAprox.tol
rfabac4a rb050cc43 16 16 "vdebuen@tol-project.org" ]]; 17 17 18 #Embed "IndexedGrid.tol"; 18 19 #Embed "DesignC2.tol"; 19 20 #Embed "DesignC2.Monomial.tol"; 20 21 #Embed "DesignC2.Spheric.tol"; 21 22 #Embed "DesignC2.Monomial.Spheric.Grid.tol"; 23 #Embed "DesignC2.Monomial.Spheric.Disconnected.tol"; 22 24 #Embed "ModelC2.tol" 23 25 -
tol_pkg/MulVarFunAprox/test/test_0001/test.tol
rfabac4a rb050cc43 14 14 15 15 //Random seed setting 16 Real PutRandomSeed( 0);16 Real PutRandomSeed(1096004697); 17 17 Real rndSeed = GetRandomSeed(0); 18 18 WriteLn("Current Random Seed = "<<rndSeed); … … 33 33 })) 34 34 }; 35 36 MulVarFunAprox::@DesignC2 design = 35 VMatrix X = generate.nodal.points(15000,2); 36 /* * / 37 MulVarFunAprox::@DesignC2 design = { 37 38 MulVarFunAprox::@DesignC2.Monomial.Spheric.Grid::New( 38 generate.nodal.points(15000,2)); 39 40 Real design::force.reference.near.nodal := True; 41 Real design::radius.margin := 0.05; 39 X, radius.margin=0.05) }; 40 /* */ 41 MulVarFunAprox::@DesignC2 design = { 42 MulVarFunAprox::@DesignC2.Monomial.Spheric.Disconnected::New( 43 VMatrix X, 44 Real radius.margin = Sqrt(2)-1, 45 Real radius = .5, 46 Real do.build.border = False)}; 47 /* */ 48 //Real design::force.reference.near.nodal := True; 49 Real design::_.radius.margin := 0.05; 42 50 Real design::max.deg := 6; 43 51 … … 87 95 VMatrix X2 = { design::generate.internal.points( 88 96 Real m = design::_.J/2, 89 Real design::force.reference.near.nodal)};97 True)}; 90 98 91 99 Real b::upgrade.cache(X2,0);
Note: See TracChangeset
for help on using the changeset viewer.
