text:
const flowerFolder = raw"C:\Users\James\.julia\dev\MushroomBody\ImageData\kaggle\natural_images\flower"
const imageFolder = raw"C:\Users\James\.julia\dev\MushroomBody\ImageData\kaggle\natural_images"

function ColorInput(arraysize::Tuple; folder=flowerFolder, normalise_to=350, stages=[0,1,0], input_bool=Bool[0,1,0], reward_bool=Bool[0,1,0], punishment_bool=Bool[0,0,0])
	img = random_image(folder)
	init = normalise(bin(get_hues(img),arraysize[1]),normalise_to)
	return ColorInput(init, stages, input_bool, reward_bool, punishment_bool, cumsum(stages),img)
end

function ColorInput(img::Image, arraysize::Tuple; normalise_to=350, stages=[0,1,0], input_bool=Bool[0,1,0], reward_bool=Bool[0,1,0], punishment_bool=Bool[0,0,0])

	ColorInput(normalise(bin(get_hues(img), arraysize[1]), normalise_to), stages, input_bool, reward_bool, punishment_bool, cumsum(stages), img)
end

function normalise(array, norm_to)

	# normalises the hues by the maximum
	array .* (norm_to / maximum(array))
end

function display_hues(hues)

	minH = 0
	maxH = 360
	hue_list = range(minH, maxH, length=length(hues))
	out = Array{HSL,1}(undef,length(hues))
	hues = normalise(hues,1)

	for i = 1:length(hues)

		out[i] = HSL(hue_list[i],hues[i],0.5)

	end
	return out
end

function many_folders_sequence(arraysize, nstim; folder=imageFolder, stages=[0,1,0], input_bool=Bool[0,1,0], reward_bool=Bool[0,1,0], punishment_bool=Bool[0,0,0])

	out = Vector{Inputs}(undef, nstim)
	multifolders = readdir(folder)
	for i = 1:nstim
		image_folder = folder * "\\" * rand(multifolders)
		out[i] = ColorInput(arraysize, folder=image_folder, stages=stages, input_bool=input_bool, reward_bool=reward_bool, punishment_bool=punishment_bool)
	end

	return InputSequence(out, RestInput(arraysize))
end

function color_sequence(arraysize, nstim; folder=flowerFolder, stages=[0,1,0], input_bool=Bool[0,1,0], reward_bool=Bool[0,1,0], punishment_bool=Bool[0,0,0])

	out = Vector{Inputs}(undef, nstim)
	for i = 1:nstim
		out[i] = ColorInput(arraysize, folder=folder, stages=stages, input_bool=input_bool, reward_bool=reward_bool, punishment_bool=punishment_bool)
	end

	return InputSequence(out, RestInput(arraysize))
end


text:
module GeoTime

using Dates, Printf

import Dates: UTInstant, Millisecond
import Dates: year,  month,  day, hour, minute, second, millisecond
import Dates: daysinmonth, daysinyear, yearmonthday, yearmonth
import Dates: monthday, len

import Base: +, -, isless, string, show, convert, reinterpret

export daysinmonth, daysinyear, yearmonthday, yearmonth, monthday

export timeencode, timedecode, datetuple

export DateTimeStandard, DateTimeJulian, DateTimeProlepticGregorian,
       DateTimeAllLeap, DateTimeNoLeap, DateTime360Day, AbstractCFDateTime

include("time.jl")

end

ext: jl


text:
# Differential equations with mixed units using ArrayPartition and no specified algorithm
using MechGlueDiffEqBase
using MechanicalUnits: @import_expand, ∙, g
@import_expand(km, N, s, m, km, kg, °, inch)
using OrdinaryDiffEq: ODEProblem
using DifferentialEquations: solve, Tsit5
using Logging: disable_logging, LogLevel, Debug, with_logger
import Logging

using Test: @inferred

# Constants we don't think we'll change ever
const x₀ = 0.0km
const y₀ = 0.0km
const ø = 15inch
const ρ = 1.225kg/m³

# Calculated constants
const Aₚᵣ = π/4 * ø^2
# Constants that we define as functions, because we may
# want to modify them later in the same scripting session.
α₀()  = 30°
v₀()  = 1050m/s

mₚ()   = 495kg
C_s() = 0.4
x´₀() = v₀() * cos(α₀())
y´₀() = v₀() * sin(α₀())

# Local tuple initial condition
u₀ = ArrayPartition([x₀, y₀], [x´₀(), y´₀()])

# Functions
v(vx, vy) = √(vx^2 + vy^2)
R(vx, vy) = 0.5∙ρ∙C_s()∙Aₚᵣ∙v(vx, vy)^2
α(vx, vy) = atan(vy, vx)
Rx(vx, vy) = R(vx, vy) * cos(α(vx, vy))
Ry(vx, vy) = R(vx, vy) * sin(α(vx, vy))

# Local tuple, i.e. the interesting degrees of freedom
# and their derivatives
function  Γ!(u´, u, p, t)
    (x, y), (x´, y´) = packout(u)
    @debug "Γ"    x    y    x´   y´   maxlog = 3
    # Calculate the acceleration for this step
    x´´ =     -Rx(x´, y´) / mₚ()
    y´´ = -1g -Ry(x´, y´) / mₚ()
    @debug "Γ Acceleration"    x´´|> g    y´´|>g    maxlog = 3
    packin!(u´, x´, y´, x´´, y´´)
end

packout(u) = u.x[1], u.x[2]
function packin!(u´, x´, y´, x´´, y´´)
    # We access the values through its field name, x (unrelated to the specific problem)
    u´.x[1][1] = x´;  u´.x[1][2] = y´
    u´.x[2][1] = x´´; u´.x[2][2] = y´´
    return u´
end

function solve_guarded(u₀; alg = Tsit5(), debug=false)
    # Test the functions
    @inferred packout(u₀)
    @inferred packout(0.01u₀/s)
    @inferred Γ!(u₀/s,u₀, nothing, nothing)
    !debug  && disable_logging(Debug)
    prob = ODEProblem( Γ!,u₀,(0.0,60)s)
    sol = if debug
        with_logger(Logging.ConsoleLogger(stderr, Logging.Debug)) do
            solve(prob, Tsit5())
        end
    else
        solve(prob, Tsit5())
    end
    # Re-enable
    disable_logging(LogLevel(Debug-1))
    sol
end
@time sol = solve_guarded(u₀, debug=true);
nothing


text:
#----------------------------------------------------------------
#
# Representation of transcripts similar to GTF
#

export Transcript,
    Gene,
    iscoding,
    exons,
    cdss,
    GTFFile,
    transcripts,
    genes,
    tabix

#abstract type AbstractGenomicFeature end
abstract type AbstractGene <: AbstractGenomicFeature end

struct Transcript <: AbstractGenomicFeature
    chrom::AbstractString
    bpStart::Int64
    bpEnd::Int64
    strand::Char
    id::AbstractString
    name::AbstractString
    exons::Vector{Tuple{Int64, Int64}}
    cdss::Vector{Tuple{Int64, Int64}}
end

name(t::Transcript) = t.name
iscoding(t::Transcript) = length(t.cdss) > 0
exons(t::Transcript) = t.exons
cdss(t::Transcript) = t.cdss

struct Gene <: AbstractGene
    chrom::AbstractString
    bpStart::Int64
    bpEnd::Int64
    strand::Char
    id::AbstractString
    name::AbstractString
    transcripts::Vector{Transcript}
end

name(g::Gene) = g.name
iscoding(g::Gene) = any(iscoding.(g.transcripts))
transcripts(g::Gene) = g.transcripts
exons(g::Gene) = union([exons(t) for t in transcripts(g)]...)
cdss(g::Gene) = union([cdss(t) for t in transcripts(g)]...)

Base.show(io::IO, g::Gene) = print(io, "Gene $(g.name) ($(g.chrom):$(g.bpStart)-$(g.bpEnd)) with $(length(g.transcripts)) transcripts")

struct GTFFile
    path::AbstractString
    io::IO
end

import GZip
GTFFile(filename::AbstractString) = GTFFile(filename, endswith(filename, ".gz") ? GZip.open(filename) : open(filename, "r"))
GTFFile(io::IO) = GTFFile("", io)

function genes(gtf::GTFFile)

    # genes, transcripts and exons might come in any order
    genes = Gene[]
    gene2transcripts = Dict{AbstractString, Vector{Transcript}}()
    transcript2exons = Dict{AbstractString, Vector{Tuple{Int64, Int64}}}()
    transcript2cdss = Dict{AbstractString, Vector{Tuple{Int64, Int64}}}()
    count = 1
    for line in eachline(gtf.io)

        count += 1
        if startswith(line, "#")
            continue
        end

        chrom, source, feature, bpStart, bpEnd, score, strand, frame, attributes = split(chomp(line), "\t")
        attributes = parse_attributes(attributes)

        if feature == "gene"
            push!(genes, Gene(chrom, parse(Int64, bpStart), parse(Int64, bpEnd), strand[1],
                              attributes["gene_id"], attributes["gene_name"], Transcript[]))
        elseif feature == "transcript"
            t = Transcript(chrom, parse(Int64, bpStart), parse(Int64, bpEnd), strand[1],
                           attributes["transcript_id"], attributes["transcript_name"],
                           Tuple{Int64, Int64}[], Tuple{Int64, Int64}[])
            geneid = attributes["gene_id"]
            ~haskey(gene2transcripts, geneid) && (gene2transcripts[geneid] = Transcript[])
            push!(gene2transcripts[geneid], t)
        elseif feature == "exon"
            e = (parse(Int64, bpStart), parse(Int64, bpEnd))
            transcriptid = attributes["transcript_id"]
            ~haskey(transcript2exons, transcriptid) && (transcript2exons[transcriptid] = Tuple{Int64, Int64}[])
            push!(transcript2exons[transcriptid], e)
        elseif feature == "CDS" || feature == "cds"
            c = (parse(Int64, bpStart), parse(Int64, bpEnd))
            transcriptid = attributes["transcript_id"]
            ~haskey(transcript2cdss, transcriptid) && (transcript2cdss[transcriptid] = Tuple{Int64, Int64}[])
            push!(transcript2cdss[transcriptid], c)
        end
    end
    close(gtf.io)
    #info("Parsed $count lines from $(gtf.path)")

    for g in genes
        for t in gene2transcripts[g.id]
            push!(g.transcripts, t)
            for e in transcript2exons[t.id]
                push!(t.exons, e)
            end
            for c in get(transcript2cdss, t.id, Tuple{Int64, Int64}[])
                push!(t.cdss, c)
            end
        end
    end
    genes
end

function transcripts(gtf::GTFFile)
    gs = genes(gtf)
    transcripts = Transcript[]
    for g in gs
        append!(transcripts, g.transcripts)
    end
    transcripts
end

function parse_attributes(attributes)
    res = Dict()
    for kv in split(attributes, ";", keepempty=false)
        k, v = split(kv, " ", keepempty=false)
        res[k] = replace(v, "\"" => "")
    end
    res
end

# get file handle to tabix stream
function tabix(filename::AbstractString, chrom::AbstractString, bpStart::Int64, bpEnd::Int64)
    cmd = `tabix $filename $chrom:$bpStart-$bpEnd`
    process = open(cmd, "r")
    process
end


text:
using Plots;
using Printf;
gr();
anim = Animation();
f(x) = 1/2*cos(8pi*(x-1/2));

function forward(u_old::Array, i::Int64)
	return (3u_old[ir[i]] + 3u_old[i] - 7u_old[il[i]] + u_old[il2[i]]) / 8;
end

function backward(u_old::Array, i::Int64)
	return (3u_old[il[i]] + 3u_old[i] - 7u_old[ir[i]] + u_old[ir2[i]]) / 8;
end

function main()
#=********Numerical Set Up*********=#
	global c = 1.0;
	global nx = 101;
	global lx = 100;
	global x = range(0.0, stop=lx, length=nx);
	global dx = lx /(nx-1);
	global u  = zeros(Float64, nx);
	@inbounds for i = 1:nx
		if 40.0 < x[i] < 60.0;
			u[i] = f(x[i]);
		end
	end
	global t = 0.0;
	global dt = 2e-2;
	global tlims = 300;
	global C = dt * c / dx;
	global ir = zeros(Int64, nx);
	global ir2 = zeros(Int64, nx);
	global il = zeros(Int64, nx);
	global il2 = zeros(Int64, nx);
	@inbounds for i = 1:nx
		ir[i] = i + 1;
		ir2[i] = i + 2;
		il[i] = i - 1;
		il2[i] = i - 2;
	end
	ir[end] = 1;
	ir2[end - 1] = 1;
	ir2[end] = 2;
	il[begin] = nx;
	il2[begin + 1] = nx;
	il2[begin] = nx - 1;
#=******Finish Numerical Set Up*************=#

	while t < tlims
		u_old = copy(u);
		##First substep
		@inbounds for i = 1:nx-1
			u[i] = u_old[i] - C * forward(u_old, i);
		end
		u[end] = u[begin];

		##Second substep
		@inbounds for i = 1:nx
			u_old[i] = 0.5*(u_old[i] + u[i]) - 0.5*C/8.0*(forward(u, i));
		end
		u_old[end] = u_old[begin];

		u = copy(u_old);

		plt = plot(x, u, xlims=(0, lx), xticks=0:lx/5:lx, ylims=(-0.2, 1.2), yticks=-0.2:0.2:1.2,
		xlabel="x", ylabel="u", color="blue", title="QUICK method by using Heun", legend=false);
		frame(anim, plt);
		t += dt;
		@printf("t = %.4f\n", t);
	end
	gif(anim, "quick.gif", fps=5);
