// test_simeq_newton3.java  solve nonlinear system of equations
//                     method: newton iteration using Jacobian
//                     use list for higher order terms
//
//  Given problem  A X = Y  where X may have terms x1, x2, x3, x4,
//                 and higher order such as: x1*x2, x1*x3, x1*x4, x2*x3, ...
//                 A sparse matrix may be used, coefficients given
//                 Y is vector of reals given
//                 independent unknowns are x1, x2, x3, x4
//
//  for testing, generate A using pseudo random numbers
//               choose x1=1.1 x2=1.2 x3=1.4 x4 =1.5, compute products
//               compute terms of Y using Y = A X
//
//  Solve by initial guess at values of x1, x2, x3, x4 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)
//
//
//  The user matrix equation for this case is: y = A * x
//  Given A and y find x
//
//   x1     x2     x3     x4     x1*x2 ... x3*x4  (may be subset of nl terms) 
//
// | A1,1   A1,2   A1,3   A1,4   A1,5 ...  A1,10  | | x1    | |y1 |
// | A2,1   A2,2   A2,3   A2,4   A2,5 ...  A2,10  | | x2    | |y2 |
// | A3,1   A3,2   A3,3   A3,4   A3,5 ...  A3,10  | | x3    | |y3 |
// | A4,1   A4,2   A4,3   A4,4   A4,5 ...  A4,10  |*| x4    |=|y4 |
//                                                  | x1*x2 |
//                                                  | x1*x3 | 
//                                                  | ...   | 
//                                                  | x3*x4 |
//
//   var1 ... n                  var3 ... nlin
//
//  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,5*x2 + A1,6*x3 + A1,7*x4   deriv Row 1 wrt x1
//  J1,2 = A1,2 + A1,5*x1 + A1,8*x3 + A1,9*x4   deriv Row 1 wrt x2
//  J1,3 = A1,3 + A1,6*x1 + A1,8*x2 + A1,10*x4  deriv Row 1 wrt x3
//  J1,4 = A1,4 + A1,7*x1 + A1,9*x2 + A1,10*x3  deriv Row 1 wrt x4
//  J1,5 = A1,5                                 deriv Row 1 wrt x1, wrt x2
//  J1,6 = A1,6                                 deriv Row 1 wrt x1, wrt x3
//  J1,7 = A1,7                                 deriv Row 1 wrt x1, wrt x4
//  J1,8 = A1,8                                 deriv Row 1 wrt x2, wrt x3
//  J1,9 = A1,9                                 deriv Row 1 wrt x2, wrt x4
//  J1,10 = A1,10                               deriv Row 1 wrt x3, wrt x4
//
//  J2,1 = A2,1 + A2,5*x2 + A2,6*x3 + A2,7*x4   deriv Row 2 wrt x1
//  J2,2 = A2,2 + A2,5*x1 + A2,8*x3 + A2,9*x4   deriv Row 2 wrt x2
//
//  etc.
//
//  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_newton3.java compiled

