text:
using SFGAnalysis
using Test
import Distributions: Normal

@testset "SFGAnalysis.jl" begin
    λ1 = 515
    λ2 = 3450
    ω1 = wl2freq(λ1)
    ω2 = wl2freq(λ2)
    ω = ω1 + ω2
    λ = freq2wl(ω)
    β1 = π/4
    β2 = π/6
    β = sfangle(ω, ω1, ω2, β1, β2)

    #TODO: are these really frequencies or angular frequencies or wavenumbers?
    @test wl2freq(λ1) ≈ 5.821212776699029e14 # source?
    @test λ ≈ 448.1084489281211  # source? but looks correct

    # Test Setup Type
    n1  = 1.0
    n11 = 1.0
    n12 = 1.0
    n2  = 1.3
    n21 = 1.3
    n22 = 1.4
    s = Setup((ω1, ω2), (β1, β2), (n1, n11, n12), (n2, n21, n22))

    # test sfgintensity
    χ2_eff = 1e-5 + 1e-6im   # completely arbitrary. change for calculated value
    @test sfgintensity(ω, β, χ2_eff) == sfgintensity(s, χ2_eff)

    # Test ninterface
    # https://doi.org/10.1103/PhysRevB.59.1263 for verification
    @test round(ninterface(1.5), digits=2) == 1.22
    @test round(ninterface(1.34), digits=2) == 1.15

    # Test susceptibility
    for t in [:yyz, :yzy, :zyy, :xxz, :xzx, :zxx, :zzz]
        susceptibility(π/4, t; pointgroup=:c3v, mode=:ss, R=1.7)
    end
    for t in [:yyz, :yzy, :zyy, :xxz, :xzx, :zxx, :zzz]
        susceptibility(π/4, t; pointgroup=:c3v, mode=:as)
    end
    @test_throws ErrorException susceptibility(π/4, :yyz)
    @test_throws ErrorException susceptibility(π/4, :yyz; pointgroup=:c3v)
    @test_throws ErrorException susceptibility(π/4, :yyz; mode=:c3v)
    @test_throws ErrorException susceptibility(π/4, :yyz; pointgroup=:c3v, mode=:c3v)

    # Test effective_susceptibility
    χ_ssp = susceptibility(π/8, :yyz; pointgroup=:c3v, mode=:as)
    χ_ppp = [susceptibility(π/8, t; pointgroup=:c3v, mode=:as) for t in [:xxz, :xzx, :zxx, :zzz]]
    χ2_eff_ssp = effective_susceptibility_ssp(s, χ_ssp)
    χ2_eff_ppp = effective_susceptibility_ppp(s, χ_ppp)

end

@testset "J. Phys. Chem. C 2007, 111, 8716-8725" begin
    # Test https://doi.org/10.1021/jp067062h
    ω1 = wl2freq(532)
    ω2 = wl2freq(3333)
    β1 = [deg2rad(62), deg2rad(37)]
    β2 = [deg2rad(53), deg2rad(51)]
    n1  = 1
    n11 = 1
    n12 = 1
    n2  = 1.33
    n21 = 1.33
    n22 = 1.33
    R = 1.7
    s = Setup((ω1, ω2), (β1[1], β2[1]), (n1, n11, n12), (n2, n21, n22))

    θ = range(0, π/2, length=91)

    # a factor is appended to the effective susceptibilities
    χ_eff_ssp_ss = effective_susceptibility.(θ, s, :ssp; pointgroup=:c3v, mode=:ss, R=R) * sec(s.Ω.β)
    χ_eff_ssp_as = effective_susceptibility.(θ, s, :ssp; pointgroup=:c3v, mode=:as) * sec(s.Ω.β)
    χ_eff_ppp_ss = effective_susceptibility.(θ, s, :ppp; pointgroup=:c3v, mode=:ss, R=R) * sec(s.Ω.β)
    χ_eff_ppp_as = effective_susceptibility.(θ, s, :ppp; pointgroup=:c3v, mode=:as) * sec(s.Ω.β)
    χ_eff_sps_ss = effective_susceptibility.(θ, s, :sps; pointgroup=:c3v, mode=:ss, R=R) * sec(s.Ω.β)
    χ_eff_sps_as = effective_susceptibility.(θ, s, :sps; pointgroup=:c3v, mode=:as) * sec(s.Ω.β)
    χ_eff_pss_ss = effective_susceptibility.(θ, s, :pss; pointgroup=:c3v, mode=:ss, R=R) * sec(s.Ω.β)
    χ_eff_pss_as = effective_susceptibility.(θ, s, :pss; pointgroup=:c3v, mode=:as) * sec(s.Ω.β)

    @test round(χ_eff_ssp_ss[1], digits=1) == 0.9
    @test round(χ_eff_ssp_ss[45], digits=1) == 0.6
    @test round(χ_eff_ssp_ss[90], digits=1) == 0.0

    @test round(χ_eff_ppp_ss[1], digits=1) == -0.1
    @test round(χ_eff_ppp_ss[45], digits=1) == 0.1
    @test round(χ_eff_ppp_ss[90], digits=1) == 0.0

    # The susceptibilities for pss polarization are the same as for
    # sps polarization. Differences in effective susceptibility are
    # solely due to differences in the fresnel factors with pss
    # polarization generally having a smaller eff. sus.
    @test abs(χ_eff_pss_ss[1]) <= abs(χ_eff_sps_ss[1])
    @test abs(χ_eff_pss_ss[45]) <= abs(χ_eff_sps_ss[45])
    @test abs(χ_eff_pss_ss[90]) <= abs(χ_eff_sps_ss[90])
    @test abs(χ_eff_pss_as[1]) <= abs(χ_eff_sps_as[1])
    @test abs(χ_eff_pss_as[45]) <= abs(χ_eff_sps_as[45])
    @test abs(χ_eff_pss_as[90]) <= abs(χ_eff_sps_as[90])
    # check sign
    @test sign(χ_eff_pss_ss[1]) == sign(χ_eff_sps_ss[1])
    @test sign(χ_eff_pss_ss[45]) == sign(χ_eff_sps_ss[45])
    @test sign(χ_eff_pss_ss[90]) == sign(χ_eff_sps_ss[90])
    @test sign(χ_eff_pss_as[1]) == sign(χ_eff_sps_as[1])
    @test sign(χ_eff_pss_as[45]) == sign(χ_eff_sps_as[45])
    @test sign(χ_eff_pss_as[90]) == sign(χ_eff_sps_as[90])

    @test round(χ_eff_sps_ss[1], digits=1) == 0.0
    @test round(χ_eff_sps_ss[45], digits=1) == -0.1
    @test round(χ_eff_sps_ss[90], digits=1) == 0.0

    @test round(χ_eff_ssp_as[1], digits=1) == 0.0
    @test round(χ_eff_ssp_as[45], digits=1) == -0.2
    @test round(χ_eff_ssp_as[90], digits=1) == 0.0

    @test round(χ_eff_ppp_as[1], digits=1) == 0.0
    @test round(χ_eff_ppp_as[45], digits=1) == 0.4
    @test round(χ_eff_ppp_as[90], digits=1) == 0.0

    @test round(χ_eff_sps_as[1], digits=1) == 0.6
    @test round(χ_eff_sps_as[45], digits=1) == 0.2
    @test round(χ_eff_sps_as[90], digits=1) == 0.0
end

#TODO: Add tests for pss

@testset "Distributions" begin
    # Test the distribution functionality
    # Do everything like in the testse "J. Phys. Chem. C 2007, 111, 8716-8725"
    # This time we just assume a very narrow distribution. This way the results
    # should be the same.
    ω1 = wl2freq(532)
    ω2 = wl2freq(3333)
    β1 = [deg2rad(62), deg2rad(37)]
    β2 = [deg2rad(53), deg2rad(51)]
    n1  = 1
    n11 = 1
    n12 = 1
    n2  = 1.33
    n21 = 1.33
    n22 = 1.33
    R = 1.7
    s = Setup((ω1, ω2), (β1[1], β2[1]), (n1, n11, n12), (n2, n21, n22))
    display(s)

    dist = [Normal(θ, π/100) for θ in range(0, π/2, length=91)]

    # a factor is appended to the effective susceptibilities
    χ_eff_ssp_ss = effective_susceptibility.(dist, s, :ssp; pointgroup=:c3v, mode=:ss, R=R) * sec(s.Ω.β)
    χ_eff_ssp_as = effective_susceptibility.(dist, s, :ssp; pointgroup=:c3v, mode=:as) * sec(s.Ω.β)
    χ_eff_ppp_ss = effective_susceptibility.(dist, s, :ppp; pointgroup=:c3v, mode=:ss, R=R) * sec(s.Ω.β)
    χ_eff_ppp_as = effective_susceptibility.(dist, s, :ppp; pointgroup=:c3v, mode=:as) * sec(s.Ω.β)

    @test round(χ_eff_ssp_ss[1], digits=1) == 0.9
    @test round(χ_eff_ssp_ss[45], digits=1) == 0.6
    @test round(χ_eff_ssp_ss[90], digits=1) == 0.0

    @test round(χ_eff_ppp_ss[1], digits=1) == -0.1
    @test round(χ_eff_ppp_ss[45], digits=1) == 0.1
    @test round(χ_eff_ppp_ss[90], digits=1) == 0.0

    @test round(χ_eff_ssp_as[1], digits=1) == 0.0
    @test round(χ_eff_ssp_as[45], digits=1) == -0.2
    @test round(χ_eff_ssp_as[90], digits=1) == 0.0

    @test round(χ_eff_ppp_as[1], digits=1) == 0.0
    @test round(χ_eff_ppp_as[45], digits=1) == 0.4
    @test round(χ_eff_ppp_as[90], digits=1) == 0.0
end

@testset "Spectra" begin
    x = range(2850, 3050, length=201)
    A = [1, 1]
    ω = [2900, 3000]
    Γ = [8, 8]
    φ_number = 0
    φ_array  = [0, 0]
    χnr_real = 0
    χnr_comp = Complex(0)
    y1 = sfspectrum.(x, Ref(A), Ref(ω), Ref(Γ))
    # φ belongs to χnr
    y2 = sfspectrum.(x, Ref(A), Ref(ω), Ref(Γ), φ_number, χnr_real)
    # φ belongs to A
    y3 = sfspectrum.(x, Ref(A), Ref(ω), Ref(Γ), Ref(φ_array),  χnr_real)
    # no phase
    y4 = sfspectrum.(x, Ref(A), Ref(ω), Ref(Γ), χnr_real)
    # all spectra should be equal
    @test y1 == y2 == y3 == y4
end


text:
function SeisSmoothStructure(m::AbstractString,d::AbstractString,param)

	adj = get(param,"adj",false)
	param["f"] = [smooth_structure]
	param["group"] = "gather"
	param["key"] = ["iaz","iang"]
	tmp1 = join(["tmp_SmoothStructure1_",string(int(rand()*100000))])
	tmp2 = join(["tmp_SmoothStructure2_",string(int(rand()*100000))])
	if (adj==false)
		SeisSort(m,tmp1,key=["iaz","iang","imy","imx"]);
		SeisProcess(tmp1,tmp2,param)
		SeisSort(tmp2,d,key=["imx","imy","iang","iaz"]);
	else
		SeisSort(d,tmp1,key=["iaz","iang","imy","imx"]);
		SeisProcess(tmp1,tmp2,param)
		SeisSort(tmp2,m,key=["imx","imy","iang","iaz"]);
	end
	SeisRemove(tmp1)
	SeisRemove(tmp2)

end

function smooth_structure(d,h,param)

	Nsmooth2 = get(param,"Nsmooth2",3)
	dip_field = get(param,"dip","dip")
	dip,h_dip = SeisRead(dip_field)
	#dip2 = atan(abs(dip))
	dip2 = atan(dip)

	nx = size(d,2)
	nz = size(d,1)
	WL = int(Nsmooth2/2)
	d2 = copy(d)
	smootharray!(d,d2,dip,dip2,nx,nz,Nsmooth2,WL)

	return d2,h;
end

function smootharray!(d::Array{Float32,2},d2::Array{Float32,2},dip::Array{Float32,2},dip2::Array{Float32,2},nx::Int64,nz::Int64,Nsmooth2::Int64,WL::Int64)

	for ix = WL : size(d,2) - WL
		for iz = WL : size(d,1) - WL
			@inbounds d2[iz,ix] = 0.
			@inbounds d2[iz,ix] = smooth1pixel(d,dip[iz,ix],cos(dip2[iz,ix]),sin(dip2[iz,ix]),ix - int(WL*cos(abs(dip2[iz,ix]))),iz - int(WL*sin(dip2[iz,ix])),nx,nz,Nsmooth2)
		end
	end

end

function smooth1pixel(m::Array{Float32,2},dip::Float32,xval::Float32,zval::Float32,ixo::Int64,izo::Int64,nx::Int64,nz::Int64,nw::Int64)

	val = 0.
	for iw = 1 : nw
		@inbounds val += m[izo + int(iw*zval),ixo + int(iw*xval)]
	end

	return val/nw

end

ext: jl

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
using Test, ModifiedHankelFunctionsOfOrderOneThird

println("\nStarting tests")

@time include("ModifiedHankelFunctionTests.jl")


text:
@assert NA | true == true
@assert isna(NA | false)
@assert isna(NA | NA)
@assert true | NA == true
@assert isna(false | NA)

@assert isna(NA & true)
@assert NA & false == false
@assert isna(NA & NA)
@assert isna(true & NA)
@assert false & NA == false

@assert any(DataVector[1, 2, NA] .== 1) == true
@assert any(DataVector[NA, 1, 2] .== 1) == true
@assert isna(any(DataVector[1, 2, NA] .== 3))
@assert any(DataVector[1, 2, 3] .== 4) == false

@assert isna(all(DataVector[1, 1, NA] .== 1))
@assert isna(all(DataVector[NA, 1, 1] .== 1))
@assert all(DataVector[1, 1, 1] .== 1) == true
@assert all(DataVector[1, 2, 1] .== 1) == false


