text:
import Lbfgsb

# 3<=m<=20 recommended
# factr =1.d+12 for low accuracy; factr =1.d+7 for moderate accuracy; factr =1.d+1 for extremely high accuracy.
# If factr =0, the test will stop the algorithm only if the objective function remains unchanged after one iteration.
function lbfgsb( z, alpha, c, lb = nothing; factr = 1e1, m=30, maxiter = 5000, iprint=-1, pgtol=1e-7)
  x       = vcat(z, alpha)
  x_      = [real(x);imag(x)]
  n       = length(x_)

  function objFunc(x)
    z, alpha = unpackComplex(x)
    z = residual(z, alpha, c)
    return z
  end

  function gradFunc!(x, g::Array)
    z, alpha = unpackComplex(x)
    y = vcat(Fgradient(z,alpha, c)...)
    g[:] = 2*[real(y);imag(y)]
    return
  end

  if (lb != nothing)
    lower = [ones(div(n,2))*-Inf;ones(div(n,4))*lb;ones(div(n,4))*-Inf] # Stores the lower bounds for every entry of x
    btype = [ones(div(n,2))*0;ones(div(n,4))*1;ones(div(n,4))*0] # Set bounds
    x_[x_ .< lower] = lower[x_ .< lower] # project into bounds
    if (!isempty(x_[x_ .< lower]))
      error("Could not project into boundaries.")
    end
    f, x, numCall, numIter, status = Lbfgsb.lbfgsb( objFunc, gradFunc!, x_, btype=btype, lb=lower, factr=factr, m=m, maxiter=maxiter, iprint=iprint, pgtol=pgtol )
  else
    f, x, numCall, numIter, status = Lbfgsb.lbfgsb( objFunc, gradFunc!, x_, m=m, factr=factr, maxiter=maxiter, iprint=iprint, pgtol=pgtol )
  end

  return unpackComplex(x)
end

function splitVec( x )
  return x[1:Int(length(x)/2)],x[Int(length(x)/2) + 1:end]
end

function unpackComplex(x)
  Re,Im = splitVec(x)
  z_re,a_re = splitVec(Re)
  z_im,a_im = splitVec(Im)

  return z_re + im .* z_im, a_re + im .* a_im
end


text:
# Implement the MutableArithmetics API for UnitAffExpr

const UnitAffExpr{U} = UnitExpression{AffExpr,U}

UnitAffExpr(aff::AffExpr, unit::U) where {U} = UnitExpression(aff, unit)

_MA.mutability(::Type{UnitAffExpr}) = _MA.IsMutable()

#  Minimal support of operate!! to allow @rewrite macro to work on linear
# constraints

_UnitAffOrVar = Union{UnitVariableRef,UnitAffExpr}
_AffOrVar = Union{AffExpr,VariableRef}
_AddSub = Union{typeof(_MA.add_mul),typeof(_MA.sub_mul)}

function Base.convert(::Type{UnitAffExpr}, uv::UnitVariableRef)
    return UnitAffExpr(convert(AffExpr, uv.variable), uv.unit)
end

function _update_expression(ua::UnitAffExpr, a::Number, x::UnitVariableRef)
    aval = Unitful.ustrip(
        Unitful.uconvert(ua.unit, a * Unitful.Quantity(1, x.unit)),
    )
    return UnitAffExpr(
        JuMP.add_to_expression!(ua.expr, aval, x.variable),
        ua.unit,
    )
end

function _update_expression(
    ua::UnitAffExpr,
    a::Unitful.Quantity,
    x::VariableRef,
)
    aval = Unitful.ustrip(Unitful.uconvert(ua.unit, a))
    return UnitAffExpr(JuMP.add_to_expression!(ua.expr, aval, x), ua.unit)
end

function _update_expression(ua::UnitAffExpr, a::Unitful.Quantity)
    aval = Unitful.ustrip(Unitful.uconvert(ua.unit, a))
    return UnitAffExpr(JuMP.add_to_expression!(ua.expr, aval), ua.unit)
end

function _update_expression(ua::UnitAffExpr, a::Number, x::UnitAffExpr)
    factor = Unitful.ustrip(
        Unitful.uconvert(ua.unit, a * Unitful.Quantity(1, x.unit)),
    )
    return UnitAffExpr(
        JuMP.add_to_expression!(ua.expr, factor, x.expr),
        ua.unit,
    )
end

function _update_expression(ua::UnitAffExpr, a::Number, b::Unitful.Quantity)
    aval = Unitful.ustrip(Unitful.uconvert(ua.unit, a * b))
    return UnitAffExpr(JuMP.add_to_expression!(ua.expr, aval), ua.unit)
end

function _update_expression(ua::UnitAffExpr, a::Unitful.Quantity, b::Number)
    aval = Unitful.ustrip(Unitful.uconvert(ua.unit, a * b))
    return UnitAffExpr(JuMP.add_to_expression!(ua.expr, aval), ua.unit)
end

function _update_expression(
    ua::UnitAffExpr,
    a::Unitful.Quantity,
    b::Unitful.Quantity,
)
    aval = Unitful.ustrip(Unitful.uconvert(ua.unit, a * b))
    return UnitAffExpr(JuMP.add_to_expression!(ua.expr, aval), ua.unit)
end

function _create_expression(
    t::_AddSub,
    z::typeof(_MA.Zero()),
    a::Number,
    x::UnitVariableRef,
)
    val = a * Unitful.Quantity(1, x.unit)
    return UnitAffExpr(
        _MA.operate!!(t, z, Unitful.ustrip(val), x.variable),
        Unitful.unit(val),
    )
end

function _create_expression(
    t::_AddSub,
    z::typeof(_MA.Zero()),
    a::Unitful.Quantity,
    x::VariableRef,
)
    return UnitAffExpr(
        _MA.operate!!(t, z, Unitful.ustrip(a), x),
        Unitful.unit(a),
    )
end

function _create_expression(
    t::_AddSub,
    z::typeof(_MA.Zero()),
    a::Unitful.Quantity,
)
    return UnitAffExpr(
        convert(AffExpr, _MA.operate!!(t, z, Unitful.ustrip(a))),
        Unitful.unit(a),
    )
end

# Two arguments
function _MA.operate!!(
    ::typeof(_MA.add_mul),
    uav::_UnitAffOrVar,
    x::UnitVariableRef,
)
    return _update_expression(convert(UnitAffExpr, uav), 1, x)
end

function _MA.operate!!(
    ::typeof(_MA.sub_mul),
    uav::_UnitAffOrVar,
    x::UnitVariableRef,
)
    return _update_expression(convert(UnitAffExpr, uav), -1, x)
end

function _MA.operate!!(
    ::typeof(_MA.add_mul),
    uav::_UnitAffOrVar,
    x::UnitAffExpr,
)
    return _update_expression(convert(UnitAffExpr, uav), 1, x)
end

function _MA.operate!!(
    ::typeof(_MA.sub_mul),
    uav::_UnitAffOrVar,
    x::UnitAffExpr,
)
    return _update_expression(convert(UnitAffExpr, uav), -1, x)
end

function _MA.operate!!(
    ::typeof(_MA.add_mul),
    uav::_UnitAffOrVar,
    a::Unitful.Quantity,
)
    return _update_expression(convert(UnitAffExpr, uav), a)
end

function _MA.operate!!(
    ::typeof(_MA.sub_mul),
    uav::_UnitAffOrVar,
    a::Unitful.Quantity,
)
    return _update_expression(convert(UnitAffExpr, uav), -a)
end

function _MA.operate!!(::_AddSub, ::_UnitAffOrVar, ::_AffOrVar)
    return error(
        "Can not combine variables with and without units in the samme expression",
    )
end

function _MA.operate!!(::_AddSub, ::_AffOrVar, ::_UnitAffOrVar)
    return error(
        "Can not combine variables with and without units in the samme expression",
    )
end

function _MA.operate!!(::_AddSub, ::_UnitAffOrVar, ::Real)
    return error("Can not combine unit expression with numbers without units")
end

function _MA.operate!!(::_AddSub, ::_AffOrVar, ::Unitful.Quantity)
    return error("Can not combine unit parameters with variables without units")
end

# Three arguments
function _MA.operate!!(
    ::typeof(_MA.add_mul),
    uav::_UnitAffOrVar,
    a::Number,
    x::_UnitAffOrVar,
)
    return _update_expression(
        convert(UnitAffExpr, uav),
        a,
        convert(UnitAffExpr, x),
    )
end

function _MA.operate!!(
    ::typeof(_MA.sub_mul),
    uav::_UnitAffOrVar,
    a::Number,
    x::_UnitAffOrVar,
)
    return _update_expression(
        convert(UnitAffExpr, uav),
        -a,
        convert(UnitAffExpr, x),
    )
end

function _MA.operate!!(
    ::typeof(_MA.add_mul),
    uav::_UnitAffOrVar,
    a::Unitful.Quantity,
    x::VariableRef,
)
    return _update_expression(convert(UnitAffExpr, uav), a, x)
end

function _MA.operate!!(
    ::typeof(_MA.sub_mul),
    uav::_UnitAffOrVar,
    a::Unitful.Quantity,
    x::VariableRef,
)
    return _update_expression(convert(UnitAffExpr, uav), -a, x)
end

function _MA.operate!!(
    ::typeof(_MA.add_mul),
    uav::_UnitAffOrVar,
    a::Number,
    b::Number,
)
    return _update_expression(convert(UnitAffExpr, uav), a, b)
end

function _MA.operate!!(
    ::typeof(_MA.sub_mul),
    uav::_UnitAffOrVar,
    a::Number,
    b::Number,
)
    return _update_expression(convert(UnitAffExpr, uav), -a, b)
end

function _MA.operate!!(
    ::typeof(_MA.add_mul),
    uav::_UnitAffOrVar,
    a::Number,
    b::Unitful.Units,
)
    return _update_expression(
        convert(UnitAffExpr, uav),
        a * Unitful.Quantity(1.0, b),
    )
end

function _MA.operate!!(
    ::typeof(_MA.sub_mul),
    uav::_UnitAffOrVar,
    a::Number,
    b::Unitful.Units,
)
    return _update_expression(
        convert(UnitAffExpr, uav),
        -a * Unitful.Quantity(1.0, b),
    )
end

function _MA.operate!!(::_AddSub, ::_UnitAffOrVar, ::Real, ::VariableRef)
    return error(
        "Can not combine variables with and without units in the samme expression",
    )
end

function _MA.operate!!(::_AddSub, ::_AffOrVar, ::Number, ::UnitVariableRef)
    return error(
        "Can not combine variables with and without units in the samme expression",
    )
end

function _MA.operate!!(
    ::_AddSub,
    ::_AffOrVar,
    ::Unitful.Quantity,
    ::VariableRef,
)
    return error(
        "Can not combine variables with and without units in the samme expression",
    )
end

# Multiple arguments
function _MA.operate!!(t::_AddSub, uav::_UnitAffOrVar, x, y, z, other_args...)
    args = (x, y, z, other_args...)
    n = length(args)
    varidx = findall(
        t ->
            (typeof(t) <: UnitAffExpr) ||
                (typeof(t) <: UnitVariableRef) ||
                (typeof(t) <: VariableRef),
        args,
    )
    var = args[varidx[1]]
    val = prod(args[i] for i in setdiff(1:n, varidx))
    return _MA.operate!!(t, uav, val, var)
end

# Zero handling
function _MA.operate!!(t::_AddSub, z::_MA.Zero, a::Unitful.Quantity)
    return _create_expression(t, z, a)
end

function _MA.operate!!(t::_AddSub, z::_MA.Zero, x::UnitVariableRef)
    return _create_expression(t, z, 1, x)
end

_MA.operate!!(::_AddSub, ::_MA.Zero, ua::UnitAffExpr) = ua

function _MA.operate!!(t::_AddSub, z::_MA.Zero, a::Number, x::UnitVariableRef)
    return _create_expression(t, z, a, x)
end

function _MA.operate!!(t::_AddSub, z::_MA.Zero, x::UnitVariableRef, a::Number)
    return _create_expression(t, z, a, x)
end

function _MA.operate!!(
    t::_AddSub,
    z::_MA.Zero,
    a::Unitful.Quantity,
    x::VariableRef,
)
    return _create_expression(t, z, a, x)
end

function _MA.operate!!(
    t::_AddSub,
    z::_MA.Zero,
    x::VariableRef,
    a::Unitful.Quantity,
)
    return _create_expression(t, z, a, x)
