text:
module ApplyTests

using Test
using Gridap.Arrays
using Gridap.Arrays: ArrayWithCounter, reset_counter!
using FillArrays

a = rand(3,2,4)
test_array(a,a)

a = rand(3,2,4)
a = CartesianIndices(a)
test_array(a,a)

a = rand(3,2)
a = CartesianIndices(a)
c = apply(-,a)
test_array(c,-a)

a = rand(12)
c = apply(-,a)
test_array(c,-a)

a = rand(12)
b = rand(12)
c = apply(-,a,b)
test_array(c,a.-b)

c = apply(Float64,-,a,b)
test_array(c,a.-b)

a = rand(0)
b = rand(0)
c = apply(-,a,b)
test_array(c,a.-b)

a = fill(rand(2,3),12)
b = rand(12)
c = apply(bcast(-),a,b)
test_array(c,[ai.-bi for (ai,bi) in zip(a,b)])

a = fill(rand(2,3),0)
b = rand(0)
c = apply(bcast(-),a,b)
test_array(c,[ai.-bi for (ai,bi) in zip(a,b)])

a = fill(rand(2,3),12)
b = rand(12)
c = apply(bcast(-),a,b)
d = apply(bcast(+),a,c)
e = apply(bcast(*),d,c)
test_array(e,[((ai.-bi).+ai).*(ai.-bi) for (ai,bi) in zip(a,b)])

a = fill(rand(Int,2,3),12)
b = fill(rand(Int,1,3),12)
c = array_caches(a,b)
i = 1
v = getitems!(c,(a,b),i)
@test c == (nothing,nothing)
@test v == (a[i],b[i])

a = fill(rand(Int,2,3),12)
b = fill(rand(Int,1,3),12)
ai = testitem(a)
@test ai == a[1]
ai, bi = testitems(a,b)
@test ai == a[1]
@test bi == b[1]

a = fill(rand(Int,2,3),0)
b = fill(1,0)
ai = testitem(a)
@test ai == Array{Int,2}(undef,0,0)
ai, bi = testitems(a,b)
@test ai == Array{Int,2}(undef,0,0)
@test bi == zero(Int)

a = fill(+,10)
x = rand(10)
y = rand(10)
v = apply(a,x,y)
r = [(xi+yi) for (xi,yi) in zip(x,y)]
test_array(v,r)
v = apply(Float64,a,x,y)
test_array(v,r)

a = Fill(bcast(+),10)
x = [rand(2,3) for i in 1:10]
y = [rand(1,3) for i in 1:10]
v = apply(a,x,y)
r = [(xi.+yi) for (xi,yi) in zip(x,y)]
test_array(v,r)

a = Fill(bcast(+),10)
x = [rand(mod(i-1,3)+1,3) for i in 1:10]
y = [rand(1,3) for i in 1:10]
v = apply(a,x,y)
r = [(xi.+yi) for (xi,yi) in zip(x,y)]
test_array(v,r)

# Test the intermediate results caching mechanism

a = ArrayWithCounter(fill(rand(2,3),12))
b = ArrayWithCounter(rand(12))
c = apply(bcast(-),a,b)
d = apply(bcast(+),a,c)
e = apply(bcast(*),d,c)
r = [ (ai.-bi).*(ai.+(ai.-bi)) for (ai,bi) in zip(a,b)]
cache = array_cache(e)
reset_counter!(a)
reset_counter!(b)
for i in 1:length(e)
  ei = getindex!(cache,e,i)
  ei = getindex!(cache,e,i)
  ei = getindex!(cache,e,i)
end

@test all(a.counter .== 2)
@test all(b.counter .== 1)

l = 10
ai = 1.0
bi = 2.0
a = Fill(ai,l)
b = Fill(bi,l)
c = apply(+,a,b)
r = map(+,a,b)
test_array(c,r)
@test isa(c,Fill)

l = 10
ai = [8 0; 0 4]
a = Fill(ai,l)
c = apply(inv,a)
r = map(inv,a)
test_array(c,r)
@test isa(c,Fill)

l = 10
ai = [8 0; 0 4]
a = fill(ai,l)
c = apply(inv,a)
r = map(inv,a)
test_array(c,r)

l = 0
ai = [8 0; 0 4]
a = fill(ai,l)
c = apply(inv,a)
r = map(inv,a)
test_array(c,r)

ai = [8 0; 0 4]
a = CompressedArray([ai,],Int[])
c = apply(inv,a)
r = map(inv,a)
test_array(c,r)

l = 10
ai = [8, 0]
a = fill(ai,l)
f(ai) = ai[2]-ai[1]
c = apply(f,a)
r = map(f,a)
test_array(c,r)

g(ai) = ai[2]-ai[1]
import Gridap.Arrays: kernel_testitem!
kernel_testitem!(c,::typeof(g),ai) = zero(eltype(ai))
l = 0
ai = [8, 0]
a = fill(ai,l)
c = apply(g,a)
r = map(g,a)
test_array(c,r)

#f = rand(10)
#a = rand(10)
#@test apply(f,a) === f
#
#f = fill(rand(4),10)
#a = rand(10)
#@test apply(f,a) === f
#
#f = fill(rand(4),10)
#g = fill(rand(4),10)
#a = rand(10)
#h = apply_all((f,g),a)
#@test h[1] === f
#@test h[2] === g

#l = 10
#ai = 1.0
#bi = 2.0
#a = fill(ai,l)
#b = fill(bi,l)
#c = apply(+,a,b)
#d = apply(c,a,b)
#@test c == d
#@test typeof(d) == typeof(c)

end # module


text:
const MAX_WRITE_SIZE = Int(typemax(UInt16))
struct BufferedOutput
    action::Function
    waiter::Task
    iob::IOBuffer
    lck::ReentrantLock

    function BufferedOutput(action::Function; delay::Float64=0.05, bytelimit::Int=1024*8, maxbuff::Int=MAX_WRITE_SIZE)
        lck = ReentrantLock()
        iob = IOBuffer(;maxsize=maxbuff)
        waiter = @async begin
            while isopen(iob)
                timedwait(delay) do
                    lock(lck) do
                        position(iob) >= bytelimit
                    end
                end
                lock(lck) do
                    if position(iob) > 0
                        action(take!(iob))
                    end
                end
            end
        end
        new(action, waiter, iob, lck)
    end
end

function write(buff::BufferedOutput, data)
    towrite = length(data)
    pos = 1
    while pos <= towrite
        (pos == 1) || sleep(0.05)
        endpos = min(pos+buff.iob.maxsize-1, towrite)
        pos += lock(buff.lck) do
            write(buff.iob, view(data, pos:endpos))
        end
    end
end

close(buff::BufferedOutput) = close(buff.iob)
function flush(buff::BufferedOutput)
    lock(buff.lck) do
        if position(buff.iob) > 0
            (buff.action)(take!(buff.iob))
        end
    end
end

writefully(io::T, iob::IOBuffer) where {T <: IO} = writefully(io, take!(iob))
function writefully(io::T, b::Vector{UInt8}, nbytes::Int=length(b)) where {T <: IO}
    pos = write(io, b)
    while pos < nbytes
        pos += write(io, view(b, pos:nbytes))
    end
    nbytes
end

ext: jl


text:
# index.jl getindex and setindex methods for NamedArray
## (c) 2013--2020 David A. van Leeuwen

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

import Base: getindex, setindex!

## AbstractArray Interface, integers have precedence over everything else
@inline Base.@propagate_inbounds function Base.getindex(n::NamedArray{T, N, AT, DT},
    i::Int) where {T, N, AT, DT}
    getindex(n.array, i)
end

@inline Base.@propagate_inbounds function Base.getindex(n::NamedArray{T, N, AT, DT},
    I::Vararg{Int, N}) where {T, N, AT, DT}
    getindex(n.array, I...)
end

@inline Base.@propagate_inbounds function Base.setindex!(n::NamedArray{T, N, AT, DT},
     v, i::Int) where {T, N, AT, DT}
     setindex!(n.array, v, i::Int)
end

@inline Base.@propagate_inbounds function Base.setindex!(n::NamedArray{T, N, AT, DT},
     v, I::Vararg{Int, N}) where {T, N, AT, DT}
     setindex!(n.array, v, I...)
end

## optional methods
Base.IndexStyle(n::NamedArray) = IndexStyle(n.array)

## Ambiguity
#getindex(n::NamedArray{T, 1, AT, DT}, i::Int64) where {T, AT, DT} = getindex(n.array, i)
#setindex!(n::NamedArray{T, 1, AT, DT}, v::Any, i::Int64) where {T, AT, DT} = setindex!(n.array, v, i)

function flattenednames(n::NamedArray)
    L = length(n) # elements in array
    cols = Array[]
    factor = 1
    for d in 1:ndims(n)
        nlevels = size(n, d)
        nrep = L ÷ (nlevels * factor)
        data = repeat(vcat([fill(x, factor) for x in names(n, d)]...), nrep)
        push!(cols, data)
        factor *= nlevels
    end
    return collect(zip(cols...))
end

## from subarray.jl
@inline Base.@propagate_inbounds getindex(n::NamedVector, ::Colon) = n
getindex(n::NamedArray, ::Colon) = NamedArray(n.array[:], (flattenednames(n),) , (tuple(dimnames(n)...),))

## special 0-dimensional case
## getindex{T}(n::NamedArray{T,0}, i::Real) = getindex(n.array, i)

getindex(n::NamedArray{T, N, AT, DT}, I::Vararg{Any,N}) where {T, N, AT, DT} = namedgetindex(n, map((d,i)->indices(d, i), n.dicts, I)...)

Base.view(n::NamedArray{T,N}, I::Vararg{Union{AbstractArray,Colon,Real},N}) where {T,N} = namedgetindex(n, map((d,i)->indices(d, i), n.dicts, I)...; useview=true)
Base.view(n::NamedArray{T,N}, I::Vararg{Any,N}) where {T,N} = namedgetindex(n, map((d,i)->indices(d, i), n.dicts, I)...; useview=true)

## indices computes numeric indices from named or other

## indices(::Associative, index) converts any type `index` to Integer

## single index
indices(dict::AbstractDict{K,V}, i::Integer) where {K, V<:Integer} = i ## integer index takes precedence
## indices(dict::Associative{K,V}, i::K) where {K<:Real, V<:Integer} = dict[i]
## indices(dict::Associative{K,V}, i::Real) where {K, V<:Integer} = Base.to_index(i)
indices(dict::AbstractDict{K,V}, i::K) where {K, V<:Integer} = dict[i]
indices(dict::AbstractDict{K,V}, i::Name{IK}) where IK <: K where {K, V<:Integer} = dict[i.name]

## ambiguity if dict key is CartesionIndex, this should never happen
# indices(dict::AbstractDict{K,V}, i::K) where {K<:CartesianIndex,V<:Integer} = dict[i]
@inline indices(dict::AbstractDict, ci::CartesianIndex) = ci
@inline indices(dict::AbstractDict{Any,V}, ci::CartesianIndex) where {V<:Integer} = ci

## multiple indices
## the following two lines are partly because of ambiguity
#indices(dict::Associative{T,V}, i::AbstractArray{T}) where {T<:Integer,V<:Integer} = [dict[k] for k in i]
#indices(dict::Associative{T,V}, i::AbstractArray{T}) where {T<:Real,V<:Integer} = [dict[k] for k in i]

indices(dict::AbstractDict{K,V}, i::AbstractArray) where {K,V<:Integer} = [indices(dict, k) for k in i]
#indices(dict::Associative{K,V}, i::AbstractArray{T}) where {T<:Integer,K,V<:Integer} = i
#indices(dict::Associative{K,V}, i::AbstractArray{K}) where {K,V<:Integer} = [dict[k] for k in i]
#indices(dict::Associative{K,V}, i::AbstractArray{Name{K}}) where {K, V<:Integer} = [dict[k.name] for k in i]
## in 0.4, we need to take care of : ourselves it seems
indices(dict::AbstractDict{K,V}, ::Colon) where {K,V<:Integer} = collect(1:length(dict))

## negation
indices(dict::AbstractDict{K,V}, i::K) where {K<:Not,V<:Integer} = dict[i]
indices(dict::AbstractDict, i::Not) = setdiff(1:length(dict), indices(dict, i.skip))

## namedgetindex collects the elements from the array, and takes care of the index names
## `index` is an integer now, or an array of integers, or a cartesianindex
## and has been computed by `indices()`

## Simple scalar indexing
@inline namedgetindex(n::NamedArray, I::Vararg{Integer,N}) where {N} = getindex(n.array, I...)

@inline dimkeepingtype(x) = false
@inline dimkeepingtype(x::AbstractArray) = true
@inline dimkeepingtype(x::AbstractRange) = true
@inline dimkeepingtype(x::BitArray) = true

## Slices etc.
@inline Base.@propagate_inbounds namedgetindex(n::NamedArray, index::CartesianIndex) = getindex(n.array, index)
Base.@propagate_inbounds function namedgetindex(n::NamedArray, index...; useview::Bool=false)
    if useview
        a = view(n.array, index...)
    else
        a = getindex(n.array, index...)
    end
    if !isa(a, AbstractArray)
        return(a)
    end
    N = length(index)
    keeping = filter(i -> dimkeepingtype(index[i]), 1:N)
    if ndims(a) < length(keeping) ## || length(dims) == 1 && ndims(n) > 1
        warn("Dropped names for ", typeof(n.array), " with index ", index)
        return a;               # number of dimension changed, this should not happen
    end
    newnames = Array{Any}(undef, ndims(a))
    newdimnames = Array{eltype(n.dimnames)}(undef, ndims(a))
    count = 0
    for d in keeping
        if ndims(index[d]) > 1
            ## take over the names of the index for this dimension
            for (name, dimname) in zip(defaultnames(index[d]), dimnames(index[d]))
                count += 1
                newnames[count] = name
                newdimnames[count] = Symbol(string(n.dimnames[d], "_", dimname))
            end
        else
            count += 1
            newnames[count] = names(n, d)[index[d]]
            newdimnames[count] = n.dimnames[d]
        end
    end
    return NamedArray(a, tuple(newnames...), tuple(newdimnames...))
