text:
export sample_points

import Distributions

"""
    sample_points(m::TriMesh, num_samples::Int=5000; returns_normals::Bool=false, eps::Number = 1e-6)

Uniformly samples `num_samples` points from the surface of TriMesh `m`.

`returns_normals` is optional keyword argument, to returns normals
from respective faces of samples. `eps` is optional keyword argument for
a small number to prevent division by zero for small surface areas.

### Examples:
```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> points = sample_points(m, 5000)
julia> points, normals = sample_points(m, 5000; returns_normals=true)
```
"""
function sample_points(
    m::TriMesh{T,R,S},
    num_samples::Int = 5000;
    eps::Number = EPS,
)::S{T,3} where {T,R,S}
    verts_padded = get_verts_padded(m)
    faces_padded = get_faces_padded(m)
    faces_areas_padded = compute_faces_areas_padded(m)

    #(Fi,B)
    #TODO: condition for probvec fails in Float32
    faces_areas_prob = Zygote.ignore() do
        faces_areas_padded = Float64.(faces_areas_padded)
        faces_areas_prob =
            faces_areas_padded ./ max.(sum(faces_areas_padded; dims = 2), eps)
        faces_areas_prob[:, end, :] +=
            reshape(max.(1 .- sum(faces_areas_prob; dims = 2), 0.0), 1, :)
        return faces_areas_prob
    end
    samples = @ignore similar(verts_padded, 3, num_samples, m.N)
    samples = Zygote.bufferfrom(samples)

    for (i, _len) in enumerate(m._faces_len)
        # Distributions package expects Array!
        probvec = cpu(faces_areas_prob[1, 1:_len, i])
        dist = Distributions.Categorical(probvec)
        sample_faces_idx = @ignore rand(dist, num_samples)
        sample_faces = faces_padded[:, sample_faces_idx, i]
        samples[:, :, i] = _sample_points(
            S,
            verts_padded[:, 1:m._verts_len[i], i],
            sample_faces,
            num_samples,
        )
    end

    return copy(samples)
end

function _sample_points(
    ::Type{S},
    verts::AbstractArray{T,2},
    sample_faces::AbstractArray{R,2},
    num_samples::Int,
) where {S,T,R}

    v1 = verts[:, sample_faces[1, :]]
    v2 = verts[:, sample_faces[2, :]]
    v3 = verts[:, sample_faces[3, :]]
    (w1, w2, w3) = S.(_rand_barycentric_coords(num_samples))
    samples = (w1 .* v1) + (w2 .* v2) + (w3 .* v3)
    return samples
end

function _rand_barycentric_coords(num_samples::Int)
    u = sqrt.(rand(Float32, 1, num_samples))
    v = rand(Float32, 1, num_samples)
    w1 = 1.0f0 .- u
    w2 = u .* (1.0f0 .- v)
    w3 = u .* v
    return (w1, w2, w3)
end

"""
    normalize!(m::TriMesh)

Normalize each mesh in TriMesh `m` with mean centered at origin and unit standard deviation
and overwrite the `m` with normalized TriMesh.

See also: [`normalize`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> normalize!(m)
```
"""
function normalize!(m::TriMesh{T,R,S}) where {T,R,S}
    verts_padded = get_verts_padded(m)
    _len = S(reshape(m._verts_len, 1, 1, :)) # move to `S` storage type
    _centroid = sum(verts_padded; dims = 2) ./ _len
    _correction = ((_centroid .^ 2) .* (m.V .- _len))
    _std =
        sqrt.(
            (sum((verts_padded .- _centroid) .^ 2; dims = 2) - _correction) ./ (_len .- 1),
        )
    verts_padded = (verts_padded .- _centroid) ./ max.(_std, EPS)
    m._verts_list = _padded_to_list(verts_padded, m._verts_len)
    return m
end

"""
    normalize(m::TriMesh)

Normalize each mesh in TriMesh `m` with mean centered at origin and unit standard deviation

See also: [`normalize!`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> m = normalize(m)
```
"""
function normalize(m::TriMesh)
    m = deepcopy(m)
    normalize!(m)
    return m
end

"""
    scale!(m::TriMesh, factor::Number)
    scale!(m::TriMesh, factor::AbstractArray{<:Number})

Scale the TriMesh `m` by scaling factor `factor`
and overwrite `m` with scaled TriMesh. If `factor`
is array of size `(3, )`, then TriMesh will be scale
by scaling factor of respective dimension.

Scaling factor `factor` (each element in case of array)
should be strictly greater than `0.0`.

See also: [`scale`](@ref)

### Examples:
```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> scale!(m, 1.0)
julia> scale!(m, [1.0, 1.0, 1.0])
```
"""
function scale!(m::TriMesh, factor::Float32)
    (factor > 0.0) || error("factor must be greater than 0.0")
    verts_packed = get_verts_packed(m)
    lmul!(factor, verts_packed)
    m._verts_packed = verts_packed
    return m
end

function scale!(m::TriMesh, factor::AbstractArray{Float32})
    (size(factor) == (3,)) ||
        error("factor must be (3, ), but instead got $(size(factor)) array")
    (factor .> 0.0) || error("factor must be greater than 0.0")
    verts_packed = get_verts_packed(m) .* reshape(factor, :, 1)
    m._verts_packed = verts_packed
    return m
end

scale!(m::TriMesh, factor) = scale!(m, Float32.(factor))

"""
    scale(m::TriMesh, factor::Number)
    scale(m::TriMesh, factor::AbstractArray{<:Number})

Scale the TriMesh `m` by scaling factor `factor`.
If `factor` is array of size `(3, )`, then TriMesh
will be scaleby scaling factor of respective dimension.

Scaling factor `factor` (each element in case of array)
should be strictly greater than `0.0`.

See also: [`scale!`](@ref)

### Examples:
```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> m = scale(m, 1.0)
julia> m = scale!(m, [1.0, 1.0, 1.0])
```
"""
function scale(m::TriMesh, factor::Union{Float32,AbstractArray{Float32}})
    m = deepcopy(m)
    scale!(m, factor)
    return m
end

scale(m::TriMesh, factor) = scale(m, Float32.(factor))

"""
    rotate!(m::TriMesh, rotmat::AbstractArray{<:Number,2})
    rotate!(m::TriMesh, rotmat::AbstractArray{<:Number,3})

Rotate the TriMesh `m` by rotation matrix `rotmat`
and overwrite `m` with rotated TriMesh.

Rotation matrix `rotmat` should be of size `(3,3)` or `(3,3,B)`
where B is the batch size of TriMesh.

See also: [`rotate`](@ref)

### Examples:
```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> rotmat = rand(3,3)
julia> rotate!(m, rotmat)
```
"""
function rotate!(m::TriMesh, rotmat::AbstractArray{Float32,2})
    size(rotmat) == (3, 3) ||
        error("rotmat must be (3, 3) array, but instead got $(size(rotmat)) array")
    verts_packed = transpose(rotmat) * get_verts_packed(m)
    m._verts_packed = verts_packed
    return m
end

function rotate!(m::TriMesh, rotmat::AbstractArray{Float32,3})
    size(rotmat) == (3, 3, m.N) ||
        error("rotmat must be (3, 3, $(m.N)) array, but instead got $(size(rotmat)) array")
    verts_padded = Flux.batched_mul(Flux.batched_transpose(rotmat), get_verts_padded(m))
    m._verts_padded = verts_padded
    return m
end

rotate!(m::TriMesh, rotmat::AbstractArray{<:Number}) = rotate!(m, Float32.(rotmat))

"""
    rotate(m::TriMesh, rotmat::AbstractArray{<:Number,2})

Rotate the TriMesh `m` by rotation matrix `rotmat`.

Rotation matrix `rotmat` should be of size `(3,3)`

See also: [`rotate!`](@ref)

### Examples:
```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> rotmat = rand(3,3)
julia> m = rotate(m, rotmat)
```
"""
function rotate(m::TriMesh, rotmat::AbstractArray{<:Number})
    m = deepcopy(m)
    rotate!(m, rotmat)
    return m
end

"""
    realign!(src::TriMesh, tgt::TriMesh)
    realign!(src::TriMesh, tgt_min::AbstractArray{<:Number,2}, tgt_max::AbstractArray{<:Number,2})

Re-Align the TriMesh `src` with the axis aligned bounding box of mesh at `index` in TriMesh `tgt`
and overwrite `src` with re-aligned TriMesh.

See also: [`realign`](@ref)

### Examples:
```jldoctest
julia> src = load_trimesh("teapot.obj")
julia> tgt = scale(src, 2.0)
julia> realign!(src, tgt)
```
"""
function realign!(
    src::TriMesh,
    tgt_min::AbstractArray{Float32,2},
    tgt_max::AbstractArray{Float32,2},
)
    verts_padded = get_verts_padded(src)
    src_min = minimum(verts_padded, dims = 2)
    src_max = maximum(verts_padded, dims = 2)
    verts_padded =
        ((verts_padded .- src_min) ./ (src_max - src_min .+ EPS)) .* (tgt_max - tgt_min) .+
        tgt_min
    src._verts_packed = _padded_to_packed(verts_padded, src._verts_len)
    return src
end

realign!(
    src::TriMesh,
    tgt_min::AbstractArray{<:Number,2},
    tgt_max::AbstractArray{<:Number,2},
) = realign!(src, Float32.(tgt_min), Float32.(tgt_max))

realign!(src::TriMesh, tgt::AbstractArray{<:Number,2}) = realign!(src, Float32.(tgt))

function realign!(src::TriMesh, tgt::AbstractArray{Float32,2})
    tgt_min = minimum(tgt, dims = 2)
    tgt_max = maximum(tgt, dims = 2)
    realign!(src, tgt_min, tgt_max)
    return src
end

realign!(src::TriMesh, tgt::TriMesh, index::Integer = 1) =
    realign!(src, Float32.(get_verts_list(tgt)[index]))

"""
    realign(src::TriMesh, tgt::TriMesh, index::Integer=1)
    realign(src::TriMesh, tgt_min::AbstractArray{<:Number,2}, tgt_max::AbstractArray{<:Number,2})

Re-Align the TriMesh `src` with the axis aligned bounding box of mesh at `index` in TriMesh `tgt`.

See also: [`realign`](@ref)

### Examples:
```jldoctest
julia> src = load_trimesh("teapot.obj")
julia> tgt = scale(src, 2.0)
julia> src = realign(src, tgt)
```
"""
function realign(
    src::TriMesh,
    tgt_min::AbstractArray{<:Number,2},
    tgt_max::AbstractArray{<:Number,2},
)
    src = deepcopy(src)
    realign!(src, tgt_min, tgt_max)
    return src
end

function realign(src::TriMesh, tgt::TriMesh, index::Integer = 1)
    src = deepcopy(src)
    realign!(src, tgt, index)
    return src
end

"""
    translate!(m::TriMesh, vector::Number)
    translate!(m::TriMesh, vector::AbstractArray{<:Number})

Translate the TriMesh `m` by translating vector `vector`
and overwrite `m` with translated TriMesh. If `vector`
is a number, then TriMesh will be translated by same number
in all dimension.

See also: [`translate`](@ref)

### Examples:
```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> translate!(m, 0.0)
julia> translate!(m, [0.0, 0.0, 0.0])
```
"""
translate!(m::TriMesh, vector::Float32) = translate!(m, fill(vector, (3,)))

function translate!(m::TriMesh, vector::AbstractArray{Float32})
    (size(vector) == (3,)) ||
        error("vector must be (3, ), but instead got $(size(vector)) array")
    verts_packed = get_verts_packed(m)
    verts_packed = verts_packed .+ reshape(vector, :, 1)
    m._verts_packed = verts_packed
    return m
end

translate!(m::TriMesh, vector) = translate!(m, Float32.(vector))

"""
    translate(m::TriMesh, vector::Number)
    translate(m::TriMesh, vector::AbstractArray{<:Number})

Translate the TriMesh `m` by translating vector `vector`.
If `vector` is a number, then TriMesh will be translated
by same number in all dimension.

See also: [`translate!`](@ref)

### Examples:
```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> m = translate(m, 0.0)
julia> m = translate(m, [0.0, 0.0, 0.0])
```
"""
function translate(m::TriMesh, args...)
    m = deepcopy(m)
    translate!(m, args...)
end

"""
    offset!(m::TriMesh, offset_verts::AbstractArray{<:Number,2})

Add offset to the vertices of the TriMesh `m` by offset vertices `offset_verts_packed`
and overwrite `m` with updated TriMesh.

See also: [`offset`](@ref)

### Examples:
```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> offset_verts = ones(get_verts_packed(m))
julia> offset!(m, offset_verts)
```
"""

