text:
## Spectral density ##

### Multivariate ###

"""
    add_sdf!(out, model::TimeSeriesModel, ω)

Add the sdf of a model to a storage vector.
"""
function add_sdf!(out, model::TimeSeriesModel, ω) # default sdf returns error
    error("sdf not yet defined for model of type $(typeof(model)).")
end

"""
    add_asdf!(out, model::TimeSeriesModel, ω, Δ)

Add the aliased sdf of a model to a storage vector.
"""
function add_asdf!(out, model::TimeSeriesModel, ω, Δ) # default method approximates the asdf (overload if known)
    for k ∈ -nalias(model):nalias(model)
        add_sdf!(out, model, ω + k * 2π / Δ)
    end
    return nothing
end

"""
    asdf!(out, model::TimeSeriesModel, ω, Δ)

Computes the asdf and overwrites a storage vector.
"""
function asdf!(out, model::TimeSeriesModel, ω, Δ)
    out .= zero(eltype(out))
    add_asdf!(out, model, ω, Δ)
    return nothing
end

"""
    asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel, freq::FreqAcvEst)

Compute the asdf for all frequencies and allocates to appropriate location in storage.
"""
function asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel, freq::FreqAcvEst)
    size(store.allocatedarray,2) == length(freq.Ωₘ) || throw(ArgumentError("size(store.allocatedarray,2) != length(Ω)."))
    @inbounds for i ∈ 1:length(freq.Ωₘ)
        @views asdf!(store.allocatedarray[:, i], model, freq.Ωₘ[i], freq.Δ)
    end
    return nothing
end
function asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel{D,T}, freq::FreqAcvEst) where {D,T<:Real} ## real values
    size(store.allocatedarray,2) == length(freq.Ωₘ) || throw(ArgumentError("size(store.allocatedarray,2) != length(Ω)."))
    m = length(freq.Ωₘ)
    startcopyind = m÷2+2
    @inbounds for i ∈ 1:startcopyind-1
        @views asdf!(store.allocatedarray[:, i], model, freq.Ωₘ[i], freq.Δ)
    end
    @inbounds for i in startcopyind:length(freq.Ωₘ)
        @views store.allocatedarray[:, i] .= conj.(store.allocatedarray[:, m-i+2])
    end
    return nothing
end
function asdf!(store::TimeSeriesModelStorageFunction, model::TimeSeriesModel)
    asdf!(store.funcmemory, model, store.encodedtime)
    return nothing
end

"""
    sdf(model::TimeSeriesModel, ω)

Compute sdf at `ω`.
"""
function sdf(model::TimeSeriesModel, ω)
    out = zeros(ComplexF64, nlowertriangle_dimension(model))
    add_sdf!(out, model, ω)
    return SHermitianCompact(SVector{nlowertriangle_dimension(model)}(out))
end

"""
    asdf(model::TimeSeriesModel, ω, Δ)

Compute asdf with sampling rate `Δ` at `ω`.
"""
function asdf(model::TimeSeriesModel, ω, Δ)
    out = zeros(ComplexF64, nlowertriangle_dimension(model))
    asdf!(out, model, ω, Δ)
    return SHermitianCompact(SVector{nlowertriangle_dimension(model)}(out))
end

### Univariate ###
function sdf(model::TimeSeriesModel{1,T}, ω) where {T} # default sdf returns error
    error("sdf not yet defined for model of type $(typeof(model)).")
end
function asdf(model::TimeSeriesModel{1,T}, ω, Δ) where {T} # default method approximates the asdf (overload if known)
    val = zero(ComplexF64)
    for k ∈ -nalias(model):nalias(model)
        val += sdf(model, ω + k * 2π / Δ)
    end
    return val
end

function asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel{1,T}, freq::FreqAcvEst) where {T}
    length(store.allocatedarray) == length(freq.Ωₘ) || throw(ArgumentError("length(store.allocatedarray) != length(Ω)."))
    @inbounds for i ∈ 1:length(freq.Ωₘ)
        store.allocatedarray[i] = asdf(model, freq.Ωₘ[i], freq.Δ)
    end
    return nothing
end
function asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel{1,T}, freq::FreqAcvEst) where {T<:Real} ## real values
    length(store.allocatedarray) == length(freq.Ωₘ) || throw(ArgumentError("length(store.allocatedarray) != length(Ω)."))
    m = length(freq.Ωₘ)
    startcopyind = m÷2+2
    @inbounds for i ∈ 1:startcopyind-1
        store.allocatedarray[i] = asdf(model, freq.Ωₘ[i], freq.Δ)
    end
    @inbounds for i ∈ startcopyind:length(freq.Ωₘ)
        store.allocatedarray[i] = store.allocatedarray[m-i+2]
    end
    return nothing
end

### Additive ###
function add_sdf!(out, model::AdditiveTimeSeriesModel, ω)
    add_sdf!(out, model.model1, ω)
    add_sdf!(out, model.model2, ω)
    return nothing
end
function add_asdf!(out, model::AdditiveTimeSeriesModel, ω, Δ)
    add_asdf!(out, model.model1, ω, Δ)
    add_asdf!(out, model.model2, ω, Δ)
    return nothing
end
function sdf(model::AdditiveTimeSeriesModel{M₁,M₂,1,T}, ω) where {M₁<:TimeSeriesModel{1,T},M₂<:TimeSeriesModel{1,T}} where {T}
    return sdf(model.model1,ω)+sdf(model.model2,ω)
end
function asdf(model::AdditiveTimeSeriesModel{M₁,M₂,1,T}, ω, Δ) where {M₁<:TimeSeriesModel{1,T},M₂<:TimeSeriesModel{1,T}} where {T}
    return asdf(model.model1,ω,Δ)+asdf(model.model2,ω,Δ)
end

extract_asdf(store::SdfStorage) = store.allocatedarray
extract_asdf(store::SdfStorageUni) = store.allocatedarray

function asdf!(store::AdditiveStorage, model::AdditiveTimeSeriesModel)
    @views asdf!(store.store1, model.model1)
    @views asdf!(store.store2, model.model2)
    extract_asdf(extract_S(store.store1)) .+= extract_asdf(extract_S(store.store2))
    return nothing
end

## Gradient of spectral density

### Multivariate ###

"""
    grad_add_sdf!(out, model::TimeSeriesModel, ω)

Add the gradient of the sdf of a model to an input storage vector.
"""
function grad_add_sdf!(out, model::TimeSeriesModel, ω) # default sdf returns error
    error("grad_add_sdf! not yet defined for model of type $(typeof(model)).")
end

"""
    grad_add_asdf!(out, model::TimeSeriesModel, ω, Δ)

Add the gradient of the aliased sdf of a model to an input storage vector.
"""
function grad_add_asdf!(out, model::TimeSeriesModel, ω, Δ) # default method approximates the asdf (overload if known)
    for k ∈ -nalias(model):nalias(model)
        grad_add_sdf!(out, model, ω + k * 2π / Δ)
    end
    return nothing
end

"""
    grad_asdf!(out, model::TimeSeriesModel, ω, Δ)

Compute the gradient of the aliased sdf of a model and store in storage vector.
"""
function grad_asdf!(out, model::TimeSeriesModel, ω, Δ)
    out .= zero(eltype(out))
    grad_add_asdf!(out, model, ω, Δ)
    return nothing
end

"""
    grad_asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel, freq::FreqAcvEst)
    grad_asdf!(store::TimeSeriesModelStorageGradient, model::TimeSeriesModel)

Compute the gradient of the asdf for all frequencies and allocate to appropriate location in storage.
"""
function grad_asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel, freq::FreqAcvEst)
    size(store.allocatedarray,3) == length(freq.Ωₘ) || throw(ArgumentError("size(store.allocatedarray,3) != length(freq.Ωₘ)."))
    @inbounds for i ∈ 1:length(freq.Ωₘ)
        @views grad_asdf!(store.allocatedarray[:, :, i], model, freq.Ωₘ[i], freq.Δ)
    end
    return nothing
end
function grad_asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel{D,T}, freq::FreqAcvEst) where {D,T<:Real} # real case
    size(store.allocatedarray,3) == length(freq.Ωₘ) || throw(ArgumentError("size(store.allocatedarray,3) != length(freq.Ωₘ)."))
    m = length(freq.Ωₘ)
    startcopyind = m÷2+2
    @inbounds for i ∈ 1:startcopyind-1
        @views grad_asdf!(store.allocatedarray[:, :, i], model, freq.Ωₘ[i], freq.Δ)
    end
    @inbounds for i in startcopyind:length(freq.Ωₘ)
        @views store.allocatedarray[:, :, i] .= conj.(store.allocatedarray[:, :, m-i+2])
    end
    return nothing
end
function grad_asdf!(store::TimeSeriesModelStorageGradient, model::TimeSeriesModel)
    grad_asdf!(store.gradmemory, model, store.encodedtime)
    return nothing
end

"""
    grad_sdf(model::TimeSeriesModel, ω)

Compute the gradient of the sdf at `ω`.
"""
function grad_sdf(model::TimeSeriesModel, ω)
    G = zeros(ComplexF64,nlowertriangle_dimension(model),npars(model))
    grad_add_sdf!(G, model, ω)
    return G
end

### Univariate ###

function grad_asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel{1,T}, freq::FreqAcvEst) where {T}
    size(store.allocatedarray,2) == length(freq.Ωₘ) || throw(ArgumentError("size(store.allocatedarray,2) != length(freq.Ωₘ)."))
    @inbounds for i ∈ 1:length(freq.Ωₘ)
        @views grad_asdf!(store.allocatedarray[:, i], model, freq.Ωₘ[i], freq.Δ)
    end
    return nothing
end
function grad_asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel{1,T}, freq::FreqAcvEst) where {T<:Real} # real case
    size(store.allocatedarray,2) == length(freq.Ωₘ) || throw(ArgumentError("size(store.allocatedarray,2) != length(freq.Ωₘ)."))
    m = length(freq.Ωₘ)
    startcopyind = m÷2+2
    @inbounds for i ∈ 1:startcopyind-1
        @views grad_asdf!(store.allocatedarray[:, i], model, freq.Ωₘ[i], freq.Δ)
    end
    @inbounds for i in startcopyind:length(freq.Ωₘ)
        @views store.allocatedarray[:, i] .= store.allocatedarray[:, m-i+2]
    end
    return nothing
end
function grad_sdf(model::TimeSeriesModel{1,T}, ω) where {T}
    G = zeros(ComplexF64,npars(model))
    grad_add_sdf!(G, model, ω)
    return G
end

### Additive ###

function grad_add_sdf!(out, model::AdditiveTimeSeriesModel, ω)
    @views grad_add_sdf!(out[1:npars(M₁)], model.model1, ω)
    @views grad_add_sdf!(out[npars(M₁)+1:end], model.model2, ω)
    return nothing
end

function grad_add_asdf!(out, model::AdditiveTimeSeriesModel, ω, Δ)
    @views grad_add_asdf!(out[1:npars(M₁)], model.model1, ω, Δ)
    @views grad_add_asdf!(out[npars(M₁)+1:end], model.model2, ω, Δ)
    return nothing
end

function grad_asdf!(store::AdditiveStorage, model::AdditiveTimeSeriesModel)
    @views grad_asdf!(store.store1, model.model1)
    @views grad_asdf!(store.store2, model.model2)
    return nothing
end

function grad_sdf(model::AdditiveTimeSeriesModel, ω)
    return hcat(grad_sdf(model.model1, ω), grad_sdf(model.model2, ω)) # not very efficient implementation but not designed to be used.
end
function grad_sdf(model::AdditiveTimeSeriesModel{M₁,M₂,1,T}, ω) where {M₁<:TimeSeriesModel{1,T},M₂<:TimeSeriesModel{1,T}} where {T}
    return vcat(grad_sdf(model.model1, ω), grad_sdf(model.model2, ω)) # not very efficient implementation but not designed to be used.
end

## Hessian of spectral density

"""
    hess_add_sdf!(out, model::TimeSeriesModel, ω)

Add the Hessian of the sdf of a model to an input storage vector.
"""
function hess_add_sdf!(out, model::TimeSeriesModel, ω) # default sdf returns error
    error("hess_add_sdf! not yet defined for model of type $(typeof(model)).")
end

"""
    hess_add_asdf!(out, model::Type{<:TimeSeriesModel}, ω, Δ)

Add the Hessian of the aliased sdf of a model to an input storage vector.
"""
function hess_add_asdf!(out, model::TimeSeriesModel, ω, Δ) # default method approximates the asdf (overload if known)
    for k ∈ -nalias(model):nalias(model)
        hess_add_sdf!(out, model, ω + k * 2π / Δ)
    end
    return nothing
end

"""
    hess_asdf!(out, model::TimeSeriesModel, ω, Δ)

Compute the Hessian of the aliased sdf of a model and store in storage vector.
"""
function hess_asdf!(out, model::TimeSeriesModel, ω, Δ)
    out .= zero(eltype(out))
    hess_add_asdf!(out, model, ω, Δ)
    return nothing
end

"""
    hess_asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel, freq::FreqAcvEst)
    hess_asdf!(store::TimeSeriesModelStorageHessian, model::TimeSeriesModel)

Compute the Hessian of the asdf for all frequencies and allocate to appropriate location in storage.
"""
function hess_asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel, freq::FreqAcvEst)
    size(store.allocatedarray,3) == length(freq.Ωₘ) || throw(ArgumentError("size(store.allocatedarray,3) != length(freq.Ωₘ)."))
    @inbounds for i ∈ 1:length(freq.Ωₘ)
        @views hess_asdf!(store.allocatedarray[:, :, i], model, freq.Ωₘ[i], freq.Δ)
    end
    return nothing
