text:
"""
```
segments                = felzenszwalb(img, k, [min_size])
index_map, num_segments = felzenszwalb(edges, num_vertices, k, [min_size])
```

Segments an image using Felzenszwalb's graph-based algorithm. The function can be used in either of two ways -

1. `segments = felzenszwalb(img, k, [min_size])`

Segments an image using Felzenszwalb's segmentation algorithm and returns the result as `SegmentedImage`. The algorithm uses
euclidean distance in color space as edge weights for the region adjacency graph.

Parameters:
-    img            = input image
-    k              = Threshold for region merging step. Larger threshold will result in bigger segments.
-    min_size       = Minimum segment size

2. `index_map, num_segments = felzenszwalb(edges, num_vertices, k, [min_size])`

Segments an image represented as Region Adjacency Graph(RAG) using Felzenszwalb's segmentation algorithm. Each pixel/region
 corresponds to a node in the graph and weights on each edge measure the dissimilarity between pixels.
The function returns the number of segments and index mapping from nodes of the RAG to segments.

Parameters:
-    edges          = Array of edges in RAG. Each edge is represented as `ImageEdge`.
-    num_vertices   = Number of vertices in RAG
-    k              = Threshold for region merging step. Larger threshold will result in bigger segments.
-    min_size       = Minimum segment size

"""
function felzenszwalb(edges::Array{ImageEdge}, num_vertices::Int, k::Real, min_size::Int = 0)

    num_edges = length(edges)
    G = IntDisjointSets(num_vertices)
    set_size = ones(num_vertices)
    threshold = fill(convert(Float64,k), num_vertices)

    sort!(edges, lt = (x,y)->(x.weight<y.weight))

    for edge in edges
        w = edge.weight
        a = find_root!(G, edge.index1)
        b = find_root!(G, edge.index2)
        if a!=b
            if w <= min(threshold[a], threshold[b])
                merged_root = union!(G, a, b)
                set_size[merged_root] = set_size[a] + set_size[b]
                threshold[merged_root] = w + k/set_size[merged_root]
            end
        end
    end

    #merge small segments
    for edge in edges
        a = find_root!(G, edge.index1)
        b = find_root!(G, edge.index2)
        if a!=b && (set_size[a] < min_size || set_size[b] < min_size)
            union!(G, a, b)
        end
    end

    segments = OrderedSet()
    for i in 1:num_vertices
        push!(segments, find_root!(G, i))
    end

    num_sets = length(segments)
    segments2index = Dict{Int, Int}()
    for (i, s) in enumerate(segments)
        segments2index[s] = i
    end

    index_map = Array{Int}(undef, num_vertices)
    for i in 1:num_vertices
        index_map[i] = segments2index[find_root!(G, i)]
    end

    return index_map, num_sets
end

meantype(::Type{T}) where T = typeof(zero(accum_type(T))/2)

function felzenszwalb(img::AbstractArray{T, 2}, k::Real, min_size::Int = 0) where T<:Union{Real,Color}

    rows, cols = size(img)
    num_vertices = rows*cols
    num_edges = 4*rows*cols - 3*rows - 3*cols + 2
    edges = Array{ImageEdge}(undef, num_edges)

    R = CartesianIndices(size(img))
    I1, Iend = first(R), last(R)
    num = 1
    for I in R
        for J in CartesianIndices(_colon(max(I1, I-I1), min(Iend, I+I1)))
            if I >= J
                continue
            end
            edges[num] = ImageEdge((I[2]-1)*rows+I[1], (J[2]-1)*rows+J[1], sqrt(sum(abs2,(img[I])-meantype(T)(img[J]))))
            num += 1
        end
    end

    index_map, num_segments = felzenszwalb(edges, num_vertices, k, min_size)

    result              = similar(img, Int)
    labels              = Array(1:num_segments)
    region_means        = Dict{Int, meantype(T)}()
    region_pix_count    = Dict{Int, Int}()

    for j in axes(img, 2)
        for i in axes(img, 1)
            result[i, j] = index_map[(j-1)*rows+i]
            region_pix_count[result[i,j]] = get(region_pix_count, result[i, j], 0) + 1
            region_means[result[i,j]] = get(region_means, result[i,j], zero(meantype(T))) + (img[i, j] - get(region_means, result[i,j], zero(meantype(T))))/region_pix_count[result[i,j]]
        end
    end

    return SegmentedImage(result, labels, region_means, region_pix_count)
end


text:
"""
find_missing(filename, opening_char, closing_char)

Locates missing brackets etc. in a Latex document.
Example: `find_missing(thesis.tex, '{', '}')
Does not work if opening and closing chars are the same, e.g., \$ \$
"""
function find_missing(filename, oc::AbstractChar, cc::AbstractChar, opens = [])
    @assert oc != cc "I do not work with matching opening and closing characters"
    (filename[end-3:end] == ".tex")  || @warn "This function was implemented to operate on .tex files. I will try anyway..."
    if !isempty(opens)
        @warn("Entering $filename with non-empty list of opened characters")
        @show opens
    end
    each_texline(filename) do lineno,line,file
        for char in line
            if char == oc
                push!(opens, (file, lineno))
            elseif char == cc
                if isempty(opens)
                    @info("Unexpected closing character found", file, lineno, line)
                    @info("Continuing with a count of 0")
                else
                    pop!(opens)
                end
            end
        end
    end
    if !isempty(opens)
        filename, lineno = opens[end]
        @info("Non matched opening character found", filename, lineno)
    end
end

ext: jl


text:
module Electric

using Modia

Pin = Model(
    v = potential,
    i = flow
)

OnePort = Model(
    p = Pin,
    n = Pin,
    equations = :[
        v = p.v - n.v
        0 = p.i + n.i
        i = p.i ],
)

ConstantVoltage = OnePort | Model(
    V = 1.0,
    equations = :[ v = V ] )

VarVoltage = OnePort | Model(
    V = Var(input=true, start=1.0),
    equations = :[v = V] )

Ground = Model(
    p = Pin,
    equations = :[ p.v = 0.0] )

end

max_stars_repo_name: johhell/EPSphasor

max_issues_repo_name: johhell/EPSphasor

max_forks_repo_name: johhell/EPSphasor


text:
using Random, Printf, Statistics
# XPU kernels
macro coords() esc(:( dx*(co1*(nx-2) + ix-0.5)*v1 + dy*(co2*(ny-2) + iy-0.5)*v2 )) end
@parallel_indices (ix,iy) function compute_1!(Yf::Data.Array, v1::Data.Number, v2::Data.Number, a::Data.Number, b::Data.Number, dx::Data.Number, dy::Data.Number, nx::Int, ny::Int, co1::Int, co2::Int)
    if (ix<=size(Yf,1) && iy<=size(Yf,2))  Yf[ix,iy] = Yf[ix,iy] + a*sin( @coords() ) + b*cos( @coords() )  end
    return
end

@parallel function compute_2!(Yf::Data.Array, c::Data.Number)
    @all(Yf) = @all(Yf)*c
    return
end

# 2D Gaussian random field with exponnential covariance
@views function grf2D_expon!(Yf::Data.Array, sf::Data.Number, cl, nh::Int, nx::Int, ny::Int, dx::Data.Number, dy::Data.Number; me=0::Int, co1=0::Int, co2=0::Int, do_reset=true)
    # Resetting the random seed if needed
    if do_reset  Random.seed!(1234)  end
    # Derived numerics
    c = sf/sqrt(nh)
    # Scalar allocations
    ϕ=0.0; k=0.0; d=0.0; θ=0.0; v1=0.0; v2=0.0; a=0.0; b=0.0
    if (me==0) println("Starting 2D RandomField generation (anisotropic exponential covariance function)...") end
    # Loop over nh harmonics
    for ih = 1:nh
        if (ih==501)  global wtime0 = Base.time()  end
        ϕ    = 2.0*pi*rand()
        # Gaussian spectrum
        flag = true
        while flag
            k  = tan(pi*0.5*rand())
            d  = (k*k)/(1.0+(k*k))
            if (rand()<d)  flag = false  end
        end
        θ    = acos(1.0-2.0*rand())
        v1   = k*sin(ϕ)*sin(θ)/cl[1]
        v2   = k*cos(ϕ)*sin(θ)/cl[2]
        a, b = randn(), randn()
        @parallel compute_1!(Yf, v1, v2, a, b, dx, dy, nx, ny, co1, co2)
    end
    @parallel compute_2!(Yf, c)
    # Performance
    wtime    = Base.time() - wtime0
    wtime_it = wtime/(nh-500)                   # Execution time per iteration [s]
    if (me==0) @printf("Total harmonic iters=%d, time=%1.3e sec \n", nh, wtime) end
    return wtime_it
end

# 2D Gaussian random field with Gaussian covariance
@views function grf2D_gauss!(Yf::Data.Array, sf::Data.Number, cl, nh::Int, k_m::Data.Number, nx::Int, ny::Int, dx::Data.Number, dy::Data.Number; me=0::Int, co1=0::Int, co2=0::Int, do_reset=true)
    # Resetting the random seed if needed
    if do_reset  Random.seed!(1234)  end
    # Derived numerics
    c      = sf/sqrt(nh)
    lf     = 2.0*cl/sqrt(pi)
    # Scalar allocations
    ϕ=0.0; k=0.0; d=0.0; θ=0.0; v1=0.0; v2=0.0; a=0.0; b=0.0
    if (me==0) println("Starting 2D RandomField generation (isotropic Gaussian covariance function)...") end
    # Loop over nh harmonics
    for ih = 1:nh
        if (ih==501)  global wtime0 = Base.time()  end
        ϕ    = 2.0*pi*rand()
        # Gaussian spectrum
        flag = true
        while flag
            k  = k_m*rand()
            d  = k*k*exp(-0.5*k*k)
            if (rand()*2.0*exp(-1.0)<d)  flag = false  end
        end
        k    = k/lf*sqrt(2.0)
        θ    = acos(1.0-2.0*rand())
        v1   = k*sin(ϕ)*sin(θ)
        v2   = k*cos(ϕ)*sin(θ)
        a, b = randn(), randn()
        @parallel compute_1!(Yf, v1, v2, a, b, dx, dy, nx, ny, co1, co2)
    end
    @parallel compute_2!(Yf, c)
    # Performance
    wtime    = Base.time() - wtime0
    wtime_it = wtime/(nh-500)                   # Execution time per iteration [s]
    if (me==0) @printf("Total harmonic iters=%d, time=%1.3e sec \n", nh, wtime) end
    return wtime_it
end


text:
"""
    struct Torus{T,CO,TR,TP,TT} <: AbstractVolumePrimitive{T,CO}

Volume primitive describing a [Torus](@ref). It is defined as all points that are
within a given radius to a circle, parallel to the `xy` plane, with constant
radius around a given origin.

## Parametric types
* `T`: Precision type.
* `CO`: Describes whether the surface belongs to the primitive.
    It can be `ClosedPrimitive`, i.e. the surface points belong to the primitive,
    or `OpenPrimitive`, i.e. the surface points do not belong to the primitive.
* `TR`: Type of `r_tube`.
    * `TR == T`: Full tube without cutout (constant radius `r_tube`).
* `TP`: Type of the azimuthial angle `φ`.
    * `TP == Nothing`: Full 2π in `φ`.
* `TT`: Type of the polar angle `θ`.
    * `TT == Nothing`: Full 2π in `θ`.

## Fields
* `r_torus::T`: Distance of the center of the `Torus` to the center of the tube (in m).
* `r_tube::TR`: Radius of the tube of the `Torus` (in m).
* `φ::TP`: Range in azimuthial angle `φ` of the `Torus`.
* `θ::TT`: Range in polar angle `θ` of the `Torus`.
* `origin::CartesianPoint{T}`: The position of the center of the `Torus`.
* `rotation::SMatrix{3,3,T,9}`: Matrix that describes a rotation of the `Torus` around its `origin`.

## Definition in Configuration File

So far, the only `Torus` implemented so far is a `FullTorus`.
A `FullTorus` is defined in the configuration file as part of the `geometry` field
of an object through the field `torus`.

Example definitions of a `FullTorus` looks like this:
```yaml
torus:
  r_torus: 10.0   # => r_torus = 10.0
  r_tube: 2       # => r_tube = 2.0
  phi:
      from: 0.0°
      to: 360.0°  # => φ = nothing
  theta:
      from: 0.0°
      to: 360.0°  # => θ = nothing
```
The fields `phi` and `theta` do not need to defined if they are full 2π.

To define a torus with inner cut-out, use [`CSGDifference`](@ref):
```yaml
difference:
  - torus:
      r_torus: 10.0   # => r_torus = 10.0
      r_tube: 2       # => r_tube = 2.0
  - torus:
      r_torus: 10.0   # => r_torus = 10.0
      r_tube: 1       # => r_tube = 1.0
```
This is a torus with `r_tube` having an inner radius of 1 and an outer radius of 2.

See also [Constructive Solid Geometry (CSG)](@ref).
"""
@with_kw struct Torus{T,CO,TR,TP,TT} <: AbstractVolumePrimitive{T,CO}
    r_torus::T = 1
    r_tube::TR = 1  # (r_tube_in, r_tube_out)
    φ::TP = nothing
    θ::TT = nothing

    origin::CartesianPoint{T} = zero(CartesianPoint{T})
    rotation::SMatrix{3,3,T,9} = one(SMatrix{3, 3, T, 9})
end

Torus{T,CO,TR,TP,TT}( t::Torus{T,CO,TR,TP,TT}; COT = CO,
            origin::CartesianPoint{T} = t.origin,
            rotation::SMatrix{3,3,T,9} = t.rotation) where {T,CO<:Union{ClosedPrimitive, OpenPrimitive},TR,TP,TT} =
    Torus{T,COT,TR,TP,TT}(t.r_torus, t.r_tube, t.φ, t.θ, origin, rotation)

const FullTorus{T,CO} = Torus{T,CO,T,Nothing,Nothing}

