text:
###############################################################################
#
# Calculate empirical gsw from the equations
#
###############################################################################
"""
    stomatal_conductance(
                model::EmpiricalStomatalModel{FT},
                leaf::Leaf{FT},
                envir::AirLayer{FT},
                β::FT
    ) where {FT<:AbstractFloat}
    stomatal_conductance(
                model::EmpiricalStomatalModel{FT},
                canopyi::CanopyLayer{FT},
                envir::AirLayer{FT},
                β::FT
    ) where {FT<:AbstractFloat}
    stomatal_conductance(
                model::EmpiricalStomatalModel{FT},
                canopyi::CanopyLayer{FT},
                envir::AirLayer{FT},
                β::FT,
                ind::Int
    ) where {FT<:AbstractFloat}

Steady state gsw from empirical approach given
- `model` [`EmpiricalStomatalModel`](@ref) type empirical model parameter set
- `leaf` [`Leaf`] type struct
- `canopyi` [`CanopyLayer`](@ref) type struct
- `envir` [`AirLayer`] type struct
- `β` Correction factor over the g1 part of an empirical model
- `ind` Nth leaf in the canopy layer
"""
function stomatal_conductance(
            model::ESMBallBerry{FT},
            leaf::Leaf{FT},
            envir::AirLayer{FT},
            β::FT
) where {FT<:AbstractFloat}
    @unpack g0, g1    = model;
    @unpack An, p_s   = leaf;
    @unpack p_atm, RH = envir;

    return g0 + g1 * RH * p_atm * FT(1e-6) * β * An / p_s
end

function stomatal_conductance(
            model::ESMBallBerry{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT},
            β::FT
) where {FT<:AbstractFloat}
    @unpack g0, g1    = model;
    @unpack An, p_s   = canopyi;
    @unpack p_atm, RH = envir;

    return g0 .+ g1 * RH * p_atm * FT(1e-6) * β .* An ./ p_s
end

function stomatal_conductance(
            model::ESMBallBerry{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT},
            β::FT,
            ind::Int
) where {FT<:AbstractFloat}
    @unpack g0, g1    = model;
    @unpack An, p_s   = canopyi;
    @unpack p_atm, RH = envir;

    return g0 + g1 * RH * p_atm * FT(1e-6) * β * An[ind] / p_s[ind]
end

function stomatal_conductance(
            model::ESMGentine{FT},
            leaf::Leaf{FT},
            envir::AirLayer{FT},
            β::FT
) where {FT<:AbstractFloat}
    @unpack g0, g1  = model;
    @unpack An, p_i = leaf;
    @unpack p_atm   = envir;

    return g0 + g1 * p_atm * FT(1e-6) * β * An / p_i
end

function stomatal_conductance(
            model::ESMGentine{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT},
            β::FT
) where {FT<:AbstractFloat}
    @unpack g0, g1  = model;
    @unpack An, p_i = canopyi;
    @unpack p_atm   = envir;

    return g0 .+ g1 * p_atm * FT(1e-6) * β .* An ./ p_i
end

function stomatal_conductance(
            model::ESMGentine{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT},
            β::FT,
            ind::Int
) where {FT<:AbstractFloat}
    @unpack g0, g1  = model;
    @unpack An, p_i = canopyi;
    @unpack p_atm   = envir;

    return g0 + g1 * p_atm * FT(1e-6) * β * An[ind] / p_i[ind]
end

function stomatal_conductance(
            model::ESMLeuning{FT},
            leaf::Leaf{FT},
            envir::AirLayer{FT},
            β::FT
) where {FT<:AbstractFloat}
    @unpack d0, g0, g1             = model;
    @unpack An, p_s, p_sat, Γ_star = leaf;
    @unpack p_atm, p_H₂O           = envir;

    return g0 + g1 * p_atm * FT(1e-6) / (1 + (p_sat - p_H₂O)/d0) *
                β * An / (p_s - Γ_star)
end

function stomatal_conductance(
            model::ESMLeuning{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT},
            β::FT
) where {FT<:AbstractFloat}
    @unpack d0, g0, g1     = model;
    @unpack An, p_s, p_sat = canopyi;
    @unpack Γ_star         = canopyi.ps;
    @unpack p_atm, p_H₂O   = envir;

    return g0 .+ g1 * p_atm * FT(1e-6) / (1 + (p_sat - p_H₂O)/d0) *
                 β .* An ./ (p_s .- Γ_star)
end

function stomatal_conductance(
            model::ESMLeuning{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT},
            β::FT,
            ind::Int
) where {FT<:AbstractFloat}
    @unpack d0, g0, g1     = model;
    @unpack An, p_s, p_sat = canopyi;
    @unpack Γ_star         = canopyi.ps;
    @unpack p_atm, p_H₂O   = envir;

    return g0 + g1 * p_atm * FT(1e-6) / (1 + (p_sat - p_H₂O)/d0) *
                β * An[ind] / (p_s[ind] - Γ_star)
end

function stomatal_conductance(
            model::ESMMedlyn{FT},
            leaf::Leaf{FT},
            envir::AirLayer{FT},
            β::FT
) where {FT<:AbstractFloat}
    @unpack g0, g1            = model;
    @unpack An, p_sat         = leaf;
    @unpack p_a, p_atm, p_H₂O = envir;
    vpd = max(FT(0.001), p_sat - p_H₂O);

    return g0 + p_atm * FT(1e-6) / p_a * (1 + g1/sqrt(vpd)) * β * An
end

function stomatal_conductance(
            model::ESMMedlyn{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT},
            β::FT
) where {FT<:AbstractFloat}
    @unpack g0, g1            = model;
    @unpack An, p_sat         = canopyi;
    @unpack p_a, p_atm, p_H₂O = envir;
    vpd = max(FT(0.001), p_sat - p_H₂O);

    return g0 .+ p_atm * FT(1e-6) / p_a * (1 + g1/sqrt(vpd)) * β .* An
end

function stomatal_conductance(
            model::ESMMedlyn{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT},
            β::FT,
            ind::Int
) where {FT<:AbstractFloat}
    @unpack g0, g1            = model;
    @unpack An, p_sat         = canopyi;
    @unpack p_a, p_atm, p_H₂O = envir;
    vpd = max(FT(0.001), p_sat - p_H₂O);

    return g0 + p_atm * FT(1e-6) / p_a * (1 + g1/sqrt(vpd)) * β * An[ind]
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using PowerModels
using .PP2PM

function run_powermodels(json_path)
    pm = PP2PM.load_pm_from_json(json_path)

    model = PP2PM.get_model(pm["pm_model"])

    solver = PP2PM.get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"],
    pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"])

    result = PowerModels.run_pf(pm, model, solver)
    # add branch flows
    update_data!(pm, result["solution"])
    flows = calc_branch_flow_ac(pm)
    update_data!(result["solution"], flows)
    return result
end

ext: jl

max_issues_repo_name: shravanchipli/pandapower

max_forks_repo_name: gdgarcia/pandapower


text:
@doc doc"""
    dimer_lattice()
    dimer_lattice(param::Dict)

Generates dimer lattice (indeed, this is not a "lattice")
"""
function dimer_lattice()
    dim = 1
    sz = [2]
    nsitetypes = 2
    nbondtypes = 1
    sites = Vector{Int}[zeros(Int,0) for i in 1:nsitetypes]
    bonds = Vector{Int}[zeros(Int,0) for i in 1:nbondtypes]
    nsites = 2
    nbonds = 1
    sitetypes = [1,2]
    bondtypes = [1]
    transvector = @compat Matrix(1.0I, dim, dim)
    coords = zeros(dim, 2)
    bond_dirs = zeros(dim, 1)
    neighborsites = [[2],[1]]
    neighborbonds = [[1],[1]]
    source = zeros(Int,1)
    target = zeros(Int,1)
    site_L2 = zeros(Int,2)
    site_L4 = zeros(Int,2)

    sites[1] = [1]
    sites[2] = [2]
    bonds[1] = [1]
    coords[1,1] = 0.0
    coords[1,2] = 1.0
    bond_dirs[1,1] = 1.0
    source[1] = 1
    target[1] = 2
    site_L2[1] = 2
    site_L2[2] = 1
    site_L4[1] = 2
    site_L4[2] = 1

    Lattice(dim, sz, nsitetypes, nbondtypes,
            sites, bonds, nsites, nbonds, sitetypes, bondtypes,
            transvector, coords, bond_dirs,
            neighborsites, neighborbonds,
            source,target, site_L2, site_L4)
end
dimer_lattice(param::Dict) = dimer_lattice()

@doc doc"""
    chain_lattice(L::Integer)
    chain_lattice(param::Dict)

Generates chain lattice with length `L = param["L"]`

`nsitetypes` is 2: all `1` (`2`) sites do not connect to another `1` (`2`) site.

`nbondtypes` is 2: `1` bond connects `2n-1` and `2n` sites and `2` bond connects `2n` and `2n+1` sites.
"""
function chain_lattice(L::Integer)
    dim = 1
    sz = [L]
    nsitetypes = 2
    nbondtypes = 2
    nsites = L
    nbonds = L
    sitetypes = zeros(Int,L)
    bondtypes = zeros(Int,L)
    transvector = @compat Matrix(1.0I, dim, dim)
    coords = zeros(dim, L)
    bond_dirs = zeros(dim, L)
    neighborsites = Vector{Int}[]
    neighborbonds = Vector{Int}[]
    source = zeros(Int,L)
    target = zeros(Int,L)
    site_L2 = zeros(Int,L)
    site_L4 = zeros(Int,L)

    ns = zeros(Int,2)
    nb = zeros(Int,2)
    for s in 1:L
        coords[1,s] = (s-1)
        bond_dirs[1,s] = 1.0
        ns[1] = mod1(s+1,L)
        ns[2] = mod1(s-1,L)
        nb[1] = s
        nb[2] = mod1(s-1,L)
        push!(neighborsites, copy(ns))
        push!(neighborbonds, copy(nb))
        source[s] = s
        target[s] = neighborsites[s][1]
        sitetypes[s] = ifelse(isodd(s),1,2)
        bondtypes[s] = ifelse(isodd(s),1,2)
        site_L2[s] = mod1(s+div(L,2),L)
        site_L4[s] = mod1(s+div(L,4),L)
    end
    sites = Vector{Int}[zeros(Int,0) for i in 1:nsitetypes]
    bonds = Vector{Int}[zeros(Int,0) for i in 1:nbondtypes]
    for s in 1:nsites
        push!(sites[sitetypes[s]], s)
    end
    for b in 1:nbonds
        push!(bonds[bondtypes[b]], b)
    end
    Lattice(dim, sz, nsitetypes, nbondtypes,
            sites, bonds, nsites, nbonds, sitetypes, bondtypes,
            transvector, coords, bond_dirs, neighborsites, neighborbonds,
            source,target, site_L2, site_L4)
end
chain_lattice(param::Dict) = chain_lattice(param["L"])

@doc doc"""
    square_lattice(L::Integer, W::Integer=L)
    square_lattice(param::Dict)

Generates square lattice with size `L=param["L"]` $\times$ `W=param["W"]`.

`nsitetypes` is 2: all `1` (`2`) sites do not connect to another `1` (`2`) site.

`nbondtypes` is 2: `1` bonds are parallel to `x` axis and `2` are parallel to `y` axis.
"""
square_lattice(L::Integer) = square_lattice(L,L)
function square_lattice(L::Integer, W::Integer)
    s2xy(s::Integer) = mod(s-1,L),div(s-1,L)
    xy2s(x::Integer, y::Integer) = mod(y,W)*L + mod(x,L) + 1

    dim = 2
    sz = [L,W]
    nsitetypes = 2
    nbondtypes = 2
    nsites = L*W
    nbonds = 2*L*W
    sitetypes = zeros(Int, nsites)
    bondtypes = zeros(Int, nbonds)
    transvector = @compat Matrix(1.0I, dim, dim)
    coords = zeros(dim, nsites)
    bond_dirs = zeros(dim, nbonds)
    neighborsites = Vector{Int}[]
    neighborbonds = Vector{Int}[]
    source = zeros(Int,nbonds)
    target = zeros(Int,nbonds)
    sitetypes = zeros(Int,nsites)
    bondtypes = zeros(Int,nbonds)
    site_L2 = zeros(Int,nsites)
    site_L4 = zeros(Int,nsites)
    L2 = div(L,2)
    L4 = div(L,4)
    W2 = div(W,2)
    W4 = div(W,4)
    ns = zeros(Int,4)
    nb = zeros(Int,4)
    @inbounds for s in 1:nsites
        x,y = s2xy(s)
        ns[1] = xy2s(x+1,y)
        ns[2] = xy2s(x,y+1)
        ns[3] = xy2s(x-1,y)
        ns[4] = xy2s(x,y-1)
        push!(neighborsites, copy(ns))
        nb[1] = 2s-1
        nb[2] = 2s
        nb[3] = 2ns[3]-1
        nb[4] = 2ns[4]
        push!(neighborbonds, copy(nb))
        source[2s-1:2s] .= s
        target[2s-1:2s] .= neighborsites[s][1:2]
        coords[1,s] = x
        coords[2,s] = y
        sitetypes[s] = ifelse(iseven(x+y),1,2)
        bond_dirs[:, 2s-1] = [1.0, 0.0]
        bond_dirs[:, 2s] = [0.0, 1.0]
        bondtypes[2s-1] = 1
        bondtypes[2s] = 2
        site_L2[s] = xy2s(x+L2, y+W2)
        site_L4[s] = xy2s(x+L4, y+W4)
    end
    sites = Vector{Int}[zeros(Int,0) for i in 1:nsitetypes]
    bonds = Vector{Int}[zeros(Int,0) for i in 1:nbondtypes]
    for s in 1:nsites
        push!(sites[sitetypes[s]], s)
    end
    for b in 1:nbonds
        push!(bonds[bondtypes[b]], b)
    end
    Lattice(dim, sz, nsitetypes, nbondtypes,
            sites, bonds, nsites, nbonds, sitetypes, bondtypes,
            transvector, coords, bond_dirs, neighborsites, neighborbonds,
            source,target, site_L2, site_L4)