end

## work out n[:A => "1", :C => "5"]
function indices(n::NamedArray, I::Pair...)
    dict = Dict(I...)
    Set(keys(dict)) ⊆ Set(n.dimnames) || error("Dimension name mismatch")
    result = Vector{Union{Int,Colon}}(undef, ndims(n))
    fill!(result, :) ## unspecified dimensions act as colon
    for (i, dim) in enumerate(n.dimnames)
        if dim in keys(dict)
            result[i] = n.dicts[i][dict[dim]]
        end
    end
    return result
end

@inline Base.@propagate_inbounds getindex(n::NamedArray, I::Pair...) = getindex(n.array, indices(n, I...)...)
## 0.6 ambiguity
@inline Base.@propagate_inbounds getindex(n::NamedVector, I::CartesianIndex{1}) = getindex(n.array, I)
@inline Base.@propagate_inbounds getindex(n::NamedArray{T,N}, I::CartesianIndex{N}) where {T,N} = getindex(n.array, I)

## Setindex is remarkably more simple than getindex.  I wonder why...

## This takes care of most cases
@inline Base.@propagate_inbounds function setindex!(n::NamedArray{T,N}, x, I::Vararg{Any,N}) where {T,N}
    II = map((d,i)->indices(d, i), n.dicts, I)
    n.array[II...] = x
end

## assignment via n[:B => "two"] = [...]
@inline Base.@propagate_inbounds setindex!(n::NamedArray, x, I::Vararg{Pair}) = setindex!(n.array, x, indices(n, I...)...)


text:
using Ipopt
using MathOptInterface

const MOI = MathOptInterface

function moi_solve(
    optimizer::MOI.AbstractOptimizer,
    nlp::ExaPF.AbstractNLPEvaluator,
)
    block_data = MOI.NLPBlockData(nlp)
    u♭, u♯ = ExaPF.bounds(nlp, ExaPF.Variables())
    u0 = ExaPF.initial(nlp)
    n = ExaPF.n_variables(nlp)
    vars = MOI.add_variables(optimizer, n)
    # Set bounds and initial values
    for i in 1:n
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(vars[i]),
            MOI.LessThan(u♯[i])
        )
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(vars[i]),
            MOI.GreaterThan(u♭[i])
        )
        MOI.set(optimizer, MOI.VariablePrimalStart(), vars[i], u0[i])
    end
    MOI.set(optimizer, MOI.NLPBlock(), block_data)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(optimizer)
    x_opt = [MOI.get(optimizer, MOI.VariablePrimal(), v) for v in vars]
    solution = (
        minimum=MOI.get(optimizer, MOI.ObjectiveValue()),
        minimizer=x_opt
    )
    MOI.empty!(optimizer)
    return solution
end

@testset "MOI wrapper" begin
    CASE57_SOLUTION = [
        1.0260825400262428,
        0.0,
        0.6,
        0.0,
        8.603379123083476,
        0.0,
        1.3982297290334753,
        1.016218932360429,
        1.009467718490475,
        1.027857273911734,
        1.06,
        0.9962215167521776,
        0.9965415804936182
    ]

    datafile = joinpath(dirname(@__FILE__), "..", "data", "case57.m")
    nlp = ExaPF.ReducedSpaceEvaluator(datafile)
    optimizer = Ipopt.Optimizer(
        print_level=0,
        limited_memory_max_history=50,
        hessian_approximation="limited-memory",
        tol=1e-2,
    )
    solution = moi_solve(optimizer, nlp)
    @test solution.minimum ≈ 3.7589338e+04
    @test solution.minimizer ≈ CASE57_SOLUTION atol=1e-6
end


text:
#here all the different input structures for the different calculations go
#TODO should we also create a config file for each job with stuff like server etc? and other config things,
#      which if not supplied could contain the default stuff?
"""
Represents a full DFT job with multiple input files and calculations.
"""
mutable struct DFJob
    id           ::Int #REVIEW: For now ID does nothing but I can see it doing things later
    name         ::String
    structure    ::AbstractStructure
    inputs       ::Vector{DFInput}
    local_dir    ::String
    server       ::String
    server_dir   ::String
    header       ::Vector{String}
    metadata     ::Dict
    function DFJob(name, structure, calculations, local_dir, server, server_dir, header = getdefault_jobheader())
        if local_dir != ""
            local_dir = local_dir
        end

        if server_dir != ""
            server_dir = server_dir
        end
        if !isabspath(local_dir)
            local_dir = abspath(local_dir)
        end
        return new(0, name, structure, calculations, local_dir, server, server_dir, header, Dict())
    end
end

#TODO implement abinit
function DFJob(job_name, local_dir, structure::AbstractStructure, calculations::Vector, common_flags...;
                    server=getdefault_server(),
                    server_dir="",
                    package=QE,
                    bin_dir=joinpath("usr","local","bin"),
                    pseudoset=:default,
                    pseudospecifier="",
                    header=getdefault_jobheader())

    @assert package==QE "Only implemented for Quantum Espresso!"

    job_calcs = DFInput[]
    if typeof(common_flags) != Dict
        common_flags = Dict(common_flags)
    end
    bin_dir = bin_dir
    for (calc, (excs, data)) in calculations
        if !isa(data, SymAnyDict)
            data = SymAnyDict(data)
        end
        calc_ = typeof(calc) == String ? Symbol(calc) : calc
        if in(calc_, [:vc_relax, :relax, :scf])
            k_points = get(data, :kpoints, [1, 1, 1, 0, 0, 0])
            k_option = :automatic
        elseif calc_ == :nscf
            k_points = kgrid(get(data, :kpoints, [1, 1, 1])[1:3]..., QE)
            k_option = :crystal
        elseif calc_ == :bands
            k_points = get(data, :kpoints, [[0., 0., 0., 1.]])
            num_k = 0.0
            for point in k_points
                num_k += point[4]
            end
            if num_k > 100.
                if !haskey(data, :flags)
                    data[:flags] = Pair{Symbol, Any}[]
                end
                push!(data[:flags], :verbosity => "high")
            end
            k_option = :crystal_b
        end
        flags  = convert(Vector{Pair{Symbol, Any}}, get(data, :flags, Pair{Symbol, Any}[]))
        if excs[2].exec == "pw.x"
            push!(flags, :calculation => "$(string(calc_))")
            datablocks = [InputData(:k_points, k_option, k_points)]
        else
            datablocks =  InputData[]
        end
        input_ = DFInput{package}(string(calc_), local_dir,
                         Dict{Symbol, Any}(),
                         datablocks, excs, true)
        setflags!(input_, common_flags..., print=false) #This could be changed to use suppressor
        setflags!(input_, flags..., print=false)
        push!(job_calcs, input_)
    end
    out = DFJob(job_name, structure, job_calcs, local_dir, server, server_dir, header)
    setatoms!(out, structure.atoms, pseudoset = pseudoset, pseudospecifier= pseudospecifier)
    if !haskey(common_flags, :ecutwfc)
        @info "No :ecutwfc specified in the flags, determining minimum cutoff from pseudos."
        setcutoffs!(out)
    end
    return out
end

function DFJob(job_name, local_dir, ciffile::String, calculations::Vector, args...; kwargs...)
    structure = Structure(ciffile, name=job_name)
    return DFJob(job_name, local_dir, structure, calculations, args... ; kwargs...)
end

function DFJob(job::DFJob, flagstoset...; cell_=copy(cell(job)), atoms_=copy(atoms(job)), name=job.name,
                                          server_dir = job.server_dir,
                                          local_dir  = job.local_dir,
                                          pseudoset  = nothing,
                                          pseudospecifier = "")
    newjob = deepcopy(job)

    setcell!(newjob, cell_)
    if pseudoset == nothing
        pseudoset, specifier = getpseudoset(job.structure.atoms[1])
        specifier = pseudospecifier == nothing ? specifier : pseudospecifier
        setatoms!(newjob, atoms_, pseudoset = pseudoset, pseudospecifier=specifier)
    else
        setatoms!(newjob, atoms_, pseudoset = pseudoset, pseudospecifier= pseudospecifier)
    end
    setserverdir!(newjob, server_dir)
    setlocaldir!(newjob, local_dir)
    newjob.name = name

    setflags!(newjob, flagstoset..., print=false)
    return newjob
end

"""
    DFJob(job_dir::String, T=Float64; job_fuzzy = "job", new_job_name=nothing, new_local_dir=nothing, server=getdefault_server(),server_dir="")

Loads and returns a local DFJob. If local_dir is not specified the job directory will be registered as the local one.
"""
function DFJob(job_dir::String, T=Float64;
                  job_fuzzy     = "job",
                  new_job_name  = "",
                  new_local_dir = nothing,
                  server        = getdefault_server(),
                  server_dir    = "")
    name, header, inputs, structure = read_job_inputs(joinpath(job_dir, searchdir(job_dir, job_fuzzy)[1]))
    j_name = isempty(new_job_name) ? name : new_job_name
    structure_name = split(j_name, "_")[1]
    structure.name = structure_name

    if new_local_dir != nothing
        return DFJob(j_name, structure, inputs, new_local_dir, server, server_dir, header)
    else
        return DFJob(j_name, structure, inputs, job_dir, server, server_dir, header)
    end
end

"""
    DFJob(server_dir::String, local_dir::String, server=getdefault_server(); job_fuzzy="*job*", new_job_name="")

Pulls a server job to local directory and then loads it. A fuzzy search for the job file will be performed and the found input files will be pulled.
"""
function DFJob(server_dir::String, local_dir::String, server = getdefault_server();
                         job_fuzzy    = "*job*",
                         new_job_name = "")

    pulljob(server, server_dir, local_dir)
    return DFJob(local_dir, server=server, server_dir=server_dir, new_job_name=new_job_name)
end

name(job) = job.name
#-------------------BEGINNING GENERAL SECTION-------------#
scriptpath(job::DFJob)       = joinpath(job.local_dir, "job.tt")
starttime(job::DFJob)        = mtime(scriptpath(job))

runslocal(job::DFJob)        = job.server        =="localhost"
structure(job::DFJob)        = job.structure
iswannierjob(job::DFJob)     = any(x->package(x) == Wannier90, inputs(job)) && any(x->isnscfcalc(x), inputs(job))
getnscfcalc(job::DFJob)      = getfirst(x -> isnscfcalc(x), inputs(job))
cell(job::DFJob)             = cell(structure(job))

input(job::DFJob, n::String) = getfirst(x -> occursin(n, name(x)), inputs(job))
inputs(job::DFJob)           = job.inputs

"""
    inputs(job::DFJob, names::Vector)

Returns an array of the inputs that match the names.
"""
inputs(job::DFJob, names::Vector, fuzzy=true) = fuzzy ? filter(x -> any(occursin.(names, name(x))), inputs(job)) : input.(job, names)
inputs(job::DFJob, n::String, fuzzy=true) = inputs(job, [n], fuzzy)
inputs(job::DFJob, package_::Package) = filter(x->package(x)==package_, inputs(job))
inpath(job::DFJob, n) = inpath(input(job,n))
outpath(job::DFJob, n) = outpath(input(job,n))

"Runs some checks on the set flags for the inputs in the job, and sets metadata (:prefix, :outdir etc) related flags to the correct ones. It also checks whether flags in the various inputs are allowed and set to the correct types."
function sanitizeflags!(job::DFJob)
    setflags!(job, :prefix => "$(job.name)", print=false)
    if iswannierjob(job)
	    nscfcalc = getnscfcalc(job)
	    if package(nscfcalc) == QE && hasflag(nscfcalc, :nbnd)
	        setflags!(job, :num_bands => nscfcalc[:nbnd], print=false)
        elseif package(nscfcalc) == Elk
	        setflags!(job, :num_bands => length(nscfcalc[:wann_bands]))
	        nscfcalc[:wann_projections] = projections_string.(unique(filter(x->!isempty(projections(x)), atoms(job))))
	        nscfcalc[:elk2wan_tasks]    = ["602", "604"]
	        nscfcalc[:wann_seedname]    = Symbol(name(job))
	        if job[:wannier_plot] == true
		        push!(nscfcalc[:elk2wan_tasks], "605")
	        end
        end
    end
    for i in filter(x -> package(x) == QE, inputs(job))
	    set_hubbard_flags!(i, job.structure)
	    set_starting_magnetization_flags!(i, job.structure)
    end
    sanitize_pseudos!(job)
    sanitizeflags!.(inputs(job))
end

function sanitize_pseudos!(job::DFJob)
	all_pseudos = pseudo.(atoms(job))
	uni_dirs    = unique(map(x->x.dir, all_pseudos))
	uni_pseudos = unique(all_pseudos)
	pseudo_dir  = length(uni_dirs) == 1 ? uni_dirs[1] : job.local_dir
	if length(uni_dirs) > 1
		@info "Found pseudos in multiple directories, copying them to job directory"
		for pseudo in uni_pseudos
			cp(path(pseudo), joinpath(job.local_dir, pseudo.name), force=true)
		end
	end
	for p in all_pseudos
		p.dir = pseudo_dir
	end
	for i in filter(x -> package(x) == QE, inputs(job))
		setflags!(i, :pseudo_dir => pseudo_dir; print=false)
	end
end