function Geometry(::Type{T}, ::Type{Torus}, dict::AbstractDict, input_units::NamedTuple, transformations::Transformations{T}) where {T}
    length_unit = input_units.length
    angle_unit = input_units.angle
    origin = get_origin(T, dict, length_unit)
    rotation = get_rotation(T, dict, angle_unit)

    r_torus = _parse_value(T, dict["r_torus"], length_unit)
    r_tube = _parse_radial_interval(T, dict["r_tube"], length_unit)
    φ = parse_φ_of_primitive(T, dict, angle_unit)
    if !(φ === nothing)
        error("Partial Torus (`φ = φ`) is not yet supported.")
    end
    θ = parse_θ_of_primitive(T, dict, angle_unit)
    if !(θ === nothing)
        error("Partial Torus (`θ = θ`) is not yet supported.")
    end
    if haskey(dict, "z")
        @warn "Deprecation warning: Field `z` for `Torus` is deprecated.
                Use instead (only) `origin` to specify the origin of the primitive.
                There might be a conflict with the possible field `origin`:
                The `z` component of the origin of the primitive is overwritten by the `z`."
        z = _parse_value(T, dict["z"], length_unit)
        origin = CartesianPoint{T}(origin[1], origin[2], z)
    end

    t = if r_tube isa Real
        Torus{T,ClosedPrimitive,typeof(r_tube),typeof(φ),typeof(θ)}(
            r_torus = r_torus, r_tube = r_tube, φ =φ, θ = θ,
            origin = origin, rotation = rotation)
    else
        Torus{T,ClosedPrimitive,typeof(r_tube[2]),typeof(φ),typeof(θ)}(
            r_torus = r_torus, r_tube = r_tube[2], φ =φ, θ = θ,
            origin = origin, rotation = rotation) -
        Torus{T,ClosedPrimitive,typeof(r_tube[1]),typeof(φ),typeof(θ)}(
            r_torus = r_torus, r_tube = r_tube[1], φ =φ, θ = θ,
            origin = origin, rotation = rotation)
    end
    transform(t, transformations)
end

function Dictionary(t::Torus{T})::OrderedDict{String, Any} where {T}
    dict = OrderedDict{String, Any}()
    dict["r_torus"] = t.r_torus
    dict["r_tube"] = t.r_tube # always a Real
    if !isnothing(t.φ) error("Partial Torus (`φ = φ`) is not yet supported.") end
    if !isnothing(t.θ) error("Partial Torus (`θ = θ`) is not yet supported.") end
    if t.origin != zero(CartesianVector{T}) dict["origin"] = t.origin end
    if t.rotation != one(SMatrix{3,3,T,9}) dict["rotation"] = Dictionary(t.rotation) end
    OrderedDict{String, Any}("torus" => dict)
end

function surfaces(t::FullTorus{T,ClosedPrimitive}) where {T}
    tm = FullTorusMantle{T,:inwards}(t.r_torus, t.r_tube, t.φ, t.θ, t.origin, t.rotation)
    (tm, )
end
function surfaces(t::FullTorus{T,OpenPrimitive}) where {T}
    tm = FullTorusMantle{T,:outwards}(t.r_torus, t.r_tube, t.φ, t.θ, t.origin, t.rotation)
    (tm, )
end

function _in(pt::CartesianPoint{T}, t::FullTorus{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    _r = hypot(hypot(pt.x, pt.y) - t.r_torus, pt.z)
    return _r <= t.r_tube + csgtol
end
function _in(pt::CartesianPoint{T}, t::FullTorus{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    _r = hypot(hypot(pt.x, pt.y) - t.r_torus, pt.z)
    return _r < t.r_tube - csgtol
end

# #Constructors
# function Torus(;r_torus = 1, r_tubeMin = 0, r_tubeMax = 1, φMin = 0, φMax = 2π, θMin = 0, θMax = 2π, z = 0)
#     T = float(promote_type(typeof.((r_torus, r_tubeMin, r_tubeMax, φMin, φMax, θMin, θMax, z))...))
#     r_tube = r_tubeMin == 0 ? T(r_tubeMax) : T(r_tubeMin)..T(r_tubeMax)
#     φ = mod(T(φMax) - T(φMin), T(2π)) == 0 ? nothing : T(φMin)..T(φMax)
#     θ = mod(T(θMax) - T(θMin), T(2π)) == 0 ? nothing : T(θMin)..T(θMax)
#     Torus( T, T(r_torus), r_tube, φ, θ, T(z))
# end

# Torus(r_torus, r_tubeMin, r_tubeMax, φMin, φMax, θMin, θMax, z) = Torus(;r_torus = r_torus, r_tubeMin = r_tubeMin, r_tubeMax = r_tubeMax, φMin = φMin, φMax = φMax, θMin = θMin, θMax = θMax, z = z)

# function Torus(r_torus::R1, r_tube::R2, z::TZ) where {R1<:Real, R2<:Real, TZ<:Real}
#     T = float(promote_type(R1, R2, TZ))
#     Torus( T, T(r_torus), T(r_tube), nothing, nothing, T(z))
# end

# function RoundChamfer(r_torus::R1, r_tube::R2, z::TZ) where {R1<:Real, R2<:Real, TZ<:Real}
#     T = float(promote_type(R1, R2, TZ))
#     Torus( T, T(r_torus), T(r_tube), nothing, T(0)..T(π/2), T(z))
# end

# # read-in
function Geometry(T::DataType, ::Type{Torus}, dict::AbstractDict, input_units::NamedTuple)
    length_unit = input_units.length
    angle_unit = input_units.angle
    r_torus::T = _parse_value(T, dict["r_torus"], length_unit)
    r_tube = _parse_radial_interval(T, dict["r_tube"], length_unit)
    φ = parse_φ_of_primitive(T, dict, angle_unit)
    θ = parse_θ_of_primitive(T, dict, angle_unit)
    z::T = ("z" in keys(dict) ? _parse_value(T, dict["z"], length_unit) : T(0))
    Torus(T, r_torus, r_tube, φ, θ, z)
end

# function Dictionary(t::Torus{T}) where {T}
#     dict = OrderedDict{String,Any}()
#     dict["r_torus"] = t.r_torus
#     dict["r_tube"] = typeof(t.r_tube) == T ? t.r_tube : OrderedDict{String,Any}("from" => t.r_tube.left, "to" => t.r_tube.right)
#     if !isnothing(t.φ) dict["phi"] = OrderedDict{String,Any}("from" => t.φ.left, "to" => t.φ.right) end
#     if !isnothing(t.θ) dict["theta"] = OrderedDict{String,Any}("from" => t.θ.left, "to" => t.θ.right) end
#     if t.z != 0 dict["z"] = t.z end
#     OrderedDict{String,Any}("torus" => dict)
# end

# get_r_tube_limits(t::Torus{T}) where {T} = _radial_endpoints(t.r_tube)
# get_φ_limits(t::Torus{T, <:Any, Nothing, <:Any}) where {T} = (T(0), T(2π), true)
# get_φ_limits(t::Torus{T, <:Any, <:AbstractInterval, <:Any}) where {T} = (t.φ.left, t.φ.right, false)

# get_θ_limits(t::Torus{T, <:Any, <:Any, Nothing}) where {T} = (T(0), T(2π), true)
# get_θ_limits(t::Torus{T, <:Any, <:Any, <:AbstractInterval}) where {T} = (t.θ.left, t.θ.right, false)

# function _is_torus_collapsed(t::Torus{T}) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     isapprox(r_tubeMin, r_tubeMax, atol = geom_atol_zero(T))
# end

# function _get_decomposed_surfaces_torus(t::Torus{T, T}) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     AbstractSurfacePrimitive[TorusMantle(t, r_tube = r_tubeMax)]
# end

# function _get_decomposed_surfaces_torus(t::Torus{T, <:AbstractInterval{T}}) where {T}
#     surfaces = AbstractSurfacePrimitive[]
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     if isapprox(r_tubeMin, r_tubeMax, atol = geom_atol_zero(T))
#         push!(surfaces, TorusMantle(t, r_tube = r_tubeMax))
#     else
#         push!(surfaces, TorusMantle(t, r_tube = r_tubeMin), TorusMantle(t, r_tube = r_tubeMax))
#     end
#     surfaces
# end

# get_decomposed_surfaces(t::Torus{T, <:Any, Nothing, Nothing}) where {T} = _get_decomposed_surfaces_torus(t)

# function get_decomposed_surfaces(t::Torus{T, <:Any, <:AbstractInterval, Nothing}) where {T}
#     φMin::T, φMax::T, _ = get_φ_limits(t)
#     surfaces = _get_decomposed_surfaces_torus(t)
#     if !_is_torus_collapsed(t)
#         push!(surfaces, ToroidalAnnulus(t, φ = φMin), ToroidalAnnulus(t, φ = φMax))
#     end
#     surfaces
# end

# function get_decomposed_surfaces(t::Torus{T, <:Any, Nothing, <:AbstractInterval}) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     θMin::T, θMax::T, _ = get_θ_limits(t)
#     surfaces = _get_decomposed_surfaces_torus(t)
#     if !_is_torus_collapsed(t)
#         if r_tubeMin == T(0) && minmax(mod.((θMin, θMax), T(2π))...) == (T(0),T(π))
#             rMin = t.r_torus - r_tubeMax
#             rMax = t.r_torus + r_tubeMax
#             r = rMin == T(0) ? T(rMax) : T(rMin)..T(rMax)
#             push!(surfaces, CylindricalAnnulus(T, r, t.φ, t.z))
#         else
#             for θ in [θMin, θMax]
#                 mod(θ, T(2π)) in [T(0),T(π)] ? push!(surfaces, CylindricalAnnulus(t, θ = θ)) : push!(surfaces, ConeMantle(t, θ = θ))
#             end
#         end
#     end
#     surfaces
# end

# function get_decomposed_surfaces(t::Torus{T, <:Any, <:AbstractInterval, <:AbstractInterval}) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     θMin::T, θMax::T, _ = get_θ_limits(t)
#     φMin::T, φMax::T, _ = get_φ_limits(t)
#     surfaces = _get_decomposed_surfaces_torus(t)
#     if !_is_torus_collapsed(t)
#         if r_tubeMin == T(0) && minmax(mod.((θMin, θMax), T(2π))...) == (T(0),T(π))
#             rMin = t.r_torus - r_tubeMax
#             rMax = t.r_torus + r_tubeMax
#             r = rMin == T(0) ? T(rMax) : T(rMin)..T(rMax)
#             push!(surfaces, CylindricalAnnulus(T, r, t.φ, t.z))
#         else
#             for θ in [θMin, θMax]
#                 mod(θ, T(2π)) in [T(0),T(π)] ? push!(surfaces, CylindricalAnnulus(t, θ = θ)) : push!(surfaces, ConeMantle(t, θ = θ))
#             end
#         end
#         push!(surfaces, ToroidalAnnulus(t, φ = φMin), ToroidalAnnulus(t, φ = φMax))
#     end
#     surfaces
# end

# function sample(t::Torus{T}, step::Real) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     θMin::T, θMax::T, _ = get_θ_limits(t)
#     φMin::T, φMax::T, _ = get_φ_limits(t)
#     samples = [
#         CylindricalPoint{T}(t.r_torus+r_tube*cos(θ),φ,r_tube*sin(θ))
#         for r_tube in r_tubeMin:step:r_tubeMax
#         for θ in (r_tube == 0 ? θMin : θMin:step/r_tube:θMax)
#         for φ in (r_tube == 0 ? φMin : φMin:step/r_tube:φMax)
#     ]
# end

# function sample(t::Torus{T}, Nsamps::NTuple{3,Int} = (2,5,3)) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     θMin::T, θMax::T, _ = get_θ_limits(t)
#     φMin::T, φMax::T, _ = get_φ_limits(t)
#     samples = [
#         CylindricalPoint{T}(t.r_torus+r_tube*cos(θ),φ,r_tube*sin(θ))
#         for r_tube in (Nsamps[1] ≤ 1 ? r_tubeMin : range(r_tubeMin, r_tubeMax, length = Nsamps[1]))
#         for θ in (Nsamps[3] ≤ 1 ? θMin : range(θMin, θMax, length = Nsamps[3]))
#         for φ in (Nsamps[2] ≤ 1 ? φMin : range(φMin, φMax, length = Nsamps[2]))
#     ]
# end


text:
export cyclicProximalPoint
@doc doc"""
    cyclicProximalPoint(M, F, proximalMaps, x)
perform a cyclic proximal point algorithm.
# Input
* `M` – a manifold $\mathcal M$
* `F` – a cost function $F\colon\mathcal M\to\mathbb R$ to minimize
* `proximalMaps` – an Array of proximal maps (`Function`s) `(λ,x) -> y` for the summands of $F$
* `x` – an initial value $x\in\mathcal M$

# Optional
the default values are given in brackets
* `evaluationOrder` – ( [`LinearEvalOrder`](@ref) ) – whether
  to use a randomly permuted sequence ([`FixedRandomEvalOrder`](@ref)), a per
  cycle permuted sequence ([`RandomEvalOrder`](@ref)) or the default linear one.
* `λ` – ( `iter -> 1/iter` ) a function returning the (square summable but not
  summable) sequence of λi
* `stoppingCriterion` – ([`stopWhenAny`](@ref)`(`[`stopAfterIteration`](@ref)`(5000),`[`stopWhenChangeLess`](@ref)`(10.0^-8))`) a [`StoppingCriterion`](@ref).

# Output
* `xOpt` – the resulting (approximately critical) point of gradientDescent
* `record` – (if activated) a String containing the stopping criterion stopping
  reason.
"""
function cyclicProximalPoint(M::Mc,
  F::Function, proximalMaps::Array{Function,N} where N, x0::MP;
  evaluationOrder::EvalOrder = LinearEvalOrder(),
  stoppingCriterion::StoppingCriterion = stopWhenAny( stopAfterIteration(5000), stopWhenChangeLess(10.0^-12)),
  λ = i -> typicalDistance(M)/2/i,
  kwargs... #decorator options
  ) where {Mc <: Manifold, MP <: MPoint}
    p = ProximalProblem(M,F,proximalMaps)
    o = CyclicProximalPointOptions(x0,stoppingCriterion,λ,evaluationOrder)

    o = decorateOptions(o; kwargs...)
    resultO = solve(p,o)
    if hasRecord(resultO)
        return getSolverResult(p,getOptions(resultO)), getRecord(resultO)
    end
    return getSolverResult(p,resultO)
end
function initializeSolver!(p::ProximalProblem, o::CyclicProximalPointOptions)
    c = length(p.proximalMaps)
    o.order = updateOrder(c,0,[1:c...],o.orderType)
end
function doSolverStep!(p::ProximalProblem, o::CyclicProximalPointOptions, iter)
    c = length(p.proximalMaps)
    λi = o.λ(iter)
    for k=o.order
        o.x = getProximalMap(p,λi,o.x,k)
    end
    o.order = updateOrder(c,iter,o.order,o.orderType)
end
getSolverResult(p::ProximalProblem, o::CyclicProximalPointOptions) = o.x
updateOrder(n,i,o,::LinearEvalOrder) = o
updateOrder(n,i,o,::RandomEvalOrder) = collect(1:n)[randperm(length(o))]
updateOrder(n,i,o,::FixedRandomEvalOrder) = (i==0) ? collect(1:n)[randperm(length(o))] : o


text:
# StdString

eval(cxxnewfile("StdString.cxx", """
    #include <memory>
    #include <string>
    #include <utility>

    static_assert(sizeof(bool) == 1, "");
    """))

################################################################################

Base.codepoint(ch::StdChar) = ch.cxx

Base.print(io::IO, ch::StdChar) = print(io, Char(ch))
Base.show(io::IO, ch::StdChar) = show(io, Char(ch))

################################################################################

function generate(::Type{StdString})
    eval(cxxfunction(FnName(:RefStdString_new, "std_string_new", libSTL),
                     FnResult(Ptr{StdString}, "std::string *", RefStdString, expr -> :(RefStdString($expr))), FnArg[],
                     "return new std::string;"))
    eval(cxxfunction(FnName(:RefStdString_new, "std_string_new_const_char_ptr_std_size_t", libSTL),
                     FnResult(Ptr{StdString}, "std::string *", RefStdString, expr -> :(RefStdString($expr))),
                     [FnArg(:ptr, Ptr{Cchar}, "ptr", "const char *", AbstractString, identity),
                      FnArg(:size, Csize_t, "size", "std::size_t", Integer, identity)], "return new std::string(ptr, size);"))

    eval(cxxfunction(FnName(:RefStdString_delete, "std_string_delete", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:str, Ptr{StdString}, "str", "std::string * restrict", RefStdString, identity)], "delete str;"))

    eval(cxxfunction(FnName(:(Base.copy), "std_string_copy", libSTL),
                     FnResult(Ptr{StdString}, "std::string *", RefStdString, expr -> :(RefStdString($expr))),
                     [FnArg(:str, Ptr{StdString}, "str", "std::string * restrict", RefStdString, identity)],
                     "return new std::string(*str);"))

    # GCStdString

    eval(cxxfunction(FnName(:GCStdString_sizeof, "std_string_sizeof", libSTL), FnResult(Csize_t, "std::size_t"), FnArg[],
                     "return sizeof(std::string);"))
    push!(selftests, () -> @assert GCStdString_sizeof() <= GCStdString_size)

    eval(cxxfunction(FnName(:GCStdString_construct, "std_string_construct", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{StdString}, "ptr", "void *", GCStdString, identity)], "new(ptr) std::string;"))
    eval(cxxfunction(FnName(:GCStdString_construct, "std_string_construct_const_char_ptr_std_size_t", libSTL),
                     FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{StdString}, "ptr", "void *", GCStdString, identity),
                      FnArg(:str, Ptr{Cchar}, "str", "const char *", AbstractString, identity),
                      FnArg(:size, Csize_t, "size", "std::size_t", Integer, identity)], "new(ptr) std::string(str, size);"))

    eval(cxxfunction(FnName(:GCStdString_destruct, "std_string_destruct", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{StdString}, "ptr", "std::string * restrict", GCStdString, identity)],
                     "ptr->~basic_string();"))

    eval(cxxfunction(FnName(:GCStdString_copy_construct, "std_string_copy_construct", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{StdString}, "ptr", "std::string * restrict", GCStdString, identity),
                      FnArg(:str, Ptr{StdString}, "str", "const std::string * restrict", GCStdString, identity)],
                     "new(ptr) std::string(*str);"))

    # SharedStdString

    eval(cxxfunction(FnName(:SharedStdString_sizeof, "std_shared_ptr_std_string_sizeof", libSTL), FnResult(Csize_t, "std::size_t"),
                     FnArg[], "return sizeof(std::shared_ptr<std::string>);"))
    push!(selftests, () -> @assert SharedStdString_sizeof() <= SharedStdString_size)

    eval(cxxfunction(FnName(:SharedStdString_construct, "std_shared_ptr_std_string_placement_new", libSTL),
                     FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "void *", SharedStdString, expr -> :(pointer_from_objref($expr)))],
                     """
                     auto res = new(ptr) std::shared_ptr<std::string>;
                     *res = std::make_shared<std::string>();
                     """))
    eval(cxxfunction(FnName(:SharedStdString_construct, "std_shared_ptr_std_string_placement_new_const_char_ptr_std_size_t",
                            libSTL), FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "void *", SharedStdString, expr -> :(pointer_from_objref($expr))),
                      FnArg(:str, Ptr{Cchar}, "str", "const char *", AbstractString, identity),
                      FnArg(:size, Csize_t, "size", "std::size_t", Integer, identity)],
                     """
                     auto res = new(ptr) std::shared_ptr<std::string>;
                     *res = std::make_shared<std::string>(str, size);
                     """))

    eval(cxxfunction(FnName(:SharedStdString_destruct, "std_shared_ptr_std_string_placement_delete", libSTL),
                     FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "std::shared_ptr<std::string> * restrict", SharedStdString,
                            expr -> :(pointer_from_objref($expr)))], "ptr->~shared_ptr();"))

    eval(cxxfunction(FnName(:SharedStdString_copy_construct, "std_shared_ptr_std_string_placement_copy", libSTL),
                     FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "void *", SharedStdString, expr -> :(pointer_from_objref($expr))),
                      FnArg(:str, Ptr{StdString}, "str", "const std::string * restrict", SharedStdString, identity)],
                     """
                     auto res = new(ptr) std::shared_ptr<std::string>;
                     *res = std::make_shared<std::string>(*str);
                     """))

    #TODO eval(cxxfunction(FnName(:SharedStdString_get, "std_shared_ptr_std_string_get", libSTL),
    #TODO                  FnResult(Ptr{StdString}, "std::string *"),
    #TODO                  [FnArg(:ptr, Ptr{Cvoid}, "ptr", "std::shared_ptr<std::string> * restrict", SharedStdString,
    #TODO                         expr -> :(pointer_from_objref($expr)))], "return ptr->get();"))
    eval(cxxfunction(FnName(:SharedStdString_get, "std_shared_ptr_std_string_get", libSTL),
                     FnResult(Ptr{StdString}, "std::string *"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "std::shared_ptr<std::string> * restrict", SharedStdString,
                            expr -> :(pointer_from_objref($expr)))], "return ptr->get();"))

    # StdString

    eval(cxxfunction(FnName(:(Base.length), "std_string_length", libSTL),
                     FnResult(Csize_t, "std::size_t", Int, expr -> :(convert(Int, $expr))),
                     [FnArg(:str, Ptr{StdString}, "str", "const std::string * restrict", StdString, identity)],
                     "return str->size();"))

    eval(cxxfunction(FnName(:(Base.getindex), "std_string_getindex", libSTL),
                     FnResult(Cchar, "char", StdChar, expr -> :(StdChar($expr))),
                     [FnArg(:str, Ptr{StdString}, "str", "const std::string * restrict", StdString, identity),
                      FnArg(:idx, Csize_t, "idx", "std::size_t", Integer, expr -> :($expr - 1))], "return (*str)[idx];"))

    eval(cxxfunction(FnName(:(Base.setindex!), "std_string_setindex_", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:str, Ptr{StdString}, "str", "std::string * restrict", StdString, identity),
                      FnArg(:elt, Cchar, "elt", "char", Union{AbstractChar,Integer,StdChar}, expr -> :(convert(StdChar, $expr))),
                      FnArg(:idx, Csize_t, "idx", "std::size_t", Integer, expr -> :($expr - 1))], "(*str)[idx] = elt;"))

    eval(cxxfunction(FnName(:(Base.:(==)), "std_string_equals", libSTL), FnResult(Bool, "bool"),
                     [FnArg(:str1, Ptr{StdString}, "str1", "const std::string * restrict", StdString, identity),
                      FnArg(:str2, Ptr{StdString}, "str2", "const std::string * restrict", StdString, identity)],
                     "return *str1 == *str2;"))
    eval(cxxfunction(FnName(:(Base.:(<)), "std_string_less", libSTL), FnResult(Bool, "bool"),
                     [FnArg(:str1, Ptr{StdString}, "str1", "const std::string * restrict", StdString, identity),
                      FnArg(:str2, Ptr{StdString}, "str2", "const std::string * restrict", StdString, identity)],
                     "return *str1 < *str2;"))
    eval(cxxfunction(FnName(:(Base.cmp), "std_string_cmp", libSTL),
                     FnResult(Cint, "int", Int, expr -> :($expr == 0 ? 0 : $expr < 0 ? -1 : 1)),
                     [FnArg(:str1, Ptr{StdString}, "str1", "const std::string * restrict", StdString, identity),
                      FnArg(:str2, Ptr{StdString}, "str2", "const std::string * restrict", StdString, identity)],
                     "return str1->compare(*str2);"))

    eval(cxxfunction(FnName(:(Base.string), "std_string_String", libSTL),
                     FnResult(Pair{Ptr{Cchar},Csize_t}, "std::pair<const char *, std::size_t>", String,
                              expr -> :(unsafe_string($expr[1], $expr[2]))),
                     [FnArg(:str, Ptr{StdString}, "str", "const std::string * restrict", StdString, identity)],
                     "return std::pair<const char *, std::size_t>(str->c_str(), str->size());"))

    return nothing
