text:
## Process domains

"Associate the domain with a simpler domain for numerical integration."
has_integration_domain(d) = hasparameterization(d)
"Return the integration domain associated `d`."
integration_domain(d) = parameterdomain(d)
"Return the map from the integration domain of `d` to `d`."
mapfrom_integration_domain(d) = mapfrom_parameterdomain(d)

# Avoid the affine map for intervals and cubes
has_integration_domain(d::Interval) = false
has_integration_domain(d::DomainSets.HyperRectangle) = false

domain_splits(domain) = false
domain_splits(domain::UnionDomain) = true
domain_splits(domain::ProductDomain) = any(map(domain_splits, components(domain)))

domain_split(domain) = (domain,)

domain_split(domain::UnionDomain) = nonoverlapping_domains(domain)
domain_split(domain::ProductDomain) = domain_split(domain, components(domain)...)

function domain_split(domain::ProductDomain, domain1::Domain, domain2::Domain)
    domains1 = domain_split(domain1)
    domains2 = domain_split(domain2)
    [productdomain(d1, d2) for d1 in domains1,d2 in domains2]
end

function domain_split(domain::ProductDomain, domain1::Domain, domain2::Domain, domain3::Domain)
    domains1 = domain_split(domain1)
    domains2 = domain_split(domain2)
    domains3 = domain_split(domain3)
    [productdomain(d1, d2, d3) for d1 in domains1, d2 in domains2, d3 in domains3]
end

function domain_split(domain::ProductDomain, domain1::Domain, domain2::Domain, domain3::Domain, domain4::Domain)
    domains1 = domain_split(domain1)
    domains2 = domain_split(domain2)
    domains3 = domain_split(domain3)
    domains4 = domain_split(domain4)
    [productdomain(d1,d2,d3,d4) for d1 in domains1, d2 in domains2, d3 in domains3, d4 in domains4]
end

"Convert a union of domains into a vector of domains without overlap"
function nonoverlapping_domains(domain::UnionDomain{T}) where {T}
    domains = Vector{Domain{T}}()
    push!(domains, component(domain, 1))
    for i in 2:ncomponents(domain)
        d = component(domain, i)
        for j=1:i-1
            d = d \ domains[j]
        end
        push!(domains, d)
    end
    domains
end

process_domain(qs, integrand, domain, measure, properties...) =
    (integrand, domain, measure, properties...)

function process_domain(qs, integrand, domain::Domain, measure::Lebesgue, properties...)
    if has_integration_domain(domain)
        paramdomain = integration_domain(domain)
        fmap = mapfrom_integration_domain(domain)
        (DomainSets.diffvolume(fmap) * (integrand ∘ fmap), paramdomain, Lebesgue{eltype(paramdomain)}(), properties...)
    else
        (integrand, domain, measure, properties...)
    end
end

integrate_domain(qs, integrand, domain::EmptySpace, measure, properties...) =
    zero_result(integrand)

integrate_domain(qs, integrand, domain::Point, measure, properties...) =
    zero_result(integrand)

# function process_domain(qs, integrand, domain::MappedDomain, measure::Lebesgue, properties...)
#     m = forward_map(domain)
#     @show m
#     process_domain(qs, (t->jacdet(m,t)) * (integrand ∘ m), superdomain(domain), measure, properties...)
# end

function process_domain(qs, integrand, domain::LowerRightTriangle{T}, measure::Lebesgue, properties...) where {T}
    a = domain.a
    b = domain.b
    # For x from a to b, and y from a to x, we set y = a+(x-a)*u, where u goes from 0 to 1. We then have dy = (x-1)du.
    d = UnitInterval{T}()
    square = (a..b) × d
    triangle_map = x -> SVector(x[1],a+x[2]*(x[1]-a))
    # TODO: change the singularity to something sensible: the diagonal is mapped to the right side of the square
    ((x->(x[1]-a)) * (integrand ∘ triangle_map), square, measure)
end

function process_domain(qs, integrand, domain::UpperRightTriangle{T}, measure::Lebesgue, properties...) where {T}
    a = domain.a
    b = domain.b
    # For x from a to b, and y from x to b, we set y = b-(b-x)*u, where u goes from 0 to 1. We then have dy = (b-x)du.
    d = UnitInterval{T}()
    square = (a..b) × d
    triangle_map = x -> SVector(x[1], b-(b-x[1])*(1-x[2]))
    ((x->(b-x[1])) * (integrand ∘ triangle_map), square, measure)
end


text:
function histogram(vector::AbstractVector; use_bin=true)

  return @vlplot(:bar, x={vector, bin=use_bin}, y="count()")

end

function histogram(data, col_name; use_bin=true)

  return data |> @vlplot(:bar, x={col_name, bin=use_bin}, y="count()")

end

ext: jl


text:
const EARTH_RADIUS = 6371.0  # km

struct Layer
    name::AbstractString
    r_min
    r_max
    density::Polynomial
end

struct EarthModel
    layers::Vector{Layer}
end

PREM = EarthModel(
    [
        Layer("Inner Core", 0.0, 1221.5, Polynomial([13.0885, 0.0, -8.8381])),
        Layer("Outer Core", 1221.5, 3480.0, Polynomial([12.5815, -1.2638, -3.6426, -5.5281])),
        Layer("Lower Mantle", 3480.0, 5701.0, Polynomial([7.9565, -6.4761, 5.5283, -3.0807])),
        Layer("Transition Zone 1", 5701.0, 5771.0, Polynomial([5.3197, -1.4836])),
        Layer("Transition Zone 2", 5771.0, 5971.0, Polynomial([11.2494, -8.0298])),
        Layer("Transition Zone 3", 5971.0, 6151.0, Polynomial([7.1089, -3.8045])),
        Layer("LVZ", 6151.0, 6291.0, Polynomial([2.6910, 0.6924])),
        Layer("LID", 6291.0, 6346.6, Polynomial([2.6910, 0.6924])),
        Layer("Crust 1", 6346.6, 6356.0, Polynomial([2.9])),
        Layer("Crust 2", 6356.0, 6368.0, Polynomial([2.6])),
        Layer("Ocean", 6368.0, 6371.0, Polynomial([1.020]))
    ]
)

function (m::EarthModel)(r)
    for layer ∈ m.layers
        if r < layer.r_max && r >= layer.r_min
            return layer.density(r / m.layers[end].r_max)
        end
    end
end

function tracklength(zenith, zposition)
"""
$(SIGNATURES)

Total path length through earth from detector position

# Arguments
- `zenith`: Zenith angle of the path with respect to the detector frame [rad]
- `zposition` Distance below the surface of the Earth (positive value) [km]
"""
    if isapprox(zenith, pi)
        return 2*EARTH_RADIUS - zposition
    elseif isapprox(zenith, 0)
        return zposition
    end

    theta = pi - zenith
    zprime = EARTH_RADIUS - zposition
    sintheta = sin(theta)
    costheta = cos(theta)
    halfcoord = sqrt(EARTH_RADIUS^2 - (zprime * sintheta)^2)
    chordsection = zprime * costheta
    if isapprox(zenith, pi/2)
        return halfcoord
    else
        return halfcoord + chordsection
    end
end

function prempath(zenith::Vector{Float64}, zposition; samples=100, discrete_densities=nothing)
"""
$(SIGNATURES)

# Arguments
- `zenith::Vector{Float64}`: Zenith angles of the paths with respect to the detector frame [rad]
- `zposition::Float64` Distance below the surface of the Earth (positive value) [km]
- `samples` The number of steps with equal distance
- `discrete_densities` List of density values to be used for discretization
"""
    map(z->prempath(z, zposition; samples=samples, discrete_densities=discrete_densities), zenith)
end

function prempath(zenith::T, zposition; samples=100, discrete_densities=nothing) where {T <: Number}
"""
$(SIGNATURES)

# Arguments
- `zenith::Float64`: Zenith angle of the path with respect to the detector frame [rad]
- `zposition::Float64` Distance below the surface of the Earth (positive value) [km]
- `samples` The number of steps with equal distance
"""
    trklen = tracklength(zenith, zposition)
    x = Array(range(0.0; stop=trklen, length=samples))
    sections = (x[2:end] - x[1:end-1])
    total_pathlen = 0.5 * (x[2:end] + x[1:end-1])
    zprime = EARTH_RADIUS - zposition
    radii = map(x -> sqrt(zprime^2 + x^2 + 2*zprime*x*cos(zenith)), total_pathlen)
    densities = PREM.(radii)
    if !isnothing(discrete_densities)
        idx = map(d->searchsortedfirst(discrete_densities,d), densities)
        densities = map(i->discrete_densities[i], idx)
    end
    Path(densities, sections)
end


text:
include("../demo_utilities.jl")
include("../../src/cg_modified.jl")
using BenchmarkTools, Profile

n = 1024
A = Diagonal(collect(LinRange(1e-3, 1, n)).^2)
b = randn(n)
bNorm = norm(b)
function trial()
    cgm = ConjGradModified(A, b)
    ChangeStorageLimit(cgm, 1024)
    TurnOffReorthgonalize(cgm)
    ResNorm = norm(cgm.r)
    iterationCount = 0
    while ResNorm > 1e-10
        ResNorm = cgm()/bNorm
        iterationCount += 1
        if ResNorm == Inf || ResNorm == NaN
            error("Resnorm is inf or nan.")
        end
    end
    println("iterationCount: $iterationCount")
end
@time trial()
display(@benchmark trial())

# function Trial2()
#     cgm = ConjGradModified(A, b)
#     SetStorageLimit(cgm, 1024)
#     # TurnOffReorthgonalize(cgm)
#     ResNorm = norm(cgm.r)
#     iterationCount = 0
#     while ResNorm > 1e-10
#         ResNorm = cgm()/bNorm
#         iterationCount += 1
#         if ResNorm == Inf || ResNorm == NaN
#             error("Resnorm is inf or nan.")
#         end
#     end
#     println("iterationCount: $iterationCount")

# end
# @time Trial2()


text:
module CUDNN

using CUDAapi
using CUDAapi: libraryPropertyType

using CUDAdrv
using CUDAdrv: CUstream

import CUDAnative

using CEnum

using ..CuArrays
using ..CuArrays: active_context, @argout, @workspace
import ..CuArrays.unsafe_free!

import NNlib

const libcudnn = Ref("libcudnn")

# core library
include("libcudnn_common.jl")
include("error.jl")
include("libcudnn.jl")

# low-level wrappers
include("util.jl")
include("base.jl")
include("tensor.jl")
include("conv.jl")
include("pooling.jl")
include("activation.jl")
include("filter.jl")
include("softmax.jl")
include("batchnorm.jl")
include("dropout.jl")
include("rnn.jl")

# high-level integrations
include("nnlib.jl")

include("compat.jl")

const _handles = Dict{CuContext,cudnnHandle_t}()
const _handle = Ref{cudnnHandle_t}(C_NULL)

function handle()
    if _handle[] == C_NULL
        CUDAnative.maybe_initialize("CUDNN")
        _handle[] = get!(_handles, active_context[]) do
            context = active_context[]
            handle = cudnnCreate()
            atexit(()->CUDAdrv.isvalid(context) && cudnnDestroy(handle))
            handle
        end
    end

    return _handle[]
end

end


text:
const Maybe{T} = Union{Nothing, T}

const Status = Dict(
    sfPlaying => :playing,
    sfPaused => :paused,
    sfStopped => :stopped,
)

function showprops(sound::AbstractSound)
    for (name, value) in pairs(getprops(sound))
        println("$name: $value")
    end
end

function timeToSamples(pos, samplerate, channels)
    samples = ((pos * samplerate * channels) + 500000) / 1e6
    return floor(Int, samples)
end


text:
module InterpolationTestUtils

using Test, Interpolations, ForwardDiff, StaticArrays
using Interpolations: degree, itpflag, bounds, lbounds, ubounds

export check_axes, check_inbounds_values, check_oob, can_eval_near_boundaries,
       check_gradient, check_hessian
export MyPair

const failstore = Ref{Any}(nothing)   # stash the inputs to failing tests here

## Property accessors
coefs(itp) = itp.coefs

boundaryconditions(itp::AbstractInterpolation) = boundaryconditions(itpflag(itp))
boundaryconditions(bs::BSpline) = degree(bs)
boundaryconditions(ni::NoInterp) = ni

ndaccessor(x, d) = x
ndaccessor(x::Tuple, d) = x[d]

haspadding(itp) = haspadding(boundaryconditions(itp))
haspadding(::Union{Constant,Linear,NoInterp}) = false
haspadding(::Quadratic{BC}) where BC = haspadding(BC())
haspadding(::Cubic{BC})     where BC = haspadding(BC())
haspadding(::BC) where BC<:Interpolations.BoundaryCondition = !(BC <: Union{Periodic,InPlace,InPlaceQ})
haspadding(bcs::Tuple) = map(haspadding, bcs)
haspadding(::NoInterp) = false

getindexib(itp, i...) = @inbounds itp[i...]
callib(itp, i...)     = @inbounds itp(i...)

⊂(r1::AbstractRange, r2::AbstractRange) = first(r2) < first(r1) < last(r2) && first(r2) < last(r1) < last(r2)

function check_axes(itp, A, isinplace=false)
    @test ndims(itp) == ndims(A)
    axsi, axsA = @inferred(axes(itp)), axes(A)
    szi, szA = @inferred(size(itp)), size(A)
    haspad = haspadding(itp)
    for d = 1:ndims(A)
        if isinplace && ndaccessor(haspad, d)
            @test axsi[d] != axsA[d] && axsi[d] ⊂ axsA[d]
            @test szi[d] < szA[d]
        else
            @test axsi[d] == axsA[d]
            @test szi[d] == szA[d]
        end
    end
    nothing
end

