text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Function to obtain the price index
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

"""
    create_price_index_hat(w_hat::Vector, τ_hat_Z::Matrix, τ_hat_F::Matrix, VA_coeff::Vector, γ::Matrix, π_Z::Matrix, π_F::Matrix, θ::Vector)

The function computes the optimal price indices (change) conditional on the economys' wage structure by minimizing the distance between successive iterations
    in the intermediate input price index.

# Arguments
- `w_hat::Vector`: N×1, country wage (change) vector w_hat.
- `τ_hat_Z::Matrix`: NS×NS, origin country-industry destination country-industry intermediate goods demand trade cost (change) matrix τ_hat_Z.
- `τ_hat_F::Matrix`: NS×N, origin country-industry destination country final goods trade cost (change) matrix τ_hat_F.
- `VA_coeff::Vector`: NS×1, country-industry value added coefficients vector VA_coeff.
- `γ::Matrix`: S×NS, country-industry intermediate input expenditure share matrix γ.
- `π_Z::Matrix`: NS×NS, origin country-industry destination country-industry intermediate goods trade share matrix π_Z.
- `π_F::Matrix`: NS×N, origin country-industry destination country final goods trade share matrix π_F.
- `θ::Vector`: NS×1, country-industry trade elasticities vector θ.

# Output
- `P_hat_Z::Matrix{Float64}`: S×NS, country-industry composite intermediate goods price index (change) per industry matrix P_hat_Z.
- `P_hat_F::Matrix{Float64}`: S×N, country-industry composite final goods price index (change) matrix P_hat_F.
- `π_hat_Z::Matrix{Float64}`: NS×NS, origin country-industry destination country-industry intermediate goods trade share (change) matrix π_hat_Z.
- `π_hat_F::Matrix{Float64}`: NS×N, origin country-industry destination country final goods trade share (change) matrix π_hat_F.
- `cost_hat::Matrix{Float64}`: NS×1, country-industry production cost (change) vector cost_hat.

# Examples
```julia-repl
julia> P_hat_Z, P_hat_F, π_hat_Z, π_hat_F, cost_hat = create_price_index_hat(w_hat, τ_hat_Z, τ_hat_F, VA_coeff, γ, π_Z, π_F, θ)
```
"""

function create_price_index_hat(w_hat::Vector, τ_hat_Z::Matrix, τ_hat_F::Matrix, VA_coeff::Vector, γ::Matrix, π_Z::Matrix, π_F::Matrix, θ::Vector)

    # initialize
    w0_hat = copy(w_hat) # N×1, update wage vector
    P0_hat_Z = ones(S, N*S) # S×NS, intermediate goods price indices

    # iteration parameters
    tolerance_p = 1e-3
    max_iteration_p = 100
    iteration_p = 0
    max_error_p = 1.0

    while max_error_p > tolerance_p && iteration_p <= max_iteration_p

        # # cost structure of the economy, equation (47)
        # cost_hat_w = [w0_hat[ceil(Int,i/S)]^VA_coeff[i] for i in 1:N*S] # NS×1, wages
        # cost_hat_p = P0_hat_Z.^γ
        # cost_hat_p = [prod(cost_hat_p[:,i]) for i in 1:N*S] # NS×1, price indices (prod is the same as sum just for multiplication)
        # cost_hat = cost_hat_w .* cost_hat_p # NS×1
        # global cost_hat = reshape(cost_hat, S, N) # S×N

        global cost_hat = ones(S,N)
        VAcoeff = reshape(VA_coeff, S, N)
        for s in 1:S
            for j in 1:N
                cost_hat[s,j] = w0_hat[j]^VAcoeff[s,j]
                for r in 1:S
                    cost_hat[s,j] = cost_hat[s,j]*P0_hat_Z[r,(j-1)*S+s]^γ[r,(j-1)*S+s]
                end
            end
        end

        # ----------
        # price index for intermediate goods from equation (48)
        global cost_hat_Z = zeros(N*S, N*S)
        θ = reshape(θ, S, N) # S×N

        for i in 1:N
            for r in 1:S
                for j in 1:N
                    for s in 1:S
                        iirr = (i-1)*S+r
                        jjss = (j-1)*S+s
                        cost_hat_Z[iirr,jjss] = (cost_hat[r,i]*τ_hat_Z[iirr,jjss])^(-θ[r,i])
                    end
                end
            end
        end

        # price index calculation
        global P_hat_Z = zeros(S, N*S)
        for r in 1:S
            for j in 1:N
                for s in 1:S
                    jjss = (j-1)*S+s
                    P_hat_Z[r,jjss] = (sum(π_Z[r:S:(N-1)*S+r,jjss] .* cost_hat_Z[r:S:(N-1)*S+r,jjss]))^(-1/θ[r,j])
                end
            end
        end

        P_hat_Z = ifelse.(isinf.(P_hat_Z), 0.0, P_hat_Z)

        # update iteration parameters
        error = abs.(P_hat_Z .- P0_hat_Z)
        max_error_p = maximum(error) # update error
        P0_hat_Z[:] = copy(P_hat_Z) # update to new price index
        iteration_p += 1 # update iteration count

        println(" - Inner loop: Iteration $iteration_p completed with error $max_error_p") # print update on progress

    end

    # ----------
    # price index for final goods from equation (49) --- not used in optimization, i.e. residual optimum (since cost is same in F as in Z!)
    cost_hat_F = zeros(N*S,N)
    for i in 1:N
        for r in 1:S
            for j in 1:N
                iirr = (i-1)*S+r
                cost_hat_F[iirr,j] = (cost_hat[r,i]*τ_hat_F[iirr,j])^(-θ[r,i])
            end
        end
    end

    global P_hat_F = zeros(S, N)
    for r in 1:S
        for j in 1:N
            P_hat_F[r,j] = (sum(π_F[r:S:(N-1)*S+r,j] .* cost_hat_F[r:S:(N-1)*S+r,j]))^(-1/θ[r,j])
        end
    end

    P_hat_F = ifelse.(isinf.(P_hat_F), 0.0, P_hat_F)

    # ----------
    # trade shares from equation (45) and (46)
    θ = vec(θ)
    θ = repeat(θ, 1, N*S)
    global π_hat_Z = cost_hat_Z ./ repeat(P_hat_Z, N) .^ (.-θ) # NS×NS, - cost_hat_Z is already to power of -θ => only P_hat_Z^-θ

    θ = θ[:, 1:N] # reduce to one NS×N again
    global π_hat_F = cost_hat_F ./ repeat(P_hat_F, N) .^ (.-θ) # NS×N,  - cost_hat_Z is already to power of -θ => only P_hat_Z^-θ

    return P_hat_Z, P_hat_F, π_hat_Z, π_hat_F, cost_hat
end


text:
# Command: Get Named Cookie
"""
    cookie(session::Session, cookie::AbstractString)::Cookie

Get Named Cookie.
"""
function cookie(session::Session, cookie::AbstractString)
    # x = "WhoIsAwesome"
    @unpack addr, id = session
    response = HTTP.get(
        "$addr/session/$id/cookie/$(escapeuri(cookie))",
        [("Content-Type" => "application/json")],
    )
    @assert response.status == 200
    JSON3.read(response.body).value
end

ext: jl


text:
#####
##### NetCDFOutputWriter tests
#####

function test_DateTime_netcdf_output(arch)
    grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
    clock = Clock(time=DateTime(2021, 1, 1))
    model = NonhydrostaticModel(architecture=arch, grid=grid, clock=clock)

    Δt = 5days + 3hours + 44.123seconds
    simulation = Simulation(model, Δt=Δt, stop_time=DateTime(2021, 2, 1))

    filepath = "test_DateTime.nc"
    simulation.output_writers[:cal] = NetCDFOutputWriter(model, fields(model), filepath=filepath, schedule=IterationInterval(1))

    run!(simulation)

    ds = NCDataset(filepath)
    @test ds["time"].attrib["units"] == "seconds since 2000-01-01 00:00:00"

    Nt = length(ds["time"])
    @test Nt == 8

    for n in 1:Nt-1
        @test ds["time"][n] == DateTime(2021, 1, 1) + (n-1) * Millisecond(1000Δt)
    end

    @test ds["time"][Nt] == DateTime(2021, 2, 1)

    close(ds)
    rm(filepath)

    return nothing
end

function test_TimeDate_netcdf_output(arch)
    grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
    clock = Clock(time=TimeDate(2021, 1, 1))
    model = NonhydrostaticModel(architecture=arch, grid=grid, clock=clock)

    Δt = 5days + 3hours + 44.123seconds
    simulation = Simulation(model, Δt=Δt, stop_time=TimeDate(2021, 2, 1))

    filepath = "test_TimeDate.nc"
    simulation.output_writers[:cal] = NetCDFOutputWriter(model, fields(model), filepath=filepath, schedule=IterationInterval(1))

    run!(simulation)

    ds = NCDataset(filepath)
    @test ds["time"].attrib["units"] == "seconds since 2000-01-01 00:00:00"

    Nt = length(ds["time"])
    @test Nt == 8

    for n in 1:Nt-1
        @test ds["time"][n] == DateTime(2021, 1, 1) + (n-1) * Millisecond(1000Δt)
    end

    @test ds["time"][Nt] == DateTime(2021, 2, 1)

    close(ds)
    rm(filepath)

    return nothing
end

function test_thermal_bubble_netcdf_output(arch)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100

    topo = (Periodic, Periodic, Bounded)
    grid = RegularRectilinearGrid(topology=topo, size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz))
    closure = IsotropicDiffusivity(ν=4e-2, κ=4e-2)
    model = NonhydrostaticModel(architecture=arch, grid=grid, closure=closure)
    simulation = Simulation(model, Δt=6, stop_iteration=10)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    CUDA.@allowscalar model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    outputs = Dict("v" => model.velocities.v,
                   "u" => model.velocities.u,
                   "w" => model.velocities.w,
                   "T" => model.tracers.T,
                   "S" => model.tracers.S)

    nc_filepath = "test_dump_$(typeof(arch)).nc"
    nc_writer = NetCDFOutputWriter(model, outputs, filepath=nc_filepath, schedule=IterationInterval(10), verbose=true)
    push!(simulation.output_writers, nc_writer)

    i_slice = 1:2:10
    j_slice = 13
    k_slice = 9:11
    field_slicer = FieldSlicer(i=i_slice, j=j_slice, k=k_slice)
    j_slice = j_slice:j_slice  # So we can correctly index with it for later tests.

    nc_sliced_filepath = "test_dump_sliced_$(typeof(arch)).nc"
    nc_sliced_writer = NetCDFOutputWriter(model, outputs, filepath=nc_sliced_filepath, schedule=IterationInterval(10),
                                          array_type=Array{Float32}, field_slicer=field_slicer, verbose=true)

    push!(simulation.output_writers, nc_sliced_writer)

    run!(simulation)

    ds3 = Dataset(nc_filepath)

    @test haskey(ds3.attrib, "date") && !isnothing(ds3.attrib["date"])
    @test haskey(ds3.attrib, "Julia") && !isnothing(ds3.attrib["Julia"])
    @test haskey(ds3.attrib, "Oceananigans") && !isnothing(ds3.attrib["Oceananigans"])
    @test haskey(ds3.attrib, "schedule") && ds3.attrib["schedule"] == "IterationInterval"
    @test haskey(ds3.attrib, "interval") && ds3.attrib["interval"] == 10
    @test haskey(ds3.attrib, "output iteration interval") && !isnothing(ds3.attrib["output iteration interval"])

    @test eltype(ds3["time"]) == eltype(model.clock.time)

    @test eltype(ds3["xC"]) == Float32
    @test eltype(ds3["xF"]) == Float32
    @test eltype(ds3["yC"]) == Float32
    @test eltype(ds3["yF"]) == Float32
    @test eltype(ds3["zC"]) == Float32
    @test eltype(ds3["zF"]) == Float32

    @test length(ds3["xC"]) == Nx
    @test length(ds3["yC"]) == Ny
    @test length(ds3["zC"]) == Nz
    @test length(ds3["xF"]) == Nx
    @test length(ds3["yF"]) == Ny
    @test length(ds3["zF"]) == Nz+1  # z is Bounded

    @test ds3["xC"][1] == grid.xC[1]
    @test ds3["xF"][1] == grid.xF[1]
    @test ds3["yC"][1] == grid.yC[1]
    @test ds3["yF"][1] == grid.yF[1]
    @test ds3["zC"][1] == grid.zC[1]
    @test ds3["zF"][1] == grid.zF[1]

    @test ds3["xC"][end] == grid.xC[Nx]
    @test ds3["xF"][end] == grid.xF[Nx]
    @test ds3["yC"][end] == grid.yC[Ny]
    @test ds3["yF"][end] == grid.yF[Ny]
    @test ds3["zC"][end] == grid.zC[Nz]
    @test ds3["zF"][end] == grid.zF[Nz+1]  # z is Bounded

    @test eltype(ds3["u"]) == Float32
    @test eltype(ds3["v"]) == Float32
    @test eltype(ds3["w"]) == Float32
    @test eltype(ds3["T"]) == Float32
    @test eltype(ds3["S"]) == Float32

    u = ds3["u"][:, :, :, end]
    v = ds3["v"][:, :, :, end]
    w = ds3["w"][:, :, :, end]
    T = ds3["T"][:, :, :, end]
    S = ds3["S"][:, :, :, end]

    close(ds3)

    @test all(u .≈ Array(interior(model.velocities.u)))
    @test all(v .≈ Array(interior(model.velocities.v)))
    @test all(w .≈ Array(interior(model.velocities.w)))
    @test all(T .≈ Array(interior(model.tracers.T)))
    @test all(S .≈ Array(interior(model.tracers.S)))

    ds2 = Dataset(nc_sliced_filepath)

    @test haskey(ds2.attrib, "date") && !isnothing(ds2.attrib["date"])
    @test haskey(ds2.attrib, "Julia") && !isnothing(ds2.attrib["Julia"])
    @test haskey(ds2.attrib, "Oceananigans") && !isnothing(ds2.attrib["Oceananigans"])
    @test haskey(ds2.attrib, "schedule") && ds2.attrib["schedule"] == "IterationInterval"
    @test haskey(ds2.attrib, "interval") && ds2.attrib["interval"] == 10
    @test haskey(ds2.attrib, "output iteration interval") && !isnothing(ds2.attrib["output iteration interval"])

    @test eltype(ds2["time"]) == eltype(model.clock.time)

    @test eltype(ds2["xC"]) == Float32
    @test eltype(ds2["xF"]) == Float32
    @test eltype(ds2["yC"]) == Float32
    @test eltype(ds2["yF"]) == Float32
    @test eltype(ds2["zC"]) == Float32
    @test eltype(ds2["zF"]) == Float32

    @test length(ds2["xC"]) == length(i_slice)
    @test length(ds2["xF"]) == length(i_slice)
    @test length(ds2["yC"]) == length(j_slice)
    @test length(ds2["yF"]) == length(j_slice)
    @test length(ds2["zC"]) == length(k_slice)
    @test length(ds2["zF"]) == length(k_slice)

    @test ds2["xC"][1] == grid.xC[i_slice[1]]
    @test ds2["xF"][1] == grid.xF[i_slice[1]]
    @test ds2["yC"][1] == grid.yC[j_slice[1]]
    @test ds2["yF"][1] == grid.yF[j_slice[1]]
    @test ds2["zC"][1] == grid.zC[k_slice[1]]
    @test ds2["zF"][1] == grid.zF[k_slice[1]]

    @test ds2["xC"][end] == grid.xC[i_slice[end]]
    @test ds2["xF"][end] == grid.xF[i_slice[end]]
    @test ds2["yC"][end] == grid.yC[j_slice[end]]
    @test ds2["yF"][end] == grid.yF[j_slice[end]]
    @test ds2["zC"][end] == grid.zC[k_slice[end]]
    @test ds2["zF"][end] == grid.zF[k_slice[end]]

    @test eltype(ds2["u"]) == Float32
    @test eltype(ds2["v"]) == Float32
    @test eltype(ds2["w"]) == Float32
    @test eltype(ds2["T"]) == Float32
    @test eltype(ds2["S"]) == Float32

    u_sliced = ds2["u"][:, :, :, end]
    v_sliced = ds2["v"][:, :, :, end]
    w_sliced = ds2["w"][:, :, :, end]
    T_sliced = ds2["T"][:, :, :, end]
    S_sliced = ds2["S"][:, :, :, end]

    close(ds2)

    @test all(u_sliced .≈ Array(interior(model.velocities.u))[i_slice, j_slice, k_slice])
    @test all(v_sliced .≈ Array(interior(model.velocities.v))[i_slice, j_slice, k_slice])
    @test all(w_sliced .≈ Array(interior(model.velocities.w))[i_slice, j_slice, k_slice])
    @test all(T_sliced .≈ Array(interior(model.tracers.T))[i_slice, j_slice, k_slice])
    @test all(S_sliced .≈ Array(interior(model.tracers.S))[i_slice, j_slice, k_slice])

    rm(nc_filepath)
    rm(nc_sliced_filepath)

    return nothing