end

generate(StdString)

Base.copy(str::GCStdString) = GCStdString(Base.copy, str)
Base.copy(str::SharedStdString) = SharedStdString(Base.copy, str)

free(str::StdString) = StdString_delete(str)

################################################################################

Base.convert(::Type{String}, str::StdString) = string(str)

Base.isempty(str::StdString) = length(str) == 0
Base.eltype(::Type{<:StdString}) = StdChar

function Base.iterate(str::StdString, pos::Int=1)
    pos > length(str) && return nothing
    return str[pos], pos + 1
end

function Base.print(io::IO, str::StdString)
    for i in 1:length(str)
        print(io, str[i])
    end
    return nothing
end
function Base.show(io::IO, str::StdString)
    print(io, '"')
    for i in 1:length(str)
        ch = str[i]
        if ch ∈ ('"', '\\')
            print(io, '\\', ch)
        elseif !isprint(ch)
            print(io, "\\u", string(Int(ch); base=16, pad=2))
        else
            print(io, ch)
        end
    end
    print(io, '"')
    return nothing
end

Base.ncodeunits(str::StdString) = length(str)
Base.codeunit(::StdString) = Cchar
Base.firstindex(::StdString) = 1
Base.lastindex(str::StdString) = length(str)
Base.thisind(str::StdString, i::Int) = Int(i)
Base.prevind(str::StdString, i::Int) = Int(i) - 1
Base.nextind(str::StdString, i::Int) = Int(i) + 1
Base.isvalid(::StdString, ::Int) = true

################################################################################

function TestAbstractTypes.generator(::Type{StdChar})
    function gen(channel::Channel{StdChar})
        while true
            put!(channel, StdChar(rand(Cchar)))
        end
    end
    return Channel{StdChar}(gen)
end

function TestAbstractTypes.generator(::Type{RefStdString})
    source = generator(String)
    function gen(channel::Channel{RefStdString})
        while true
            put!(channel, RefStdString(take!(source)))
        end
    end
    return Channel{RefStdString}(gen)
end

function TestAbstractTypes.generator(::Type{GCStdString})
    source = generator(String)
    function gen(channel::Channel{GCStdString})
        while true
            put!(channel, GCStdString(take!(source)))
        end
    end
    return Channel{GCStdString}(gen)
end

function TestAbstractTypes.generator(::Type{SharedStdString})
    source = generator(String)
    function gen(channel::Channel{SharedStdString})
        while true
            put!(channel, SharedStdString(take!(source)))
        end
    end
    return Channel{SharedStdString}(gen)
end


