text:
function newunit(color=(1,0,0),loc=(0,0,Inf),unitspec::Dict=Dict())
	ir=haskey(unitspec,:ir) ? unitspec[:ir] : 2
	pl=haskey(unitspec,:pl) ? unitspec[:pl] : [2]
	passover=haskey(unitspec,:passover) ? unitspec[:passover] : false
	passoverself=haskey(unitspec,:passoverself) ? unitspec[:passoverself] : true
	inclusive=haskey(unitspec,:inclusive) ? unitspec[:inclusive] : false
	groundlevel=haskey(unitspec,:groundlevel) ? unitspec[:groundlevel] : true
	live=haskey(unitspec,:live!) ? unitspec[:live!] : spreadlife!
	baselife=haskey(unitspec,:baselife) ? unitspec[:baselife] : 6
	harvest=haskey(unitspec,:harvest) ? unitspec[:harvest] : unitharvest #(game,unit)->[0.0,0,0,0,0]
	name=haskey(unitspec,:name) ? unitspec[:name] : ""
	costfun=haskey(unitspec,:costfun) ? unitspec[:costfun] : standard_costfun
	canspawn=haskey(unitspec,:canspawn) ? unitspec[:canspawn] : false
	graphic=haskey(unitspec,:graphic) ? unitspec[:graphic] : []
	extra=haskey(unitspec,:extra) ? unitspec[:extra] : Nothing
	return Unit(color,ir,pl,passover,passoverself,inclusive,loc,groundlevel,live,baselife,harvest,name, costfun,canspawn,false,graphic,extra)
end
function newunit(dic::Dict)
	unitspec=units[dic[:name]]
	return newunit(dic[:color],dic[:loc],unitspec)
end
function newgroup(unit::Unit)
	body=[unit]
	units=[unit]
	spawns=Unit[]
	if unit.canspawn
		push!(spawns,unit)
	end
	group=Group(spawns,body,units,[0.0,0,0,0,0],false)
	#group.lifmap=lifmap(group)
	return group
end

function newgroup(units::Array{Unit})
	spawns=Unit[]
	for u in units
		if u.canspawn
			push!(spawns,u)
		end
	end
	body=Unit[]
	if !isempty(spawns)
		body=getcellgroup(spawns[1])
	end
	group=Group(spawns,body,units,[0.0,0,0,0,0],false)
	#group.lifmap=lifmap(group)
	return group
end
function newgroup(spawns::Array,units::Array) #deprecated
	return Group(spawns,Unit[],units,[0.0,0,0,0,0],false)
end
function newgroup(spawns::Array,body::Array,units::Array)
	return Group(spawns,body,units,[0.0,0,0,0,0],false)
end

function newboard(shells=6,initlocs=[(0,0,2)],sizemod=5,size=30,offsetx=0,offsety=0,bgcolor=(0,0,0),gridcolor=(1/2,1/2,1/2),grid=0,c=@GtkCanvas())
	if grid==0
		grid=makegrid(shells,initlocs)
	end
	board=Board(shells,initlocs,grid,c, sizemod,size,offsetx,offsety,0,0,bgcolor,gridcolor)
	return board
end

function newgame(name=string(round(Integer,time())),boardparams=[];unitparams=["standard"],map=Dict(),unit=0,color=(1,0,0),colors=colorsets[1],colind=1,colmax=3,colock=false,delete=false,sequence=[newunit((1,0,0),(-1,1,2),units["queen"]),newunit((0,1,0),(0,-1,2),units["queen"]),newunit((0,0,1),(1,0,2),units["queen"])],board=0,printscore=false,points=[0.0,0,0,0,0],season=0,win=0,window=(900,700),autoharvest=true)
#sequence=[newunit((1,0,0),(-1,1,2),units["queen"]),newunit((0,1,0),(0,-1,2),units["queen"]),newunit((0,0,1),(1,0,2),units["queen"])] #make arbitrary sequences to start games with: newgame(sequence=yoursequence). For example tunnelracing where storages in the corners have to be connected but only the first person to cross across the board can do it without tunneling
	if board==0
		board=newboard(boardparams...)
	end
	for loc in board.grid
		map[loc]=0
	end
	game=Game(name,map,Group[],Unit[],Unit[],unitparams,color,colors,colind,colmax, colock,delete,sequence,board,printscore,points,season,win,window,0,0,Dict(),autoharvest)
	if win==0
		box=GtkBox(:h)
		savebtn=GtkButton("Save")
		loadbtn=GtkButton("Load")
		nameentry=GtkEntry()
		set_gtk_property!(nameentry,:text,game.name)
		scorelabel=GtkLabel("")#pointslabel(game))
		newslabel=GtkLabel("")#infolabel(game))
		passbtn=GtkButton("Next color")
		backbtn=GtkButton("Prev color")
		#colabel=GtkLabel(string(game.color))
		#autoharvestcheck = GtkCheckButton("Autoharvest")
		#set_gtk_property!(autoharvestcheck,:active,game.autoharvest)
		colockcheck=GtkCheckButton("Lock color")
		set_gtk_property!(colockcheck,:active,game.colock)
		scalemaxfac=100
		zoomscale = GtkScale(false, 1:scalemaxfac*10)
		zlabel=GtkLabel("Zoom")
		zadj=Gtk.Adjustment(zoomscale)
		set_gtk_property!(zadj,:value,game.board.sizemod*10)
		omax=scalemaxfac*30
		xoscale = GtkScale(false, -omax:omax)
		xlabel=GtkLabel("Pan x")
		xadj=Gtk.Adjustment(xoscale)
		set_gtk_property!(xadj,:value,game.board.panx)
		yoscale = GtkScale(false, -omax:omax)
		ylabel=GtkLabel("Pan y")
		yadj=Gtk.Adjustment(yoscale)
		set_gtk_property!(yadj,:value,game.board.pany)
		spexpx=GtkSpinButton(-1000:1000)
		Gtk.G_.value(spexpx,0)
		spexpy=GtkSpinButton(-1000:1000)
		Gtk.G_.value(spexpy,0)
		spexpshell=GtkSpinButton(0:3000)
		Gtk.G_.value(spexpshell,game.board.shells+1)
		xexplabel=GtkLabel("X")
		yexplabel=GtkLabel("Y")
		shellexplabel=GtkLabel("radius")
		placebtn=GtkButton("Place unit at")
		expbtn=GtkButton("Expand board at (X,Y)")
		centerbtn=GtkButton("Center board on (X,Y)")
		deletecheck=GtkCheckButton("Enable deletion")
		set_gtk_property!(deletecheck,:active,game.delete)
		clabel1=GtkLabel("Place")
		clabel2=GtkLabel("units")
		withlabel=GtkLabel("with")
		unitscombo=GtkComboBoxText()
		staind=0;staindset=false
		for c in keys(units)
			push!(unitscombo,c)
			if !staindset && c=="standard"
				staindset=true
			elseif !staindset
				staind+=1
			end
		end
		set_gtk_property!(unitscombo,:active,staind)
		g=GtkGrid()
		g[1,1]=savebtn
		g[2,1]=nameentry
		g[3,1]=loadbtn
		g[1,2]=scorelabel
		g[2,2]=newslabel
		g[2,4]=deletecheck
		g[1,3]=passbtn
		g[1,4]=backbtn
		g[2,3]=colockcheck
		g[1,5]=clabel1
		g[3,5]=clabel2
		g[2,5]=unitscombo
		g[1,6]=zlabel
		g[1,7]=xlabel
		g[1,8]=ylabel
		g[2,6]=zoomscale
		g[2,7]=xoscale
		g[2,8]=yoscale
		g[2,9]=placebtn
		g[1,10]=xexplabel
		g[1,11]=yexplabel
		g[2,10]=spexpx
		g[2,11]=spexpy
		g[1,13]=shellexplabel
		g[2,13]=spexpshell
		g[3,12]=withlabel
		g[2,12]=expbtn
		g[2,14]=centerbtn
		push!(box,game.board.c)
		push!(box,g)
		game.g=g
		game.gui[:scorelabel]=scorelabel
		game.gui[:newslabel]=newslabel
		game.gui[:yoscale]=yoscale
		game.gui[:xoscale]=xoscale
		game.gui[:zadj]=zadj
		game.gui[:xadj]=xadj
		game.gui[:yadj]=yadj
		game.gui[:colockcheck]=colockcheck
		game.gui[:deletecheck]=deletecheck
		set_gtk_property!(box,:expand,game.board.c,true)
		game.win=GtkWindow(box,"Weilianqi $name",window[1],window[2])
		showall(game.win)
		id = signal_connect(savebtn, "clicked") do widget
			game.name=nameentry.text[String]
			save(game)
		end
		id = signal_connect(loadbtn, "clicked") do widget
			@sigatom loadgame(nameentry.text[String])
		end
		id = signal_connect(passbtn, "clicked") do widget
			pass!(game,colockcheck.active[Bool])
		end
		id = signal_connect(backbtn, "clicked") do widget
			pass!(game,colockcheck.active[Bool],true)
		end
		signal_connect(unitscombo, "changed") do widget, others...
			unitname=Gtk.bytestring( GAccessor.active_text(unitscombo) )
			game.unitparams[end]=unitname
		end
		id = signal_connect(zoomscale, "value-changed") do widget
			wval=Gtk.G_.value(widget)
			game.board.sizemod=wval/10+exp(wval/100)
			x=Gtk.G_.value(spexpx)
			y=Gtk.G_.value(spexpy)
			center(game,(x,y))
		end
		id = signal_connect(xoscale, "value-changed") do widget
			game.board.panx=-Gtk.G_.value(widget)*10
			drawboard(game)
		end
		id = signal_connect(yoscale, "value-changed") do widget
			game.board.pany=-Gtk.G_.value(widget)*10
			drawboard(game)
		end
		id = signal_connect(colockcheck, "clicked") do widget
			game.colock=widget.active[Bool]
		end
		id = signal_connect(deletecheck, "clicked") do widget
			game.delete=widget.active[Bool]
		end
		id = signal_connect(expbtn, "clicked") do widget
			x=Gtk.G_.value(spexpx)
			y=Gtk.G_.value(spexpy)
			r=Gtk.G_.value(spexpshell)
			remain=expandboard!(game,Integer(r),[(x,y,2)])
		end
		id = signal_connect(centerbtn, "clicked") do widget
			x=Gtk.G_.value(spexpx)
			y=Gtk.G_.value(spexpy)
			center(game,(x,y))
		end
		id = signal_connect(placebtn, "clicked") do widget
			x=Gtk.G_.value(spexpx)
			y=Gtk.G_.value(spexpy)
			nu=newunit(game.color,(x,y,2),units[game.unitparams[end]])
			placeunit!(game,nu)
			sync!(game)
			drawboard(game)
		end
	end
	@guarded function drawsignal(widget)
		ctx=getgc(widget)
		h=height(widget)
		w=width(widget)
		game.window=(w,h)
		game.board.size=game.window[2]/(game.board.shells*game.board.sizemod)
		drawboard(game,ctx,w,h)
	end
	draw(drawsignal,game.board.c)
	game.board.c.mouse.button1press = @guarded (widget, event) -> begin
		ctx=getgc(widget)
		h=height(game.board.c)
		w=width(game.board.c)
		nu=newunit(game.color,0,units[game.unitparams[end]])
		size=game.board.size
		offx=game.board.offsetx+game.board.panx
		offy=game.board.offsety+game.board.pany
		q,r=pixel_to_hex(event.x-w/2-offx,event.y-h/2-offy,size)
		maindiff=abs(round(q)-q)+abs(round(r)-r)
		qup,rup=pixel_to_hex(event.x-w/2+size*cos(pi/6)-offx,event.y-h/2+sin(pi/6)*size-offy,size)
		updiff=abs(round(qup)-qup)+abs(round(rup)-rup)
		qdown,rdown=pixel_to_hex(event.x-w/2+size*cos(pi/6)-offx,event.y-h/2-sin(pi/6)*size-offy,size)
		downdiff=abs(round(qdown)-qdown)+abs(round(rdown)-rdown)
		if length(nu.pl)==1
			best=(nu.pl[1]+1)%3+1
		else
			best=findmin([maindiff,updiff,downdiff])[2]
		end
		hex=[(round(Int,q),round(Int,r),2),(round(Int,qup),round(Int,rup),3),(round(Int,qdown),round(Int,rdown),1)][best]
		nu.loc=hex
		exists=in(hex,keys(game.map))
		if exists
			if game.delete==true && game.map[hex]!=0 && isa(game.map[hex],Unit)
				removeunit!(game,game.map[hex])
			elseif game.map[hex]==0 && placeable(game,nu)
				placeunit!(game,nu)
				if !game.colock
					game.colind=game.colind%game.colmax+1
					game.color=game.colors[game.colind]
				end
			end
			sync!(game)
			drawboard(game,ctx,w,h)
			reveal(widget)
		end
		sync!(game) #sync twice to correct corruptions of reality, usually works
	end
	placeseq!(game)
	sync!(game)
	show(game.board.c)
	return game
end


text:
module Manipulate

export manipulate!

using RigidBodyTreeInspector
using DrakeVisualizer
using RigidBodyDynamics
using RigidBodyDynamics: Bounds, position_bounds, lower, upper
using InteractBase: slider, Widget, observe, vbox
using DataStructures: OrderedDict

function remove_infs(b::Bounds, default=Float64(π))
    Bounds(isfinite(lower(b)) ? lower(b) : -default,
           isfinite(upper(b)) ? upper(b) : default)
end

slider_range(joint::Joint) = remove_infs.(position_bounds(joint))
function slider_range(joint::Joint{T, <: QuaternionFloating}) where {T}
    defaults = [1., 1, 1, 1, 10, 10, 10]
    remove_infs.(position_bounds(joint), defaults)
end

slider_labels(joint::Joint) = [string("q", i) for i in 1:num_positions(joint)]
slider_labels(joint::Joint{T, <:QuaternionFloating}) where {T} = ["rw", "rx", "ry", "rz", "x", "y", "z"]

