text:
export armijo, wolfe, armijo_wolfe, shamanskii_stop, goldstein

"""
armijo: check if a step size is admissible according to the Armijo criterion.

Armijo criterion: f(x + θd) - f(x) < τ₀ θ ∇f(x+θd)d

Note: ht, h₀ and g₀ are required in the LSAtT
"""
function armijo(h      :: Any, #LineModel
                h_at_t :: LSAtT;
                τ₀     :: Float64 = 0.01,
                kwargs...)

    if (h_at_t.ht == nothing) || (h_at_t.h₀ == nothing) || (h_at_t.g₀ == nothing)
     #printstyled("Warning: Nothing entries in the State.\n", color = :red)
     return throw(error("Nothing entries in the State. ht, h₀ and g₀ are mandatory."))
    else

     fact = -0.8
     Eps = 1e-10
     hgoal = h_at_t.ht - h_at_t.h₀ - h_at_t.g₀ * h_at_t.x * τ₀
     # Armijo = (h_at_t.ht <= hgoal)# || ((h_at_t.ht <= h_at_t.h₀ + Eps * abs(h_at_t.h₀)) & (h_at_t.gt <= fact * h_at_t.g₀))
	 # Armijo_HZ =
	 # positive = h_at_t.x > 0.0   # positive step
     return max(hgoal, 0.0)
    end
end

"""
wolfe: check if a step size is admissible according to the Wolfe criterion.

Strong Wolfe criterion: |∇f(x+θd)| < τ₁||∇f(x)||.

Note: gt and g₀ are required in the LSAtT
"""
function wolfe(h      :: Any, #LineModel,
               h_at_t :: LSAtT;
               τ₁     :: Float64 = 0.99,
               kwargs...)

    if (h_at_t.g₀ == nothing) || (h_at_t.gt == nothing)
     return throw(error("Nothing entries in the State."))
    else

     #wolfe = (τ₁ .* h_at_t.g₀) - (abs(h_at_t.gt))
     wolfe = abs(h_at_t.gt) - τ₁ * abs(h_at_t.g₀)
	 #positive = h_at_t.x > 0.0   # positive step
     return max(wolfe, 0.0)
    end
end

"""
armijo_wolfe: check if a step size is admissible according to the Armijo and Wolfe criteria.

Note: ht, h₀, gt and g₀ are required in the LSAtT
"""
function armijo_wolfe(h      :: Any, #LineModel,
                      h_at_t :: LSAtT;
                      τ₀     :: Float64 = 0.01,
                      τ₁     :: Float64 = 0.99,
                      kwargs...)

   if (h_at_t.ht == nothing) || (h_at_t.gt == nothing) || (h_at_t.h₀ == nothing) || (h_at_t.g₀ == nothing)
    #printstyled("Warning: Nothing entries in the State.\n", color = :red)
    return throw(error("Nothing entries in the State. ht, h₀, gt and g₀ are mandatory."))
   else
    wolfe  = abs(h_at_t.gt) - τ₁ * abs(h_at_t.g₀)
    armijo = h_at_t.ht - h_at_t.h₀ - h_at_t.g₀ * h_at_t.x * τ₀
    return max(armijo, wolfe, 0.0)
   end
end

"""
shamanskii_stop: check if a step size is admissible according to the "Shamanskii" criteria.
This criteria was proposed in:
Lampariello, F., & Sciandrone, M. (2001).
Global convergence technique for the Newton method with periodic Hessian evaluation.
Journal of optimization theory and applications, 111(2), 341-358.

Note: * h.d accessible (specific LineModel)
      * ht, h₀ are required in the LSAtT
"""
function shamanskii_stop(h      :: Any, #LineModel,
                         h_at_t :: LSAtT;
                         γ      :: Float64 = 1.0e-09,
                         kwargs...)

    admissible = h_at_t.ht - h_at_t.h₀ - γ * (h_at_t.x)^3 * norm(h.d)^3
    return max(admissible, 0.0)
end

"""
goldstein: check if a step size is admissible according to the Goldstein criteria.

Note: ht, h₀ and g₀ are required in the LSAtT
"""
function goldstein(h      :: Any, #LineModel,
                   h_at_t :: LSAtT;
                   τ₀     :: Float64 = 0.0001,
                   τ₁     :: Float64 = 0.9999,
                   kwargs...)

 	goldstein = max(h_at_t.h₀ + h_at_t.x * (1 - τ₀) * h_at_t.g₀ - h_at_t.ht,
                    h_at_t.ht - (h_at_t.h₀ + h_at_t.x *  τ₀ * h_at_t.g₀))
 	# positive = h_at_t.x > 0.0   # positive step
 	return max(goldstein, 0.0) #&& positive
end


text:
# ------------------------------------------------------------
# Reconstruction of impulse responses, i.e., a single coefficient in the Haar
# basis from equidistant frequency samples

# ------------------------------------------------------------
# Reconstruction of scaling function

# Input:
# M: The number of frequency samples
# epsislon: The distance between two consecutive samples
# N: The number of wavelet coefficients recovered

M = 8
epsilon = 0.25
J = floor(Int, log2(M))

# Sampling points
xi = [-M:M;]*epsilon

# Change of basis matrix
T = freq2Haar(xi, J)

# Fourier observations
b = FourHaarScaling(xi, J, 3)

# Scaling function coefficients: 1 in the 4th entry, 0 elsewhere
y = pinv(T)*b

# ------------------------------------------------------------
# Reconstruction of wavelet function

J2 = [0:2;]

# Change of basis matrix
T2 = freq2Haar(xi, J2)

# Fourier observations
b2 = FourHaarWavelet(xi, 2, 1)
#b2 = FourHaarScaling(xi, 1, 0)

# Wavelet coefficients: 1 in the 5th entry, 0 elsewhere
y2 = pinv(T2)*b2

ext: jl


text:
module Vmrs
using CUDA
using DocStringExtensions
using Adapt

export AbstractVmr, Vmr, VmrGM, init_vmr, get_vmr

"""
    AbstractVmr{FT}
"""
abstract type AbstractVmr{FT <: AbstractFloat} end
"""
    VmrGM{FT,FTA1D,FTA2D} <: AbstractVmr{FT}

Volume mixing ratios for various gases in the atmosphere. This struct can be used
when only H₂O and O₃ concentrations vary spatially and a global mean is used
for all other gases.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct VmrGM{FT <: AbstractFloat, FTA1D <: AbstractArray{FT, 1}, FTA2D <: AbstractArray{FT, 2}} <: AbstractVmr{FT}
    "volume mixing ratio of H₂O"
    vmr_h2o::FTA2D
    "volume mixing ratio of Ozone"
    vmr_o3::FTA2D
    "volume mixing ratio of all other gases, which are independent of location and column"
    vmr::FTA1D
end
VmrGM(vmr_h2o, vmr_o3, vmr) = VmrGM{eltype(vmr_h2o), typeof(vmr), typeof(vmr_h2o)}(vmr_h2o, vmr_o3, vmr)
Adapt.@adapt_structure VmrGM

function VmrGM(
    ::Type{FT},
    ::Type{DA},
    vmr_h2o::FTA2D,
    vmr_o3::FTA2D,
    vmr::FTA1D,
) where {FT <: AbstractFloat, FTA1D <: AbstractArray{FT, 1}, FTA2D <: AbstractArray{FT, 2}, DA}
    return VmrGM{FT, typeof(vmr), typeof(vmr_h2o)}(vmr_h2o, vmr_o3, vmr)
end

"""
    Vmr{FT,FTA3D} <: AbstractVmr{FT}

Volume mixing ratios for various gases in the atmosphere. This struct can be used
when concentrations vary spatially for all gases.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct Vmr{FT <: AbstractFloat, FTA3D <: AbstractArray{FT, 3}} <: AbstractVmr{FT}
    "volume mixing ratio of all gases as a function of location and column"
    vmr::FTA3D
end
Vmr(vmr) = Vmr{eltype(vmr), typeof(vmr)}(vmr)
Adapt.@adapt_structure Vmr
"""
    get_vmr(
        vmr::VmrGM{FT},
        ig::Int,
        ilay::Int,
        icol::Int,
    ) where {FT<:AbstractFloat}

Obtain volume mixing ratio of gas `ig` for layer `ilay` of column `icol`.
"""
@inline function get_vmr(vmr::VmrGM{FT}, ig::Int, ilay::Int, icol::Int) where {FT <: AbstractFloat}
    if ig == 0
        return FT(1)
    elseif ig == 1 # h2o / h2o_foreign / h2o_self-continua
        return @inbounds vmr.vmr_h2o[ilay, icol]
    elseif ig == 3 # ozone
        return @inbounds vmr.vmr_o3[ilay, icol]
    else # other gases
        return @inbounds vmr.vmr[ig]
    end
end

"""
    get_vmr(
        vmr::Vmr{FT},
        ig::Int,
        ilay::Int,
        icol::Int,
    ) where {FT<:AbstractFloat}

Obtain volume mixing ratio of gas `ig` for layer `ilay` of column `icol`.
"""
@inline function get_vmr(vmr::Vmr{FT}, ig::Int, ilay::Int, icol::Int) where {FT <: AbstractFloat}
    if ig == 0
        return FT(1)
    else
        return @inbounds vmr.vmr[ilay, icol, ig]
    end
end

"""
    init_vmr(
        ngas,
        nlay,
        ncol,
        ::Type{FT},
        ::Type{DA};
        gm::Bool = true,
    ) where {FT<:AbstractFloat,DA}

Initialize the Vmr struct.
"""
function init_vmr(
    ngas::Int,
    nlay::Int,
    ncol::Int,
    ::Type{FT},
    ::Type{DA};
    gm::Bool = true,
) where {FT <: AbstractFloat, DA}
    if gm
        FTA1D = DA{FT, 1}
        FTA2D = DA{FT, 2}
        return VmrGM{FT, FTA1D, FTA2D}(FTA2D(zeros(nlay, ncol)), FTA2D(zeros(nlay, ncol)), FTA1D(zeros(ngas)))
    else
        FTA3D = DA{FT, 3}
        return Vmr{FT, FTA3D}(FTA3D(zeros(nlay, ncol, ngas)))
    end
end

end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
import Base: iszero, zero, *

# Define addition of zeros.
for T in [
    :FiniteZeroKernel,
    :LhsFiniteZeroCrossKernel,
    :RhsFiniteZeroCrossKernel,
    :FiniteZeroCrossKernel,
]
    @eval +(x::$T, x′::$T) = zero(x)
end


text:
/**
 * This class models a cloud storage API.
 */
class Main {

    static int xTESTOUTPUT = 7;
    static int yTESTOUTPUT = 7;
    static int zTESTOUTPUT = 7;