end
function hess_asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel{D,T}, freq::FreqAcvEst) where {D,T<:Real} # real case
    size(store.allocatedarray,3) == length(freq.Ωₘ) || throw(ArgumentError("size(store.allocatedarray,3) != length(freq.Ωₘ)."))
    m = length(freq.Ωₘ)
    startcopyind = m÷2+2
    @inbounds for i ∈ 1:startcopyind-1
        @views grad_asdf!(store.allocatedarray[:, :, i], model, freq.Ωₘ[i], freq.Δ)
    end
    @inbounds for i in startcopyind:length(freq.Ωₘ)
        @views store.allocatedarray[:, :, i] .= conj.(store.allocatedarray[:, :, m-i+2])
    end
    return nothing
end
function hess_asdf!(store::TimeSeriesModelStorageHessian, model::TimeSeriesModel)
    hess_asdf!(store.hessmemory, model, store.encodedtime)
    return nothing
end

"""
    hess_sdf(model::TimeSeriesModel, ω)

Compute the Hessian of the sdf at `ω`.
"""
function hess_sdf(model::TimeSeriesModel, ω)
    H = zeros(ComplexF64,nlowertriangle_dimension(model),nlowertriangle_parameter(model))
    hess_add_sdf!(H, model, ω)
    return H
end

### Univariate ###

function hess_asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel{1,T}, freq::FreqAcvEst) where {T}
    size(store.allocatedarray,2) == length(freq.Ωₘ) || throw(ArgumentError("size(store.allocatedarray,2) != length(freq.Ωₘ)."))
    @inbounds for i ∈ 1:length(freq.Ωₘ)
        @views hess_asdf!(store.allocatedarray[:, i], model, freq.Ωₘ[i], freq.Δ)
    end
    return nothing
end
function hess_asdf!(store::TimeSeriesModelStorage, model::TimeSeriesModel{1,T}, freq::FreqAcvEst) where {T<:Real} # real case
    size(store.allocatedarray,2) == length(freq.Ωₘ) || throw(ArgumentError("size(store.allocatedarray,2) != length(freq.Ωₘ)."))
    m = length(freq.Ωₘ)
    startcopyind = m÷2+2
    @inbounds for i ∈ 1:startcopyind-1
        @views grad_asdf!(store.allocatedarray[:, i], model, freq.Ωₘ[i], freq.Δ)
    end
    @inbounds for i in startcopyind:length(freq.Ωₘ)
        @views store.allocatedarray[:, i] .= store.allocatedarray[:, m-i+2]
    end
    return nothing
end
function hess_sdf(model::TimeSeriesModel{1,T}, ω) where {T}
    H = zeros(ComplexF64,nlowertriangle_parameter(model))
    hess_add_sdf!(H, model, ω)
    return H
end

### Additive ###

function hess_asdf!(store::AdditiveStorage, model::AdditiveTimeSeriesModel)
    @views hess_asdf!(store.store1, model.model1)
    @views hess_asdf!(store.store2, model.model2)
    return nothing
end

function hess_sdf(model::AdditiveTimeSeriesModel, ω)
    return hcat(hess_sdf(model.model1, ω), hess_sdf(model.model2, ω)) # not very efficient implementation but not designed to be used.
end
function hess_sdf(model::AdditiveTimeSeriesModel{M₁,M₂,1,T}, ω) where {M₁<:TimeSeriesModel{1,T},M₂<:TimeSeriesModel{1,T}} where {T}
    return vcat(hess_sdf(model.model1, ω), hess_sdf(model.model2, ω)) # not very efficient implementation but not designed to be used.
end


text:
abstract type RandomGraphType end
struct ErdosRenyi <: RandomGraphType end
struct RandomGeometric <: RandomGraphType end

# --------------------------------------------------------------------------- #
#						       Random Graph Models
# --------------------------------------------------------------------------- #

