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