    public static void main(String[] args) {

	// the following is dangerous data that should be erased sometime.
	condition c = false;

	int dangerous = ({L /c H})args.length;

	// get a key pair
	int pub = (PUBKEY(L /c H){L})42;
	int priv = (PRIVKEY(L /c H){L /c H})42;

	// let's encrypt it, and put the result in the cloud.

	Cloud cld = new Cloud();
	xTESTOUTPUT = CryptoLibrary.encrypt(pub, dangerous); // Should be L

	cld.put(CryptoLibrary.encrypt(pub, dangerous));

	// set the condition.
	c = true;

	yTESTOUTPUT = cld.get(); // should be L, since the ciphertext is still low security.

	// at this point, priv should no longer be accessible at PRIVKEY.
	zTESTOUTPUT = CryptoLibrary.decrypt(priv, cld.get()); // the level should be H, since it can no longer be at L /c H, since c was set.
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


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

using Markdown
using InteractiveUtils

# ╔═╡ dfc4afc0-69ad-11eb-0ece-c7507e610108
using GMT

# ╔═╡ 832779c0-6a37-11eb-21ab-79e35f04da4e
md"### Plot Landsat 8 images

The data used in this example is the band 4 (red channel) of a Landsat 8 scene. Those are relatively big images (~116 MB) so we will download it first and take that into account when the results of the commands bellow do not show instantly."

# ╔═╡ c7c5e480-6a38-11eb-2273-7dc4ddd09418
I = gmtread("/vsicurl/http://landsat-pds.s3.amazonaws.com/c1/L8/037/034/LC08_L1TP_037034_20160712_20170221_01_T1/LC08_L1TP_037034_20160712_20170221_01_T1_B4.TIF");

# ╔═╡ d07b0c8e-6a38-11eb-3371-69a77e57dae1
md"and have a look at what we got"

# ╔═╡ e0e5ee10-6a38-11eb-0b67-fb9fbf7cbdb4
imshow(I)

# ╔═╡ e905b3f0-6a38-11eb-3354-072a5d9152e0
md"_Ah, nice ... but it's so dark that we can't really see much!_

Well that is true and has one explanation. Modern satellite data is acquired with sensors of 12 bits or more but the data is stored in variables with 16 bits. This means that many of those data bits will not be used. But on screens we are limited to 8 bits per channel so we must scale the full 16 bits range [0 65535] to the [0 255] interval and in this process we tend to have much more dark pixels.

To see this better, let's look at the image's histogram.
"

# ╔═╡ 32c60de2-6a3a-11eb-1cfd-1b74d0ad957b
histogram(I, auto=true, bin=20, show=true)

# ╔═╡ 3a707cb0-6a3a-11eb-3975-b1b091d05e2f
md"We have used here the option **auto**=*true* that will try to guess where the data in histogram plot starts and ~ ends. It did behave well and we will use those numbers to do an operation that is called *histogram stretch* that consists in picking only part of the histogram and stretch it to [0 255]. And while at it we van visually observe that the limit [6000 24000] seems slightly better than the automatic one. Note that in fact we have data to the 40000 DN (Digital Number) but they are very few and at the end we must choose a balance to show *almos all* DNs and not making the image too dark. Reducing the higher value to 23400 would have made the image sligthy lighter."

# ╔═╡ 32e2c830-6a3b-11eb-2da3-cd08989b96a5
imshow(I, stretch=[6000 25000])

# ╔═╡ Cell order:
# ╠═dfc4afc0-69ad-11eb-0ece-c7507e610108
# ╠═a02af2c0-69ad-11eb-0cc5-57d1f41ddac5
# ╟─832779c0-6a37-11eb-21ab-79e35f04da4e
# ╠═c7c5e480-6a38-11eb-2273-7dc4ddd09418
# ╟─d07b0c8e-6a38-11eb-3371-69a77e57dae1
# ╠═e0e5ee10-6a38-11eb-0b67-fb9fbf7cbdb4
# ╟─e905b3f0-6a38-11eb-3354-072a5d9152e0
# ╠═32c60de2-6a3a-11eb-1cfd-1b74d0ad957b
# ╟─3a707cb0-6a3a-11eb-3975-b1b091d05e2f
# ╠═32e2c830-6a3b-11eb-2da3-cd08989b96a5


text:
include("../src/qspice.jl")

using QSpice.Gates, QSpice.State, QSpice.BitOps

function pow2between(n1, n2)
    for i = n1:n2
        if (i & (i - 1)) == 0
            return i
        end
    end
    error("Cannot find appropriate power of 2 in the given interval")
end

# Two random integers have a ~60% chance of being coprimes, so just resamping
# until we find one is a viable algorithm
function randcoprime(n)
    cprm = rand(3:n - 1)
    while gcd(cprm, n) != 1
        cprm = rand(3:n - 1)
    end
    return cprm
end

function approxrational(n::Real, limit::Int)
    m = [1 0; 0 1]

    current = n
    cint = floor(Int, current)

    while m[2,1] * cint + m[2, 2] <= limit
        m = m * [cint 1; 1 0]
        if current == float(cint)
            break
        end
        current = 1.0 / (current - cint)
        cint = floor(Int, current)
    end

    e1 = n - m[1,1] / m[2,1]

    cint = floor(Int, (limit - m[2,2]) / m[2,1])
    n2 = m[1,1] * cint + m[1,2]
    d2 = m[2,1] * cint + m[2,2]
    e2 = n - n2 / d2

    if abs(e1) <= abs(e2)
        return m[1,1], m[2,1]
    else
        return n2, d2
    end
end

function factorize(n)
    if isprime(n)
        return [n]
    end

    if iseven(n)
        return [2, factorize(n / 2)...]
    end

    q = pow2between(n * n, 2 * n * n - 1)
    x = randcoprime(n)

    factor = 1
    while factor == 1
        m = findperiod(n, q, x)
        r = approxrational(m / q, q)[2]
        factor = max(gcd(x^div(r, 2) + 1, n) % n,
                     gcd(x^div(r, 2) - 1, n) % n)
    end

    return [factor, div(n, factor)]
end

function findperiod(n, q, x)
    # First initialize register 1 with an equally weighed superposition of integers
    # from 0 to q - 1
    reg1 = fromstates([QUBIT0 for _ in 1:floor(Int, log2(q))]...)
    for i = 1:reg1.bits
        reg1 = hadamard(reg1, i)
    end

    # This part is not yet done with actual quantum gates/functions
    bits = ceil(Int, log2(n))
    reg2 = QuantumState([0.0 + 0im for _ in 1:2^bits], bits)
    for a = 0:q - 1
        reg2[powermod(x, a, n) + 1] += reg1[a + 1]
    end
    n = sqrt(sumabs2(reg2.vector))
    reg2.vector = reg2.vector .* (1.0 / n)

    reg = Gates.join(reg1, reg2)
    post, k = partialmeasure(reg, reg1.bits + 1:reg.bits...)

    collapsed, m = measure(qft(post))
    return todecimal(m)
end

max_stars_repo_name: leios/quantum

max_issues_repo_name: leios/quantum

max_forks_repo_name: leios/quantum


text:
using Test, LinearAlgebra, MkCell

# testing data
include("primtive_cells.jl")
include("paper/lattice_data.jl")

# tests
include("paper/volume.jl")
include("paper/angles.jl")
include("paper/rmax_sphere.jl")
include("paper/dmin.jl")
include("optimization.jl")


text:
abstract Population
abstract PopulationWithFitness <: Population

type FloatVectorPopulation <: PopulationWithFitness
  # The population is a matrix of floats, each row being an individual.
  individuals::Array{Float64, 2}

  # The fitnesses is a matrix of floats, each row being the fitness for one
  # individual, and each column corresponding to one objective/goal.
  fitness::Array{Float64, 2}
  fitness_scheme::FitnessScheme

  # A population always saves the top/best individuals seen during a search,
  # together with their fitness values.
  top::Array{Float64, 2}
  top_fitness::Array{Float64, 2}
  top_size::Int

  function FloatVectorPopulation(size = 100, dimensions = 1,
    scheme = float_vector_scheme_min(), numTopIndividuals = 10)
    inds = rand(size, dimensions)
    fs = ones(size, dimensions) * scheme.worst_fitness # Bug! Need not be the same num of objectives as there are dimensions...
    new(inds, fs, scheme, inds[1:numTopIndividuals,:], fs[1:numTopIndividuals,:], numTopIndividuals)
  end
end

best(pop::FloatVectorPopulation) = pop.top[1]
bestfitness(pop::FloatVectorPopulation) = pop.top_fitness[1]

# Update the toplist if this individual is good enough.
function update_toplist!(candidate, candidateFitness, population::PopulationWithFitness)
  fscheme = population.fitness_scheme
  if isbetter(candidateFitness, top_fitness[end], scheme)
    i = length(population.top)
    fs = population.top_fitness
    while(i >= 1 & isbetter(candidateFitness, fs[i], fscheme))
      i -= 1
    end
    population.top = vcat(population.top[1:i,:], candidate, population.top[i:end-1,:])
    population.top_fitness = vcat(population.top_fitness[1:i,:], candidateFitness, population.top_fitness[i:end-1,:])
  end
end


text:
module ModiaInterface

export Object3D, Scene, Visual, Solid
export Box, Beam, Cylinder, Sphere, Ellipsoid
export Cone, Capsule, GearWheel, Grid, SpringShape
export CoordinateSystem, FileMesh
export Font, TextShape
export VisualMaterial
export MassProperties, MassPropertiesFromShape
export MassPropertiesFromShapeAndMass
export UniformGravityField
export RefPath, ptpJointSpace, scheduleReferenceMotion
export calculateRobotMovement
export getRefPathPosition, getRefPathInitPosition, getVariables
export multibodyResiduals
export Fix
export Revolute, RevoluteWithFlange
export Prismatic, PrismaticWithFlange
export J123, J132, FreeMotion

export buildModia3D

import Modia3D
using ModiaLang

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")
include("$(ModiaLang.path)/models/Translational.jl")

include("buildModia3D.jl")
include("model3D.jl")

end


text:
if abspath(PROGRAM_FILE) == @__FILE__
    println("OK")
    quit(1)
end


text:
using ImageCore, Colors, FixedPointNumbers, ColorVectorSpace
using Base.Test

@testset "Deprecated" begin
    @testset "constructors" begin
        # grayim
        for a in (0x00:0x18, collect(0x00:0x18))
            B = reshape(a, 5, 5)
            img = grayim(B)
            @test colorspace(img) == "Gray"
            @test colordim(img) == 0
            @test img[5,2] === Gray{U8}(9/255)
        end
        for a in (0x00:0x7c, collect(0x00:0x7c))
            B = reshape(a, 5, 5, 5)
            img = grayim(B)
            @test colorspace(img) == "Gray"
            @test colordim(img) == 0
            @test img[5,2,3] === Gray{U8}(59/255)
        end
        for a in (0x0000:0x0018, collect(0x0000:0x0018))
            B = reshape(a, 5, 5)
            img = grayim(B)
            @test colorspace(img) == "Gray"
            @test colordim(img) == 0
            @test img[5,2] === Gray{U16}(9/65535)
        end
        for a in (0x0000:0x007c, collect(0x0000:0x007c))
            B = reshape(a, 5, 5, 5)
            img = grayim(B)
            @test colorspace(img) == "Gray"
            @test colordim(img) == 0
            @test img[5,2,3] === Gray{U16}(59/65535)
        end
        # colorim
        for a in (0x00:0x4a, collect(0x00:0x4a))
            C = reshape(a, 3, 5, 5)
            img = colorim(C)
            @test colorspace(img) == "RGB"
            @test colordim(img) == 0
            @test img[1,1] === RGB{U8}(0,1/255,2/255)
            C = reshape(a, 5, 5, 3)
            img = colorim(C)
            @test colorspace(img) == "RGB"
            @test colordim(img) == 0
            @test img[1,1] === RGB{U8}(0,25/255,50/255)
            @test_throws ErrorException colorim(reshape(a, 5, 3, 5))
        end
        for a in (0x0000:0x004a, collect(0x0000:0x004a))
            C = reshape(a, 3, 5, 5)
            img = colorim(C)
            @test colorspace(img) == "RGB"
            @test colordim(img) == 0
            @test img[1,1] === RGB{U16}(0,1/65535,2/65535)
            C = reshape(a, 5, 5, 3)
            img = colorim(C)
            @test colorspace(img) == "RGB"
            @test colordim(img) == 0
            @test img[1,1] === RGB{U16}(0,25/65535,50/65535)
            @test_throws ErrorException colorim(reshape(a, 5, 3, 5))
        end
        for a in (0x00:0x63, collect(0x00:0x63))
            for (S,T) in (("RGBA", RGBA), ("ARGB",ARGB))
                C = reshape(a, 4, 5, 5)
                img = colorim(C, S)
                @test colorspace(img) == S
                @test colordim(img) == 0
                @test img[1,1] === T{U8}(0,1/255,2/255,3/255)
                C = reshape(a, 5, 5, 4)
                img = colorim(C, S)
                @test colorspace(img) == S
                @test colordim(img) == 0
                @test img[1,1] === T{U8}(0,25/255,50/255,75/255)
            end
            @test_throws ErrorException colorim(reshape(a, 4, 5, 5), "RRRR")
        end
        # ambiguous sizes
        for a in (0x00:0x2c, collect(0x00:0x2c))
            C = reshape(a, 3, 5, 3)
            @test_throws ErrorException colorim(C)
        end
        for a in (0x00:0x4f, collect(0x00:0x4f))
            C = reshape(a, 4, 5, 4)
            @test_throws ErrorException colorim(C)
        end
    end
    @testset "traits" begin
        for (B,S) in ((rand(UInt16(1):UInt16(20), 3, 5),"Gray"),
                      (rand(Gray{Float32}, 3, 5),"Gray"),
                      (rand(RGB{Float16}, 3, 5),"RGB"),
                      (bitrand(3, 5),"Binary"),
                      (rand(UInt32, 3, 5),"RGB24"))
            T = eltype(B)
            @test data(B) === B
            @test colorspace(B) == S
            @test colordim(B) == 0
            @test timedim(B) == 0
            @test spatialorder(B) == [:y,:x]
            @test isdirect(B)
            @test limits(B) == (T==Bool ? (0,1) : (zero(T),one(T)))
            @test storageorder(B) == [:y,:x]
            @test ncolorelem(B) == (T <: Colorant ? length(T) : 1)
            assert2d(B)
            assert_scalar_color(B)
            @test isyfirst(B)
            @test !isxfirst(B)
            assert_yfirst(B)
            @test_throws ErrorException assert_xfirst(B)
            @test spatialproperties(B) == String[]
            @test_throws ErrorException permutedims(B, ("x", "y"))
            @test_throws ErrorException permutedims(B, (:x, :y))
        end
        for (B,S) in ((rand(UInt16(1):UInt16(20), 5),"Gray"),
                      (rand(Gray{Float32}, 5),"Gray"),
                      (rand(RGB{Float16}, 5),"RGB"),
                      (bitrand(5),"Binary"))
            @test colorspace(B) == S
        end
        for (B,S) in ((rand(Gray{Float32}, 5, 5, 5),"Gray"),
                      (rand(RGB{Float16}, 5, 5, 5),"RGB"),
                      (bitrand(5, 5, 5),"Binary"))
            @test colorspace(B) == S
        end
        @test colorspace(rand(Float32, 5, 5, 3)) == "RGB"
        @test_throws ErrorException colorspace(rand(Float32, 5, 5, 5))
    end
    @testset "Reinterpret, convert, separate, raw" begin
        # some of these are redundant with convert_reinterpret.jl, but
        # these are "classic" tests so it's good to make sure they
        # pass
        a = RGB{Float64}[RGB(1,1,0)]
        af = reinterpret(Float64, a)
        @test vec(af) == [1.0,1.0,0.0]
        @test size(af) == (3,1)
        @test_throws DimensionMismatch reinterpret(Float32, a)
        anew = reinterpret(RGB, af)
        @test anew == a
        anew = reinterpret(RGB, vec(af))
        @test anew[1] == a[1]
        @test ndims(anew) == 0
        anew = reinterpret(RGB{Float64}, af)
        @test anew == a
        @test_throws DimensionMismatch reinterpret(RGB{Float32}, af)
        Au8 = rand(0x00:0xff, 3, 5, 4)
        A8 = reinterpret(UFixed8, Au8)
        rawrgb8 = reinterpret(RGB, A8)
        @test eltype(rawrgb8) == RGB{UFixed8}
        @test reinterpret(UFixed8, rawrgb8) == A8
        @test reinterpret(UInt8, rawrgb8) == Au8
        rawrgb32 = convert(Array{RGB{Float32}}, rawrgb8)
        @test eltype(rawrgb32) == RGB{Float32}
        @test ufixed8(rawrgb32) == rawrgb8
        @test reinterpret(UFixed8, rawrgb8) == A8
        # imrgb8 = convert(Image, rawrgb8)
        # @test spatialorder(imrgb8) == ImageCore.yx
        # cvt = convert(Image, imrgb8)
        # @test cvt == imrgb8 && typeof(cvt) == typeof(imrgb8)
        # cvt = convert(Image{RGB{UFixed8}}, imrgb8)
        # @test cvt == imrgb8 && typeof(cvt) == typeof(imrgb8)
        imrgb8 = rawrgb8
        im8 = reinterpret(UFixed8, imrgb8)
        @test data(im8) == A8
        @test reinterpret(UInt8, imrgb8) == Au8
        @test reinterpret(RGB, im8) == imrgb8
        ims8 = separate(imrgb8)
        @test colorspace(ims8) == "RGB"
        # cvt = convert(Image, ims8)
        # @test cvt == ims8 && typeof(cvt) == typeof(ims8)
        # @test (cvt = convert(Image{UFixed8}, ims8)) == ims8 && typeof(cvt) == typeof(ims8)
        # @test (cvt = separate(ims8)) == ims8 && typeof(cvt) == typeof(ims8)
        # imrgb8_2 = convert(Image{RGB}, ims8)
        # @test isa(imrgb8_2, Image{RGB{UFixed8}})
        A = reinterpret(UFixed8, UInt8[1 2; 3 4])
        # imgray = convert(Image{Gray{UFixed8}}, A)
        imgray = convert(Array{Gray{U8}}, A)
        @test spatialorder(imgray) == ImageCore.yx
        @test data(imgray) == reinterpret(Gray{UFixed8}, [0x01 0x02; 0x03 0x04])
        # @test eltype(convert(Image{HSV{Float32}}, imrgb8)) == HSV{Float32}
        # @test eltype(convert(Image{HSV}, float32(imrgb8))) == HSV{Float32}
        # local img = Image(reinterpret(Gray{UFixed16}, rand(UInt16, 5, 5)))
        # imgs = subim(img, :, :)
        img = reinterpret(Gray{UFixed16}, rand(UInt16, 5, 5))
        imgs = view(img, :, :)
        @test isa(minimum(imgs), Gray{UFixed16})  # change in behavior from Images issue#232
        imgdata = rand(UInt16, 5, 5)
        # img = Image(reinterpret(Gray{UFixed16}, imgdata))
        img = reinterpret(Gray{UFixed16}, imgdata)
        @test all(raw(img) .== imgdata)
        # @test size(raw(Image(rawrgb8))) == (3,5,4)
        @test all(raw(imgdata) .== imgdata)
    end
end

nothing


text:
"""

# OptimizeAlgorithm types

### Types defined
```julia
* Lbfgs::OptimizeAlgorithm   : Euclidean manifold with unit metric
* Bfgs::OptimizeAlgorithm    : Euclidean manifold with unit metric
* Newton::OptimizeAlgorithm  : Euclidean manifold with unit metric
```
"""
abstract type OptimizeAlgorithm end

"""

# Lbfgs type and constructor

Settings for method=Lbfgs() in Optimize().

### Method
```julia
Lbfgs(;init_alpha=0.001, tol_obj=1e-8, tol_grad=1e-8, tol_param=1e-8, history_size=5)
```
### Optional arguments
```julia
* `init_alpha::Float64`        : Linear search step size for first iteration
* `tol_obj::Float64`           : Convergence tolerance for objective function
* `tol_rel_obj::Float64`       : Relative change tolerance in objective function
* `tol_grad::Float64`          : Convergence tolerance on norm of gradient
* `tol_rel_grad::Float64`      : Relative change tolerance on norm of gradient
* `tol_param::Float64`         : Convergence tolerance on parameter values
* `history_size::Int`          : No of update vectors to use in Hessian approx
```

### Related help
```julia
?OptimizeMethod               : List of available optimize methods
?Optimize                      : Optimize arguments
```
"""
mutable struct Lbfgs <: OptimizeAlgorithm
  init_alpha::Float64
  tol_obj::Float64
  tol_rel_obj::Float64
  tol_grad::Float64
  tol_rel_grad::Float64
  tol_param::Float64
  history_size::Int64
end
Lbfgs(;init_alpha=0.001, tol_obj=1e-8,  tol_rel_obj=1e4,
  tol_grad=1e-8, tol_rel_grad=1e7, tol_param=1e-8, history_size=5) =
    Lbfgs(init_alpha, tol_obj, tol_rel_obj, tol_grad, tol_rel_grad,
      tol_param, history_size)

"""

# Bfgs type and constructor

Settings for method=Bfgs() in Optimize().

### Method
```julia
Bfgs(;init_alpha=0.001, tol_obj=1e-8, tol_rel_obj=1e4,
  tol_grad=1e-8, tol_rel_grad=1e7, tol_param=1e-8)
```
### Optional arguments
```julia
* `init_alpha::Float64`        : Linear search step size for first iteration
* `tol_obj::Float64`           : Convergence tolerance for objective function
* `tol_rel_obj::Float64`       : Relative change tolerance in objective function
* `tol_grad::Float64`          : Convergence tolerance on norm of gradient
* `tol_rel_grad::Float64`      : Relative change tolerance on norm of gradient
* `tol_param::Float64`         : Convergence tolerance on parameter values
```

### Related help
```julia
?OptimizeMethod                : List of available optimize methods
?Optimize                      : Optimize arguments
```
"""
mutable struct Bfgs <: OptimizeAlgorithm
  init_alpha::Float64
  tol_obj::Float64
  tol_rel_obj::Float64
  tol_grad::Float64
  tol_rel_grad::Float64
  tol_param::Float64
end
Bfgs(;init_alpha=0.001, tol_obj=1e-8, tol_rel_obj=1e4,
  tol_grad=1e-8, tol_rel_grad=1e7, tol_param=1e-8) =
    Bfgs(init_alpha, tol_obj, tol_rel_obj, tol_grad, tol_rel_grad, tol_param)

"""

# Newton type and constructor

Settings for method=Newton() in Optimize().

### Method
```julia
Newton()
```
### Related help
```julia
?OptimizeMethod                : List of available optimize methods
?Optimize                      : Optimize arguments
```
"""
mutable struct Newton <: OptimizeAlgorithm
end

"""

# Optimize type and constructor

Settings for Optimize top level method.

### Method
```julia
Optimize(;
  method=Lbfgs(),
  iter=2000,
  save_iterations=false
)
```
### Optional arguments
```julia
* `method::OptimizeMethod`      : Optimization algorithm
* `iter::Int`                   : Total number of iterations
* `save_iterations::Bool`       : Stream optimization progress to output
```

### Related help
```julia
?Stanmodel                      : Create a StanModel
?OptimizeAlgorithm              : Available algorithms
```
"""
mutable struct Optimize <: Method
  method::OptimizeAlgorithm
  iter::Int64
  save_iterations::Bool
end
Optimize(;method::OptimizeAlgorithm=Lbfgs(), iter::Number=2000,
  save_iterations::Bool=false) =
  Optimize(method, iter, save_iterations)
Optimize(method::OptimizeAlgorithm) = Optimize(method, 2000, false)

function optimize_show(io::IO, o::Optimize, compact::Bool)
  if compact
    println("Optimize($(o.method), $(o.iter), $(o.save_iterations))")
  else
    println("  method =                  Optimize()")
    if isa(o.method, Lbfgs)
      println("    algorithm =               Lbfgs()")
      println("      init_alpha =              ", o.method.init_alpha)
      println("      tol_obj =                 ", o.method.tol_obj)
      println("      tol_grad =                ", o.method.tol_grad)
      println("      tol_param =               ", o.method.tol_param)
      println("      history_size =            ", o.method.history_size)
    elseif isa(o.method, Bfgs)
      println("    algorithm =               Bfgs()")
      println("      init_alpha =              ", o.method.init_alpha)
      println("      tol_obj =                 ", o.method.tol_obj)
      println("      tol_grad =                ", o.method.tol_grad)
      println("      tol_param =               ", o.method.tol_param)
    else
      println("    algorithm =               Newton()")
    end
    println("    iterations =              ", o.iter)
    println("    save_iterations =         ", o.save_iterations)
  end
end

show(io::IO, o::Optimize) = optimize_show(io, o, false)


text:
## tape compilation

## CONVERSION TO EXGRAPH

make_name(op::AbstractOp) = Symbol("%$(getid(getvar(op)))")
make_name(var::TAny) = Symbol("%$(getid(var))")

to_exnode(op::Input) = ExNode{:input}(make_name(op), make_name(op); val=getvalue(op))
to_exnode(op::Constant) = ExNode{:constant}(make_name(op), getvalue(op); val=getvalue(op))
to_exnode(op::Assign) = ExNode{:(=)}(make_name(op), make_name(op.src); val=getvalue(op))

function to_exnode(op::Call)
    arg_names = map(make_name, op.args)
    # fns = maybe_to_symbol(op.fn)
    if isempty(op.kwargs)
        ex = Expr(:call, op.fn, arg_names...)
    else
        ex = Expr(:call, op.fn, Espresso.make_kw_params(op.kwargs), arg_names...)
    end
    return ExNode{:call}(make_name(op), ex; val=getvalue(op))
end

function to_exnode(op::Bcast)
    arg_names = map(make_name, op.args)
    #  fns = maybe_to_symbol(op.fn)
    ex = Expr(:., op.fn, Expr(:tuple, arg_names...))
    ExNode{:bcast}(make_name(op), ex; val=getvalue(op))
end

function to_exgraph(tape::Tape)
    g = ExGraph()
    for op in tape
        push!(g, to_exnode(op))
    end
    return g
end

## CODE GENERATION

"""
Rewrite :(Z = X * Y) into :(mul!(Z, X, Y)),  but only if X and Y are arrays
"""
function rewrite_mul(tape::Tape, ex::Expr)
    @assert ex.head == :block
    new_ex = Expr(:block)
    for subex in ex.args
        if (matchingex(:(_R = $*(_X, _Y)), subex)
            && tape[subex.args[1]].val isa AbstractArray
            && all(tape[subex.args[2].args[i]].val isa AbstractArray for i=2:3))
            new_subex = rewrite(subex, :(_Z = $*(_X, _Y)), :($mul!(_Z, _X, _Y)))
            push!(new_ex.args, new_subex)
        else
            push!(new_ex.args, subex)
        end
    end
    return new_ex
end

"""
Generate initialization block of function code
"""
function generate_prologue(tape::Tape)
    ex = Expr(:body)
    for op in tape
        name = make_name(op)
        push!(ex.args, :($name = $(op.var).val))
    end
    return ex
end

"""
Main part of generated function code
"""
function generate_body(tape::Tape)
    ret_var_ids = vcat(values(tape.derivs) |> collect, [tape.resultid])
    ret_var_names = [Symbol("%$id") for id in ret_var_ids]
    exg = to_exgraph(tape)
    # exg = Espresso.fuse_broadcasting(exg)  # TODO: this breaks test_examples.jl
    exg = Espresso.fuse_assigned(exg; outvars=ret_var_names)
    exg = Espresso.eliminate_common(exg)
    ex = Espresso.to_expr_kw(exg)
    # ex = rewrite_all(ex, INPLACE_RULES; phs=Set([:X, :Y, :Z]))
    ex = rewrite_mul(tape, ex)
    return ex
end

"""
Generate last block of function code
"""
function generate_epilogue(tape::Tape)
    ex = Expr(:body)
    for op in tape
        name = make_name(op)
        push!(ex.args, :($(op.var).val = $name))
    end
    return ex
end

function generate_function_expr(tape::Tape)
    fn_ex = :(function $(gensym("tape_fn"))() end)
    fn_ex_body = fn_ex.args[2]
    prologue = generate_prologue(tape)
    body = generate_body(tape)
    epilogue = generate_epilogue(tape)
    # abusing LineNumberNode to insert comments
    push!(fn_ex_body.args, LineNumberNode(0, "prologue"))
    push!(fn_ex_body.args, prologue.args...)
    push!(fn_ex_body.args, LineNumberNode(0, "body"))
    push!(fn_ex_body.args, body.args...)
    push!(fn_ex_body.args, LineNumberNode(0, "epilogue"))
    push!(fn_ex_body.args, epilogue.args...)
    return fn_ex
end

## COMPILATION

function compile(tape::Tape)
    fn_ex = generate_function_expr(tape)
    return Core.eval(@__MODULE__, fn_ex)
end

function compile!(tape::Tape)
    tape.compiled = compile(tape)
end

function rerecord_inputs!(tape::Tape, args...)
    minitape = Tape()
    targs = make_tracked_args(minitape, args...)
    for i=1:length(minitape)
        val = getvalue(minitape[i])
        setvalue!(tape[i], val)
    end
end

function play!(tape::Tape, args...; use_compiled=true)
    rerecord_inputs!(tape, args...)
    if use_compiled && tape.compiled != nothing
        Base.invokelatest(tape.compiled)
    else
        for op in tape
            exec!(tape, op)
        end
    end
end


text:
module TestSynthesis
using CustomTest
using Base.Test
using ControlSystems

P = tf(1,[1,1])
C = tf([1,1],[1,0])
L = P*C
Lsys = ss(L)

B = [1]
A = [1,1]
R = [1,1]
S = [1]
T = [1]

@test isapprox(minreal(feedback(P,C),1e-5), minreal(feedback(L),1e-5))
@test isapprox(numpoly(minreal(feedback(L),1e-5))[1].a, numpoly(tf(1,[1,1]))[1].a)# This test is ugly, but numerical stability is poor for minreal
@test feedback2dof(B,A,R,S,T) == tf(B.*T, conv(A,R) + [0;0;conv(B,S)])
@test feedback2dof(P,R,S,T) == tf(B.*T, conv(A,R) + [0;0;conv(B,S)])
@test isapprox(pole(minreal(tf(feedback(Lsys)),1e-5)) , pole(minreal(feedback(L),1e-5)), atol=1e-5)

@test_throws ErrorException feedback(ss(1),ss(1))
@test_throws ErrorException feedback(ss(eye(2), ones(2,2), ones(1,2),0))

end


text:
using Test
using Pat

@testset "Basic Construction" begin
    cb = Pat.CelestialBody(1.3, 1)
    @test Pat.gravitational_parameter(cb) == 1.3
    @test Pat.naif_id(cb) == 1

    @test Pat.gravitational_parameter(Pat.CelestialBody(12.3)) == 12.3
    @test Pat.mass(Pat.CelestialBody(6.6743e-20)) == 1
end


text:
using BitBasis
using YaoArrayRegister: SDDiagonal, ArrayReg, staticize
using TupleTools
using StaticArrays
using NiLang.AD: GVar

export i_instruct!, VecStack, incstack!, stack4reg

"""
    VecStack{T}

A vector stack, its data is initialized as an empty Matrix.
It has a stack `top` that points to one column of the matrix.
The `VecStack` then behaves likes a normal vector.
One can use `incstack!` to increase the stack top.
"""
struct VecStack{T}
	data::Matrix{T}
	top::Int
end

Base.getindex(x::VecStack, i::Int) = @inbounds x.data[i,x.top]
Base.setindex!(x::VecStack, val, i::Int) = @inbounds setindex!(x.data, val, i, x.top)
Base.copy(vs::VecStack) = VecStack(copy(vs.data), vs.top)

"""
    stack4reg(reg, n)

Build a `VecStack` of size (`n` × the register size).
"""
function stack4reg(reg::ArrayReg{1,T}, n::Int) where T<:Tropical
	VecStack(ones(T, length(reg.state), n), 0)
end

"""
    incstack!(vs::VecStack)

Increase the stack size for 1.
"""
@i function incstack!(vs::VecStack)
	@safe @assert size(vs.data, 2) > vs.top
	INC(vs.top)
end

@i @inline function unsafe_store!(neginf::Tropical{T}, x::Tropical{T}) where T
    neginf *= x
end

"""
    i_instruct!(state, U0, locs, clocs, cvals, REG_STACK)

Reversible instruction function.

* `state`is a register with tropical numbers as elements.
* `U0` is the gate matrix, also with tropical numbers as elements.
* `locs` is the locations (tuple of ints) to apply gate `U0`.
* `clocs` is the locations of control bits.
* `cvals` is the values of control bits, should be a tuple of `0` and `1`.
* `REG_STACK` is the stack for caching intermediate state.
    If U0 is dense, one-register-sized space is required,
    if `U0` is diagonal, it can be empty due to the reversibility.
    For stack constructor, see `stack4reg` for details.
"""
@i function i_instruct!(state::StridedVector{T},
        U0::AbstractMatrix{T},
        locs::NTuple{M, Int},
        clocs::NTuple{C, Int},
        cvals::NTuple{C, Int},
		REG_STACK::VecStack{T}) where {T<:Tropical,C, M}
    @routine @invcheckoff begin
        nbit ← log2i(length(state))
        if (!all(TupleTools.diff(locs).>0), ~)
            @safe error("Sorry, `TropicalYao` only supports the location specification in increasing order. File an issue if you need it!")
        end
        MM ← size(U0, 1)
        locked_bits ← [clocs..., locs...]
        locked_vals ← [cvals..., zeros(Int, M)...]
        locs_raw ← [i+1 for i in itercontrol(nbit, setdiff(1:nbit, locs), zeros(Int, nbit-M))] |> staticize
        configs ← itercontrol(nbit, locked_bits, locked_vals)
    end
    loop_kernel(state, configs, U0, locs_raw, REG_STACK)
    ~@routine
end

@i function i_instruct!(state::StridedVector{T},
		val::Val{:X},
        locs::NTuple{M, Int},
        clocs::NTuple{C, Int},
        cvals::NTuple{C, Int},
		REG_STACK::VecStack{T}) where {T<:Tropical,C, M}
	@routine @invcheckoff begin
	    configs ← itercontrol(log2dim1(state), [clocs..., locs[1]], [cvals..., 0])
	    mask ← bmask(locs...)
	end
	@invcheckoff for i=1:length(configs)
		@invcheckoff b ← configs[i]
        @inbounds NiLang.SWAP(state[b+1], state[flip(b, mask) + 1])
		@invcheckoff b → configs[i]
    end
	~@routine
end

@i function loop_kernel(state::StridedVector{<:Tropical}, configs, U::SDDiagonal, locs_raw, REG_STACK)
    @invcheckoff @inbounds for i=1:length(configs)
        x ← configs[i]
        for i=1:length(U.diag)
            state[x+locs_raw[i]] *= U.diag[i]
        end
        x → configs[i]
    end
end

_scache(u::MMatrix{N,N,T}) where {N,T} = MVector{N}(ones(T, N))
_scache(u::AbstractMatrix{T}) where {T} = MVector{size(u,1)}(ones(T, size(u, 1)))

@i function loop_kernel(state::StridedVector{T}, configs, U::AbstractMatrix, locs_raw, REG_STACK) where {T<:Tropical}
    @invcheckoff begin
        @routine begin
            branch_keeper ← zeros(Bool, size(U, 2))
            scache ← _scache(U)
        end
		incstack!(REG_STACK)
        @inbounds for i=1:length(configs)
			@routine begin
            	x ← configs[i]
				for l=1:length(locs_raw)
					scache[l] *= state[x + locs_raw[l]]
				end
			end
			for l=1:size(U,1)
				el ← zero(T)
				@routine for k=1:size(U,2)
					scache[k] *= U[l,k]
					if (el.n < scache[k].n, branch_keeper[k])
						FLIP(branch_keeper[k])
						NiLang.SWAP(el, scache[k])
					end
				end
                unsafe_store!(REG_STACK[x+locs_raw[l]], el)
				~@routine
				el → zero(T)
			end
			~@routine
        end
		for i=1:length(state)
        	NiLang.SWAP(state[i], REG_STACK[i])
		end
        ~@routine
    end
end

@i function loop_kernel(state::StridedVector{T}, configs, U::MMatrix{2,2}, locs_raw, REG_STACK) where T<:Tropical
    @invcheckoff begin
		@routine @inbounds begin
			a ← one(T)
			b ← one(T)
			c ← one(T)
			d ← one(T)
			a *= U[1]
			b *= U[2]
			c *= U[3]
			d *= U[4]
		end
		incstack!(REG_STACK)
        @inbounds for i=1:length(configs)
            @routine begin
                x ← configs[i]
                is1 ← x + locs_raw[1]
                is2 ← x + locs_raw[2]
                s1 ← one(T)
                s2 ← one(T)
            end
			NiLang.SWAP(s1, state[is1])
			NiLang.SWAP(s2, state[is2])
			@routine begin
				s3 ← one(T)
				s4 ← one(T)
				s3 *= s1 * b
				s4 *= s2 * d
				s1 *= a
				s2 *= c
				bk1 ← s1 > s2
				bk2 ← s3 > s4
			end
			if (bk1, ~)
				state[is1] *= s1
			else
				state[is1] *= s2
			end
			if (bk2, ~)
				state[is2] *= s3
			else
				state[is2] *= s4
			end
			~@routine
        	NiLang.SWAP(s1, REG_STACK[is1])
        	NiLang.SWAP(s2, REG_STACK[is2])
            ~@routine
    	end
		~@routine
	end
end

NiLang.AD.GVar(x::ArrayReg{B}) where B = ArrayReg{B}(GVar(x.state))
NiLang.AD.GVar(x::VecStack) where B = VecStack(GVar(x.data), x.top)


text:
"""
    App(handler)

calls handler with the session and the http request object.
f is expected to return a valid DOM object,
so something renderable by jsrender, e.g. `DOM.div`.
"""
struct App
    handler::Function
    function App(handler::Function)
        if hasmethod(handler, Tuple{Session, HTTP.Request})
            return new(handler)
        elseif hasmethod(handler, Tuple{Session})
            return new((session, request) -> handler(session))
        elseif hasmethod(handler, Tuple{HTTP.Request})
            return new((session, request) -> handler(request))
        elseif hasmethod(handler, Tuple{})
           return new((session, request) -> handler())
        else
            error("""
            Handler function must have the following signature:
                handler() -> DOM
                handler(session::Session) -> DOM
                handler(request::Request) -> DOM
                handler(session, request) -> DOM
            """)
        end
    end
    function App(dom_object)
        return new((s, r)-> dom_object)
    end
end

"""
The string part of JSCode.
"""
struct JSString
    source::String
end

"""
Javascript code that supports interpolation of Julia Objects.
Construction of JSCode via string macro:
```julia
jsc = js"console.log(\$(some_julia_variable))"
```
This will decompose into:
```julia
jsc.source == [JSString("console.log("), some_julia_variable, JSString("\"")]
```
"""
struct JSCode
    source::Vector{Union{JSString, Any}}
end

"""
Represent an asset stored at an URL.
We try to always have online & local files for assets
If one gives an online resource, it will be downloaded, to host it locally.
"""
struct Asset
    media_type::Symbol
    # We try to always have online & local files for assets
    # If you only give an online resource, we will download it
    # to also be able to host it locally
    online_path::String
    local_path::String
    onload::Union{Nothing, JSCode}
end

"""
Encapsulates frontend dependencies. Can be used in the following way:

```Julia
const noUiSlider = Dependency(
    :noUiSlider,
    # js & css dependencies are supported
    [
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.js",
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.css"
    ]
)
# use the dependency on the frontend:
evaljs(session, js"\$(noUiSlider).some_function(...)")
```
jsrender will make sure that all dependencies get loaded.
"""
struct Dependency
    name::Symbol # The JS Module name that will get loaded
    assets::Vector{Asset}
    # The global -> Function name, JSCode -> the actual function code!
    functions::Dict{Symbol, JSCode}
end

"""
    UrlSerializer
Struct used to encode how an url is rendered
Fields:
```julia
# uses assetserver?
assetserver::Bool
# if assetserver == false, we move all assets into asset_folder
# for someone else to serve them!
asset_folder::Union{Nothing, String}

absolute::Bool
# Used to prepend if absolute == true
content_delivery_url::String
```
"""
struct UrlSerializer
    # uses assetserver?
    assetserver::Bool
    # if assetserver == false, we move all assets into asset_folder
    # for someone else to serve them!
    asset_folder::Union{Nothing, String}

    absolute::Bool
    # Used to prepend if absolute == true
    content_delivery_url::String
    # Inlines all content directly into the html
    # Makes all above options obsolete
    inline_all::Bool
end

function UrlSerializer(;
        proxy = JSSERVE_CONFIGURATION.content_delivery_url[],
        assetserver=true, asset_folder=nothing, absolute=proxy!="",
        inline_all=false
    )
    return UrlSerializer(
        assetserver, asset_folder, absolute, proxy, inline_all
    )
end

struct JSException <: Exception
    exception::String
    message::String
    stacktrace::Vector{String}
end

"""
Creates a Julia exception from data passed to us by the frondend!
"""
function JSException(js_data::AbstractDict)
    stacktrace = String[]
    if js_data["stacktrace"] !== nothing
        for line in split(js_data["stacktrace"], "\n")
            push!(stacktrace, replace(line, ASSET_URL_REGEX => replace_url))
        end
    end
    return JSException(js_data["exception"], js_data["message"], stacktrace)
end

function Base.show(io::IO, exception::JSException)
    println(io, "An exception was thrown in JS: $(exception.exception)")
    println(io, "Additional message: $(exception.message)")
    println(io, "Stack trace:")
    for line in exception.stacktrace
        println(io, "    ", line)
    end
end

"""
A web session with a user
"""
struct Session
    # IOBuffer for testing
    connection::Base.RefValue{Union{Nothing, WebSocket, IOBuffer}}
    # Bool -> if already registered with Frontend
    observables::Dict{String, Tuple{Bool, Observable}}
    message_queue::Vector{Dict{Symbol, Any}}
    dependencies::Set{Asset}
    on_document_load::Vector{JSCode}
    id::String
    js_fully_loaded::Channel{Bool}
    on_websocket_ready::Any
    url_serializer::UrlSerializer
    # Should be checkd on js_fully_loaded to see if an error occured
    init_error::Ref{Union{Nothing, JSException}}
    js_comm::Observable{Union{Nothing, Dict{String, Any}}}
    on_close::Observable{Bool}
    deregister_callbacks::Vector{Observables.ObserverFunction}
    unique_object_cache::Dict{String, WeakRef}
end

struct Routes
    table::Vector{Pair{Any, Any}}
end

"""
The application one serves
"""
struct Server
    url::String
    port::Int
    sessions::Dict{String, Session}
    server_task::Ref{Task}
    server_connection::Ref{TCPServer}
    routes::Routes
    websocket_routes::Routes
end


text:
using Images, MXNet, JLD

nnet = mx.load_checkpoint("weights/inception-v3/InceptionV3-FE", 0, mx.FeedForward; context = mx.gpu());

PROJECT_IMG_FOLDER = joinpath("data", "book-covers");
DATASETS = ["train", "test"];
BATCH_SIZE = 100;

for dataset_name in DATASETS

    categories = Int[]

    current_folder = joinpath(PROJECT_IMG_FOLDER, dataset_name)
    file_name = joinpath(PROJECT_IMG_FOLDER, "bookcover30-labels-$dataset_name.txt")
    lines = readlines(file_name);

    results = zeros(Float16, (2048, size(lines, 1)))

    for i=1:BATCH_SIZE:size(lines, 1)

        println("BATCH: $i")

        mx_data = mx.zeros((299, 299, 3, BATCH_SIZE));

        for j=1:BATCH_SIZE

            if j % 50 == 0
                println("\tIMG: $j")
            end

            img = imresize(load(joinpath(current_folder, split(lines[i + j - 1])[1])), (299, 299));
            img = permutedims(Float16.(channelview(RGB.(img))), (3, 2, 1));

            img[:, :, :] *= 256.0;
            img[:, :, :] -= 128.;
            img[:, :, :] /= 128.;

            img = reshape(img, 299, 299, 3, 1);
            mx_data[j:j] = img;
        end

        data_provider = mx.ArrayDataProvider(:data => mx_data);
        results[:, i:i+(BATCH_SIZE-1)] = mx.predict(nnet, data_provider)
    end

    JLD.save(joinpath(PROJECT_IMG_FOLDER, "dump-$dataset_name-inception.jld"), "data", results)
end


text:
using KernelDensity
using PyPlot

KernelDensity.kde(t::BenchmarkTools.Trial; kwargs...) = kde(t.times; kwargs...)
PyPlot.plot(t::BenchmarkTools.Trial; kwargs...) = scatter(0:length(t)-1, t.times; kwargs...)
PyPlot.plot(k::KernelDensity.UnivariateKDE; kwargs...) = PyPlot.plot(k.x, k.density; kwargs...)

function plotkde(trials::Vector; kwargs...)
    for t in trials
        plotkde(t; kwargs...)
    end
end

function plotkde(trial::BenchmarkTools.Trial; cut = 0.1, bandwidth = nothing, kwargs...)
    trial = BenchmarkTools.trim(trial, cut)
    k = bandwidth == nothing ? kde(trial) : kde(trial; bandwidth = bandwidth)
    plot(k; kwargs...)
end


text: println(sum(filter(x -> mod(x, 3) == 0 || mod(x, 5) == 0, 1:999)))

max_stars_repo_name: seaneshbaugh/rosetta-euler

max_issues_repo_name: seaneshbaugh/rosetta-euler

max_forks_repo_name: seaneshbaugh/rosetta-euler


text:
## MPI Init
using MPIClusterManagers, Distributed
import MPI
MPI.Init()
rank = MPI.Comm_rank(MPI.COMM_WORLD)
size = MPI.Comm_size(MPI.COMM_WORLD)
if rank == 0
    @info "size is $size"
end
manager = MPIClusterManagers.start_main_loop(MPI_TRANSPORT_ALL)
@info "there are $(nworkers()) workers"
@everywhere @info gethostname()
include("01_read_data.jl")
include("02_set_options.jl")
## MPI checks
# split into sequential iterations of parallel soundings
nsoundings = length(soundings)
ncores = nworkers()
@assert mod(ncores+1,nchainspersounding+1) == 0
@assert mod(ppn, nchainspersounding+1) == 0
nparallelsoundings = Int((ncores+1)/(nchainspersounding+1))
nsequentialiters = ceil(Int, nsoundings/nparallelsoundings)
@info "will require $nsequentialiters iterations of $nparallelsoundings"
## set up McMC
@everywhere using Distributed
@everywhere using HiQGA.transD_GP
## do the parallel soundings
@info "starting"
transD_GP.SkyTEM1DInversion.loopacrosssoundings(soundings, opt;
                    nsequentialiters   = nsequentialiters,
                    nparallelsoundings = nparallelsoundings,
                    zfixed             = zfixed,
                    ρfixed             = ρfixed,
                    zstart             = zstart,
                    extendfrac         = extendfrac,
                    dz                 = dz,
                    ρbg                = ρbg,
                    nlayers            = nlayers,
                    ntimesperdecade    = ntimesperdecade,
                    nfreqsperdecade    = nfreqsperdecade,
                    Tmax               = Tmax,
                    nsamples           = nsamples,
                    nchainsatone       = nchainsatone,
                    nchainspersounding = nchainspersounding,
                    useML              = useML)

MPIClusterManagers.stop_main_loop(manager)
rmprocs(workers())
exit()

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
order = ([
    "Nuclear",
    "Coal",
    "Hydro",
    "Gas_CC",
    "Gas_CT",
    "Storage",
    "Oil_ST",
    "Oil_CT",
    "Sync_Cond",
    "Wind",
    "Solar",
    "CSP",
    "curtailment",
])
function _get_iterator(sys::PSY.System, results::IS.Results)
    iterators = []
    for (k, v) in results.variables
        if "$k"[1:2] == "P_"
            datatype = (split("$k", "P_")[2])
            if datatype == "ThermalStandard"
                iterators =
                    vcat(iterators, collect(PSY.get_components(PSY.ThermalStandard, sys)))
            elseif datatype == "RenewableDispatch"
                iterators =
                    vcat(iterators, collect(PSY.get_components(PSY.RenewableDispatch, sys)))
            end
        end
    end
    iterators_sorted = Dict{Any, Any}()
    for iterator in iterators
        name = iterator.name
        iterators_sorted[name] = []
        if isdefined(iterator.tech, :fuel)
            iterators_sorted[name] = vcat(
                iterators_sorted[name],
                (
                    NamedTuple{(:primemover, :fuel)},
                    ((iterator.tech.primemover), (iterator.tech.fuel)),
                ),
            )
        else
            iterators_sorted[name] = vcat(
                iterators_sorted[name],
                (NamedTuple{(:primemover, :fuel)}, (iterator.tech.primemover, nothing)),
            )
        end
    end
    return iterators_sorted
end

"""
    generators = make_fuel_dictionary(system::PSY.System, results::IS.Results)

This function makes a dictionary of fuel type and the generators associated.

# Arguments
- `c_sys5_re::PSY.System`: the system that is used to create the results
- `results::IS.Results`: results

# Key Words
- `categories::Dict{String, NamedTuple}`: if stacking by a different category is desired

# Example
results = solve_op_model!(OpModel)
generators = make_fuel_dictionary(c_sys5_re, results)

"""
function make_fuel_dictionary(sys::PSY.System, res::IS.Results; kwargs...)

    categories = Dict()
    categories["Solar"] = NamedTuple{(:primemover, :fuel)}, (PSY.PVe, nothing)
    categories["Wind"] = NamedTuple{(:primemover, :fuel)}, (PSY.WT, nothing)
    categories["Oil_CT"] =
        NamedTuple{(:primemover, :fuel)}, (PSY.CT, PSY.DISTILLATE_FUEL_OIL)
    categories["Oil_ST"] =
        NamedTuple{(:primemover, :fuel)}, (PSY.ST, PSY.DISTILLATE_FUEL_OIL)
    categories["Storage"] = NamedTuple{(:primemover, :fuel)}, (PSY.BA, nothing)
    categories["Gas_CT"] = NamedTuple{(:primemover, :fuel)}, (PSY.CT, PSY.NATURAL_GAS)
    categories["Gas_CC"] = NamedTuple{(:primemover, :fuel)}, (PSY.CC, PSY.NATURAL_GAS)
    categories["Hydro"] = NamedTuple{(:primemover, :fuel)}, (PSY.HY, nothing)
    categories["Coal"] = NamedTuple{(:primemover, :fuel)}, (PSY.ST, PSY.COAL)
    categories["Nuclear"] = NamedTuple{(:primemover, :fuel)}, (PSY.ST, PSY.NUCLEAR)
    categories = get(kwargs, :categories, categories)
    iterators = _get_iterator(sys, res)
    generators = Dict()

    for (category, fuel_type) in categories
        generators["$category"] = []
        for (name, fuels) in iterators
            for fuel in fuels
                if fuel == fuel_type
                    generators["$category"] = vcat(generators["$category"], name)
                end
            end
        end
        if isempty(generators["$category"])
            delete!(generators, "$category")
        end
    end
    return generators
end

function _aggregate_data(res::IS.Results, generators::Dict)
    All_var = DataFrames.DataFrame()
    var_names = collect(keys(res.variables))
    for i in 1:length(var_names)
        All_var = hcat(All_var, res.variables[var_names[i]], makeunique = true)
    end
    fuel_dataframes = Dict()

    for (k, v) in generators
        generator_df = DataFrames.DataFrame()
        for l in v
            generator_df = hcat(generator_df, All_var[:, Symbol("$(l)")], makeunique = true)
        end
        fuel_dataframes[k] = generator_df
    end

    return fuel_dataframes
end
"""
    stack = get_stacked_aggregation_data(res::IS.Results, generators::Dict)

This function aggregates the data into a struct type StackedGeneration
so that the results can be plotted using the StackedGeneration recipe.

# Example
```julia
using Plots
gr()
generators = make_fuel_dictionary(res, system)
stack = get_stacked_aggregation_data(res, generators)
plot(stack)
```
*OR*
```julia
using Plots
gr()
fuel_plot(res, system)
```
"""
function get_stacked_aggregation_data(res::IS.Results, generators::Dict)
    # order at the top
    category_aggs = _aggregate_data(res, generators)
    time_range = res.time_stamp[!, :Range]
    labels = collect(keys(category_aggs))
    new_labels = []

    for fuel in order
        for label in labels
            if label == fuel
                new_labels = vcat(new_labels, label)
            end
        end
    end
    variable = category_aggs[(new_labels[1])]
    data_matrix = sum(Matrix(variable), dims = 2)
    legend = [string.(new_labels)[1]]
    for i in 2:length(new_labels)
        variable = category_aggs[(new_labels[i])]
        legend = hcat(legend, string.(new_labels)[i])
        data_matrix = hcat(data_matrix, sum(Matrix(variable), dims = 2))
    end
    return StackedGeneration(time_range, data_matrix, legend)
end
"""
    bar = get_bar_aggregation_data(results::IS.Results, generators::Dict)

This function aggregates the data into a struct type StackedGeneration
so that the results can be plotted using the StackedGeneration recipe.

# Example
```julia
using Plots
gr()
generators = make_fuel_dictionary(res, system)
bar = get_bar_aggregation_data(res, generators)
plot(bar)
```
*OR*
```julia
using Plots
gr()
fuel_plot(res, system)
```
"""
function get_bar_aggregation_data(res::IS.Results, generators::Dict)
    category_aggs = _aggregate_data(res, generators)
    time_range = res.time_stamp[!, :Range]
    labels = collect(keys(category_aggs))
    new_labels = []
    for fuel in order
        for label in labels
            if label == fuel
                new_labels = vcat(new_labels, label)
            end
        end
    end
    variable = category_aggs[(new_labels[1])]
    data_matrix = sum(Matrix(variable), dims = 2)
    legend = [string.(new_labels)[1]]
    for i in 2:length(new_labels)
        variable = category_aggs[(new_labels[i])]
        data_matrix = hcat(data_matrix, sum(Matrix(variable), dims = 2))
        legend = hcat(legend, string.(new_labels)[i])
    end
    bar_data = sum(data_matrix, dims = 1)
    return BarGeneration(time_range, bar_data, legend)
end


text:
include("infdecomp.jl")
using InfDecomp

list_of_functions = [
                     :rdn,
                     :unq,
                     :xor,
                     :and,
                     :rdnxor,
                     :rdnunqxor,
                     :xorand,
#                     :discr_normal
                     ]

function sample(fun::Symbol, Pr_noise::Float64, sz)
    if fun==:rdn # X=Y=Z
        @assert sz ≥ 2  "Need size ≥ 2 for rdn"
        y=rand(0:sz-1)
        z=y
        x=y
        if rand() ≤ Pr_noise
            x = mod(x + rand([+1,-1])  ,sz)
        end
        return (x,y,z)
    elseif fun==:unq # (Y,Z),Y,Z
        @assert sz ≥ 2  "Need size ≥ 2 for unq"
        y=rand(0:sz-1)
        z=rand(0:sz-1)
        x=(y,z)
        if rand() ≤ Pr_noise
            while x==(y,z)
                x = (  mod(y+rand([+1,0,-1]),sz), mod(z+rand([+1,0,-1]),sz)  )
            end
        end
        return ( x,y,z )
    elseif fun==:xor # (Y xor Z),Y,Z
        @assert ispow2(sz)  "Need size to be a power of 2 for xor"
        y=rand(0:sz-1)
        z=rand(0:sz-1)
        x= y $ z
        if rand() ≤ Pr_noise
            l = Int(round(log2(sz)))
            b = rand(0:l-1)
            N = 2^b
            x $= N
        end
        return ( x,y,z )
    elseif fun==:and # (Y and Z),Y,Z
        @assert ispow2(sz)  "Need size to be a power of 2 for and"
        y=rand(0:sz-1)
        z=rand(0:sz-1)
        x= y & z
        if rand() ≤ Pr_noise
            l = Int(round(log2(sz)))
            b = rand(0:l-1)
            N = 2^b
            x $= N
        end
        return ( x,y,z )
    elseif fun==:rdnxor # (  (U xor V,W),  (U,W), (V,W) )
        @assert ispow2(sz)  "Need size to be a power of 2 for rdnxor"
        u=rand(0:sz-1)
        v=rand(0:sz-1)
        w=rand(0:sz-1)
        y = (u,w)
        z = (u,w)
        x = ( u $ v , w)
        if rand() ≤ Pr_noise
            l = Int(round(log2(sz)))
            b = rand(0:l-1)
            N = 2^b
            r = rand()
            if r ≤ 2/5
                x = ( u $ v $ N, w)
            elseif r ≤ 4/5
                x = ( u $ v, w $ N)
            else
                b = rand(0:l-1)
                M = 2^b
                x = ( u $ v $ M, w $ N)
            end
        end
        return ( x,y,z )
    elseif fun==:xorand
        @assert ispow2(sz)  "Need size to be a power of 2 for and"
        y=rand(0:sz-1)
        z=rand(0:sz-1)
        x= (y$z, y&z)
        if rand() ≤ Pr_noise
            l = Int(round(log2(sz)))
            b = rand(0:l-1)
            N = 2^b
            if rand() ≤ 1/2
                x = (y$z $ N, y&z)
            else
                x = (y$z, y&z $ N)
            end
        end
        return ( x,y,z )
    elseif fun==:rdnunqxor # (
        @assert typeof(sz) <: Tuple  "Need size to be tuple for rdnunqxor"
        @assert length(sz) == 2      "Need size to be pair for rdnunqxor"
        @assert ispow2(sz[1])        "Need sizes to be powers of 2 for rdnunqxor"
        @assert sz[1]>=4             "Need 1st size to be ≥ 2 for rdnunqxor"
        l  = Int(round(log2(sz[1])))
        l₁ = Int(floor(l/2))
        l₂ = Int(ceil(l/2))
        s₁ = 2^l₁
        s₂ = 2^l₂
        u₁ = rand(0:s₁-1)
        u₂ = rand(0:s₂-1)
        v₁ = rand(0:s₁-1)
        v₂ = rand(0:s₂-1)
        w  = rand(0:sz[2]-1)
        x = (u₁$v₁,u₂,v₂,w)
        y = (u₁,u₂,w)
        z = (v₁,v₂,w)
        if rand() ≤ Pr_noise
            r = rand()
            if r ≤ 1/5
                u₁ = rand(0:s₁-1)
            elseif r ≤ 2/5
                u₂ = rand(0:s₂-1)
            elseif r ≤ 3/5
                v₁ = rand(0:s₁-1)
            elseif r ≤ 4/5
                v₂ = rand(0:s₂-1)
            else
                w  = rand(0:sz[2]-1)
            end
            x = (u₁$v₁,u₂,v₂,w)
        end
        return ( x,y,z )
    # elseif fun==:discr_normal
    #     blah
    else
        error("sample(): Unknown function, $fun")
    end
    ;
end

function make_PDF(fun::Symbol, n_samples::Int, Pr_noise::Float64, sz, FLOAT::DataType=Float64)
    t = sample(fun,0.,sz)
    T = typeof(t)
    @assert n_samples >= 1
    @assert T <: Tuple "Something went wrong, I didn't get a tuple!"
    @assert length(t) == 3 "Something went wrong, I didn't get a 3-tuple!"

    pdf =  Dict{T,FLOAT}()

    for iter = 1:n_samples
        t = sample(fun,Pr_noise,sz)
        pdf[t] = get!(pdf,t,FLOAT(0)) + FLOAT(1)
    end

    for t in keys(pdf)
        pdf[t] /= n_samples
    end

    return pdf
    ;
end # make_PDF()

function test()
    for f in list_of_functions
        if f ≠ :rdnunqxor
            q = make_PDF(f,100,.05,16)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,8)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,4)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,2)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,16)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,8)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,4)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,2)
            d,e = InfDecomp.create_stuff(q)
        else
            q = make_PDF(f,100,.05,(16,3))
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,(8,3))
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,(4,3))
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,(16,5))
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,(8,5))
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,(4,5))
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(16,3))
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(8,3))
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(4,3))
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(16,5))
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(8,5))
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(4,5))
            d,e = InfDecomp.create_stuff(q)
        end
        if f ≠ :rdnunqxor
            q = make_PDF(f,100,.05,16,BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,8,BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,4,BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,2,BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,16,BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,8,BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,4,BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,2,BigFloat)
            d,e = InfDecomp.create_stuff(q)
        else
            q = make_PDF(f,100,.05,(16,3),BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,(8,3),BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,(4,3),BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,(16,5),BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,(8,5),BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.05,(4,5),BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(16,3),BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(8,3),BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(4,3),BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(16,5),BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(8,5),BigFloat)
            d,e = InfDecomp.create_stuff(q)
            q = make_PDF(f,100,.95,(4,5),BigFloat)
            d,e = InfDecomp.create_stuff(q)
        end
    end
