text:
import Base: rand,
             ∈

export AbstractHPolygon,
       an_element,
       isredundant,
       remove_redundant_constraints!,
       addconstraint!,
       vertices_list,
       constraints_list,
       isbounded

# This constant marks the threshold for the number of constraints of a polygon
# above which we use a binary search to find the relevant constraint in a
# support vector query.
#
# NOTE: The value must be strictly greater than 2.
const BINARY_SEARCH_THRESHOLD = 10

"""
    AbstractHPolygon{N<:Real} <: AbstractPolygon{N}

Abstract type for polygons in H-representation (i.e., constraints).

### Notes

All subtypes must satisfy the invariant that constraints are sorted
counter-clockwise.

Every concrete `AbstractHPolygon` must have the following fields:
- `constraints::Vector{LinearConstraint{N, AbstractVector{N}}}` -- the constraints

New subtypes should be added to the `convert` method in order to be convertible.

```jldoctest; setup = :(using LazySets: subtypes)
julia> subtypes(AbstractHPolygon)
2-element Array{Any,1}:
 HPolygon
 HPolygonOpt
```
"""
abstract type AbstractHPolygon{N<:Real} <: AbstractPolygon{N} end

isconvextype(::Type{<:AbstractHPolygon}) = true

# --- AbstractPolygon interface functions ---

"""
    tovrep(P::AbstractHPolygon{N}) where {N<:Real}

Build a vertex representation of the given polygon.

### Input

- `P` -- polygon in constraint representation

### Output

The same polygon but in vertex representation, a `VPolygon`.
"""
function tovrep(P::AbstractHPolygon{N}) where {N<:Real}
    return VPolygon(vertices_list(P; apply_convex_hull=false))
end

"""
    tohrep(P::HPOLYGON) where {HPOLYGON<:AbstractHPolygon}

Build a contraint representation of the given polygon.

### Input

- `P` -- polygon in constraint representation

### Output

The identity, i.e., the same polygon instance.
"""
function tohrep(P::HPOLYGON) where {HPOLYGON<:AbstractHPolygon}
    return P
end

# --- AbstractPolytope interface functions ---

"""
    vertices_list(P::AbstractHPolygon{N};
                  apply_convex_hull::Bool=true,
                  check_feasibility::Bool=true) where {N<:Real}

Return the list of vertices of a polygon in constraint representation.

### Input

- `P`                 -- polygon in constraint representation
- `apply_convex_hull` -- (optional, default: `true`) flag to post-process the
                         intersection of constraints with a convex hull
- `check_feasibility` -- (optional, default: `true`) flag to check whether the
                         polygon was empty (required for correctness in case of
                         empty polygons)

### Output

List of vertices.

### Algorithm

We compute each vertex as the intersection of consecutive lines defined by the
half-spaces.
If `check_feasibility` is active, we then check if the constraints of the
polygon were actually feasible (i.e., they pointed in the *right* direction).
For this we compute the *average* of all vertices and check membership in each
constraint.
"""
function vertices_list(P::AbstractHPolygon{N};
                       apply_convex_hull::Bool=true,
                       check_feasibility::Bool=true) where {N<:Real}
    n = length(P.constraints)
    points = Vector{Vector{N}}(undef, n)
    if n == 0
        return points
    end
    @inbounds for i in 1:n-1
        cap = intersection(Line(P.constraints[i]), Line(P.constraints[i+1]))
        if cap isa EmptySet
            return Vector{Vector{N}}()
        else
            points[i] = element(cap)
        end
    end
    cap = intersection(Line(P.constraints[n]), Line(P.constraints[1]))
    if cap isa EmptySet
        return Vector{Vector{N}}()
    else
        points[n] = element(cap)
    end

    # check if polygon was empty
    if check_feasibility
        avg = sum(points) / length(points)
        if avg ∉ P
            return Vector{Vector{N}}(undef, 0)
        end
    end

    return apply_convex_hull ? convex_hull(points) : points
end

"""
    constraints_list(P::AbstractHPolygon{N}) where {N<:Real}

Return the list of constraints defining a polygon in H-representation.

### Input

- `P` -- polygon in H-representation

### Output

The list of constraints of the polygon.
The implementation guarantees that the constraints are sorted counter-clockwise.
"""
function constraints_list(P::AbstractHPolygon{N}) where {N<:Real}
    return P.constraints
end

# --- LazySet interface functions ---

"""
    an_element(P::AbstractHPolygon{N}) where {N<:Real}

Return some element of a polygon in constraint representation.

### Input

- `P` -- polygon in constraint representation

### Output

A vertex of the polygon in constraint representation (the first one in the order
of the constraints).
"""
function an_element(P::AbstractHPolygon{N}) where {N<:Real}
    @assert length(P.constraints) >= 2 "polygon has less than two constraints"
    return element(intersection(Line(P.constraints[1]),
                                Line(P.constraints[2])))
end

"""
    ∈(x::AbstractVector{N}, P::AbstractHPolygon{N}) where {N<:Real}

Check whether a given 2D point is contained in a polygon in constraint
representation.

### Input

- `x` -- two-dimensional point/vector
- `P` -- polygon in constraint representation

### Output

`true` iff ``x ∈ P``.

### Algorithm

This implementation checks if the point lies on the outside of each edge.
"""
function ∈(x::AbstractVector{N}, P::AbstractHPolygon{N}) where {N<:Real}
    @assert length(x) == 2

    for c in P.constraints
        if !_leq(dot(c.a, x), c.b)
            return false
        end
    end
    return true
end

"""
    rand(::Type{HPOLYGON}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing,
         [num_constraints]::Int=-1) where {HPOLYGON<:AbstractHPolygon}

Create a random polygon in constraint representation.

### Input

- `HPOLYGON`        -- type for dispatch
- `N`               -- (optional, default: `Float64`) numeric type
- `dim`             -- (optional, default: 2) dimension
- `rng`             -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`            -- (optional, default: `nothing`) seed for reseeding
- `num_constraints` -- (optional, default: `-1`) number of constraints of the
                       polygon (must be 3 or bigger; see comment below)

### Output

A random polygon in constraint representation.

### Algorithm

We create a random polygon in vertex representation and convert it to constraint
representation.
See [`rand(::Type{VPolygon})`](@ref).
For non-flat polygons the number of vertices and the number of constraints are
identical.
"""
function rand(::Type{HPOLYGON};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing,
              num_constraints::Int=-1) where {HPOLYGON<:AbstractHPolygon}
    @assert dim == 2 "cannot create a random $HPOLYGON of dimension $dim"
    @assert num_constraints < 0 || num_constraints >= 3 "cannot construct a " *
        "random $HPOLYGON with only $num_constraints constraints"
    rng = reseed(rng, seed)
    vpolygon = rand(VPolygon; N=N, dim=dim, rng=rng, seed=seed,
                    num_vertices=num_constraints)
    return convert(HPOLYGON, vpolygon)
end

# --- common AbstractHPolygon functions ---

"""
    isredundant(cmid::LinearConstraint{N}, cright::LinearConstraint{N},
                cleft::LinearConstraint{N}) where {N<:Real}

Check whether a linear constraint is redundant wrt. two surrounding constraints.

### Input

- `cmid`   -- linear constraint of concern
- `cright` -- linear constraint to the right (clockwise turn)
- `cleft`  -- linear constraint to the left (counter-clockwise turn)

### Output

`true` iff the constraint is redundant.

### Algorithm

We first check whether the angle between the surrounding constraints is < 180°,
which is a necessary condition (unless the direction is identical to one of the
other two constraints).
If so, we next check if the angle is 0°, in which case the constraint `cmid` is
redundant unless it is strictly tighter than the other two constraints.
If the angle is strictly between 0° and 180°, the constraint `cmid` is redundant
if and only if the vertex defined by the other two constraints lies inside the
set defined by `cmid`.
"""
function isredundant(cmid::LinearConstraint{N},
                     cright::LinearConstraint{N},
                     cleft::LinearConstraint{N}) where {N<:Real}
    samedir_check = false
    # determine angle between surrounding constraints
    if !is_right_turn(cright.a, cleft.a)
        # angle is > 180°
        samedir_check = true
    elseif is_right_turn(cleft.a, cright.a)
        # angle is 0° or 180°
        if samedir(cright.a, cleft.a)[1]
            # angle is 0°
            if samedir(cright.a, cmid.a)[1] && samedir(cleft.a, cmid.a)[1]
                # all three constraints have the same direction
                # constraint is redundant unless it is tighter than the others
                return !is_tighter_same_dir_2D(cmid, cright, strict=true) &&
                       !is_tighter_same_dir_2D(cmid, cleft, strict=true)
            else
                # surrounding constraints have the same direction but the
                # central constraint does not -> corner case with just three
                # unbounding constraints (usually occurs during incremental
                # addition of constraints)
                return false
            end
        else
            # angle is 180°
            samedir_check = true
        end
    end
    # check if the constraint has the same direction as one of the two
    if samedir(cright.a, cmid.a)[1]
        return !is_tighter_same_dir_2D(cmid, cright, strict=true)
    elseif samedir(cleft.a, cmid.a)[1]
        return !is_tighter_same_dir_2D(cmid, cleft, strict=true)
    elseif samedir_check
        # not the same direction => constraint is not redundant
        return false
    end
    cap = intersection(Line(cright), Line(cleft))
    @assert cap isa Singleton
    return cap ⊆ cmid
end

"""
    remove_redundant_constraints!(P::AbstractHPolygon)

Remove all redundant constraints of a polygon in constraint representation.

### Input

- `P` -- polygon in constraint representation

### Output

The same polygon with all redundant constraints removed.

### Notes

Since we only consider bounded polygons and a polygon needs at least three
constraints to be bounded, we stop removing redundant constraints if there are
three or less constraints left.
This means that for non-bounded polygons the result may be unexpected.

### Algorithm

We go through all consecutive triples of constraints and check if the one in the
middle is redundant.
For this we assume that the constraints are sorted.
"""
function remove_redundant_constraints!(P::AbstractHPolygon)
    C = P.constraints
    i = 1
    go_on = true
    cright = C[length(C)] # define initial cright here
    while length(C) >= 3 && go_on
        cmid = C[i]
        if i < length(C)
            cleft = C[i+1]
        elseif i == length(C)
            cleft = C[1]
            go_on = false
        end
        if isredundant(cmid, cright, cleft)
            deleteat!(C, i)
        else
            i += 1
            cright = C[i-1] # update cright (updates less often than cmid/cleft)
        end
    end
    return P
end

"""
    addconstraint!(P::AbstractHPolygon{N},
                   constraint::LinearConstraint{N};
                   [linear_search]::Bool=(length(P.constraints) <
                                          BINARY_SEARCH_THRESHOLD),
                   [prune]::Bool=true) where {N<:Real}

Add a linear constraint to a polygon in constraint representation, keeping the
constraints sorted by their normal directions.

### Input

- `P`             -- polygon in constraint representation
- `constraint`    -- linear constraint to add
- `linear_search` -- (optional, default: `length(constraints) <
                     BINARY_SEARCH_THRESHOLD`) flag to choose between linear
                     and binary search
- `prune`         -- (optional, default: `true`) flag for removing redundant
                     constraints in the end
"""
function addconstraint!(P::AbstractHPolygon{N},
                        constraint::LinearConstraint{N};
                        linear_search::Bool=(length(P.constraints) <
                                             BINARY_SEARCH_THRESHOLD),
                        prune::Bool=true) where {N<:Real}
    return addconstraint!(P.constraints, constraint,
                          linear_search=linear_search, prune=prune)
end