function check_inbounds_values(itp, A)
    i = first(eachindex(itp))
    @test A[i] ≈ @inferred(itp[i]) == @inferred(itp[Tuple(i)...]) ≈ @inferred(itp(i)) ≈ @inferred(itp(float.(Tuple(i))...))
    for i in eachindex(itp)
        @test A[i] ≈ itp[i] == itp[Tuple(i)...] ≈ itp(i) ≈ itp(float.(Tuple(i))...)
    end
    cb = Base.JLOptions().check_bounds == 1
    if ndims(itp) == 1
        for i in eachindex(itp)
            @test itp[i,1] ≈ A[i]   # used in the AbstractArray display infrastructure
            @test_throws BoundsError itp[i,2]
            # if cb
            #     @test_throws BoundsError getindexib(itp, i, 2)
            #     @test_throws BoundsError callib(itp, i, 2)
            # else
            #     @test getindexib(itp, i, 2) ≈ A[i]
            #     @test callib(itp, i, 2) ≈ A[i]
            # end
        end
    end
    nothing
end

function check_oob(itp)
    widen(r) = first(r)-1:last(r)+1
    indsi = axes(itp)
    indsci = CartesianIndices(indsi)
    for i in CartesianIndices(widen.(indsi))
        i ∈ indsci && continue
        @test_throws BoundsError itp[i]
    end
    nothing
end

function can_eval_near_boundaries(itp::AbstractInterpolation{T,1}) where T
    l, u = bounds(itp, 1)
    @test isfinite(itp(l+0.1))
    @test isfinite(itp(u-0.1))
    @test_throws BoundsError itp(l-0.1)
    @test_throws BoundsError itp(u+0.1)
end

function can_eval_near_boundaries(itp::AbstractInterpolation)
    l, u = Float64.(lbounds(itp)), Float64.(ubounds(itp))
    for d = 1:ndims(itp)
        nearl = substitute(l, d, l[d]+0.1)
        # @show summary(itp) nearl
        @test isfinite(itp(nearl...))
        outl = substitute(l, d, l[d]-0.1)
        @test_throws BoundsError itp(outl...)
        nearu = substitute(u, d, u[d]-0.1)
        # @show nearu
        @test isfinite(itp(nearu...))
        outu = substitute(u, d, u[d]+0.1)
        @test_throws BoundsError itp(outu...)
    end
end

function substitute(default::NTuple{N,T}, d::Integer, val::T) where {T,N}
    ntuple(i->ifelse(i==d, val, default[i]), Val(N))
end

# Generate a grid of points [1.0, 1.3333, 1.6667, 2.0, 2.3333, ...] along each coordinate
thirds(axs) = Iterators.product(_thirds(axs...)...)

_thirds(a, axs...) =
    (sort(Float64[a; (first(a):last(a)-1) .+ 1/3; (first(a)+1:last(a)) .- 1/3]), _thirds(axs...)...)
_thirds() = ()

function check_gradient(itp::AbstractInterpolation, gtmp)
    val(x) = itp(Tuple(x)...)
    g!(gstore, x) = ForwardDiff.gradient!(gstore, val, x)
    gtmp2 = similar(gtmp)
    i = first(thirds(axes(itp)))
    @inferred(Interpolations.gradient(itp, i...))
    for i in thirds(axes(itp))
        @test Interpolations.gradient(itp, i...) ≈ g!(gtmp, SVector(i))
        @test Interpolations.gradient!(gtmp2, itp, i...) ≈ gtmp
    end
end

function check_hessian(itp::AbstractInterpolation, htmp)
    val(x) = itp(Tuple(x)...)
    h!(hstore, x) = ForwardDiff.hessian!(hstore, val, x)
    htmp2 = similar(htmp)
    i = first(thirds(axes(itp)))
    @inferred(Interpolations.hessian(itp, i...))
    for i in thirds(axes(itp))
        @test Interpolations.hessian(itp, i...) ≈ h!(htmp, SVector(i))
        @test Interpolations.hessian!(htmp2, itp, i...) ≈ htmp
    end
end

## A type used for multi-valued tests
import Base: +, -, *, /, ≈

struct MyPair{T}
    first::T
    second::T
end

# Here's the interface your type must define
(+)(p1::MyPair, p2::MyPair) = MyPair(p1.first+p2.first, p1.second+p2.second)
(-)(p1::MyPair, p2::MyPair) = MyPair(p1.first-p2.first, p1.second-p2.second)
(*)(n::Number, p::MyPair) = MyPair(n*p.first, n*p.second)
(*)(p::MyPair, n::Number) = n*p
(/)(p::MyPair, n::Number) = MyPair(p.first/n, p.second/n)
Base.zero(::Type{MyPair{T}}) where {T} = MyPair(zero(T),zero(T))
Base.promote_rule(::Type{MyPair{T1}}, ::Type{T2}) where {T1,T2<:Number} = MyPair{promote_type(T1,T2)}
≈(p1::MyPair, p2::MyPair) = (p1.first ≈ p2.first) & (p1.second ≈ p2.second)

end


text:
using Viznet, Test

@testset "mathjax" begin
    @test isa(mathjax(0.5,0.4,0.1,0.2,raw"\frac{x}{2}"), Viznet.MathJax)
    @test isa(mathjax(rand(2),rand(2),rand(2),rand(2),[raw"\frac{x}{2}",raw"\frac{x}{3}"]), Viznet.MathJax)
end


text:
module Solve

export marginalCost!, bertrandSolve!, partialBertrand, collectEquil!
using TypeDef, Helper, Tools, Compute, FOCs, UserFileIO
using NLsolve, Printf

function bertrandSolve!(
  eq  ::Array{equilType,1},
  dat ::datType,
  dim ::dimType,
  env::envType,
  mat ::matType,
  par ::parType,
  set ::setType,
  sim ::simType,
  indx::indxType,
  mask::maskType,
  vInt::Bool=false,
  nEq ::Int64=16)

  if any(eq[1].fail) # Skip if marginal cost solution failed
    @printf("Terminating bertrandSolve!\n")
    return
  end

  if vInt==false
    game2stages!(eq,dat,dim,env,mat,par,set,sim,indx,mask,nEq)
  else
    game1stage!(eq,dat,dim,env,mat,par,set,sim,indx,mask,nEq)
  end
end

function partialBertrand2(
  eq  ::Array{equilType,1},
  dat ::datType,
  dim ::dimType,
  env::envType,
  mat ::matType,
  par ::parType,
  set ::setType,
  sim ::simType,
  indx::indxType,
  mask::maskType,
  vInt::Bool=false)

  peq = deepcopy(eq[14]) # Market without Fight brands
  @. peq.mcost = eq[1].mcost
  for j=1:dim.product
    if j!=10
      peq.price[j] = eq[16].price[j] * peq.pdline[j]
    end
  end

  sim = allocateSim(dim)
  for t=1:dim.qtr
    setupIndex!(sim,indx,peq.pdline,dim,mat,mask,t)
    for j=1:dim.product
      if j!=10
        indx.anyActive[j] = false
      end
    end
    mc = [peq.mcost[1][:,t], peq.mcost[2][:,t]]
    if any(indx.anyActive)
      p0 = peq.price[indx.anyActive,t]
      pj = peq.price[:,t]
      wp = peq.wprice[:,t]
      foc2(x) = foc_retailPartial(x,pj,wp,mc,dim,mat,par,sim,indx,mask)
      stage2 = nlsolve(foc2,p0)
      peq.price[indx.anyActive,t] = stage2.zero
    end
    collectEquil!(peq,t,vInt,peq.mcost,dim,mat,par,sim,indx,mask)
  end

  return peq
end

function partialBertrand(
  eq  ::Array{equilType,1},
  dat ::datType,
  dim ::dimType,
  env::envType,
  mat ::matType,
  par ::parType,
  set ::setType,
  sim ::simType,
  indx::indxType,
  mask::maskType,
  vInt::Bool=false)

  peq = [deepcopy(eq[1]) for i=1:3]
  @. peq[1].pdline = eq[ 3].pdline  # Market without Free
  @. peq[2].pdline = eq[16].pdline  # Market without Free & Fight brands
  @. peq[3].pdline = eq[14].pdline  # Market without Fight brands

  for i=1:3
    for j=1:dim.product
      peq[i].price[j] = peq[i].price[j] * peq[i].pdline[j]
    end
  end

  for i=1:2
    for t=1:dim.qtr
      setupIndex!(sim,indx,peq[i].pdline,dim,mat,mask,t)
      collectEquil!(peq[i],t,vInt,peq[1].mcost,dim,mat,par,sim,indx,mask)
    end
  end

  i=3
  sim = allocateSim(dim)
  for t=1:dim.qtr
    setupIndex!(sim,indx,peq[i].pdline,dim,mat,mask,t)
    for j=1:dim.product
      if j!=10
        indx.anyActive[j] = false
      end
    end
    mc = [peq[1].mcost[1][:,t], peq[1].mcost[2][:,t]]
    if any(indx.anyActive)
      p0 = peq[1].price[indx.anyActive,t]
      pj = peq[1].price[:,t]
      wp = peq[1].wprice[:,t]
      foc2(x) = foc_retailPartial(x,pj,wp,mc,dim,mat,par,sim,indx,mask)
      stage2 = nlsolve(foc2,p0)
      peq[i].price[indx.anyActive,t] = stage2.zero
    end
    collectEquil!(peq[i],t,vInt,peq[1].mcost,dim,mat,par,sim,indx,mask)
  end

  return peq
end

function game2stages!(
  eq  ::Array{equilType,1},
  dat ::datType,
  dim ::dimType,
  env::envType,
  mat ::matType,
  par ::parType,
  set ::setType,
  sim ::simType,
  indx::indxType,
  mask::maskType,
  nEq ::Int64=16)

  for i=1:nEq
    eq[i].fail .= false
    tic = time()
    for t=1:dim.qtr
      setupIndex!(sim,indx,eq[i].pdline,dim,mat,mask,t)
      wprice0 = eq[1].wprice[indx.mvno,t]
      # mcost   = eq[1].mcost[:,t]
      mcost   = [eq[1].mcost[1][:,t], eq[1].mcost[2][:,t]]
      f0      = zeros(sum(indx.anyActive[mask.mvno]))
      fail    = eq[i].fail[t:t]
      foc1(f,x) = foc_whsale!(f,x,sim.p,mcost,fail,dim,mat,par,sim,indx,mask)
      df = OnceDifferentiable(foc1,wprice0,f0,autodiff=:finiteforward)
      stage1 = nlsolve(df,wprice0,iterations=100)
      iter = 1
      while (checkConvg(iter,stage1,fail[1])==false) & (iter<=2)
        fail .= false
        @printf("Re-solving the foc_whsale for %d times....\n",iter+1)
        stage1 = nlsolve(foc1,wprice0+2randn(length(wprice0)),iterations=100)
        iter += 1
      end
      eq[i].fail[t] = (checkConvg(iter,stage1,fail[1])==false)
      if eq[i].fail[t]
        @printf("Solving foc_whsale failed %d times.\n",iter)
        @printf("Exiting bertrandSolve!.\n")
        break
      end
      eq[i].wprice[indx.mvno,t] = stage1.zero

      foc2(x) = foc_retail(x,mcost,eq[i].wprice[:,t],dim,mat,par,sim,indx,mask)
      p0 = view(sim.p,indx.anyActive)
      stage2 = nlsolve(foc2,p0)

      eq[i].price[indx.anyActive,t] = stage2.zero
      collectEquil!(eq[i],t,false,eq[1].mcost,dim,mat,par,sim,indx,mask)
      # @printf("Equil %2d Qtr %2d complete in %.1f seconds.\n",i,t,time()-tic)
    end
    if any(eq[i].fail)
      @printf("Simulating equil %2d failed.\n",i)
      break
    else
      @printf("Equil %2d complete in %.1f seconds.\n",i,time()-tic)
    end
  end # end of i loop
end

function game1stage!(
  eq  ::Array{equilType,1},
  dat ::datType,
  dim ::dimType,
  env::envType,
  mat ::matType,
  par ::parType,
  set ::setType,
  sim ::simType,
  indx::indxType,
  mask::maskType,
  nEq ::Int64=16)

  for i=1:nEq
    eq[i].fail .= false
    # tic = time()
    for t=1:dim.qtr
      setupIndex!(sim,indx,eq[i].pdline,dim,mat,mask,t)
      # mcost   = eq[1].mcost[:,t]
      mcost   = [eq[1].mcost[1][:,t], eq[1].mcost[2][:,t]]
      p0      = sim.p[indx.anyActive]
      f0      = zeros(sum(indx.anyActive))
      fail    = eq[i].fail[t:t]
      foc(f,x) = foc_retailV!(f,x,mcost,dim,mat,par,sim,indx,mask)
      df = OnceDifferentiable(foc,p0,f0,autodiff=:finiteforward)
      stage1 = nlsolve(df,p0,iterations=100)
      iter = 1
      while (checkConvg(iter,stage1,fail[1])==false) & (iter<=2)
        fail .= false
        @printf("Re-solving the foc_retailV for %d times....\n",iter+1)
        stage1 = nlsolve(foc,p0+2randn(length(p0)),iterations=100)
        iter += 1
      end
      eq[i].fail[t] = (checkConvg(iter,stage1,fail[1])==false)
      if eq[i].fail[t]
        @printf("Solving foc_retailV failed %d times.\n",iter)
        @printf("Exiting bertrandSolve!.\n")
        break
      end
      eq[i].price[indx.anyActive,t] = stage1.zero

      D, csplus = demandOnly(sim,eq[i].price[:,t],dim,mat,par,indx)
      eq[i].demand[:,t] = D
      eq[i].csplus[t] = csplus
      eq[i].markup[:,t,2] = eq[i].price[:,t] - mcost[2]
      @. eq[i].profit[:,t,2] = 3D * eq[i].markup[:,t,2]
      # @printf("Equil %2d Qtr %2d complete in %.1f seconds.\n",i,t,time()-tic)
    end
    if any(eq[i].fail)
      @printf("Simulating equil %2d failed.\n",i)
      break
    else
      # @printf("Equil %2d complete in %.1f seconds.\n",i,time()-tic)
    end
  end # end of i loop

end

function marginalCost!(
  eq  ::Array{equilType,1},
  dat ::datType,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  par ::parType,
  set ::setType,
  sim ::simType,
  indx::indxType,
  mask::maskType,
  vInt::Bool=false)

  for t=1:dim.qtr
    if vInt==false
      mc!(eq,t,dim,env,mat,par,sim,indx,mask)
    else
      mcvInt!(eq,t,dim,env,mat,par,sim,indx,mask)
    end
  end

end

function mc!(
  eq  ::Array{equilType,1},
  t   ::Int64,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  par ::parType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  tic = time()
  setupIndex!(sim,indx,eq[1].pdline,dim,mat,mask,t)
  demand!(sim,sim.p,dim,mat,par,indx)
  wprice = whsaleprice(sim.p,dim,mat,sim,indx)

  mcostfree = wprice[mask.free]
  mcostinc0 = wprice[indx.anyActive .& .!mask.free]

  foc(x) = foc_mcost(x,mcostfree,sim.p,wprice,dim,mat,par,sim,indx,mask)
  result = nlsolve(foc,mcostinc0,autodiff=:finiteforward)
  mcost  = mergeMcosts(result.zero,mcostfree,dim,indx,mask)

  fail = eq[1].fail[t:t]
  foc1(f,x) = foc_whsale!(f,x,sim.p,mcost,fail,dim,mat,par,sim,indx,mask)
  Op = nlsolve(foc1,wprice[indx.mvno],autodiff=:finiteforward)
  d = maximum(abs.(Op.zero - wprice[indx.mvno]))
  if (d<1e-6) & converged(Op)
    @printf("Test passed for market %d in %.1f sec.\n",t,time()-tic)
    if env.mc_mvno==0      # MVNO's marginal cost = 0
      # eq[1].mcost[1][:,t] = mcost
      mcostmvno = zeros(sum(indx.mvno))
    elseif env.mc_mvno==1  # MVNO's marginal cost = 10% of wholesale price
      mcostmvno = 0.1*wprice[indx.mvno]
    elseif env.mc_mvno==2  # MVNO's marginal cost = 5% of retail price
      mcostmvno = 0.05*sim.p[indx.mvno]
    end
    eq[1].mcost[2][indx.mno, t] = mcost[2][indx.mno]
    eq[1].mcost[2][indx.mvno,t] = mcostmvno
    eq[1].mcost[1][indx.mvno,t] = mcost[1][indx.mvno] - mcostmvno
    eq[1].wprice[indx.mvno,t]   = wprice[indx.mvno]- mcostmvno
  else
    @printf("Test failed in mc! Markups not matching\n")
    eq[1].fail .= true
    return
  end
end

function mcvInt!(
  eq  ::Array{equilType,1},
  t   ::Int64,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  par ::parType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  tic = time()
  setupIndex!(sim,indx,eq[1].pdline,dim,mat,mask,t)
  demandD1only!(sim,sim.p,dim,mat,par,indx)
  mcost  = [zeros(dim.product) for t=1:2]
  any = indx.anyActive
  Ds = sim.∂D∂p' .* mat.host
  mcost[2][any] = sim.p[any] + Ds[any,any]\sim.D[any]

  # Test validity
  foc(f,x) = foc_retailV!(f,x,mcost,dim,mat,par,sim,indx,mask)
  Op = nlsolve(foc,sim.p[indx.anyActive],autodiff=:finiteforward)
  d = maximum(abs.(Op.zero - sim.p[indx.anyActive]))
  if (d<1e-6) & converged(Op)
    # @printf("Test passed for market %d in %.1f sec.\n",t,time()-tic)
    eq[1].mcost[2][:,t] = mcost[2]
  else
    @printf("Test failed in mcvInt! Markups not matching\n")
    eq[1].fail .= true
    return
  end
end

function collectEquil!(
  eq  ::equilType,
  t   ::Int64,
  vInt::Bool,
  mc  ::Vector{Matrix{Float64}},
  dim ::dimType,
  mat ::matType,
  par ::parType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  p  = eq.price[:,t]
  w  = eq.wprice[:,t]

  D, csplus = demandOnly(sim,p,dim,mat,par,indx)

  for i=1:dim.product
    if indx.anyActive[i]
      if mask.mno[i]
        eq.markup[i,t,1] = 0
        eq.markup[i,t,2] = p[i] - mc[2][i,t]
      elseif mask.mvno[i]
        if vInt==false
          eq.markup[i,t,1] = w[i] - mc[1][i,t]
          eq.markup[i,t,2] = p[i] - w[i] - mc[2][i,t]
        else
          eq.markup[i,t,2] = p[i] - mc[2][i,t]
        end
      end
      eq.demand[i,t]   = D[i]
      eq.profit[i,t,1] = 3D[i] * eq.markup[i,t,1]
      eq.profit[i,t,2] = 3D[i] * eq.markup[i,t,2]
    end
  end
  eq.csplus[t] = csplus

  # mRetail = price - mcost
  # mRetail[mask.mvno] = price[mask.mvno] - wprice[mask.mvno]
  # mWhsale = wprice - mcost
  # mWhsale[mask.mno] .= 0
end

end

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
@testset "Cone" begin

cone = Cone()
@test cone.id == 1
@test cone.h == 1.0
@test GRAYTR.can_intersect(cone)

cone = Cone(2)
@test cone.id == 2
@test cone.h == 1.0

# Create with transformation
cone = Cone(GRAYTR.scaling(2.0, 1.0, 1.0))
@test GRAYTR.obj_bounds(cone) == GRAYTR.BoundingBox(Point3(-1, -1, 0), Point3(1, 1, 1))
@test GRAYTR.world_bounds(cone) == GRAYTR.BoundingBox(Point3(-2, -1, 0), Point3(2, 1, 1))

cone = Cone(2, GRAYTR.scaling(2.0, 1.0, 1.0))
@test cone.id == 2

# Transform after creation
T = GRAYTR.scaling(1.0, 2.0, 1.0)
cone = T(Cone())
@test GRAYTR.obj_bounds(cone) == GRAYTR.BoundingBox(Point3(-1, -1, 0), Point3(1, 1, 1))
@test GRAYTR.world_bounds(cone) == GRAYTR.BoundingBox(Point3(-1, -2, 0), Point3(1, 2, 1))

@testset "Intersections" begin

    # intersectP
    D = Cone(translation(0, 0, -1))

    R = GRAYTR.Ray(Point3(0.5, 0, 1), Vector3(0, 0, -1))
    @test GRAYTR.intersectP(R, D)

    R = GRAYTR.Ray(Point3(0.5, 0, 1), Vector3(0, 0, -1), 10, Inf, 1)
    @test !GRAYTR.intersectP(R, D)

    R = GRAYTR.Ray(Point3(0, 2, 1), Vector3(0, 0, -1))
    @test !GRAYTR.intersectP(R, D)

    R = GRAYTR.Ray(Point3(0.9*sqrt(2), 0.9*sqrt(2), 1), Vector3(0, 0, -1))
    @test !GRAYTR.intersectP(R, D)

    R = GRAYTR.Ray(Point3(0.5, 0, -2), Vector3(0, 0, 1))
    @test GRAYTR.intersectP(R, D)

    R = GRAYTR.Ray(Point3(0.0, 0, -0.2), Vector3(0, 1, 0))
    @test GRAYTR.intersectP(R, D)

    R = GRAYTR.Ray(Point3(2.0, 0, 0.2), Vector3(-1, 0, 0))
    @test !GRAYTR.intersectP(R, D)

    # intersect

    R = GRAYTR.Ray(Point3(0.5, 0, 1), Vector3(0, 0, -1))
    dg, t, e = GRAYTR.shape_intersect(R, D)
    @test dg.p ≈ Point3(0.5, 0, -0.5)
    @test dg.n ≈ normalize(Normal3(1, 0, 1))
    @test t ≈ 1.5

    R = GRAYTR.Ray(Point3(0.5, 0, -3), Vector3(0, 0, 1))
    dg, t, e = GRAYTR.shape_intersect(R, D)
    @test dg.p ≈ Point3(0.5, 0, -0.5)
    @test dg.n ≈ normalize(Normal3(-1.0, 0, -1.0))
    @test t ≈ 2.5

    R = GRAYTR.Ray(Point3(0.9*sqrt(2), 0.9*sqrt(2), 2), Vector3(0, 0, -1))
    dg, t, e = GRAYTR.shape_intersect(R, D)
    @test dg == nothing
    @test isnan(t)
    @test isnan(e)

    R = GRAYTR.Ray(Point3(1, 0, -0.5), Vector3(-1, 0, 0))
    dg, t, e = GRAYTR.shape_intersect(R, D)
    @test dg.p ≈ Point3(0.5, 0.0, -0.5)
    @test dg.n ≈ normalize(Normal3(1.0, 0.0, 1.0))
    @test t ≈ 0.5

    R = GRAYTR.Ray(Point3(0, 0, 1), Vector3(0, 0, -1), 10, Inf, 1)
    dg, t, e = GRAYTR.shape_intersect(R, D)
    @test dg == nothing
    @test isnan(t)
    @test isnan(e)

    R = GRAYTR.Ray(Point3(2.0, 0, 0.2), Vector3(-1, 0, 0))
    dg, t, e = GRAYTR.shape_intersect(R, D)
    @test dg == nothing
    @test isnan(t)
    @test isnan(e)

end # testset "intersections"

@testset "Create from point to point" begin

    cone = GRAYTR.cone_between_points(1, Point3(0,0,0), 1.0, Point3(0,0,1), 0.5)
    @test cone.h ≈ 0.5
    @test !any(isnan, cone.obj_to_world.M)
    @test !any(isnan, cone.world_to_obj.M)
    @test !any(isnan, cone.obj_to_world.MInv)
    @test !any(isnan, cone.world_to_obj.MInv)
    @test cone.obj_to_world(Point3(0,0,0)) ≈ Point3(0,0,0)
    @test cone.obj_to_world(Point3(0,0,0.5)) ≈ Point3(0,0,1)

    cone = GRAYTR.cone_between_points(1, Point3(0,0,0), 0.5, Point3(0,0,1), 1.0)
    @test cone.h ≈ 0.5
    @test !any(isnan, cone.obj_to_world.M)
    @test !any(isnan, cone.world_to_obj.M)
    @test !any(isnan, cone.obj_to_world.MInv)
    @test !any(isnan, cone.world_to_obj.MInv)
    @test cone.obj_to_world(Point3(0,0,0)) ≈ Point3(0,0,1)
    @test isapprox(cone.obj_to_world(Point3(0,0,0.5)), Point3(0,0,0) ; atol=1e-14)

    cone = GRAYTR.cone_between_points(1, Point3(0,0,0), 1.0, Point3(0,0,-1), 0.5)
    @test cone.h ≈ 0.5
    @test !any(isnan, cone.obj_to_world.M)
    @test !any(isnan, cone.world_to_obj.M)
    @test !any(isnan, cone.obj_to_world.MInv)
    @test !any(isnan, cone.world_to_obj.MInv)
    @test cone.obj_to_world(Point3(0,0,0)) ≈ Point3(0,0,0)
    @test isapprox(cone.obj_to_world(Point3(0,0,0.5)), Point3(0,0,-1) ; atol=1e-14)

    cone = GRAYTR.cone_between_points(1, Point3(0,0,0), 0.5, Point3(0,0,-1), 1.0)
    @test cone.h ≈ 0.5
    @test !any(isnan, cone.obj_to_world.M)
    @test !any(isnan, cone.world_to_obj.M)
    @test !any(isnan, cone.obj_to_world.MInv)
    @test !any(isnan, cone.world_to_obj.MInv)
    @test isapprox(cone.obj_to_world(Point3(0,0,0)), Point3(0,0,-1) ; atol=1e-15)
    @test isapprox(cone.obj_to_world(Point3(0,0,0.5)), Point3(0,0,0) ; atol=1e-15)

    P1 = Point3(1, 0, 0)
    P2 = Point3(3, 0, 0)
    r1 = 2.0
    r2 = 1.0

    cone = GRAYTR.cone_between_points(1, P1, r1, P2, r2)
    @test cone.h ≈ 0.5
    @test isapprox(cone.obj_to_world(Point3(0,0,0)), Point3(1,0,0) ; atol=1e-15)
    @test isapprox(cone.obj_to_world(Point3(0,0,0.5)), Point3(3,0,0) ; atol=1e-15)

    # Intersections from Y direction
    R = GRAYTR.Ray(Point3(0.0, -5.0, 0.0), Vector3(0, 1, 0))
    @test !GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(2.0, -5.0, 0.0), Vector3(0, 1, 0))
    @test GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(2.0, -5.0, 1.0), Vector3(0, 1, 0))
    @test GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(2.0, -5.0, -1.0), Vector3(0, 1, 0))
    @test GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(4.0, -5.0, 0.0), Vector3(0, 1, 0))
    @test !GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(2.0, -5.0, 2.0), Vector3(0, 1, 0))
    @test !GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(2.0, -5.0, -2.0), Vector3(0, 1, 0))
    @test !GRAYTR.intersectP(R, cone)

    # Intersections from Z direction
    R = GRAYTR.Ray(Point3(0.0, 0.0, 5.0), Vector3(0, 0, -1))
    @test !GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(2.0, 0.0, 5.0), Vector3(0, 0, -1))
    @test GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(2.0, -1.0, 5.0), Vector3(0, 0, -1))
    @test GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(2.0, -1.0, 5.0), Vector3(0, 0, -1))
    @test GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(4.0, 0.0, 5.0), Vector3(0, 0, -1))
    @test !GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(2.0, 2.0, 5.0), Vector3(0, 0, -1))
    @test !GRAYTR.intersectP(R, cone)

    R = GRAYTR.Ray(Point3(2.0, -2.0, 5.0), Vector3(0, 0, -1))
    @test !GRAYTR.intersectP(R, cone)

end

end # main Cone testset

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using DimensionalData, Test

A = zeros(X(4.0:7.0), Y(10.0:12.0))
di = DimIndices(A)

ci = CartesianIndices(A)
@test val.(collect(di)) == Tuple.(collect(ci))

di[4, 3] == (X(4), Y(3))
@test di[X(1)] == [(Y(1),), (Y(2),), (Y(3),)]
@test map(ds -> A[ds...] + 2, di) == fill(2.0, 4, 3)
@test map(ds -> A[ds...], di[X(At(7.0))]) == [fill(0.0, 4) for i in 1:3]
@test_throws ArgumentError DimIndices(zeros(2, 2))

di = DimKeys(A)
di[4, 3] == (X(At(7.0)), Y(At(12.0)))
@test di[X(1)] == [(Y(At(10.0)),), (Y(At(11.0)),), (Y(At(12.0)),)]
@test map(ds -> A[ds...] + 2, di) == fill(2.0, 4, 3)
@test map(ds -> A[ds...], di[X(At(7.0))]) == [fill(0.0, 4) for i in 1:3]
@test_throws ArgumentError DimIndices(zeros(2, 2))


text:
################################################################################
"""
$(TYPEDEF)
Default Configuration for Metropolis sampler.

# Fields
$(TYPEDFIELDS)
"""
struct DiagnosticsMetropolis{T<:AbstractFloat} <: MCMCKernelDiagnostics
    "Discretization size"
    ϵ::T
    function DiagnosticsMetropolis(ϵ::T) where {T<:AbstractFloat}
        ArgCheck.@argcheck ϵ > 0.0 "Discretization size has to be positive"
        return new{T}(ϵ)
    end
end

function results(
    diagnostics::AbstractVector{DiagnosticsMetropolis{F}},
    Ndigits::Integer,
    quantiles::Vector{T},
) where {F<:AbstractFloat,T<:Real}
    println(
        "Metropolis sampler has average (final) stepsize of ",
        round.(
            mean(diagnostics[iter].ϵ for iter in eachindex(diagnostics)); digits=Ndigits
        ),
        " (",
        round.(diagnostics[end].ϵ; digits=Ndigits),
        ").",
    )
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Show relevant diagnostic results.

# Examples
```julia
```

"""
function generate_showvalues(diagnostics::D) where {D<:DiagnosticsMetropolis}
    return function showvalues()
        return (:Metropolis, "diagnostics"), (:ϵ, diagnostics.ϵ)
    end
end

################################################################################
#export
export DiagnosticsMetropolis, generate_showvalues


text:
const ERROR_MSG = "Pretend Error"

struct MinderUser end

function hear(s::Scene{MinderUser}, ::Val{:die!})
    invite!(s, MinderUser())
    error(ERROR_MSG)
end

function hear(s::Scene{MinderUser}, ::Leave!)
    # Should be called by the minder before Stage has chance
    @test isnothing(Actors.stage_ref(s).state.time_to_leave)
    leave!(s)
end

struct MinderPlay end

function hear(s::Scene{MinderPlay}, ::Genesis!)
    user = invite!(s, MinderUser())

    say(s, user, Val(:die!))

    msg = expect(s, Died!)
    @test msg.corpse.task.exception.msg == ERROR_MSG
    leave!(s)
end

@testset LuvvyTestSet expect=2 "Minders" begin
    play!(MinderPlay())
end


text:
b=3;a=@equ a=$b #:a ≖ 3
@test a==Equation(:a,3)
P=@equ P=Ten($(map(x->pi^x,1:3)),i)
@test floor(P.rhs&@equ(i=1))==3
arr=@equs(e=$ℯ, pi=$pi, M=$(rand(3,2)))
@test arr[1].rhs == ℯ
@test arr[2].rhs == pi

@test Ten(:I,[:i,:i])&@equ(I=[1 0;0 1])==2
@test Ten(:A,[:i,:i])&@equ(A=[:a 0;0 :b])==:a+:b
@test Ten(:A,:i)*Ten(:B,:j)&@equs(A=[1,2,3],B=[3,2,1], j=i)==10
@test Ten(:A,[:j,:i,:i])*Ten(:B,:j)&@equs(A=ones(3,3,3), B=[1,2,3])==18
@test (Alt([:i,:j,:k])*Ten([:a1,:a2,:a3],:j)*Ten([:b1,:b2,:b3],:k)&@equ i=1 )==simplify(:a2*:b3+-1.0*:a3*:b2)

x=@equ x=a*b^sqrt(y)+c/d
@test (x&@equs(a=3, b=2, y=9, c=8, d=4)).rhs==3*2^sqrt(9)+8/4

tri=@equ c^2=a^2+b^2
@test sqrt(tri)==simplify(:c ≖ Equations.Sqrt(:a*:a+:b*:b))

energy=@equ E=m*c^2
c=@equ c=299792458
m=@equ m=3*n
n=@equ n=9
@test (energy&c&m&n).rhs == 3*9*299792458^2

r=(Der(:x^:n,:x)-Der(-0.1*:x^:m,:x)+1/:a*Der(:a*sqrt(:x),:x))&relations["Der"]
@test r==simplify(:n*Equations.Pow(:x,:n+-1)+0.1*:m*Equations.Pow(:x,:m+-1)+0.5*Equations.Pow(:x,-0.5))

relation=@equ Log(:a,:a)=1
@test Log(:e)&relation ==1
@test Log(9,9)&relation ==1
@test Log(:x+:y,:x+:y)&relation ==1

eq=3*:x^2-5*:x+1.5 ≖ 0
meq=eq&quadratic
@test evaluate(eq.lhs,Dict(:x=>meq[1].rhs))==0

rel=@equ Oneable(a)*x*z=y
@test :q*:r&rel ==:y

eq=3*:x^2-5*:x+1.5 ≖ 0
meq=eq&quadratic
@assert evaluate(eq.lhs,Dict(:x=>meq[1].rhs))==0
f1(ex::Expression)=ex[1][1]
f2(fac::Factor)=3*fac
@test (:a+:b)&[f1,f2]==3*:a
f3(eq::Equation)=eq'
f4(eq::Equation)=sqrt(eq)
@test (@equ(a=b^2)&[f3,f4]).lhs==:b

l=U(:l,:meter);t=U(:t,:second);v=l/t
@test v==U(:l*Equations.╱(:t),:meter*Equations.╱(:second))

rule=Der(:a*:x,:x)≖:a #equivalent to Equation(Der(:a*:x,:x),:a)
ex=Der(3*:x,:x)
m=matches(ex,rule)[1] #equivalent to ex&rule
@test m==ex&rule==3

eq=Equation(:x*:z+:y)
@test eq.rhs==0
@test !isempty(matches(eq))
eq=Equation(:x^2,9)
@test matches(eq,Sqrt)[1].rhs==3

meq=matches(:x^2+:a*:x≖0,Div)[1]
try
	evaluate(meq,Dict(:x=>0))
	@test false
catch er
	@test true
end


text:
export Rendering

module Rendering

using ..Ranger

include("colors.jl")
include("state.jl")

include("raster_font.jl")
include("vector_font.jl")
include("render_context.jl")

end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
"""
Plot the raw signal and its envelope
"""
@userplot EnvelopePlot
@recipe function f(e::EnvelopePlot)
    x,y = e.args
    z = abs.(hilbert(y))
    zeroline = zeros(length(x))

    legend --> false
    grid --> false
    linecolor := [:grey :black :black]
    linewidth := [0.5 2 0.7]
    linestyle := [:solid :solid :dash]

    dpi --> 300
    size --> (800,600)

    @series begin
        seriestype := :line
        x, [y z zeroline]
    end
end


text:
struct StepAssertFailure <: Exception
    assertion::String
end

macro expect(ex)
    # Create a human readable form of the expectation.
    if length(ex.args) > 0 && ex.args[1] isa Expr && ex.args[1].head == :globalref
        operator = ex.args[1].args[2]
        readable_ex = copy(ex)
        readable_ex.args[1] = operator
    else
        readable_ex = ex
    end

    quote
        # Check the expectation.
        if !($(esc(ex)))
            throw(StepAssertFailure($(string(readable_ex))))
        end
    end
end


text:
### Reference:
### Heikki Haario, Eero Saksman and Johanna Tamminen
### An Adaptive Metropolis Algorithm
### Bernoulli, 2001, 7 (2), pp 223-242

### Abstract AM state

abstract AMState{F<:VariateForm} <: MHSamplerState{F}

### AM state subtypes

## MuvAMState holds the internal state ("local variables") of the AM sampler for multivariate parameters

type MuvAMState <: AMState{Multivariate}
  pstate::ParameterState{Continuous, Multivariate} # Parameter state used internally by AM
  tune::MCTunerState
  ratio::Real # Acceptance ratio
  lastmean::RealVector
  secondlastmean::RealVector
  C::RealMatrix
  cholC::RealLowerTriangular
  count::Integer

  function MuvAMState(
    pstate::ParameterState{Continuous, Multivariate},
    tune::MCTunerState,
    ratio::Real,
    lastmean::RealVector,
    secondlastmean::RealVector,
    C::RealMatrix,
    cholC::RealLowerTriangular,
    count::Integer
  )
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    new(pstate, tune, ratio, lastmean, secondlastmean, C, cholC, count)
  end
end

MuvAMState(
  pstate::ParameterState{Continuous, Multivariate},
  tune::MCTunerState=BasicMCTune(),
  lastmean::RealVector=Array(eltype(pstate), pstate.size),
  C::RealMatrix=Array(eltype(pstate), pstate.size, pstate.size),
  cholC::RealLowerTriangular=chol(C, Val{:L})
) =
  MuvAMState(pstate, tune, NaN, lastmean, Array(eltype(pstate), pstate.size), C, cholC, 0)

### Adaptive Metropolis (AM) sampler

immutable AM <: MHSampler
  C0::RealMatrix # Initial covariance of proposal according to best prior knowledge
  t0::Integer # Number of initial Monte Carlo steps for which C0 will be used
  sd::Real # Constant that depends only on the dimension of the parameter space
  ε::Real # Non-negative constant with relatively small value

  function AM(C0::RealMatrix, t0::Integer, sd::Real, ε::Real)
    @assert 0 < t0 "Number of initial Monte Carlo steps for which C0 will be used must be positive"
    @assert 0 < sd "Constant sd must be positive"
    @assert 0 <= ε "Constant ε must be non-negative"
    new(C0, t0, sd, ε)
  end
end

AM(C0::RealMatrix, d::Integer; t0::Real=3, ε::Real=0.05) = AM(C0, t0, abs2(2.4)/d, ε)

AM(C0::RealVector, d::Integer; t0::Real=3, ε::Real=0.05) = AM(diagm(C0), d, t0=t0, ε=ε)

AM(C0::Real=1., d::Integer=1; t0::Real=3, ε::Real=0.05) = AM(fill(C0, d), d, t0=t0, ε=ε)

### Initialize AM sampler

## Initialize parameter state

function initialize!{F<:VariateForm}(
  pstate::ParameterState{Continuous, F},
  parameter::Parameter{Continuous, F},
  sampler::AM
)
  parameter.logtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"
end

## Initialize AM state

sampler_state(
  parameter::Parameter{Continuous, Multivariate},
  sampler::AM,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, Multivariate},
  vstate::VariableStateVector
) =
  MuvAMState(
    generate_empty(pstate),
    tuner_state(parameter, sampler, tuner),
    copy(pstate.value),
    copy(sampler.C0),
    RealLowerTriangular(Array(eltype(pstate), pstate.size, pstate.size))
  )

## Reset parameter state

function reset!(
  pstate::ParameterState{Continuous, Multivariate},
  x::RealVector,
  parameter::Parameter{Continuous, Multivariate},
  sampler::AM
)
  pstate.value = copy(x)
  parameter.logtarget!(pstate)
end

## Reset sampler state

function reset!(
  sstate::MuvAMState,
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::AM,
  tuner::MCTuner
)
  reset!(sstate.tune, sampler, tuner)
  sstate.lastmean = copy(pstate.value)
  sstate.C = copy(sampler.C0)
  sstate.cholC = chol(sstate.C, Val{:L})
  sstate.count = 0
end

function covariance!(
  C::RealMatrix,
  lastC::RealMatrix,
  k::Integer,
  x::RealVector,
  lastmean::RealVector,
  secondlastmean::RealVector,
  d::Integer,
  sd::Real,
  ε::Real
)
  C[:, :] = zeros(d, d)
  BLAS.ger!(1.0, x, x, C)
  BLAS.ger!(1.0, -(k+1)*lastmean, lastmean, C)
  BLAS.ger!(1.0, k*secondlastmean, secondlastmean, C)
  C[:, :] = ((k-1)*lastC+sd*(C+ε*eye(d)))/k
end

Base.show(io::IO, sampler::AM) = print(io, "AM sampler: t0 = $(sampler.t0), sd = $(sampler.sd), ε = $(sampler.ε)")

Base.writemime(io::IO, ::MIME"text/plain", sampler::AM) = show(io, sampler)


text:
module ExtrusionPolytopesTests

using Test
using InplaceArrays.Helpers
using InplaceArrays.TensorValues
using InplaceArrays.Arrays
using InplaceArrays.Fields
using InplaceArrays.ReferenceFEs

@test num_faces(VERTEX) == 1
@test num_dims(VERTEX) == 0

p = Polytope(HEX_AXIS, HEX_AXIS)
test_polytope(p,optional=true)

@test get_faces(p,0,0) == [[1], [2], [3], [4]]
@test get_faces(p,1,0) == [[1, 2], [3, 4], [1, 3], [2, 4]]
@test get_faces(p,2,0) == [[1, 2, 3, 4]]
@test get_faces(p,0,1) == [[1, 3], [1, 4], [2, 3], [2, 4]]
@test get_faces(p,1,1) == [[1], [2], [3], [4]]
@test get_faces(p,2,1) == [[1, 2, 3, 4]]
@test get_faces(p,0,2) == [[1], [1], [1], [1]]
@test get_faces(p,1,2) == [[1], [1], [1], [1]]
@test get_faces(p,2,2) == [[1]]
@test get_facedims(p) == [0,0,0,0,1,1,1,1,2]
@test get_offsets(p) == [0,4,8]
@test num_facets(p) == 4
@test num_edges(p) == 4

r = Point{2,Float64}[(1, 0), (1, 0), (0, 1), (0, 1)]
@test get_edge_tangents(p) == r

r = [
  [1, 2, 3, 4], [1, 3, 2, 4], [2, 1, 4, 3], [2, 4, 1, 3],
  [3, 1, 4, 2], [3, 4, 1, 2], [4, 2, 3, 1], [4, 3, 2, 1]]
@test get_vertex_permutations(p) == r

r = Point{2,Float64}[(0, -1), (0, 1), (-1, 0), (1, 0)]
@test get_facet_normals(p) == r

@test num_faces(p) == 9
@test num_dims(p) == 2

p = Polytope(TET_AXIS, TET_AXIS, TET_AXIS)
test_polytope(p,optional=true)

@test num_faces(p) == 15
@test num_dims(p) == 3

x = Point{3,Float64}[(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)]
@test get_vertex_coordinates(p) == x

p = SEGMENT
test_polytope(p,optional=true)
@test get_vertex_coordinates(p) == VectorValue{1,Float64}[(0),(1)]
@test get_edge_tangents(p) == VectorValue{1,Float64}[(1)]
@test get_facet_normals(p) == VectorValue{1,Float64}[(-1),(1)]

p = VERTEX
test_polytope(p,optional=true)
@test get_vertex_coordinates(p) == VectorValue{0,Float64}[()]
@test get_edge_tangents(p) == VectorValue{0,Float64}[]
@test get_facet_normals(p) == VectorValue{0,Float64}[]

test_polytope(TRI,optional=true)
test_polytope(QUAD,optional=true)
test_polytope(TET,optional=true)
test_polytope(HEX,optional=true)

@test num_facets(SEGMENT) == 2
@test num_facets(TRI) == 3
@test num_facets(QUAD) == 4
@test num_facets(TET) == 4
@test num_facets(HEX) == 6

@test num_vertices(PYRAMID) == 5
# TODO There is a problem with PYRAMID
#@test num_edges(PYRAMID) == 8
#@test num_facets(PYRAMID) == 5

end # module


text:
include("networks/Network-Visualization.jl")
include("LiftAlgorithm.jl")

function lift_visualization!(fgp::FixedGraphPlot, lift_alg::LiftAlgorithm)
    _lift_nodecolors!(fgp, lift_alg)
    _lift_edges_visualization!(fgp, lift_alg)
    return fgp
end

function _lift_nodecolors!(fgp::FixedGraphPlot, lift_alg::LiftAlgorithm)
    map!(s -> isempty(s) ? colorant"white" : colorant"snow4", fgp.nodecolors, lift_alg.request_sets)
    for inq_node ∈ lift_alg.inquirer
        fgp.nodecolors[inq_node] = colorant"navajowhite"
    end
    token_pos = lift_alg.token.pos
    if token_pos != 0
        fgp.nodecolors[token_pos] = colorant"crimson"
    end
    return fgp
end

function _lift_edges_visualization!(fgp::FixedGraphPlot, lift_alg::LiftAlgorithm)
    fill!(fgp.edgecolors, colorant"green")
    fill!(fgp.edgelabels, ' ')

    for letter in all_letters(lift_alg.nw)
        edge = fgp.edgeindices[letter.e]
        if letter.m == token
            fgp.edgelabels[edge] = 'T'
            fgp.edgecolors[edge] = colorant"crimson"
        else
            fgp.edgelabels[edge] = 'R'
        end
    end
    return fgp
end

max_stars_repo_name: tomtuamnuq/dist-algorithms

max_issues_repo_name: tomtuamnuq/dist-algorithms

max_forks_repo_name: tomtuamnuq/dist-algorithms


text:
using DataFrames
using Query
using PyPlot
using LaTeXStrings

include("schema.jl")
include("plotting.jl")

const LEGEND_FIGURE_SIZE = (3.2, 3.5)
const PLOT_FIGURE_SIZE = (3, 2.5)
const COLOR_LIKELIHOOD_WEIGHTING_ONE_PARTICLE = "white"
const COLOR_GOLD_STANDARD = "gold"
const ANNOTATION_FONTSIZE = 11
const TITLE_FONTSIZE=13

function annotate_num_clusters_plot(data::DataFrame)

    # draw filled circle for likelihood weighting with one particle
    likelihood_weighting_one_particle = @from i in data begin
            @where ((i.proposal_name == PRIOR_PROPOSAL_NAME) &&
                    (i.num_sweeps == 0) &&
                    (i.num_particles == 1))
            @select {i.num_particles, i.num_clusters_estimate}
            @collect DataFrame
        end
    add_filled_circle(likelihood_weighting_one_particle[COL_NUM_PARTICLES],
                      likelihood_weighting_one_particle[COL_NUM_CLUSTERS_ESTIMATE],
                      COLOR_LIKELIHOOD_WEIGHTING_ONE_PARTICLE)

    # draw filled circle for gold-standard
    gold_standard = @from i in data begin
            @where ((i.proposal_name == OPTIMAL_PROPOSAL_NAME) &&
                    (i.num_sweeps == 4) &&
                    (i.num_particles == 100))
            @select {i.num_particles, i.num_clusters_estimate}
            @collect DataFrame
        end
    add_filled_circle(gold_standard[COL_NUM_PARTICLES],
                      gold_standard[COL_NUM_CLUSTERS_ESTIMATE],
                      COLOR_GOLD_STANDARD)

    # draw arrow and text box for likelihood weighting with one particle
    head = (0.07, 0.69)
    base = (0.35, 0.87)
    draw_arrow(base[1], base[2], head[1] - base[1], head[2] - base[2])
    text(0.35, 0.9, "Appears accurate", ANNOTATION_FONTSIZE)

    # draw horizontal dashed line for gold standard
    xlim = get_xlim(data)
    gold_standard_value = gold_standard[COL_NUM_CLUSTERS_ESTIMATE]
    plt[:plot](xlim, [gold_standard_value, gold_standard_value], "k--")
end

function plot_num_clusters(data::DataFrame, fname::String, to_plot::Vector{Tuple{String,Int}})
    plt[:figure](figsize=PLOT_FIGURE_SIZE)
    for (proposal_name, num_sweeps) in to_plot
        result = @from i in data begin
            @where ((i.proposal_name == proposal_name) &&
                    (i.num_sweeps == num_sweeps))
            @select {i.num_particles, i.num_clusters_estimate, i.num_clusters_stderr}
            @collect DataFrame
        end
        plt[:plot](result[COL_NUM_PARTICLES], result[COL_NUM_CLUSTERS_ESTIMATE])
        lower = result[COL_NUM_CLUSTERS_ESTIMATE] - result[COL_NUM_CLUSTERS_STDERR]
        upper = result[COL_NUM_CLUSTERS_ESTIMATE] + result[COL_NUM_CLUSTERS_STDERR]
        plt[:fill_between](result[COL_NUM_PARTICLES], lower, upper, alpha=0.3)
    end

    annotate_num_clusters_plot(data)

    ax = plt[:gca]()
    ax[:set_xscale]("log")
    ax[:set_xlim](get_xlim(data))
    plt[:ylabel]("Average number\nof clusters")
    plt[:xlabel]("Number of particles")
    plt[:title]("Heuristic diagnostic", fontsize=TITLE_FONTSIZE)
    plt[:tight_layout](pad=0)
    plt[:savefig](fname)
end

function annotate_aide_plot(data::DataFrame)

    # draw filled circle for likelihood weighting with one particle
    likelihood_weighting_one_particle = @from i in data begin
            @where ((i.proposal_name == PRIOR_PROPOSAL_NAME) &&
                    (i.num_sweeps == 0) &&
                    (i.num_particles == 1))
            @select {i.num_particles, i.aide_estimate}
            @collect DataFrame
        end
    add_filled_circle(likelihood_weighting_one_particle[COL_NUM_PARTICLES],
                      likelihood_weighting_one_particle[COL_AIDE_ESTIMATE],
                      COLOR_LIKELIHOOD_WEIGHTING_ONE_PARTICLE)

    # draw arrow and text box for likelihood weighting with one particle
    head = (0.07, 0.93)
    base = (0.23, 0.9)
    draw_arrow(base[1], base[2], head[1] - base[1], head[2] - base[2])
    text(0.25, 0.7, "Likelihood weighting\nwith 1 particle appears\nleast accurate", ANNOTATION_FONTSIZE)

    # draw horizontal dashed line for zero KL divergence
    xlim = get_xlim(data)
    plt[:plot](xlim, [0., 0.], "k--")
end

function plot_aide(data::DataFrame, fname::String, to_plot::Vector{Tuple{String,Int}})
    plt[:figure](figsize=PLOT_FIGURE_SIZE)
    for (proposal_name, num_sweeps) in to_plot
        result = @from i in data begin
            @where ((i.proposal_name == proposal_name) &&
                    (i.num_sweeps == num_sweeps))
            @select {i.num_particles, i.aide_estimate, i.aide_stderr}
            @collect DataFrame
        end
        plt[:plot](result[COL_NUM_PARTICLES], result[COL_AIDE_ESTIMATE])
        lower = result[COL_AIDE_ESTIMATE] - result[COL_AIDE_STDERR]
        upper = result[COL_AIDE_ESTIMATE] + result[COL_AIDE_STDERR]
        plt[:fill_between](result[COL_NUM_PARTICLES], lower, upper, alpha=0.3)
    end

    annotate_aide_plot(data)

    ax = plt[:gca]()
    ax[:set_xscale]("log")
    ax[:set_xlim](get_xlim(data))
    plt[:ylabel]("nats")
    plt[:xlabel]("Number of particles")
    plt[:title]("AIDE estimates", fontsize=TITLE_FONTSIZE)
    plt[:tight_layout](pad=0)
    plt[:savefig](fname)
end

function plot_legend(data::DataFrame, fname::String)
    plt[:figure](figsize=LEGEND_FIGURE_SIZE)
    plt[:plot]([], [], label="SMC, prior proposal\n0 rejuvenation sweeps")
    plt[:plot]([], [], label="SMC, optimal proposal\n0 rejuvenation sweeps")
    plt[:plot]([], [], label="SMC, optimal proposal\n4 rejuvenation sweeps")
    plt[:scatter]([], [], color="gold", s=70, edgecolor="black", label="Gold-standard")
    plt[:scatter]([], [], color="white", s=70, edgecolor="black", label="Likelihood-weighting\n(1 particle)")
    plt[:legend](fontsize=13, ncol=1, frameon=false, borderpad=0, borderaxespad=0, loc="lower left", labelspacing=2)
    plt[:axis]("off")
    plt[:tight_layout](pad=0)
    plt[:savefig](fname)
end

data_dir = "data"

plot_dir = "plots"
if !Base.Filesystem.isdir(plot_dir)
    Base.Filesystem.mkdir(plot_dir)
end

data = readtable("$data_dir/sampler_data.csv")

# which proposals, and which number of sweeps to plot, as a function of number
# of particles
to_plot = [
    (PRIOR_PROPOSAL_NAME, 0),
    (OPTIMAL_PROPOSAL_NAME, 0),
    (OPTIMAL_PROPOSAL_NAME, 4)]
plot_aide(data, "$plot_dir/aide.pdf", to_plot)
plot_num_clusters(data, "$plot_dir/num_clusters.pdf", to_plot)
plot_legend(data, "$plot_dir/legend.pdf")


text:
using Documenter, RealNeuralNetworks
using RealNeuralNetworks.Neurons
using RealNeuralNetworks.Neurons.Segments
using RealNeuralNetworks.Neurons.Segments.Synapses
using RealNeuralNetworks.SWCs

makedocs(
	modules=[RealNeuralNetworks, Neurons, Segments, Synapses, SWCs],
	sitename="RealNeuralNetworks.jl",
	authors="Jingpeng Wu",
	format=:html,
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Getting Started" => "man/getting_started.md",
            "man/examples.md",
            "man/references.md",
        ],
        "Library" => Any[
            "Public" => "lib/public.md",
            hide("Internals" => "lib/internals.md", Any[
                "lib/internals/Manifests.md",
                "lib/internals/NodeNets.md",
                "lib/internals/PointArrays.md",
            ])
        ],
        "man/contributing.md",
    ],
    # use clean URLs, unless built as a "local" build
    html_prettyurls = !("local" in ARGS),
    html_canonical = "https://seung-lab.github.io/RealNeuralNetworks.jl/latest",
)

deploydocs(
	repo="github.com/seung-lab/RealNeuralNetworks.jl",
	target="build",
	julia="0.7",
	deps=nothing,
	make=nothing
)


text:
using Images

"""
    _substitution(image, keys, type; path_for_result="./encrypted.png")

Performs substitution encryption/decryption on a given image with the given keys.

See [`substitution_encryption`](@ref) and [`substitution_decryption`](@ref) for more details.

# Arguments
- `image::Array{RGB{N0f8},2}`: A loaded image.
- `keys::Array{Int64, 1}`: Keys for encryption.
- `type::Symbol`: Can be `:encrypt` or `:decrypt`.
- `path_for_result::String`: The path for storing the encrypted image.
- `inplace::Boolean`: Perform substitution on the provided image.
"""
function _substitution(
    image::Union{String,Array{RGB{N0f8},2}},
    keys::Vector{Int64},
    type::Symbol;
    path_for_result::String="./encrypted.png",
    inplace=false,
)

    if typeof(image) == String
        image = load(image)
    end

    # Generating dimensions of the image
    height = size(image)[1]
    width = size(image)[2]

    if length(keys) != height * width
        throw(ArgumentError("Number of keys must be equal to height * width of image."))
    end

    # generate a copy if not inplace
    ~inplace && (image = copy(image))

    if type == :encrypt
        @info "ENCRYPTING"
    else
        @info "DECRYPTING"
    end

    # reshape keys for broadcasting
    keys = reshape(keys, height, width)

    # substitute all pixels in one go
    @. image = _substitute_pixel(image, keys)

    if type == :encrypt
        @info "ENCRYPTED"
    else
        @info "DECRYPTED"
    end

    save(path_for_result, image)
    image
end

"""
    _substitute_pixel(pixel::RGB, key::Int64)

Returns the pixel after XORing the R, G, and B values with the key.
Specifically developed to return an `Array` (or the complete image)
of XORed RGB values in one go.

See [`_substitution`](@ref) for more details.

# Arguments
- `pixel::RGB`: Pixel value with r, g, and b components.
- `key::Int64`: The key.

# Returns
- `pixel::RGB`: Substituted pixel.
"""
_substitute_pixel(pixel::RGB, key::Int64) = RGB(
    (trunc(Int, pixel.r * 255) ⊻ key) / 255,
    (trunc(Int, pixel.g * 255) ⊻ key) / 255,
    (trunc(Int, pixel.b * 255) ⊻ key) / 255
)

"""
    substitution_encryption(image, keys; path_for_result="./encrypted.png")

Performs substitution encryption on a given image with the given keys.

# Algorithm
Iterates simulataneously over each pixel and key, and XORs the pixel value
(all R, G, and B) with the given key. Hence, the order of the keys matter.

# Arguments
- `image::Array{RGB{N0f8},2}`: A loaded image.
- `keys::Array{Int64, 1}`: Keys for encryption.
- `path_for_result::String`: The path for storing the encrypted image.

# Returns
- `image::Array{RGB{N0f8}, 2}`: Encrypted image.

# Example
```jldoctest
julia> using Images

julia> img = load("../test_images/camera.jfif")
225×225 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 ⋮                             ⋱
 RGB{N0f8}(0.4,0.4,0.4)           RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.447,0.447,0.447)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.427,0.427,0.427)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.451,0.451,0.451)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.51,0.51,0.51)     …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.537,0.537,0.537)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.412,0.412,0.412)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.149,0.149,0.149)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.0,0.0,0.0)           RGB{N0f8}(0.0,0.0,0.0)

julia> height, width = size(img)
(225, 225)

julia> keys = logistic_key(0.01, 3.97, height * width)
50625-element Vector{Int64}:
   0
  44
   7
  26
  14
 224
  16
 250
 162
 211
   ⋮
  72
 255
  80
  84
  21
 154
 197
  82
 147

julia> substitution_encryption(img, keys)
[ Info: ENCRYPTING
[ Info: ENCRYPTED
225×225 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.839,0.839,0.839)
 RGB{N0f8}(0.718,0.718,0.718)     RGB{N0f8}(0.361,0.361,0.361)
 RGB{N0f8}(0.612,0.612,0.612)     RGB{N0f8}(0.055,0.055,0.055)
 RGB{N0f8}(0.506,0.506,0.506)     RGB{N0f8}(0.58,0.58,0.58)
 RGB{N0f8}(0.584,0.584,0.584)     RGB{N0f8}(0.631,0.631,0.631)
 RGB{N0f8}(0.482,0.482,0.482)  …  RGB{N0f8}(0.639,0.639,0.639)
 RGB{N0f8}(0.545,0.545,0.545)     RGB{N0f8}(0.196,0.196,0.196)
 RGB{N0f8}(0.38,0.38,0.38)        RGB{N0f8}(0.047,0.047,0.047)
 RGB{N0f8}(0.224,0.224,0.224)     RGB{N0f8}(0.847,0.847,0.847)
 RGB{N0f8}(0.282,0.282,0.282)     RGB{N0f8}(0.827,0.827,0.827)
 ⋮                             ⋱
 RGB{N0f8}(0.573,0.573,0.573)     RGB{N0f8}(0.282,0.282,0.282)
 RGB{N0f8}(0.361,0.361,0.361)     RGB{N0f8}(1.0,1.0,1.0)
 RGB{N0f8}(0.847,0.847,0.847)     RGB{N0f8}(0.314,0.314,0.314)
 RGB{N0f8}(0.349,0.349,0.349)     RGB{N0f8}(0.329,0.329,0.329)
 RGB{N0f8}(0.984,0.984,0.984)  …  RGB{N0f8}(0.082,0.082,0.082)
 RGB{N0f8}(0.031,0.031,0.031)     RGB{N0f8}(0.604,0.604,0.604)
 RGB{N0f8}(0.071,0.071,0.071)     RGB{N0f8}(0.773,0.773,0.773)
 RGB{N0f8}(0.765,0.765,0.765)     RGB{N0f8}(0.322,0.322,0.322)
 RGB{N0f8}(0.902,0.902,0.902)     RGB{N0f8}(0.576,0.576,0.576)
```
"""
substitution_encryption(
    image::Array{RGB{N0f8},2},
    keys::Vector{Int64};
    path_for_result::String="./encrypted.png"
) = _substitution(image, keys, :encrypt; path_for_result=path_for_result)

"""
    substitution_encryption!(image, keys; path_for_result="./encrypted.png")

Performs substitution encryption on a given image with the given keys.

# Algorithm
Iterates simulataneously over each pixel and key, and XORs the pixel value
(all R, G, and B) with the given key. Hence, the order of the keys matter.

# Arguments
- `image::Array{RGB{N0f8},2}`: A loaded image.
- `keys::Array{Int64, 1}`: Keys for encryption.
- `path_for_result::String`: The path for storing the encrypted image.

# Returns
- `image::Array{RGB{N0f8}, 2}`: Encrypted image.

# Example
```jldoctest
julia> using Images

julia> img = load("../test_images/camera.jfif")
225×225 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 ⋮                             ⋱
 RGB{N0f8}(0.4,0.4,0.4)           RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.447,0.447,0.447)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.427,0.427,0.427)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.451,0.451,0.451)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.51,0.51,0.51)     …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.537,0.537,0.537)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.412,0.412,0.412)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.149,0.149,0.149)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.0,0.0,0.0)           RGB{N0f8}(0.0,0.0,0.0)

julia> height, width = size(img)
(225, 225)

julia> keys = logistic_key(0.01, 3.97, height * width)
50625-element Vector{Int64}:
   0
  44
   7
  26
  14
 224
  16
 250
 162
 211
   ⋮
  72
 255
  80
  84
  21
 154
 197
  82
 147

julia> substitution_encryption!(img, keys)
[ Info: ENCRYPTING
[ Info: ENCRYPTED
225×225 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.839,0.839,0.839)
 RGB{N0f8}(0.718,0.718,0.718)     RGB{N0f8}(0.361,0.361,0.361)
 RGB{N0f8}(0.612,0.612,0.612)     RGB{N0f8}(0.055,0.055,0.055)
 RGB{N0f8}(0.506,0.506,0.506)     RGB{N0f8}(0.58,0.58,0.58)
 RGB{N0f8}(0.584,0.584,0.584)     RGB{N0f8}(0.631,0.631,0.631)
 RGB{N0f8}(0.482,0.482,0.482)  …  RGB{N0f8}(0.639,0.639,0.639)
 RGB{N0f8}(0.545,0.545,0.545)     RGB{N0f8}(0.196,0.196,0.196)
 RGB{N0f8}(0.38,0.38,0.38)        RGB{N0f8}(0.047,0.047,0.047)
 RGB{N0f8}(0.224,0.224,0.224)     RGB{N0f8}(0.847,0.847,0.847)
 RGB{N0f8}(0.282,0.282,0.282)     RGB{N0f8}(0.827,0.827,0.827)
 ⋮                             ⋱
 RGB{N0f8}(0.573,0.573,0.573)     RGB{N0f8}(0.282,0.282,0.282)
 RGB{N0f8}(0.361,0.361,0.361)     RGB{N0f8}(1.0,1.0,1.0)
 RGB{N0f8}(0.847,0.847,0.847)     RGB{N0f8}(0.314,0.314,0.314)
 RGB{N0f8}(0.349,0.349,0.349)     RGB{N0f8}(0.329,0.329,0.329)
 RGB{N0f8}(0.984,0.984,0.984)  …  RGB{N0f8}(0.082,0.082,0.082)
 RGB{N0f8}(0.031,0.031,0.031)     RGB{N0f8}(0.604,0.604,0.604)
 RGB{N0f8}(0.071,0.071,0.071)     RGB{N0f8}(0.773,0.773,0.773)
 RGB{N0f8}(0.765,0.765,0.765)     RGB{N0f8}(0.322,0.322,0.322)
 RGB{N0f8}(0.902,0.902,0.902)     RGB{N0f8}(0.576,0.576,0.576)

julia> img  # inplace
225×225 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.839,0.839,0.839)
 RGB{N0f8}(0.718,0.718,0.718)     RGB{N0f8}(0.361,0.361,0.361)
 RGB{N0f8}(0.612,0.612,0.612)     RGB{N0f8}(0.055,0.055,0.055)
 RGB{N0f8}(0.506,0.506,0.506)     RGB{N0f8}(0.58,0.58,0.58)
 RGB{N0f8}(0.584,0.584,0.584)     RGB{N0f8}(0.631,0.631,0.631)
 RGB{N0f8}(0.482,0.482,0.482)  …  RGB{N0f8}(0.639,0.639,0.639)
 RGB{N0f8}(0.545,0.545,0.545)     RGB{N0f8}(0.196,0.196,0.196)
 RGB{N0f8}(0.38,0.38,0.38)        RGB{N0f8}(0.047,0.047,0.047)
 RGB{N0f8}(0.224,0.224,0.224)     RGB{N0f8}(0.847,0.847,0.847)
 RGB{N0f8}(0.282,0.282,0.282)     RGB{N0f8}(0.827,0.827,0.827)
 ⋮                             ⋱
 RGB{N0f8}(0.573,0.573,0.573)     RGB{N0f8}(0.282,0.282,0.282)
 RGB{N0f8}(0.361,0.361,0.361)     RGB{N0f8}(1.0,1.0,1.0)
 RGB{N0f8}(0.847,0.847,0.847)     RGB{N0f8}(0.314,0.314,0.314)
 RGB{N0f8}(0.349,0.349,0.349)     RGB{N0f8}(0.329,0.329,0.329)
 RGB{N0f8}(0.984,0.984,0.984)  …  RGB{N0f8}(0.082,0.082,0.082)
 RGB{N0f8}(0.031,0.031,0.031)     RGB{N0f8}(0.604,0.604,0.604)
 RGB{N0f8}(0.071,0.071,0.071)     RGB{N0f8}(0.773,0.773,0.773)
 RGB{N0f8}(0.765,0.765,0.765)     RGB{N0f8}(0.322,0.322,0.322)
 RGB{N0f8}(0.902,0.902,0.902)     RGB{N0f8}(0.576,0.576,0.576)```
"""
substitution_encryption!(
    image::Array{RGB{N0f8},2},
    keys::Vector{Int64};
    path_for_result::String="./encrypted.png"
) = _substitution(image, keys, :encrypt; path_for_result=path_for_result, inplace=true)

"""
    substitution_decryption(image, keys; path_for_result="./decrypted.png")

Performs substitution decryption on a given image with the given keys.

# Algorithm
Iterates simulataneously over each pixel and key, and XORs the pixel value
(all R, G, and B) with the given key. Hence, the keys provided must be the same
as the ones provided during encryption.

# Arguments
- `image::Union{String,Array{RGB{N0f8},2}}`: The path to the image or the loaded image to be decrypted.
- `keys::Array{Int64, 1}`: Keys for decryption.
- `path_for_result::String`: The path for storing the decrypted image.

# Returns
- `image::Array{RGB{N0f8}, 2}`: Decrypted image.

# Example
```jldoctest
julia> using Images

julia> img = load("../test_images/camera.jfif")
225×225 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 ⋮                             ⋱
 RGB{N0f8}(0.4,0.4,0.4)           RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.447,0.447,0.447)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.427,0.427,0.427)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.451,0.451,0.451)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.51,0.51,0.51)     …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.537,0.537,0.537)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.412,0.412,0.412)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.149,0.149,0.149)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.0,0.0,0.0)           RGB{N0f8}(0.0,0.0,0.0)

julia> height, width = size(img)
(225, 225)

julia> keys = logistic_key(0.01, 3.97, height * width)
50625-element Vector{Int64}:
   0
  44
   7
  26
  14
 224
  16
 250
 162
 211
   ⋮
  72
 255
  80
  84
  21
 154
 197
  82
 147

julia> substitution_decryption(img, keys)
[ Info: DECRYPTING
[ Info: DECRYPTED
225×225 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.839,0.839,0.839)
 RGB{N0f8}(0.718,0.718,0.718)     RGB{N0f8}(0.361,0.361,0.361)
 RGB{N0f8}(0.612,0.612,0.612)     RGB{N0f8}(0.055,0.055,0.055)
 RGB{N0f8}(0.506,0.506,0.506)     RGB{N0f8}(0.58,0.58,0.58)
 RGB{N0f8}(0.584,0.584,0.584)     RGB{N0f8}(0.631,0.631,0.631)
 RGB{N0f8}(0.482,0.482,0.482)  …  RGB{N0f8}(0.639,0.639,0.639)
 RGB{N0f8}(0.545,0.545,0.545)     RGB{N0f8}(0.196,0.196,0.196)
 RGB{N0f8}(0.38,0.38,0.38)        RGB{N0f8}(0.047,0.047,0.047)
 RGB{N0f8}(0.224,0.224,0.224)     RGB{N0f8}(0.847,0.847,0.847)
 RGB{N0f8}(0.282,0.282,0.282)     RGB{N0f8}(0.827,0.827,0.827)
 ⋮                             ⋱
 RGB{N0f8}(0.573,0.573,0.573)     RGB{N0f8}(0.282,0.282,0.282)
 RGB{N0f8}(0.361,0.361,0.361)     RGB{N0f8}(1.0,1.0,1.0)
 RGB{N0f8}(0.847,0.847,0.847)     RGB{N0f8}(0.314,0.314,0.314)
 RGB{N0f8}(0.349,0.349,0.349)     RGB{N0f8}(0.329,0.329,0.329)
 RGB{N0f8}(0.984,0.984,0.984)  …  RGB{N0f8}(0.082,0.082,0.082)
 RGB{N0f8}(0.031,0.031,0.031)     RGB{N0f8}(0.604,0.604,0.604)
 RGB{N0f8}(0.071,0.071,0.071)     RGB{N0f8}(0.773,0.773,0.773)
 RGB{N0f8}(0.765,0.765,0.765)     RGB{N0f8}(0.322,0.322,0.322)
 RGB{N0f8}(0.902,0.902,0.902)     RGB{N0f8}(0.576,0.576,0.576)
```
"""
substitution_decryption(
    image::Union{String,Array{RGB{N0f8},2}},
    keys::Vector{Int64};
    path_for_result::String="./decrypted.png",
) = _substitution(image, keys, :decrypt; path_for_result=path_for_result)

"""
    substitution_decryption!(image, keys; path_for_result="./decrypted.png")

