text:
# # HK (Hegselmann and Krause) opinion dynamics model

# This example showcases
# * How to do synchronous updating of Agent properties
#   (also know as [Synchronous update schedule](http://jmckalex.org/compass/syn-and-asynch-expl.html)).
#   In a Synchronous update schedule changes made to an agent are not seen by
#   other agents until the next step, see also [Wilensky 2015, p.286](https://mitpress.mit.edu/books/introduction-agent-based-modeling)).
# * How to terminate the system evolution on demand according to a boolean function.
# * How to terminate the system evolution according to what happened on the _previous_ step.

# ## Model overview

# This is an implementation of a simple version of the
# [Hegselmann and Krause (2002)](http://jasss.soc.surrey.ac.uk/5/3/2.html) model.
# It is a model of opinion formation with the question: which
# parameters' values lead to consensus, polarization or fragmentation?
# It models interacting **groups** of agents (as opposed to interacting pairs, typical in
# the literature) in which it is assumed that if an agent disagrees too much with
# the opinion of a source of influence, the source can no longer influence the
# agent's opinion. There is then a "bound of confidence". The model shows that the
# systemic configuration is heavily dependent on this parameter's value.

# The model has the following components:

# - A set of n Agents with opinions xᵢ in the range [0,1] as attribute
# - A parameter ϵ called "bound" in (0, 0.3]
# - The update rule: at each step every agent adopts the mean of the opinions which are within
#   the confidence bound ( |xᵢ - xⱼ| ≤ ϵ).

# It is also available from the `Models` module as [`Models.hk`](@ref).

# ## Core structures
# We start by defining the Agent type and initializing the model.
# The Agent type has two fields so that we can implement the synchronous update.
using Agents
using Statistics: mean

mutable struct HKAgent <: AbstractAgent
    id::Int
    old_opinion::Float64
    new_opinion::Float64
    previous_opinon::Float64
end

# There is a reason the agent has three fields that are "the same".
# The `old_opinion` is used for the synchronous agent update, since we require access
# to a property's value at the start of the step and the end of the step.
# The `previous_opinion` is the opinion of the agent in the _previous_ step,
# as the model termination requires access to a property's value
# at the end of the previous step, and the end of the current step.

# We could, alternatively, make the three opinions a single field with vector value.

function hk_model(; numagents = 100, ϵ = 0.2)
    model = ABM(HKAgent, scheduler = fastest, properties = Dict(:ϵ => ϵ))
    for i in 1:numagents
        o = rand()
        add_agent!(model, o, o, -1)
    end
    return model
end

model = hk_model()

# Add some helper functions for the update rule. As there is a filter in
# the rule we implement it outside the `agent_step!` method. Notice that the filter
# is applied to the `:old_opinion` field.
function boundfilter(agent, model)
    filter(
        j -> abs(agent.old_opinion - j) < model.ϵ,
        [a.old_opinion for a in allagents(model)],
    )
end
nothing # hide

# Now we implement the `agent_step!`
function agent_step!(agent, model)
    agent.previous_opinon = agent.old_opinion
    agent.new_opinion = mean(boundfilter(agent, model))
end
nothing # hide

# and `model_step!`
function model_step!(model)
    for a in allagents(model)
        a.old_opinion = a.new_opinion
    end
end
nothing # hide

# From this implementation we see that to implement synchronous scheduling
# we define an Agent type with `old` and `new` fields for attributes that
# are changed via the synchronous update. In `agent_step!` we use the `old` field
# then, after updating all the agents `new` fields, we use the `model_step!`
# to update the model for the next iteration.

# ## Running the model
# The parameter of interest is now `:new_opinion`, so we assign
# it to variable `adata` and pass it to the `run!` method
# to be collected in a DataFrame.

# In addition, we want to run the model only until all agents have converged to an opinion.
# From the documentation of [`step!`](@ref) one can see that instead of specifying the
# amount of steps we can specify a function instead.
function terminate(model, s)
    if any(
        !isapprox(a.previous_opinon, a.new_opinion; rtol = 1e-12) for a in allagents(model)
    )
        return false
    else
        return true
    end
end

step!(model, agent_step!, model_step!, terminate)
model[1]

# Alright, let's wrap everything in a function and do some data collection using [`run!`](@ref).

function model_run(; kwargs...)
    model = hk_model(; kwargs...)
    agent_data, _ = run!(model, agent_step!, model_step!, terminate; adata = [:new_opinion])
    return agent_data
end

data = model_run(numagents = 100)
data[(end - 19):end, :]

# Notice that here we didn't speciy when to collect data, so this is done at every step.
# Instead, we could collect data only at the final step, by re-using the same
# function for the `when` argument:

model = hk_model()
agent_data, _ = run!(
    model,
    agent_step!,
    model_step!,
    terminate;
    adata = [:new_opinion],
    when = terminate,
)
agent_data

# Finally we run three scenarios, collect the data and plot it.
using Plots
using Random # hide
gr() # hide

plotsim(data, ϵ) = plot(
    data.step,
    data.new_opinion,
    leg = false,
    group = data.id,
    title = "epsilon = $(ϵ)",
)

Random.seed!(42) # hide

plt001, plt015, plt03 =
    map(e -> (model_run(ϵ = e), e) |> t -> plotsim(t[1], t[2]), [0.05, 0.15, 0.3])

plot(plt001, plt015, plt03, layout = (3, 1))


text:
"""
equivalent to [0; cumsum(x)]
"""
function zero_cumsum(x::Array{T, 1}) where T
    res = Array{T, 1}(undef, length(x) + 1)

    res[1] = zero(T)
    for i in 1:length(x)
        res[i + 1] = res[i] + x[i]
    end

    return res
end

function zero_cumsum(f::Function, x::Array{T, 1}) where T
    res = Array{T, 1}(undef, length(x) + 1)

    res[1] = zero(T)
    for i in 1:length(x)
        res[i + 1] = res[i] + f(x[i])
    end

    return res
end

"""
    NormalMeanSegment(data[, σ = 1.0])

Constructs a function which calculates the cost of fitting a Normal distribution
with unknown mean, and known standard deviation `σ` for a specified segment of `data`.
By default, `sigma` is set to `1.0`.

See also: [`NormalVarSegment`](@ref), [`NormalMeanVarSegment`](@ref), [`ExponentialSegment`](@ref), [`PoissonSegment`](@ref), [`GammaRateSegment`](@ref), [`GammaShapeSegment`](@ref),
[`Nonparametric`](@ref), [`PELT`](@ref), [`BS`](@ref)

# Returns
* `cost::Function`: Function which takes indices (s, t) where
                    0 ≤ s < t < n and n is the length of the time series
                    and returns the cost of fitting distribution to the segment [s+1, ..., t]

# Example
```julia-repl
n = 1000
λ = 100
μ, σ = Normal(0.0, 10.0), 1.0
data, cps = @changepoint_sampler n λ Normal(μ, σ)
cost = NormalMeanSegment(data, σ = 1.0)
cps, cost = BS(cost, n)
```
"""
function NormalMeanSegment(data::Array{Float64}, σ::Real = 1.0)
    cd = zero_cumsum(data)
    cd_2 = zero_cumsum(x -> x^2, data)
    cost(s::Int64, t::Int64) =
        ( cd_2[t+1] - cd_2[s+1] - abs2(cd[t+1] - cd[s+1]) / (t-s) ) / (σ^2)
    return cost
end

"""
    NormalVarSegment(data, μ)

Constructs a segment cost function for `data` assuming it follows a Normal distribution with changing variance
and fixed mean `μ`.

See also: [`NormalMeanSegment`](@ref)
"""
function NormalVarSegment(data::Array{Float64}, μ::Real)
    ss = zero_cumsum(x -> (x - μ)^2, data)
    cost(s::Int64, t::Int64) = (t-s) * log( (ss[t+1] - ss[s+1])/(t-s) )
end

"""
    NormalMeanVarSegment(data)

Constructs a segment cost function for `data` assuming it follows a Normal distribution with changing mean and variance.

See also: [`NormalMeanSegment`](@ref)
"""
function NormalMeanVarSegment(data::Array{Float64})
    cd = zero_cumsum(data)
    cd_2 = zero_cumsum(x -> x^2, data)
    function cost(s::Int64, t::Int64)
      mu = (cd[t + 1] - cd[s + 1]) / (t - s)
      sig = ( cd_2[t+1] - cd_2[s+1] ) / (t-s) - mu^2
      return (t-s)*log(sig) + (cd_2[t+1] - cd_2[s+1])/sig - 2*(cd[t+1] - cd[s+1])*mu/sig + ((t-s)*mu^2)/sig
    end
    return cost
end

"""
    ExponentialSegment(data)

Constructs a segment cost function for `data` assuming it follows an Exponential distribution with changing mean.

See also: [`NormalMeanSegment`](@ref)
"""
function ExponentialSegment(data::Array{Float64})
    cd = zero_cumsum(data)
    cost(s::Int64, t::Int64) = -*(t-s) * ( log(t-s) - log(cd[t+1] - cd[s+1]))
    return cost
end

"""
    PoissonSegment(data)

Constructs a segment cost function for `data` assuming it follows Poisson distribution with changing mean.

See also: [`NormalMeanSegment`](@ref)
"""
function PoissonSegment(data::Array{Float64})
    cd = zero_cumsum(data)
    cost(s::Int64, t::Int64) = -2*(cd[t+1]-cd[s+1]) * ( log(cd[t+1]-cd[s+1]) - log(t-s) - 1 )
    return cost
end

"""
    GammaShapeSegment(data, beta)

Constructs a segment cost function for `data` assuming it follows a Gamma distribution with changing shape parameter
and fixed rate parameter `beta`.

See also: [`NormalMeanSegment`](@ref)
"""
function GammaShapeSegment(data::Array{Float64}, beta::Float64)
    lcd = zero_cumsum(log, data)
    function cost(s::Int64, t::Int64)
        alpha_hat = invdigamma( log(beta) + (lcd[t+1] - lcd[s+1])/(t-s) )
        cost = (t-s)*( alpha_hat*log(beta) - lgamma(alpha_hat) ) + (alpha_hat-1)*(lcd[t+1]-lcd[s+1])
        return cost*-2
    end
end

"""
    GammaRateSegment(data, alpha)

Constructs a segment cost function for `data` assuming it follows a Gamma distribution with changing rate parameter
and fixed shape parameter `alpha`.

See also: [`NormalMeanSegment`](@ref)
"""
function GammaRateSegment(data::Array{Float64}, alpha::Float64)
    cd = zero_cumsum(data)
    function cost(s::Int64, t::Int64)
        beta_hat = ( alpha * (t-s) )/( cd[t+1] - cd[s+1] )
        cost = (t-s)*alpha*log(beta_hat) - beta_hat*( cd[t+1] - cd[s+1] )
        return cost*-2
    end
end

"""
    NonparametricSegment(data, k)

Constructs a segment cost function for `data` using the Nonparametric model with
parameter `k`.

See also: [`NormalMeanSegment`](@ref)
"""
function NonparametricSegment(data::Array{Float64}, k::Int64)
    n = length(data)
    if (k > n)
        k = n
    end
    x = sort(data)
    yk = -1 .+ (2 .* collect(1:k) ./ k .- 1/k)
    c = -log(2*n-1)
    pk = (1 .+ exp.(c .* yk)) .^ -1.0
    j = floor(Int64, (n-1)*pk[1]+1)
    Q = vcat(0,cumsum(data .< x[j]) .+ 0.5 .* cumsum(data.==x[j]))
    for i in 2:k
        j = convert(Int64,floor((n-1)*pk[i]+1))
        Q =vcat(Q,(vcat(0,cumsum(data .< x[j]) + 0.5*cumsum(data.==x[j]))))
    end

    function cost(s::Int64, t::Int64)
        Fkl = zeros(Float64, k)
        for i in 1:k
            Fkl[i] = ((Q[(t+1)+(i-1)*(n+1)]) -  (Q[(s+1)+(i-1)*(n+1)]))/(t-s)
        end
        cost = 0
        for i in 1:k
            if (Fkl[i] > 0 && Fkl[i] < 1)
                cost  = cost + (t-s)*(Fkl[i]*log(Fkl[i]) + (1-Fkl[i])*log(1-Fkl[i]))
            end
        end
        return 2*c*cost/k
    end
end

"""
    OLSSegment(data)

Create a segment cost function for piecewise linear regressions, fitted using OLS (assuming Normally distributed errors).

See also: [`NormalMeanSegment`](@ref)
"""
function OLSSegment(data::Array{Float64})
    function cost(s::Int64, t::Int64)
	      if t-s > 1
            y = data[s:t]
	          m = length(y)
	          x = 1:m
	          a = ( 2*(2*m+1)*sum(y) - 6*sum(x.*y) )/( m*(m-1) )
	          b = ( 12*sum(x.*y) - 6*(m+1)*sum(y) )/( m*(m-1)*(m+1) )
	          sig = sum( (y-a-b*x).^2 )/m
	          return m/2 * ( log(sig) + 1 )
	      else
	          return Inf
	      end
    end
end

ext: jl


text:
using Base.Test

using IntegerMask

@testset "mask" begin

# copy
inds = [2, 3]
m1 = Mask{4}(inds)
m2 = copy(m1)
push!(inds, 4)

@test m1.indices == [2, 3, 4]
@test m2.indices == [2, 3]

# length
@test length(Mask{10}()) == 0
@test length(Mask{3}(1, 2)) == 2

# equality
@test Mask{3}(1, 2) == Mask{3}(1, 2)
@test Mask{3}(1, 2) != Mask{3}(3)

# getindex array[mask]
@test_throws AssertionError (1:10)[Mask{9}(1)]
@test (1:10)[Mask{10}(1, 3, 7)] == [1, 3, 7]

# getindex mask[mask]
@test Mask{4}(1, 3, 4)[Mask{4}()] == Mask{0}()
@test Mask{4}(1, 3, 4)[Mask{4}(2, 3)] == Mask{2}(2)
@test Mask{4}(1, 3, 4)[Mask{4}(1, 2, 3)] == Mask{3}(1, 3)

# complement
@test Mask{3}(1, 2)' == Mask{3}(3)
@test Mask{5}(3)' == Mask{5}(1, 2, 4, 5)
@test Mask{4}(1, 2, 3, 4)' == Mask{4}()

# union
@test Mask{2}() ∪ Mask{2}() == Mask{2}()
@test Mask{2}() ∪ Mask{2}(1) == Mask{2}(1)
@test Mask{2}(1) ∪ Mask{2}() == Mask{2}(1)
@test Mask{4}(1, 3) ∪ Mask{4}(2, 3) == Mask{4}(1, 2, 3)
@test Mask{4}(2, 4) ∪ Mask{4}(3, 4) == Mask{4}(2, 3, 4)