function sanitize_magnetization!(job::DFJob)
    if !any(x->package(x) == QE, inputs(job))
        return
    end
    magnetic_ats = filter(ismagnetic, atoms(job))
    magnetic_elements = map(element, magnetic_ats)
    for e in magnetic_elements
        magnetizations = Vec3[]
        dftus          = DFTU[]
        for a in magnetic_ats
            if e == element(a)
                magid = findfirst(x -> x == magnetization(a), magnetizations)
                dftuid = findfirst(x-> x == dftu(a), dftus)
                if magid === nothing
                    push!(magnetizations, magnetization(a))
                    tid1 = length(magnetizations)
                else
                    tid1 = magid
                end
                if dftuid === nothing
                    push!(dftus, dftu(a))
                    tid2 = length(dftus)
                else
                    tid2 = dftuid
                end
                atid = max(tid1, tid2)
                a.name = Symbol(string(e.symbol)*"$atid")
            end
        end
    end
end

function sanitize_projections!(job::DFJob)
    if !any(x->!isempty(projections(x)), atoms(job))
        return
    end
    uats = unique(atoms(job))
    projs = unique([name(at) => [p.orb.name for p in projections(at)] for at in uats])
    setprojections!(job, projs...)
end

"Checks the last created output file for a certain job."
function runninginput(job::DFJob)
    @assert job.server == "localhost" "Intended use for now is locally."
    t = mtime(scriptpath(job))
    for i in reverse(inputs(job))
        p = outpath(i)
        if ispath(p) && mtime(p) > t
            return i
        end
    end
end

"Finds the input corresponding to the name and returns the full output path."
outpath(job::DFJob, n::String) = outpath(input(job,n))

Base.joinpath(job::DFJob, n::AbstractString) = joinpath(job.local_dir, n)

runslocal_assert(job::DFJob) =
    @assert runslocal(job) "This only works if the job runs on `localhost`."

function Base.pop!(job::DFJob, name::String)
    i = findfirst(x -> x.name == name, job.inputs)
    if i === nothing
        error("No input with name $name found.")
    end
    out = job.inputs[i]
    deleteat!(job.inputs, i)
    return out
end


text:
module ColorfulCodeGen

using InteractiveUtils: gen_call_with_extracted_types_and_kwargs,
    code_lowered, code_typed, code_warntype, code_llvm, code_native
using Markdown: MD
function open_reader(cmd, proc_out)
    process = open(cmd, proc_out; write=true)
    return (process.in, process)
end

using Base.Meta: show_sexpr

mutable struct PygmentizeConfig
    command::Cmd
    options::Dict{Any, Cmd}
end

code_intel(args...) = code_native(args...; syntax=:intel)

const PYGMENTIZE = PygmentizeConfig(
    `pygmentize -f terminal256`,
    Dict(
        code_lowered  => `-l julia`,
        code_typed    => `-l julia`,
        code_warntype => `-l julia`,
        code_llvm     => `-l llvm`,
        code_native   => `-l asm`,
        code_intel    => `-l nasm`,
        show_sexpr    => `-l scheme`,
        Expr          => `-l julia`,
        MD            => `-l md`,
    ),
)

struct ColoredCode
    text::String
    optionkey
end

function Base.show(io::IO, ::MIME"text/plain", cc::ColoredCode)
    if get(io, :compact, false) || get(io, :typeinfo, Any) === ColoredCode
        print(io, "ColoredCode(")
        text = repr(cc.text)
        if get(io, :limit, true) && length(text) > 20
            print(io, first(text, 20))
            printstyled(io, '…'; color=:light_black)
            print(io, '"')
        else
            print(io, text)
        end
        print(io, ", ", cc.optionkey, ")")
        return
    end
    highlighter = get_command(cc.optionkey)
    proc_in, process = open_reader(highlighter, io)
    try
        print(proc_in, cc.text)
    finally
        close(proc_in)
        wait(process)
    end
    return
end

function get_command(key, config::PygmentizeConfig = PYGMENTIZE)
    option = get(config.options, key, nothing)
    option === nothing && error("Highlighter not found for: ", key)
    return `$(config.command) $option`
end

mimefor(_) = MIME"text/plain"()
mimefor(::MD) = MIME"text/markdown"()

function highlight(x)
    text = sprint(show, mimefor(x), x)
    return ColoredCode(text, typeof(x))
end

const _with_io = (:code_warntype, :code_llvm, :code_native, :code_intel)
const _no_io = (:code_typed, :code_lowered)

for fname in (_with_io..., :show_sexpr)
    colored_name = Symbol("c$fname")
    @eval begin
        @doc $("""
            $(colored_name)([io,] args...; kwargs...)

        Colored version of `$(fname)([io,] args...; kwargs...)`.
        """)
        function $(colored_name)(args...; kwargs...)
            text = sprint() do io
                $fname(io, args...; kwargs...)
            end
            return ColoredCode(text, $fname)
        end
    end
end

for fname in _no_io
    colored_name = Symbol("c$fname")
    @eval begin
        @doc $("""
            $(colored_name)([io,] args...; kwargs...)

        Colored version of `$(fname)(args...; kwargs...)`.
        """)
        function $(colored_name)(args...; kwargs...)
            results = $fname(args...; kwargs...)
            result = length(results) == 1 ? results[1] : results
            # ^- from interactiveutil.jl
            text = sprint(show, MIME("text/plain"), result)
            return ColoredCode(text, $fname)
        end
    end
end

for fname in (_with_io..., _no_io...)
    colored_name = Symbol("c$fname")
    @eval begin
        $(colored_name)(io::IO, args...; kwargs...) =
            show(io, MIME"text/plain"(), $(colored_name)(args...; kwargs...))
        macro ($colored_name)(ex0...)
            quoted = $(Expr(:quote, colored_name))
            gen_call_with_extracted_types_and_kwargs(__module__, quoted, ex0)
        end

        export $colored_name
        export $(Symbol("@$colored_name"))
    end
end

dehygiene(expr::Expr) = Expr(expr.head, dehygiene.(expr.args)...)
dehygiene(x::Any) = x  # LineNumberNode etc.

function dehygiene(sym::Symbol)
    str = string(sym)
    if startswith(str, "#")
        return Symbol(replace(str, "#" => "😷"))
    else
        return sym
    end
end

"""
    @cmacroexpand ex

Print syntax highlighted expression of what `@macroexpand ex` returns.
"""
macro cmacroexpand(ex)
    return :(highlight(dehygiene(macroexpand($__module__, $(QuoteNode(ex))))))
end

export @cmacroexpand

macro cmacroexpand1(ex)
    :(highlight(dehygiene(macroexpand($__module__, $(QuoteNode(ex));
                                      recursive = false))))
end

export @cmacroexpand1

export @cshow_sexpr

"""
    @cshow_sexpr expression

Print S-expression of Julia `expression`; it's a colored version of
`Meta.show_sexpr(:(expression))`.
"""
macro cshow_sexpr(ex)
    :($cshow_sexpr($(QuoteNode(ex))))
end

export @clowered

"""
    @clowered expression

Print lowered form of `expression`; it's a colored version of
`Meta.lower(Main, :(expression))`.
"""
macro clowered(ex)
    quote
        let ex = Meta.lower($__module__, $(QuoteNode(ex)))
            if ex isa $Expr && length(ex.args) == 1 &&
                    ex.args[1] isa $(Core.CodeInfo)
                ex = ex.args[1]
            end
            $highlight(ex)
        end
    end
end

end # module


text:
module KernelGoodnessOfFit

# exports
export
    # kernels
    Kernel, kernel, GaussianRBF,

    # ksd
    u_p,
    Sᵤ_estimate,
    KSDTest,

    # bootstrap
    degenerate_bootstrap,

    # Finite Set Stein Discrepancy
    ξ,
    fssd,
    compute_Ξ,
    τ_from_Ξ,
    Σₚ,
    σ²_H₁,
    fssd_H₁_opt_factor,
    FSSDTest,
    FSSDopt,
    FSSDrand,

    # core
    perform,
    pvalue

include("core.jl")
include("bootstrap.jl")
include("kernels.jl")
include("ksd.jl")
include("fssd.jl")

end # module


text:
using Peacock
function Chern_number(solver::Peacock.FDFD.Solver, polar::Polarisation; Nkx = 4, Nky = 4, bands = 1:4)
    a1 = solver.basis.a1
    a2 = solver.basis.a2
    Px = 1     # lattice constant along x
    Py = 1     # lattice constant along y
    N_max = bands[end]  # maximum band number ***

    Pkx = 2pi/Px/a2[2]
    Pky = 2pi/Py/a2[2]
    deltakx = Pkx/Nkx
    deltaky = Pky/Nky
    NK = Nkx*Nky
    kcx = zeros(NK, 4)
    kcy = zeros(NK, 4)
    fieldTot = zeros(N_max, 1)
    fieldtemp = zeros(ComplexF64, N_max, 1)

    for m=1:Nkx
        for n=1:Nky

            index = (n-1)*Nkx+m
            kcxg = [(m-1)*deltakx,(m)*deltakx,(m)*deltakx,(m-1)*deltakx]
            kcyg = [(n-1)*deltaky,(n-1)*deltaky,(n)*deltaky,(n)*deltaky]
            kcx[index, :] = kcxg * a2[2] + kcyg * a1[2]
            kcy[index, :] = kcyg * a1[1] - a2[1]*kcxg  # Rotation coordinates

            modes = Peacock.FDFD.solve(solver, [kcx[index, 1], kcy[index, 1]], polar, bands = 1:N_max)
            field1 = [get_field_FDFD(modes[i].data, modes[i].basis, k0=modes[i].k0) for i in 1:N_max]

            modes = Peacock.FDFD.solve(solver, [kcx[index, 2], kcy[index, 2]], polar, bands = 1:N_max)
            field2 = [get_field_FDFD(modes[i].data, modes[i].basis, k0=modes[i].k0) for i in 1:N_max]

            modes = Peacock.FDFD.solve(solver, [kcx[index, 3], kcy[index, 3]], polar, bands = 1:N_max)
            field3 = [get_field_FDFD(modes[i].data, modes[i].basis, k0=modes[i].k0) for i in 1:N_max]

            modes = Peacock.FDFD.solve(solver, [kcx[index, 4], kcy[index, 4]], polar, bands = 1:N_max)
            field4 = [get_field_FDFD(modes[i].data, modes[i].basis, k0=modes[i].k0) for i in 1:N_max]

            weighting = transpose(modes[1].weighting)

            for mm = 1:N_max
                temp = field1[mm][:]' .* weighting * field2[mm][:] / abs.(field1[mm][:]' .* weighting * field2[mm][:]) *
                    field2[mm][:]' .* weighting * field3[mm][:] / abs.(field2[mm][:]' .* weighting * field3[mm][:]) *
                    field3[mm][:]' .* weighting * field4[mm][:] / abs.(field3[mm][:]' .* weighting * field4[mm][:]) *
                    field4[mm][:]' .* weighting * field1[mm][:] / abs.(field4[mm][:]' .* weighting * field1[mm][:])

                temp = field1[mm][:]' .* weighting * field2[mm][:] * field2[mm][:]' .* weighting * field3[mm][:] * field3[mm][:]' .* weighting * field4[mm][:] * field4[mm][:]' .* weighting * field1[mm][:] #不归一化对结果没影响
                fieldtemp[mm] = temp[1]
                fieldTot[mm] = fieldTot[mm] + imag(log(fieldtemp[mm]))
            end
        end
    end
    return round.(fieldTot[:]/2/pi)
end

function Chern_number(epf::Function, muf::Function, a1::Array{<:Real,1}, a2::Array{<:Real,1}, d1::Real, d2::Real, polar::Polarisation; Nkx = 4, Nky = 4, bands=1:4)
    geometry = Geometry(epf, muf, a1, a2, d1, d2, polar)
    solver = Peacock.FDFD.Solver(geometry, 2*[d1, d2])
    return Chern_number(solver, polar, Nkx=Nkx, Nky=Nky, bands = bands)
end


text:
using BQCESubroutine
using LinearAlgebra
using YaoLocations: YaoLocations, plain, Locations, CtrlLocations, merge_locations
using Test

@testset "broutine" begin
    include("broutine.jl")
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to printing a table to a string.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Print table to string" begin
    data     = rand(10,4)
    header_m = [1 2 3 4]
    header_v = [1 2 3 4]

    expected = sprint(pretty_table, data)
    result   = pretty_table(String, data)
    @test result == expected

    expected = sprint((io,data)->pretty_table(io, data, header_m), data)
    result   = pretty_table(String, data, header_m)
    @test result == expected

    expected = sprint((io,data)->pretty_table(io, data, header_v), data)
    result   = pretty_table(String, data, header_v)
    @test result == expected

    dict = Dict(:a => 1,
                :b => 2,
                :c => 3,
                :d => 4)

    expected = sprint(pretty_table, data)
    result   = pretty_table(String, data)
    @test result == expected
end


text:
mutable struct Cat{T} where T <: AbstractString
    id::UInt
    name::T
    food::AbstractVector{T}
end
function FindID(Winter::Cat{T}) where T <: AbstractString
    return Winter.id
end

max_stars_repo_name: Rratic/JuliaGuide

max_issues_repo_name: Rratic/JuliaGuide

max_forks_repo_name: Rratic/JuliaGuide


text:
const DSL_STATIC_ANNOTATION = :static
const DSL_ARG_GRAD_ANNOTATION = :grad
const DSL_RET_GRAD_ANNOTATION = :grad
const DSL_TRACK_DIFFS_ANNOTATION = :diffs
const DSL_NO_JULIA_CACHE_ANNOTATION = :nojuliacache

struct Argument
    name::Symbol
    typ::Union{Symbol,Expr}
    annotations::Set{Symbol}
end

Argument(name, typ) = Argument(name, typ, Set{Symbol}())

