# test_simeq_newton5.rb  solve nonlinear system of equations
#                     method: newton iteration using jacobian
#                     use list for higher order terms
#
#  a maximal term can have up to power of 3 divided by power of 2
#  in this implementation 
#  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)
#
#  The matrix equation could be:
#
#   x1    x2    x3   x1*x2 x3*x2^2 x2/x3 x1^3/x3^2 ... 
#   var1  var1  var1  var2   var3     varil  vari2
#   a1    a2    a3    a4     a5       a6     a7    coefficients for each equation
#   n=3 variables x1, x2, x3
#   nlin=4 nonlinear terms  power 2, power 3, p0123/power 1  p0123/power 2
#   thus 3+4= 7 equations in 7 unknowns  (index 1 will be index 0 in code)
#
#
# | 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 |
# | ...                                       ...| |            | |   |
#       zero based variable numbers  x1 is 0, x2 is 1, x3 is 2, -1 means none
# var1    0   1   2   0   2   1   0
# var2   -1  -1  -1   1   1  -1   0
# var3   -1  -1  -1  -1   1  -1   0
# vari1  -1  -1  -1  -1  -1   2   2
# vari2-  1  -1  -1  -1  -1  -1   2
#
#  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
#  ...
#  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_newton5.rb compiled
#
# test_simeq_newton5.rb

require_relative 'Simeq_newton5'

def build_a(n, nlin, var, var2, var3, vari1, vari2)
  # using random numbers, may get different results on another run
  puts "build_a runnung"
  for i in 0...n  # all subscripts zero based, one less than comments
    for j in 0...(n+nlin)
      a[i][j] = Math.random()
    end # j
  end # i
  puts "first a generated "
  for i in 0...n
    for j in 0...(n+nlin)
      puts "a[#{i}][#{j}]=#{a[i][j]}"
    end # j
  end # i
  puts " "
end # build_a

def build_y(n, nlin, a, x_soln, y, var1, var2, var3, vari1, vari2)
  puts "build_y running"
  for i in n...nlin
    x_soln[i] = 1.0 # build nonlinear term values 
    if var1[i]>=0
      x_soln[i]  =  x_soln[var1[i]]
    end # if
    if var2[i]>=0
      x_soln[i] = x_soln[i] * x_soln[var2[i]]
    end # if
    if var3[i]>=0
      x_soln[i] = x_soln[i] * x_soln[var3[i]]
    end # if
    if vari1[i]>=0
      x_soln[i] = x_soln[i] / x_soln[vari1[i]]
    end # if
    if vari2[i]>=0
      x_soln[i] = x_soln[i] / x_soln[vari2[i]]
    end # if
  end # i
  # debug print
  for i in 0...(n+nlin)
    puts "x_soln[#{i}]=#{x_soln[i]}"
  end # i
  puts " "

  # set up y
  for i in 0...n
    y[i] = 0.0
    for j in 0...(n+nlin)
      y[i] = y[i] + a[i][j]*x_soln[j]
    end # j
  end # i
  # debug print
  for i in 0...n
    puts "y[#{i}]=#{y[i]}"
  end # i
  puts " "
end # build_y

def print_eqn(n, nlin, a, x, y, var1, var2, var3, vari1, vari2, x_soln)
  xname = ["x0", "x1", "x2", "x3", "x4", "x5"] # all linear terms
  puts "print_eqn running"
  puts "solve system of equations a * x = y for x"
  puts "the equations are for i=0,#{(n-1)}"
  for i in 0...n
    puts "a[i][#{i}]*#{xname[i]} + "
  end # i
  for i in n...(n+nlin)
    if var1[i]>=n || var2[i]>=n || var3[i]>=n || vari1[i]>=n || vari2[i]>=n
      puts "var error #{var1[i]} #{var2[i]} #{var3[i]} #{vari1[i]} #{vari2[i]}"
      return
    else
      # puts "a[#{i}][#{i}]= #{a[i][i]}"
      if var1[i]>=0
        puts " * #{xname[var1[i]]}"
      end # if
      if var2[i]>=0
        puts " * #{xname[var2[i]]}"
      end # if
      if var3[i]>=0
        puts " * #{xname[var3[i]]}"
      end # if
      if vari1[i]>=0 || vari2[i]>=0
        puts " /("
      end # if
      if vari1[i]>=0
        puts "#{xname[vari1[i]]}"
      end # if
      if vari2[i]>=0
        puts " * #{xname[vari2[i]]}"
      end # if
      if vari1[i]>=0 || vari2[i]>=0
        puts ")"
      end # if
      if i==n+nlin-1
        puts " "
      else
        puts " + "
      end # if
    end # if
    puts "  = y[i]"
    puts " "
  end # i
  
  puts "desired solution, may not be unique"
  for i in 0...n
      puts "x_soln[#{i}]=#{x_soln[i]}"
  end # i
  puts " "

  puts "initial guess"
  for i in 0...n
    puts "x[#{i}]=#{x[i]}"
  end # i
  puts " "
end # print_eqn

def print_results(n, nlin, a, x, y, x_soln)
  puts "print_results running, Returned solution vs expected solution"
  for i in 0...n
    puts "x[#{i}]=#{x[i]}  err=#{(x[i]-x_soln[i])}"
  end # i
  err = 0.0
  for i in 0...n
    aerr = 0.0
    for j in 0...(n+nlin)
      aerr = aerr + a[i][j]*x[j]
    end # j
    err = err + (aerr-y[i]).abs
  end # i
  puts "Returned solution in given equation, sum of errors=#{err}"
  puts " "
end # print_results

puts "test_simeq_newton5.rb running"
#  Simeq_newton5.new.simeq_newton5(n, nlin, a, y, vec1, vec2, vec3,veci1,veci2)
xname = ["x0", "x1", "x2", "x3", "x4", "x5"] # all linear terms
# build test case 1
n = 2   # total number of terms
nlin = 2 # linear unknowns, required, even if zero coef
# 2 equations with 4 terms 
a = Array.new(n){Array.new(n+nlin)}
y = Array.new(n)
x = Array.new(n+nlin)
x_soln = Array.new(n+nlin)
# build_a(n, nlin, a, y, var1, var2, var3, vari1, vari2, x) not used
#
# a1*x1 + a2*x2 + a3*x1^2 + a4*x2^2 = y1  equation
#
#        x1  x2  x1^2  x2^2
var1 =  [ 0,  1,    0,    1]  # zero based subscript
var2 =  [-1, -1,    0,    1]  # e.g. last is x3*x4
var3 =  [-1, -1,   -1,   -1]  # possible cube or 3 term
vari1 = [-1, -1,   -1,   -1] # zero based subscript
vari2 = [-1, -1,   -1,   -1] # e.g. last is x3*x4
puts "test case 1, n=#{n}, nlin=#{nlin}"
puts "var1   #{var1[0]}   #{var1[0]}   #{var1[2]}   #{var1[3]}"
puts "var2   #{var2[0]}   #{var2[0]}   #{var2[2]}   #{var2[3]}"
puts "var3   #{var3[0]}   #{var3[0]}   #{var3[2]}   #{var3[3]}"
puts "vari1  #{vari1[0]}   #{vari1[0]}   #{vari1[2]}   #{vari1[3]}"
puts "vari2  #{vari2[0]}   #{vari2[0]}   #{vari2[2]}   #{vari2[3]}"
puts " "
for i in 0...n
  for j in 0...(n+nlin)
    a[i][j] = 1.0
  end # j
end # j
a[0][3] = 0.0 # x + y + x^2 = 4
a[1][2] = 0.0 # x + y + y^2 = 7
for i in 0...n
  for j in 0...(n+nlin)
    puts "a[#{i}][#{j}]=#{a[i][j]}"
  end # j
end # i
puts " "

# choose x's
# arbitrary test solution if random a matrix
x_soln[0] = 1.0
x_soln[1] = 2.0
x_soln[2] = 1.0
x_soln[3] = 4.0

build_y(n, nlin, a, x_soln, y, var1, var2, var3, vari1, vari2)

# initial guess
x[0] = 0.5
x[1] = 0.5

print_eqn(n, nlin, a, x, y, var1, var2, var3, vari1, vari2, x_soln)
   
Simeq_newton5.new.simeq(n, nlin, a, y, var1, var2, var3, vari1, vari2, x)

print_results(n, nlin, a, x, y, x_soln)
puts "test case 1 finished"