# intersect
@test Mask{2}() ∩ Mask{2}() == Mask{2}()
@test Mask{2}() ∩ Mask{2}(1) == Mask{2}()
@test Mask{2}(1) ∩ Mask{2}() == Mask{2}()
@test Mask{4}(1, 3) ∩ Mask{4}(2, 3) == Mask{4}(3)
@test Mask{4}(2, 3) ∩ Mask{4}(1, 3) == Mask{4}(3)
@test Mask{4}(2, 4) ∩ Mask{4}(3, 4) == Mask{4}(4)
@test Mask{6}(1, 3, 5, 6) ∩ Mask{6}(2, 3, 4) == Mask{6}(3)
@test Mask{6}(2, 3, 4) ∩ Mask{6}(1, 3, 5, 6) == Mask{6}(3)

# setdiff
@test Mask{2}() \ Mask{2}() == Mask{2}()
@test Mask{2}() \ Mask{2}(1) == Mask{2}()
@test Mask{2}(1) \ Mask{2}() == Mask{2}(1)
@test Mask{4}(1, 3) \ Mask{4}(2, 3) == Mask{4}(1)
@test Mask{4}(2, 4) \ Mask{4}(3, 4) == Mask{4}(2)
@test Mask{4}(2, 3, 4) \ Mask{4}(3, 4) == Mask{4}(2)

# in
@test 1 ∈ Mask{2}(1)
@test 2 ∉ Mask{2}(1)

# sum
@test_throws AssertionError sum([1, 2], Mask{1}())

c = sum(1:4, Mask{4}())
@test isa(c, Int)
@test c == 0

c = sum([1., 2., 3., 4.], Mask{4}())
@test isa(c, Float64)
@test c == 0

c = sum(complex([1., 2., 3., 4.]), Mask{4}())
@test isa(c, Complex128)
@test c == 0

@test sum(1:4, Mask{4}(2, 4)) == 6

# prod
@test_throws AssertionError prod([1, 2], Mask{1}())

c = prod(1:4, Mask{4}())
@test isa(c, Int)
@test c == 1

c = prod([1., 2., 3., 4.], Mask{4}())
@test isa(c, Float64)
@test c == 1

c = prod(complex([1., 2., 3., 4.]), Mask{4}())
@test isa(c, Complex128)
@test c == 1

@test prod(1:4, Mask{4}(2, 4)) == 8

# reduce
@test_throws AssertionError reduce(+, [1, 2], Mask{1}())
@test_throws AssertionError reduce(+, [], Mask{0}())

c = reduce(*, 1:4, Mask{4}(2, 4))
@test isa(c, Int)
@test c == 8

c = reduce(+, [[1,3], [5,2], [6,3], [2, -9]], Mask{4}(1, 3))
@test isa(c, Vector{Int})
@test c == [7, 6]

end # testset

max_issues_repo_name: bastikr/IntegerMask

max_forks_repo_name: bastikr/IntegerMask


text:
"""
    TypeB(var_1::TypeC)

Lorem ipsum dolor sit amet.
"""
struct TypeB

  var_1::TypeC

end


text:
module SpatialMaths

export AbstractPoint, AbstractFrame, AbstractPose
export Point, Frame, Pose
export Point2, Frame2, Pose2
export rot2, trot2, transl2, se2, so2
export Point3, Frame3, Pose3
export rotx, roty, rotz, trotx, troty, trotz, transl3, se3, so3
export transl

include("spatialmaths/pose.jl")
include("spatialmaths/pose_graph.jl")
include("spatialmaths/point.jl")
include("spatialmaths/frame.jl")
include("spatialmaths/transforms.jl")
include("spatialmaths/spatial_types.jl")

end # (sub)module SpatialMaths


text:
@test SpecificV(0.5, 1500.0) ≈ 0.138_455_090E1
@test SpecificV(0.5u"MPa", 1500.0u"K") ≈ 0.138_455_090E1u"m^3/kg"

@test SpecificV(30.0, 1500.0) ≈ 0.230_761_299E-1
@test SpecificV(30.0u"MPa", 1500.0u"K") ≈ 0.230_761_299E-1u"m^3/kg"

@test SpecificV(30.0, 2000.0) ≈ 0.311_385_219E-1
@test SpecificV(30.0u"MPa", 2000.0u"K") ≈ 0.311_385_219E-1u"m^3/kg"

@test SpecificH(0.5, 1500.0) ≈ 0.521_976_855E4
@test SpecificH(0.5u"MPa", 1500.0u"K") ≈ 0.521_976_855E4u"kJ/kg"

@test SpecificH(30.0, 1500.0) ≈ 0.516_723_514E4
@test SpecificH(30.0u"MPa", 1500.0u"K") ≈ 0.516_723_514E4u"kJ/kg"

@test SpecificH(30.0, 2000.0) ≈ 0.657_122_604E4
@test SpecificH(30.0u"MPa", 2000.0u"K") ≈ 0.657_122_604E4u"kJ/kg"

@test SpecificU(0.5, 1500.0) ≈ 0.452_749_310E4
@test SpecificU(0.5u"MPa", 1500.0u"K") ≈ 0.452_749_310E4u"kJ/kg"

@test SpecificU(30.0, 1500.0) ≈ 0.447_495_124E4
@test SpecificU(30.0u"MPa", 1500.0u"K") ≈ 0.447_495_124E4u"kJ/kg"

@test SpecificU(30.0, 2000.0) ≈ 0.563_707_038E4
@test SpecificU(30.0u"MPa", 2000.0u"K") ≈ 0.563_707_038E4u"kJ/kg"

@test SpecificS(0.5, 1500.0) ≈ 0.965_408_875E1
@test SpecificS(0.5u"MPa", 1500.0u"K") ≈ 0.965_408_875E1u"kJ/kg/K"

@test SpecificS(30.0, 1500.0) ≈ 0.772_970_133E1
@test SpecificS(30.0u"MPa", 1500.0u"K") ≈ 0.772_970_133E1u"kJ/kg/K"

@test SpecificS(30.0, 2000.0) ≈ 0.853_640_523E1
@test SpecificS(30.0u"MPa", 2000.0u"K") ≈ 0.853_640_523E1u"kJ/kg/K"

@test SpecificCP(0.5, 1500.0) ≈ 0.261_609_445E1
@test SpecificCP(0.5u"MPa", 1500.0u"K") ≈ 0.261_609_445E1u"kJ/kg/K"

@test SpecificCP(30.0, 1500.0) ≈ 0.272_724_317E1
@test SpecificCP(30.0u"MPa", 1500.0u"K") ≈ 0.272_724_317E1u"kJ/kg/K"

@test SpecificCP(30.0, 2000.0) ≈ 0.288_569_882E1
@test SpecificCP(30.0u"MPa", 2000.0u"K") ≈ 0.288_569_882E1u"kJ/kg/K"

@test SpeedOfSound(0.5, 1500.0) ≈ 0.917_068_690E3
@test SpeedOfSound(0.5u"MPa", 1500.0u"K") ≈ 0.917_068_690E3u"m/s"

@test SpeedOfSound(30.0, 1500.0) ≈ 0.928_548_002E3
@test SpeedOfSound(30.0u"MPa", 1500.0u"K") ≈ 0.928_548_002E3u"m/s"

@test SpeedOfSound(30.0, 2000.0) ≈ 0.106_736_948E4
@test SpeedOfSound(30.0u"MPa", 2000.0u"K") ≈ 0.106_736_948E4u"m/s"


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

using Markdown
using InteractiveUtils

# ╔═╡ 4b5b292a-699e-11eb-0b9f-7daa58c4eeba
begin
	using DataFrames
	using TimeSeries
	using Dates
	using StatsPlots
end

# ╔═╡ 8ab35a9c-69a6-11eb-2112-1733984fe260
md"""
TODO: pp 50 from the book
Try working with the RData format instead.
"""

# ╔═╡ b92c9ec8-699f-11eb-2416-dd71aa4bf2d9
gr();

# ╔═╡ 574c6094-69a2-11eb-35d9-df19dc196d0d
full_df = DataFrame();

# ╔═╡ ae93f34e-69a2-11eb-216d-53292bcbd85a
# files = Dict("cmort" => "cmort.dat",
# 	"temp" => "temp.dat",
# 	"part" => "particle.dat"
# );

# ╔═╡ 07e736d8-69a6-11eb-09c5-ab3114753f2b
files = Dict(:cmort => "cmort.dat");

# ╔═╡ 653a7668-699e-11eb-022e-d79092879b0c
function read_files_to_df!(df, k, v)
	base_path = "../data"
	open(joinpath(base_path, v), "r") do io
		df.k = read(io)
	end
end;

# ╔═╡ e313c112-69a2-11eb-1a37-b9695c591827
read_files_to_df!(full_df, "cmort", get(files, :cmort, 1))

# ╔═╡ 552adfe4-69a6-11eb-1d17-4dcde0b9ffe2
first(full_df, 3)

# ╔═╡ 71961502-699e-11eb-3626-d504ba310af6
dates = Date(1970, 1, 1):Day(7):Date(1979, 10);

# ╔═╡ af740eb0-699e-11eb-2ef7-f5785d24db12
cmort_ar = TimeArray(dates[1:end-1], cmort[!, :cmort], [:cmort]);

# ╔═╡ e3071582-699f-11eb-1657-792d26702790
df_cmort = DataFrame(cmort_ar);

# ╔═╡ 26f05f74-69a0-11eb-32f7-97e305b6aba1
@df df_cmort plot(:timestamp, :cmort, leg=false)

# ╔═╡ 3727b36c-69a0-11eb-15e2-237d37b7be34

# ╔═╡ 679228c0-699e-11eb-1095-0f981944a5f1

# ╔═╡ 677c9046-699e-11eb-20b9-e37da734c1a3

# ╔═╡ 6763604e-699e-11eb-1f20-b716b1e3a4d1

# ╔═╡ 6749ae56-699e-11eb-06c1-5f35812d1b06

# ╔═╡ 672ef066-699e-11eb-3ab5-b972365504fe

# ╔═╡ 6715ed6e-699e-11eb-365e-cb335f1e8a73

# ╔═╡ 66fb044a-699e-11eb-300d-fd661918f5bd

# ╔═╡ 66e4225c-699e-11eb-132f-47a9bd831a87

# ╔═╡ 66ca2744-699e-11eb-2412-fb201a7c3e98

# ╔═╡ 66ae1838-699e-11eb-14db-c9fe4e1b5b45

# ╔═╡ 66969334-699e-11eb-05f9-498acac4c604

# ╔═╡ 66772d0a-699e-11eb-1a76-ed5e301953dd

# ╔═╡ 6657f7fa-699e-11eb-1992-8521c50bd37b

# ╔═╡ Cell order:
# ╠═8ab35a9c-69a6-11eb-2112-1733984fe260
# ╠═4b5b292a-699e-11eb-0b9f-7daa58c4eeba
# ╠═b92c9ec8-699f-11eb-2416-dd71aa4bf2d9
# ╠═574c6094-69a2-11eb-35d9-df19dc196d0d
# ╠═ae93f34e-69a2-11eb-216d-53292bcbd85a
# ╠═07e736d8-69a6-11eb-09c5-ab3114753f2b
# ╠═653a7668-699e-11eb-022e-d79092879b0c
# ╠═e313c112-69a2-11eb-1a37-b9695c591827
# ╠═552adfe4-69a6-11eb-1d17-4dcde0b9ffe2
# ╠═71961502-699e-11eb-3626-d504ba310af6
# ╠═af740eb0-699e-11eb-2ef7-f5785d24db12
# ╠═e3071582-699f-11eb-1657-792d26702790
# ╠═26f05f74-69a0-11eb-32f7-97e305b6aba1
# ╠═3727b36c-69a0-11eb-15e2-237d37b7be34
# ╠═679228c0-699e-11eb-1095-0f981944a5f1
# ╠═677c9046-699e-11eb-20b9-e37da734c1a3
# ╠═6763604e-699e-11eb-1f20-b716b1e3a4d1
# ╠═6749ae56-699e-11eb-06c1-5f35812d1b06
# ╠═672ef066-699e-11eb-3ab5-b972365504fe
# ╠═6715ed6e-699e-11eb-365e-cb335f1e8a73
# ╠═66fb044a-699e-11eb-300d-fd661918f5bd
# ╠═66e4225c-699e-11eb-132f-47a9bd831a87
# ╠═66ca2744-699e-11eb-2412-fb201a7c3e98
# ╠═66ae1838-699e-11eb-14db-c9fe4e1b5b45
# ╠═66969334-699e-11eb-05f9-498acac4c604
# ╠═66772d0a-699e-11eb-1a76-ed5e301953dd
# ╠═6657f7fa-699e-11eb-1992-8521c50bd37b

max_issues_repo_name: BioLogs/TimeSeriesConcepts


text:
using YaoZX
using Test

@testset "push gate" begin
    include("push_gate.jl")
end

@testset "to qbir" begin
    include("toqbir.jl")
end

@testset "simplification" begin
    include("simplification.jl")
end


text:
function randomized_first_improvement(tuning_run::Run,
                                      channel::RemoteChannel;
                                      cutoff::Integer          = 10_000,
                                      threshold::AbstractFloat = 0.6)
    function iterate(tuning_run::Run)
        if rand() <= threshold
            return random_walk(tuning_run)
        else
            return first_improvement(tuning_run, cutoff = cutoff)
        end
    end

    technique(tuning_run,
              channel,
              iterate,
              name = "Randomized First Improvement")
end


text:
# In this example, we use EKI to tune the closure parameters of a HydrostaticFreeSurfaceModel
# with a CATKEBasedVerticalDiffusivity closure in order to align the predictions of the model
# to those of a high-resolution LES data generated in LESbrary.jl. Here `predictions` refers to the
# 1-D profiles of temperature, velocity, and turbulent kinetic energy horizontally averaged over a
# 3-D physical domain.

pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "projects", "OceanBoundaryLayerParameterizations", "src"))

image_dir = joinpath(@__DIR__, "quick_calibrate")
mkpath(image_dir)

using Oceananigans
using OceanTurbulenceParameterEstimation
using OceanBoundaryLayerParameterizations

# CATKE parameters involved in setting field diffusivities
@free_parameters StabilityFnParameters CᴷRiʷ CᴷRiᶜ Cᴷu⁻ Cᴷuʳ Cᴷc⁻ Cᴷcʳ Cᴷe⁻ Cᴷeʳ

parameters = Parameters(
    RelevantParameters = CATKEParametersRiDependent,  # Parameters that are used in CATKE
    ParametersToOptimize = StabilityFnParameters      # Subset of RelevantParameters that we want to optimize
)