end

; # EOF

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
module MathTeXEngine
# Adapted from matplotlib mathtext engine
# https://github.com/matplotlib/matplotlib/blob/master/lib/matplotlib/_mathtext.py

using AbstractTrees
using Automa
using FreeTypeAbstraction
using LaTeXStrings

using Automa.RegExp: @re_str
using DataStructures: Stack
using GeometryBasics: Point2f, Rect2f
using REPL.REPLCompletions: latex_symbols
using RelocatableFolders

import FreeTypeAbstraction:
    ascender, boundingbox, descender, get_extent, hadvance, inkheight, inkwidth,
    height_insensitive_boundingbox, leftinkbound, rightinkbound,
    topinkbound, bottominkbound

const re = Automa.RegExp

export TeXExpr, texparse
export TeXElement, TeXChar, VLine, HLine, generate_tex_elements
export get_font, get_fontpath

# Reexport from LaTeXStrings
export @L_str

include("parser/texexpr.jl")
include("parser/commands_data.jl")
include("parser/commands_registration.jl")
include("parser/parser.jl")

include("engine/computer_modern_data.jl")
include("engine/fonts.jl")
include("engine/layout_context.jl")
include("engine/texelements.jl")
include("engine/layout.jl")

end # module


text:
using DifferentialEquations
using LSODA
using JLD2
using ArgParse

parser = ArgParseSettings(allow_ambiguous_opts=false)
@add_arg_table! parser begin
    "output"
        help = "output JLD2"
        required = true
    "form"
        help = "noise form: add or mul"
        required = true
    "initstate"
        help = "initial state: hi or lo"
        required = true
    "signals"
        help = "signal values"
        required = true
    "--cells"
        help = "number of cells to simulate"
        default = 500
        arg_type = Int
    "--noises"
        help = "noise values"
        default = "0.5:0.5:1.5"
    "--stochtime"
        help = "stochastic simulation endpoint"
        arg_type = Float64
        default = 200.0
end
args = parse_args(parser)

function model_geneosc(du, u, p, t)
    aM, aP, aF, bM, bP, bF, kf, kb, σM, σP, σF, τ = p

    M, P, F, DA, DR = u[1], u[2], u[3], u[4], u[5]

    du[1] = aM * DA - bM * M + u[6]
    du[2] = aP * M - bP * P + u[7]
    du[3] = aF * P - bF * F - kf * F * DA + kb * DR + u[8]
    du[5] = kf * F * DA - kb * DR - bF * DR
    du[4] = -du[5]

    du[6:8] = -u[6:8]/τ
end

function stoch_add(du,u,p,t)
    aM, aP, aF, bM, bP, bF, kf, kb, σM, σP, σF, τ = p
    du[1:5] .= 0
    du[6] = sqrt(2.0*σM^2.0/τ)
    du[7] = sqrt(2.0*σP^2.0/τ)
    du[8] = sqrt(2.0*σF^2.0/τ)