end

function test_thermal_bubble_netcdf_output_with_halos(arch)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100

    topo = (Periodic, Periodic, Bounded)
    grid = RegularRectilinearGrid(topology=topo, size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz))
    closure = IsotropicDiffusivity(ν=4e-2, κ=4e-2)
    model = NonhydrostaticModel(architecture=arch, grid=grid, closure=closure)
    simulation = Simulation(model, Δt=6, stop_iteration=10)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    CUDA.@allowscalar model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    nc_filepath = "test_dump_with_halos_$(typeof(arch)).nc"
    nc_writer = NetCDFOutputWriter(model, merge(model.velocities, model.tracers),
                                   filepath=nc_filepath,
                                   schedule=IterationInterval(10),
                                   field_slicer=FieldSlicer(with_halos=true))

    push!(simulation.output_writers, nc_writer)

    run!(simulation)

    ds = Dataset(nc_filepath)

    @test haskey(ds.attrib, "date") && !isnothing(ds.attrib["date"])
    @test haskey(ds.attrib, "Julia") && !isnothing(ds.attrib["Julia"])
    @test haskey(ds.attrib, "Oceananigans") && !isnothing(ds.attrib["Oceananigans"])
    @test haskey(ds.attrib, "schedule") && ds.attrib["schedule"] == "IterationInterval"
    @test haskey(ds.attrib, "interval") && ds.attrib["interval"] == 10
    @test haskey(ds.attrib, "output iteration interval") && !isnothing(ds.attrib["output iteration interval"])

    @test eltype(ds["time"]) == eltype(model.clock.time)

    @test eltype(ds["xC"]) == Float32
    @test eltype(ds["xF"]) == Float32
    @test eltype(ds["yC"]) == Float32
    @test eltype(ds["yF"]) == Float32
    @test eltype(ds["zC"]) == Float32
    @test eltype(ds["zF"]) == Float32

    Hx, Hy, Hz = grid.Hx, grid.Hy, grid.Hz
    @test length(ds["xC"]) == Nx+2Hx
    @test length(ds["yC"]) == Ny+2Hy
    @test length(ds["zC"]) == Nz+2Hz
    @test length(ds["xF"]) == Nx+2Hx
    @test length(ds["yF"]) == Ny+2Hy
    @test length(ds["zF"]) == Nz+2Hz+1  # z is Bounded

    @test ds["xC"][1] == grid.xC[1-Hx]
    @test ds["xF"][1] == grid.xF[1-Hx]
    @test ds["yC"][1] == grid.yC[1-Hy]
    @test ds["yF"][1] == grid.yF[1-Hy]
    @test ds["zC"][1] == grid.zC[1-Hz]
    @test ds["zF"][1] == grid.zF[1-Hz]

    @test ds["xC"][end] == grid.xC[Nx+Hx]
    @test ds["xF"][end] == grid.xF[Nx+Hx]
    @test ds["yC"][end] == grid.yC[Ny+Hy]
    @test ds["yF"][end] == grid.yF[Ny+Hy]
    @test ds["zC"][end] == grid.zC[Nz+Hz]
    @test ds["zF"][end] == grid.zF[Nz+Hz+1]  # z is Bounded

    @test eltype(ds["u"]) == Float32
    @test eltype(ds["v"]) == Float32
    @test eltype(ds["w"]) == Float32
    @test eltype(ds["T"]) == Float32
    @test eltype(ds["S"]) == Float32

    u = ds["u"][:, :, :, end]
    v = ds["v"][:, :, :, end]
    w = ds["w"][:, :, :, end]
    T = ds["T"][:, :, :, end]
    S = ds["S"][:, :, :, end]

    close(ds)

    @test all(u .≈ Array(model.velocities.u.data.parent))
    @test all(v .≈ Array(model.velocities.v.data.parent))
    @test all(w .≈ Array(model.velocities.w.data.parent))
    @test all(T .≈ Array(model.tracers.T.data.parent))
    @test all(S .≈ Array(model.tracers.S.data.parent))

    rm(nc_filepath)

    return nothing
end

function test_netcdf_function_output(arch)
    N = 16
    L = 1
    Δt = 1.25
    iters = 3

    grid = RegularRectilinearGrid(size=(N, N, N), extent=(L, 2L, 3L))
    model = NonhydrostaticModel(architecture=arch, grid=grid)
    simulation = Simulation(model, Δt=Δt, stop_iteration=iters)
    grid = model.grid

    # Define scalar, vector, and 2D slice outputs
    f(model) = model.clock.time^2

    g(model) = model.clock.time .* exp.(znodes(Center, grid))

    h(model) = model.clock.time .* (   sin.(xnodes(Center, grid, reshape=true)[:, :, 1])
                                    .* cos.(ynodes(Face, grid, reshape=true)[:, :, 1]))

    outputs = (scalar=f, profile=g, slice=h)
    dims = (scalar=(), profile=("zC",), slice=("xC", "yC"))

    output_attributes = (
        scalar = (longname="Some scalar", units="bananas"),
        profile = (longname="Some vertical profile", units="watermelons"),
        slice = (longname="Some slice", units="mushrooms")
    )

    global_attributes = (location="Bay of Fundy", onions=7)

    nc_filepath = "test_function_outputs_$(typeof(arch)).nc"

    simulation.output_writers[:food] =
        NetCDFOutputWriter(model, outputs; filepath=nc_filepath,
                           schedule=TimeInterval(Δt), dimensions=dims, array_type=Array{Float64}, verbose=true,
                           global_attributes=global_attributes, output_attributes=output_attributes)

    run!(simulation)

    ds = Dataset(nc_filepath, "r")

    @test haskey(ds.attrib, "date") && !isnothing(ds.attrib["date"])
    @test haskey(ds.attrib, "Julia") && !isnothing(ds.attrib["Julia"])
    @test haskey(ds.attrib, "Oceananigans") && !isnothing(ds.attrib["Oceananigans"])
    @test haskey(ds.attrib, "schedule") && !isnothing(ds.attrib["schedule"])
    @test haskey(ds.attrib, "interval") && !isnothing(ds.attrib["interval"])
    @test haskey(ds.attrib, "output time interval") && !isnothing(ds.attrib["output time interval"])

    @test eltype(ds["time"]) == eltype(model.clock.time)

    @test eltype(ds["xC"]) == Float64
    @test eltype(ds["xF"]) == Float64
    @test eltype(ds["yC"]) == Float64
    @test eltype(ds["yF"]) == Float64
    @test eltype(ds["zC"]) == Float64
    @test eltype(ds["zF"]) == Float64

    @test length(ds["xC"]) == N
    @test length(ds["yC"]) == N
    @test length(ds["zC"]) == N
    @test length(ds["xF"]) == N
    @test length(ds["yF"]) == N
    @test length(ds["zF"]) == N+1  # z is Bounded

    @test ds["xC"][1] == grid.xC[1]
    @test ds["xF"][1] == grid.xF[1]
    @test ds["yC"][1] == grid.yC[1]
    @test ds["yF"][1] == grid.yF[1]
    @test ds["zC"][1] == grid.zC[1]
    @test ds["zF"][1] == grid.zF[1]

    @test ds["xC"][end] == grid.xC[N]
    @test ds["yC"][end] == grid.yC[N]
    @test ds["zC"][end] == grid.zC[N]
    @test ds["xF"][end] == grid.xF[N]
    @test ds["yF"][end] == grid.yF[N]
    @test ds["zF"][end] == grid.zF[N+1]  # z is Bounded

    @test ds.attrib["location"] == "Bay of Fundy"
    @test ds.attrib["onions"] == 7

    @test eltype(ds["scalar"]) == Float64
    @test eltype(ds["profile"]) == Float64
    @test eltype(ds["slice"]) == Float64

    @test length(ds["time"]) == iters+1
    @test ds["time"][:] == [n*Δt for n in 0:iters]

    @test length(ds["scalar"]) == iters+1
    @test ds["scalar"].attrib["longname"] == "Some scalar"
    @test ds["scalar"].attrib["units"] == "bananas"
    @test ds["scalar"][:] == [(n*Δt)^2 for n in 0:iters]
    @test dimnames(ds["scalar"]) == ("time",)

    @test ds["profile"].attrib["longname"] == "Some vertical profile"
    @test ds["profile"].attrib["units"] == "watermelons"
    @test size(ds["profile"]) == (N, iters+1)
    @test dimnames(ds["profile"]) == ("zC", "time")

    for n in 0:iters
        @test ds["profile"][:, n+1] == n*Δt .* exp.(znodes(Center, grid))
    end

    @test ds["slice"].attrib["longname"] == "Some slice"
    @test ds["slice"].attrib["units"] == "mushrooms"
    @test size(ds["slice"]) == (N, N, iters+1)
    @test dimnames(ds["slice"]) == ("xC", "yC", "time")

    for n in 0:iters
        @test ds["slice"][:, :, n+1] == n*Δt .* (   sin.(xnodes(Center, grid, reshape=true)[:, :, 1])
                                                 .* cos.(ynodes(Face, grid, reshape=true)[:, :, 1]))
    end

    close(ds)

    #####
    ##### Take 1 more time step and test that appending to a NetCDF file works
    #####

    iters += 1
    simulation = Simulation(model, Δt=Δt, stop_iteration=iters)

    simulation.output_writers[:food] =
        NetCDFOutputWriter(model, outputs; filepath=nc_filepath, mode="a",
                           schedule=IterationInterval(1), array_type=Array{Float64}, dimensions=dims, verbose=true,
                           global_attributes=global_attributes, output_attributes=output_attributes)

    run!(simulation)

    ds = Dataset(nc_filepath, "r")

    @test length(ds["time"]) == iters+1
    @test length(ds["scalar"]) == iters+1
    @test size(ds["profile"]) == (N, iters+1)
    @test size(ds["slice"]) == (N, N, iters+1)

    @test ds["time"][:] == [n*Δt for n in 0:iters]
    @test ds["scalar"][:] == [(n*Δt)^2 for n in 0:iters]

    for n in 0:iters
        @test ds["profile"][:, n+1] == n*Δt .* exp.(znodes(Center, grid))
        @test ds["slice"][:, :, n+1] == n*Δt .* (   sin.(xnodes(Center, grid, reshape=true)[:, :, 1])
                                                 .* cos.(ynodes(Face, grid, reshape=true)[:, :, 1]))
    end

    close(ds)

    rm(nc_filepath)

    return nothing
end