function erdos_renyi(n,p)
    A = sprand(n,n,p)
    U = triu(A,1)
    is,js,_ = findnz(max.(U,U'))
    return sparse(is,js,ones(length(is)),n,n)
end

function erdos_renyi(n;degreedist=LogNormal(log(5),1))

    println("made it")
    # form the edges for sparse
    ei = Int[]
    ej = Int[]

    for i=1:n
      deg = ceil(Int,minimum((rand(degreedist),n-1)))
      neighbors = sample(1:n,deg+1, replace = false)

      for j in neighbors
        if i != j
          push!(ei,i)
          push!(ej,j)
        end
      end
    end

    A = sparse(ei,ej,1.0,n,n)

    return max.(A,A')
end

#--                             David's graph code                          --#

function random_geometric_graph(n,k)
  xy = rand(2,n)
  T = BallTree(xy)
  idxs = knn(T, xy, minimum((k,n)))[1]
  # form the edges for sparse
  ei = Int[]
  ej = Int[]
  for i=1:n
    for j=idxs[i]
      if i != j
        push!(ei,i)
        push!(ej,j)
      end
    end
  end
  A = sparse(ei,ej,1.0,n,n)

  return max.(A,A')
end

function spatial_network(n::Integer, d::Integer; degreedist=LogNormal(log(5),1))
  xy, ei, ej = spatial_graph_edges(n, d;degreedist=degreedist)
  A = sparse(ei,ej,1.0,n,n)
  return max.(A,A')
end

function spatial_graph_edges(n::Integer,d::Integer;degreedist=LogNormal(log(5),1))
  xy = rand(d,n)
  T = BallTree(xy)
  # form the edges for sparse
  ei = Int[]
  ej = Int[]
  for i=1:n
    deg = ceil(Int,minimum((rand(degreedist),n-2)))
    idxs, dists = knn(T, xy[:,i], deg+1)
    for j in idxs
      if i != j
        push!(ei,i)
        push!(ej,j)
      end
    end
  end
  return xy, ei, ej
end

ext: jl

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
using OrdinaryDiffEq, DiffEqCallbacks, Test

function some_dynamics(u,p,t)
    du = zeros(length(u))
    for i in 1:length(u)
        for j in 1:length(u)
            du[i] += p[i,j]*u[j]*u[i]
        end
    end
    return du
end

p = rand(4,4)
startp = copy(p)

u0 = zeros(length(p[1,:])) .+ 0.1
tspan = (0.,1.)
prob = ODEProblem(some_dynamics,u0,tspan,p)
cb = PresetTimeCallback(0.5,integrator -> integrator.p .= rand(4,4))
sol = solve(prob,Tsit5(),callback = cb)
@test 0.5 ∈ sol.t
@test p != startp

p = rand(4,4)
startp = copy(p)

prob = ODEProblem(some_dynamics,u0,tspan,p)
cb = PresetTimeCallback([0.3,0.6],integrator -> integrator.p .= rand(4,4))
sol = solve(prob,Tsit5(),callback = cb)
@test 0.3 ∈ sol.t
@test 0.6 ∈ sol.t
@test p != startp

notcalled = true
prob = ODEProblem(some_dynamics,u0,tspan,p)
cb = PresetTimeCallback([1.2],integrator -> notcalled = false)
sol = solve(prob,Tsit5(),callback = cb)
@test notcalled

cb = PresetTimeCallback([1.2],integrator -> begin
    global notcalled
    notcalled = false
    end, filter_tstops = false)
sol = solve(prob,Tsit5(),callback = cb)
@test !notcalled

notcalled = true
prob = ODEProblem(some_dynamics,u0,(1.0,0.0),p)
cb = PresetTimeCallback([-0.2],integrator -> notcalled = false)
sol = solve(prob,Tsit5(),callback = cb)
@test notcalled

cb = PresetTimeCallback([-0.2],integrator -> begin
    global notcalled
    notcalled = false
    end, filter_tstops = false)
sol = solve(prob,Tsit5(),callback = cb)
@test !notcalled


text:
include("init-namedarrays.jl")

@test names(n) == Array[["one", "two"], ["a", "b", "c", "d"]]
@test defaultnames(n.array) == Array[["1", "2"], ["1", "2", "3", "4"]]

@test dimnames(n) == [:A, :B]
@test dimnames(n.array) == [:A, :B]

@test dimnames(n.array, 1) == :A
@test dimnames(n.array, 2) == :B

dn1 = ["一", "二"]
dn2 = ["一", "二", "三", "四"]

@inferred setnames!(n, dn1, 1)
for (i, zh) in enumerate(dn2)
	@inferred setnames!(n, zh, 2, i)
end

@test names(n) == Array[dn1, dn2]

@inferred setdimnames!(n, ("thing1", :thing2))
@test dimnames(n) == ["thing1", :thing2]
setdimnames!(n, ["magnificent", 7]) ## can't infer mixed types
@test dimnames(n) == ["magnificent", 7]

println(n)

@test_throws TypeError setnames!(n, [:a, :v], 1)
@test_throws TypeError setnames!(n, :a, 1, 1)
@test_throws DimensionMismatch setnames!(n, ["a"], 1)
for i in [0,5], j in 1:2
	@test_throws BoundsError setnames!(n, "a", j, i)
end


text:
using Phylodendron2

HLINE = '─'
VLINE = '│'
TOP_ELBOW = '┌'
THREEWAY = '├'
BOTTOM_ELBOW = '└'
NODE_SYMBOL = '○'

txt = "((((🐇:1,Whale🐳:2):1,(Snake:1,🐢:2):1):6,🐠:1):1,🐌:5);"
txt = "(((A:4,B:1,C:2,D:3):1,E:2):1,(F:1,G:2):4);"
tree = read_newick(text=txt)[1]

function branch2edge(
	p::Node, q::Node, anc::Int, i::Int, edgematrix, labels, branchlengths, nanc
	)
	i += 1		# Last node index
	j = i		# Index of this node
	edgematrix = vcat(edgematrix, [anc j])
	push!(labels, label(q))
	push!(branchlengths, brlength(p, q))
	push!(nanc, nanc[anc] + 1)

	for link in q.links
		link.to == p && continue
		i, edgematrix = branch2edge(
			q, link.to, j, i, edgematrix, labels, branchlengths, nanc
			)
	end

	return i, edgematrix
end

function branch2edge(tree::Tree)
	i = 1
	anc = 1
	edgematrix = Matrix{Int}(undef, 0, 2)
	nodelabels = String[]			# Node labels
	branchlengths = Union{Float64, Nothing}[]		# Branch lengths
	nanc = [0]						# Number of ancestors of the node

	push!(nodelabels, label(tree.origin))
	push!(branchlengths, 0.0)

	for link in tree.origin.links
		i, edgematrix =
			branch2edge(tree.origin, link.to, anc, i, edgematrix, nodelabels, branchlengths, nanc)
	end

	return (
		edgematrix = edgematrix,
		nodelabels = nodelabels,
		branchlengths = branchlengths,
		nanc = nanc
		)
end

function children_edgematrix(edgematrix, i)
	edges = findall(==(i), edgematrix[:, 1])

	return edgematrix[edges, 2]
end

function compute_ypos(edgematrix, n)
	ypos = zeros(Int, n)
	tips = [i ∉ edgematrix[:, 1] for i in 1:n]

	y = 1
	for i in eachindex(ypos)
		if tips[i]
			ypos[i] = y
			y += 2
		end
	end

	for i in n:-1:1
		chld = children_edgematrix(edgematrix, i)
		if ! isempty(chld)
			y = floor(Int, sum(extrema(ypos[chld])) / 2)
			ypos[i] = y
		end
	end

	return ypos
end

function compute_xpos(edgematrix, brlen, n, s)
	nedge = size(edgematrix, 1)

	xposanc = ones(Int, n)		# Position of the ancestor node
	xposdesc = ones(Int, n)		# Position of the descendant node

	for i in 1:nedge		# start on 2nd index to skip the root
		anc = edgematrix[i, 1]
		desc = edgematrix[i, 2]
		xposanc[desc] = xposdesc[anc]
		brdesc = floor(Int, (isnothing(brlen[desc]) ? 0 : brlen[desc]) * s)
		xposdesc[desc] = xposanc[desc] + brdesc + 1
	end

	return xposanc, xposdesc
end

function paint_canvas!(A::Matrix{Char}, edgematrix, ypos, xposanc, xposdesc, nodelabels)
	for i in eachindex(ypos)
		x0, x1, y = xposanc[i], xposdesc[i], ypos[i]
		x1 == 0 && continue

		# Draw vertical lines of the fork
		chld = children_edgematrix(edgematrix, i)
		if length(chld) > 1
			ytop, ybottom = extrema(ypos[chld])
			A[ytop+1:ybottom-1, x1] .= VLINE
			A[ypos[chld], x1] .= THREEWAY
			A[ytop, x1] = TOP_ELBOW
			A[ybottom, x1] = BOTTOM_ELBOW
		end

		# Draw node
		A[y, x1] = NODE_SYMBOL

		thislab = nodelabels[i]
		if ! isempty(thislab)
			for j in eachindex(thislab)
				A[y, x1 + 1 + j] = thislab[j]
			end
		end

		# Draw branch
		x0 += 1
		x1 -= 1
		if x1 ≥ x0
			A[y, x0:x1] .= HLINE
		end
	end
end

function print_canvas(io::IO, A::Array{Char})
	nline, ncol = size(A)
	for i in 1:nline
		for j in 1:ncol
			print(io, A[i, j])
		end
		print(io, '\n')
	end

	return nothing
end

function textplot(io::IO, tree::Tree)
	edgematrix, nodelabels, brlen, nanc = branch2edge(tree)

	n = length(nodelabels)

	ypos = compute_ypos(edgematrix, n)

	canvas_width = displaysize(io)[2]
	minimum_widths = nanc .+ length.(nodelabels) .+ 2
	minimum_width = maximum(minimum_widths)

	minimum_width > canvas_width && @error "the tree is too big to be displayed here"

	# Compute scaling factor
	s = (canvas_width  .- minimum_widths) ./ brlen |> minimum
	xposanc, xposdesc = compute_xpos(edgematrix, brlen, n, s)

	canvas_height = maximum(ypos)
	canvas = fill(' ', canvas_width, canvas_height);

	paint_canvas!(canvas, edgematrix, ypos, xposanc, xposdesc, nodelabels);
	print_canvas(stdout, canvas)

	return nothing
end

textplot(tree::Tree) = textplot(stdout, tree)

textplot(tree)


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
using GenerativeAD.Models: anomaly_score, generalized_anomaly_score_gpu
using BSON
using PyCall
using StatsBase: fit!, predict, sample

using Flux
using MLDataPattern
using OrderedCollections

s = ArgParseSettings()
@add_arg_table! s begin
	"max_seed"
		default = 1
		arg_type = Int
		help = "maximum number of seeds to run through"
	"dataset"
		default = "iris"
		arg_type = String
		help = "dataset"
	"sampling"
		default = "random"
		arg_type = String
		help = "sampling of hyperparameters"
	"contamination"
		arg_type = Float64
		help = "contamination rate of training data"
		default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, max_seed, sampling, contamination = parsed_args

modelname ="GANomaly"

function sample_params()
	parameters_rng = (
		hdim   		= 2 .^(4:9),
		zdim   		= 2 .^(1:8),
		nlayers   	= 3:4,
		activation  = ["relu", "swish", "tanh"],
		lr   		= 10f0 .^ (-4:-3),
		decay   	= 0f0:0.1:0.5,
		batch_size  = 2 .^ (5:7),
		w1 			= [1,10:10:90 ...],
		w2 			= [1,10:10:90 ...],
		w3 			= [1,10:10:90 ...],
		init_seed   = 1:Int(1e8),
	)

	return (;zip(keys(parameters_rng), map(x->sample(x, 1)[1], parameters_rng))...)
end

function create_space()
	pyReal = pyimport("skopt.space")["Real"]
	pyInt = pyimport("skopt.space")["Integer"]
	pyCat = pyimport("skopt.space")["Categorical"]

	parameters_rng = (
		hdim   		= pyInt(4, 9,						          name="log2_hdim"),
		zdim   		= pyInt(1, 8,						          name="log2_zdim"),
		nlayers   	= pyInt(3, 4,						          name="nlayers"),
		activation  = pyCat(categories=["relu", "swish", "tanh"], name="activation"),
		lr   		= pyReal(1f-5, 1f-3, prior="log-uniform",	  name="lr"),
		#decay   	= pyReal(1f-5, 5f-1, prior="log-uniform", 	  name="decay"),
		batch_size  = pyInt(5, 7, 								  name="log2_batch_size"),
		w1			= pyInt(1, 90, 								  name="w1"),
		w2			= pyInt(1, 90, 								  name="w2"),
		w3			= pyInt(1, 90, 								  name="w3"),
	)
end

function fit(data, parameters)
	weights = [parameters[:w1],parameters[:w2],parameters[:w3]] # change :weights -> :w1, :w2, :w3 because of bayes optimization
	default_params = (iters=10000, check_every=30, patience=10, weights=weights,)

	all_parameters = merge(parameters, default_params)

	# define models (Generator, Discriminator)
	generator, discriminator = GenerativeAD.Models.tabular_ganomaly_constructor(all_parameters)

	# define optimiser
	try
		global info, fit_t, _, _, _ = @timed fit!(generator|>gpu, discriminator|>gpu, data, all_parameters)
	catch e
		println("Error caught.")
		return (fit_t = NaN, model = nothing, history = nothing, n_parameters = NaN), []
	end

	training_info = (
		fit_t = fit_t,
		model = (info[2]|>cpu, info[3]|>cpu),
		history = info[1], # losses through time
		npars = info[4], # number of parameters
		iters = info[5] # optim iterations of model
		)

	return training_info, [(x -> GenerativeAD.Models.anomaly_score(generator|>cpu, x; dims=1)[:], parameters)]
end

#_________________________________________________________________________________________________

try_counter = 0
max_tries = 10*max_seed
cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
sampling_string = sampling == "bayes" ? "_bayes" : ""
prefix = "experiments$(sampling_string)/tabular$(cont_string)"
dataset_folder = datadir("$(prefix)/$(modelname)/$(dataset)")
while try_counter < max_tries
	if sampling == "bayes"
		parameters = GenerativeAD.bayes_params(
								create_space(),
								dataset_folder,
								sample_params; add_model_seed=true)
	else
		parameters = sample_params()
	end

    for seed in 1:max_seed
		savepath = joinpath(dataset_folder, "seed=$(seed)")
		mkpath(savepath)

		# get data
		data = GenerativeAD.load_data(dataset, seed=seed, contamination=contamination)

		# update parameter
		parameters = merge(parameters, (idim=size(data[1][1],1), ))
		# here, check if a model with the same parameters was already tested
		if GenerativeAD.check_params(savepath, parameters)
			#(X_train,_), (X_val, y_val), (X_test, y_test) = data
			training_info, results = fit(data, parameters)
			# saving model separately
			if training_info.model !== nothing
				tagsave(joinpath(savepath, savename("model", parameters, "bson", digits=5)), Dict("model"=>training_info.model), safe = true)
				training_info = merge(training_info, (model = nothing,))
			end
			save_entries = merge(training_info, (modelname = modelname, seed = seed, dataset = dataset, contamination = contamination))
			# now loop over all anomaly score funs
			all_scores = [GenerativeAD.experiment(result..., data, savepath; save_entries...) for result in results]
			if sampling == "bayes" && length(all_scores) > 0
				@info("Updating cache with $(length(all_scores)) results.")
				GenerativeAD.update_bayes_cache(dataset_folder,
						all_scores; ignore=Set([:init_seed, :L, :idim]))
			end
			global try_counter = max_tries + 1
		else
			@info "Model already present, sampling new hyperparameters..."
			global try_counter += 1
		end
	end
end
(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing


text:
include("scaler.jl")
using DataFrames

export fill_missing_aggregate
export fill_missing_scalar

"""
_given a dataframe and optionally columns, return the dataframe with missing values filled in with column's mean, median, etc._
#### parameters:
    data : DataFrame
        dataframe to min max scale
    columns : Array{Symbol}
        an array of column names to scale, default is all numeric
    method : Function
        [mean, median, mode], or any desired function that operates on an array and returns a single value
#### returns: DataFrame
    a new dataframe with desired columns altered
"""
function fill_missing_aggregate(data::DataFrame, columns::Array=[]; method)
    tmp = copy(data)
    numeric_cols = numeric_columns(tmp)
    iter_columns = isempty(columns) ? numeric_cols :
                   intersect(numeric_cols, Symbol.(columns))
    for col in iter_columns
        tmp[!, Symbol(col)] = recode(tmp[!, Symbol(col)], missing => method(skipmissing(tmp[!, Symbol(col)])))
    end
    return tmp
end

"""
_given a dataframe and optionally columns, return the dataframe with missing values filled in with provided scalar_
#### parameters:
    data : DataFrame
        dataframe to min max scale
    columns : Array{Symbol}
        an array of column names to scale, default is all numeric
    scalar : Any
        any desired value to replace missings, must be of type of column(s)
#### returns: DataFrame
    a new dataframe with desired columns altered
"""
function fill_missing_scalar(data::DataFrame, columns::Array=[]; scalar::Any)
    tmp = copy(data)
    iter_columns = isempty(columns) ? names(tmp) :
                   intersect(names(tmp), Symbol.(columns))
    for col in iter_columns
        try
            tmp[!, Symbol(col)] = recode(tmp[!, Symbol(col)], missing => scalar)
        catch e
            throw(ArgumentError("Scalar must match variable type of ALL columns."))
        end
    end
    return tmp
end


text:
using BenchmarkTools
using MAVs, NLOptControl, VehicleModels

# change defaults
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 2.50
BenchmarkTools.DEFAULT_PARAMETERS.samples = 1

const SUITE = BenchmarkGroup()

#case_names = ["s1","s2","s3","s4","s5","s6"]
SUITE["test"] = BenchmarkGroup(["case","solver"])
planner_name = "RTPP"
vehicle_name = "hmmwv"
models =[:ThreeDOFv2, :KinematicBicycle2]
model = models[1]

for case_name in ("s1","s2")
    for solver in (:KNITRO,:Ipopt)
        c = load(open(string(Pkg.dir("MAVs"),"/config/planner/",planner_name,".yaml")))
        c["vehicle"] = load(open(string(Pkg.dir("MAVs"),"/config/vehicle/",vehicle_name,".yaml")))
        c["goal"] = load(open(string(Pkg.dir("MAVs"),"/config/case/",case_name,".yaml")))["goal"]
        c["X0"] = load(open(string(Pkg.dir("MAVs"),"/config/case/",case_name,".yaml")))["X0"]
        c["obstacle"] = load(open(string(Pkg.dir("MAVs"),"/config/case/",case_name,".yaml")))["obstacle"]
        setConfig(c, "misc";(:N=>15),(:model=>model),(:solver=>solver), (:integrationScheme=>:trapezoidal))
        fixYAML(c)   # fix messed up data types
        n = initializeAutonomousControl(c);
        SUITE["test"][case_name,string(solver)] = @benchmarkable simMPC!($n;updateFunction=updateAutoParams!,checkFunction=checkCrash)
    end
end

run(SUITE)


text:
"""
    FixedRateDescent

Descent method which minimizes the objective function in the direction
of antigradient at each step.
"""
mutable struct FixedRateDescent{T<:AbstractFloat,V<:AbstractVector{T}} <: CoreMethod
    x::V
    g::V
    α::T
end

@inline gradientvec(M::FixedRateDescent) = M.g
@inline argumentvec(M::FixedRateDescent) = M.x
@inline step_origin(M::FixedRateDescent) = M.x

function FixedRateDescent(x::AbstractVector, α::Real)
    F = float(eltype(x))
    FixedRateDescent(similar(x, F), similar(x, F), convert(F, α))
end

FixedRateDescent(x::AbstractVector) = FixedRateDescent(x, 0.01)

function init!(M::FixedRateDescent{T}, optfn!, x0; kw...) where {T}
    optfn!(x0, zero(T), x0)
    return
end

@inline reset!(::FixedRateDescent) = nothing

function reset!(M::FixedRateDescent, x0, α = M.α)
    if length(M.x) != length(x0)
        foreach((M.x, M.g)) do v
            resize!(v, length(x0))
        end
    end
    M.α = α
    return
end

@inline function callfn!(M::FixedRateDescent, fdf, x, α, d)
    __update_arg!(M, x, α, d)
    y, g = fdf(M.x, M.g)
    __update_grad!(M, g)
    return y, g
end

function step!(M::FixedRateDescent, optfn!; constrain_step = infstep)
    d = rmul!(M.g, -1)
    maxstep = constrain_step(M.x, d)
    s = M.α <= maxstep ? M.α : maxstep / 2
    optfn!(M.x, s, d)
    return s
end

@inline function __update_arg!(M::FixedRateDescent, x, α, d)
    map!(M.x, d, x) do a, b
        muladd(α, a, b)
    end
    return
end

@inline function __update_arg!(M::FixedRateDescent, x)
    if x !== M.x
        copy!(M.x, x)
    end
    return
end

@inline function __update_grad!(M::FixedRateDescent, g)
    if M.g !== g
        copy!(M.g, g)
    end
    return
end

"""
    MomentumDescent

Descent method which minimizes the objective function in the direction
of antigradient at each step.
"""
mutable struct MomentumDescent{T<:AbstractFloat,V<:AbstractVector{T}} <: CoreMethod
    x::V
    g::V
    v::V
    learn_rate::T
    decay_rate::T
end

@inline gradientvec(M::MomentumDescent) = M.g
@inline argumentvec(M::MomentumDescent) = M.x
@inline step_origin(M::MomentumDescent) = M.x

function MomentumDescent(
    x::AbstractVector;
    learn_rate::Real = 0.01, decay_rate::Real = 0.9
)
    F = float(eltype(x))
    return MomentumDescent(
        similar(x, F),
        similar(x, F),
        similar(x, F),
        convert(F, learn_rate),
        convert(F, decay_rate),
    )
end

function init!(M::MomentumDescent{T}, optfn!, x0; kw...) where {T}
    optfn!(x0, zero(T), x0)
    fill!(M.v, zero(T))
    return M
end

@inline function reset!(M::MomentumDescent)
    fill!(M.v, false)
    return M
end

function reset!(
    M::MomentumDescent,
    x0,
    learn_rate = M.learn_rate,
    decay_rate = M.decay_rate,
)
    n = length(x0)
    if length(M.x) != n
        for v in (M.x, M.g, M.v)
            resize!(v, n)
        end
    end
    M.v .= false
    M.x .= x0
    M.learn_rate = learn_rate
    M.decay_rate = decay_rate
    return M
end

@inline function callfn!(M::MomentumDescent, fdf, x, α, d)
    __update_arg!(M, x, α, d)
    y, g = fdf(M.x, M.g)
    __update_grad!(M, g)
    return y, g
end

function step!(M::MomentumDescent, optfn!; constrain_step = infstep)
    map!(M.v, M.v, M.g) do v, g
        M.decay_rate * v - M.learn_rate * g
    end
    maxstep = constrain_step(M.x, M.v)
    if maxstep <= 1
        M.v *= maxstep / 2
    end
    optfn!(M.x, one(maxstep), M.v)
    return M.learn_rate
end

@inline function __update_arg!(M::MomentumDescent, x, α, d)
    map!(M.x, d, x) do a, b
        muladd(α, a, b)
    end
    return
end

@inline function __update_arg!(M::MomentumDescent, x)
    if x !== M.x
        copy!(M.x, x)
    end
    return
end

@inline function __update_grad!(M::MomentumDescent, g)
    if M.g !== g
        copy!(M.g, g)
    end
    return
end

"""
    NesterovMomentum

Descent method which minimizes the objective function in the direction
of antigradient at each step.
"""
mutable struct NesterovMomentum{T<:AbstractFloat,V<:AbstractVector{T}} <: CoreMethod
    x::V
    g::V
    v::V
    α::T # learning rate
    β::T # decay rate
end

@inline gradientvec(M::NesterovMomentum) = M.g
@inline argumentvec(M::NesterovMomentum) = M.x
@inline step_origin(M::NesterovMomentum) = M.x

function NesterovMomentum(x::AbstractVector; learn_rate::Real=0.01, decay_rate::Real=0.9)
    F = float(eltype(x))

    # α and β are parameters of "half step" v <- βv - α∇f
    # applied twice, it must yield v <- decay_rate × v - learn_rate × ∇f
    α = learn_rate
    β = decay_rate
    return NesterovMomentum(
        similar(x, F),
        similar(x, F),
        similar(x, F),
        convert(F, α),
        convert(F, β),
    )
end

function init!(M::NesterovMomentum{T}, optfn!, x0; kw...) where {T}
    optfn!(x0, zero(T), x0)
    fill!(M.v, zero(T))
    return M
end

@inline function reset!(M::NesterovMomentum)
    fill!(M.v, false)
    return M
end

function reset!(M::NesterovMomentum, x0, learn_rate=M.learn_rate, decay_rate=M.decay_rate)
    n = length(x0)
    if length(M.x) != n
        for v in (M.x, M.g, M.v)
            resize!(v, n)
        end
    end
    M.α = learn_rate / (1 + decay_rate)
    M.β = sqrt(decay_rate)
    return M
end

@inline function callfn!(M::NesterovMomentum, fdf, x, α, d)
    __update_arg!(M, x, α, d)
    y, g = fdf(M.x, M.g)
    __update_grad!(M, g)
    return y, g
end

function step!(M::NesterovMomentum, optfn!; constrain_step = infstep)
    α, β = M.α, M.β
    M.v *= β
    maxstep = constrain_step(M.x, M.v)
    if maxstep <= 1
        M.v *= maxstep / 2
    end
    optfn!(M.x, one(maxstep), M.v)
    M.v .-= α * M.g
    M.x .-= α * M.g
    return M.α
end

@inline function __update_arg!(M::NesterovMomentum, x, α, d)
    map!(M.x, d, x) do a, b
        muladd(α, a, b)
    end
    return
end

@inline function __update_arg!(M::NesterovMomentum, x)
    if x !== M.x
        copy!(M.x, x)
    end
    return
end

@inline function __update_grad!(M::NesterovMomentum, g)
    if M.g !== g
        copy!(M.g, g)
    end
    return
end


text:
using LaserTag
using Plots
using POMDPToolbox
using Reel
using ProgressMeter
using PmapProgressMeter
using ParticleFilters
using ContinuousPOMDPTreeSearchExperiments
using Plots
using QMDP
using JLD
# using DESPOT
using BasicPOMCP
using ARDESPOT

file_contents = readstring(@__FILE__())

@everywhere begin
    using POMDPs
    using POMDPToolbox
    using ContinuousPOMDPTreeSearchExperiments
    using DiscreteValueIteration
    using ParticleFilters
    using POMCPOW
    using LaserTag
    using QMDP
    # using DESPOT
    using BasicPOMCP
    using ARDESPOT

    N = 1000
    n = 1_000_000
    P = typeof(gen_lasertag(rng=MersenneTwister(18)))

    solvers = Dict{String, Union{Policy, Solver}}(

        #=
        "pomcpow" => begin
            # ro = MoveTowards()
            solver = POMCPOWSolver(tree_queries=1_000_000_000, #500_000
                                   criterion=MaxUCB(60.0),
                                   final_criterion=MaxTries(),
                                   max_depth=90,
                                   max_time=1.0,
                                   enable_action_pw=false,
                                   # k_action=4.0,
                                   # alpha_action=1/8,
                                   k_observation=2.0,
                                   alpha_observation=1/20,
                                   estimate_value=FOValue(ValueIterationSolver()),
                                   check_repeat_act=false,
                                   check_repeat_obs=false,
                                   init_N=InevitableInit(),
                                   init_V=InevitableInit(),
                                   rng=MersenneTwister(13)
                                  )
            solver
        end,
        =#

        # "move_towards_sampled" => MoveTowardsSampled(MersenneTwister(17)),

        # "ml" => OptimalMLSolver(ValueIterationSolver()),

        # "be" => BestExpectedSolver(ValueIterationSolver()),

        # "random" => RandomSolver(rng=MersenneTwister(4)),

        "despot" => begin
            DESPOTSolver(lambda=0.01,
                         K=500,
                         max_trials=1_000_000,
                         T_max=1.0,
                         bounds=LaserBounds{P}(),
                         default_action=NoGapTag(),
                         bounds_warnings=false,
                         random_source=MemorizingSource(500, 90, MersenneTwister(5), min_reserve=10)
                         rng=MersenneTwister(4))
        end,

        # "qmdp" => QMDPSolver(max_iterations=1000),

        #=
        "pomcp" => POMCPSolver(tree_queries=n,
                                   c=20.0,
                                   max_depth=100,
                                   estimate_value=FOValue(ValueIterationSolver()),
                                   rng=MersenneTwister(13)
                                  )
        =#
    )
end

@show N
@show solvers["despot"].lambda
@show solvers["despot"].K
@show solvers["despot"].T_max
@show solvers["pomcpow"].max_time

rdict = Dict{String, Any}()
for (k,sol) in solvers
    prog = Progress(N, desc="Simulating...")
    @show k
    rewards = pmap(prog, 1:N) do i
    # rewards = map(1:N) do i
        pomdp = gen_lasertag(rng=MersenneTwister(i+700_000))
        if isa(sol,Solver)
            p = solve(deepcopy(sol), pomdp)
        else
            p = sol
        end
        hr = HistoryRecorder(max_steps=100, rng=MersenneTwister(i))
        up_rng = MersenneTwister(i+140_000)
        up = ObsAdaptiveParticleFilter(deepcopy(pomdp), LowVarianceResampler(100_000), 0.05, up_rng)
        hist = simulate(hr, pomdp, p, up)
        discounted_reward(hist)
    end
    @show mean(rewards)
    @show std(rewards)/sqrt(N)
    rdict[k] = rewards
end

filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "data", "laser_discrete_run_$(Dates.format(now(), "E_d_u_HH_MM")).jld")
println("saving to $filename...")
@save(filename, solvers, rdict, file_contents)
println("done.")


text:
using ProximalOperators

A = [  1.0  -2.0   3.0  -4.0  5.0;
       2.0  -1.0   0.0  -1.0  3.0;
      -1.0   0.0   4.0  -3.0  2.0;
      -1.0  -1.0  -1.0   1.0  3.0]
b = [1.0, 2.0, 3.0, 4.0]

m, n = size(A)

f = Translate(LogisticLoss(ones(m)), -b)
lam = 0.1
g = NormL1(lam)

x_star = [0, 0, 2.114635341704963e-01, 0, 2.845881348733116e+00]

# Nonfast/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fs=f, As=A, g=g, tol=1e-6, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 1700
println(sol)

# Fast/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fs=f, As=A, g=g, tol=1e-6, adaptive=true, fast=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 500
println(sol)

# ZeroFPR/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fs=f, As=A, g=g, tol=1e-6, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 25
println(sol)

# PANOC/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.PANOC(x0; fs=f, As=A, g=g, tol=1e-6, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 35
println(sol)


text:
using LinearAlgebra, Random
using Statistics, Distributions

"""
@article{Lewandowski2009a,
    Author = {Daniel Lewandowski and Dorota Kurowicka and Harry Joe},
    Journal = {Journal of Multivariate Analysis},
    Number = {9},
    Pages = {1989 - 2001},
    Title = {Generating random correlation matrices based on vines and extended onion method},
    Volume = {100},
    Year = {2009}}
https://doi.org/10.1016/j.jmva.2009.04.008
"""
function onion(n; η = 1.0)
    S = Symmetric(ones(n, n))
    β = η + (n - 2) / 2
    S.data[1, 2] = 2 * rand(Beta(β, β)) - 1
    for k = 2:n-1
        β -= 0.5
        y = rand(Beta(k / 2, β))
        u = normalize!(Random.randn(k))
        w = sqrt(y) * u
        F = cholesky(Symmetric(S[1:k, 1:k]))
        S.data[1:k, k+1] = F.L * w
    end
    return S
end

function randncm(n; seed = 0, γ = 0.0, p = 0.5, gaussian_noise = false)
    rng = Random.seed!(seed)

    # Random target matrix
    U = onion(n)

    # Random sparse H
    if gaussian_noise
        # Hij ∈ [1, 10] for Hij ≠ 0
        Htmp = [rand() < p ? 1.0 + 9.0 * rand() : 0.0 for i = 1:n, j = 1:n]
    else
        Htmp = [rand() < p ? rand() : 0.0 for i = 1:n, j = 1:n]
    end
    H = Symmetric(triu(Htmp, 1) + I)

    # Random noise
    if gaussian_noise
        # G = U + E./H
        G = copy(U)
        for j = 1:n
            for i = 1:j-1
                Hij = H.data[i, j]
                if Hij > 0.0
                    G.data[i, j] += γ * Random.randn() / Hij
                else
                    G.data[i, j] = 0.0
                end
            end
        end
    else
        Etmp = 2 * rand(n, n) .- 1
        E = Symmetric(triu(Etmp, 1) + I)
        Gtmp = (1 - γ) .* U .+ γ .* E
        G = Symmetric(triu(Gtmp, 1) + I)
    end

    Random.seed!(rng)

    return U, G, H
end

max_stars_repo_name: krislock/IR-FISTA

max_issues_repo_name: krislock/IR-FISTA

max_forks_repo_name: krislock/IR-FISTA


text:
using MAT
using jInv.Mesh
using jInv.LinearSolvers
using jInv.ForwardShare
using jInv.InverseSolve
using jInv.Utils
using OED

model   = 1;
IPmodes = [0;1;2;3];
lambdas = logspace(-4.,2,20);
nrep    = 4;

r       = 48
l       = 2
pad     = 4
nex     = 20
matfile = "data/OED_shepp.mat"
resfile = "results/plotObjectiveB/OED_shepp"
matfile = matread(matfile)
n       = Int(matfile["n"])
L       = .1*speye(n*n)
Xtrue   = matfile["x_true"]

Ce = matfile["Ce"]
ce = Ce*Xtrue
Ci = matfile["Ci"]
ci = matfile["ci"]

if nex!==-1
	idx = randperm(size(Xtrue,2))[1:nex]
	Xtrue = Xtrue[:,idx]
	ci   = ci[:,idx]
	ce   = ce[:,idx]
end

Mimg   = getRegularMesh([-1 1 -1 1.],[n,n])
domain = copy(Mimg.domain);
domain[1:2:end] -= pad*Mimg.h
domain[2:2:end] += pad*Mimg.h
Mgrid = getRegularMesh(domain,[n+2*pad,n+2*pad])
wGrid = getRegularMesh([0 180 0 180],[32 32]);

# divide problems among available workers
ids = mod(0:size(Xtrue,2)-1,nworkers())+1
OEDparam = Array{OEDTypeBParam}(nworkers())
for k=1:nworkers()
	OEDparam[k] = getOEDTypeBParam(r,l,Xtrue[:,ids.==k],Mimg,Mgrid)
	OEDparam[k].noiseLevel = 1e0
end
ipDict  = Dict(0=>:unconstr, 1=>:eqconstr, 2=>:nonnegconstr, 3=>:boxconstr)

pFor = 1;
for k=3:3
	IPmode  = ipDict[IPmodes[k]]
	outfile = @sprintf "%s-%s.mat" resfile IPmode
	outf    = matread(outfile)
    F       = outf["F"]
	F[F.==0] = NaN
	w1,w2 = getCellCenteredAxes(wGrid)
	i,j     = ind2sub((length(w1),length(w2)), indmin(F))
	wOpt = [w1[i];w2[j]]

	println("wOpt = $wOpt")
	if IPmode==:unconstr
		pFor = Array{UnconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
			pFor[j] = UnconstrainedOEDParam(OEDparam[j],L,zeros(0,0),[])
		end
	elseif IPmode==:eqconstr
		if isempty(Ce) || isempty(ce)
			error("equality constraints must be provided. To run unconstrained OED, use IPmode=:unconstr.")
		end
		pFor = Array{EqconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
			pFor[j] = EqconstrainedOEDParam(OEDparam[j],L,Ce,ce[:,ids.==j],zeros(0,0),[])
		end
	elseif IPmode==:nonnegconstr
		I = speye(size(Xtrue,1))
		pFor = Array{IneqconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
				pFor[j] = getIneqconstrainedOEDParam(OEDparam[j],L,Ce,ce[:,ids.==j],I,0*OEDparam[j].Ftrue,
				          ftolIPQP=1e-3,mutolIPQP=1e-3,outIPQP=1)
		end

	elseif IPmode==:boxconstr
		if isempty(Ci) || isempty(ci)
			error("inequality constraints must be provided. To run equality constrained OED, use IPmode=:eqconstr.")
		end
		pFor = Array{IneqconstrainedOEDParam}(nworkers())

		for j=1:nworkers()
		pFor[j] = 	getIneqconstrainedOEDParam(OEDparam,L,Ce,ce[:,ids.==j],Ci,ci[:,ids.==j],
		ftolIPQP=1e-3,mutolIPQP=1e-3,outIPQP=0)
		end
	else
		error("Invalid IPmode = $IPmode")
	end
	display(pFor[1])

	F    = zeros(length(lambdas),nrep)
	dF   = zeros(length(lambdas),nrep)
	Dc    = zeros(size(Xtrue,1),size(Xtrue,2),length(lambdas))

	for j=1:nrep
	for k=1:length(lambdas)
	    println("lambda=$(lambdas[k])")
			if nworkers()>1
				pMis  = Array{RemoteChannel}(nworkers())
				for j2=1:nworkers()
					pFor[j2].L = sqrt(lambdas[k])*speye(n*n)
					pMis[j2] = initRemoteChannel(getMisfitParam,workers()[j2], pFor[j2],ones(size(OEDparam[j2].Ftrue)),Xtrue[:,ids.==j2],SSDFun)
				end
		        Dcr,F[k,j],dFk,d2F,pMis,tMis = computeMisfit(wOpt,pMis);
				dF[k,j] = norm(dFk)
				for j2=1:nworkers()
					Dc[:,ids.==j2,k] = fetch(Dcr[j2])
				end
			else
				pFor[1].L = sqrt(lambdas[k])*speye(n*n)
				pMis  = getMisfitParam(pFor[1],ones(size(Xtrue)),Xtrue,SSDFun);
		        Dc[:,:,k],F[k,j],dFk,d2F,pMis,tMis = computeMisfit(wOpt,pMis);
				dF[k,j] = norm(dFk)
			end

	end
end
	outfile2 = @sprintf "%s-lambdas-%s.mat" resfile IPmode
	matwrite(outfile2,Dict("F"=>F,"lambdas"=>lambdas,"Dc"=>Dc,"dF"=>dF,"wOpt"=>wOpt))
end


text:
"""This module implements a LALR(1) Parser
"""
# Author: Erez Shinan (2017)
# Email : erezshin@gmail.com

# Adapted to Julia by James Hester (2019)
# Email : james.r.hester@gmail.com

export feed_token!

struct ParseConf
    parse_table::ParseTable
    start_state::Int64
    end_state::Int64
    states::Array{Dict{String,Tuple{Symbol,Any}},1}
    callbacks::Dict{Rule,Function}
    start::String
end

struct _Parser
    parse_table::ParseTable
    callbacks::Dict{Rule,Function}
    debug::Bool
end

_Parser(parse_table,callbacks;debug=false) = begin
    _Parser(parse_table,callbacks,debug)
end

ParseConf(parse_table,callbacks,start) = begin
    ParseConf(parse_table,parse_table.start_states[start],
              parse_table.end_states[start],
              parse_table.states,
              callbacks,
              start)
end

struct LALRParser
    _parse_table::ParseTable
    parser_conf::ParserConf
    parser::_Parser
    debug::Bool
end

LALRParser(parser_conf;debug=false) = begin
    analysis = LALR_Analyzer(parser_conf,debug=debug)
    compute_lalr!(analysis)
    callbacks = parser_conf.callbacks
    parser = _Parser(analysis.parse_table,callbacks)
    LALRParser(analysis.parse_table,parser_conf,parser,debug)
end

parse(l::LALRParser,args...;kwargs...) = begin
    #println("Parsing with $args, $kwargs")
    parse(l.parser,args...,kwargs...)
end

struct ParserState
    parse_conf::ParseConf
    state_stack::Stack{Int64}
    value_stack::Stack{Any}
end

ParserState(parse_conf;state_stack=nothing,value_stack=nothing) = begin
    if state_stack === nothing
        state_stack = Stack{Int64}()
        push!(state_stack,parse_conf.start_state)
    end

    ParserState(parse_conf, state_stack,
                value_stack === nothing ? Stack{Any}() : value_stack
                )
end

position(ps::ParserState) = first(ps.state_stack)

feed_token!(ps::ParserState,token;is_end=false) = begin
    state_stack = ps.state_stack
    value_stack = ps.value_stack
    states = ps.parse_conf.states
    end_state = ps.parse_conf.end_state
    callbacks = ps.parse_conf.callbacks
    while true
        arg = nothing  #create outside scope of try
        action = nothing # ditto
        state = first(state_stack)
        #println("State $state, current token $(token.type_) ($token)")
        try
            #println("Possibles: $(keys(states[state]))\n")
            action,arg = states[state][token.type_]
        catch e
            if e isa KeyError
                expected = [s for s in keys(states[state]) if all(x->isuppercase(x),s)]
                throw(UnexpectedToken(token, expected, state=state, puppet=nothing))
            else
                rethrow(e)
            end
        end
        @assert arg != end_state
        if action == :shift
            @assert !is_end
            #println("Shift!")
            push!(state_stack,arg)
            push!(value_stack,token)
            return
        end
        size = length(arg.expansion)::Int
        #println("Preparing to pop $size values")
        if size > 0
            s = grab!(value_stack,size)
            for i = 1:size
                pop!(state_stack)
            end
        else
            s = []
        end
        value = callbacks[arg](s)

        #println("Reduced, value now $value")
        _action, new_state = states[first(state_stack)][arg.origin.name]::Tuple{Symbol,Int64}
        @assert _action == :shift
        push!(state_stack,new_state)
        push!(value_stack,value)
        if is_end && first(state_stack) == end_state
            return first(value_stack)
        end
    end
end

parse(p::_Parser,lexer::LexerThread,start::String; value_stack = nothing, state_stack=nothing) = begin
    parse_conf = ParseConf(p.parse_table,p.callbacks,start)
    parser_state = ParserState(parse_conf,state_stack=state_stack,value_stack=value_stack)
    return parse_from_state(p,lexer,parser_state)
end

parse_from_state(p::_Parser,lexer,state::ParserState) = begin
    # Main LALR-parser loop
    lp = LexerParser(lexer,state)
    try
        token = nothing
        for token in lp
            #println("Seen $(token.type_): $token")
            feed_token!(state,token)
        end
        token = !(token===nothing) ? new_borrow_pos("\$END","",token) : Token("\$END","",
                                                                              pos_in_stream=0,
                                                                              line=1,
                                                                              column=1)
        return feed_token!(state,token,is_end=true)
    catch e
        if e isa UnexpectedInput
            println("Unexpected input: $e")
            try
                e.puppet = ParserPuppet(p,state,lexer)
            catch f
                if !(f isa UndefVarError) #parser puppet might not be defined
                    throw(f)
                end
            end
            rethrow(e)
        else
            if p.debug
                println("")
                println("STATE STACK DUMP")
                println("----------------")
                for (i,s) in enum(state.state_stack)
                    println("$i)$s")
                end
                println("")
            end
            rethrow(e)
        end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""
    DefsOkFalse(;
    )

"""
mutable struct DefsOkFalse <: SwaggerModel

    function DefsOkFalse(;)
        o = new()
        o
    end
end # type DefsOkFalse

const _property_map_DefsOkFalse = Dict{Symbol,Symbol}()
const _property_types_DefsOkFalse = Dict{Symbol,String}()
Base.propertynames(::Type{ DefsOkFalse }) = collect(keys(_property_map_DefsOkFalse))
Swagger.property_type(::Type{ DefsOkFalse }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DefsOkFalse[name]))}
Swagger.field_name(::Type{ DefsOkFalse }, property_name::Symbol) =  _property_map_DefsOkFalse[property_name]

function check_required(o::DefsOkFalse)
    true
end

function validate_property(::Type{ DefsOkFalse }, name::Symbol, val)
end


text:
using AtomsBase
using DFTK
using InteratomicPotentials
using StaticArrays
using Test
using Unitful
using UnitfulAtomic

@testset "DFTK -> InteratomicPotentials" begin
    functionals = [:lda_x, :lda_c_pw]
    scf_kwargs = Dict(:damping => 0.7, :tol => 1e-4)
    potential = DFTKPotential(5u"hartree", [1, 1, 1]; functionals, scf_kwargs)

    particles = [
        :Ar => [21.0, 21.0, 21.0]u"bohr",
        :Ar => [7.0, 21.0, 21.0]u"bohr",
        :Ar => [21.0, 7.0, 21.0]u"bohr",
        :Ar => [7.0, 7.0, 21.0]u"bohr",
        :Ar => [21.0, 21.0, 7.0]u"bohr",
        :Ar => [7.0, 21.0, 7.0]u"bohr",
        :Ar => [21.0, 7.0, 7.0]u"bohr",
        :Ar => [7.0, 7.0, 7.0]u"bohr"
    ]
    box = [[28.0, 0.0, 0.0], [0.0, 28.0, 0.0], [0.0, 0.0, 28.0]]u"bohr"
    system = attach_psp(periodic_system(particles, box); functional="lda")

    eandf = energy_and_force(system, potential)
    @test eandf.e isa Unitful.Energy
    @test eandf.f isa AbstractVector{<:SVector{3,<:Unitful.Force}}

    @test haskey(potential.scf_kwargs, :ψ)
    @test haskey(potential.scf_kwargs, :ρ)
end


text:
# lazy outer products of Fields, which comes up alot in automatic differentiation
*(x::Field, y::Adjoint{<:Any, <:Field}) = OuterProdOp(x, y.parent)

# this does basis promotion, unlike Zygote's default for AbstractArrays
Zygote.accum(a::Field, b::Field) = a+b
# this may create a LazyBinaryOp, unlike Zygote's
Zygote.accum(a::LinOp, b::LinOp) = a+b

## Fields

# we have to define several custom adjoints which account for the automatic
# basis conversion which CMBLensing does. part of the reason we have to do this
# explicilty is because of:
# https://discourse.julialang.org/t/how-to-deal-with-zygote-sometimes-pirating-its-own-adjoints-with-worse-ones

# ℝᴺ -> ℝ¹
@adjoint sum(f::Field{B}) where {B} = sum(f), Δ -> (Δ*one(f),)
@adjoint norm(f::Field) = Zygote.pullback(f->sqrt(dot(f,f)), f)
@adjoint dot(f::Field{B1}, g::Field{B2}) where {B1,B2} = dot(f,g), Δ -> (Δ*B1(g), Δ*B2(f))
@adjoint *(f::Adjoint{<:Any,<:Field}, g::Field) = Zygote.pullback((f,g)->dot(f',g),f,g)
# ℝᴺˣᴺ -> ℝ¹
@adjoint logdet(L::ParamDependentOp, θ) = Zygote._pullback(θ->logdet(L(;θ...)), θ) # dont need to take out offset here like in
@adjoint logdet(L::DiagOp) = logdet(L), Δ -> (Δ * pinv(L)',) # this was wrong in Zygote before

# basis conversion
@adjoint (::Type{B})(f::Field{B′}) where {B<:Basis, B′} = B(f), Δ -> (B′(Δ),)

# algebra
@adjoint +(f::Field{B1}, g::Field{B2}) where {B1,B2} = f+g, Δ -> (B1(Δ), B2(Δ))
@adjoint *(a::Real, f::Field{B}) where {B} = a*f, Δ -> (f'*Δ, B(Δ*a))
@adjoint *(a::Real, L::DiagOp) = a*L, Δ -> (tr(L'*Δ), a*Δ) # need to use trace here since it unfolds the diagonal

# operators
@adjoint *(D::DiagOp{<:Field{B}}, v::Field{B′}) where {B,B′} = D*v, Δ->(B(Δ)*B(v)', B′(D'*Δ))
@adjoint \(D::DiagOp{<:Field{B}}, v::Field{B′}) where {B,B′} = begin
    z = D \ v
    function back(Δ)
        v̄ = D' \ Δ
        -B(v̄)*B(z)', B′(v̄)
    end
    z, back
end
@adjoint *(∇::DiagOp{<:∇diag}, f::Field{B}) where {B} = ∇*f, Δ->(nothing, B(∇'*Δ))
# this makes it so we only have to define adjoints for L*f, and the f'*L adjoint just uses that
@adjoint *(f::Adjoint{<:Any,<:Field}, D::Diagonal) = Zygote.pullback((f,D)->(D'*f')', f, D)

## FieldVectors

# following two definitions are almost definitely not totally right w.r.t
# putting stuff in the correct basis, although they're working for everything
# I've needed thus far

@adjoint function *(x::FieldOrOpRowVector, y::FieldVector)
    z = x * y
    # when x is a vector of Fields
    back(Δ) = (Δ * y', x' * Δ)
    # when x is a vector of Diagonals. in this case, Δ * basis(Δ)(y)' create an
    # OuterProdOp in the same basis as the Diagonals in x
    back(Δ::Field{B}) where {B} = (Δ * basis(Δ)(y)'), (x' * Δ)
    z, back
end

@adjoint function *(A::FieldOrOpMatrix, x::FieldOrOpVector)
    z = A * x
    back(Δ::FieldVector) = (B=basis(eltype(z)); (B(Δ) * B(x)', A' * Δ))
    back(Δ::FieldOrOpVector) = (Δ * x', A' * Δ) # not sure if anything special needed here?
    back(Δ) = (Δ * x', A' * Δ)
    z, back
end

# don't know why Zygote's default adjoint for this breaks in various ways but this is simple enough
@adjoint +(I::UniformScaling, L::Union{LinOp, FieldOrOpMatrix}) = I+L, Δ->(nothing, Δ)

# Zygote/lib/array.jl:311 would suggest this should be:
#    M⁻¹, Δ->(-M⁻¹' * Δ * M⁻¹' + (- M * M⁻¹ * Δ' * M⁻¹ * M⁻¹' + Δ' * M⁻¹ * M⁻¹') + (M⁻¹' * M⁻¹ * Δ' - M⁻¹' * M⁻¹ * Δ' * M⁻¹ * M),)
# I haven't derived their version, but numerically the one gives the right answer where as their doesn't...
@adjoint function pinv(L::Union{LinOp, FieldOrOpMatrix})
    L⁻¹ = pinv(L)
    L⁻¹, Δ->(-L⁻¹' * Δ * L⁻¹',)
end

@adjoint sqrt(L::DiagOp) = (z=sqrt(L);), Δ -> ((pinv(z)/2)'*Δ,)

# some stuff which arguably belongs in Zygote or ChainRules
# see also: https://github.com/FluxML/Zygote.jl/issues/316

@adjoint broadcasted(::typeof(\), x ::Numeric, y::Numeric) =
    broadcast(\, x, y), Δ -> (nothing, unbroadcast(x, @. -Δ*y/x^2), unbroadcast(y, @. Δ/x))

@adjoint (::Type{SA})(tup) where {SA<:SArray} = SA(tup), Δ->(tuple(Δ...),)

# workaround for https://github.com/FluxML/Zygote.jl/issues/686
if versionof(Zygote) > v"0.4.15"
    Zygote._zero(xs::StaticArray, T) = SizedArray{Tuple{size(xs)...},Union{T,Nothing}}(map(_->nothing, xs))
end

# functions with no gradient which Zygote would otherwise fail on

@nograd fieldinfo

# finite difference Hessian using Zygote gradients
function hessian(f, xs::Vector; ε=1f-3)
    hcat(finite_difference(xs->vcat(gradient(f,xs)[1]...),xs,ε=ε)...)
end

function finite_difference(f, xs::Vector; ε=1f-3, progress=false)
    @showprogress (progress ? 1 : Inf) map(1:length(xs)) do i
        xs₊ = copy(xs); xs₊[i] += ε
        xs₋ = copy(xs); xs₋[i] -= ε
        (f(xs₊) .- f(xs₋)) ./ (2ε)
    end
end


text:
using VegaLite
using DataFrames

## histograms by group

df= DataFrame(group=rand(0:1, 200))
df[:x] = df[:group]*2 + randn(size(df,1))

data_values(df) +
    mark_bar() +
    encoding_x_quant(:x; bin=Dict(:maxbins=>15)) +
    encoding_y_quant(:*, aggregate="count") +
    encoding_row_nominal(:group) +
    encoding_color_nominal(:group)


text:
export TimelessInstantModel

"""
    abstract type TimelessInstantModel{I,O} <: InstantModel{I,O}

An InstantModel in which the sfunc made does not depend on time.
Must implement a version of `make_initial` that does not take the current time as argument.
Note that `make_initial` can be defined to take keyword arguments, so the sfunc created
need not be exactly the same every time.
"""

abstract type TimelessInstantModel{I,O} <: InstantModel{I,O} end

make_initial(m::TimelessInstantModel,t;keys...) = make_initial(m;keys...)


text:
push!(LOAD_PATH,"..")
using Documenter, PlanerTransformer

const PAGES = [
  "Home" => "index.md",
  "Installation" => "install.md",
  "Introduction" => "introduction.md",
  "Design API" => "design_api.md",
  "Data Entry API" => "data_entry_api.md",
  "Internal API" => "internal_api.md",
  "API index" => "APIindex.md",
]

makedocs(
  modules = [PlanerTransformer],
  doctest   = "doctest" in ARGS,
  linkcheck = "linkcheck" in ARGS,
  format    = "pdf" in ARGS ? :latex : :html,
  build = "site",
  sitename = "PlanerTransformer",
  authors = "Chris Stook",
  pages = PAGES,
  html_prettyurls = "deploy" in ARGS,
)

if "deploy" in ARGS
  fake_travis = "C:/Users/Chris/fake_travis_PlanerTransformer.jl"
  if isfile(fake_travis)
    include(fake_travis)
  end
  deploydocs(
    repo = "github.com/cstook/PlanerTransformer.jl.git",
    target = "site",
    branch = "gh-pages",
    latest = "master",
    osname = "linux",
    julia  = "1.0",
    deps = nothing,
    make = nothing,
  )
end


text:
function  ueberpruefePara(PX1,PX2,theta,a,b,ParameterdarstellungX1,ParameterdarstellungX2)

t = 0.6221;
Z = 3;
S = 1;
x = 0.7291;

   if ParameterdarstellungX1(t,Z,S) == PX1[Z,S]+t*cosd(theta[S])
        if ParameterdarstellungX2(t,Z,S) == PX2[Z,S]+t*sind(theta[S])
            println("Eure Parameterdarstellung für die Strahlen ist korrekt!");
        elseif ParameterdarstellungX2(t,Z,S) == PX2[Z,S]+t*cosd(theta[S])
            println("Der Teil für die X1 Koordinate ist korrekt. Schaut euch für die X2 Koordinaten noch einmal an, wie der Kosinus und Sinus definiert sind.")
        else
            println("Der Teil für die X1 Koordinate ist korrekt. Überprüft erneut eure Formel für die X2 Koordinate.")
        end
    elseif ParameterdarstellungX1(t,Z,S) == PX1[Z,S]+t*sind(theta[S])
        if ParameterdarstellungX2(t,Z,S) == PX2[Z,S]+t*sind(theta[S])
            println("Der Teil für die X2 Koordinate ist korrekt. Schaut euch für die X1 Koordinaten noch einmal an, wie der Kosinus und Sinus definiert sind.")
        else
            println("Schaut euch erneut an, wie der Sinus und Kosinus definiert sind und überprüft erneut eure Formeln.")
        end
    else
        println("Schaut euch erneut an, wie der Sinus und Kosinus definiert sind und überprüft erneut eure Formeln.")
   end
end


text:
# 118 Bus Model based on the IEEE 118 bus case data in http://motor.ece.iit.edu/data/JEAS_IEEE118.doc
# Wind data and modifications partly based on Roald et al's paper "Optimal power flow with wind power control and limited expected risk of overloads"

# OPTIONS for caseID:
# wind_penetration_level_1: 25% wind penetration
# wind_penetration_level_2: 50% wind penetration
# wind_penetration_level_3: 75% wind penetration
# wind_penetration_level_4: 100% wind penetration
# wind_penetration_level_5: 125% wind penetration

const caseID = "wind_penetration_level_1"

# include basic model parameters
include("params.jl")

#*========== BASIC TOPOLOGY ===========

# number of buses, generators (regular and wind), loads
const numBuses = 118
const numGen = 79
const numRegGen = 54
const numWindGen = 25
const numLoads = 91

# map regular generators to buses
# from column 2 of Table 1 in the reference
const regGenToBus = [4,6,8,10,12,15,18,19,24,25,26,27,31,32,34,36,40,42,46,49,54,55,56,59,61,62,65,66,69,70,72,73,74,76,77,80,82,85,87,89,90,91,92,99,100,103,104,105,107,110,111,112,113,116]

# map wind generators to buses
# from Roald et al's spreadsheet
const windGenToBus = [3,6,9,18,25,29,33,35,41,43,46,56,60,61,70,78,83,85,92,93,96,107,113,117,118]

# map loads to buses
# from column 1 of Table 6 in the reference
const loadToBus = [1,2,3,4,6,7,11,12,13,14,15,16,17,18,19,20,21,22,23,27,28,29,31,32,33,34,35,36,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,62,66,67,70,74,75,76,77,78,79,80,82,83,84,85,86,88,90,92,93,94,95,96,97,98,100,101,102,103,104,105,106,107,108,109,110,112,114,115,117,118]

# lines between buses
# from columns 2 and 3 of Table 3 in the reference
const lines = ((1,2), (1,3), (4,5), (3,5), (5,6), (6,7), (8,9), (8,5), (9,10), (4,11), (5,11), (11,12), (2,12), (3,12), (7,12), (11,13), (12,14), (13,15), (14,15), (12,16), (15,17), (16,17), (17,18), (18,19), (19,20), (15,19), (20,21), (21,22), (22,23), (23,24), (23,25), (26,25), (25,27), (27,28), (28,29), (30,17), (8,30), (26,30), (17,31), (29,31), (23,32), (31,32), (27,32), (15,33), (19,34), (35,36), (35,37), (33,37), (34,36), (34,37), (38,37), (37,39), (37,40), (30,38), (39,40), (40,41), (40,42), (41,42), (43,44), (34,43), (44,45), (45,46), (46,47), (46,48), (47,49), (42,49), (42,49), (45,49), (48,49), (49,50), (49,51), (51,52), (52,53), (53,54), (49,54), (49,54), (54,55), (54,56), (55,56), (56,57), (50,57), (56,58), (51,58), (54,59), (56,59), (56,59), (55,59), (59,60), (59,61), (60,61), (60,62), (61,62), (63,59), (63,64), (64,61), (38,65), (64,65), (49,66), (49,66), (62,66), (62,67), (65,66), (66,67), (65,68), (47,69), (49,69), (68,69), (69,70), (24,70), (70,71), (24,72), (71,72), (71,73), (70,74), (70,75), (69,75), (74,75), (76,77), (69,77), (75,77), (77,78), (78,79), (77,80), (77,80), (79,80), (68,81), (81,80), (77,82), (82,83), (83,84), (83,85), (84,85), (85,86), (86,87), (85,88), (85,89), (88,89), (89,90), (89,90), (90,91), (89,92), (89,92), (91,92), (92,93), (92,94), (93,94), (94,95), (80,96), (82,96), (94,96), (80,97), (80,98), (80,99), (92,100), (94,100), (95,96), (96,97), (98,100), (99,100), (100,101), (92,102), (101,102), (100,103), (100,104), (103,104), (103,105), (100,106), (104,105), (105,106), (105,107), (105,108), (106,107), (108,109), (103,110), (109,110), (110,111), (110,112), (17,113), (32,113), (32,114), (27,115), (114,115), (68,116), (12,117), (75,118), (76,118))

#*========== GENERATOR BOUNDS ===========

# lower and upper power bounds for regular generators (pmin and pmax)
const regGenMinPow = zeros(Float64,numRegGen)
# from column 6 of Table 1 in the reference
# change upper limit of generator#49 to be large enough to ensure relatively complete recourse holds
const regGenMaxPow = [30.0,30.0,30.0,300.0,300.0,30.0,100.0,30.0,30.0,300.0,350.0,30.0,30.0,100.0,30.0,100.0,30.0,30.0,100.0,250.0,250.0,100.0,100.0,200.0,200.0,100.0,420.0,420.0,300.0,80.0,30.0,30.0,20.0,100.0,100.0,300.0,100.0,30.0,300.0,200.0,20.0,50.0,300.0,300.0,300.0,20.0,100.0,100.0,10000.0,50.0,100.0,100.0,100.0,50.0]

# lower power bounds for wind generators
const windGenMinPow = zeros(Float64,numWindGen)

#*========== COST COEFFICIENTS ===========

# regular generator costs including cost of reserves
# only use linear cost coefficients
# from column 4 of Table 1 in the reference
const regGenCost = [26.2438,26.2438,26.2438,12.8875,12.8875,26.2438,17.8200,26.2438,26.2438,12.8875,10.7600,26.2438,26.2438,17.8200,26.2438,17.8200,26.2438,26.2438,17.8200,12.3299,12.3299,17.8200,17.8200,13.2900,13.2900,17.8200,8.3391,8.3391,12.8875,15.4708,26.2438,26.2438,37.6968,17.8200,17.8200,12.8875,17.8200,26.2438,10.7600,12.8875,37.6968,22.9423,12.8875,12.8875,12.8875,37.6968,17.8200,17.8200,37.6968,22.9423,17.8200,17.8200,17.8200,22.9423]

# scale reserves cost by suitable factor
const regGenPosResCost = reservesCostFactor*regGenCost
const regGenNegResCost = 1.0*regGenPosResCost

# wind generator costs including cost of reserves (assume zero marginal costs)
const windGenCost = zeros(Float64,numWindGen)

# set cost of reserves for each wind generator to be the minimum reserves cost of all regular generators
const minRegGenPosResCost = minimum(regGenPosResCost)
const minRegGenNegResCost = minimum(regGenNegResCost)

const windGenPosResCost = minRegGenPosResCost*windReservesCostFactor*ones(Float64,numWindGen)
const windGenNegResCost = minRegGenNegResCost*windReservesCostFactor*ones(Float64,numWindGen)

#*========== RESERVES INFORMATION ===========

# which regular generators are ready to provide reserves?
const regGenRes = trues(numRegGen)

# which regular generators are ALWAYS READY to provide reserves?
# make sure that the generator with large enough max generation limit always provides reserves to ensure relatively complete recourse holds
# we assume that all generators provide guaranteed reserves with a minimum participation factor min_alpha so that numerical issues don't crop up while running the stochastic approximation method
const regGenGuarRes = trues(numRegGen)

# which wind generators are ready to provide reserves?
const windGenRes = trues(numWindGen)

# which wind generators are ALWAYS READY to provide reserves?
const windGenGuarRes = trues(numWindGen)

# minimum participation factors for generators ALWAYS READY to provide reserves
# reset this min participation factor to not restrict the space of solutions too much, if necessary
const min_alpha = 0.001
const numGenProvRes = sum(windGenGuarRes) + sum(regGenGuarRes)
if(numGenProvRes > 10)
	min_alpha = 1.0/(100.0*numGenProvRes)
end

#*========== LINE PARAMETERS ===========

# these line parameters are assumed to be in the same order as "lines"
# line susceptances (beta)
# from column 6 of Table 3 in the reference
const beta = 1./[0.09990,0.04240,0.00798,0.10800,0.05400,0.02080,0.03050,0.02670,0.03220,0.06880,0.06820,0.01960,0.06160,0.16000,0.03400,0.07310,0.07070,0.24440,0.19500,0.08340,0.04370,0.18010,0.05050,0.04930,0.11700,0.03940,0.08490,0.09700,0.15900,0.04920,0.08000,0.03820,0.16300,0.08550,0.09430,0.03880,0.05040,0.08600,0.15630,0.03310,0.11530,0.09850,0.07550,0.12440,0.24700,0.01020,0.04970,0.14200,0.02680,0.00940,0.03750,0.10600,0.16800,0.05400,0.06050,0.04870,0.18300,0.13500,0.24540,0.16810,0.09010,0.13560,0.12700,0.18900,0.06250,0.32300,0.32300,0.18600,0.05050,0.07520,0.13700,0.05880,0.16350,0.12200,0.28900,0.29100,0.07070,0.00955,0.01510,0.09660,0.13400,0.09660,0.07190,0.22930,0.25100,0.23900,0.21580,0.14500,0.15000,0.01350,0.05610,0.03760,0.03860,0.02000,0.02680,0.09860,0.03020,0.09190,0.09190,0.21800,0.11700,0.03700,0.10150,0.01600,0.27780,0.32400,0.03700,0.12700,0.41150,0.03550,0.19600,0.18000,0.04540,0.13230,0.14100,0.12200,0.04060,0.14800,0.10100,0.19990,0.01240,0.02440,0.04850,0.10500,0.07040,0.02020,0.03700,0.08530,0.03665,0.13200,0.14800,0.06410,0.12300,0.20740,0.10200,0.17300,0.07120,0.18800,0.09970,0.08360,0.05050,0.15810,0.12720,0.08480,0.15800,0.07320,0.04340,0.18200,0.05300,0.08690,0.09340,0.10800,0.20600,0.29500,0.05800,0.05470,0.08850,0.17900,0.08130,0.12620,0.05590,0.11200,0.05250,0.20400,0.15840,0.16250,0.22900,0.03780,0.05470,0.18300,0.07030,0.18300,0.02880,0.18130,0.07620,0.07550,0.06400,0.03010,0.20300,0.06120,0.07410,0.01040,0.00405,0.14000,0.04810,0.05440]

# line flow limits
const lineFlowFactor = 0.75
# from column 8 of Table 3 in the reference
const flow_max = [175.0,175.0,500.0,175.0,175.0,175.0,500.0,500.0,500.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,500.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,500.0,500.0,500.0,175.0,175.0,500.0,175.0,500.0,175.0,175.0,140.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,500.0,500.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,500.0,175.0,175.0,500.0,500.0,500.0,500.0,500.0,500.0,500.0,175.0,175.0,500.0,175.0,500.0,175.0,175.0,500.0,500.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,500.0,175.0,175.0,175.0,175.0,175.0,175.0,500.0,500.0,175.0,500.0,500.0,200.0,200.0,175.0,175.0,175.0,500.0,500.0,175.0,175.0,500.0,500.0,500.0,175.0,500.0,500.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,200.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,500.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,175.0,500.0,175.0,175.0,175.0,500.0,175.0,175.0,175.0]*lineFlowFactor

#*========== LOAD PARAMETERS ===========

# power demands/loads
const demandFactor = 1.5

# from column 2 of Table 6 in the reference
const meanDemand = [54.14,21.23,41.40,31.85,55.20,20.17,74.31,49.89,36.09,14.86,95.54,26.54,11.68,63.69,47.77,19.11,14.86,10.62,7.43,65.82,18.05,25.48,45.65,62.63,24.42,62.63,35.03,32.91,27.00,20.00,37.00,37.00,18.00,16.00,53.00,28.00,34.00,20.00,87.00,17.00,17.00,18.00,23.00,113.00,63.00,84.00,12.00,12.00,277.00,78.00,77.00,39.00,28.00,66.00,68.00,47.00,68.00,61.00,71.00,39.00,130.00,54.00,20.00,11.00,24.00,21.00,48.00,78.00,65.00,12.00,30.00,42.00,38.00,15.00,34.00,37.00,22.00,5.00,23.00,38.00,31.00,43.00,28.00,2.00,8.00,39.00,25.00,8.49,23.35,21.23,33.00]*demandFactor

#*========== WIND PARAMETERS ===========

# Wind prediction data
const windFactor = 1.0

# from Roald et al's spreadsheet
const windPredMean = [109.0,70.0,200.0,147.0,102.0,148.0,70.0,60.0,96.0,80.0,105.0,113.0,80.0,84.0,59.0,146.0,69.0,250.0,118.0,89.0,127.0,76.0,72.0,129.0,129.0]

const sumMeanWind = sum(windPredMean)
const sumMeanDemand = sum(meanDemand)

# set average wind output depending on the wind penetration level
if(caseID == "wind_penetration_level_1")
	windFactor = 0.25*sumMeanDemand/sumMeanWind
elseif(caseID == "wind_penetration_level_2")
	windFactor = 0.5*sumMeanDemand/sumMeanWind
elseif(caseID == "wind_penetration_level_3")
	windFactor = 0.75*sumMeanDemand/sumMeanWind
elseif(caseID == "wind_penetration_level_4")
	windFactor = 1.0*sumMeanDemand/sumMeanWind
elseif(caseID == "wind_penetration_level_5")
	windFactor = 1.25*sumMeanDemand/sumMeanWind
end

windPredMean *= windFactor

#*========== UNCERTAIN PARAMETERS ===========

# standard deviation of loads/demands. keep this small enough to be realistic
const demandStdevFactor = 0.05
const demandStdev = demandStdevFactor*meanDemand

# randomly generated correlation matrix for the demands
include("118bus_demandCorrMat.jl")

# standard deviation of wind power output
const windPredStdevFactor = 0.1
const windPredStdev = windPredStdevFactor*windPredMean

# randomly generated correlation matrix for the wind power outputs
include("118bus_windCorrMat.jl")

#*========== SCENARIO GENERATION ===========

# store the Cholesky factorization of the demand covariance matrix for numerical efficiency
const demandCovMat = Diagonal(demandStdev)*demandCorrMat*Diagonal(demandStdev)
const demandCovMat_chol = (cholfact(Hermitian(demandCovMat)))[:L]

# store the Cholesky factorization of the wind covariance matrix for numerical efficiency
const windCovMat = Diagonal(windPredStdev)*windCorrMat*Diagonal(windPredStdev)
const windCovMat_chol = (cholfact(Hermitian(windCovMat)))[:L]

# tailored implementation of multivariate normal distribution for demand scenarios
function generateDemandScenarios(numSamples::Int64)

	omega_d = zeros(Float64,numLoads,numSamples)

	omega_d_tmp = rand(Normal(0.0,1.0),numLoads,numSamples)
	for samp = 1:numSamples
		omega_d[:,samp] = demandCovMat_chol*omega_d_tmp[:,samp] + meanDemand
	end

	return omega_d
end

# tailored implementation of multivariate normal distribution for wind scenarios
function generateWindScenarios(numSamples::Int64)

	omega_w = zeros(Float64,numWindGen,numSamples)

	omega_w_tmp = rand(Normal(0.0,1.0),numWindGen,numSamples)
	for samp = 1:numSamples
		omega_w[:,samp] = windCovMat_chol*omega_w_tmp[:,samp] + windPredMean
	end

	return omega_w
end

# Generate demand and wind scenarios
function generateScenarios(numSamples::Int64)

	omega_d = generateDemandScenarios(numSamples)
	omega_w = generateWindScenarios(numSamples)

	return omega_d, omega_w
end

max_stars_repo_name: rohitkannan/DCOPF-reserve-saturation

max_issues_repo_name: rohitkannan/DCOPF-reserve-saturation

max_forks_repo_name: rohitkannan/DCOPF-reserve-saturation


text:
"""
    C_iio_context_clone(context)

Duplicate a pre-existing IIO context.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure

# Returns
- On success, A pointer to an iio_context structure
- On failure, throws an error if the assertions are enabled, or NULL otherwise.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga1815e7c39b9a69aa11cf948b0433df01)
"""
function C_iio_context_clone(context::Ptr{iio_context})
    @assert_null_pointer clone = ccall(
        (:iio_context_clone, libIIO),
        Ptr{iio_context}, (Ptr{iio_context},),
        context
    );
    return clone;
end

"""
    C_iio_context_destroy(context)

Destroy the given context.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure

# NOTE
After that function, the iio_context pointer shall be invalid.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga75de8dae515c539818e52b408830d3ba)
"""
function C_iio_context_destroy(context::Ptr{iio_context})
    ccall(
        (:iio_context_destroy, libIIO),
        Cvoid, (Ptr{iio_context},),
        context
    );
end

"""
    C_iio_context_find_device(context, name)

Try to find a device structure by its name of ID.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure
- `name::String` : A NULL-terminated string corresponding to the name or the ID of the device to search for

# Returns
- On success, a pointer to an iio_device structure
- If the name or ID does not correspond to any known device, an error is thrown if the assertions are enabled, or NULL otherwise.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#gade1dadfb5bc3c3b236add67f803c50c3)
"""
function C_iio_context_find_device(context::Ptr{iio_context}, name::String)
    @assert_null_pointer device = ccall(
        (:iio_context_find_device, libIIO),
        Ptr{iio_device}, (Ptr{iio_context}, Cstring),
        context, name
    );
    return device;
end

"""
    C_iio_context_get_attr(context, index)

Retrieve the name and value of a context-specific attribute.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure
- `index::UInt32` : The index corresponding to the attribute

# Returns
- On success, `(0, name::String, value::String)` is returned.
- On error, `(errno, "", "")` is returned, where errno is a negative code.

Introduced in version 0.9.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga477dfddaefe0acda401f600247e13fc7)
"""
function C_iio_context_get_attr(context::Ptr{iio_context}, index::UInt32)
    name, value = Ref{Cstring}(), Ref{Cstring}();
    ret = ccall(
        (:iio_context_get_attr, libIIO),
        Cint, (Ptr{iio_context}, Cuint, Ptr{Cstring}, Ptr{Cstring}),
        context, index, name, value
    );
    return (ret == 0) ? (ret, Base.unsafe_string(name[]), Base.unsafe_string(value[])) : (ret, "", "");
end

"""
    C_iio_context_get_attr_value(context, name)

Retrieve the value of a context-specific attribute.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure
- `name::String` : The name of the context attribute to read

Returns
- On success, a NULL-terminated string.
- If the name does not correspond to any attribute and the assertions are enabled, throws an error.
- If the name does not correspond to any attribute and the assertions are disabled, returns an empty string.

Introduced in version 0.9.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga6394d108d425e4a6ed28d00c0e93d6ed)
"""
function C_iio_context_get_attr_value(context::Ptr{iio_context}, name::String)
    @assert_Cstring value = ccall(
        (:iio_context_get_attr_value, libIIO),
        Cstring, (Ptr{iio_context}, Cstring),
        context, name
    );
    return value != C_NULL ? Base.unsafe_string(value) : "";
end

"""
    C_iio_context_get_attrs_count(context)

Get the number of context-specific attributes.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure

# Returns
- The number of context-specific attributes

Introduced in version 0.9.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga91e0c4ed91d760b411d4cbea28c993da)
"""
function C_iio_context_get_attrs_count(context::Ptr{iio_context})
    return ccall(
        (:iio_context_get_attrs_count, libIIO),
        Cuint, (Ptr{iio_context},),
        context
    );
end

"""
    C_iio_context_get_description(context)

Get a description of the given context.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure

# Returns
- A NULL-terminated string

# NOTE
The returned string will contain human-readable information about the current context.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga5591da0927887e88be4ef7d670cb60a9)
"""
function C_iio_context_get_description(context::Ptr{iio_context})
    return Base.unsafe_string(ccall(
        (:iio_context_get_description, libIIO),
        Cstring, (Ptr{iio_context},),
        context
    ));
end

"""
    C_iio_context_get_device(context, index)

Get the device present at the given index.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure
- `index::UInt32` : The index corresponding to the device

# Returns
- On success, a pointer to an iio_device structure
- If the index is invalid and the assertions are enabled, an error is thrown.
- If the index is invalid and the assertions are disabled, NULL is returned.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga3f2813ff34bf96c7c85dd05909f1c709)
"""
function C_iio_context_get_device(context::Ptr{iio_context}, index::UInt32)
    @assert_null_pointer device = ccall(
        (:iio_context_get_device, libIIO),
        Ptr{iio_device}, (Ptr{iio_context}, Cuint),
        context, index
    );
    return device;
end

"""
    C_iio_context_get_devices_count(context)

Enumerate the devices found in the given context.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure

# Returns
- The number of devices found

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#gab4fc2a93fd5824f3c9e06aa81e8097d1)
"""
function C_iio_context_get_devices_count(context::Ptr{iio_context})
    return ccall(
        (:iio_context_get_devices_count, libIIO),
        Cuint, (Ptr{iio_context},),
        context
    );
end

"""
    C_iio_context_get_name(context)

Get the name of the given context.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure

# Returns
- A NULL-terminated string

# NOTE
The returned string will be local, xml or network when the context has been created with the local, xml and network backends respectively.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#gafed8e036873ad6f70c3db92c7136ad31)
"""
function C_iio_context_get_name(context::Ptr{iio_context})
    return Base.unsafe_string(ccall(
        (:iio_context_get_name, libIIO),
        Cstring, (Ptr{iio_context},),
        context
    ));
end

"""
    C_iio_context_get_version(context)

Get the version of the backend in use.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure

# Returns
- `ret::Int` : 0 if no errors, negative error code otherwise
- `major::Int` : The major version
- `minor::Int` : The minor version
- `git_tag::String` : The git tag

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga342bf90d946e7ed3815372db22c4d3a6)
"""
function C_iio_context_get_version(context::Ptr{iio_context})
    major, minor, git_tag = Ref{Cuint}(0), Ref{Cuint}(0), zeros(UInt8, 8);
    ret = ccall(
        (:iio_context_get_version, libIIO),
        Cint, (Ptr{iio_context}, Ref{Cuint}, Ref{Cuint}, Ptr{Cchar}),
        context, major, minor, pointer(git_tag)
    );
    return ret, Int(major[]), Int(minor[]), toString(git_tag);
end

"""
    C_iio_context_get_xml(context)