two_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/2DaySuite"
four_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/4DaySuite"
six_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/6DaySuite"

two_day_suite = TwoDaySuite(two_day_suite_dir)
four_day_suite = FourDaySuite(four_day_suite_dir)
six_day_suite = SixDaySuite(six_day_suite_dir)

# InverseProblem represents the model, data, loss function, and parameters
# calibration = InverseProblem(two_day_suite, # "Truth data" for model calibration
#                              parameters;   # Model parameters
#                              # Loss function parameters
#                              relative_weights = relative_weight_options["all_but_e"],
#                              # Model (hyper)parameters
#                              architecture = GPU(),
#                              ensemble_size = 10,
#                              Nz = 16,
#                              Δt = 30.0)

calibration = InverseProblem(two_day_suite, parameters; relative_weights = relative_weight_options["all_but_e"],
                             architecture = GPU(), ensemble_size = 10, Δt = 30.0)

validation = InverseProblem(four_day_suite, calibration; Nz = 64);

# Loss on default parameters
l0 = calibration()

# Example parameters
θ = calibration.default_parameters

# Loss on parameters θ.
# θ can be
#   1. a vector
#   2. a FreeParameters object
#   3. a vector of parameter vectors (one for each ensemble member)
#   4. or a vector of FreeParameter objects (one for each ensemble member)
# If (1) or (2), the ensemble members are redundant and the loss is computed for just the one parameter set.
lθ = calibration(θ)

# Output files/figures
directory = joinpath(pwd(), "quick_calibrate")

# Run the model forward and store the solution
output = model_time_series(calibration, θ)

# Run the model forward with parameters θ and visualize the solution compared to the truth
visualize!(output; filename = "visualize_default_parameters.png")

#=

# Use EKI to calibrate the model parameters
calibration_algorithm = EKI(noise_level = 10^(-2.0),
                            N_iter = 15,
                            stds_within_bounds = 0.6,
                            informed_priors = false)

best_parameters = calibrate(calibration; algorithm = calibration_algorithm)

loss = calibration(best_parameters)

=#

# Runs `visualize!` and records a summary of the calibration results in a `result.txt` file.
visualize_and_save!(calibration, validation, θ, directory)


text:
"""
	abstract type AbstractCell{T,S,U}

The basic type of J4PR. Most of the types are its subtypes to ensure that operations
between various types of `Cells` are possible, especially when creating `Pipes` that
can have arbitrary levels of nesting.
"""
abstract type AbstractCell{T,S,U} end

#############################
# Operators [AbstractCells] #
#############################

# Unary operators for field access
+(ac::T where T<:AbstractCell) = ac.x
-(ac::T where T<:AbstractCell) = ac.y
~(ac::T where T<:AbstractCell) = dump(ac)

# Define any new aliases used throughout subsequent code
const PTuple{T} = Tuple{Vararg{<:T}}


text:
export sgiter

"""
    λ,v = sgiter([eltype],nep::NEP,j::Integer;[λ_min,][λ_max,][λ,][errmeasure,][tol,][maxit,][logger,][eigsolvertype::Type,])

Finds the `j`-th eigenvalue of the NEP using safeguarded iteration, with eigenvalue numbering
according to min-max theory. The method only works for Hermitian problems, and the eigenvalues
are assumed to be real. If an interval [`λ_min`,`λ_max`] is given, then the Rayleigh functional
is assumed to be unique on the interval. If no interval is given, then the minimum solution is
always taken. The method requires the computation of (all) eigenvalues of a matrix. The `eigsolvertype`
is a `Type` that specifies which eigevalue solver is used inside the algorithm.

See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> nep = nep_gallery("real_quadratic");
julia> λ,v = sgiter(nep, 1, λ_min = -10, λ_max = 0,  λ = -10, maxit = 100);
julia> minimum(svdvals(compute_Mder(nep,λ)))
0.0
julia> norm(v)
1.0
```

# References
* V. Mehrmann and H. Voss, Nonlinear eigenvalue problems: a challenge for modern eigenvalue methods, GAMM‐Mitteilungen 27.2 (2004): 121-152.
* H. Voss and B. Werner, Solving sparse nonlinear eigenvalue problems. Technical Report 82/4, Inst. f. Angew. Mathematik, Universität Hamburg, 1982.
* B. Werner. Das Spektrum von Operatorenscharen mit verallgemeinerten Rayleighquotienten. PhD thesis, Fachbereich Mathematik, Universität Hamburg, 1970

"""
sgiter(nep::NEP, j::Integer; params...) = sgiter(ComplexF64, nep, j; params...)
function sgiter(::Type{T},
                   nep::NEP,
                   j::Integer;
                   λ_min::Real = NaN,
                   λ_max::Real = NaN,
                   λ::Number = zero(real(T)),
                   errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
                   tol::Real = eps(real(T)) * 100,
                   maxit::Integer = 100,
                   inner_solver=@default_compute_rf_inner_solver(nep),
                   logger = 0,
                   eigsolvertype::Type = DefaultEigSolver
                   ) where {T<:Number}

    @parse_logger_param!(logger)

    n = size(nep,1)
    if (j > n) || (j <= 0)
       error("j must be a number between 1 and size(nep) = ", n, ". You have selected j = ", j, ".")
    end
    real_T = real(T)
    λ_min::real_T = real_T(λ_min)
    λ_max::real_T = real_T(λ_max)
    if !( (isnan(λ_min) && isnan(λ_max)) || (!isnan(λ_min) && !isnan(λ_max)) )
        error("A proper interval is not chosen.") #Both should be either given or not given
    end
    if (λ_max < λ_min)
       error("The interval cannot be empty, λ_max >= λ_min is required. Supplied interval [λ_min,λ_max] = [", λ_min, ",", λ_max, "].")
    end
    λ::real_T = real_T(λ)
    if ( (λ < λ_min) || (λ > λ_max) )
       error("The starting guess is outside the interval.")
    end

    v::Array{T,1} = zeros(T,n)
    err = 0

    for k = 1:maxit
       eig_solver = eigsolvertype(compute_Mder(nep, λ, 0))
       v[:] = compute_jth_eigenvector(eig_solver, nep, λ, j)
       λ_vec = compute_rf(real_T, nep, v, inner_solver, TOL = tol/10)
       push_info!(logger,2,"compute_rf: $λ_vec")
       λ = choose_correct_eigenvalue_from_rf(λ_vec, λ_min, λ_max)
       err = estimate_error(errmeasure,λ, v)
       push_iteration_info!(logger,k,err=err,λ=λ);
       if (err < tol)
           return (λ, v)
       end
    end

    msg = "Number of iterations exceeded. maxit = $(maxit)."
    throw(NoConvergenceException(λ, v, err ,msg))
end

function compute_jth_eigenvector(eig_solver, nep::NEP, λ, j)
    n = size(nep,1)
    Λ, V = eig_solve(eig_solver, nev = n)
    if (n>1)
        p = sortperm(Λ)
        return vec(V[:,p[j]])
    else
        return V
    end
end

function choose_correct_eigenvalue_from_rf(λ_vec, λ_min, λ_max)
    if ( isnan(λ_min) && isnan(λ_max) )
        return minimum(λ_vec)
    else
        idxes = (λ_vec .<= λ_max) .& (λ_vec .>= λ_min)
        if sum(idxes) > 1
            error("Multiple values of λ found in the interval.")
        end
        if sum(idxes) == 0
            error("No λ found in the prescribed interval.")
        end
        idx = findall(idxes)[1]::Int #This vector is only 1 element
        return λ_vec[idx]
    end
end


text:
module data_repeatrate
function Onetree(df,x1,x2,y1,y2)
    dfm=df[(df.X.>=x1).&(df.X.<=x2).&(df.Y.>=y1).&(df.Y.<=y2).&(df.Close.>0), :]
    return dfm, dfm.Close
end
function Areafruit(df_onetree)
    avgx=sum(df_onetree.X.*df_onetree.Close)/size(df_onetree)[1];
    avgy=sum(df_onetree.Y.*df_onetree.Close)/size(df_onetree)[1];
    index=findmax(df_onetree.Close)[2];
    modex=df_onetree.X[index];#近似的树的中心
    modey=df_onetree.Y[index];
    #xmax=findmax(df.X)[1];
    #ymax=findmax(df.Y)[1];
    distance_x=df_onetree.X.-modex;
    distance_y=df_onetree.Y.-modey;
    R=(distance_x.^2+distance_y.^2).^0.5;
    df_onetree.distance = R;
    c1=df_onetree[df_onetree.distance.<=1,:];
    c11=findmax(c1.Close)[1];
    c2=df_onetree[(df_onetree.distance.<=3).&(df_onetree.distance.>1),:];
    c22=sum(c2.Close)/size(c2)[1]*8;
    c3=df_onetree[(df_onetree.distance.<=5).&(df_onetree.distance.>3),:];
    c33=sum(c3.Close)/size(c3)[1]*16;

    d=sum((df_onetree.X.-avgx).^2+(df_onetree.Y.-avgy).^2)/size(df_onetree)[1];
    #repeatratio=(acos(d/2)/pi-(d*sqrt(1-d^2/4))/2)/pi
    r=sqrt(findmax((df_onetree.X.-avgx).^2+(df_onetree.Y.-avgy).^2)[1]);
    number_fruit=c11+c22+c33;
    #return modex,modey,number_fruit
    return number_fruit
end
end


text:
@testset "Path Normalization" begin
    mixed_path = ".\\\\path///to\\a\\place/..\\///"
    norm_path = "path/to/a"
    @test Zarr.normalize_path(mixed_path) == norm_path
    @test Zarr.DirectoryStore(mixed_path).folder == norm_path
    @test Zarr.normalize_path("/") == "/"
    @test Zarr.normalize_path("/a/") == "/a"
    @test Zarr.normalize_path("/path/to/a") == "/path/to/a"
end

"""
Function to test the interface of AbstractStore. Every complete implementation should pass this test.
"""
function test_store_common(ds)
  @test !Zarr.is_zgroup(ds)
  ds[".zgroup"]=rand(UInt8,50)
  @test Zarr.is_zgroup(ds)
  @test !Zarr.is_zarray(ds)

  @test Zarr.zname(ds)=="foo"
  @test isempty(Zarr.subdirs(ds))
  @test sort(collect(Zarr.keys(ds)))==[".zgroup"]

  #Create a subgroup
  snew = Zarr.newsub(ds,"bar")
  @test !Zarr.is_zarray(ds)
  snew[".zarray"] = rand(UInt8,50)
  @test Zarr.is_zgroup(ds)
  #Test getindex and setindex
  data = rand(UInt8,50)
  snew["0.0.0"] = data
  @test snew["0.0.0"]==data
  @test Zarr.storagesize(snew)==50
  @test Zarr.isinitialized(snew,"0.0.0")
  @test !Zarr.isinitialized(snew,"0.0.1")
  Zarr.writeattrs(snew,Dict("a"=>"b"))
  @test Zarr.getattrs(snew)==Dict("a"=>"b")
  snew2 = Zarr.getsub(ds,"bar")
  @test Zarr.getattrs(snew2)==Dict("a"=>"b")
  @test snew2["0.0.0"]==data
  delete!(snew2,"0.0.0")
  @test !Zarr.isinitialized(snew2,"0.0.0")
  snew["0.0.0"] = data
end

@testset "DirectoryStore" begin
  A = fill(1.0, 30, 20)
  chunks = (5,10)
  metadata = Zarr.Metadata(A, chunks; fill_value=-1.5)
  p = tempname()
  mkpath(joinpath(p,"foo"))
  ds = Zarr.DirectoryStore(joinpath(p,"foo"))
  test_store_common(ds)
  @test isdir(joinpath(p,"foo"))
  @test isfile(joinpath(p,"foo",".zgroup"))
  @test isdir(joinpath(p,"foo","bar"))
  @test isfile(joinpath(p,"foo","bar","0.0.0"))
  @test isfile(joinpath(p,"foo","bar",".zarray"))
  @test Zarr.path(ds)==joinpath(p,"foo")
end

@testset "DictStore" begin
  A = fill(1.0, 30, 20)
  chunks = (5,10)
  metadata = Zarr.Metadata(A, chunks; fill_value=-1.5)
  ds = Zarr.DictStore("foo")
  test_store_common(ds)
  @test ds.name == "foo"
  @test haskey(ds.a,".zgroup")
  @test ds.subdirs["bar"].a[".zattrs"]==UInt8[0x7b, 0x22, 0x61, 0x22, 0x3a, 0x22, 0x62, 0x22, 0x7d]
  @test sort(collect(keys(ds.subdirs["bar"].a)))==[".zarray", ".zattrs", "0.0.0"]
  @test isempty(keys(ds.subdirs["bar"].subdirs))
end

# @testset "AWS S3 Storage" begin
#     bucket = "zarr-demo"
#     store = "store/foo/"
#     region = "eu-west-2"
#     S3 = S3Store(bucket, store, region)
#     @test storagesize(S3) == 69
#     @test Zarr.zname(S3) == "foo"
#     @test Zarr.is_zgroup(S3) == true
#     S3group = zopen(S3)
#     @test Zarr.zname(S3group) == "foo"
#     S3Array = S3group.groups["bar"].arrays["baz"]
#     @test Zarr.zname(S3Array) == "baz"
#     @test eltype(S3Array) == Zarr.ASCIIChar
#     @test storagesize(S3Array) == 69
#     @test String(S3Array[:]) == "Hello from the cloud!"
# end


text:
using Test, NiLangCore
using NiLangCore: SymbolTable, allocate!, deallocate!, operate!, swapvars!, variable_analysis_ex

@testset "variable analysis" begin
    st = SymbolTable()
    # allocate! : not exist
    allocate!(st, :x)
    allocate!(st, :y)
    @test st.existing == [:x, :y]
    # allocate! : existing
    @test_throws InvertibilityError allocate!(st, :x)
    @test st.existing == [:x, :y]
    # deallocate! : not exist
    @test_throws InvertibilityError deallocate!(st, :z)
    # deallocate! : existing
    deallocate!(st, :y)
    @test st.existing == [:x]
    @test st.deallocated == [:y]
    # deallocate! : deallocated
    @test_throws InvertibilityError deallocate!(st, :y)
    # operate! : deallocated
    @test_throws InvertibilityError operate!(st, :y)
    # allocate! : deallocated
    allocate!(st, :y)
    @test st.existing == [:x, :y]
    @test st.deallocated == []
    # operate! : not exist
    operate!(st, :j)
    @test st.unclassified == [:j]
    # operate! : existing
    operate!(st, :y)
    @test st.unclassified == [:j]
    # allocate! unclassified
    @test_throws InvertibilityError allocate!(st, :j)
    # operate! : unclassified
    operate!(st, :j)
    @test st.unclassified == [:j]
    # deallocate! : unclassified
    @test_throws InvertibilityError deallocate!(st, :j)

    # swap both existing
    swapvars!(st, :j, :x)
    @test st.unclassified == [:x]
    @test st.existing == [:j, :y]

    # swap existing - nonexisting
    swapvars!(st, :j, :k)
    @test st.unclassified == [:x, :j]
    @test st.existing == [:k, :y]

    # swap nonexisting - existing
    swapvars!(st, :o, :x)
    @test st.unclassified == [:o, :j, :x]
    @test st.existing == [:k, :y]

    # swap both not existing
    swapvars!(st, :m, :n)
    @test st.unclassified == [:o, :j, :x, :m, :n]

    # push and pop variables