function offset!(m::TriMesh, offset_verts_packed::AbstractArray{Float32,2})
    verts_packed = get_verts_packed(m)
    (size(offset_verts_packed) == size(verts_packed)) ||
        error("mesh and offset_verts size mismatch")
    verts_packed += offset_verts_packed
    m._verts_packed = verts_packed
    return m
end

offset!(m::TriMesh, offset_verts_packed::AbstractArray{<:Number,2}) =
    offset!(m, Float32.(offset_verts_packed))

"""
    offset(m::TriMesh, offset_verts_packed::AbstractArray{<:Number,2})

Add offset to the vertices of the TriMesh `m` by offset vertices `offset_verts_packed`.

See also: [`offset!`](@ref)

### Examples:
```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> offset_verts = ones(get_verts_packed(m))
julia> m = offset(m, offset_verts)
```
"""
function offset(m::TriMesh, offset_verts_packed::AbstractArray{<:Number,2})
    m = deepcopy(m)
    offset!(m, offset_verts_packed)
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() # HIDE
let
    odeSol(x,y) = Point(-x, 2y) # x'(t) = -x, y'(t) = 2y
    fig = Figure(resolution = (700, 450), fontsize = 18, font = "sans")
    ax = Axis(fig, xlabel = "x", ylabel = "y", backgroundcolor = :black)
    stplt = streamplot!(ax, odeSol, -2..4, -2..2, colormap = Reverse(:plasma),
        gridsize= (32,32), arrow_size = 10)
    fig[1, 1] = ax
    fig
    save(joinpath(@__DIR__, "output", "streamplot.png"), fig, px_per_unit = 2) # HIDE
end
using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

ext: jl

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using CSV, JSON, DataFrames, JLD, MAT, FileIO, JLD2
include("aggreg.jl")

"""Equivalent annual cost (EAC)
# Arguments
* `cost::Real`: Net present value of a project.
* `n::Integer`: Number of payments.
* `r::Real`: Interest rate.
"""
function equivalent_annual_cost(cost::Real, n::Integer, r::Real)
    factor = if iszero(r) n else (1 - (1 + r)^(-n)) / r end
    return cost / factor
end

"""Loads parameter values for an instance from CSV and JSON files. Reads the following files from `instances_path`.
- `indices.json` with fields `G`, `G_r`, `N`, `L`, `T`, `S`
- `constants.json` with fields `kappa`, `C`, `C_bar`, `r`
- `nodes/` -- Time clustered data from the nodes with fields `Dem_Inc`, `Load_mod`, `Max_Load`, `Avail_Win`, `Avail_Sol`
  - `1.csv`
  - `2.csv`
  - ...
- `technology.csv` with fields `cost`, `lifetime`, `M`, `fuel_cost_1`, `fuel_cost_2`, `r_minus`, `r_plus`
- `transmission.csv` with fields `M`, `cost`, `dist`, `lifetime`, `C`, `B`
- `storage.csv` with fields `xi`, `cost`, `lifetime`, `C`, `b0_1, ..., b0_n`
# Arguments
- `DataInput_path::AbstractString`: Path to the instance directory.
"""
function Params(DataInput_path::AbstractString, instances_path_ftr::AbstractString)
    # Load indexes and constant parameters
    indices = JSON.parsefile(joinpath(instances_path_ftr, "indices.json"))

    # Load indices. Convert JSON values to right types.
    G = indices["G"] |> Array{Int}
    G_r = indices["G_r"] |> Array{Int}
    N = indices["N"] |> Array{Int}
    L = indices["L"] |> Array{Array{Int}}
    L_ind = indices["L_ind"] |> Array{Int}
    S = indices["S"] |> Array{Int}
    H = indices["H"] |> Array{Int}
    T = 1:indices["T"] |> Array{Int}

    # Load constant parameters
    constants = JSON.parsefile(joinpath(DataInput_path, "constants.json"))
    κ = constants["kappa"]
    μ = constants["mu"]
    C = constants["C"]
    C̄ = constants["C_bar"]
    C_E = constants["C_E"]
    interest_rate = constants["r"]
    R_E = constants["R_E"]

    # Load time clustered parameters and region
    D_nt = zeros(length(N), length(T))
    A_gnt = ones(length(G), length(N), length(T))
    AH_nt = zeros(length(N), length(T))
    AR_nt = zeros(length(N), length(T))

    # Load generation parameters (per node / per generation technology)
    Gmin_gn = zeros(length(G), length(N))
    Gmax_gn = zeros(length(G), length(N))
    gen_capacity = CSV.File(joinpath(instances_path_ftr, "gen_capacity.csv")) |> DataFrame

    for n in N
        # Load node values from CSV files.
        nodes = CSV.File(joinpath(instances_path_ftr, "nodes", "$n.csv")) |> DataFrame
        for g in G
            line_search = findfirst((gen_capacity.gen_tech .== g) .& (gen_capacity.node .== n))
            Gmin_gn[g,n] = gen_capacity.gcap_min[line_search]
            Gmax_gn[g,n] = gen_capacity.gcap_max[line_search]
        end
        D_nt[n, :] = round.(nodes.Demand[T]; digits = 5)
        A_gnt[1, n, :] = round.(nodes.Avail_Wind_On[T]; digits = 5)
        A_gnt[2, n, :] = round.(nodes.Avail_Wind_Off[T]; digits = 5)
        A_gnt[3, n, :] = round.(nodes.Avail_Sol[T]; digits = 5)
        A_gnt[A_gnt .< 0.001] .= 0
        AH_nt[n,:] = round.(nodes.Hyd_In[T]; digits = 0)
        AR_nt[n,:] = round.(nodes.HydRoR_In[T]; digits = 0)
    end

    # Nodes specifications
    region_n = Array{AbstractString, 1}(undef, length(N))
    max_dem_n = Array{Float64, 1}(undef, length(N))

    nodes_specs = CSV.File(joinpath(instances_path_ftr, "nodes_specs.csv")) |> DataFrame

    for n in N
        region_n[n] =  nodes_specs.Name[n]
        max_dem_n[n] = nodes_specs.Max_Demand[n]
    end

    # Load weights
    clust_weights = CSV.File(joinpath(instances_path_ftr, "weights.csv")) |> DataFrame
    τ_t = clust_weights.Weights[T]

    # Load technology parameters
    gen_technology = joinpath(instances_path_ftr, "gen_technology.csv") |>
        CSV.File |> DataFrame
    I_g = equivalent_annual_cost.(gen_technology.investment_cost .* 1000, gen_technology.lifetime,
                                  interest_rate) |> Array{AbstractFloat, 1}
    M_g = gen_technology.fixedOM .* 1000 |> Array{AbstractFloat, 1}
    C_g = gen_technology.fuel_cost ./ gen_technology.efficiency .+ gen_technology.varOM |> Array{AbstractFloat, 1}
    e_g = gen_technology.efficiency |> Array{AbstractFloat, 1}
    E_g = gen_technology.emissions |> Array{AbstractFloat, 1}
    r⁻_g = gen_technology.r_minus |> Array{AbstractFloat, 1}
    r⁺_g = gen_technology.r_plus |> Array{AbstractFloat, 1}
    technology_g = gen_technology.name

    # Load transmission parameters
    I_l = zeros(length(L_ind))
    M_l = zeros(length(L_ind))
    C_l = zeros(length(L_ind))
    B_l = zeros(length(L_ind))
    e_l = zeros(length(L_ind))
    Tmin_l = zeros(length(L_ind))
    Tmax_l = zeros(length(L_ind))
    transmission = joinpath(instances_path_ftr, "transmission.csv") |>
        CSV.File |> DataFrame
    I_l = equivalent_annual_cost.(transmission.cost[1] .* transmission.dist .+ transmission.converter_cost[1],
                                  transmission.lifetime[1]|>Int64, interest_rate) |> Array{AbstractFloat, 1}
    M_l = transmission.M .* I_l |> Array{AbstractFloat, 1}
    C_l = transmission.C |> Array{AbstractFloat, 1}
    B_l = transmission.B |> Array{AbstractFloat, 1}
    e_l = transmission.efficiency |> Array{AbstractFloat, 1}
    Tmin_l = transmission.tcap_min |> Array{AbstractFloat, 1}
    Tmax_l = transmission.tcap_max |> Array{AbstractFloat, 1}

    # Load storage parameters
    ξ_s = zeros(length(S))
    I_s = zeros(length(S))
    C_s = zeros(length(S))
    Smin_sn = zeros(length(S),length(N))
    Smax_sn = zeros(length(S),length(N))
    storage = joinpath(instances_path_ftr, "storage.csv") |>
        CSV.File |> DataFrame
    sto_capacity = joinpath(instances_path_ftr, "sto_capacity.csv") |>
        CSV.File |> DataFrame
    for s in S
        ξ_s = storage.xi |> Array{AbstractFloat, 1}
        I_s = equivalent_annual_cost.(storage.cost .* 1000, storage.lifetime, interest_rate) |> Array{AbstractFloat, 1}
        C_s = storage.C |> Array{AbstractFloat, 1}
        for n in N
            line_search = findfirst((sto_capacity.s .== s) .& (sto_capacity.node .== n))
            Smin_sn[s,n] = sto_capacity.scap_min[line_search]
            Smax_sn[s,n] = sto_capacity.scap_max[line_search]
        end
    end
    Smin_sn = Smin_sn  |> Array{AbstractFloat, 2}
    Smax_sn = Smax_sn  |> Array{AbstractFloat, 2}

    # Load hydro capacity and technology parameters
    Wmax_hn = zeros(length(H),length(N))
    Wmin_hn = zeros(length(H),length(N))
    Hmin_hn = zeros(length(H),length(N))
    Hmax_hn = zeros(length(H),length(N))
    Fmin_n = zeros(length(N))
    HRmax_n = zeros(length(N))

    hydro_capacity = joinpath(instances_path_ftr, "hydro_capacity.csv") |> CSV.File |> DataFrame
    for h in H, n in N
        line_search = findfirst((hydro_capacity.hydro_tech .== h) .& (hydro_capacity.node .== n))
        Hmin_hn[h,n] = hydro_capacity.hcap_min[line_search]
        Hmax_hn[h,n] = hydro_capacity.hcap_max[line_search]
        Wmin_hn[h,n] = hydro_capacity.wcap_min[line_search]
        Wmax_hn[h,n] = hydro_capacity.wcap_max[line_search]
    end

    hydro = joinpath(instances_path_ftr, "hydro.csv") |> CSV.File |> DataFrame
    HRmax_n[1:length(N)] = hydro.HydroRoR[1:length(N)] |> Array{AbstractFloat, 1}
    Fmin_n[1:length(N)] = hydro.hyd_flow_min[1:length(N)] |> Array{AbstractFloat, 1}

    hydro_technology = joinpath(instances_path_ftr, "hydro_technology.csv") |> CSV.File |> DataFrame;
    I_h = equivalent_annual_cost.(hydro_technology.investment_cost .* 1000, hydro_technology.lifetime,
                                    interest_rate) |> Vector{Float64}
    M_h = hydro_technology.fixedOM .* 1000 |> Vector{Float64}
    C_h = hydro_technology.fuel_cost ./ hydro_technology.efficiency .+ hydro_technology.varOM |> Vector{Float64}
    e_h = hydro_technology.efficiency |> Vector{Float64}
    E_h = hydro_technology.emissions |> Vector{Float64}
    r⁻_h = hydro_technology.r_minus |> Vector{Float64}
    r⁺_h = hydro_technology.r_plus |> Vector{Float64}

    # Return Params struct
    Params(
        region_n, max_dem_n, technology_g, G, G_r, N, L, L_ind, T, S, H, κ, μ, C, C̄, C_E, R_E, τ_t, Gmin_gn, Gmax_gn, A_gnt, D_nt, I_g, M_g, C_g,
        e_g, E_g, r⁻_g, r⁺_g, I_l, M_l, C_l, B_l, e_l, Tmin_l, Tmax_l, ξ_s, I_s, C_s, Smin_sn, Smax_sn,
        Wmax_hn, Wmin_hn, Hmax_hn, Hmin_hn, HRmax_n, Fmin_n, AH_nt, AR_nt,
        I_h, M_h, C_h, e_h, E_h, r⁻_h, r⁺_h)
end

"""Perform function Params() and save objects (i.e., A_gnt, AH_nt, etc.) as a dictionary.
# Arguments
- `constants_path::AbstractString`: Path to the constants JSON file
- `instances_path_ftr::AbstractString`: Path to the instance directory
- `output_path::AbstractString`: Path where the dictionary with the parameters is saved.
"""
function perform_Params(constants_path::AbstractString, instances_path_ftr::AbstractString, params_path::AbstractString)
    # Create parameters object
    parameters = Params(constants_path, instances_path_ftr)

    # Create parameters dictionary folder
    mkpath(params_path)

    # Saving parameters into a dictionary
    parameters_dict = Dict(fieldnames(typeof(parameters)) .|> String .=> [getfield(parameters,i) for i in fieldnames(typeof(parameters))])
    JLD2.save(joinpath(params_path, "parameters.jld2"), parameters_dict)