Obtain a XML representation of the given context.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure

# Returns
- A NULL-terminated string

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga86554706f17faf41e61e3295dc68a70c)
"""
function C_iio_context_get_xml(context::Ptr{iio_context})
    return Base.unsafe_string(ccall(
        (:iio_context_get_xml, libIIO),
        Cstring, (Ptr{iio_context},),
        context
    ));
end

"""
    C_iio_context_set_timeout(context, timeout_ms)

Set a timeout for I/O operations.

# Parameters
- `context::Ptr{iio_context}` : A pointer to an iio_context structure
- `timeout_ms::UInt32` : A positive integer representing the time in milliseconds after which a timeout occurs. A value of 0 is used to specify that no timeout should occur.

# Returns
- On success, 0 is returned
- On error, a negative errno code is returned

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#gaba3f4c4f9f885f41a6c0b9ac79b7f28d)
"""
function C_iio_context_set_timeout(context::Ptr{iio_context}, timeout_ms::UInt32)
    return ccall(
        (:iio_context_set_timeout, libIIO),
        Cint, (Ptr{iio_context}, Cuint),
        context, timeout_ms
    );
end

"""
    C_iio_create_context_from_uri(uri)

Create a context from a URI description.

# Parameters
- `uri::String` : A URI describing the context location

