using Statistics
"""
   linear1D(data,xgrid, x) 

linear interpolation of data(x)

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- x: x
"""

@inline function linear1D(data, xgrid, x)

    xarray = xgrid.grid

    xi0,xi1 = 0,0
    if(x<=xarray[firstindex(xgrid)])
        xi0=1
        xi1=2
    elseif(x>=xarray[lastindex(xgrid)])
        xi0=lastindex(xgrid)-1
        xi1=xi0+1
    else
        xi0=floor(xgrid,x)
        xi1=xi0+1
    end

    dx0, dx1 = x - xarray[xi0], xarray[xi1] - x

    d0, d1 = data[xi0], data[xi1]

    g = d0 * dx1 + d1 * dx0

    gx = g / (dx0 + dx1) 
    return gx
end


"""
   linear2D(data, xgrid, ygrid, x, y) 

linear interpolation of data(x, y)

#Arguments:
- xgrid: one-dimensional grid of x
- ygrid: one-dimensional grid of y
- data: two-dimensional array of data
- x: x
- y: y
"""
@inline function linear2D(data, xgrid, ygrid, x, y)

    xarray, yarray = xgrid.grid, ygrid.grid

    # if (
    #     x <= xarray[firstindex(xgrid)] ||
    #     x >= xarray[lastindex(xgrid)] ||
    #     y <= yarray[firstindex(ygrid)] ||
    #     y >= yarray[lastindex(ygrid)]
    # )
    #     return 0.0
    # end
    # xi0, yi0 = floor(xgrid, x), floor(ygrid, y)
    # xi1, yi1 = xi0 + 1, yi0 + 1

    xi0,xi1,yi0,yi1 = 0,0,0,0
    if(x<=xarray[firstindex(xgrid)])
        xi0=1
        xi1=2
    elseif(x>=xarray[lastindex(xgrid)])
        xi0=lastindex(xgrid)-1
        xi1=xi0+1
    else
        xi0=floor(xgrid,x)
        xi1=xi0+1
    end

    if(y<=yarray[firstindex(ygrid)])
        yi0=1
        yi1=2
    elseif(y>=yarray[lastindex(ygrid)])
        yi0=lastindex(ygrid)-1
        yi1=yi0+1
    else
        yi0=floor(ygrid,y)
        yi1=yi0+1
    end

    dx0, dx1 = x - xarray[xi0], xarray[xi1] - x
    dy0, dy1 = y - yarray[yi0], yarray[yi1] - y

    d00, d01 = data[xi0, yi0], data[xi0, yi1]
    d10, d11 = data[xi1, yi0], data[xi1, yi1]

    g0 = data[xi0, yi0] * dx1 + data[xi1, yi0] * dx0
    g1 = data[xi0, yi1] * dx1 + data[xi1, yi1] * dx0

    gx = (g0 * dy1 + g1 * dy0) / (dx0 + dx1) / (dy0 + dy1)
    return gx
end


"""
testInterpolation1D(func, grid1, grid2)

Test interpolation on grid1 with function func, compare results with func(grid2).
Return max deviation and std of deviation.

#Arguments:
- func: function needs to be interpolated
- grid1: grid to be tested
- grid2: grid to be interpolated
"""
function testInterpolation1D(func, grid1, grid2, isRelativeError=false, find=false)
    data1 = [func(grid1[i]) for i in 1:(grid1.size)]

    data2 = [func(grid2[i]) for i in 1:(grid2.size)]
    data2_int = [linear1D(data1, grid1, grid2[i]) for i in 1:(grid2.size)]

    d = data2-data2_int
    if isRelativeError
        d = abs.(d ./ maximum(  abs.(data2) ))
    else
        d = abs.( d )
    end
    # d_max = maximum( d )
    d_max, d_max_index = findmax(d)
    d_std = std(d)
    if find
        return d_max, d_std, d_max_index
    else
        return d_max, d_std
    end
end

# function testInterpolation1D_rel(func, grid1, grid2)
#     data1 = [func(grid1[i]) for i in 1:(grid1.size)]

#     data2 = [func(grid2[i]) for i in 1:(grid2.size)]
#     data2_int = [linear1D(data1, grid1, grid2[i]) for i in 1:(grid2.size)]

    #     d_rel = abs.((data2-data2_int) ./ data2)
#     d_max = maximum(d_rel)
#     d_std = std(d_rel)
#     return d_max, d_std
# end

"""
optimizeUniLog(generator, para, MN, func)

Return a pair of parameters (M, N) which give optimized unilog grid
generated by generator that minimize absolute error when interpolate func.

#Arguments:
- generator: function generator(para, M, N) that generates unilog grid. 
- para: parameters needed by generator
- MN: the result is constrained by (M+1)*N<MN
- func: the function to be interpolated
"""
function optimizeUniLog(generator, para, MN, func)
    Ns = floor(Int, sqrt(MN))
    fine_grid = generator(para, 2Ns, 2Ns)
    f(m, n) = testInterpolation1D(func, generator(para, m, n), fine_grid, true)[2]

    M, N = Ns-1, Ns
    d_curr = f(M, N)

    Na, Nb = floor(Int, sqrt(MN)), floor(Int, sqrt(MN))
    while Na > 1
        if floor(Int, MN/(Nb+1))==Na
            # maximize Nb for current Na
            Nb = Nb + 1
            continue
        end

        # test (Na-1, Nb) and (Nb-1, Na)
        d1 = f(Na-1, Nb)
        d2 = f(Nb-1, Na)

        if d1<=d2 && d1<d_curr
            M, N = Na-1, Nb
            d_curr = d1
        elseif d2<d1 && d2<d_curr
            M, N = Nb-1, Na
            d_curr = d2
        end

        Nb = Nb + 1
        Na = floor(Int, MN/Nb)
    end

    return M, N, d_curr
end
