Free cookie consent management tool by TermsFeed Policy Generator

source: branches/MathNetNumerics-Exploration-2789/Testing/Form.cs @ 15428

Last change on this file since 15428 was 15348, checked in by gkronber, 7 years ago

#2789 worked on NLR with constraints
File size: 13.5 KB
Line 
1using System;
2using System.Linq;
3using System.Windows.Forms;
4using System.Windows.Forms.DataVisualization.Charting;
5using MathNet.Numerics.Distributions;
6
7namespace Testing {
8  public partial class TestForm : Form {
9    public TestForm() {
10      InitializeComponent();
11
12      //Interpolate();
13      MinAdaptiveGradientSampling();
14    }
15
16    // E P p v T Tmax
17    double[,] data = new double[,]
18    {
19        {0.8368, 1.3357, 0.8023, 12.3839, 53.7792, 201.9808},
20        {0.5197, 1.4065, 0.5648, 17.5391, 53.0914, 162.0023},
21        {0.8406, 1.4284, 0.5660, 17.5387, 52.9766, 205.0202},
22        {1.1622, 1.4325, 0.5670, 17.5386, 53.2497, 234.6500},
23        {0.5208, 0.8740, 0.3042, 17.5388, 52.7885, 148.3409},
24        {0.8411, 0.8848, 0.3054, 17.5390, 52.7723, 172.0350},
25        {1.1628, 0.8866, 0.3060, 17.5396, 53.1371, 186.4282},
26        {0.8472, 1.9431, 0.8010, 17.5390, 52.8836, 221.6862},
27        {1.1679, 1.9221, 0.8020, 17.5391, 53.0837, 258.8506},
28        {1.5070, 1.9104, 0.8031, 17.5384, 53.3303, 286.1709},
29        {0.5085, 1.3330, 0.8003, 12.3837, 53.2506, 156.4186},
30        {0.8374, 1.3556, 0.8015, 12.3837, 53.1234, 200.7889},
31        {1.4964, 1.3682, 0.8036, 12.3838, 53.6585, 252.7940},
32        {0.5039, 0.7592, 0.5658, 9.5036, 53.6265, 149.1639},
33        {0.8291, 0.7754, 0.5671, 9.5032, 53.5168, 182.6419},
34        {1.4895, 0.7920, 0.5689, 9.5036, 53.9511, 219.4112},
35        {0.8373, 1.3749, 0.8015, 12.3835, 53.5531, 201.8539},
36        {0.5038, 0.7309, 0.8014, 6.7081, 53.7301, 152.8836},
37        {0.8278, 0.7412, 0.8031, 6.7077, 53.8957, 192.8725},
38        {1.0551, 0.7459, 0.8041, 6.7081, 54.1965, 213.8598},
39        {0.8374, 1.3772, 0.8011, 12.3836, 53.4341, 202.3147},
40        {0.5108, 1.2544, 1.4784, 6.7082, 53.7971, 162.4619},
41        {0.8344, 1.2516, 1.4843, 6.7078, 53.9234, 215.3835},
42        {1.0614, 1.2461, 1.4868, 6.7079, 54.0830, 244.6310},
43        {0.8370, 1.3392, 0.8009, 12.3842, 53.4277, 200.7755},
44        {0.5167, 1.6944, 2.0975, 6.7078, 53.8257, 166.3573},
45        {0.8405, 1.6901, 2.1060, 6.7078, 53.8974, 224.3230},
46        {1.0674, 1.6778, 2.1095, 6.7082, 54.1107, 257.8381},
47        {0.8371, 1.3172, 0.8010, 12.3839, 53.4347, 199.6043},
48        {0.8329, 1.2290, 0.2155, 24.8521, 52.6891, 175.7116},
49        {1.1958, 0.9784, 0.2163, 24.8525, 52.8742, 197.7365},
50        {1.5166, 0.9721, 0.2168, 24.8520, 53.0470, 213.4699},
51        {0.8383, 1.4018, 0.8009, 12.3835, 53.5269, 201.9488},
52        {0.8453, 1.5948, 0.3972, 24.8515, 52.8218, 205.8914},
53        {1.1740, 1.5958, 0.3981, 24.8518, 52.7543, 236.6865},
54        {1.4963, 1.6219, 0.3988, 24.8524, 53.1004, 260.2091},
55        {0.8383, 1.4094, 0.8008, 12.3838, 53.6704, 202.8565},
56        {0.8538, 2.1021, 0.5648, 24.8523, 52.6869, 221.3082},
57        {1.1752, 2.0947, 0.5656, 24.8519, 52.8542, 260.5521},
58        {1.4978, 2.0919, 0.5661, 24.8529, 53.2062, 291.1189},
59        {0.8383, 1.3991, 0.8004, 12.3842, 53.5893, 201.8672},
60        {0.8423, 1.8295, 1.4793, 9.5039, 53.7088, 219.3284},
61        {1.1670, 1.8139, 1.4836, 9.5038, 53.8615, 262.5016},
62        {1.5025, 1.7922, 1.4866, 9.5035, 54.1248, 297.1750},
63        {0.8377, 1.3349, 0.8005, 12.3838, 53.6432, 198.6240},
64        {0.8371, 1.2622, 0.8137, 12.3844, 101.9483, 240.1088},
65        {0.5196, 1.3274, 0.5751, 17.5397, 102.0049, 202.1158},
66        {0.8409, 1.3472, 0.5760, 17.5391, 100.8094, 239.7598},
67        {1.1625, 1.3534, 0.5763, 17.5390, 100.9805, 265.2706},
68        {0.5225, 0.8233, 0.3109, 17.5394, 100.0307, 186.0851},
69        {0.8433, 0.8225, 0.3116, 17.5395, 99.9612, 207.5332},
70        {1.1640, 0.8288, 0.3120, 17.5393, 100.3567, 222.7550},
71        {0.8473, 1.8966, 0.8123, 17.5390, 100.1643, 254.6133},
72        {1.1681, 1.8691, 0.8132, 17.5392, 100.3375, 289.9712},
73        {1.5075, 1.8382, 0.8136, 17.5392, 100.7590, 315.6024},
74        {0.5081, 1.2802, 0.8109, 12.3841, 100.6839, 196.4331},
75        {0.8371, 1.2978, 0.8123, 12.3839, 100.6797, 238.7059},
76        {1.4964, 1.2975, 0.8145, 12.3835, 101.5454, 286.6940},
77        {0.5038, 0.7123, 0.5749, 9.5036, 101.3295, 185.5991},
78        {0.8289, 0.7297, 0.5757, 9.5036, 101.3065, 214.1675},
79        {1.4891, 0.7365, 0.5764, 9.5035, 102.1764, 244.6417},
80        {0.8374, 1.3175, 0.8121, 12.3840, 101.3808, 240.0907},
81        {0.5037, 0.7071, 0.8116, 6.7080, 102.1794, 191.0549},
82        {0.8278, 0.7135, 0.8142, 6.7079, 102.5190, 225.8846},
83        {1.0549, 0.7150, 0.8156, 6.7079, 102.9568, 243.3557},
84        {0.8371, 1.3160, 0.8119, 12.3841, 101.1813, 240.1163},
85        {0.5110, 1.2456, 1.5031, 6.7075, 102.2212, 200.6899},
86        {0.8345, 1.2408, 1.5095, 6.7075, 102.5058, 249.4029},
87        {1.0615, 1.2257, 1.5114, 6.7077, 103.0208, 276.3802},
88        {0.8372, 1.2756, 0.8115, 12.3839, 101.2421, 239.1648},
89        {0.5168, 1.6720, 2.1301, 6.7078, 102.2741, 204.8973},
90        {0.8403, 1.6587, 2.1390, 6.7079, 102.5263, 258.4755},
91        {1.0668, 1.6367, 2.1399, 6.7078, 102.9995, 291.2107},
92        {0.8368, 1.2379, 0.8111, 12.3839, 101.3221, 238.3200},
93        {0.8191, 1.3475, 0.2200, 24.8519, 99.4185, 213.5561},
94        {1.1952, 0.9615, 0.2205, 24.8520, 99.7065, 237.2100},
95        {1.5198, 0.9392, 0.2206, 24.8526, 100.2240, 254.1504},
96        {0.8373, 1.3301, 0.8111, 12.3842, 101.4987, 239.6880},
97        {0.8401, 1.5435, 0.4047, 24.8529, 99.6446, 242.0860},
98        {1.1737, 1.5140, 0.4051, 24.8520, 99.8763, 270.9338},
99        {1.4966, 1.5101, 0.4053, 24.8527, 100.3614, 295.7137},
100        {0.8370, 1.3303, 0.8109, 12.3839, 101.5883, 239.4128},
101        {0.8515, 2.0246, 0.5744, 24.8527, 99.7521, 256.5944},
102        {1.1740, 2.0306, 0.5749, 24.8528, 99.9245, 293.3329},
103        {1.4970, 2.0055, 0.5749, 24.8522, 100.4382, 320.3913},
104        {0.8373, 1.3181, 0.8106, 12.3846, 101.6315, 238.0053},
105        {0.8417, 1.7571, 1.5057, 9.5039, 102.0292, 251.9787},
106        {1.1665, 1.7413, 1.5092, 9.5036, 102.2874, 292.1148},
107        {1.5021, 1.7017, 1.5109, 9.5040, 102.7524, 325.0641},
108        {0.8365, 1.2463, 0.8105, 12.3839, 101.9477, 236.4320},
109        {0.8370, 1.2600, 0.8101, 12.3843, 104.2213, 240.013},
110    };
111
112
113
114    //private void MinBleic() {
115    //
116    //  double[] w = new double[14];
117    //  double[] optW;
118    //
119    //  double[,] constraints = new double[,]
120    //  {
121    //  // 0  1  2  3  4  5  6  7  8  9 10 11 12 13  y
122    //    {0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0 }
123    //  };
124    //  int[] constraintType = new int[]
125    //  {
126    //    -1       
127    //  };
128    //
129    //  alglib.minbleicstate state;
130    //  alglib.minbleicreport report;
131    //
132    //  alglib.minbleiccreate(w, out state);
133    //  alglib.minbleicsetcond(state, 0, 1E-10, 0.0, 0);
134    //  alglib.minbleicsetlc(state, constraints, constraintType);
135    //  alglib.minbleicsetgradientcheck(state, 0.001);
136    //  alglib.minbleicsetxrep(state, true);
137    //  alglib.minbleicoptimize(state, Grad, Rep, null);
138    //  alglib.minbleicresults(state, out optW, out report);
139    //
140    //}
141
142    private AutoDiff.Term f;
143
144    private void MinAdaptiveGradientSampling() {
145
146      double[] w = new double[14];
147      for (int i = 0; i < w.Length; i++) w[i] = 1.0;
148      double[] optW;
149
150      double epsx = 0.00001;
151      double radius = 100;
152      double rho = 10.0;
153      int maxits = 0;
154      int numEqConstraints = 0;
155      int numIneqConstraints = 6;
156
157
158      alglib.minnsstate state;
159      alglib.minnsreport report;
160
161      alglib.minnscreate(w, out state);
162      alglib.minnssetnlc(state, numEqConstraints, numIneqConstraints);
163      alglib.minnssetalgoags(state, radius, rho);
164      alglib.minnssetcond(state, 0, 0);
165      alglib.minnssetxrep(state, true);
166      alglib.minnsoptimize(state, Jac, Rep, null);
167      alglib.minnsresults(state, out optW, out report);
168
169    }
170
171
172    private void Rep(double[] w, double func, object o) {
173      Console.WriteLine(Math.Sqrt(func / data.GetLength(0)));
174    }
175
176
177    private void Jac(double[] w, double[] fi, double[,] jac, object o) {
178      Array.Clear(fi, 0, fi.Length);
179      Array.Clear(jac, 0, jac.Length);
180
181      // objective function
182      double f = 0.0;
183      double[] fgrad = new double[w.Length];
184      Grad(w, ref f, fgrad, o);
185      fi[0] = f;
186      Copy(fgrad, jac, 0);
187
188      // d Tmax / d Tmin > 0 for all valid inputs
189      fi[1] = -dTmaxOverdE(3, 4, w, fgrad);
190      Copy(fgrad, jac, 1);
191      fi[2] = -dTmaxOverdE(3, 0.1, w, fgrad);
192      Copy(fgrad, jac, 2);
193      fi[3] = -dTmaxOverdE(0.1, 4, w, fgrad);
194      Copy(fgrad, jac, 3);
195      fi[4] = -dTmaxOverdE(0.1, 0.1, w, fgrad);
196      Copy(fgrad, jac, 4);
197
198      // d Tmax / d p > 0 for all valid inputs
199      fi[5] = -dTmaxOverdp(3, 3, 0.1, 50, w, fgrad);
200      Copy(fgrad, jac, 5);
201      fi[6] = -dTmaxOverdp(3, 0.5, 0.1, 50, w, fgrad);
202      Copy(fgrad, jac, 6);
203
204      // d Tmax / d v > 0 for all valid inputs
205      // d Tmax / d E > 0 for all valid inputs
206
207    }
208
209    private void Copy(double[] from, double[,] to, int col)
210    {
211      for (int i = 0; i < from.Length; i++) {
212        to[col, i] = from[i];
213      }
214    }
215
216
217    private void Grad(double[] w, ref double func, double[] grad, object o) {
218      //
219      double ssq = 0.0;
220      for (int i = 0; i < grad.Length; i++) grad[i] = 0.0;
221      var fgrad = new double[grad.Length];
222      for (int r = 0; r < data.GetLength(0); r++) {
223        var pred = Tmax(data[r, 0], data[r, 2], data[r, 3], data[r, 4], w, fgrad);
224        var target = data[r, 5];
225        var err = pred - target;
226        ssq += err * err;
227
228        for (int i = 0; i < grad.Length; i++) {
229          grad[i] += -2.0 * target * fgrad[i];
230          grad[i] += 2.0 * pred * fgrad[i];
231        }
232      }
233
234      func = ssq;
235    }
236
237    private double Tmax(double E, double p, double v, double Tmin, double[] w, double[] grad) {
238      grad[0] = Tmin;
239      grad[1] = (v);
240      grad[2] = (p);
241      grad[3] = (E);
242      grad[4] = (Math.Log((p)) * v);
243      grad[5] = (Math.Log((p)) * E);
244      grad[6] = (v * p);
245      grad[7] = (Tmin * p);
246      grad[8] = (E * p);
247      grad[9] = (v * E);
248      grad[10] = (Tmin * E);
249      grad[11] = (E * E);
250      grad[12] = (E * 1 / p);
251      grad[13] = (1 / p);
252
253      return
254        w[0] * (Tmin) +
255        w[1] * (v) +
256        w[2] * (p) +
257        w[3] * (E) +
258
259        w[4] * (Math.Log((p)) * v) +
260        w[5] * (Math.Log((p)) * E) +
261
262        w[6] * (v * p) +
263        w[7] * (Tmin * p) +
264        w[8] * (E * p) +
265
266        w[9] * (v * E) +
267        w[10] * (Tmin * E) +
268        w[11] * (E * E) +
269
270        w[12] * (E * 1 / p) +
271
272        w[13] * (1 / p);
273    }
274
275    private double dTmaxOverdp(double E, double p, double v, double Tmin, double[] w, double[] grad) {
276      grad[0] = 0;
277      grad[1] = 0;
278      grad[2] = 1;
279      grad[3] = 0;
280      grad[4] = (v * 1/p);
281      grad[5] = (E * 1/p);
282      grad[6] = v ;
283      grad[7] = Tmin;
284      grad[8] = E;
285      grad[9] = 0;
286      grad[10] = 0;
287      grad[11] = 0;
288      grad[12] = -(E * 1 / (p*p));
289      grad[13] = -(1 / (p*p));
290
291      return
292        w[2] +
293
294        w[4] * (v * 1/p) +
295        w[5] * (E * 1/p) +
296
297        w[6] * (v ) +
298        w[7] * (Tmin ) +
299        w[8] * (E) +
300
301        w[12] * -(E * 1 / (p*p)) +
302
303        w[13] * -(1 / (p*p));
304    }
305
306    private double dTmaxOverdE(double E, double p, double[] w, double[] grad) {
307      grad[0] = 1;
308      grad[1] = 0;
309      grad[2] = 0;
310      grad[3] = 0;
311      grad[4] = 0;
312      grad[5] = 0;
313      grad[6] = 0;
314      grad[7] = p;
315      grad[8] = 0;
316      grad[9] = 0;
317      grad[10] = E;
318      grad[11] = 0;
319      grad[12] = 0;
320      grad[13] = 0;
321
322      return
323        w[0] +
324        w[7] * (p) +
325        w[10] * (E);
326    }
327
328    private void Interpolate() {
329      chart.Series.Clear();
330
331
332      for (int i = 0; i < 10; i++) {
333        var s = new Series();
334        s.ChartType = SeriesChartType.FastLine;
335
336        var noise = new Normal(0, 0.05);
337        var x = new double[] {
338          0,
339          1,
340          2,
341          3,
342          4
343        };
344        var eps = noise.Samples().Take(x.Length).ToArray();
345        var y = new double[]
346        {
347          16 + eps[0],
348          9 + eps[1],
349          4 + eps[2],
350          1 + eps[3],
351          0 + eps[4]
352      };
353        var spline = MathNet.Numerics.Interpolation.CubicSpline.InterpolateNatural(x, y
354          //, SplineBoundaryCondition.SecondDerivative, 0, SplineBoundaryCondition.FirstDerivative, 0
355          );
356        for (double xi = 0.0; xi < 10; xi += 0.1) {
357          s.Points.AddXY(xi, spline.Interpolate(xi));
358          Console.WriteLine("{0:N2} {1:N2}", xi, spline.Interpolate(xi));
359        }
360        chart.Series.Add(s);
361      }
362    }
363
364    // public void InterpolateF()
365    // {
366    //   // extrapolating function
367    //   // 6 unknowns
368    //   // f(x,y)     = xy + x² + y² + x + y + 1
369    //   // df(x,y)/dx =              +2x + y + 1
370    //   // df(x,y)/dy =              + x +2y + 1
371    //
372    //   double[,] m =
373    //   {
374    //     // known f
375    //     { 0, 1, 0, 1, 0, 1 }, //  f(1,0)
376    //     { 0, 0, 1, 0, 1, 1 }, //  f(0,1)
377    //     { 1, 1, 1, 1, 1, 1 }, //  f(1,1)
378    //     // known df
379    //     { 0, 0, 0, 2, 0, 1 }, // df/dx(1,0)
380    //     { 0, 0, 0, 0, 1, 1 }, // df/dx(0,1)
381    //     { 0, 0, 0, 1, 0, 1 }, // df/dy(1,0)
382    //     { 0, 0, 0, 0, 2, 1 }, // df/dy(0,1)
383    //   }
384    //
385    // }
386
387    private double F(double x, double y) {
388      return x * x + y * y + x * y + x + y + 1;
389    }
390
391    private double dFdx(double x, double y) {
392      return 2 * x + 1;
393    }
394
395    private double dFdy(double x, double y) {
396      return 2 * y + x + 1;
397    }
398  }
399}
Note: See TracBrowser for help on using the repository browser.