text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

const OBJC_NEW_PROPERTIES = 1

# Skipping MacroDefinition: FGAPI __attribute__ ( ( visibility ( "default" ) ) )

const FG_VERSION = "0.9.2"
const FG_VERSION_MAJOR = 0
const FG_VERSION_MINOR = 9
const FG_VERSION_PATCH = 2
const FG_API_VERSION = FG_API_VERSION_CURRENT

const fg_window = Ptr{Void}
const fg_font = Ptr{Void}
const fg_chart = Ptr{Void}
const fg_image = Ptr{Void}
const fg_histogram = Ptr{Void}
const fg_plot = Ptr{Void}
const fg_surface = Ptr{Void}
const fg_vector_field = Ptr{Void}

# begin enum ANONYMOUS_1
# const ANONYMOUS_1 = UInt32
const FG_ERR_NONE = (UInt32)(0)
const FG_ERR_SIZE = (UInt32)(1001)
const FG_ERR_INVALID_TYPE = (UInt32)(1002)
const FG_ERR_INVALID_ARG = (UInt32)(1003)
const FG_ERR_GL_ERROR = (UInt32)(2001)
const FG_ERR_FREETYPE_ERROR = (UInt32)(3001)
const FG_ERR_FILE_NOT_FOUND = (UInt32)(4001)
const FG_ERR_NOT_SUPPORTED = (UInt32)(5001)
const FG_ERR_NOT_CONFIGURED = (UInt32)(5002)
const FG_ERR_FONTCONFIG_ERROR = (UInt32)(6001)
const FG_ERR_FREEIMAGE_UNKNOWN_FORMAT = (UInt32)(7001)
const FG_ERR_FREEIMAGE_BAD_ALLOC = (UInt32)(7002)
const FG_ERR_FREEIMAGE_SAVE_FAILED = (UInt32)(7003)
const FG_ERR_INTERNAL = (UInt32)(9001)
const FG_ERR_RUNTIME = (UInt32)(9002)
const FG_ERR_UNKNOWN = (UInt32)(9003)
# end enum ANONYMOUS_1

# begin enum fg_err
const fg_err = UInt32
const FG_ERR_NONE = (UInt32)(0)
const FG_ERR_SIZE = (UInt32)(1001)
const FG_ERR_INVALID_TYPE = (UInt32)(1002)
const FG_ERR_INVALID_ARG = (UInt32)(1003)
const FG_ERR_GL_ERROR = (UInt32)(2001)
const FG_ERR_FREETYPE_ERROR = (UInt32)(3001)
const FG_ERR_FILE_NOT_FOUND = (UInt32)(4001)
const FG_ERR_NOT_SUPPORTED = (UInt32)(5001)
const FG_ERR_NOT_CONFIGURED = (UInt32)(5002)
const FG_ERR_FONTCONFIG_ERROR = (UInt32)(6001)
const FG_ERR_FREEIMAGE_UNKNOWN_FORMAT = (UInt32)(7001)
const FG_ERR_FREEIMAGE_BAD_ALLOC = (UInt32)(7002)
const FG_ERR_FREEIMAGE_SAVE_FAILED = (UInt32)(7003)
const FG_ERR_INTERNAL = (UInt32)(9001)
const FG_ERR_RUNTIME = (UInt32)(9002)
const FG_ERR_UNKNOWN = (UInt32)(9003)
# end enum fg_err

# begin enum ANONYMOUS_2
# const ANONYMOUS_2 = UInt32
const FG_GRAYSCALE = (UInt32)(100)
const FG_RG = (UInt32)(200)
const FG_RGB = (UInt32)(300)
const FG_BGR = (UInt32)(301)
const FG_RGBA = (UInt32)(400)
const FG_BGRA = (UInt32)(401)
# end enum ANONYMOUS_2

# begin enum fg_channel_format
const fg_channel_format = UInt32
const FG_GRAYSCALE = (UInt32)(100)
const FG_RG = (UInt32)(200)
const FG_RGB = (UInt32)(300)
const FG_BGR = (UInt32)(301)
const FG_RGBA = (UInt32)(400)
const FG_BGRA = (UInt32)(401)
# end enum fg_channel_format

# begin enum ANONYMOUS_3
# const ANONYMOUS_3 = UInt32
const FG_CHART_2D = (UInt32)(2)
const FG_CHART_3D = (UInt32)(3)
# end enum ANONYMOUS_3

# begin enum fg_chart_type
const fg_chart_type = UInt32
const FG_CHART_2D = (UInt32)(2)
const FG_CHART_3D = (UInt32)(3)
# end enum fg_chart_type

# begin enum ANONYMOUS_4
# const ANONYMOUS_4 = UInt32
const FG_COLOR_MAP_DEFAULT = (UInt32)(0)
const FG_COLOR_MAP_SPECTRUM = (UInt32)(1)
const FG_COLOR_MAP_RAINBOW = (UInt32)(2)
const FG_COLOR_MAP_RED = (UInt32)(3)
const FG_COLOR_MAP_MOOD = (UInt32)(4)
const FG_COLOR_MAP_HEAT = (UInt32)(5)
const FG_COLOR_MAP_BLUE = (UInt32)(6)
const FG_COLOR_MAP_INFERNO = (UInt32)(7)
const FG_COLOR_MAP_MAGMA = (UInt32)(8)
const FG_COLOR_MAP_PLASMA = (UInt32)(9)
const FG_COLOR_MAP_VIRIDIS = (UInt32)(10)
# end enum ANONYMOUS_4

# begin enum fg_color_map
const fg_color_map = UInt32
const FG_COLOR_MAP_DEFAULT = (UInt32)(0)
const FG_COLOR_MAP_SPECTRUM = (UInt32)(1)
const FG_COLOR_MAP_RAINBOW = (UInt32)(2)
const FG_COLOR_MAP_RED = (UInt32)(3)
const FG_COLOR_MAP_MOOD = (UInt32)(4)
const FG_COLOR_MAP_HEAT = (UInt32)(5)
const FG_COLOR_MAP_BLUE = (UInt32)(6)
const FG_COLOR_MAP_INFERNO = (UInt32)(7)
const FG_COLOR_MAP_MAGMA = (UInt32)(8)
const FG_COLOR_MAP_PLASMA = (UInt32)(9)
const FG_COLOR_MAP_VIRIDIS = (UInt32)(10)
# end enum fg_color_map

# begin enum ANONYMOUS_5
# const ANONYMOUS_5 = UInt32
const FG_RED = (UInt32)(0x00000000ff0000ff)
const FG_GREEN = (UInt32)(16711935)
const FG_BLUE = (UInt32)(65535)
const FG_YELLOW = (UInt32)(0x00000000ffff00ff)
const FG_CYAN = (UInt32)(16777215)
const FG_MAGENTA = (UInt32)(0x00000000ff00ffff)
const FG_WHITE = (UInt32)(0x00000000ffffffff)
const FG_BLACK = (UInt32)(255)
# end enum ANONYMOUS_5

# begin enum fg_color
const fg_color = UInt32
const FG_RED = (UInt32)(0x00000000ff0000ff)
const FG_GREEN = (UInt32)(16711935)
const FG_BLUE = (UInt32)(65535)
const FG_YELLOW = (UInt32)(0x00000000ffff00ff)
const FG_CYAN = (UInt32)(16777215)
const FG_MAGENTA = (UInt32)(0x00000000ff00ffff)
const FG_WHITE = (UInt32)(0x00000000ffffffff)
const FG_BLACK = (UInt32)(255)
# end enum fg_color

# begin enum ANONYMOUS_6
# const ANONYMOUS_6 = UInt32
const FG_INT8 = (UInt32)(0)
const FG_UINT8 = (UInt32)(1)
const FG_INT32 = (UInt32)(2)
const FG_UINT32 = (UInt32)(3)
const FG_FLOAT32 = (UInt32)(4)
const FG_INT16 = (UInt32)(5)
const FG_UINT16 = (UInt32)(6)
# end enum ANONYMOUS_6

# begin enum fg_dtype
const fg_dtype = UInt32
const FG_INT8 = (UInt32)(0)
const FG_UINT8 = (UInt32)(1)
const FG_INT32 = (UInt32)(2)
const FG_UINT32 = (UInt32)(3)
const FG_FLOAT32 = (UInt32)(4)
const FG_INT16 = (UInt32)(5)
const FG_UINT16 = (UInt32)(6)
# end enum fg_dtype

# begin enum ANONYMOUS_7
# const ANONYMOUS_7 = UInt32
const FG_PLOT_LINE = (UInt32)(0)
const FG_PLOT_SCATTER = (UInt32)(1)
const FG_PLOT_SURFACE = (UInt32)(2)
# end enum ANONYMOUS_7

# begin enum fg_plot_type
const fg_plot_type = UInt32
const FG_PLOT_LINE = (UInt32)(0)
const FG_PLOT_SCATTER = (UInt32)(1)
const FG_PLOT_SURFACE = (UInt32)(2)
# end enum fg_plot_type

# begin enum ANONYMOUS_8
# const ANONYMOUS_8 = UInt32
const FG_MARKER_NONE = (UInt32)(0)
const FG_MARKER_POINT = (UInt32)(1)
const FG_MARKER_CIRCLE = (UInt32)(2)
const FG_MARKER_SQUARE = (UInt32)(3)
const FG_MARKER_TRIANGLE = (UInt32)(4)
const FG_MARKER_CROSS = (UInt32)(5)
const FG_MARKER_PLUS = (UInt32)(6)
const FG_MARKER_STAR = (UInt32)(7)
# end enum ANONYMOUS_8

# begin enum fg_marker_type
const fg_marker_type = UInt32
const FG_MARKER_NONE = (UInt32)(0)
const FG_MARKER_POINT = (UInt32)(1)
const FG_MARKER_CIRCLE = (UInt32)(2)
const FG_MARKER_SQUARE = (UInt32)(3)
const FG_MARKER_TRIANGLE = (UInt32)(4)
const FG_MARKER_CROSS = (UInt32)(5)
const FG_MARKER_PLUS = (UInt32)(6)
const FG_MARKER_STAR = (UInt32)(7)
# end enum fg_marker_type


text:
export scale_weighting, freq_weighting, weighting

function weighting(x,condition,params)
  condition == :scales ? scale_weighting(x,params[:W_m_σ],params[:W_m_c]) :
  condition == :freqs ? freq_weighting(x,params[:W_m_σ],params[:W_m_c]) :
  condition == :track ? track_weighting(x,params[:W_m_σ_t],params[:W_m_σ_ϕ],
                                        params[:W_m_σ_N],params[:W_m_c]) :
  error("No condition named $condition.")
end