end
function square_lattice(param::Dict)
    L = param["L"]
    W = get(param, "W", L)
    return square_lattice(L, W)
end

@doc doc"""
    triangular_lattice(L::Integer, W::Integer=L)
    triangular_lattice(param::Dict)

Generates triangular lattice with size `L=param["L"]` $\times$ `W=param["W"]`.

`nsitetypes` is 3: all `1`, `2`, and `3` sites do not connect to another `1`, `2`, and `3` site, respectively.

`nbondtypes` is 3: `1`, `2`, and `3` bonds make an angle of 0, 60, and 120 degree with `x` axis, respectively.
"""
triangular_lattice(L::Integer) = triangular_lattice(L,L)
function triangular_lattice(L::Integer, W::Integer)
    s2xy(s::Integer) = mod(s-1,L),div(s-1,L)
    xy2s(x::Integer, y::Integer) = mod(y,W)*L + mod(x,L) + 1
    ex = [1.0, 0.0]
    ey = [0.5, 0.5sqrt(3)]

    dim = 2
    sz = [L,W]
    nsitetypes = 3
    nbondtypes = 3
    nsites = L*W
    nbonds = 3nsites
    sitetypes = zeros(Int,nsites)
    bondtypes = zeros(Int,nbonds)
    transvector = [1.0 0.5; 0.0 0.5sqrt(3)]
    coords = zeros(dim, nsites)
    bond_dirs = zeros(dim, nbonds)
    neighborsites = Vector{Int}[]
    neighborbonds = Vector{Int}[]
    source = zeros(Int,nbonds)
    target = zeros(Int,nbonds)
    site_L2 = zeros(Int, nsites)
    site_L4 = zeros(Int, nsites)
    L2 = div(L,2)
    L4 = div(L,4)
    W2 = div(W,2)
    W4 = div(W,4)
    ns = zeros(Int,6)
    nb = zeros(Int,6)
    @inbounds for s in 1:nsites
        x,y = s2xy(s)
        ns[1] = xy2s(x+1,y)
        ns[2] = xy2s(x,y+1)
        ns[3] = xy2s(x-1,y+1)
        ns[4] = xy2s(x-1,y)
        ns[5] = xy2s(x,y-1)
        ns[6] = xy2s(x+1,y-1)
        push!(neighborsites, copy(ns))
        nb[1] = 3s-2
        nb[2] = 3s-1
        nb[3] = 3s
        nb[4] = 3ns[4]-2
        nb[5] = 3ns[5]-1
        nb[6] = 3ns[6]
        push!(neighborbonds, copy(nb))
        source[3s-2:3s] .= s
        target[3s-2:3s] .= neighborsites[s][1:3]
        coords[1,s] = x
        coords[2,s] = y
        if sitetypes[s]==0
            sitetypes[s] = 1
        else
            if sitetypes[neighbors[1,s]]==0
                sitetypes[neighbors[1,s]]=mod1(sitetypes[s]+1,3)
            end
            if sitetypes[neighbors[2,s]]==0
                sitetypes[neighbors[2,s]]=mod1(sitetypes[s]+2,3)
            end
        end
        bond_dirs[:, 3s-2] .= [1.0,0.0]
        bond_dirs[:, 3s-1] .= [0.0,1.0]
        bond_dirs[:, 3s] .= [-1.0,1.0]
        bondtypes[3s-2] = 1
        bondtypes[3s-1] = 2
        bondtypes[3s] = 3
        site_L2[s] = xy2s(x+L2, y+W2)
        site_L4[s] = xy2s(x+L4, y+W4)
    end
    sites = Vector{Int}[zeros(Int,0) for i in 1:nsitetypes]
    bonds = Vector{Int}[zeros(Int,0) for i in 1:nbondtypes]
    for s in 1:nsites
        push!(sites[sitetypes[s]], s)
    end
    for b in 1:nbonds
        push!(bonds[bondtypes[b]], b)
    end
    Lattice(dim, sz, nsitetypes, nbondtypes,
            sites, bonds, nsites, nbonds, sitetypes, bondtypes,
            transvector, coords, bond_dirs, neighborsites, neighborbonds,
            source,target, site_L2, site_L4)
end
function triangular_lattice(param::Dict)
    L = param["L"]
    W = get(param, "W", L)
    return triangular_lattice(L, W)
end

@doc doc"""
    cubic_lattice(L::Integer, W::Integer=L, H::Integer=W)
    cubic_lattice(param::Dict)

Generates cubic lattice with size `L=param["L"]` $\times$ `W=param["W"]` $\times$ `H=param["H"]`.

`nsitetypes` is 2: all `1` (`2`) sites do not connect to another `1` (`2`) site.

`nbondtypes` is 2: `1`, `2`, and `3` bonds are parallel to `x`, `y`, and `z` axis, respectively.
"""
cubic_lattice(L::Integer) = cubic_lattice(L,L,L)
cubic_lattice(L::Integer, W::Integer) = cubic_lattice(L,W,W)
function cubic_lattice(L::Integer, W::Integer, H::Integer)
    function s2xyz(s::Integer)
        z = div(s-1, L*W)
        xy = mod(s-1, L*W)
        x = mod(xy, L)
        y = div(xy, L)
        return x, y, z
    end
    xyz2s(x::Integer, y::Integer, z::Integer) = mod(z,H)*W*L + mod(y,W)*L + mod(x,L) + 1

    dim = 3
    sz = [L,W,H]
    nsitetypes = 2
    nbondtypes = 3
    nsites = L*W*H
    nbonds = 3nsites
    sitetypes = zeros(Int,nsites)
    bondtypes = zeros(Int,nbonds)
    transvector = @compat Matrix(1.0I, dim, dim)
    coords = zeros(dim, nsites)
    bond_dirs = zeros(dim, nbonds)
    neighborsites = Vector{Int}[]
    neighborbonds = Vector{Int}[]
    source = zeros(Int,nbonds)
    target = zeros(Int,nbonds)
    site_L2 = zeros(Int,nsites)
    site_L4 = zeros(Int,nsites)
    L2 = div(L,2)
    L4 = div(L,4)
    W2 = div(W,2)
    W4 = div(W,4)
    H2 = div(H,2)
    H4 = div(H,4)
    ns = zeros(Int,6)
    nb = zeros(Int,6)
    @inbounds for s in 1:nsites
        x,y,z = s2xyz(s)
        ns[1] = xyz2s(x+1,y,z)
        ns[2] = xyz2s(x,y+1,z)
        ns[3] = xyz2s(x,y,z+1)
        ns[4] = xyz2s(x-1,y,z)
        ns[5] = xyz2s(x,y-1,z)
        ns[6] = xyz2s(x,y,z-1)
        push!(neighborsites, copy(ns))
        nb[1] = 3s-2
        nb[2] = 3s-1
        nb[3] = 3s
        nb[4] = 3ns[4]-2
        nb[5] = 3ns[5]-1
        nb[6] = 3ns[6]
        push!(neighborbonds, copy(nb))
        source[3s-2:3s] .= s
        target[3s-2:3s] .= neighborsites[s][1:3]
        coords[1,s] = x
        coords[2,s] = y
        coords[3,s] = z
        sitetypes[s] = ifelse(iseven(x+y+z),1,2)
        bondtypes[3s-2] = 1
        bondtypes[3s-1] = 2
        bondtypes[3s] = 3
        bond_dirs[:, 3s-2] = [1.0, 0.0, 0.0]
        bond_dirs[:, 3s-1] = [0.0, 1.0, 0.0]
        bond_dirs[:, 3s-0] = [0.0, 0.0, 1.0]
        site_L2[s] = xyz2s(x+L2,y+W2,z+H2)
        site_L4[s] = xyz2s(x+L4,y+W4,z+H4)
    end
    sites = Vector{Int}[zeros(Int,0) for i in 1:nsitetypes]
    bonds = Vector{Int}[zeros(Int,0) for i in 1:nbondtypes]
    for s in 1:nsites
        push!(sites[sitetypes[s]], s)
    end
    for b in 1:nbonds
        push!(bonds[bondtypes[b]], b)
    end
    Lattice(dim, sz, nsitetypes, nbondtypes,
            sites, bonds, nsites, nbonds, sitetypes, bondtypes,
            transvector, coords, bond_dirs,
            neighborsites, neighborbonds,
            source,target, site_L2, site_L4)
end
function cubic_lattice(param::Dict)
    L = param["L"]
    W = get(param, "W", L)
    H = get(param, "H", W)
    return cubic_lattice(L, W, H)
end

@doc doc"""
    fully_connected_lattice(N::Integer)
    fully_connected_lattice(param::Dict)

Generates `N=param["N"]` site fully connected lattice

Both `nsitetypes` and `nbondtypes` are 1.
"""
function fully_connected_lattice(N::Integer)
    dim = 1
    nbonds = div(N*(N-1),2)
    nsitetypes = 1
    nbondtypes = 1
    sites = [collect(1:N)]
    bonds = [collect(1:nbonds)]
    sitetypes = ones(N)
    bondtypes = ones(nbonds)
    neighborsites = [vcat(1:s-1, s+1:N) for s in 1:N]
    neighborbonds = Vector{Int}[zeros(Int,0) for s in 1:N]
    source = zeros(Int,nbonds)
    target = zeros(Int,nbonds)
    ib = 0
    for i in 1:N
        for j in (i+1):N
            ib += 1
            source[ib] = i
            target[ib] = j
            push!(neighborbonds[i], ib)
            push!(neighborbonds[j], ib)
        end
    end

    lat = Lattice(dim, [N], nsitetypes, nbondtypes, sites, bonds, N, nbonds, sitetypes, bondtypes,
            zeros(1,1), zeros(1,N), zeros(1,nbonds), neighborsites, neighborbonds, source, target, zeros(Int,N), zeros(Int,N))
    return lat
end
fully_connected_lattice(param::Dict) = fully_connected_lattice(param["N"])


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	load_inputs(setup::Dict,path::AbstractString)

Loads various data inputs from multiple input .csv files in path directory and stores variables in a Dict (dictionary) object for use in model() function

inputs:
setup - dict object containing setup parameters
path - string path to working directory

