## Copyright 2021 Ian Jauslin
## 
## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License.
## You may obtain a copy of the License at
## 
##     http://www.apache.org/licenses/LICENSE-2.0
## 
## Unless required by applicable law or agreed to in writing, software
## distributed under the License is distributed on an "AS IS" BASIS,
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
## See the License for the specific language governing permissions and
## limitations under the License.

# linear interpolation: given vectors x,y, compute a linear interpolation for y(x0)
# assume x is ordered
@everywhere function linear_interpolation(x0,x,y)
  # if x0 is beyond all x's, then return the corresponding boundary value.
  if x0>x[length(x)]
    return y[length(y)]
  elseif x0<x[1]
    return y[1]
  end

  # find bracketing interval
  i=bracket(x0,x,1,length(x))

  # interpolate
  return y[i]+(y[i+1]-y[i])*(x0-x[i])/(x[i+1]-x[i])
end
@everywhere function bracket(x0,x,a,b)
  i=floor(Int64,(a+b)/2)
  if x0<x[i]
    return bracket(x0,x,a,i)
  elseif x0>x[i+1]
    return bracket(x0,x,i+1,b)
  else
    return i
  end
end


# polynomial interpolation of a family of points
@everywhere function poly_interpolation(x,y)
  # init for recursion
  rec=Array{Polynomial{Float64}}(undef,length(x))
  for i in 1:length(x)
    rec[i]=Polynomial([1.])
  end

  # compute \prod (x-x_i)
  poly_interpolation_rec(rec,x,1,length(x))

  # sum terms together
  out=0.
  for i in 1:length(y)
    out+=rec[i]/rec[i](x[i])*y[i]
  end

  return out
end
# recursive helper function
@everywhere function poly_interpolation_rec(out,x,a,b)
  if a==b
    return
  end
  # mid point
  mid=floor(Int64,(a+b)/2)
  # multiply left and right of mid
  right=Polynomial([1.])
  for i in mid+1:b
    right*=Polynomial([-x[i],1.])
  end
  left=Polynomial([1.])
  for i in a:mid
    left*=Polynomial([-x[i],1.])
  end
  
  # multiply into left and right
  for i in a:mid
    out[i]*=right
  end
  for i in mid+1:b
    out[i]*=left
  end

  # recurse
  poly_interpolation_rec(out,x,a,mid)
  poly_interpolation_rec(out,x,mid+1,b)

  return
end
## the following does the same, but has complexity N^2, the function above has N*log(N)
#@everywhere function poly_interpolation(x,y)
#  out=Polynomial([0.])
#  for i in 1:length(x)
#    prod=Polynomial([1.])
#    for j in 1:length(x)
#      if j!=i
#	prod*=Polynomial([-x[j]/(x[i]-x[j]),1/(x[i]-x[j])])
#      end
#    end
#    out+=prod*y[i]
#  end
#
#  return out
#end