function test_netcdf_time_averaging(arch)
    topo = (Periodic, Periodic, Periodic)
    domain = (x=(0, 1), y=(0, 1), z=(0, 1))
    grid = RegularRectilinearGrid(topology=topo, size=(4, 4, 4); domain...)

    λ1(x, y, z) = x + (1 - y)^2 + tanh(z)
    λ2(x, y, z) = x + (1 - y)^2 + tanh(4z)

    Fc1(x, y, z, t, c1) = - λ1(x, y, z) * c1
    Fc2(x, y, z, t, c2) = - λ2(x, y, z) * c2

    c1_forcing = Forcing(Fc1, field_dependencies=:c1)
    c2_forcing = Forcing(Fc2, field_dependencies=:c2)

    model = NonhydrostaticModel(
                grid = grid,
        architecture = arch,
         timestepper = :RungeKutta3,
             tracers = (:c1, :c2),
             forcing = (c1=c1_forcing, c2=c2_forcing),
            coriolis = nothing,
            buoyancy = nothing,
             closure = nothing
    )

    set!(model, c1=1, c2=1)

    Δt = 1/512  # Nice floating-point number
    simulation = Simulation(model, Δt=Δt, stop_time=50Δt)

    ∫c1_dxdy = AveragedField(model.tracers.c1, dims=(1, 2))
    ∫c2_dxdy = AveragedField(model.tracers.c2, dims=(1, 2))

    nc_outputs = Dict("c1" => ∫c1_dxdy, "c2" => ∫c2_dxdy)
    nc_dimensions = Dict("c1" => ("zC",), "c2" => ("zC",))

    horizontal_average_nc_filepath = "decay_averaged_field_test.nc"
    simulation.output_writers[:horizontal_average] =
        NetCDFOutputWriter(model, nc_outputs, filepath=horizontal_average_nc_filepath, schedule=TimeInterval(10Δt),
                           dimensions=nc_dimensions, array_type=Array{Float64}, verbose=true)

    multiple_time_average_nc_filepath = "decay_windowed_time_average_test.nc"
    single_time_average_nc_filepath = "single_decay_windowed_time_average_test.nc"
    window = 6Δt
    stride = 2

    single_nc_output = Dict("c1" => ∫c1_dxdy)
    single_nc_dimension = Dict("c1" => ("zC",))

    simulation.output_writers[:single_output_time_average] =
        NetCDFOutputWriter(model, single_nc_output, filepath=single_time_average_nc_filepath, array_type=Array{Float64},
                           schedule=AveragedTimeInterval(10Δt, window=window, stride=stride),
                           dimensions=single_nc_dimension, verbose=true)

    simulation.output_writers[:multiple_output_time_average] =
        NetCDFOutputWriter(model, nc_outputs, filepath=multiple_time_average_nc_filepath, array_type=Array{Float64},
                           schedule=AveragedTimeInterval(10Δt, window=window, stride=stride),
                           dimensions=nc_dimensions, verbose=true)

    run!(simulation)

    ##### For each λ, horizontal average should evaluate to
    #####
    #####     c̄(z, t) = ∫₀¹ ∫₀¹ exp{- λ(x, y, z) * t} dx dy
    #####             = 1 / (Nx*Ny) * Σᵢ₌₁ᴺˣ Σⱼ₌₁ᴺʸ exp{- λ(i, j, k) * t}
    #####
    ##### which we can compute analytically.

    ds = NCDataset(horizontal_average_nc_filepath)

    Nx, Ny, Nz = size(grid)
    xs, ys, zs = nodes(model.tracers.c1)

    c̄1(z, t) = 1 / (Nx * Ny) * sum(exp(-λ1(x, y, z) * t) for x in xs for y in ys)
    c̄2(z, t) = 1 / (Nx * Ny) * sum(exp(-λ2(x, y, z) * t) for x in xs for y in ys)

    for (n, t) in enumerate(ds["time"])
        @test ds["c1"][:, n] ≈ c̄1.(zs, t)
        @test ds["c2"][:, n] ≈ c̄2.(zs, t)
    end

    close(ds)

    # Compute time averages...
    c̄1(ts) = 1/length(ts) * sum(c̄1.(zs, t) for t in ts)
    c̄2(ts) = 1/length(ts) * sum(c̄2.(zs, t) for t in ts)

    #####
    ##### Test strided windowed time average against analytic solution
    ##### for *single* NetCDF output
    #####

    single_ds = NCDataset(single_time_average_nc_filepath)

    attribute_names = ("schedule", "interval", "output time interval",
                       "time_averaging_window", "time averaging window",
                       "time_averaging_stride", "time averaging stride")

    for name in attribute_names
        @test haskey(single_ds.attrib, name) && !isnothing(single_ds.attrib[name])
    end

    window_size = Int(window/Δt)

    @info "    Testing time-averaging of a single NetCDF output [$(typeof(arch))]..."

    for (n, t) in enumerate(single_ds["time"][2:end])
        averaging_times = [t - n*Δt for n in 0:stride:window_size-1]
        @test single_ds["c1"][:, n+1] ≈ c̄1(averaging_times)
    end

    close(single_ds)

    #####
    ##### Test strided windowed time average against analytic solution
    ##### for *multiple* NetCDF outputs
    #####

    ds = NCDataset(multiple_time_average_nc_filepath)

    @info "    Testing time-averaging of multiple NetCDF outputs [$(typeof(arch))]..."

    for (n, t) in enumerate(ds["time"][2:end])
        averaging_times = [t - n*Δt for n in 0:stride:window_size-1]
        @test ds["c1"][:, n+1] ≈ c̄1(averaging_times)
        @test ds["c2"][:, n+1] ≈ c̄2(averaging_times)
    end

    close(ds)

    rm(horizontal_average_nc_filepath)
    rm(multiple_time_average_nc_filepath)

    return nothing
end

function test_netcdf_output_alignment(arch)
    grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(architecture=arch, grid=grid)
    simulation = Simulation(model, Δt=0.2, stop_time=40)

    test_filename1 = "test_output_alignment1.nc"
    simulation.output_writers[:stuff] =
        NetCDFOutputWriter(model, model.velocities, filepath=test_filename1,
                           schedule=TimeInterval(7.3))

    test_filename2 = "test_output_alignment2.nc"
    simulation.output_writers[:something] =
        NetCDFOutputWriter(model, model.tracers, filepath=test_filename2,
                           schedule=TimeInterval(3.0))

    run!(simulation)

    Dataset(test_filename1, "r") do ds
        @test all(ds["time"] .== 0:7.3:40)
    end

    Dataset(test_filename2, "r") do ds
        @test all(ds["time"] .== 0:3.0:40)
    end

    rm(test_filename1)
    rm(test_filename2)

    return nothing
end

function test_netcdf_vertically_stretched_grid_output(arch)
    Nx = Ny = 8
    Nz = 16
    zF = [k^2 for k in 0:Nz]
    grid = VerticallyStretchedRectilinearGrid(architecture=arch, size=(Nx, Ny, Nz), x=(0, 1), y=(-π, π), z_faces=zF)

    model = NonhydrostaticModel(architecture=arch, grid=grid)

    Δt = 1.25
    iters = 3
    simulation = Simulation(model, Δt=Δt, stop_iteration=iters)

    nc_filepath = "test_netcdf_vertically_stretched_grid_output_$(typeof(arch)).nc"

    simulation.output_writers[:fields] =
        NetCDFOutputWriter(model, merge(model.velocities, model.tracers),
                             filepath = nc_filepath,
                             schedule = IterationInterval(1),
                           array_type = Array{Float64},
                              verbose = true)

    run!(simulation)

    grid = model.grid
    ds = NCDataset(nc_filepath)

    @test length(ds["xC"]) == Nx
    @test length(ds["yC"]) == Ny
    @test length(ds["zC"]) == Nz
    @test length(ds["xF"]) == Nx
    @test length(ds["yF"]) == Ny
    @test length(ds["zF"]) == Nz+1  # z is Bounded

    @test ds["xC"][1] == grid.xᶜᵃᵃ[1]
    @test ds["xF"][1] == grid.xᶠᵃᵃ[1]
    @test ds["yC"][1] == grid.yᵃᶜᵃ[1]
    @test ds["yF"][1] == grid.yᵃᶠᵃ[1]

    @test CUDA.@allowscalar ds["zC"][1] == grid.zᵃᵃᶜ[1]
    @test CUDA.@allowscalar ds["zF"][1] == grid.zᵃᵃᶠ[1]

    @test ds["xC"][end] == grid.xᶜᵃᵃ[Nx]
    @test ds["xF"][end] == grid.xᶠᵃᵃ[Nx]
    @test ds["yC"][end] == grid.yᵃᶜᵃ[Ny]
    @test ds["yF"][end] == grid.yᵃᶠᵃ[Ny]

    @test CUDA.@allowscalar  ds["zC"][end] == grid.zᵃᵃᶜ[Nz]
    @test CUDA.@allowscalar  ds["zF"][end] == grid.zᵃᵃᶠ[Nz+1]  # z is Bounded

    close(ds)
    rm(nc_filepath)

    return nothing
end

using Oceananigans.Models.HydrostaticFreeSurfaceModels: VectorInvariant

function test_netcdf_regular_lat_lon_grid_output(arch)
    Nx = Ny = Nz = 16
    grid = RegularLatitudeLongitudeGrid(size=(Nx, Ny, Nz), longitude=(-180, 180), latitude=(-80, 80), z=(-100, 0))
    model = HydrostaticFreeSurfaceModel(architecture=arch, momentum_advection = VectorInvariant(), grid=grid)

    Δt = 1.25
    iters = 3
    simulation = Simulation(model, Δt=Δt, stop_iteration=iters)

    nc_filepath = "test_netcdf_regular_lat_lon_grid_output_$(typeof(arch)).nc"

    simulation.output_writers[:fields] =
        NetCDFOutputWriter(model, merge(model.velocities, model.tracers),
                             filepath = nc_filepath,
                             schedule = IterationInterval(1),
                           array_type = Array{Float64},
                              verbose = true)

    run!(simulation)

    grid = model.grid
    ds = NCDataset(nc_filepath)

    @test length(ds["xC"]) == Nx
    @test length(ds["yC"]) == Ny
    @test length(ds["zC"]) == Nz
    @test length(ds["xF"]) == Nx
    @test length(ds["yF"]) == Ny+1  # y is Bounded
    @test length(ds["zF"]) == Nz+1  # z is Bounded

    @test ds["xC"][1] == grid.λᶜᵃᵃ[1]
    @test ds["xF"][1] == grid.λᶠᵃᵃ[1]
    @test ds["yC"][1] == grid.φᵃᶜᵃ[1]
    @test ds["yF"][1] == grid.φᵃᶠᵃ[1]
    @test ds["zC"][1] == grid.zᵃᵃᶜ[1]
    @test ds["zF"][1] == grid.zᵃᵃᶠ[1]

    @test ds["xC"][end] == grid.λᶜᵃᵃ[Nx]
    @test ds["xF"][end] == grid.λᶠᵃᵃ[Nx]
    @test ds["yC"][end] == grid.φᵃᶜᵃ[Ny]
    @test ds["yF"][end] == grid.φᵃᶠᵃ[Ny+1]  # y is Bounded
    @test ds["zC"][end] == grid.zᵃᵃᶜ[Nz]
    @test ds["zF"][end] == grid.zᵃᵃᶠ[Nz+1]  # z is Bounded

    close(ds)
    rm(nc_filepath)

    return nothing
end

for arch in archs
    @testset "NetCDF output writer [$(typeof(arch))]" begin
        @info "  Testing NetCDF output writer [$(typeof(arch))]..."
        test_DateTime_netcdf_output(arch)
        test_TimeDate_netcdf_output(arch)
        test_thermal_bubble_netcdf_output(arch)
        test_thermal_bubble_netcdf_output_with_halos(arch)
        test_netcdf_function_output(arch)
        test_netcdf_output_alignment(arch)
        test_netcdf_time_averaging(arch)
        test_netcdf_vertically_stretched_grid_output(arch)
        test_netcdf_regular_lat_lon_grid_output(arch)
    end
end


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "../.."))

using Oceananigans
using Plots, LinearAlgebra, Distributions, JLD2
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity
using OceanTurbulenceParameterEstimation

include("utils/lesbrary_paths.jl")
include("utils/one_dimensional_ensemble_model.jl")
include("utils/parameters.jl")
include("utils/visualize_profile_predictions.jl")

# two_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/2DaySuite"
# four_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/4DaySuite"
# six_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/6DaySuite"

# two_day_suite = TwoDaySuite(two_day_suite_dir)
# four_day_suite = FourDaySuite(four_day_suite_dir)
# six_day_suite = SixDaySuite(six_day_suite_dir)

# calibration = InverseProblem(two_day_suite, parameters; relative_weights = relative_weight_options["all_but_e"],
#                              architecture = GPU(), ensemble_size = 10, Δt = 30.0)

# validation = InverseProblem(four_day_suite, calibration; Nz = 32);

#####
##### Set up ensemble model
#####

## NEED TO IMPLEMENT COARSE-GRAINING
##
##

lesbrary_directory = "/Users/adelinehillier/Desktop/dev/"

observations = TwoDaySuite(lesbrary_directory; first_iteration = 13, last_iteration = nothing, normalize = ZScore, Nz = 128)

parameter_set = CATKEParametersRiDependent
closure = closure_with_parameter_set(CATKEVerticalDiffusivity(Float64;), parameter_set)

ensemble_model = OneDimensionalEnsembleModel(observations;
    architecture = CPU(),
    ensemble_size = 20,
    closure = closure
)

ensemble_simulation = Simulation(ensemble_model; Δt = 10seconds, stop_time = 2days)

pop!(ensemble_simulation.diagnostics, :nan_checker)

#####
##### Build free parameters
#####

free_parameter_names = keys(parameter_set.defaults)
free_parameter_means = collect(values(parameter_set.defaults))
priors = NamedTuple(pname => ConstrainedNormal(0.0, 1.0, bounds(pname) .* 0.5...) for pname in free_parameter_names)

free_parameters = FreeParameters(priors)

#####
##### Build the Inverse Problem
#####

calibration = InverseProblem(observations, ensemble_simulation, free_parameters);

# #####
# ##### Calibrate
# #####

iterations = 5
# eki = EnsembleKalmanInversion(calibration; noise_covariance = 1e-2)
# params = iterate!(eki; iterations = iterations)

directory =
# visualize!(calibration, params;
#     field_names = [:u, :v, :b, :e],
#     directory = @__DIR__,
#     filename = "perfect_model_visual_calibrated.png"
# )
# @show params


text:
module Reader

export

    readdata!,
    readdata

const NROW = 32
const NCOL = 32
const NCHAN = 3
const NBYTE = NROW * NCOL * NCHAN + 2 # "+ 2" for label

function readnext!(buffer::Array{UInt8}, io::IO)
    c = Int(read(io, UInt8))
    f = Int(read(io, UInt8))
    read!(io, buffer)
    buffer, c, f
end

function readdata!(buffer::Array{UInt8}, io::IO, index::Integer)
    seek(io, (index - 1) * NBYTE)
    readnext!(buffer, io)
end

function readdata(io::IO, nobs::Int, index::Integer)
    buffer = Array{UInt8}(undef, NROW, NCOL, NCHAN)
    readdata!(buffer, io, index)
end

function readdata(io::IO, nobs::Int)
    X = Array{UInt8}(undef, NROW, NCOL, NCHAN, nobs)
    C = Array{Int}(undef, nobs)
    F = Array{Int}(undef, nobs)
    buffer = Array{UInt8}(undef, NROW, NCOL, NCHAN)
    @inbounds for index in 1:nobs
        _, tc, tf = readnext!(buffer, io)
        copyto!(view(X,:,:,:,index), buffer)
        C[index] = tc
        F[index] = tf
    end
    X, C, F
end

function readdata(file::AbstractString, nobs::Int, index::Integer)
    open(file, "r") do io
        readdata(io, nobs, index)
    end::Tuple{Array{UInt8,3},Int,Int}
end

function readdata(file::AbstractString, nobs::Int)
    open(file, "r") do io
        readdata(io, nobs)
    end::Tuple{Array{UInt8,4},Vector{Int},Vector{Int}}
end

end


text:
module TreeReducer
    export tmapreduce

    function tmapreduce(f::Function, op::Function, workers::Array{Int}=procs())
        futures = [@spawnat w f() for w in workers]

        return mapreduce(fetch, op, futures)
    end
end # module

max_stars_repo_name: gueltro/TreeReducer

max_issues_repo_name: gueltro/TreeReducer

max_forks_repo_name: gueltro/TreeReducer


text:
# basic graph tests for trivial cases

using Test
using Statistics
using IncrementalInference

@testset "Test basic single variable graph with one prior..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))

tree, smt, hist = solveTree!(fg)

# check mean and covariance
@test (getKDE(fg, :x0) |> getKDEMean .|> abs)[1] < 0.5
@test 0.4 < Statistics.cov( getPoints(getKDE(fg, :x0))[1,:] ) < 1.8

end

@testset "Test basic single variable graph with one prior offset by 1000..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(1000.0,1.0)))

tree, smt, hist = solveTree!(fg)

# check mean and covariance
@test abs((getKDE(fg, :x0) |> getKDEMean)[1]-1000) < 0.5
@test 0.4 < Statistics.cov( getPoints(getKDE(fg, :x0))[1,:] ) < 1.8