end

@testset "variable analysis" begin
    st = SymbolTable([:x, :y], [], [])
    ex = :((x,y) ↔ (a, b))
    variable_analysis_ex(ex, st)
    @test st.existing == [:a, :b]
    @test st.unclassified == [:x, :y]
    st = SymbolTable([:x, :y], [], [])
    ex = :((x,y) ↔ b)
    variable_analysis_ex(ex, st)
    @test st.existing == [:b]
    @test st.unclassified == [:x, :y]
    ex = :(b ↔ (x,y))
    variable_analysis_ex(ex, st)
    @test st.existing == [:x, :y]
    @test st.unclassified == [:b]

    st = SymbolTable([:x, :y], [], [])
    ex = :(b ↔ x)
    variable_analysis_ex(ex, st)
    @test st.existing == [:b, :y]
    @test st.unclassified == [:x]

    st = SymbolTable([], [], [])
    ex = :(b ↔ (x, y))
    variable_analysis_ex(ex, st)
    @test st.existing == []
    @test st.unclassified == [:b, :x, :y]
end


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@bench @testset "cd player native" begin
    nep = nep_gallery("nlevp_native_cd_player")
    n = size(nep, 1)

    (Λ,V)=polyeig(nep)

    (λ,v)=quasinewton(nep,λ=Λ[1],v=V[:,1],logger=displaylevel, armijo_factor=0.5,armijo_max=10, tol=1e-10)
    verify_lambdas(1, nep, λ, v, 1e-10)

    (λ,v)=resinv(nep,λ=Λ[2],v=V[:,2],logger=displaylevel, armijo_factor=0.5,armijo_max=10, tol=1e-10)
    verify_lambdas(1, nep, λ, v, 1e-10)

    @testset "Errors thrown" begin
        @test_throws MethodError nep_gallery("nlevp_native_cd_player", 15)
        @test_throws MethodError nep_gallery("nlevp_native_cd_player", t=15)
    end

end


text:
struct AdamIterable{T, F, S, R}
    w0::T
    f::F
    stepsize::S
    beta1::R
    beta2::R
    epsilon::R
end

AdamIterable(w0, f; stepsize=1e-3, beta1=0.9, beta2=0.999, epsilon=1e-8) = AdamIterable(w0, f, to_iterator(stepsize), beta1, beta2, epsilon)

Base.IteratorSize(::Type{<:AdamIterable}) = Base.IsInfinite()

mutable struct AdamState{T, G, S, R}
    w::T
    f_w
    grad_f_w::G
    stepsize_iterator::S
    m::T
    v::T
    beta1t::R
    beta2t::R
end

function Base.iterate(iter::AdamIterable)
    w = copy(iter.w0)
    grad_f_w, f_w = gradient(iter.f, w)
    state = AdamState(
        w, f_w, grad_f_w, Iterators.Stateful(iter.stepsize),
        zero(w), zero(w), iter.beta1, iter.beta2
    )
    return IterationOutput(state.w, state.f_w, state.grad_f_w), state
end

function Base.iterate(iter::AdamIterable, state::AdamState)
    state.m .= iter.beta1 .* state.m .+ (1 - iter.beta1) .* state.grad_f_w
    state.v .= iter.beta2 .* state.v .+ (1 - iter.beta2) .* state.grad_f_w .^ 2
    stepsize = popfirst!(state.stepsize_iterator)
    alpha_t = stepsize * sqrt(1 - state.beta2t)/(1 - state.beta1t)
    state.w .-= alpha_t .* state.m ./ (sqrt.(state.v) .+ iter.epsilon)
    state.beta1t *= iter.beta1
    state.beta2t *= iter.beta2
    state.grad_f_w, state.f_w = gradient(iter.f, state.w)
    return IterationOutput(state.w, state.f_w, state.grad_f_w), state
end

struct Adam
    kwargs
    Adam(; kwargs...) = new(kwargs)
end

(alg::Adam)(args...) = AdamIterable(args...; alg.kwargs...)


text:
module ProposalFunctions

export swapTwo, widescaleRandomNoise

function swapTwo(x)
    n = length(x)
    index1 = rand(1:n)
    index2 = rand(1:n)
    new_x = x
    new_x[index1], new_x[index2] = x[index2], x[index1]
    new_x
end

function widescaleRandomNoise(x::Vector{Float64})
    n = length(x)
    expon = rand()*5.0 - 2.0
    displacement = randn(n)
    new_x = x + displacement*2^expon
end

function widescaleRandomNoise(x::Float64)
    expon = rand()*5.0
    displacement = randn()
    new_x = x + displacement * 2^expon
end

end

max_stars_repo_name: jacob-ruth/pt-optimization

max_issues_repo_name: jacob-ruth/pt-optimization

max_forks_repo_name: jacob-ruth/pt-optimization


text:
module RigidToricGraph

using Oscar
using Oscar.Polymake
# Write your package code here.

function get_threefaces(HD, bad_flag)
    threefaceics = @Polymake.convert_to Array{Int} Polymake.graph.nodes_of_rank(HD, 3)
    threefaceics = Polymake.to_one_based_indexing(threefaceics)
    F = @Polymake.convert_to Array{Set{Int}} HD.FACES
    F = Polymake.to_one_based_indexing(F)
    threefaces = []
    for l in threefaceics
        tf = F[l]
        if !bad_flag &&length(tf)==3  # https://xkcd.com/927
            push!(threefaces, tf)
        elseif bad_flag && length(tf)>3
            push!(threefaces, tf)
        end
    end
    return threefaces
end

function get_twofaces(HD)
    twofaceics = @Polymake.convert_to Array{Int} Polymake.graph.nodes_of_rank(HD, 2)
    twofaceics = Polymake.to_one_based_indexing(twofaceics)
    F = @Polymake.convert_to Array{Set{Int}} HD.FACES
    F = Polymake.to_one_based_indexing(F)
    twofaces = []
    for l in twofaceics
        tf = F[l]
        push!(twofaces, tf)
        end
    return twofaces
end

function crosscut_skeleton(Cdual, HD, Defdegree)
    skeleton = []
    good_faces = get_threefaces(HD, false)
    append!(good_faces, get_twofaces(HD))
    for gtf in good_faces
        is_good = true
        for r in gtf
            prod = dot(rays(Cdual)[r], Defdegree)
            if prod < 1
                is_good = false
                break
            end
        end
        if length(gtf) == 2
            if findfirst(s->issubset(gtf, s), skeleton) != nothing
                is_good = false
            end
        end
        if is_good
            push!(skeleton, gtf)
        end
    end
    return skeleton
end

export
    crosscut_skeleton,
    get_threefaces,
    get_twofaces

end


text:
# =========================================================================== #
# Ode Problem
# =========================================================================== #

struct OdeProblem
    # add mass matrix
    # add event_fcn
    # add output_fcn
    ode_fcn
    time_interval::TimeInterval
    init_cond::Vector{Float64}
end

function n_dim(ode_problem::OdeProblem)
    return length(ode_problem.init_cond)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using DataStructures, Plots
#include("generator.jl")
#include("line_tools.jl")

#########################################################################
# ------------------- TYPES FOR THE BO ALGORITHM  -----------------------
#########################################################################

import Base:<

struct Event
    point::Array{Float64, 1}
    nature::String
    assets::Union{RootedLine, Tuple{RootedLine, RootedLine}}
end

function <(e::Event, f::Event)
    return e.point < f.point
end

#########################################################################
# ------------------- PREPROCESSING  ------------------------------------
#########################################################################

function insert_endpoints_of_one_line!(X::AVLTree{Event}, l::RootedLine, r::Real)
    (startpoint, endpoint) = intersection_with_circle(l, r)
    startevent = Event(startpoint, "startpoint", l)
    endevent = Event(endpoint, "endpoint", l)
    insert!(X, startevent) ; insert!(X, endevent)
end

function initialize_event_queue(H::Array{Line, 1}, r::Real)
    X = AVLTree{Event}()
    t = SweepLine(-r-1)
    R = [RootedLine(l, t) for l in H]
    for l in R
        insert_endpoints_of_one_line!(X, l, r)
    end
    return R, X, t
end

#########################################################################
# ------------------- HANDLING EVENTS  ----------------------------------
#########################################################################

function check_intersection_and_update!(a::Real, X::AVLTree{Event}, list, l1::RootedLine, l2::RootedLine, r::Real)
    """
    Checks if l1 and l2 intersect to the right of the sweep line at t and inside the disk of
    radius r. If yes, checks if the intersection point had already been seen. If no,
    registers this point in the list and add the event to the X-structure.
    """
    if !haskey(list, Set([l1, l2])) #do nothing
        point::Array{Float64} = l1 ∩ l2
        if normof(point) < r && point[1] > a
            parents = point_in_line(l1, a)<point_in_line(l2, a) ? (l1, l2) : (l2, l1)
            e = Event(point, "intersection", parents)
            insert!(X, e) ; list[Set([l1, l2])] = point
        end
    end
end

function handle_startpoint!(t::SweepLine, X::AVLTree{Event}, ℒ::AVLTree{RootedLine}, list, e::Event, r::Real)
    point = e.point
    ℓ = e.assets
    t.r = point[1] #update the sweep line
    insert!(ℒ, ℓ)
    line_index = sorted_rank(ℒ, ℓ)

    if ℒ.count>1
        if line_index==1
            successor = ℒ[line_index+1]
            check_intersection_and_update!(t.r, X, list, ℓ, successor, r)
        elseif line_index==ℒ.count
            predecessor = ℒ[line_index-1]
            check_intersection_and_update!(t.r, X, list, ℓ, predecessor, r)
        else
            successor = ℒ[line_index+1]
            predecessor = ℒ[line_index-1]
            check_intersection_and_update!(t.r, X, list, ℓ, successor, r)
            check_intersection_and_update!(t.r, X, list, ℓ, predecessor, r)
        end
    end
end

function handle_endpoint!(t::SweepLine, X::AVLTree{Event}, ℒ::AVLTree{RootedLine}, list, e::Event, r::Real)
    point = e.point
    ℓ = e.assets
    t.r = point[1]    #update the sweep line
    line_index = sorted_rank(ℒ, ℓ)

    if line_index!=1 && line_index!=ℒ.count
        successor = ℒ[line_index+1]
        predecessor = ℒ[line_index-1]
        check_intersection_and_update!(t.r, X, list, successor, predecessor, r)
    end
    delete!(ℒ, ℓ)
end

function handle_intersection!(t::SweepLine, X::AVLTree{Event}, ℒ::AVLTree{RootedLine}, list, e::Event, r::Real)
    point = e.point
    (lower, upper) = e.assets
    (lower_rank, upper_rank) = [sorted_rank(ℒ, lower), sorted_rank(ℒ, upper)] #simplify ?

    if lower_rank!=1
        check_intersection_and_update!(point[1], X, list, upper, ℒ[lower_rank-1], r)
    end

    if upper_rank!=ℒ.count
        check_intersection_and_update!(point[1], X, list, lower, ℒ[upper_rank+1], r)
    end

    delete!(ℒ, lower) ; delete!(ℒ, upper)
    next_future_event = X[1]
    mid = 0.5 * (point[1] + next_future_event.point[1])
    t.r = mid
    insert!(ℒ, upper) ; insert!(ℒ, lower)

end

#########################################################################
# -------------------MAIN FUNCTION --------------------------------------
#########################################################################

function find_all_intersections(H::Array{Line} ; r::Real = 1.)

    @assert r>0 throw(ArgumentError("Radius r must me >0."))

    R, X, t = initialize_event_queue(H, r)
    list = Dict{Set{RootedLine}, Array{Float64}}()
    Y = AVLTree{RootedLine}()

    while X.count > 0

        event = X[1]
        delete!(X, event)

        if event.nature == "startpoint"
            handle_startpoint!(t, X, Y, list, event, r)
        elseif event.nature == "endpoint"
            handle_endpoint!(t, X, Y, list, event, r)
        else
            handle_intersection!(t, X, Y, list, event, r)
        end
    end

    return list
end

function intersections_naive(H::Array{Line}; r=1.)
    """
    Loops all over the pairs of distinct lines,
    yielding a O(n^2) algorithm.
    """
    s = size(H)[1]
    k = Int(s*(s-1)/2)
    res = Set{Array{Float64, 1}}()
    for i in 1:s
        for j in 1:(i-1)
            p = H[i] ∩ H[j]
            if normof(p)<r
                push!(res, p)
            end
        end
    end

    return res
end

# H = [Line(0.7544959394550933, 0.15515333752706018)
# Line(1.6832151948996787, 0.04401847448027585)
# Line(5.698726857470466, 0.07601306086085136)
# Line(4.232076668938423, 0.3089918936601266)]

# #find_all_intersections(H[1:2])

# H = hyperplanes_poisson(1, 5)

# include("generator.jl")
# @time a = find_all_intersections(H)
# @time b = intersections_naive(H)
# length(a)==length(b)

max_stars_repo_name: SimonCoste/LineIntersections

max_issues_repo_name: SimonCoste/LineIntersections

max_forks_repo_name: SimonCoste/LineIntersections


text:
# This code shows how visualisations created with [Luxor.jl](https://github.com/JuliaGraphics/Luxor.jl) can be
# rendered easily with MiniFB.

# Add the Luxor, Colors and MiniFB packages to your environment before running this code.
using Luxor, MiniFB, Colors

# [Guilloché](https://en.wikipedia.org/wiki/Guilloch%C3%A9) is a decorative technique in which a very precise, intricate and repetitive pattern is mechanically
# engraved into an underlying material via engine turning. The following struct contains the metadata used to render the pattern.
mutable struct Guilloche
    f1::Float64
    f2::Float64
    f3::Float64
    f4::Float64
    # damping factors: increase to make the energy losses faster
    xd1::Float64
    xd2::Float64
    yd1::Float64
    yd2::Float64
    # phase factors
    xp1::Float64
    xp2::Float64
    yp3::Float64
    yp4::Float64
