/***********************************************************************
 * Xpress Optimizer Examples
 * =========================
 *
 *  Polygon_vecmapdelta.cs
 *`````````````````````
 *
 *  (c) 2017 Fair Isaac Corporation
 *
 *
 * Polygon example: maximise the area of an N sided polygon
 *
 * *** Demonstrating using a vector map (R^k->R) userfunction calculating its own derivatives ***
 *
 * Variables:
 *
 *  rho   : 0..N-1   ! Distance of vertex from the base point
 *  theta : 0..N-1   ! Angle from x-axis
 *
 * Objective:
 *   (sum (i in 1..N-2) (rho(i)*rho(i-1)*sin(theta(i)-theta(i-1)))) * 0.5
 *
 * Constraints:
 *   Vertices in increasing degree order:
 *      theta(i) >= theta(i-1) +.01  : i = 1..N-2
 *   Boundary conditions:
 *      theta(N-1) <= Pi
 *      0.1 <= rho(i) <= 1   : i = 0..N-2
 *   Third side of all triangles <= 1
 *      rho(i)^2 + rho(j)^2 - rho(i)*rho(j)*2*cos(theta(j)-theta(i)) <= 1    : i in 0..N-3, j in i..N-2
***********************************************************************/

using System;
using System.IO;
using Optimizer;
using XSLPdnet;

namespace XSLPExamples
{
  class XSLPExample
  {

    // userfunction of type "vecmapdelta", implementing (x0,x1) -> sin(x0-x1) and (x0,x1) -> cos(x0-x1) alongside its partials
    public double mySinDifference(double[] x, double[] deltas, double[] partialDerivatives)
    {
      if (deltas != null) // Delta may be used as a suggestion for a finite difference step size
                          // however it also indicates if a partial is requested, saving on effort in case only an evaluation is needed
      {
        if (deltas[0] != 0.0)
        {
          partialDerivatives[0] = Math.Cos(x[0] - x[1]);
        }
        if (deltas[1] != 0.0)
        {
          partialDerivatives[1] = -Math.Cos(x[0] - x[1]);
        }
      }
      return (Math.Sin(x[0] - x[1]));
    }
    public double myCosDifference(double[] x, double[] deltas, double[] partialDerivatives)
    {
      if (deltas != null) // Delta may be used as a suggestion for a finite difference step size
                          // however it also indicates if a partial is requested, saving on effort in case only an evaluation is needed
      {
        if (deltas[0] != 0.0)
        {
          partialDerivatives[0] = -Math.Sin(x[0] - x[1]);
        }
        if (deltas[1] != 0.0)
        {
          partialDerivatives[1] = Math.Sin(x[0] - x[1]);
        }
      }
      return (Math.Cos(x[0] - x[1]));
    }

    public void RunExp1()
    {
      /* Example 1 loads a simple problem, and setups the SLP coefficient using a string. */
      XPRS.Init("");
      try
      {
        XSLP sProb = new XSLP();
        try
        {

          // Tell Optimizer to call OptimizerMsg whenever a message is output
          sProb.XPRS.AddMessageCallback(this.OptimizerMsg);

          System.Console.WriteLine("####################");
          System.Console.WriteLine("# Polygon example  #");
          System.Console.WriteLine("####################");
          System.Console.WriteLine(sProb.getBanner());

          // Number of side of the Polygon
          int nSide = 5;

          /* Create rows */
          int nRow = nSide - 2 + (nSide - 1) * (nSide - 2) / 2 + 1;
          double[] RHS = new double[nRow];
          for (int i = 0; i < nRow; i++) RHS[i] = 0;

          char[] RowType = new char[nRow + 1];
          string[] RowNames = new string[nRow];
          nRow = 0;

          // Objective constraint
          RowType[nRow] = 'E'; /* OBJEQ */
          RowNames[nRow] = "OBJEQ";
          nRow++;

          // Vertices in increasing degree order
          for (int i = 1; i < nSide - 1; i++)
          {
            RowType[nRow] = 'G'; /* T2T1 .. T4T3 */
            RowNames[nRow] = "T" + (i + 1) + "T" + i;
            nRow++;
            RHS[i] = 0.001;
          }

          // Third side of all triangles <= 1
          for (int i = 1; i < nSide - 1; i++)
          {
            for (int j = i + 1; j < nSide; j++)
            {
              RowType[nRow] = 'L';
              RHS[nRow] = 1.0;
              RowNames[nRow] = "V" + i + "V" + j;
              nRow++;
            }
          }
          RowType[nRow] = '\0';

          /* Columns and linear coefficients*/
          int nCol = (nSide - 1) * 2 + 2;
          int nElement = 0;
          double[] OBJ = new double[nCol];
          double[] Lower = new double[nCol];
          double[] Upper = new double[nCol];
          for (int i = 0; i < nCol; i++)
          {
            OBJ[i] = 0;           /* objective function */
            Lower[i] = 0;         /* lower bound normally zero */
            Upper[i] = XPRS.PLUSINFINITY; /* upper bound infinity */
          }

          // Arrays to load linear coefficients
          int[] ColStart = new int[nRow + 1];
          double[] Element = new double[10 * nCol];
          int[] RowIndex = new int[10 * nCol];
          string[] ColNames = new string[nCol];

          /* OBJX */
          nCol = 0;
          ColStart[nCol] = nElement;
          OBJ[nCol] = 1.0;
          ColNames[nCol] = "OBJX";
          Lower[nCol++] = XPRS.MINUSINFINITY; /* free column */
          Element[nElement] = -1.0;
          RowIndex[nElement++] = 0;

          /* THETA1 - THETA 4 */
          int iRow = 0;
          for (int i = 1; i < nSide; i++)
          {
            ColNames[nCol] = "THETA" + i;
            ColStart[nCol++] = nElement;
            if (i < nSide - 1)
            {
              Element[nElement] = -1;
              RowIndex[nElement++] = iRow + 1;
            }
            if (i > 1)
            {
              Element[nElement] = 1;
              RowIndex[nElement++] = iRow;
            }
            iRow++;
          }

          Upper[nCol - 1] = 3.1415926;

          /* Rho's */
          for (int i = 1; i < nSide; i++)
          {
            Lower[nCol] = 0.01;   /* lower bound */
            Upper[nCol] = 1;
            ColNames[nCol] = "RHO" + i;
            ColStart[nCol++] = nElement;
          }
          ColStart[nCol] = nElement;

          sProb.XPRS.LoadLP("Polygon", nCol, nRow, RowType, RHS, null, OBJ, ColStart, null, RowIndex, Element, Lower, Upper);

          sProb.XPRS.AddNames(1, RowNames, 0, nRow - 1);
          sProb.XPRS.AddNames(2, ColNames, 0, nCol - 1);

          /* Add the user function*/
          sProb.AddUserFunction("mySinDifference", 2 /* number if inputs arguments */, 0, mySinDifference);
          sProb.AddUserFunction("myCosDifference", 2 /* number if inputs arguments */, 0, myCosDifference);

          /* Build up nonlinear coefficients */

          /* Area */
          string CoefBuffer = "0.5 * ( ";
          for (int i = 1; i < nSide - 1; i++)
          {
            if (i > 1)
            {
              CoefBuffer += " + ";
            }
            // Expression using userfunction "mySinDifference" in place of the builtin sine function
            CoefBuffer += "RHO" + (i + 1) + " * RHO" + i + " * mySinDifference ( THETA" + (i + 1) + " , THETA" + i + " )";
          }
          CoefBuffer += " )";
          sProb.ChgStringFormula(0, CoefBuffer);

          /* Distances */
          for (int i = 1; i < nSide - 1; i++)
          {
            for (int j = i + 1; j < nSide; j++)
            {
              // Expression using userfunction "myCosDifference" in place of the builtin cosine function
              CoefBuffer = "RHO" + j + " ^ 2 + RHO" + i + " ^ 2 - 2 * RHO" + j + " * RHO" + i + " * myCosDifference ( THETA" + j + " , THETA" + i + " )";
              sProb.ChgStringFormula(iRow, CoefBuffer);
              iRow++;
            }
          }

          sProb.Maxim("");
        }
        finally
        {
          /* Cleanup */
          sProb.destroyProb();
        }
      }
      finally
      {
        XPRS.Free();
      }
    }

    public void OptimizerMsg(XPRSprob prob, object data, string message, int len, int msglvl)
    {
      if (message != null)
      {
        Console.WriteLine(message);
      }
    }


    public static void Main(string[] args)
    {
      XSLPExample e = new XSLPExample();
      e.RunExp1();
      Console.ReadLine();
    }
  }
}