text:
using Oceananigans.Grids: total_length, topology

using OffsetArrays: OffsetArray

default_indices(n) = Tuple(Colon() for i=1:n)

#####
##### Creating offset arrays for field data by dispatching on architecture.
#####

"""
Return a range of indices for a field located at either cell `Center`s or `Face`s along a
grid dimension which is `Periodic`, or cell `Center`s for a grid dimension which is `Bounded`.
The dimension has length `N` and `H` halo points.
"""
offset_indices(loc, topo, N, H=0) = 1 - H : N + H

"""
Return a range of indices for a field located at cell `Face`s along a grid dimension which
is `Bounded` and has length `N` and with halo points `H`.
"""
offset_indices(::Type{Face}, ::Type{<:BoundedTopology}, N, H=0) = 1 - H : N + H + 1

"""
Return a range of indices for a field along a 'reduced' dimension.
"""
offset_indices(::Type{Nothing}, topo, N, H=0) = 1 : 1

offset_indices(L, T, N, H, ::Colon) = offset_indices(L, T, N, H)
offset_indices(L, T, N, H, r::UnitRange) = r
offset_indices(::Type{Nothing}, T, N, H, ::UnitRange) = 1:1

function offset_data(underlying_data::AbstractArray, loc, topo, N, H, indices=default_indices(length(loc)))
    ii = offset_indices.(loc, topo, N, H, indices)
    # Add extra indices for arrays of higher dimension than loc, topo, etc.
    extra_ii = Tuple(axes(underlying_data, d) for d in length(ii)+1:ndims(underlying_data))
    return OffsetArray(underlying_data, ii..., extra_ii...)
end

"""
    offset_data(underlying_data, grid::AbstractGrid, loc)

Returns an `OffsetArray` that maps to `underlying_data` in memory,
with offset indices appropriate for the `data` of a field on
a `grid` of `size(grid)` and located at `loc`.
"""
offset_data(underlying_data::AbstractArray, grid::AbstractGrid, loc, indices=default_indices(length(loc))) =
    offset_data(underlying_data, loc, topology(grid), size(grid), halo_size(grid), indices)

"""
    new_data(FT, grid, loc, indices)

Returns an `OffsetArray` of zeros of float type `FT` on `arch`itecture,
with indices corresponding to a field on a `grid` of `size(grid)` and located at `loc`.
"""
function new_data(FT::DataType, grid::AbstractGrid, loc, indices=default_indices(length(loc)))
    arch = architecture(grid)
    Tx, Ty, Tz = total_size(loc, grid, indices)
    underlying_data = zeros(FT, arch, Tx, Ty, Tz)
    return offset_data(underlying_data, grid, loc, indices)
end

new_data(grid, loc, indices=default_indices) = new_data(eltype(grid), grid, loc, indices)