end

function stoch_mul(du,u,p,t)
    aM, aP, aF, bM, bP, bF, kf, kb, σM, σP, σF, τ = p
    du[1:5] .= 0
    du[6] = sqrt(2.0*σM^2.0/τ) * u[1]
    du[7] = sqrt(2.0*σP^2.0/τ) * u[2]
    du[8] = sqrt(2.0*σF^2.0/τ) * u[3]
end

p_kim_add = [
    15.1745, # aM
    1, # aP
    1, # aF
    1, # bM
    1, # bP
    1, # bF
    200, # kf
    50, # kb
    20, # σM
    20, # σP
    2, # σF
    0.1, # τ
]
p_kim_mul = copy(p_kim_add)
p_kim_mul[9:11] .= 1. # σ parameters

function filter_sols(sols)
    nvars, ntps, nsols = size(sols)
    good_sols = []
    for i in 1:nsols
        if isassigned(sols[i][1,:], ntps) == true
            push!(good_sols, i)
        end
    end
    println("Obtained $(length(good_sols))/$nsols converged solutions.")
    sols[:, :, good_sols]
end

function sim_geneosc_ic(model, model_stoch, u0, pars; popsize=1000, selectR::String="lo", saveat=0.1)
    selectR == "lo" || selectR == "hi" || error("invalid selectR")
    prob_ode = ODEProblem(model, u0, (0., 200.), pars);
    prob_sde = SDEProblem(model, model_stoch, u0, tspan, pars);
    range = LinRange(0, 500, 20)
    init_conds = [(m, p) for m in range for p in range]
    function prob_func_cells_init(prob, i, repeat)
        m, p = init_conds[i]
        remake(prob, u0=[m; p; u0[3:end]])
    end
    prob_cells_init = EnsembleProblem(prob_ode, prob_func=prob_func_cells_init)
    @time sim_init = solve(prob_cells_init, lsoda(), saveat=saveat, dt=0.01, reltol=1e-6, abstol=1e-6, trajectories=length(init_conds));
    sol_init = sim_init[:, 400:end, :]
    if selectR == "lo"
        sol_init[isnan.(sol_init)] .= Inf
        idsel = argmin(sol_init[1, :, :])
    else
        sol_init[isnan.(sol_init)] .= -Inf
        idsel = argmax(sol_init[1, :, :])
    end
    u0_s = sol_init[:, idsel[1], idsel[2]]
    println("Starting stochastic simulation with init cond ", u0_s)
    function prob_func_cells(prob, i, repeat)
        remake(prob, u0=u0_s)
    end
    prob_cells = EnsembleProblem(prob_sde, prob_func=prob_func_cells)
    @time sim = solve(prob_cells, STrapezoid(), saveat=saveat, maxiter=1e12, dt=0.0001, trajectories=popsize);
    filter_sols(sim)
end

stochs = Dict("add" => stoch_add, "mul" => stoch_mul)
psets = Dict("add" => p_kim_add, "mul" => p_kim_mul)

stoch = stochs[args["form"]]
u0 = [0., 0., 0., 164.75, 0., 0., 0., 0.]
p_basal = psets[args["form"]]
tspan = (0., args["stochtime"])
noises = eval(Meta.parse(args["noises"]))
signals = eval(Meta.parse(args["signals"]))

all_results = map(noises) do noise
    results = map(signals) do signal
        p = copy(p_basal)
        p[1] = signal
        p[9:11] *= noise
        signal => sim_geneosc_ic(model_geneosc, stoch, u0, p; popsize=args["cells"], selectR=args["initstate"])
    end
    noise => Dict(results)
end
save(args["output"], "results", Dict(all_results))


text:
"""
    Sphere{T,N} <: ParametricSurface{T,N}

Spherical surface centered at the origin.

```julia
Sphere(radius::T = 1.0; interface::NullOrFresnel{T} = nullinterface(T))
```
"""
struct Sphere{T,N} <: ParametricSurface{T,N}
    radius::T
    interface::NullOrFresnel{T}

    function Sphere(radius::T; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real}
        @assert !isnan(radius)
        @assert radius > zero(T)
        return new{T,3}(radius, interface)
    end
end
export Sphere

interface(a::Sphere{T,N}) where {T<:Real,N} = a.interface
radius(a::Sphere{T}) where {T<:Real} = a.radius

uvrange(::Type{Sphere{T,N}}) where {T<:Real,N} = ((-T(π), T(π)), (zero(T), T(π)))

onsurface(sph::Sphere{T,3}, x::T, y::T, z::T) where {T<:Real} = samepoint(x^2 + y^2 + z^2, radius(sph)^2)

inside(sph::Sphere{T,3}, x::T, y::T, z::T) where {T<:Real} = x^2 + y^2 + z^2 - radius(sph)^2 < zero(T)

point(sph::Sphere{T,3}, ϕ::T, θ::T) where {T<:Real} = SVector{3,T}(radius(sph) * sin(θ) * cos(-ϕ), radius(sph) * sin(θ) * sin(-ϕ), radius(sph) * cos(θ))

normal(::Sphere{T,3}, ϕ::T, θ::T) where {T<:Real} = SVector{3,T}(sin(θ) * cos(-ϕ), sin(θ) * sin(-ϕ), cos(θ))

partials(sph::Sphere{T,3}, ϕ::T, θ::T) where {T<:Real} = SVector{3,T}(radius(sph) * sin(θ) * -sin(-ϕ), radius(sph) * sin(θ) * cos(-ϕ), 0.0), SVector{3,T}(radius(sph) * cos(θ) * cos(-ϕ), radius(sph) * cos(θ) * sin(-ϕ), radius(sph) * -sin(θ))
function uv(::Sphere{T,3}, x::T, y::T, z::T) where {T<:Real}
    # avoid divide by zero for ForwardDiff
    ϕ = -NaNsafeatan(y, x)
    if x == zero(T) && y == zero(T)
        θ = zero(T)
    else
        θ = NaNsafeatan(sqrt(x^2 + y^2), z)
    end
    return SVector{2,T}(ϕ, θ)
end

# Assumes the ray has been transformed into the canonical sphere coordinate frame which has the vertical axis passing through (0,0,0) and aligned with the z axis.
function surfaceintersection(sph::Sphere{T,N}, r::AbstractRay{T,N}) where {T<:Real,N}
    ox, oy, oz = origin(r)
    dx, dy, dz = direction(r)
    rad = radius(sph)

    a = dx^2 + dy^2 + dz^2
    b = 2 * (ox * dx + oy * dy + oz * dz)
    c = (ox^2 + oy^2 + oz^2) - rad^2

    temp = quadraticroots(a, b, c)

    if temp === nothing
        return EmptyInterval(T) # no intersection with sphere and ray not contained entirely in sphere
    end

    t1, t2 = temp

    if isapprox(t1, t2, rtol = 1e-12, atol = 2 * eps(T))
        return EmptyInterval(T) # single root which indicates a tangent ray sphere intersection
    end

    if t1 > zero(T)
        pt1 = point(r, t1)
    else
        pt1 = nothing
    end

    if t2 > zero(T)
        pt2 = point(r, t2)
    else
        pt2 = nothing
    end

    let int1 = nothing, int2 = nothing
        if pt1 !== nothing
            θ, ρ = uv(sph, pt1)
            int1 = Intersection(t1, pt1, pt1, θ, ρ, interface(sph))
        end

        if pt2 !== nothing
            θ, ρ = uv(sph, pt2)
            int2 = Intersection(t2, pt2, pt2, θ, ρ, interface(sph))
        end

        if int1 !== nothing && int2 !== nothing
            if t1 <= t2
                return Interval(int1, int2)
            else
                return Interval(int2, int1)
            end
        elseif int1 !== nothing
            return rayorigininterval(int1)
        elseif int2 !== nothing
            return rayorigininterval(int2)
        else
            return EmptyInterval(T)
        end
    end
end

BoundingBox(sph::Sphere{T,3}) where {T<:Real} = BoundingBox(-radius(sph), radius(sph), -radius(sph), radius(sph), -radius(sph), radius(sph))


text:
@doc raw"""
    UnweightedParticleEnsemble{T}

An ensemble of `N` particles, each of dimension `n`.
"""
struct MultivariateGaussian{S, TM, TP} <: ParametricRepresentation{Vector{S}, S}
    mean::TM
    cov::TP
    MultivariateGaussian(mean::AbstractVector{S}, cov::AbstractMatrix{S}) where S = new{S, typeof(mean), typeof(cov)}(mean, cov)
end

mean(gauss::MultivariateGaussian) = gauss.mean
cov(gauss::MultivariateGaussian)  = gauss.cov
var(gauss::MultivariateGaussian)  = LinearAlgebra.diag(gauss.cov)

parameter_vector(gauss::MultivariateGaussian) = vcat(mean(gauss), vec(cov(gauss)))

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
@testset "Misc tests: show, recipes, inferrability" begin

    # test Base.show
    rd = RefElemData(Tri(),N=3)
    md = MeshData(uniform_mesh(Tri(),1)...,rd)
    @test (@capture_out Base.show(stdout,MIME"text/plain"(),rd)) == "RefElemData for a degree 3 Polynomial() approximation on Tri() element."
    @test (@capture_out Base.show(stdout,rd)) == "RefElemData{N=3,Polynomial(),Tri()}."
    @test (@capture_out Base.show(stdout,MIME"text/plain"(),md)) == "MeshData of dimension 2 with 2 elements"
    @test (@capture_out Base.show(stdout,md)) == "MeshData{2}"

    # test recipes
    # see https://discourse.julialang.org/t/how-to-test-plot-recipes/2648/6?u=jlchan
    meshIO = triangulate_domain(Scramjet())
    recipe = RecipesBase.apply_recipe(Dict{Symbol, Any}(),BoundaryTagPlotter(meshIO));
    @test getfield(recipe[1],1)[:label]=="1"
    @test any(isnan.(getfield(recipe[1],2)[1]))

    recipe = RecipesBase.apply_recipe(Dict{Symbol,Any}(),VertexMeshPlotter(meshIO))
    @test getfield(recipe[1],1)[:legend] == false
    @test getfield(recipe[1],1)[:aspect_ratio] == 1
    @test getfield(recipe[1],1)[:linecolor] == :black
    @test any(isnan.(getfield(recipe[1],2)[1]))

    rd = RefElemData(Tri(),2)
    md = MeshData(triangulateIO_to_VXYEToV(meshIO)...,rd)
    recipe = RecipesBase.apply_recipe(Dict{Symbol,Any}(),MeshPlotter(rd,md))
    @test getfield(recipe[1],1)[:legend] == false
    @test getfield(recipe[1],1)[:aspect_ratio] == 1
    @test getfield(recipe[1],1)[:linecolor] == :black
    @test any(isnan.(getfield(recipe[1],2)[1]))

    rd = RefElemData(Hex(),N=1)
    md = MeshData(uniform_mesh(Hex(),1)...,rd)

    # wrap inferrability test in function
    # https://discourse.julialang.org/t/unexpected-type-instability-with-getproperty-but-not-setproperty/26975/15?u=jlchan
    function foo(rd::RefElemData)
        rd.Nfq,rd.Np,rd.Nq
    end
    function foo(md::MeshData)
        md.VX,md.VY,md.VZ,md.num_elements,md.K
    end
    @inferred foo(rd)
    @inferred foo(md)
end


text:
using Flux, CuArrays, Metalhead, Images
using Metalhead: trainimgs
using Flux: onehotbatch, argmax, crossentropy, throttle, @epochs
using Base.Iterators: partition
using BSON: @save

struct Bottleneck
  layer
end

Flux.treelike(Bottleneck)

Bottleneck(in_planes, growth_rate) = Bottleneck(
                                           Chain(BatchNorm(in_planes, relu),
                                           Conv((1, 1), in_planes=>4growth_rate),
                                           BatchNorm(4growth_rate, relu),
                                           Conv((3, 3), 4growth_rate=>growth_rate, pad = (1, 1))))

(b::Bottleneck)(x) = cat(3, b.layer(x), x)

Transition(chs::Pair{<:Int, <:Int}) = Chain(BatchNorm(chs[1], relu),
                                            Conv((1, 1), chs),
                                            x -> meanpool(x, (2, 2)))

function _make_dense_layers(block, in_planes, growth_rate, nblock)
  local layers = []
  for i in 1:nblock
    push!(layers, block(in_planes, growth_rate))
    in_planes += growth_rate
  end
  Chain(layers...)
end

function DenseNet(nblocks; block = Bottleneck, growth_rate = 12, reduction = 0.5, num_classes = 10)
  num_planes = 2growth_rate
  layers = []
  push!(layers, Conv((3, 3), 3=>num_planes, pad = (1, 1)))
  push!(layers, BatchNorm(num_planes, relu))

  for i in 1:3
    push!(layers, _make_dense_layers(block, num_planes, growth_rate, nblocks[i]))
    num_planes += nblocks[i] * growth_rate
    out_planes = Int(floor(num_planes * reduction))
    push!(layers, Transition(num_planes=>out_planes))
    num_planes = out_planes
  end

  push!(layers, _make_dense_layers(block, num_planes, growth_rate, nblocks[4]))
  num_planes += nblocks[4] * growth_rate
  push!(layers, BatchNorm(num_planes, relu))

  Chain(layers..., x -> meanpool(x, (4, 4)),
        x -> reshape(x, :, size(x, 4)),
        Dense(num_planes, num_classes), softmax)
end

densenet = DenseNet([6, 12, 24, 16]) |> gpu

im_mean = reshape([0.5071, 0.4867, 0.4408], 1, 1, 3)
im_std = reshape([0.2675, 0.2565, 0.2761], 1, 1, 3)

im_to_arr(x) = (permutedims(float.(channelview(x)), [3, 2, 1]) .- im_mean) ./ im_std

X = trainimgs(CIFAR10)
imgs = [im_to_arr(X[i].img) for i in 1:50000]
labels = onehotbatch([X[i].ground_truth.class for i in 1:50000],1:10)
train = [(cat(4, imgs[i]...), labels[:,i]) for i in partition(1:50000, 32)]

loss(x, y) = crossentropy(x, y)

accuracy(x, y) = mean(argmax(x, 1:10) .== argmax(y, 1:10))

opt = ADAM(params(densenet))

@epochs 5 begin
  for data in train
    forward = (densenet(data[1] |> gpu), data[2] |> gpu)
    l = loss(forward[1], forward[2])
    println("Loss = $l and Accuracy = $(accuracy(forward[1], forward[2]))")
    Flux.back!(l)
    opt()
  end
  densenet = densenet |> cpu
  @save "cifar10_checkpoint.bson" densenet
  densenet = densenet |> gpu
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "casacore"
version = v"3.4.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/casacore/casacore.git", "a735bf5f31ea012ac4f7f7378a1f89c1fc136a06"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/casacore

# Apply upstream patch that fixes i386 builds
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/0001-Fix-compilation-of-RefTable-on-i386-platform-1129.patch

# Apply strerror patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/0002-Use-overloading-to-detect-strerror-variant.patch

mkdir build && cd build
CMAKE_FLAGS=()

# Release build for best performance
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=Release)

# Install things into $prefix
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})

# Set the toolchain
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})

# Disable OpenMP, it'll probably interfere with Julia's threads
CMAKE_FLAGS+=(-DUSE_OPENMP=OFF)

# Don't build python hooks
CMAKE_FLAGS+=(-DBUILD_PYTHON=no)