function sliders(joint::Joint, values=zeros(num_positions(joint));
                 bounds=slider_range(joint),
                 labels=slider_labels(joint),
                 resolution=0.01, prefix="")
    map(bounds, labels, values) do b, label, value
        slider(lower(b):resolution:upper(b),
               value=value,
               label=string(prefix, label))
    end
end

function combined_observable(joint::Joint, sliders::AbstractVector)
    map(observe.(sliders)...) do args...
        q = vcat(args...)
        normalize_configuration!(q, joint)
        q
    end
end

function widget(joint::Joint{T, <:Fixed}, args...) where T
    Widget{:rbd_joint}()
end

function widget(joint::Joint, initial_value=zeros(num_positions(joint)); prefix=string(joint, '.'))
    s = sliders(joint, initial_value, prefix=prefix)
    keys = Symbol.(slider_labels(joint))
    w = Widget{:rbd_joint}(OrderedDict(zip(keys, s)))
    w.output = combined_observable(joint, s)
    w.layout = x -> vbox(s...)
    w
end

"""
    manipulate!(callback::Function, state::MechanismState)

Create Interact sliders to manipulate the state of the mechanism, and call
callback(state) each time a slider is changed. This mutates the state in-place.
"""
function manipulate!(callback::Function, state::MechanismState)
    joint_list = joints(state.mechanism)
    widgets = widget.(joint_list, configuration.(state, joint_list))
    keys = Symbol.(joint_list)
    w = Widget{:rbd_manipulator}(OrderedDict(zip(keys, widgets)))
    w.output = map(observe.(widgets)...) do signals...
        for i in 1:length(joint_list)
            if num_positions(joint_list[i]) > 0
                set_configuration!(state, joint_list[i], signals[i])
            end
        end
        callback(state)
    end
    w.layout = x -> vbox(widgets...)
    w
end

"""
    manipulate(callback::Function, mechanism::Mechanism)

Create Interact sliders to manipulate the state of the mechanism, and call
callback(state) each time a slider is changed. This constructs a new
MechanismState object to mutate internally.
"""
manipulate(callback::Function, mechanism::Mechanism) =
    manipulate(callback, MechanismState{Float64}(mechanism))

end

ext: jl


text:
"""
    AbstractBitstream

Inherit from this type to create a custom bitstream type.
"""
abstract type AbstractBitstream <: Number end

"""
    push!(s::AbstractBitstream, b)

Push a bit(s) `b` onto bitstream `s`.

Fields:
- `s::AbstractBitstream`: the bitstream object
- `b`: the bit(s) to push onto the stream
"""
function Base.push!(s::AbstractBitstream, bs::Vector)
    for b in bs
        push!(s, b)
    end
end
Base.push!(s::VecOrMat{<:AbstractBitstream}, bs::VecOrMat) = push!.(s, bs)

"""
    pop!(s::AbstractBitstream)

Pop a bit from bitstream `s`.

Fields:
- `s::AbstractBitstream`: the bitstream object
"""
Base.pop!(s::VecOrMat{<:AbstractBitstream}) = pop!.(s)

"""
    observe(s::AbstractBitstream)

Examine the most recent bit added to the stream without removing it.

Fields:
- `s::AbstractBitstream`: the bitstream object
"""
observe(s::VecOrMat{<:AbstractBitstream}) = observe.(s)

"""
    length(s::AbstractBitstream)

Return the number of bits in `s`.
"""
Base.length(s::AbstractBitstream) = length(bits(s))

Base.iterate(s::AbstractBitstream, state...) = iterate(bits(s), state...)

Base.eltype(s::AbstractBitstream) = eltype(bits(s))

Base.getindex(s::AbstractBitstream, i::Integer) = Iterators.take(bits(s), i) |> collect |> last
Base.firstindex(s::AbstractBitstream) = 1
Base.lastindex(s::AbstractBitstream) = length(s)


text:
module DAparalpipTests

using Test
using LightPropagation: fluence_DA_semiinf_CW, fluence_DA_semiinf_TD
using LightPropagation: fluence_DA_slab_CW, fluence_DA_slab_TD
using LightPropagation: fluence_DA_paralpip_CW, fluence_DA_paralpip_TD
using LightPropagation: flux_DA_semiinf_TD, flux_DA_paralpip_TD, flux_DA_slab_TD

# Test parallelepied against semi-inf and slabs

# CW

# test against semi-inf
@test fluence_DA_paralpip_CW(0.1, 10.0, n_ext = 1.0, n_med = 1.0, rd = [24.0, 25.0, 0.0], rs = [25.0,25.0], L = [50.0,50.0,50.0], xs = 10) ≈ 0.024035998288332954
@test fluence_DA_paralpip_CW(0.1, 10.0, n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0,5.0], L = [10.0,10.0,10.0]) ≈  fluence_DA_semiinf_CW(1.0, 0.1, 10.0,  n_med = 1.0, n_ext = 1.0, z = 0.0)
@test fluence_DA_paralpip_CW(0.23, 25.0, n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0,5.0], L = [10.0,10.0,10.0]) ≈  fluence_DA_semiinf_CW(1.0, 0.23, 25.0,  n_med = 1.0, n_ext = 1.0, z = 0.0)
@test fluence_DA_paralpip_CW(0.4, 42.0, n_ext = 1.0, n_med = 1.0, rd = [3.0, 5.0, 0.0], rs = [5.0,5.0], L = [10.0,10.0,10.0]) ≈  fluence_DA_semiinf_CW(2.0, 0.4, 42.0,  n_med = 1.0, n_ext = 1.0, z = 0.0)
@test fluence_DA_paralpip_CW(0.1, 30.0, n_ext = 1.2, n_med = 1.4, rd = [3.5, 5.0, 0.0], rs = [5.0,5.0], L = [10.0,10.0,10.0], xs = 20) ≈  fluence_DA_semiinf_CW(1.5, 0.1, 30.0,  n_med = 1.4, n_ext = 1.2, z = 0.0)

# test against slab
@test fluence_DA_paralpip_CW(0.1, 10.0, n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0,5.0], L = [10.0,10.0,10.0]) ≈  fluence_DA_slab_CW(1.0, 0.1, 10.0, s = 10.0, n_med = 1.0, n_ext = 1.0, z = 0.0)
@test fluence_DA_paralpip_CW(0.1, 10.0, n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0,5.0], L = [10.0,10.0,0.5]) ≈  fluence_DA_slab_CW(1.0, 0.1, 10.0, s = 0.5, n_med = 1.0, n_ext = 1.0, z = 0.0)
@test fluence_DA_paralpip_CW(0.2, 21.0, n_ext = 1.1, n_med = 1.3, rd = [3.8, 5.0, 0.0], rs = [5.0,5.0], L = [10.0,10.0,0.2]) ≈  fluence_DA_slab_CW(1.2, 0.2, 21.0, s = 0.2, n_med = 1.3, n_ext = 1.1, z = 0.0)
@test fluence_DA_paralpip_CW(0.2, 21.0, n_ext = 1.1, n_med = 1.3, rd = [3.8, 5.0, 0.5], rs = [5.0,5.0], L = [10.0,10.0,1.0]) ≈  fluence_DA_slab_CW(1.2, 0.2, 21.0, s = 1.0, n_med = 1.3, n_ext = 1.1, z = 0.5)

# TD
t = 0.1:0.5:8

# compare against semi-inf
@test fluence_DA_paralpip_TD(t, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, rd = [24.0, 25.0, 0.0], rs = [25.0,25.0], L = [50.0,50.0,50.0], xs = 10) ≈  fluence_DA_semiinf_TD(t, 1.0, 0.1, 10.0,  n_med = 1.0, n_ext = 1.0, z = 0.0)
@test fluence_DA_paralpip_TD(t, 0.1, 70.0, n_ext = 1.0, n_med = 1.0, rd = [2.0, 5.0, 0.0], rs = [5.0,5.0], L = [10.0,10.0,10.0]) ≈  fluence_DA_semiinf_TD(t, 3.0, 0.1, 70.0,  n_med = 1.0, n_ext = 1.0, z = 0.0)

# compare against slabs
@test fluence_DA_paralpip_TD(t, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0, 5.0], L = [10.0, 10.0, 10.0], xs = 10) ≈  fluence_DA_slab_TD(t, 1.0, 0.1, 10.0, s = 10.0, n_med = 1.0, n_ext = 1.0, z = 0.0)
@test fluence_DA_paralpip_TD(t, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0, 5.0], L = [10.0, 10.0, 0.5], xs = 10) ≈  fluence_DA_slab_TD(t, 1.0, 0.1, 10.0, s = 0.5, n_med = 1.0, n_ext = 1.0, z = 0.0)

@test fluence_DA_paralpip_TD(t, 0.2, 26.0, n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0, 5.0], L = [10.0, 10.0, 0.2], xs = 10) ≈  fluence_DA_slab_TD(t, 1.0, 0.2, 26.0, s = 0.2, n_med = 1.0, n_ext = 1.0, z = 0.0)
@test fluence_DA_paralpip_TD(t, 0.4, 32.0, n_ext = 1.0, n_med = 1.0, rd = [3.5, 5.0, 0.1], rs = [5.0, 5.0], L = [10.0, 10.0, 0.2], xs = 10) ≈  fluence_DA_slab_TD(t, 1.5, 0.4, 32.0, s = 0.2, n_med = 1.0, n_ext = 1.0, z = 0.1)

# flux against semi-inf (reflectance)
@test flux_DA_paralpip_TD(t, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, rd = [24.0, 25.0, 0.0], rs = [25.0,25.0], L = [50.0,50.0,50.0], xs = 10) ≈  flux_DA_semiinf_TD(t, 1.0, 0.1, 10.0,  n_med = 1.0, n_ext = 1.0)
@test flux_DA_paralpip_TD(t, 0.2, 12.0, n_ext = 1.2, n_med = 1.0, rd = [23.0, 25.0, 0.0], rs = [25.0,25.0], L = [50.0,50.0,50.0], xs = 10) ≈  flux_DA_semiinf_TD(t, 2.0, 0.2, 12.0,  n_med = 1.0, n_ext = 1.2)

# flux against slabs (transmittance)
@test flux_DA_paralpip_TD(t, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, rd = [5.0, 5.0, 1.0], rs = [5.0, 5.0], L = [10.0, 10.0, 1.0], xs = 10) ≈  flux_DA_slab_TD(t, 0.0, 0.1, 10.0, s = 1.0, n_med = 1.0, n_ext = 1.0, z = 1.0)

end # module

max_issues_repo_name: heltonmc/LightPropagation


text:
using Test

import LinearAlgebra.I
import LinearAlgebra.pinv
import LinearAlgebra.lu

include("../../../lib/KTBC.jl")
import .KTBC.CreateKTBC

# 8
K₅, T₅ = CreateKTBC(5)

6 * (inv(T₅) - inv(K₅))

6 * inv(K₅)

# 9
K₄, T₄, B₄, C₄  = CreateKTBC(4)
f = [
  1
  -1
  1
  -1
]
e = [
  1
  1
  1
  1
]

u = pinv(C₄) * f

function roundWithLosts(value, lost=10^(-5))
    round.(Int, lost * value)
end

# left division operator: \ - action its just like inv(a) * b
@test roundWithLosts(C₄ \ e) == roundWithLosts(inv(C₄) * e)
@test roundWithLosts(C₄ \ f) == roundWithLosts(inv(C₄) * f)

# 10
H₃ = [
  2 -1 0
  -1 2 -1
  0 -1 1
]

invH₃ = round.(Int, inv(H₃))

I₃ = convert(Array{Int}, I(3))

J₃ = reverse(I₃; dims=2)

@test J₃ * T₃ * J₃ == H₃

@test invH₃ == J₃ * round.(Int, inv(T₃)) * J₃

U₃ = createU(3)
S₃ = createS(3)

@test U₃ * U₃' == H₃

@test S₃' * S₃ == invH₃

# 11

U₃ = createS(3)

SE = [
 0.0  0.0  1.0
 0.0  1.0  1.0
 1.0  1.0  1.0
]

NW = [
 1.0  1.0  1.0
 1.0  1.0  0.0
 1.0  0.0  0.0
]

L = [
 1.0  0.0  0.0
 1.0  1.0  0.0
 1.0  1.0  1.0
]

@test J₃ * U₃ == SE

@test U₃ * J₃ == NW

@test J₃ * U₃ * J₃ == L

@test SE * NW == [
 1.0  0.0  0.0
 2.0  1.0  0.0
 3.0  2.0  1.0
]

# 12, 13
L, U = lu(C₄)

@test round.(Int, U * 6) == [
  12  -6   0  -6
  0   9  -6  -3
  0   0   8  -8
  0   0   0   0
]

U = 6 \ round.(Int, U * 6)

# Only non-zero rows
U = U[vec(mapslices(row -> any(row .!= 0), U', dims=1)), :]
Uᵀ = U'

UUᵀ = U * Uᵀ

UUᵀ = 36 \ round.(Int, 36 * U * Uᵀ)
U⁺ = Uᵀ * inv(UUᵀ)

L = C₄ * U⁺

# Proof
@test C₄ == round.(Int, L * U)

# Or
Lₐ = C₄ * pinv(U)
@test C₄ == round.(Int, Lₐ * U)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApimachineryPkgApisMetaV1StatusCause <: SwaggerModel
    field::Any # spec type: Union{ Nothing, String } # spec name: field
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason

    function IoK8sApimachineryPkgApisMetaV1StatusCause(;field=nothing, message=nothing, reason=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusCause, Symbol("field"), field)
        setfield!(o, Symbol("field"), field)
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusCause, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusCause, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1StatusCause

const _property_map_IoK8sApimachineryPkgApisMetaV1StatusCause = Dict{Symbol,Symbol}(Symbol("field")=>Symbol("field"), Symbol("message")=>Symbol("message"), Symbol("reason")=>Symbol("reason"))
const _property_types_IoK8sApimachineryPkgApisMetaV1StatusCause = Dict{Symbol,String}(Symbol("field")=>"String", Symbol("message")=>"String", Symbol("reason")=>"String")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1StatusCause }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1StatusCause))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1StatusCause }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1StatusCause[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1StatusCause }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1StatusCause[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1StatusCause)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1StatusCause }, name::Symbol, val)
end