end

"""Read time-dependent parameters (i.e., T, D_nt, A_gnt, AH_nt, AR_nt) and create a new parameters struct.
# Arguments
- `params_path_ftr::AbstractString`: Path to the FTR parameters
- `instance_path_clust::AbstractString`: Path to the instance folder.
"""
function change_time_parameters(params_path_ftr::AbstractString, instance_path_clust::AbstractString; nosun::Bool = false)
    # Reading FTR parameters (assuming perform_Params was ran and the dictionary parameters.jld2 is available)
    ParamsDict = load(joinpath(params_path_ftr,"parameters.jld2"))

    # Declare representative periods
    rep_periods = JSON.parsefile(joinpath(instance_path_clust,"rep_periods.json"))
    T = 1:rep_periods["T"]
    ParamsDict["T"] = T

    # Load weights
    clust_weights = CSV.File(joinpath(instance_path_clust, "weights.csv")) |> DataFrame
    τ_t = clust_weights.Weights[T]
    ParamsDict["τ_t"] = τ_t

    # Read time-dependent parameters
    D_nt = zeros(length(ParamsDict["N"]), length(T))
    A_gnt = ones(length(ParamsDict["G"]), length(ParamsDict["N"]), length(T))
    AH_nt = zeros(length(ParamsDict["N"]), length(T))
    AR_nt = zeros(length(ParamsDict["N"]), length(T))

    for n in ParamsDict["N"]
            # Load node values from CSV files.
            nodes = CSV.File(joinpath(instance_path_clust,"nodes","$n.csv")) |> DataFrame
            D_nt[n, :] = round.(nodes.Demand[T]; digits = 5)
            A_gnt[1, n, :] = round.(nodes.Avail_Wind_On[T]; digits = 5)
            A_gnt[2, n, :] = round.(nodes.Avail_Wind_Off[T]; digits = 5)
            A_gnt[3, n, :] = round.(nodes.Avail_Sol[T]; digits = 5)
            A_gnt[A_gnt .< 0.001] .= 0
            AH_nt[n,:] = round.(nodes.Hyd_In[T]; digits = 0)
            AR_nt[n,:] = round.(nodes.HydRoR_In[T]; digits = 0)
    end

    ParamsDict["A_gnt"] = A_gnt
    ParamsDict["D_nt"] = D_nt
    ParamsDict["AH_nt"] = AH_nt
    ParamsDict["AR_nt"] = AR_nt

    if nosun
        gsun = 3
        ParamsDict["A_gnt"][gsun,:,:] .= 0
    end

    # Forming parameters struct
    parameters = Params(
        ParamsDict["region_n"], ParamsDict["max_dem_n"], ParamsDict["technology_g"], ParamsDict["G"], ParamsDict["G_r"], ParamsDict["N"], ParamsDict["L"], ParamsDict["L_ind"],
        ParamsDict["T"], ParamsDict["S"], ParamsDict["H"], ParamsDict["κ"], ParamsDict["μ"], ParamsDict["C"], ParamsDict["C̄"], ParamsDict["C_E"], ParamsDict["R_E"], ParamsDict["τ_t"],
        ParamsDict["Gmin_gn"], ParamsDict["Gmax_gn"], ParamsDict["A_gnt"], ParamsDict["D_nt"], ParamsDict["I_g"], ParamsDict["M_g"], ParamsDict["C_g"], ParamsDict["e_g"], ParamsDict["E_g"],
        ParamsDict["r⁻_g"], ParamsDict["r⁺_g"], ParamsDict["I_l"], ParamsDict["M_l"], ParamsDict["C_l"], ParamsDict["B_l"], ParamsDict["e_l"], ParamsDict["Tmin_l"], ParamsDict["Tmax_l"],
        ParamsDict["ξ_s"], ParamsDict["I_s"], ParamsDict["C_s"], ParamsDict["Smin_sn"], ParamsDict["Smax_sn"], ParamsDict["Wmax_hn"], ParamsDict["Wmin_hn"], ParamsDict["Hmax_hn"],
        ParamsDict["Hmin_hn"], ParamsDict["HRmax_n"], ParamsDict["Fmin_n"], ParamsDict["AH_nt"], ParamsDict["AR_nt"], ParamsDict["I_h"], ParamsDict["M_h"], ParamsDict["C_h"],
        ParamsDict["e_h"], ParamsDict["E_h"], ParamsDict["r⁻_h"], ParamsDict["r⁺_h"]
    );
    return parameters
end

"""Save object into JSON file.
# Arguments
- `object`
- `output_path::AbstractString`: Full filepath, e.g., `path.json`.
"""
function save_json(object, filepath::AbstractString)
    open(filepath, "w") do io
        JSON.print(io, object)
    end
end

flatten(x::Array{<:Array, 1}) = Iterators.flatten(x)|> collect|> flatten
flatten(x::Array{<:Number, 1}) = x
flatten(x::AffExpr) = x
flatten(x::VecOrMat{AffExpr}) = x[:]
flatten(x::Array{<:Any}) = x[:]

shape(x::Array{<:Array, 1}) = vcat(shape(first(x)), [length(x)])
shape(x::Array{<:Number, 1}) = [length(x)]

mdim(x::Array{<:Array, 1}) = reshape(flatten(x), shape(x)...)
mdim(x::Array{<:Number, 1}) = x

transform(x::Array, t::Type{Array{T, N}}) where T <: Number where N = mdim(x)
transform(x::Array, t::Type{Array{T, 1}}) where T <: Array = x
transform(x::Number, t::Type{T}) where T <: Number = x
transform(x::Array, t::Type{Array{AbstractString, 1}}) = x

function convert_type(::Type{Array{T, N}}) where T <: Number where N
    t = T
    for i in 1:N
        t = Array{t, 1}
    end
    return t
end
convert_type(::Type{Array{T, 1}}) where T <: Array = Array{T, 1}
convert_type(t::Type{T}) where T <: Number = t
convert_type(::Type{Array{AbstractString, 1}}) = Array{AbstractString, 1}

"""Load values to type from JSON file.
# Arguments
- `type`
- `filepath::AbstractString`
"""
function load_json(type, filepath::AbstractString)
    objects = JSON.parsefile(filepath)
    fields = []
    for (s, t) in zip(fieldnames(type), fieldtypes(type))
        push!(fields,
              objects[string(s)] |>
              v -> convert(convert_type(t), v) |>
              v -> transform(v, t))
    end
    type(fields...)
end

#Replace NaN's with 0
function replace_nans!(array::Array{Float64, N}) where N
    for i in eachindex(array)
        if isnan(array[i])
            array[i] = zero(1)
        end
    end
end

"""
read_clusters(ClustersDataPath::AbstractString, k::Int)
Read clusters from a julia dictionary (.jld2) with the structure ClustInstance.
# Output:
- _ClustInstance: Instance of clustering with weights, centroids, etc.
- _SeriesInstance: Instance of series used to form the clusters.
"""
function read_clusters(ClustersDataPath::AbstractString, k::Int)
    _ClustUpdate = load(joinpath(ClustersDataPath,"clust_out.jld2"))
    _SeriesUpdate = load(joinpath(ClustersDataPath,"series_out.jld2"));

    return _ClustUpdate[string(k)], _SeriesUpdate[string(k)]
end

"""
unpack_clusters_features(_ClustInstance::ClustInstance)
# Output:
- k_cent: centroids.
- weights: respective centroids' weights.
- series_clust: mapping between series and clusters.
"""
function unpack_clusters_features(_ClustInstance::ClustInstance)
    @assert _ClustInstance.nclusters == _ClustInstance.series_clust[end] "Number of clusters inconsistent in the cluster instance"

    k_cent = _ClustInstance.k_cent
    weights = _ClustInstance.weights
    series_clust = _ClustInstance.series_clust

    return k_cent, weights, series_clust
end


text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.MPIStateArrays
using CLIMA.DGBalanceLawDiscretizations
using Printf
using LinearAlgebra
using Logging

const _nauxstate = 7
@inline function auxiliary_state_initialization!(aux, x, y, z, dim)
    @inbounds begin
        aux[1] = x
        aux[2] = y
        aux[3] = z
        if dim == 2
            aux[4] = x * y + z * y
            aux[5] = 2 * x * y + sin(x) * y^2 / 2 - (z - 1)^2 * y^3 / 3
        else
            aux[4] = x * z + z^2 / 2
            aux[5] = 2 * x * z + sin(x) * y * z - (1 + (z - 1)^3) * y^2 / 3
        end
    end
end
@inline function integral_knl(val, Q, aux)
    @inbounds begin
        x, y, z = aux[1], aux[2], aux[3]
        val[1] = x + z
        val[2] = 2 * x + sin(x) * y - (z - 1)^2 * y^2
    end
end

using Test
function run(mpicomm, dim, ArrayType, Ne, N, FT)

    brickrange = ntuple(j -> range(FT(0); length = Ne[j] + 1, stop = 3), dim)
    topl = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = ntuple(j -> true, dim),
    )

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    spacedisc = DGBalanceLaw(
        grid = grid,
        length_state_vector = 0,
        flux! = (x...) -> (),
        numerical_flux! = (x...) -> (),
        auxiliary_state_length = _nauxstate,
        auxiliary_state_initialization! = (x...) ->
            auxiliary_state_initialization!(x..., dim),
    )

    Q = MPIStateArray(spacedisc)

    DGBalanceLawDiscretizations.indefinite_stack_integral!(
        spacedisc,
        integral_knl,
        Q,
        (6, 7),
    )

    # Wrapping in Array ensure both GPU and CPU code use same approx
    @test Array(spacedisc.auxstate.data[:, 4, :]) ≈
          Array(spacedisc.auxstate.data[:, 6, :])
    @test Array(spacedisc.auxstate.data[:, 5, :]) ≈
          Array(spacedisc.auxstate.data[:, 7, :])
end

let
    CLIMA.init()
    ArrayTypes = (CLIMA.array_type(),)

    mpicomm = MPI.COMM_WORLD

    numelem = (5, 5, 5)
    lvls = 1

    polynomialorder = 4

    @testset "$(@__FILE__)" for ArrayType in ArrayTypes
        for FT in (Float64,) #Float32)
            for dim in 2:3
                err = zeros(FT, lvls)
                for l in 1:lvls
                    @info (ArrayType, FT, dim)
                    run(
                        mpicomm,
                        dim,
                        ArrayType,
                        ntuple(j -> 2^(l - 1) * numelem[j], dim),
                        polynomialorder,
                        FT,
                    )
                end
            end
        end
    end
end

nothing

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
module TimeFrames

import Base: range, +, -, *

using Dates

export TimeFrame, Boundary
export YearBegin, YearEnd
export MonthBegin, MonthEnd
# export Millisecond, Second, Minute, Hour, Day, Week
export NoTimeFrame
export apply, range
export Begin, End
export @tf_str

abstract type TimeFrame end

@enum(Boundary,
    UndefBoundary = 0,  # Undefined boundary
    Begin = 1,  # begin of interval
    End   = 2,  # end of interval
)

#T should be Dates.TimePeriod

abstract type AbstractPeriodFrame <: TimeFrame end
abstract type AbstractTimePeriodFrame <: AbstractPeriodFrame end
abstract type AbstractDatePeriodFrame <: AbstractPeriodFrame end

struct TimePeriodFrame{T <: Dates.TimePeriod} <: AbstractTimePeriodFrame
    period::T
    boundary::Boundary
end
TimePeriodFrame{T}(; boundary=Begin::Boundary) where T<:Dates.TimePeriod = TimePeriodFrame(T(1), boundary)
TimePeriodFrame{T}(n::Integer; boundary=Begin::Boundary) where T<:Dates.TimePeriod = TimePeriodFrame(T(n), boundary)

struct DatePeriodFrame{T <: Dates.DatePeriod} <: AbstractDatePeriodFrame
    period::T
    boundary::Boundary
end
DatePeriodFrame{T}(; boundary=Begin::Boundary) where T<:Dates.DatePeriod = DatePeriodFrame(T(1), boundary)
DatePeriodFrame{T}(n::Integer; boundary=Begin::Boundary) where T<:Dates.DatePeriod = DatePeriodFrame(T(n), boundary)

struct NoTimeFrame <: TimeFrame
    NoTimeFrame(args...; kwargs...) = new()
end
TimeFrame() = NoTimeFrame()