function parse_annotations(annotations_expr)
    annotations = Set{Symbol}()
    if isa(annotations_expr, Symbol)
        push!(annotations, annotations_expr)
    elseif isa(annotations_expr, Expr) && annotations_expr.head == :tuple
        for annotation in annotations_expr.args
            if !isa(annotation, Symbol)
                error("syntax error in annotations_expr at $annotation")
            else
                push!(annotations, annotation)
            end
        end
    else
        error("syntax error in annotations at $annotations")
    end
    annotations
end

function parse_arg(expr)
    if isa(expr, Symbol)
        # x
        arg = Argument(expr, :Any)
    elseif isa(expr, Expr) && expr.head == :(::)
        # x::Int
        arg = Argument(expr.args[1], expr.args[2])
    elseif isa(expr, Expr) && expr.head == :call
        # (grad,foo)(x::Int)
        annotations_expr = expr.args[1]
        sub_arg = parse_arg(expr.args[2])
        annotations = parse_annotations(annotations_expr)
        arg = Argument(sub_arg.name, sub_arg.typ, annotations)
    else
        dump(expr)
        error("syntax error in gen function argument at $expr")
    end
    arg
end

include("dynamic.jl")
include("static.jl")

function parse_gen_function(ast, annotations)
    if ast.head != :function
        error("syntax error at $ast in $(ast.head)")
    end
    if length(ast.args) != 2
        error("syntax error at $ast in $(ast.args)")
    end
    signature = ast.args[1]
    body = ast.args[2]
    if signature.head == :(::)
        (call_signature, return_type) = signature.args
    elseif signature.head == :call
        (call_signature, return_type) = (signature, :Any)
    else
        error("syntax error at $(signature)")
    end
    name = call_signature.args[1]
    args = map(parse_arg, call_signature.args[2:end])
    static = DSL_STATIC_ANNOTATION in annotations
    if static
        make_static_gen_function(name, args, body, return_type, annotations)
    else
        make_dynamic_gen_function(name, args, body, return_type, annotations)
    end
end

macro gen(annotations_expr, ast)

    # parse the annotations
    annotations = parse_annotations(annotations_expr)

    # parse the function definition
    parse_gen_function(ast, annotations)
end

macro gen(ast)
    parse_gen_function(ast, Set{Symbol}())
end

max_stars_repo_name: dendisuhubdy/Gen

max_issues_repo_name: kumsh/Gen

max_forks_repo_name: kumsh/Gen


text:
module Interpolation
using CLIMA
using MPI
import GaussQuadrature
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.Mesh.Geometry
using CLIMA.Mesh.Elements
using LinearAlgebra
using StaticArrays
using NCDatasets
#-------------------
using CUDAnative
using GPUifyLoops

export InterpolationBrick, write_interpolated_data,
       InterpolationCubedSphere, interpolate_local!
#--------------------------------------------------------
"""
    InterpolationBrick(grid::DiscontinuousSpectralElementGrid{FT}, xbnd::Array{FT,2}, xres) where FT <: AbstractFloat

This interpolation structure and the corresponding functions works for a brick, where stretching/compression happens only along the x, y & z axis.
Here x1 = X1(ξ1); x2 = X2(ξ2); x3 = X3(ξ3)

# input for the inner constructor
 - `grid` DiscontinousSpectralElementGrid
 - `xbnd` Domain boundaries in x1, x2 and x3 directions
 - `xres` Resolution of the interpolation grid in x1, x2 and x3 directions
"""
struct InterpolationBrick{FT  <:AbstractFloat,
                           T  <:Int,
                         FTV  <:AbstractVector{FT},
                         FTVD <:AbstractVector{FT},
                          TVD <:AbstractVector{T},
                        FTA2  <:Array{FT,2},
                        UI8AD <:AbstractArray{UInt8,2},
                        UI16V <: AbstractVector{UInt16},
                       UI16VD <:AbstractVector{UInt16},
                         I32V <: AbstractVector{Int32}}
    Nel::T         # # of elements
    Np::T
    Npl::T
    poly_order::T

    xbnd::FTA2     # domain bounds, [2(min/max),ndim]
    xres::FTV      # resolutions in x1, x2, x3 directions for the uniform grid

    x1g::FTV # 1D interpolation grid in x1 direction
    x2g::FTV
    x3g::FTV

    ξ1::FTVD # unique ξ1 coordinates of interpolation points within each element
    ξ2::FTVD # unique ξ2 coordinates of interpolation points within each element
    ξ3::FTVD # unique ξ3 coordinates of interpolation points within each element

    flg::UI8AD # flags when ξ1/ξ2/ξ3 interpolation point matches with a GLL point

    fac::FTVD # normalization factor

    x1i::UI16VD # interpolation grid indices of interpolation points within each element
    x2i::UI16VD
    x3i::UI16VD

    offset::TVD  # offsets for each element for v

    m1_r::FTVD  # GLL points
    m1_w::FTVD  # GLL weights
    wb::FTVD    # Barycentric weights
    # MPI setup for gathering interpolated variable on proc # 0
    Nel_all::I32V

    x1i_all::UI16V
    x2i_all::UI16V
    x3i_all::UI16V
#--------------------------------------------------------
    function InterpolationBrick(grid::DiscontinuousSpectralElementGrid{FT}, xbnd::Array{FT,2}, xres) where FT <: AbstractFloat
        mpicomm = grid.topology.mpicomm
        pid     = MPI.Comm_rank(mpicomm)
        npr     = MPI.Comm_size(mpicomm)

        DA     = arraytype(grid)                    # device array
        device = arraytype(grid) <: Array ? CPU() : CUDA()

        poly_order = polynomialorder(grid)
        qm1 = poly_order + 1
        ndim = 3
        toler = 4 * eps(FT) # tolerance

        x1g = collect(range(xbnd[1,1], xbnd[2,1], step=xres[1])); n1g = length(x1g)
        x2g = collect(range(xbnd[1,2], xbnd[2,2], step=xres[2])); n2g = length(x2g)
        x3g = collect(range(xbnd[1,3], xbnd[2,3], step=xres[3])); n3g = length(x3g)

        Np = n1g*n2g*n3g
        marker = BitArray{3}(undef,n1g,n2g,n3g); fill!(marker,true)
        #-----------------------------------------------------------------------------------
        Nel       = length(grid.topology.realelems) # # of elements on local process
        offset    = Vector{Int}(undef,Nel+1) # offsets for the interpolated variable
        n123      = zeros(Int,    ndim)        # # of unique ξ1, ξ2, ξ3 points in each cell
        xsten     = zeros(Int, 2, ndim)        # x1, x2, x3 start and end for each brick element
        xbndl     = zeros(FT,2, ndim)        # x1,x2,x3 limits (min,max) for each brick element

        ξ1 = map( i -> zeros(FT,i), zeros(Int,Nel))
        ξ2 = map( i -> zeros(FT,i), zeros(Int,Nel))
        ξ3 = map( i -> zeros(FT,i), zeros(Int,Nel))

        x1i = map( i -> zeros(UInt16,i), zeros(UInt16,Nel))
        x2i = map( i -> zeros(UInt16,i), zeros(UInt16,Nel))
        x3i = map( i -> zeros(UInt16,i), zeros(UInt16,Nel))

        x  = map( i -> zeros(FT,ndim,i), zeros(Int,Nel)) # interpolation grid points embedded in each cell

        offset[1] = 0
        #-----------------------------------------------------------------------------------
        for el in 1:Nel
            #-------------------------------------
            for (xg,dim) in zip((x1g, x2g, x3g), 1:ndim)
                xbndl[1,dim], xbndl[2,dim] = extrema(grid.topology.elemtocoord[dim,:,el])

                st = findfirst( xg .≥ xbndl[1,dim].-toler )
                if st ≠ nothing
                    if xg[st] > (xbndl[2,dim] + toler)
                        st = nothing
                    end
                end

                if st ≠ nothing
                    xsten[1,dim] = st
                    xsten[2,dim] = findlast( temp -> temp .≤ xbndl[2,dim].+toler, xg )
                	n123[dim] = xsten[2,dim] - xsten[1,dim] + 1
                else
                    n123[dim] = 0
                end
            end

            if prod(n123) > 0
                for k in xsten[1,3]:xsten[2,3], j in xsten[1,2]:xsten[2,2], i in xsten[1,1]:xsten[2,1]
                    if marker[i, j, k]
						push!(ξ1[el],  2 * ( x1g[i] - xbndl[1,1] ) / (xbndl[2,1]-xbndl[1,1]) -  1  )
						push!(ξ2[el],  2 * ( x2g[j] - xbndl[1,2] ) / (xbndl[2,2]-xbndl[1,2]) -  1  )
						push!(ξ3[el],  2 * ( x3g[k] - xbndl[1,3] ) / (xbndl[2,3]-xbndl[1,3]) -  1  )

                        push!(x1i[el], UInt16(i)); push!(x2i[el], UInt16(j)); push!(x3i[el], UInt16(k))
                        marker[i, j, k] = false
                    end
                end
                offset[el+1] = offset[el] + length(ξ1[el])
            end
    	        #-------------------------------------
        end # el loop
        #-----------------------------------------------------------------------------------
        m1_r, m1_w = GaussQuadrature.legendre(FT,qm1,GaussQuadrature.both)
        wb = Elements.baryweights(m1_r)

        Npl = offset[end]

        ξ1_d  = Array{FT}(undef,Npl);     ξ2_d  = Array{FT}(undef,Npl);     ξ3_d  = Array{FT}(undef,Npl)
        x1i_d = Array{UInt16}(undef,Npl); x2i_d = Array{UInt16}(undef,Npl); x3i_d = Array{UInt16}(undef,Npl)
        fac_d = zeros(FT,Npl);            flg_d = zeros(UInt8,3,Npl)

        for i in 1:Nel
            ctr = 1
            for j in offset[i]+1:offset[i+1]
                ξ1_d[j]  = ξ1[i][ctr]
                ξ2_d[j]  = ξ2[i][ctr]
                ξ3_d[j]  = ξ3[i][ctr]
                x1i_d[j] = x1i[i][ctr]
                x2i_d[j] = x2i[i][ctr]
                x3i_d[j] = x3i[i][ctr]
                #-----setting up interpolation
                fac1 = FT(0); fac2 = FT(0); fac3 = FT(0)
				for ib in 1:qm1

                    if abs(m1_r[ib] - ξ1_d[j]) < toler
                        @inbounds flg_d[1,j] = UInt8(ib)
                    else
                        @inbounds fac1 += wb[ib] / (ξ1_d[j]-m1_r[ib])
                    end

                    if abs(m1_r[ib] - ξ2_d[j]) < toler
                        @inbounds flg_d[2,j] = UInt8(ib)
                    else
                        @inbounds fac2 += wb[ib] / (ξ2_d[j]-m1_r[ib])
                    end

                    if abs(m1_r[ib] - ξ3_d[j]) < toler
                        @inbounds flg_d[3,j] = UInt8(ib)
                    else
                        @inbounds fac3 += wb[ib] / (ξ3_d[j]-m1_r[ib])
                    end
                end

                flg_d[1,j] ≠ UInt8(0) && (fac1 = FT(1))
                flg_d[2,j] ≠ UInt8(0) && (fac2 = FT(1))
                flg_d[3,j] ≠ UInt8(0) && (fac3 = FT(1))

                fac_d[j] = FT(1) / (fac1 * fac2 * fac3)

                #-----------------------------
                ctr += 1
            end
        end
		#----MPI setup for gathering data on proc # 0----------------------
        root = 0
        Nel_all        = zeros(Int32,npr)
        Nel_all[pid+1] = Npl

        MPI.Allreduce!(Nel_all, +, mpicomm)

        if pid ≠ root
            x1i_all = zeros(UInt16,0); x2i_all = zeros(UInt16,0); x3i_all = zeros(UInt16,0)
        else
            x1i_all = Array{UInt16}(undef,sum(Nel_all))
            x2i_all = Array{UInt16}(undef,sum(Nel_all))
            x3i_all = Array{UInt16}(undef,sum(Nel_all))
    	end

        MPI.Gatherv!(x1i_d, x1i_all, Nel_all, root, mpicomm)
        MPI.Gatherv!(x2i_d, x2i_all, Nel_all, root, mpicomm)
        MPI.Gatherv!(x3i_d, x3i_all, Nel_all, root, mpicomm)
        #-----------------------------------------------------
        if device==CUDA()
            ξ1_d  = DA(ξ1_d);  ξ2_d  = DA(ξ2_d);  ξ3_d   = DA(ξ3_d)
            x1i_d = DA(x1i_d); x2i_d = DA(x2i_d); x3i_d  = DA(x3i_d)
            flg_d = DA(flg_d); fac_d = DA(fac_d); offset = DA(offset)
            m1_r  = DA(m1_r);  m1_w  = DA(m1_w);  wb     = DA(wb)
        end
        return new{FT,
                   Int,
                   typeof(x1g),
                   typeof(ξ1_d),
                   typeof(offset),
                   typeof(xbnd),
                   typeof(flg_d),
                   typeof(x1i_all),
                   typeof(x1i_d),
                   typeof(Nel_all)}(Nel, Np, Npl, poly_order, xbnd, xres, x1g, x2g, x3g, ξ1_d, ξ2_d, ξ3_d,
                                   flg_d, fac_d, x1i_d, x2i_d, x3i_d, offset, m1_r, m1_w, wb, Nel_all, x1i_all, x2i_all, x3i_all)

    end
