Add 1D interpolation.

Author: sriharshakandala

.dev/Project.toml

@@ -0,0 +1,2 @@
+[deps]
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"

.dev/climaformat.jl

@@ -0,0 +1,83 @@
+#!/usr/bin/env julia
+#
+# This is an adapted version of format.jl from JuliaFormatter with the
+# following license:
+#
+#    MIT License
+#    Copyright (c) 2019 Dominique Luna
+#
+#    Permission is hereby granted, free of charge, to any person obtaining a
+#    copy of this software and associated documentation files (the
+#    "Software"), to deal in the Software without restriction, including
+#    without limitation the rights to use, copy, modify, merge, publish,
+#    distribute, sublicense, and/or sell copies of the Software, and to permit
+#    persons to whom the Software is furnished to do so, subject to the
+#    following conditions:
+#
+#    The above copyright notice and this permission notice shall be included
+#    in all copies or substantial portions of the Software.
+#
+#    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+#    OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+#    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
+#    NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+#    DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+#    OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+#    USE OR OTHER DEALINGS IN THE SOFTWARE.
+#
+using Pkg
+Pkg.activate(@__DIR__)
+Pkg.instantiate()
+
+using JuliaFormatter
+
+help = """
+Usage: climaformat.jl [flags] [FILE/PATH]...
+
+Formats the given julia files using the CLIMA formatting options.  If paths
+are given it will format the julia files in the paths. Otherwise, it will
+format all changed julia files.
+
+    -v, --verbose
+        Print the name of the files being formatted with relevant details.
+
+    -h, --help
+        Print this message.
+"""
+
+function parse_opts!(args::Vector{String})
+    i = 1
+    opts = Dict{Symbol,Union{Int,Bool}}()
+    while i ≤ length(args)
+        arg = args[i]
+        if arg[1] != '-'
+            i += 1
+            continue
+        end
+        if arg == "-v" || arg == "--verbose"
+            opt = :verbose
+        elseif arg == "-h" || arg == "--help"
+            opt = :help
+        else
+            error("invalid option $arg")
+        end
+        if opt in (:verbose, :help)
+            opts[opt] = true
+            deleteat!(args, i)
+        end
+    end
+    return opts
+end
+
+opts = parse_opts!(ARGS)
+if haskey(opts, :help)
+    write(stdout, help)
+    exit(0)
+end
+if isempty(ARGS)
+    filenames = readlines(`git ls-files "*.jl"`)
+else
+    filenames = ARGS
+end
+
+format(filenames; opts...)

Project.toml

@@ -0,0 +1,4 @@
+name = "ClimaInterpolations"
+uuid = "dd0f122e-fa3b-47f3-bcf0-93bbc60d885e"
+authors = ["sriharshakandala <[email protected]>"]
+version = "0.1.0"

src/ClimaInterpolations.jl

@@ -0,0 +1,5 @@
+module ClimaInterpolations
+
+include("Interpolation1D.jl")
+
+end

src/Interpolation1D.jl