# Returns
- On success, a pointer to a iio_context structure
- On failure, if the assertions are enabled, an error is thrown.
- On failure, if the assertions are disabled, NULL is returned.

# NOTE
The following URIs are supported based on compile time backend support:

- Local backend, "local:":
    Does not have an address part. For example "local:"

- XML backend, "xml:"
    Requires a path to the XML file for the address part. For example "xml:/home/user/file.xml"

- Network backend, "ip:"
    Requires a hostname, IPv4, or IPv6 to connect to a specific running IIO Daemon or no address part for automatic discovery when library is compiled with ZeroConf support. For example "ip:192.168.2.1", or "ip:localhost", or "ip:" or "ip:plutosdr.local"

- USB backend, "usb:"
    When more than one usb device is attached, requires bus, address, and interface parts separated with a dot. For example "usb:3.32.5". Where there is only one USB device attached, the shorthand "usb:" can be used.

- Serial backend, "serial:" requires :
    + a port (/dev/ttyUSB0),
    + baud_rate (default 115200)
    + serial port configuration
        - data bits (5 6 7 8 9)
        - parity ('n' none, 'o' odd, 'e' even, 'm' mark, 's' space)
        - stop bits (1 2)
        - flow control ('\0' none, 'x' Xon Xoff, 'r' RTSCTS, 'd' DTRDSR)
    + For example "serial:/dev/ttyUSB0,115200" or "serial:/dev/ttyUSB0,115200,8n1"

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#gafdcee40508700fa395370b6c636e16fe)
"""
function C_iio_create_context_from_uri(uri::String)
    @assert_null_pointer context = ccall(
        (:iio_create_context_from_uri, libIIO),
        Ptr{iio_context}, (Cstring,),
        uri
    );
    return context;
end

"""
    C_iio_create_default_context()