text:
module Mainrun
    using Dates
    using InteractiveUtils
    import ..LTK_universe:Universe,show_parameters,make_WdagWmatrix,calc_Action,set_β!,set_βs!,get_β,
                            Wilsonloops_actions,calc_looptrvalues,calc_trainingdata,calc_looptrvalues_site
    import ..Actions:Setup_Gauge_action,Setup_Fermi_action,GaugeActionParam_autogenerator
    import ..Measurements:calc_plaquette,measure_correlator,Measurement,calc_polyakovloop,measure_chiral_cond,calc_topological_charge,
                measurements,Measurement_set
    import  ..MD:md_initialize!,MD_parameters_standard,md!,metropolis_update!,construct_MD_parameters
    import ..System_parameters:Params,print_parameters,parameterloading,Params_set#,parameterloading2
    import ..Print_config:write_config
    import ..Smearing:gradientflow!,calc_fatlink_APE,calc_stout
    import ..ILDG_format:ILDG,load_gaugefield,load_gaugefield!,save_binarydata
    import ..Bridge_format:load_BridgeText!,save_textdata
    import ..Heatbath:heatbath!,overrelaxation!
    import ..Wilsonloops:make_plaq,make_loopforactions,make_plaqloops,make_rectloops,make_polyakovloops
    import ..IOmodule:saveU,loadU,loadU!
    import ..SLMC:SLMC_data,show_effbeta,update_slmcdata!
    import ..Gaugefields:calc_GaugeAction

    import ..Actions:GaugeActionParam_standard,
                    GaugeActionParam,
                    GaugeActionParam_autogenerator
    import ..Verbose_print:println_verbose1,println_verbose2,Verbose_1

    import ..System_parameters:system,actions,md,cg,wilson,staggered,measurement

    function run_LQCD(filenamein::String)
        #if isdemo
        #    params_set = Params_set(Demo.system,Demo.actions,Demo.md,Demo.cg,Demo.wilson,Demo.staggered,Demo.measurement)
        #else
        filename = filenamein
        include(abspath(filename))
        params_set = Params_set(system,actions,md,cg,wilson,staggered,measurement)
        #end

        plaq = run_LQCD(params_set)
    end

    #=
    function run_LQCD()
        parameters = parameterloading()

        load_gaugefield()

        univ = Universe(parameters)

        run_LQCD!(univ,parameters)
        return
    end
    =#

    function run_LQCD(parameters::Params)
        univ = Universe(parameters)
        plaq = run_LQCD!(univ,parameters)

        return plaq
    end

    function run_LQCD(params_set::Params_set)
        parameters = parameterloading(params_set)
        univ = Universe(parameters)
        plaq = run_LQCD!(univ,parameters)

        return  plaq
    end

    function run_LQCD!(univ::Universe,parameters::Params)
        verbose = univ.kind_of_verboselevel
        #println("# ",pwd())
        #println("# ",Dates.now())
        println_verbose1(verbose,"# ",pwd())
        println_verbose1(verbose,"# ",Dates.now())
        versioninfo(verbose)

        #show_parameters(univ)

        mdparams = construct_MD_parameters(parameters)

        measset = Measurement_set(univ,parameters.measuredir,measurement_methods=parameters.measurement_methods)

        #if isdemo
        #run_demo!(parameters,univ,measset)
        #else
        plaq = run_core!(parameters,univ,mdparams,measset)
        #end

    end

    function checkfile(filename,list)
        for name in list

        end
    end

    function run_init_Fileloading!(parameters,univ,mdparams,meas,verbose)
        ildg = nothing
        println_verbose1(verbose,"load U from ",parameters.loadU_dir)
        if parameters.loadU_format == "JLD"
            datatype = "JLD"
            filename_load =  filter(f -> contains(f,".jld2"),readdir("./$(parameters.loadU_dir)"))
            #filename_load =  filter(f -> contains(f,".jld"),readdir("./$(parameters.loadU_dir)"))
        elseif parameters.loadU_format == "ILDG"
            datatype = "ILDG"
            filename_load =  filter(f -> contains(f,"ildg"),readdir("./$(parameters.loadU_dir)"))
        elseif parameters.loadU_format == "BridgeText"
            datatype = "BridgeText"
            filename_load =  filter(f -> contains(f,"txt"),readdir("./$(parameters.loadU_dir)"))
        else
            error("loadU_format should be JLD or ILDG")
        end

        if parameters.loadU_fromfile
            data = readlines("./$(parameters.loadU_dir)/$(parameters.loadU_filename)")
            datafromlist = String[]

            for (i,datai) in enumerate(data)
                havesharp = findfirst('#',datai)
                #println("havesharp: $havesharp")
                if havesharp == 1
                    datau = ""
                elseif havesharp != nothing
                    datau = datai[begin:havesharp-1]
                else
                    datau = datai
                end

                datau = strip(datau)

                if datatype == "JLD"
                    doescontains = contains(datau,".jld2")
                    #doescontains = contains(datau,".jld")
                elseif datatype == "ILDG"
                    doescontains = contains(datau,".ildg")
                elseif datatype == "BridgeText"
                    doescontains = contains(datau,".txt")
                end

                #println(doescontains)
                if doescontains
                    if isfile("./$(parameters.loadU_dir)/$(datau)")
                        push!(datafromlist,datau)
                    else
                        println("Warning! $(datau) does not exist in $(parameters.loadU_dir)! We ignore this.")
                    end
                    #println(datau)

                end
            end
            filename_load = datafromlist
            #println(datafromlist)

        end

        #exit()

        #filename = filter(f -> isfile(f), readdir("./$(parameters.loadU_dir)"))
        #println(filename)
        numfiles = length(filename_load)
        println_verbose1(verbose,"Num of files = $numfiles")
        for file in filename_load
            println_verbose1(verbose,"$file")
        end
        #println("Num of files = $numfiles")
        Nsteps = numfiles-1
        filename_i = filename_load[1]
        if parameters.loadU_format == "JLD"
            loadU!(parameters.loadU_dir*"/"*filename_i,univ.U)
        elseif parameters.loadU_format == "ILDG"
            ildg = ILDG(parameters.loadU_dir*"/"*filename_i)
            i = 1
            load_gaugefield!(univ.U,i,ildg,parameters.L,parameters.NC)
        elseif parameters.loadU_format == "BridgeText"
            load_BridgeText!(parameters.loadU_dir*"/"*filename_i,univ.U,parameters.L,parameters.NC)
        end
        return Nsteps,numfiles,filename_load,ildg

    end

    function print_trainingdata(trs,Sg,Sf)
        for i=1:length(trs)
            print(real(trs[i]),"\t",imag(trs[i]),"\t")
        end
        println(Sg,"\t",Sf,"\t #trainingdata")
    end

    function print_trainingdata(fp,trs,Sg,Sf)
        for i=1:length(trs)
            print(fp,real(trs[i]),"\t",imag(trs[i]),"\t")
        end
        println(fp,Sg,"\t",Sf,"\t #trainingdata")
    end

    function print_trainingdata_site(fp,trs,Sg,Sf)
        return

        numloops,NX,NY,NZ,NT = size(trs)
        println(fp,Sg,"\t",Sf," $numloops $NX $NY $NZ $NT")
        for it = 1:NT
            for iz = 1:NZ
                for iy = 1:NY
                    for ix=1:NX
                        print(fp,"$ix $iy $iz $it ")
                        for i=1:numloops
                            print(fp,real(trs[i,ix,iy,iz,it]),"\t",imag(trs[i,ix,iy,iz,it]),"\t")
                        end
                        println(fp,"\t")

                    end
                end
            end
        end
        #println(fp,Sg,"\t",Sf,"\t #trainingdata_site")

    end

    function run_core!(parameters,univ,mdparams,meas)

        # If an algorithm uses fermion integration (trlog(D+m)),
        # it is flagged
        if parameters.update_method == "IntegratedHMC" ||
                    parameters.update_method == "SLHMC" ||
                    parameters.update_method == "IntegratedHB" ||
                    parameters.update_method == "SLMC"
            isIntegratedFermion = true
        else
            isIntegratedFermion = false
        end
        verbose = univ.kind_of_verboselevel

        if parameters.integratedFermionAction
            loopactions = Wilsonloops_actions(univ)
            trainingfp = open(parameters.training_data_name,"w")
            trainingfp2 = open("trainset.txt","w")
            couplinglist = loopactions.couplinglist
            print(trainingfp,"# ")
            for i=1:loopactions.numloops
                print(trainingfp,"Re($(loopactions.couplinglist[i])) Im($(loopactions.couplinglist[i])) ")
            end
            println(trainingfp,"Sg Sf")
            #println(trainingfp,"#Re(plaq) Im(plaq) Re(rect) Im(rect) Re(polyx) Im(polyy) Re(polyz) Im(polyz) Re(polyt) Im(polyt) Sg Sf")
        end

        Nsteps = parameters.Nsteps
        if parameters.update_method == "Fileloading"
            Nsteps,numfiles,filename_load,ildg = run_init_Fileloading!(parameters,univ,mdparams,meas,verbose)
        elseif parameters.update_method == "SLHMC" || parameters.update_method == "SLMC"
            #slmc_data = SLMC_data(1,univ.NC)
            if typeof(univ.gparam) == GaugeActionParam_autogenerator
                slmc_data = SLMC_data(length(univ.gparam.couplinglist),univ.NC)
            else
                slmc_data = SLMC_data(1,univ.NC)
            end
        end

        numaccepts = 0

        if (parameters.Nthermalization ≤ 0 && parameters.initialtrj == 1) || parameters.update_method == "Fileloading"
            plaq,poly = measurements(0,univ.U,univ,meas;verbose = univ.kind_of_verboselevel) # check consistency of preparation.
        end

        if parameters.integratedFermionAction
            trs,Sg,Sf = calc_trainingdata(loopactions,univ)
            print_trainingdata(trs,Sg,Sf)
            print_trainingdata(trainingfp,trs,Sg,Sf)
            trs_site = calc_looptrvalues_site(loopactions,univ)
            print_trainingdata_site(trainingfp2,trs_site,Sg,Sf)
        end

        if parameters.saveU_format != nothing && parameters.update_method != "Fileloading"
            itrj = 0
            itrjstring = lpad(itrj,8,"0")
            itrjsavecount = 0
            println_verbose1(verbose,"save gaugefields U every $(parameters.saveU_every) trajectory")
            #println("save gaugefields U every $(parameters.saveU_every) trajectory")
        end

        if isIntegratedFermion
            Sfold = nothing
        end

        for itrj=parameters.initialtrj:Nsteps
            println("# itrj = $itrj")
            # Update for different updaters
            # HMC: Hybrid Monte-Carlo
            if parameters.update_method == "HMC"
                Hold = md_initialize!(univ)

                @time Hnew = md!(univ,mdparams)

                if parameters.integratedFermionAction
                    trs,Sg,Sf = calc_trainingdata(loopactions,univ)
                    print_trainingdata(trs,Sg,Sf)
                    print_trainingdata(trainingfp,trs,Sg,Sf)
                    trs_site = calc_looptrvalues_site(loopactions,univ)
                    print_trainingdata_site(trainingfp2,trs_site,Sg,Sf)
                end

                accept = metropolis_update!(univ,Hold,Hnew)
                numaccepts += ifelse(accept,1,0)
                #println("Acceptance $numaccepts/$itrj : $(round(numaccepts*100/itrj)) %")
                println_verbose1(verbose,"Acceptance $numaccepts/$itrj : $(round(numaccepts*100/itrj)) %")

            # Heatbath
            elseif parameters.update_method == "Heatbath"
                @time heatbath!(univ)
                if parameters.useOR
                    for ior=1:parameters.numOR
                        @time overrelaxation!(univ)
                    end
                end

                if parameters.integratedFermionAction
                    trs,Sg,Sf = calc_trainingdata(loopactions,univ)
                    print_trainingdata(trs,Sg,Sf)
                    print_trainingdata(trainingfp,trs,Sg,Sf)
                    trs_site = calc_looptrvalues_site(loopactions,univ)
                    print_trainingdata_site(trainingfp2,trs_site,Sg,Sf)
                end

            # Integrated HMC
            # HMC with S = -tr(log(D+m)), instead of the pseudo-fermins
            elseif parameters.update_method == "IntegratedHMC"
                Sgold = md_initialize!(univ)

                @time Sgnew,Sfnew,Sgold,Sfold = md!(univ,Sfold,Sgold,mdparams)
                Sold = Sgold + Sfold
                Snew = Sgnew + Sfnew

                if parameters.integratedFermionAction
                    trs,Sg,Sf = calc_trainingdata(loopactions,univ)
                    print_trainingdata(trs,Sg,Sf)
                    print_trainingdata(trainingfp,trs,Sg,Sf)
                    trs_site = calc_looptrvalues_site(loopactions,univ)
                    print_trainingdata_site(trainingfp2,trs_site,Sg,Sf)
                end

                accept = metropolis_update!(univ,Sold,Snew)
                numaccepts += ifelse(accept,1,0)
                #println("Acceptance $numaccepts/$itrj : $(round(numaccepts*100/itrj)) %")
                println_verbose1(verbose,"Acceptance $numaccepts/$itrj : $(round(numaccepts*100/itrj)) %")

                Sfold = ifelse(accept,Sfnew,Sfold)

            # File loading
            # This actually loads files, and performs measurements
            elseif parameters.update_method == "Fileloading"
                filename_i = filename_load[itrj+1]
                if parameters.loadU_format == "JLD"
                    loadU!(parameters.loadU_dir*"/"*filename_i,univ.U)
                elseif parameters.loadU_format == "ILDG"
                    ildg = ILDG(parameters.loadU_dir*"/"*filename_i)
                    i = 1
                    load_gaugefield!(univ.U,i,ildg,parameters.L,parameters.NC)
                elseif parameters.loadU_format == "BridgeText"
                    load_BridgeText!(parameters.loadU_dir*"/"*filename_i,univ.U,parameters.L,parameters.NC)
                end

                if parameters.integratedFermionAction
                    trs,Sg,Sf = calc_trainingdata(loopactions,univ)
                    print_trainingdata(trs,Sg,Sf)
                    print_trainingdata(trainingfp,trs,Sg,Sf)
                    trs_site = calc_looptrvalues_site(loopactions,univ)
                    print_trainingdata_site(trainingfp2,trs_site,Sg,Sf)
                end

                #loadU!(parameters.loadU_dir*"/"*filename_i,univ.U)

            # SLHMC: Self-learing Hybrid Monte-Carlo
            # This uses gluonic effective action for MD
            elseif parameters.update_method == "SLHMC"

                Sgold = md_initialize!(univ)

                @time Sgnew,Sfnew,Sgold,Sfold,plaq = md!(univ,Sfold,Sgold,mdparams)
                Sold = Sgold + Sfold
                Snew = Sgnew + Sfnew

            elseif parameters.update_method == "SLMC"
                Sgold = md_initialize!(univ)

                @time Sgnew,Sfnew,Sgeffnew,Sgold,Sfold,Sgeffold = md!(univ,Sfold,Sgold,mdparams)
                Sold = Sgold + Sfold -Sgeffold
                Snew = Sgnew + Sfnew -Sgeffnew

            elseif parameters.update_method == "IntegratedHB"
                Sgold = md_initialize!(univ)

                @time Sgnew,Sfnew,Sgeffnew,Sgold,Sfold,Sgeffold = md!(univ,Sfold,Sgold,mdparams)
                Sold = Sgold + Sfold -Sgeffold
                Snew = Sgnew + Sfnew -Sgeffnew

                if parameters.integratedFermionAction
                    trs,Sg,Sf = calc_trainingdata(loopactions,univ)
                    print_trainingdata(trs,Sg,Sf)
                    print_trainingdata(trainingfp,trs,Sg,Sf)
                    trs_site = calc_looptrvalues_site(loopactions,univ)
                    print_trainingdata_site(trainingfp2,trs_site,Sg,Sf)
                end

                accept = metropolis_update!(univ,Sold,Snew)
                numaccepts += ifelse(accept,1,0)
                #println("Acceptance $numaccepts/$itrj : $(round(numaccepts*100/itrj)) %")
                println_verbose1(verbose,"Acceptance $numaccepts/$itrj : $(round(numaccepts*100/itrj)) %")

                Sfold = ifelse(accept,Sfnew,Sfold)
            end

            if parameters.update_method == "SLHMC" || parameters.update_method == "SLMC"
                S2,plaqetc = calc_GaugeAction(univ)
                plaq = plaqetc
                #S2,plaq = calc_Action(univ)

                if typeof(plaqetc) == Float64
                    plaqetc = [plaq]
                end
                #println(plaq)
                #exit()
                Sf = Sfnew

                outputdata = univ.gparam.β*plaqetc[1]*slmc_data.factor + Sf

                update_slmcdata!(slmc_data,plaqetc,outputdata)
                βeffs,Econst,IsSucs = show_effbeta(slmc_data,univ.gparam)
                println("βeffs = ",βeffs)

                println_verbose1(verbose,"#S = ",outputdata)
                println_verbose1(verbose,"#Estimated Seff = ",Econst + sum(βeffs[:].*plaqetc[:])*slmc_data.factor)
                if IsSucs && itrj ≥ parameters.firstlearn
                    mdparams.βeff = βeffs[:]
                end

                if parameters.integratedFermionAction
                    trs,Sg,Sf = calc_trainingdata(loopactions,univ)
                    print_trainingdata(trs,Sg,Sf)
                    print_trainingdata(trainingfp,trs,Sg,Sf)
                    trs_site = calc_looptrvalues_site(loopactions,univ)
                    print_trainingdata_site(trainingfp2,trs_site,Sg,Sf)
                end

                accept = metropolis_update!(univ,Sold,Snew)
                numaccepts += ifelse(accept,1,0)
                println_verbose1(verbose,"Acceptance $numaccepts/$itrj : $(round(numaccepts*100/itrj)) %")
                #println("Acceptance $numaccepts/$itrj : $(round(numaccepts*100/itrj)) %")
                Sfold = ifelse(accept,Sfnew,Sfold)
            end# update end

            if (itrj ≥ parameters.Nthermalization) || parameters.update_method == "Fileloading"
                plaq,poly = measurements(itrj,univ.U,univ,meas;verbose = univ.kind_of_verboselevel)
            end

            if itrj % parameters.saveU_every == 0 && parameters.saveU_format != nothing && parameters.update_method != "Fileloading"
                itrjsavecount += 1

                itrjstring = lpad(itrj,8,"0")
                #itrjstring = lpad(itrjsavecount,8,"0")
                if parameters.saveU_format == "JLD"
                    #filename = parameters.saveU_dir*"/conf_$(itrjstring).jld"
                    filename = parameters.saveU_dir*"/conf_$(itrjstring).jld2"
                    saveU(filename,univ.U)
                elseif parameters.saveU_format == "ILDG"
                    filename = parameters.saveU_dir*"/conf_$(itrjstring).ildg"
                    save_binarydata(univ.U,filename)
                elseif parameters.saveU_format == "BridgeText"
                    filename = parameters.saveU_dir*"/conf_$(itrjstring).txt"
                    save_textdata(univ.U,filename)
                else
                    error("$(parameters.saveU_format) is not supported")
                end
            end

            println_verbose1(verbose,"-------------------------------------")
            #println("-------------------------------------")
            flush(stdout)
            flush(verbose)
            if parameters.integratedFermionAction
                flush(trainingfp)
                flush(trainingfp2)
            end
        end
        if parameters.integratedFermionAction
            close(trainingfp)
        end
        return plaq
    end