end

function _MA.operate!!(t::_AddSub, zero::_MA.Zero, x, y, z, other_args...)
    args = (x, y, z, other_args...)
    n = length(args)
    varidx = findall(
        t ->
            (typeof(t) <: UnitAffExpr) ||
                (typeof(t) <: UnitVariableRef) ||
                (typeof(t) <: VariableRef),
        args,
    )
    var = args[varidx[1]]
    val = prod(args[i] for i in setdiff(1:n, varidx))
    return _MA.operate!!(t, zero, val, var)
end

ext: jl


text:
@defcomp sealevelrise begin

    ThermSLR    = Variable(index=[time])    # Path of SLR from thermal expansion
    GSICSLR     = Variable(index=[time])    # Path of SLR from G&SIC
    GISSLR      = Variable(index=[time])    # Path of SLR from GIS
    AISSLR      = Variable(index=[time])    # Path of SLR from AIS
    TotSLR      = Variable(index=[time])    # Path of total SLR

    TATM        = Parameter(index=[time])   # Path of atmospheric temperature anomalies (from climate dynamics?)

    therm0      = Parameter()   # Initial SLR from thermal expansion (meters above 2000)
    gsic0       = Parameter()   # Initial SLR from G&SIC
    gis0        = Parameter()   # Initial SLR from GIS
    ais0        = Parameter()   # Initial SLR from AIS

    therm_asym  = Parameter()   # Asymptotic rise from thermal expansion
    gsic_asym   = Parameter()   # Asymptotic rise from G&SIC
    gis_asym    = Parameter()   # Asymptotic rise from GIS
    ais_asym    = Parameter()   # Asymptotic rise from AIS

    thermrate   = Parameter()   # Rate of thermal expansion
    gsicrate    = Parameter()   # Rate of G&SIC
    gisrate     = Parameter()   # Rate of GIS
    aisrate     = Parameter()   # Rate of AIS

    slrthreshold = Parameter()  # Temperature threshold for AIS

    function run_timestep(p, v, d, t)
        if is_first(t)
            v.ThermSLR[t] = p.therm0
            v.GSICSLR[t]  = p.gsic0
            v.GISSLR[t]   = p.gis0
            v.AISSLR[t]   = p.ais0
        else

            v.ThermSLR[t] = v.ThermSLR[t - 1] + p.thermrate * p.TATM[t]
            v.GSICSLR[t]  = v.GSICSLR[t - 1] + p.gsicrate * (p.gsic_asym - v.GSICSLR[t - 1]) * p.TATM[t]
            v.GISSLR[t]   = v.GISSLR[t - 1]  + p.gisrate  * (p.gis_asym  - v.GISSLR[t - 1])  * p.TATM[t]
            v.AISSLR[t]   = 0

            if p.TATM[t] > p.slrthreshold
                v.AISSLR[t] = v.AISSLR[t - 1] + p.aisrate * (p.ais_asym - v.AISSLR[t - 1]) * p.TATM[t]
            end
        end

        v.TotSLR[t] = v.ThermSLR[t] + v.GSICSLR[t] + v.GISSLR[t] + v.AISSLR[t]
    end
end


text:
using SFML

socket = SocketTCP()

connect(socket, "127.0.0.1", 53000, 10000)
println("Connected")

packet = Packet()
write(packet, "First String")
write(packet, "Second String")

send(socket, packet)


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module NullableArrays

using Compat
using Compat.view
using Reexport
@reexport using Base.Cartesian

export NullableArray,
       NullableVector,
       NullableMatrix,

       # Macros

       # Methods
       dropnull,
       anynull,
       allnull,
       nullify!,
       padnull!,
       padnull

include("typedefs.jl")
include("constructors.jl")
include("primitives.jl")
include("indexing.jl")
include("map.jl")
include("nullablevector.jl")
include("operators.jl")
include("broadcast.jl")
include("reduce.jl")
include("show.jl")
include("subarray.jl")

end


text:
using Documenter, TextModels

makedocs(
    modules = [TextModels],
    sitename = "TextAnalysis",
    format = Documenter.HTML(
    ),
    pages = [
        "Home" => "index.md",
        "Conditional Random Fields" => "crf.md",
        "ULMFiT" => "ULMFiT.md",
        "Named Entity Recognition" => "ner.md",
        "Tagging Schemes" => "tagging.md",
        "Sentiment Analyzer" => "sentiment.md",
        "API References" => "APIReference.md"
    ],
)


text:
using FastCap

# Generate the two wire shapes
np = 4
nl = 5
ew = 0.2
c1 = Cube("bar1",edgewidth=ew,size=(1,nl,1),panels=(np,(np-2)*nl+2,np))
c2 = Cube("bar2",edgewidth=ew,size=(nl,1,1),panels=((np-2)*nl+2,np,np))

# Generate the geometries
geoms = [
    Conductor("cond1",1,(1,0,0),c1),
    Conductor("cond2",1,(3,0,0),c1),
    Conductor("cond3",1,(0,1,2),c2),
    Conductor("cond4",1,(0,3,2),c2)
    ]

# View
#linedrawing(geoms,"examples\\buscrossing.pdf",x="5")

# Simulate!
res1 = fastcapwr(geoms,o=3,t=0.001)
res2 = fastcap2(geoms,o=3,t=0.001)
res3 = fastercap(geoms,a=0.005)

max_stars_repo_name: jebej/FastCap

max_issues_repo_name: jebej/FastCap

max_forks_repo_name: jebej/FastCap


text:
import GitHub

const BRANCH_OR_COMMIT = Union{GitHub.Branch, GitHub.Commit}

struct AlwaysAssertionError
    msg::String
end

function always_assert(condition::Bool, msg::String)
    if !condition
        throw(AlwaysAssertionError(msg))
    end
    return nothing
end

function _get_branch(repo;
                     api::GitHub.GitHubAPI,
                     auth::GitHub.Authorization,
                     branch)
    return GitHub.branch(api, repo, branch; auth = auth)
end

function _get_commit(repo;
                     api::GitHub.GitHubAPI,
                     auth::GitHub.Authorization,
                     commit)
    return GitHub.commit(api, repo, commit; auth = auth)
end

function _get_commit_or_branch(repo;
                               api::GitHub.GitHubAPI,
                               auth::GitHub.Authorization,
                               branch = nothing,
                               commit = nothing)
    if (branch isa Nothing) && (commit isa Nothing)
        msg = string(
            "You must provide either ",
            "the branch or the commit",
        )
        throw(ArgumentError(msg))
    elseif (!(branch isa Nothing)) && (!(commit isa Nothing))
        msg = string(
            "You cannot provide both ",
            "a branch and a commit. ",
            "You must provide exactly one.",
        )
        throw(ArgumentError(msg))
    elseif (!(branch isa Nothing))
        return _get_branch(repo; api=api, auth=auth, branch=branch)
    else
        always_assert(!(commit isa Nothing), "!(commit isa Nothing)")
        return _get_commit(repo; api = api, auth = auth, commit = commit)
    end
end

function get_statuses(repo;
                      api::GitHub.GitHubAPI,
                      auth::GitHub.Authorization,
                      branch,
                      commit)
    branch_or_commit = _get_commit_or_branch(
        repo;
        api = api,
        auth = auth,
        branch = branch,
        commit = commit,
    )
    return get_statuses(repo, branch_or_commit; api, auth)
end

function get_statuses(repo,
                      branch_or_commit::BRANCH_OR_COMMIT;
                      api::GitHub.GitHubAPI,
                      auth::GitHub.Authorization)
    combined_status = GitHub.status(
        api,
        repo,
        branch_or_commit;
        auth = auth,
    )
    combined_status_statuses = combined_status.statuses
    statuses, _ = GitHub.statuses(
        api,
        repo,
        branch_or_commit;
        auth = auth,
    )
    status_contexts = vcat(
        [x.context for x in combined_status_statuses],
        [x.context for x in statuses],
    )
    unique!(status_contexts)
    sort!(status_contexts)
    return status_contexts
end

function _convert_to_commit(repo,
                            branch::GitHub.Branch;
                            api::GitHub.GitHubAPI,
                            auth::GitHub.Authorization)
    commit = GitHub.commit(api, repo, branch; auth = auth)
    return commit
end

function _convert_to_commit(repo,
                            commit::GitHub.Commit;
                            api::GitHub.GitHubAPI,
                            auth::GitHub.Authorization)
    return commit
end

function _get_commit_sha_from_commit(commit::GitHub.Commit)
    return commit.sha
end

function get_check_runs(repo;
                        api::GitHub.GitHubAPI,
                        auth::GitHub.Authorization,
                        branch,
                        commit)
    return nothing
end

function get_check_runs(repo,
                        branch_or_commit::BRANCH_OR_COMMIT;
                        api::GitHub.GitHubAPI,
                        auth::GitHub.Authorization)
    commit = _convert_to_commit(
        repo,
        branch_or_commit;
        api = api,
        auth = auth,
    )
    commit_sha = _get_commit_sha_from_commit(commit)
    endpoint = "/repos/$(repo.full_name)/commits/$(commit_sha)/check-runs"
    check_runs = GitHub.gh_get_json(
        api,
        endpoint;
        auth = auth,
    )
    check_run_names = [x["name"] for x in check_runs["check_runs"]]
    unique!(check_run_names)
    sort!(check_run_names)
    return check_run_names
end

function get_statuses_and_check_runs(repo;
                                     api::GitHub.GitHubAPI,
                                     auth::GitHub.Authorization,
                                     branch,
                                     commit)
    status_contexts_and_check_run_names = get_statuses_and_check_runs(
        repo,
        branch_or_commit;
        api = api,
        auth = auth)
    unique!(status_contexts_and_check_run_names)
    sort!(status_contexts_and_check_run_names)
    return status_contexts_and_check_run_names
end

function get_statuses_and_check_runs(repo,
                                     branch_or_commit::BRANCH_OR_COMMIT;
                                     api::GitHub.GitHubAPI,
                                     auth::GitHub.Authorization)

    status_contexts = get_statuses(
        repo,
        branch_or_commit;
        api = api,
        auth = auth,
    )
    check_run_names = get_check_runs(
        repo,
        branch_or_commit;
        api = api,
        auth = auth,
    )
    status_contexts_and_check_run_names = vcat(
        status_contexts,
        check_run_names,
    )
    unique!(status_contexts_and_check_run_names)
    sort!(status_contexts_and_check_run_names)
    return status_contexts_and_check_run_names
end

function _show(v::AbstractVector)
    for (i, x) ∈ enumerate(v)
        println(stdout, "$(i). $(x)")
    end
end

function main()
    api = GitHub.DEFAULT_API
    auth = GitHub.AnonymousAuth()
    repo = GitHub.repo(api, "JuliaGPU/KernelAbstractions.jl"; auth = auth);
    # branch_name = "master"
    # branch_name = "staging"
    branch_name = "trying"
    branch = _get_commit_or_branch(repo; api = api, auth = auth, branch = branch_name);
    result = get_statuses_and_check_runs(repo, branch; api = api, auth = auth)
    @info "" branch_name
    show(result)
    _show(result)
    return result
end

main()


text:
using Viznet: canvas
import Viznet
using Compose: CurvePrimitive, Form
using YaoBlocks
using BitBasis

export CircuitStyles, CircuitGrid, circuit_canvas, vizcircuit