text:
mutable struct Interval
    index::Int
    xmin::Float64
    xmax::Float64
    label::AbstractString
end

mutable struct Tier
    num::Int
    class::AbstractString
    name::AbstractString
    size::Int
    xmin::Float64
    xmax::Float64
    contents::Vector{Interval}
end

duration(interval::Interval) = interval.xmax - interval.xmin
duration(tier::Tier) = tier.xmax - tier.xmin

TextGrid = Vector{Tier}

function Base.show(io::IO, ::MIME"text/plain", tg::TextGrid)
    xmin, xmax = round.(extrema(tg), digits = 3)
    println("$(length(tg))-tier TextGrid ($(xmin)-$(xmax)s)")

    if length(tg) > 1
        print("Tiers:")
    elseif length(tg) == 1
        print("Tier:")
    end

    for tier in tg
        print(" ($(tier.num)) $(tier.name)")
    end
    line_width = displaysize(io)[2] < 80 ? (displaysize(io)[2] - 5) : 75
    length(tg) > 0 && print_tiers(io, tg, line_width)
end

function print_tiers(io::IO, tg::TextGrid, width = 75)
    println(io)
    println(io, "   +" * '-'^width * "+")

    if get(io, :limit, false) && length(tg) > 6
        for tier in tg[1:3]
            print_indiv_tier(io, tier, width)
        end
        println("⋮" * ' '^(3 + floor(Int, width/2)), "⋮")
        for tier in tg[end-1:end]
            print_indiv_tier(io, tier, width)
        end
    else
        for tier in tg
            print_indiv_tier(io, tier, width)
        end
    end
    print(io, "   +" * '-'^width * "+")
end

function print_indiv_tier(io::IO, tier::Tier, width = 75)
    print(io, rpad(tier.num, 3, ' ') * "|")

    chunk_width = duration(tier) / width
    chunk = fill(' ', width)
    if length(tier.contents) > 1
        in_chunk = unique([ceil(Int, (x.xmin - tier.xmin) / chunk_width) for x in tier.contents[2:end]])
        chunk[in_chunk] .= '|'
    end
    print(io, String(chunk))

    println(io, "|")
end


text:
prog = "1 + 1"
ex1 = Meta.parse(prog)
print(typeof(ex1))

max_stars_repo_name: zqsheng/snippet

max_issues_repo_name: zqsheng/snippet

max_forks_repo_name: zqsheng/snippet


text:
using Interact
using Blink
#using Plots
using WAV
using PortAudio
include("ColorCode.jl")

function plotBelief(belief::Belief)
  plot(OrderedDict(reverse([(string(k),v) for (k,v) in belief.b],dims=1)),
       plot_title="Belief",
       ylabel="Key",
       xlabel="Probability",
       seriestype=:bar,
       #bar_width=0.6,
       orientation=:horizontal,
       xlims=(0,1),
       yticks=:all,
       legend=false)
       #size=(1100,500))
  gui()
end

struct KeyboardFrontEnd
  letters::OrderedDict{Symbol,Widget{:latex,String}}
  commString::Widget{:latex,String}
end

function layoutKeyboard(keyboard::KeyboardFrontEnd)
  letterObjs = collect(values(keyboard.letters))
  return pad(1em,vbox(pad(0.5em,keyboard.commString),
                      vbox(pad(0.5em,hbox(map(x->pad(["right"],0.6em,x),letterObjs[1:9])...)),
                           pad(0.5em,hbox(map(x->pad(["right"],0.6em,x),letterObjs[10:18])...)),
                           pad(0.5em,hbox(map(x->pad(["right"],0.6em,x),letterObjs[19:26])...)),
                           pad(0.5em,hbox(map(x->pad(["right"],0.6em,x),letterObjs[27:end])...)))))
end

function renderAssignment(keyboard::KeyboardFrontEnd, assignment::Dict{Symbol,Int})
  colors = ["red" "blue"]
  foreach(x->keyboard.letters[x[1]][] = "\\fbox{\\color{$(colors[x[2]])} \\Huge \\texttt{$(keyboardStrings[x[1]])}}",assignment)
end

buttons = [button(x,style=Dict(:backgroundColor=>x)) for x in ["red" "blue"]]

keyboard = KeyboardFrontEnd(OrderedDict([sym => latex(str) for (sym,str) in keyboardStrings]),latex("\\Large\\textrm{|}"))

lmModel = LM.getModel()
lmState = LM.getStartState(lmModel)
prior = getPrior(lmModel,lmState)
commString = Ref("")
belief = Belief(prior,99,1)
history = BeliefHistory()
certaintyThreshold = 0.95

assignment = Dict([(k,1) for k in keys(keyboardStrings)])
changeAssignment(belief,assignment)
renderAssignment(keyboard,assignment)

#plotBelief(belief)

popSound,fs = wavread("pop.wav");
audioStream = PortAudioStream(0,2,adjust_channels=true)
write(audioStream,popSound) #stops lag for some reason does it? idk

function buttonCallback(button,commString)
  updateBelief(belief,button,assignment)
  newCommString = chooseLetter(belief,commString[],certaintyThreshold,history,lmModel,lmState)
  if length(newCommString) != length(commString[])
    @async write(audioStream,popSound)
  end
  commString[] = newCommString
  changeAssignment(belief,assignment)
  keyboard.commString[] = "\\Large\\textrm{$(replace(commString[]," "=>"\\  "))|}"

  #plotBelief(belief)

  renderAssignment(keyboard,assignment)
end

buttonCallbacks = [on(_->buttonCallback(n,commString),button) for (button,n) in zip(buttons,1:length(buttons))]

w = Window()
title(w,"ColorCode")
body!(w,vbox(layoutKeyboard(keyboard),pad(["left"],6em,hbox(pad(1em,buttons[1]), pad(1em,buttons[2])))))
js(w,Blink.JSString("""document.onkeydown = function (e) {Blink.msg("press",e.keyCode)}; """))
handle(w,"press") do key
  if key == 37
    buttonCallback(1,commString)
  end
  if key == 39
    buttonCallback(2,commString)

  end
end

wait()
LM.releaseModel(lmModel)
LM.releaseState(lmModel)

max_stars_repo_name: mrdaly/ColorCode

max_issues_repo_name: mrdaly/ColorCode

max_forks_repo_name: mrdaly/ColorCode


text:
#inspired by flatwhite syntax theme for Atom
const sand_palette = [
    colorant"#6494ED", # blue
    colorant"#73C990", # green
    colorant"#E2C08D", # brown
    colorant"#FF6347", # red
    colorant"#2E2C29", # dark,
    colorant"#4B4844" # medium
]

sand_bg = colorant"#F7F3EE"

_themes[:sand] = PlotTheme(
    bg = sand_bg,
    bginside = colorant"#E2DCD4",
    fg = colorant"#CBBFAF",
    fgtext = colorant"#725B61",
    fgguide = colorant"#725B61",
    fglegend = colorant"#725B61",
    palette = expand_palette(sand_bg, sand_palette),
    colorgradient = :dense
)


text:
using ChromNet

# project_groupgm
IC = [1 0.1 0.4 0.2; 0.1 1 0.3 0.5; 0.4 0.3 1 0.6; 0.2 0.5 0.6 1]
header = ["A1", "A2", "B1", "B2"]
groups = build_groups(inv(IC), header)
G,headerG = build_groupgm(IC, header, groups)

ans = [
    1.0 0.9 0.6 0.6;
    0.9 1.0 0.9 0.9;
    0.6 0.9 1.0 0.6;
    0.6 0.9 0.6 1.0
]
@test all((ChromNetPaper.project_groupgm(G, header, groups, x->true) .- ans) .< 1e-8)

ans = [
    1.0 0.1 0.6 0.6;
    0.1 1.0 0.8 0.8;
    0.6 0.8 1.0 0.6;
    0.6 0.8 0.6 1.0
]
@test all((ChromNetPaper.project_groupgm(G, header, groups, x->x[1] < 0.1) .- ans) .< 1e-8)

# run with default group filter
ChromNetPaper.project_groupgm(G, header, groups)


text:
"""
    count_nucleotides(strand)

The frequency of each nucleotide within `strand` as a dictionary.

Invalid strands raise a `DomainError`.

"""
function count_nucleotides(strand)
    d = Dict(n => 0 for n ∈ ['A','C','G','T'])
    for nucleotide ∈ strand
        if nucleotide ∉ keys(d)
            throw(DomainError("$nucleotide","Invalid nucleotide"))
        end
        d[nucleotide] += 1
    end
    return d
end

function count_nucleotides2(strand)
    ca = cc = cg = ct = 0
    for nuc in strand
        if nuc == 'A'
            ca += 1
        elseif nuc == 'C'
            cc += 1
        elseif nuc == 'G'
            cg += 1
        elseif nuc == 'T'
            ct += 1
        else
            throw(DomainError(nuc, "only A, C, G and T are valid nucleotides"))
        end
    end
    return Dict('A' => ca, 'C' => cc, 'G' => cg, 'T' => ct)
end

function count_nucleotides3(strand)
    counts = Dict((base => count(==(base), strand)) for base in "ACGT")
    if sum(values(counts)) != length(strand)
        throw(DomainError(strand, "only A, C, G and T are valid nucleotides"))
    end
    counts
end

function count_nucleotides4(strand)
    utf8 = transcode(UInt8, strand)
    # println(utf8)
    counts = zeros(Int, 256)
    @inbounds for byte in utf8
        counts[byte + 1] += 1
    end
    # println(counts)
    result = Dict(base => counts[Int(base) + 1] for base in "ACGT")
    if sum(values(result)) != length(utf8)
        throw(DomainError(strand, "only A, C, G and T are valid nucleotides"))
    end
    result
end

# println(count_nucleotides4("ACGTACGTACGT"))

max_stars_repo_name: vladflore/exercism

max_issues_repo_name: vladflore/exercism

max_forks_repo_name: vladflore/exercism


text: abstract type Boundary{N} end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
datapath = joinpath(dirname(pathof(Lasso)),"..","test","data")
plotspath = joinpath(dirname(@__FILE__), "plots")
mkpath(plotspath)

Random.seed!(243214)
@testset "plot GammaLassoPath's" begin
    @testset "$family" for (family, dist, link) in (("gaussian", Normal(), IdentityLink()), ("binomial", Binomial(), LogitLink()), ("poisson", Poisson(), LogLink()))
        data = CSV.read(joinpath(datapath,"gamlr.$family.data.csv"))
        y = convert(Vector{Float64},data[1])
        X = convert(Matrix{Float64},data[2:end])
        (n,p) = size(X)
        @testset "γ=$γ" for γ in (0, 2, 10)
            @testset "x=$x" for x in (:segment, :λ, :logλ)
                fitname = "gamma$γ.pf1"
                if !isfile(joinpath(datapath,"gamlr.$family.$fitname.params.csv"))
                    # file names in older Lasso packages
                    fitname = "gamma$γ"
                end
                # get gamlr.R params and estimates
                params = CSV.read(joinpath(datapath,"gamlr.$family.$fitname.params.csv"))
                fittable = CSV.read(joinpath(datapath,"gamlr.$family.$fitname.fit.csv"))
                gcoefs = convert(Matrix{Float64},CSV.read(joinpath(datapath,"gamlr.$family.$fitname.coefs.csv")))
                family = params[1,Symbol("fit.family")]
                γ=params[1,Symbol("fit.gamma")]

                # fit julia version
                glp = fit(GammaLassoPath, X, y, dist, link; γ=γ, λminratio=0.001)

                # test plots
                p = plot(glp;x=x)
                filename = joinpath(plotspath,"$family.$fitname.$x.svg")
                savefig(p, filename)
                @test isfile(filename)
            end
        end
    end
end


text:
###########################################################

"""
$(TYPEDSIGNATURES)

Calculate the area of each individual Beta function for a given protein concentration.
The SDS concentration axis needs a unit conversion which is handled by the `volumescal`
arguments.

# Arguments

- `betamix`: structure of type `BetaMix2D` containing the parameters of the Beta functions,
            the modes, confidence parameters and amplitudes and protein name
- `protcon`: protein concentration value (y axis) at which to perform the calculation of the area
- `volumescal`: scaling factor in μl to convert from mM to mole, instrument dependent. The default
               is 203.0μl.
"""
function area_enthalpy(betamix::BetaMix2D,protcon::Real ; volumescal::Real=203.0)

    # @assert betamix.betpar.ymin<=protcon<=betamix.betpar.ymax

    # scale the SDS concentration from mM to mole 10^(-3)*203*10^(-6) =2.03*10^(-7).
    # the 203 µl is specific for certain instruments
    scalfact = volumescal*10^(-9)

    ##------------------------------------
    ## Define a new ScaledBeta2DParams with the new units
    betpar_inp = deepcopy(betamix.betpar)
    betpar_scal = ScaledBeta2DParams(modefuny=betpar_inp.modefuny,konfuny=betpar_inp.konfuny,
                                     ampfuny=betpar_inp.ampfuny,a=betpar_inp.a*scalfact,b=betpar_inp.b*scalfact,
                                     ymin=betpar_inp.ymin,ymax=betpar_inp.ymax)

    ## Define a new BetaMix2D with the new units
    modkonamp_scal = copy(betamix.modkonamp)
    # scale the mode in terms of SDS concentration:
    modkonamp_scal[1:2,:] .*= scalfact
    betamix_scal = BetaMix2D(betpar_scal,modkonamp_scal,betamix.protein)
    ##------------------------------------

    ncomp = size(betamix_scal.modkonamp,2)
    integr = zeros(ncomp)
    err = zeros(ncomp)
    for c=1:ncomp

        ## betamix_scal  _scal!!!
        bp = betamix_scal.betpar
        mcur = betamix_scal.modkonamp[:,c]

        @assert length(mcur)==bp.nummodpar

        ##---------------------------------------------------------
        # (val,err) = hquadrature(f::Function, xmin::Real, xmax::Real;
        #                     reltol=1e-8, abstol=0, maxevals=0)
        abeta = bp.a
        bbeta = bp.b
        ## get the values of model parameters at y=ycur
        mode,kon,amp = getmodparbeta(bp,mcur,protcon)
        #@show c,mode,kon,amp,protcon

        integr[c],err[c] = hquadrature(x->scaledbeta(mode,kon,abeta,bbeta,amp,x),abeta,bbeta,
                                       reltol=1e-8,abstol=0,maxevals=0)

        # ## gaussian quadrature check
        # Npts = 500
        # qtpts,weights = gauss(Npts,abeta,bbeta)
        # nodes = [scaledbeta(mode,kon,abeta,bbeta,amp,x) for x in qtpts]
        # gaussquad = dot( weights, nodes )
        # @show gaussquad-integr[c]

    end

    return integr,err
end

###########################################################

"""
$(TYPEDSIGNATURES)

Calculate the area of each individual Beta function for a set of  protein concentrations.

# Arguments

- `betamix`: structure of type `BetaMix2D` containing the parameters of the Beta functions,
            the modes, confidence parameters and amplitudes and protein name
- `protcosn`: array of protein concentrations at which to calculate the areas
- `volumescal`: scaling factor in μl to convert from mM to mole, instrument dependent. The default
               is 203.0μl.
"""
function areasvsprotcon(betamix::BetaMix2D,protcons::Vector{Float64},outdir::String,protein::String;
                        volumescal::Real=203.0)

    ##-------------------------------------------------------
    N=length(protcons)
    ncomp = size(betamix.modkonamp,2)
    areas = zeros(N,ncomp)
    errar = zeros(N,ncomp)
    for i=1:N
        areas[i,:],errar[i,:] = area_enthalpy(betamix,protcons[i],volumescal=volumescal)
    end

    ##-------------------------------------------------------
    ## linear fit to area curves
    angcoe = Vector{Float64}(undef,ncomp)
    intercept = Vector{Float64}(undef,ncomp)
    resstdev = Vector{Float64}(undef,ncomp)
    for i=1:ncomp
        angcoe[i],intercept[i],resstdev[i] = lssqregr([areas[:,i] protcons])
    end
    linfitres = [angcoe intercept]

    ##-------------------------------------------------------
    # save in JLD2 format
    outfile1 = joinpath(outdir,protein*"_ITCareasprotcon.jld2")
    println("Saving areas data in JLD2 to $outfile1\n")
    jldopen(outfile1, "w") do fl
        fl["linearfit"] = linfitres
        fl["areas"] = areas
        fl["errar"] = errar
        fl["resstdev"] = resstdev
        fl["protconareas"] = protcons
        fl["volumescal"] = volumescal
    end

    # ##-------------------------------------------------------
    # # save in plain text format
    # outdata = [areas protcons]
    # outfile2 = joinpath(outdir,protein*"_ITCareasprotcon.dat")
    # println("Saving results in text file to $outfile2\n")
    # header = "# Binding isotherm results for protein $protein, columns: freeSDS and Nbound"
    # ofl = open(outfile2,"w")
    # writedlm(ofl,[header])
    # writedlm(ofl,outdata)
    # close(ofl)

    return areas,errar,linfitres,resstdev
end

###########################################################

"""
$(TYPEDSIGNATURES)

Calculate the volume of each single Beta function within given bounds of protein concentration.
The SDS concentration axis needs a unit conversion which is handled by the `volumescal`
arguments.

# Arguments

- `betamix`: structure of type `BetaMix2D` containing the parameters of the Beta functions,
            the modes, confidence parameters and amplitudes and protein name
- `minprotcon`: lower bound of protein concentration value (y axis) to perform the integral
- `maxprotcon`: upper bound of protein concentration value (y axis) to perform the integral
- `volumescal`=203.0: scaling factor in μl to convert from mM to mole, instrument dependent. The default
               is 203.0μl.
"""
function volume_enthalpy(betamix::BetaMix2D,minprotcon::Real,maxprotcon::Real ; volumescal::Real=203.0)

    @assert betamix.betpar.ymin<=minprotcon<=betamix.betpar.ymax
    @assert betamix.betpar.ymin<=maxprotcon<=betamix.betpar.ymax
    @assert minprotcon<maxprotcon

    # scale the SDS concentration from mM to mole 10^(-3)*203*10^(-6) =2.03*10^(-7).
    # the 203 µl is specific for certain instruments
    scalfact = volumescal*10^(-9)

    ##------------------------------------
    ## Define a new ScaledBeta2DParams with the new units
    betpar_inp = deepcopy(betamix.betpar)
    betpar_scal = ScaledBeta2DParams(modefuny=betpar_inp.modefuny,konfuny=betpar_inp.konfuny,
                                     ampfuny=betpar_inp.ampfuny,a=betpar_inp.a*scalfact,b=betpar_inp.b*scalfact,
                                     ymin=betpar_inp.ymin,ymax=betpar_inp.ymax)

    ## Define a new BetaMix2D with the new units
    modkonamp_scal = copy(betamix.modkonamp)
    # scale the mode in terms of SDS concentration:
    modkonamp_scal[1:2,:] .*= scalfact
    betamix_scal = BetaMix2D(betpar_scal,modkonamp_scal,betamix.protein)
    ##------------------------------------

    ##---------------------------------------------------
    function betaval(xy::Vector{<:Real})::Real
        ## get the values of model parameters at y=ycur
        x=xy[1]
        y=xy[2]
        mode,kon,amp = getmodparbeta(bp,mcur,y)
        ## calculate value of beta function
        bval = scaledbeta(mode,kon,bp.a,bp.b,amp,x)
        return bval
    end
    ##---------------------------------------------------

    ##---------------------------------------------------------
    # (val,err) = hquadrature(f::Function, xmin::Real, xmax::Real;
    #                     reltol=1e-8, abstol=0, maxevals=0)
    bp = betamix_scal.betpar
    xmin = [bp.a,minprotcon]
    xmax = [bp.b,maxprotcon]

    ncomp = size(betamix_scal.modkonamp,2)
    integr = zeros(ncomp)
    err = zeros(ncomp)
    mcur=nothing # given the closure of betaval() and the private scope of for loop
    for c=1:ncomp
        mcur = betamix_scal.modkonamp[:,c]
        integr[c],err[c] = hcubature(betaval,xmin,xmax,reltol=1e-8,abstol=0,maxevals=0)
    end

    return integr,err
end

############################################################


text:
######################################################################
# greedyp.jl
#
# like greedy.jl, but includes penalty for being too close to target
######################################################################
mutable struct GreedyP <: Policy
    actions::Vector{Action}
    lambda::Float64
    d::Float64
end
function GreedyP(max_step::Real, n::Int, headings=0; stay=true)
    GreedyP(make_action_list(max_step, n, headings, stay=stay), 10.0, 15.0)
end
function GreedyP(x::Vehicle, n::Int, headings=0; stay=true)
    GreedyP(make_action_list(x.max_step, n, headings, stay=stay), 10.0,15.0)
end

# loop over all actions.
# The one with smallest expected entropy is best
# We know this depends on sensor
function action(m::SearchDomain, x::Vehicle, o, f::AbstractFilter, p::GreedyP)
    best_score = -Inf
    best_a = (0.0, 0.0, 0.0)
    best_mi = 0.0
    best_p = 0.0

    # propagate target dynamics forward
    f2 = deepcopy(f)
    predict!(f2)

    for a in p.actions
        # find out where a will take you
        xp = new_pose(m, x, a)
        # compute best mutual information
        mi = mutual_information(f2, xp)

        # compute penalty for being too close
        penalty = p.lambda * nmac_penalty(f2, xp, p.d)

        score = mi - penalty

        if score > best_score
            best_mi = mi
            best_p = penalty
            best_score = score
            best_a = a
        end
    end

    #println("best_mi = ", best_mi, ", best_p = ", best_p)
    #c = covariance(f2)
    #ev = eigvals(c)
    #e1 = ev[1]
    #e2 = ev[2]
    #println("c = ", covariance(f))
    #println("e1 = ", e1, ", e2 = ", e2)
    #println("ce = ", cheap_entropy(f, 200.0, 40))
    return best_a
end

function nmac_penalty(pf::PF, xp::Pose, d::Float64)
    nmac_penalty(pf._b, xp, d, pf.n)
end
function nmac_penalty(b::ParticleCollection, xp::Pose, d::Float64, n::Int)
    penalty = 0.0
    for i = 1:n
        tp = particle(b,i)
        dx = tp[1] - xp[1]
        dy = tp[2] - xp[2]

        if dx*dx + dy*dy < d*d
            penalty += 1.0
        end
    end
    return penalty / n
end


text:
##
## Find an approximate gcd (agcd) for two polynomials
## The `gcd` function can exactly handle Poly{Int} if Rational{BigInt} is used:
## p = poly([1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,3,3])
## q = Poly(map(u->big(u)//1, coeffs(p)), p.var)
## gcd(p, polyder(p)) # Poly(5.960464477539063e-8)
## gcd(q, polyder(q)) # degree 17, as expected
## But what to do with Poly{Float}? There floating point issues will be a problem,
## This create `agcd` for approximate gcd:
## u,v,w,err = agcd(p, polyder(p))
## degree(u)  # 17
## roots(v)   # 1,2,3

## reimplement Zeng's gcd code for approximate gcds

## Special matrices
function cauchy_matrix{T}(p::Poly{T}, k::Integer)
    n = degree(p) + 1
    out = zeros(T, n + k - 1, k)
    for i in 1:k
        out[(1:n) + (i-1), i] = rcoeffs(p)
    end
    out
end

function sylvester_matrix(p::Poly, q::Poly; k::Int=0)
    @assert k >= 0
    n,m = degree(p), degree(q)
    if n < m
        p,q = q,p
        n,m = m,n
    end

    del = n - m
    i,j = k + del, k
    if k == 0
        i,j = n,m
    end
    hcat(cauchy_matrix(q, i), cauchy_matrix(p, j))
end

## preconditioning code
## taken from https://who.rocq.inria.fr/Jan.Elias/pdf/je_masterthesis.pdf
function geometric_mean(a::Vector, epsilon=Base.eps())
    a = filter(x -> abs(x) > epsilon, a)
    n = length(a)
    prod(abs(a) .^ (1/n))
end

function ratio(p,q, atol=Base.eps(), rtol=Base.eps())
    is_nonzero(x) = !isapprox(x, 0; rtol=rtol, atol=atol)
    as = abs(filter(is_nonzero, p.a))
    length(as) == 0 && return Inf

    bs = abs(filter(is_nonzero, q.a))
    length(bs) == 0 && return Inf

    max(maximum(as), maximum(bs)) / min(minimum(as), minimum(bs))
end

## find alpha, gamma that minimize ratio of max coefficients to min coefficients, for getting zeros
## 1.12 writes this as a linear programming problem, we just ...
function precondition(p::Poly,q::Poly)
    alphas = [2.0^i for i in -5:5]
    phis = [2.0^i for i in -5:5]
    out = (1,1)
    m = ratio(p,q)

    for α in alphas
        for ϕ in phis
            r = ratio(polyval(p, ϕ * variable(p)), α * polyval(q, ϕ * variable(q)))
            if r < m
                out = (α, ϕ)
            end
        end
    end

    α, ϕ = out

    p = polyval(p, ϕ * variable(p))
    q = α * polyval(q, ϕ * variable(q))

    p = p * (1/geometric_mean(coeffs(p)))
    q = q * (1/geometric_mean(coeffs(q)))

    p,q, ϕ, α

end

## converge on right singular vector of A associated with singular value sigma (not a double value)
## we return if sigma < tol; delta \appro 0;
## how to @assert that sigma_2 > sigma_1?
function lemma24(A::Matrix; θ::Real=1e-8)
    const MAXSTEPS = 100

    Q,R = Base.qr(A)
    if rank(R) < size(R)[2]
        λs, vs = eig(R)
        _, ind = findmin(abs(λs))
        return(λs[ind], vs[:,ind])
    end

    x = rand(size(A)[2]); x/norm(x,2)   # use random initial guess
    σ, σ1 = 1e8, Inf

    function update(x)
        y = conj(R)' \ x
        z = R \ y
        x = z/norm(z,2)
        sigma = norm(R * x, 2)  # y/norm(z,2)
        (x, minimum(abs(sigma)))
    end

    ## how long do we update? Size is only issue, not accuracy so we iterate until we stop changing much
    for _ in 1:MAXSTEPS
        x, σ1 = update(x)
        ((σ1 < θ) | (abs((σ - σ1) / σ1) < 1.1)) && break
        σ = σ1
    end

    return(σ1, x)
end

"""
Return k, u,v,w where k reveals rank; u*v \approx p; u*w \approx q; v & w coprime

Following Zeng, section 4, this could be made more efficient.
In Zeng, theta=1e-8. Here agcd uses 1e-12, which is an improvement on some...
"""
function reveal_rank(p::Poly, q::Poly, theta=1e-8)
    n,m = map(degree, (p, q))
    if n < m
        p,q = q,p
        n,m=m,n
    end

    for k in 1:m
        A = sylvester_matrix(p,q,k=k)
        psi, x = lemma24(A)
        if abs(psi) < norm(p,2) * theta #  norm(A,2)?
            ## degree u = n - k; degree(v) = k
            v = monic(Poly(reverse(x[1:(length(x)-k)]), p.var))
            w = monic(Poly(reverse(x[(length(x)-k+1):end]), p.var))
            u = monic(Poly(reverse(cauchy_matrix(v, degree(p) - degree(v) + 1) \ rcoeffs(p)), p.var))
            return (k, u, v,w)
        end
    end
    return (n+m, Poly(ones(eltype(coeffs(p)),1)), p, q)
end

## Lemma 4.1, (25)
## solve weighted least squares problem W*(Ax - b) = 0
function weighted_least_square(A, b, w)
    W = diagm(w)
    (W * A) \ (W * b)
end

## Jacobian F(u,v,w) = [p,p'] is J(u,v,w)
function JF{T}(u::Poly{T}, v::Poly, w::Poly)
    j, k, l= degree(u), degree(v), degree(w)
    n, m = j + k, j + l

    a = cauchy_matrix(v, n+1-k)
    b = cauchy_matrix(u, n+1-j)
    c = zeros(T, n+1, m+1-j)

    d = cauchy_matrix(w, m+1-l)
    e = zeros(T, m+1, n+1-j)
    f = cauchy_matrix(u, m+1-j)

    A = hcat(a,b,c)
    B = hcat(d,e,f)
    vcat(A,B)