end

@testset "Test basic single variable graph with two identical priors..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))

tree, smt, hist = solveTree!(fg)

# check mean and covariance
@test (getKDE(fg, :x0) |> getKDEMean .|> abs)[1] < 0.3
# should be sqrt(0.5) = 0.7, but lands near 0.6 instead -- computation is too confident.
@test 0.3 < Statistics.cov( getPoints(getKDE(fg, :x0))[1,:] ) < 1.0

end

@testset "Test basic single variable graph with three identical priors..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))

tree, smt, hist = solveTree!(fg)

# check mean and covariance
@test (getKDE(fg, :x0) |> getKDEMean .|> abs)[1] < 0.4
# should be sqrt(1/3) = 0.577, but lands near 0.35 instead -- computation is too confident.
@test 0.1 < Statistics.cov( getPoints(getKDE(fg, :x0))[1,:] ) < 0.75

end

@testset "Test basic single variable graph with two priors at + and - 1..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(-1.0,1.0)))
addFactor!(fg, [:x0;], Prior(Normal(+1.0,1.0)))

tree, smt, hist = solveTree!(fg)

# check mean and covariance -- should be zero
@test (getKDE(fg, :x0) |> getKDEMean .|> abs)[1] < 0.8
# should be sqrt(1/2) = 0.707 -- computation results nearer 0.7.
@test 0.2 < Statistics.cov( getPoints(getKDE(fg, :x0))[1,:] ) < 1

end

@testset "Test basic single variable graph with two priors at + and - 1, offset by -1000..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(-1.0-1000,1.0)))
addFactor!(fg, [:x0;], Prior(Normal(+1.0-1000,1.0)))

tree, smt, hist = solveTree!(fg)

# check mean and covariance -- should be zero
@test abs((getKDE(fg, :x0) |> getKDEMean)[1] + 1000) < 0.6
# should be sqrt(1/2) = 0.707 -- computation results nearer 0.7.
@test 0.2 < Statistics.cov( getPoints(getKDE(fg, :x0))[1,:] ) < 1

end

@testset "Test basic two variable graph with two identical priors and weak connection..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))
addFactor!(fg, [:x1;], Prior(Normal(0.0,1.0)))
addFactor!(fg, [:x0;:x1;], LinearConditional(Normal(0.0,10.0)))

tree, smt, hist = solveTree!(fg)

# check mean and covariance -- should be zero
@test (getKDE(fg, :x0) |> getKDEMean .|> abs)[1] < 0.6
@test (getKDE(fg, :x1) |> getKDEMean .|> abs)[1] < 0.6

@test 0.4 < Statistics.cov( getPoints(getKDE(fg, :x0))[1,:] ) < 2.2
@test 0.4 < Statistics.cov( getPoints(getKDE(fg, :x1))[1,:] ) < 2.3

end

@testset "Test basic two variable graph with two separated priors and weak connection..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(-1.0,1.0)))
addFactor!(fg, [:x1;], Prior(Normal(+1.0,1.0)))
addFactor!(fg, [:x0;:x1;], LinearConditional(Normal(0.0,10.0)))

tree, smt, hist = solveTree!(fg)

# check mean and covariance -- should be near each prior
@test abs((getKDE(fg, :x0) |> getKDEMean)[1]+1) < 0.5
@test abs((getKDE(fg, :x1) |> getKDEMean)[1]-1) < 0.5

@test 0.4 < Statistics.cov( getPoints(getKDE(fg, :x0))[1,:] ) < 1.8
@test 0.4 < Statistics.cov( getPoints(getKDE(fg, :x1))[1,:] ) < 1.8

end

@testset "Test basic two variable graph with two separated priors and strong connection..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(-1.0,1.0)))
addFactor!(fg, [:x2;], Prior(Normal(+1.0,1.0)))
addFactor!(fg, [:x0;:x1;], LinearConditional(Normal(0.0,1.0)))
addFactor!(fg, [:x1;:x2;], LinearConditional(Normal(0.0,1.0)))

tree, smt, hist = solveTree!(fg)

# check mean and covariance -- should between two priors somewhere
@test abs((getKDE(fg, :x0) |> getKDEMean)[1] + 1) < 0.8
@test abs((getKDE(fg, :x1) |> getKDEMean)[1]) < 0.9
@test abs((getKDE(fg, :x2) |> getKDEMean)[1] - 1) < 0.9

@test 0.3 < Statistics.cov( getPoints(getKDE(fg, :x0))[1,:] ) < 1.8
@test 0.3 < Statistics.cov( getPoints(getKDE(fg, :x1))[1,:] ) < 1.8
@test 0.3 < Statistics.cov( getPoints(getKDE(fg, :x2))[1,:] ) < 1.8

end

@testset "Test basic five variable graph with two separated priors and nominal connection..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addVariable!(fg, :x3, ContinuousScalar)
addVariable!(fg, :x4, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(-3.0,1.0)))
addFactor!(fg, [:x4;], Prior(Normal(+3.0,1.0)))
addFactor!(fg, [:x0;:x1;], LinearConditional(Normal(0.0,1.0)))
addFactor!(fg, [:x1;:x2;], LinearConditional(Normal(0.0,1.0)))
addFactor!(fg, [:x2;:x3;], LinearConditional(Normal(0.0,1.0)))
addFactor!(fg, [:x3;:x4;], LinearConditional(Normal(0.0,1.0)))

# drawGraph(fg, show=true)

tree, smt, hist = solveTree!(fg)

# using KernelDensityEstimatePlotting
# plotKDE((x->getKDE(fg,x)).([:x0;:x1;:x2;:x3;:x4]))
# using Gadfly, Cairo, Fontconfig
# drawTree(tree,show=true,imgs=true)

# check mean and covariance -- should be zero
X0 = (getKDE(fg, :x0) |> getKDEMean)[1]
X1 = (getKDE(fg, :x1) |> getKDEMean)[1]
X2 = (getKDE(fg, :x2) |> getKDEMean)[1]
X3 = (getKDE(fg, :x3) |> getKDEMean)[1]
X4 = (getKDE(fg, :x4) |> getKDEMean)[1]

@test X0 < X1 < X2 < X3 < X4

@test abs(X0+X4) < 2.2
@test abs(X1+X3) < 2.2
@test abs(X2) < 2.2

@test 0.2 < Statistics.cov( getPoints(getKDE(fg, :x0))[1,:] ) < 2.3
@test 0.2 < Statistics.cov( getPoints(getKDE(fg, :x1))[1,:] ) < 2.3
@test 0.2 < Statistics.cov( getPoints(getKDE(fg, :x2))[1,:] ) < 2.3
@test 0.2 < Statistics.cov( getPoints(getKDE(fg, :x3))[1,:] ) < 2.3
@test 0.2 < Statistics.cov( getPoints(getKDE(fg, :x4))[1,:] ) < 2.3

end

#


text:
"""
Tools for computing the formal solution of the radiative transfer equation.
"""

"""
Compute intensity by "brute force" trapezoidal integration. Doesn't work well in many
cases, its use is discouraged.
"""
function calc_intensity_brute_force(
    distance::Array{<:Unitful.Length, 1},
    extinction::Array{<:PerLength, 1},
    source_function::Array{<:UnitsIntensity_λ, 1}
)
    @assert distance[2] > distance[1] "Distance must be monotonically increasing"
    # Since integration functions don't work with units,
    # need to ensure quantities are in compatible units
    dist = ustrip.(distance .|> u"m")
    ext = ustrip.(extinction .|> u"m^-1")
    source = ustrip.(source_function .|> u"kW / (m^2 * nm)")
    τ = cumul_integrate(dist, ext, TrapezoidalFast())
    return integrate(τ, source .* exp.(-τ))u"kW / (m^2 * nm)"
end

"""
    piecewise_1D_nn(
        z::Array{<:Unitful.Length, 1},
        α::Array{<:PerLength, 1},
        source_function::Array{<:Unitful.Quantity, 1};
        to_end::Bool=false,
        initial_condition=:source,
    )

Compute piecewise integration of the radiative transfer equation,
assuming nearest-neighbour integration of the source function, for a given
height `z`, extinction `α` and `source_function`. The optional
keyword argument `to_end` defines the direction of the integration:
if `false` (default) will start to integrate intensity from the last
element to the first, and if `true` will integrate from the first
element to the last. `initial_condition` can take two values: `:zero` for
no radiation, or `:source` (default) to take the source function at the
starting point.
"""
function piecewise_1D_nn(
    z::Array{<:Unitful.Length, 1},
    α::Array{<:PerLength, 1},
    source_function::Array{<:Unitful.Quantity, 1};
    to_end::Bool=false,
    initial_condition=:source
)
    ndep = length(z)
    if to_end
        start = 1
        incr = 1
        depth_range = 2:ndep
    else
        start = ndep
        incr = -1
        depth_range = ndep-1:-1:1
    end
    intensity = similar(source_function)
    if initial_condition == :source
        intensity[start] = source_function[start]
    elseif initial_condition == :zero
        intensity[start] *= 0
    else
        error("NotImplemented initial condition $initial_condition")
    end
    for i in depth_range
        Δτ = abs(z[i] - z[i-incr]) * (α[i] + α[i-incr]) / 2
        w, _ = _w2(Δτ)
        intensity[i] = ((1 - w)*intensity[i-incr] +
                         w * (source_function[i] + source_function[i-incr]) / 2)
    end
    return intensity
end

"""
    piecewise_1D_linear(
        z::Array{<:Unitful.Length, 1},
        α::Array{<:PerLength, 1},
        source_function::Array{<:Unitful.Quantity, 1};
        to_end::Bool=false,
        initial_condition=:source,
    )

Compute piecewise integration of the radiative transfer equation,
assuming linear integration of the source function, for a given
height `z`, extinction `α` and `source_function`. The optional
keyword argument `to_end` defines the direction of the integration:
if `false` (default) will start to integrate intensity from the last
element to the first, and if `true` will integrate from the first
element to the last. `initial_condition` can take two values: `:zero` for
no radiation, or `:source` (default) to take the source function at the
starting point.
"""
function piecewise_1D_linear(
    z::Array{<:Unitful.Length, 1},
    α::Array{<:PerLength, 1},
    source_function::Array{<:Unitful.Quantity, 1};
    to_end::Bool=false,
    initial_condition=:source,
)
    ndep = length(z)
    if to_end
        start = 1
        incr = 1
        depth_range = 2:ndep
    else
        start = ndep
        incr = -1
        depth_range = ndep-1:-1:1
    end
    intensity = similar(source_function)
    if initial_condition == :source
        intensity[start] = source_function[start]
    elseif initial_condition == :zero
        intensity[start] *= 0
    else
        error("NotImplemented initial condition $initial_condition")
    end
    for i in depth_range
        Δτ = abs(z[i] - z[i-incr]) * (α[i] + α[i-incr]) / 2
        ΔS = (source_function[i-incr] - source_function[i]) / Δτ
        w1, w2 = _w2(Δτ)
        intensity[i] = (1 - w1)*intensity[i-incr] + w1*source_function[i] + w2*ΔS
    end
    return intensity
end

"""
Computes weights for linear integration of source function,
approximating `exp(-Δτ)` for very small and very large values of `Δτ`.
"""
function _w2(Δτ::T) where T <: AbstractFloat
    if Δτ < 5e-4
        w1 = Δτ * (1 - Δτ / 2)
        w2 = Δτ^2 * (0.5f0 - Δτ / 3)
    elseif Δτ > 50
        w1 = w2 = one(T)
    else
        expΔτ = exp(-Δτ)
        w1 = 1 - expΔτ
        w2 = w1 - Δτ * expΔτ
    end
    return w1, w2
end


text:
const Void = Cvoid
const NULL = C_NULL

const GMT_SESSION_NORMAL   = 0   # Typical mode to GMT_Create_Session
const GMT_SESSION_NOEXIT   = 1   # Call return and not exit when error
const GMT_SESSION_EXTERNAL = 2   # Called by an external API (e.g., Matlab, Julia, Python).
const GMT_SESSION_COLMAJOR = 4   # External API uses column-major formats (e.g., Julai, MATLAB, Fortran). [Row-major format]
const GMT_SESSION_LOGERRORS = 8
const GMT_SESSION_RUNMODE  = 16	 # If set enable GMT's modern runmode. [Classic]
const GMT_SESSION_BEGIN    = 32  # Begin a new session. [Sets modern mode]
const GMT_SESSION_END      = 64  # End a session. [Ends modern mode]
const GMT_SESSION_CLEAR    = 128 # Clear session files/directories
const GMT_SESSION_FIGURE   = 256 # Add a figure to the session queue. [Modern mode only]
# begin enum GMT_enum_runmode
const GMT_CLASSIC        = 0     # Select Classic GMT behavior with -O -K -R -J
const GMT_MODERN         = 1     # Select Modern behavior where -O -K are disabled and -R -J optional if possible
# begin enum GMT_enum_type
const GMT_CHAR = 0
const GMT_UCHAR = 1
const GMT_SHORT = 2
const GMT_USHORT = 3
const GMT_INT = 4
const GMT_UINT = 5
const GMT_LONG = 6
const GMT_ULONG = 7
const GMT_FLOAT = 8
const GMT_DOUBLE = 9
const GMT_TEXT = 10
const GMT_DATETIME = 11
const GMT_N_TYPES = 12
const GMT_VIA_CHAR = 100         # int8_t, 1-byte signed integer type */
const GMT_VIA_UCHAR = 200        # uint8_t, 1-byte unsigned integer type */
const GMT_VIA_SHORT = 300        # int16_t, 2-byte signed integer type */
const GMT_VIA_USHORT = 400       # uint16_t, 2-byte unsigned integer type */
const GMT_VIA_INT = 500          # int32_t, 4-byte signed integer type */
const GMT_VIA_UINT = 600         # uint32_t, 4-byte unsigned integer type */
const GMT_VIA_LONG = 700         # int64_t, 8-byte signed integer type */
const GMT_VIA_ULONG = 800        # uint64_t, 8-byte unsigned integer type */
const GMT_VIA_FLOAT = 900        # 4-byte data float type */
const GMT_VIA_DOUBLE = 1000
# end enum GMT_enum_type
# begin enum GMT_enum_opt
const GMT_OPT_INFILE = 60
# end enum GMT_enum_opt

# begin enum GMT_enum_method
const GMT_IS_OUTPUT = 1024
# end enum GMT_enum_method

# begin enum GMT_enum_via
const GMT_VIA_NONE = 0
const GMT_VIA_MODULE_INPUT = 64
const GMT_VIA_VECTOR = 128
const GMT_VIA_MATRIX = 256
# end enum GMT_enum_via

# GMT_enum_container
const GMT_CONTAINER_AND_DATA = 0    # Create|Read|write both container and the data array
const GMT_CONTAINER_ONLY     = 1    # Create|read|write the container but no data array
const GMT_DATA_ONLY          = 2    # Create|Read|write the container's array only
const GMT_WITH_STRINGS       = 32   # Allocate string array also [DATASET, MATRIX, VECTOR only]
const GMT_NO_STRINGS         = 0