module CircuitStyles
    using Compose
    import Viznet
    const r = Ref(0.2)
	const lw = Ref(1pt)
	const textsize = Ref(16pt)
	const paramtextsize = Ref(10pt)
	const fontfamily = Ref("Helvetica Neue")
    const linecolor = Ref("#000000")
    const gate_bgcolor = Ref("#FFFFFF")
    const textcolor = Ref("#000000")
    const scale = Ref(1.0)
    getlw() = lw[] * scale[]
    gettextsize() = textsize[] * scale[]
    getparamtextsize() = paramtextsize[] * scale[]
    G() = compose(context(), rectangle(-r[], -r[], 2*r[], 2*r[]), fill(gate_bgcolor[]), stroke(linecolor[]), linewidth(getlw()))
    C() = compose(context(), circle(0.0, 0.0, r[]/3), fill(linecolor[]), linewidth(0))
    NC() = compose(context(), circle(0.0, 0.0, r[]/3), fill(gate_bgcolor[]), stroke(linecolor[]), linewidth(getlw()))
    X() = compose(context(), xgon(0.0, 0.0, r[], 4), fill(linecolor[]), linewidth(0))
	NOT() = compose(context(),
                    (context(), circle(0.0, 0.0, r[]), stroke(linecolor[]), linewidth(getlw()), fill("transparent")),
                    (context(), polygon([(-r[], 0.0), (r[], 0.0)]), stroke(linecolor[]), linewidth(getlw())),
                    (context(), polygon([(0.0, -r[]), (0.0, r[])]), stroke(linecolor[]), linewidth(getlw()))
               )
    WG() = compose(context(), rectangle(-1.5*r[], -r[], 3*r[], 2*r[]), fill(gate_bgcolor[]), stroke(linecolor[]), linewidth(getlw()))
    MULTIGATE(h) = compose(context(), rectangle(-1.5*r[], -(h/2+r[]), 3*r[], (h+2*r[])), fill(gate_bgcolor[]), stroke(linecolor[]), linewidth(getlw()))
    LINE() = compose(context(), line(), stroke(linecolor[]), linewidth(getlw()))
    TEXT() = compose(context(), text(0.0, 0.0, "", hcenter, vcenter), fontsize(gettextsize()), fill(textcolor[]), font(fontfamily[]))
    PARAMTEXT() = compose(context(), text(0.0, 0.0, "", hcenter, vcenter), fontsize(getparamtextsize()), fill(textcolor[]), font(fontfamily[]))
    MEASURE() = compose(context(),
                        rectangle(-r[], -r[], 2*r[], 2*r[]), fill(gate_bgcolor[]), stroke(linecolor[]), linewidth(getlw()),
                        compose(context(), curve((-0.8*r[], 0.5*r[]), (-0.8*r[], -0.6*r[]), (0.8*r[], -0.6*r[]), (0.8*r[], 0.5*r[])), stroke(linecolor[]), linewidth(getlw())),
                        compose(context(), line([(0.0, 0.5*r[]), (0.7*r[], -0.4*r[])]), stroke(linecolor[]), linewidth(getlw())),
        begin
            ns = Viznet.nodestyle(:triangle, fill(linecolor[]); r=0.1*r[], θ=atan(0.7, 0.9))
            Viznet.inner_most_containers(ns) do c
                Viznet.update_locs!(c.form_children, [(0.7*r[], -0.4*r[])])
            end
            ns
        end
        )
	function setlw(_lw)
		lw[] = _lw
	end
	function setr(_r)
		r[] = _r
	end
	function settextsize(_textsize)
		textsize[] = _textsize
	end
	function setparamtextsize(_paramtextsize)
		paramtextsize[] = _paramtextsize
	end
	function setfontfamily(_fontfamily)
		fontfamily[] = _fontfamily
	end
end

struct CircuitGrid
    frontier::Vector{Int}
	w_depth::Float64
	w_line::Float64
end

nline(c::CircuitGrid) = length(c.frontier)
depth(c::CircuitGrid) = frontier(c, 1, nline(c))
Base.getindex(c::CircuitGrid, i, j) = (c.w_depth*i, c.w_line*j)
Base.typed_vcat(c::CircuitGrid, ij1, ij2) = (c[ij1...], c[ij2...])

function CircuitGrid(nline::Int; w_depth=1.0, w_line=1.0)
    CircuitGrid(zeros(Int, nline), w_depth, w_line)
end

function frontier(c::CircuitGrid, args...)
    maximum(i->c.frontier[i], min(args..., nline(c)):max(args..., 1))
end

function _draw!(c::CircuitGrid, loc_brush_texts)
    isempty(loc_brush_texts) && return
    # a loc can be a integer, or a range
    locs = Iterators.flatten(getindex.(loc_brush_texts, 1)) |> collect
    i = frontier(c, locs...) + 1
	local jpre
	loc_brush_texts = sort(loc_brush_texts, by=x->first(x[1]))
    for (k, (j, b, txt)) in enumerate(loc_brush_texts)
        length(j) == 0 && continue
        jmid = (minimum(j)+maximum(j))/2
		b >> c[i, jmid]
		if length(txt) >= 3
			CircuitStyles.PARAMTEXT() >> (c[i, jmid], txt)
		elseif length(txt) >= 1
			CircuitStyles.TEXT() >> (c[i, jmid], txt)
		end
		if k!=1
			CircuitStyles.LINE() >> c[(i, jmid); (i, jpre)]
		end
		jpre = jmid
    end

	jmin, jmax = min(locs..., nline(c)), max(locs..., 1)
	for j = jmin:jmax
		CircuitStyles.LINE() >> c[(i, j); (c.frontier[j], j)]
		c.frontier[j] = i
	end
end

function _draw_continuous_multiqubit!(c::CircuitGrid, loc_text)
    (start, stop), txt = loc_text
    stop-start<0 && return
    b = CircuitStyles.MULTIGATE((stop-start) * c.w_line)
    i = frontier(c, start:stop...) + 1
    j = (stop+start)/2

    b >> c[i, j]
    if length(txt) >= 3
        CircuitStyles.PARAMTEXT() >> (c[i, j], txt)
    elseif length(txt) >= 1
        CircuitStyles.TEXT() >> (c[i, j], txt)
    end
	for j = start:stop
		CircuitStyles.LINE() >> c[(i, j); (c.frontier[j], j)]
		c.frontier[j] = i
	end
end

function initialize!(c::CircuitGrid; starting_texts, starting_offset)
	starting_texts !== nothing && for j=1:nline(c)
        CircuitStyles.TEXT() >> (c[starting_offset, j], string(starting_texts[j]))
	end
end

function finalize!(c::CircuitGrid; show_ending_bar, ending_offset, ending_texts)
    i = frontier(c, 1, nline(c)) + 1
	for j=1:nline(c)
		show_ending_bar && CircuitStyles.LINE() >> c[(i, j-0.2); (i, j+0.2)]
		CircuitStyles.LINE() >> c[(i, j); (c.frontier[j], j)]
        ending_texts !== nothing && CircuitStyles.TEXT() >> (c[i+ending_offset, j], string(ending_texts[j]))
	end
	c.frontier .= i
end

# elementary
function draw!(c::CircuitGrid, b::AbstractBlock, address, controls)
    error("block type $(typeof(b)) does not support visualization.")
end

function draw!(c::CircuitGrid, p::PrimitiveBlock{M}, address, controls) where M
	bts = length(controls)>=1 ? get_cbrush_texts(p) : get_brush_texts(p)
	_draw!(c, [controls..., [(address[i], bts[i]...) for i=occupied_locs(p)]...])
end

function draw!(c::CircuitGrid, p::Scale, address, controls)
	fp = YaoBlocks.factor(p)
	if !(abs(fp) ≈ 1)
		error("can not visualize non-phase factor.")
	end
	draw!(c, YaoBlocks.phase(angle(fp)), [first(address)], controls)
	draw!(c, p.content, address, controls)
end

# composite
function draw!(c::CircuitGrid, p::ChainBlock{N}, address, controls) where N
	draw!.(Ref(c), subblocks(p), Ref(address), Ref(controls))
end

function draw!(c::CircuitGrid, p::PutBlock{N,M}, address, controls) where {N,M}
	locs = [address[i] for i in p.locs]
	draw!(c, p.content, locs, controls)
end

function draw!(c::CircuitGrid, m::YaoBlocks.Measure{N}, address, controls) where {N,M}
    if m.postprocess isa RemoveMeasured
        error("can not visualize post-processing: `RemoveMeasured`.")
    end
    if !(m.operator isa ComputationalBasis)
        error("can not visualize measure blocks for operators")
    end
    locs = m.locations isa AllLocs ? collect(1:N) : [address[i] for i in m.locations]
    for (i, loc) in enumerate(locs)
        _draw!(c, [(loc, CircuitStyles.MEASURE(), "")])
        if m.postprocess isa ResetTo
            val = readbit(m.postprocess.x, i)
            _draw!(c, [(loc, CircuitStyles.G(), val == 1 ? "P₁" : "P₀")])
        end
    end
end

function draw!(c::CircuitGrid, cb::ControlBlock{N,GT,C}, address, controls) where {N,GT,C}
    ctrl_locs = [address[i] for i in cb.ctrl_locs]
    locs = [address[i] for i in cb.locs]
	mycontrols = [(loc, (bit == 1 ? CircuitStyles.C() : CircuitStyles.NC()), "") for (loc, bit)=zip(ctrl_locs, cb.ctrl_config)]
	draw!(c, cb.content, locs, [controls..., mycontrols...])
end

function draw!(c::CircuitGrid, cb::LabelBlock{GT,N}, address, controls) where {N,GT}
    length(address) == 0 && return
    is_continuous_chunk(address) || error("address not continuous in a block marked as continous.")
	_draw!(c, [controls..., (minimum(address):maximum(address), CircuitStyles.MULTIGATE((length(address)-1)*c.w_line), cb.name)])
end

for (GATE, SYM) in [(:XGate, :Rx), (:YGate, :Ry), (:ZGate, :Rz)]
	@eval get_brush_texts(b::RotationGate{1,T,<:$GATE}) where T = [(CircuitStyles.WG(), "$($(SYM))($(pretty_angle(b.theta)))")]
end

pretty_angle(theta) = string(theta)
function pretty_angle(theta::AbstractFloat)
	c = ZXCalculus.continued_fraction(theta/π, 10)
	if c.den < 100
		res = if c.num == 1
			"π"
		elseif c.num==0
			"0"
		elseif c.num==-1
			"-π"
		else
			"$(c.num)π"
		end
		if c.den != 1
			res *= "/$(c.den)"
		end
		res
	else
		"$(round(theta; digits=2))"
	end
end

get_brush_texts(b::ConstGate.CNOTGate) = [(CircuitStyles.C(), ""), (CircuitStyles.X(), "")]
get_brush_texts(b::ConstGate.CZGate) = [(CircuitStyles.C(), ""), (CircuitStyles.C(), "")]
get_brush_texts(b::ConstGate.ToffoliGate) = [(CircuitStyles.C(), ""), (CircuitStyles.C(), ""), (CircuitStyles.X(), "")]
get_brush_texts(b::ConstGate.SdagGate) = [(CircuitStyles.G(), "S'")]
get_brush_texts(b::ConstGate.TdagGate) = [(CircuitStyles.G(), "T'")]
get_brush_texts(b::ConstGate.PuGate) = [(CircuitStyles.G(), "P+")]
get_brush_texts(b::ConstGate.PdGate) = [(CircuitStyles.G(), "P-")]
get_brush_texts(b::ConstGate.P0Gate) = [(CircuitStyles.G(), "P₀")]
get_brush_texts(b::ConstGate.P1Gate) = [(CircuitStyles.G(), "P₁")]
get_brush_texts(b::ConstGate.I2Gate) = []
get_brush_texts(b::SWAPGate) = [(CircuitStyles.X(), ""), (CircuitStyles.X(), "")]
get_brush_texts(b::PrimitiveBlock{M}) where M = fill((CircuitStyles.G(), ""), M)
get_brush_texts(b::PrimitiveBlock{1}) = [(CircuitStyles.G(), "")]
get_brush_texts(b::ShiftGate) = [(CircuitStyles.WG(), "ϕ($(pretty_angle(b.theta)))")]
get_brush_texts(b::PhaseGate) = [(CircuitStyles.WG(), "^$(pretty_angle(b.theta))")]
function get_brush_texts(b::T) where T<:ConstantGate{1}
    namestr = string(T.name.name)
    if endswith(namestr, "Gate")
        namestr = namestr[1:end-4]
    end
    [(CircuitStyles.G(), namestr)]
end

get_cbrush_texts(b::PrimitiveBlock) = get_brush_texts(b)
get_cbrush_texts(b::XGate) = [(CircuitStyles.NOT(), "")]
get_cbrush_texts(b::ZGate) = [(CircuitStyles.C(), "")]

# front end
plot(blk::AbstractBlock; kwargs...) = vizcircuit(blk; kwargs...)
function vizcircuit(blk::AbstractBlock; w_depth=0.85, w_line=0.75, scale=1.0, show_ending_bar=false, starting_texts=nothing, starting_offset=-0.3, ending_texts=nothing, ending_offset=0.3)
    CircuitStyles.scale[] = scale
    img = circuit_canvas(nqubits(blk); w_depth=w_depth, w_line=w_line, show_ending_bar=show_ending_bar, starting_texts=starting_texts, starting_offset=starting_offset, ending_texts=ending_texts, ending_offset=ending_offset) do c
		basicstyle(blk) >> c
	end
    CircuitStyles.scale[] = 1.0
    return img |> rescale(scale)