text:
module AlphaShapes
    import BlackBoxOptim.bboptimize, BlackBoxOptim.best_candidate
    import Distances.pairwise, Distances.Euclidean
    import LinearAlgebra.det, LinearAlgebra.inv, LinearAlgebra.norm

    using MiniQhull

    export AlphaShape, AlphaShapeVolume

    """
        CayleyMenger

    Compute the Cayley-Menger matrix from squared
    distances dij^2.

    https://mathworld.wolfram.com/Cayley-MengerDeterminant.html

    For a 2-d triangle the following example would symbolically
    represent the matrix for points x1, x2, and x3 where dij^2
    is the square Euclidean distance between i and j
    Example
```julia-repl
julia>:([0    1      1      1
         1    0      d12^2  d13^2
         1    d21^2  0      d23^2
         1    d32^2  d32^2    0])
:([0 1 1 1; 1 0 d12 ^ 2 d13 ^ 2; 1 d21 ^ 2 0 d23 ^ 2; 1 d32 ^ 2 d32 ^ 2 0])
```
    """
    function CayleyMenger(points::Array{Float64,2})::Array{Float64}
        d = pairwise(Euclidean(),points')
        n = size(points,1)
        CM = ones(n+1,n+1)
        CM[2:end,2:end] = d.^2.0
        CM[end,end] = 0.0
        CM[1,1] = 0.0
        return CM
    end

    """
        SimplexVolume

    Calculate the volume of a simplex using the
    Cayley Menger matrix

    https://westy31.home.xs4all.nl/Circumsphere/ncircumsphere.htm
    """
    function SimplexVolume(points::Array{Float64,2})::Float64

        CM = CayleyMenger(points)
        n = size(CM,1)-2
        return real(sqrt(Complex(((-1.0)^(n+1))/((2.0^n)*factorial(n)^2.0)*det(CM))))
    end

    """
        SimplexCircumSphere

    Find the centre and radius of the circumsphere of a simplex
    https://westy31.home.xs4all.nl/Circumsphere/ncircumsphere.htm
    """
    function SimplexCircumSphere(points::Array{Float64,2})::Tuple{Array{Float64,1},Float64}
        CM = CayleyMenger(points)
        cminv = inv(CM)
        R = real(sqrt(Complex(cminv[1,1]/(-2.0))))
        # convert for barycentric to cartesian coordinates
        λ = cminv[1,2:end]
        c = (λ'*points)[1,:]
        return c,R
    end

    """
        SimplexCircumSphere

    Find the centre and radius of the circumsphere of a simplex
    https://westy31.home.xs4all.nl/Circumsphere/ncircumsphere.htm
    """
    function SimplexCircumRadiusSquared(points::Array{Float64,2})::Float64
        CM = CayleyMenger(points)
        cminv = inv(CM)
        return cminv[1,1]/(-2.0)
    end

    function VertexInTriangle(x::Array{Float64,1},T::Array{Float64,2})::Bool
        if x == T[1,:] || x == T[2,:] || x == T[3,:]
            return true
        else
            return false
        end
    end
    function VertexInTriangulation(x::Array{Float64,1},T::Array{Float64,3})::Bool
        for i in 1:size(T,1)
            if VertexInTriangle(x,T[i,:,:])
                return true
            end
        end
        return false
    end
    function AllPointsInAlphaShape(X::Array{Float64,2},A::Array{Float64,3})::Bool
        for i in 1:size(X,1)
            if VertexInTriangulation(X[i,:],A) == false
                return false
            end
        end
        return true
    end

    """
        AlphaShapeVolume(A::Array{Float64,3})

    return the sum of volumes of all simplices in the alpha shapes A
    """
    function AlphaShapeVolume(A::Array{Float64,3})::Float64
        area = 0.0
        for t in 1:size(A,1)
            area += SimplexVolume(A[t,:,:])
        end
        return area
    end

    """
        GetDelaunayTriangulation(points::Array{Float64,2})::Array{Float64,3}

    Wrap MiniQhull.jl's delaunay to get a delaunay triangualation in any
    dimension
    """
    function GetDelaunayTriangulation(points::Array{Float64,2},indices=false)
        tess = delaunay(permutedims(points,(2,1)))
        Triangles = zeros(size(tess,2),size(tess,1),size(tess,1)-1)
        for i in 1:size(tess,2)
            for j in 1:size(tess,1)
                Triangles[i,j,:] = points[tess[j,i],:]
            end
        end
	if indices
	    return Triangles,tess
	end
        return Triangles
    end

    """
        FindAlpha(X::Array{Float64,2};
            search::Tuple{Float64,Float64}=(0.0, 10.0),
            MaxSteps::Int=100)::Float64

    Use BlackBocOptim to find the optimal alpha value (volume minimiser which cotains all
    points from the input).
    """
    function FindAlpha(X::Array{Float64,2};
        search::Tuple{Float64,Float64}=(0.0, 10.0),
        MaxSteps::Int=100)::Float64
        objective = function(α)
            α = α[1]
            A = AlphaShape(X,α=α);
            t = AllPointsInAlphaShape(X,A)
            o =  AlphaShapeVolume(A)
            # minimise the area but if not all points are
            # included then set to extreme value
            if t
                return o
            else
                return Inf
            end
        end
        res = bboptimize(objective; SearchRange = search, NumDimensions = 1,MaxSteps=MaxSteps);
        return best_candidate(res)[1]
    end

    """
        AlphaShape(X::Array{Float64,2};α::Union{Nothing,Float64}=nothing,
            search::Tuple{Float64,Float64}=(0.0, 10.0),
            MaxSteps::Int=100)::Array{Float64,3}

    Find the alpha shape corresponding to the 2D array of points
    X: [npoints,2], and the α value α.

    If α == nothing then a search over the range of values search is done
    to find the best value. The optimisation objective is the alpha shape
    area (minimise) subject to all points in X being included in the
    alpha shape triangulation.
    """
    function AlphaShape(X::Array{Float64,2};α::Union{Nothing,Float64}=nothing,
        search::Tuple{Float64,Float64}=(0.0, 10.0),
        MaxSteps::Int=100)::Array{Float64,3}
        T = GetDelaunayTriangulation(X)
        if α == nothing
            println("Finding the optimum α value...\n")
            α = FindAlpha(X;search=search,MaxSteps=MaxSteps)
        end
        α2 = α^2.0
        A = [SimplexCircumRadiusSquared(T[i,:,:]) < α2 for i in 1:size(T,1)]
        return T[A.==1,:,:]
    end

	include("utils.jl")
	include("Density.jl")
end # module AlphaShapes

ext: jl

max_stars_repo_name: harveydevereux/CUDA-Whirligigs

max_issues_repo_name: harveydevereux/CUDA-Whirligigs

max_forks_repo_name: harveydevereux/CUDA-Whirligigs


text:
function test_to_vec(x::T; check_inferred=true) where {T}
    check_inferred && @inferred FiniteDifferences.to_vec(x)
    x_vec, back = FiniteDifferences.to_vec(x)
    @test x_vec isa Vector
    @test all(s -> s isa Real, x_vec)
    check_inferred && @inferred back(x_vec)
    @test x == back(x_vec)
    return nothing
end

@testset "chainrulesutils" begin
    y = GBMatrix(sprand(10, 10, 0.5))
    test_to_vec(y)
    v = GBVector(sprand(10, 0.5))
    test_to_vec(v)
end


text:
# Load _g1 from neuron1.jld which is a 3D z-aligned box holding one neuron

function gloss1(par)            # gaussian, fixed a0
    # Minimum: 2.544195 at [1.0279446 1.04854 0.47514126 1.6086882 1.5215918 3.7772746 0.06005877]
    global _g1,_g2
    isdefined(:_g1) || error("Define target _g1")
    isdefined(:_g2) || (_g2 = similar(_g1))
    a0 = 0.03
    (xi,yi,zi,rx,ry,rz,ai)=par
    for x=1:41,y=1:41,z=1:20
        _g2[x,y,z] = a0 + ai * exp(-rx*(xi-x/20)^2 - ry*(yi-y/20)^2 - rz*(zi-z/20)^2)
    end
    vecnorm(broadcast!(-, _g2, _g2, _g1))
end

function gloss2(par)            # gaussian a0 optimized
    # Minimum: 2.318399 at [1.0074097f0,1.0313389f0,0.4666347f0,4.5696244f0,4.1360598f0,11.593927f0,0.07910681f0,0.04300439f0]
    global _g1,_g2
    isdefined(:_g1) || error("Define target _g1")
    isdefined(:_g2) || (_g2 = similar(_g1))
    (xi,yi,zi,rx,ry,rz,ai,a0)=par
    for x=1:41,y=1:41,z=1:20
        _g2[x,y,z] = a0 + ai * exp(-rx*(xi-x/20)^2 - ry*(yi-y/20)^2 - rz*(zi-z/20)^2)
    end
    vecnorm(broadcast!(-, _g2, _g2, _g1))
end

function gloss3(par)            # exponential fixed a0
    # Minimum: 2.554641 at [1.0148009f0,1.0381743f0,0.46588328f0,1.3417006f0,1.2868823f0,1.6828519f0,0.09428701f0]
    global _g1,_g2
    isdefined(:_g1) || error("Define target _g1")
    isdefined(:_g2) || (_g2 = similar(_g1))
    (xi,yi,zi,rx,ry,rz,ai)=par
    a0 = 0.03
    for x=1:41,y=1:41,z=1:20
        _g2[x,y,z] = a0 + ai * exp(-rx*abs(xi-x/20) - ry*abs(yi-y/20) - rz*abs(zi-z/20))
    end
    vecnorm(broadcast!(-, _g2, _g2, _g1))
end

function gloss4(par)            # exponential a0 optimized
    # Minimum: 2.455714 at [1.0124047f0,1.0313879f0,0.4658681f0,2.5158482f0,2.3875716f0,3.5193841f0,0.1388294f0,0.041223526f0]
    global _g1,_g2
    isdefined(:_g1) || error("Define target _g1")
    isdefined(:_g2) || (_g2 = similar(_g1))
    (xi,yi,zi,rx,ry,rz,ai,a0)=par
    for x=1:41,y=1:41,z=1:20
        _g2[x,y,z] = a0 + ai * exp(-rx*abs(xi-x/20) - ry*abs(yi-y/20) - rz*abs(zi-z/20))
    end
    vecnorm(broadcast!(-, _g2, _g2, _g1))
end

function gloss5(par)            # exp(-r^p) with fixed a0
    # Minimum: 2.525750 at [1.020521f0,1.0451639f0,0.46920383f0,1.4722148f0,1.402376f0,2.522549f0,0.07197181f0,1.4944962f0]
    global _g1,_g2
    isdefined(:_g1) || error("Define target _g1")
    isdefined(:_g2) || (_g2 = similar(_g1))
    (xi,yi,zi,rx,ry,rz,ai,p)=par
    a0 = 0.03
    for x=1:41,y=1:41,z=1:20
        _g2[x,y,z] = a0 + ai * exp(-rx*abs(xi-x/20)^p - ry*abs(yi-y/20)^p - rz*abs(zi-z/20)^p)
    end
    vecnorm(broadcast!(-, _g2, _g2, _g1))
end

function gloss6(par)            # exp(-r^p) with a0 optimized
    # Minimum: 2.309404 at [1.0071558f0,1.0304579f0,0.46796566f0,6.200063f0,5.4875193f0,20.799503f0,0.06861037f0,0.04346131f0,2.4839554f0]
    global _g1,_g2
    isdefined(:_g1) || error("Define target _g1")
    isdefined(:_g2) || (_g2 = similar(_g1))
    (xi,yi,zi,rx,ry,rz,ai,a0,p)=par
    for x=1:41,y=1:41,z=1:20
        _g2[x,y,z] = a0 + ai * exp(-rx*abs(xi-x/20)^p - ry*abs(yi-y/20)^p - rz*abs(zi-z/20)^p)
    end
    vecnorm(broadcast!(-, _g2, _g2, _g1))
end

function gloss7(par)            # 1/(1+r^p) with fixed a0
    # Minimum: 2.351727 at [1.0172533f0,1.0372907f0,0.46628773f0,8.777429f0,8.041203f0,20.366737f0,0.11055235f0,1.799353f0]
    global _g1,_g2
    isdefined(:_g1) || error("Define target _g1")
    isdefined(:_g2) || (_g2 = similar(_g1))
    (xi,yi,zi,rx,ry,rz,ai,p)=abs(par)
    a0 = 0.03
    for x=1:41,y=1:41,z=1:20
        _g2[x,y,z] = a0 + ai / (1 + (rx*(xi-x/20)^2 + ry*(yi-y/20)^2 + rz*(zi-z/20)^2)^(p/2))
    end
    vecnorm(broadcast!(-, _g2, _g2, _g1))
end

function gloss8(par)            # 1/(1+r^p) with optimized a0
    # Minimum: 2.318999 at [1.0090305f0,1.0319712f0,0.4678052f0,5.4503393f0,4.9420156f0,13.582409f0,0.07206944f0,0.04085429f0,3.5884995f0]
    global _g1,_g2
    isdefined(:_g1) || error("Define target _g1")
    isdefined(:_g2) || (_g2 = similar(_g1))
    (xi,yi,zi,rx,ry,rz,ai,a0,p)=abs(par)
    for x=1:41,y=1:41,z=1:20
        _g2[x,y,z] = a0 + ai / (1 + (rx*(xi-x/20)^2 + ry*(yi-y/20)^2 + rz*(zi-z/20)^2)^(p/2))
    end
    vecnorm(broadcast!(-, _g2, _g2, _g1))
end

function gloss9(par)            # 1/(1+r^2) with fixed a0
    # Minimum: 2.355760 at [1.0174382f0,1.0378766f0,0.46776757f0,6.482576f0,5.9682374f0,15.343796f0,0.09882537f0]
    global _g1,_g2
    isdefined(:_g1) || error("Define target _g1")
    isdefined(:_g2) || (_g2 = similar(_g1))
    (xi,yi,zi,rx,ry,rz,ai)=abs(par)
    a0 = 0.03; p = 2
    for x=1:41,y=1:41,z=1:20
        _g2[x,y,z] = a0 + ai / (1 + (rx*(xi-x/20)^2 + ry*(yi-y/20)^2 + rz*(zi-z/20)^2)^(p/2))
    end
    vecnorm(broadcast!(-, _g2, _g2, _g1))
end

function gloss10(par)            # 1/(1+r) with fixed a0
    # Minimum: 2.457183 at [1.0198932f0,1.0372633f0,0.46870422f0,559.0978f0,656.0105f0,1288.5681f0,0.42343268f0]
    global _g1,_g2
    isdefined(:_g1) || error("Define target _g1")
    isdefined(:_g2) || (_g2 = similar(_g1))
    (xi,yi,zi,rx,ry,rz,ai)=abs(par)
    a0 = 0.03; p = 1
    for x=1:41,y=1:41,z=1:20
        _g2[x,y,z] = a0 + ai / (1 + (rx*(xi-x/20)^2 + ry*(yi-y/20)^2 + rz*(zi-z/20)^2)^(p/2))
    end
    vecnorm(broadcast!(-, _g2, _g2, _g1))
end

# optimize(gloss, Float32[1.0,1.0,0.5,1.0,1.0,1.0,.1]; show_trace=true, show_every=10, method=LBFGS())

max_stars_repo_name: JuliaTagBot/NeuronTracking

max_issues_repo_name: JuliaTagBot/NeuronTracking

max_forks_repo_name: JuliaTagBot/NeuronTracking


text:
# Point mass configuration ####
function regular_polygon(n; radius = 1, θ_0 = 0., center=[0., 0.])
    function create_point(i)
        θ = 2*π*i/n
        point = [center[1] + radius*cos(θ + θ_0),
                 center[2] + radius*sin(θ + θ_0)]
        return point
    end

    points = Array{Float64,2}(undef, 2, n)
    for i in 1:n
        points[1, i], points[2, i] = create_point(i)
    end

    return points
end

# Spring configuration ####
function complete_graph(n)
    graph = ones(Float64, n, n)

    for i in 1:n
        graph[i, i] = 0.
    end

    return graph
end

function loop_graph(n, connections = [1,])
    graph = zeros(Float64, n, n)

    # Ensure that neighbours that are i away in other direction are connected
    complementary_connections = [n - i for i in connections]
    append!(connections, complementary_connections)

    for i in 1:n, j in (i+1):n
        if (j - i) ∈ connections
            graph[i, j], graph[j, i] = 1., 1.
        end
    end

    return graph
end


text:
export PermutationGroup, perm, circshift!, swap!

swap!(p::Nemo.Generic.perm, i::Integer, j::Integer) = (t = p[i]; p[i] = p[j]; p[j] = t)
circshift!(p::Nemo.Generic.perm, k::Integer) = (t = copy(p.d); circshift!(p.d, t, k))

getindex(x::AbstractArray, p::Nemo.Generic.perm) = x[p.d]


text:
using UnaCorda
using Test

@testset "UnaCorda.jl" begin
    # Write your tests here.
end


text:
NIScope_Horizontal_Dict = Dict()

# Nr of acquisition points
function set_nr_pts(obj::T;npts=1e6) where {T<:NIScope}
	#NISCOPE_ATTR_HORZ_MIN_NUM_PTS=1250009; # ((1000000 + 250000) + 9)
	status = niScope_SetAttributeViInt32(obj, "", NISCOPE_ATTR_HORZ_MIN_NUM_PTS, Int32(npts))
end
# Sampling rate
function set_srate(obj::T;sr=1e6) where {T<:NIScope}
	#NISCOPE_ATTR_MIN_SAMPLE_RATE=1150009; # ((1000000 + 150000) + 9)
	status = niScope_SetAttributeViReal64(obj, "", NISCOPE_ATTR_MIN_SAMPLE_RATE, Float64(sr))
end
# Acquisition duration
function set_hduration(obj::T;d=1e-3) where {T<:NIScope}
	error("not available yet")
end
# Horizontal offset origin
function set_hoffs_div(obj::T;href::Int=1) where {T<:NIScope}
	#NISCOPE_ATTR_HORZ_RECORD_REF_POSITION=1250011; # ((1000000 + 250000) + 11)
	status = niScope_SetAttributeViReal64(obj, "", NISCOPE_ATTR_HORZ_RECORD_REF_POSITION, Float64(href*10))
end


text:
include("1.2.1.1 (a+b x+c x^2)^p.jl")
include("1.2.1.2 (d+e x)^m (a+b x+c x^2)^p.jl")
include("1.2.1.3 (d+e x)^m (f+g x) (a+b x+c x^2)^p.jl")
include("1.2.1.4 (d+e x)^m (f+g x)^n (a+b x+c x^2)^p.jl")
include("1.2.1.5 (a+b x+c x^2)^p (d+e x+f x^2)^q.jl")
include("1.2.1.6 (g+h x)^m (a+b x+c x^2)^p (d+e x+f x^2)^q.jl")
include("1.2.1.7 (a+b x+c x^2)^p (d+e x+f x^2)^q (A+B x+C x^2).jl")
include("1.2.1.8 P(x) (a+b x+c x^2)^p.jl")
include("1.2.1.9 P(x) (d+e x)^m (a+b x+c x^2)^p.jl")

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.1 Quadratic/.jl


text:
immutable FourierLayer <: Mocha.Layer
    bottoms::Vector{Symbol}
    flags::UInt32
    name::AbstractString
    pathkeys::Vector{PathKey}
    timelimit::Float64
    tops::Vector{Symbol}
end

function FourierLayer( ;
        bottoms::Vector{Symbol} = Symbol[],
        flags::UInt32 = FFTW.ESTIMATE,
        name::AbstractString = "fourier",
        pathkeys::Vector{PathKey} = PathKey[],
        timelimit::Float64 = Inf,
        tops::Vector{Symbol} = Symbol[])
    @assert length(bottoms) == 1
    @assert length(tops) == 1
    FourierLayer(bottoms, flags, name, pathkeys, timelimit, tops)
end

Mocha.can_do_bp(::FourierLayer) = false
Mocha.has_neuron(::FourierLayer) = false
Mocha.has_param(::FourierLayer) = false


text:
include("euler/euler.jl")
using .Primes: prime_numbers
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

compute(n::Int)::Int = sum(prime_numbers(n))

compute(10)

compute(2_000_000)

@benchmark compute(2_000_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using NeuralQuantum, QuantumOpticsBase
using LinearAlgebra, SparseArrays
using Test

N = 4
g = 0.7
V = 2.0

hilb = HomogeneousSpin(N)
hilbq = SpinBasis(1//2)^N

ops = []
opsq = []
H = LocalOperator(hilb)
Hq = DenseOperator(hilbq)
for i=1:N
    global H += g/2.0 * sigmax(hilb, i)
    global Hq += g/2.0 * embed(hilbq, i, sigmax(SpinBasis(1//2)))
    push!(ops, sigmam(hilb, i))
    push!(opsq, embed(hilbq, i, sigmam(SpinBasis(1//2))))
end

for i=1:N
    global H  += V/4.0 * sigmaz(hilb, i) * sigmaz(hilb, mod(i, N)+1)
    global Hq += V/4.0 * embed(hilbq, i, sigmaz(SpinBasis(1//2))) * embed(hilbq, mod(i, N)+1, sigmaz(SpinBasis(1//2)))
end

@test Hq.data ≈ NeuralQuantum.to_matrix(H)

@testset "Matrix mapping" begin
    lh = liouvillian(H, [])
    lhq = liouvillian(Hq, [])
    @test NeuralQuantum.to_matrix(lh) ≈ lhq.data

    l0 = liouvillian(ops)
    l0q = liouvillian(DenseOperator(hilbq), opsq)
    @test NeuralQuantum.to_matrix(l0) ≈ l0q.data

    liouv = liouvillian(H, ops)
    liouvq = liouvillian(Hq, opsq)
    @test NeuralQuantum.to_matrix(liouv) ≈ liouvq.data
end


text:
include(joinpath("../src/FingerTrees.jl"))

if VERSION < v"0.7"
    using Base.Test
else
    using Test, Random
end
function updown(E)
    ft = FingerTrees.EmptyFT{Char}()
    for i in 'A':E
        ft = FingerTrees.conjr(ft,i)
        println(ft)
    end
    for i in 'A':E
            k , ft = FingerTrees.splitl(ft)
            println(k, "\t", ft)
    end
end
ft = updown('E')
#ft = updown('Z')

function randomft(N, start = 1, verb = false)
    ft = FingerTrees.EmptyFT{Int}()
    b = bitrand(N)
    l = sum(b) + start - 1
    u = l+1
    i = 1
    for i in 1:N
        if b[i]
            ft = FingerTrees.conjl(l, ft)
            l -= 1
        else
            ft = FingerTrees.conjr(ft,u)
            u += 1
        end
        verb &&     println(ft)

    end
    ft
end

randomft(12, true)

function torture(N, verb=false)
    ft = randomft(N, 1, verb)

    # checks integrity
    for i in 1:N
        @test ft[i] == i
    end

    b = bitrand(N)
    l = 1
    u = N
    for i in 1:N
        if b[i]
            k, ft = FingerTrees.splitl(ft)
            @test k == l
            l += 1
        else
            ft, k = FingerTrees.splitr(ft)
            @test k == u
            u -= 1
        end
verb &&     println(k, " ",i, ft)
    end
    @assert(FingerTrees.isempty(ft))
    n = N
    ft = FingerTrees.concat(randomft(n), randomft(n,n+1))
    FingerTrees.traverse((x,l)->@test(isa(x, Int)), ft)
#=    j = 1
    for i in ft
        @test j==i
        j += 1
    end
=#
    i = rand(1:N)
    println("split $N at $i")
    a, j, b = FingerTrees.split(randomft(N), i)
    for k in 1:i-1
        @test a[k] == k
    end
    @test i==j
    for k in i+1:N
        @test b[k-i] == k
    end

end

torture(3); torture(10); torture(100);
@time for i in 1:100; torture(3); torture(10); torture(100); end
println("done")
ft = FingerTrees.concat(randomft(10), randomft(10,11));


text:
module InformationInequalities
using Combinatorics
using Plots
using LinearAlgebra
using Graphs
using GraphRecipes
using LaTeXStrings
# Write your package code here.

include("utils/utils.jl")
include("entropy/entropies.jl")
include("mutual_information/mi_measures.jl")
include("i_measures/imeasures.jl")
include("visualization/ploting.jl")

export ConditionalEntropyList
export ConditionalMutualInformationList
export ElementalMeasures
export Elemental2Canonical
export Elemental2Canonical_H
export Elemental2Canonical_MI
export InformationExpressionToCanonical
export InformationMeasureToCanonical
export LinearInformationExpressionToCanonical

export entropic_matrix
export find_matrixG
export find_entropic_vector
export unique_entropy_vector
export entropy_vector
export find_entropic_vector
export find_entropic_vector_int

export minimal_EIM_list_canonical
export order_entropic
export order_entropic_expression

export plotEntropyTree
export GeometryConeGamma2
export elementsGamma2

export order_string
export find_subset
export order_entropic1
export numEIM

export plotInformationExpression
export plotIE

export simplify
export simplifyH
export entropic_terms
export minimalE

end


text:
m = [1+1im 2+2im 3+3im;
     4+4im 5+5im 6+6im;
     7+7im 8+8im 9+9im]

v = [1.0,2.0,3.0]

qm = QuArray(m)
qv = QuArray(v)

@assert rawcoeffs(qm*qm) == m*m
@assert rawcoeffs(qm'*qm') == m*m
@assert rawcoeffs(qm'*qm) == m'*m
@assert rawcoeffs(qm*qm') == m*m'

@assert rawcoeffs(qm*qv) == m*v
@assert rawcoeffs(qv'*qm') == m*v
@assert rawcoeffs(qm'*qv) == m'*v
@assert rawcoeffs(qv'*qm) == m'*v

@assert qv'*qv == dot(v,v)
@assert rawcoeffs(qv*qv') == v*v'

@assert qv'*qm*qv == first(v'*m*v)

@assert rawcoeffs(2*im*qv) == 2*im*v
@assert rawcoeffs(qv/2) == v/2
@assert rawcoeffs(scale(2,qv)) == 2*v
@assert coeffs(qm' * -im) == -im * m'
@assert norm(qv) == norm(v)
@test_approx_eq rawcoeffs(normalize(qv)) rawcoeffs(qv)/norm(v)

# a simple test of the `==` operator
@assert qm*3im == scale!(3im, copy(qm))


text:
using Markdown
using MarkdownStyles
using Test

macro renders(x, s)
    esc(quote
        @test sprint(show, "text/html", $x) == $s
    end)
end

macro rendersas(x, y)
    esc(quote
        @renders $x sprint(show, "text/html", $y)
    end)
end
macro rendersas(x, y, class)
    esc(quote
        @renders $x replace(sprint(show, "text/html", $y), "
" => "
")
    end)
end

@testset "MarkdownStyles.jl" begin
    # Types and singletons
    @test Style(:example) isa Style{:example}
    @test Style(:example) === Style(:example)
    @test Style(:example)(md"Hello") isa MarkdownStyles.Styled

    blocks = (
        header=md"# header",
        code=md"""
        ```julia
        println("code")
        ```
        """,
        paragraph=md"paragraph",
        blockquote=md"> blockquote",
        footnote=md"""
        text [^1]

        [^1]: footnote
        """,
        admonition=md"""
        !!! info
            admonition
        """,
        list=md"""- list""",
        ruler=md"""---""",
    )

    inlines = (
        code=md"`code`",
        bold=md"**bold**",
        italic=md"*italic*",
        image=md"![image](https://picsum.photos/200/200)",
        link=md"[link](https://julialang.com)",
    )

    extras = (
        table=md"""
        | a | b | c | d |
        |--:|:--|:-:|---|
        | 0 | 1 | 2 | 3 |
        """,
    )

    # New styles behave the same as Markdown.jl by default
    @testset "NoStyle" begin
        @testset "Block elements" for x ∈ blocks
            @rendersas NoStyle(x) x
        end
        @testset "Inline elements" for x ∈ inlines
            @rendersas NoStyle(x) x
        end
        @testset "Extra elements" for x ∈ extras
            @rendersas NoStyle(x) x
        end
    end

    @testset "Bulma" begin
        @testset "Block elements" for (k, x) ∈ pairs(blocks)
            if k ∉ (:admonition,)
                @rendersas Bulma(x) x "content"
            end
        end
        @testset "Inline elements" for (k, x) ∈ pairs(inlines)
            if k ∉ (:image,)
                @rendersas Bulma(x) x "content"
            end
        end
        @testset "Extra elements" for (k, x) ∈ pairs(extras)
            if k ∉ (:table,)
                @rendersas Bulma(x) x "content"
            end
        end

        @renders Bulma(blocks.admonition) """

Info

admonition
\n

\n
"""
        @renders Bulma(inlines.image) """

\n
"""
        @renders Bulma(extras.table) """

a	b	c	d

0	1	2	3

\n
"""
    end
end


text:
using Dash, DashCoreComponents, DashHtmlComponents, Printf

function LoadExampleCode(filename, wd = nothing)
  example_file_as_string = read(filename, String)

  wd = joinpath(pwd(), isnothing(wd) ? "" : wd )
  return cd(wd) do

    parsed_exprs = ParseExampleExpr(Base.parse_input_line(example_file_as_string; filename = filename))
    return (
        layout = html_div(
            className = "example-container",
            children = eval(parsed_exprs.layout),
            style = Dict("marginBottom" => "10px"),
        ),
        source_code = html_div(
            children = dcc_markdown(
                @sprintf("```julia\n%s\n```", example_file_as_string)
            ),
            className = "code-container",
            style = Dict("borderLeft" => "thin lightgrey solid"),
        ),
        callback! = eval(parsed_exprs.callback!)
    )
  end

end

function ParseExampleExpr(expr)
  layout = nothing
  callback_funcs = Expr(:block)
  for code_block in expr.args
    if code_block isa Expr
        #skip app = ....
        if code_block.head == :(=) && code_block.args[1] == :app
          continue
        end
        #skip run_server
        if code_block.head == :call && code_block.args[1] == :run_server
          continue
        end

        #skip and save layout
        if code_block.head == :(=) && code_block.args[1] == :(app.layout)
          layout = code_block.args[2]
          continue
        end

        #skip and append to list callback
        if is_callback(code_block)
          push!(callback_funcs.args, code_block)
          continue
        end

        eval(code_block)

    end
  end
  return (layout = layout, callback! = Expr(:function, :(app,), callback_funcs))
end

is_callbackcall(expr) = expr.head == :call && expr.args[1] == :callback!

function is_callback(expr)
  is_callbackcall(expr) && return true #callback!((x)->x, ....)
  (expr.head == :do && is_callbackcall(expr.args[1])) && return true #do syntax
  return false
end

function LoadAndDisplayComponent(example_string)
  return html_div(
    (
      html_div(
        children = dcc_markdown(@printf("```Julia\n%s```", example_string)),
        className = "code-container",
        style = Dict("borderLeft" => "thin lightgrey solid")
      ),
      html_div(
        children = include_string(Main, example_string),
        className = "example-container",
        style = Dict("marginBottom" => "10px", "overflow-x" => "initial")
      ),
    )
  )
end

function LoadAndDisplayComponent2(example_string)
  return html_div(
    (
      html_div(
        children = dcc_markdown(@printf("```Julia\n%s```", example_string)),
        className = "code-container",
        style = Dict("borderLeft" => "thin lightgrey solid")
      ),
      html_div(
        children = include_string(Main, example_string),
        className = "example-container",
        style = Dict("marginBottom" => "10px", "padding-bottom" => "30px")
      ),
    )
  )
end

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
function weighted_mean(numbers, weights)
    tmp_weights = []
    tmp_number = []
    for i in 1:length(weights)
        if ismissing(numbers[i])
            continue
        else
            push!(tmp_weights, weights[i])
            push!(tmp_number, numbers[i])
        end
    end
    weights = tmp_weights
    numbers = tmp_number
    if length(weights) == 0
        return 0
    end
    if sum(tmp_weights) == 0
        return 0 # Such a pixel is background noise for us
    end

    wm = 0
    for i in 1:length(tmp_number)
        wm = wm + tmp_number[i]*tmp_weights[i]
    end
    return wm/sum(tmp_weights)
end


text:
using TropicalTensors
using HDF5

function tensor_labels(neighbors)
	k=0
	nbs = copy(neighbors)
	labels = zeros(Int, size(nbs)...)
	for j=1:size(nbs, 2)
		for i=1:size(nbs, 1)
			j2 = nbs[i, j] + 1
			if j2 != 0
				k += 1
				labels[i, j] = k
				nbs[i, j] = -1
				i2 = indexin(j-1, nbs[:,j2])[]
				labels[i2, j2] = k
				nbs[i2, j2] = -1
			end
		end
	end
	@assert count(labels .== 0) == 0 && length(unique(labels)) == length(labels)÷2
	labels
end

function build_tree(order)
	N = size(order, 2)+1
	ids = collect(Any, 1:N)
	for i=1:N-1
		ta, tb = order[:,i]
		ids[ta+1] = ContractionTree(ids[ta+1], ids[tb+1])
		ids[tb+1] = nothing
	end
	ids[1]
end


text:
module Info

export AbstractCTInfo, CTInfo

import ..Utils: yaml_repr

abstract type AbstractCTInfo end

struct CTInfo{NT<:NamedTuple} <: AbstractCTInfo
    value::NT
end

CTInfo(pairs::Pair...) = CTInfo((; pairs...))
CTInfo(; kwargs...) = CTInfo((; kwargs...))

Base.getproperty(ctinfo::CTInfo, p::Symbol) =
    p == :value ? getfield(ctinfo, :value) :
    getproperty(getfield(ctinfo, :value), p)

Base.eltype(::Type{CTInfo{NT}}) where {NT} = eltype(NT)

Base.length(ctinfo::CTInfo) = length(ctinfo.value)
Base.iterate(ctinfo::CTInfo) = iterate(ctinfo.value)
Base.iterate(ctinfo::CTInfo, state) = iterate(ctinfo.value, state)

Base.keys(ctinfo::CTInfo) = keys(ctinfo.value)
Base.values(ctinfo::CTInfo) = values(ctinfo.value)

function Base.merge(a::CTInfo, others::CTInfo...)
    CTInfo(merge(a.value, getproperty.(others, :value)...))
end

yaml_repr(ctinfo::CTInfo) = yaml_repr(ctinfo.value)

end # module


text:
@testset "CCIPCA basic tests" begin
    o = CCIPCA(2, 4)
    @test OnlineStats.outdim(o) == 2
    @test OnlineStats.indim(o)  == 4
    @test size(o)               == (4, 2)
    @test length(o) == lastindex(o) == 2

    # eigen-vectors and eigen-values are 0 before we fit any values:
    @test o[1] == zeros(Float64, 4)
    @test o[2] == zeros(Float64, 4)
    @test OnlineStats.eigenvalue(o, 1) == 0.0
    @test OnlineStats.eigenvalue(o, 2) == 0.0

    # relativevariances is also a vector of zeros before we fit any values:
    @test OnlineStats.relativevariances(o) == zeros(Float64, 2)

    # first vector added goes straight into the projection matrix:
    u1 = rand(4)
    fit!(o, u1)
    @test o[1] == u1/norm(u1)
    @test o[2] == zeros(Float64, 4)

    # We can get eigen-values individually or in array:
    @test length(OnlineStats.eigenvalues(o)) == 2
    @test OnlineStats.eigenvalues(o)[1] == OnlineStats.eigenvalue(o, 1)
    @test OnlineStats.eigenvalues(o)[2] == OnlineStats.eigenvalue(o, 2)

    # errors if asking for eigen-values or vectors outside of outdim range:
    @test_throws AssertionError OnlineStats.eigenvalue(o,  3)
    @test_throws AssertionError OnlineStats.eigenvalue(o,  0)
    @test_throws AssertionError OnlineStats.eigenvalue(o, -1)

    # errors if trying to fit vector of wrong length:
    @test_throws AssertionError fit!(o, rand(3))
    @test_throws AssertionError fit!(o, rand(5))
end # @testset "CCIPCA basic tests" begin

@testset "CCIPCA with no indim to constructor" begin
    o = CCIPCA(3)
    @test OnlineStats.outdim(o) == 3
    @test OnlineStats.indim(o)  == 1
    @test size(o)               == (1, 3)
    @test length(o) == lastindex(o) == 3

    # It will now be set on first fit instead
    fit!(o, rand(5))
    @test OnlineStats.outdim(o) == 3
    @test OnlineStats.indim(o)  == 5
    @test size(o)               == (5, 3)

    # But, obviously, we cannot change the indim again
    @test_throws AssertionError fit!(o, rand(6))
    @test_throws AssertionError fit!(o, rand(4))
end

@testset "Differential test #1 with onlinePCA R package" begin
    # Differential testing with onlinePCA package in R:
    #   lambda <- c(0.0, 0.0)
    #   U <- matrix(rep.int(0, 2*4), nrow = 4, ncol = 2)
    #   library(onlinePCA)
    o = CCIPCA(2, 4)

    #   u1 <- c(1.0, 2.0, 3.0, 4.0)
    u1 = [1.0, 2.0, 3.0, 4.0]
    #   r1 <- ccipca(lambda, U, u1, 0, q=2, l=0);
    fit!(o, u1)
    #   r1$values[1]  # => 5.477226
    @test OnlineStats.eigenvalues(o)[1] ≈ 5.477225575
    ev1 = OnlineStats.eigenvector(o, 1)
    #   r1$vectors[1] # => 0.1825742
    @test isapprox(ev1[1], 0.1825742; atol = 1e-7)
    #   r1$vectors[2] # => 0.3651484
    @test isapprox(ev1[2], 0.3651484; atol = 1e-7)
    #   r1$vectors[3] # => 0.5477226
    @test isapprox(ev1[3], 0.5477226; atol = 1e-7)
    #   r1$vectors[4] # => 0.7302967
    @test isapprox(ev1[4], 0.7302967; atol = 1e-7)
    # Since we have only added one value so far, the 2nd eigenvector is 0
    @test OnlineStats.eigenvector(o, 2) == zeros(Float64, 4)

    #   u2 <- c(3.1, 2.7, 5.6, 3.0)
    u2 = [3.1, 2.7, 5.6, 3.0]
    #   xbar <- (u1 + u2)/2
    #   r2 <- ccipca(r1$values, r1$vectors, u2, 1, q=2, l=0, center = xbar);
    fit!(o, u2)
    #   r2$values[1]  # => 3.011238
    @test isapprox(OnlineStats.eigenvalues(o)[1], 3.011238; atol = 1e-6)
    ev1 = OnlineStats.eigenvector(o, 1)
    #   r2$vectors[1,1] # => 0.2822287
    @test isapprox(ev1[1], 0.2822287; atol = 1e-6)
    #   r2$vectors[2,1] # => 0.3708174
    @test isapprox(ev1[2], 0.3708174; atol = 1e-6)
    #   r1$vectors[3,1] # => 0.6419809
    @test isapprox(ev1[3], 0.6419809; atol = 1e-6)
    #   r1$vectors[4,1] # => 0.6088530
    @test isapprox(ev1[4], 0.6088530; atol = 1e-6)

    #   r2$values[2]  # => 1.500179
    @test isapprox(OnlineStats.eigenvalues(o)[2], 1.500179; atol = 1e-6)
    ev2 = OnlineStats.eigenvector(o, 2)
    #   r2$vectors[1,2] # => 0.520012308
    @test isapprox(ev2[1], 0.520012308; atol = 1e-6)
    #   r2$vectors[2,2] # => -0.003068536
    @test isapprox(ev2[2], -0.003068536; atol = 1e-6)
    #   r1$vectors[3,2] # => 0.457338448
    @test isapprox(ev2[3], 0.457338448; atol = 1e-6)
    #   r1$vectors[4,2] # => -0.721400948
    @test isapprox(ev2[4], -0.721400948; atol = 1e-6)
end # @testset "Differential test with onlinePCA R package" begin

function setup_ccipca_randomly(indimrange = 4:100, maxn = 200;
    scalerange = 0.01:0.001:42.42, fit = true)
    indim = rand(indimrange)
    outdim = rand(1:min(5, indim-1))
    o = CCIPCA(outdim, indim)
    us = Vector{Float64}[]
    for _ in 1:rand(1:maxn)
        u = rand(scalerange, indim)
        push!(us, u)
        if fit
            fit!(o, u)
        end
    end
    return o, us
end

@testset "fit!+transform == fittransform! and callable (random testing)" begin
    function test_f!tr_equals_ft!r(o1, o2, u)
        fit!(o1, u)
        uproj1 = OnlineStats.transform(o1, u)
        urec1  = OnlineStats.reconstruct(o1, uproj1)
        uproj2 = OnlineStats.fittransform!(o2, u)
        urec2 = OnlineStats.reconstruct(o2, uproj2)

        @test uproj1 == uproj2
        @test urec1 == urec2

        # calling the ccipca object is same as transform
        uproj1b = o1(u)
        @test uproj1b == uproj1
    end

    # Random testing with different indims and outdims
    for _ in 1:10
        o1, us = setup_ccipca_randomly(4:100, 30; fit = false)
        o2 = CCIPCA(OnlineStats.outdim(o1), OnlineStats.indim(o1))
        for u in us
            test_f!tr_equals_ft!r(o1, o2, u)
        end
    end
end

@testset "sort! and relativevariances (random testing)" begin
    for _ in 1:10
        o, us = setup_ccipca_randomly(4:42, 25)
        lambdapre = deepcopy(o.lambda)
        Upre = deepcopy(o.U)
        sort!(o)

        # Ensure eigen-values and -vectors are correctly sorted:
        for i in 1:(OnlineStats.outdim(o)-1)
            evi = OnlineStats.eigenvalue(o, i)
            # earlier eigenvalues are larger than equal later ones
            @test evi >= OnlineStats.eigenvalue(o, i+1)
            # current eigenvector is same as where current eigenvalue was previously
            idx = findfirst(ev -> ev == evi, lambdapre)
            @test OnlineStats.eigenvector(o, i) == Upre[:, idx]
        end

        # Ensure relativevariances is also sorted and has valid ranges
        v = OnlineStats.relativevariances(o)
        @test length(v) == OnlineStats.outdim(o)
        @test isapprox(sum(v), 1.0; atol=1e-4)
        for i in 1:(OnlineStats.outdim(o)-1)
            @test v[i] >= v[i+1]
            @test 0.0 <= v[i] <= 1.0
        end
    end
end


text:
# This is a demo script that shows how to use Celeste on a small
# cluster without a shared file system. It serves files over HTTP
# using nginx. It is appropriate for applications where network
# bandwidth significantly exceeds disk bandwith and for scales
# of < 1000 ranks.

# Note: This example assumes that SDSS data is available (in SDSS format) in
# /sdss on the master (the node from which this script is invoked). Different
# settings can be made in the simple_nginx.conf/iosettings.yml file.

usage = """
infer_boxes_http.jl <hostfile> <rcf_nsrcs> <boxfile>
"""

if length(ARGS) != 3
    println(usage)
    exit(-1)
end

hostfile = ARGS[1]
np = length(collect(eachline(hostfile)))

# Add 100 workers and initialize MPI between them
using MPI
manager=MPIManager(np=np, mpirun_cmd=`mpiexec -x JULIA_NUM_THREADS=8 -v -mca plm_rsh_no_tree_spawn 1 -np $np --map-by node --hostfile $hostfile`)
addprocs(manager)

# Load Celeste
using Celeste
@everywhere ENV["CELESTE_RANKS_PER_NODE"]=1
@everywhere ENV["CELESTE_THREADS_PER_CORE"]=1
here = @__DIR__
!isdefined(:CelesteMultiNode) && @everywhere include(joinpath($here,"..","..","src","multinode_run.jl"))
!isdefined(:CelesteHTTPIO) && @everywhere include(joinpath($here,"httpio.jl"))
include(joinpath(@__DIR__,"..","..","bin","binutil.jl"))

# Load Input settings
strategy = Celeste.read_settings_file(joinpath(@__DIR__,"iosettings.yml"))
all_rcfs, all_rcf_nsrcs = parse_rcfs_nsrcs(joinpath(ARGS[2]))
boxes, boxes_rcf_idxs = parse_boxfile(joinpath(ARGS[3]))

# Launch an nginx process for serving large static files
if isfile("/tmp/celeste-nginx.pid")
    run(pipeline(`kill -QUIT $(readstring("/tmp/celeste-nginx.pid")))`, stderr=DevNull))
end
spawn(`nginx -p $(@__DIR__) -c simple_nginx.conf`)

hostname = strip(readstring(`hostname`))

# Tell each worker to run Celeste
futures = map(workers()) do worker
    remotecall(Celeste.ParallelRun.infer_boxes, worker, CelesteMultiNode.DtreeStrategy(),
        all_rcfs, all_rcf_nsrcs, [boxes], [boxes_rcf_idxs],
        CelesteHTTPIO.HTTPStrategy("$hostname:10400", strategy), true, joinpath(ENV["HOME"],"output"))
end

# If we're not at the REPL, wait for completion before exiting
Base.isinteractive() || foreach(fetch, futures)


text:
function variable_dcbranch_current(pm::_PM.AbstractPowerModel; kwargs...)
end

"variable: `vdcm[i]` for `i` in `dcbus`es"
function variable_dcgrid_voltage_magnitude(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    vdcm = _PM.var(pm, nw)[:vdcm] = JuMP.JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :busdc)], base_name="$(nw)_vdcm",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :busdc, i), "Vdc", 1.0)
    )

    if bounded
        for (i, busdc) in _PM.ref(pm, nw, :busdc)
            JuMP.set_lower_bound(vdcm[i],  busdc["Vdcmin"])
            JuMP.set_upper_bound(vdcm[i],  busdc["Vdcmax"])
        end
    end

    report && _IM.sol_component_value(pm, nw, :busdc, :vm, _PM.ids(pm, nw, :busdc), vdcm)
end

function variable_dcgrid_voltage_magnitude(pm::_PM.AbstractLPACModel; nw::Int=pm.cnw, bounded = true, report::Bool=true)
    phivdcm = _PM.var(pm, nw)[:phi_vdcm] = JuMP.JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :busdc)], base_name="$(nw)_phi_vdcm",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :busdc, i), "Vdc")
    )

    if bounded
        for (i, busdc) in _PM.ref(pm, nw, :busdc)
            JuMP.set_lower_bound(phivdcm[i],  busdc["Vdcmin"] - 1)
            JuMP.set_upper_bound(phivdcm[i],  busdc["Vdcmax"] - 1)
        end
    end

    report && _IM.sol_component_value(pm, nw, :busdc, :phivdcm, _PM.ids(pm, nw, :busdc), phivdcm)

end
"variable: `vdcm[i]` for `i` in `dcbus`es"
function variable_dcgrid_voltage_magnitude_sqr(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    wdc = _PM.var(pm, nw)[:wdc] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :busdc)], base_name="$(nw)_wdc",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :busdc, i), "Vdc", 1.0)^2
    )
    wdcr = _PM.var(pm, nw)[:wdcr] = JuMP.@variable(pm.model,
    [(i,j) in _PM.ids(pm, nw, :buspairsdc)], base_name="$(nw)_wdcr",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :busdc, i), "Vdc", 1.0)^2
    )

    if bounded
        for (i, busdc) in _PM.ref(pm, nw, :busdc)
            JuMP.set_lower_bound(wdc[i],  busdc["Vdcmin"]^2)
            JuMP.set_upper_bound(wdc[i],  busdc["Vdcmax"]^2)
        end
        for (bp, buspairdc) in _PM.ref(pm, nw, :buspairsdc)
            JuMP.set_lower_bound(wdcr[bp],  0)
            JuMP.set_upper_bound(wdcr[bp],  buspairdc["vm_fr_max"] * buspairdc["vm_to_max"])
        end
    end

    report && _IM.sol_component_value(pm, nw, :busdc, :wdc, _PM.ids(pm, nw, :busdc), wdc)
end

"variable: `p_dcgrid[l,i,j]` for `(l,i,j)` in `arcs_dcgrid`"
function variable_active_dcbranch_flow(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    p = _PM.var(pm, nw)[:p_dcgrid] = JuMP.@variable(pm.model,
    [(l,i,j) in _PM.ref(pm, nw, :arcs_dcgrid)], base_name="$(nw)_pdcgrid",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :branchdc, l), "p_start", 1.0)
    )

    if bounded
        for arc in _PM.ref(pm, nw, :arcs_dcgrid)
            l,i,j = arc
            JuMP.set_lower_bound(p[arc], -_PM.ref(pm, nw, :branchdc, l)["rateA"])
            JuMP.set_upper_bound(p[arc],  _PM.ref(pm, nw, :branchdc, l)["rateA"])
        end
    end

    report && _IM.sol_component_value_edge(pm, nw, :branchdc, :pf, :pt, _PM.ref(pm, nw, :arcs_dcgrid_from), _PM.ref(pm, nw, :arcs_dcgrid_to), p)
end

"variable: `ccm_dcgrid[l]` for `(l)` in `branchdc`"
function variable_dcbranch_current_sqr(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    vpu = 0.8;
    cc = _PM.var(pm, nw)[:ccm_dcgrid] = JuMP.@variable(pm.model,
    [l in _PM.ids(pm, nw, :branchdc)], base_name="$(nw)_ccm_dcgrid",
    start = (_PM.comp_start_value(_PM.ref(pm, nw, :branchdc, l), "p_start", 0.0) / vpu)^2
    )
    if bounded
        for (l, branchdc) in _PM.ref(pm, nw, :branchdc)
            JuMP.set_lower_bound(cc[l], 0)
            JuMP.set_upper_bound(cc[l], (branchdc["rateA"] / vpu)^2)
        end
    end

    report && _IM.sol_component_value(pm, nw, :branchdc, :ccm, _PM.ids(pm, nw, :branchdc), cc)
end

####################### TNEP variables ################################

function variable_dcbranch_current_ne(pm::_PM.AbstractPowerModel; kwargs...)
end

"variable: `vdcm[i]` for `i` in `dcbus`es"
function variable_dcgrid_voltage_magnitude_ne(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    vdcm_ne = _PM.var(pm, nw)[:vdcm_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :busdc_ne)], base_name="$(nw)_vdcm_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :busdc_ne, i), "Vdc",  1.0)
    )
    if bounded
        for (i, busdc) in _PM.ref(pm, nw, :busdc_ne)
            JuMP.set_lower_bound(vdcm_ne[i],  busdc["Vdcmin"])
            JuMP.set_upper_bound(vdcm_ne[i],  busdc["Vdcmax"])
        end
    end

    report && _IM.sol_component_value(pm, nw, :busdc_ne, :vm, _PM.ids(pm, nw, :busdc_ne), vdcm_ne)
end

function variable_dcgrid_voltage_magnitude_ne(pm::_PM.AbstractLPACModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
        phivdcm_ne = _PM.var(pm, nw)[:phi_vdcm_ne] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :busdc_ne)], base_name="$(nw)_phi_vdcm_ne",
        lower_bound = _PM.ref(pm, nw, :busdc_ne, i, "Vdcmin") - 1,
        upper_bound = _PM.ref(pm, nw, :busdc_ne, i, "Vdcmax") - 1,
        start = _PM.comp_start_value(_PM.ref(pm, nw, :busdc_ne, i), "Vdc")
        )
        if bounded
            for (i, busdc) in _PM.ref(pm, nw, :busdc_ne)
                JuMP.set_lower_bound(phivdcm_ne[i],  busdc["Vdcmin"] - 1)
                JuMP.set_upper_bound(phivdcm_ne[i],  busdc["Vdcmax"] -1 )
            end
        end
        report && _IM.sol_component_value(pm, nw, :busdc_ne, :phivdcm_ne, _PM.ids(pm, nw, :busdc_ne), phivdcm_ne)

#TODO
# think about creating an arc/dict with branchdc_ne pointing to both existing and new buses. Then update limits with corresponding buses.
        phivdcm_fr_ne = _PM.var(pm, nw)[:phi_vdcm_fr] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :branchdc_ne)], base_name="$(nw)_phi_vdcm_fr",
        start = 0
        )
        if bounded
            for (i, branchdc) in _PM.ref(pm, nw, :branchdc_ne)
                JuMP.set_lower_bound(phivdcm_fr_ne[i],  -0.2)
                JuMP.set_upper_bound(phivdcm_fr_ne[i],  0.2 )
            end
        end
        report && _IM.sol_component_value(pm, nw, :branchdc_ne, :phivdcm_fr, _PM.ids(pm, nw, :branchdc_ne), phivdcm_fr_ne)

        phivdcm_to_ne = _PM.var(pm, nw)[:phi_vdcm_to] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :branchdc_ne)], base_name="$(nw)_phi_vdcm_to",
        start = 0
        )

        if bounded
            for (i, branchdc) in _PM.ref(pm, nw, :branchdc_ne)
                JuMP.set_lower_bound(phivdcm_to_ne[i],  -0.2)
                JuMP.set_upper_bound(phivdcm_to_ne[i],  0.2 )
            end
        end
        report && _IM.sol_component_value(pm, nw, :branchdc_ne, :phivdcm_to, _PM.ids(pm, nw, :branchdc_ne), phivdcm_to_ne)
end

"variable: `vdcm[i]` for `i` in `dcbus`es"
function variable_dcgrid_voltage_magnitude_sqr_ne(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    bi_bp = Dict([(i, (b["fbusdc"], b["tbusdc"])) for (i,b) in _PM.ref(pm, nw, :branchdc_ne)])
    bus_vdcmax = merge(Dict([(b,bus["Vdcmax"]) for (b,bus) in _PM.ref(pm, nw, :busdc)]),
    Dict([(b,bus["Vdcmax"]) for (b,bus) in _PM.ref(pm, nw, :busdc_ne)]))
    bus_vdcmin = merge(Dict([(b,bus["Vdcmin"]) for (b,bus) in _PM.ref(pm, nw, :busdc)]),
    Dict([(b,bus["Vdcmin"]) for (b,bus) in _PM.ref(pm, nw, :busdc_ne)]))
         # display(_PM.ids(pm, nw, :buspairsdc_ne))
        wdc_ne = _PM.var(pm, nw)[:wdc_ne] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :busdc_ne)], base_name="$(nw)_wdc_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :busdc_ne, i), "Vdc",  1.0)^2,
        )
        wdcr_ne = _PM.var(pm, nw)[:wdcr_ne] = JuMP.@variable(pm.model,
        [l in _PM.ids(pm, nw, :branchdc_ne)], base_name="$(nw)_wdcr_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :busdc_ne, bi_bp[l][1]), "Vdc",  1.0)^2,
        )
        if bounded
            for (i, busdc) in _PM.ref(pm, nw, :busdc_ne)
                JuMP.set_lower_bound(wdc_ne[i],  busdc["Vdcmin"]^2)
                JuMP.set_upper_bound(wdc_ne[i],  busdc["Vdcmax"]^2)
            end
            for (br, branchdc) in _PM.ref(pm, nw, :branchdc_ne)
                JuMP.set_lower_bound(wdcr_ne[br],  0)
                JuMP.set_upper_bound(wdcr_ne[br],  bus_vdcmax[bi_bp[br][1]] * bus_vdcmax[bi_bp[br][2]])
            end
        end
        report && _IM.sol_component_value(pm, nw, :busdc_ne, :wdc_ne, _PM.ids(pm, nw, :busdc_ne), wdc_ne)
        report && _IM.sol_component_value(pm, nw, :branchdc_ne, :wdcr_ne, _PM.ids(pm, nw, :branchdc_ne), wdcr_ne)
end

function variable_dcgrid_voltage_magnitude_sqr_du(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true) # this has to to every branch, different than its counterpart(Wdc_fr) since two candidate branches can be connected to same node and two duplicate variables will be needed
    bi_bp = Dict([(i, (b["fbusdc"], b["tbusdc"])) for (i,b) in _PM.ref(pm, nw, :branchdc_ne)])
    wdc_fr_ne = _PM.var(pm, nw)[:wdc_du_fr] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :branchdc_ne)], base_name="$(nw)_wdc_du_fr",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :busdc_ne, bi_bp[i][1]), "Vdc",  1.0)^2,
    )
    wdc_to_ne = _PM.var(pm, nw)[:wdc_du_to] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :branchdc_ne)], base_name="$(nw)_wdc_du_to",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :busdc_ne, bi_bp[i][1]), "Vdc",  1.0)^2,
    )
    #TODO replace wdc_du_fr and wdc_du_to with wdc_fr and wdc_to make make it consistent with PM, there multiplication is defined by wr - real and wi- imag
    wdcr_frto_ne = _PM.var(pm, nw)[:wdcr_du] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :branchdc_ne)], base_name="$(nw)_wdcr_du",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :busdc_ne, bi_bp[i][1]), "Vdc",  1.0)^2,
    )

    if bounded
        for (i, branchdc) in _PM.ref(pm, nw, :branchdc_ne)
            JuMP.set_lower_bound(wdc_fr_ne[i],  0)
            JuMP.set_upper_bound(wdc_fr_ne[i],  1.21)
            JuMP.set_lower_bound(wdc_to_ne[i],  0)
            JuMP.set_upper_bound(wdc_to_ne[i],  1.21)
            JuMP.set_lower_bound(wdcr_frto_ne[i],  0)
            JuMP.set_upper_bound(wdcr_frto_ne[i],  1.21)
        end
    end
    report && _IM.sol_component_value(pm, nw, :branchdc_ne, :wdc_du_fr, _PM.ids(pm, nw, :branchdc_ne), wdc_fr_ne)
    report && _IM.sol_component_value(pm, nw, :branchdc_ne, :wdc_du_to, _PM.ids(pm, nw, :branchdc_ne), wdc_to_ne)
    report && _IM.sol_component_value(pm, nw, :branchdc_ne, :wdcr_du, _PM.ids(pm, nw, :branchdc_ne), wdcr_frto_ne)
end

"variable: `p_dcgrid[l,i,j]` for `(l,i,j)` in `arcs_dcgrid`"
function variable_active_dcbranch_flow_ne(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
        p = _PM.var(pm, nw)[:p_dcgrid_ne]= JuMP.@variable(pm.model,
        [(l,i,j) in _PM.ref(pm, nw, :arcs_dcgrid_ne)], base_name="$(nw)_pdcgrid_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :branchdc_ne, l), "p_start",  1.0)
        )

    if bounded
        for arc in _PM.ref(pm, nw, :arcs_dcgrid_ne)
            l,i,j = arc
            JuMP.set_lower_bound(p[arc], -_PM.ref(pm, nw, :branchdc_ne, l)["rateA"])
            JuMP.set_upper_bound(p[arc],  _PM.ref(pm, nw, :branchdc_ne, l)["rateA"])
        end
    end

    report && _IM.sol_component_value_edge(pm, nw, :branchdc_ne, :pf, :pt, _PM.ref(pm, nw, :arcs_dcgrid_from_ne), _PM.ref(pm, nw, :arcs_dcgrid_to_ne), p)
end

"variable: `ccm_dcgrid[l]` for `(l)` in `branchdc`"
function variable_dcbranch_current_sqr_ne(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    vpu = 0.8
    cc= _PM.var(pm, nw)[:ccm_dcgrid_ne] = JuMP.@variable(pm.model,
    [l in _PM.ids(pm, nw, :branchdc_ne)], base_name="$(nw)_ccm_dcgrid_ne",
    start = (_PM.comp_start_value(_PM.ref(pm, nw, :branchdc_ne, l), "p_start",  0.0) / vpu)^2
    )

    if bounded
        for (l, branchdc) in _PM.ref(pm, nw, :branchdc_ne)
            JuMP.set_lower_bound(cc[l], 0)
            JuMP.set_upper_bound(cc[l], (branchdc["rateA"] / vpu)^2)
        end
    end

    report && _IM.sol_component_value(pm, nw, :branchdc_ne, :ccm, _PM.ids(pm, nw, :branchdc), cc)
end

"variable: `0 <= convdc_ne[c] <= 1` for `c` in `candidate converters"
function variable_branch_ne(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, relax::Bool=false, report::Bool=true)
    if !relax
        Z_dc_branch_ne = _PM.var(pm, nw)[:branchdc_ne] = JuMP.@variable(pm.model, #branch_ne is also name in PowerModels, branchdc_ne is candidate branches
        [l in _PM.ids(pm, nw, :branchdc_ne)], base_name="$(nw)_branch_ne",
        binary = true,
        start = _PM.comp_start_value(_PM.ref(pm, nw, :branchdc_ne, l), "convdc_tnep_start",  0.0)
        )
    else
        Z_dc_branch_ne = _PM.var(pm, nw)[:branchdc_ne] = JuMP.@variable(pm.model, #branch_ne is also name in PowerModels, branchdc_ne is candidate branches
        [l in _PM.ids(pm, nw, :branchdc_ne)], base_name="$(nw)_branch_ne",
        lower_bound = 0,
        upper_bound = 1,
        start = _PM.comp_start_value(_PM.ref(pm, nw, :branchdc_ne, l), "convdc_tnep_start",  0.0)
        )
    end
    report && _IM.sol_component_value(pm, nw, :branchdc_ne, :isbuilt, _PM.ids(pm, nw, :branchdc_ne), Z_dc_branch_ne)
end


text:
###################################################################################################
#                                        Production Conditions
###################################################################################################
function can_attend(actr)
    c1(actr) = !isempty(actr.visual_location.buffer)
    c2(actr) = !actr.visual.state.busy
    return c1, c2
end

function can_encode(actr)
    c1(actr) = !actr.visual.state.empty
    c2(actr) = !actr.imaginal.state.busy
    c3(actr) = actr.imaginal.state.empty
    c4(actr) = actr.visual.buffer[1].slots.text != "start test"
    c5(actr) = actr.goal.buffer[1].slots.goal != :test
    return c1, c2, c3, c4, c5
end

function can_start(actr)
    c1(actr) = !actr.visual.state.empty
    c2(actr) = actr.visual.buffer[1].slots.text == "start test"
    return c1, c2
end

function can_retrieve(actr)
    c1(actr) = !actr.declarative.state.busy
    c2(actr) = !actr.visual.state.empty
    c3(actr) = actr.goal.buffer[1].slots.goal == :test
    return c1, c2, c3
end

function can_respond_yes(actr)
    c1(actr) = !actr.declarative.state.empty
    c2(actr) = !actr.motor.state.busy
    return c1, c2
end

function can_respond_no(actr)
    c1(actr) = actr.declarative.state.error
    c2(actr) = !actr.motor.state.busy
    return c1, c2
end
###################################################################################################
#                                        Production Actions
###################################################################################################
function encode_action(actr, task)
    buffer = actr.visual.buffer
    word = buffer[1].slots.text
    chunk = Chunk(;word)
    clear_buffer!(actr.visual)
    # modify chunk
    tΔ =  encoding!(actr, chunk)
    register!(actr, clear_buffer!, after, tΔ , actr, actr.imaginal)
    return nothing
end

function attend_action(actr, args...)
    buffer = actr.visual_location.buffer
    chunk = deepcopy(buffer[1])
    clear_buffer!(actr.visual_location)
    attending!(actr, chunk)
    return nothing
end

function start_action(actr, task)
    clear_buffer!(actr.visual)
    actr.goal.buffer[1].slots = (goal=:test,)
    register!(actr, ()->(), now; description="start test phase")
    return nothing
end

function retrieve_word(actr, task)
    word = actr.visual.buffer[1].slots.text
    clear_buffer!(actr.visual)
    retrieving!(actr; word)
end

function respond_yes(actr, task)
    clear_buffer!(actr.declarative)
    key = "y"
    responding!(actr, task, key)
    return nothing
end

function respond_no(actr, task)
    clear_buffer!(actr.declarative)
    key = "n"
    responding!(actr, task, key)
    return nothing
end


text:
# ---
# title: 1163. Last Substring in Lexicographical Order
# id: problem1163
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Suffix Array
# link: <https://leetcode.com/problems/last-substring-in-lexicographical-order/description/>
# hidden: true
# ---
#
# Given a string `s`, return the last substring of `s` in lexicographical order.
#
#
#
# **Example 1:**
#
#
#
#     Input: "abab"
#     Output: "bab"
#     Explanation: The substrings are ["a", "ab", "aba", "abab", "b", "ba", "bab"]. The lexicographically maximum substring is "bab".
#
#
# **Example 2:**
#
#
#
#     Input: "leetcode"
#     Output: "tcode"
#
#
#
#
# **Note:**
#
#   1. `1 <= s.length <= 4 * 10^5`
#   2. s contains only lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
struct Checker{is_supervised,
                  supports_weights,
                  input_scitype,
                    target_scitype} end

function Base.getproperty(::Checker{is_supervised,
                                    supports_weights,
                                    input_scitype,
                                    target_scitype},
                          field::Symbol) where {is_supervised,
                                                supports_weights,
                                                input_scitype,
                                                target_scitype}
    if field === :is_supervised
        return is_supervised
    elseif field === :supports_weights
        return supports_weights
    elseif field === :input_scitype
        return input_scitype
    elseif field === :target_scitype
        return target_scitype
    else
        throw(ArgumentError("Unsupported property. "))
    end
end

Base.propertynames(::Checker) =
    (:is_supervised, :supports_weights, :input_scitype, :target_scitype)

function _as_named_tuple(s::Checker)
    names = propertynames(s)
    NamedTuple{names}(Tuple(getproperty(s, p) for p in names))
end

# function Base.show(io::IO, ::MIME"text/plain", S::Checker)
#     show(io, MIME("text/plain"), _as_named_tuple(S))
# end

matching(X)       = Checker{false,false,scitype(X),missing}()
matching(X, y)    = Checker{true,false,scitype(X),scitype(y)}()
matching(X, y, w) = Checker{true,true,scitype(X),scitype(y)}()


text:
#! format: off

abstract type AbstractRenewableFormulation <: AbstractDeviceFormulation end
abstract type AbstractRenewableDispatchFormulation <: AbstractRenewableFormulation end
struct RenewableFullDispatch <: AbstractRenewableDispatchFormulation end
struct RenewableConstantPowerFactor <: AbstractRenewableDispatchFormulation end

get_variable_sign(_, ::Type{<:PSY.RenewableGen}, ::AbstractRenewableFormulation) = 1.0
########################### ActivePowerVariable, RenewableGen #################################

get_variable_binary(::ActivePowerVariable, ::Type{<:PSY.RenewableGen}, ::AbstractRenewableFormulation) = false

get_variable_expression_name(::ActivePowerVariable, ::Type{<:PSY.RenewableGen}) = :nodal_balance_active

get_variable_lower_bound(::ActivePowerVariable, d::PSY.RenewableGen, ::AbstractRenewableFormulation) = 0.0
get_variable_upper_bound(::ActivePowerVariable, d::PSY.RenewableGen, ::AbstractRenewableFormulation) = PSY.get_max_active_power(d)

########################### ReactivePowerVariable, RenewableGen #################################

get_variable_binary(::ReactivePowerVariable, ::Type{<:PSY.RenewableGen}, ::AbstractRenewableFormulation) = false

get_variable_expression_name(::ReactivePowerVariable, ::Type{<:PSY.RenewableGen}) = :nodal_balance_reactive

#! format: on

####################################### Reactive Power constraint_infos #########################
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractDeviceFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.RenewableGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ReactivePowerVariable,
                T,
            ),
            variable_name = make_variable_name(ReactivePowerVariable, T),
            limits_func = x -> PSY.get_reactive_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    ::Type{T},
    ::Type{<:RenewableConstantPowerFactor},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.RenewableGen}
    return DeviceRangeConstraintSpec(;
        custom_optimization_container_func = custom_reactive_power_constraints!,
    )
end

function custom_reactive_power_constraints!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::Type{RenewableConstantPowerFactor},
) where {T <: PSY.RenewableGen}
    names = [PSY.get_name(d) for d in devices]
    time_steps = model_time_steps(optimization_container)
    p_var = get_variable(optimization_container, ACTIVE_POWER, T)
    q_var = get_variable(optimization_container, REACTIVE_POWER, T)
    constraint_val = JuMPConstraintArray(undef, names, time_steps)
    assign_constraint!(optimization_container, REACTIVE_RANGE, T, constraint_val)
    for t in time_steps, d in devices
        name = PSY.get_name(d)
        pf = sin(acos(PSY.get_power_factor(d)))
        constraint_val[name, t] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            q_var[name, t] == p_var[name, t] * pf
        )
    end
    return
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractRenewableDispatchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.RenewableGen}
    if !use_parameters && !use_forecasts
        return DeviceRangeConstraintSpec(;
            range_constraint_spec = RangeConstraintSpec(;
                constraint_name = make_constraint_name(
                    RangeConstraint,
                    ActivePowerVariable,
                    T,
                ),
                variable_name = make_variable_name(ActivePowerVariable, T),
                limits_func = x -> (min = 0.0, max = PSY.get_active_power(x)),
                constraint_func = device_range!,
                constraint_struct = DeviceRangeConstraintInfo,
            ),
        )
    end

    return DeviceRangeConstraintSpec(;
        timeseries_range_constraint_spec = TimeSeriesConstraintSpec(;
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            parameter_name = use_parameters ? ACTIVE_POWER : nothing,
            forecast_label = "max_active_power",
            multiplier_func = x -> PSY.get_max_active_power(x),
            constraint_func = use_parameters ? device_timeseries_param_ub! :
                              device_timeseries_ub!,
        ),
    )
end

########################## Addition to the nodal balances ##################################
function NodalExpressionSpec(
    ::Type{T},
    ::Type{<:PM.AbstractPowerModel},
    use_forecasts::Bool,
) where {T <: PSY.RenewableGen}
    return NodalExpressionSpec(
        "max_active_power",
        REACTIVE_POWER,
        use_forecasts ? x -> PSY.get_max_reactive_power(x) : x -> PSY.get_reactive_power(x),
        1.0,
        T,
    )
end

function NodalExpressionSpec(
    ::Type{T},
    ::Type{<:PM.AbstractActivePowerModel},
    use_forecasts::Bool,
) where {T <: PSY.RenewableGen}
    return NodalExpressionSpec(
        "max_active_power",
        ACTIVE_POWER,
        use_forecasts ? x -> PSY.get_max_active_power(x) : x -> PSY.get_active_power(x),
        1.0,
        T,
    )
end

##################################### renewable generation cost ############################
function AddCostSpec(
    ::Type{T},
    ::Type{U},
    ::OptimizationContainer,
) where {T <: PSY.RenewableDispatch, U <: AbstractRenewableDispatchFormulation}
    # TODO: remove once cost_function is required
    cost_function = x -> (x === nothing ? 1.0 : PSY.get_variable(x))
    return AddCostSpec(;
        variable_type = ActivePowerVariable,
        component_type = T,
        variable_cost = cost_function,
        multiplier = OBJECTIVE_FUNCTION_NEGATIVE,
    )
end


text:
module Noise
using PoissonRandom
using Random
using ImageCore

"""
    complex_copy(x)

If `x` is real, it returns `x + im * x`
If `x` is complex, it returns `x`,
"""
complex_copy(a::T) where T = (a + 1im * a)::Complex{T}
complex_copy(a::Complex{T}) where T= a::Complex{T}

 # clipping a single value
function clip_v(x)
    return max(0, min(1, x))
end

function clip_v(x, minv, maxv)
    return max(minv, min(maxv, x))
end

 # get the maximum single value of a RGB image
function max_rgb(X::AbstractArray{RGB{T}}) where T
    max_v = -Inf
    for i in eachindex(X)
        max_v = max(red(X[i]), max_v)
        max_v = max(green(X[i]), max_v)
        max_v = max(blue(X[i]), max_v)
    end
    return max_v
end

function apply_noise!(pixel_f, noise_f, X::Union{AbstractArray{Gray{T}}, AbstractArray{RGB{T}},
            AbstractArray{T}}, clip) where T

    f(x) = pixel_f(x, noise_f(x))
     # defining the core functions inside here, gives 6x speed improvement
     # core functions non clipping
    function core_f(x::RGB)
        return RGB(f(red(x)),
                   f(green(x)),
                   f(blue(x)))
    end

    function core_f(x::Gray)
        return Gray(f(gray(x)))
    end

    function core_f(x)
        return f(x)
    end

    function core_f_clip(x::RGB)
        return RGB(clip_v(f(red(x))),
                   clip_v(f(green(x))),
                   clip_v(f(blue(x))))
    end

    function core_f_clip(x::Gray)
        return Gray(clip_v(f(gray(x))))
    end

    function core_f_clip(x)
        return clip_v(f(x))
    end
    # if normed clip the values to be in [0, 1]
    if T <: Normed || clip
        @inbounds for i in eachindex(X)
            X[i] = core_f_clip(X[i])
        end
    else
        @inbounds for i in eachindex(X)
            X[i] = core_f(X[i])
        end
    end
    return X
end

 # function which use exactly the same noise for each color channel of a pixel
function apply_noise_chn!(pixel_f, noise_f, X::AbstractArray{RGB{T}}, clip) where T
    if T <: Normed || clip
        @inbounds for i in eachindex(X)
            a = X[i]
            n = noise_f(red(a)*0)
            X[i] = RGB(clip_v(pixel_f(red(a), n)),
                       clip_v(pixel_f(green(a), n)),
                       clip_v(pixel_f(blue(a), n)))
        end
    else
        @inbounds for i in eachindex(X)
            a = X[i]
            n = noise_f(red(a)*0)
            X[i] = RGB(pixel_f(red(a), n),
                       pixel_f(green(a), n),
                       pixel_f(blue(a), n))
        end
    end
    return X
end

include("poisson.jl")
include("salt_pepper.jl")
include("white_noise_additive.jl")
include("quantization.jl")
include("multiplicative_noise.jl")

end

max_issues_repo_name: roflmaostc/Noise

max_forks_repo_name: roflmaostc/Noise


text:
"""
    max_compression_overhead()

The maximum overhead during compression in bytes.
"""
max_compression_overhead() = Lib.BLOSC_MAX_OVERHEAD

"""
    max_compressed_size(src)

The maximum size after compression in bytes.
"""
max_compressed_size(src) = sizeof(src) + max_compression_overhead()

"""
    max_compress_buffer(src)

Creates a buffer that is guaranteed to contain compressed data
"""
make_compress_buffer(src) = Vector{UInt8}(undef, max_compressed_size(src))

"""
    unsafe_compress!([ctx = CompressionContext(T; kwargs...)], dest::Ptr{T}, dest_size, src::Ptr{T}, src_size; kwargs...)

Compress `src_size` bytes from `src` buffer, into `dest` buffer (with `dest_size` bytes size) using context `ctx` with no checks, return size of the data written to dest

The `unsafe` prefix on this function indicates that no validation is performed on the pointers dest and src to ensure that they are valid. Incorrect usage may corrupt or segfault your program, in the same manner as C.
"""
function unsafe_compress!(ctx::CompressionContext, dest::Ptr{UInt8}, dest_size, src::Ptr{T}, src_size) where {T}
    return Lib.blosc2_compress_ctx(ctx.context, src, src_size, dest, dest_size)
end

unsafe_compress!(dest::Ptr{UInt8}, dest_size, src::Ptr{T}, src_size; kwargs...) where {T} =
    unsafe_compress!(CompressionContext(T; kwargs...), dest, dest_size, src, src_size)

"""
    compress!([ctx = CompressionContext(T; kwargs...)],
         dest::Vector{UInt8}, src::Vector{T},
         [dest_offset = 1], [src_range = 1:length(src)]
        ; kwargs...
        )

Compress elements from  `src_range` of `src` into `dest` buffer starting from `dest_offset` (1-indexed) using context `ctx`
Return size of the data written to dest
"""
function compress!(ctx::CompressionContext, dest::Vector{UInt8}, src::Vector{T}, dest_offset = 1,  src_range::UnitRange{<:Integer} = 1:length(src)) where {T}
    @boundscheck begin
        checkbounds(dest, dest_offset)
        checkbounds(src, src_range)
    end
    dest_size = length(dest) - dest_offset + 1
    src_size = length(src_range) * sizeof(T)
    result = GC.@preserve dest src begin
        unsafe_compress!(ctx, pointer(dest, dest_offset), dest_size, pointer(src, first(src_range)), src_size)
    end
    result < 0 && error("Compression error. Error code: $result")
    return result
end

compress!(dest::Vector{UInt8}, src::Vector{T}, dest_offset = 1,  src_range::UnitRange{<:Integer} = 1:length(src); kwargs...) where {T} =
    compress!(CompressionContext(T; kwargs...), dest, src, dest_offset, src_range)

"""
    compress([ctx = CompressionContext(T; kwargs...)], src::Array{T}, [src_range = 1:length(src)]; kwargs...)::Vector{UInt8}

Return `Vector{UInt8}` consisting of `src_range` of `src` in compressed form using context `ctx`
"""
function compress(ctx::CompressionContext, src::Array{T}, src_range::UnitRange{<:Integer} = 1:length(src)) where {T}
    buffer = make_compress_buffer(src)
    sz = compress!(ctx, buffer, src, 1, src_range)
    resize!(buffer, sz)
    return buffer
end
compress(src::Array{T}, src_range::UnitRange{<:Integer} = 1:length(src); kwargs...) where {T} = compress(CompressionContext(T; kwargs...), src, src_range)


text:
for F in (:acos, :acosd, :acosh, :acot, :acotd, :acoth, :acsc, :acscd, :acsch,
          :asec, :asecd, :asech, :asin, :asind, :asinh, :atan, :atan2, :atand,
          :atanh, :beta, :cbrt, :cos,:cosc, :cosd, :cosh, :cospi, :cot, :cotd,
          :coth, :csc, :cscd, :csch, :deg2rad, :exp, :exp10, :exp2, :expm1,
          :exponent, :frexp, :gamma, :lfact, :lgamma, :log, :log10, :log1p, :log2,
          :mod2pi, :modf, :rad2deg, :rem2pi, :sec, :secd, :sech, :significand,
          :sin, :sinc, :sincos, :sind, :sinh, :sinpi, :sqrt, :tan, :tand, :tanh)
    @eval begin
        $F(x::T) where {T<:MarkableFloat64} = T($F(Float64(x)))
        $F(x::T) where {T<:MarkableFloat32} = T($F(Float32(x)))
        $F(x::T) where {T<:MarkableFloat16} = T($F(Float16(x)))
    end
end


text:
handle("editor.julia.hints") do editor, data
  cur = cursor(data["cursor"]["line"], data["cursor"]["col"])
  code = data["code"]
  mod = Jewel.getmodule(code, cur, filemod = data["module"])
  completions = allcompletions(code, cur, mod = mod, file = data["path"])
  if completions == nothing
    raise(editor, "editor.julia.hints.update",
          @d(:hints => c()))
  else
    raise(editor, "editor.julia.hints.update",
          @d(:hints => completions[:hints],
             :pattern => completions[:pattern].pattern))
  end
end


text:
using RDatasets, DataArrays, Gadfly

set_default_plot_size(6inch, 3inch)

plot(dataset("car", "Womenlf"), x=:HIncome, y=:Region, Geom.histogram2d)


text:
E = SNN.HH(;N = 1)
E.I = [0.001]

SNN.monitor(E, [:v])
SNN.sim!([E], []; dt = 0.01ms, duration = 100ms)


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

import Base:*
"""
  `*(A::WFST,B::WFST)`

Returns `C`, the concatenation/product of `A` and `B`. If `A` converts the sequence `a_i` to `a_o` with weight `wa` and `B` converts the sequence `b_i` to `b_o` with weight `wb` then C` converts the sequence `[a_i,b_i]` to `[a_o,b_o]` with weight `wa*wb`.

```julia
julia> sym = ["a","b","c","d","α","β","γ","δ"];

julia> A = WFST(sym);

julia> add_arc!(A,1=>2,"a"=>"α");

julia> add_arc!(A,2=>3,"b"=>"β");

julia> initial!(A,1); final!(A,3,5)
WFST #states: 3, #arcs: 2, #isym: 8, #osym: 8
|1/0.0f0|
a:α/0.0f0 → (2)
(2)
b:β/0.0f0 → (3)
((3/5.0f0))

julia> B = WFST(sym);

julia> add_arc!(B,1=>2,"c"=>"γ");

julia> add_arc!(B,2=>3,"d"=>"δ");

julia> initial!(B,1); final!(B,3,2)
WFST #states: 3, #arcs: 2, #isym: 8, #osym: 8
|1/0.0f0|
c:γ/0.0f0 → (2)
(2)
d:δ/0.0f0 → (3)
((3/2.0f0))

julia> C = A*B
WFST #states: 6, #arcs: 5, #isym: 8, #osym: 8
|1/0.0f0|
a:α/0.0f0 → (2)
(2)
b:β/0.0f0 → (3)
(3)
ϵ:ϵ/5.0f0 → (4)
(4)
c:γ/0.0f0 → (5)
(5)
d:δ/0.0f0 → (6)
((6/2.0f0))

julia> A(["a","b"])
(["α", "β"], 5.0f0)

julia> B(["c","d"])
(["γ", "δ"], 2.0f0)

julia> C(["a","b","c","d"])
(["α", "β", "γ", "δ"], 7.0f0)

```
"""
function *(A::WFST,B::WFST)
  isym = get_isym(A)
  osym = get_osym(A)
  W = typeofweight(A)
  if isym != get_isym(B)
    throw(ErrorException("WFSTs have different input symbol tables"))
  end
  if osym != get_osym(B)
    throw(ErrorException("WFSTs have different output symbol tables"))
  end
  if W != typeofweight(B)
    throw(ErrorException("WFSTs have different type of weigths $W and $(typeofweight(B))"))
  end
  C = deepcopy(A)
  offset = length(C)
  add_states!(C,length(B))
  for (q,arcs) in enumerate(get_states(B))
    for arc in arcs
      ilab, olab = get_ilabel(arc), get_olabel(arc)
      w, n = get_weight(arc), get_nextstate(arc)
      push!(C[q+offset],Arc(ilab,olab,w,n+offset))
    end
  end
  final_A = get_final(A; single=false)
  initial_B = get_initial(B; single=false)
  for (qf,wf) in final_A
    for (qi,wi) in initial_B
      push!(C[qf],Arc(0,0,wf*wi,qi+offset))
    end
    rm_final!(C,qf)
  end
  for (q,w) in get_final(B; single=false)
    final!(C,q+offset,w)
  end
  return C
end


text:
using Assets, Currencies, Instruments, FixedPointDecimals
# using Currencies: currency, symbol, unit, code, name

@cash USD, EUR, JPY, JOD, CNY

using Test

# Check that some basic currencies have been loaded correctly
# (this part should really be in Currencies.jl, not here)
currencies = ((USD, :USD, 2, 840, "US Dollar"),
              (EUR, :EUR, 2, 978, "Euro"),
              (JPY, :JPY, 0, 392, "Yen"),
              (JOD, :JOD, 3, 400, "Jordanian Dinar"),
              (CNY, :CNY, 2, 156, "Yuan Renminbi"))

@testset "Basic currencies" begin
    for (pos, s, u, c, n) in currencies
        ccy = currency(pos)
        @test symbol(ccy) == s
        @test unit(ccy) == u
        @test name(ccy) == n
        @test code(ccy) == c
    end
end

@testset "All currencies" begin
    for sym in Currencies.allsymbols()
        ccy = Currency{sym}
        ct = cash(sym)
        @test ct == cash(ccy)
        @test ct == cash(sym)
        @test currency(ct) == ccy
        @test symbol(ct) == symbol(ccy)
        @test unit(ct) == unit(ccy)
        @test code(ct) == code(ccy)
        @test name(ct) == name(ccy)

        pos = Position(ct,1)
        pt = typeof(pos)
        @test currency(pt) == currency(ct)
        @test currency(pt) == ccy
        @test 1pt == pos
        @test pos * 1 == pos
        @test 1pos + 1pos == Position(ct,2)
        @test 1pos - 1pos == Position(ct,0)
        @test 20pos / 4pos == FixedDecimal{Int,unit(ct)}(5)
        @test 20pos / 4 == Position(ct,5)
    end
end


text:
"""
`run(doc::WeaveDoc; doctype = :auto, plotlib="Gadfly",
        out_path=:doc, fig_path = "figures", fig_ext = nothing,
        cache_path = "cache", cache = :off)`

Run code chunks and capture output from parsed document.

* `doctype`: :auto = set based on file extension or specify one of the supported formats.
  See `list_out_formats()`
* `plotlib`: `"PyPlot"`, `"Gadfly"`, or `"Winston"`
* `out_path`: Path where the output is generated. Can be: `:doc`: Path of the source document, `:pwd`: Julia working directory,
  `"somepath"`: Path as a AbstractString e.g `"/home/mpastell/weaveout"`
* `fig_path`: where figures will be generated, relative to out_path
* `fig_ext`: Extension for saved figures e.g. `".pdf"`, `".png"`. Default setting depends on `doctype`.
* `cache_path`: where of cached output will be saved.
* `cache`: controls caching of code: `:off` = no caching, `:all` = cache everything,
 `:user` = cache based on chunk options, `:refresh`, run all code chunks and save new cache.

**Note:** Run command from terminal and not using IJulia, Juno or ESS, they tend to mess with capturing output.
"""
function Base.run(doc::WeaveDoc; doctype = :auto, plotlib=:auto,
        out_path=:doc, fig_path = "figures", fig_ext = nothing,
        cache_path = "cache", cache = :off)
    #cache :all, :user, :off, :refresh

    doc.cwd = get_cwd(doc, out_path)
    doctype == :auto && (doctype = detect_doctype(doc.source))
    doc.doctype = doctype
    doc.format = formats[doctype]

    set_rc_params(doc.format.formatdict, fig_path, fig_ext)

    #New sandbox for each document
    sandbox = "ReportSandBox$(rcParams[:doc_number])"
    eval(parse("module $sandbox\nend"))
    SandBox = eval(parse(sandbox))
    rcParams[:doc_number] += 1

    if haskey(doc.format.formatdict, :mimetypes)
      mimetypes = doc.format.formatdict[:mimetypes]
    else
      mimetypes = default_mime_types
    end

    #Reset plotting
    rcParams[:plotlib_set] = false
    plotlib == :auto || init_plotting(plotlib)

    report = Report(doc.cwd, doc.basename, doc.format.formatdict, mimetypes)
    pushdisplay(report)

    if cache != :off && cache != :refresh
        cached = read_cache(doc, cache_path)
        cached == nothing && info("No cached results found, running code")
    else
        cached = nothing
    end

    executed = Any[]
    n = length(doc.chunks)

    for i = 1:n
        chunk = doc.chunks[i]

        if typeof(chunk) == CodeChunk
            options = merge(rcParams[:chunk_defaults], chunk.options)
            merge!(chunk.options, options)

        end

        restore = (cache ==:user && typeof(chunk) == CodeChunk && chunk.options[:cache])

        if cached != nothing && (cache == :all || restore)
            result_chunks = restore_chunk(chunk, cached)
        else

        result_chunks = run_chunk(chunk, report, SandBox)
        end

        executed = [executed; result_chunks]
    end

    popdisplay(report)

    #Clear variables from used sandbox
    clear_sandbox(SandBox)
    doc.chunks = executed

    if cache != :off
        write_cache(doc, cache_path)
    end

    return doc
end

"""Detect the output format based on file extension"""
function detect_doctype(source::AbstractString)
  ext = lowercase(splitext(source)[2])
  ext == ".jl" && return "md2html"
  contains(ext, "md") && return "md2html"
  contains(ext, "rst") && return "rst"
  contains(ext, "tex") && return "texminted"
  contains(ext, "txt") && return "asciidoc"

  return "pandoc"
end

function run_chunk(chunk::CodeChunk, report::Report, SandBox::Module)
    result_chunk = eval_chunk(chunk, report, SandBox)
end

function run_chunk(chunk::DocChunk, report::Report, SandBox::Module)
    return chunk
end

function reset_report(report::Report)
    report.cur_result = ""
    report.figures = AbstractString[]
    report.term_state = :text
end

function run_code(chunk::CodeChunk, report::Report, SandBox::Module)
    expressions = parse_input(chunk.content)
    N = length(expressions)
    #@show expressions
    result_no = 1
    results = ChunkOutput[ ]

    for (str_expr, expr) = expressions
        reset_report(report)
        lastline = (result_no == N)
        rcParams[:plotlib_set] || detect_plotlib(chunk) #Try to autodetect plotting library
        (obj, out) = capture_output(expr, SandBox, chunk.options[:term],
                      chunk.options[:display], rcParams[:plotlib], lastline)
        figures = report.figures #Captured figures
        result = ChunkOutput(str_expr, out, report.cur_result, report.rich_output, figures)
        report.rich_output = ""
        push!(results, result)
        result_no += 1
    end

    #Save figures only in the end of chunk for PyPlot
    if rcParams[:plotlib] == "PyPlot"
        savefigs_pyplot(report::Report)
    end

    return results
end

function capture_output(expr, SandBox::Module, term, disp, plotlib,
                        lastline)
    oldSTDOUT = STDOUT
    out = nothing
    obj = nothing
    rw, wr = redirect_stdout()
    try
        obj = eval(SandBox, expr)
        if term || disp
            obj != nothing && display(obj)
        elseif typeof(expr) == Symbol
            display(obj)
        elseif plotlib == "Gadfly" && typeof(obj) == Gadfly.Plot
            obj != nothing && display(obj)
        #This shows images and lone variables, result can
        #Handle last line sepately
        elseif lastline && obj != nothing
          #elseif mimewritable("image/png", obj) && expr.head == :call
          if expr.head == :call
            display(obj)
          end
        end
    finally
        redirect_stdout(oldSTDOUT)
        close(wr)
        out = @compat readstring(rw)
        close(rw)
    end
    return (obj, out)
end

#Parse chunk input to array of expressions
function parse_input(input::AbstractString)
    parsed = Tuple{AbstractString, Any}[]
    n = length(input)
    pos = 2 #The first character is extra line end
    while pos ≤ n
        oldpos = pos
        code,  pos = parse(input, pos)
        push!(parsed, (input[oldpos:pos-1] , code ))
    end
    parsed
end

function eval_chunk(chunk::CodeChunk, report::Report, SandBox::Module)
    info("Weaving chunk $(chunk.number) from line $(chunk.start_line)")

    if !chunk.options[:eval]
        chunk.output = ""
        chunk.options[:fig] = false
        return chunk
    end

    #Run preexecute_hooks
    for hook in preexecute_hooks
      chunk = hook(chunk)
    end

    report.fignum = 1
    report.cur_chunk = chunk

    if haskey(report.formatdict, :out_width) && chunk.options[:out_width] == nothing
        chunk.options[:out_width] = report.formatdict[:out_width]
    end

    chunk.result = run_code(chunk, report, SandBox)

    #Run post_execute chunks
    for hook in postexecute_hooks
      chunk = hook(chunk)
    end

    if chunk.options[:term]
        chunks = collect_results(chunk, TermResult())
    elseif chunk.options[:hold]
        chunks = collect_results(chunk, CollectResult())
    else
        chunks = collect_results(chunk, ScriptResult())
    end

      #else
     #   chunk.options[:fig] && (chunk.figures = copy(report.figures))
    #end

    chunks
end

#function eval_chunk(chunk::DocChunk, report::Report, SandBox)
#    chunk
#end

#Set all variables to nothing
function clear_sandbox(SandBox::Module)
    for name = names(SandBox, true)
        if name != :eval && name != names(SandBox)[1]
            try eval(SandBox, parse(AbstractString(AbstractString(name), "=nothing"))) end
        end
    end
end

function get_figname(report::Report, chunk; fignum = nothing, ext = nothing)
    figpath = joinpath(report.cwd, chunk.options[:fig_path])
    isdir(figpath) || mkpath(figpath)
    ext == nothing && (ext = chunk.options[:fig_ext])
    fignum == nothing && (fignum = report.fignum)

    chunkid = (chunk.options[:name] == nothing) ? chunk.number : chunk.options[:name]
    full_name = joinpath(report.cwd, chunk.options[:fig_path],
    "$(report.basename)_$(chunkid)_$(fignum)$ext")
    rel_name = "$(chunk.options[:fig_path])/$(report.basename)_$(chunkid)_$(fignum)$ext" #Relative path is used in output
    return full_name, rel_name
end

function init_plotting(plotlib)
    rcParams[:plotlib_set] = true
    if plotlib == nothing
        rcParams[:plotlib] = nothing
    else
        l_plotlib = lowercase(plotlib)
        rcParams[:chunk_defaults][:fig] = true
        if l_plotlib == "winston"
            eval(parse("""include(Pkg.dir("Weave","src","winston.jl"))"""))
            rcParams[:plotlib] = "Winston"
        elseif l_plotlib == "pyplot"
            eval(parse("""include(Pkg.dir("Weave","src","pyplot.jl"))"""))
            rcParams[:plotlib] = "PyPlot"
        elseif l_plotlib == "plots"
            eval(parse("""include(Pkg.dir("Weave","src","plots.jl"))"""))
            rcParams[:plotlib] = "Plots"
        elseif l_plotlib == "gadfly"
            eval(parse("""include(Pkg.dir("Weave","src","gadfly.jl"))"""))
            rcParams[:plotlib] = "Gadfly"
      end
    end
    return true
end

function get_cwd(doc::WeaveDoc, out_path)
    #Set the output directory
    if out_path == :doc
        cwd = doc.path
    elseif out_path == :pwd
        cwd = pwd()
    else
        #If there is no extension, use as path
        splitted = splitext(out_path)
        if splitted[2] == ""
            cwd = expanduser(out_path)
        else
            cwd = splitdir(expanduser(out_path))[1]
        end
    end
    return cwd
end

"""Get output file name based on out_path"""
function get_outname(out_path::Symbol, doc::WeaveDoc; ext = nothing)
    ext == nothing && (ext = doc.format.formatdict[:extension])
    outname = "$(doc.cwd)/$(doc.basename).$ext"
end

"""Get output file name based on out_path"""
function get_outname(out_path::AbstractString, doc::WeaveDoc; ext = nothing)
    ext == nothing && (ext = doc.format.formatdict[:extension])
    splitted = splitext(out_path)
    if (splitted[2]) == ""
        outname = "$(doc.cwd)/$(doc.basename).$ext"
    else
        outname = expanduser(out_path)
    end
end

function set_rc_params(formatdict, fig_path, fig_ext)
    if fig_ext == nothing
        rcParams[:chunk_defaults][:fig_ext] = formatdict[:fig_ext]
        docParams[:fig_ext] = formatdict[:fig_ext]
    else
        rcParams[:chunk_defaults][:fig_ext] = fig_ext
        docParams[:fig_ext] = fig_ext
    end
        rcParams[:chunk_defaults][:fig_path] = fig_path
        docParams[:fig_path] = fig_path
    return nothing
