text:
# HalidePerovskites
# Uses the newly forked off FeynmanKadanoffOsakaHellwarth.jl package
# Codes by Jarvist Moore Frost, 2017
# These codes were developed with Julia 0.5.0, and requires the Optim and Plots packages.

# This file, when run under Julia, should regenerate all data (and plots) associated with Arxiv paper:
# https://arxiv.org/abs/1704.05404
# Polaron mobility in halide perovskites
# Jarvist Moore Frost
# (Submitted on 18 Apr 2017 [v1])

push!(LOAD_PATH,"../src/") # load module from local directory

using PolaronMobility
using PlotPolaron # Plots dependency

##### load in library routines... #####
#"Rob-MAPCl",  Plot figures with Plots, which defaults to Pyplot backend
using Plots
pyplot()
default(grid=false) # No silly dotted grid lines
default(size=(400,250)) # A good small size for two-column EPS output

#default(size=(800,600)) # Nice size for small-ish PNGs for slides

# Physical constants
const hbar = const ħ = 1.05457162825e-34;          # kg m2 / s
const eV = const q = const ElectronVolt = 1.602176487e-19;                         # kg m2 / s2
const me=MassElectron = 9.10938188e-31;                          # kg
const Boltzmann = const kB =  1.3806504e-23;                  # kg m2 / K s2
const ε_0 = 8.854E-12 #Units: C2N−1m−2, permittivity of free space

" Copy and pasted out of a Jupyter notebook; this calculates 'alpha' parameters
for various materials, as a comparison to the literature used when figuring out
the oft-quoted units. "
function checkalpha()
	println(" Alpha-parameter, Cross check 'feynmanalpha()' fn vs. literature values.\n")
	print("\t NaCl Frohlich paper α=",feynmanalpha(2.3, 5.6, (4.9E13/(2*pi)), 1.0))
    println(" should be ~about 5 (Feynman1955)")
	print("\t CdTe  α=",feynmanalpha(7.1,   10.4,  5.08E12, 0.095))
    println(" Stone 0.39 / Devreese 0.29 ")
	print("\t GaAs  α=",feynmanalpha(10.89, 12.9,  8.46E12, 0.063))
    println(" Devreese 0.068 ")

    println()
    println("Guess at PCBM: 4.0, 6.0 ; α=",feynmanalpha(4.0,6.0, 1E12, 50))
    println("MAPI:")
    println("MAPI  4.5, 24.1, 9THz ; α=",feynmanalpha(4.5,   24.1,  9.0E12,    0.12))
    println("MAPI  4.5, 24.1, 2.25THz - 75 cm^-1 ; α=",feynmanalpha(4.5,   24.1,  2.25E12,    0.12))
    println("MAPI  6.0, 25.7, 9THz ; α=",feynmanalpha(6.0,   25.7,  9.0E12,    0.12))
    println("MAPI  6.0, 36.0, 9THz ; α=",feynmanalpha(6.0,   36,  9.0E12,    0.12))
    println("MAPI  6.0, 36.0, 1THz ; α=",feynmanalpha(6.0,   36,  1.0E12,    0.12))
end
checkalpha()

#####
# Call simulation

# CsSnX3 X={Cl,Br,I}
# L. Huang, W.Lambrecht - PRB 88, 165203 (2013)
# Dielectric consts, from TABLE VII
# Effective masses from TABLE VI, mh*
const cm1=2.997e10 # cm-1 to Herz

function CsSn()
    CsSnCl3=polaronmobility(10:10:1000, 4.80, 29.4,  243cm1, 0.140, figures=true) # alpha= 1.386311
    savepolaron("CsSnCl3",CsSnCl3)

    CsSnBr3=polaronmobility(10:10:1000, 5.35, 32.4,  183cm1, 0.082, figures=true) # alpha= 1.094468
    savepolaron("CsSnBr3",CsSnBr3)

    CsSnI3=polaronmobility(10:10:1000, 6.05, 48.2,  152cm1, 0.069, figures=true) # alpha= 1.020355
    savepolaron("CsSnI3",CsSnI3)

end

CsSn() # generate all that lovely tin data

# Ts, βreds, Kμs, Hμs, FHIPμs, vs, ws, ks, Ms, As, Bs, Cs, Fs, Taus
# CsSnBr3.dat:300 0.877382 511.513 356.358 874.499 8.13686 7.31511 12.6976 0.23729 -3.78992 2.15682 0.84436 0.788735 0.185396
# CsSnCl3.dat:300 1.165048 212.034 147.36 272.993 6.50772 5.58522 11.1557 0.357614 -3.5556 2.42611 0.947833 0.181657 0.143969
# CsSnI3.dat:300 0.728755 703.502 487.085 1448.02 9.64022 8.80508 15.4044 0.198691 -4.01438 2.18761 0.859401 0.967371 0.207864

# So for CsSnSi3, a hole-mobility of 487(Hellwarth)-703(Kadanoff). Compares to expt. 400 (transport) -  585 (Hall effect). Ace!

#Ts,Kμs, Hμs, FHIPμs, ks, Ms, As, Bs, Cs, Fs, Taus
#effectivemass=0.12 # the bare-electron band effective-mass.
# --> 0.12 for electrons and 0.15 for holes, in MAPI. See 2014 PRB.
# MAPI  4.5, 24.1, 2.25THz - 75 cm^-1 ; α=
# MAPI  4.5, 24.1, 2.25THz - 75 cm^-1 ; α=
MAPIe=polaronmobility(10:10:1000, 4.5, 24.1, 2.25E12, 0.12)
plotpolaron("MAPI-electron", MAPIe)
savepolaron("MAPI-electron",MAPIe)

MAPIh=polaronmobility(10:10:1000, 4.5, 24.1, 2.25E12, 0.15)
plotpolaron("MAPI-hole", MAPIh)
savepolaron("MAPI-hole",MAPIh)

# PCBM: 4.0, 5.0, 2.25Thz, effective-mass=1.0
#polaronmobility(4.0, 5.0, 2.25E12, 1.00)

# CsPbI3: 3.82, 8.27, 2.83 THz (highest freq mode in cubic)
#    ^-- v. quick Kmesh=3x3x3 PBESol Calc, JMF 2017-04-12. Files: jarvist@titanium:~/phonopy-work/2017-03-Scott-Distortions/1000-Dielectric/3x3x3-kmesh/
# Horribly non-converged!
#                         Kmesh=6x6x6; 7.2/12.1
# Kmesh=9x9x9x, Ediff=10^-9;           6.1/12.0, 2.57 THz
function CsPbI()
    CsPbI=polaronmobility(10:10:1000, 6.1,6.1+12.0, 2.57E12, 0.12)
    plotpolaron("CsPbI3-electron",CsPbI)
    savepolaron("CsPbI3-electron",CsPbI)
end
CsPbI() # generate CsPbI3 data for later plotting

function SendnerCrosscheck()
    const cm1=2.997e10 # cm-1 to Herz

    # Rob's / Sendner's paper - values extracted form IR measures
    # https://doi.org/10.1039%2Fc6mh00275g

    # These data mainly from Private communication.
    # It is not well described in the paper but they
    # used one the Hellwarth effective-mode method to reduce the observed IR
    # oscillators down to a single mode. Rob provided these values by email (as
    # below). The Effetive Masses and dielectric constants are from the paper.
    # Combined, this reproduces their mobilities, and the internal w and
    # v parameters (again, data from email from Rob):
    #  Paper: "For omega we used: MAPbI/Br/Cl = 112.9/149.4/214.0"
    RobMAPI =polaronmobility(10:10:400, 5.0, 33.5, 112.9*cm1, 0.104)
    RobMAPBr=polaronmobility(10:10:400, 4.7, 32.3,149.4*cm1, 0.117)
    RobMAPCl=polaronmobility(10:10:400, 4.0, 29.8, 214.0*cm1, 0.2)

    savepolaron("Rob-MAPI",RobMAPI)
    savepolaron("Rob-MAPBr",RobMAPBr)
    savepolaron("Rob-MAPCl",RobMAPCl)

    plot(RobMAPI.T,RobMAPI.Hμ,label="(Rob's values) MAPI",markersize=2,marker=:uptriangle,ylim=(0,400))
    plot!(RobMAPBr.T,RobMAPBr.Hμ,label="(Rob's values) MAPBr",markersize=2,marker=:diamond)
    plot!(RobMAPCl.T,RobMAPCl.Hμ,label="(Rob's values) MAPCl",markersize=2,marker=:diamond)
    savefig("Rob-comparison.png")
#    savefig("Rob-comparison.eps")
end

SendnerCrosscheck() # cross-check data for Rob

#####
## Expt. data to compare against
# Milot/Herz 2015 Time-Resolved-Microwave-Conductivity mobilities
# Data from table in SI of: DOI: 10.1002/adfm.201502340
# Absolute values possibly dodge due to unknown yield of charge carriers; but hopefully trend A.OK!
Milot= [
8 184
40 321
80 143
120 62
140 40
160 52
180 44
205 41
230 39
265 26
295 35
310 24
320 24
330 19
340 16
355 15
]

# IV estimated mobilities (?!) from large single crystals, assumed ambient T
# Nature Communications 6, Article number: 7586 (2015)
# doi:10.1038/ncomms8586
Saidaminov =
[ 300 67.2 ]

#Semonin2016,
#  doi = {10.1021/acs.jpclett.6b01308},
Semonin =
[ 300 115 ] # +- 15 cm^2/Vs, holes+electrons

#####
## Calculated mobilities vs. expt
plot(Milot[:,1],Milot[:,2],label="Milot T-dep TRMC Polycrystal",marker=:hexagon,markersize=3,
xlab="Temperature (K)",ylab="Mobility (cm\$^2\$/Vs)", ylims=(0,400) )
plot!(Saidaminov[:,1],Saidaminov[:,2],label="",markersize=6,marker=:utriangle)
plot!(Semonin[:,1],Semonin[:,2],label="",markersize=6,marker=:square)
#plot!(Saidaminov[:,1],Saidaminov[:,2],label="Saidaminov JV Single Crystal", markersize=6,marker=:utriangle)
#plot!(Semonin[:,1],Semonin[:,2],label="Semonin Single Crystal TRMC", markersize=6,marker=:hexagon)

#plot!(MAPIe.T,MAPIe.Kμs,label="(electrons) Kadanoff Polaron mobility",marker=2)
#plot!(Ts,hHμs,label="Calculated (hole) mobility",markersize=3,marker=:dtriangle)

plot!(MAPIe.T,MAPIe.Hμ,label="Calculated (electron) mobility",markersize=3,marker=:dtriangle)
plot!(MAPIh.T,MAPIh.Hμ,label="Calculated (hole) mobility",markersize=3,marker=:utriangle)

savefig("MAPI-eh-mobility-calculated-experimental.png")
savefig("MAPI-eh-mobility-calculated-experimental.eps")

plot!(ylims=(),yscale=:log10)
plot!(ylims=(),xscale=:log10) # on log-log axes; powerlaw becomes straight
plot!(x->1e5*x^(-3/2),20:1000,yscale=:log10) # Trend line for the T^-3/2 powerlaw

savefig("MAPI-eh-mobility-calculated-experimental-log10.png")
savefig("MAPI-eh-mobility-calculated-experimental-log10.eps")

# OK; fit some stuff.
using LsqFit # Least squares fit, orig from Optim package.

function fitPowerLaw()
    model(x,p)=p[1]*x.^(p[2])
    p0=[1.0,-3.0/2]

    # Data we're fitting against
    plot(MAPIe.T,MAPIe.Hμ,label="Calculated (electron) mobility",markersize=3,marker=:dtriangle,ylims=(1,400))

    # They're nasty 'Any' data types; force to Float64
    T=convert(Array{Float64},MAPIe.T)
    μ=convert(Array{Float64},MAPIe.Hμ)

    # Ta da! Least squares fit.
    fit=curve_fit(model,T,μ,p0)
    p=fit.param
    #plot!(T,model(T,p))

    # Fit just from 100K onwards
    fit=curve_fit(model,T[10:end],μ[10:end],p0)
    p=fit.param
    plot!(T[10:end],model(T[10:end],p), label="Calc [10:end] Powerlaw: $p")

    # Fit all from 50K onwards - very dubious of quality of fit
    #fit=curve_fit(model,T[5:end],μ[5:end],p0)
    #p=fit.param
    #plot!(T[5:end],model(T[5:end],p), label="Calc [5:end] Powerlaw: $p")

    # TRMC Polycrystal data for comparison
    plot!(Milot[:,1],Milot[:,2],label="Milot T-dep TRMC Polycrystal",marker=:hexagon,markersize=3,
    xlab="Temperature (K)",ylab="Mobility (cm\$^2\$/Vs)", ylims=(0,400) )

    # Fit all data points except for spurious first one
    MT=convert(Array{Float64},Milot[2:end,1])
    Mu=convert(Array{Float64},Milot[2:end,2])
    fit=curve_fit(model,MT,Mu,p)
    p=fit.param
    plot!(T[5:end], model(T[5:end],p), label="Milot Powerlaw Fit 2:end: $p")

    # Fit only sensible data points > 120 K and up
    MT=convert(Array{Float64},Milot[4:end,1])
    Mu=convert(Array{Float64},Milot[4:end,2])
    fit=curve_fit(model,MT,Mu,p)
    p=fit.param
    plot!(T[5:end], model(T[5:end],p), label="Milot Powerlaw Fit 4:end: $p")

    savefig("MAPI-fit-lin.png")
    savefig("MAPI-fit-lin.pdf")

    plot!(ylims=(1,), yscale=:log10)# log-log; power-law --> straight line
    plot!(xlims=(1,), xscale=:log10)
    savefig("MAPI-fit-log.png")
    plot!(legend=false) # no key
    savefig("MAPI-fit-log.pdf")