end
## compute F(u,v,w) - [p, p'] = [u*v, u*w] - [p, p']
Fmp(p,q,u,v,w) = [rcoeffs(u*v); rcoeffs(u*w)] - [rcoeffs(p); rcoeffs(q)]
## error in estimate for p=u*v,q=u*w for some weights
residual_error(p,q,u,v,w, wts=ones(degree(p) + degree(q) + 2)) = norm(Fmp(p,q,u,v,w) .* wts, 2)
function agcd_update(p, q, u, v, w, wts)
    m,n = map(degree, (u,v))
    A = JF(u, v, w)
    b = Fmp(p,q,u,v,w)
    inc = weighted_least_square(A, b, wts)

    x = vcat(rcoeffs(u), rcoeffs(v), rcoeffs(w))
    x = x - inc

    u = Poly(reverse(x[1:(1+m)]), p.var)
    v = Poly(reverse(x[(m+2):(m+n+2)]), p.var)
    w = Poly(reverse(x[(m+2+n+1):end]), p.var)

    err = residual_error(p,q,u,v,w, wts)
    (u, v, w, err)
end

"""

Find an approximate GCD for polynomials `p` and `q` using an algorithm of [Zeng](http://www.ams.org/journals/mcom/2005-74-250/S0025-5718-04-01692-8/home.html).

Returns u,v,w, err where:

* `u*v \approx monic(p)`
* `u*w \approx monic(q)`
* The total residual error in these approximations is bounded by `err`.

Further,
* `v` and `w` should share no common roots (`u` is a gcd of `u*v` and `u*w`)
* `roots(v)` should exhaust unique values of `roots(p)`.

The tolerances are:

* theta: passed to the reveal_rank function. In Zeng 1e-8 is used. Here 1e-12 seems to work better?
* \rho: if the residual error does not come below this mark, then we use the initial guess

"""
function agcd{T,S}(p::Poly{T}, q::Poly{S}=p';
              theta = 1e-12,    # reveal_rank tolerance. (1e-8 in paper, this seems better?)
              ρ::Real = 1e-10   # residual tolerance
              )

    n, m = map(degree,(p,q))
    if m > n
        p,q=q,p
    end

    if m == 0
        return (Poly(ones(1)), p, q, 0)
    end

    p0,q0 = map(copy, (p,q))
    p,q, phi, alpha = precondition(p,q)

    k, u, v, w = reveal_rank(p, q, theta)
    u0,v0,w0 = map(copy, (u,v,w))

    m,n,k = map(degree, (u, v, w))
    wts =  map(pj -> 1/max(1, abs(pj)), vcat(rcoeffs(p), rcoeffs(q)))

    ## iterate to solve
    ## uvw_j+1 = uvw_j - Jw[F - p]

    err0, err1 = Inf, residual_error(p,q,u,v,w, wts)
    for ctr in 1:20
        try
            u1, v1, w1, err1 = agcd_update(p, q, u, v, w, wts)
            if err1 < err0
                err0, u, v, w = err1, u1, v1, w1
            else
                break
            end
        catch err
            break                       # possibly singular
        end
    end

    if err0 > norm(sylvester_matrix(p,q), 2) * ρ
        u,v,w = map(monic, (u0, v0, w0))  ## failed to converge, so we return the initial guess
    end

    x = (1/phi) * variable(p) # reverse preconditioning
    u,v,w = map(monic, (polyval(u, x), polyval(v, x), polyval(w, x)))

    (u,v,w, residual_error(p0,q0,u,v,w))
end


text:
module Physt

import Base.<<, Base.deepcopy, Base.copy
export h1, update, update!, find_bin, 📊

abstract type Binning{AxisType} end
abstract type Histogram end

struct RegularBinning{AxisType} <: Binning{AxisType}
    nbins::Int64
end

struct Histogram1D{AxisType, ValueType} <: Histogram where AxisType
    edges::Array{AxisType, 1}
    values::Array{ValueType, 1}
end

function copy(old::Histogram1D{AxisType, ValueType}) where AxisType where ValueType
    new = Histogram1D{AxisType, ValueType}(copy(old.edges), copy(old.values))
end

deepcopy(histogram::Histogram) = copy(histogram)

function fit(binning::RegularBinning{AxisType}, entries::AbstractArray)::Array{AxisType, 1} where AxisType
    min = convert(AxisType, minimum(entries))
    max = convert(AxisType, maximum(entries))
    step = (max - min) / binning.nbins
    return collect(range(min, length=binning.nbins + 1, step=step))
end

@inline function find_bin(edges::AbstractArray{AxisType}, value::AxisType) where AxisType
    rightedge = findfirst(e -> e > value, edges)
    if rightedge == nothing
        # Last edge is inclusive
        if value == edges[end]
            return length(edges) - 1
        end
    elseif rightedge > 1
        return rightedge - 1
    end
    return nothing
end

function find_bin(histogram::Histogram1D, value)
    return find_bin(histogram.edges, value)
end

function transform(edges::AbstractArray{AxisType}, entries::AbstractArray{AxisType})::Array{Int64, 1} where AxisType
    # Least effective of all implementations
    nbins = length(edges) - 1
    result = zeros(Int64, nbins)
    for value in entries
        bin = find_bin(edges, value)
        if bin !== nothing
            result[bin] += 1
        end
    end
    return result
end

function fit_transform(binning::Binning, entries::AbstractArray)
    edges = fit(binning, entries)
    return edges, transform(edges, entries)
end

function update(histogram::Histogram, any)
    new = copy(histogram)
    update!(new, any)
    return new
end

function update!(histogram::Histogram1D{AxisType}, entry) where AxisType
    entry_converted = convert(AxisType, entry)
    bin = find_bin(histogram.edges, entry_converted)
    if bin != nothing
        histogram.values[bin] += 1
    end
    return histogram
end

function update!(histogram::Histogram1D{AxisType}, entries::AbstractArray) where AxisType
    new_values = transform(histogram.edges, entries)
    histogram.values .+= new_values
    return histogram
end

<<(histogram::Histogram, what_to_update) = update(histogram, what_to_update)

function h1(entries::AbstractArray{AxisType}, nbins::Integer=10) where AxisType
    binning = RegularBinning{AxisType}(nbins)
    edges, values = fit_transform(binning, entries)
    return Histogram1D{AxisType, Int64}(edges, values)
end

function h1(entries::AbstractArray{AxisType}, edges::AbstractArray{AxisType}) where AxisType
    values = transform(edges, entries)
    return Histogram1D{AxisType, Int64}(edges, values)
end

📊 = h1  # Because we can...

end # module


text:
# we also need this before defining a module
using DotOverload

@dotted module PlotTest

using DotOverload
using PyCall
using PyPlot

DotOverload.getMember(t::PyObject, k) = t[k]
DotOverload.setMember!(t::PyObject, k, v) = t[k] = v
DotOverload.getMember(t::Figure, k) = t[k]
DotOverload.setMember!(t::Figure, k, v) = t[k] = v

plt.ioff()

fig, ax = plt.subplots()
line, = ax.plot([],[])

plt.ion()

function init()
  fig.show()
  fig.canvas.draw()
  plt.axis([0, 2*π, -2, 2])
end

function doStuff()
  x = linspace(0,2*pi,1000);
  y = sin(3*x + 4*cos(2*x))

  line.set_xdata(x)
  line.set_ydata(y)

  ax.draw_artist(line)
  fig.canvas.update()
end

end #module

# usage:
# using PlotTest
# PlotTest.init()
# PlotTest.doStuff()


text:
mutable struct DiagonalPreconditioner{T, S<:AbstractVector{T}} <: AbstractPreconditioner
    D::S
end
function EmptyDiagonalPreconditioner(A::AbstractMatrix)
    D = zeros(eltype(A), size(A,1))
    return DiagonalPreconditioner{eltype(D), typeof(D)}(D)
end
function DiagonalPreconditioner(A::AbstractMatrix)
    D = diag(A, 0)
    return DiagonalPreconditioner{eltype(D), typeof(D)}(D)
end
function diag!(D, A::AbstractMatrix)
    length(D) == size(A,1) == size(A,2) || throw("D and A sizes are not compatible.")
    @inbounds @simd for i in 1:length(D)
        D[i] = A[i,i]
    end
    return
end
function UpdatePreconditioner!(D::DiagonalPreconditioner, K::AbstractMatrix)
    diag!(D.D, K)
    return D
end
function ldiv!(y, C::DiagonalPreconditioner, b)
    @inbounds @simd for i in 1:length(C.D)
        y[i,:] .= view(b, i, :) ./ C.D[i]
    end
    return y
end
function (\)(C::DiagonalPreconditioner, b)
    y = zeros(b)
    @inbounds @simd for i in 1:length(C.D)
        y[i,:] .= view(b, i, :) / C.D[i]
    end
    return y
end


text:
#Factorials and elementary coefficients

export
    derangement,
    partialderangement,
    factorial,
    subfactorial,
    doublefactorial,
    hyperfactorial,
    multifactorial,
    gamma,
    primorial,
    multinomial

# TODO: This should really live in Base, otherwise it's type piracy
"""
    factorial(n, k)

Compute ``n!/k!``.
"""
function Base.factorial(n::T, k::T) where T<:Integer
    if k < 0 || n < 0 || k > n
        throw(DomainError((n, k), "n and k must be nonnegative with k ≤ n"))
    end
    f = one(T)
    while n > k
        f = Base.checked_mul(f, n)
        n -= 1
    end
    return f
end
Base.factorial(n::Integer, k::Integer) = factorial(promote(n, k)...)

"""
    derangement(n)

Compute the number of permutations of `n` with no fixed points, also known as the
subfactorial. An alias `subfactorial` for this function is provided for convenience.
"""
function derangement(sn::Integer)
    n = BigInt(sn)
    return numerator(factorial(n) * sum([(-1)^k // factorial(k) for k = 0:n]))
end
const subfactorial = derangement

function doublefactorial(n::Integer)
    if n < 0
        throw(DomainError(n, "n must be nonnegative"))
    end
    z = Ref{BigInt}(0)
    ccall((:__gmpz_2fac_ui, :libgmp), Cvoid, (Ref{BigInt}, UInt), z, UInt(n))
    return z[]
end

"""
    partialderangement(n, k)

Compute the number of permutations of `n` with exactly k fixed points.
"""
function partialderangement(n::Integer, k::Integer)
    if n < 0
        throw(DomainError(n))
    end
    if k < 0 || k > n
        throw(DomainError(k))
    end
    a = BigInt(n)
    b = BigInt(k)
    return subfactorial(a - b) * binomial(a, b)
end

# Hyperfactorial
hyperfactorial(n::Integer) = prod(i->i^i, BigInt(2):n)

function multifactorial(n::Integer, m::Integer)
    if n < 0
        throw(DomainError(n, "n must be nonnegative"))
    end
    z = Ref{BigInt}(0)
    ccall((:__gmpz_mfac_uiui, :libgmp), Cvoid, (Ref{BigInt}, UInt, UInt), z, UInt(n), UInt(m))
    return z[]
end

function primorial(n::Integer)
    if n < 0
        throw(DomainError(n, "n must be nonnegative"))
    end
    z = Ref{BigInt}(0)
    ccall((:__gmpz_primorial_ui, :libgmp), Cvoid, (Ref{BigInt}, UInt), z, UInt(n))
    return z[]
end

"""
    multinomial(k...)

Multinomial coefficient where `n = sum(k)`.
"""
function multinomial(k...)
    s = 0
    result = 1
    @inbounds for i in k
        s += i
        result *= binomial(s, i)
    end
    result
end


text:
mutable struct FIXMessageIterator
    data::Vector{UInt8}
    messages::OrderedDict{Int64, OrderedDict{Int64, String}}
    tag::Vector{Char}
    value::Vector{Char}
    state::Bool
    message_id::Int64
    function FIXMessageIterator(data::Vector{UInt8})
        return new(data, OrderedDict{Int64, OrderedDict{Int64, String}}(), Char[], Char[], true, 0)
    end
end

length(this::FIXMessageIterator) = length(this.data)
function start(this::FIXMessageIterator)::Int64
    return 1
end

function next(this::FIXMessageIterator, idx::Int64)::Tuple{FIXMessageIterator, Int64}
    x = this.data[idx]
    if x == 0x01
        int_tag = parse(Int64, String(this.tag))
        str_value = String(this.value)
        !haskey(this.messages, this.message_id) && (this.messages[this.message_id] = OrderedDict{Int64, String}())
        this.messages[this.message_id][int_tag] = str_value
        resize!(this.tag, 0)
        resize!(this.value, 0)
        if int_tag == 10
            this.message_id += 1
        end
        this.state = true
    elseif x == 0x3d
        this.state = !this.state
    else
        if this.state
            push!(this.tag, Char(x))
        else
            push!(this.value, Char(x))
        end
    end

    return (this, idx + 1)
end

function done(this::FIXMessageIterator, idx::Int64)
    return length(this.data) > idx
end

function fixparse(data::Vector{UInt8})::OrderedDict{Int64, OrderedDict{Int64, String}}
    f = FIXMessageIterator(data)
    for k in eachindex(data)
        next(f, k)
    end
    return f.messages
end

function fixconvert(tags::Dict{Int64, String}, msg::OrderedDict{Int64, String})::OrderedDict{String, String}
    out = OrderedDict{String, String}()
    for (k, v) in msg
        out[tags[k]] = v
    end
    return out
end


text:
push!(LOAD_PATH,"../src/")
push!(LOAD_PATH,"../docs/")

using Documenter, ISA

makedocs(;
    modules=[ISA],
    authors="Steven Sandoval, Hasan Al-Shammari, Mamta Dalal",
    sitename="ISA.jl",
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Instantaneous Spectral Analysis" =>
                    ["Cannonical Triplets and Component Sets"   => "basics/triplets.md",
                     "AM--FM Components"     => "basics/components.md",
                     "AM--FM Models"         => "basics/models.md",
                     "Instantaneous Spectra" => "basics/spectra.md",],
        "Basic Components" =>
                    ["Simple Harmonic Components"     => "basicComponents/SHC.md",
                     "AM Components"     => "basicComponents/AM.md",
                     "FM Components"     => "basicComponents/FM.md",
                     "AM--FM Components"     => "basicComponents/AMFM.md",],
        "Basic Signal Models" =>
                    ["Fourier Series"     => "basicModels/Fourier.md",
                     "AM--FM Models"     => "basicModels/AMFM.md",],
        "Numerical ISA" =>
                    ["Numerical AM--FM Components"     => "numerical/components.md",
                     "Numerical AM--FM Models"     => "numerical/models.md",
                     "Numerical Triplets and Sets"     => "numerical/tripletsSets.md",
                     "Numerical Instantaneous Spectra"     => "numerical/spectra.md",],
        "Types/Methods/Functions" => "typesMethodsFunctions.md",
        "Utility Functions" => "utilityFunctions.md",
        "Plotting Recipes" => "recipes.md",
        "Cite" => "cite.md",
    ],
    repo="https://github.com/NMSU-ISA/ISA.jl/blob/{commit}{path}#{line}",
    assets=String[],
)

deploydocs(;
    repo="github.com/NMSU-ISA/ISA",
    devurl = "dev"
)


text:
#background-anomaly separation escalator video with Minkowski decomposition
using Distributed
@everywhere using SetIntersectionProjection
using MAT
using LinearAlgebra
using Statistics

ENV["MPLBACKEND"]="qt5agg"
using PyPlot

#define computational grid structure
@everywhere mutable struct compgrid
  d :: Tuple
  n :: Tuple
end

#select working precision
TF=Float32 #Float32

#load escalator video) (Mablab files for compatibility with matlab only solvers for comparison...)
#load escalator video) (Mablab files for compatibility with matlab only solvers for comparison...)
if isfile("escalator_data.mat") == true
  file     = matopen("escalator_data.mat")
  #file = matopen(joinpath(dirname(pathof(SetIntersectionProjection)), "../examples/Data/escalator_data.mat"))
else
  println("downloading escalator video from http://cvxr.com/tfocs/demos/rpca/escalator_data.mat")
  run(`wget http://cvxr.com/tfocs/demos/rpca/escalator_data.mat`)
  #error("download escalator video from http://cvxr.com/tfocs/demos/rpca/escalator_data.mat")
end
mtrue    = read(file, "X")
n1       = convert(Integer,read(file, "m"))
n2       = convert(Integer,read(file, "n"))
m_mat    = convert(Array{TF,2},mtrue)
m_tensor = convert(Array{TF,3},reshape(mtrue,n1,n2,Integer(200)))

#preprocess data
for i=1:size(m_tensor,3)
  m_tensor[:,:,i].=m_tensor[:,:,i].-median(m_tensor[:,:,i])
end

#split data tensor in a training and evaluation data
m_train      = m_tensor[:,:,181:200] #last 20 frames do not have people -> this should be the background so we learn constraints from it
m_evaluation = m_tensor[:,:,1:180]
m_est        = zeros(TF,size(m_evaluation)) #allocate results

#computational grid for the training images (all images assumed to be on the same grid here)
comp_grid_eval  = compgrid((1f0,1f0,1f0),(size(m_evaluation,1),size(m_evaluation,2), size(m_evaluation,3)))
comp_grid_train = compgrid((1f0,1f0,1f0),(size(m_train,1),size(m_train,2), size(m_train,3)))

#initialize constraint dictionaries
constraint_background = Vector{SetIntersectionProjection.set_definitions}()
constraint_anomaly    = Vector{SetIntersectionProjection.set_definitions}()
constraint_sum        = Vector{SetIntersectionProjection.set_definitions}()

#observe and define bound constraints for background
  #bounds:
  m_min     = minimum(m_train,dims=3)
  m_min     = vec(repeat(m_min[:,:,:],outer=[1,1,comp_grid_eval.n[3]])) .- 40.0
  m_max     = maximum(m_train,dims=3)
  m_max     = vec(repeat(m_max[:,:,:],outer=[1,1,comp_grid_eval.n[3]])) .+ 40.0
  set_type  = "bounds"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint_background, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

  #subspace constraints on each z-slice of the tensor
  m_min    = 0.0
  m_max    = 0.0
  set_type = "subspace"
  TD_OP     = "identity"
  A = m_train
  A = reshape(A,size(A,1)*size(A,2),size(A,3))
  F = svd(A)
  app_mode  = ("slice","z")
  custom_TD_OP = (F.U,true)
  push!(constraint_background, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#set bound constraints for sum
  m_min     = minimum(m_tensor) .- 10f0
  m_max     = maximum(m_tensor) .+ 10f0
  set_type  = "bounds"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint_sum, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#define constraints for anomaly
  # #set bound constraints for anomaly (this should be bounds-sum - bounds-background)
  # m_min     = minimum(m_tensor) .- 10f0 .- constraint_background[1].min .- 20.0
  # m_max     = maximum(m_tensor) .- 10f0 .- constraint_background[1].max .+ 20.0
  # set_type  = "bounds"
  # TD_OP     = "identity"
  # app_mode  = ("matrix","")
  # custom_TD_OP = ([],false)
  # push!(constraint_anomaly, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

  #set maximum of nonzero pixels in anomaly component for every time-frame
  m_min     = 0
  m_max     = Integer(round(0.25*comp_grid_eval.n[1]*0.25*comp_grid_eval.n[2]))
  set_type  = "cardinality"
  TD_OP     = "identity"
  app_mode  = ("slice","z")
  custom_TD_OP = ([],false)
  push!(constraint_anomaly, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

  #set maximum of nonzero jumps (x) in pixel values in anomaly component for every time-frame
  m_min     = 0
  m_max     = 4*12*10
  set_type  = "cardinality"
  TD_OP     = "D_x"
  app_mode  = ("slice","z")
  custom_TD_OP = ([],false)
  push!(constraint_anomaly, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

  #set maximum of nonzero jumps (xy in pixel values in anomaly component for every time-frame
  m_min     = 0
  m_max     = 22*2*10
  set_type  = "cardinality"
  TD_OP     = "D_y"
  app_mode  = ("slice","z")
  custom_TD_OP = ([],false)
  push!(constraint_anomaly, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#define solver options
options=PARSDMM_options()
options.FL=Float32
options=default_PARSDMM_options(options,options.FL)
options.adjust_gamma           = true
options.adjust_rho             = true
options.adjust_feasibility_rho = true
options.Blas_active            = true
options.maxit                  = 100
options.feas_tol               = 0.001
options.obj_tol                = 0.001
options.evol_rel_tol           = 0.00001

set_zero_subnormals(true)
BLAS.set_num_threads(4)
options.parallel         = false
options.feasibility_only = false
options.zero_ini_guess   = true

#set up constraints, precompute some things and define projector
(P_sub_sum,TD_OP_sum,set_Prop_sum) = setup_constraints(constraint_sum,comp_grid_eval,options.FL)
(P_sub_c1,TD_OP_c1,set_Prop_c1)    = setup_constraints(constraint_background,comp_grid_eval,options.FL)
(P_sub_c2,TD_OP_c2,set_Prop_c2)    = setup_constraints(constraint_anomaly,comp_grid_eval,options.FL)

options.Minkowski = true
(TD_OP,set_Prop,AtA,l,y) = PARSDMM_precompute_distribute_Minkowski(TD_OP_c1,TD_OP_c2,TD_OP_sum,set_Prop_c1,set_Prop_c2,set_Prop_sum,comp_grid_eval,options)

P_sub=deepcopy(P_sub_c1);
append!(P_sub,P_sub_c2)
append!(P_sub,P_sub_sum)

options.rho_ini = [1.0f0]
options.rho_ini = ones(length(TD_OP))*options.rho_ini[1]
for i=1:length(options.rho_ini)
  if set_Prop.ncvx[i]==true
    options.rho_ini[i]=0.1f0
  end
end

#compute generalized Minkowski decomposition
@time (v_decomposition,log_PARSDMM,dummy1,dummy2) = PARSDMM(vec(m_evaluation),AtA,TD_OP,set_Prop,P_sub,comp_grid_eval,options);
@time (v_decomposition,log_PARSDMM,dummy1,dummy2) = PARSDMM(vec(m_evaluation),AtA,TD_OP,set_Prop,P_sub,comp_grid_eval,options);

figure();
subplot(3, 3, 3);semilogy(log_PARSDMM.r_pri)          ;title("r primal")
subplot(3, 3, 4);semilogy(log_PARSDMM.r_dual)         ;title("r dual")
subplot(3, 3, 1);semilogy(log_PARSDMM.obj)            ;title(L"$ \frac{1}{2} || \mathbf{m}-\mathbf{x} ||_2^2 $")
subplot(3, 3, 2);semilogy(log_PARSDMM.set_feasibility);title("TD feasibility violation")
subplot(3, 3, 5);plot(log_PARSDMM.cg_it)              ;title("nr. of CG iterations")
subplot(3, 3, 6);semilogy(log_PARSDMM.cg_relres)      ;title("CG rel. res.")
subplot(3, 3, 7);semilogy(log_PARSDMM.rho)            ;title("rho")
subplot(3, 3, 8);plot(log_PARSDMM.gamma)              ;title("gamma")
subplot(3, 3, 9);semilogy(log_PARSDMM.evol_x)         ;title("x evolution")
tight_layout()

#plot
m_est_background = v_decomposition[1:prod(comp_grid_eval.n)]
m_est_anomaly    = v_decomposition[1+prod(comp_grid_eval.n):end]
m_est_sum        = m_est_background + m_est_anomaly

m_est_background = reshape(m_est_background,comp_grid_eval.n)
m_est_anomaly    = reshape(m_est_anomaly,comp_grid_eval.n)
m_est_sum        = reshape(m_est_sum,comp_grid_eval.n)

figure(figsize=(10,9))
for i=1:4
  subplot(4,3,(i-1)*3+1);imshow(m_est_background[:,:,i*25],cmap="gray");axis("off");colorbar();title(string("Background est, ",i*25)) #title("training image", fontsize=10)
  subplot(4,3,(i-1)*3+2);imshow(m_est_anomaly[:,:,i*25],cmap="gray");axis("off");colorbar();title(string("Anomaly est, ",i*25))
  #subplot(4,4,i+8);imshow(m_est_sum[:,:,i*25],cmap="gray");axis("off");colorbar();;title(string("sum, ",i*25))
  subplot(4,3,(i-1)*3+3);imshow(m_evaluation[:,:,i*25],cmap="gray");axis("off");colorbar();;title(string("Original, ",i*25))
end
tight_layout()
savefig("escalator_decomposition.png")

########### Plot results from another method #############
#run example in MatLab: https://github.com/stephenbeckr/fastRPCA/blob/master/demos/demo_escalatorVideo.m
#figures below are the result of "split_SPCP" in the above script
file     = matopen("L.mat")
L = read(file,"L") #low rank component
file     = matopen("S.mat")
S = read(file,"S") #sparse component

L = reshape(L,comp_grid_eval.n[1],comp_grid_eval.n[2],200)
S = reshape(S,comp_grid_eval.n[1],comp_grid_eval.n[2],200)

figure(figsize=(10,9))
for i=1:4
  subplot(4,3,(i-1)*3+1);imshow(L[:,:,i*25],cmap="gray");axis("off");colorbar();title(string("Background sSPCP, ",i*25)) #title("training image", fontsize=10)
  subplot(4,3,(i-1)*3+2);imshow(S[:,:,i*25],cmap="gray");axis("off");colorbar();title(string("Anomaly sSPCP, ",i*25))
  #subplot(4,4,i+8);imshow(m_est_sum[:,:,i*25],cmap="gray");axis("off");colorbar();;title(string("sum, ",i*25))
  subplot(4,3,(i-1)*3+3);imshow(m_evaluation[:,:,i*25],cmap="gray");axis("off");colorbar();;title(string("Original, ",i*25))
end
tight_layout()
savefig("escalator_decomposition_splitSPCP.png")

#Plot results of both methods in one figure
figure(figsize=(10,10))
for i=1:4
  subplot(5,4,i);imshow(m_evaluation[:,:,i*25],cmap="gray");axis("off");;title(string("Original, ",i*25))
  subplot(5,4,i+4);imshow(m_est_background[:,:,i*25],cmap="gray");axis("off");title(string("Background gMinkowski, ",i*25)) #title("training image", fontsize=10)
  subplot(5,4,i+8);imshow(m_est_anomaly[:,:,i*25],cmap="gray");axis("off");title(string("Anomaly gMinkowski, ",i*25))
  subplot(5,4,i+12);imshow(L[:,:,i*25],cmap="gray");axis("off");title(string("Background sSPCP, ",i*25)) #title("training image", fontsize=10)
  subplot(5,4,i+16);imshow(S[:,:,i*25],cmap="gray");axis("off");title(string("Anomaly sSPCP, ",i*25))
end
tight_layout()
savefig("escalator_decomposition_allresults.png")


text:
@doc """ data flow of whole project
""" ->
function dataflow()

    # make data
    if debug
        #get dict sample_label label_sample
        @info("processing label")
        @time readlabel()

        @info("processing training data")
        @time X,lbs = train_data()

        @info("shuffle and split data")
        data = hcat(X, lbs)
        data = shuffle(data)
        @time train,val,test = split_data(data)
        JLD.save(train_fl, "train", train)
        JLD.save(val_fl,   "val",   val)
        JLD.save(test_fl,  "test",  test)
    end
    train = JLD.load(train_fl, "train")
    val = JLD.load(val_fl, "val")
    test = JLD.load(test_fl, "test")
    # tain model
    @info("train data with XGBoost")
    model_xgboost(train,val,test,100_000)

end

@doc """ submission with the result model
""" ->
function submission()

end


text:
obs(env::CartPoleEnv, ::Nothing) = obs(env, zeros(4))
obs(env::CartPoleEnv, (x, x̄, θ, θ̄)) =
    Dict("x" => x, "theta"=>θ)

struct CartPoleDrawParams <: AbstractDrawParams
    screen_height::UInt32
    screen_width::UInt32
    world_width::Float32
    scale::Float32
    carty::UInt32
    pole_width::Float32
    pole_length::Float32
    cart_width::Float32
    cart_height::Float32
end

CartPoleDrawParams() =
    CartPoleDrawParams(
        400,       # screen_height
        600,       # screen_width
        48f-1,     # world_width
        600/48f-1, # scale (screen_width / world_width)
        100,       # carty
        10f0,      # pole_width
        600/48f-1, # pole_length (scale)
        50f0,      # cart_width
        30f0       # cart_height
    )

Ctx(env::CartPoleEnv, mode::Symbol = :human_window) = Ctx(env, Val(mode))

function Ctx(::CartPoleEnv, ::Val{:human_pane})
    draw_params = CartPoleDrawParams()
    viewer = CairoRGBSurface(draw_params.screen_width, draw_params.screen_height)

    CairoCtx(draw_params, viewer)
end

@init @require Gtk="4c0ca9eb-093a-5379-98c5-f87ac0bbbf44" function Ctx(::CartPoleEnv, ::Val{:human_window})
    draw_params = CartPoleDrawParams()
    viewer = CairoRGBSurface(draw_params.screen_width, draw_params.screen_height)

    canvas = @GtkCanvas()
    canvas.backcc = CairoContext(viewer)
    win = GtkWindow(canvas, "CartPole",
            draw_params.screen_width, draw_params.screen_height; resizable=false)
    show(canvas)
    visible(win, false)
    signal_connect(win, "delete-event") do widget, event
        ccall((:gtk_widget_hide_on_delete, Gtk.libgtk), Bool, (Ptr{GObject},), win)
    end

    GtkCtx(draw_params, canvas, win)
end

function Ctx(::CartPoleEnv, ::Val{:rgb})
    draw_params = CartPoleDrawParams()
    viewer = CairoRGBSurface(draw_params.screen_width, draw_params.screen_height)

    RGBCtx(draw_params, viewer)
end

function drawcanvas!(env::CartPoleEnv, viewer::CairoContext, params::CartPoleDrawParams)
    # Background
    set_source_rgb(viewer, 1, 1, 1)
    rectangle(viewer, 0, 0, params.screen_width, params.screen_height)
    fill(viewer)

    # Track
    set_source_rgb(viewer, 0, 0, 0)
    move_to(viewer, 0, params.screen_height - params.carty)
    line_to(viewer, params.screen_width, params.screen_height - params.carty)
    set_line_width(viewer, 0.5)
    stroke(viewer)

    # Cart
    cartx = env.state[1] * params.scale + params.screen_width/2f0 # MIDDLE OF CART
    translation_dist = Pair(cartx, params.screen_height - params.carty - params.cart_height/2f0)
    translate(viewer, translation_dist.first, translation_dist.second) # first = x, second = y
    set_source_rgb(viewer, 0, 0, 0)
    rectangle(viewer, 0, 0, params.cart_width, params.cart_height)
    fill(viewer)

    # Undoing translation
    translate(viewer, -translation_dist.first, -translation_dist.second)

    # Pole
    translation_dist = Pair(cartx + params.cart_width/2f0, params.screen_height - params.carty - params.cart_height/8f0)
    # Translating from origin to perform rotation; first = x, second = y
    translate(viewer, translation_dist.first, translation_dist.second)
    # translate(viewer, params.cart_width/2, params.cart_height/2)
    rotate(viewer, env.state[3])
    set_source_rgb(viewer, 8f-1, 6f-1, 4f-1)
    rectangle(viewer, -params.pole_width/2, 0, params.pole_width, -params.pole_length)
    fill(viewer)

    #Axle
    set_source_rgb(viewer, 5f-1, 5f-1, 8f-1)
    circle(viewer, 0, 0, 5)
    fill(viewer)

    # Undoing translations and rotations
    rotate(viewer, -env.state[3])
    translate(viewer, -translation_dist.first, -translation_dist.second)
end


text:
quarter_date_to_number_out = DSGE.quarter_date_to_number(DSGE.quartertodate("2019-Q1"))
quarter_date_to_number_ans = 2019.00
quarter_number_to_date_out = DSGE.quarter_number_to_date(quarter_date_to_number_ans)
quarter_number_to_date_ans = DSGE.quartertodate("2019-Q1")
dates = [DSGE.quartertodate("2019-Q1"),DSGE.quartertodate("2019-Q2"),DSGE.quartertodate("2019-Q3"),DSGE.quartertodate("2019-Q4")]
date_to_float_out = date_to_float(dates, 2)
date_to_float_ans = [2019.25 2019.5 2019.75 2020.;
          2019.25 2019.5 2019.75 2020.]
m = AnSchorfheide(testing = true)
transform_list = [DSGE.loggrowthtopct_annualized_percapita, DSGE.loggrowthtopct_annualized,
                  DSGE.loggrowthtopct, DSGE.loggrowthtopct_percapita,
                  DSGE.logleveltopct_annualized_percapita, DSGE.logleveltopct_annualized,
                  DSGE.quartertoannual, DSGE.loggrowthtopct, identity]
untrans_ans = ["Q/Q Log Growth Rate", "Q/Q Log Growth Rate", "Q/Q Log Growth Rate",
               "Q/Q Log Growth Rate", "Log Level",  "Log Level", "Percent Q/Q",
               "Q/Q Log Growth Rate", ""]
q4_ans = ["Percent 4Q Growth", "Percent 4Q Growth", "Percent 4Q Growth",
          "Percent 4Q Growth", "Percent 4Q Growth",  "Percent 4Q Growth",
          "Percent Annualized", "Percent 4Q Growth", ""]
other_ans = ["Percent Q/Q Annualized", "Percent Q/Q Annualized", "Percent Q/Q Annualized",
          "Percent Q/Q Annualized", "Percent Q/Q Annualized",  "Percent Q/Q Annualized",
          "Percent Annualized", "Percent Q/Q Annualized", ""]

@testset "Check utility functions for plotting" begin
    @test "Real GDP Growth" == DSGE.describe_series(m, :obs_gdp, :obs)
    @test "Output Growth" == DSGE.describe_series(m, :y_t, :pseudo)
    @test "π_t" == DSGE.describe_series(m, :π_t, :states)
    @test "pi_t" == DSGE.describe_series(m, :π_t, :states, detexify = true)
    @test "z" == DSGE.describe_series(m, :z_sh, :shocks)
    for (i,tf) in enumerate(transform_list)
        m.observable_mappings[:obs_gdp].rev_transform = tf
        m.pseudo_observable_mappings[:y_t].rev_transform = tf
        @test untrans_ans[i] == DSGE.series_ylabel(m, :obs_gdp, :obs, untrans = true)
        @test untrans_ans[i] == DSGE.series_ylabel(m, :y_t, :pseudo, untrans = true)
        @test q4_ans[i] == DSGE.series_ylabel(m, :obs_gdp, :obs, fourquarter = true)
        @test q4_ans[i] == DSGE.series_ylabel(m, :y_t, :pseudo, fourquarter = true)
        @test other_ans[i] == DSGE.series_ylabel(m, :obs_gdp, :obs)
        @test other_ans[i] == DSGE.series_ylabel(m, :y_t, :pseudo)
    end
    m.observable_mappings[:obs_gdp].rev_transform = sum
    @test_throws ErrorException DSGE.series_ylabel(m, :obs_gdp, :obs)
    @test_throws ErrorException DSGE.series_ylabel(m, :obs_gdp, :test)
    @test "Standard Deviations" == DSGE.series_ylabel(m, :test, :stdshocks)
    @test isempty(DSGE.series_ylabel(m, :test, :states))
    @test isempty(DSGE.series_ylabel(m, :test, :shocks))
    @test quarter_date_to_number_out == quarter_date_to_number_ans
    @test quarter_number_to_date_out == quarter_number_to_date_ans
    @test date_to_float_out == date_to_float_ans
end

nothing


text:
include("./CapcacityMarket.jl")
include("./capacity_data.jl")
include("./market_operations.jl")


text:
"""
    generate_image(size, porosity, distribution)

Generate an image of overlapping `n`-dimensional balls with radii
sampled from `distribution`. The image is generated with specified
`size` and `porosity`.
"""
function generate_image(size,
                        porosity     :: AbstractFloat,
                        distribution :: Distribution)
    image     = zeros(Bool, size) |> BitArray
    cur_pores = prod(size)
    pores     = porosity * cur_pores

    indices    = CartesianIndices(image)
    fidx, lidx = first(indices), last(indices)
    uidx       = oneunit(fidx)

    while true
        radius  = rand(distribution)
        iradius = radius |> ceil |> Int
        center  = rand(indices)

        from = max(center - iradius*uidx, fidx)
        to   = min(center + iradius*uidx, lidx)

        old_pores = count(iszero, image[from:to])

        for idx in from:to
            r = sqrt(sum(x^2 for x in Tuple(idx - center)))
            if r <= radius
                image[idx] = 1
            end
        end

        new_pores = count(iszero, image[from:to])
        cur_pores += new_pores - old_pores

        if cur_pores < pores
            return image
        end
    end
end


text:
# abstract model definition
"""
Abstract model.
"""
abstract type Model end

# Model interface: mandatory

"""
    lattice(model)

Returns the lattice of a given model.
"""
lattice(::T) where {T <: Model} = error("lattice() not implemented for $T")

# Optional

"""
    parameters(model)

Collects relevant parametrs of a model into a named tuple.
"""
@inline parameters(::Model) = NamedTuple{}()

# See configurations.jl - compression of configurations
# compress, decompress


text:
Base.size(A::StandardArray) = size(A.data)
Base.getindex(A::StandardArray, I::Int) = A.data[I]
Base.getindex(A::StandardArray{T,N}, I::Vararg{Int,N}) where {T,N} = A.data[I...]
Base.setindex(A::StandardArray, filler, I::Int) = A.data[I] = filler
Base.setindex(A::StandardArray{T,N}, filler, I::Vararg{Int,N}) where {T,N} = A.data[I...] = filler

"""
Input Types, they are an array to be inputed, and hold info about how long to last etc.
"""

"""
Input type constructors
"""
AbstractInput(construct::Type{T}, I::Vararg) where {T<:AbstractInput} = construct(I...)
AbstractInput(construct::Type{T}, A::Array; stages=[0,1,0], input_bool=Bool[0,1,0], reward_bool=Bool[0,1,0], punishment_bool=Bool[0,0,0]) where {T<:AbstractInput} = construct(A, stages=stages, input_bool=input_bool, reward_bool=reward_bool, punishment_bool=punishment_bool)
AbstractInput(construct::Type{T}, nn::Tuple; stages=[0,1,0], input_bool=Bool[0,1,0], reward_bool=Bool[0,1,0], punishment_bool=Bool[0,0,0]) where {T<:AbstractInput} = construct(nn, stages=stages, input_bool=input_bool, reward_bool=reward_bool, punishment_bool=punishment_bool)

SparseInput(A::Array; stages=[0,1,0], input_bool=Bool[0,1,0], reward_bool=Bool[0,1,0], punishment_bool=Bool[0,0,0]) = SparseInput(A, stages, input_bool, reward_bool=reward_bool, punishment_bool=punishment_bool)
SparseInput(A::Array, stages::Vector{Int}, input_bool::Vector{Bool}, reward_bool::Vector{Bool}, punishment_bool::Vector{Bool}) = SparseInput(A, stages, input_bool, reward_bool, punishment_bool, cumsum(stages))
function SparseInput(arraysize::Tuple; density=0.1, filler=300, stages=[0,1,0], input_bool=Bool[0,1,0], reward_bool=Bool[0,1,0], punishment_bool=Bool[0,0,0]) # actual density is density^n. n=ndims

	init = fill(0.0, arraysize)
	fillsize = round.(Int, arraysize .* density)
	init[map(x -> 1:x, fillsize)...] .= filler * ones(fillsize)
	shuffle!(init)
	return SparseInput(init, stages, input_bool, reward_bool, punishment_bool)
end

SparseRandInput(A::Array; stages=[0,1,0], input_bool=Bool[0,1,0], reward_bool=Bool[0,1,0], punishment_bool=Bool[0,0,0]) = SparseRandInput(A, stages, input_bool, reward_bool, punishment_bool)
SparseRandInput(A::Array, stages::Vector{Int}, input_bool::Vector{Bool}, reward_bool::Vector{Bool}, punishment_bool::Vector{Bool}) = SparseRandInput(A, stages, input_bool, reward_bool, punishment_bool, cumsum(stages))
function SparseRandInput(arraysize::Tuple; density=0.1, mean_filler=300, filler_range=50, stages=[0,1,1], input_bool=Bool[0,1,0], reward_bool=Bool[0,1,0], punishment_bool=Bool[0,0,0])

	out = fill(0.0, arraysize)
	fillsize = round.(Int, arraysize .* density)
	out[map(x -> 1:x, fillsize)...] .= mean_filler .* ones(fillsize) .+ filler_range .* randn(fillsize)
	shuffle!(out)
	return SparseRandInput(out, stages, input_bool, reward_bool, punishment_bool)
end

RestInput(A::Array; stages=[0,1,0], input_bool=Bool[0,1,0], reward_bool=Bool[0,1,0], punishment_bool=Bool[0,0,0]) = RestInput(A, stages, input_bool, reward_bool, punishment_bool)
RestInput(A::Array, stages::Vector{Int}, input_bool::Vector{Bool}, reward_bool::Vector{Bool}, punishment_bool::Vector{Bool}) = RestInput(A, stages, input_bool, reward_bool, punishment_bool, cumsum(stages))
RestInput(nn::Tuple; stages=[0,1,0], input_bool=Bool[0,1,0], reward_bool=Bool[0,0,0], punishment_bool=Bool[0,0,0]) = RestInput(zeros(nn), stages, input_bool, reward_bool, punishment_bool) #always no DA

"""
Input Sequence
"""
abstract type AbstractSequence{T} <: StandardArray{T,1} end
struct InputSequence <: AbstractSequence{Inputs}
	data::Vector{Inputs}
	inputdurations::Vector{Int}
	default::RestInput
end
InputSequence(inputarray::Vector{Inputs}, default::RestInput) = InputSequence(inputarray, inputdurations(inputarray), default)
InputSequence(inputarray::Vector{<:Inputs}) = InputSequence(inputarray, inputdurations(inputarray), RestInput((length(inputarray[1]),)))

"""
functions for input types
"""
duration(input::AbstractInput) = sum(input.stages)
duration(inseq::InputSequence) = sum(inseq.inputdurations)
inputdurations(inseq::InputSequence) = [duration(it) for it in inseq] #? does this need to be here??
inputdurations(inseq::Vector{<:Inputs}) = [duration(it) for it in inseq]
time_index(inseq::InputSequence, t) = findfirst(cumsum(inseq.inputdurations) .>= t)

function constructinputsequence(nn, typelist; 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, length(typelist))
	typlen = length(typelist)

	for i in 1:typlen
		out[i] = AbstractInput(typelist[i], nn, stages=stages, input_bool=input_bool, reward_bool=reward_bool, punishment_bool=punishment_bool)
	end

	default = RestInput(nn)
	return InputSequence(out, default)
end

emptyinputsequence() = InputSequence(Vector{Inputs}(undef,0))

"""
returns the index (or stage) that the current time corresponds to
"""
function get_input(inseq::InputSequence, t)

	if t > duration(inseq) || t < 1
		error("no input available for that time")
	end

	nowind = time_index(inseq, t)
	tᵢ = nowind == 1 ? t : t - cumsum(inseq.inputdurations)[nowind-1]

	stage = get_stage(inseq[nowind], tᵢ)
	return inseq[nowind], inseq[nowind].input_bool[stage], inseq[nowind].reward_bool[stage], inseq[nowind].punishment_bool[stage]
end

function get_input(sensory::AbstractInput, t)
	stage = get_stage(sensory, t)
	return sensory, sensory.input_bool[stage], sensory.reward_bool[stage], sensory.punishment_bool[stage]
end

get_stage(input::AbstractInput, t) = findfirst(input.cum_stages .>= t)

"""
wrapper for sensory and reward functions
"""
function inputandreward!(t, input, sensory, τ, da_on, da)

	current_sensory, give_input, give_reward, give_punishment = get_input(sensory,t)
	input .= give_input ? current_sensory.data : sensory.default.data
	reward_BA = give_reward ? da_on : 0.0
	punishment_BA = give_punishment ? da_on : 0.0
	da.reward = update_da(da.reward, reward_BA, τ)
	da.punishment = update_da(da.punishment, punishment_BA, τ)

	return da
end

#input_or_default(sensory, give_input)

# prime exameple of where over-defininng input types is actually a hindrence.
#=function inputandreward!(t, input, sensory::AbstractInput, τ, da_on; da=0)
	give_input, give_da = get_input(sensory, t)
	input .= give_input ? sensory.data : sensory.default.data
	BA = give_da ? da_on : 0
	da = update_da(da, BA, τ)

	return(BA,da)
end
=#


text:
"""
    unicodesymbol2string(us::Char)

Return an ascii representation (LaTeX name for a unicode symbol) for the input character. If the input character is already an ASCII or it is not listed in `REPL.REPLCompletions.latex_symbols`, the input string is returned.

```julia-repl
julia> unicodesymbol2string("β")
"beta"
julia> unicodesymbol2string("b")
"b"
```
"""
function unicodesymbol2string(us::Char)
    us = string(us)
    return isascii(us) ? us : isempty(symbol_latex(us)) ? us : symbol_latex(us)[2:end]
end

"""
    reverse_dict(d::Dict)

Return a Dict instance with each `key => val` pair in `d` reversed as `val => key` pair. Same values in `d` will be overwirten by the last occurd pair.
"""
reverse_dict(d::Dict) = Dict(v => k for (k, v) in d)

"""
    _sort_tuple2(t)

Sort two-element Tuple to increase order. Example

```julia-REPL
julia> _sort_tuple2((2, 1))
(1, 2)
julia> _sort_tuple2((1, 2))
(1, 2)
```
"""
function _sort_tuple2(t)
    a, b = t
    return a > b ? (b, a) : (a, b)
end

"""
    retrieve(dict, key_of_interest, output=[])

Return all values matched the `key_of_interest` in a nested Dict.
"""
function retrieve(dict, key_of_interest, output=[])
    for (key, value) in dict
        if key == key_of_interest
            push!(output, value)
        end
        if value isa AbstractDict
            retrieve(value, key_of_interest, output)
        end
    end
    return output
end


text:
module PostgreSQL

import LibPQ, DBInterface

struct Connection <: DBInterface.Connection
    inner::LibPQ.Connection
end

struct Stmt <: DBInterface.Statement
    inner::LibPQ.Statement
end

DBInterface.connect(::Type{Connection}, args...; kw...) = Connection(LibPQ.Connection(args...; kw...))

DBInterface.close!(conn::Connection) = LibPQ.close(conn.inner)

function DBInterface.prepare(conn::Connection, sql::AbstractString)
    # The definition of the sql is different for `DBInterface` and `libPQ`.

    sql_new = sql
    i = 1
    while true
        r = match(r"(:[\w]+)|([\?])", sql_new)
        if isnothing(r)
            break
        end
        sql_new = sql_new[1:r.offset - 1] * "\$$i" * sql_new[r.offset + length(r.match):end]
        i += 1
    end

    return Stmt(LibPQ.prepare(conn.inner, sql_new))
end

DBInterface.execute(stmt::Stmt, params=()) = LibPQ.execute(stmt.inner, params)

end


text:
__precompile__(true)
module ImageMagick

using Compat
import Compat.String

using FixedPointNumbers, ColorTypes, Images, ColorVectorSpace
using FileIO: DataFormat, @format_str, Stream, File, filename, stream

export MagickWand
export constituteimage
export exportimagepixels!
export getblob
export getimagealphachannel
export getimagecolorspace
export getimagedepth
export getnumberimages
export importimagepixels
export readblob
export readimage
export resetiterator
export setimagecolorspace
export setimagecompression
export setimagecompressionquality
export setimageformat
export writeimage
export image2wand

include("libmagickwand.jl")

typealias AbstractGray{T} Color{T, 1}

const is_little_endian = ENDIAN_BOM == 0x04030201

# Image / Video formats

image_formats = [
    format"BMP",
    format"AVI",
    format"CRW",
    format"CUR",
    format"DCX",
    format"DOT",
    format"EPS",
    format"GIF",
    format"HDR",
    format"ICO",
    format"INFO",
    format"JP2",
    format"JPEG",
    format"PCX",
    format"PDB",
    format"PDF",
    format"PGM",
    format"PNG",
    format"PSD",
    format"RGB",
    format"TIFF",
    format"WMF",
    format"WPG",
    format"TGA"
]

load{T <: DataFormat}(imagefile::File{T}, args...; key_args...) = load_(filename(imagefile), args...; key_args...)
load(filename::AbstractString, args...; key_args...) = load_(filename, args...; key_args...)
save{T <: DataFormat}(imagefile::File{T}, args...; key_args...) = save_(filename(imagefile), args...; key_args...)
save(filename::AbstractString, args...; key_args...) = save_(filename, args...; key_args...)

load{T <: DataFormat}(imgstream::Stream{T}, args...; key_args...) = load_(stream(imgstream), args...; key_args...)
load(imgstream::IO, args...; key_args...) = load_(imgstream, args...; key_args...)
save{T <: DataFormat}(imgstream::Stream{T}, args...; key_args...) = save_(imgstream, args...; key_args...)

const ufixedtype = Dict(10=>UFixed10, 12=>UFixed12, 14=>UFixed14, 16=>UFixed16)

readblob(data::Vector{UInt8}) = load_(data)

function load_(file::Union{AbstractString,IO,Vector{UInt8}}; ImageType=Image, extraprop="", extrapropertynames=false)
    wand = MagickWand()
    readimage(wand, file)
    resetiterator(wand)

    if extrapropertynames
        return(getimageproperties(wand, "*"))
    end

    # Determine what we need to know about the image format
    sz = size(wand)
    n = getnumberimages(wand)
    if n > 1
        sz = tuple(sz..., n)
    end

    imtype      = getimagetype(wand)
    havealpha   = getimagealphachannel(wand)
    cs          = getimagecolorspace(wand)
    if imtype == "GrayscaleType" || imtype == "GrayscaleMatteType"
        cs = "Gray"
    end

    depth = getimagechanneldepth(wand, DefaultChannels)
    if depth <= 8
        T = UFixed8     # always use 8-bit for 8-bit and less
    else
        T = ufixedtype[2*((depth+1)>>1)]  # always use an even # of bits (see issue 242#issuecomment-68845157)
    end

    channelorder = cs
    if havealpha
        if channelorder == "sRGB" || channelorder == "RGB"
            if is_little_endian
                T, channelorder = BGRA{T}, "BGRA"
            else
                T, channelorder = ARGB{T}, "ARGB"
            end
        elseif channelorder == "Gray"
            T, channelorder = GrayA{T}, "IA"
        else
            error("Cannot parse colorspace $channelorder")
        end
    else
        if channelorder == "sRGB" || channelorder == "RGB"
            T, channelorder = RGB{T}, "RGB"
        elseif channelorder == "Gray"
            T, channelorder = Gray{T}, "I"
        else
            error("Cannot parse colorspace $channelorder")
        end
    end
    # Allocate the buffer and get the pixel data
    buf = Array(T, sz...)
    exportimagepixels!(buf, wand, cs, channelorder)

    prop = Dict{Compat.UTF8String, Any}()
    orient = getimageproperty(wand, "exif:Orientation", false)
    if haskey(orientation_dict, orient)
        prop["spatialorder"] = orientation_dict[orient]
    else
        warn("orientation $orient not yet supported")
        prop["spatialorder"] = ["x", "y"]
    end
    n > 1 && (prop["timedim"] = ndims(buf))
    prop["colorspace"] = cs

    if extraprop != ""
        for extra in [extraprop;]
            prop[extra] = getimageproperty(wand,extra)
        end
    end

    ImageType(buf, prop)
end

function save_(filename::AbstractString, img, permute_horizontal=true; mapi = mapinfo(img), quality = nothing)
    wand = image2wand(img, mapi, quality, permute_horizontal)
    writeimage(wand, filename)
end

# This differs from `save_` for files because this is primarily used
# by IJulia, and we want to restrict large images to make display faster.
function save_(s::Stream, img, permute_horizontal=true; mapi = Images.mapinfo_writemime(img), quality = nothing)
    wand = image2wand(img, mapi, quality, permute_horizontal)
    blob = getblob(wand, formatstring(s))
    write(stream(s), blob)
end

function image2wand(img, mapi, quality, permute_horizontal=true)
    local imgw
    try
        imgw = map(mapi, img)
    catch
        warn("Mapping to the storage type failed; perhaps your data had out-of-range values?\nTry `map(Images.Clamp01NaN(img), img)` to clamp values to a valid range.")
        rethrow()
    end
    permute_horizontal && (imgw = permutedims_horizontal(imgw))
    have_color = colordim(imgw)!=0
    if ndims(imgw) > 3+have_color
        error("At most 3 dimensions are supported")
    end
    wand = MagickWand()
    if haskey(img, "colorspace")
        cs = img["colorspace"]
    else
        cs = colorspace(imgw)
    end
    if in(cs, ("RGB", "RGBA", "ARGB", "BGRA", "ABGR"))
        cs = libversion > v"6.7.5" ? "sRGB" : "RGB"
    end
    channelorder = colorspace(imgw)
    if channelorder == "Gray"
        channelorder = "I"
    elseif channelorder == "GrayA"
        channelorder = "IA"
    elseif channelorder == "AGray"
        channelorder = "AI"
    end
    tmp = to_explicit(to_contiguous(data(imgw)))
    constituteimage(tmp, wand, cs, channelorder)
    if quality != nothing
        setimagecompressionquality(wand, quality)
    end
    resetiterator(wand)
    wand
end

formatstring{S}(s::Stream{DataFormat{S}}) = string(S)

# ImageMagick mapinfo client. Converts to RGB and uses UFixed.
mapinfo(img::AbstractArray{Bool}) = MapNone{UFixed8}()
mapinfo{T<:UFixed}(img::AbstractArray{T}) = MapNone{T}()
mapinfo{T<:AbstractFloat}(img::AbstractArray{T}) = MapNone{UFixed8}()
let handled = Set()
for ACV in (Color, AbstractRGB)
    for CV in subtypes(ACV)
        (length(CV.parameters) == 1 && !(CV.abstract)) || continue
        CVnew = CV<:AbstractGray ? Gray : RGB
        @eval mapinfo{T<:UFixed}(img::AbstractArray{$CV{T}}) = MapNone{$CVnew{T}}()
        @eval mapinfo{CV<:$CV}(img::AbstractArray{CV}) = MapNone{$CVnew{UFixed8}}()
        CVnew = CV<:AbstractGray ? Gray : BGR
        AC, CA       = alphacolor(CV), coloralpha(CV)
        if AC in handled
            continue
        end
        push!(handled, AC)
        ACnew, CAnew = alphacolor(CVnew), coloralpha(CVnew)
        @eval begin
            mapinfo{T<:UFixed}(img::AbstractArray{$AC{T}}) = MapNone{$ACnew{T}}()
            mapinfo{P<:$AC}(img::AbstractArray{P}) = MapNone{$ACnew{UFixed8}}()
            mapinfo{T<:UFixed}(img::AbstractArray{$CA{T}}) = MapNone{$CAnew{T}}()
            mapinfo{P<:$CA}(img::AbstractArray{P}) = MapNone{$CAnew{UFixed8}}()
        end
    end
end
end
mapinfo(img::AbstractArray{RGB24}) = MapNone{RGB{UFixed8}}()
mapinfo(img::AbstractArray{ARGB32}) = MapNone{BGRA{UFixed8}}()

# Make the data contiguous in memory, this is necessary for
# imagemagick since it doesn't handle stride.
to_contiguous(A::Array) = A
to_contiguous(A::AbstractArray) = copy(A)
to_contiguous(A::SubArray) = copy(A)

to_explicit(A::Image) = to_explicit(data(A))
to_explicit(A::AbstractArray) = to_explicit(copy(A))

to_explicit{T<:UFixed}(A::Array{T}) = reinterpret(FixedPointNumbers.rawtype(T), A)

to_explicit{T<:UFixed}(A::Array{RGB{T}}) = reinterpret(FixedPointNumbers.rawtype(T), A, tuple(3, size(A)...))
to_explicit{T<:AbstractFloat}(A::Array{RGB{T}}) = to_explicit(map(ClampMinMax(RGB{UFixed8}, zero(RGB{T}), one(RGB{T})), A))
to_explicit{T<:UFixed}(A::Array{Gray{T}}) = reinterpret(FixedPointNumbers.rawtype(T), A, size(A))
to_explicit{T<:AbstractFloat}(A::Array{Gray{T}}) = to_explicit(map(ClampMinMax(Gray{UFixed8}, zero(Gray{T}), one(Gray{T})), A))

to_explicit{T<:UFixed}(A::Array{GrayA{T}}) = reinterpret(FixedPointNumbers.rawtype(T), A)
to_explicit{T<:AbstractFloat}(A::Array{GrayA{T}}) = to_explicit(map(ClampMinMax(GrayA{UFixed8}, zero(GrayA{T}), one(GrayA{T})), A))

to_explicit{T<:UFixed}(A::Array{BGRA{T}}) = reinterpret(FixedPointNumbers.rawtype(T), A, tuple(4, size(A)...))
to_explicit{T<:AbstractFloat}(A::Array{BGRA{T}}) = to_explicit(map(ClampMinMax(BGRA{UFixed8}, zero(BGRA{T}), one(BGRA{T})), A))
to_explicit{T<:UFixed}(A::Array{RGBA{T}}) = reinterpret(FixedPointNumbers.rawtype(T), A, tuple(4, size(A)...))
to_explicit{T<:AbstractFloat}(A::Array{RGBA{T}}) = to_explicit(map(ClampMinMax(RGBA{UFixed8}, zero(RGBA{T}), one(RGBA{T})), A))

# Permute to a color, horizontal, vertical, ... storage order (with time always last)
function permutation_horizontal(img)
    cd = colordim(img)
    td = timedim(img)
    p = spatialpermutation(["x", "y"], img)
    if cd != 0
        p[p .>= cd] += 1
        insert!(p, 1, cd)
    end
    if td != 0
        push!(p, td)
    end
    p
end

permutedims_horizontal(img) = permutedims(img, permutation_horizontal(img))

end # module


text:
#! /usr/bin/julia

# Rosetta Code, Five weekends

isdefined(:Date) || using Dates

const wday = Dates.Fri
const lo = 1900
const hi = 2100
const showres = 5

mfive = recur(Date(lo, 1):Month(1):Date(hi, 12)) do m
    Dates.daysinmonth(m) == 31 && Dates.dayofweek(m) == wday
end

println("Considering the years from ", lo, " to ", hi, ".\n")
println("There are ", length(mfive), " months having 5 3-day weekends.")

println("The first ", showres, " such months are:")
for m in mfive[1:showres]
    println("    ", Dates.monthname(m), " ", Dates.year(m))
end

println("\nThe last ", showres, " such months are:")
for m in mfive[end-showres+1:end]
    println("    ", Dates.monthname(m), " ", Dates.year(m))
end

print("\nThere are ", length(filter(y -> !(y in year(mfive)), lo:hi)))
println(" years that have no such months.")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using aI_Calculations
using Test

@testset "aI_Calculations.jl" begin
    @test Decibels(10) == 20
    @test Decibels(123) == 41.79810222878796
    @test Decibels(10, 100) == (20, 40)
    @test Decibels(10, 1_000, 100) == (20, 60, 40)
    @test RMS(0, 100, 200) == 99.74968671630002
    @test RMS(-200, 100, 10) == 114.0175425099138
end


text:
using Distributed
addprocs(16)
@everywhere include("get_objective_solenoid.jl")
get_Javg_vs_s(1)

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
"""
`DisplacementField <: ImageTransformation`

Description
============

A displacement field is an image-size agnostic specification of a
elastic deformation for images with two spatial dimensions and one
color channel. The term "image-size agnostic" means that once
defined, the deformation can be applied to arbitrarily sized images.

As a programmer one can think of a displacement field as a 2d grid
of small 2d vectors placed on each position of said grid.
The positions themself need not be equally spaced. However, both
(x and y) coordinates for each of these positions must be within
the interval [0, 1]. This way the positions are on a relative scale
and can be streched to any image width and height.
The positions are meant to denote reference points of the input image.
Each position has a vector placed on it to specify where the pixel
located on the reference point of the input image should end up on
the output image (so a translation vector basically).

The grid usually has a lot less of these reference points than the
image has pixels. The two reason for this are performance on the
one hand, and also that grids are intended to be used with arbitrary
sized images.
As a consequence of this all the pixel inbetween reference points
need to be interpolated in some form. On way to do that is using
piecewise affine transformations, which is how a DisplacementMesh
does it. See the documentation of DisplacementMesh for more
information.

Usage
======

    DisplacementField(x, y, delta_X, delta_Y)

    DisplacementField(delta_X, delta_Y)

Fields
=======

- **`x`** : The horizontal component of the position for the
displacement vectors in the field. If not specified then
a uniformly spaced grid will be generated automatically.
Note that this must be a one dimensional array (a vector).

- **`y`** : The vertical component of the position for the
displacement vectors in the field. If not specified then
a uniformly spaced grid will be generated automatically.
Note that this must be a one dimensional array (a vector).

- **`delta_X`** : The horizontal components of each of the
displacement vectors themselves. They essentially specify how
much each of the references point should be translated along the
horizontal axis.
Note that this must be a two dimensional array (a matrix).

- **`delta_Y`** : The vertical components of each of the
displacement vectors themselves. They essentially specify how
much each of the references point should be translated along the
vertical axis.
Note that this must be a two dimensional array (a matrix).

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Details
========

Another way that a displacement field can be thought about is
in terms of a discrete two-dimensional vector field of a domain
D ⊂ U × V, where U = {i ∈ R | 0 ≤ i ≤ 1}, and V = {j ∈ R | 0 ≤ j ≤ 1}.
Note that U and V are distinct sets, and that vectors along the
first the dimension share the second dimension and vice versa.
This way later triangulation is trivial, because the position
points will allways be layed out in terms of rectangles.

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Augmentor
    using TestImages

    # load an example image
    img = testimage("lena")

    # Randomly distort the image using a 10x15 smoothed displacement field
    field = gaussian_displacement(10, 15)
    @assert typeof(field) <: DisplacementField
    img_new = transform(field, img)

see also
=========

`uniform_displacement`, and `gaussian_displacement` for generating
random displacement fields

`RandomDisplacement`, and `SmoothedDisplacement` for
generating new random displacement fields on every use

`ImageTransformation`, `transform`
"""
immutable DisplacementField <: ImageTransformation
    x::Vector{Float64}
    y::Vector{Float64}
    delta_X::Matrix{Float64}
    delta_Y::Matrix{Float64}

    function DisplacementField(x::AbstractVector{Float64}, y::AbstractVector{Float64}, delta_X::Matrix{Float64}, delta_Y::Matrix{Float64})
        @assert length(x) == size(delta_X, 2) == size(delta_Y, 2)
        @assert length(y) == size(delta_X, 1) == size(delta_Y, 1)
        @assert issorted(x) && issorted(y)
        @assert minimum(x) >= 0 && maximum(x) <= 1
        @assert minimum(y) >= 0 && maximum(y) <= 1
        new(convert(Vector{Float64}, x), convert(Vector{Float64},y), delta_X, delta_Y)
    end

    function DisplacementField(delta_X::Matrix{Float64}, delta_Y::Matrix{Float64})
        x = linspace(0, 1, size(delta_X,2))
        y = linspace(0, 1, size(delta_X,1))
        DisplacementField(x, y, delta_X, delta_Y)
    end
end

function Base.show(io::IO, df::DisplacementField)
    print(io, "DisplacementField (width: $(length(df.x)), height: $(length(df.y)))")
end

Base.showcompact(io::IO, df::DisplacementField) = print(io, "Displace with $(length(df.x))x$(length(df.y)) displacement field.")
multiplier(::DisplacementField) = 1 # not true, basically infinity...

@recipe function plot(df::DisplacementField, img::Image)
    @series begin
        seriestype := :image
        img
    end
    df, size(img, "x"), size(img, "y")
end

@recipe function plot(df::DisplacementField, width = 1, height = 2)
    h, w = size(df.delta_X)

    yflip := true
    seriestype := :quiver

    quiver := vec([(Float64(df.delta_X[i,j] * width), Float64(df.delta_Y[i,j] * height)) for i=1:h, j=1:w])

    vec([(df.x[j] * width, df.y[i] * height) for i=1:h, j=1:w])
end

function transform{T<:AbstractImage}(df::DisplacementField, img::T)
    dm = DisplacementMesh(df, img)
    result = _transform(dm, img)::T
    _log_operation!(result, df)::T
end

# uniform displacement

"""
`uniform_displacement(gridwidth, gridheight[, scale=0.2, static_border=true, normalize=true])` → `DisplacementField`

Generates a `DisplacementField` by placing uniformly random
generated displacement vectors on an equally space grid of size
`gridheight × gridwidth`.

Note that the resulting `DisplacementField` will remain static,
which means that it will always perform the same displacement.
If one wants to apply a different random displacement on every use
then `RandomDisplacement` should be used instead.

see `RandomDisplacement` for more information on the parameters.
"""
function uniform_displacement(gridwidth::Int, gridheight::Int, scale::Real, static_border::Bool = true, normalize::Bool = true)
    delta_X = _1d_uniform_displacement(gridwidth, gridheight, Float64(scale), static_border, normalize)
    delta_Y = _1d_uniform_displacement(gridwidth, gridheight, Float64(scale), static_border, normalize)
    DisplacementField(delta_X, delta_Y)
end

function uniform_displacement(gridwidth::Int, gridheight::Int; scale = .2, static_border = true, normalize = true)
    uniform_displacement(gridwidth, gridheight, scale, static_border, normalize)
end

function _1d_uniform_displacement(gridwidth::Int, gridheight::Int, scale::Float64, static_border::Bool, normalize::Bool)
    A = if static_border
        @assert gridwidth > 2 && gridheight > 2
        A_t = rand(gridheight, gridwidth)
        _set_border!(A_t, .5)
    else
        @assert gridwidth > 0 && gridheight > 0
        rand(gridheight, gridwidth)
    end::Matrix{Float64}
    broadcast!(*, A, A, 2.)
    broadcast!(-, A, A, 1.)
    if normalize
        broadcast!(*, A, A, scale / norm(A))
    end
    A
end

# gaussian displacement fields

"""
`gaussian_displacement(gridwidth, gridheight[, scale=0.2, sigma=2, iterations=1, static_border=true, normalize=true])` → `DisplacementField`

Generates a `DisplacementField` by placing uniformly random
generated displacement vectors on an equally space grid of size
`gridheight × gridwidth`. In addition, this vectorfield is smoothed
using a gaussian filter with parameter `sigma`. This will result
in a less chaotic displacement field and be much more similar to
an elastic distortion.

Note that the resulting `DisplacementField` will remain static,
which means that it will always perform the same displacement.
If one wants to apply a different random displacement on every use
then `SmoothedRandomDisplacement` should be used instead.

see `SmoothedRandomDisplacement` for more information on the parameters.
"""
function gaussian_displacement(gridwidth::Int, gridheight::Int, scale::Real, sigma::Vector{Float64}, iterations::Int = 1, static_border::Bool = true, normalize::Bool = true)
    delta_X = _1d_gaussian_displacement(gridwidth, gridheight, Float64(scale), sigma, iterations, static_border, normalize)
    delta_Y = _1d_gaussian_displacement(gridwidth, gridheight, Float64(scale), sigma, iterations, static_border, normalize)
    DisplacementField(delta_X, delta_Y)
end

function gaussian_displacement(gridwidth::Int, gridheight::Int, scale::Real, sigma::Real, iterations::Int = 1, static_border::Bool = true, normalize::Bool = true)
    gaussian_displacement(gridwidth, gridheight, Float64(scale), [Float64(sigma), Float64(sigma)], iterations, static_border, normalize)
end

function gaussian_displacement(gridwidth::Int, gridheight::Int; scale = .2, sigma = 2, iterations = 1, static_border = true, normalize = true)
    gaussian_displacement(gridwidth, gridheight, scale, sigma, iterations, static_border, normalize)
end

function _1d_gaussian_displacement(gridwidth::Int, gridheight::Int, scale::Float64, sigma::Vector{Float64}, iterations::Int, static_border::Bool, normalize::Bool)
    @assert length(sigma) == 2
    @assert all(sigma .> 0)
    @assert iterations > 0
    A = if static_border
        A_t = _1d_uniform_displacement(gridwidth, gridheight, 1., false, false)
        _set_border!(A_t, 0.)
        for iter = 1:iterations
            Images.imfilter_gaussian_no_nans!(A_t, sigma)
            _set_border!(A_t, 0.)
        end
        A_t
    else
        A_t = _1d_uniform_displacement(gridwidth, gridheight, 1., false, false)
        for iter = 1:iterations
            Images.imfilter_gaussian_no_nans!(A_t, sigma)
        end
        A_t
    end::Matrix{Float64}
    if normalize
        broadcast!(*, A, A, scale / norm(A))
    end
    A
end

function _set_border!{T}(A::Matrix{T}, val::T)
    h, w = size(A)
    @inbounds for i = 1:h, j = [1,w]
        A[i,j] = val
    end
    @inbounds for i = [1,h], j = 1:w
        A[i,j] = val
    end
    A
end