returns: Dict (dictionary) object containing all data inputs
"""
function load_inputs(setup::Dict,path::AbstractString)

	data_directory = chop(replace(path, pwd() => ""), head = 1, tail = 0)

	## Read input files
	println("Reading Input CSV Files")
	## Declare Dict (dictionary) object used to store parameters
	inputs = Dict()
	# Read input data about power network topology, operating and expansion attributes
    if isfile(joinpath(path,"Network.csv"))
		inputs, network_var = load_network_data(setup, path, inputs)
	else
		inputs["Z"] = 1
		inputs["L"] = 0
	end

	# Read temporal-resolved load data, and clustering information if relevant
	inputs = load_load_data(setup, path, inputs)
	# Read fuel cost data, including time-varying fuel costs
	inputs, cost_fuel, CO2_fuel = load_fuels_data(setup, path, inputs)
	# Read in generator/resource related inputs
	inputs = load_generators_data(setup, path, inputs, cost_fuel, CO2_fuel)
	# Read in generator/resource availability profiles
	inputs = load_generators_variability(setup, path, inputs)

	if setup["CapacityReserveMargin"]==1
		inputs = load_cap_reserve_margin(setup, path, inputs)
		if inputs["Z"] >1
			inputs = load_cap_reserve_margin_trans(setup, inputs,network_var)
		end
	end

	# Read in general configuration parameters for reserves (resource-specific reserve parameters are read in generators_data())
	if setup["Reserves"]==1
		inputs = load_reserves(setup, path, inputs)
	end

	if setup["MinCapReq"] == 1
		inputs = load_minimum_capacity_requirement(path, inputs, setup)
	end

	if setup["EnergyShareRequirement"]==1
		inputs = load_energy_share_requirement(setup, path, inputs)
	end

	if setup["CO2Cap"] >= 1
		inputs = load_co2_cap(setup, path, inputs)
	end

	# Read in mapping of modeled periods to representative periods
	if setup["OperationWrapping"]==1 && !isempty(inputs["STOR_LONG_DURATION"]) && (isfile(data_directory*"/Period_map.csv") || isfile(joinpath(data_directory,joinpath(setup["TimeDomainReductionFolder"],"Period_map.csv")))) # Use Time Domain Reduced data for GenX)
		inputs = load_period_map(setup, path, inputs)
	end

	println("CSV Files Successfully Read In From $path")

	return inputs
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
open("crud.txt", "w") do io
    write(io, "Hello world!")
end;

open(f->read(f, String), "crud.txt")

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming


text:
"""
Compare an on-disk table to a table schema.
"""
module ondisk_table

using CSV
using DataFrames

using ..types
using ..common

function compare(tableschema::TableSchema, input_data_file::String, output_data_file::String,
                 input_issues_file::String, output_issues_file::String, sorted_by_primarykey::Bool)
    # Init
    tablename     = tableschema.name
    outdata       = init_outdata(tableschema, input_data_file)
    issues_in     = init_issues(tableschema)  # Issues for indata
    issues_out    = init_issues(tableschema)  # Issues for outdata
    pk_colnames   = tableschema.primarykey
    pk_ncols      = length(pk_colnames)
    primarykey    = missings(String, pk_ncols)  # Stringified primary key
    pkvalues_in   = Set{String}()  # Values of the primary key
    pk_issues_in  = pk_ncols == 1 ? issues_in[:columnissues][pk_colnames[1]]  : Dict{Symbol, Int}()
    i_outdata     = 0
    nconstraints  = length(tableschema.intrarow_constraints)
    pk_colnames_set   = Set(pk_colnames)
    colname2colschema = tableschema.colname2colschema
    uniquevalues_in   = Dict(colname => Set{colschema.datatype}() for (colname, colschema) in colname2colschema if colschema.isunique==true)
    uniquevalues_out  = Dict(colname => Set{colschema.datatype}() for (colname, colschema) in colname2colschema if colschema.isunique==true)

    nr            = 0  # Total number of rows in the output data
    n_outdata     = size(outdata, 1)
    delim_outdata = output_data_file[(end - 2):end] == "csv" ? "," : "\t"
    delim_iniss   = input_issues_file[(end - 2):end] == "csv" ? "," : "\t"
    delim_outiss  = output_issues_file[(end - 2):end] == "csv" ? "," : "\t"
    quotechar     = nothing  # In some files values are delimited and quoted. E.g., line = "\"v1\", \"v2\", ...".
    colissues_in  = issues_in[:columnissues]
    colissues_out = issues_out[:columnissues]
    CSV.write(output_data_file, init_outdata(tableschema, 0); delim=delim_outdata)  # Write column headers to disk
    csvrows = CSV.Rows(input_data_file; reusebuffer=true)
    for inputrow in csvrows
        # Parse inputrow into outputrow according to ColumnSchema
        i_outdata += 1
        outputrow  = outdata[i_outdata, :]
        parserow!(outputrow, inputrow, colname2colschema)

        # Assess input row
        if pk_ncols == 1
            pk_n_missing   = pk_issues_in[:n_missing]    # Number of earlier rows (excluding the curent row) with missing primary key
            pk_n_notunique = pk_issues_in[:n_notunique]  # Number of earlier rows (excluding the curent row) with duplicated primary key
        end
        assess_row!(colissues_in, outputrow, colname2colschema, uniquevalues_in)  # Note: outputrow used because inputrow contains only Strings
        nconstraints > 0 && test_intrarow_constraints!(issues_in[:intrarow_constraints], tableschema, outputrow)
        if pk_ncols == 1
            pk_incomplete, pk_duplicated = assess_singlecolumn_primarykey!(issues_in, pk_issues_in, pk_n_missing, pk_n_notunique)
        else
            pk_incomplete, pk_duplicated = assess_multicolumn_primarykey!(issues_in, primarykey, pk_colnames, pkvalues_in, sorted_by_primarykey, outputrow)
        end

        # Assess output row
        # For speed, avoid testing value_is_valid directly. Instead reuse assessment of input.
        # Testing intra-row constraints is unnecessary because either outputrow hasn't changed or the tests return early due to missingness
        pk_has_changed = false
        for (colname, colschema) in colname2colschema
            val = outputrow[colname]
            ci  = colissues_out[colname]
            if colissues_in[colname][:n_invalid] == ci[:n_invalid]  # input value (=output value) is valid...no change to outputrow
                if ismissing(val)
                    if colschema.isrequired
                        ci[:n_missing] += 1
                    end
                else
                    if colschema.isunique
                        uniquevals = uniquevalues_out[colname]
                        if in(val, uniquevals)
                            ci[:n_notunique] += 1
                        else
                            push!(uniquevals, val)
                        end
                    end
                end
            else  # input value (=output value) is invalid...set to missing and report as missing, not as invalid
                @inbounds outputrow[colname] = missing
                pk_has_changed = pk_ncols > 1 && !pk_has_changed && in(colname, pk_colnames_set)
                ci[:n_invalid] += 1  # Required for comparison between input and output. Reset to 0 below.
                if colschema.isrequired
                    ci[:n_missing] += 1
                end
            end
        end
        if pk_has_changed || pk_incomplete
            issues_out[:primarykey_incomplete] += 1
        elseif pk_duplicated
            issues_out[:primarykey_duplicates] += 1  # Output pk is valid and equals input pk, which is duplicated (and therefore complete)
        end

        # If outdata is full append it to output_data_file
        i_outdata != n_outdata && continue
        CSV.write(output_data_file, outdata; append=true, delim=delim_outdata)
        nr        += n_outdata
        i_outdata  = 0  # Reset the row number
    end
    if i_outdata != 0
        CSV.write(output_data_file, view(outdata, 1:i_outdata, :); append=true, delim=delim_outdata)
        nr += i_outdata
    end

    # Column-level checks
    datacols_match_schemacols!(issues_in, tableschema, Set(csvrows.names))  # By construction this issue doesn't exist for outdata

    # Format result
    for (colname, colissues) in issues_out[:columnissues]
        colissues[:n_invalid] = 0  # Invalid values have been set to missing in the output data (which are then reported as missing rather than invalid)
    end
    issues_in  = construct_issues_table(issues_in,  tableschema, nr)
    issues_out = construct_issues_table(issues_out, tableschema, nr)
    CSV.write(input_issues_file,  issues_in;  delim=delim_iniss)
    CSV.write(output_issues_file, issues_out; delim=delim_outiss)
end

end


text:
# Julia wrapper for header: /opt/local/include/netcdf.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0
macro c(ret_type, func, arg_types, lib)
    local args_in = Any[ symbol(string('a',x)) for x in 1:length(arg_types.args) ]
    quote
        function $(esc(func))($(args_in...))
            retval=ccall( Libdl.dlsym(libnc,$(Expr(:quote,func))), $ret_type, $arg_types, $(args_in...) )
            if retval!=NC_NOERR
                haskey(error_description,retval) ? error("NetCDF library error: $(error_description[retval])") : error("NetCDF library error: $(retval)")
            end
        end
    end
end

using BinDeps
@BinDeps.load_dependencies [:libnetcdf]
const libnc = Libdl.dlopen(libnetcdf[1][2])

const NC_NAT = 0
const NC_BYTE = 1
const NC_CHAR = 2
const NC_SHORT = 3
const NC_INT = 4
const NC_LONG =
const NC_FLOAT = 5
const NC_DOUBLE = 6
const NC_UBYTE = 7
const NC_USHORT = 8
const NC_UINT = 9
const NC_INT64 = 10
const NC_UINT64 = 11
const NC_STRING = 12
const NC_MAX_ATOMIC_TYPE =
const NC_VLEN = 13
const NC_OPAQUE = 14
const NC_ENUM = 15
const NC_COMPOUND = 16
const NC_FIRSTUSERTYPEID = 32

const NC_FILL = 0
const NC_NOFILL = 0x100
const NC_NOWRITE = 0x0000
const NC_WRITE = 0x0001
const NC_CLOBBER = 0x0000
const NC_NOCLOBBER = 0x0004
const NC_DISKLESS = 0x0008
const NC_MMAP = 0x0010
const NC_CLASSIC_MODEL = 0x0100
const NC_64BIT_OFFSET = 0x0200
const NC_LOCK = 0x0400
const NC_SHARE = 0x0800
const NC_NETCDF4 = 0x1000
const NC_MPIIO = 0x2000
const NC_MPIPOSIX = 0x4000
const NC_PNETCDF = 0x8000
const NC_FORMAT_CLASSIC = (1)
const NC_FORMAT_64BIT = (2)
const NC_FORMAT_NETCDF4 = (3)
const NC_FORMAT_NETCDF4_CLASSIC = (4)
const NC_SIZEHINT_DEFAULT = 0
# Skipping MacroDefinition: NC_ALIGN_CHUNK((size_t)(-1))
const NC_UNLIMITED = 0
const NC_GLOBAL = -1
const NC_MAX_DIMS = 1024
const NC_MAX_ATTRS = 8192
const NC_MAX_VARS = 8192
const NC_MAX_NAME = 256
const NC_MAX_VAR_DIMS = 1024
const NC_MAX_HDF4_NAME = 64
const NC_ENDIAN_NATIVE = 0
const NC_ENDIAN_LITTLE = 1
const NC_ENDIAN_BIG = 2
const NC_CHUNKED = 0
const NC_CONTIGUOUS = 1
const NC_NOCHECKSUM = 0
const NC_FLETCHER32 = 1
const NC_NOSHUFFLE = 0
const NC_SHUFFLE = 1
# Skipping MacroDefinition: NC_ISSYSERR(err)((err)>0)
const NC_NOERR = 0
const NC2_ERR = (-1)
const NC_EBADID = (-33)
const NC_ENFILE = (-34)
const NC_EEXIST = (-35)
const NC_EINVAL = (-36)
const NC_EPERM = (-37)
const NC_ENOTINDEFINE = (-38)
const NC_EINDEFINE = (-39)
const NC_EINVALCOORDS = (-40)
const NC_EMAXDIMS = (-41)
const NC_ENAMEINUSE = (-42)
const NC_ENOTATT = (-43)
const NC_EMAXATTS = (-44)
const NC_EBADTYPE = (-45)
const NC_EBADDIM = (-46)
const NC_EUNLIMPOS = (-47)
const NC_EMAXVARS = (-48)
const NC_ENOTVAR = (-49)
const NC_EGLOBAL = (-50)
const NC_ENOTNC = (-51)
const NC_ESTS = (-52)
const NC_EMAXNAME = (-53)
const NC_EUNLIMIT = (-54)
const NC_ENORECVARS = (-55)
const NC_ECHAR = (-56)
const NC_EEDGE = (-57)
const NC_ESTRIDE = (-58)
const NC_EBADNAME = (-59)
const NC_ERANGE = (-60)
const NC_ENOMEM = (-61)
const NC_EVARSIZE = (-62)
const NC_EDIMSIZE = (-63)
const NC_ETRUNC = (-64)
const NC_EAXISTYPE = (-65)
const NC_EDAP = (-66)
const NC_ECURL = (-67)
const NC_EIO = (-68)
const NC_ENODATA = (-69)
const NC_EDAPSVC = (-70)
const NC_EDAS = (-71)
const NC_EDDS = (-72)
const NC_EDATADDS = (-73)
const NC_EDAPURL = (-74)
const NC_EDAPCONSTRAINT = (-75)
const NC_ETRANSLATION = (-76)
const NC4_FIRST_ERROR = (-100)
const NC_EHDFERR = (-101)
const NC_ECANTREAD = (-102)
const NC_ECANTWRITE = (-103)
const NC_ECANTCREATE = (-104)
const NC_EFILEMETA = (-105)
const NC_EDIMMETA = (-106)
const NC_EATTMETA = (-107)
const NC_EVARMETA = (-108)
const NC_ENOCOMPOUND = (-109)
const NC_EATTEXISTS = (-110)
const NC_ENOTNC4 = (-111)
const NC_ESTRICTNC3 = (-112)
const NC_ENOTNC3 = (-113)
const NC_ENOPAR = (-114)
const NC_EPARINIT = (-115)
const NC_EBADGRPID = (-116)
const NC_EBADTYPID = (-117)
const NC_ETYPDEFINED = (-118)
const NC_EBADFIELD = (-119)
const NC_EBADCLASS = (-120)
const NC_EMAPTYPE = (-121)
const NC_ELATEFILL = (-122)
const NC_ELATEDEF = (-123)
const NC_EDIMSCALE = (-124)
const NC_ENOGRP = (-125)
const NC_ESTORAGE = (-126)
const NC_EBADCHUNK = (-127)
const NC_ENOTBUILT = (-128)
const NC_EDISKLESS = (-129)
const NC4_LAST_ERROR = (-129)
const DIM_WITHOUT_VARIABLE = "This is a netCDF dimension but not a netCDF variable."
const NC_HAVE_NEW_CHUNKING_API = 1
# Skipping MacroDefinition: NC_EURL(NC_EDAPURL)
# Skipping MacroDefinition: NC_ECONSTRAINT(NC_EDAPCONSTRAINT)
const NC_INDEPENDENT = 0
const NC_COLLECTIVE = 1
# Skipping MacroDefinition: NC_COMPOUND_OFFSET(S,M)(offsetof(S,M))
# Skipping MacroDefinition: ncvarcpy(ncid_in,varid,ncid_out)ncvarcopy((ncid_in),(varid),(ncid_out))
# Skipping MacroDefinition: nc_set_log_level(e)

typealias nc_type Cint
typealias nclong Cint

@c Ptr{UInt8} nc_inq_libvers () libnetcdf
@c Ptr{UInt8} nc_strerror (Cint,) libnetcdf
@c Cint nc__create (Ptr{UInt8}, Cint, Csize_t, Ptr{Csize_t}, Ptr{Cint}) libnetcdf
@c Cint nc_create (Ptr{UInt8}, Cint, Ptr{Cint}) libnetcdf
@c Cint nc__open (Ptr{UInt8}, Cint, Ptr{Csize_t}, Ptr{Cint}) libnetcdf
@c Cint nc_open (Ptr{UInt8}, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_inq_path (Cint, Ptr{Csize_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_var_par_access (Cint, Cint, Cint) libnetcdf
@c Cint nc_inq_ncid (Cint, Ptr{UInt8}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_grps (Cint, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_grpname (Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_inq_grpname_full (Cint, Ptr{Csize_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_inq_grpname_len (Cint, Ptr{Csize_t}) libnetcdf
@c Cint nc_inq_grp_parent (Cint, Ptr{Cint}) libnetcdf
@c Cint nc_inq_grp_ncid (Cint, Ptr{UInt8}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_grp_full_ncid (Cint, Ptr{UInt8}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_varids (Cint, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_dimids (Cint, Ptr{Cint}, Ptr{Cint}, Cint) libnetcdf
@c Cint nc_inq_typeids (Cint, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_type_equal (Cint, nc_type, Cint, nc_type, Ptr{Cint}) libnetcdf
@c Cint nc_def_grp (Cint, Ptr{UInt8}, Ptr{Cint}) libnetcdf
@c Cint nc_def_compound (Cint, Csize_t, Ptr{UInt8}, Ptr{nc_type}) libnetcdf
@c Cint nc_insert_compound (Cint, nc_type, Ptr{UInt8}, Csize_t, nc_type) libnetcdf
@c Cint nc_insert_array_compound (Cint, nc_type, Ptr{UInt8}, Csize_t, nc_type, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_inq_type (Cint, nc_type, Ptr{UInt8}, Ptr{Csize_t}) libnetcdf
@c Cint nc_inq_typeid (Cint, Ptr{UInt8}, Ptr{nc_type}) libnetcdf
@c Cint nc_inq_compound (Cint, nc_type, Ptr{UInt8}, Ptr{Csize_t}, Ptr{Csize_t}) libnetcdf
@c Cint nc_inq_compound_name (Cint, nc_type, Ptr{UInt8}) libnetcdf
@c Cint nc_inq_compound_size (Cint, nc_type, Ptr{Csize_t}) libnetcdf
@c Cint nc_inq_compound_nfields (Cint, nc_type, Ptr{Csize_t}) libnetcdf
@c Cint nc_inq_compound_field (Cint, nc_type, Cint, Ptr{UInt8}, Ptr{Csize_t}, Ptr{nc_type}, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_compound_fieldname (Cint, nc_type, Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_inq_compound_fieldindex (Cint, nc_type, Ptr{UInt8}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_compound_fieldoffset (Cint, nc_type, Cint, Ptr{Csize_t}) libnetcdf
@c Cint nc_inq_compound_fieldtype (Cint, nc_type, Cint, Ptr{nc_type}) libnetcdf
@c Cint nc_inq_compound_fieldndims (Cint, nc_type, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_inq_compound_fielddim_sizes (Cint, nc_type, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_def_vlen (Cint, Ptr{UInt8}, nc_type, Ptr{nc_type}) libnetcdf
@c Cint nc_inq_vlen (Cint, nc_type, Ptr{UInt8}, Ptr{Csize_t}, Ptr{nc_type}) libnetcdf
@c Cint nc_free_vlen (Ptr{nc_vlen_t},) libnetcdf
@c Cint nc_free_vlens (Csize_t, Ptr{nc_vlen_t}) libnetcdf
@c Cint nc_put_vlen_element (Cint, Cint, Ptr{Void}, Csize_t, Ptr{Void}) libnetcdf
@c Cint nc_get_vlen_element (Cint, Cint, Ptr{Void}, Ptr{Csize_t}, Ptr{Void}) libnetcdf
@c Cint nc_free_string (Csize_t, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_inq_user_type (Cint, nc_type, Ptr{UInt8}, Ptr{Csize_t}, Ptr{nc_type}, Ptr{Csize_t}, Ptr{Cint}) libnetcdf
@c Cint nc_put_att (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{Void}) libnetcdf
@c Cint nc_get_att (Cint, Cint, Ptr{UInt8}, Ptr{Void}) libnetcdf
@c Cint nc_def_enum (Cint, nc_type, Ptr{UInt8}, Ptr{nc_type}) libnetcdf
@c Cint nc_insert_enum (Cint, nc_type, Ptr{UInt8}, Ptr{Void}) libnetcdf
@c Cint nc_inq_enum (Cint, nc_type, Ptr{UInt8}, Ptr{nc_type}, Ptr{Csize_t}, Ptr{Csize_t}) libnetcdf
@c Cint nc_inq_enum_member (Cint, nc_type, Cint, Ptr{UInt8}, Ptr{Void}) libnetcdf
@c Cint nc_inq_enum_ident (Cint, nc_type, Clonglong, Ptr{UInt8}) libnetcdf
@c Cint nc_def_opaque (Cint, Csize_t, Ptr{UInt8}, Ptr{nc_type}) libnetcdf
@c Cint nc_inq_opaque (Cint, nc_type, Ptr{UInt8}, Ptr{Csize_t}) libnetcdf
@c Cint nc_put_var (Cint, Cint, Ptr{Void}) libnetcdf
@c Cint nc_get_var (Cint, Cint, Ptr{Void}) libnetcdf
@c Cint nc_put_var1 (Cint, Cint, Ptr{Csize_t}, Ptr{Void}) libnetcdf
@c Cint nc_get_var1 (Cint, Cint, Ptr{Csize_t}, Ptr{Void}) libnetcdf
@c Cint nc_put_vara (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Void}) libnetcdf
@c Cint nc_get_vara (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Void}) libnetcdf
@c Cint nc_put_vars (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Void}) libnetcdf
@c Cint nc_get_vars (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Void}) libnetcdf
@c Cint nc_put_varm (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Void}) libnetcdf
@c Cint nc_get_varm (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Void}) libnetcdf
@c Cint nc_def_var_deflate (Cint, Cint, Cint, Cint, Cint) libnetcdf
@c Cint nc_inq_var_deflate (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_var_szip (Cint, Cint, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint nc_def_var_fletcher32 (Cint, Cint, Cint) libnetcdf
@c Cint nc_inq_var_fletcher32 (Cint, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_def_var_chunking (Cint, Cint, Cint, Ptr{Csize_t}) libnetcdf
@c Cint nc_inq_var_chunking (Cint, Cint, Ptr{Cint}, Ptr{Csize_t}) libnetcdf
@c Cint nc_def_var_fill (Cint, Cint, Cint, Ptr{Void}) libnetcdf
@c Cint nc_inq_var_fill (Cint, Cint, Ptr{Cint}, Ptr{Void}) libnetcdf
@c Cint nc_def_var_endian (Cint, Cint, Cint) libnetcdf
@c Cint nc_inq_var_endian (Cint, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_set_fill (Cint, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_set_default_format (Cint, Ptr{Cint}) libnetcdf
@c Cint nc_set_chunk_cache (Csize_t, Csize_t, Cfloat) libnetcdf
@c Cint nc_get_chunk_cache (Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cfloat}) libnetcdf
@c Cint nc_set_var_chunk_cache (Cint, Cint, Csize_t, Csize_t, Cfloat) libnetcdf
@c Cint nc_get_var_chunk_cache (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cfloat}) libnetcdf
@c Cint nc_redef (Cint,) libnetcdf
@c Cint nc__enddef (Cint, Csize_t, Csize_t, Csize_t, Csize_t) libnetcdf
@c Cint nc_enddef (Cint,) libnetcdf
@c Cint nc_sync (Cint,) libnetcdf
@c Cint nc_abort (Cint,) libnetcdf
@c Cint nc_close (Cint,) libnetcdf
@c Cint nc_inq (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_ndims (Cint, Ptr{Cint}) libnetcdf
@c Cint nc_inq_nvars (Cint, Ptr{Cint}) libnetcdf
@c Cint nc_inq_natts (Cint, Ptr{Cint}) libnetcdf
@c Cint nc_inq_unlimdim (Cint, Ptr{Cint}) libnetcdf
@c Cint nc_inq_unlimdims (Cint, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_format (Cint, Ptr{Cint}) libnetcdf
@c Cint nc_def_dim (Cint, Ptr{UInt8}, Csize_t, Ptr{Cint}) libnetcdf
@c Cint nc_inq_dimid (Cint, Ptr{UInt8}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_dim (Cint, Cint, Ptr{UInt8}, Ptr{Csize_t}) libnetcdf
@c Cint nc_inq_dimname (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_inq_dimlen (Cint, Cint, Ptr{Csize_t}) libnetcdf
@c Cint nc_rename_dim (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_inq_att (Cint, Cint, Ptr{UInt8}, Ptr{nc_type}, Ptr{Csize_t}) libnetcdf
@c Cint nc_inq_attid (Cint, Cint, Ptr{UInt8}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_atttype (Cint, Cint, Ptr{UInt8}, Ptr{nc_type}) libnetcdf
@c Cint nc_inq_attlen (Cint, Cint, Ptr{UInt8}, Ptr{Csize_t}) libnetcdf
@c Cint nc_inq_attname (Cint, Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_copy_att (Cint, Cint, Ptr{UInt8}, Cint, Cint) libnetcdf
@c Cint nc_rename_att (Cint, Cint, Ptr{UInt8}, Ptr{UInt8}) libnetcdf
@c Cint nc_del_att (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_put_att_text (Cint, Cint, Ptr{UInt8}, Csize_t, Ptr{UInt8}) libnetcdf
@c Cint nc_get_att_text (Cint, Cint, Ptr{UInt8}, Ptr{UInt8}) libnetcdf
@c Cint nc_put_att_uchar (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_att_uchar (Cint, Cint, Ptr{UInt8}, Ptr{Cuchar}) libnetcdf
@c Cint nc_put_att_schar (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{Int8}) libnetcdf
@c Cint nc_get_att_schar (Cint, Cint, Ptr{UInt8}, Ptr{Int8}) libnetcdf
@c Cint nc_put_att_short (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{Int16}) libnetcdf
@c Cint nc_get_att_short (Cint, Cint, Ptr{UInt8}, Ptr{Int16}) libnetcdf
@c Cint nc_put_att_int (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{Cint}) libnetcdf
@c Cint nc_get_att_int (Cint, Cint, Ptr{UInt8}, Ptr{Cint}) libnetcdf
@c Cint nc_put_att_long (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{Clong}) libnetcdf
@c Cint nc_get_att_long (Cint, Cint, Ptr{UInt8}, Ptr{Clong}) libnetcdf
@c Cint nc_put_att_float (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{Cfloat}) libnetcdf
@c Cint nc_get_att_float (Cint, Cint, Ptr{UInt8}, Ptr{Cfloat}) libnetcdf
@c Cint nc_put_att_double (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{Cdouble}) libnetcdf
@c Cint nc_get_att_double (Cint, Cint, Ptr{UInt8}, Ptr{Cdouble}) libnetcdf
@c Cint nc_put_att_ushort (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{UInt16}) libnetcdf
@c Cint nc_get_att_ushort (Cint, Cint, Ptr{UInt8}, Ptr{UInt16}) libnetcdf
@c Cint nc_put_att_uint (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{UInt32}) libnetcdf
@c Cint nc_get_att_uint (Cint, Cint, Ptr{UInt8}, Ptr{UInt32}) libnetcdf
@c Cint nc_put_att_longlong (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{Clonglong}) libnetcdf
@c Cint nc_get_att_longlong (Cint, Cint, Ptr{UInt8}, Ptr{Clonglong}) libnetcdf
@c Cint nc_put_att_ulonglong (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{Culonglong}) libnetcdf
@c Cint nc_get_att_ulonglong (Cint, Cint, Ptr{UInt8}, Ptr{Culonglong}) libnetcdf
@c Cint nc_put_att_string (Cint, Cint, Ptr{UInt8}, Csize_t, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_get_att_string (Cint, Cint, Ptr{UInt8}, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_def_var (Cint, Ptr{UInt8}, nc_type, Cint, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_var (Cint, Cint, Ptr{UInt8}, Ptr{nc_type}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_varid (Cint, Ptr{UInt8}, Ptr{Cint}) libnetcdf
@c Cint nc_inq_varname (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_inq_vartype (Cint, Cint, Ptr{nc_type}) libnetcdf
@c Cint nc_inq_varndims (Cint, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_inq_vardimid (Cint, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_inq_varnatts (Cint, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_rename_var (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_copy_var (Cint, Cint, Cint) libnetcdf
@c Cint nc_put_var1_text (Cint, Cint, Ptr{Csize_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_get_var1_text (Cint, Cint, Ptr{Csize_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_put_var1_uchar (Cint, Cint, Ptr{Csize_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_var1_uchar (Cint, Cint, Ptr{Csize_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_put_var1_schar (Cint, Cint, Ptr{Csize_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_get_var1_schar (Cint, Cint, Ptr{Csize_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_put_var1_short (Cint, Cint, Ptr{Csize_t}, Ptr{Int16}) libnetcdf
@c Cint nc_get_var1_short (Cint, Cint, Ptr{Csize_t}, Ptr{Int16}) libnetcdf
@c Cint nc_put_var1_int (Cint, Cint, Ptr{Csize_t}, Ptr{Cint}) libnetcdf
@c Cint nc_get_var1_int (Cint, Cint, Ptr{Csize_t}, Ptr{Cint}) libnetcdf
@c Cint nc_put_var1_long (Cint, Cint, Ptr{Csize_t}, Ptr{Clong}) libnetcdf
@c Cint nc_get_var1_long (Cint, Cint, Ptr{Csize_t}, Ptr{Clong}) libnetcdf
@c Cint nc_put_var1_float (Cint, Cint, Ptr{Csize_t}, Ptr{Cfloat}) libnetcdf
@c Cint nc_get_var1_float (Cint, Cint, Ptr{Csize_t}, Ptr{Cfloat}) libnetcdf
@c Cint nc_put_var1_double (Cint, Cint, Ptr{Csize_t}, Ptr{Cdouble}) libnetcdf
@c Cint nc_get_var1_double (Cint, Cint, Ptr{Csize_t}, Ptr{Cdouble}) libnetcdf
@c Cint nc_put_var1_ushort (Cint, Cint, Ptr{Csize_t}, Ptr{UInt16}) libnetcdf
@c Cint nc_get_var1_ushort (Cint, Cint, Ptr{Csize_t}, Ptr{UInt16}) libnetcdf
@c Cint nc_put_var1_uint (Cint, Cint, Ptr{Csize_t}, Ptr{UInt32}) libnetcdf
@c Cint nc_get_var1_uint (Cint, Cint, Ptr{Csize_t}, Ptr{UInt32}) libnetcdf
@c Cint nc_put_var1_longlong (Cint, Cint, Ptr{Csize_t}, Ptr{Clonglong}) libnetcdf
@c Cint nc_get_var1_longlong (Cint, Cint, Ptr{Csize_t}, Ptr{Clonglong}) libnetcdf
@c Cint nc_put_var1_ulonglong (Cint, Cint, Ptr{Csize_t}, Ptr{Culonglong}) libnetcdf
@c Cint nc_get_var1_ulonglong (Cint, Cint, Ptr{Csize_t}, Ptr{Culonglong}) libnetcdf
@c Cint nc_put_var1_string (Cint, Cint, Ptr{Csize_t}, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_get_var1_string (Cint, Cint, Ptr{Csize_t}, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_put_vara_text (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_get_vara_text (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_put_vara_uchar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_vara_uchar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_put_vara_schar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Int8}) libnetcdf
@c Cint nc_get_vara_schar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Int8}) libnetcdf
@c Cint nc_put_vara_short (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Int16}) libnetcdf
@c Cint nc_get_vara_short (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Int16}) libnetcdf
@c Cint nc_put_vara_int (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cint}) libnetcdf
@c Cint nc_get_vara_int (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cint}) libnetcdf
@c Cint nc_put_vara_long (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Clong}) libnetcdf
@c Cint nc_get_vara_long (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Clong}) libnetcdf
@c Cint nc_put_vara_float (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cfloat}) libnetcdf
@c Cint nc_get_vara_float (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cfloat}) libnetcdf
@c Cint nc_put_vara_double (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cdouble}) libnetcdf
@c Cint nc_get_vara_double (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cdouble}) libnetcdf
@c Cint nc_put_vara_ushort (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{UInt16}) libnetcdf
@c Cint nc_get_vara_ushort (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{UInt16}) libnetcdf
@c Cint nc_put_vara_uint (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{UInt32}) libnetcdf
@c Cint nc_get_vara_uint (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{UInt32}) libnetcdf
@c Cint nc_put_vara_longlong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Clonglong}) libnetcdf
@c Cint nc_get_vara_longlong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Clonglong}) libnetcdf
@c Cint nc_put_vara_ulonglong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Culonglong}) libnetcdf
@c Cint nc_get_vara_ulonglong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Culonglong}) libnetcdf
@c Cint nc_put_vara_string (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_get_vara_string (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_put_vars_text (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_get_vars_text (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_put_vars_uchar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_vars_uchar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_put_vars_schar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_get_vars_schar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_put_vars_short (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Int16}) libnetcdf
@c Cint nc_get_vars_short (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Int16}) libnetcdf
@c Cint nc_put_vars_int (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Cint}) libnetcdf
@c Cint nc_get_vars_int (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Cint}) libnetcdf
@c Cint nc_put_vars_long (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Clong}) libnetcdf
@c Cint nc_get_vars_long (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Clong}) libnetcdf
@c Cint nc_put_vars_float (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Cfloat}) libnetcdf
@c Cint nc_get_vars_float (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Cfloat}) libnetcdf
@c Cint nc_put_vars_double (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Cdouble}) libnetcdf
@c Cint nc_get_vars_double (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Cdouble}) libnetcdf
@c Cint nc_put_vars_ushort (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{UInt16}) libnetcdf
@c Cint nc_get_vars_ushort (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{UInt16}) libnetcdf
@c Cint nc_put_vars_uint (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{UInt32}) libnetcdf
@c Cint nc_get_vars_uint (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{UInt32}) libnetcdf
@c Cint nc_put_vars_longlong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Clonglong}) libnetcdf
@c Cint nc_get_vars_longlong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Clonglong}) libnetcdf
@c Cint nc_put_vars_ulonglong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Culonglong}) libnetcdf
@c Cint nc_get_vars_ulonglong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Culonglong}) libnetcdf
@c Cint nc_put_vars_string (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_get_vars_string (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_put_varm_text (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_get_varm_text (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_put_varm_uchar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_varm_uchar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_put_varm_schar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_get_varm_schar (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{UInt8}) libnetcdf
@c Cint nc_put_varm_short (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Int16}) libnetcdf
@c Cint nc_get_varm_short (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Int16}) libnetcdf
@c Cint nc_put_varm_int (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Cint}) libnetcdf
@c Cint nc_get_varm_int (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Cint}) libnetcdf
@c Cint nc_put_varm_long (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Clong}) libnetcdf
@c Cint nc_get_varm_long (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Clong}) libnetcdf
@c Cint nc_put_varm_float (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Cfloat}) libnetcdf
@c Cint nc_get_varm_float (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Cfloat}) libnetcdf
@c Cint nc_put_varm_double (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Cdouble}) libnetcdf
@c Cint nc_get_varm_double (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Cdouble}) libnetcdf
@c Cint nc_put_varm_ushort (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{UInt16}) libnetcdf
@c Cint nc_get_varm_ushort (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{UInt16}) libnetcdf
@c Cint nc_put_varm_uint (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{UInt32}) libnetcdf
@c Cint nc_get_varm_uint (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{UInt32}) libnetcdf
@c Cint nc_put_varm_longlong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Clonglong}) libnetcdf
@c Cint nc_get_varm_longlong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Clonglong}) libnetcdf
@c Cint nc_put_varm_ulonglong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Culonglong}) libnetcdf
@c Cint nc_get_varm_ulonglong (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Culonglong}) libnetcdf
@c Cint nc_put_varm_string (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_get_varm_string (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_put_var_text (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_get_var_text (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_put_var_uchar (Cint, Cint, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_var_uchar (Cint, Cint, Ptr{Cuchar}) libnetcdf
@c Cint nc_put_var_schar (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_get_var_schar (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint nc_put_var_short (Cint, Cint, Ptr{Int16}) libnetcdf
@c Cint nc_get_var_short (Cint, Cint, Ptr{Int16}) libnetcdf
@c Cint nc_put_var_int (Cint, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_get_var_int (Cint, Cint, Ptr{Cint}) libnetcdf
@c Cint nc_put_var_long (Cint, Cint, Ptr{Clong}) libnetcdf
@c Cint nc_get_var_long (Cint, Cint, Ptr{Clong}) libnetcdf
@c Cint nc_put_var_float (Cint, Cint, Ptr{Cfloat}) libnetcdf
@c Cint nc_get_var_float (Cint, Cint, Ptr{Cfloat}) libnetcdf
@c Cint nc_put_var_double (Cint, Cint, Ptr{Cdouble}) libnetcdf
@c Cint nc_get_var_double (Cint, Cint, Ptr{Cdouble}) libnetcdf
@c Cint nc_put_var_ushort (Cint, Cint, Ptr{UInt16}) libnetcdf
@c Cint nc_get_var_ushort (Cint, Cint, Ptr{UInt16}) libnetcdf
@c Cint nc_put_var_uint (Cint, Cint, Ptr{UInt32}) libnetcdf
@c Cint nc_get_var_uint (Cint, Cint, Ptr{UInt32}) libnetcdf
@c Cint nc_put_var_longlong (Cint, Cint, Ptr{Clonglong}) libnetcdf
@c Cint nc_get_var_longlong (Cint, Cint, Ptr{Clonglong}) libnetcdf
@c Cint nc_put_var_ulonglong (Cint, Cint, Ptr{Culonglong}) libnetcdf
@c Cint nc_get_var_ulonglong (Cint, Cint, Ptr{Culonglong}) libnetcdf
@c Cint nc_put_var_string (Cint, Cint, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_get_var_string (Cint, Cint, Ptr{Ptr{UInt8}}) libnetcdf
@c Cint nc_put_att_ubyte (Cint, Cint, Ptr{UInt8}, nc_type, Csize_t, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_att_ubyte (Cint, Cint, Ptr{UInt8}, Ptr{Cuchar}) libnetcdf
@c Cint nc_put_var1_ubyte (Cint, Cint, Ptr{Csize_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_var1_ubyte (Cint, Cint, Ptr{Csize_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_put_vara_ubyte (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_vara_ubyte (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_put_vars_ubyte (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_vars_ubyte (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_put_varm_ubyte (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_varm_ubyte (Cint, Cint, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{ptrdiff_t}, Ptr{ptrdiff_t}, Ptr{Cuchar}) libnetcdf
@c Cint nc_put_var_ubyte (Cint, Cint, Ptr{Cuchar}) libnetcdf
@c Cint nc_get_var_ubyte (Cint, Cint, Ptr{Cuchar}) libnetcdf
@c Cint nc_show_metadata (Cint,) libnetcdf
@c Cint nc__create_mp (Ptr{UInt8}, Cint, Csize_t, Cint, Ptr{Csize_t}, Ptr{Cint}) libnetcdf
@c Cint nc__open_mp (Ptr{UInt8}, Cint, Cint, Ptr{Csize_t}, Ptr{Cint}) libnetcdf
@c Cint nc_delete (Ptr{UInt8},) libnetcdf
@c Cint nc_delete_mp (Ptr{UInt8}, Cint) libnetcdf
@c Cint nc_set_base_pe (Cint, Cint) libnetcdf
@c Cint nc_inq_base_pe (Cint, Ptr{Cint}) libnetcdf
@c Cint nctypelen (nc_type,) libnetcdf
@c None nc_advise (Ptr{UInt8}, Cint, Ptr{UInt8}) libnetcdf
@c Cint nccreate (Ptr{UInt8}, Cint) libnetcdf
@c Cint ncopen (Ptr{UInt8}, Cint) libnetcdf
@c Cint ncsetfill (Cint, Cint) libnetcdf
@c Cint ncredef (Cint,) libnetcdf
@c Cint ncendef (Cint,) libnetcdf
@c Cint ncsync (Cint,) libnetcdf
@c Cint ncabort (Cint,) libnetcdf
@c Cint ncclose (Cint,) libnetcdf
@c Cint ncinquire (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint ncdimdef (Cint, Ptr{UInt8}, Clong) libnetcdf
@c Cint ncdimid (Cint, Ptr{UInt8}) libnetcdf
@c Cint ncdiminq (Cint, Cint, Ptr{UInt8}, Ptr{Clong}) libnetcdf
@c Cint ncdimrename (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint ncattput (Cint, Cint, Ptr{UInt8}, nc_type, Cint, Ptr{Void}) libnetcdf
@c Cint ncattinq (Cint, Cint, Ptr{UInt8}, Ptr{nc_type}, Ptr{Cint}) libnetcdf
@c Cint ncattget (Cint, Cint, Ptr{UInt8}, Ptr{Void}) libnetcdf
@c Cint ncattcopy (Cint, Cint, Ptr{UInt8}, Cint, Cint) libnetcdf
@c Cint ncattname (Cint, Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint ncattrename (Cint, Cint, Ptr{UInt8}, Ptr{UInt8}) libnetcdf
@c Cint ncattdel (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint ncvardef (Cint, Ptr{UInt8}, nc_type, Cint, Ptr{Cint}) libnetcdf
@c Cint ncvarid (Cint, Ptr{UInt8}) libnetcdf
@c Cint ncvarinq (Cint, Cint, Ptr{UInt8}, Ptr{nc_type}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}) libnetcdf
@c Cint ncvarput1 (Cint, Cint, Ptr{Clong}, Ptr{Void}) libnetcdf
@c Cint ncvarget1 (Cint, Cint, Ptr{Clong}, Ptr{Void}) libnetcdf
@c Cint ncvarput (Cint, Cint, Ptr{Clong}, Ptr{Clong}, Ptr{Void}) libnetcdf
@c Cint ncvarget (Cint, Cint, Ptr{Clong}, Ptr{Clong}, Ptr{Void}) libnetcdf
@c Cint ncvarputs (Cint, Cint, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Void}) libnetcdf
@c Cint ncvargets (Cint, Cint, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Void}) libnetcdf
@c Cint ncvarputg (Cint, Cint, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Void}) libnetcdf
@c Cint ncvargetg (Cint, Cint, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Void}) libnetcdf
@c Cint ncvarrename (Cint, Cint, Ptr{UInt8}) libnetcdf
@c Cint ncrecinq (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Clong}) libnetcdf
@c Cint ncrecget (Cint, Clong, Ptr{Ptr{Void}}) libnetcdf
@c Cint ncrecput (Cint, Clong, Ptr{Ptr{Void}}) libnetcdf


text:
using Optim
using Statistics

  """
      objFunctionMSE( p::Vector{Float64}; predMat::Matrix{Float64}, t::Vector{Float64} )

      evaluate MeanSquaredError under given params.
      #Arguments
      - `p::Vector{Float64}`:       initial solution.
      - `predMat::Matrix{Float64}`: each column represents predictions of one model.
      - `t::Vector{Float64}`:       ground truth labels.
      #Return
      - `Float64`:                  MeanSquaredError.
  """
  function objFunctionMSE( p::Vector{Float64}; predMat::Matrix{Float64}, t::Vector{Float64} )
    return mean( lossFunctions.MSE( predMat * p, t ) );
  end

  """
      objFunctionHinge( p::Vector{Float64}; predMat::Matrix{Float64}, t::Vector{Float64} )

      evaluate MeanSquaredError under given params.
      #Arguments
      - `p::Vector{Float64}`:       initial solution.
      - `predMat::Matrix{Float64}`: each column represents predictions of one model.
      - `t::Vector{Float64}`:       ground truth labels.
      #Return
      - `Float64`:                  hingeLoss.
  """
  function objFunctionHinge( p::Vector{Float64}; predMat::Matrix{Float64}, t::Vector{Float64} )
    return mean( lossFunctions.hingeLoss( predMat * p, t ) );
  end

  """
      directOptimNaiveMSE( predMat::Matrix{Float64}, t::Vector{Float64}, p::Vector{Float64} )

      compute refined solution _for_ mixing parameter p.
      #Arguments
      - `predMat::Matrix{Float64}`: each column is the prediction _of_ one hypothesis.
      - `t::Vector{Float64}`:       label vector.
      - `p::Vector{Float64}`:       initial solution for mixing coefficients.
      #Return
      - `Vector{Float64}`:          improved initial solution.
  """
  function directOptimNaiveMSE( predMat::Matrix{Float64}, t::Vector{Float64}, p::Vector{Float64} )
    @assert size( t, 1 ) == size( t, 1 )
    numberSamples = size( predMat, 1 );
    width         = size( predMat, 2 );
    f(p)          = objFunctionMSE( p,predMat=predMat, t=t )
    resultObject  = optimize( f, p );
    Optim.minimizer( resultObject )
  end

  """
      directOptimHinge( predMat::Matrix{Float64}, t::Vector{Float64}, p::Vector{Float64} )

      compute refined solution _for_ mixing parameter p.
      #Arguments
      - `predMat::Matrix{Float64}`: each column is the prediction _of_ one hypothesis.
      - `t::Vector{Float64}`:       label vector.
      - `p::Vector{Float64}`:       initial solution for mixing coefficients.
      #Return
      - `Vector{Float64}`:          improved initial solution.
  """
  function directOptimHinge( predMat::Matrix{Float64}, t::Vector{Float64}, p::Vector{Float64} )
    @assert size( t, 1 ) == size( t, 1 )
    numberSamples = size( predMat, 1 );
    width         = size( predMat, 2 );
    f(p)          = objFunctionHinge( p,predMat=predMat, t=t );
    resultObject  = optimize( f, p );
    Optim.minimizer( resultObject )
  end


text:
using LinearAlgebra
# using SpecialFunctions

"""
    isdependent returns true if eq is dependent on x