# Base.hash(tf::TimePeriodFrame, h::UInt) = hash(tf.period, hash(tf.boundary))
# Base.:(==)(tf1::TimePeriodFrame, tf2::TimePeriodFrame) = hash(tf1) == hash(tf2)

function _period_step(::Type{Date})
    Dates.Day(1)
end

const period_step = Dates.Millisecond(1)

function _period_step(::Type{DateTime})
    period_step
end

#struct Microsecond <: AbstractTimePeriodFrame
#    period::Dates.TimePeriod
#    boundary::Boundary
#end
#Microsecond() = Microsecond(Dates.Microsecond(1), Begin)
#Microsecond(n::Integer) = Microsecond(Dates.Microsecond(n), Begin)

struct Millisecond <: AbstractTimePeriodFrame
    period::Dates.Millisecond
    boundary::Boundary
end
Millisecond() = Millisecond(Dates.Millisecond(1), Begin)
Millisecond(n::Integer) = Millisecond(Dates.Millisecond(n), Begin)

struct Second <: AbstractTimePeriodFrame
    period::Dates.Second
    boundary::Boundary
end
Second() = Second(Dates.Second(1), Begin)
Second(n::Integer) = Second(Dates.Second(n), Begin)

struct Minute <: AbstractTimePeriodFrame
    period::Dates.Minute
    boundary::Boundary
end
Minute() = Minute(Dates.Minute(1), Begin)
Minute(n::Integer) = Minute(Dates.Minute(n), Begin)

struct Hour <: AbstractTimePeriodFrame
    period::Dates.Hour
    boundary::Boundary
end
Hour() = Hour(Dates.Hour(1), Begin)
Hour(n::Integer) = Hour(Dates.Hour(n), Begin)

struct Day <: AbstractDatePeriodFrame
    period::Dates.Day
    boundary::Boundary
end
Day() = Day(Dates.Day(1), Begin)
Day(n::Integer) = Day(Dates.Day(n), Begin)

struct Week <: AbstractDatePeriodFrame
    period::Dates.Week
    boundary::Boundary
end
Week() = Week(Dates.Week(1), Begin)
Week(n::Integer) = Week(Dates.Week(n), Begin)

struct MonthEnd <: AbstractDatePeriodFrame
    period::Dates.Month
    boundary::Boundary
end
MonthEnd() = MonthEnd(Dates.Month(1), End)
MonthEnd(n::Integer) = MonthEnd(Dates.Month(n), End)

struct MonthBegin <: AbstractDatePeriodFrame
    period::Dates.Month
    boundary::Boundary
end
MonthBegin() = MonthBegin(Dates.Month(1), Begin)
MonthBegin(n::Integer) = MonthBegin(Dates.Month(n), Begin)

struct YearEnd <: AbstractDatePeriodFrame
    period::Dates.Year
    boundary::Boundary
end
YearEnd() = YearEnd(Dates.Year(1), End)
YearEnd(n::Integer) = YearEnd(Dates.Year(n), End)

struct YearBegin <: AbstractDatePeriodFrame
    period::Dates.Year
    boundary::Boundary
end
YearBegin() = YearBegin(Dates.Year(1), Begin)
YearBegin(n::Integer) = YearBegin(Dates.Year(n), Begin)

const _D_STR2TIMEFRAME = Dict(
    "A"=>YearEnd,
    "AS"=>YearBegin,
    "M"=>MonthEnd,
    "MS"=>MonthBegin,
    "W"=>Week,
    "D"=>Day,
    "H"=>Hour,
    "T"=>Minute,
    "S"=>Second,
    "L"=>Millisecond,
    #"U"=>Microsecond,
    ""=>NoTimeFrame
)
# Reverse key/value
const _D_TIMEFRAME2STR = Dict{DataType,String}()
for (key, typ) in _D_STR2TIMEFRAME
    _D_TIMEFRAME2STR[typ] = key
end
# Additional shortcuts
const _D_STR2TIMEFRAME_ADDITIONAL = Dict(
    "MIN"=>Minute,
)
for (key, value) in _D_STR2TIMEFRAME_ADDITIONAL
    _D_STR2TIMEFRAME[key] = value
end

# To string
function String(tf::TimeFrame)
    s_tf = _D_TIMEFRAME2STR[typeof(tf)]
    if tf.period.value == 1
        s_tf
    else
        "$(tf.period.value)$(s_tf)"
    end
end

# Parse
function TimeFrame(s::String; boundary=UndefBoundary)
    freq_pattern = join(keys(_D_STR2TIMEFRAME), "|")
    #pattern = r"^([\d]*)([Y|M|W|D|H|T])$"
    pattern = Regex("^([\\d]*)([$freq_pattern]*)\$", "i")
    m = match(pattern, s)
    if m == nothing
        error("Can't parse '$s' to TimeFrame")
    else
        s_freq = uppercase(m[2])
        tf_typ = _D_STR2TIMEFRAME[s_freq]
        if m[1] != ""
            value = parse(Int, m[1])
        else
            value = 1
        end
        tf = tf_typ(value)
        if boundary != UndefBoundary
            tf.boundary = boundary
        end
        tf
    end
end

# grouper
function dt_grouper(tf::AbstractPeriodFrame)
    dt -> _d_f_boundary[tf.boundary](dt, tf.period)
end

function dt_grouper(tf::AbstractPeriodFrame, t::Type)
    if tf.boundary == Begin
        dt -> _d_f_boundary[tf.boundary](dt, tf.period)
    elseif tf.boundary == End
        dt -> _d_f_boundary[tf.boundary](dt, tf.period) - _period_step(t)
    else
        error("Unsupported boundary $(tf.boundary)")
    end
end

struct CustomTimeFrame <: TimeFrame
    f_group::Function
end

function TimeFrame(f_group::Function)
    CustomTimeFrame(f_group)
end

function TimeFrame(td::Dates.TimePeriod; boundary=Begin::Boundary)
    T = typeof(td)
    TimePeriodFrame{T}(td.value, boundary=boundary)
end

function TimeFrame(td::Dates.DatePeriod; boundary=Begin::Boundary)
    T = typeof(td)
    DatePeriodFrame{T}(td.value, boundary=boundary)
end

function dt_grouper(tf::CustomTimeFrame, ::Type)
    tf.f_group
end

const _d_f_boundary = Dict(
    Begin::Boundary => floor,
    End::Boundary => ceil
)

function apply(tf::TimeFrame, dt)
    dt_grouper(tf, typeof(dt))(dt)
end

function tonext(tf::TimeFrame, dt::Dates.TimeType; same=false)
    dt2 = apply(tf, dt)
    if dt2 < dt
        dt2 + tf
    else
        if !same && dt2 == dt
            dt2 + tf
        else
            dt2
        end
    end
end

# range
function range(dt1::Dates.TimeType, tf::AbstractPeriodFrame, dt2::Dates.TimeType; apply_tf=true)
    td = _period_step(typeof(dt2))
    if apply_tf
        apply(tf, dt1):tf.period:apply(tf, dt2-td)
    else
        dt1:tf.period:dt2-td
    end
end

function range(dt1::Dates.TimeType, tf::AbstractPeriodFrame, len::Integer)
    range(dt1, step=tf.period, length=len)
end

function range(tf::AbstractPeriodFrame, dt2::Dates.TimeType, len::Integer)
    range(dt2 - len * tf.period, step=tf.period, length=len)
end

range(dt1::DateTime, tf::NoTimeFrame, dt2::DateTime) = [dt1]

macro tf_str(tf)
  :(TimeFrame($tf))
end

promote_timetype(::Type{DateTime}, ::Type) = DateTime

promote_timetype(::Type{Date}, ::Type)              = Date
promote_timetype(::Type{Date}, ::Type{Hour})        = DateTime
promote_timetype(::Type{Date}, ::Type{Minute})      = DateTime
promote_timetype(::Type{Date}, ::Type{Second})      = DateTime
promote_timetype(::Type{Date}, ::Type{Millisecond}) = DateTime

promote_timetype(::Type{Dates.Time}, ::Type)             = Dates.Time
promote_timetype(::Type{Dates.Time}, ::Type{YearBegin})  = throw(InexactError(:none, Any, nothing))
promote_timetype(::Type{Dates.Time}, ::Type{YearEnd})    = throw(InexactError(:none, Any, nothing))
promote_timetype(::Type{Dates.Time}, ::Type{MonthBegin}) = throw(InexactError(:none, Any, nothing))
promote_timetype(::Type{Dates.Time}, ::Type{MonthEnd})   = throw(InexactError(:none, Any, nothing))
promote_timetype(::Type{Dates.Time}, ::Type{Week})       = throw(InexactError(:none, Any, nothing))
promote_timetype(::Type{Dates.Time}, ::Type{Day})        = throw(InexactError(:none, Any, nothing))

+(t::T, tf::TF) where {T<:Dates.TimeType, TF<:TimeFrame} =
  convert(promote_timetype(T, TF), t) + tf.period

+(tf::TimeFrame, t::TimeType) = t + tf

-(t::T, tf::TF) where {T<:Dates.TimeType, TF<:TimeFrame} =
  convert(promote_timetype(T, TF), t) - tf.period

*(tf::AbstractPeriodFrame, n::Int) = typeof(tf)(tf.period * n, tf.boundary)
*(n::Int, tf::AbstractPeriodFrame) = *(tf, n)

end # module


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

name = "LCIO_Julia_Wrapper"
version = v"0.8.1"

# Collection of sources required to build LCIOWrapBuilder
sources = [
	GitSource("https://github.com/jstrube/LCIO_Julia_Wrapper.git", "0cebd2e94e739e086ee49b2fcce4839786a39b0a"),
	ArchiveSource("https://julialang-s3.julialang.org/bin/linux/x64/1.3/julia-1.3.1-linux-x86_64.tar.gz", "faa707c8343780a6fe5eaf13490355e8190acf8e2c189b9e7ecbddb0fa2643ad"; unpack_target="julia-x86_64-linux-gnu"),
        ArchiveSource("https://github.com/Gnimuc/JuliaBuilder/releases/download/v1.3.0/julia-1.3.0-x86_64-apple-darwin14.tar.gz", "f2e5359f03314656c06e2a0a28a497f62e78f027dbe7f5155a5710b4914439b1"; unpack_target="julia-x86_64-apple-darwin14"),
]