@@ -0,0 +1,104 @@
+module Interpolation1D
+
+include("stencils1d.jl")
+
+"""
+    interpolate(xpts, fval, xtarg)
+
+Interpolate `fval` defined on `xpts` at `xtarg`.
+Reference:
+1). Berrut, J., and Trefethen, L. N., 2004. Barycentric Lagrange Interpolation, SIAM REVIEW
+Vol. 46, No. 3, pp. 501–517.
+DOI. 10.1137/S0036144502417715
+"""
+Base.@pure function interpolate(
+    xpts::AbstractArray{FT,1},
+    fval::AbstractArray{FT,1},
+    xtarg::FT,
+) where {FT<:AbstractFloat}
+    length(xpts) == 1 && return fval[1] # flat extrapolation
+    l, tot, mi, tol = FT(1), -FT(0), 0, (10 * eps(FT))
+    @inbounds begin
+        for i in eachindex(xpts)
+            dx = xtarg - xpts[i]
+            if abs(dx) < tol
+                mi = i
+                break
+            else
+                l *= dx
+            end
+        end
+        if mi == 0
+            for i in eachindex(xpts)
+                xi = xpts[i]
+                wi = FT(1)
+                for j in eachindex(xpts)
+                    if j ≠ i
+                        wi *= (xi - xpts[j])
+                    end
+                end
+                tot += fval[i] * l / ((xtarg - xi) * wi)
+            end
+        else
+            tot = fval[mi]
+        end
+    end
+    return tot
+end
+
+"""
+    interpolate1d!(
+        ftarget::AbstractArray{FT,N},
+        xsource::AbstractArray{FT,N},
+        xtarget::AbstractArray{FT,N},
+        fsource::AbstractArray{FT,N},
+        order,
+        extrapolate = Flat(),
+    ) where {FT,N}
+
+Interpolate `fsource`, defined on grid `xsource`, onto the `xtarget` grid.
+Here the source grid `xsource` is an N-dimensional array of columns.
+The first dimension is assumed to be the column dimension. 
+Each column can have a different grid.
+"""
+function interpolate1d!(
+    ftarget::AbstractArray{FT,N},
+    xsource::AbstractArray{FT,N},
+    xtarget::AbstractArray{FT,N},
+    fsource::AbstractArray{FT,N},
+    order,
+    extrapolate = Flat(),
+) where {FT,N}
+    @assert N ≥ 1
+    (nsource, coldims_source...) = size(xsource)
+    (ntarget, coldims_target...) = size(xtarget)
+    @assert coldims_source == coldims_target
+    @assert size(ftarget) == size(xtarget)
+    @assert size(fsource) == size(xsource)
+    coldims = coldims_source
+    colcidxs = CartesianIndices(coldims)
+    @inbounds begin
+        for colcidx in colcidxs
+            colidx = Tuple(colcidx)
+            first = 1
+            for i1 = 1:ntarget
+                (st, en) = get_stencil(
+                    order,
+                    view(xsource, :, colidx...),
+                    xtarget[i1, colidx...],
+                    first = first,
+                    extrapolate = extrapolate,
+                )
+                first = st
+                ftarget[i1, colidx...] = interpolate(
+                    view(xsource, st:en, colidx...),
+                    view(fsource, st:en, colidx...),
+                    xtarget[i1, colidx...],
+                )
+            end
+        end
+    end
+    return nothing
+end
+
+end

src/stencils1d.jl

@@ -0,0 +1,32 @@
+
+abstract type Order1D end
+abstract type Extrapolate1D end
+
+struct Linear <: Order1D end
+
+struct Flat <: Extrapolate1D end
+struct LinearExtrapolation <: Extrapolate1D end
+
+function get_stencil(alg::Linear, xsource, xtarg; first = 1, extrapolate = Flat())
+    last = length(xsource)
+    if xtarg < xsource[1] # extrapolation
+        # Linear extrapolation
+        extrapolate == LinearExtrapolation() && return (1, 2)
+        # Flat extrapolation by default
+        return (1, 1)
+    end
+    if xtarg > xsource[last] # extrapolation
+        # Linear extrapolation
+        extrapolate == LinearExtrapolation() && return (last - 1, last)
+        # Flat extrapolation by default
+        return (last, last)
+    end
+    st = first
+    for i = first:last
+        if xtarg ≤ xsource[i]
+            st = max(i - 1, 1)
+            break
+        end
+    end
+    return (st, st < last ? st + 1 : st)
+end

test/Project.toml

@@ -0,0 +1,4 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+ClimaInterpolations = "dd0f122e-fa3b-47f3-bcf0-93bbc60d885e"

test/interpolation1D.jl

@@ -0,0 +1,130 @@
+import ClimaInterpolations.Interpolation1D:
+    interpolate, Linear, get_stencil, interpolate1d!, Flat, LinearExtrapolation
+using Test
+
+
+function get_uniform_column_grids(
+    ::Type{FT},
+    xmin,
+    xmax,
+    xmintarg,
+    xmaxtarg,
+    nsource,
+    ntarget,
+) where {FT}
+    return (
+        Vector{FT}(range(xmin, xmax, length = nsource)),
+        Vector{FT}(range(xmintarg, xmaxtarg, length = ntarget)),
+    )
+end
+
+function test_single_column(
+    ::Type{FT},
+    xmin,
+    xmax,
+    nsource,
+    ntarget;
+    toler,
+    xmintarg = xmin,
+    xmaxtarg = xmax,
+    extrapolation = Flat(),
+) where {FT}
+    @testset "1D linear interpolation on single column with $FT" begin
+        xsource, xtarget =
+            get_uniform_column_grids(FT, xmin, xmax, xmintarg, xmaxtarg, nsource, ntarget)
+        fsource = sin.(xsource) # function defined on source grid
+        ftarget = zeros(FT, ntarget) # allocated function on target grid
+        interpolate1d!(ftarget, xsource, xtarget, fsource, Linear(), extrapolation)
+        diff = maximum(
+            abs.(ftarget .- sin.(xtarget)) .* (xtarget .≤ xmax) .* (xtarget .≥ xmin),
+        )
+        @test diff ≤ toler
+        # test extrapolation
+        if xmintarg < xmin || xmaxtarg > xmax
+            if extrapolation == Flat()
+                left_boundary_pass = true
+                right_boundary_pass = true
+                for i = 1:length(xtarget)
+                    if xtarget[i] < xmin
+                        left_boundary_pass = ftarget[i] == fsource[1]
+                    end
+                    if xtarget[i] > xmax
+                        right_boundary_pass = ftarget[i] == fsource[end]
+                    end
+                end
+                @testset "testing Flat extrapolation" begin
+                    @test left_boundary_pass
+                    @test right_boundary_pass
+                end
+            end
+        end
+
+    end
+    return nothing
+end
+
+function test_multiple_columns(
+    ::Type{FT},
+    xmin,
+    xmax,
+    nsource,
+    ntarget,
+    nlon,
+    nlat;
+    toler,
+    xmintarg = xmin,
+    xmaxtarg = xmax,
+    extrapolation = Flat(),
+) where {FT}
+    @testset "1D linear interpolation on multiple columns with $FT" begin
+        xsource, xtarget =
+            get_uniform_column_grids(FT, xmin, xmax, xmintarg, xmaxtarg, nsource, ntarget)
+
+        xsourcecols = repeat(xsource, 1, nlon, nlat)
+        xtargetcols = repeat(xtarget, 1, nlon, nlat)
+        fsourcecols = sin.(xsourcecols)
+        ftargetcols = zeros(FT, ntarget, nlon, nlat)
+
+        interpolate1d!(
+            ftargetcols,
+            xsourcecols,
+            xtargetcols,
+            fsourcecols,
+            Linear(),
+            extrapolation,
+        )
+        diff = maximum(abs.(ftargetcols .- sin.(xtargetcols))[:])
+        @test diff ≤ toler
+    end
+    return nothing
+end
+
+# single column linear interpolation tests without extrapolation
+test_single_column(Float32, 0.0, 2 * π, 150, 200, toler = 0.0003)
+test_single_column(Float64, 0.0, 2 * π, 150, 200, toler = 0.0003)
+# single column linear interpolation tests with Flat extrapolation
+test_single_column(
+    Float32,
+    0.0,
+    2 * π,
+    150,
+    200,
+    toler = 0.0003,
+    xmintarg = -1.0,
+    xmaxtarg = 2 * π + 1.0,
+    extrapolation = Flat(),
+)
+test_single_column(
+    Float64,
+    0.0,
+    2 * π,
+    150,
+    200,
+    toler = 0.0003,
+    xmintarg = -1.0,
+    xmaxtarg = 2 * π + 1.0,
+    extrapolation = Flat(),
+)
+# multiple column liner interpolation tests without extrapolation
+test_multiple_columns(Float32, 0.0, 2 * π, 150, 200, 2560, 1280, toler = 0.0003)
+test_multiple_columns(Float64, 0.0, 2 * π, 150, 200, 2560, 1280, toler = 0.0003)