#--------------------------------------------------------
end # struct InterpolationBrick
#--------------------------------------------------------
"""
    interpolate_local!(intrp_brck::InterpolationBrick{FT}, sv::AbstractArray{FT}, v::AbstractArray{FT}) where {FT <: AbstractFloat}

This interpolation function works for a brick, where stretching/compression happens only along the x, y & z axis.
Here x1 = X1(ξ1); x2 = X2(ξ2); x3 = X3(ξ3)

# input
 - `intrp_brck` Initialized InterpolationBrick structure
 - `sv` State Array consisting of various variables on the discontinuous Galerkin grid
 - `v`  Interpolated variables
"""
function interpolate_local!(intrp_brck::InterpolationBrick{FT}, sv::AbstractArray{FT}, v::AbstractArray{FT}) where {FT <: AbstractFloat}
    #-------------- ----------------------------------------------------------------------------
    offset = intrp_brck.offset; m1_r = intrp_brck.m1_r; wb = intrp_brck.wb
    ξ1 = intrp_brck.ξ1; ξ2 = intrp_brck.ξ2; ξ3 = intrp_brck.ξ3
    flg = intrp_brck.flg; fac = intrp_brck.fac;
    #------------------------------------------------------------------------------------------
    qm1   = length(m1_r)
    Nel   = length(offset) - 1
    nvars = size(sv,2)

    device = typeof(sv) <: Array ? CPU() : CUDA()

    if device==CPU()
        #------------------------------------------------------------------------------------------
        Nel = length(offset) - 1

        vout    = zeros(FT,nvars)
        vout_ii = zeros(FT,nvars)
        vout_ij = zeros(FT,nvars)

        for el in 1:Nel #-----for each element elno
            np  = offset[el+1] - offset[el]
            off = offset[el]

            for i in 1:np # interpolating point-by-point
                ξ1l = ξ1[off+i]; ξ2l = ξ2[off+i]; ξ3l = ξ3[off+i]

                f1 = flg[1,off+i]; f2 = flg[2,off+i]; f3 = flg[3,off+i]
                fc = fac[off+i]

                vout .= 0.0
                f3 == 0 ? (ikloop = 1:qm1) : (ikloop = f3:f3)

                for ik in ikloop
                    #--------------------------------------------
                    vout_ij .= 0.0
                    f2 == 0 ? (ijloop = 1:qm1) : (ijloop = f2:f2)
                    for ij in ijloop #1:qm1
                        #----------------------------------------------------------------------------
                        vout_ii .= 0.0

                        if f1 == 0
                            for ii in 1:qm1
                                for vari in 1:nvars
	                                @inbounds vout_ii[vari] += sv[ii + (ij-1)*qm1 + (ik-1)*qm1*qm1, vari, el] *  wb[ii] / (ξ1l-m1_r[ii])#phir[ii]
                                end
                            end
                        else
                            for vari in 1:nvars
	                            @inbounds vout_ii[vari] = sv[f1 + (ij-1)*qm1 + (ik-1)*qm1*qm1, vari, el]
                            end
                        end
                        if f2==0
                            for vari in 1:nvars
	                            @inbounds vout_ij[vari] += vout_ii[vari] * wb[ij] / (ξ2l-m1_r[ij])#phis[ij]
                            end
                        else
                            for vari in 1:nvars
	                            @inbounds vout_ij[vari] = vout_ii[vari]
                            end
                        end
                        #----------------------------------------------------------------------------
                    end
                    if f3==0
                        for vari in 1:nvars
                            @inbounds vout[vari] += vout_ij[vari] * wb[ik] / (ξ3l-m1_r[ik])#phit[ik]
                        end
                    else
                        for vari in 1:nvars
                            @inbounds vout[vari] = vout_ij[vari]
                        end
                    end
                    #--------------------------------------------
                end
                for vari in 1:nvars
	                @inbounds v[off + i, vari] = vout[vari] * fc
				end
            end
        end
        #------------------------------------------------------------------------------------------
    else
        @cuda threads=(qm1,qm1) blocks=(Nel,nvars) shmem=qm1*(qm1+2)*sizeof(FT) interpolate_brick_CUDA!(offset, m1_r, wb, ξ1, ξ2, ξ3, flg, fac, sv, v)
    end
end
#--------------------------------------------------------

#--------------------------------------------------------
function interpolate_brick_CUDA!(offset::AbstractArray{T,1},    m1_r::AbstractArray{FT,1},   wb::AbstractArray{FT,1},
                                     ξ1::AbstractArray{FT,1},     ξ2::AbstractArray{FT,1},   ξ3::AbstractArray{FT,1},
                                    flg::AbstractArray{UInt8,2}, fac::AbstractArray{FT,1},
                                     sv::AbstractArray{FT},        v::AbstractArray{FT}) where {T <: Int, FT <: AbstractFloat}

    tj     = threadIdx().x; tk = threadIdx().y; # thread ids
    el     = blockIdx().x                       # assigning one element per block
    st_idx = blockIdx().y
    qm1    = length(m1_r)
    #--------creating views for shared memory
    shm_FT = @cuDynamicSharedMem(FT, (qm1,qm1+2))

    vout_jk = view(shm_FT,:,1:qm1)
    wb_sh   = view(shm_FT,:,qm1+1)
    m1_r_sh = view(shm_FT,:,qm1+2)
    #------loading shared memory-----------------------------
    if tk==1
        wb_sh[tj]   = wb[tj]
        m1_r_sh[tj] = m1_r[tj]
    end
    sync_threads()
    #-----------------------------------
    np  = offset[el+1] - offset[el]
    off = offset[el]

    for i in 1:np # interpolating point-by-point
        #-----------------------------------
        ξ1l = ξ1[off+i]; ξ2l = ξ2[off+i]; ξ3l = ξ3[off+i]

        f1 = flg[1,off+i]; f2 = flg[2,off+i]; f3 = flg[3,off+i]
        fc = fac[off+i]

        if f1==0 # applying phir
            @inbounds vout_jk[tj,tk] = sv[1 + (tj-1)*qm1 + (tk-1)*qm1*qm1, st_idx, el] *  wb_sh[1] / (ξ1l-m1_r_sh[1])
            for ii in 2:qm1
                @inbounds vout_jk[tj,tk] += sv[ii + (tj-1)*qm1 + (tk-1)*qm1*qm1, st_idx, el] *  wb_sh[ii] / (ξ1l-m1_r_sh[ii])
            end
        else
            @inbounds vout_jk[tj,tk] = sv[f1 + (tj-1)*qm1 + (tk-1)*qm1*qm1, st_idx, el]
        end

        if f2==0 # applying phis
            @inbounds vout_jk[tj,tk] *= (wb_sh[tj] / (ξ2l-m1_r_sh[tj]))
        end
        sync_threads()

        if tj==1 # reduction
            if f2==0
                for ij in 2:qm1
                    @inbounds vout_jk[1,tk]  += vout_jk[ij,tk]
                end
            else
                if f2 ≠ 1
                    @inbounds vout_jk[1,tk] = vout_jk[f2,tk]
                end
            end

            if f3==0 # applying phit
                @inbounds vout_jk[1,tk] *= (wb_sh[tk] / (ξ3l-m1_r_sh[tk]))
            end
        end
        sync_threads()

        if tj==1 && tk==1 # reduction
            if f3==0
                for ik in 2:qm1
                    @inbounds vout_jk[1,1] += vout_jk[1,ik]
                end
            else
                if f3 ≠ 1
                    @inbounds vout_jk[1,1] = vout_jk[1,f3]
                end
            end
            @inbounds v[off + i, st_idx] = vout_jk[1,1] * fc
        end

            #-----------------------------------
    end
    #-----------------------------------

    return nothing
