text:
module SpatQSims

using FisherySim
using Distributions
using NeutralLandscapes
using Plots
using Random
using HDF5
using CSV
using StructArrays
using Arrow

import FisherySim: simulate, domain
import Plots: plot

include("realizations.jl")
export
    generate_realizations,
    save_realizations,
    load_realization,
    load_movement,
    # save_realizations
    load_spatial_eq,
    load_popstate,
    # save_realizations,
    load_cathcability

include("prep_sims.jl")
export
    prep_sims

include("scale_devs.jl")
export
    scale_devs,
    transform_log_qdevs

include("simprep.jl")
export
    get_realization,
    SpatQSimPrep,
    movement,
    init_pop,
    log_qdevs,
    realization,
    prep_file

include("simspec.jl")
export
    SpatQSimSpec,
    # realization,
    # prep_file,
    domain,
    pop_dynamics,
    survey_targeting,
    comm_targeting,
    base_catchability,
    survey_catchability,
    comm_catchability,
    tweedie_shape,
    tweedie_dispersion,
    n_years,
    survey_vessel,
    comm_vessel,
    fleet

include("qdevscaling.jl")
export
    QDevScalingSpec,
    sim_value,
    comm_qdev_scale

include("sharedq.jl")
export
    SharedQSpec,
    sim_value,
    survey_catchability

include("prefintensity.jl")
export
    PrefIntensitySpec,
    sim_value,
    comm_targeting

include("densdepq.jl")
export
    DensityDependentQSpec,
    sim_value,
    comm_catchability

include("habq.jl")
export
    HabQSpec,
    HabQPrep,
    edge_hab,
    general_hab_distr,
    rocky_hab_distr,
    generate_habitat,
    general_hab_pref,
    rocky_hab_pref_gen,
    generate_habitat,
    move_rate,
    make_moveop,
    sim_value,
    survey_catchability,
    comm_catchability,
    prepare,
    setup,
    habitat,
    domain,
    movement,
    init_pop,
    prepare_habitat,
    get_habitat

include("sim_pars.jl")
export
    sim_values,
    sim_value_idx

include("files.jl")
export
    file_paths,
    simstudy_dir,
    simstudy_prefix,
    prep_path

include("simsetup.jl")
export
    SpatQSimSetup,
    simulate

include("run_sims.jl")
export
    run_sims

include("simresult.jl")
export
    SpatQSimResult,
    # realization,
    file_paths,
    save,
    save_pop_hdf5,
    save_pop_csv,
    save_catch_csv

include("plots.jl")
export
    catch_by_year,
    plot_gif,
    plot

# function run_simulation(scenario::Symbol,
#                         realization::Integer,
#                         prep_fn = "prep.h5")
#     # Don't need habitat for simulation
#     _, movement, init_pop, catchability_devs =
#         get_realization(realization, prep_fn)

#     # 100×100 gridded domain for all simulations
#     Ω = GriddedFisheryDomain()

#     #-Population dynamics-------------------------------------------------------
#     r = 0.06
#     K = 100.0

#     # K defined above to get spatial equilibrium
#     schaefer = Schaefer(r, K)

#     #-Targeting-----------------------------------------------------------------
#     # Commercial targeting depends on population, so isn't constructed prior to
#     # simulation.
#     survey_stations = vec(LinearIndices(Ω.n)[2:4:98, 5:10:95])
#     survey_targeting = FixedTargeting(survey_stations)
#     # Commercial preference is based on spatial distribution of biomass
#     comm_targeting = DynamicPreferentialTargeting(init_pop.P, identity)

#     #-Catchability--------------------------------------------------------------
#     survey_q_base = 0.2
#     comm_q_base = 0.2
#     if scenario == :pref
#         survey_catchability = Catchability(survey_q_base)
#         comm_catchability = Catchability(comm_q_base)
#     elseif scenario == :spat
#         survey_catchability = Catchability(survey_q_base)
#         comm_catchability = Catchability(comm_q_base .* catchability_devs)
#         comm_targeting = RandomTargeting()
#     elseif scenario == :combo
#         survey_catchability = Catchability(survey_q_base)
#         comm_catchability = Catchability(comm_q_base .* catchability_devs)
#     else
#         @error "scenario must be one of :pref, :spat, or :combo."
#     end

#     #-Vessels-------------------------------------------------------------------
#     # Tweedie power parameter
#     ξ = 1.84
#     # Tweedie variance parameter
#     ϕ = 1.2
#     survey_vessel = Vessel(survey_targeting,
#                            survey_catchability,
#                            ξ, ϕ)
#     comm_vessel = Vessel(comm_targeting,
#                          comm_catchability,
#                          ξ, ϕ)

#     #-Fleet---------------------------------------------------------------------
#     fleet = Fleet([survey_vessel, comm_vessel],
#                   [length(survey_vessel.target.locations),
#                    2_500])

#     #-Simulate------------------------------------------------------------------
#     simulate(init_pop, fleet, movement, schaefer, Ω, 25)
# end

# function run_scenarios(repl::Integer, flnm = "prep.h5")
#     Pdict = Dict{Symbol, Vector{PopState}}()
#     Cdict = Dict{Symbol, Vector{Catch}}()
#     for sc in [:pref, :spat, :combo]
#         P, C = run_simulation(sc, repl, flnm)
#         push!(Pdict, sc => P)
#         push!(Cdict, sc => C)
#     end
#     Pdict, Cdict
# end

# function save_scenarios(repl::Integer, Pdict, Cdict)
#     mkpath("repl_" * string(repl, pad = 2))
#     cd("repl_" * string(repl, pad = 2)) do
#         # Save population states to HDF5 file
#         flnm = "pop_" * string(repl, pad = 2) * ".h5"
#         dom_size = size(Pdict[:pref][1])
#         time_size = length(Pdict[:pref])
#         h5open(flnm, "w") do fid
#             for sc in keys(Pdict)
#                 sc_grp = g_create(fid, string(sc))
#                 sc_pop = create_dataset(sc_grp, "popstate", datatype(Float64),
#                                   dataspace(dom_size..., time_size))
#                 for yr in eachindex(Pdict[sc])
#                     sc_pop[:, :, yr] = Pdict[sc][yr].P
#                 end
#             end
#         end

#         for sc in keys(Pdict)
#             popstate = (pop = sum.(Pdict[sc]),)
#             CSV.write("popstate_" * string(repl, pad = 2) *
#                       "_" * string(sc) * ".csv",
#                       StructArray(popstate))
#         end

#         # Save catch records to CSV file
#         for sc in keys(Cdict)
#             flnm = "catch_" * string(repl, pad = 2) *
#                   "_" * string(sc) * ".csv"
#             CSV.write(flnm, StructArray(Cdict[sc]))
#         end
#     end
#     nothing
# end

# function get_coordref(Ω)
#     [(loc_idx = i, s1 = s1, s2 = s2)
#         for (i, (s1, s2)) in enumerate(vec(Ω.locs))]
# end

# function write_coordref(Ω)
#     coordref = get_coordref(Ω)
#     CSV.write("coordref.csv", coordref)
# end

end # module


text:
# This is a test file to make codecov happy, technically all of the
# tests here are already done / integrated within other tests.

@testset "strings" begin
    st = "blah"

    @test J.str(st) == "blah"

    sst = SubString("blahblah", 1:4)
    @test sst == "blah"
    @test J.str(sst) == "blahblah"

    sst = SubString("blah✅💕and etcσ⭒ but ∃⫙∀ done", 1:27)
    @test J.to(sst) == 27

    s = "aabccabcdefabcg"
    for m ∈ eachmatch(r"abc", s)
        @test s[J.matchrange(m)] == "abc"
    end
end

@testset "ocblock" begin
    st = "This is a block  and done"
    τ = J.find_tokens(st, J.MD_TOKENS, J.MD_1C_TOKENS)
    ocb = J.OCBlock(:COMMENT, (τ[1]=>τ[2]))
    @test J.otok(ocb) == τ[1]
    @test J.ctok(ocb) == τ[2]
end