"""
    addconstraint!(constraints::Vector{LC},
                   new_constraint::LinearConstraint{N};
                   [linear_search]::Bool=(length(P.constraints) <
                                          BINARY_SEARCH_THRESHOLD),
                   [prune]::Bool=true
                  ) where {N<:Real, LC<:LinearConstraint{N}}

Add a linear constraint to a sorted vector of constrains, keeping the
constraints sorted by their normal directions.

### Input

- `constraints`    -- vector of linear constraintspolygon in constraint
                      representation
- `new_constraint` -- linear constraint to add
- `linear_search`  -- (optional, default: `length(constraints) <
                      BINARY_SEARCH_THRESHOLD`) flag to choose between linear
                      and binary search
- `prune`          -- (optional, default: `true`) flag for removing redundant
                      constraints in the end

### Algorithm

If `prune` is active, we check if the new constraint is redundant.
If the constraint is not redundant, we perform the same check to the left and to
the right until we find the first constraint that is not redundant.
"""
function addconstraint!(constraints::Vector{LC},
                        new_constraint::LinearConstraint{N};
                        linear_search::Bool=(length(constraints) <
                                             BINARY_SEARCH_THRESHOLD),
                        prune::Bool=true
                       ) where {N<:Real, LC<:LinearConstraint{N}}
    m = length(constraints)
    k = m
    if k > 0
        d = new_constraint.a
        if d <= constraints[1].a
            k = 0
        elseif linear_search
            # linear search
            while d <= constraints[k].a
                k -= 1
            end
        else
            # binary search
            k = binary_search_constraints(
                d, constraints, k, 1 + div(k, 2), choose_lower=true)
        end
    end

    # here constraints[k] <= new_constraint <= constraints[(k%m)+1]
    if prune && m >= 2
        # check if new constraint is redundant
        k += 1
        below = k == 1 ? m : k - 1
        above = k == m + 1 ? 1 : k
        if isredundant(new_constraint, constraints[below], constraints[above])
            return nothing
        end
        # insert new constraint
        insert!(constraints, k, new_constraint)
        m += 1
        # check if old constraints below became redundant
        while m > 2
            center = k == 1 ? m : k - 1
            below = center == 1 ? m : center - 1
            if isredundant(constraints[center], constraints[below],
                           new_constraint)
                deleteat!(constraints, center)
                if center < k
                    k -= 1
                end
                m -= 1
            else
                break
            end
        end
        # check if old constraints above became redundant
        while m > 2
            center = k == m ? 1 : k + 1
            above = center == m ? 1 : center + 1
            if isredundant(constraints[center], new_constraint,
                            constraints[above])
                deleteat!(constraints, center)
                if center < k
                    k -= 1
                end
                m -= 1
            else
                break
            end
        end
    else
        insert!(constraints, k+1, new_constraint)
    end
    return nothing
end

"""
    binary_search_constraints(d::AbstractVector{N},
                              constraints::Vector{<:LinearConstraint{N}},
                              n::Int,
                              k::Int;
                              [choose_lower]::Bool=false) where {N}

Performs a binary search in the constraints.

### Input

- `d`            -- direction
- `constraints`  -- constraints
- `n`            -- number of constraints
- `k`            -- start index
- `choose_lower` -- (optional, default: `false`) flag for choosing the lower
                    index (see the 'Output' section)

### Output

In the default setting, the result is the smallest index `k` such that
`d <= constraints[k]`, or `n+1` if no such `k` exists.
If the `choose_lower` flag is set, the result is the largest index `k` such
that `constraints[k] < d`, which is equivalent to being `k-1` in the normal
setting.
"""
function binary_search_constraints(d::AbstractVector{N},
                                   constraints::Vector{<:LinearConstraint{N}},
                                   n::Int,
                                   k::Int;
                                   choose_lower::Bool=false) where {N}
    lower = 1
    upper = n+1
    while lower + 1 < upper
        if constraints[k].a <= d
            lower = k
        else
            upper = k
        end
        k = lower + div(upper - lower, 2)
    end
    if choose_lower
        return lower
    else
        if lower == 1 && !(constraints[1].a <= d)
            # special case for index 1
            return 1
        end
        return upper
    end
end

"""
    isbounded(P::AbstractHPolygon, [use_type_assumption]::Bool=true)

Determine whether a polygon in constraint representation is bounded.

### Input

- `P`                   -- polygon in constraint representation
- `use_type_assumption` -- (optional, default: `true`) flag for ignoring the
                           type assumption that polygons are bounded

### Output

`true` if `use_type_assumption` is activated.
Otherwise, `true` iff `P` is bounded.

### Algorithm

If `!use_type_assumption`, we convert `P` to an `HPolyhedron` `P2` and then use
`isbounded(P2)`.
"""
function isbounded(P::AbstractHPolygon, use_type_assumption::Bool=true)
    if use_type_assumption
        return true
    end
    if isempty(constraints_list(P))
        return false
    end
    return isbounded_unit_dimensions(P)
end


text:
function ∑(coll)
    isempty(coll) ? 0. : sum(coll)
end
function propersubsets(coll)
    (S for S in subsets(coll) if 0 < length(S) < length(coll))
end

ext: jl


text:
function smooth(y, win_len=11, win_method=2)
# Based on http://www.simulkade.com/posts/2015-05-07-how-to-smoothen-noisy-data.html
# This function requires the DSP package to be installed
# 1: flat
# 2: hanning
# 3: hamming ...
if win_len%2==0
  win_len+=1 # only use odd numbers
end
if win_method == 1
  w=ones(win_len)
elseif win_method==2
  w=DSP.hanning(win_len)
elseif win_method==3
  w=DSP.hamming(win_len)
end

if win_len<3
  return y
elseif length(y)<win_len
  return y
else
        y_new = [2*y[1].-reverse(y[1:win_len],1); y[:]; 2*y[end].-reverse(y[end-win_len:end],1)]
  y_smooth = conv(y_new, w/sum(w))
  ind = floor(Int, 1.5*win_len)
  return y_smooth[1+ind:end-ind-1]
end

end # end of function


text:
d2 = readdlm("./tsp/distances_usa.txt")
c2 = readdlm("./tsp/cities_usa.txt")

# Some basic unit tests
# If a Test fails, there is an error in your implementation
# Passing the tests does not mean your implementation is correct!
@testset begin
    _s1 = find_subtours(
        [
            [0.0 1.0];
            [1.0 0.0]
        ]
    )  # no sub-tour should be returned here

    _s2 = find_subtours(
        [
            [0.0 1.0 0.0 0.0];
            [1.0 0.0 0.0 0.0];
            [0.0 0.0 0.0 1.0];
            [0.0 0.0 1.0 0.0]
        ]
    )  # two sub-tours should be returned here

    @test length(_s1) == 0
    @test length(_s2) == 2
end;

sol2 = solve_tsp2(50, d2, GLPK.Optimizer)
RoadTrip(50, c2, sol2)
![](tsp/solution_usa.png)


text:
""" PIV PARAMETER OBJECT """
mutable struct PIVParameters
         corr::Union{ZNCC,FFT}          # cross-correlation algorithm, "zncc" or "fft"
    interArea::Tuple{Int64,Int64,Int64} # interrogation area size(s) in pixels
 searchMargin::Tuple{Int64,Int64,Int64} # search margin(s) in pixels
      overlap::Tuple{Int64,Int64,Int64} # pixel distance(s) between contiguous interrogation areas
        mpass::Int64    # multi-pass depth
        width::Int64    # exclusion radius around 1st peak for finding 2nd peak
         peak::String   # one of "", "gaussian", "centroid"
     sigNoise::String   # one of "", "peak2peak"
      filtFun::Function # Filtering function to skip cross-correlating certain interrogation areas (e.g. background)
    threshold::Float64  # skip cross correlation if filtFun( interrogation area ) < threshold
        mpFun::Function # scaling function for IA, SM and ovp during multipass
       mpReps::Array{Int64,1} # number of times each multipass iteration is repeated. Not implemented.
        units::Tuple{DU,DU,DU,TU} # spatial and temporal units of the analyzed data
         args::Tuple{III,III,III,I,I,S,S,Function,F,Function,Array{I,1}}
end

PP = PIVParameters;

# constructor
PIVParameters( corr::Union{ZNCC,FFT}, IA::UnI, SM::UnI, ovp::UnI, mpass::I, width::I, peak::S,
               sn::S, ffn::Function, th::F, mpfn::Function, mpReps::Array{I,1}, uns::UNTS
             ) = PIVParameters( corr, toIII(IA), toIII(SM), toIII(ovp), mpass, width, peak, sn,
                                ffn, th, mpfn, mpReps, uns, ( toIII(IA), toIII(SM), toIII(ovp),
                                mpass, width, peak, sn, ffn, th, mpfn, mpReps )
                              )

function parseCorr( corr::String )
    lowercase(corr) == "zncc" &&  return ZNCC();
    lowercase(corr) == "fft"  &&  return  FFT();
end

# interrogation area size, search margin and overlap can be specified for each dimension,
# or a sigle value can be applied to all dimensions.
toIII( input::Int64 ) = ( input, input, input )
toIII( input::Tuple{Int64,Int64} ) = ( input[1], input[2], 0 )
toIII( input::Tuple{Int64,Int64, Int64} ) = input

# constructor with default values
function setPIVParameters(; interArea::UnI=32, searchMargin::UnI=0, overlap::UnI=0,
                            mpass::I=1, width::I=2, corr::S="zncc",
                            peak::S="gaussian", sigNoise::S="",
                            filtFun=(x)->maxval(x), threshold::F=-1.0,
                            mpFun=(x,mp)->(x*mp), mpReps::Array{I,1}=[1,],
                            units=(1m,1m,1m,1second) )
    return PIVParameters( parseCorr(corr), interArea, searchMargin, overlap, mpass,
                          width, peak, sigNoise, filtFun, threshold, mpFun, mpReps, units )
end

function checkPIVParameters( dsize, N, params )
    !( N !== 2  ||  N !== 3 ) && error("Your input data isn't 2D nor 3D.")
    IA = params.interArea;
    if N == 2
         IA = params.interArea[1:2]
    end
    any( params.interArea .< 2 ) && error("interArea = 1 makes little sense.")
    any( IA .> dsize ) && error("interArea is bigger than the data")
    any( params.overlap .>= params.interArea ) && error("Overlap must be < IA size")

    if size(params.mpReps,1) !== params.mpass
        updateParameters!( params, :mpReps, ones( Int64, params.mpass ) )
    end
end

""" SYNTHETIC PARTICLES PARAMETER OBJECT """
mutable struct synthParameters
      n::Int64   # total number of particles, ignored when dens > 0
   dens::Int64   # particle density ( nºparticles/IA )
      d::Int64   # depth. If depth > 1, generates a volume
      w::Int64   # width of generated image/volume
      h::Int64   # height of generated image/volume
     i0::Float64 # max intensity of particles
     dt::Float64 # standard deviation of gaussian around each particles
     th::Float64 # standard deviation of gaussian around lightsheet
    err::Float64 # determines rendering radius around each pixel
      z::Int64   # light-sheet position. Default z = 0
  noise::Int64   # varible for controlling noise. Not really implemented.
    rad::Int64   # overwrites the rendering radius computed from err.
   mode::Bool    # switches between ( x, y, z ) or ( y, x, z ) coordinates
   args::Tuple{I,I,I,I,I,F,F,F,F,I,I,I,Bool}
end

SP = synthParameters;

# constructor
synthParameters( n::Int64, dens::Int64, d::Int64, w::Int64, h::Int64, i0::Float64, dt::Float64,
                 th::Float64, err::Float64, z::Int64, noise::Int64, rad::Int64, mode::Bool
               ) = synthParameters( n, dens, d, w, h, i0, dt, th, err, z, noise, rad, mode,
                                   (n, dens, d, w, h, i0, dt, th, err, z, noise, rad, mode)
                                  )

# constructor with default values
function setSynthParameters(; n=0::I, dens=0::I, d=1::I, w=300::I, h=300::I,
                              z=0::I, i0=255.0::F, dt=3.0::F, err=0.1::F,
                              th=10.0::F, noise=0::I, rad=0::I, mode="whd"::S )

    return synthParameters( n, dens, d, w, h, i0, dt, th, err, z, noise, rad, mode=="whd" )
