source: branches/2789_MathNetNumerics-Exploration/Testing/Form.cs @ 16421

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;
    }
  }
}