text:
import GPLinearODEMaker.powers_of_negative_one

"""
se_kernel function created by kernel_coder(). Requires 1 hyperparameters. Likely created using se_kernel_base() as an input.
Use with include("src/kernels/se_kernel.jl").
hyperparameters == ["λ"]
"""
function se_kernel(
    hyperparameters::Vector{<:Real},
    δ::Real,
    dorder::Vector{<:Integer};
    shift_ind::Integer=0)

    include_shift = shift_ind!=0
    @assert length(hyperparameters)==1+Int(include_shift) "hyperparameters is the wrong length"
    dorder_len = 3 + Int(include_shift)
    @assert length(dorder)==dorder_len "dorder is the wrong length"
    @assert maximum(dorder) < 3 "No more than two derivatives for each time or hyperparameter can be calculated"

    dorder2 = dorder[2]
    dorder[2] += dorder[1]

    if include_shift
        dorder[2] += dorder[shift_ind+2]
        δ += hyperparameters[shift_ind]
        hyperparameters_view = view(hyperparameters, 1:2 .!= shift_ind)
        dorder_view = view(dorder, 2:dorder_len)
        dorder_view = view(dorder_view, 1:(dorder_len-1) .!= (shift_ind+1))
    else
        hyperparameters_view = hyperparameters
        dorder_view = view(dorder, 2:dorder_len)
    end

    λ = hyperparameters_view[1]

    if dorder_view==[6, 2]
        func = -630*exp((-1/2)*δ^2/λ^2)/λ^8 + 3465*exp((-1/2)*δ^2/λ^2)*δ^2/λ^10 - 2520*exp((-1/2)*δ^2/λ^2)*δ^4/λ^12 + 546*exp((-1/2)*δ^2/λ^2)*δ^6/λ^14 - 42*exp((-1/2)*δ^2/λ^2)*δ^8/λ^16 + exp((-1/2)*δ^2/λ^2)*δ^10/λ^18
    end

    if dorder_view==[5, 2]
        func = -630*exp((-1/2)*δ^2/λ^2)*δ/λ^8 + 945*exp((-1/2)*δ^2/λ^2)*δ^3/λ^10 - 315*exp((-1/2)*δ^2/λ^2)*δ^5/λ^12 + 33*exp((-1/2)*δ^2/λ^2)*δ^7/λ^14 - exp((-1/2)*δ^2/λ^2)*δ^9/λ^16
    end

    if dorder_view==[4, 2]
        func = 60*exp((-1/2)*δ^2/λ^2)/λ^6 - 285*exp((-1/2)*δ^2/λ^2)*δ^2/λ^8 + 165*exp((-1/2)*δ^2/λ^2)*δ^4/λ^10 - 25*exp((-1/2)*δ^2/λ^2)*δ^6/λ^12 + exp((-1/2)*δ^2/λ^2)*δ^8/λ^14
    end

    if dorder_view==[3, 2]
        func = 60*exp((-1/2)*δ^2/λ^2)*δ/λ^6 - 75*exp((-1/2)*δ^2/λ^2)*δ^3/λ^8 + 18*exp((-1/2)*δ^2/λ^2)*δ^5/λ^10 - exp((-1/2)*δ^2/λ^2)*δ^7/λ^12
    end

    if dorder_view==[2, 2]
        func = -6*exp((-1/2)*δ^2/λ^2)/λ^4 + 27*exp((-1/2)*δ^2/λ^2)*δ^2/λ^6 - 12*exp((-1/2)*δ^2/λ^2)*δ^4/λ^8 + exp((-1/2)*δ^2/λ^2)*δ^6/λ^10
    end

    if dorder_view==[1, 2]
        func = -6*exp((-1/2)*δ^2/λ^2)*δ/λ^4 + 7*exp((-1/2)*δ^2/λ^2)*δ^3/λ^6 - exp((-1/2)*δ^2/λ^2)*δ^5/λ^8
    end

    if dorder_view==[0, 2]
        func = -3*exp((-1/2)*δ^2/λ^2)*δ^2/λ^4 + exp((-1/2)*δ^2/λ^2)*δ^4/λ^6
    end

    if dorder_view==[6, 1]
        func = 90*exp((-1/2)*δ^2/λ^2)/λ^7 - 375*exp((-1/2)*δ^2/λ^2)*δ^2/λ^9 + 195*exp((-1/2)*δ^2/λ^2)*δ^4/λ^11 - 27*exp((-1/2)*δ^2/λ^2)*δ^6/λ^13 + exp((-1/2)*δ^2/λ^2)*δ^8/λ^15
    end

    if dorder_view==[5, 1]
        func = 90*exp((-1/2)*δ^2/λ^2)*δ/λ^7 - 95*exp((-1/2)*δ^2/λ^2)*δ^3/λ^9 + 20*exp((-1/2)*δ^2/λ^2)*δ^5/λ^11 - exp((-1/2)*δ^2/λ^2)*δ^7/λ^13
    end

    if dorder_view==[4, 1]
        func = -12*exp((-1/2)*δ^2/λ^2)/λ^5 + 39*exp((-1/2)*δ^2/λ^2)*δ^2/λ^7 - 14*exp((-1/2)*δ^2/λ^2)*δ^4/λ^9 + exp((-1/2)*δ^2/λ^2)*δ^6/λ^11
    end

    if dorder_view==[3, 1]
        func = -12*exp((-1/2)*δ^2/λ^2)*δ/λ^5 + 9*exp((-1/2)*δ^2/λ^2)*δ^3/λ^7 - exp((-1/2)*δ^2/λ^2)*δ^5/λ^9
    end

    if dorder_view==[2, 1]
        func = 2*exp((-1/2)*δ^2/λ^2)/λ^3 - 5*exp((-1/2)*δ^2/λ^2)*δ^2/λ^5 + exp((-1/2)*δ^2/λ^2)*δ^4/λ^7
    end

    if dorder_view==[1, 1]
        func = 2*exp((-1/2)*δ^2/λ^2)*δ/λ^3 - exp((-1/2)*δ^2/λ^2)*δ^3/λ^5
    end

    if dorder_view==[0, 1]
        func = exp((-1/2)*δ^2/λ^2)*δ^2/λ^3
    end

    if dorder_view==[6, 0]
        func = -15*exp((-1/2)*δ^2/λ^2)/λ^6 + 45*exp((-1/2)*δ^2/λ^2)*δ^2/λ^8 - 15*exp((-1/2)*δ^2/λ^2)*δ^4/λ^10 + exp((-1/2)*δ^2/λ^2)*δ^6/λ^12
    end

    if dorder_view==[5, 0]
        func = -15*exp((-1/2)*δ^2/λ^2)*δ/λ^6 + 10*exp((-1/2)*δ^2/λ^2)*δ^3/λ^8 - exp((-1/2)*δ^2/λ^2)*δ^5/λ^10
    end

    if dorder_view==[4, 0]
        func = 3*exp((-1/2)*δ^2/λ^2)/λ^4 - 6*exp((-1/2)*δ^2/λ^2)*δ^2/λ^6 + exp((-1/2)*δ^2/λ^2)*δ^4/λ^8
    end

    if dorder_view==[3, 0]
        func = 3*exp((-1/2)*δ^2/λ^2)*δ/λ^4 - exp((-1/2)*δ^2/λ^2)*δ^3/λ^6
    end

    if dorder_view==[2, 0]
        func = -exp((-1/2)*δ^2/λ^2)/λ^2 + exp((-1/2)*δ^2/λ^2)*δ^2/λ^4
    end

    if dorder_view==[1, 0]
        func = -exp((-1/2)*δ^2/λ^2)*δ/λ^2
    end

    if dorder_view==[0, 0]
        func = exp((-1/2)*δ^2/λ^2)
    end

    dorder[2] = dorder2  # resetting dorder[2]
    return powers_of_negative_one(dorder2) * float(func)  # correcting for amount of t2 derivatives

end

return se_kernel, 1  # the function handle and the number of kernel hyperparameters


text:
# BigBed File Format
# ==================
#
# The BigWig/BigBed format is documented in
#   Kent, W. James, et al. "BigWig and BigBed: # enabling browsing of large
#   distributed datasets." Bioinformatics 26.17 (2010): 2204-2207.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module BigBed

import Automa
import Automa.RegExp: @re_str
import BufferedStreams
import ColorTypes
import FixedPointNumbers: N0f8
import GenomicFeatures: GenomicFeatures, BBI, BED, Interval
import Libz
using BioCore

include("reader.jl")
include("writer.jl")
include("record.jl")
include("overlap.jl")

end

ext: jl


text:
function isreducible(m::Module, v)
    inputs = map(x -> x.name, getinputs(m.dfg, v))

    return all(x -> haskey(m.parameters, x), inputs)
end

function reducenode!(m::Module, netlist::Netlist, v)
    # determine constant
    args = map(x -> m.parameters[x.name], getinputs(m.dfg, v))
    op = getoperator(m.dfg, v)
    value = eval(:($op($(args...))))

    # double check result
    isreal(value) || error("""
        Tried to reduce node, but received non-numeric result: $value
            args: $args
            op: $op
        """)

    # delete node and propagate constant
    children = outneighbors(m.dfg, v)
    outputs = getoutputs(m.dfg, v)
    for child in children
        for output in outputs
            inputs = getinputs(m.dfg, child)
            for (i, input) in enumerate(inputs)
                if getname(input) ∈ getname.(outputs)
                    inputs[i] = Variable(Symbol(value), Symbol(typeof(value)))
                end
            end
            set_prop!(m.dfg, child, :inputs, inputs)
        end
    end
    rem_vertex!(m.dfg, v)

    # add constant to netlist
    update!(netlist, Net(name = string(value), class = :constant, size = (1, 1)))

    return m
end

function constantreduction!(m::Module, netlist::Netlist)
    # if there are no nodes to reduce then return
    nodes = filter(v -> isreducible(m, v), getroots(m.dfg))
    isempty(nodes) && return m

    # reduce one node then try again
    reducenode!(m, netlist, nodes[1])
    constantreduction!(m, netlist)

    return m
end


text:
@doc raw"""
    AbstractMultinomialDoublyStochastic{N} <: AbstractEmbeddedManifold{ℝ, DefaultIsometricEmbeddingType}

A common type for manifolds that are doubly stochastic, for example by direct constraint
[`MultinomialDoubleStochastic`](@ref) or by symmetry [`MultinomialSymmetric`](@ref),
as long as they are also modeled as [`DefaultIsometricEmbeddingType`](@ref)
[`AbstractEmbeddedManifold`](@ref)s.
"""
abstract type AbstractMultinomialDoublyStochastic{N} <:
              AbstractEmbeddedManifold{ℝ,DefaultIsometricEmbeddingType} end

@doc raw"""
    MultinomialDoublyStochastic{n} <: AbstractMultinomialDoublyStochastic{N}

The set of doubly stochastic multinomial matrices consists of all $n×n$ matrices with
stochastic columns and rows, i.e.
````math
\begin{aligned}
\mathcal{DP}(n) \coloneqq \bigl\{p ∈ ℝ^{n×n}\ \big|\ &p_{i,j} > 0 \text{ for all } i=1,…,n, j=1,…,m,\\
& p\mathbf{1}_n = p^{\mathrm{T}}\mathbf{1}_n = \mathbf{1}_n
\bigr\},
\end{aligned}
````

where $\mathbf{1}_n$ is the vector of length $n$ containing ones.

The tangent space can be written as

````math
T_p\mathcal{DP}(n) \coloneqq \bigl\{
X ∈ ℝ^{n×n}\ \big|\ X = X^{\mathrm{T}} \text{ and }
X\mathbf{1}_n = X^{\mathrm{T}}\mathbf{1}_n = \mathbf{0}_n
\bigr\},
````

where $\mathbf{0}_n$ is the vector of length $n$ containing zeros.

More details can be found in Section III[^DouikHassibi2019].

# Constructor

    MultinomialDoubleStochastic(n)

Generate the manifold of matrices $\mathbb R^{n×n}$ that are doubly stochastic and symmetric.

[^DouikHassibi2019]:
    > A. Douik, B. Hassibi:
    > Manifold Optimization Over the Set of Doubly Stochastic Matrices: A Second-Order Geometry,
    > IEEE Transactions on Signal Processing 67(22), pp. 5761–5774, 2019.
    > doi: [10.1109/tsp.2019.2946024](http://doi.org/10.1109/tsp.2019.2946024),
    > arXiv: [1802.02628](https://arxiv.org/abs/1802.02628).
"""
struct MultinomialDoubleStochastic{N} <: AbstractMultinomialDoublyStochastic{N} end

function MultinomialDoubleStochastic(n::Int)
    return MultinomialDoubleStochastic{n}()
end

@doc raw"""
    check_manifold_point(M::MultinomialDoubleStochastic, p)

Checks whether `p` is a valid point on the [`MultinomialDoubleStochastic`](@ref)`(n)` `M`,
i.e. is a  matrix with positive entries whose rows and columns sum to one.
"""
function check_manifold_point(M::MultinomialDoubleStochastic{n}, p; kwargs...) where {n}
    mpv =
        invoke(check_manifold_point, Tuple{supertype(typeof(M)),typeof(p)}, M, p; kwargs...)
    mpv === nothing || return mpv
    # positivity and columns are checked in the embedding, we further check
    r = sum(p, dims = 2)
    if !isapprox(norm(r - ones(n, 1)), 0.0; kwargs...)
        return DomainError(
            r,
            "The point $(p) does not lie on $M, since its rows do not sum up to one.",
        )
    end
    return nothing