public class test_simeq_newton3
{
  public test_simeq_newton3()
       // build (double A[][], double Y[], int var1[], int var2[], double X[])
  {
    int n = 4;    // total number of unknown variables  rows of A matrix
    int nlin = 6; // nonlinear unknowns, products var2
    int ntot = n + nlin;    // columns of A matrix
    int var1[] = { 0, 1, 2, 3, 0, 0, 0, 1, 1, 2}; // zero based subscript
    int var2[] = {-1,-1,-1,-1, 1, 2, 3, 2, 3, 3}; // e.g. last is x3*x4
    int var3[] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}; // possible cube or 3 term
    System.out.println("n="+n+"  nlin="+nlin);
    System.out.println("x0*x1, x0*x2, x0*x3, x1*x2, x1*x3, x2*x3");
    // could have int var4[] etc for higher powers
    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];
    String Xname[] = {"X0", "X1", "X2", "X3"}; // all linear terms
    String Xname1[] = {"X1", "X2", "X3", "X4"}; // all linear terms
    double start, now;

    System.out.println("test_simeq_newton3.java running");
    start = System.currentTimeMillis()/1000.0;

    // build first test case
    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("A["+i+"]["+j+"]= "+A[i][j]);
      }
    }
    System.out.println("first A generated random");
    // choose x1=1.1 x2=1.2 x3=1.4 x4 =1.5, compute x1*x2, x3*x4
    // 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;
    System.out.println("solve system of equations A X = Y for X = ");
    for(int i=0; i<n; i++) System.out.print(Xname1[i]+", ");
    System.out.println(" ");
    for(int i=n; i<n+nlin; i++)
    {
      if(var1[i]>=n || var2[i]>=n || var3[i]>=n)
	  System.out.println("var error "+var1[i]+" "+var2[i]+" "+var3[i]);
      System.out.print(Xname1[var1[i]]);
      if(var2[i]>=0) System.out.print("*"+Xname1[var2[i]]);
      if(var3[i]>=0) System.out.println("*"+Xname1[var3[i]]);
      if(i==n-1) System.out.println("\n");
      else       System.out.println(",");
    }

    for(int i=n; i<n+nlin; i++)
    {
      X_soln[i] = X_soln[var1[i]]; // at least one required
      if(var2[i]>=0) X_soln[i] *= X_soln[var2[i]];
      if(var3[i]>=0) X_soln[i] *= X_soln[var3[i]];
    }
    // debug print
    for(int i=0; i<n; i++)
    {
      System.out.println("X_soln["+i+"]="+X_soln[i]);
    }
    System.out.println(" ");

    // 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(" ");
    // initial guess
    for(int i=0; i<n; i++) X[i] = 1.0;
   
    new simeq_newton3(n, nlin, A, Y, var1, var2, var3, X,
                      10, 1.0, 1.0e-6, 3 );

    System.out.println("Returned solution");

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

    // second  test case
    // int var1[] = { 0, 1, 2, 3, 0, 1 ,2, 3, 0, 2}; //  zero based subscript
    // int var2[] = {-1,-1,-1,-1, 0, 1, 2, 3, 1, 3}; // e.g. last is x3*x4
    // int var3[] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};
    //                            --squared--
    var1[4]=0; var2[4]=0; // var3 initialized to -1
    var1[5]=1; var2[5]=1;
    var1[6]=2; var2[6]=2;
    var1[7]=3; var2[7]=3;
    var1[8]=0; var2[8]=1;
    var1[9]=2; var2[9]=3;

    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("second A generated ");
    // choose x1=1.1 x2=1.2 x3=1.4 x4 =1.5, compute x1*x2, x3*x4
    // 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;
    System.out.println("solve system of equations A X = Y for X = ");
    for(int i=0; i<n; i++) System.out.print(Xname[i]+", ");
    System.out.println(" ");
    for(int i=n; i<n+nlin; i++)
    {
      if(var1[i]>=n || var2[i]>=n || var3[i]>=n)
	  System.out.println("var error "+var1[i]+" "+var2[i]+" "+var3[i]);
      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(i==n-1) System.out.println("\n");
      else       System.out.println(",");
    }

    for(int i=n; i<n+nlin; i++)
    {
      X_soln[i] = X_soln[var1[i]]; // at least one required
      if(var2[i]>=0) X_soln[i] *= X_soln[var2[i]];
      if(var3[i]>=0) X_soln[i] *= X_soln[var3[i]];
    }
    // debug print
    for(int i=0; i<n+nlin; i++)
    {
      System.out.println("X_soln["+i+"]="+X_soln[i]);
    }
    System.out.println(" ");

    // 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(" ");
    // initial guess
    for(int i=0; i<n; i++) X[i] = 1.0;
   
    new simeq_newton3(n, nlin, A, Y, var1, var2, var3, X,
                      10, 1.0, 1.0e-6, 3 );

    System.out.println("Returned solution");

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

    // third  test case
    // int var1[] = { 0, 1, 2, 3, 0, 1 ,2, 3, 0, 1}; //  zero based subscript
    // int var2[] = {-1,-1,-1,-1, 0, 1, 2, 3, 1, 2}; //  e.g. last is
    // int var3[] = {-1,-1,-1,-1, 0, 1, 2, 3, 2, 3}; //  x2*x3*x4
    //                            ---cubed---
    var1[4]=0; var2[4]=0; var3[4]=0;
    var1[5]=1; var2[5]=1; var3[5]=1;
    var1[6]=2; var2[6]=2; var3[6]=2;
    var1[7]=3; var2[7]=3; var3[7]=3;
    var1[8]=0; var2[8]=1; var3[8]=2;
    var1[9]=1; var2[9]=2; var3[9]=3;

    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("third A generated ");
    // choose x1=1.1 x2=1.2 x3=1.4 x4 =1.5, compute x1*x2, x3*x4
    // 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;
    System.out.println("solve system of equations A X = Y for X = ");
    for(int i=0; i<n; i++) System.out.print(Xname[i]+", ");
    System.out.println(" ");
    for(int i=n; i<n+nlin; i++)
    {
      if(var1[i]>=n || var2[i]>=n || var3[i]>=n)
	  System.out.println("var error "+var1[i]+" "+var2[i]+" "+var3[i]);
      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(i==n-1) System.out.println("\n");
      else       System.out.println(",");
    }

    for(int i=nlin; i<n; i++)
    {
      X_soln[i] = X_soln[var1[i]]; // at least one required
      if(var2[i]>=0) X_soln[i] *= X_soln[var2[i]];
      if(var3[i]>=0) X_soln[i] *= X_soln[var3[i]];
    }
    // debug print
    for(int i=0; i<n; i++)
    {
      System.out.println("X_soln["+i+"]="+X_soln[i]);
    }
    System.out.println(" ");

    // set up Y
    for(int i=0; i<n; i++)
    {
      Y[i] = 0.0;
      for(int j=0; j<n; 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(" ");
    // initial guess
    for(int i=0; i<n; i++) X[i] = 0.5;
   
    new simeq_newton3(n, nlin, A, Y, var1, var2, var3, X,
                      10, 1.0, 1.0e-6, 2 );

    System.out.println("Returned solution");

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

    System.out.println("Fourth test case ");
    n = 2;
    nlin = 2;
    double A1[][] = new double[n][n+nlin];
    var1[2]=0; var2[2]=0; var3[2]=-1;
    var1[3]=1; var2[3]=1; var3[5]=-1;

    A1[0][0] = 1.0;
    A1[0][1] = 1.0;
    A1[0][2] = 1.0; // x[0]^2 
    A1[0][3] = 0.0; // x[1]^2 
    A1[1][0] = 1.0;
    A1[1][1] = 1.0;
    A1[1][2] = 0.0;
    A1[1][3] = 1.0;

    System.out.println("fourth A generated ");
    // choose x1=1 x2=1  compute x1*x1, ... 
    X_soln[0] = 1.0;
    X_soln[1] = 2.0;
    System.out.println("solve system of equations A X = Y for X = ");
    for(int i=0; i<n; i++) System.out.print(Xname[i]+", ");
    System.out.println(" ");
    for(int i=n; i<n+nlin; i++)
    {
      if(var1[i]>=n || var2[i]>=n || var3[i]>=n)
	  System.out.println("var error "+var1[i]+" "+var2[i]+" "+var3[i]);
      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(i==n-1) System.out.println("\n");
      else       System.out.println(",");
    }

    for(int i=n; i<n+nlin; i++)
    {
      X_soln[i] = X_soln[var1[i]]; // at least one required 
      if(var2[i]>=0) X_soln[i] *= X_soln[var2[i]];
      if(var3[i]>=0) X_soln[i] *= X_soln[var3[i]];
    }
    for(int i=0; i<n+nlin; i++)
    {
      System.out.println("X_soln["+i+"]="+X_soln[i]);
    }
    System.out.println(" ");

    // set up Y 
    for(int i=0; i<n; i++)
    {
      Y[i] = 0.0;
      for(int j=0; j<n+nlin; j++)
      {
        Y[i] += A1[i][j]*X_soln[j];
      }
    }
    for(int i=0; i<n; i++)
    {
      System.out.println("Y["+i+"]="+Y[i]);
    }
    System.out.println(" ");
    // initial guess 
    X[0] = 0.5;
    X[1] = 0.5;

    new simeq_newton3(n, nlin, A1, Y, var1, var2, var3, X,
                      10, 1.0, 1.0e-6, 1 );

    System.out.println("Returned solution ");

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

    now = System.currentTimeMillis()/1000.0;
    System.out.println("test_simeq_newton3 finished in "+(now-start)+" seconds");
  }

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