# GMT_enum_read
const GMT_READ_NORMAL = 0	# Normal read mode [Default]
const GMT_READ_DATA   = 1	# Read ASCII data record and return double array
const GMT_READ_TEXT   = 2	# Read ASCII data record and return text string
const GMT_READ_MIXED  = 3   # Read ASCII data record and return double array but tolerate conversion errors
const GMT_READ_FILEBREAK = 4

# begin enum GMT_enum_family
const GMT_IS_DATASET = 0
const GMT_IS_GRID = 1
const GMT_IS_IMAGE = 2
const GMT_IS_CPT = 3				# To be removed whem >= GMT5.3 only
const GMT_IS_PALETTE = 3
const GMT_IS_POSTSCRIPT = 4
const GMT_IS_TEXTSET = 5
const GMT_IS_MATRIX = 6 - (GMTver >= 6)		# FCK, I hate this
const GMT_IS_VECTOR = 7 - (GMTver >= 6)
# begin enum GMT_enum_comment
const GMT_COMMENT_IS_TEXT = 0
# end enum GMT_enum_comment
# begin enum GMT_api_err_enum
const GMT_NOTSET = -1
const GMT_NOERROR = 0
# end enum GMT_api_err_enum
# begin enum GMT_module_enum
const GMT_MODULE_USAGE	  = -6
const GMT_MODULE_SYNOPSIS = -5
const GMT_MODULE_LIST     = -4
const GMT_MODULE_EXIST    = -3
const GMT_MODULE_PURPOSE  = -2
const GMT_MODULE_OPT      = -1
const GMT_MODULE_CMD      = 0
# end enum GMT_module_enum
# begin enum GMT_io_enum
const GMT_IN = 0
const GMT_OUT = 1
const GMT_ERR = 2
# end enum GMT_io_enum
# begin enum GMT_enum_dimensions
const GMT_X = 0
const GMT_Y = 1
const GMT_Z = 2
# end enum GMT_enum_dimensions

const GMT_ALLOC_EXTERNALLY = 0    # Allocated outside of GMT: We cannot reallocate or free this memory
const GMT_ALLOC_INTERNALLY = 1    # Allocated by GMT: We may reallocate as needed and free when no longer needed
# begin enum GMT_enum_write
const GMT_STRICT_CONVERSION = 1024
const GMT_LAX_CONVERSION = 2048
# end enum GMT_enum_write
# begin enum GMT_enum_verbose
const GMT_MSG_QUIET = 0
const GMT_MSG_NORMAL = 1
const GMT_MSG_TICTOC = 2
const GMT_MSG_COMPAT = 3
const GMT_MSG_VERBOSE = 4
const GMT_MSG_LONG_VERBOSE = 5
const GMT_MSG_DEBUG = 6
# end enum GMT_enum_verbose
# begin enum GMT_enum_reg
const GMT_GRID_NODE_REG = 0
const GMT_GRID_PIXEL_REG = 1
# end enum GMT_enum_reg
# begin enum GMT_enum_gridindex
const GMT_XLO = 0
const GMT_XHI = 1
const GMT_YLO = 2
const GMT_YHI = 3
const GMT_ZLO = 4
const GMT_ZHI = 5
# end enum GMT_enum_gridindex
# begin enum GMT_enum_dimindex
const GMT_TBL = 0
const GMT_SEG = 1
const GMT_ROW = 2
const GMT_COL = 3
# end enum GMT_enum_dimindex

# begin enum GMT_enum_gridio
const GMT_GRID_ALL = 0
# end enum GMT_enum_gridio

const DOUBLE_CLASS = 1
const SINGLE_CLASS = 2
const INT64_CLASS  = 3
const UINT64_CLASS = 4
const INT32_CLASS  = 5
const UINT32_CLASS = 6
const INT16_CLASS  = 7
const UINT16_CLASS = 8
const INT8_CLASS   = 9
const UINT8_CLASS  = 10

# begin enum GMT_enum_fmt
if (GMTver < 6.0)
	const GMT_IS_ROW_FORMAT = 0
	const GMT_IS_COL_FORMAT = 1
else
	const GMT_IS_ROW_FORMAT = 1
	const GMT_IS_COL_FORMAT = 2
end
# end enum GMT_enum_fmt

# begin enum GMT_enum_geometry
const GMT_IS_POINT = 1
const GMT_IS_LINE = 2
const GMT_IS_POLY = 4
const GMT_IS_PLP = 7
const GMT_IS_SURFACE = 8
const GMT_IS_NONE = 16
# end enum GMT_enum_geometry

# begin enum GMT_enum_color
const GMT_RGB = 0
const GMT_CMYK = 1
const GMT_HSV = 2
# end enum GMT_enum_color

# GMT_enum_cptflags
const GMT_CPT_NO_BNF     = 1
const GMT_CPT_EXTEND_BNF = 2
const GMT_CPT_HINGED     = 4
const GMT_CPT_TIME       = 8

struct GMT_OPTION			# Structure for a single GMT command option
	option::UInt8				# 1-char command line -<option> (e.g. D in -D) identifying the option (* if file)
	arg::Ptr{UInt8}				# If not NULL, contains the argument for this option
	next::Ptr{GMT_OPTION}
	previous::Ptr{GMT_OPTION}
end

mutable struct GMT_GRID_HEADER_v6
	n_columns::UInt32
	n_rows::UInt32
	registration::UInt32
	wesn::NTuple{4,Cdouble}
	z_min::Cdouble
	z_max::Cdouble
	inc::NTuple{2,Cdouble}
	z_scale_factor::Cdouble
	z_add_offset::Cdouble
	x_unit::NTuple{80,UInt8}
	y_unit::NTuple{80,UInt8}
	z_unit::NTuple{80,UInt8}
	title::NTuple{80,UInt8}
	command::NTuple{320,UInt8}
	remark::NTuple{160,UInt8}
	# Items not stored in the data file for grids but explicitly used in macros computing node numbers
	nm::Csize_t 			# Number of data items in this grid (n_columns * n_rows) [padding is excluded]
	size::Csize_t 			# Actual number of items required to hold this grid (= mx * my)
	bits::UInt32
	complex_mode::UInt32
	_type::UInt32
	n_bands::UInt32
	mx::UInt32 				# Actual dimensions of the grid in memory, allowing for the padding
	my::UInt32
	pad::NTuple{4,UInt32}
	mem_layout::NTuple{4,UInt8}
	nan_value::Cfloat
	xy_off::Cdouble
	ProjRefPROJ4::Ptr{UInt8}
	ProjRefWKT::Ptr{UInt8}
	ProjRefEPSG::Cint
	hidden::Ptr{Cvoid}		# Lower-level information for GMT use only
end

mutable struct GMT_GRID_HEADER_v5
	n_columns::UInt32
	n_rows::UInt32
	registration::UInt32
	wesn::NTuple{4,Cdouble}
	z_min::Cdouble
	z_max::Cdouble
	inc::NTuple{2,Cdouble}
	z_scale_factor::Cdouble
	z_add_offset::Cdouble
	x_unit::NTuple{80,UInt8}
	y_unit::NTuple{80,UInt8}
	z_unit::NTuple{80,UInt8}
	title::NTuple{80,UInt8}
	command::NTuple{320,UInt8}
	remark::NTuple{160,UInt8}
	# Variables "hidden" from the API. This section is flexible and considered private
	_type::UInt32
	bits::UInt32
	complex_mode::UInt32
	mx::UInt32 				# Actual dimensions of the grid in memory, allowing for the padding
	my::UInt32
	nm::Csize_t 			# Number of data items in this grid (n_columns * n_rows) [padding is excluded]
	size::Csize_t 			# Actual number of items required to hold this grid (= mx * my)
	n_alloc::Csize_t
	trendmode::UInt32
	arrangement::UInt32
	n_bands::UInt32
	pad::NTuple{4,UInt32}
	BC::NTuple{4,UInt32}
	grdtype::UInt32
	reset_pad::UInt32
	name::NTuple{256,UInt8}
	varname::NTuple{80,UInt8}
	ProjRefPROJ4::Ptr{UInt8}
	ProjRefWKT::Ptr{UInt8}
	row_order::Cint
	z_id::Cint
	ncid::Cint
	xy_dim::NTuple{2,Cint}
	t_index::NTuple{3,Csize_t}
	data_offset::Csize_t
	stride::UInt32
	nan_value::Cfloat
	xy_off::Cdouble
	r_inc::NTuple{2,Cdouble}
	flags::NTuple{4,UInt8}
	pocket::Ptr{UInt8}
	mem_layout::NTuple{4,UInt8}
	bcr_threshold::Cdouble
	has_NaNs::UInt32
	bcr_interpolant::UInt32
	bcr_n::UInt32
	nxp::UInt32
	nyp::UInt32
	no_BC::UInt32
	gn::UInt32
	gs::UInt32
	is_netcdf4::UInt32
	z_chunksize::NTuple{2,Csize_t}
	z_shuffle::UInt32
	z_deflate_level::UInt32
	z_scale_autoadust::UInt32
	z_offset_autoadust::UInt32
	xy_adjust::NTuple{2,UInt32}
	xy_mode::NTuple{2,UInt32}
	xy_unit::NTuple{2,UInt32}
	xy_unit_to_meter::NTuple{2,Cdouble}
	index_function::Ptr{Cvoid}
end

if (GMTver < 6.0)
	const GMT_GRID_HEADER = GMT_GRID_HEADER_v5
else
	const GMT_GRID_HEADER = GMT_GRID_HEADER_v6
end

mutable struct GMT_GRID_v5
	header::Ptr{GMT_GRID_HEADER}
	data::Ptr{Cfloat}
	id::UInt32
	alloc_level::UInt32
	alloc_mode::UInt32
	x::Ptr{Cdouble}
	y::Ptr{Cdouble}
	extra::Ptr{Cvoid}
end

mutable struct GMT_GRID_v6
	header::Ptr{GMT_GRID_HEADER}
	data::Ptr{Cfloat}
	x::Ptr{Cdouble}
	y::Ptr{Cdouble}
	hidden::Ptr{Cvoid}
end

if (GMTver < 6.0)
	const GMT_GRID = GMT_GRID_v5
else
	const GMT_GRID = GMT_GRID_v6
end

if (GMTver < 6.0)
	struct GMT_OGR
		geometry::UInt32
		n_aspatial::UInt32
		region::Ptr{UInt8}
		proj::NTuple{4,Ptr{UInt8}}
		_type::Ptr{UInt32}
		name::Ptr{Ptr{UInt8}}
		pol_mode::UInt32
		tvalue::Ptr{Ptr{UInt8}}
		dvalue::Ptr{Cdouble}
	end
	struct GMT_OGR_SEG
		pol_mode::UInt32
		n_aspatial::UInt32
		tvalue::Ptr{Ptr{UInt8}}
		dvalue::Ptr{Cdouble}
	end

	mutable struct GMT_DATASEGMENT
		n_rows::UInt64
		n_columns::UInt64
		min::Ptr{Cdouble}
		max::Ptr{Cdouble}
		data::Ptr{Ptr{Cdouble}}
		label::Ptr{UInt8}
		header::Ptr{UInt8}
		text::Ptr{Ptr{UInt8}}
		mode::UInt32
		pol_mode::UInt32
		id::UInt64
		n_alloc::Csize_t
		range::Cint
		pole::Cint
		dist::Cdouble
		lat_limit::Cdouble
		ogr::Ptr{GMT_OGR_SEG}
		next::Ptr{GMT_DATASEGMENT}
		alloc_mode::UInt32
		file::NTuple{2,Ptr{UInt8}}
	end

	mutable struct GMT_DATATABLE
		n_headers::UInt32
		n_columns::UInt64
		n_segments::UInt64
		n_records::UInt64
		min::Ptr{Cdouble}
		max::Ptr{Cdouble}
		header::Ptr{Ptr{UInt8}}
		segment::Ptr{Ptr{GMT_DATASEGMENT}}
		id::UInt64
		n_alloc::Csize_t
		mode::UInt32
		ogr::Ptr{GMT_OGR}
		file::NTuple{2,Ptr{UInt8}}
	end

	mutable struct GMT_DATASET
		n_tables::UInt64
		n_columns::UInt64
		n_segments::UInt64
		n_records::UInt64
		min::Ptr{Cdouble}
		max::Ptr{Cdouble}
		table::Ptr{Ptr{GMT_DATATABLE}}
		id::UInt64
		n_alloc::Csize_t
		dim::NTuple{4,UInt64}
		geometry::UInt32
		alloc_level::UInt32
		io_mode::UInt32
		alloc_mode::UInt32
		file::NTuple{2,Ptr{UInt8}}
	end

else			# GMT6

	mutable struct GMT_DATASEGMENT
		n_rows::UInt64
		n_columns::UInt64
		min::Ptr{Cdouble}
		max::Ptr{Cdouble}
		data::Ptr{Ptr{Cdouble}}
		label::Ptr{UInt8}
		header::Ptr{UInt8}
		text::Ptr{Ptr{UInt8}}
		hidden::Ptr{Cvoid}
	end

	mutable struct GMT_DATATABLE
		n_headers::UInt32
		n_columns::UInt64
		n_segments::UInt64
		n_records::UInt64
		min::Ptr{Cdouble}
		max::Ptr{Cdouble}
		header::Ptr{Ptr{UInt8}}
		segment::Ptr{Ptr{GMT_DATASEGMENT}}
		hidden::Ptr{Cvoid}
	end

	mutable struct GMT_DATASET
		n_tables::UInt64
		n_columns::UInt64
		n_segments::UInt64
		n_records::UInt64
		min::Ptr{Cdouble}
		max::Ptr{Cdouble}
		table::Ptr{Ptr{GMT_DATATABLE}}
		type_::UInt32
		geometry::UInt32
		ProjRefPROJ4::Ptr{UInt8}
		ProjRefWKT::Ptr{UInt8}
		ProjRefEPSG::Cint
		hidden::Ptr{Cvoid}
	end
end

mutable struct GMT_TEXTSEGMENT
	n_rows::UInt64
	data::Ptr{Ptr{UInt8}}
	label::Ptr{UInt8}
	header::Ptr{UInt8}
	id::UInt64
	mode::UInt32
	n_alloc::Csize_t
	file::NTuple{2,Ptr{UInt8}}
	tvalue::Ptr{Ptr{UInt8}}
end
mutable struct GMT_TEXTTABLE
	n_headers::UInt32
	n_segments::UInt64
	n_records::UInt64
	header::Ptr{Ptr{UInt8}}
	segment::Ptr{Ptr{GMT_TEXTSEGMENT}}
	id::UInt64
	n_alloc::Csize_t
	mode::UInt32
	file::NTuple{2,Ptr{UInt8}}