end
@doc raw"""
    check_tangent_vector(M::MultinomialDoubleStochastic p, X; check_base_point = true, kwargs...)

Checks whether `X` is a valid tangent vector to `p` on the [`MultinomialDoubleStochastic`](@ref) `M`.
This means, that `p` is valid, that `X` is of correct dimension and sums to zero along any
column or row.

The optional parameter `check_base_point` indicates, whether to call
[`check_manifold_point`](@ref check_manifold_point(::MultinomialDoubleStochastic, ::Any))  for `p`.
"""
function check_tangent_vector(
    M::MultinomialDoubleStochastic{n},
    p,
    X;
    check_base_point = true,
    kwargs...,
) where {n}
    if check_base_point
        mpe = check_manifold_point(M, p; kwargs...)
        mpe === nothing || return mpe
    end
    mpv = invoke(
        check_tangent_vector,
        Tuple{supertype(typeof(M)),typeof(p),typeof(X)},
        M,
        p,
        X;
        check_base_point = false, # already checked above
        kwargs...,
    )
    mpv === nothing || return mpv
    # columns are checked in the embedding, we further check
    r = sum(X, dims = 2) # check for stochastic rows
    if !isapprox(norm(r), 0.0; kwargs...)
        return DomainError(
            r,
            "The matrix $(X) is not a tangent vector to $(p) on $(M), since its rows do not sum up to zero.",
        )
    end
    return nothing
end

function decorated_manifold(::MultinomialDoubleStochastic{N}) where {N}
    return MultinomialMatrices(N, N)
end

@doc raw"""
    manifold_dimension(M::MultinomialDoubleStochastic{n}) where {n}

returns the dimension of the [`MultinomialDoubleStochastic`](@ref) manifold
namely
````math
\operatorname{dim}_{\mathcal{DP}(n)} = (n-1)^2.
````
"""
@generated function manifold_dimension(::MultinomialDoubleStochastic{n}) where {n}
    return (n - 1)^2
end

@doc raw"""
    project(M::MultinomialDoubleStochastic{n}, p, Y) where {n}

Project `Y` onto the tangent space at `p` on the [`MultinomialDoubleStochastic`](@ref) `M`, return the result in `X`.
The formula reads
````math
    \operatorname{proj}_p(Y) = Y - (α\mathbf{1}_n^{\mathrm{T}} + \mathbf{1}_nβ^{\mathrm{T}}) ⊙ p,
````
where $⊙$ denotes the Hadamard or elementwise product and $\mathbb{1}_n$ is the vector of length $n$ containing ones.
The two vectors $α,β ∈ ℝ^{n×n}$ are computed as a solution (typically using the left pseudo inverse) of
````math
    \begin{pmatrix} I_n & p\\p^{\mathrm{T}} & I_n \end{pmatrix}
    \begin{pmatrix} α\\ β\end{pmatrix}
    =
    \begin{pmatrix} Y\mathbf{1}\\Y^{\mathrm{T}}\mathbf{1}\end{pmatrix},
````
where $I_n$ is the $n×n$ unit matrix and $\mathbf{1}_n$ is the vector of length $n$ containing ones.

"""
project(::MultinomialDoubleStochastic, ::Any, ::Any)

function project!(::MultinomialDoubleStochastic{n}, X, p, Y) where {n}
    ζ = [I p; p I] \ [sum(Y, dims = 2); sum(Y, dims = 1)'] # Formula (25) from 1802.02628
    return X .= Y .- (repeat(ζ[1:n], 1, 3) .+ repeat(ζ[(n + 1):end]', 3, 1)) .* p
end

@doc raw"""
    project(
        M::AbstractMultinomialDoublyStochastic,
        p;
        maxiter = 100,
        tolerance = eps(eltype(p))
    )

project a matrix `p` with positive entries applying Sinkhorn's algorithm.
Note that this projct method – different from the usual case, accepts keywords.
"""
function project(M::AbstractMultinomialDoublyStochastic, p; kwargs...)
    q = allocate_result(M, project, p)
    project!(M, q, p; kwargs...)
    return q
end

function project!(
    ::AbstractMultinomialDoublyStochastic{n},
    q,
    p;
    maxiter = 100,
    tolerance = eps(eltype(p)),
) where {n}
    any(p .<= 0) &&
        throw(DomainError("The matrix $p can not be projected, since it has nonpositive entries."))
    iter = 0
    d1 = sum(p, dims = 1)
    d2 = 1 ./ (p * d1')
    row = d2' * p
    gap = 2 * tolerance
    while iter < maxiter && (gap >= tolerance)
        iter += 1
        row .= d2' * p
        gap = maximum(abs.(row .* d1 .- 1))
        d1 .= 1 ./ row
        d2 .= 1 ./ (p * d1')
    end
    q .= p .* (d2 * d1)
    return q
end

@generated function representation_size(::MultinomialDoubleStochastic{n}) where {n}
    return (n, n)
end

@doc raw"""
    retract(M::MultinomialDoubleStochastic, p, X, ::ProjectionRetraction)

compute a projection based retraction by projecting $p\odot\exp(X⨸p)$ back onto the manifold,
where $⊙,⨸$ are elementwise multiplication and division, respectively. Similarly, $\exp$
refers to the elementwise exponentiation.
"""
retract(::MultinomialDoubleStochastic, ::Any, ::Any, ::ProjectionRetraction)

function retract!(M::MultinomialDoubleStochastic, q, p, X, ::ProjectionRetraction)
    return project!(M, q, p .* exp.(X ./ p))
end

"""
    vector_transport_to(M::MultinomialDoubleStochastic, p, X, q)

transport the tangent vector `X` at `p` to `q` by projecting it onto the tangent space
at `q`.
"""
vector_transport_to(
    ::MultinomialDoubleStochastic,
    ::Any,
    ::Any,
    ::Any,
    ::ProjectionTransport,
)

function vector_transport_to!(
    M::MultinomialDoubleStochastic,
    Y,
    p,
    X,
    q,
    ::ProjectionTransport,
)
    project!(M, Y, q, X)
    return Y
end

function Base.show(io::IO, ::MultinomialDoubleStochastic{n}) where {n}
    return print(io, "MultinomialDoubleStochastic($(n))")
end


text:
@testset "CIR approximate filtering tests" begin
    Random.seed!(4)

    X = DualOptimalFiltering.generate_CIR_trajectory(range(0, stop = 2, length = 20), 3, 3., 0.5, 1);
    Y = map(λ -> rand(Poisson(λ),10), X);
    data = Dict(zip(range(0, stop = 2, length = 20), Y))

    Λ_of_t, wms_of_t, θ_of_t = DualOptimalFiltering.filter_CIR_keep_fixed_number(3., 0.5, 1.,1., data, 10);

    @test length(keys(Λ_of_t)) == 20
    @test length(keys(wms_of_t)) == 20

    Λ_of_t, wms_of_t = DualOptimalFiltering.filter_CIR_keep_above_threshold(3., 0.5, 1.,1., data, 0.0001)
    @test length(keys(Λ_of_t)) == 20
    @test length(keys(wms_of_t)) == 20

    Λ_of_t, wms_of_t = DualOptimalFiltering.filter_CIR_fixed_fraction(3., 0.5, 1.,1., data, .99)
    @test length(keys(Λ_of_t)) == 20
    @test length(keys(wms_of_t)) == 20

    fixed_number = 10
    function prune_keeping_fixed_number(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = DualOptimalFiltering.keep_fixed_number_of_weights(Λ_of_t, wms_of_t, fixed_number)
        return Λ_of_t_kept, DualOptimalFiltering.normalise(wms_of_t_kept)
    end

end;

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
using PerceptualColourMaps
using Base.Test

include("test_cmap.jl")
include("test_applycolourmaps.jl")
include("test_relief.jl")
include("test_utilities.jl")


text:
using PyPlot

## Draw Graphs
Graph=digraph[1,2]
## Add all nodes to appropriate place
for n=2:Nodes
    Graph=addedge[Graph,rainbow[n,1]+1,n+1]
end
p=plot(Graph,"Layout','force','NodeLabel",{})
## Color the nodes appropriately
highlight[p,1,"NodeColor','k"]
labelnode[p,1,"SEED"]
for n=2:length(rainbow)+1
    if rainbow[n-1,2]==1
        NODECOLOR='r'
    elseif rainbow[n-1,2]==2
        NODECOLOR=[1,0.5,0]
    elseif rainbow[n-1,2]==3
        NODECOLOR='y'
    elseif rainbow[n-1,2]==4
        NODECOLOR='g'
    elseif rainbow[n-1,2]==5
        NODECOLOR='b'
    elseif rainbow[n-1,2]==6
        NODECOLOR=[.75,0.25,.9]
    end
    highlight[p,n,"NodeColor",NODECOLOR]
    labelnode[p,n,n-1]
end

max_stars_repo_name: DesignEngrLab/ColorGraph

max_issues_repo_name: DesignEngrLab/ColorGraph

max_forks_repo_name: DesignEngrLab/ColorGraph


text:
using Oceananigans.Operators: assumed_field_location

#####
##### Default boundary conditions
#####

struct DefaultPrognosticFieldBoundaryCondition end

default_prognostic_bc(::Grids.Periodic, loc)       = PeriodicBoundaryCondition()
default_prognostic_bc(::Bounded,        ::Center)  = NoFluxBoundaryCondition()
default_prognostic_bc(::Bounded,        ::Face)    = ImpenetrableBoundaryCondition()

default_prognostic_bc(::Bounded,        ::Nothing) = nothing
default_prognostic_bc(::Flat,           ::Nothing) = nothing
default_prognostic_bc(::Grids.Periodic, ::Nothing) = nothing
default_prognostic_bc(::Flat, loc) = nothing

default_auxiliary_bc(topo, loc) = default_prognostic_bc(topo, loc)
default_auxiliary_bc(::Bounded, ::Face) = nothing

#####
##### Field boundary conditions
#####

mutable struct FieldBoundaryConditions{W, E, S, N, B, T, I}
        west :: W
        east :: E
       south :: S
       north :: N
      bottom :: B
         top :: T
    immersed :: I
end

function FieldBoundaryConditions(indices::Tuple, west, east, south, north, bottom, top, immersed)
    # Turn bcs in windowed dimensions into nothing
    west, east   = window_boundary_conditions(indices[1], west, east)
    south, north = window_boundary_conditions(indices[2], south, north)
    bottom, top  = window_boundary_conditions(indices[3], bottom, top)
    return FieldBoundaryConditions(west, east, south, north, bottom, top, immersed)
end

FieldBoundaryConditions(indices::Tuple, bcs::FieldBoundaryConditions) =
    FieldBoundaryConditions(indices, (getproperty(bcs, side) for side in fieldnames(FieldBoundaryConditions))...)

window_boundary_conditions(::Colon, left, right) = left, right
window_boundary_conditions(::UnitRange, left, right) = nothing, nothing

"""
    FieldBoundaryConditions(; kwargs...)

Return a template for boundary conditions on prognostic fields.

Keyword arguments
=================

Keyword arguments specify boundary conditions on the 7 possible boundaries:

  - `west`, left end point in the `x`-direction where `i=1`
  - `east`, right end point in the `x`-direction where `i=grid.Nx`
  - `south`, left end point in the `y`-direction where `j=1`
  - `north`, right end point in the `y`-direction where `j=grid.Ny`
  - `bottom`, right end point in the `z`-direction where `k=1`
  - `top`, right end point in the `z`-direction where `k=grid.Nz`
  - `immersed`, boundary between solid and fluid for immersed boundaries (experimental support only)

If a boundary condition is unspecified, the default for prognostic fields
and the topology in the boundary-normal direction is used:

  - `PeriodicBoundaryCondition` for `Periodic` directions
  - `NoFluxBoundaryCondition` for `Bounded` directions and `Centered`-located fields
  - `ImpenetrableBoundaryCondition` for `Bounded` directions and `Face`-located fields
  - `nothing` for `Flat` directions and/or `Nothing`-located fields
"""
function FieldBoundaryConditions(;     west = DefaultPrognosticFieldBoundaryCondition(),
                                       east = DefaultPrognosticFieldBoundaryCondition(),
                                      south = DefaultPrognosticFieldBoundaryCondition(),
                                      north = DefaultPrognosticFieldBoundaryCondition(),
                                     bottom = DefaultPrognosticFieldBoundaryCondition(),
                                        top = DefaultPrognosticFieldBoundaryCondition(),
                                   immersed = NoFluxBoundaryCondition())

   return FieldBoundaryConditions(west, east, south, north, bottom, top, immersed)
end

"""
    FieldBoundaryConditions(grid, location; kwargs...)

Return boundary conditions for auxiliary fields (fields whose values are
derived from a model's prognostic fields) on `grid` and at `location`.

Keyword arguments
=================

Keyword arguments specify boundary conditions on the 6 possible boundaries:

  - `west`, left end point in the `x`-direction where `i=1`
  - `east`, right end point in the `x`-direction where `i=grid.Nx`
  - `south`, left end point in the `y`-direction where `j=1`
  - `north`, right end point in the `y`-direction where `j=grid.Ny`
  - `bottom`, right end point in the `z`-direction where `k=1`
  - `top`, right end point in the `z`-direction where `k=grid.Nz`

If a boundary condition is unspecified, the default for auxiliary fields
and the topology in the boundary-normal direction is used:

  - `PeriodicBoundaryCondition` for `Periodic` directions
  - `GradientBoundaryCondition(0)` for `Bounded` directions and `Centered`-located fields
  - `nothing` for `Bounded` directions and `Face`-located fields
  - `nothing` for `Flat` directions and/or `Nothing`-located fields)
"""
function FieldBoundaryConditions(grid, loc, indices=(:, :, :);
                                 west     = default_auxiliary_bc(topology(grid, 1)(), loc[1]()),
                                 east     = default_auxiliary_bc(topology(grid, 1)(), loc[1]()),
                                 south    = default_auxiliary_bc(topology(grid, 2)(), loc[2]()),
                                 north    = default_auxiliary_bc(topology(grid, 2)(), loc[2]()),
                                 bottom   = default_auxiliary_bc(topology(grid, 3)(), loc[3]()),
                                 top      = default_auxiliary_bc(topology(grid, 3)(), loc[3]()),
                                 immersed = NoFluxBoundaryCondition())

    return FieldBoundaryConditions(indices, west, east, south, north, bottom, top, immersed)
end

#####
##### Boundary condition "regularization"
#####

regularize_boundary_condition(::DefaultPrognosticFieldBoundaryCondition, topo, loc, dim, args...) =
    default_prognostic_bc(topo[dim](), loc[dim]())

regularize_boundary_condition(bc, args...) = bc # fallback

"""
Compute default boundary conditions and attach field locations to ContinuousBoundaryFunction
boundary conditions for prognostic model field boundary conditions.

!!! warn "No support for `ContinuousBoundaryFunction` for immersed boundary conditions"
    Do not regularize immersed boundary conditions.

    Currently, there is no support `ContinuousBoundaryFunction` for immersed boundary
    conditions.
"""
function regularize_field_boundary_conditions(bcs::FieldBoundaryConditions,
                                              grid::AbstractGrid,
                                              field_name::Symbol,
                                              prognostic_names=nothing)
    topo = topology(grid)
    loc = assumed_field_location(field_name)

    west     = regularize_boundary_condition(bcs.west,   topo, loc, 1, 1,       prognostic_names)
    east     = regularize_boundary_condition(bcs.east,   topo, loc, 1, grid.Nx, prognostic_names)
    south    = regularize_boundary_condition(bcs.south,  topo, loc, 2, 1,       prognostic_names)
    north    = regularize_boundary_condition(bcs.north,  topo, loc, 2, grid.Ny, prognostic_names)
    bottom   = regularize_boundary_condition(bcs.bottom, topo, loc, 3, 1,       prognostic_names)
    top      = regularize_boundary_condition(bcs.top,    topo, loc, 3, grid.Nz, prognostic_names)

    # Eventually we could envision supporting ContinuousForcing-style boundary conditions
    # for the immersed boundary condition, which would benefit from regularization.
    # But for now we don't regularize the immersed boundary condition.
    immersed = bcs.immersed

    return FieldBoundaryConditions(west, east, south, north, bottom, top, immersed)
end

# For nested NamedTuples of boundary conditions (eg diffusivity boundary conditions)
function regularize_field_boundary_conditions(boundary_conditions::NamedTuple,
                                              grid::AbstractGrid,
                                              group_name::Symbol,
                                              prognostic_names=nothing)

    return NamedTuple(field_name => regularize_field_boundary_conditions(field_bcs, grid, field_name, prognostic_names)
                      for (field_name, field_bcs) in pairs(boundary_conditions))
end

regularize_field_boundary_conditions(::Missing,
                                     grid::AbstractGrid,
                                     field_name::Symbol,
                                     prognostic_names=nothing) = missing

#####
##### Outer interface for model constructors
#####

regularize_field_boundary_conditions(boundary_conditions::NamedTuple, grid::AbstractGrid, prognostic_names::Tuple) =
    NamedTuple(field_name => regularize_field_boundary_conditions(field_bcs, grid, field_name, prognostic_names)
               for (field_name, field_bcs) in pairs(boundary_conditions))


text:
export birdcageSensitivity

"""
    birdcageSensitivity(N::Int64, ncoils::Int64, relative_radius::Float64)

Computes the sensitivity maps for each coils that are arranged
in a birdcage manner.
"""
function birdcageSensitivity(N::Int64,
                             ncoils::Int64,
                             relative_radius::Float64)

  out = zeros(ComplexF64, N, N, 1, ncoils)
  for c=0:ncoils-1
    coilx = relative_radius*cos(c*(2*pi/ncoils))
    coily = relative_radius*sin(c*(2*pi/ncoils))
    coil_phase = -c*(2*pi/ncoils)

    for y=0:N-1
      y_co = (y - (N/2))/(N/2) - coily
      for x=0:N-1
        x_co = (x - (N/2))/(N/2) - coilx
        rr = sqrt(x_co^2+y_co^2)
        phi = atan(x_co, -y_co) + coil_phase
        out[x+1,y+1, 1, c+1] = 1/(rr) * exp(1im*phi)
      end
    end
  end
 return out
end


text:
"""
Convert Julia arrays to and from PyTorch `Tensor`s. Dimension ordering.

What do arrays/tensors look like?

- In Julia, the batch dimension is the LAST, in PyTorch is FIRST.
- In Julia, the channel dimensions comes AFTER the spatial dimensions,
    in PyTorch it comes BEFORE.

Example permutations:

- Batch of images: HWCB -> BCHW (4,3,1,2)
- Batch of onehot vectors: CB -> BC (2,1)

"""

abstract type ArrayFormat end

struct BatchFormat{F<:ArrayFormat} <: ArrayFormat
    format::F
end

BatchFormat() = BatchFormat(SameFormat())

struct SameFormat <: ArrayFormat end
struct ImageFormat <: ArrayFormat end

tensorpermutation(f::BatchFormat, n::Int) = (n, (tensorpermutation(f.format, n-1))...)
tensorpermutation(::SameFormat, n::Int) = ntuple(identity, n)
tensorpermutation(::ImageFormat, n::Int) = (n, ntuple(identity, n-1)...)

defaultformat(::AbstractArray) = BatchFormat()
defaultformat(::AbstractArray{T, 4}) where T = BatchFormat(ImageFormat())

function totensor(a::AbstractArray, format = defaultformat(a); kwargs...)
    return torch.tensor(permutedims(a, tensorpermutation(format, ndims(a))); kwargs...)
end

totensor(t::PyObject, args...; kwargs...) = t

function fromtensor(a::PyObject, format = defaultformat(a))
    return pycall(a.cpu().numpy, PyArray)
end

@testset "tensorpermutation" begin
    @test tensorpermutation(BatchFormat(), 4) == (4, 1, 2, 3)
    @test tensorpermutation(SameFormat(), 4) == (1, 2, 3, 4)
    @test tensorpermutation(ImageFormat(), 3) == (3, 1, 2)
end


text:
using ECOS
const ECOS_PRIMAL_FACTORY =
    MOI.OptimizerWithAttributes(ECOS.Optimizer, MOI.Silent() => true)
const ECOS_DUAL_FACTORY = dual_optimizer(ECOS_PRIMAL_FACTORY)
const ECOS_PRIMAL_OPT = MOI.instantiate(ECOS_PRIMAL_FACTORY)
const ECOS_DUAL_OPT = MOI.instantiate(ECOS_DUAL_FACTORY)

# push!(dual_linear_optimizer, DualOptimizer(ECOS.Optimizer(verbose = 0)))

@testset "ECOS conic Problems" begin
    @testset "ECOS Conic linear problems" begin
        list_of_soc_problems = [
            conic_linear1_test,
            conic_linear2_test,
            conic_linear3_test,
            conic_linear4_test,
        ]
        test_strong_duality(list_of_soc_problems, ECOS_PRIMAL_FACTORY)
    end
    @testset "ECOS SOC problems" begin
        list_of_soc_problems =
            [soc1_test, soc2_test, soc3_test, soc4_test, soc5_test, soc6_test]
        test_strong_duality(list_of_soc_problems, ECOS_PRIMAL_FACTORY)
    end
    @testset "ECOS RotatedSOC problems" begin
        list_of_rsoc_problems = [rsoc1_test, rsoc2_test, rsoc3_test, rsoc4_test]
        test_strong_duality(list_of_rsoc_problems, ECOS_PRIMAL_FACTORY)
    end
end


text:
struct VenusBarycenter <: abstractBarycenter end

naif_id(::Type{VenusBarycenter}) = 2
parent(::Type{VenusBarycenter}) = SSB


text:
using Test
using QCMaterial

@testset "computation.apply_qubit_op" begin
    n_qubit = 2
    # Prepare |Psi> = |00>
    state = QulacsQuantumState(n_qubit)
    set_computational_basis!(state, 0b00)

    # Apply c^dagger_1 to |Psi>, which will yield |10>.
    op = jordan_wigner(FermionOperator("2^"))

    # Prepare |phi> = |01>
    state0_bra = QulacsQuantumState(n_qubit)
    set_computational_basis!(state0_bra, 0b01)

    # Fit <01|U(theta)^dagger c_2^dagger |00>

    circuit_bra = uccgsd(n_qubit, orbital_rot=true, conserv_Sz_singles=false)

    update_circuit_param!(circuit_bra, rand(num_theta(circuit_bra)))
    #update_circuit_param!(circuit_bra, rand(size(circuit_bra.theta_offset)[1]))
    squared_norm = apply_qubit_op!(op, state, circuit_bra, state0_bra)
    println("squared_norm=",squared_norm)

    # Verify the result
    #@test≈ 1.0
    @test isapprox( abs(squared_norm), 1.0, rtol=1e-3)
end

#@testset "computation.apply_qubit_op2" begin
    #n_qubit = 2
    ## Prepare |Psi> = (|01> + |10>)/sqrt(2)
    #state = QulacsQuantumState(n_qubit)
    #set_computational_basis!(state, 0b01)
#
    #state2 = QulacsQuantumState(n_qubit)
    #set_computational_basis!(state, 0b10)
#
    ##state + state2
#
#end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
function approximate_p(::Type{InferenceTestUtils}, arg0::jdouble, arg1::jint, arg2::jint)
    return jcall(InferenceTestUtils, "approximateP", jdouble, (jdouble, jint, jint), arg0, arg1, arg2)
end

function chi_square(::Type{InferenceTestUtils}, arg0::Vector{Vector{jlong}})
    return jcall(InferenceTestUtils, "chiSquare", jdouble, (Vector{Vector{jlong}},), arg0)
end

function chi_square(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jlong})
    return jcall(InferenceTestUtils, "chiSquare", jdouble, (Vector{jdouble}, Vector{jlong}), arg0, arg1)
end

function chi_square_data_sets_comparison(::Type{InferenceTestUtils}, arg0::Vector{jlong}, arg1::Vector{jlong})
    return jcall(InferenceTestUtils, "chiSquareDataSetsComparison", jdouble, (Vector{jlong}, Vector{jlong}), arg0, arg1)
end

function chi_square_test(::Type{InferenceTestUtils}, arg0::Vector{Vector{jlong}})
    return jcall(InferenceTestUtils, "chiSquareTest", jdouble, (Vector{Vector{jlong}},), arg0)
end

function chi_square_test(::Type{InferenceTestUtils}, arg0::Vector{Vector{jlong}}, arg1::jdouble)
    return jcall(InferenceTestUtils, "chiSquareTest", jboolean, (Vector{Vector{jlong}}, jdouble), arg0, arg1)
end

function chi_square_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jlong})
    return jcall(InferenceTestUtils, "chiSquareTest", jdouble, (Vector{jdouble}, Vector{jlong}), arg0, arg1)
end

function chi_square_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jlong}, arg2::jdouble)
    return jcall(InferenceTestUtils, "chiSquareTest", jboolean, (Vector{jdouble}, Vector{jlong}, jdouble), arg0, arg1, arg2)
end

function chi_square_test_data_sets_comparison(::Type{InferenceTestUtils}, arg0::Vector{jlong}, arg1::Vector{jlong})
    return jcall(InferenceTestUtils, "chiSquareTestDataSetsComparison", jdouble, (Vector{jlong}, Vector{jlong}), arg0, arg1)
end

function chi_square_test_data_sets_comparison(::Type{InferenceTestUtils}, arg0::Vector{jlong}, arg1::Vector{jlong}, arg2::jdouble)
    return jcall(InferenceTestUtils, "chiSquareTestDataSetsComparison", jboolean, (Vector{jlong}, Vector{jlong}, jdouble), arg0, arg1, arg2)
end

function exact_p(::Type{InferenceTestUtils}, arg0::jdouble, arg1::jint, arg2::jint, arg3::jboolean)
    return jcall(InferenceTestUtils, "exactP", jdouble, (jdouble, jint, jint, jboolean), arg0, arg1, arg2, arg3)
end

function g(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jlong})
    return jcall(InferenceTestUtils, "g", jdouble, (Vector{jdouble}, Vector{jlong}), arg0, arg1)
end

function g_data_sets_comparison(::Type{InferenceTestUtils}, arg0::Vector{jlong}, arg1::Vector{jlong})
    return jcall(InferenceTestUtils, "gDataSetsComparison", jdouble, (Vector{jlong}, Vector{jlong}), arg0, arg1)
end

function g_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jlong})
    return jcall(InferenceTestUtils, "gTest", jdouble, (Vector{jdouble}, Vector{jlong}), arg0, arg1)
end

function g_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jlong}, arg2::jdouble)
    return jcall(InferenceTestUtils, "gTest", jboolean, (Vector{jdouble}, Vector{jlong}, jdouble), arg0, arg1, arg2)
end

function g_test_data_sets_comparison(::Type{InferenceTestUtils}, arg0::Vector{jlong}, arg1::Vector{jlong})
    return jcall(InferenceTestUtils, "gTestDataSetsComparison", jdouble, (Vector{jlong}, Vector{jlong}), arg0, arg1)
end

function g_test_data_sets_comparison(::Type{InferenceTestUtils}, arg0::Vector{jlong}, arg1::Vector{jlong}, arg2::jdouble)
    return jcall(InferenceTestUtils, "gTestDataSetsComparison", jboolean, (Vector{jlong}, Vector{jlong}, jdouble), arg0, arg1, arg2)
end

function g_test_intrinsic(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jlong})
    return jcall(InferenceTestUtils, "gTestIntrinsic", jdouble, (Vector{jdouble}, Vector{jlong}), arg0, arg1)
end

function homoscedastic_t(::Type{InferenceTestUtils}, arg0::StatisticalSummary, arg1::StatisticalSummary)
    return jcall(InferenceTestUtils, "homoscedasticT", jdouble, (StatisticalSummary, StatisticalSummary), arg0, arg1)
end

function homoscedastic_t(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "homoscedasticT", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function homoscedastic_t_test(::Type{InferenceTestUtils}, arg0::StatisticalSummary, arg1::StatisticalSummary)
    return jcall(InferenceTestUtils, "homoscedasticTTest", jdouble, (StatisticalSummary, StatisticalSummary), arg0, arg1)
end

function homoscedastic_t_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "homoscedasticTTest", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function homoscedastic_t_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(InferenceTestUtils, "homoscedasticTTest", jboolean, (Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function kolmogorov_smirnov_statistic(::Type{InferenceTestUtils}, arg0::RealDistribution, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "kolmogorovSmirnovStatistic", jdouble, (RealDistribution, Vector{jdouble}), arg0, arg1)
end

function kolmogorov_smirnov_statistic(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "kolmogorovSmirnovStatistic", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function kolmogorov_smirnov_test(::Type{InferenceTestUtils}, arg0::RealDistribution, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "kolmogorovSmirnovTest", jdouble, (RealDistribution, Vector{jdouble}), arg0, arg1)
end

function kolmogorov_smirnov_test(::Type{InferenceTestUtils}, arg0::RealDistribution, arg1::Vector{jdouble}, arg2::jboolean)
    return jcall(InferenceTestUtils, "kolmogorovSmirnovTest", jdouble, (RealDistribution, Vector{jdouble}, jboolean), arg0, arg1, arg2)
end

function kolmogorov_smirnov_test(::Type{InferenceTestUtils}, arg0::RealDistribution, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(InferenceTestUtils, "kolmogorovSmirnovTest", jboolean, (RealDistribution, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function kolmogorov_smirnov_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "kolmogorovSmirnovTest", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function kolmogorov_smirnov_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jboolean)
    return jcall(InferenceTestUtils, "kolmogorovSmirnovTest", jdouble, (Vector{jdouble}, Vector{jdouble}, jboolean), arg0, arg1, arg2)
end

function one_way_anova_f_value(::Type{InferenceTestUtils}, arg0::Collection)
    return jcall(InferenceTestUtils, "oneWayAnovaFValue", jdouble, (Collection,), arg0)
end

function one_way_anova_p_value(::Type{InferenceTestUtils}, arg0::Collection)
    return jcall(InferenceTestUtils, "oneWayAnovaPValue", jdouble, (Collection,), arg0)
end

function one_way_anova_test(::Type{InferenceTestUtils}, arg0::Collection, arg1::jdouble)
    return jcall(InferenceTestUtils, "oneWayAnovaTest", jboolean, (Collection, jdouble), arg0, arg1)
end

function paired_t(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "pairedT", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function paired_t_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "pairedTTest", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function paired_t_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(InferenceTestUtils, "pairedTTest", jboolean, (Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function root_log_likelihood_ratio(::Type{InferenceTestUtils}, arg0::jlong, arg1::jlong, arg2::jlong, arg3::jlong)
    return jcall(InferenceTestUtils, "rootLogLikelihoodRatio", jdouble, (jlong, jlong, jlong, jlong), arg0, arg1, arg2, arg3)
end

function t(::Type{InferenceTestUtils}, arg0::StatisticalSummary, arg1::StatisticalSummary)
    return jcall(InferenceTestUtils, "t", jdouble, (StatisticalSummary, StatisticalSummary), arg0, arg1)
end

function t(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "t", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function t(::Type{InferenceTestUtils}, arg0::jdouble, arg1::StatisticalSummary)
    return jcall(InferenceTestUtils, "t", jdouble, (jdouble, StatisticalSummary), arg0, arg1)
end

function t(::Type{InferenceTestUtils}, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "t", jdouble, (jdouble, Vector{jdouble}), arg0, arg1)
end

function t_test(::Type{InferenceTestUtils}, arg0::StatisticalSummary, arg1::StatisticalSummary)
    return jcall(InferenceTestUtils, "tTest", jdouble, (StatisticalSummary, StatisticalSummary), arg0, arg1)
end

function t_test(::Type{InferenceTestUtils}, arg0::StatisticalSummary, arg1::StatisticalSummary, arg2::jdouble)
    return jcall(InferenceTestUtils, "tTest", jboolean, (StatisticalSummary, StatisticalSummary, jdouble), arg0, arg1, arg2)
end

function t_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "tTest", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function t_test(::Type{InferenceTestUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(InferenceTestUtils, "tTest", jboolean, (Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function t_test(::Type{InferenceTestUtils}, arg0::jdouble, arg1::StatisticalSummary)
    return jcall(InferenceTestUtils, "tTest", jdouble, (jdouble, StatisticalSummary), arg0, arg1)
end

function t_test(::Type{InferenceTestUtils}, arg0::jdouble, arg1::StatisticalSummary, arg2::jdouble)
    return jcall(InferenceTestUtils, "tTest", jboolean, (jdouble, StatisticalSummary, jdouble), arg0, arg1, arg2)
end

function t_test(::Type{InferenceTestUtils}, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(InferenceTestUtils, "tTest", jdouble, (jdouble, Vector{jdouble}), arg0, arg1)
end

function t_test(::Type{InferenceTestUtils}, arg0::jdouble, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(InferenceTestUtils, "tTest", jboolean, (jdouble, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end


text:
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization

@testset "Immersed boundaries with hydrostatic free surface models" begin
    @info "Testing immersed boundaries vertical integrals"

    for arch in archs
        Nx = 5
        Ny = 5

        underlying_grid = RectilinearGrid(arch,
                                          size = (Nx, Ny, 3),
                                          extent = (Nx, Ny, 3),
                                          topology = (Periodic, Periodic, Bounded))

        # B for bathymetry
        B = [-3. for i=1:Nx, j=1:Ny ]
        B[2:Nx-1,2:Ny-1] .= [-2. for i=2:Nx-1, j=2:Ny-1 ]
        B[3:Nx-2,3:Ny-2] .= [-1. for i=3:Nx-2, j=3:Ny-2 ]

        grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(B))

        model = HydrostaticFreeSurfaceModel(grid = grid,
                                            free_surface = ImplicitFreeSurface(),
                                            tracer_advection = WENO5(),
                                            buoyancy = nothing,
                                            tracers = nothing,
                                            closure = nothing)

        x_ref = [0.0  0.0  0.0  0.0  0.0  0.0  0.0
                 0.0  3.0  3.0  3.0  3.0  3.0  0.0
                 0.0  3.0  2.0  2.0  2.0  2.0  0.0
                 0.0  3.0  2.0  1.0  1.0  2.0  0.0
                 0.0  3.0  2.0  2.0  2.0  2.0  0.0
                 0.0  3.0  3.0  3.0  3.0  3.0  0.0
                 0.0  0.0  0.0  0.0  0.0  0.0  0.0]'

        y_ref = [0.0  0.0  0.0  0.0  0.0  0.0  0.0
                 0.0  3.0  3.0  3.0  3.0  3.0  0.0
                 0.0  3.0  2.0  2.0  2.0  3.0  0.0
                 0.0  3.0  2.0  1.0  2.0  3.0  0.0
                 0.0  3.0  2.0  1.0  2.0  3.0  0.0
                 0.0  3.0  2.0  2.0  2.0  3.0  0.0
                 0.0  0.0  0.0  0.0  0.0  0.0  0.0]'

        fs = model.free_surface
        vertically_integrated_lateral_areas = fs.implicit_step_solver.vertically_integrated_lateral_areas

        ∫Axᶠᶜᶜ = vertically_integrated_lateral_areas.xᶠᶜᶜ
        ∫Ayᶜᶠᶜ = vertically_integrated_lateral_areas.yᶜᶠᶜ

        ∫Axᶠᶜᶜ = Array(parent(∫Axᶠᶜᶜ))
        ∫Ayᶜᶠᶜ = Array(parent(∫Ayᶜᶠᶜ))

        Ax_ok = ∫Axᶠᶜᶜ[:, :, 1] ≈ x_ref
        Ay_ok = ∫Ayᶜᶠᶜ[:, :, 1] ≈ y_ref

        @test (Ax_ok & Ay_ok)
    end
end


text:
using Pseudospectra, Test

# Toeplitz matrix from S&P, Ch. 16
# Need a large domain for lower bound computation.
@testset "Matrix powers" begin
    opts = Dict{Symbol,Any}(:npts => 40,:ax => [-110,110,-110,110])
    A = diagm(1 => 100.0*ones(4)) + diagm(-2 => 1e-4*ones(3))
    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)
    @test iscomputed(ps_data)

    mtxpowersplot(gs,ps_data,20,lbmethod=:combined)
end


text:
using Bokeh

# more or less the simplest plotting example
x = linspace(0, 2pi, 10)
y = sin(x)#[sin(x) cos(x)]
# by default autoopen is set to true so `plot` will try to open a browser
# with the plot file (which by default is "bokeh_plot.html")
Bokeh.debug(true)
# plot(y, "k--")
Bokeh.height(400)
Bokeh.width(600)
hold(true)
autoopen(false)
ys = [sin(x) cos(x)]
plot(ys)
showplot()


text:
#=
Function declarations for PCRE2 library API

Copyright 2018 Gandalf Software, Inc., Scott P. Jones, and contributors to pcre2.h
(based on pcre2.h (copyright University of Cambridge), see PCRE2_LICENSE)
Licensed under MIT License, see LICENSE.md
=#

struct _UCHAR end

## Type definitions for PCRE bindings

const UCharP          = Ptr{_UCHAR} # This is really Ptr{UInt8},Ptr{UInt16},Ptr{UInt32}
const StrP            = UCharP
const VoidP           = Ptr{Cvoid}
const GeneralContextP = Ptr{Cvoid}
const CompileContextP = Ptr{Cvoid}
const MatchContextP   = Ptr{Cvoid}
const ConvertContextP = Ptr{Cvoid}
const MatchDataP      = Ptr{Cvoid}
const CodeP           = Ptr{Cvoid}
const JitStackP       = Ptr{Cvoid}
const SizeP           = Ptr{Csize_t}
const SizeRef         = Ref{Csize_t}
const MatchOptions    = UInt32
const CompileOptions  = UInt32

const funclist =
(
 (:config, Cint, (CONFIG, VoidP)),

## Functions for manipulating contexts

# GENERAL CONTEXT FUNCTIONS
 (:general_context_copy,      GeneralContextP, (GeneralContextP,)),
 #(:general_context_create, GeneralContextP,
 #    (VoidP (*)(Csize_t, VoidP), Cvoid (*)(VoidP, VoidP), VoidP)),
 (:general_context_free,      Cvoid,           (GeneralContextP,)),

# COMPILE CONTEXT FUNCTIONS
 (:compile_context_copy,      CompileContextP, (CompileContextP,)),
 (:compile_context_create,    CompileContextP, (GeneralContextP,)),
 (:compile_context_free,      Cvoid,           (CompileContextP,)),
 (:set_bsr,                   Cint,            (CompileContextP, UInt32)),
 (:set_character_tables,      Cint,            (CompileContextP, Ptr{UInt8})),
 (:set_compile_extra_options, Cint,            (CompileContextP, UInt32)),
 (:set_max_pattern_length,    Cint,            (CompileContextP, Csize_t)),
 (:set_newline,               Cint,            (CompileContextP, UInt32)),
 (:set_parens_nest_limit,     Cint,            (CompileContextP, UInt32)),
 #(:set_compile_recursion_guard, Cint, (CompileContextP, Cint (*)(UInt32, VoidP), VoidP)

# MATCH CONTEXT FUNCTIONS
 (:match_context_copy,        MatchContextP,   (MatchContextP,)),
 (:match_context_create,      MatchContextP,   (GeneralContextP,)),
 (:match_context_free,        Cvoid,           (MatchContextP,)),
 #(:set_callout, Cint, (MatchContextP, Cint (*)(callout_block *, VoidP), VoidP)),
 (:set_depth_limit,           Cint,            (MatchContextP, UInt32)),
 (:set_heap_limit,            Cint,            (MatchContextP, UInt32)),
 (:set_match_limit,           Cint,            (MatchContextP, UInt32)),
 (:set_offset_limit,          Cint,            (MatchContextP, Csize_t)),
 (:set_recursion_limit,       Cint,            (MatchContextP, UInt32)),
 #(:set_recursion_memory_management, Cint,
 # (MatchContextP, VoidP (*)(Csize_t, VoidP), Cvoid (*)(VoidP, VoidP), VoidP)),

# CONVERT CONTEXT FUNCTIONS
 (:convert_context_copy,      ConvertContextP, (ConvertContextP,)),
 (:convert_context_create,    ConvertContextP, (GeneralContextP,)),
 (:convert_context_free,      Cvoid,           (ConvertContextP,)),
 (:set_glob_escape,           Cint,            (ConvertContextP, UInt32)),
 (:set_glob_separator,        Cint,            (ConvertContextP, UInt32)),

## Functions concerned with compiling a pattern to PCRE internal code

# COMPILE FUNCTIONS
 (:compile, CodeP, (StrP, Csize_t, CompileOptions, Ref{Cint}, SizeRef, CompileContextP)),
 (:code_free,                 Cvoid,           (CodeP,)),
 (:code_copy,                 CodeP,           (CodeP,)),
 (:code_copy_with_tables,     CodeP,           (CodeP,)),

## Functions that give information about a compiled pattern

# PATTERN INFO FUNCTIONS
 (:pattern_info,              Cint,            (CodeP, INFO, VoidP)),
#(:callout_enumerate, Cint, (CodeP, Cint (*)(callout_enumerate_block *, VoidP), VoidP)),

## Functions for running a match and inspecting the result.

# MATCH FUNCTIONS
 (:match_data_create,          MatchDataP,     (UInt32, GeneralContextP)),
 (:match_data_create_from_pattern, MatchDataP, (CodeP, GeneralContextP)),
 (:dfa_match,                  Cint,
  (CodeP, StrP, Csize_t, Csize_t, MatchOptions, MatchDataP, MatchContextP, Ptr{Cint}, Csize_t)),
 (:match,                      Cint,
  (CodeP, StrP, Csize_t, Csize_t, MatchOptions, MatchDataP, MatchContextP)),
 (:match_data_free,            Cvoid,          (MatchDataP,)),
 (:get_mark,                   StrP,           (MatchDataP,)),
 (:get_ovector_count,          UInt32,         (MatchDataP,)),
 (:get_ovector_pointer,        SizeP,          (MatchDataP,)),
 (:get_startchar,              Csize_t,        (MatchDataP,)),

## Convenience functions for handling matched substrings

# SUBSTRING FUNCTIONS
 (:substring_copy_byname,      Cint,           (MatchDataP, StrP, UCharP, SizeRef)),
 (:substring_copy_bynumber,    Cint,           (MatchDataP, UInt32, UCharP, SizeRef)),
 (:substring_free,             Cvoid,          (UCharP,)),
 (:substring_get_byname,       Cint,           (MatchDataP, StrP, Ref{UCharP}, SizeRef)),
 (:substring_get_bynumber,     Cint,           (MatchDataP, UInt32, Ref{UCharP}, SizeRef)),
 (:substring_length_byname,    Cint,           (MatchDataP, StrP, SizeRef)),
 (:substring_length_bynumber,  Cint,           (MatchDataP, UInt32, SizeRef)),
 (:substring_nametable_scan,   Cint,           (CodeP, StrP, Ptr{StrP}, Ptr{StrP})),
 (:substring_number_from_name, Cint,           (CodeP, StrP)),
 (:substring_list_free,        Cvoid,          (Ptr{StrP},)),
 (:substring_list_get,         Cint,           (MatchDataP, Ptr{Ptr{UCharP}}, Ptr{SizeP})),

## Functions for serializing / deserializing compiled patterns

# SERIALIZE FUNCTIONS
 (:serialize_encode,           Int32,
    (Ptr{CodeP}, Int32, Ptr{Ptr{UInt8}}, SizeP, GeneralContextP)),
 (:serialize_decode,           Int32,          (Ptr{CodeP}, Int32, Ptr{UInt8}, GeneralContextP)),
 (:serialize_get_number_of_codes, Int32,       (Ptr{UInt8},)),
 (:serialize_free,             Cvoid,          (Ptr{UInt8},)),

## Convenience function for match + substitute.

 (:substitute,                 Cint,
    (CodeP, StrP, Csize_t, Csize_t, MatchOptions, MatchDataP, MatchContextP,
     StrP, Csize_t, UCharP, SizeP)),

## Functions for converting pattern source strings

# CONVERT FUNCTIONS
 (:pattern_convert,            Cint, (StrP, Csize_t, UInt32, Ptr{UCharP}, SizeP, ConvertContextP)),
 (:converted_pattern_free,     Cvoid,          (UCharP,)),

## Functions for JIT processing

# JIT FUNCTIONS
 (:jit_compile,                Cint,           (CodeP, UInt32)),
 (:jit_match,                  Cint,
    (CodeP, StrP, Csize_t, Csize_t, UInt32, MatchDataP, MatchContextP)),
 (:jit_free_unused_memory,     Cvoid,          (GeneralContextP,)),
 (:jit_stack_create,           JitStackP,      (Csize_t, Csize_t, GeneralContextP)),
 (:jit_stack_assign,           Cvoid,          (MatchContextP, VoidP, VoidP)), # jit_callback
 (:jit_stack_free,             Cvoid,          (JitStackP,)),

## Other miscellaneous functions

 (:get_error_message,          Cint,           (Cint, UCharP, Csize_t)),
 (:maketables,                 Ptr{UInt8},     (GeneralContextP,))
)


text:
push!(LOAD_PATH,"../src/")
module Selects
using DataFrames
export Select
  function Select(Tabla::DataFrame,query::Vector{String}, nombres_de_columnas::Vector{String} ,Formato::String)::DataFrame
    #las columnas no se pueden repetir
    #
    #caso base por si las listas CSV no tienen header
    #
    #Checar si la columna esta repetida o no
    checker = Dict()

    # Procesar query
    for (indice,nombre_columna) in enumerate(query) # supongan que "1" es el nombre de la column      a
      idx = 0
      if !(nombre_columna in nombres_de_columnas)
        try
          idx = parse(Int32,nombre_columna)
        catch
          error("$nombre_columna no esta en al archivo CSV")
        end
        if idx > length(nombres_de_columnas)
          error("$idx esta fuera del rango de columnas: rango = 1-$(length(nombres_de_columnas))")
        else
          query[indice] = nombres_de_columnas[idx]
        end
      end
    end
    for nombre_columna in query
      #Checar si la columna esta en la tabla
      if !(nombre_columna in nombres_de_columnas)
        error("$x no se encuentra en el archivo CSV")
      end
      #Checar si la columna esta repetida
      if haskey(checker , nombre_columna)
        error("Columna $nombre_columna repetida")
      else
        push!(checker, nombre_columna => 1)
      end
    end
    query_inversa= Vector{String}()
    for nombre_columna in nombres_de_columnas
      if !(nombre_columna in query)
        push!(query_inversa,nombre_columna)
      end
    end
    q1 = select(Tabla, Not(query_inversa))
    # print(q1)
    return q1
  end
end

max_stars_repo_name: MarleenMs/Covid

max_issues_repo_name: MarleenMs/Covid

max_forks_repo_name: MarleenMs/Covid


text:
### A Pluto.jl notebook ###
# v0.17.2

using Markdown
using InteractiveUtils

using Pkg, DrWatson

begin
    using Distributions
    using StatsPlots
    using StatsBase
    using KernelDensity
    using StatisticalRethinking
end

md"##### Setting default attributes for plots."
default(labels=false)

md"### Snippet 3.1"

begin
    Pr_Positive_Vampire = 0.95
    Pr_Positive_Mortal = 0.01
    Pr_Vampire = 0.001
    tmp = Pr_Positive_Vampire * Pr_Vampire
    Pr_Positive = tmp + Pr_Positive_Mortal * (1 - Pr_Vampire)
    Pr_Vampire_Positive = tmp / Pr_Positive
    Pr_Vampire_Positive
end

md"### Snippet 3.2"

begin
    size = 1000
    p_grid = range(0, 1; length=size)
    prob_p = repeat([1.0], size);
    prob_data = [pdf(Binomial(9, p), 6) for p in p_grid];
    posterior = prob_data .* prob_p
    posterior /= sum(posterior);
end

md"### Snippet 3.3"

begin
    samples_count = 10_000
    cat = Categorical(posterior);
    indices = rand(cat, samples_count)
    samples = p_grid[indices];
end

md"### Snippet 3.4"

scatter(samples; alpha=0.2)

md"### Snippet 3.5"

density(samples)

md"### Snippet 3.6"

sum(posterior[p_grid .< 0.5])

md"### Snippet 3.7"

sum(samples .< 0.5) / samples_count

md"### Snippet 3.8"

sum(@. (samples > 0.5) & (samples < 0.75)) / samples_count

md"### Snippet 3.9"

quantile(samples, 0.8)

md"### Snippet 3.10"

quantile(samples, [0.1, 0.9])

md"### Snippet 3.11"

let
    size = 1000
    p_grid = range(0, 1; length=size)
    prob_p = repeat([1.0], size);
    prob_data = [pdf(Binomial(3, p), 3) for p in p_grid];
    posterior = prob_data .* prob_p
    posterior /= sum(posterior)

    samples_count = 10_000
    cat = Categorical(posterior);
    samples = p_grid[rand(cat, samples_count)];
end;

md"### Snippet 3.12"

percentile(samples, [25, 75])

md"### Snippet 3.13"

hpdi(samples, alpha=0.5)

md"### Snippet 3.14"

p_grid[argmax(posterior)]

md"### Snippet 3.15"

k = kde(samples, bandwidth=0.01)
k.x[argmax(k.density)]

md"### Snippet 3.16"

mean(samples), median(samples)

md"### Snippet 3.17"

sum(@. posterior * abs(0.5 - p_grid))

md"### Snippet 3.18"

loss = map(d -> sum(@. posterior * abs(d - p_grid)), p_grid);

md"### Snippet 3.19"

p_grid[argmin(loss)]

md"### Snippet 3.20"

[pdf(Binomial(2, 0.7), n) for n ∈ 0:2]

md"### Snippet 3.21"

rand(Binomial(2, 0.7))

md"### Snippet 3.22"

s = rand(Binomial(2, 0.7), 10)
println(s)

md"### Snippet 3.23"

let
    dummy_w = rand(Binomial(2, 0.7), 100_000);
    proportions(dummy_w)  # or counts(dummy_w)/100000
end

md"### Snippet 3.24"

let
    dummy_w = rand(Binomial(9, 0.7), 100_000);
    histogram(dummy_w; xlabel="dummy water count", ylabel="Frequency")
end

md"### Snippet 3.25"

w = rand(Binomial(9, 0.6), 10_000);

md"### Snippet 3.26"

w = [rand(Binomial(9, p)) for p in samples];

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╠═1f4452e6-57af-11ec-3c3a-9b6f651fda27
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
## Exercise 5-3
## Fermat’s Last Theorem says that there are no positive integers a, b, and c such that

## a^n+b^n=c^n
## for any value of n greater than 2.

## Write a function named checkfermat that takes four parameters—a, b, c and n—and checks to see if Fermat’s theorem holds. If n is greater than 2 and a^n + b^n == c^n the program should print, “Holy smokes, Fermat was wrong!” Otherwise the program should print, “No, that doesn’t work.”
println("Ans 1: ")

function check_fermats_theorem(a::Int, b::Int, c::Int, n::Int)
    if a^n + b^n == c^n && n > 2
        println("Holy smokes, Fermat was wrong!")
    else
        println("No, that doesn’t work.")
    end
end

check_fermats_theorem(3, 5, 3, 3)

## 2. Write a function that prompts the user to input values for a, b, c and n, converts them to integers, and uses checkfermat to check whether they violate Fermat’s theorem.
println("Ans 2: ")

function check_fermats_theorem2()
    print("Enter value for a: ")
    a = parse(Int64, readline())
    print("Enter value for b: ")
    b = parse(Int64, readline())
    print("Enter value for c: ")
    c = parse(Int64, readline())
    print("Enter value for n: ")
    n = parse(Int64, readline())

    check_fermats_theorem(a, b, c, n)
end

check_fermats_theorem2()

println("End.")


text:
include("../../src/ToyModelsMD.jl")
using .ToyModelsMD

mb = MullerBrown()
langevin = Langevin(1.0,1.0)
system = ThermostattedSystem(mb,langevin)

dt = 0.0001
tf = 0.0015
nsteps = convert(Int64, round(tf/dt))

system.x = [0.63, 0.03]
integrator = StochasticEuler(system,mb)
runtraj!(nsteps,dt,integrator)


text:
# reverse-mode evaluation of an expression tree

# assumes partials_storage is already updated
# dense gradient output, assumes initialized to zero
# if subexpressions are present, must run reverse_eval on subexpression tapes afterwards
function reverse_eval(reverse_storage::AbstractVector{T}, partials_storage::AbstractVector{T},
                      nd::Vector{NodeData},adj) where {T}

    @assert length(reverse_storage) >= length(nd)
    @assert length(partials_storage) >= length(nd)

    # nd is already in order such that parents always appear before children
    # so a forward pass through nd is a backwards pass through the tree

    # reverse_storage[k] is the partial derivative of the output with respect to
    # the value of node k
    reverse_storage[1] = one(T)

    for k in 2:length(nd)
        @inbounds nod = nd[k]
        if nod.nodetype == VALUE || nod.nodetype == LOGIC || nod.nodetype == COMPARISON || nod.nodetype == PARAMETER
            continue
        end
        @inbounds rev_parent = reverse_storage[nod.parent]
        @inbounds partial = partials_storage[k]
        @inbounds reverse_storage[k] = ifelse(rev_parent == 0.0 && !isfinite(partial), rev_parent, rev_parent*partial)
        #@inbounds reverse_storage[k] = reverse_storage[nod.parent]*partials_storage[k]
    end
    #@show storage

    nothing

end

export reverse_eval

# assume we've already run the reverse pass, now just extract the answer
# given the scaling value
function reverse_extract(output::AbstractVector{T}, reverse_storage::AbstractVector{T}, nd::Vector{NodeData},
                         adj, subexpression_output, scale_value::T) where {T}

    @assert length(reverse_storage) >= length(nd)

    for k in 1:length(nd)
        @inbounds nod = nd[k]
        if nod.nodetype == VARIABLE
            @inbounds output[nod.index] += scale_value*reverse_storage[k]
        elseif nod.nodetype == SUBEXPRESSION
            @inbounds subexpression_output[nod.index] += scale_value*reverse_storage[k]
        end
    end
    #@show storage

    nothing

end

export reverse_extract

# Compute directional derivatives of the reverse pass, goes with forward_eval_ϵ
# to compute hessian-vector products.
function reverse_eval_ϵ(output_ϵ::AbstractVector{ForwardDiff.Partials{N,T}},
                        reverse_storage::AbstractVector{T}, reverse_storage_ϵ, partials_storage::AbstractVector{T},
                        partials_storage_ϵ, nd::Vector{NodeData}, adj, subexpression_output,
                        subexpression_output_ϵ, scale_value::T, scale_value_ϵ::ForwardDiff.Partials{N,T}) where {N,T}

    @assert length(reverse_storage_ϵ) >= length(nd)
    @assert length(partials_storage_ϵ) >= length(nd)

    if nd[1].nodetype == VARIABLE
        @inbounds output_ϵ[nd[1].index] += scale_value_ϵ
        return
    elseif nd[1].nodetype == SUBEXPRESSION
        @inbounds subexpression_output[nd[1].index] += scale_value*reverse_storage[1]
        @inbounds subexpression_output_ϵ[nd[1].index] += scale_value_ϵ
        return
    end

    reverse_storage_ϵ[1] = scale_value_ϵ

    for k in 2:length(nd)
        @inbounds nod = nd[k]
        if nod.nodetype == VALUE || nod.nodetype == LOGIC || nod.nodetype == COMPARISON || nod.nodetype == PARAMETER
            continue
        end
        # compute the value of reverse_storage[k]
        @inbounds parentval = scale_value*reverse_storage[nod.parent]
        @inbounds parentval_ϵ = reverse_storage_ϵ[nod.parent]
        @inbounds partial = partials_storage[k]
        @inbounds partial_ϵ = partials_storage_ϵ[k]

        #reverse_storage_ϵ[k] = parentval*partial_ϵ + partial*parentval_ϵ
        if !isfinite(partial) && parentval == 0.0
            reverse_storage_ϵ[k] = zero(ForwardDiff.Partials{N,T})
        else
            reverse_storage_ϵ[k] = ForwardDiff._mul_partials(partial_ϵ,parentval_ϵ,parentval,partial)
        end

        if nod.nodetype == VARIABLE
            @inbounds output_ϵ[nod.index] += reverse_storage_ϵ[k]
        elseif nod.nodetype == SUBEXPRESSION
            @inbounds subexpression_output[nod.index] += scale_value*reverse_storage[k]
            @inbounds subexpression_output_ϵ[nod.index] += reverse_storage_ϵ[k]
        end
    end
    #@show storage

    nothing

end

export reverse_eval_ϵ


text:
"""
```julia
nanstd(A; dims=:, mean=nothing, corrected=true)
```
Compute the variance of all non-`NaN` elements in `A`, optionally over dimensions specified by `dims`.
As `Statistics.var`, but ignoring `NaN`s.

A precomputed `mean` may optionally be provided, which results in a somewhat faster
calculation. If `corrected` is `true`, then _Bessel's correction_ is applied, such
that the sum is divided by `n-1` rather than `n`.

As an alternative to `dims`, `nanstd` also supports the `dim` keyword, which
behaves identically to `dims`, but also drops any singleton dimensions that have
been reduced over (as is the convention in some other languages).

## Examples
```julia
julia> using NaNStatistics

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> nanstd(A, dims=1)
1×2 Matrix{Float64}:
 1.41421  1.41421

julia> nanstd(A, dims=2)
2×1 Matrix{Float64}:
 0.7071067811865476
 0.7071067811865476
```
"""
nanstd(A; dims=:, dim=:, mean=nothing, corrected=true) = sqrt!(__nanvar(mean, corrected, A, dims, dim))
export nanstd

sqrt!(x::Number) = sqrt(x)
function sqrt!(A::AbstractArray)
    @turbo for i ∈ eachindex(A)
        A[i] = sqrt(A[i])
    end
    return A
end


text:
module LibPQ

export status, reset!, execute, clear, fetch!, prepare,
    num_columns, num_rows, num_params, num_affected_rows

using Compat.Dates
using DocStringExtensions
using Decimals
using DataStreams
using Base.Iterators: zip, product
using IterTools: imap
using LayerDicts
using Memento
using Missings
using OffsetArrays
using TimeZones
using Compat: Compat, @__MODULE__, AbstractDict, Cvoid, axes, @compat, lastindex, replace,
              something, undef
import Compat.Distributed: clear!

const Parameter = Union{String, Missing}
const LOGGER = getlogger(@__MODULE__)

function __init__()
    Memento.register(LOGGER)
end

# Docstring template for types using DocStringExtensions
@template TYPES =
    """
        $(TYPEDEF)

    $(DOCSTRING)

    ## Fields:

    $(FIELDS)
    """

include(joinpath(@__DIR__, "utils.jl"))

module libpq_c
    export Oid

    using Compat: Cvoid

    include(joinpath(@__DIR__, "..", "deps", "deps.jl"))

    function __init__()
        check_deps()
    end

    include(joinpath(@__DIR__, "headers", "libpq-fe.jl"))
end

using .libpq_c

include("typemaps.jl")

"""
    const LIBPQ_TYPE_MAP::PQTypeMap

The [`PQTypeMap`](@ref) containing LibPQ-level type mappings for LibPQ.jl.
Adding type mappings to this constant will override the default type mappings for all code
using LibPQ.jl.
"""
const LIBPQ_TYPE_MAP = PQTypeMap()

"""
    const LIBPQ_CONVERSIONS::PQConversions

The [`PQConversions`](@ref) containing LibPQ-level conversion functions for LibPQ.jl.
Adding conversions to this constant will override the default conversions for all code using
LibPQ.jl.
"""
const LIBPQ_CONVERSIONS = PQConversions()

### CONNECTIONS BEGIN
show_option(str::String) = string(replace(str, [' ', '\\'] => s -> "\\$s"))
show_option(bool::Bool) = ifelse(bool, 't', 'f')
show_option(num::Real) = num

# values containing spaces may not work correctly on PostgreSQL versions before 9.6
const CONNECTION_OPTION_DEFAULTS = Dict{String, String}(
    "DateStyle" => "ISO,YMD",
    "IntervalStyle" => "iso_8601",
    "TimeZone" => "UTC",
)

const CONNECTION_PARAMETER_DEFAULTS = Dict{String, String}(
    "client_encoding" => "UTF8",
    "application_name" => "LibPQ.jl",
    "options" => join(
        ("-c $k=$(show_option(v))" for (k, v) in CONNECTION_OPTION_DEFAULTS),
        " ",
    ),
)

"A connection to a PostgreSQL database."
mutable struct Connection
    "A pointer to a libpq PGconn object (C_NULL if closed)"
    conn::Ptr{libpq_c.PGconn}

    "libpq client encoding (string encoding of returned data)"
    encoding::String

    "Integer counter for generating connection-level unique identifiers"
    uid_counter::UInt

    "Connection-level type correspondence map"
    type_map::PQTypeMap

    "Connection-level conversion functions"
    func_map::PQConversions

    "True if the connection is closed and the PGconn object has been cleaned up"
    closed::Bool

    function Connection(
        conn::Ptr,
        closed=false;
        type_map::AbstractDict=PQTypeMap(),
        conversions::AbstractDict=PQConversions(),
    )
        return new(conn, "UTF8", 0, PQTypeMap(type_map), PQConversions(conversions), closed)
    end
end

"""
    handle_new_connection(jl_conn::Connection; throw_error=true) -> Connection

Check status and handle errors for newly-created connections.
Also set the client encoding ([23.3. Character Set Support](https://www.postgresql.org/docs/10/static/multibyte.html))
to `jl_conn.encoding`.

If `throw_error` is `true`, an error will be thrown if the connection's status is
`CONNECTION_BAD` and the PGconn object will be cleaned up.
Otherwise, a warning will be shown and the user should call `close` or `reset!` on the
returned `Connection`.
"""
function handle_new_connection(jl_conn::Connection; throw_error::Bool=true)
    if status(jl_conn) == libpq_c.CONNECTION_BAD
        err = error_message(jl_conn)

        if throw_error
            close(jl_conn)
            error(LOGGER, err)
        else
            warn(LOGGER, err)
        end
    else
        # if connection is successful, set client_encoding
        reset_encoding!(jl_conn)
    end

    return jl_conn
end

"""
    Connection(
        str::AbstractString;
        throw_error::Bool=true,
        type_map::AbstractDict=LibPQ.PQTypeMap(),
        conversions::AbstractDict=LibPQ.PQConversions(),
    ) -> Connection

Create a `Connection` from a connection string as specified in the PostgreSQL
documentation ([33.1.1. Connection Strings](https://www.postgresql.org/docs/10/static/libpq-connect.html#LIBPQ-CONNSTRING)).

For information on the `type_map` and `conversions` arguments, see [Type Conversions](@ref typeconv).

See [`handle_new_connection`](@ref) for information on the `throw_error` argument.
"""
function Connection(str::AbstractString; throw_error::Bool=true, kwargs...)
    ci_array = conninfo(str)

    keywords = String[]
    values = String[]

    for (k, v) in CONNECTION_PARAMETER_DEFAULTS
        push!(keywords, k)
        push!(values, v)
    end

    for ci in ci_array
        if !ismissing(ci.val)
            push!(keywords, ci.keyword)
            push!(values, ci.val)
        end
    end

    # Make the connection
    jl_conn = Connection(libpq_c.PQconnectdbParams(keywords, values, false); kwargs...)

    # If password needed and not entered, prompt the user
    if libpq_c.PQconnectionNeedsPassword(jl_conn.conn) == 1
        push!(keywords, "password")
        user = unsafe_string(libpq_c.PQuser(jl_conn.conn))
        prompt = "Enter password for PostgreSQL user $user:"
        push!(values, Base.getpass(prompt))
        return handle_new_connection(
            Connection(libpq_c.PQconnectdbParams(keywords, values, false); kwargs...);
            throw_error=throw_error,
        )
    else
        return handle_new_connection(
            jl_conn;
            throw_error=throw_error,
        )
    end

end

"""
    Connection(f, args...; kwargs...) -> Connection

A utility method to support `do` syntax.
Constructs the `Connection`, calls `f` on it, then closes it.
"""
function Connection(f::Base.Callable, args...; kwargs...)
    jl_conn = Connection(args...; kwargs...)

    try
        return f(jl_conn)
    finally
        close(jl_conn)
    end
end

"""
    server_version(jl_conn::Connection) -> VersionNumber

Get the PostgreSQL version of the server.

See [33.2. Connection Status Functions](https://www.postgresql.org/docs/10/static/libpq-status.html#LIBPQ-PQSERVERVERSION)
for information on the integer returned by `PQserverVersion` that is parsed by this
function.

See [`@pqv_str`](@ref) for information on how this packages represents PostgreSQL version
numbers.
"""
function server_version(jl_conn::Connection)
    version_int = libpq_c.PQserverVersion(jl_conn.conn)

    first_major = version_int ÷ 10000
    version = if first_major >= 10
        # new style (only major-minor)
        return VersionNumber(first_major, 0, version_int % 100)
    else
        # old style (major-major-minor)
        return VersionNumber(first_major, (version_int % 10000) ÷ 100, version_int % 100)
    end

    return version
end

"""
    @pqv_str -> VersionNumber

Parse a PostgreSQL version.

!!! note

    As of version 10.0, PostgreSQL moved from a three-part version number (where the first
    two parts represent the major version and the third represents the minor version) to a
    two-part major-minor version number.
    In LibPQ.jl, we represent this using the first two `VersionNumber` components as the
    major version and the third as the minor version.

    ## Examples

    ```jldoctest
    julia> using LibPQ: @pqv_str

    julia> pqv"10.1" == v"10.0.1"
    true

    julia> pqv"9.2.5" == v"9.2.5"
    true
    ```
"""
macro pqv_str(str)
    _pqv_str(str)
end

function _pqv_str(str)
    splitted = map(x -> parse(Int, x)::Int, split(str, '.'))

    if isempty(splitted)
        throw(ArgumentError("PostgreSQL version must contain at least one integer"))
    end

    version = if splitted[1] >= 10
        if length(splitted) == 1
            VersionNumber(splitted[1])
        elseif length(splitted) == 2
            VersionNumber(splitted[1], 0, splitted[2])
        else
            throw(ArgumentError(
                "PostgreSQL versions only have two components starting at version 10"
            ))
        end
    else
        if length(splitted) > 3
            throw(ArgumentError(
                "PostgreSQL versions cannot have more than three components"
            ))
        end

        VersionNumber(splitted...)
    end

    return version
end

"""
    encoding(jl_conn::Connection) -> String

Return the client encoding name for the current connection (see
[Table 23.1. PostgreSQL Character Sets](https://www.postgresql.org/docs/10/static/multibyte.html#CHARSET-TABLE)
for possible values).

Currently all Julia connections are set to use `UTF8` as this makes conversion to and from
`String` straighforward.

See also: [`set_encoding!`](@ref), [`reset_encoding!`](@ref)
"""
function encoding(jl_conn::Connection)
    encoding_id::Cint = libpq_c.PQclientEncoding(jl_conn.conn)

    if encoding_id == -1
        error(LOGGER, "libpq could not retrieve the connection's client encoding")
    end

    return unsafe_string(libpq_c.pg_encoding_to_char(encoding_id))
end

"""
    set_encoding!(jl_conn::Connection, encoding::String)

Set the client encoding for the current connection (see
[Table 23.1. PostgreSQL Character Sets](https://www.postgresql.org/docs/10/static/multibyte.html#CHARSET-TABLE)
for possible values).

Currently all Julia connections are set to use `UTF8` as this makes conversion to and from
`String` straighforward.
Other encodings are not explicitly handled by this package and will probably be very buggy.

See also: [`encoding`](@ref), [`reset_encoding!`](@ref)
"""
function set_encoding!(jl_conn::Connection, encoding::String)
    status = libpq_c.PQsetClientEncoding(jl_conn.conn, encoding)

    if status == -1
        error(LOGGER, "libpq could not set the connection's client encoding to $encoding")
    else
        jl_conn.encoding = encoding
    end

    return nothing
end

"""
    reset_encoding!(jl_conn::Connection, encoding::String)

Reset the client encoding for the current connection to `jl_conn.encoding`.

See also: [`encoding`](@ref), [`set_encoding!`](@ref)
"""
function reset_encoding!(jl_conn::Connection)
    set_encoding!(jl_conn, jl_conn.encoding)
end

"""
    unique_id(jl_conn::Connection, prefix::AbstractString="") -> String

Return a valid PostgreSQL identifier that is unique for the current connection.
This is mostly used to create names for prepared statements.
"""
function unique_id(jl_conn::Connection, prefix::AbstractString="")
    id_number, jl_conn.uid_counter = jl_conn.uid_counter, jl_conn.uid_counter + 1

    return "__libpq_$(prefix)_$(id_number)__"
end

"""
    status(jl_conn::Connection) -> libpq_c.ConnStatusType

Return the status of the PostgreSQL database connection according to libpq.
Only `CONNECTION_OK` and `CONNECTION_BAD` are valid for blocking connections, and only
blocking connections are supported right now.

See also: [`error_message`](@ref)
"""
status(jl_conn::Connection) = libpq_c.PQstatus(jl_conn.conn)

"""
    transaction_status(jl_conn::Connection) -> libpq_c.PGTransactionStatusType

Return the PostgreSQL database server's current in-transaction status for the connection.
See [](https://www.postgresql.org/docs/10/static/libpq-status.html#LIBPQ-PQTRANSACTIONSTATUS)
for information on the meaning of the possible return values.
"""
transaction_status(jl_conn::Connection) = libpq_c.PQtransactionStatus(jl_conn.conn)

"""
    close(jl_conn::Connection)

Close the PostgreSQL database connection and free the memory used by the `PGconn` object.
This function calls [`PQfinish`](https://www.postgresql.org/docs/10/static/libpq-connect.html#LIBPQ-PQFINISH),
but only if `jl_conn.closed` is `false`, to avoid a double-free.
"""
function Base.close(jl_conn::Connection)
    if !jl_conn.closed
        libpq_c.PQfinish(jl_conn.conn)
    end

    jl_conn.closed = true
    jl_conn.conn = C_NULL
    return nothing
end

"""
    isopen(jl_conn::Connection) -> Bool

Check whether a connection is open.
"""
Base.isopen(jl_conn::Connection) = !jl_conn.closed

"""
    reset!(jl_conn::Connection; throw_error=true)

Reset the communication to the PostgreSQL server.
The `PGconn` object will be recreated using identical connection parameters.

See [`handle_new_connection`](@ref) for information on the `throw_error` argument.

!!! note

    This function can be called on a connection with status `CONNECTION_BAD`, for example,
    but cannot be called on a connection that has been closed.
"""
function reset!(jl_conn::Connection; throw_error::Bool=true)
    if jl_conn.closed
        error(LOGGER, "Cannot reset a connection that has been closed")
    end

    libpq_c.PQreset(jl_conn.conn)
    handle_new_connection(jl_conn; throw_error=throw_error)

    return nothing
end

"""
    error_message(jl_conn::Connection) -> String

Return the error message most recently generated by an operation on the connection.
Includes a trailing newline.
"""
error_message(jl_conn::Connection) = unsafe_string(libpq_c.PQerrorMessage(jl_conn.conn))

"""
Indicator for how to display a PostgreSQL connection option (`PQconninfoOption`).

Possible values are:

* `Normal` (libpq: ""): display as is
* `Password` (libpq: "*"): hide the value of this field
* `Debug` (libpq: "D"): don't show by default
"""
@enum ConninfoDisplay Normal Password Debug

"""
    parse(::Type{ConninfoDisplay}, str::AbstractString) -> ConninfoDisplay

Parse a `ConninfoDisplay` from a string. See [`ConninfoDisplay`](@ref).
"""
function Base.parse(::Type{ConninfoDisplay}, str::AbstractString)::ConninfoDisplay
    if length(str) < 1
        Normal
    elseif first(str) == '*'
        Password
    elseif first(str) == 'D'
        Debug
    else
        error(LOGGER, "Unexpected dispchar '$str' in PQconninfoOption")
    end
end

"A Julia representation of a PostgreSQL connection option (`PQconninfoOption`)."
struct ConnectionOption
    "The name of the option"
    keyword::String

    "The name of the fallback environment variable for this option"
    envvar::Union{String, Missing}

    "The PostgreSQL compiled-in default for this option"
    compiled::Union{String, Missing}

    "The value of the option if set"
    val::Union{String, Missing}

    "The label of the option for display"
    label::String

    "Indicator for how to display the option (see [`ConninfoDisplay`](@ref))"
    disptype::ConninfoDisplay

    "The size of field to provide for entry of the option value (not used here)"
    dispsize::Int
end

"""
    ConnectionOption(pq_opt::libpq_c.PQconninfoOption) -> ConnectionOption

Construct a `ConnectionOption` from a `libpg_c.PQconninfoOption`.
"""
function ConnectionOption(pq_opt::libpq_c.PQconninfoOption)
    ConnectionOption(
        unsafe_string(pq_opt.keyword),
        unsafe_string_or_null(pq_opt.envvar),
        unsafe_string_or_null(pq_opt.compiled),
        unsafe_string_or_null(pq_opt.val),
        unsafe_string(pq_opt.label),
        parse(ConninfoDisplay, unsafe_string(pq_opt.dispchar)),
        pq_opt.dispsize,
    )
end

"""
    conninfo(jl_conn::Connection) -> Vector{ConnectionOption}

Get all connection options for a connection.
"""
function conninfo(jl_conn::Connection)
    ci_ptr = libpq_c.PQconninfo(jl_conn.conn)

    if ci_ptr == C_NULL
        error(LOGGER, "libpq could not allocate memory for connection info")
    end

    ci_array = conninfo(ci_ptr)
    libpq_c.PQconninfoFree(ci_ptr)
    return ci_array
end

function conninfo(ci_ptr::Ptr{libpq_c.PQconninfoOption})
    ci_array = Vector{ConnectionOption}()

    # ci_ptr is an array of PQconninfoOptions terminated by a PQconninfoOption with the
    # keyword field set to C_NULL
    ci_opt_idx = 1
    ci_opt = unsafe_load(ci_ptr, ci_opt_idx)
    while ci_opt.keyword != C_NULL
        push!(ci_array, ConnectionOption(ci_opt))

        ci_opt_idx += 1
        ci_opt = unsafe_load(ci_ptr, ci_opt_idx)
    end

    return ci_array
end

"""
    conninfo(str::AbstractString) -> Vector{ConnectionOption}

Parse connection options from a connection string (either a URI or key-value pairs).
"""
function conninfo(str::AbstractString)
    err_ref = Ref{Ptr{UInt8}}(C_NULL)
    ci_ptr = libpq_c.PQconninfoParse(str, err_ref)

    if ci_ptr == C_NULL && err_ref[] == C_NULL
        error(LOGGER, "libpq could not allocate memory for connection info")
    end

    if err_ref[] != C_NULL
        err_msg = unsafe_string(err_ref[])
        libpq_c.PQfreemem(err_ref[])
        error(err_msg)
    end

    ci_array = conninfo(ci_ptr)
    libpq_c.PQconninfoFree(ci_ptr)
    return ci_array
end

"""
    show(io::IO, jl_conn::Connection)

Display a [`Connection`](@ref) by showing the connection status and each connection option.
"""
function Base.show(io::IO, jl_conn::Connection)
    if jl_conn.closed
        print(io, "PostgreSQL connection (closed)")
        return nothing
    end

    print(io, "PostgreSQL connection ($(status(jl_conn))) with parameters:")
    for ci_opt in conninfo(jl_conn)
        if !ismissing(ci_opt.val) && ci_opt.disptype != Debug
            print(io, "\n  ", ci_opt.keyword, " = ")

            if ci_opt.disptype == Password
                print(io, "*" ^ ci_opt.dispsize)
            else
                print(io, ci_opt.val)
            end
        end
    end
end

### CONNECTIONS END

### RESULTS BEGIN

"A result from a PostgreSQL database query"
mutable struct Result <: Data.Source
    "A pointer to a libpq PGresult object (C_NULL if cleared)"
    result::Ptr{libpq_c.PGresult}

    "PostgreSQL Oids for each column in the result"
    column_oids::Vector{Oid}

    "Julia types for each column in the result"
    column_types::Vector{Type}

    "Whether to expect NULL for each column (whether output data can have `missing`)"
    not_null::Vector{Bool}

    "Conversions from PostgreSQL data to Julia types for each column in the result"
    column_funcs::Vector{Base.Callable}

    # TODO: attach encoding per https://wiki.postgresql.org/wiki/Driver_development#Result_object_and_client_encoding
    function Result(
        result::Ptr{libpq_c.PGresult},
        jl_conn::Connection;
        column_types::AbstractDict=ColumnTypeMap(),
        type_map::AbstractDict=PQTypeMap(),
        conversions::AbstractDict=PQConversions(),
        not_null=false,
    )
        jl_result = new(result)

        column_type_map = ColumnTypeMap()
        for (k, v) in column_types
            column_type_map[column_number(jl_result, k)] = v
        end

        type_lookup = LayerDict(
            PQTypeMap(type_map),
            jl_conn.type_map,
            LIBPQ_TYPE_MAP,
            _DEFAULT_TYPE_MAP,
        )

        func_lookup = LayerDict(
            PQConversions(conversions),
            jl_conn.func_map,
            LIBPQ_CONVERSIONS,
            _DEFAULT_CONVERSIONS,
            _FALLBACK_CONVERSION,
        )

        jl_result.column_oids = col_oids = map(1:num_columns(jl_result)) do col_num
            libpq_c.PQftype(jl_result.result, col_num - 1)
        end

        jl_result.column_types = col_types = collect(Type, imap(enumerate(col_oids)) do itr
            col_num, col_oid = itr
            get(column_type_map, col_num) do
                get(type_lookup, col_oid, String)
            end
        end)

        jl_result.column_funcs = collect(Base.Callable, imap(col_oids, col_types) do oid, typ
            func_lookup[(oid, typ)]
        end)

        # figure out which columns the user says may contain nulls
        if not_null isa Bool
            jl_result.not_null = fill(not_null, size(col_types))
        elseif not_null isa AbstractArray
            if eltype(not_null) === Bool
                if length(not_null) != length(col_types)
                    throw(ArgumentError(
                        "The length of keyword argument not_null, when an array, must be equal to the number of columns"
                    ))
                end

                jl_result.not_null = not_null
            else
                # assume array of column names
                jl_result.not_null = fill(false, size(col_types))

                for col_name in not_null
                    col_num = column_number(jl_result, col_name)
                    jl_result.not_null[col_num] = true
                end
            end
        else
            throw(ArgumentError(
                "Unsupported type $(typeof(not_null)) for keyword argument not_null"
            ))
        end

        # NOTE: The @compat annotation is necessary for 0.6 support, where the function
        # argument to `finalizer` comes last instead of first.
        @compat finalizer(close, jl_result)

        return jl_result
    end
end

"""
    show(io::IO, jl_result::Result)

Show a PostgreSQL result and whether it has been cleared.
"""
function Base.show(io::IO, jl_result::Result)
    print(io, "PostgreSQL result")

    if !isopen(jl_result)
        print(io, " (cleared)")
    end
end

"""
    status(jl_result::Result) -> libpq_c.ExecStatusType

Return the status of a result's corresponding database query according to libpq.
Only `CONNECTION_OK` and `CONNECTION_BAD` are valid for blocking connections, and only
blocking connections are supported right now.

See also: [`error_message`](@ref)
"""
status(jl_result::Result) = libpq_c.PQresultStatus(jl_result.result)

"""
    error_message(jl_result::Result) -> String

Return the error message associated with the result, or an empty string if there was no
error.
Includes a trailing newline.
"""
function error_message(jl_result::Result)
    unsafe_string(libpq_c.PQresultErrorMessage(jl_result.result))
end

"""
    close(jl_result::Result)

Clean up the memory used by the `PGresult` object.
The `Result` will no longer be usable.
"""
function Base.close(jl_result::Result)
    ptr, jl_result.result = jl_result.result, C_NULL
    if ptr != C_NULL
        libpq_c.PQclear(ptr)
    end
    return nothing
end

"""
    isopen(jl_result::Result)

Determine whether the given `Result` has been `close`d, i.e. whether the memory
associated with the underlying `PGresult` object has been cleared.
"""
Base.isopen(jl_result::Result) = jl_result.result != C_NULL

"""
    handle_result(jl_result::Result; throw_error::Bool=true) -> Result

Check status and handle errors for newly-created result objects.

If `throw_error` is `true`, throw an error and clear the result if the query results in a
fatal error or unreadable response.
Otherwise a warning is shown.

Also print an info message about the result.
"""
function handle_result(jl_result::Result; throw_error::Bool=true)
    err_msg = error_message(jl_result)
    result_status = status(jl_result)

    if result_status in (libpq_c.PGRES_BAD_RESPONSE, libpq_c.PGRES_FATAL_ERROR)
        if throw_error
            close(jl_result)
            error(LOGGER, err_msg)
        else
            warn(LOGGER, err_msg)
        end
    else
        if result_status == libpq_c.PGRES_NONFATAL_ERROR
            warn(LOGGER, err_msg)
        end

        info(LOGGER, unsafe_string(libpq_c.PQcmdStatus(jl_result.result)))
    end

    return jl_result
end

"""
    execute(
        {jl_conn::Connection, query::AbstractString | stmt::Statement},
        [parameters::AbstractVector,]
        throw_error::Bool=true,
        column_types::AbstractDict=ColumnTypeMap(),
        type_map::AbstractDict=LibPQ.PQTypeMap(),
        conversions::AbstractDict=LibPQ.PQConversions(),
    ) -> Result

Run a query on the PostgreSQL database and return a `Result`.
If `throw_error` is `true`, throw an error and clear the result if the query results in a
fatal error or unreadable response.

The query may be passed as `Connection` and `AbstractString` (SQL) arguments, or as a
`Statement`.

`execute` optionally takes a `parameters` vector which passes query parameters as strings to
PostgreSQL.

`column_types` accepts type overrides for columns in the result which take priority over
those in `type_map`.
For information on the `column_types`, `type_map`, and `conversions` arguments, see
[Type Conversions](@ref typeconv).
"""
function execute end

function execute(
    jl_conn::Connection,
    query::AbstractString;
    throw_error::Bool=true,
    kwargs...
)
    return handle_result(
        Result(libpq_c.PQexec(jl_conn.conn, query), jl_conn; kwargs...);
        throw_error=throw_error,
    )
end

function execute(
    jl_conn::Connection,
    query::AbstractString,
    parameters::AbstractVector;
    throw_error::Bool=true,
    kwargs...
)
    num_params = length(parameters)
    string_params = string_parameters(parameters)

    return handle_result(
        Result(libpq_c.PQexecParams(
            jl_conn.conn,
            query,
            num_params,
            C_NULL,  # set paramTypes to C_NULL to have the server infer a type
            parameter_pointers(string_params),
            C_NULL,  # paramLengths is ignored for text format parameters
            zeros(Cint, num_params),  # all parameters in text format
            zero(Cint),  # return result in text format
        ), jl_conn; kwargs...);
        throw_error=throw_error,
    )
end

"""
    string_parameters(parameters::AbstractVector) -> Vector{Union{String, Missing}}

Convert parameters to strings which can be passed to libpq, propagating `missing`.
"""
function string_parameters end

string_parameters(parameters::AbstractVector{<:Parameter}) = parameters

# vector which can't contain missing
string_parameters(parameters::AbstractVector) = map(string, parameters)

# vector which might contain missings
function string_parameters(parameters::AbstractVector{>:Missing})
    collect(
        Union{String, Missing},
        imap(parameters) do parameter
            ismissing(parameter) ? missing : string(parameter)
        end
    )
end

"""
    parameter_pointers(parameters::AbstractVector{<:Parameter}) -> Vector{Ptr{UInt8}}

Given a vector of parameters, returns a vector of pointers to either the string bytes in the
original or `C_NULL` if the element is `missing`.
"""
function parameter_pointers(parameters::AbstractVector{<:Parameter})
    pointers = Vector{Ptr{UInt8}}(undef, length(parameters))

    map!(pointers, parameters) do parameter
        ismissing(parameter) ? C_NULL : pointer(parameter)
    end

    return pointers
end

"""
    num_params(jl_result::Result) -> Int

Return the number of parameters in a prepared statement.
If this result did not come from the description of a prepared statement, return 0.
"""
function num_params(jl_result::Result)::Int
    # todo: check cleared?
    libpq_c.PQnparams(jl_result.result)
end

"""
    num_rows(jl_result::Result) -> Int

Return the number of rows in the query result.
This will be 0 if the query would never return data.
"""
function num_rows(jl_result::Result)::Int
    # todo: check cleared?
    libpq_c.PQntuples(jl_result.result)
end

"""
    num_affected_rows(jl_result::Result) -> Int

Return the number of rows affected by the command returning the result.
This is useful for counting the rows affected by operations such as INSERT,
UPDATE and DELETE that do not return rows but affect them.
This will be 0 if the query does not affect any row.
"""
function num_affected_rows(jl_result::Result)::Int
    # todo: check cleared?
    str = unsafe_string(libpq_c.PQcmdTuples(jl_result.result))
    if isempty(str)
        throw(ArgumentError("Result generated by an incompatible command"))
    else
        return parse(Int, str)
    end
end

"""
    num_columns(jl_result::Result) -> Int

Return the number of columns in the query result.
This will be 0 if the query would never return data.
"""
function num_columns(jl_result::Result)::Int
    # todo: check cleared?
    libpq_c.PQnfields(jl_result.result)
end

"""
    column_name(jl_result::Result, column_number::Integer) -> String

Return the name of the column at index `column_number` (1-based).
"""
function column_name(jl_result::Result, column_number::Integer)
    # todo: check cleared?
    unsafe_string(libpq_c.PQfname(jl_result.result, column_number - 1))
end

"""
    column_names(jl_result::Result) -> Vector{String}

Return the names of all the columns in the query result.
"""
function column_names(jl_result::Result)
    return [column_name(jl_result, i) for i in 1:num_columns(jl_result)]
end

"""
    column_number(jl_result::Result, column_name::Union{AbstractString, Symbol}) -> Int

Return the index (1-based) of the column named `column_name`.
"""
function column_number(jl_result::Result, column_name::Union{AbstractString, Symbol})::Int
    # todo: check cleared?
    return libpq_c.PQfnumber(jl_result.result, String(column_name)) + 1
end

"""
    column_number(jl_result::Result, column_idx::Integer) -> Int

Return the index of the column if it is valid, or error.
"""
function column_number(jl_result::Result, column_idx::Integer)::Int
    if !checkindex(Bool, 1:num_columns(jl_result), column_idx)
        throw(BoundsError(column_names(jl_result), column_idx))
    end

    return column_idx
end

"""
    column_oids(jl_result::Result) -> Vector{LibPQ.Oid}

Return the PostgreSQL oids for each column in the result.
"""
column_oids(jl_result::Result) = jl_result.column_oids

"""
    column_types(jl_result::Result) -> Vector{Type}

Return the corresponding Julia types for each column in the result.
"""
column_types(jl_result::Result) = jl_result.column_types

### RESULTS END

### PREPARE BEGIN

"A PostgreSQL prepared statement"
struct Statement
    """
    A `Connection` for which this statement is valid.
    It may become invalid if the connection is reset.
    """
    jl_conn::Connection

    "An autogenerated neame for the prepared statement (using [`unique_id`](@ref)"
    name::String

    "The query string of the prepared statement"
    query::String

    "A `Result` containing a description of the prepared statement"
    description::Result

    "The number of parameters accepted by this statement according to `description`"
    num_params::Int
end

"""
    prepare(jl_conn::Connection, query::AbstractString) -> Statement

Create a prepared statement on the PostgreSQL server using libpq.
The statement is given an generated unique name using [`unique_id`](@ref).

!!! note

    Currently the statement is not explicitly deallocated, but it is deallocated at the end
    of session per the [PostgreSQL documentation on DEALLOCATE](https://www.postgresql.org/docs/10/static/sql-deallocate.html).
"""
function prepare(jl_conn::Connection, query::AbstractString)
    uid = unique_id(jl_conn, "stmt")

    jl_result = handle_result(
        Result(libpq_c.PQprepare(
            jl_conn.conn,
            uid,
            query,
            0,  # infer all parameters from the query string
            C_NULL,
        ), jl_conn);
        throw_error=true,
    )

    close(jl_result)

    description = handle_result(
        Result(libpq_c.PQdescribePrepared(
            jl_conn.conn,
            uid,
        ), jl_conn);
        throw_error=true,
    )

    Statement(jl_conn, uid, query, description, num_params(description))
end

"""
    show(io::IO, jl_result::Statement)

Show a PostgreSQL prepared statement and its query.
"""
function Base.show(io::IO, stmt::Statement)
    print(
        io,
        "PostgreSQL prepared statement named ",
        stmt.name,
        " with query ",
        stmt.query,
    )
end

"""
    num_params(stmt::Statement) -> Int

Return the number of parameters in the prepared statement.
"""
num_params(stmt::Statement) = num_params(stmt.description)

"""
    num_columns(stmt::Statement) -> Int

Return the number of columns that would be returned by executing the prepared statement.
"""
num_columns(stmt::Statement) = num_columns(stmt.description)

"""
    column_name(stmt::Statement, column_number::Integer) -> String

Return the name of the column at index `column_number` (1-based) that would be returned by
executing the prepared statement.
"""
function column_name(stmt::Statement, column_number::Integer)
    column_name(stmt.description, column_number)
end

"""
    column_names(stmt::Statement) -> Vector{String}

Return the names of all the columns in the query result that would be returned by executing
the prepared statement.
"""
column_names(stmt::Statement) = column_names(stmt.description)

"""
    column_number(stmt::Statement, column_name::AbstractString) -> Int

Return the index (1-based) of the column named `column_name` that would be returned by
executing the prepared statement.
"""
function column_number(stmt::Statement, column_name::AbstractString)
    column_number(stmt.description, column_name)
end

function execute(
    stmt::Statement,
    parameters::AbstractVector;
    throw_error::Bool=true,
    kwargs...
)
    num_params = length(parameters)
    string_params = string_parameters(parameters)

    return handle_result(
        Result(libpq_c.PQexecPrepared(
            stmt.jl_conn.conn,
            stmt.name,
            num_params,
            parameter_pointers(string_params),
            C_NULL,  # paramLengths is ignored for text format parameters
            zeros(Cint, num_params),  # all parameters in text format
            zero(Cint),  # return result in text format
        ), stmt.jl_conn; kwargs...);
        throw_error=throw_error,
    )
end

function execute(
    stmt::Statement;
    throw_error::Bool=true,
    kwargs...
)
    return handle_result(
        Result(libpq_c.PQexecPrepared(
            stmt.jl_conn.conn,
            stmt.name,
            0,  # no parameters
            C_NULL,
            C_NULL,  # paramLengths is ignored for text format parameters
            C_NULL,  # all parameters in text format
            zero(Cint),  # return result in text format
        ), stmt.jl_conn; kwargs...);
        throw_error=throw_error,
    )
end

### PREPARE END

include("parsing.jl")
include("datastreams.jl")

Base.@deprecate clear!(jl_result::Result) close(jl_result)

end


text: include("string.jl")