end

function circuit_canvas(f, nline::Int; w_depth=0.85, w_line=0.75, show_ending_bar=false, starting_texts=nothing, starting_offset=-0.3, ending_texts=nothing, ending_offset=0.3)
	c = CircuitGrid(nline; w_depth=w_depth, w_line=w_line)
	g = canvas() do
        initialize!(c; starting_texts=starting_texts, starting_offset=starting_offset)
        f(c)
        finalize!(c; show_ending_bar=show_ending_bar, ending_texts=ending_texts, ending_offset=ending_offset)
	end
	a, b = (depth(c)+1)*w_depth, nline*w_line
	Compose.set_default_graphic_size(a*2.5*cm, b*2.5*cm)
	compose(context(0.5/a, -0.5/b, 1/a, 1/b), g)
end

Base.:>>(blk::AbstractBlock{N}, c::CircuitGrid) where N = draw!(c, blk, collect(1:N), [])
Base.:>>(blk::Function, c::CircuitGrid) = blk(nline(c)) >> c

function rescale(factor)
	a, b = Compose.default_graphic_width, Compose.default_graphic_height
	Compose.set_default_graphic_size(a*factor, b*factor)
	graph -> compose(context(), graph)
end

vizcircuit(; kwargs...) = c->vizcircuit(c; kwargs...)

function basicstyle(blk::AbstractBlock)
	YaoBlocks.Optimise.simplify(blk, rules=[YaoBlocks.Optimise.to_basictypes])
end


text:
__precompile__()

module MarketTechnicals

using Reexport
using StatsBase

@reexport using TimeSeries

export sma, ema, kama, env,
       bollingerbands, truerange, atr, keltnerbands, chaikinvolatility, donchianchannels,
       obv, vwap, adl,
       doji,
       rsi, macd, cci, roc, adx, stochasticoscillator, chaikinoscillator, aroon,
       floorpivots, woodiespivots,
       typical

include("candlesticks.jl")
include("levels.jl")
include("movingaverages.jl")
include("momentum.jl")
include("utilities.jl")
include("volatility.jl")
include("volume.jl")

# for user customization
# FIXME: using @__DIR__ while we upgrade to julia 0.6+
RC_FILE = joinpath(dirname(@__FILE__), ".rc.jl")
if !ispath(RC_FILE)
    touch(RC_FILE)
end

include(RC_FILE)

end


text:
using DataFrames
using CSV
using Plots
#%%
data_k1 = open("qho_res.dat")
df_k1=CSV.read(data_k1, header=false, delim=" ")
r_k1=df_k1[1]
phi_k1=df_k1[2]

#%%
plot(df_k1[1], df_k1[2], label="res", dpi=600)
plot!(grid=true, legend=:topleft)
xlabel!("E")
ylabel!("residual")
savefig("res.png")
# xaxis!(:log10)
# yaxis!(:log10)
#%%
2+2

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
if !@isdefined TEST
    using Revise
end
using Soss
using Mitosis
using Random, Test, LinearAlgebra, Statistics

# Define some vectors and matrices
ξ0 = [1., 0.]
x = [1.2, 0.1]
P0 = Matrix(1.0I, 2, 2)

Φ = [0.8 0.5; -0.1 0.8]
β = [0.1, 0.2]
Q = [0.2 0.0; 0.0 1.0]

H = [1.0 0.0]
R = Matrix(1.0I, 1, 1)

# and a nonlinear vector function
f(x) = [0.8 0.5; -0.1 0.8]*[atan(x[1]), atan(x[2])] + [0.1, 0.2]

# Define some transition kernels.

# We define the equivalent of the Soss model
m = @model ξ0 begin
    a ~ MvNormal(ξ0, P0) # priortransition
    b ~ MvNormal(H*a, R) # partialobservation
    c ~ MvNormal(f(a), Q) # nonlineartransition
    d ~ MvNormal(H*c, R) # partial observations
    e ~ MvNormal(c, P0) # full observation
    # return b, d, e
end

# We use AffineMap, LinearMap and ConstantMap
# For example
@test AffineMap(Φ, β)(x) == Φ*x + β

# Prior
# a ~ N(ξ0, P0)
k_0a = Mitosis.kernel(Gaussian; μ=LinearMap(I(2)), Σ=ConstantMap(P0))

# Partial observation
# b ~ N(H*a, R)
k_ab = Mitosis.kernel(Gaussian; μ=LinearMap(H), Σ=ConstantMap(R))

# Nonlinear transition and linear approximation
# c ~ N(f(a), Q)
k_ac = Mitosis.kernel(Gaussian; μ=f, Σ=ConstantMap(Q))
k_ac_approx = Mitosis.kernel(Gaussian; μ=AffineMap(Φ, β), Σ=ConstantMap(Q))

# d ~ N(H*c, R)
k_cd = Mitosis.kernel(Gaussian; μ=LinearMap(H), Σ=ConstantMap(R))

# e ~ N(c, P0)
k_ce = Mitosis.kernel(Gaussian; μ=LinearMap(I(2)), Σ=ConstantMap(P0))

# And a deterministic copy kernel
cp2 = Mitosis.Copy{2}() # copy kernel
@test cp2(1.1) == (1.1, 1.1)

# Forward sample a Bayes net # rand(m(ξ0=ξ0))
a = rand(k_0a(ξ0))
b = rand(k_ab(a))
c = rand(k_ac(a)) # forward model with nonlinear transition
d = rand(k_cd(c))
e = rand(k_ce(c))

# We actually want to write down the
# forward sample with explicit copies
# so every state is only used as input of a single kernel
a = rand(k_0a(ξ0))
a1,a2 = cp2(a)
b = rand(k_ab(a1))
c = rand(k_ac(a2)) # forward model with nonlinear transition
c1,c2 = cp2(c)
d = rand(k_cd(c1))
e = rand(k_ce(c2))

# observations are b d e (all leaves)
mc2, pc2 = backward(BFFG(), k_ce, e; unfused=true)
mc1, pc1 = backward(BFFG(), k_cd, d; unfused=true)
mc, pc = backward(BFFG(), cp2, pc1, pc2; unfused=true)
ma2, pa2 = backward(BFFG(), k_ac_approx, pc; unfused=true)
ma1, pa1 = backward(BFFG(), k_ab, b; unfused=true)
ma, pa = backward(BFFG(), cp2, pa1, pa2; unfused=true)
m0, p0 = backward(BFFG(), k_0a, pa) # not a child of a copy

samples = []
lls = []
for i in 1:1000
    a = rand(forward(BFFG(), k_0a, m0, weighted(ξ0)))
    a1, a2 = rand(forward(BFFG(), cp2, ma, a))
    b = rand(forward(BFFG(), k_ab, ma1, a1))
    c = rand(forward(BFFG(), (k_ac, k_ac_approx), ma2, a2)) # forward model with nonlinear transition
    c1, c2 = rand(forward(BFFG(), cp2, mc, c))
    d = rand(forward(BFFG(), k_cd, mc1, c1))
    e = rand(forward(BFFG(), k_ce, mc2, c2))

    # Sample and weight (sum over the weights of the leafs.)
    push!(samples, (a=a.x,b=b.x,c=c.x,d=d.x,e=e.x))
    push!(lls, b.ll + d.ll + e.ll)

end


text:
module STBOutput

using DataFrames: DataFrame, DataFrameRow, Not, select!

using PovertyAndInequalityMeasures

using CSV

using StatsBase
using ScottishTaxBenefitModel
using .Definitions
using .Utils
using .GeneralTaxComponents:
    WEEKS_PER_YEAR

using .Results:
    BenefitUnitResult,
    HouseholdResult,
    IndividualResult

using .STBIncomes

using .ModelHousehold

using .FRSHouseholdGetter:
    get_slot_for_household,
    get_slot_for_person

using .RunSettings

using .SimplePovertyCounts:
    GroupPoverty,
    calc_child_poverty

export
    add_to_frames!,
    dump_frames,
    initialise_frames,
    summarise_frames,
    make_poverty_line,
    make_gain_lose,
    hdiff,
    idiff,
    irdiff