end
mutable struct GMT_TEXTSET
	n_tables::UInt64
	n_segments::UInt64
	n_records::UInt64
	table::Ptr{Ptr{GMT_TEXTTABLE}}
	id::UInt64
	n_alloc::Csize_t
	geometry::UInt32
	alloc_level::UInt32
	io_mode::UInt32
	alloc_mode::UInt32
	file::NTuple{2,Ptr{UInt8}}
end
struct GMT_FILL
	rgb::NTuple{4,Cdouble}
	f_rgb::NTuple{4,Cdouble}
	b_rgb::NTuple{4,Cdouble}
	use_pattern::Bool
	pattern_no::Int32
	dpi::UInt32
	pattern::NTuple{256,UInt8}		# was char pattern[GMT_BUFSIZ];
end
if (GMTver < 6.0)
	struct GMT_LUT
		z_low::Cdouble
		z_high::Cdouble
		i_dz::Cdouble
		rgb_low::NTuple{4,Cdouble}
		rgb_high::NTuple{4,Cdouble}
		rgb_diff::NTuple{4,Cdouble}
		hsv_low::NTuple{4,Cdouble}
		hsv_high::NTuple{4,Cdouble}
		hsv_diff::NTuple{4,Cdouble}
		annot::UInt32
		skip::UInt32
		fill::Ptr{GMT_FILL}
		label::Ptr{UInt8}
	end
else
	struct GMT_LUT
		z_low::Cdouble
		z_high::Cdouble
		i_dz::Cdouble
		rgb_low::NTuple{4,Cdouble}
		rgb_high::NTuple{4,Cdouble}
		rgb_diff::NTuple{4,Cdouble}
		hsv_low::NTuple{4,Cdouble}
		hsv_high::NTuple{4,Cdouble}
		hsv_diff::NTuple{4,Cdouble}
		annot::UInt32
		skip::UInt32
		fill::Ptr{GMT_FILL}
		label::Ptr{UInt8}
		key::Ptr{UInt8}
	end
end
struct GMT_BFN
	rgb::NTuple{4,Cdouble}
	hsv::NTuple{4,Cdouble}
	skip::UInt32
	fill::Ptr{GMT_FILL}
end
mutable struct GMT_PALETTE_v5
	n_headers::UInt32
	n_colors::UInt32
	mode::UInt32
	data::Ptr{GMT_LUT}
	bfn::NTuple{3,GMT_BFN}
	header::Ptr{Ptr{UInt8}}
	id::UInt64
	alloc_mode::UInt32
	alloc_level::UInt32
	auto_scale::UInt32
	model::UInt32
	is_wrapping::UInt32
	is_gray::UInt32
	is_bw::UInt32
	is_continuous::UInt32
	has_pattern::UInt32
	has_hinge::UInt32
	has_range::UInt32
	skip::UInt32
	categorical::UInt32
	z_adjust::NTuple{2,UInt32}
	z_mode::NTuple{2,UInt32}
	z_unit::NTuple{2,UInt32}
	z_unit_to_meter::NTuple{2,Cdouble}
	minmax::NTuple{2,Cdouble}
	hinge::Cdouble
	wrap_length::Cdouble
end
mutable struct GMT_PALETTE_v6
	data::Ptr{GMT_LUT}
	bfn::NTuple{3,GMT_BFN}
	n_headers::UInt32
	n_colors::UInt32
	mode::UInt32
	model::UInt32
	is_wrapping::UInt32
	is_gray::UInt32
	is_bw::UInt32
	is_continuous::UInt32
	has_pattern::UInt32
	has_hinge::UInt32
	has_range::UInt32
	categorical::UInt32
	minmax::NTuple{2,Cdouble}
	hinge::Cdouble
	wrap_length::Cdouble
	header::Ptr{Ptr{UInt8}}
	hidden::Ptr{Cvoid}
end

if (GMTver < 6.0)
	const GMT_PALETTE = GMT_PALETTE_v5
else
	const GMT_PALETTE = GMT_PALETTE_v6
end

mutable struct GMT_IMAGE_v5
	_type::UInt32
	colormap::Ptr{Cint}
	n_indexed_colors::Cint
	header::Ptr{GMT_GRID_HEADER}
	data::Ptr{Cuchar}
	alpha::Ptr{Cuchar}
	id::UInt64
	alloc_level::UInt32
	alloc_mode::UInt32
	color_interp::Ptr{UInt8}
	x::Ptr{Cdouble}
	y::Ptr{Cdouble}
end

mutable struct GMT_IMAGE_v6
	_type::UInt32
	colormap::Ptr{Cint}
	n_indexed_colors::Cint
	header::Ptr{GMT_GRID_HEADER}
	data::Ptr{Cuchar}
	alpha::Ptr{Cuchar}
	color_interp::Ptr{UInt8}
	x::Ptr{Cdouble}
	y::Ptr{Cdouble}
	hidden::Ptr{Cvoid}
end

if (GMTver < 6.0)
	const GMT_IMAGE = GMT_IMAGE_v5
else
	const GMT_IMAGE = GMT_IMAGE_v6
end

mutable struct GMT_POSTSCRIPT_v5
	n_alloc::Csize_t
	n_bytes::Csize_t
	mode::UInt32
	n_headers::UInt32
	data::Ptr{UInt8}
	header::Ptr{Ptr{UInt8}}
	id::UInt64
	alloc_level::UInt32
	alloc_mode::UInt32
end
mutable struct GMT_POSTSCRIPT_v6
	n_bytes::Csize_t
	mode::UInt32
	n_headers::UInt32
	data::Ptr{UInt8}
	header::Ptr{Ptr{UInt8}}
	hidden::Ptr{Cvoid}
end

if (GMTver < 6.0)
	const GMT_POSTSCRIPT = GMT_POSTSCRIPT_v5
else
	const GMT_POSTSCRIPT = GMT_POSTSCRIPT_v6
end

struct GMT_UNIVECTOR
	uc1::Ptr{UInt8}
	sc1::Ptr{Int8}
	ui2::Ptr{UInt16}
	si2::Ptr{Int16}
	ui4::Ptr{UInt32}
	si4::Ptr{Int32}
	ui8::Ptr{UInt64}
	si8::Ptr{Int64}
	f4::Ptr{Float32}
	f8::Ptr{Float64}
end

struct GMT_VECTOR_v6
	n_columns::UInt64
	n_rows::UInt64
	n_headers::UInt32
	registration::UInt32
	_type::Ptr{UInt32}
	range::NTuple{2,Cdouble}
	#data::Ptr{GMT_UNIVECTOR}
	data::Ptr{Ptr{Cvoid}}
	text::Ptr{Ptr{UInt8}};			# Pointer to optional array of strings [NULL]
	command::NTuple{320,UInt8}
	remark::NTuple{160,UInt8}
	header::Ptr{Ptr{UInt8}};		# Array with all Vector header records, if any)
	ProjRefPROJ4::Ptr{UInt8}
	ProjRefWKT::Ptr{UInt8}
	ProjRefEPSG::Cint
	hidden::Ptr{Cvoid}
end

struct GMT_VECTOR_v5
	n_columns::UInt64
	n_rows::UInt64
	registration::UInt32
	_type::Ptr{UInt32}
	#data::Ptr{GMT_UNIVECTOR}
	data::Ptr{Ptr{Cvoid}}
	range::NTuple{2,Cdouble}
	command::NTuple{320,UInt8}
	remark::NTuple{160,UInt8}
	id::UInt64
	alloc_level::UInt32
	alloc_mode::UInt32
end

if (GMTver < 6.0)
	const GMT_VECTOR = GMT_VECTOR_v5
else
	const GMT_VECTOR = GMT_VECTOR_v6
end

mutable struct GMT_MATRIX_v6
	n_rows::UInt64
	n_columns::UInt64
	n_layers::UInt64
	n_headers::UInt32
	shape::UInt32
	registration::UInt32
	dim::Csize_t
	size::Csize_t
	_type::UInt32
	range::NTuple{6,Cdouble}
	inc::NTuple{3,Cdouble}
#	data::GMT_UNIVECTOR
	data::Ptr{Cvoid}
	text::Ptr{Ptr{UInt8}}
	command::NTuple{320,UInt8}
	remark::NTuple{160,UInt8}
	header::Ptr{Ptr{UInt8}};		# Array with all Matrix header records, if any)
	ProjRefPROJ4::Ptr{UInt8}
	ProjRefWKT::Ptr{UInt8}
	ProjRefEPSG::Cint
	hidden::Ptr{Cvoid}
end

mutable struct GMT_MATRIX_v5
	n_rows::UInt64
	n_columns::UInt64
	n_layers::UInt64
	shape::UInt32
	registration::UInt32
	dim::Csize_t
	size::Csize_t
	_type::UInt32
	range::NTuple{6,Cdouble}
#	data::GMT_UNIVECTOR
	data::Ptr{Cvoid}
	command::NTuple{320,UInt8}
	remark::NTuple{160,UInt8}
	id::UInt64
	alloc_level::UInt32
	alloc_mode::UInt32
end

if (GMTver < 6.0)
	const GMT_MATRIX = GMT_MATRIX_v5
else
	const GMT_MATRIX = GMT_MATRIX_v6
end

mutable struct GMT_RESOURCE
	family::UInt32          # GMT data family, i.e., GMT_IS_DATASET, GMT_IS_GRID, etc.
	geometry::UInt32        # One of the recognized GMT geometries
	direction::UInt32       # Either GMT_IN or GMT_OUT
	option::Ptr{GMT_OPTION} # Pointer to the corresponding module option
	name::NTuple{16,UInt8} # Object ID returned by GMT_Register_IO
	pos::Cint               # Corresponding index into external object in|out arrays
	mode::Cint              # Either primary (0) or secondary (1) resource
	object::Ptr{Cvoid}       # Pointer to the actual GMT object
end

struct GMTAPI_DATA_OBJECT
	# Information for each input or output data entity, including information
	# needed while reading/writing from a table (file or array)
	rec::UInt64                 # Current rec to read [GMT_DATASET and GMT_TEXTSET to/from MATRIX/VECTOR only]
	n_rows::UInt64              # Number or rows in this array [GMT_DATASET and GMT_TEXTSET to/from MATRIX/VETOR only]
	n_columns::UInt64			# Number of columns to process in this dataset [GMT_DATASET only]
	n_expected_fields::UInt64	# Number of expected columns for this dataset [GMT_DATASET only]
	n_alloc::Csize_t			# Number of items allocated so far if writing to memory
	ID::UInt32					# Unique identifier which is >= 0
	alloc_level::UInt32			# Nested module level when object was allocated
	status::UInt32				# 0 when first registered, 1 after reading/writing has started, 2 when finished
	selected::Cint				# true if requested by current module, false otherwise
	close_file::Cint			# true if we opened source as a file and thus need to close it when done
	region::Cint				# true if wesn was passed, false otherwise
	no_longer_owner::Cint		# true if the data pointed to by the object was passed on to another object
	messenger::Cint				# true for output objects passed from the outside to receive data from GMT. If true we destroy data pointer before writing
	alloc_mode::UInt32			# GMT_ALLOCATED_{BY_GMT|EXTERNALLY}
	direction::UInt32			# GMT_IN or GMT_OUT
	family::UInt32				# One of GMT_IS_{DATASET|TEXTSET|CPT|IMAGE|GRID|MATRIX|VECTOR|COORD}
	actual_family::UInt32		# May be GMT_IS_MATRIX|VECTOR when one of the others are created via those
	method::UInt32              # One of GMT_IS_{FILE,STREAM,FDESC,DUPLICATE,REFERENCE} or sum with enum GMT_enum_via (GMT_VIA_{NONE,VECTOR,MATRIX,OUTPUT}); using unsigned type because sum exceeds enum GMT_enum_method
	geometry::UInt32			# One of GMT_IS_{POINT|LINE|POLY|PLP|SURFACE|NONE}
	wesn::NTuple{4,Cdouble}	# Grid domain limits
	resource::Ptr{Cvoid}			# Points to registered filename, memory location, etc., where data can be obtained from with GMT_Get_Data.
	data::Ptr{Cvoid}				# Points to GMT object that was read from a resource
	#FILE *fp;					# Pointer to source/destination stream [For rec-by-rec procession, NULL if memory location]
	fp::Ptr{Cvoid}				# Pointer to source/destination stream [For rec-by-rec procession, NULL if memory location]
	filename::Ptr{UInt8}		# Filename, stream, of file handle (otherwise NULL)
	#Cvoid *(*import) (struct GMT_CTRL *, FILE *, uint64_t *, int *);	# Pointer to input function (for DATASET/TEXTSET only)
	ifun::Ptr{Cvoid} 			# Pointer to input function (for DATASET/TEXTSET only)
	# Start of temporary variables for API debug - They are only set when building GMT with /DEBUG
	G::Ptr{Cvoid}				# struct GMT_GRID *G;
	D::Ptr{Cvoid}				# struct GMT_DATASET *D;
	T::Ptr{Cvoid}				# struct GMT_TEXTSET *T;
	C::Ptr{Cvoid}				# struct GMT_PALETTE *C;
	M::Ptr{Cvoid}				# struct GMT_MATRIX *M;
	V::Ptr{Cvoid}				# struct GMT_VECTOR *V;
	I::Ptr{Cvoid}				# struct GMT_IMAGE *I;
end

struct Gmt_libinfo
	name::Ptr{UInt8}	# Library tag name [without leading "lib" and extension], e.g. "gmt", "gmtsuppl" */
	path::Ptr{UInt8}	# Full path to library as given in GMT_CUSTOM_LIBS */
	skip::Ptr{Bool}		# true if we tried to open it and it was not available the first time */
	handle::Ptr{Cvoid}	# Handle to the shared library, returned by dlopen or dlopen_special */
end

