text:
# DifferentialEquations
using DifferentialEquations
f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

# DiffEqSensitivity
using DiffEqSensitivity, OrdinaryDiffEq, Zygote
function fiip(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end
p = [1.5,1.0,3.0,1.0]; u0 = [1.0;1.0]
prob = ODEProblem(fiip,u0,(0.0,10.0),p)
sol = concrete_solve(prob,Tsit5())

max_stars_repo_name: michaeldeistler/diffeqtorch

max_issues_repo_name: michaeldeistler/diffeqtorch

max_forks_repo_name: michaeldeistler/diffeqtorch


text:
# code for FastTransforms.jl until that package is compatible with Polynomials v2.0.0

## Use FastTransform for conversion, as  possible
## FastTransforms.kind2string.(0:15)
# Legendre<--Chebyshev
function _convert(
    ::Type{Q},
    p::P,
) where {
    T<:AbstractFloat,
    Q<:Union{Legendre,OrthonormalLegendre},
    P<:Union{Chebyshev{T},OrthonormalChebyshev{T}},
}
    ps = coeffs(p)
    Q(cheb2leg(ps, normcheb=isorthonormal(P), normleg=isorthonormal(Q)))
end

# Chebyshev<--Legendre
function _convert(
    ::Type{Q},
    p::P,
) where {
    T<:AbstractFloat,
    Q<:Union{Chebyshev,OrthonormalChebyshev},
    P<:Union{Legendre{T},OrthonormalLegendre{T}},
}
    ps = coeffs(p)
    Q(leg2cheb(ps, normleg=isorthonormal(P), normcheb=isorthonormal(Q)))
end

# ultraspherical<--ultraspherical
function _convert(
    ::Type{Q},
    p::P,
) where {
    β,
    α,
    T<:AbstractFloat,
    Q<:Union{Gegenbauer{β},OrthonormalGegenbauer{β}},
    P<:Union{Gegenbauer{α,T},OrthonormalGegenbauer{α,T}},
}
    ps = coeffs(p)
    Q(ultra2ultra(ps, α, β, norm1=isorthonormal(P), norm2=isorthonormal(Q)))
end

# Jacobi<--Jacobi
function _convert(
    ::Type{Q},
    p::P,
) where {
    γ,
    δ,
    α,
    β,
    T<:AbstractFloat,
    Q<:Union{Jacobi{γ,δ},OrthonormalJacobi{γ,δ}},
    P<:Union{Jacobi{α,β,T},OrthonormalJacobi{α,β,T}},
}
    ps = coeffs(p)
    Q(jac2jac(ps, α, β, γ, δ, norm1=isorthonormal(P), norm2=isorthonormal(Q)))
end

# Laguerre<--Laguerre
function _convert(
    ::Type{Q},
    p::P,
) where {
    α,
    β,
    T<:AbstractFloat,
    Q<:Union{Laguerre{β},OrthonormalLaguerre{β}},
    P<:Union{Laguerre{α,T},OrthonormalLaguerre{α,T}},
}
    ps = coeffs(p)
    Q(lag2lag(ps, α, β, norm1=isorthonormal(P), norm2=isorthonormal(Q)))
end

# Jacobi<--ultraspherical
function _convert(
    ::Type{Q},
    p::P,
) where {
    γ,
    δ,
    α,
    T<:AbstractFloat,
    Q<:Union{Jacobi{γ,δ},OrthonormalJacobi{γ,δ}},
    P<:Union{Gegenbauer{α,T},OrthonormalGegenbauer{α,T}},
}
    ps = coeffs(p)
    Q(ultra2jac(ps, α, γ, δ, normultra=isorthonormal(P), normjac=isorthonormal(Q)))
end

# ultraspherical<--Jacobi
function _convert(
    ::Type{Q},
    p::P,
) where {
    α,
    γ,
    δ,
    T<:AbstractFloat,
    Q<:Union{Gegenbauer{α,T},OrthonormalGegenbauer{α,T}},
    P<:Union{Jacobi{γ,δ},OrthonormalJacobi{γ,δ}},
}
    ps = coeffs(p)
    Q(jac2ultra(ps, γ, δ, α, normjac=isorthonormal(P), normultra=isorthonormal(Q)))
end

# Jacobi<--Chebyshev
function _convert(
    ::Type{Q},
    p::P,
) where {
    γ,
    δ,
    T<:AbstractFloat,
    Q<:Union{Jacobi{γ,δ},OrthonormalJacobi{γ,δ}},
    P<:Union{Chebyshev{T},OrthonormalChebyshev{T}},
}
    ps = coeffs(p)
    Q(cheb2jac(ps, γ, δ, normcheb=isorthonormal(P), normjac=isorthonormal(Q)))
end

# Chebyshev<--Jacobi
function _convert(
    ::Type{Q},
    p::P,
) where {
    γ,
    δ,
    T<:AbstractFloat,
    Q<:Union{Chebyshev,OrthonormalChebyshev},
    P<:Union{Jacobi{γ,δ,T},OrthonormalJacobi{γ,δ,T}},
}
    ps = coeffs(p)
    Q(jac2cheb(ps, γ, δ, normjac=isorthonormal(P), normcheb=isorthonormal(Q)))
end

# ultraspherical<--Chebyshev
function _convert(
    ::Type{Q},
    p::P,
) where {
    α,
    T<:AbstractFloat,
    Q<:Union{Gegenbauer{α},OrthonormalGegenbauer{α}},
    P<:Union{Chebyshev{T},OrthonormalChebyshev{T}},
}
    ps = coeffs(p)
    Q(cheb2ultra(ps, α, normcheb=isorthonormal(P), normultra=isorthonormal(Q)))
end

# Chebyshev<--ultraspherical
function _convert(
    ::Type{Q},
    p::P,
) where {
    α,
    T<:AbstractFloat,
    Q<:Union{Chebyshev,OrthonormalChebyshev},
    P<:Union{Gegenbauer{α,T},OrthonormalGegenbauer{α,T}},
}
    ps = coeffs(p)
    Q(ultra2cheb(ps, α, normultra=isorthonormal(P), normcheb=isorthonormal(Q)))
end

ext: jl


text:
module GtkServer

export gtk_init, gtk_send, gtk_exit

struct Pipes
    pin  :: Pipe
    pout :: Pipe
end

# Internal Gtk-Server state
struct State
    pipes    :: Pipes
    process  :: Union{Base.Process, Nothing}
    nullmode :: Bool
end

function gtk_init(; nullmode = false, logfile::String = "")
    # Start gtk-server and define pipes
    pin = Pipe()
    pout = Pipe()
    try
        success(`gtk-server`)
        # set up gtk-server command
        if logfile != ""
            cmd = `gtk-server -stdin -log=$logfile`
        else
            cmd = `gtk-server -stdin`
        end
        # execute gtk-server and set up communtb vert term julitb vert term juliacation pipes
        gtkproc = run(pipeline(cmd, stdin = pin, stdout = pout), wait = false)
        process_running(gtkproc) || error("There was a problem starting gtk-server")
        close(pout.in)
        close(pin.out)
        Base.start_reading(pout.out)
        global state = State(Pipes(pin, pout), gtkproc, false)
    catch
        @warn "gtk-server is not available on this platform;\nGtkServer.jl is operating in null mode."
        global state = State(Pipes(pin,pout), nothing, true)
    end
    return state
end

function gtk_send(s::AbstractString)::String
    global state

    if !state.nullmode
        write(state.pipes.pin, s)
        flush(state.pipes.pin)
        yield()
        response = String(readavailable(state.pipes.pout))
    else
        @warn "GtkServer is in null mode; no operation has been performed."
        response = ""
    end
    return response
end

function gtk_exit()::String
    global state

    if !state.nullmode
        response = gtk_send("gtk_server_exit")
        close(state.pipes.pin)
        close(state.pipes.pout)
    else
        @warn "GtkServer is in null mode; no operation has been performed."
        response = ""
    end
    return response
end

end # module


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

mutable struct ExpressRouteCircuitPropertiesFormat <: SwaggerModel
    allowClassicOperations::Any # spec type: Union{ Nothing, Bool } # spec name: allowClassicOperations
    circuitProvisioningState::Any # spec type: Union{ Nothing, String } # spec name: circuitProvisioningState
    serviceProviderProvisioningState::Any # spec type: Union{ Nothing, String } # spec name: serviceProviderProvisioningState
    authorizations::Any # spec type: Union{ Nothing, Vector{ExpressRouteCircuitAuthorization} } # spec name: authorizations
    peerings::Any # spec type: Union{ Nothing, Vector{ExpressRouteCircuitPeering} } # spec name: peerings
    serviceKey::Any # spec type: Union{ Nothing, String } # spec name: serviceKey
    serviceProviderNotes::Any # spec type: Union{ Nothing, String } # spec name: serviceProviderNotes
    serviceProviderProperties::Any # spec type: Union{ Nothing, ExpressRouteCircuitServiceProviderProperties } # spec name: serviceProviderProperties
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    gatewayManagerEtag::Any # spec type: Union{ Nothing, String } # spec name: gatewayManagerEtag

    function ExpressRouteCircuitPropertiesFormat(;allowClassicOperations=nothing, circuitProvisioningState=nothing, serviceProviderProvisioningState=nothing, authorizations=nothing, peerings=nothing, serviceKey=nothing, serviceProviderNotes=nothing, serviceProviderProperties=nothing, provisioningState=nothing, gatewayManagerEtag=nothing)
        o = new()
        validate_property(ExpressRouteCircuitPropertiesFormat, Symbol("allowClassicOperations"), allowClassicOperations)
        setfield!(o, Symbol("allowClassicOperations"), allowClassicOperations)
        validate_property(ExpressRouteCircuitPropertiesFormat, Symbol("circuitProvisioningState"), circuitProvisioningState)
        setfield!(o, Symbol("circuitProvisioningState"), circuitProvisioningState)
        validate_property(ExpressRouteCircuitPropertiesFormat, Symbol("serviceProviderProvisioningState"), serviceProviderProvisioningState)
        setfield!(o, Symbol("serviceProviderProvisioningState"), serviceProviderProvisioningState)
        validate_property(ExpressRouteCircuitPropertiesFormat, Symbol("authorizations"), authorizations)
        setfield!(o, Symbol("authorizations"), authorizations)
        validate_property(ExpressRouteCircuitPropertiesFormat, Symbol("peerings"), peerings)
        setfield!(o, Symbol("peerings"), peerings)
        validate_property(ExpressRouteCircuitPropertiesFormat, Symbol("serviceKey"), serviceKey)
        setfield!(o, Symbol("serviceKey"), serviceKey)
        validate_property(ExpressRouteCircuitPropertiesFormat, Symbol("serviceProviderNotes"), serviceProviderNotes)
        setfield!(o, Symbol("serviceProviderNotes"), serviceProviderNotes)
        validate_property(ExpressRouteCircuitPropertiesFormat, Symbol("serviceProviderProperties"), serviceProviderProperties)
        setfield!(o, Symbol("serviceProviderProperties"), serviceProviderProperties)
        validate_property(ExpressRouteCircuitPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(ExpressRouteCircuitPropertiesFormat, Symbol("gatewayManagerEtag"), gatewayManagerEtag)
        setfield!(o, Symbol("gatewayManagerEtag"), gatewayManagerEtag)
        o
    end
end # type ExpressRouteCircuitPropertiesFormat

const _property_map_ExpressRouteCircuitPropertiesFormat = Dict{Symbol,Symbol}(Symbol("allowClassicOperations")=>Symbol("allowClassicOperations"), Symbol("circuitProvisioningState")=>Symbol("circuitProvisioningState"), Symbol("serviceProviderProvisioningState")=>Symbol("serviceProviderProvisioningState"), Symbol("authorizations")=>Symbol("authorizations"), Symbol("peerings")=>Symbol("peerings"), Symbol("serviceKey")=>Symbol("serviceKey"), Symbol("serviceProviderNotes")=>Symbol("serviceProviderNotes"), Symbol("serviceProviderProperties")=>Symbol("serviceProviderProperties"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("gatewayManagerEtag")=>Symbol("gatewayManagerEtag"))
const _property_types_ExpressRouteCircuitPropertiesFormat = Dict{Symbol,String}(Symbol("allowClassicOperations")=>"Bool", Symbol("circuitProvisioningState")=>"String", Symbol("serviceProviderProvisioningState")=>"String", Symbol("authorizations")=>"Vector{ExpressRouteCircuitAuthorization}", Symbol("peerings")=>"Vector{ExpressRouteCircuitPeering}", Symbol("serviceKey")=>"String", Symbol("serviceProviderNotes")=>"String", Symbol("serviceProviderProperties")=>"ExpressRouteCircuitServiceProviderProperties", Symbol("provisioningState")=>"String", Symbol("gatewayManagerEtag")=>"String")
Base.propertynames(::Type{ ExpressRouteCircuitPropertiesFormat }) = collect(keys(_property_map_ExpressRouteCircuitPropertiesFormat))
Swagger.property_type(::Type{ ExpressRouteCircuitPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ExpressRouteCircuitPropertiesFormat[name]))}
Swagger.field_name(::Type{ ExpressRouteCircuitPropertiesFormat }, property_name::Symbol) =  _property_map_ExpressRouteCircuitPropertiesFormat[property_name]

const _allowed_ExpressRouteCircuitPropertiesFormat_serviceProviderProvisioningState = ["NotProvisioned", "Provisioning", "Provisioned", "Deprovisioning"]

function check_required(o::ExpressRouteCircuitPropertiesFormat)
    true
end

function validate_property(::Type{ ExpressRouteCircuitPropertiesFormat }, name::Symbol, val)
    if name === Symbol("serviceProviderProvisioningState")
        Swagger.validate_param(name, "ExpressRouteCircuitPropertiesFormat", :enum, val, _allowed_ExpressRouteCircuitPropertiesFormat_serviceProviderProvisioningState)
    end
end


text:
# compare stack depth to trampoline depth

function stack(n, m)
    if n > m
        n
    else
        if n % 1000 == 0
#            println("stack $n")
        end
        n + stack(n+1, m)
    end
end

stack(0, 10)
@time println(stack(0, 100_000))
# stack limit is somewhere around 100,000 (certainly less than 200,000)

abstract Msg

type Call<:Msg
    before::Function
    after::Function
    value::Int
end

type Return<:Msg
    value::Int
end

function inc(n, m)
    if n > m
        Return(n)
    else
        if n % 1000 == 0
#            println("trampoline $n")
        end
        Call(inc, (x, m) -> Return(n+x), n+1)
    end
end

function sum(n, m)
    Return(n)
end
function trampoline(n, m)
    stack = Function[inc]
    while length(stack) > 0
        f = pop!(stack)
        msg = f(n, m)
        if isa(msg, Call)
            push!(stack, msg.after)
            push!(stack, msg.before)
        end
        n = msg.value
    end
    n
end

trampoline(0, 10)
@time println(trampoline(0, 100_000))
@time println(trampoline(0, 1_000_000))
@time println(trampoline(0, 10_000_000))


text:
module Cn2An

number_cn2an_map = Dict(
    "零" => 0,
    "一" => 1,
    "二" => 2,
    "三" => 3,
    "四" => 4,
    "五" => 5,
    "六" => 6,
    "七" => 7,
    "八" => 8,
    "九" => 9
)

unit_cn2an_map = Dict(
    "十" => 10,
    "百" => 100,
    "千" => 1000,
    "万" => 10000,
    "亿" => 100000000
)

function cn2an(inputs::String)
    output = 0
    unit = 1
    ten_thousand_unit = 1
    num = 0

    for (index, cn_num) in enumerate(reverse(inputs))
        cn_num = string(cn_num)

        if haskey(number_cn2an_map, cn_num)
            num = number_cn2an_map[cn_num]
            output = output + num * unit

        elseif haskey(unit_cn2an_map, cn_num)
            unit = unit_cn2an_map[cn_num]

            if unit % 10000 == 0
                if unit > ten_thousand_unit
                    ten_thousand_unit = unit
                else
                    ten_thousand_unit = unit * ten_thousand_unit
                    unit = ten_thousand_unit
                end
            end

            if unit < ten_thousand_unit
                unit = ten_thousand_unit * unit
            end

        else
            throw(UndefVarError(cn_num))
        end

    end

    return output
end

number_an2cn_map = Dict(
    0 => "零",
    1 => "一",
    2 => "二",
    3 => "三",
    4 => "四",
    5 => "五",
    6 => "六",
    7 => "七",
    8 => "八",
    9 => "九"
)

unit_an2cn_order_list = ["", "十", "百", "千", "万", "十", "百", "千", "亿", "十", "百", "千", "万", "十", "百", "千"]

function an2cn(inputs::Int)
    inputs = string(inputs)
    len_input = length(inputs)
    output = ""
    for (index, an_num) in enumerate(inputs)
        an_num = parse(Int, an_num)
        if an_num > 0
            output = output * number_an2cn_map[an_num] * unit_an2cn_order_list[len_input - index + 1]
        else
            if (len_input - index) % 4 == 0
                output = output * number_an2cn_map[an_num] * unit_an2cn_order_list[len_input - index + 1]
            end

            if index > 1 && string(output[end]) != "零"
                output = output * number_an2cn_map[an_num]
            end
        end
    end
    output = replace(output, r"零零" => "零")
    output = replace(output, r"零万" => "万")
    output = replace(output, r"零亿" => "亿")
    output = replace(output, r"^零" => "")
    output = replace(output, r"零$" => "")

    return output
end

end # module


text:
"""
Convert from geodetic coordinates to local East, North, Up (ENU) coordinates.
"""
function geodetic2enu(lon, lat, h, lon0, lat0, h0)
	#
	D1 = mapproject([lon lat h], geod2ecef=true);
	D0 = mapproject([lon0 lat0 h0], geod2ecef=true);
	d  = D1[1].data .- D0[1].data;
	xEast, yNorth, zUp = ecef2enuv(view(d, :, 1), view(d, :, 2), view(d, :, 3), lon0, lat0)
end

function ecef2enuv(u, v, w, lon0, lat0)
	# Adapted from Octave, that adapted from anonymous contributor that probably adapted from Matlab
	t     =  cosd(lon0) .* u + sind(lon0) .* v;
	east  = -sind(lon0) .* u + cosd(lon0) .* v;
	up    =  cosd(lat0) .* t + sind(lat0) .* w;
	north = -sind(lat0) .* t + cosd(lat0) .* w;
	return east, north, up
end


text:
using RasterDataSources, URIs, Test, Dates
using RasterDataSources: rastername, rasterpath, zipurl, zipname, zippath

@testset "WorldClim Climate" begin

    zip_url = URI(scheme="https", host="biogeo.ucdavis.edu", path="/data/worldclim/v2.1/base/wc2.1_10m_wind.zip")
    @test zipurl(WorldClim{Climate}, :wind; res="10m") == zip_url
    @test zipname(WorldClim{Climate}, :wind; res="10m") == "wc2.1_10m_wind.zip"

    raster_file = joinpath(ENV["RASTERDATASOURCES_PATH"], "WorldClim", "Climate", "wind", "wc2.1_10m_wind_01.tif")
    @test rasterpath(WorldClim{Climate}, :wind; month=1, res="10m") == raster_file
    @test getraster(WorldClim{Climate}, (:wind,); month=1, res="10m") == (wind=raster_file,)
    @test getraster(WorldClim{Climate}, :wind; month=1:1, res="10m") == [raster_file]
    @test getraster(WorldClim{Climate}, (:wind,); month=1:1, res="10m") == [(; wind=raster_file)]
    @test getraster(WorldClim{Climate}, [:wind]; month=1:1, res="10m") == [(; wind=raster_file)]
    @test isfile(raster_file)
end

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
# About the construction of a grid
#
# - The domain boundaries are always primal grid planes.
#
# - The dual grid planes are always mid-sections between adjacent primal grid planes.
#
# - In each Cartesian axis, the negative- and positive-end boundary conditions are always
# chosen the same.  For example, if the negative end is Bloch, the positive end must be also
# Bloch.  Similarly, if the negative end is PEC, the positive end is also PEC.  (More on
# this below.)
#
# - The primal and dual grid planes are not necessarily associated with the E- and H-field
# planes.  Which of the primal and dual grid planes to use for the E- and H-fields is
# decided by the choice of the boundary condition.
#
#       - The PEC and PMC symmetry boundary conditions are realized by putting the relavant
#       fields on the domain boundary as the tangential components.  Therefore, we use the
#       primal grid planes as the E-field planes for the PEC boundary condition, and as the
#       H-field planes for the PMC boundary condition.
#
#       - On the other hand, for the Bloch boundary condition we can use the primal grid
#       planes as either of the E- and H-field planes.  Note that Maxwell's equations are
#       ∇_h × μ⁻¹ ∇_e × E - ω² ε E = -i ω J.  Once we decide that in the w-axis direciton
#       (w = x,y,z) we are going to use the primal grid planes as the E-field planes, we can
#       simply use the forward difference (which is for primal fields) as ∂/∂w in
#       constructing ∇_e.
#
# - The notion of U, V, and ξ, η is now inappropriate.  I used to use U to indicate the
# primal field, but there is no such thing like the primal field, because the E-field, for
# example, does not have to be completely on the primal grid planes.  The notion makes sense
# only in one Cartesian direction.

# About the symmetry boundary condition
#
# - The symmetry boundary is used to implement PEC and PMC on the boundaries of the domain.
# The symmetry boundary condition zeros the tangential fields to the boundary.  Because the
# domain boundaries are always primal grid planes, this means that PEC is implemented by
# using the primal grid planes as the E-tangential planes and the dual grid plane as the
# H-tangential planes.
#
# - An important property of the symmetry boundary plane is that it simultaneously zeros the
# tangential field and normal field defined on the boundary.  Here, the two fields are of
# different kinds.  For example, if the primal grid planes are E-tangential planes, on the
# boundary planes tangential E-fields and normal H-fields are defined.  The symmetry
# boundary condition zeros both kinds of the fields.
#
# The fact that symmetry boundary condition zeros *any* fields defined on the boundary
# greatly simplify the implementation of the effect of the boundary condition in
# differential operators.  For example, for the PEC boundary condition, the forward
# derivative operator subtracts fields between primal grid planes, but these fields can be
# either tangential E-fields or normal H-fields.  We have to create the derivative operators
# to handle the boundary conditions properly (such that they zero the fields that must be
# zero), but this property allows us to use the same derivative operator for both E- and H-
# fields.
#
# - You never need to mix the Bloch boundary condition and symmetry boundary condition like
# PEC and PMC, because by using the symmetry boundary condition on both ends of the domain
# your domain is automatically repeated to infinity.  (In other words, symmetry boundary
# condition is a way to simulate a periodic structure whose unit cell is mirror-symmetric.
# Note that inside a mirror-symmetric structure, the field distribution can be mirror-
# symmetric or antisymmetric.)
#
# In some cases you may want to put PEC on one end and PMC on the other end of the domain to
# simulate a half of a mirror-symmetric unit cell of a periodic structure where the field
# distribution is symmetric at the center plane of the unhalved unit cell but antisymmetric
# at the boundaries of the unit cell.  At first glance, supporting such a mixed symmetry
# boundary condition would be difficult to support, because PEC and PMC are realized by
# zeroing the tangential E- and H-fields.  Implementing PEC on the negative end and PMC on
# the positive end means therefore, e.g., to start the domain with a primal grid plane
# (containing tangential E-fields) and end the domain with a dual grid plane (containing
# tangential H-fields), which requires putting only a half grid on the positive end and thus
# not ideal.  In fact, we can still keep a whole integral number of grid cells to implement
# the PMC on the negative end and PEC on the positive end by using the cool trick I used
# in FD3D.  (We cannot implement PEC on the negative end and PMC on the positive end even
# with that trick.)  However, supporting such a mixed boundary condition makes the code
# unnecesasrily convoluted, so I will not support it for now.  (Consider supporting it by
# passing something like an "ismixed" flag.)

# About matrix bandwidth vs. boundary condition
#
# Sparsity pattern of matrix with Bloch boundary condition:
# ⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙
# ⠀⠈⠻⣦⡀
# ⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦
#
# Note the nonzero entries at the top-right and bottom-left corners.
#
# Sparsity pattern of matrix with symmetry boundary condition:
# ⠻⣦⡀
# ⠀⠈⠻⣦⡀
# ⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀
# ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦
#
# No nonzero entries at those corners.  Hence, the symmetry boundary condition produces much
# narrower bandwidth, which leads less fill-ins in matrix factors.

export Grid  # types
export isproper_blochphase, lghost  # functions

# Below, the transformation τ refers to the transformation of a ghost point back to the
# non-ghost point that the material parameter at the ghost point comes from.  For the Bloch
# boundary condition, τ transforms the ghost point to the other side of the domain; for the
# symmetry boundary condition, τ transforms the ghost point to its mirror symmetry position
# with respect to the boundary.
#
# All the fields of Ghosted include ghost points.  Therefore, their lengths are N.+1, rather
# than N.  (Even thought σ seems to deal with reflection of ghost points, it doesn't; it has
# N entries, so it is defined in Grid rather than in Ghosted.)
struct Ghosted{K}
    l::Tuple2{NTuple{K,VecFloat}}  # l[PRIM][k] = primal vertex locations in k-direction with ghost points
    τl::Tuple2{NTuple{K,VecFloat}}  # τl[PRIM][k] = primal vertex locations in k-direction with transformed ghost points into domain; all τl are in domain
    τind::Tuple2{NTuple{K,VecInt}}  # τind[PRIM][k] = indices of Ghosted.l corresponding to transformed points in to domain by boundary conditions
    ∆τ::Tuple2{NTuple{K,VecFloat}}  # ∆τ[PRIM][k] = amount of "shift" to get points into domain by Bloch boundary conditions (ignoring transformation by reflection by symmetry boundary, which is handled by Grid.σ)
end

# Consider storing lvxlbounds, which is
# ([(ldual₀,ldual₁), (ldual₁,ldual₂), ..., (ldualₙ₋₁,ldualₙ)],
#  [(lprim₁,lprim₂), ..., (lprimₙ₋₁,lprimₙ), (lprimₙ,lprimₙ₊₁)]
# The purpose of this field is easy retrieval of voxel bounds surrounding a voxel center.
# For example, for a voxel center l[gt][i], the surrounding voxel bounds are lvxlbounds[gt][i],
# regardless of the value of gt.  You don't need to adjust i by ±1.
#
# Another option.  To get the voxel bounds for voxel centers at l[gt], take l[alter(gt)],
# and prepend or append the ghost point appropriately.  Let the resulting array lvxlbounds.
# Then, for a voxel center l[gt][i], the voxel bounds are lvxlbounds[[i,i+1]], regardless of
# the value of gt.
struct Grid{K}
    N::SVec{K,Int}  # N[k] = number of grid cells in k-direction
    L::SVec{K,Float}  # L[k] = length of grid (domain) in k-direction
    l::Tuple2{NTuple{K,VecFloat}}  # l[PRIM][k] = primal vertex locations in k-direction
    ∆l::Tuple2{NTuple{K,VecFloat}}  # ∆l[PRIM][k] = (∆l at primal vertices in w) == diff(l[DUAL][k] including ghost point)
    isbloch::SVec{K,Bool}  # isbloch[k]: true if boundary condition in k-direction is Bloch
    σ::Tuple2{NTuple{K,VecBool}}  # false only for non-ghost points exactly on symmetry boundary (= first point in primal grid)
    bounds::Tuple2{SVec{K,Float}}  # bounds[NEG][k] = boundary of domain at (-) end in k-direction
    ghosted::Ghosted{K}  # data related to ghost points
end

# Constructor taking non-static vectors.
Grid(lprim::NTuple{K,AbsVecReal}, isbloch::AbsVecBool) where {K} = Grid(lprim, SVec{K}(isbloch))

Base.size(grid::Grid) = grid.N.data
Base.length(grid::Grid) = prod(size(grid))

# Constructor calling the inner constructor.
function Grid(lprim::NTuple{K,AbsVecReal},  # primal grid plane locations, including both domain boundaries
              isbloch::SVec{K,Bool}) where {K}
    all(issorted.(lprim)) || throw(ArgumentError("all entry vectors of lprim = $(lprim) must be sorted."))

    # For array inputs, create separate copies.
    lprim = collect.(float.(lprim))  # collect instead of deepcopy to turn range (e.g., 1.0:10.0) into vector

    ldual = movingavg.(lprim)  # NTuple{K,VecFloat}

    lbound = SVec(lprim)  # SVec{K,Vector{<:Real}}; to make broadcast and map on lprim to produce SVec
    bounds = (getindex.(lbound,1), getindex.(lbound,lastindex.(lbound)))  # Tuple2{SVec{K,<:Real}}
    L = bounds[nP] - bounds[nN]  # SVec{K,<:Real}

    # Prepare lprim and ldual to calculate ∆lprim and ∆ldual.
    # Note that ∆lprim is not ∆'s or lprim, but ∆'s defined at lprim locations.
    for k = 1:K
        # To calculate ∆'s at lprim, we need one more ldual outside the negative boundary.
        ldual₀ = isbloch[k] ? ldual[k][end]-L[k] : 2bounds[nN][k]-ldual[k][1]
        prepend!(ldual[k], ldual₀)  # length(ldual[k]) = N[k]+1
    end

    # Below, (∆lprim, ∆ldual) is not (diff(lprim), diff(ldual)), but swapped.
    ∆lprim = diff.(ldual)  # NTuple{K,VecFloat}; length(∆lprim[k]) = N[k]
    ∆ldual = diff.(lprim)  # NTuple{K,VecFloat}; length(∆ldual[k]) = N[k]
    ∆l = (∆lprim, ∆ldual)

    # Set N (number of grid cells along the axis).
    @assert length.(∆l[nPR]) == length.(∆l[nDL])  # lprim, ldual, ∆lprim, ∆ldual have the same length
    N = SVec(length.(∆l[nPR]))  # SInt{3}

    # Find the locations of the ghost points transformed into the domain by the boundary
    # conditions.  Note that regardless of the boundary condition, the primal ghost point is
    # on the positive side (so it is lprim[end]) and the dual ghost point is on the negative
    # side (so it is ldual[1]).
    # Make sure these locations are where the objects assigned to the ghost points are found.

    # Construct an instance of Ghosted.  Note that lprim and ldual include ghost points now.
    τind = (map(n->collect(1:n+1), N.data), map(n->collect(1:n+1), N.data))  # Tuple23{VecInt}
    τindg_prim = .!isbloch.*N .+ 1  # SInt{K}: N+1 for symmetry; 1 for Bloch
    τindg_dual = isbloch.*N .+ 1 .+ .!isbloch  # SInt{K}: 2 for symmetry; N+1 for Bloch
    τindg = (τindg_prim, τindg_dual)  # Tuple2{SInt{K}}

    τl = (deepcopy(lprim), deepcopy(ldual))

    ∆τ = (zeros.((N.+1).data), zeros.((N.+1).data))  # Tuple2{NTuple{K,VecFloat}}

    for k = 1:K
        τind[nPR][k][end] = τindg[nPR][k]  # N+1 for symmetry; 1 for Bloch
        τind[nDL][k][1] = τindg[nDL][k]  # 2 for symmetry; N+1 for Bloch

        # Below, getindex() does not return an index, but the value at the given index.
        τl[nPR][k][end] = getindex(lprim[k], τind[nPR][k][end])  # lprim[k][τind[nPR][k][end]]
        τl[nDL][k][1] = getindex(ldual[k], τind[nDL][k][1])  # ldual[k][τind[nDL][k][1]]

        ∆τ[nPR][k][end] = -L[k] * isbloch[k]
        ∆τ[nDL][k][1] = L[k] * isbloch[k]
    end

    ghosted = Ghosted{K}((deepcopy(lprim), deepcopy(ldual)), τl, τind, ∆τ)

    # Construct σ, which is false on non-Bloch (= symmetry) boundaries.
    # Note that σ's are vectors of length-N and therefore do not include ghost points.
    σprim = ones.(Bool, N.data)  # NTuple{K,VecBool}; all filled with 'True'
    σdual = ones.(Bool, N.data)  # NTuple{K,VecBool}; all filled with 'True'
    for k = 1:K
        # Primal grid points are on the symmetry boundary if they are the first (non-ghost)
        # points and the symmetry boundary is used.
        σprim[k][1] = isbloch[k]

        # Dual grid points are never on any domain boundary, so no change on σdual.
    end
    σ = (σprim, σdual)

    # Eliminate ghost points.  Regardless of the boundary type, the primal ghost point
    # is on the positive side (so lprim[end] is dropped) and the dual ghost point is on
    # the negative side (so ldual[1] is dropped).
    pop!.(lprim)
    popfirst!.(ldual)
    @assert length.(lprim) == length.(ldual) == N.data  # lprim, ldual, ∆lprim, ∆ldual have the same length
    l = (lprim, ldual)

    return Grid{K}(N, L, l, ∆l, isbloch, σ, bounds, ghosted)
end

# To-do
# When needed, implement `section` (or `cross_section`) that returns Grid{K₁} from Grid{K₀}
# for K₁ < K₀.  Provide the directions (Axis) across which the cross section in taken.  Use
# the inner constructor of Grid to construct the cross section.

# Additional functions
Base.ndims(::Grid{K}) where {K} = K
Base.in(l::SVec{K}, g::Grid{K}) where {K} = all(g.bounds[nN] .≤ l .≤ g.bounds[nP])

# For the Bloch boundary, check if e⁻ⁱᵏᴸ is indeed a phase factor (i.e., amplitude = 1).
# For the non-Bloch (= symmetry) boundary, check if e⁻ⁱᵏᴸ = 1.
#
# e⁻ⁱᵏᴸ is the accumulated phase from the negative-end boundary to the positive-end one.  L
# is the distance between the boundaries (= domain size).  It is exp(-ikL), not exp(+ikL),
# because we assume the exp(+iωt) time dependence and thus the position dependence is
# exp(-ik⋅r).
isproper_blochphase(e⁻ⁱᵏᴸ::Number, isbloch::Bool) = (isbloch && abs(e⁻ⁱᵏᴸ)==1) || (!isbloch && e⁻ⁱᵏᴸ==1)
isproper_blochphase(e⁻ⁱᵏᴸ::SVec{K,<:Number}, isbloch::SVec{K,Bool}) where {K} =
    all(ntuple(k->isproper_blochphase(isbloch[k], e⁻ⁱᵏᴸ[k]), Val(K)))


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

function permute_states(fst::WFST{W,A,D},perm::Vector{D}) where {W,A,D,I}
  if length(perm) != length(fst)
    throw(ErrorException("Number of states of WFST $(lenth(fst)) not equal to permutation length  $(length(perm))"))
  end
  isym=get_isym(fst)
  osym=get_osym(fst)
  iperm = invperm(perm)
  initials = Dict(iperm[k] => fst.initials[k] for k in keys(fst.initials))
  finals = Dict(iperm[k] => fst.finals[k] for k in keys(fst.finals))
  states = [Vector{A}() for i in 1:size(fst,1)]
  for (j,i) in enumerate(perm)
    for arc in fst[i]
      push!(states[j],
            Arc{W,D}(get_ilabel(arc),get_olabel(arc),
                     get_weight(arc),iperm[get_nextstate(arc)])
           )
    end
  end
  return WFST(states,initials,finals,isym,osym)
end


text:
using Test
using BCVUMPS
using ADBCVUMPS
using ADBCVUMPS:diaglocal

@testset "hamiltonianmodels" for Ni = [1,2,3], Nj = [1,2,3]
    @test Ising(Ni,Nj) isa HamiltonianModel
    @test TFIsing(Ni,Nj,1.0) isa HamiltonianModel
    @test Heisenberg(Ni,Nj) isa HamiltonianModel
    @test diaglocal(Ni,Nj,[1.,-1]) isa HamiltonianModel
end


text:
###############################################################################
"""
Provides utilities for removal of radial distortion, by predicting the true u,v
coordinates of a pixel from the distorted coordinates, and the intrinsic
calibration parameters of the camera.

Also provides utilities to simulate radial distoriton.
"""
module RemoveImageDistortion
###############################################################################

###############################################################################
import ExtrinsicCameraCalibration: IntrinsicCalibrationParameters
# using Statistics
###############################################################################

###############################################################################
export RemoveRadialDistortion, SimulateRadialDistortion
###############################################################################

###############################################################################
function RemoveRadialDistortion(u_d::Array, v_d::Array, ι::IntrinsicCalibrationParameters)
    return RemoveRadialDistortion(u_d::Array, v_d::Array, ι.radial_distortion_poly_degrees,
                                  ι.radial_distortion_poly_coeffs, ι.u₀, ι.v₀)
end
function RemoveRadialDistortion(u_d::Array, v_d::Array, Degrees, Coeffs, u₀, v₀)
    # Convert the *observed* pixel to polar coordinates:
    ρ_d = abs.((u_d .- u₀) .+ (v_d .- v₀)im)
    Θ_d = angle.((u_d .- u₀) .+ (v_d .- v₀)im)
    # Calculate the Δρ correction needed at each point:
    Δρ = Polyvalc(ρ_d, Coeffs, Degrees)
    Δu = cos.(Θ_d) .* Δρ
    Δv = sin.(Θ_d) .* Δρ
    # Apply the correction to ρ and convert back to pixel coordinates:
    u_p = u_d - Δu
    v_p = v_d - Δv
    return u_p,v_p
end
###############################################################################

###############################################################################
function SimulateRadialDistortion(uv_p::Array, ι::IntrinsicCalibrationParameters)
    return SimulateRadialDistortion(u_p, v_p, ι)
end
function SimulateRadialDistortion(u_p::Array, v_p::Array, ι::IntrinsicCalibrationParameters)
    return SimulateRadialDistortion(u_p, v_p,
                                    ι.radial_distortion_poly_degrees,
                                    ι.radial_distortion_poly_coeffs,
                                    ι.u₀, ι.v₀)
end
function SimulateRadialDistortion(u_p::Array, v_p::Array, Degrees, Coeffs, u₀, v₀)
    ρ_p = abs.((u_p .- u₀) .+ (v_p .- v₀)im)
    Θ_p = angle.((u_p .- u₀) .+ (v_p .- v₀)im)
    # Calculate the Δρ correction needed at each point:
    Δρ = Polyvalc(ρ_p, Coeffs, Degrees)
    Δu = cos.(Θ_p) .* Δρ
    Δv = sin.(Θ_p) .* Δρ
    # Apply the correction to ρ and convert back to pixel coordinates:
    u_d = u_p + Δu
    v_d = v_p + Δv
    return u_d,v_d
end
###############################################################################

###############################################################################
# Polyfitc is similar to polyfit, but allows choosing powers (5,3,1, etc.)
# (based on Matlab file exchange Polyfitc calculations):
function Polyfitc(X,Y,Degrees)
    lsqM = ones(length(X), length(Degrees))
    for n = 1:length(Degrees)
        lsqM[:, n] = X.^Degrees[n]
    end
    Coeffs = (lsqM \ Y)'
    return Coeffs[:]
end
# ----------
function Polyvalc(X, Coeffs, Degrees)
    val(x,Coeffs,Degrees) = sum( [ Coeffs[i]*x^(Degrees[i]) for i in 1:length(Coeffs) ] )
    Y = [ val(x,Coeffs,Degrees) for x in X ]
end
###############################################################################

###############################################################################
###############################################################################
###############################################################################
end                             # end of module
###############################################################################
###############################################################################
###############################################################################

max_stars_repo_name: saschweitzer/camera-calibration

max_issues_repo_name: saschweitzer/camera-calibration

max_forks_repo_name: saschweitzer/camera-calibration


text:
function main()
    result = 0
    for i in 2:1_000_000
        if i == reduce(+, map((x) -> x^5, digits(i)))
            result += i
        end
    end
    @printf("sum of terms = %d\n", result)
end

main()


text:
function didyoumean(A,x)
    d = Float64[ compare(Levenshtein(), lowercase(a), lowercase(x)) for a in A ]
    i = indmax(d)
    A[i]
end

function _addmissing(df::AbstractDataFrame, cols::AbstractVector{Symbol})
    df = copy(df)
    for c in cols
        c in names(df) || (df[c] = missing)
    end
    df
end

function _vcat(args::AbstractDataFrame...)
    cols = unique(vcat(names.(args)...))
    vcat( [ _addmissing(df,cols) for df in args ]... )
end

function printtable(io::IO, T::DataFrame)
    A = fill("", size(T,1)+1, 11) # +1 for header
    P = Vector{Function}(undef, 11) # determines left or right padding
    S = fill("", 11)

    P[ 1]=rpad; S[ 1]=" "; A[:, 1]=vcat("Lag", map(string, T[:Team]))
    P[ 2]=lpad; S[ 2]=" "; A[:, 2]=vcat("Säs", map(string, T[:Seasons]))
    P[ 3]=lpad; S[ 3]=" "; A[:, 3]=vcat("S",   map(string, T[:GamesPlayed]))
    P[ 4]=lpad; S[ 4]=" "; A[:, 4]=vcat("V",   map(string, T[:Wins]))
    P[ 5]=lpad; S[ 5]=" "; A[:, 5]=vcat("O",   map(string, T[:Draws]))
    P[ 6]=lpad; S[ 6]=" "; A[:, 6]=vcat("F",   map(string, T[:Losses]))
    P[ 7]=lpad; S[ 7]="-"; A[:, 7]=vcat("GM",  map(string, T[:GoalsFor]))
    P[ 8]=rpad; S[ 8]=" "; A[:, 8]=vcat("IM",  map(string, T[:GoalsAgainst]))
    P[ 9]=lpad; S[ 9]=" "; A[:, 9]=vcat("D",   map(string, T[:GoalDifference]))
    P[10]=lpad; S[10]=" "; A[:,10]=vcat("P",   map(string, T[:Points]))
    P[11]=lpad; S[11]="";  A[:,11]=vcat("PS",  map(x->@sprintf("%.2f",x), T[:PointAverage]))

    # padding
    for (i,f) in enumerate(P)
        A[:,i] = f.(A[:,i], maximum(length.(A[:,i])))
    end

    for i=1:size(A,1)
        for j=1:size(A,2)
            print(io, A[i,j], S[j])
        end
        println(io)
    end
end
printtable(T::DataFrame) = printtable(stdout, T)

function printgames(io::IO, G::DataFrame)
    size(G,1)==0 && return
    A = fill("", size(G,1), 6)
    P = Vector{Function}(undef,6) # determines left or right padding
    S = fill("",6)

    P[ 1]=rpad; S[ 1]=" ";   A[:, 1]=map(string, G[:Date])
    P[ 2]=rpad; S[ 2]=" - "; A[:, 2]=map(string, G[:HomeTeam])
    P[ 3]=rpad; S[ 3]=" ";   A[:, 3]=map(string, G[:AwayTeam])
    P[ 4]=lpad; S[ 4]="-";   A[:, 4]=map(string, G[:HomeScore])
    P[ 5]=rpad; S[ 5]="";    A[:, 5]=map(string, G[:AwayScore])
    P[ 6]=rpad; S[ 6]="";    A[:, 6].="" # for extra notes

    # handle unusual extra fields
    if :Outcome in names(G)
        for i=1:size(A,1)
            o = G[i,:Outcome]
            if !ismissing(o)
                o=="HomeVictory" && (A[i,6] = string(A[i,6], ' ', G[i,:HomeTeam], " tilldömd segern."))
                o=="AwayVictory" && (A[i,6] = string(A[i,6], ' ', G[i,:AwayTeam], " tilldömd segern."))
            end
        end
    end
    if :PointsDeductedHome in names(G)
        for i=1:size(A,1)
            p = G[i,:PointsDeductedHome]
            ismissing(p) || (A[i,6] = string(A[i,6], ' ', p, " poängs avdrag för ", G[i,:HomeTeam], '.'))
        end
    end
    if :PointsDeductedAway in names(G)
        for i=1:size(A,1)
            p = G[i,:PointsDeductedAway]
            ismissing(p) || (A[i,6] = string(A[i,6], ' ', p, " poängs avdrag för ", G[i,:AwayTeam], '.'))
        end
    end

    # padding
    for (i,f) in enumerate(P)
        A[:,i] = f.(A[:,i], maximum(length.(A[:,i])))
    end

    for i=1:size(A,1)
        for j=1:size(A,2)
            print(io, A[i,j], S[j])
        end
        println(io)
    end
end
printgames(G::DataFrame) = printgames(stdout, G)


text:
module StayingPutSim

    using DataFrames
    using StatFiles
    using IterableTables
    using IteratorInterfaceExtensions
    using TableTraits
    using CSV
    using GLM
    # using Plots
    using StatsModels
    using StatsBase
    using Statistics
    using Query
    using DataFramesMeta

    using TableTraits

    # local imports
    using Utils
    using GlobalDecls

    #import CareData: Carer, YP, create_base_datasets, DataSettings, get_yp, make_yp_frame,
#        make_carer_frame,  DataSettings, CPIINDEX, CarerOutcomes, uprate,  AFC_SURVEY_YEAR,
#        EmploymentStatus

    using CareData
    using FosterParameters
    using LAModelData: get_pay_class_from_region_code
    export do_one_calc, override_outcomes!, do_basic_calc, track_series

    function do_basic_calc(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes

        fee1 = ismissing(council_data[:sp_fee_1]) ? 0.0 : council_data[:sp_fee_1]
        allow1 = ismissing(council_data[:sp_allowance_1]) ? 0.0 : council_data[:sp_allowance_1]
        # fixme do something with these
        fee2 = ismissing(council_data[:sp_fee_2]) ? 0.0 : council_data[:sp_fee_2]
        allow2 = ismissing(council_data[:sp_allowance_2]) ? 0.0 : council_data[:sp_allowance_2]
        fee3 = ismissing(council_data[:sp_fee_3]) ? 0.0 : council_data[:sp_fee_3]
        allow3 = ismissing(council_data[:sp_allowance_3]) ? 0.0 : council_data[:sp_allowance_3]
        fee4 = ismissing(council_data[:sp_fee_3]) ? 0.0 : council_data[:sp_fee_4]

        fee1 = uprate( fee1, year, AFC_SURVEY_YEAR )
        allow1 = uprate( allow1, year, AFC_SURVEY_YEAR )
        @assert (fee1+allow1) > 0.0
        CarerOutcomes(
            fee          = fee1,
            allowance    = allow1,
            housing_cont = 0.0,
            other_cont   = 0.0 )
    end # do_basic_calc

    FUNCS = Dict()

    # so taper(10,[1.0,0.5,0.25]) => 10, 5, 2.5
    function taper( amount :: Float64, rates :: Vector{Float64} ) :: Vector{Float64}
        n = size( rates )[1]
        out .*= rates*amount
        out
    end

    function darlington(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes ## return type illegal here
        println("Darlington ($ccode)")
        # allowance £246.47 a week plus fee until two weeks after a young person completes their 13th year of education.
        # Afterwards it reduces to the normal rate of £160.00 a week plus £20.00 from the young person.
        # since 13th yearis 18, assume 160
        allow = uprate(160.00,year,AFC_SURVEY_YEAR)
        contributions = uprate(20.00,year,AFC_SURVEY_YEAR)
        # bands A,A+topup,B,C
        fees = uprate.([0, 75.0,90,100,150],year,AFC_SURVEY_YEAR)
        fee = fees[carer.skill_level]
        CarerOutcomes(
            fee          = fee,
            allowance    = allow,
            housing_cont = 0.0,
            other_cont   = contributions )
        # CarerOutcomes( payments, contributions, 0.0, max(0.0,payments+contributions ), fee )
    end

    function sefton(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        hcontrib = uprate( 50.0, year, AFC_SURVEY_YEAR )
        if yp.employment_status == NEET
            outcomes.housing_cont = min( hcontrib, yp.housing_costs )
            outcomes.other_cont = uprate( 12.0, year, AFC_SURVEY_YEAR )
        end
        outcomes
    end

    function oldham_style_fee(
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        fees         :: Vector ) :: Real
        ufees = uprate.( fees, year, AFC_SURVEY_YEAR  )
        fee = ufees[carer.skill_level]
        if yp.age == 19
            fee /= 2
        elseif yp.age == 20
            fee = 0.0
        end
        fee
    end

    function random_oldham_fee(
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        fees         :: Vector ) :: Real
        fee = oldham_style_fee( year, carer, yp, fees )
        prop = 0.86 # see AFC note
        if rand() > prop
            fee = 0.0
        end
        fee
    end

    """
    Year one –Full payment allowance and fee
    Year 2- Full allowance plus half the skill payment
    Year 3- Full allowance
    """
    function oldham(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        fees = uprate.( [0, 78.38, 158.76, 237.23, 340.76], AFC_SURVEY_YEAR  )
        fee = fees[carer.skill_level]
        if yp.age == 19
            fee /= 2
        elseif yp.age == 20
            fee = 0.0
        end
        outcomes.fee = fee
        outcomes
    end

    function rochdale(
            ccode        :: AbstractString,
            year         :: Integer,
            carer        :: Carer,
            yp           :: YP,
            council_data :: DataFrameRow,
            params       :: Params ) :: CarerOutcomes

        # This sounds mich the same as Sefron, so ..
        outcomes = sefton( ccode, year, carer, yp, council_data, params )
        # The ‘Staying Put’ Allowance therefore is based on the former foster placement rate, which would be
        # applied on the young person’s 18th birthday minus the pocket money and clothing allowance element
        # as this should be replaced by the young person’s welfare benefit claim. The ‘Staying Put’ allowance is funded
        # from a number of sources including Housing and other Benefits, Young Person’s Contribution,
        # Personalised Budget and Local Authority. The ‘Staying Put’ Carer will receive the
        # boarding out rate minus pocket money and clothing element on a weekly basis. They will not receive contributions
        # in respect of Christmas or Festival payments, Birthday payment or Holiday payment in order to ensure equity with other care leavers.
        # so.. assume 1/2 comes from other sources GKS
        # outcomes.payments_from_la /= 2.0 #
        outcomes
    end

    function tameside(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        # Staying Put is paid at the same rate as supported
        # lodging which is £177.76 plus £22 contribution from the
        # young person Total £199.76
        outcomes.other_cont = uprate( 22.0, year,  AFC_SURVEY_YEAR )
        outcomes
    end

    function wirral(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        paylim = uprate( 50.0, year,  AFC_SURVEY_YEAR )
        outcomes.other_cont = max(0.0, 0.5*(yp.earnings - paylim))
        outcomes
    end

    function derby(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        outcomes.other_cont = uprate( 12.0, year, AFC_SURVEY_YEAR )
        # 241.72 less YP & HB payments
        if yp.employment_status == NEET
            outcomes.housing_cont = yp.housing_costs
        end
        outcomes
    end

    """
    Staying Put is based on the rate we are paying when a young person reaches
    18 with a deduction of £23 which the young person is expected to
    contribute and any housing benefit that is payable.
    """
    function derbyshire(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        outcomes.other_cont = uprate( 23.0, year, AFC_SURVEY_YEAR )
        if yp.employment_status == NEET
            outcomes.housing_cont = yp.housing_costs
        end
        outcomes
    end

    """
    The weekly allowance is £80.00 Per week with an additional £86.30 Per week (£166.30) housing
    benefit depending on circumstances the £86.30  may be paid by the Local Authority.
    (essentially if a young person cannot claim housing benefit). This payment remains the same for any period that a Young person is on the Staying Put scheme
    """
    function leicester(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        if yp.employment_status == NEET
            outcomes.housing_cont = yp.housing_costs
        end
        outcomes
    end

    """
        We do not operate a fixed payment level / band. Calculations are based on
        what the foster carers would have been getting as a carer minus housing
        benefit for their local area and a contribution by young person (£30),
    """
    function lincolnshire(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        outcomes.other_cont = uprate( 30.0, year, AFC_SURVEY_YEAR )
        if yp.employment_status == NEET
            outcomes.housing_cont = yp.housing_costs
        end
        outcomes
    end

    """
     £118.00 per week. £210 (£118.00 + £92.00) if the young person does not receive housing benefits.
    """
    function coventry(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        hcontrib = uprate( 92.0, year, AFC_SURVEY_YEAR )
        if yp.employment_status == NEET
            outcomes.housing_cont = min( hcontrib, yp.housing_costs )
        end
        outcomes
    end

    """
    Basic 195.00. Hybrid 245.00. High Needs 295.00.
    """
    function shropshire(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        levels = uprate.([195.0, 245.0, 295.0 ], year, AFC_SURVEY_YEAR )
        allowance = assignrand( [0.333,0.666,1.0], levels )
        CarerOutcomes(
            fee          = 0.0,
            allowance    = allowance,
            housing_cont = 0.0,
            other_cont   = 0.0 )
    end

    """
    216.42. Comprising: Contribution from the Young Person's wages/JSA. 25.13.
    100% of the Young Person's housing benefit entitlement (variable) with the balance met by Children's Social Services.
    191.29

    """
    function solihull(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        othcontrib = uprate( 25.13, year, AFC_SURVEY_YEAR )
        if yp.employment_status == NEET
            outcomes.housing_cont = yp.housing_costs
        end
        outcomes.other_cont = othcontrib
        outcomes
    end

    """
    100 + c. 60 housing benefit
    """
    function stokeontrent(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        hcontrib = uprate( 60.0, year, AFC_SURVEY_YEAR )
        if yp.employment_status == NEET
            outcomes.housing_cont = min( hcontrib, yp.housing_costs )
        end
        outcomes
    end

    """
    Staying Put is a fixed rate of £232 per week, and does not change each year.
    The young person should apply for Housing Benefit (unless the carer is a relative or
    claiming Housing Benefit in their own right).  The amount awarded to the yp in Housing Benefit is then deducted from the £232.00.
    """
    function walsall(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        if yp.employment_status == NEET
            outcomes.housing_cont = yp.housing_costs
        end
        outcomes
    end

    """
     Warwickshire County council pays an allowance of £180.00 per week for 'Staying Put' carers,
     but the Young Person pays a contribution if they are working.
    """
    function warwickshire(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        if yp.earnings > 0.0
            outcomes.other_cont = yp.earnings * 0.5 # FIXME completely made up number
        end
        outcomes
    end

    """
    Payments are static irrespective of year one, two or three. Maximum care payment is either £255.00 or £331.50 depending on level of Foster Carer.
    It is broken down into three components totalling a maximum care payment:
    1. Local Housing Allowance for accommodating Authority e.g. Wolverhampton is £86.30
    2. Contribution from young person - £25.00.
    3. Care Payment from Wolverhampton to total up to maximum care payment.
    """
    function wolverhampton(
            ccode        :: AbstractString,
            year         :: Integer,
            carer        :: Carer,
            yp           :: YP,
            council_data :: DataFrameRow,
            params       :: Params ) :: CarerOutcomes
            outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
            if carer.skill_level > 3
                outcomes.allowance = uprate( 331.50, year, AFC_SURVEY_YEAR )
            end
            othcontrib = uprate( 25.00, year, AFC_SURVEY_YEAR )
            if yp.employment_status == NEET
                outcomes.housing_cont = yp.housing_costs
            end
            outcomes.other_cont = othcontrib
            outcomes
    end

    """
    £137.50 per week if young person is receiving housing benefit and £197.50 if young
    person is not receiving housing benefit.
    This is a flat weekly rate during the staying put agreement.
    """
    function worcestershire(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        hcontrib = uprate( 60.0, year, AFC_SURVEY_YEAR )
        if yp.employment_status == NEET
            outcomes.housing_cont = min( hcontrib, yp.housing_costs )
        end
        outcomes
    end

    """
    We are supporting five  tier 1, one tier 2, and three tier 3.  we also have 2 with no tier group.#
    For the first 4 weeks we pay the full fostering  payment while the client is arranging benefits.
    After 4 weeks the staying put payments are made up of housing benefit approx. £111.00 per week,
    The carer receives their Tier payment, The  carer receives £25.00 for food & £20.00 utilities.
    FIXME is this what that means?
    """
    function luton(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        if yp.employment_status == NEET
            outcomes.housing_cont = yp.housing_costs
        end
        outcomes
    end

    """
        # Standard Allowance£271.08
        # Pocket Money-£12.00
        # Clothing Allowance-£19.39
        # Savings-£20.00
        # TOTAL£219.69
        FIXME WTF
    """
    function bedford(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        outcomes
    end

    function bracknellforest(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        payments = uprate.([160.0, 229, 298, 367, 506.0 ], year, AFC_SURVEY_YEAR )
        allowance = payments[carer.skill_level]
        CarerOutcomes(
            fee=0,
            allowance=allowance,
            housing_cont = 0,
            other_cont= 0 )
    end

    """
      16+ fostering allowance  = £242.08 per week - £34.41.(PA) = £207.67, Housing Benefit = £70 per week (for example), CSD contribution is £207.67 - £70 = £137.67 per week.
    """
    function hampshire(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        othcontrib = uprate( 34.41, year, AFC_SURVEY_YEAR )
        if yp.employment_status == NEET
            outcomes.housing_cont = yp.housing_costs
        end
        outcomes.other_cont = othcontrib
        outcomes
    end

    function medway(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        allowances = uprate.([142.45, 194.11, 226.52 ], year, AFC_SURVEY_YEAR )
        allowance = assignrand( [0.333,0.666,1.0], allowances )
        CarerOutcomes(
            fee=0,
            allowance=allowance,
            housing_cont = 0,
            other_cont= 0 )
    end

    """
    Each carer will receive £200 per week (Staying Put Rate -2017-18) Some may receive £180 per week and the young person makes a £20 contribution towards this weekly.
    """
    function swindon(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        othcontrib = uprate( 20.0, year, AFC_SURVEY_YEAR )
        outcomes.other_cont = othcontrib
        outcomes
    end

    """
      16+ fostering allowance  = £242.08 per week - £34.41.(PA) = £207.67, Housing Benefit = £70 per week (for example), CSD contribution is £207.67 - £70 = £137.67 per week.
    """
    function wiltshire(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        othcontrib = uprate( 20.0, year, AFC_SURVEY_YEAR )
        outcomes.other_cont = othcontrib
        outcomes
    end

    ## FIXME random band fee took middle
    FUNCS["XXHerefordshire"]=function( cohort :: Vector, data :: DataFrameRow, params :: Params )
        frame = do_basic_calc( "Herefordshire", cohort, data, params )
        # £118.00 per week. £210 (£118.00 + £92.00) if the young person does not receive housing benefits.
    end

    function westsussex(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        othcontrib = uprate( 10.0, year, AFC_SURVEY_YEAR )
        outcomes.other_cont = othcontrib
        outcomes
    end

    function bathandnesomerset(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        othcontrib = uprate( 13.0, year, AFC_SURVEY_YEAR )
        outcomes.other_cont = othcontrib
        outcomes
    end

    ## FIXME come back & see if I can make some sense out of Wigan
    FUNCS["XXWigan"]=function( cohort :: Vector, data :: DataFrameRow, params :: Params )
        frame = do_basic_calc( "Wigan", cohort, data, params )
        # Year 1. We have 2 rates within the first year 1 – Standard Staying Put Rate of £170.00
        # per week 2 – Providers are paid their current Fostering rate until the end of the Academic Year – these fees range between £307 and £586 per week dependant on what Tier they are on as a foster carer.
        # Years 2 & 3. We pay all Staying Put arrangements standard fee of £170.00 per week.

        fostered_by_age = round.( Integer, cohort )
        num_fostered = sum( fostered_by_age )
        frame = makeframe( "Wigan", data[:ccode], data[:rcode], num_fostered )
        y_1 = fostered_by_age[1]
        after_y_1 = sum( fostered_by_age[2:params.num_years] )
        for i in (y_1+1):params.num_years
            frame.payments[i] = 170.0;
        end
        frame
    end

    """
    #support costs for Project carers if high support needs over the age of 18#
    #£66.70 per week rent - payable by Housing Benefit (or Staying Put budget
    # if HB claim not possible)
    # £30.00 per week contribution from the young person for food and utilities. 17 on standard level and 1 on enhanced
    """
    function dorset(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        hcontrib = uprate( 66.70, year, AFC_SURVEY_YEAR )
        othcontrib = uprate( 30.0, year, AFC_SURVEY_YEAR )
        if yp.employment_status == NEET
            outcomes.housing_cont = min( hcontrib, yp.housing_costs )
            outcomes.other_cont = outcomes.housing_cont
        end
        outcomes
    end

    function somerset(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        if yp.employment_status == NEET
            outcomes.other_cont = yp.benefits
        end
        outcomes
    end

    function byage(
        yp           :: YP,
        payments     :: Vector,
        year         :: Integer
         ) :: CarerOutcomes
         payments = uprate.(payments, year, AFC_SURVEY_YEAR )
         i = yp.age - 17
         allowance = payments[i]
         CarerOutcomes(
            fee=0,allowance=allowance,housing_cont =0, other_cont=0 )
    end

    """
    Year 1 (18-19 yrs) - £250
    Year 2 (19-20 yrs) - £230
    Year 3 (20-21 yrs) - £180
    """
    function cambridgeshire(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        allowances = uprate.( [250, 230, 180 ], year, AFC_SURVEY_YEAR )
        byage( yp, allowances, year )
    end

    """
    For carers accredited at Levels 1 to 4, the carer’s total payment in respect
    of the placement will remain in line with their accredited status (i.e. as if
    they were caring for a looked after child aged 17).
    Payments to carers accredited at Level 5 will be capped at the Level 4 rate.
    NCC will pay the carer the appropriate accreditation level minus the rent
    liability for the young person (which is set at the one bedroom
    self-contained Local Housing Allowance) and £10pw utilities.
    100pw retainer while at uni
    """
    function norfolk(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes

        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        othcontrib = uprate( 10.0, year, AFC_SURVEY_YEAR )
        if yp.employment_status == NEET
            outcomes.housing_cont = yp.housing_costs
        end
        outcomes.other_cont = othcontrib
        outcomes
    end

    """
    The Staying Put Policy matches the current arrangements for foster carers
    for children aged 16+, including the skills element but minus pocket
    money, clothing and holiday (£303.10). Staying Put providers receive the
    full skills payment in the first year of a Staying Put arrangement, 75% in
    the second year and 50% in year three.
    """
    function buckinghamshire(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        othcontrib = uprate( 10.0, year, AFC_SURVEY_YEAR )
        outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
        if yp.age == 19
            outcomes.allowance *= 0.75
        elseif yp.age == 20
            outcomes.allowance *= 0.5
        end
        outcomes
    end

    """
    No allowances are paid directly to the former foster carer, a subsistence
    allowance is paid directly to the young person

    Former foster carers are paid on the relevant Band fee that they were
    previously on or remain on as foster carers,
    (see table in response to Q3) there is no differentiation in subsequent years,
    and they continue to be paid at the same rate throughout the Staying Put
    arrangement. 7 carers.
    """
    function ealing(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
        fee = oldham_style_fee( year, carer, yp, [0, 78.38, 158.76, 237.23, 340.76] )
        CarerOutcomes( fee=fee, allowance=0.0, housing_cont =0, other_cont=0)
    end

    """
    The Carer would receive the amount/tier level as per the attached that they were on previous to becoming Staying Put minus the below income:-
    Income support.
    Housing Benefit.
    Wages.
    £11.00 utility bills.
    The Staying Put Policy is completed over 2 Stages not in 3 year stages.  Stage 1 they receive allowance and fee.Stage 2 minus the fee.
    """
    function suffolk(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes
    outcomes = do_basic_calc( ccode, year, carer, yp, council_data, params )
    othcontrib = uprate( 11.0, year, AFC_SURVEY_YEAR )
    if yp.employment_status == NEET
        outcomes.housing_cont = yp.housing_costs
        outcomes.other_cont = yp.benefits
    end
    outcomes
end

    # TODO simple assign by age West Beds 248.67 198.67 148.67 ??
    # TODO Great Matrix Suffolk
    # TODO Brighton
    # TODO Kent
    # TODO Bristol
    # TODO Bournemouth
    # TODO North Somerset
    # TODO Plymouth
    # TODO POOLE academic year
    # TODO Somerset Fee structure

    FUNCS["E06000005"]=darlington
    FUNCS["E08000005"]=rochdale
    FUNCS["E08000014"]=sefton
    FUNCS["E08000008"]=tameside
    FUNCS["E08000015"]=wirral
    FUNCS["E06000015"]=derby
    FUNCS["E10000007"]=derbyshire
    FUNCS["E06000016"]=leicester
    FUNCS["E10000019"]=lincolnshire
    FUNCS["E08000026"]=coventry
    FUNCS["E06000051"]=shropshire
    FUNCS["E08000029"]=solihull
    FUNCS["E06000021"]=stokeontrent
    FUNCS["E08000030"]=walsall
    FUNCS["E10000031"]=warwickshire
    FUNCS["E08000031"]=wolverhampton
    FUNCS["E10000034"]=worcestershire
    FUNCS["E06000032"]=luton
    FUNCS["E06000055"]=bedford
    FUNCS["E06000036"]=bracknellforest
    FUNCS["E10000014"]=hampshire
    FUNCS["E06000035"]=medway
    FUNCS["E06000030"]=swindon
    FUNCS["E06000054"]=wiltshire
    FUNCS["E10000032"]=westsussex
    FUNCS["E06000022"]=bathandnesomerset
    FUNCS["E10000009"]=dorset
    FUNCS["E10000027"]=somerset
    FUNCS["E10000003"]=cambridgeshire
    FUNCS["E08000004"]=oldham
    FUNCS["E10000020"]=norfolk
    FUNCS["E10000002"]=buckinghamshire

    FUNCS["E10000029"]= suffolk
    FUNCS["E09000009"]= ealing

    # @enum ContributionType actual_contribution no_contribution all_people benefits_only flat_rate
    # @enum PaymentType actual_payment min_payment age_16_17
    # @enum UpratingType no_uprating by_benefits by_cpi

    # now override
    function override_outcomes!(
        outcomes     :: CarerOutcomes,
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes

        ## FIXME these two things overlap horribly
        if params.yp_contrib_type == no_contribution
                outcomes.housing_cont = 0.0
                outcomes.other_cont = 0.0
        elseif params.yp_contrib_type == flat_rate
                outcomes.other_cont = params.yp_contrib
                outcomes.housing_cont = 0.0
        elseif params.yp_contrib_type == benefits_only
                if yp.benefits == 0.0
                    outcomes.other_cont = 0.0
                    outcomes.housing_cont = 0.0
                end
        elseif params.yp_contrib_type == hb_only
            if yp.benefits == 0.0
                outcomes.other_cont = 0.0
                # outcomes.housing_cont = 0.0
            end

        end

        if params.contrib_hb == no_contribution
            outcomes.housing_cont = 0.0
        elseif params.contrib_hb == all_people
            outcomes.housing_cont = yp.housing_costs
        elseif params.contrib_hb == benefits_only
            if yp.benefits != 0.0
                outcomes.housing_cont = yp.housing_costs
            end
        end
        if params.payment == min_payment
            reg = get_pay_class_from_region_code( council_data.rcode )
            minp = getAllowances( reg, year ).age_16_17
            outcomes.allowance = max( outcomes.allowance, minp )
        end
        fsize = size(params.fee)[1]
        if fsize > 0
            @assert fsize == 5
            fee = params.fee[carer.skill_level]
            outcomes.fee = max( fee, outcomes.fee )
        end
        tsize = size(params.taper)[1]
        if tsize > 0
            @assert tsize == 3
            p = yp.age - 17
            outcomes.allowance *= params.taper[p]
            outcomes.fee *= params.taper[p]
            outcomes.housing_cont *= params.taper[p]
            outcomes.other_cont *= params.taper[p]
        end
        delete_fee = rand() < params.prop_fees_deleted
        if( delete_fee )
            outcomes.fee = 0
        end
        outcomes
    end

    function do_one_calc(
        ccode        :: AbstractString,
        year         :: Integer,
        carer        :: Carer,
        yp           :: YP,
        council_data :: DataFrameRow,
        params       :: Params ) :: CarerOutcomes# now override

        if ccode in keys( FUNCS )
            outcomes = FUNCS[ccode](
                ccode,
                year,
                carer,
                yp,
                council_data,
                params )
        else
            outcomes = do_basic_calc(
                ccode,
                year,
                carer,
                yp,
                council_data,
                params )
        end
        override_outcomes!(
            outcomes,
            ccode,
            year,
            carer,
            yp,
            council_data,
            params )
        fillcomponents!( outcomes )
        outcomes
    end # do_one_calc

    """
    Given some number `base` and a series of (e.g.) populations `series`,
    return a vector which starts at base and moves in proportion to `series`
    useful for (e.g.) forecasting some grant allocation. the `base_period` is the 1st element of
    series to track from, base year allows the base to represent some period
    in the past, relative to AFC_SURVEY_YEAR

    so: base = 100 series=[1,2,3] => 100,200,300

    """
    function track_series(
        base    :: Real,
        series  :: Vector,
        base_period :: Integer = 1,
        base_year   :: Integer = AFC_SURVEY_YEAR
        ) :: Vector

        b = series[base_period]
        y = base_year
        n = size( series )[1]
        out = zeros(n)
        for i in 1:n
            x = base * series[i]/b
            out[i]=uprate( x, y, base_year )
            y += 1
        end
        out
    end

end # module Staying Put Sim


text:
# This file was generated, do not modify it. # hide
tuned_forest = TunedModel(model=forest,
                          tuning=Grid(resolution=3),
                          resampling=CV(nfolds=6, rng=StableRNG(32)),
                          ranges=[r_sf, r_bf],
                          measure=rms)
m = machine(tuned_forest, X, y)
e = evaluate!(m, resampling=Holdout(fraction_train=0.8),
              measure=[rms, rmslp1])
e |> pprint


text:
function set_up_operator(p::Param)

    # Variables
    rmin        = p.rmin
    rmax        = p.rmax
    rnodes      = p.rnodes

    thnodes     = p.thnodes

    #Create position vectors
    dr = (rmax - rmin)/(rnodes-1)
    r = [rmin + (i-1)*dr for i in 1:rnodes]

    dth = 2*pi/(thnodes)
    theta = [ (i-1)*dth for i in 1:thnodes]

    #Create Differential Operators
    Dr  =   Diff_Operator_2D(1 ,1, [dr,dth] , [rnodes,thnodes])
    Drr =   Diff_Operator_2D(2 ,1, [dr,dth] , [rnodes,thnodes])

    Dp  =  Diff_Operator_2D_Periodic(1 ,2, [dr,dth] , [rnodes,thnodes])
    Dpp =  Diff_Operator_2D_Periodic(2 ,2, [dr,dth] , [rnodes,thnodes])

    one_over_r2 = spzeros(rnodes*thnodes,rnodes*thnodes)

    for i in 1:rnodes*thnodes
        one_over_r2[i,i] = (1/r[1 + (i-1)%rnodes])^2
    end

    RHS_Opp = Drr + one_over_r2*Dpp + (0.25)*one_over_r2
    #sparse(Drr + (0.25)*one_over_r2*I  + one_over_r2*Dpp);

    av1 = zeros(rnodes*thnodes)
    av2 = zeros(rnodes*thnodes)

    #Operator structure
    return Opps(r, theta , av1 , av2  , Dr, Dp, RHS_Opp)

end


text:
#==============================================================================#
# Point-like particle
#==============================================================================#

"""
	mutable struct Point <: AbstractSpecies

Species corresponding to a single point-like particle composed of a single monomer type.
"""
mutable struct Point <: AbstractSpecies
	mids     :: Vector{Int}
	Nref     :: Float64

	# SCF fields
	Q        :: Float64
	density  :: Dict{Int,FieldGrid{Float64}}
	system   :: Option{FieldSystem}

	function Point(mon::Monomer)
		return new([mon.id], mon.vol, 0.0, Dict(), nothing)
	end
end

#==============================================================================#

Base.show(io::IO, point::Point) = @printf(io, "Point(mid = %d)", point.mids[1])

monomer_fraction(point::Point, mid::Integer) = hasmonomer(point, mid) ? 1.0 : 0.0

function setup!(point::Point, sys::FieldSystem)
	point.system = sys
	point.Q = 0.0
	point.density[point.mids[1]] = zeros(Float64, sys.dims)
end

function density!(point::Point)
	@assert !isnothing(point.system)
	sys = point.system
	mon = sys.monomers[point.mids[1]]
	omega = sys.fields[mon.id]

	# Partition function via Boltzmann weight
	Q = 0.0
	@simd for i in eachindex(omega)
		@inbounds Q += exp(-mon.vol * omega[i])
	end
	point.Q = Q / ngrid(sys) # Normalize by number of grid points

	# Update point density field
	@. point.density[mon.id] = exp(-mon.vol * omega) / point.Q

	return nothing
end

function scfstress(point::Point)
	@assert !isnothing(point.system)
	return zeros(nparams(point.system.cell))
end


text:
## sparse matrix multiplication

function (*){TvA,TiA,TvB,TiB}(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB})
    Tv = promote_type(TvA, TvB)
    Ti = promote_type(TiA, TiB)
    A  = convert(SparseMatrixCSC{Tv,Ti}, A)
    B  = convert(SparseMatrixCSC{Tv,Ti}, B)
    A * B
end

(*){TvA,TiA}(A::SparseMatrixCSC{TvA,TiA}, X::BitArray{1}) = invoke(*, (SparseMatrixCSC, AbstractVector), A, X)
# In matrix-vector multiplication, the correct orientation of the vector is assumed.
function A_mul_B!(α::Number, A::SparseMatrixCSC, x::AbstractVector, β::Number, y::AbstractVector)
    A.n == length(x) || throw(DimensionMismatch(""))
    A.m == length(y) || throw(DimensionMismatch(""))
    if β != 1
        β != 0 ? scale!(y,β) : fill!(y,zero(eltype(y)))
    end
    nzv = A.nzval
    rv = A.rowval
    for col = 1 : A.n
        αx = α*x[col]
        @inbounds for k = A.colptr[col] : (A.colptr[col+1]-1)
            y[rv[k]] += nzv[k]*αx
        end
    end
    y
end
A_mul_B!(y::AbstractVector, A::SparseMatrixCSC, x::AbstractVector) = A_mul_B!(one(eltype(x)), A, x, zero(eltype(y)), y)

function *{TA,S,Tx}(A::SparseMatrixCSC{TA,S}, x::AbstractVector{Tx})
    T = promote_type(TA,Tx)
    A_mul_B!(one(T), A, x, zero(T), Array(T, A.m))
end

function Ac_mul_B!(α::Number, A::SparseMatrixCSC, x::AbstractVector, β::Number, y::AbstractVector)
    A.n == length(y) || throw(DimensionMismatch(""))
    A.m == length(x) || throw(DimensionMismatch(""))
    nzv = A.nzval
    rv = A.rowval
    zro = zero(eltype(y))
    @inbounds begin
        for i = 1 : A.n
            if β != 1
                if β != 0
                    y[i] *= β
                else
                    y[i] = zro
                end
            end
            tmp = zero(eltype(y))
            for j = A.colptr[i] : (A.colptr[i+1]-1)
                tmp += conj(nzv[j])*x[rv[j]]
            end
            y[i] += α*tmp
        end
    end
    y
end
Ac_mul_B!(y::AbstractVector, A::SparseMatrixCSC, x::AbstractVector) = Ac_mul_B!(one(eltype(x)), A, x, zero(eltype(y)), y)
function At_mul_B!(α::Number, A::SparseMatrixCSC, x::AbstractVector, β::Number, y::AbstractVector)
    A.n == length(y) || throw(DimensionMismatch(""))
    A.m == length(x) || throw(DimensionMismatch(""))
    nzv = A.nzval
    rv = A.rowval
    zro = zero(eltype(y))
    @inbounds begin
        for i = 1 : A.n
            if β != 1
                if β != 0
                    y[i] *= β
                else
                    y[i] = zro
                end
            end
            tmp = zero(eltype(y))
            for j = A.colptr[i] : (A.colptr[i+1]-1)
                tmp += nzv[j]*x[rv[j]]
            end
            y[i] += α*tmp
        end
    end
    y
end
At_mul_B!(y::AbstractVector, A::SparseMatrixCSC, x::AbstractVector) = At_mul_B!(one(eltype(x)), A, x, zero(eltype(y)), y)
function Ac_mul_B{TA,S,Tx}(A::SparseMatrixCSC{TA,S}, x::AbstractVector{Tx})
    T = promote_type(TA, Tx)
    Ac_mul_B!(one(T), A, x, zero(T), Array(T, A.n))
end
function At_mul_B{TA,S,Tx}(A::SparseMatrixCSC{TA,S}, x::AbstractVector{Tx})
    T = promote_type(TA, Tx)
    At_mul_B!(one(T), A, x, zero(T), Array(T, A.n))
end

*(X::BitArray{1}, A::SparseMatrixCSC) = invoke(*, (AbstractVector, SparseMatrixCSC), X, A)
# In vector-matrix multiplication, the correct orientation of the vector is assumed.
# XXX: this is wrong (i.e. not what Arrays would do)!!
function *{T1,T2}(X::AbstractVector{T1}, A::SparseMatrixCSC{T2})
    A.m==length(X) || throw(DimensionMismatch(""))
    Y = zeros(promote_type(T1,T2), A.n)
    nzv = A.nzval
    rv = A.rowval
    for col =1:A.n, k=A.colptr[col]:(A.colptr[col+1]-1)
        Y[col] += X[rv[k]] * nzv[k]
    end
    Y
end

*{TvA,TiA}(A::SparseMatrixCSC{TvA,TiA}, X::BitArray{2}) = invoke(*, (SparseMatrixCSC, AbstractMatrix), A, X)
function (*){TvA,TiA,TX}(A::SparseMatrixCSC{TvA,TiA}, X::AbstractMatrix{TX})
    mX, nX = size(X)
    A.n==mX || throw(DimensionMismatch(""))
    Y = zeros(promote_type(TvA,TX), A.m, nX)
    nzv = A.nzval
    rv = A.rowval
    colptr = A.colptr
    for multivec_col=1:nX, col=1:A.n
        Xc = X[col, multivec_col]
        @inbounds for k = colptr[col] : (colptr[col+1]-1)
           Y[rv[k], multivec_col] += nzv[k] * Xc
        end
    end
    Y
end

*{TvA,TiA}(X::BitArray{2}, A::SparseMatrixCSC{TvA,TiA}) = invoke(*, (AbstractMatrix, SparseMatrixCSC), X, A)
# TODO: Tridiagonal * Sparse should be implemented more efficiently
*{TX,TvA,TiA}(X::Tridiagonal{TX}, A::SparseMatrixCSC{TvA,TiA}) = invoke(*, (Tridiagonal, AbstractMatrix), X, A)
*{TvA,TiA}(X::Triangular, A::SparseMatrixCSC{TvA,TiA}) = full(X)*A
function *{TX,TvA,TiA}(X::AbstractMatrix{TX}, A::SparseMatrixCSC{TvA,TiA})
    mX, nX = size(X)
    nX == A.m || throw(DimensionMismatch(""))
    Y = zeros(promote_type(TX,TvA), mX, A.n)
    for multivec_row=1:mX, col=1:A.n, k=A.colptr[col]:(A.colptr[col+1]-1)
        Y[multivec_row, col] += X[multivec_row, A.rowval[k]] * A.nzval[k]
    end
    Y
end

# Sparse matrix multiplication as described in [Gustavson, 1978]:
# http://www.cse.iitb.ac.in/graphics/~anand/website/include/papers/matrix/fast_matrix_mul.pdf
function *{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti})
    mA, nA = size(A)
    mB, nB = size(B)
    nA==mB || throw(DimensionMismatch(""))

    colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval
    colptrB = B.colptr; rowvalB = B.rowval; nzvalB = B.nzval
    # TODO: Need better estimation of result space
    nnzC = min(mA*nB, length(nzvalA) + length(nzvalB))
    colptrC = Array(Ti, nB+1)
    rowvalC = Array(Ti, nnzC)
    nzvalC = Array(Tv, nnzC)

    @inbounds begin
        ip = 1
        xb = zeros(Ti, mA)
        x  = zeros(Tv, mA)
        for i in 1:nB
            if ip + mA - 1 > nnzC
                resize!(rowvalC, nnzC + max(nnzC,mA))
                resize!(nzvalC, nnzC + max(nnzC,mA))
                nnzC = length(nzvalC)
            end
            colptrC[i] = ip
            for jp in colptrB[i]:(colptrB[i+1] - 1)
                nzB = nzvalB[jp]
                j = rowvalB[jp]
                for kp in colptrA[j]:(colptrA[j+1] - 1)
                    nzC = nzvalA[kp] * nzB
                    k = rowvalA[kp]
                    if xb[k] != i
                        rowvalC[ip] = k
                        ip += 1
                        xb[k] = i
                        x[k] = nzC
                    else
                        x[k] += nzC
                    end
                end
            end
            for vp in colptrC[i]:(ip - 1)
                nzvalC[vp] = x[rowvalC[vp]]
            end
        end
        colptrC[nB+1] = ip
    end

    deleteat!(rowvalC, colptrC[end]:length(rowvalC))
    deleteat!(nzvalC, colptrC[end]:length(nzvalC))

    # The Gustavson algorithm does not guarantee the product to have sorted row indices.
    Cunsorted = SparseMatrixCSC(mA, nB, colptrC, rowvalC, nzvalC)
    Ct = Cunsorted.'
    Ctt = Base.SparseMatrix.transpose!(Ct, SparseMatrixCSC(mA, nB, colptrC, rowvalC, nzvalC))
end

## solvers
function A_ldiv_B!(A::SparseMatrixCSC, b::AbstractVecOrMat)
    if iseltype(b, Complex); A = complex(A); end

    if istril(A)
        # TODO: Fix diagonal case. Diagonal(A.nzval) needs to handle
        # the case where there are zeros on the diagonal and error out.
        # It also does not work in the complex case. VBS.
        #if istriu(A); return A_ldiv_B!(Diagonal(A.nzval), b); end
        return fwdTriSolve!(A, b)
    end
    if istriu(A); return bwdTriSolve!(A, b); end
    return A_ldiv_B!(lufact(A),b)
end

function fwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat)
# forward substitution for CSC matrices
    n = length(B)
    if isa(B, Vector)
        nrowB = n
        ncolB = 1
    else
        nrowB, ncolB = size(B)
    end
    ncol = chksquare(A)
    if nrowB != ncol
        throw(DimensionMismatch("A is $(ncol)X$(ncol) and B has length $(n)"))
    end

    aa = A.nzval
    ja = A.rowval
    ia = A.colptr

    joff = 0
    for k = 1:ncolB
        for j = 1:(nrowB-1)
            jb = joff + j
            i1 = ia[j]
            i2 = ia[j+1]-1
            B[jb] /= aa[i1]
            bj = B[jb]
            for i = i1+1:i2
                B[joff+ja[i]] -= bj*aa[i]
            end
        end
        joff += nrowB
        B[joff] /= aa[end]
    end
    return B
end

function bwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat)
# backward substitution for CSC matrices
    n = length(B)
    if isa(B, Vector)
        nrowB = n
        ncolB = 1
    else
        nrowB, ncolB = size(B)
    end
    ncol = chksquare(A)
    if nrowB != ncol throw(DimensionMismatch("A is $(ncol)X$(ncol) and B has length $(n)")) end

    aa = A.nzval
    ja = A.rowval
    ia = A.colptr

    joff = 0
    for k = 1:ncolB
        for j = nrowB:-1:2
            jb = joff + j
            i1 = ia[j]
            i2 = ia[j+1]-1
            B[jb] /= aa[i2]
            bj = B[jb]
            for i = i2-1:-1:i1
                B[joff+ja[i]] -= bj*aa[i]
            end
        end
        B[joff+1] /= aa[1]
        joff += nrowB
    end
   return B
end

## triu, tril

function triu{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, k::Integer)
    m,n = size(S)
    colptr = Array(Ti, n+1)
    nnz = 0
    for col = 1 : min(max(k+1,1), n+1)
        colptr[col] = 1
    end
    for col = max(k+1,1) : n
        for c1 = S.colptr[col] : S.colptr[col+1]-1
            S.rowval[c1] > col - k && break
            nnz += 1
        end
        colptr[col+1] = nnz+1
    end
    rowval = Array(Ti, nnz)
    nzval = Array(Tv, nnz)
    A = SparseMatrixCSC{Tv,Ti}(m, n, colptr, rowval, nzval)
    for col = max(k+1,1) : n
        c1 = S.colptr[col]
        for c2 = A.colptr[col] : A.colptr[col+1]-1
            A.rowval[c2] = S.rowval[c1]
            A.nzval[c2] = S.nzval[c1]
            c1 += 1
        end
    end
    A
end

function tril{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, k::Integer)
    m,n = size(S)
    colptr = Array(Ti, n+1)
    nnz = 0
    colptr[1] = 1
    for col = 1 : min(n, m+k)
        l1 = S.colptr[col+1]-1
        for c1 = 0 : (l1 - S.colptr[col])
            S.rowval[l1 - c1] < col - k && break
            nnz += 1
        end
        colptr[col+1] = nnz+1
    end
    for col = max(min(n, m+k)+2,1) : n+1
        colptr[col] = nnz+1
    end
    rowval = Array(Ti, nnz)
    nzval = Array(Tv, nnz)
    A = SparseMatrixCSC{Tv,Ti}(m, n, colptr, rowval, nzval)
    for col = 1 : min(n, m+k)
        c1 = S.colptr[col+1]-1
        l2 = A.colptr[col+1]-1
        for c2 = 0 : l2 - A.colptr[col]
            A.rowval[l2 - c2] = S.rowval[c1]
            A.nzval[l2 - c2] = S.nzval[c1]
            c1 -= 1
        end
    end
    A
end

## diff

function sparse_diff1{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti})
    m,n = size(S)
    m > 1 || return SparseMatrixCSC{Tv,Ti}(0, n, ones(n+1), Ti[], Tv[])
    colptr = Array(Ti, n+1)
    numnz = 2 * nnz(S) # upper bound; will shrink later
    rowval = Array(Ti, numnz)
    nzval = Array(Tv, numnz)
    numnz = 0
    colptr[1] = 1
    for col = 1 : n
        last_row = 0
        last_val = 0
        for k = S.colptr[col] : S.colptr[col+1]-1
            row = S.rowval[k]
            val = S.nzval[k]
            if row > 1
                if row == last_row + 1
                    nzval[numnz] += val
                    nzval[numnz]==zero(Tv) && (numnz -= 1)
                else
                    numnz += 1
                    rowval[numnz] = row - 1
                    nzval[numnz] = val
                end
            end
            if row < m
                numnz += 1
                rowval[numnz] = row
                nzval[numnz] = -val
            end
            last_row = row
            last_val = val
        end
        colptr[col+1] = numnz+1
    end
    deleteat!(rowval, numnz+1:length(rowval))
    deleteat!(nzval, numnz+1:length(nzval))
    return SparseMatrixCSC{Tv,Ti}(m-1, n, colptr, rowval, nzval)
end

function sparse_diff2{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti})

    m,n = size(a)
    colptr = Array(Ti, max(n,1))
    numnz = 2 * nnz(a) # upper bound; will shrink later
    rowval = Array(Ti, numnz)
    nzval = Array(Tv, numnz)

    z = zero(Tv)

    colptr_a = a.colptr
    rowval_a = a.rowval
    nzval_a = a.nzval

    ptrS = 1
    colptr[1] = 1

    n == 0 && return SparseMatrixCSC{Tv,Ti}(m, n, colptr, rowval, nzval)

    startA = colptr_a[1]
    stopA = colptr_a[2]

    rA = startA : stopA - 1
    rowvalA = rowval_a[rA]
    nzvalA = nzval_a[rA]
    lA = stopA - startA

    for col = 1:n-1
        startB, stopB = startA, stopA
        startA = colptr_a[col+1]
        stopA = colptr_a[col+2]

        rowvalB = rowvalA
        nzvalB = nzvalA
        lB = lA

        rA = startA : stopA - 1
        rowvalA = rowval_a[rA]
        nzvalA = nzval_a[rA]
        lA = stopA - startA

        ptrB = 1
        ptrA = 1

        while ptrA <= lA && ptrB <= lB
            rowA = rowvalA[ptrA]
            rowB = rowvalB[ptrB]
            if rowA < rowB
                rowval[ptrS] = rowA
                nzval[ptrS] = nzvalA[ptrA]
                ptrS += 1
                ptrA += 1
            elseif rowB < rowA
                rowval[ptrS] = rowB
                nzval[ptrS] = -nzvalB[ptrB]
                ptrS += 1
                ptrB += 1
            else
                res = nzvalA[ptrA] - nzvalB[ptrB]
                if res != z
                    rowval[ptrS] = rowA
                    nzval[ptrS] = res
                    ptrS += 1
                end
                ptrA += 1
                ptrB += 1
            end
        end

        while ptrA <= lA
            rowval[ptrS] = rowvalA[ptrA]
            nzval[ptrS] = nzvalA[ptrA]
            ptrS += 1
            ptrA += 1
        end

        while ptrB <= lB
            rowval[ptrS] = rowvalB[ptrB]
            nzval[ptrS] = -nzvalB[ptrB]
            ptrS += 1
            ptrB += 1
        end

        colptr[col+1] = ptrS
    end
    deleteat!(rowval, ptrS:length(rowval))
    deleteat!(nzval, ptrS:length(nzval))
    return SparseMatrixCSC{Tv,Ti}(m, n-1, colptr, rowval, nzval)
end

diff(a::SparseMatrixCSC, dim::Integer)= dim==1 ? sparse_diff1(a) : sparse_diff2(a)

## norm and rank
vecnorm(A::SparseMatrixCSC, p::Real=2) = vecnorm(A.nzval, p)

function norm(A::SparseMatrixCSC,p::Real=2)
    m, n = size(A)
    if m == 0 || n == 0 || isempty(A)
        return float(real(zero(eltype(A))))
    elseif m == 1 || n == 1
        # TODO: compute more efficiently using A.nzval directly
        return norm(full(A), p)
    else
        Tnorm = typeof(float(real(zero(eltype(A)))))
        Tsum = promote_type(Float64,Tnorm)
        if p==1
            nA::Tsum = 0
            for j=1:n
                colSum::Tsum = 0
                for i = A.colptr[j]:A.colptr[j+1]-1
                    colSum += abs(A.nzval[i])
                end
                nA = max(nA, colSum)
            end
            return convert(Tnorm, nA)
        elseif p==2
            throw(ArgumentError("2-norm not yet implemented for sparse matrices. Try norm(full(A)) or norm(A, p) where p=1 or Inf."))
        elseif p==Inf
            rowSum = zeros(Tsum,m)
            for i=1:length(A.nzval)
                rowSum[A.rowval[i]] += abs(A.nzval[i])
            end
            return convert(Tnorm, maximum(rowSum))
        end
    end
    throw(ArgumentError("invalid p-norm p=$p. Valid: 1, Inf"))
end

# TODO

# kron

function kron{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti}, b::SparseMatrixCSC{Tv,Ti})
    numnzA = nnz(a)
    numnzB = nnz(b)

    numnz = numnzA * numnzB

    mA,nA = size(a)
    mB,nB = size(b)

    m,n = mA*mB, nA*nB

    colptr = Array(Ti, n+1)
    rowval = Array(Ti, numnz)
    nzval = Array(Tv, numnz)

    colptr[1] = 1

    colptrA = a.colptr
    colptrB = b.colptr
    rowvalA = a.rowval
    rowvalB = b.rowval
    nzvalA = a.nzval
    nzvalB = b.nzval

    col = 1

    @inbounds for j = 1:nA
        startA = colptrA[j]
        stopA = colptrA[j+1]-1
        lA = stopA - startA + 1

        for i = 1:nB
            startB = colptrB[i]
            stopB = colptrB[i+1]-1
            lB = stopB - startB + 1

            ptr_range = (1:lB) + (colptr[col]-1)

            colptr[col+1] = colptr[col] + lA * lB
            col += 1

            for ptrA = startA : stopA
                ptrB = startB
                for ptr = ptr_range
                    rowval[ptr] = (rowvalA[ptrA]-1)*mB + rowvalB[ptrB]
                    nzval[ptr] = nzvalA[ptrA] * nzvalB[ptrB]
                    ptrB += 1
                end
                ptr_range += lB
            end
        end
    end
    SparseMatrixCSC(m, n, colptr, rowval, nzval)
end

function kron{Tv1,Ti1,Tv2,Ti2}(A::SparseMatrixCSC{Tv1,Ti1}, B::SparseMatrixCSC{Tv2,Ti2})
    Tv_res = promote_type(Tv1, Tv2)
    Ti_res = promote_type(Ti1, Ti2)
    A = convert(SparseMatrixCSC{Tv_res,Ti_res}, A)
    B = convert(SparseMatrixCSC{Tv_res,Ti_res}, B)
    return kron(A,B)
end

kron(A::SparseMatrixCSC, B::VecOrMat) = kron(A, sparse(B))
kron(A::VecOrMat, B::SparseMatrixCSC) = kron(sparse(A), B)

## det, inv, cond

inv(A::SparseMatrixCSC) = error("The inverse of a sparse matrix can often be dense and can cause the computer to run out of memory. If you are sure you have enough memory, please convert your matrix to a dense matrix.")

# TODO

## scale methods

# multiply by diagonal matrix as vector
function scale!{Tv,Ti}(C::SparseMatrixCSC{Tv,Ti}, A::SparseMatrixCSC, b::Vector)
    m, n = size(A)
    (n==length(b) && size(A)==size(C)) || throw(DimensionMismatch(""))
    numnz = nnz(A)
    C.colptr = convert(Array{Ti}, A.colptr)
    C.rowval = convert(Array{Ti}, A.rowval)
    C.nzval = Array(Tv, numnz)
    for col = 1:n, p = A.colptr[col]:(A.colptr[col+1]-1)
        C.nzval[p] = A.nzval[p] * b[col]
    end
    C
end

function scale!{Tv,Ti}(C::SparseMatrixCSC{Tv,Ti}, b::Vector, A::SparseMatrixCSC)
    m, n = size(A)
    (n==length(b) && size(A)==size(C)) || throw(DimensionMismatch(""))
    numnz = nnz(A)
    C.colptr = convert(Array{Ti}, A.colptr)
    C.rowval = convert(Array{Ti}, A.rowval)
    C.nzval = Array(Tv, numnz)
    for col = 1:n, p = A.colptr[col]:(A.colptr[col+1]-1)
        C.nzval[p] = A.nzval[p] * b[A.rowval[p]]
    end
    C
end

scale{Tv,Ti,T}(A::SparseMatrixCSC{Tv,Ti}, b::Vector{T}) =
    scale!(SparseMatrixCSC(size(A,1),size(A,2),Ti[],Ti[],promote_type(Tv,T)[]), A, b)

scale{T,Tv,Ti}(b::Vector{T}, A::SparseMatrixCSC{Tv,Ti}) =
    scale!(SparseMatrixCSC(size(A,1),size(A,2),Ti[],Ti[],promote_type(Tv,T)[]), b, A)

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# Autogenerated wrapper script for LicenseSpring_jll for x86_64-apple-darwin
export libLicenseSpring

using LibCURL_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("LicenseSpring")
JLLWrappers.@declare_library_product(libLicenseSpring, "@rpath/libLicenseSpring.dylib")
function __init__()
    JLLWrappers.@generate_init_header(LibCURL_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libLicenseSpring,
        "lib/libLicenseSpring.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using StructuralCausalModels, Test

d_str = "dag{k1 -> k2;k2 -> y;v -> x2;w -> k1;x1 -> v;x1 -> w;x2 -> y;x3 -> w;x3 -> y}"

dag = DAG("test_desc", d_str)

to_dagitty(dag) |> display

fname = joinpath(ProjDir, "test_descendents_03.dot")
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

display(dag)
display(dag.a)
println(dag.a[:x1, :w])
println(dag.e[:x1, :w])
println()

bs = basis_set(dag)
bs |> display

@show adjustment_sets(dag, :x2, :y)