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,84 @@
+using ClimaInterpolations
+using ClimaInterpolations.Interpolation1D
+import .Interpolation1D:
+    interpolate, Linear, get_stencil, interpolate1d!, Flat, LinearExtrapolation
+using Statistics
+using BenchmarkTools
+using CairoMakie
+using Test
+
+FT = Float64
+
+xmin, xmax = FT(0), FT(2 * π)
+nsource, ntarget = 150, 200
+#nsource, ntarget = 15, 20
+xsource = Vector{FT}(range(xmin, xmax, length = nsource))
+#xtarget = Vector{FT}(range(xmin-1, xmax+1, length = ntarget))
+xtarget = Vector{FT}(range(xmin, xmax, length = ntarget))
+
+fsource = sin.(xsource) # function defined on source grid
+ftarget = zeros(FT, ntarget) # allocated function on target grid
+
+#extrapolation = LinearExtrapolation()
+extrapolation = Flat()
+@time interpolate1d!(ftarget, xsource, xtarget, fsource, Linear(), extrapolation)
+@time interpolate1d!(ftarget, xsource, xtarget, fsource, Linear(), extrapolation)
+@allocated interpolate1d!(ftarget, xsource, xtarget, fsource, Linear(), extrapolation)
+
+bm_trial = @benchmark interpolate1d!(
+    $ftarget,
+    $xsource,
+    $xtarget,
+    $fsource,
+    $Linear(),
+    $extrapolation,
+)
+
+@show (nsource, ntarget)
+@show bm_trial
+@show Statistics.median(bm_trial)
+
+fig, ax, plt = scatterlines(xsource, fsource, color = :blue, marker = :star4)
+scatterlines!(ax, xtarget, ftarget, color = :green, marker = :circle)
+save("test.pdf", fig)
+# multiple 1D columns
+nlon, nlat = 2560, 1280
+#nlon, nlat = 5, 7
+xsourcecols = repeat(xsource, 1, nlon, nlat)
+xtargetcols = repeat(xtarget, 1, nlon, nlat)
+fsourcecols = repeat(fsource, 1, nlon, nlat)
+ftargetcols = zeros(ntarget, nlon, nlat)
+
+@time interpolate1d!(
+    ftargetcols,
+    xsourcecols,
+    xtargetcols,
+    fsourcecols,
+    Linear(),
+    extrapolation,
+)
+@time interpolate1d!(
+    ftargetcols,
+    xsourcecols,
+    xtargetcols,
+    fsourcecols,
+    Linear(),
+    extrapolation,
+)
+
+bm_trial_mcols = @benchmark interpolate1d!(
+    $ftargetcols,
+    $xsourcecols,
+    $xtargetcols,
+    $fsourcecols,
+    $Linear(),
+    $extrapolation,
+)
+
+@show bm_trial_mcols
+@show Statistics.median(bm_trial_mcols)
+
+@testset "testing interpolate1d! on multiple columns" begin
+    @test ftargetcols[:, 1, 1] ≈ ftarget[:]
+    @test ftargetcols[:, 25, 87] ≈ ftarget[:]
+end