#=
mutable struct GMTAPI_CTRL
	# Master controller which holds all GMT API related information at run-time for a single session.
	# Users can run several GMT sessions concurrently; each session requires its own structure.
	# Use GMTAPI_Create_Session to initialize a new session and GMTAPI_Destroy_Session to end it.

	current_rec::NTuple{2,UInt64}	# Current record number >= 0 in the combined virtual dataset (in and out)
	n_objects::UInt32			# Number of currently active input and output data objects
	unique_ID::UInt32			# Used to create unique IDs for duration of session
	session_ID::UInt32			# ID of this session
	unique_var_ID::UInt32		# Used to create unique object IDs (grid,dataset, etc) for duration of session
	current_item::NTuple{2,Cint}	# Array number of current dataset being processed (in and out)
	pad::UInt32					# Session default for number of rows/cols padding for grids [2]
	external::UInt32			# 1 if called via external API (MATLAB, Python) [0]
	runmode::UInt32				# nonzero for GMT modern runmode [0 = classic]
	shape::Cint                 # GMT_IS_COL_FORMAT (1) if column-major (MATLAB, Fortran), GMT_IS_ROW_FORMAT (0) if row-major
	leave_grid_scaled::UInt32	# 1 if we dont want to unpack a grid after we packed it for writing [0]
	n_cores::UInt32             # Number of available cores on this system
	verbose::UInt32             # Used until GMT is set up
	registered::NTuple{2,Bool}	# true if at least one source/destination has been registered (in and out)
	io_enabled::NTuple{2,Bool}	# true if access has been allowed (in and out)
	module_input::Bool          # true when we are about to read inputs to the module (command line) */
	n_objects_alloc::Csize_t	# Allocation counter for data objects
	error::Int32				# Error code from latest API call [GMT_OK]
	last_error::Int32			# Error code from previous API call [GMT_OK]
	shelf::Int32				# Place to pass hidden values within API
	io_mode::NTuple{2,UInt32}	# 1 if access as set, 0 if record-by-record
	PPID::Cint                  # The Process ID of the parent (e.g., shell) or the external caller
	#GMT::Ptr{GMT_CTRL}			# Key structure with low-level GMT internal parameters
	GMT::Ptr{Cvoid}				# Maybe one day. Till than just keep it as Cvoid
	object::Ptr{Ptr{GMTAPI_DATA_OBJECT}}	# List of registered data objects
	session_tag::Ptr{UInt8}		# Name tag for this session (or NULL)
	tmp_dir::Ptr{UInt8}         # System tmp_dir (NULL if not found)
	gwf_dir::Ptr{UInt8}         # GMT WorkFlow dir (NULL if not running in modern mode)
	internal::Bool				# true if session was initiated by gmt.c
	deep_debug::Bool			# temporary for debugging
	#int (*print_func) (FILE *, const char *);	# Pointer to fprintf function (may be reset by external APIs like MEX)
	pf::Ptr{Cvoid}				# Don't know what to put here, so ley it be *void
	do_not_exit::UInt32			# 0 by default, mieaning it is OK to call exit  (may be reset by external APIs like MEX to call return instead)
	lib::Ptr{Gmt_libinfo}		# List of shared libs to consider
	n_shared_libs::UInt32		# How many in lib
end
=#


text:
using CoarseGrainingPrism3D
println("data")
@time data = PrAOpti()
@show length(data[1])

println("prism eff")
#@time dataEff = PrismEff(data,1)
#@show length(dataEff[1])

println("test 2-2 move on initial prism")
println("one time")
# @time Move2_2(data,12,12,[10,7,3,11])
# @time dataIniPrismMove2_2 = Move2_2(data,12,12,[10,7,3,11])
# @show length(dataIniPrismMove2_2[2])
# @show sum(numchop(sort(data[1])-sort(dataIniPrismMove2_2[1])).==0)

println("3-1 move test")
#Move3_1testIniPrismtoPyra = Move3_1(data,12,[1,4,8,9,5,11])
#ExtractPyraFromSVD = FullsvdPrismTruncated(data,1)
#StoreExtractPyraFromSVDAmp = Array{Real}(undef,0)
#StoreExtractPyraFromSVDIndex = Array{Int64}(undef,0)
#for i in 1:length(ExtractPyraFromSVD[2])
#    for j in 1:length(ExtractPyraFromSVD[2][i])
#        push!(StoreExtractPyraFromSVDAmp,ExtractPyraFromSVD[2][i][j])
        #push!(StoreExtractPyraFromSVDIndex,ExtractPyraFromSVD[4][i][j])
#    end
#end
#@show (length(Move3_1testIniPrismtoPyra[1]),length(Move3_1testIniPrismtoPyra[2]))
#@show (length(StoreExtractPyraFromSVDAmp),length(StoreExtractPyraFromSVDIndex))

#@show sum(numchop(sort(Move3_1testIniPrismtoPyra[1])-sort(StoreExtractPyraFromSVDAmp)).==0)
#@show sum(numchop(sort(Move3_1testIniPrismtoPyra[2])-sort(StoreExtractPyraFromSVDIndex)).==0)
#
#println("time test for move")
#println("one order")
#@time Move2_2FirstTriangle = Move2_2(dataEff,18,16,[2,10,6,13])
#@time Move2_2_BigDiag = Move2_2(Move2_2FirstTriangle,18,18,[10,11,14,15])
#@time Move3_1FirstTriangle = Move3_1(Move2_2_BigDiag,18,[10,13,14,18,16,7])
#@show length(Move3_1FirstTriangle[2])
#@time Move2_2SecondTriangle = Move2_2(Move3_1FirstTriangle,15,14,[4,11,8,12])
#@time Move3_1SecondTriangle = Move3_1(Move2_2SecondTriangle,15,[10,11,12,15,3,14])
#@show length(Move3_1SecondTriangle[2])
println("another order")
#@time Move2_2_BigDiag = Move2_2(dataEff,18,18,[10,11,14,15])
#@time Move2_2FirstTriangle = Move2_2(Move2_2_BigDiag,18,16,[2,10,6,13])
#@time Move3_1FirstTriangle = Move3_1(Move2_2FirstTriangle,18,[10,13,14,18,16,7])
#@show length(Move3_1FirstTriangle[2])
#@time Move2_2SecondTriangle = Move2_2(Move3_1FirstTriangle,15,14,[4,11,8,12])
#@time Move3_1SecondTriangle = Move3_1(Move2_2SecondTriangle,15,[10,11,12,15,3,14])
#@show length(Move3_1SecondTriangle[2])

#
println("test Embedding")
@time DataEmbPrismEff = EmbPrismEff(data,1)
# @show length(unique(DataEmbPrismEff[2]))
# @show length(DataEmbPrismEff[1])
# @show sum(DataEmbPrismEff[1].!=0)

println("comparing after embedding and initial prism")
# @show sum(numchop(sort(DataEmbPrismEff[1])-sort(data[1])).==0)
# @show sum((sort(DataEmbPrismEff[2])-sort(data[2])).==0)


text:
# ---
# title: Fourier transform of Moisan periodic image component
# id: demo_perfft2
# cover: assets/perfft2.png
# author: Peter Kovesi
# date: 2018-10-26
# ---

# The function `perfft2()` implements Moisan's "Periodic plus Smooth Image
# Decomposition" which decomposes an image into two components
#
#         img = p + s
#
# where `s` is the 'smooth' component with mean 0 and `p` is the 'periodic' component
# which has no sharp discontinuities when one moves cyclically across the image
# boundaries.
#
# This decomposition is very useful when one wants to obtain an FFT of an image
# with minimal artifacts introduced from the boundary discontinuities. The image
# `p` gathers most of the image information but avoids periodization artifacts.
#
# Reference:
# L. Moisan, "Periodic plus Smooth Image Decomposition", Journal of
# Mathematical Imaging and Vision, vol 39:2, pp. 161-179, 2011.

using Images
using FFTW
using ImagePhaseCongruency
using ImageContrastAdjustment
using TestImages

img = Float64.(Gray.(testimage("lena")))

IMG = fft(img)               # 'Standard' fft
(P, S, p, s) = perfft2(img)  # 'Periodic' fft

mosaic(
    adjust_histogram(Gray.(p), LinearStretching()),
    adjust_histogram(s, LinearStretching()),
    ## Note the vertical and horizontal cross in
    ## the spectrum induced by the non-periodic edges.
    adjust_histogram(log.(abs.(fftshift(IMG)) .+ 1), LinearStretching()),
    ## Note the clean spectrum because p is periodic.
    adjust_histogram(log.(abs.(fftshift(P)) .+ 1), LinearStretching());
    nrow=2, rowmajor=true
)
# Top 1) left: periodic component 2) right: smooth component
#
# Bottom 3) left: spectrum of standard FFT 4) right: spectrum of periodic component

# save cover image #src
isdir("assets") || mkdir("assets") #src
cover = Gray.(adjust_histogram(log.(abs.(fftshift(P)) .+ 1), LinearStretching())) #src
save(joinpath("assets", "perfft2.png"), cover) #src

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
module MCMC

using Distributions
using StatsFuns
import Dates
using ProgressBars
import LinearAlgebra

include("misc.jl")
include("TuningParam.jl")
include("MvTuner.jl")
include("metropolis.jl")
include("gibbs.jl")

end # module

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    API implementation for J2 osculating orbit propagator.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

get_epoch(orbp::OrbitPropagatorJ2osc) = orbp.j2d.epoch

function get_mean_elements(orbp::OrbitPropagatorJ2osc)
    j2d = orbp.j2oscd.j2d
    orb = KeplerianElements(j2d.epoch + j2d.Δt/86400,
                            j2d.al_k * j2d.j2_gc.R0,
                            j2d.e_k,
                            j2d.i_k,
                            j2d.Ω_k,
                            j2d.ω_k,
                            j2d.f_k)

    return orb
end

"""
    init_orbit_propagator(Val(:J2osc), epoch::Number, a_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, f_0::Number, dn_o2::Number = 0, ddn_o6::Number = 0; j2_gc::J2_GravCte{T} = j2_gc_egm08) where T
    init_orbit_propagator(Val(:J2osc), orb_0::Orbit, dn_o2::Number = 0, ddn_o6::Number = 0; j2_gc::J2_GravCte = j2_gc_egm08)

Initialize the J2 osculating orbit propagator.

# Args

* `epoch`: Initial orbit epoch [Julian Day].
* `a_0`: Initial mean semi-major axis [m].
* `e_0`: Initial mean eccentricity.
* `i_0`: Initial mean inclination [rad].
* `Ω_0`: Initial mean right ascension of the ascending node [rad].
* `ω_0`: Initial mean argument of perigee [rad].
* `f_0`: Initial mean true anomaly [rad].
* `dn_o2`: (OPTIONAL) First time derivative of mean motion divided by 2
           \\[rad/s²] (**Default** = 0).
* `ddn_o6`: (OPTIONAL) Second time derivative of mean motion divided by 6
            \\[rad/s³] (**Default** = 0).
* `orb_0`: Instance of the structure `KeplerianElements` with the initial mean
           orbital elements [SI].

## Keywords

* `j2_gc`: (OPTIONAL) J2 orbit propagator gravitational constants
           (**Default** = `j2_gc_egm08`).

"""
function init_orbit_propagator(::Val{:J2osc}, epoch::Number, a_0::Number,
                               e_0::Number, i_0::Number, Ω_0::Number,
                               ω_0::Number, f_0::Number, dn_o2::Number = 0,
                               ddn_o6::Number = 0;
                               j2_gc::J2_GravCte{T} = j2_gc_egm08) where T

    # Create the new Two Body propagator structure.
    j2oscd = j2osc_init(epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0, dn_o2, ddn_o6;
                        j2_gc = j2_gc)

    # Create and return the orbit propagator structure.
    return OrbitPropagatorJ2osc(j2oscd)
end

function init_orbit_propagator(::Val{:J2osc},
                               orb_0::Orbit,
                               dn_o2::Number = 0,
                               ddn_o6::Number = 0;
                               j2_gc::J2_GravCte = j2_gc_egm08)
    # Convert the orbit representation to Keplerian elements.
    k_0 = convert(KeplerianElements, orb_0)

    return init_orbit_propagator(Val(:J2osc), k_0.t, k_0.a, k_0.e, k_0.i, k_0.Ω,
                                 k_0.ω, k_0.f, dn_o2, ddn_o6; j2_gc = j2_gc)
end

function propagate!(orbp::OrbitPropagatorJ2osc{T}, t::Number) where T
    # Auxiliary variables.
    j2oscd = orbp.j2oscd

    # Propagate the orbit.
    r_i, v_i = j2osc!(j2oscd, t)

    # Return.
    return r_i, v_i
end

function step!(orbp::OrbitPropagatorJ2osc, Δt::Number)
    # Auxiliary variables.
    j2oscd = orbp.j2oscd

    # Propagate the orbit.
    r_i, v_i = j2osc!(j2oscd, j2oscd.Δt + Δt)

    # Return the information about the step.
    return r_i, v_i
end


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

"""
Unit tests for the package installer/uninstaller functions in `TestTools`.
"""

# --- Imports

# Standard library
using Test

# External packages
using Suppressor

# Local modules
using TestTools: TestTools
using TestTools.jltest: EnhancedTestSet

# --- Tests

@testset EnhancedTestSet "TestTools: install(), uninstall()" begin

    # --- Preparations

    # Construct path to test installation directory
    bin_dir = abspath(joinpath(@__DIR__, "testing-bin-dir"))

    # Cache common variables
    jltest_cmd = "jltest"
    jlcoverage_cmd = "jlcoverage"
    jlcodestyle_cmd = "jlcodestyle"
    if Sys.iswindows()
        jltest_cmd = "$jltest_cmd.cmd"
        jlcoverage_cmd = "$jlcoverage_cmd.cmd"
        jlcodestyle_cmd = "$jlcodestyle_cmd.cmd"
    end

    # --- install() tests

    jltest_exec_path = Base.contractuser(joinpath(bin_dir, jltest_cmd))
    jlcoverage_exec_path = Base.contractuser(joinpath(bin_dir, jlcoverage_cmd))
    jlcodestyle_exec_path = Base.contractuser(joinpath(bin_dir, jlcodestyle_cmd))

    expected_output_install = """
        [ Info: Installed $jltest_cmd to `$jltest_exec_path`.
        [ Info: Installed $jlcoverage_cmd to `$jlcoverage_exec_path`.
        [ Info: Installed $jlcodestyle_cmd to `$jlcodestyle_exec_path`.
        ┌ Info: Make sure that `$bin_dir` is in PATH, or manually add a
        └ symlink from a directory in PATH to the installed program file.
        """

    # Remove existing install directory
    if isdir(bin_dir)
        rm(bin_dir; force=true, recursive=true)
    end

    # Case: fresh installation
    output = @capture_err begin
        TestTools.install(; bin_dir=bin_dir)
    end

    @test output == expected_output_install

    # Case: attempt to reinstall with force=false
    local error
    try
        TestTools.install(; bin_dir=bin_dir)
    catch error
    end

    expected_error =
        "File `$jltest_exec_path` already exists. " *
        "Use `TestTools.install(force=true)` to overwrite."
    @test error.msg == expected_error

    # Case: attempt to reinstall with force=true
    output = @capture_err begin
        TestTools.install(; bin_dir=bin_dir, force=true)
    end

    @test output == expected_output_install

    # --- uninstall() tests

    # Case: uninstall
    output = @capture_err begin
        TestTools.uninstall(; bin_dir=bin_dir)
    end

    expected_output = """
        [ Info: Uninstalled `$jltest_exec_path`.
        [ Info: Uninstalled `$jlcoverage_exec_path`.
        [ Info: Uninstalled `$jlcodestyle_exec_path`.
        """
    @test output == expected_output

    # --- Clean up

    rm(bin_dir)
end

@testset EnhancedTestSet "TestTools: install_cli(), uninstall_cli(): invalid arguments" begin

    # --- install_cli() tests

    # Case: invalid `cli` arg
    local error
    try
        TestTools.install_cli("invalid-cli")
    catch error
    end

    expected_error = "Invalid `cli`: invalid-cli"
    @test error.msg == expected_error

    # --- uninstall_cli() tests

    # Case: invalid `cli` arg
    local error
    try
        TestTools.uninstall_cli("invalid-cli")
    catch error
    end

    expected_error = "Invalid `cli`: invalid-cli"
    @test error.msg == expected_error