# Bash recipe for building across all platforms
script = raw"""
case "$target" in
	arm-linux-gnueabihf|x86_64-linux-gnu)
        Julia_PREFIX=${WORKSPACE}/srcdir/julia-$target/julia-1.3.1
        ;;
    x86_64-apple-darwin14|x86_64-w64-mingw32)
        Julia_PREFIX=${WORKSPACE}/srcdir/julia-$target/juliabin
        ;;
esac

cd ${WORKSPACE}/srcdir/LCIO_Julia_Wrapper/
mkdir build && cd build
cmake -DJulia_PREFIX=${Julia_PREFIX} -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_FIND_ROOT_PATH=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
make install
install_license $WORKSPACE/srcdir/LCIO_Julia_Wrapper/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = Platform[
    Linux(:x86_64, libc=:glibc),
    MacOS(:x86_64)
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("liblciowrap", :lciowrap)
]

# Dependencies that must be installed before this package can be built
dependencies = [
	Dependency("libcxxwrap_julia_jll"),
	Dependency("LCIO_jll")
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"7")

max_stars_repo_name: tbeason/Yggdrasil

max_issues_repo_name: tbeason/Yggdrasil

max_forks_repo_name: tbeason/Yggdrasil


text:
function query(node::TrapezoidalSearchNode, edge::E)
    node.condition[1] == 0 && return node.data
    p = edge.p
    if node.condition[1] == 1
        if p.x < node.condition[2].p.x
            # println("Point left: ", node.condition[2].p)
            return query(node.left, edge)
        else
            # println("Point right: ", node.condition[2].p)
            return query(node.right, edge)
        end
    else
        side = sideof(node.condition[2], p)
        if side == 1
            # println("Edge left: ", node.condition[2])
            return query(node.left, edge)
        elseif side == -1
            # println("Edge right: ", node.condition[2])
            return query(node.right, edge)
        else
            side = sideof(node.condition[2], edge.q)
            if side == 1
                # println("Edge left (2): ", node.condition[2])
                return query(node.left, edge)
            elseif side == -1
                # println("Edge right (2): ", node.condition[2])
                return query(node.right, edge)
            else
                # TODO: propagate error from here
                # println("You've got to be kidding...")
                return node.data
            end
        end
    end
end

function query(node::TrapezoidalSearchNode, p::V)
    node.condition[1] == 0 && return node.data
    if node.condition[1] == 1
        if p.x < node.condition[2].p.x
            # println("Point left: ", node.condition[2].p)
            return query(node.left, p)
        else
            # println("Point right: ", node.condition[2].p)
            return query(node.right, p)
        end
    else
        side = sideof(node.condition[2], p)
        if side == 1
            # println("Edge left: ", node.condition[2])
            return query(node.left, p)
        else
            # println("Edge right: ", node.condition[2])
            return query(node.right, p)
        end
    end
end

function buildnode(leftp, rightp, top, bottom, face)
    t = Trapezoid(leftp, rightp, top, bottom, face)
    node = TrapezoidalSearchNode(t)
    t.node = node

    return t
end

function buildnode(t)
    node = TrapezoidalSearchNode(t)
    t.node = node

    return t
end

similar(node::TrapezoidalSearchNode{T, T1}) where {T, T1} = TrapezoidalSearchNode{T, T1}(nothing, nothing, (0, zero(E{T}), nothing))

function findadjacent(s, t)
    # If we have only one adjacent trapezoid, return it
    t.rb === nothing && return t.rt
    t.rt === nothing && return t.rb

    # We have two adjacent trapezoids, must choose one of them
    p = lp(s)
    if t.rightp == p
        return sideof(s, t.rb.top.q) == 1 ? t.rb : t.rt
    end

    # General position
    return sideof(s, t.rightp) == 1 ? t.rb : t.rt
end

function initialsplit!(t, s)
    p = lp(s)
    node = t.node

    if p.x > t.leftp.x && p.x < t.rightp.x
        # Transforming node to internal and adding type 1 condition
        t.node = nothing
        node.data = nothing
        node.condition = (1, E(p, p))

        # Split initial trapezoid and add results to the left and right nodes
        # Fake right trapezoid, just to simplify things
        t2 = buildnode(p, t.rightp, t.top, t.bottom, t.face)
        t2.rb = t.rb
        t2.rt = t.rt
        node.right = t2.node

        t.rightp = p
        t = buildnode(t)
        node.left = t.node

        t2.lt = t
        t2.lb = t
        t = t2
    end

    topt = buildnode(p, p, t.top, edge(s), lhs(s))
    bott = buildnode(p, p, edge(s), t.bottom, rhs(s))
    if t.lt !== nothing
        t.lt.rt = topt
        topt.lt = t.lt
    end
    if t.lb !== nothing
        t.lb.rb = bott
        bott.lb = t.lb
    end

    return t, topt, bott
end

# TODO: next couple of functions is redundant
function topsplit!(t, s, topt)
    rightp = t.rightp.x < rp(s).x ? t.rightp : rp(s)

    if topt.top == t.top
        topt.rightp = rightp
    else
        topt2 = buildnode(t.leftp, rightp, t.top, edge(s), lhs(s))
        topt.rb = topt2
        topt2.lb = topt
        topt = topt2

        topt.lt = t.lt
        if topt.lt !== nothing
            topt.lt.rt = topt
        end
    end

    topt.rt = t.rt
    if topt.rt !== nothing
        topt.rt.lt = topt
    end

    return topt
end

function bottomsplit!(t, s, bott)
    rightp = t.rightp.x < rp(s).x ? t.rightp : rp(s)

    if bott.bottom == t.bottom
        bott.rightp = rightp
    else
        bott2 = buildnode(t.leftp, rightp, edge(s), t.bottom, rhs(s))
        bott.rt = bott2
        bott2.lt = bott
        bott = bott2

        bott.lb = t.lb
        if bott.lb != nothing
            bott.lb.rb = bott
        end
    end

    bott.rb = t.rb
    if bott.rb !== nothing
        bott.rb.lb = bott
    end

    return bott
end

function finalsplit!(t, s)
    q = rp(s)
    node = t.node
    node.data = nothing
    t.node = nothing

    if t.rightp.x <= q.x
        return true, findadjacent(s, t)
    end

    if t.rt !== nothing
        t.rt.lt = t
    end
    if t.rb !== nothing
        t.rb.lb = t
    end

    # Introduce rightmost node and split trapezoid vertically
    t.leftp = q
    t = buildnode(t)
    node.left.data.rt = t
    t.lt = node.left.data

    node.right.data.rb = t
    t.lb = node.right.data

    t2 = buildnode(node.left.data)
    node.left.left = t2.node
    node.left.right = node.right
    node.left.condition = node.condition
    node.left.data = nothing

    node.right = t.node
    node.condition = (1, E(q, q))

    return false, t
end

function update!(tree::TrapezoidalSearchNode, s)
    t = query(tree, edge(s))

    t, topt, bott = initialsplit!(t, s)

    splitting = true
    while splitting
        # we are building upper and lower division of the current trapezoid
        topt = topsplit!(t, s, topt)
        bott = bottomsplit!(t, s, bott)

        t.node.left = topt.node
        t.node.right = bott.node
        t.node.condition = (2, edge(s))

        splitting, t = finalsplit!(t, s)
    end
end

function buildsearch(s, lb, rt)
    bb = buildnode(lb, rt, E(V(lb.x, rt.y), rt), E(lb, V(rt.x, lb.y)), -1)
    root = bb.node
    for item in s
        update!(root, item)
    end

    return root
end

function buildsearch(s)
    bb = buildnode(boundingbox(s))
    root = bb.node
    for item in s
        update!(root, item)
    end

    return root
end


text:
include("dependencies_for_runtests.jl")
include("dependencies_for_poisson_solvers.jl")

#####
##### Run pressure solver tests 2
#####

@testset "Poisson solvers 2" begin
    @info "Testing Poisson solvers (vertically stretched grid)..."

    # Vertically stretched topologies to test.
    vs_topos = [
        (Periodic, Periodic, Bounded),
        (Periodic, Bounded,  Bounded),
        (Bounded,  Periodic, Bounded),
        (Bounded,  Bounded,  Bounded),
        (Flat,     Bounded,  Bounded),
        (Flat,     Periodic, Bounded),
        (Bounded,  Flat,     Bounded),
        (Periodic, Flat,     Bounded)
    ]

    for arch in archs, topo in vs_topos
        @testset "Vertically stretched Poisson solver [FACR, $(typeof(arch)), $topo]" begin
            @info "  Testing vertically stretched Poisson solver [FACR, $(typeof(arch)), $topo]..."

            @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 8, 8, 1:8)
            @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 7, 7, 1:7)
            @test vertically_stretched_poisson_solver_correct_answer(Float32, arch, topo, 8, 8, 1:8)

            zF_even = [1, 2, 4, 7, 11, 16, 22, 29, 37]      # Nz = 8
            zF_odd  = [1, 2, 4, 7, 11, 16, 22, 29, 37, 51]  # Nz = 9

            for zF in [zF_even, zF_odd]
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 8,  8, zF)
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 16, 8, zF)
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 8, 16, zF)
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 8, 11, zF)
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 5,  8, zF)
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 7, 13, zF)
            end
        end
    end
end


text:
@testset "Power series integration" begin

    T, t = Nemo.PowerSeriesRing(Nemo.QQ, 300, "t"; model=:capped_absolute)

    @test ps_integrate(zero(T)) == zero(T)

    @test ps_integrate(one(T)) == gen(T)

    @test ps_integrate(gen(T)) == gen(T)^2 * (1 // 2)

    @test ps_integrate(t^3 - t^4 * (1 // 8) + 5 * t^10) == t^4 * (1 // 4) - t^5 * (1 // 40) + t^11 * (5 // 11)

    @test ps_integrate( (1 + t)^1000 ) == truncate((1 + t)^1001 * (1 // 1001) - 1 // 1001, 300)

end


text:
"""
Given a PowerModels data dictionary, following the InfrastructureModels
multi-infrastructure conventions, use PowerModels.jl to solve the
AC OPF problem and return the primal and dual variable values, and
whether the AC OPF solver converged to an Optimal solution.
"""
function run_ac_opf(network_data::Dict; print_level=0, from_py=false,
				    solver=JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL))
	if from_py
		network_data = fix_network_data_dict(network_data)
	end

	pm = PM.instantiate_model(network_data, PM.ACPPowerModel, PM.build_opf)
	JuMP.set_optimizer(pm.model, solver)
	if print_level < 2
		JuMP.set_silent(pm.model)
	end
	results = PM.optimize_model!(pm)

	# Extracting Results
	bus_ids = sort(collect(ids(pm, :bus)))
	bus_ids_str = string.(bus_ids)
	bus_nums = collect(1:length(bus_ids))
	bus_num_ref = Dict(zip(bus_ids_str, bus_nums))

	num_buses = length(bus_ids)

	net_gen = network_data["gen"]
	gen_buses = sort([gen["gen_bus"] for gen in values(net_gen)])
	gen_bus_nums = [bus_num_ref[bus_idx] for bus_idx in string.(gen_buses)]
	num_gens = length(net_gen)

	not_slack(i) = net_bus[string(net_gen[i]["gen_bus"])]["bus_type"] != 3 # bus_type 3: Slack bus # UNUSED
	get_index(i) = net_gen[i]["index"] # UNUSED

	# Gather primal values
	primals = Dict()
	vm = zeros(num_buses)
	va = zeros(num_buses)
	res_bus = results["solution"]["bus"]
	for bus_idx in bus_ids_str
		bus_num = bus_num_ref[bus_idx]
		vm[bus_num] = res_bus[bus_idx]["vm"]
		va[bus_num] = res_bus[bus_idx]["va"]
	end
	v = vm.*exp.(1im*deg2rad.(va))
	vmg = vm[gen_bus_nums]

	primals["vm_bus"] = Array(vm')
	primals["va_bus"] = Array(va')
	primals["vm_gen"] = Array(vmg')
	primals["v_bus"] = Array(transpose(v))

	pg = zeros(num_gens)
	qg = zeros(num_gens)
	res_gen = results["solution"]["gen"]
	for gen_idx in keys(res_gen)
		gen_num = net_gen[gen_idx]["index"]
		pg[gen_num] = res_gen[gen_idx]["pg"]
		qg[gen_num] = res_gen[gen_idx]["qg"]
	end

	primals["p_gen"] = Array(pg')
	primals["q_gen"] = Array(qg')

	net_branch = network_data["branch"]
	num_branches = length(net_branch)
	p_to = zeros(num_branches)
	q_to = zeros(num_branches)
	p_fr = zeros(num_branches)
	q_fr = zeros(num_branches)
	res_branch = results["solution"]["branch"]
	for (idx, branch) in res_branch
		index = network_data["branch"][idx]["index"]
		p_to[index] = branch["pt"]
		q_to[index] = branch["qt"]
		p_fr[index] = branch["pf"]
		q_fr[index] = branch["qf"]
	end

	primals["p_fr"] = Array(p_fr')
	primals["p_to"] = Array(p_to')
	primals["q_fr"] = Array(q_fr')
	primals["q_to"] = Array(q_to')
	primals["total_cost"] = Array([results["objective"]])

	if print_level > 1
		println(string("OPF: ",results["termination_status"]))
	end
	converged = results["termination_status"] in SOLVED_STATUSES

	# Gather dual values
	vm = var(pm, :vm)
	pg = var(pm, :pg)
	qg = var(pm, :qg)
	p = var(pm, :p)
	q = var(pm, :q)

	#TASK: Determine if I need to sort these voltage values?
	v_min = [JuMP.dual(JuMP.LowerBoundRef(vm[i])) for i in bus_ids]
	v_max = [JuMP.dual(JuMP.UpperBoundRef(vm[i])) for i in bus_ids]

	pg_min = zeros(num_gens)
	pg_max = zeros(num_gens)
	qg_min = zeros(num_gens)
	qg_max = zeros(num_gens)
	for gen in values(net_gen)
		gen_idx = gen["index"]
		pg_min[gen_idx] = JuMP.dual(JuMP.LowerBoundRef(pg[gen_idx]))
		pg_max[gen_idx] = JuMP.dual(JuMP.UpperBoundRef(pg[gen_idx]))
		qg_min[gen_idx] = JuMP.dual(JuMP.LowerBoundRef(qg[gen_idx]))
		qg_max[gen_idx] = JuMP.dual(JuMP.UpperBoundRef(qg[gen_idx]))
	end

	p_to_max = zeros(num_branches)
	q_to_max = zeros(num_branches)
	p_fr_max = zeros(num_branches)
	q_fr_max = zeros(num_branches)
	if haskey(collect(values(net_branch))[1], "rate_a")
	for branch in values(net_branch)
		idx = branch["index"]
		fr = branch["f_bus"]
		to = branch["t_bus"]
		to_key = (idx, fr, to)
		fr_key = (idx, to, fr)
		p_to = p[to_key]
		q_to = q[to_key]
		p_fr = p[fr_key]
		q_fr = q[fr_key]

		p_to_max[idx] = JuMP.dual(JuMP.UpperBoundRef(p_to))
		q_to_max[idx] = JuMP.dual(JuMP.UpperBoundRef(q_to))
		p_fr_max[idx] = JuMP.dual(JuMP.UpperBoundRef(p_fr))
		q_fr_max[idx] = JuMP.dual(JuMP.UpperBoundRef(q_fr))
	end
	end

	duals = Dict("v_min" => Array(v_min'),
				 "v_max" => Array(v_max'),
				 "pg_min" => Array(pg_min'),
				 "pg_max" => Array(pg_max'),
				 "qg_min" => Array(qg_min'),
				 "qg_max" => Array(qg_max'),
				 "pto_max" => Array(p_to_max'),
				 "qto_max" => Array(q_to_max'),
				 "pfr_max" => Array(p_fr_max'),
				 "qfr_max" => Array(q_fr_max'),
				 )

	results = (primals, duals)
	return results, converged
end

"Calculates and returns the admittance matrix for the given network data dictionary"
function admittance_matrix(network_data::Dict)
	network_data = fix_network_data_dict(network_data)
	return PM.calc_admittance_matrix(network_data).matrix
end

"""
	Sets the PowerModels data dictionaries types to Dict{String,Any} as Python strips the
	string type which is required by PowerModels functions.