# count of the aggregates added to the income_frame - total benefits and so on
const EXTRA_INC_COLS = 10

    #=
    struct FrameStarts
        hh :: Integer
        bu :: Integer
        pers :: Integer
        income :: Integer
    end
    =#

    function make_household_results_frame( n :: Int ) :: DataFrame
        make_household_results_frame( Float64, n )
    end

    function make_household_results_frame( RT :: DataType, n :: Int ) :: DataFrame
        DataFrame(
            hid       = zeros( BigInt, n ),
            sequence  = zeros( Int, n ),
            data_year  = zeros( Int, n ),

            weight    = zeros(RT,n),
            weighted_people = zeros(RT,n),
            hh_type   = zeros( Int, n ),
            num_people = zeros( Int, n ),
            tenure    = fill( Missing_Tenure_Type, n ),
            region    = fill( Missing_Standard_Region, n ),
            gross_decile = zeros( Int, n ),
            bhc_net_income = zeros(RT,n),
            ahc_net_income = zeros(RT,n),
            # eq_scale = zeros(RT,n),
            eq_bhc_net_income = zeros(RT,n),
            eq_ahc_net_income = zeros(RT,n), # etc.
            income_taxes = zeros(RT,n),
            means_tested_benefits = zeros(RT,n),
            other_benefits = zeros(RT,n),
            scottish_income_tax = zeros(RT,n),
            num_children = zeros(RT,n)
        )
    end

    function hdiff( df1 :: DataFrame, df2 :: DataFrame )
        p1 = index_of_field( df1, "bhc_net_income")
        p2 = index_of_field( df1, "num_children")-1
        return df_diff( df1, df2, p1, p2 )
    end

    function make_bu_results_frame( n :: Int ) :: DataFrame
        return make_bu_results_frame( Float64, n )
    end

    function make_bu_results_frame( RT :: DataType, n :: Int ) :: DataFrame
        DataFrame(
            hid       = zeros(BigInt,n),
            buno      = zeros( Int, n ),
            data_year = zeros( Int, n ),
            weight    = zeros(RT,n),
            bu_type   = zeros( Int, n ),
            tenure    = zeros( Int, n ),
            region    = zeros( Int, n ),
            gross_decile = zeros( Int, n ),
            net_income = zeros(RT,n),

            income_taxes = zeros(RT,n),
            means_tested_benefits = zeros(RT,n),
            other_benefits = zeros(RT,n)
        )
    end

    function make_individual_results_frame( n :: Int ) :: DataFrame
        make_individual_results_frame( Float64, n )
    end

    function make_incomes_frame( RT :: DataType, n :: Int; id = 1 ) :: DataFrame
        frame :: DataFrame = create_incomes_dataframe( RT, n )
        # extra calculated fields
        frame.employers_ni = zeros( n )
        frame.scottish_income_tax = zeros( n ) # i.e. excluding savings & dividends
        frame.total_benefits = zeros( n )
        frame.legacy_mtbs  = zeros( n )
        frame.means_tested_bens = zeros( n )
        frame.non_means_tested_bens = zeros( n )
        frame.sickness_illness = zeros( n )
        frame.scottish_benefits = zeros( n )
        frame.pension_relief_at_source = zeros( n )
        frame.net_cost = zeros( n )
        # add some crosstab fields ...
        frame.id = fill( id, n )
        frame.data_year = zeros( Int, n )
        frame.sex = fill(Missing_Sex,n)
        frame.ethnic_group = fill(Missing_Ethnic_Group,n)
        frame.is_child = fill( false, n )
        frame.age_band  = zeros(Int,n)
        frame.employment_status = fill(Missing_ILO_Employment,n)
        frame.tenure    = fill( Missing_Tenure_Type, n )
        frame.region    = fill( Missing_Standard_Region, n )
        frame.gross_decile = zeros( Int, n )
        frame.council = fill( Symbol( "No_Council"), n)

        # ... and so on
        return frame
    end

    """
    Change in the numeric fields of incomes
    """
    function idiff( d1 :: DataFrame, d2 :: DataFrame )
        # the numeric fields are between "income_tax" and "id"
        p1 = index_of_field( d1, "income_tax" )
        p2 = index_of_field( d1, "id")-1
        return df_diff( d1, d2, p1, p2 )
    end

    function make_individual_results_frame( RT :: DataType, n :: Int ) :: DataFrame
       DataFrame(
         hid = zeros(BigInt,n),
         pid = zeros(BigInt,n),
         data_year = zeros(Int,n),
         weight = zeros(RT,n),
         sex = fill(Missing_Sex,n),
         ethnic_group = fill(Missing_Ethnic_Group,n),
         is_child = fill( false, n ),
         age_band  = zeros(Int,n),
         employment_status = fill(Missing_ILO_Employment,n),
         # ... and so on

         income_taxes = zeros(RT,n),
         means_tested_benefits = zeros(RT,n),
         other_benefits = zeros(RT,n),

         income_tax = zeros(RT,n),
         it_non_savings = zeros(RT,n),
         it_savings = zeros(RT,n),
         it_dividends = zeros(RT,n),
         it_pension_relief_at_source = zeros(RT,n),
         ni_above_lower_earnings_limit = fill( false, n ),
         ni_total_ni = zeros(RT,n),
         ni_class_1_primary = zeros(RT,n),
         ni_class_1_secondary = zeros(RT,n),
         ni_class_2  = zeros(RT,n),
         ni_class_3  = zeros(RT,n),
         ni_class_4  = zeros(RT,n),
         assumed_gross_wage = Vector{Union{Real,Missing}}(missing, n),
         metr = Vector{Union{Real,Missing}}(missing, n),
         tax_credit = zeros(RT,n),
         replacement_rate = Vector{Union{Real,Missing}}(missing, n))
    end

    function irdiff( d1 :: DataFrame, d2 :: DataFrame )
        # the numeric fields are between "income_tax" and "id"
        p1 = index_of_field( d1, "income_taxes" )
        p2 = size(d1)[2]
        return df_diff( d1, d2, p1, p2 )
    end

    function initialise_frames( T::DataType, settings :: Settings, num_systems :: Integer  ) :: NamedTuple
        indiv = []
        bu = []
        hh = []
        income = []
        for s in 1:num_systems
            push!(indiv, make_individual_results_frame( T, settings.num_people ))
            push!(bu, make_bu_results_frame( T, settings.num_people )) # overstates but we don't actually know this at the start
            push!(hh, make_household_results_frame( T, settings.num_households ))
            push!(income, make_incomes_frame( T, settings.num_people )) # overstates but we don't actually know this at the start
        end
        (hh=hh, bu=bu, indiv=indiv, income=income)
    end

    function fill_hh_frame_row!( hr :: DataFrameRow, hh :: Household, hres :: HouseholdResult )
        nps =  num_people(hh)
        hr.hid = hh.hid
        hr.sequence = hh.sequence
        hr.weight = hh.weight
        hr.data_year = hh.data_year
        hr.weighted_people = hh.weight*nps
        hr.num_people = nps
        hr.hh_type = -1
        hr.tenure = hh.tenure
        hr.region = hh.region
        hr.gross_decile = -1
        hr.income_taxes = isum(hres.income, INCOME_TAXES )
        hr.means_tested_benefits = isum( hres.income, MEANS_TESTED_BENS )
        hr.other_benefits = isum( hres.income, NON_MEANS_TESTED_BENS )
        hr.bhc_net_income = hres.bhc_net_income
        hr.ahc_net_income = hres.ahc_net_income
        # hr.eq_scale = hres.eq_scale
        hr.eq_bhc_net_income = hres.eq_bhc_net_income
        hr.eq_ahc_net_income = hres.eq_ahc_net_income
        hr.num_children = num_children( hh )
    end

    function fill_bu_frame_row!(
        br :: DataFrameRow,
        hh :: Household,
        bres :: BenefitUnitResult )

        # ...

    end

    function summarise_inc_frame( incd :: DataFrame ) :: DataFrame
        nrows = 80
        out = make_incomes_frame(Float64, nrows)
        out.label = fill("",nrows)

        # labels
        out[1,:label]="Grant Total £p.a"
        out[2,:label]="Counts"
        row = 3
        for em in instances( ILO_Employment )
            out[row,:label]="$em - £p.a"
            row += 1
            out[row,:label]="$em - counts"
            row += 1
        end
        for em in instances( Tenure_Type )
            out[row,:label]="$em - £p.a"
            row += 1
            out[row,:label]="$em - counts"
            row += 1
        end
        for a in 1:17
            em = age_str(a)
            out[row,:label]="$em - £p.a"
            row += 1
            out[row,:label]="$em - counts"
            row += 1
        end

        col = 3
        # FIXME is there something subtly wrong with the weighting here?
        for i in 1:(INC_ARRAY_SIZE+EXTRA_INC_COLS)
            col += 1
            # println( "on column $col")
            out[1,col] = sum( WEEKS_PER_YEAR .* incd[:,col] .* incd[:,:weight] ) # £mn
            out[2,col] = sum((incd[:,col] .> 0) .* incd[:,:weight]) # counts
            row = 3
            for em in instances( ILO_Employment )
                selected = incd.employment_status.==em
                out[row,col] = sum( WEEKS_PER_YEAR .* incd[selected,col] .* incd[selected,:weight] ) # £mn
                row += 1
                out[row,col] = sum((incd[selected,col] .> 0) .* incd[selected,:weight]) # Counts
                row += 1
            end
            for em in instances( Tenure_Type )
                selected = incd.tenure .== em
                out[row,col] =sum(  WEEKS_PER_YEAR .* incd[selected,col] .* incd[selected,:weight] ) # £mn
                row += 1
                out[row,col] = sum((incd[selected,col] .> 0) .* incd[selected,:weight]) # Counts
                row += 1
            end
            for a in 1:17
                selected = incd.age_band .== a
                out[row,col] = sum( WEEKS_PER_YEAR .* incd[selected,col] .* incd[selected,:weight] ) # £mn
                row += 1
                out[row,col] = sum((incd[selected,col] .> 0) .* incd[selected,:weight]) # Counts
                row += 1
            end

        end
        select!(out, Not([:pid,:hid,:weight])) # clear out pids
        return out
    end

    function fill_inc_frame_row!(
        ir :: DataFrameRow,
        hh :: Household,
        pers :: Person,
        pres :: IndividualResult,
        from_child_record :: Bool )

        STBIncomes.fill_inc_frame_row!(
            ir, pers.pid, hh.hid, hh.weight, pres.income )
        # some aggregate income fields
        ir.income_tax -=  pres.it.pension_relief_at_source
        ir.employers_ni = pres.ni.class_1_secondary

        ## FIXME the pension_relief thing might not be quite right
        ir.scottish_income_tax = pres.it.non_savings_tax - pres.it.pension_relief_at_source

        ir.total_benefits = isum( pres.income, BENEFITS )
        ir.legacy_mtbs  = isum( pres.income, LEGACY_MTBS )
        ir.means_tested_bens = isum( pres.income, MEANS_TESTED_BENS )
        ir.non_means_tested_bens  = isum( pres.income, NON_MEANS_TESTED_BENS )
        ir.sickness_illness = isum( pres.income, SICKNESS_ILLNESS )
        ir.scottish_benefits = isum( pres.income, SCOTTISH_BENEFITS )
        ir.pension_relief_at_source = pres.it.pension_relief_at_source

        ir.tenure = hh.tenure
        ir.data_year = hh.data_year
        ir.region = hh.region
        ir.gross_decile = -1
        ir.council = hh.council
        ir.employment_status = pers.employment_status
        ir.age_band = age_range( pers.age )
        ir.is_child = from_child_record
        ir.employers_ni = pres.ni.class_1_secondary
        ir.net_cost = isum( pres.income, NET_COST ) + ir.pension_relief_at_source

    end

    function fill_pers_frame_row!(
        pr :: DataFrameRow,
        hh :: Household,
        pers :: Person,
        pres :: IndividualResult,
        from_child_record :: Bool )
        pr.hid = hh.hid
        pr.pid = pers.pid
        pr.data_year = hh.data_year
        pr.weight = hh.weight
        pr.sex = pers.sex
        pr.age_band  = age_range( pers.age )
        pr.employment_status = pers.employment_status
        pr.ethnic_group = pers.ethnic_group
        pr.is_child = from_child_record

        pr.income_taxes = isum( pres.income, INCOME_TAXES)
        pr.means_tested_benefits = isum( pres.income, MEANS_TESTED_BENS )
        pr.other_benefits = isum( pres.income, NON_MEANS_TESTED_BENS )

        pr.income_tax = pres.income[INCOME_TAX]
        pr.it_non_savings = pres.it.non_savings_tax
        pr.it_savings = pres.it.savings_tax
        pr.it_dividends = pres.it.dividends_tax
        pr.it_pension_relief_at_source = pres.it.pension_relief_at_source

        pr.ni_above_lower_earnings_limit = pres.ni.above_lower_earnings_limit
        pr.ni_total_ni = pres.income[NATIONAL_INSURANCE]
        pr.ni_class_1_primary = pres.ni.class_1_primary
        pr.ni_class_1_secondary = pres.ni.class_1_secondary
        pr.ni_class_2  = pres.ni.class_2
        pr.ni_class_3  = pres.ni.class_3
        pr.ni_class_4  = pres.ni.class_4
        pr.assumed_gross_wage = pres.ni.assumed_gross_wage
        if pres.metr != -12345.0 # missing indicator
            pr.metr = pres.metr
        end
        pr.replacement_rate = pres.replacement_rate
    end

    #
    # fill the rows in the output dataframes for this hhld
    #
    function add_to_frames!(
        frames :: NamedTuple,
        hh     :: Household,
        hres   :: HouseholdResult,
        sysno  :: Integer,
        num_systems :: Integer  )

        hfno = get_slot_for_household( hh.hid, hh.data_year )
        fill_hh_frame_row!(
            frames.hh[sysno][hfno, :], hh, hres)
        nbus = length(hres.bus)
        np = length( hh.people )
        bus = get_benefit_units( hh )
        pfbu = 0
        for buno in 1:nbus
            # this won't work at the moment & isn't used
            # fill_bu_frame_row!( frames.bu[sysno][bfno,:], hh, hres.bus[buno])
            for( pid, pers ) in bus[buno].people
                pfno = get_slot_for_person( pid, hh.data_year )
                # pfbu += 1
                from_child_record = pid in bus[buno].children
                fill_pers_frame_row!(
                    frames.indiv[sysno][pfno,:],
                    hh,
                    pers,
                    hres.bus[buno].pers[pid],
                    from_child_record )
                fill_inc_frame_row!(
                    frames.income[sysno][pfno,:],
                    hh,
                    pers,
                    hres.bus[buno].pers[pid],
                    from_child_record )

            end # person loop
        end # bu loop
    end

    """
    got to be a better way...
    """
    function income_measure_as_sym( i :: IneqIncomeMeasure)::Symbol
        return if i == bhc_net_income
            :bhc_net_income
        elseif i == ahc_net_income
            :ahc_net_income
        elseif i == eq_ahc_net_income
            :eq_ahc_net_income
        elseif i == eq_bhc_net_income
            :eq_bhc_net_income
        end
    end

    function eq_income_measure( i :: IneqIncomeMeasure )
        return i in [bhc_net_income, eq_bhc_net_income ] ? eq_bhc_net_income : eq_ahc_net_income
    end

    function make_poverty_line( hhs :: DataFrame, settings :: Settings ) :: Real
        income = income_measure_as_sym( settings.ineq_income_measure )
        deciles = PovertyAndInequalityMeasures.binify( hhs, 10, :weighted_people, income )
        return deciles[5,3]*(2.0/3.0)
    end

    function make_gain_lose( post :: DataFrame, pre :: DataFrame, measure :: Symbol )::NamedTuple
        pre_inc = pre[:,measure]
        post_inc = post[:,measure]
        n = size(post_inc)[1]
        diff = zeros(n)
        # so check for 0 pre income
        for i in 1:n
            if pre_inc[i] != 0
                diff[i] = (pre_inc[i]-post_inc[i])/pre_inc[i]
            end
        end
        # println( "diff=$diff n=$n")
        popn = sum(post[:,:weighted_people])
        gainers = (diff.>=0.01).*post[:,:weighted_people]
        sg = sum(gainers)
        losers = (diff.<= -0.01).*post[:,:weighted_people]
        sl = sum(losers)
        nc = (abs.(diff).< 0.01).*post[:,:weighted_people]
        sn = sum(nc)
        @assert sg+sl+sn ≈ popn "sg=$sg + sl=$sl + sn=$sn !≈ popn=$popn"
        return ( gainers=sg, losers=sl,nc=sn, popn = popn )
    end

    function make_gain_lose( post :: DataFrame, pre :: DataFrame, settings :: Settings )::NamedTuple
        income = income_measure_as_sym( settings.ineq_income_measure )
        return make_gain_lose( post, pre, income )
    end

    function metrs_to_hist( indiv :: DataFrame ) :: NamedTuple
        # these 2 convoluted lines make this draw only
        # over the non-missing (children, retired)
        p = collect(keys(skipmissing( indiv.metr )))
        indp = indiv[p,[:metr, :weight]] # just non missing
        # so .. <=0, >0 <=10, >10<=20 and so on
        # skip near-infinite mrs mwhen averaging
        mmtr = mean( indp[(indp.metr.<150),:].metr, Weights(indp[(indp.metr.<150),:].weight))
        hist = fit( Histogram, indp.metr, Weights( indp.weight ), [-Inf, 0.00001, 10.0, 20.0, 30.0, 50.0, 80.0, 100.0, Inf], closed=:right )
        return ( mean=mmtr, hist=hist)
    end

    function summarise_frames(
        frames :: NamedTuple,
        settings :: Settings ) :: NamedTuple
        ns = size( frames.indiv )[1] # num systems
        income_summary = []
        gain_lose = []
        poverty = []
        inequality = []
        deciles = []
        quantiles = []
        metrs = []
        poverty_line = []
        child_poverty = []
        income_measure = income_measure_as_sym( settings.ineq_income_measure )

        for sysno in 1:ns
            push!( metrs, metrs_to_hist( frames.indiv[sysno] ))
            push!(income_summary,
                summarise_inc_frame(frames.income[sysno]))
            push!( deciles,
                PovertyAndInequalityMeasures.binify(
                    frames.hh[sysno],
                    10,
                    :weighted_people,
                    income_measure ))
            push!( quantiles,
                PovertyAndInequalityMeasures.binify(
                    frames.hh[sysno],
                    50,
                    :weighted_people,
                    income_measure  ))

            ineq = make_inequality(
                frames.hh[sysno],
                :weighted_people,
                income_measure  )
            push!( inequality, ineq )
            push!(
                poverty,
                PovertyAndInequalityMeasures.make_poverty(
                    frames.hh[sysno],
                    settings.poverty_line,
                    settings.growth,
                    :weighted_people,
                    income_measure ))
            push!( poverty_line, settings.poverty_line )
            cp = calc_child_poverty(
                settings.poverty_line,
                frames.hh[sysno],
                measure=income_measure
            )
            push!( child_poverty, cp )
        end
        return (
            quantiles=quantiles,
            deciles = deciles,
            income_summary = income_summary,
            poverty=poverty,
            inequality=inequality,
            metrs = metrs,
            child_poverty=child_poverty,
            poverty_line = settings.poverty_line )
    end

    ## FIXME eventually, move this to DrWatson
    function dump_frames(
        settings :: Settings,
        frames :: NamedTuple )
        ns = size( frames.indiv )[1] # num systems
        fbase = basiccensor(settings.run_name)
        mkpath(settings.output_dir)
        for fno in 1:ns
            fname = "$(settings.output_dir)/$(fbase)_$(fno)_hh.csv"
            CSV.write( fname, frames.hh[fno] )
            fname = "$(settings.output_dir)/$(fbase)_$(fno)_bu.csv"
            CSV.write( fname, frames.bu[fno] )
            fname = "$(settings.output_dir)/$(fbase)_$(fno)_pers.csv"
            CSV.write( fname, frames.indiv[fno] )
            fname = "$(settings.output_dir)/$(fbase)_$(fno)_income.csv"
            CSV.write( fname, frames.income[fno] )
            income_summary = summarise_inc_frame(frames.income[fno])
            fname = "$(settings.output_dir)/$(fbase)_$(fno)_income-summary.csv"
            CSV.write( fname, income_summary )

        end
    end