end

fitPowerLaw()

#plot!(RobMAPI.T,RobMAPI.Hμ,label="(Rob's values) MAPI",markersize=2,marker=:rect)
#savefig("MAPI-eh-mobility-calculated-experimental-Rob.png")
#savefig("MAPI-eh-mobility-calculated-experimental-Rob.eps")

println("That's me!")


text:
using ProgressMeter, Kinetic

begin
    # space
    x0 = -1.5
    x1 = 1.5
    y0 = -1.5
    y1 = 1.5
    nx = 100
    ny = 100
    dx = (x1 - x0) / nx
    dy = (y1 - y0) / ny

    pspace = PSpace2D(x0, x1, nx, y0, y1, ny)

    # time
    tEnd = 1.0
    cfl = 0.95

    # quadrature
    quadratureorder = 6
    points, triangulation = octa_quadrature(quadratureorder)
    weights = Kinetic.quadrature_weights(points, triangulation)
    nq = size(points, 1)

    # particle
    SigmaS = 1 * ones(ny + 4, nx + 4)
    SigmaA = 0 * ones(ny + 4, nx + 4)
    SigmaT = SigmaS + SigmaA
end

# initial distribution
phi = zeros(nq, nx, ny)
s2 = 0.03^2
flr = 1e-4
init_field(x, y) = max(flr, 1.0 / (4.0 * pi * s2) * exp(-(x^2 + y^2) / 4.0 / s2))
for j = 1:nx
    for i = 1:ny
        y = y0 + dy / 2 + (i - 3) * dy
        x = x0 + dx / 2 + (j - 3) * dx
        for q = 1:nq
            phi[q, i, j] = init_field(x, y) / 4.0 / π
        end
    end
end

dt = cfl / 2 * (dx * dy) / (dx + dy)
global t = 0.0

flux1 = zeros(nq, nx + 1, ny)
flux2 = zeros(nq, nx, ny + 1)

@showprogress for iter = 1:50
    for i = 2:nx, j = 1:ny
        tmp = @view flux1[:, i, j]
        flux_kfvs!(tmp, phi[:, i-1, j], phi[:, i, j], points[:, 1], dt)
    end
    for i = 1:nx, j = 2:ny
        tmp = @view flux2[:, i, j]
        flux_kfvs!(tmp, phi[:, i, j-1], phi[:, i, j], points[:, 2], dt)
    end

    for j = 1:ny, i = 1:nx
        integral = discrete_moments(phi[:, i, j], weights)
        integral *= 1.0 / 4.0 / pi

        for q = 1:nq
            phi[q, i, j] =
                phi[q, i, j] +
                (flux1[q, i, j] - flux1[q, i+1, j]) / dx +
                (flux2[q, i, j] - flux2[q, i, j+1]) / dy +
                (integral - phi[q, i, j]) * dt
        end
    end

    global t += dt
end

ρ = zeros(nx, ny)
for i = 1:nx, j = 1:ny
    ρ[i, j] = discrete_moments(phi[:, i, j], weights)
end

using Plots
contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], ρ[:, :])

ext: jl


text:
module Singularity

include("generate.jl")
include("building.jl")
include("reproduce.jl")

using FileIO
using DrWatson