end

# A function that draws a single frame of the Guilloché, using Luxor.jl functions
function drawguilloche(g::Guilloche, radius=20;
        npoints=100, limit=1000, stepby=0.025, colorfunction=(t) -> setopacity(1.0))
    function nextposition(t)
        x = exp(-g.xd1 * t) * (radius * sin((g.f1 * t) + g.xp1)) + exp(-g.xd2 * t) * (radius * sin((g.f2 * t) + g.xp2))
        y = exp(-g.yd1 * t) * (radius * cos((g.f3 * t) + g.yp3)) + exp(-g.yd2 * t) * (radius * cos((g.f4 * t) + g.yp4))
        return (x, y)
    end
    counter = 0
    move(nextposition(0)...)
    for t in 0:stepby:npoints
        colorfunction(t)
        line(nextposition(t)...)
        counter += 1
        if counter > limit
            counter = 0
            strokepath()
            colorfunction(t)
            move(nextposition(t)...)
        end
    end
    strokepath()
end

# The main method sets up the Guilloche struct, opens the window for display, and animates the frames using a while loop.
function main()
    WIDTH  = 600
    HEIGHT = 600
    window = mfb_open_ex("Luxor -> MiniFB", WIDTH, HEIGHT, MiniFB.WF_RESIZABLE)
    y = 0.1
    g2 = Guilloche(2, 1.5, 3, 1,
        0.0005, 0.0005, 0.0005, 0.0005,
        .1, .1, .1, .1)
    t1 = time_ns()
    while true
# The Luxor `@imagematrix` macro creates a 2D matrix of Color objects in ARGB32 format, which is exactly what MiniFB requires
        buffer = @imagematrix begin
                    background("grey10")
                    fontsize(20)
                    sethue("orange")
                    g2.xd1 = y/70
                    g2.xd2 = y/80
                    g2.yd1 = y/90
                    g2.yd2 = y/100
                    g2.xp1 = 8y
                    g2.xp2 = 4y
                    g2.yp3 = 1y
                    g2.yp4 = 3y
                    drawguilloche(g2, 150,
                        npoints=300,
                        stepby=0.05,
                        limit=10,
                        colorfunction = (t) -> sethue(HSL(rescale(t, 0, 300, 0, 359), 0.9, 0.7)))
                    t2 = time_ns()
                    text(string("fps: ", round(1E9/(t2-t1), digits=2)),
                        boxtopleft(BoundingBox() * 0.9))
                end 600 600
        y = y + 0.001
        (y > 1) && (y = 0.1)
# The image matrix created by Luxor is, in standard Julia convention, column first. The buffer
# expected by MiniFB is row first. The call to `permutedims` makes that change.
        state = mfb_update(window, permutedims(buffer, (2, 1)))
        t1 = t2
# Exit when the user closes the window.
        if state != MiniFB.STATE_OK
            break
        end
    end
    mfb_close(window)
end

# Call the main function. This produces really attractive animations -- I highly reccomend you run this yourself!
# The static image below shows you one frame, but does not do justice to the whole animation.
main()

# ![](../assets/guilloche.png)

# _Thanks to [`@Cormullion`](https://github.com/cormullion) for this example_


text:
struct ArrowVisualizer{V<:AbstractVisualizer}
    vis::V
end

function setobject!(vis::ArrowVisualizer, material::AbstractMaterial=defaultmaterial();
        shaft_material::AbstractMaterial=material,
        head_material::AbstractMaterial=material)
    settransform!(vis, zero(Point{3, Float64}), zero(Vec{3, Float64}))
    shaft = Cylinder(zero(Point{3, Float64}), Point(0.0, 0.0, 1.0), 1.0)
    setobject!(vis.vis[:shaft], shaft, shaft_material)
    head = Cone(zero(Point{3, Float64}), Point(0.0, 0.0, 1.0), 1.0)
    setobject!(vis.vis[:head], head, head_material)
    vis
end

function Base.show(io::IO, v::ArrowVisualizer)
    print(io, "MeshCat ArrowVisualizer with path $(v.vis.path)/arrow")
end

function settransform!(vis::ArrowVisualizer, base::Point{3}, vec::Vec{3};
        shaft_radius=0.01,
        max_head_radius=2*shaft_radius,
        max_head_length=max_head_radius)
    vec_length = norm(vec)
    T = typeof(vec_length)
    rotation = if vec_length > eps(T)
        rotation_between(SVector(0., 0., 1.), vec)
    else
        one(Quat{Float64})
    end |> LinearMap

    vis_tform = Translation(base) ∘ rotation
    settransform!(vis.vis, vis_tform)

    shaft_length = max(vec_length - max_head_length, 0)
    shaft_scaling_diag = SVector(shaft_radius, shaft_radius, shaft_length)
    if iszero(shaft_length)
        # This case is necessary to ensure that the shaft
        # completely disappears in animations.
        shaft_scaling_diag = zero(shaft_scaling_diag)
    end
    shaft_scaling = LinearMap(Diagonal(shaft_scaling_diag))
    settransform!(vis.vis[:shaft], shaft_scaling)

    head_length = vec_length - shaft_length
    head_radius = max_head_radius * head_length / max_head_length
    head_scaling = LinearMap(Diagonal(SVector(head_radius, head_radius, head_length)))
    head_tform = Translation(shaft_length * Vec(0, 0, 1)) ∘ head_scaling
    # head_tform = Translation(base) ∘ rotation ∘ Translation(shaft_length * Vec(0, 0, 1)) ∘ head_scaling
    settransform!(vis.vis[:head], head_tform)

    vis
end

function settransform!(vis::ArrowVisualizer, base::Point{3}, tip::Point{3}; kwargs...)
    settransform!(vis, base, Vec(tip - base))
end


text:
module Polytopes

using Gridap
using Gridap.Helpers
using StaticArrays
using Base.Cartesian
using LinearAlgebra

using Combinatorics

export Polytope
export NFace
export HEX_AXIS, TET_AXIS
export num_nfaces
export face_normals

export space_dim
export dim
export anchor
export extrusion
export nfaces
export nfaces_dim
export nface_connections
export nf_nfs
export nf_dim

export num_nfaces
export nface_ref_polytopes
export generate_admissible_permutations
export equidistant_interior_nodes_coordinates
export vertices_coordinates
export face_normals
export edge_tangents
export ref_nface_polytope
export nfaces_vertices

# Module constants

const HEX_AXIS = 1
const TET_AXIS = 2

# Concrete implementations

"""
n-face of the polytope, i.e., any polytope of lower dimension `N` representing
its boundary and the polytope itself (for `N` equal to the space dimension `D`)
"""
struct NFace{D}
  anchor::Point{D,Int}
  extrusion::Point{D,Int}
end

"""
Aggregation of all n-faces that compose the polytope boundary and the polytope
itself, the classification of n-faces with respect to their dimension and type
"""
struct Polytope{D}
  extrusion::Point{D,Int}
  nfaces::Vector{NFace}
  nf_nfs::Vector{Vector{Int64}}
  nf_dim::Vector{Vector{UnitRange{Int64}}}
end

"""
Constructs a `Polytope` given the type of extrusion, i.e., HEX_AXIS (=1) for "hex" extrusion
and TET_AXIS (=2) for "tet" extrusion
"""
function Polytope(extrusion::Vararg{Int,N}) where N
  return Polytope(Point{N,Int}(extrusion))
end

function Polytope(extrusion::NTuple{N,Int}) where N
  return Polytope(extrusion...)
end

function Polytope(extrusion::Point{D,Int}) where D
  zerop = Point{D,Int}(zeros(Int64, D))
  pol_nfs_dim = _polytopenfaces(zerop, extrusion)
  pol_nfs = pol_nfs_dim[1]
  pol_dim = pol_nfs_dim[2]
  nfs_id = Dict(nf => i for (i, nf) in enumerate(pol_nfs))
  nf_nfs_dim = _polytopemesh(pol_nfs, nfs_id)
  nf_nfs = nf_nfs_dim[1]
  nf_dim = nf_nfs_dim[2]
  Polytope{D}(extrusion, pol_nfs, nf_nfs, nf_dim)
end

"""
Provides the dimension of the environment space in which the n-face is defined
"""
space_dim(::NFace{D}) where {D} = D

anchor(nf::NFace) = nf.anchor

"""
Provides the dimension of the n-face, i.e., n
"""
dim(nf::NFace) = _nfdim(nf.extrusion)

extrusion(nf::NFace) = nf.extrusion

dim(::Polytope{D}) where {D} = D

extrusion(p::Polytope) = p.extrusion

nfaces(p::Polytope) = p.nfaces

nfaces_dim(p::Polytope, d) = p.nf_dim[end][d+1]

"""
# It provides for every df-face in the polytope all its dt-faces
# We use dim-wise numbering, i.e., we start numbering from 1 at every dim
"""
nface_connections(p::Polytope, dfrom, dto) = _dimfrom_fs_dimto_fs(p, dfrom, dto)

# @santiagobadia : I would prefer not to make public the following ones
"""
Provides the label of the n-faces that are a subset (on the boundary) or
equal to a given n-face. We use a global numbering for n-faces in the returned
array. n-faces are sorted by increasing dimension.
"""
nf_nfs(p::Polytope) = p.nf_nfs

"""
Provides for the result of the `nf_nfs` method, i.e., n-faces in an n-face,
the range of n-faces for every dimension lower or equal than the n-face dimension.
"""
nf_dim(p::Polytope) = p.nf_dim

"""
Returns number of nfaces of dimension dim
"""
function num_nfaces(polytope::Polytope, dim::Integer)
  n = length(polytope.nf_dim)
  k = 0
  for nface = 1:n
    d = length(polytope.nf_dim[nface]) - 1
    if d == dim
      k += 1
    end
  end
  k
end

"""
Returns number of nfaces
"""
num_nfaces(polytope::Polytope) = length(polytope.nfaces)

"""
# Returns an array with the reference polytopes for all n-faces (undef for vertices)
"""
function nface_ref_polytopes(p::Polytope)
  function _eliminate_zeros(a)
    b = Int[]
    for m in a
      if (m != 0)
        push!(b, m)
      end
    end
    return Tuple(b)
  end
  nf_ref_p = Vector{Polytope}(undef, length(p.nfaces))
  ref_nf_ps = Polytope[]
  v = _vertex()
  for (i_nf, nf) in enumerate(p.nfaces)
    r_ext = _eliminate_zeros(nf.extrusion)
    if r_ext != ()
      k = 0
      for (i_p, ref_p) in enumerate(ref_nf_ps)
        if r_ext == ref_p.extrusion
          k = i_p
          nf_ref_p[i_nf] = ref_p
        end
      end
      if k == 0
        ref_p = Polytope(r_ext)
        push!(ref_nf_ps, ref_p)
        k = length(ref_nf_ps) + 1
        nf_ref_p[i_nf] = ref_p
      end
    else
        nf_ref_p[i_nf] = v
    end
  end
  return nf_ref_p
end

"""
It generates all the admissible permutations of nodes that lead to an
admissible polytope
"""
function generate_admissible_permutations(p::Polytope)
  p_dims = length(p.extrusion)
  p_vs = Gridap.Polytopes._dimfrom_fs_dimto_fs(p, p_dims, 0)
  vs = p.nfaces[p_vs...]
  num_vs = length(vs)
  ext = p.extrusion
  l = [i for i = 1:num_vs]
  # @santiagobadia : Here we have to decide how we want this info stored
  permuted_polytopes = Vector{Int}[]
  for c in Combinatorics.permutations(l, p_dims + 1)
    admissible_polytope = true
    c1 = vs[c[1]].anchor
    for j = 2:p_dims+1
      c2 = vs[c[j]].anchor
      if (!_are_nodes_connected(c1, c2, ext))
        admissible_polytope = false
      end
    end
    if (admissible_polytope)
      push!(permuted_polytopes, c)
    end
  end
  return permuted_polytopes
end

"""
It generates the set of nodes (its coordinates) in the interior of an n-face,
for a given order. The node coordinates are the ones for a equispace case.
"""
function equidistant_interior_nodes_coordinates(p::Polytope{D}, order::Int) where D
  _order = order*ones(Int64,dim(p))
  equidistant_interior_nodes_coordinates(p, _order)
end

function equidistant_interior_nodes_coordinates(p::Polytope{D}, order::Vector{Int}) where D

  if (all(extrusion(p).array .== HEX_AXIS) || all(order .== order[1]))
    ns = _interior_nodes_int_coords(p, order)
    return ns_float = _interior_nodes_int_to_real_coords(ns, order)
  else
    error("One can consider anisotropic orders on n-cubes only")
  end
end

"""
It generates the list of coordinates of all vertices in the polytope. It is
assumed that the polytope has the bounding box [0,1]**dim
"""
function vertices_coordinates(p::Polytope{D}) where D
  vs = _dimfrom_fs_dimto_fs(p, D, 0)[1]
  vcs = Point{D,Float64}[]
  for i = 1:length(vs)
    cs = convert(Vector{Float64}, [p.nfaces[vs[i]].anchor...])
    push!(vcs, cs)
  end
  return vcs
end

"""
It generates the outwards normals of the facets of a polytope. It returns two
arrays, the first one being the outward normal and the second one the orientation.
"""
function face_normals(p::Polytope{D}) where D
  nf_vs = _dimfrom_fs_dimto_fs(p, D - 1, 0)
  vs = vertices_coordinates(p)
  f_ns = Point{D,Float64}[]
  f_os = Int[]
  for i_f = 1:length(p.nf_dim[end][end-1])
    n, f_o = _facet_normal(p, nf_vs, vs, i_f)
    push!(f_ns, Point{D,Float64}(n))
    push!(f_os, f_o)
  end
  return f_ns, f_os
end

"""
It generates the tangent vectors for polytope edges.
"""
function edge_tangents(p::Polytope{D}) where D
  ed_vs = nfaces_vertices(p,1)
  return ts = [(t = vs[2]-vs[1])/norm(t) for vs in ed_vs ]
end

"""
# Returns the reference polytope for n-faces of a given dimension.
"""
function ref_nface_polytope(p,nf_dim)
  nfs = nfaces_dim(p,nf_dim)
  fps = nface_ref_polytopes(p)[nfs]
  @assert(all(extrusion.(fps) .== extrusion(fps[1])), "All n-faces must be of the same type")
  return fps[1]
end

"""
# Return the n-faces vertices coordinates array for a given n-face dimension
"""
function nfaces_vertices(p,d)
  nc = num_nfaces(p,d)
  verts = vertices_coordinates(p)
  faces_vs = nface_connections(p,d,0)
  fvs = Gridap.CellValuesGallery.CellValueFromArray(faces_vs)
  vs = Gridap.CellValuesGallery.CellValueFromArray(verts)
  cfvs = Gridap.CellValuesGallery.CellVectorFromLocalToGlobal(fvs,vs)