end
#--------------------------------------------------------
"""
    InterpolationCubedSphere(grid::DiscontinuousSpectralElementGrid, vert_range::AbstractArray{FT}, nhor::Int, lat_res::FT, long_res::FT, rad_res::FT) where {FT <: AbstractFloat}

This interpolation structure and the corresponding functions works for a cubed sphere topology. The data is interpolated along a lat/long/rad grid.

# input for the inner constructor
 - `grid` DiscontinousSpectralElementGrid
 - `vert_range` vertex range along the radial coordinate
 - `lat_res` Resolution of the interpolation grid along the latitude coordinate in radians
 - `long_res` Resolution of the interpolation grid along the longitude coordinate in radians
 - `rad_res` Resolution of the interpolation grid along the radial coordinate
"""
struct InterpolationCubedSphere{FT <: AbstractFloat,
                                 T <: Int,
                               FTV <: AbstractVector{FT},
                              FTVD <: AbstractVector{FT},
                               TVD <: AbstractVector{T},
                             UI8AD <: AbstractArray{UInt8,2},
                             UI16V <: AbstractVector{UInt16},
                            UI16VD <: AbstractVector{UInt16},
                              I32V <: AbstractVector{Int32}}

    Nel::T
    Np::T             # # of points on the interpolation grid
    Npl::T            # # of interpolation points on the local process
    poly_order::T

    lat_min::FT;  long_min::FT;  rad_min::FT; # domain bounds, min
    lat_max::FT;  long_max::FT;  rad_max::FT; # domain bounds, max
    lat_res::FT;  long_res::FT;  rad_res::FT; # respective resolutions for the uniform grid

    n_rad::T; n_lat::T; n_long::T;             # of lat, long & rad grid locations

    rad_grd::FTV; lat_grd::FTV; long_grd::FTV; # rad, lat & long locations of interpolation grid

    ξ1::FTVD #Device array containing ξ1 coordinates of interpolation points within each element
    ξ2::FTVD #Device array containing ξ2 coordinates of interpolation points within each element
    ξ3::FTVD #Device array containing ξ3 coordinates of interpolation points within each element

    flg::UI8AD # flags when ξ1/ξ2/ξ3 interpolation point matches with a GLL point

    fac::FTVD # normalization factor

    radi::UI16VD  # rad coordinates of interpolation points within each element
    lati::UI16VD  # lat coordinates of interpolation points within each element
    longi::UI16VD # long coordinates of interpolation points within each element

    offset::TVD  # offsets for each element for v

    m1_r::FTVD   # GLL points
    m1_w::FTVD   # GLL weights
    wb::FTVD     # Barycentric weights

    # MPI setup for gathering interpolated variable on proc # 0
    Nel_all::I32V

    radi_all::UI16V
    lati_all::UI16V
    longi_all::UI16V

  #--------------------------------------------------------
    function InterpolationCubedSphere(grid::DiscontinuousSpectralElementGrid, vert_range::AbstractArray{FT}, nhor::Int, lat_res::FT, long_res::FT, rad_res::FT) where {FT <: AbstractFloat}
        mpicomm = MPI.COMM_WORLD
        pid     = MPI.Comm_rank(mpicomm)
        npr     = MPI.Comm_size(mpicomm)

        DA = arraytype(grid)                    # device array
        device = arraytype(grid) <: Array ? CPU() : CUDA()
        poly_order = polynomialorder(grid)
        qm1    = poly_order + 1
        toler1 = FT(eps(FT) * vert_range[1] * 2.0) # tolerance for unwarp function
        toler2 = FT(eps(FT) * 4.0)                 # tolerance
        toler3 = FT(eps(FT) * vert_range[1] * 10.0) # tolerance for Newton-Raphson

        lat_min,   lat_max = FT(-90.0), FT(90.0)            # inclination/zeinth angle range
        long_min, long_max = FT(-180.0), FT(180.0) 		    # azimuthal angle range
        rad_min,   rad_max = vert_range[1], vert_range[end] # radius range

        Nel = length(grid.topology.realelems) # # of local elements on the local process

        nvert_range = length(vert_range)
        nvert = nvert_range - 1              # # of elements in vertical direction
        Nel_glob = nvert * nhor * nhor * 6

        nblck = nhor * nhor * nvert
        Δh = 2 / nhor                               # horizontal grid spacing in unwarped grid

        lat_grd  = range(lat_min,  lat_max,  step=lat_res)
        long_grd = range(long_min, long_max, step=long_res)
        rad_grd  = range(rad_min,  rad_max,  step=rad_res)

        n_lat, n_long, n_rad = Int(length(lat_grd)), Int(length(long_grd)), Int(length(rad_grd))

        Np = n_lat * n_long * n_rad

        uw_grd = zeros(FT, 3)
        diffv  = zeros(FT, 3)
        ξ      = zeros(FT, 3)
        #----------------------------------------------
        glob_ord = grid.topology.origsendorder # to account for reordering of elements after the partitioning process

        glob_elem_no = zeros(Int, nvert*length(glob_ord))

        for i in 1:length(glob_ord), j in 1:nvert
            glob_elem_no[j + (i-1)*nvert] = (glob_ord[i] - 1)*nvert + j
        end
        glob_to_loc = Dict( glob_elem_no[i] => Int(i) for i in 1:Nel ) # using dictionary for speedup

        ξ1, ξ2, ξ3        = map( i -> zeros(FT,i), zeros(Int,Nel)), map( i -> zeros(FT,i), zeros(Int,Nel)), map( i -> zeros(FT,i), zeros(Int,Nel))

        radi, lati, longi = map( i -> zeros(UInt16,i), zeros(UInt16,Nel)), map( i -> zeros(UInt16,i), zeros(UInt16,Nel)), map( i -> zeros(UInt16,i), zeros(UInt16,Nel))

        offset_d = zeros(Int,Nel+1)

        #--------------------------------
        for i in 1:n_rad
            #--------------------------------
            rad = rad_grd[i]
            if rad ≤ vert_range[1]       # accounting for minor rounding errors from unwarp function at boundaries
                vert_range[1] - rad < toler1 ? l_nrm = 1 :  error("fatal error, rad lower than inner radius: ", vert_range[1] - rad," $x1_grd /// $x2_grd //// $x3_grd" )
            elseif rad ≥ vert_range[end] # accounting for minor rounding errors from unwarp function at boundaries
                rad - vert_range[end] < toler1 ? l_nrm = nvert : error("fatal error, rad greater than outer radius")
            else                         # normal scenario
                for l in 2:nvert_range
                    if vert_range[l] - rad > FT(0)
                        l_nrm = l - 1
                        break
                    end
                end
            end
            #--------------------------------
            for j in 1:n_lat
                @inbounds x3_grd = rad * sind(lat_grd[j])
                for k in 1:n_long
                    @inbounds x1_grd = rad * cosd(lat_grd[j]) * cosd(long_grd[k]) # inclination -> latitude; azimuthal -> longitude.
                    @inbounds x2_grd = rad * cosd(lat_grd[j]) * sind(long_grd[k]) # inclination -> latitude; azimuthal -> longitude.

                    uw_grd[1], uw_grd[2], uw_grd[3] = Topologies.cubedshellunwarp(x1_grd, x2_grd, x3_grd) # unwarping from sphere to cubed shell
        		    #--------------------------------
                    x1_uw2_grd = uw_grd[1] / rad # unwrapping cubed shell on to a 2D grid (in 3D space, -1 to 1 cube)
                    x2_uw2_grd = uw_grd[2] / rad
                    x3_uw2_grd = uw_grd[3] / rad
                    #--------------------------------
                    if     abs(x1_uw2_grd + 1) < toler2 # face 1 (x1 == -1 plane)
                        l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
                        l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(l_nrm + (nhor-l2)*nvert + (l3-1)*nvert*nhor)
                    elseif abs(x2_uw2_grd + 1) < toler2 # face 2 (x2 == -1 plane)
                        l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
                        l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(l_nrm + (l1-1)*nvert + (l3-1)*nvert*nhor + nblck*1)
                    elseif abs(x1_uw2_grd - 1) < toler2 # face 3 (x1 == +1 plane)
                        l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
                        l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(l_nrm + (l2-1)*nvert + (l3-1)*nvert*nhor + nblck*2 )
                    elseif abs(x3_uw2_grd - 1) < toler2 # face 4 (x3 == +1 plane)
                        l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
                        l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(l_nrm + (l1-1)*nvert + (l2-1)*nvert*nhor + nblck*3)
                    elseif abs(x2_uw2_grd - 1) < toler2 # face 5 (x2 == +1 plane)
                        l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
                        l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(l_nrm + (l1-1)*nvert + (nhor-l3)*nvert*nhor + nblck*4 )
                    elseif abs(x3_uw2_grd + 1) < toler2 # face 6 (x3 == -1 plane)
                        l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
                        l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(l_nrm + (l1-1)*nvert + (nhor-l2)*nvert*nhor + nblck*5)
                    else
                        error("error: unwrapped grid does not lie on any of the 6 faces")
                    end
    	        	#--------------------------------
                    el_loc = get(glob_to_loc,el_glob,nothing)
                    if el_loc ≠ nothing # computing inner coordinates for local elements
                        invert_trilear_mapping_hex!(view(grid.topology.elemtocoord,1,:,el_loc),
                                                    view(grid.topology.elemtocoord,2,:,el_loc),
                                                    view(grid.topology.elemtocoord,3,:,el_loc), uw_grd, diffv, toler3, ξ)
                        push!(ξ1[el_loc],ξ[1])
                        push!(ξ2[el_loc],ξ[2])
                        push!(ξ3[el_loc],ξ[3])
                        push!(radi[el_loc],  UInt16(i))
                        push!(lati[el_loc],  UInt16(j))
                        push!(longi[el_loc], UInt16(k))
                        offset_d[el_loc+1] += 1
                    end
                #--------------------------------
                end
            end
        end

        for i in 2:Nel+1
            @inbounds offset_d[i] += offset_d[i-1]
        end

        Npl = offset_d[Nel+1]

        v = Vector{FT}(undef,offset_d[Nel+1]) # Allocating storage for interpolation variable

        ξ1_d = Vector{FT}(undef,Npl); ξ2_d = Vector{FT}(undef,Npl); ξ3_d = Vector{FT}(undef,Npl)

        flg_d = zeros(UInt8,3,Npl); fac_d  = ones(FT,Npl)

        rad_d  = Vector{UInt16}(undef,Npl); lat_d  = Vector{UInt16}(undef,Npl); long_d = Vector{UInt16}(undef,Npl)

        m1_r, m1_w = GaussQuadrature.legendre(FT,qm1,GaussQuadrature.both)
        wb = Elements.baryweights(m1_r)
        for i in 1:Nel
            ctr = 1
            for j in offset_d[i]+1:offset_d[i+1]
                @inbounds ξ1_d[j]   = ξ1[i][ctr]
                @inbounds ξ2_d[j]   = ξ2[i][ctr]
                @inbounds ξ3_d[j]   = ξ3[i][ctr]
                @inbounds rad_d[j]  = radi[i][ctr]
                @inbounds lat_d[j]  = lati[i][ctr]
                @inbounds long_d[j] = longi[i][ctr]
                #-----setting up interpolation
                fac1 = FT(0); fac2 = FT(0); fac3 = FT(0)
                for ib in 1:qm1
                    if abs(m1_r[ib] - ξ1_d[j]) < toler2
                        @inbounds flg_d[1,j] = UInt8(ib)
                    else
                        @inbounds fac1 += wb[ib] / (ξ1_d[j]-m1_r[ib])
                    end

                    if abs(m1_r[ib] - ξ2_d[j]) < toler2
                        @inbounds flg_d[2,j] = UInt8(ib)
                    else
                        @inbounds fac2 += wb[ib] / (ξ2_d[j]-m1_r[ib])
                    end

                    if abs(m1_r[ib] - ξ3_d[j]) < toler2
                        @inbounds flg_d[3,j] = UInt8(ib)
                    else
                        @inbounds fac3 += wb[ib] / (ξ3_d[j]-m1_r[ib])
                    end
                end

                flg_d[1,j] ≠ 0 && (fac1 = FT(1))
                flg_d[2,j] ≠ 0 && (fac2 = FT(1))
                flg_d[3,j] ≠ 0 && (fac3 = FT(1))

                fac_d[j] = FT(1) / (fac1 * fac2 * fac3)
                #-----------------------------
                ctr += 1
            end
        end
        #----MPI setup for gathering data on proc # 0----------------------
        root = 0
        Nel_all        = zeros(Int32,npr)
        Nel_all[pid+1] = Int32(Npl)

        MPI.Allreduce!(Nel_all, +, mpicomm)

        if pid ≠ root
            radi_all = zeros(UInt16,0); lati_all = zeros(UInt16,0); longi_all = zeros(UInt16,0)
        else
            radi_all  = Array{UInt16}(undef,sum(Nel_all))
            lati_all  = Array{UInt16}(undef,sum(Nel_all))
            longi_all = Array{UInt16}(undef,sum(Nel_all))
        end

        MPI.Gatherv!(rad_d,  radi_all, Nel_all, root, mpicomm)
        MPI.Gatherv!(lat_d,  lati_all, Nel_all, root, mpicomm)
        MPI.Gatherv!(long_d,longi_all, Nel_all, root, mpicomm)
		#-----------------------------------------------------
        if device==CUDA()
            ξ1_d  = DA(ξ1_d);  ξ2_d  = DA(ξ2_d); ξ3_d = DA(ξ3_d)

            flg_d = DA(flg_d); fac_d = DA(fac_d)

            rad_d = DA(rad_d); lat_d = DA(lat_d); long_d = DA(long_d)

            m1_r = DA(m1_r); m1_w = DA(m1_w); wb = DA(wb);

            offset_d = DA(offset_d);
        end

        return new{FT,
                   Int,
                   typeof(rad_grd),
                   typeof(ξ1_d),
                   typeof(offset_d),
                   typeof(flg_d),
                   typeof(radi_all),
                   typeof(rad_d),
                   typeof(Nel_all)}(Nel, Np, Npl, poly_order, lat_min, long_min, rad_min, lat_max, long_max, rad_max, lat_res, long_res, rad_res,
                    n_rad, n_lat, n_long, rad_grd, lat_grd, long_grd, ξ1_d, ξ2_d, ξ3_d, flg_d, fac_d, rad_d, lat_d, long_d, offset_d, m1_r, m1_w, wb,
                    Nel_all, radi_all, lati_all, longi_all)
    #-----------------------------------------------------------------------------------
    end # Inner constructor function InterpolationCubedSphere
#-----------------------------------------------------------------------------------
end # structure InterpolationCubedSphere
#--------------------------------------------------------
"""
    invert_trilear_mapping_hex!(X1::AbstractArray{FT,1}, X2::AbstractArray{FT,1}, X3::AbstractArray{FT,1},
                                      x::AbstractArray{FT,1}, d::AbstractArray{FT,1}, tol::FT, ξ::AbstractArray{FT,1}) where FT <: AbstractFloat

This function computes ξ = (ξ1,ξ2,ξ3) given x = (x1,x2,x3) and the (8) vertex coordinates of a Hexahedron. Newton-Raphson method is used
# input
 - `X1` X1 coordinates of the (8) vertices of the hexahedron
 - `X2` X2 coordinates of the (8) vertices of the hexahedron
 - `X3` X3 coordinates of the (8) vertices of the hexahedron
 - `x` (x1,x2,x3) coordinates of the point
 - `d` (x1,x2,x3) coordinates, temporary storage
# output
 - `ξ` (ξ1,ξ2,ξ3) coordinates of the point
"""
function invert_trilear_mapping_hex!(X1::AbstractArray{FT,1}, X2::AbstractArray{FT,1}, X3::AbstractArray{FT,1},
                                      x::AbstractArray{FT,1}, d::AbstractArray{FT,1}, tol::FT, ξ::AbstractArray{FT,1}) where FT <: AbstractFloat
    max_it = 10     # maximum # of iterations
    ξ     .= FT(0) #zeros(FT,3,1) # initial guess => cell centroid
    trilinear_map_minus_x!(ξ, X1, X2, X3,x,d)
    err = sqrt(d[1]*d[1] + d[2]*d[2] + d[3]*d[3])
    ctr = 0
    #---Newton-Raphson iterations---------------------------
    while err > tol
        trilinear_map_IJac_x_vec!(ξ, X1, X2, X3, d)
        ξ .-= d
        trilinear_map_minus_x!(ξ, X1, X2, X3,x,d)
        err = sqrt(d[1]*d[1] + d[2]*d[2] + d[3]*d[3]) #norm(d)
        ctr += 1
        if ctr > max_it
            error("invert_trilinear_mapping_hex: Newton-Raphson not converging to desired tolerance after max_it = ", max_it," iterations; err = ", err,"; toler = ", tol)
        end
    end
    #-------------------------------------------------------
    clamp!(ξ,FT(-1),FT(1))
    return nothing
end
#--------------------------------------------------------
function trilinear_map_minus_x!(ξ::AbstractArray{FT,1}, x1v::AbstractArray{FT,1}, x2v::AbstractArray{FT,1}, x3v::AbstractArray{FT,1},
                                x::AbstractArray{FT,1}, d::AbstractArray{FT,1}) where FT <: AbstractFloat
    p1 = 1 + ξ[1]; p2 = 1 + ξ[2]; p3 = 1 + ξ[3];
    m1 = 1 - ξ[1]; m2 = 1 - ξ[2]; m3 = 1 - ξ[3];

    d[1] = ( m1 * ( m2 * ( m3 * x1v[1] + p3 * x1v[5] ) +  p2 * ( m3 * x1v[3] + p3 * x1v[7] ) ) +
             p1 * ( m2 * ( m3 * x1v[2] + p3 * x1v[6] ) +  p2 * ( m3 * x1v[4] + p3 * x1v[8] ) ) )/8.0 - x[1]

    d[2] = ( m1 * ( m2 * ( m3 * x2v[1] + p3 * x2v[5] ) +  p2 * ( m3 * x2v[3] + p3 * x2v[7] ) ) +
             p1 * ( m2 * ( m3 * x2v[2] + p3 * x2v[6] ) +  p2 * ( m3 * x2v[4] + p3 * x2v[8] ) ) )/8.0 - x[2]

    d[3] = ( m1 * ( m2 * ( m3 * x3v[1] + p3 * x3v[5] ) +  p2 * ( m3 * x3v[3] + p3 * x3v[7] ) ) +
             p1 * ( m2 * ( m3 * x3v[2] + p3 * x3v[6] ) +  p2 * ( m3 * x3v[4] + p3 * x3v[8] ) ) )/8.0 - x[3]

    return nothing
