// Poly.scala    polyfit and polyval

object Poly
{
  var S = Simeq

  // polyval  simple Horners evaluation y = p[0] + p[1]*x + ...p[pwr]x^pwr
  def polyval(pwr:Int, p:Array[Double], x:Double):Double=
  {
    var y = p(pwr)*x
    var j = 0
    for(i <- 0 until pwr-1)
    {
      j = pwr-1-i
      y = (p(j)+y)*x
    }
    return y + p(0)
  }

  // polyfit    polyfit(pwr, n, x, y, p)   p computed
  // i=0, i<n     y(i) = p(0) + p(1)*x(i) + p(2)*x(i)^2 + ... x(i)^pwr
  // The problem is stated as follows :
  // Given measured data for values of y based on values of x. e.g.
  //
  //  i    y_actual         x
  // ---   --------       -----
  //  0     32.5           1.0
  //  1      7.2           2.0
  //  2      6.9           3.2
  //  3      8.4           2.7
  //  4      9.6           1.4
  //
  // Find p() such that  
  // y_approximate =  p(0) + p(1)*x(i) + p(2)*x(i)^2 ... +p(pwr)*x(i)^pwr
  // and such that the sum of (y_actual - y_approximate) squared is minimized.
  // To find the coefficients  p  solve the linear system of equations
  //
  // |SUM(1*1)   SUM(1*X)   SUM(1*x^2)   SUM(1*x^3)  | |p(0)| |SUM(1*y)  |
  // |SUM(x*1)   SUM(x*x)   SUM(x*x^2)   SUM(x*x^3)  |x|p(1)|=|SUM(x*y)  |
  // |SUM(x^2*1) SUM(x^2*x) SUM(x^2*x^2) SUM(x^2*x^3)| |p(2)| |SUM(x^2*y)|
  // |SUM(x^3*1) SUM(x^3*x) SUM(x^3*x^2) SUM(x^3*x^3)| |p(3)| |SUM(x^3*y)|
  //   solve  |A|*|p|=|Y| using simeq
  

  def polyfit(pwr:Int, n:Int, x:Array[Double], y:Array[Double]):Array[Double]=
  {
    var pn = pwr+1;
    var xpwr = Array.ofDim[Double](pn)
    var A =    Array.ofDim[Double](pn,pn)
    var Y =    Array.ofDim[Double](pn)
    var p =    Array.ofDim[Double](pn)

    for(i <- 0 until pn) // zero A and Y
    {
      for(j <- 0 until pn)
      {
        A(i)(j) = 0.0
      }
      Y(i) = 0.0
    }
    for(k <- 0 until n)
    {
      xpwr(0) = 1.0
      for(i <- 1 until pn) xpwr(i) = xpwr(i-1)*x(k)
      for(i <- 0 until pn)
      {
        for(j <- 0 until pn)
        {
          A(i)(j) = A(i)(j) + xpwr(i)*xpwr(j)
        }
        Y(i) = Y(i) + y(k)*xpwr(i)
      }
    }

    p = S.solve(A, Y)
    return p
  } // end polyfit

} // end Poly.scala