end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    CurveData(data, x, y, z, ...)

Spatial `data` georeferenced with parametric coordinates `x`, `y`, `z`, ...
The `data` argument is a dictionary mapping variable names to Julia arrays
with the actual data.

See also: [`Curve`](@ref)
"""
struct CurveData{T,N} <: AbstractData{T,N}
  data::OrderedDict{Symbol,<:AbstractArray}
  domain::Curve{T,N}

  function CurveData{T,N}(data, domain) where {N,T}
    nvals = [length(array) for array in values(data)]
    @assert all(nvals .== npoints(domain)) "data and domain must have the same number of points"
    new(data, domain)
  end
end

CurveData(data::OrderedDict{Symbol,<:AbstractArray},
          coords::AbstractMatrix{T}) where {T} =
  CurveData{T,size(coords,1)}(data, Curve(coords))

CurveData(data::OrderedDict{Symbol,<:AbstractArray},
          coordarrays::Vararg{<:AbstractVector{T},N}) where {N,T} =
  CurveData{T,N}(data, Curve(coordarrays...))

# ------------
# IO methods
# ------------
function Base.show(io::IO, geodata::CurveData{T,N}) where {N,T}
  npts = npoints(geodata)
  print(io, "$npts CurveData{$T,$N}")
end


text:
mutable struct VoidUpdater <: Updater end
POMDPs.update(::VoidUpdater, ::B, ::Any, ::Any, b=nothing) where B = nothing

mutable struct RandomPolicy{P <: Union{MDP, POMDP}} <: Policy
    rng::AbstractRNG
    problem::P
end
RandomPolicy(problem::Union{POMDP,MDP};
            rng=Random.GLOBAL_RNG) = RandomPolicy(rng, problem)

function POMDPs.action(policy::RandomPolicy, s)
    return rand(policy.rng, actions(policy.problem, s))
end

mutable struct RandomSolver <: Solver
    rng::AbstractRNG
end

RandomSolver(;rng=Base.GLOBAL_RNG) = RandomSolver(rng)
POMDPs.solve(solver::RandomSolver, problem::P) where {P<:Union{POMDP,MDP}} = RandomPolicy(solver.rng, problem)

let
    rng = MersenneTwister(7)

    mdp = GridWorld()
    s = initialstate(mdp, rng)
    a = rand(rng, actions(mdp))
    @inferred generate_sri(mdp, s, a, rng)

    pomdp = TigerPOMDP()
    s = initialstate(pomdp, rng)
    a = rand(rng, actions(pomdp))
    @inferred generate_sori(pomdp, s, a, rng)

    up = VoidUpdater()
    policy = RandomPolicy(rng, pomdp)
    @inferred action_info(policy, s)

    solver = RandomSolver(rng=rng)
    policy, sinfo = solve_info(solver, pomdp)
    @test isa(sinfo, Nothing)

    d = initialstate_distribution(pomdp)
    b = initialize_belief(up, d)
    a = action(policy, b)
    sp, o = generate_so(pomdp, rand(rng, d), a, rng)
    @inferred update_info(up, b, a, o)
end


text:
module TestFind

using Test
using LazyGroupBy

@testset begin
    @test findall.(<(5), grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5])) ==
          Dict(false => [1, 7, 9], true => [3, 4, 8])
    @test findfirst.(<(5), grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5])) ==
          Dict(false => 1, true => 3)
    @test findlast.(<(5), grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5])) ==
          Dict(false => 9, true => 8)
end

end  # module


text:
@doc raw"""
# Lotka-Volterra model in 2D

```math
\begin{aligned}
L (q, \dot{q}) &= \bigg( q_2 + \frac{\log q_2}{q_1} \bigg) \, \dot{q_1} + q_1 \, \dot{q_2} - H(q) , \\
H(q) &= a_1 \, q_1 + a_2 \, q_2 + b_1 \, \log q_1 + b_2 \, \log q_2
\end{aligned}
```

"""
module LotkaVolterra2d

    export lotka_volterra_2d_dg_gauge

    ϑ₁(t, q) = q[2] + log(q[2]) / q[1]
    ϑ₂(t, q) = q[1]

    dϑ₁dx₁(t, q) = - log(q[2]) / q[1]^2
    dϑ₁dx₂(t, q) = 1 + 1 / (q[1] * q[2])

    dϑ₂dx₁(t, q) = one(eltype(q))
    dϑ₂dx₂(t, q) = zero(eltype(q))

    include("lotka_volterra_2d_common.jl")
    include("lotka_volterra_2d_equations.jl")

    function d²ϑ₁d₁d₁(t, q)
        + 2 * log(q[2]) / q[1]^3
    end

    function d²ϑ₁d₁d₂(t, q)
        - 1 / (q[1]^2 * q[2])
    end

    function d²ϑ₁d₂d₁(t, q)
        - 1 / (q[1]^2 * q[2])
    end

    function d²ϑ₁d₂d₂(t, q)
        - 1 / (q[1] * q[2]^2)
    end

    function d²ϑ₂d₁d₁(t, q)
        zero(eltype(q))
    end

    function d²ϑ₂d₁d₂(t, q)
        zero(eltype(q))
    end

    function d²ϑ₂d₂d₁(t, q)
        zero(eltype(q))
    end

    function d²ϑ₂d₂d₂(t, q)
        zero(eltype(q))
    end

    function g̅₁(t, q, v)
        d²ϑ₁d₁d₁(t,q) * q[1] * v[1] + d²ϑ₁d₂d₁(t,q) * q[1] * v[2] + d²ϑ₂d₁d₁(t,q) * q[2] * v[1] + d²ϑ₂d₂d₁(t,q) * q[2] * v[2]
    end

    function g̅₂(t, q, v)
        d²ϑ₁d₁d₂(t,q) * q[1] * v[1] + d²ϑ₁d₂d₂(t,q) * q[1] * v[2] + d²ϑ₂d₁d₂(t,q) * q[2] * v[1] + d²ϑ₂d₂d₂(t,q) * q[2] * v[2]
    end

    function lotka_volterra_2d_ϑ_κ(κ, t, q, v, Θ, params)
        Θ[1] = (1-κ) * ϑ₁(t,q) - κ * f₁(t,q,q)
        Θ[2] = (1-κ) * ϑ₂(t,q) - κ * f₂(t,q,q)
        nothing
    end

    function lotka_volterra_2d_f_κ(κ::Real, t::Real, q::Vector, v::Vector, f::Vector, params)
        f[1] = (1-κ) * f₁(t,q,v) - κ * (g₁(t,q,v) + g̅₁(t,q,v)) - dHd₁(t, q, params)
        f[2] = (1-κ) * f₂(t,q,v) - κ * (g₂(t,q,v) + g̅₂(t,q,v)) - dHd₂(t, q, params)
        nothing
    end

    function lotka_volterra_2d_g_κ(κ::Real, t::Real, q::Vector, v::Vector, g::Vector, params)
        g[1] = (1-κ) * f₁(t,q,v) - κ * (g₁(t,q,v) + g̅₁(t,q,v))
        g[2] = (1-κ) * f₂(t,q,v) - κ * (g₂(t,q,v) + g̅₂(t,q,v))
        nothing
    end

    # function lotka_volterra_2d_g(κ::Real, t::Real, q::Vector, v::Vector, g::Vector)
    #     g[1] = (1-κ) * g₁(t,q,v) - κ * g̅₁(t,q,v) - κ * f₁(t,q,v)
    #     g[2] = (1-κ) * g₂(t,q,v) - κ * g̅₂(t,q,v) - κ * f₂(t,q,v)
    #     nothing
    # end

    function lotka_volterra_2d_dg_gauge(q₀=q₀, p₀=ϑ(0, q₀), κ=0; params=parameters)
        lotka_volterra_2d_ϑ = (t, q, v, p, params) -> lotka_volterra_2d_ϑ_κ(κ, t, q, v, p, params)
        lotka_volterra_2d_f = (t, q, v, f, params) -> lotka_volterra_2d_f_κ(κ, t, q, v, f, params)
        lotka_volterra_2d_g = (t, q, λ, g, params) -> lotka_volterra_2d_g_κ(κ, t, q, λ, g, params)

        IODE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
             lotka_volterra_2d_g, q₀, p₀;
             parameters=params, v̄=lotka_volterra_2d_v)
    end

end


text:
# Hyperbolic network layer from Lensink et al. (2019)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export HyperbolicLayer

"""
    HyperbolicLayer(nx, ny, n_in, batchsize, kernel, stride, pad; action="same", α=1f0, hidden_factor=1)

or

    HyperbolicLayer(W, b, nx, ny, batchsize, stride, pad; action="same", α=1f0)

Create an invertible hyperbolic coupling layer.

*Input*:

 - `nx`, `ny`, `n_in`, `batchsize`: Dimensions of input tensor

 - `kernel`, `stride`, `pad`: Kernel size, stride and padding of the convolutional operator

 - `action`: String that defines whether layer keeps the number of channels fixed (`"same"`),
    increases it by a factor of 4 (`"up"`) or decreased it by a factor of 4 (`"down"`)

 - `W`, `b`: Convolutional weight and bias. `W` has dimensions of `(kernel, kernel, n_in, n_in)`.
   `b` has dimensions of `n_in`.

 - `α`: Step size for second time derivative. Default is 1.

 - `hidden_factor`: Increase the no. of channels by `hidden_factor` in the forward convolution.
    After applying the transpose convolution, the dimensions are back to the input dimensions.

