text:
"""
    JSONWorksheet

construct 'Array{OrderedDict, 1}' for each row from Worksheet

# Constructors
```julia
JSONWorksheet("Example.xlsx", "Sheet1")
JSONWorksheet("Example.xlsx", 1)

```
# Arguments
- `row_oriented` : if 'true'(the default) it will look for colum names in '1:1', if `false` it will look for colum names in 'A:A'
- `start_line` : starting index of position of columnname.
- `squeeze` : squeezes all rows of Worksheet to a singe row.
- `delim` : a String or Regrex that of deliminator for converting single cell to array.

"""
mutable struct JSONWorksheet
    xlsxpath::String
    pointer::Array{Pointer,1}
    data::Array{T,1} where T
    sheetname::String
end
function JSONWorksheet(xlsxpath, sheet, arr;
                        delim=";", squeeze=false)
    arr = dropemptyrange(arr)
    @assert !isempty(arr) "'$(xlsxpath)!$(sheet)' don't have valid column names, try change optional argument'start_line'"

    pointer = _column_to_pointer.(arr[1, :])
    real_keys = map(el -> el.tokens, pointer)
    # TODO more robust key validity check
    if !allunique(real_keys)
        throw(AssertionError("column names must be unique, check for duplication $(arr[1, :])"))
    end

    if squeeze
        data = squeezerow_to_jsonarray(arr, pointer, delim)
    else
        data = eachrow_to_jsonarray(arr, pointer, delim)
    end
    JSONWorksheet(normpath(xlsxpath), pointer, data, String(sheet))
end
function JSONWorksheet(xf::XLSX.XLSXFile, sheet;
                       start_line=1,
                       row_oriented=true,
                       delim=";", squeeze=false)
    ws = isa(sheet, Symbol) ? xf[String(sheet)] : xf[sheet]
    sheet = ws.name
    # orientation handling
    ws = begin
        rg = XLSX.get_dimension(ws)
        if row_oriented
            rg = XLSX.CellRange(XLSX.CellRef(start_line, rg.start.column_number), rg.stop)
            dt = ws[rg]
        else
            rg = XLSX.CellRange(XLSX.CellRef(rg.start.row_number, start_line), rg.stop)
            dt = permutedims(ws[rg])
        end
    end

    JSONWorksheet(xf.filepath, sheet, ws; delim=delim, squeeze=squeeze)
end
function JSONWorksheet(xlsxpath, sheet; kwargs...)
    xf = XLSX.readxlsx(xlsxpath)
    x = JSONWorksheet(xf, sheet; kwargs...)
    close(xf)
    return x
end

function eachrow_to_jsonarray(data::Array{T,2}, pointers, delim) where T
    json = Array{OrderedDict,1}(undef, size(data, 1) - 1)
    Threads.@threads for i in 1:length(json)
        json[i] = row_to_jsonarray(data[i + 1, :], pointers, delim)
    end
    return json
end

function row_to_jsonarray(row, pointers, delim)
    x = OrderedDict{String,Any}()
    for (col, p) in enumerate(pointers)
        x[p] = collect_cell(p, row[col], delim)
    end
    return x
end

function squeezerow_to_jsonarray(data::Array{T,2}, pointers, delim) where T
    arr_pointer = map(p -> begin
        U = Vector{eltype(p)}; Pointer{U}(p.tokens)
    end, pointers)

    squzzed_json = OrderedDict{String, Any}()
    @inbounds for (col, p) in enumerate(pointers)
        val = map(i -> collect_cell(p, data[i + 1, :][col], delim), 1:size(data, 1) - 1)
        squzzed_json[arr_pointer[col]] = val
    end
    return [squzzed_json]
end

@inline function dropemptyrange(arr::Array{T,2}) where T
    cols = falses(size(arr, 2))
    @inbounds for c in 1:size(arr, 2)
        # There must be a column name, or it's a commet line
        if !ismissing(arr[1, c])
            for r in 1:size(arr, 1)
                if !ismissing(arr[r, c])
                    cols[c] = true
                    break
                end
            end
        end
    end

    arr = arr[:, cols]
    rows = falses(size(arr, 1))
    @inbounds for r in 1:size(arr, 1)
        for c in 1:size(arr, 2)
            if !ismissing(arr[r, c])
                rows[r] = true
                break
            end
        end
    end
    return arr[rows, :]
end

function collect_cell(p::Pointer{T}, cell, delim) where T
    if ismissing(cell)
        val = JSONPointer._null_value(p)
    else
        if T <: AbstractArray
            if isa(cell, AbstractString)
                val = split(cell, delim)
                isempty(val[end]) && pop!(val)
                if eltype(T) <: Real
                    val = parse.(eltype(T), val)
                elseif eltype(T) <: AbstractString
                    val = string.(val)
                end
            else
                val = cell
                if eltype(T) <: Real
                    if isa(cell, AbstractString)
                        val = parse(eltype(T), cell)
                    end
                elseif eltype(T) <: AbstractString
                    if !isa(cell, AbstractString)
                        val = string(cell)
                    end
                end
                val = convert(T, [val])
            end
        else
            val = cell
        end
    end
    return val
end

data(jws::JSONWorksheet) = getfield(jws, :data)
xlsxpath(jws::JSONWorksheet) = getfield(jws, :xlsxpath)
sheetnames(jws::JSONWorksheet) = getfield(jws, :sheetname)
Base.keys(jws::JSONWorksheet) = jws.pointer
function Base.haskey(jws::JSONWorksheet, key::Pointer)
    t = key.tokens
    for el in getfield.(keys(jws), :tokens)
        if el == key.tokens
            return true
        elseif length(el) > length(t)
            if el[1:length(t)] == t
                return true
            end
        end
    end
    return false
end

Base.iterate(jws::JSONWorksheet) = iterate(data(jws))
Base.iterate(jws::JSONWorksheet, i) = iterate(data(jws), i)

Base.size(jws::JSONWorksheet) = (length(jws.data), length(jws.pointer))
function Base.size(jws::JSONWorksheet, d)
    d == 1 ? length(jws.data) :
    d == 2 ? length(jws.pointer) : throw(DimensionMismatch("only 2 dimensions of `JSONWorksheets` object are measurable"))
end
Base.length(jws::JSONWorksheet) = length(data(jws))

StructTypes.StructType(::JSONWorksheet) = StructTypes.ArrayType()

########################################################################
##
## getindex() definitions
##
########################################################################
Base.getindex(jws::JSONWorksheet, i) = getindex(jws.data, i)
Base.getindex(jws::JSONWorksheet, ::Colon, ::Colon) = getindex(jws, eachindex(jws.data), eachindex(jws.pointer))
Base.getindex(jws::JSONWorksheet, row_ind, ::Colon) = getindex(jws, row_ind, eachindex(jws.pointer))

Base.firstindex(jws::JSONWorksheet) = firstindex(jws.data)
Base.lastindex(jws::JSONWorksheet) = lastindex(jws.data)
function Base.lastindex(jws::JSONWorksheet, i::Integer)
    i == 1 ? lastindex(jws.data) :
    i == 2 ? lastindex(jws.pointer) :
    throw(DimensionMismatch("JSONWorksheet only has two dimensions"))
end

function Base.getindex(jws::JSONWorksheet, row_ind::Integer, col_ind::Integer)
    p = keys(jws)[col_ind]

    jws[row_ind, p]
end
function Base.getindex(jws::JSONWorksheet, row_ind::Integer, col_ind::AbstractArray)
    pointers = keys(jws)[col_ind]

    permutedims(map(p -> jws[row_ind, p], pointers))
end
@inline function Base.getindex(jws::JSONWorksheet, row_inds::AbstractArray, col_ind::AbstractArray)
    pointers = keys(jws)[col_ind]
    rows = jws[row_inds]

    # v = vcat(map(el -> jws[el, col_ind], row_inds)...)
    v = Array{Any,2}(undef, length(rows), length(pointers))
    @inbounds for (r, _row) in enumerate(rows)
        for (c, _col) in enumerate(pointers)
            v[r, c] = if haskey(_row, _col)
                _row[_col]
            else
                missing
            end
        end
    end

    return v
end

function Base.getindex(jws::JSONWorksheet, row_ind::Integer, col_ind::Pointer)
    row = jws[row_ind]

    return row[col_ind]
end
@inline function Base.getindex(jws::JSONWorksheet, row_inds, p::Pointer)
    map(row -> row[p], jws[row_inds])
end
@inline function Base.getindex(jws::JSONWorksheet, row_inds, col_ind::Integer)
    p = keys(jws)[col_ind]

    getindex(jws, row_inds, p)
end

function Base.setindex!(jws::JSONWorksheet, value::Vector, p::Pointer)
    if length(jws) != length(value)
        throw(ArgumentError("New column must have the same length as old columns"))
    end
    @inbounds for (i, row) in enumerate(jws)
        row[p] = value[i]
    end
    if !haskey(jws, p)
        push!(jws.pointer, p)
    end
    return jws
end
function Base.setindex!(jws::JSONWorksheet, value, i::Integer, p::Pointer)
    jws[i][p] = value
end

"""
    merge(a::JSONWorksheet, b::JSONWorksheet, bykey::AbstractString)

Construct a merged JSONWorksheet from the given JSONWorksheets.
If the same Pointer is present in another collection, the value for that key will be the
value it has in the last collection listed.
"""
function Base.merge(a::JSONWorksheet, b::JSONWorksheet, key::AbstractString)
    merge(a::JSONWorksheet, b::JSONWorksheet, Pointer(key))
end
function Base.merge(a::JSONWorksheet, b::JSONWorksheet, key::Pointer)
    @assert haskey(a, key) "$key is not found in the JSONWorksheet(\"$(a.sheetname)\")"
    @assert haskey(b, key) "$key is not found in the JSONWorksheet(\"$(b.sheetname)\")"

    pointers = unique([a.pointer; b.pointer])

    keyvalues_a = map(el -> el[key], a.data)
    keyvalues_b = map(el -> el[key], b.data)
    ind = indexin(keyvalues_b, keyvalues_a)

    data = deepcopy(a.data)
    for (i, _b) in enumerate(b.data)
        j = ind[i]
        if isnothing(j)
            _a = deepcopy(_b)
            for p in a.pointer
                _a[p] = JSONPointer._null_value(p)
            end
            push!(data, _a)
        else
            _a = data[j]
        end
        for p in b.pointer
            _a[p] = _b[p]
        end
    end
    JSONWorksheet(a.xlsxpath, pointers, data, a.sheetname)
end
function Base.append!(a::JSONWorksheet, b::JSONWorksheet)
    ak = map(el -> el.tokens, keys(a))
    bk = map(el -> el.tokens, keys(b))

    if sort(ak) != sort(bk)
        throw(AssertionError("""Column names must be same for append!
         $(setdiff(collect(ak), collect(bk)))"""))
    end

    append!(a.data, b.data)
end

function Base.sort!(jws::JSONWorksheet, key; kwargs...)
    sort!(jws, Pointer(key); kwargs...)
end
function Base.sort!(jws::JSONWorksheet, pointer::Pointer; kwargs...)
    sorted_idx = sortperm(map(el -> el[pointer], data(jws)); kwargs...)
    jws.data = data(jws)[sorted_idx]
    return jws
end

function Base.summary(io::IO, jws::JSONWorksheet)
    @printf("%d×%d %s - %s!%s\n", size(jws)..., "JSONWorksheet",
        basename(xlsxpath(jws)), sheetnames(jws))
end
function Base.show(io::IO, jws::JSONWorksheet)
    summary(io, jws)
    # TODO truncate based on screen size
    print(io, "row 1 => ")
    print(io, JSON3.pretty(jws[1]))
    if length(jws) > 1
        print("...")
        print(io, "row $(length(jws)) => ")
        print(io, JSON3.write(jws[end]))
    end
end


text:
@testset "latex-wspace" begin
    s = raw"""
    \newcommand{\hello}{hello}
    A\hello B
    """ |> fd2html_td
    @test isapproxstr(s, "
Ahello B
")
    s = raw"""
    \newcommand{\eqa}[1]{\begin{eqnarray}#1\end{eqnarray}}
    A\eqa{B}C
    \eqa{
        D
    }E
    """ |> fd2html_td
    @test isapproxstr(s, raw"""

            A\[\begin{array}{c} B\end{array}\]C
            \[\begin{array}{c} D\end{array}\]E

""")

    s = raw"""
    \newcommand{\eqa}[1]{\begin{eqnarray}#1\end{eqnarray}}

    \eqa{A\\
        D
    }

    E
    """ |> fd2html_td
    @test isapproxstr(s, raw"""
        \[\begin{array}{c} A\\
        D\end{array}\]

E

        """)

    s = raw"""
    @def indented_code = false
    \newcommand{\eqa}[1]{\begin{eqnarray}#1\end{eqnarray}}

    \eqa{A\\
        D}

    E""" |> fd2html_td
    @test isapproxstr(s, raw"""
        \[\begin{array}{c} A\\
        D\end{array}\]

E

           """)
end

@testset "latex-wspmath" begin
    s = raw"""
    \newcommand{\esp}{\quad\!\!}
    $$A\esp=B$$
    """ |> fd2html_td
    @test isapproxstr(s, raw"\[A\quad\!\!=B\]")
end

@testset "code-wspace" begin
    s = raw"""
    A
    ```
    C
        B
            E
    D
    ```
    """ |> fd2html_td
    @test isapproxstr(s, """
A

C\n    B\n        E\nD
""")
end

@testset "auto html esc" begin
    s = raw"""
    Blah

    ```html
    <div class="foo">Blah</div>
    ```

    End
    """ |> fd2html_td
    @test isapproxstr(s, """

        Blah

<div class="foo">Blah</div>

End
""")
end

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "PCRE2"
version = v"10.36"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://ftp.pcre.org/pub/pcre/pcre2-$(version.major).$(version.minor).tar.gz",
                  "b95ddb9414f91a967a887d69617059fb672b914f56fa3d613812c1ee8e8a1a37"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/pcre2-*/

# Update configure scripts
update_configure_scripts

# Force optimization
export CFLAGS="${CFLAGS} -O3"

# Apply patches
atomic_patch -d src/sljit -p2 ${WORKSPACE}/srcdir/patches/sljit-apple-silicon-support.patch
atomic_patch -d src/sljit -p2 ${WORKSPACE}/srcdir/patches/sljit-nomprotect.patch

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --disable-static \
    --enable-jit \
    --enable-pcre2-16 \
    --enable-pcre2-32

make -j${nproc}
make install

# On windows we need libcpre2-8.dll as well
if [[ ${target} == *mingw* ]]; then
    ln -s libpcre2-8-0.dll  ${libdir}/libpcre2-8.dll
    ln -s libpcre2-16-0.dll ${libdir}/libpcre2-16.dll
    ln -s libpcre2-32-0.dll ${libdir}/libpcre2-32.dll
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libpcre2-8", :libpcre2_8),
    LibraryProduct("libpcre2-16", :libpcre2_16),
    LibraryProduct("libpcre2-32", :libpcre2_32)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    AbstractMaxEntropyModel

An abstract type for a MaxEntropyModel. Each model has one or more structural constraints
that are fixed while the rest of the network is completely random.
"""
abstract type AbstractMaxEntropyModel end

"""
    UBCM{T,N} <: AbstractMaxEntropyModel where {T<:Real, N<:UInt}

# Undirected Binary Configuration Model (UBCM)
Maximum entropy model with a fixed degree sequence. Uses the model where ``x_i = e^{-\\theta_i}``

See also: [AbstractMaxEntropyModel](@ref)
"""
mutable struct UBCM{T,N} <: AbstractMaxEntropyModel where {T<:Real, N<:UInt}
    idx::Vector{N}
    κ::Vector{T}
    f::Vector{T}
    method::Symbol
    x0::Vector{T} # initial guess
    F::Vector{T}  # buffer for solution => used?
    f!::Function  # function that will need to be solved
    xs::Vector{T} # solution vector
    x::IndirectArray{T,1, N, Vector{N}, Vector{T}}  # solution vector expanded
end

Base.show(io::IO, model::UBCM{T,N}) where {T,N} = print(io, """UBCM{$(T)} model ($(length(model.κ) < length(model.idx)  ? "$(round((1- length(model.κ)/length(model.idx)) * 100,digits=2))% compression in $(N)" : "uncompressed"))""")

"""
	UBCM(k::Vector) where {T<:Int}

Generate a UBCM model from degree vector of a graph.
### Input
- `k` degree vector of the graph
- keyword arguments:
    - `compact` (optional, default: `true`) flag to indicate compression usage (as the same constraint leads to the same value)
	- `P` (optional, default: `Float64`) datatype used for storage
	- `method` (optional, default: `:newton`) method to uses to solve the system of equations; should be one of `:newton` or `fixedpoint`
	- `initial` (optional, default: `:degree`) vector of initial values or symbol indication what method tu use for initial guess

### Output
`UBCM` object

### Examples
```Julia
UBCM([1;2;2;3])
````

### See also
[`UBCM`](@ref)
"""
function UBCM(k::Vector{T}; compact::Bool=true,
							P::DataType=Float64,
							method::Symbol=:newton,
							initial::Union{Symbol, Vector}=:degree, kwargs...) where {T<:Int}
	# check input
	method ∈ [:newton, :fixedpoint] ? nothing : throw(ArgumentError("'$(method)' is not a valid solution method"))

	# convert datatypes
	k = convert.(P, k)

	# compress or not
	if compact
		idxtype = length(unique(k)) < typemax(UInt8) ? UInt8 : length(unique(k)) < typemax(UInt16) ? UInt16 : UInt32 # limited to 4.29e9 nodes
		idx, κ = IndirectArray{idxtype}(k).index, IndirectArray{idxtype}(k).values
		f = countmap(k, ones(P, size(k))) # using weight vector for proper type inference
		f = [f[v] for v in κ] # addition
	else
		κ = k
		idxtype = length(κ) < typemax(UInt8) ? UInt8 : length(κ) < typemax(UInt16) ? UInt16 : UInt32 # limited to 4.29e9 nodes
		idx = collect(one(idxtype):idxtype(length(κ)))
		f = ones(P,length(κ)) #     Dict(v => one(P) for v in κ)
	end

	# initial vector - need more work (current conditions are unstable)
	if isa(initial, Vector)
		@warn "THIS IS NOT WORKING!\n$(initial)\n$(initial[idx])\n$(κ)"
		length(initial[idx]) == length(κ) ? nothing : throw(DimensionMismatch("Length of initial vector $(length(initial)) does not match the length of κ $(length(κ))"))
		x0 = P.(initial[idx])
	elseif isa(initial, Symbol)
		initial ∈ [:random, :uniform, :degree, :nodes, :edges] ? nothing : throw(ArgumentError("'$(initial)' is not a valid initial argument method in combination with `$(method)` as a solution method"))
		# should be evaluated more
		if isequal(initial, :random)      		# random in  [0,1]
			x0 = rand(P, length(κ))
		elseif isequal(initial, :uniform) 		# uniform unique random value in [0,1]
			x0 = rand(P) .* ones(P, length(κ))
		elseif isequal(initial, :degree) 		# proportional to degree in [0,1]
			x0 = κ ./ maximum(κ)
		elseif isequal(initial, :nodes) 		# proportional to sqrt of number or nodes
			x0 = κ ./ P.(sqrt(sum(f)))
		elseif isequal(initial, :edges)			# proportional to sqrt of number of edges
			x0 = κ ./ P.(sqrt(sum(f .* κ)/2))
		end
	else
		throw(ArgumentError("`initial` should be either a Vector or a Symbol"))
	end

	# functions to compute
	if isequal(method, :newton)
		f! = (F::Vector, x::Vector) -> UBCM_∂ℒ_∂x!(F, x, κ, f)
	elseif isequal(method, :fixedpoint)
		f! = (F::Vector, x::Vector) -> UBCM_∂ℒ_∂x_it!(F, x, κ, f)
	end

	# buffers (required?)
	F = similar(x0) # modify in place method
	xs = similar(x0) # solution value

	# outresult with indirect indexing for further use
	x = IndirectArray(idx, xs)

	return UBCM(idx, κ, f, method, x0, F, f!, xs, x)
end

"""
	UBCM(G::T; kwargs...) where T<:AbstractGraph

Generate UBCM model directly from a graph.

### See also:
[`UBCM`](@ref)
"""
function UBCM(G::T; kwargs...) where T<:AbstractGraph
	if is_directed(G)
		@warn "Using a directed graph ($(G)) with an undirected network model"
	end

	return  UBCM(degree(G); kwargs...)
end

"""
	fffff(x)

helper function for UBCM gradient computation
"""
fffff(x) = x / (1 + x)

"""
	UBCM_∂ℒ_∂x!(F::Vector, x::Vector, κ::Vector, f::Vector)

Gradient of the likelihood function of the UBCM model using the ``x_{i}`` formulation.

F value of gradients
x value of parameters
κ value of (reduced) degree vector
f frequency associated with each value in κ

### See also:
[`UBCM`](@ref)

"""
function UBCM_∂ℒ_∂x!(F::Vector, x::Vector, k::Vector,f::Vector)
	@tturbo for i in eachindex(x)
		F[i] = -fffff(x[i]*x[i]) - k[i]
		for j in eachindex(x)
			F[i] += f[j] * fffff(x[i]*x[j])
		end
	end

	return F
end

"""
	UBCM_∂ℒ_∂x_it!(F::Vector, x::Vector, κ::Vector, f::Vector)

Iterative gradient of the likelihood function of the UBCM model using the `x_i` formulation.

F value of function
x value of parameters
κ value of (reduce) degree vector
f frequency associated with each value in κ

### See also:
[`UBCM`](@ref)

"""
function UBCM_∂ℒ_∂x_it!(F::Vector{T}, x::Vector{T}, κ::Vector{T}, f::Vector{T}) where {T}
    @tturbo for i in eachindex(x)
		fx = -x[i] / (1 + x[i] * x[i])
		for j in eachindex(x)
			fx += f[j] * x[j] / (1 + x[j] * x[i])
		end

		F[i] = κ[i] / fx
	end

	return F
end

"""
	solve!model::UBCM; , kwargs...)

Solve the equations of a UBCM model to obtain the maximum likelihood parameters.

### See also [`UBCM`](@ref)
"""
function solve!(model::UBCM; kwargs...)
	if isequal(model.method, :newton)
		res = nlsolve(model.f!, model.x0; autodiff=:forward, kwargs...)
	elseif isequal(model.method, :fixedpoint)
		res = fixedpoint(model.f!, model.x0; kwargs...)
	end

	if res.f_converged || res.x_converged
		model.x = res.zero[model.idx]
	else
		@warn "Method `$(model.method)` did not converge. Verify initial conditions or try another method."
	end

	return res
end

"""
    DBCM{T,N} <: AbstractMaxEntropyModel where {T<:Real, N<:UInt}

# Directed Binary Configuration Model (DBCM)
Maximum entropy model with a fixed in- and outdegree sequence.
"""
mutable struct DBCM{T,N} <: AbstractMaxEntropyModel where {T<:Real, N<:UInt}

end

"""
	UBCM_∂ℒ_∂x!(F::Vector, x::Vector, κ::Vector, f::Vector)

Gradient of the likelihood function of the UBCM model using the ``x{i}`` formulation. Used to generate a function
`(F, x) -> UBCM_∂ℒ_∂x!(F, x, κ, f)` by a `UBCM` instance.

F value of gradients
x value of parameters
κ value of (reduced degree vector)
f frequency associated with each value in κ

## see also:
[`UBCM`](@ref)

"""
function UBCM_∂ℒ_∂θ!(F::Vector, x::Vector, k::Vector,f::Vector)
	@tturbo for i in eachindex(x)
		F[i] = -fffff(x[i]*x[i]) - k[i]
		for j in eachindex(x)
			F[i] += f[j] * fffff(x[i]*x[j])
		end
	end
	return F
end


text:
using Test
using AtomicGraphNets
using ChemistryFeaturization
using Flux

using Pkg

# define a test model input
dummyfzn = GraphNodeFeaturization(["X" for i = 1:40], nbins = 1) # cheeky way to get a matrix of all 1's
ag = AtomGraph(Float64.([0 1; 1 0]), ["H", "O"])
input = featurize(ag, dummyfzn)

@testset "CGCNN" begin
    # build a model with deterministic initialization of weights
    in_fea_len = 40
    conv_fea_len = 20
    pool_fea_len = 10
    model = build_CGCNN(
        in_fea_len,
        atom_conv_feature_length = conv_fea_len,
        pooled_feature_length = pool_fea_len,
        num_hidden_layers = 2,
        initW = (d...) -> fill(0.1, d...),
    )

    @testset "initialization" begin
        @test length(model) == 5
        @test size(model[1].convweight) ==
              size(model[1].selfweight) ==
              (conv_fea_len, in_fea_len)
        @test size(model[2].convweight) ==
              size(model[2].selfweight) ==
              (conv_fea_len, conv_fea_len)
        @test size(model[4].weight) == (pool_fea_len, pool_fea_len)
        @test size(model[5].weight) == (1, pool_fea_len)
    end

    @testset "forward pass" begin
        lapl, output1 = model[1](input)
        int_mat =
            model[2].convweight * output1 * lapl +
            model[2].selfweight * output1 +
            hcat([model[2].bias for i = 1:size(output1, 2)]...)
        @test all(isapprox.(model[1:2](input)[2], zeros(Float64, 20, 2), atol = 2e-3))
        @test model(input)[1] ≈ 0.693 atol = 1e-3
    end

    @testset "backward pass" begin
        loss(x, y) = Flux.Losses.mse(model(x), y)
        data = zip([input], [1.0]) # high-quality data...
        opt = Descent(0.1) # should be deterministic
        Flux.train!(loss, params(model), data, opt)
        @test model(input)[1] ≈ 1.051 atol = 1e-3
    end
end

@testset "SGCNN" begin
    # could probably do with some more detailed tests here but better something than nothing for now
    model = build_SGCNN(40, initW = (d...) -> fill(0.1, d...))

    @testset "initialization" begin
        @test length(model) == 5
        @test model[1].connection == vcat
        @test length.(model[1].layers) == (3, 3)
    end

    @testset "forward pass" begin
        @test model((input, input))[1] ≈ 45.331 atol = 1e-3
    end

    @testset "backward pass" begin
        loss(x, y) = Flux.Losses.mse(model(x), y)
        data = zip([(input, input)], [1.0]) # high-quality data...
        opt = Descent(0.1) # should be deterministic
        Flux.train!(loss, params(model), data, opt)
        @test model((input, input))[1] ≈ -75.985 atol = 1e-3
    end
end

@testset "DEQ" begin
    # build a model with deterministic initialization of weights
    in_fea_len = 40
    pool_fea_len = 10
    model = build_CGCNN_DEQ(
        in_fea_len,
        pooled_feature_length = pool_fea_len,
        num_hidden_layers = 2,
        initW = (d...) -> fill(0.1, d...),
    )

    @testset "initialization" begin
        @test length(model) == 4
        @test size(model[1].conv.convweight) ==
              size(model[1].conv.selfweight) ==
              (in_fea_len, in_fea_len)
        @test size(model[3].weight) == (pool_fea_len, pool_fea_len)
        @test size(model[4].weight) == (1, pool_fea_len)
    end

    @testset "forward pass" begin
        @test model(input)[1] ≈ 0.693147 atol = 1e-6
    end

    @testset "backward pass" begin
        loss(x, y) = Flux.Losses.mse(model(x), y)
        data = zip([input], [1.0]) # high-quality data...
        opt = Descent(0.1) # should be deterministic
        Flux.train!(loss, params(model), data, opt)
        @test model(input)[1] ≈ 1.051563 atol = 1e-6
    end
end


text:
is_skewsymm(M) = M' == -M

col_negs(R::Matrix) = [sum(R[:,j].<0) for j=1:size(R)[2]]
row_negs(R::Matrix) = [sum(R[i,:].<0) for i=1:size(R)[1]]
negs(M) = (row_negs(M),col_negs(M))

straight_cols(R) = all(sort(col_negs(R)).==Vector(0:size(R)[2]-1))
straight_rows(R) = all(sort(row_negs(R)).==Vector(0:size(R)[1]-1))
is_straight(R) = straight_cols(R) && straight_rows(R)

crossing_cols(R) = any(sort(col_negs(R)).>Vector(0:size(R)[2]-1))
crossing_rows(R) = any(sort(row_negs(R)).>Vector(0:size(R)[1]-1))
is_crossing(R) = crossing_cols(R) || crossing_rows(R)

sum_above_cols(R) = sum(sort(col_negs(R)).>Vector(0:size(R)[2]-1))
sum_above_rows(R) = sum(sort(row_negs(R)).>Vector(0:size(R)[1]-1))

function upper_sum(A::Matrix)
    S, (H,W) = 0, size(A)
    @assert W == H
    for i=1:H, j=i+1:W  S += A[i,j] end
    return S
end

is_min_reached(R,T) = upper_sum(R[T,T]) == sum(R[R.<0])

max_stars_repo_name: vaa-msu/ranking-tests

max_issues_repo_name: vaa-msu/ranking-tests

max_forks_repo_name: vaa-msu/ranking-tests


text:
println("Start")
# Naddprocs = 1
# if Naddprocs > 1
#     addprocs(Naddprocs)
#
#     println("Load module")
#     @everywhere push!(LOAD_PATH, ".")
#     @everywhere using RobustOptimization, ScikitLearn
#     using Ipopt, PyPlot
#     @everywhere @sk_import model_selection: train_test_split
#     @everywhere @sk_import metrics: r2_score
#
#     @everywhere nbfeatures = 13
#     solver = IpoptSolver(print_level = 2)
#
#     println("Create date")
#     @everywhere df_train, df_test = create_data("housing_scale", nbfeatures, 0.6)
# else
println("Load module")
push!(LOAD_PATH, ".")
using RobustOptimization
using ScikitLearn
using Ipopt, PyPlot
@sk_import model_selection: train_test_split
@sk_import metrics: r2_score

nbfeatures = 13
solver = IpoptSolver(print_level = 2)

println("Create data")
df_train, df_test = create_data("housing_scale", nbfeatures, 0.1, train_test_split)

println("Define model")
N = size(df_train)[1]
ϵ = 0.1
verbosity = 100
itmax = 5000
sample = 4

ambiguity = "wasserstein"
robustModel = RobustModel(N, nbfeatures, ϵ, ambiguity, LinearRegression())
α = 1/norm(robustModel.descent_direction)
projParams = ProjParams(Int(1e6), 1e-5, sample, para_proj=Sequential(), para_inter=Sequential())
optParams = OptParams(itmax, 1e-7, α, verbosity = verbosity);

println("Solve normal")
xnormal = normal_opt(df_train, solver)
y_pred = pred(df_test, xnormal, robustModel.regressionModel)
y_true = df_test[:,end]
println(" ")
println(" ")
println("xnormal = ", xnormal)
println("rnormal = ", r2_score(y_true, y_pred))
println(" ")
println(" ")

println("Solve wasserstein")
x0 = init_proj(df_train, robustModel, projParams)
println("fin_init")
xalg, yalg, dmwas, mem, mini = run_algo(x0, df_train, robustModel, optParams, projParams)
xrobust1 = getsolution(xalg, ambiguity, nbfeatures)
xrobust2 = getsolution(mem, ambiguity, nbfeatures)
y_pred1 = pred(df_test, xrobust1, robustModel.regressionModel)
y_pred2 = pred(df_test, xrobust2, robustModel.regressionModel)
y_true = df_test[:,end]
r1 = r2_score(y_true, y_pred1)
r2 = r2_score(y_true, y_pred2)
if r1 > r2
    println(" ")
    println(" ")
    println("xrobwas = ", xrobust1)
    println("r1 = ", r2_score(y_true, y_pred1))
    println(" ")
    println(" ")
else
    println(" ")
    println(" ")
    println("xrobwas = ", xrobust2)
    println("r2 = ", r2_score(y_true, y_pred2))
    println(" ")
    println(" ")
end
println("Fin wasserstein")

ambiguity = "KLdivergence"
robustModel = RobustModel(N, nbfeatures, ϵ, ambiguity, LinearRegression())
α = 1/norm(robustModel.descent_direction)
projParams = ProjParams(Int(1e6), 1e-5, sample, para_proj=Sequential(), para_inter=Sequential())
optParams = OptParams(itmax, 1e-7, α, verbosity = verbosity);

println("Solve KL")
x0 = init_proj(df_train, robustModel, projParams)
println("fin_init")
xalg, yalg, dmKL, mem, mini = run_algo(x0, df_train, robustModel, optParams, projParams)
xrobust1 = getsolution(xalg, ambiguity, nbfeatures)
xrobust2 = getsolution(mem, ambiguity, nbfeatures)
y_pred1 = pred(df_test, xrobust1, robustModel.regressionModel)
y_pred2 = pred(df_test, xrobust2, robustModel.regressionModel)
y_true = df_test[:,end]
r1 = r2_score(y_true, y_pred1)
r2 = r2_score(y_true, y_pred2)
if r1 > r2
    println(" ")
    println(" ")
    println("xrobKL = ", xrobust1)
    println("r1 = ", r2_score(y_true, y_pred1))
    println(" ")
    println(" ")
else
    println(" ")
    println(" ")
    println("xrobKl = ", xrobust2)
    println("r2 = ", r2_score(y_true, y_pred2))
    println(" ")
    println(" ")
end

println("Fin KL")

println("dmwas = ", dmwas)
println("dmKl = ", dmKL)

max_stars_repo_name: Henri-Gerard/robox

max_issues_repo_name: Henri-Gerard/robox

max_forks_repo_name: Henri-Gerard/robox


text:
@testset "Linear disjoint" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  _K, _ = NumberField([x^2 - 2, x^2 - 3], "a", cached = false)
  K, _ = simple_extension(_K)
  L, b = NumberField(x^2 - 2, "b", cached = false)
  @test !islinearly_disjoint(K, L)

  M, c = NumberField(x^2 - 3, "c", cached = false)
  @test islinearly_disjoint(L, M)
end

@testset "Random" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x^32 + 2, "a")

  b = @inferred rand([a], -10:10)
  @test b isa nf_elem
  @test_throws ErrorException rand(nf_elem[], -10:10)

  b = @inferred rand(basis(K), 1:100, 10)
  @test count(!iszero, (coeff(b, i) for i in 0:31)) <= 10
  @test_throws ErrorException rand(nf_elem[], -10:10, 5)
  @test_throws ErrorException rand([a, a^2], -10:10, -10)
  @test_throws ErrorException rand(basis(K), -10:10, 100)

  @inferred rand!(b, basis(K), 1:100, 20)
  @test count(!iszero, (coeff(b, i) for i in 0:31)) <= 20
  @test_throws ErrorException rand!(b, basis(K), 1:100, 120)
  @test_throws ErrorException rand!(b, basis(K), 1:100, -100)

  @inferred rand!(b, basis(K), 1:100)
  @test_throws ErrorException rand!(b, nf_elem[], 1:100)
end

@testset "Polynomial" begin
  Qx, x = QQ["x"]
  K, _a = NumberField(x^3 - 3*x - 1, "a")
  Kt, t = K["t"]
  f = t^4+(-28*_a^2 + 26*_a + 124)*t^2+(81*_a^2 + 504*_a + 936)
  @test @inferred isirreducible(f)

  f = x^3-7*x^2+6*x-1
  K, a = number_field(f, "a")
  Ky, y = K["y"]
  h = y^3+(15037//140*a^2 - 109//40*a - 915//14)*y^2+(16375724527//78400*a^2 - 993527643//4900*a + 393774209//11200)*y+(107296943419277//878080*a^2 - 2594040461688323//21952000*a + 17784340885567//878080)
  @assert isirreducible(h)
end

@testset "Is integral" begin
  Qx, x = FlintQQ["x"]
  f = x^2 + 1
  K, a = number_field(f, "a")

  @test Hecke.isintegral(a) == true
  @test Hecke.isintegral(fmpq(1, 2)*a) == false

  g = x^3 + 3
  L, b = number_field([f, g], "b")

  @test Hecke.isintegral(b[1]) == true
  @test Hecke.isintegral(fmpq(1, 2)*b[1]) == false
end

@testset "Compositum" begin
  Qx, x = FlintQQ["x"]
  f = x^2 + 1
  K, a = number_field(f, "a")
  L, b = number_field(x^2-3, "b")
  C, mK, mL = compositum(K, L)

  @test iszero(K.pol(mK(gen(K))))
  @test iszero(L.pol(mL(gen(L))))
end

@testset "PolyFactor" begin
  Zx, x = FlintZZ["x"]
  k, a = number_field(swinnerton_dyer(3, x))
  kt, t = k["t"]

  g = swinnerton_dyer(8, x)
  @test length(factor((t^2-a)*(t^3-a-1))) == 2 #Trager
  @test length(factor((t^2-a)*(t^3-a-1)*(t+a^2+1)*(t^5+t+a))) == 4 #Zassenhaus
  @test length(factor(change_base_ring(k, g))) == 8 # van Hoeij
  @test length(factor(t)) == 1
  @test length(factor(t^10)) == 1

  # Now the same with non-nice defining polynomial

  for i in 1:10
    n = rand(1:10)
    d = rand(1:10)
    k, a = number_field(n//d * change_base_ring(FlintQQ, swinnerton_dyer(3, x)))
    kt, t = k["t"]
    g = swinnerton_dyer(8, x)
    @test length(factor((t^2-a)*(t^3-a-1))) == 2 #Trager
    @test length(factor((t^2-a)*(t^3-a-1)*(t+a^2+1)*(t^5+t+a))) == 4 #Zassenhaus
    @test length(factor(change_base_ring(k, g))) == 8 # van Hoeij
    @test length(factor(t)) == 1
    @test length(factor(t^10)) == 1
  end

  K, a = NumberField(x - 1, "a")
  Kt, t = K["t"]
  f = t^5 -3 * t^4 - 104 * t^3 + 312 * t^2 + 400*t -1200
  @test length(factor(f)) == 5
  @test length(factor(f*t)) == 6

  for i in 1:10
    n = rand(1:10)
    d = rand(1:10)
    K, a = NumberField(n//d * change_base_ring(FlintQQ, x - 1), "a")
    Kt, t = K["t"]
    f = t^5 -3 * t^4 - 104 * t^3 + 312 * t^2 + 400*t -1200
    @test length(factor(f)) == 5
    @test length(factor(f*t)) == 6

    K, a = NumberField(change_base_ring(FlintQQ, x) - n//d, "a")
    Kt, t = K["t"]
    f = t^5 -3 * t^4 - 104 * t^3 + 312 * t^2 + 400*t -1200
    @test length(factor(f)) == 5
    @test length(factor(f*t)) == 6
  end

  #Tommys
  K, a = number_field(x^2 - x - 4)
  Ky, y = K["y"]
  f = y^16+(39)*y^14+(449)*y^12+(1794)*y^10+(2830)*y^8+(1794)*y^6+(449)*y^4+(39)*y^2+(1)
  @test length(factor(f)) == 2

  # Dan's example #76
  K, a = number_field(x^2 + x + 1, cached = false);
  Kt, t = K["t"]
  g = t^4-t^3-t+(1)
  @test length(factor(g)) == 3

  f = x^4 + 24*x^2+28
  K, a = number_field(f, "a", cached = false);
  g = x^8-1
  @test length(factor(g, K)) == 4

  K, a = number_field(x^3 - 2, "a")
  Kx, x = PolynomialRing(K, "x")
  fa = factor(x^4+(1//2*a^2 + 1*a + 2)*x^3+(a^2 + 2*a + 2)*x^2+(1//2*a^2 + 1*a + 2)*x+1)
  @test length(fa) == 3
end

@testset "Root computation" begin
  Qx, x = QQ["x"]
  f = x^3-39*x-65
  K, a = number_field(f, "a")
  r = @inferred roots(f, K)
  @test length(r) == 3
  @test all(iszero, (f(b) for b in r))
end

@testset "Norm of factored element" begin
  Qx, x = QQ["x"]
  f = x^3 - 2
  K, a = number_field(f, "a")
  Kt, t = K["t"]
  L, b = number_field(x^6 - 6*x^4 - 4*x^3 + 12*x^2 - 24*x - 4)
  h = hom(K, L, -12//155*b^5 - 9//310*b^4 + 16//31*b^3 + 78//155*b^2 - 76//155*b + 182//155)
  z = -3*b^5 - 4*b^3 + 3*b^2 + 6*b + 2
  @test norm(h, z) == 1620*a^2 - 3342*a + 1120
  @test evaluate(norm(h, FacElem(z))) == evaluate(FacElem(K, Dict(1620*a^2 - 3342*a + 1120 => 1)))
  @test evaluate(norm(h, FacElem(L))) == evaluate(FacElem(K))
end


text:
module TestIAEA
using Test
using PhaseSpaceIO
using PhaseSpaceIO: position
using PhaseSpaceIO.Testing
using Accessors
using PhaseSpaceIO: ptype, read_particle, write_particle, CompressedIAEAParticle
using StaticArrays
using LinearAlgebra
using CoordinateTransformations

@testset "CompressedIAEAParticle" begin
    path = testdatapath("some_file.IAEAphsp")
    p = first(PhspVector(path))

    # does not throw
    show(IOBuffer(), p)
    show(IOBuffer(), [p])
    p2 = IAEAParticle(p)
    @test p isa CompressedIAEAParticle
    @test p2 isa IAEAParticle
    @test propertynames(p) == propertynames(p2)
    for prop in propertynames(p)
        @test getproperty(p, prop) == getproperty(p2, prop)
    end
    @test p == p2
    @test [p] == [p2]
    pos = @SVector randn(3)
    dir = normalize(@SVector randn(3))
    @test (@set position(p) = pos) == (@set position(p2) = pos)
    @test (@set direction(p) = dir) == (@set direction(p2) = dir)

    trafo = AffineMap(@SMatrix(randn(3,3)), @SVector(randn(3)))
    @test trafo(p) == trafo(p2)
end

@testset "IAEA save load with many constants" begin
    mktempdir() do dir
        h = IAEAHeader{0,0}((x=1, y=2, z=3, u=0, v=1, weight=5))
        p_truth = IAEAParticle(x=1,y=2,z=3,u=0,v=1,w=0,weight=5, typ=photon, E=6)
        path = IAEAPath(joinpath(dir, "phsp"))
        iaea_writer(path, h) do w
            write(w, p_truth)
        end

        r_reloaded, ps_reloaded = phsp_iterator(path) do iter
            iter.header.record_contents, collect(iter)
        end
        @test r_reloaded == h.record_contents
        @test length(ps_reloaded) == 1
        @test first(ps_reloaded) == p_truth
        @test eltype(ps_reloaded) <: IAEAParticle

        v = PhspVector(path)
        @test length(v) == 1
        @test eltype(v) <: CompressedIAEAParticle
        @test first(v) == p_truth
    end
end

@testset "similar_header" begin
    path = testdatapath("some_file.IAEAphsp")
    h = phsp_iterator(PhaseSpaceIO.similar_header, path)
    @test typeof(h) == typeof(IAEAHeader{0,1}())
    ps = phsp_iterator(collect, path)
    h  = PhaseSpaceIO.similar_header(ps)
    @test typeof(h) == typeof(IAEAHeader{0,1}())
end

@testset "finalizer" begin
    function write_sloppy(path, ps, r)
        w = iaea_writer(path, r)
        for p in ps
            write(w, p)
        end
    end
    path = IAEAPath(tempname())
    ps = [arbitrary(IAEAParticle{0,1})]
    r = IAEAHeader{0,1}()
    write_sloppy(path, ps, r)
    GC.gc()
    ps2 = phsp_iterator(collect, path)
    @test length(ps) == length(ps2) == 1
    @test first(ps) ≈ first(ps2)
    rm(path)
end

@testset "read write single particle" begin
    h = IAEAHeader{0,1}()
    P = ptype(h)
    @test P == IAEAParticle{0,1}
    p_ref = P(photon,
        1.0f0, 2.0f0,
        3.0f0, 4.0f0, 5.0f0,
        0.53259337f0, 0.3302265f0, -0.7792912f0,
        true, (), (13,))

    path = testdatapath("some_file.IAEAphsp")
    ps = iaea_iterator(collect, path)
    @test length(ps) == 1
    @test first(ps) == p_ref
    iaea_iterator(path) do iter
        @test iter.header.attributes[:INSTITUTION] == "hello world"
        @test iter.header.attributes[:ORIG_HISTORIES] == "18446744073709551615"
    end

    io = IOBuffer()
    write_particle(io, p_ref, h)
    seekstart(io)
    p = @inferred read_particle(io, h)
    @test p === p_ref
    @test eof(io)
end

@testset "test IAEAPhspIterator" begin
    path = testdatapath("some_file.IAEAphsp")
    phsp = iaea_iterator(path)
    @test length(phsp) == 1
    @test eltype(phsp) === IAEAParticle{0,1}
    @test collect(phsp) == collect(phsp)
    @test length(collect(phsp)) == 1
    close(phsp)
end

function test_header_contents(path)
    s = read(path, String)
    for key in [
        :IAEA_INDEX,
        :TITLE,
        :FILE_TYPE,
        :CHECKSUM,
        :RECORD_LENGTH,
        :BYTE_ORDER,
        :ORIG_HISTORIES,
        :PARTICLES,
        :TRANSPORT_PARAMETERS,
        :MACHINE_TYPE,
        :MONTE_CARLO_CODE_VERSION,
        :GLOBAL_PHOTON_ENERGY_CUTOFF,
        :GLOBAL_PARTICLE_ENERGY_CUTOFF,
        :COORDINATE_SYSTEM_DESCRIPTION,
        :BEAM_NAME,
        :FIELD_SIZE,
        :NOMINAL_SSD,
        :MC_INPUT_FILENAME,
        :VARIANCE_REDUCTION_TECHNIQUES,
        :INITIAL_SOURCE_DESCRIPTION,
        :PUBLISHED_REFERENCE,
        :AUTHORS,
        :INSTITUTION,
        :LINK_VALIDATION,
        :ADDITIONAL_NOTES,]
        @test occursin(string(key), s)
    end
end

@testset "test iaea_iterator iaea_writer" begin
    for _ in 1:100
        f = rand(1:3)
        i = rand(1:3)
        n = rand(1:1000)
        ps = [arbitrary(IAEAParticle{f,i}) for _ in 1:n]
        r = IAEAHeader{f,i}()
        dir = tempname()
        mkpath(dir)
        path = IAEAPath(joinpath(dir, "hello"))
        iaea_writer(path,r) do w
            for p in ps
                write(w, p)
            end
        end
        @test ispath(path.header)
        @test ispath(path.phsp)
        test_header_contents(path.header)

        ps_reload = iaea_iterator(collect, path)
        @test all(ps_reload .≈ ps)
        @test eltype(ps_reload) === eltype(ps)
        rm(path)

        @test !ispath(path.header)
        @test !ispath(path.phsp)
    end
end

end #module


text: RpcApi.Statistics(150,16,681659,900,600,ASCIIString["orgMutationProbs: [100,300,95,0,1,1,1]","orgCloneMutation: 1","ORGANISM_MAX_MUTATIONS_ON_CLONE: 100","orgRainMutationPeriod: 500","ORGANISM_MAX_MUTATION_PERIOD: 10000","ORGANISM_MUTATION_AMOUNT: 2","ORGANISM_MAX_MUTATION_AMOUNT: 100","orgStartAmount: 250","orgStartEnergy: 5000","ORGANISM_MAX_ENERGY: 4294967295","orgEnergySpendPeriod: 200","orgEnergySpendAmount: 1","orgRemoveWeakPeriod: 500","orgRemoveWeakAmount: 5","orgClonePeriod: 5","orgStartColor: 1","codeFuncParamAmount: 2","worldWidth: 900","worldHeight: 600","WORLD_FRAME_DELAY: 0","worldIps: 16","worldMaxOrgs: 400","worldMinOrgs: 50","worldStartEnergyDots: 10000","worldStartEnergyInDot: 500","worldEnergyCheckPercent: 0.3","worldEnergyCheckPeriod: 500","worldZoom: 1","backupPeriod: 1","backupAmount: 7","STAT_WIDTH: 650","STAT_HEIGHT: 500","STAT_FRAME_DELAY: 5","conStepEnergySpendPercent: 20","conServerPort: 2000","conFastServerPort: 2001","conServerIp: 127.0.0.1","conLeftServerPort: 0","conLeftServerIp: 127.0.0.1","conRightServerPort: 0","conRightServerIp: 127.0.0.1","conUpServerPort: 0","conUpServerIp: 127.0.0.1","conDownServerPort: 0","conDownServerIp: 127.0.0.1"],RpcApi.SimpleOrganism(0x00000000000a6648,:(function (o::Creature.Organism,)

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
struct LinearCache{TA,Tb,Tu,Tp,Talg,Tc,Tl,Tr,Ttol}
    A::TA
    b::Tb
    u::Tu
    p::Tp
    alg::Talg
    cacheval::Tc  # store alg cache here
    isfresh::Bool # false => cacheval is set wrt A, true => update cacheval wrt A
    Pl::Tl        # preconditioners
    Pr::Tr
    abstol::Ttol
    reltol::Ttol
    maxiters::Int
    verbose::Bool
end

"""
$(SIGNATURES)
"""
function set_A(cache::LinearCache, A)
    @set! cache.A = A
    @set! cache.isfresh = true
    return cache
end

"""
$(SIGNATURES)
"""
function set_b(cache::LinearCache, b)
    @set! cache.b = b
    return cache
end

"""
$(SIGNATURES)
"""
function set_u(cache::LinearCache, u)
    @set! cache.u = u
    return cache
end

"""
$(SIGNATURES)
"""
function set_p(cache::LinearCache, p)
    @set! cache.p = p
#   @set! cache.isfresh = true
    return cache
end

"""
$(SIGNATURES)
"""
function set_prec(cache, Pl, Pr)
    @set! cache.Pl = Pl
    @set! cache.Pr = Pr
    return cache
end

function set_cacheval(cache::LinearCache, alg_cache)
    if cache.isfresh
        @set! cache.cacheval = alg_cache
        @set! cache.isfresh = false
    end
    return cache
end

init_cacheval(alg::SciMLLinearSolveAlgorithm, args...) = nothing

SciMLBase.init(prob::LinearProblem, args...; kwargs...) = SciMLBase.init(prob,nothing,args...;kwargs...)

default_tol(::Type{T}) where T = √(eps(T))
default_tol(::Type{Complex{T}}) where T = √(eps(T))
default_tol(::Type{<:Rational}) = 0
default_tol(::Type{<:Integer}) = 0
default_tol(::Type{Any}) = 0

function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorithm,Nothing}, args...;
                        alias_A = false, alias_b = false,
                        abstol=default_tol(eltype(prob.A)),
                        reltol=default_tol(eltype(prob.A)),
                        maxiters=length(prob.b),
                        verbose=false,
                        Pl = Identity(),
                        Pr = Identity(),
                        kwargs...,
                       )
    @unpack A, b, u0, p = prob

    u0 = if u0 !== nothing
        u0
    else
        u0 = similar(b, size(A, 2))
        fill!(u0,false)
    end

    cacheval = init_cacheval(alg, A, b, u0, Pl, Pr, maxiters, abstol, reltol, verbose)
    isfresh = true
    Tc = typeof(cacheval)

    A = alias_A ? A : deepcopy(A)
    b = alias_b ? b : deepcopy(b)

    cache = LinearCache{
        typeof(A),
        typeof(b),
        typeof(u0),
        typeof(p),
        typeof(alg),
        Tc,
        typeof(Pl),
        typeof(Pr),
        typeof(reltol),
    }(
        A,
        b,
        u0,
        p,
        alg,
        cacheval,
        isfresh,
        Pl,
        Pr,
        abstol,
        reltol,
        maxiters,
        verbose,
    )
    return cache
end

SciMLBase.solve(prob::LinearProblem, args...; kwargs...) = solve(init(prob, nothing, args...; kwargs...))

SciMLBase.solve(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorithm,Nothing},
                args...; kwargs...) = solve(init(prob, alg, args...; kwargs...))

SciMLBase.solve(cache::LinearCache, args...; kwargs...) =
    solve(cache, cache.alg, args...; kwargs...)


text:
import JSON.print
import JSON.Parser.parse
using NetCDF
using NCDatasets
using DataStructures

indent = 4 # figure out how to read user defined indent

"""
    ncmeta_addvarattr(filename_ncfile::String, filename_json::String)

Read user-definied JSON file with defined global meta-data information and apply to existing netcdf file
"""
function ncmeta_addvarattr(var :: String, filename_ncfile::String, filename_json::String)

    # Open the JSON file and get the global attributes
    fio = open(filename_json,"r")
    varattrinput = parse(fio,dicttype=OrderedDict)
    close(fio)

    # Read the netcdf files and update the global data
    ncds = NCDataset(filename_ncfile,"a")
    for (attrname,attrdesc) in varattrinput
        ncds[var].attrib[attrname] = attrdesc
    end
    close(ncds)
    return nothing
end

"""
    ncmeta_addglobalattr(filename_ncfile::String, filename_json::String)

Read user-definied JSON file with defined global meta-data information and apply to existing netcdf file
"""
function ncmeta_addglobalattr(filename_ncfile::String, filename_json::String)

    # Open the JSON file and get the global attributes
    fio = open(filename_json,"r")
    globalinput = parse(fio,dicttype=OrderedDict)
    close(fio)

    # Read the netcdf files and update the global data
    ncds = NCDataset(filename_ncfile,"a")
    for (attrname,attrdesc) in globalinput
        ncds.attrib[attrname] = attrdesc
    end
    close(ncds)
    return nothing
end

"""
    ncmeta_globaltemplate(filename_template::String, filename_ncfile::String="")

Define a JSON template file with necessary meta-date fields for global attributes.  If `filename_ncfile` is provided, grab the global attributes already included and add any missing default global attributes.
"""
function ncmeta_globaltemplate(filename_template::String, filename_ncfile::String="")

    # If a netcdf file is given, read the attribute meta-data into an OrderedDict
    if filename_ncfile != ""
        ncds = NCDataset(filename_ncfile)
        glbattr = OrderedDict(ncds.attrib)
        tempglobal = merge(globaltemplate(),glbattr) # Merge with the import holding priority
    else
        tempglobal = globaltemplate()
    end
    close(ncds)
    # TO DO: Get the variable attribute metadata and merge with global

    # Write the template to a file
    open(filename_template,"w") do fio
        print(fio, tempglobal, indent)
    end
    return nothing
end

"""
    ncmeta_vartemplate(filename_template::String, filename_ncfile::String="")

Define a JSON template file with necessary meta-date fields for global attributes.  If `filename_ncfile` is provided, grab the global attributes already included and add any missing default global attributes.
"""
function ncmeta_vartemplate(filename_template::String)

    # Write the template to a file
    open(filename_template,"w") do fio
        print(fio, vartemplate(), indent)
    end
    return nothing
end

"""
    ncmeta_cfcheck(filename :: String)

Check the user-defined JSON file for CF conformance (min requirements?)
"""
function ncmeta_cfcheck(filename :: String)
    return nothing
end

"""
    ncmeta_checkvar(filename :: String)

Check that the user-defined variables have attributes
"""

function ncmeta_checkvar(filename :: String)

    # Check that the file has variablesd
    ncds = NCDataset(filename)
    if isempty(keys(ncds))
        error("The dataset in ", filename, "does not contain any variables")
    end

    # Report which keys if any have variables
    for (varname,var) in ncds
        println("Checking ", varname)

        # Check if variable has any attributes
        if isempty(var.attrib)
            println(varname," is missing attribute data")
        end

        # TO DO: cross reference against standard variable metadata?

    end
    return nothing

end

"""
    globaltemplate()

Function that defines the minimum set of global attributes as a dict
"""
function globaltemplate()

    # Define the minimum global attribute dict
    globaldict = OrderedDict("comment" => "",
                             "conventions" => "",
                             "history" => "",
                             "institution" => "",
                             "references" => "",
                             "source" => "",
                             "title" => ""
                             )
    return globaldict
end

"""
    vartemplate()

Function that defines the minimum set of var attributes as a dict
"""
function vartemplate()

    # Define the minimum variable attribute dict
    vardict = OrderedDict("units" => "",
                   "long_name" => "",
                   "mode" => ""
                  )
    return vardict
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    add_financial_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")

Adds the Financial results to the dictionary passed back from `run_reopt` using the solved model `m` and the `REoptInputs` for node `_n`.
Note: the node number is an empty string if evaluating a single `Site`.

Financial results:
- `lcc` Optimal lifecycle cost
- `lifecycle_capital_costs_plus_om` Capital cost for all technologies plus present value of operations and maintenance over anlaysis period
- `lifecycle_capital_costs` Net capital costs for all technologies, in present value, including replacement costs and incentives.
"""
function add_financial_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
    r = Dict{String, Any}()
    r["lcc"] = value(m[Symbol("Costs"*_n)]) + 0.0001 * value(m[Symbol("MinChargeAdder"*_n)])
    r["lifecycle_om_costs_before_tax"] = value(m[Symbol("TotalPerUnitSizeOMCosts"*_n)] +
                                           m[Symbol("TotalPerUnitProdOMCosts"*_n)])
    r["year_one_om_costs_before_tax"] = r["lifecycle_om_costs_before_tax"] / (p.pwf_om * p.third_party_factor)
    r["lifecycle_om_costs_after_tax"] = r["lifecycle_om_costs_before_tax"] * (1 - p.s.financial.owner_tax_pct)
    r["year_one_om_costs_after_tax"] = r["lifecycle_om_costs_after_tax"] / (p.pwf_om * p.third_party_factor)

    r["lifecycle_capital_costs_plus_om"] = value(m[Symbol("TotalTechCapCosts"*_n)] + m[Symbol("TotalStorageCapCosts"*_n)]) +
        r["lifecycle_om_costs_after_tax"]
    r["lifecycle_capital_costs"] = value(m[Symbol("TotalTechCapCosts"*_n)] + m[Symbol("TotalStorageCapCosts"*_n)])
    r["initial_capital_costs"] = initial_capex(m, p; _n=_n)
    r["initial_capital_costs_after_incentives"] = initial_capex_after_incentives(m, p, r["lifecycle_capital_costs"]; _n=_n)

    future_replacement_cost, present_replacement_cost = replacement_costs_future_and_present(m, p; _n=_n)
    r["replacement_costs"] = future_replacement_cost
    r["om_and_replacement_present_cost_after_tax"] = present_replacement_cost + r["lifecycle_om_costs_after_tax"]
    r["developer_om_and_replacement_present_cost_after_tax"] = r["om_and_replacement_present_cost_after_tax"] /
        p.third_party_factor

    if !isempty(p.techs.fuel_burning)
        r["lifecycle_fuel_costs_after_tax"] = value(m[:TotalFuelCosts]) * (1 - p.s.financial.offtaker_tax_pct)
    end

    d["Financial"] = Dict(k => round(v, digits=2) for (k,v) in r)
    nothing
end

"""
    initial_capex(m::JuMP.AbstractModel, p::REoptInputs; _n="")

Calculate and return the up-front capital costs for all technologies, in present value, excluding replacement costs and
incentives.
"""
function initial_capex(m::JuMP.AbstractModel, p::REoptInputs; _n="")
    initial_capex = 0

    if !isempty(p.techs.gen) && isempty(_n)  # generators not included in multinode model
        initial_capex += p.s.generator.installed_cost_per_kw * value.(m[Symbol("dvPurchaseSize"*_n)])["Generator"]
    end

    if !isempty(p.techs.pv)
        for pv in p.s.pvs
            initial_capex += pv.installed_cost_per_kw * value.(m[Symbol("dvPurchaseSize"*_n)])[pv.name]
        end
    end

    for b in p.s.storage.types
        if p.s.storage.max_kw[b] > 0
            initial_capex += p.s.storage.raw_inputs[b].installed_cost_per_kw * value.(m[Symbol("dvStoragePower"*_n)])[b] +
                p.s.storage.raw_inputs[b].installed_cost_per_kwh * value.(m[Symbol("dvStorageEnergy"*_n)])[b]
        end
    end

    if "Wind" in p.techs.all
        initial_capex += p.s.wind.installed_cost_per_kw * value.(m[Symbol("dvPurchaseSize"*_n)])["Wind"]
    end

    if "CHP" in p.techs.all
        # CHP.installed_cost_per_kw is now a list with potentially > 1 elements
        cost_list = p.s.chp.installed_cost_per_kw
        size_list = p.s.chp.tech_sizes_for_cost_curve
        chp_size = value.(m[Symbol("dvPurchaseSize"*_n)])["CHP"]
        if typeof(cost_list) == Vector{Float64}
            if chp_size <= size_list[1]
                initial_capex += chp_size * cost_list[1]  # Currently not handling non-zero cost ($) for 0 kW size input
            elseif chp_size > size_list[end]
                initial_capex += chp_size * cost_list[end]
            else
                for s in 2:length(size_list)
                    if (chp_size > size_list[s-1]) && (chp_size <= size_list[s])
                        slope = (cost_list[s] * size_list[s] - cost_list[s-1] * size_list[s-1]) /
                                (size_list[s] - size_list[s-1])
                        initial_capex += cost_list[s-1] * size_list[s-1] + (chp_size - size_list[s-1]) * slope
                    end
                end
            end
        else
            initial_capex += cost_list * chp_size
        #Add supplementary firing capital cost
        # chp_supp_firing_size = self.nested_outputs["Scenario"]["Site"][tech].get("size_supplementary_firing_kw")
        # chp_supp_firing_cost = self.inputs[tech].get("supplementary_firing_capital_cost_per_kw") or 0
        # initial_capex += chp_supp_firing_size * chp_supp_firing_cost
        end
    end

    # TODO thermal tech costs

    return initial_capex
end

"""
    initial_capex_after_incentives(m::JuMP.AbstractModel, p::REoptInputs, lifecycle_capital_costs::Float64; _n="")

The lifecycle_capital_costs output is the initial capex after incentives, except it includes the battery
replacement cost in present value. So we calculate the initial_capex_after_incentives as lifecycle_capital_costs
minus the battery replacement cost in present value.
Note that the owner_discount_pct and owner_tax_pct are set to the offtaker_discount_pct and offtaker_tax_pct
respectively when third_party_ownership is False.
"""
function initial_capex_after_incentives(m::JuMP.AbstractModel, p::REoptInputs, lifecycle_capital_costs::Float64; _n="")
    initial_capex_after_incentives = lifecycle_capital_costs / p.third_party_factor

    for b in p.s.storage.types

        if !(:inverter_replacement_year in fieldnames(typeof(p.s.storage.raw_inputs[b])))
            continue
        end

        pwf_inverter = 1 / ((1 + p.s.financial.owner_discount_pct)^p.s.storage.raw_inputs[b].inverter_replacement_year)

        pwf_storage  = 1 / ((1 + p.s.financial.owner_discount_pct)^p.s.storage.raw_inputs[b].battery_replacement_year)

        inverter_future_cost = p.s.storage.raw_inputs[b].replace_cost_per_kw * value.(m[Symbol("dvStoragePower"*_n)])[b]

        storage_future_cost = p.s.storage.raw_inputs[b].replace_cost_per_kwh * value.(m[Symbol("dvStorageEnergy"*_n)])[b]

        # NOTE these initial costs include the tax benefit available to commercial entities
        initial_capex_after_incentives -= inverter_future_cost * pwf_inverter * (1 - p.s.financial.owner_tax_pct)
        initial_capex_after_incentives -= storage_future_cost  * pwf_storage  * (1 - p.s.financial.owner_tax_pct)
    end

    return initial_capex_after_incentives
end

"""
    replacement_costs_future_and_present(m::JuMP.AbstractModel, p::REoptInputs; _n="")

returns two values: the future and present costs of replacing all storage systems
"""
function replacement_costs_future_and_present(m::JuMP.AbstractModel, p::REoptInputs; _n="")
    future_cost = 0
    present_cost = 0
    for b in p.s.storage.types

        if !(:inverter_replacement_year in fieldnames(typeof(p.s.storage.raw_inputs[b])))
            continue
        end
        future_cost_inverter = p.s.storage.raw_inputs[b].replace_cost_per_kw * value.(m[Symbol("dvStoragePower"*_n)])[b]
        future_cost_storage = p.s.storage.raw_inputs[b].replace_cost_per_kwh * value.(m[Symbol("dvStorageEnergy"*_n)])[b]
        future_cost += future_cost_inverter + future_cost_storage

        present_cost += future_cost_inverter * (1 - p.s.financial.owner_tax_pct) /
            ((1 + p.s.financial.owner_discount_pct)^p.s.storage.raw_inputs[b].inverter_replacement_year)
        present_cost += future_cost_storage * (1 - p.s.financial.owner_tax_pct) /
            ((1 + p.s.financial.owner_discount_pct)^p.s.storage.raw_inputs[b].battery_replacement_year)
    end
    return future_cost, present_cost
end

"""
    calculate_lcoe(p::REoptInputs, tech_results::Dict, tech::AbstractTech)

The Levelized Cost of Energy (LCOE) is calculated as annualized costs (capital and O+M translated to current value)
divided by annuale energy output.
"""
function calculate_lcoe(p::REoptInputs, tech_results::Dict, tech::AbstractTech)
    existing_kw = :existing_kw in fieldnames(typeof(tech)) ? tech.existing_kw : 0.0
    new_kw = get(tech_results, "size_kw", 0) - existing_kw # new capacity
    if new_kw == 0
        return 0.0
    end

    years = p.s.financial.analysis_years # length of financial life
    # TODO is most of this calculated in proforma metrics?
    if p.s.financial.third_party_ownership
        discount_pct = p.s.financial.owner_discount_pct
        federal_tax_pct = p.s.financial.owner_tax_pct
    else
        discount_pct = p.s.financial.offtaker_discount_pct
        federal_tax_pct = p.s.financial.offtaker_tax_pct
    end

    capital_costs = new_kw * tech.installed_cost_per_kw # pre-incentive capital costs

    annual_om = new_kw * tech.om_cost_per_kw # NPV of O&M charges escalated over financial life

    om_series = [annual_om * (1+p.s.financial.om_cost_escalation_pct)^yr for yr in 1:years]
    npv_om = sum([om * (1.0/(1.0+discount_pct))^yr for (yr, om) in enumerate(om_series)])

    #Incentives as calculated in the spreadsheet, note utility incentives are applied before state incentives
    utility_ibi = min(capital_costs * tech.utility_ibi_pct, tech.utility_ibi_max)
    utility_cbi = min(new_kw * tech.utility_rebate_per_kw, tech.utility_rebate_max)
    state_ibi = min((capital_costs - utility_ibi - utility_cbi) * tech.state_ibi_pct, tech.state_ibi_max)
    state_cbi = min(new_kw * tech.state_rebate_per_kw, tech.state_rebate_max)
    federal_cbi = new_kw * tech.federal_rebate_per_kw
    ibi = utility_ibi + state_ibi  #total investment-based incentives
    cbi = utility_cbi + federal_cbi + state_cbi #total capacity-based incentives

    #calculate energy in the BAU case, used twice later on
    existing_energy_bau = get(tech_results, "year_one_energy_produced_kwh_bau", 0)

    #calculate the value of the production-based incentive stream
    npv_pbi = 0
    year_one_energy_produced = get(tech_results, "year_one_energy_produced_kwh", 0)
    degradation_pct = :degradation_pct in fieldnames(typeof(tech)) ? tech.degradation_pct : 0.0
    if tech.production_incentive_max_benefit > 0
        for yr in 1:years
            if yr < tech.production_incentive_years
                degredation_pct = (1- degradation_pct)^yr
                base_pbi = minimum([tech.production_incentive_per_kwh *
                    (year_one_energy_produced - existing_energy_bau) * degredation_pct,
                    tech.production_incentive_max_benefit * degredation_pct
                ])
                npv_pbi += base_pbi * (1.0/(1.0+discount_pct))^(yr+1)
            end
        end
    end

    npv_federal_itc = 0
    federal_itc_basis = capital_costs - state_ibi - utility_ibi - state_cbi - utility_cbi - federal_cbi
    federal_itc_amount = tech.federal_itc_pct * federal_itc_basis
    npv_federal_itc = federal_itc_amount * (1.0/(1.0+discount_pct))

    depreciation_schedule = zeros(years)
    if tech.macrs_option_years in [5,7]
        if tech.macrs_option_years == 5
            schedule = p.s.financial.macrs_five_year
        elseif tech.macrs_option_years == 7
            schedule = p.s.financial.macrs_seven_year
        end
        macrs_bonus_basis = federal_itc_basis - (federal_itc_basis * tech.federal_itc_pct * tech.macrs_itc_reduction)
        macrs_basis = macrs_bonus_basis * (1 - tech.macrs_bonus_pct)
        for (i,r) in enumerate(schedule)
            if i-1 < length(depreciation_schedule)
                depreciation_schedule[i] = macrs_basis * r
            end
        end
        depreciation_schedule[1] += tech.macrs_bonus_pct * macrs_bonus_basis
    end

    tax_deductions = (om_series + depreciation_schedule) * federal_tax_pct
    npv_tax_deductions = sum([i* (1.0/(1.0+discount_pct))^yr for (yr,i) in enumerate(tax_deductions)])

    #we only care about the energy produced by new capacity in LCOE calcs
    annual_energy = year_one_energy_produced - existing_energy_bau
    npv_annual_energy = sum([annual_energy * (1.0/(1.0+discount_pct))^yr *
        (1- degradation_pct)^(yr-1) for yr in 1:years])

    #LCOE is calculated as annualized costs divided by annualized energy
    lcoe = (capital_costs + npv_om - npv_pbi - cbi - ibi - npv_federal_itc - npv_tax_deductions ) / npv_annual_energy

    return round(lcoe, digits=4)
end


text:
# This is code for computing a transit model integrated over
# a time step, giving the fluence in units of time (since
# flux is normalized to unity).

include("transit_poly_struct.jl")

# First the version without derivatives, which tracks the number of depths:
function integrate_timestep_track(param::Array{T,1},trans::Transit_Struct{T},t1::T,t2::T,tol::T,maxdepth::Int64) where {T <: Real}

  function solver(param::Array{T,1},time::T) where {T <: Real}
  # This predicts the impact parameter as a function of time in
  # terms of the transit_struct.
  # For now, transit is approximated as a straight line, where params are given as:
  # param[1] = t0 (mid-transit time; units of days/JD)
  # param[2] = v/R_* (sky velocity in terms of radius of the star - units are day^{-1}
  # param[3] = b (impact parameter at time t0)
  return sqrt(param[3]^2+(param[2]*(time-param[1]))^2)
  end

  # Now, adopt Dan Foreman-Mackey's recursive integration approach:
  function integrate(param::Array{T,1}, f1::T, f2::T, t1::T, t2::T, depth::Int64) where {T <: Real}
  tmid = 0.5*(t1+t2);
  bmid = solver(param,tmid)
  trans.b = bmid
  fmid = transit_poly_d!(trans)
  fapprox = 0.5*(f1+f2)
  d = abs(fmid-fapprox)
  if d > tol && depth < maxdepth
     a,dmaxa,da = integrate(param,f1,fmid,t1,tmid,depth+1)
     b,dmaxb,db = integrate(param,fmid,f2,tmid,t2,depth+1)
     return a+b,maximum([dmaxa,dmaxb]),maximum([da,db])
  end
  return 0.5*(fapprox+fmid)*(t2-t1),depth,d
  end

  b1 = solver(param,t1)
  trans.b = b1
  f1 = transit_poly_d!(trans)
  b2 = solver(param,t2)
  trans.b = b2
  f2 = transit_poly_d!(trans)
  fint,dmax,diff = integrate(param,f1,f2,t1,t2,0)
return fint,dmax,diff
end

# First the version without derivatives, which doesn't track depth:
function integrate_timestep(param::Array{T,1},trans::Transit_Struct{T},t1::T,t2::T,tol::T,maxdepth::Int64) where {T <: Real}

  function solver(param::Array{T,1},time::T) where {T <: Real}
  # This predicts the impact parameter as a function of time in
  # terms of the transit_struct.
  # For now, transit is approximated as a straight line, where params are given as:
  # param[1] = t0 (mid-transit time; units of days/JD)
  # param[2] = v/R_* (sky velocity in terms of radius of the star - units are day^{-1}
  # param[3] = b (impact parameter at time t0)
  return sqrt(param[3]^2+(param[2]*(time-param[1]))^2)
  end

  # Now, adopt Dan Foreman-Mackey's recursive integration approach:
  function integrate(param::Array{T,1}, f1::T, f2::T, t1::T, t2::T, depth::Int64) where {T <: Real}
  tmid = 0.5*(t1+t2);
  bmid = solver(param,tmid)
  trans.b = bmid
  fmid = transit_poly_d!(trans)
  fapprox = 0.5*(f1+f2)
  d = abs(fmid-fapprox)
  if d > tol && depth < maxdepth
     a = integrate(param,f1,fmid,t1,tmid,depth+1)
     b = integrate(param,fmid,f2,tmid,t2,depth+1)
     return a+b
  end
  return 0.5*(fapprox+fmid)*(t2-t1)
  end

  b1 = solver(param,t1)
  trans.b = b1
  f1 = transit_poly_d!(trans)
  b2 = solver(param,t2)
  trans.b = b2
  f2 = transit_poly_d!(trans)
  fint = integrate(param,f1,f2,t1,t2,0)
return fint
end

# Now the version with derivatives:
function integrate_timestep_gradient(param::Array{T,1},trans::Transit_Struct{T},t1::T,t2::T,tol::T,maxdepth::Int64) where {T <: Real}

  function fill_flux!(tmid::T,fmid0::T,trans_mid::Transit_Struct{T}) where {T <: Real}
  fmid = Array{T}(6+trans_mid.n)
  fmid[1] = fmid0   # flux
  fmid[2] = trans_mid.dfdrb[1]  # r derivative
  fmid[3] = trans_mid.dfdrb[2]/trans_mid.b*param[2]^2*(param[1]-tmid)  # t0 derivative
  fmid[4] = trans_mid.dfdrb[2]*param[2]/trans_mid.b*(tmid-param[1])^2  # v derivative
  fmid[5] = trans_mid.dfdrb[2]*param[3]/trans_mid.b  # b0 derivative
  @inbounds for i=1:trans_mid.n+1
    fmid[5+i] = trans_mid.dfdg[i]  # d_i derivatives
  end
  return fmid
  end

  function solver(time::T) where {T <: Real}
  # This predicts the impact parameter as a function of time in
  # terms of the transit_struct.
  # For now, transit is approximated as a straight line, where params are given as:
  # param[1] = t0 (mid-transit time; units of days/JD)
  # param[2] = v/R_* (sky velocity in terms of radius of the star - units are day^{-1}
  # param[3] = b (impact parameter at time t0)
  return sqrt(param[3]^2+(param[2]*(time-param[1]))^2)
  end

  # Now, adopt Dan Foreman-Mackey's recursive integration approach:
  function integrate!(f1::Array{T,1}, f2::Array{T,1}, t1::T, t2::T, depth::Int64) where {T <: Real}
  tmid = 0.5*(t1+t2)
  trans.b = solver(tmid)
  fmid0 = transit_poly_d!(trans)
  fmid = fill_flux!(tmid,fmid0,trans)
#  fmid = [fmid0;trans.dfdrb[1];trans.dfdrb[2]/trans.b*param[2]^2*(param[1]-tmid);
#     trans.dfdrb[2]*param[2]/trans.b*(tmid-param[1])^2;trans.dfdrb[2]*param[3]/trans.b;trans.dfdg]
  # Extrapolate the flux at the boundaries to the midpoint:
  # Compute: flux, df/dr, df/dt0, df/dv, df/db0, df/d d_n:
  fapprox = Array{T}(6+trans.n)
  @inbounds for i=1:6+trans.n
    fapprox[i] = 0.5*(f1[i]+f2[i])
  end
#  println("depth: ",depth," b: ",trans.b," f1: ",f1," fmid: ",fmid," f2: ",f2)
  # Tolerance is set by the flux term only:
#  d = abs(fmid0-fapprox[1])
  d = maximum(abs,fmid-fapprox)
  if d > tol && depth < maxdepth
     # Compute: flux, df/dr, df/dt0, df/dv, df/db0, df/d d_n:
    a = integrate!(f1,fmid,t1,tmid,depth+1)
    b = integrate!(fmid,f2,tmid,t2,depth+1)
    return a+b
  end
  return 0.5*(fapprox+fmid)*(t2-t1)
#  return fapprox
  end

  trans.b = solver(t1)
  f10 = transit_poly_d!(trans)
  f1 = fill_flux!(t1,f10,trans)
  trans.b = solver(t2)
  f20 = transit_poly_d!(trans)
  f2 = fill_flux!(t2,f20,trans)
#  f2 = [f20;trans.dfdrb[1];trans.dfdrb[2]/trans.b*param[2]^2*(param[1]-t2);
#     trans.dfdrb[2]*param[2]/trans.b*(t2-param[1])^2;trans.dfdrb[2]*param[3]/trans.b;trans.dfdg]
  fint = integrate!(f1,f2,t1,t2,0)
return fint
end


text:
using Random: AbstractRNG
using Distributions: Distributions, rand, logpdf, ContinuousMultivariateDistribution, MixtureModel, MvNormal

struct Ring{T<:AbstractFloat} <: ContinuousMultivariateDistribution
    mixturemodel
end

function Ring(n_clusters::Int, distance::T, var::T) where {T<:AbstractFloat}
    pi_typed = convert(T, pi)
    cluster_idcs = collect(0:n_clusters-1)
    base_angle = pi_typed * 2 / n_clusters
    angle = (base_angle .* cluster_idcs) .- pi_typed / 2
    mean = [distance * cos.(angle) distance * sin.(angle)]'
    return Ring{T}(MixtureModel([MvNormal(mean[:,i], var) for i in 1:n_clusters]))
end

Distributions.rand(rng::AbstractRNG, ring::Ring{T}, n::Int) where {T<:AbstractFloat} =
    convert.(T, rand(rng, ring.mixturemodel, n))
Distributions.logpdf(dataset::Ring, X::AbstractArray{<:AbstractFloat,2}) = logpdf(dataset.mixturemodel, X)

struct RingDataset{D} <: AbstractDataset{D}
    X
    Xt
end

n_display(d::RingDataset) = 2_000

function RingDataset(
    n_data::Int,
    z_angle::T1=Float32(pi / 3),
    z_std::T1=1f-1,
    n_clusters::Int=10,
    distance::T2=25f-1,
    var::T2=25f-2;
    seed::Int=1,
    test_ratio=1/6,
    n_test::Int=ratio2num(n_data, test_ratio),
    rng=MersenneTwister(seed),
) where {T<:AbstractFloat, T1<:Union{Nothing, T}, T2<:T}
    ring = Ring(n_clusters, distance, var)
    Tz(angle) = [cos(angle) 0 sin(angle); 0 1 0; -sin(angle) 0 cos(angle)]
    make3d(X) = Tz(z_angle) * cat(X, randn(T1, 1, size(X, 2)) * z_std; dims=1)
    X, Xt = rand(rng, ring, n_data), rand(rng, ring, n_test)
    if isnothing(z_angle) || isnothing(z_std)
        @info "Oh you just get the 2D ring dataset" n_data=n_data n_test=n_test n_clusters=n_clusters distance=distance var=var
        return RingDataset{2}(X, Xt)
    else
        X, Xt = make3d.((X, Xt))
        @info "Oh you just get the 3D ring dataset" n_data=n_data n_test=n_test n_clusters=n_clusters distance=distance var=var z_angle=z_angle z_std=z_std
        return RingDataset{3}(X, Xt)
    end
end

RingDataset{2}(n_data::Int=60_000, args...; kwargs...) = RingDataset(n_data, nothing, nothing, args...; kwargs...)
TwoDimRingDataset = RingDataset{2}

RingDataset{3}(n_data::Int=60_000, args...; kwargs...) = RingDataset(n_data, args...; kwargs...)
ThreeDimRingDataset = RingDataset{3}


text:
module MadNLPGPU

import CUDA: CUBLAS, CUSOLVER, CuVector, CuMatrix, toolkit_version, R_64F, has_cuda
import MadNLP:
    @kwdef, Logger, @debug, @warn, @error,
    AbstractOptions, AbstractLinearSolver, set_options!, MadNLPLapackCPU,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    introduce, factorize!, solve!, improve!, is_inertia, inertia, tril_to_full!

if has_cuda()
    include("lapackgpu.jl")
    export MadNLPLapackGPU
end

end # module


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
# import ._ProtoBuf_Top_.google

struct __enum_EncodingType <: ProtoEnum
    NONE::Int32
    UTF8::Int32
    UTF16::Int32
    UTF32::Int32
    __enum_EncodingType() = new(0,1,2,3)
end #struct __enum_EncodingType
const EncodingType = __enum_EncodingType()

struct __enum_Document_Type <: ProtoEnum
    TYPE_UNSPECIFIED::Int32
    PLAIN_TEXT::Int32
    HTML::Int32
    __enum_Document_Type() = new(0,1,2)
end #struct __enum_Document_Type
const Document_Type = __enum_Document_Type()

mutable struct Document <: ProtoType
    _type::Int32
    content::AbstractString
    gcs_content_uri::AbstractString
    language::AbstractString
    Document(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Document
const __oneofs_Document = Int[0,1,1,0]
const __oneof_names_Document = [Symbol("source")]
meta(t::Type{Document}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, __oneofs_Document, __oneof_names_Document, ProtoBuf.DEF_FIELD_TYPES)

mutable struct Sentiment <: ProtoType
    magnitude::Float32
    score::Float32
    Sentiment(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Sentiment
const __fnum_Sentiment = Int[2,3]
meta(t::Type{Sentiment}) = meta(t, ProtoBuf.DEF_REQ, __fnum_Sentiment, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

struct __enum_PartOfSpeech_Tag <: ProtoEnum
    UNKNOWN::Int32
    ADJ::Int32
    ADP::Int32
    ADV::Int32
    CONJ::Int32
    DET::Int32
    NOUN::Int32
    NUM::Int32
    PRON::Int32
    PRT::Int32
    PUNCT::Int32
    VERB::Int32
    X::Int32
    AFFIX::Int32
    __enum_PartOfSpeech_Tag() = new(0,1,2,3,4,5,6,7,8,9,10,11,12,13)
end #struct __enum_PartOfSpeech_Tag
const PartOfSpeech_Tag = __enum_PartOfSpeech_Tag()

struct __enum_PartOfSpeech_Aspect <: ProtoEnum
    ASPECT_UNKNOWN::Int32
    PERFECTIVE::Int32
    IMPERFECTIVE::Int32
    PROGRESSIVE::Int32
    __enum_PartOfSpeech_Aspect() = new(0,1,2,3)
end #struct __enum_PartOfSpeech_Aspect
const PartOfSpeech_Aspect = __enum_PartOfSpeech_Aspect()

struct __enum_PartOfSpeech_Case <: ProtoEnum
    CASE_UNKNOWN::Int32
    ACCUSATIVE::Int32
    ADVERBIAL::Int32
    COMPLEMENTIVE::Int32
    DATIVE::Int32
    GENITIVE::Int32
    INSTRUMENTAL::Int32
    LOCATIVE::Int32
    NOMINATIVE::Int32
    OBLIQUE::Int32
    PARTITIVE::Int32
    PREPOSITIONAL::Int32
    REFLEXIVE_CASE::Int32
    RELATIVE_CASE::Int32
    VOCATIVE::Int32
    __enum_PartOfSpeech_Case() = new(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14)
end #struct __enum_PartOfSpeech_Case
const PartOfSpeech_Case = __enum_PartOfSpeech_Case()

struct __enum_PartOfSpeech_Form <: ProtoEnum
    FORM_UNKNOWN::Int32
    ADNOMIAL::Int32
    AUXILIARY::Int32
    COMPLEMENTIZER::Int32
    FINAL_ENDING::Int32
    GERUND::Int32
    REALIS::Int32
    IRREALIS::Int32
    SHORT::Int32
    LONG::Int32
    ORDER::Int32
    SPECIFIC::Int32
    __enum_PartOfSpeech_Form() = new(0,1,2,3,4,5,6,7,8,9,10,11)
end #struct __enum_PartOfSpeech_Form
const PartOfSpeech_Form = __enum_PartOfSpeech_Form()

struct __enum_PartOfSpeech_Gender <: ProtoEnum
    GENDER_UNKNOWN::Int32
    FEMININE::Int32
    MASCULINE::Int32
    NEUTER::Int32
    __enum_PartOfSpeech_Gender() = new(0,1,2,3)
end #struct __enum_PartOfSpeech_Gender
const PartOfSpeech_Gender = __enum_PartOfSpeech_Gender()

struct __enum_PartOfSpeech_Mood <: ProtoEnum
    MOOD_UNKNOWN::Int32
    CONDITIONAL_MOOD::Int32
    IMPERATIVE::Int32
    INDICATIVE::Int32
    INTERROGATIVE::Int32
    JUSSIVE::Int32
    SUBJUNCTIVE::Int32
    __enum_PartOfSpeech_Mood() = new(0,1,2,3,4,5,6)
end #struct __enum_PartOfSpeech_Mood
const PartOfSpeech_Mood = __enum_PartOfSpeech_Mood()

struct __enum_PartOfSpeech_Number <: ProtoEnum
    NUMBER_UNKNOWN::Int32
    SINGULAR::Int32
    PLURAL::Int32
    DUAL::Int32
    __enum_PartOfSpeech_Number() = new(0,1,2,3)
end #struct __enum_PartOfSpeech_Number
const PartOfSpeech_Number = __enum_PartOfSpeech_Number()

struct __enum_PartOfSpeech_Person <: ProtoEnum
    PERSON_UNKNOWN::Int32
    FIRST::Int32
    SECOND::Int32
    THIRD::Int32
    REFLEXIVE_PERSON::Int32
    __enum_PartOfSpeech_Person() = new(0,1,2,3,4)
end #struct __enum_PartOfSpeech_Person
const PartOfSpeech_Person = __enum_PartOfSpeech_Person()

struct __enum_PartOfSpeech_Proper <: ProtoEnum
    PROPER_UNKNOWN::Int32
    PROPER::Int32
    NOT_PROPER::Int32
    __enum_PartOfSpeech_Proper() = new(0,1,2)
end #struct __enum_PartOfSpeech_Proper
const PartOfSpeech_Proper = __enum_PartOfSpeech_Proper()

struct __enum_PartOfSpeech_Reciprocity <: ProtoEnum
    RECIPROCITY_UNKNOWN::Int32
    RECIPROCAL::Int32
    NON_RECIPROCAL::Int32
    __enum_PartOfSpeech_Reciprocity() = new(0,1,2)
end #struct __enum_PartOfSpeech_Reciprocity
const PartOfSpeech_Reciprocity = __enum_PartOfSpeech_Reciprocity()

struct __enum_PartOfSpeech_Tense <: ProtoEnum
    TENSE_UNKNOWN::Int32
    CONDITIONAL_TENSE::Int32
    FUTURE::Int32
    PAST::Int32
    PRESENT::Int32
    IMPERFECT::Int32
    PLUPERFECT::Int32
    __enum_PartOfSpeech_Tense() = new(0,1,2,3,4,5,6)
end #struct __enum_PartOfSpeech_Tense
const PartOfSpeech_Tense = __enum_PartOfSpeech_Tense()

struct __enum_PartOfSpeech_Voice <: ProtoEnum
    VOICE_UNKNOWN::Int32
    ACTIVE::Int32
    CAUSATIVE::Int32
    PASSIVE::Int32
    __enum_PartOfSpeech_Voice() = new(0,1,2,3)
end #struct __enum_PartOfSpeech_Voice
const PartOfSpeech_Voice = __enum_PartOfSpeech_Voice()

mutable struct PartOfSpeech <: ProtoType
    tag::Int32
    aspect::Int32
    case::Int32
    form::Int32
    gender::Int32
    mood::Int32
    number::Int32
    person::Int32
    proper::Int32
    reciprocity::Int32
    tense::Int32
    voice::Int32
    PartOfSpeech(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct PartOfSpeech

struct __enum_DependencyEdge_Label <: ProtoEnum
    UNKNOWN::Int32
    ABBREV::Int32
    ACOMP::Int32
    ADVCL::Int32
    ADVMOD::Int32
    AMOD::Int32
    APPOS::Int32
    ATTR::Int32
    AUX::Int32
    AUXPASS::Int32
    CC::Int32
    CCOMP::Int32
    CONJ::Int32
    CSUBJ::Int32
    CSUBJPASS::Int32
    DEP::Int32
    DET::Int32
    DISCOURSE::Int32
    DOBJ::Int32
    EXPL::Int32
    GOESWITH::Int32
    IOBJ::Int32
    MARK::Int32
    MWE::Int32
    MWV::Int32
    NEG::Int32
    NN::Int32
    NPADVMOD::Int32
    NSUBJ::Int32
    NSUBJPASS::Int32
    NUM::Int32
    NUMBER::Int32
    P::Int32
    PARATAXIS::Int32
    PARTMOD::Int32
    PCOMP::Int32
    POBJ::Int32
    POSS::Int32
    POSTNEG::Int32
    PRECOMP::Int32
    PRECONJ::Int32
    PREDET::Int32
    PREF::Int32
    PREP::Int32
    PRONL::Int32
    PRT::Int32
    PS::Int32
    QUANTMOD::Int32
    RCMOD::Int32
    RCMODREL::Int32
    RDROP::Int32
    REF::Int32
    REMNANT::Int32
    REPARANDUM::Int32
    ROOT::Int32
    SNUM::Int32
    SUFF::Int32
    TMOD::Int32
    TOPIC::Int32
    VMOD::Int32
    VOCATIVE::Int32
    XCOMP::Int32
    SUFFIX::Int32
    TITLE::Int32
    ADVPHMOD::Int32
    AUXCAUS::Int32
    AUXVV::Int32
    DTMOD::Int32
    FOREIGN::Int32
    KW::Int32
    LIST::Int32
    NOMC::Int32
    NOMCSUBJ::Int32
    NOMCSUBJPASS::Int32
    NUMC::Int32
    COP::Int32
    DISLOCATED::Int32
    ASP::Int32
    GMOD::Int32
    GOBJ::Int32
    INFMOD::Int32
    MES::Int32
    NCOMP::Int32
    __enum_DependencyEdge_Label() = new(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82)
end #struct __enum_DependencyEdge_Label
const DependencyEdge_Label = __enum_DependencyEdge_Label()

mutable struct DependencyEdge <: ProtoType
    head_token_index::Int32
    label::Int32
    DependencyEdge(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DependencyEdge

mutable struct TextSpan <: ProtoType
    content::AbstractString
    begin_offset::Int32
    TextSpan(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TextSpan

mutable struct Token <: ProtoType
    text::TextSpan
    part_of_speech::PartOfSpeech
    dependency_edge::DependencyEdge
    lemma::AbstractString
    Token(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Token

mutable struct Sentence <: ProtoType
    text::TextSpan
    sentiment::Sentiment
    Sentence(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Sentence

struct __enum_EntityMention_Type <: ProtoEnum
    TYPE_UNKNOWN::Int32
    PROPER::Int32
    COMMON::Int32
    __enum_EntityMention_Type() = new(0,1,2)
end #struct __enum_EntityMention_Type
const EntityMention_Type = __enum_EntityMention_Type()

mutable struct EntityMention <: ProtoType
    text::TextSpan
    _type::Int32
    sentiment::Sentiment
    EntityMention(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EntityMention

struct __enum_Entity_Type <: ProtoEnum
    UNKNOWN::Int32
    PERSON::Int32
    LOCATION::Int32
    ORGANIZATION::Int32
    EVENT::Int32
    WORK_OF_ART::Int32
    CONSUMER_GOOD::Int32
    OTHER::Int32
    PHONE_NUMBER::Int32
    ADDRESS::Int32
    DATE::Int32
    NUMBER::Int32
    PRICE::Int32
    __enum_Entity_Type() = new(0,1,2,3,4,5,6,7,9,10,11,12,13)
end #struct __enum_Entity_Type
const Entity_Type = __enum_Entity_Type()

mutable struct Entity_MetadataEntry <: ProtoType
    key::AbstractString
    value::AbstractString
    Entity_MetadataEntry(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Entity_MetadataEntry (mapentry)

mutable struct Entity <: ProtoType
    name::AbstractString
    _type::Int32
    metadata::Base.Dict{AbstractString,AbstractString} # map entry
    salience::Float32
    mentions::Base.Vector{EntityMention}
    sentiment::Sentiment
    Entity(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Entity

mutable struct ClassificationCategory <: ProtoType
    name::AbstractString
    confidence::Float32
    ClassificationCategory(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ClassificationCategory

mutable struct AnalyzeSentimentRequest <: ProtoType
    document::Document
    encoding_type::Int32
    AnalyzeSentimentRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AnalyzeSentimentRequest

mutable struct AnalyzeSentimentResponse <: ProtoType
    document_sentiment::Sentiment
    language::AbstractString
    sentences::Base.Vector{Sentence}
    AnalyzeSentimentResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AnalyzeSentimentResponse

mutable struct AnalyzeEntitySentimentRequest <: ProtoType
    document::Document
    encoding_type::Int32
    AnalyzeEntitySentimentRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AnalyzeEntitySentimentRequest

mutable struct AnalyzeEntitySentimentResponse <: ProtoType
    entities::Base.Vector{Entity}
    language::AbstractString
    AnalyzeEntitySentimentResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AnalyzeEntitySentimentResponse

mutable struct AnalyzeEntitiesRequest <: ProtoType
    document::Document
    encoding_type::Int32
    AnalyzeEntitiesRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AnalyzeEntitiesRequest

mutable struct AnalyzeEntitiesResponse <: ProtoType
    entities::Base.Vector{Entity}
    language::AbstractString
    AnalyzeEntitiesResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AnalyzeEntitiesResponse

mutable struct AnalyzeSyntaxRequest <: ProtoType
    document::Document
    encoding_type::Int32
    AnalyzeSyntaxRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AnalyzeSyntaxRequest

mutable struct AnalyzeSyntaxResponse <: ProtoType
    sentences::Base.Vector{Sentence}
    tokens::Base.Vector{Token}
    language::AbstractString
    AnalyzeSyntaxResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AnalyzeSyntaxResponse

mutable struct ClassifyTextRequest <: ProtoType
    document::Document
    ClassifyTextRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ClassifyTextRequest

mutable struct ClassifyTextResponse <: ProtoType
    categories::Base.Vector{ClassificationCategory}
    ClassifyTextResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ClassifyTextResponse

mutable struct AnnotateTextRequest_Features <: ProtoType
    extract_syntax::Bool
    extract_entities::Bool
    extract_document_sentiment::Bool
    extract_entity_sentiment::Bool
    classify_text::Bool
    AnnotateTextRequest_Features(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AnnotateTextRequest_Features
const __fnum_AnnotateTextRequest_Features = Int[1,2,3,4,6]
meta(t::Type{AnnotateTextRequest_Features}) = meta(t, ProtoBuf.DEF_REQ, __fnum_AnnotateTextRequest_Features, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct AnnotateTextRequest <: ProtoType
    document::Document
    features::AnnotateTextRequest_Features
    encoding_type::Int32
    AnnotateTextRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AnnotateTextRequest

mutable struct AnnotateTextResponse <: ProtoType
    sentences::Base.Vector{Sentence}
    tokens::Base.Vector{Token}
    entities::Base.Vector{Entity}
    document_sentiment::Sentiment
    language::AbstractString
    categories::Base.Vector{ClassificationCategory}
    AnnotateTextResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AnnotateTextResponse

export EncodingType, Document_Type, Document, Sentence, Entity_Type, Entity_MetadataEntry, Entity, Token, Sentiment, PartOfSpeech_Tag, PartOfSpeech_Aspect, PartOfSpeech_Case, PartOfSpeech_Form, PartOfSpeech_Gender, PartOfSpeech_Mood, PartOfSpeech_Number, PartOfSpeech_Person, PartOfSpeech_Proper, PartOfSpeech_Reciprocity, PartOfSpeech_Tense, PartOfSpeech_Voice, PartOfSpeech, DependencyEdge_Label, DependencyEdge, EntityMention_Type, EntityMention, TextSpan, ClassificationCategory, AnalyzeSentimentRequest, AnalyzeSentimentResponse, AnalyzeEntitySentimentRequest, AnalyzeEntitySentimentResponse, AnalyzeEntitiesRequest, AnalyzeEntitiesResponse, AnalyzeSyntaxRequest, AnalyzeSyntaxResponse, ClassifyTextRequest, ClassifyTextResponse, AnnotateTextRequest_Features, AnnotateTextRequest, AnnotateTextResponse
# mapentries: "Entity_MetadataEntry" => ("AbstractString", "AbstractString")

max_stars_repo_name: sivakon/grpc-julia-gcp-example

max_issues_repo_name: sivakon/grpc-julia-gcp-example

max_forks_repo_name: sivakon/grpc-julia-gcp-example


text:
@doc """ sam_index_load will call sam_index_load2 and search the string(fn,".bai")
         for the index file
""" ->
function sam_index_load(sios::SamIOStream,fname::AbstractString)
    fp = sios.handle
    fn = convert(Ptr{Cchar}, pointer(fname))
    pidx = sam_index_load(fp,fn)
end

@doc """ sam_index_build will always build .bai index file
""" ->
function sam_index_build(fname::AbstractString)
    fn = convert(Ptr{Cchar}, pointer(fname))
    sam_index_build(fn,0)
end

@doc """ query reads in a region
""" ->
function sam_itr_querys(pidx,phdr::Ptr{Header},region::ASCIIString)
    region = string(region,'\0')
    region = convert(Ptr{Cchar},pointer(region))
    ret = sam_itr_querys(pidx,phdr,region)
    if ret == C_NULL
        info("sam_itr_querys return NULL")
    end
    ret
end

@doc """ produce record in the itr
""" ->
function sam_itr_next!(handle::Ptr{Void}, itr::Ptr{Void}, r::Ptr{Record})
    phtsf = convert(Ptr{HTSFile}, handle)
    htsfile = unsafe_load(phtsf) # may change unsigned bits to signed
#    @show htsfile.lineno
#    @show htsfile.line
#    @show bytestring(htsfile.fn)
    bgzf = htsfile.bgzf
    flag = sam_itr_next(bgzf,itr,r,handle)
    flag >=0 ? true : false
end


text:
using Images, SIUnits, Base.Test

img = load("test.imagine")
@test ndims(img) == 4
@test size(img) == (5,7,3,4)
@test timedim(img) == 4
@test spatialorder(img) == ["x", "l", "z"]
@test pixelspacing(img)[1:2] == [0.71e-6Meter, 0.71e-6Meter]
@test_approx_eq pixelspacing(img)[3].val (2/3)*1e-4

bn = tempname()
ifn = string(bn, ".imagine")
cfn = string(bn, ".cam")
A = rand(Float32,2,3,4,5)
open(cfn, "w") do io
    write(io, A)
end
ImagineFormat.save_header(ifn, "test.imagine", A)
img2 = load(ifn)
@test eltype(img2) == Float32
@test data(img2) == A
rm(ifn)
rm(cfn)

using ImagineFormat
io = IOBuffer()
imagine2nrrd(io, img["imagineheader"])
str = takebuf_string(io)
@test str == "NRRD0001\ntype: uint16\ndimension: 4\nsizes: 5 7 3 4\nkinds: space space space time\nencoding: raw\nendian: little\n"


text:
using OmegaGrav
using Dierckx
using Tables, CSV
# %% Compute the matter power spectrum using CLASS
include("compute_pk_class.jl")
# %% Compute Ωgrav at seven redshifts
redshift = [0, 0.3, 0.5, 0.7, 1, 1.3, 1.5]
nred = length(redshift)
Ωgpk = zeros(nred)
Ωglin = zeros(nred)
Ωghalo = zeros(nred)
Ωgh12 = zeros(nred)
Ωgh13 = zeros(nred)
Ωgh14 = zeros(nred)
# ΩghaloT10 = zeros(nred)
# Ωgh12T10 = zeros(nred)
# Ωgh13T10 = zeros(nred)
# Ωgh14T10 = zeros(nred)
for ired = 1:nred
   z = redshift[ired]
   # Define functions to return linear and non-linear power spectra
   # Note: The CLASS code takes wavenumbers in units of 1/Mpc (no h) and
   # return power spectra in units of Mpc^3 (no 1/h^3).
   pkcb_class(kovh) = cosmo.pk_cb_lin(kovh * h0, z) * h0^3
   pknl_class(kovh) = cosmo.pk(kovh * h0, z) * h0^3
   pklin_class(kovh) = cosmo.pk_lin(kovh * h0, z) * h0^3
   # Spline interpolate in log(k)
   lnk = log(1e-4):0.05:log(100)
   pkcb = Spline1D(lnk, pkcb_class.(exp.(lnk)))
   pknl = Spline1D(lnk, pknl_class.(exp.(lnk)))
   pklin = Spline1D(lnk, pklin_class.(exp.(lnk)))
   # %% Compute Ωgrav from non-linear total matter P(k)
   Ωgpk[ired] = ograv_pk(x -> pknl(log(x)), z, Ωm)
   # %% Compute Ωgrav from linear total matter P(k)
   Ωglin[ired] = ograv_pk(x -> pklin(log(x)), z, Ωm)
   # %% Compute Ωgrav from Halos, excluding the neutrino contribution
   Ωghalo[ired] = ograv_halo(x -> pkcb(log(x)), z, Ωm, Ωcb)
   Ωgh12[ired] = ograv_halo(x -> pkcb(log(x)), z, Ωm, Ωcb, Mmin = 1e12)
   Ωgh13[ired] = ograv_halo(x -> pkcb(log(x)), z, Ωm, Ωcb, Mmin = 1e13)
   Ωgh14[ired] = ograv_halo(x -> pkcb(log(x)), z, Ωm, Ωcb, Mmin = 1e14)
   # %% Compute Ωgrav from Halos, excluding the neutrino contribution, with Tinker10MF
   # ΩghaloT10[ired] = ograv_halo(x -> pkcb(log(x)), z, Ωm, Ωcb, t10MF = true)
   # Ωgh12T10[ired] = ograv_halo(x -> pkcb(log(x)), z, Ωm, Ωcb, Mmin = 1e12, t10MF = true)
   # Ωgh13T10[ired] = ograv_halo(x -> pkcb(log(x)), z, Ωm, Ωcb, Mmin = 1e13, t10MF = true)
   # Ωgh14T10[ired] = ograv_halo(x -> pkcb(log(x)), z, Ωm, Ωcb, Mmin = 1e14, t10MF = true)
end

#%% Save table data to table1.csv
fb = params["omega_b"] / (params["omega_cdm"] + params["omega_b"])
t = Tables.table([redshift Ωglin Ωgpk Ωghalo Ωgh12 Ωgh13 Ωgh14 2 * fb * Ωghalo / 3 2 * fb * (Ωgpk - Ωglin) / 3])
header = [
   "z",
   "Omega_W/2_pklin",
   "Omega_W/2_pknl",
   "Omega_W/2_halo",
   "Omega_W/2_halo_Mgt1e12",
   "Omega_W/2_halo_Mgt1e13",
   "Omega_W/2_halo_Mgt1e14",
   "Omega_W/2_halo_times_twothirdfb",
   "Omega_W/2_pkdiff_times_twothirdfb",
]
CSV.write("table1.csv", t, header = header)

# %% Clean CLASS (the equivalent of the struct_free() in the `main`
# of CLASS. This step is primordial when running in a loop over different
# cosmologies, as you will saturate your memory very fast if you ommit
# it.
cosmo.struct_cleanup()
# If you want to change completely the cosmology, you should also
# clean the arguments, otherwise, if you are simply running on a loop
# of different values for the same parameters, this step is not needed
cosmo.empty()


text:
# Autogenerated wrapper script for eccodes_jll for x86_64-apple-darwin
export eccodes

using libpng_jll
using OpenJpeg_jll
JLLWrappers.@generate_wrapper_header("eccodes")
JLLWrappers.@declare_library_product(eccodes, "@rpath/libeccodes.dylib")
function __init__()
    JLLWrappers.@generate_init_header(libpng_jll, OpenJpeg_jll)
    JLLWrappers.@init_library_product(
        eccodes,
        "lib/libeccodes.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
## Model.jl --- copy model for protein structures

export extractframe
export JAtom, JResidue, JConnectivity, StructureFrame
export getatom, setatom!

"""
    JAtom{T<:AbstractFloat}

Represents an atom in a probed frame. Contains basic information such as its
`name` and `type` as given in the PDB (input file). Mass and charge are
obtained from Chemfiles inference given atom type.

**Note**: charge of an atom is `e * jatom.charge` with `e` as elementary charge.
"""
struct JAtom{T<:AbstractFloat}   ## make this mutable??
    name::Symbol
    type::String
    mass::T
    charge::T
end

"""
    JResidue

A residue in a probed frame. Contains the residue `name`. Holds a list of *indices*
of atoms in the `StructureFrame`. Other properties include the `chainid` specifying
the chain in which the residue belongs and whether the residue is included in the
described protein structure.
"""
struct JResidue
    name::Symbol
    chainid::Union{String,Missing}
    at_dict::Dict{Symbol,UInt64}
    standard_pdb::Union{Bool,Nothing}
end

"""
    JConnectivity

Topology of atoms in the probed frame. Contains bonds,
angles, dihedrals, and improper dihedrals. The four connectivity
types are represented by `NTuple{N,UInt64}` containing indices of atoms in the
probed frame.
"""
struct JConnectivity
    bonds::Vector{NTuple{2,UInt64}}
    angles::Vector{NTuple{3,UInt64}}
    dihedrals::Vector{NTuple{4,UInt64}}
    impropers::Vector{NTuple{4,UInt64}}
end

"""
    StructureFrame

Overall model for protein structure. Contains information on the atomic
coordinates, atom information, and the residues. Is obtained from using
[`extractframe`](@ref) on a Chemfiles Frame object.
"""
struct StructureFrame
    step::UInt64                    ## timestep probed (only one for PDB)
    at_pos::Vector{SVector{3,T}} where T <: AbstractFloat
    at_list::Vector{JAtom}
    res_list::Vector{JResidue}
end

"""
    extractframe(fr::Chemfiles.Frame)

Extract `Chemfiles.Frame` to get objects of structs defined in the package.
Note that not all properties will be inherited/converted.
"""
function extractframe(fr::Chemfiles.Frame; sort=false)
    step_ = step(fr)    # frame step
    sz = size(fr)

    ## get topology and get residue/connectivity
    top = Chemfiles.Topology(fr)

    ## get positions and list of atoms; assert equal size
    pos = Chemfiles.positions(fr)
    pos_ = map(SVector{3}, Iterators.partition(pos, 3))

    ### get residues where atoms are and create from Chemfiles data
    at_list_ = [JAtom(Chemfiles.Atom(fr, i-1)) for i in 1:sz]

    @assert length(pos_) == length(at_list_) "atom and positions list sizes differ"

    # get residues with atom name dictionary
    res_list_ = map(1:Chemfiles.count_residues(top)) do i
        JResidue(Chemfiles.Residue(top, i-1), at_list_)
    end

    if sort
        sort!(res_list_, lt=(x,y)->isless(x.standard_pdb,y.standard_pdb), rev=true)
    end

    return StructureFrame(step_, pos_, at_list_, res_list_)
end

function JAtom(atom::Chemfiles.Atom)
    name_ = Symbol(Chemfiles.name(atom))
    type_ = Chemfiles.type(atom)
    mass_ = Chemfiles.mass(atom)
    charge_ = Chemfiles.charge(atom)
    return JAtom(name_, type_ , mass_, charge_)
end

function JResidue(res::Chemfiles.Residue, reflist)
    plist = Chemfiles.list_properties(res)

    ## atom name dictionary based on reflist (frame atom list)
    at_dict = Dict(reflist[i+1].name=>i+1 for i in Chemfiles.atoms(res))

    name_ = Symbol(Chemfiles.name(res))
    chainid_ = let kw = "chainid"
        ifelse(kw in plist, Chemfiles.property(res, kw), missing)
    end
    pdb_ = let kw = "is_standard_pdb"
        ifelse(kw in plist, Chemfiles.property(res, kw), nothing)
    end
    return JResidue(name_, chainid_, at_dict, pdb_)
end

"""
    getatom(residue, name::Symbol) -> UInt64

Get the _index_ of an atom in the parent frame of the residue that has the
name provided.
"""
getatom(res, s) = res.at_dict[s]

"""
    setatom!(residue, atom::JAtom, idx)

Sets an atom to a residue.

**Note**: This does not check if the atom with a given `.name` is already in
the residue and instead mutates it.
"""
setatom!(res, atom::JAtom, idx) = push!(res.at_dict, atom.name=>idx)

### extension models
include("extensions/FrameTools.jl")
include("extensions/HBond.jl")