end

# Helpers

_nfdim(a::Point{D,Int}) where D = sum([a[i] > 0 ? 1 : 0 for i = 1:D])

# Generates the array of n-faces of a polytope
function _polytopenfaces(anchor::Point{D,Int}, extrusion::Point{D,Int}) where D
  dnf = _nfdim(extrusion)
  zerop = Point{D,Int}(zeros(Int64, D))
  nf_nfs = []
  nf_nfs = _nfaceboundary!(anchor, zerop, extrusion, true, nf_nfs)
  [sort!(nf_nfs, by = x -> x.anchor[i]) for i = 1:length(extrusion)]
  [sort!(nf_nfs, by = x -> x.extrusion[i]) for i = 1:length(extrusion)]
  [sort!(nf_nfs, by = x -> sum(x.extrusion))]
  numnfs = length(nf_nfs)
  nfsdim = [_nfdim(nf_nfs[i].extrusion) for i = 1:numnfs]
  dimnfs = Array{UnitRange{Int64},1}(undef, dnf + 1)
  dim = 0
  i = 1
  for iface = 1:numnfs
    if (nfsdim[iface] > dim)
      # global dim; # global i
      dim += 1
      dimnfs[dim] = i:iface-1
      i = iface
    end
  end
  dimnfs[dnf+1] = numnfs:numnfs
  return [nf_nfs, dimnfs]
end

# Provides for all n-faces of a polytope the d-faces for 0 <= d <n on its
# boundary (e.g., given a face, it provides the set of edges and corners on its
# boundary) using the global n-face numbering of the base polytope
function _polytopemesh(nfaces::Vector{NFace{D}}, nfaceid::Dict) where D
  num_nfs = length(nfaces)
  nfnfs = Vector{Vector{Int64}}(undef, num_nfs)
  nfnfs_dim = Vector{Vector{UnitRange{Int64}}}(undef, num_nfs)
  for (inf, nf) in enumerate(nfaces)
    nfs_dim_nf = _polytopenfaces(nf.anchor, nf.extrusion)
    nf_nfs = nfs_dim_nf[1]
    dimnfs = nfs_dim_nf[2]
    nfnfs[inf] = [get(nfaceid, nf, nf) for nf in nf_nfs]
    nfnfs_dim[inf] = dimnfs
  end
  return [nfnfs, nfnfs_dim]
end

# Generates the list of n-face of a polytope the d-faces for 0 <= d <n on its
# boundary
function _nfaceboundary!(
  anchor::Point{D,Int},
  extrusion::Point{D,Int},
  extend::Point{D,Int},
  isanchor::Bool,
  list
) where D
  newext = extend
  list = [list..., NFace{D}(anchor, extrusion)]
  for i = 1:D
    curex = newext[i]
    if (curex > 0) # Perform extension
      func1 = (j -> j == i ? 0 : newext[j])
      newext = Point{D,Int}([func1(i) for i = 1:D])
      func2 = (j -> j == i ? 1 : 0)
      edim = Point{D,Int}([func2(i) for i = 1:D])
      func3 = (j -> j >= i ? anchor[j] : 0)
      tetp = Point{D,Int}([func3(i) for i = 1:D]) + edim
      if (curex == 1) # Quad extension
        list = _nfaceboundary!(anchor + edim, extrusion, newext, false, list)
      elseif (isanchor)
        list = _nfaceboundary!(tetp, extrusion, newext, false, list)
      end
      list = _nfaceboundary!(
        anchor,
        extrusion + edim * curex,
        newext,
        false,
        list
      )
    end
  end
  return list
end

function _dimfrom_fs_dimto_fs(p::Polytope, dim_from::Int, dim_to::Int)
  @assert dim_to <= dim_from
  dim_from += 1
  dim_to += 1
  dffs_r = p.nf_dim[end][dim_from]
  dffs_dtfs = Vector{Vector{Int}}(undef, dffs_r[end] - dffs_r[1] + 1)
  offs = p.nf_dim[end][dim_to][1] - 1
  for (i_dff, dff) in enumerate(dffs_r)
    dff_nfs = p.nf_nfs[dff]
    dff_dtfs_r = p.nf_dim[dff][dim_to]
    dff_dtfs = dff_nfs[dff_dtfs_r]
    dffs_dtfs[i_dff] = dff_dtfs .- offs
    # @santiagobadia : With or without offset ?
  end
  return dffs_dtfs
end

# Auxiliary function that determines whether two nodes are connected
function _are_nodes_connected(c1, c2, ext)
  sp_dims = length(c1)
  d = zeros(length(c1))
  for i = 1:length(d)
    d[i] = c2[i] - c1[i]
  end
  dn = sum(d .* d)
  connected = false
  if (dn == 1)
    connected = true
  else
    k = 0
    for j = 1:sp_dims
      if (ext[j] == 2)
        if (c1[j] == 1 || c2[j] == 1)
          k = j
        end
      end
    end
    for l = 1:k
      d[l] = 0
    end
    dn = sum(d .* d)
    if (dn == 0)
      connected = true
    end
  end
  return connected
end

# It generates the set of nodes (its coordinates) in the interior of an n-face,
# for a given order. The node coordinates are `Int` and from 0 to `order` per
# direction
function _interior_nodes_int_coords(p::Polytope{D}, order) where D
  ext = p.extrusion
  _ord = [order...]
  verts = Point{D,Int}[]
  coor = zeros(Int, D)
  _generate_nodes!(D, p.extrusion, _ord, coor, verts)
  return verts
end

# Auxiliary private recursive function to implement _interior_nodes_int_coords
function _generate_nodes!(dim, ext, order, coor, verts)
  ncoo = copy(coor)
  nord = copy(order)
  for i = 1:order[dim]-1
    ncoo[dim] = i
    if dim > 1
      if (ext[dim] == TET_AXIS)
        nord .-= 1
      end
      _generate_nodes!(dim - 1, ext, nord, ncoo, verts)
    else
      push!(verts, Point(ncoo...))
    end
  end
end

# Transforms the int coordinates to float coordinates
function _interior_nodes_int_to_real_coords(nodes, order)
  if length(nodes) > 0
    dim = length(nodes[1])
    cs_float = Point{dim,Float64}[]
    cs = zeros(Float64, dim)
    for cs_int in nodes
      for i = 1:dim
        cs[i] = cs_int[i] / order[i]
      end
      push!(cs_float, Point{dim,Float64}(cs))
    end
  else
    cs_float = Point{0,Float64}[]
  end
  return cs_float
end

function _facet_normal(p::Polytope{D}, nf_vs, vs, i_f) where D
  if (length(p.extrusion) > 1)
    v = Float64[]
    for i = 2:length(nf_vs[i_f])
      vi = vs[nf_vs[i_f][i]] - vs[nf_vs[i_f][1]]
      push!(v, vi...)
    end
    n = nullspace(transpose(reshape(v, D, length(nf_vs[i_f]) - 1)))
    n = n .* 1 / sqrt(dot(n, n))
    ext_v = _vertex_not_in_facet(p, i_f, nf_vs)
    v3 = vs[nf_vs[i_f][1]] - vs[ext_v]
    f_or = 1
    if dot(v3, n) < 0.0
      n *= -1
      f_or = -1
    end
  elseif (length(p.extrusion) == 1)
    ext_v = _vertex_not_in_facet(p, i_f, nf_vs)
    n = vs[nf_vs[i_f][1]] - vs[ext_v]
    n = n .* 1 / dot(n, n)
    f_or = 1
  else
    error("O-dim polytopes do not have properly define outward facet normals")
  end
  return n, f_or
end

function _vertex_not_in_facet(p, i_f, nf_vs)
  for i in p.nf_dim[end][1]
    is_in_f = false
    for j in nf_vs[i_f]
      if i == j
        is_in_f = true
        break
      end
    end
    if !is_in_f
      return i
      break
    end
  end
end

# Generates a zero-dim polytope (vertex)
function _vertex()
  ext = ()
  nfdim = [[1:1]]
  nfnfs = [[1]]
  nfanc = Point{0,Int}()
  nf = NFace{0}(nfanc,nfanc)
  nfs = [nf]
  return Polytope{0}(ext, nfs, nfnfs, nfdim)
end

end # module Polytopes


text:
using FractionalDiffEq, Plots
fun(x, y) = 1-y
prob = SingleTermFODEProblem(fun, 0.5, 0, (0, 5))
sol = solve(prob, 0.01, PECE())
tspan = collect(0:0.01:5)
#plot(result.t, result.u)
plot(sol, linewidth=2, legend=:bottomright)


text:
abstract type SamplingOperator <: AbstractOperator
end

dest_space(op::SamplingOperator) = Span(dest(op))

function apply(op::SamplingOperator, f::AbstractVector; options...)
	@warn "This is a funny function to use"
	@assert length(f) == size(op,2)
	f
end

"""
A `GridSampling` is an operator that maps a function to its samples
in a grid.
"""
struct GridSampling <: SamplingOperator
    src     ::  FunctionSpace
    dest    ::  GridBasis
end

GridSampling(grid::AbstractGrid{S}, ::Type{T} = subeltype(S)) where {S,T} =
    GridSampling(GenericFunctionSpace{S,T}(), GridBasis{T}(grid))

GridSampling(gridbasis::GridBasis{S,T}) where {S,T} =
	GridSampling(GenericFunctionSpace{eltype(grid(gridbasis)),T}(), gridbasis)

GridSampling(m::DiscreteWeight, ::Type{T} = numeltype(m)) where T = GridSampling(points(m), T)

grid(op::GridSampling) = grid(dest(op))

(op::GridSampling)(f) = apply(op, f)

dest(op::GridSampling) = op.dest

src_space(op::GridSampling) = op.src

apply!(result, op::GridSampling, f; options...) = sample!(result, grid(op), f; options...)

"Sample the function f on the given grid."
sample(g::AbstractGrid, f, T = numtype(g)) = sample!(zeros(T, size(g)), g, f)

broadcast(f::Function, grid::AbstractGrid) = sample(grid, f)

component(op::GridSampling, i) = GridSampling(component(dest(op),i))
components(op::GridSampling) = map( s -> GridSampling(s), components(dest(op)))

factors(op::GridSampling) = components(op)

function tensorproduct(op1::GridSampling, op2::GridSampling)
	T = promote_type(coefficienttype(dest(op1)),coefficienttype(dest(op2)))
	grid1 = grid(op1)
	grid2 = grid(op2)
	GridSampling(grid1 × grid2, T)
end

function tensorproduct(op1::GridSampling, op2::AbstractOperator)
	@assert dest(op2) isa GridBasis
	@assert ncomponents(op2) == 2
	tensorproduct(IdentityOperator(dest(op1)), component(op2,2)) * tensorproduct(op1, component(op2,1))
end

function tensorproduct(ops::GridSampling...)
	T = promote_type(map(t->coefficienttype(dest(t)), ops)...)
	g = ProductGrid(map(grid, ops)...)
	GridSampling(g, T)
end

function tensorproduct(op1::AbstractOperator, op2::AbstractOperator)
	@assert dest(op1) isa GridBasis
	@assert dest(op2) isa GridBasis
	@assert ncomponents(op1) == 2
	@assert ncomponents(op2) == 2
	T = promote_type(coefficienttype(dest(op1)),coefficienttype(dest(op2)))
	grid1 = grid(dest(op1))
	grid2 = grid(dest(op2))
	g = grid1 × grid2
	tensorproduct(component(op1,2),component(op2,2)) * tensorproduct(component(op1,1),component(op2,1))
end

function tensorproduct(op1::AbstractOperator, op2::AbstractOperator, op3::AbstractOperator)
	@assert dest(op1) isa GridBasis
	@assert dest(op2) isa GridBasis
	@assert dest(op3) isa GridBasis
	# TODO: generalize to longer operators
	@assert ncomponents(op1) == 2
	@assert ncomponents(op2) == 2
	@assert ncomponents(op3) == 2
	tensorproduct(component(op1,2),component(op2,2),component(op3,2)) * tensorproduct(component(op1,1),component(op2,1),component(op3,1))
end

# We don't want to assume that f can be called with a vector argument.
# In order to avoid the overhead of splatting, we capture a number of special cases
call_function_with_vector(f, x::Number) = f(x)
call_function_with_vector(f, x::SVector{1}) = f(x[1])
call_function_with_vector(f, x::SVector{2}) = f(x[1], x[2])
call_function_with_vector(f, x::SVector{3}) = f(x[1], x[2], x[3])
call_function_with_vector(f, x::SVector{4}) = f(x[1], x[2], x[3], x[4])
call_function_with_vector(f, x::SVector{N}) where {N} = f(x...)
call_function_with_vector(f, x::AbstractVector) = f(x...)

function sample!(result, grid, f; options...)
    for i in eachindex(grid)
		result[i] = call_function_with_vector(f, grid[i])
	end
	result
end

show(io::IO, mime::MIME"text/plain", op::GridSampling) = composite_show(io, mime, op)
Display.object_parentheses(op::GridSampling) = false
Display.stencil_parentheses(op::GridSampling) = false
Display.displaystencil(op::GridSampling) = ["GridSampling(", src_space(op), ", ", grid(op), ")"]

"""
A `ProjectionSampling` is an operator that maps a function to its inner products
with a projection basis.
"""
struct ProjectionSampling <: SamplingOperator
    dict		::  Dictionary
	measure		::	Measure
	src_space	::	FunctionSpace
end

const AnalysisOperator = ProjectionSampling

space(dict::Dictionary) = space(measure(dict))

ProjectionSampling(dict::Dictionary) = ProjectionSampling(dict, measure(dict), Span(dict))

ProjectionSampling(dict::Dictionary, measure::Measure) = ProjectionSampling(dict, measure, space(measure))

ProjectionSampling(dict::Dictionary, space::FunctionSpace) = ProjectionSampling(dict, measure(space), space)

dictionary(op::ProjectionSampling) = op.dict

dest(op::ProjectionSampling) = dictionary(op)

src_space(op::ProjectionSampling) = op.src_space

measure(op::ProjectionSampling) = op.measure

apply!(result, op::ProjectionSampling, f; options...) = project!(result, f, dictionary(op), measure(op); options...)

function project!(result, f, dict::Dictionary, measure::Measure; options...)
    for i in eachindex(result)
		result[i] = innerproduct(dict[i], f, measure; options...)
	end
	result
end