@testset "isexactly" begin
    steps, b, λ = J.isexactly("

ext: jl


text:
# Boiler Plate "using" in seperate file
include(pwd() * "/scripts/dependencies.jl")
include(pwd() * "/scripts/loop_names.jl")

Qlist = [10, 1, 5]
τlist = [0.2, 0.01, 0.05, 0.1]

## Test first
Qlist = [10]
τlist = [0.2]

for Q in Qlist, τ in τlist
    simulationdescription = generate_descriptor(Q, τ)
    println("Currentely doing ", simulationdescription)
end

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
##### Mass tendencies

#####
##### First order fluxes
#####

function flux(::Mass, ::Advect, atmos, args)
    return args.state.ρu
end

#####
##### Second order fluxes
#####

function flux(::Mass, ::MoistureDiffusion, atmos, args)
    @unpack state, diffusive = args
    @unpack D_t = args.precomputed.turbulence
    d_q_tot = (-D_t) .* diffusive.moisture.∇q_tot
    return d_q_tot * state.ρ
end

#####
##### Sources
#####

function source(::Mass, s::Subsidence, m, args)
    @unpack state, aux, diffusive = args
    z = altitude(m, aux)
    w_sub = subsidence_velocity(s, z)
    k̂ = vertical_unit_vector(m, aux)
    return -state.ρ * w_sub * dot(k̂, diffusive.moisture.∇q_tot)
end

function source(::Mass, s::RemovePrecipitation, m, args)
    @unpack state = args
    @unpack ts = args.precomputed
    if has_condensate(ts)
        nt = remove_precipitation_sources(s, m, args)
        return nt.S_ρ_qt
    else
        FT = eltype(state)
        return FT(0)
    end
end

function source(::Mass, s::WarmRain_1M, m, args)
    @unpack cache = args.precomputed.precipitation
    return cache.S_ρ_qt
end

function source(::Mass, s::RainSnow_1M, m, args)
    @unpack cache = args.precomputed.precipitation
    return cache.S_ρ_qt
end


text:
module journal

"""
```Accounting jounal
aka "General Ledger"
to keep tabs on day-to-day monetary activities (to value)
(that we have to keep)
ample Transaction:
Date   Dr account amount
                         Cr Account amount

-In each transaction there is:
1.date of an entry (with time if entries are more frequent, i.e. day-trading more frequent than investor moe frequent than a wage-based worker )
2.'Title' of this  journal entry - Account logistics: what happens, who withraws from which account
3.Comment (optional) ,Explination, elaboration on this account

the amount is written twice in the same transaction:
1. Debit side( on the Right) where the debit account name is written [plus the amount (default:\$ or other currency)]
2. credit side (on the Left) where the credit account namee is written [plus the amount (default:\$ or other currency)  ]
3.auto-checking for match -might be needed  [TODO: reactive journaling with gui]
(though it's hardcoded to yeild the same result)

the transaction is also written
2. in dr account's book
3. in cr account's book
not mentioning the T-accounts for either the:
4. dr account
5. cr account

thus, each transaction could literally get translated into (at least) to be
written down in 5 different seperate entities (structures)
```

"""
using Test
#module journal
#get drAccount, crAccount
@test function transaction(dr = 1, cr = 1, drAccount, crAccount, amount)

  return
end
end

#=

function read()

end

function write()

end

function truncate()

end

function delete()

end

end
=#

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
abstract type WSRuleModel <: MixingRule end

struct WSRule{γ} <: WSRuleModel
    components::Array{String,1}
    activity::γ
end

@registermodel WSRule
export WSRule
function WSRule(components::Vector{String}; activity = Wilson, userlocations::Vector{String}=String[],activity_userlocations::Vector{String}=String[], verbose::Bool=false)
    init_activity = activity(components;userlocations = activity_userlocations,verbose)

    model = WSRule(components, init_activity)
    return model
end

function mixing_rule(model::RKModel,V,T,z,mixing_model::WSRuleModel,α,a,b,c)
    n = sum(z)
    invn2 = (one(n)/n)^2
    num = sum(z[i]*(b[i,i]-a[i,i]*α[i]/(R̄*T)) for i ∈ @comps)/n
    den = 1 - (sum(z[i]*a[i,i]*α[i]/b[i,i]/(n*R̄*T) for i ∈ @comps)-excess_gibbs_free_energy(mixing_model.activity,1e5,T,z)/(n*R̄*T)/log(2))
    c̄ = dot(z,c)/n
    b̄  = num/den
    ā  = R̄*T*(b̄-num)
    return ā,b̄,c̄
end

function mixing_rule(model::PRModel,V,T,z,mixing_model::WSRuleModel,α,a,b,c)
    λ = 1/(2*√(2))*log((2+√(2))/(2-√(2)))
    n = sum(z)
    invn2 = (one(n)/n)^2
    num = sum(z[i]*(b[i,i]-a[i,i]*α[i]/(R̄*T)) for i ∈ @comps)/n
    den = 1 - (sum(z[i]*a[i,i]*α[i]/b[i,i]/(n*R̄*T) for i ∈ @comps)-excess_gibbs_free_energy(mixing_model.activity,1e5,T,z)/(n*R̄*T)/λ)
    c̄ = dot(z,c)/n
    b̄  = num/den
    ā  = R̄*T*(b̄-num)
    return ā,b̄,c̄
end


text:
"""
    masked_crc32c(data)
"""
function masked_crc32c(data)
    x = UInt32(CRC32c.crc32c(data))
    return UInt32(((x >> 15) | UInt32(x << 17)) + 0xa282ead8)
end


text:
#!/usr/bin/env julia

using MAT
using AliasedSeabed

# Creates an aliased seabed mask for the given frequency,
# outputing a MAT file with a "mask" variable.
#
# E.g. ./asb.jl 38 myfile.mat
# Output is myfile-asb18.mat

function main(args)
    freq = args[1]

    filenames = args[2:end]

    for filename in filenames
        a = matread(filename)
        m = asbmask(a["Sv$freq"], a["ntheta$freq"], a["nphi$freq"])

        out ="$(splitext(basename(filename))[1])-asb$freq.mat"

        matwrite(out, Dict("mask" =>  convert(Array{Bool}, m)))
    end
end

main(ARGS)


text:
type SDF_MYD_EP <: SDFDynamicSystem
    omg::Float64
    zta::Float64
    kpa::Float64
    phi::Float64
    alpha::Float64
    u0::Float64
    v0::Float64
    hv::AbstractArray{Float64,1}

    SDF_MYD_EP(omg::Float64, zta::Float64, kpa::Float64, phi::Float64, alpha::Float64) = new(omg, zta, kpa, phi, alpha, 0.0, 0.0, zeros(5))
    SDF_MYD_EP(omg::Float64, zta::Float64, kpa::Float64, phi::Float64) = new(omg, zta, kpa, phi, 0.0, 0.0, 0.0, zeros(5))
end

function fd(s::SDF_MYD_EP, u::Float64; update::Bool=false)
    up, dup, fp, kp, status = s.hv[1:5]
    k0 = s.kpa*s.omg*s.omg
    k1 = s.alpha*k0
    fy = s.phi*g
    uy = fy/k0
    du = u - up
    f_try = fp+du*kp
    bup = fy+k1*(u-uy)
    bdn = -fy+k1*(u+uy)

    if du == 0.0
        kc = kp
        f = fp
    elseif status == 0
        if du > 0.0
            if f_try > bup
                f = bup
                kc = k1
                status = 1
            else
                kc = k0
                f = f_try
            end
        else
            if f_try < bdn
                f = bdn
                kc = k1
                status = 1
            else
                kc = k0
                f = f_try
            end
        end
    elseif status == 1
        if dup*du > 0
            f = f_try
            kc = k1
        else
            kc = k0
            f = fp + kc*du
            status = 0
            if f > bup
                f = bup
                kc = k1
                status = 1
            elseif f < bdn
                f = bdn
                kc = k1
                status = 1
            end
        end
    end
    if update
        s.hv[1:5] = [u, du, f, kc, status]
    end
    f
end

function systemEquation(s::SDF_MYD_EP, p::Function)
    function dydt(y::AbstractArray{Float64,1}, t::Float64)
        u = y[1]
        v = y[2]
        omg2 = s.omg*s.omg
        dvdt = p(t)-2.0*s.zta*s.omg*v-fd(s,u)-omg2*u
        [v; dvdt]
    end
    dydt
end

function response(s::SDF_MYD_EP, p::Function, t::AbstractArray{Float64,1})
    y0 = [s.u0; s.v0]

    f = systemEquation(s, p)

    updateHistVar = (y) -> fd(s, y[1]; update=true)
    y = ode4hv(f, y0, t, updateHistVar)

    u = y[:,1]
    v = y[:,2]
    fdz = 2.0*s.zta*s.omg*v
    fdp = y[:,3]
    a = map(p,t)-fdz-fdp-s.omg*s.omg*u

    [u v a fdz fdp]
end

max_stars_repo_name: Panchatantra/DSR

max_issues_repo_name: Panchatantra/DSR

max_forks_repo_name: Panchatantra/DSR


text:
export splitdims

Base.zeros(x::Number) = zero(x)
Base.ones(x::Number) = one(x)

function pad(x::Array, N::Int)
    nd = ndims(x)
    if nd == N
        x
    elseif nd < N
        dims = fill(1,N)
        for i = 1:ndims(x)
            dims[i] = size(x,i)
        end
        reshape(x, dims...)
    else
        # TODO
        throw("Not implemented yet.")
    end
end

function size3d(x::Array, dim::Int)
    dim == 0 && return (1, length(x), 1)
    dim1, dim2, dim3 = 1, size(x,dim), 1
    for i = 1:dim-1
        dim1 *= size(x, i)
    end
    for i = dim+1:ndims(x)
        dim3 *= size(x, i)
    end
    (dim1, dim2, dim3)
end

"""
    redim(x, n, [pad])
"""
function redim{T,N}(x::Array{T,N}, n::Int; pad=0)
    n == N && return x
    dims = ntuple(n) do i
        1 <= i-pad <= N ? size(x,i-pad) : 1
    end
    reshape(x, dims)
end

#=
import Base: normalize, normalize!
function normalize{T,N}(x::Array{T,N})
    z = mapreducedim(v -> v*v, +, x, N-1)
    @inbounds @simd for i = 1:length(z)
        z[i] = 1 / sqrt(z[i]+1e-9)
    end
    x .* z
end
function normalize!{T,N}(x::Array{T,N})
    z = mapreducedim(v -> v*v, +, x, N-1)
    @inbounds @simd for i = 1:length(z)
        z[i] = 1 / sqrt(z[i]+1e-9)
    end
    x .*= z
    x
end
=#

# Workaround a lack of optimization in gcc
#const exp_cst1 = 2139095040.f0
#const exp_cst2 = 0.f0

#=
@inline function exp_approx(val::Float32)
val2 = 12102203.1615614f0 * val + 1065353216.f0
val3 = val2 < exp_cst1 ? val2 : exp_cst1
val4 = val3 > exp_cst2 ? val3 : exp_cst2
val4i = floor(Int32, val4)
xu = val4i & 0x7F800000
xu2 = (val4i & 0x7FFFFF) | 0x3F800000
b = reinterpret(Float32, Int32(xu2))
xuf = reinterpret(Float32, Int32(xu))
xuf * (0.510397365625862338668154f0 + b *
(0.310670891004095530771135f0 + b *
(0.168143436463395944830000f0 + b *
(-2.88093587581985443087955f-3 + b *
1.3671023382430374383648148f-2))))
end
=#

#=
export fastexp!, normalexp!
const FASTEXP_F32 = Libdl.dlsym(library, :fastexp)
const NORMALEXP_F32 = Libdl.dlsym(library, :normalexp)
function fastexp!{T}(x::Vector{T}, y::Vector{T})
ccall(FASTEXP_F32, Void, (Ptr{T}, Ptr{T}, Cint), x, y, length(x))
end
function normalexp!{T}(x::Vector{T}, y::Vector{T})
ccall(NORMALEXP_F32, Void, (Ptr{T}, Ptr{T}, Cint), x, y, length(x))
end
=#
#=
macro fastmap(f, T, src)
quote
local src = $(esc(src))
local f = $(esc(f))
local T = $(esc(T))
dest = Array(T, length(src))
for i = 1:length(src)
dest[i] = f(src[i])
end
dest
end
end
=#


text:
## Comparison of Optim nnls and this package

using NonNegLeastSquares
import Optim

r = Optim.nnls(randn(5,5),randn(5))
b = nonneg_lsq(randn(5,5),randn(5))

A = randn(2000,2000)
b = randn(2000)

@time r = Optim.nnls(A,b)
@time b = nonneg_lsq(A,b) # roughly 2x faster


text:
using ApproxFunSingularities, ApproxFunOrthogonalPolynomials, ApproxFunBase, InfiniteArrays, BlockArrays, BlockBandedMatrices, IntervalSets, SpecialFunctions, LinearAlgebra, Random, Test
import ApproxFunBase: HeavisideSpace, PointSpace, ArraySpace, DiracSpace, PiecewiseSegment, testbandedoperator, UnionDomain, testtransforms, testfunctional, testbandedblockbandedoperator,
                        resizedata!, CachedOperator, RaggedMatrix

@testset "Sqrt" begin
    x=Fun(identity);
    @time @test sqrt(cos(π/2*x))(.1) ≈ sqrt(cos(.1π/2))

    x=Fun(identity,-2..2)
    @time u=sqrt(4-x^2)/(2π)

    @test u(.1) ≈ sqrt(4-0.1^2)/(2π)
    @test sum(u) ≈ 1

    #this call threw an error, which we check
    @test length(values(u)) == 1

    f = Fun(x->x*cot(π*x/2))
    x = Fun(identity)
    u = Fun(JacobiWeight(1.,1.,ChebyshevInterval()), (f/(1-x^2)).coefficients)
    @test 1/(0.1*cot(π*.1/2)) ≈ (1/u)(.1)

    @test (x/u)(.1) ≈ tan(π*.1/2)

    f=Fun(x->exp(-x^2),Line(0.,0.,-.5,-.5),400)
    @time @test sum(f) ≈ sqrt(π)

    f=Fun(x->exp(x)/sqrt(1-x.^2),JacobiWeight(-.5,-.5))
    @test f(.1) ≈ (x->exp(x)/sqrt(1-x.^2))(.1)

    @time @test norm(Fun(exp,Legendre(0..1))+sqrt(Fun(0..1))) ≈ 2.491141949903508

    @testset "sampling Chebyshev" begin
        x=Fun(identity)
        f = exp(x)/sqrt(1-x^2)
        @time g = cumsum(f)
        @test abs(g(-1)) ≤ 1E-15
        @test g'(0.1) ≈ f(0.1)
    end

    @testset "Complex domains sqrt" begin
        a=1+10*im; b=2-6*im
        d = IntervalCurve(Fun(x->1+a*x+b*x^2))

        x=Fun(d)
        w=sqrt(abs(leftendpoint(d)-x))*sqrt(abs(rightendpoint(d)-x))

        @test sum(w/(x-2.))/(2π*im) ≈ (-4.722196879007759+2.347910413861846im)
        @test linesum(w*log(abs(x-2.)))/π ≈ (88.5579588360686)

        a=Arc(0.,1.,0.,π/2)
        ζ=Fun(identity,a)
        f=Fun(exp,a)*sqrt(abs((ζ-1)*(ζ-im)))
    end

    @time @test norm(Fun(exp,Legendre(0..1))+sqrt(Fun(0..1))) ≈ 2.491141949903508
end

@testset "JacobiWeight" begin
    @testset "Sub-operator re-view bug" begin
        D = Derivative(Chebyshev())
        S = view(D[:, 2:end], Block.(3:4), Block.(2:4))
        @test parent(S) == D
        @test parentindices(S) == (3:4,2:4)
        @test bandwidths(S)  == (-2,2)

        DS=JacobiWeight(1,1,Jacobi(1,1))
        D=Derivative(DS)[2:end,:]
        @test domainspace(D) == DS | (1:∞)
        testbandedoperator(D)
    end

    @testset "Multiplication functions" begin
        x = Fun()
        M = Multiplication(x, JacobiWeight(0,0,Chebyshev()))
        @test exp(M).f == Multiplication(exp(x), Chebyshev()).f
    end

    @testset "Derivative" begin
        S = JacobiWeight(-1.,-1.,Chebyshev(0..1))

        # Checks bug in Derivative(S)
        @test typeof(ConstantSpace(0..1)) <: Space{ClosedInterval{Int},Float64}

        D=Derivative(S)
        f=Fun(S,Fun(exp,0..1).coefficients)
        x=0.1
        @test f(x) ≈ exp(x)*x^(-1)*(1-x)^(-1)/4
        @test (D*f)(x) ≈ -exp(x)*(1+(x-3)*x)/(4*(x-1)^2*x^2)

        S=JacobiWeight(-1.,0.,Chebyshev(0..1))
        D=Derivative(S)

        f=Fun(S,Fun(exp,0..1).coefficients)
        x=.1
        @test f(x) ≈ exp(x)*x^(-1)/2
        @test (D*f)(x) ≈ exp(x)*(x-1)/(2x^2)
    end

    @testset "Jacobi singularity" begin
        x = Fun(identity)
        f = exp(x)/(1-x.^2)

        @test f(.1) ≈ exp(.1)/(1-.1^2)
        f = exp(x)/(1-x.^2).^1
        @test f(.1) ≈ exp(.1)/(1-.1^2)
        f = exp(x)/(1-x.^2).^1.0
        @test f(.1) ≈ exp(.1)/(1-.1^2)

        ## 1/f with poles
        x=Fun(identity)
        f=sin(10x)
        g=1/f

        @test g(.123) ≈ csc(10*.123)
    end

    @testset "Jacobi conversions" begin
        S1,S2=JacobiWeight(3.,1.,Jacobi(1.,1.)),JacobiWeight(1.,1.,Jacobi(0.,1.))
        f=Fun(S1,[1,2,3.])
        C=Conversion(S1,S2)
        Cf=C*f
        @test Cf(0.1) ≈ f(0.1)

        S1,S2=JacobiWeight(3.,2.,Jacobi(1.,1.)),JacobiWeight(1.,1.,Jacobi(0.,0.))
        f=Fun(S1,[1,2,3.])
        C=Conversion(S1,S2)
        Cf=C*f
        @test Cf(0.1) ≈ f(0.1)

        C=Conversion(JacobiWeight(2,0,Jacobi(2,-2)),Legendre())
        @test C[1,1]==4/3 && C[1,2]==-2 && C[2,1]==2 && C[2,2]==-2.4
    end

    @testset "Array Conversion" begin
        a = ArraySpace(JacobiWeight(1/2,1/2, Chebyshev()), 2)
        b = ArraySpace(JacobiWeight(1/2,1/2, Ultraspherical(1)), 2)
        C = Conversion(a, b)

        f = Fun(a, rand(10))
        @test f(0.1) ≈ (C*f)(0.1)

        a = ArraySpace(JacobiWeight(1/2,1/2, Chebyshev()), 2,3)
        b = ArraySpace(JacobiWeight(1/2,1/2, Ultraspherical(1)), 2,3)
        C = Conversion(a, b)

        f = Fun(a, rand(10))
        @test f(0.1) ≈ (C*f)(0.1)
    end

    @testset "Equivalent spaces" begin
        @test norm(Fun(cos,Chebyshev)-Fun(cos,Jacobi(-0.5,-0.5)))<100eps()
        @test norm(Fun(cos,Chebyshev)-Fun(cos,JacobiWeight(0,0)))<100eps()
        @test norm(Fun(cos,Jacobi(-0.5,-0.5))-Fun(cos,JacobiWeight(0,0))) < 100eps()
        @test norm(Fun(cos,Chebyshev)-Fun(cos,JacobiWeight(0,0,Jacobi(-0.5,-0.5))))<100eps()
        @test norm(Fun(cos,Jacobi(-0.5,-0.5))-Fun(cos,JacobiWeight(0,0,Jacobi(-0.5,-0.5))))<100eps()
    end
end

@testset "Ray and Line" begin
    @testset "broadcast!" begin
        f = Fun()
        @test_throws ArgumentError (f .= Fun(Line()))
    end
    @testset "Ray" begin
        @test Inf in Ray()   # this was a bug
        @test Space(0..Inf) == Chebyshev(Ray())
        f=Fun(x->exp(-x),0..Inf)
        @test f(0.1) ≈ exp(-0.1)
        @test f'(.1) ≈ -f(.1)

        x=Fun(identity,Ray())
        f=exp(-x)
        u=integrate(f)
        @test (u(1.)-u(0)-1) ≈ -f(1)

        x=Fun(identity,Ray())
        f=x^(-0.123)*exp(-x)
        @test integrate(f)'(1.) ≈ f(1.)

        @test ≈(sum(Fun(sech,0..Inf)),sum(Fun(sech,0..40));atol=1000000eps())
        @test Line() ∪ Ray() == Line()
        @test Line() ∩ Ray() == Ray()

        f=Fun(sech,Line())
        @test Fun(f,Ray())(2.0) ≈ sech(2.0)
        @test Fun(f,Ray(0.,π))(-2.0) ≈ sech(-2.0)
        @test Fun(sech,Ray(0.,π))(-2.0) ≈ sech(-2.0)
    end

    @testset "Ei (Exp Integral)" begin
        y=Fun(Ray())
        q=integrate(exp(-y)/y)
        @test (q-last(q))(2.) ≈ (-0.04890051070806113)

        ## Line
        f=Fun(x->exp(-x^2),Line())

        @test f'(0.1) ≈ -2*0.1exp(-0.1^2)
        @test (Derivative()*f)(0.1) ≈ -2*0.1exp(-0.1^2)
    end
end

@testset "LogWeight" begin
    x=Fun(identity,-1..1)
    f=exp(x+1)-1
    @test log(f)(0.1) ≈ log(f(0.1))

    x=Fun(identity,0..1)
    f=exp(x)-1
    @test log(f)(0.1) ≈ log(f(0.1))

    x=Fun(identity,0..1)
    @test Fun(exp(x)/x-1/x,Chebyshev)(0.1) ≈ (exp(0.1)-1)/0.1

    x=Fun(identity,0..1)
    f=1/x
    p=integrate(f)
    @test (p-p(1.))(0.5) ≈ log(0.5)

    f=1/(1-x)
    p=integrate(f)
    @test (p-p(0.))(0.5) ≈ -log(1-0.5)

    @testset "#393" begin
        x = Fun(0..1)
        @time f = exp(x)*sqrt(x)*log(1-x)
        @test f(0.1) ≈ exp(0.1)*sqrt(0.1)*log(1-0.1)
    end
end

@testset "DiracDelta and PointSpace" begin
    a,b=DiracDelta(0.),DiracDelta(1.)
    f=Fun(exp)
    g=a+0.2b+f
    @test components(g)[2](0.) ≈ 1.
    @test g(.1) ≈ exp(.1)
    @test sum(g) ≈ (sum(f)+1.2)

    #Checks prevoius bug
    δ=DiracDelta()
    x=Fun()
    w=sqrt(1-x^2)
    @test (w+δ)(0.1) ≈ w(0.1)
    @test sum(w+δ) ≈ sum(w)+1

    ## PointSpace
    f=Fun(x->(x-0.1),PointSpace([0,0.1,1]))
    g = f + Fun(2..3)
    @test f(0.0) ≈ g(0.0) ≈ -0.1
    @test f(0.1) ≈ g(0.1) ≈ 0.0
    @test f(1.0) ≈ g(1.0) ≈ 0.9

    @test g(2.3) ≈ 2.3

    h = a + Fun(2..3)

    # for some reason this test is broken only on Travis
    @test_skip g/h ≈ f/a + Fun(1,2..3)
end

@testset "Multiple roots" begin
    x=Fun(identity,-1..1)

    @test (1/x^2)(0.1) ≈ 100.
    @test (1/x^2)(-0.1) ≈ 100.

    fc=x*(1+x)^2
    @time @test (1/fc)(0.1) ≈ 1/fc(0.1)

    fc=x*(1-x)^2
    @test (1/fc)(0.1) ≈ 1/fc(0.1)
end

@testset "special function singularities" begin
    x=Fun(0..1)
    @time @test erf(sqrt(x))(0.1) ≈ erf(sqrt(0.1))
    @time @test erfc(sqrt(x))(0.1) ≈ erfc(sqrt(0.1))

    ## roots of log(abs(x-y))
    x=Fun(-2..(-1))
    @test space(abs(x)) == Chebyshev(-2 .. (-1))

    @test roots(abs(x+1.2)) ≈ [-1.2]

    f = sign(x+1.2)
    @test space(f) isa PiecewiseSpace
    @test f(-1.4) == -1
    @test f(-1.1) == 1

    f=abs(x+1.2)
    @test abs(f)(-1.3) ≈ f(-1.3)
    @test abs(f)(-1.1) ≈ f(-1.1)
    @test norm(abs(f)-f)<10eps()

    @test norm(sign(f)-Fun(1,space(f)))<10eps()

    @test log(f)(-1.3) ≈ log(abs(-1.3+1.2))
    @test log(f)(-1.1) ≈ log(abs(-1.1+1.2))
end

@testset "SumSpace" begin
    @testset "SumSpace bug" begin
        dsp=JacobiWeight(1.,0.,Jacobi(1.,0.,0..1))⊕JacobiWeight(0.5,0.,Jacobi(0.5,-0.5,0..1))
        rsp=Legendre(0..1)⊕JacobiWeight(0.5,0.,Jacobi(0.5,0.5,0..1))

        C=Conversion(dsp,rsp)

        f=Fun(dsp,[1.,2.,3.,4.,5.])
        @test f(0.1) ≈ (C*f)(0.1)
    end

    @testset "Triple SumSpace" begin
        x=Fun()
        w=log(1-x)+sqrt(1-x^2)
        @time f=w+x
        @test f(0.1) ≈ (w(0.1)+0.1)
        @test (w+1)(0.1) ≈ (w(0.1)+1)
        @test (w+x+1)(0.1) ≈ (w(0.1)+1.1)
        @test ((w+x)+1)(0.1) ≈ (w(0.1)+1.1)
    end

    @testset "Previoius segfault" begin
        x=Fun(identity,-1..1)
        @time f=x+sin(2x)*sqrt(1-x^2)
        @test f(0.1) ≈ 0.1+sin(2*0.1)*sqrt(1-0.1^2)
    end

    @testset "Multiple piecewisespace" begin
        x=Fun(identity,-3 .. -2)+Fun(identity,2..3)
        w=sqrt(9-x^2)
        f=w+Fun()
        @test (f+w)(2.5) ≈ 2w(2.5)
        @test (f+w)(.5) ≈ f(.5)
    end

    @testset "cancellation conversion" begin
        x=Fun(0..1)
        f=exp(x)-1
        Fun(f,JacobiWeight(1.,0.,0..1))
    end

    @testset "Union of ChebyshevDirichlet" begin
        dom = UnionDomain(0..1, 2..3)
        @test components(union(JacobiWeight.(-0.5,-0.5,ChebyshevDirichlet{1,1}.(components(dom)))...)) ==
            (JacobiWeight.(-0.5,-0.5,ChebyshevDirichlet{1,1}.(components(dom)))...,)
    end

    @testset "Ultraspherical special functions" begin
        x = Fun(Ultraspherical(2,0..1))
        sqrt(x)(0.1) ≈ sqrt(0.1)

        f = Fun(x->x*exp(x),Ultraspherical(1,0..1))
        sqrt(f(0.1)) ≈ sqrt(f)(0.1)
    end

    @testset "one for SumSpace" begin
        S = Jacobi(0,1) ⊕ JacobiWeight(1/3,0,Jacobi(1/3,2/3)) ⊕ JacobiWeight(2/3,0,Jacobi(2/3,1/3))
        o = ones(S)
        @test o(0.5) ≈ 1
    end

    @testset "Complex piecewise" begin
        x = Fun(identity, Segment(im,0) ∪ Segment(0,1))

        @test sqrt(1-x)(0.2im) ≈ sqrt(1-0.2im)
        @test sqrt(1-x)(0.2) ≈ sqrt(1-0.2)

        w=2/(sqrt(1-x)*sqrt(1+im*x))
        @test w(0.2im) ≈ 2/(sqrt(1-0.2im)*sqrt(1+im*(0.2im)))
        @test w(0.2) ≈ 2/(sqrt(1-0.2)*sqrt(1+im*(0.2)))
    end
end

@testset "WeightedJacobi" begin
    m=20
    @time testtransforms(JacobiWeight(0.,m,Jacobi(0.,2m+1)))
    f=Fun(x->((1-x)/2).^m.*exp(x),JacobiWeight(0.,m,Jacobi(0.,2m+1)))
    @test abs(f(.1)-(x->((1-x)/2).^m.*exp(x))(.1))<10eps()

    m=10
    @time f=Fun(x->besselj(m,m*(1-x)),JacobiWeight(0.,m,Jacobi(0.,2m+1)))
    @test f(0.) ≈ besselj(m,m)

    @testset "Conversion" begin
        testtransforms(Jacobi(-0.5,-0.5))
        @test norm(Fun(Fun(exp),Jacobi(-.5,-.5))-Fun(exp,Jacobi(-.5,-.5))) < 300eps()

        x=Fun(identity)
        ri=0.5/(1-x)
        @test ((1-x)/2 .* Fun(exp,JacobiWeight(0.,0.,Jacobi(0.,1.))))(.1) ≈ (1-.1)./2*exp(.1)

        @test ((1-x)/2 .* Fun(exp,JacobiWeight(0.,0.,Jacobi(0.,1.))))(.1) ≈ (1-.1)./2*exp(.1)

        @test (ri*Fun(exp,JacobiWeight(0.,0.,Jacobi(0.,1.))))(.1) ≈ .5/(1-.1)*exp(.1)
    end

    @testset "Derivative" begin
        D=Derivative(Jacobi(0.,1.,Segment(1.,0.)))
        @time testbandedoperator(D)

        S=JacobiWeight(0.,0.,Jacobi(0.,1.,Segment(1.,0.)))
        D=Derivative(S)
        testbandedoperator(D)

        f=Fun(exp,domainspace(D))
        @test (D*f-f).coefficients|>norm < eps(100000.)
        @test (f'-f).coefficients|>norm < eps(100000.)
        @test (D^2*f-f).coefficients|>norm < eps(100000000.)
        @test (D*(D*f)-f).coefficients|>norm < eps(100000000.)

        S=JacobiWeight(1,1,Ultraspherical(1))

        f=Fun(S,[1.,2.,3.])
        @test (Derivative(S,2)*f)(0.1) ≈ f''(0.1)
    end

    @testset "special derivative" begin
        x=Fun()
        f=exp(x)*sqrt(1-x^2)
        D=Derivative(WeightedJacobi(.5,.5))

        testtransforms(WeightedJacobi(.5,.5))
        testbandedoperator(D)

        @time g=(D*Fun(f,domainspace(D)))
        @test f'(0.1) ≈ g(0.1)
    end

    @testset "==" begin
        @test WeightedJacobi(0.1,0.2) == WeightedJacobi(0.1+eps(),0.2)
    end

    @testset "subspace bug" begin
        f=Fun(WeightedJacobi(0.1,0.2),rand(10))  # convert to Legendre expansion

        g=(f|(2:∞))

        @test ApproxFunBase.coefficients(g.coefficients,space(g),ApproxFunBase.canonicalspace(g))[1] ==0.
        @test norm((Fun(g,space(f))|(2:∞)-g).coefficients) < 10eps()
    end

    @testset "JacobiWeight cumsum bug Issue #557" begin
        x = Fun(0.0..1.0)
        ν = 2
        @time f = x^(ν/2-1) * exp(-x/2) # 0.05s
        @test cumsum(f)' ≈ f
        @test cumsum(f)(1.0) ≈ 0.7869386805747332 # Mathematic

        x = Fun(Ray())
        ν = 2
        @time f = x^(ν/2-1) * exp(-x/2) # 0.05s
        @test cumsum(f)' ≈ f
        @test cumsum(f)(1.0) ≈ 0.7869386805747332
    end

    @testset "Definite integral" begin
        for S in (WeightedJacobi(0,0), JacobiWeight(0,0, Legendre(1.1..2.3)), Legendre())
            B = DefiniteIntegral(S)
            testfunctional(B)
            @test ApproxFunBase.rowstop(B,1) == 1
            B[1] == arclength(domain(S))
            f = Fun(exp, S)
            B*f == sum(Fun(exp,domain(S)))
        end
    end
end

@testset "WeightedLaguerre" begin
    @testset "WeightedLaguerre cumsum" begin
        α = 2.7
        f = Fun(WeightedLaguerre(α), [1.0]);
        f = Fun(f, JacobiWeight(α,0,Chebyshev(Ray())));
        g = integrate(f)
        g(3.0) - cumsum(Fun(x -> f(x), 0..6))(3.0)
    end

    @testset "Log with squareroot singularities" begin
        a = 1.0; b = 2.0+im
        d = Segment(a,b)
        z = Fun(d)

        f = real(exp(z) / (sqrt(z-a)*sqrt(b-z)))
        S=space(f)
        x=4.0+2im;
        @test linesum(f*log(abs(x-z))) ≈ 13.740676344264614
    end
end

@testset "Hermite Integration" begin
    @test_throws ArgumentError integrate(Fun(GaussWeight(Hermite(2),1), [0.0,1.0]))

    w = Fun(GaussWeight(Hermite(2), 0), [1.0,2.0,3.0])
    g = integrate(w)
    g̃ = Fun(Hermite(2), [0.0, 0.5, 0.5, 0.5])
    @test g(0.1) == g̃(0.1)

    w = Fun(GaussWeight(), Float64[])
    g = integrate(w)
    @test g(0.1) == 0.0

    w = Fun(GaussWeight(), [1.0])
    g = integrate(w)
    @test_skip w̃ = Fun(w, -7..7)
    w̃ = Fun( x-> w(x), -7..7)
    g̃ = cumsum(w̃)
    @test g(3) - g(-7) ≈ g̃(3)

    w = Fun(GaussWeight(), Float64[1.0])
    g = integrate(w)
    @test_skip w̃ = Fun(w, -7..7)
    w̃ = Fun(x -> w(x), -7..7)
    g̃ = cumsum(w̃)
    @test g(3) - g(-7) ≈ g̃(3)

    w = Fun(GaussWeight(Hermite(2), 2), Float64[1.0])
    g = integrate(w)
    @test_skip w̃ = Fun(w, -7..7)
    w̃ = Fun(x -> w(x), -7..7)
    g̃ = cumsum(w̃)
    @test g(3) - g(-7) ≈ g̃(3)

    w = Fun(GaussWeight(), Float64[0.0, 1.0])
    g = integrate(w)
    @test_skip w̃ = Fun(w, -7..7)
    w̃ = Fun(x -> w(x), -7..7)
    g̃ = cumsum(w̃)
    @test g(3) - g(-7) ≈ g̃(3)

    w = Fun(GaussWeight(Hermite(2), 2), Float64[0.0, 1.0])
    g = integrate(w)
    @test_skip w̃ = Fun(w, -7..7)
    w̃ = Fun(x -> w(x), -7..7)
    g̃ = cumsum(w̃)
    @test g(3) - g(-7) ≈ g̃(3)
end

@testset "Laguerre" begin
    @testset "Integration" begin
        @test_throws ArgumentError integrate(Fun(LaguerreWeight(1, Laguerre(2)), [1.0, 2.0]))

        g = integrate(Fun(WeightedLaguerre(), []))
        @test g(0.1) == 0.0

        α = 2.8
        f = Fun(WeightedLaguerre(α), [1.0])
        g = integrate(f)
        @test g(3.0) - cumsum(Fun(x -> f(x), 0..6))(3.0) ≈ g(0.0)

        α = 2
        f = Fun(WeightedLaguerre(α), [1.0])
        g = integrate(f)
        @test g(3.0) - cumsum(Fun(x -> f(x), 0..6))(3.0) ≈ g(0.0)

        α = 2.8
        f = Fun(WeightedLaguerre(α), [0.0, 1.0])
        g = integrate(f)
        f̃ = Fun(x -> f(x), 0 .. 100)
        g̃ = integrate(f̃)
        g̃ = g̃ - last(g̃)
        @test g(3.0) ≈ g̃(3.0)

        α = 2
        f = Fun(WeightedLaguerre(α), [1.0, 1.0])
        g = integrate(f)
        @test g(3.0) - cumsum(Fun(x -> f(x), 0..6))(3.0) ≈ g(0.0)

        α = 5
        f = Fun(WeightedLaguerre(α), [1.0])
        g = integrate(f)
        @test g(3.0) - cumsum(Fun(x -> f(x), 0..6))(3.0) ≈ g(0.0)

        α = 5
        f = Fun(WeightedLaguerre(α), [0.0, 1.0])
        g = integrate(f)
        f̃ = Fun(x -> f(x), 0 .. 100)
        g̃ = integrate(f̃)
        g̃ = g̃ - last(g̃)
        @test g(3.0) ≈ g̃(3.0)

        α = 5
        f = Fun(WeightedLaguerre(α), [1.0, 1.0])
        g = integrate(f)
        @test g(3.0) - cumsum(Fun(x -> f(x), 0..6))(3.0) ≈ g(0.0)
    end

    @testset "Correct domain" begin
        w = Fun(WeightedLaguerre(0.5),[1.0])
        h = cumsum(w)
        @test domain(h) == Ray()
    end
end

include("IntegralEquationsTest.jl")

@testset "PDE" begin
    @testset "Laplacian" begin
        S=WeightedJacobi(1,1)^2
        L=Laplacian(S)
        testbandedblockbandedoperator(L)
    end
    @testset "Zero Dirichlet" begin
        S = JacobiWeight(1.,1.,Jacobi(1.,1.))^2
        Δ = Laplacian(S)

        testbandedblockbandedoperator(Δ)

        u = Fun((x,y)->sin(π*x)*sin(π*y),S)
        f = -2π^2*u

        F = qr(Δ)
        ApproxFunBase.resizedata!(F,:,1000)
        @time v=F\f
        @test norm((u-v).coefficients)<100eps()

        F=qr(Δ)
        ApproxFunBase.resizedata!(F.R_cache,:,100)
        ApproxFunBase.resizedata!(F.R_cache,:,1000)
        @time v=F \ f
        @test norm((u-v).coefficients)<100eps()

        F=qr(Δ)
        @time v=F\f
        @test norm((u-v).coefficients)<100eps()
    end
    @testset "Gaussian zero Dirichlet Poisson" begin
        S=JacobiWeight(1.,1.,Jacobi(1.,1.))^2
        Δ=Laplacian(S)

        f=Fun((x,y)->exp(-10(x+.2)^2-20(y-.1)^2),rangespace(Δ))  #default is [-1,1]^2
        @time v = \(Δ,f;tolerance=1E-14)
        @test norm((Δ*v-f).coefficients) < 1E-14

    end
    @testset "check we dispatch correctly to get fast build" begin
        S = JacobiWeight(1.,1.,Jacobi(1.,1.))^2
        Δ = Laplacian(S)
        @time S = view(Δ.op.ops[1].ops[1].op,Block.(1:40), Block.(1:40))
        @test typeof(S.parent.domaintensorizer) == ApproxFunBase.Trivial2DTensorizer
    end

    @testset "Operator resize" begin
        S=ChebyshevDirichlet()^2
        B=Dirichlet(S)
        f = Fun((x,y)->exp(x)*sin(y),S)
        @test norm((Fun((x,y)->exp(x)*sin(y),∂(domain(S))) - B*f).coefficients) < 100eps()

        S=JacobiWeight(1.,1.,Jacobi(1.,1.))^2
        Δ=Laplacian(S)

        @test cache(Δ)[1:100,1:100]  ≈ Δ[1:100,1:100]
        @test cache(Δ;padding=true)[1:100,1:100]  ≈ Δ[1:100,1:100]

        @test cache(Δ)[5:100,7:100]  ≈ Δ[5:100,7:100]
        @test cache(Δ;padding=true)[5:100,7:100]  ≈ Δ[5:100,7:100]

        # Check that QR is growing correctly
        for col in (1,2,3,10,11,40)
            QR=qr(Δ)
            resizedata!(QR.R_cache,:,col+100)
            resizedata!(QR,:,col)
            QR2=qr!(CachedOperator(RaggedMatrix,Δ;padding=true))
            resizedata!(QR2.R_cache,:,QR.ncols+100)
            resizedata!(QR2,:,QR.ncols)
            n=min(size(QR.H,1),size(QR2.H,1))
            @test QR.H[1:n,1:col] ≈ QR2.H[1:n,1:col]
            @test QR.R_cache[1:col,1:col] ≈ QR2.R_cache[1:col,1:col]
            @test QR.R_cache[1:col+10,1:col+10] ≈ QR2.R_cache[1:col+10,1:col+10]
        end

        QRR=qr(Δ)
        QR2=qr!(CachedOperator(RaggedMatrix,Δ;padding=true))
        for col in (80,200)
            resizedata!(QRR,:,col)
            resizedata!(QR2,:,QRR.ncols)
            n=min(size(QRR.H,1),size(QR2.H,1))
            @test QRR.H[1:n,1:col] ≈ QR2.H[1:n,1:col]
            @test QRR.R_cache[1:col,1:col] ≈ QR2.R_cache[1:col,1:col]
            @test QRR.R_cache[1:col+10,1:col+10] ≈ QR2.R_cache[1:col+10,1:col+10]
        end

        # this checks a bug
        QRR=qr(Δ)
        resizedata!(QRR,:,548)
        resizedata!(QRR,:,430)

        u=Fun((x,y)->sin(π*x)*sin(π*y),S)
        f=-2π^2*u

        QRR=qr(Δ)
        v=QRR\f
        @test norm((u-v).coefficients)<100eps()

        v=Δ\f
        @test norm((u-v).coefficients)<100eps()

        f=Fun((x,y)->exp(-10(x+.2)^2-20(y-.1)^2),rangespace(Δ))  #default is [-1,1]^2
        @time v=\(Δ,f;tolerance=1E-14)
        @test norm((Δ*v-f).coefficients)<1E-14

        KO=Δ.op.ops[1].ops[1].op

        M=BandedBlockBandedMatrix(view(KO,1:4,1:4))
        @test norm(BandedBlockBandedMatrix(view(KO,1:4,2:4))-M[:,2:4]) < 10eps()
        @test norm(BandedBlockBandedMatrix(view(KO,1:4,3:4))-M[:,3:4]) < 10eps()
    end
end

@testset "Cauchy fun" begin
    f = Fun((x,y)->1/(2π*(x^2+y^2+1)^(3/2)),Line()^2)
    @test f(0.1,0.2) ≈ 1/(2π*(0.1^2+0.2^2+1)^(3/2))

    #TODO: improve tolerance
    f = LowRankFun((x,y)->1/(2π*(x^2+y^2+1)^(3/2)),JacobiWeight(2.,2.,Line())^2)
    @test ≈(f(0.1,0.2),1/(2π*(0.1^2+0.2^2+1)^(3/2));atol=1E-4)
end

@testset "Bessel" begin
    @time for ν in (1.,0.5,2.,3.5)
        println("        ν = $ν")
        S=JacobiWeight(-ν,0.,Chebyshev(0..1))
        D=Derivative(S)
        x=Fun(identity,domain(S))
        L=(x^2)*D^2+x*D+(x^2-ν^2);
        u=\([rdirichlet(S);rneumann(S);L],[bessely(ν,1.),.5*(bessely(ν-1.,1.)-bessely(ν+1.,1.)),0];
                    tolerance=1E-10)
        @test ≈(u(.1),bessely(ν,.1);atol=eps(1000000.)*max(abs(u(.1)),1))
        u=Fun(x->bessely(ν,x),S)
        @test ≈(u(.1),bessely(ν,.1);atol=eps(10000.)*max(abs(u(.1)),1))
        u=Fun(x->besselj(ν,x),S)
        @test ≈(u(.1),besselj(ν,.1);atol=eps(10000.)*max(abs(u(.1)),1))
    end

    @time for ν in (1.,0.5,0.123,3.5)
        println("        ν = $ν")
        S=JacobiWeight(ν,0.,Chebyshev(0..1))
        D=Derivative(S)
        x=Fun(identity,domain(S))
        L=(x^2)*D^2+x*D+(x^2-ν^2);

        u=\([rdirichlet(S);rneumann(S);L],[besselj(ν,1.),.5*(besselj(ν-1.,1.)-besselj(ν+1.,1.)),0];
                    tolerance=1E-10)
        @test ≈(u(.1),besselj(ν,.1);atol=eps(1000000.)*max(abs(u(.1)),1))
        u=Fun(x->besselj(ν,x),S)
        @test ≈(u(.1),besselj(ν,.1);atol=eps(10000.)*max(abs(u(.1)),1))
    end
end

@testset "Speed test" begin
    S = JacobiWeight(1.,1.,Jacobi(1.,1.))^2
    Δ = Laplacian(S)

    f = Fun((x,y)->sin(π*x)*sin(π*y),S)

    QR1=qr(Δ)
    ApproxFunBase.resizedata!(QR1,:,400)
        \(QR1,f; tolerance=1E-10)
    QR1=qr(Δ)
        @time Δ[Block.(1:40), Block.(1:40)]
        @time ApproxFunBase.resizedata!(QR1,:,400)
        @time \(QR1,f; tolerance=1E-10)
    println("Laplace Dirichlet: should be ~0.015, 0.015, 0.001")
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    Dates

The `Dates` module provides `Date`, `DateTime`, `Time` types, and related functions.

The types are not aware of time zones, based on UT seconds
(86400 seconds a day, avoiding leap seconds), and
use the proleptic Gregorian calendar, as specified in [ISO 8601](https://en.wikipedia.org/wiki/ISO_8601).
For time zone functionality, see the TimeZones.jl package.

```jldoctest
julia> dt = DateTime(2017,12,31,23,59,59,999)
2017-12-31T23:59:59.999

julia> d1 = Date(Dates.Month(12), Dates.Year(2017))
2017-12-01

julia> d2 = Date("2017-12-31", Dates.DateFormat("y-m-d"))
2017-12-31

julia> Dates.yearmonthday(d2)
(2017, 12, 31)

julia> d2-d1
30 days
```

Please see the manual section on [`Date`](@ref) and [`DateTime`](@ref)
for more information.
"""
module Dates

import ..Base: ==, div, fld, mod, rem, gcd, lcm, +, -, *, /, %
import ..Base.broadcast

using Base.Iterators

include("types.jl")
include("periods.jl")
include("accessors.jl")
include("query.jl")
include("arithmetic.jl")
include("conversions.jl")
include("ranges.jl")
include("adjusters.jl")
include("rounding.jl")
include("io.jl")
include("parse.jl")

export Period, DatePeriod, TimePeriod,
       Year, Month, Week, Day, Hour, Minute, Second, Millisecond,
       Microsecond, Nanosecond,
       TimeZone, UTC, TimeType, DateTime, Date, Time,
       # periods.jl
       canonicalize,
       # accessors.jl
       yearmonthday, yearmonth, monthday, year, month, week, day,
       hour, minute, second, millisecond, dayofmonth,
       microsecond, nanosecond,
       # query.jl
       dayofweek, isleapyear, daysinmonth, daysinyear, dayofyear, dayname, dayabbr,
       dayofweekofmonth, daysofweekinmonth, monthname, monthabbr,
       quarterofyear, dayofquarter,
       Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, Sunday,
       Mon, Tue, Wed, Thu, Fri, Sat, Sun,
       January, February, March, April, May, June,
       July, August, September, October, November, December,
       Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec,
       # conversions.jl
       unix2datetime, datetime2unix, now, today,
       rata2datetime, datetime2rata, julian2datetime, datetime2julian,
       # adjusters.jl
       firstdayofweek, lastdayofweek,
       firstdayofmonth, lastdayofmonth,
       firstdayofyear, lastdayofyear,
       firstdayofquarter, lastdayofquarter,
       adjust, tonext, toprev, tofirst, tolast,
       # io.jl
       ISODateTimeFormat, ISODateFormat, ISOTimeFormat, DateFormat, RFC1123Format, @dateformat_str

end # module

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
using RandomNumbers
using DataFrames
using Plots, StatsPlots

@inline @fastmath function randbn(n,p,rng)
    q = 1.0 - p
    s = p/q
    a = (n+1)*s
    r = exp(n*log(q))
    x = 0
    u = rand(rng)
    while true
        if (u < r)
            return x
        end
        u -= r
        x += 1
        r *= (a/x)-s
    end
end

@inline @fastmath function sir(u, parms, rng)
    (S, I, R, Y) = u
    (β, γ, ι, N, δt) = parms
    λ = β * (I + ι) / N
    ifrac = 1.0 - exp(-λ * δt)
    rfrac = 1.0 - exp(-γ * δt)
    infection = randbn(S, ifrac, rng)
    recovery = randbn(I, rfrac, rng)
    return (S - infection, I + infection - recovery, R + recovery, Y + infection)
end

function simulate(r)
    parms = (0.1, 0.05, 0.01, 1000.0, 0.1)
    tf = 200
    t = 0:0.1:tf
    tl = length(t)
    S = zeros(tl)
    I = zeros(tl)
    R = zeros(tl)
    Y = zeros(tl)
    u0 = (999, 1, 0, 0)
    (S[1],I[1],R[1],Y[1]) = u0
    u = u0
    for i in 2:tl
        u = sir(u, parms, r)
        (S[i],I[i],R[i],Y[i]) = u
    end
    return DataFrame(Time=t,S=S,I=I,R=R,Y=Y)
end

seed = 42
r = Xorshifts.Xorshift128Plus(seed);

sir_out = simulate(r);

head(sir_out)

# Plot
@df sir_out plot(:Time, [:S :I :R], colour = [:red :green :blue], xlabel="Time",ylabel="Number")

max_stars_repo_name: epimodels/epicookbook

max_issues_repo_name: epimodels/epicookbook

max_forks_repo_name: epimodels/epicookbook


text:
"""
    SOCtoRSOCBridge{T} <: Bridges.Variable.AbstractBridge

Same transformation as [`MOI.Bridges.Constraint.SOCRBridge`](@ref).
"""
struct SOCtoRSOCBridge{T} <: AbstractBridge
    variables::Vector{MOI.VariableIndex}
    constraint::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.RotatedSecondOrderCone}
end
function bridge_constrained_variable(
    ::Type{SOCtoRSOCBridge{T}}, model::MOI.ModelLike,
    set::MOI.SecondOrderCone) where T
    variables, constraint = MOI.add_constrained_variables(
        model, MOI.RotatedSecondOrderCone(MOI.dimension(set)))
    return SOCtoRSOCBridge{T}(variables, constraint)
end

function supports_constrained_variable(
    ::Type{<:SOCtoRSOCBridge}, ::Type{MOI.SecondOrderCone})
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:SOCtoRSOCBridge})
    return [(MOI.RotatedSecondOrderCone,)]
end
function MOIB.added_constraint_types(::Type{<:SOCtoRSOCBridge})
    return Tuple{DataType, DataType}[]
end

# Attributes, Bridge acting as a model
function MOI.get(bridge::SOCtoRSOCBridge, ::MOI.NumberOfVariables)
    return length(bridge.variables)
end
function MOI.get(bridge::SOCtoRSOCBridge, ::MOI.ListOfVariableIndices)
    return bridge.variables
end
function MOI.get(bridge::SOCtoRSOCBridge,
                 ::MOI.NumberOfConstraints{MOI.VectorOfVariables,
                                           MOI.RotatedSecondOrderCone})
    return 1
end
function MOI.get(bridge::SOCtoRSOCBridge,
                 ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,
                                               MOI.RotatedSecondOrderCone})
    return [bridge.constraint]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::SOCtoRSOCBridge)
    MOI.delete(model, bridge.variables)
end

# Attributes, Bridge acting as a constraint

function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::SOCtoRSOCBridge{T}) where T
    return MOI.SecondOrderCone(length(bridge.variables))
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal, MOI.ConstraintDual},
    bridge::SOCtoRSOCBridge)
    return rotate_result(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.VariablePrimal,
                 bridge::SOCtoRSOCBridge, i::IndexInVector)
    return rotate_result(model, attr, bridge.variables, i)
end

function MOIB.bridged_function(bridge::SOCtoRSOCBridge{T}, i::IndexInVector) where T
    return rotate_bridged_function(T, bridge.variables, i)
end
function unbridged_map(bridge::SOCtoRSOCBridge{T}, vis::Vector{MOI.VariableIndex}) where T
    return rotate_unbridged_map(T, bridge.variables, vis)
end


text:
const ALL_DIFFERENT = Dict(
    :concept => all_unique,
)


text:
using ProgressLogging
import Base: size, getindex, setindex!, similar, eltype, deleteat!, filter, union, intersect, convert, promote_rule, +, \

abstract type AbstractFeatureSet <: AbstractVector{Function} end
export AbstractFeatureSet

"""
    FeatureSet(methods, [names, keywords, descriptions])
    FeatureSet(features::Vector{T}) where {T <: AbstractFeature}

Construct a `FeatureSet` from `methods` (a vector of functions) and optionally provide `names` as a vector of symbols, `keywords` as a vector of vectors of strings and `descriptions` as a vector of strings.
A `FeatureSet` can be called on a time series vector or matrix `X` (with time series occupying columns) to return a `FeatureArray` of feature values.
Subsets of a `FeatureSet` `𝒇` can be obtained by indexing with feature names (as symbols) or the regular linear and logical indices.
`FeatureSet`s also support simple set operations defined for arrays, such as unions and intersections, as well as convenient syntax for concatenation (`+`) and set differencing (`\\`).
Note that two features are considered the same if and only if their names are equal.

# Examples
```julia
𝒇 = FeatureSet([sum, length], [:sum, :length], [["distribution"], ["sampling"]], ["∑x¹", "∑x⁰"])
X = randn(100, 2) # 2 time series, 100 samples long
F = 𝒇(X)

# Joining feature sets
𝒇₁ = FeatureSet([x->min(x...), x->max(x...)], [:min, :max], [["distribution"], ["distribution"]], ["minimum", "maximum"])
𝒈₁ = 𝒇 + 𝒇₁
G = 𝒈₁(X)

# Intersecting feature sets, where features are identified exclusively by their names
𝒇₂ = FeatureSet(x->prod, :sum, ["distributions"], "∏x")
𝒈₂ = 𝒇 ∩ 𝒇₂ # The intersection of two feature sets, both with their own :sum
G = 𝒈₂(X) # The intersection contains the :sum of the first argument to ∩; 𝒇
```
"""
struct FeatureSet <: AbstractFeatureSet
    features::Vector{AbstractFeature}
    FeatureSet(features::Vector{T}) where {T <: AbstractFeature} = new(features)
end

FeatureSet( methods::AbstractVector,
            names=Symbol.(methods),
            keywords=fill([], length(methods)),
            descriptions=fill("", length(methods))) =
            FeatureSet(Feature.(methods, names, keywords, descriptions))

FeatureSet( methods::Function,
            names=Symbol(methods),
            keywords=[],
            descriptions="") =
            FeatureSet([Feature(methods, names, keywords, descriptions)])

FeatureSet(f::AbstractFeature) = FeatureSet([f])

export FeatureSet

getfeatures(𝒇::AbstractFeatureSet) = 𝒇.features
getmethods(𝒇::AbstractFeatureSet)  = getmethod.(𝒇)
getnames(𝒇::AbstractFeatureSet)  = getname.(𝒇)
getkeywords(𝒇::AbstractFeatureSet)  = getkeywords.(𝒇)
getdescriptions(𝒇::AbstractFeatureSet)  = getdescription.(𝒇)
export getfeatures, getmethods, getnames, getkeywords, getdescriptions

size(𝒇::AbstractFeatureSet) = size(getfeatures(𝒇))

getindex(𝒇::AbstractFeatureSet, i::Int) = getfeatures(𝒇)[i]
getindex(𝒇::AbstractFeatureSet, I) = FeatureSet(getfeatures(𝒇)[I])

function getindex(𝒇::AbstractFeatureSet, 𝐟::Vector{Symbol})
    i = [findfirst(x -> x == f, getnames(𝒇)) for f ∈ 𝐟]
    getindex(𝒇, i)
end

function getindex(𝒇::AbstractFeatureSet, f::Symbol)
    i = findfirst(x -> x == f, getnames(𝒇))
    getindex(𝒇, i)
end

function setindex!(𝒇::AbstractFeatureSet, f, i::Int)
    setindex!(𝒇.features, f, i)
    ()
end

IndexStyle(::AbstractFeatureSet) = IndexLinear()
eltype(::AbstractFeatureSet) = AbstractFeature

similar(::AbstractFeatureSet, ::Type{S}, dims::Dims) where {S} = FeatureSet(Vector{AbstractFeature}(undef, dims[1]))

deleteat!(𝒇::AbstractFeatureSet, args...) = deleteat!(𝒇.features, args...)

filter(f, 𝒇::AbstractFeatureSet) = FeatureSet(filter(f, getfeatures(𝒇)))

(+)(𝒇::AbstractFeatureSet, 𝒇′::AbstractFeatureSet) = FeatureSet(
                    [vcat(g(𝒇), g(𝒇′)) for g ∈ [getfeatures,
                                                getnames,
                                                getkeywords,
                                                getdescriptions]]...)
(\)(𝒇::AbstractFeatureSet, 𝒇′::AbstractFeatureSet) = setdiff(𝒇, 𝒇′)

# Allow operations between FeatureSet and Feature by converting the Feature
for p ∈ [:+, :\, :union, :intersect]
    eval(quote
        ($p)(𝒇::AbstractFeatureSet, f::AbstractFeature) = ($p)(𝒇, FeatureSet(f))
        ($p)(f::AbstractFeature, 𝒇::AbstractFeatureSet) = ($p)(FeatureSet(f), 𝒇)
    end)
end

(𝒇::AbstractFeatureSet)(x::AbstractVector) = FeatureVector([𝑓(x) for 𝑓 ∈ 𝒇], 𝒇)

function (𝒇::AbstractFeatureSet)(X::AbstractArray)
    F = Array{Float64}(undef, (length(𝒇), size(X)[2:end]...))
    threadlog = 0
    threadmax = prod(size(F)[2:end])/Threads.nthreads()
    @withprogress name="catch22" begin
        Threads.@threads for i ∈ CartesianIndices(size(F)[2:end])
            F[:, Tuple(i)...] = vec(𝒇(X[:, Tuple(i)...]))
            Threads.threadid() == 1 && (threadlog += 1)%50 == 0 && @logprogress threadlog/threadmax
        end
    end
    FeatureArray(F, 𝒇)
end

(𝒇::AbstractFeatureSet)(X::AbstractDimArray) = FeatureArray(𝒇(Array(X)), (Dim{:feature}(getnames(𝒇)), dims(X)[2:end]...))

(𝒇::AbstractFeatureSet)(x, f::Symbol) = 𝒇[f](x)


text:
using DiffEqFlux, Flux, OrdinaryDiffEq, Test, Optim, DiffEqSensitivity

x = Float32[0.8; 0.8]
tspan = (0.0f0,10.0f0)

ann = Chain(Dense(2,10,tanh), Dense(10,1))
p = Float32[-2.0,1.1]
p2,re = Flux.destructure(ann)
_p = [p;p2]
θ = [x;_p]

function dudt2_(u,p,t)
    x, y = u
    [(re(p[3:end])(u)[1]),p[1]*y + p[2]*x]
end

prob = ODEProblem(dudt2_,x,tspan,_p)
solve(prob,Tsit5())

function predict_rd(θ)
  Array(solve(prob,Tsit5(),u0=θ[1:2],p=θ[3:end],abstol=1e-7,reltol=1e-5,sensealg=TrackerAdjoint()))
end
loss_rd(p) = sum(abs2,x-1 for x in predict_rd(p))
l = loss_rd(θ)

cb = function (θ,l)
  println(l)
  #display(plot(solve(remake(prob,u0=Flux.data(_x),p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6)))
  false
end

# Display the ODE with the current parameter values.
cb(θ,l)

loss1 = loss_rd(θ)
res = DiffEqFlux.sciml_train(loss_rd, θ, BFGS(initial_stepnorm = 0.01), cb = cb)
loss2 = res.minimum
@test 10loss2 < loss1

## Partial Neural Adjoint

u0 = Float32[0.8; 0.8]
tspan = (0.0f0,25.0f0)

ann = Chain(Dense(2,10,tanh), Dense(10,1))

p1,re = Flux.destructure(ann)
p2 = Float32[-2.0,1.1]
p3 = [p1;p2]
θ = [u0;p3]

function dudt_(du,u,p,t)
    x, y = u
    du[1] = re(p[1:41])(u)[1]
    du[2] = p[end-1]*y + p[end]*x
end
prob = ODEProblem(dudt_,u0,tspan,p3)
solve(prob,Tsit5(),abstol=1e-8,reltol=1e-6)

function predict_adjoint(θ)
  Array(solve(prob,Tsit5(),u0=θ[1:2],p=θ[3:end],saveat=0.0:1:25.0))
end
loss_adjoint(θ) = sum(abs2,x-1 for x in predict_adjoint(θ))
l = loss_adjoint(θ)

cb = function (θ,l)
  println(l)
  #display(plot(solve(remake(prob,p=Flux.data(p3),u0=Flux.data(u0)),Tsit5(),saveat=0.1),ylim=(0,6)))
  false
end

# Display the ODE with the current parameter values.
cb(θ,l)

loss1 = loss_adjoint(θ)
res1 = DiffEqFlux.sciml_train(loss_adjoint, θ, ADAM(0.01), cb = cb, maxiters = 100)
res = DiffEqFlux.sciml_train(loss_adjoint, res1.minimizer, BFGS(initial_stepnorm = 0.01), cb = cb)
loss2 = res.minimum
@test 10loss2 < loss1


text:
module DatasetsCF

# package code goes here
using Persa
using MLBase: labelencode, labelmap
using DataFrames
using BinDeps
using CSV

include("movielens.jl")

end # module


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function testtable(fname)
    currentwidth = 800
    currentheight = 800
    Drawing(currentwidth, currentheight, fname)
    origin()
    background("white")
    setline(0.5)

    # test single cell tables
    t = Table(1, 1, 0, 0)
    for n in t
        pt = n[1]
        i = n[2]
    end

    t = Table(
        collect(range(30, stop=85, length=10)),
        collect(range(20, stop=75, length=10)))
    for n in t
        pt = n[1]
        i = n[2]
    end

    sethue("red")
    circle.(t[:, 3], 10, :fill)
    sethue("blue")
    circle.(t[3, :], 9, :fill)
    sethue("green")
    circle(t[4, 5], 8, :fill)

    for n in 1:20
        randomhue()
        circle.(t[rand(1:10), :], 7, :fill)
    end
    for n in 1:30
        randomhue()
        circle.(t[:, rand(1:10)], 6, :fill)
    end
    setline(2)
    highlightcells(t, 1:3:length(t), offset=2)

    sethue("black")
    fontsize(10)
    for (pt, n) in t
        label(string(n), :n, t[n], offset=15)
    end

    @test finish() == true
end

fname = "table.png"
testtable(fname)
println("table test: output in $(fname)")
println("...finished table tests: output in $(fname)")


text:
using Printf

const STREAM_ARRAY_SIZE = 1_000_000_000
const OFFSET = 0
const NTIMES = 100
const scalar = 3.0
const label = ["Copy:\t", "Scale:\t", "Add:\t", "Triad:\t"]
const bytes = [
    2 * sizeof(Float64) * STREAM_ARRAY_SIZE,
    2 * sizeof(Float64) * STREAM_ARRAY_SIZE,
    3 * sizeof(Float64) * STREAM_ARRAY_SIZE,
    3 * sizeof(Float64) * STREAM_ARRAY_SIZE
]
const sep = "----------------------------------------------"

# only print 3 decimal places per default
Base.show(io::IO, f::Float64) = @printf(io, "%.3f", f)
# 1. setup data + arrays
# (2. give overview over expected time )
# 3. do benchmark
#   1. copy STREAM_ARRAY_SIZE values from array to array
#   2. do scalar multiplication on array
#   3. add STREAM_ARRAY_SIZE values from two arrays together and store in third
#   4. do triad: for all STREAM_ARRAY_SIZE values: a[j] = b[j] + scalar * c[j]
# 4. Calculate avg, min, max times
# 5. Print best bandwidth and times for each catagory

function main()
    # initialization
    a = Array{Float64, 1}(undef, STREAM_ARRAY_SIZE + OFFSET)
    b = Array{Float64, 1}(undef, STREAM_ARRAY_SIZE + OFFSET)
    c = Array{Float64, 1}(undef, STREAM_ARRAY_SIZE + OFFSET)

    arraySize = sizeof(Float64) * STREAM_ARRAY_SIZE / 1024 / 1024

    println("STREAM.jl")
    println(sep)
    println("Array size = ", STREAM_ARRAY_SIZE, " (elements) Offset = ", OFFSET, " (elements)")
    println("Memory per array = ", arraySize, " MiB (= ", arraySize / 1024, " GiB)")
    println("Total memory = ", arraySize * 3, " MiB")
    println("Each kernel will be executed ", NTIMES, " times")
    println(" The *best* time for each kernel (excluding the first run)")
    println(" will be used to compute the reported bandwidth.")
    println(sep)

    Threads.@threads for i in 1:STREAM_ARRAY_SIZE
        a[i] = 1.0
        b[i] = 2.0
        c[i] = 0.0
    end

    times = zeros(Float64, 4, NTIMES)

    avgtime = zeros(4)
    maxtime = zeros(4)
    mintime = fill(typemax(Float64), 4)

    println("Using ", Threads.nthreads(), " threads")

    if ((q = quantum()) >= 1)
        println("Your clock granularity/precision appears to be ", q, "ns")
    else
        println("Your clock granularity appears to be less than 1ns")
        q = 1
    end

    # part of the timer code
    tmp = time_ns()
    Threads.@threads for i in 1:STREAM_ARRAY_SIZE
        a[i] = 2.0 * a[i]
    end
    tmp = (time_ns() - tmp) * 1e-3

    println("Each test below will take on the order of ", tmp, " microseconds")
    println("    (= ", floor(tmp/q), " clock ticks)")
    println("Increase the size of the arrays if this shows that")
    println(" you are not getting at least 20,000 clock ticks per test.")
    println(sep)

    for k in 1:NTIMES
	    times[1,k] = time_ns()
        Threads.@threads for j in 1:STREAM_ARRAY_SIZE
            c[j] = a[j]
        end
        times[1,k] = time_ns() - times[1,k]

	    times[2,k] = time_ns()
        Threads.@threads for j in 1:STREAM_ARRAY_SIZE
            b[j] = scalar * c[j]
        end
        times[2,k] = time_ns() - times[2,k]

	    times[3,k] = time_ns()
	    Threads.@threads for j in 1:STREAM_ARRAY_SIZE
            c[j] = a[j] + b[j]
        end
        times[3,k] = time_ns() - times[3,k]

        times[4,k] = time_ns()
	    Threads.@threads for j in 1:STREAM_ARRAY_SIZE
            a[j] = b[j] + scalar * c[j]
        end
        times[4,k] = time_ns() - times[4,k]
	    #tmp = @timed Threads.@threads for j in 1:STREAM_ARRAY_SIZE
        #    c[j] = a[j]
        #end
        #times[1,k] = tmp[2]
    end

    # Summarize the times
    # the first iteration is only for warmup
    for k in 2:NTIMES
        for j in 1:4
            # times in nanoseconds -> seconds
            times[j,k] = times[j,k] * 1e-9

            avgtime[j] = avgtime[j] + times[j,k]
            mintime[j] = min(mintime[j], times[j,k])
            maxtime[j] = max(maxtime[j], times[j,k])

            #tmp = @timed Threads.@threads for j in 1:STREAM_ARRAY_SIZE
            #    b[j] = scalar * c[j]
            #end
            #times[2,k] = tmp[2]

            #tmp = @timed Threads.@threads for j in 1:STREAM_ARRAY_SIZE
            #    c[j] = a[j] + b[j]
            #end
            #times[3,k] = tmp[2]

            #tmp = @timed Threads.@threads for j in 1:STREAM_ARRAY_SIZE
            #    a[j] = b[j] + scalar * c[j]
            #end
            #times[4,k] = tmp[2]
        end
    end

    println("Function   Best Rate MB/s  Avg time    Min time    Max time")
    for j in 1:4
        avgtime[j] = avgtime[j]/(NTIMES - 1)
        Printf.@printf "%s%12.1f  %11.6f  %11.6f  %11.6f\n" label[j] 1.0e-6 * bytes[j]/mintime[j] avgtime[j] mintime[j] maxtime[j]
    end
    checkResults(a, b, c)

	open("jl.results", "w") do f
        	# print all times to disk
	        Printf.@printf f "JCopy;JScale;JAdd;JTriad\n"
		for k in 2:NTIMES
			Printf.@printf f "%11.6f;%11.6f;%11.6f;%11.6f\n" times[1,k] times[2,k] times[3,k] times[4,k]
		end
	end
end

function quantum()
    M = 20
    timesfound = Array{Float64, 1}(undef, M)
    t1::UInt64 = 0
    t2::UInt64 = 0

    for i in 1:M
	t1 = time_ns()
        while ((t2 = time_ns()) - t1) < 1e3
	end
	timesfound[i] = (t1 = t2)
    end
    minDelta = 1000000000
    for i in 2:M
	delta = floor(timesfound[i] - timesfound[i-1])
	minDelta = min(minDelta, max(delta, 0))
    end
    return minDelta
end

function checkResults(a::Array{Float64,1}, b::Array{Float64,1}, c::Array{Float64,1})
    aj::Float64 = 1.0
    bj::Float64 = 2.0
    cj::Float64 = 0.0

    # part of timing
    aj = 2.0 * aj

    for k in 1:NTIMES
        cj = aj
        bj = scalar * cj
        cj = aj + bj
        aj = bj + scalar * cj
    end

    Printf.@printf "aj = %lf, bj = %lf, cj = %lf\n" aj bj cj

    aSumErr = 0.0
    bSumErr = 0.0
    cSumErr = 0.0
    for j in 1:STREAM_ARRAY_SIZE
        aSumErr += abs(a[j] - aj)
        bSumErr += abs(b[j] - bj)
        cSumErr += abs(c[j] - cj)
    end

    aAvgErr = aSumErr / STREAM_ARRAY_SIZE
    bAvgErr = bSumErr / STREAM_ARRAY_SIZE
    cAvgErr = cSumErr / STREAM_ARRAY_SIZE

    epsilon = 1.e-13

    err = 0
    if (abs(aAvgErr/aj) > epsilon)
        err += 1
        println("Validation of array a failed with rate: ", abs(aAvgErr/aj))
    end
    if (abs(bAvgErr/bj) > epsilon)
        err += 1
        println("Validation of array b failed with rate: ", abs(bAvgErr/bj))
    end
    if (abs(cAvgErr/cj) > epsilon)
        err += 1
        println("Validation of array c failed with rate: ", abs(cAvgErr/cj))
    end
    if (err == 0)
        Printf.@printf "Solution Validates: avg error less than %e on all three arrays\n" epsilon
    end
end

main()


text:
using CSV,DataFrames,Statistics,StatsBase
using Random
using ScottishTaxBenefitModel
using .Utils: coarse_match

#
# Scripts for creating a merged FRS (recipient) and SHS (donor) datase
# see the following:
#
# * match routine `` in [.Utils](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/blob/master/src/Utils.jl)
# * [blog post](https://stb-blog.virtual-worlds.scot/articles/2021/02/07/matching-notes.html) (includes some references);
# * test cases [matching_tests.jl]()https://github.com/grahamstark/ScottishTaxBenefitModel.jl/blob/master/test/matching_tests.jl
#
const DIR="/mnt/data/"
const SCOTLAND=299999999
const MAX_REPEATS = 200
#
# useful lists for searching for cases with no matches
#
const bestmatches = [:sernum, :datayear, :shelter_1,:tenure_1,:acctype_1,:singlepar_1,:numadults_1,:numkids_1,:empstathigh_1,:sochigh_1,:agehigh_1,:ethnichigh_1]
const worstmatches = [:sernum, :datayear, :shelter_3,:tenure_3,:acctype_3,:singlepar_3,:numadults_3,:numkids_3,:empstathigh_3,:sochigh_3,:agehigh_3,:ethnichigh_3]
#
# things that *have* to match
const critmatches = [:sernum, :datayear, :shelter_3,:singlepar_3,:numadults_2,:numkids_2,:empstathigh_3,:agehigh_3]
const critmatche1 = [:sernum, :datayear, :shelter_1,:singlepar_1,:numadults_2,:numkids_2,:empstathigh_1,:agehigh_1]

function init_council_counter( shs_councils :: DataFrame ) :: Dict{String,Float64}
    counts = Dict()
    for c in eachrow( shs_councils )
        counts[c.Code] = 0.0
    end
    counts
end

"""
Add a matched council to each model frs household record
this is chosen randomly from all the matched records, using
a the calculated inverse sample frequency of that council and
the quality of the match (high q's -> more likely, obvs.)
"""
function add_council_to_frs!(
    ;
    mhh   :: DataFrame,
    recip :: DataFrame,
    shs_councils :: DataFrame,
    inv_freqs :: Dict )

    n_rows = size( mhh )[1]
    mhh[!,:council] = fill("  ", n_rows)
    mhh[!,:nhs_board] = fill("  ", n_rows)

    counts =  init_council_counter( shs_councils )
    println( "counts $counts ")
    for r in eachrow(recip)
        println( "on r" )
        cnhs = []
        weights = []
        println( "r.sernum $(r.sernum) datayear $(r.datayear)" )
        dy = r.datayear+2000
        mpp = ((mhh.hid .== r.sernum).&(mhh.data_year .== dy ))
        mppsize = count( mpp )
        @assert mppsize in [0,1]
        println( "mppsize $mppsize" )
        if mppsize == 1
            # mra = mra[1,:] # convert output frs to row
            for n in 1:MAX_REPEATS
                idkey = Symbol( "shs_uniqidnew_$(n)" )
                ykey = Symbol( "shs_datayear_$(n)" )
                qkey = qkey = Symbol( "shs_quality_$(n)" )
                ckey = Symbol( "shs_council_$(n)" )
                nkey = Symbol( "shs_nhs_board_$(n)" )
                y = r[ykey]
                id = r[idkey]
                qual = r[qkey]
                nhs = r[nkey]
                ccode = r[ckey]
                # println( "qual=$qual" )
                if r[ykey] == 0 # no more match candidates
                    break
                end
                weight = inv_freqs[ccode]/qual
                # println( "got weight $weight, ccode $ccode " )
                push!( weights,  weight )
                push!( cnhs, (council=ccode,nhs=nhs))
            end # best to worst matches
            probs = ProbabilityWeights(weights/sum(weights)) # to a probability
            println("probs $(probs)" )
            println( "councils $(cnhs)" )
            ccode = sample( cnhs, probs )
            counts[ccode.council] += 1
            mhh[mpp,:council] .= ccode.council
            mhh[mpp,:nhs_board] .= ccode.nhs
            println( "selected $ccode" )
        end # frs model record exists
    end # loop round recipients
    return counts
end

function count_councils( shs_all_years :: DataFrame, shs_councils  :: DataFrame ) :: Dict
    ccounts = Dict()
    for c in eachrow( shs_councils )
        ccounts[c.Code] = 0.0
    end
    println( ccounts )
    for s in eachrow(shs_all_years)
        if ! ismissing(s.council)
            cc = shs_councils[(shs_councils.SHS_ID .== s.council),:Code]
            cc = cc[1]
            println(cc)
            ccounts[cc] += 1
        end
    end
    return ccounts
end

"""
Retrofit council and nhs region keys to the recipient dataset
"""
function add_in_las_to_recip!(
    recip         :: DataFrame,
    shs_all_years :: DataFrame,
    shs_councils  :: DataFrame )
    n_rows = size( recip )[1]
    for i in 1:MAX_REPEATS
       ckey = Symbol( "shs_council_$(i)" )
       nkey = Symbol( "shs_nhs_board_$(i)" )
       recip[!,ckey] = fill("  ",n_rows)
       recip[!,nkey] = fill("  ",n_rows)
    end
    for r in eachrow( recip )
        println( "on row $(r.sernum) $(r.datayear)" )
        for n in 1:MAX_REPEATS
            idkey = Symbol( "shs_uniqidnew_$(n)" )
            ykey = Symbol( "shs_datayear_$(n)" )
            ckey = Symbol( "shs_council_$(n)" )
            nkey = Symbol( "shs_nhs_board_$(n)" )
            y = r[ykey]
            id = r[idkey]
            if y > 0 # a match
                shsrow=shs_all_years[((shs_all_years.uniqidnew.==id).&(shs_all_years.datayear.==y)),:][1,:]
                cc = shs_councils[(shs_councils.SHS_ID .== shsrow.council),:Code][1]
                r[ckey] = cc
                r[nkey] = shsrow.hlthbd2014
            end
        end
    end
end

function make_highest_hh_income_marker( hhad :: DataFrame ) :: Vector{Bool}
	sns = unique(hhad.sernum)
	n = size( hhad )[1]
	highest_income = fill( false, n )
	pers_count = 0
	for sn in sns
		ads = hhad[ (hhad.sernum .== sn),:]
		max_inc_in_hh = -9999999999999999.99
		richest_in_hh = -1   # pos in this hh of richest member
		richest_array_count = -1 # pos in the overall array of richest hh member
		adno = 0
		for ad in eachrow(ads)
			adno += 1
			pers_count += 1
			inc = ad.indinc === missing ? 0.0 : Float64( ad.indinc )
			if inc > max_inc_in_hh
				max_inc_in_hh = inc
				richest_array_count = pers_count
				richest_in_hh = adno
			end
		end
		if richest_in_hh > 2
			println( "sernum=$sn max_inc_in_hh=$max_inc_in_hh richest_in_hh=$richest_in_hh " )
		end
		@assert richest_in_hh > 0
		highest_income[richest_array_count] = true
	end
	return highest_income
end

function loadshs( year::Int )::DataFrame
	year -= 2000
	ystr = "$(year)$(year+1)"
	fname = "$(DIR)/shs/$(ystr)/tab/shs20$(year)_social_public.tab"
	println( "loading '$fname'" )
	shs = CSV.File( fname;
	    missingstrings=["NA",""],
	    types=Dict(:uniqidnew => String,:UNIQIDNEW=>String)) |> DataFrame
	lcnames = Symbol.(lowercase.(string.(names(shs))))
    rename!(shs,lcnames)
    shs[!,:datayear] .= year
    return shs
end

function loadfrs( year::Int, fname :: String ) :: DataFrame
	ystr = "$(year)$(year+1)"
	fname = "$(DIR)/frs/$(year)/tab/$(fname).tab"
	println( "loading '$fname'" )
	frs = CSV.File( fname; missingstrings=["-1"] ) |> DataFrame
	lcnames = Symbol.(lowercase.(string.(names(frs))))
    rename!(frs,lcnames)
    frs[!,:datayear] .= year-2000
    return frs
end

"""
Load FRS hhld datasets, and then make Scotland only subsets. with just the adult with highest income
merged in.
FIXME years for FRS 2015:2019, etc; years for
"""
function create_frs( years :: UnitRange ) :: DataFrame
    n = length(years)
    frs=Array{DataFrame}(undef,n)
    i = 0
    for year in years
        i += 1
        hh = loadfrs( year, "househol" )
        shh = hh[(hh.gvtregn .== SCOTLAND),:]
        ad = loadfrs( year, "adult" )
        shhad=innerjoin( ad, shh, on=:sernum; makeunique=true )
        highest_income = make_highest_hh_income_marker( shhad )
        frs[i] = shhad[highest_income,:]
    end
    return vcat( frs...; cols=:intersect )
end

"""
Stack scottish household surveys.
"""
function create_shs( years :: UnitRange ) :: DataFrame
    n = length(years)
    shs = Array{DataFrame}(undef,n)
    i = 0
    for year in years
        i += 1
        shs[i] = loadshs(year)
    end
    return vcat( shs...; cols=:intersect )
end

# SHS Tenure
# Pos. = 56	Variable = tenure	Variable label = Tenure - SHS but non-harmonised version
# This variable is    numeric, the SPSS measurement level is NOMINAL
# 	Value label information for tenure
# 	Value = 1.0	Label = Owned outright
# 	Value = 2.0	Label = Buying with help of loan/mortgage
# 	Value = 3.0	Label = Rent – LA
# 	Value = 4.0	Label = Rent - HA, Co-op
# 	Value = 5.0	Label = Rent - private landlord
# 	Value = 6.0	Label = Other
# 	Value = 999998.0	Label = Don't know
# 	Value = 999999.0	Label = Refused
#
# # FRS TENURE
#
# os. = 258	Variable = TENTYP2	Variable label = Tenure type
# This variable is    numeric, the SPSS measurement level is NOMINAL
# SPSS user missing values = -9.0 thru -1.0
# 	Value label information for TENTYP2
# 	Value = 1.0	Label = LA / New Town / NIHE / Council rented
# 	Value = 2.0	Label = Housing Association / Co-Op / Trust rented
# 	Value = 3.0	Label = Other private rented unfurnished
# 	Value = 4.0	Label = Other private rented furnished
# 	Value = 5.0	Label = Owned with a mortgage (includes part rent / part own)
# 	Value = 6.0	Label = Owned outright
# 	Value = 7.0	Label = Rent-free
# 	Value = 8.0	Label = Squats

# harmonised tenure
#
# 1 -> OO
# 2 -> Mortgaged
# 3 -> LA/Council Rented
# 4 -> HA Rented
# 5 -> Private Rented
# 6 -> Other

# coarsened tenure
# 1 -> Owned
# 2 -> Rented
# 3 -> Other

# _3
# 1 -> Any

function shs_tenuremap( tenure :: Union{Int,Missing} ) :: Vector{Int}
    out = fill( 0, 3 )
    if ismissing( tenure ) # || tenure >= 5
        return out;
    end
    out[3] = 1
    if tenure == 1 # OO
        out[1] = 1
        out[2] = 1
    elseif tenure == 2
        out[1] = 2
        out[2] = 1
    elseif tenure == 3
        out[1] = 3
        out[2] = 2
    elseif tenure == 4
        out[1] = 4
        out[2] = 2
    elseif tenure == 5
        out[1] = 5
        out[2] = 2
    elseif tenure in [6, 7, 999998, 999998] # 7 is poss coding error
        out[1] = 6
        out[2] = 3
    else
        @assert false "unmatched tenure $tenure";
    end
    out[3] = out[2]
    return out
end

#
# level 1 -> actual number of people
# level 2
#   0 -> 0 (for kids)
#   1 -> 1
#   2 -> 2
#   3 -> 3:5
#   4 -> > 5
# level 3
#  1 adult -> 0
#  2 adults -> 1
#  > 2 adults -> 2
#  0 child -> 0
#  > 0 child -> 1
#
function total_people( n :: Union{Int,Missing}, def :: Int, is_child :: Bool ) :: Vector{Int}
    out = fill( def, 3 )
    if ismissing( n )
        return out
    end
    out[1] = n
    if n == 0
        out[2] = 0
    elseif n == 1
        out[2] = 1
    elseif n == 2
        out[2] = 2
    elseif n in 3:5
        out[2] = 3
    else
        out[2] = 4
    end
    if is_child # any children
       out[3] = out[2] > 0 ? 1 : 0
    else
       @assert out[2] > 0 "no adults"
       if out[2] == 1
            out[3] == 0
       elseif out[2] == 2
            out[3] = 1
       else
           out[3] = 2
       end
    end
    return out
end

"""
1. age (max 80)
2. age 5 year bands
3. age 20 year bands
"""
function age( age  :: Union{Int,Missing} )  :: Vector{Int}
    out = fill( 0, 3 )
    if ismissing( age )
        return out
    end
    age = min( 80, age )
    out[1] = age
    out[2] = Int(trunc(age/5))
    out[3] = Int(trunc(age/20))
    return out
end

function frs_tenuremap( tentyp2 :: Union{Int,Missing} ) :: Vector{Int}
    out = fill( 0, 3 )
    out[3] = 1
    if ismissing( tentyp2 )
        out[1] = 6
        out[2] = 3
    elseif tentyp2 == 1
        out[1] = 3
        out[2] = 2
    elseif tentyp2 == 2
        out[1] = 4
        out[2] = 2
    elseif tentyp2 in [3,4]
        out[1] = 5
        out[2] = 2
    elseif tentyp2 == 5
        out[1] = 2
        out[2] = 1
    elseif tentyp2 == 6
        out[1] = 1
        out[2] = 1
    elseif tentyp2 in 7:8
        out[1] = 6
        out[2] = 3
    else
        @assert false "unmatched tentyp2 $tentyp2";
    end
    return out
end

function to_i( i :: Union{Int,Missing} ) :: Int
    if ismissing(i)
        return -1
    end
    return i
end

function is_sp(i::Integer)
    i in 9:11
end

function setone( i :: Union{Int,Missing}, decider::Function ) :: Vector{Int}
    i = to_i( i )
    o = decider(i)
    return fill( o, 3 )
end

function setone( i :: Union{Int,Missing}, target :: Int = 1 ) :: Vector{Int}
    i = to_i( i )
    o = (i == target )
    return fill( o, 3 )
end

#Pos. = 58	Variable = hb1	Variable label = hb1 - Is the household's accommodation...
# SHS Accomodation
#
# This variable is  numeric, the SPSS measurement level is NOMINAL
# SPSS user missing values = -1.0 thru None
# 	Value label information for hb1
# 	Value = 1.0	Label = House or bungalow
# 	Value = 2.0	Label = A flat, maisonette or apartment (including )
# 	Value = 3.0	Label = Other, including room(s), caravan/mobile homes
#
# Pos. = 59	Variable = hb2	Variable label = hb2 - Is it...
# This variable is  numeric, the SPSS measurement level is SCALE
# SPSS user missing values = -1.0 thru None
# 	Value label information for hb2
# 	Value = 1.0	Label = Detached
# 	Value = 2.0	Label = Semi-detached
# 	Value = 3.0	Label = or terraced/end of terrace?
#
#
#
# FRS TYPEACC
#  1     | Whole house/bungalow, detached
#  2     | Whole house/bungalow, semi-detache
#  3     | Whole house/bungalow, terraced
#  4     | Purpose-built flat or maisonette
#  5     | Converted house/building
#  6     | Caravan/Mobile home or Houseboat
#  7     | Other

"""
level 1
1.     | Whole house/bungalow, detached
2.     | Whole house/bungalow, semi-detache
3.     | Whole house/bungalow, terraced
4.     | Purpose-built flat or maisonette/Converted house/building
5.     | Caravan/Mobile home or Houseboat, Other
level 2
1.     House
2.     Flat
3.     Other
level 3
1. any

"""
function frs_btype( typeacc :: Union{Missing,Int} ) :: Vector{Int}
    out = fill( 1, 3 )
    if ismissing(typeacc)
        return out
    end
    if typeacc < 5
        out[1] = typeacc
    elseif typeacc == 5
        out[1] = 4
    elseif typeacc in 6:7
        out[1] = 5
    else
        @assert false "typeac $typacc not recognised"
    end
    if typeacc in 1:3
        out[2] = 1
    elseif typeacc in 4:5
        out[2] = 2
    elseif typeacc in 6:7
        out[2] = 3
    end
    return out
end

function shs_btype( hb1 :: Union{Missing,Int}, hb2 :: Union{Missing,Int} ) :: Vector{Int}
    out = fill( -886, 3 )
    if ismissing(hb1)
        return fill( -991, 3 )
    end
    if hb1 == 1
        if hb2 == 1
            out[1] = 1
        elseif hb2 == 2
            out[1] = 2
        elseif hb2 == 3
            out[1] = 3
        else
            @assert false "hb2=$hb2"
        end
    elseif hb1 == 2
        out[1] = 4
    elseif hb1 == 3
        out[1] = 5
    else
        @assert false "unrecognised hb1 $hb1"
    end
    out[2] = hb1
    out[3]=out[2]
    return out
end

#

# shelter
#
# accsup1 shs 1=yes 2=no
# shelter frs 1=yes 2=no
#
# frs CTREB
# shs hh57_08 == 1 (? maybe 2 bus)
#
# hhsize shs
# totads shs   adulth frs
# totkids shs  /depchldh frs/comp
# totpeeps (hhsize?)
#
#
# single parent hhtype_new == 3  shs hhcomps in 9:11 frs
#
# accom type TYPEACC hb1

function assign!( df :: DataFrame, name :: Symbol, vals )
    n = size(vals[1])[1]
    m = size( df )[1]
    ET = eltype(vals[1])
    # println( vals )
    for i in 1:n
        sym = Symbol("$(String(name))_$(i)")
        println( "writing $sym $i" )
        df[!,sym] = zeros(ET,m)
        for j in 1:m
            df[j,sym] = vals[j][i]
        end
    end
end

#   SHS hihecon
# 	Value = 1.0	Label = A - Self employed
# 	Value = 2.0	Label = B - Employed full time
# 	Value = 3.0	Label = C - Employed part time
# 	Value = 4.0	Label = D - Looking after the home or family
# 	Value = 5.0	Label = E - Permanently retired from work
# 	Value = 6.0	Label = F - Unemployed and seeking work
# 	Value = 7.0	Label = G - At school
# 	Value = 8.0	Label = H - In further / higher education
# 	Value = 9.0	Label = I - Gov't work or training scheme
# 	Value = 10.0	Label = J - Permanently sick or disabled
# 	Value = 11.0	Label = K - Unable to work because of short-term illness or injury
# 	Value = 12.0	Label = L - Pre school / Not yet at school
# 	Value = 13.0	Label = Other (specify)
#
#
# 	FRS EMPSTATI
#   Value = 1.0	Label = Full-time Employee
# 	Value = 2.0	Label = Part-time Employee
# 	Value = 3.0	Label = Full-time Self-Employed
# 	Value = 4.0	Label = Part-time Self-Employed
# 	Value = 5.0	Label = Unemployed
# 	Value = 6.0	Label = Retired
# 	Value = 7.0	Label = Student
# 	Value = 8.0	Label = Looking after family/home
# 	Value = 9.0	Label = Permanently sick/disabled
# 	Value = 10.0	Label = Temporarily sick/injured
# 	Value = 11.0	Label = Other Inactive

# Level 1
# 1 Full-time Employee
# 2 Part-time Employee
# 3 Self Employed
# 4 Unemployed
# 5 Retired
# 6 Student
# 7 Looking after the home or family
# 8 Permanently sick/disabled
# 9 Temporarily sick/injured
# 10 Other Inactive

# level 2
# 1 employed
# 2 self-employed
# 3 retired
# 4 student
# 5 other, inc unemployed

# level3
# 1 working
# 2 retired
# 3 not working

function frs_empstat( empstati :: Int ) :: Vector{Int}
    out = fill( 0, 3 )
    out[1] = empstati
    if empstati > 3
        out[1] -= 1
    end
    if out[1] in 1:2
        out[2] = 1
    elseif out[1] in 3
        out[2] = 2
    elseif out[1] in 5
        out[2] = 3
    elseif out[1] in 6
        out[2] = 4
    elseif out[1] in [4,7,8,9,10]
        out[2] = 5
    else
        @assert false "empstati $empstati"
    end
    if out[1] in 1:3
        out[3] = 1
    elseif out[1] in 5
        out[3] = 2
    else
        out[3] = 3
    end
    return out
end

function shs_empstat( hihecon :: Union{Missing,Int} ) :: Vector{Int}
    out = fill( 0, 3 )
    if ismissing(hihecon) # value 14 not documented
        return fill( -990, 3 )
    end
    if hihecon > 13
        return fill( -989, 3 )
    end
    if  hihecon == 1
        out[1] = 3
    elseif hihecon == 2
        out[1] = 1
    elseif hihecon == 3
        out[1] = 2
    elseif hihecon == 4
        out[1] = 7
    elseif hihecon == 5
        out[1] = 5
    elseif hihecon == 6
        out[1] = 4
    elseif hihecon in [7,9,13]
        out[1] = 10
    elseif hihecon == 8
        out[1] = 6
    elseif hihecon == 10
        out[1] = 8
    elseif hihecon == 11
        out[1] = 9
    else
        @assert false "hihecon $hihecon"
    end
    if out[1] in 1:2
        out[2] = 1
    elseif out[1] in 3
        out[2] = 2
    elseif out[1] in 5
        out[2] = 3
    elseif out[1] in 6
        out[2] = 4
    elseif out[1] in [4,7,8,9,10]
        out[2] = 5
    else
        @assert false "out[1] $(out[1])"
    end
    if out[1] in 1:3
        out[3] = 1
    elseif out[1] in 4
        out[3] = 2
    else
        out[3] = 3
    end
    return out
end

#
#
# FRS ethgr3
#
# Pos. = 512	Variable = ETHGR3	Variable label = Ethnicity of Adult (harmonised version)
# This variable is  numeric, the SPSS measurement level is NOMINAL
# SPSS user missing values = -9.0 thru -1.0
# 	Value label information for ETHGR3
# 	Value = 1.0	Label = White
# 	Value = 2.0	Label = Mixed/ Multiple ethnic groups
# 	Value = 3.0	Label = Asian/ Asian British
# 	Value = 4.0	Label = Black/ African/ Caribbean/ Black British
# 	Value = 5.0	Label = Other ethnic group
#
# SHS  hih_eth2012
#
# Pos. = 35	Variable = HIH_ETH2012	Variable label = Ethnic origin of HiH
# This variable is  numeric, the SPSS measurement level is SCALE
# 	Value label information for HIH_ETH2012
# 	Value = 1.0	Label = White
# 	Value = 2.0	Label = Minority ethnic groups
# 	Value = 3.0	Label = Don't know
# 	Value = 4.0	Label = Refused

#
# out
# level1
# 1-> White
# 2-> Non White
# level2
# 1-> Any
# level3
# 1-> Any

function frs_ethnic( ethgr3 :: Int ) :: Vector{Int}
    out = fill( 1, 3 )
    out[1] = ethgr3 == 1 ? 1 : 2
    return out
end

function shs_ethnic( hih_eth2012 :: Union{Missing,Int} ) :: Vector{Int}
    out = fill( 1, 3 )
    if ismissing(hih_eth2012) # value 14 not documented
        return fill( -986, 3 )
    end
    if hih_eth2012 > 2
        return out
    end
    out[1] = hih_eth2012 == 1 ? 1 : 2
    return out
end

#
# annual hh net income
#  mean(collect(skipmissing(shs_all_years.annetinc)))
# SHS mean: 27,046.418224914407 median 23,000.0
# mean(collect(skipmissing(frs_all_years_scot_he.hhinc*52)))
# FRS mean 37,704.65 median 28,184.0

# SHS SOC hihsoc
# Pos. = 1410	Variable = hihsoc	Variable label = HIH Social Occupational Classification
# This variable is  numeric, the SPSS measurement level is NOMINAL
# 	Value label information for hihsoc
# 	Value = 1.0	Label = MANAGERS, DIRECTORS AND SENIOR OFFICIALS
# 	Value = 2.0	Label = PROFESSIONAL OCCUPATIONS
# 	Value = 3.0	Label = ASSOCIATE PROFESSIONAL AND TECHNICAL OCCUPATIONS
# 	Value = 4.0	Label = ADMINISTRATIVE AND SECRETARIAL OCCUPATIONS
# 	Value = 5.0	Label = SKILLED TRADES OCCUPATIONS
# 	Value = 6.0	Label = CARING, LEISURE AND OTHER SERVICE OCCUPATIONS
# 	Value = 7.0	Label = SALES AND CUSTOMER SERVICE OCCUPATIONS
# 	Value = 8.0	Label = PROCESS, PLANT AND MACHINE OPERATIVES
# 	Value = 9.0	Label = ELEMENTARY OCCUPATIONS
# 	Value = -9.0	Label = NO INFORMATION
#
# FRS
# Pos. = 425	Variable = SOC2010	Variable label = Standard Occupational Classification
# This variable is  numeric, the SPSS measurement level is NOMINAL
# SPSS user missing values = -9.0 thru -1.0
# 	Value label information for SOC2010
# 	Value = 0.0	Label = Undefined
# 	Value = 1000.0	Label = Managers Directors & Senior Officials
# 	Value = 2000.0	Label = Professional Occupations
# 	Value = 3000.0	Label = Associate Prof. & Technical Occupations
# 	Value = 4000.0	Label = Admin & Secretarial Occupations
# 	Value = 5000.0	Label = Skilled Trades Occupations
# 	Value = 6000.0	Label = Caring leisure and other service occupations
# 	Value = 7000.0	Label = Sales & Customer Service
# 	Value = 8000.0	Label = Process, Plant & Machine Operatives
# 	Value = 9000.0	Label = Elementary Occupations
#
# level 1
# as hisoc if 1:9 else 0
# level2
# 0 -> undefined
# 1 -> 1..2
# 2 -> 3,4,7
# 3 -> 5,8
# 4 -> 6,9
# level2
# 1 -> all
# 11,000 undefined occs in shs, 263 in FRS
#
function map_social( soc :: Int ) :: Vector{Int}
    out = fill(0,3)
    if ! (soc in 1:9)
        return [0,0,1]
    end
    out[1] = soc
    if soc in 1:2
        out[2] = 1
    elseif soc in [3,4,7]
        out[2] = 2
    elseif soc in [5,8]
         out[2] = 3
    elseif soc in [6,9]
         out[2] = 4
    else
        @assert false "soc=$soc"
    end
    out[3]=1
    return out
end

function bedrooms( rooms :: Union{Missing,Int} ) :: Vector{Int}
    rooms = min(6, rooms )
    out = fill(0,3)
    if (ismissing(rooms) || (rooms == 0 ))
        return [0,0, 1]
    end
    out = fill(0,3)
    out[1] = rooms
    out[2] = min( rooms, 3)
    out[3] = rooms == 1 ? 1 : 2
    return out
end

function shs_map_social( hihsoc :: Union{Missing,Int} ) :: Vector{Int}
    if ismissing(hihsoc)
        return [0,0,1]
    end
    return map_social( hihsoc )
end

function frs_map_social( soc2010 )
    return map_social( Int(soc2010/1000))
end

function data_year( dy :: Int ) :: Vector{Int}
    return [dy,1,1]
end

function create_donor_and_recip_datasets( frs_all_years_scot_he :: DataFrame, shs_all_years::DataFrame )
    #
    # Create donor  (SHS) and recipient (FRS) datasets.
    #
    donor = DataFrame( datayear=shs_all_years.datayear, uniqidnew=shs_all_years.uniqidnew )
    recip = DataFrame( datayear=frs_all_years_scot_he.datayear, sernum=frs_all_years_scot_he.sernum )

    # add placeholders for N matches

    n_rows = size( recip )[1]
    for i in 1:MAX_REPEATS
       idkey = Symbol( "shs_uniqidnew_$(i)" )
       ykey = Symbol( "shs_datayear_$(i)" )
       qkey = Symbol( "shs_quality_$(i)" )
       recip[!,idkey] = fill("  ",n_rows)
       recip[!,ykey] = fill(0,n_rows)
       recip[!,qkey] = fill(0,n_rows)
    end

    assign!( recip, :shelter, setone.( frs_all_years_scot_he.shelter ))
    assign!( recip, :tenure, frs_tenuremap.(frs_all_years_scot_he.tentyp2))
    assign!( recip, :singlepar, setone.(frs_all_years_scot_he.hhcomps, is_sp))
    assign!( recip, :numadults, total_people.( frs_all_years_scot_he.adulth, -777, false ))
    assign!( recip, :numkids, total_people.( frs_all_years_scot_he.depchldh, -776, true ))
    assign!( recip, :acctype, frs_btype.( frs_all_years_scot_he.typeacc ))
    assign!( recip, :agehigh, age.( frs_all_years_scot_he.age80 ))
    assign!( recip, :empstathigh, frs_empstat.( frs_all_years_scot_he.empstati ))
    assign!( recip, :ethnichigh, frs_ethnic.( frs_all_years_scot_he.ethgr3 ))
    assign!( recip, :sochigh, frs_map_social.( frs_all_years_scot_he.soc2010 ))
    assign!( recip, :datayear, data_year.( frs_all_years_scot_he.datayear ))
    assign!( recip, :bedrooms, bedrooms.( frs_all_years_scot_he.bedroom6 ))

    assign!( donor, :shelter, setone.( shs_all_years.accsup1 ))
    assign!( donor, :tenure, shs_tenuremap.(shs_all_years.tenure))
    assign!( donor, :singlepar, setone.( shs_all_years.hhtype_new, 3 ))
    assign!( donor, :numadults, total_people.( shs_all_years.totads, -888, false ))

    #
    # fixme: check distinction between `totkids` (max 4 in the data) and `numkids`
    # nothing much in docs
    #
    assign!( donor, :numkids, total_people.( shs_all_years.numkids, -887, true ))
    assign!( donor, :acctype, shs_btype.( shs_all_years.hb1, shs_all_years.hb2))
    assign!( donor, :agehigh, age.( shs_all_years.hihage ))
    assign!( donor, :empstathigh, shs_empstat.( shs_all_years.hihecon ))
    assign!( donor, :ethnichigh, shs_ethnic.( shs_all_years.hih_eth2012 ))
    assign!( donor, :sochigh, shs_map_social.( shs_all_years.hihsoc ))
    assign!( donor, :datayear, data_year.( shs_all_years.datayear ))
    assign!( donor, :bedrooms, bedrooms.( shs_all_years.hc4 ))
    ( donor, recip )
end

function print_matches( matches )
    tm = sum( matches.matches)
    println( "mot matches= $tm" )
    cc=StatsBase.counts(matches.quality)
    for i in 1:size(cc)[1]
        if cc[i] !== 0
            println( "$(i-10) = $(cc[i])")
        end
    end
end

struct MatchRec
    quality   :: Int
    uniqidnew :: String
    datayear  :: Int
end

"""
Sort the matches by match quality and then randomise within blocks of matches with the
same match quality.
"""
function shuffle_blocks( a :: Vector ) :: Vector
    sort!( a, by = x -> x.quality )
    out = fill( a[1], 0 )
    block = fill( a[1], 0)
    n = size(a)[1]
    last = a[1]
    for i in 1:n
        if a[i].quality == last.quality
            push!( block, a[i] )
        else
            # found a new qual; shuffle the block, append it to output & start a new block
            last = a[i]
            shuffle!( block )
            out = vcat( out, block )
            resize!( block, 0 )
            push!( block, a[i] )
        end
    end
    # handle whatever's left
    if size(block)[1] > 0
        shuffle!( block )
        out = vcat( out, block )
    end
    return out
end

"""
Add values for columns `shs_uniqidnew_n`, `shs_datayear_n` and `shs_quality_n` for one row of
the recipient dataset, for up to MAX_REPEATS matches with best in 1, next in 2 and so on.
"""
function load_matches_to_recip!( r :: DataFrameRow, matches :: NamedTuple, donor :: DataFrame )
    num_matches = count(matches.matches)
    # matches is a list of all the donor records that match - get a count of trues
    println( "num_matches=$num_matches")
    if num_matches == 0
        println("Zero matches for frs $(r.sernum) $(r.datayear) ")
        return 0;
    end
    # println( "matches $(matches.matches) ")
    nq = size( matches.quality )[1]
    nd = size( donor )[1]
    @assert nq == nd " nq = $nq nd = $nd num_matches = $num_matches"
    v = fill(MatchRec(0,"",0),0)
    for i in 1:nq
        if matches.matches[i]
            dr = donor[i,:]
            push!( v, MatchRec( matches.quality[i], dr.uniqidnew, dr.datayear ))
        end
    end
    nv = size( v )
    # println( "size v = $nv ")
    vs = shuffle_blocks( v )
    svs = size( vs )[1]
    @assert num_matches == svs "num_matches=$num_matches svs = $svs "
    num_matches = min( MAX_REPEATS, num_matches )
    for i in 1:num_matches
       idkey = Symbol( "shs_uniqidnew_$(i)" )
       ykey = Symbol( "shs_datayear_$(i)" )
       qkey = Symbol( "shs_quality_$(i)" )
       r[idkey] = vs[i].uniqidnew
       r[ykey] = vs[i].datayear
       r[qkey] = vs[i].quality
       if i > 1
            @assert vs[i].quality >= vs[i-1].quality
       end
    end
    return num_matches
end

function match_up_unmatched!( recip :: DataFrame, donor :: DataFrame, fields :: Vector{Symbol}, targets :: Vector )
    n = size( targets )[1]
    nr = size(recip)[1]
    @assert n == nr " n = $n nr = $nr "
    s = 0
    for i in 1:n #iter
        if targets[i]
            s += 1
            rs = recip[i,:]
            matches = Utils.coarse_match(
                rs,
                donor,
                fields,
                3 )
            load_matches_to_recip!( rs, matches, donor )
       end
    end
    return s
end

function make_initial_match!( recip :: DataFrame, donor :: DataFrame )
    #
    # initial match across all FRS households.
    #
    targets = [:shelter,:tenure,:acctype,:bedrooms,:singlepar,:numadults,:numkids,:empstathigh,:sochigh,:agehigh,:ethnichigh] #  drop :datayear since year 15 doesn't overlap
    all = fill( true, size( recip )[1] ) # match every frs household
    match_up_unmatched!( recip, donor, targets, all )
end


text:
#ДАНО: Робот - рядом с горизонтальной перегородкой (под ней), бесконечно продолжающейся в обе стороны,
#в которой имеется проход шириной в одну клетку.
#РЕЗУЛЬТАТ: Робот - в клетке под проходом

function glavnaia(r::Robot)
    num_steps=1
    side=Ost
    while isborder(r,Nord)==true
        check_space(r,side,num_steps)
        side=inverse(side)
        num_steps+=1
    end
end

function check_space(r::Robot, side::HorizonSide, num_steps::Int) #двигает робота на num шагов и проверяет наличие прохода сверху
    for _ in 1:num_steps
        move!(r,side)
        if isborder(r,Nord)==false
            break
        end
    end
end

inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4)) # Возвращает направление, противоположное заданному

max_stars_repo_name: Aminulik/mirea-progs

max_issues_repo_name: Aminulik/mirea-progs

max_forks_repo_name: Aminulik/mirea-progs


text:
using FANN
using Base.Test

net = ANN([1, 1], [:linear])

# train on some data
X = [0.0 1.0]
Y = [0.0 1.0]
train!(net, DataSet(X, Y), max_epochs=10000, epochs_between_reports=0)

# save and load
savenet(net, "net")
net2 = loadnet("net")
rm("net")

# should be the same prediction
@test norm(predict(net, X) - predict(net2, X)) == 0.0


text:
"Parse PDDL domain description."
parse_domain(expr::Vector) =
    parse_domain(expr, GenericDomain)
parse_domain(str::AbstractString) =
    parse_domain(parse_string(str))
parse_domain(expr::Vector, domain_type::Type) =
    domain_type(parse_description(:domain, expr)...)
parse_domain(str::AbstractString, domain_type::Type) =
    parse_domain(parse_string(str), domain_type)
top_level_parsers[:domain] = parse_domain

"Parse domain requirements."
function parse_requirements(expr::Vector)
    reqs = Dict{Symbol,Bool}(e.name => true for e in expr[2:end])
    reqs = merge(DEFAULT_REQUIREMENTS, reqs)
    unchecked = [k for (k, v) in reqs if v == true]
    while length(unchecked) > 0
        req = pop!(unchecked)
        deps = get(IMPLIED_REQUIREMENTS, req, Symbol[])
        if length(deps) == 0 continue end
        reqs = merge(reqs, Dict{Symbol,Bool}(d => true for d in deps))
        append!(unchecked, deps)
    end
    return reqs
end
parse_requirements(expr::Nothing) = copy(DEFAULT_REQUIREMENTS)
head_field_parsers[:domain][:requirements] = parse_requirements

"Parse type hierarchy."
function parse_types(expr::Vector)
    @assert (expr[1].name == :types) ":types keyword is missing."
    types = Dict{Symbol,Vector{Symbol}}(:object => Symbol[])
    all_subtypes = Set{Symbol}()
    accum = Symbol[]
    is_supertype = false
    for e in expr[2:end]
        if e == :-
            is_supertype = true
            continue
        end
        subtypes = get!(types, e, Symbol[])
        if is_supertype
            append!(subtypes, accum)
            union!(all_subtypes, accum)
            accum = Symbol[]
            is_supertype = false
        else
            push!(accum, e)
        end
    end
    maxtypes = setdiff(keys(types), all_subtypes, [:object])
    append!(types[:object], collect(maxtypes))
    return (typetree=types,)
end
parse_types(expr::Nothing) = Dict{Symbol,Vector{Symbol}}(:object => Symbol[])
head_field_parsers[:domain][:types] = parse_types

"Parse constants in a planning domain."
function parse_constants(expr::Vector)
    @assert (expr[1].name == :constants) ":constants keyword is missing."
    objs, types = parse_typed_consts(expr[2:end])
    types = Dict{Const,Symbol}(o => t for (o, t) in zip(objs, types))
    return (constants=objs, constypes=types)
end
parse_constants(::Nothing) =
    (constants=Const[], constypes=Dict{Const,Symbol}())
head_field_parsers[:domain][:constants] = parse_constants

"Parse predicate list."
function parse_predicates(expr::Vector)
    @assert (expr[1].name == :predicates) ":predicates keyword is missing."
    preds = Dict{Symbol,Signature}()
    declarations, types = parse_typed_declarations(expr[2:end], :boolean)
    for ((term, argtypes), type) in zip(declarations, types)
        if type !== :boolean error("Predicate $term is not boolean.") end
        preds[term.name] = Signature(term, argtypes, type)
    end
    return preds
end
parse_predicates(::Nothing) = Dict{Symbol,Signature}()
head_field_parsers[:domain][:predicates] = parse_predicates

"Parse list of function (i.e. fluent) declarations."
function parse_functions(expr::Vector)
    @assert (expr[1].name == :functions) ":functions keyword is missing."
    funcs = Dict{Symbol,Signature}()
    declarations, types = parse_typed_declarations(expr[2:end], :numeric)
    for ((term, argtypes), type) in zip(declarations, types)
        funcs[term.name] = Signature(term, argtypes, type)
    end
    return funcs
end
parse_functions(::Nothing) = Dict{Symbol,Signature}()
head_field_parsers[:domain][:functions] = parse_functions

"Parse axioms (a.k.a. derived predicates)."
function parse_axiom(expr::Vector)
    @assert (expr[1].name in [:axiom, :derived]) ":derived keyword is missing."
    head = parse_formula(expr[2])
    body = parse_formula(expr[3])
    return Clause(head, Term[body])
end
body_field_parsers[:domain][:axiom] = parse_axiom

"Parse axioms (a.k.a. derived predicates)."
parse_derived(expr::Vector) = parse_axiom(expr)
body_field_parsers[:domain][:derived] = parse_derived

"Parse action definition."
function parse_action(expr::Vector)
    args = Dict(expr[i].name => expr[i+1] for i in 1:2:length(expr))
    @assert (:action in keys(args)) ":action keyword is missing"
    name = args[:action]
    params, types = parse_typed_vars(get(args, :parameters, []))
    precondition = parse_formula(get(args, :precondition, []))
    effect = parse_formula(args[:effect])
    cont_prec = parse_formula(args[:continuous_precondition])
    cont_eff = parse_formula(args[:continuous_effect])
    dynamics = parse_formula(args[:dynamics])
    return GenericAction(name, params, types, precondition, effect, cont_prec,
                    cont_eff, dynamics)
end
body_field_parsers[:domain][:action] = parse_action