function scale_weight(cort,σ,c)
  s = log.(ustrip.(uconvert.(cycoct,scales(cort))))
  W = @. c*(1 - exp(-(s - s')^2 / (σ*log(2))^2))

  W
end

function scale_weighting(cort,σ,c)
  W = scale_weight(cort,σ,c)
  function helper(x::AbstractArray{T,2}) where T
    reshape(W*vec(sum(x,2)),:,1)
  end

  function helper(x::AbstractArray{T,3}) where T
    reshape(W*vec(sum(x,(1,3))),1,:,1)
  end

  function helper(x::AbstractArray{T,1}) where T
    W*x
  end

  helper
end

function track_weighting(tracks,σ_t,σ_p,σ_N,c)
  tcs, sds, Ns = zip(axisvalues(AxisArrays.axes(tracks,Axis{:params}()))[1]...)

  N = collect(log.(ustrip.(Ns))) # prior strength
  sd = collect(sds) # prior standard deviation
  t = collect(log.(ustrip.(uconvert.(s,tcs)))) # prior time constant

  W = @. c*(1 - exp(-(sd - sd')^2 / (σ_p*log(2))^2 +
                    # the `t` term is current always zero (see settings.toml)
                    -(t - t')^2 / (σ_t*log(2))^2 +
                    -(N - N')^2 / (σ_N*log(2))^2))

  x -> W*x
end

function freq_weighting(spect,σ,c)
  f = log.(ustrip.(uconvert.(Hz,frequencies(spect))))
  W = @. c*(1 - exp(-(f - f')^2 / (σ*log(2))^2))

  x -> W*x
end

ext: jl

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
using HDF5, FileIO, Test

fn = tempname() * ".h5"

hfile = h5open(fn, "w")
hfile["A"] = 1.0
hfile["B"] = [1,2,3]
create_group(hfile, "G")
hfile["G/A"] = collect(-3:4)
create_group(hfile, "G1/G2")
hfile["G1/G2/A"] = "hello"
close(hfile);

# test loader
data = Dict("A" => 1.0, "B"=> [1,2,3], "G/A"=>collect(-3:4), "G1/G2/A"=>"hello")
@test load(fn) == data
@test load(fn, "A") == 1.0
@test load(fn, "A","B") == (1.0, [1,2,3])
@test load(fn, "G/A") == collect(-3:4)

rm(fn)

# test saver
save(fn, data)
@test load(fn) == data
@test load(fn, "A") == 1.0
fr = h5open(fn, "r")
read(fr, "A") == 1.0
close(fr)

rm(fn)


text:
module Airdata

using LinearAlgebra
using StaticArrays
using UnPack

using Flight.Plotting
using Flight.Attitude

using Flight.Geodesy
using Flight.Atmosphere
using Flight.Kinematics

using Flight.Atmosphere: T_std, p_std, ρ_std, γ

import Flight.Plotting: plots

export get_airflow_angles, get_wind_axes
export AirData

#compute airflow angles at frame c from the c-frame aerodynamic velocity
@inline function get_airflow_angles(v_wOc_c::AbstractVector{<:Real})::Tuple{Float64, Float64}
    if norm(v_wOc_c) < 1
        α = β = 0.0
    else
        α = atan(v_wOc_c[3], v_wOc_c[1])
        β = atan(v_wOc_c[2], √(v_wOc_c[1]^2 + v_wOc_c[3]^2))
    end
    return (α, β)
end

#compute c-axes to wind axes rotation from the c-frame aerodynamic velocity
@inline function get_wind_axes(v_wOc_c::AbstractVector{<:Real})
    α, β = get_airflow_angles(v_wOc_c)
    get_wind_axes(α, β)
end

#compute c-axes to wind axes rotation from the c-frame airflow angles
@inline function get_wind_axes(α::Real, β::Real)
    q_cw = Ry(-α) ∘ Rz(β)
    return q_cw
end

struct AirData
    v_ew_n::SVector{3,Float64} #wind velocity, NED axes
    v_ew_b::SVector{3,Float64} #wind-relative, airframe axes
    v_eOb_b::SVector{3,Float64} #airframe velocity vector
    v_wOb_b::SVector{3,Float64} #airframe aerodynamic velocity vector
    α_b::Float64 #airframe angle of attack
    β_b::Float64 #airframe angle of sideslip
    T::Float64 #static temperature
    p::Float64 #static pressure
    ρ::Float64 #density
    a::Float64 #speed of sound
    μ::Float64 #dynamic viscosity
    M::Float64 #Mach number
    Tt::Float64 #total temperature
    pt::Float64 #total pressure
    Δp::Float64 #impact pressure
    q::Float64 #dynamic pressure
    TAS::Float64 #true airspeed
    EAS::Float64 #equivalent airspeed
    CAS::Float64 #calibrated airspeed
end

function AirData()
    AirData(KinData(), AtmosphericData())
end

function AirData(kin::KinData, atm_sys::AtmosphericSystem)
    #the AtmosphericData constructor accepts any Geographic subtype, but it's
    #likely that ISA SL conditions and wind will be expressed in LatLon
    atm_data = AtmosphericData(atm_sys, Geographic(kin.pos.ϕ_λ, kin.pos.h_o))
    AirData(kin, atm_data)
end

function AirData(kin::KinData, atm_data::AtmosphericData)

    v_eOb_b = kin.vel.v_eOb_b
    v_ew_n = atm_data.wind.v_ew_n
    v_ew_b = kin.pos.q_nb'(v_ew_n)
    v_wOb_b = v_eOb_b - v_ew_b
    α_b, β_b = get_airflow_angles(v_wOb_b)

    @unpack T, p, ρ, a, μ = atm_data.isa_
    TAS = norm(v_wOb_b)
    M = TAS / a
    Tt = T * (1 + (γ - 1)/2 * M^2)
    pt = p * (Tt/T)^(γ/(γ-1))
    Δp = pt - p
    q = 1/2 * ρ * TAS^2

    EAS = TAS * √(ρ / ρ_std)
    CAS = √(2γ/(γ-1) * p_std/ρ_std * ( (1 + q/p_std)^((γ-1)/γ) - 1) )

    AirData(v_ew_n, v_ew_b, v_eOb_b, v_wOb_b, α_b, β_b,
            T, p, ρ, a, μ, M, Tt, pt, Δp, q, TAS, EAS, CAS)

end

function plots(t, data::AbstractVector{<:AirData}; mode, save_path, kwargs...)

    @unpack v_ew_n, v_ew_b, v_eOb_b, v_wOb_b, α_b, β_b, a, μ, ρ, TAS, EAS, CAS, M,
            p, pt, T, Tt, Δp, q = StructArray(data)

    pd = Dict{String, Plots.Plot}()

    pd["01_v_ew_n"] = thplot(t, v_ew_n;
        plot_title = "Velocity (Wind / ECEF) [NED]",
        label = ["North" "East" "Down"],
        ylabel = [L"$v_{ew}^{N} \ (m/s)$" L"$v_{ew}^{E} \ (m/s)$" L"$v_{ew}^{D} \ (m/s)$"],
        th_split = :h,
        kwargs...)

    pd["02_v_ew_b"] = thplot(t, v_ew_b;
        plot_title = "Velocity (Wind / ECEF) [Airframe]",
        ylabel = [L"$v_{ew}^{x_b} \ (m/s)$" L"$v_{ew}^{y_b} \ (m/s)$" L"$v_{ew}^{z_b} \ (m/s)$"],
        th_split = :h,
        kwargs...)

    pd["03_v_eOb_b"] = thplot(t, v_eOb_b;
        plot_title = "Velocity (Airframe / ECEF) [Airframe]",
        ylabel = [L"$v_{eb}^{x_b} \ (m/s)$" L"$v_{eb}^{y_b} \ (m/s)$" L"$v_{eb}^{z_b} \ (m/s)$"],
        th_split = :h,
        kwargs...)

    pd["04_v_wOb_b"] = thplot(t, v_wOb_b;
        plot_title = "Velocity (Airframe / Wind) [Airframe]",
        ylabel = [L"$v_{eb}^{x_b} \ (m/s)$" L"$v_{eb}^{y_b} \ (m/s)$" L"$v_{eb}^{z_b} \ (m/s)$"],
        th_split = :h,
        kwargs...)

    splt_α = thplot(t, rad2deg.(α_b);
        title = "Angle of Attack", ylabel = L"$\alpha \ (deg)$",
        label = "", kwargs...)

    splt_β = thplot(t, rad2deg.(β_b);
        title = "Angle of Sideslip", ylabel = L"$\beta \ (deg)$",
        label = "", kwargs...)

    pd["05_α_β"] = plot(splt_α, splt_β;
        plot_title = "Airflow Angles [Airframe]",
        layout = (1,2),
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    splt_a = thplot(t, a;
        title = "Speed of Sound", ylabel = L"$a \ (m/s)$",
        label = "", kwargs...)

    splt_ρ = thplot(t, ρ;
        title = "Density", ylabel = L"$\rho \ (kg/m^3)$",
        label = "", kwargs...)

    splt_μ = thplot(t, μ;
        title = "Dynamic Viscosity", ylabel = L"$\mu \ (Pa \ s)$",
        label = "", kwargs...)

    pd["06_ρ_a"] = plot(splt_ρ, splt_a, splt_μ;
        plot_title = "Freestream Properties",
        layout = (1,3),
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    splt_T = thplot(t, hcat(T, Tt);
        title = "Temperature",
        label = ["Static"  "Total"],
        ylabel = L"$T \ (K)$",
        th_split = :none, kwargs...)

    splt_p = thplot(t, hcat(p, pt)/1000;
        title = "Pressure",
        label = ["Static"  "Total"],
        ylabel = L"$p \ (kPa)$",
        th_split = :none, kwargs...)

    pd["07_T_p"] = plot(splt_T, splt_p;
        plot_title = "Freestream Properties",
        layout = (1,2),
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    splt_airspeed = thplot(t, hcat(TAS,EAS,CAS);
        title = "Airspeed",
        label = ["True" "Equivalent" "Calibrated"],
        ylabel = L"$v \ (m/s)$",
        th_split = :none, kwargs...)

    splt_Mach = thplot(t, M;
        title = "Mach", ylabel = L"M",
        label = "", kwargs...)

    splt_q = thplot(t, q/1000;
        title = "Dynamic Pressure", ylabel = L"$q \ (kPa)$",
        label = "", kwargs...)

    l3 = @layout [a{0.5w} [b; c{0.5h}]]
    pd["08_airspeed_M_q"] = plot(splt_airspeed, splt_Mach, splt_q;
        layout = l3,
        plot_title = "Freestream Properties",
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    save_plots(pd; save_path)

end

end #module


text:
function part1(input)
    start = map_cave(input)
    n = Ref(0)
    find_paths(start, n, false)
    return n[]
end

function part2(input)
    start = map_cave(input)
    n = Ref(0)
    find_paths(start, n, true)
    return n[]
end

mutable struct Cave
    name::Symbol
    small::Bool
    visited::Bool
    neighbors::Vector{Cave}
end

Cave(name) = Cave(Symbol(name), all(islowercase, collect(name)), false, Cave[])

function map_cave(input)
    caves = Dict{String, Cave}()
    for line in eachline(input)
        cave1, cave2 = get_cave.(Ref(caves), split(line, "-"))
        push!(cave1.neighbors, cave2)
        push!(cave2.neighbors, cave1)
    end
    return caves["start"]
end

get_cave(caves, name) = get!(() -> Cave(name), caves, name)

function find_paths(cave, n, allow_second_visit)
    if cave.name == :end
        n[] += 1
        return
    end
    second_visit = false
    if cave.small && cave.visited
        if !allow_second_visit || cave.name == :start
            return
        end
        second_visit = true
        allow_second_visit = false
    end
    cave.visited = true
    for neighbor in cave.neighbors
        find_paths(neighbor, n, allow_second_visit)
    end
    if !second_visit
        cave.visited = false
    end
end


text:
"""
Sensitivity vector with iterative methods
"""

"""
stsenguess(Q::MatT, ::Type{Tv} = Float64)::Vector{Tv}

Get a vector which guesses the stationary vector of CTMC.
This is used for the iterative methods solving the stationary vector.
"""

function stsenguess(Q::MatT, ::Type{Tv} = Float64)::Vector{Tv} where {Tv,MatT}
    m, n = size(Q)
    @assert m == n
    fill(Tv(0), n)
end

"""
dtcmstcheck(Q::MatT, pis::Vector{Tv}; eps = Tv(1.0e-8))
ctmcstcheck(Q::MatT, pis::Vector{Tv}; eps = Tv(1.0e-8))

Check whether a given vector pis is a stationary vector of P ro Q
"""

function dtmcstcheck(P::MatT, pis::Vector{Tv}; eps = Tv(1.0e-8)) where {Tv,MatT}
    v = P' * pis - pis
    maximum(abs.(v)) < eps
end

function ctmcstcheck(Q::MatT, pis::Vector{Tv}; eps = Tv(1.0e-8)) where {Tv,MatT}
    v = Q' * pis
    maximum(abs.(v)) < eps
end

"""
stsengs(Q::SparseCSC{Tv,Ti}, pis::Vector{Tv}, b::Vector{Tv}:
  x0::Vector{Tv}=stsenguess(Q,Tv), maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6))

Get a sensitivity vector for stationary vector of CTMC.

Parameters:
- Q: CTMC Kernal
- pis: A stationary vector, i.e., pis * Q = 0
- b: A vector. In the case of first derivative, b = pis * Qdash where Qdash is the first derivate of Q. It may change when we want to the high-order derivatives.
- x0: Initial vector for iteration
- maxiter: The maximum number of iteration. The algorithm stops when the number of iteration becomes maxiter.
- steps: The number of steps to check the convergence
- rtol: the tolerance error. When the relative errors of two successive vectors with steps attains rtol, the algorithm stops.
Return value:
A tuple of
- x: sensitivity vector
- conv: A boolean whether the algorithm converges or not
- iter: The number of iterations
- rerror: The relative error when the algorithm stops
"""

function stsengs(Q::SparseMatrixCSC{Tv,Ti}, pis::Vector{Tv}, b::Vector{Tv};
    x0::Vector{Tv}=stsenguess(Q), maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6)) where {Tv,Ti}
    stsengs(SparseCSC(Q), pis, b, x0=x0, maxiter=maxiter, steps=steps, rtol=rtol)
end

function stsengs(Q::SparseCSC{Tv,Ti}, pis::Vector{Tv}, b::Vector{Tv};
    x0::Vector{Tv}=stsenguess(Q), maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6)) where {Tv,Ti}
    m, n = size(Q)
    @assert m == n
    @assert ctmcstcheck(Q, pis)
    x = copy(x0)
    iter = 0
    conv = false
    rerror::Tv = Tv(0)
    prevx = similar(x)
    while true
        prevx .= x
        for i in 1:steps
            gsstep!(x, Q, b, alpha=-Tv(1))
            @axpy(-sum(x), pis, x)
        end
        # rerror = maximum(abs.((x - prevx) ./ x))
        rerror = maximum(abs.(x - prevx)) / maximum(x)
        iter += steps
        if rerror < rtol
            conv = true
            break
        end
        if iter >= maxiter
            break
        end
    end
    return x, conv, iter, rerror
end

"""
stsengs(Q::SparseCSC{Tv,Ti}, pis::Vector{Tv}, b::Vector{Tv}:
  x0::Vector{Tv}=stsenguess(Q,Tv), maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6))

Get a sensitivity vector for stationary vector of CTMC.
For instance, the first derivative of stationary vector with power method

  s_n = (s_{n-1} P + pis * Pdash)(I - 1 pis)

Parameters:
- P: DTMC Kernal
- pis: A stationary vector, i.e., pis * P = pis
- b: A vector. In the case of first derivative, b = pis * Pdash where Pdash is the first derivate of P. It may change when we want to the high-order derivatives.
- x0: Initial vector for iteration
- maxiter: The maximum number of iteration. The algorithm stops when the number of iteration becomes maxiter.
- steps: The number of steps to check the convergence
- rtol: the tolerance error. When the relative errors of two successive vectors with steps attains rtol, the algorithm stops.
Return value:
A tuple of
- x: sensitivity vector
- conv: A boolean whether the algorithm converges or not
- iter: The number of iterations
- rerror: The relative error when the algorithm stops
"""

function stsenpower(P::AbstractMatrix{Tv}, pis::Vector{Tv}, b::Vector{Tv};
    x0::Vector{Tv}=stsenguess(P,Tv), maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6)) where {Tv}
    m, n = size(P)
    @assert m == n
    @assert dtmcstcheck(P, pis)
    Pdash = P'
    x = copy(x0)
    iter = 0
    conv = false
    rerror::Tv = 0
    prevx = similar(x)
    while true
        prevx .= x
        for i in 1:steps
            x = Pdash * x + b
            @axpy(-sum(x), pis, x)
        end
        # rerror = maximum(abs.((x - prevx) ./ x))
        rerror = maximum(abs.(x - prevx)) / maximum(x)
        iter += steps
        if rerror < rtol
            conv = true
            break
        end
        if iter >= maxiter
            break
        end
    end
    return x, conv, iter, rerror
end


text:
# joints
function complementarity(mechanism, joint::JointConstraint{T,N,Nc};
    scaling::Bool=false) where {T,N,Nc}

    c = []
    for (i, element) in enumerate((joint.translational, joint.rotational))
        λi = joint.impulses[2][joint_impulse_index(joint, i)]
        si, γi = split_impulses(element, λi)
        push!(c, si .* γi)
    end

    return vcat(c...)
end

# contacts
complementarity(mechanism, contact::ContactConstraint;  scaling=false) = contact.impulses[2] .* contact.impulses_dual[2]

function complementarity(mechanism, contact::ContactConstraint{T,N,Nc,Cs,N½};
    scaling::Bool = false) where {T,N,Nc,Cs<:NonlinearContact{T,N},N½}
    γ = contact.impulses[2]
    s = contact.impulses_dual[2]
    return vcat(γ[1] * s[1], cone_product(γ[@SVector [2,3,4]], s[@SVector [2,3,4]]))
end

complementarityμ(mechanism, contact::ContactConstraint; scaling=false) = complementarity(mechanism, contact, scaling=scaling) - mechanism.μ * neutral_vector(contact.model)


text:
include("base.jl")

########################################
########################################
#########      INDIVIDUAL      #########
########################################
########################################
mutable struct Individual
    genotype::Array{MGene}
    # fitness::function
    model::Chain
    Individual() = new()

end

function construct_model!(ind::Individual, input_channels::Int64, output_layer::Dense)
    layers = []
    # last_output_channels = input_channels
    last_gene = MConv((40, 40, 4), (40, 40, 4))
    last_pheno = Nothing

    for g in ind.genotype
        if typeof(g) <: MBlockGene

            ph = get_phenotype(last_gene, g)
            last_pheno = Chain(ph...)
            last_gene = g

            push!(layers, last_pheno)
        else
            if typeof(g) <: MSkipConnection
                ph = get_phenotype(last_gene, g, last_pheno)
                pop!(layers)
                push!(layers, ph)
            else
                ph = get_phenotype(last_gene, g)
                push!(layers, ph)
            end
            last_pheno = ph
            last_gene = g
            # println(ph)
        end
    end

    ### TODO: placeholder
    ex_data = randn(Float32, 40, 40, 4, 1)
    dense_in_size = size(Chain(layers..., x -> reshape(x, :, size(x, 4)))(ex_data))[1]
    ###

    ind.model = Chain(layers..., x -> reshape(x, :, size(x, 4)), Dense(dense_in_size, 1, relu; initb = ones))
end

max_stars_repo_name: mahmoud-al-najar/neuroevolution

max_issues_repo_name: mahmoud-al-najar/neuroevolution

max_forks_repo_name: mahmoud-al-najar/neuroevolution


text:
using JuliaWebAPI
using Base64
using GZip

const CONTENT_DISPOSITION_TEMPLATE = "attachment; filename="
const FILE_DOWNLOAD_HDR = Dict{String,String}("Content-Type" => "application/octet-stream", "Content-Disposition" => "")

function filebytes(filename)
    open(filename, "r") do fp
        buff = Array{UInt8}(filesize(filename))
        return read!(fp, buff)
    end
end

function servefile(filename; zipped=false)
    zipped = (zipped == "true")
    FILE_DOWNLOAD_HDR["Content-Disposition"] = string(CONTENT_DISPOSITION_TEMPLATE, filename, zipped?".gz":"")
    open(filename, "r") do fp
        buff = filebytes(filename)
        zipped || return buff

        # gzip file before returning
        gzfname = tempname()
        gzfile = GZip.open(gzfname, "w")
        try
            write(gzfile, buff)
        catch ex
            rethrow(ex)
        finally
            close(gzfile)
            buff = filebytes(gzfname)
            rm(gzfname)
        end
        buff
    end
end

function listfiles()
    iob = IOBuffer()
    println(iob, """
                    Upload a file:
                    <form method="POST" enctype="multipart/form-data" action="/savefile">
                        File: <input type="file" name="upfile">

                        Name: <input type="text" name="filename">

                        <input type="submit" value="Submit">
                    </form>

                    Click on a file to download:

""")
    for fname in readdir()
        println(iob, "
$fname | zipped
")
    end
    println(iob, "
")
    String(take!(iob))
end

function savefile(; upfile=nothing, filename=nothing)
    fname = String(base64decode(filename))
    data = base64decode(upfile)
    savefile(fname, data)
end

function savefile(filename::String, upfile::Vector{UInt8})
    filename = joinpath(pwd(), String(filename))
    open(filename, "w") do f
        write(f, upfile)
    end
    listfiles()
end

const REGISTERED_APIS = [
        (listfiles, false),
        (servefile, false, FILE_DOWNLOAD_HDR),
        (savefile, false),
    ]

process(JuliaWebAPI.create_responder(REGISTERED_APIS, "tcp://127.0.0.1:9999", true, ""))


text:
"""Struct holding the prognostic spectral variables."""
struct PrognosticVariables{NF<:AbstractFloat}
    vor         ::Array{Complex{NF},3}      # Vorticity of horizontal wind field
    div         ::Array{Complex{NF},3}      # Divergence of horizontal wind field
    temp        ::Array{Complex{NF},3}      # Absolute temperature [K]
    pres_surf   ::Array{Complex{NF},2}      # Logarithm of surface pressure [log(Pa)]
    humid       ::Array{Complex{NF},3}      # Specific humidity [g/kg]
end

"""Initialize prognostic variables from rest or restart from file."""
function initial_conditions(    P::Parameters,      # Parameter struct
                                B::Boundaries,      # Boundaries struct
                                G::GeoSpectral)     # GeoSpectral struct

    @unpack initial_conditions = P

    if initial_conditions == :rest
        ProgVars = initialize_from_rest(P,B,G)
    elseif initial_conditions == :restart
        ProgVars = initialize_from_file(P,B,G)
    else
        throw(error("Incorrect initialization option, $initial_conditions given."))
    end

    return ProgVars
end

"""Initialize a PrognosticVariables struct for an atmosphere at rest. No winds,
hence zero vorticity and divergence, but temperature, pressure and humidity are
initialised """
function initialize_from_rest(  P::Parameters,
                                B::Boundaries,
                                G::GeoSpectral)

    @unpack nlev = G.geometry
    @unpack mx, nx = G.spectral

    # conversion to type NF later when creating a PrognosticVariables struct
    vor         = zeros(Complex{Float64},mx,nx,nlev)  # vorticity
    div         = zeros(Complex{Float64},mx,nx,nlev)  # divergence
    temp        = zeros(Complex{Float64},mx,nx,nlev)  # absolute Temperature
    pres_surf   = zeros(Complex{Float64},mx,nx)       # logarithm of surface pressure
    humid       = zeros(Complex{Float64},mx,nx,nlev)  # specific humidity

    initialize_temperature!(temp,P,B,G)                     # temperature from lapse rates
    pres_surf_grid = initialize_pressure!(pres_surf,P,B,G)  # pressure from temperature profile
    initialize_humidity!(humid,pres_surf_grid,P,G)          # specific humidity from pressure

    # conversion to NF happens here implicitly
    return PrognosticVariables{P.NF}(vor,div,temp,pres_surf,humid)
end

"""Initialize spectral temperature from surface absolute temperature and constant
lapse rate (troposphere) and zero lapse rate (stratosphere)."""
function initialize_temperature!(   temp::AbstractArray{Complex{NF},3}, # spectral temperature in 3D
                                    P::Parameters,                      # Parameters struct
                                    B::Boundaries{NF},                  # Boundaries struct
                                    G::GeoSpectral{NF}                  # Geospectral struct
                                    ) where {NF<:AbstractFloat}         # number format NF

    mx,nx,nlev = size(temp)
    @unpack geopot_surf = B     # spectral surface geopotential [m²/s²]

    # temp_ref:     Reference absolute T [K] at surface z = 0, constant lapse rate
    # temp_top:     Reference absolute T in the stratosphere [K], lapse rate = 0
    # lapse_rate:   Reference temperature lapse rate -dT/dz [K/km]
    # gravity:      Gravitational acceleration [m/s^2]
    # R:            Specific gas constant for dry air [J/kg/K]
    @unpack temp_ref, temp_top, lapse_rate, gravity, R = P
    @unpack n_stratosphere_levels = P

    lapse_rate_scaled = lapse_rate/gravity/1000     # Lapse rate scaled by gravity [K/m / (m²/s²)]

    temp_surf = -lapse_rate_scaled*geopot_surf      # spectral surface air temperature from orography and lapse rate
    temp_surf[1,1] += √2*temp_ref                   # adjust mean value (spectral coefficient 1,1) with temp_ref

    # Stratosphere, set the first spectral coefficient (=mean value)
    # in uppermost levels (default: k=1,2) for lapse rate = 0
    for k in 1:n_stratosphere_levels
        temp[1,1,k] = √2*temp_top
    end

    # Temperature at tropospheric levels
    @unpack σ_levels_full = G.geometry

    for k in n_stratosphere_levels+1:nlev
        for j in 1:nx
            for i in 1:mx
                temp[i,j,k] = temp_surf[i,j]*σ_levels_full[k]^(R*lapse_rate_scaled)
            end
        end
    end
end

"""Initialize the logarithm of surface pressure `logp0` consistent with temperature profile."""
function initialize_pressure!(  pres_surf::AbstractArray{Complex{NF},2},    # logarithm of surface pressure
                                P::Parameters,                              # Parameters struct
                                B::Boundaries{NF},                          # Boundaries struct
                                G::GeoSpectral{NF}                          # Geospectral struct
                                ) where {NF<:AbstractFloat}                 # number format NF

    @unpack nlon, nlat = P

    # temp_ref:     Reference absolute T [K] at surface z = 0, constant lapse rate
    # temp_top:     Reference absolute T in the stratosphere [K], lapse rate = 0
    # lapse_rate:   Reference temperature lapse rate -dT/dz [K/km]
    # gravity:      Gravitational acceleration [m/s^2]
    # R:            Specific gas constant for dry air [J/kg/K]
    # pres_ref:     Reference surface pressure [hPa]
    @unpack temp_ref, temp_top, lapse_rate, gravity, pres_ref, R = P
    @unpack geopot_surf = B                     # spectral surface geopotential
    geopot_surf_grid = gridded(geopot_surf,G)   # convert to grid-point space

    lapse_rate_scaled = lapse_rate/gravity/1000  # Lapse rate scaled by gravity [K/m / (m²/s²)]
    log_pres_ref = log(pres_ref)                # logarithm of reference surface pressure
    pres_surf_grid = zeros(nlon, nlat)          # logarithm of surface pressure by grid point

    for j in 1:nlat
        for i in 1:nlon
            pres_surf_grid[i,j] = log_pres_ref +
                log(1 - lapse_rate_scaled*geopot_surf_grid[i,j]/temp_ref)/(R*lapse_rate_scaled)
        end
    end

    pres_surf .= spectral(pres_surf_grid,G) # convert to spectral space
    spectral_truncation!(pres_surf,G)       # truncate in spectral space

    return pres_surf_grid                   # return grid for use in initialize_humidity!
end

"""Initialize specific humidity in spectral space."""
function initialize_humidity!(  humid::AbstractArray{Complex{NF},3},    # spectral specific humidity
                                pres_surf_grid::AbstractArray{NF,2},    # logarithm of surface pressure (grid space)
                                P::Parameters,                          # Parameters struct
                                G::GeoSpectral{NF}                      # Geospectral struct
                                ) where {NF<:AbstractFloat}             # number format NF

    mx,nx,nlev = size(humid)
    @unpack nlon, nlat, n_stratosphere_levels = P
    @unpack σ_levels_full = G.geometry

    # reference saturation water vapour pressure [Pa]
    # relative humidity reference [1]
    @unpack water_pres_ref, relhumid_ref = P
    humid_ref = relhumid_ref*0.622*water_pres_ref   # reference specific humidity [Pa]

    @unpack scale_height, scale_height_humid = P            # scale height [km], scale height for spec humidity [km]
    scale_height_ratio = scale_height/scale_height_humid    # ratio of scale heights [1]

    # Specific humidity at the surface (grid space)
    humid_surf_grid = humid_ref*exp.(scale_height_ratio*pres_surf_grid)

    # Convert to spectral space and truncate
    humid_surf = spectral(humid_surf_grid,G)
    spectral_truncation!(humid_surf,G)

    # stratospheric humidity zero
    fill!(view(humid,:,:,1:n_stratosphere_levels),0)

    # Specific humidity at tropospheric levels
    for k in n_stratosphere_levels+1:nlev
        for j in 1:nx
            for i in 1:mx
                humid[i,j,k] = humid_surf[i,j]*σ_levels_full[k]^scale_height_ratio
            end
        end
    end
end


text:
import URIParser: URI
import Glob: glob

# Generic constructor which will create the appropriate
# implementation based on the host platform.
function Path(str::AbstractString)
    return @unix ? PosixPath(str) : WindowsPath(str)
end

# non-standard string literal
macro p_str(path)
    Path(path)
end

cwd() = Path(pwd())
home() = Path(homedir())

#=
Path Modifiers
===============================================
The following are methods for working with and extracting
path components
=#
hasparent(path::AbstractPath) = length(parts(path)) > 1

Base.parent(path::AbstractPath) = parents(path)[end]

function parents{T<:AbstractPath}(path::T)
    if hasparent(path)
        return map(1:length(parts(path))-1) do i
            T(parts(path)[1:i])
        end
    else
        error("$(string(path)) has no parents")
    end
end

function Base.joinpath{T<:AbstractPath}(pieces::T...)
    all_parts = ASCIIString[]

    for p in pieces
        push!(all_parts, parts(p)...)
    end

    return T(tuple(all_parts...))
end

Base.basename(path::AbstractPath) = parts(path)[end]

function filename(path::AbstractPath)
    name = basename(path)
    return split(name, '.')[1]
end

function extension(path::AbstractPath)
    name = basename(path)

    tokenized = split(name, '.')
    if length(tokenized) > 1
        return tokenized[end]
    else
        return ""
    end
end

function extensions(path::AbstractPath)
    name = basename(path)

    tokenized = split(name, '.')
    if length(tokenized) > 1
        return tokenized[2:end]
    else
        return []
    end
end

exists(path::AbstractPath) = ispath(string(path))
Base.real(path::AbstractPath) = Path(realpath(string(path)))

function Base.norm{T<:AbstractPath}(path::T)
    p = parts(path)
    result = ASCIIString[]
    rem = length(p)
    count = 0
    del = 0

    while count < length(p)
        str = p[end-count]

        if str == ".."
            del += 1
        elseif str != "."
            if del == 0
                push!(result, str)
            else
                del -= 1
            end
        end

        rem -= 1
        count += 1
    end

    return T(tuple(fill("..", del)..., reverse(result)...))
end

function Base.abs(path::AbstractPath)
    result = expanduser(path)

    if isabs(result)
        return norm(result)
    else
        return norm(joinpath(cwd(), result))
    end
end

function relative{T<:AbstractPath}(path::T, start::T=T("."))
    curdir = "."
    pardir = ".."

    p = parts(abs(path))
    s = parts(abs(start))

    p == s && return curdir

    i = 0
    while i < min(length(p), length(s))
        i += 1
        if p[i] != s[i]
            i -= 1
            break
        end
    end

    pathpart = p[i+1:findlast(x -> !isempty(x), p)]
    #pathpart = join(path_arr[i+1:findlast(x -> !isempty(x), path_arr)], path_separator)
    prefix_num = findlast(x -> !isempty(x), s) - i - 1
    if prefix_num >= 0
        relpath_ = isempty(pathpart) ?
            tuple(fill(pardir, prefix_num + 1)...) :
            tuple(fill(pardir, prefix_num + 1)..., pathpart...)
    else
        relpath_ = pathpart
    end
    return isempty(relpath_) ? T(curdir) : T(relpath_)
end

function glob{T<:AbstractPath}(path::T, pattern)
    matches = glob(pattern, string(path))

    return map(
        i -> T(i),
        matches
    )
end

function uri(path::AbstractPath)
    if isempty(root(path))
        error("$(string(path)) is not an absolute path")
    end

    uri_str = "file://$(string(path))"

    return URI(uri_str)
end

#=
The following a descriptive methods for paths
built around stat
=#
Base.stat(path::AbstractPath) = Status(stat(string(path)))
Base.lstat(path::AbstractPath) = Status(lstat(string(path)))
mode(path::AbstractPath) = stat(path).mode
Base.size(path::AbstractPath) = stat(path).size
modified(path::AbstractPath) = stat(path).mtime
created(path::AbstractPath) = stat(path).ctime
Base.isdir(path::AbstractPath) = isdir(mode(path))
Base.isfile(path::AbstractPath) = isfile(mode(path))
Base.islink(path::AbstractPath) = islink(mode(path))
Base.issocket(path::AbstractPath) = issocket(mode(path))
Base.isfifo(path::AbstractPath) = issocket(mode(path))
Base.ischardev(path::AbstractPath) = ischardev(mode(path))
Base.isblockdev(path::AbstractPath) = isblockdev(mode(path))

function Base.isexecutable(path::AbstractPath)
    s = stat(path)
    usr = User()

    return isexecutable(s.mode, :ALL) || isexecutable(s.mode, :OTHER) ||
        ( usr.uid == s.user.uid && isexecutable(s.mode, :USER) ) ||
        ( usr.gid == s.group.gid && isexecutable(s.mode, :GROUP) )
end

function Base.iswritable(path::AbstractPath)
    s = stat(path)
    usr = User()

    return iswritable(s.mode, :ALL) || iswritable(s.mode, :OTHER) ||
        ( usr.uid == s.user.uid && iswritable(s.mode, :USER) ) ||
        ( usr.gid == s.group.gid && iswritable(s.mode, :GROUP) )
end

function Base.isreadable(path::AbstractPath)
    s = stat(path)
    usr = User()

    return isreadable(s.mode, :ALL) || isreadable(s.mode, :OTHER) ||
        ( usr.uid == s.user.uid && isreadable(s.mode, :USER) ) ||
        ( usr.gid == s.group.gid && isreadable(s.mode, :GROUP) )
end

function Base.ismount(path::AbstractPath)
    isdir(path) || return false
    s1 = lstat(path)
    # Symbolic links cannot be mount points
    islink(s1) && return false
    s2 = lstat(parent(path))
    # If a directory and its parent are on different devices,  then the
    # directory must be a mount point
    (s1.device != s2.device) && return true
    (s1.inode == s2.inode) && return true
    false
end

#=
Path Operations
===============================================
The following are methods for actually manipulating the
filesystem.

NOTE: Currently, we are just wrapping base julia functions,
but in the future we'll likely be handling platform specific
code in the implementation instances.
=#

Base.cd(path::AbstractPath) = cd(string(path))
function Base.cd(fn::Function, dir::AbstractPath)
    old = cwd()
    try
        cd(dir)
        fn()
   finally
        cd(old)
    end
end

function Base.mkdir(path::AbstractPath; mode=0o777, recursive=false, exist_ok=false)
    if exists(path) && !exist_ok
        error("$path already exists.")
    else
        if !hasparent(path) || exists(parent(path))
            mkdir(string(path), mode)
        elseif hasparent(path) && !exists(parent(path)) && recursive
            mkdir(parent(path); mode=mode, recursive=recursive, exist_ok=exist_ok)
        else
            error(
                "The parent of $path does not exist. " *
                "Pass recursive=true to create it."
            )
        end
    end
end

function Base.symlink(src::AbstractPath, dest::AbstractPath; exist_ok=false, overwrite=false)
    if exists(src)
        if exists(dest) && exist_ok && overwrite
            remove(dest, recursive=true)
        end

        if !exists(dest)
            symlink(string(src), string(dest))
        elseif !exist_ok
            error("$dest already exists.")
        end
    else
        error("$src is not a valid path")
    end
end

function Base.copy(src::AbstractPath, dest::AbstractPath; recursive=false, exist_ok=false, overwrite=false, symlinks=false)
    if exists(src)
        if exists(dest) && exist_ok && overwrite
            remove(dest, recursive=true)
        end

        if !exists(dest)
            if hasparent(dest) && recursive
                mkdir(parent(dest); recursive=recursive, exist_ok=true)
            end

            cp(src, dest; follow_symlinks=symlinks)
        elseif !exist_ok
            error("$dest already exists.")
        end
    else
        error("$src is not a valid path")
    end
end

function move(src::AbstractPath, dest::AbstractPath; recursive=false, exist_ok=false, overwrite=false)
    if exists(src)
        if exists(dest) && exist_ok && overwrite
            remove(dest, recursive=true)
        end

        if !exists(dest)
            # If the destination is has missing parents
            # and parents is true then we'll create the necessary parent
            # directories.
            if hasparent(dest) && recursive
                mkdir(parent(dest); recursive=recursive, exist_ok=true)
            end

            mv(string(src), string(dest))
        elseif !exist_ok
            error("$dest already exists.")
        end
    else
        error("$src is not a valid path")
    end
end

function Base.cp(src::AbstractPath, dest::AbstractPath; remove_destination::Bool=false, follow_symlinks::Bool=false)
    cp(string(src), string(dest); remove_destination=remove_destination, follow_symlinks=follow_symlinks)
end

remove(path::AbstractPath; recursive=false) = rm(string(path); recursive=recursive)
Base.touch(path::AbstractPath) = touch(string(path))

tmpname() = Path(tempname())
tmpdir() = Path(tempdir())

function mktmp(parent::AbstractPath=Path(tempdir()))
    path, io = mktemp(string(parent))
    return Path(path), io
end

mktmpdir(parent::AbstractPath=tmpdir()) = Path(mktempdir(string(parent)))

function mktmp(fn::Function, parent=tmpdir())
    (tmp_path, tmp_io) = mktmp(parent)
    try
        fn(tmp_path, tmp_io)
    finally
        close(tmp_io)
        remove(tmp_path)
    end
end

function mktmpdir(fn::Function, parent=tmpdir())
    tmpdir = mktmpdir(parent)
    try
        fn(tmpdir)
    finally
        remove(tmpdir, recursive=true)
    end
end

function chown(path::AbstractPath, user::AbstractString, group::AbstractString; recursive=false)
    @unix_only begin
        chown_cmd = ByteString["chown"]
        if recursive
            push!(chown_cmd, "-R")
        end
        append!(chown_cmd, ByteString["$(user):$(group)", string(path)])

        run(Cmd(chown_cmd))
    end

    @windows_only error("chown is currently not supported on windows.")
end

function Base.chmod(path::AbstractPath, mode::Mode; recursive=false)
    chmod_path = string(path)
    chmod_mode = raw(mode)

    if isdir(path) && recursive
        for p in glob(path, "*")
            chmod(chmod_path, chmod_mode; recursive=recursive)
        end
    end

    chmod(string(path), raw(mode))
end

function Base.chmod(path::AbstractPath; user::UInt8=0o0, group::UInt8=0o0, other::UInt8=0o0, recursive=false)
    chmod(path, Mode(user=user, group=group, other=other); recursive=recursive)
end

function Base.chmod(path::AbstractPath, symbolic_mode::AbstractString; recursive=false)
    who_char = ['u', 'g', 'o']
    who_actual = [:USER, :GROUP, :OTHER]
    act_char = ['+', '-', '=']
    perm_char = ['r', 'w', 'x']
    perm_actual = [READ, WRITE, EXEC]
    unsupported_perm_char = ['s', 't', 'X', 'u', 'g', 'o']

    tokenized = split(symbolic_mode, act_char)
    if length(tokenized) != 2
        error("Invalid symbolic string expected format <who><action><perm>.")
    end

    who_raw = tokenized[1]
    perm_raw = tokenized[2]

    who = [:ALL]
    perm = 0o0

    for i in 1:3
        if who_char[i] in who_raw
            push!(who, who_actual[i])
        end
    end

    for i in 1:3
        if perm_char[i] in perm_raw
            perm += perm_actual[i]
        end
    end

    for x in unsupported_perm_char
        if x in perm_raw
            error("$x is currently an unsupported permission char for symbolic modes.")
        end
    end

    m = mode(path)
    new_m = Mode(perm, who...)

    if '+' in symbolic_mode
        chmod(path, m + new_m; recursive=recursive)
    elseif '-' in symbolic_mode
        chmod(path, m - new_m; recursive=recursive)
    elseif '=' in symbolic_mode
        chmod(path, new_m; recursive=recursive)
    else
        error("No valid action found in symbolic mode string.")
    end
end

Base.read(path::AbstractPath) = open(readall, string(path))

function Base.write(path::AbstractPath, content::AbstractString)
    open(string(path), "w") do f
        write(f, content)
    end
end


text:
using CairoMakie, Makie.GeometryBasics
n = 20
function ngonShape(h, k, r, n)
    Polygon([Point2f(h .+ r*sin.(m*2π/n),k .+ r*cos.(m*2π/n)) for m in 1:n])
end
polysCentric = [ngonShape(0, 0, 3/i^1.5, i) for i in 3:n]
polysCircular = [ngonShape(√2/2*sin(θ), √2/2*cos(θ), 0.15/√idx, idx + 2)
    for (idx, θ) in enumerate(LinRange(0,2π*(1 -1/(n-2)), n-2))]
cmap = to_colormap(:linear_protanopic_deuteranopic_kbw_5_98_c40_n256)[3:end]

with_theme(theme_black()) do
    fig, ax, = poly(polysCentric; color = 1:n-2, colormap = cmap,
        axis = (;aspect = DataAspect()), figure = (;resolution = (600,400)))
    poly!(polysCircular; color = 1:n-2, colormap = cmap)
    hidedecorations!(ax; grid = false)
    hidespines!(ax)
    save("polygons.svg", fig)
    display(fig)
end


text:
using Individual.SchemaBase
using Individual.SchemaEvents
using Catlab.ACSetInterface
using Test

@testset "basic model construction with integer initial states works" begin

    initial_states = [1,2,3,1,2,3]
    state_labels = ["S", "I", "R"]

    SIR = SchedulingIBM{String, Int64, String, Vector{Function}}()
    initialize_states(SIR, initial_states, state_labels)

    @test nparts(SIR, :Person) == length(initial_states)
    @test nparts(SIR, :State) == length(state_labels)
    @test subpart(SIR, :statelabel) == state_labels
    @test subpart(SIR, :state) == initial_states
    @test subpart(SIR, :state_update) == zeros(Int64, length(initial_states))

    @test npeople(SIR) == length(initial_states)
    @test npeople(SIR, "S") == count(x -> x == 1, initial_states)
    @test npeople(SIR, "I") == count(x -> x == 2, initial_states)
    @test npeople(SIR, "R") == count(x -> x == 3, initial_states)
    @test nstate(SIR) == length(state_labels)
    @test statelabel(SIR) == state_labels
    @test get_index_state(SIR) == 1:length(initial_states)
    @test get_index_state(SIR, "S") == findall(x -> x == 1, initial_states)
    @test get_index_state(SIR, "I") == findall(x -> x == 2, initial_states)
    @test get_index_state(SIR, "R") == findall(x -> x == 3, initial_states)

end

@testset "basic model construction with string initial states works" begin

    initial_states = ["S", "I", "R", "S", "I", "R"]
    state_labels = ["S", "I", "R"]

    SIR = SchedulingIBM{String, Int64, String, Vector{Function}}()
    initialize_states(SIR, initial_states, state_labels)

    @test nparts(SIR, :Person) == length(initial_states)
    @test nparts(SIR, :State) == length(state_labels)
    @test subpart(SIR, :statelabel) == state_labels
    @test subpart(SIR, :state) == indexin(initial_states, state_labels)
    @test subpart(SIR, :state_update) == zeros(Int64, length(initial_states))

    @test npeople(SIR) == length(initial_states)
    @test npeople(SIR, "S") == count(x -> x == "S", initial_states)
    @test npeople(SIR, "I") == count(x -> x == "I", initial_states)
    @test npeople(SIR, "R") == count(x -> x == "R", initial_states)
    @test nstate(SIR) == length(state_labels)
    @test statelabel(SIR) == state_labels
    @test get_index_state(SIR) == 1:length(initial_states)
    @test get_index_state(SIR, "S") == findall(x -> x == "S", initial_states)
    @test get_index_state(SIR, "I") == findall(x -> x == "I", initial_states)
    @test get_index_state(SIR, "R") == findall(x -> x == "R", initial_states)

end

@testset "construct model errors properly" begin

    initial_states = [1,2,3,1,2,3,4]
    state_labels = ["S", "I", "R"]

    SIR = SchedulingIBM{String, Int64, String, Vector{Function}}()
    @test_throws ArgumentError initialize_states(SIR, initial_states, state_labels)

    initial_states = ["S", "I", "R", "S", "I", "R", "X"]
    @test_throws ArgumentError initialize_states(SIR, initial_states, state_labels)

end

@testset "events can be added to model" begin

    initial_states = ["S", "I", "R", "S", "I", "R"]
    state_labels = ["S", "I", "R"]

    SIR = SchedulingIBM{String, Int64, String, Vector{Function}}()
    initialize_states(SIR, initial_states, state_labels)

    function listener1(target, t::Int)
        queue_state_update(SIR, target, "R")
    end

    add_event(SIR, "Recovery", listener1)

    @test length(subpart(SIR, 1, :eventlistener)) == 1
    @test subpart(SIR, 1, :eventlistener) isa Vector{Function}
    @test subpart(SIR, 1, :eventlabel) == "Recovery"

    function listener2(target, t::Int)
        queue_state_update(SIR, target, "I")
    end

    add_event(SIR, "Infection", listener2)

    @test length(subpart(SIR, 2, :eventlistener)) == 1
    @test subpart(SIR, 2, :eventlistener) isa Vector{Function}
    @test subpart(SIR, 2, :eventlabel) == "Infection"

    function listener3(target, t::Int)
        queue_state_update(SIR, target[1], "I")
    end

    function listener4(target, t::Int)
        queue_state_update(SIR, target[2], "I")
    end

    add_event(SIR, "Infection1", listener3, listener4)

    @test length(subpart(SIR, 3, :eventlistener)) == 2
    @test subpart(SIR, 3, :eventlistener) isa Vector{Function}
    @test subpart(SIR, 3, :eventlabel) == "Infection1"

end

@testset "events can be scheduled, canceled, and fired" begin

    initial_states = ["S", "I", "R", "S", "I", "R"]
    state_labels = ["S", "I", "R", "D"]

    SIR = SchedulingIBM{String, Int64, String, Vector{Function}}()
    initialize_states(SIR, initial_states, state_labels)

    function listener1(target, t::Int)
        queue_state_update(SIR, target, "R")
    end

    add_event(SIR, "Recovery", listener1)

    function listener2(target, t::Int)
        queue_state_update(SIR, target, "D")
    end

    add_event(SIR, "Death", listener2)

    schedule_event(SIR, 2:4, 3, "Recovery")

    @test get_scheduled(SIR, "Recovery") == 2:4
    @test subpart(SIR, :delay) == fill(3,3)
    @test subpart(SIR, :scheduled_to_event) == fill(1,3)

    clear_schedule(SIR, 3)
    @test get_scheduled(SIR, "Recovery") == [2,4]

    schedule_event(SIR, [1,5], [2,4], "Death")

    # tick 1
    event_tick(SIR)
    event_process(SIR, 1)
    apply_state_updates(SIR)
    @test subpart(SIR, :state) == indexin(initial_states, state_labels)

    # tick 2
    event_tick(SIR)
    event_process(SIR, 1)
    apply_state_updates(SIR)
    @test subpart(SIR, :state) == indexin(["D", "I", "R", "S", "I", "R"], state_labels)

    # tick 3
    event_tick(SIR)
    event_process(SIR, 1)
    apply_state_updates(SIR)
    @test subpart(SIR, :state) == indexin(["D", "R", "R", "R", "I", "R"], state_labels)

    # tick 4
    event_tick(SIR)
    event_process(SIR, 1)
    apply_state_updates(SIR)
    @test subpart(SIR, :state) == indexin(["D", "R", "R", "R", "D", "R"], state_labels)

end

@testset "events with multiple listeners work" begin

    initial_states = ["S", "S", "S", "S", "S", "S"]
    state_labels = ["S", "I", "R"]

    SIR = SchedulingIBM{String, Int64, String, Vector{Function}}()
    initialize_states(SIR, initial_states, state_labels)

    function listener1(target, t::Int)
        queue_state_update(SIR, target[1], "I")
    end

    function listener2(target, t::Int)
        queue_state_update(SIR, target[2], "I")
    end

    add_event(SIR, "Infection", listener1, listener2)

    schedule_event(SIR, [1,5,3], [0,0,0], "Infection")

    event_process(SIR, 1)
    apply_state_updates(SIR)

    @test subpart(SIR, :state) == indexin(["I", "S", "S", "S", "I", "S"], state_labels)

end

@testset "can clear all events" begin

    initial_states = ["S", "I", "R", "S", "I", "R"]
    state_labels = ["S", "I", "R", "D"]

    SIR = SchedulingIBM{String, Int64, String, Vector{Function}}()
    initialize_states(SIR, initial_states, state_labels)

    function listener1(target, t::Int)
        queue_state_update(SIR, target, "R")
    end

    add_event(SIR, "Recovery", listener1)

    schedule_event(SIR, 2:4, 3, "Recovery")

    @test get_scheduled(SIR, "Recovery") == 2:4

    clear_schedule(SIR)

    @test length(get_scheduled(SIR, "Recovery")) == 0

end


text:
"""
    LinearMixedModel

Linear mixed-effects model representation

## Fields

* `formula`: the formula for the model
* `reterms`: a `Vector{ReMat{T}}` of random-effects terms.
* `feterms`: a `Vector{FeMat{T}}` of the fixed-effects model matrix and the response
* `sqrtwts`: vector of square roots of the case weights.  Can be empty.
* `A`: an `nt × nt` symmetric `BlockMatrix` of matrices representing `hcat(Z,X,y)'hcat(Z,X,y)`
* `L`: a `nt × nt` `BlockMatrix` - the lower Cholesky factor of `Λ'AΛ+I`
* `optsum`: an [`OptSummary`](@ref) object

## Properties

* `θ` or `theta`: the covariance parameter vector used to form λ
* `β` or `beta`: the fixed-effects coefficient vector
* `λ` or `lambda`: a vector of lower triangular matrices repeated on the diagonal blocks of `Λ`
* `σ` or `sigma`: current value of the standard deviation of the per-observation noise
* `b`: random effects on the original scale, as a vector of matrices
* `u`: random effects on the orthogonal scale, as a vector of matrices
* `lowerbd`: lower bounds on the elements of θ
* `X`: the fixed-effects model matrix
* `y`: the response vector
"""
struct LinearMixedModel{T <: AbstractFloat} <: MixedModel{T}
    formula::FormulaTerm
    reterms::Vector{ReMat{T}}
    feterms::Vector{FeMat{T}}
    sqrtwts::Vector{T}
    A::BlockMatrix{T}            # cross-product blocks
    L::BlockMatrix{T}
    optsum::OptSummary{T}
end
LinearMixedModel(f::FormulaTerm, tbl;
                 contrasts = Dict{Symbol,Any}(),
                 wts = []) =
    LinearMixedModel(f::FormulaTerm, Tables.columntable(tbl),
                     contrasts = contrasts,
                     wts = wts)
function LinearMixedModel(f::FormulaTerm, tbl::Tables.ColumnTable;
                          contrasts = Dict{Symbol,Any}(),
                          wts = [])
    # TODO: perform missing_omit() after apply_schema() when improved
    # missing support is in a StatsModels release
    tbl, _ = StatsModels.missing_omit(tbl, f)
    form = apply_schema(f, schema(f, tbl, contrasts), LinearMixedModel)
    # tbl, _ = StatsModels.missing_omit(tbl, form)

    y, Xs = modelcols(form, tbl)

    y = reshape(float(y), (:, 1)) # y as a floating-point matrix
    T = eltype(y)

    reterms = ReMat{T}[]
    feterms = FeMat{T}[]
    for (i,x) in enumerate(Xs)
        if isa(x, ReMat{T})
            push!(reterms, x)
        else
            cnames = coefnames(form.rhs[i])
            push!(feterms, FeMat(x, isa(cnames, String) ? [cnames] : collect(cnames)))
        end
    end
    push!(feterms, FeMat(y, [""]))

    # detect and combine RE terms with the same grouping var
    if length(reterms) > 1
        reterms = amalgamate(reterms)
    end

    sort!(reterms, by=nranef, rev=true)

    # create A and L
    terms = vcat(reterms, feterms)
    k = length(terms)
    sz = append!(size.(reterms, 2), rank.(feterms))
    A = BlockArray(undef_blocks, AbstractMatrix{T}, sz, sz)
    L = BlockArray(undef_blocks, AbstractMatrix{T}, sz, sz)
    for j in 1:k
        for i in j:k
            Lij = L[Block(i,j)] = densify(terms[i]'terms[j])
            A[Block(i,j)] = deepcopy(isa(Lij, BlockedSparse) ? Lij.cscmat : Lij)
        end
    end
    for i in 2:length(reterms) # check for fill-in due to non-nested grouping factors
        ci = reterms[i]
        for j in 1:(i - 1)
            cj = reterms[j]
            if !isnested(cj, ci)
                for l in i:k
                    L[Block(l, i)] = Matrix(L[Block(l, i)])
                end
                break
            end
        end
    end
    lbd = foldl(vcat, lowerbd(c) for c in reterms)
    θ = foldl(vcat, getθ(c) for c in reterms)
    optsum = OptSummary(θ, lbd, :LN_BOBYQA, ftol_rel = T(1.0e-12), ftol_abs = T(1.0e-8))
    fill!(optsum.xtol_abs, 1.0e-10)
    LinearMixedModel(form, reterms, feterms, sqrt.(convert(Vector{T}, wts)), A, L, optsum)
end
fit(::Type{LinearMixedModel}, f::FormulaTerm, tbl;
    wts = [],
    contrasts = Dict{Symbol,Any}(),
    verbose::Bool = false,
    REML::Bool = false) =
    fit(LinearMixedModel, f, Tables.columntable(tbl),
        wts = wts, contrasts = contrasts, verbose = verbose, REML = REML)
fit(::Type{LinearMixedModel},
    f::FormulaTerm,
    tbl::Tables.ColumnTable;
    wts = wts, contrasts = contrasts, verbose = verbose, REML = REML) =
    fit!(LinearMixedModel(f, tbl,
                          contrasts = contrasts,
                          wts = wts),
                          verbose = verbose,
                          REML = REML)

StatsBase.coef(m::MixedModel) = fixef(m, false)

function βs(m::LinearMixedModel)
    fetrm = first(m.feterms)
    (; (k => v for (k,v) in zip(Symbol.(fetrm.cnames), fixef(m)))...)
end

StatsBase.coefnames(m::LinearMixedModel) = first(m.feterms).cnames

function StatsBase.coeftable(m::MixedModel)
    co = coef(m)
    se = stderror(m)
    z = co ./ se
    pvalue = ccdf.(Chisq(1), abs2.(z))
    CoefTable(hcat(co, se, z, pvalue), ["Estimate", "Std.Error", "z value", "P(>|z|)"],
        first(m.feterms).cnames, 4)
end

"""
    cond(m::MixedModel)

Return a vector of condition numbers of the λ matrices for the random-effects terms
"""
LinearAlgebra.cond(m::MixedModel) = cond.(m.λ)

"""
    condVar(m::LinearMixedModel)

Return the conditional variances matrices of the random effects.

The random effects are returned by `ranef` as a vector of length `k`,
where `k` is the number of random effects terms.  The `i`th element
is a matrix of size `vᵢ × ℓᵢ`  where `vᵢ` is the size of the
vector-valued random effects for each of the `ℓᵢ` levels of the grouping
factor.  Technically those values are the modes of the conditional
distribution of the random effects given the observed data.

This function returns an array of `k` three dimensional arrays,
where the `i`th array is of size `vᵢ × vᵢ × ℓᵢ`.  These are the
diagonal blocks from the conditional variance-covariance matrix,

    s² Λ(Λ'Z'ZΛ + I)⁻¹Λ'
"""
function condVar(m::LinearMixedModel{T}) where {T}
    retrms = m.reterms
    t1 = first(retrms)
    L11 = m.L[Block(1, 1)]
    if !isone(length(retrms)) || !isa(L11, Diagonal{T, Vector{T}})
        throw(ArgumentError("code for vector-valued r.e. or more than one term not yet written"))
    end
    ll = first(t1.λ)
    Ld = L11.diag
    Array{T, 3}[reshape(abs2.(ll ./ Ld) .* varest(m), (1, 1, length(Ld)))]
end

"""
    describeblocks(io::IO, m::MixedModel)
    describeblocks(m::MixedModel)

Describe the types and sizes of the blocks in the lower triangle of `m.A` and `m.L`.
"""
function describeblocks(io::IO, m::LinearMixedModel)
    A = m.A
    L = m.L
    for i in 1:BlockArrays.nblocks(A, 2), j in 1:i
        println(io, i, ",", j, ": ", typeof(A[Block(i, j)]), " ",
            BlockArrays.blocksize(A, (i, j)), " ", typeof(L[Block(i, j)]))
    end
end
describeblocks(m::MixedModel) = describeblocks(stdout, m)

StatsBase.deviance(m::MixedModel) = objective(m)

GLM.dispersion(m::LinearMixedModel, sqr::Bool=false) = sqr ? varest(m) : sdest(m)

GLM.dispersion_parameter(m::LinearMixedModel) = true

StatsBase.dof(m::LinearMixedModel) = size(m)[2] + sum(nθ, m.reterms) + 1

function StatsBase.dof_residual(m::LinearMixedModel)::Int
    (n, p, q, k) = size(m)
    n - m.optsum.REML * p
end

"""
    feL(m::MixedModel)

Return the lower Cholesky factor for the fixed-effects parameters, as an `LowerTriangular`
`p × p` matrix.
"""
feL(m::LinearMixedModel) = LowerTriangular(m.L.blocks[end - 1, end - 1])

"""
    fit!(m::LinearMixedModel[; verbose::Bool=false, REML::Bool=false])

Optimize the objective of a `LinearMixedModel`.  When `verbose` is `true` the values of the
objective and the parameters are printed on stdout at each function evaluation.
"""
function fit!(m::LinearMixedModel{T}; verbose::Bool=false, REML::Bool=false) where {T}
    optsum = m.optsum
    opt = Opt(optsum)
    feval = 0
    optsum.REML = REML
    function obj(x, g)
        isempty(g) || error("gradient not defined")
        feval += 1
        val = objective(updateL!(setθ!(m, x)))
        feval == 1 && (optsum.finitial = val)
        verbose && println("f_", feval, ": ", round(val, digits=5), " ", x)
        val
    end
    NLopt.min_objective!(opt, obj)
    fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial))
    ## check if small non-negative parameter values can be set to zero
    xmin_ = copy(xmin)
    lb = optsum.lowerbd
    for i in eachindex(xmin_)
        if iszero(lb[i]) && zero(T) < xmin_[i] < T(0.001)
            xmin_[i] = zero(T)
        end
    end
    if xmin_ ≠ xmin
        if (zeroobj = obj(xmin_, T[])) ≤ (fmin + 1.e-5)
            fmin = zeroobj
            copyto!(xmin, xmin_)
        end
    end
    ## ensure that the parameter values saved in m are xmin
    updateL!(setθ!(m, xmin))

    optsum.feval = feval
    optsum.final = xmin
    optsum.fmin = fmin
    optsum.returnvalue = ret
    ret == :ROUNDOFF_LIMITED && @warn("NLopt was roundoff limited")
    if ret ∈ [:FAILURE, :INVALID_ARGS, :OUT_OF_MEMORY, :FORCED_STOP, :MAXFEVAL_REACHED]
        @warn("NLopt optimization failure: $ret")
    end
    m
end

function fitted!(v::AbstractArray{T}, m::LinearMixedModel{T}) where {T}
    ## FIXME: Create and use `effects(m) -> β, b` w/o calculating β twice
    vv = vec(v)
    mul!(vv, first(m.feterms), fixef(m))
    for (rt, bb) in zip(m.reterms, ranef(m))
        unscaledre!(vv, rt, bb)
    end
    v
end

StatsBase.fitted(m::LinearMixedModel{T}) where {T} = fitted!(Vector{T}(undef, nobs(m)), m)

"""
    fixef!(v::Vector{T}, m::LinearMixedModel{T})

Overwrite `v` with the pivoted and, possibly, truncated fixed-effects coefficients of model `m`
"""
fixef!(v::AbstractVector{T}, m::LinearMixedModel{T}) where {T} =
    ldiv!(feL(m)', copyto!(v, m.L.blocks[end, end - 1]))

"""
    fixef(m::MixedModel, permuted=true)

Return the fixed-effects parameter vector estimate of `m`.

If `permuted` is `true` the vector elements are permuted according to
`m.trms[end - 1].piv` and truncated to the rank of that term.
"""
function fixef(m::LinearMixedModel{T}, permuted=true) where {T}
    val = ldiv!(feL(m)', vec(copy(m.L.blocks[end, end-1])))
    if !permuted
        Xtrm = first(m.feterms)
        piv = Xtrm.piv
        p = length(piv)
        if Xtrm.rank < p
            val = copyto!(fill(-zero(T), p), val)
        end
        invpermute!(val, piv)
    end
    val
end

"""
    fnames(m::MixedModel)

Return the names of the grouping factors for the random-effects terms.
"""
fnames(m::MixedModel) = ((Symbol(tr.trm.sym) for tr in m.reterms)...,)

"""
    getθ(m::LinearMixedModel)

Return the current covariance parameter vector.
"""
getθ(m::LinearMixedModel{T}) where {T} = foldl(vcat, getθ.(m.reterms))

function getθ!(v::AbstractVector{T}, m::LinearMixedModel{T}) where {T}
    k = 0
    for t in m.reterms
        nt = nθ(t)
        getθ!(view(v, (k+1):(k+nt)), t)
        k += nt
    end
    v
end

function Base.getproperty(m::LinearMixedModel, s::Symbol)
    if s == :θ || s == :theta
        getθ(m)
    elseif s == :β || s == :beta
        coef(m)
    elseif s == :βs || s == :betas
        βs(m)
    elseif s == :λ || s == :lambda
        getproperty.(m.reterms, :λ)
    elseif s == :σ || s == :sigma
        sdest(m)
    elseif s == :σs || s == :sigmas
        σs(m)
    elseif s == :σρs || s == :sigmarhos
        σρs(m)
    elseif s == :b
        ranef(m)
    elseif s == :objective
        objective(m)
    elseif s == :pvalues
        ccdf.(Chisq(1), abs2.(coef(m) ./ stderror(m)))
    elseif s == :stderror
        stderror(m)
    elseif s == :u
        ranef(m, uscale = true)
    elseif s == :lowerbd
        m.optsum.lowerbd
    elseif s == :X
        modelmatrix(m)
    elseif s == :y
        vec(last(m.feterms).x)
    elseif s == :rePCA
        rePCA(m)
    else
        getfield(m, s)
    end
end

"""
    issingular(m::LinearMixedModel, θ=m.θ)

Test whether the model `m` is singular if the parameter vector is `θ`.
"""
issingular(m::LinearMixedModel, θ::AbstractVector=m.θ) = any(isapprox.(m.optsum.lowerbd, θ))

function StatsBase.loglikelihood(m::LinearMixedModel)
    if m.optsum.REML
        throw(ArgumentError("loglikelihood not available for models fit by REML"))
    end
    -objective(m) / 2
end

lowerbd(m::LinearMixedModel) = m.optsum.lowerbd

"""
    likelihoodratiotest(m::LinearMixedModel...)

Likeihood ratio test applied to a set of nested models.

Note that nesting of the models is not checked.  It is incumbent on the user to check this.
"""
function likelihoodratiotest(m::LinearMixedModel...)
    m = collect(m)   # change the tuple to an array
    dofs = dof.(m)
    ord = sortperm(dofs)
    dofs = dofs[ord]
    devs = deviance.(m)[ord]
    dofdiffs = diff(dofs)
    devdiffs = .-(diff(devs))
    pvals = ccdf.(Chisq.(dofdiffs), devdiffs)
    (models=(dof=dofs, deviance=devs), tests=(dofdiff=dofdiffs, deviancediff=devdiffs, p_values=pvals))
end

function StatsBase.modelmatrix(m::LinearMixedModel)
    fetrm = first(m.feterms)
    if fetrm.rank == size(fetrm, 2)
        fetrm.x
    else
        fetrm.x[:, invperm(fetrm.piv)]
    end
end

nθ(m::LinearMixedModel) = sum(nθ, m.reterms)

StatsBase.nobs(m::LinearMixedModel) = first(size(m))

"""
    objective(m::LinearMixedModel)

Return negative twice the log-likelihood of model `m`
"""
function objective(m::LinearMixedModel)
    wts = m.sqrtwts
    logdet(m) + dof_residual(m)*(1 + log2π + log(varest(m))) - (isempty(wts) ? 0 : 2sum(log, wts))
end

StatsBase.predict(m::LinearMixedModel) = fitted(m)

Base.propertynames(m::LinearMixedModel, private=false) =
    (:formula, :sqrtwts, :A, :L, :optsum, :θ, :theta, :β, :beta, :λ, :lambda, :stderror,
     :σ, :sigma, :σs, :sigmas, :b, :u, :lowerbd, :X, :y, :rePCA, :reterms, :feterms,
     :objective, :pvalues)

"""
    pwrss(m::LinearMixedModel)

The penalized, weighted residual sum-of-squares.
"""
pwrss(m::LinearMixedModel) = abs2(sqrtpwrss(m))

"""
    ranef!(v::Vector{Matrix{T}}, m::MixedModel{T}, β, uscale::Bool) where {T}

Overwrite `v` with the conditional modes of the random effects for `m`.

If `uscale` is `true` the random effects are on the spherical (i.e. `u`) scale, otherwise
on the original scale
"""
function ranef!(v::Vector, m::LinearMixedModel{T}, β::AbstractArray{T}, uscale::Bool) where {T}
    (k = length(v)) == length(m.reterms) || throw(DimensionMismatch(""))
    L = m.L
    for j in 1:k
        mulαβ!(vec(copyto!(v[j], L[Block(BlockArrays.nblocks(L, 2), j)])),
            L[Block(k + 1, j)]', β, -one(T), one(T))
    end
    for i in k: -1 :1
        Lii = L[Block(i, i)]
        vi = vec(v[i])
        ldiv!(adjoint(isa(Lii, Diagonal) ? Lii : LowerTriangular(Lii)), vi)
        for j in 1:(i - 1)
            mulαβ!(vec(v[j]), L[Block(i, j)]', vi, -one(T), one(T))
        end
    end
    if !uscale
        for (t, vv) in zip(m.reterms, v)
            lmul!(t.λ, vv)
        end
    end
    v
end

ranef!(v::Vector, m::LinearMixedModel, uscale::Bool) = ranef!(v, m, fixef(m), uscale)

"""
    ranef(m::LinearMixedModel; uscale=false) #, named=true)

Return, as a `Vector{Vector{T}}` (`Vector{NamedVector{T}}` if `named=true`),
the conditional modes of the random effects in model `m`.

If `uscale` is `true` the random effects are on the spherical (i.e. `u`) scale, otherwise on
the original scale.
"""
function ranef(m::LinearMixedModel{T}; uscale=false, named=false) where {T}
    v = [Matrix{T}(undef, size(t.z, 1), nlevs(t)) for t in m.reterms]
    ranef!(v, m, uscale)
    named || return v
    vnmd = map(NamedArray, v)
    for (trm, vnm) in zip(m.reterms, vnmd)
        setnames!(vnm, trm.cnames, 1)
        setnames!(vnm, string.(trm.trm.contrasts.levels), 2)
    end
    vnmd
end

LinearAlgebra.rank(m::LinearMixedModel) = first(m.feterms).rank

function rePCA(m::LinearMixedModel{T}) where {T}
    re = m.reterms
    nt = length(re)
    tup = ntuple(i -> normalized_variance_cumsum(re[i].λ), nt)
    NamedTuple{ntuple(i -> re[i].trm.sym, nt), typeof(tup)}(tup)
end

"""
    reevaluateAend!(m::LinearMixedModel)

Reevaluate the last column of `m.A` from `m.feterms`.  This function should be called
after updating the response, `m.feterms[end]`.
"""
function reevaluateAend!(m::LinearMixedModel)
    A = m.A
    ftrms = m.feterms
    trmn = reweight!(last(ftrms), m.sqrtwts)
    nblk = BlockArrays.nblocks(A, 1)
    for (j, trm) in enumerate(vcat(m.reterms, ftrms))
        mul!(A[Block(nblk, j)], trmn', trm)
    end
    m
end

"""
    refit!(m::LinearMixedModel[, y::Vector])

Refit the model `m` after installing response `y`.

If `y` is omitted the current response vector is used.
"""
refit!(m::LinearMixedModel) = fit!(reevaluateAend!(m))

function refit!(m::LinearMixedModel, y)
    resp = last(m.feterms)
    length(y) == size(resp, 1) || throw(DimensionMismatch(""))
    copyto!(resp, y)
    refit!(m)
end

StatsBase.residuals(m::LinearMixedModel) = response(m) .- fitted(m)

StatsBase.response(m::LinearMixedModel) = vec(last(m.feterms).x)

function reweight!(m::LinearMixedModel, weights)
    sqrtwts = map!(sqrt, m.sqrtwts, weights)
    reweight!.(m.feterms, Ref(sqrtwts))
    reweight!.(m.reterms, Ref(sqrtwts))
    updateA!(m)
    updateL!(m)
end

"""
    sdest(m::LinearMixedModel)

Return the estimate of σ, the standard deviation of the per-observation noise.
"""
sdest(m::LinearMixedModel) = √varest(m)

"""
    setθ!{T}(m::LinearMixedModel{T}, v::Vector{T})

Install `v` as the θ parameters in `m`.
"""
function setθ!(m::LinearMixedModel, v)
    offset = 0
    for trm in m.reterms
        k = nθ(trm)
        setθ!(trm, view(v, (1:k) .+ offset))
        offset += k
    end
    m
end

Base.setproperty!(m::LinearMixedModel, s::Symbol, y) = s == :θ ? setθ!(m, y) : setfield!(m, s, y)

function Base.show(io::IO, m::LinearMixedModel)
    if m.optsum.feval < 0
        @warn("Model has not been fit")
        return nothing
    end
    n, p, q, k = size(m)
    REML = m.optsum.REML
    println(io, "Linear mixed model fit by ", REML ? "REML" : "maximum likelihood")
    println(io, " ", m.formula)
    oo = objective(m)
    if REML
        println(io, " REML criterion at convergence: ", oo)
    else
        nums = showoff([-oo/ 2, oo, aic(m), bic(m)])
        fieldwd = max(maximum(textwidth.(nums)) + 1, 11)
        for label in [" logLik", "-2 logLik", "AIC", "BIC"]
            print(io, rpad(lpad(label, (fieldwd + textwidth(label)) >> 1), fieldwd))
        end
        println(io)
        print.(Ref(io), lpad.(nums, fieldwd))
        println(io)
    end
    println(io)

    show(io,VarCorr(m))

    print(io," Number of obs: $n; levels of grouping factors: ")
    join(io, nlevs.(m.reterms), ", ")
    println(io)
    println(io,"\n  Fixed-effects parameters:")
    show(io,coeftable(m))
end

function σs(m::LinearMixedModel)
    σ = sdest(m)
    NamedTuple{fnames(m)}(((σs(t, σ) for t in m.reterms)...,))
end

function σρs(m::LinearMixedModel)
    σ = sdest(m)
    NamedTuple{fnames(m)}(((σρs(t, σ) for t in m.reterms)...,))
end

function Base.size(m::LinearMixedModel)
    n, p = size(first(m.feterms))
    n, p, sum(size.(m.reterms, 2)), length(m.reterms)
end

"""
    sqrtpwrss(m::LinearMixedModel)

Return the square root of the penalized, weighted residual sum-of-squares (pwrss).

This value is the contents of the `1 × 1` bottom right block of `m.L`
"""
sqrtpwrss(m::LinearMixedModel) = first(m.L.blocks[end, end])

"""
    std(m::MixedModel)

Return the estimated standard deviations of the random effects as a `Vector{Vector{T}}`.
"""
function Statistics.std(m::LinearMixedModel)
    rl = rowlengths.(m.reterms)
    s = sdest(m)
    isfinite(s) ? rmul!(push!(rl, [1.]), s) : rl
end

"""
    updateA!(m::LinearMixedModel)

Update the cross-product array, `m.A`, using `m.reterms` and `m.feterms`

This is usually done after a reweight! operation.
"""
function updateA!(m::LinearMixedModel)
    terms = vcat(m.reterms, m.feterms)
    k = length(terms)
    A = m.A
    for j in 1:k
        for i in j:k
            mul!(A[Block(i, j)], terms[i]', terms[j])
        end
    end
    m
end

"""
    updateL!(m::LinearMixedModel)

Update the blocked lower Cholesky factor, `m.L`, from `m.A` and `m.reterms` (used for λ only)

This is the crucial step in evaluating the objective, given a new parameter value.
"""
function updateL!(m::LinearMixedModel{T}) where {T}
    A = m.A
    L = m.L
    k = BlockArrays.nblocks(A, 2)
    for j in 1:k                         # copy lower triangle of A to L
        for i in j:BlockArrays.nblocks(A, 1)
            copyto!(L[Block(i, j)], A[Block(i, j)])
        end
    end
    for (j, cj) in enumerate(m.reterms)  # pre- and post-multiply by Λ, add I to diagonal
        scaleinflate!(L[Block(j, j)], cj)
        for i in (j+1):k         # postmultiply column by Λ
            rmulΛ!(L[Block(i, j)], cj)
        end
        for jj in 1:(j-1)        # premultiply row by Λ'
            lmulΛ!(cj', L[Block(j, jj)])
        end
    end
    for j in 1:k                         # blocked Cholesky
        Ljj = L[Block(j, j)]
        LjjH = isa(Ljj, Diagonal) ? Ljj : Hermitian(Ljj, :L)
        for jj in 1:(j - 1)
            rankUpdate!(LjjH, L[Block(j, jj)], -one(T))
        end
        cholUnblocked!(Ljj, Val{:L})
        LjjT = isa(Ljj, Diagonal) ? Ljj : LowerTriangular(Ljj)
        for i in (j + 1):k
            Lij = L[Block(i, j)]
            for jj in 1:(j - 1)
                mulαβ!(Lij, L[Block(i, jj)], L[Block(j, jj)]', -one(T), one(T))
            end
            rdiv!(Lij, LjjT')
        end
    end
    m
end

"""
    varest(m::LinearMixedModel)

Returns the estimate of σ², the variance of the conditional distribution of Y given B.
"""
varest(m::LinearMixedModel) = pwrss(m) / dof_residual(m)

function StatsBase.vcov(m::LinearMixedModel{T}) where {T}
    Xtrm = first(m.feterms)
    iperm = invperm(Xtrm.piv)
    p = length(iperm)
    r = Xtrm.rank
    Linv = inv(feL(m))
    permvcov = varest(m) * (Linv'Linv)
    if p == Xtrm.rank
        permvcov[iperm, iperm]
    else
        covmat = fill(zero(T)/zero(T), (p, p))
        for j in 1:r, i in 1:r
            covmat[i,j] = permvcov[i, j]
        end
        covmat[iperm, iperm]
    end
end

"""
    zerocorr!(m::LinearMixedModel[, trmnms::Vector{Symbol}])

Rewrite the random effects specification for the grouping factors in `trmnms` to zero correlation parameter.

The default for `trmnms` is all the names of random-effects terms.

A random effects term is in the zero correlation parameter configuration when the off-diagonal elements of
λ are all zero - hence there are no correlation parameters in that term being estimated.
"""
function zerocorr!(m::LinearMixedModel{T}, trmns) where {T}
    reterms = m.reterms
    for trm in reterms
        if fname(trm) in trmns
            zerocorr!(trm)
        end
    end
    optsum = m.optsum
    optsum.lowerbd = foldl(vcat, lowerbd(c) for c in reterms)
    optsum.initial = foldl(vcat, getθ(c) for c in reterms)
    optsum.final = copy(optsum.initial)
    optsum.xtol_abs = fill!(copy(optsum.initial), 1.0e-10)
    optsum.initial_step = T[]

    # the model is no longer fitted
    optsum.feval == -1

    m
end

zerocorr!(m::LinearMixedModel) = zerocorr!(m, fnames(m))


text:
using Documenter, ColorSchemes

makedocs(
  modules = [ColorSchemes],
  format = :html,
  sitename = "ColorSchemes",
  pages    = Any[
    "Introduction"      => "index.md",
    "Basic usage"       => "basics.md",
    "Finding colors"    => "inverse.md",
    "Plotting"          => "plotting.md",
    "Images"            => "images.md",
    "Index"             => "functionindex.md"
    ]
  )

deploydocs(
    repo = "github.com/JuliaGraphics/ColorSchemes.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
## Look at musical artists
using DelimitedFiles
using SparseArrays
names = readlines("artistsim.names")
data = readdlm("artistsim.smat")
##
data[2:end,1] .+= 1
data[2:end,2] .+= 1
##
A = sparse(Int.(data[2:end,1]),
           Int.(data[2:end,2]),
           Int.(data[2:end,3]),
           Int(data[1,1]),Int(data[1,1]))
## Show the data
using Random
p = randperm(size(data,1))
display(Int.(data[p,:]))
##
names[34261],names[14195]
##
names[11990],names[16868]
## The data is not symmetric.
sum(abs.(A-A'))
## Show the artists with most in-coming similarity
using Plots
plot(vec(sum(A,dims=2)))
##
using Printf
d = vec(sum(A,dims=1))
p = sortperm(d, rev=true)
for i=1:25
    @printf("%2i - %s\n", i, names[p[i]])
end
##
using MatrixNetworks
x = pagerank(A,0.85)
p = sortperm(x,rev=true)
##
for i=1:25
    @printf("%2i - %s\n", i, names[p[i]])
end


text:
#################### Raftery and Lewis Diagnostic ####################
"""
    rafterydiag(x::Vector{T}; q::Real=0.025, r::Real=0.005,
                      s::Real=0.95, eps::Real=0.001,
                      range::AbstractRange=1:1:length(x)) where {T<:Real}
"""
function rafterydiag(x::Vector{T}; q::Real=0.025, r::Real=0.005,
                      s::Real=0.95, eps::Real=0.001,
                      range::AbstractRange=1:1:length(x)) where {T<:Real}
  nx = length(x)
  phi = sqrt(2.0) * erfinv(s)
  nmin = ceil(Int, q * (1.0 - q) * (phi / r)^2)
  if nmin > nx
    @warn("At least $nmin samples are needed for specified q, r, and s")
    kthin = burnin = total = NaN
  else
    dichot = Int[(x .<= quantile(x, q))...]
    kthin = 0
    bic = 1.0
    local test, ntest
    while bic >= 0.0
      kthin += 1
      test = dichot[1:kthin:nx]
      ntest = length(test)
      temp = test[1:(ntest - 2)] + 2 * test[2:(ntest - 1)] + 4 * test[3:ntest]
      trantest = reshape(counts(temp, 0:7), 2, 2, 2)
      g2 = 0.0
      for i1 in 1:2, i2 in 1:2, i3 in 1:2
        tt = trantest[i1, i2, i3]
        if tt > 0
          fitted = sum(trantest[:, i2, i3]) * sum(trantest[i1, i2, :]) /
                   sum(trantest[:, i2, :])
          g2 += 2.0 * tt * log(tt / fitted)
        end
      end
      bic = g2 - 2.0 * log(ntest - 2.0)
    end
    tranfinal = counts(test[1:(ntest - 1)] + 2 * test[2:ntest], 0:3)
    alpha = tranfinal[3] / (tranfinal[1] + tranfinal[3])
    beta = tranfinal[2] / (tranfinal[2] + tranfinal[4])
    kthin *= step(range)
    m = log(eps * (alpha + beta) / max(alpha, beta)) /
        log(abs(1.0 - alpha - beta))
    burnin = kthin * ceil(m) + range.start - 1
    n = ((2.0 - alpha - beta) * alpha * beta * phi^2) /
        (r^2 * (alpha + beta)^3)
    keep = kthin * ceil(n)
    total = burnin + keep
  end
  [kthin, burnin, total, nmin, total / nmin]
end
"""
    rafterydiag(c::AbstractChains; q::Real=0.025, r::Real=0.005,
                     s::Real=0.95, eps::Real=0.001)

Compute the convergence diagnostic of Raftery and Lewis for MCMC sampler output. The diagnostic is designed to determine the number of autocorrelated samples required to estimate a specified quantile ``\\theta_q``, such that ``\\Pr(\\theta \\le \\theta_q) = q``, within a desired degree of accuracy. In particular, if ``\\hat{\\theta}_q`` is the estimand and ``\\Pr(\\theta \\le \\hat{\\theta}_q) = \\hat{P}_q`` the estimated cumulative probability, then accuracy is specified in terms of r and s, where ``\\Pr(q - r < \\hat{P}_q < q + r) = s``. Thinning may be employed in the calculation of the diagnostic to satisfy its underlying assumptions. However, users may not want to apply the same (or any) thinning when estimating posterior summary statistics because doing so results in a loss of information. Accordingly, sample sizes estimated by the diagnostic tend to be conservative (too large).

Returns a `ChainSummary` type object with parameters contained in the rows of the value field, and thinning intervals employed, numbers of samples to discard as burn-in sequences, total numbers ``(N)`` to burn-in and retain, numbers of independent samples that would be needed ``(Nmin)``, and dependence factors ``(N / Nmin)`` in the columns. Results are chain-specific.

* `c` : sampler output on which to perform calculations.

* `q` : posterior quantile of interest.

* `r` : margin of error for estimated cumulative probabilities.

* `s` : probability for the margin of error.

* `eps` : tolerance within which the probabilities of transitioning from initial to retained iterations are within the equilibrium probabilities for the chain. This argument determines the number of samples to discard as a burn-in sequence and is typically left at its default value.

* `x` : vector on which to perform calculations.

* `range` : ??
"""
function rafterydiag(c::AbstractChains; q::Real=0.025, r::Real=0.005,
                     s::Real=0.95, eps::Real=0.001)
  _, p, m = size(c.value)
  vals = Array{Float64}(undef, p, 5, m)
  for j in 1:p, k in 1:m
    vals[j, :, k] = rafterydiag(c.value[:, j, k], q=q, r=r, s=s, eps=eps,
                                range=c.range)
  end
  hdr = header(c) * "\nRaftery and Lewis Diagnostic:\n" *
        "Quantile (q) = $q\nAccuracy (r) = $r\nProbability (s) = $s\n"
  ChainSummary(vals, c.names, ["Thinning", "Burn-in", "Total", "Nmin",
                               "Dependence Factor"], hdr)
end


text:
using Dash, DashHtmlComponents, DashCoreComponents

app = dash()
markdown_text = "
### Dash and Markdown

Dash apps can be written in Markdown.
Dash uses the [CommonMark](http://commonmark.org/)
specification of Markdown.
Check out their [60 Second Markdown Tutorial](http://commonmark.org/help/)
if this is your first introduction to Markdown!
"

app.layout = html_div() do
    dcc_markdown(markdown_text)
end

run_server(app, "0.0.0.0", debug=true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
@testset "device arrays" begin

@testset "constructors" begin
    # inner constructors
    let
        dp = CUDAnative.DevicePtr{Int,AS.Generic}(0)
        CuDeviceArray{Int,1,AS.Generic}((1,), dp)
    end

    # outer constructors
    for I in [Int32,Int64]
        a = I(1)
        b = I(2)

        dp = CUDAnative.DevicePtr{I,AS.Generic}(0)

        # not parameterized
        CuDeviceArray(b, dp)
        CuDeviceArray((b,), dp)
        CuDeviceArray((b,a), dp)

        # partially parameterized
        CuDeviceArray{I}(b, dp)
        CuDeviceArray{I}((b,), dp)
        CuDeviceArray{I}((a,b), dp)
        CuDeviceArray{I,1}(b, dp)
        CuDeviceArray{I,1}((b,), dp)
        @test_throws MethodError CuDeviceArray{I,1}((a,b), dp)
        @test_throws MethodError CuDeviceArray{I,2}(b, dp)
        @test_throws MethodError CuDeviceArray{I,2}((b,), dp)
        CuDeviceArray{I,2}((a,b), dp)

        # fully parameterized
        CuDeviceArray{I,1,AS.Generic}(b, dp)
        CuDeviceArray{I,1,AS.Generic}((b,), dp)
        @test_throws MethodError CuDeviceArray{I,1,AS.Generic}((a,b), dp)
        @test_throws MethodError CuDeviceArray{I,1,AS.Shared}((a,b), dp)
        @test_throws MethodError CuDeviceArray{I,2,AS.Generic}(b, dp)
        @test_throws MethodError CuDeviceArray{I,2,AS.Generic}((b,), dp)
        CuDeviceArray{I,2,AS.Generic}((a,b), dp)

        # type aliases
        CuDeviceVector{I}(b, dp)
        CuDeviceMatrix{I}((a,b), dp)
    end
end

@testset "basics" begin     # argument passing, get and setindex, length
    dims = (16, 16)
    len = prod(dims)

    function kernel(input::CuDeviceArray{Float32}, output::CuDeviceArray{Float32})
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x

        if i <= length(input)
            output[i] = Float64(input[i])   # force conversion upon setindex!
        end

        return
    end

    input = round.(rand(Float32, dims) * 100)

    input_dev = CuTestArray(input)
    output_dev = CuTestArray(input)

    @cuda threads=len kernel(input_dev, output_dev)
    output = Array(output_dev)
    @test input ≈ output
end

@testset "iteration" begin     # argument passing, get and setindex, length
    dims = (16, 16)
    function kernel(input::CuDeviceArray{T}, output::CuDeviceArray{T}) where {T}
        acc = zero(T)
        for elem in input
            acc += elem
        end
        output[1] = acc
        return
    end

    input = round.(rand(Float32, dims) * 100)

    input_dev = CuTestArray(input)
    output_dev = CuTestArray(Float32[0])

    @cuda kernel(input_dev, output_dev)
    output = Array(output_dev)
    @test sum(input) ≈ output[1]
end

@testset "bounds checking" begin
    @testset "#313" begin
        function kernel(dest)
            dest[1] = 1
            nothing
        end
        tt = Tuple{SubArray{Float64,2,CuDeviceArray{Float64,2,AS.Global},
                            Tuple{UnitRange{Int64},UnitRange{Int64}},false}}

        ir = sprint(io->CUDAnative.code_llvm(io, kernel, tt))
        @test !occursin("jl_invoke", ir)
        CUDAnative.code_ptx(devnull, kernel, tt)
    end
end

@testset "views" begin
    function kernel(array)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x

        _sub = view(array, 2:length(array)-1)
        if i <= length(_sub)
            _sub[i] = i
        end

        return
    end

    array = zeros(Int64, 100)
    array_dev = CuTestArray(array)

    sub = view(array, 2:length(array)-1)
    for i in 1:length(sub)
        sub[i] = i
    end

    @cuda threads=100 kernel(array_dev)
    @test array == Array(array_dev)
end

@testset "non-Int index to unsafe_load" begin
    function load_index(a)
        return a[UInt64(1)]
    end

    a = [1]
    p = pointer(a)
    dp = Base.bitcast(CUDAnative.DevicePtr{eltype(p), AS.Generic}, p)
    da = CUDAnative.CuDeviceArray(1, dp)
    load_index(da)
end

@testset "ldg" begin
    function kernel(a, b, i)
        b[i] = ldg(a, i)
        return
    end

    buf = IOBuffer()

    a = CuTestArray([0])
    b = CuTestArray([0])
    @device_code_ptx io=buf @cuda kernel(a, b, 1)
    @test Array(a) == Array(b)

    asm = String(take!(copy(buf)))
    @test occursin("ld.global.nc", asm)
end

end


text:
using Plots
using Test: @test

macro isplot(ex) # @isplot macro to streamline tests
    :(@test $(esc(ex)) isa Plots.Plot)
end


text: _is_incremental(::Type{CPLEX.Optimizer}) = false

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: MatthewStuber/EAGO


text:
using Test
using Yao
using QuAlgorithmZoo

"""
    gaussian_pdf(x, μ::Real, σ::Real)

gaussian probability density function.
"""
function gaussian_pdf(x, μ::Real, σ::Real)
    pl = @. 1 / sqrt(2pi * σ^2) * exp(-(x - μ)^2 / (2 * σ^2))
    pl / sum(pl)
end

@testset "qcbm" begin
    # problem setup
    n = 6
    depth = 6

    N = 1<<n
    kernel = rbf_kernel(0:N-1, 0.25)
    pg = gaussian_pdf(1:N, N/2-0.5, N/4)
    circuit = random_diff_circuit(n, depth, pair_ring(n)) |> autodiff(:QC)
    dispatch!(circuit, :random)
    qcbm = QCBM(circuit, kernel, pg)

    # training
    niter = 100
    optim = Adam(lr=0.1)
    for info in QCBMGo!(qcbm, optim, niter) end
    @test qcbm |> loss < 1e-4
end


text:
using Documenter
using MarketTechnicals

makedocs(
    format = :html,
    sitename = "MarketTechnicals.jl",
    pages = [
        "index.md",
        "getting_started.md",
        "ma.md",
        "levels.md",
        "momentum.md",
        "volatility.md",
        "volume.md",
        "candlesticks.md",
        "utils.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaQuant/MarketTechnicals.jl.git",
    julia  = "0.6",
    latest = "master",
    target = "build",
    deps = nothing,  # we use the `format = :html`, without `mkdocs`
    make = nothing,  # we use the `format = :html`, without `mkdocs`
)


text:
using Test
using CLIMAParameters: AbstractEarthParameterSet
using CLIMAParameters.SubgridScale

@testset "SubgridScale" begin
    struct EarthParameterSet <: AbstractEarthParameterSet end
    ps = EarthParameterSet()
    @test von_karman_const(ps) ≈ 0.4
end

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
##############################################################################
#
# Fallback methods, usually overridden for specific distributions
#
##############################################################################

binaryentropy(d::Distribution) = entropy(d) / log(2)

function cor(d::MultivariateDistribution)
    R = copy(d.cov)
    m, n = size(R)
    for j in 1:n
        for i in 1:n
            R[i, j] = d.cov[i, j] / sqrt(d.cov[i, i] * d.cov[j, j])
        end
    end
    return R
end

cov(d::MultivariateDistribution) = var(d)

ccdf(d::UnivariateDistribution, q::Real) = 1.0 - cdf(d, q)

cquantile(d::UnivariateDistribution, p::Real) = quantile(d, 1.0 - p)

function deviance{M<:Real,Y<:Real,W<:Real}(d::Distribution,
                                           mu::AbstractArray{M},
                                           y::AbstractArray{Y},
                                           wt::AbstractArray{W})
    promote_shape(size(mu), promote_shape(size(y), size(wt)))
    ans = 0.0
    for i in 1:length(y)
        ans += wt[i] * logpdf(d, mu[i], y[i])
    end
    return -2.0 * ans
end

function devresid(d::Distribution, y::Real, mu::Real, wt::Real)
    return -2.0 * wt * logpdf(d, mu, y)
end

function devresid{Y<:Real,M<:Real,W<:Real}(d::Distribution,
                                           y::AbstractArray{Y},
                                           mu::AbstractArray{M},
                                           wt::AbstractArray{W})
    R = Array(Float64, promote_shape(size(y), promote_shape(size(mu), size(wt))))
    for i in 1:length(mu)
        R[i] = devresid(d, y[i], mu[i], wt[i])
    end
    return R
end

function devresid(d::Distribution, y::Vector{Float64},
                  mu::Vector{Float64}, wt::Vector{Float64})
    return [devresid(d, y[i], mu[i], wt[i]) for i in 1:length(y)]
end

function insupport(d::Distribution, x::AbstractArray)
    for e in x
        if !insupport(d, e)
            return false
        end
    end
    return true
end

invlogccdf(d::Distribution, lp::Real) = quantile(d, -expm1(lp))

invlogcdf(d::Distribution, lp::Real) = quantile(d, exp(lp))

# kurtosis returns excess kurtosis by default
# proper kurtosis can be returned with correction = false
function kurtosis(d::Distribution, correction::Bool)
    if correction
        return kurtosis(d)
    else
        return kurtosis(d) + 3.0
    end
end
excess(d::Distribution) = kurtosis(d)
excess_kurtosis(d::Distribution) = kurtosis(d)
proper_kurtosis(d::Distribution) = kurtosis(d, false)

logcdf(d::Distribution, q::Real) = log(cdf(d,q))

logccdf(d::Distribution, q::Real) = log(ccdf(d,q))

logpdf(d::Distribution, x::Real) = log(pdf(d,x))

function logpdf!(r::AbstractArray, d::UnivariateDistribution, x::AbstractArray)
    if size(x) != size(r)
        throw(ArgumentError("Inconsistent array dimensions."))
    end
    for i in 1:length(x)
        r[i] = logpdf(d, x[i])
    end
end

function logpdf(d::MultivariateDistribution, x::AbstractMatrix)
    n::Int = size(x, 2)
    r = Array(Float64, n)
    for i in 1:n
        r[i] = logpdf(d, x[:, i])
    end
    r
end

# function logpdf{T <: Real}(d::Dirichlet, x::Matrix{T})
#     r = Array(Float64, size(x, 2))
#     logpdf!(r, d, x)
#     return r
# end

function logpdf!(r::AbstractArray, d::MultivariateDistribution, x::AbstractMatrix)
    n::Int = size(x, 2)
    if length(r) != n
        throw(ArgumentError("Inconsistent array dimensions."))
    end
    for i = 1:n
        r[i] = logpdf(d, x[:, i])
    end
end

logpmf(d::DiscreteDistribution, args::Any...) = logpdf(d, args...)

function logpmf!(r::AbstractArray, d::DiscreteDistribution, args::Any...)
    return logpdf!(r, d, args...)
end

function mustart{Y<:Real,W<:Real}(d::Distribution,
                                  y::AbstractArray{Y},
                                  wt::AbstractArray{W})
    M = Array(Float64, promote_shape(size(y), size(wt)))
    for i in 1:length(M)
        M[i] = mustart(d, y[i], wt[i])
    end
    return M
end

pmf(d::DiscreteDistribution, args::Any...) = pdf(d, args...)

function pdf(d::MultivariateDistribution, X::AbstractMatrix)
    r = Array(Float64, size(X, 2))
    return pdf!(r, d, X)
end

function pdf!(r::AbstractArray, d::MultivariateDistribution, X::AbstractMatrix)
    for i in 1:size(X, 2)
        r[i] = pdf(d, X[:, i])
    end
    return r
end

function rand!(d::UnivariateDistribution, A::Array)
    for i in 1:length(A)
        A[i] = rand(d)
    end
    return A
end

function rand(d::ContinuousDistribution, dims::Dims)
    return rand!(d, Array(Float64, dims))
end

function rand(d::DiscreteDistribution, dims::Dims)
    return rand!(d, Array(Int, dims))
end

function rand(d::NonMatrixDistribution, dims::Integer...)
    return rand(d, map(int, dims))
end

function rand(d::ContinuousMultivariateDistribution, dims::Integer)
    return rand!(d, Array(Float64, length(mean(d)), dims))
end

function rand(d::DiscreteMultivariateDistribution, dims::Integer)
    return rand!(d, Array(Int, length(mean(d)), dims))
end

function rand(d::MatrixDistribution, dims::Integer)
    return rand!(d, Array(Matrix{Float64}, int(dims)))
end

function rand!(d::MultivariateDistribution, X::Matrix)
    k = length(mean(d))
    m, n = size(X)
    if m != k
        error("Wrong dimensions")
    end
    for i in 1:n
        X[:, i] = rand(d)
    end
    return X
end

function rand!(d::MatrixDistribution, X::Array{Matrix{Float64}})
    for i in 1:length(X)
        X[i] = rand(d)
    end
    return X
end

function sprand(m::Integer, n::Integer, density::Real, d::Distribution)
    return sprand(m, n, density, n -> rand(d, n))
end

std(d::Distribution) = sqrt(var(d))

function var{M <: Real}(d::Distribution, mu::AbstractArray{M})
    V = similar(mu, Float64)
    for i in 1:length(mu)
        V[i] = var(d, mu[i])
    end
    return V
end

# Vectorize methods
for f in (:pdf, :logpdf, :cdf, :logcdf, :ccdf, :logccdf, :quantile,
          :cquantile, :invlogcdf, :invlogccdf)
    @eval begin
        function ($f)(d::UnivariateDistribution, x::AbstractArray)
            res = Array(Float64, size(x))
            for i in 1:length(res)
                res[i] = ($f)(d, x[i])
            end
            return res
        end
    end
end


text:
#!/usr/bin/env julia

#!/usr/bin/env julia
using FASTX
using BioAlignments

targets = readlines(ARGS[1])

# Read sequence dictionary for HLA
seqFile = "data/hla_gen.fa"
type2seq = Dict()
open(FASTA.Reader,seqFile) do reader
	for record in reader
		# HLA:HLA14787_DQA1*01:01:01:03_5667_bp
		# HLA:HLA17305_DQA1*01:01:01:05_6493_bp
		name = FASTA.identifier(record)
		FASTA.identifier

		#type2seq[name]=FASTA.sequence(record)
	end
end

#type2seqID=Dict()
#for (name,seq) in name2seq

#end

fields = split(name,":")
5-element Vector{SubString{String}}:
 "HLA"
 "HLA14787_DQA1*01"
 "01"
 "01"
 "03_5667_bp"

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
function rrsync(source, destination)
    run(`rsync $source $destination`)
end

function mnfsrsync(source::String, destination, user="data")
    rrsync(mnfs2rsync(source, destination, user)...)
end

function mnfsrsync(source::Array{T,1}, destination, user="data") where T
	files = cat(2,mnfs2rsync.(source, [destination])...)
	tname = tempname()
	writedlm(tname,files[1,:])
	run(`rsync --progress --no-R --files-from $tname / $(files[2,1])`)
end

function mnfs2rsync(source::String, destination, user="data")
    ssource = split(source, '/')
    sdest = split(destination, '/')
    if (ssource[2] == "mnfs") & (sdest[2] == "mnfs")
        error("Both source & destination paths cannot be remote")
    elseif ssource[2] == "mnfs"
        mnfspath = ssource; rev=true
    elseif sdest[2] == "mnfs"
        mnfspath = sdest; rev=false
    else
        warn("Both local paths")
        return [source, destination]
    end
    if isfile("/etc/auto.nfs")
        autonfs = readdlm("/etc/auto.nfs", use_mmap=false)
        server = split(autonfs[autonfs[:,1].==mnfspath[3],3][1],':')[1]
        server = split(server, '.')[1]
    else
        error("Need /etc/auto.nfs to resolve server location")
    end
    localpath = "/$(join(mnfspath[3:end], '/'))"
    rsyncpath = "$user@$server:/$localpath"
    if rev
        return [rsyncpath, destination]
    else
        return [source, rsyncpath]
    end
end


text:
type FakeDriveAction <: DriveAction end
type FakeActionContext <: ActionContext end

let
    trajdata = get_test_trajdata()
    roadway = trajdata.roadway
    veh = get_vehicle(trajdata, 1, 1)

    @test_throws ErrorException length(FakeDriveAction)
    @test_throws ErrorException convert(FakeDriveAction, Float64[])
    @test_throws ErrorException convert(Vector{Float64}, FakeDriveAction())
    @test_throws ErrorException copy!(Float64[], FakeDriveAction())
    @test_throws ErrorException propagate(veh, FakeDriveAction(), FakeActionContext(), roadway)

    let
        a = AccelTurnrate(0.1,0.2)
        @test a == convert(AccelTurnrate, [0.1,0.2])
        @test convert(Vector{Float64}, a) == [0.1,0.2]
        @test copy!([NaN, NaN], AccelTurnrate(0.1,0.2)) == [0.1,0.2]

        s = propagate(veh, AccelTurnrate(0.0,0.0), IntegratedContinuous(1.0,1), roadway)
        @test isapprox(s.posG.x, veh.state.v)
        @test isapprox(s.posG.y, 0.0)
        @test isapprox(s.posG.θ, 0.0)
    end

    let
        a = AccelDesang(0.1,0.2)
        @test a == convert(AccelDesang, [0.1,0.2])
        @test convert(Vector{Float64}, a) == [0.1,0.2]
        @test copy!([NaN, NaN], AccelDesang(0.1,0.2)) == [0.1,0.2]

        s = propagate(veh, AccelDesang(0.0,0.0), IntegratedContinuous(1.0,1), roadway)
        @test isapprox(s.posG.x, veh.state.v)
        @test isapprox(s.posG.y, 0.0)
        @test isapprox(s.posG.θ, 0.0)
    end

    let
        a = NextState(VehicleState())
        @test VehicleState() == propagate(veh, a, ContextFree(), roadway)
    end
end