text:
module STBParameters

    #
    # This module models all the parameters (tax rates, benefit levels, ages, etc.) needed to model
    # the Scottish/UK Tax and Benefit System. Also functions to convert everything to weekly amounts.
    # The default values that everything is initialised to are 2019/20 values, which shows how long I've been at this.
    # At the bottom are functions that can read other sets of parameters from files. The 20/21 and later
    # parameters should be loaded that way rather than by changing the defaults, as the 19/20 system is used
    # heavily in the unit tests.
    #
    using Dates
    using Dates: Date, now, TimeType, Year
    using TimeSeries
    using StaticArrays
    using Parameters
    using DataFrames,CSV

    using ScottishTaxBenefitModel
    using .GeneralTaxComponents: RateBands, WEEKS_PER_YEAR, WEEKS_PER_MONTH
    using .Definitions
    using .Utils
    using .TimeSeriesUtils: fy, fy_array
    using .STBIncomes

    # FIXME make this ordered list
    export IncomeTaxSys, NationalInsuranceSys, TaxBenefitSystem, SavingsCredit
    export WorkingTaxCredit, SavingsCredit, IncomeRules, MinimumWage, PersonalAllowances
    export weeklyise!, annualise!, AgeLimits, HoursLimits, LegacyMeansTestedBenefitSystem
    export HousingBenefits, HousingRestrictions, Premia, ChildTaxCredit, LocalTaxes, CouncilTax
    export state_pension_age, reached_state_pension_age, load_file, load_file!
    export BRMA, loadBRMAs, DEFAULT_BRMA_2021
    export AttendanceAllowance, ChildBenefit, DisabilityLivingAllowance
    export CarersAllowance, PersonalIndependencePayment, ContributoryESA
    export WidowsPensions, BereavementSupport, RetirementPension, JobSeekersAllowance
    export NonMeansTestedSys, MaternityAllowance, ChildLimits
    export BenefitCapSys

    const MCA_DATE = Date(1935,4,6) # fixme make this a parameter

    ## TODO Use Unitful to have currency weekly monthly annual counts as annotations
    # using Unitful

    function Default_Fuel_Dict_2020_21(t::Type):: Fuel_Type_Dict
        Fuel_Type_Dict{t}(
            Missing_Fuel_Type=>0.1,
            No_Fuel=>0.1,
            Other=>0.1,
            Dont_know=>0.1,
            Petrol=>0.25, # dunno
            Diesel=>0.37,
            Hybrid_use_a_combination_of_petrol_and_electricity=>0.16,
            Electric=>0.02,
            LPG=>0.02,
            Biofuel_eg_E85_fuel=>0.02 )
    end

    @with_kw mutable struct AttendanceAllowance{RT<:Real}
        abolished :: Bool = false
        higher :: RT = 87.65
        lower :: RT = 58.70
        extra_people :: RT = 0
        candidates = Set{OneIndex}()
        slot :: Incomes = ATTENDANCE_ALLOWANCE
    end

    @with_kw mutable struct ChildBenefit{RT<:Real}
        abolished :: Bool = false
        first_child :: RT = 20.70
        other_children :: RT = 13.70
        high_income_thresh :: RT = 50_000.0
        withdrawal = 1/100
        guardians_allowance :: RT = 17.20
    end

    @with_kw mutable struct DisabilityLivingAllowance{RT<:Real}
        abolished :: Bool = false
        care_high :: RT = 87.65
        care_middle  :: RT = 58.70
        care_low  :: RT = 23.20
        mob_high  :: RT = 61.20
        mob_low  :: RT = 23.20
        extra_people :: RT = 0
        candidates = Set{OneIndex}()
        care_slot :: Incomes = DLA_SELF_CARE
        mob_slot :: Incomes = DLA_MOBILITY
    end

    @with_kw mutable struct CarersAllowance{RT<:Real}
        abolished :: Bool = false
        allowance :: RT = 66.15
        scottish_supplement :: RT = 231.40 # per 6 months
        hours :: Int = 35
        gainful_employment_min :: RT = 123.0
        earnings = [SELF_EMPLOYMENT_INCOME,WAGES]
        deductions = [INCOME_TAX,NATIONAL_INSURANCE]
        extra_people :: RT = 0
        candidates = Set{OneIndex}()
        slot :: Incomes = CARERS_ALLOWANCE
    end

    @with_kw mutable struct PersonalIndependencePayment{RT<:Real}
        abolished :: Bool = false
        dl_standard :: RT = 58.70
        dl_enhanced :: RT = 87.85
        mobility_standard :: RT = 23.20
        mobility_enhanced :: RT = 61.20
        extra_people :: RT = 0 # FIXME this is the same for mob/daily living
        mobility_candidates = Set{OneIndex}()
        dl_candidates = Set{OneIndex}()
        care_slot :: Incomes = PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING
        mob_slot :: Incomes = PERSONAL_INDEPENDENCE_PAYMENT_MOBILITY
    end

    @with_kw mutable struct ContributoryESA{RT<:Real}
        abolished :: Bool = false
        assessment_u25 :: RT = 57.90
        assessment_25p :: RT = 73.10
        main           :: RT = 73.10
        work           :: RT = 29.05
        support        :: RT = 38.55
    end

    @with_kw mutable struct JobSeekersAllowance{RT}
        abolished :: Bool = false
        u25 :: RT = 57.90
        o24 :: RT = 73.10
    end

    #=
    @with_kw IncapacityBenefit{RT<:Real}
        upens_lower :: RT = 84.65
        upens_higher :: RT = 100.20
        upens_adult_dep :: RT = 50.80
        opens_lower :: RT = 107.65
        opens_higher :: RT = 112.25
        opens_adult_dep :: RT = 62.75

    end
    =#

    @with_kw mutable struct RetirementPension{RT}
        abolished :: Bool = false
        new_state_pension :: RT = 168.60
        # pension_start_date = Date( 2016, 04, 06 )
        cat_a     :: RT = 129.20
        cat_b     :: RT = 129.20
        cat_b_survivor :: RT = 77.45
        cat_d     :: RT = 77.45
    end

    @with_kw mutable struct BereavementSupport{RT}
        abolished :: Bool = false
        # higher effectively just means 'with children';
        lump_sum_higher :: RT = 3_500 # convert to weekly
        lump_sum_lower  :: RT = 2_500
        higher :: RT = 350 # monthly
        lower  :: RT = 100 # monthly
        deaths_after = Date( 2017, 04, 06 ) # noway of using this, but, key
    end

    @with_kw mutable struct WidowsPensions{RT}
        abolished :: Bool = false
        industrial_higher :: RT = 129.20
        industrial_lower :: RT = 38.76
        standard_rate :: RT = 119.90
        parent :: RT = 119.90
        ages = collect(54:-1:45)
        age_amounts = Vector{RT}([111.51,103.11,94.72,86.33,77.94,69.54,61.15,52.76,44.36,35.97])
    end

    @with_kw mutable struct MaternityAllowance{RT}
        abolished :: Bool = false
        rate :: RT = 148.68

    end

    @with_kw mutable struct BenefitCapSys{RT<:Real}
        abolished :: Bool = false
        outside_london_single :: RT = 257.69
        outside_london_couple :: RT = 384.62
        # not really needed, but anyway ..
        inside_london_single :: RT = 296.35
        inside_london_couple  :: RT = 442.31
        uc_incomes_limit :: RT = 542.88
    end

    function weeklyise!( bc :: BenefitCapSys; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        bc.uc_incomes_limit /= wpm
    end
    #
    # initial version - will be progressively replaced
    # with actual calculations based on disability, hours caring etc.
    @with_kw mutable struct NonMeansTestedSys{RT<:Real}
        attendance_allowance = AttendanceAllowance{RT}()
        child_benefit = ChildBenefit{RT}()
        dla = DisabilityLivingAllowance{RT}()
        carers = CarersAllowance{RT}()
        pip = PersonalIndependencePayment{RT}()
        esa = ContributoryESA{RT}()
        # ?? TODO, maybe incapacity = IncapacityBenefit{RT}()
        jsa = JobSeekersAllowance{RT}()
        pensions = RetirementPension{RT}()
        bereavement = BereavementSupport{RT}()
        widows_pension = WidowsPensions{RT}()
        # young carer grant
        maternity = MaternityAllowance{RT}()
        smp :: RT = 148.68
        # not modelled SDA,Incapacity which we just wrap into
        # ESA
    end

    function weeklyise!( nmt :: NonMeansTestedSys; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        nmt.bereavement.lump_sum_higher /= wpy
        nmt.bereavement.lump_sum_lower /= wpy
        nmt.bereavement.lower /= wpm
        nmt.bereavement.higher /= wpm
        nmt.carers.scottish_supplement /= (wpy/2) # kinda sorta
        nmt.child_benefit.high_income_thresh /= wpy
        # this is unintuitive, but the weekly amount of CB
        # is withdrawn by 1% of the annual excess of income,
        # and the whole model is weekly, so ...
        nmt.child_benefit.withdrawal *= wpy
    end

    @with_kw mutable struct IncomeTaxSys{RT<:Real}
        abolished :: Bool = false
        non_savings_rates :: RateBands{RT} =  [19.0,20.0,21.0,41.0,46.0]
        non_savings_thresholds :: RateBands{RT} =  [2_049.0, 12_444.0, 30_930.0, 150_000.0]
        non_savings_basic_rate :: Int = 2 # above this counts as higher rate

        savings_rates  :: RateBands{RT} =  [0.0, 20.0, 40.0, 45.0]
        savings_thresholds  :: RateBands{RT} =  [5_000.0, 37_500.0, 150_000.0]
        savings_basic_rate :: Int = 2 # above this counts as higher rate

        dividend_rates :: RateBands{RT} =  [0.0, 7.5,32.5,38.1]
        dividend_thresholds :: RateBands{RT} =  [2_000.0, 37_500.0, 150_000.0]
        dividend_basic_rate :: Int = 2 # above this counts as higher rate

        personal_allowance :: RT          = 12_500.00
        personal_allowance_income_limit :: RT = 100_000.00
        personal_allowance_withdrawal_rate  :: RT= 50.0
        blind_persons_allowance    :: RT  = 2_450.00

        married_couples_allowance   :: RT = 8_915.00
        mca_minimum                 :: RT = 3_450.00
        mca_income_maximum          :: RT = 29_600.00
        mca_credit_rate             :: RT = 10.0
        mca_withdrawal_rate         :: RT= 50.0

        marriage_allowance          :: RT = 1_250.00
        personal_savings_allowance  :: RT = 1_000.00

        # FIXME better to have it straight from
        # the book with charges per CO2 range
        # and the data being an estimate of CO2 per type
        company_car_charge_by_CO2_emissions :: Fuel_Type_Dict{RT} = Default_Fuel_Dict_2020_21(RT)
        fuel_imputation  :: RT = 24_100.00

        #
        # pensions
        #
        pension_contrib_basic_amount :: RT = 3_600.00
        pension_contrib_annual_allowance :: RT = 40_000.00
        pension_contrib_annual_minimum :: RT = 10_000.00
        pension_contrib_threshold_income :: RT = 150_000.00
        pension_contrib_withdrawal_rate :: RT = 50.0

        non_savings_income = NON_SAVINGS_INCOME
        all_taxable = ALL_TAXABLE_INCOME
        savings_income = SAVINGS_INCOME
        dividend_income = DIVIDEND_INCOME

        mca_date = MCA_DATE

    end

    function annualise!( it :: IncomeTaxSys; wpy = WEEKS_PER_YEAR, wpm=WEEKS_PER_MONTH )
        it.non_savings_rates .*= 100.0
        it.savings_rates .*= 100.0
        it.dividend_rates .*= 100.0
        it.personal_allowance_withdrawal_rate *= 100.0
        it.non_savings_thresholds .*= wpy
        it.savings_thresholds .*= wpy
        it.dividend_thresholds .*= wpy
        it.personal_allowance *= wpy
        it.blind_persons_allowance *= wpy
        it.married_couples_allowance *= wpy
        it.mca_minimum *= wpy
        it.marriage_allowance *= wpy
        it.personal_savings_allowance *= wpy
        it.pension_contrib_basic_amount *= wpy

        it.mca_income_maximum       *= wpy
        it.mca_credit_rate             *= 100.0
        it.mca_withdrawal_rate         *= 100.0
        for k in it.company_car_charge_by_CO2_emissions
            it.company_car_charge_by_CO2_emissions[k.first] *= wpy
        end
        it.pension_contrib_basic_amount *= wpy
        it.pension_contrib_annual_allowance *= wpy
        it.pension_contrib_annual_minimum *= wpy
        it.pension_contrib_threshold_income *= wpy
        it.pension_contrib_withdrawal_rate *= 100.0
    end

    function weeklyise!( it :: IncomeTaxSys; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )

        it.non_savings_rates ./= 100.0
        it.savings_rates ./= 100.0
        it.dividend_rates ./= 100.0
        it.personal_allowance_withdrawal_rate /= 100.0
        it.non_savings_thresholds ./= wpy
        it.savings_thresholds ./= wpy
        it.dividend_thresholds ./= wpy
        it.personal_allowance /= wpy
        it.blind_persons_allowance /= wpy
        it.married_couples_allowance /= wpy
        it.mca_minimum /= wpy
        it.marriage_allowance /= wpy
        it.personal_savings_allowance /= wpy
        it.mca_income_maximum       /= wpy
        it.mca_credit_rate             /= 100.0
        it.mca_withdrawal_rate         /= 100.0
        for (k,v) in it.company_car_charge_by_CO2_emissions
            it.company_car_charge_by_CO2_emissions[k] /= wpy
        end
        it.pension_contrib_basic_amount /= wpy
        it.pension_contrib_annual_allowance /= wpy
        it.pension_contrib_annual_minimum /= wpy
        it.pension_contrib_threshold_income /= wpy
        it.pension_contrib_withdrawal_rate /= 100.0
    end

    """
    Note this is *very* approximate as the pension
    age actually increases monthly; see https://en.wikipedia.org/wiki/State_Pension_(United_Kingdom)
    """
    function pension_ages()::TimeArray
        n = 2050-2010+1
        dates = Vector{Date}(undef,n)
        males=zeros(Int,n)
        females=zeros(Int,n)
        i = 0
        #
        # we could actually just use this if-else as the function lookup,
        # except we might want to parameterise this. Note if we do and have different m/f ages we have
        # to change the ModelHousehold count-by-age functions.
        #
        for y in 2010:1:2050
            i += 1
            dates[i] = Date( y, 04, 06 )
            if y < 2012
                males[i] = 65
                females[i] = 60
            elseif y < 2014
                males[i] = 65
                females[i] = 61
            elseif y < 2016
                males[i] = 65
                females[i] = 62
            elseif y < 2018
                males[i] = 65
                females[i] = 63
            elseif y < 2020
                males[i] = 65
                females[i] = 64
            elseif y < 2022
                males[i] = 65
                females[i] = 65
            elseif y < 2024
                males[i] = 66
                females[i] = 66
            elseif y < 2046
                males[i] = 67
                females[i] = 67
            else
                males[i] = 68
                females[i] = 68
            end
        end # loop
        data=(dates=dates,females=females,males=males)
        ts = TimeArray(data,timestamp=:dates)
        ts
    end

    @with_kw mutable struct AgeLimits
        state_pension_ages = pension_ages(); # fixme can't serialise using json3
        savings_credit_to_new_state_pension :: Date = Date( 2016, 04, 06 )
    end

    function state_pension_age( limits :: AgeLimits, sex :: Sex, when :: Integer )::Integer
        y = fy( when )
        # temp temp fixme
        limits.state_pension_ages = pension_ages()
        n = size( limits.state_pension_ages )[1]
        if y < timestamp( limits.state_pension_ages[1] )
            py = sex == Male ? limits.state_pension_ages[1].males : limits.state_pension_ages[1].females
        elseif y > timestamp( limits.state_pension_ages[n] )
            py = sex == Male ? limits.state_pension_ages[n].males : limits.state_pension_ages[n].females
        else
            py = sex == Male ? limits.state_pension_ages[y].males : limits.state_pension_ages[y].females
        end
        return values(py)[1]
    end

    function state_pension_age( limits :: AgeLimits, sex :: Sex, when :: DateTime = now() )::Integer
        return state_pension_age( limits, sex, Dates.year( when ))
    end

    function reached_state_pension_age(
        limits :: AgeLimits,
        age  :: Int,
        sex  :: Sex,
        when :: Integer ) :: Bool

        return Utils.age_then( age, when ) >= state_pension_age( limits, sex, when )
    end

    function reached_state_pension_age(
        limits :: AgeLimits,
        age  :: Int,
        sex  :: Sex,
        when :: DateTime = now()) :: Bool
        return reached_state_pension_age( limits, age, sex, Dates.year( when ))
    end

    function reached_state_pension_age(
        limits :: AgeLimits,
        age  :: Int,
        sex  :: Sex,
        when :: Date ) :: Bool
        return reached_state_pension_age( limits, age, sex, Dates.year( when ))
    end

    @with_kw mutable struct NationalInsuranceSys{RT<:Real}
        abolished :: Bool = false
        primary_class_1_rates :: RateBands{RT} = [0.0, 0.0, 12.0, 2.0 ]
        primary_class_1_bands :: RateBands{RT} = [118.0, 166.0, 962.0, 9999999999999.9]
        secondary_class_1_rates :: RateBands{RT} = [0.0, 13.8, 13.8 ]
        secondary_class_1_bands :: RateBands{RT} = [166.0, 962.0, 99999999999999.9 ]
        state_pension_age :: Int = 66; # fixme move
        class_2_threshold ::RT = 6_365.0;
        class_2_rate ::RT = 3.00;
        class_4_rates :: RateBands{RT} = [0.0, 9.0, 2.0 ]
        class_4_bands :: RateBands{RT} = [8_632.0, 50_000.0, 99999999999999.9 ]
        class_1_income = IncludedItems([WAGES],[PENSION_CONTRIBUTIONS_EMPLOYER])
        class_4_income = [SELF_EMPLOYMENT_INCOME]
            ## some modelling of u21s and u25s in apprentiships here..
        # gross_to_net_lookup = BudgetConstraint(undef,0)
    end

    function weeklyise!( ni :: NationalInsuranceSys; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        ni.primary_class_1_rates ./= 100.0
        ni.secondary_class_1_rates ./= 100.0
        ni.class_2_threshold /= wpy
        ni.class_4_rates ./= 100.0
        ni.class_4_bands ./= wpy
    end

    @with_kw mutable struct PersonalAllowances{ RT<:Real }
        age_18_24 :: RT = 57.90
        age_25_and_over :: RT = 73.10 # FIXME rename this allowance !! Applies to u25 ESA people too (CPAG21/22 p165)
        age_18_and_in_work_activity :: RT = 73.10
        over_pension_age :: RT = 181.00
        lone_parent :: RT = 73.10
        lone_parent_over_pension_age :: RT = 181.00
        couple_both_under_18 :: RT = 87.50 # this isn't quite right; see cpag p336
        couple_both_over_18 :: RT = 114.85
        couple_over_pension_age :: RT = 270.60
        couple_one_over_18_high :: RT = 114.85
        couple_one_over_18_med :: RT = 173.10
        pa_couple_one_over_18_low :: RT = 57.90
        child :: RT = 66.90
        pc_mig_single :: RT = 167.25
        pc_mig_couple :: RT = 255.25
        pc_child :: RT = 53.34 # FIXME also something messy for older 1st children
    end

    @with_kw mutable struct Premia{ RT<:Real }
        family :: RT = 17.45
        family_lone_parent = 22.00 # FIXME this is not used??
        disabled_child :: RT = 64.19
        carer_single :: RT = 36.85
        carer_couple :: RT = 73.70 # FIXME is this used?
        disability_single :: RT = 34.35
        disability_couple :: RT = 48.95
        enhanced_disability_child :: RT = 26.04
        enhanced_disability_single :: RT = 16.80
        enhanced_disability_couple :: RT = 24.10
        severe_disability_single :: RT = 65.85
        severe_disability_couple :: RT = 131.70
        pensioner_is :: RT = 140.40
        disability_premium_qualifying_benefits = [
            DLA_SELF_CARE,
            DLA_MOBILITY,
            PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING,
            PERSONAL_INDEPENDENCE_PAYMENT_MOBILITY,
            ATTENDANCE_ALLOWANCE,
            SEVERE_DISABILITY_ALLOWANCE,
            INCAPACITY_BENEFIT,
            SCOTTISH_DISABILITY_ASSISTANCE_CHILDREN_DAILY_LIVING,
            SCOTTISH_DISABILITY_ASSISTANCE_CHILDREN_MOBILITY,
            SCOTTISH_DISABILITY_ASSISTANCE_OLDER_PEOPLE,
            SCOTTISH_DISABILITY_ASSISTANCE_WORKING_AGE_DAILY_LIVING,
            SCOTTISH_DISABILITY_ASSISTANCE_WORKING_AGE_MOBILITY ]
        enhanced_disability_premium_qualifying_benefits = [
            PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING,
            DLA_SELF_CARE,
            ATTENDANCE_ALLOWANCE,
            SCOTTISH_DISABILITY_ASSISTANCE_WORKING_AGE_DAILY_LIVING,
            SCOTTISH_DISABILITY_ASSISTANCE_CHILDREN_DAILY_LIVING,
            SCOTTISH_DISABILITY_ASSISTANCE_OLDER_PEOPLE
        ]
    end

    @with_kw mutable struct WorkingTaxCredit{ RT<:Real }
        ## PA
        abolished :: Bool = false
        basic :: RT = 1_960.00
        lone_parent :: RT = 2_010.00
        couple  :: RT = 2_010.00
        hours_ge_30 :: RT = 810.00
        disability :: RT = 3_165.00
        severe_disability :: RT = 1_365.00
        age_50_plus  :: RT = 1_365.00 # discontinued 2012 - not modelled
        age_50_plus_30_hrs :: RT = 2_030.00 # discontinued - not modelled
        childcare_max_2_plus_children :: RT  = 300.0 # pw
        childcare_max_1_child :: RT  = 175.0
        childcare_proportion :: RT = 70.0 # pct
        taper :: RT = 41.0
        threshold :: RT = 6_420.0
        non_earnings_minima :: RT = 300.0
        # incomes :: Incomes_Dict = make_all_taxable( RT )
    end

    function weeklyise!( wtc :: WorkingTaxCredit; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        wtc.basic /= wpy
        wtc.lone_parent /= wpy
        wtc.couple /= wpy
        wtc.hours_ge_30 /= wpy
        wtc.disability /= wpy
        wtc.severe_disability /= wpy
        wtc.age_50_plus /= wpy
        wtc.age_50_plus_30_hrs /= wpy
        wtc.childcare_proportion /= 100.0
        wtc.non_earnings_minima /= wpy
        wtc.threshold /= wpy
        wtc.taper /= 100.0
    end

    @with_kw mutable struct ScottishChildPayment{ RT<:Real }
        # the guidance is really ambigious about weeks/months
        # just jam on weeks
        amount :: RT = 10.0
        maximum_age :: Int = 5
        qualifying_benefits = [
            CHILD_TAX_CREDIT,
            UNIVERSAL_CREDIT,
            INCOME_SUPPORT,
            WORKING_TAX_CREDIT,
            PENSION_CREDIT,
            NON_CONTRIB_JOBSEEKERS_ALLOWANCE,
            NON_CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE
        ]
    end

    @with_kw mutable struct ChildTaxCredit{ RT<:Real }
        abolished :: Bool = false
        family :: RT = 545.00
        child  :: RT = 2_780.00
        disability :: RT = 3_355.00
        severe_disability :: RT = 1_360.00
        threshold :: RT = 16_105.00
    end

    function weeklyise!( ctc :: ChildTaxCredit; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        ctc.family /= wpy
        ctc.child /= wpy
        ctc.disability /= wpy
        ctc.severe_disability /= wpy
        ctc.threshold/= wpy
    end

    @with_kw mutable struct HoursLimits
        lower :: Int = 16
        med   :: Int = 24
        higher :: Int = 30
    end

    @with_kw mutable struct IncomeRules{RT<:Real}
        permitted_work :: RT= 131.50
        lone_parent_hb :: RT = 25.00
        high :: RT = 20.0
        low_couple :: RT = 10.0
        low_single :: RT = 5.0
        hb_additional :: RT = 17.10
        childcare_max_1 :: RT = 175.00
        childcare_max_2 :: RT = 300.00
        incomes     = LEGACY_MT_INCOME
        hb_incomes  = LEGACY_HB_INCOME
        pc_incomes  = LEGACY_PC_INCOME
        sc_incomes  = LEGACY_SAVINGS_CREDIT_INCOME
        capital_min :: RT = 6_000.0
        capital_max :: RT = 16_000.0
        pc_capital_min :: RT = 10_000.0
        pc_capital_max :: RT = 99999999999999.9
        pensioner_capital_min :: RT = 10_000.0
        pensioner_capital_max :: RT = 16_000.0

        capital_tariff :: RT = 250 # £1pw per 250
        pensioner_tariff :: RT = 500
    end

    @with_kw mutable struct MinimumWage{RT<:Real}
        abolished :: Bool = false
        ages = [16,18,21,25]
        # 19/20
        wage_per_hour :: Vector{RT} = [4.35, 6.15, 7.70, 8.21]
        # 20/1
        # wage_per_hour :: Vector{RT} = [4.55, 6.45, 8.20, 8.72]
        # apprentice_rate :: RT = 4.15;
        apprentice_rate :: RT = 3.90;
    end

    function get_minimum_wage( mwsys :: MinimumWage, age :: Int ):: Real
        p = 0
        if age < mwsys.ages[1]
            return zero(eltype( mwsys.wage_per_hour))
        end
        for r in mwsys.ages
            p += 1
            if age <= r
                break
            end
        end
        return mwsys.wage_per_hour[p]
    end

    function default_band_ds( RT :: Type ) :: Dict
        return Dict{Symbol,RT}(
            :S12000033  =>  1_377.30,
            :S12000034  =>  1_300.81,
            :S12000041  =>  1_206.54,
            :S12000035  =>  1_367.73,
            :S12000036  =>  1_338.59,
            :S12000005  =>  1_304.63,
            :S12000006  =>  1_222.63,
            :S12000042  =>  1_379.00,
            :S12000008  =>  1_375.35,
            :S12000045  =>  1_308.98,
            :S12000010  =>  1_302.62,
            :S12000011  =>  1_289.96,
            :S12000014  =>  1_225.58,
            :S12000047  =>  1_280.80,
            :S12000049  =>  1_386.00,
            :S12000017  =>  1_332.33,
            :S12000018  =>  1_331.84,
            :S12000019  =>  1_409.00,
            :S12000020  =>  1_322.87,
            :S12000013  =>  1_193.49,
            :S12000021  =>  1_342.69,
            :S12000050  =>  1_221.25,
            :S12000023  =>  1_208.48,
            :S12000048  =>  1_318.00,
            :S12000038  =>  1_315.42,
            :S12000026  =>  1_253.91,
            :S12000027  =>  1_206.33,
            :S12000028  =>  1_344.96,
            :S12000029  =>  1_203.00,
            :S12000030  =>  1_344.28,
            :S12000039  =>  1_293.55,
            :S12000040  =>  1_276.42 )
    end

    function default_ct_ratios(RT)
        return Dict{CT_Band,RT}(
            Band_A=>240/360,
            Band_B=>280/360,
            Band_C=>320/360,
            Band_D=>360/360,
            Band_E=>473/360,
            Band_F=>585/360,
            Band_G=>705/360,
            Band_H=>882/360,
            Band_I=>-1, # wales only
            Household_not_valued_separately => 0.0 ) # see CT note
    end

    @with_kw mutable struct CouncilTax{RT<:Real}
        abolished :: Bool = false
        band_d :: Dict{Symbol,RT} = default_band_ds(RT)
        relativities :: Dict{CT_Band,RT} = default_ct_ratios(RT)
        single_person_discount :: RT = 25.0
        # TODO see CT note on disabled discounts
    end

    @with_kw mutable struct LocalTaxes{RT<:Real}
        ct = CouncilTax{RT}()
        # other possible local taxes go here
    end

    function weeklyise!( lt :: LocalTaxes; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        for (c,v) in lt.ct.band_d
            lt.ct.band_d[c] /= wpy
        end
        lt.ct.single_person_discount /= 100.0
    end

    @with_kw mutable struct SavingsCredit{RT<:Real}
        abolished :: Bool = false
        withdrawal_rate :: RT = 60.0
        threshold_single :: RT = 144.38
        threshold_couple :: RT =229.67
        max_single :: RT = 13.73
        max_couple :: RT = 15.35
        available_till = Date( 2016, 04, 06 )
    end

    function weeklyise!( sc :: SavingsCredit; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        sc.withdrawal_rate /= 100.0
    end

     @with_kw mutable struct HousingBenefits{RT<:Real}
        abolished :: Bool = false
        taper :: RT = 65.0
        passported_bens = DEFAULT_PASSPORTED_BENS
        ndd_deductions :: RateBands{RT} =  [15.60,35.85,49.20,80.55,91.70,100.65]
        ndd_incomes :: RateBands{RT} =  [143.0,209.0,271.0,363.0,451.0,99999999999999.9]
     end

     function weeklyise!( hb :: HousingBenefits; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEARs )
        hb.taper /= 100.0
     end

    @with_kw mutable struct LegacyMeansTestedBenefitSystem{RT<:Real}
        # CPAG 2019/bur.pers[pid].20 p335
        abolished :: Bool = false
        ## FIXME we can't turn off pension credit individually here..
        premia :: Premia = Premia{RT}()
        allowances :: PersonalAllowances = PersonalAllowances{RT}()
        income_rules :: IncomeRules = IncomeRules{RT}()
        # FIXME why do we need a seperate copy of HoursLimits here?
        hours_limits :: HoursLimits = HoursLimits()
        savings_credit :: SavingsCredit = SavingsCredit{RT}()
        working_tax_credit = WorkingTaxCredit{RT}()
        child_tax_credit = ChildTaxCredit{RT}()
        hb = HousingBenefits{RT}()
        ctr = HousingBenefits{RT}( false, 20.0, DEFAULT_PASSPORTED_BENS, RateBands{RT}[], RateBands{RT}[])
    end

    struct BRMA{N,T}
        name :: String
        code :: Symbol
        bedrooms :: SVector{N,T}
        room :: T
    end

    function loadBRMAs( N :: Int, T :: Type, file :: String  ) :: Dict{Symbol,BRMA{N,T}}
        bd = CSV.File( file ) |> DataFrame
        println( bd )
        # FIXME infer N from bd
        dict = Dict{ Symbol, BRMA{ N, T }}()
        for r in eachrow( bd )
            println( r.bname )
            obd = BRMA( String(r.bname), Symbol( r.bcode ), SVector{N,T}([r.bed_1,r.bed_2,r.bed_3,r.bed_4]), T(r.room) )
            dict[ Symbol( r.bcode ) ] = obd
        end
        dict
    end

    const DEFAULT_BRMA_2021 = joinpath(MODEL_DATA_DIR,"local", "lha_rates_scotland_2020_21.csv")

    @with_kw mutable struct HousingRestrictions{RT<:Real}
        abolished :: Bool = false
        # Temp till we figure this stuff out
        maximum_rooms :: Int = 4
        rooms_rent_reduction = SVector{2,RT}(14, 25)
        # FIXME!!! load this somewhere
        brmas = loadBRMAs( 4, RT, DEFAULT_BRMA_2021 )
        single_room_age = 35 # FIXME expand
    end

    function weeklyise!( hr :: HousingRestrictions; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        hr.rooms_rent_reduction /= 100.0
    end

    @with_kw mutable struct ChildLimits
        policy_start  ::  TimeType = Date( 2017, 4, 6 )
        max_children  :: Int = 2
    end

    @with_kw mutable struct UniversalCreditSys{RT<:Real}
        abolished :: Bool = false
        threshold :: RT = 2_500.0 ## NOT USED
        age_18_24 :: RT = 251.77
        age_25_and_over :: RT = 317.82

        couple_both_under_25 :: RT = 395.20
        couple_oldest_25_plus :: RT = 498.89

        first_child  :: RT = 277.08
        subsequent_child :: RT = 231.67
        disabled_child_lower :: RT = 126.11
        disabled_child_higher :: RT = 392.08
        limited_capcacity_for_work_activity:: RT = 336.20
        carer ::  RT = 160.20

        ndd :: RT = 73.89

        childcare_max_2_plus_children :: RT  = 1_108.04 # pm
        childcare_max_1_child :: RT  = 646.35
        childcare_proportion :: RT = 85.0 # pct

        minimum_income_floor_hours :: RT = 35*WEEKS_PER_MONTH # FIXME this jams on the DWP 4.35 weeks per month which we may not want in tests

        work_allowance_w_housing :: RT = 287.0
        work_allowance_no_housing :: RT = 503.0
        other_income = UC_OTHER_INCOME
        earned_income :: IncludedItems = UC_EARNED_INCOME
        capital_min :: RT = 6_000.0
        capital_max :: RT = 16_000.0
        # £1 *per week* ≆ 4.35 pm FIXME make 4.35 wpm?
        capital_tariff :: RT = 250.0/4.35
        taper  :: RT= 63.0
        ctr_taper  :: RT = 20.0 # not really part of UC, I suppose, but still...
    end

    function weeklyise!( uc :: UniversalCreditSys; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        uc.threshold /= wpm
        uc.age_18_24  /= wpm
        uc.age_25_and_over  /= wpm
        uc.capital_tariff *= 4.35 # !!! wpm
        uc.couple_both_under_25  /= wpm
        uc.couple_oldest_25_plus  /= wpm
        uc.minimum_income_floor_hours /= wpm
        uc.first_child   /= wpm
        uc.subsequent_child  /= wpm
        uc.disabled_child_lower  /= wpm
        uc.disabled_child_higher  /= wpm
        uc.limited_capcacity_for_work_activity /= wpm
        uc.carer  /= wpm
        uc.ndd /= wpm
        uc.childcare_max_2_plus_children  /= wpm
        uc.childcare_max_1_child  /= wpm
        uc.childcare_proportion  /= 100.0
        uc.taper /= 100.0
        uc.ctr_taper /= 100.0
        uc.work_allowance_w_housing /= wpm
        uc.work_allowance_no_housing /= wpm

    end
    @with_kw mutable struct UBISys{RT}
        # numbers from "horizon 3; 2019 values"
        abolished :: Bool = true
        adult_amount :: RT = 4_800.0
        child_amount :: RT = 3_000.0
        universal_pension :: RT = 8_780.0
        adult_age :: Int = 17
        retirement_age :: Int = 66
    end

    function weeklyise!( ubi :: UBISys; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        ubi.adult_amount /= wpy
        ubi.child_amount /= wpy
        ubi.universal_pension /= wpy
    end

    @with_kw mutable struct TaxBenefitSystem{RT<:Real}
        name :: String = "Scotland 2919/20"
        it   = IncomeTaxSys{RT}()
        ni   = NationalInsuranceSys{RT}()
        lmt  = LegacyMeansTestedBenefitSystem{RT}()
        uc   = UniversalCreditSys{RT}()
        scottish_child_payment = ScottishChildPayment{RT}()
        age_limits = AgeLimits()
        # just a copy of standard ft/pt hours; mt benefits may have their own copy
        hours_limits :: HoursLimits = HoursLimits()
        child_limits :: ChildLimits = ChildLimits()
        minwage = MinimumWage{RT}()
        hr = HousingRestrictions{RT}() # fixme better name
        loctax = LocalTaxes{RT}() # fixme better name
        nmt_bens = NonMeansTestedSys{RT}()
        bencap = BenefitCapSys{RT}()
        ubi = UBISys{RT}()
    end

    function weeklyise!( lmt :: LegacyMeansTestedBenefitSystem; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        println( "weeklyise lmt wpm = $wpm wpy=$wpy")
        weeklyise!( lmt.working_tax_credit; wpm=wpm, wpy=wpy )
        weeklyise!( lmt.child_tax_credit; wpm=wpm, wpy=wpy )
        weeklyise!( lmt.savings_credit; wpm=wpm, wpy=wpy)
        weeklyise!( lmt.hb; wpm=wpm, wpy=wpy )
        weeklyise!( lmt.ctr; wpm=wpm, wpy=wpy )
    end

    function weeklyise!( tb :: TaxBenefitSystem; wpm=WEEKS_PER_MONTH, wpy=WEEKS_PER_YEAR )
        println( "weeklyise tb wpm = $wpm wpy=$wpy")

        weeklyise!( tb.it; wpm=wpm, wpy=wpy )
        weeklyise!( tb.ni; wpm=wpm, wpy=wpy )
        weeklyise!( tb.lmt; wpm=wpm, wpy=wpy )
        weeklyise!( tb.hr; wpm=wpm, wpy=wpy )
        weeklyise!( tb.loctax; wpm=wpm, wpy=wpy )
        weeklyise!( tb.nmt_bens; wpm=wpm, wpy=wpy )
        weeklyise!( tb.uc; wpm=wpm, wpy=wpy )
        weeklyise!( tb.bencap; wpm=wpm, wpy=wpy )
        weeklyise!( tb.ubi; wpm=wpm, wpy=wpy )
    end

   """
   Load a file `filename` and use it to create a modified version
   of the default parameter system. The file should contain
   entries like `sys.it.personal_allowance=999` but can also contain
   arbitrary code.Note: probably not thread safe: poss global variable?
   The file is executed as julia code but doesn't actually need
   a `.jl` extension.
   """
   function load_file( sysname :: AbstractString, RT :: Type = Float64 ) :: TaxBenefitSystem
        begin
            T = RT
            sys = TaxBenefitSystem{T}()
            global sys
            global T
            include( sysname )
        end
        return sys
   end

   """
   Load a file `filename` and use it to modify the given parameter system. The file should contain
   entries like `sys.it.personal_allowance=999` but can also contain
   arbitrary code.Note: probably not thread safe: poss global variable?
   The file is executed as julia code but doesn't actually need
   a `.jl` extension.
   """
   function load_file!( psys :: TaxBenefitSystem{RT}, sysname :: AbstractString ) where RT
            begin
                # I have no idea why scoping
                # like this works ..
                T = RT
                sys = psys
                global sys
                global T
                include( sysname )
            end
   end

end # module


text:
using MLJ

function asymptotics(estimator, problem_fn, rng, Ns)
    abs_mean_rel_errors = []
    abs_vars = []
    ATE = 0
    for n in Ns
        estimates_at_n = []
        for i in 1:10
            T, W, y, ATE = problem_fn(rng; n=n)
            mach = machine(estimator, T, W, y)
            fit!(mach, verbosity=0)
            push!(estimates_at_n, briefreport(mach)[1].estimate)
        end
        push!(abs_mean_rel_errors, 100mean([abs((x-ATE)/ATE) for x in estimates_at_n]))
        push!(abs_vars, var(estimates_at_n))
    end
    return abs_mean_rel_errors, abs_vars
end

max_stars_repo_name: olivierlabayle/TMLE

max_issues_repo_name: olivierlabayle/TMLE

max_forks_repo_name: olivierlabayle/TMLE


text:
function ekepler(m::T,ecc::T) where {T <: Real}
KEPLER_TOL = sqrt(eps(m))
if m != zero(T)
#real*8 e,e0,eps,m,ms,pi2,f0,f1,f2,f3,d1,d2,d3
# This routine solves Kepler's equation for E as a function of (e,M)
# using the procedure outlined in Murray & Dermott:
  pi2=2pi
  ms=mod(m,pi2)
  d3 =one(T)
  de0=ecc*0.85*sign(ms)
  de1=2*de0
  de2=3*de0
  iter = 0
  ITMAX = 20
#  while abs(d3) > KEPLER_TOL
  while true
    de2 = de1
    de1 = de0
    f3=ecc*cos(de0+ms)
    f2=ecc*sin(de0+ms)
#    f1=1.0-f3
#    f0=de0-f2
#    d1=-f0/f1
#    d2=-f0/(f1+0.5*d1*f2)
#    d3=-f0/(f1+d2*0.5*(f2+d2*f3/3.))
#    de0 += d3
    de0 = (f2-de1*f3)/(1-f3)
    iter += 1
    if iter >= ITMAX || de0 == de1 || de0 == de2
      break
    end
  end
  if iter >= ITMAX && !(T == BigFloat && minimum([abs(de0-de1),abs(de0-de2)]) < eps(one(T)))
    println("iterations in ekepler: ",iter," de0: ",de0," de1-de0: ",de1-de0," de2-de0: ",de2-de0)
  end
  ekep=de0+m
else
  ekep = zero(T)
end
return ekep::typeof(m)
end

function ekepler2(m::T,ecc::T) where {T <: Real}
KEPLER_TOL = sqrt(eps(m))
if m != 0.0
#real*8 e,e0,eps,m,ms,pi2,f0,f1,f2,f3,d1,d2,d3
# This routine solves Kepler's equation for E as a function of (e,M)
# using the procedure outlined in Murray & Dermott:
  pi2=2.0*pi
  ms=mod(m,pi2)
  d3 =one(T)
  e0=ms+ecc*0.85*sin(ms)/abs(sin(ms))
  e1=2*e0
  e2=3*e0
#  while abs(d3) > KEPLER_TOL
  ITMAX = 50
  while true
    e2 = e1
    e1 = e0
    f3=ecc*cos(e0)
    f2=ecc*sin(e0)
    f1=1.0-f3
    f0=e0-ms-f2
    d1=-f0/f1
    d2=-f0/(f1+0.5*d1*f2)
    d3=-f0/(f1+d2*0.5*(f2+d2*f3/3.))
    e0=e0+d3
    if iter >= ITMAX || e0 == e1 || e0 == e2
      break
    end
  end
  if iter >= ITMAX
    println("iterations in ekepler2: ",iter," de0: ",de0," de1-de0: ",de1-de0," de2-de0: ",de2-de0)
  end
  ekep=e0+m-ms
else
  ekep = 0.0
end
return ekep::typeof(m)
end

function kepler(m,ecc)
@assert(ecc >= 0.0)
@assert(ecc <= 1.0)
f=m
if ecc > 0
  ekep=ekepler(m,ecc)
#  println(m-ekep+ecc*sin(ekep))
#  f=2.0*atan(sqrt((1.0+ecc)/(1.0-ecc))*tan(0.5*ekep))
#  f=2.0*atan2(sqrt(1.0+ecc)*sin(0.5*ekep),sqrt(1.0-ecc)*cos(0.5*ekep))
  f=2.0*atan(sqrt(1.0+ecc)*sin(0.5*ekep),sqrt(1.0-ecc)*cos(0.5*ekep))
end
return f
end


text:
import Statistics: mean
using Printf
using DelimitedFiles
using Base.Threads
using MultivariateStats
using Statistics
using Flux
"""
    accuracy(x, y, modell)

calculate euclidean distances between points in y and y' for all samples
"""
function accuracy(x, y, modell, dimensions)
  testmode!(modell)
  resp = modell(x)
  testmode!(modell, false)
  resp = relu.(resp)
  mean_accs = zeros(convert(Int32,size(y,1)/dimensions), size(y,2))
  counter = 1
  for i in 1:dimensions:size(y)[1]
    coordsy = y[i:i+(dimensions-1), :]
    coordsresp = resp[i:i+(dimensions-1), :]
    dist = (coordsy .- coordsresp).^2
    dists = zeros(size(dist,2), 1)
    for j in 1:size(dist)[2]
      dists[j, 1] = sqrt(sum(dist[:,j]))
    end
    mean_accs[counter, :] = dists'[:,1]
    counter += 1
  end
  return mean_accs
end

"""
    print_accuracy(x, y, modell)

calculate and print euclidean distances between points in y and y' for all samples
"""
function print_accuracy(x, y, modell, dimensions)
  testmode!(modell)
  resp = modell(x)
  testmode!(modell, false)
  mean_accs = []
  counter = 1
  for i in 1:dimensions:size(y)[1]
    coordsy = y[i:i+(dimensions-1), :]
    coordsresp = resp[i:i+(dimensions-1), :]
    dist = (coordsy .- coordsresp).^2
    dists = []
    for j in 1:size(dist)[2]
      push!(dists, sqrt(sum(dist[:,j])))
    end
    println("mean distance on point ", counter, " is ", mean(dists))
    push!(mean_accs, dists)
    counter += 1
  end
  return mean_accs
end

"""
    cost_whole_data(x, y)

calculates cost (as defined by function 'cost') over an entire 5d tensor
containing 3d image data.
"""
function cost_whole_data_3D(x, y, cost)
  costs = []
  for i in 1:size(y)[2]
    push!(costs, cost(x[:,:,:,:,i:i], y[:,i:i]))
  end
  return sum(costs)
end

"""
    cost_whole_data_2D(x, y)

calculates cost (as defined by function 'cost') over an entire 4d tensor
containing 2d image data.
"""
function cost_whole_data_2D(x, y, cost)
  costs = []
  for i in 1:4:size(y)[2]
    if i+3>size(y, 2)
      push!(costs, cost(x[:,:,:,i:end], y[:,i:end]))
    else
      push!(costs, cost(x[:,:,:,i:i+3], y[:,i:i+3]))
    end
  end
  return sum(costs)
end

"""
    avg_accuracy_per_point(modelo, x, y, dims)

calculates average, min and max distance between predicted and actual landmarks.
"""
function avg_accuracy_per_point(modelo, x, y, dims)
  devs = []
  for i in 1:size(y)[2]
    if size(devs)[1] == 0
      devs = accuracy(x[:,:,:,i:i], y[:,i:i], modelo, dims)
    else
      devs = hcat(devs, accuracy(x[:,:,:,i:i], y[:,i:i], modelo, dims))
    end
  end
  max_acc = []
  min_acc = []
  avg_accs = []
  for i in 1:size(devs)[1]
    push!(avg_accs, median(devs[i, :]))
    push!(min_acc, minimum(devs[i, :]))
    push!(max_acc, maximum(devs[i, :]))
  end
  return avg_accs, max_acc, min_acc
end

"""
    predict_single(model, x)

Returns the output of 'model' in testing mode given input x and applies the relu
function to each output. Takes a single sample.
"""
function predict_single(model, x)
  testmode!(model)
  resp = model(x)
  resp = relu.(resp)
  testmode!(model, false)
  return resp
end

"""
    predict_set(X, model)

Returns the output of 'model' in testing mode given input x and applies the relu
function to each output. Takes multiple samples.
"""
function predict_set(X, model)
  resps = []
  testmode!(model)
  for i in 1:size(X, 4)
    if i == 1
      resps = model(X[:,:,:,i:i])
    else
      resps = cat(resps, model(X[:,:,:,i:i]), dims=2)
    end
  end
  testmode!(model, false)
  return relu.(resps)
end

"""
    save_vols_to_folder(folder, vols, names)

Save a 4D tensor containing 3D volumes to a folder, each in its own
sub folder,  in .tif format.
"""
function save_vols_to_folder(folder, vols, names)
  for i in 1:length(names)
    dir = joinpath(folder, names[i])
    mkdir(dir)
    for img in 1:size(vols, 3)
      name = joinpath(dir, string(img, ".tif"))
      Images.save(name, Images.colorview(Images.Gray, vols[:,:,img,i]))
    end
  end
end

"""
    array_to_lm_file(output_path, coordinates)

Create an avizo landmark file (output_path) with one
set of landmarks. The coordinates to these landmarks need to be passed
in the form of a one-dimensional array as returned by a network when one
volume is given as the input.
"""
function array_to_lm_file(output_path, coordinates)
  io = open(output_path, "w")
  num_lms = 0
  temp = []
  for i in 1:3:length(coordinates)
    push!(temp, string(@sprintf("%.15e ", coordinates[i]),@sprintf("%.15e ", coordinates[i+1]),
    @sprintf("%.15e ", coordinates[i+2])))
    num_lms += 1
  end
  write(io, string("""# Avizo 3D ASCII 3.0

  define Markers $(num_lms)

  Parameters {
  \tNumSets 1,
  \tContentType "LandmarkSet"
  }

  Markers { float[3] Coordinates } @1

  # Data section follows
  @1
  """))
  close(io)
  io = open(output_path, "a")
  writedlm(io, temp,'\n')
  close(io)
end

max_stars_repo_name: imi-diener/AutoLM

max_issues_repo_name: imi-diener/AutoLM

max_forks_repo_name: imi-diener/AutoLM


text:
# ============================================================================ #
# Stack
# ============================================================================ #

abstract type AbstractStack
end

mutable struct Stack{T} <: AbstractStack
    ptr::Int64 # points to the last valid element
    capacity::Int64
    data::Vector{T}

    function Stack{T}(start_size::Int) where {T}
        ptr      = 0
        capacity = start_size < 0 ? 0 : start_size
        data     = Vector{T}(undef, capacity)
        new(ptr, capacity, data)
    end
end
function Stack{T}() where {T}
    start_size = 0
    return Stack{T}(start_size)
end

function next_index(stack::S, index::Int) where {S <: AbstractStack}
    n = index + 1
    if n > stack.capacity
        n = 1
    end
    return n
end

function prev_index(stack::S, index::Int) where {S <: AbstractStack}
    n = index - 1
    if n < 1
        n = stack.capacity
    end
    return n
end

function is_capacity_reached(
    stack::S,
    index = stack.ptr) where {S <: AbstractStack}
    return stack.capacity == 0 || stack.capacity < index
end

function grow_stack!(stack::Stack{T}) where {T}
    growth_factor   = 2
    new_capa        = convert(Int, ceil(stack.capacity * growth_factor))
    stack.capacity  = stack.capacity == 0 ? 1 : new_capa
    new_data        = similar( stack.data, T, stack.capacity )
    copyto!(new_data, stack.data)
    stack.data      = new_data;
    return stack
end

function Base.push!(stack::Stack{T}, new_element::T) where {T}
    stack.ptr += 1;
    if is_capacity_reached(stack)
        grow_stack!(stack)
    end
    stack.data[stack.ptr] = new_element;
    return stack
end

function Base.pop!(stack::Stack{T})::T where {T}
    if isempty(stack)
        underflow_error("empty stack")
    end
    stack.ptr -= 1;
    return stack.data[stack.ptr+1]
end

function Base.peek(stack::Stack{T})::T where {T}
    return stack.data[stack.ptr]
end

function Base.peek(stack::Stack{T}, n::Int = 0)::T where {T}
    return stack.data[stack.ptr - n]
end

function Base.isempty(stack::S) where {S <: AbstractStack}
    return length(stack) <= 0
end

function Base.length(stack::S) where {S <: AbstractStack}
    return stack.ptr
end

function Base.size(stack::S) where {S <: AbstractStack}
    return (stack.ptr,)
end

function capacity(stack::S) where {S <: AbstractStack}
    return stack.capacity
end

function last_valid_index(stack::S) where {S <: AbstractStack}
    return stack.ptr
end

function Base.getindex(stack::Stack{T}, idx::Int) where {T}
    return stack.data[idx]
end

function Base.getindex(stack::Stack{T}, I::Vararg{Int, N}) where {T, N}
    return stack.data[I]
end

function Base.eltype(::Stack{T}) where {T}
    return T
end

function shrink_to_fit!(stack::Stack{T}) where {T}
	vec = similar(stack.data, stack.ptr)
	copyto!(vec, 1, stack.data, 1, stack.ptr)
	stack.data = vec
	stack.capacity = stack.ptr
	return stack
end

function copy(stack::Stack{T}) where {T}
    new_stack           = Stack{T}()
    new_stack.data      = copy(stack.data)
    new_stack.capacity  = stack.capacity
    new_stack.ptr       = stack.ptr
    return new_stack
end

function merge!(stack::Stack{T}, vec::AbstractArray{T}) where {T}
    L = length(vec)
    L_needed = stack.ptr + L
    if L_needed > stack.capacity
        # allocate new array with the given size
        new_data = Vector{T}(undef, L_needed)
        copyto!(new_data, 1, stack.data, stack.ptr)
        copyto!(new_data, stack.ptr+1, vec, 1, L)
    else
        # we have enough space to fit in the new array
        copyto!(stack.data, stack.ptr+1, vec, 1, L)
    end
    return stack
end

function Base.show(io::IO, stack::Stack{T}) where {T}
    println("Stack{$(T)} with properties:")
    println("\t capacity:\t$(stack.capacity)")
    println("\t ptr:\t\t$(stack.ptr)")
    print("\t data:\t\tVector{$(T)}")
    print( view(stack.data, 1:stack.ptr) )
    print("\n\n")
end

mutable struct MatrixStack{T} <: AbstractStack
    ptr::Int64
    capacity::Int64
    data::Matrix{T} # the columns are the items

    function MatrixStack{T}(vec_size::Int, start_size::Int) where {T}
        ptr = 0
        capacity = start_size < 0 ? 0 : start_size
        data     = Matrix{T}(undef, vec_size, capacity)
        new(ptr, capacity, data)
    end
end
function MatrixStack(vec_size::Int, start_size::Int)
    return MatrixStack{Float64}(vec_size, start_size)
end

function grow_stack!(stack::MatrixStack{T}) where {T}
    growth_factor   = 2
    new_capa        = convert(Int, ceil(stack.capacity * growth_factor))
    stack.capacity  = stack.capacity == 0 ? 1 : new_capa

    new_data        = similar( stack.data, size(stack.data,1), stack.capacity )

    copyto!(new_data, stack.data)
    stack.data      = new_data;
    return stack
end

function Base.push!(stack::MatrixStack{T}, new_element) where {T}
    stack.ptr += 1;
    if is_capacity_reached(stack)
        grow_stack!(stack)
    end
    stack.data[:, stack.ptr] = new_element;
    return stack
end

function Base.pop!(stack::MatrixStack{T}) where {T}
    if isempty(stack)
        underflow_error("empty stack")
    end
    stack.ptr -= 1;
    return stack.data[:, stack.ptr+1]
end

function Base.peek(stack::MatrixStack{T}) where {T}
    v = view(stack.data, 1:size(stack.data,1), stack.ptr)
    return v
end

function Base.peek(stack::MatrixStack{T}, n::Int = 0) where {T}
    v = view(stack.data, 1:size(stack.data,1), stack.ptr - n)
    return v
end

function Base.getindex(stack::MatrixStack{T}, idx::Int) where {T}
    return stack.data[idx]
end

function Base.getindex(stack::MatrixStack{T}, I::Vararg{Int, N}) where {T, N}
    return stack.data[I]
end

function shrink_to_fit!(stack::MatrixStack{T}) where {T}
    mat = similar(stack.data, size(stack.data, 1), stack.ptr)
	copyto!( mat, 1, stack.data, 1, stack.ptr * size(stack.data, 1) )
	stack.data     = mat
	stack.capacity = stack.ptr
	return stack
end

# function Base.show(io::IO, stack::MatrixStack{T}) where {T}
# end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
const FILE = Void # not sure if this works
const time_t = Void # not sure if this works


text:
using InfluxDB, Statistics
using Test

container_id = nothing
try
    # First, build the test server
    success(`docker build -t influxdbjl_test_server influxdb_test_instance`)

    # Next, run it in the background
    global container_id = chomp(String(read(`docker run -d --rm -p 8086:8086 influxdbjl_test_server`)))
    @info("Started testing docker influxdb container, running in container $(container_id)")

    include("querying.jl")
    include("databases.jl")
    include("writing.jl")
    include("escaping.jl")
finally
    if container_id !== nothing
        @info("Stopping docker influxdb container $(container_id)")
        run(`docker stop $(container_id)`)
    end
end


text:
using Random
using Zygote

# Copied from Zygote's test code
function ngradient(f, xs::AbstractArray...)
  grads = zero.(xs)
  for (x, Δ) in zip(xs, grads), i in 1:length(x)
    δ = sqrt(eps())
    tmp = x[i]
    x[i] = tmp - δ/2
    y1 = f(xs...)
    x[i] = tmp + δ/2
    y2 = f(xs...)
    x[i] = tmp
    Δ[i] = (y2-y1)/δ
  end
  return grads
end

function gradcheck(f, xs...)
    res = all(isapprox.(ngradient(f, xs...),
                   gradient(f, xs...), rtol = 1e-3, atol = 1e-3))
    res
end

gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(sin.(f(xs...))), xs...)
gradtest(f, dims...) = gradtest(f, rand.(Float64, dims)...)

function gradtestnoerror(f, xs...)
    try
        gradient(f, xs...)
        return true
    catch e
        throw(e)
        return false
    end
end

Random.seed!(0)

X, Y, N, W, R, B = 2, 2, 3, 8, 2, 4
controut = 16

# Dimension map
d = Dict(
    :kr => (W, R, B),
    :βr => (R, B),
    :kw => (W, 1, B),
    :βw => (1, B),
    :wr => (N, R, B),
    :ww => (N, 1, B),
    :u => (N, B),
    :f => (1, B),
    :g => (1, B),
    :π => (3, R, B),
    :L => (N, N, B),
    :M => (N, W, B),
    :controllerout => (Y, B),
    :readvec => (R*W, B),
    :Wr => (Y, R*W, B),
    :ξ => (DNC.outputsize(R, N, W, X, Y)-Y, B),
    :erase => (W, B),
   )

tr = DNC.inputmappings(d[:ξ][1], R, W)
inputs = DNC.split_ξ(rand(d[:ξ]...), tr)

@testset "Addressing" begin
    @test gradtest(DNC.contentaddress, d[:kr], d[:M], d[:βr])
    @test gradtest(DNC.contentaddress, d[:kw], d[:M], d[:βw])
    @test gradtest(DNC.memoryretention, d[:wr], d[:f])
    @test gradtest(DNC.usage, d[:u], d[:ww], d[:wr], d[:f])
    @test gradtest(DNC.writeweight, d[:ww], d[:u], d[:g], d[:g])
    @test gradtest(DNC.precedenceweight, d[:u], d[:ww])
    @test gradtest(DNC.forwardweight, d[:L], d[:wr])
    @test gradtest(DNC.backwardweight, d[:L], d[:wr])
    @test gradtest(DNC.readweight, d[:wr], d[:wr], d[:wr], d[:π])
end

@testset "Utils" begin
    @test gradtest(x -> DNC.oneplus.(x), [-10.0, -1, 0, 1, 10])
    @test gradtest(DNC.calcoutput, d[:controllerout], d[:readvec], d[:Wr])
    @test_broken gradtest(x->DNC.split_ξ(x, tr).kr, d[:ξ])
end

M = rand(Float32, d[:M]...)
L = rand(Float32, d[:L]...)
wr = rand(Float32, d[:wr]...)
u = rand(Float32, d[:u]...)
ma = DNC.MemoryAccess(d[:ξ][1], N, W, R, B)

@testset "Access" begin
    @test gradtest(DNC.eraseandadd, d[:M], d[:ww], d[:erase], d[:erase])
    @test gradtest(DNC.readmem, d[:M], d[:wr])
    # readweights and writeweights use inputs-tuple as argument, so gradtest won't work.
    @test gradtestnoerror((m, i, u) ->sum(DNC.writeweights(m, i, u)), M, inputs, u)
    @test gradtestnoerror((m, i, l, wr) -> sum(DNC.readweights(m, i, l, wr)), M, inputs, L, wr)
    @test gradtestnoerror(x->sum(ma(x)), rand(Float32, d[:ξ]...))
end

model = Dnc(X, Y, controut, N, W, R, B)

@testset "Computer" begin
    @test gradtestnoerror(x->sum(model(x)), rand(Float32, X, B))
end


text:
function updateBound!(u::Array{Float64,2}, size_total_x, size_total_y, neighbors, comm,
                       me, xs, ys, xe, ye, xcell, ycell, nproc)

    mep1 = me + 1

    #assume, to start with, that this process is not going to receive anything
    rreq = Dict{String, MPI.Request}(
                                "N" => MPI.REQUEST_NULL,
                                "S" => MPI.REQUEST_NULL,
                                "E" => MPI.REQUEST_NULL,
                                "W" => MPI.REQUEST_NULL
                                )
    recv = Dict{String, Array{Float64,1}}()
    ghost_boundaries = Dict{String, Any}(
                                "N" => (xe[mep1]+1, ys[mep1]:ye[mep1]),
                                "S" => (xs[mep1]-1, ys[mep1]:ye[mep1]),
                                "E" => (xs[mep1]:xe[mep1], ye[mep1]+1),
                                "W" => (xs[mep1]:xe[mep1], ys[mep1]-1)
                                )
    is_receiving = Dict{String, Bool}("N" => false, "S" => false, "E" => false, "W" => false)

    #send
    neighbors["N"] >=0 && MPI.Isend(u[xe[mep1], ys[mep1]:ye[mep1]], neighbors["N"], me + 40, comm)
    neighbors["S"] >=0 && MPI.Isend(u[xs[mep1], ys[mep1]:ye[mep1]], neighbors["S"], me + 50, comm)
    neighbors["E"] >=0 && MPI.Isend(u[xs[mep1]:xe[mep1], ye[mep1]], neighbors["E"], me + 60, comm)
    neighbors["W"] >=0 && MPI.Isend(u[xs[mep1]:xe[mep1], ys[mep1]], neighbors["W"], me + 70, comm)

    #receive
    if (neighbors["N"] >= 0)
        recv["N"] = Array{Float64,1}(undef, ycell)
        is_receiving["N"] = true
        rreq["N"] = MPI.Irecv!(recv["N"], neighbors["N"], neighbors["N"] + 50, comm)
    end
    if (neighbors["S"] >= 0)
        recv["S"] = Array{Float64,1}(undef, ycell)
        is_receiving["S"] = true
        rreq["S"] = MPI.Irecv!(recv["S"], neighbors["S"], neighbors["S"] + 40, comm)
    end
    if (neighbors["E"] >= 0)
        recv["E"] = Array{Float64,1}(undef, xcell)
        is_receiving["E"] = true
        rreq["E"] = MPI.Irecv!(recv["E"], neighbors["E"], neighbors["E"] + 70, comm)
    end
    if (neighbors["W"] >= 0)
        recv["W"] = Array{Float64,1}(undef, xcell)
        is_receiving["W"] = true
        rreq["W"] = MPI.Irecv!(recv["W"], neighbors["W"], neighbors["W"] + 60, comm)
    end

    MPI.Waitall!([rreq[k] for k in keys(rreq)])
    for (k, v) in is_receiving
        if v
            u[ghost_boundaries[k][1], ghost_boundaries[k][2]] = recv[k]
        end
    end
end

function computeNext!(u0::Array{Float64,2}, u::Array{Float64,2},
    size_total_x::Int, size_total_y::Int, dt::Float64, hx::Float64, hy::Float64,
    me::Int, xs::Array{Int,1}, ys::Array{Int,1}, xe::Array{Int,1}, ye::Array{Int,1},
    nproc::Int, k0::Float64)

    # The stencil of the explicit operator for the heat equation
    # on a regular rectangular grid using a five point finite difference
    # scheme in space is :
    #
    # |                                    weightx * u[i-1][j]                                    |
    # |                                                                                |
    # | weighty * u[i][j-1]   (diagx * weightx + diagy * weighty) * u[i][j]   weighty * u[i][j+1] |
    # |                                                                                |
    # |                                    weightx * u[i+1][j]                                    |

    mep1 = me + 1

    diagx = Float64(-2.0 + hx * hx / (2 * k0 * dt))
    diagy = Float64(-2.0 + hy * hy / (2 * k0 * dt))
    weightx = Float64(k0 * dt / (hx * hx))
    weighty = Float64(k0 * dt / (hy * hy))

    #Perform an explicit update on the points within the domain.
    #Optimization : inner loop on column index (second index) since
    #Julia is row major
    diff = 0.0
    for i = xs[mep1]:xe[mep1]
        for j = ys[mep1]:ye[mep1]
            u[i,j] = weightx * (u0[i-1,j] + u0[i+1,j] + u0[i,j]*diagx) +
                weighty * (u0[i,j-1] + u0[i,j+1] + u0[i,j]*diagy)
        end
    end

   #Compute the difference into domain for convergence.
   #Update the value u0(i,j).
   #Optimization : inner loop on column index (second index) since
   #Julia is row major
   diff = 0.0
   for j = ys[mep1]:ye[mep1]
     for i = xs[mep1]:xe[mep1]
        ldiff = u0[i,j] - u[i,j]
        diff = diff + ldiff * ldiff
        u0[i,j] = u[i,j]
     end
   end

   return diff
end


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>20.0, :alg=>GLGM06(δ=0.01)),
    "X0"   => Hyperrectangle(low=[zeros(64); fill(0.001, 16); fill(-0.002, 4)],
                            high=[zeros(64); fill(0.0015, 16); fill(-0.0015, 4)]),
    "info" => Dict(
        "name" => "pde",
        "dim" => 84,
        "linear" => true,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using HDF5, JLD
using PyPlot

file = "/Users/deguignet/Documents/Julia/imagerec.jld"

admmst = load(file,"admmst")
x = admmst.x
lastiter = load(file, "algost.lastiter")
nfreq = load(file, "algost.nfreq")

tol1 = load(file, "toolst.tol1")
tol2 = load(file, "toolst.tol2")
tol3 = load(file, "toolst.tol3")
tol4 = load(file, "toolst.tol4")
tol5 = load(file, "toolst.tol5")
mydata = load(file, "skyst.mydata")
sky = load(file, "skyst.sky")

nu = load(file, "psfst.nu")
nu0 = load(file, "psfst.nu0")

d = size(x)[1]
if isodd(size(x)[1]) == true
    d12 = round((size(x)[1])/2)
else d12 = d/2
end

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

include("alpha_rec.jl")

#############################################################################
#############################################################################
######################################   Objet reconstruit
# figure(1)
# clf()
# title("Reconstructed object")
# e = 0
# for z in [1 5 10 15]
#     e += 1
#     subplot(2,2,e)
#     a = nu[z]/1e9
#     axis("off")
#     imshow(x[:,:,z],vmin=0,vmax=0.9)
#     title("v = $a[z] GHz")
# end
# subplots_adjust(bottom=0.1, right=0.8, top=0.95)
# cax = axes([0.85, 0.1, 0.025, 0.8])
# colorbar(cax=cax)
#############################################################################
#############################################################################
######################################  RMSE
# figure(2)
# clf()
# errr = zeros(Float64,15)
# for z = 1:nfreq
#     errr[z] = sqrt(sum((x[:,:,z] - sky[:,:,z]).^2)/sum(sky[:,:,z].^2))
# end
# # plot(nu./1e9,err[350,:]',marker="o",linewidth=2)
# plot(nu./1e9,errr[:],marker="o",linewidth=2)
# xlabel(L"\nu \; (GHz)",fontsize=18)
# ylabel(L"Relative \, rmse",fontsize=16)
# axis([1, 1.8, 0.1, 0.2])
#############################################################################
#############################################################################
######################################  cste / alpha / beta
# figure(3)
# clf()
# title("Cste Reconstruction")
# subplot(121)
# colorbar(imshow(-alpharec[:,:,1]))
# subplot(122)
# colorbar(imshow(-alpharecsky[:,:,1]))
#
# figure(4)
# clf()
# title("Alpha Reconstruction")
# subplot(121)
# colorbar(imshow(-alpharec[:,:,2]))
# subplot(122)
# colorbar(imshow(-alpharecsky[:,:,2]))
#
# figure(5)
# clf()
# title("Beta Reconstruction")
# subplot(121)
# colorbar(imshow(-alpharec[:,:,3]))
# subplot(122)
# colorbar(imshow(-alpharecsky[:,:,3]))
############################################# plot alpha
#############################################
figure(6)
clf()
m = d12-40
n = d12+40
p = linspace(m,n,int64(n-m+1))
q = alpharecsky[d12,m:n,2]'
r = alpharec[d12,m:n,2]'
dirty = alpharecdata[d12,m:n,2]'
plot(p,q,linewidth=2)
plot(p,r,linewidth=2)
plot(p,dirty,linewidth=2)
legend( ("Sky","Reconstructed sky", "Dirty image"), loc="upper left")
xlabel(L"pixel",fontsize=18)
ylabel(L"Spectral \, index \, \alpha",fontsize=18)
############################################# plot beta
#############################################
figure(7)
clf()
p = linspace(d12-40,d12+40,81)
q = alpharecsky[d12,d12-40:d12+40,3]'
r = alpharec[d12,d12-40:d12+40,3]'
dirty = alpharecdata[d12,d12-40:d12+40,3]'
plot(p,q)
plot(p,r)
plot(p,dirty)
legend( ("Sky","Reconstructed sky", "Dirty image"), loc="lower center")
xlabel("pix")
ylabel(L"Spectral index $\beta$")
############################################# plot spectre
#############################################
figure(8)
clf()

subplot(1,2,1)
    p = nu./nu0
    i = d12
    j = d12-20

    q = 1*(squeeze(squeeze(x[i,j,:],1),1))
    d = (squeeze(squeeze(sky[i,j,:],1),1))

    plot(log10(p),log10(d),linewidth=2)
    plot(log10(p),log10(q),linewidth=2,linewidth=2)

    label1 = "Sky spectrum"
    label2 = "Reconstructed spectrum"

    legend( (label1, label2), loc="upper right")
    # xlabel(L"$log(\nu/\nu0)$,fontsize=16)
    # xlabel(L"$log(\nu/\nu0)$",fontsize=16)

subplot(1,2,2)
    p = nu./nu0
    i = d12
    j = d12+20

    q = 1.*(squeeze(squeeze(x[i,j,:],1),1))
    d = (squeeze(squeeze(sky[i,j,:],1),1))

    plot(log10(p),log10(d),linewidth=2)
    plot(log10(p),log10(q),linewidth=2)

    # plot(log10(p),log10(s))

    label1 = "Sky spectrum"
    label2 = "Reconstructed spectrum"

    legend( (label1, label2), loc="upper left")
    # xlabel(L"$log(\nu/\nu0)$,fontsize=16)
    # xlabel(L"$log(\nu/\nu0)$",fontsize=16)
#############################################
#############################################

#############################################################################
#############################################################################
######################################  sky / dirty / rec_obj
# figure(9)
# clf()
# title("multiplot")
# e = 1
# for z in [1 5 10 15]
#     subplot(3,4,e)
#     axis("off")
#     a = nu[z]/1e9
#     imshow(sky[:,:,z],vmin=0,vmax=1.3)
#     title(string(L"$\nu \, = \,$","$a ", L"\, GHz"),fontsize=14)
#     e += 1
# end
#
# e = 5
# for z in [1 5 10 15]
#     subplot(3,4,e)
#     axis("off")
#     imshow(mydata[:,:,z]./maximum(mydata).*1.2,vmin=0,vmax=1.3)
#
#     e += 1
# end
# e = 9
# for z in [1 5 10 15]
#     subplot(3,4,e)
#     axis("off")
#     imshow(x[:,:,z],vmin=0,vmax=1.3)
#
#     e += 1
# end
#
# subplots_adjust(bottom=0.1, right=0.8, top=0.9)
# cax = axes([0.85, 0.1, 0.025, 0.8])
# colorbar(cax=cax)
#############################################################################
#############################################################################
######################################  Reconstruction error
# figure(10)
# clf()
# e = 1
# for z in [1 5 10 15]
#     subplot(2,2,e)
#     axis("off")
#     a = nu[z]/1e9
#     # imshow(sqrt(abs(sky[:,:,z]-x[:,:,z])),vmin = 0, vmax = round(maximum(sqrt(abs(sky[:,:,:]-x[:,:,:]))),1))
#     imshow(((sky[:,:,z]-x[:,:,z])),vmin = 0, vmax = 0.1)
#     title(string(L"$\nu \, = \,$","$a ", L"\, GHz"))
#     e += 1
# end
#
# subplots_adjust(bottom=0.1, right=0.8, top=0.9)
# cax = axes([0.85, 0.1, 0.025, 0.8])
# colorbar(cax=cax)
#############################################################################
#############################################################################
######################################  sky / dirty / rec_obj / error
# figure(11)
# clf()
# title("multiplot")
# e = 1
# for z in [1 5 10 15]
#     subplot(4,4,e)
#     axis("off")
#     a = nu[z]/1e9
#     imshow(sky[:,:,z],vmin=0,vmax=1.3)
#     title(string(L"$\nu \, = \,$","$a ", L"\, GHz"),fontsize=14)
#     e += 1
# end
# e = 5
# for z in [1 5 10 15]
#     subplot(4,4,e)
#     axis("off")
#     imshow(mydata[:,:,z]./maximum(mydata).*1.2,vmin=0,vmax=1.3)
#     e += 1
# end
# e = 9
# for z in [1 5 10 15]
#     subplot(4,4,e)
#     axis("off")
#     imshow(x[:,:,z],vmin=0,vmax=1.3)
#     e += 1
# end
#
# subplots_adjust(bottom=0.1, right=0.8, top=0.9)
# cax = axes([0.85, 0.3, 0.025, 0.6])
# colorbar(cax=cax)
#
#
# e = 13
# for z in [1 5 10 15]
#     subplot(4,4,e)
#     axis("off")
#     imshow(((sky[:,:,z]-x[:,:,z])),vmin = 0, vmax = 0.08)
#     e += 1
# end
#
# subplots_adjust(bottom=0.1, right=0.8, top=0.9)
# cax = axes([0.85, 0.1, 0.025, 0.17])
# colorbar(cax=cax)
# #############################################################################
# #############################################################################
# ######################################  sky / dirty / error
# figure(12)
# clf()
# title("multiplot")
# e = 1
# for z in [1 5 10 15]
#     subplot(3,4,e)
#     axis("off")
#     a = nu[z]/1e9
#     imshow(sky[:,:,z],vmin=0,vmax=1.3)
#     title(string(L"$\nu \, = \,$","$a ", L"\, GHz"),fontsize=14)
#     e += 1
# end
# e = 5
# for z in [1 5 10 15]
#     subplot(3,4,e)
#     axis("off")
#     imshow(mydata[:,:,z]./maximum(mydata).*1.2,vmin=0,vmax=1.3)
#     e += 1
# end
#
# subplots_adjust(bottom=0.1, right=0.8, top=0.9)
# cax = axes([0.85, 0.41, 0.025, 0.46])
# colorbar(cax=cax)
#
# e = 9
# for z in [1 5 10 15]
#     subplot(3,4,e)
#     axis("off")
#     imshow(sqrt(abs(sky[:,:,z]-x[:,:,z])),vmin = 0, vmax = round(maximum(sqrt(abs(sky[:,:,:]-x[:,:,:]))),1))
#     e += 1
# end
# #vmax = round(maximum(sky[:,:,:]-x[:,:,:]),1)
# subplots_adjust(bottom=0.1, right=0.8, top=0.9)
# cax = axes([0.85, 0.13, 0.025, 0.18])
# colorbar(cax=cax)
#############################################################################
#############################################################################


text:
"""
    concordance(class, var, tie)

Computes concordant, tied and discordant pairs.\\
`class` can be either a BitVector or a 2 level categorical target variable in which case
`true` is defined by the last value in sorted sequence.\\
`var` is a Vector of predictor, same length as `class`,\\
`tie` (optional) can be a number (default is 1e-6) that defines a tied region,
or it can be a function that when called with a scalar value will return a tuple
of lower bound and upper bound of a tied region,
useful when you want to do percentage tied region for instance.

Pair-wise comparison between class 1 with class 0 values are made as follows:
- class 1 value > class 0 value is Concordant
- class 1 value ≈ class 0 value (within `tie`) is Tied
- class 1 value < class 0 value is Discordant

Returns:
- concordant, number of concordant comparisons
- tied, number of tied comparisons
- discordant, number of discordant comparisons
- auroc, or C, is (Concordant + 0.5Tied) / Total comparisons; same as numeric integration of ROC curve
- gini, 2C-1, also known as Somer's D, is (Concordant - Discordant) / Total comparisons

Concordance calculation is the same as numeric integration of the ROC curve,
but it allows for fuzzy tied regions which can be useful.

Note:
- Goodman-Kruskal Gamma is (Concordant - Discordant) / (Concordant + Discordant)
- Kendall's Tau is (Concordant - Discordant) / (0.5 x Total count x (Total count - 1))
"""
function concordance(class::BitVector, var::Vector, tie = 1e-6)
    n = length(class)
    n == length(var) || throw(ArgumentError("class and var should be the same length"))
    n1 = sum(class)
    n1 == 0 && throw(ArgumentError("there are no class 1"))
    n0 = n - n1
    n0 == 0 && throw(ArgumentError("there are no class 0"))

    c1 = sort!(var[class])              # sorted var of class 1
    c0 = sort!(var[.!class])            # sorted var of class 0

    conc = 0
    tied = 0
    l_ix = 1
    u_ix = 1
    @inbounds for v in c1               # loop over this usually smaller array
        l_v = v - tie
        u_v = v + tie

        while l_ix <= n0                # first index within [l_v, u_v] window
            l_v <= c0[l_ix] && break
            l_ix += 1                   # could exit as n0 + 1
        end

        while u_ix <= n0                # first index beyond [l_v, u_v] window
            u_v < c0[u_ix] && break
            u_ix += 1                   # could exit as n0 + 1
        end

        conc += (l_ix - 1)              # lower than l_ix are Concordant
        tied += (u_ix - l_ix)           # within window are tied
    end

    pairs = n1 * n0
    disc = pairs - conc - tied
    auc = (conc + 0.5tied) / pairs

    (conc = conc, tied = tied, disc = disc, auc = auc, gini = 2auc - 1)
end

function concordance(class::BitVector, var::Vector, tie::Function)
    n = length(class)
    n == length(var) || throw(ArgumentError("class and var should be the same length"))
    n1 = sum(class)
    n1 == 0 && throw(ArgumentError("there are no class 1"))
    n0 = n - n1
    n0 == 0 && throw(ArgumentError("there are no class 0"))

    c1 = sort!(var[class])              # sorted var of class 1
    c0 = sort!(var[.!class])            # sorted var of class 0

    conc = 0
    tied = 0
    l_ix = 1
    u_ix = 1
    @inbounds for v in c1               # loop over this usually smaller array
        l_v, u_v = tie(v)               # tie scalar function returning lower & upper bound

        while l_ix <= n0                # first index within [l_v, u_v] window
            l_v <= c0[l_ix] && break
            l_ix += 1                   # could exit as n0 + 1
        end

        while u_ix <= n0                # first index beyond [l_v, u_v] window
            u_v < c0[u_ix] && break
            u_ix += 1                   # could exit as n0 + 1
        end

        conc += (l_ix - 1)              # lower than l_ix are Concordant
        tied += (u_ix - l_ix)           # within window are tied
    end

    pairs = n1 * n0
    disc = pairs - conc - tied
    auc = (conc + 0.5tied) / pairs

    (conc = conc, tied = tied, disc = disc, auc = auc, gini = 2auc - 1)
end

function concordance(class::Vector, var::Vector, tie = 1e-6)
    uc = sort!(unique(class))
    length(uc) == 2 || throw(ArgumentError("class must be 2 levels"))

    concordance(class .== uc[2], var, tie)
end

roc = concordance                       # synonym


text:
@testset "utilities - polyeval" begin
# Testing the polynomial evaluation macro:
p = randn(6)
x = 0.33
@test @evalpoly(x, p[1], p[2]) ≈ Psychro.polyeval(x, p, 2) atol=1e-13
@test @evalpoly(x, p[1], p[2], p[3]) ≈ Psychro.polyeval(x, p, 3) atol=1e-13
@test @evalpoly(x, p[1], p[2], p[3], p[4]) ≈ Psychro.polyeval(x, p, 4) atol=1e-13
@test @evalpoly(x, p[1], p[2], p[3], p[4], p[5]) ≈ Psychro.polyeval(x, p, 5) atol=1e-13
@test @evalpoly(x, p[1], p[2], p[3], p[4], p[5], p[6]) ≈ Psychro.polyeval(x, p, 6) atol=1e-13
end


text:
abstract type AbstractAgent end

function action(agent, env, ϵ) end
function update!(agent) end
function run_episode!(agent, env) end
function evaluate(agent, env, n = 100) end
function train!(agent, env, max_episodes = Inf) end


text:
module Terrain

using StaticArrays

using Flight.Geodesy
using Flight.Utils

#TerrainModel does not belong to the Aircraft itself. it must be defined
#separately, and passed as an external data source. the same goes for
#AtmosphereModel. there must be a level above the Aircraft, which will typically
#be the simulation scheduler, that defines both the environmental models and all
#the aircraft participating in the simulation. this may be a block based
#simulation engine or a custom made one.

#the AtmosphericModel contained in the Aircraft should behave as a gateway to
#the actual AtmosphericModel, through which the evolving atmospheric model can
#be queried by the Aircraft for the values at its current location. the
#TerrainModel can be handled similarly, because even if it is not evolving in
#time, it may be an arbitrarily complex terrain database which must serve all
#vehicles. initially, these Model "clients" will be simply references to the
#models themselves, because these will be constant, simple and shared with no
#one else

#an alternative approach, instead of querying the models in a client-server
#architecture, would be simply to provide access references to the environment
#model data and methods and store them in the external data source fields of the
#Aircraft System being simulated. then, the Aircraft can locally evaluate those
#methods for their particular position. in that case, the exchanged information
#is not the product of model evaluation at the Aircraft location, but the data
#required to perform those evaluations within the Aircraft model.

#both scenarios could be realized with an Observer model, with any vehicle in
#the simulation subscribing to the AtmosphericModel and the TerrainModel. each
#time these update, they notify all their subscribers so they can also be
#evolved by the scheduler. if this communication occurs through Channels, each
#Model (Atmospheric, Terrain, Vehicle...) can run on separate Tasks, and each of
#these in a different thread. Thread safety.

export AbstractTerrain, DummyTerrain, HorizontalTerrain
export TerrainData, SurfaceType
export get_terrain_data

@enum SurfaceType DryTarmac WetTarmac IcyTarmac

struct TerrainData
    altitude::Altitude{Orthometric}
    normal::SVector{3,Float64} #NED components, inward pointing
    surface::SurfaceType
end

TerrainData(; altitude = AltOrth(0),
                normal = SVector{3,Float64}(0,0,1),
                surface = DryTarmac) =
    TerrainData(altitude, SVector{3,Float64}(normal), surface)

######################## AbstractTerrain ##########################

abstract type AbstractTerrain end
get_terrain_data(::AbstractTerrain, args...) = throw(MethodError(get_terrain_data, args))

struct DummyTerrain <: AbstractTerrain end

struct HorizontalTerrain <: AbstractTerrain
    altitude::Altitude{Orthometric}
    surface::SurfaceType
end

HorizontalTerrain(; altitude = AltOrth(0), surface = DryTarmac) =
    HorizontalTerrain(altitude, surface)

get_terrain_data(trn::HorizontalTerrain, ::Abstract2DLocation) =
    TerrainData(trn.altitude, SVector{3,Float64}(0,0,1), trn.surface)

end #module


text:
module StatisticalRethinkingBase

end # module


text:
precompile( inject, (TwObj{TwTreeData}, Symbol ) )
precompile( inject, (TwObj{TwTreeData}, ASCIIString ) )
precompile( inject, (TwObj{TwDfTableData}, Symbol ) )
precompile( inject, (TwObj{TwDfTableData}, ASCIIString ) )
precompile( inject, (TwObj{TwEntryData}, Symbol ) )
precompile( inject, (TwObj{TwEntryData}, ASCIIString ) )
precompile( inject, (TwObj{TwScreenData}, Symbol ) )
precompile( inject, (TwObj{TwScreenData}, ASCIIString ) )