end

main();


text:
#-------------------------------------------------------------------------------
# Decomposition
#-------------------------------------------------------------------------------
export proj_to_irrep, block_decomposition
"""
Transform matrix from a (reducible) representation to a given irrep.

For a representation D(g) and irrep d(g), find S that
    S⁺⋅D(g)⋅S = d(g)
"""
function proj_to_irrep(
    rep::AbstractVector{<:AbstractMatrix},
    irep::AbstractVector{<:AbstractMatrix};
    R::Bool=false
)
    N, D = length(rep), size(irep[1], 1)
    vecs = map(1:D) do i
        Pii = proj_operator(rep, [conj(m[i,i]) for m in irep])
        e, v = eigen(Hermitian(Pii))
        v[:, e .> 1e-7]
    end
    for i = 2:D
        Pi1 = proj_operator(rep, [conj(m[i,1]) for m in irep])
        U = vecs[i]' * Pi1 * vecs[1] * D / N
        #@assert norm(U*U'-I) < 1e-7 "Not unitary."
        vecs[i] = vecs[i] * U
    end
    spaces = map(1:size(vecs[1], 2)) do i
        hcat((vecs[j][:, i] for j=1:D)...)
    end
    if R && sum(norm.(imag.(rep))) > 1e-7
        [sqrt(2)*[real(s) imag(s)] for s in spaces]
    else
        spaces
    end
end
#-------------------------------------------------------------------------------
"""
Find a list of projection Sᵢ that
    Sᵢ⁺⋅D(g)⋅Sᵢ = dᵢ(g),
where {dᵢ} is a list of irreducible representations as an input.
"""
function proj_to_irrep(
    rep::AbstractVector{<:AbstractMatrix},
    ireps::AbstractVector{<:AbstractVector{<:AbstractMatrix}};
    R::Bool=false
)
    vcat([proj_to_irrep(rep, irep, R=R) for irep in ireps]...)
end
#-------------------------------------------------------------------------------
function block_decomposition(
    rep::AbstractVector{<:AbstractMatrix},
    group::AbstractFiniteGroup;
    R::Bool=true
)
    ct = charactertable(group)
    row = single_complex_row(ct)
    ireps =  irreps(group, ct[row, :])
    proj_to_irrep(rep, ireps, R=R)
end
#-------------------------------------------------------------------------------
function block_decomposition(
    rep::AbstractVector{<:AbstractMatrix},
    group::String;
    R::Bool=true
)
    block_decomposition(rep, pointgroup(group), R=R)
end
#-------------------------------------------------------------------------------
function block_decomposition(rep::AbstractVector{<:AbstractMatrix}; R::Bool=true)
    n = length(rep)
    multab = Matrix{Integer}(undef, n, n)
    for i=1:n, j=1:n
        gij = rep[i]*rep[j]
        gij_index = 0
        for k = 1:n
            if norm(rep[k] - gij) < 1e-10
                @assert iszero(gij_index) "Not faithful representation. An explicit `FiniteGroup` should be provided."
                gij_index = k
            end
        end
        @assert gij_index > 0 "The input `rep` is not a representation."
        multab[i,j] = gij_index
    end
    group = FiniteGroup(multab)
    block_decomposition(rep, group, R=R)
end


text:
module Test_cplyap

using LinearAlgebra
using MatrixEquations
using Test

@testset "Testing positive continuous Lyapunov equation solvers" begin

n = 30
m = 5
p = 3
Ty = Float64

@testset "Positive continuous Lyapunov equations" begin

reltol = eps(float(100))
a = -1; b = 2im; @time u = plyapc(a,b)
@test abs(a*u*u'+u*u'*a'+b*b') < reltol

reltol = eps(float(100f0))
a = -1f0; b = 2f0im; @time u = plyapc(a,b)
@test abs(a*u*u'+u*u'*a'+b*b') < reltol

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ar = ar-2*norm(ar)*Matrix(I,n,n)
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
ac = ac-2*norm(ac)*Matrix(I,n,n)
br = rand(Ty,n,m)
brw = rand(Ty,n,n+m)
bc = br+im*rand(Ty,n,m)
bcw = brw+im*rand(Ty,n,n+m)
cr = rand(Ty,p,n)
crt = rand(Ty,n+p,n)
cc = cr+im*rand(Ty,p,n)
cct = crt+im*rand(Ty,n+p,n)
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

@time u = plyapc(ar,br);
x = u*u'; @test norm(ar*x+x*ar'+br*br')/norm(x)/norm(ar) < reltol

@time u = plyapc(ar,0*br);
x = u*u'; @test norm(ar*x+x*ar') < reltol

@time u = plyapc(ar,brw);
x = u*u'; @test norm(ar*x+x*ar'+brw*brw')/norm(x)/norm(ar) < reltol

@time u = plyapc(ar,I,br);
x = u*u'; @test norm(ar*x+x*ar'+br*br')/norm(x)/norm(ar) < reltol

@time u = plyapc(ar',cr');
x = u'*u; @test norm(ar'*x+x*ar+cr'*cr)/norm(x)/norm(ar) < reltol

@time u = plyapc(ar',crt');
x = u'*u; @test norm(ar'*x+x*ar+crt'*crt)/norm(x)/norm(ar) < reltol

@time u = plyapc(ac,bc);
x = u*u'; @test norm(ac*x+x*ac'+bc*bc')/norm(x)/norm(ac) < reltol

@time u = plyapc(ac,0*bc);
x = u*u'; @test norm(ac*x+x*ac') < reltol

@time u = plyapc(ac,bcw);
x = u*u'; @test norm(ac*x+x*ac'+bcw*bcw')/norm(x)/norm(ac) < reltol

@time u = plyapc(ac',cc');
x = u'*u; @test norm(ac'*x+x*ac+cc'*cc)/norm(x)/norm(ac) < reltol

@time u = plyapc(ac',cct');
x = u'*u; @test norm(ac'*x+x*ac+cct'*cct)/norm(x)/norm(ac) < reltol

@time u = plyapc(ac',I,cc');
x = u'*u; @test norm(ac'*x+x*ac+cc'*cc)/norm(x)/norm(ac) < reltol

@time u = plyapc(ac,br);
x = u*u'; @test norm(ac*x+x*ac'+br*br')/norm(x)/norm(ac) < reltol

@time u = plyapc(ac',cr');
x = u'*u; @test norm(ac'*x+x*ac+cr'*cr)/norm(x)/norm(ac) < reltol

@time u = plyapc(ar,bc);
x = u*u'; @test norm(ar*x+x*ar'+bc*bc')/norm(x)/norm(ar) < reltol

@time u = plyapc(ar',cc');
x = u'*u; @test norm(ar'*x+x*ar+cc'*cc)/norm(x)/norm(ar)  < reltol

end
end

@testset "Positive continuous generalized Lyapunov equations" begin

reltol = eps(float(100))
a = -1; ee = 4; b = 2im; @time u = plyapc(a,ee,b)
@test abs(a*u*u'*ee'+ee*u*u'*a'+b*b') < reltol

reltol = eps(float(100f0))
a = -1f0; ee = 4f0; b = 2f0im; @time u = plyapc(a,ee,b)
@test abs(a*u*u'*ee'+ee*u*u'*a'+b*b')  < reltol

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ar = ar-2*norm(ar)*Matrix(I,n,n)
er = rand(n,n)
ar = er*ar
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
ac = ac-2*norm(ac)*Matrix(I,n,n)
ec = er+im*rand(Ty,n,n)
ac = ec*ac
br = rand(Ty,n,m)
brw = rand(Ty,n,n+m)
bc = br+im*rand(Ty,n,m)
bcw = brw+im*rand(Ty,n,n+m)
cr = rand(Ty,p,n)
crt = rand(Ty,n+p,n)
cc = cr+im*rand(Ty,p,n)
cct = crt+im*rand(Ty,n+p,n)
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

@time u = plyapc(ar,er,br);
x = u*u'; @test norm(ar*x*er'+er*x*ar'+br*br')/norm(x)/norm(ar) < reltol

@time u = plyapc(ar,er,0*br);
x = u*u'; @test norm(ar*x*er'+er*x*ar') < reltol

@time u = plyapc(ar,er,brw);
x = u*u'; @test norm(ar*x*er'+er*x*ar'+brw*brw')/norm(x)/norm(ar) < reltol

@time u = plyapc(ar',er',cr');
x = u'*u; @test norm(ar'*x*er+er'*x*ar+cr'*cr)/norm(x)/norm(ar) < reltol

@time u = plyapc(ar',er',crt');
x = u'*u; @test norm(ar'*x*er+er'*x*ar+crt'*crt)/norm(x)/norm(ar) < reltol

@time u = plyapc(ac,ec,bc);
x = u*u'; @test norm(ac*x*ec'+ec*x*ac'+bc*bc')/norm(x)/norm(ac)/norm(ec) < reltol

@time u = plyapc(ac,ec,0*bc);
x = u*u'; @test norm(ac*x*ec'+ec*x*ac') < reltol

@time u = plyapc(ac,ec,bcw);
x = u*u'; @test norm(ac*x*ec'+ec*x*ac'+bcw*bcw')/norm(x)/norm(ac)/norm(ec) < reltol

@time u = plyapc(ac',ec',cc');
x = u'*u; @test norm(ac'*x*ec+ec'*x*ac+cc'*cc)/norm(x)/norm(ac) < reltol

@time u = plyapc(ac',ec',cct');
x = u'*u; @test norm(ac'*x*ec+ec'*x*ac+cct'*cct)/norm(x)/norm(ac) < reltol

@time u = plyapc(ar,er,bc);
x = u*u'; @test norm(ar*x*er'+er*x*ar'+bc*bc')/norm(x)/norm(ar) < reltol

@time u = plyapc(ar',er',cc');
x = u'*u; @test norm(ar'*x*er+er'*x*ar+cc'*cc)/norm(x)/norm(ar) < reltol

@time u = plyapc(ac,ec,br);
x = u*u'; @test norm(ac*x*ec'+ec*x*ac'+br*br')/norm(x)/norm(ac)/norm(ec) < reltol

@time u = plyapc(ac',ec',cr');
x = u'*u; @test norm(ac'*x*ec+ec'*x*ac+cr'*cr)/norm(x)/norm(ac) < reltol

end
end

@testset "Positive 2x2 continuous Lyapunov equations" begin

A = [-1.1 1.; -1. -1.]
E = [1. 1.; 0. 1.]
R = UpperTriangular([1. 1.; 0. 1.])
reltol = eps(float(100))

U = copy(R)
β, α = MatrixEquations.plyap2!(A, U, adj = true, disc = false)
X = U'*U; @test norm(A'*X+X*A+R'*R)/max(1,norm(X))/norm(A) < reltol &&
                norm(β*U-U*A)/max(1,norm(U))/norm(A) < reltol &&
                norm(α*U - R)/max(1,norm(R)) < reltol

U = copy(R)
β, α = MatrixEquations.plyap2!(A, U, adj = false, disc = false)
X = U*U'; @test norm(A*X+X*A'+R*R')/max(1,norm(X))/norm(A) < reltol &&
                norm(U*β-A*U)/max(1,norm(U))/norm(A) < reltol &&
                norm(U*α - R)/max(1,norm(R)) < reltol

U = copy(R)
β, α = MatrixEquations.pglyap2!(A, E, U, adj = true, disc = false)
X = U'*U; @test norm(A'*X*E+E'*X*A+R'*R)/max(1,norm(X))/norm(A) < reltol &&
                norm(β*U*E-U*A)/max(1,norm(U))/norm(A) < reltol &&
                norm(α*U*E - R)/max(1,norm(R)) < reltol

U = copy(R)
β, α = MatrixEquations.pglyap2!(A, E, U, adj = false, disc = false)
X = U*U'; @test norm(A*X*E'+E*X*A'+R*R')/max(1,norm(X))/norm(A) < reltol &&
                norm(A*U-E*U*β)/max(1,norm(U))/norm(A) < reltol &&
                norm(E*U*α - R)/max(1,norm(R)) < reltol

U = copy(R)
Q = qr(rand(2,2)).Q; A1 = Q*A; E1 = Q*E;
β, α = MatrixEquations.pglyap2!(A1, E1, U, adj = false, disc = false)
X = U*U'; @test norm(A1*X*E1'+E1*X*A1'+R*R')/max(1,norm(X))/norm(A) < reltol &&
                norm(A1*U-E1*U*β)/max(1,norm(U))/norm(A) < reltol &&
                norm(E1*U*α - R)/max(1,norm(R)) < reltol

U = copy(0*R)
β, α = MatrixEquations.pglyap2!(A, E, U, adj = false, disc = false)
X = U*U'; @test norm(A*X*E'+E*X*A')/max(1,norm(X))/norm(A) < reltol &&
                norm(A*U-E*U*β)/max(1,norm(U))/norm(A) < reltol &&
                norm(E*U*α - 0*R)/max(1,norm(R)) < reltol

A = [-1.1 1.; -1. -1.]/10; A = convert(Matrix{Float32},A)
E = [1. 1.; 0. 1.]; E = convert(Matrix{Float32},E)
R = [1. 1.; 0. 1.]; R = UpperTriangular(convert(Matrix{Float32},R))
reltol = eps(100f0)

U = copy(R)
β, α = MatrixEquations.plyap2!(A, U, adj = true, disc = false)
X = U'*U; @test norm(A'*X+X*A+R'*R)/max(1,norm(X))/norm(A) < reltol &&
                norm(β*U-U*A)/max(1,norm(U))/norm(A) < reltol &&
                norm(α*U - R)/max(1,norm(R)) < reltol