Create a context from local or remote IIO devices.

# Returns
- On success, A pointer to an iio_context structure
- On failure, if the assertions are enabled, an error is thrown
- On failure, if the assertions are disabled, NULL is returned

# NOTE
This function will create a network context if the IIOD_REMOTE environment variable is set to the hostname where the IIOD server runs.
If set to an empty string, the server will be discovered using ZeroConf. If the environment variable is not set, a local context will be created instead.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga21076125f817a680e0c01d4a490f0416)
"""
function C_iio_create_default_context()
    @assert_null_pointer context = ccall(
        (:iio_create_default_context, libIIO),
        Ptr{iio_context}, ()
    );
    return context;
end

"""
    C_iio_create_local_context()

Create a context from local IIO devices (Linux only)

# Returns
- On success, A pointer to an iio_context structure
- On failure, if the assertions are enabled, an error is thrown.
- On failure, if the assertions are disabled, NULL is returned.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#gaf31acec2d0f9f498870cc52a1d49783e)
"""
function C_iio_create_local_context()
    @assert_null_pointer context = ccall(
        (:iio_create_local_context, libIIO),
        Ptr{iio_context}, ()
    );
    return context;
end

"""
    C_iio_create_network_context(host)

Create a context from the network.

# Parameters
- `host::String` : Hostname, IPv4 or IPv6 address where the IIO Daemon is running

