using System; using System.Linq; using System.Windows.Forms; using System.Windows.Forms.DataVisualization.Charting; using MathNet.Numerics.Distributions; namespace Testing { public partial class TestForm : Form { public TestForm() { InitializeComponent(); //Interpolate(); MinAdaptiveGradientSampling(); } // E P p v T Tmax double[,] data = new double[,] { {0.8368, 1.3357, 0.8023, 12.3839, 53.7792, 201.9808}, {0.5197, 1.4065, 0.5648, 17.5391, 53.0914, 162.0023}, {0.8406, 1.4284, 0.5660, 17.5387, 52.9766, 205.0202}, {1.1622, 1.4325, 0.5670, 17.5386, 53.2497, 234.6500}, {0.5208, 0.8740, 0.3042, 17.5388, 52.7885, 148.3409}, {0.8411, 0.8848, 0.3054, 17.5390, 52.7723, 172.0350}, {1.1628, 0.8866, 0.3060, 17.5396, 53.1371, 186.4282}, {0.8472, 1.9431, 0.8010, 17.5390, 52.8836, 221.6862}, {1.1679, 1.9221, 0.8020, 17.5391, 53.0837, 258.8506}, {1.5070, 1.9104, 0.8031, 17.5384, 53.3303, 286.1709}, {0.5085, 1.3330, 0.8003, 12.3837, 53.2506, 156.4186}, {0.8374, 1.3556, 0.8015, 12.3837, 53.1234, 200.7889}, {1.4964, 1.3682, 0.8036, 12.3838, 53.6585, 252.7940}, {0.5039, 0.7592, 0.5658, 9.5036, 53.6265, 149.1639}, {0.8291, 0.7754, 0.5671, 9.5032, 53.5168, 182.6419}, {1.4895, 0.7920, 0.5689, 9.5036, 53.9511, 219.4112}, {0.8373, 1.3749, 0.8015, 12.3835, 53.5531, 201.8539}, {0.5038, 0.7309, 0.8014, 6.7081, 53.7301, 152.8836}, {0.8278, 0.7412, 0.8031, 6.7077, 53.8957, 192.8725}, {1.0551, 0.7459, 0.8041, 6.7081, 54.1965, 213.8598}, {0.8374, 1.3772, 0.8011, 12.3836, 53.4341, 202.3147}, {0.5108, 1.2544, 1.4784, 6.7082, 53.7971, 162.4619}, {0.8344, 1.2516, 1.4843, 6.7078, 53.9234, 215.3835}, {1.0614, 1.2461, 1.4868, 6.7079, 54.0830, 244.6310}, {0.8370, 1.3392, 0.8009, 12.3842, 53.4277, 200.7755}, {0.5167, 1.6944, 2.0975, 6.7078, 53.8257, 166.3573}, {0.8405, 1.6901, 2.1060, 6.7078, 53.8974, 224.3230}, {1.0674, 1.6778, 2.1095, 6.7082, 54.1107, 257.8381}, {0.8371, 1.3172, 0.8010, 12.3839, 53.4347, 199.6043}, {0.8329, 1.2290, 0.2155, 24.8521, 52.6891, 175.7116}, {1.1958, 0.9784, 0.2163, 24.8525, 52.8742, 197.7365}, {1.5166, 0.9721, 0.2168, 24.8520, 53.0470, 213.4699}, {0.8383, 1.4018, 0.8009, 12.3835, 53.5269, 201.9488}, {0.8453, 1.5948, 0.3972, 24.8515, 52.8218, 205.8914}, {1.1740, 1.5958, 0.3981, 24.8518, 52.7543, 236.6865}, {1.4963, 1.6219, 0.3988, 24.8524, 53.1004, 260.2091}, {0.8383, 1.4094, 0.8008, 12.3838, 53.6704, 202.8565}, {0.8538, 2.1021, 0.5648, 24.8523, 52.6869, 221.3082}, {1.1752, 2.0947, 0.5656, 24.8519, 52.8542, 260.5521}, {1.4978, 2.0919, 0.5661, 24.8529, 53.2062, 291.1189}, {0.8383, 1.3991, 0.8004, 12.3842, 53.5893, 201.8672}, {0.8423, 1.8295, 1.4793, 9.5039, 53.7088, 219.3284}, {1.1670, 1.8139, 1.4836, 9.5038, 53.8615, 262.5016}, {1.5025, 1.7922, 1.4866, 9.5035, 54.1248, 297.1750}, {0.8377, 1.3349, 0.8005, 12.3838, 53.6432, 198.6240}, {0.8371, 1.2622, 0.8137, 12.3844, 101.9483, 240.1088}, {0.5196, 1.3274, 0.5751, 17.5397, 102.0049, 202.1158}, {0.8409, 1.3472, 0.5760, 17.5391, 100.8094, 239.7598}, {1.1625, 1.3534, 0.5763, 17.5390, 100.9805, 265.2706}, {0.5225, 0.8233, 0.3109, 17.5394, 100.0307, 186.0851}, {0.8433, 0.8225, 0.3116, 17.5395, 99.9612, 207.5332}, {1.1640, 0.8288, 0.3120, 17.5393, 100.3567, 222.7550}, {0.8473, 1.8966, 0.8123, 17.5390, 100.1643, 254.6133}, {1.1681, 1.8691, 0.8132, 17.5392, 100.3375, 289.9712}, {1.5075, 1.8382, 0.8136, 17.5392, 100.7590, 315.6024}, {0.5081, 1.2802, 0.8109, 12.3841, 100.6839, 196.4331}, {0.8371, 1.2978, 0.8123, 12.3839, 100.6797, 238.7059}, {1.4964, 1.2975, 0.8145, 12.3835, 101.5454, 286.6940}, {0.5038, 0.7123, 0.5749, 9.5036, 101.3295, 185.5991}, {0.8289, 0.7297, 0.5757, 9.5036, 101.3065, 214.1675}, {1.4891, 0.7365, 0.5764, 9.5035, 102.1764, 244.6417}, {0.8374, 1.3175, 0.8121, 12.3840, 101.3808, 240.0907}, {0.5037, 0.7071, 0.8116, 6.7080, 102.1794, 191.0549}, {0.8278, 0.7135, 0.8142, 6.7079, 102.5190, 225.8846}, {1.0549, 0.7150, 0.8156, 6.7079, 102.9568, 243.3557}, {0.8371, 1.3160, 0.8119, 12.3841, 101.1813, 240.1163}, {0.5110, 1.2456, 1.5031, 6.7075, 102.2212, 200.6899}, {0.8345, 1.2408, 1.5095, 6.7075, 102.5058, 249.4029}, {1.0615, 1.2257, 1.5114, 6.7077, 103.0208, 276.3802}, {0.8372, 1.2756, 0.8115, 12.3839, 101.2421, 239.1648}, {0.5168, 1.6720, 2.1301, 6.7078, 102.2741, 204.8973}, {0.8403, 1.6587, 2.1390, 6.7079, 102.5263, 258.4755}, {1.0668, 1.6367, 2.1399, 6.7078, 102.9995, 291.2107}, {0.8368, 1.2379, 0.8111, 12.3839, 101.3221, 238.3200}, {0.8191, 1.3475, 0.2200, 24.8519, 99.4185, 213.5561}, {1.1952, 0.9615, 0.2205, 24.8520, 99.7065, 237.2100}, {1.5198, 0.9392, 0.2206, 24.8526, 100.2240, 254.1504}, {0.8373, 1.3301, 0.8111, 12.3842, 101.4987, 239.6880}, {0.8401, 1.5435, 0.4047, 24.8529, 99.6446, 242.0860}, {1.1737, 1.5140, 0.4051, 24.8520, 99.8763, 270.9338}, {1.4966, 1.5101, 0.4053, 24.8527, 100.3614, 295.7137}, {0.8370, 1.3303, 0.8109, 12.3839, 101.5883, 239.4128}, {0.8515, 2.0246, 0.5744, 24.8527, 99.7521, 256.5944}, {1.1740, 2.0306, 0.5749, 24.8528, 99.9245, 293.3329}, {1.4970, 2.0055, 0.5749, 24.8522, 100.4382, 320.3913}, {0.8373, 1.3181, 0.8106, 12.3846, 101.6315, 238.0053}, {0.8417, 1.7571, 1.5057, 9.5039, 102.0292, 251.9787}, {1.1665, 1.7413, 1.5092, 9.5036, 102.2874, 292.1148}, {1.5021, 1.7017, 1.5109, 9.5040, 102.7524, 325.0641}, {0.8365, 1.2463, 0.8105, 12.3839, 101.9477, 236.4320}, {0.8370, 1.2600, 0.8101, 12.3843, 104.2213, 240.013}, }; //private void MinBleic() { // // double[] w = new double[14]; // double[] optW; // // double[,] constraints = new double[,] // { // // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 y // {0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0 } // }; // int[] constraintType = new int[] // { // -1 // }; // // alglib.minbleicstate state; // alglib.minbleicreport report; // // alglib.minbleiccreate(w, out state); // alglib.minbleicsetcond(state, 0, 1E-10, 0.0, 0); // alglib.minbleicsetlc(state, constraints, constraintType); // alglib.minbleicsetgradientcheck(state, 0.001); // alglib.minbleicsetxrep(state, true); // alglib.minbleicoptimize(state, Grad, Rep, null); // alglib.minbleicresults(state, out optW, out report); // //} private AutoDiff.Term f; private void MinAdaptiveGradientSampling() { double[] w = new double[14]; for (int i = 0; i < w.Length; i++) w[i] = 1.0; double[] optW; double epsx = 0.00001; double radius = 100; double rho = 10.0; int maxits = 0; int numEqConstraints = 0; int numIneqConstraints = 6; alglib.minnsstate state; alglib.minnsreport report; alglib.minnscreate(w, out state); alglib.minnssetnlc(state, numEqConstraints, numIneqConstraints); alglib.minnssetalgoags(state, radius, rho); alglib.minnssetcond(state, 0, 0); alglib.minnssetxrep(state, true); alglib.minnsoptimize(state, Jac, Rep, null); alglib.minnsresults(state, out optW, out report); } private void Rep(double[] w, double func, object o) { Console.WriteLine(Math.Sqrt(func / data.GetLength(0))); } private void Jac(double[] w, double[] fi, double[,] jac, object o) { Array.Clear(fi, 0, fi.Length); Array.Clear(jac, 0, jac.Length); // objective function double f = 0.0; double[] fgrad = new double[w.Length]; Grad(w, ref f, fgrad, o); fi[0] = f; Copy(fgrad, jac, 0); // d Tmax / d Tmin > 0 for all valid inputs fi[1] = -dTmaxOverdE(3, 4, w, fgrad); Copy(fgrad, jac, 1); fi[2] = -dTmaxOverdE(3, 0.1, w, fgrad); Copy(fgrad, jac, 2); fi[3] = -dTmaxOverdE(0.1, 4, w, fgrad); Copy(fgrad, jac, 3); fi[4] = -dTmaxOverdE(0.1, 0.1, w, fgrad); Copy(fgrad, jac, 4); // d Tmax / d p > 0 for all valid inputs fi[5] = -dTmaxOverdp(3, 3, 0.1, 50, w, fgrad); Copy(fgrad, jac, 5); fi[6] = -dTmaxOverdp(3, 0.5, 0.1, 50, w, fgrad); Copy(fgrad, jac, 6); // d Tmax / d v > 0 for all valid inputs // d Tmax / d E > 0 for all valid inputs } private void Copy(double[] from, double[,] to, int col) { for (int i = 0; i < from.Length; i++) { to[col, i] = from[i]; } } private void Grad(double[] w, ref double func, double[] grad, object o) { // double ssq = 0.0; for (int i = 0; i < grad.Length; i++) grad[i] = 0.0; var fgrad = new double[grad.Length]; for (int r = 0; r < data.GetLength(0); r++) { var pred = Tmax(data[r, 0], data[r, 2], data[r, 3], data[r, 4], w, fgrad); var target = data[r, 5]; var err = pred - target; ssq += err * err; for (int i = 0; i < grad.Length; i++) { grad[i] += -2.0 * target * fgrad[i]; grad[i] += 2.0 * pred * fgrad[i]; } } func = ssq; } private double Tmax(double E, double p, double v, double Tmin, double[] w, double[] grad) { grad[0] = Tmin; grad[1] = (v); grad[2] = (p); grad[3] = (E); grad[4] = (Math.Log((p)) * v); grad[5] = (Math.Log((p)) * E); grad[6] = (v * p); grad[7] = (Tmin * p); grad[8] = (E * p); grad[9] = (v * E); grad[10] = (Tmin * E); grad[11] = (E * E); grad[12] = (E * 1 / p); grad[13] = (1 / p); return w[0] * (Tmin) + w[1] * (v) + w[2] * (p) + w[3] * (E) + w[4] * (Math.Log((p)) * v) + w[5] * (Math.Log((p)) * E) + w[6] * (v * p) + w[7] * (Tmin * p) + w[8] * (E * p) + w[9] * (v * E) + w[10] * (Tmin * E) + w[11] * (E * E) + w[12] * (E * 1 / p) + w[13] * (1 / p); } private double dTmaxOverdp(double E, double p, double v, double Tmin, double[] w, double[] grad) { grad[0] = 0; grad[1] = 0; grad[2] = 1; grad[3] = 0; grad[4] = (v * 1/p); grad[5] = (E * 1/p); grad[6] = v ; grad[7] = Tmin; grad[8] = E; grad[9] = 0; grad[10] = 0; grad[11] = 0; grad[12] = -(E * 1 / (p*p)); grad[13] = -(1 / (p*p)); return w[2] + w[4] * (v * 1/p) + w[5] * (E * 1/p) + w[6] * (v ) + w[7] * (Tmin ) + w[8] * (E) + w[12] * -(E * 1 / (p*p)) + w[13] * -(1 / (p*p)); } private double dTmaxOverdE(double E, double p, double[] w, double[] grad) { grad[0] = 1; grad[1] = 0; grad[2] = 0; grad[3] = 0; grad[4] = 0; grad[5] = 0; grad[6] = 0; grad[7] = p; grad[8] = 0; grad[9] = 0; grad[10] = E; grad[11] = 0; grad[12] = 0; grad[13] = 0; return w[0] + w[7] * (p) + w[10] * (E); } private void Interpolate() { chart.Series.Clear(); for (int i = 0; i < 10; i++) { var s = new Series(); s.ChartType = SeriesChartType.FastLine; var noise = new Normal(0, 0.05); var x = new double[] { 0, 1, 2, 3, 4 }; var eps = noise.Samples().Take(x.Length).ToArray(); var y = new double[] { 16 + eps[0], 9 + eps[1], 4 + eps[2], 1 + eps[3], 0 + eps[4] }; var spline = MathNet.Numerics.Interpolation.CubicSpline.InterpolateNatural(x, y //, SplineBoundaryCondition.SecondDerivative, 0, SplineBoundaryCondition.FirstDerivative, 0 ); for (double xi = 0.0; xi < 10; xi += 0.1) { s.Points.AddXY(xi, spline.Interpolate(xi)); Console.WriteLine("{0:N2} {1:N2}", xi, spline.Interpolate(xi)); } chart.Series.Add(s); } } // public void InterpolateF() // { // // extrapolating function // // 6 unknowns // // f(x,y) = xy + x² + y² + x + y + 1 // // df(x,y)/dx = +2x + y + 1 // // df(x,y)/dy = + x +2y + 1 // // double[,] m = // { // // known f // { 0, 1, 0, 1, 0, 1 }, // f(1,0) // { 0, 0, 1, 0, 1, 1 }, // f(0,1) // { 1, 1, 1, 1, 1, 1 }, // f(1,1) // // known df // { 0, 0, 0, 2, 0, 1 }, // df/dx(1,0) // { 0, 0, 0, 0, 1, 1 }, // df/dx(0,1) // { 0, 0, 0, 1, 0, 1 }, // df/dy(1,0) // { 0, 0, 0, 0, 2, 1 }, // df/dy(0,1) // } // // } private double F(double x, double y) { return x * x + y * y + x * y + x + y + 1; } private double dFdx(double x, double y) { return 2 * x + 1; } private double dFdy(double x, double y) { return 2 * y + x + 1; } } }