end


text:
using SpecialFunctions
using Statistics
using FFTW

function fMultiSinGen(N::Integer,
                      P::Integer,
                      M::Integer;
                      fMin::Float64=0.99,
                      fMax::Float64=1.00,
                      fs::Float64=1.0,
                      type_signal::String="full",
                      nGroup::Integer=3,
                      uStd::Float64=1.)
    """
    generates a zero-mean random phase multisine with std = 1
    INPUT
    options.N: number of points per period
    options.P: number of periods
    options.M: number of realizations
    options.fMin: minimum excited frequency
    options.fMax: maximum escited frequency
    options.fs: sample frequency
    options.type: "full", "odd", "oddrandom"

    OPTIONAL
    options.nGroup: in case of oddrandom, 1 out of nGroup odd lines is
                     discarded. Default = 3
    options.std: std of the generated signals. Default = 1

    OUTPUT
    u: NPxM record of the generated signals
    lines: excited frequency lines -> 1 = dc, 2 = fs/N

    copyright:
    Maarten Schoukens
    Vrije Universiteit Brussel, Brussels Belgium
    10/05/2017

    translated from Matlab to Julia by
    Wouter Kouw
    TU Eindhoven, Eindhoven, Netherlands
    22/01/2021

    This work is licensed under a
    Creative Commons Attribution-NonCommercial 4.0 International License
    (CC BY-NC 4.0)
    https://creativecommons.org/licenses/by-nc/4.0/
    """

    # Lines selection - select which frequencies to excite
    f0 = fs/N
    linesMin = Int64(ceil(fMin / f0) + 1)
    linesMax = Int64(floor(fMax / f0) + 1)
    lines = linesMin:linesMax

    # Remove DC component
    if lines[1] == 1; lines = lines[2:end]; end

    if type_signal == "full"
        # do nothing
    elseif type_signal == "odd"

        # remove even lines - odd indices
        if Bool(mod(lines[1],2)) # lines(1) is odd
            lines = lines[2:2:end]
        else
            lines = lines[1:2:end]
        end

    elseif type_signal == "oddrandom"

        # remove even lines - odd indices
        if Bool(mod(lines[1],2)) # lines(1) is odd
            lines = lines[2:2:end]
        else
            lines = lines[1:2:end]
        end

        # remove 1 out of nGroup lines
        nLines = length(lines)
        nRemove = floor(nLines / nGroup)
        removeInd = rand(1:nGroup, [1 nRemove])
        removeInd = removeInd + nGroup*[0:nRemove-1]
        lines = lines(!removeInd)
    end
    nLines = length(lines)

    # multisine generation - frequency domain implementation
    U = zeros(ComplexF64, N,M)

    # excite the selected frequencies
    U[lines,:] = exp.(2im*pi*rand(nLines,M))

    # go to time domain
    u = real(ifft(U))

    # rescale to obtain desired rms std
    u = uStd * u ./ std(u[:,1])

    # generate P periods
    u = repeat(u, outer=(P,1))

    return u, lines
end


text:
struct TransformsRow{T, F}
    row::T
    funcs::F
end

Base.getproperty(row::TransformsRow, nm::Symbol) = (getfunc(row, getfield(row, 2), nm))(getproperty(getfield(row, 1), nm))
Base.propertynames(row::TransformsRow) = propertynames(getfield(row, 1))

struct Transforms{C, T, F}
    source::T
    funcs::F # NamedTuple of columnname=>transform function
end

Base.propertynames(t::Transforms{true}) = propertynames(getfield(t, 1))
Base.getproperty(t::Transforms{true}, nm::Symbol) = Base.map(getfunc(t, getfield(t, 2), nm), getproperty(getfield(t, 1), nm))

transform(funcs) = x->transform(x, funcs)
transform(; kw...) = transform(kw.data)
function transform(src::T, funcs::F) where {T, F}
    C = columnaccess(T)
    x = C ? columns(src) : rows(src)
    return Transforms{C, typeof(x), F}(x, funcs)
end

getfunc(row, nt::NamedTuple, nm) = get(nt, nm, identity)
getfunc(row, d::Dict{String, <:Base.Callable}, nm) = get(d, String(nm), identity)
getfunc(row, d::Dict{Symbol, <:Base.Callable}, nm) = get(d, nm, identity)
getfunc(row, d::Dict{Int, <:Base.Callable}, nm) = get(d, findfirst(isequal(nm), propertynames(row)), identity)

istable(::Type{<:Transforms}) = true
rowaccess(::Type{Transforms{C, T, F}}) where {C, T, F} = !C
rows(t::Transforms{false}) = t
columnaccess(::Type{Transforms{C, T, F}}) where {C, T, F} = C
columns(t::Transforms{true}) = t
# avoid relying on inference here and just let sinks figure things out
schema(t::Transforms) = nothing
IteratorInterfaceExtensions.isiterable(x::Transforms) = true
IteratorInterfaceExtensions.getiterator(x::Transforms) = datavaluerows(columntable(x))

Base.IteratorSize(::Type{Transforms{false, T, F}}) where {T, F} = Base.IteratorSize(T)
Base.length(t::Transforms{false}) = length(getfield(t, 1))
Base.eltype(t::Transforms{false, T, F}) where {T, F} = TransformsRow{eltype(getfield(t, 1)), F}

@inline function Base.iterate(t::Transforms{false}, st=())
    state = iterate(getfield(t, 1), st...)
    state === nothing && return nothing
    return TransformsRow(state[1], getfield(t, 2)), (state[2],)
end

# select
struct Select{T, columnaccess, names}
    source::T
end

select(names::Symbol...) = x->select(x, names...)
select(names::String...) = x->select(x, Base.map(Symbol, names)...)
select(inds::Integer...) = x->select(x, Base.map(Int, inds)...)

function select(x::T, names...) where {T}
    colaccess = columnaccess(T)
    r = colaccess ? columns(x) : rows(x)
    return Select{typeof(r), colaccess, names}(r)
end

istable(::Type{<:Select}) = true
IteratorInterfaceExtensions.isiterable(x::Select) = true
IteratorInterfaceExtensions.getiterator(x::Select) = datavaluerows(x)

Base.@pure function typesubset(::Schema{names, types}, nms::NTuple{N, Symbol}) where {names, types, N}
    return Tuple{Any[columntype(names, types, nm) for nm in nms]...}
end

Base.@pure function typesubset(::Schema{names, types}, inds::NTuple{N, Int}) where {names, types, N}
    return Tuple{Any[fieldtype(types, i) for i in inds]...}
end

typesubset(::Schema{names, types}, ::Tuple{}) where {names, types} = Tuple{}

namesubset(::Schema{names, types}, nms::NTuple{N, Symbol}) where {names, types, N} = nms
Base.@pure namesubset(::Schema{names, T}, inds::NTuple{N, Int}) where {names, T, N} = ntuple(i -> names[inds[i]], N)
namesubset(::Schema{names, types}, ::Tuple{}) where {names, types} = ()
namesubset(names, nms::NTuple{N, Symbol}) where {N} = nms
namesubset(names, inds::NTuple{N, Int}) where {N} = ntuple(i -> names[inds[i]], N)
namesubset(names, ::Tuple{}) = ()

function schema(s::Select{T, columnaccess, names}) where {T, columnaccess, names}
    sch = schema(getfield(s, 1))
    sch === nothing && return nothing
    return Schema(namesubset(sch, names), typesubset(sch, names))
end

# columns: make Select property-accessible
Base.getproperty(s::Select{T, true, names}, nm::Symbol) where {T, names} = getproperty(getfield(s, 1), nm)
Base.propertynames(s::Select{T, true, names}) where {T, names} = namesubset(propertynames(getfield(s, 1)), names)
columnaccess(::Type{Select{T, C, names}}) where {T, C, names} = C
columns(s::Select{T, true, names}) where {T, names} = s

# rows: implement Iterator interface
Base.IteratorSize(::Type{Select{T, false, names}}) where {T, names} = Base.IteratorSize(T)
Base.length(s::Select{T, false, names}) where {T, names} = length(getfield(s, 1))
Base.IteratorEltype(::Type{Select{T, false, names}}) where {T, names} = Base.IteratorEltype(T)
Base.eltype(s::Select{T, false, names}) where {T, names} = SelectRow{eltype(getfield(s, 1)), names}
rowaccess(::Type{Select{T, columnaccess, names}}) where {T, columnaccess, names} = !columnaccess
rows(s::Select{T, false, names}) where {T, names} = s

# we need to iterate a "row view" in case the underlying source has unknown schema
# to ensure each iterated row only has `names` propertynames
struct SelectRow{T, names}
    row::T
end

Base.getproperty(row::SelectRow, nm::Symbol) = getproperty(getfield(row, 1), nm)

getprops(row, nms::NTuple{N, Symbol}) where {N} = nms
getprops(row, inds::NTuple{N, Int}) where {N} = ntuple(i->propertynames(getfield(row, 1))[inds[i]], N)
getprops(row, ::Tuple{}) = ()

Base.propertynames(row::SelectRow{T, names}) where {T, names} = getprops(row, names)

@inline function Base.iterate(s::Select{T, false, names}) where {T, names}
    state = iterate(getfield(s, 1))
    state === nothing && return nothing
    row, st = state
    return SelectRow{typeof(row), names}(row), st
end

@inline function Base.iterate(s::Select{T, false, names}, st) where {T, names}
    state = iterate(getfield(s, 1), st)
    state === nothing && return nothing
    row, st = state
    return SelectRow{typeof(row), names}(row), st
end


text:
using Compose, Plots, Printf
import Cairo, Fontconfig
using Graphs: AbstractGraph, Edge
using GraphPlot: spring_layout, gplot
using Colors: Colorant, FixedPointNumbers, RGB

struct FixedGraphPlot{T<:Integer}
    # fixed graph plot to change colors and edgelables
    graph::AbstractGraph{T}
    locs_x::Vector{Float64}
    locs_y::Vector{Float64}
    NODESIZE::Float64
    nodelabels::Union{Vector{T},Vector{String}}
    nodecolors::Vector{RGB{FixedPointNumbers.N0f8}}
    edgeindices::Dict{Edge,T} # not compatbile with flip_edge! - call reset_edgeindices after an edge has changed
    edgelabels::Vector{Char}
    edgecolors::Vector{RGB{FixedPointNumbers.N0f8}}

    function FixedGraphPlot{T}(graph::AbstractGraph{T}; layout = spring_layout(graph), NODESIZE = 0.25 / sqrt(nv(graph)), nodelabels = 1:nv(graph)) where {T<:Integer}
        locs_x, locs_y = layout
        nodelabels = collect(nodelabels)
        nodecolors = fill(colorant"white", nv(graph))
        edgeindices = Dict{Edge,T}()
        for (i, e) in enumerate(edges(graph))
            edgeindices[e] = i
            edgeindices[reverse(e)] = i
        end
        edgelabels = fill(' ', ne(graph))
        edgecolors = fill(colorant"purple", ne(graph))
        new(graph, locs_x, locs_y, NODESIZE, nodelabels, nodecolors, edgeindices, edgelabels, edgecolors)
    end
end

function reset_edgeindices!(fgp::FixedGraphPlot, clear = false)
    if clear
        empty!(fgp.edgeindices)
    end
    for (i, e) in enumerate(edges(fgp.graph))
        fgp.edgeindices[e] = i
        fgp.edgeindices[reverse(e)] = i
    end
    return fgp
end

function frame(anim::Animation, gp::Context, width = 20cm, height = 12cm; dpi = 96)
    i = length(anim.frames) + 1
    filename = @sprintf("%06d.png", i)
    draw(PNG(joinpath(anim.dir, filename), width, height, dpi = dpi), gp)
    push!(anim.frames, filename)
end

function graph_context(fgp::FixedGraphPlot)
    return gplot(fgp.graph, fgp.locs_x, fgp.locs_y, NODESIZE = fgp.NODESIZE, nodelabel = fgp.nodelabels, nodefillc = fgp.nodecolors, edgelabel = fgp.edgelabels, edgestrokec = fgp.edgecolors)
end

max_stars_repo_name: tomtuamnuq/dist-algorithms

max_issues_repo_name: tomtuamnuq/dist-algorithms

max_forks_repo_name: tomtuamnuq/dist-algorithms


text:
@testset "Ising Model" begin
    IsingSpin = MonteCarlo.IsingSpin

    # constructors
    m = IsingModel(L=8,dims=2);
    @test m.L == 8 && m.dims == 2
    @test ndims(m) == 2
    @test MonteCarlo.nsites(m) == 64
    @test typeof(m) == MonteCarlo.IsingModel{MonteCarlo.SquareLattice}
    m = IsingModel(dims=1, L=10);
    @test typeof(m) == MonteCarlo.IsingModel{MonteCarlo.Chain}
    @test m.L == 10 && m.dims == 1
    @test MonteCarlo.nsites(m) == 10
    @test ndims(m) == 1
    d = Dict(:dims=>3, :L=>3)
    m = IsingModel(d)
    @test typeof(m) == MonteCarlo.IsingModel{MonteCarlo.CubicLattice{Array{Int64,3}}}
    @test m.L == 3 && m.dims == 3

    # energy, general
    @test MonteCarlo.energy(MC(m, beta=1), m, reshape(IsingSpin.(1:27), (3,3,3))) == -18333.0
    # energy, square lattice
    m = IsingModel(dims=2, L=8)
    mc = MC(m, beta=1)
    conf = reshape(IsingSpin.(1:64), (8,8))
    conff = deepcopy(conf)
    @test MonteCarlo.energy(mc, m, conf) == -164320.0

    # rand, conftype
    Random.seed!(123)
    @test MonteCarlo.rand(MC, m) == IsingSpin[-1 -1 1 -1 1 1 1 -1; -1 1 1 -1 -1 1 1 1; -1 1 1 1 1 1 -1 1; -1 1 1 1 -1 -1 1 1; 1 1 1 -1 1 -1 -1 -1; 1 -1 -1 -1 1 1 -1 1; -1 -1 1 1 1 1 -1 1; 1 -1 1 -1 -1 -1 1 1]

    # propose, accept
    @test MonteCarlo.propose_local(mc, m, 13, conff) == (1352.0, nothing)
    @test conf == conff
    @test MonteCarlo.accept_local!(mc, m, 13, conff, 2, 1352.0) == nothing
    conff[13] *= -1
    @test conf == conff

    # observables
    obs = MonteCarlo.prepare_observables(mc, m);
    @test typeof(obs) == Dict{String,MonteCarloObservable.Observable}
    @test MonteCarlo.measure_observables!(mc, m, obs, conf) == nothing
    @test MonteCarlo.finish_observables!(mc, m, obs) == nothing
end

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC