text:
# Copyright 2020 Digital Domain 3.0
#
# Licensed under the Apache License, Version 2.0 (the "Apache License")
# with the following modification; you may not use this file except in
# compliance with the Apache License and the following modification to it:
# Section 6. Trademarks. is deleted and replaced with:
#
# 6. Trademarks. This License does not grant permission to use the trade
#    names, trademarks, service marks, or product names of the Licensor
#    and its affiliates, except as required to comply with Section 4(c) of
#    the License and to reproduce the content of the NOTICE file.
#
# You may obtain a copy of the Apache License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the Apache License with the above modification is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Apache License for the specific
# language governing permissions and limitations under the Apache License.

"""
Extension of Blades into a vector space.

KVectors are containers for Blades of the same grade.  New Blade elements are added with the + operator
"""

export
KVector,
NullKVector,
sort_basis,
isnull,
coords,
prune,
normalize_safe,
basis_1vector,
reciprocal_frame

include("Blades.jl")
using StaticArrays
using LinearAlgebra

"Container of Blades with the same grade and element type"
struct KVector{T,K,N} <: AbstractArray{T,1}
  k::SVector{N,Blade{T,K}}
end

Base.getindex(B::BT, i::T) where {BT<:KVector, T<:Integer} = B.k[i]
Base.firstindex(k::K) where {K<:KVector} = firstindex(k.k)
Base.lastindex(k::K) where {K<:KVector} = lastindex(k.k)

const NullKVector{T,K} = KVector{T,K,0}

NullKVector{T,K}() where {T,K} = KVector{T,K,0}(SVector{0,T}())
isnull(b::B) where {B<:NullKVector} = true
isnull(b::B) where {B<:KVector} = false

Base.first(k::K) where {K<:KVector} = isempty(k) ? nothing : k.k[1]
Base.length(k::K) where {K<:KVector} = length(k.k)
Base.size(k::K) where {K<:KVector} = (length(k.k),)
Base.isempty(k::K) where {K<:KVector} = isempty(k.k)
Base.iterate(A::K) where {K<:KVector} = isempty(A) ? nothing : (A[1], 2)
Base.iterate(A::K, i::T) where {K<:KVector, T} = length(A) < i ? nothing : (A[i], i+1)
Base.eltype(::Type{K}) where {T,P,K<:KVector{T,P}} = Blade{T,P}

Base.conj(k::K) where {K<:KVector} = K(map(conj, k))

KVector{T}(v::Blade{T,K}) where {T,K} = KVector{T,K,1}(SArray{Tuple{1},Blade{T,K},1,1}(v))
KVector(v::Blade{T,K}) where {T,K} = KVector{T}(v)
Base.convert(::Type{B}, k::K) where {T,N,B<:KVector, K<:Blade{T,N}} = B(k)
KVector(kv::V) where {T, KN, K<:Blade{T,KN}, V<:AbstractVector{K}} = KVector{T,KN,length(kv)}(kv)
KVector(kv::KV) where KV<:KVector = kv

KVector(uv::C, 𝐼::Blade{KT,N}) where {T, N, KT, C<:SVector{N,T}} = KVector(uv.*(basis_1blades(𝐼)))

KVector(s::T) where T<:Real = s

function KVector(uv::C, ns::Module) where {T, N, C<:SVector{N,T}}
  KVector(uv.*(basis_1blades(pseudoscalar(ns))))
end

function KVector(uv::C, ns::Bool = false) where {T, N, C<:SVector{N,T}}
  @warn "resolving Blade types in top-level module"
  KVector(uv, Main)
end

KVector(uv::A, ns::MN=false) where {T<:Real, A<:Vector{T}, MN<:Union{Bool, Module}} = KVector(SVector{length(uv)}(uv), ns)

KVector(uv::A, 𝐼::Blade{KT,N}) where
  {T<:Real, A<:Vector{T}, N, KT} =
  KVector(SVector{length(uv)}(uv), 𝐼)

KVector(uv, ::Type{I}) where {I<:Blade} = KVector(uv, 1I)

vector(k::KA) where {K<:KVector, KA<:AbstractVector{K}} = SVector{3}.(coords.(k))

Base.one(::Type{KVector{T,K,N}}) where {T,K,N} = KVector(normalize(one(T) .* (basis_kblades(dual(1),K))))
Base.one(k::K) where K<:KVector = one(K)

Base.iszero(b::B) where {T,K,N,B<:KVector{T,K,0}} = true
Base.iszero(b::B) where {T,K,N,B<:KVector{T,K,N}} = iszero(sum((x->x*x).((k->k.x).(b.k))))

Base.:+(b::KA, c::KB) where {T<:Number,K, KA<:Blade{T,K}, KB<:Blade{T,K}} =
  KVector{T,K,2}(SArray{Tuple{2},Blade{T,K},1,2}(b,c))
Base.:-(b::KA, c::KB) where {T<:Number,K, KA<:Blade{T,K}, KB<:Blade{T,K}} =
  KVector{T,K,2}(SArray{Tuple{2},Blade{T,K},1,2}(b,-c))

function Base.:+(b::V, k::KVector{TK,K,N}) where {TK,TB,K, V<:Blade{TB,K},N}
  T = promote_type(TK, TB)
  i = findfirst(v->typeof(v)==V, k.k)
  if i == nothing
    #KVector{T,K,N+1}(vcat(Blade{T,K}[T(scalar(b))*untype(b)],k.k))
    KVector{T,K,N+1}(vcat(Blade{T,K}[b],k.k))
  else
    KVector(setindex(k.k,k.k[i]+b,i))
  end
end

Base.:+(k::KVector{TK,K,N}, b::V) where {TK,TB,K, V<:Blade{TB,K},N} = b+k
Base.:+(b::B, k::V) where {T,K,B<:KVector{T,K},V<:Blade{T,K}} = k+b
Base.:+(b::KVector{T,K}, c::KVector{T,K}) where {T,K} = reduce( (bc,v)->v+bc, b.k; init=c )
Base.:-(b::B, k::V) where {T,K,B<:KVector{T,K},V<:Blade{T,K}} = -k+b
Base.:-(k::V, b::B) where {T,K,B<:KVector{T,K},V<:Blade{T,K}} = k+(-b)
Base.:-(a::B) where {T, B<:KVector{T}} = -one(T)*a
Base.:-(b::KVector{T,K}, c::KVector{T,K}) where {T,K} = b+(-c)

"apply reverse to all Blades in this KVector"
Base.reverse(b::B) where B<:KVector = B(reverse.(b.k))

"""
    grade(a)

the grade k of a k-vector.  roughly speaking, the dimension of the subspace spanned by the k-vector
"""
grade(b::KVector{T,K}) where {T,K} = K
grade(b::KVector{T,K}, i::Integer) where {T,K} = i==K ? b : zero(T)

"apply dual to all Blades in this KVector"
dual(b::KVector{T,K}) where {T,K} = K==grade(pseudoscalar(first(b))) ? b.k[1].x : KVector(dual.(b))
dual(b::KVector{T,K,0}) where {T,K} = b
"apply ⟂ to all Blades in this KVector"
⟂(b::KVector{T,K})  where {T,K} = KVector((⟂).(b))
⟂(b::KVector{T,K,0})  where {T,K} = b
"apply ! to all Blades in this KVector"
Base.:!(b::B) where {B<:KVector} = dual(b)

#rc(b::KVector{T,K}) where {T,K} = K==grade(pseudoscalar(first(b))) ? b.k[1].x : KVector(rc.(b), parentmodule(typeof(first(b))))
rc(b::KVector{T,K}) where {T,K} = K==grade(pseudoscalar(first(b))) ? b.k[1].x : sum(rc.(b))
rc(b::KVector{T,K,0}) where {T,K} = b
#lc(b::KVector{T,K}) where {T,K} = K==grade(pseudoscalar(first(b))) ? b.k[1].x : KVector(lc.(b), parentmodule(typeof(first(b))))
lc(b::KVector{T,K}) where {T,K} = K==grade(pseudoscalar(first(b))) ? b.k[1].x : sum(lc.(b))
lc(b::KVector{T,K,0}) where {T,K} = b

Base.zero(b::Type{B}) where {T,K,N,B<:KVector{T,K,N}} = NullKVector{T,K}()
Base.zero(b::KVector) = zero(typeof(b))

Base.:*(nb::NullKVector{T,K}, nb2::NullKVector{T,K}) where {T,K} = nb
Base.:*(nb::NullKVector{T,K}, b::KVector{T,K}) where {T,K} = nb
Base.:*(b::KVector{T,K}, nb::NullKVector{T,K}) where {T,K} = nb*b
Base.:+(nb::NullKVector{T,K}, nb2::NullKVector{T,K}) where {T,K} = nb
Base.:+(nb::NullKVector{T,K}, b::KVector{T,K}) where {T,K} = b
Base.:+(b::KVector{T,K}, nb::NullKVector{T,K}) where {T,K} = nb+b

Base.:*(s::T, b::B) where {T<:Real, B<:KVector} = B(s*b.k)
Base.:*(b::B, s::T) where {T<:Real, B<:KVector} = s*b
Base.:/(b::B, s::T) where {T<:Real, B<:KVector} = b*(one(T)/s)

pseudoscalar(k::K) where {K<:KVector} = pseudoscalar(first(k))

#==
"""
    inv(k)

Inverse of a KVector when it exists.  Left inverse.  k⋅inv(k) = 1
"""
Base.inv(k::K) where {K<:KVector} = reverse(k)/(k⋅reverse(k))
==#

"""
    ⋆(k)

Hodge star operator mapping k to it's Hodge dual.
"""
⋆(k::K) where {K<:KVector} = mapreduce(⋆,+,k)

"""
    sort_basis(b)

sort blades by bases indices in ascending order within the k-vector
"""
sort_basis(B::BT) where {BT<:KVector} = BT(sort(Vector(B.k); by=subspace∘typeof))
sort_basis(B::BT) where {BT<:Blade} = B

Base.:(==)(B::BT, B2::BT2) where {BT<:Union{KVector, Blade}, BT2<:Union{KVector, Blade}} = iszero(B-B2)

Base.promote_rule(::Type{B}, ::Type{K}) where {B<:KVector,K<:Blade} = KVector

prune(k::KVector, epsi = eps()) = sum(Iterators.filter(x->abs(scalar(x)) > epsi, k))
prune(b::Blade, epsi = eps()) = abs(scalar(b)) > epsi ? b : scalar(b)

"""
    ∧(a,b)

Wedge product between two k-vectors
"""
function ∧(A::KV,B::KV2) where {T,K,N,KV<:KVector{T,K,N}, K2,N2,KV2<:KVector{T,K2,N2}}
  AwB = Blade{T, K+K2}[]
  sizehint!(AwB, N*N2)
  for aᵢ in A.k
    for bⱼ in B.k
      awb = aᵢ∧bⱼ
      if !iszero(awb)
        push!(AwB, awb)
      end
    end
  end
  sum(AwB)
end

"""
    ∧(a,b)

Wedge product between blade and k-vector
"""
function ∧(a::BT,B::KV2) where {T,K,BT<:Blade{T,K}, K2,N2,KV2<:KVector{T,K2,N2}}
  AwB = Blade{T, K+K2}[]
  sizehint!(AwB, N2)
  for bⱼ in B.k
    awb = a∧bⱼ
    if !iszero(awb)
      push!(AwB, awb)
    end
  end
  sum(AwB)
end

∧(B::M,a::N) where {M<:KVector,N<:Blade} = -a∧B
∧(s::T,B::N) where {T<:Real,N<:KVector} = s*B
∧(B::N,s::T) where {T<:Real,N<:KVector} = s*B

"""
    gram(a, u)

The Gram matrix, a length(a)xlength(a) matrix where each entry is the dot product aᵢ⋅uⱼ

Entries are sorted by acending basis indices from a and u.

Example:
for length(a) = 2, gram( a, u ) = [ a₁⋅u₁ a₁⋅u₂ ;
                                    a₂⋅u₁ a₂⋅u₂ ]
"""
function gram(a::K,u::L) where {K<:KVector, L<:KVector}
  a = sort_basis(prune(a))
  u = sort_basis(prune(u))
  n = length(a)
  if n != length(u)
    zero(eltype(eltype(a)))
  else
    [aᵢ⋅uⱼ for uⱼ in u for aᵢ in a] |> A->(reshape(A,n,n))
  end
end

import LinearAlgebra: det

"Determinant of Gram matrix formed from k-vectors a and u"
det(a::K,u::L) where {K<:KVector, L<:KVector} = det(gram(a,u))

import LinearAlgebra: norm,normalize,norm_sqr

norm_sqr(k::K) where {K<:KVector} = k⋅reverse(k)
norm(k::K) where {K<:KVector} = sqrt(norm_sqr(k))
normalize(k::K) where {K<:KVector} = k/norm(k)
normalize_safe(k::K) where {K<:KVector} = (kk = norm_sqr(k); abs(kk) > 0 ? k/norm(k) : k)

import Base: cos

"""
    cos(a, b)

cosine of angle between blades a and be spanned by the k-vectors a and b
"""
cos( a::A, b::B ) where { A<:Union{KVector, Blade}, B<:Union{KVector, Blade}} = norm(a⋅b)/(norm(a)*norm(b))

basis_1blades( k::K ) where {K<:KVector} = basis_1blades(first(k))

"""A KVector formed from all basis 1-vectors set to magnitude 1.0. Useful as a test direction"""
basis_1vector( k::K ) where {K<:KVector} = KVector(one.(basis_1blades(k)))

"""
    coords(k)

coordinate scalar values for ordered k-vector basis of k.
"""
function coords( k::KV ) where {T,K,KV<:KVector{T,K}}
  be = basis_kblades(first(k),K)
  kb = untype.(k)
  [ (i = findfirst(kbᵢ->beᵢ==kbᵢ, kb); isnothing(i) ? zero(T) : scalar(k[i])) for beᵢ in be ]
end

coords( b::Blade ) = coords(KVector(b))

factor( k::K ) where K<:KVector = factor.(k)

outermorphism(L, k::K) where K<:KVector = mapreduce(b->outermorphism(L, b), +, k)

Base.in(be::BK, bs::BK2) where {BK<:Union{Blade,KVector}, BK2<:Union{Blade,KVector}} = iszero(be∧bs)

fieldtype(k::K) where {F<:Number, K<:KVector{F}} = F

function Base.isapprox(k::K, l::K2; kwargs...) where {K<:KVector, K2<:KVector}
  mapreduce( (b,c)->isapprox(b,c; kwargs...), (acc,e)->acc && e, (sort_basis∘prune)(k), (sort_basis∘prune)(l))
end

"""
    reciprocal_frame(fᵢ)

given a set of 1-blades and/or  1-vectors defining a frame of reference,

return a set of reciprocal frame vectors where fᵢ⋅fʲ = δᵢʲ

only valid for metric vectors
"""
function reciprocal_frame(fᵢ::Vector{K}) where {K<:Union{KVector, Blade}}
  Fn⁻¹ = inv(reduce(∧, fᵢ))

  map(1:length(fᵢ)) do i
    ((-1)^(i-1)) * reduce(∧, vcat(fᵢ[1:i-1], fᵢ[i+1:end])) * Fn⁻¹
  end
end


text:
# This file was formerly a part of Julia. License is MIT: https://julialang.org/license

include("../perfutil.jl")

rpath(filename::AbstractString) = joinpath(@__DIR__, filename)

include("binary_trees.jl")
@timeit binary_trees(10) "binary_trees" "Allocate and deallocate many many binary trees"

include("fannkuch.jl")
@timeit fannkuch(7) "fannkuch" "Indexed-access to tiny integer-sequence"

include("fasta.jl")
@timeit fasta(100) "fasta" "Generate and write random DNA sequences"

include("k_nucleotide.jl")
@timeit k_nucleotide(rpath("knucleotide-input.txt")) "k_nucleotide" "Hashtable update and k-nucleotide strings"

include("mandelbrot.jl")
@timeit mandelbrot(200, rpath("mandelbrot-output-julia.txt")) "mandelbrot" "Generate Mandelbrot set portable bitmap file"

include("meteor_contest.jl")
@timeit meteor_contest() "meteor_contest" "Search for solutions to shape packing puzzle"

include("nbody.jl")
using NBody
@timeit NBody.nbody() "nbody" "Double-precision N-body simulation"

include("nbody_vec.jl")
using NBodyVec
@timeit NBodyVec.nbody_vec() "nbody_vec" "A vectorized double-precision N-body simulation"

include("pidigits.jl")
@assert pidigits(1000) == 9216420198
@timeit pidigits(1000) "pidigits" "Streaming arbitrary-precision arithmetic"

include("regex_dna.jl")
@timeit regex_dna(rpath("regexdna-input.txt")) "regex_dna" "Match DNA 8-mers and substitute nucleotides for IUB codes"

include("revcomp.jl")
@timeit revcomp(rpath("revcomp-input.txt")) "revcomp" "Read DNA sequences - write their reverse-complement"

include("spectralnorm.jl")
@timeit spectralnorm() "spectralnorm" "Eigenvalue using the power method"

maxrss("shootout")

ext: jl

max_stars_repo_name: JuliaTagBot/Benchmarks

max_issues_repo_name: JuliaArchive/Benchmarks

max_forks_repo_name: JuliaArchive/Benchmarks


text:
module ParetoSmooth
using Requires
using DocStringExtensions

function __init__()

    @require MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d" begin
        include("MCMCChainsHelpers.jl")
        @require Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0" begin
            include("TuringHelpers.jl")
        end
    end

end

include("AbstractCV.jl")
include("ESS.jl")
include("GPD.jl")
include("InternalHelpers.jl")
include("ImportanceSampling.jl")
include("LeaveOneOut.jl")
include("ModelComparison.jl")
include("NaiveLPD.jl")
include("PublicHelpers.jl")

end


text:
#### Tested moist thermodynamic profiles

#=
This file contains functions to compute all of the
thermodynamic _states_ that Thermodynamics is
tested with in runtests.jl
=#

"""
    unpack_fields(_struct, syms...)

Unpack struct properties `syms`
from struct `_struct`

# Example
```julia
julia> struct Foo;a;b;c;end

julia> f = Foo(1,2,3)
Foo(1, 2, 3)

julia> @unpack_fields f a c; @show a c
a = 1
c = 3
```
"""
macro unpack_fields(_struct, syms...)
    thunk = Expr(:block)
    for sym in syms
        push!(
            thunk.args,
            :($(esc(sym)) = getproperty($(esc(_struct)), $(QuoteNode(sym)))),
        )
    end
    push!(thunk.args, nothing)
    return thunk
end

"""
    ProfileSet

A set of profiles used to test Thermodynamics.
"""
struct ProfileSet{FT}
    z::Array{FT}                        # Altitude
    T::Array{FT}                        # Temperature
    p::Array{FT}                        # Pressure
    RS::Array{FT}                       # Relative humidity
    e_int::Array{FT}                    # Internal energy
    ρ::Array{FT}                        # Density
    θ_liq_ice::Array{FT}                # Potential temperature
    q_tot::Array{FT}                    # Total specific humidity
    q_liq::Array{FT}                    # Liquid specific humidity
    q_ice::Array{FT}                    # Ice specific humidity
    q_pt::Array{PhasePartition{FT}}     # Phase partition
    RH::Array{FT}                       # Relative humidity
    SS::Array{FT}                       # Super saturation
end

"""
    input_config(
        FT;
        n=50,
        n_RS1=10,
        n_RS2=20,
        T_min=FT(150),
        T_surface=FT(350)
    ) where {FT}

Return input arguments to construct profiles
"""
function input_config(
    FT;
    n = 50,
    n_RS1 = 10,
    n_RS2 = 20,
    T_surface = FT(350),
    T_min = FT(150),
)
    n_RS = n_RS1 + n_RS2
    z_range = range(FT(0), stop = FT(2.5e4), length = n)
    relative_sat1 = range(FT(0), stop = FT(1), length = n_RS1)
    relative_sat2 = range(FT(1), stop = FT(1.02), length = n_RS2)
    relative_sat = [relative_sat1..., relative_sat2...]
    return z_range, relative_sat, T_surface, T_min
end

"""
    shared_profiles(
        param_set::AbstractParameterSet,
        z_range::AbstractArray,
        relative_sat::AbstractArray,
        T_surface::FT,
        T_min::FT,
    ) where {FT}

Compute profiles shared across `PhaseDry`,
`PhaseEquil` and `PhaseNonEquil` thermodynamic
states, including:
 - `z` altitude
 - `T_virt` virtual temperature
 - `p` pressure
 - `RS` relative saturation
"""
function shared_profiles(
    param_set::AbstractParameterSet,
    z_range::AbstractArray,
    relative_sat::AbstractArray,
    T_surface::FT,
    T_min::FT,
) where {FT}
    n_RS = length(relative_sat)
    n = length(z_range)
    T_virt = Array{FT}(undef, n * n_RS)
    p = Array{FT}(undef, n * n_RS)
    RS = Array{FT}(undef, n * n_RS)
    z = Array{FT}(undef, n * n_RS)
    linear_indices = LinearIndices((1:n, 1:n_RS))
    profile = DecayingTemperatureProfile{FT}(param_set, T_surface, T_min)
    for i in linear_indices.indices[1]
        for j in linear_indices.indices[2]
            k = linear_indices[i, j]
            z[k] = z_range[i]
            T_virt[k], p[k] = profile(param_set, z[k])
            RS[k] = relative_sat[j]
        end
    end
    return z, T_virt, p, RS
end

####
#### PhaseDry
####

"""
    PhaseDryProfiles(param_set, ::Type{FT})

Returns a `ProfileSet` used to test dry thermodynamic states.
"""
function PhaseDryProfiles(
    param_set::AbstractParameterSet,
    ::Type{FT},
) where {FT}

    z_range, relative_sat, T_surface, T_min = input_config(FT)
    z, T_virt, p, RS =
        shared_profiles(param_set, z_range, relative_sat, T_surface, T_min)
    _R_d::FT = R_d(param_set)
    T = T_virt
    ρ = p ./ (_R_d .* T)

    # Additional variables
    phase_type = PhaseDry
    q_tot = zeros(FT, length(RS))
    q_pt = PhasePartition_equil.(Ref(param_set), T, ρ, q_tot, Ref(phase_type))
    e_int = internal_energy.(Ref(param_set), T, q_pt)
    θ_liq_ice = liquid_ice_pottemp.(Ref(param_set), T, ρ, q_pt)
    q_liq = getproperty.(q_pt, :liq)
    q_ice = getproperty.(q_pt, :ice)
    RH = relative_humidity.(Ref(param_set), T, p, e_int, Ref(phase_type), q_pt)

    # TODO: Update this once a super saturation method exists
    # SS = super_saturation.(Ref(param_set), T, p, e_int, Ref(phase_type), q_pt)
    SS = zeros(FT, length(RH))

    return ProfileSet(
        z,
        T,
        p,
        RS,
        e_int,
        ρ,
        θ_liq_ice,
        q_tot,
        q_liq,
        q_ice,
        q_pt,
        RH,
        SS,
    )
end

####
#### PhaseEquil
####

"""
    PhaseEquilProfiles(param_set, ::Type{FT})

Returns a `ProfileSet` used to test moist states in thermodynamic equilibrium.
"""
function PhaseEquilProfiles(
    param_set::AbstractParameterSet,
    ::Type{FT},
) where {FT}

    z_range, relative_sat, T_surface, T_min = input_config(FT)

    z, T_virt, p, RS =
        shared_profiles(param_set, z_range, relative_sat, T_surface, T_min)
    _R_d::FT = R_d(param_set)
    T = T_virt
    ρ = p ./ (_R_d .* T)

    # Additional variables
    phase_type = PhaseEquil
    q_sat = q_vap_saturation.(Ref(param_set), T, ρ, Ref(phase_type))
    q_tot = min.(RS .* q_sat, FT(1))
    q_pt = PhasePartition_equil.(Ref(param_set), T, ρ, q_tot, Ref(phase_type))
    e_int = internal_energy.(Ref(param_set), T, q_pt)
    θ_liq_ice = liquid_ice_pottemp.(Ref(param_set), T, ρ, q_pt)
    q_liq = getproperty.(q_pt, :liq)
    q_ice = getproperty.(q_pt, :ice)
    RH = relative_humidity.(Ref(param_set), T, p, e_int, Ref(phase_type), q_pt)

    # TODO: Update this once a super saturation method exists
    # SS = super_saturation.(Ref(param_set), T, p, e_int, Ref(phase_type), q_pt)
    SS = zeros(FT, length(RH))

    return ProfileSet(
        z,
        T,
        p,
        RS,
        e_int,
        ρ,
        θ_liq_ice,
        q_tot,
        q_liq,
        q_ice,
        q_pt,
        RH,
        SS,
    )
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
########## Broadcast GBP messages: Factor to variable ##########
function messageFactorVariableBroadcast(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateFactor
        Mrow = 0.0; Vrow = 0.0
        for j in gbp.graph.rowptr[i]
            Mrow += (gbp.graph.coefficient[j] * gbp.inference.meanVariableFactor[j])
            Vrow += (gbp.graph.coefficient[j]^2 * gbp.inference.varianceVariableFactor[j])
        end
        for j in gbp.graph.rowptr[i]
            q = gbp.graph.toVariable.nzval[j]
            gbp.inference.meanFactorVariable[q] = (gbp.graph.meanIndirect[j] - Mrow) / gbp.graph.coefficient[j] + gbp.inference.meanVariableFactor[j]
            gbp.inference.varianceFactorVariable[q] = (gbp.graph.varianceIndirect[j] + Vrow - gbp.graph.coefficient[j]^2 * gbp.inference.varianceVariableFactor[j]) / (gbp.graph.coefficient[j]^2)
        end
    end
end

########## Broadcast GBP means: Factor to variable ##########
function meanFactorVariableBroadcast(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateFactor
        Mrow = 0.0
        for j in gbp.graph.rowptr[i]
            Mrow += (gbp.graph.coefficient[j] * gbp.inference.meanVariableFactor[j])
        end
        for j in gbp.graph.rowptr[i]
            gbp.inference.meanFactorVariable[gbp.graph.toVariable.nzval[j]] = (gbp.graph.meanIndirect[j] - Mrow) / gbp.graph.coefficient[j] + gbp.inference.meanVariableFactor[j]
        end
    end
end

########## Broadcast GBP variances: Factor to variable ##########
function varianceFactorVariableBroadcast(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateFactor
        Vrow = 0.0
        for j in gbp.graph.rowptr[i]
            Vrow += (gbp.graph.coefficient[j]^2 * gbp.inference.varianceVariableFactor[j])
        end
        for j in gbp.graph.rowptr[i]
            gbp.inference.varianceFactorVariable[gbp.graph.toVariable.nzval[j]] = (gbp.graph.varianceIndirect[j] + Vrow - gbp.graph.coefficient[j]^2 * gbp.inference.varianceVariableFactor[j]) / (gbp.graph.coefficient[j]^2)
        end
    end
end

########## Broadcast GBP damp messages: Factor to variable ##########
function messageDampFactorVariableBroadcast(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateFactor
        Mrow = 0.0; Vrow = 0.0
        for j in gbp.graph.rowptr[i]
            Mrow += (gbp.graph.coefficient[j] * gbp.inference.meanVariableFactor[j])
            Vrow += (gbp.graph.coefficient[j]^2 * gbp.inference.varianceVariableFactor[j])
        end
        for j in gbp.graph.rowptr[i]
            q = gbp.graph.toVariable.nzval[j]
            gbp.inference.meanFactorVariable[q] = gbp.graph.alphaNew[j] * ((gbp.graph.meanIndirect[j] - Mrow) / gbp.graph.coefficient[j] + gbp.inference.meanVariableFactor[j]) + gbp.graph.alphaOld[j] * gbp.inference.meanFactorVariable[q]
            gbp.inference.varianceFactorVariable[q] = (gbp.graph.varianceIndirect[j] + Vrow - gbp.graph.coefficient[j]^2 * gbp.inference.varianceVariableFactor[j]) / (gbp.graph.coefficient[j]^2)
        end
    end
end

########## Broadcast GBP damp means: Factor to variable ##########
function meanDampFactorVariableBroadcast(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateFactor
        Mrow = 0.0
        for j in gbp.graph.rowptr[i]
            Mrow += (gbp.graph.coefficient[j] * gbp.inference.meanVariableFactor[j])
        end
        for j in gbp.graph.rowptr[i]
            q = gbp.graph.toVariable.nzval[j]
            gbp.inference.meanFactorVariable[q] = gbp.graph.alphaNew[j] * ((gbp.graph.meanIndirect[j] - Mrow) / gbp.graph.coefficient[j] + gbp.inference.meanVariableFactor[j]) + gbp.graph.alphaOld[j] * gbp.inference.meanFactorVariable[q]
        end
    end
end

########## Broadcast GBP messages: Variable to factor ##########
function messageVariableFactorBroadcast(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateVariable
        Mcol = gbp.graph.meanDirect[i]; Wcol = gbp.graph.weightDirect[i]
        for j in gbp.graph.colptr[i]
            Mcol += gbp.inference.meanFactorVariable[j] / gbp.inference.varianceFactorVariable[j]
            Wcol += 1 / gbp.inference.varianceFactorVariable[j]
        end
        for j in gbp.graph.colptr[i]
            q = gbp.graph.toFactor.nzval[j]
            gbp.inference.varianceVariableFactor[q] = 1 / (Wcol - 1 / gbp.inference.varianceFactorVariable[j])
            gbp.inference.meanVariableFactor[q] = (Mcol - gbp.inference.meanFactorVariable[j] / gbp.inference.varianceFactorVariable[j]) * gbp.inference.varianceVariableFactor[q]
        end
    end
end

########## Broadcast GBP means: Variable to factor ##########
function meanVariableFactorBroadcast(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateVariable
        Mcol = gbp.graph.meanDirect[i]
        for j in gbp.graph.colptr[i]
            Mcol += gbp.inference.meanFactorVariable[j] / gbp.inference.varianceFactorVariable[j]
        end
        for j in gbp.graph.colptr[i]
            q = gbp.graph.toFactor.nzval[j]
            gbp.inference.meanVariableFactor[q] = (Mcol - gbp.inference.meanFactorVariable[j] / gbp.inference.varianceFactorVariable[j]) * gbp.inference.varianceVariableFactor[q]
        end
    end
end

########## Broadcast GBP variances: Variable to factor ##########
function varianceVariableFactorBroadcast(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateVariable
        Wcol = gbp.graph.weightDirect[i]
        for j in gbp.graph.colptr[i]
            Wcol += 1 / gbp.inference.varianceFactorVariable[j]
        end
        for j in gbp.graph.colptr[i]
            gbp.inference.varianceVariableFactor[gbp.graph.toFactor.nzval[j]] = 1 / (Wcol - 1 / gbp.inference.varianceFactorVariable[j])
        end
    end
end


text:
# ---
# title: 65. Valid Number
# id: problem65
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, String
# link: <https://leetcode.com/problems/valid-number/description/>
# hidden: true
# ---
#
# Validate if a given string can be interpreted as a decimal number.
#
# Some examples:
# `"0"` => `true`
# `" 0.1 "` => `true`
# `"abc"` => `false`
# `"1 a"` => `false`
# `"2e10"` => `true`
# `" -90e3   "` => `true`
# `" 1e"` => `false`
# `"e3"` => `false`
# `" 6e-1"` => `true`
# `" 99e2.5 "` => `false`
# `"53.5e93"` => `true`
# `" --6 "` => `false`
# `"-+3"` => `false`
# `"95a54e53"` => `false`
#
# **Note:** It is intended for the problem statement to be ambiguous. It would
# be best if you gathered all requirements up front before implementing one.
# However, here is a list of characters that can be in a valid decimal number:
#
#   * Numbers 0-9
#   * Exponent - "e"
#   * Positive/negative sign - "+"/"-"
#   * Decimal point - "."
#
# Of course, the context of these characters also matters in the input.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "0"
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: s = "3"
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 20`
#   * `s` consists of only English letters, digits, space `' '`, plus `'+'`, minus `'-'`, or dot `'.'`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
## just takes realpart of operator
struct ReOperator{O,T} <: Operator{T}
    op::O
end

ReOperator(op)=ReOperator{typeof(op),Float64}(op)
convert(::Type{Operator{T}},R::ReOperator) where {T} = ReOperator{typeof(R.op),T}(R.op)

@wrapperstructure ReOperator
@wrapperspaces ReOperator

getindex(RI::ReOperator{O,T},k::Integer,j::Integer) where {O,T} =
    convert(T,real(RI.op[k,j]))

choosedomainspace(R::ReOperator,sp::Space) = choosedomainspace(R.op,sp)
for OP in (:promotedomainspace,:promoterangespace)
    @eval begin
        $OP(R::ReOperator,sp::UnsetSpace) = ReOperator($OP(R.op,sp))
        $OP(R::ReOperator,sp::Space) = ReOperator($OP(R.op,sp))
    end
end

# TODO: can't do this because UnsetSpace might change type
#real{T<:Real}(op::Operator{T})=op
real(op::Operator) = ReOperator(op)


text:
using Iris
using Test

@testset "Points" begin
    @testset "1D Cartesian" begin
        p = Point(1)
        @test_throws ArgumentError Point{Cartesian}(Inf)
        @test_throws ArgumentError Point{Cartesian}(NaN)
        @test 3*p == Point{Cartesian}(3)
        @test p*3 == Point{Cartesian}(3)
        @test p+p≈2p
        @test p-p≈0p
        @test ndims(p)==1
        @test p[1] ≈ 1
        @test firstindex(p)==1
        @test p[end] ≈ 1
        @test length(p)==1
        @test size(p)==(1,)
        @test eltype(p)==Float64
        @test p.x ≈ 1
        @test p.y ≈ 0
        @test p.z ≈ 0
        @test p.s ≈ 1
        @test p.r ≈ hypot(1)
        @test p.ϕ ≈ 0
        @test p.θ ≈ π/2
    end
    @testset "2D Cartesian" begin
        p = Point(1,2.)
        @test_throws ArgumentError Point{Cartesian}(Inf,2.)
        @test_throws ArgumentError Point{Cartesian}(1,NaN)
        @test 3*p == Point{Cartesian}(3,6)
        @test p*3 == Point{Cartesian}(3,6)
        @test p+p≈2p
        @test p-p≈0p
        @test ndims(p)==2
        @test p[2]==2
        @test firstindex(p)==1
        @test p[end]==2
        @test length(p)==2
        @test size(p)==(2,)
        @test eltype(p)==Float64
        @test p.x==1
        @test p.y==2
        @test p.r==hypot(1,2)

        p = Point(3,4)
        @test p.x ≈ 3
        @test p.y ≈ 4
        @test p.z ≈ 0
        @test p.s ≈ 5
        @test p.r ≈ 5
        @test p.ϕ ≈ 0.92729521
        @test p.θ ≈ π/2
    end
    @testset "3D Cartesian" begin
        p = Point(1,2,3.)
        @test_throws ArgumentError Point{Cartesian}(Inf,2.,3)
        @test_throws ArgumentError Point{Cartesian}(1,NaN,3)
        @test 3*p == Point{Cartesian}(3,6,9)
        @test p*3 == Point{Cartesian}(3,6,9)
        @test p+p≈2p
        @test p-p≈0p
        @test ndims(p)==3
        @test p[2]==2
        @test firstindex(p)==1
        @test p[end]==3
        @test length(p)==3
        @test size(p)==(3,)
        @test eltype(p)==Float64

        p = Point{Cartesian}(1,sqrt(8),4)
        @test p.x ≈ 1
        @test p.y ≈ sqrt(8)
        @test p.z ≈ 4
        @test p.s ≈ 3
        @test p.r ≈ 5
        @test p.ϕ ≈ 1.2309594173407747
        @test p.θ ≈ 0.6435011087932843
    end
    @testset "4D Cartesian" begin
        p = Point(1,2,3,4.)
        @test_throws ArgumentError Point{Cartesian}(Inf,2.,3,4.)
        @test_throws ArgumentError Point{Cartesian}(1,NaN,3,4.)
        @test 3*p == Point{Cartesian}(3,6,9,12)
        @test p*3 == Point{Cartesian}(3,6,9,12)
        @test p+p≈2p
        @test p-p≈0p
        @test ndims(p)==4
        @test p[2]==2
        @test firstindex(p)==1
        @test p[end]==4
        @test length(p)==4
        @test size(p)==(4,)
        @test eltype(p)==Float64

        p = Point{Cartesian}(1,sqrt(8),4,12)
        @test p.x ≈ 1
        @test p.y ≈ sqrt(8)
        @test p.z ≈ 4
        @test p.s ≈ 3
        @test p.r ≈ 13
        @test_throws ErrorException p.ϕ
        @test_throws ErrorException p.θ
    end
    @testset "1D Polar" begin
        p = Point{Polar}(1)
        @test_throws ArgumentError Point{Polar}(Inf)
        @test_throws ArgumentError Point{Polar}(NaN)
        @test_throws ArgumentError Polar(-1)
        @test 3*p == Point{Polar}(3)
        @test p*3 == Point{Polar}(3)
        @test p+p≈2p
        @test p-p≈0p
        @test ndims(p)==1
        @test firstindex(p)==1
        @test p[end]==1
        @test length(p)==1
        @test size(p)==(1,)
        @test eltype(p)==Float64

        p = Point{Polar}(1)
        @test p.x ≈ 1
        @test p.y ≈ 0
        @test p.z ≈ 0
        @test p.s ≈ 1
        @test p.r ≈ 1
        @test p.ϕ ≈ 0
        @test p.θ ≈ π/2
    end
    @testset "2D Polar" begin
        p = Point{Polar}(1,2.)
        @test_throws ArgumentError Point{Polar}(Inf,2.)
        @test_throws ArgumentError Point{Polar}(NaN,1)
        @test_throws ArgumentError Polar(-1,2.)
        @test 3*p == Point{Polar}(3,2)
        @test p*3 == Point{Polar}(3,2)
        @test p+p≈2p
        @test p-p≈0p
        @test ndims(p)==2
        @test p[2]==2
        @test firstindex(p)==1
        @test p[end]==2
        @test length(p)==2
        @test size(p)==(2,)
        @test eltype(p)==Float64

        p = Point{Polar}(1,π/3)
        @test p.x ≈ 1/2
        @test p.y ≈ sqrt(3)/2
        @test p.z ≈ 0
        @test p.s ≈ 1
        @test p.r ≈ 1
        @test p.ϕ ≈ π/3
        @test p.θ ≈ π/2
    end
    @testset "3D Polar" begin
        p = Point{Polar}(1,2,3.)
        @test_throws ArgumentError Point{Polar}(Inf,2.,NaN)
        @test_throws ArgumentError Point{Polar}(NaN,1,3.)
        @test_throws ArgumentError Point{Polar}(-1,2.,3.)
        @test 3*p == Point{Polar}(3,2,9)
        @test p*3 == Point{Polar}(3,2,9)
        @test p+p≈2p
        @test p-p≈0p
        @test ndims(p)==3
        @test p[2]==2
        @test firstindex(p)==1
        @test p[end]==3
        @test length(p)==3
        @test size(p)==(3,)
        @test eltype(p)==Float64

        p = Point{Polar}(3,π/3,4)
        @test p.x ≈ 3/2
        @test p.y ≈ 3*sqrt(3)/2
        @test p.z ≈ 4
        @test p.s ≈ 3
        @test p.r ≈ 5
        @test p.ϕ ≈ π/3
        @test p.θ ≈ 0.6435011087932843
    end
    @testset "1D Spherical" begin
        p = Point{Spherical}(1)
        @test_throws ArgumentError Point{Spherical}(NaN)
        @test_throws ArgumentError Point{Spherical}(Inf)
        @test_throws ArgumentError Point{Spherical}(-1)
        @test 3*p ≈ Point{Spherical}(3)
        @test p*3 ≈ Point{Spherical}(3)
        @test p+p≈2p
        @test p-p≈0p
        @test ndims(p)==1
        @test firstindex(p)==1
        @test p[end]==1
        @test length(p)==1
        @test size(p)==(1,)
        @test eltype(p)==Float64

        p = Point{Spherical}(3)
        @test p.x ≈ 3
        @test p.y ≈ 0
        @test p.z ≈ 0
        @test p.s ≈ 3
        @test p.r ≈ 3
        @test p.ϕ ≈ 0
        @test p.θ ≈ 0
    end
    @testset "2D Spherical" begin
        p = Point{Spherical}(1,π/3)
        @test_throws ArgumentError Point{Spherical}(NaN,1)
        @test_throws ArgumentError Point{Spherical}(Inf,1)
        @test_throws ArgumentError Point{Spherical}(-1,1)
        @test 3*p ≈ Point{Spherical}(3,π/3)
        @test p*3 ≈ Point{Spherical}(3,π/3)
        @test p+p≈2p
        @test p-p≈0p
        @test ndims(p)==2
        @test firstindex(p)==1
        @test p[end]==π/3
        @test length(p)==2
        @test size(p)==(2,)
        @test eltype(p)==Float64

        p = Point{Spherical}(3,π/3)
        @test p.x ≈ 3*sqrt(3)/2
        @test p.y ≈ 0
        @test p.z ≈ 3/2
        @test p.s ≈ p.x
        @test p.r ≈ 3
        @test p.ϕ ≈ 0
        @test p.θ ≈ π/3
    end
    @testset "3D Spherical" begin
        p = Point{Spherical}(1,2,3)
        @test_throws ArgumentError Point{Spherical}(NaN,1,2)
        @test_throws ArgumentError Point{Spherical}(Inf,1,NaN)
        @test_throws ArgumentError Point{Spherical}(-1,1,1)
        @test 3*p == Point{Spherical}(3,2,3)
        @test p*3 == Point{Spherical}(3,2,3)
        @test p+p≈2p
        @test p-p≈0p
        @test ndims(p)==3
        @test p[2]==2
        @test firstindex(p)==1
        @test p[end]==3
        @test length(p)==3
        @test size(p)==(3,)
        @test eltype(p)==Float64

        p = Point{Spherical}(3,π/3,π/4)
        @test p.x ≈ 3*(sqrt(3)/2)*(1/sqrt(2))
        @test p.y ≈ 3*(sqrt(3)/2)*(1/sqrt(2))
        @test p.z ≈ 3/2
        @test p.s ≈ 3*sqrt(3)/2
        @test p.r ≈ 3
        @test p.ϕ ≈ π/4
        @test p.θ ≈ π/3
    end
end


text:
# Julia wrapper for header: gdal_utils.h
# Automatically generated using Clang.jl

"""
    GDALInfoOptionsNew(char ** papszArgv,
                       GDALInfoOptionsForBinary * psOptionsForBinary) -> GDALInfoOptions *

Allocates a GDALInfoOptions struct.

### Parameters
* **papszArgv**: NULL terminated list of options (potentially including filename and open options too), or NULL. The accepted options are the ones of the gdalinfo utility.
* **psOptionsForBinary**: (output) may be NULL (and should generally be NULL), otherwise (gdalinfo_bin.cpp use case) must be allocated with GDALInfoOptionsForBinaryNew() prior to this function. Will be filled with potentially present filename, open options, subdataset number...

### Returns
pointer to the allocated GDALInfoOptions struct. Must be freed with GDALInfoOptionsFree().
"""
function gdalinfooptionsnew(papszArgv, psOptionsForBinary)
    aftercare(ccall((:GDALInfoOptionsNew, libgdal), Ptr{GDALInfoOptions}, (Ptr{Cstring}, Ptr{GDALInfoOptionsForBinary}), papszArgv, psOptionsForBinary))
end

"""
    GDALInfoOptionsFree(GDALInfoOptions * psOptions) -> void

Frees the GDALInfoOptions struct.

### Parameters
* **psOptions**: the options struct for GDALInfo().
"""
function gdalinfooptionsfree(psOptions)
    aftercare(ccall((:GDALInfoOptionsFree, libgdal), Cvoid, (Ptr{GDALInfoOptions},), psOptions))
end

"""
    GDALInfo(GDALDatasetH hDataset,
             const GDALInfoOptions * psOptions) -> char *

Lists various information about a GDAL supported raster dataset.

### Parameters
* **hDataset**: the dataset handle.
* **psOptions**: the options structure returned by GDALInfoOptionsNew() or NULL.

### Returns
string corresponding to the information about the raster dataset (must be freed with CPLFree()), or NULL in case of error.
"""
function gdalinfo(hDataset, psOptions)
    aftercare(ccall((:GDALInfo, libgdal), Cstring, (GDALDatasetH, Ptr{GDALInfoOptions}), hDataset, psOptions), true)
end

"""
    GDALTranslateOptionsNew(char ** papszArgv,
                            GDALTranslateOptionsForBinary * psOptionsForBinary) -> GDALTranslateOptions *

Allocates a GDALTranslateOptions struct.

### Parameters
* **papszArgv**: NULL terminated list of options (potentially including filename and open options too), or NULL. The accepted options are the ones of the gdal_translate utility.
* **psOptionsForBinary**: (output) may be NULL (and should generally be NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with GDALTranslateOptionsForBinaryNew() prior to this function. Will be filled with potentially present filename, open options,...

### Returns
pointer to the allocated GDALTranslateOptions struct. Must be freed with GDALTranslateOptionsFree().
"""
function gdaltranslateoptionsnew(papszArgv, psOptionsForBinary)
    aftercare(ccall((:GDALTranslateOptionsNew, libgdal), Ptr{GDALTranslateOptions}, (Ptr{Cstring}, Ptr{GDALTranslateOptionsForBinary}), papszArgv, psOptionsForBinary))
end

"""
    GDALTranslateOptionsFree(GDALTranslateOptions * psOptions) -> void

Frees the GDALTranslateOptions struct.

### Parameters
* **psOptions**: the options struct for GDALTranslate().
"""
function gdaltranslateoptionsfree(psOptions)
    aftercare(ccall((:GDALTranslateOptionsFree, libgdal), Cvoid, (Ptr{GDALTranslateOptions},), psOptions))
end

"""
    GDALTranslateOptionsSetProgress(GDALTranslateOptions * psOptions,
                                    GDALProgressFunc pfnProgress,
                                    void * pProgressData) -> void

Set a progress function.

### Parameters
* **psOptions**: the options struct for GDALTranslate().
* **pfnProgress**: the progress callback.
* **pProgressData**: the user data for the progress callback.
"""
function gdaltranslateoptionssetprogress(psOptions, pfnProgress, pProgressData)
    aftercare(ccall((:GDALTranslateOptionsSetProgress, libgdal), Cvoid, (Ptr{GDALTranslateOptions}, GDALProgressFunc, Ptr{Cvoid}), psOptions, pfnProgress, pProgressData))
end

"""
    GDALTranslate(const char * pszDest,
                  GDALDatasetH hSrcDataset,
                  const GDALTranslateOptions * psOptionsIn,
                  int * pbUsageError) -> GDALDatasetH

Converts raster data between different formats.

### Parameters
* **pszDest**: the destination dataset path.
* **hSrcDataset**: the source dataset handle.
* **psOptionsIn**: the options struct returned by GDALTranslateOptionsNew() or NULL.
* **pbUsageError**: the pointer to int variable to determine any usage error has occurred or NULL.

### Returns
the output dataset (new dataset that must be closed using GDALClose()) or NULL in case of error.
"""
function gdaltranslate(pszDestFilename, hSrcDataset, psOptions, pbUsageError)
    aftercare(ccall((:GDALTranslate, libgdal), GDALDatasetH, (Cstring, GDALDatasetH, Ptr{GDALTranslateOptions}, Ptr{Cint}), pszDestFilename, hSrcDataset, psOptions, pbUsageError))
end

"""
    GDALWarpAppOptionsNew(char ** papszArgv,
                          GDALWarpAppOptionsForBinary * psOptionsForBinary) -> GDALWarpAppOptions *

Allocates a GDALWarpAppOptions struct.

### Parameters
* **papszArgv**: NULL terminated list of options (potentially including filename and open options too), or NULL. The accepted options are the ones of the gdalwarp utility.
* **psOptionsForBinary**: (output) may be NULL (and should generally be NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with GDALWarpAppOptionsForBinaryNew() prior to this function. Will be filled with potentially present filename, open options,...

### Returns
pointer to the allocated GDALWarpAppOptions struct. Must be freed with GDALWarpAppOptionsFree().
"""
function gdalwarpappoptionsnew(papszArgv, psOptionsForBinary)
    aftercare(ccall((:GDALWarpAppOptionsNew, libgdal), Ptr{GDALWarpAppOptions}, (Ptr{Cstring}, Ptr{GDALWarpAppOptionsForBinary}), papszArgv, psOptionsForBinary))
end

"""
    GDALWarpAppOptionsFree(GDALWarpAppOptions * psOptions) -> void

Frees the GDALWarpAppOptions struct.

### Parameters
* **psOptions**: the options struct for GDALWarp().
"""
function gdalwarpappoptionsfree(psOptions)
    aftercare(ccall((:GDALWarpAppOptionsFree, libgdal), Cvoid, (Ptr{GDALWarpAppOptions},), psOptions))
end

"""
    GDALWarpAppOptionsSetProgress(GDALWarpAppOptions * psOptions,
                                  GDALProgressFunc pfnProgress,
                                  void * pProgressData) -> void

Set a progress function.

### Parameters
* **psOptions**: the options struct for GDALWarp().
* **pfnProgress**: the progress callback.
* **pProgressData**: the user data for the progress callback.
"""
function gdalwarpappoptionssetprogress(psOptions, pfnProgress, pProgressData)
    aftercare(ccall((:GDALWarpAppOptionsSetProgress, libgdal), Cvoid, (Ptr{GDALWarpAppOptions}, GDALProgressFunc, Ptr{Cvoid}), psOptions, pfnProgress, pProgressData))
end

"""
    GDALWarpAppOptionsSetQuiet(GDALWarpAppOptions * psOptions,
                               int bQuiet) -> void

Set a progress function.

### Parameters
* **psOptions**: the options struct for GDALWarp().
* **bQuiet**: whether GDALWarp() should emit messages on stdout.
"""
function gdalwarpappoptionssetquiet(psOptions, bQuiet)
    aftercare(ccall((:GDALWarpAppOptionsSetQuiet, libgdal), Cvoid, (Ptr{GDALWarpAppOptions}, Cint), psOptions, bQuiet))
end

"""
    GDALWarpAppOptionsSetWarpOption(GDALWarpAppOptions * psOptions,
                                    const char * pszKey,
                                    const char * pszValue) -> void

Set a warp option.

### Parameters
* **psOptions**: the options struct for GDALWarp().
* **pszKey**: key.
* **pszValue**: value.
"""
function gdalwarpappoptionssetwarpoption(psOptions, pszKey, pszValue)
    aftercare(ccall((:GDALWarpAppOptionsSetWarpOption, libgdal), Cvoid, (Ptr{GDALWarpAppOptions}, Cstring, Cstring), psOptions, pszKey, pszValue))
end

"""
    GDALWarp(const char * pszDest,
             GDALDatasetH hDstDS,
             int nSrcCount,
             GDALDatasetH * pahSrcDS,
             const GDALWarpAppOptions * psOptionsIn,
             int * pbUsageError) -> GDALDatasetH

Image reprojection and warping function.

### Parameters
* **pszDest**: the destination dataset path or NULL.
* **hDstDS**: the destination dataset or NULL.
* **nSrcCount**: the number of input datasets.
* **pahSrcDS**: the list of input datasets.
* **psOptionsIn**: the options struct returned by GDALWarpAppOptionsNew() or NULL.
* **pbUsageError**: the pointer to int variable to determine any usage error has occurred.

### Returns
the output dataset (new dataset that must be closed using GDALClose(), or hDstDS if not NULL) or NULL in case of error.
"""
function gdalwarp(pszDest, hDstDS, nSrcCount, pahSrcDS, psOptions, pbUsageError)
    aftercare(ccall((:GDALWarp, libgdal), GDALDatasetH, (Cstring, GDALDatasetH, Cint, Ptr{GDALDatasetH}, Ptr{GDALWarpAppOptions}, Ptr{Cint}), pszDest, hDstDS, nSrcCount, pahSrcDS, psOptions, pbUsageError))
end

"""
    GDALVectorTranslateOptionsNew(char ** papszArgv,
                                  GDALVectorTranslateOptionsForBinary * psOptionsForBinary) -> GDALVectorTranslateOptions *

allocates a GDALVectorTranslateOptions struct.

### Parameters
* **papszArgv**: NULL terminated list of options (potentially including filename and open options too), or NULL. The accepted options are the ones of the ogr2ogr utility.
* **psOptionsForBinary**: (output) may be NULL (and should generally be NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with GDALVectorTranslateOptionsForBinaryNew() prior to this function. Will be filled with potentially present filename, open options,...

### Returns
pointer to the allocated GDALVectorTranslateOptions struct. Must be freed with GDALVectorTranslateOptionsFree().
"""
function gdalvectortranslateoptionsnew(papszArgv, psOptionsForBinary)
    aftercare(ccall((:GDALVectorTranslateOptionsNew, libgdal), Ptr{GDALVectorTranslateOptions}, (Ptr{Cstring}, Ptr{GDALVectorTranslateOptionsForBinary}), papszArgv, psOptionsForBinary))
end

"""
    GDALVectorTranslateOptionsFree(GDALVectorTranslateOptions * psOptions) -> void

Frees the GDALVectorTranslateOptions struct.

### Parameters
* **psOptions**: the options struct for GDALVectorTranslate().
"""
function gdalvectortranslateoptionsfree(psOptions)
    aftercare(ccall((:GDALVectorTranslateOptionsFree, libgdal), Cvoid, (Ptr{GDALVectorTranslateOptions},), psOptions))
end

"""
    GDALVectorTranslateOptionsSetProgress(GDALVectorTranslateOptions * psOptions,
                                          GDALProgressFunc pfnProgress,
                                          void * pProgressData) -> void

Set a progress function.

### Parameters
* **psOptions**: the options struct for GDALVectorTranslate().
* **pfnProgress**: the progress callback.
* **pProgressData**: the user data for the progress callback.
"""
function gdalvectortranslateoptionssetprogress(psOptions, pfnProgress, pProgressData)
    aftercare(ccall((:GDALVectorTranslateOptionsSetProgress, libgdal), Cvoid, (Ptr{GDALVectorTranslateOptions}, GDALProgressFunc, Ptr{Cvoid}), psOptions, pfnProgress, pProgressData))
end

"""
    GDALVectorTranslate(const char * pszDest,
                        GDALDatasetH hDstDS,
                        int nSrcCount,
                        GDALDatasetH * pahSrcDS,
                        const GDALVectorTranslateOptions * psOptionsIn,
                        int * pbUsageError) -> GDALDatasetH

Converts vector data between file formats.

### Parameters
* **pszDest**: the destination dataset path or NULL.
* **hDstDS**: the destination dataset or NULL.
* **nSrcCount**: the number of input datasets (only 1 supported currently)
* **pahSrcDS**: the list of input datasets.
* **psOptionsIn**: the options struct returned by GDALVectorTranslateOptionsNew() or NULL.
* **pbUsageError**: the pointer to int variable to determine any usage error has occurred

### Returns
the output dataset (new dataset that must be closed using GDALClose(), or hDstDS is not NULL) or NULL in case of error.
"""
function gdalvectortranslate(pszDest, hDstDS, nSrcCount, pahSrcDS, psOptions, pbUsageError)
    aftercare(ccall((:GDALVectorTranslate, libgdal), GDALDatasetH, (Cstring, GDALDatasetH, Cint, Ptr{GDALDatasetH}, Ptr{GDALVectorTranslateOptions}, Ptr{Cint}), pszDest, hDstDS, nSrcCount, pahSrcDS, psOptions, pbUsageError))
end

"""
    GDALDEMProcessingOptionsNew(char ** papszArgv,
                                GDALDEMProcessingOptionsForBinary * psOptionsForBinary) -> GDALDEMProcessingOptions *

Allocates a GDALDEMProcessingOptions struct.

### Parameters
* **papszArgv**: NULL terminated list of options (potentially including filename and open options too), or NULL. The accepted options are the ones of the gdaldem utility.
* **psOptionsForBinary**: (output) may be NULL (and should generally be NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with GDALDEMProcessingOptionsForBinaryNew() prior to this function. Will be filled with potentially present filename, open options,...

### Returns
pointer to the allocated GDALDEMProcessingOptions struct. Must be freed with GDALDEMProcessingOptionsFree().
"""
function gdaldemprocessingoptionsnew(papszArgv, psOptionsForBinary)
    aftercare(ccall((:GDALDEMProcessingOptionsNew, libgdal), Ptr{GDALDEMProcessingOptions}, (Ptr{Cstring}, Ptr{GDALDEMProcessingOptionsForBinary}), papszArgv, psOptionsForBinary))
end

"""
    GDALDEMProcessingOptionsFree(GDALDEMProcessingOptions * psOptions) -> void

Frees the GDALDEMProcessingOptions struct.

### Parameters
* **psOptions**: the options struct for GDALDEMProcessing().
"""
function gdaldemprocessingoptionsfree(psOptions)
    aftercare(ccall((:GDALDEMProcessingOptionsFree, libgdal), Cvoid, (Ptr{GDALDEMProcessingOptions},), psOptions))
end

"""
    GDALDEMProcessingOptionsSetProgress(GDALDEMProcessingOptions * psOptions,
                                        GDALProgressFunc pfnProgress,
                                        void * pProgressData) -> void

Set a progress function.

### Parameters
* **psOptions**: the options struct for GDALDEMProcessing().
* **pfnProgress**: the progress callback.
* **pProgressData**: the user data for the progress callback.
"""
function gdaldemprocessingoptionssetprogress(psOptions, pfnProgress, pProgressData)
    aftercare(ccall((:GDALDEMProcessingOptionsSetProgress, libgdal), Cvoid, (Ptr{GDALDEMProcessingOptions}, GDALProgressFunc, Ptr{Cvoid}), psOptions, pfnProgress, pProgressData))
end

"""
    GDALDEMProcessing(const char * pszDest,
                      GDALDatasetH hSrcDataset,
                      const char * pszProcessing,
                      const char * pszColorFilename,
                      const GDALDEMProcessingOptions * psOptionsIn,
                      int * pbUsageError) -> GDALDatasetH

Apply a DEM processing.

### Parameters
* **pszDest**: the destination dataset path.
* **hSrcDataset**: the source dataset handle.
* **pszProcessing**: the processing to apply (one of "hillshade", "slope", "aspect", "color-relief", "TRI", "TPI", "Roughness")
* **pszColorFilename**: color file (mandatory for "color-relief" processing, should be NULL otherwise)
* **psOptionsIn**: the options struct returned by GDALDEMProcessingOptionsNew() or NULL.
* **pbUsageError**: the pointer to int variable to determine any usage error has occurred or NULL.

### Returns
the output dataset (new dataset that must be closed using GDALClose()) or NULL in case of error.
"""
function gdaldemprocessing(pszDestFilename, hSrcDataset, pszProcessing, pszColorFilename, psOptions, pbUsageError)
    aftercare(ccall((:GDALDEMProcessing, libgdal), GDALDatasetH, (Cstring, GDALDatasetH, Cstring, Cstring, Ptr{GDALDEMProcessingOptions}, Ptr{Cint}), pszDestFilename, hSrcDataset, pszProcessing, pszColorFilename, psOptions, pbUsageError))
end

"""
    GDALNearblackOptionsNew(char ** papszArgv,
                            GDALNearblackOptionsForBinary * psOptionsForBinary) -> GDALNearblackOptions *

Allocates a GDALNearblackOptions struct.

### Parameters
* **papszArgv**: NULL terminated list of options (potentially including filename and open options too), or NULL. The accepted options are the ones of the nearblack utility.
* **psOptionsForBinary**: (output) may be NULL (and should generally be NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with GDALNearblackOptionsForBinaryNew() prior to this function. Will be filled with potentially present filename, open options,...

### Returns
pointer to the allocated GDALNearblackOptions struct. Must be freed with GDALNearblackOptionsFree().
"""
function gdalnearblackoptionsnew(papszArgv, psOptionsForBinary)
    aftercare(ccall((:GDALNearblackOptionsNew, libgdal), Ptr{GDALNearblackOptions}, (Ptr{Cstring}, Ptr{GDALNearblackOptionsForBinary}), papszArgv, psOptionsForBinary))
end

"""
    GDALNearblackOptionsFree(GDALNearblackOptions * psOptions) -> void

Frees the GDALNearblackOptions struct.

### Parameters
* **psOptions**: the options struct for GDALNearblack().
"""
function gdalnearblackoptionsfree(psOptions)
    aftercare(ccall((:GDALNearblackOptionsFree, libgdal), Cvoid, (Ptr{GDALNearblackOptions},), psOptions))
end

"""
    GDALNearblackOptionsSetProgress(GDALNearblackOptions * psOptions,
                                    GDALProgressFunc pfnProgress,
                                    void * pProgressData) -> void

Set a progress function.

### Parameters
* **psOptions**: the options struct for GDALNearblack().
* **pfnProgress**: the progress callback.
* **pProgressData**: the user data for the progress callback.
"""
function gdalnearblackoptionssetprogress(psOptions, pfnProgress, pProgressData)
    aftercare(ccall((:GDALNearblackOptionsSetProgress, libgdal), Cvoid, (Ptr{GDALNearblackOptions}, GDALProgressFunc, Ptr{Cvoid}), psOptions, pfnProgress, pProgressData))
end

"""
    GDALNearblack(const char * pszDest,
                  GDALDatasetH hDstDS,
                  GDALDatasetH hSrcDataset,
                  const GDALNearblackOptions * psOptionsIn,
                  int * pbUsageError) -> GDALDatasetH

Convert nearly black/white borders to exact value.

### Parameters
* **pszDest**: the destination dataset path or NULL.
* **hDstDS**: the destination dataset or NULL. Might be equal to hSrcDataset.
* **hSrcDataset**: the source dataset handle.
* **psOptionsIn**: the options struct returned by GDALNearblackOptionsNew() or NULL.
* **pbUsageError**: the pointer to int variable to determine any usage error has occurred or NULL.

### Returns
the output dataset (new dataset that must be closed using GDALClose(), or hDstDS is not NULL) or NULL in case of error.
"""
function gdalnearblack(pszDest, hDstDS, hSrcDS, psOptions, pbUsageError)
    aftercare(ccall((:GDALNearblack, libgdal), GDALDatasetH, (Cstring, GDALDatasetH, GDALDatasetH, Ptr{GDALNearblackOptions}, Ptr{Cint}), pszDest, hDstDS, hSrcDS, psOptions, pbUsageError))
end

"""
    GDALGridOptionsNew(char ** papszArgv,
                       GDALGridOptionsForBinary * psOptionsForBinary) -> GDALGridOptions *

Allocates a GDALGridOptions struct.

### Parameters
* **papszArgv**: NULL terminated list of options (potentially including filename and open options too), or NULL. The accepted options are the ones of the gdal_translate utility.
* **psOptionsForBinary**: (output) may be NULL (and should generally be NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with GDALGridOptionsForBinaryNew() prior to this function. Will be filled with potentially present filename, open options,...

### Returns
pointer to the allocated GDALGridOptions struct. Must be freed with GDALGridOptionsFree().
"""
function gdalgridoptionsnew(papszArgv, psOptionsForBinary)
    aftercare(ccall((:GDALGridOptionsNew, libgdal), Ptr{GDALGridOptions}, (Ptr{Cstring}, Ptr{GDALGridOptionsForBinary}), papszArgv, psOptionsForBinary))
end

"""
    GDALGridOptionsFree(GDALGridOptions * psOptions) -> void

Frees the GDALGridOptions struct.

### Parameters
* **psOptions**: the options struct for GDALGrid().
"""
function gdalgridoptionsfree(psOptions)
    aftercare(ccall((:GDALGridOptionsFree, libgdal), Cvoid, (Ptr{GDALGridOptions},), psOptions))
end

"""
    GDALGridOptionsSetProgress(GDALGridOptions * psOptions,
                               GDALProgressFunc pfnProgress,
                               void * pProgressData) -> void

Set a progress function.

### Parameters
* **psOptions**: the options struct for GDALGrid().
* **pfnProgress**: the progress callback.
* **pProgressData**: the user data for the progress callback.
"""
function gdalgridoptionssetprogress(psOptions, pfnProgress, pProgressData)
    aftercare(ccall((:GDALGridOptionsSetProgress, libgdal), Cvoid, (Ptr{GDALGridOptions}, GDALProgressFunc, Ptr{Cvoid}), psOptions, pfnProgress, pProgressData))
end

"""
    GDALGrid(const char * pszDest,
             GDALDatasetH hSrcDataset,
             const GDALGridOptions * psOptionsIn,
             int * pbUsageError) -> GDALDatasetH

Create raster from the scattered data.

### Parameters
* **pszDest**: the destination dataset path.
* **hSrcDataset**: the source dataset handle.
* **psOptionsIn**: the options struct returned by GDALGridOptionsNew() or NULL.
* **pbUsageError**: the pointer to int variable to determine any usage error has occurred or NULL.

### Returns
the output dataset (new dataset that must be closed using GDALClose()) or NULL in case of error.
"""
function gdalgrid(pszDest, hSrcDS, psOptions, pbUsageError)
    aftercare(ccall((:GDALGrid, libgdal), GDALDatasetH, (Cstring, GDALDatasetH, Ptr{GDALGridOptions}, Ptr{Cint}), pszDest, hSrcDS, psOptions, pbUsageError))
end

"""
    GDALRasterizeOptionsNew(char ** papszArgv,
                            GDALRasterizeOptionsForBinary * psOptionsForBinary) -> GDALRasterizeOptions *

Allocates a GDALRasterizeOptions struct.

### Parameters
* **papszArgv**: NULL terminated list of options (potentially including filename and open options too), or NULL. The accepted options are the ones of the gdal_rasterize utility.
* **psOptionsForBinary**: (output) may be NULL (and should generally be NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with GDALRasterizeOptionsForBinaryNew() prior to this function. Will be filled with potentially present filename, open options,...

### Returns
pointer to the allocated GDALRasterizeOptions struct. Must be freed with GDALRasterizeOptionsFree().
"""
function gdalrasterizeoptionsnew(papszArgv, psOptionsForBinary)
    aftercare(ccall((:GDALRasterizeOptionsNew, libgdal), Ptr{GDALRasterizeOptions}, (Ptr{Cstring}, Ptr{GDALRasterizeOptionsForBinary}), papszArgv, psOptionsForBinary))
end

"""
    GDALRasterizeOptionsFree(GDALRasterizeOptions * psOptions) -> void

Frees the GDALRasterizeOptions struct.

### Parameters
* **psOptions**: the options struct for GDALRasterize().
"""
function gdalrasterizeoptionsfree(psOptions)
    aftercare(ccall((:GDALRasterizeOptionsFree, libgdal), Cvoid, (Ptr{GDALRasterizeOptions},), psOptions))
end

"""
    GDALRasterizeOptionsSetProgress(GDALRasterizeOptions * psOptions,
                                    GDALProgressFunc pfnProgress,
                                    void * pProgressData) -> void

Set a progress function.

### Parameters
* **psOptions**: the options struct for GDALRasterize().
* **pfnProgress**: the progress callback.
* **pProgressData**: the user data for the progress callback.
"""
function gdalrasterizeoptionssetprogress(psOptions, pfnProgress, pProgressData)
    aftercare(ccall((:GDALRasterizeOptionsSetProgress, libgdal), Cvoid, (Ptr{GDALRasterizeOptions}, GDALProgressFunc, Ptr{Cvoid}), psOptions, pfnProgress, pProgressData))
end

"""
    GDALRasterize(const char * pszDest,
                  GDALDatasetH hDstDS,
                  GDALDatasetH hSrcDataset,
                  const GDALRasterizeOptions * psOptionsIn,
                  int * pbUsageError) -> GDALDatasetH

Burns vector geometries into a raster.

### Parameters
* **pszDest**: the destination dataset path or NULL.
* **hDstDS**: the destination dataset or NULL.
* **hSrcDataset**: the source dataset handle.
* **psOptionsIn**: the options struct returned by GDALRasterizeOptionsNew() or NULL.
* **pbUsageError**: the pointer to int variable to determine any usage error has occurred or NULL.

### Returns
the output dataset (new dataset that must be closed using GDALClose(), or hDstDS is not NULL) or NULL in case of error.
"""
function gdalrasterize(pszDest, hDstDS, hSrcDS, psOptions, pbUsageError)
    aftercare(ccall((:GDALRasterize, libgdal), GDALDatasetH, (Cstring, GDALDatasetH, GDALDatasetH, Ptr{GDALRasterizeOptions}, Ptr{Cint}), pszDest, hDstDS, hSrcDS, psOptions, pbUsageError))
end

"""
    GDALBuildVRTOptionsNew(char ** papszArgv,
                           GDALBuildVRTOptionsForBinary * psOptionsForBinary) -> GDALBuildVRTOptions *

Allocates a GDALBuildVRTOptions struct.

### Parameters
* **papszArgv**: NULL terminated list of options (potentially including filename and open options too), or NULL. The accepted options are the ones of the gdalbuildvrt utility.
* **psOptionsForBinary**: (output) may be NULL (and should generally be NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with GDALBuildVRTOptionsForBinaryNew() prior to this function. Will be filled with potentially present filename, open options,...

### Returns
pointer to the allocated GDALBuildVRTOptions struct. Must be freed with GDALBuildVRTOptionsFree().
"""
function gdalbuildvrtoptionsnew(papszArgv, psOptionsForBinary)
    aftercare(ccall((:GDALBuildVRTOptionsNew, libgdal), Ptr{GDALBuildVRTOptions}, (Ptr{Cstring}, Ptr{GDALBuildVRTOptionsForBinary}), papszArgv, psOptionsForBinary))
end

"""
    GDALBuildVRTOptionsFree(GDALBuildVRTOptions * psOptions) -> void

Frees the GDALBuildVRTOptions struct.

### Parameters
* **psOptions**: the options struct for GDALBuildVRT().
"""
function gdalbuildvrtoptionsfree(psOptions)
    aftercare(ccall((:GDALBuildVRTOptionsFree, libgdal), Cvoid, (Ptr{GDALBuildVRTOptions},), psOptions))
end

"""
    GDALBuildVRTOptionsSetProgress(GDALBuildVRTOptions * psOptions,
                                   GDALProgressFunc pfnProgress,
                                   void * pProgressData) -> void

Set a progress function.

### Parameters
* **psOptions**: the options struct for GDALBuildVRT().
* **pfnProgress**: the progress callback.
* **pProgressData**: the user data for the progress callback.
"""
function gdalbuildvrtoptionssetprogress(psOptions, pfnProgress, pProgressData)
    aftercare(ccall((:GDALBuildVRTOptionsSetProgress, libgdal), Cvoid, (Ptr{GDALBuildVRTOptions}, GDALProgressFunc, Ptr{Cvoid}), psOptions, pfnProgress, pProgressData))
end

"""
    GDALBuildVRT(const char * pszDest,
                 int nSrcCount,
                 GDALDatasetH * pahSrcDS,
                 const char *const * papszSrcDSNames,
                 const GDALBuildVRTOptions * psOptionsIn,
                 int * pbUsageError) -> GDALDatasetH

Build a VRT from a list of datasets.

### Parameters
* **pszDest**: the destination dataset path.
* **nSrcCount**: the number of input datasets.
* **pahSrcDS**: the list of input datasets (or NULL, exclusive with papszSrcDSNames)
* **papszSrcDSNames**: the list of input dataset names (or NULL, exclusive with pahSrcDS)
* **psOptionsIn**: the options struct returned by GDALBuildVRTOptionsNew() or NULL.
* **pbUsageError**: the pointer to int variable to determine any usage error has occurred.

### Returns
the output dataset (new dataset that must be closed using GDALClose()) or NULL in case of error.
"""
function gdalbuildvrt(pszDest, nSrcCount, pahSrcDS, papszSrcDSNames, psOptions, pbUsageError)
    aftercare(ccall((:GDALBuildVRT, libgdal), GDALDatasetH, (Cstring, Cint, Ptr{GDALDatasetH}, Ptr{Cstring}, Ptr{GDALBuildVRTOptions}, Ptr{Cint}), pszDest, nSrcCount, pahSrcDS, papszSrcDSNames, psOptions, pbUsageError))
end


text:
"""
**Heist.jl** is a Julia package that wraps art's gallery.
The name comes from James Stapleton's similar package for Python.

Note that you must run the Julia executable that is compatable with Root.

You must have the gallery package and dependencies set up in your environment.
You must have your data product packages and dependencies set up in your environment.
"""
module Heist

using Cxx, ROOT  # ROOT is needed for TFile

# Load header paths
addHeaderDir(ENV["GALLERY_INC"], kind=C_System)
addHeaderDir(ENV["CANVAS_INC"], kind=C_System)
addHeaderDir(ENV["CLHEP_INC"], kind=C_System)
addHeaderDir(ENV["BOOST_INC"], kind=C_System)
addHeaderDir(ENV["CETLIB_INC"], kind=C_System)
addHeaderDir(ENV["CETLIB_EXCEPT_INC"], kind=C_System)
addHeaderDir(ENV["FHICLCPP_INC"], kind=C_System)
addHeaderDir(ENV["MESSAGEFACILITY_INC"], kind=C_System)

addHeaderDir(ENV["G4INCLUDE"], kind=C_System)

addHeaderDir(ENV["GM2RINGSIM_INC"], kind=C_System)
addHeaderDir(ENV["GM2DATAPRODUCTS_INC"], kind=C_System)

# Load gallery and dependent libraries
Libdl.dlopen("libgallery", Libdl.RTLD_GLOBAL)
Libdl.dlopen("libcanvas", Libdl.RTLD_GLOBAL)
Libdl.dlopen("libCLHEP", Libdl.RTLD_GLOBAL)

# Load some headers that we'll need
cxx"""
   #include "canvas/Utilities/InputTag.h"
   #include "gallery/Event.h"

   #include <vector>
"""

# files.jl
export
    HeistFiles

# recordSpec.jl
export
    HeistRecordSpec

# event.jl
export
    HeistEvent,
    nextEvent,
    atEnd

# events.jl
export
    HeistEvents

# getRecord.jl
export
    getRecord

include("files.jl")
include("recordSpec.jl")
include("event.jl")
include("events.jl")
include("getRecord.jl")

end # module


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the API function to invert rotations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/inv_rotations.jl
# ============================

# Functions: inv_rotation
# -----------------------

@testset "Invert rotations (Float64)" begin
    T = Float64

    # DCM
    # ==========================================================================

    # Create a random DCM.
    D = create_rotation_matrix(_rand_ang(T), :Z) *
        create_rotation_matrix(_rand_ang(T), :Y) *
        create_rotation_matrix(_rand_ang(T), :X)

    Di = inv_rotation(D)
    @test eltype(Di) === T

    Die = inv(D)
    @test Di ≈ Die

    # Quaternion
    # ==========================================================================

    # Create a random quaternion.
    q = Quaternion(_rand_ang(T), _rand_ang(T), _rand_ang(T), _rand_ang(T))
    q = q / norm(q)

    qi = inv_rotation(q)
    @test eltype(qi) === T

    qie = inv(q)
    @test qi ≈ qie

    # Euler angle and axis
    # ==========================================================================

    # Create a random Euler angle and axis.
    v = @SVector rand(T, 3)
    a = _rand_ang(T)
    av = EulerAngleAxis(a, v)

    avi = inv_rotation(av)
    @test eltype(avi) === T

    avie = inv(av)
    @test avi ≈ avie

    # Euler angles
    # ==========================================================================

    # Create random Euler angles.
    ea = EulerAngles(_rand_ang(T), _rand_ang(T), _rand_ang(T), rand(valid_rot_seqs))

    eai = inv_rotation(ea)
    @test eltype(eai) === T

    eaie = inv(ea)
    @test eai ≈ eaie
end

@testset "Invert rotations (Float32)" begin
    T = Float32

    # DCM
    # ==========================================================================

    # Create a random DCM.
    D = create_rotation_matrix(_rand_ang(T), :Z) *
        create_rotation_matrix(_rand_ang(T), :Y) *
        create_rotation_matrix(_rand_ang(T), :X)

    Di = inv_rotation(D)
    @test eltype(Di) === T

    Die = inv(D)
    @test Di ≈ Die

    # Quaternion
    # ==========================================================================

    # Create a random quaternion.
    q = Quaternion(_rand_ang(T), _rand_ang(T), _rand_ang(T), _rand_ang(T))
    q = q / norm(q)

    qi = inv_rotation(q)
    @test eltype(qi) === T

    qie = inv(q)
    @test qi ≈ qie

    # Euler angle and axis
    # ==========================================================================

    # Create a random Euler angle and axis.
    v = @SVector rand(T, 3)
    a = _rand_ang(T)
    av = EulerAngleAxis(a, v)

    avi = inv_rotation(av)
    @test eltype(avi) === T

    avie = inv(av)
    @test avi ≈ avie

    # Euler angles
    # ==========================================================================

    # Create random Euler angles.
    ea = EulerAngles(_rand_ang(T), _rand_ang(T), _rand_ang(T), rand(valid_rot_seqs))

    eai = inv_rotation(ea)
    @test eltype(eai) === T

    eaie = inv(ea)
    @test eai ≈ eaie
end


text:
using Test

@testset "test folder/one" begin
    @test true
    @test true
    @test false
end

# @testset "test folder/err" begin
#     error("this is an error")
# end


text:
using LinearAlgebra: norm, dot

function gdoptimize(f, g!, fg!, x0::AbstractArray{T}, nmf_flag,
                    maxiter::Int = 1,
                    g_rtol::T = sqrt(eps(T)), g_atol::T = eps(T)) where T <: Number

	linesearch = BackTracking(order=3)
#	linesearch = StrongWolfe()

	x = copy(x0)
	xx = copy(x0)
	gvec = similar(x)
	g!(gvec, x)
	fx = f(x)

	gnorm = norm(gvec)
	gtol = max(g_rtol*gnorm, g_atol)

	# Univariate line search functions
	function ϕ(α)
		for i in eachindex(xx)
			xx[i] = x[i] + α * s[i]
		end
		nmf_flag && prox!(xx)
		f(xx)
	end
	function dϕ(α)
		g!(gvec, x)
		return dot(gvec, s)
	end
	function ϕdϕ(α)
		for i in eachindex(xx)
			xx[i] = x[i] + α * s[i]
		end
		nmf_flag && prox!(xx)
		phi = fg!(nothing, xx)
		fg!(gvec, x)
		dphi = dot(gvec, s)
		return (phi, dphi)
	end

	s = similar(gvec) # Step direction

	iter = 0
	while iter < maxiter && gnorm > gtol
		iter += 1
		s .= -gvec

		dϕ_0 = dot(s, gvec)
		α, fx = linesearch(ϕ, dϕ, ϕdϕ, 1.0, fx, dϕ_0)

		for i in eachindex(x)
			x[i] = x[i] + α * s[i]
		end
		nmf_flag && prox!(x)
		g!(gvec, x)
		gnorm = norm(gvec)
	end

	copyto!(x0, x)

	return fx
end


text:
module ResTime
  using PDBTools
  using FortranFiles
  using Printf
  using StructTypes
  using StaticArrays
  using CellListMap
  using LinearAlgebra
  using LsqFit
  using Polynomials
  using EasyFit

  export distance
  export correlation
  export corr_bench
  export autocorr
  export autocorr_cell
  export fitting

  include("Selection.jl")
  include("FileOperations.jl")
  include("viewmol.jl")
  include("distance.jl")

  include("Trajectory.jl")
  include("ChemFiles.jl")
  include("NamdDCD.jl")
  include("PDBTraj.jl")

  include("sample_space.jl")
  include("survivor_prob.jl")
  include("print_file.jl")
  include("events.jl")
  include("domain_eval.jl")
  include("correlation.jl")
  include("autocorr.jl")
  include("autocorr_CellList.jl")
  # include("fitting.jl")

end


text:
using CSV, StatsBase, DataFrames

include("common.jl")

APPROACHES_TO_COMPARE = [
    "grb-$method-$num_layers" for
    method in ("nn", "cda", "standard-nmdt", "tightened-nmdt"),
    num_layers in (2, 4, 6, 8)
]

if length(ARGS) != 1 || !ispath(ARGS[1])
    error("Invalid csv path.")
end
const CSV_PATH = ARGS[1]

_shifted_geomean(vals, shift) = StatsBase.geomean(vals .+ shift) - shift

df = CSV.read(CSV_PATH, DataFrame)

# Figure out the best primal objective value available in dataset
# for each instance.
df[!, :best_primal] .= 0.0
for i = 1:size(df, 1)
    row = df[i, :]
    to_slice = df[df.instance.==row.instance, :true_objective]
    best_primal_value = minimum(filter(!isnan, to_slice))
    cutoffs = unique(df[df.instance.==row.instance, :cutoff])
    @assert length(cutoffs) == 1
    cutoff = cutoffs[1]
    @assert !isnan(best_primal_value)
    best_primal_value = min(best_primal_value, cutoff)
    df.best_primal[i] = best_primal_value
end

# Add a column to record MIP gap, with respect to best primal solution in data set.
function _mip_gap(row)
    return 100 * abs(row.dual_bound - row.best_primal) / abs(row.best_primal)
end
df[!, :gap] .= 0.0
for i = 1:size(df, 1)
    df.gap[i] = _mip_gap(df[i, :])
end

# Split off problem classes based on how easy they were. We will do separate analysis
# on each problem class.
PROBLEM_CLASS_DFS = Dict()
PROBLEM_CLASS_DFS["solved"] = filter(df) do row
    instance = row.instance
    slice = df[df.instance.==row.instance, :]
    return all(slice.termination_status .== "OPTIMAL")
end

PROBLEM_CLASS_DFS["contested"] = filter(df) do row
    instance = row.instance
    slice = df[df.instance.==row.instance, :]
    return any(slice.termination_status .== "OPTIMAL") &&
           !all(slice.termination_status .== "OPTIMAL")
end

PROBLEM_CLASS_DFS["unsolved"] = filter(df) do row
    instance = row.instance
    slice = df[df.instance.==row.instance, :]
    return !any(slice.termination_status .== "OPTIMAL")
end

# 1. Solved instances
for (class, subset_df) in PROBLEM_CLASS_DFS
    mask = map(t -> t in APPROACHES_TO_COMPARE, subset_df.approach)
    subset_df = subset_df[mask, :]
    total_instance_count = length(unique(subset_df.instance))
    println("Class: $class ($total_instance_count instances)")
    println("="^(30))

    winners = Dict(approach => 0 for approach in APPROACHES_TO_COMPARE)
    fails = Dict(approach => 0 for approach in APPROACHES_TO_COMPARE)
    best_bound = Dict(approach => 0 for approach in APPROACHES_TO_COMPARE)
    for instance in unique(subset_df.instance)
        slice = subset_df[subset_df.instance.==instance, :]
        if Set(unique(slice.approach)) ⊉ Set(APPROACHES_TO_COMPARE)
            @show Set(unique(slice.approach))
            @show Set(APPROACHES_TO_COMPARE)
            error("Incomplete experiment in table.")
        end
        if size(slice, 1) == 0
            error("Unexpected case, aborting.")
        end
        best_solve_time = minimum(slice.solve_time_sec)
        best_gap = minimum(slice.gap)
        all_time_out = best_solve_time >= 0.99 * TIME_LIMIT
        # NOTE: This excludes ties from the tally!
        if size(slice[slice.solve_time_sec.==best_solve_time, :], 1) != 1
            @show slice
            @warn("Need to handle ties")
        end
        for approach in APPROACHES_TO_COMPARE
            sub_slice = slice[slice.approach.==approach, :]
            solve_times = sub_slice.solve_time_sec
            @assert size(solve_times, 1) == 1
            gaps = sub_slice.gap
            solve_time = minimum(solve_times)
            gap = minimum(gaps)
            best_primal = minimum(sub_slice.best_primal)
            solver_time_out = solve_time >= 0.99 * TIME_LIMIT
            if solver_time_out
                fails[approach] += 1
            end
            if all_time_out
                if gap == best_gap
                    winners[approach] += 1
                end
            else
                if solve_time == best_solve_time
                    winners[approach] += 1
                    # break
                end
            end
            # Either attain best bound, or reach Gurobi's optimality cutoff,
            # which is 1e-4 * |primal_bound|
            # Note that the gaps are recorded in the table as percentages,
            # so we must multiply through by 100.
            if gap == best_gap || gap < 1e-4 * 100
                best_bound[approach] += 1
            end
        end
    end

    time_shift = minimum(subset_df.solve_time_sec)
    gap_shift = max(1e-4, minimum(subset_df.gap))

    for approach in APPROACHES_TO_COMPARE
        slice = subset_df[subset_df.approach.==approach, :]
        println("$approach")
        println("-"^(length(approach) + 4))
        println(
            "  * solve time: ",
            _shifted_geomean(slice.solve_time_sec, time_shift),
            " sec",
        )
        println("  * MIP gap:    ", _shifted_geomean(slice.gap, gap_shift))
        num_winners = winners[approach]
        num_fails = fails[approach]
        num_best_bound = best_bound[approach]
        println("  * Winners:    ", num_winners, " / ", total_instance_count)
        println("  * Fails:      ", num_fails, " / ", total_instance_count)
        println("  * Best bound: ", num_best_bound, " / ", total_instance_count)
    end
end


text:
#**************************************************************************************
# test_GeoEfficiency.jl package
# ============================== part of the GeoEfficiency.jl package.
#
#
#
#**************************************************************************************

# about() test is in the runtests file


text:
#=
1D advection using FV
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("FVadvection1d");

useLog("FVadvection1dlog", level=3)

# Configuration setup
domain(1)
solverType(FV)
timeStepper(RK4)

# Mesh
n = 60 # number of elements
mesh(LINEMESH, elsperdim=n, bids=2)

# Set the order for flux
# finiteVolumeOrder(3);

# Variables and BCs
u = variable("u", location=CELL)
v = variable("v", location=CELL)
w = variable("w", location=CELL)
# boundary(u, 1, FLUX, "t<0.2 ? sin(pi*t/0.2)^2 : 0")
# boundary(u, 2, NO_BC)
# boundary(v, 1, FLUX, "t<0.2 ? sin(pi*t/0.2)^2 : 0")
# boundary(v, 2, NO_BC)
# boundary(w, 1, FLUX, "t<0.2 ? sin(pi*t/0.2)^2 : 0")
# boundary(w, 2, NO_BC)

boundary(u, 1, FLUX, "t<0.2 ? sin(pi*t/0.2)^2 : 0")
boundary(v, 1, DIRICHLET, "t<0.2 ? sin(pi*t/0.2)^2 : 0")
boundary(w, 1, DIRICHLET, 1)
boundary(u, 2, NO_BC)
boundary(v, 2, NO_BC)
boundary(w, 2, DIRICHLET, 0)

# Time interval and initial condition
T = 0.5;
timeInterval(T)
initial(u, "0")
initial(v, "0")
initial(w, "0")

# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)/V) = int(S dx) - int(F ds)
coefficient("a", 1) # advection velocity
# The "upwind" function applies upwinding to the term (a.n)*f with flow velocity a.
# The optional third parameter is for tuning. Default upwind = 0, central = 1. Choose something between these.
flux([u, v, w], ["upwind(a,u)", "upwind(a,v)", "upwind(a,w)"])
# Note that there is no source() for this problem.

# @exportCode("fvad1dcode") # uncomment to export generated code to a file
# @importCode("advec1dcode") # uncomment to import code from a file

solve([u,v,w])

finalize_finch()

##### Uncomment below to compare to exact solution

# # The exact solution with constant velocity v
# a = 1;
# n = size(Finch.fv_info.cellCenters,2);
# exact = zeros(n);
# x = Finch.fv_info.cellCenters[:]
# for i=1:n
#     xt = x[i] - a*T;
#     if xt < 0
#         exact[i] = xt > -0.2 ? sin(pi*-xt/0.2)^2 : 0
#     end
# end

# using Plots
# pyplot();
# display(plot([x x x x], [exact u.values[:] v.values[:] w.values[:]], markershape=:circle, label=["exact" "u" "v" "w"]))

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using RecipesBase

# Plot a 2D IntervalBox:
@recipe function f(xx::IntervalBox{2}) #; customcolor = :green, alpha=0.5)

    (x, y) = xx

    alpha --> 0.5
    seriestype := :shape

    x = [x.lo, x.hi, x.hi, x.lo]
    y = [y.lo, y.lo, y.hi, y.hi]

    x, y
end

# Plot a vector of 2D IntervalBoxes:
@recipe function f(v::Vector{T}) where T<:IntervalBox{2}

    seriestype := :shape

    xs = Float64[]
    ys = Float64[]

    # build up coords:  # (alternative: use @series)
    for xx in v
        (x, y) = xx

        # use NaNs to separate
        append!(xs, [x.lo, x.hi, x.hi, x.lo, NaN])
        append!(ys, [y.lo, y.lo, y.hi, y.hi, NaN])

    end

    alpha --> 0.5

    #x = xs
    #y = ys

    #x, y

    xs, ys

end


text:
using LinearAlgebra
using jInv.Mesh
using Test

print("   test getInterpolationMatrix (RegularMesh) ... ")
domain = zeros(6); domain[2:2:end] = rand(3).+1; domain[1:2:end] = -rand(3).-1
n1    = [2; 3; 4]
n2    = 2*n1
for dim=2:3
	print(" dim=",dim)
	Mc    = getRegularMesh(domain[1:2*dim],n1[1:dim])
	Mf    = getRegularMesh(domain[1:2*dim],n2[1:dim])
	Pcf   = getInterpolationMatrix(Mc,Mf)
	Pfc   = getInterpolationMatrix(Mf,Mc)
	Xc    = getCellCenteredGrid(Mc)
	Xf    = getCellCenteredGrid(Mf)

	@test isapprox(sum(Pfc,dims=2),ones(size(Pfc,1)))
	@test norm(Xc-Pfc*Xf)/norm(Xc) < 1e-13
end
print(" passed\n")
print("   test getInterpolationMatrix (TensorMesh) ... ")
x0    = randn(3)
h1    = rand(4)
h2    = rand(6)
h3    = rand(8)
H1    = h1[1:2:end]+h1[2:2:end]
H2    = h2[1:2:end]+h2[2:2:end]
H3    = h3[1:2:end]+h3[2:2:end]

Mct    = getTensorMesh3D(H1,H2,H3,x0)
Mft    = getTensorMesh3D(h1,h2,h3,x0)
Pcft   = getInterpolationMatrix(Mct,Mft)
Pfct   = getInterpolationMatrix(Mft,Mct)
Xc    = getCellCenteredGrid(Mct)
Xf    = getCellCenteredGrid(Mft)

@test isapprox(sum(Pfct,dims=2),ones(size(Pfct,1)))
@test norm(Xc-Pfct*Xf)/norm(Xc) < 1e-13
print("passed\n")


text:
include("../../MaximinOPF/src/MaximinOPF.jl")
using PowerModels
using JuMP
using Mosek
using MosekTools
PowerModels.silence()

testcase = Dict(
	"file" => "data/case30.m",
 	"name" => "case30K4",
 	"attack_budget" => 4,
 	"inactive_indices" => [],
 	"protected_indices" => []
	)

pm_data = PowerModels.parse_file(testcase["file"])
pm_data["attacker_budget"] = testcase["attack_budget"] ###Adding another key and entry
pm_data["inactive_branches"] = testcase["inactive_indices"] ###Adding another key and entry
pm_data["protected_branches"] = testcase["protected_indices"] ###Adding another key and entry

soc_solver=with_optimizer(Mosek.Optimizer,MSK_IPAR_LOG=0)

for pm_form in [SOCWRConicPowerModel] ### MaximinOPF.conic_supported_pm
    if !(pm_form in MaximinOPF.sdp_pm)
        start_time = time_ns()
        println("Formulating and solving the form ",pm_form)
        maxmin_model=MaximinOPF.MaximinOPFModel(pm_data, pm_form; enforce_int=true)
        JuMP.set_optimizer(maxmin_model,soc_solver)
        JuMP.optimize!(maxmin_model)
        branch_ids=sort(collect(pm_data["undecided_branches"]))
        println("Optimal value using powerform ", pm_form, " is: ",JuMP.objective_value(maxmin_model), " with status ",JuMP.termination_status(maxmin_model))
        x_soln_str = ""
        x_soln = Dict{Int64,Float64}()
        for l in branch_ids
            x_var = variable_by_name(maxmin_model,"x[$l]_1")
            x_val = JuMP.value(x_var)
            if x_val > 1.0-1.0e-8
	            x_val = 1
                x_soln_str = string(x_soln_str," $l")
            elseif x_val < 1.0e-8
	            x_val = 0
            end
            x_soln[l] = x_val
        end
        end_time = time_ns()
        runtime = (end_time-start_time)/1e9
        println("Final best solution: ",x_soln_str," with value ",JuMP.objective_value(maxmin_model))
        println("Runtime: ",runtime)
    end
end


text:
#
println("Hb) Tests of the StrongField module to calculate energy and momentum distributions of photoelectron in SFA.")
#
## setDefaults("print summary: open", "zzz-HHG.sum")

#Pulse shape integrals: Properties of pulses

using DelimitedFiles

if  true
    asfSettings  = AsfSettings(AsfSettings(), generateScf=true)
    grid         = Radial.Grid(Radial.Grid(false), hp=5.0e-2)
    nuclearModel = Nuclear.Model(3.)
    # Compute the initial and final levels for the SFA computation
    wa           = Atomic.Computation(Atomic.Computation(), name="Li ground-state levels", grid=grid, nuclearModel=nuclearModel, configs=[Configuration("[He] 2s")],  asfSettings=asfSettings )
    wb           = perform(wa, output=true)
    initialLevel = wb["multiplet:"].levels[1]
    #
    wa           = Atomic.Computation(Atomic.Computation(), name="He ground-state level", grid=grid, nuclearModel=nuclearModel, configs=[Configuration("[He]")],  asfSettings=asfSettings )
    wb           = perform(wa, output=true)
    finalLevel   = wb["multiplet:"].levels[1]
    #
    #omega        = convertUnits("energy: from wavelength [nm] to atomic", 800.)
    omega        = 0.056939228883964354555047293615643866360187530517578125 #Since we use different conversion factors (precision!) in Mathematica and JAC
    intensity    = convertUnits("intensity: from W/cm^2 to atomic", 1.0e14)
    #A0           = Pulse.computeFieldAmplitude(intensity, omega)
    A0           = 0.4015398801828808927893987856805324554443359375 #Since we use different conversion factors (precision!) in Mathematica and JAC

    observable   = StrongField.SfaEnergyDistribution(pi/2, pi/2, 10, 6*omega) #Arguments are: theta, phi, number of energies, maximum energy
    #observable    = StrongField.SfaMomentumDistribution(pi/2, 360, 5, 10*omega) #Arguments are: theta, number of azimuthal angles, number of energies, masimum energy

    beam         = Pulse.PlaneWaveBeam(A0, omega, 0.)  # (A0, omega, cep)
    #envelope     = Pulse.InfiniteEnvelope()
    envelope     = Pulse.SinSquaredEnvelope(14)
    #polarization = Basics.RightCircular()
    polarization = Basics.RightElliptical(0.5)
    volkov       = StrongField.CoulombVolkov(1)  ## FreeVolkov(), CoulombVolkov(Z) with the charge of the ion Z, DistortedVolkov()

    sfaSettings  = StrongField.Settings()
end
#elseif true

if true
    wa = StrongField.Computation(observable::StrongField.AbstractSFAObservable, nuclearModel, grid, initialLevel, finalLevel, beam, envelope, polarization, volkov, sfaSettings)
    wb = StrongField.perform(wa, output=true)

    distribution = wb["energy distribution"]

    writedlm("Probabilities.csv",distribution.probabilities)

    nothing
end

if true
    StrongField.plot(wa,wb,"omega","linear")
end

#Convergence test for pulse shape integrals
if false
    Tp = 2;
    envelope = Pulse.GaussianEnvelope(Tp * 110.3489708296542008205005779602025351176903820678791359884044538)
    FPlus = ComplexF64[]

    thetap = pi/2
    phip = 0.

    Ep = 0.1138784578
    initialEn = initialLevel.energy

    orderGH = Int[]

    for j = 1:500
        push!(orderGH, j*100)
        push!(FPlus, Pulse.envelopeVolkovIntegral( true, envelope, beam, polarization, thetap, phip, 0.1138784578, initialEn, j*100  ) )
    end

    writedlm("ReFPlus.csv",real(FPlus))
    writedlm("ImFPlus.csv",imag(FPlus))

end


text:
using Test
using SumOfSquares
using DynamicPolynomials

function bivariate_quadratic_test(optimizer,
                                  config::MOIT.TestConfig,
                                  cone::SumOfSquares.PolyJuMP.PolynomialSet)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @variable(model, α)

    @polyvar x
    cref = @constraint(model, x^2 + α*x + 1 in cone)

    # See https://github.com/JuliaOpt/MathOptInterface.jl/issues/676
    @objective(model, Max, α + 1)
    optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ 3.0 atol=atol rtol=rtol

    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test value(α) ≈ 2.0 atol=atol rtol=rtol

    p = gram_matrix(cref)
    @test getmat(p) ≈ ones(2, 2) atol=atol rtol=rtol
    @test p.x == [x, 1]

    @test dual_status(model) == MOI.FEASIBLE_POINT
    μ = dual(cref)
    @test μ isa AbstractMeasure{Float64}
    @test length(moments(μ)) == 3
    a = moment_value.(moments(μ))
    @test a[2] ≈ -1.0  atol=atol rtol=rtol
    @test a[1] + a[3] ≈ 2.0  atol=atol rtol=rtol
    @test monomial.(moments(μ)) == [x^2, x, 1]

    ν = moment_matrix(cref)
    @test getmat(ν) ≈ [a[1] a[2]
                       a[2] a[3]] atol=atol rtol=rtol
    @test ν.x == [x, 1]

    S = SumOfSquares.SOSPolynomialSet{
        SumOfSquares.FullSpace, typeof(cone), SumOfSquares.MonomialBasis,
        Monomial{true},MonomialVector{true},Tuple{}
    }
    @test list_of_constraint_types(model) == [(Vector{AffExpr}, S)]
    test_delete_bridge(
        model, cref, 1,
        ((MOI.VectorOfVariables, MOI.Nonnegatives, 0),
         (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
         (MOI.VectorAffineFunction{Float64},
          SumOfSquares.PolyJuMP.ZeroPolynomialSet{
              SumOfSquares.FullSpace, SumOfSquares.MonomialBasis,
              Monomial{true}, MonomialVector{true}},
          0)))
end
sos_bivariate_quadratic_test(optimizer, config)   = bivariate_quadratic_test(optimizer, config, SOSCone())
sdsos_bivariate_quadratic_test(optimizer, config) = bivariate_quadratic_test(optimizer, config, SDSOSCone())
dsos_bivariate_quadratic_test(optimizer, config)  = bivariate_quadratic_test(optimizer, config, DSOSCone())