// test_simeq_newton5.java  solve nonlinear system of equations
//                          method: newton iteration using Jacobian
//                          use list for higher order terms
//
//  may have 1 to about 10 variables, X1, X2, X3, ... 
//  a maximal term can have up to power of 3 divided by power of 2 for each var
//  n is number of variables, nlin is number of non linear terms
//  
//  in this implementation (if we had symbolic differentiation, no limit)
//  a term could be  x1 , x1*x3 , x2/x5 , x1*x2*x3/(x4*x5) ,  x2^3/x4^2 
//
//  for testing, generate A using pseudo random numbers
//               choose x1=1.1 x2=1.2 x3=1.4 x4  compute products
//               and reciprocals
//               compute terms of Y using Y = A X
//
//  Solve by initial guess at values of x1, x2, x3 computing products
//    X_next = X_initial - J_initial^-1 * (A *  X_initial - Y)
//    in general X_next = X_prev - (J_prev^-1 * (A * X_prev - Y))*b
//                                 where 0 < b < 1, often 0.5, for stability
//
//  solved when  abs sum each row A * X_next -Y < epsilon
//
//  It may stall, stop if abs(X_next-X_prev)<epsilon
//  It may diverge, stop, indicate no solution
//                        (or try a different initial guess)
//  It may oscillate, stop, indicate no solution
//                        (or try a different initial guess)
//
//  Having unknown values for x1, x2, x3
//  Solve x1 + x2 + x3 + x1*x2 + x3*x2*x2 + x2/x3 + x1*x1*x1/(x3*x3) = y
//     there are 7 terms, thus you need 7 equations, constants A and y i=1..7
//     A[i,1]*x1 + A[i,2]*x2 + A[i,3]*x3 + A[i,4]*x1*x2 +
//     A[i,5]*x3*x2*x2 + A[i,6]*x2/x3 + A[i,7]*x1*x1*x1/(x3*x3)=y[i]
//
//  The matrix equation becomes: A * X_soln = Y
//
//   x1    x2    x3   x1*x2 x3*x2^2 x2/x3 x1^3/x3^2  
//
// | A1,1  A1,2  A1,3  A1,4   A1,5  A1,6  A1,7 | | x1         | |y1 |
// | A2,1  A2,2  A2,3  A2,4   A2,5  A2,6  A2,7 | | x2         | |y2 |
// | A3,1  A3,2  A3,3  A3,4   A3,5  A3,6  A3,7 | | x3         | |y3 |
// | A4,1  A4,2  A4,3  A4,4   A4,5  A4,6  A4,7 |*| x1*x2      | |y4 |
// | A5,1  A5,2  A5,3  A5,4   A5,5  A5,6  A5,7 | | x3*x2^2    |=|y5 |
// | A6,1  A6,2  A6,3  A6,4   A6,5  A6,6  A6,7 | | x2/x3      | |y6 |
// | A7,1  A7,2  A7,3  A7,4   A7,5  A7,6  A7,7 | | x1^3/x3^2  | |y7 |
//   A constants                          unknown terms, soln    y constants
//
//  The Jacobian, J is the numerically computed derivatives based on A, X_prev
//  The derivative coefficients may be computed once based on A and the
//  unknown variables. The numeric value plugs in the X_prev values.
//
//  J1,1 = A1,1 + A1,4*x2 + 3*A1,7*x1^2/x3^2               deriv Row 1 wrt x1
//  J1,2 = A1,2 + A1,4*x1 + 2*A1,5*x3*x2 + A1,6/x3         deriv Row 1 wrt x2
//  J1,3 = A1,3 + A1,5*x2^2 + -A1,6*x2/x3^2 + -2*A1,7*x1^2/x3^3  
//                                                         deriv Row 1 wrt x3
//  J2,1 = A2,1 + A2,4*x2 + 3*A2,7*x1^2/x3^2               deriv Row 2 wrt x1
//  J2,2 = A2,2 + A2,4*x1 + 2*A2,5*x3*x2 + A2,6/x3         deriv Row 2 wrt x2
//  J2,3 = A1,3 + A2,5*x2^2 + -A2,6*x2/x3^2 + -2*A2,7*x1^2/x3^3  
//                                                         deriv Row 2 wrt x3
//  J3,1 = A3,1 + A3,4*x2 + 3*A3,7*x1^2/x3^2               deriv Row 3 wrt x1
//  J3,2 = A3,2 + A3,4*x1 + 2*A3,5*x3*x2 + A3,6/x3         deriv Row 3 wrt x2
//  J3,3 = A3,3 + A3,5*x2^2 + -A3,6*x2/x3^2 + -2*A3,7*x1^2/x3^3  
//                                                         deriv Row 3 wrt x3
//
//  The non linear equation is given by arrays, 7 terms, subscript 0..6
//   index = 0   1   2     3          4      5         6
//  var1     1   2   3     1          3      2         1
//  var2    -1  -1  -1     1          2     -1         1  -1 is unused variable
//  var3    -1  -1  -1    -1          2     -1         1
//  vari1   -1  -1  -1    -1         -1      3         3
//  vari2   -1  -1  -1    -1         -1     -1         3
//          x1  x2  x3  x1*x2  x3*x2*x2  x2/x3  x1*x1*x1/x3*x3   term
//
//  Code or a data structure must be available to know the equation
//  of the entries in the Y vector. Symbolic computation of
//  derivatives is assumed to be available, when needed.
//

