1 2printf( "Starting the npsol test...\n" ); 3 4assert = strip (function (t) 5{ 6 if (!test(t)) 7 { 8 message ("...failed.\a"); 9 exception (); 10 } 11}); 12 13# This code attempts to find the polygon with largest area among all 14# polygons with N sides and diameter no greater than one. For N odd 15# the answer is a regular polygon, but that's not the case when N is 16# even. (For quadrilaterals, though, the regular figure has the same 17# area as the non-regular solution.) 18 19# More information, and some optimization program benchmarks, may be 20# found at http://www-unix.mcs.anl.gov/~more/cops/bcops/polygon.html. 21 22# I stuck this here mainly because I thought it was cool, the npsol 23# performance seemed to compare pretty well with the benchmarked 24# codes, and because I'm a packrat. 25 26# The design variables are stored in a vector as pairs of distance 27# and angle. The final vertex is taken at the origin. 28 29dv2r = function (dv) 30{ 31 dv = vector (dv); 32 return dv [1:dv.ne:2]; 33}; 34 35dv2theta = function (dv) 36{ 37 dv = vector (dv); 38 return dv [2:dv.ne:2]; 39}; 40 41objf = function (dv) 42{ 43 local (r; theta; r1; r2); 44 45 r = dv2r (dv); 46 theta = dv2theta (dv); 47 48 r1 = r[1:r.ne-1]; 49 r2 = r[2:r.ne]; 50 51 return - 0.5 * sum (r1 @ r2 @ sin (diff (theta))); 52}; 53 54objgrd = function (dv) 55{ 56 local (r; theta; r1; r2; dt; st; ct; n; g); 57 58 dv = vector (dv); 59 r = dv2r (dv); 60 theta = dv2theta (dv); 61 62 r1 = r[1:r.ne-1]; 63 r2 = r[2:r.ne]; 64 dt = diff (theta); 65 st = sin (dt); 66 ct = cos (dt); 67 68 n = dv.ne; 69 g = fill (n; 0.0); 70 g[1:n-3:2] = r2 @ st; 71 g[3:n:2] += r1 @ st; 72 g[2:n-2:2] = - r1 @ r2 @ ct; 73 g[4:n:2] += r1 @ r2 @ ct; 74 75 return -g/2; 76}; 77 78conf = function (dv) 79{ 80 local (r; theta; n; nc; i; j; k; c); 81 82 r = dv2r (dv); 83 n = r.ne; 84 theta = dv2theta (dv); 85 86 nc = n*(n-1)/2 + (n-1); 87 c = fill (nc; 0.0); 88 89 k = 0; 90 for (i in seq(n-1)) 91 { 92 for (j in i+1:n) 93 { 94 c[k+=1] = 1.0 - r[i]^2 - r[j]^2 + 2.0*r[i]*r[j]*cos(theta[j]-theta[i]); 95 } 96 } 97 98 for (i in seq(n-1)) 99 { 100 c[k+=1] = theta[i+1] - theta[i]; 101 } 102 103 return c; 104}; 105 106congrd = function (dv) 107{ 108 local (r; theta; n; nc; i; j; k; c); 109 110 r = dv2r (dv); 111 n = r.ne; 112 theta = dv2theta (dv); 113 114 nc = n*(n-1)/2 + (n-1); 115 c = fill (nc,dv.ne; 0.0); 116 117 k = 0; 118 for (i in seq(n-1)) 119 { 120 for (j in i+1:n) 121 { 122 k += 1; 123 c[k;2*i-1] = 2.0 * (r[j]*cos(theta[j]-theta[i]) - r[i]); 124 c[k;2*j-1] = 2.0 * (r[i]*cos(theta[j]-theta[i]) - r[j]); 125 c[k;2*i] = 2.0*r[i]*r[j]*sin(theta[j]-theta[i]); 126 c[k;2*j] = -2.0*r[i]*r[j]*sin(theta[j]-theta[i]); 127 } 128 } 129 130 for (i in seq(n-1)) 131 { 132 k += 1; 133 c[k;2*i] = -1.0; 134 c[k;2*(i+1)] = 1.0; 135 } 136 137 return c; 138}; 139 140crd = function (dv) 141{ 142 local (r; theta); 143 144 r = dv2r (dv); 145 theta = dv2theta (dv); 146 147 return 0.0, label (r@sin(theta); r@cos(theta)), 0.0; 148}; 149 150sides = function (dv) 151{ 152 local (s; n); 153 154 dv = vector (dv); 155 n = dv.ne; 156 s = fill (n,2; 0.0); 157 158 s [1:n; 1] = 0.0; 159 s [1:n:2; 2] = 1.0; 160 s [2:n:2; 2] = acos (-1); 161 162 return s; 163}; 164 165bounds = function (dv) 166{ 167 local (b; n); 168 169 n = conf(dv).ne; 170 b = fill (n,2; 0.0); 171 172 b [;1] = 0.0; 173 b [;2] = 1e10; 174 175 return b; 176}; 177 178start = function (n) 179{ 180 local (i; pi; z; s); 181 182 i = sqrt (-1); 183 pi = acos (-1); 184 185 z = (exp (i*(2*pi*seq(n)/n-pi/2)) + i) / 2; 186 187 s = fill (2*n; 0.0); 188 s [1:2*n:2] = abs (z); 189 s [2:2*n:2] = atan2 (imag (z); real (z)); 190 191 return s[1:s.ne-2]; 192}; 193 194go = function (n) 195{ 196 local (s); 197 198 n = vector (n); 199 if (n.ne == 1) 200 { 201 s = start (n[1]); 202 elseif (n.ne%2) 203 message ("start vector must have even length"); 204 exception (); 205 elseif (n.ne >= 4) 206 s = n; 207 else 208 message ("start vector too short"); 209 exception (); 210 } 211 212 return npsol ({objf; objgrd}; s; 213 {side_constraints = sides(s); 214 nonlinear_constraints = {values = {conf; congrd}; 215 bounds = bounds(s); 216 }}; 217 {derivative_level = 3; verify_level = -1}); 218}; 219 220nudge = function (n) 221{ 222 local (i; pi; z; s); 223 224 i = sqrt (-1); 225 pi = acos (-1); 226 227 z = (exp (i*(2*pi*sort(rand(n))-pi/2)) + i) / 2; 228 229 s = fill (2*n; 0.0); 230 s [1:2*n:2] = abs (z); 231 s [2:2*n:2] = atan2 (imag (z); real (z)); 232 233 return s[1:s.ne-2]; 234}; 235 236tries = function (n; N) 237{ 238 local (i; r; rr; v); 239 240 r = go (n); 241 v = r.objective; 242 243 for (i in 1:N) 244 { 245 rr = go (nudge (n)); 246 if (rr.objective < v) 247 { 248 r = rr; 249 v = r.objective; 250 } 251 } 252 253 return r; 254}; 255 256poly = function (n) 257{ 258 local (r); 259 r = tries (n; 5); 260 plot (crd (adjust (r.solution)); crd (start (50)); 1); 261 return -r.objective; 262}; 263 264adjust = function (dv) 265{ 266 dv[2:dv.ne:2] += (acos(1-2*min(dv[1],1)^2)/2 - dv[2] + 267 acos(-1)-acos(1-2*min(dv[dv.ne-1],1)^2)/2 - dv[dv.ne]) / 2; 268 return dv; 269}; 270 271 272if (npsol) 273{ 274 275 npsol_test_fn = function( x ) 276 { 277 local( x1; x2 ); 278 x1 = scalar( x[1] ); x2 = scalar( x[2] ); 279 return x1^4 - 4*x1*x2 + x2^4; 280 }; 281 282 N = npsol( npsol_test_fn; 10,10; {}; { deriv = 0 } ); 283 284 assert (norm( N.solution - (1,1) ) < 1e-5); 285 286 assert (abs (tries(5; 10).objective + 0.6571638901) < 1e-7); 287 288 printf ("...passed.\n"); 289 290else 291 292 printf ("...skipped.\n"); 293} 294