end
#--------------------------------------------------------
function trilinear_map_IJac_x_vec!(ξ::AbstractArray{FT,1}, x1v::AbstractArray{FT,1}, x2v::AbstractArray{FT,1}, x3v::AbstractArray{FT,1},
                                   v::AbstractArray{FT,1}) where FT <: AbstractFloat
    p1 = 1 + ξ[1]; p2 = 1 + ξ[2]; p3 = 1 + ξ[3];
    m1 = 1 - ξ[1]; m2 = 1 - ξ[2]; m3 = 1 - ξ[3];
    #-------------------------------------------------------------------------
    Jac11 = ( m2 * ( m3 * (x1v[2] - x1v[1]) + p3 * (x1v[6] - x1v[5]) ) +
              p2 * ( m3 * (x1v[4] - x1v[3]) + p3 * (x1v[8] - x1v[7]) ) )/ 8.0

    Jac12 = ( m1 * ( m3 * (x1v[3] - x1v[1]) + p3 * (x1v[7] - x1v[5]) ) +
              p1 * ( m3 * (x1v[4] - x1v[2]) + p3 * (x1v[8] - x1v[6]) ) )/ 8.0

    Jac13 = ( m1 * ( m2 * (x1v[5] - x1v[1]) + p2 * (x1v[7] - x1v[3]) ) +
              p1 * ( m2 * (x1v[6] - x1v[2]) + p2 * (x1v[8] - x1v[4]) ) )/ 8.0
    #-------------------------------------------------------------------------
    Jac21 = ( m2 * ( m3 * (x2v[2] - x2v[1]) + p3 * (x2v[6] - x2v[5]) ) +
              p2 * ( m3 * (x2v[4] - x2v[3]) + p3 * (x2v[8] - x2v[7]) ) )/ 8.0

    Jac22 = ( m1 * ( m3 * (x2v[3] - x2v[1]) + p3 * (x2v[7] - x2v[5]) ) +
              p1 * ( m3 * (x2v[4] - x2v[2]) + p3 * (x2v[8] - x2v[6]) ) )/ 8.0

    Jac23 = ( m1 * ( m2 * (x2v[5] - x2v[1]) + p2 * (x2v[7] - x2v[3]) ) +
              p1 * ( m2 * (x2v[6] - x2v[2]) + p2 * (x2v[8] - x2v[4]) ) )/ 8.0
    #-------------------------------------------------------------------------
    Jac31 = ( m2 * ( m3 * (x3v[2] - x3v[1]) + p3 * (x3v[6] - x3v[5]) ) +
              p2 * ( m3 * (x3v[4] - x3v[3]) + p3 * (x3v[8] - x3v[7]) ) )/ 8.0

    Jac32 = ( m1 * ( m3 * (x3v[3] - x3v[1]) + p3 * (x3v[7] - x3v[5]) ) +
              p1 * ( m3 * (x3v[4] - x3v[2]) + p3 * (x3v[8] - x3v[6]) ) )/ 8.0

    Jac33 = ( m1 * ( m2 * (x3v[5] - x3v[1]) + p2 * (x3v[7] - x3v[3]) ) +
              p1 * ( m2 * (x3v[6] - x3v[2]) + p2 * (x3v[8] - x3v[4]) ) )/ 8.0
    #-------------------------------------------------------------------------
    # computing cofactor matrix
    C11 =  Jac22*Jac33 - Jac23*Jac32; C12 = -Jac21*Jac33 + Jac23*Jac31; C13 =  Jac21*Jac32 - Jac22*Jac31
    C21 = -Jac12*Jac33 + Jac13*Jac32; C22 =  Jac11*Jac33 - Jac13*Jac31; C23 = -Jac11*Jac32 + Jac12*Jac31
    C31 =  Jac12*Jac23 - Jac13*Jac22; C32 = -Jac11*Jac23 + Jac13*Jac21; C33 =  Jac11*Jac22 - Jac12*Jac21;

    # computing determinant
    det = Jac11*C11 + Jac12*C12 + Jac13*C13

    Jac11 = (C11*v[1] + C21*v[2] + C31*v[3]) / det
    Jac21 = (C12*v[1] + C22*v[2] + C32*v[3]) / det
    Jac31 = (C13*v[1] + C23*v[2] + C33*v[3]) / det

    v[1] = Jac11; v[2] = Jac21; v[3] = Jac31

    return nothing
end
#--------------------------------------------------------
"""
    interpolate_local!(intrp_cs::InterpolationCubedSphere{FT}, sv::AbstractArray{FT}, v::AbstractArray{FT}) where {FT <: AbstractFloat}

This interpolation function works for cubed spherical shell geometry.

# input
 - `intrp_cs` Initialized cubed sphere structure
 - `sv` State Array consisting of various variables on the discontinuous Galerkin grid
 - `v`  Interpolated variables
"""
function interpolate_local!(intrp_cs::InterpolationCubedSphere{FT}, sv::AbstractArray{FT}, v::AbstractArray{FT}) where {FT <: AbstractFloat}
    #------------------------------------------------------------------------------------------
    offset = intrp_cs.offset; m1_r = intrp_cs.m1_r; wb = intrp_cs.wb;
    ξ1 = intrp_cs.ξ1; ξ2 = intrp_cs.ξ2; ξ3 = intrp_cs.ξ3
    flg = intrp_cs.flg; fac = intrp_cs.fac
    #------------------------------------------------------------------------------------------
    qm1   = length(m1_r)
    nvars = size(sv,2)
    Nel   = length(offset) - 1

    device = typeof(sv) <: Array ? CPU() : CUDA()

    if device==CPU()
        #------------------------------------------------------------------------------------------
        Nel   = length(offset) - 1

        vout    = zeros(FT,nvars) #FT(0)
        vout_ii = zeros(FT,nvars) #FT(0)
        vout_ij = zeros(FT,nvars) #FT(0)

        for el in 1:Nel #-----for each element elno
            np  = offset[el+1] - offset[el]
            off = offset[el]

            for i in 1:np # interpolating point-by-point
                ξ1l = ξ1[off+i]; ξ2l = ξ2[off+i]; ξ3l = ξ3[off+i]

                f1 = flg[1,off+i]; f2 = flg[2,off+i]; f3 = flg[3,off+i]
                fc = fac[off+i]

                vout .= 0.0
                f3 == 0 ? (ikloop = 1:qm1) : (ikloop = f3:f3)

                for ik in ikloop
                    #--------------------------------------------
                    vout_ij .= 0.0
                    f2 == 0 ? (ijloop = 1:qm1) : (ijloop = f2:f2)
                    for ij in ijloop #1:qm1
                        #----------------------------------------------------------------------------
                        vout_ii .= 0.0

                        if f1 == 0
                            for ii in 1:qm1
                                for vari in 1:nvars
                                    @inbounds vout_ii[vari] += sv[ii + (ij-1)*qm1 + (ik-1)*qm1*qm1, vari, el] *  wb[ii] / (ξ1l-m1_r[ii])#phir[ii]
                                end
                            end
                        else
                            for vari in 1:nvars
                                @inbounds vout_ii[vari] = sv[f1 + (ij-1)*qm1 + (ik-1)*qm1*qm1, vari, el]
                            end
                        end
                        if f2==0
                            for vari in 1:nvars
                            	@inbounds vout_ij[vari] += vout_ii[vari] * wb[ij] / (ξ2l-m1_r[ij])#phis[ij]
                            end
                        else
                            for vari in 1:nvars
                                @inbounds vout_ij[vari] = vout_ii[vari]
                            end
                        end
                        #----------------------------------------------------------------------------
                    end
                    if f3==0
                        for vari in 1:nvars
                            @inbounds vout[vari] += vout_ij[vari] * wb[ik] / (ξ3l-m1_r[ik])#phit[ik]
                        end
                    else
                        for vari in 1:nvars
                            @inbounds vout[vari] = vout_ij[vari]
                        end
                    end
                    #--------------------------------------------
                end
                for vari in 1:nvars
                    @inbounds v[off + i,vari] = vout[vari] * fc
                end
            end
        end
        #------------------------------------------------------------------------------------------
    else
        #------------------------------------------------------------------------------------------
        @cuda threads=(qm1,qm1) blocks=(Nel,nvars) shmem=qm1*(qm1+2)*sizeof(FT) interpolate_cubed_sphere_CUDA!(offset, m1_r, wb, ξ1, ξ2, ξ3, flg, fac, sv, v)
        #------------------------------------------------------------------------------------------
    end
    return nothing
end
#--------------------------------------------------------
function interpolate_cubed_sphere_CUDA!(offset::AbstractArray{T,1}, m1_r::AbstractArray{FT,1}, wb::AbstractArray{FT,1},
                                        ξ1::AbstractArray{FT,1}, ξ2::AbstractArray{FT,1}, ξ3::AbstractArray{FT,1},
                                        flg::AbstractArray{UInt8,2}, fac::AbstractArray{FT,1},
                                        sv::AbstractArray{FT}, v::AbstractArray{FT}) where {T <: Integer, FT <: AbstractFloat}

    tj    = threadIdx().x; tk = threadIdx().y; # thread ids
    el    = blockIdx().x                       # assigning one element per block
    st_no = blockIdx().y

    qm1 = length(m1_r)
    nvars = size(sv,2)
    #--------creating views for shared memory
    shm_FT = @cuDynamicSharedMem(FT, (qm1,qm1+2))

    vout_jk = view(shm_FT,:,1:qm1)
    wb_sh   = view(shm_FT,:,qm1+1)
    m1_r_sh = view(shm_FT,:,qm1+2)
    #------loading shared memory-----------------------------
    if tk==1
        wb_sh[tj]   = wb[tj]
        m1_r_sh[tj] = m1_r[tj]
    end
    sync_threads()
    #-----------------------------------
    np  = offset[el+1] - offset[el]
    off = offset[el]

    for i in 1:np # interpolating point-by-point
        #-----------------------------------
        ξ1l = ξ1[off+i]; ξ2l = ξ2[off+i]; ξ3l = ξ3[off+i]

        f1 = flg[1,off+i]; f2 = flg[2,off+i]; f3 = flg[3,off+i]
        fc = fac[off+i]

        if f1==0 # applying phir
            @inbounds vout_jk[tj,tk] = sv[1 + (tj-1)*qm1 + (tk-1)*qm1*qm1, st_no, el] *  wb_sh[1] / (ξ1l-m1_r_sh[1])
            for ii in 2:qm1
                @inbounds vout_jk[tj,tk] += sv[ii + (tj-1)*qm1 + (tk-1)*qm1*qm1, st_no, el] *  wb_sh[ii] / (ξ1l-m1_r_sh[ii])
            end
        else
            @inbounds vout_jk[tj,tk] = sv[f1 + (tj-1)*qm1 + (tk-1)*qm1*qm1, st_no, el]
        end

        if f2==0 # applying phis
            @inbounds vout_jk[tj,tk] *= (wb_sh[tj] / (ξ2l-m1_r_sh[tj]))
        end
        sync_threads()

        if tj==1 # reduction
            if f2==0
                for ij in 2:qm1
                    @inbounds vout_jk[1,tk]  += vout_jk[ij,tk]
                end
            else
                if f2 ≠ 1
                    @inbounds vout_jk[1,tk] = vout_jk[f2,tk]
                end
            end

            if f3==0 # applying phit
                @inbounds vout_jk[1,tk] *= (wb_sh[tk] / (ξ3l-m1_r_sh[tk]))
            end
        end
        sync_threads()

        if tj==1 && tk==1 # reduction
            if f3==0
                for ik in 2:qm1
                    @inbounds vout_jk[1,1] += vout_jk[1,ik]
                end
            else
                if f3 ≠ 1
                    @inbounds vout_jk[1,1] = vout_jk[1,f3]
                end
            end
            @inbounds v[off + i, st_no] = vout_jk[1,1] * fc
        end
        #-----------------------------------
    end
    #-----------------------------------
    return nothing
end
#--------------------------------------------------------
"""
    write_interpolated_data(intrp_cs::InterpolationCubedSphere{FT}, iv::AbstractArray{FT}, varnames, filename) where {FT <: AbstractFloat}

This interpolation function gathers interpolated data onto process # 0 and writed data to a NetCDF file.

# input
 - `intrp_cs` Initialized cubed sphere interpolation structure
 - `iv` interpolated variables
 - `varnames` Tuple of Interpolated variable name strings
 - `filename` Filename of the NetCDF file to be written
"""
function write_interpolated_data(intrp_cs::InterpolationCubedSphere{FT}, iv::AbstractArray{FT}, varnames, filename) where {FT <: AbstractFloat}

    device  = CLIMA.array_type() <: Array ? CPU() : CUDA()
    DA      = CLIMA.array_type()                    # device array

    mpicomm = MPI.COMM_WORLD
    pid     = MPI.Comm_rank(mpicomm)
    npr     = MPI.Comm_size(mpicomm)
    root    = 0

    nrad    = length(intrp_cs.rad_grd)
    nlat    = length(intrp_cs.lat_grd)
    nlong   = length(intrp_cs.long_grd)
    nvars   = size(iv,2)

    if length(varnames) ≠ nvars
        error("# of varnames $(length(varnames)) do not match with number of variables $nvars ")
    end
    # Gathering data onto process id 0
    v = deepcopy(Array(iv))

    pid==0 ? v_all = Array{FT}(undef,length(intrp_cs.radi_all),nvars) : v_all = Array{FT}(undef,0,nvars)

    for vari in 1:nvars
        MPI.Gatherv!(view(v,:,vari), view(v_all,:,vari), intrp_cs.Nel_all, root, mpicomm)
    end

    if pid==0
        svi   = Array{FT}(undef,nrad,nlat,nlong,nvars)
        radi  = intrp_cs.radi_all
        lati  = intrp_cs.lati_all
        longi = intrp_cs.longi_all

        for i in 1:length(radi)
            for vari in 1:nvars
                svi[radi[i], lati[i], longi[i], vari] = v_all[i,vari]
            end
        end
        #-------writing data to NetCDF file
        ds = Dataset(filename,"c") # creating a netCDF file

        ds.dim["rad"]  = nrad    # defining dimensions
        ds.dim["lat"]  = nlat
        ds.dim["long"] = nlong

        rad_grd  = defVar(ds, "rad",  intrp_cs.rad_grd,  ("rad",))
        lat_grd  = defVar(ds, "lat",  intrp_cs.lat_grd,  ("lat",))
        long_grd = defVar(ds, "long", intrp_cs.long_grd, ("long",))

        for vari in 1:nvars
            sym = Symbol(varnames[vari])
            @eval ($sym = defVar($ds,$(varnames[vari]), $(svi[:,:,:,vari]), ("rad", "lat","long")))
        end

        close(ds)
        #--------------------------------------------------------------------------------------------------
    else
        svi = Array{FT}(undef,0,0,0,nvars)
	end
    return svi
