// test_nuderiv.swift
// file nuderiv.swift included below, should be import

print("test nuderiv.test")
var npoint = 7
var x = [0.1, 0.2, 0.4, 0.6, 0.8, 0.9, 1.1]
var y = [Double](repeating:0,count:npoint)   // function
var yp = [Double](repeating:0,count:npoint)  // first deriv
var ypp = [Double](repeating:0,count:npoint) // second deriv
var cx = [Double](repeating:0,count:npoint)
var cxx = [Double](repeating:0,count:npoint)
var deriv = 0.0
var order = 1

for i in 0..<npoint {
  y[i] = x[i]*x[i]*x[i] + 2.0*x[i]*x[i]
  yp[i] = 3.0*x[i]*x[i] + 4.0*x[i]
  ypp[i] = 6.0*x[i] + 4.0
}
// check deriv at all points i
order = 1
for i in 0..<npoint {
  print("order=\(order), npoint=\(npoint), point=\(i), x")
  print("x=\(x)")
  cx = nuderiv(order, npoint, i, x)
  print("cx=\(cx)")
  deriv = 0.0
  for j in 0..<npoint {
    deriv = deriv + cx[j]*y[j]
  }
  print("deriv=\(deriv), should be=\(yp[i]) \n")
}
// check deriv at all points i
order = 2
for i in 0..<npoint {
  print("order=\(order), npoint=\(npoint), point=\(i), x")
  print("x=\(x)")
  cxx = nuderiv(order, npoint, i, x)
  print("cxx=\(cxx)")
  deriv = 0.0
  for j in 0..<npoint {
    deriv = deriv + cxx[j]*y[j]
  }
  print("deriv=\(deriv), should be=\(ypp[i]) \n")
}
print("test_nuderiv.swift finished")

// nuderiv.swift  non uniformly spaced derivative coefficients 
//    given  x0, x1, x2, x3  non uniformly spaced
//    find the coefficients of y0, y1, y2, y3
//    in order to compute the derivatives:
//      y'(x0) = c00*y0 + c01*y1 + c02*y2 + c03*y3
//      y'(x1) = c10*y0 + c11*y1 + c12*y2 + c13*y3
//      y'(x2) = c20*y0 + c21*y1 + c22*y2 + c23*y3
//      y'(x3) = c30*y0 + c31*y1 + c32*y2 + c33*y3
//
//  Method:  fit data to  y(x)    = a + b*x + c*x^2 + d*x^3
//                        y'(x)   =     b   + 2*c*x + 3*d*x^2
//                        y''(x)  =           2*c   + 6*d*x
//                        y'''(x) =                   6*d
//
//           with y0, y1, y2 and y3 variables:
//
//                              y0  y1  y2  y3 
//  form:  | 1  x0  x0^2  x0^3   1   0   0   0 | = a
//         | 1  x1  x1^2  x1^3   0   1   0   0 | = b
//         | 1  x2  x2^2  x2^3   0   0   1   0 | = c
//         | 1  x3  x3^2  x3^3   0   0   0   1 | = d
//
// reduce  | 1   0     0     0  a0  a1  a2  a3 | = a
//         | 0   1     0     0  b0  b1  b2  b3 | = b
//         | 0   0     1     0  c0  c1  c2  c3 | = c
//         | 0   0     0     1  d0  d1  d2  d3 | = d
//
//                              this is just the inverse
//
//   y'(x0) =  b0*y0        + b1*y1        + b2*y2        + b3*y3        +
//             2*c0*y0*x0   + 2*c1*y1*x0   + 2*c2*y2*x0   + 2*c3*y3*x0   +
//             3*d0*y0*x0^2 + 3*d1*y1*x0^2 + 3*d2*y2*x0^2 + 3*d3*y3*x0^2
//
//   or  c00 = b0 + 2*c0*x0 + 3*d0*x0^2
//       c01 = b1 + 2*c1*x0 + 3*d1*x0^2
//       c02 = b2 + 2*c2*x0 + 3*d2*x0^2
//       c03 = b3 + 2*c3*x0 + 3*d3*x0^2
//
//   y'(x0) = c00*y0 + c01*y1 + c02*y2 +c03*y3
//
//   y'(x1) =  b0*y0        + b1*y1        + b2*y2        + b3*y3        +
//             2*c0*y0*x1   + 2*c1*y1*x1   + 2*c2*y2*x1   + 2*c3*y3*x1   +
//             3*d0*y0*x1^2 + 3*d1*y1*x1^2 + 3*d2*y2*x1^2 + 3*d3*y3*x1^2
//
//   or  c10 = b0 + 2*c0*x1 + 3*d0*x1^2
//       c11 = b1 + 2*c1*x1 + 3*d1*x1^2
//       c12 = b2 + 2*c2*x1 + 3*d2*x1^2
//       c13 = b3 + 2*c3*x1 + 3*d3*x1^2
//
//       cij = bj + 2*cj*xi + 3*dj*xi^2
//
//
//   y'(x1) = c10*y0 + c11*y1 + c12*y2 +c13*y3
//
//   y''(x0) = 2*c0*y0    + 2*c1*y1    + 2*c2*y2    + 2*c3*y3    +
//             6*d0*y0*x0 + 6*d1*y1*x0 + 6*d2*y2*x0 + 6*d3*y3*x0
//
//   or  c00 = 2*c0 + 6*d0*x0
//       c01 = 2*c1 + 6*d1*x0
//       c02 = 2*c2 + 6*d2*x0
//       c03 = 2*c3 + 6*d3*x0
//
//   y'''(x0) = 6*d0*y0 + 6*d1*y1 + 6*d2*y2 + 6*d3*y3
//
//   or  c00 = 6*d0    independent of x
//       c01 = 6*d1
//       c02 = 6*d2
//       c03 = 6*d3
//
// uses  matinv  from  matrix.swift, inserted below

func nuderiv(_ order:Int,_ npoint:Int,_ point:Int,_ x:[Double])->[Double]{
  let n = npoint
  var B = [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
  var fct = [Double](repeating:0,count:n)
  var c = [Double](repeating:0,count:n)
  var pwr = 0.0

  for i in 0..<n {   // build matrix that will be inverted 
    pwr = 1.0
    for j in 0..<n {
      B[i][j] = pwr
      pwr = pwr*x[i]
    }
  }
  let A = matinv(B)  // invert the matrix 
  // form derivative for order and point 
  for i in 0..<n { // compute factors for order 
    fct[i] = 1.0
  }
  for j in 0..<order {
    for i in 0..<n {
      fct[i] = fct[i]*(Double)(i-j)
    }
  }
  for i in 0..<n {
    c[i] = 0.0
    pwr = 1.0
    for j in order..<n {
      c[i] = c[i] + fct[j]*A[j][i]*pwr
      pwr = pwr * x[point]
    }
  }
  return c
} // end nuderiv 


// matinv.swift     matrix invert  should be imported
func matinv(_ a: [[Double]]) -> [[Double]] {
  let n = a.count
  var inv = [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
  for i in 0..<n {
    for j in 0..<n {
      inv[i][j] = a[i][j]
    }
  }

  var row = [Int](repeating:0,count:n)
  var col = [Int](repeating:0,count:n)
  var temp = [Double](repeating:0,count:n)
  var hold = 0
  var I_pivot = 0
  var J_pivot = 0
  var pivot = 0.0
  var abs_pivot = 0.0

  // set up row and column interchange vectors
  for k in 0..<n {
    row[k] = k
    col[k] = k
  }
  // begin main reduction loop
  for k in 0..<n {
    // find largest element for pivot
    pivot = inv[row[k]][col[k]]
    I_pivot = k
    J_pivot = k
    for i in k..<n {
      for j in k..<n {
        abs_pivot = abs(pivot)
        if abs(inv[row[i]][col[j]]) > abs_pivot {
          I_pivot = i
          J_pivot = j
          pivot = inv[row[i]][col[j]]
        }
      }
    }
    if abs(pivot) < 1.0E-10 {
      print("Matrix is singular !")
      return [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
    }
    hold = row[k]
    row[k] = row[I_pivot]
    row[I_pivot] = hold
    hold = col[k]
    col[k] = col[J_pivot]
    col[J_pivot] = hold
    // reduce about pivot
    inv[row[k]][col[k]] = 1.0 / pivot
    for j in 0..<n {
      if j != k {
        inv[row[k]][col[j]] = inv[row[k]][col[j]] * inv[row[k]][col[k]]
      }
    }
    // inner reduction loop
    for i in 0..<n {
      if k != i {
        for j in 0..<n {
          if k != j {
            inv[row[i]][col[j]] = inv[row[i]][col[j]] - inv[row[i]][col[k]] * inv[row[k]][col[j]]
          }
        }
        inv[row[i]][col [k]] = -inv[row[i]][col[k]] * inv[row[k]][col[k]]
      }
    }
  }
  // end main reduction loop

  // unscramble rows
  for j in 0..<n {
    for i in 0..<n {
      temp[col[i]] = inv[row[i]][j]
    }
    for i in 0..<n {
      inv[i][j] = temp[i]
    }
  }
  // unscramble columns
  for i in 0..<n {
    for j in 0..<n {
      temp[row[j]] = inv[i][col[j]]
    }
    for j in 0..<n {
      inv[i][j] = temp[j]
    }
  }
  return inv
} // end matinv   that should be imported

// end nuderiv.swift

// end test_nuderiv.swift