# Stock FFTW3 has threads already baked in
CMAKE_FLAGS+=(-DFFTW3_DISABLE_THREADS=yes)

# linux builds need to link against rt (they don't do this for some reason)
if [[ ${target} == *-linux-* ]]; then
    CMAKE_FLAGS+=(-DCASACORE_ARCH_LIBS="-lrt")
fi

# Congire
cmake ${CMAKE_FLAGS[@]} ..

# Make and install
make -j${nproc}
make install
exit
"""

# Exclude windows, casacore needs pread and pwrite, which are POSIX-only
# Also exclude FreeBSD, lots of upstream bugs that seem not worth fixing ourself
platforms = supported_platforms(exclude=(platform)-> Sys.iswindows(platform) || Sys.isfreebsd(platform))
# Deal with the fact that we have std::string values, which causes issues across the gcc 4/5 boundary
platforms = expand_cxxstring_abis(platforms)
# expand gfortran versions as well
platforms = expand_gfortran_versions(platforms)
# Exclude all musl builds, upstream doesn't care and they use a few glibc-specific features
filter!(p -> libc(p) != "musl", platforms)
# Exclude powerpc64le-linux-gnu-libgfortran3, something is very wonky with that particular combination
# Is it worth fixing?
filter!(p -> !(arch(p) == "powerpc64le" && libgfortran_version(p) == v"3"), platforms)

# The products that we will ensure are always built
products = Product[
    LibraryProduct("libcasa_casa", :libcasa_casa),
    LibraryProduct("libcasa_coordinates", :libcasa_coordinates),
    LibraryProduct("libcasa_derivedmscal", :libcasa_derivedmscal),
    LibraryProduct("libcasa_fits", :libcasa_fits),
    LibraryProduct("libcasa_images", :libcasa_images),
    LibraryProduct("libcasa_lattices", :libcasa_lattices),
    LibraryProduct("libcasa_meas", :libcasa_meas),
    LibraryProduct("libcasa_measures", :libcasa_measures),
    LibraryProduct("libcasa_mirlib", :libcasa_mirlib),
    LibraryProduct("libcasa_msfits", :libcasa_msfits),
    LibraryProduct("libcasa_ms", :libcasa_ms),
    LibraryProduct("libcasa_scimath_f", :libcasa_scimath_f),
    LibraryProduct("libcasa_scimath", :libcasa_scimath),
    LibraryProduct("libcasa_tables", :libcasa_tables),
    ExecutableProduct("findmeastable", :findmeastable),
    ExecutableProduct("fits2table", :fits2table),
    ExecutableProduct("image2fits", :image2fits),
    ExecutableProduct("imagecalc", :imagecalc),
    ExecutableProduct("imageregrid", :imageregrid),
    ExecutableProduct("imageslice", :imageslice),
    ExecutableProduct("lsmf", :lsmf),
    ExecutableProduct("measuresdata", :measuresdata),
    ExecutableProduct("measuresdata-update", :measuresdata_update),
    ExecutableProduct("ms2uvfits", :ms2uvfits),
    ExecutableProduct("msselect", :msselect),
    ExecutableProduct("readms", :readms),
    ExecutableProduct("showtableinfo", :showtableinfo),
    ExecutableProduct("showtablelock", :showtablelock),
    ExecutableProduct("tablefromascii", :tablefromascii),
    ExecutableProduct("taql", :taql),
    ExecutableProduct("tomf", :tomf),
    ExecutableProduct("writems", :writems),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="FFTW_jll", uuid="f5851436-0d7a-5f13-b9de-f02708fd171a"); compat = "3.3.10")
    Dependency(PackageSpec(name="CFITSIO_jll", uuid="b3e40c51-02ae-5482-8a39-3ace5868dcf4"); compat = "4.0.0")
    Dependency(PackageSpec(name="WCS_jll", uuid="550c8279-ae0e-5d1b-948f-937f2608a23e"); compat = "7.7.0")
    Dependency(PackageSpec(name="Readline_jll", uuid="05236dd9-4125-5232-aa7c-9ec0c9b2c25a"); compat = "8.1.1")
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.7", preferred_gcc_version = v"5.2")

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


text:
using DiffEqFlux, Flux
""" CTRNN1 implements a continuous time recurrent neural net
    with 1 hidden layer:
        dxdt = -x/τ + Wσ(Ax + Bu + μ) + ν

    Available types:
        * CTRNN1Nominal             dxdt = -x/τ + Wσ(Ax + Bu + μ) + ν
        * CTRNN1Canonical           dxdt = -x/τ + W(σ(Ax + Bu + μ) - σ(μ))
        * CTRNN1Zero                dxdt = -x/τ + Wσ(Ax + Bu)
"""
abstract type CTRNN1 <: Ode end

""" dxdt = -x/τ + Wσ(Ax + Bu + μ) + ν
    CTRNN1Nominal """
struct CTRNN1Nominal{F,FU,F2} <: CTRNN1
    state_dim::Int
    hidden_dim::Int
    input_dim::Int
    paramsum::Int
    index_dict::Dict{String,UnitRange{Int}}
    σ::F
    σ_u::FU
    initial_params::F2
    function CTRNN1Nominal(
        state_dim::Int,
        hidden_dim::Int,
        input_dim::Int;
        σ = relu,
        initW = Flux.glorot_uniform,
        initb = Flux.zeros,
        σ_u = identity,
    )
        # compute indices
        i = 1
        index_dict = Dict{String,UnitRange{Int}}()
        index_dict["logτ"] = i:i
        i += 1
        index_dict["W"] = i:i+state_dim*hidden_dim-1
        i += state_dim * hidden_dim
        index_dict["A"] = i:i+hidden_dim*state_dim-1
        i += hidden_dim * state_dim
        index_dict["B"] = i:i+hidden_dim*input_dim-1
        i += hidden_dim * input_dim
        index_dict["μ"] = i:i+hidden_dim-1
        i += hidden_dim
        index_dict["ν"] = i:i+state_dim-1
        i += state_dim

        # initialize initial parameters
        initial_params() = vcat(
            Flux.zeros(1),
            vec(initW(state_dim, hidden_dim)),
            vec(initW(hidden_dim, state_dim)),
            vec(initW(hidden_dim, input_dim)),
            vec(initb(hidden_dim)),
            vec(initb(state_dim)),
        )

        new{typeof(σ),typeof(σ_u),typeof(initial_params)}(
            state_dim,
            hidden_dim,
            input_dim,
            i - 1,
            index_dict,
            σ,
            σ_u,
            initial_params,
        )
    end
end

function unpack_params(f::CTRNN1Nominal, p)
    logτ = view(p, f.index_dict["logτ"])
    W = reshape(view(p, f.index_dict["W"]), f.state_dim, f.hidden_dim)
    A = reshape(view(p, f.index_dict["A"]), f.hidden_dim, f.state_dim)
    B = reshape(view(p, f.index_dict["B"]), f.hidden_dim, f.input_dim)
    μ = view(p, f.index_dict["μ"])
    ν = view(p, f.index_dict["ν"])
    return logτ, W, A, B, μ, ν
end

function (f::CTRNN1Nominal)(xu, p)
    x, u = extract_inputs(f, xu)
    logτ, W, A, B, μ, ν = unpack_params(f, p)
    return -x .* exp.(-logτ) .+ (W * f.σ.((A * x .+ B * f.σ_u.(u)) .+ μ) .+ ν)
end

""" dxdt = -x/τ + W(σ(Ax + Bu + μ) - σ(μ))
    CTRNN1Canonical """
struct CTRNN1Canonical{F,FU,F2} <: CTRNN1
    state_dim::Int
    hidden_dim::Int
    input_dim::Int
    paramsum::Int
    index_dict::Dict{String,UnitRange{Int}}
    σ::F
    σ_u::FU
    initial_params::F2
    function CTRNN1Canonical(
        state_dim::Int,
        hidden_dim::Int,
        input_dim::Int;
        σ = relu,
        initW = Flux.glorot_uniform,
        initb = Flux.zeros,
        σ_u = identity,
    )
        # compute indices
        i = 1
        index_dict = Dict{String,UnitRange{Int}}()
        index_dict["logτ"] = i:i
        i += 1
        index_dict["W"] = i:i+state_dim*hidden_dim-1
        i += state_dim * hidden_dim
        index_dict["A"] = i:i+hidden_dim*state_dim-1
        i += hidden_dim * state_dim
        index_dict["B"] = i:i+hidden_dim*input_dim-1
        i += hidden_dim * input_dim
        index_dict["μ"] = i:i+hidden_dim-1
        i += hidden_dim

        # initialize initial parameters
        function initial_params()
            W = initW(state_dim, hidden_dim)
            A = initW(hidden_dim, state_dim)
            return vcat(
                Flux.zeros(1),
                vec(W),
                vec(A / maximum(0.5 * eigen(A * W + W' * A').values)),
                vec(initW(hidden_dim, input_dim)),
                vec(initb(hidden_dim)),
            )
        end
        new{typeof(σ),typeof(σ_u),typeof(initial_params)}(
            state_dim,
            hidden_dim,
            input_dim,
            i - 1,
            index_dict,
            σ,
            σ_u,
            initial_params,
        )
    end
end

function unpack_params(f::CTRNN1Canonical, p)
    logτ = view(p, f.index_dict["logτ"])
    W = reshape(view(p, f.index_dict["W"]), f.state_dim, f.hidden_dim)
    A = reshape(view(p, f.index_dict["A"]), f.hidden_dim, f.state_dim)
    B = reshape(view(p, f.index_dict["B"]), f.hidden_dim, f.input_dim)
    μ = view(p, f.index_dict["μ"])
    ν = -W * f.σ.(μ)
    return logτ, W, A, B, μ, ν
end

function (f::CTRNN1Canonical)(xu, p)
    x, u = extract_inputs(f, xu)
    logτ = view(p, f.index_dict["logτ"])
    W = reshape(view(p, f.index_dict["W"]), f.state_dim, f.hidden_dim)
    A = reshape(view(p, f.index_dict["A"]), f.hidden_dim, f.state_dim)
    B = reshape(view(p, f.index_dict["B"]), f.hidden_dim, f.input_dim)
    μ = view(p, f.index_dict["μ"])
    return -x .* exp.(-logτ) .+ W * (f.σ.(A * x .+ B * f.σ_u.(u) .+ μ) .- f.σ.(μ))
end

""" dxdt = -x/τ + W(σ(Ax + Bu + μ) - σ(μ))
    CTRNN1Zero """
struct CTRNN1Zero{F,FU,F2} <: CTRNN1
    state_dim::Int
    hidden_dim::Int
    input_dim::Int
    paramsum::Int
    index_dict::Dict{String,UnitRange{Int}}
    σ::F
    σ_u::FU
    initial_params::F2
    function CTRNN1Zero(
        state_dim::Int,
        hidden_dim::Int,
        input_dim::Int;
        σ = tanh,
        initW = Flux.glorot_uniform,
        initb = Flux.zeros,
        σ_u = identity,
    )
        # compute indices
        i = 1
        index_dict = Dict{String,UnitRange{Int}}()
        index_dict["logτ"] = i:i
        i += 1
        index_dict["W"] = i:i+state_dim*hidden_dim-1
        i += state_dim * hidden_dim
        index_dict["A"] = i:i+hidden_dim*state_dim-1
        i += hidden_dim * state_dim
        index_dict["B"] = i:i+hidden_dim*input_dim-1
        i += hidden_dim * input_dim

        # initialize initial parameters
        initial_params() = vcat(
            Flux.zeros(1),
            vec(initW(state_dim, hidden_dim)),
            vec(initW(hidden_dim, state_dim)),
            vec(initW(hidden_dim, input_dim)),
        )

        new{typeof(σ),typeof(σ_u),typeof(initial_params)}(
            state_dim,
            hidden_dim,
            input_dim,
            i - 1,
            index_dict,
            σ,
            σ_u,
            initial_params,
        )
    end
end

function unpack_params(f::CTRNN1Zero, p)
    logτ = view(p, f.index_dict["logτ"])
    W = reshape(view(p, f.index_dict["W"]), f.state_dim, f.hidden_dim)
    A = reshape(view(p, f.index_dict["A"]), f.hidden_dim, f.state_dim)
    B = reshape(view(p, f.index_dict["B"]), f.hidden_dim, f.input_dim)
    μ = zeros(f.hidden_dim)
    ν = zeros(f.state_dim)
    return logτ, W, A, B, μ, ν
end

function (f::CTRNN1Zero)(xu, p)
    x, u = extract_inputs(f, xu)
    logτ = view(p, f.index_dict["logτ"])
    W = reshape(view(p, f.index_dict["W"]), f.state_dim, f.hidden_dim)
    A = reshape(view(p, f.index_dict["A"]), f.hidden_dim, f.state_dim)
    B = reshape(view(p, f.index_dict["B"]), f.hidden_dim, f.input_dim)
    return -x .* exp.(-logτ) .+ W * f.σ.(A * x .+ B * f.σ_u.(u))
end

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
using SampleChainsContinuous
using Test

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