end
#--------------------------------------------------------
"""
    write_interpolated_data(intrp_brck::InterpolationBrick{FT}, iv::AbstractArray{FT}, varnames, filename) where {FT <:AbstractFloat}

This interpolation function gathers interpolated data onto process # 0 and writed data to a NetCDF file.

# input
 - `intrp_brck` Initialized brick interpolation structure
 - `iv` interpolated variables
 - `varnames` Tuple of Interpolated variable name strings
 - `filename` Filename of the NetCDF file to be written
# output
 - `svi` full interpolated variables on process # 0
"""
function write_interpolated_data(intrp_brck::InterpolationBrick{FT}, iv::AbstractArray{FT}, varnames, filename) where {FT <: AbstractFloat}

    device  = CLIMA.array_type() <: Array ? CPU() : CUDA()
    DA      = CLIMA.array_type()                    # device array

    mpicomm = MPI.COMM_WORLD
    pid     = MPI.Comm_rank(mpicomm)
    npr     = MPI.Comm_size(mpicomm)
    root    = 0

    nx1    = length(intrp_brck.x1g)
    nx2    = length(intrp_brck.x2g)
    nx3    = length(intrp_brck.x3g)
    nvars  = size(iv,2)

    if length(varnames) ≠ nvars
        error("# of varnames $(length(varnames)) do not match with number of variables $nvars ")
    end

    v = deepcopy(Array(iv))

    if pid==0
        v_all = Array{FT}(undef,length(intrp_brck.x1i_all),nvars)
    else
        v_all = Array{FT}(undef,0,nvars)
    end
    #----------------------
    for vari in 1:nvars
        MPI.Gatherv!(view(v,:,vari), view(v_all,:,vari), intrp_brck.Nel_all, root, mpicomm)
    end
    #----------------------
    if pid==0
        #------------------------
        svi = Array{FT}(undef,nx1,nx2,nx3,nvars)
        x1i = intrp_brck.x1i_all
        x2i = intrp_brck.x2i_all
        x3i = intrp_brck.x3i_all

        for vari in 1:nvars
            for i in 1:length(x1i)
                svi[x1i[i], x2i[i], x3i[i], vari] = v_all[i,vari]
            end
        end
        #-------writing data to NetCDF file
        ds = Dataset(filename,"c") # creating a netCDF file

        ds.dim["x1"] = nx1    # defining dimensions
        ds.dim["x2"] = nx2
        ds.dim["x3"] = nx3

        x1_grd = defVar(ds, "x1", intrp_brck.x1g, ("x1",))
        x2_grd = defVar(ds, "x2", intrp_brck.x2g, ("x2",))
        x3_grd = defVar(ds, "x3", intrp_brck.x3g, ("x3",))

        for vari in 1:nvars
            sym = Symbol(varnames[vari])
            @eval ($sym = defVar($ds,$(varnames[vari]), $(svi[:,:,:,vari]), ("x1","x2","x3")))
        end

        close(ds)
        #--------------------------------------------------------------------------------------------------
    else
        svi = Array{FT}(undef,0,0,0,0)
    end
    return svi
end
#--------------------------------------------------------
end # module interploation

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
using DecisionMakingUtils
using Test

@testset "ZeroOneNormalization Test" begin
    k = 3
    T = Float32
    low = -collect(T, 1:k)
    high = collect(T, 1:k)
    ranges = hcat(low, high)
    n = ZeroOneNormalization(ranges)
    x = zeros(T, k)
    buff = zero(x)
    y = n(x)
    @test typeof(y) == typeof(x)
    @test size(y) == size(x)

    @test all(y .≥ 0.0)
    @test all(y .≤ 1.0)
    y = n(low)
    @test all(y .≥ 0.0)
    @test all(y .≤ 1.0)
    y = n(high)
    @test all(y .≥ 0.0)
    @test all(y .≤ 1.0)
    n(buff, high)
    @test all(isapprox.(buff, y))
end

@testset "PosNegNormalization Test" begin
    k = 3
    T = Float32
    low = -collect(T, 1:k)
    high = collect(T, 1:k)
    ranges = hcat(low, high)
    n = PosNegNormalization(ranges)

    x = zeros(T, k)
    buff = zero(x)
    y = n(x)
    @test typeof(y) == typeof(x)
    @test size(y) == size(x)
    @test all(y .≥ -1.0)
    @test all(y .≤ 1.0)
    y = n(low)
    println(low)
    @test all(y .≥ -1.0)
    @test all(y .≤ 1.0)
    y = n(high)
    @test all(y .≥ -1.0)
    @test all(y .≤ 1.0)
    n(buff, high)
    @test all(isapprox.(buff, y))
end

@testset "GaussianNormalization Test" begin
    k = 3
    T = Float32
    μ = collect(T, 1:k)
    σ = collect(T, 1:k)
    n = GaussianNormalization(μ, σ)

    x = zeros(T, k)
    buff = zero(x)
    y = n(x)
    @test typeof(y) == typeof(x)
    @test size(y) == size(x)
    println(y)
    @test all(isapprox.(y, -1.0))
    y = n(μ)
    @test all(isapprox.(y, 0.0))
    n(buff, μ)
    @test all(isapprox.(buff, y))
end


text:
# Pairwise Alignment
# ==================
#
# Pairwise alignment type.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
Pairwise alignment
"""
type PairwiseAlignment{S1,S2}
    a::AlignedSequence{S1}
    b::S2
end

Base.start(aln::PairwiseAlignment) = (2, 1)
Base.done(aln::PairwiseAlignment, ij) = ij[1] > endof(aln.a.aln.anchors)
function Base.next(aln::PairwiseAlignment, ij)
    i, j = ij
    anchors = aln.a.aln.anchors
    anchor = anchors[i]
    seq = aln.a.seq
    ref = aln.b
    seqpos = anchors[i-1].seqpos
    refpos = anchors[i-1].refpos

    if ismatchop(anchor.op)
        x = seq[seqpos + j]
        y = ref[refpos + j]
    elseif isinsertop(anchor.op)
        x = seq[seqpos + j]
        y = Bio.Seq.gap(eltype(ref))
    elseif isdeleteop(anchor.op)
        x = Bio.Seq.gap(eltype(seq))
        y = ref[refpos + j]
    else
        @assert false
    end

    if ismatchop(anchor.op) || isinsertop(anchor.op)
        if j < anchor.seqpos - seqpos
            j += 1
        else
            i += 1
            j = 1
        end
    else
        if j < anchor.refpos - refpos
            j += 1
        else
            i += 1
            j = 1
        end
    end

    return (x, y), (i, j)
end

Base.length(aln::PairwiseAlignment) = count_aligned(aln)
Base.eltype{S1,S2}(::Type{PairwiseAlignment{S1,S2}}) = Tuple{eltype(S1),eltype(S2)}

"""
    count(aln::PairwiseAlignment, target::Operation)

Count the number of positions where the `target` operation is applied.
"""
function Base.count(aln::PairwiseAlignment, target::Operation)
    anchors = aln.a.aln.anchors
    n = 0
    for i in 2:endof(anchors)
        op = anchors[i].op
        if op == target
            if ismatchop(op) || isinsertop(op)
                n += anchors[i].seqpos - anchors[i-1].seqpos
            elseif isdeleteop(op)
                n += anchors[i].refpos - anchors[i-1].refpos
            end
        end
    end
    return n
end

"""
    count_matches(aln)

Count the number of matching positions.
"""
count_matches(aln::PairwiseAlignment) = count(aln, OP_SEQ_MATCH)

"""
    count_mismatches(aln)

Count the number of mismatching positions.
"""
count_mismatches(aln::PairwiseAlignment) = count(aln, OP_SEQ_MISMATCH)

"""
    count_insertions(aln)

Count the number of inserting positions.
"""
count_insertions(aln::PairwiseAlignment) = count(aln, OP_INSERT)

"""
    count_deletions(aln)

Count the number of deleting positions.
"""
count_deletions(aln::PairwiseAlignment) = count(aln, OP_DELETE)

"""
    count_aligned(aln)

Count the number of aligned positions.
"""
function count_aligned(aln::PairwiseAlignment)
    anchors = aln.a.aln.anchors
    n = 0
    for i in 2:endof(anchors)
        op = anchors[i].op
        if ismatchop(op) || isinsertop(op)
            n += anchors[i].seqpos - anchors[i-1].seqpos
        elseif isdeleteop(op)
            n += anchors[i].refpos - anchors[i-1].refpos
        end
    end
    return n
end

function ref2seq(aln::PairwiseAlignment, i::Integer)
    return ref2seq(aln.a, i)
end

function seq2ref(aln::PairwiseAlignment, i::Integer)
    return seq2ref(aln.a, i)
end

# Printers
# --------

function Base.show{S1,S2}(io::IO, aln::PairwiseAlignment{S1,S2})
    println(io, "PairwiseAlignment{", S1, ",", S2, "}:")
    print(io, aln)
end

function Base.print(io::IO, aln::PairwiseAlignment, width::Integer=60)
    seq = aln.a.seq
    ref = aln.b
    anchors = aln.a.aln.anchors
    # width of position numbers
    posw = ndigits(max(anchors[end].seqpos, anchors[end].refpos)) + 1

    i = 0
    seqpos = anchors[1].seqpos
    refpos = anchors[1].refpos
    seqbuf = IOBuffer()
    refbuf = IOBuffer()
    matbuf = IOBuffer()
    s = start(aln)
    while !done(aln, s)
        (x, y), s = next(aln, s)

        i += 1
        if x != Bio.Seq.gap(eltype(seq))
            seqpos += 1
        end
        if y != Bio.Seq.gap(eltype(ref))
            refpos += 1
        end

        if i % width == 1
            print(seqbuf, "  seq:", lpad(seqpos, posw), ' ')
            print(refbuf, "  ref:", lpad(refpos, posw), ' ')
            print(matbuf, " "^(posw + 7))
        end

        print(seqbuf, x)
        print(refbuf, y)
        print(matbuf, x == y ? '|' : ' ')

        if i % width == 0
            print(seqbuf, lpad(seqpos, posw))
            print(refbuf, lpad(refpos, posw))
            print(matbuf)

            println(io, takebuf_string(seqbuf))
            println(io, takebuf_string(matbuf))
            println(io, takebuf_string(refbuf))

            if !done(aln, s)
                println(io)
                seek(seqbuf, 0)
                seek(matbuf, 0)
                seek(refbuf, 0)
            end
        end
    end

    if i % width != 0
        print(seqbuf, lpad(seqpos, posw))
        print(refbuf, lpad(refpos, posw))
        print(matbuf)

        println(io, takebuf_string(seqbuf))
        println(io, takebuf_string(matbuf))
        println(io, takebuf_string(refbuf))
    end
end


text:
# Following RCP8.5 CO2e concentrations
# Raw data at https://www.iiasa.ac.at/web-apps/tnt/RcpDb/dsd?Action=htmlpage&page=compare
#
# See below link for 2020 initial condition:
# https://www.eea.europa.eu/data-and-maps/indicators/atmospheric-greenhouse-gas-concentrations-6/assessment-1
function ramp_emissions(t, q0::Float64, n::Float64, t1::Float64, t2::Float64)
    t0 = t[1]
    Δt0 = t1 - t0
    Δt1 = t2 - t1
    q = zeros(size(t))
    increase_idx = (t .<= t1)
    decrease_idx = ((t .> t1) .& (t .<= t2))
    q[increase_idx] .= q0 * (1. .+ (n-1) .*(t[increase_idx] .- t0)/Δt0)
    q[decrease_idx] .= n * q0 * (t2 .- t[decrease_idx])/Δt1
    q[t .> t2] .= 0.
    return q
end
function ramp_emissions(dom::Domain)
    return ramp_emissions(t(dom), 7.5, 3., 2100., 2150.)
end

emissions(q, M) = q .* (1. .- M)
function emissions(m::ClimateModel; M=false)
    return emissions(
        m.economics.baseline_emissions,
        m.controls.mitigate .* (1. .- .~past_mask(m) * ~M)
    )
end

function effective_emissions(r, q, M, R)
    return r*(emissions(q, M) .- q[1]*R)
end
function effective_emissions(m; M=false, R=false)
    return effective_emissions(
        m.physics.r,
        m.economics.baseline_emissions,
        m.controls.mitigate .* (1. .- .~past_mask(m) * ~M),
        m.controls.remove .* (1. .- .~past_mask(m) * ~R)
    )
end

function c(c0, effective_emissions, dt)
    return c0 .+ cumsum(effective_emissions * dt)
end

function c(m; M=false, R=false)
    return c(
        m.physics.c0,
        effective_emissions(m, M=M, R=R),
        m.domain.dt
    )
end


text:
@connector function MagneticPort(;name, V_m_start=0.0, Phi_start=0.0)
    @variables V_m(t)=V_m_start # [Wb] Magnetic potential at the port
    @variables Phi(t)=Phi_start [connect=Flow] # [A] Magnetic flux flowing into the port"
    ODESystem(Equation[], t, [V_m, Phi], []; name=name)
end
Base.@doc "Port for a Magnetic system." MagneticPort

"""
Positive magnetic port
"""
const PositiveMagneticPort = MagneticPort

"""
Negative magnetic port
"""
const NegativeMagneticPort = MagneticPort

"""
    TwoPort(;name, V_m_start=0.0, Phi_start=0.0)

Partial component with magnetic potential difference between two magnetic ports p and n and magnetic flux Phi from p to n.

# Parameters:
- `V_m_start`: Initial magnetic potential difference between both ports
- `Phi_start`: Initial magnetic flux from port_p to port_n
"""
function TwoPort(;name, V_m_start=0.0, Phi_start=0.0)
    @named port_p = PositiveMagneticPort()
    @named port_n = NegativeMagneticPort()
    @variables V_m(t)=V_m_start Phi(t)=Phi_start
    eqs = [
        V_m ~ port_p.V_m - port_n.V_m
        Phi ~ port_p.Phi
        0 ~ port_p.Phi + port_n.Phi
    ]
    ODESystem(eqs, t, [V_m, Phi], [], systems=[port_p, port_n], name=name)
end