"""
function fix_network_data_dict(network_data)
	# Dicts & subDicts passed from python have type Dict{Any,Any}
	network_data = convert(Dict{String,Any}, network_data)
	for key in keys(network_data)  # Need Dict{String,Any}
	if typeof(network_data[key]) == Dict{Any,Any}
		network_data[key] = convert(Dict{String,Any}, network_data[key])
		for key2 in keys(network_data[key])
		if typeof(network_data[key][key2]) == Dict{Any,Any}
			network_data[key][key2] = convert(Dict{String,Any}, network_data[key][key2])
		end
		end
	end
	end

	return network_data
end


text:
abstract type Graph{V,E} end

function add!(g::Graph{V,E}, u::Int64, v::Int64) where {V,E<:Edge}
    add!(g, E(u,v))
end

"""

This is why cities cannot have the same name.
"""
function add!(g::Graph{V,E}, u::V, v::V) where {V,E<:Edge}
    add!(g, index_of(g, u), index_of(g, v))
end

function add!(g::Graph{V,E}, edge::Edge) where {V,E<:Edge}
    edges = edges_of(g, edge[1]) #g.edges_lists[edge.u]
    push!(edges, edge)

    edges = edges_of(g, edge[2]) #g.edges_lists[edge.v]
    push!(edges, reverse(edge))

    nothing
end

"""
    Base.append!(g::Graph{V,E}, vertices) where {V,E<:Edge}

Append vertices to a graph.  Creates empty edge lists for each vertex.
"""
function Base.append!(g::Graph{V,E}, vertices) where {V,E<:Edge} # is <:Edge necessary, could be more generic later
    append!(g.vertices, vertices)

    # One list of edges for each vertex
    for i in 1:length(vertices)
        push!(g.edges_lists, Vector{Vector{E}}())
    end
end

"""
    edges_of(g::Graph{V,E}, v::V) where {V,E<:Edge}

Get the edges of a node (vertex).

# Arguments
- `g`: Graph
- `v`: Node (vertex)

"""
function edges_of(g::Graph{V,E}, v::V) where {V,E<:Edge}
    index = index_of(g, v)
    return edges_of(g, index)
end

"""
    edges_of(g::Graph{V,E}, index::Int64)  where {V,E<:Edge}

Get the edges of a node (vertex).

# Arguments
- `g`: Graph
- `index`: Index of the node (vertex)

"""
function edges_of(g::Graph{V,E}, index::Int64)  where {V,E<:Edge}
    return g.edges_lists[index]
end

"""

This is why cities cannot have the same name.
"""
function index_of(g::Graph{V,E}, v::V) where {V,E<:Edge}
    findfirst(x->x==v, g.vertices)
end

function neighbor_of(g::Graph{V,E}, v::V) where {V,E<:Edge}
    index = index_of(g, v)
    return neighbor_of(g, index)
end

function neighbor_of(g::Graph{V,E}, index::Int64) where {V,E<:Edge}
    edges = g.edges_lists[index]

    neighbors = Vector{V}()
    for edge in edges
        push!(neighbors, vertex_at(g, edge[2]))
    end

    return neighbors
end

function vertex_at(g::Graph{V,E}, index::Int64) where {V,E<:Edge}
    return g.vertices[index]
end

function vertex_at(g::Graph{V,E}, v::V) where {V,E<:Edge}
    return g.vertices[index_of(v)]
end

function Base.show(io::IO, g::Graph)
    for v in g.vertices
        print(io, "$v -> [")
        neighbors = neighbor_of(g,v)
        for neighbor in neighbor_of(g,v)
            print(io, "$neighbor")

            if neighbor != neighbors[end]
                print(io, ", ")
            end
        end
        print(io, "]")

        #print(io, "$v -> $(neighbor_of(g, v))")

        println(io)
    end
end

"""
    UnweightedGraph{V,E<:Edge} <: Graph{V,E}

An unweighted graph that uses any type of node and `Edge`
"""
struct UnweightedGraph{V,E<:Edge} <: Graph{V,E}
    vertices::Vector{V}            # list of vertices
    edges_lists::Vector{Vector{E}} # edges_lists[i] -> list of edges that connects to vertices[i]

    function UnweightedGraph{V,E}() where {V,E<:Edge} # is E<:Edge necessary
        new(Vector{V}(), Vector{Vector{E}}())
    end
end

function UnweightedGraph{V,E}(v::Vector{V}) where {V,E<:Edge}
    me = UnweightedGraph{V,E}()
    append!(me, v)

    return me
end

"""
    WeightedGraph{V,E<:WeightedEdge} <: Graph{V,E}

A weighted graph that uses any type of node with `WeightedEdge`
"""
struct WeightedGraph{V,E<:WeightedEdge} <: Graph{V,E}
    vertices::Vector{V}            # list of vertices
    edges_lists::Vector{Vector{E}} # edges_lists[i] -> list of edges that connects to vertices[i]

    function WeightedGraph{V,E}() where {V,E<:WeightedEdge} # is E<:Edge necessary
        new(Vector{V}(), Vector{Vector{E}}())
    end
end

function WeightedGraph{V,E}(v::Vector{V}) where {V,E<:WeightedEdge}
    me = WeightedGraph{V,E}()
    append!(me, v)

    return me
end

function add!(g::WeightedGraph{V,E}, u::Int64, v::Int64, w) where {V,E<:WeightedEdge}
    add!(g, E(u,v,w))
end

function add!(g::WeightedGraph{V,E}, u::V, v::V, w) where {V,E<:WeightedEdge}
    add!(g, index_of(g, u), index_of(g, v), w)
end
function Base.show(io::IO, g::WeightedGraph)
    for v in g.vertices
        print(io, "$v -> [")
        edges = edges_of(g, v)
        for edge in edges
            print(io, "($(vertex_at(g, edge[2])), $(weight(edge)))")
            if edge != edges[end]
                print(io, ", ")
            end
        end
        # neighbors = neighbor_of(g,v)
        # for neighbor in neighbor_of(g,v)
        #     print(io, "($neighbor, )")
        #     if neighbor != neighbors[end]
        #         print(io, ", ")
        #     end
        # end
        print(io, "]")

        #print(io, "$v -> $(neighbor_of(g, v))")

        println(io)
    end
end

function Base.print(io::IO, g::WeightedGraph, edges::Vector{WeightedEdge})
    for edge in edges
        println(io, "$(vertex_at(g, edge[1])) $(weight(edge)) > $(vertex_at(g, edge[2]))")
    end
end

function next_points(grid, here)
    locations = Vector{Tuple{Int64,Int64}}()

    hr = here[1] # row index
    hc = here[2] # col index
    rcount = size(grid)[1]
    ccount = size(grid)[2]

    nr = hr + 1 # next row index
    nc = hc
    if (nr <= rcount) && (grid[nr,nc] != BLOCKED) # () is not necessary, but make it clear
        push!(locations, (nr,nc))
    end

    nr = hr - 1
    nc = hc
    if (1 <= nr) && (grid[nr,nc] != BLOCKED)
        push!(locations, (nr,nc))
    end

    nr = hr
    nc = hc + 1
    if (nc <= ccount) && (grid[nr,nc] != BLOCKED)
        push!(locations, (nr,nc))
    end

    nr = hr
    nc = hc - 1
    if (1 <= nc) && (grid[nr,nc] != BLOCKED)
        push!(locations, (nr,nc))
    end

    return locations
end

"""
Deprecated, use next_points
"""
successors(grid, pt) = next_points(grid, pt)

function is_goal(goal, here)
    #return (goal[1] == here[1]) && (goal[2] == here[2])
    #@show goal, here
    return goal == here
end


text:
type TURI
    value::AbstractString
end

function ==(u1::TURI, u2::TURI)
    (u1.value == u2.value)
end

function Base.convert(::Type{URIParser.URI}, turi::TURI)
    URIParser.URI(turi.value)
end

function Base.string(turi::TURI)
    turi.value
end


text:
module TrialDocs

include("Random/new_functions.jl")

export hello, domath, test_args, test_args_kw

"""
    hello(who::String)

Return "Hello, `who`".
"""
hello(who::String) = "Hello, $who"

"""
    domath(x::Number)