U = copy(R)
β, α = MatrixEquations.plyap2!(A, U, adj = false, disc = false)
X = U*U'; @test norm(A*X+X*A'+R*R')/max(1,norm(X))/norm(A) < reltol &&
                norm(U*β-A*U)/max(1,norm(U))/norm(A) < reltol &&
                norm(U*α - R)/max(1,norm(R)) < reltol

U = copy(R)
β, α = MatrixEquations.pglyap2!(A, E, U, adj = true, disc = false)
X = U'*U; @test norm(A'*X*E+E'*X*A+R'*R)/max(1,norm(X))/norm(A) < reltol &&
                norm(β*U*E-U*A)/max(1,norm(U))/norm(A) < reltol &&
                norm(α*U*E - R)/max(1,norm(R)) < reltol

U = copy(R)
β, α = MatrixEquations.pglyap2!(A, E, U, adj = false, disc = false)
X = U*U'; @test norm(A*X*E'+E*X*A'+R*R')/max(1,norm(X))/norm(A) < reltol &&
                norm(A*U-E*U*β)/max(1,norm(U))/norm(A) < reltol &&
                norm(E*U*α - R)/max(1,norm(R)) < reltol

U = copy(0*R)
β, α = MatrixEquations.pglyap2!(A, E, U, adj = false, disc = false)
X = U*U'; @test norm(A*X*E'+E*X*A')/max(1,norm(X))/norm(A) < reltol &&
                norm(A*U-E*U*β)/max(1,norm(U))/norm(A) < reltol &&
                norm(E*U*α)/max(1,norm(R)) < reltol

end

@testset "Continuous positive Lyapunov equations - Schur form" begin

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ar = ar-2*norm(ar)*Matrix(I,n,n)
as,  = schur(ar)
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
ac = ac-2*norm(ac)*Matrix(I,n,n)
acs,  = schur(ac)
br = rand(Ty,n,m)
brw = rand(Ty,n,n+m)
bc = br+im*rand(Ty,n,m)
cr = rand(Ty,p,n)
crt = rand(Ty,n+p,n)
cc = cr+im*rand(Ty,p,n)
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

@time u = plyaps(as,br);
x = u*u'; @test norm(as*x+x*as'+br*br')/norm(x)/norm(as) < reltol

@time u = plyaps(as,brw);
x = u*u'; @test norm(as*x+x*as'+brw*brw')/norm(x)/norm(as) < reltol

@time u = plyaps(as,I,br);
x = u*u'; @test norm(as*x+x*as'+br*br')/norm(x)/norm(as) < reltol

@time u = plyaps(as',cr');
x = u'*u; @test norm(as'*x+x*as+cr'*cr)/norm(x)/norm(as) < reltol

@time u = plyaps(as',crt');
x = u'*u; @test norm(as'*x+x*as+crt'*crt)/norm(x)/norm(as) < reltol

@time u = plyaps(acs,bc);
x = u*u'; @test norm(acs*x+x*acs'+bc*bc')/norm(x)/norm(as) < reltol

@time u = plyaps(acs',cc');
x = u'*u; @test norm(acs'*x+x*acs+cc'*cc)/norm(x)/norm(as) < reltol

@time u = plyaps(acs',I,cc');
x = u'*u; @test norm(acs'*x+x*acs+cc'*cc)/norm(x)/norm(as) < reltol

F = UpperTriangular(rand(Ty,n,n))
R = copy(F)
@time plyapcs!(as,R,adj = true);
x = R'*R; @test norm(as'*x+x*as+F'*F)/norm(x)/norm(as) < reltol

F = UpperTriangular(rand(Ty,n,n))
R = copy(F)
@time plyapcs!(as,I,R,adj = true);
x = R'*R; @test norm(as'*x+x*as+F'*F)/norm(x)/norm(as) < reltol

F = UpperTriangular(rand(Ty,n,n))
R = copy(F)
@time plyapcs!(as,R,adj = false);
x = R*R'; @test norm(as*x+x*as'+F*F')/norm(x)/norm(as) < reltol

F = UpperTriangular(rand(Ty,n,n)+im*rand(Ty,n,n))
R = copy(F)
@time plyapcs!(acs,R,adj = true);
x = R'*R; @test norm(acs'*x+x*acs+F'*F)/norm(x)/norm(acs) < reltol

F = UpperTriangular(rand(Ty,n,n)+im*rand(Ty,n,n))
R = copy(F)
@time plyapcs!(acs,R,adj = false);
x = R*R'; @test norm(acs*x+x*acs'+F*F')/norm(x)/norm(acs) < reltol

F = UpperTriangular(rand(Ty,n,n)+im*rand(Ty,n,n))
R = copy(F)
@time plyapcs!(acs,I,R,adj = true);
x = R'*R; @test norm(acs'*x+x*acs+F'*F)/norm(x)/norm(acs) < reltol

er = rand(Ty,n,n)
ar = er*ar
as, es = schur(ar,er)
ec = er+im*rand(Ty,n,n)
ac = ec*ac
acs, ecs = schur(ac,ec)

@time u = plyaps(as,es,br);
x = u*u'; @test norm(as*x*es'+es*x*as'+br*br')/norm(x)/norm(as) < reltol

@time u = plyaps(as,es,brw);
x = u*u'; @test norm(as*x*es'+es*x*as'+brw*brw')/norm(x)/norm(as) < reltol

@time u = plyaps(as',es',cr');
x = u'*u; @test norm(as'*x*es+es'*x*as+cr'*cr)/norm(x)/norm(as) < reltol

@time u = plyaps(as',es',crt');
x = u'*u; @test norm(as'*x*es+es'*x*as+crt'*crt)/norm(x)/norm(as) < reltol

@time u = plyaps(acs,ecs,bc);
x = u*u'; @test norm(acs*x*ecs'+ecs*x*acs'+bc*bc')/norm(x)/norm(as) < reltol

@time u = plyaps(acs',ecs',cc');
x = u'*u; @test norm(acs'*x*ecs+ecs'*x*acs+cc'*cc)/norm(x)/norm(as) < reltol

F = UpperTriangular(rand(Ty,n,n))
R = copy(F)
@time plyapcs!(as,es,R,adj = true);
x = R'*R; @test norm(as'*x*es+es'*x*as+F'*F)/norm(x)/norm(as) < reltol  #fails

F = UpperTriangular(rand(Ty,n,n))
R = copy(0*F)
@time plyapcs!(as,es,R,adj = true);
x = R'*R; @test norm(as'*x*es+es'*x*as) < reltol

F = UpperTriangular(rand(Ty,n,n))
R = copy(F)
@time plyapcs!(as,es,R,adj = false);
x = R*R'; @test norm(as*x*es'+es*x*as'+F*F')/norm(x)/norm(as) < reltol

F = UpperTriangular(rand(Ty,n,n)+im*rand(Ty,n,n))
R = copy(F)
@time plyapcs!(acs,ecs,R,adj = false);
x = R*R'; @test norm(acs*x*ecs'+ecs*x*acs'+F*F')/norm(x)/norm(acs) < reltol

F = UpperTriangular(rand(Ty,n,n)+im*rand(Ty,n,n))
R = copy(0*F)
@time plyapcs!(acs,ecs,R,adj = false);
x = R*R'; @test norm(acs*x*ecs'+ecs*x*acs') < reltol

F = UpperTriangular(rand(Ty,n,n)+im*rand(Ty,n,n))
R = copy(F)
@time plyapcs!(acs,ecs,R,adj = true);
x = R'*R; @test norm(acs'*x*ecs+ecs'*x*acs+F'*F)/norm(x)/norm(acs) < reltol

end

end
end

end


text:
__precompile__(true)

module Binary_BB_Trees

using Printf
using GeometryTypes: HomogenousMesh
using Rotations
using LinearAlgebra
using StaticArrays
using DataStructures
using Statistics
using NearestNeighbors
using NumericalTricks

include("box_types.jl")
include("mesh.jl")
include("util.jl")
include("vector_cache.jl")
include("tree_types.jl")
include("bb_intersection.jl")
include("blob_types.jl")
include("top_down.jl")
include("obb_construction.jl")
include("mesh_create_rot_sym.jl")
include("mesh_create_swept.jl")
include("extensions.jl")

export
    # box_types.jl
    BoundingBox,
    OBB,

    # mesh.jl
    Tri,
    Tet,
    eMesh,
    as_tet_eMesh,
    as_tri_eMesh,
    vertex_pos_for_tri_ind,
    vertex_pos_for_tet_ind,
    n_point,
    n_tri,
    n_tet,
    get_tri,
    get_tet,
    get_point,
    verify_mesh,
    eMesh_transform!,
    invert!,
    mesh_inplace_rekey!,
    mesh_remove_unused_points!,
    delete_triangles!,
    sub_div_mesh,
    mesh_repair!,
    output_eMesh_half_plane,
    output_eMesh_sphere,
    output_eMesh_box,

    # util.jl
    calc_min_max,
    calc_obb,
    crop_mesh,
    sort_so_big_ϵ_last,

    # vector_cache.jl
    VectorCache,
    expand!,
    returnNext,
    addCacheItem!,

    # tree_types.jl
    bin_BB_Tree,
    TT_Cache,
    update_TT_Cache!,
    is_leaf,
    is_not_leaf,
    treeDepth,
    leafNumber,
    tree_tree_intersect,
    extractData,

    # bb_intersection.jl
    BB_BB_intersect,

    # blob_types.jl
    blob,
    eMesh_to_tree,

    # top_down.jl

    # obb_construction.jl
    fit_tet_obb,
    fit_tri_obb,

    # mesh_create_rot_sym.jl
    obj_from_point_sequence,

    # mesh_create_swept.jl
    create_swept_mesh,
    f_swept_circle,
    f_swept_triv

    # extensions.jl

end # module


text:
using Revise; using BenchmarkTools; using FunctionWrappers:FunctionWrapper;

RheoFloat = Float64

struct TT
   G::FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}}
   J::FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}}

   Ga::FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}}
   Ja::FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}}
end

function ftest()
   # test base case
   p=(:a1, :a2)
   unpack_expr = Meta.parse(string(join(string.(p), ","), ", = a"))
   G = :(a1*t + a2)
   J = :(a2*t + a1)
   @eval return(TT(
      ((t, a) -> ($unpack_expr; $G)) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},
      ((t, a) -> ($unpack_expr; $J)) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},

      ((ta, a) -> ($unpack_expr; [$G for t in ta])) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
      ((ta, a) -> ($unpack_expr; [$J for t in ta])) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
      ))
end

macro mtest()
   # test base case
   p=(:a1, :a2)
   unpack_expr = Meta.parse(string(join(string.(p), ","), ", = a"))
   G = :(a1*t + a2)
   J = :(a2*t + a1)
   return quote TT(
      ((t, a) -> ($unpack_expr; $G)) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},
      ((t, a) -> ($unpack_expr; $J)) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},

      ((ta, a) -> ($unpack_expr; [$G for t in ta])) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
      ((ta, a) -> ($unpack_expr; [$J for t in ta])) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
     ) end
end

params = [1.0, 1.0];
testvector = Vector(0.0:0.01:100.0)

t1 = ftest()
t2 = @mtest

#println("Function")
@btime t1.Ga($testvector, $params) # 8.739 μs (3 allocations: 78.30 KiB)
@btime t2.Ga($testvector, $params) # 8.437 μs (3 allocations: 78.30 KiB)

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
"""
    bezier(line::LineString; resolution::Real=10000, sharpness::Real=0.85)

Take a LineString and returns a curved version
by applying a [Bezier spline](http://en.wikipedia.org/wiki/B%C3%A9zier_spline)
algorithm.
"""
function bezier(line::LineString; resolution::Real=10000., sharpness::Real=0.85)

    points::Vector{Point} = []

    for p in line.coordinates
        push!(points, Point(p))
    end

    spline = Spline(points=points, duration=resolution, sharpness=sharpness)

    coords = []
    for i in 1:10:spline.duration
        pos = Splines.pos(i)

        floor(i / 100) % 2 == 0 && push!(coords, [pos.coordinates[1], pos.coordinates[2]])

    end

    return LineString(coords)

end


text:
export tfullcontract, fullcontract

function fullcontract(contractor, region::T, tol = default_tolerance, maxiters = 100) where T
    for i in 1:maxiters
        region2 = region .∩ contractor(region)
        if diam(region2) < diam(region)
            region = region2
        else
            break
        end
        maximum(diam.(region)) < tol && break
    end
    return region
end

tfullcontract(contractor, regions, tol = default_tolerance) =
    qmap(r -> fullcontract(contractor, r, tol), regions)


text:
# some utility functions.
using SeisIO, SeisNoise, Printf, DelimitedFiles,Logging

export readfilelist, writesac_rich,compile_stationinfo,update_siteinfo!

const sac_nul_f = -12345.0f0
const sac_nul_i = Int32(-12345)
const sac_nul_start = 0x2d
const sac_nul_Int8 = UInt8[0x31, 0x32, 0x33, 0x34, 0x35]
const μs = 1.0e-6

function readfilelist(infilelist::String)
    outlist = String[];
    open(infilelist) do f
        for ln in eachline(f)
          push!(outlist,ln)
        end
    end
    return(outlist)
end

function compile_stationinfo(siteinfofile::String;separator::AbstractChar=' ',header::Bool=false)
    siteinfotemp=readdlm(siteinfofile,',',header=header)
    siteinfodict=Dict()
    prevsiteid=""
    for i = 1:size(siteinfotemp[1],1)
        siteid="$(siteinfotemp[1][i,1])"*"."*"$(siteinfotemp[1][i,2])"
        if siteid == prevsiteid
            continue
        end
      siteinfodict["$(siteinfotemp[1][i,1])"*"."*"$(siteinfotemp[1][i,2])"] = siteinfotemp[1][i,:]
      prevsiteid=siteid
    end
    return(siteinfodict)
end
#utility function
function update_siteinfo!(S::SeisData,siteid::String,sitedict::Dict)
    for i = 1:size(S)[1]
        SC=SeisChannel()
        SC=S[i]
        SC.id=siteid
        SC.name=siteid
        siteid_temp=join(split(siteid,".")[1:2],".")
        setfield!(SC.loc,:lat,Float64(sitedict[siteid_temp][3]))
        setfield!(SC.loc,:lon,Float64(sitedict[siteid_temp][4]))
        setfield!(SC.loc,:el,Float64(sitedict[siteid_temp][5]))
        # = C.loc
        # TO-DO get misc values
        SC.misc["stla"] = sitedict[siteid_temp][3]
        SC.misc["stlo"] = sitedict[siteid_temp][4]
        SC.misc["stel"] = sitedict[siteid_temp][5]
        S[i]=SC
    end
    return nothing
end

# this function and fill_sac_rich(), and writesac() were modified from those in SeisIO
# I modified it to save correlation data, which has two station information in the name.
function write_sac_file(fname::String, fv::Array{Float32,1}, iv::Array{Int32,1}, cv::Array{UInt8,1}, x::Array{Float32,1};
                        t=[Float32(0)]::Array{Float32,1}, ts=true::Bool)
  f = open(fname, "w")
  write(f, fv)
  write(f, iv)
  write(f, cv)
  write(f, x)
  if ts
    write(f, t)
  end
  close(f)
  return
end
#rich header information saved in misc field. It handles CorrData headers. If not CorrData, it does the same thing as fill_sac() in SeisIO package
function fill_sac_rich(S::GphysChannel, ts::Bool, leven::Bool;dtype::DataType=CorrData)
  fv = sac_nul_f.*ones(Float32, 70)
  iv = sac_nul_i.*ones(Int32, 40)
  cv = repeat(vcat(sac_nul_Int8, sac_nul_start, 0x20, 0x20), 24)
  cv[17:22] .= 0x20

  # Ints
  T = getfield(S, :t)
  tt = Int32[Base.parse(Int32, i) for i in split(string(u2d(T[1,2]*μs)),r"[\.\:T\-]")]
  length(tt) == 6 && append!(tt, zero(Int32))
  y = tt[1]
  j = Int32(md2j(y, tt[2], tt[3]))
  iv[1:6] = prepend!(tt[4:7], [y, j])
  iv[7] = Int32(6)
  iv[10] = Int32(length(S.x))
  iv[16] = Int32(ts ? 4 : 1)
  iv[36] = Int32(leven ? 1 : 0)

  # Floats
  dt = 1.0/S.fs
  fv[1] = Float32(dt)
  fv[4] = Float32(S.gain)

  y_s = string(y); y_s="0"^(4-length(y_s))*y_s
  j_s = string(j); j_s="0"^(3-length(j_s))*j_s
  h_s = string(tt[4]); h_s="0"^(2-length(h_s))*h_s
  m_s = string(tt[5]); m_s="0"^(2-length(m_s))*m_s
  s_s = string(tt[6]); s_s="0"^(2-length(s_s))*s_s
  ms_s = string(tt[7]); ms_s="0"^(3-length(ms_s))*ms_s

  if dtype == CorrData
    #get station and event locations from misc fiels, assuming they are saved.
    fv[36]=Float32(S.misc["evla"])
    fv[37]=Float32(S.misc["evlo"])
    fv[39]=Float32(S.misc["evdp"])
    fv[32]=Float32(S.misc["stla"])
    fv[33]=Float32(S.misc["stlo"])
    fv[34]=Float32(S.misc["stel"])
    fv[6]=Float32(S.misc["b"])
    fv[7]=Float32(S.misc["e"])

    # Chars (ugh...)
    id_temp = String.(split(S.name,'.')) #8 entries for CorrData naming convention: two stations.
    evname = id_temp[1]*"."*id_temp[2]*"."*id_temp[3]*"."*id_temp[4] #
    # println(evname)
    id = [evname; id_temp[5:8]]
    #change component name to CorrData.comp
    id[5] = S.misc["comp"]
    #make the first station as event
    ci = [9, 169, 1, 25, 161]
    Lc = [16, 8, 8, 8, 8]
    for i = 1:length(id)
      if !isempty(id[i])
        L_max = Lc[i]
        si = ci[i]
        ei = ci[i] + L_max - 1
        s = codeunits(id[i])
        Ls = length(s)
        L = min(Ls, L_max)
        copyto!(cv, si, s, 1, L)
        if L < L_max
          cv[si+L:ei] .= 0x20
        end
      end
    end
    #other character headers
    kh= String[]
    push!(kh,S.misc["corr_type"])
    push!(kh,S.misc["whitened"])
    push!(kh,S.misc["time_norm"])

    ci = [137, 145, 153]
    Lc = [8, 8, 8]
    for i = 1:length(kh)
      if !isempty(kh[i])
        L_max = Lc[i]
        si = ci[i]
        ei = ci[i] + L_max - 1
        s = codeunits(kh[i])
        Ls = length(s)
        L = min(Ls, L_max)
        copyto!(cv, si, s, 1, L)
        if L < L_max
          cv[si+L:ei] .= 0x20
        end
      end
    end
    #other header values.
    fv[2]=Float32(S.misc["depmin"])
    fv[3]=Float32(S.misc["depmax"])
    fv[57]=Float32(S.misc["depmen"])
    fv[41]=Float32(S.misc["freqmin"])
    fv[42]=Float32(S.misc["freqmax"])
    fv[43]=Float32(S.misc["cc_len"])
    fv[44]=Float32(S.misc["cc_step"])

    # Assign a filename
    fname = join([y_s, j_s,h_s, m_s, s_s, ms_s, id[1], id[2], id[3], id[4],id_temp[8],id[5],"sac"],'.')
  else
    fv[6] = rem(T[1,2], 1000)*1.0f-3
    fv[7] = Float32(dt*length(S.x) + sum(T[2:end,2])*μs)
    if !isempty(S.loc)
      loc = getfield(S, :loc)
      # println(typeof(loc))
      if typeof(loc) == GeoLoc
        # println("getting locations.")
        fv[32] = Float32(getfield(loc, :lat))
        fv[33] = Float32(getfield(loc, :lon))
        fv[34] = Float32(getfield(loc, :el))
        fv[58] = Float32(getfield(loc, :az))
        fv[59] = Float32(getfield(loc, :inc))
      end
    # else
    #   @warn "loc field is empty"
    end

    # Chars (ugh...)
    id_temp = String.(split(S.id,'.')) #8 entries for CorrData naming convention: two stations.
    id = id_temp[1:4]
    #make the first station as event
    ci = [169, 1, 25, 161]
    Lc = [8, 8, 8, 8]
    for i = 1:length(id)
      if !isempty(id[i])
        L_max = Lc[i]
        si = ci[i]
        ei = ci[i] + L_max - 1
        s = codeunits(id[i])
        Ls = length(s)
        L = min(Ls, L_max)
        copyto!(cv, si, s, 1, L)
        if L < L_max
          cv[si+L:ei] .= 0x20
        end
      end
    end

    # Assign a filename
    fname = join([y_s, j_s,h_s, m_s, s_s, ms_s, id[1], id[2], id[3], id[4],"sac"],'.')
  end

  return (fv, iv, cv, fname)
end

"""
    writesac(S::Union{SeisData}[; ts=false, v=0])

Write all data in SeisData structure `S` to auto-generated SAC files.
"""
function writesac_rich(S::GphysData; datatype::DataType=CorrData, outdir::String="", outfile::String="", ts::Bool=false, v::Int64=0)
  if ts
    ift = Int32(4); leven = false
  else
    ift = Int32(1); leven = true
  end
  tdata = Array{Float32}(undef, 0)
  N     = S.n
  for i = 1:N
    T = getindex(S, i)
    b = T.t[1,2]
    dt = 1.0/T.fs
    (fv, iv, cv, fname) = fill_sac_rich(T, ts, leven,dtype=datatype)

    # Data
    x = eltype(T.x) == Float32 ? getfield(T, :x) : map(Float32, T.x)
    ts && (tdata = map(Float32, μs*(t_expand(T.t, dt) .- b)))

    # Write to file
    # println(fname)
    if isempty(outfile)
      if N > 1
        fname2 =fname[1:end-4]*"_"*string(i)*".sac"
      else
        fname2 = fname
      end
    else
      fname2 = outfile
    end
    # println(fname2)
    if length(outdir) > 0
      mkpath(outdir)
      write_sac_file(outdir*"/"*fname2, fv, iv, cv, x, t=tdata, ts=ts)
    else
      write_sac_file(fname2, fv, iv, cv, x, t=tdata, ts=ts)
    end
    v > 0  && @printf(stdout, "%s: Wrote file %s from channel %i\n", string(now()), fname2, i)
  end
end


text:
"""
Data types and functions support `Configuration`s, which describe the participants
in a Paxos cluster.
"""
module Configurations

export isquorum,
  quorumSize,
  NoQuorumAvailableException,
  Configuration,
  configuration,
  memberAddress,
  memberAddresses,
  addMember,
  removeMember

using ..Nodes

struct OrdinaryConfiguration
  """
  Dictionary mapping IDs of members to their addresses
  """
  members::Dict{NodeID,Any}
end

"""
A transition configuration (inspired by Raft) supports rollout of a
new configuration, with potentially differences in membership
"""
struct TransitionConfiguration
  old::OrdinaryConfiguration
  new::OrdinaryConfiguration
end

"""
A configuration is the general type for maintaining a list of menbers
and addreses where they can be reached over transport
"""
const Configuration = Union{OrdinaryConfiguration,TransitionConfiguration}

"""
Return the address for the specfied node, or error if no such address
"""
function memberAddress(cfg::OrdinaryConfiguration, memberID::NodeID)
  cfg.members[memberID]
end

function memberAddress(cfg::TransitionConfiguration, memberID::NodeID)
  if haskey(cfg.old, memberID)
    memberAddress(cfg.old, memberID)
  else
    memberAddress(cfg.new, memberID)
  end
end

"""
Return the addresses of all membes in the configuration
"""
function memberAddresses(cfg::OrdinaryConfiguration)
  Set(keys(cfg.members))
end

function memberAddresses(cfg::TransitionConfiguration)
  union(memberAddresses(cfg.old), memberAddresses(cfg.new))
end

"""
Return the minimum size for a majority of members in the configuration
"""
function quorumSize(cfg::OrdinaryConfiguration)
  quorumSize(length(cfg.members))
end

"""
For `Integer` values, return the required number of participants for a majority
"""
function quorumSize(len::Integer)
  convert(Integer, if iseven(len)
    (len / 2) + 1
  else
    ceil(len / 2)
  end)
end

"""
Exception thrown when no quorum is available to respond to Paxos messages
"""
struct NoQuorumAvailableException <: Exception end

"""
Return true if the list of members represents a quorum for the configuration
"""
function isquorum(cfg::OrdinaryConfiguration, members::Set{NodeID})
  count = 0
  for member in members
    if haskey(cfg.members, member)
      count += 1
    end
  end
  count >= quorumSize(cfg)
end

"""
During transitions, eturn true if the list of members represents a quorum in
both old and new configurations
"""
function isquorum(cfg::TransitionConfiguration, members::Set{NodeID})
  isquorum(cfg.old, members) && isquorum(cfg.new, members)
end

function Base.isempty(cfg::OrdinaryConfiguration)
  isempty(cfg.members)
end

function Base.isempty(cfg::TransitionConfiguration)
  isempty(cfg.old) && isempty(cfg.new)
end

"""
Create a new, empty `OrdinaryConfiguration`
"""
function configuration(members::Dict{NodeID,Any} = Dict{NodeID,Any}())
  OrdinaryConfiguration(members)
end

function addMember(cfg::OrdinaryConfiguration, id::NodeID, address)
  TransitionConfiguration(cfg, configuration(merge(cfg.members, Dict(id => address))))
end

function addMember(cfg::TransitionConfiguration, id::NodeID, address)
  TransitionConfiguration(
    cfg.old,
    configuration(merge(cfg.new.members, Dict(id => address))),
  )
end

function removeMember(cfg::OrdinaryConfiguration, id::NodeID)
  TransitionConfiguration(cfg, configuration(delete!(Dict(cfg.members), id)))
end

function removeMember(cfg::TransitionConfiguration, id::NodeID)
  TransitionConfiguration(cfg, configuration(delete!(Dict(cfg.new.members), id)))
end

end


text:
abstract type AbstractStochasticProcess end
abstract type AbstractMarkovChain <: AbstractStochasticProcess end

"""
    MarkovChain

A discrete-time, discrete-space, finite-space Markov chain.

##### Fields
- `transitionmatrix::Matrix{Float64}` : stochastic matrix in which element `[i,j]` denotes
the probability that the Markov chain transitions from state `i` to state `j`.
- `initialdist::Vector{Float64}` : stochastic vector in which element `i` denotes the
probability that the Markov chain's initial state is `i`.

##### Complexity (for the inner constructor)
- Time complexity: `O(n^2)`, where `n` denotes the number of rows (or columns) in a square
matrix.
"""
struct MarkovChain <: AbstractMarkovChain
    transitionmatrix::Matrix{Float64}
    initialdist::Vector{Float64}
    # ensure that the transition matrix and the initial distribution are stochastic
    function MarkovChain(transitionmatrix::AbstractMatrix{<:Real},
                         initialdist::AbstractVector{<:Real})
        size(transitionmatrix)[1] != length(initialdist) &&
            ArgumentError("Uh-oh! The transition matrix must have the same number of rows as
                          the initial distribution has elements.") |> throw
        !isstochastic(initialdist) || !isstochastic(transitionmatrix) &&
            ArgumentError("Oops! The transition matrix and the initial distribution must be
                          stochastic.") |> throw
        return new(transitionmatrix, initialdist)
    end
end

"""
    Base.show(io::IO, markovchain::MarkovChain) -> nothing

Show a `MarkovChain` in a human-friendly manner.
"""
function Base.show(io::IO, markovchain::MarkovChain)
    println(io, typeof(markovchain))
    println(io, " Initial distribution:")
    println(io, "  ", markovchain.initialdist)
    println(io, " Transition matrix:")
    println(io, "  ", markovchain.transitionmatrix)
    return nothing
end

"""
    MarkovChain(transitionmatrix::AbstractMatrix{<:Real},
                initialdist::AbstractMatrix{<:Real}) -> MarkovChain

Construct a `MarkovChain` when the initial distribution `initialdist` is specified as either
a one-row or a one-column matrix.

##### Complexity
- Time complexity: `O(n^2)`, where `n` denotes the number of rows (or columns) in a square
matrix.
"""
function MarkovChain(transitionmatrix::AbstractMatrix{<:Real},
                     initialdist::AbstractMatrix{<:Real})
    nrows, ncols = size(initialdist)[1], size(initialdist)[2]
    if nrows == 1
        # make initialdist a vector
        initialdist = initialdist[1,:]
    elseif ncols == 1
        # make initialdist a vector
        initialdist = initialdist[:,1]
    else
        throw(ArgumentError("The initial distribution must be a vector."))
    end
    # multiple dispatch
    return MarkovChain(transitionmatrix, initialdist)
end

"""
    isstochastic(matrix::AbstractMatrix, tol::Float64 = 1e-5) -> Bool

Decide if `matrix` is stochastic.

`tol` is the tolerance associated with the row sums. We consider the row sums to be valid if
and only if they are in the interval `[1-tol, 1+tol]`. This is important when the element
type is `Float64`, for example.

##### Complexity
- Time complexity: `O(n^2)`, where `n` denotes the number of rows (or columns) in a square
matrix.
"""
function isstochastic(matrix::AbstractMatrix, tol::Float64 = 1e-5)
    nrows, ncols = size(matrix)[1], size(matrix)[2]
    # if matrix has only one row
    if nrows == 1
        vec = matrix[1,:]::Vector
        # multiple dispatch
        return isstochastic(vec, tol)
    # if matrix has only one column
    elseif ncols == 1
        vec = matrix[:,1]::Vector
        # multiple dispatch
        return isstochastic(vec, tol)
    # if matrix is non-square
    elseif nrows != ncols
        return false
    end
    # check if elements are in [0,1]
    for element in matrix
        if element < 0 || element > 1
            return false
        end
    end
    # check row sums
    for row in eachrow(matrix)
        rowsum = sum(row)
        if rowsum < 1 - tol || rowsum > 1 + tol
            return false
        end
    end
    return true  # this line is reached only if all above checks pass
end

"""
    isstochastic(vector::AbstractVector, tol::Float64 = 1e-5) -> Bool

Decide if `vector` is stochastic.

`tol` is the tolerance associated with the sum over `vector`'s elements. We consider the sum
to be valid if and only if it is in the interval `[1-tol, 1+tol]`. This is important when
the element type is `Float64`, for example.

##### Complexity
- Time complexity: `O(n)`.
"""
function isstochastic(vector::AbstractVector, tol::Float64 = 1e-5)
    # check if elements are in [0,1]
    for element in vector
        if element < 0 || element > 1
            return false
        end
    end
    vecsum = sum(vector)
    if vecsum < 1 - tol || vecsum > 1 + tol
        return false
    end
    return true  # this line is reached only if all above checks pass
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
# included in module Wut

struct ScalarEncoder

	w::Number            # width of the contiguous 1 bits, must be odd
	minval::Number       # minimum value of input signal
	maxval::Number       # maximum value of input signal (strictly less when periodic = true)
	periodic::Bool       # input values wrap around and maxval input is mapped onto minval
	n::Number            # output bit width (must be greater than w)
	radius::Number       # input value difference > radius do not overlap bits
	resolution::Number   # input value difference >= resolution do not give identical bits
	name::AbstractString # optional descriptive string
	verbosity::Number
	clipInput::Bool      # Non-periodic inputs are clipped to minval / maxval
	forced::Bool         # skip some safety checks
    padding::Number

    function ScalarEncoder(;w=0,minval=0,maxval=0,periodic=false,
				n=0,radius=0,resolution=1,name="",
				verbosity=0,clipInput=false,forced=false,padding=0)
        res=(maxval-minval)/n
		new(w,minval,maxval,periodic,n,w*res,res,name,verbosity,clipInput,forced,periodic?0:(w-1)/2)

end

end
export ScalarEncoder

#struct BitPat
#	e::Unsigned
#	b::BitArray{1}
#	BitPat(e,b) = new(e,BitArray(e))
#end
#export BitPat

function encode(e::ScalarEncoder, n::Number)
	encodeIntoArray(e, n, BitPat(e.n,n))
end
export encode

#function Base.show(io::IO, m::BitPat)
#	print(io,"[")
#	for (i,v) in enumerate(m.b)
#		v ? print(io," 1") : print(io," 0")
#	end
#	print(io," ]")
#end

function getWidth(e::ScalarEncoder)
	e.n
end
export getWidth

function getDescription(e::ScalarEncoder)
	e.name
end
export getDescription

function getBucketIndices(e::ScalarEncoder)

end
export getBucketIndices

function encodeIntoArray(e::ScalarEncoder,n::Number,b::BitPat;learn=true)
	fill!(b.b,false)
	t1 = e.periodic ? n % e.maxval : n
	t2 = t1 > e.maxval ? e.maxval : t1
	t = t2 < e.minval ? e.minval : t2

    c = convert(Int64,round((((t - e.minval) + e.resolution/2) / e.resolution )))# + e.padding))

#	i=convert(Int64,round(t*e.n/e.maxval))
    d=convert(Int64,(e.w-1)/2)
    for i = c-d : c+d
	  if i <= e.n && i > 0
        b.b[i]=true
      end
    end

	b
end
export encodeIntoArray

function decode(e::ScalarEncoder)

end
export decode

function getBucketValues(e::ScalarEncoder)

end
export getBucketValues

function getBucketInfo(e::ScalarEncoder)

end
export getBucketInfo

function topDownCompute(e::ScalarEncoder)

end
export topDownCompute

function closenessScores(e::ScalarEncoder)

end
export closenessScores

# end # module


text:
#=
Functions that are derived from Base or other packages
=#
function Base.show(io::IO, m::MonitorProperties)
    println(io, "name: ", m.name)
    println(io, "physicalsize: ",  m.physicalsize[1], "x", m.physicalsize[2])
    println(io, "resolution: ", m.videomode.width, "x", m.videomode.height)
    println(io, "dpi: ", m.dpi[1], "x", m.dpi[2])
end
function Base.show(io::IO, m::Screen)
    println(io, "name: ", m.name)
    println(io, "children: ", length(m.children))
    println(io, "Inputs:")
    for (key, value) in m.inputs
        println(io, "  ", key, " => ", typeof(value))
    end
end

"""
mouse position is in coorditnates relative to `screen`
"""
function GeometryTypes.isinside(screen::Screen, mpos)
    isinside(zeroposition(value(screen.area)), mpos...) || return false
    for s in screen.children
        # if inside any children, it's not inside screen
        isinside(value(s.area), mpos...) && return false
    end
    true
end

"""
Args: `screens_mpos` -> Tuple{Vector{Screen}, Vec{2,T}}
"""
function isoutside(screens_mouseposition)
    screens, mpos = screens_mouseposition
    for screen in screens
        isinside(screen, mpos) && return false
    end
    true
end

"""
Returns the monitor resolution of the primary monitor.
"""
function primarymonitorresolution()
    props = MonitorProperties(GLFW.GetPrimaryMonitor())
    w,h = props.videomode.width, props.videomode.height
    Vec(Int(w),Int(h))
end

"""
Create a new rectangle with x,y == 0,0 while taking the widths from the original
Rectangle
"""
zeroposition{T}(r::SimpleRectangle{T}) = SimpleRectangle(zero(T), zero(T), r.w, r.h)


text:
"""
    DRR_JJnnp(n,np,a,j,IntTable_K=IntTable_K_serial,tabSMARes=tabSMARes_serial)

Computation of DRR_JJnn, contribution of the (n,np) resonance number to DRR.

Using log-sampling in the integration for better interpolation.

# Arguments
- `n   ::Int64`  : first resonance number of the orbit.
- `np  ::Int64`  : second resonance number of the orbit.
- `a   ::Float64`: semi-major axis of the orbit.
- `j   ::Float64`: reduced angular momentum of the orbit.
- `IntTable_K=IntTable_K_serial:: IntTable`: Structure which contains the tables needed for the computation of this integral.
- `tabSMARes=tabSMARes_serial`: Static container for the resonant lines, contains the edge of the domain [amin,amax] of the resonant line.
"""
function DRR_JJnnp(n::Int64,np::Int64,a::Float64,j::Float64,IntTable_K::IntTable=IntTable_K_serial,tabSMARes::Array{Float64,2}=tabSMARes_serial)
    #####
    omega = n*nup(a,j) # Value of the resonant frequency
    #####
    ResPoints = tabSMARes!(n,np,a,j,tabSMARes) # Trying to compute the resonant line
    #####
    if (ResPoints == 0)
        return 0.0 # If the resonant line is empty, we return 0.0 for the diffusion coefficient
    end
    #####
    # First resonance line
    apmin1, apmax1 = tabSMARes[1,1], tabSMARes[1,2] # Boundary of the SMA-region where the resonant line lies
    dalphap1 = log(apmax1/apmin1)/(nbResPoints) # Step distance in log-space

    if (ResPoints == 2) # If second resonance line
        apmin2, apmax2 = tabSMARes[2,1], tabSMARes[2,2] # Boundary of the SMA-region where the resonant line lies
        dalphap2 = log(apmax2/apmin2)/(nbResPoints) # Step distance in log-space
    end
    #####
    res = 0.0 # Initialisation of the result
    #####
    for alphap1 = range(log(apmin1),length=nbResPoints,log(apmax1)-dalphap1) # Uniform logarithmic spacing of the SMA !! ATTENTION, we do not go through the last element, as we will sum in the middle of the interval
        ap = exp(alphap1 + 0.5*dalphap1) # Central position of the SMA considered
        jp = getjRes(np,ap,omega) # Finding the location of the resonant jp
        res += dalphap1*ap*FtotBath(ap,jp)*AnnpSQ(n,np,a,j,ap,jp,IntTable_K)/(abs(dnupdj(ap,jp))) # Local contribution !! ATTENTION, not to forget the ap
    end
    if (ResPoints == 2) # If two resonance lines
        for alphap2 = range(log(apmin2),length=nbResPoints,log(apmax2)-dalphap2) # Uniform logarithmic spacing of the SMA !! ATTENTION, we do not go through the last element, as we will sum in the middle of the interval
            ap = exp(alphap2 + 0.5*dalphap2) # Central position of the SMA considered
            jp = getjRes(np,ap,omega) # Finding the location of the resonant jp
            res += dalphap2*ap*FtotBath(ap,jp)*AnnpSQ(n,np,a,j,ap,jp,IntTable_K)/(abs(dnupdj(ap,jp))) # Local contribution !! ATTENTION, not to forget the ap
        end
    end

    #####
    res *= (n^(2))/(abs(np)) # Adding the prefactor in n^2/|np|
    return res # Output
end

"""
    DRR_JJ(a,j,IntTable_K=IntTable_K_serial,tabSMARes=tabSMARes_serial)

Computation of DRR_JJ the scalar resonant diffusion coefficient in JJ-coordinates.

# Arguments
- `a   ::Float64`: semi-major axis of the orbit.
- `j   ::Float64`: reduced angular momentum of the orbit.
- `IntTable_K=IntTable_K_serial:: IntTable`: Structure which contains the tables needed for the computation of this integral.
- `tabSMARes=tabSMARes_serial`: Static container for the resonant lines, contains the edge of the domain [amin,amax] of the resonant line.

"""
function DRR_JJ(a::Float64,j::Float64,IntTable_K::IntTable=IntTable_K_serial,tabSMARes::Array{Float64,2}=tabSMARes_serial)
    res = 0.0 # Initialisation of the result
    #####
    for iPair=1:nbResPairs # Loop over the resonant pairs
        n, np = tabResPairs[1,iPair], tabResPairs[2,iPair] # Value of the considered resonance (n,np)
        res += DRR_JJnnp(n,np,a,j,IntTable_K,tabSMARes) # Adding the contribution from the resonance pair (n,np)
    end
    #####
    res *= 4.0*pi*G^(2) # Adding the prefactor 4piG^2
    return res # Output
end

"""
    DRR_jj(a,j,IntTable_K=IntTable_K_serial,tabSMARes=tabSMARes_serial)

Computation of DRR_jj the (normalized) scalar resonant diffusion coefficient in jj-coordinates.

# Arguments
- `a   ::Float64`: semi-major axis of the orbit.
- `j   ::Float64`: reduced angular momentum of the orbit.
- `IntTable_K=IntTable_K_serial:: IntTable`: Structure which contains the tables needed for the computation of this integral.
- `tabSMARes=tabSMARes_serial`: Static container for the resonant lines, contains the edge of the domain [amin,amax] of the resonant line.

"""
function DRR_jj(a::Float64,j::Float64,IntTable_K::IntTable=IntTable_K_serial,tabSMARes::Array{Float64,2}=tabSMARes_serial)
    return DRR_JJ(a,j,IntTable_K,tabSMARes)/((Jc(a))^(2)) # Output
end


text:
export PowerMeasure
using Base: eltype

"""
    struct PowerMeasure{M,N}
        μ::M
        size::NTuple{N,Int}
    end

A power measure is a product of a measure with itself. The number of elements in
the product determines the dimensionality of the resulting support.

Note that power measures are only well-defined for integer powers.

The nth power of a measure μ can be written μ^x.
"""
struct PowerMeasure{M,N}
    μ::M
    size::NTuple{N,Int}
end

function Base.:^(μ::Measure{X}, n::Integer) where {X}
    components = ntuple(i -> μ, n)
    ProductMeasure(components...)
end


text:
using GLVisualize, GeometryTypes, FileIO
using GLAbstraction, Colors, Reactive
if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(0f0:0.1:1f0)
end

description = """
Example how to animate a 2D barplot.
"""

primitive = SimpleRectangle(0f0, -0.5f0, 1f0, 1f0)
scalars = rand(10f0:0.01f0:200f0, 200)

function interpolate(a, scalars, t)
    [ae + ((be - ae) * t) for (ae, be) in zip(a, scalars)]
end
t = const_lift(*, timesignal, 10f0)
v0 = (scalars, rand(10f0:0.01f0:200f0, 200), scalars)
interpolated = foldp(v0, t) do v0_v1_ip, td
    v0, v1, ip = v0_v1_ip
    pol = td % 1
    if isapprox(pol, 0.0)
        v0 = v1
        v1 = map(x-> rand(linspace(-50f0, 60f0, 100)), v0)
    end
    v0, v1, interpolate(v0, v1, pol)
end
b_sig = map(last, interpolated)
bars = visualize(
    (RECTANGLE, b_sig),
    ranges = linspace(0, 600, 10),
    color_norm = Vec2f0(-40, 200),
    color_map = GLVisualize.default(Vector{RGBA})
)

_view(bars, window, camera = :orthographic_pixel)
if !isdefined(:runtests)
    renderloop(window)
end


text:
using ArgParse
using BioSequences
using CodecZlib
using FASTX

## Command Line Interface ##

cli = ArgParseSettings()
@add_arg_table cli begin
	"--window"
		help = "window size"
		arg_type = Int
		default = 100
	"fasta"
		help = "fasta file, required"
		required = true
end
arg = parse_args(ARGS, cli)

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
module FieldArithmetic

export Ring, RingElem, CRing, CRingElem, Field, FieldElem, Zmod, ZmodElem, PolyRing, PolyRingElem,
QuotientRing, QuotientRingElem, Fpsquare, order, squareroot, cuberoot, cuberootOfOne, allcuberoots,
allElements

##########################################
#         For overwriting
##########################################

import Base: +,-,*,/,^,==,inv,show,zeros,gcd,gcdx,iszero,isone,zero,one
import Primes: isprime

##########################################
# Abstract types for defining domains
##########################################

abstract type Ring end

abstract type RingElem end

abstract type CRing <: Ring end

abstract type CRingElem <: RingElem end

##########################################
#        Generic helper functions
##########################################

function sameRing(x::RingElem, y::RingElem)
    return x.ring == y.ring
end

function zeros(R::T, len::Integer) where T <: Ring
    return [zero(R) for i in 1:len]
end

##########################################
#
#         Ring of Integers mod n
#
##########################################

struct Zmod{T <: Integer} <: CRing
    n::T
end

struct ZmodElem{T <: Integer} <: CRingElem
    d::T
    ring::Zmod{T}
end

# Ring-Constructor
function Zmod(n::T) where T <: Integer
    return Zmod{T}(n)
end

# Element-Constructor
function (R::Zmod{T})(a::Integer) where T <: Integer
    a = T(a)
    n = R.n
    d = a%n
    if d < 0
        d += n
    end
    return ZmodElem{T}(d, R)
end

###### Printing ######
function show(io::IO, R::Zmod{T}) where T <: Integer
    print(io, "Z_$(R.n)")
end

function show(io::IO, x::ZmodElem{T}) where T <: Integer
    print(io, "$(x.d)")
end

###### Helper functions ######

function zero(R::Zmod{T}) where T <: Integer
    return R(T(0))
end

function one(R::Zmod{T}) where T <: Integer
    return R(T(1))
end

function iszero(x::ZmodElem{T}) where T <: Integer
    return x.d == 0
end

function isone(x::ZmodElem{T}) where T <: Integer
    return x.d == 1
end

function order(R::Zmod{T}) where T <: Integer
    return R.n
end

function allElements(R::Zmod{T}) where T <: Integer
    return [R(i) for i in 0:R.n-1]
end

###### Basic operators #######
function -(x::ZmodElem{T}) where T <: Integer
    if x.d == 0
        return x
    else
        R = x.ring
        return R(R.n - x.d)
    end
end

function +(x::ZmodElem{T}, y::ZmodElem{T}) where T <: Integer
    if !sameRing(x, y)
        error("Cannot add elements of different rings")
    end
    z = x.d + y.d
    R = x.ring
    return R(z)
end

function -(x::ZmodElem{T}, y::ZmodElem{T}) where T <: Integer
    return x+(-y)
end

function *(x::ZmodElem{T}, y::ZmodElem{T}) where T <: Integer
    if !sameRing(x, y)
        error("Cannot multiply elements of different rings")
    end
    z = x.d * y.d
    R = x.ring
    return R(z)
end

function *(x::Integer, y::ZmodElem{T}) where T <: Integer
    return y.ring(x)*y
end

function *(x::ZmodElem{T}, y::Integer) where T <: Integer
    return x.ring(y)*x
end

function inv(x::ZmodElem{T}) where T <: Integer
    R = x.ring
    g, a, b = gcdx(x.d, R.n)
    g != 1 && error("Element $x is not invertible in $(x.ring)")
    return R(a)
end

#exponensiation uses the simple square and multiply algorithm
function ^(x::ZmodElem{T}, y::Integer) where T <: Integer
    if x.d == 0 || y == 1
        return x
    end
    R = x.ring
    if y == 0
        return R(T(1))
    end
    if y < 0
        x = inv(x)
        y = -y
    end
    # Should include some reduction of y modulo order if order is known
    bit = T(1) << (ndigits(y, base = 2) - 1)
    z = x
    bit >>= 1
    while bit != 0
        z = z*z
        if (bit & y) != 0
            z *= x
        end
        bit >>=1
    end
    return z
end

function ==(x::ZmodElem{T}, y::ZmodElem{T}) where T <: Integer
    return sameRing(x, y) && (x.d == y.d)
end

function ==(R1::Zmod{T}, R2::Zmod{T}) where T <: Integer
    return R1.n == R2.n
end

##########################################
#
#         Polynomial Ring
#
##########################################

struct PolyRing{T <: CRing} <: CRing
    parent::T
    variable::String
end

struct PolyRingElem{T <: CRingElem, S <: CRing} <: CRingElem
    d::Array{T}
    ring::PolyRing{S}
end

# Ring-Constructor
function PolyRing(R::T,var::String) where T <: CRing
    return PolyRing{T}(R,var)
end

function PolyRing(R::T) where T <: CRing
    return PolyRing{T}(R,"X")
end

# Element-Constructor
function (R::PolyRing{S})(a::Array{T}) where {T <: CRingElem, S <: CRing}
    while length(a) > 0
        if !(iszero(a[length(a)]))
            return PolyRingElem{T,S}(a, R)
        end
        pop!(a)
    end
    return PolyRingElem{T,S}([zero(R.parent)], R)
end

function (R::PolyRing{Zmod{T}})(a::Array{S}) where {T <: Integer, S <: Integer}
    Z = R.parent
    a = [Z(i) for i in a]
    return R(a)
end

###### Printing ######
function show(io::IO, R::PolyRing{T}) where T <: CRing
    print(io, "$(R.parent)[$(R.variable)]")
end

function show(io::IO, x::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    s = "$(x.d[1])"
    l = length(x.d)
    R = x.ring
    if l > 1
        s *= " + $(x.d[2])*$(R.variable)"
        if l > 2
            for i in 3:l
                s *= " + $(x.d[i])*$(R.variable)^$(i-1)"
            end
        end
    end
    print(io, s)
end

###### Useful helper functions #######
function zero(R::PolyRing{T}) where T <: CRing
    return R([zero(R.parent)])
end

function one(R::PolyRing{T}) where T <: CRing
    return R([one(R.parent)])
end

function iszero(a::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    return a==zero(a.ring)
end

function isone(a::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    return a==one(a.ring)
end

function degree(a::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    if iszero(a)
        return -1
    else
        return length(a.d)-1
    end
end

function LC(a::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    return a.d[length(a.d)]
end

function monic(a::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    x = inv(LC(a))
    return a*x
end

###### Basic operators #######
function -(a::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    d = deepcopy(a.d)
    l = length(d)
    R = a.ring
    z = zeros(R.parent, l)
    for i in 1:l
        z[i] = -d[i]
    end
    return R(z)
end

function +(a::PolyRingElem{T,S}, b::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    if !sameRing(a, b)
        error("Cannot add elements from different rings")
    end
    R = a.ring
    d1 = deepcopy(a.d)
    d2 = deepcopy(b.d)
    lena = length(d1)
    lenb = length(d2)
    lenz = max(lena,lenb)
    lenmin = min(lena,lenb)
    z = zeros(R.parent, lenz)
    for i in 1:lenmin
        z[i] = d1[i] + d2[i]
    end
    if lena > lenmin
        for i in lenmin+1:lena
            z[i] = d1[i]
        end
    end
    if lenb > lenmin
        for i in lenmin+1:lenb
            z[i] = d2[i]
        end
    end
    return R(z)
end

function -(a::PolyRingElem{T,S}, b::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    return a + (-b)
end

#Might be subject to more effective implementations later...
function *(a::PolyRingElem{T,S}, b::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    if !sameRing(a, b)
        error("Cannot multiply elements from different rings")
    end
    R = a.ring
    d1 = deepcopy(a.d)
    d2 = deepcopy(b.d)
    lena = length(d1)
    lenb = length(d2)
    if lena == 0 || lenb == 0
        return R([zero(R.parent)])
    end
    lenz = lena + lenb - 1
    z = zeros(R.parent, lenz)
    k = lenz-lena
    for i in 1:k
        push!(d1,zero(R.parent))
    end
    k = lenz-lenb
    for i in 1:k
        push!(d2,zero(R.parent))
    end
    for i in 1:lenz
        s = zero(R.parent)
        for j in 1:i
            s += d1[j] * d2[i+1-j]
        end
        z[i] = s
    end
    return R(z)
end

function *(a::PolyRingElem{T,S}, b::T) where {T <: CRingElem, S <: CRing}
    R = a.ring
    z = deepcopy(a.d)
    z = [k*b for k in z]
    return R(z)
end

function *(b::PolyRingElem{T,S}, a::T) where {T <: CRingElem, S <: CRing}
    R = b.ring
    z = deepcopy(b.d)
    z = [k*a for k in z]
    return R(z)
end

function *(n::U, a::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing, U <:Integer}
    if n < 0
        a = -a
        n = -n
    end
    if n == 1
        return a
    end
    Rpoly = a.ring
    R = Rpoly.parent
    if n == 0
        return zero(Rpoly)
    end
    bit = U(1) << (ndigits(n, base = 2) - 1)
    z = a
    bit >>= 1
    while bit != 0
        z = z+z
        if (bit & n) != 0
            z += a
        end
        bit >>=1
    end
    return z
end

function *(a::PolyRingElem{T,S}, n::U) where {T <: CRingElem, S <: CRing, U <:Integer}
    return n*a
end

#Again, just using square and multiply algorithm
function ^(a::PolyRingElem{T,S}, y::U) where {T <: CRingElem, S <: CRing, U <: Integer}
    if y < 0
        #Will implement invert function latert
        error("Not yet supported")
    end
    if y == 1
        return x
    end
    Rpoly = a.ring
    R = Rpoly.parent
    if y == 0
        return one(Rpoly)
    end
    # Should include some reduction of y modulo order if order is known
    bit = U(1) << (ndigits(y, base = 2) - 1)
    z = a
    bit >>= 1
    while bit != 0
        z = z*z
        if (bit & y) != 0
            z *= a
        end
        bit >>=1
    end
    return z
end

function ==(a::PolyRingElem{T,S}, b::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    return sameRing(a, b) && (a.d == b.d)
end

####### GCD (Very much needed for constructing quotient rings) #######
function gcd(a::PolyRingElem{T,S}, b::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    if !sameRing(a, b)
        error("Elements must be elements of same ring")
    end
    if degree(a) < degree(b)
        return gcd(b,a)
    end
    if iszero(a) && iszero(b)
        return zero(a.ring)
    end
    while degree(b) > -1
        q,r = divmod(a,b)
        a = b
        b = r
    end
    return monic(a)
end

function gcdx(a::PolyRingElem{T,S}, b::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    if !sameRing(a, b)
        error("Elements must be elements of same ring")
    end
    if degree(a) < degree(b)
        (a, y1, x1) = gcdx(b,a)
        return (a, x1, y1)
    end
    Rpol = a.ring
    if iszero(b)
        factor = inv(LC(a))
        return (a*factor, one(Rpol)*factor, zero(Rpol))
    end
    x0, x1 = zero(Rpol), one(Rpol)
    y1, y0 = zero(Rpol), one(Rpol)
    while degree(b) > -1
        q, r = divmod(a,b)
        x = x1 - q*x0
        y = y1 - q*y0
        a, b = b, r
        x0, x1, y0, y1 = x, x0, y, y0
    end
    factor = inv(LC(a))
    return (a*factor, x1*factor, y1*factor)
end

#gcdx subroutine
function divmod(a::PolyRingElem{T,S}, b::PolyRingElem{T,S}) where {T <: CRingElem, S <: CRing}
    if !sameRing(a, b)
        error("Elements must be elements of same ring")
    end
    Rpol = a.ring
    R = Rpol.parent
    lcdiv = LC(b)
    quotient = zero(Rpol)
    remainder = a
    while degree(remainder) >= degree(b)
        z = zeros(R, degree(remainder)-degree(b))
        push!(z,LC(remainder)*(inv(lcdiv)))
        z = Rpol(z)
        quotient += z
        remainder -= z*b
    end
    return quotient, remainder
end

##########################################
#
#         Quotient Ring
#
##########################################

struct QuotientRing{T <: CRing, S <: CRingElem} <: CRing
    parent::T
    modulo::S
end

struct QuotientRingElem{T <: CRing, S <: CRingElem} <: CRingElem
    d::S
    ring::QuotientRing{T,S}
end

# Ring-Constructor
function QuotientRing(R::T, m::S) where {T <: CRing, S <: CRingElem}
    return QuotientRing{T,S}(R,m)
end

# Element-Constructor
function (R::QuotientRing{T,S})(a::S) where {T <: CRing, S <: CRingElem}
    m = R.modulo
    q,r = divmod(a,m)
    return QuotientRingElem{T,S}(r, R)
end

function (R::QuotientRing{PolyRing{T},PolyRingElem{S,T}})(a::Array{S}) where {T <: CRing, S <: CRingElem}
    a = R.parent(a)
    m = R.modulo
    q,r = divmod(a,m)
    return QuotientRingElem{PolyRing{T},PolyRingElem{S,T}}(r, R)
end

###### Printing ######
function show(io::IO, R::QuotientRing{T,S}) where {T <: CRing, S <: CRingElem}
    print(io, "$(R.parent)/<$(R.modulo)>")
end

function show(io::IO, a::QuotientRingElem{T,S}) where {T <: CRing, S <: CRingElem}
    print(io, "$(a.d)")
end

###### Helper functions #######

function zero(R::QuotientRing{T,S}) where {T <: CRing, S <: CRingElem}
    return R(zero(R.parent))
end

function one(R::QuotientRing{T,S}) where {T <: CRing, S <: CRingElem}
    return R(one(R.parent))
end

function iszero(a::QuotientRingElem{T,S}) where {T <: CRing, S <: CRingElem}
    return iszero(a.d)
end

function isone(a::QuotientRingElem{T,S}) where {T <: CRing, S <: CRingElem}
    return isone(a.d)
end

function order(R::QuotientRing{PolyRing{T},PolyRingElem{S,T}}) where {T <: CRing, S <: CRingElem}
    return order(R.parent.parent)^degree(R.modulo)
end

###### Basic operators #######
function -(a::QuotientRingElem{T,S}) where {T <: CRing, S <: CRingElem}
    if iszero(a)
        return a
    else
        R = a.ring
        return R(-a.d)
    end
end

function +(a::QuotientRingElem{T,S}, b::QuotientRingElem{T,S}) where {T <: CRing, S <: CRingElem}
    if !sameRing(a, b)
        error("Cannot add elements of different rings")
    end
    z = a.d + b.d
    R = a.ring
    return R(z)
end

function -(a::QuotientRingElem{T,S}, b::QuotientRingElem{T,S}) where {T <: CRing, S <: CRingElem}
    return a+(-b)
end

function *(a::QuotientRingElem{T,S}, b::QuotientRingElem{T,S}) where {T <: CRing, S <: CRingElem}
    if !sameRing(a, b)
        error("Cannot multiply elements of different rings")
    end
    z = a.d * b.d
    R = a.ring
    return R(z)
end

function *(n::Integer, a::QuotientRingElem{T,S}) where {T <: CRing, S <: CRingElem}
    z = n*a.d
    R = a.ring
    return R(z)
end

function *(a::QuotientRingElem{T,S}, n::Integer) where {T <: CRing, S <: CRingElem}
    return n*a
end

function inv(x::QuotientRingElem{T,S}) where {T <: CRing, S <: CRingElem}
    R = x.ring
    g, a, b = gcdx(x.d, R.modulo)
    (!isone(g)) && error("Element $x is not invertible in $(x.ring)")
    return R(a)
end

#yet again exponensiation uses the simple square and multiply algorithm
function ^(x::QuotientRingElem{T,S}, y::U) where {T <: CRing, S <: CRingElem, U <: Integer}
    if iszero(x) || y == 1
        return x
    end
    R = x.ring
    if y == 0
        return one(R)
    end
    if y < 0
        x = inv(x)
        y = -y
    end
    # Should include some reduction of y modulo order if order is known
    bit = U(1) << (ndigits(y, base = 2) - 1)
    z = x
    bit >>= 1
    while bit != 0
        z = z*z
        if (bit & y) != 0
            z *= x
        end
        bit >>=1
    end
    return z
end

function ==(a::QuotientRingElem{T,S}, b::QuotientRingElem{T,S}) where {T <: CRing, S <: CRingElem}
    return sameRing(a, b) && (a.d == b.d)
end

##########################################
#
#      Finite field of order p^2
#
##########################################
# Just a special case of the construction above for easy usage
# If p = 3 (mod 4), then -1 is a quadratic non-residue and x^2 + 1 is irredusible
function Fpsquare(p::T) where T <: Integer
    Zp = Zmod(p)
    ZpPol = PolyRing(Zp, "i")
    if !isprime(p)
        error("$p is not a prime!")
    end
    if p%4 != 3
        error("Only supported for p ≡ 3 (mod 4)")
    end
    m = ZpPol([Zp(T(1)),Zp(T(0)),Zp(T(1))])
    return QuotientRing(ZpPol, m)
end

# Element-Constructor
function (R::QuotientRing{PolyRing{Zmod{T}},PolyRingElem{ZmodElem{T},Zmod{T}}})(a::Array{S}) where {T <: Integer, S <: Integer}
    a = [T(x) for x in a]
    a = R.parent([R.parent.parent(i) for i in a])
    m = R.modulo
    q,r = divmod(a,m)
    return QuotientRingElem{PolyRing{Zmod{T}},PolyRingElem{ZmodElem{T},Zmod{T}}}(r, R)
end

#################################
#
# Dangerous or specific functions
#
#################################

function squareroot(a::ZmodElem{T}) where {T <: Integer}
    return a.ring(TShanks(a.d, a.ring.n)[1])
end

# This is just for the special case, Zp[X]/<f> where f is like x^2+1...
function cuberootOfOne(R::QuotientRing{PolyRing{Zmod{T}},PolyRingElem{ZmodElem{T},Zmod{T}}}) where {T <: Integer, S <: Integer}
    #idea is that the non trivial cube roots of 1 are -1±sqrt(-3)/2 as x^3-1 = (x-1)(x^2+x+1)
    Z = R.parent.parent
    if isone(Z(-3)^((Z.n-1)÷2))
        d = squareroot(Z(-3))
        u = R([d])
    elseif isone(Z(3)^((Z.n-1)÷2))
        d = squareroot(Z(3))
        u = R([zero(Z), d])
    else
        error("Probably no cuberoot of 1 in the field")
    end
    id = one(R)
    z = inv(2*id)
    return (-id+u)*z
end

# Make sure that the ring is actually a field before using this... Further 3|(q-1)
function cuberoot(a::CRingElem)
    q = order(a.ring)
    if q%9 != 7
        error("Currently only supported for q ≡ 7 (mod 9)")
    end
    #general euler criterion
    k = (q-1)÷3
    if !isone(a^k)
        error("$a is not a cube in $(a.ring)")
    end
    k = (q+2)÷9
    return a^k
end

function allcuberoots(a::CRingElem)
    w = cuberootOfOne(a.ring)
    k = cuberoot(a)
    return (k, k*w, k*w^2)
end

function allElements(R::QuotientRing{PolyRing{Zmod{T}},PolyRingElem{ZmodElem{T},Zmod{T}}}) where T <: Integer
    deg = degree(R.modulo)
    if deg != 2
        println("Too many!")
    end
    Z = R.parent.parent
    l = QuotientRingElem{PolyRing{Zmod{T}},PolyRingElem{ZmodElem{T},Zmod{T}}}[]
    for i1 in allElements(Z)
        for i2 in allElements(Z)
            push!(l,R([i1,i2]))
        end
    end
    return l
end

#################################
#
#            Other
#
#################################

legendre(a, p) = powermod(a, (p - 1) ÷ 2, p)

function TShanks(n::T, p::T) where T <: Union{Int, Int128, BigInt}
    legendre(n, p) != 1 && throw(ArgumentError("$n not a square (mod $p)"))
    local q::T = p - one(p)
    local s::T = 0
    while iszero(q % 2)
        q ÷= 2
        s += one(s)
    end
    if s == one(s)
        r = powermod(n, (p + 1) >> 2, p)
        return r, p - r
    end
    local z::T
    for z in 2:(p - 1)
        p - 1 == legendre(z, p) && break
    end
    local c::T = powermod(z, q, p)
    local r::T = powermod(n, (q + 1) >> 1, p)
    local t::T = powermod(n, q, p)
    local m::T = s
    local t2::T = zero(p)
    while !iszero((t - 1) % p)
        t2 = (t * t) % p
        local i::T
        for i in Base.OneTo(m)
            iszero((t2 - 1) % p) && break
            t2 = (t2 * t2) % p
        end
        b = powermod(c, 1 << (m - i - 1), p)
        r = (r * b) % p
        c = (b * b) % p
        t = (t * c) % p
        m = i
    end
    return r, p - r
end

end


text:
using .DataFrames
import .DataFrames: DataFrame

"""

# convert_a3d

# Convert the output file(s) created by cmdstan to a DataFrame.

$(SIGNATURES)

"""
function convert_a3d(a3d_array, cnames, ::Val{:dataframe})
  # Inital DataFrame
  df = DataFrame(a3d_array[:, :, 1], Symbol.(cnames))

  # Append the other chains
  for j in 2:size(a3d_array, 3)
    df = vcat(df, DataFrame(a3d_array[:, :, j], Symbol.(cnames)))
  end
  df
end

"""

# convert_a3d

# Convert the output file(s) created by cmdstan to a Vector{DataFrame).

$(SIGNATURES)

"""
function convert_a3d(a3d_array, cnames, ::Val{:dataframes})

  dfa = Vector{DataFrame}(undef, size(a3d_array, 3))
  for j in 1:size(a3d_array, 3)
    dfa[j] = DataFrame(a3d_array[:, :, j], Symbol.(cnames))
  end

  dfa
end

"""

DataFrame()

# Block Stan named parameters, e.g. b.1, b.2, ... in a DataFrame.

$(SIGNATURES)

Example:

df_log_lik = DataFrame(m601s_df, :log_lik)
log_lik = Matrix(df_log_lik)

"""
function DataFrame(df::DataFrame, sym::Union{Symbol, String})
    n = string.(names(df))
    syms = string(sym)
    sel = String[]
    for (i, s) in enumerate(n)
        if length(s) > length(syms) && syms == n[i][1:length(syms)] &&
            n[i][length(syms)+1] in ['[', '.', '_']
            append!(sel, [n[i]])
        end
    end
    length(sel) == 0 && error("$syms not in $n")
    df[:, sel]
end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 1, Zadanie 1.1.
#

# assert the program’s arguments
if length(ARGS) == 0
    println("error: provide the floating point precision (16, 32 or 64 bits)")
    exit(1)
end

fl_precision_string = ARGS[1]

if fl_precision_string == "16"
    fl = Float16
elseif fl_precision_string == "32"
    fl = Float32
elseif fl_precision_string == "64"
    fl = Float64
else
    println("error: provide one of [16, 32, 64]")
    exit(1)
end
# assert end

# start finding the machine epsilon at 0.5
mach_eps = fl(1/2)

# iteratively find the machine epsilon
for i in range(1, length=100)
    # if adding even smaller eps would not change the value break the loop
    if fl(1) + fl(mach_eps/2) == fl(1)
        break
    end
    # otherwise continue
    global mach_eps = fl(mach_eps/2)
end

# print the output
println("iterative method = ", mach_eps, "; eps(fl) = ", eps(fl), "; equal = ", (eps(fl) == mach_eps))

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
# # Linear regression

using DynamicHMCModels

ProjDir = @__DIR__
cd(ProjDir)

# Import the dataset.

# ### snippet 5.4

df = DataFrame(CSV.read(joinpath("..", "..", "data", "WaffleDivorce.csv"), delim=';'))

mean_ma = mean(df[:, :Marriage])
df[!, :Marriage_s] = convert(Vector{Float64},
  (df[:, :Marriage]) .- mean_ma)/std(df[:, :Marriage]);

mean_mam = mean(df[:, :MedianAgeMarriage])
df[!, :MedianAgeMarriage_s] = convert(Vector{Float64},
  (df[:, :MedianAgeMarriage]) .- mean_mam)/std(df[:, :MedianAgeMarriage]);

# Model ``y ∼ Xβ + ϵ``, where ``ϵ ∼ N(0, σ²)`` IID. Student on σ

Base.@kwdef mutable struct WaffleDivorce{Ty <: AbstractVector, Tx <: AbstractMatrix}
    "Observations."
    y::Ty
    "Covariates"
    x::Tx
end

# Write a function to return a properly dimensioned transformation.

function make_transformation(model::WaffleDivorce)
  as((β = as(Array, size(model.x, 2)), σ = asℝ₊))
end

# Instantiate the model with data and inits.

x = hcat(ones(size(df, 1)), df[:, :Marriage_s], df[:, :MedianAgeMarriage_s]);
model = WaffleDivorce(;y=df[:, :Divorce], x=x);

# Make the type callable with the parameters *as a single argument*.

function (model::WaffleDivorce)(θ)
    @unpack y, x = model   # extract the data
    @unpack β, σ = θ            # works on the named tuple too
    ll = 0.0
    ll += logpdf(Normal(10, 10), x[1])
    ll += logpdf(Normal(0, 1), x[2])
    ll += logpdf(Normal(0, 1), x[3])
    ll += logpdf(TDist(1.0), σ)
    ll += loglikelihood(Normal(0, σ), y .- x*β)
    ll
end

println()
model((β = [1.0, 2.0, 3.0], σ = 1.0)) |> display
println()

# Instantiate the model with data and inits.

println()
model((β = [1.0, 2.0, 3.0], σ = 1.0))
println()

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P);

# Tune and sample.

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

println()
DynamicHMC.Diagnostics.EBFMI(results.tree_statistics) |> display
println()

DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics) |> display
println()

a3d = Array{Float64, 3}(undef, 1000, 4, 1);
for j in 1:1
  for i in 1:1000
    a3d[i, 1:3, j] = values(posterior[i].β)
    a3d[i, 4, j] = values(posterior[i].σ)
  end
end

pnames = ["β[1]", "β[2]", "β[3]", "σ"]
sections =   Dict(
  :parameters => pnames,
)
chns = create_mcmcchains(a3d, pnames, sections, start=1);
chns = set_names(chns, Dict("β[1]" => "α", "β[2]" => "β[1]", "β[3]" => "β[2]"))

stan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
          Mean        SD       Naive SE       MCSE      ESS
    a  9.69137275 0.21507432 0.0034006235 0.0038501180 1000
   bA -1.12184710 0.29039965 0.0045916216 0.0053055477 1000
   bM -0.12106472 0.28705400 0.0045387223 0.0051444688 1000
sigma  1.52326545 0.16272599 0.0025729239 0.0034436330 1000

Quantiles:
         2.5%       25.0%      50.0%      75.0%       97.5%
    a  9.2694878  9.5497650  9.6906850  9.83227750 10.11643500
   bA -1.6852295 -1.3167700 -1.1254650 -0.92889225 -0.53389157
   bM -0.6889247 -0.3151695 -0.1231065  0.07218513  0.45527243
sigma  1.2421182  1.4125950  1.5107700  1.61579000  1.89891925
";

describe(chns)

# end of m4.5d.jl


text:
using CodeToString
using Base.Test

try
    f(a, b) = a + b
    code_llvm_to_bytestring(f, (Int64, Int64))
    @code_llvm_to_bytestring f(1,2)
    @test typeof(code_llvm_to_bytestring(f, (Int64, Int64))) == Array{UInt8,1}
    @test typeof(@code_llvm_to_bytestring f(1,2)) == Array{UInt8,1}
    code_native_to_bytestring(f, (Int64, Int64))
    @code_native_to_bytestring f(1,2)
    @test typeof(code_native_to_bytestring(f, (Int64, Int64))) == Array{UInt8,1}
    @test typeof(@code_native_to_bytestring f(1,2)) == Array{UInt8,1}
    code_warntype_to_bytestring(f, (Int64, Int64))
    @code_warntype_to_bytestring f(1,2)
    @test typeof(code_warntype_to_bytestring(f, (Int64, Int64))) == Array{UInt8,1}
    @test typeof(@code_warntype_to_bytestring f(1,2)) == Array{UInt8,1}
catch
    @test false
end

@test true


text:
include("address_manipulate.jl")
include("optimise.jl")

print_blocktree() = print_subtypetree(AbstractBlock)


text:
module SingleDispatchArrays
export SingleDispatchVector

"Stores all subtypes for a given base type as Basetype => Type{Tuple{Subtypes...}}."
const subtype_dict = Dict()

get_subtypes(basetype) = get(subtype_dict, basetype, Tuple{})
function add_subtype(basetype, new_type)
    old_types = fieldtypes(get_subtypes(basetype))
    if !(new_type in old_types)
        subtype_dict[basetype] = Tuple{old_types..., new_type}
    end
end

"""
    SingleDispatchVector{BaseType} <: AbstractVector
Vector of non-homogenous elements with fast static dispatch table.
!!! note
A full array interface has not been implemented, please use [`index_value_pair`](@ref) to manually manipulate the underlying vector.
"""
mutable struct SingleDispatchVector{T} <: AbstractVector{T}
    vec::Vector{Pair{UInt, T}}
    # Cache of `get_subtypes(basetype)`, only updated when inserting elements of uncached type.
    subtypes::Type
end

SingleDispatchVector{T}() where T = SingleDispatchVector{T}([], get_subtypes(T))

function SingleDispatchVector{T}(vec::AbstractVector{T}) where T
    subtypes = get_subtypes(T)
    pairs = map(vec) do v
        Pair(type_index(subtypes, v), v)
    end
    SingleDispatchVector{T}(pairs, subtypes)
end

basetype(::SingleDispatchVector{T}) where T = T

function generate_dispatch_table(types)
    root_expr = :(if false nothing end)
    cur = root_expr
    for (i, t) in enumerate(types)
        branch = Expr(:elseif, :(ind == $i), :(f(value::$t)))
        push!(cur.args, branch)
        cur = branch
    end
    fallback = :(error("Element $value has invalid type index $ind for basetype $(basetype(a))"))
    push!(cur.args, fallback)
    return root_expr
end

@generated function Base.foreach(f, a::SingleDispatchVector, subtypes::Type{T}) where T
    types = fieldtypes(T)
    dispatch_table = generate_dispatch_table(types)
    return quote
        for (ind, value) in a.vec
            $dispatch_table
        end
    end
end

# Helper function to pass the signature
Base.foreach(f, a::SingleDispatchVector) = foreach(f, a, a.subtypes)

# Too complicated for compiler to fully constant propagation, so the index is directly generated.
@generated function type_index(::Type{Ts}, ::T) where {Ts,T}
    ind = 0
    for (i, t) in enumerate(fieldtypes(Ts))
        if t == T
            ind = i
        end
    end

    return ind
end

function index_value_pair(a::SingleDispatchVector, value)
    # Fast path
    ind = type_index(a.subtypes, value)
    iszero(ind) || return Pair(ind, value)
    # Update subtype set and try again
    a.subtypes = get_subtypes(basetype(a))
    ind = type_index(a.subtypes, value)
    if iszero(ind)
        error("Type $(typeof(value)) is not a recorded subtype of $(basetype(a))")
    else
        Pair(ind, value)
    end
end

# Array interface (very incomplete)
Base.setindex!(a::SingleDispatchVector, val, key) = Base.setindex!(a.vec, index_value_pair(a, val), key)
Base.push!(a::SingleDispatchVector, val) = Base.push!(a.vec, index_value_pair(a, val))
Base.getindex(a::SingleDispatchVector, i) = Base.getindex(a.vec, i)|>last
Base.size(a::SingleDispatchVector) = Base.size(a.vec)
Base.resize!(a::SingleDispatchVector, n) = Base.resize!(a.vec, n)

end # module


text:
struct MbplsrAvg
    fm
    w::Vector{Float64}
end

function mbplsr_avg(X, Y, weights = ones(size(X[1], 1)); nlv, w = nothing)
    nbl = length(X)
    isnothing(w) ? w = ones(nbl) / nbl : nothing
    fm = list(nbl)
    @inbounds for i = 1:nbl
        fm[i] = plskern(X[i], Y, weights; nlv = nlv[i])
    end
    MbplsrAvg(fm, w)
end

"""
    predict(object::MbplsrAvg, X)
Compute Y-predictions from a fitted model.
* `object` : The fitted model.
* `X` : A list (vector) of X-data for which predictions are computed.
"""
function predict(object::MbplsrAvg, X)
    nbl = length(object.fm)
    pred = object.w[1] * predict(object.fm[1], X[1]).pred
    if nbl > 1
        @inbounds for i = 2:nbl
            pred .+= object.w[i] * predict(object.fm[i], X[i]).pred
        end
    end
    (pred = pred,)
end


text:
## Adapted from  doc/libsundials-serial-dev/examples/ida/serial/idaRoberts_dns.c and
##               sundialsTB/ida/examples_ser/midasRoberts_dns.m

## /*
##  * -----------------------------------------------------------------
##  * $Revision: 1.2 $
##  * $Date: 2009/01/21 21:46:40 $
##  * -----------------------------------------------------------------
##  * Programmer(s): Allan Taylor, Alan Hindmarsh and
##  *                Radu Serban @ LLNL
##  * -----------------------------------------------------------------
##  * This simple example problem for IDA, due to Robertson,
##  * is from chemical kinetics, and consists of the following three
##  * equations:
##  *
##  *      dy1/dt = -.04*y1 + 1.e4*y2*y3
##  *      dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
##  *         0   = y1 + y2 + y3 - 1
##  *
##  * on the interval from t = 0.0 to t = 4.e10, with initial
##  * conditions: y1 = 1, y2 = y3 = 0.
##  *
##  * While integrating the system, we also use the rootfinding
##  * feature to find the points at which y1 = 1e-4 or at which
##  * y3 = 0.01.
##  *
##  * The problem is solved with IDA using IDADENSE for the linear
##  * solver, with a user-supplied Jacobian. Output is printed at
##  * t = .4, 4, 40, ..., 4e10.
##  * -----------------------------------------------------------------
##  */

using Sundials

## Define the system residual function.
function resrob(tres, yy_nv, yp_nv, rr_nv, user_data)
    yy = convert(Vector, yy_nv)
    yp = convert(Vector, yp_nv)
    rr = convert(Vector, rr_nv)
    rr[1]  = -0.04*yy[1] + 1.0e4*yy[2]*yy[3]
    rr[2]  = -rr[1] - 3.0e7*yy[2]*yy[2] - yp[2]
    rr[1] -=  yp[1]
    rr[3]  =  yy[1] + yy[2] + yy[3] - 1.0
    return Sundials.IDA_SUCCESS
end

## Root function routine. Compute functions g_i(t,y) for i = 0,1.
function grob(t, yy_nv, yp_nv, gout_ptr, user_data)
    yy = convert(Vector, yy_nv)
    gout = Sundials.asarray(gout_ptr, (2,))
    gout[1] = yy[1] - 0.0001
    gout[2] = yy[3] - 0.01
    return Sundials.IDA_SUCCESS
end

## Define the Jacobian function. BROKEN - JJ is wrong
function jacrob(Neq, tt, cj, yy, yp, resvec,
                JJ, user_data, tempv1, tempv2, tempv3)
    JJ = pointer_to_array(convert(Ptr{Float64}, JJ), (3,3))
    JJ[1,1] = -0.04 - cj
    JJ[2,1] = 0.04
    JJ[3,1] = 1.0
    JJ[1,2] = 1.0e4*yy[3]
    JJ[2,2] = -1.0e4*yy[3] - 6.0e7*yy[2] - cj
    JJ[3,2] = 1.0
    JJ[1,3] = 1.0e4*yy[2]
    JJ[2,3] = -1.0e4*yy[2]
    JJ[3,3] = 1.0
    return Sundials.IDA_SUCCESS
end

neq = 3
nout = 12
t0 = 0.0
yy0 = [1.0,0.0,0.0]
yp0 = [-0.04,0.04,0.0]
rtol = 1e-4
avtol = [1e-8, 1e-14, 1e-6]
tout1 = 0.4

mem = Sundials.IDACreate()
Sundials.@checkflag Sundials.IDAInit(mem, resrob, t0, yy0, yp0)
Sundials.@checkflag Sundials.IDASVtolerances(mem, rtol, avtol)

## Call IDARootInit to specify the root function grob with 2 components
Sundials.@checkflag Sundials.IDARootInit(mem, 2, grob)

## Call IDADense and set up the linear solver.
A = Sundials.SUNDenseMatrix(length(y0),length(y0))
LS = Sundials.SUNDenseLinearSolver(y0,A)
Sundials.@checkflag Sundials.IDADlsSetLinearSolver(mem, LS, A)

iout = 0
tout = tout1
tret = [1.0]

while iout < nout
    yy = similar(yy0)
    yp = similar(yp0)
    retval = Sundials.IDASolve(mem, tout, tret, yy, yp, Sundials.IDA_NORMAL)
    println("T=", tout, ", Y=", yy)
    if retval == Sundials.IDA_ROOT_RETURN
        rootsfound = zeros(Cint, 2)
        Sundials.@checkflag Sundials.IDAGetRootInfo(mem, rootsfound)
        println("roots=", rootsfound)
    elseif retval == Sundials.IDA_SUCCESS
        global iout += 1
        global tout *= 10.0
    end
end

Sundials.SUNLinSolFree_Dense(LS)
Sundials.SUNMatDestroy_Dense(A)


text:
using JSON

macro test_str(s)
    return quote
        str = @interpolate $s
        # Do what you want to do with interpolated string here.
        sprint(print,str^3)
    end
end

function test1(x)
    test"Non-standard string literal has interpolation: $x\n"
end

println("************")
println("Test Result:")
println("************")
println("Call test1(2):")
println(test1(2))
println("Call test1(Dict(:a=>2)):")
println(test1(Dict(:a=>2)))