Performs substitution decryption on a given image with the given keys.

# Algorithm
Iterates simulataneously over each pixel and key, and XORs the pixel value
(all R, G, and B) with the given key. Hence, the keys provided must be the same
as the ones provided during encryption.

# Arguments
- `image::Union{String,Array{RGB{N0f8},2}}`: The path to the image or the loaded image to be decrypted.
- `keys::Array{Int64, 1}`: Keys for decryption.
- `path_for_result::String`: The path for storing the decrypted image.

# Returns
- `image::Array{RGB{N0f8}, 2}`: Decrypted image.

# Example
```jldoctest
julia> using Images

julia> img = load("../test_images/camera.jfif")
225×225 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.608,0.608,0.608)     RGB{N0f8}(0.0,0.0,0.0)
 ⋮                             ⋱
 RGB{N0f8}(0.4,0.4,0.4)           RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.447,0.447,0.447)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.427,0.427,0.427)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.451,0.451,0.451)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.51,0.51,0.51)     …  RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.537,0.537,0.537)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.412,0.412,0.412)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.149,0.149,0.149)     RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.0,0.0,0.0)           RGB{N0f8}(0.0,0.0,0.0)

julia> height, width = size(img)
(225, 225)

julia> keys = logistic_key(0.01, 3.97, height * width)
50625-element Vector{Int64}:
   0
  44
   7
  26
  14
 224
  16
 250
 162
 211
   ⋮
  72
 255
  80
  84
  21
 154
 197
  82
 147

julia> substitution_decryption!(img, keys)
[ Info: DECRYPTING
[ Info: DECRYPTED
225×225 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.839,0.839,0.839)
 RGB{N0f8}(0.718,0.718,0.718)     RGB{N0f8}(0.361,0.361,0.361)
 RGB{N0f8}(0.612,0.612,0.612)     RGB{N0f8}(0.055,0.055,0.055)
 RGB{N0f8}(0.506,0.506,0.506)     RGB{N0f8}(0.58,0.58,0.58)
 RGB{N0f8}(0.584,0.584,0.584)     RGB{N0f8}(0.631,0.631,0.631)
 RGB{N0f8}(0.482,0.482,0.482)  …  RGB{N0f8}(0.639,0.639,0.639)
 RGB{N0f8}(0.545,0.545,0.545)     RGB{N0f8}(0.196,0.196,0.196)
 RGB{N0f8}(0.38,0.38,0.38)        RGB{N0f8}(0.047,0.047,0.047)
 RGB{N0f8}(0.224,0.224,0.224)     RGB{N0f8}(0.847,0.847,0.847)
 RGB{N0f8}(0.282,0.282,0.282)     RGB{N0f8}(0.827,0.827,0.827)
 ⋮                             ⋱
 RGB{N0f8}(0.573,0.573,0.573)     RGB{N0f8}(0.282,0.282,0.282)
 RGB{N0f8}(0.361,0.361,0.361)     RGB{N0f8}(1.0,1.0,1.0)
 RGB{N0f8}(0.847,0.847,0.847)     RGB{N0f8}(0.314,0.314,0.314)
 RGB{N0f8}(0.349,0.349,0.349)     RGB{N0f8}(0.329,0.329,0.329)
 RGB{N0f8}(0.984,0.984,0.984)  …  RGB{N0f8}(0.082,0.082,0.082)
 RGB{N0f8}(0.031,0.031,0.031)     RGB{N0f8}(0.604,0.604,0.604)
 RGB{N0f8}(0.071,0.071,0.071)     RGB{N0f8}(0.773,0.773,0.773)
 RGB{N0f8}(0.765,0.765,0.765)     RGB{N0f8}(0.322,0.322,0.322)
 RGB{N0f8}(0.902,0.902,0.902)     RGB{N0f8}(0.576,0.576,0.576)

julia> img  # inplace
225×225 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(0.608,0.608,0.608)  …  RGB{N0f8}(0.839,0.839,0.839)
 RGB{N0f8}(0.718,0.718,0.718)     RGB{N0f8}(0.361,0.361,0.361)
 RGB{N0f8}(0.612,0.612,0.612)     RGB{N0f8}(0.055,0.055,0.055)
 RGB{N0f8}(0.506,0.506,0.506)     RGB{N0f8}(0.58,0.58,0.58)
 RGB{N0f8}(0.584,0.584,0.584)     RGB{N0f8}(0.631,0.631,0.631)
 RGB{N0f8}(0.482,0.482,0.482)  …  RGB{N0f8}(0.639,0.639,0.639)
 RGB{N0f8}(0.545,0.545,0.545)     RGB{N0f8}(0.196,0.196,0.196)
 RGB{N0f8}(0.38,0.38,0.38)        RGB{N0f8}(0.047,0.047,0.047)
 RGB{N0f8}(0.224,0.224,0.224)     RGB{N0f8}(0.847,0.847,0.847)
 RGB{N0f8}(0.282,0.282,0.282)     RGB{N0f8}(0.827,0.827,0.827)
 ⋮                             ⋱
 RGB{N0f8}(0.573,0.573,0.573)     RGB{N0f8}(0.282,0.282,0.282)
 RGB{N0f8}(0.361,0.361,0.361)     RGB{N0f8}(1.0,1.0,1.0)
 RGB{N0f8}(0.847,0.847,0.847)     RGB{N0f8}(0.314,0.314,0.314)
 RGB{N0f8}(0.349,0.349,0.349)     RGB{N0f8}(0.329,0.329,0.329)
 RGB{N0f8}(0.984,0.984,0.984)  …  RGB{N0f8}(0.082,0.082,0.082)
 RGB{N0f8}(0.031,0.031,0.031)     RGB{N0f8}(0.604,0.604,0.604)
 RGB{N0f8}(0.071,0.071,0.071)     RGB{N0f8}(0.773,0.773,0.773)
 RGB{N0f8}(0.765,0.765,0.765)     RGB{N0f8}(0.322,0.322,0.322)
 RGB{N0f8}(0.902,0.902,0.902)     RGB{N0f8}(0.576,0.576,0.576)
```
"""
substitution_decryption!(
    image::Array{RGB{N0f8},2},
    keys::Vector{Int64};
    path_for_result::String="./decrypted.png",
) = _substitution(image, keys, :decrypt; path_for_result=path_for_result, inplace=true)


text:
## Solve the first linear system
x,y,z = [0 1/2 1/4; 1 1 2; 2 1 1.5] \ [ 2,6,7]
## Solve the reduced carbs linear system
x,y,z = [0 1/2 1/4; 1 1 2; 2 1 1.5] \ [ 2.5,6.5,5]
## Solve the LP
using JuMP
using Clp
##
model = Model(with_optimizer(Clp.Optimizer))
@variable(model, 0 <= x)
@variable(model, 0 <= y)
@variable(model, 0 <= z)
@objective(model, Min, x+y+z)
@constraint(model,
  [0 1/2 1/4; 1 1 2; 1 1 1.5]*[x,y,z] .>= [2.5,6.5,5])
optimize!(model)
##
value(x), value(y), value(z)


text:
# Google cloud storage resource

import .PoorGCloudAuth: GCSClient
import .PoorGCloudAuth

struct GCSResource <: RemoteResource
    gsurl::String
    client::GCSClient
    size::Int64
    function GCSResource(gsurl, client)
        size = PoorGCloudAuth.fetch_object_size(gsurl, client)
        resource = new(gsurl, client, size)
        resource
    end
end

const GCSFileIO = RemoteResourceIO{GCSResource}

function RemoteResourceIO{GCSResource}(gsurl::AbstractString; kwargs...)
    RemoteResourceIO(GCSResource(gsurl); kwargs...)
end

function RemoteResourceIO{GCSResource}(gsurl::AbstractString, client::GCSClient; kwargs...)
    RemoteResourceIO(GCSResource(gsurl, client); kwargs...)
end

function RemoteResourceIO{GCSResource}(client::GCSClient; kwargs...)
    gsurl -> GCSResource(gsurl, client; kwargs...)
end

function GCSResource(gsurl::AbstractString)
    GCSResource(gsurl, GCSClient())
end

function resc_length(x::GCSResource)
    x.size
end

function resc_fetch(x::GCSResource, range::UnitRange{<:Integer})
    PoorGCloudAuth.fetch_gcs(x.gsurl, x.client; range = range)
end


text:
function _create_juliadefault()
    cs = ColorScheme()
    symbol!(cs, Crayon(bold = true))
    comment!(cs, Crayon(bold = true))
    string!(cs, Crayon(bold = true))
    call!(cs, Crayon(bold = true))
    op!(cs, Crayon(bold = true))
    keyword!(cs, Crayon(bold = true))
    text!(cs, Crayon(bold = true))
    macro!(cs, Crayon(bold = true))
    function_def!(cs, Crayon(bold = true))
    error!(cs, Crayon(bold = true))
    argdef!(cs, Crayon(bold = true))
    number!(cs, Crayon(bold = true))
    return cs
end

# Try to represent the Monokai colorscheme.
function _create_monokai()
    cs = ColorScheme()
    symbol!(cs, Crayon(foreground = :magenta))
    comment!(cs, Crayon(foreground = :dark_gray))
    string!(cs, Crayon(foreground = :yellow))
    call!(cs, Crayon(foreground = :cyan))
    op!(cs, Crayon(foreground = :light_red))
    keyword!(cs, Crayon(foreground = :light_red))
    text!(cs, Crayon(foreground = :default))
    macro!(cs, Crayon(foreground = :cyan))
    function_def!(cs, Crayon(foreground = :green))
    error!(cs, Crayon(foreground = :default))
    argdef!(cs, Crayon(foreground = :cyan))
    number!(cs, Crayon(foreground = :magenta))
    return cs
end

function _create_monokai_256()
    cs = ColorScheme()
    symbol!(cs, Crayon(foreground = 141)) # purpleish
    comment!(cs, Crayon(foreground = 60)) # greyish
    string!(cs, Crayon(foreground = 208)) # beigish
    call!(cs, Crayon(foreground = 81)) # cyanish
    op!(cs, Crayon(foreground = 197)) # redish
    keyword!(cs, Crayon(foreground = 197)) # redish
    text!(cs, Crayon(foreground = :default))
    macro!(cs, Crayon(foreground = 208)) # cyanish
    function_def!(cs, Crayon(foreground = 148))
    error!(cs, Crayon(foreground = :default))
    argdef!(cs, Crayon(foreground = 81))  # cyanish
    number!(cs, Crayon(foreground = 141)) # purpleish
    return cs
end

function _create_monokai_24()
    cs = ColorScheme()
    symbol!(cs, Crayon(foreground = (174,129,255)))
    comment!(cs, Crayon(foreground = (89, 89, 89)))
    string!(cs, Crayon(foreground = (253,151,31)))
    call!(cs, Crayon(foreground = (102,217,239)))
    op!(cs, Crayon(foreground = (249, 38, 114)))
    keyword!(cs, Crayon(foreground = (249, 38, 114)))
    text!(cs, Crayon(foreground = :default))
    macro!(cs, Crayon(foreground = (253,151,31)))
    function_def!(cs, Crayon(foreground = (166,226,42)))
    error!(cs, Crayon(foreground = :default))
    argdef!(cs, Crayon(foreground = (102,217,239)))
    number!(cs, Crayon(foreground = (174,129,255)))
    return cs
end

function _create_boxymonokai_256()
    cs = ColorScheme()
    symbol!(cs, Crayon(foreground = 148))
    comment!(cs, Crayon(foreground = 95))
    string!(cs, Crayon(foreground = 148))
    call!(cs, Crayon(foreground = 81))
    op!(cs, Crayon(foreground = 158))
    keyword!(cs, Crayon(foreground = 141))
    text!(cs, Crayon(foreground = :default))
    macro!(cs, Crayon(foreground = 81))
    function_def!(cs, Crayon(foreground = 81))
    error!(cs, Crayon(foreground = :default))
    argdef!(cs, Crayon(foreground = 186))
    number!(cs, Crayon(foreground = 208))
    return cs
end

function _create_tomorrow_night_bright()
    cs = ColorScheme()
    symbol!(cs, Crayon(foreground = 185))
    comment!(cs, Crayon(foreground = 246))
    string!(cs, Crayon(foreground = 185))
    call!(cs, Crayon(foreground = 73))
    op!(cs, Crayon(foreground = 231))
    keyword!(cs, Crayon(foreground = 140))
    text!(cs, Crayon(foreground = :default))
    macro!(cs, Crayon(foreground = 73))
    function_def!(cs, Crayon(foreground = 110))
    error!(cs, Crayon(foreground = :default))
    argdef!(cs, Crayon(foreground = 255)) # nothing special added here
    number!(cs, Crayon(foreground = 208))
    return cs
end

function _create_tomorrow_night_bright_24()
    cs = ColorScheme()
    symbol!(cs, Crayon(foreground = (185, 202, 74)))
    comment!(cs, Crayon(foreground = (150, 152, 150)))
    string!(cs, Crayon(foreground = (185, 202, 74)))
    call!(cs, Crayon(foreground = (112, 192, 177)))
    op!(cs, Crayon(foreground = (255, 255, 255)))
    keyword!(cs, Crayon(foreground = (195, 151, 216)))
    text!(cs, Crayon(foreground = :default))
    macro!(cs, Crayon(foreground = (112, 192, 177)))
    function_def!(cs, Crayon(foreground = (122, 166, 218)))
    error!(cs, Crayon(foreground = :default))
    argdef!(cs, Crayon(foreground = 255)) # nothing special added here
    number!(cs, Crayon(foreground = (231, 140, 69)))
    return cs
end

function _create_tomorrow_24()
    cs = ColorScheme()
    symbol!(cs, Crayon(foreground = (113, 140, 0)))
    comment!(cs, Crayon(foreground = (142, 144, 140)))
    string!(cs, Crayon(foreground = (113, 140, 0)))
    call!(cs, Crayon(foreground = (62, 153, 159)))
    op!(cs, Crayon(foreground = :default))
    keyword!(cs, Crayon(foreground = (137, 89, 168)))
    text!(cs, Crayon(foreground = :default))
    macro!(cs, Crayon(foreground = (62, 153, 159)))
    function_def!(cs, Crayon(foreground = (66, 113, 174)))
    error!(cs, Crayon(foreground = :default))
    argdef!(cs, Crayon(foreground = :default)) # nothing special added here
    number!(cs, Crayon(foreground = (245, 135, 31)))
    return cs
end

function _create_tomorrow_256()
    cs = ColorScheme()
    symbol!(cs, Crayon(foreground = 64))
    comment!(cs, Crayon(foreground = 102))
    string!(cs, Crayon(foreground = 64))
    call!(cs, Crayon(foreground = 67))
    op!(cs, Crayon(foreground = :default))
    keyword!(cs, Crayon(foreground = 97))
    text!(cs, Crayon(foreground = :default))
    macro!(cs, Crayon(foreground = 67))
    function_def!(cs, Crayon(foreground = 61))
    error!(cs, Crayon(foreground = :default))
    argdef!(cs, Crayon(foreground = :default)) # nothing special added here
    number!(cs, Crayon(foreground = 208))
    return cs
end

function _create_distinguished()
    cs = ColorScheme()
    symbol!(cs, Crayon(foreground = 66, bold=true))
    comment!(cs, Crayon(foreground = 243))
    string!(cs, Crayon(foreground = 143))
    call!(cs, Crayon(foreground = :default))
    op!(cs, Crayon(foreground = 180))
    keyword!(cs, Crayon(foreground = 173))
    text!(cs, Crayon(foreground = :default))
    macro!(cs, Crayon(foreground = 247))
    function_def!(cs, Crayon(foreground = :default))
    error!(cs, Crayon(foreground = :default))
    argdef!(cs, Crayon(foreground = 67))
    number!(cs, Crayon(foreground = 173))
    return cs
end


text:
# syntax: proto2
mutable struct Version <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Version(; kwargs...)
        obj = new(meta(Version), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Version
const __meta_Version = Ref{ProtoMeta}()
function meta(::Type{Version})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Version)
            __meta_Version[] = target = ProtoMeta(Version)
            allflds = Pair{Symbol,Union{Type,String}}[:major => Int32, :minor => Int32, :patch => Int32, :suffix => AbstractString]
            meta(target, Version, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Version[]
    end
end
function Base.getproperty(obj::Version, name::Symbol)
    if name === :major
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :minor
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :patch
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :suffix
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct CodeGeneratorRequest <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CodeGeneratorRequest(; kwargs...)
        obj = new(meta(CodeGeneratorRequest), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CodeGeneratorRequest
const __meta_CodeGeneratorRequest = Ref{ProtoMeta}()
function meta(::Type{CodeGeneratorRequest})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CodeGeneratorRequest)
            __meta_CodeGeneratorRequest[] = target = ProtoMeta(CodeGeneratorRequest)
            fnum = Int[1,2,15,3]
            allflds = Pair{Symbol,Union{Type,String}}[:file_to_generate => Base.Vector{AbstractString}, :parameter => AbstractString, :proto_file => Base.Vector{ProtoBuf.GoogleProtoBuf.FileDescriptorProto}, :compiler_version => Version]
            meta(target, CodeGeneratorRequest, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CodeGeneratorRequest[]
    end
end
function Base.getproperty(obj::CodeGeneratorRequest, name::Symbol)
    if name === :file_to_generate
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    elseif name === :parameter
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :proto_file
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ProtoBuf.GoogleProtoBuf.FileDescriptorProto}
    elseif name === :compiler_version
        return (obj.__protobuf_jl_internal_values[name])::Version
    else
        getfield(obj, name)
    end
end

mutable struct CodeGeneratorResponse_File <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CodeGeneratorResponse_File(; kwargs...)
        obj = new(meta(CodeGeneratorResponse_File), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CodeGeneratorResponse_File
const __meta_CodeGeneratorResponse_File = Ref{ProtoMeta}()
function meta(::Type{CodeGeneratorResponse_File})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CodeGeneratorResponse_File)
            __meta_CodeGeneratorResponse_File[] = target = ProtoMeta(CodeGeneratorResponse_File)
            fnum = Int[1,2,15]
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :insertion_point => AbstractString, :content => AbstractString]
            meta(target, CodeGeneratorResponse_File, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CodeGeneratorResponse_File[]
    end
end
function Base.getproperty(obj::CodeGeneratorResponse_File, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :insertion_point
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :content
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct CodeGeneratorResponse <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CodeGeneratorResponse(; kwargs...)
        obj = new(meta(CodeGeneratorResponse), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CodeGeneratorResponse
const __meta_CodeGeneratorResponse = Ref{ProtoMeta}()
function meta(::Type{CodeGeneratorResponse})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CodeGeneratorResponse)
            __meta_CodeGeneratorResponse[] = target = ProtoMeta(CodeGeneratorResponse)
            fnum = Int[1,15]
            allflds = Pair{Symbol,Union{Type,String}}[:error => AbstractString, :file => Base.Vector{CodeGeneratorResponse_File}]
            meta(target, CodeGeneratorResponse, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CodeGeneratorResponse[]
    end
end
function Base.getproperty(obj::CodeGeneratorResponse, name::Symbol)
    if name === :error
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :file
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{CodeGeneratorResponse_File}
    else
        getfield(obj, name)
    end
end

export Version, CodeGeneratorRequest, CodeGeneratorResponse_File, CodeGeneratorResponse


text:
using NMF
using Pkg
pkg"activate ."
using NNSparseCoding; const N = NNSparseCoding;
using Test
using Plots;
theme(:dark);
using BenchmarkTools

function scatty(M1,M2)
  scatter(M1[:], M2[:]; leg = false, ratio = 1)
  return plot!([extrema(M1)...].* 1.1, identity; linewidth = 3)
end

##

n = 80
d = n
k = 300
A,S = N.initialize_mats(n,d,k)
# let's have a sparse S
map!(s-> rand() < 0.9 ? 0.0 : 3s ,S,S)
X = A*S
lambda=0.1
nsteps=10_000
stepsize=0.2
optthresh=1E-4
o = N.nnsc_optimize(X,d,lambda ;
  step_zero=stepsize, nstepmax=nsteps,
  opt_threshold=optthresh,verbose=true)

extrema(o.AS .- X)
##
# X2,res = let res = nnmf(X,d)
#   (res.W * res.H,res)
# end

(W,H) = let (W,H) = NMF.nndsvd(X,d ;variant=:a),
  alg = NMF.CoordinateDescent{Float64}(
      l₁ratio=0.8,α=10.0,regularization=:both)
  res = NMF.solve!(alg,X,W,H)
  (res.W,res.H)
end
X2=W*H

scatty(X2,X)
scatty(o.AS,X)

histogram(o.S[:])
histogram(H[:])

heatmap(o.S)
heatmap(H)


text:
# ---
# title: 815. Bus Routes
# id: problem815
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/bus-routes/description/>
# hidden: true
# ---
#
# We have a list of bus routes. Each `routes[i]` is a bus route that the i-th
# bus repeats forever. For example if `routes[0] = [1, 5, 7]`, this means that
# the first bus (0-th indexed) travels in the sequence 1->5->7->1->5->7->1->...
# forever.
#
# We start at bus stop `S` (initially not on a bus), and we want to go to bus
# stop `T`. Travelling by buses only, what is the least number of buses we must
# take to reach our destination? Return -1 if it is not possible.
#
#
#
#     **Example:**
#     Input:
#     routes = [[1, 2, 7], [3, 6, 7]]
#     S = 1
#     T = 6
#     Output: 2
#     Explanation:
#     The best strategy is take the first bus to the bus stop 7, then take the second bus to the bus stop 6.
#
#
#
#
# **Constraints:**
#
#   * `1 <= routes.length <= 500`.
#   * `1 <= routes[i].length <= 10^5`.
#   * `0 <= routes[i][j] < 10 ^ 6`.
#
#
## @lc code=start
using LeetCode

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


text:
module Chips

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export chip

register_normal_element("q__chip", context = @__MODULE__)

function chip(args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)

  wrap() do
    q__chip(args...; attributes([kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

end


text:
"""
This function defines the constraints for the water level (or state of charge)
for the Hydro Reservoir.
"""
function energy_target_constraint!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, S},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.Component, S <: AbstractDeviceFormulation}
    key = ICKey(EnergyLevel, T)
    parameters = model_has_parameters(optimization_container)
    use_forecast_data = model_uses_forecasts(optimization_container)
    time_steps = model_time_steps(optimization_container)
    target_forecast_label = "storage_target"
    constraint_infos_target = Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    if use_forecast_data
        for (ix, d) in enumerate(devices)
            ts_vector_target =
                get_time_series(optimization_container, d, target_forecast_label)
            constraint_info_target = DeviceTimeSeriesConstraintInfo(
                d,
                x -> get_target_multiplier(x),
                ts_vector_target,
            )
            constraint_infos_target[ix] = constraint_info_target
        end
    else
        for (ix, d) in enumerate(devices)
            ts_vector_target =
                length(time_steps) == 1 ? [PSY.get_storage_target(d)] :
                vcat(zeros(time_steps[end - 1]), PSY.get_storage_target(d))
            constraint_info_target = DeviceTimeSeriesConstraintInfo(
                d,
                x -> get_target_multiplier(x),
                ts_vector_target,
            )
            constraint_infos_target[ix] = constraint_info_target
        end
    end

    if parameters
        energy_target_param!(
            optimization_container,
            constraint_infos_target,
            make_constraint_name(ENERGY_TARGET, T),
            (
                make_variable_name(ENERGY, T),
                make_variable_name(ENERGY_SHORTAGE, T),
                make_variable_name(ENERGY_SURPLUS, T),
            ),
            UpdateRef{T}(TARGET, target_forecast_label),
        )
    else
        energy_target!(
            optimization_container,
            constraint_infos_target,
            make_constraint_name(ENERGY_TARGET, T),
            (
                make_variable_name(ENERGY, T),
                make_variable_name(ENERGY_SHORTAGE, T),
                make_variable_name(ENERGY_SURPLUS, T),
            ),
        )
    end
    return
end


text:
### Define the target distribution and its gradient
using DiffResults: GradientResult
import DiffResults
using ForwardDiff: gradient!

### Build up a HMC sampler to draw samples
using AdvancedHMC

export advancedHMC

"""
    advancedHMC(m, data, N = 1000; n_adapts = 1000)