function __init__()
    if !Sys.islinux()
        error("This package currently requires singularity to be installed,
        which is only available for linux")
        # TODO add check for singularity (even further, add singularity as binary)
    end
end

end # module


text:
# ============================================================
# Uncertainty Quantification Methods
# ============================================================

"""
Calculate collocation -> polynomial chaos

"""
function ran_chaos(ran::AbstractArray{<:AbstractFloat,1}, uq::UQ1D)
    chaos = zeros(eltype(ran), uq.nr + 1)
    for j = 1:uq.nr+1
        chaos[j] =
            sum(@. uq.op.quad.weights * ran * uq.phiRan[:, j]) /
            (uq.t2Product[j-1, j-1] + 1.e-7)
    end

    return chaos
end

function ran_chaos(ran::AbstractArray{<:AbstractFloat,1}, uq::UQ2D)
    chaos = zeros(eltype(ran), uq.nm + 1)
    for j in eachindex(chaos)
        chaos[j] =
            sum(@. uq.weights * ran * uq.phiRan[:, j]) / (uq.t2Product[j-1, j-1] + 1.e-7)
    end

    return chaos
end

function ran_chaos(ran::AbstractArray{<:AbstractFloat,1}, op::AbstractOrthoPoly)

    phiRan = evaluate(collect(0:op.deg), op.quad.nodes, op)
    t2 = Tensor(2, op)

    chaos = zeros(eltype(ran), op.deg + 1)
    for j = 1:op.deg+1
        chaos[j] =
            sum(@. op.quad.weights * ran * phiRan[:, j]) / (t2.get([j - 1, j - 1]) + 1.e-7)
    end

    return chaos

end

function ran_chaos(uRan::AbstractArray{<:AbstractFloat,2}, idx::Integer, uq::AbstractUQ)

    if idx == 1

        uChaos = zeros(uq.nr + 1, axes(uRan, 2))
        for j in axes(uChaos, 2)
            uChaos[:, j] .= ran_chaos(uRan[:, j], uq)
        end

    elseif idx == 2

        uChaos = zeros(axes(uRan, 1), uq.nr + 1)
        for i in axes(uChaos, 1)
            uChaos[i, :] .= ran_chaos(uRan[i, :], uq)
        end

    end

    return uChaos

end

function ran_chaos(uRan::AbstractArray{<:AbstractFloat,3}, idx::Integer, uq::AbstractUQ)

    if idx == 1

        uChaos = zeros(uq.nr + 1, axes(uRan, 2), axes(uRan, 3))
        for k in axes(uChaos, 3)
            for j in axes(uChaos, 2)
                uChaos[:, j, k] .= ran_chaos(uRan[:, j, k], uq)
            end
        end

    elseif idx == 2

        uChaos = zeros(axes(uRan, 1), uq.nr + 1, axes(uRan, 3))
        for k in axes(uChaos, 3)
            for i in axes(uChaos, 1)
                uChaos[i, :, k] .= ran_chaos(uRan[i, :, k], uq)
            end
        end

    elseif idx == 3

        uChaos = zeros(uq.nr + 1, axes(uRan, 2), axes(uRan, 3))
        for k in axes(uChaos, 3)
            for j in axes(uChaos, 2)
                uChaos[:, j, k] .= ran_chaos(uRan[:, j, k], uq)
            end
        end

    end

    return uChaos

end

function ran_chaos(uRan::AbstractArray{<:AbstractFloat,4}, idx::Integer, uq::AbstractUQ)

    if idx == 1

        uChaos = zeros(uq.nr + 1, axes(uRan, 2), axes(uRan, 3), axes(uRan, 4))
        for l in axes(uChaos, 4), k in axes(uChaos, 3), j in axes(uChaos, 2)
            uChaos[:, j, k, l] .= ran_chaos(uRan[:, j, k, l], uq)
        end

    elseif idx == 2

        uChaos = zeros(axes(uRan, 1), uq.nr + 1, axes(uRan, 3), axes(uRan, 4))
        for l in axes(uChaos, 4), k in axes(uChaos, 3), i in axes(uChaos, 1)
            uChaos[i, :, k, l] .= ran_chaos(uRan[i, :, k, l], uq)
        end

    elseif idx == 3

        uChaos = zeros(axes(uRan, 1), axes(uRan, 2), uq.nr + 1, axes(uRan, 4))
        for l in axes(uChaos, 4), j in axes(uChaos, 2), i in axes(uChaos, 1)
            uChaos[i, j, :, l] .= ran_chaos(uRan[i, j, :, l], uq)
        end

    elseif idx == 4

        uChaos = zeros(axes(uRan, 1), axes(uRan, 2), axes(uRan, 3), uq.nr + 1)
        for k in axes(uChaos, 3), j in axes(uChaos, 2), i in axes(uChaos, 1)
            uChaos[i, j, k, :] .= ran_chaos(uRan[i, j, k, :], uq)
        end

    end

    return uChaos

end

"""
Calculate polynomial chaos -> collocation

"""
chaos_ran(chaos::AbstractArray{<:AbstractFloat,1}, uq::UQ1D) =
    evaluatePCE(chaos, uq.op.quad.nodes, uq.op)

chaos_ran(chaos::AbstractArray{<:AbstractFloat,1}, uq::UQ2D) =
    evaluatePCE(chaos, uq.points, uq.op)

chaos_ran(chaos::AbstractArray{<:AbstractFloat,1}, op::AbstractOrthoPoly) =
    evaluatePCE(chaos, op.quad.nodes, op)

function chaos_ran(uChaos::AbstractArray{Float64,2}, idx::Integer, uq::AbstractUQ)

    if idx == 1

        uRan = zeros(uq.op.quad.Nquad, axes(uChaos, 2))
        for j in axes(uRan, 2)
            uRan[:, j] .= chaos_ran(uChaos[:, j], uq)
        end

    elseif idx == 2

        uRan = zeros(axes(uChaos, 1), uq.op.quad.Nquad)
        for i in axes(uRan, 1)
            uRan[i, :] .= chaos_ran(uChaos[i, :], uq)
        end

    end

    return uRan

end

function chaos_ran(uChaos::AbstractArray{Float64,3}, idx::Integer, uq::AbstractUQ)

    if idx == 1

        uRan = zeros(uq.op.quad.Nquad, axes(uChaos, 2), axes(uChaos, 3))
        for k in axes(uRan, 3)
            for j in axes(uRan, 2)
                uRan[:, j, k] .= chaos_ran(uChaos[:, j, k], uq)
            end
        end

    elseif idx == 2

        uRan = zeros(axes(uChaos, 1), uq.op.quad.Nquad, axes(uChaos, 3))
        for k in axes(uRan, 3)
            for i in axes(uRan, 1)
                uRan[i, :, k] .= chaos_ran(uChaos[i, :, k], uq)
            end
        end

    elseif idx == 3

        uRan = zeros(axes(uChaos, 1), axes(uChaos, 2), uq.op.quad.Nquad)
        for j in axes(uRan, 2)
            for i in axes(uRan, 1)
                uRan[i, j, :] .= chaos_ran(uChaos[i, j, :], uq)
            end
        end

    end

    return uRan

end

"""
Calculate λ -> T in polynomial chaos

"""
function lambda_tchaos(lambdaChaos::Array{<:AbstractFloat,1}, mass::Real, uq::AbstractUQ)
    lambdaRan = chaos_ran(lambdaChaos, uq)
    TRan = mass ./ lambdaRan
    TChaos = ran_chaos(TRan, uq)

    return TChaos
end

"""
Calculate T -> λ in polynomial chaos

"""
function t_lambdachaos(TChaos::Array{<:AbstractFloat,1}, mass::Real, uq::AbstractUQ)
    TRan = chaos_ran(TChaos, uq)
    lambdaRan = mass ./ TRan
    lambdaChaos = ran_chaos(lambdaRan, uq)

    return lambdaChaos

end


text:
import VideoIO
import Plots
import AbstractPlotting

import Makie
using Images, CoordinateTransformations, OffsetArrays

function main_makie()
    f = VideoIO.opencamera()
    img = read(f)
    tfm = recenter(RotMatrix(pi/2), Images.center(img))
    img=parent(warp(img,tfm))
    scene = Makie.Scene(resolution = size(img))
    makieimg = Makie.image!(
        scene,
        img,
        show_axis = false,
        scale_plot = true
    )[end]
    display(scene)

    while !eof(f)
        # OffsetArrays -> Array
        img=parent(warp(read(f),tfm))
        @info size(img)
        makieimg[1] = img
        sleep(1/f.framerate)
    end
end

function main_plots()
    f = VideoIO.opencamera()
    img = read(f)
    p=Plots.plot(img) |> display
    for i in 1:100
        @info i
        read!(f,img)
        Plots.plot(img, title="$i") |> display
    end
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Distributions
using Random
using Plots
#set parameters for data generation

μ=0
σ=1
#Sample size, number of simulations, number and size of hidden layers
n=100000
numsim=100
nHidden = [6,15,20,15,6]
#list for keeping track of simulations
βlist=[]
βbiaslist=[]
βnnbabysitlist=[]
βnnvanlist=[]
ValErrBS=[]
ValErrVan=[]
misslist=[]
#Neural Net Parameters and commands
maxIter = 10000
include("NeuralNet.jl")
include("sgdfunc.jl")

for j in 1:numsim

#Draw random varaibles
normdist=Normal(μ,σ)
X1=rand(normdist,n)
X2=rand(normdist,n)
ϵ1=rand(Normal(0,1),n)
ϵ2=rand(Normal(0,10),n)
#form dependent variables
X3=X1+X2+ϵ1
y=2*X1-3*X2+4*X3+ϵ2

#Check OLS estimates for all data
X=hcat(X1,X2,X3)
β=inv(X'*X)*X'*y
push!(βlist,β)

#Systematic Removal
indexl= [] #Contains list of indices for nonmissing data
indexm=[] #Contains list of indices for missing data
for i in 1:n
    k=rand(Binomial(1,.5))
    if (y[i]>1 && ϵ2[i]<5) && k==true #removal conditions
        push!(indexm,i)
    else
        push!(indexl,i)
    end
end

#Number of missing
push!(misslist,length(indexm))

#Get Bias β
Xmiss=hcat(X1[indexl],X2[indexl],X3[indexl])
ymiss=y[indexl]
βbias=inv(Xmiss'*Xmiss)*Xmiss'*ymiss
push!(βbiaslist,βbias)

#Create training and validation sets (70% training/30% validation)
numtrain=convert(Int,round(.7*length(indexl)))
xtrain=hcat(X1[indexl[1:numtrain]],X2[indexl[1:numtrain]])
x3train=(X3[indexl[1:numtrain]])
xvalid=hcat(X1[indexl[(numtrain+1):length(indexl)]],X2[indexl[(numtrain+1):length(indexl)]])
x3valid=(X3[indexl[numtrain+1:length(indexl)]])

#Number of training examples
ntrain = size(xtrain,1)
#add bias variables
xtrain = [ones(ntrain,1) xtrain]
#number of features
dtrain = size(xtrain,2)

#get number of total weights
nParams = NeuralNet_nParams(dtrain,nHidden)

#Run neural network to get weights with Vanilla SGD
bestvanval=Inf
wbestvan=randn(nParams)
for α in [.1, .01,.001,.0001,.00001]
w,validvan=VanillaSGD(NeuralNet_backprop,NeuralNet_predict,maxIter,nHidden, nParams,xtrain,x3train,xvalid,x3valid, α ,10)
if validvan < bestvanval
wbestvan=w
bestvanval=validvan
end
end
push!(ValErrVan,bestvanval)
#predict on missing data
X3nnvan=NeuralNet_predict(wbestvan,hcat(ones(length(indexm)),X1[indexm],X2[indexm]),nHidden)
#form new X and y
Xnnvan=vcat(hcat(X1[indexl],X2[indexl],X3[indexl]),hcat(X1[indexm],X2[indexm],X3nnvan))
ynnvan=vcat(y[indexl],y[indexm])

#Get vanilla SGD OLS estimates
βnnvan=inv(Xnnvan'*Xnnvan)*Xnnvan'*ynnvan
push!(βnnvanlist,βnnvan)

#Run neural network to get weights with SGDbabysitter
w, validBB = SGDBabysitter(NeuralNet_backprop, NeuralNet_predict,
maxIter, nHidden, nParams, xtrain, x3train, xvalid, x3valid,500)
push!(ValErrBS,validBB)
#predict on missing data
X3nnbabysit=NeuralNet_predict(w,hcat(ones(length(indexm)),X1[indexm],X2[indexm]),nHidden)
#form new X and y
Xnnbabysit=vcat(hcat(X1[indexl],X2[indexl],X3[indexl]),hcat(X1[indexm],X2[indexm],X3nnbabysit))
ynnbabysit=vcat(y[indexl],y[indexm])

#Get SGDBabysitter OLS estimates
βnnbabysit=inv(Xnnbabysit'*Xnnbabysit)*Xnnbabysit'*ynnbabysit
push!(βnnbabysitlist,βnnbabysit)

end
#Final report of all OLS estimates
print("All Data Coefficients: ", sum(βlist)/numsim, " Coefficients with Missing Data: ", sum(βbiaslist)/numsim," Coefficients with NN Imputation (SGDBabysitter): ",sum(βnnbabysitlist)/numsim," Coefficients with NN Imputation (VanillaSGD): " ,sum(βnnvanlist)/numsim," Average Number of Missing Variables: ", sum(misslist)/numsim, " ")
print("Average Validation Error SGDBabysitter: ", sum(ValErrBS)/numsim, " Average Validation Error VanillaSGD ", sum(ValErrVan)/numsim, " ", "Average Test Error SGDBabysitter: ", sum((X3nnbabysit-X3[indexm]).^2)/numsim, " Average Validation Error VanillaSGD ", sum((X3nnvan-X3[indexm]).^2)/numsim,)


text:
function OneVsAll(X,y;λs=0,Norms=false)
   m=size(X,1)
   n=size(X,2)
   NLabels=size(unique(y),1)

    Θ=zeros(NLabels,n+1)

    for i=1:NLabels
       Cost, grad, avg, st=LogCostGradient(X,(y.==i).*1.0,λ=λs, Norm=Norms)
       θi=zeros(n+1)
       r=optimize(Cost,grad,θi,LBFGS(),inplace = false)
        Θ[i,:]=r.minimizer

    end

    return Θ
end


text:
using SymbolicUtils
using SymbolicUtils: Term
using SpecialFunctions
using Test
import IfElse: ifelse
import IfElse

using SymbolicUtils: showraw, Symbolic

function rand_input(T)
    if T == Bool
        return rand(Bool)
    elseif T <: Integer
        x = rand(-100:100)
        while iszero(x)
            x = rand(-100:100)
        end
        return x
    elseif T == Rational
        return Rational(rand_input(Int), rand(1:50)) # no 0 denominator tests yet!
    elseif T == Real
        return rand_input(rand([Int, Float64, Rational]))
    elseif T == Complex
        return rand_input(Real) + rand_input(Real) * im
    elseif T == Number
        # more real than complex
        return rand_input(rand([Real, Real, Real, Complex{Float64}]))
    elseif !isabstracttype(T)
        return rand(T)
    else
        error("Don't know how to generate input of type $T")
    end
end

rand_input(i::Symbolic{T}) where {T} = rand_input(T)

const num_spec = let
    @syms a b::Real c::Integer d::Float64 e::Rational f

    leaf_funcs = [()->rand_input(Real),
                  ()->rand_input(Complex),
                  ()->rand([a b c d e f]),
                  ()->rand([a b c d e f])]

    binops = SymbolicUtils.diadic
    nopow  = filter(x->x!==(^), binops)
    twoargfns = vcat(nopow, (x,y)->x isa Union{Int, Rational, Complex{<:Rational}} ? x * y : x^y)
    fns = vcat(1 .=> vcat(SymbolicUtils.monadic, [one, zero]),
               2 .=> vcat(twoargfns, fill(+, 5), [-,-], fill(*, 5)),
    3 .=> [+, *])

    (leaves=leaf_funcs, funcs=fns, input=rand_input)
end

const bool_spec = let
    @syms b::Bool x::Real y::Real z::Complex

    bool_leaf_funcs = [()->rand(Bool),
                       ()->rand([b, (x, b) => ((x > 0) | b), (x,)=>(x < 0), (y,z) => y==z, (y, z) => y!=z])]

    fns = vcat(1 .=> [(!), (~)],
               2 .=> [(|), (&), xor],
               3 .=> [ifelse]) # cond will still stay in bool by condtruction

    (leaves=bool_leaf_funcs,
     funcs=fns,
     input=rand_input
    )
end

function gen_rand_expr(inputs;
                       spec=num_spec,
                       leaf_prob=0.92,
                       depth=0,
                       min_depth=1,
                       max_depth=5)

    if depth > max_depth  || (min_depth <= depth && rand() < leaf_prob)
        leaf = rand(spec.leaves)()
        if leaf isa SymbolicUtils.Sym
            push!(inputs, leaf)
        elseif leaf isa Pair
            foreach(i->push!(inputs, i), leaf[1])
            return leaf[2]
        end
        return leaf
    end
    arity, f = rand(spec.funcs)
    args = [gen_rand_expr(inputs,
                          spec=spec,
                          leaf_prob=leaf_prob,
                          depth=depth+1,
                          min_depth=min_depth,
                          max_depth=max_depth) for i in 1:arity]
    try
        return f(args...)
    catch err
        if err isa DomainError || err isa DivideError || err isa MethodError ||
            err isa SymbolicUtils.SpecialFunctions.AmosException
            return gen_rand_expr(inputs,
                                 spec=spec,
                                 leaf_prob=leaf_prob,
                                 depth=depth,
                                 min_depth=min_depth,
                                 max_depth=max_depth)
        else
            @show f
            @show arity
            @show args
            rethrow(err)
        end
    end
end

struct Errored
    err
end

function fuzz_test(ntrials, spec, simplify=simplify;kwargs...)
    inputs = Set()
    expr = gen_rand_expr(inputs; spec=spec, kwargs...)
    inputs = collect(inputs)
    code = try
        SymbolicUtils.Code.toexpr(expr)
    catch err
        @show expr
        rethrow(err)
    end
    unsimplifiedstr = """
    function $(tuple(inputs...))
        $(sprint(io->print(io, code)))
    end
    """

    simplifiedstr = """
    function $(tuple(inputs...))
        $(sprint(io->showraw(io, simplify(expr))))
    end
    """
    f = include_string(Main, unsimplifiedstr)
    g = include_string(Main, simplifiedstr)

    for i=1:ntrials
        args = [spec.input(i) for i in inputs]
        unsimplified = try
            Base.invokelatest(f, args...)
        catch err
            Errored(err)
        end

        simplified  = try
            Base.invokelatest(g, args...)
        catch err
            Errored(err)
        end
        if unsimplified isa Errored
            if !(simplified isa Errored)
                @test_skip false
                @goto print_err
            end
            @test true
        elseif simplified isa Errored
            if !(unsimplified isa Errored)
                @test_skip false
                @goto print_err
            end
            @test true
        elseif isnan(unsimplified)
            if !isnan(simplified)
                @test_skip false
                @goto print_err
            end
            @test true
        elseif isnan(simplified)
            if !isnan(unsimplified)
                @test_skip false
                @goto print_err
            end
            @test true
        else
            if !(unsimplified ≈ simplified)
                @test_skip false
                @goto print_err
            end
            @test true
        end
        continue

        @label print_err
        println("""Test failed for expression
                    $(sprint(io->showraw(io, expr))) = $unsimplified
                Simplified:
                    $(sprint(io->showraw(io, simplify(expr)))) = $simplified
                Inputs:
                    $inputs = $args
                """)
    end
end

leaves = [(@syms a b c d e g)..., a^3, a^-2, b^2, b^(-1), big.((3//5, 0//2, -1//2, 1//2, 1//2+2im))...]
function gen_expr(lvl=5)
    if lvl == 0
        x = rand(leaves)
        x, x
    elseif rand() < 0.5
        f = rand((+, *))
        n = rand(1:5)
        args = [gen_expr(lvl-1) for i in 1:n]

        Term{Number}(f, first.(args)), f(last.(args)...)
    else
        f = rand((-,/))
        l = gen_expr(lvl-1)
        r = gen_expr(lvl-1)
        if f === (/) && r[2] isa Number && iszero(r[2])
            return gen_expr(lvl)
        end
        args = [l, r]

        Term{Number}(f, first.(args)), f(last.(args)...)
    end
end

test_dict = Dict{Any, Rational{BigInt}}(a=>1,b=>-1,c=>2,d=>-2,e=>5//3,g=>-2//3)
function fuzz_addmulpow(lvl, d=test_dict)
    l, r = gen_expr()
    rl = try
        substitute(l, d)
    catch err
        err
    end
    rr = try
        substitute(r, d)
    catch err
        err
    end

    if !(rl isa Number) && rr isa Number
        return # lhs errored, rhs did not
    end
    if rl isa Number || rr isa Number
        if isequal(rl, rr)
            @test true
        else
            println("Weird bug here:")
            @show r l
            @show rl rr
            @test false
        end
    end
end


text:
import LinearAlgebra: det, logdet, diagm, Diagonal, diag

const ∇ScalarDiag = Diagonal{<:∇Scalar}

@explicit_intercepts Diagonal Tuple{∇AbstractMatrix}
function ∇(
    ::Type{Diagonal},
    ::Type{Arg{1}},
    p,
    Y::∇ScalarDiag,
    Ȳ::∇ScalarDiag,
    X::∇AbstractMatrix,
)
    X̄ = zeroslike(X)
    copyto!(view(X̄, diagind(X)), Ȳ.diag)
    return X̄
end
function ∇(
    X̄::∇AbstractMatrix,
    ::Type{Diagonal},
    ::Type{Arg{1}},
    p,
    Y::∇ScalarDiag,
    Ȳ::∇ScalarDiag,
    X::∇AbstractMatrix,
)
    X̄_diag = view(X̄, diagind(X̄))
    broadcast!(+, X̄_diag, X̄_diag, Ȳ.diag)
    return X̄
end

# NOTE: diagm can't go through the @explicit_intercepts machinery directly because as of
# Julia 0.7, its methods are not sufficiently straightforward; we need to dispatch on one
# of the parameters in the parametric type of diagm's one argument. However, we can cheat
# a little bit and use an internal helper function _diagm that has simple methods that
# dispatch to diagm when no arguments are Nodes, and we'll extend diagm to dispatch to
# _diagm when it receives arguments that are nodes. _diagm can go through the intercepts
# machinery, so it knows how to deal.

# TODO: Possibly we should overload `Pair` so that it constructs a `Node{Pair}` then this
# would hit sensitivities that we have defined via ChainRules.

_diagm(x::∇AbstractVector, k::Integer=0) = diagm(k => x)
LinearAlgebra.diagm(x::Pair{<:Integer, <:Node{<:∇AbstractVector}}) = _diagm(last(x), first(x))

@explicit_intercepts _diagm Tuple{∇AbstractVector}
function ∇(
    ::typeof(_diagm),
    ::Type{Arg{1}},
    p,
    Y::∇AbstractMatrix,
    Ȳ::∇AbstractMatrix,
    x::∇AbstractVector,
)
    return copyto!(similar(x), view(Ȳ, diagind(Ȳ)))
end
function ∇(
    x̄::∇AbstractVector,
    ::typeof(_diagm),
    ::Type{Arg{1}},
    p,
    Y::∇AbstractMatrix,
    Ȳ::∇AbstractMatrix,
    x::∇AbstractVector,
)
    return broadcast!(+, x̄, x̄, view(Ȳ, diagind(Ȳ)))
end
@explicit_intercepts _diagm Tuple{∇AbstractVector, Integer} [true, false]
function ∇(
    ::typeof(_diagm),
    ::Type{Arg{1}},
    p,
    Y::∇AbstractMatrix,
    Ȳ::∇AbstractMatrix,
    x::∇AbstractVector,
    k::Integer,
)
    return copyto!(similar(x), view(Ȳ, diagind(Ȳ, k)))
end
function ∇(
    x̄::∇AbstractVector,
    ::typeof(_diagm),
    ::Type{Arg{1}},
    p,
    Y::∇AbstractMatrix,
    Ȳ::∇AbstractMatrix,
    x::∇AbstractVector,
    k::Integer,
)
    return broadcast!(+, x̄, x̄, view(Ȳ, diagind(Ȳ, k)))
end


text:
# The old-gdemo model.
@model gdemo(x, y) = begin
  s ~ InverseGamma(2, 3)
  m ~ Normal(0, sqrt(s))
  x ~ Normal(m, sqrt(s))
  y ~ Normal(m, sqrt(s))
  return s, m
end

@model gdemo_d() = begin
  s ~ InverseGamma(2, 3)
  m ~ Normal(0, sqrt(s))
  1.5 ~ Normal(m, sqrt(s))
  2.0 ~ Normal(m, sqrt(s))
  return s, m
end

gdemo_default = gdemo_d()

@model MoGtest(D) = begin
    mu1 ~ Normal(1, 1)
    mu2 ~ Normal(4, 1)
    z1 ~ Categorical(2)
    if z1 == 1
        D[1] ~ Normal(mu1, 1)
    else
        D[1] ~ Normal(mu2, 1)
    end
    z2 ~ Categorical(2)
    if z2 == 1
        D[2] ~ Normal(mu1, 1)
    else
        D[2] ~ Normal(mu2, 1)
    end
    z3 ~ Categorical(2)
    if z3 == 1
        D[3] ~ Normal(mu1, 1)
    else
        D[3] ~ Normal(mu2, 1)
    end
    z4 ~ Categorical(2)
    if z4 == 1
        D[4] ~ Normal(mu1, 1)
    else
        D[4] ~ Normal(mu2, 1)
    end
    z1, z2, z3, z4, mu1, mu2
end

MoGtest_default = MoGtest([1.0 1.0 4.0 4.0])

# Declare empty model to make the Sampler constructor work.
@model empty_model() = begin x = 1; end

# A collection of models for which the mean-of-means for the posterior should
# be same.
@model function gdemo1(x = 10 * ones(2), ::Type{TV} = Vector{Float64}) where {TV}
    # `dot_assume` and `observe`
    m = TV(undef, length(x))
    m .~ Normal()
    x ~ MvNormal(m, 0.5 * ones(length(x)))
end

@model function gdemo2(x = 10 * ones(2), ::Type{TV} = Vector{Float64}) where {TV}
    # `assume` with indexing and `observe`
    m = TV(undef, length(x))
    for i in eachindex(m)
        m[i] ~ Normal()
    end
    x ~ MvNormal(m, 0.5 * ones(length(x)))
end

@model function gdemo3(x = 10 * ones(2))
    # Multivariate `assume` and `observe`
    m ~ MvNormal(length(x), 1.0)
    x ~ MvNormal(m, 0.5 * ones(length(x)))
end

@model function gdemo4(x = 10 * ones(2), ::Type{TV} = Vector{Float64}) where {TV}
    # `dot_assume` and `observe` with indexing
    m = TV(undef, length(x))
    m .~ Normal()
    for i in eachindex(x)
        x[i] ~ Normal(m[i], 0.5)
    end
end

# Using vector of `length` 1 here so the posterior of `m` is the same
# as the others.
@model function gdemo5(x = 10 * ones(1))
    # `assume` and `dot_observe`
    m ~ Normal()
    x .~ Normal(m, 0.5)
end

@model function gdemo6()
    # `assume` and literal `observe`
    m ~ MvNormal(2, 1.0)
    [10.0, 10.0] ~ MvNormal(m, 0.5 * ones(2))
end

@model function gdemo7(::Type{TV} = Vector{Float64}) where {TV}
    # `dot_assume` and literal `observe` with indexing
    m = TV(undef, 2)
    m .~ Normal()
    for i in eachindex(m)
        10.0 ~ Normal(m[i], 0.5)
    end
end

@model function gdemo8()
    # `assume` and literal `dot_observe`
    m ~ Normal()
    [10.0, ] .~ Normal(m, 0.5)
end

@model function _prior_dot_assume(::Type{TV} = Vector{Float64}) where {TV}
    m = TV(undef, 2)
    m .~ Normal()

    return m
end

@model function gdemo9()
    # Submodel prior
    m = @submodel _prior_dot_assume()
    for i in eachindex(m)
        10.0 ~ Normal(m[i], 0.5)
    end
end

@model function _likelihood_dot_observe(m, x)
    x ~ MvNormal(m, 0.5 * ones(length(m)))
end

@model function gdemo10(x = 10 * ones(2), ::Type{TV} = Vector{Float64}) where {TV}
    m = TV(undef, length(x))
    m .~ Normal()

    # Submodel likelihood
    @submodel _likelihood_dot_observe(m, x)
end

const gdemo_models = (gdemo1(), gdemo2(), gdemo3(), gdemo4(), gdemo5(), gdemo6(), gdemo7(), gdemo8(), gdemo9(), gdemo10())


text:
using .LinearMaps, .IterativeSolvers
export CGIndirectKKTSolver, MINRESIndirectKKTSolver
mutable struct IndirectReducedKKTSolver{TP, TA, T} <: AbstractKKTSolver
    m::Integer
    n::Integer
    P::TP
    A::TA
    σ::T
    ρ::Vector{T}
    tol_constant::T
    tol_exponent::T
    # Memory and statistics for cg
    solver_type::Symbol
    previous_solution::Vector{T}
    iteration_counter::Int
    multiplications::Vector{Int}
    tmp_n::Vector{T} # n-dimensional used to avoid memory allocation
    tmp_m::Vector{T} # m-dimensional used to avoid memory allocation

    function IndirectReducedKKTSolver(P::TP, A::TA, σ::T, ρ;
        solver_type::Symbol=:CG, tol_constant::T=T(1.0), tol_exponent::T=T(1.5)
        ) where {TP, TA, T}
        m, n = size(A)
        if isa(ρ, T)
            ρ = ρ*ones(T, m)
        end
        @assert eltype(P) == eltype(A) == T "Inconsistent element types."
        @assert solver_type == :CG || solver_type == :MINRES "Solver symbol must be either :CG or :MINRES"
        new{TP, TA, T}(m, n, P, A, σ, ρ,
            tol_constant, tol_exponent,
            solver_type,
            zeros(T, n), 1, zeros(Int, 0), zeros(T, n), zeros(T, m))
    end
end

function solve!(S::IndirectReducedKKTSolver, y::AbstractVector{T}, x::AbstractVector{T}) where {T}
    # Solves the (KKT) linear system
    # | P + σI     A'  |  |y1|  =  |x1|
    # | A        -I/ρ  |  |y2|  =  |x2|
    # x1,y1: R^n, x2/y2: R^m
    # where [y1; y2] := y, [x1; x2] := x

    # In particular we perform cg/minres on the reduced system
    # (P + σΙ + Α'ρΑ)y1 = x1 + A'ρx2
    # And then recover y2 as
    # y2 = ρ(Ay1 - x2)

    x1 = view(x, 1:S.n); y1 = view(y, 1:S.n)
    x2 = view(x, S.n+1:S.n+S.m); y2 = view(y, S.n+1:S.n+S.m)

    # Form right-hand side for cg/minres: y1 = x1 + A'ρx2
    @. y2 = S.ρ*x2
    mul!(y1, S.A', y2)
    y1 .+= x1

    push!(S.multiplications, 0)
    function reduced_mul!(y::AbstractVector{T}, x::AbstractVector{T}) where {T}
        # y = (P + σΙ + Α'ρΑ)x
        mul!(S.tmp_m, S.A, x)
        @. S.tmp_m .*= S.ρ
        mul!(S.tmp_n, S.A', S.tmp_m)
        axpy!(S.σ, x, S.tmp_n)
        mul!(y, S.P, x)
        axpy!(one(T), S.tmp_n, y)
        S.multiplications[end] += 1
        return y
    end
    L = LinearMap(reduced_mul!, S.n; ismutating=true, issymmetric=true)
    if S.solver_type == :CG
        cg!(S.previous_solution, L, y1, tol=get_tolerance(S)/norm(y1))
    elseif S.solver_type == :MINRES
        init_residual = norm(L*S.previous_solution - y1)
        minres!(S.previous_solution, L, y1, tol=get_tolerance(S)/init_residual)
    end
    # Sanity check for tolerance
    # might not always hold for MINRES, as its termination criterion is approximate, (see its documentation)
    # @assert get_tolerance(S) > norm(L*S.previous_solution - y1)
    copyto!(y1, S.previous_solution)

    # y2 = Ay1 - x2
    mul!(y2, S.A, y1)
    axpy!(-one(T), x2, y2)
    @. y2 .*= S.ρ

    S.iteration_counter += 1

    return y
end

mutable struct IndirectKKTSolver{TP, TA, T} <: AbstractKKTSolver
    m::Integer
    n::Integer
    P::TP
    A::TA
    σ::T
    ρ::Vector{T}
    tol_constant::T
    tol_exponent::T
    # Memory and statistics for cg
    solver_type::Symbol
    previous_solution::Vector{T}
    iteration_counter::Int
    multiplications::Vector{Int}
    tmp_n::Vector{T} # n-dimensional used to avoid memory allocation
    tmp_m::Vector{T} # m-dimensional used to avoid memory allocation

    function IndirectKKTSolver(P::TP, A::TA, σ::T, ρ; solver_type::Symbol=:MINRES,
        tol_constant::T=T(1.0), tol_exponent::T=T(1.5)
        ) where {TP, TA, T}
        m, n = size(A)
        if isa(ρ, T)
            ρ = ρ*ones(T, m)
        end
        @assert eltype(P) == eltype(A) == T "Inconsistent element types."
        @assert solver_type == :MINRES "Solver symbol must be :MINRES"
        new{TP, TA, T}(m, n, P, A, σ, ρ,
            tol_constant, tol_exponent,
            solver_type,
            zeros(T, m + n), 1, zeros(Int, 0), zeros(T, n), zeros(T, m))
    end
end

function solve!(S::IndirectKKTSolver, y::AbstractVector{T}, x::AbstractVector{T}) where {T}
    # Solves the (KKT) linear system
    # | P + σI     A'  |  |y1|  =  |x1|
    # | A        -I/ρ  |  |y2|  =  |x2|
    # x1,y1: R^n, x2/y2: R^m
    # where [y1; y2] := y, [x1; x2] := x
    push!(S.multiplications, 0)
    function kkt_mul!(y::AbstractVector{T}, x::AbstractVector{T}) where {T}
        # Performs
        # |y1| = | P + σI     A'  |  |x1|
        # |y2| = | A        -I/ρ  |  |x2|
        x1 = view(x, 1:S.n); y1 = view(y, 1:S.n)
        x2 = view(x, S.n+1:S.n+S.m); y2 = view(y, S.n+1:S.n+S.m)

        mul!(S.tmp_n, S.A', x2)
        axpy!(S.σ, x1, S.tmp_n)
        mul!(y1, S.P, x1)
        axpy!(one(T), S.tmp_n, y1)

        @. y2 = -x2/S.ρ
        mul!(S.tmp_m, S.A, x1)
        axpy!(one(T), S.tmp_m, y2)

        S.multiplications[end] += 1
        return y
    end
    L = LinearMap(kkt_mul!, S.n + S.m; ismutating=true, issymmetric=true)
    if S.solver_type == :MINRES
        init_residual = norm(L*S.previous_solution - x)
        minres!(S.previous_solution, L, x, tol=get_tolerance(S)/init_residual)
    end
    # Sanity check for tolerance
    # might not always hold for MINRES, as its termination criterion is approximate, (see its documentation)
    # @assert get_tolerance(S) > norm(L*S.previous_solution - y)
    copyto!(y, S.previous_solution)

    S.iteration_counter += 1

    return y
end

function update_rho!(S::Union{IndirectReducedKKTSolver, IndirectKKTSolver}, ρ)
    S.ρ .= ρ
end

function get_tolerance(S::Union{IndirectReducedKKTSolver, IndirectKKTSolver})
    return S.tol_constant/S.iteration_counter^S.tol_exponent
end

# For convenience we define CG and Minres wrappers around IndirectKKTSolver and IndirectReducedKKTSolver
struct CGIndirectKKTSolver{TP, TA, T} <: AbstractKKTSolver
    indirect_kktsolver::IndirectReducedKKTSolver{TP, TA, T}
    function CGIndirectKKTSolver(P::TP, A::TA, σ::T, ρ; tol_constant::T=T(1.0), tol_exponent::T=T(1.5)) where {TP, TA, T}
        new{TP, TA, T}(IndirectReducedKKTSolver(P, A, σ, ρ; solver_type = :CG, tol_constant = tol_constant, tol_exponent = tol_exponent))
    end
end

struct MINRESIndirectKKTSolver{TP, TA, T} <: AbstractKKTSolver
    indirect_kktsolver::IndirectKKTSolver{TP, TA, T}
    function MINRESIndirectKKTSolver(P::TP, A::TA, σ::T, ρ; tol_constant::T=T(1.0), tol_exponent::T=T(1.5)) where {TP, TA, T}
        new{TP, TA, T}(IndirectKKTSolver(P, A, σ, ρ; solver_type = :MINRES, tol_constant = tol_constant, tol_exponent = tol_exponent))
    end
end

update_rho!(S::Union{CGIndirectKKTSolver, MINRESIndirectKKTSolver}, ρ) = update_rho!(S.indirect_kktsolver, ρ)
get_tolerance(S::Union{CGIndirectKKTSolver, MINRESIndirectKKTSolver}, ρ) = get_tolerance(S.indirect_kkt_solver, ρ)
solve!(S::Union{CGIndirectKKTSolver, MINRESIndirectKKTSolver}, y::AbstractVector{T}, x::AbstractVector{T}) where {T} = solve!(S.indirect_kktsolver, y, x)


text:
#     These functions deal with specified websocket protocols.
#      Included in server_functions.jl

function ws_test_protocol(ws::WebSockets.WebSocket)
    id = cid = ""
    while isopen(ws)
        id, msg = wsmsg_listen(ws)
        cid = "ws_test_protocol $(id)\t"
        # WEBSOCKETS_SUBPROTOCOL[id] = ws
        sendmsg = "Special welcome. Our ref.: $(id)"
        clog(cid,  " Sending: $sendmsg\n")
        writeto(id, sendmsg)
        clog(cid, "Received: \n", :yellow, :bold, msg, "\n")
        if startswith(msg, "Ping me!")
            sendmsg = "No ping for you. Our ref.: $(id). Will close on next received message."
            clog(cid,  " Sending: $sendmsg\n")
            writeto(id, sendmsg)
        elseif startswith(msg, "ws2 echo: No ping for you.")
            clog(cid,  "Closing from server\n")
            close(ws)
        end
    end
    clog(cid, "Websocket on subprotocol $(id) was closed; exits handler.\n")
    nothing
end

function ws_test_binary(ws::WebSockets.WebSocket)
    id = cid = ""
    width = Int16(0)
    height = Int16(0)
    byteperpixel = Int16(4)
    ub = 0
    while isopen(ws)
        id, data = wsmsg_listen_binary(ws)
        cid = "ws_test_binary $(id)\t"
        clog(cid,  " Handler is running \n")
        ub = length(data)
        if ub == 6
            # Don't do this like this at home perhaps. Dealing with net endianness by reversing the order of bytes AND arguments..
            byteperpixel, height, width = reinterpret(Int16, reverse(data))
            clog(cid, "Received 6 bytes = 3xInt16 custom header: \n\t\t\t\t\t"
                    , :yellow, " height: ", :bold, height
                    , :normal, :yellow, " width: ", :bold, width
                    , :normal, :yellow, " byteperpixel: ", byteperpixel, "\n")
        elseif ub == Int(width) * height * byteperpixel
            clog(cid, "Received image data for manipulation, ", div(length(data), 1024),  " kB\n")
            tic()
            ired =  1:4:ub
            igreen =  2:4:ub
            iblue =  3:4:ub
            ialpha = 4:4:ub
            itpquarter = 1:div(ub,4)
            fill!(view(data, ialpha), 0xee)
            fill!(view(data, iblue), 0xff)
            fill!(view(data, intersect(itpquarter, ired)), 0x88)
            rand!(view(data, ired), [0x00, 0x88, 0xff])
            hedata = UInt16.([height, width, 4])
            clog(id, "Manipulated image and made ready in $(toq()*1000) ms. \n\t\tSending header, then image.\n")
            writeto(id,  reinterpret(UInt8, hedata))
            writeto(id,  data)
        else
            clog(cid, "Received unexpected data: \t", :yellow, :bold, length(data), "-element ", typeof(data), "\n")
            clog(cid, "expected length: ",  Int(width) * height * byteperpixel
                    , :yellow, " height: ", :bold, height
                    , :normal, :yellow, " width: ", :bold, width
                    , :normal, :yellow, " byteperpixel: ", byteperpixel, "\n")
        end
    end
    clog(cid, "Websocket on subprotocol $(id) was closed; exits handler.\n")
    nothing
end
"""
Listens for the next binary websocket message
"""
function wsmsg_listen_binary(ws::WebSockets.WebSocket)
    id = cid = ""
    data = UInt8[]
    try
        # Holding here. Cleanup code could work with InterruptException and Base.throwto.
        # This is not considered necessary for the test (and doesn't really work well).
        # So, this may continue running after the test is finished.
        id = String(read(ws))
        data , t, allocbytes = @timed read(ws)
        cid = "wsmsg_listen_binary $(id)\t"
        clog(cid,  " Binary listening...\n")

        if !haskey(WEBSOCKETS,id)
            clog(:red,"Pushing WS...\n")
            WEBSOCKETS[id] = ws
            WEBSOCKETS_BINARY[id] = ws
            RECEIVED_WS_MSGS_LENGTH[id] = Vector{Int64}()
            RECEIVED_WS_MSGS_TIME[id] = Vector{Float64}()
            RECEIVED_WS_MSGS_ALLOCATED[id] = Vector{Int64}()
        end
        push!(RECEIVED_WS_MSGS_LENGTH[id], length(data))
        push!(RECEIVED_WS_MSGS_TIME[id], t)
        push!(RECEIVED_WS_MSGS_ALLOCATED[id], allocbytes)
    catch e
        if typeof(e) == WebSockets.WebSocketClosedError
            clog(cid, :green, " Websocket was or is being closed.\n")
        else
            clog(cid, "Caught exception..\n", "\t\t", :red, e, "\n")
            if isopen(ws)
                clog(cid, "Closing\n")
                close(ws)
            end
        end
    end
    return id, data
end

return nothing


text:
using SymbolicRegression
using Bijectors
I = collect(0:0.1:10)
f(x) = 0.2exp(-0.01*x)*x
Y = f.(I)

options = SymbolicRegression.Options(
    binary_operators=(+, *, /, -),
    unary_operators=(exp,),
    npopulations=20
)

hallOfFame = EquationSearch(I', Y, niterations=5, options=options, numprocs=4)
dominating = calculateParetoFrontier(I, Y, hallOfFame, options)
eqn = node_to_symbolic(dominating[end].tree, options)


text:
include("Helpers.jl")

include("CachedArrays.jl")

include("CachedSubVectors.jl")

include("CachedValues.jl")


text:
using Documenter
using MixedModelsExtras

makedocs(
    root = joinpath(dirname(pathof(MixedModelsExtras)), "..", "docs"),
    sitename = "MixedModelsExtras",
    doctest = true,
    pages = [
        "index.md",
    ],
)

deploydocs(repo = "github.com/palday/MixedModelsExtras.jl.git", push_preview = true)


text:
using TextPlots

plot(x -> cos(x), -1:1)
#> cos(x)
#>        1 ⡤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠊⠉⠉⠉⠉⠉⠉⠑⠒⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠒⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠌⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠡⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⡐⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢂⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⡠⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠡⡀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⢐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⡀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠠⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠄⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⢀⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢂⠀⠀⠀⢸
#>          ⡇⠀⢀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢂⠀⠀⢸
#>          ⡇⢀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢂⠀⢸
#>          ⡇⡐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢂⢸
#>  0.54030 ⠓⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠚
#>         -1                                                            1

plot([x -> cos(x), x -> cos(x + pi)], 0:5)
#> cos(x), cos(x + pi)
#>        1 ⡤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤
#>          ⡇⠉⠉⠉⠒⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠔⠊⠉⠉⠉⠉⠉⠑⠒⠤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠉⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⠤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠌⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢄⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⡐⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠡⡀⠀⠀⠀⠀⠀⢀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠢⡀⠀⠀⠠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢂⠀⠀⠀⡠⠁⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢄⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⣀⠌⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠊⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠉⢂⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠔⠁⠀⠀⠐⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠌⠀⠀⠀⠑⡀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠡⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡐⠁⠀⠀⠀⠀⠀⠈⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠊⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⣀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⠤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⣀⣀⣀⠤⠒⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⠢⢄⣀⣀⣀⣀⣀⡠⠤⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>       -1 ⠓⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠚
#>          0                                                            5

plot([1, 3, 5, 7], [13, 11, 9, 7])
#> scatter plot
#>       13 ⡤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤
#>          ⡇⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⢸
#>          ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
#>        7 ⠓⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠚
#>          1                                                            7

max_stars_repo_name: floswald/JuliaByExample

max_issues_repo_name: floswald/JuliaByExample

max_forks_repo_name: floswald/JuliaByExample


text:
using BCTRNN
using DiffEqSensitivity
using OrdinaryDiffEq
import DiffEqFlux: FastChain, FastDense
import Flux: ClipValue, ADAM

# Not in Project.toml
using Plots
gr()

include("half_cheetah_data_loader.jl")

function train_cheetah(epochs, solver=nothing; sensealg=nothing,
  T=Float32, model_size=5, batchsize=1, seq_len=32, normalise=true,
  kwargs...)

  train_dl, test_dl, _, _ = get_2d_dl(T; batchsize, seq_len, normalise=true)
  @show size(first(train_dl)[1])
  @show size(first(train_dl)[1][1])

  f_in = 17
  f_out = 17
  n_neurons = model_size
  n_sens = n_neurons
  n_out = n_neurons

  model = FastChain(BCTRNN.Mapper(f_in),
                    BCTRNN.LTC(f_in, n_neurons, solver, sensealg; n_sens, n_out),
                    FastDense(n_out, f_out))

  cb = BCTRNN.MyCallback(T; ecb=mycb, nepochs=epochs, nsamples=length(train_dl))
  #opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.02))
  opt = BCTRNN.ClampBoundOptim(BCTRNN.get_bounds(model,T)..., ClipValue(T(0.8)), ADAM(T(0.005)))
  BCTRNN.optimize(model, BCTRNN.loss_seq, cb, opt, train_dl, epochs, T), model
end
#1173.351351 seconds (1.02 G allocations: 65.414 GiB, 1.82% gc time, 0.51% compilation time)
train_cheetah(30, Tsit5(); sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)), model_size=8, batchsize=10, abstol=1e-4, reltol=1e-3
)
train_cheetah(30, AutoTsit5(Rosenbrock23(autodiff=false)); sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)), model_size=6, batchsize=10, abstol=1e-4, reltol=1e-3
)


text:
immutable GradientState
    vol::Vol
    multiplierL1::Float64
    multiplierL2::Float64
end

abstract Layer

outputSize(self::Layer) = size(self.outputActivation.w)
inputSize(self::Layer) = size(self.inputActivation.w)

## +----------+
## |          |
## |          |
## |          |
## |          |
## |          |
## +----------+
## |..........|
## |..........|
## |..........<------------+  Input Activation
## |..........|            |
## |..........|            |
## +----------+            |
##                         |
##                         |
## +----------+            |
## |..........|            |
## |..........|            |
## |..........>------------+  Input Activation
## |..........|
## |..........|
## +----------+
## |          |
## |          |
## |          <------------+  Input Activation
## |          |            |
## |          |            |
## +----------+            |
##                         |
##                         |
##                         |
## +----------+            |
## |          |            |
## |          |            |
## |          >------------+  Output Activation
## |          |
## |          |
## +----------+
## |..........|
## |..........|
## |..........>------------+  Activation from lower levels
## |..........|            |
## |..........|            |
## +----------+            |


text:
include("test_utils.jl")
include("test_models.jl")
include("test_optimizers.jl")
include("test_fine_tuning.jl")


text:
import Pkg
Pkg.instantiate() # downloads all dependencies for the current project

Pkg.add("JSON")
Pkg.add("HTTP")
Pkg.add("JSONTables")

using HTTP, JSON, JSONTables

# treat a json object of arrays or array of objects as a "table"
jtable = jsontable(json_source)

# turn json table into DataFrame
df = DataFrame(jtable)

# turn DataFrame back into json object of arrays
objecttable(df)

# turn DataFrame back into json array of objects
arraytable(df)


text:
# test_util.jl

function test_util()
    @testset "$(rpad("Index reversal",80))" begin
        test_reverse([1,2,3])
        test_reverse(OffsetArray([1,2,3], -2:0))
    end
end

function test_reverse(a)
    b = reverse_indices(a)
    for i in eachindex(a)
        @test a[i] == b[-i]
    end
end


text:
using PAPI, Test

include("eventset.jl")
#include("counting.jl")
include("serialization.jl")


text:
#=
    OLD OUT DATED CODE
=#

struct Disk_old
    identifier::String
    color::String
    location::Tuple{Float64, Float64}
    radius::Float64
    parameters::Vector{Disk_old}
end

"""
    add_disk(rt::Disk_old, parent, child = [])

The basic tree creation for disk operads. Specifying the new Disk_old, the parent identifier, and any potential children.
"""
function add_disk(rt::Disk_old, parent, child = [])
    found = false

    if parent == rt.identifier
        push!(rt.parameters, child)
        return (rt, true)
    end

    for i in rt.parameters
        found_m = false
        (i, found_m) = add_disk(i, parent, child)
        if found_m == true
            found = true
        end
    end

    return rt, found
end

"""
    draw_disk(to_draw, origin, dash)
Basics for drawing a disk. THe origin is what the disk will be relative to, and dash will be used for circles with outlines
"""
function draw_disk(to_draw, origin, dash)
    center = (origin[1] + (to_draw.location[1] * cosd(to_draw.location[2])), origin[2] + (to_draw.location[1] * sind(to_draw.location[2])))
    return compose(context(), circle(center[1], center[2], to_draw.radius), fill(to_draw.color))
    #if dash == 0
    #    circle = compose(context(), circle(center[1], center[2], to_draw.radius), fill(to_draw.color))
    #elseif dash == 1
    #    circle = compse(context(), arc(center[1], center[2], to_draw.radius, range(0, step=pi/4, length=8), range(pi/8, step=pi/4, length = 8)), fill(to_draw.color))
    #else
    #    circle = 0;
    #end
end

max_stars_repo_path: src/Old Code/Disk_Old.jl

max_issues_repo_path: src/Old Code/Disk_Old.jl

max_forks_repo_path: src/Old Code/Disk_Old.jl


text:
using Test
using MARY_fRG.Basics
using Printf

@testset "测试三角形" begin
    tri1 = Triangle([Point2D(1., 0.), Point2D(0., 1), Point2D(0., 0.)])
    println(tri1)
    println(Point2D(1., 0.))
    @test typeof(tri1) <: AbstractPolygon
    tri2 = Triangle(Point2D(1., 0.), Point2D(0., 1), Point2D(0., 0.))
    println(tri2)
    @test area(tri1) == 0.5
    @test area(tri1) == tri1.area
end

@testset "测试直角三角形" begin
    tri1 = RtTriangle([Point2D(0., 0.), Point2D(0., 1), Point2D(1., 0.)])
    @test area(tri1) == 0.5
    @test area(tri1) == tri1.area
    tri2 = RtTriangle(Point2D(1., 1.), Point2D(0., 0), Point2D(2., 0.))
    @test area(tri2) == 1.0
    @test area(tri2) == tri2.area
end

@testset "测试等边三角形" begin
    tri1 = EqTriangle([Point2D(0., 0.), Point2D(1., sqrt(3)), Point2D(2., 0.)])
    @test area(tri1) == sqrt(3)
    @test area(tri1) == tri1.area
    tri2 = EqTriangle(Point2D(1., 1.), Point2D(2., 1+sqrt(3)), Point2D(3., 1.))
    @test area(tri2) == sqrt(3)
    @test area(tri2) == tri2.area
end


text:
# Skew-normal simulation example for c-posterior for mixtures
module Suite
using Distributions
using PyPlot
using HDF5, JLD
draw() = (pause(0.001); get_current_fig_manager()[:window][:raise_]())

# Settings
ns = [20,100,500,2000,10000]  # sample sizes n to use
alphas = [Inf,100]  # robustification params alpha to use
nreps = 5  # number of times to run the simulation
n_samples = 10^5  # number of MCMC iterations
nburn = int(n_samples/10)  # number of iterations to discard as burn-in
from_file = false  # load previous results from file

# Model parameters and sampler code
include("setup.jl")

nns = length(ns)
for (i_a,alpha) in enumerate(alphas)

    if from_file
        k_posteriors = load("k_posteriors-alpha=$alpha.jld","k_posteriors")
    else
        k_posteriors = zeros(m,nreps,nns)
        for (i_n,n) in enumerate(ns)
            for rep in 1:nreps
                srand(n+rep) # Reset RNG
                x = [(rand()<p1 ? skewrnd(1,l1,s1,a1) : skewrnd(1,l2,s2,a2))[1]::Float64 for j=1:n] # Sample data

                # Run sampler
                zeta = (1/n)/(1/n + 1/alpha)
                p,theta,k_r,v_r,art,arv,m_r,l_r = sampler(x,n_samples,m,c,sigma,zeta)

                # Compute posterior on k
                edges,counts = hist(k_r[nburn+1:end],[0:m])
                k_posteriors[:,rep,i_n] = counts/(n_samples-nburn)
            end
        end
        save("k_posteriors-alpha=$alpha.jld","k_posteriors",k_posteriors)
    end

    # Prior on k: p(k) \propto Binomial(k|m,lambda/m)I(k>0).
    pk = pdf(Binomial(m,lambda/m),1:m); pk = pk/sum(pk)

    # Plot
    figure(i_a,figsize=(8,3.2)); clf(); hold(true)
    subplots_adjust(bottom=0.2)
    colors = "bgyrm"
    shapes = "ds^v*"
    plot(1:m,pk,"ko-",label="prior",markersize=8,linewidth=2)
    for (i_n,n) in enumerate(ns)
        plot(1:m,vec(mean(k_posteriors[:,:,i_n],2)),"$(colors[i_n])$(shapes[i_n])-",label="n = $n",markersize=8,linewidth=2)
    end
    xlim(0,m)
    xticks(1:m)
    if alpha<Inf; title("\$\\mathrm{Coarsened\\,posterior}\$",fontsize=17)
    else; title("\$\\mathrm{Standard\\,posterior}\$",fontsize=17)
    end
    xlabel("\$k \\,\\,\\mathrm{(\\#\\,of\\,components)}\$",fontsize=17)
    ylabel("\$\\Pi(k|\\mathrm{data})\$",fontsize=15)
    legend(loc="upper right",numpoints=1)
    draw()
    savefig("skewsuite-alpha=$alpha.png",dpi=150)

end

end # module

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
### Training utilities

export training!

function training!(
	T::InvertibleNetwork,
	fval::Array{Float32, 1},
	loss::Function, Ztrain::Array{Float32, 4},
	nepochs::Int64, batchsize::Int64,
	opt::Optimiser;
	verbose::Bool = true, verbose_b::Bool = false,
	randinput::Bool = false,
	gradclip::Bool = false)

	trainsize = size(Ztrain, 4)
	nbatches = Int64(round(trainsize/batchsize))

	# Loop over iterations
	for n = 1:nepochs

	    # Shuffling data
	    idx_train_perm = randperm(trainsize)

	    # Loop over batches
	    for idx_batch = 1:nbatches

			# Selecting batch
	        idx = idx_train_perm[(idx_batch-1)*batchsize+1:idx_batch*batchsize]
	        Z = Ztrain[:, :, :, idx]
			randinput && (Z += 1f-2*randn(Float32, size(Z)))

	        # Evaluate objective and backpropagating
	        l = loss(Z)
	        append!(fval, l)

	        # Update network
			for p in get_params(T)
				gradclip && (p.grad ./= (norm(p.grad)+1f-10))
				update!(opt, p.data, p.grad)
			end

	        # Clear gradients
	        clear_grad!(T)

		    # Print msg
		    verbose_b && println("         Batch: [", idx_batch, "/", nbatches, "] --- fval = ", fval[end])

	    end

	    # Print msg
	    verbose && println("Epoch: [", n, "/", nepochs, "] --- fval = ", fval[end])

	end

	return fval

end


text:
# Implements invariants that combine other invariants
# - `AllInvariant` passes only if all child invariants are true
# - `AnyInvariant` passes only if any one child invariant is true
# - `SequenceInvariant` runs multiple invariants in an order and fails
#    early if one fails

abstract type InvariantList <: AbstractInvariant end

title(invs::InvariantList) = invs.title
description(invs::InvariantList) = invs.description

AbstractTrees.children(invs::InvariantList) = invs.invariants

# ## AllInvariant

struct AllInvariant{I<:AbstractInvariant} <: InvariantList
    invariants::Vector{I}
    title::String
    description::Union{Nothing, String}
    shortcircuit::Bool
end

AllInvariant(
    invariants, title::String;
    description = nothing,
    shortcircuit = true,
    kwargs...
) = invariant(AllInvariant(invariants, title, description, shortcircuit); kwargs...)

# TODO: check for errors thrown and check invariants individually
# TODO: remove SequenceInvariant and add `shortcircuit` kwarg to `AllInvariant`

function satisfies(invs::AllInvariant, input)
    results = []
    keepchecking = true
    for inv in invs.invariants
        if !keepchecking
            push!(results, missing)
            continue
        else
            res = try
                satisfies(inv, input)
            catch e
                "Unexpected error while checking invariant: $e"
            end
            push!(results, res)
            if !isnothing(res) && invs.shortcircuit
                keepchecking = false
            end
        end

    end
    return all(isnothing, results) ? nothing : results
end

function errormessage(io::IO, invs::AllInvariant, msgs)
    __combinator_errormessage(io, invs, msgs, map(__getmarker, msgs),
        faint("All invariants listed below should be satisfied:\n\n "))
end

__getmarker(::Nothing) = PASS
__getmarker(::Missing) = UNKNOWN
__getmarker(_) = FAIL

# ## AnyInvariant

struct AnyInvariant{I<:AbstractInvariant} <: InvariantList
    invariants::Vector{I}
    title::String
    description::Union{Nothing, String}
end

AnyInvariant(
    invariants, title::String;
    description = nothing,
    kwargs...
) = invariant(AnyInvariant(invariants, title, description); kwargs...)

function satisfies(invs::AnyInvariant, input)
    results = []

    for inv in invs.invariants
        res = satisfies(inv, input)
        push!(results, res)
        if isnothing(res)
            return nothing
        end
    end
    return results
end

function errormessage(io::IO, invs::AnyInvariant, msgs)
    __combinator_errormessage(io, invs, msgs, map(m -> isnothing(m) ? PASS : FAIL, msgs),
        faint("At least one of the invariants listed below should be satisfied:\n\n"))
end

# ## Helpers for printing children

const PASS = "\e[32m✔\e[0m\e[2m"
const FAIL = "\e[31m⨯\e[0m"
const UNKNOWN = "\e[33m?\e[0m"

function __combinator_errormessage(io::IO, invs, msgs, markers, msg)
    showdescription(io, invs)
    print(io, msg)
    for (inv, message, marker) in zip(invs.invariants, msgs, markers)
        __errormessage_child(io, inv; marker, message)
    end
end

function __errormessage_child(
        io,
        inv;
        marker = "o",
        message=nothing,
        indent = "    "
        )
    print(io, marker, " ")
    showtitle(io, inv)
    print(io, "\e[0m")
    println(io)

    if !(isnothing(message) || ismissing(message))
        ioi = WrapIO(io; indent)
        errormessage(ioi, inv, message)
    end
end

@testset "Combinators" begin
    i = exampleinvariant()
    @test invariant([i, i], "test") isa AllInvariant
    @test invariant([i, i], "test", :all) isa AllInvariant
    @test invariant([i, i], "test", :any) isa AnyInvariant
end


text:
using Images
using Colors
using HDF5
using FixedPointNumbers

function safe(name)
	result = replace(name, ":", " ")
	result = replace(result, "\\", " ")
	result = replace(result, "/", " ")
	result = replace(result, "*", " ")
	result = replace(result, "?", " ")
	result = replace(result, "|", " ")
	result = replace(result, "<", " ")
	result = replace(result, ">", " ")
	result = replace(result, "\"", " ")
	result = replace(result, r"[^\x20-\x7e]", "")
	result = strip(result)
	return result
end

function urlsafe(title)
	result = replace(title, r"[^\x20-\x7e]", "")
	result = strip(result)
	result = replace(result, " ", "+")
	result = replace(result, ",", "")
	result = replace(result, "!", "")
	result = replace(result, "'", "")
	result = replace(result, "-", "")
	return result
end

function process(list, movies, X, name)
	@printf "Start %s\n" name
	numImages = size(list)[1] * 3

	h5 = h5open("$(name).hdf5", "w")
	dset_data = d_create(h5, "data", datatype(Float32), dataspace(90, 120, 3, numImages))
	dset_label = d_create(h5, "label", datatype(Float32), dataspace(size(X)[2], numImages))

	idx = 0
	for i in list
		if idx % 1000 == 0
			@printf "%i/%i\t%s\n" idx numImages now()
		end

		index = findfirst(movies[:, 1] .== i)
		title = movies[index, 2]

		if isfile("movies/$(safe(title))/1.jpg") &&
			isfile("movies/$(safe(title))/2.jpg") &&
			isfile("movies/$(safe(title))/3.jpg")
			for j = 1:3
				img = load("movies/$(safe(title))/$(j).jpg")
				data = convert(Array, separate(img))
				if size(data)[1] < 90 || size(data)[2] < 120 || size(data)[3] < 3
					@show title
					continue
				else
					data = convert(Array, separate(img))[1:90, 1:120, 1:3]
				end
				#data = Array{UInt8}(round(data * 255))
				idx += 1
				h5["label"][1:size(X)[2], idx] = X[index, :]
				h5["data"][map(x -> 1:x, size(data))..., idx] = data
			end
		end
	end

	close(h5)
	@printf "End %s\n" name
end

function main(dataPath)
	lines = parse(readall(pipeline(`cat $[dataPath]movies.csv`, `wc -l`)))
	@time movies = readdlm(dataPath * "movies.csv", ',', header=true, dims=(lines, 3))
	movies = movies[1]

	lines = parse(readall(pipeline(`cat X-results.csv`, `wc -l`)))
	@time X = readdlm("X-results.csv", ',', header=false, dims=(lines, 18))

	function exclude(x)
		index = x
		title = movies[index, 2]
		if isfile("movies/$(safe(title))/1.jpg") &&
			isfile("movies/$(safe(title))/2.jpg") &&
			isfile("movies/$(safe(title))/3.jpg")
			return true
		else
			return false
		end
	end

	validMovies = filter(exclude, movies[:, 1])
	order = shuffle(collect(1:size(validMovies)[1]))
	validMovies = validMovies[order, :]

	trainEnd = Int(floor(size(validMovies)[1] * 0.8))
	trainMovies = validMovies[1:trainEnd, :]
	crossEnd = Int(floor(size(validMovies)[1] * 0.9))
	crossMovies = validMovies[trainEnd+1 : crossEnd, :]
	testMovies = validMovies[crossEnd+1 : end, :]

	@printf "Finished loading\n"

	process(trainMovies, movies, X, "train")
	process(crossMovies, movies, X, "cross")
	process(testMovies, movies, X, "test")

	@printf "Done :)\n"
end

main("../ml-latest/fixed/")


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export AccuracyLayer
type AccuracyNeuron <: Neuron
    value    :: Float32
    inputs   :: Vector{Float32}
    AccuracyNeuron() = new(0.0, [])
end

type AccuracyEnsemble <: NormalizationEnsemble
    name         :: Symbol
    neurons      :: Vector{AccuracyNeuron}
    connections  :: Vector{Connection}
    phase        :: Phase
    net_subgroup :: Cint
end

function AccuracyEnsemble(net::Net, name::Symbol)
    ens = AccuracyEnsemble(name, [], [], Test, 1)
    add_ensemble(net, ens)
    ens
end

function get_forward_args(ens::AccuracyEnsemble)
    return [symbol(ens.name,:value)]
end

function get_backward_args(ens::AccuracyEnsemble)
    return []
end

function init(ensemble::AccuracyEnsemble, net::Net)
    set_buffer(net, symbol(ensemble.name,:value), Array(Float32, 1))
end

function AccuracyLayer(name::Symbol, net::Net, input_ensemble::AbstractEnsemble,
                       label_ensemble::AbstractEnsemble)
    accuracy = AccuracyEnsemble(net, name)
    add_connections(net, input_ensemble, accuracy, (i) -> (1:length(input_ensemble),))
    add_connections(net, label_ensemble, accuracy, (i) -> (1:1,))
    accuracy
end

@output AccuracyEnsemble function forward(acc::Array, input::Array, label::Array)
    accuracy = 0.0
    for n in 1:size(input, 2)
        max_idx = 1
        max_val = input[1, n]
        for i in 2:size(input, 1)
            if input[i, n] > max_val
                max_val = input[i, n]
                max_idx = i
            end
        end
        if max_idx == round(Int, label[1, n]) + 1
            accuracy += 1.0
        end
    end
    acc[1] = accuracy / size(input, 2)
    return 0
end


text:
using DatasetManager
using Test

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


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module Flatbuf

using ..FlatBuffers

include("Schema.jl")
include("File.jl")
include("Message.jl")

end # module

max_forks_repo_name: XpressAI/arrow


text:
function custom_gate_cx(args, outs...; qasm_cgc=nothing, N=nothing)
    if length(outs) == 2
        append!(qasm_cgc, [circuit_gate((outs[2]), X, (outs[1]))])
    else
        append!(qasm_cgc, circuit_gate((outs[end]), X, (outs[1:end-1]...,)))
    end
end

function custom_gate_cy(args, outs...; qasm_cgc=nothing, N=nothing)
    if length(outs) == 2
        append!(qasm_cgc, [circuit_gate((outs[2]), Y, (outs[1]))])
    else
        append!(qasm_cgc, circuit_gate((outs[end]), Y, (outs[1:end-1]...,)))
    end
end

function custom_gate_cz(args, outs...; qasm_cgc=nothing, N=nothing)
    if length(outs) == 2
        append!(qasm_cgc, [circuit_gate((outs[2]), Z, (outs[1]))])
    else
        append!(qasm_cgc, circuit_gate((outs[end]), Z, (outs[1:end-1]...,)))
    end
end

function custom_gate_u(args, outs...; qasm_cgc=nothing, N=nothing)
    if length(outs) == 1
        append!(qasm_cgc, [circuit_gate(outs[1], RzGate(args[1])),
                  circuit_gate(outs[1], RyGate(args[2])),
                  circuit_gate(outs[1], RzGate(args[3]))
                ])
    else
        append!(qasm_cgc, [circuit_gate(outs[1], RzGate(args[1], ccntrl)),
                  circuit_gate(outs[1], RyGate(args[2]), ccntrl),
                  circuit_gate(outs[1], RzGate(args[3]), ccntrl)
                ])
    end
end


text:
using MappedArrays
using Base.Test

a = [1,4,9,16]

b = mappedarray(sqrt, a)
@test parent(b) === a
@test eltype(b) == Float64
@test @inferred(getindex(b, 1)) == 1
@test b[2] == 2
@test b[3] == 3
@test b[4] == 4
@test_throws ErrorException b[3] = 0
@test isa(eachindex(b), AbstractUnitRange)

c = mappedarray((sqrt, x->x*x), a)
@test parent(c) === a
@test @inferred(getindex(c, 1)) == 1
@test c[2] == 2
@test c[3] == 3
@test c[4] == 4
c[3] = 2
@test a[3] == 4
@test_throws InexactError c[3] = 2.2  # because the backing array is Array{Int}
@test isa(eachindex(c), AbstractUnitRange)

sb = similar(b)
@test isa(sb, Array{Float64})
@test size(sb) == size(b)


text:
"""
    SqExponential(lnℓ::AbstractFloat) <: RadialKernel{SqEuclidean}

The squared exponential kernel (aka the Gaussian kernel) with length scale
``\\exp(2 \\ln\\ell) = \\ell^2 > 0``.

```math
k_{\\ell}(x, y) = k_{\\ell}(\\lVert x - y\\rVert_2^2) =
\\exp\\bigg\\lbrace -\\frac{\\lVert x - y\\rVert_2^2}{2\\ell^2}\\bigg\\rbrace
```

External links
* [RBF Kernel on Wikipedia](https://en.wikipedia.org/wiki/Radial_basis_function_kernel)

"""
mutable struct SqExponential{F<:AbstractFloat} <: RadialKernel{SqEuclidean}
    dist::SqEuclidean
    lnℓ::F
    twoℓ²::F
    function SqExponential(dist::SqEuclidean, lnℓ::AbstractFloat)
        new{typeof(lnℓ)}(dist, lnℓ, exp(2 * lnℓ + float(logtwo)))
    end
end

SqExponential(lnℓ::AbstractFloat) = SqExponential(SqEuclidean(), lnℓ)

@inline (k::SqExponential)(τ::AbstractFloat) = exp(-τ / k.twoℓ²)

Gaussian = SqExponential

numparams(::SqExponential) = (1,)
paramtypes(::SqExponential{F}) where F = (F,)
params(k::SqExponential) = (lnℓ = k.lnℓ,)
function setparams!(k::SqExponential{F}, lnℓ::F) where F
    k.lnℓ = lnℓ
    k.twoℓ² = exp(2 * lnℓ + float(logtwo))
end


text:
# Assuming dev, net, feats, ndeps are predefined
eqparse(a,b)=(isequal(a[1],b[1]), isequal(sortcols(vcat(a[2],a[3])), sortcols(vcat(b[2],b[3]))))
gnet = copy(net, :gpu)
gnet[end].f = KUnet.logp

isdefined(:g0) || @date g0=KUparser.gparse(dev,gnet,feats,ndeps)
isdefined(:g1) || @date g1=KUparser.gparse(dev,gnet,feats,ndeps,1700)
@show eqparse(g1,g0)
isdefined(:g2) || @date g2=KUparser.gparse(dev,gnet,feats,ndeps,1700,20)
@show eqparse(g2,g0)
isdefined(:g3) || @date g3=KUparser.gparse(dev,gnet,feats,ndeps,100)
@show eqparse(g3,g0)
isdefined(:g4) || @date g4=KUparser.gparse(dev,gnet,feats,ndeps,100,20)
@show eqparse(g4,g0)
isdefined(:g100) || @date g100=KUparser.gparse(dev[1:100],gnet,feats,ndeps,100,20)

if false

isdefined(:b0) || @date b0=KUparser.bparse(dev,gnet,feats,ndeps,1)
@show eqparse(b0,g0)
isdefined(:b1) || @date b1=KUparser.bparse(dev,gnet,feats,ndeps,1,1700)
@show eqparse(b1,g0)
isdefined(:b2) || @date b2=KUparser.bparse(dev,gnet,feats,ndeps,1,1700,20)
@show eqparse(b2,g0)
isdefined(:b3) || @date b3=KUparser.bparse(dev,gnet,feats,ndeps,1,100)
@show eqparse(b3,g0)
isdefined(:b4) || @date b4=KUparser.bparse(dev,gnet,feats,ndeps,1,100,20)
@show eqparse(b4,g0)

isdefined(:b01) || @date b01=KUparser.bparse(dev[1:100],gnet,feats,ndeps,10)
@show eqparse(b01,g100)
isdefined(:b11) || @date b11=KUparser.bparse(dev[1:100],gnet,feats,ndeps,10,100)
@show eqparse(b11,b01)
isdefined(:b51) || @date b51=KUparser.bparse(dev[1:100],gnet,feats,ndeps,10,5)
@show eqparse(b51,b01)
isdefined(:b21) || @date b21=KUparser.bparse(dev[1:100],gnet,feats,ndeps,10,5,20)
@show eqparse(b21,b01)
isdefined(:b31) || @date b31=KUparser.bparse(dev[1:100],gnet,feats,ndeps,10,1)
@show eqparse(b31,b01)
isdefined(:b41) || @date b41=KUparser.bparse(dev[1:100],gnet,feats,ndeps,10,1,20)
@show eqparse(b41,b01)

isdefined(:b02) || @date b02=KUparser.bparse(dev[1:100],gnet,feats,ndeps,100)
@show eqparse(b02,g100)
@show eqparse(b02,b01)
isdefined(:b12) || @date b12=KUparser.bparse(dev[1:100],gnet,feats,ndeps,100,100)
@show eqparse(b12,b02)
isdefined(:b52) || @date b52=KUparser.bparse(dev[1:100],gnet,feats,ndeps,100,5)
@show eqparse(b52,b02)
isdefined(:b22) || @date b22=KUparser.bparse(dev[1:100],gnet,feats,ndeps,100,5,20)
@show eqparse(b22,b02)
isdefined(:b32) || @date b32=KUparser.bparse(dev[1:100],gnet,feats,ndeps,100,1)
@show eqparse(b32,b02)
isdefined(:b42) || @date b42=KUparser.bparse(dev[1:100],gnet,feats,ndeps,100,1,20)
@show eqparse(b42,b02)

end


text:
const badevs = Ref{Int}(0)
using Printf
function schurtest(A::Matrix{T}, tol; normal=false) where {T<:Complex}
    n = size(A,1)
    ulp = eps(real(T))
    if real(T) <: BlasFloat
        S = GenericSchur.gschur(A)
    else
        S = schur(A)
    end
    # test 1: S.T is upper triangular
    @test all(tril(S.T,-1) .== 0)
    # test 2: norm(A - S.Z * S.T * S.Z') / (n * norm(A) * ulp) < tol
    @test norm(A - S.Z * S.T * S.Z') / (n * norm(A) * ulp) < tol
    # test 3: S.Z is unitary: norm(I - S.Z * S.Z') / (n * ulp) < tol
    @test norm(I - S.Z * S.Z') / (n * ulp) < tol
    # test 4: S.values are e.v. of T
    @test all(csort(S.values) .== csort(diag(S.T)))

    # It is tempting to check eigenvalues against LAPACK eigvals(A)
    # for suitable types, but the comparison is misleading without
    # detailed analysis because many of the test matrices are severly
    # non-normal.

    if normal
        # verify that Schur diagonalizes normal matrices
        @test norm(triu(S.T,1)) / (n * norm(A) * ulp) < tol
    end
    # CHECKME: a few hard cases fail; is overflow protection sufficient?
    fom = norm(A) / (floatmin(real(T))/eps(real(T)))
    FOM = norm(A) / (floatmax(real(T))*eps(real(T)))
    if (fom > 25) && (FOM < 0.04)
        As = A
        Ss = S
        vtol = tol
    else
        vtol = (T == Complex{Float16} ? 200 : 20)
        # test work-arounds
        if fom <= 25
            As = A * (floatmax(real(T))*eps(real(T)))
            Ss = GenericSchur.gschur(As)
        elseif FOM >= 0.04
            As =  A * (floatmin(real(T))/eps(real(T)))
            Ss = GenericSchur.gschur(As)
        end
    end
    VR = eigvecs(Ss)
    @test norm(As * VR - VR * diagm(0 => Ss.values)) / (n * norm(As) * ulp) < vtol
    VL = eigvecs(Ss, left=true)
    @test norm(As' * VL - VL * diagm(0 => conj.(Ss.values))) / (n * norm(As) * ulp) < vtol
    if normal
        # check orthonormality of vectors where appropriate
        @test norm(VR*VR'-I) / (n * ulp) < vtol
        @test norm(VL*VL'-I) / (n * ulp) < vtol
    end
end

"""
generate a random unitary matrix

Just something good enough for present purposes.
Unitary matrices are normal, so Schur decomposition is diagonal for them.
(This is not only another test, but sometimes actually useful.)
"""
function randu(::Type{T},n) where {T<:Complex}
    if T ∈ [Float16,Float32,Float64]
        A = randn(T,n,n)
        F = qr(A)
        return Matrix(F.Q)
    else
        # don't have normal deviates for other types, but need appropriate QR
        A = randn(ComplexF64,n,n)
        F = qr(convert.(T,A))
        return Matrix(F.Q)
    end
end

Random.seed!(1234)

for T in [Complex{BigFloat},Complex{Float16}]
@testset "group $T" begin
unfl = floatmin(real(T))
ovfl = one(real(T)) / unfl
ulp = eps(real(T))
ulpinv = one(real(T)) / ulp
rtulp = sqrt(ulp)
rtulpi = one(real(T)) / rtulp

ens = [4,32]
@testset "random general" begin
    for n in ens
        A = rand(T,n,n)
        schurtest(A,10)
    end
end
@testset "random unitary" begin
    for n in ens
        A = randu(T,n)
        schurtest(A,10,normal=true)
    end
end

end # group testset
end # type loop (non-BlasFloat)

@testset "Godunov (complex)" begin
    setprecision(BigFloat, 80) do
        A,v = godunov(Complex{BigFloat})
        S = schur(A)
        @test isapprox(csort(S.values),v,atol=2e-5)
        vnew = eigvals(A)
        @test isapprox(csort(vnew),v,atol=2e-5)
    end
end

for T in [ComplexF64, ComplexF32]
@testset "group $T" begin

unfl = floatmin(real(T))
ovfl = one(real(T)) / unfl
ulp = eps(real(T))
ulpinv = one(real(T)) / ulp
rtulp = sqrt(ulp)
rtulpi = one(real(T)) / rtulp

# sizes of matrices to test
ens = [4,32]

itypemax = 10
# On author's dev system, tests pass with tol=5 for most and 6 for all,
# but generalization is risky.
# (LAPACK uses 20, FWIW.)
tols = fill(20,itypemax)

simrconds = [1.0, rtulpi, 0.0]
magns = [1.0, ovfl*ulp, unfl*ulpinv]

# test matrix from zdrves:
# conds=0 means index 3 for simrconds (i.e. irrelevant)
# type magn mode conds
# 1 1 0 0
# 2 1 0 0
# 3 1 0 0
# 4 1 4 0
# 4 1 3 0
# 4 1 1 0
# 4 2 4 0
# 4 3 4 0
# 6 1 4 1
# 6 1 3 1
# 6 1 1 1
# 6 1 5 1
# 6 1 4 2
# 6 1 3 2
# 6 1 1 2
# 6 1 5 2
# 6 2 5 2
# 6 3 5 2
# 9 1 4 0
# 9 2 3 0
# 9 3 1 0

# Other stuff in zdrves: pointless or not even used
# 1- zero
# 2- identity
# 4- diagonal, ev specified
# 5- symmetric, ev specified, not used in zdrves

@testset "Jordan" begin
    itype=3
    kmagn = [1]
    nj = length(kmagn)
    for n in ens
        for j in 1:nj
            anorm = magns[kmagn[j]]
            A = diagm(1 => ones(T,n-1)) + anorm * I
            schurtest(A,tols[itype])
        end
    end
end

@testset "general, ev specified" begin
    itype=6
    kmagn =  [1,1,1,1,1,1,1,1,2,3]
    kmode =  [4,3,1,5,4,3,1,5,5,5]
    kconds = [1,1,1,1,2,2,2,2,2,2]
    nj = length(kmode)
    for n in ens
        for j in 1:nj
#               CALL ZLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE,
#     $                      'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
#     $                      A, LDA, WORK( 2*N+1 ), IINFO )
            A = zeros(T,n,n)
            simrcond = simrconds[kconds[j]]
            imode = kmode[j]
            anorm = magns[kmagn[j]]
            rcond = ulpinv
            verbosity[] > 1 &&
                println("latme: n=$n $anorm, $imode, $rcond, $simrcond")
            latme!(A,anorm,imode,rcond,simrcond)
            schurtest(A,tols[itype])
        end
    end
end

@testset "diagonal, random ev" begin
    itype=7
    kmode = [6]
    kmagn = [1]
    nj = length(kmagn)
    for n in ens
        for j in 1:nj
     #           CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
            A = zeros(T,n,n)
            imode = kmode[j]
            anorm = magns[kmagn[j]]
            rcond = 1.0
            latmr!(A,anorm,imode,rcond,kl=0,ku=0)
            schurtest(A,tols[itype])
        end
    end
end

@testset "symmetric, random ev" begin
    itype=8
    kmagn = [1]
    kmode = [6]
    nj = length(kmagn)
    for n in ens
        for j in 1:nj
     # note sym='H'
     #           CALL ZLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
            A = zeros(T,n,n)
            imode = kmode[j]
            anorm = magns[kmagn[j]]
            rcond = 1.0
            latmr!(A,anorm,imode,rcond,sym='H',kl=n,ku=n)
            schurtest(A,tols[itype])
        end
    end
end

@testset "general, random ev" begin
    itype=9
    kmagn = [1,2,3]
    # kmode = [4,3,1] # data statement is a lie
    kmode = [6,6,6]
    nj = length(kmagn)
    for n in ens
        for j in 1:nj
     #           CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
            # zdrves zeros first 2 rows/cols & last row w/o explanation
            # presumably to test balancing, which we don't do, so skip it.
            A = zeros(T,n,n)
            imode = kmode[j]
            anorm = magns[kmagn[j]]
            rcond = 1.0
            (verbosity[] > 1) &&
                println("latmr: n=$n $anorm, $imode, $rcond")
            latmr!(A,anorm,imode,rcond)
            schurtest(A,tols[itype])
        end
    end
end

@testset "triangular, random ev" begin
    # trivial, not even run in zdrves
    itype=10
    kconds = [1]
    kmagn = [1]
    kmode = [1]
    nj = length(kmagn)
    for n in ens
        for j in 1:nj
     # note: kl=0
     #           CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
#
            A = zeros(T,n,n)
            imode = 6
            anorm = magns[kmagn[j]]
            rcond = 1.0
            latmr!(A,anorm,imode,rcond,kl=0)
            schurtest(A,tols[itype])
        end
    end
end

end # group testset
end # type loop (BlasFloat)


text:
## from base/boot.jl:
#
# immutable UTF8String <: String
#     data::Array{Uint8,1}
# end
#

## basic UTF-8 decoding & iteration ##

const utf8_offset = [
    0x00000000, 0x00003080,
    0x000e2080, 0x03c82080,
    0xfa082080, 0x82082080,
]

const utf8_trailing = [
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, 3,3,3,3,3,3,3,3,4,4,4,4,5,5,5,5,
]

is_utf8_start(byte::Uint8) = ((byte&0xc0)!=0x80)

## required core functionality ##

function endof(s::UTF8String)
    d = s.data
    i = length(d)
    i == 0 && return i
    while !is_utf8_start(d[i])
        i -= 1
    end
    i
end
length(s::UTF8String) = int(ccall(:u8_strlen, Csize_t, (Ptr{Uint8},), s.data))

function getindex(s::UTF8String, i::Int)
    # potentially faster version
    # d = s.data
    # a::Uint32 = d[i]
    # if a < 0x80; return char(a); end
    # #if a&0xc0==0x80; return '\ufffd'; end
    # b::Uint32 = a<<6 + d[i+1]
    # if a < 0xe0; return char(b - 0x00003080); end
    # c::Uint32 = b<<6 + d[i+2]
    # if a < 0xf0; return char(c - 0x000e2080); end
    # return char(c<<6 + d[i+3] - 0x03c82080)

    d = s.data
    b = d[i]
    if !is_utf8_start(b)
        j = i-1
        while 0 < j && !is_utf8_start(d[j])
            j -= 1
        end
        if 0 < j && i <= j+utf8_trailing[d[j]+1] <= length(d)
            # b is a continuation byte of a valid UTF-8 character
            error("invalid UTF-8 character index")
        end
        return '\ufffd'
    end
    trailing = utf8_trailing[b+1]
    if length(d) < i + trailing
        return '\ufffd'
    end
    c::Uint32 = 0
    for j = 1:trailing+1
        c <<= 6
        c += d[i]
        i += 1
    end
    c -= utf8_offset[trailing+1]
    char(c)
end

# this is a trick to allow inlining and tuple elision
next(s::UTF8String, i::Int) = (s[i], i+1+utf8_trailing[s.data[i]+1])

function first_utf8_byte(c::Char)
    c < 0x80    ? uint8(c)            :
    c < 0x800   ? uint8((c>>6 )|0xc0) :
    c < 0x10000 ? uint8((c>>12)|0xe0) :
                  uint8((c>>18)|0xf0)
end

## overload methods for efficiency ##

sizeof(s::UTF8String) = sizeof(s.data)

isvalid(s::UTF8String, i::Integer) =
    (1 <= i <= endof(s.data)) && is_utf8_start(s.data[i])

const empty_utf8 = UTF8String(Uint8[])

function getindex(s::UTF8String, r::UnitRange{Int})
    isempty(r) && return empty_utf8
    i, j = first(r), last(r)
    d = s.data
    if !is_utf8_start(d[i])
        i = nextind(s,i)
    end
    if j > endof(s)
        throw(BoundsError())
    end
    j = nextind(s,j)-1
    UTF8String(d[i:j])
end

function search(s::UTF8String, c::Char, i::Integer)
    if c < 0x80 return search(s.data, uint8(c), i) end
    while true
        i = search(s.data, first_utf8_byte(c), i)
        if i==0 || s[i]==c return i end
        i = next(s,i)[2]
    end
end

function rsearch(s::UTF8String, c::Char, i::Integer)
    if c < 0x80 return rsearch(s.data, uint8(c), i) end
    b = first_utf8_byte(c)
    while true
        i = rsearch(s.data, b, i)
        if i==0 || s[i]==c return i end
        i = prevind(s,i)
    end
end

function string(a::ByteString...)
    if length(a) == 1
        return a[1]::UTF8String
    end
    # ^^ at least one must be UTF-8 or the ASCII-only method would get called
    data = Array(Uint8,0)
    for d in a
        append!(data,d.data)
    end
    UTF8String(data)
end

function reverse(s::UTF8String)
    out = similar(s.data)
    if ccall(:u8_reverse, Cint, (Ptr{Uint8}, Ptr{Uint8}, Csize_t),
             out, s.data, length(out)) == 1
        error("invalid UTF-8 data")
    end
    UTF8String(out)
end

## outputing UTF-8 strings ##

print(io::IO, s::UTF8String) = (write(io, s.data);nothing)
write(io::IO, s::UTF8String) = write(io, s.data)

## transcoding to UTF-8 ##

utf8(x) = convert(UTF8String, x)
convert(::Type{UTF8String}, s::UTF8String) = s
convert(::Type{UTF8String}, s::ASCIIString) = UTF8String(s.data)
convert(::Type{UTF8String}, a::Array{Uint8,1}) = is_valid_utf8(a) ? UTF8String(a) : error("invalid UTF-8 sequence")
function convert(::Type{UTF8String}, a::Array{Uint8,1}, invalids_as::String)
    l = length(a)
    idx = 1
    iscopy = false
    while idx <= l
        if is_utf8_start(a[idx])
            nextidx = idx+1+utf8_trailing[a[idx]+1]
            (nextidx <= (l+1)) && (idx = nextidx; continue)
        end
        !iscopy && (a = copy(a); iscopy = true)
        endn = idx
        while endn <= l
            is_utf8_start(a[endn]) && break
            endn += 1
        end
        (endn > idx) && (endn -= 1)
        splice!(a, idx:endn, invalids_as.data)
        l = length(a)
    end
    UTF8String(a)
end
convert(::Type{UTF8String}, s::String) = utf8(bytestring(s))

# The last case is the replacement character 0xfffd (3 bytes)
utf8sizeof(c::Char) = c < 0x80 ? 1 : c < 0x800 ? 2 : c < 0x10000 ? 3 : c < 0x110000 ? 4 : 3

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
canreflect(::CallInfo) = false
isambiguous(::CallInfo) = false
reflect(::CallInfo; optimize=true, params=current_params()) = nothing

# Call could be resolved to a singular MI
struct MICallInfo <: CallInfo
    mi::MethodInstance
    rt
end
canreflect(::MICallInfo) = true
reflect(mi::MICallInfo; optimize=true, params=current_params()) =
    reflect(mi.mi, optimize=optimize, params=params)

struct BuiltinCallInfo <: CallInfo
    types
    rt
end

struct MultiCallInfo <: CallInfo
    sig
    rt
    callinfos::Vector{CallInfo}
end
isambiguous(::MultiCallInfo) = true
Base.collect(mci::MultiCallInfo) = mci.callinfos

struct GeneratedCallInfo <: CallInfo
    sig
    rt
end

struct FailedCallInfo <: CallInfo
    sig
    rt
end

function process_invoke(::Consumer, ref::Reflection, id, c)
    @assert c.head === :invoke
    rt = ref.CI.ssavaluetypes[id]
    return Callsite(id, MICallInfo(c.args[1], rt))
end

function process_call(::Consumer, ref::Reflection, id, c)
    @assert c.head === :call

    rt = ref.CI.ssavaluetypes[id]
    sig = map(arg -> widenconst(argextype(arg, ref.CI, ref.sptypes, ref.slottypes)), c.args)

    # Look through _apply
    ok = true
    while sig[1] === typeof(Core._apply)
        new_sig = Any[types[2]]
        for t in sig[3:end]
            if !(t <: Tuple) || t isa Union
                ok = false
                break
            end
            append!(new_sig, t.parameters)
        end
        ok || break
        sig = new_sig
    end
    if !ok
        # when does this happen?
        @error "Failed to look through _apply" id c ref
        return nothing
    end

    # Filter out builtin functions and intrinsic function
    if sig[1] <: Core.Builtin || sig[1] <: Core.IntrinsicFunction
        return id, BuiltinCallInfo(sig, rt)
    end
    return id, callinfo(sig, rt, ref)
end

function callinfo(sig, rt, ref)
    methds = Base._methods_by_ftype(sig, -1, ref.world)
    (methds === false || length(methds) < 1) && return FailedCallInfo(sig, rt)
    callinfos = CallInfo[]
    for x in methds
        atypes = x[1]
        sparams = x[2]
        meth = x[3]
        if isdefined(meth, :generator) && !Base.may_invoke_generator(meth, atypes, sparams)
            push!(callinfos, GeneratedCallInfo(sig, rt))
        else
            mi = code_for_method(meth, atypes, sparams, ref.world)
            if mi === nothing
                push!(callinfos, FailedCallInfo(sig, rt))
            end
            push!(callinfos, MICallInfo(mi, rt))
        end
    end

    @assert length(callinfos) != 0
    length(callinfos) == 1 && return first(callinfos)
    return MultiCallInfo(sig, rt, callinfos)
end