end

""" TRANSFORMATION PARAMETERS """
mutable struct transformParameters{N,M}
      kind::Symbol            # currently ony :translation is supported
     means::NTuple{N,Float64} # mean value of transform in x, y [ and z ] if N = 3
      vars::NTuple{N,Float64} # variance component in covariance matrix
 cov_ratio::NTuple{M,Float64} # covariance components in covariance matrix. -1 < 1
      args::Tuple{Symbol,NTuple{N,Float64},NTuple{N,Float64},NTuple{M,Float64}}
end

# constructor
transformParameters( kind::Symbol, means::NTuple{N,Float64}, vars::NTuple{N,Float64},
                     covs::NTuple{M,Float64}
                   ) where {N,M} = transformParameters{N,M}( kind, means, vars, covs,
                                                            (kind, means, vars, covs)
                                                           )

TP2D = transformParameters{2,1}
TP3D = transformParameters{3,3}
TP   = Union{TP2D,TP3D}

# constructor with default values
function setTransformParameters(; kind=:translation, means=(1.0, 1.0), vars=(eps(Float64),
								  eps(Float64)), cov_ratio=(0.0,) )

    return transformParameters( kind, means, vars, cov_ratio );
end

to2DTransform( t2D::TP2D ) = t2D

function to2DTransform( t3D::TP3D )
    means     = (  t3D.means[1], t3D.means[2] );
    vars      = (  t3D.vars[1] ,  t3D.vars[2] );
    cov_ratio = (  t3D.cov_ratio[1], );

    return transformParameters{2,1}( t3D.kind, means, vars, cov_ratio );
end

function to3DTransform( t2D::TP2D )
    means     = (   t2D.means[1]  ,   t2D.means[2]  , (t2D.means[1] + t2D.means[2])/2 );
    vars      = (    t2D.vars[1]  ,    t2D.vars[2]  , ( t2D.vars[1] +  t2D.vars[2])/2 );
    cov_ratio = ( t2D.cov_ratio[1], t2D.cov_ratio[1], t2D.cov_ratio[1] );

    return transformParameters{3,3}( t2D.kind, means, vars, cov_ratio );
end

to3DTransform( t3D::TP3D ) = t3D

""" EVALUATION PARAMETER OBJECT """
mutable struct metaParameters
 metaloop::Int64    # determines the number of iterations where the changing variable goes from its "min" to "max" values
 variable::Symbol   # specifies what variable will be changed each iteration.
      min::Float64  # minimum value of changing variable
      max::Float64  # maximum value of changing variable
  repeats::Int64    # number of measurements for computing bias and error
     args::Tuple{Int64,Symbol,Float64,Float64,Int64}
end

MP = metaParameters;

metaParameters( metaloop::Int64, variable::Symbol, min::Float64, max::Float64, repeats::Int64
              ) = metaParameters( metaloop, variable, min, max, repeats,
                                 (metaloop, variable, min, max, repeats)
                                )

"""
    Creates a set of parameters for directing PIV evaluation. Accepted parameters are:
        variable (  Symbol ): a symbol indicating what variable should be evaluated. Can be any parameters
                              in synthParameters, PIVParameters or transformParameters.
        metaloop (  Int64  ): number of iterations
            min  ( Float64 ): min value of changing variable
            max  ( Float64 ): max value of changing variable
        repeats  (  Int64  ): number of repetitions
"""
function setMetaParameters( ; metaloop=0::I, var=Symbol("")::Symbol, min=0.0::F,  max=0.0::F, repeats=1::I )

    return metaParameters( metaloop, var, min,  max, repeats )
end

# function to update parameters, and also update the .args tuple, which was created to pass parameters
# by tuple deconstruction "parameter.args..."

function updateParameters!( params, field, newVal )
    hasfield = false;
    fieldidx =     0;
    for e in fieldnames( typeof(params) )
        fieldidx += 1;
        if ( e == field )
            hasfield = true;
            setproperty!( params, e, newVal )
            break
        end
    end

    ( !hasfield ) && return nothing;

    if typeof(params) == PIVParameters
        fieldidx -= 1 # ignore corr parameter
    end
    newargs = [ arg for arg in params.args ]
    newargs[ fieldidx ] = newVal;
    setproperty!( params, :args, Tuple(newargs) )
end


text:
using CUDA

@inline function CheckBlockBounds(offset, block, grids)
        x = offset.x + block.x
        if x < 1 || x > grids
                return true
        end
        x = offset.y + block.y
        if x < 1 || x > grids
                return true
        end
        x = offset.z + block.z
        if x < 1 || x > grids
                return true
        end
        return false
end

@inline function CalcDistance(query, point)
    tempdist = Float32(0)
    for d = 1:length(query)
        @inbounds tempdist += CUDA.pow(query[d] - point[d], 2)
    end
    tempdist = CUDA.sqrt(tempdist)
    return tempdist
end

function cuGridKnnSimpleNoTypesFunction(devPoints, devQueries,
        Pointsperblock, Queriesperblock,
        IntPointsperblock, IntQueriesperblock,
        devDistances, devNeighbours,
        dimensions, grids, offset)

    if CheckBlockBounds(offset, blockIdx(), grids)
        return nothing
    end
    q_bid = ((blockIdx().z - 1)*gridDim().y*gridDim().x
            + (blockIdx().y - 1)*gridDim().x
            + blockIdx().x)
    p_bid = ((blockIdx().z + offset.z- 1)*gridDim().y*gridDim().x
            + (blockIdx().y + offset.y - 1)*gridDim().x
            + blockIdx().x + offset.x)
    tid = threadIdx().x + (threadIdx().y-1)*blockDim().x
    stride = (blockDim().x)*(blockDim().y)

    SharedPoints = @cuDynamicSharedMem(Float32, (dimensions, stride))
    SharedQueries = @cuDynamicSharedMem(Float32, (dimensions, stride), dimensions*stride*sizeof(Float32))

    @inbounds startPoints = IntPointsperblock[p_bid]
    @inbounds startQueries = IntQueriesperblock[q_bid]
    @inbounds totalPoints = Pointsperblock[p_bid]
    @inbounds totalQueries = Queriesperblock[q_bid]
    @inbounds Points = @view devPoints[(startPoints+1):(startPoints+totalPoints), :]
    @inbounds Queries = @view devQueries[(startQueries+1):(startQueries+totalQueries), :]
    @inbounds Distances = @view devDistances[(startQueries+1):(startQueries+totalQueries)]
    @inbounds Neighbours = @view devNeighbours[(startQueries+1):(startQueries+totalQueries)]
    @inbounds query = @view SharedQueries[:, tid]
    dist = Float32(0)
    nb = Int32(0)
    for q in 0:stride:(totalQueries-1)
        if tid + q <= totalQueries
            for d = 1:dimensions
                @inbounds  query[d] = Queries[tid+q, d]
            end
            @inbounds dist = Distances[tid+q]
            @inbounds nb = Neighbours[tid+q]
        end
        for p in 0:stride:(totalPoints-1)

            CUDA.sync_threads()
            if tid + p <= totalPoints
                for d in 1:dimensions
                    @inbounds SharedPoints[d, tid] = Points[tid+p, d]
                end
            end
            CUDA.sync_threads()
            bounds::Int32 = CUDA.min(stride, totalPoints-p)
            for i in 1:bounds
                @inbounds point = @view SharedPoints[:, (i+tid-2)%bounds+1]
                tempdist = CalcDistance(query, point)
                if tempdist < dist
                    dist = tempdist
                    nb = Int32(startPoints + p + (i+tid-2)%bounds+1)
                end
            end
        end
        if tid + q <= totalQueries
            @inbounds Distances[tid+q] = dist
            @inbounds Neighbours[tid+q] = nb
        end
    end
    return nothing
end