*Output*:

 - `HL`: Invertible hyperbolic coupling layer

 *Usage:*

 - Forward mode: `X_curr, X_new = HL.forward(X_prev, X_curr)`

 - Inverse mode: `X_prev, X_curr = HL.inverse(X_curr, X_new)`

 - Backward mode: `ΔX_prev, ΔX_curr, X_prev, X_curr = HL.backward(ΔX_curr, ΔX_new, X_curr, X_new)`

 *Trainable parameters:*

 - `HL.W`: Convolutional kernel

 - `HL.b`: Bias

 See also: [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct HyperbolicLayer <: NeuralNetLayer
    W::Parameter
    b::Parameter
    α::Float32
    cdims::DenseConvDims
    action::String
end

@Flux.functor HyperbolicLayer

# Constructor 2D
function HyperbolicLayer(nx::Int64, ny::Int64, n_in::Int64, batchsize::Int64, kernel::Int64,
    stride::Int64, pad::Int64; action="same", α=1f0, hidden_factor=1)

    # Set ouput/hidden dimensions
    if action == "same"
        n_out = n_in
    elseif action == "up"
        n_out = Int(n_in/4)
        nx = Int(nx*2)
        ny = Int(ny*2)
    elseif action == "down"
        n_out = Int(n_in*4)
        nx = Int(nx/2)
        ny = Int(ny/2)
    end
    n_hidden = n_out*hidden_factor

    W = Parameter(glorot_uniform(kernel, kernel, n_out, n_hidden))
    b = Parameter(zeros(Float32, n_hidden))

    cdims = DenseConvDims((nx, ny, n_out, batchsize), (kernel, kernel, n_out, n_hidden);
        stride=(stride, stride), padding=(pad, pad))

    return HyperbolicLayer(W, b, α, cdims, action)
end

# Constructor 3D
function HyperbolicLayer(nx::Int64, ny::Int64, nz::Int64, n_in::Int64, batchsize::Int64, kernel::Int64,
    stride::Int64, pad::Int64; action="same", α=1f0, hidden_factor=1)

    # Set ouput/hidden dimensions
    if action == "same"
        n_out = n_in
    elseif action == "up"
        n_out = Int(n_in/8)
        nx = Int(nx*2)
        ny = Int(ny*2)
        nz = Int(nz*2)
    elseif action == "down"
        n_out = Int(n_in*8)
        nx = Int(nx/2)
        ny = Int(ny/2)
        nz = Int(nz/2)
    end
    n_hidden = n_out*hidden_factor

    W = Parameter(glorot_uniform(kernel, kernel, kernel, n_out, n_hidden))
    b = Parameter(zeros(Float32, n_hidden))

    cdims = DenseConvDims((nx, ny, nz, n_out, batchsize), (kernel, kernel, kernel, n_out, n_hidden);
        stride=(stride, stride, stride), padding=(pad, pad, pad))

    return HyperbolicLayer(W, b, α, cdims, action)
end

# Constructor for given weights 2D
function HyperbolicLayer(W::AbstractArray{Float32, 4}, b::AbstractArray{Float32, 1}, nx::Int64, ny::Int64,
    batchsize::Int64, stride::Int64, pad::Int64; action="same", α=1f0)

    kernel, n_in, n_hidden = size(W)[2:4]

    # Set ouput/hidden dimensions
    if action == "same"
        n_out = n_in
    elseif action == "up"
        n_out = Int(n_in/4)
        nx = Int(nx*2)
        ny = Int(ny*2)
    elseif action == "down"
        n_out = Int(n_in*4)
        nx = Int(nx/2)
        ny = Int(ny/2)
    end

    W = Parameter(W)
    b = Parameter(b)

    cdims = DenseConvDims((nx, ny, n_out, batchsize), (kernel, kernel, n_out, n_hidden);
        stride=(stride, stride), padding=(pad, pad))

    return HyperbolicLayer(W, b, α, cdims, action)
end

# Constructor for given weights 3D
function HyperbolicLayer(W::AbstractArray{Float32, 5}, b::AbstractArray{Float32, 1}, nx::Int64, ny::Int64, nz::Int64,
    batchsize::Int64, stride::Int64, pad::Int64; action="same", α=1f0)

    kernel, n_in, n_hidden = size(W)[3:5]

    # Set ouput/hidden dimensions
    if action == "same"
        n_out = n_in
    elseif action == "up"
        n_out = Int(n_in/8)
        nx = Int(nx*2)
        ny = Int(ny*2)
        nz = Int(nz*2)
    elseif action == "down"
        n_out = Int(n_in*8)
        nx = Int(nx/2)
        ny = Int(ny/2)
        nz = Int(nz/2)
    end

    W = Parameter(W)
    b = Parameter(b)

    cdims = DenseConvDims((nx, ny, nz, n_out, batchsize), (kernel, kernel, kernel, n_out, n_hidden);
        stride=(stride, stride, stride), padding=(pad, pad, pad))

    return HyperbolicLayer(W, b, α, cdims, action)
end

# Forward pass
function forward(X_prev_in, X_curr_in, HL::HyperbolicLayer)

    # Change dimensions
    if HL.action == "same"
        X_prev = identity(X_prev_in)
        X_curr = identity(X_curr_in)
    elseif HL.action == "up"
        X_prev = wavelet_unsqueeze(X_prev_in)
        X_curr = wavelet_unsqueeze(X_curr_in)
    elseif HL.action == "down"
        X_prev = wavelet_squeeze(X_prev_in)
        X_curr = wavelet_squeeze(X_curr_in)
    else
        throw("Specified operation not defined.")
    end

    # Symmetric convolution w/ relu activation
    if length(size(X_curr)) == 4
        X_conv = conv(X_curr, HL.W.data, HL.cdims) .+ reshape(HL.b.data, 1, 1, :, 1)
    else
        X_conv = conv(X_curr, HL.W.data, HL.cdims) .+ reshape(HL.b.data, 1, 1, 1, :, 1)
    end
    X_relu = ReLU(X_conv)
    X_convT = -∇conv_data(X_relu, HL.W.data, HL.cdims)

    # Update
    X_new = 2f0*X_curr - X_prev + HL.α*X_convT

    return X_curr, X_new
end

# Inverse pass
function inverse(X_curr, X_new, HL::HyperbolicLayer; save=false)

    # Symmetric convolution w/ relu activation
    if length(size(X_curr)) == 4
        X_conv = conv(X_curr, HL.W.data, HL.cdims) .+ reshape(HL.b.data, 1, 1, :, 1)
    else
        X_conv = conv(X_curr, HL.W.data, HL.cdims) .+ reshape(HL.b.data, 1, 1, 1, :, 1)
    end
    X_relu = ReLU(X_conv)
    X_convT = -∇conv_data(X_relu, HL.W.data, HL.cdims)

    # Update
    X_prev = 2*X_curr - X_new + HL.α*X_convT

    # Change dimensions
    if HL.action == "same"
        X_prev_in = identity(X_prev)
        X_curr_in = identity(X_curr)
    elseif HL.action == "down"
        X_prev_in = wavelet_unsqueeze(X_prev)
        X_curr_in = wavelet_unsqueeze(X_curr)
    elseif HL.action == "up"
        X_prev_in = wavelet_squeeze(X_prev)
        X_curr_in = wavelet_squeeze(X_curr)
    else
        throw("Specified operation not defined.")
    end

    if save == false
        return X_prev_in, X_curr_in
    else
        return X_prev_in, X_curr_in, X_conv, X_relu
    end
end

# Backward pass
function backward(ΔX_curr, ΔX_new, X_curr, X_new, HL::HyperbolicLayer)

    # Recompute forward states
    X_prev_in, X_curr_in, X_conv, X_relu = inverse(X_curr, X_new, HL; save=true)

    # Backpropagate data residual and compute gradients
    ΔX_convT = copy(ΔX_new)
    ΔX_relu = -HL.α*conv(ΔX_convT, HL.W.data, HL.cdims)
    ΔW = -HL.α*∇conv_filter(ΔX_convT, X_relu, HL.cdims)

    ΔX_conv = ReLUgrad(ΔX_relu, X_conv)
    ΔX_curr += ∇conv_data(ΔX_conv, HL.W.data, HL.cdims)
    ΔW += ∇conv_filter(X_curr, ΔX_conv, HL.cdims)
    if length(size(X_curr)) == 4
        Δb = sum(ΔX_conv; dims=[1,2,4])[1,1,:,1]
    else
        Δb = sum(ΔX_conv; dims=[1,2,3,5])[1,1,1,:,1]
    end
    ΔX_curr += 2f0*ΔX_new
    ΔX_prev = -ΔX_new

    # Set gradients
    HL.W.grad = ΔW
    HL.b.grad = Δb

    # Change dimensions
    if HL.action == "same"
        ΔX_prev_in = identity(ΔX_prev)
        ΔX_curr_in = identity(ΔX_curr)
    elseif HL.action == "down"
        ΔX_prev_in = wavelet_unsqueeze(ΔX_prev)
        ΔX_curr_in = wavelet_unsqueeze(ΔX_curr)
    elseif HL.action == "up"
        ΔX_prev_in = wavelet_squeeze(ΔX_prev)
        ΔX_curr_in = wavelet_squeeze(ΔX_curr)
    else
        throw("Specified operation not defined.")
    end

    return ΔX_prev_in, ΔX_curr_in, X_prev_in, X_curr_in
end

# Clear gradients
function clear_grad!(HL::HyperbolicLayer)
    HL.W.grad = nothing
    HL.b.grad = nothing
end

# Put parameters
function put_params!(HL::HyperbolicLayer, Params::Array{Any,1})
    HL.W.data = Params[1].data
    HL.b.data = Params[2].data
end

# Get parameters
get_params(HL::HyperbolicLayer) = [HL.W, HL.b]


text:
# this file is not used in the code. it contains all the testcases but we now use getProblemSpecificParameters instead
# you can however move a testcase from here to the getProblemSpecificParameters file but have to remove the KS. everywhere

function setTestCaseInitialCondition!(KS::KineticSolver,phi::Array{Float64,3})

    if KS.ictype=="linesource"
        KS.SigmaS = 1*ones(KS.ny,KS.nx);
        SigmaA = 0*ones(KS.ny,KS.nx);
        KS.SigmaT = KS.SigmaS+SigmaA;
        KS.x0 = -1.5
        KS.x1 =  1.5
        KS.y0 = -1.5
        KS.y1 =  1.5

        setInitialConditionsGaussian(KS,phi)
    elseif KS.ictype=="checkerboard"
        KS.x0 = 0.0;
        KS.x1 = 7.0;
        KS.y0 = 0.0;
        KS.y1 = 7.0;
        KS.dx = (KS.x1-KS.x0)/(KS.nx-4)
        KS.dy = (KS.y1-KS.y0)/(KS.ny-4)
        KS.SigmaS = zeros(KS.ny,KS.nx);
        KS.SigmaT = zeros(KS.ny,KS.nx);
        phi[:] = 1e-10;
        KS.tEnd = 3.2
        KS.dt = KS.cfl/2*(KS.dx*KS.dy)/(KS.dx+KS.dy)

        KS.nt = ceil(KS.tEnd/KS.dt);

       # KS.rangex = collect(1:KS.nx)+2 # two ghost cells left and right
       # KS.rangey = collect(1:KS.ny)+2 # same
        for j=KS.rangex
            for i=KS.rangey
                y = KS.y0+KS.dy/2 + (i-3)*KS.dy
                x = KS.x0+KS.dx/2 + (j-3)*KS.dx
                cx = ceil(x);
                cy = ceil(y);
                #cx = x;
                #cy = y;
                g = float(ceil((cx+cy)/2.)*2.==(cx+cy))*((1.<cx)&&(cx<7.)&&(1.<cy)&&((cy-2.*abs(cx-4.0))<4.0));
                KS.SigmaS[i,j] = (1-g)*1+g*0;
                KS.SigmaT[i,j] = (1-g)*0+g*10 + KS.SigmaS[i,j]

                if ((3<x)&&(x<4)&&(3<y)&&(y<4))
                        for q=KS.rangequad
                                  KS.Source[q,i,j] = 1.0;#*KS.dx*KS.dy;#*4.0*pi/KS.nquadpoints;
                          end
                    end
            end
        end
    elseif KS.ictype=="checkerboard2"
        KS.x0 = 0.0;
        KS.x1 = 7.0;
        KS.y0 = 0.0;
        KS.y1 = 7.0;
        KS.dx = (KS.x1-KS.x0)/(KS.nx-4)
        KS.dy = (KS.y1-KS.y0)/(KS.ny-4)
        KS.SigmaS = zeros(KS.ny,KS.nx);
        KS.SigmaT = zeros(KS.ny,KS.nx);
        phi[:] = 1e-12;
        KS.tEnd = 3.2
        KS.dt = KS.cfl/2*(KS.dx*KS.dy)/(KS.dx+KS.dy)

        KS.nt = ceil(KS.tEnd/KS.dt);

       # KS.rangex = collect(1:KS.nx)+2 # two ghost cells left and right
       # KS.rangey = collect(1:KS.ny)+2 # same
        for j=KS.rangex
            for i=KS.rangey
                y = KS.y0+KS.dy/2 + (i-3)*KS.dy
                x = KS.x0+KS.dx/2 + (j-3)*KS.dx
                cx = ceil(x);
                cy = ceil(y);
                #cx = x;
                #cy = y;
                g = float(ceil((cx+cy)/2.)*2.==(cx+cy))*((1.<cx)&&(cx<7.)&&(1.<cy)&&((cy-2.*abs(cx-4.0))<4.0));
                KS.SigmaS[i,j] = (1-g)*1+g*0;
                KS.SigmaT[i,j] = (1-g)*0+g*10 + KS.SigmaS[i,j]

                if ((3<x)&&(x<4)&&(3<y)&&(y<4))
                        for q=KS.rangequad
                                  KS.Source[q,i,j] = 1.0;#KS.nquadpoints;#*KS.dx*KS.dy;#*4.0*pi/KS.nquadpoints;
                        end
                        KS.SigmaS[i,j] = 0.0;
                        KS.SigmaT[i,j] = 10.0;
                    end
            end
        end
    elseif KS.ictype=="beam"
        KS.x0 = 0.0;
        KS.x1 = 7.0;
        KS.y0 = 0.0;
        KS.y1 = 7.0;
        KS.dx = (KS.x1-KS.x0)/(KS.nx-4)
        KS.dy = (KS.y1-KS.y0)/(KS.ny-4)
        KS.SigmaS = zeros(KS.ny,KS.nx);
        KS.SigmaT = zeros(KS.ny,KS.nx);
        phi[:] = 1e-8;
        KS.tEnd = 4.0
        KS.dt = KS.cfl/2*(KS.dx*KS.dy)/(KS.dx+KS.dy)

        KS.nt = ceil(KS.tEnd/KS.dt);

        tmpMax = -10.0
        tmpMid = -5.0;
        qMax = -1;
        qMid = -1;
        for q=KS.rangequad
            if KS.Q.pointsxyz[q,1]>= tmpMax
                tmpMid = tmpMax;
                tmpMax = KS.Q.pointsxyz[q,1];
                qMid = qMax;
                qMax = q;
            elseif KS.Q.pointsxyz[q,1]>= tmpMid
                tmpMid = KS.Q.pointsxyz[q,1];
                qMid = q;
            end
        end
        #println(KS.Q.pointsxyz[qMid,1])
        #println(KS.Q.pointsxyz[qMax,1])
        #error("Max directions are ",qMax, " and ", qMid);

       # KS.rangex = collect(1:KS.nx)+2 # two ghost cells left and right
       # KS.rangey = collect(1:KS.ny)+2 # same
        for j=KS.rangex
            for i=KS.rangey
                y = KS.y0+KS.dy/2 + (i-3)*KS.dy
                x = KS.x0+KS.dx/2 + (j-3)*KS.dx
                cx = ceil(x);
                cy = ceil(y);
                #cx = x;
                #cy = y;
                #g = float(ceil((cx+cy)/2.)*2.==(cx+cy))*((1.<cx)&&(cx<7.)&&(1.<cy)&&((cy-2.*abs(cx-4.0))<4.0));
                KS.SigmaS[i,j] = 0.0#(1-g)*1+g*0;
                KS.SigmaT[i,j] = 0.1#(1-g)*0+g*10 + KS.SigmaS[i,j]

                if (3.45<x && x< 3.55 &&i==20)
                    for q=KS.rangequad
                        #if q==qMid || q == qMax
                        if KS.Q.pointsxyz[q,1] > 0.7
                        #if KS.Q.pointsxyz[q,1]>=0
                        #KS.Source[q,i,j] = 10.0;#KS.nquadpoints;#*KS.dx*KS.dy;#*4.0*pi/KS.nquadpoints;
                        phi[q,i,j] = 10.0;
                        end
                    end
                end
            end
        end
    elseif KS.ictype=="circle"
        KS.x0 = -4.0;
        KS.x1 = 4.0;
        KS.y0 = -4.0;
        KS.y1 = 4.0;
        KS.dx = (KS.x1-KS.x0)/(KS.nx-4)
        KS.dy = (KS.y1-KS.y0)/(KS.ny-4)
        KS.SigmaS = zeros(KS.ny,KS.nx);
        KS.SigmaT = zeros(KS.ny,KS.nx);
        phi[:] = 1e-5;
        KS.tEnd = 40
        KS.dt = KS.cfl/2*(KS.dx*KS.dy)/(KS.dx+KS.dy)

        KS.nt = ceil(KS.tEnd/KS.dt);

       # KS.rangex = collect(1:KS.nx)+2 # two ghost cells left and right
       # KS.rangey = collect(1:KS.ny)+2 # same
        for j=KS.rangex
            for i=KS.rangey
                y = KS.y0+KS.dy/2 + (i-3)*KS.dy
                x = KS.x0+KS.dx/2 + (j-3)*KS.dx
                cx = ceil(x);
                cy = ceil(y);
                #cx = x;
                #cy = y;
                KS.SigmaS[i,j] = 1# 1.81-0.09056
                KS.SigmaT[i,j] = 1.5#1.81
                r = 1.0;
                dr = 0.2; # width of the circle
                if (x*x+y*y)>(r-dr) && (x*x+y*y)<(r+dr)
                        for q=KS.rangequad
                                  KS.Source[q,i,j] = 1.0;#KS.nquadpoints;#*KS.dx*KS.dy;#*4.0*pi/KS.nquadpoints;
                        end
                    end
            end
        end

    elseif KS.ictype=="spiral"
        KS.x0 = -4.0;
        KS.x1 = 4.0;
        KS.y0 = -4.0;
        KS.y1 = 4.0;
        KS.dx = (KS.x1-KS.x0)/(KS.nx-4)
        KS.dy = (KS.y1-KS.y0)/(KS.ny-4)
        KS.SigmaS = zeros(KS.ny,KS.nx);
        KS.SigmaT = zeros(KS.ny,KS.nx);
        phi[:] = 1e-5;
        KS.tEnd = 4
        KS.dt = KS.cfl/2*(KS.dx*KS.dy)/(KS.dx+KS.dy)

        KS.nt = ceil(KS.tEnd/KS.dt);

       # KS.rangex = collect(1:KS.nx)+2 # two ghost cells left and right
       # KS.rangey = collect(1:KS.ny)+2 # same

        for j=KS.rangex
            for i=KS.rangey
                y = KS.y0+KS.dy/2 + (i-3)*KS.dy
                x = KS.x0+KS.dx/2 + (j-3)*KS.dx
                cx = ceil(x);
                cy = ceil(y);
                #cx = x;
                #cy = y;.
                KS.SigmaS[i,j] = 1
                KS.SigmaT[i,j] = 0

               # if sqrt(x*x+y*y)>sqrt(theta)-dphi && sqrt(x*x+y*y)<sqrt(theta)+dphi;
               #     for q=KS.rangequad
               #       KS.Source[q,i,j] = 1.0;#KS.nquadpoints;#*KS.dx*KS.dy;#*4.0*pi/KS.nquadpoints;
               #     end
               # end

            end
        end

        for theta = 0:0.01:3*pi
            for r = theta/3:0.01:theta/3+0.1;
            for j=KS.rangex
                for i=KS.rangey
                    y = KS.y0+KS.dy/2 + (i-3)*KS.dy
                    x = KS.x0+KS.dx/2 + (j-3)*KS.dx

                    xx = cos(theta)*r/2
                    yy = sin(theta)*r/2

                    inside = (x<=xx && xx<=x+KS.dx) && (y<=yy && yy<=y+KS.dy)
                    if inside
                        for q=KS.rangequad
                          KS.Source[q,i,j] = 1.0;#KS.nquadpoints;#*KS.dx*KS.dy;#*4.0*pi/KS.nquadpoints;
                        end
                    end
                end
            end
            end
        end

    end

    nicewrite("SigmaS.txt",KS.SigmaS)
    nicewrite("SigmaT.txt",KS.SigmaT)

end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
using Wand
using Base.Test

using PyCall
@pyimport wand.color as color
Color = color.Color

@test Color("white") == Color("#ffffff")


text:
include("../functions/Functions.jl")
include("../evaluate/Evaluate.jl")
include("../opt/SAG.jl")
include("../data/data.jl")

X, y = getFullData()

# init initial point x_0
x0 = zeros(size(X,2),1)

lambda = 5

iter = 250
step =0.01

f = f_logreg(X,y,x0,lambda)
g = sub_g_logreg()

#f = f_svm(X,y,x0,lambda)
#g = sub_g_svm()

@time (theta, strings, iter) = SAG(X,y,x0,f,g,1e-8,lambda,iter,step,size(X,1))

grad = g_logreg(X,y,theta,lambda)(theta)
grad_norm = norm(grad, Inf)
@printf("Inf-Norm of gradient is %f.\n",grad_norm)
println(evaluate(X,y,theta,"logreg",false))

# Lippschitz = largest eigenvalue of Hessian*Hessian'

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
function precompile_adfem()
    PWD = pwd()
    DEPS = joinpath(@__DIR__, "..", "deps")
    change_directory(DEPS)
    ADCME.precompile()

    install_adept()
    install_mfem()
    install_had()

    change_directory(joinpath(DEPS, "build"))
    require_file("CMakeCache.txt") do
        ADCME.cmake()
    end
    require_library("adfem") do
        ADCME.make()
    end

    change_directory(joinpath(DEPS, "MFEM", "build"))
    require_file("CMakeCache.txt") do
        ADCME.cmake()
    end
    require_library("admfem") do
        ADCME.make()
    end

    change_directory(joinpath(DEPS, "MFEM3", "build"))
    require_file("CMakeCache.txt") do
        ADCME.cmake()
    end
    require_library("admfem") do
        ADCME.make()
    end

    PIP = get_pip()
    try
        run(`$(PIP) install pyvista`)
    catch
        @warn "pyvista installation was not successful. 3D plots functionalities are disabled."
    end
    change_directory(PWD)

end


text:
import CSV
import DataFrames; const df = DataFrames
import Plots; const pl = Plots
import ColorSchemes; const cs = ColorSchemes
pl.plotly()

path = "/home/tobias/doc/students/Bachelor_Michelle-Tavernaro/Results/27.11.18/winkel03/"
angles = 135:5:225
data = df.DataFrame[]
for (i, angle) in enumerate(angles)
    push!(data, CSV.read(path*string(angle)*".csv"))
end

oben_whk = map(x -> x[5,4], data)
mitte_whk = map(x -> x[9,4], data)
unten_whk = map(x -> x[13,4], data)

p = pl.plot([oben_whk])
pl.gui(p)


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

using Markdown
using InteractiveUtils

# ╔═╡ dd7b857f-4e75-4bc8-b904-f94f9b220959
using Test

# ╔═╡ d5cea183-ac32-4c5a-b1df-0cebe4421536
using PlutoUI

# ╔═╡ 4104c8bc-4e8b-4988-88c3-05699b9829fb
md"FinancialReportModel"

# ╔═╡ 4b6fef6e-f6ab-4d17-b801-82785324bc75
md"# What is this for? "

# ╔═╡ d6827019-3284-471c-983f-124b1dcf4bb5
md"""
## Financial modelling
There is a broad segment of financial modelling that tries to calculate future state of a company under different circumstances. It is known as a *financial model* or *a 3-way financial model* or *a 3 statement financial model*. It can be used for business valuation, asessment of business robustness, etc. [More] (https://corporatefinanceinstitute.com/resources/knowledge/modeling/types-of-financial-models/)

Conceptually the model is quite straitforward: given a number of assumptions and a starting balance sheet, one applies accounting rules (reporting, more precisely) and calculates future state of accounts.

The model has descrete time (years, quarters), it is deterministic, dependences are known beforehand. So the only goal is to simulate.
"""

# ╔═╡ 1ce2d04c-6c93-42b0-9fe2-4ddccd891674
md"""

## Pluto vs Excel
Usually this type of calculations is done in Excel, due to its ubuquity, and responsiveness. It is quite easy to make what-if analysis. Pluto, on the other hand, can be a substitution for this interface.

However, Excel is akin to a list of paper, one can draw there anything. It is very error prone. The modeller has to be very disciplined in order to keep trust to the calculation results.

There are modelling standards that try to reduce the modelling risk by imposing restrictions on the modeller. One have to adhere to the special worksheet structure, include checks for model integrity, include checks for business logic (e.g. negative cash), etc.

One rule reads as \"one row -- one formula\". In this case the number of unique formulas reduces substantially (degrees of magnitude) and the modeller (or an auditor) can easier keep track of what is behind those walls of figures.

Basicly these formulas represent relationships betweeen accounts and one can model only those. It is quite similar to ODEs. This module tries to do exacly this.

It is interesting that so far the module is Pluto agnostic.

"""

# ╔═╡ 29468f52-5366-4c16-a667-fc3df7b25b3c
md"""
## Why not the ...
### ... DataFrames
I suspect that it will be tricky to implement reactivity for objects inside a df. Another point is the less dependencies, the better.

### ... Modeling Toolkit & DifferentialEquations
I don't know, too ignorant for this.
On one hand, the model looks like a DE with discrete time.
On another hand, DE might be an overkill.
"""

# ╔═╡ 3a47d90f-7e51-42d3-81b4-3ff7d86bf154
md"
# TODO
* build own DAG
* calculate based on DAG <- can be solved by reactivity
* check if DAG is acyclic
* adding variable together with its initial value?
* function for restrictions?
* If Accounts have the same units, they can be added.
* reactivity. Via Observables (?)
* make errors readable. It is now hard to understand what and where goes wrong
* error if DAG refers to a missing
* error if parameters used in formula with wrong syntax: scalar in place for vecrot and vice versa
* Base.show() for model, account
* prettytables? intercept if for Pluto (mime::MIME\"text/html\")"

# ╔═╡ b2802ff6-de7b-4b3b-883c-79d53896da0b
md"# Module"

# ╔═╡ 5adf2e27-7437-42a9-84cb-1375f06862ae
begin
	mutable struct Account
		acc::Array{Union{Number, Missing}, 1}
		unit::Symbol
	end

	Account(len::Number, unit::Symbol) = Account(fill(missing, len), unit)
end

# ╔═╡ 1a04ad02-41ad-4675-a887-ec8e27928987
begin
	struct Finmodel
		n_periods::Integer
		variables::Dict{Symbol, Account}
		parameters
		rules
	end

	Finmodel(n_periods::Integer) = Finmodel(n_periods, Dict(), Dict(), [])
end

# ╔═╡ e24bafaf-b7f8-470b-8dbb-c21f25642800
Base.length(acc::Account) = length(acc.acc)

# ╔═╡ 0607acd3-e09c-4440-90de-931653d92769
function add_variables!(model::Finmodel, variables)
	for v in variables
		# :unit symbol is a temp stab
		model.variables[v] = Account(model.n_periods+1, :unit)
	end
end

# ╔═╡ 12c0a5c6-1eb7-4ec9-a707-cff8c058f2e8
function set_initials!(model::Finmodel, init)
	for (var, val) in init
		model[var][1] = val
	end
end

# ╔═╡ 8cebec2d-586e-4835-bf2c-dc0715cbdfc4
# The name of model variable is unknown at runtime.
# It depends upon the user. This macro passes model name as a string.

macro rules(model, args)
	modelname = string(model)
	return :(set_rules!($(esc(model)), $modelname, $(esc(args))))
end

# ╔═╡ 5097e2d3-7cb2-4de3-a7de-a63bd5d2c69d
# overwrites existing parameters
function set_parameters!(model::Finmodel, params)
	while length(model.parameters)>0
		pop!(model.parameters)
	end

	for (param, val) in params
		@assert length(val)==1 || length(val)==model.n_periods+1
		model.parameters[param] = val
	end
end

# ╔═╡ e8e659ee-5cf6-42ef-b8b2-13a863ffe02b
function parse_formula_on_register(modelname, formula)
	# :var[smth] reference
	replacement = SubstitutionString(modelname * s"[\1][_$time_\2]")
	formula = replace(formula, r"(:\w+)\[([\+-]?\d+)]" => replacement)
	# :parameter reference
	replacement = SubstitutionString(modelname * s"[\1]")
	replace(formula, r"(?<!\[)(:\w+)(?!\])" => replacement)
end

# ╔═╡ c3045fa3-4035-43b1-bff9-1c96dab85326
# overwrites existing rules
function set_rules!(model::Finmodel, modelname, rules)
	while length(model.rules)>0
		pop!(model.rules)
	end

	for (var, rule) in rules
			push!(model.rules,(var, parse_formula_on_register(modelname, rule)))
	end
end

# ╔═╡ febaa730-6616-400a-b451-00d01732542f
function parse_formula_on_calc(formula, t)
	replace(formula, r"\_\$time\_" => "$t")
end

# ╔═╡ 0038633e-3fda-4743-98cd-bb861ca9ebb3
# common [] interface for both vars and parameters
function Base.getindex(m::Finmodel, name)
	if haskey(m.variables, name)
		return m.variables[name]
	elseif haskey(m.parameters, name)
		return m.parameters[name]
	else
		throw(KeyError(name))
	end
end

# ╔═╡ fe7f64de-a0b9-42fd-8cd4-bc4d6f5eedc6
function Base.getindex(m::Finmodel, name, idx::Number)
	if haskey(m.variables, name)
		return m.variables[name][idx]
	elseif haskey(m.parameters, name)
		return m.parameters[name][idx]
	else
		throw(KeyError(name))
	end
end

# ╔═╡ 4a94087d-1953-4409-ac97-44738bc21afc
function Base.getindex(acc::Account, idx)
	return acc.acc[idx]
end

# ╔═╡ acb16d3c-c3f6-4241-994f-1b48967fc816
function Base.setindex!(acc::Account, rhs, idx)
	return acc.acc[idx] =  rhs
end

# ╔═╡ 594c2cbe-0929-43dc-b736-36b6d21d22f6
# cycles over timesteps and calculates var value
# TODO there is a problem if the rules are recorded out of order.
# E.g. for rules "b = a * 2, aₜ=aₜ-1 + 5". calculate! will try to compute b first
# and will result in an error.

function calculate!(m::Finmodel)
	for t in 2:m.n_periods+1
		for (var, rule) in m.rules
			formula = parse_formula_on_calc(rule, t)
			# TODO: somehow catch refs to missing in order to produce
			# meaningful error text
			try
				m.variables[var][t] = eval(Meta.parse(formula))
			finally
				# @show formula
			end
		end
	end
end

# ╔═╡ a31bed2b-f6f4-4777-83c7-cd1d9572a14f

# ╔═╡ ff014c1f-fbff-4db3-aa86-62668af8b71d
md"# Tests"

# ╔═╡ 87ec27f9-a05b-4cfa-95db-7be344f0f868
@testset "Integration test" begin

	# TODO get rid of global. eval() works only for global context
	n = 40
    global model1 = Finmodel(n)
	@test model1.n_periods==n

	@testset "Add variables" begin
		add_variables!(model1,[:a, :b])
		@test isnothing(model1[:a]) == false
		@test isnothing(model1[:b]) == false
		@test length(model1[:a]) == n + 1
	end

	a_init = 10
	b_init = 30
	@testset "Set initials" begin
		set_initials!(model1, [:a => a_init, :b => b_init])
		@test model1[:a][1] == a_init
		@test model1[:b][1] == b_init
		@test ismissing(model1[:a][2])
	end

	@testset "Set parameter" begin
		set_parameters!(model1, [:c => 5])
		@test model1[:c] == 5
	end

	@testset "Formula parsing" begin
		@test eval(Meta.parse("model1[:a][1]*4")) == a_init * 4
		test_formula1 = parse_formula_on_register("model1", ":a[-1]+:b[-1]")
		parsed_f = parse_formula_on_calc(test_formula1, 2)
		@test eval(Meta.parse(parsed_f)) == a_init * 4
	end

	@testset "Setting rules" begin
		set_rules!(model1, "model1", [:a => ":a[-1]+1", :b => ":a[+0]*3"])
		@test !isnothing(model1.rules[1])
		@test length(model1.rules) == 2
	end

	@testset "Calculations" begin
		calculate!(model1)
		@test model1[:a][41] == a_init + n * 1
		@test model1[:b][41] == (a_init + n * 1) * 3
	end
end

# ╔═╡ 9b837bf5-572f-4f8d-84c7-e958f951fc75
model1

# ╔═╡ 62a698bb-dbec-42b9-aa29-e940f79d7611
@testset "Parameters" begin

	# TODO get rid of global. eval() works only for global context
	n2 = 40
    global model2 = Finmodel(n2)
	add_variables!(model2,[:a, :b])

	@testset "Vector parameter" begin
		set_parameters!(model2, [:d => 8, :e => repeat([6], 41)])
		@test length(model2[:e]) == n2 + 1
	end

	# test for using parameter
	# scalar prm instead of vector
	# vector prm instead of scalar
	# should a scalar prm broadcasted to a vector prm?
	# vector prm of length that deviates from the number of periods
end

# ╔═╡ 35442681-d3b4-48c4-b0da-d83f6dfdca93
md"# Elaborated example (WIP)"

# ╔═╡ 49344309-23b0-4edb-b572-3271d1f5bccd
fm = Finmodel(12)

# ╔═╡ 39c1c3f0-4329-440f-be63-4a753bfea510
set_parameters!(fm, [
		:production => fill(10, 13),
		:price => 8,
		:cost => 4,
		:depreciation => 2,
		:incometaxrate => 0.2
		])

# ╔═╡ 9ab3aeb4-3db1-4655-a984-fd5c21fdae48
add_variables!(fm, [
		:revenue, :cogs, :ebitda, :operatingp, :incometax, :netp,
		:cash, :fixedassets, :retained])

# ╔═╡ 3b0a6701-7102-4798-b364-b3989d6f2762
set_initials!(fm, [
		:cash => 0,
		:fixedassets => 100,
		:retained => 100])

# ╔═╡ 8c389a79-5330-4d5a-a21e-acb587b9ab65
@rules fm [
		:revenue => ":production[+0] * :price",
		:cogs => ":production[+0] * :cost",
		:operatingp => ":revenue[+0] - :cogs[+0] - :depreciation",
		:incometax => ":operatingp[+0] * :incometaxrate",
		:netp => ":operatingp[+0] - :incometax[+0]",

		:cash => ":cash[-1] + :netp[+0]",
		:retained => ":retained[-1] + :netp[+0]",
		:fixedassets => ":fixedassets[-1]"
		]

# ╔═╡ cd7ddb24-de63-4c9f-85c9-aa6f7720f452
fm

# ╔═╡ bfb866fd-58a1-44d5-b5eb-8856a7a30566
calculate!(fm)

# ╔═╡ 8bb1e62a-541b-48d7-9b2e-de196f6fd080
fm

# ╔═╡ 64ed55b4-adc0-47c9-8b73-490d8af9cb1f
fm[:production]

# ╔═╡ e2d4aaef-c74a-4733-a0cb-8161b8f1159c
fm[:retained].acc .≈ fm[:fixedassets].acc + fm[:cash].acc

# ╔═╡ 705e87b8-5c65-43c1-9de5-413c6b967c61
PlutoUI.TableOfContents(aside = true)

# ╔═╡ Cell order:
# ╟─4104c8bc-4e8b-4988-88c3-05699b9829fb
# ╟─4b6fef6e-f6ab-4d17-b801-82785324bc75
# ╟─d6827019-3284-471c-983f-124b1dcf4bb5
# ╟─1ce2d04c-6c93-42b0-9fe2-4ddccd891674
# ╟─29468f52-5366-4c16-a667-fc3df7b25b3c
# ╟─3a47d90f-7e51-42d3-81b4-3ff7d86bf154
# ╟─b2802ff6-de7b-4b3b-883c-79d53896da0b
# ╠═1a04ad02-41ad-4675-a887-ec8e27928987
# ╠═5adf2e27-7437-42a9-84cb-1375f06862ae
# ╠═e24bafaf-b7f8-470b-8dbb-c21f25642800
# ╠═0607acd3-e09c-4440-90de-931653d92769
# ╠═12c0a5c6-1eb7-4ec9-a707-cff8c058f2e8
# ╠═8cebec2d-586e-4835-bf2c-dc0715cbdfc4
# ╠═c3045fa3-4035-43b1-bff9-1c96dab85326
# ╠═5097e2d3-7cb2-4de3-a7de-a63bd5d2c69d
# ╠═e8e659ee-5cf6-42ef-b8b2-13a863ffe02b
# ╠═febaa730-6616-400a-b451-00d01732542f
# ╠═0038633e-3fda-4743-98cd-bb861ca9ebb3
# ╠═fe7f64de-a0b9-42fd-8cd4-bc4d6f5eedc6
# ╠═4a94087d-1953-4409-ac97-44738bc21afc
# ╠═acb16d3c-c3f6-4241-994f-1b48967fc816
# ╠═594c2cbe-0929-43dc-b736-36b6d21d22f6
# ╠═a31bed2b-f6f4-4777-83c7-cd1d9572a14f
# ╟─ff014c1f-fbff-4db3-aa86-62668af8b71d
# ╠═dd7b857f-4e75-4bc8-b904-f94f9b220959
# ╠═87ec27f9-a05b-4cfa-95db-7be344f0f868
# ╠═9b837bf5-572f-4f8d-84c7-e958f951fc75
# ╠═62a698bb-dbec-42b9-aa29-e940f79d7611
# ╟─35442681-d3b4-48c4-b0da-d83f6dfdca93
# ╠═49344309-23b0-4edb-b572-3271d1f5bccd
# ╠═39c1c3f0-4329-440f-be63-4a753bfea510
# ╠═9ab3aeb4-3db1-4655-a984-fd5c21fdae48
# ╠═3b0a6701-7102-4798-b364-b3989d6f2762
# ╠═8c389a79-5330-4d5a-a21e-acb587b9ab65
# ╠═cd7ddb24-de63-4c9f-85c9-aa6f7720f452
# ╠═bfb866fd-58a1-44d5-b5eb-8856a7a30566
# ╠═8bb1e62a-541b-48d7-9b2e-de196f6fd080
# ╠═64ed55b4-adc0-47c9-8b73-490d8af9cb1f
# ╠═e2d4aaef-c74a-4733-a0cb-8161b8f1159c
# ╟─d5cea183-ac32-4c5a-b1df-0cebe4421536
# ╟─705e87b8-5c65-43c1-9de5-413c6b967c61

max_stars_repo_name: elgator/FinancialReportModel

max_issues_repo_name: hdlopeza/FinancialReportModel

max_forks_repo_name: hdlopeza/FinancialReportModel


text:
import Clang
import HiGHS_jll

const LIBHIGHS_HEADERS = [
    joinpath(HiGHS_jll.artifact_dir, "include", "interfaces", "highs_c_api.h"),
    # joinpath(HiGHS_jll.artifact_dir, "include", "lp_data", "HConst.h"),
    joinpath(HiGHS_jll.artifact_dir, "include", "util", "HighsInt.h"),
]

const LIBHIGHS_INCLUDES =
    [joinpath(HiGHS_jll.artifact_dir, "include", "util", "HighsInt.h")]

const GEN_DIR = joinpath(@__DIR__, "..", "src", "gen")
if !isdir(GEN_DIR)
    mkdir(GEN_DIR)
end

clang_args = ["-I", joinpath(HiGHS_jll.artifact_dir, "include")]
for header in Clang.find_std_headers()
    push!(clang_args, "-I" * header)
end

wc = Clang.init(
    headers = LIBHIGHS_HEADERS,
    output_file = joinpath(GEN_DIR, "libhighs_api.jl"),
    common_file = joinpath(GEN_DIR, "libhighs_common.jl"),
    clang_includes = LIBHIGHS_INCLUDES,
    clang_args = clang_args,
    header_wrapped = (root, current) -> root == current,
    header_library = x -> "libhighs",
    clang_diagnostics = true,
)

run(wc)

function manual_corrections()
    format_off = "#! format: off\n\n"
    for file in ["libhighs_api.jl"]
        filename = joinpath(GEN_DIR, file)
        content = read(filename, String)
        write(filename, format_off * content)
    end
end

manual_corrections()

rm(joinpath(GEN_DIR, "LibTemplate.jl"))
rm(joinpath(GEN_DIR, "ctypes.jl"))


text:
using TextMining
using FactCheck
using Compat

reload("feature_vector_tests.jl")
reload("cluster_tests.jl")
reload("data_set_tests.jl")
reload("clustering_tests.jl")
reload("distribution_tests.jl")
exitstatus()

max_stars_repo_name: drawkward/SLU-TextMining

max_issues_repo_name: drawkward/SLU-TextMining

max_forks_repo_name: drawkward/SLU-TextMining


text:
win = GtkWindow("A new window")
g = GtkGrid()
a = GtkEntry()  # a widget for entering text
set_gtk_property!(a, :text, "This is Gtk!")
b = GtkCheckButton("Check me!")
c = GtkScale(false, 0:10)     # a slider

# Now let's place these graphical elements into the Grid:
g[1,1] = a    # Cartesian coordinates, g[x,y]
g[2,1] = b
g[1:2,2] = c  # spans both columns
set_gtk_property!(g, :column_homogeneous, true)
set_gtk_property!(g, :column_spacing, 15)  # introduce a 15-pixel gap between columns
push!(win, g)
showall(win)

max_stars_repo_name: ashwani-rathee/Chatbots

max_issues_repo_name: ashwani-rathee/Chatbots

max_forks_repo_name: ashwani-rathee/Chatbots


text:
function solve(case, D, X)
    earliest = D
    for i in Iterators.reverse(X)
        earliest = div(earliest, i) * i
    end

    println("Case #$(case): $(earliest)")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        _, D = parse.(Int, split(readline(stdin)))
        X = parse.(Int, split(readline(stdin)))
        solve(case, D, X)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive