text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

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

#"pseudo-types" of aspheres
"""
AphericSurfaces polynomial evaluation is optimized for no terms `CONIC`, odd terms `ODD`, even terms `EVEN` or any `ODDEVEN`
"""
@enum AsphericSurfaceType CONIC ODD EVEN ODDEVEN

"""
    AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N}

Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter,.
`T` is the datatype, `N` is the dimensionality, `Q` is the number of aspheric terms, and `M` is the type of aspheric polynomial.

```julia
AsphericSurface(semidiameter; radius, conic, aspherics=nothing, normradius = semidiameter)
```

The surface is centered at the origin and treated as being the cap of an infinite cylinder, thus creating a true half-space.
Outside of `0 <= ρ <= 1` the height of the surface is not necessarily well defined, so NaN may be returned.

`aspherics`  aspherics should be vectors containing tuples of the form (i, v) where i is the polynomial power of the aspheric term.
An empty vector is not permitted. Use `nothing` instead.

The sag is defined by the equation

```math
z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir^{i}
```

where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, and ``k = \\texttt{conic}`` .

The function checks if the aspheric terms are missing, even, odd or both and uses an efficient polynomial evaluation strategy.

"""
struct AsphericSurface{T, N, Q, M} <: ParametricSurface{T, N}
    semidiameter::T
    curvature::T
    conic::T
    aspherics::SVector{Q, T}
    normradius::T
    boundingcylinder::Cylinder{T, N}

    function AsphericSurface(
        M::AsphericSurfaceType,
        semidiameter::T,
        curvature::T,
        conic::T,
        aspherics::SVector{Q, T},
        normradius::T,
        boundingcylinder,
    ) where {T <: Real, Q}
        return new{T, 3, Q, M}(
            semidiameter,
            curvature,
            conic,
            aspherics,
            normradius,
            boundingcylinder,
        )
    end

    function AsphericSurface(
        semidiameter::T;
        radius::T = typemax(T),
        conic::T = zero(T),
        aspherics::Union{Nothing, Vector{Tuple{Int, T}}} = nothing,
        normradius::T = semidiameter,
    ) where {T <: Real}
        @assert semidiameter > 0
        @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic)
        @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)"
        @assert aspherics != []

        acs = []
        if aspherics === nothing
            M = CONIC
        else
            asphericTerms = [i for (i,) in aspherics]
            minAsphericTerm = minimum(asphericTerms)
            maxAsphericTerm = maximum(asphericTerms)
            @assert minAsphericTerm > 0 "Aspheric Terms must be Order 1 or higher ($minAsphericTerm)"
            acs = zeros(T, maxAsphericTerm)
            for (i, k) in aspherics
                acs[i] = k
            end
            odd = any([isodd(i) && a != zero(T) for (i, a) in enumerate(acs)])
            even = any([iseven(i) && a != zero(T) for (i, a) in enumerate(acs)])
            if odd && even
                M = ODDEVEN
            elseif even
                M = EVEN
                acs = [acs[2i] for i in 1:(maxAsphericTerm ÷ 2)]
            elseif odd
                M = ODD
                acs = [acs[2i - 1] for i in 1:((maxAsphericTerm + 1) ÷ 2)]
            else #there are no nonzero aspherics terms in the list
                M = CONIC
                acs = []
            end
        end
        Q = length(acs)
        surf = new{T, 3, Q, M}(
            semidiameter,
            1 / radius,
            conic,
            acs,
            normradius,
            Cylinder(semidiameter, interface = opaqueinterface(T)),
        )
        return surf
    end
end

"""
    EvenAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter)

Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter.

`aspherics` should be an array of the even coefficients of the aspheric polynomial starting with A2

"""
function EvenAsphericSurface(
    semidiameter::T,
    curvature::T,
    conic::T,
    aspherics::Vector{T};
    normradius::T = semidiameter,
) where {T <: Real}
    Q = length(aspherics)
    return AsphericSurface(
        EVEN,
        semidiameter,
        curvature,
        conic,
        SVector{Q, T}(aspherics),
        normradius,
        Cylinder(semidiameter, interface = opaqueinterface(T)),
    )
end

"""
    OddAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter)

Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter.

`aspherics`  should be an array of the odd coefficients of the aspheric polynomial starting with A1

"""
function OddAsphericSurface(
    semidiameter::T,
    curvature::T,
    conic::T,
    aspherics::Vector{T};
    normradius::T = semidiameter,
) where {T <: Real}
    Q = length(aspherics)
    return AsphericSurface(
        ODD,
        semidiameter,
        curvature,
        conic,
        SVector{Q, T}(aspherics),
        normradius,
        Cylinder(semidiameter, interface = opaqueinterface(T)),
    )
end

"""
    OddEvenAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter)

Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter.

`aspherics` should be an array of the both odd and even coefficients of the aspheric polynomial starting with A1

"""
function OddEvenAsphericSurface(
    semidiameter::T,
    curvature::T,
    conic::T,
    aspherics::Vector{T};
    normradius::T = semidiameter,
) where {T <: Real}
    Q = length(aspherics)
    return AsphericSurface(
        ODDEVEN,
        semidiameter,
        curvature,
        conic,
        SVector{Q, T}(aspherics),
        normradius,
        Cylinder(semidiameter, interface = opaqueinterface(T)),
    )
end

#don't have a function for CONIC as it would better to directly solve for the surface intersection instead of the rootfinding of a ParametricSurface

"""
    asphericType(surf::AsphericSurface)

Query the polynomial type of `asp.  Returns CONIC, ODD, EVEN, or ODDEVEN. CONIC corresponds to no aspheric terms, ODD
means it only has odd aspheric terms, EVEN means only even aspheric terms and ODDEVEN means both even and odd terms.

This function is to enable proper interpretation of `surf.aspherics` by any optimization routines that directly query the aspheric coefficients.
"""
asphericType(z::AsphericSurface{T, 3, Q, M}) where {T <: Real, Q, M} = M

function uvrange(::AsphericSurface{T, N}) where {T <: Real, N}
    return ((zero(T), one(T)), (-T(π), T(π)))
end # ρ and ϕ

semidiameter(z::AsphericSurface{T}) where {T <: Real} = z.semidiameter
halfsizeu(z::AsphericSurface{T}) where {T <: Real} = semidiameter(z)
halfsizev(z::AsphericSurface{T}) where {T <: Real} = semidiameter(z)

boundingobj(z::AsphericSurface{T}) where {T <: Real} = z.boundingcylinder

function prod_step(
    z::AsphericSurface{T, N, Q, ODDEVEN},
    r,
    r2,
) where {T <: Real, N, Q}
    return r, r
end
function prod_step(
    z::AsphericSurface{T, N, Q, ODD},
    r,
    r2,
) where {T <: Real, N, Q}
    return r, r2
end
function prod_step(
    z::AsphericSurface{T, N, Q, EVEN},
    r,
    r2,
) where {T <: Real, N, Q}
    return r2, r2
end

function point(
    z::AsphericSurface{T, 3, Q, M},
    ρ::T,
    ϕ::T,
)::SVector{3, T} where {T <: Real, Q, M}
    rad = z.semidiameter
    r = ρ * rad
    r2 = r^2
    t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2
    if t < zero(T)
        return SVector{3, T}(NaN, NaN, NaN)
    end
    h = z.curvature * r2 / (one(T) + sqrt(t))
    # sum aspheric
    if M != CONIC
        prod, step = prod_step(z, r, r2)  #multiple dispatch on M
        asp, rest = Iterators.peel(z.aspherics)
        h += asp * prod
        for asp in rest
            prod *= step
            h += asp * prod
        end
    end
    return SVector{3, T}(r * cos(ϕ), r * sin(ϕ), h)
end

function partial_prod_step(
    z::AsphericSurface{T, 3, Q, EVEN},
    r::T,
    r2::T,
) where {T <: Real, Q}
    return r, r2, 2:2:(2Q)
end
function partial_prod_step(
    z::AsphericSurface{T, 3, Q, ODD},
    r::T,
    r2::T,
) where {T <: Real, Q}
    return one(T), r2, 1:2:(2Q - 1)
end
function partial_prod_step(
    z::AsphericSurface{T, 3, Q, ODDEVEN},
    r::T,
    r2::T,
) where {T <: Real, Q}
    return one(T), r, 1:1:Q
end

function partials(
    z::AsphericSurface{T, 3, Q, M},
    ρ::T,
    ϕ::T,
)::Tuple{SVector{3, T}, SVector{3, T}} where {T <: Real, Q, M}
    rad = z.semidiameter
    r = ρ * rad
    r2 = r * r
    t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2
    if t < zero(T)
        return SVector{3, T}(NaN, NaN, NaN), SVector{3, T}(NaN, NaN, NaN)
    end
    dhdρ = rad * z.curvature * r * sqrt(t) / t
    # sum aspherics partial
    if M != CONIC
        prod, step, mIter = partial_prod_step(z, r, r2)

        ((m, asp), rest) = Iterators.peel(zip(mIter, z.aspherics))
        dhdρ += rad * asp * 2 * prod #first term m=1*2 and prod = r
        for (m, asp) in rest
            prod *= step
            dhdρ += rad * m * asp * prod
        end
    end
    dhdϕ = zero(T)
    cosϕ = cos(ϕ)
    sinϕ = sin(ϕ)
    return SVector{3, T}(rad * cosϕ, rad * sinϕ, dhdρ),
    SVector{3, T}(r * -sinϕ, r * cosϕ, dhdϕ)
end

function normal(
    z::AsphericSurface{T, 3},
    ρ::T,
    ϕ::T,
)::SVector{3, T} where {T <: Real}
    du, dv = partials(z, ρ, ϕ)
    if ρ == zero(T) && norm(dv) == zero(T)
        # in cases where there is no δϕ at ρ = 0 (i.e. anything which is rotationally symetric)
        # then we get some big problems, hardcoding this case solves the problems
        return SVector{3, T}(0, 0, 1)
    end
    return normalize(cross(du, dv))
end

function uv(z::AsphericSurface{T, 3}, p::SVector{3, T}) where {T <: Real}
    # avoid divide by zero for ForwardDiff
    ϕ = NaNsafeatan(p[2], p[1])
    if p[1] == zero(T) && p[2] == zero(T)
        ρ = zero(T)
    else
        ρ = sqrt(p[1]^2 + p[2]^2) / semidiameter(z)
    end
    return SVector{2, T}(ρ, ϕ)
end

function onsurface(
    surf::AsphericSurface{T, 3},
    p::SVector{3, T},
) where {T <: Real}
    ρ, ϕ = uv(surf, p)
    if ρ > one(T)
        return false
    else
        surfpoint = point(surf, ρ, ϕ)
        return samepoint(p[3], surfpoint[3])
    end
end

function inside(surf::AsphericSurface{T, 3}, p::SVector{3, T}) where {T <: Real}
    ρ, ϕ = uv(surf, p)
    if ρ > one(T)
        return false
    else
        surfpoint = point(surf, ρ, ϕ)
        return p[3] < surfpoint[3]
    end
end

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

# Assumes the ray has been transformed into the canonical coordinate frame which has the vertical axis passing through (0,0,0) and aligned with the z axis.
function surfaceintersection(
    surf::AcceleratedParametricSurface{T, 3, AsphericSurface{T, 3}},
    r::AbstractRay{T, 3},
) where {T <: Real}
    cylint = surfaceintersection(surf.surface.boundingcylinder, r)
    if cylint isa EmptyInterval{T}
        return EmptyInterval(T)
    else
        if doesintersect(surf.triangles_bbox, r) ||
           inside(surf.triangles_bbox, origin(r))
            surfint = triangulatedintersection(surf, r)
            if !(surfint isa EmptyInterval{T})
                return intervalintersection(cylint, surfint)
            end
        end
        # hasn't hit the surface
        if lower(cylint) isa RayOrigin{T} && upper(cylint) isa Infinity{T}
            if inside(surf.surface, origin(r))
                return Interval(RayOrigin(T), Infinity(T))
            else
                return EmptyInterval(T)
            end
            # otherwise check that the intersection is underneath the surface
        else
            p = point(closestintersection(cylint, false))
            ρ, ϕ = uv(surf, p)
            surfpoint = point(surf.surface, ρ, ϕ)
            if p[3] < surfpoint[3]
                return cylint # TODO!! UV (and interface) issues?
            else
                return EmptyInterval(T)
            end
        end
    end
end

function AcceleratedParametricSurface(
    surf::T,
    numsamples::Int = 17;
    interface::NullOrFresnel{S} = NullInterface(S),
) where {S <: Real, N, T <: AsphericSurface{S, N}}
    # Zernike uses ρ, ϕ uv space so need to modify extension of triangulation
    a = AcceleratedParametricSurface(
        surf,
        triangulate(surf, numsamples, true, false, true, false),
        interface = interface,
    )
    emptytrianglepool!(S)
    return a
end

function asphericSag(surf::AsphericSurface{T, 3, Q, EVEN}) where {T <: Real, Q}
    amin = sum(
        k < zero(T) ? k * surf.semidiameter^(2m) : zero(T) for
        (m, k) in enumerate(surf.aspherics)
    )
    amax = sum(
        k > zero(T) ? k * surf.semidiameter^(2m) : zero(T) for
        (m, k) in enumerate(surf.aspherics)
    )
    return amin, amax
end

function asphericSag(surf::AsphericSurface{T, 3, Q, ODD}) where {T <: Real, Q}
    amin = sum(
        k < zero(T) ? k * surf.semidiameter^(2m - 1) : zero(T) for
        (m, k) in enumerate(surf.aspherics)
    )
    amax = sum(
        k > zero(T) ? k * surf.semidiameter^(2m - 1) : zero(T) for
        (m, k) in enumerate(surf.aspherics)
    )
    return amin, amax
end

function asphericSag(
    surf::AsphericSurface{T, 3, Q, ODDEVEN},
) where {T <: Real, Q}
    amin = sum(
        k < zero(T) ? k * surf.semidiameter^(m) : zero(T) for
        (m, k) in enumerate(surf.aspherics)
    )
    amax = sum(
        k > zero(T) ? k * surf.semidiameter^(m) : zero(T) for
        (m, k) in enumerate(surf.aspherics)
    )
    return amin, amax
end

function asphericSag(surf::AsphericSurface{T, 3, Q, CONIC}) where {T <: Real, Q}
    amin = zero(T)
    amax = zero(T)
    return amin, amax
end

function BoundingBox(surf::AsphericSurface{T, 3}) where {T <: Real}
    xmin = -semidiameter(surf)
    xmax = semidiameter(surf)
    ymin = -semidiameter(surf)
    ymax = semidiameter(surf)
    # aspherics can be more comlpicated, so just take sum of all negative and all positives
    amin, amax = asphericSag(surf)
    zmin = amin
    zmax = amax
    # curvature only goes one way
    q = one(T) - (one(T) + surf.conic) * surf.curvature^2 * surf.semidiameter^2
    if q < zero(T)
        throw(
            ErrorException(
                "The surface is invalid, no bounding box can be constructed",
            ),
        )
    end
    hmax = surf.curvature * surf.semidiameter^2 / (one(T) + sqrt(q))
    if hmax > zero(T)
        zmax += hmax
    else
        zmin += hmax
    end
    return BoundingBox(xmin, xmax, ymin, ymax, zmin, zmax)
end


text: # import JuliaWebAPI

ext: jl


text:
struct Camera{T <: AbstractFloat}
	origin::Vec3{T}
	lower_left_corner::Vec3{T}
	horizontal::Vec3{T}
	vertical::Vec3{T}
	u::Vec3{T}
	v::Vec3{T}
	w::Vec3{T}
	lens_radius::T
end

"""
	Args:
		vfov: vertical field-of-view in degrees
		aspect_ratio: horizontal/vertical ratio of pixels
      aperture: if 0 - no depth-of-field
"""
function default_camera(lookfrom::Vec3{T}=(SA{T}[0,0,0]),
						lookat::Vec3{T}=(SA{T}[0,0,-1]),
						vup::Vec3{T}=(SA{T}[0,1,0]),
						vfov::T=T(90), aspect_ratio::T=T(16/9),
						aperture::T=T(0), focus_dist::T=T(1)) where T
	viewport_height = T(2) * tand(vfov/T(2))
	viewport_width = aspect_ratio * viewport_height

	w = normalize(lookfrom - lookat)
	u = normalize(vup × w)
	v = w × u

	origin = lookfrom
	horizontal = focus_dist * viewport_width * u
	vertical = focus_dist * viewport_height * v
	lower_left_corner = origin - horizontal/T(2) - vertical/T(2) - focus_dist*w
	lens_radius = aperture/T(2)
	Camera{T}(origin, lower_left_corner, horizontal, vertical, u, v, w, lens_radius)
end

default_camera(lookfrom, lookat, vup, vfov, aspect_ratio, aperture, focus_dist; elem_type::Type{T}) where T =
	default_camera(Vec3{T}(lookfrom), Vec3{T}(lookat), Vec3{T}(vup),
		T(vfov), T(aspect_ratio), T(aperture), T(focus_dist)
	)

@inline @fastmath function get_ray(c::Camera{T}, s::T, t::T) where T
	rd = SVector{2,T}(c.lens_radius * random_vec2_in_disk(T))
	offset = c.u * rd.x + c.v * rd.y #offset = c.u * rd.x + c.v * rd.y
    Ray(c.origin + offset, normalize(c.lower_left_corner + s*c.horizontal +
									 t*c.vertical - c.origin - offset))
end


text:
#File: maxSharpe.jl
#Author:
#Date: 29-June-2020

#OBJECTIVE FUNCTION: MAXIMISE Sharpe Ratio
function maxSharpe(W,R,rf,VCV)
###########################################################################
#INPUTS:
#R     - returns vector
#rf    - risk-free rate
#VCV   - VCV matrix
###########################################################################
    R = R-rf;
    fval = -W'*R/sqrt(W'*VCV*W);
    return fval
end


text:
using Merly
using Test
using HTTP
using JSON

ip = "127.0.0.1"
port = 8086

@test webserverpath("website") == joinpath(pwd(),"website")
@test File("index.html") == read( open(joinpath(pwd(),"index.html")), String)

useCORS(AllowOrigins = "*", AllowHeaders = "*" , AllowMethods = "GET,POST,PUT,DELETE", MaxAge = "80")

headersalways(["X-PINGOTHER" => "pingpong"])

notfound("""<!DOCTYPE html>

              <title>Not found</title>

404, Not found

              """)

notfound("notfound.html")

webserverfiles("*")

u = 1

function tojson(data::String)
   return JSON.parse(data)
end

formats["application/json"] = tojson

@page "/" HTTP.Response(200,"Hello World!")

@page "/hola/:usr" HTTP.Response(200,string("Hello ",request.params["usr"],"!"))

@page "/mifile" HTTP.Response(200, File("index.html"))

@page "/get1" (;u=u) HTTP.Response(200,string("Get1 ",u," !"))

@route GET "/get/:data1" (;u=u) begin
  u = u +1
  HTTP.Response(200, string(u ,request.params["data1"]))
end

@route GET "/get/data/hola" begin
  HTTP.Response(200,"get this back: data")
end

@route GET "/regex/(\\w+\\d+)" begin
  return HTTP.Response(200, string("datos ",request.params["2"]))# $(request.params[1])"
end

@route POST "/post" begin
  HTTP.Response(200,"I did something!")
end

@route POST|PUT|DELETE "/" begin
  println("query: ",request.query)
  println("body: ",request.body)

  HTTP.Response(200
          , HTTP.mkheaders(["Content-Type" => "text/plain"])
          , body="I did something2!")
end

Get("/data", (request,HTTP)->begin

  println("params: ",request.params)
  println("query: ",request.query)

  HTTP.Response(200
          , HTTP.mkheaders(["Content-Type" => "text/plain"])
          , body=string(u,"data ", get(request.query,"hola","")))

end)

Post("/data", (request,HTTP)-> begin
  println("params: ",request.params)
  println("query: ",request.query)
  println("body: ",request.body)

  global u="bye"

  HTTP.Response(200
          , HTTP.mkheaders(["Content-Type" => "text/plain"])
          , body=string("I did something! ", request.body["query"]))

end)

 Put("/data", (request,HTTP) -> begin
     HTTP.response(200, "test text put")
 end)
 Delete("/data", (request,HTTP) -> begin
     HTTP.response(200, "test text delete")
 end)
 Connect("/data", (request,HTTP) -> begin
     HTTP.response(200, "test text Connect")
 end)
 Trace("/data", (request,HTTP) -> begin
     HTTP.response(200, "test text Trace")
 end)
 Head("/data", (request,HTTP) -> begin
     HTTP.response(200, "test text head")
 end)
 Patch("/data", (request,HTTP) -> begin
     HTTP.response(200, "test text Patch")
 end)

Get("/test1/:usr",
  (request, HTTP) -> begin
        HTTP.Response(200,string("test1 ",request.params["usr"],"!"))
    end
)

Get("/test2/:usr",
    (result(;u=u) = (request, HTTP)-> begin
          u= u+1
          HTTP.Response(200,string("test2 ",u,request.params["usr"]," !"))
        end)()
)

@route GET "/test3/" (;u=u) begin
  u= u+1
  HTTP.Response(200,string("get this back:", u))
end

function authenticate(request, HTTP)
  isAuthenticated = false

  if (request.params["status"] === "authenticated")
    isAuthenticated = true
  end

  return request, HTTP, isAuthenticated
end

Get("/verify/:status",

  (result(;middleware=authenticate) = (request, HTTP)-> begin

  myfunction = (request, HTTP, isAuthenticated)-> begin

      if (isAuthenticated == false )
          return  HTTP.Response(403,string("Unauthenticated. Please signup!"))
      end
          return  HTTP.Response(200,string("verify !"))
      end

  return myfunction(middleware(request,HTTP)...)

  end)()

)

 @async start(host = ip, port = port)
 sleep(2)

myjson = Dict("query"=>"data")

my_headers = HTTP.mkheaders(["Content-Type" => "application/json"])

  global r = HTTP.get("http://$(ip):$(port)/")
 @test r.status == 200
 @test String(r.body) == "Hello World!"

 global r = HTTP.get("http://$(ip):$(port)/hola/usuario")
@test r.status == 200
@test String(r.body) == "Hello usuario!"

 global r = HTTP.get("http://$(ip):$(port)/get/testdata")
@test r.status == 200
@test String(r.body) == "2testdata"

 global r = HTTP.get("http://$(ip):$(port)/get/data")
@test r.status == 200
@test String(r.body) == "3data"

 global r = HTTP.get("http://$(ip):$(port)/get/data/hola")
@test r.status == 200
@test String(r.body) == "get this back: data"

r= HTTP.get("http://$(ip):$(port)/regex/test1")
@test String(r.body) == "datos test1"

 global r = HTTP.get("http://$(ip):$(port)/data?hola=1")
@test r.status == 200
@test String(r.body) == "1data 1"

myjson = Dict("query"=>"data")
my_headers = HTTP.mkheaders(["Accept" => "application/json","Content-Type" => "application/json"])
 global r = HTTP.post("http://$(ip):$(port)/data",my_headers,JSON.json(myjson))
@test r.status == 200
@test String(r.body) == "I did something! data"

 global r = HTTP.get("http://$(ip):$(port)/data")
@test r.status == 200
@test String(r.body) == "byedata "

 global r = HTTP.get("http://$(ip):$(port)/test1/neomatrix")
@test r.status == 200
@test String(r.body) == "test1 neomatrix!"

 global r = HTTP.get("http://$(ip):$(port)/test2/neomatrixcode")
@test r.status == 200
@test String(r.body) == "test2 2neomatrixcode !"

 global r = HTTP.get("http://$(ip):$(port)/get1")
@test r.status == 200
@test String(r.body) == "Get1 1 !"

 global r = HTTP.get("http://$(ip):$(port)/test3")
@test r.status == 200
@test String(r.body) == "get this back:2"

my_headers = HTTP.mkheaders(["Accept" => "application/json"])
 global r = HTTP.post("http://$(ip):$(port)/post",my_headers,JSON.json(myjson))
@test r.status == 200
@test String(r.body) == "I did something!"

my_headers = HTTP.mkheaders(["Accept" => "application/json","Content-Type" => "application/json"])
 global r = HTTP.post("http://$(ip):$(port)/",my_headers,JSON.json(myjson))
@test r.status == 200
@test String(r.body) == "I did something2!"

 global r = HTTP.put("http://$(ip):$(port)/",my_headers,JSON.json(myjson))
@test r.status == 200
@test String(r.body) == "I did something2!"

 global r = HTTP.delete("http://$(ip):$(port)/")
@test r.status == 200
@test String(r.body) == "I did something2!"
@test r.headers == Pair{SubString{String},SubString{String}}["Content-Type" => "text/plain", "Access-Control-Allow-Origin" => "*", "X-PINGOTHER" => "pingpong", "Transfer-Encoding" => "chunked"]

try
   global r = HTTP.get("http://$(ip):$(port)/nada")
catch e
  @test String(e.response.body) == File("notfound.html")
end

 global r = HTTP.get("http://$(ip):$(port)/prueba.txt")
@test r.status == 200
@test String(r.body) == "probando webserver"
@test r.headers == Pair{SubString{String},SubString{String}}["Content-Type" => "text/plain", "Access-Control-Allow-Origin" => "*", "X-PINGOTHER" => "pingpong", "Transfer-Encoding" => "chunked"]

global r= HTTP.get("http://$(ip):$(port)/testfolder/testfile.txt")
@test String(r.body) == "test text"

global r= HTTP.get("http://$(ip):$(port)/mifile")
@test r.status == 200

if Sys.iswindows()
    @test String(r.body) ==  "<!DOCTYPE html>\r\n\r\n\r\n\t<meta charset=\"UTF-8\">\r\n\t<title>Document</title>\r\n\r\n\r\n
hola
\r\n\r\n"
else
   @test String(r.body) ==  "<!DOCTYPE html>\n\n\n\t<meta charset=\"UTF-8\">\n\t<title>Document</title>\n\n\n
hola
\n\n"
end

global r= HTTP.get("http://$(ip):$(port)/verify/authenticated")
@test r.status == 200
@test String(r.body) == "verify !"

try
   global r= HTTP.get("http://$(ip):$(port)/verify/unauthenticated")
catch e
  @test e.response.status == 403
  @test String(e.response.body) == "Unauthenticated. Please signup!"
end


text:
#===============================================================================

  Code for investment under uncertainty by Heterogenous firms

  Translated from matlab code from Ben Moll's website:
      www.princeton.edu/~moll/HACTproject.htm

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

using LinearAlgebra, SparseArrays, Plots

# Set up model parameters
ρ = 0.05  #discount rate
α = 0.5 # Cobb-douglas production function
Θ = 2.7 # quadratic adjustment cost
δ = 0.025 # depreciation rate

# Parameters for the Ornstein-Uhlenbeck process for z
var = 0.026^2
zmean = exp(var/2)
corr = 0.859
θ = -log(corr)
σ_sq = 2*θ*var

# Set up the gridspace for z
J = 40
zmin = zmean*0.6
zmax=zmean*1.4
z = LinRange(zmin,zmax,J)
dz = (zmax-zmin)/(J-1)
dz_sq = dz^2

# Set up gridspace for capital
H =100
kmin=1
kmax = 100
k = LinRange(kmin,kmax, H)
dk = (kmax-kmin)/(H-1)

kk = k*ones(1,J)
zz = (z*ones(1,H))'

# Set up terms for drift and variance of the UO process that depend on zz
μ = (-θ*log.(z) .+ σ_sq/2).*z
Σ_sq = σ_sq*(z.^2)

# Production function
F = zz.*kk.^α

# Parameters for our loop
maxit = 20
crit=10^(-6)
Δ = 1000

# Empty spaces for saving things
Vkb, Vkf, c = [zeros(H,J) for i in 1:3]
dist =[]
x_all=[]

# Create the matrix that summarizes the evolution of our system in the z dimension
 yy = (-Σ_sq/dz_sq + min.(μ,0)/dz - max.(μ,0)/dz)
 χ = Σ_sq/(2*dz_sq)- min.(μ,0)/dz
 ζ = max.(μ,0)/dz + Σ_sq/(2*dz_sq)

 # Define the Diagonals of this matrix
 updiag_z = zeros(H,1)
 	for j = 1:J
		global updiag_z =[updiag_z; repeat([ζ[j]], H, 1)]
	end

 centerdiag_z=repeat([χ[1]+yy[1]],H,1)
	for j = 2:J-1
		global centerdiag_z = [centerdiag_z; repeat([yy[j]], H, 1)]
	end
 centerdiag_z=[centerdiag_z; repeat([yy[J]+ζ[J]], H, 1)]

lowdiag_z = repeat([χ[2]], H, 1)
	for j=3:J
		global lowdiag_z = [lowdiag_z; repeat([χ[j]],H,1)]
	end
# spdiags in Matlab allows for automatic trimming/adding of zeros
    # spdiagm does not do this

Bswitch =  spdiagm(0=>centerdiag_z[:], -H=>lowdiag_z[:], H=>updiag_z[H+1:end-H])

# Inital guess for v
v0 = (F-δ.*kk-0.5*Θ*δ^2 .*kk)/ρ
global v =v0

for n in 1:maxit
	global V=v
    #forward difference
    Vkf[1:H-1,:] = (V[2:H,:]-V[1:H-1,:])/dk
    Vkf[H,:] = (1+Θ*δ)*ones(1,J)#will never be used, but impose state constraint a<=amax just in case
    #backward difference
    Vkb[2:H,:] = (V[2:H,:]-V[1:H-1,:])/dk
    Vkb[1,:] = (1+Θ*δ)*ones(1,J) #state constraint boundary condition

    #I_concave = Vab .> Vaf #indicator whether value function is concave (problems arise if this is not the case)

    # investment and savings with forward difference
    xf = (Vkf.-1)/Θ.*kk
    ssf = xf-δ.*kk
	Hf = F -xf -0.5*Θ*(xf./kk).^2 .*kk + Vkf.*ssf
    # consumption and savings with backward difference
	xb = (Vkb.-1)/Θ.*kk
    ssb = xb-δ.*kk
	Hb = F -xb -0.5*Θ*(xb./kk).^2 .*kk + Vkb.*ssb
    # investment at steady state
    x0 = δ.*kk

    #= dV_upwind makes a choice of forward or backward differences based on
     the sign of the drift    =#
	Ieither = (1 .- (ssf.>0)) .* (1 .- (ssb.<0))
	Iunique = (ssb.<0).*(1 .- (ssf.>0)) + (1 .- (ssb.<0)).*(ssf.>0)
	Iboth = (ssb.<0).*(ssf.>0)
    If= Iunique.*(ssf .> 0) + Iboth.*(Hf.>=Hb)  #positive drift → forward difference
    Ib = Iunique.*(ssb .< 0) + Iboth.*(Hb.>=Hf) #negative drift → backward difference
    I0 = Ieither  #at steady state

	x = xf.*If + xb.*Ib + x0.*I0
    profits = F-x-0.5*Θ*(x./kk).^2. .*kk

    # CONSTRUCT MATRIX
    X = -ssb.*Ib/dk
    Y = - ssf.*If/dk  + ssb.*Ib/dk
    Z = ssf.*If/dk

	# Define the Diagonals of this transition matrix
	updiag = 0
       for j = 1:J
           updiag =[updiag; Z[1:H-1,j]; 0]
       end
    updiag =(updiag[:])

    centerdiag=reshape(Y, H*J, 1)
    centerdiag = (centerdiag[:]) # for tuples

   lowdiag = X[2:H, 1]
       for j = 2:J
           lowdiag = [lowdiag; 0; X[2:H,j]]
       end
   lowdiag=(lowdiag)

   # spdiags in Matlab allows for automatic trimming/adding of zeros
       # spdiagm does not do this
  	AA = sparse(Diagonal(centerdiag))+ [zeros(1, H*J); sparse(Diagonal(lowdiag)) zeros(H*J-1,1)] + sparse(Diagonal(updiag))[2:end, 1:(H*J)] # trim first element
	A = AA + Bswitch

    B = (1/Δ + ρ)*sparse(I, H*J, H*J) - A

    profits_stacked = reshape(profits, H*J,1)
    V_stacked = reshape(V, H*J,1)

    b = profits_stacked + V_stacked/Δ
    V_stacked = B\b #SOLVE SYSTEM OF EQUATIONS

    global V = reshape(V_stacked, H,J)

    Vchange = V-v
	global v = V

	push!(dist, findmax(abs.(Vchange))[1])
	push!(x_all,x)

	if dist[n]<crit
		println("Value function converged iteration = $(n)")
		break
	end
end

plot(k,v, legend=false)

plot(dist, legend=false)

kdot=x_all[end]-δ.*kk
plot(kk,kdot,label="")
plot!(kk,zeros(H,1), line=:dash, color=:black, label="")


text:
###############################################################################
#
#  Isogenies
#
#
#  (C) 2022 Jeroen Hanselman
#
###############################################################################

export Isogeny

export isogeny_from_kernel, domain, codomain, degree, image

mutable struct Isogeny
  domain::EllCrv
  codomain::EllCrv
  coordx::RingElem
  coordy::RingElem
  kernel::Vector{EllCrvPt}
  kernel_polynomial::RingElem
  degree::Int

  function Isogeny(E,psi)
    r = new()
    r.domain = E
    r.kernel_polynomial = psi
    r.codomain, r.coordx, r.coordy, r.degree = isogeny_kernel_polynomial(E,psi)
    return r
  end

end

@doc Markdown.doc"""
    isogeny_from_kernel(E::EllCrv, psi::RingElem) -> Isogeny

Compute the isogeny $\phi$: $E$ -> $E'$ where the kernel of $\phi$ contains
the points whose $x$-coordinates are roots of the input polynomial $\psi$.
"""
function isogeny_from_kernel(E::EllCrv, psi::RingElem)
  return Isogeny(E,psi)
end

@doc Markdown.doc"""
    domain(phi::Isogeny) -> EllCrv

Return the domain of the isogeny $\phi$.
"""
function domain(f::Isogeny)
  return f.domain
end

@doc Markdown.doc"""
    domain(phi::Isogeny) -> EllCrv

Return the codomain of the isogeny $\phi$.
"""
function codomain(f::Isogeny)
  return f.codomain
end

@doc Markdown.doc"""
    degree(phi::Isogeny) -> Int

Return the degree of the isogeny $\phi$.
"""
function degree(f::Isogeny)
  return f.degree
end

@doc Markdown.doc"""
    image(phi::Isogeny, P::EllCrvPt) -> EllCrvPt

Return $\phi(P)$.
"""
function image(phi::Isogeny, P::EllCrvPt)
  @assert domain(phi) == parent(P)
  x = P.coordx
  y = P.coordy

  if evaluate(phi.kernel_polynomial,x) == 0
    return infinity(codomain(phi))
  end

  phix = evaluate(evaluate(phi.coordx,y), x)
  phiy = evaluate(evaluate(phi.coordy,y), x)

  return codomain(phi)([phix, phiy])
end

function isogeny_kernel_polynomial(E, psi)

#TODO: Test if polynomial defines a possible kernel, i.e. it needs to be separable, defined over the base field and its roots (in the alg. closure)
#need to be x-coordinates of torsion points on the curve

  R = parent(psi)
  psi = numerator(psi//R(leading_coefficient(psi)))

  psi_G = gcd(division_polynomial_univariate(E, 2, gen(R))[1], psi)
  # force this polynomial to be monic
  psi_G = numerator(psi_G//R(leading_coefficient(psi_G)))

  if (0 != degree(psi_G))
    psi_quo = numerator((psi//psi_G))
      if (degree(psi_quo)!=0)
          throw(DomainError(psi, "In order to use Kohel's algorithm, a kernel of even degree needs to be contained in the two torsion."))
      end
    phi, omega, v, w, n, d = even_kernel_polynomial(E,psi)
  else
    phi, omega, v, w, n, d = odd_kernel_polynomial(E,psi)
  end
  Ry = parent(omega)
  phi = Ry(phi)
  psi = Ry(psi)
  return compute_codomain(E,v,w), phi//psi^2, omega//psi^3, d
end

function even_kernel_polynomial(E, psi_G)
  R = parent(psi_G)
  x = gen(R)
  psi_2tor = division_polynomial_univariate(E,2, x)[1]
  if  mod(psi_2tor,psi_G) != 0
    throw(DomainError(psi_G,"Does not define a finite subgroup of the elliptic curve."))
  end
  n = degree(psi_G)
  d = n+1
  K = base_field(E)
  char = characteristic(K)

  Kxy,y = PolynomialRing(R,"y")

  a1,a2,a3,a4,a6 = a_invars(E)
  b2, b4, b6 = E.b_invars

  if n == 1
    x0 = constant_coefficient(-psi_G)

    # determine y0
    if (2 == characteristic(base_field(E)))
      y0 = sqrt(x0^3 + a2*x0^2 + a4*x0 + a6)
    else
      y0 = -(a1*x0 + a3)//2
      v = (3*x0^2 + 2*a2*x0 + a4 - a1*y0)
      w = x0*v

      phi = (x*psi_G + v)*psi_G
      omega = (y*psi_G^2 - v*(a1*psi_G + (y - y0)))*psi_G
    end
  elseif (3 == n)
    s = coefficients(psi_G)
    s1 = -s[n-1]
    s2 = s[n-2]
    s3 = -s[n-3]

    dpsi_G = derivative(psi_G)
    ddpsi_G = derivative(dpsi_G)

    phi = (dpsi_G^2) + (-2*ddpsi_G + (4*x - s1))*psi_G
    dphi = derivative(phi)

    psi_2 = 2*y + a1*x + a3

    omega = (psi_2*(dphi*psi_G - phi*dpsi_G) - (a1*phi + a3*psi_G)*psi_G)//2

    phi = phi*psi_G
    omega = omega*psi_G

    v = 3*(s1^2 - 2*s2) + b2*s1//2 + 3*b4//2
    w = 3*(s1^3 - 3*s1*s2 + 3*s3) + b2*(s1^2 - 2*s2)//2 + b4*s1//2
  else
     throw(DomainError(psi_G, "Input polynomial must be of degree 1 or 3."))
  end
return phi, omega, v, w, n, d
end

function odd_kernel_polynomial(E, psi)
  n = degree(psi)
  d = 2*n+1

  char = characteristic(base_field(E))

  a1, a2, a3, a4, a6 = a_invars(E)
  b2, b4, b6 = b_invars(E)

  R = parent(psi)
  x = gen(R)

  Rxy,y = PolynomialRing(R,"y")
  psi_2 = 2*y + a1*x + a3

  psicoeffs = coefficients(psi)

  s1 = s2 = s3 = 0

  if (1 <= n)
    s1 = -psicoeffs[n-1]
  end
  if (2 <= n)
    s2 = psicoeffs[n-2]
  end
  if (3 <= n)
    s3 = -psicoeffs[n-3]
  end

  #Variables for computing the codomain
  v = 6*(s1^2 - 2*s2) + b2*s1 + n*b4
  w = 10*(s1^3 - 3*s1*s2 + 3*s3) + 2*b2*(s1^2 - 2*s2) + 3*b4*s1 + n*b6

  dpsi = derivative(psi)
  ddpsi = derivative(dpsi)

  phi = (4*x^3 + b2*x^2 + 2*b4*x + b6)*(dpsi^2 - ddpsi*psi) - (6*x^2 + b2*x + b4)*dpsi*psi + (d*x - 2*s1)*psi^2
  dphi = derivative(phi)

  if (char!=2)
    omega = numerator(dphi*psi*psi_2//2) - phi*dpsi*psi_2 -numerator((a1*phi + a3*psi^2)*psi//2)
  else
  #Characteristic 2 case:
    if (0 < n)
      s1 = -psicoeffs[n-1]
    else
      s1 = 0
    end
    ddpsi = 0
    temp = 1
    for j in (0:n-1)
      ddpsi = ddpsi + binomial(j+2,2)*psicoeffs[j+2]*temp
      temp = x*temp
    end

    dddpsi = 0
    temp = 1
    for j in (0:n-2)
      dddpsi = dddpsi + (3*binomial(j+3,3))*psicoeffs[(j+3)]*temp
      temp = x*temp
    end

    omega = dphi*psi*y - phi*dpsi*psi_2 + ((a1*x + a3)*(psi_2^2)*(ddpsi*dpsi-dddpsi*psi) + (a1*psi_2^2 - 3*(a1*x + a3)*(6*x^2 + b2*x + b4))*dpsi*psi +(a1*x^3 + 3*a3*x^2 + (2*a2*a3 - a1*a4)*x + (a3*a4 - 2*a1*a6))*dpsi^2 +(-(3*a1*x^2 + 6*a3*x + (-a1*a4 + 2*a2*a3)) + (a1*x + a3)*(d*x - 2*s1) )*dpsi*psi + (a1*s1 + a3*n)*psi^2)*psi
  end
  return phi, omega, v, w, n, d
end

function compute_codomain(E::EllCrv, v, w)
  a1, a2, a3, a4, a6 = a_invars(E)
  newa4 = a4 - 5*v
  newa6 = a6 - (a1^2 + 4*a2)*v - 7*w
  return EllipticCurve([a1, a2, a3, newa4, newa6])
end


text:
module GmshTools

using Libdl

deps_jl = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(deps_jl)
  s = """
  Package GmshTools not installed properly.
  Run Pkg.build(\"GridapGmsh\"), restart Julia and try again.
  """
  error(s)
end

include(deps_jl)

if GMSH_FOUND
  include(gmsh_jl)
  # Hack taken from MPI.jl
  function __init__()
    @static if Sys.isunix()
      Libdl.dlopen(gmsh.lib, Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL)
    end
  end
else
    s = """
    Gmsh not found in system paths.
    Install Gmsh or export path to Gmsh and rebuild the project.
    Run Pkg.build(\"GmshTools\"), restart Julia and try again.
    """
    @warn s
end

export
    gmsh,
    match_tuple,
    @addField,
    @gmsh_do,
    @gmsh_open

include("macros.jl")

end


text:
using Test
using UUIDs
using Preferences
set_preferences!(
    UUID("b0ab02a7-8576-43f7-aa76-eaa7c3897c54"),
    "color"=>Dict("fieldname"=>"blue"),
    force=true,
)

using GarishPrint
@test GarishPrint.ColorPreference().fieldname === :blue


text:
using DiffusionMCMC
using ExtensibleMCMC
using Test

const dMCMC = DiffusionMCMC
const eMCMC = ExtensibleMCMC

@testset "schedule.jl" begin
    schedule_params = (
        num_mcmc_iter = 10,
        num_params = 4,
        exclude_params = [(1,3:8), (2,4:2:10), (4,6:10)]
    )
    ∅ = nothing
    t,f = true, false

    schedule = eMCMC.MCMCSchedule(
        schedule_params...;
        start=(
            prev_mcmciter=∅,
            prev_pidx=∅,
            mcmciter=1,
            pidx=1,
            same_layout=(t,t,t)
        ),
        backend=DiffusionMCMCBackend(),
        extra_info=(
            layout_types=(
                (1,1,1),
                (1,2,2),
                (1,1,3),
                (1,2,3),
            ),
            no_blocking=false,
            blocking_layout=:not_needed,
        )
    )

    expected = (
        (∅,∅,1,1,(t,t,t)),  (1,1,1,2,(t,f,f)),  (1,2,1,3,(t,f,f)),  (1,3,1,4,(t,f,t)),
        (1,4,2,1,(t,f,f)),  (2,1,2,2,(t,f,f)),  (2,2,2,3,(t,f,f)),  (2,3,2,4,(t,f,t)),
                            (2,4,3,2,(t,t,f)),  (3,2,3,3,(t,f,f)),  (3,3,3,4,(t,f,t)),
                                                (3,4,4,3,(t,f,t)),  (4,3,4,4,(t,f,t)),
                            (4,4,5,2,(t,t,f)),  (5,2,5,3,(t,f,f)),  (5,3,5,4,(t,f,t)),
                                                (5,4,6,3,(t,f,t)),
                            (6,3,7,2,(t,f,f)),  (7,2,7,3,(t,f,f)),
                                                (7,3,8,3,(t,t,t)),
        (8,3,9,1,(t,t,f)),  (9,1,9,2,(t,f,f)),  (9,2,9,3,(t,f,f)),
        (9,3,10,1,(t,t,f)),                     (10,1,10,3,(t,t,f)),
    )
    for (i,s) in enumerate(schedule)
        @test expected[i] == tuple(s...)
    end
end

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


text:
import Base: copy

export Entity, release!, copy

immutable Entity{T}
    dna::DNA
    body::T
end

function Entity{T}(dna::DNA, body::T)
    Entity{T}(dna, body)
end

@generated function release!{T}(entity::Entity{T})
    if length(methods(release!, (T, ))) > 0
        :( release!(entity.body); release!(entity.dna) )
    else
        :( release!(entity.dna) )
    end
end

@generated function copy{T}(entity::Entity{T})
    if length(methods(copy, (T, ))) > 0
        :( Entity(copy(entity.dna), copy(entity.body)) )
    elseif length(methods(deepcopy, (T, ))) > 0
        :( Entity(copy(entity.dna), deepcopy(entity.body)) )
    else
        :( Entity(copy(entity.dna), entity.body) )
    end
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# ---
# title: step2line
# id: demo_step2line
# cover: assets/step2line.png
# author: Peter Kovesi
# date: 2018-10-26
# ---

# The `step2line()` function generates a phase congruent test image where angle at which the
# congruency occurs is interpolated from 0 at the top of the image to pi/2 at the bottom.
# This produces an interpolation of feature type from step edge to line.  The point being
# that phase congruency at any angle produces a feature and the angle at which the
# congruency occurs defines the feature type. Gradient based edge detectors will only
# correctly mark the step-like feature towards the top of the image and incorrectly mark two
# features towards the bottom of the image whereas phase congruency will correctly mark a
# single feature from top to bottom.  In general, natural images contain a roughly uniform
# distribution of the full continuum of feature types from step to line.

using Images
using ImagePhaseCongruency

img1 = step2line(ampexponent=-1)
## note the softer features
img2 = step2line(ampexponent=-1.5)

## Compute phase congruency on the `step2line` image using default parameters
(pc,) = phasecongmono(step2line(ampexponent = -1))

fimg = imfilter(step2line(ampexponent = -1), KernelFactors.gaussian((2, 2)))
(gx, gy) = imgradients(fimg, KernelFactors.ando3)
∇img = sqrt.(gx.^2 + gy.^2)

mosaicview(
    adjust_histogram(Gray.(img1), LinearStretching()),
    adjust_histogram(img2, LinearStretching()),
    adjust_histogram(pc, LinearStretching()),
    adjust_histogram(∇img, LinearStretching()),
    nrow=2, rowmajor=true
)

# save cover image #src
isdir("assets") || mkdir("assets") #src
save(joinpath("assets", "step2line.png"), adjust_histogram(Gray.(img1), LinearStretching())) #src

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
"Dummy function and forcing default."
@inline zero_func(args...) = 0

"""
    Forcing(Fu, Fv, Fw, FF, FS)

    Forcing(; Fu=zero_func, Fv=zero_func, Fw=zero_func, FT=zero_func, FS=zero_func)

Construct a `Forcing` to specify functions that force `u`, `v`, `w`, `T`, and `S`.
Forcing functions default to `zero_func`, which does nothing.

Forcing functions have the following function signature:
    f(grid::Grid, u::A, v::A, w::A, T::A, S::A, i::Int, j::Int, k::Int)
where A <: AbstractArray, e.g. Array or CuArray.
"""
struct Forcing{Tu,Tv,Tw,TT,TS}
    u::Tu
    v::Tv
    w::Tw
    T::TT
    S::TS
    function Forcing(Fu, Fv, Fw, FT, FS)
        Fu = Fu === nothing ? zero_func : Fu
        Fv = Fv === nothing ? zero_func : Fv
        Fw = Fw === nothing ? zero_func : Fw
        FT = FT === nothing ? zero_func : FT
        FS = FS === nothing ? zero_func : FS
        new{typeof(Fu),typeof(Fv),typeof(Fw),typeof(FT),typeof(FS)}(Fu, Fv, Fw, FT, FS)
    end
end

Forcing(; Fu=nothing, Fv=nothing, Fw=nothing, FT=nothing, FS=nothing) = Forcing(Fu, Fv, Fw, FT, FS)


text:
using JLD
using Distributions
using Random

N = 1000    #population size
time_step = 10.0 # time step length
N_communities = 1 # how many communites there are in the population
community_probs = [1.0] # probability of being in each community
community_contact_rate = [1.0]
density_dependent_fecundity = 0.0007
average_worm_lifespan = 5.7
max_age = 100
initial_worms = 0
initial_miracidia_days = 3
initial_miracidia = 50000*N/1000
init_env_cercariae = 50000*N/1000
worm_stages = 1
egg_production_distribution = "NegBin"
female_factor = 1
male_factor = 1
age_contact_rates = [0.1094983,
0.4363393,
0.4291753,
0.02498713]
c1 = age_contact_rates[1]
c2 = age_contact_rates[2]
c3 = age_contact_rates[3]
c4 = age_contact_rates[4]
ages_for_contacts = [4,9,15,100]
contact_rate_by_age_array = fill(0.09, trunc(Int,(max_age+1)))
mda_adherence = 1
mda_access = 1
female_factor = 1
male_factor = 1
birth_rate = 28*time_step/(1000*365)
human_cercariae_prop = 1
predis_aggregation = 0.24
cercariae_survival = 0.5
miracidia_survival = 0.5
miracidia_maturity = 24
death_prob_by_age = [0.0656, 0.0093, 0.003, 0.0023, 0.0027, 0.0038, 0.0044, 0.0048, 0.0053,
                     0.0065, 0.0088, 0.0106, 0.0144, 0.021, 0.0333, 0.0529, 0.0851, 0.1366, 0.2183, 0.2998 , 0.3698, 1]

ages_for_death = [1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60,
                   65, 70, 75, 80, 85, 90, 95, 100, 110];
r = 0.03
kato_katz_par = 0.87
heavy_burden_threshold = 16
vaccine_effectiveness  = 0.86
drug_effectiveness = 0.86
spec_ages = [7639, 7082, 6524, 5674, 4725, 4147, 3928, 3362,
              2636, 1970, 1468, 1166, 943, 718, 455, 244]
ages_per_index = 5
record_frequency = 1/24;

scenario = "moderate adult";
filename = "eq_runs.jld";
use_kato_katz = 1
# main parameters that we change
# contact_rate = 0.007068083
predis_aggregation = 	0.08489151
max_fecundity = 5
max_fec_contact_rate_product = 1/10
contact_rate = max_fec_contact_rate_product/max_fecundity

M0 = 1000
rate_acquired_immunity = 0
human_larvae_maturity_time = 40
egg_sample_size = 1/100
last_uptake = 0

egg_multiplier = 100

sd_decrease = 1

#
# pars = Parameters(N, time_step, N_communities, community_probs, community_contact_rate,
#         density_dependent_fecundity, average_worm_lifespan,
#         max_age, initial_worms, initial_miracidia, initial_miracidia_days, init_env_cercariae,
#         worm_stages, contact_rate, max_fec_contact_rate_product, max_fecundity, age_contact_rates,
#         ages_for_contacts, contact_rate_by_age_array, mda_adherence, mda_access,  female_factor, male_factor, miracidia_maturity,
#         birth_rate, human_cercariae_prop, predis_aggregation, cercariae_survival, miracidia_survival,
#         death_prob_by_age, ages_for_death, r, vaccine_effectiveness, drug_effectiveness,
#         spec_ages, ages_per_index, record_frequency, use_kato_katz, kato_katz_par, heavy_burden_threshold,
#         rate_acquired_immunity, M0, human_larvae_maturity_time, egg_sample_size)
# pars = make_age_contact_rate_array(pars, scenario, [],[]);


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetNetworkProfile <: SwaggerModel
    healthProbe::Any # spec type: Union{ Nothing, ApiEntityReference } # spec name: healthProbe
    networkInterfaceConfigurations::Any # spec type: Union{ Nothing, Vector{VirtualMachineScaleSetNetworkConfiguration} } # spec name: networkInterfaceConfigurations

    function VirtualMachineScaleSetNetworkProfile(;healthProbe=nothing, networkInterfaceConfigurations=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetNetworkProfile, Symbol("healthProbe"), healthProbe)
        setfield!(o, Symbol("healthProbe"), healthProbe)
        validate_property(VirtualMachineScaleSetNetworkProfile, Symbol("networkInterfaceConfigurations"), networkInterfaceConfigurations)
        setfield!(o, Symbol("networkInterfaceConfigurations"), networkInterfaceConfigurations)
        o
    end
end # type VirtualMachineScaleSetNetworkProfile

const _property_map_VirtualMachineScaleSetNetworkProfile = Dict{Symbol,Symbol}(Symbol("healthProbe")=>Symbol("healthProbe"), Symbol("networkInterfaceConfigurations")=>Symbol("networkInterfaceConfigurations"))
const _property_types_VirtualMachineScaleSetNetworkProfile = Dict{Symbol,String}(Symbol("healthProbe")=>"ApiEntityReference", Symbol("networkInterfaceConfigurations")=>"Vector{VirtualMachineScaleSetNetworkConfiguration}")
Base.propertynames(::Type{ VirtualMachineScaleSetNetworkProfile }) = collect(keys(_property_map_VirtualMachineScaleSetNetworkProfile))
Swagger.property_type(::Type{ VirtualMachineScaleSetNetworkProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetNetworkProfile[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetNetworkProfile }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetNetworkProfile[property_name]

function check_required(o::VirtualMachineScaleSetNetworkProfile)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetNetworkProfile }, name::Symbol, val)
end


text: module PredictMDTestModuleA end


text:
using Test
using Unfold

firbase = firbasis((-1,1),10)

# test optional call
@test firbase.kernel == firbasis((-1,1),10).kernel

@test typeof(firbase.name) <: String
# test basics of basisfunction
@test length(collect(firbase.colnames)) == 21
@test unique(firbase.kernel(1))==[1.0,0.0]

# test length consistency
@test length(firbase.colnames) == length(firbase.times) == size(firbase.kernel(3.1))[2]

# testing the non-sampling rate samples
@test firbase.kernel(0.5)[1:3,1:3] == [0.5 0.0 0.0; 0.5 0.5 0.0; 0.0 0.5 0.5]

boldbase = hrfbasis(2.,name="test")

@test boldbase.name == "test"
@test boldbase.kernel(0) == boldbase.kernel(1)
@test boldbase.kernel(0.1) != boldbase.kernel(1) # testing fractional kernel generation
@test findmax(boldbase.kernel(0.3))[2][1] == 4


text:
"
    Truncate a given state using svd
"
@with_kw struct SvdCut <: Algorithm
    trschemes = [notrunc()]
end

function changebonds(state::Union{FiniteMPS{T},MPSComoving{T}},alg::SvdCut) where T<: GenericMPSTensor{Sp,N} where {Sp,N} # made it work for GenericMPSTensor
    for i in length(state)-1:-1:1
        a = state.AC[i]
        b = _permute_tail(state.AR[i+1]);

        (U,S,V) = tsvd(a*b,trunc=alg.trschemes[mod1(i,end)],alg=TensorKit.SVD());

        state.AC[i] = (U,complex(S));
        state.AC[i+1] = (complex(S),_permute_front(V));
    end

    return state
end

function changebonds(state::InfiniteMPS,alg::SvdCut)
    nAL = copy(state.AL); #not actually left orthonormalized

    for i in 1:length(state)
        (U,S,V) = tsvd(state.CR[i],trunc=alg.trschemes[mod1(i,end)],alg=TensorKit.SVD());
        @tensor nAL[i][-1 -2;-3]:=nAL[i][-1,-2,1]*U[1,-3]
        @tensor nAL[i+1][-1 -2;-3]:=conj(U[1,-1])*nAL[i+1][1,-2,-3]
    end

    return InfiniteMPS(nAL)
end

function changebonds(state,H,alg::SvdCut,pars=params(state,H))
    state = changebonds(state,alg);
    return state,pars;
end


text:
@cdef :PyNumber_Add PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_Subtract PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_Multiply PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_MatrixMultiply PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_FloorDivide PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_TrueDivide PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_Remainder PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_Divmod PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_Power PyPtr (PyPtr, PyPtr, PyPtr)
@cdef :PyNumber_Negative PyPtr (PyPtr,)
@cdef :PyNumber_Positive PyPtr (PyPtr,)
@cdef :PyNumber_Absolute PyPtr (PyPtr,)
@cdef :PyNumber_Invert PyPtr (PyPtr,)
@cdef :PyNumber_Lshift PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_Rshift PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_And PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_Xor PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_Or PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlaceAdd PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlaceSubtract PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlaceMultiply PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlaceMatrixMultiply PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlaceFloorDivide PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlaceTrueDivide PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlaceRemainder PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlacePower PyPtr (PyPtr, PyPtr, PyPtr)
@cdef :PyNumber_InPlaceLshift PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlaceRshift PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlaceAnd PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlaceXor PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_InPlaceOr PyPtr (PyPtr, PyPtr)
@cdef :PyNumber_Long PyPtr (PyPtr,)
@cdef :PyNumber_Float PyPtr (PyPtr,)
@cdef :PyNumber_Index PyPtr (PyPtr,)

for n in [:Number, :Complex, :Real, :Rational, :Integral]
    p = Symbol(:Py, n, :ABC)
    t = Symbol(p, :_Type)
    tr = Symbol(p, :__ref)
    c = Symbol(p, :_Check)
    @eval const $tr = Ref(PyNULL)
    @eval $t(doimport::Bool = true) = begin
        ptr = $tr[]
        isnull(ptr) || return ptr
        a = doimport ? PyImport_ImportModule("numbers") : PyImport_GetModule("numbers")
        isnull(a) && return a
        b = PyObject_GetAttrString(a, $(string(n)))
        Py_DecRef(a)
        isnull(b) && return b
        $tr[] = b
    end
    @eval $c(o) = begin
        t = $t(false)
        isnull(t) && return (PyErr_IsSet() ? Cint(-1) : Cint(0))
        PyObject_IsInstance(o, t)
    end
end


text:
using MathematicalSystems, LazySets, TypedPolynomials
using SparseArrays, LinearAlgebra
using InteractiveUtils: subtypes
using Test

# fix namespace conflicts with LazySets
using MathematicalSystems: LinearMap, AffineMap, ResetMap

const SECOND_ORDER_CTYPES = [SecondOrderAffineContinuousSystem,
                             SecondOrderConstrainedAffineControlContinuousSystem,
                             SecondOrderConstrainedLinearControlContinuousSystem,
                             SecondOrderLinearContinuousSystem]

@testset "Utilities" begin include("utilities.jl") end
@testset "Continuous Systems" begin include("continuous.jl") end
@testset "Discrete Systems" begin include("discrete.jl") end
@testset "Identity Multiple" begin include("identity.jl") end
@testset "Inputs" begin include("inputs.jl") end
@testset "Maps" begin include("maps.jl") end
@testset "Outputs" begin include("outputs.jl") end
@testset "Succesor of discrete system" begin include("successor.jl") end
@testset "Vector field of continuous system" begin include("vector_field.jl") end
@testset "Discretize continuous system" begin include("discretize.jl") end
@testset "@sytem macro" begin include("@system.jl") end
@testset "@ivp macro" begin include("@ivp.jl") end
@testset "Other macros" begin include("macros.jl") end


text:
using CairoMakie, Random, Distributions
Random.seed!(1234)
d = Normal()
b = Binomial(15, 0.7)
n = 350
with_theme(theme_dark()) do
    fig, = scatter(rand(d,n), rand(b, n);  markersize = 12*abs.(rand(d, n)),
        color = tuple.(:dodgerblue, rand(n)), strokewidth = 0.5,
        strokecolor = :white, axis = (; xlabel = "x", ylabel = "y"),
        figure = (; resolution = (600,400)))
    save("scatter_bubbles.png", fig, px_per_unit = 2)
    display(fig)
end


text:
# This is just a bit of code to set up a temporary package environment since not everyone
# will have the same environment set up; this is largely unnecessary when working only on
# your own code.
# See https://bkamins.github.io/julialang/2020/06/28/automatic-project-environments.html
using Pkg: Pkg
# Change to a newly created temporary directory
cd(mktempdir()) do
    Pkg.activate(".")  # activate a new project in the temporary directory
    Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true  # don't update the registry; on Windows it's slow...
    # Add the required packages
    Pkg.add([Pkg.PackageSpec(name="StochasticDiffEq"),
        Pkg.PackageSpec(name="OrdinaryDiffEq"),
        Pkg.PackageSpec(name="StaticArrays"),
        Pkg.PackageSpec(name="Plots"),])
end

# Load required packages
using StochasticDiffEq  # SDE solvers
using OrdinaryDiffEq  # ODE solvers
using StaticArrays  # SVectors for speed
using Plots  # Plotting

# Note on StaticArrays: Ordinary vectors (like Matlab) are created with square brackets,
# e.g., [1, 2, 3], and are good for most purposes. However, since they are variable size the
# compiler cannot optimise as much as it could if they were a fixed size. StaticArrays
# provides fixed-size arrays which allow the compiler to do lots of optimisations and so
# often work much faster for small arrays (typically where there are fewer than 20 elements
# in the array).

"""
    NTMD

A nonlinear tuned-mass-damper system with stochastic forcing of the form

```math
m₁x₁″ + c₁(x₁′ - x₂′) + k₁x₁ + k₂(x₁ - x₂) + μ(x₁ - x₂)^3 = f
m₂x₂″ + c₂(x₂′ - x₁′) + k₂(x₂ - x₁) + μ(x₂ - x₁)^3 = 0
```

where

```math
f = F₀cos(ωt) + ξ(t)
```

and `ξ(t)` is a Brownian motion with variance `σ`.

"""
Base.@kwdef mutable struct NTMD
    m₁::Float64 = 1.0
    c₁::Float64 = 0.02
    k₁::Float64 = 1.0
    m₂::Float64 = 0.5
    c₂::Float64 = 0.05
    k₂::Float64 = 0.95
    μ::Float64 = 0.1
    σ::Float64 = 0.5
    F₀::Float64 = 1.0
    ω::Float64 = 1.0
end

function ntmd_dt(u, p::NTMD, t)
    # Deterministic part
    x₁ = u[1]
    x₁′ = u[2]
    x₂ = u[3]
    x₂′ = u[4]
    f = p.F₀*cos(p.ω*t)
    return SVector(
        x₁′,
        (f - p.c₁*(x₁′ - x₂′) - p.k₁*x₁ - p.k₂*(x₁ - x₂) - p.μ*(x₁ - x₂)^3)/p.m₁,
        x₂′,
        (-p.c₂*(x₂′ - x₁′) - p.k₂*(x₂ - x₁) - p.μ*(x₂ - x₁)^3)/p.m₂,
    )
end

function ntmd_dW(u, p::NTMD, t)
    # Stochastic part
    return SVector(0.0, p.σ/p.m₁, 0.0, 0.0)
end

function ntmd_integrate_sde(u₀, p::NTMD, t)
    # Solve the full SDE
    sde = SDEProblem(ntmd_dt, ntmd_dW, SVector{4, Float64}(u₀), (0.0, Float64(t)), p)
    return solve(sde, SOSRA())
end

function ntmd_integrate_ode(u₀, p::NTMD, t)
    # Solve only the deterministic part of the SDE
    ode = ODEProblem(ntmd_dt, SVector{4, Float64}(u₀), (0.0, Float64(t)), p)
    return solve(ode, Tsit5())
end

function ntmd_example()
    p = NTMD(ω=1.1)  # NTMD with default parameters except ω
    sol = ntmd_integrate_ode(ones(4), p, 2π/p.ω*1000)  # integrate to remove the transient
    u₀ = sol[:, end]  # use the end point for more simulations
    sol = ntmd_integrate_ode(u₀, p, 2π/p.ω*20)  # integrate the deterministic part only (for comparison purposes)
    sol_sde = ntmd_integrate_sde(u₀, p, 2π/p.ω*20)  # integrate the full stochastic differential equation
    plot(sol.t, sol[1, :])
    plot!(sol_sde.t, sol_sde[1, :])
    gui()
    return sol_sde
end

max_stars_repo_name: dawbarton/julia-de-examples

max_issues_repo_name: dawbarton/julia-de-examples

max_forks_repo_name: dawbarton/julia-de-examples


text:
struct Searching <: AbstractPokemonScene
    frame_time::Float64
end
feature_image_name(::Type{Searching}, ::PokemonBattle) = "searching"


text:
#=
There are many attribute-label elements.
The common usage is that 'label' usually be read as annotation-label.

Attribute-labels do not have associated graphics elements. Since <graphics> are
optional for annotation-labels they share the same implementation.

Unknown tags get parsed by unclaimed_label.  Annotation-labels usually have
known tags and dedicated dictonary keys. Pnml-node-elements put unregistered children
into the :labels collection.  It can include annotations and attributes.

Because any tag not present in the tagmap are processed by `unclaimed_label`
it is not necessary to define a parse method unless valididation, documentation,
or additional processing is desired. Some are defined here anyway.
=#

"""
$(TYPEDSIGNATURES)

Return [`Declaration`](@ref) label of 'net' and 'page' nodes.
"""
function parse_declaration(node, pntd; kw...)
    nn = nodename(node)
    nn == "declaration" || error("element name wrong: $nn")
    d = pnml_label_defaults(node, :tag=>Symbol(nn))

    foreach(elements(node)) do child
        @match nodename(child) begin
            "structure" => (d[:structure] = decl_structure(child, pntd; kw...))
            _ => parse_pnml_label_common!(d, child, pntd; kw...)
         end
    end
    Declaration(d[:structure], ObjectCommon(d))
end

# <declaration><structure><declarations><namedsort id="weight" name="Weight">...
# optional,     required,  zero or more
function decl_structure(node, pntd; kw...)
    nn = nodename(node)
    nn == "structure" || error("element name wrong: $nn")
    #TODO warn if more than one?
    declarations = getfirst("declarations", node)
    isnothing(declarations) ? AbstractDeclaration[] :
                            parse_declarations(declarations, pntd; kw...)
end

"""
$(TYPEDSIGNATURES)

Return an Vector{[`AbstractDeclaration`](@ref)} subtype,
"""
function parse_declarations(node, pntd; kw...)
    nn = nodename(node)
    nn == "declarations" || error("element name wrong: $nn")

    v = AbstractDeclaration[]
    foreach(elements(node)) do child
        @match nodename(child) begin
            "namedsort" => push!(v, parse_namedsort(child, pntd; kw...))
            "namedoperator" => push!(v, parse_namedoperator(child, pntd; kw...))
            "variabledecl" => push!(v, parse_variabledecl(child, pntd; kw...))
            _ =>  push!(v, parse_unknowndecl(child, pntd; kw...))
        end
    end
    return v
end

"""
$(TYPEDSIGNATURES)
"""
function parse_namedsort(node, pntd; kw...)
    nn = nodename(node)
    nn == "namedsort" || error("element name wrong: $nn")
    EzXML.haskey(node, "id") || throw(MissingIDException(nn, node))
    EzXML.haskey(node, "name") || throw(MalformedException("$nn missing name attribute", node))

    def = parse_sort(firstelement(node), pntd; kw...)
    NamedSort(register_id!(kw[:reg], node["id"]), node["name"], def)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_namedoperator(node, pntd; kw...)
    nn = nodename(node)
    nn == "namedoperator" || error("element name wrong: $nn")
    EzXML.haskey(node, "id") || throw(MissingIDException(nn, node))
    EzXML.haskey(node, "name") || throw(MalformedException("$(nn) missing name attribute", node))

    @warn "namedoperator under test"
    # <parameter> holds zero or more VariableDeclaration
    def = parse_sort(getfirst("def", node), pntd; kw...)
    parameters = parse_variabledecl.(elements(getfirst("parameter", node)), Ref(pntd); kw...)
    NamedOperator(register_id!(kw[:reg], node["id"]), node["name"], def)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_variabledecl(node, pntd; kw...)
    nn = nodename(node)
    nn == "variabledecl" || error("element name wrong: $nn")
    EzXML.haskey(node, "id") || throw(MissingIDException(nn, node))
    EzXML.haskey(node, "name") || throw(MalformedException("$nn missing name attribute", node))
    # Assert only 1 element
    sort = parse_sort(firstelement(node), pntd; kw...)

    VariableDeclaration(Symbol(node["id"]), node["name"], sort) #TODO register id?
end

"""
$(TYPEDSIGNATURES)
"""
function parse_unknowndecl(node, pntd; kw...)
    nn = nodename(node)
    @info("unknown declaration: $nn")
    EzXML.haskey(node, "id") || throw(MissingIDException(nn, node))
    EzXML.haskey(node, "name") || throw(MalformedException("$nn missing name attribute", node))

    content = anyelement.(elements(node), Ref(pntd); kw...) #TODO Turn children into?
    UnknownDeclaration(Symbol(node["id"]), node["name"], nn, content)
end

#------------------------
"""
$(TYPEDSIGNATURES)

Defines the "sort" of tokens held by the place and semantics of the marking.
NB: The "type" of a place is different from the "type" of a net or "pntd".

"""
function parse_type(node, pntd; kwargs...)
    nn = nodename(node)
    nn == "type" || error("element name wrong: $nn")
    anyelement(node, pntd; kwargs...) #TODO implement sort type
end

"""
$(TYPEDSIGNATURES)

Sorts are found within a <structure> element.
"""
function parse_sort(node, pntd; kw...)
    nn = nodename(node)
    # Builtin
    sort =  nn == "bool" ? anyelement(node, pntd; kw...) :
            nn == "finiteenumeration" ? anyelement(node, pntd; kw...) :
            nn == "finiteintrange" ? anyelement(node, pntd; kw...) :
            nn == "cyclicenumeration"  ? anyelement(node, pntd; kw...) :
            nn == "dot" ? anyelement(node, pntd; kw...) :
            # Also do these.
            nn == "mulitsetsort" ? anyelement(node, pntd; kw...) :
            nn == "productsort" ? anyelement(node, pntd; kw...) :
            nn == "usersort" ? anyelement(node, pntd; kw...) : nothing
    #TODO partition
    isnothing(sort) && error("$nn is not a known sort")
    return sort
end
# BuiltInSort
# MultisetSort
# ProductSort ordered list of sorts
# UserSort

# NamedSort id, name

"""
$(TYPEDSIGNATURES)
"""
function parse_usersort(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "usersort" || error("element name wrong: $nn")
    EzXML.haskey(node, "declaration") || throw(MalformedException("$(nn) missing declaration attribute", node))
    UserSort(anyelement(node, pntd; kwargs...))
end

"""
$(TYPEDSIGNATURES)

There will be no node <term>.
Instead it is the interpertation of the child of some <structure> elements.
"""
function parse_term(node, pntd; kwargs...)
    nn = nodename(node)
    #TODO Validate that it is a kind of term? How? nn == "term" || error("element name wrong: $nn")
    Term(unclaimed_label(node, pntd; kwargs...))
end

#! TODO Variable is one kind of term.
#! TODO Operator is another kind of term.

"""
$(TYPEDSIGNATURES)
"""
function parse_and(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "and" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_arbitraryoperator(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "arbitraryoperator" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_arbitrarysort(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "arbitrarysort" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_bool(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "bool" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_booleanconstant(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "booleanconstant" || error("element name wrong: $nn")
    EzXML.haskey(node, "declaration") || throw(MalformedException("$(nn) missing declaration attribute", node))
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_equality(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "equality" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_imply(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "imply" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_inequality(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "inequality" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_mulitsetsort(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "mulitsetsort" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_not(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "not" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_or(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "or" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_productsort(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "productsort" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_tuple(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "tuple" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_unparsed(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "unparsed" || error("element name wrong: $nn")
    PnmlLabel(node; kwargs...)
end

"""
$(TYPEDSIGNATURES)
"""
function parse_useroperator(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "useroperator" || error("element name wrong: $nn")
    EzXML.haskey(node, "declaration") || throw(MalformedException("$(nn) missing declaration attribute", node))
    UserOperator(Symbol(node["declaration"]))
end

"""
$(TYPEDSIGNATURES)
"""
function parse_variable(node, pntd; kwargs...)
    @debug node
    nn = nodename(node)
    nn == "variable" || error("element name wrong: $nn")
    # The 'primer' UML2 uses variableDecl
    EzXML.haskey(node, "refvariable") || throw(MalformedException("$(nn) missing refvariable attribute", node))
    Variable(Symbol(node["refvariable"]))
end


text:
#All these examples of path are in 4 stages

using Distributions, Random
rng = MersenneTwister(01012019);

"""
	gaussian_path1D()

Returns a '4x1' dimensional array of Gaussian random walk
"""
function gaussian_path1D()
    return vcat(0.0, cumsum(randn(rng, 3, 1), dims = 1)) #4 stages
end

"""
	gaussian_path2D()

Returns a '4x2' dimensional array of Gaussian random walk
"""
function gaussian_path2D()
    gsmatrix = randn(rng, 4, 2) * [1.0 0.0 ; 0.9 0.3] #will create an (dimension x nstages) matrix
    gsmatrix[1,:] .= 0.0
    return cumsum(gsmatrix .+ [1.0 0.0], dims = 1)
end

"""
	running_maximum1D()

Returns a '4x1' dimensional array of Running Maximum process.
"""
function running_maximum1D()
    rmatrix = vcat(0.0, cumsum(randn(rng, 3, 1), dims = 1))
    for i = 2 : 4
        rmatrix[i] = max.(rmatrix[i-1], rmatrix[i])
    end
    return rmatrix
end

"""
	running_maximum2D()

Returns a '4x2' dimensional array of Running Maximum process.
"""
function running_maximum2D()
    rmatrix = vcat(0.0, cumsum(randn(rng, 3, 1), dims = 1))
    rmatrix2D = zeros(4, 2)
    rmatrix2D[:,1] .= vec(rmatrix)
    for j = 2 : 2
        for i = 2 : 4
            rmatrix2D[i,j] = max.(rmatrix[i-1], rmatrix[i])
        end
    end
    return rmatrix2D * [1.0 0.0; 0.9 0.3]
end

"""
	path()

Returns a sample of stock prices following the a simple random random process.
"""
function path()
    return  100 .+ 50 * vcat(0.0, cumsum(randn(rng, 3, 1), dims = 1))
end


text:
using ArrayFire
using Base.Test

#Basic math
a = rand(Float32, 10, 10)
ad = AFArray(a)
@test sumabs2(Array(ad + 2) - (a + 2)) < 1e-6
@test sumabs2(Array(ad .+ 2) - (a .+ 2)) < 1e-6
@test sumabs2(Array(2 + ad) - (2 + a)) < 1e-6
@test sumabs2(Array(2 .+ ad) - (2 .+ a)) < 1e-6
@test sumabs2(Array(ad - 2) - (a - 2)) < 1e-6
@test sumabs2(Array(ad .- 2) - (a .- 2)) < 1e-6
@test sumabs2(Array(2 - ad) - (2 - a)) < 1e-6
@test sumabs2(Array(2 .- ad) - (2 .- a)) < 1e-6
@test sumabs2(Array(ad * 2) - (a * 2)) < 1e-6
@test sumabs2(Array(ad .* 2) - (a .* 2)) < 1e-6
@test sumabs2(Array(2 * ad) - (2 * a)) < 1e-6
@test sumabs2(Array(2 .* ad) - (2 .* a)) < 1e-6
@test sumabs2(Array(ad / 2) - (a / 2)) < 1e-6
@test sumabs2(Array(ad ./ 2) - (a ./ 2)) < 1e-6
@test sumabs2(Array(2 ./ ad) - (2 ./ a)) < 1e-6
@test sumabs2(Array(ad .^ 2) - (a .^ 2)) < 1e-6
@test sumabs2(Array(2 .^ ad) - (2 .^ a)) < 1e-6

#Trig functions
@test sumabs2(Array(sin(ad)) - sin(a)) < 1e-6
@test sumabs2(Array(cos(ad)) - cos(a)) < 1e-6
@test sumabs2(Array(tan(ad)) - tan(a)) < 1e-6
@test sumabs2(Array(sinh(ad)) - sinh(a)) < 1e-6
@test sumabs2(Array(cosh(ad)) - cosh(a)) < 1e-6
@test sumabs2(Array(tanh(ad)) - tanh(a)) < 1e-6
@test sumabs2(Array(atan2(ad, ad)) - atan2(a, a)) < 1e-6
@test sumabs2(Array(atan2(ad, 2)) - atan2(a, 2)) < 1e-6
@test sumabs2(Array(atan2(2, ad)) - atan2(2, a)) < 1e-6

#Measures
@test minimum(a) == minimum(ad)
@test maximum(a) == maximum(ad)
@test sumabs2(Array(maximum(ad,1)) - maximum(a,1)) < 1e-6
@test sumabs2(Array(maximum(ad,2)) - maximum(a,2)) < 1e-6
@test sumabs2(Array(minimum(ad,1)) - minimum(a,1)) < 1e-6
@test sumabs2(Array(minimum(ad,2)) - minimum(a,2)) < 1e-6
@test sumabs2(Array(max(ad,0.5f0)) - max(a,0.5f0)) < 1e-6
@test sumabs2(Array(min(ad,0.5f0)) - min(a,0.5f0)) < 1e-6
@test_approx_eq  mean(ad) mean(a)
@test sumabs2(Array(mean(ad,1)) - mean(a,1)) < 1e-6
@test sumabs2(Array(mean(ad,2)) - mean(a,2)) < 1e-6
@test median(ad) == median(a)
@test sumabs2(Array(median(ad,1)) - median(a,1)) < 1e-6
@test_approx_eq var(ad) var(a)

#Linalg
@test sumabs2(Array(ad') - a') < 1e-6
ld, ud, pd = lu(ad)
l, u, p = lu(a)
@test sumabs(Array(ld) - l) < 1e-5
@test sumabs(Array(ud) - u) < 1e-5
@test sumabs(Array(pd) - p) < 1e-5
@test sumabs2(Array(chol(ad*ad')) - chol(a*a')) < 1e-5
@test sumabs2(Array(ctranspose(chol(ad*ad'))) - ctranspose(chol(a*a'))) < 1e-5
ud, sd, vtd = svd(ad)
u, s, v = svd(a)
@test sumabs(Array(ud) - u) < 1e-4
@test sumabs(Array(sd) - s) < 1e-4
@test sumabs(Array(vtd') - v) < 1e-4

#Complex numbers
@test Array(complex(ad,ad)) == complex(a,a)
@test Array(real(complex(ad,ad))) == real(complex(a,a))
@test Array(imag(complex(a,a))) == imag(complex(a,a))


text:
using Test

include("Transform/dct.jl")
include("BM3DDenoise.jl")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule faust_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("faust")
JLLWrappers.@generate_main_file("faust", UUID("ee02ff01-307c-5f43-b4c2-04ffc4645b0f"))
end  # module faust_jll


text:
import Registrator: RegServer.parse_submission_string

using Test

@test parse_submission_string("register()") == ("register", [], Dict{Symbol,String}())
@test parse_submission_string("register(qux)") == ("register", ["qux"], Dict{Symbol,String}())
@test parse_submission_string("register(qux, baz)") == ("register", ["qux", "baz"], Dict{Symbol,String}())

@test parse_submission_string("approved()") == ("approved", [], Dict{Symbol,String}())

@test parse_submission_string("register(target=qux)") == ("register", String[], Dict(:target=>"qux"))
@test parse_submission_string("register(target=qux, branch=baz)") == ("register", String[], Dict(:branch=>"baz",:target=>"qux"))
@test parse_submission_string("register(qux, baz, target=foo, branch=bar)") == ("register", ["qux", "baz"], Dict(:branch=>"bar",:target=>"foo"))


text:
# Install (Julia 0.6.4)
Pkg.add("JuMP")
Pkg.add("MacroTools")
Pkg.add("Parameters")
Pkg.add("Reexport")
Pkg.add("BenchmarkTools")
Pkg.add("ProgressMeter")
Pkg.add("Plots")
Pkg.add("PyPlot") # Requires matplotlib
Pkg.add("PlotRecipes")

# Solvers (Gurobi requires license)
Pkg.add("Clp")
Pkg.add("Gurobi")

Pkg.clone("https://github.com/martinbiel/TraitDispatch.jl.git")
Pkg.clone("https://github.com/martinbiel/StochasticPrograms.jl.git")
Pkg.clone("https://github.com/martinbiel/HydroModels.jl.git")
Pkg.clone("https://github.com/martinbiel/LShapedSolvers.jl.git")

using JuMP
using BenchmarkTools
using Clp
using Gurobi
using StochasticPrograms
using HydroModels
using LShapedSolvers

max_stars_repo_name: martinbiel/paw-atm-reproducibility

max_issues_repo_name: martinbiel/paw-atm-reproducibility

max_forks_repo_name: martinbiel/paw-atm-reproducibility


text:
module Day17

export get_inputs, get_solution1, get_solution2

using OffsetArrays

## Input getting
function get_inputs()
    test_input1 = test_input2 = read_input(IOBuffer(
        """
        .#.
        ..#
        ###
        """
    ))
    test_output1 = 112
    test_output2 = 848
    data = read_input(joinpath(@__DIR__, "input.txt"))
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

function read_input(io)
    char_matrix = permutedims(reduce(hcat, collect.(readlines(io))), (2,1))
    active_matrix = char_matrix .== '#'
    return active_matrix
end

## Solution functions
# Allocate an OffsetArray for the space the active cells will grow into
function make_space(active_matrix, dims=3, steps=6)
    init_dim_size = (axes(active_matrix)..., ntuple(d->0, dims-2)...)
    dim_size = grow(init_dim_size, steps+1)

    # This will make an OffsetArray with indices that reach back into the negatives
    space = zeros(Bool, dim_size...)
    for ij in CartesianIndices(active_matrix)
        space[Tuple(ij)..., ntuple(d->0, dims-2)...] = active_matrix[ij]
    end
    return space
end

# Adjacent indices of (0,0,0,...)
adjacent_indices(dims) = (
    Tuple(ijk)
    for ijk in CartesianIndices(ntuple(i->-1:1, dims))
    if !all(Tuple(ijk) .== 0)
)
# Get the adjacent cells
get_adjacent(data, idx, adj_inds) = (data[(adj_idx .+ idx)...] for adj_idx in adj_inds)
get_adjacent(data, idx::CartesianIndex, adj_inds) = get_adjacent(data, idx.I, adj_inds)

# Step the simulation by iterating through the bounding_box axes
function step!(space, temp, bounding_box)
    adj_inds = adjacent_indices(ndims(space))
    for idx in CartesianIndices(bounding_box)
        count_adj = count(get_adjacent(space, idx, adj_inds))
        temp[idx] = if space[idx] == true
            count_adj==2 || count_adj==3
        else
            count_adj==3
        end
    end
    space .= temp
end

# Grow a bounding box of array axes
grow(box, n=1) = map(range -> first(range)-n:last(range)+n, box)

# Part 1
function get_solution1(data, dims=3, steps=6)
    space = make_space(data, dims, steps)
    temp = copy(space)
    bounding_box = grow((axes(data)..., ntuple(d->0, dims-2)...))
    for _ in 1:6
        step!(space, temp, bounding_box)
        bounding_box = grow(bounding_box)
    end
    return count(space)
end

# Part 2
get_solution2(data) = get_solution1(data, 4)

end


text:
"""
   `interp_chunk_to_grid_gp_temporal!( flux_out, var_out, chunk_of_spectrum, wavelengths )`
Return spectra interpolated onto a grid of points using linear interpolation.
# Arguments:
- flux_out: (results stored into this array)
- var_out: (results stored into this array)
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
- boost_factor: divide wavelengths by boost_factor
Optional Arguments:
- Filter: Vector with pre-allocated workspace (if length>=1)
# Returns
- flux_out
"""
function interp_chunk_to_shifted_grid_gp_temporal!( flux_out::AA1, var_out::AA2, chunk::AC, grid::AR, boost_factor::Real; use_logx::Bool = true,  use_logy::Bool = false, smooth_factor::Real=1 ) where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, AC<:AbstractChunkOfSpectrum, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    @assert size(flux_out) == size(var_out)
    @assert size(flux_out) == size(grid)
    y = eltype(chunk.flux) == Float64  ? chunk.flux : convert.(Float64,chunk.flux)
    var = eltype(chunk.var) == Float64  ? chunk.var : convert.(Float64,chunk.var)
    # TODO: Could reduce allocations with Lazy Arrays?
    mean_y = mean(y)
    if !use_logy
        y_to_use .= y_to_use./mean_y .- 1.0
        var_to_use .= var_to_use ./ mean_y^2
    else
        var_to_use .= var_to_use ./ y_to_use.^2
        y_to_use .= log.(y_to_use).-log(mean_y)
    end

    flux_out .= TemporalGPInterpolation.predict_mean(chunk.λ./boost_factor, y_to_use, grid, sigmasq_obs = var, use_logx=use_logx, use_logy=use_logy, smooth_factor=smooth_factor )
    # TODO: Update var_out to actually use the right GP
    var_out .= TemporalGPInterpolation.predict_mean(chunk.λ./boost_factor, var, grid, sigmasq_obs = var, use_logx=use_logx, use_logy=use_logy, smooth_factor=smooth_factor )
    if !use_logy
        flux_out .= mean_y.*(flux_out.+1.0)
        var_out .*= mean_y^2
    else
        flux_out .= exp.(flux_out.+log(mean_y))
        var_out .*= flux_out.^2
    end
    return flux_out
end

"""
   `interp_chunk_to_shifted_grid_gp_temporal( chunk_of_spectrum, wavelengths, boost_factor )
Return spectra interpolated onto a grid of points using linear interpolation.
# Arguments:
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
- boost_factor: divide wavelengths by boost_factor
Optional Arguments:
- Filter: Vector with pre-allocated workspace (if length>=1)
# Returns
- flux_out
"""
function interp_chunk_to_shifted_grid_gp_temporal( chunk::AC, grid::AR, boost_factor::Real; use_logx::Bool = true,  use_logy::Bool = false,smooth_factor::Real=1 ) where {  AC<:AbstractChunkOfSpectrum, T2<:Real, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    flux_out = Array{Float64,1}(undef,length(grid))
    var_out = Array{Float64,1}(undef,length(grid))
    interp_chunk_to_grid_gp!(flux_out, var_out, chunk, grid, boost_factor=boost_factor, use_logx=use_logx, use_logy=use_logy, smooth_factor=smooth_factor)
    return (flux=flux_out, var=var_out)
end

"""
   `interp_chunk_to_grid_gp_temporal( chunk_of_spectrum, wavelengths )`
Return spectra interpolated onto a grid of points using linear interpolation.
# Arguments:
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
# Returns
- flux_out
"""
function interp_chunk_to_grid_gp_temporal!( flux_out::AA1, var_out::AA2, chunk::AC, grid::AR ; use_logx::Bool = true,  use_logy::Bool = false, smooth_factor::Real=1 ) where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, AC<:AbstractChunkOfSpectrum, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    @assert size(flux_out) == size(var_out)
    @assert size(flux_out) == size(grid)
    y_to_use = eltype(chunk.flux) == Float64  ? copy(chunk.flux) : convert.(Float64,chunk.flux)
    var_to_use = eltype(chunk.var) == Float64  ? copy(chunk.var) : convert.(Float64,chunk.var)

    # TODO: Could reduce allocations with Lazy Arrays?
    mean_y = mean(y_to_use)
    if !use_logy
        y_to_use .= y_to_use./mean_y .- 1.0
        var_to_use .= var_to_use ./ mean_y^2
    else
        var_to_use .= var_to_use ./ y_to_use.^2
        y_to_use .= log.(y_to_use).-log(mean_y)
    end

    flux_out .= TemporalGPInterpolation.predict_mean(chunk.λ, y_to_use, grid, sigmasq_obs = var_to_use, use_logx=use_logx, use_logy=use_logy, smooth_factor=smooth_factor )
    # TODO: Update var_out to actually use the right GP
    var_out .= TemporalGPInterpolation.predict_mean(chunk.λ, var_to_use, grid, sigmasq_obs = var_to_use, use_logx=use_logx, use_logy=use_logy, smooth_factor=smooth_factor )
    if !use_logy
        flux_out .= mean_y.*(flux_out.+1.0)
        var_out .*= mean_y^2
    else
        flux_out .= exp.(flux_out.+log(mean_y))
        var_out .*= flux_out.^2
    end
    return flux_out
end

"""
   `interp_chunk_to_grid_gp_temporal( chunk_of_spectrum, wavelengths )`
Return spectra interpolated onto a grid of points using linear interpolation.
# Arguments:
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
# Returns
- flux_out
"""
function interp_chunk_to_grid_gp_temporal( chunk::AC, grid::AR ; use_logx::Bool = true,  use_logy::Bool = false, smooth_factor::Real=1 ) where {  AC<:AbstractChunkOfSpectrum, T2<:Real, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    flux_out = Array{Float64,1}(undef,length(grid))
    var_out = Array{Float64,1}(undef,length(grid))
    interp_chunk_to_grid_gp!(flux_out, var_out, chunk, grid, use_logx=use_logx, use_logy=use_logy, smooth_factor=smooth_factor)
    return (flux=flux_out, var=var_out)
end


text:
using Redux
using CImGui

include(joinpath(@__DIR__, "..", "counter", "Counter.jl"))
using .Counter

include("../Renderer.jl")
using .Renderer

# middlewares
function add_logging(dispatch)
    new_dispatch = (s::Store, a::AbstractAction) -> begin
        @info "dispatching action:" a
        result = dispatch(s, a)
        @info "next state:" get_state(s)
        return result
    end
    return new_dispatch
end

function more_logging(dispatch)
    new_dispatch = (s::Store, a::AbstractAction) -> begin
        result = dispatch(s, a)
        @info """------------
                 reducer: $(s.reducer)
                 current_state: $(s.current_state)
                 current_listeners: $(s.current_listeners)
                 next_listeners: $(s.next_listeners)
                 is_dispatching: $(s.is_dispatching)
              """
        return result
    end
    return new_dispatch
end

const store = create_store(Counter.counter, Counter.State(0), [add_logging, more_logging])

function counter_ui(store)
    flag = CImGui.ImGuiWindowFlags_NoTitleBar |
           CImGui.ImGuiWindowFlags_NoResize |
           CImGui.ImGuiWindowFlags_AlwaysAutoResize |
           CImGui.ImGuiWindowFlags_NoSavedSettings |
           CImGui.ImGuiWindowFlags_NoFocusOnAppearing |
           CImGui.ImGuiWindowFlags_NoNav
    CImGui.Begin("Counter", Ref(true), flag)
        spacing = CImGui.GetStyle().ItemInnerSpacing.x
        CImGui.PushButtonRepeat(true)
        CImGui.ArrowButton("##left", CImGui.ImGuiDir_Left) && dispatch!(store, Counter.DECREMENT)
        CImGui.SameLine(0.0, spacing)
        CImGui.ArrowButton("##right", CImGui.ImGuiDir_Right) && dispatch!(store, Counter.INCREMENT)
        CImGui.PopButtonRepeat()
        CImGui.SameLine()
        value = get_state(store).counter
        CImGui.Text("$value")
    CImGui.End()
end

Renderer.render(()->counter_ui(store), width=180, height=50, title="App: Counter")


text:
module Dep442B

using Requires

export check442B

check442B() = true

function link_442A()
    @debug "Loading 442A support into 442B"
    include("support_442A.jl")
end

function __init__()
    @require Dep442A="76238f47-ed95-4e4a-a4d9-95a3fb1630ea" link_442A()
end

end # module


text:
println(0)
println(1)
a = 1
b = a + a
println(a)
println(b)
while a < 1000000
    global a = b + a
    global b = b + a
    println(a)
    println(b)
end

max_stars_repo_name: GingerIndustries/Fibonacci

max_issues_repo_name: GingerIndustries/Fibonacci

max_forks_repo_name: GingerIndustries/Fibonacci


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

"""
    DiscretizationMethod

A method for discretizing geometries into meshes.
"""
abstract type DiscretizationMethod end

"""
    discretize(geometry, method)

Discretize `geometry` with discretization `method`.
"""
function discretize end

"""
    BoundaryDiscretizationMethod

A method for discretizing geometries based on their boundary.
"""
abstract type BoundaryDiscretizationMethod <: DiscretizationMethod end

"""
    discretizewithin(boundary, method)

Discretize geometry within `boundary` with boundary discretization `method`.
"""
function discretizewithin end

discretize(geometry, method::BoundaryDiscretizationMethod) =
  discretizewithin(boundary(geometry), method)

function discretize(polygon::Polygon{Dim,T}, method::BoundaryDiscretizationMethod) where {Dim,T}
  # build bridges in case the polygon has holes,
  # i.e. reduce to a single outer boundary
  chain, dups = bridge(unique(polygon), width=2atol(T))

  # discretize using outer boundary
  mesh = discretizewithin(chain, method)

  if isempty(dups)
    # nothing to be done, return mesh
    mesh
  else
    # remove duplicate vertices
    points = vertices(mesh)
    for (i, j) in dups
      points[i] = centroid(Segment(points[i], points[j]))
    end
    repeated = sort(last.(dups))
    deleteat!(points, repeated)

    # adjust connectivities
    elems = elements(topology(mesh))
    twin  = Dict(reverse.(dups))
    rrep  = reverse(repeated)
    einds = map(elems) do elem
      inds = indices(elem)
      [get(twin, ind, ind) for ind in inds]
    end
    for inds in einds
      for r in rrep
        for i in 1:length(inds)
          inds[i] > r && (inds[i] -= 1)
        end
      end
    end
    connec = connect.(Tuple.(einds))

    # return mesh without duplicates
    SimpleMesh(points, connec)
  end
end

discretize(multi::Multi, method::BoundaryDiscretizationMethod) =
  mapreduce(geom -> discretize(geom, method), merge, multi)

function discretizewithin(chain::Chain{3}, method::BoundaryDiscretizationMethod)
  # collect vertices to get rid of static containers
  points = vertices(chain) |> collect

  # project points on 2D plane of maximum variance
  projected = proj2D(points)

  # discretize within 2D chain with given method
  chain2D = Chain([projected; first(projected)])
  mesh    = discretizewithin(chain2D, method)

  # return mesh with original points
  SimpleMesh(points, topology(mesh))
end

"""
    triangulate(object)

Triangulate `object` of parametric dimension 2 into
triangles using an appropriate triangulation method.
"""
function triangulate end

triangulate(box::Box{2}) = discretize(box, FanTriangulation())

triangulate(tri::Triangle) = discretize(tri, FanTriangulation())

triangulate(quad::Quadrangle) = discretize(quad, FanTriangulation())

triangulate(ngon::Ngon) = discretize(ngon, Dehn1899())

triangulate(poly::PolyArea) = discretize(poly, FIST())

triangulate(multi::Multi) = mapreduce(triangulate, merge, multi)

function triangulate(mesh::Mesh)
  points = vertices(mesh)
  elems  = elements(mesh)
  topo   = topology(mesh)
  connec = elements(topo)

  # initialize vector of global indices
  ginds = Vector{Int}[]

  # triangulate each element and append global indices
  for (e, c) in zip(elems, connec)
    # triangulate single element
    mesh′   = triangulate(e)
    topo′   = topology(mesh′)
    connec′ = elements(topo′)

    # global indices
    inds = indices(c)

    # convert from local to global indices
    einds = [[inds[i] for i in indices(c′)] for c′ in connec′]

    # save global indices
    append!(ginds, einds)
  end

  # new connectivities
  newconnec = connect.(Tuple.(ginds), Triangle)

  SimpleMesh(points, newconnec)
end

triangulate(sphere::Sphere{3}) =
  discretize(sphere, RegularDiscretization(100)) |> triangulate

triangulate(ball::Ball{2}) =
  discretize(ball, RegularDiscretization(100)) |> triangulate

# ----------------
# IMPLEMENTATIONS
# ----------------

include("discretization/fan.jl")
include("discretization/regular.jl")
include("discretization/fist.jl")
include("discretization/dehn.jl")