function cuda_knn_simple_no_types_function(OrderedPoints, OrderedQueries,
    PointsPerBlock, QueriesPerBlock,
    IntegralPointsPerBlock, IntegralQueriesPerBlock,
    numOfPoints, numOfQueries, numOfGrids, dimensions)

    ## Data transfer
    devPoints = CuArray(OrderedPoints)
    devQueries = CuArray(OrderedQueries)
    devPointsPerBlock = CuArray(PointsPerBlock)
    devQueriesPerBlock = CuArray(QueriesPerBlock)
    devIntegralPointsPerBlock = CuArray(IntegralPointsPerBlock)
    devIntegralQueriesPerBlock = CuArray(IntegralQueriesPerBlock)
    devRes = CUDA.fill(Float32(100), numOfQueries)
    devNeighbours = CUDA.zeros(Int32, numOfQueries)
    ## Config
    # dimensions = Int32(dimensions)
    dimensions = (dimensions)
    grids = (numOfGrids)
    offset = (x=(0),y=(0),z=(0))
    thread_groups = 2
    shmem=2*thread_groups*32*dimensions*sizeof(Float32)
    ## Running for main block.
    @cuda(blocks=(grids,grids,grids), threads=(32,thread_groups), shmem=shmem,
      cuGridKnnSimpleNoTypesFunction(devPoints, devQueries,
        devPointsPerBlock, devQueriesPerBlock,
        devIntegralPointsPerBlock, devIntegralQueriesPerBlock,
        devRes, devNeighbours, dimensions, grids, (x=0,y=0,z=0)))
    #Running for neighbour blocks.
    for x = -1:1
        for y = -1:1
            for z = -1:1
                if x == 0 && y == 0 && z == 0
                    continue
                end
                offset = (x=x,y=y,z=z)
                @cuda(blocks=(grids,grids,grids), threads=(32,thread_groups), shmem=shmem,
                  cuGridKnnSimpleNoTypesFunction(devPoints, devQueries,
                    devPointsPerBlock, devQueriesPerBlock,
                    devIntegralPointsPerBlock, devIntegralQueriesPerBlock,
                    devRes, devNeighbours, dimensions, grids, offset))

            end
        end
    end

    return Array(devNeighbours), Array(devRes)

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type SumpTank
	SumpTank()=begin
		PP=outers.PP
		NComp=outers.NComp
		[
			:(LI = Levelpercent_Initial),
			:(OutletLiquid.T = Temperature_Initial),
			:(OutletLiquid.z(1:NComp - 1) = Composition_Initial(1:NComp - 1)/sum(Composition_Initial)),
		],
		[
			"Initial level Percent","Initial Outlet Liquid Temperature","Initial Outlet Liquid Composition Normalized",
		],
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of components",
				:Lower=>1
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Mol Weight",
				:Hidden=>true
			)),(NComp)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Pi value",
				:Default=>3.141593,
				:Hidden=>true,
				:Symbol=>"\\pi"
			)),
			acceleration(Dict{Symbol,Any}(
				:Brief=>"Gravity Acceleration",
				:Default=>9.81,
				:Hidden=>true
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Valid=>["elliptical","hemispherical","flat"],
				:Default=>"flat"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Vessel diameter",
				:Symbol=>"D_{i}"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Side length of the cylinder shell",
				:Symbol=>"L_{vessel}"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Elliptical Head Total Volume",
				:Hidden=>true,
				:Symbol=>"V_{head}^{elliptical}"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Hemispherical Head Total Volume",
				:Hidden=>true,
				:Symbol=>"V_{head}^{hemispherical}"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Cylinder Total Volume",
				:Hidden=>true,
				:Symbol=>"V_{cylinder}"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Vessel radius",
				:Hidden=>true,
				:Symbol=>"R_{cylinder}"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Initial liquid height in Percent",
				:Default=>0.70
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Initial Liquid Temperature",
				:Default=>330
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Initial Composition",
				:Default=>0.10
			)),(NComp)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Feed Stream",
				:PosX=>0.22,
				:PosY=>0,
				:Symbol=>"_{in}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream",
				:PosX=>0.50,
				:PosY=>1,
				:Symbol=>"_{out}^{Liquid}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream",
				:PosX=>1,
				:PosY=>0.20,
				:Symbol=>"_{out}^{Vapour}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream",
				:PosX=>0.68,
				:PosY=>0,
				:Symbol=>"_{out}^{Vapour}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Heat Duty",
				:Protected=>false,
				:Symbol=>"Q_{in}"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Vessel total volume",
				:Protected=>true,
				:Symbol=>"V_{total}"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Vessel volume content",
				:Protected=>true,
				:Symbol=>"V_{filled}"
			)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar Holdup in the Vessel",
				:Protected=>true
			)),(NComp)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar liquid holdup",
				:Protected=>true
			)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar vapour holdup",
				:Protected=>true
			)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Total Energy Holdup in the Vessel",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Volume",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar volume",
				:Protected=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"liquid height",
				:Protected=>true
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Vessel cylinder shell Cross section area",
				:Hidden=>true,
				:Symbol=>"A_{cross}"
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Equilibrium pressure on the liquid surface",
				:Protected=>true,
				:Symbol=>"P_{eq}"
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Static head at the bottom of the tank",
				:Protected=>true,
				:Symbol=>"P_{static}^{Liquid}"
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Level Indicator",
				:PosX=>1,
				:PosY=>0.73,
				:Protected=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Temperature Indicator",
				:PosX=>1,
				:PosY=>0.40,
				:Protected=>true
			)),
			[
				:(Across = 0.25*(pi*Diameter^2)),
				:(Vtotal = Vhead_elliptical + Vcylinder),
				:(Vfilled = 0.25*pi*(((Diameter*Level)/(0.25*Diameter))^2)*(0.25*Diameter-Level/3)),
				:(Vfilled = 0.25*pi*(Diameter^2)*(Level - 0.25*Diameter/3)),
				:(Vtotal = Vhead_hemispherical + Vcylinder),
				:(Vfilled = 0.25*pi*(Level^2)*(2*Diameter-4*Level/3)),
				:(Vfilled = 0.25*pi*((2/3)*((0.5*Diameter)^3) - (0.25*(Diameter)^3) + Level*Diameter^2)),
				:(Vtotal = Vcylinder),
				:(Vfilled = Across*Level),
				:(diff(TotalHoldup) = InletLiquid.F*InletLiquid.z + InletVapour.F*InletVapour.z- OutletLiquid.F*OutletLiquid.z - OutletVapour.F*OutletVapour.z),
				:(diff(E) = InletLiquid.F*InletLiquid.h + InletVapour.F*InletVapour.h - OutletLiquid.F*OutletLiquid.h - OutletVapour.F*OutletVapour.h + InletQ),
				:(TotalHoldup = LiquidHoldup*OutletLiquid.z + VapourHoldup*OutletVapour.z),
				:(E = LiquidHoldup*OutletLiquid.h + VapourHoldup*OutletVapour.h - OutletLiquid.P*Vtotal),
				:(sum(OutletLiquid.z)=1.0),
				:(sum(OutletLiquid.z)=sum(OutletVapour.z)),
				:(vL = PP.LiquidVolume(OutletLiquid.T, Peq, OutletLiquid.z)),
				:(vV = PP.VapourVolume(OutletVapour.T, Peq, OutletVapour.z)),
				:(PP.LiquidFugacityCoefficient(OutletLiquid.T, Peq, OutletLiquid.z)*OutletLiquid.z = PP.VapourFugacityCoefficient(OutletVapour.T, Peq, OutletVapour.z)*OutletVapour.z),
				:(OutletVapour.T = OutletLiquid.T),
				:(OutletVapour.P = Peq),
				:(Pstatic = PP.LiquidDensity(OutletLiquid.T, Peq, OutletLiquid.z) * g * Level),
				:(OutletLiquid.P = Peq + Pstatic),
				:(Vtotal = LiquidHoldup * vL + VapourHoldup * vV),
				:(LI*Vtotal= Vfilled),
				:(TI * "K" = OutletVapour.T),
				:(LiquidHoldup * vL = Vfilled),
			],
			[
				"Vessel Cross Section Area","Vessel Total Volume","Vessel Filled Volume","Vessel Filled Volume","Vessel Total Volume","Vessel Filled Volume","Vessel Filled Volume","Vessel Total Volume","Vessel Filled Volume","Component Molar Balance","Energy Balance","Molar Holdup","Energy Holdup","Mol fraction normalisation","Mol fraction normalisation","Liquid Volume","Vapour Volume","Chemical Equilibrium","Thermal Equilibrium","Mechanical Equilibrium for the Vapour Phase","Static Head","Mechanical Equilibrium for the Liquid Phase","Geometry Constraint","Level indicator","Temperature indicator","Liquid Level",
			],
			[:PP,:NComp,:Mw,:pi,:g,:Heads,:Diameter,:Lenght,:Vhead_elliptical,:Vhead_hemispherical,:Vcylinder,:radius,:Levelpercent_Initial,:Temperature_Initial,:Composition_Initial,],
			[:InletLiquid,:OutletLiquid,:InletVapour,:OutletVapour,:InletQ,:Vtotal,:Vfilled,:TotalHoldup,:LiquidHoldup,:VapourHoldup,:E,:vL,:vV,:Level,:Across,:Peq,:Pstatic,:LI,:TI,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Mw::Array{molweight}
	pi::positive
	g::acceleration
	Heads::DanaSwitcher
	Diameter::length
	Lenght::length
	Vhead_elliptical::volume
	Vhead_hemispherical::volume
	Vcylinder::volume
	radius::length
	Levelpercent_Initial::positive
	Temperature_Initial::temperature
	Composition_Initial::Array{fraction}
	InletLiquid::stream
	OutletLiquid::liquid_stream
	InletVapour::stream
	OutletVapour::vapour_stream
	InletQ::power
	Vtotal::volume
	Vfilled::volume
	TotalHoldup::Array{mol}
	LiquidHoldup::mol
	VapourHoldup::mol
	E::energy
	vL::volume_mol
	vV::volume_mol
	Level::length
	Across::area
	Peq::pressure
	Pstatic::pressure
	LI::control_signal
	TI::control_signal
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export SumpTank
function set(in::SumpTank)
	Mw=PP.MolecularWeight()
	 g = 9.81 * "m/(s^2)"
	 Vhead_elliptical = (pi*Diameter^3)/12
	 Vhead_hemispherical = (pi*Diameter^3)/6
	 Vcylinder = 0.25*(pi*Diameter^2)*Lenght
	 radius = 0.5*Diameter

end
function initial(in::SumpTank)
	addEquation(1)
	addEquation(2)
	addEquation(3)
end
function setEquationFlow(in::SumpTank)
	addEquation(1)
	let switch=Heads
		if switch=="elliptical"
			addEquation(2)
			if Level < 0.25*Diameter
				addEquation(3)
			else
				addEquation(4)
			end
		elseif switch=="hemispherical"
			addEquation(5)
			if Level < 0.5*Diameter
				addEquation(6)
			else
				addEquation(7)
			end
		elseif switch=="flat"
			addEquation(8)
			addEquation(9)
		end
	end
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
	addEquation(22)
	#"Pressure Drop"
#	Pdrop = min([InletLiquid.P,InletVapour.P]) - OutletLiquid.P;
#	#OutletLiquid.P  = InletLiquid.P - Pdrop;
#
addEquation(23)
	addEquation(24)
	addEquation(25)
	addEquation(26)
end
function atributes(in::SumpTank,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/SumpTank"
	fields[:Brief]="Model of a tank with 2 material inlet streams and with thermodynamic equilibrium."
	fields[:Info]="== ASSUMPTIONS ==
* perfect mixing of both phases;
* thermodynamics equilibrium.

== SET ==
*Orientation: vessel position - vertical or horizontal;
*Heads (bottom and top heads are identical)
**elliptical: 2:1 elliptical heads (25% of vessel diameter);
**hemispherical: hemispherical heads (50% of vessel diameter);
**flat: flat heads (0% of vessel diameter);
*Diameter: Vessel diameter;
*Lenght: Side length of the cylinder shell;

== SPECIFY ==
* the Inlet stream;
* the outlet flows: OutletVapour.F and OutletLiquid.F;
* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model).

== OPTIONAL ==
* the SumpTank model has three control ports
** TI OutletLiquid Temperature Indicator;
** PI OutletLiquid Pressure Indicator;
** LI Level Indicator;

== INITIAL CONDITIONS ==
* Initial_Temperature :  the Tank temperature (OutletLiquid.T);
* Levelpercent_Initial : the Tank liquid level in percent (LI);
* Initial_Composition : (NoComps) OutletLiquid compositions.
"
	drive!(fields,_)
	return fields
end
SumpTank(_::Dict{Symbol,Any})=begin
	newModel=SumpTank()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
export check_sstate

"""
    check_sstate(model; <options>)

Run a diagnostic test to determine if the steady state solution stored within
the given model object is indeed a solution of the steady state system of
equations.

Return the number of steady state equations that are violated by the current
steady state solution. If `verbose=true`, also display diagnostic information in
the form of listing all bad equations and their residuals.

### Options
Standard options (default values from `model.options`)
  * `verbose`
  * `tol` - an equation is considered satisfied if its residual, in absolute
    value, is smaller than this number.

"""
function check_sstate(model::Model; verbose::Bool = model.options.verbose, tol::Float64 = model.options.tol)
    R, J = global_SS_RJ(model.sstate.values, model)
    bad_eqn = 0
    bad_eqn_str = ""
    for (ind, res) ∈ sort(collect(enumerate(R)), lt = (a, b)->abs(a[2]) > abs(b[2]))
        if abs(res) > tol
            bad_eqn += 1
            if verbose
                eqn_str = ModelBaseEcon.geteqn(ind, model.sstate)
                res_str = @sprintf "%- 10g" res
                bad_eqn_str *= "\nE$ind  res=$(res_str)  $(eqn_str)"
            end
        end
    end
    if bad_eqn > 0 && verbose
        s = bad_eqn > 1 ? "equations are" : "equation is"
        b = bad_eqn > 1 ? "$(bad_eqn) " : ""
        @warn "The following $(b)steady state $(s) not satisfied: $(bad_eqn_str)"
    end
    if bad_eqn == 0 && verbose
        @info "All steady state equations are satisfied."
    end
    return bad_eqn
end
@assert precompile(check_sstate, (Model,))

export diagnose_sstate
"""
    diagnose_sstate([point,] model)

Run diagnostics on the steady state of the given model. If `point` is not given,
then we check the steady state solution stored inside the given model.

Return a tuple of "bad" equations and "bad" variables.

The set of "bad" equations is one that is inconsistent, i.e. there is no
solution. This might happen if the system is overdetermined.

The set of "bad" variables contains variables that cannot be solved uniquely.
This might happen if the system is underdetermined. In this case, try adding
steady state constraints until you get a unique solution. See
[`@steadystate`](@ref ModelBaseEcon.@steadystate) in ModelBaseEcon.

!!! warning "Internal function"
    The output from this function may be difficult to read.<br>
    Call [`check_sstate`](@ref) instead.
"""
@inline diagnose_sstate(model::Model) = diagnose_sstate(model.sstate.values, model)
function diagnose_sstate(point::AbstractVector{Float64}, model::Model)
    tol = model.options.tol
    rr, jj = global_SS_RJ(point, model)
    neqns, nvars = size(jj)
    ff = qr([jj Matrix(I, neqns, neqns); zeros(neqns, nvars) Matrix(I, neqns, neqns)], Val(true))
    sol = Vector{Float64}(undef, neqns + nvars)
    rj = rank(ff.R)
    rhs = ff.Q' * vcat(rr, zeros(neqns))
    sol[ff.p[1:rj]] = ff.R[1:rj,1:rj] \ rhs[1:rj]
    bad = falses(neqns)
    bad_res = abs.(sol[nvars + 1:end])
    bad[bad_res .> tol] .= true
    bad = sort(findall(bad), lt = (a, b)->bad_res[a] > bad_res[b])
    bad_eqn = tuple(("E$(be)  residual=$(bad_res[be])  $(geteqn(be, model.sstate))" for be in bad)...)
    vsyms = Symbol[ModelBaseEcon.ss_symbol(model.sstate, vi) for vi = 1:nvars]
    bad_var_mask = falses(size(vsyms))
    bad_var_mask[ff.p[rj + 1:end]] .= true
    if model.flags.ssZeroSlope
        bad_var_mask[2:2:end] .= false
    end
    for (vi, v) ∈ enumerate(model.allvars)
        if issteady(v)
            bad_var_mask[2vi] = false
        end
        if isshock(v)
            bad_var_mask[2vi-1] = false
            bad_var_mask[2vi] = false
        end
    end
    # bad_var = vsyms[filter(i->i <= nvars, findall(bad_var_mask))]
    bad_var = vsyms[bad_var_mask]
    return bad_eqn, bad_var
end
@assert precompile(diagnose_sstate, (Vector{Float64}, Model))


text:
# This script is licensed under the MIT License (MIT).

using Weave

convert_doc("legend-julia-software-tutorial.jmd", "legend-julia-software-tutorial.jl")
convert_doc("legend-julia-software-tutorial.jmd", "legend-julia-software-tutorial.ipynb")

max_stars_repo_name: fhagemann/legend-julia-tutorial

max_issues_repo_name: fhagemann/legend-julia-tutorial

max_forks_repo_name: fhagemann/legend-julia-tutorial


text:
function plot_hardi(directions, values)
    facets, points = convexhull(directions)
    for (i, v) ∈ enumerate(values)
        points[:, i] .*= v
    end
    mesh(points', facets', color = values, shading = false)
end


text:
module Day23

using LinearAlgebra

inp = Dict()

for line in eachline("day23.txt")
    m = match(r"pos=<(\-?\d+(\.\d+)?),(\-?\d+(\.\d+)?),(\-?\d+(\.\d+)?)>, r=(\d+)", line)
    m = [parse(Int, i) for i in m.captures if i !== nothing]
    inp[[m[1], m[2], m[3]]] = m[4] # pos => r
end

function solvea()
    r, coord = findmax(inp)
    count(x -> norm(x - coord, 1) <= r, keys(inp))
end

function solveb()
    #unfinished
    tmpmat = hcat(keys(inp)...)
    extrema(tmpmat; dims = 2) |> Iterators.product
end

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
println("------------------------------------")
println("Fusion Trees")
println("------------------------------------")
ti = time()
@timedtestset "Fusion trees for $(TensorKit.type_repr(I))" for I in sectorlist
    Istr = TensorKit.type_repr(I)
    N = 5
    out = ntuple(n->randsector(I), N)
    isdual = ntuple(n->rand(Bool), N)
    in = rand(collect(⊗(out...)))
    numtrees = count(n->true, fusiontrees(out, in, isdual))
    while !(0 < numtrees < 30)
        out = ntuple(n->randsector(I), N)
        in = rand(collect(⊗(out...)))
        numtrees = count(n->true, fusiontrees(out, in, isdual))
    end
    it = @constinferred fusiontrees(out, in, isdual)
    @constinferred Nothing iterate(it)
    f = @constinferred first(it)
    @testset "Fusion tree $I: printing" begin
        @test eval(Meta.parse(sprint(show,f))) == f
    end
    @testset "Fusion tree $Istr: braiding" begin
        for in = ⊗(out...)
            for i = 1:N-1
                for f in fusiontrees(out, in, isdual)
                    d1 = @constinferred TK.artin_braid(f, i)
                    @test norm(values(d1)) ≈ 1
                    d2 = empty(d1)
                    for (f1, coeff1) in d1
                        for (f2,coeff2) in TK.artin_braid(f1, i; inv = true)
                            d2[f2] = get(d2, f2, zero(coeff1)) + coeff2*coeff1
                        end
                    end
                    for (f2, coeff2) in d2
                        if f2 == f
                            @test coeff2 ≈ 1
                        else
                            @test isapprox(coeff2, 0; atol = 1e-12, rtol = 1e-12)
                        end
                    end
                end
            end
        end

        f = rand(collect(it))
        d1 = TK.artin_braid(f, 2)
        d2 = empty(d1)
        for (f1, coeff1) in d1
            for (f2,coeff2) in TK.artin_braid(f1, 3)
                d2[f2] = get(d2, f2, zero(coeff1)) + coeff2*coeff1
            end
        end
        d1 = d2
        d2 = empty(d1)
        for (f1, coeff1) in d1
            for (f2,coeff2) in TK.artin_braid(f1, 3; inv = true)
                d2[f2] = get(d2, f2, zero(coeff1)) + coeff2*coeff1
            end
        end
        d1 = d2
        d2 = empty(d1)
        for (f1, coeff1) in d1
            for (f2,coeff2) in TK.artin_braid(f1, 2; inv = true)
                d2[f2] = get(d2, f2, zero(coeff1)) + coeff2*coeff1
            end
        end
        d1 = d2
        for (f1, coeff1) in d1
            if f1 == f
                @test coeff1 ≈ 1
            else
                @test isapprox(coeff1, 0; atol = 1e-12, rtol = 1e-12)
            end
        end
    end
    @testset "Fusion tree $Istr: braiding and permuting" begin
        p = tuple(randperm(N)...,)
        ip = invperm(p)

        levels = ntuple(identity, N)
        d = @constinferred TK.braid(f, levels, p)
        d2 = Dict{typeof(f), valtype(d)}()
        levels2 = p
        for (f2, coeff) in d
            for (f1,coeff2) in TK.braid(f2, levels2, ip)
                d2[f1] = get(d2, f1, zero(coeff)) + coeff2*coeff
            end
        end
        for (f1, coeff2) in d2
            if f1 == f
                @test coeff2 ≈ 1
            else
                @test isapprox(coeff2, 0; atol = 1e-12, rtol = 1e-12)
            end
        end

        if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
            Af = convert(Array, f)
            Afp = permutedims(Af, (p..., N+1))
            Afp2 = zero(Afp)
            for (f1, coeff) in d
                Afp2 .+= coeff .* convert(Array, f1)
            end
            @test Afp ≈ Afp2
        end
    end
    @testset "Fusion tree $Istr: insertat" begin
        N = 4
        out2 = ntuple(n->randsector(I), N)
        in2 = rand(collect(⊗(out2...)))
        isdual2 = ntuple(n->rand(Bool), N)
        f2 = rand(collect(fusiontrees(out2, in2, isdual2)))
        for i = 1:N
            out1 = ntuple(n->randsector(I), N)
            out1 = Base.setindex(out1, in2, i)
            in1 = rand(collect(⊗(out1...)))
            isdual1 = ntuple(n->rand(Bool), N)
            isdual1 = Base.setindex(isdual1, false, i)
            f1 = rand(collect(fusiontrees(out1, in1, isdual1)))

            trees = @constinferred TK.insertat(f1, i, f2)
            @test norm(values(trees)) ≈ 1

            f1a, f1b = @constinferred TK.split(f1, $i)
            @test length(TK.insertat(f1b, 1, f1a)) == 1
            @test first(TK.insertat(f1b, 1, f1a)) == (f1 => 1)

            levels = ntuple(identity, N)
            gen = Base.Generator(TK.braid(f1, levels, (i, (1:i-1)..., (i+1:N)...))) do (t, c)
                (t′, c′) = first(TK.insertat(t, 1, f2))
                @test c′ == one(c′)
                return t′ => c
            end
            trees2 = Dict(gen)
            trees3 = empty(trees2)
            p = ((N+1:N+i-1)..., (1:N)..., (N+i:2N-1)...)
            levels = ((i:N+i-1)..., (1:i-1)..., (i+N:2N-1)...)
            for (t,coeff) in trees2
                for (t′,coeff′) in TK.braid(t, levels, p)
                    trees3[t′] = get(trees3, t′, zero(coeff′)) + coeff*coeff′
                end
            end
            for (t, coeff) in trees3
                coeff′ = get(trees, t, zero(coeff))
                @test isapprox(coeff′, coeff; atol = 1e-12, rtol = 1e-12)
            end

            if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
                Af1 = convert(Array, f1)
                Af2 = convert(Array, f2)
                Af = TensorOperations.tensorcontract(Af1, [1:i-1; -1; N-1 .+ (i+1:N+1)],
                                                     Af2, [i-1 .+ (1:N); -1], 1:2N)
                Af′ = zero(Af)
                for (f, coeff) in trees
                    Af′ .+= coeff .* convert(Array, f)
                end
                @test isapprox(Af, Af′; atol = 1e-12, rtol = 1e-12)
            end
        end
    end
    @testset "Fusion tree $Istr: merging" begin
        N = 3
        out1 = ntuple(n->randsector(I), N)
        in1 = rand(collect(⊗(out1...)))
        f1 = rand(collect(fusiontrees(out1, in1)))
        out2 = ntuple(n->randsector(I), N)
        in2 = rand(collect(⊗(out2...)))
        f2 = rand(collect(fusiontrees(out2, in2)))

        @constinferred TK.merge(f1, f2, first(in1 ⊗ in2), 1)
        if !(FusionStyle(I) isa DegenerateNonAbelian)
            @constinferred TK.merge(f1, f2, first(in1 ⊗ in2), nothing)
            @constinferred TK.merge(f1, f2, first(in1 ⊗ in2))
        end
        @test dim(in1)*dim(in2) ≈ sum(abs2(coeff)*dim(c) for c in in1 ⊗ in2
                                                    for μ in 1:Nsymbol(in1, in2, c)
                                                    for (f,coeff) in TK.merge(f1, f2, c, μ))

        for c in in1 ⊗ in2
            R = Rsymbol(in1, in2, c)
            for μ = 1:Nsymbol(in1, in2, c)
                μ′ = FusionStyle(I) isa DegenerateNonAbelian ? μ : nothing
                trees1 = TK.merge(f1, f2, c, μ′)

                # test merge and braid interplay
                trees2 = Dict{keytype(trees1), complex(valtype(trees1))}()
                trees3 = Dict{keytype(trees1), complex(valtype(trees1))}()
                for ν = 1:Nsymbol(in2, in1, c)
                    ν′ = FusionStyle(I) isa DegenerateNonAbelian ? ν : nothing
                    for (t, coeff) in TK.merge(f2, f1, c, ν′)
                        trees2[t] = get(trees2, t, zero(valtype(trees2))) + coeff*R[μ,ν]
                    end
                end
                perm = ( (N.+(1:N))... , (1:N)...)
                levels = ntuple(identity, 2*N)
                for (t, coeff) in trees1
                    for (t′, coeff′) in TK.braid(t, levels, perm)
                        trees3[t′] = get(trees3, t′, zero(valtype(trees3))) + coeff*coeff′
                    end
                end
                for (t, coeff) in trees3
                    coeff′ = get(trees2, t, zero(coeff))
                    @test isapprox(coeff, coeff′; atol = 1e-12, rtol = 1e-12)
                end

                # test via conversion
                if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
                    Af1 = convert(Array, f1)
                    Af2 = convert(Array, f2)
                    Af0 = convert(Array, FusionTree((f1.coupled, f2.coupled), c, (false, false), (), (μ,)))
                    _Af = TensorOperations.tensorcontract(Af1, [1:N;-1],
                                                            Af0, [-1;N+1;N+2], 1:N+2)
                    Af = TensorOperations.tensorcontract(Af2, [N .+ (1:N); -1],
                                                            _Af, [1:N; -1; 2N+1], 1:2N+1)
                    Af′ = zero(Af)
                    for (f, coeff) in trees1
                        Af′ .+= coeff .* convert(Array, f)
                    end
                    @test Af ≈ Af′
                end
            end
        end
    end

    if I <: ProductSector
        N = 3
    else
        N = 4
    end
    out = ntuple(n->randsector(I), N)
    numtrees = count(n->true, fusiontrees((out..., map(dual, out)...)))
    while !(0 < numtrees < 100)
        out = ntuple(n->randsector(I), N)
        numtrees = count(n->true, fusiontrees((out..., map(dual, out)...)))
    end
    incoming = rand(collect(⊗(out...)))
    f1 = rand(collect(fusiontrees(out, incoming, ntuple(n->rand(Bool), N))))
    f2 = rand(collect(fusiontrees(out[randperm(N)], incoming, ntuple(n->rand(Bool), N))))

    @testset "Double fusion tree $Istr: repartioning" begin
        for n = 0:2*N
            d = @constinferred TK.repartition(f1, f2, $n)
            @test dim(incoming) ≈ sum(abs2(coef)*dim(f1.coupled) for ((f1,f2), coef) in d)
            d2 = Dict{typeof((f1,f2)), valtype(d)}()
            for ((f1′,f2′),coeff) in d
                for ((f1′′,f2′′),coeff2) in TK.repartition(f1′,f2′, N)
                    d2[(f1′′,f2′′)] = get(d2, (f1′′,f2′′), zero(coeff)) + coeff2*coeff
                end
            end
            for ((f1′,f2′), coeff2) in d2
                if f1 == f1′ && f2 == f2′
                    @test coeff2 ≈ 1
                else
                    @test isapprox(coeff2, 0; atol = 1e-12, rtol = 1e-12)
                end
            end
            if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
                Af1 = convert(Array, f1)
                Af2 = permutedims(convert(Array, f2), [N:-1:1; N+1])
                sz1 = size(Af1)
                sz2 = size(Af2)
                d1 = prod(sz1[1:end-1])
                d2 = prod(sz2[1:end-1])
                dc = sz1[end]
                A = reshape(reshape(Af1, (d1, dc)) * reshape(Af2, (d2, dc))',
                            (sz1[1:end-1]..., sz2[1:end-1]...))
                A2 = zero(A)
                for ((f1′,f2′), coeff) in d
                    Af1′ = convert(Array, f1′)
                    Af2′ = permutedims(convert(Array, f2′), [(2N-n):-1:1; 2N-n+1])
                    sz1′ = size(Af1′)
                    sz2′ = size(Af2′)
                    d1′ = prod(sz1′[1:end-1])
                    d2′ = prod(sz2′[1:end-1])
                    dc′ = sz1′[end]
                    A2 += coeff*reshape(reshape(Af1′,(d1′,dc′)) * reshape(Af2′,(d2′,dc′))',
                                        (sz1′[1:end-1]..., sz2′[1:end-1]...))
                end
                @test A ≈ A2
            end
        end
    end
    @testset "Double fusion tree $Istr: permutation" begin
        if BraidingStyle(I) isa SymmetricBraiding
            for n = 0:2N
                p = (randperm(2*N)...,)
                p1, p2 = p[1:n], p[n+1:2N]
                ip = invperm(p)
                ip1, ip2 = ip[1:N], ip[N+1:2N]

                d = @constinferred TensorKit.permute(f1, f2, p1, p2)
                @test dim(incoming) ≈ sum(abs2(coef)*dim(f1.coupled) for ((f1,f2), coef) in d)
                d2 = Dict{typeof((f1,f2)), valtype(d)}()
                for ((f1′,f2′), coeff) in d
                    d′ = TensorKit.permute(f1′,f2′, ip1, ip2)
                    for ((f1′′,f2′′), coeff2) in d′
                        d2[(f1′′,f2′′)] = get(d2, (f1′′,f2′′), zero(coeff)) + coeff2*coeff
                    end
                end
                for ((f1′,f2′), coeff2) in d2
                    if f1 == f1′ && f2 == f2′
                        @test coeff2 ≈ 1
                    else
                        @test abs(coeff2) < 1e-12
                    end
                end

                if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
                    Af1 = convert(Array, f1)
                    Af2 = convert(Array, f2)
                    sz1 = size(Af1)
                    sz2 = size(Af2)
                    d1 = prod(sz1[1:end-1])
                    d2 = prod(sz2[1:end-1])
                    dc = sz1[end]
                    A = reshape(reshape(Af1, (d1, dc)) * reshape(Af2, (d2, dc))',
                                (sz1[1:end-1]..., sz2[1:end-1]...))
                    Ap = permutedims(A, (p1..., p2...));
                    A2 = zero(Ap)
                    for ((f1′,f2′), coeff) in d
                        Af1′ = convert(Array, f1′)
                        Af2′ = convert(Array, f2′)
                        sz1′ = size(Af1′)
                        sz2′ = size(Af2′)
                        d1′ = prod(sz1′[1:end-1])
                        d2′ = prod(sz2′[1:end-1])
                        dc′ = sz1′[end]
                        A2 += coeff*reshape(reshape(Af1′,(d1′,dc′)) *
                                            reshape(Af2′,(d2′,dc′))',
                                            (sz1′[1:end-1]..., sz2′[1:end-1]...))
                    end
                    @test Ap ≈ A2
                end
            end
        end
    end
    @testset "Double fusion tree $Istr: transposition" begin
        for n = 0:2N
            i0 = rand(1:2N)
            p = mod1.(i0 .+ (1:2N), 2N)
            ip = mod1.(-i0 .+ (1:2N), 2N)
            p′ = tuple(getindex.(Ref(vcat(1:N, 2N:-1:N+1)), p)...)
            p1, p2 = p′[1:n], p′[2N:-1:n+1]
            ip′ = tuple(getindex.(Ref(vcat(1:n, 2N:-1:n+1)), ip)...)
            ip1, ip2 = ip′[1:N], ip′[2N:-1:N+1]

            d = @constinferred transpose(f1, f2, p1, p2)
            @test dim(incoming) ≈ sum(abs2(coef)*dim(f1.coupled) for ((f1,f2), coef) in d)
            d2 = Dict{typeof((f1,f2)), valtype(d)}()
            for ((f1′,f2′), coeff) in d
                d′ = transpose(f1′,f2′, ip1, ip2)
                for ((f1′′,f2′′), coeff2) in d′
                    d2[(f1′′,f2′′)] = get(d2, (f1′′,f2′′), zero(coeff)) + coeff2*coeff
                end
            end
            for ((f1′,f2′), coeff2) in d2
                if f1 == f1′ && f2 == f2′
                    @test coeff2 ≈ 1
                else
                    @test abs(coeff2) < 1e-12
                end
            end

            if BraidingStyle(I) isa Bosonic
                d3 = permute(f1, f2, p1, p2)
                for (f1′,f2′) in union(keys(d), keys(d3))
                    coeff1 = get(d, (f1′,f2′), zero(valtype(d)))
                    coeff3 = get(d3, (f1′,f2′), zero(valtype(d3)))
                    @test isapprox(coeff1, coeff3; atol = 1e-12)
                end
            end

            if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
                Af1 = convert(Array, f1)
                Af2 = convert(Array, f2)
                sz1 = size(Af1)
                sz2 = size(Af2)
                d1 = prod(sz1[1:end-1])
                d2 = prod(sz2[1:end-1])
                dc = sz1[end]
                A = reshape(reshape(Af1, (d1, dc)) * reshape(Af2, (d2, dc))',
                            (sz1[1:end-1]..., sz2[1:end-1]...))
                Ap = permutedims(A, (p1..., p2...));
                A2 = zero(Ap)
                for ((f1′,f2′), coeff) in d
                    Af1′ = convert(Array, f1′)
                    Af2′ = convert(Array, f2′)
                    sz1′ = size(Af1′)
                    sz2′ = size(Af2′)
                    d1′ = prod(sz1′[1:end-1])
                    d2′ = prod(sz2′[1:end-1])
                    dc′ = sz1′[end]
                    A2 += coeff*reshape(reshape(Af1′,(d1′,dc′)) *
                                        reshape(Af2′,(d2′,dc′))',
                                        (sz1′[1:end-1]..., sz2′[1:end-1]...))
                end
                @test Ap ≈ A2
            end
        end
    end
end
tf = time()
printstyled("Finished fusion tree tests in ",
            string(round(tf-ti; sigdigits=3)),
            " seconds."; bold = true, color = Base.info_color())
println()


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _mag2flux(mag::T, zero_point::T, ABwave::T) where {T<:AbstractFloat}
    if isnan(ABwave)
        return exp10(-0.4*(mag + zero_point))
    else
        return exp10(-0.4*(mag + 2.406 + 5*log10(ABwave)))
    end
end

"""
    mag2flux(mag[, zero_point, ABwave=number]) -> flux

### Purpose ###

Convert from magnitudes to flux expressed in erg/(s cm² Å).

### Explanation ###

This is the reverse of `flux2mag`.

### Arguments ###

* `mag`: the magnitude to be converted in flux.
* `zero_point`: the zero point level of the magnitude.  If not supplied then defaults to
 21.1 (Code et al 1976).  Ignored if the `ABwave` keyword is supplied
* `ABwave` (optional numeric keyword): wavelength, in Angstroms.  If supplied, then the
 input `mag` is assumed to contain Oke AB magnitudes (Oke & Gunn 1983, ApJ, 266, 713;
 http://adsabs.harvard.edu/abs/1983ApJ...266..713O).

### Output ###

The flux.

If the `ABwave` keyword is set, then the flux is given by the expression

\$\$\\text{flux} = 10^{-0.4(\\text{mag} +2.406 + 4\\log_{10}(\\text{ABwave}))}\$\$

Otherwise the flux is given by

\$\$\\text{flux} =  10^{-0.4(\\text{mag} + \\text{zero point})}\$\$

### Example ###

```jldoctest
julia> using AstroLib

julia> mag2flux(8.3)
1.7378008287493692e-12

julia> mag2flux(8.3, 12)
7.58577575029182e-9

julia> mag2flux(8.3, ABwave=12)
3.6244115683017193e-7
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
mag2flux(mag::Real, zero_point::Real=21.1; ABwave::Real=NaN) =
    _mag2flux(promote(float(mag), float(zero_point), float(ABwave))...,)


text:
using Test
using ObservableCollector

output = []

@test begin
    output = []
    obs = @observations begin
        @every 2 "O2" (s,t)->t
    end
    for t in 1:10
        obs(output,0,t)
    end
    output == [(name=:O2, val=2j) for j in 1:5]
end

@test begin
    output = []
    obs = @observations begin
        @at 2 "AT" (s,t)->t
    end
    for t in 1:10
        obs(output,0,t)
    end
    output == [(name=:AT, val=2)]
end

@test begin
    output = []
    obs = @observations begin
        @condition (s,t)->s+t==5 "sumtofive" (s,t)->t
    end
    for t in 1:10
        obs(output,1,t)
    end
    output == [(name=:sumtofive, val=4)]
end

@test begin
    output = []
    obs = @observations begin
        @condition (s,t)->t>=3 begin
            "A" --> (s,t)->s
            "T" --> (s,t)->t
        end
    end
    for t in 1:4
        obs(output,1,t)
    end
    output == [(name=:A, val=1),(name=:T, val=3),(name=:A, val=1),(name=:T, val=4)]
end


text:
@testset "IO" begin
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x^2 - 2;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [2, 0, 0, 0, 1, 0, 0, 0, 16]);
  gens = [[2, 0, 0], [2, 0, 0], [2, 2, 0], [a, a, 0], [a + 2, 1, 1//4], [3*a + 3, 3//2*a, 3//8*a]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test Hecke.to_hecke_string(L) isa String
  @test Hecke.to_magma_string(L) isa String

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^4-x^3-4*x^2+4*x+1
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2+(a^3 - 1*a^2 - 4*a + 5)
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 3, 3, [(1), 0, 0, 0, (1), 0, 0, 0, (1)])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [(1), 0, 0]), map(E, [0, (1), 0]), map(E, [0, 0, (1)])]
  L = hermitian_lattice(E, generators = gens, gram_ambient_space = D)
  @test Hecke.to_hecke_string(L) isa String
  @test Hecke.to_magma_string(L) isa String
end


text:
using Quadrature
using Cubature, Cuba
using Test

singledim_algs = [
        QuadGKJL(),
        HCubatureJL(),
        VEGAS(),
        CubatureJLh(),
        CubatureJLp(),
        CubaVegas(),
        CubaSUAVE()
        ]

multidim_algs = [
        HCubatureJL(),
        VEGAS(),
        CubatureJLh(),
        CubatureJLp(),
        CubaVegas(),
        CubaSUAVE(),
        CubaDivonne(),
        CubaCuhre()
        ]

singledim_batch_algs = [
        VEGAS(),
        CubatureJLh(),
        CubatureJLp(),
        CubaVegas(),
        CubaSUAVE()
        ]

multidim_batch_algs = [
        VEGAS(),
        CubatureJLh(),
        CubatureJLp(),
        CubaVegas(),
        CubaSUAVE(),
        CubaDivonne(),
        CubaCuhre()
        ]

integrands = [
              (x,p) -> 1,
              (x,p) -> sum(x),
              (x,p) -> sum(sin.(x))
             ]

iip_integrands = [
                  (dx,x,p) -> (dx .= 1),
                  (dx,x,p) -> (dx .= sum(x)),
                  (dx,x,p) -> (dx .= sum(sin.(x)))
              ]

batch_f(f) = (pts,p) -> begin
  fevals = zeros(size(pts, 2))
  for i = 1:size(pts, 2)
     x = pts[:,i]
     fevals[i] = f(x,p)
  end
  fevals
end

batch_iip_f(f) = (fevals,pts,p) -> begin
  for i = 1:size(pts, 2)
     x = pts[:,i]
     fevals[i] = f(x,p)
  end
  nothing
end

prob = QuadratureProblem(integrands[1],ones(2),3ones(2))
_sol = solve(prob,HCubatureJL())
solve(prob,CubatureJLh(),reltol=1e-3,abstol=1e-3)

@testset "Standard Single Dimension Integrands" begin
    for i in 1:length(integrands)
        prob = QuadratureProblem(integrands[i],1,3)
        _sol = solve(prob,HCubatureJL())
        for alg in singledim_algs
            @info "Dimension = 1, Alg = $alg, Integrand = $i"
            sol = solve(prob,alg,reltol=1e-3,abstol=1e-3)
            @test sol.u ≈ _sol.u rtol = 1e-2
        end
    end
end

@testset "Standard Integrands" begin
    for i in 1:length(integrands)
        for dim = 1:5
            prob = QuadratureProblem(integrands[i],ones(dim),3ones(dim))
            _sol = solve(prob,CubatureJLh())
            for alg in multidim_algs
                if dim == 1 && (alg isa CubaDivonne || alg isa CubaCuhre) ||
                    dim > 3 && alg isa VEGAS
                    continue
                end
                @info "Dimension = $dim, Alg = $alg, Integrand = $i"
                sol = solve(prob,alg,reltol=1e-3,abstol=1e-3)
                @test sol.u ≈ _sol.u rtol = 1e-2
            end
        end
    end
end

@testset "In-place Standard Integrands" begin
    for i in 1:length(iip_integrands)
        for dim = 1:5
            prob = QuadratureProblem(iip_integrands[i],ones(dim),3ones(dim))
            _sol = solve(prob,CubatureJLh())
            for alg in multidim_algs
                if dim == 1 && (alg isa CubaDivonne || alg isa CubaCuhre) ||
                    dim > 3 && alg isa VEGAS # Large VEGAS omitted because it's slow!
                    continue
                end
                @info "Dimension = $dim, Alg = $alg, Integrand = $i"
                sol = solve(prob,alg,reltol=1e-3,abstol=1e-3)
                @test sol.u ≈ _sol.u rtol = 1e-2
            end
        end
    end
end

@testset "Batched Single Dimension Integrands" begin
    for i in 1:length(integrands)
        _prob = QuadratureProblem(integrands[i],1,3)
        _sol = solve(_prob,CubatureJLh())
        prob = QuadratureProblem(batch_f(integrands[i]),1,3,batch=10)
        for alg in singledim_batch_algs
            @info "Dimension = 1, Alg = $alg, Integrand = $i"
            sol = solve(prob,alg,reltol=1e-3,abstol=1e-3)
            if !isapprox(_sol.u,sol.u,rtol = 1e-2)
                @info "$(alg) failed"
                @show sol.u,_sol.u
                @test_broken sol.u ≈ _sol.u rtol = 1e-2
            else
                @test sol.u ≈ _sol.u rtol = 1e-2
            end
        end
    end
end

@testset "Batched Standard Integrands" begin
    for i in 1:length(integrands)
        for dim = 1:5
            _prob = QuadratureProblem(integrands[i],ones(dim),3ones(dim))
            _sol = solve(_prob,CubatureJLh())
            prob = QuadratureProblem(batch_f(integrands[i]),ones(dim),3ones(dim),batch=10)
            for alg in multidim_batch_algs
                if dim == 1 && (alg isa CubaDivonne || alg isa CubaCuhre) ||
                   dim > 3 && alg isa VEGAS # Large VEGAS omitted because it's slow!
                    continue
                end
                @info "Dimension = $dim, Alg = $alg, Integrand = $i"
                sol = solve(prob,alg,reltol=1e-3,abstol=1e-3)
                @test sol.u ≈ _sol.u rtol = 1e-2
            end
        end
    end
end

@testset "In-Place Batched Standard Integrands" begin
    for i in 1:length(iip_integrands)
        for dim = 1:5
            _prob = QuadratureProblem(iip_integrands[i],ones(dim),3ones(dim))
            _sol = solve(_prob,CubatureJLh())
            prob = QuadratureProblem(batch_iip_f(integrands[i]),ones(dim),3ones(dim),batch=10)
            for alg in multidim_batch_algs
                if dim == 1 && (alg isa CubaDivonne || alg isa CubaCuhre)
                    continue
                end
                @info "Dimension = $dim, Alg = $alg, Integrand = $i"
                if alg isa VEGAS
                    @test_broken sol = solve(prob,alg,reltol=1e-3,abstol=1e-3).retcode == :Success
                else
                    sol = solve(prob,alg,reltol=1e-3,abstol=1e-3)
                    @test sol.u ≈ _sol.u rtol = 1e-2
                end
            end
        end
    end
end


text:
using Base.Threads

function softmax!(out::AbstractVecOrMat{T}, xs::AbstractVecOrMat{T}) where T<:AbstractFloat
  @threads for j = 1:size(xs, 2)
    @inbounds begin
      # out[end, :] .= maximum(xs, 1)
      out[end, j] = xs[end, j]
      for i = 1:size(xs, 1)
        out[end, j] = max(out[end, j], xs[i, j])
      end
      # out .= exp(xs .- out[end, :])
      for i = 1:size(out, 1)
        out[i, j] = exp(xs[i, j] - out[end, j])
      end
      # out ./= sum(out, 1)
      s = zero(eltype(out))
      for i = 1:size(out, 1)
        s += out[i, j]
      end
      for i = 1:size(out, 1)
        out[i, j] /= s
      end
    end
  end
  return out
end

softmax!(xs) = softmax!(xs, xs)
softmax(xs) = softmax!(similar(xs), xs)

function ∇softmax!(out::AbstractVecOrMat, Δ::AbstractVecOrMat, xs::AbstractVecOrMat)
  sf = softmax(xs)
  out .= sf .* (Δ .- sum(Δ .*sf, dims = 1))
end

∇softmax!(Δ, xs) = ∇softmax!(Δ, Δ, xs)
∇softmax(Δ, xs) = ∇softmax!(similar(Δ), Δ, xs)

"""
    softmax(xs) = exp.(xs) ./ sum(exp.(xs))

[Softmax](https://en.wikipedia.org/wiki/Softmax_function) takes
log-probabilities (any real vector) and returns a probability distribution that
sums to 1.

If given a matrix it will treat it as a batch of vectors, with each column
independent.

    julia> softmax([1,2,3.])
    3-element Array{Float64,1}:
      0.0900306
      0.244728
      0.665241
"""
softmax


text:
mutable struct Context

end

ytest(u,i,C)=mean(u)+i
ytest()="Y"
xtest(u,x,j,i,C)=sum(u)
xtest()="U"
dc(u,x,xc,j,i,C)=3*u[j]+x[1]-xc[1]
dc()="u"

function du(u,x,xc,j,i,C)
  Y=x[1]
  U=xc[1]
end

c=(x=[xtest],f=(k=1.0,r=0.1),deriv=dc)

C=(x=[ytest],c=(c,c,c),fixed=(e=2.0,))

C.c.x[1]()

function dudt(du,u,C)
  #each instance i test
  i=1
  x=zeros(1)
    [x[q]=C.x[q](u,i,C) for q in eachindex(C.x)]
    for j in 1:length(C.c)
        xc=[x(u,x,j,i,C) for x in C.c[i].x]
        du[j]=c.deriv(u,x,xc,j,i,C)
        use splatting with an id vector of variable names!
    end
  #end
  return du
end

u=rand(3)
du=similar(u)
dudt(du,u,C)

#get argument names:
e=Base.uncompressed_ast(collect(methods(f))[1]).slotnames
#Get function name
typeof(f).name.mt.name

thisIsAFunction_τ(a,b,g)=a+b+g
d=thisIsAFunction_τ
d(1,2,3)
string(typeof(d).name.mt.name)
#=
first register all functions
create a hashtable for id's to varnames
create nice labels & variable from function names
create idlist for all functions
=#


text:
"""
TLeakyIntegrator: Integration with a constant leak
(No leakage of background non-currently experienced s-a pairs)

NOT USED IN PAPER!
"""
struct TLeakyIntegrator <: TNs1a0s0
    nr_states::Int
    nr_actions::Int
    etaleak::Float64
    lowerbound::Float64
    SPE::Array{Float64, 1}
    Nsa::Array{Float64, 2}
    Ns1a0s0::Array{Dict{Tuple{Int, Int}, Float64}, 1}
    terminalstates::Array{Int,1}
end
function TLeakyIntegrator(; nr_states = 10, nr_actions = 4, etaleak = .9, lowerbound = eps())
    SPE = [0.]
    Nsa = zeros(nr_actions, nr_states) .+ nr_states*eps()
    Ns1a0s0 = [Dict{Tuple{Int, Int}, Float64}() for _ in 1:nr_states]
    [Ns1a0s0[sprime][(a, s)] = eps() for sprime in 1:nr_states for a in 1:nr_actions for s in 1:nr_states]
    TLeakyIntegrator(nr_states, nr_actions, etaleak, lowerbound, SPE, Nsa, Ns1a0s0, Int[])
end
export TLeakyIntegrator
""" X[t] = etaleak * X[t-1] + etaleak * I[.] """
function updatet!(learnerT::TLeakyIntegrator, s0, a0, s1, done)
    calcSPE!(learnerT, s0, a0, s1)
    leaka0s0!(learnerT, s0, a0, s1, done)
    computeterminalNs1a0s0!(learnerT, s1, done)
end
export updatet!
function leaka0s0!(learnerT::TLeakyIntegrator, s0, a0, s1, done)
    learnerT.Nsa[a0, s0] *= learnerT.etaleak # Discount transition
    learnerT.Nsa[a0, s0] += learnerT.etaleak # Increase observed transition
    learnerT.Nsa[a0, s0] += (1. - learnerT.etaleak) * learnerT.nr_states * learnerT.lowerbound # Bound it

    nextstates = [s for s in 1:learnerT.nr_states if haskey(learnerT.Ns1a0s0[s],(a0,s0))]
    for sprime in nextstates
        learnerT.Ns1a0s0[sprime][(a0, s0)] *= learnerT.etaleak # Discount all outgoing transitions
        learnerT.Ns1a0s0[sprime][(a0, s0)] += (1. - learnerT.etaleak) * learnerT.lowerbound # Bound it
    end
    if haskey(learnerT.Ns1a0s0[s1], (a0, s0))
        learnerT.Ns1a0s0[s1][(a0, s0)] += learnerT.etaleak # Increase observed
    else

        learnerT.Ns1a0s0[s1][(a0, s0)] = learnerT.etaleak
        learnerT.Ns1a0s0[s1][(a0, s0)] += (1. - learnerT.etaleak) * learnerT.lowerbound # Bound it
    end
end
function computeterminalNs1a0s0!(learnerT::TLeakyIntegrator, s1, done)
    if done
        if !in(s1, learnerT.terminalstates)
            push!(learnerT.terminalstates, s1)
            for a in 1:learnerT.nr_actions
                for s in 1:learnerT.nr_states
                    if s == s1
                        learnerT.Nsa[a, s1] = 1.
                        learnerT.Ns1a0s0[s][(a, s1)] = 1.
                    else
                        learnerT.Ns1a0s0[s][(a, s1)] = 0.
                    end
                end
            end
        end
    end
end
function calcSPE!(learnerT::TLeakyIntegrator, s0, a0, s1)
    Psasprime = 0.
    if learnerT.Nsa[a0, s0]!=0
        if haskey(learnerT.Ns1a0s0[s1], (a0, s0))
            Psasprime = learnerT.Ns1a0s0[s1][(a0, s0)]/learnerT.Nsa[a0, s0]
        end
    end
    learnerT.SPE[1] = 1. - Psasprime
end


text:
using ToyFHE
using OffsetArrays
using Test

# For now use the BFV parameter estimator and just copy it over.
# We need to figure out the right way to do this in CKKS
bfv_params = BFVParams(
    7, # plaintext modulus
    ; eval_mult_count = 3
)

params = CKKSParams(bfv_params.ℛ, bfv_params.relin_window,
    bfv_params.σ)

scale = 2^40
Tscale = FixedRational{scale}
Tscalesq = FixedRational{Int128(scale)^2}

plain = CKKSEncoding{Tscale}(zero(plaintext_space(params)))
plain .= OffsetArray(LinRange(0.0, 1.0, 2048), 0:2047)

re = convert(ToyFHE.NTT.RingElement, plain)

# TODO: Can we do better here by being more careful with the FFT?
@test real.(collect(CKKSEncoding{Tscalesq}(re*re))) ≈ LinRange(0.0, 1.0, 2048).^2 atol=10^-4

# First test the encoder in isolation without the extra noise from encryption

kp = keygen(params)
c = encrypt(kp, plain)

@test real.(collect(decrypt(kp, c))) ≈ LinRange(0.0, 1.0, 2048) atol=10^-4
@test real.(collect(decrypt(kp, c*c))) ≈ LinRange(0.0, 1.0, 2048).^2 atol=10^-4


text:
push!(LOAD_PATH,"../src/")
using Documenter, GL2017Replication

makedocs(modules = [GL2017Replication], sitename = "GL2017Replication.jl"
        ,pages = [
                "Home" => "index.md",
                "Steady State" => "steadystate.md",
                "Transition" => "transition.md",
                "Output" => "output.md",
                "Functions" => "functions.md",
                ]
                )

deploydocs(repo = "github.com/ChristopherSchang/GL2017Replication.jl.git", devbranch = "main")


text:
"""
    countBinaryOnes(ind::Individual)

documentation
"""
function countBinaryOnes(genotype::GeneticJulia.BinaryGenotype)
    genotype = genotype._representation
    tam = length(genotype)
    count=0
    for i=1:tam
        if genotype[i]==1
            count+=1
        end
    end
    return count
end # function

function nNodes(canonicalgp::GeneticJulia.GPGenotype, gpInfo)
    length(canonicalgp._representation)
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
module AirfoilMeshers

using Statistics, Printf, LinearAlgebra
using StaticArrays
using Dierckx
using TimerOutputs
using FortranFiles
struct AirfoilMesh
    topo::Symbol
    grid::Matrix{MVector{2,Float64}}
    ncon::Int
end

const to = TimerOutput()

include("tridiag/solver.jl")
include("surface/surface.jl")
include("hyperbolic/hyperbolic.jl")
include("output/output.jl")

export SurfaceDistribution, TanhSpacing, EqualSpacing
export LeadingNode, TrailingNode, SurfaceNode
export surface_mesh
export WakeInfo, create_wake, join_wake
export StepSize, UserStepSize, SpecifiedStepSize, ExpStepSize, ConstStepSize, MinStepSize, FilterStepSize
export ConstantX,ConstantY,ConstantAngle,FreeBoundary,SymmetryBoundary,LineBoundary
export HyperbolicMesher
export gen_mesh, save_mesh

end # module


text:
using Test, Mill, HierarchicalUtils

@testset "replacein" begin
    metadata = fill("metadata", 4)
    an1 = ArrayNode(rand(3,4))
    b = BagNode(an1, [1:4, 0:-1], metadata)
    an2 = ArrayNode(randn(5, 4))
    wb = WeightedBagNode(an2, [1:2,3:4], rand(1:4, 4), metadata)
    pn = ProductNode((b=b,wb=wb))
    an3 = ArrayNode(rand(10, 2))
    x1 = ProductNode((pn, an3))

    t = list_traversal(x1)
    for t in list_traversal(x1)
        x2 = deepcopy(x1)
        x3 = replacein(x2, x2[t], x1[t])
        @test x3[t] === x1[t]
    end

    m1 = reflectinmodel(x1)
    for t in list_traversal(m1)
        m2 = deepcopy(m1)
        m3 = replacein(m2, m2[t], m1[t])
        @test m3[t] === m1[t]
    end

end

@testset "findin" begin
    metadata = fill("metadata", 4)
    an1 = ArrayNode(rand(3,4))
    b = BagNode(an1, [1:4, 0:-1], metadata)
    an2 = ArrayNode(randn(5, 4))
    wb = WeightedBagNode(an2, [1:2,3:4], rand(1:4, 4), metadata)
    pn = ProductNode((b=b,wb=wb))
    an3 = ArrayNode(rand(10, 2))
    x1 = ProductNode((a = pn, b = an3))

    t = list_traversal(x1)
    for t in list_traversal(x1)
        l = findin(x1, x1[t])
        @test get(x1,l) === x1[t]
    end

    m1 = reflectinmodel(x1)
    for t in list_traversal(m1)
        l = findin(m1,  m1[t])
        @test get(m1, l) === m1[t]
    end

end


text:
"""
Orthogonalize!(V, hv, vv, orth; verbose=false)

Orthogonalize the Krylov vectors using your (my) choice of
methods. Anything other than classical Gram-Schmidt twice (cgs2) is
likely to become an undocumented and UNSUPPORTED option. Methods other
than cgs2 are mostly for CI for the linear solver.

DO NOT use anything other than "cgs2" with Anderson acceleration.
"""
function Orthogonalize!(V, hv, vv, orth = "cgs2"; verbose=false)
    orthopts = ["mgs1", "mgs2", "cgs1", "cgs2"]
    orth in orthopts || error("Impossible orth spec in Orthogonalize!")
    if orth == "mgs1"
        mgs!(V, hv, vv; verbose=verbose)
    elseif orth == "mgs2"
        mgs!(V, hv, vv, "twice"; verbose=verbose)
    elseif orth == "cgs1"
        cgs!(V, hv, vv, "once"; verbose=verbose)
    else
        cgs!(V, hv, vv, "twice"; verbose=verbose)
    end
end

"""
mgs!(V, hv, vv, orth; verbose=false)
"""
function mgs!(V, hv, vv, orth = "once"; verbose=false)
    k = length(hv) - 1
    normin = norm(vv)
    #p=copy(vv)
    @views for j = 1:k
        p = vec(V[:, j])
        hv[j] = p' * vv
        vv .-= hv[j] * p
    end
    hv[k+1] = norm(vv)
    if (normin + 0.001 * hv[k+1] == normin) && (orth == "twice")
        @views for j = 1:k
            p = vec(V[:, j])
            hr = p' * vv
            hv[j] += hr
            vv .-= hr * p
        end
        hv[k+1] = norm(vv)
    end
    nv = hv[k+1]
    #
    # Watch out for happy breakdown
    #
    #if hv[k+1] != 0
    #@views vv .= vv/hv[k+1]
    (nv != 0) || (verbose && (println("breakdown in mgs1")))
    if nv != 0
        vv ./= nv
    end
end

"""
cgs!(V, hv, vv, orth="twice"; verbose=false)

Classical Gram-Schmidt.
"""
function cgs!(V, hv, vv, orth="twice"; verbose=false)
    #
    #   no BLAS. Seems faster than BLAS since 1.6 and allocates
    #   far less memory.
    #
    k = length(hv)
    T = eltype(V)
    onep = T(1.0)
    zerop = T(0.0)
    @views rk = hv[1:k-1]
    pk = zeros(T, size(rk))
    qk = vv
    Qkm = V
    # Orthogonalize
    rk .+= Qkm' * qk
#    qk .-= Qkm * rk
    mul!(qk, Qkm, rk, -1.0, 1.0)
    if orth == "twice"
        # Orthogonalize again
        pk .= Qkm' * qk
#        qk .-= Qkm * pk
        mul!(qk, Qkm, pk, -1.0, 1.0)
        rk .+= pk
    end
    # Keep track of what you did.
    nqk = norm(qk)
    (nqk != 0) || (verbose && (println("breakdown in cgs")))
    (nqk > 0.0) && (qk ./= nqk)
    hv[k] = nqk
end


text:
import Base: abs2, angle, convert, exponent
export BigComplex, abs2, precision, conj, atan2, angle

function show(io::IO, a::BigComplex)
  print(io, a.re, " + ", a.im, "im")
end

function +(a::BigComplex, b::BigComplex)
  return BigComplex(a.re + b.re, a.im + b.im)
end

function -(a::BigComplex, b::BigComplex)
  return BigComplex(a.re - b.re, a.im - b.im)
end

function -(a::BigComplex)
  return BigComplex(-a.re, -a.im)
end

function *(a::BigComplex, b::BigComplex)
  return BigComplex(a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re)
end

function +(a::BigComplex, b::BigFloat)
  return BigComplex(a.re + b, a.im)
end

function *(a::BigComplex, b::BigFloat)
  return BigComplex(a.re * b, a.im * b)
end

function *(a::BigFloat, b::BigComplex)
  return BigComplex(b.re * a, b.im * a)
end

function /(a::BigComplex, b::BigComplex)
  return a*conj(b)/abs2(b)
end

function conj(a::BigComplex)
  return BigComplex(a.re, -a.im)
end

function abs2(a::BigComplex)
  return a.re^2 + a.im^2
end

function angle(a::BigComplex)
  return atan2(imag(a), real(a))
end

function /(a::BigComplex, b::BigFloat)
  return BigComplex(a.re/b, a.im/b)
end

function precision(a::BigComplex)
  return min(Base.precision(a.re), Base.precision(a.im))
end

function real(a::BigComplex)
  return a.re
end

function imag(a::BigComplex)
  return a.im
end


text:
# constrained to one -liners for this one:
println("The answer to Euler Problem 6 is $(sum(1:100)^2 - sum(map(x->x^2, 1:100)))")

max_stars_repo_path: Problem 6.jl

max_stars_repo_name: andrew-christianson/Polyglot-Euler

max_issues_repo_path: Problem 6.jl

max_issues_repo_name: andrew-christianson/Polyglot-Euler

max_forks_repo_path: Problem 6.jl

max_forks_repo_name: andrew-christianson/Polyglot-Euler