// needs  simeq_newton5.java compiled

public class test_simeq_newton5
{
  simeq_newton5 SN = new simeq_newton5(); // default
  String Xname[] = {"X0", "X1", "X2", "X3", "X4", "X5"}; // all linear terms
  double A[][];    // = new double[n][n+nlin];
  double Y[];      // = new double[n];
  double X[];      // = new double[n+nlin];
  double X_soln[]; // = new double[n+nlin];

  public test_simeq_newton5()
    // build (double A[][], double Y[], int var1[], int var2[], double X[])
  {
    System.out.println("test_simeq_newton5.java running");
  { // local block
    int n = 2;   // total number of terms
    int nlin = 2; // linear unknowns, required, even if zero coef
    System.out.println("test case 1, n="+n+", nlin="+nlin);
    System.out.println("equation X[0] + X[1] + X[0]*X[0] + X[1]*X[1] = Y");
    //              --n-- --nlin--
    int var1[]  = { 0,  1,  0,  1}; // zero based subscript
    int var2[]  = {-1, -1,  0,  1}; // e.g. last is x3*x4
    int var3[]  = {-1, -1, -1, -1}; // possible cube or 3 term
    int vari1[] = {-1, -1, -1, -1}; // zero based subscript
    int vari2[] = {-1, -1, -1, -1}; // e.g. last is x3*x4

    A = new double[n][n+nlin];
    Y = new double[n];
    X = new double[n+nlin];
    X_soln = new double[n+nlin];

    // build test case 1
    // build_A(n, nlin, var1, var2, var3, vari1, vari2);
    for(int i=0; i<n; i++)
      for(int j=0; j<n+nlin; j++)
        A[i][j] = 1.0;

    A[0][3] = 0.0; // x + y + x^2 = 4
    A[1][2] = 0.0; // x + y + y^2 = 7
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n+nlin; j++)
      {
	System.out.println("A["+i+"]["+j+"]="+A[i][j]);
      }
    }
    System.out.println(" ");

    // choose X's
    // arbitrary test solution for random A matrix
    X_soln[0] = 1.0;
    X_soln[1] = 2.0;

    build_Y(n, nlin, var1, var2, var3, vari1, vari2);

    // initial guess
    X[0] = 0.5;
    X[1] = 0.5;

    print_eqn(n, nlin, var1, var2, var3, vari1, vari2);
   
    SN.simeq(n, nlin, A, Y, var1, var2, var3, vari1, vari2, X);

    print_results(n, nlin);