Return `x + 5`.
"""
domath(x::Number) = x + 5

end # module

max_stars_repo_name: sandreza/TrialDocs

max_issues_repo_name: sandreza/TrialDocs

max_forks_repo_name: sandreza/TrialDocs


text:
using Gadfly

set_default_plot_size(8inch,3.3inch)

# Issue #1357

p1 = plot(xmin=[.25], xmax=[.75], Geom.band, color=[colorant"red"])
p2 = plot(xmin=[.25], xmax=[.75], Geom.band, color=[colorant"red"],
    layer(Geom.rect, xmin=[0], xmax=[1], ymin=[0], ymax=[1], color=[colorant"blue"]))
p3 = plot(x=[0,10], y=[0,1], Geom.blank,
    layer(xintercept=[1.2inch], Geom.vline(color="black")) )

hstack(p1, p2, p3)


text:
iscompact(d::TypedEndpointsInterval{:closed,:closed}) = true
iscompact(d::TypedEndpointsInterval) = false

isinterval(d::Domain) = false
isinterval(d::AbstractInterval) = true

hash(d::AbstractInterval, h::UInt) =
    hashrec(isleftopen(d), isrightopen(d), leftendpoint(d), rightendpoint(d), h)

Display.object_parentheses(::AbstractInterval) = true

approx_indomain(x, d::AbstractInterval, tolerance) =
    (x >= leftendpoint(d)-tolerance) && (x <= rightendpoint(d)+tolerance)

approx_indomain(x, d::TypedEndpointsInterval{:closed,:closed}, tolerance) =
    (x >= leftendpoint(d)-tolerance) && (x <= rightendpoint(d)+tolerance)
approx_indomain(x, d::TypedEndpointsInterval{:open,:closed}, tolerance) =
    (x > leftendpoint(d)-tolerance) && (x <= rightendpoint(d)+tolerance)
approx_indomain(x, d::TypedEndpointsInterval{:closed,:open}, tolerance) =
    (x >= leftendpoint(d)-tolerance) && (x < rightendpoint(d)+tolerance)
approx_indomain(x, d::TypedEndpointsInterval{:open,:open}, tolerance) =
    (x > leftendpoint(d)-tolerance) && (x < rightendpoint(d)+tolerance)

function point_in_domain(d::AbstractInterval)
    isempty(d) && throw(BoundsError())
    mean(d)
end

# For an interval of integers, try to find an integer point
function point_in_domain(d::AbstractInterval{T}) where {T<:Integer}
    isempty(d) && throw(BoundsError())
    x = round(T, mean(d))
    if x ∈ d
        x
    else
        # the mean is not inside the interval when rounded, this means
        # we have to choose an endpoint
        if isleftclosed(d)
            leftendpoint(d)
        elseif isrightclosed(d)
            rightendpoint(d)
        else
            # the interval is open and contains no integers in its interior
            throw(BoundsError())
        end
    end
end

isapprox(d1::AbstractInterval, d2::AbstractInterval) =
    isapprox(leftendpoint(d1), leftendpoint(d2); atol=default_tolerance(d1)) &&
    isapprox(rightendpoint(d1), rightendpoint(d2); atol=default_tolerance(d1))

boundary(d::AbstractInterval) = Point(leftendpoint(d)) ∪ Point(rightendpoint(d))
corners(d::AbstractInterval) = [leftendpoint(d), rightendpoint(d)]

normal(d::AbstractInterval, x) = (abs(minimum(d)-x) < abs(maximum(d)-x)) ? -one(eltype(d)) : one(eltype(d))

distance_to(d::AbstractInterval, x) = x ∈ d ? zero(eltype(d)) : min(abs(x-supremum(d)), abs(x-infimum(d)))

boundingbox(d::AbstractInterval) = d

volume(d::AbstractInterval) = width(d)

similar_interval(d::AbstractInterval, a, b) = similar_interval(d, promote(a, b)...)

## Some special intervals follow, e.g.:
# - the unit interval [0,1]
# - the 'Chebyshev' interval [-1,1]
# - ...
# Unlike a generic interval, these specific intervals have no data, and only one
# type parameter (T).

"""
The abstract type `FixedInterval` is the supertype of intervals with endpoints
determined by the type, rather than field values. Examples include `UnitInterval`
and `ChebyshevInterval`.
"""
abstract type FixedInterval{L,R,T} <: TypedEndpointsInterval{L,R,T} end
const ClosedFixedInterval{T} = FixedInterval{:closed,:closed,T}

closure(d::AbstractInterval) = ClosedInterval(endpoints(d)...)
closure(d::ClosedInterval) = d

interior(d::AbstractInterval) = OpenInterval(endpoints(d)...)

"""
Return an interval that is similar to the given interval, but with endpoints
`a` and `b` instead.
"""
similar_interval(d::ClosedFixedInterval{T}, a::S, b::S) where {S,T} =
    ClosedInterval{promote_type(float(T),S)}(a, b)

"The closed unit interval [0,1]."
struct UnitInterval{T} <: ClosedFixedInterval{T} end

UnitInterval() = UnitInterval{Float64}()

endpoints(d::UnitInterval{T}) where {T} = (zero(T), one(T))

similardomain(::UnitInterval, ::Type{T}) where {T} = UnitInterval{T}()

"The closed interval [-1,1]."
struct ChebyshevInterval{T} <: ClosedFixedInterval{T}
end

ChebyshevInterval() = ChebyshevInterval{Float64}()

endpoints(d::ChebyshevInterval{T}) where {T} = (-one(T),one(T))

similardomain(::ChebyshevInterval, ::Type{T}) where {T} = ChebyshevInterval{T}()

-(d::ChebyshevInterval) = d

"Map the interval [a,b] to the interval [c,d]."
interval_map(a, b, c, d) = AffineMap((d-c)/(b-a), c - a*(d-c)/(b-a))

canonicaldomain(d::ClosedInterval{T}) where {T} = ChebyshevInterval{float(T)}()
canonicaldomain(d::Interval{L,R,T}) where {L,R,T} = Interval{L,R,float(T)}(-1, 1)
canonicaldomain(d::FixedInterval) = d

mapfrom_canonical(d::Interval{L,R,T}) where {L,R,T} = interval_map(-one(T), one(T), endpoints(d)...)
mapfrom_canonical(d::FixedInterval{L,R,T}) where {L,R,T} = StaticIdentityMap{T}()

mapto(d1::D, d2::D) where {D <: FixedInterval} = identitymap(d1)
mapto(d1::D1, d2::D2) where {D1 <: FixedInterval,D2 <: FixedInterval} =
    interval_map(infimum(d1), supremum(d1), infimum(d2), supremum(d2))

"The half-open positive halfline `[0,∞)`."
struct HalfLine{T} <: FixedInterval{:closed,:open,T} end
HalfLine() = HalfLine{Float64}()

endpoints(d::HalfLine{T}) where {T} = (zero(T), T(Inf))

boundary(d::HalfLine) = Point(leftendpoint(d))

similardomain(::HalfLine, ::Type{T}) where {T} = HalfLine{T}()

# intercept and simplify the definition of IntervalSets
in(x, d::HalfLine) = x >= 0

approx_indomain(x, d::HalfLine, tolerance) = x >= -tolerance

function similar_interval(d::HalfLine{T}, a::S, b::S) where {T,S}
    @assert a == 0
    @assert isinf(b) && b > 0
    HalfLine{promote_type(float(T),S)}()
end

point_in_domain(d::HalfLine) = zero(eltype(d))

"The open negative halfline `(-∞,0)`."
struct NegativeHalfLine{T} <: FixedInterval{:open,:open,T} end
NegativeHalfLine() = NegativeHalfLine{Float64}()

similardomain(::NegativeHalfLine, ::Type{T}) where {T} = NegativeHalfLine{T}()

endpoints(d::NegativeHalfLine{T}) where {T} = (-T(Inf), zero(T))

boundary(d::NegativeHalfLine) = Point(rightendpoint(d))

# intercept and simplify the definition of IntervalSets
in(x, d::NegativeHalfLine) = x < 0

approx_indomain(x, d::NegativeHalfLine, tolerance) = x < tolerance

function similar_interval(d::NegativeHalfLine{T}, a::S, b::S) where {S,T}
    @assert isinf(a) && a < 0
    @assert b == 0
    NegativeHalfLine{promote_type(S,float(T))}()
end

point_in_domain(d::NegativeHalfLine) = -one(eltype(d))

similar_interval(d::Interval{L,R,T}, a::S, b::S) where {L,R,T,S} =
    Interval{L,R,promote_type(float(T),S)}(a, b)

#########################################
# A few set operations with known result
#########################################

# We define an exhaustive list of combinations of the four fixed intervals
# combined above in the routines 'intersect', 'union' and 'setdiff' where
# the output is known explicitly.

# Override the definition of Intervals.jl for FixedInterval's defined here
union(d1::FixedInterval, d2::FixedInterval) = uniondomain(d1, d2)
intersect(d1::FixedInterval, d2::FixedInterval) = intersectdomain(d1, d2)

uniondomain(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval) =
    uniondomain(promote(d1,d2)...)

# type-unstable union of intervals. This function differs from `union` in
# IntervalSets.jl because that one throws an error if the intervals do not overlap
function uniondomain(d1::TypedEndpointsInterval{L1,R1,T}, d2::TypedEndpointsInterval{L2,R2,T}) where {L1,R1,L2,R2,T}
    a1 = leftendpoint(d1)
    b1 = rightendpoint(d1)
    a2 = leftendpoint(d2)
    b2 = rightendpoint(d2)

    # are they empty?
    isempty(d1) && return d2
    isempty(d2) && return d1
    # are they equal?
    d1 == d2 && return d1
    # does one lie within the other?
    a2 > a1 && b2 < b1 && return d1
    a1 > a2 && b1 < b2 && return d2
    # are they disjoint?
    b2 < a1 && return UnionDomain(d1, d2)   # return a UnionDomain for disjoint intervals
    b1 < a2 && return UnionDomain(d1, d2)
    # at this stage they must overlap
    if a1 < a2
        a = a1
        L = L1
    elseif a1 == a2
        a = a1
        L = (L1==L2==:open) ? :open : :closed
    else
        a = a2
        L = L2
    end
    if b1 > b2
        b = b1
        R = R1
    elseif b1 == b2
        b = b1
        R = (R1==R2==:open) ? :open : :closed
    else
        b = b2
        R = R2
    end
    Interval{L,R,T}(a, b)
end

function intersectdomain(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval)
    # go back to the definition of IntervalSets.jl
    d = intersect(d1, d2)
    if isempty(d)
        # but avoid returning an interval like 2..1
        EmptySpace{eltype(d)}()
    elseif leftendpoint(d)==rightendpoint(d)
        # and avoid an interval like 1..1
        Point(leftendpoint(d))
    else
        d
    end
end

==(d1::TypedEndpointsInterval, d2::Point) =
    isclosedset(d1) && (leftendpoint(d1)==rightendpoint(d1)==d2.x)

# intersectdomain(d1::AbstractInterval, d2::AbstractInterval) = intersect(d1, d2)
# uniondomain(d1::AbstractInterval, d2::AbstractInterval) = union(d1, d2)

# Since fixed intervals are fully determined by their type,
# the result of intersect, union or setdiff is always known for two
# domains of the same type.
intersectdomain(d1::D, d2::D) where {D <: FixedInterval} = d1
uniondomain(d1::D, d2::D) where {D <: FixedInterval} = d1
setdiffdomain(d1::D, d2::D) where {D <: FixedInterval} = EmptySpace{eltype(D)}()

# [0,1] ∩ [-1,1] = [0,1]
intersectdomain(d1::UnitInterval{T}, d2::ChebyshevInterval{T}) where {T} = UnitInterval{T}()
intersectdomain(d1::ChebyshevInterval{T}, d2::UnitInterval{T}) where {T} = UnitInterval{T}()
# [0,1] ∩ [0,∞) = [0,1]
intersectdomain(d1::UnitInterval{T}, d2::HalfLine{T}) where {T} = UnitInterval{T}()
intersectdomain(d1::HalfLine{T}, d2::UnitInterval{T}) where {T} = UnitInterval{T}()
# [0,1] ∩ (-∞,0) = {}
intersectdomain(d1::UnitInterval{T}, d2::NegativeHalfLine{T}) where {T} = EmptySpace{T}()
intersectdomain(d1::NegativeHalfLine{T}, d2::UnitInterval{T}) where {T} = EmptySpace{T}()
# [-1,1] ∩ [0,∞) = [0,1]
intersectdomain(d1::ChebyshevInterval{T}, d2::HalfLine{T}) where {T} = UnitInterval{T}()
intersectdomain(d1::HalfLine{T}, d2::ChebyshevInterval{T}) where {T} = UnitInterval{T}()
# [0,∞) ∩ (-∞,0) = {}
intersectdomain(d1::HalfLine{T}, d2::NegativeHalfLine{T}) where {T} = EmptySpace{T}()
intersectdomain(d1::NegativeHalfLine{T}, d2::HalfLine{T}) where {T} = EmptySpace{T}()

# [0,1] ∪ [-1,1] = [-1,1]
uniondomain(d1::UnitInterval{T}, d2::ChebyshevInterval{T}) where {T} = ChebyshevInterval{T}()
uniondomain(d1::ChebyshevInterval{T}, d2::UnitInterval{T}) where {T} = ChebyshevInterval{T}()
# [0,1] ∪ [0,∞) = [0,∞)
uniondomain(d1::UnitInterval{T}, d2::HalfLine{T}) where {T} = HalfLine{T}()
uniondomain(d1::HalfLine{T}, d2::UnitInterval{T}) where {T} = HalfLine{T}()

# (-∞,0) ∪ [0,∞) = (-∞,∞)
# Note: T<:real to ensure that FullSpace{T} is not larger than intended.
uniondomain(d1::NegativeHalfLine{T}, d2::HalfLine{T}) where {T<:Real} = FullSpace{T}()
uniondomain(d1::HalfLine{T}, d2::NegativeHalfLine{T}) where {T<:Real} = FullSpace{T}()

# [0,1] ∖ [-1,1] = {}
setdiffdomain(d1::UnitInterval{T}, d2::ChebyshevInterval{T}) where {T} = EmptySpace{T}()
# [0,1] ∖ [0,∞) = {}
setdiffdomain(d1::UnitInterval{T}, d2::HalfLine{T}) where {T} = EmptySpace{T}()
# [0,1] ∖ (-∞,0) = [0,1]
setdiffdomain(d1::UnitInterval{T}, d2::NegativeHalfLine{T}) where {T} = UnitInterval{T}()
# [-1,1] ∖ (-∞,0) = [0,1]
setdiffdomain(d1::ChebyshevInterval{T}, d2::NegativeHalfLine{T}) where {T} = UnitInterval{T}()
# [0,∞) ∖ (-∞,0) = [0,∞)
setdiffdomain(d1::HalfLine{T}, d2::NegativeHalfLine{T}) where {T} = HalfLine{T}()
# (-∞,0) ∖ [0,1] = (-∞,0)
setdiffdomain(d1::NegativeHalfLine{T}, d2::UnitInterval{T}) where {T} = NegativeHalfLine{T}()
# (-∞,0) ∖ [0,∞) = (-∞,0)
setdiffdomain(d1::NegativeHalfLine{T}, d2::HalfLine{T}) where {T} = NegativeHalfLine{T}()

#################################
# Conversions between intervals
#################################
# Avoid depcrecated warning: Warning: Constructors no longer fall back to `convert`.
# example: A constructor `AbstractInterval{Float64}(::IntervalSets.ClosedInterval{Float64})` should be defined instead.

function convert(::Type{UnitInterval{T}}, d::AbstractInterval) where {T}
    endpoints(d) == (0,1) || throw(InexactError(:convert,UnitInterval,d))
    UnitInterval{T}()
end

function convert(::Type{ChebyshevInterval{T}}, d::AbstractInterval) where {T}
    endpoints(d) == (-1,1)|| throw(InexactError(:convert,ChebyshevInterval,d))
    ChebyshevInterval{T}()
end

UnitInterval{T}(d::AbstractInterval) where T = convert(UnitInterval{T}, d)
ChebyshevInterval{T}(d::AbstractInterval) where T = convert(ChebyshevInterval{T}, d)

function show(io::IO, d::ChebyshevInterval)
    print(io, Interval(d))
    print(io, " (Chebyshev)")
end
function show(io::IO, d::UnitInterval)
    print(io, Interval(d))
    print(io, " (Unit)")
end
function show(io::IO, d::HalfLine)
    print(io, Interval(d))
    print(io, " (HalfLine)")
end

########################
# Arithmetic operations
########################

function setdiffdomain(d1::AbstractInterval{T}, d2::AbstractInterval{T}) where T
    a1 = leftendpoint(d1)
    b1 = rightendpoint(d1)
    a2 = leftendpoint(d2)
    b2 = rightendpoint(d2)

    isempty(d1) && return d1
    isempty(d2) && return d1
    b1 < a2 && return d1
    a1 < a2 ≤ b1 ≤ b2 && return (a1..a2)
    a1 < a2 ≤ b2 < b1 && return uniondomain(a1..a2, b2..b1)
    a2 ≤ a1 < b2 < b1 && return (b2..b1)
    a2 ≤ a1 ≤ b1 ≤ b2 && return EmptySpace{T}()

    @assert b2 ≤ a1
    d1
end

setdiffdomain(d1::AbstractInterval, d2::AbstractInterval) =
    setdiffdomain(promote(d1,d2)...)

switch_open_closed(d::AbstractInterval) = d
switch_open_closed(d::Interval{L,R,T}) where {L,R,T} =
    Interval{R,L,T}(leftendpoint(d),rightendpoint(d))

function map_domain(map::AbstractAffineMap{<:Number}, domain::AbstractInterval)
    le = map(leftendpoint(domain))
    re = map(rightendpoint(domain))
    if le<re
        similar_interval(domain,le,re)
    else
        similar_interval(switch_open_closed(domain),re,le)
    end
end

mapped_domain(invmap::AbstractAffineMap{<:Number}, domain::AbstractInterval) =
    map_domain(inverse(invmap), domain)

# Preserve maps when acting on the half line and negative halfline, because
# the map may not always preserve the domain type
map_domain(map::AbstractAffineMap{<:Number}, domain::HalfLine) = MappedDomain(inverse(map), domain)
map_domain(map::AbstractAffineMap{<:Number}, domain::NegativeHalfLine) = MappedDomain(inverse(map), domain)


text:
using EconPDEs

struct BoltonChenWangModel
  r::Float64
  δ::Float64
  A::Float64
  σ::Float64
  θ::Float64
  λ::Float64
  l::Float64
  γ::Float64
  ϕ::Float64
end

function BoltonChenWangModel(r = 0.06, δ = 0.10, A = 0.18, σ = 0.09, θ = 1.5, λ = 0.01, l = 0.9, γ = 0.06, ϕ = 0.01)
  BoltonChenWang(r, δ, A, σ, θ, λ, l , γ, ϕ)
end

function initialize_stategrid(m::BoltonChenWangModel; n = 100)
  OrderedDict(:w => range(0.0, stop = 0.3, length = n))
end

function initialize_y(m::BoltonChenWangModel, stategrid::OrderedDict)
    OrderedDict(:v =>  stategrid[:w])
end

function (m::BoltonChenWangModel)(state::NamedTuple, y::NamedTuple)
  r = m.r; δ = m.δ; A = m.A; σ = m.σ; θ = m.θ; λ = m.λ; l = m.l; γ = m.γ; ϕ = m.ϕ
  w = state.w
  v, vw, vww = y.v, y.vw, y.vww
  i = 1 / θ * (v / vw - w - 1)
  vt = (i - δ) * (v - vw * w) + ((r - λ) * w + A - i - θ * i^2 / 2) * vw + σ^2 / 2 * vww - r * v
  μw = (r - λ) * w + A - i - θ * i^2 / 2 - (i - δ) * w
  return (vt,), (μw, ), (v = v, vw = vw, vww = vww, w = w)
end

m = BoltonChenWangModel()
stategrid = initialize_stategrid(m)
y0 = initialize_y(m, stategrid)
y, result, distance = pdesolve(m, stategrid, y0; bc = OrderedDict(:vw => (1.5, 1.0)))

#========================================================================================

Iterate on boundary conditions until the solution satisfies the conditions given in Bolton Chen Wang (2009)

========================================================================================#

function f(m, x, stategrid, y)
  y, result, distance = pdesolve(m, stategrid, y; bc = OrderedDict(:vw => (x[1], x[2])))
  u1, u2, w = result[:v], result[:vw], result[:w]
  out = zeros(4)
  mi = argmin(abs.(1 + m.γ .- u2))
  out[1] = u2[mi] - 1 - m.γ
  out[2] =  u1[mi] - m.ϕ - (1 + m.γ) * w[mi] - u1[1]
  wi = argmin(abs.(1 .-u2))
  out[3] = u2[wi] - 1
  i = 1 / m.θ * (u1[wi] - w[wi] - 1)
  out[4] = m.r * u1[wi] - (i - m.δ) * (u1[wi] - w[wi]) - ((m.r - m.λ) * w[wi] + m.A - i - m.θ * i^2 / 2)
  return out
end
# using LeastSquaresOptim
# newsol = optimize(x -> f(m, x, stategrid, y0), [1.0,  1.0], Dogleg())
# y, result, distance = pdesolve(m, stategrid, y0; bc = OrderedDict(:vw => (newsol.minimizer[1],newsol.minimizer[2])))


text:
using Test

using Algorithm101
const ALGO = Algorithm101

using DataStructure101
const DS = DataStructure101

@testset "DFS" begin

    g = DS.Graph(10, :di_matrix)

    DS.relate!(g, 10, 9)
    DS.relate!(g, 10, 8)

    DS.relate!(g, 9, 7)
    DS.relate!(g, 9, 6)
    DS.relate!(g, 8, 5)
    DS.relate!(g, 8, 4)

    DS.relate!(g, 7, 3)
    DS.relate!(g, 7, 2)
    DS.relate!(g, 6, 1)

    @test ALGO.dfs(g, 10) == [10, 8, 4, 5, 9, 6, 1, 7, 2, 3]

    n = 100
    city = DS.Graph(100, :w_di_matrix, random=true)
    @test ALGO.dfs(city, 1) == collect(1:100)

end


text:
abstract type AbstractDim{dim} end

struct Dim{dim} <: AbstractDim{dim} end # generic
struct UniaxialStrain{dim} <: AbstractDim{dim} end # 1D uniaxial strain state
struct UniaxialStress{dim} <: AbstractDim{dim} end # 1D uniaxial stress state
struct PlaneStrain{dim} <: AbstractDim{dim} end # 2D plane strain state
struct PlaneStress{dim} <: AbstractDim{dim} end # 2D plane stress state
# constructors without dim parameter
UniaxialStrain() = UniaxialStrain{1}()
UniaxialStress() = UniaxialStress{1}()
PlaneStrain() = PlaneStrain{2}()
PlaneStress() = PlaneStress{2}()

reduce_dim(A::Tensor{1,3}, ::AbstractDim{dim}) where dim = Tensor{1,dim}(i->A[i])
reduce_dim(A::Tensor{2,3}, ::AbstractDim{dim}) where dim = Tensor{2,dim}((i,j)->A[i,j])
reduce_dim(A::Tensor{4,3}, ::AbstractDim{dim}) where dim = Tensor{4,dim}((i,j,k,l)->A[i,j,k,l])
reduce_dim(A::SymmetricTensor{2,3}, ::AbstractDim{dim}) where dim = SymmetricTensor{2,dim}((i,j)->A[i,j])
reduce_dim(A::SymmetricTensor{4,3}, ::AbstractDim{dim}) where dim = SymmetricTensor{4,dim}((i,j,k,l)->A[i,j,k,l])

increase_dim(A::Tensor{1,dim,T}) where {dim, T} = Tensor{1,3}(i->(i <= dim ? A[i] : zero(T)))
increase_dim(A::Tensor{2,dim,T}) where {dim, T} = Tensor{2,3}((i,j)->(i <= dim && j <= dim ? A[i,j] : zero(T)))
increase_dim(A::SymmetricTensor{2,dim,T}) where {dim, T} = SymmetricTensor{2,3}((i,j)->(i <= dim && j <= dim ? A[i,j] : zero(T)))

# restricted strain states
function material_response(
    dim::Union{Dim{d}, UniaxialStrain{d}, PlaneStrain{d}},
    m::AbstractMaterial,
    Δε::AbstractTensor{2,d,T},
    state::AbstractMaterialState,
    Δt = nothing;
    cache = nothing,
    options = Dict{Symbol, Any}(),
    ) where {d,T}

    Δε_3D = increase_dim(Δε)
    σ, dσdε, material_state = material_response(m, Δε_3D, state, Δt; cache=cache, options=options)

    return reduce_dim(σ, dim), reduce_dim(dσdε, dim), material_state
end

# restricted stress states
function material_response(
    dim::Union{UniaxialStress{d}, PlaneStress{d}},
    m::AbstractMaterial,
    Δε::AbstractTensor{2,d,T},
    state::AbstractMaterialState,
    Δt = nothing;
    cache =  get_cache(m),
    options = Dict{Symbol, Any}(),
    ) where {d, T}

    tol = get(options, :plane_stress_tol, 1e-8)
    max_iter = get(options, :plane_stress_max_iter, 10)
    converged = false

    Δε_3D = increase_dim(Δε)

    zero_idxs = get_zero_indices(dim, Δε_3D)
    nonzero_idxs = get_nonzero_indices(dim, Δε_3D)

    Δε_voigt = cache.x_f
    fill!(Δε_voigt, 0.0)
    σ_mandel = cache.F
    J = cache.DF

    for _ in 1:max_iter
        σ, ∂σ∂ε, temp_state = material_response(m, Δε_3D, state, Δt; cache=cache, options=options)
        tomandel!(σ_mandel, σ)
        if norm(view(σ_mandel, zero_idxs)) < tol
            converged = true
            ∂σ∂ε_2D = fromvoigt(SymmetricTensor{4,d}, inv(inv(tovoigt(∂σ∂ε))[nonzero_idxs, nonzero_idxs]))
            return reduce_dim(σ, dim), ∂σ∂ε_2D, temp_state
        end
        tomandel!(J, ∂σ∂ε)
        fill!(Δε_voigt, 0.0)
        Δε_voigt[zero_idxs] = -view(J, zero_idxs, zero_idxs) \ view(σ_mandel, zero_idxs)
        Δε_correction = frommandel(SymmetricTensor{2,3}, Δε_voigt)
        Δε_3D = Δε_3D + Δε_correction
    end
    error("Not converged. Could not find plane stress state.")
end

# out of plane / axis components for restricted stress cases
get_zero_indices(::PlaneStress{2}, ::SymmetricTensor{2,3}) = [3, 4, 5] # for voigt/mandel format, do not use on tensor data!
get_nonzero_indices(::PlaneStress{2}, ::SymmetricTensor{2,3}) = [1, 2, 6] # for voigt/mandel format, do not use on tensor data!
get_zero_indices(::PlaneStress{2}, ::Tensor{2,3}) = [3, 4, 5, 7, 8] # for voigt/mandel format, do not use on tensor data!
get_nonzero_indices(::PlaneStress{2}, ::Tensor{2,3}) = [1, 2, 6, 9] # for voigt/mandel format, do not use on tensor data!

get_zero_indices(::UniaxialStress{1}, ::SymmetricTensor{2,3}) = collect(2:6) # for voigt/mandel format, do not use on tensor data!
get_nonzero_indices(::UniaxialStress{1}, ::SymmetricTensor{2,3}) = [1] # for voigt/mandel format, do not use on tensor data!
get_zero_indices(::UniaxialStress{1}, ::Tensor{2,3}) = collect(2:9) # for voigt/mandel format, do not use on tensor data!
get_nonzero_indices(::UniaxialStress{1}, ::Tensor{2,3}) = [1] # for voigt/mandel format, do not use on tensor data!

# fallback in case there is no cache defined
struct PlaneStressCache{TF, TDF, TX}
    F::TF
    DF::TDF
    x_f::TX
end

# generic fallback, Materials without field σ need to define it
get_stress_type(state::AbstractMaterialState) = typeof(state.σ)

# fallback for optional cache argument
function get_cache(m::AbstractMaterial)
    state = initial_material_state(m)
    stress_type = get_stress_type(state)
    T = eltype(stress_type)
    M = Tensors.n_components(Tensors.get_base(stress_type))
    cache = PlaneStressCache(Vector{T}(undef,M), Matrix{T}(undef,M,M), Vector{T}(undef,M))
    return cache
end