"""
isdependent(eq, x) = !isequal(expand_derivatives(Differential(x)(eq)), 0)

Base.signbit(z::Complex{T}) where T<:Number = signbit(real(z))

# this is the main heurisctic used to find the test fragments
function generate_basis(eq, x, h=[])
    Δeq = expand_derivatives(Differential(x)(eq))
    kers = expand(eq + Δeq)
    # return [one(x); candidates(kers, x); h]
    return [one(x); candidates(kers, x); h]
end

"""
    candidates returns a list of candidate expressions to form the integration
    basis
"""
candidates(eq, x) = isdependent(eq,x) ? [eq] : []

candidates(eq::Num, x) = candidates(value(eq), x)

# the candidates of an Add is the union of the candidates of the terms
# ∫ Σᵢ fᵢ(x) dx = Σᵢ ∫ fᵢ(x) dx
candidates(eq::SymbolicUtils.Add, x) = unique(∪([candidates(t,x) for t in arguments(eq)]...))

# the candidates of a Mul is the outer product of the candidates of the terms
# d(uv)/dx = u dv/dx + v + du/dx
function candidates(eq::SymbolicUtils.Mul, x)
    terms = [candidates(q,x) for q in arguments(eq)]
    n = length(terms)

    l = Any[one(x)]

    for j = 1:n
        m = length(l)
        for t in terms[j]
            for k = 1:m
                push!(l, l[k]*t)
            end
        end
    end

    unique(l[2:end])    # removing the initial 1
end

# the candidates of a Pow encode different integration rules
function candidates(eq::SymbolicUtils.Pow, x)
    if !isdependent(eq,x) return [one(x)] end

    p = arguments(eq)[1]    # eq = p ^ k
    k = arguments(eq)[2]

    if k < 0 && k ≈ round(k)
        return candidate_pow_minus(p, k, x)
    elseif k ≈ 0.5 || k ≈ -0.5
        # ∫ √f(x) dx = ... + c * log(df/dx + √f) if deg(f) == 2
        Δ = expand_derivatives(Differential(x)(p))
        # return [[p^k, p^(k+1)]; log(abs(0.5*Δ + sqrt(p)))]
        return [[p^k, p^(k+1)]; log(0.5*Δ + sqrt(p))]
    end

    # ∫ p^k dp = c * p^(k+1)
    return [p^k, p^(k+1)]
end

function candidate_pow_minus(p, k, x)
    q = to_poly(p, x)   # q is a DynamicPolynomials Polynomial
                        # the reason is to be able to use find_roots
    if q == nothing return [p^k, p^(k+1), log(p)] end

    q = 0 + 1*q         # take care of monomials and constants
    r, s = find_roots(q, var(q))
    s = s[1:2:end]
    r = nice_parameter.(r)
    s = nice_parameter.(s)
    # s = Complex.(nice_parameters(real.(s)), nice_parameters(imag(s)))

    # ∫ 1 / ((x-z₁)(x-z₂)) dx = ... + c₁ * log(x-z₁) + c₂ * log(x-z₂)
    # q = sum(log(x - u) for u in r; init=zero(x)) +
    #     sum(atan((x - real(u))/imag(u)) for u in s; init=zero(x)) +
    #     sum(log(x^2 - 2*real(u)*x + abs2(u)) for u in s; init=zero(x))

    q₁ = [[log(x - u) for u in r];
          [atan((x - real(u))/imag(u)) for u in s];
          [log(x^2 - 2*real(u)*x + abs2(u)) for u in s]
         ]

    q₂ = [[(x - u)^-1 for u in r];
          [(x^2 - 2*real(u)*x + abs2(u))^-1 for u in s];
          [x*(x^2 - 2*real(u)*x + abs2(u))^-1 for u in s]
         ]

    # return [[p^k, p^(k+1)]; candidates(q₁, x)]
    if k ≈ -1
        return [[p^k]; q₁; q₂]
    else
        return [[p^k, p^(k+1)]; q₁; q₂]
    end
end

###############################################################################

"""
    integrate is the main entry point

    input:
    ------
    eq: a Symbolics expression to integrate
    abstol: the desired tolerance
    num_steps: the number of different steps with expanding basis to be tried
    num_trials: the number of trials in each step (no changes to the basis)
    lo and hi: the range used to generate random values of x (the independent variable)
    show_basis: if true, the basis is printed

    output:
    -------
    solved, unsolved

    a pair of expressions, solved is the solved integral and unsolved is the residual unsolved
    portion of the input
"""
function integrate(eq; kwargs...)
    x = var(eq)
    integrate(eq, x; kwargs...)
end

function integrate(eq, x; abstol=1e-6, num_steps=2, num_trials=10, radius=5.0,
                   show_basis=false, opt = STLSQ(exp.(-10:1:0)), bypass=false,
                   attempt_ratio=5, symbolic=false, bypart=true, max_basis=110,
                   verbose=false)
    eq = expand(eq)

    # eq is a constant
    if x == nothing
        @syms 𝑥
        return 𝑥 * eq, 0, 0
    end

    # check if eq is a rational function
    # if so, we perform a partial-fraction decomposition first (the first part of the Hermite's method)

    # q = to_rational(eq, x)
    # if q != nothing
    #     eq = q
    # end

    s₁, u₁, ϵ = integrate_sum(eq, x; bypass, abstol, num_trials, num_steps,
                              radius, show_basis, opt, attempt_ratio, symbolic, max_basis, verbose)

    if isequal(u₁, 0) || !bypart
        return s₁, u₁, ϵ
    else
        s₂, u₂, ϵ = try_integration_by_parts(u₁, x; abstol, num_trials, num_steps,
                                             radius, show_basis, opt, attempt_ratio,
                                             symbolic, max_basis, verbose)
        return s₁ + s₂, u₂, ϵ
    end
end

"""
    ∫ Σᵢ fᵢ(x) dx = Σᵢ ∫ fᵢ(x) dx
"""
function integrate_sum(eq::SymbolicUtils.Add, x; bypass=false, kwargs...)
    if bypass
        return integrate_term(eq, x; kwargs...)
    else
        solved = 0
        unsolved = 0
        ϵ₀ = 0

        for p in arguments(eq)
            s, u, ϵ = integrate_term(p, x; kwargs...)
            solved += s
            unsolved += u
            ϵ₀ = max(ϵ₀, ϵ)
        end

        return solved, unsolved, ϵ₀
    end
end

function integrate_sum(eq, x; kwargs...)
    integrate_term(eq, x; kwargs...)
end

function accept_solution(eq, x, sol; abstol=1e-6)
    try
        Δ = substitute(expand_derivatives(Differential(x)(sol)-eq), Dict(x => Complex(rand())))
        return abs(Δ) < abstol
    catch e
        #
    end
    return false
end

function integrate_term(eq, x; kwargs...)
    args = Dict(kwargs)
    abstol, num_steps, num_trials, show_basis, symbolic, verbose, max_basis, radius =
        args[:abstol], args[:num_steps], args[:num_trials], args[:show_basis],
        args[:symbolic], args[:verbose], args[:max_basis], args[:radius]

    # note that the order of the operations is important!
    # first, collecing hints, then applying transformation rules, and finally finding the basis.
    h = collect_hints(eq, x)
    eq = apply_integration_rules(eq)
    basis = generate_basis(eq, x, h)

    if verbose printstyled("|β| = ", length(basis), ". "; color=:yellow) end
    if length(basis) > max_basis return 0, eq, Inf end

    D = Differential(x)
    ϵ₀ = Inf
    y₀ = 0

    for i = 1:num_steps
        basis = unique([basis; basis*x])
        Δbasis = [expand_derivatives(D(f)) for f in basis]
        if show_basis println(basis); println(Δbasis) end

        if symbolic
            y, ϵ = try_symbolic(Float64, eq, x, basis, Δbasis; kwargs...)
            if !isequal(y, 0) && accept_solution(eq, x, y; abstol)
                if verbose printstyled("$i, symbolic\n"; color=:yellow) end
                return y, 0, 0
            end
        end

        for j = 1:num_trials
            y, ϵ = try_integrate(Complex{Float64}, eq, x, basis, Δbasis, radius*rand(); kwargs...)
            if ϵ < abstol && accept_solution(eq, x, y; abstol)
                if verbose printstyled("$i, $j\n"; color=:yellow) end
                return y, 0, ϵ
            else
                ϵ₀ = min(ϵ, ϵ₀)
                y₀ = y
            end
        end
    end

    if accept_solution(eq, x, y₀; abstol)
        if verbose printstyled("rescue\n"; color=:yellow) end
        return y₀, 0, ϵ₀
    else
        return 0, eq, ϵ₀
    end
end

rms(x) = sqrt(sum(x.^2) / length(x))

"""
    the core of the randomized parameter-fitting algorithm

    `try_integrate` tries to find a linear combination of the basis, whose
    derivative is equal to eq

    output
    -------
    integral, error
"""
function try_integrate(T, eq, x, basis, Δbasis, radius; kwargs...)
    args = Dict(kwargs)
    abstol, opt, attempt_ratio = args[:abstol], args[:opt], args[:attempt_ratio]

    n = length(basis)
    # A is an nxn matrix holding the values of the fragments at n random points
    # b hold the value of the input function at those points
    A = zeros(T, (n, n))
    b = zeros(T, n)

    i = 1
    k = 1

    while i <= n
        if T isa Complex
            x₀ = T(radius*randn()*cis(2π*rand()))
        else
            x₀ = T(radius*rand())
        end
        d = Dict(x => x₀)
        try
            for j = 1:n
                A[i, j] = T(substitute(Δbasis[j], d))
            end
            b[i] = T(substitute(eq, d))
            i += 1
        catch e
            println("basis matrix error: ", e)
        end
        if k > attempt_ratio*n return nothing, 1e6 end
        k += 1
    end

    # find a linearly independent subset of the basis
    l = find_independent_subset(A; abstol)
    A, b, basis, Δbasis, n = A[l,l], b[l], basis[l], Δbasis[l], sum(l)

    if det(A) ≈ 0 return nothing, 1e6 end

    coefs = ones(T, n)
    for j = 1:n
        coefs[j] = coef(Δbasis[j], x)
        A[:,j] /= coefs[j]
    end

    # q₀ = A \ b
    q₀ = Optimize.init(opt, A, b)
    @views Optimize.sparse_regression!(q₀, A, permutedims(b)', opt, maxiter = 1000)
    ϵ = rms(A * q₀ - b)
    q = nice_parameter.(q₀ ./ coefs)
    sum(q[i]*basis[i] for i = 1:length(basis) if q[i] != 0; init=zero(x)), abs(ϵ)
end

function find_independent_basis(T, eq, x, Δbasis; kwargs...)
    args = Dict(kwargs)
    abstol, opt, radius, attempt_ratio = args[:abstol], args[:opt], args[:radius], args[:attempt_ratio]

    n = length(Δbasis)
    A = zeros(T, (n, n))

    i = 1
    k = 1

    while i <= n
        if T isa Complex
            x₀ = T(radius*randn()*cis(2π*rand()))
        else
            x₀ = T(radius*rand())
        end
        d = Dict(x => x₀)
        try
            for j = 1:n
                A[i, j] = T(substitute(Δbasis[j], d))
            end
            i += 1
        catch e
            println(e)
        end
        if k > attempt_ratio*n return nothing, 1e6 end
        k += 1
    end

    # find a linearly independent subset of the basis
    find_independent_subset(A; abstol)
end

"""
    returns a list of the indices of a linearly independent subset of the columns of A
"""
function find_independent_subset(A; abstol=1e-3)
    Q, R = qr(A)
    abs.(diag(R)) .> abstol
end

function find_independent_subset2(A; abstol=1e-3)
    n = size(A, 1)
    l = BitVector(undef, n)
    for i = 1:n
        l[i] = 1
        if abs(det(A[l,l])) < abstol
            l[i] = 0
        end
    end
    l
end

"""
    converts float to int or small rational numbers
"""
function nice_parameters(p; abstol=1e-3)
    c = lcm(collect(1:10)...)
    n = length(p)
    q = Array{Any}(undef, n)
    for i = 1:n
        den = 1
        while den < 10
            if abs(round(p[i]*den) - p[i]*den) < abstol
                a = round(Int, p[i]*den) // den
                q[i] = (denominator(a) == 1 ? numerator(a) : a)
                den = 10
            else
                q[i] = Float64(p[i])
            end
            den += 1
        end
    end
    q
end

function nice_parameter(u::T; abstol=1e-3, M=10) where T<:Real
    c = lcm(collect(1:M)...)
    for den = 1:M
        try
            if abs(round(u*den) - u*den) < abstol
                a = round(Int, u*den) // den
                return (denominator(a) == 1 ? numerator(a) : a)
            end
        catch e
        end
    end
    return u
end

function nice_parameter(u::Complex{T}; abstol=1e-3, M=10) where T<:Real
    α = nice_parameter(real(u))
    β = nice_parameter(imag(u))
    return β ≈ 0 ? α : Complex(α, β)
end

function try_integration_by_parts(eq, x; kwargs...)
    f = factors(eq, x)
    if length(f) <= 2 return zero(x), eq, Inf end

    D = Differential(x)
    ϵ₀ = Inf

    for u in f
        v′ = eq / u
        if !is_number(u) && !is_number(v′)
            v, r, ϵ = integrate_term(v′, x; kwargs...)
            if isequal(r, 0)
                uv = expand_derivatives(v*D(u))
                s, r, ϵ = integrate_sum(uv, x; kwargs...)
                if isequal(r, 0)
                    return expand(u*v - s), 0, ϵ
                else
                    zero(x), eq, ϵ
                    # ϵ₀ = min(ϵ, ϵ₀)
                end
            end
        end
    end

    return zero(x), eq, ϵ₀
end

var_index(v) = istree(v) && Symbol(operation(v)) == :getindex ? arguments(v)[2] : -1

function try_symbolic(T, eq, x, basis, Δbasis; kwargs...)
    n = length(basis)
    @syms θ[1:n]
    D = Differential(x)
    Δeq = expand(sum(θ[j]*Δbasis[j] for j=1:n) - eq)

    terms = collect_terms(Δeq, x)
    eqs = collect(values(terms))
    # eqs = filter(p->length(get_variables(p)) > 0, eqs)

    sol = solve_symbolic(eqs)

    for i = 1:n
        if !haskey(sol, θ[i])
            sol[θ[i]] = 0
        end
    end

    p = substitute(expand(sum(θ[j]*basis[j] for j=1:n)), sol)
    p, 0
end

mutable struct Fragment
    eq
    lhs
end

function solve_symbolic(eqs)
    n = length(eqs)
    solved = Set()
    unfinished = true
    frags = Fragment[]
    k = 1

    while unfinished
        unfinished = false
        for eq in eqs
            δf = [v for v in get_variables(eq) if v ∉ solved]
            if length(δf) == 1
                push!(solved, δf[1])
                push!(frags, Fragment(eq, δf[1]))
                unfinished = true
            end
        end
    end

    sol = Dict()
    for f in frags
        eq = substitute(f.eq, sol)
        u = Symbolics.solve_for(eq ~ 0, f.lhs)
        sol[f.lhs] = nice_parameter(u)
    end

    sys = []
    vars = Set()

    for eq in eqs
        δf = [v for v in get_variables(eq) if v ∉ solved]
        if length(δf) > 1
            for v in δf
                push!(vars, v)
            end
            q = substitute(eq, sol)
            push!(sys, q)
        end
    end

    sys = unique(sys)
    vars = [v for v in vars]

    if !isempty(vars) && length(vars) == length(sys)
        try
            vals = nice_parameters.(Symbolics.solve_for(sys .~ 0, vars))
            # vals = Symbolics.solve_for(sys .~ 0, vars)
            for (v,u) in zip(vars, vals)
                sol[v] = u
            end
        catch e
            # println("from symbolic: ", e)
        end
    end

    sol
end

########################## Transformation Rules ###############################

trig_rule1 = @rule tan(~x) => sin(~x) / cos(~x)
trig_rule2 = @rule sec(~x) => one(~x) / cos(~x)
# trig_rule2 = @rule sec(~x) => (tan(~x)/cos(~x) + 1/cos(~x)^2) / (tan(~x) + 1/cos(~x))
trig_rule3 = @rule csc(~x) => one(~x) / sin(~x)
# trig_rule3 = @rule csc(~x) => (cot(~x)/sin(~x) + 1/sin(~x)^2) / (cot(~x) + 1/sin(~x))
trig_rule4 = @rule cot(~x) => cos(~x) / sin(~x)

trig_rules = [trig_rule1, trig_rule2, trig_rule3, trig_rule4]

hyper_rule1 = @rule tanh(~x) => sinh(~x) / cosh(~x)
hyper_rule2 = @rule sech(~x) => one(~x) / cosh(~x)
hyper_rule3 = @rule csch(~x) => one(~x) / sinh(~x)
hyper_rule4 = @rule coth(~x) => cosh(~x) / sinh(~x)

hyper_rules = [hyper_rule1, hyper_rule2, hyper_rule3, hyper_rule4]

misc_rule1 = @rule sqrt(~x) => ^(~x, 0.5)

misc_rules = [misc_rule1]

int_rules = [trig_rules; hyper_rules; misc_rules]
# int_rules = misc_rules

apply_integration_rules(eq) = Fixpoint(Prewalk(PassThrough(Chain(int_rules))))(value(eq))

function U(u...)
    u = map(x -> x isa AbstractArray ? x : [], u)
    return union(u...)
end

hints(eq::SymbolicUtils.Add, x, h) = map(t->hints(t,x,h), arguments(eq))
hints(eq::SymbolicUtils.Mul, x, h) = map(t->hints(t,x,h), arguments(eq))
hints(eq::SymbolicUtils.Pow, x, h) = hints(arguments(eq)[1],x,h)

function hints(eq::SymbolicUtils.Term, x, h)
    s = Symbol(operation(eq))
    u = arguments(eq)[1]

    if s == :sec
        push!(h, log(1/cos(u) + sin(u)/cos(u)))
    elseif s == :csc
        push!(h, log(1/sin(u) - cos(u)/sin(u)))
    elseif s == :tan
        push!(h, log(cos(u)))
    elseif s == :cot
        push!(h, log(sin(u)))
    elseif s == :tanh
        push!(h, log(cosh(u)))
    end
end

function hints(eq, x, h)
end

function collect_hints(eq, x)
    h = []
    hints(eq, x, h)
    h
end

########################## Convert to a Polynomial? #############################

@polyvar 𝑦

is_number(x::T) where T<:Integer = true
is_number(x::T) where T<:Real = true
is_number(x::T) where T<:Complex = true
is_number(x::T) where T<:Rational = true
is_number(x) = false

function is_poly(eq, x)
    p = collect_powers(value(eq), x)
    all(is_number, values(p))
end

function to_poly(p::SymbolicUtils.Add, x)
    l = [to_poly(t,x) for t in arguments(p)]
    if any(x->x==nothing, l) return nothing end
    sum(l; init=0)
end

function to_poly(p::SymbolicUtils.Mul, x)
    l = [to_poly(t,x) for t in arguments(p)]
    if any(x->x==nothing, l) return nothing end
    prod(l; init=1)
end

function to_poly(p::SymbolicUtils.Pow, x)
    if !isequal(arguments(p)[1], x) return nothing end
    return 𝑦^arguments(p)[2]
end

to_poly(p::SymbolicUtils.Sym, x) = isequal(p, x) ? 𝑦 : nothing
to_poly(p::SymbolicUtils.Term, x) = nothing
to_poly(p, x) = is_number(p) ? p : nothing

function to_rational(p::SymbolicUtils.Mul, x)
    P = one(x)
    Q = one(x)

    for t in arguments(p)
        y = to_poly(t, x)
        if y != nothing
            P *= y
        else
            y = to_poly(1/t, x)
            if y != nothing
                Q *= y
            else
                return nothing
            end
        end
    end

    if is_number(P) || is_number(Q) return nothing end
    return sym(factor_rational(P, Q), x)
end

to_rational(p, x) = nothing

is_multiple_x(p::SymbolicUtils.Add, x) = all(is_multiple_x(t,x) for t in arguments(p))
is_multiple_x(p::SymbolicUtils.Mul, x) = any(is_multiple_x(t,x) for t in arguments(p))
is_multiple_x(p::SymbolicUtils.Pow, x) = isequal(arguments(p)[1], x) && arguments(p)[2] >= 1
is_multiple_x(p::SymbolicUtils.Sym, x) = isequal(p, x)
is_multiple_x(p, x) = false

factors(eq, x) = isdependent(eq, x) ? [one(x), eq] : [one(x)]

function factors(eq::SymbolicUtils.Pow, x)
    p, k = arguments(eq)
    [p^(i*sign(k)) for i=0:abs(k)]
end

function factors(eq::SymbolicUtils.Mul, x)
    terms = [factors(q,x) for q in arguments(eq)]
    n = length(terms)

    l = Any[one(x)]

    for j = 1:n
        m = length(l)
        for t in terms[j]
            for k = 1:m
                push!(l, l[k]*t)
            end
        end
    end

    unique(l)
end

coef(eq::SymbolicUtils.Mul, x) = prod(t for t in arguments(eq) if !isdependent(t,x); init=1)
coef(eq::SymbolicUtils.Add, x) = minimum(abs(coef(t,x)) for t in arguments(eq))
coef(eq, x) = 1

function collect_terms(eq::SymbolicUtils.Add, x)
    d = Dict{Any,Any}(1 => 0)
    for t in arguments(eq)
        if isdependent(t, x)
            for s in collect_terms(t, x)
                v, c = first(s), last(s)
                if haskey(d, v)
                    d[v] += c
                else
                    d[v] = c
                end
            end
        else
            d[1] += t
        end
    end
    d
end

function collect_terms(eq, x)
    c = coef(eq, x)
    return Dict(eq/c => c)
end

##############################################################################

@syms x β

"""
    a list of basic standard integral tests
    based on http://integral-table.com/ with modifications
"""
basic_integrals = [
# Basic Forms
    1,
    x^2,
    4x^3,
# Integrals of Rational Functions
    1/x,
    1/(2x + 5),
    1/(x + 1)^2,
    (x + 3)^3,
    x*(x - 2)^4,
    1/(1 + x^2),
    1/(9 + x^2),
    x/(4 + x^2),
    x^2/(16 + x^2),
    x^3/(1 + x^2),
    1/(x^2 - 5x + 6),
    1/(x^2 + x + 1),
    x/(x + 4)^2,
    x/(x^2 + x + 1),
# Integrals with Roots
    sqrt(x - 2),
    1 / sqrt(x - 1),
    1 / sqrt(x + 1),
    1 / sqrt(4 - x),
    x * sqrt(x - 3),
    sqrt(2x + 5),
    (3x - 1)^1.5,
    x / sqrt(x - 1),
    x / sqrt(x + 1),
    sqrt(x / (4 - x)),
    sqrt(x / (4 + x)),
    x * sqrt(2x + 3),
    sqrt(x*(x+2)),
    sqrt(x^3*(x+3)),
    sqrt(x^2 + 4),
    sqrt(x^2 - 4),
    sqrt(4 - x^2),
    x * sqrt(x^2 + 9),
    x * sqrt(x^2 - 9),
    1 / sqrt(x^2 + 4),
    1 / sqrt(x^2 - 4),
    1 / sqrt(4 - x^2),
    x / sqrt(x^2 + 4),
    x / sqrt(x^2 - 4),
    x / sqrt(4 - x^2),
    x^2 / sqrt(x^2 + 4),
    x^2 / sqrt(x^2 - 4),
    sqrt(x^2 - 5x + 6),
    x * sqrt(x^2 - 5x + 6),
    1 / sqrt(x^2 - 5x + 6),
    1 / (4 + x^2)^1.5,
# Integrals with Logarithms
    log(x),
    x * log(x),
    x^2 * log(x),
    log(2x) / x,
    log(x) / x^2,
    log(2x + 1),
    log(x^2 + 4),
    log(x^2 - 4),
    log(x^2 - 5x + 6),
    x * log(x + 2),
    x * log(9 - 4x^2),
    log(x)^2,
    log(x)^3,
    x * log(x)^2,
    x^2 * log(x)^2,
# Integrals with Exponentials
    exp(x),
    sqrt(x) * exp(x),
    x * exp(x),
    x * exp(3x),
    x^2 * exp(x),
    x^2 * exp(5x),
    x^3 * exp(x),
    x^3 * exp(2x),
    exp(x^2),
    x * exp(x^2),
# Integrals with Trigonometric Functions
    sin(4x),
    sin(x)^2,
    sin(x)^3,
    cos(3x),
    cos(x)^2,
    cos(2x)^3,
    sin(x) * cos(x),
    sin(3x) * cos(5x),
    sin(x)^2 * cos(x),
    sin(3x)^2 * cos(x),
    sin(x) * cos(x)^2,
    sin(x) * cos(5x)^2,
    sin(x)^2 * cos(x),
    sin(x)^2 * cos(x)^2,
    sin(4x)^2 * cos(4x^2),
    tan(x),
    tan(7x),
    tan(x)^2,
    tan(x)^3,
    sec(x),
    sec(x) * tan(x),
    sec(x)^2 * tan(x),
    csc(x),
    sec(x) * csc(x),
# Products of Trigonometric Functions and Monomials
    x * cos(x),
    x * cos(3x),
    x^2 * cos(x),
    x^2 * cos(5x),
    x * sin(x),
    x * sin(3x),
    x^2 * sin(x),
    x^2 * sin(5x),
    x * cos(x)^2,
    x * sin(x)^2,
    x * tan(x)^2,
    x * sec(x)^2,
    x^3 * sin(x),
    x^4 * cos(2x),
    sin(x)^2 * cos(x)^3,
# Products of Trigonometric Functions and Exponentials
    exp(x) * sin(x),
    exp(3x) * sin(2x),
    exp(x) * cos(x),
    exp(2x) * cos(7x),
    x * exp(x) * sin(x),
    x * exp(x) * cos(x),
# Integrals of Hyperbolic Functions
    cosh(x),
    exp(x) * cosh(x),
    sinh(3x),
    exp(2x) * sinh(3x),
    tanh(x),
    exp(x) * tanh(x),
    cos(x) * cosh(x),
    cos(x) * sinh(x),
    sin(x) * cosh(x),
    sin(x) * sinh(x),
    sinh(x) * cosh(x),
    sinh(3x) * cosh(5x),
# Misc
    exp(x) / (1 + exp(x)),
    cos(exp(x)) * sin(exp(x)) * exp(x),
    cos(exp(x))^2 * sin(exp(x)) * exp(x),
    1 / (x*log(x)),
    (log(x - 1) + (x - 1)^-1) * log(x),
    1 / (exp(x) - 1),
    1 / (exp(x) + 5),
    sqrt(x)*log(x),
    log(log(x)) / x,
    x^3 * exp(x^2),
    sin(log(x)),
    x * cos(x) * exp(x),
    log(x - 1)^2,
    1 / (exp(2x) - 1),
    exp(x) / (exp(2x) - 1),
    x / (exp(2x) - 1),
# derivative-divide examples (Lamangna 7.10.2)
    exp(x) * exp(exp(x)),
    exp(sqrt(x)) / sqrt(x),
    log(log(x)) / (x*log(x)),
    log(cos(x)) * tan(x),
# rothstein-Trager examples (Lamangna 7.10.9)
    1 / (x^3 - x),
    1 / (x^3 + 1),
    1 / (x^2 - 8),
    (x + 1) / (x^2 + 1),
    x / (x^4 - 4),
    x^3 / (x^4 + 1),
    1 / (x^4 + 1),
# bypass = true
    β,      # turn of bypass = true
    exp(x)/x - exp(x)/x^2,
    cos(x)/x - sin(x)/x^2,
    1/log(x) - 1/log(x)^2,
]

function test_integrals(; symbolic=false, num_trials=3, verbose=false, bypart=true)
    bypass = false
    misses = []
    D = Differential(x)

    for eq in basic_integrals
        if isequal(eq, β)
            printstyled("**** bypass on ****\n"; color=:red)
            bypass = true
        else
            printstyled(eq, " =>\t"; color=:green)
            solved, unsolved = integrate(eq; bypass, symbolic, num_trials, verbose, bypart)
            printstyled(solved; color=:white)
            if isequal(unsolved, 0)
                println()
            else
                printstyled(" + ∫ ", unsolved, '\n'; color=:red)
                push!(misses, eq)
            end
        end
    end

    n = length(misses)
    if n > 0 println("**** missess (n=$n) *****") end
    for eq in misses
        printstyled(eq, '\n'; color=:red)
    end
end

##################### Special Functions ######################################

# """
#     logarithmic integral
# """
# function li(x; n=10)
#     z = log(abs(x))
#     s = sum(z^k / (factorial(k) * k) for k = 1:n)
#     return SpecialFunctions.γ + log(z) + s
# end


text:
function contract!(tiJaB_d_temp3, tiJaB_i, Wmnij)
  rvir = axes(tiJaB_d_temp3, 4)
  nvir = last(rvir)
  rocc = axes(tiJaB_d_temp3, 1)
  @inbounds @fastmath for b in rvir, a = b:nvir, j in rocc, i = j:last(rocc)
    temp = zero(eltype(tiJaB_i))
    for n in rocc, m in rocc
      temp += tiJaB_i[m, n, a, b] * Wmnij[m, n, i, j]
    end
    tiJaB_d_temp3[i, j, a, b] = temp
    tiJaB_d_temp3[j, i, b, a] = temp
  end
end

function contracttest1!(tiJaB_d_temp3, tiJaB_i, Wmnij)
  rvir = axes(tiJaB_d_temp3, 4)
  nvir = last(rvir)
  rocc = axes(tiJaB_d_temp3, 1)
  for b in rvir, j in rocc
    @turbo for a = b:nvir, i = j:last(rocc)
      temp = zero(eltype(tiJaB_i))
      for n in rocc, m in rocc
        temp += tiJaB_i[m, n, a, b] * Wmnij[m, n, i, j]
      end
      tiJaB_d_temp3[i, j, a, b] = temp
      tiJaB_d_temp3[j, i, b, a] = temp
    end
  end
  tiJaB_d_temp3
end
function contracttest2!(tiJaB_d_temp3, tiJaB_i, Wmnij)
  rvir = axes(tiJaB_d_temp3, 4)
  nvir = last(rvir)
  rocc = axes(tiJaB_d_temp3, 1)
  for b in rvir, a = b:nvir, j in rocc
    @turbo for i = j:last(rocc)
      temp = zero(eltype(tiJaB_i))
      for n in rocc, m in rocc
        temp += tiJaB_i[m, n, a, b] * Wmnij[m, n, i, j]
      end
      tiJaB_d_temp3[i, j, a, b] = temp
      tiJaB_d_temp3[j, i, b, a] = temp
    end
  end
  tiJaB_d_temp3
end

@testset "Tensors" begin
  @show Float64, @__LINE__
  LA, LIM = 31, 23
  A = rand(LIM, LIM, LA, LA)
  B = rand(LIM, LIM, LIM, LIM)

  C1 = fill(-999.999, LIM, LIM, LA, LA)
  C2 = fill(-999.999, LIM, LIM, LA, LA)
  C3 = fill(-999.999, LIM, LIM, LA, LA)
  # C1 = Array{Float64}(undef, LIM, LIM, LA, LA);
  # C2 = similar(C1); C3 = similar(C1);

  contract!(C1, A, B)
  @test C1 ≈ contracttest1!(C2, A, B)
  @test C1 ≈ contracttest2!(C3, A, B)

  Apermute = PermutedDimsArray(permutedims(A, (2, 4, 1, 3)), (3, 1, 4, 2))
  Bpermute = PermutedDimsArray(permutedims(B, (2, 1, 4, 3)), (2, 1, 4, 3))
  @test C1 ≈ contracttest1!(C2, Apermute, Bpermute)
  @test C1 ≈ contracttest2!(C3, Apermute, Bpermute)
end


text:
using SymEngine
import YaoBase: rot_mat

rot_mat(::Type{T}, ::Val{:Rx}, theta::Basic) where T =
    Basic[cos(theta/2) -im * sin(theta/2); -im * sin(theta/2) cos(theta/2)]
rot_mat(::Type{T}, ::Val{:Ry}, theta::Basic) where T =
    Basic[cos(theta/2) -sin(theta/2); sin(theta/2) cos(theta/2)]
rot_mat(::Type{T}, ::Val{:Rz}, theta::Basic) where T =
    Diagonal(Basic[exp(-im*theta/2), exp(im*theta/2)])
rot_mat(::Type{T}, ::Val{:CPHASE}, theta::Basic) where T =
    Diagonal(Basic[1, 1, 1, exp(im*theta)])
rot_mat(::Type{T}, ::Val{:PSWAP}, theta::Basic) where T =
    rot_mat(Basic, Const.SWAP, theta)

for G in [:Rx, :Ry, :Rz, :CPHASE]
    # forward single gates
    @eval function YaoBase.instruct!(state::AbstractVecOrMat{T}, g::Val{$(QuoteNode(G))},
                            locs::Union{Int, NTuple{N3,Int}},
                            control_locs::NTuple{N1, Int},
                            control_bits::NTuple{N2, Int}, theta::Basic) where {T, N1, N2, N3}
        m = rot_mat(T, g, theta)
        instruct!(state, m, locs, control_locs, control_bits)
    end
end

# forward single gates
@eval function YaoBase.instruct!(state::AbstractVecOrMat{T}, g::Val,
                        locs::Union{Int, NTuple{N1, Int}}, theta::Basic) where {T, N1}
    instruct!(state, g, locs, (), (), theta)
end

@eval function YaoBase.instruct!(state::AbstractVecOrMat{T}, g::Val{:PSWAP},
                        locs::Tuple{Int, Int}, theta::Basic) where {T, N1}
    instruct!(state, g, locs, (), (), theta)
end


text:
type ZVSample
  samples::Matrix{Float64}
  gradients::Matrix{Float64}
  diagnostics::Dict

  function ZVSample(samples::Matrix{Float64}, gradients::Matrix{Float64}, diagnostics::Dict)
    @assert size(samples) == size(gradients) "Samples and gradients must be of the same size"
    new(samples, gradients, diagnostics)
  end
end

ZVSample(samples::Matrix{Float64}, gradients::Matrix{Float64}) = ZVSample(samples, gradients, Dict())

ZVSample(chain::MCMCChain) = ZVSample(chain.samples, chain.gradients, chain.diagnostics)

type ZVMean
  means::Matrix{Float64}
  diagnostics::Dict
end

ZVMean(means::Matrix{Float64}) = ZVMean(means, Dict())


text:
using Test
using Petri
using LabelledArrays
using OrdinaryDiffEq
using SteadyStateDiffEq
using StochasticDiffEq
using DiffEqJump
using Random

Random.seed!(1234);

@testset "Generation of ODE formulas" begin
    sir = Petri.Model([:S,:I,:R],[(LVector(S=1,I=1), LVector(I=2)),
                                  (LVector(I=1),     LVector(R=1))])

    sirf = vectorfield(sir)
    rates = sirf(LVector(S=0.0, I=0.0, R=0.0), LVector(S=100.0,I=1.0,R=0.0), [0.5/101,0.25],0.0)
    @test rates.S ≈ -.495 atol=1e-2
    @test rates.I ≈ .24505 atol=1e-2
    @test rates.R ≈ .25 atol=1e-2
end

@testset "Generation of ODE solutions" begin
    @testset "SIR Solving" begin
        sir = Petri.Model([:S,:I,:R],[(LVector(S=1,I=1), LVector(I=2)),
                                      (LVector(I=1),     LVector(R=1))])
        u0 = LVector(S=990.0,I=10.0,R=0.0)
        p = [(u,t)->0.5/sum(u), 0.25]
        prob = ODEProblem(sir,u0,(0.0,40.0),p)
        sol = OrdinaryDiffEq.solve(prob,Tsit5())
        @test sum(sol[end]) ≈ 1000 atol=1
        @test sol[end].S ≈ 210 atol=1
        @test sol[end].I ≈ 14.5 atol=1
        @test sol[end].R ≈ 775.5 atol=1
    end
end

@testset "Generation of SDE solutions" begin
    @testset "SIR Solving" begin
        sir = Petri.Model([:S,:I,:R],[(LVector(S=1,I=1), LVector(I=2)),
                                      (LVector(I=1),     LVector(R=1))])
        u0 = LVector(S=990.0,I=10.0,R=0.0)
        p = [(u,t)->(0.5/sum(u)), 0.25]
        prob,cb = SDEProblem(sir,u0,(0.0,40.0),p)
        sol = StochasticDiffEq.solve(prob,SRA1(),callback=cb)
        @test sum(sol[end]) ≈ 1000 atol=1
        @test sol[end].S ≈ 201 atol=1
        @test sol[end].I ≈ 13 atol=1
        @test sol[end].R ≈ 786 atol=1
    end
end

@testset "Generation of Jump solutions" begin
    @testset "SIR Solving" begin
        sir = Petri.Model([:S,:I,:R],[(LVector(S=1,I=1), LVector(I=2)),
                                      (LVector(I=1),     LVector(R=1))])
        u0 = LVector(S=990.0,I=10.0,R=0.0)
        p = [0.5/sum(u0), 0.25]
        prob = JumpProblem(sir,u0,(0.0,40.0),p)
        sol = DiffEqJump.solve(prob,SSAStepper())
        @test sum(sol[end]) == 1000
        @test sol[end].S == 236
        @test sol[end].I == 45
        @test sol[end].R == 719
    end
end

@testset "Generation of SteadyState solutions" begin
    @testset "SIR Steady State Solving" begin
        sir = Petri.Model([:S,:I,:R],[(LVector(S=1,I=1), LVector(I=2)),
                                      (LVector(I=1),     LVector(R=1))])
        u0 = LVector(S=990.0,I=0.0,R=0.0)
        p = [(u,t)->0.5/sum(u), 0.25]
        prob = SteadyStateProblem(sir,u0,(0.0,40.0),p)
        sol = OrdinaryDiffEq.solve(prob,DynamicSS(Tsit5()))
        @test sum(sol.u) ≈ 990 atol=1e-2
        @test sol.u.S ≈ 990 atol=1e-2
        @test sol.u.I ≈ 0 atol=1e-2
        @test sol.u.R ≈ 0 atol=1e-2
    end
end


text:
using DelayedKalmanFilter
using DelimitedFiles
using Test

@testset "ODE tests" begin
    steady_state = calculate_steady_state_of_ode(10000.,5.,log(2)/30,log(2)/90,1.,1.);

    @test length(steady_state) == 2
    @test steady_state[1] ≈ 41.41785876041631
    @test steady_state[2] ≈ 5377.800549410291
end #testset