Draw `N` samples from the posterior distribution of parameters defined in Soss model `m`, conditional on `data`. Samples are drawn using Hamiltonial Monte Carlo (HMC) from the `advancedHMC.jl` package.

## Keywords
*  `n_adapts = 1000`: The number of interations used to set HMC parameters.

Returns a tuple of length 2:
1. Samples from the posterior distribution of parameters.
2. Sample summary statistics.

## Example

```
m = @model x begin
    β ~ Normal()
    yhat = β .* x
    y ~ For(eachindex(x)) do j
        Normal(yhat[j], 2.0)
    end
end

x = randn(10);
truth = [0.61, -0.34, -1.74];

post = advancedHMC(m(x=x), (y=truth,));
E_β = mean(post[1])[1]

println("Posterior mean β: " * string(round(E_β, digits=2)))
```

"""
function advancedHMC(m :: ConditionalModel{A,B}, _data, N = 1000;
                                                         n_adapts  = 1000) where {A,B}

    ℓ(pars) = logdensity_def(m, merge(pars, _data))

    t = as(m,_data)

    function f(x)
        (θ, logjac) = transform_and_logjac(t,x)
        ℓ(θ) + logjac
    end

    function ∂f(x)
        res = GradientResult(x)
        gradient!(res, f, x)
        return (DiffResults.value(res), DiffResults.gradient(res))
    end

    # Sampling parameter settings

    D = t.dimension
    # Draw a random starting points
    initial_θ = randn(D)

    # Define metric space, Hamiltonian, sampling method and adaptor
    metric = DiagEuclideanMetric(D)
    hamiltonian = Hamiltonian(metric, f, ∂f)
    prop = AdvancedHMC.NUTS(Leapfrog(find_good_stepsize(hamiltonian, initial_θ)))
    initial_ϵ = find_good_stepsize(hamiltonian, initial_θ)
    integrator = Leapfrog(initial_ϵ)
    adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))

    # Draw samples via simulating Hamiltonian dynamics
    # - `samples` will store the samples
    # - `stats` will store statistics for each sample
    samples, stats = simulate(hamiltonian, prop, initial_θ, N, adaptor, n_adapts; progress=false, verbose=false)

end


text:
using Try

import ..TryExperimental: @tryable, istryable
using ..TryExperimental: TryExperimental
include("tryable.jl")

macro exported_function(name::Symbol)
    m = @__MODULE__
    quote
        $m.@tryable($name)
        export $name
    end |> esc
end