    System.out.println("test 1 finished \n");
  } // end local block 1

  { // local block 2
    int n = 3;   // total number of terms
    int nlin = 4; // linear unknowns, required, even if zero coef
    System.out.println("test case 2, n="+n+", nlin="+nlin);
    //             x1  x2  x3 x1*x2 x3*x2^2
    //                                  x2/x3
    //                                     x1^3/x3^2 ... 
    int var1[]  = { 0,  1,  2,  0,  2,  1, 0}; // zero based subscript
    int var2[]  = {-1, -1, -1,  1,  1, -1, 0}; // e.g. last is x3*x4
    int var3[]  = {-1, -1, -1, -1,  1, -1, 0}; // possible cube or 3 term
    int vari1[] = {-1, -1, -1, -1, -1,  2, 2}; // zero based subscript
    int vari2[] = {-1, -1, -1, -1, -1, -1, 2}; // e.g. last is x3*x4

    A = new double[n][n+nlin];
    Y = new double[n];
    X = new double[n+nlin];
    X_soln = new double[n+nlin];

    // build test case 2
    build_A(n, nlin, var1, var2, var3, vari1, vari2);

    // choose x1=1.1 x2=1.2 x3=1.4  compute x1*x2, etc
    // arbitrary test solution for random A matrix
    X_soln[0] = 1.1;
    X_soln[1] = 1.2;
    X_soln[2] = 1.4;

    build_Y(n, nlin, var1, var2, var3, vari1, vari2);

    // initial guess
    for(int i=0; i<n; i++) X[i] = 1.0;

    print_eqn(n, nlin, var1, var2, var3, vari1, vari2);
   
    SN.simeq(n, nlin, A, Y, var1, var2, var3, vari1, vari2, X);

    print_results(n, nlin);

    System.out.println("test 2 finished \n");
  } // end local block 2

  { // local block 3
    int n = 4;   // total number of terms
    int nlin = 3; // linear unknowns, required, even if zero coef
    System.out.println("test case 3, n="+n+", nlin="+nlin);
    int var1[]  = { 0,  1,  2,  3,  0,  0,  3}; // zero based subscript
    int var2[]  = {-1, -1, -1, -1,  1,  0,  3}; // e.g. last is x3*x4
    int var3[]  = {-1, -1, -1, -1, -1,  2,  3}; // possible cube or 3 term
    int vari1[] = {-1, -1, -1, -1, -1, -1, -1}; // zero based subscript
    int vari2[] = {-1, -1, -1, -1, -1, -1, -1}; // e.g. last is x3*x4

    A = new double[n][n+nlin];
    Y = new double[n];
    X = new double[n+nlin];
    X_soln = new double[n+nlin];

    // build test case 3
    build_A(n, nlin, var1, var2, var3, vari1, vari2);

    // choose x1=1.1 x2=1.2 x3=1.4  compute x1*x2, etc
    // arbitrary test solution for random A matrix
    X_soln[0] = 1.1;
    X_soln[1] = 1.2;
    X_soln[2] = 1.4;
    X_soln[3] = 1.5;

    build_Y(n, nlin, var1, var2, var3, vari1, vari2);

    // initial guess
    for(int i=0; i<n; i++) X[i] = 1.0;

    print_eqn(n, nlin, var1, var2, var3, vari1, vari2);
   
    SN.simeq(n, nlin, A, Y, var1, var2, var3, vari1, vari2, X);

    print_results(n, nlin);

    System.out.println("test 3 finished \n");
  } // end local block 3

  { // local block 4
    int n = 4;   // total number of terms
    int nlin = 6; // linear unknowns, required, even if zero coef
    System.out.println("test case 4, n="+n+", nlin="+nlin);
    int var1[]  = { 0,  1,  2,  3,  0,  1,  2,  3,  0,  1};
    int var2[]  = {-1, -1, -1, -1,  0,  1,  2,  3,  1,  2};
    int var3[]  = {-1, -1, -1, -1,  0,  1,  2,  3,  2,  3};
    int vari1[] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1};
    int vari2[] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1};

    A = new double[n][n+nlin];
    Y = new double[n];
    X = new double[n+nlin];
    X_soln = new double[n+nlin];

    // build test case 4
    build_A(n, nlin, var1, var2, var3, vari1, vari2);

    // choose x1=1.1 x2=1.2 x3=1.4  compute x1*x2, etc
    // arbitrary test solution for random A matrix
    X_soln[0] = 1.1;
    X_soln[1] = 1.2;
    X_soln[2] = 1.4;
    X_soln[3] = 1.5;

    build_Y(n, nlin, var1, var2, var3, vari1, vari2);

    // initial guess
    X[0] = 1.0;
    X[1] = 1.0;
    X[2] = 1.0;
    X[3] = 1.0;

    print_eqn(n, nlin, var1, var2, var3, vari1, vari2);
   
    SN.simeq(n, nlin, A, Y, var1, var2, var3, vari1, vari2, X);

    print_results(n, nlin);

    System.out.println("test 4 finished \n");
  } // end local block 4

  { // local block 5
    int n = 5;   // total number of terms
    int nlin = 15; // linear unknowns, required, even if zero coef
    System.out.println("test case 5, n="+n+", nlin="+nlin);
    int var1[]  = { 0,  1,  2,  3,  4,  1,  1,  1,  2,  2,  3,
                    1,  1,  1,  2,  2,  2,  3,  3,  3};
    int var2[]  = {-1, -1, -1, -1, -1,  1,  2,  3,  2,  3,  3,
                    1,  1,  1,  2,  2,  2,  3,  3,  3};
    int var3[]  = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                    1,  2,  3,  1,  2,  3,  1,  2,  3};
    int vari1[] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                   -1, -1, -1, -1, -1, -1, -1, -1, -1};
    int vari2[] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                   -1, -1, -1, -1, -1, -1, -1, -1, -1};

    A = new double[n][n+nlin];
    Y = new double[n];
    X = new double[n+nlin];
    X_soln = new double[n+nlin];

    // build test case 5
    build_A(n, nlin, var1, var2, var3, vari1, vari2);

    // choose x1=1.1 x2=1.2 x3=1.4  compute x1*x2, etc
    // arbitrary test solution for random A matrix
    X_soln[0] = 1.1;
    X_soln[1] = 1.2;
    X_soln[2] = 1.4;
    X_soln[3] = 1.5;
    X_soln[4] = 1.7;

    build_Y(n, nlin, var1, var2, var3, vari1, vari2);

    // initial guess
    X[0] = 1.0;
    X[1] = 1.0;
    X[2] = 1.0;
    X[3] = 1.0;
    X[4] = 1.0;

    print_eqn(n, nlin, var1, var2, var3, vari1, vari2);
   
    SN.simeq(n, nlin, A, Y, var1, var2, var3, vari1, vari2, X);

    System.out.println("solution may not be unique");
    print_results(n, nlin);

    // initial guess
    X[0] = 1.1;
    X[1] = 1.2;
    X[2] = 1.4;
    X[3] = 1.5;
    X[4] = 1.7;
   
    SN.simeq(n, nlin, A, Y, var1, var2, var3, vari1, vari2, X);

    System.out.println("give solution as initial guess, check");
    print_results(n, nlin);

    System.out.println("test 5 finished \n");
  } // end local block 5

  { // local block 6
    simeq_newton5 SN24 = new simeq_newton5(1.0e-7, 1.0, 24, false);
    int n = 5;   // total number of terms
    int nlin = 1; // linear unknowns, required, even if zero coef
    System.out.println("test case 4, n="+n+", nlin="+nlin);
    int var1[]  = { 0,  1,  2,  3,  4,  0};
    int var2[]  = {-1, -1, -1, -1, -1,  1};
    int var3[]  = {-1, -1, -1, -1, -1,  2};
    int vari1[] = {-1, -1, -1, -1, -1,  3};
    int vari2[] = {-1, -1, -1, -1, -1,  4};

    A = new double[n][n+nlin];
    Y = new double[n];
    X = new double[n+nlin];
    X_soln = new double[n+nlin];

    // build test case 6
    for(int i=0; i<n; i++)
      for(int j=0; j<n+nlin; j++)
	  A[i][j] = 1.0/(double)(i+j+1);

    // choose x1=1.1 x2=1.2 x3=1.4  compute x1*x2, etc
    // arbitrary test solution for random A matrix
    X_soln[0] = 1.1;
    X_soln[1] = 1.2;
    X_soln[2] = 1.4;
    X_soln[3] = 1.5;
    X_soln[4] = 1.9;

    build_Y(n, nlin, var1, var2, var3, vari1, vari2);

    // initial guess
    for(int i=0; i<n; i++) X[i] = 2.0;

    print_eqn(n, nlin, var1, var2, var3, vari1, vari2);
   
    SN24.simeq(n, nlin, A, Y, var1, var2, var3, vari1, vari2, X);

    print_results(n, nlin);

    System.out.println("test 6 finished \n");
  } // end local block 6

  } // end test_simeq_newton5 constructor


  void build_A(int n, int nlin, int var1[], int var2[],
               int var3[], int vari1[], int vari2[])
  { // using random numbers, may get different results on another run
    for(int i=0; i<n; i++) // all subscripts zero based, one less than comments
    {
      for(int j=0; j<n+nlin; j++)
      {
	  A[i][j] = Math.random();
      }
    }
    System.out.println("first A generated ");
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n+nlin; j++)
      {
	System.out.println("A["+i+"]["+j+"]="+A[i][j]);
      }
    }
    System.out.println(" ");
  } // end build_A

  void build_Y(int n, int nlin, int var1[], int var2[],
               int var3[], int vari1[], int vari2[])
  {
    for(int i=n; i<n+nlin; i++)
    {
      X_soln[i] = 1.0; /* build nonlinear term values */
      if(var1[i]>=0) X_soln[i]  =  X_soln[var1[i]];
      if(var2[i]>=0) X_soln[i]  *= X_soln[var2[i]];
      if(var3[i]>=0) X_soln[i]  *= X_soln[var3[i]];
      if(vari1[i]>=0) X_soln[i] /= X_soln[vari1[i]];
      if(vari2[i]>=0) X_soln[i] /= X_soln[vari2[i]];
    }

    // set up Y
    for(int i=0; i<n; i++)
    {
      Y[i] = 0.0;
      for(int j=0; j<n+nlin; j++)
      {
	Y[i] += A[i][j]*X_soln[j];
      }
    }
    // debug print
    for(int i=0; i<n; i++)
    {
      System.out.println("Y["+i+"]="+Y[i]);
    }
    System.out.println(" ");
  } // end build_Y

  void print_eqn(int n, int nlin, int var1[], int var2[],
                 int var3[], int vari1[], int vari2[])
  {
    System.out.println("solve system of equations A X = Y for X");
    System.out.println("the equations are for i=0,"+(n-1));
    for(int i=0; i<n; i++) System.out.println("A[i]["+i+"]*"+Xname[i]+"+");
    for(int i=n; i<n+nlin; i++)
    {
      if(var1[i]>=n || var2[i]>=n || var3[i]>=n ||
         vari1[i]>=n || vari2[i]>=n)
	  System.out.println("var error "+var1[i]+" "+var2[i]+" "+var3[i]+
                             " "+vari1[i]+" "+vari2[i]);
      System.out.print("A[i]["+i+"]");
      if(var1[i]>=0) System.out.print("*"+Xname[var1[i]]);
      if(var2[i]>=0) System.out.print("*"+Xname[var2[i]]);
      if(var3[i]>=0) System.out.print("*"+Xname[var3[i]]);
      if(vari1[i]>=0 || vari2[i]>=0) System.out.print("/(");
      if(vari1[i]>=0) System.out.print(Xname[vari1[i]]);
      if(vari2[i]>=0) System.out.print("*"+Xname[vari2[i]]);
      if(vari1[i]>=0 || vari2[i]>=0) System.out.print(")");
      if(i==n+nlin-1) System.out.println(" ");
      else       System.out.println("+");
    }
    System.out.println("  = Y[i]");
    System.out.println(" ");

    System.out.println("desired solution, may not be unique");
    for(int i=0; i<n; i++)
    {
      System.out.println("X_soln["+i+"]="+X_soln[i]);
    }
    System.out.println(" ");

    System.out.println("initial guess");
    for(int i=0; i<n; i++)
    {
      System.out.println("X["+i+"]="+X[i]);
    }
    System.out.println(" ");

  } // end print_eqn

  void print_results(int n, int nlin)
  {
    System.out.println("Returned solution vs expected solution");

    for(int i=0; i<n; i++) System.out.println("X["+i+"]="+X[i]+
	"       err="+(X[i]-X_soln[i]));

    double err = 0.0;
    for(int i=0; i<n; i++)
    {
      double aerr = 0.0;
      for(int j=0; j<n+nlin; j++)
      {
	aerr += A[i][j]*X[j];
      }
      err += Math.abs(aerr-Y[i]);
    }
    System.out.println("Returned solution in given equation, sum of errors="+
                       err);
    System.out.println(" ");
  } // end print_results

  public static void main (String[] args)
  {
      new test_simeq_newton5();
  }
} // end class test_simeq_newton5