print(`PABLO: A small Maple package`): print(`written by Aaron Robertson and Pablo Parrilo.`): print(`It accompanies an article by`): print(`Parrilo, Robertson, and Saracino`): print(``): print(`To verify the result in Section 3.1`): print(`use the procedure Run16().`): print(``): print(`To verify the results in Section 4`): print(`use the procedure FindMin()`): print(``): with(student): with(combinat): LargeRect16:=proc(r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16) local S: S:=(r1+r2+r3+r4+r5+r6+r7+r8)*(1/4-r5-r6-r7-r8)+(1/2-r1-r2-r3-r4-r5-r6-r7-r8)*(r5+r6+r7+r8): S:=S+(r9+r10+r11+r12+r13+r14+r15+r16)*(1/4-r9-r10-r11-r12): S:=S+(1/2-r9-r10-r11-r12-r13-r14-r15-r16)*(r9+r10+r11+r12): S: end: MidRect16:=proc(w,x,y,z,a,b): (w+x+y+z)*(1/8-a-b)+(1/4-w-x-y-z)*(a+b): end: SmallRect16:=proc(x,y,z): (x+y)*(1/16-z)+(1/8-x-y)*z: end: TinySquares16:=proc(x,y): x*(1/16-y)+y*(1/16-x): end: binarrys16:=proc(n) #outputs all binary sequences of length n local i,j,k,S,B,C,c: S:=[]: B:=[]: for k from 1 to n do B:=[op(B),[seq(0,i=1..k-1),1/16,seq(0,i=1..n-k)]]: od: for j from 2 to nops(B) do C:=choose(B,j): for i from 1 to nops(C) do c:=C[i]: S:=[op(S),add(c[k],k=1..j)]: od: od: [[seq(0,i=1..n)],op(B),op(S)]: end: Total16:=proc(r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16) local A: A:=LargeRect16(r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16): #large rects A:=A+MidRect16(r1,r2,r3,r4,r3,r4)+MidRect16(r5,r6,r7,r8,r9,r10): #medium rects A:=A+MidRect16(r9,r10,r11,r12,r7,r8)+MidRect16(r13,r14,r15,r16,r13,r14): # A:=A+SmallRect16(r1,r2,r2)+SmallRect16(r3,r4,r9)+SmallRect16(r5,r6,r4): #small rects A:=A+SmallRect16(r7,r8,r11)+SmallRect16(r9,r10,r6): # A:=A+SmallRect16(r11,r12,r13)+SmallRect16(r13,r14,r8): # A:=A+SmallRect16(r15,r16,r15): # #A:=A+r1*(1/16-r1)+r16*(1/16-r16)+1/256: # end triangles A:=A+TinySquares16(r1,r1)+TinySquares16(r16,r16)+2/1024: # A:=A+14/1024: #int triangles A:=A+TinySquares16(r3,r2)+TinySquares16(r5,r3)+TinySquares16(r7,r4): # A:=A+TinySquares16(r9,r5)+TinySquares16(r11,r6)+TinySquares16(r13,r7): # A:=A+TinySquares16(r15,r8): # A:=A+TinySquares16(r2,r9)+TinySquares16(r4,r10)+TinySquares16(r6,r11): # A:=A+TinySquares16(r8,r12)+TinySquares16(r10,r13)+TinySquares16(r12,r14): # A:=A+TinySquares16(r14,r15): # #A:=A+Rects12(r3,r4,r2)+Rects12(r5,r6,r3)+Rects12(r7,r8,r4): # #A:=A+Rects12(r9,r10,r5)+Rects12(r11,r12,r6): # #A:=A+Rects12(r1,r2,r7)+Rects12(r3,r4,r8)+Rects12(r5,r6,r9): # #A:=A+Rects12(r7,r8,r10)+Rects12(r9,r10,r11): # A: end: Run16:=proc() #if r1 is an endpoint, may assume r1=0 local Vold,tag,restrt,St,i,m,B,b,q,expr2,C,CC,flg,flg2,j,S,R,RR,T,V,Vars,Sols,expr,r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15,i16,f,Z,z,vv,s,flgcnt,flg4,RRR,fcheck,parm: b:=0: Sols:=[]: R:=[r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16]: restrt:=0: T:=[1,0,1/16]: #0 means no reds, 1/16 means all reds, 1 means keep as variable #St:=[1,1,1,3,3,3,3,1,1,3,3,3]: #restrt:=1: #tag:=3: St:=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]: for i1 from St[1] to 2 do for i2 from St[2] to 3 do for i3 from St[3] to 3 do for i4 from St[4] to 3 do for i5 from St[5] to 3 do for i6 from St[6] to 3 do for i7 from St[7] to 3 do for i8 from St[8] to 3 do for i9 from St[9] to 3 do for i10 from St[10] to 3 do for i11 from St[11] to 3 do for i12 from St[12] to 3 do for i13 from St[13] to 3 do for i14 from St[14] to 3 do for i15 from St[15] to 3 do for i16 from St[16] to 3 do if restrt=1 then St:=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]: restrt:=0: fi: flg4:=0: V:=[T[i1],T[i2],T[i3],T[i4],T[i5],T[i6],T[i7],T[i8],T[i9],T[i10],T[i11],T[i12],T[i13],T[i14],T[i15],T[i16]]: R:=[r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16]: for j from 1 to 16 do if V[j]<>1/16 then R[j]:=R[j]*V[j] else R[j]:=1/16: fi: od: #below used to track restarts if V[tag]=0 and Vold[tag]=1 then print(R): tag:=tag-1: fi: Vars:={op(R)} minus {0,1/16}: expr:=Total16(op(R)): flg:=0: flg2:=0: C:={0} union extrema(expr,{},Vars,'s'): flg:=0: S:=[op(op(s))]: #print(sis,S,R): if C={0} then # no critical point, so max at endpoint to be checked further in algorithm flg:=1: else # check if critical point is in bounds (necessarily only one) flgcnt:=0: for j from 1 to nops(S) do # if all numeric then can compare subsequent critical points to # one that we know is valid if is(rhs(S[j]),numeric)=true then flgcnt:=flgcnt+1: if is(rhs(S[j])<0)=true or is(rhs(S[j])>1/16)=true then j:=nops(S): flg:=1: fi: else flg2:=1: #print(flag2,R): fi: od: r1:='r1': r2:='r2': r3:='r3': r4:='r4': r5:='r5': r6:='r6': r7:='r7': r8:='r8': r9:='r9': r10:='r10': r11:='r11': r12:='r12': r13:='r13': r14:='r14': r15:='r15': r16:='r16': fi: if flg2=1 then # solution is parametric RR:=[r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16]: # need to check if solution is possibly in bounds (flg4=0 if in bounds) flg4:=0: for m from 1 to 16 do if R[m] = 0 or R[m] = 1/16 then RR[m]:=R[m]: else RRR:={}: for q from 1 to nops(S) do parm:=rhs(S[q]): if lhs(S[q])=R[m] then RR[m]:=parm: #print(HEREisPARM,parm): if is(parm,numeric)=false then if has(parm,r1) then RRR:=RRR union {r1}:fi: if has(parm,r2) then RRR:=RRR union {r2}:fi: if has(parm,r3) then RRR:=RRR union {r3}:fi: if has(parm,r4) then RRR:=RRR union {r4}:fi: if has(parm,r5) then RRR:=RRR union {r5}:fi: if has(parm,r6) then RRR:=RRR union {r6}:fi: if has(parm,r7) then RRR:=RRR union {r7}:fi: if has(parm,r8) then RRR:=RRR union {r8}:fi: if has(parm,r9) then RRR:=RRR union {r9}:fi: if has(parm,r10) then RRR:=RRR union {r10}:fi: if has(parm,r11) then RRR:=RRR union {r11}:fi: if has(parm,r12) then RRR:=RRR union {r12}:fi: if has(parm,r13) then RRR:=RRR union {r13}:fi: if has(parm,r14) then RRR:=RRR union {r14}:fi: if has(parm,r15) then RRR:=RRR union {r15}:fi: if has(parm,r16) then RRR:=RRR union {r16}:fi: fcheck:=unapply(parm,op(RRR)): if nops(RRR)=1 then if (fcheck(0)<0 and fcheck(1/16)<0) then flg4:=1: q:=nops(S): m:=16: elif (fcheck(0)>1/16 and fcheck(1/16)>1/16) then flg4:=1: q:=nops(S): m:=16: fi: fi: # if nops(RRR)>1 let it go, since there are not many like this fi: #if is(parm,numeric)=false q:=nops(S): fi: #if lhs(S[q])=R[m] then RR[m]:=parm: od: # q fi: # if R[m] = 0 or R[m] = 1/16 od: # m # if solution is in bounds, then add to set of critical points if flg4=0 then expr2:=expand(Total16(op(RR))): C:=C union {expr2}: #expr2 must be numerical fi: # flg4=0 (in bounds) fi: # if flg2=1 (parametric) if flg=0 and flg4=0 then if max(op(C)) > b then #.2716 print(1/2*(3/8-evalf(max(op(C)))),1/2*(3/8-max(op(C))),S,R): fi: if b Heaviside(1-x+2*y)*Heaviside(1+x-2*y); q := (x,y) -> p(x,y)+p(y,x); inti := (a,b,c,d) -> int(int(q(x,y),x=a..b),y=c..d): np := nops(part)-1; tot := 0 : for j from 1 to np do for k from 1 to np do tot := tot + (-1)^(j+k) *inti(part[j],part[j+1],part[k],part[k+1]); od: od: return tot; end proc: # This is the approximate partition obtained numerically # part := [-1,-.898,-.8762,-.7738,-.6386,-.4234,0,.4234,.6386,.7738,.8762,.898,1]: # Let's get the integral as a function of parameters # Here the structure is from the numerically obtained partitions FindMin:=proc() local p1,p2,p3,p4,p5,LOC,k,A,eqs,bmat,popt,part,costpart: part := [-1,-p1,-p2,-p3,-p4,-p5,0,p5,p4,p3,p2,p1,1]; costpart := computecost(part): # Now, let's get a Taylor expansion, around the points computed numerically... LOC:=mtaylor(costpart,{p1=8980/10000,p2=8762/10000,p3=7738/10000,p4=6386/10000,p5=4234/10000},4); # Great, it's a quadratic function # Let's compute the minimum, finding the gradient first... eqs:=map(k->diff(LOC,k),[p1,p2,p3,p4,p5]); # Let's get the matrix, and verify it's PSD.. A:=linalg[genmatrix](eqs,[p1,p2,p3,p4,p5],'bmat'); print(`The eigenvalues are:`): print(linalg[eigenvalues](evalf(evalm(A)))); print(`So it is PSD.`): # Here's the global minimum.... print(``): print(`Now, here's the minimum:`): popt:=linalg[linsolve](A,bmat); end: