Subject | : Re: Packages/optimize.n 1.42 |
Date | : 2010/04/16(Fri) 11:24:27 |
Contributor | : Akio Morita |
DownhillSimplex[]関数の能力が、デグレードしているように見えます 目標関数として正定値な二次形式(下に凸な放物面)を使った場合の例は以下の通り N = 1と 2のケースで明らかな劣化を確認しました (* from Packages/optimize.n 1.42 *) SimplexNEW=Class[{},{},{CurrentSimplex,Status="Stop", limit,n,lo,tol,per,minq,exp,cont,quick,exl,cl,rr={1,-1}}, Amoebar[p_,f_,r_,kr_]:=Module[ {pr,fac1=(1-r)/(Length[p]-1),fr,tp=Take[p[[,2]],rr]}, ! {pr,fac1=(1-r)/(L),fr,tp=Take[p[[,2]],rr]}, tp=(Plus@@tp)/Length[tp]; pr=(1-r)*tp+p[[kr,2]]*r; {pr,fr=f[pr]+ (Plus@@MapThread[(#<#2||#>#3)&,{pr,limit[[1]],limit[[2]]}])* Abs[p[[-1,1]]-p[[1,1]]]*10, If[fr<p[[kr,1]], Sort[ReplacePart[p,{fr,pr},kr]], p]}]; Amoeba[p_,f_]:=Module[ {pr,pr1,pl,fr,fr1,prr,n,kr,fr0,pr0,pt}, Do[ Status=StandardForm["Reflect "//kr]; {pt,fr,pr}=Amoebar[p,f,-1,-kr]; (* Print["dh1 ",{fr,p[[exl,1]]}]; *) If[kr==1,fr0=fr;pr0=pr]; If[fr <= p[[exl,1]], Status="Extend"; {pt,fr1,pr1}=Amoebar[p,f,-exp,-kr]; Return[If[fr1 <= fr, pr1, pr]]], (* If[fr1<=pr[[-1,1]], Sort[ReplacePart[pr,{fr1,pt},-1]], pr]]]], *) {kr,cont}]; If[cont<>1,fr=fr0;pr=pr0]; If[Length[p]>2 && fr < p[[-1-cl,1]], pr, If[fr < p[[-1,1]], Status="Contract+", Status="Contract-"]; {pt,fr1,pr1}=Amoebar[pr,f,0.5,-1]; If[fr1 < pr[[-1-cl,1]], pr1, pl=pr[[1,2]]; prr=Thread[(pl+Thread[Drop[pr,1][[,2]]])/2]; prr=(MapThread[(If[#2<#<#3,(#>#4)*(#3-#2)+#2,#])&, {#,pl-minq,pl+minq,pl}]&/@prr); Do[ Status=StandardForm["Shrink "//n]; fv=f[prr[[n-1]]]; pr[[n]]={fv,prr[[n-1]]}; If[quick && (fv<pr[[n-1,1]] || fv>pr[[-1,1]]), Break[]], {n,Length[pr],2,-1}]; If[per<>1, pr[[1]]={pr[[1,1]]*per+f[pl]*(1-per),pl}]; Sort[pr]]]]; DownhillSimplex[p0_,f_,opt___]:=Module[ {def={VariableRange:>{Table[Max[],{Length[p0[[1,2]]]}], Table[Min[],{Length[p0[[1,2]]]}]}, MaxIteration:>Max[100,10*Length[p0]], ExtendLevel->1,ContractLevel->1,ReflectRange->{1,-1}, Output->6,MinimumSize->0,ExpansionRatio->2,TryReflect->1, Tolerance->1e-6,Persistence->1,QuickShrink->0}, lo,p=Sort[p0],nc=Floor[(PageWidth-3)/22]}, {limit,n,lo,tol,per,minq,exp,cont,quick,exl,cl,rr}= {VariableRange,MaxIteration,Output, Tolerance,Persistence,MinimumSize,ExpansionRatio,TryReflect,QuickShrink, ExtendLevel,ContractLevel,ReflectRange}/.{opt}/.def; CurrentSimplex=p; exl=Restrict[exl,1,Length[p]]; cont=Restrict[cont,1,Length[p]]; Check[ Catch[ Do[CurrentSimplex=p=Amoeba[p,f]; If[lo > 0, StandardForm[ $FORM='10.6'; Write[lo,Drop[p,{nc+1,Length[p0]-nc}][[,1]]]]]; If[(p[[-1,1]]-p[[1,1]])/(1e-100+Abs[p[[-1,1]]]+Abs[p[[1,1]]]) < tol, Break[]],{n}]],]; Status="Stop"; p]; ]; (* from Packages/optimize.n 1.41 *) SimplexOLD=Class[{},{},{CurrentSimplex,Status="Stop"}, Amoebar[p_,f_,r_,limit_,kr_]:=Module[ {pr,fac1=(1-r)/(Length[p]-1),fr}, pr=Plus@@p[[,2]]*fac1+p[[kr,2]]*(r-fac1); {pr,fr=f[pr]+ (Plus@@MapThread[(#<#2||#>#3)&,{pr,limit[[1]],limit[[2]]}])* Abs[p[[-1,1]]-p[[1,1]]]*10, If[fr<p[[kr,1]], Sort[ReplacePart[p,{fr,pr},kr]], p]}]; Amoeba[p_,f_,limit_,per_,minq_,exp_,cont_,quick_,exl_]:=Module[ {pr,pr1,pl,fr,fr1,prr,n,kr,fr0,pr0,pt}, Do[ Status=StandardForm["Reflect "//kr]; {pt,fr,pr}=Amoebar[p,f,-1,limit,-kr]; (* Print["dh1 ",{fr,p[[exl,1]]}]; *) If[kr==1,fr0=fr;pr0=pr]; If[fr <= p[[exl,1]], Status="Extend"; {pt,fr1,pr1}=Amoebar[p,f,-exp,limit,-kr]; Return[If[fr1 <= fr, pr1, pr]]], (* If[fr1<=pr[[-1,1]], Sort[ReplacePart[pr,{fr1,pt},-1]], pr]]]], *) {kr,cont}]; If[cont<>1,fr=fr0;pr=pr0]; If[Length[p]>2 && fr < p[[-2,1]], pr, If[fr < p[[-1,1]], Status="Contract+", Status="Contract-"]; {pt,fr1,pr1}=Amoebar[pr,f,0.5,limit,-1]; If[fr1 < pr[[-1,1]], pr1, pl=pr[[1,2]]; prr=Thread[(pl+Thread[Drop[pr,1][[,2]]])/2]; prr=(MapThread[(If[#2<#<#3,(#>#4)*(#3-#2)+#2,#])&, {#,pl-minq,pl+minq,pl}]&/@prr); Do[ Status=StandardForm["Shrink "//n]; fv=f[prr[[n-1]]]; pr[[n]]={fv,prr[[n-1]]}; If[quick && (fv<pr[[n-1,1]] || fv>pr[[-1,1]]), Break[]], {n,Length[pr],2,-1}]; If[per<>1, pr[[1]]={pr[[1,1]]*per+f[pl]*(1-per),pl}]; Sort[pr]]]]; DownhillSimplex[p0_,f_,opt___]:=Module[ {def={VariableRange:>{Table[Max[],{Length[p0[[1,2]]]}], Table[Min[],{Length[p0[[1,2]]]}]}, MaxIteration:>Max[100,10*Length[p0]], ExtendLevel->1, Output->6,MinimumSize->0,ExpansionRatio->2,TryReflect->1, Tolerance->1e-6,Persistence->1,QuickShrink->0}, limit,n,tol,per,minq,exp,cont,quick,exl, lo,p=Sort[p0],nc=Floor[(PageWidth-3)/22]}, {limit,n,lo,tol,per,minq,exp,cont,quick,exl}= {VariableRange,MaxIteration,Output, Tolerance,Persistence,MinimumSize,ExpansionRatio,TryReflect,QuickShrink, ExtendLevel}/.{opt}/.def; CurrentSimplex=p; exl=Restrict[exl,1,Length[p]]; cont=Restrict[cont,1,Length[p]]; Check[ Catch[ Do[CurrentSimplex=p=Amoeba[p,f,limit,per,minq,exp,cont,quick,exl]; If[lo > 0, StandardForm[ $FORM='10.6'; Write[lo,Drop[p,{nc+1,Length[p0]-nc}][[,1]]]]]; If[(p[[-1,1]]-p[[1,1]])/(1e-100+Abs[p[[-1,1]]]+Abs[p[[1,1]]]) < tol, Break[]],{n}]],]; Status="Stop"; p]; ]; The$SimplexOLD = SimplexOLD[]; The$SimplexNEW = SimplexNEW[]; func = With[{l = Range[Length[#]]}, Plus@@(Reverse[l] * (# + l)^2)]&; f1 = func[{#1}]&; f2 = func[{#1, #2}]&; N = 1; size = 1e-1; x0 = Table[0, {N}]; p = {func[x0 + #], x0 + #}&/@Append[size * IdentityMatrix[N], Table[0, {N}]]; l = Transpose[Table[2 * N * {-1, 1}, {N}]]; opt = {VariableRange->l, Tolerance->1e-12, MaxIteration->1000}; rO = The$SimplexOLD@DownhillSimplex[p, func, Null@@opt]; rN = The$SimplexNEW@DownhillSimplex[p, func, Null@@opt]; Print["\nOLD ->\n", rO]; Print["\nNEW ->\n", rN]; Switch[N, 1, Plot[f1[z], {z, -2, 2}]; TkWait[], 2, ListContourPlot[Table[func[{i, j}], {i, -4, 4, 0.1}, {j, -4, 4, 0.1}], MeshRange->{{-4, 4}, {-4, 4}}, Contours->30]; TkWait[], _, (* Nothing *)];