end # STBOutput

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using Arrow, DataFrames
using Statistics, Distributions, StatsBase
#using Roots
using Optim, Roots
using PyPlot, Printf, PyCall

include("variables.jl")
include("plot_styles.jl")
using .Vars: classes

measures = DataFrame()
basedir = "/work/sschaub/JuliaForHEP/final_plotting/"
tex = true

for i in 1:10, feature in filter(x -> endswith(x, "_$i") && isdir(joinpath(basedir, x)), readdir(basedir))
    output_dir = joinpath(basedir, "$feature/")
    feature = replace(feature, Regex("_$i\$") => "")
    @show feature, i

    all_features = DataFrame(Arrow.Table(joinpath(output_dir, "all_features.arrow")))
    transform!(all_features,
        AsTable(r"output_predicted_.+") => ByRow(x -> Symbol(split(string(argmax(x)), '_')[end])) => :prediction,
    )

    for ((kind,), df) in pairs(groupby(all_features, :kind))
        #kind === :validation || continue
        split = (train=.75 * .8, test=.25, validation=.75 * .2)[kind]

        hists = DataFrame()
        for ((node, true_class), df) in pairs(groupby(df, [:prediction, :output_expected], sort=true))
            data = df[!, Symbol(:output_predicted_, node)]
            weights = Weights(df.weights ./ split)# .* (300 / 41.5))
            bins = range(.25, 1; length=16)
            values = fit(Histogram, data, weights, bins).weights
            append!(hists, DataFrame(; bins=bins[1:end-1], values, node, true_class))
        end
        @show combine(groupby(hists, [:true_class, :node]), :values => sum)
        hists = groupby(hists, [:bins, :node])
        hists = combine(hists) do x
            sig = x.true_class .=== :ttH
            #(signal=round(Int, sum(x[sig, :values])), bg=round(Int, sum(x[Not(sig), :values])))
            (signal=sum(x[sig, :values]), bg=sum(x[Not(sig), :values]))
        end
        #filter!([:signal, :bg] => (s, b) -> s + b > 0, hists)
        k = round.(Int, hists.signal .+ hists.bg)
        λ(μ) = μ .* hists.signal .+ hists.bg
        NLL(μ) = -sum(logpdf.(Poisson.(λ(μ)), k))
        μ_exp = Optim.minimizer(optimize(NLL, 0., 2.))
        z1, z2 = find_zeros(μ -> 2(NLL(μ) - NLL(μ_exp)) - 1, 0, 2)
        push!(measures, (; feature, kind, σ_μ = (z2 - z1) / 2, i))

        μ = range(.5, 1.5; length=100)
        t = NLL.(μ) .- NLL(μ_exp)

        fig, ax = subplots()
        ax.plot(μ, t)
        ax.axhline(.5; color=:grey, ls="--", lw=1)
        ax.axvline.([z1, z2]; color=:grey, ls="--", lw=1)
        fig.suptitle("Likelihood Fit ttH")
        tex && ax.set_title("\\verb|$feature|"; fontsize=16)
        ax.set_xlabel(L"\mu")
        ax.set_xlim(.5, 1.5)
        #tex && ax.set_ylabel(L"-2\log(\mathcal L(\mu) / \mathcal L(\langle \mu \rangle))")
        tex && ax.set_ylabel(L"\mathrm{NLL}(\mu) - \mathrm{NLL}(\langle \mu \rangle)")
        ax.legend([@sprintf("\$\\mu = %.3f^{+%.3f}_{-%.3f}\$", μ_exp, μ_exp-z1, z2-μ_exp)]; loc="upper center")
        annotate_cms(ax)
        display(fig)
        fig.savefig(joinpath(output_dir, "likelihood_ttH_$(kind)_15bins.pdf"))

        fig1, axs1 = subplots(ncols=2, nrows=2, figsize=(16, 14))
        foreach(
            pairs(groupby(hists, :node, sort=true)),
            axs1,
        ) do ((node,), df), ax1
            ax1.set_title("$node Node")
            #ax1.step(df.bins, df.bg; label="background")
            ax1.step(df.bins, df.signal; label="ttH Signal")
            ax1.step(df.bins, df.bg + df.signal; label="Total")
            ax1.legend()
            ax1.set_xlabel("P($node)")
            ax1.set_yscale(:log)
            #ax1.set_ylim(minimum(filter(!=(0), df.bg)) * .8, maximum(df.bg + df.signal) * 2)
            ax1.set_ylim(.05, maximum(df.bg + df.signal) * 2)
            ax1.set_ylabel("Events/Bin")
            annotate_cms(ax1)
        end
        tex && fig1.suptitle("Binning \\verb|$feature|")
        fig1.tight_layout()
        display(fig1)
        fig1.savefig(joinpath(output_dir, "binning_ttH_$(kind)_15bins.pdf"))
    end
end

#replace!(measures.feature, "lola+scalars11" => "lola+all scalars", "scalars2" => "only scalars")
Arrow.write(
    joinpath(basedir, "stddevs_mu_15bins.arrow"),
    measures,
)

begin
idx = sortperm(measures.feature; by=x -> findfirst(==(x), ["lola+" .* ["none", "tier3", "tier2+3", "all"]; "jets_as_scalars+all"; "none+all"]))
measures = measures[idx, :]
fig, ax = subplots(figsize=(11, 11))
marker = (train=:v, test=:P, validation=:o)
foreach(pairs(groupby(measures, :kind))) do ((kind,), df)
    kind === :validation || return
    df = combine(groupby(df, :feature), (:σ_μ .=> (mean, std))...)
    x = axes(df, 1)
    ax.errorbar(
        x, df.σ_μ_mean; yerr=df.σ_μ_std,
        color=:blue#=[:red; fill(:blue, length(x)-3); :green; :orange]=#, marker=marker[kind],
        lw=0, elinewidth=1, capsize=8, markersize=8,
        label=kind,
    )
end

x = unique(measures.feature)
ax.set_xticks(eachindex(x))
ax.set_xticklabels(string.("\\verb|", x, "|"); rotation=90)
ax.legend(; loc="upper right", fontsize=16, frameon=true)
ax.set_ylabel(L"\sigma_\mu")
fig.suptitle("Standard Deviations of \$\\mu\$ for Training with LoLa + X")
annotate_cms(ax)
fig.tight_layout()
fig.savefig(joinpath(basedir, "stddevs_mu_15bins.pdf"))
display(fig)
end


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

mutable struct IoK8sApiCoreV1FCVolumeSource <: SwaggerModel
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    lun::Any # spec type: Union{ Nothing, Int32 } # spec name: lun
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    targetWWNs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: targetWWNs
    wwids::Any # spec type: Union{ Nothing, Vector{String} } # spec name: wwids

    function IoK8sApiCoreV1FCVolumeSource(;fsType=nothing, lun=nothing, readOnly=nothing, targetWWNs=nothing, wwids=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1FCVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1FCVolumeSource, Symbol("lun"), lun)
        setfield!(o, Symbol("lun"), lun)
        validate_property(IoK8sApiCoreV1FCVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1FCVolumeSource, Symbol("targetWWNs"), targetWWNs)
        setfield!(o, Symbol("targetWWNs"), targetWWNs)
        validate_property(IoK8sApiCoreV1FCVolumeSource, Symbol("wwids"), wwids)
        setfield!(o, Symbol("wwids"), wwids)
        o
    end
end # type IoK8sApiCoreV1FCVolumeSource

const _property_map_IoK8sApiCoreV1FCVolumeSource = Dict{Symbol,Symbol}(Symbol("fsType")=>Symbol("fsType"), Symbol("lun")=>Symbol("lun"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("targetWWNs")=>Symbol("targetWWNs"), Symbol("wwids")=>Symbol("wwids"))
const _property_types_IoK8sApiCoreV1FCVolumeSource = Dict{Symbol,String}(Symbol("fsType")=>"String", Symbol("lun")=>"Int32", Symbol("readOnly")=>"Bool", Symbol("targetWWNs")=>"Vector{String}", Symbol("wwids")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1FCVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1FCVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1FCVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1FCVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1FCVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1FCVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1FCVolumeSource)
    true
end

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


text:
export PARAMETER_DEFAULT, PARAMETER_UNSTABLE_BUT_FAST, PARAMETER_STABLE_BUT_SLOW

using MathOptInterface
MOI = MathOptInterface
const MOIU = MOI.Utilities
const AFFEQ = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Cdouble}, MOI.EqualTo{Cdouble}}

mutable struct Optimizer <: MOI.AbstractOptimizer
    objconstant::Cdouble
    objsign::Int
    blockdims::Vector{Int}
    varmap::Vector{Tuple{Int, Int, Int}} # Variable Index vi -> blk, i, j
    b::Vector{Cdouble}
    problem::Union{Nothing, SDPAProblem}
    solve_time::Float64
    silent::Bool
    options::Dict{Symbol, Any}
    function Optimizer(; kwargs...)
		optimizer = new(
            zero(Cdouble), 1, Int[], Tuple{Int, Int, Int}[], Cdouble[],
            nothing, NaN, false, Dict{Symbol, Any}())
		for (key, value) in kwargs
			MOI.set(optimizer, MOI.RawParameter(key), value)
		end
		return optimizer
    end
end

varmap(optimizer::Optimizer, vi::MOI.VariableIndex) = optimizer.varmap[vi.value]

function MOI.supports(optimizer::Optimizer, param::MOI.RawParameter)
	return param.name in keys(SET_PARAM)
end
function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value)
	if !MOI.supports(optimizer, param)
		throw(MOI.UnsupportedAttribute(param))
	end
	optimizer.options[param.name] = value
end
function MOI.get(optimizer::Optimizer, param::MOI.RawParameter)
	# TODO: This gives a poor error message if the name of the parameter is invalid.
	return optimizer.options[param.name]
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
	optimizer.silent = value
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

MOI.get(::Optimizer, ::MOI.SolverName) = "SDPA"

# See https://www.researchgate.net/publication/247456489_SDPA_SemiDefinite_Programming_Algorithm_User's_Manual_-_Version_600
# "SDPA (SemiDefinite Programming Algorithm) User's Manual — Version 6.00" Section 6.2
const RAW_STATUS = Dict(
    noINFO     => "The iteration has exceeded the maxIteration and stopped with no informationon the primal feasibility and the dual feasibility.",
    pdOPT      => "The normal termination yielding both primal and dual approximate optimal solutions.",
    pFEAS      => "The primal problem got feasible but the iteration has exceeded the maxIteration and stopped.",
    dFEAS      => "The dual problem got feasible but the iteration has exceeded the maxIteration and stopped.",
    pdFEAS     => "Both primal problem and the dual problem got feasible, but the iterationhas exceeded the maxIteration and stopped.",
    pdINF      => "At least one of the primal problem and the dual problem is expected to be infeasible.",
    pFEAS_dINF => "The primal problem has become feasible but the dual problem is expected to be infeasible.",
    pINF_dFEAS => "The dual problem has become feasible but the primal problem is expected to be infeasible.",
    pUNBD      => "The primal problem is expected to be unbounded.",
    dUNBD      => "The dual problem is expected to be unbounded.")

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    if optimizer.problem === nothing
        return "`MOI.optimize!` not called."
    end
	return RAW_STATUS[getPhaseValue(optimizer.problem)]
end
function MOI.get(optimizer::Optimizer, ::MOI.SolveTime)
	return optimizer.solve_time
end

function MOI.is_empty(optimizer::Optimizer)
    return iszero(optimizer.objconstant) &&
        optimizer.objsign == 1 &&
        isempty(optimizer.blockdims) &&
        isempty(optimizer.varmap) &&
        isempty(optimizer.b)
end
function MOI.empty!(optimizer::Optimizer)
    optimizer.objconstant = zero(Cdouble)
    optimizer.objsign = 1
    empty!(optimizer.blockdims)
    empty!(optimizer.varmap)
    empty!(optimizer.b)
    optimizer.problem = nothing
end

function MOI.supports(
    optimizer::Optimizer,
    ::Union{MOI.ObjectiveSense,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Cdouble}}})
    return true
end

function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.VectorOfVariables}, ::Type{MOI.Reals})
    return false
end
const SupportedSets = Union{MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle}
function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.VectorOfVariables},
    ::Type{<:SupportedSets})
    return true
end
function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.ScalarAffineFunction{Cdouble}},
    ::Type{MOI.EqualTo{Cdouble}})
    return true
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end
MOIU.supports_allocate_load(::Optimizer, copy_names::Bool) = !copy_names

function MOIU.allocate(optimizer::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    # To be sure that it is done before load(optimizer, ::ObjectiveFunction, ...), we do it in allocate
    optimizer.objsign = sense == MOI.MIN_SENSE ? -1 : 1
end
function MOIU.allocate(::Optimizer, ::MOI.ObjectiveFunction, ::MOI.ScalarAffineFunction) end

function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense) end
# Loads objective coefficient α * vi
function load_objective_term!(optimizer::Optimizer, α, vi::MOI.VariableIndex)
    blk, i, j = varmap(optimizer, vi)
    coef = optimizer.objsign * α
    if i != j
        coef /= 2
    end
    # in SDP format, it is max and in MPB Conic format it is min
    inputElement(optimizer.problem, 0, blk, i, j, float(coef), false)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction, f::MOI.ScalarAffineFunction)
    obj = MOIU.canonical(f)
    optimizer.objconstant = f.constant
    for t in obj.terms
        if !iszero(t.coefficient)
            load_objective_term!(optimizer, t.coefficient, t.variable_index)
        end
    end
end

function new_block(optimizer::Optimizer, set::MOI.Nonnegatives)
    push!(optimizer.blockdims, -MOI.dimension(set))
    blk = length(optimizer.blockdims)
    for i in 1:MOI.dimension(set)
        push!(optimizer.varmap, (blk, i, i))
    end
end

function new_block(optimizer::Optimizer, set::MOI.PositiveSemidefiniteConeTriangle)
    push!(optimizer.blockdims, set.side_dimension)
    blk = length(optimizer.blockdims)
    for i in 1:set.side_dimension
        for j in 1:i
            push!(optimizer.varmap, (blk, i, j))
        end
    end
end

function MOIU.allocate_constrained_variables(optimizer::Optimizer,
                                             set::SupportedSets)
    offset = length(optimizer.varmap)
    new_block(optimizer, set)
    ci = MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}(offset + 1)
    return [MOI.VariableIndex(i) for i in offset .+ (1:MOI.dimension(set))], ci
end

function MOIU.load_constrained_variables(
    optimizer::Optimizer, vis::Vector{MOI.VariableIndex},
    ci::MOI.ConstraintIndex{MOI.VectorOfVariables},
    set::SupportedSets)
end

function MOIU.load_variables(optimizer::Optimizer, nvars)
    @assert nvars == length(optimizer.varmap)
    dummy = isempty(optimizer.b)
    if dummy
        optimizer.b = [one(Cdouble)]
        optimizer.blockdims = [optimizer.blockdims; -1]
    end
    optimizer.problem = SDPAProblem()
    setParameterType(optimizer.problem, PARAMETER_DEFAULT)
	# TODO Take `silent` into account here
    setparameters!(optimizer.problem, optimizer.options)
    inputConstraintNumber(optimizer.problem, length(optimizer.b))
    inputBlockNumber(optimizer.problem, length(optimizer.blockdims))
    for (i, blkdim) in enumerate(optimizer.blockdims)
        inputBlockSize(optimizer.problem, i, blkdim)
        inputBlockType(optimizer.problem, i, blkdim < 0 ? LP : SDP)
    end
    initializeUpperTriangleSpace(optimizer.problem)
    for i in eachindex(optimizer.b)
        inputCVec(optimizer.problem, i, optimizer.b[i])
    end
    if dummy
        inputElement(optimizer.problem, 1, length(optimizer.blockdims), 1, 1, one(Cdouble), false)
    end
end

function MOIU.allocate_constraint(optimizer::Optimizer,
                                  func::MOI.ScalarAffineFunction{Cdouble},
                                  set::MOI.EqualTo{Cdouble})
    push!(optimizer.b, MOI.constant(set))
    return AFFEQ(length(optimizer.b))
end

function MOIU.load_constraint(m::Optimizer, ci::AFFEQ,
                              f::MOI.ScalarAffineFunction, s::MOI.EqualTo)
    if !iszero(MOI.constant(f))
        throw(MOI.ScalarFunctionConstantNotZero{
            Cdouble, MOI.ScalarAffineFunction{Cdouble}, MOI.EqualTo{Cdouble}}(
                MOI.constant(f)))
    end
    f = MOIU.canonical(f) # sum terms with same variables and same outputindex
    for t in f.terms
        if !iszero(t.coefficient)
            blk, i, j = varmap(m, t.variable_index)
            coef = t.coefficient
            if i != j
                coef /= 2
            end
            inputElement(m.problem, ci.value, blk, i, j, float(coef), false)
        end
    end
end

function MOI.optimize!(m::Optimizer)
	start_time = time()
    SDPA.initializeUpperTriangle(m.problem, false)
    SDPA.initializeSolve(m.problem)
    SDPA.solve(m.problem)
    m.solve_time = time() - start_time
end

function MOI.get(m::Optimizer, ::MOI.TerminationStatus)
    if m.problem === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    end
    status = getPhaseValue(m.problem)
    if status == noINFO
        return MOI.ITERATION_LIMIT
    elseif status == pFEAS
        return MOI.SLOW_PROGRESS
    elseif status == dFEAS
        return MOI.SLOW_PROGRESS
    elseif status == pdFEAS
        return MOI.OPTIMAL
    elseif status == pdINF
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif status == pFEAS_dINF
        return MOI.DUAL_INFEASIBLE
    elseif status == pINF_dFEAS
        return MOI.INFEASIBLE
    elseif status == pdOPT
        return MOI.OPTIMAL
    elseif status == pUNBD
        return MOI.DUAL_INFEASIBLE
    elseif status == dUNBD
        return MOI.INFEASIBLE
    end
end

function MOI.get(m::Optimizer, attr::MOI.PrimalStatus)
    if attr.N > MOI.get(m, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    status = getPhaseValue(m.problem)
    if status == noINFO
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == pFEAS
        return MOI.FEASIBLE_POINT
    elseif status == dFEAS
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == pdFEAS
        return MOI.FEASIBLE_POINT
    elseif status == pdINF
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == pFEAS_dINF
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == pINF_dFEAS
        return MOI.INFEASIBLE_POINT
    elseif status == pdOPT
        return MOI.FEASIBLE_POINT
    elseif status == pUNBD
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == dUNBD
        return MOI.INFEASIBLE_POINT
    end
end

function MOI.get(m::Optimizer, attr::MOI.DualStatus)
    if attr.N > MOI.get(m, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    status = getPhaseValue(m.problem)
    if status == noINFO
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == pFEAS
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == dFEAS
        return MOI.FEASIBLE_POINT
    elseif status == pdFEAS
        return MOI.FEASIBLE_POINT
    elseif status == pdINF
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == pFEAS_dINF
        return MOI.INFEASIBLE_POINT
    elseif status == pINF_dFEAS
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == pdOPT
        return MOI.FEASIBLE_POINT
    elseif status == pUNBD
        return MOI.INFEASIBLE_POINT
    elseif status == dUNBD
        return MOI.INFEASIBILITY_CERTIFICATE
    end
end

MOI.get(m::Optimizer, ::MOI.ResultCount) = m.problem === nothing ? 0 : 1
function MOI.get(m::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(m, attr)
    return m.objsign * getPrimalObj(m.problem) + m.objconstant
end
function MOI.get(m::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(m, attr)
    return m.objsign * getDualObj(m.problem) + m.objconstant
end
struct PrimalSolutionMatrix <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::PrimalSolutionMatrix) = true
MOI.get(optimizer::Optimizer, ::PrimalSolutionMatrix) = PrimalSolution(optimizer.problem)

struct DualSolutionVector <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::DualSolutionVector) = true
function MOI.get(optimizer::Optimizer, ::DualSolutionVector)
    return unsafe_wrap(Array, getResultXVec(optimizer.problem), getConstraintNumber(optimizer.problem))
end

struct DualSlackMatrix <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::DualSlackMatrix) = true
MOI.get(optimizer::Optimizer, ::DualSlackMatrix) = VarDualSolution(optimizer.problem)

function block(optimizer::Optimizer, ci::MOI.ConstraintIndex{MOI.VectorOfVariables})
    return optimizer.varmap[ci.value][1]
end
function dimension(optimizer::Optimizer, ci::MOI.ConstraintIndex{MOI.VectorOfVariables})
    blockdim = optimizer.blockdims[block(optimizer, ci)]
    if blockdim < 0
        return -blockdim
    else
        return MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(blockdim))
    end
end
function vectorize_block(M, blk::Integer, s::Type{MOI.Nonnegatives})
    return diag(block(M, blk))
end
function vectorize_block(M::AbstractMatrix{Cdouble}, blk::Integer, s::Type{MOI.PositiveSemidefiniteConeTriangle}) where T
    B = block(M, blk)
    d = LinearAlgebra.checksquare(B)
    n = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(d))
    v = Vector{Cdouble}(undef, n)
    k = 0
    for j in 1:d
        for i in 1:j
            k += 1
            v[k] = B[i, j]
        end
    end
    @assert k == n
    return v
end

function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex)
    MOI.check_result_index_bounds(optimizer, attr)
    blk, i, j = varmap(optimizer, vi)
    return block(MOI.get(optimizer, PrimalSolutionMatrix()), blk)[i, j]
end

function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
    MOI.check_result_index_bounds(optimizer, attr)
    return vectorize_block(MOI.get(optimizer, PrimalSolutionMatrix()), block(optimizer, ci), S)
end
function MOI.get(m::Optimizer, attr::MOI.ConstraintPrimal, ci::AFFEQ)
    MOI.check_result_index_bounds(m, attr)
    return m.b[ci.value]
end

function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
    MOI.check_result_index_bounds(optimizer, attr)
    return vectorize_block(MOI.get(optimizer, DualSlackMatrix()), block(optimizer, ci), S)
end
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual, ci::AFFEQ)
    MOI.check_result_index_bounds(optimizer, attr)
    return -MOI.get(optimizer, DualSolutionVector())[ci.value]
end


text:
export UpperConfidenceBound, EpsilonGreedySelector, WeightedSample, AlternateSelector

using StatsBase:weights, sample

include("epsilon_greedy.jl")
include("upper_confidence_bound.jl")
include("weighted_sample.jl")
include("alternate.jl")


text:
using ModelingStructs: maybe_parent, indices, recursive_length

@testset "Utilities" begin
    @testset "maybe_parent" begin
        a = [1,2,3,4]
        b = view(a, 1:3)
        c = view(b, 2)
        @test maybe_parent(b) === a
        @test maybe_parent(c) === a
    end

    @testset "indices" begin
        lengths = [1,3,6,1]
        inds = indices(lengths)
        @test inds == (1:1, 2:4, 5:10, 11:11)
        @test last(last(inds)) == sum(lengths)
    end

    @testset "recursive_length" begin
        single = 1
        vect = [1, 2, 3, 4]
        nest_vect = [1, 2, [3, 4]]
        nt = (a=1, b=2)
        nest_nt = (a=single, b=nt, c=[nt, nt], d=(a=vect, b=nest_vect), e=(a=nt, b=nest_vect))
        @test recursive_length(single) == 1
        @test recursive_length(vect) == 4
        @test recursive_length(nest_vect) == 4
        @test recursive_length(nt) == 2
        @test recursive_length(nest_nt) == 21
    end
end

max_stars_repo_name: jonniedie/NamedViewVectors


text:
using EmpiricalBayes
using Test
using Distributions

true_dist = MixtureModel([ Normal(-0.3,.5), Normal(1.05,.5)])
marginal_grid = collect(range(-6, stop=6, length=1001));

marginal_h = marginal_grid[2] - marginal_grid[1] # should really replace by range

prior_grid = collect(range(-3, stop=3, length=121));
d_true = NormalConvolutionProblem(true_dist, marginal_grid);

ds = MixingNormalConvolutionProblem(Normal, 0.2, prior_grid, marginal_grid);

m = 200

cb_cb = ComteButucea(MarginalDensityTarget(0.0), m, marginal_grid)
cb_kdecalib = KDECalibrator(SincKernel, m, marginal_grid)

max_bias_cb = check_bias(cb_cb, ds; maximization=true)

max_bias_cb_kdecalib = check_bias(cb_kdecalib, ds; maximization=true)

@test max_bias_cb ≈ max_bias_cb_kdecalib atol=0.0001


text:
# Fourier Transform and visualization of which

using Plots
# idk why I just like to load plot function first
# - weired isn't it
gr()
# rotation factor

W_rotate(N::Int)=exp(-2*im*π/N)

W_rotate.(1:10)
plot(1:100,W_rotate.(1:100))#takes long if runned for the 1st time
plot(W_rotate.(1:100))#

#Amplitute can be calcd with abs(x::Complex)
abs(1+im)

# DFT and visualization

function dft(x::Array)
    n = length(x)
    r = Array{Any,1}(missing,n)
    #Note how sum is used to aggregate results after broadcasting
    for i in 1:n
        r[i]=sum(x.*exp(1).^(-(i-1)*collect(0:(n-1)).*(2*π*im/n)))
    end
    return(r)
end

x=1:.1:2π
y =sin.(x)
plot(x,y)
plot(1:length(x),dft(y))
plot(1:length(x),abs.(dft(y)))
plot(1:length(x),real.(dft(y)))
# FFT and visualization of FFT

function idft(x::Array)
    n = length(x)
    r = Array{Any,1}(missing,n)
    #Note how sum is used to aggregate results after broadcasting
    for i in 1:n
        r[i]=sum(x.*exp(1).^(+(i-1)*collect(0:(n-1)).*(2*π*im/n)))#change - to+
    end
    return(r)
end

#below inspired by
#https://www.mathworks.com/matlabcentral/answers/3187-fft-recursive-code-problem

function p2_fft_rec(f)
  n = length(f)
  if (n == 1)
    return(f)
  else
    a= p2_fft_rec(f[1:2:n])
    b= p2_fft_rec(f[2:2:n]).*exp(1).^(-2/n*π*im.*collect(0:1:(n/2)-1))
    return append!(a+b,a-b)
  end
end

function p2_fft(x)
  n = length(x)
  N = 2^(ceil(log2(n)))
  x = append!(x,repeat([0],Int64(N - n)))
  return(p2_fft_rec(x)[1:n])
end

#p2_fft(collect(1:4))
#inspired by
#https://in.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html

function p2_psd_fft(x)
  N = length(x)
  xdft = p2_fft(x)[1:Int(N/2+1)]
  psdx = (1/N/N) * abs.(xdft).^2
  psdx[2:end] = 2*psdx[2:end]
  freq = 0:1:(N/2)
  return(Dict("freq"=>freq,"psdx"=>psdx,"psd_est"=>10*log10.(psdx)))
end

max_stars_repo_path: Fourier Transform and power spectral.jl

max_issues_repo_path: Fourier Transform and power spectral.jl

max_forks_repo_path: Fourier Transform and power spectral.jl


text:
# 6.3s.jl

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';');
scale!(df, [:kcal_per_g, :perc_fat, :perc_lactose])

stan6_3 = "
data{
    int <lower=1> N;              // Sample size
    vector[N] K;
    vector[N] F;
}
parameters{
    real a;
    real bF;
    real<lower=0> sigma;
}
model{
    vector[N] mu;
    sigma ~ exponential( 1 );
    a ~ normal( 0 , 0.2 );
    bF ~ normal( 0 , 0.5 );
    mu = a + bF * F;
    K ~ normal( mu , sigma );
}
";

# Define the SampleModel and set the output format to :mcmcchains.

begin
    m6_3s = SampleModel("m6.3", stan6_3);
    m6_3_data = Dict("N" => size(df, 1), "F" => df.perc_fat_s, "K" => df.kcal_per_g_s);
    rc6_3s = stan_sample(m6_3s, data=m6_3_data);
end

if success(rc6_3s)
    part6_3s = read_samples(m6_3s, :particles)
    part6_3s |> display
end

# End of 6.3s.jl


text:
module Mods

import Base.isequal, Base.==, Base.+, Base.-, Base.*
import Base.inv, Base./, Base.^
import Base.hash, Base.show

export Mod
export isequal, ==, +, -, *
export is_invertible, inv, /, ^
export hash, CRT

"""
`Mod(v,m)` creates a modular number in mod `m` with value `v%m`.
`Mod(m)` is equivalent to `Mod(0,m)`.
"""
struct Mod{T<:Integer}
    val::T
    mod::T
    function Mod(a::S, m::T) where {S<:Integer, T<:Integer}
        if m < 1
            error("Modulus must be at least 1")
        end

        a,m = promote(a,m)
        return new{typeof(a)}(mod(a,m),m)
    end
end

Mod(m::Integer) = Mod(0,m)

function hash(x::Mod, h::UInt64= UInt64(0))
    v = BigInt(x.val)
    m = BigInt(x.mod)
    return hash(v,hash(m,h))
end

# Test for equality
isequal(x::Mod, y::Mod) = x.mod==y.mod && x.val==y.val
==(x::Mod,y::Mod) = isequal(x,y)

function modcheck(x::Mod, y::Mod)
    if x.mod != y.mod
        error("Cannot operate on two Mod objects with different moduli")
    end
    true
end

# Easy arithmetic
function +(x::Mod, y::Mod)
    modcheck(x,y)
    s,flag = Base.add_with_overflow(x.val,y.val)
    if !flag
        return Mod(x.val+y.val, x.mod)
    end
    s = widen(x.val) + widen(y.val)    # add with added precision
    s = mod(s,x.mod)                   # reduce by modulus
    return Mod(oftype(x.mod,s),x.mod)  # return with proper type
end

function -(x::Mod)
    return Mod(-x.val, x.mod)
end

-(x::Mod,y::Mod) = x + (-y)

function *(x::Mod, y::Mod)
    modcheck(x,y)
    p,flag = Base.mul_with_overflow(x.val,y.val)
    if !flag
        return Mod(x.val*y.val, x.mod)
    end
    p = widemul(x.val, y.val)         # multipy with added precision
    p = mod(p,x.mod)                  # reduce by the modulus
    return Mod(oftype(x.mod,p),x.mod) # return with proper type
end

# Division stuff
"""
`is_invertible(x::Mod)` determines if `x` is invertible.
"""
is_invertible(x::Mod) = return gcd(x.val,x.mod)==1 && x.mod>1

"""
`inv(x::Mod)` gives the multiplicative inverse of `x`.
This may be abbreviated by `x'`.
"""
function inv(x::Mod)
    (g, v, ignore) = gcdx(x.val, x.mod)
    if g != 1 || x.mod==1
        error(x, " is not invertible")
    end
    return Mod(v, x.mod)
end

# Typing shortcut for inv(x)
# adjoint(x::Mod)

function /(x::Mod, y::Mod)
    modcheck(x,y)
    return x * inv(y)
end

function ^(x::Mod, k::Integer)
    if k>0
        return Mod(powermod(x.val, k, x.mod), x.mod)
    end
    if k==0
        T = typeof(x.val)
        return Mod(one(T), x.mod)
    end
    y = inv(x)
    return y^(-k)
end

# Operations with Integers

+(x::Mod, k::Integer) = Mod(k,x.mod)+x
+(k::Integer, x::Mod) = x+k

-(x::Mod, k::Integer) = x + (-k)
-(k::Integer, x::Mod) = (-x) + k

*(x::Mod, k::Integer) = Mod(k,x.mod) * x
*(k::Integer, x::Mod) = x*k

/(x::Mod, k::Integer) = x / Mod(k, x.mod)
/(k::Integer, x::Mod) = Mod(k, x.mod) / x

# Comparison with Integers

isequal(x::Mod, k::Integer) = mod(k,x.mod) == x.val
isequal(k::Integer, x::Mod) = isequal(x,k)
==(x::Mod, k::Integer) = isequal(x,k)
==(k::Integer, x::Mod) = isequal(x,k)

# Chinese remainder theorem functions

# private helper function
function CRT_work(x::Mod, y::Mod)
    n = x.mod
    m = y.mod
    if gcd(n,m) != 1
        error("Moduli must be pairwise relatively prime")
    end

    a = x.val
    b = y.val

    k = inv(Mod(n,m)) * (b-a)

    z = a + k.val*n

    return Mod(z, n*m)
end

# public interface
"""
`CRT(m1,m2,...)`: Chinese Remainder Theorem
```
julia> CRT( Mod(4,11), Mod(8,14) )
Mods.Mod(92,154)

julia> 92%11
4

julia> 92%14
8
```
"""
function CRT(mtuple::Mod...)
    n = length(mtuple)
    if n == 0
        return Mod(1)
    end

    result::Mod = mtuple[1]

    for k=2:n
        result = CRT_work(result,mtuple[k])
    end

    return result
end

show(io::IO, m::Mod) = print(io,"Mod(",m.val,",",m.mod,")")

end # end of module Mods