end

function collect_results(chunk::CodeChunk, fmt::ScriptResult)
    content = ""
    result_no = 1
    result_chunks = CodeChunk[ ]
    for r = chunk.result
        #Check if there is any output from chunk
        if strip(r.stdout) == "" && isempty(r.figures)  && strip(r.rich_output) == ""
            content *= r.code
        else
            content = "\n" * content * r.code
            rchunk = CodeChunk(content, chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
            content = ""
            rchunk.result_no = result_no
            result_no *=1
            rchunk.figures = r.figures
            rchunk.output = r.stdout * r.displayed
            rchunk.rich_output = r.rich_output
            push!(result_chunks, rchunk)
        end
    end
    if content != ""
        startswith(content, "\n") || (content = "\n" * content)
        rchunk = CodeChunk(content, chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
        push!(result_chunks, rchunk)
    end

    return result_chunks
end

function collect_results(chunk::CodeChunk, fmt::TermResult)
    output = ""
    prompt = chunk.options[:prompt]
    result_no = 1
    result_chunks = CodeChunk[ ]
    for r = chunk.result
        output *= prompt * r.code
        output *= r.displayed * r.stdout
        if !isempty(r.figures)
            rchunk = CodeChunk("", chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
            rchunk.output = output
            output = ""
            rchunk.figures = r.figures
            push!(result_chunks, rchunk)
        end
    end
    if output != ""
         rchunk = CodeChunk("", chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
         rchunk.output = output
        push!(result_chunks, rchunk)
    end

    return result_chunks
end

function collect_results(chunk::CodeChunk, fmt::CollectResult)
    result_no = 1
    for r =chunk.result
        chunk.output *=  r.stdout
        chunk.rich_output *= r.rich_output
        chunk.figures = [chunk.figures; r.figures]
    end
    return [chunk]
end

function detect_plotlib(chunk::CodeChunk)
  if isdefined(:Plots)
    init_plotting("Plots")
    #Need to set size before plots are created
    plots_set_size(chunk)
    return
  end
  isdefined(:PyPlot) && init_plotting("PyPlot") && return
  isdefined(:Gadfly) && init_plotting("Gadfly") && return
  isdefined(:Winston) && init_plotting("Winston") && return
end


text:
const PNG_BYTES_TO_CHECK = 8

const PNG_COLOR_MASK_PALETTE = 1
const PNG_COLOR_MASK_COLOR = 2
const PNG_COLOR_MASK_ALPHA = 4

const PNG_COLOR_TYPE_GRAY = 0
const PNG_COLOR_TYPE_RGB = PNG_COLOR_MASK_COLOR        # 2
const PNG_COLOR_TYPE_PALETTE = PNG_COLOR_MASK_COLOR | PNG_COLOR_MASK_PALETTE # 3
const PNG_COLOR_TYPE_GRAY_ALPHA = PNG_COLOR_MASK_ALPHA # 4
const PNG_COLOR_TYPE_RGB_ALPHA = PNG_COLOR_MASK_COLOR | PNG_COLOR_MASK_ALPHA # 6
const PNG_COLOR_TYPE_RGBA = PNG_COLOR_TYPE_RGB_ALPHA   # 6

function png_error_handler(::Ptr{Cvoid}, msg::Cstring)
    error("Png error: $msg")
end

function png_warn_handler(::Ptr{Cvoid}, msg::Cstring)
    warn("Png warn: $msg")
end
const png_error_fn = @cfunction(png_error_handler, Cvoid, (Ptr{Cvoid}, Cstring))
const png_warn_fn = @cfunction(png_warn_handler, Cvoid, (Ptr{Cvoid}, Cstring))

# Returns the libpng version string
function get_libpng_version()
    ver = ccall((:png_access_version_number, libpng), Cuint, ())
    # Version in the format of xxyyzz, where x=major, yy=minor, z=release
    # But on the major version the first x is excluded if 0.
    dig = digits(ver)[end:-1:1]
    length(dig) == 5 ? prepend!(dig, 0) : error("Unknown libpng version: $ver")

    @inbounds major = dig[1] == 0 ? string(dig[2]) : string(dig[1], dig[2])
    @inbounds minor = dig[3] == 0 ? string(dig[4]) : string(dig[3], dig[4])
    @inbounds release = dig[5] == 0 ? string(dig[6]) : string(dig[5], dig[6])

    ver_string = "$major.$minor.$release"
end

const PNG_LIBPNG_VER_STRING = get_libpng_version()

function open_png(filename::String)
    fp = ccall(:fopen, Ptr{Cvoid}, (Cstring, Cstring), filename, "rb")
    fp == C_NULL && error("Failed to open $filename")

    header = zeros(UInt8, PNG_BYTES_TO_CHECK)
    header_size = ccall(:fread, Csize_t, (Ptr{UInt8}, Cint, Cint, Ptr{Cvoid}), header, 1, PNG_BYTES_TO_CHECK, fp)
    header_size != 8 && error("Failed to read header from $filename")

    is_png = ccall((:png_sig_cmp, libpng), Cint, (Ptr{UInt8}, Csize_t, Csize_t), header, 0, PNG_BYTES_TO_CHECK)
    is_png != 0 && error("File $filename is not a png file")

    return fp
end

function create_read_struct()
    png_ptr = ccall((:png_create_read_struct, libpng), Ptr{Cvoid}, (Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
                    PNG_LIBPNG_VER_STRING, C_NULL, png_error_fn, png_warn_fn)
    png_ptr == C_NULL && error("Failed to create png read struct")
    return png_ptr
end

function create_info_struct(png_ptr)
    info_ptr = ccall((:png_create_info_struct, libpng), Ptr{Cvoid}, (Ptr{Cvoid},), png_ptr)
    info_ptr == C_NULL && error("Failed to create png info struct")
    return info_ptr
end

function png_init_io(png_ptr::Ptr{Cvoid}, fp::Ptr{Cvoid})
    ccall((:png_init_io, libpng), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), png_ptr, fp)
end

function png_set_sig_bytes(png_ptr::Ptr{Cvoid})
    ccall((:png_set_sig_bytes, libpng), Cvoid, (Ptr{Cvoid}, Cint), png_ptr, PNG_BYTES_TO_CHECK)
end

function png_read_png(png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid}, transforms::Int)
    ccall((:png_read_png, libpng), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Ptr{Cvoid}), png_ptr, info_ptr, transforms, C_NULL)
end

function png_destroy_read_struct(png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid})
    png_ptr_ptr = Ref{Ptr{Cvoid}}(png_ptr)
    info_ptr_ptr = Ref{Ptr{Cvoid}}(info_ptr)
    ccall((:png_destroy_read_struct, libpng), Cvoid, (Ref{Ptr{Cvoid}}, Ref{Ptr{Cvoid}}, Ptr{Ptr{Cvoid}}), png_ptr_ptr, info_ptr_ptr, C_NULL)
end

function close_png(fp::Ptr{Cvoid})
    ccall(:fclose, Cint, (Ptr{Cvoid},), fp)
end

function png_get_image_width(png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid})
    ccall((:png_get_image_width, libpng), UInt32, (Ptr{Cvoid}, Ptr{Cvoid}), png_ptr, info_ptr)
end

function png_get_image_height(png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid})
    ccall((:png_get_image_height, libpng), UInt32, (Ptr{Cvoid}, Ptr{Cvoid}), png_ptr, info_ptr)
end

function png_get_color_type(png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid})
    ccall((:png_get_color_type, libpng), UInt32, (Ptr{Cvoid}, Ptr{Cvoid}), png_ptr, info_ptr)
end

function png_get_bit_depth(png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid})
    ccall((:png_get_bit_depth, libpng), UInt32, (Ptr{Cvoid}, Ptr{Cvoid}), png_ptr, info_ptr)
end

function png_get_channels(png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid})
    ccall((:png_get_channels, libpng), UInt8, (Ptr{Cvoid}, Ptr{Cvoid}), png_ptr, info_ptr)
end

# Write functions

function png_create_write_struct(png_error_fn::Ptr{Cvoid}, png_warn_fn::Ptr{Cvoid})
    png_ptr = ccall((:png_create_write_struct, libpng), Ptr{Cvoid}, (Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
                    PNG_LIBPNG_VER_STRING, C_NULL, png_error_fn, png_warn_fn)
    png_ptr == C_NULL && error("Failed to create png write struct")
    return png_ptr
end

function png_create_info_struct(png_ptr::Ptr{Cvoid})
    info_ptr = ccall((:png_create_info_struct, libpng), Ptr{Cvoid}, (Ptr{Cvoid},), png_ptr)
    info_ptr == C_NULL && error("Failed to create png info struct")
    return info_ptr
end

function png_write_info(png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid})
    ccall((:png_write_info, libpng), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), png_ptr, info_ptr)
end

function png_destroy_write_struct(png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid})
    png_ptr_ptr = Ref{Ptr{Cvoid}}(png_ptr)
    info_ptr_ptr = Ref{Ptr{Cvoid}}(info_ptr)
    ccall((:png_destroy_write_struct, libpng), Cvoid, (Ref{Ptr{Cvoid}}, Ref{Ptr{Cvoid}}), png_ptr_ptr, info_ptr_ptr)
end

max_stars_repo_name: FugroRoames/LibPNG

max_issues_repo_name: FugroRoames/LibPNG

max_forks_repo_name: FugroRoames/LibPNG


text:
module MPIReco
  using RegularizedLeastSquares
  using Reexport
  @reexport using RegularizedLeastSquares
  @reexport using ImageUtils
  @reexport using MPIFiles
  const shape = MPIFiles.shape
  @reexport using DSP
  using ProgressMeter

  using LinearAlgebra
  using SparseArrays
  using Distributed
  using DistributedArrays
  # using TensorDecompositions
  using IniFile
  import LinearAlgebra: ldiv!, \

  # ML Stuff
  import Distributions
  using Flux
  using CUDA
  using BSON

  include("Utils.jl")
  include("MultiContrast.jl")
  include("RecoParameters.jl")
  include("SystemMatrix/SystemMatrix.jl")
  include("Weighting.jl")
  include("Background.jl")
  include("TemporalRegularization/TemporalRegularization.jl")
  include("Reconstruction.jl")
  include("MultiPatch.jl")
  include("MotionCompensation/MotionCompensation.jl")
  include("MachineLearning/MachineLearning.jl")
end # module


text:
# Test the solid state table. Mostly data entry / lookups work.

println(SST[AtoZ["B"]])

@test SST[AtoZ["Be"]].kF == 1.94
@test SST[AtoZ["As"]].ϵp == 7.91
@test SST[AtoZ["I"]].ϵs == 19.42


text:
using CUDA
import MPI
using ImplicitGlobalGrid
using ParallelStencil
using ParallelStencil.FiniteDifferences3D

include("part1_utils.jl")

@parallel function diffusion_3D_step_τ!(Ht, Hτ, dHdτ, dt, dτ, qx, qy, qz, dx, dy, dz, D)
    @all(qx) = D * @d_xi(Hτ) / dx
    @all(qy) = D * @d_yi(Hτ) / dy
    @all(qz) = D * @d_zi(Hτ) / dz

    @all(dHdτ) =
        -(@inn(Hτ) - @inn(Ht)) / dt + (@d_xa(qx) / dx + @d_ya(qy) / dy + @d_za(qz) / dz)
    @inn(Hτ) = @inn(Hτ) + @all(dHdτ) * dτ
    return nothing
end

@views function diffusion_3D_array_programming(; nx, ny, nz, do_vis = false, verbose=true, init_and_finalize_MPI=!isinteractive())
    # physics
    lx, ly, lz = 10.0, 10.0, 10.0
    D = 1.0
    ttot = 1.0

    # derived numerics
    me, dims, _nprocs, coords, comm_cart =
        init_global_grid(nx, ny, nz; init_MPI = init_and_finalize_MPI, quiet=!verbose)
    dx, dy, dz = lx / nx_g(), ly / ny_g(), lz / nz_g()

    # bind MPI ranks to GPUs
    if USE_GPU
        select_device()
    end

    total_N = prod(dims) * nx * ny * nz

    # numerics
    dt = 0.2
    dτ = min(dx, dy, dz)^2 ./ D / 8.1
    nt = cld(ttot, dt)
    tol = 1e-8
    iter_max = 1e5

    # array allocation
    center = [lx / 2, ly / 2, lz / 2]
    qx = @zeros(nx - 1, ny - 2, nz - 2)
    qy = @zeros(nx - 2, ny - 1, nz - 2)
    qz = @zeros(nx - 2, ny - 2, nz - 1)
    Ht = @zeros(nx, ny, nz)
    Ht .= init_local_gaussian(center, dx, dy, dz, Ht)
    apply_boundary_conditions!(Ht, coords, dims)
    Hτ = copy(Ht)

    dHdt = @zeros(nx - 2, ny - 2, nz - 2)
    H_g = zeros(nx * dims[1], ny * dims[2], nz * dims[3])

    t = 0.0
    iter_outer = 0
    iter_total = 0
    while t < ttot

        iter_inner = 0
        err = 2 * tol
        while err > tol && iter_inner < iter_max
            @parallel diffusion_3D_step_τ!(Ht, Hτ, dHdt, dt, dτ, qx, qy, qz, dx, dy, dz, D)
            update_halo!(Hτ)
            err = dist_norm_L2(dHdt * dt, comm_cart) / sqrt(total_N)
            iter_inner += 1
        end
        if verbose && me == 0
            if err <= tol
                println("Converged after $iter_inner iterations.")
            else
                println("Couldn't converge within $iter_max iterations.")
            end
        end
        iter_total += iter_inner
        iter_outer += 1
        t += dt
        Ht .= Hτ

    end

    X_g = LinRange(0 + dx/2, lx - dx/2, nx * dims[1])

    gather!(Array(Ht), H_g)

    finalize_global_grid(; finalize_MPI = init_and_finalize_MPI)

    return X_g, H_g
end