# Returns
- On success, A pointer to an iio_context structure
- On failure, if the assertions are enabled, an error is thrown.
- On failure, if the assertions are disabled, NULL is returned.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga8adf2ef4d2b62aa34201469cd7049617)
"""
function C_iio_create_network_context(host::String)
    @assert_null_pointer context = ccall(
        (:iio_create_network_context, libIIO),
        Ptr{iio_context}, (Cstring,),
        host
    );
    return context;
end

"""
    C_iio_create_xml_context(xml_file)

Create a context from a XML file.

# Parameters
- `xml_file::String` : Path to the XML file to open

# Returns
- On success, A pointer to an iio_context structure
- On failure, if the assertions are enabled, an error is thrown.
- On failure, if the assertions are disabled, NULL is returned.

# NOTE
The format of the XML must comply to the one returned by `iio_context_get_xml`.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#ga9925a84e596c3003e30b1cdd2b65d029)
"""
function C_iio_create_xml_context(xml_file::String)
    @assert_null_pointer context = ccall(
        (:iio_create_xml_context, libIIO),
        Ptr{iio_context}, (Cstring,),
        xml_file
    );
    return context;
end

"""
    C_iio_create_xml_context_mem(xml, length)

Create a context from XML data in memory.

# Parameters
- `xml::String` : Pointer to the XML data in memory
- `length::UInt` : Length of the XML string in memory (excluding the final \0)