show(io::IO, mime::MIME"text/plain", op::ProjectionSampling) = composite_show(io, mime, op)
Display.object_parentheses(op::ProjectionSampling) = false
Display.stencil_parentheses(op::ProjectionSampling) = false
Display.displaystencil(op::ProjectionSampling) = ["ProjectionSampling(", dictionary(op), ", ", measure(op), ")"]


text:
module MaxHelpingHandSnowCustomColors

using ..Ahorn, Maple

@mapdef Effect "MaxHelpingHand/SnowCustomColors" SnowCustomColors(only::String="*", exclude::String="", colors::String="FF0000,00FF00,0000FF")

placements = SnowCustomColors

function Ahorn.canFgBg(effect::SnowCustomColors)
    return true, true
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
function selectionsort!(arr::Vector{<:Real})
    len = length(arr)
    if len < 2 return arr end
    for i in 1:len-1
        lmin, j = findmin(arr[i+1:end])
        if lmin < arr[i]
            arr[i+j] = arr[i]
            arr[i] = lmin
        end
    end
    return arr
end

v = rand(-10:10, 10)
println("# unordered: $v\n -> ordered: ", selectionsort!(v))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#####
##### Convinient aliases
#####

const AG  = AbstractGrid
const RCG = RegularCartesianGrid

#####
##### Base interpolation operators
#####

@inline ℑxᶜᵃᵃ(i, j, k, grid::AG{FT}, u) where FT = @inbounds FT(0.5) * (u[i,   j, k] + u[i+1, j, k])
@inline ℑxᶠᵃᵃ(i, j, k, grid::AG{FT}, c) where FT = @inbounds FT(0.5) * (c[i-1, j, k] + c[i,   j, k])

@inline ℑyᵃᶜᵃ(i, j, k, grid::AG{FT}, v) where FT = @inbounds FT(0.5) * (v[i, j,   k] + v[i,  j+1, k])
@inline ℑyᵃᶠᵃ(i, j, k, grid::AG{FT}, c) where FT = @inbounds FT(0.5) * (c[i, j-1, k] + c[i,  j,   k])

@inline ℑzᵃᵃᶜ(i, j, k, grid::AG{FT}, w) where FT = @inbounds FT(0.5) * (w[i, j,   k] + w[i, j, k+1])
@inline ℑzᵃᵃᶠ(i, j, k, grid::AG{FT}, c) where FT = @inbounds FT(0.5) * (c[i, j, k-1] + c[i, j,   k])

#####
##### Interpolation operators acting on functions
#####

@inline ℑxᶜᵃᵃ(i, j, k, grid::AG{FT}, f::F, args...) where {FT, F<:Function} = FT(0.5) * (f(i,   j, k, grid, args...) + f(i+1, j, k, grid, args...))
@inline ℑxᶠᵃᵃ(i, j, k, grid::AG{FT}, f::F, args...) where {FT, F<:Function} = FT(0.5) * (f(i-1, j, k, grid, args...) + f(i,   j, k, grid, args...))

@inline ℑyᵃᶜᵃ(i, j, k, grid::AG{FT}, f::F, args...) where {FT, F<:Function} = FT(0.5) * (f(i, j,   k, grid, args...) + f(i, j+1, k, grid, args...))
@inline ℑyᵃᶠᵃ(i, j, k, grid::AG{FT}, f::F, args...) where {FT, F<:Function} = FT(0.5) * (f(i, j-1, k, grid, args...) + f(i, j,   k, grid, args...))

@inline ℑzᵃᵃᶜ(i, j, k, grid::AG{FT}, f::F, args...) where {FT, F<:Function} = FT(0.5) * (f(i, j, k,   grid, args...) + f(i, j, k+1, grid, args...))
@inline ℑzᵃᵃᶠ(i, j, k, grid::AG{FT}, f::F, args...) where {FT, F<:Function} = FT(0.5) * (f(i, j, k-1, grid, args...) + f(i, j, k,   grid, args...))

#####
##### Convenience operators for "interpolating constants"
#####

@inline ℑxᶠᵃᵃ(i, j, k, grid, f::Number, args...) = f
@inline ℑxᶜᵃᵃ(i, j, k, grid, f::Number, args...) = f
@inline ℑyᵃᶠᵃ(i, j, k, grid, f::Number, args...) = f
@inline ℑyᵃᶜᵃ(i, j, k, grid, f::Number, args...) = f
@inline ℑzᵃᵃᶠ(i, j, k, grid, f::Number, args...) = f
@inline ℑzᵃᵃᶜ(i, j, k, grid, f::Number, args...) = f

#####
##### Double interpolation
#####

@inline ℑxyᶜᶜᵃ(i, j, k, grid, f, args...) = ℑyᵃᶜᵃ(i, j, k, grid, ℑxᶜᵃᵃ, f, args...)
@inline ℑxyᶠᶜᵃ(i, j, k, grid, f, args...) = ℑyᵃᶜᵃ(i, j, k, grid, ℑxᶠᵃᵃ, f, args...)
@inline ℑxyᶠᶠᵃ(i, j, k, grid, f, args...) = ℑyᵃᶠᵃ(i, j, k, grid, ℑxᶠᵃᵃ, f, args...)
@inline ℑxyᶜᶠᵃ(i, j, k, grid, f, args...) = ℑyᵃᶠᵃ(i, j, k, grid, ℑxᶜᵃᵃ, f, args...)
@inline ℑxzᶜᵃᶜ(i, j, k, grid, f, args...) = ℑzᵃᵃᶜ(i, j, k, grid, ℑxᶜᵃᵃ, f, args...)
@inline ℑxzᶠᵃᶜ(i, j, k, grid, f, args...) = ℑzᵃᵃᶜ(i, j, k, grid, ℑxᶠᵃᵃ, f, args...)
@inline ℑxzᶠᵃᶠ(i, j, k, grid, f, args...) = ℑzᵃᵃᶠ(i, j, k, grid, ℑxᶠᵃᵃ, f, args...)
@inline ℑxzᶜᵃᶠ(i, j, k, grid, f, args...) = ℑzᵃᵃᶠ(i, j, k, grid, ℑxᶜᵃᵃ, f, args...)
@inline ℑyzᵃᶜᶜ(i, j, k, grid, f, args...) = ℑzᵃᵃᶜ(i, j, k, grid, ℑyᵃᶜᵃ, f, args...)
@inline ℑyzᵃᶠᶜ(i, j, k, grid, f, args...) = ℑzᵃᵃᶜ(i, j, k, grid, ℑyᵃᶠᵃ, f, args...)
@inline ℑyzᵃᶠᶠ(i, j, k, grid, f, args...) = ℑzᵃᵃᶠ(i, j, k, grid, ℑyᵃᶠᵃ, f, args...)
@inline ℑyzᵃᶜᶠ(i, j, k, grid, f, args...) = ℑzᵃᵃᶠ(i, j, k, grid, ℑyᵃᶜᵃ, f, args...)

#####
##### Triple interpolation
#####

@inline ℑxyzᶜᶜᶠ(i, j, k, grid, f, args...) = ℑxᶜᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, ℑzᵃᵃᶠ, f, args...)
@inline ℑxyzᶜᶠᶜ(i, j, k, grid, f, args...) = ℑxᶜᵃᵃ(i, j, k, grid, ℑyᵃᶠᵃ, ℑzᵃᵃᶜ, f, args...)
@inline ℑxyzᶠᶜᶜ(i, j, k, grid, f, args...) = ℑxᶠᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, ℑzᵃᵃᶜ, f, args...)

@inline ℑxyzᶜᶠᶠ(i, j, k, grid, f, args...) = ℑxᶜᵃᵃ(i, j, k, grid, ℑyᵃᶠᵃ, ℑzᵃᵃᶠ, f, args...)
@inline ℑxyzᶠᶜᶠ(i, j, k, grid, f, args...) = ℑxᶠᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, ℑzᵃᵃᶠ, f, args...)
@inline ℑxyzᶠᶠᶜ(i, j, k, grid, f, args...) = ℑxᶠᵃᵃ(i, j, k, grid, ℑyᵃᶠᵃ, ℑzᵃᵃᶜ, f, args...)

#####
##### Leftover
#####

@inline ℑxᶠᵃᵃ(i, j, k, grid::AG{FT}, c, args...) where FT = @inbounds FT(0.5) * (c[i-1, j, k] + c[i, j, k])
@inline ℑyᵃᶠᵃ(i, j, k, grid::AG{FT}, c, args...) where FT = @inbounds FT(0.5) * (c[i, j-1, k] + c[i, j, k])
@inline ℑzᵃᵃᶠ(i, j, k, grid::AG{FT}, c, args...) where FT = @inbounds FT(0.5) * (c[i, j, k-1] + c[i, j, k])


text:
type VideoMode
    width::UInt32
    height::UInt32
    bits_per_pixel::UInt32

    function VideoMode(width::Integer, height::Integer, bits_per_pixel::Integer)
        return new(width, height, bits_per_pixel)
    end

    function VideoMode(width::Integer, height::Integer)
        return new(width, height, 32)
    end
end

function get_desktop_mode()
    return ccall((:sfVideoMode_getDesktopMode, libcsfml_window), VideoMode, ())
end

function is_valid(mode::VideoMode)
    return ccall((:sfVideoMode_isValid, libcsfml_window), Bool, (VideoMode,), mode)
end


text:
#include("Crystal.jl")

"""
    module DFToutMod

Deal with energy, force, stress, band structure, etc, from a DFT calculation

Generic to different DFT codes.
"""
module DFToutMod

#using LinearAlgebra
using ..CrystalMod:crystal
using ..CrystalMod:makecrys
using ..CrystalMod:orbital_index
using Printf

using ..ThreeBodyTB:convert_energy
using ..ThreeBodyTB:convert_force
using ..ThreeBodyTB:convert_stress
using ..ThreeBodyTB:global_energy_units

##using Formatting

export dftout
export makedftout
#export test
export test2
#export crystal
#export makecrys
#export generate_supercell
#export generate_random
#export write_poscar
#export write_efs

#holds data
"""
    mutable struct bandstructure

Band structure. Has

- `nbnd::Int` Number of bands
- `nks::Int` Number of k-points
- `nelec::Float64` Number of electrons
- `efermi::Float64` Fermi energy
- `kpts::Array{Float64,2}` List of k-points, `nks × 3` array in BZ crystal units.
- `kweights::Array{Float64,1}` Weights of k-points (`nks`).
- `kgrid::Array{Int,1}` Equivalent gamma centered Monkhorst-Pack k-grid dimensions, if applies.
- `eigs::Array{Float64,2}` Eigenvalues. `nks × nbnd`

"""
mutable struct bandstructure
    nbnd::Int
    nks::Int
    nelec::Float64
    efermi::Float64
    kpts::Array{Float64,2}
    kweights::Array{Float64,1}
    kgrid::Array{Int,1}
    eigs::Array{Float64,2}
end

"""
    mutable struct dftout

DFT output struct. Has

- `crys::crystal`  crystal structure
- `energy::Float64` the actual DFT energy, depends on pseudopotentials.
- `energy_smear::Float64` smearing energy
- `forces::Array{Float64,2}` forces Ryd / a.u.
- `stress::Array{Float64,2}` stress  Ryd / (a.u.)^3
- `bandstruct::bandstructure` See bandstructure struct
- `hasband::Bool` does this object have a band structure, usually `true`
- `hasham::Bool` not used, always `false`
- `prefix::String`  A string has a name used to find output files.
- `outdir::String` Directly loaded from
- `tot_charge::Float64` If charge of unit cell is nonzero
- `atomize_energy::Float64` Atomization energy, relative to non-spin-polarized atoms.
"""
mutable struct dftout

    crys::crystal
    energy::Float64
    energy_smear::Float64
    forces::Array{Float64,2}
    stress::Array{Float64,2}
    bandstruct::bandstructure
    hasband::Bool
    hasham::Bool
    prefix::String
    outdir::String
    tot_charge::Float64
    atomize_energy::Float64

end

"""
    function get_atomize_energy(d::dftout)

Return atomization energy in current energy units
"""
function get_atomize_energy(d::dftout)

    return convert_energy(d.atomize_energy)

end

#printing
Base.show(io::IO, d::dftout) = begin
#    println(io)
#    for i in 1:3
#        @printf(io, "A%.1i=     % .5f  % .5f  % .5f\n", i, d.crys.A[i,1],  d.crys.A[i,2],  d.crys.A[i,3])
#    end
    println(io, d.crys)

    println(io)
    println("Energy (direct from dft) ", d.energy, " Ryd")
    println(io)
    println(io,"STRUCTURE =====================| FORCES =================")
    println(io)

    forces = convert_force(d.forces)

    for i in 1:d.crys.nat
        @printf(io, "%-3s  % .5f  % .5f  % .5f | % .5f % .5f % .5f\n", d.crys.types[i], d.crys.coords[i,1], d.crys.coords[i,2], d.crys.coords[i,3], forces[i,1], forces[i,2], forces[i,3])
    end
    println(io,)
    println(io,"=========================================================")
    println(io)
    println(io,"STRESS ==================")

    stress = convert_stress(d.stress)

    @printf(io, "\n")
    for i in 1:3
        @printf(io, "% .5f  % .5f  % .5f\n", stress[i,1], stress[i,2], stress[i,3])
    end
    println(io,)
    println(io, "Atomization energy (no spin): ", convert_energy(d.atomize_energy), " $global_energy_units ")
    println(io,"=========================")
    println(io,"Has bandstructure: ", d.hasband,"; Has hamiltonian: ", d.hasham)
    println(io)

end

#print bandstructure
Base.show(io::IO, d::bandstructure) = begin
    println(io,"nbnd = ", d.nbnd, "; nkpts = ", d.nks)
    println(io,"nelec = ", d.nelec, "; efermi = ", d.efermi)

    println(io)

#    pst = "{:d} {: f} {: f} {: f} , {: f},  "
#    pst = "%-3s % .5f % .5f % .5f , % .5f, "
#    println(pst)

    for i = 1:min(d.nks, 200)
        t = @sprintf("%-3s", i)
        k = @sprintf(" [ % .5f % .5f % .5f ]", d.kpts[i,1],d.kpts[i,2],d.kpts[i,3]  )
        kw = @sprintf(" % .5f ",  d.kweights[i])
        v = " "
        for j = 1:d.nbnd
            v = v * @sprintf(" % .5f",d.eigs[i,j])
        end
        println(io,t,kw,k,v)

    end
    if d.nks > 200
        println(io, "... [truncated]")
    end
    println(io)

end