# Returns
- On success, A pointer to an iio_context structure
- On failure, if the assertions are enabled, an error is thrown.
- On failure, if the assertions are disabled, NULL is returned.

# NOTE
The format of the XML must comply to the one returned by `iio_context_get_xml`

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Context.html#gabaa848ca554af5723a44b9b7fd0ba6a3)
"""
function C_iio_create_xml_context_mem(xml::String, length::UInt)
    @assert_null_pointer context = ccall(
        (:iio_create_xml_context_mem, libIIO),
        Ptr{iio_context}, (Cstring, Csize_t),
        xml, length
    );
end


text:
using Test
using Distributed
using LinearAlgebra
using SymmetricTensors
using Cumulants
using CumulantsFeatures
using Combinatorics
using Distributed
using Random
using CumulantsUpdates
import CumulantsFeatures: reduceband, greedestep, unfoldsym, hosvdstep, greedesearchdata, mev, mormbased, hosvdapprox
import CumulantsFeatures: updatemoments

include("test_select.jl")
include("test_detectors.jl")
include("symten2mattest.jl")

addprocs(3)
@everywhere using CumulantsFeatures
@everywhere using SymmetricTensors
@everywhere using Cumulants
@everywhere import SymmetricTensors: getblock
include("paralleltests.jl")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AdditionalCapabilities <: SwaggerModel
    ultraSSDEnabled::Any # spec type: Union{ Nothing, Bool } # spec name: ultraSSDEnabled

    function AdditionalCapabilities(;ultraSSDEnabled=nothing)
        o = new()
        validate_property(AdditionalCapabilities, Symbol("ultraSSDEnabled"), ultraSSDEnabled)
        setfield!(o, Symbol("ultraSSDEnabled"), ultraSSDEnabled)
        o
    end
end # type AdditionalCapabilities

const _property_map_AdditionalCapabilities = Dict{Symbol,Symbol}(Symbol("ultraSSDEnabled")=>Symbol("ultraSSDEnabled"))
const _property_types_AdditionalCapabilities = Dict{Symbol,String}(Symbol("ultraSSDEnabled")=>"Bool")
Base.propertynames(::Type{ AdditionalCapabilities }) = collect(keys(_property_map_AdditionalCapabilities))
Swagger.property_type(::Type{ AdditionalCapabilities }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AdditionalCapabilities[name]))}
Swagger.field_name(::Type{ AdditionalCapabilities }, property_name::Symbol) =  _property_map_AdditionalCapabilities[property_name]

function check_required(o::AdditionalCapabilities)
    true
end

function validate_property(::Type{ AdditionalCapabilities }, name::Symbol, val)
end


text:
using HalfIntegers, Test

# write your own tests here
@test Half(3) == 1.5
@test Half(3) + Half(1) == 2
@test Half(2) + 1 == 2
@test 1 + Half(2) == 2
@test Half(2) + 1.0 == 2
@test 1.0 + Half(2) == 2


text:
# 3D Grid examples
# ===============
#

# ## Quadrilateral
function quadrilateral(;hx=0.25, hy=0.2, hz=0.1)
    X=collect(0:hx:1)
    Y=collect(0:hy:1)
    Z=collect(0:hz:1)
    simplexgrid(X,Y,Z)
end
# ![](quadrilateral.svg)

function mask_bedges()
    grid    = quadrilateral(hx=0.25, hy=0.25, hz=0.25)

    bedgemask!(grid, [0.0, 0.0, 0.0], [0.0, 0.0, 1.0], 1)
    bedgemask!(grid, [0.0, 0.0, 0.0], [0.0, 1.0, 0.0], 2)
    bedgemask!(grid, [0.0, 1.0, 0.0], [0.0, 1.0, 1.0], 3)
    bedgemask!(grid, [0.0, 0.0, 1.0], [0.0, 1.0, 1.0], 4)
    bedgemask!(grid, [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], 5)

    true
end


text:
#!/usr/bin/env julia

include("LVAChannels.jl")

using JSON
using SIUnits.ShortUnits
using LVAChannels

radius = 20µm/2
area = 4*pi*radius^2
current = -0.025nA

stim = Stim(20ms, 150ms, current/area)
ts, vs, m, d, h = run_lva(300ms, param=LVAParam(vrest=-65mV), stim=stim, sample_dt=0.025ms)

trace = Dict(
    :name => "membrane voltage",
    :sim => "numeric",
    :model => "test_kinlva",
    :units => "mV",
    :data => Dict(
        :time => map(t->t/ms, ts),
        Symbol("soma.mid") => map(v->v/mV, vs)
    )
)

state = Dict(
    :name => "mechanisms state",
    :sim => "numeric",
    :model => "kinlva",
    :units => "1",
    :data => Dict(
        :time => map(t->t/ms, ts),
        Symbol("m") => m,
        Symbol("d") => d,
        Symbol("h") => h
    )
)
println(JSON.json([trace, state]))

max_stars_repo_name: kabicm/arbor

max_issues_repo_name: kabicm/arbor

max_forks_repo_name: kabicm/arbor


text:
include("./src/required.jl")
include("./src/functions.jl")

# Testing the imputation

# Creating the CSV file for the results
header_list = DataFrame(Rank = [], Alpha = [], Suspected = [], Unlikely = [], Delta = []);
CSV.write("Results.csv", header_list, header = true, delim =';');

# Setting the parameters used for the simulation
α_list = [[0.0, 0.0, 0.0, 1.0], [0.0, 1.0, 1.0, 0.0], [0.0, 1.0, 1.0, 1.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]];
rank_list = [1, 2, 3, 4, 5];

# Creating the list of suspected and unlikely hosts found in the new data
suspected_newData = ["Artibeus jamaicensis","Desmodus rotundus", "Hipposideros larvatus","Hipposideros pomona", "Megaerops kusnotoi", "Myonycteris angolensis", "Myotis pequinius", "Nanonycteris veldkampii", "Nycteris macrotis", "Pipistrellus deserti", "Plecotus auritus", "Pteropus lylei", "Scotophilus heathii", "Scotophilus kuhlii"];
unlikely_newData = ["Carollia sowelli", "Hipposideros gigas", "Hipposideros lekaguli", "Macroglossus minimus", "Myotis horsfieldii", "Pipistrellus coromandra", "Tadarida teniotis"];

# Interating the simulations with all parameters
for i in 1:length(α_list)[1]
    for j in 1:length(rank_list)
        (matrix, hosts, viruses, df) = buildInteractionMatrix();
        crossValidation(0.0, calculateInitialeValues(matrix, α_list[i]), α_list[i], rank_list[j], suspected_newData, unlikely_newData);
    end
end


text:
module Model

# define the type of data that expect to use
# and the user info that require

import Base: ==

using StructTypes

export Album, User

mutable struct Album
    id::Int64 # service-managed
    userid::Int64 # service-managed
    name::String
    artist::String
    year::Int64
    timespicked::Int64 # service-managed
    songs::Vector{String}
end

==(x::Album, y::Album) = x.id == y.id
Album() = Album(0, 0, "", "", 0, 0, String[])
Album(name, artist, year, songs) = Album(0, 0, name, artist, year, 0, songs)
StructTypes.StructType(::Type{Album}) = StructTypes.Mutable()
StructTypes.idproperty(::Type{Album}) = :id

mutable struct User
    id::Int64 # service-managed
    username::String
    password::String
end

==(x::User, y::User) = x.id == y.id
User() = User(0, "", "")
User(username::String, password::String) = User(0, username, password)
User(id::Int64, username::String) = User(id, username, "")
StructTypes.StructType(::Type{User}) = StructTypes.Mutable()
StructTypes.idproperty(::Type{User}) = :id

end # module

max_stars_repo_name: vanhp/MicroServiceJulia

max_issues_repo_name: vanhp/MicroServiceJulia

max_forks_repo_name: vanhp/MicroServiceJulia


text:
##
cd(@__DIR__())
using IterTools: groupby
using Underscores
@_ sum(foldl(∪, Set.(_))|>length, groupby(isempty, eachline("input.txt")))

## PARTE 2
@_ sum(foldl(∩, Set.(_))|>length, groupby(isempty, eachline("input.txt")))