"""
    function makebs(nelec::Number, efermi::Number,  kpoints, kweights,kgrid, vals)

Constructor for `bandstructure`.
"""
function makebs(nelec::Number, efermi::Number,  kpoints, kweights,kgrid, vals)
    try
        kpoints = convert(Array{Float64,2},kpoints)
        kweights = convert(Array{Float64,1}, kweights)
        vals = convert(Array{Float64,2}, vals)
    catch
        println(typeof(kpoints), typeof(kweights), typeof(vals))
        error("error make bs type conversion")
    end
    nks = size(kpoints,1)
    if nks != size(kweights,1) || nks != size(vals,1)
        error("kpoints kweights eigs don't match")
    end
    nbnd = size(vals,2)

    K2 = zeros(Float64,size(kpoints))
    for i in 1:nks
        K2[i,:] = Int64.(round.(kpoints[i,:] * 100000000))/100000000
    end

    return bandstructure(nbnd, nks, nelec, efermi, K2, kweights, kgrid, vals)

end

function make_empty_bs()

    return bandstructure(1, 1, 1, 0.0, zeros(1,3), zeros(3), [1, 1, 1], zeros(1,1))

end

"""
    function makedftout(crys::crystal, energy::Number, energy_smear::Number,  forces, stress, bandstruct=missing; prefix="PREFIX", outdir="TMPDIR", tot_charge=0.0)

Constructor for dftout. Usually called by function that loads DFT output files, not called directly.
"""
function makedftout(crys::crystal, energy::Number, energy_smear::Number,  forces, stress, bandstruct=missing; prefix="PREFIX", outdir="TMPDIR", tot_charge=0.0)
"""
Creates a struct with the desired data
"""
    try
        forces = convert(Array{Float64,2}, forces)
        stress = convert(Array{Float64,2}, stress)
        energy = convert(Float64, energy)
        energy_smear = convert(Float64, energy_smear)
    catch
        println(typeof(energy), typeof(forces), typeof(stress))
        error("Error converting types")
        return
    end

    nat = crys.nat
    if (nat,3) != size(forces)
        println( forces)
        error("Error forces")

        return
    end
    if size(stress) != (3,3)
        println( stress)
        error("Error stress", size(stress))
        return
    end

    ind2orb, orb2ind, etotal_atoms, nval =  orbital_index(crys)
    atomize_energy = energy - etotal_atoms

    if ismissing(bandstruct)
        return dftout(crys, energy, energy_smear, forces, stress, make_empty_bs() , false, false, prefix, outdir, tot_charge, atomize_energy) #, missing, False, False)
    else
        return dftout(crys, energy, energy_smear,forces, stress, bandstruct, true, false, prefix, outdir, tot_charge, atomize_energy) #, missing, False, False)
    end
end

#alterate version, will make the crys for you.
"""
    function makedftout(A, pos, types, energy::Number,energy_smear::Number,  forces, stress, bandstruct=missing; prefix="PREFIX", outdir="TMPDIR", tot_charge=0.0)
"""
function makedftout(A, pos, types, energy::Number,energy_smear::Number,  forces, stress, bandstruct=missing; prefix="PREFIX", outdir="TMPDIR", tot_charge=0.0)
    c = makecrys(A,pos,types, units="Bohr")
    if ismissing(bandstruct)
        return makedftout(c, energy, energy_smear, forces,stress, prefix=prefix, outdir=outdir, tot_charge=0.0)
    else
        return makedftout(c, energy, energy_smear, forces,stress, bandstruct,prefix=prefix, outdir=outdir, tot_charge=0.0)
    end
end

end #end module

#using .DFToutMod
#Base.:(==)(x::dftout, y::dftout) = (x.crys == y.crys && x.energy == y.energy && x.forces == y.forces && x.stress ==y.stress)


text:
"""
    parse_map!()
Get id, width and height from readline, parse the rest using parse_turn!
"""
function parse_map!()
    id = parse(Int, readline())
    width, height = parse.(Int, split(readline()))
    return GameMap(id, width, height, parse_turn!()...)
end
update_map(game::Game) = GameMap(game.id, game.width, game.height, parse_turn!()...)
start_game(game::Game) = print(game.botname, "\n")

"""
    parse_turn(parse_turn(turn::Vector{String}))
Parse a turn message received from halite.
See the [starter-kit-reference][1]. This [post][2] is also really helpful.
[1]: https://halite.io/learn-programming-challenge/downloads-and-starter-kits/create-new-starter-kit
[2]: https://forums.halite.io/t/halite-comms-protocol-explained/411
"""
function parse_turn!(turn::Vector{SubString{String}})
    players = Dict([parse_player!(turn) for _ in 1:parse(Int,shift!(turn))])
    planets = Dict([parse_planet!(turn) for _ in 1:parse(Int,shift!(turn))])
    @assert isempty(turn)
    return players, planets
end
parse_turn!() = parse_turn!(split(readline()))

"""
    parse_player!(tokens)
Parse a vector of `2+n` tokens as a `Player` where `n` is the number of owned ships.
Since constructors fall back to `convert` methods, parse all strings to `Float64` or `Int`.
"""
function parse_player!(tokens)
    id = parse(Int, shift!(tokens))
    nr_ships = parse(Int,shift!(tokens))
    ships = Dict([parse_ships!(id, tokens) for _ in 1:nr_ships])
    return id => Player(id, ships)
end

"""
    parse_ships!(owner_id, tokens)
Parse a vector of `10` tokens as a `Ship`.
Since constructors fall back to `convert` methods, parse all strings to `Float64` or `Int`.
"""
function parse_ships!(owner_id, tokens)
    id = parse(Int, shift!(tokens))
    id => Ship(owner_id, id, parse.(Float64, splice!(tokens, 1:9))...)
end

"""
    parse_planet!(tokens)
Parse a vector of `11+n` tokens as a `Planet` where `n` is the number of docked ships.
Since constructors fall back to `convert` methods, parse all strings to `Float64` or `Int`.
"""
function parse_planet!(tokens)
    id = parse(Int, shift!(tokens))
    parameters = parse.(Float64,splice!(tokens, 1:9))
    nr_docked_ships = parse(Int,shift!(tokens))
    ship_ids = parse.(Int,splice!(tokens, 1:nr_docked_ships))
    return id => Planet(id, parameters..., ship_ids)
end

# Constructors ommitting deprecated parameters
Ship(owner_id, id, x, y, hp, xv, yv, docked, docked_planet, progress, cooldown) =
Ship(owner_id, id, x, y, hp,    Int(docked), docked_planet, progress, cooldown)

Planet(id, x, y, hp, radius, docks, current_prod, rm_prod,  owned, owner_id, docked_ships_ids) =
Planet(id, x, y, hp, radius, docks, current_prod,           owned, owner_id, docked_ships_ids)

max_stars_repo_name: Cerbes/Halite-II

max_issues_repo_name: Cerbes/Halite-II

max_forks_repo_name: Cerbes/Halite-II


text:
Initial parent is IU RFSGJABGOLYWF XSMFXNIABKT, its fitness is 1
After 10 generations, the new parent is M TBINGJ IGOYSYV KIAM WIAXEL and its fitness is 13
After 20 generations, the new parent is MRTBINKJ IT  SYO KT Z WEAIEL and its fitness is 18
After 30 generations, the new parent is MGTHINKJ IT ISYLMKJ A WEASEL and its fitness is 23
After 40 generations, the new parent is MBTHINKS IT ISYLIKE A WEASEL and its fitness is 26
After 49 generations, the parent evolved into the target METHINKS IT IS LIKE A WEASEL


text:
# WINDOW clause.

mutable struct WindowClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}
    list::Vector{SQLClause}

    WindowClause(;
                 over = nothing,
                 list) =
        new(over, list)
end

WindowClause(list...; over = nothing) =
    WindowClause(over = over, list = SQLClause[list...])

"""
    WINDOW(; over = nothing, list)
    WINDOW(list...; over = nothing)

A `WINDOW` clause.

# Examples

```jldoctest
julia> c = FROM(:person) |>
           WINDOW(:w1 => PARTITION(:year_of_birth),
                  :w2 => :w1 |> PARTITION(order_by = [:month_of_birth, :day_of_birth])) |>
           SELECT(:person_id, AGG("ROW_NUMBER", over = :w2));

julia> print(render(c))
SELECT "person_id", (ROW_NUMBER() OVER ("w2"))
FROM "person"
WINDOW "w1" AS (PARTITION BY "year_of_birth"), "w2" AS ("w1" ORDER BY "month_of_birth", "day_of_birth")
```
"""
WINDOW(args...; kws...) =
    WindowClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(WINDOW), pats::Vector{Any}) =
    dissect(scr, WindowClause, pats)

function PrettyPrinting.quoteof(c::WindowClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(WINDOW))
    if isempty(c.list)
        push!(ex.args, Expr(:kw, :list, Expr(:vect)))
    else
        append!(ex.args, quoteof(c.list, ctx))
    end
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::WindowClause, c′) =
    WindowClause(over = rebase(c.over, c′), list = c.list)


text:
include("selinv.jl")

"""
    selinv(F::SimplicialLDLt) -> B::SparseMatrixCSC

Compute the entries of the inverse of `A = F.L*F.D*F.Lt` contained in the
sparsity pattern of the (incomplete) LDLt factorization `F`.
"""
function selinv(F::SimplicialLDLt)
    Fp,Fi,Fv,conj = F.colptr,F.rowval,F.nzval,F.conj
    Bv = selinv_from_ldlt(Fp,Fi,Fv; conj=conj)
    return SparseMatrixCSC(size(F)...,Fp,Fi,Bv)
end

max_stars_repo_name: ettersi/SparseFactorizations

max_issues_repo_name: ettersi/SparseFactorizations

max_forks_repo_name: ettersi/SparseFactorizations


text:
module ChromNetPaper

globalDataDir = normpath(joinpath(dirname(Base.source_path()), "..", "data"))

include("evaluate.jl")
include("groupgm.jl")
include("drive.jl")

end # module


text:
function periodic_orbit(
    T::Type, sys::ODESystem, n_mesh::Integer, n_coll::Integer; name=gensym(:PeriodicOrbit)
)
    collsys = collocation(T, sys, n_mesh, n_coll; name)
    return push!(collsys.equations, collsys.t₀ ~ 0)
end

function periodic_orbit(sys::ODESystem, args...; kwargs...)
    return periodic_orbit(Float64, sys, args...; kwargs...)
end


text:
using Nemo, LinearAlgebra, DelimitedFiles

# using Revise
using SymbolicCRN

M = [63 90 36 -45 -18
    -7 -10 -4 5 2
    28 40 16 -20 -8
    -56 -80 -32 40 16
    7 10 4 -5 -2]

K = [2 5 0 0 0 0
    0 0 1 0 0 1
    0 0 0 5 1 0
    0 7 2 4 0 0
    7 0 0 0 2 5]

R, vrs = PolynomialRing(Nemo.ZZ, vcat(["k$i" for i in 1:5], ["x$i" for i in 1:4]));

@test all(integermultiple(raysof(SymbolicCRN.cone_positiveorthant(5))) .== LinearAlgebra.I(5))
@test integermultiple(raysof(cone_positivenullspace([0 0; 0 0]))) == [1 0; 0 1]
@test raysof(cone_positivenullspace([1 -1; 2 -2])) == transpose([1 1])
Knew = raysof(cone_positivenullspace(M))
@test all(M * Knew .== 0)
@test all(map(>=(0), Knew))
@test K == integermultiple(Knew)
p = prod([vrs[i] - i for i in 1:7])
@test integermultiple(raysof(outernormalcone(p, 5))) == [
    0 1 0 0 0 0 0
    1 0 0 0 0 0 0
    0 0 1 0 0 0 0
    0 0 0 1 0 0 0
    0 0 0 0 0 0 -1
    0 0 0 0 1 0 0
    0 0 0 0 0 1 0
    0 0 0 0 0 0 0
    0 0 0 0 0 0 0
]

# TODO;
# verticesof(polytope) = permutedims(Rational.(polymake_dehomog(polytope.VERTICES)))


text:
using Test
using ExproniconLite

@testset "name_only" begin
    @test name_only(:(x::Int)) == :x
    @test name_only(:(T <: Int)) == :T
    @test name_only(:(Foo{T} where T)) == :Foo
    @test name_only(:(Foo{T})) == :Foo
    @test_throws ErrorException name_only(Expr(:fake))
end

@testset "rm_lineinfo" begin
    ex = quote
        1 + 1
        2 + 2
    end

    @test rm_lineinfo(ex) == Expr(:block, :(1 + 1), :(2 + 2))

    ex = quote
        Base.@kwdef mutable struct D
            field1::Union{ID, Missing, Nothing} = nothing
        end
        StructTypes.StructType(::Type{D}) = begin
            StructTypes.Mutable()
        end
        StructTypes.omitempties(::Type{D}) = begin
            true
        end
    end

    @test rm_lineinfo(ex).args[1].args[end] == rm_lineinfo(:(mutable struct D
        field1::Union{ID, Missing, Nothing} = nothing
    end))
    @test rm_lineinfo(ex).args[2] == rm_lineinfo(:(StructTypes.StructType(::Type{D}) = begin
        StructTypes.Mutable()
    end))
    @test rm_lineinfo(ex).args[3] == rm_lineinfo(:(StructTypes.omitempties(::Type{D}) = begin
        true
    end))
end

@testset "flatten_blocks" begin
    ex = quote
        1 + 1
        begin
            2 + 2
        end
    end

    @test rm_lineinfo(flatten_blocks(ex)) == Expr(:block, :(1+1), :(2+2))
end

@testset "rm_annotations" begin
    ex = quote
        x :: Int
        begin
            y::Float64
        end
    end

    @test rm_lineinfo(rm_annotations(ex)) == quote
        x
        begin
            y
        end
    end|>rm_lineinfo

    ex = :(sin(::Float64; x::Int=2))
    ex = rm_annotations(ex)
    @test ex.head === :call
    @test ex.args[1] === :sin
    @test ex.args[2].head === :parameters
    @test ex.args[2].args[1] === :x
    @test ex.args[3] isa Symbol
end

@testset "prettify" begin
    ex = quote
        x :: Int
        begin
            y::Float64
        end
    end

    @test prettify(ex) == quote
        x::Int
        y::Float64
    end|>rm_lineinfo
end

global_x = 2

@testset "eval_interp" begin
    ex = Expr(:call, :+, Expr(:$, :global_x), 1)
    @test eval_interp(Main, ex) == :(2 + 1)
end

@testset "eval_literal" begin
    ex = :(for i in 1:10
        1 + 1
    end)
    @test rm_lineinfo(eval_literal(Main, ex)) == rm_lineinfo(:(for i in $(1:10)
        2
    end))
end

@testset "subtitute" begin
    @test subtitute(:(x + 1), :x=>:y) == :(y + 1)
    @test subtitute(:(for i in 1:10;x += i;end), :x => :y) == :(for i in 1:10;y += i;end)
end