text:
function algpoly(plist::Vector{T}, f, n) where {T}
    @debug "Algorithm [algpoly]" plist f
    R = parent(n)
    # Construct polynomials qj that arise from treating the shift operator
	# as the difference operator plus the identity operator.
    r = length(plist) - 1
    qlist = [sum(binomial(i, j) * plist[i+1] for i in j:r) for j in 0:r]

    # Find all candidates for the degree bound on output polynomial.
    b = maximum([(Nemo.degree(q) - (j - 1)) for (j, q) in enumerate(qlist)])
    alpha = sum(lead(q)*fallfactorial(n, j-1) for (j, q) in enumerate(qlist) if Nemo.degree(q) - (j-1) == b)
    roots = [-Nemo.coeff(f, 0) for (f, _) in Nemo.factor(R(alpha))]
    roots = [r for r in roots if r >= 0 && isone(denominator(r))]

    first = Nemo.degree(f) - b
    second = -b - 1
    third = maximum([roots; 0])
    d = convert(Int, numerator(max(first, second, third)))
    @debug "Degree bound" first second third

    # Use method of undetermined coefficients to find the output polynomial.
    R, varlist = PolynomialRing(base_ring(n), ["x$i" for i in 1:d+1])
    arg = change_base_ring(R, n)
    genpoly = sum(v*arg^(i-1) for (i,v) in enumerate(varlist))
    _plist = map(x->change_base_ring(R, x), plist)
    @debug "Generic polynomial" genpoly varlist

    poly = sum(genpoly(arg+(i - 1)) * p for (i, p) in enumerate(_plist))
    poly -= change_base_ring(R, f)
    @debug "Poly" poly

    coefficients = [Nemo.coeff(poly, i) for i in 0:length(poly)]
    A = [Nemo.coeff(cf, v) for cf in coefficients, v in varlist]
    A = MatrixSpace(base_ring(R), size(A)...)(A)
    rank, ker = kernel(A)
    solutions = typeof(n)[]
    S = base_ring(n)
    for j in 1:rank
        col = Array(ker[:, j])
        s = sum(Nemo.coeff(genpoly, i)(col...)*n^i for i in 0:length(genpoly))
        push!(solutions, s)
    end
    @debug "Algorithm [algpoly] return" solutions
    return solutions
end

function monic_factors(p::PolyElem)
    fac = Nemo.factor(p)
    u = Nemo.coeff(parent(p)(unit(fac)), 0)
    isempty(fac) && return u, [one(p)]
    fs = [(lead(f), div(f, parent(p)(lead(f)))) for (f, m) in fac for i in 1:m]
    return u*prod(map(first, fs)), map(last, fs)
end

function all_monic_factors(p::PolyElem)
    unit, fs = monic_factors(p)
    fs = [one(p); fs]
    unique(prod(s) for s in Combinatorics.powerset(fs, 1))
end

_prod(itr) = reduce(*, itr, init=1)

function alghyper(plist::Vector{T}, n; all_solutions::Bool=false) where {T}
    @debug "Algorithm [alghyper]" plist
    d = length(plist)

    alist = all_monic_factors(plist[1](n))
    blist = all_monic_factors(plist[end](n-(d-2)))

    @debug "All monic factors" alist blist

    solutions = FracElem{typeof(n)}[]
    for a in alist
        for b in blist
            @debug "Monic factors" a b

            ps = [plist[i+1] * _prod(a(n+j) for j in 0:i-1) * _prod(b(n+j) for j in i:d-2) for i in 0:d-1]

            m = maximum(Nemo.degree(p) for p in ps)
            alphas = map(x->Nemo.coeff(x, m), ps)
            zpoly = sum(c*n^(i-1) for (i, c) in enumerate(alphas))
            zs = roots(zpoly)
            zs = [z for z in zs if !iszero(z)]
            @debug "Nonzero roots" alphas zpoly zs
            for z in zs
                _ps = [z^(i-1)*p for (i, p) in enumerate(ps)]
                @debug "Coefficients auxiliary recurrence" _ps

                _sols = algpoly(_ps, zero(n), n)
                @debug "Solutions of algpoly" _sols

                for c in _sols
                    s = z * (a(n) // b(n)) * (c(n+1) // c(n))
                    !all_solutions && return s
                    push!(solutions, s)
                end
            end
        end
    end
    @debug "Algorithm [alghyper] return" unique(solutions) solutions
    !all_solutions && return nothing
    return unique(solutions)
end

function gosper(t::Seq{T}, arg) where {T <: FieldElem}
    @debug "[Gosper]" t arg
    r = div(t(arg+1), t(arg))
    @assert isrational(r)
    r = first(coeffs(r))
    c1, fac1 = monic_factors(numerator(r))
    c2, fac2 = monic_factors(denominator(r))
    if isempty(fac1) && isempty(fac2)
        @debug "Rational coefficients" c1 c2
        return parent(t)(c1//c2*arg)
    end
    R = parent(arg)
    f = R(_prod(fac1))
    g = R(_prod(fac2))
    @debug "Gosper start" t r f g
    @assert gcd(f, g) == 1

    _R, _h = PolynomialRing(base_ring(R), "h")
    _arg = change_base_ring(_R, arg)
    _res = resultant(f(_arg), g(_arg+_h))
    @debug "Resultant" _res

    hs = roots(_res)
    hs = [numerator(r) for r in hs if r >= 0 && isone(denominator(r))]
    @debug "Non-negative integer roots of resultant" hs

    sz, type = length(hs), typeof(arg)
    ss = Array{type}(undef, sz)
    ps = Array{type}(undef, sz+1)
    qs = Array{type}(undef, sz+1)
    ps[1] = f
    qs[1] = g
    for j in 1:length(hs)
        ss[j] = gcd(ps[j](arg), qs[j](arg+hs[j]))
        ps[j+1] = div(ps[j], ss[j])
        qs[j+1] = div(qs[j], ss[j](arg-hs[j]))
    end
    a = c1//c2 * ps[end]
    b = qs[end]
    c = R(_prod(ss[i](arg-j) for i in 1:length(hs) for j in 1:Int(hs[i])))
    @debug "Coefficients for auxiliary equation" a b c ss

    if degree(a) != degree(b) || lead(a) != lead(b)
        D = [degree(c) - max(degree(a), degree(b))]
    else
        d = degree(a)
        if iszero(d)
            D = [degree(c) - d + 1]
        else
            A = Nemo.coeff(a, d-1)
            B = Nemo.coeff(b(arg-1), d-1)
            dd = (B-A)/lead(a)
            if dd >= 0 && isone(denominator(dd))
                D = [degree(c) - d + 1, convert(Int, numerator(dd))]
            else
                D = [degree(c) - d + 1]
            end
        end
    end
    filter!(x->x>=0, D)
    @debug "Degree set" D
    isempty(D) && return nothing
    d = maximum(D)

    # Use method of undetermined coefficients
    _S, _vs = PolynomialRing(base_ring(R), ["v$i" for i = 0:d])
    _arg = change_base_ring(_S, arg)
    _a = change_base_ring(_S, a)
    _b = change_base_ring(_S, b)
    aux = sum(_vs[i+1]*_arg^i for i = 0:d)
    poly = _a(_arg)*aux(_arg+1) - _b(_arg-1)*aux(_arg)
    @debug "Solution template" aux poly

    deg = degree(poly)
    cs = [Nemo.coeff(poly, i) for i in 0:deg]
    mA = [Nemo.coeff(c, v) for c in cs, v in _vs]
    mB = [Nemo.coeff(c, i) for i in 0:deg]
    @debug "Solve" mA mB

    mM = hcat(mA, mB)
    mM = MatrixSpace(base_ring(arg), size(mM)...)(mM)
    r, mM = rref(mM)

    xs = zeros(base_ring(arg), d+1)
    nc = ncols(mM)
    for i in 1:r
        I = findfirst(isone, Array(mM[i, 1:nc-1]))
        xs[I[2]] = mM[i, nc]
    end
    @debug "Method of undetermined coefficients" r mM xs
    iszero(xs) && return nothing

    x = sum(xs[i+1]*arg^i for i = 0:d)
    z = parent(t)(b(arg-1)*x(arg)//c(arg))*t
    @debug "Solution from gosper" z x b c
    @assert z(arg+1) - z(arg) == t
    return z
end

issolution(s::Seq, cfs::Vector{<:FracElem}, arg::PolyElem) = iszero(sum(c * s(arg + i - 1) for (i, c) in enumerate(cfs)))

function petkovsek(cfs::Vector{<:FracElem}, arg::PolyElem, level::Int=0)
    @debug "[Petkovsek]" cfs arg
    # Get rid of denominator
    m = reduce(lcm, map(denominator, cfs))
    polys = [numerator(c * m) for c in cfs]

    ls = map(monic_factors, polys)
    cs = map(first, ls)
    g = reduce(gcd, cs)
    cs = map(x->div(x, g), cs)
    polys = [parent(arg)(c*p) for (c, p) in zip(cs, map(_prod∘last, ls))]

    g = reduce(gcd, polys)
    polys = map(x->div(x, g), polys)
    S = SequenceRing(parent(arg))
    sz = length(polys)
    if sz == 2
        # Recurrence is hypergeometric
        return [S((-polys[1]//polys[2])(arg), arg+1)(arg-1)]
    end

    # Get first hypergeometric solution
    hgterm = alghyper(polys, arg)
    if isnothing(hgterm)
        @debug "No solution from Hyper"
        return nothing
    end
    s_n = S(hgterm(arg), arg+1)(arg-1)
    @debug "Hyper" hgterm s_n
    @assert issolution(s_n, cfs, arg)

    # Take y(n) = z(n)s(n) for unknown z(n), and plug into original recurrence.
    # Then reduce order by taking u(n) = z(n+1) - z(n)
    R, zs = PolynomialRing(S, ["z$i" for i in 0:sz-1])
    poly = sum(s_n(arg+(i-1))*S(cfs[i])*zs[i] for i in 1:sz)

    bs = [Nemo.coeff(poly, z) for z in zs]
    for i in reverse(2:length(bs)-1)
        bs[i] = bs[i] + bs[i+1]
    end
    bs = bs[2:end]
    g = reduce(gcd, bs)
    bs = map(x->div(x, g), bs)
    @debug "Coefficients for order-reduced recurrence" g bs map(isrational, bs) terms(bs[end])
    @assert all(isrational(b) for b in bs)
    # Recursively solve for u(n)
    u_basis = petkovsek(map(first∘coeffs, bs), arg, level+1)
    @assert all(iszero(sum(c * u(arg + i - 1) for (i, c) in enumerate(bs))) for u in u_basis)
    # Use Gosper to resolve antidifference, i.e. find z(n) satisfying u(n) = z(n+1) - z(n)
    z_basis = [gosper(u, arg) for u in u_basis]
    if any(map(isnothing, z_basis))
        @debug "No solution from Gosper"
        return nothing
    end
    @debug "Solutions from Gosper" z_basis
    ss = [s_n * z for z in z_basis]
    @assert all(issolution(s, cfs, arg) for s in ss)
    [s_n; ss]
end


text:
# returns solution of part 1
function x()
	f = open("14.txt")
	andmask = 2^36-1
	ormask = 0
	arr = Dict{Int64, Int64}()
	for i = 1:549
		l = readline(f)
		if l[1:6] == "mask ="
			msk = l[8:end]
			andmask = parse(Int, replace(msk, "X" => "1"), base=2)
			ormask = parse(Int, replace(msk, "X" => "0"), base=2)
		elseif l[1:3] == "mem"
			i = findfirst("=", l).start
			ind = parse(Int, l[5:i-3])
			val = parse(Int, l[i+2:end])
			val &= andmask
			val |= ormask
			arr[ind] = val
		end
	end
	s = 0
	for (k,v) in arr
		s += v
	end
	s
end

# returns solution of part 2
function y()
	f = open("14.txt")
	ormasks = []
	andmask = 2^36-1
	arr = Dict{Int64, Int64}()
	for i = 1:549
		l = readline(f)
		if l[1:6] == "mask ="
			msk = l[8:end]
			xs = count(c->c=='X',msk)
			ormasks = []
			for i=0:(1<<xs-1)
				mask = collect(msk)
				andmask2 = collect(repeat('1', 36))
				j = 0
				for k = length(mask):-1:1
					if mask[k] == 'X'
						if i & (1 << j) > 0
							mask[k] = '1'
						else
							mask[k] = '0'
						end
						andmask2[k] = '0'
						j += 1
					end
				end
				push!(ormasks, parse(Int, String(mask), base=2))
				andmask = parse(Int, String(andmask2), base=2)
			end
		elseif l[1:3] == "mem"
			i = findfirst("=", l).start
			ind = parse(Int, l[5:i-3])
			val = parse(Int, l[i+2:end])
			for mask in ormasks
				ind2 = ind & andmask
				ind3 = ind2 | mask
				arr[ind3] = val
			end
		end
	end
	s = 0
	for (k,v) in arr
		s += v
	end
	s
end

ext: jl


text:
@testset "Bitmask flags" begin
    test(BitmaskDefinition, bitmask_by_name, :VkQueryPipelineStatisticFlagBits, :(
        @bitmask_flag QueryPipelineStatisticFlag::UInt32 begin
            QUERY_PIPELINE_STATISTIC_INPUT_ASSEMBLY_VERTICES_BIT = 1
            QUERY_PIPELINE_STATISTIC_INPUT_ASSEMBLY_PRIMITIVES_BIT = 2
            QUERY_PIPELINE_STATISTIC_VERTEX_SHADER_INVOCATIONS_BIT = 4
            QUERY_PIPELINE_STATISTIC_GEOMETRY_SHADER_INVOCATIONS_BIT = 8
            QUERY_PIPELINE_STATISTIC_GEOMETRY_SHADER_PRIMITIVES_BIT = 16
            QUERY_PIPELINE_STATISTIC_CLIPPING_INVOCATIONS_BIT = 32
            QUERY_PIPELINE_STATISTIC_CLIPPING_PRIMITIVES_BIT = 64
            QUERY_PIPELINE_STATISTIC_FRAGMENT_SHADER_INVOCATIONS_BIT = 128
            QUERY_PIPELINE_STATISTIC_TESSELLATION_CONTROL_SHADER_PATCHES_BIT = 256
            QUERY_PIPELINE_STATISTIC_TESSELLATION_EVALUATION_SHADER_INVOCATIONS_BIT = 512
            QUERY_PIPELINE_STATISTIC_COMPUTE_SHADER_INVOCATIONS_BIT = 1024
        end
    ))

    test(BitmaskDefinition, bitmask_by_name, :VkSparseMemoryBindFlagBits, :(
        @bitmask_flag SparseMemoryBindFlag::UInt32 begin
            SPARSE_MEMORY_BIND_METADATA_BIT = 1
        end
    ))
end


text:
#=
Heru Handika
18 September 2020
Euler's number implementation in julia
=#

using Printf

function calculate_factorial(n)
    if n == 0 || n == 1
        return 1
    else
        return n * calculate_factorial(n - 1)
    end
end

function main()
    euler_number::Float64 = 1.0
    for i = 1:10
        euler_number = euler_number + 1.0 / calculate_factorial(i)
    end
    @printf("Euler number: %.8f\n", euler_number)

end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

max_stars_repo_name: hhandika/stats-from-scratch

max_issues_repo_name: hhandika/stats-from-scratch

max_forks_repo_name: hhandika/stats-from-scratch


text:
# Structs
""" Struct containing simulation settings."""
struct Settings
    coupled::Bool  # choose coupled or uncoupled Norton Equivalent model
    harmonics::Array  # only uneven harmonics, must contain "1", example: [1,3,5,7]
    base_frequency::Int  # in Hz, default: 50
    base_voltage::Number  # in V, default: 230
    base_power::Number  # in W, default: 1000

    # NR algorithm parameters
    thresh_f::Number  # default: 1e-6
    thresh_h::Number  # default: 1e-4
    max_iter_f::Int  # default: 30
    max_iter_h::Int  # default: 50

    # voltage start values [pu]
    v_f::Number  # default: 1
    ϕ_f::Number  # default: 0
    v_h::Number  # default: 0.1
    ϕ_h::Number  # default: 0

    # calculated parameters
    K::Int  # number of harmonics (without fundamental)
    K1::Int # number of harmonics (including fundamental)
    base_current::Number
    base_admittance::Number
    base_impedance::Number
end

# Functions
"""
    init_settings(coupled, harmonics; ...)

Initialize settings struct.

Default values:
- base_frequency=50
- base_voltage=230
- base_power=1000
- thresh_f = 1e-6
- max_iter_f = 30
- thresh_h=1e-4
- max_iter_h=50
- v_f=1
- ϕ_f = 0
- v_h=0.1
- ϕ_h=0
"""
function init_settings(coupled, harmonics;
    base_frequency=50, base_voltage=230, base_power=1000,
    thresh_f = 1e-6, max_iter_f = 30, thresh_h=1e-4, max_iter_h=50,
    v_f=1, ϕ_f = 0, v_h=0.1, ϕ_h=0)

    K = length(harmonics) - 1
    K1 = length(harmonics)

    # derived pu base values
    base_current = base_power/base_voltage
    base_admittance = base_current/base_voltage
    base_impedance = 1/base_admittance

    Settings(
        coupled, harmonics, base_frequency, base_voltage, base_power,
        thresh_f, thresh_h, max_iter_f, max_iter_h,
        v_f, ϕ_f, v_h, ϕ_h,
        K, K1, base_current, base_admittance, base_impedance)
end


text:
module DOT

using GraphIO.ParserCombinator.Parsers
using LightGraphs
using LightGraphs: AbstractGraphFormat

import LightGraphs: loadgraph, loadgraphs

export DOTFormat

struct DOTFormat <: AbstractGraphFormat end
# TODO: implement save

function _dot_read_one_graph(pg::Parsers.DOT.Graph)
    isdir = pg.directed
    nvg = length(Parsers.DOT.nodes(pg))
    nodedict = Dict(zip(collect(Parsers.DOT.nodes(pg)), 1:nvg))
    if isdir
        g = LightGraphs.DiGraph(nvg)
    else
        g = LightGraphs.Graph(nvg)
    end
    for es in Parsers.DOT.edges(pg)
        s = nodedict[es[1]]
        d = nodedict[es[2]]
        add_edge!(g, s, d)
    end
    return g
end

function loaddot(io::IO, gname::String)
    p = Parsers.DOT.parse_dot(read(io, String))
    for pg in p
        isdir = pg.directed
        possname = isdir ? Parsers.DOT.StringID("digraph") : Parsers.DOT.StringID("graph")
        name = get(pg.id, possname).id
        name == gname && return _dot_read_one_graph(pg)
    end
    error("Graph $gname not found")
end

function loaddot_mult(io::IO)
    p = Parsers.DOT.parse_dot(read(io, String))

    graphs = Dict{String,AbstractGraph}()

    for pg in p
        isdir = pg.directed
        possname = isdir ? Parsers.DOT.StringID("digraph") : Parsers.DOT.StringID("graph")
        name = get(pg.id, possname).id
        graphs[name] = _dot_read_one_graph(pg)
    end
    return graphs
end

loadgraph(io::IO, gname::String, ::DOTFormat) = loaddot(io, gname)
loadgraphs(io::IO, ::DOTFormat) = loaddot_mult(io)

end #module


text:
### FieldTuple types

# FieldTuple is a thin wrapper around a Tuple or NamedTuple holding some Fields
# and behaving like a Field itself
struct FieldTuple{FS<:Union{Tuple,NamedTuple},B,T} <: Field{B,T}
    fs :: FS
    function FieldTuple(fs::FS) where {FS<:Union{Tuple,NamedTuple}}
        B = BasisProd{Tuple{map(basis,values(fs))...}}
        T = promote_type(map(eltype,values(fs))...)
        new{FS,B,T}(fs)
    end
end
# FieldTuple(args...) or FieldTuple(;kwargs...) calls the inner constructor
# which takes a single Tuple/NamedTuple:
(::Type{FT})(;kwargs...) where {FT<:FieldTuple} = FT((;kwargs...))
(::Type{FT})(f1::Field,f2::Field,fs::Field...) where {FT<:FieldTuple} = FT((f1,f2,fs...))

### printing
getindex(f::FieldTuple,::Colon) = vcat(getindex.(values(f.fs),:)...)[:]
getindex(D::DiagOp{<:FieldTuple}, i::Int, j::Int) = (i==j) ? D.diag[:][i] : diagzero(D, i, j)
typealias_def(::Type{<:FieldTuple{NamedTuple{Names,FS},T}}) where {Names,FS<:Tuple,T} =
    "Field-($(join(map(string,Names),",")))-$FS"
typealias_def(::Type{<:FieldTuple{FS,T}}) where {FS<:Tuple,T} =
    "Field-$(tuple_type_len(FS))-$FS"

### array interface
size(f::FieldTuple) = (mapreduce(length, +, f.fs, init=0),)
copy(f::FieldTuple) = FieldTuple(map(copy,f.fs))
iterate(ft::FieldTuple, args...) = iterate(ft.fs, args...)
getindex(f::FieldTuple, i::Union{Int,UnitRange}) = getindex(f.fs, i)
fill!(ft::FieldTuple, x) = (map(f->fill!(f,x), ft.fs); ft)
adapt_structure(to, f::FieldTuple) = FieldTuple(map(f->adapt(to,f),f.fs))
similar(ft::FieldTuple) = FieldTuple(map(similar,ft.fs))
similar(ft::FieldTuple, ::Type{T}) where {T<:Number} = FieldTuple(map(f->similar(f,T),ft.fs))
similar(ft::FieldTuple, Nbatch::Int) = FieldTuple(map(f->similar(f,Nbatch),ft.fs))
sum(f::FieldTuple; dims=:) = dims == (:) ? sum(sum, f.fs) : error("sum(::FieldTuple, dims=$dims not supported")

### broadcasting
# see base_fields.jl for more explanation of all these pieces, its the
# exact same principle
struct FieldTupleStyle{S,Names} <: AbstractArrayStyle{1} end
function BroadcastStyle(::Type{<:FieldTuple{TS}}) where {TS<:Tuple}
    FieldTupleStyle{Tuple{map_tupleargs(typeof∘BroadcastStyle,TS)...}, Nothing}()
end
function BroadcastStyle(::Type{<:FieldTuple{NamedTuple{Names,TS}}}) where {Names,TS<:Tuple}
    FieldTupleStyle{Tuple{map_tupleargs(typeof∘BroadcastStyle,TS)...}, Names}()
end
function BroadcastStyle(::FieldTupleStyle{S₁,Names}, ::FieldTupleStyle{S₂,Names}) where {S₁,S₂,Names}
    FieldTupleStyle{Tuple{map_tupleargs((s₁,s₂)->typeof(result_style(s₁(),s₂())), S₁, S₂)...}, Names}()
end
BroadcastStyle(S::FieldTupleStyle, ::DefaultArrayStyle{0}) = S

@generated function materialize(bc::Broadcasted{FieldTupleStyle{S,Names}}) where {S,Names}
    wrapper = Names == Nothing ? :tuple : :(NamedTuple{$Names})
    exprs = map_tupleargs(S, tuple(1:tuple_type_len(S)...)) do Sᵢ, i
        :(materialize(convert(Broadcasted{$Sᵢ}, preprocess(($(S.parameters[i])(),FieldTupleComponent{$i}()), bc))))
    end
    :(FieldTuple($wrapper($(exprs...))))
end
@generated function materialize!(dst::FieldTuple, bc::Broadcasted{FieldTupleStyle{S,Names}}) where {S,Names}
    exprs = map_tupleargs(S, tuple(1:tuple_type_len(S)...)) do Sᵢ, i
        :(materialize!(dst.fs[$i], convert(Broadcasted{$Sᵢ}, preprocess(($(S.parameters[i])(),FieldTupleComponent{$i}()), bc))))
    end
    :(begin $(exprs...) end; dst)
end

struct FieldTupleComponent{i} end

preprocess(::Tuple{<:Any,FieldTupleComponent{i}}, ft::FieldTuple) where {i} = ft.fs[i]

preprocess(dest::Tuple{FieldTupleStyle{S},<:Any}, bc::Broadcasted) where {S} =
    broadcasted(S(), bc.f, preprocess_args(dest, bc.args)...)

### promotion
function promote(ft1::FieldTuple, ft2::FieldTuple)
    fts = map(promote, ft1.fs, ft2.fs)
    FieldTuple(map(first,fts)), FieldTuple(map(last,fts))
end

### conversion
Basis(ft::FieldTuple) = ft
(::Type{B})(ft::FieldTuple) where {B<:Basis}     = FieldTuple(map(B, ft.fs))
(::Type{B})(ft::FieldTuple) where {B<:Basislike} = FieldTuple(map(B, ft.fs))
(::Type{B})(ft::FieldTuple{<:Tuple}) where {Bs,B<:BasisProd{Bs}} =
    FieldTuple(map_tupleargs((B,f)->B(f), Bs, ft.fs))
(::Type{B})(ft::FieldTuple{<:NamedTuple{Names}}) where {Names,Bs,B<:BasisProd{Bs}} =
    FieldTuple(NamedTuple{Names}(map_tupleargs((B,f)->B(f), Bs, values(ft.fs))))

### properties
getproperty(f::FieldTuple, s::Symbol) = getproperty(f, Val(s))
getproperty(f::FieldTuple, ::Val{:fs}) = getfield(f, :fs)
getproperty(f::FieldTuple, ::Val{s}) where {s} = getproperty(getfield(f,:fs), s)
propertynames(f::FieldTuple) = (:fs, propertynames(f.fs)...)

### simulation
white_noise(ξ::FieldTuple, rng::AbstractRNG) = FieldTuple(map(f -> white_noise(f, rng), ξ.fs))

### Diagonal-ops
# need a method specific for FieldTuple since we don't carry around
# the basis in a way that works with the default implementation
(*)(D::DiagOp{<:FieldTuple}, f::FieldTuple) = FieldTuple(map((d,f)->Diagonal(d)*f, D.diag.fs, f.fs))
(\)(D::DiagOp{<:FieldTuple}, f::FieldTuple) = FieldTuple(map((d,f)->Diagonal(d)\f, D.diag.fs, f.fs))

# # promote before recursing for these
dot(a::FieldTuple, b::FieldTuple) = reduce(+, map(dot, getfield.(promote(a,b),:fs)...), init=0)
hash(ft::FieldTuple, h::UInt64) = foldr(hash, (typeof(ft), ft.fs), init=h)

# logdet & trace
logdet(L::Diagonal{<:Union{Real,Complex}, <:FieldTuple}) = reduce(+, map(logdet∘Diagonal, L.diag.fs), init=0)
tr(L::Diagonal{<:Union{Real,Complex}, <:FieldTuple}) = reduce(+, map(tr∘Diagonal, L.diag.fs), init=0)

# misc
batch_length(ft::FieldTuple) = only(unique(map(batch_length, ft.fs)))
batch_index(ft::FieldTuple, I) = FieldTuple(map(f -> batch_index(f, I), ft.fs))
function global_rng_for(::Type{<:FieldTuple{<:Union{FS,NamedTuple{Names,FS}}}}) where {Names,FS<:Tuple}
    only(unique(map_tupleargs(global_rng_for, FS)))
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 ApplicationSecurityGroupPropertiesFormat <: SwaggerModel
    resourceGuid::Any # spec type: Union{ Nothing, String } # spec name: resourceGuid
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function ApplicationSecurityGroupPropertiesFormat(;resourceGuid=nothing, provisioningState=nothing)
        o = new()
        validate_property(ApplicationSecurityGroupPropertiesFormat, Symbol("resourceGuid"), resourceGuid)
        setfield!(o, Symbol("resourceGuid"), resourceGuid)
        validate_property(ApplicationSecurityGroupPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ApplicationSecurityGroupPropertiesFormat

const _property_map_ApplicationSecurityGroupPropertiesFormat = Dict{Symbol,Symbol}(Symbol("resourceGuid")=>Symbol("resourceGuid"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ApplicationSecurityGroupPropertiesFormat = Dict{Symbol,String}(Symbol("resourceGuid")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ ApplicationSecurityGroupPropertiesFormat }) = collect(keys(_property_map_ApplicationSecurityGroupPropertiesFormat))
Swagger.property_type(::Type{ ApplicationSecurityGroupPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationSecurityGroupPropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationSecurityGroupPropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationSecurityGroupPropertiesFormat[property_name]

function check_required(o::ApplicationSecurityGroupPropertiesFormat)
    true
end

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


text:
"""
Ports

Utility structure

Information on scanned ports (<ports> and <extraports>)
"""
Base.@kwdef mutable struct Ports <:Marsh
ports::Vector{Port}             = Vector{Port}()
extraports::Vector{ExtraPort}   = Vector{ExtraPort}()
end
StructTypes.names(::Type{Ports}) = ((:ports, :port),)

# Ease iteration over Ports
function Base.iterate(ps::Ports, state=1)
if state + 1 > length(ps.ports)
    return nothing
end
extra = state+1 > length(ps.extraports) ? nothing : ps.extraports[state+1]
return ([ps.ports[state+1], extra], state+1)
end


text:
module IbanGen

include("common.jl")
include("check_digits.jl")
include("bban.jl")

export is_supported_country, supported_countries
export iban, iban_random, build
export ValidationException

"""
    iban_random(
        CountryCode::Maybe{String}=nothing,
        BankCode::Maybe{String}=nothing,
        AccountNumber::Maybe{String}=nothing,
        BranchCode::Maybe{String}=nothing,
        NationalCheckDigit::Maybe{String}=nothing,
        AccountType::Maybe{String}=nothing,
        OwnerAccountType::Maybe{String}=nothing,
        IdentificationNumber::Maybe{String}=nothing
    )::Dict{String,String}

Generate a random IBAN subject to the provided attributes. For an attributes that is not provided, a
random value will be used according to the rules of the (provided or generated) country. Attributes
that are not defined for a country are ignored.

## Example
```julia
julia> iban_random()
Dict{String,String} with 6 entries:
  "CountryCode"   => "GR"
  "BranchCode"    => "7500"
  "AccountNumber" => "1HRB7OApF5ABTOYH"
  "value"         => "GR8410975001HRB7OApF5ABTOYH"
  "BankCode"      => "109"
  "CheckDigits"   => "84"

julia> iban_random(CountryCode = "GR", BankCode = "109")
Dict{String,String} with 6 entries:
  "CountryCode"   => "GR"
  "BranchCode"    => "2170"
  "AccountNumber" => "24wO2qBgz1ROP82L"
  "value"         => "GR26109217024wO2qBgz1ROP82L"
  "BankCode"      => "109"
  "CheckDigits"   => "26"

```
"""
function iban_random(;
    CountryCode::MaybeString=nothing,
    BankCode::MaybeString=nothing,
    AccountNumber::MaybeString=nothing,
    BranchCode::MaybeString=nothing,
    NationalCheckDigit::MaybeString=nothing,
    AccountType::MaybeString=nothing,
    OwnerAccountType::MaybeString=nothing,
    IdentificationNumber::MaybeString=nothing
)::Dict{String,String}

    partup = (;
        Symbol("IbanGen.BankCode") => BankCode,
        Symbol("IbanGen.AccountNumber") => AccountNumber,
        Symbol("IbanGen.BranchCode") => BranchCode,
        Symbol("IbanGen.NationalCheckDigit") => NationalCheckDigit,
        Symbol("IbanGen.AccountType") => AccountType,
        Symbol("IbanGen.OwnerAccountType") => OwnerAccountType,
        Symbol("IbanGen.IdentificationNumber") => IdentificationNumber
    )

    country_code = CountryCode === nothing ? rand(supported_countries()) : CountryCode

    theiban = parse_bban!(country_code, partup, true)

    theiban.check_digits = calculate_check_digit(
        theiban.bban_str,
        theiban.country_code
    )

    theiban.iban_str =
        theiban.country_code * theiban.check_digits * theiban.bban_str

    _as_dict(theiban)
end

"""
    iban(
        CountryCode::String,
        BankCode::String,
        AccountNumber::String,
        BranchCode::Maybe{String}=nothing,
        NationalCheckDigit::Maybe{String}=nothing,
        AccountType::Maybe{String}=nothing,
        OwnerAccountType::Maybe{String}=nothing,
        IdentificationNumber::Maybe{String}=nothing,
    )::Dict{String,String}

 Generate an IBAN based on the provided parameters.

## Example
```julia
julia> iban(CountryCode = "GB", BankCode = "NWBK", BranchCode = "601613",AccountNumber = "31926819")
Dict{String,String} with 6 entries:
  "CountryCode"   => "GB"
  "BranchCode"    => "601613"
  "AccountNumber" => "31926819"
  "value"         => "GB29NWBK60161331926819"
  "BankCode"      => "NWBK"
  "CheckDigits"   => "29"

```
"""
function iban(;
    CountryCode::String,
    BankCode::String,
    AccountNumber::String,
    BranchCode::MaybeString=nothing,
    NationalCheckDigit::MaybeString=nothing,
    AccountType::MaybeString=nothing,
    OwnerAccountType::MaybeString=nothing,
    IdentificationNumber::MaybeString=nothing
)::Dict{String,String}

    partup = (;
        Symbol("IbanGen.BankCode") => BankCode,
        Symbol("IbanGen.AccountNumber") => AccountNumber,
        Symbol("IbanGen.BranchCode") => BranchCode,
        Symbol("IbanGen.NationalCheckDigit") => NationalCheckDigit,
        Symbol("IbanGen.AccountType") => AccountType,
        Symbol("IbanGen.OwnerAccountType") => OwnerAccountType,
        Symbol("IbanGen.IdentificationNumber") => IdentificationNumber
    )

    theiban = parse_bban!(CountryCode, partup)

    theiban.check_digits = calculate_check_digit(
        theiban.bban_str,
        theiban.country_code
    )

    theiban.iban_str =
        theiban.country_code * theiban.check_digits * theiban.bban_str

    _as_dict(theiban)
end

"""
    iban(iban_str::String)::Dict{String,String}

Validated the provided string and parse it into an IBAN structure.

## Example
```julia
julia> iban("BR9700360305000010009795493P1")
Dict{String,String} with 8 entries:
  "CountryCode"      => "BR"
  "CheckDigits"      => "97"
  "BranchCode"       => "00001"
  "AccountType"      => "P"
  "AccountNumber"    => "0009795493"
  "value"            => "BR9700360305000010009795493P1"
  "BankCode"         => "00360305"
  "OwnerAccountType" => "1"

```
"""
function iban(iban_str::String)::Dict{String,String}

    ensure(
        length(iban_str) >= 4,
        iban_str, "value too short"
    )

    country_code = SubString(iban_str, 1, 2)
    bban_structure = for_country(country_code)
    ensure(
        bban_structure !== nothing,
        country_code, "country code not supported"
    )

    bban_str = SubString(iban_str, 5)
    ensure(
        _length(bban_structure) == length(bban_str),
        bban_str, "unexpected BBAN length, expected: $(_length(bban_structure))"
    )

    dict = Dict(slice[1] => bban_str[slice[2]] for slice in _slicing(bban_structure))
    theiban = parse_bban!(country_code, dict)

    provided = SubString(iban_str, 3, 4)
    expected = calculate_check_digit(bban_str, country_code)
    ensure(
        expected == provided,
        provided, "invalid check digits, expected: $(expected)"
    )

    theiban.check_digits = expected
    theiban.iban_str = iban_str

    _as_dict(theiban)
end

struct IbanEntry
    value::AbstractString
    entry_type::EntryType
    function IbanEntry(value, entry_type)
        ensure(
            !isnothing(value),
            value, "Value not provided [$(typeof(entry_type))]"
        )
        ensure(
            length(value) == entry_type.spec.length,
            value, "Unexpected length $(length(value)), expected: $(entry_type.spec.length) [$(typeof(entry_type))]"
        )
        ensure(
            issubset(value, entry_type.spec.alphabet),
            value, "invalid characters [$(typeof(entry_type))]"
        )
        new(value, entry_type)
    end
end

mutable struct TheIban
    country_code::String
    check_digits::String
    bban::Vector{IbanEntry}
    bban_str::MaybeString
    iban_str::MaybeString

    function TheIban(country_code)
        ensure(
            is_supported_country(country_code),
            country_code, "country code not supported"
        )
        new(country_code, "00", [], nothing, nothing)
    end
end

#=
#############################################################
Converts the internal IBAN representation to a `Dict`
#############################################################
=#
function _as_dict(theiban::TheIban)::Dict{String,String}
    result = Dict(stringify(entry) => entry.value for entry in values(theiban.bban))
    result["CountryCode"] = theiban.country_code
    result["CheckDigits"] = theiban.check_digits
    result["value"] = theiban.iban_str
    result
end

end


text:
using Printf: @sprintf

function StatsBase.dof(m::GeneralizedLinearMixedModel)
    length(m.β) + length(m.θ) + GLM.dispersion_parameter(m.resp.d)
end

StatsBase.fitted(m::GeneralizedLinearMixedModel) = m.resp.mu

function fixef(m::GeneralizedLinearMixedModel{T}, permuted=true) where {T}
    permuted && return m.β
    Xtrm = m.LMM.trms[end - 1]
    iperm = invperm(Xtrm.piv)
    p = length(iperm)
    r = Xtrm.rank
    r == p ? m.β[iperm] : copyto!(fill(-zero(T), p), m.β)[iperm]
end

function GeneralizedLinearMixedModel(f::Formula, fr::AbstractDataFrame,
         d::Distribution, l::Link; wt=[], offset=[], contrasts = Dict())
    if d == Binomial() && isempty(wt)
        d = Bernoulli()
    end
    LMM = LinearMixedModel(f, fr; weights = wt, contrasts=contrasts, rdist=d)
    X = LMM.trms[end - 1].x
    T = eltype(X)
    y = copy(model_response(LMM))
    if isempty(wt)
        LMM = LinearMixedModel(LMM.formula, LMM.trms, fill!(similar(y), 1), LMM.A, LMM.L, LMM.optsum)
    end
    updateL!(setθ!(LMM, getθ(LMM)))
            # fit a glm to the fixed-effects only - awkward syntax is to by-pass a test
    gl = isempty(wt) ? glm(X, y, d, l) : glm(X, y, d, l, wts=wt)
    β = coef(gl)
    u = [fill(zero(T), vsize(t), nlevs(t)) for t in reterms(LMM)]
    vv = length(u) == 1 ? vec(u[1]) : T[]

    res = GeneralizedLinearMixedModel(LMM, β, copy(β), getθ(LMM), copy.(u), u,
        zero.(u), gl.rr, similar(y), oftype(y, wt), similar(vv),
        similar(vv), similar(vv), similar(vv))
    setβθ!(res, vcat(coef(gl), getθ(LMM)))
    deviance!(res, 1)
    res
end

GeneralizedLinearMixedModel(f::Formula, fr::AbstractDataFrame, d::Distribution;
    wt=[], offset=[], contrasts=Dict()) = GeneralizedLinearMixedModel(f, fr, d,
    GLM.canonicallink(d), wt=wt, offset=offset, contrasts=contrasts)

StatsBase.fit(::Type{GeneralizedLinearMixedModel},
              f::Formula,
              fr::AbstractDataFrame,
              d::Distribution,
              l::Link = GLM.canonicallink(d);
              wt=[],
              offset=[],
              contrasts=Dict(),
              verbose=false,
              fast=false,
              nAGQ=1) =
    fit!(GeneralizedLinearMixedModel(f, fr, d, l, wt=wt, offset=offset, contrasts=contrasts),
         verbose=verbose, fast=fast, nAGQ=nAGQ)

"""
    deviance(m::GeneralizedLinearMixedModel{T}, nAGQ=1)::T where {T}

Return the deviance of `m` evaluated by adaptive Gauss-Hermite quadrature

If the distribution `D` does not have a scale parameter the Laplace approximation
is defined as the squared length of the conditional modes, `u`, plus the determinant
of `Λ'Z'WZΛ + I`, plus the sum of the squared deviance residuals.
"""
function StatsBase.deviance(m::GeneralizedLinearMixedModel{T}, nAGQ=1) where {T}
    nAGQ == 1 && return T(sum(m.resp.devresid) + logdet(m) + sum(u -> sum(abs2, u), m.u))
    u = vec(m.u[1])
    u₀ = vec(m.u₀[1])
    copyto!(u₀, u)
    ra = RaggedArray(m.resp.devresid, m.LMM.trms[1].refs)
    devc0 = sum!(map!(abs2, m.devc0, u), ra)  # the deviance components at z = 0
    sd = map!(inv, m.sd, m.LMM.L.data[Block(1,1)].diag)
    mult = fill!(m.mult, 0)
    devc = m.devc
    for (z, w) in GHnorm(nAGQ)
        if !iszero(w)
            if iszero(z)  # devc == devc0 in this case
                mult .+= w
            else
                @. u = u₀ + z * sd
                updateη!(m)
                sum!(map!(abs2, devc, u), ra)
                @. mult += exp((abs2(z) + devc0 - devc)/2)*w
            end
        end
    end
    copyto!(u, u₀)
    updateη!(m)
    sum(devc0) - 2 * (sum(log, mult) + sum(log, sd))
end

"""
    deviance!(m::GeneralizedLinearMixedModel, nAGQ=1)

Update `m.η`, `m.μ`, etc., install the working response and working weights in
`m.LMM`, update `m.LMM.A` and `m.LMM.R`, then evaluate the [`deviance`](@ref).
"""
function deviance!(m::GeneralizedLinearMixedModel, nAGQ=1)
    updateη!(m)
    GLM.wrkresp!(vec(m.LMM.trms[end].x), m.resp)
    reweight!(m.LMM, m.resp.wrkwt)
    deviance(m, nAGQ)
end

function StatsBase.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T}
    accum = zero(T)
    D = Distribution(m.resp)
    if D <: Binomial
        for (μ, y, n) in zip(m.resp.mu, m.resp.y, m.wt)
            accum += logpdf(D(round(Int, n), μ), round(Int, y * n))
        end
    else
        for (μ, y) in zip(m.resp.mu, m.resp.y)
            accum += logpdf(D(μ), y)
        end
    end
    accum - (mapreduce(u -> sum(abs2, u), + , m.u) + logdet(m)) / 2
end

function lowerbd(m::GeneralizedLinearMixedModel)
    lb = lowerbd(m.LMM)
    vcat(fill(eltype(lb)(-Inf), size(m.β)), lb)
end

StatsBase.nobs(m::GeneralizedLinearMixedModel) = length(m.η)

StatsBase.predict(m::GeneralizedLinearMixedModel) = fitted(m)

updatedevresid!(r::GLM.GlmResp, η::AbstractVector) = updateμ!(r, η)

fastlogitdevres(η, y) = 2log1p(exp(iszero(y) ? η : -η))

function updatedevresid!(r::GLM.GlmResp{V,<:Bernoulli,LogitLink}, η::V) where V<:AbstractVector{<:AbstractFloat}
    map!(fastlogitdevres, r.devresid, η, r.y)
    r
end

"""
    updateη!(m::GeneralizedLinearMixedModel)

Update the linear predictor, `m.η`, from the offset and the `B`-scale random effects.
"""
function updateη!(m::GeneralizedLinearMixedModel)
    η, b, u = m.η, m.b, m.u
    trms = m.LMM.trms
    mul!(η, trms[end - 1].x, m.β)
    for i in eachindex(b)
        unscaledre!(η, trms[i], mul!(b[i], trms[i].Λ, u[i]))
    end
    updateμ!(m.resp, η)
    m
end

average(a::T, b::T) where {T<:AbstractFloat} = (a + b) / 2

"""
    pirls!(m::GeneralizedLinearMixedModel)

Use Penalized Iteratively Reweighted Least Squares (PIRLS) to determine the conditional
modes of the random effects.

When `varyβ` is true both `u` and `β` are optimized with PIRLS.  Otherwise only `u` is
optimized and `β` is held fixed.

Passing `verbose = true` provides verbose output of the iterations.
"""
function pirls!(m::GeneralizedLinearMixedModel{T}, varyβ=false, verbose=false) where {T}
    iter, maxiter, obj = 0, 100, T(-Inf)
    u₀ = m.u₀
    u = m.u
    β = m.β
    β₀ = m.β₀
    lm = m.LMM
    for j in eachindex(u)         # start from u all zeros
        copyto!(u₀[j], fill!(u[j], 0))
    end
    varyβ && copyto!(β₀, β)
    obj₀ = deviance!(m) * 1.0001
    if verbose
        print("varyβ = ", varyβ, ", obj₀ = ", obj₀)
        if varyβ
            print(", β =")
            show(β)
        end
        println()
    end

    while iter < maxiter
        iter += 1
        varyβ && ldiv!(adjoint(feL(m)), copyto!(β, lm.L.data.blocks[end, end - 1]))
        ranef!(u, m.LMM, β, true) # solve for new values of u
        obj = deviance!(m)        # update GLM vecs and evaluate Laplace approx
        verbose && println(lpad(iter,4), ": ", obj)
        nhalf = 0
        while obj > obj₀
            nhalf += 1
            if nhalf > 10
                if iter < 2
                    throw(ErrorException("number of averaging steps > 10"))
                end
                break
            end
            for i in eachindex(u)
                map!(average, u[i], u[i], u₀[i])
            end
            varyβ && map!(average, β, β, β₀)
            obj = deviance!(m)
            verbose && println(lpad(nhalf, 8), ", ", obj)
        end
        if isapprox(obj, obj₀; atol = 0.00001)
            break
        end
        copyto!.(u₀, u)
        copyto!(β₀, β)
        obj₀ = obj
    end
    m
end

"""
    setβθ!(m::GeneralizedLinearMixedModel, v)

Set the parameter vector, `:βθ`, of `m` to `v`.

`βθ` is the concatenation of the fixed-effects, `β`, and the covariance parameter, `θ`.
"""
function setβθ!(m::GeneralizedLinearMixedModel, v)
    setβ!(m, v)
    setθ!(m, view(v, (length(m.β) + 1) : length(v)))
end

function setβ!(m::GeneralizedLinearMixedModel, v)
    β = m.β
    copyto!(β, view(v, 1:length(β)))
    m
end

function setθ!(m::GeneralizedLinearMixedModel, v)
    setθ!(m.LMM, copyto!(m.θ, v))
    m
end

sdest(m::GeneralizedLinearMixedModel{T}) where {T} = convert(T, NaN)

"""
    fit!(m::GeneralizedLinearMixedModel[, verbose = false, fast = false, nAGQ=1])

Optimize the objective function for `m`.

When `fast` is `true` a potentially much faster but slightly less accurate algorithm, in
which `pirls!` optimizes both the random effects and the fixed-effects parameters,
is used.
"""
function StatsBase.fit!(m::GeneralizedLinearMixedModel{T};
                        verbose::Bool=false, fast::Bool=false, nAGQ=1) where {T}
    β = m.β
    lm = m.LMM
    optsum = lm.optsum
    if !fast
        fit!(m, verbose=verbose, fast=true, nAGQ=nAGQ)
        optsum.lowerbd = vcat(fill!(similar(β), T(-Inf)), optsum.lowerbd)
        optsum.initial = vcat(β, m.θ)
        optsum.final = copy(optsum.initial)
        optsum.initial_step = vcat(stderror(m) ./ 3, min.(T(0.05), m.θ ./ 4))
    end
    setpar! = fast ? setθ! : setβθ!
    feval = 0
    function obj(x, g)
        isempty(g) || error("gradient not defined for this model")
        feval += 1
        val = deviance(pirls!(setpar!(m, x), fast, verbose), nAGQ)
        feval == 1 && (optsum.finitial = val)
        if verbose
            println("f", lpad(feval,5), ": ", val, ", ", x)
        end
        val
    end
    opt = Opt(optsum)
    NLopt.min_objective!(opt, obj)
    fmin, xmin, ret = NLopt.optimize(opt, copyto!(optsum.final, optsum.initial))
    ## check if very small parameter values bounded below by zero can be set to zero
    xmin_ = copy(xmin)
    for i in eachindex(xmin_)
        if iszero(optsum.lowerbd[i]) && zero(T) < xmin_[i] < T(0.001)
            xmin_[i] = zero(T)
        end
    end
    if xmin ≠ xmin_
        if (zeroobj = obj(xmin_, T[])) ≤ (fmin + 1.e-5)
            fmin = zeroobj
            copyto!(xmin, xmin_)
        end
    end
    ## ensure that the parameter values saved in m are xmin
    pirls!(setpar!(m, xmin), fast, verbose)
    optsum.nAGQ = nAGQ
    optsum.feval = feval
    optsum.final = xmin
    optsum.fmin = fmin
    optsum.returnvalue = ret
    ret == :ROUNDOFF_LIMITED && @warn("NLopt was roundoff limited")
    if ret ∈ [:FAILURE, :INVALID_ARGS, :OUT_OF_MEMORY, :FORCED_STOP, :MAXEVAL_REACHED]
        @warn("NLopt optimization failure: $ret")
    end
    m
end

function Base.show(io::IO, m::GeneralizedLinearMixedModel)
    nAGQ = m.LMM.optsum.nAGQ
    println(io, "Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = $nAGQ)")
    println(io, "  ", m.LMM.formula)
    println(io, "  Distribution: ", Distribution(m.resp))
    println(io, "  Link: ", Link(m.resp), "\n")
    println(io, string("  Deviance: ", @sprintf("%.4f", deviance(m, nAGQ))), "\n")

    show(io,VarCorr(m))
    gl = grplevels(m.LMM)
    print(io, "\n Number of obs: ", length(m.η), "; levels of grouping factors: ", gl[1])
    for l in gl[2:end]
        print(io, ", ", l)
    end
    println(io)
    println(io, "\nFixed-effects parameters:")
    show(io, coeftable(m))
end

varest(m::GeneralizedLinearMixedModel{T}) where {T} = one(T)


text:
"""
# Description

Rewrite a pair of expressions, `x` and `y`, into a single expression in which
both are the operands of `x || y` in a world in which `||` uses three-valued
logic.

# Arguments

1. `x::Any`: An expression.
2. `y::Any`: An expression.

# Return Values

1. `e′::Any`: A new expression that computes the three-valued logic value of
    `x || y` with proper short-circuiting.

# Examples

```
julia> tvl_or(:a, :b)
quote
    #= REPL[1]:28 =#
    let var"##254" = a
        #= REPL[1]:29 =#
        if var"##254" === true
            #= REPL[1]:30 =#
            true
        elseif #= REPL[1]:31 =# var"##254" === false
            #= REPL[1]:32 =#
            b
        else
            #= REPL[1]:34 =#
            if b === true
                #= REPL[1]:35 =#
                true
            else
                #= REPL[1]:37 =#
                missing
            end
        end
    end
end
```
"""
function tvl_or(@nospecialize(x::Any), @nospecialize(y::Any))
    t1 = gensym()
    quote
        let $t1 = $x
            if $t1 === true
                true
            elseif $t1 === false
                $y
            else
                if $y === true
                    true
                else
                    missing
                end
            end
        end
    end
end

"""
# Description

Rewrite a pair of expressions, `x` and `y`, into a single expression in which
both are the operands of `x && y` in a world in which `&&` uses three-valued
logic.

# Arguments

1. `x::Any`: An expression.
2. `y::Any`: An expression.

# Return Values

1. `e′::Any`: A new expression that computes the three-valued logic value of
    `x && y` with proper short-circuiting.

# Examples

```
julia> tvl_and(:a, :b)
quote
    #= REPL[4]:4 =#
    let var"##255" = a
        #= REPL[4]:5 =#
        if var"##255" === false
            #= REPL[4]:6 =#
            false
        elseif #= REPL[4]:7 =# var"##255" === true
            #= REPL[4]:8 =#
            b
        else
            #= REPL[4]:10 =#
            if b === false
                #= REPL[4]:11 =#
                false
            else
                #= REPL[4]:13 =#
                missing
            end
        end
    end
end
```
"""
function tvl_and(@nospecialize(x), @nospecialize(y))
    t1 = gensym()
    quote
        let $t1 = $x
            if $t1 === false
                false
            elseif $t1 === true
                $y
            else
                if $y === false
                    false
                else
                    missing
                end
            end
        end
    end
end

"""
# Description

Rewrite an arbitary expression, `e`, to use TVL at its top-level without
recursively descending into its components.

# Arguments

1. `e::Any`: An expression.

# Return Values

1. `e′::Any`: A new expression in which `||` or `&&` are replaced with their
    three-valued logic counterparts.

# Examples

```
julia> rewrite_tvl(:(a || b))
quote
    #= REPL[1]:28 =#
    let var"##256" = a
        #= REPL[1]:29 =#
        if var"##256" === true
            #= REPL[1]:30 =#
            true
        elseif #= REPL[1]:31 =# var"##256" === false
            #= REPL[1]:32 =#
            b
        else
            #= REPL[1]:34 =#
            if b === true
                #= REPL[1]:35 =#
                true
            else
                #= REPL[1]:37 =#
                missing
            end
        end
    end
end
```
"""
function rewrite_tvl(@nospecialize(e::Any))
    if isa(e, Expr) && e.head == :||
        tvl_or(e.args[1], e.args[2])
    elseif isa(e, Expr) && e.head == :&&
        tvl_and(e.args[1], e.args[2])
    else
        e
    end
end

"""
# Description

Recursively rewrite an arbitary expression, `e`, to use TVL.

# Arguments

1. `e::Any`: An expression.

# Return Values

1. `e′::Any`: A new expression in which three-valued logic has been substituted
    throughout the entire expression tree.

# Examples

```
julia> tvl(:(a || b))
quote
    #= REPL[1]:28 =#
    let var"##259" = a
        #= REPL[1]:29 =#
        if var"##259" === true
            #= REPL[1]:30 =#
            true
        elseif #= REPL[1]:31 =# var"##259" === false
            #= REPL[1]:32 =#
            b
        else
            #= REPL[1]:34 =#
            if b === true
                #= REPL[1]:35 =#
                true
            else
                #= REPL[1]:37 =#
                missing
            end
        end
    end
end
```
"""
function tvl(@nospecialize(e::Any))
    postwalk(rewrite_tvl, e)
end

"""
# Description

A macro that rewrites the input to use three-valued logic.

# Arguments

1. `e::Any`: An expression.

# Return Values

1. `e′::Any`: A new expression in which three-valued logic is used.

# Examples

```
julia> @tvl(false || missing)
missing

julia> @tvl(true || missing)
true
```
"""
macro tvl(e)
    esc(tvl(e))
end


text:
@testset "err_if_given_invalid_K" begin
    @test err_if_given_invalid_K(5) == 5
    @test_throws ArgumentError err_if_given_invalid_K(0)
    @test_throws ArgumentError err_if_given_invalid_K(-5)
end

@testset "replace!" begin
    y = [1, 2, 2, 2, 3]
    z1 = 1:5
    z2 = 0:3
    r1 = -2:2
    r2 = 2:5
    @test_throws DimensionMismatch _replace!(y, z1, r2)
    @test_throws DimensionMismatch _replace!(y, z2, r1)
    @test _replace!(deepcopy(y), z1, r1) == Base.replace!(deepcopy(y), (z1 .=> r1)...)
    @test _replace!(deepcopy(y), z2, r2) == Base.replace!(deepcopy(y), (z2 .=> r2)...)
end

@testset "combine_weights" begin
    rng = StableRNG(100)
    knn_weights = rand(rng, 5, 3)
    sample_weights = rand(rng, 10)
    v = [sort!(rand(rng, 1:5, 3)) for _ in 1:5]
    idxs_matrix = NeighIndsMatrix(v)
    @test combine_weights(idxs_matrix, knn_weights, nothing) === knn_weights
    @test combine_weights(idxs_matrix, nothing, nothing) == Fill(1, size(idxs_matrix))
    @test combine_weights(idxs_matrix, nothing, sample_weights) == sample_weights[idxs_matrix]
    @test combine_weights(
        idxs_matrix,
        knn_weights,
        sample_weights
    ) == (knn_weights .* sample_weights[idxs_matrix])

end

@testset "_sum" begin
    rng = StableRNG(15)
    A = rand(rng, 5, 3)
    fill_value = 2
    fill_dims = (5, 3)
    x = Fill(fill_value, fill_dims)

    @test _sum(A, dims=:)  == sum(A, dims=:)
    @test _sum(A, dims=1)  == sum(A, dims=1)
    @test _sum(A, dims=2)  == sum(A, dims=2)
    @test _sum(A, dims=3)  == sum(A, dims=3)
    @test  _sum(x, dims=:) == fill_value * (*(fill_dims...))
    @test _sum(x, dims=1) == Fill(fill_value * fill_dims[1], (1, fill_dims[2]))
    @test _sum(x, dims=2) == Fill(fill_value * fill_dims[2], (fill_dims[1], 1))
    @test _sum(x, dims=3) === x
end

@testset "scale" begin
    rng = StableRNG(25)
    A = rand(rng, 5, 3)
    B = rand(rng, Int, 5, 3)

    # Create 4 copies due to mutation
    A1 = deepcopy(A)
    A2 = deepcopy(A)
    A3 = deepcopy(A)
    A4 = deepcopy(A)

    @test scale(A1, dims=:) === (A1 ./= sum(A1, dims=:))
    @test scale(A2, dims=1) === (A2 ./= sum(A2, dims=1))
    @test scale(A3, dims=2) === (A3 ./= sum(A3, dims=2))
    @test scale(A4, dims=3) === (A4 ./= sum(A4, dims=3))

    @test scale(B, dims=:) == (B ./ sum(B, dims=:))
    @test scale(B, dims=1) == (B ./ sum(B, dims=1))
    @test scale(B, dims=2) == (B ./ sum(B, dims=2))
    @test scale(B, dims=3) == (B ./ sum(B, dims=3))
    @test !(scale(B, dims=3) === (B ./ sum(B, dims=3)))
end

@testset "matrices" begin
    # create matrices to test on.
    rng = StableRNG(17)
    n, p = (5, 3)
    vdists = [rand(rng, 3) for _ in 1:5]
    vidxs = [rand(rng, UInt, 3) for _ in 1:5]
    vr1 = rand(rng, p)
    r1matrix = RankOneMatrix(vr1, n)
    idxs_matrix = NeighIndsMatrix(vidxs)
    dists_matrix = NeighDistsMatrix(vdists)

    @test idxs_matrix isa NeighIndsMatrix
    @test dists_matrix isa NeighDistsMatrix
    @test r1matrix isa RankOneMatrix

    # check size
    @test size(idxs_matrix) == (n, p)
    @test size(dists_matrix) == (n, p)
    @test size(r1matrix) == (n, p)

    # check `getindex` and `setindex!`
    i = rand(rng, 1:n)
    j = rand(rng, 1:p)

    @test idxs_matrix[i, j] == vidxs[i][j]
    idxs_matrix[i, j] = 0
    @test idxs_matrix[i, j] == 0
    @test view(idxs_matrix, i, :) == vidxs[i]
    @test_throws BoundsError view(idxs_matrix, 6, :)
    @test_throws BoundsError idxs_matrix[i, 4]
    @test_throws BoundsError (idxs_matrix[i, 4] = 0)

    @test dists_matrix[i, j] == vdists[i][j]
    dists_matrix[i, j] = 0
    @test dists_matrix[i, j] == 0
    @test view(dists_matrix, i, :) == vdists[i]
    @test_throws BoundsError view(dists_matrix, 6, :)
    @test_throws BoundsError dists_matrix[i, 4]
    @test_throws BoundsError (dists_matrix[i, 4] = 0)

    @test r1matrix[i, j] == vr1[j]
    @test_throws BoundsError r1matrix[i, 4]
    @test_throws BoundsError r1matrix[8, j]

    # check `Base.unaliascopy`
    unaliascopy_idxs = Base.unaliascopy(idxs_matrix)
    @test unaliascopy_idxs == NeighIndsMatrix(
        [copy(@inbounds(vidxs[i])) for i in eachindex(vidxs)]
    )
    @test !(unaliascopy_idxs === idxs_matrix)

    unaliascopy_dists = Base.unaliascopy(dists_matrix)
    @test unaliascopy_dists == NeighDistsMatrix(
        [copy(@inbounds(vdists[i])) for i in eachindex(vdists)]
    )
    @test !(unaliascopy_dists === dists_matrix)

    unaliascopy_r1 = Base.unaliascopy(r1matrix)
    @test unaliascopy_r1 == RankOneMatrix(copy(vr1), n)
    @test !(unaliascopy_r1 === r1matrix)

    # check `IndexStyle`
    @test Base.IndexStyle(idxs_matrix) === Base.IndexCartesian()
    @test Base.IndexStyle(dists_matrix) === Base.IndexCartesian()
    @test Base.IndexStyle(r1matrix) === Base.IndexCartesian()
end

@testset "check_one_based_indexing" begin
    rng = StableRNG(559)
    x = rand(rng, 5)
    X = rand(rng, 5, 3)
    x1 = OffsetArray(x, 0:4)
    X1 = OffsetArray(X, 0:4, 0:2)

    @test check_onebased_indexing("input", x) === nothing
    @test check_onebased_indexing("input", X) === nothing
    @test_throws ArgumentError check_onebased_indexing("input", x1)
    @test_throws ArgumentError check_onebased_indexing("input", X1)
end


text:
# Locally Linear Embedding (LLE)
# ------------------------
# Nonlinear dimensionality reduction by locally linear embedding,
# Roweis, S. & Saul, L., Science 290:2323 (2000)

"""

    LLE{NN <: AbstractNearestNeighbors, T <: Real} <: AbstractDimensionalityReduction

The `LLE` type represents a locally linear embedding model constructed for `T` type data constructed with a help of the `NN` nearest neighbor algorithm.
"""
struct LLE{NN <: AbstractNearestNeighbors, T <: Real} <: AbstractDimensionalityReduction
    nearestneighbors::NN
    component::AbstractVector{Int}
    λ::AbstractVector{T}
    proj::Projection{T}
end

## properties
outdim(R::LLE) = size(R.proj, 1)
eigvals(R::LLE) = R.λ
neighbors(R::LLE) = R.nearestneighbors.k
vertices(R::LLE) = R.component

## show
summary(io::IO, R::LLE) = print(io, "LLE(outdim = $(outdim(R)), neighbors = $(neighbors(R)))")

## interface functions
"""
    fit(LLE, data; k=12, maxoutdim=2, nntype=BruteForce, tol=1e-5)

Fit a locally linear embedding model to `data`.

# Arguments
* `data`: a matrix of observations. Each column of `data` is an observation.

# Keyword arguments
* `k`: a number of nearest neighbors for construction of local subspace representation
* `maxoutdim`: a dimension of the reduced space.
* `nntype`: a nearest neighbor construction class (derived from `AbstractNearestNeighbors`)
* `tol`: an algorithm regularization tolerance

# Examples
```julia
M = fit(LLE, rand(3,100)) # construct LLE model
R = transform(M)          # perform dimensionality reduction
```
"""
function fit(::Type{LLE}, X::AbstractMatrix{T};
             k::Int=12, maxoutdim::Int=2, nntype=BruteForce, tol::Real=1e-5) where {T<:Real}
    # Construct NN graph
    NN = fit(nntype, X, k)
    D, E = knn(NN, X)
    _, C = largest_component(SimpleWeightedGraph(adjmat(D,E)))

    X = @view X[:, C]
    n = length(C)

    # Correct indexes of neighbors if more then one connected component
    Ec = E
    if size(E,2) != n
        R = Dict(zip(C, collect(1:n)))
        Ec = zeros(Int,k,n)
        for i in 1 : n
            Ec[:,i] = map(j->get(R,j,C[i]), E[:,C[i]])
        end
    end

    if k > maxoutdim
        @warn("k > maxoutdim: regularization will be used")
    else
        tol = 0
    end

    # Reconstruct weights and compute embedding:
    # M = (I - w)'(I - w) = I - w'I - Iw + w'w
    M = spdiagm(0 => fill(one(T), n))
    Ones = fill(one(T), k, 1)
    for i in 1 : n
        J = Ec[:,i]
        Z = view(X, :, J) .- view(X, :, i)
        G = transpose(Z)*Z
        G += I * tol # regularize
        w = vec(G \ Ones)
        w ./= sum(w)
        ww = w*transpose(w)
        for (l, j) in enumerate(J)
            M[i,j] -= w[l]
            M[j,i] -= w[l]
            for (m, jj) in enumerate(J)
                M[j,jj] = ww[l,m]
            end
        end
    end

    λ, V = decompose(M, maxoutdim)
    return LLE{nntype, T}(NN, C, λ, rmul!(transpose(V), sqrt(n)))
end

"""
    transform(R::LLE)

Transforms the data fitted to the LLE model `R` into a reduced space representation.
"""
transform(R::LLE) = R.proj


text:
using PrettyTables
function PrintEstimationResults(results, names)
    results = [names[:] results]
    clabels = ["parameter", "estimate", "st. err", "t-stat", "p-value"]
    h10 = Highlighter((results,i,j) -> (j==5) && (results[i,5] < 0.1), crayon"yellow");
    h5 = Highlighter((results,i,j) -> (j==5) && (results[i,5] < 0.05), crayon"blue");
    h1 = Highlighter((results,i,j) -> (j==5) && (results[i,5] < 0.01), crayon"green");
    pretty_table(results; header=clabels, formatters=ft_printf("%12.5f"), highlighters=(h1, h5, h10))
end

max_stars_repo_name: mcreel/EconometricsNotes

max_issues_repo_name: mcreel/EconometricsNotes

max_forks_repo_name: mcreel/EconometricsNotes


text:
module Splines2

export ns, ns_, bs, bs_, is, is_, ms, ms_

using OffsetArrays
using LinearAlgebra
using Statistics

# utility functions
zeroArray(a :: Array{T}) where T<:Real = OffsetArray(a,(size(a).*0).-1)

# type definitions
abstract type AbstractSplineBasis{T<:Real} end
# Note: we have used OffsetArray for converting from C code
mutable struct SplineBasis{T<:Real} <: AbstractSplineBasis{T}
    order::Int32 # order of the spline
    ordm1::Int32 # order -1 (3 for cubic splines
    nknots::Int32 # number of knots
    curs::Int32 # current position in knots vector
    boundary::Int32
    ncoef::Int32 # number of coefficients
    ldel::OffsetArray{T,1} # differences from knots on the left
    rdel::OffsetArray{T,1} # differences from knots on the right
    knots::OffsetArray{T,1} # knot vector
    a::OffsetArray{T,1} # scratch vector
end
mutable struct BSplineBasis{T<:Real} <: AbstractSplineBasis{T}
    spline_basis::SplineBasis{T}
    boundary_knots::Tuple{T,T}
    interior_knots::Union{Array{T,1}, Nothing}
    intercept::Bool
    df::Int32
end
mutable struct NSplineBasis{T<:Real} <: AbstractSplineBasis{T}
    b_spline_basis::BSplineBasis{T}
    qmat::Array{T,2}
    tl0::Array{T,1}
    tl1::Array{T,1}
    tr0::Array{T,1}
    tr1::Array{T,1}
end

# constructors
function SplineBasis(knots :: Array{T,1}, order :: Int32 = 4) where T<:Real
    ordm1 = order -1
    SplineBasis(order,
                ordm1,
                length(knots),
                0,
                0,
                length(knots)-order,
                zeroArray(zeros(T,ordm1)),
                zeroArray(zeros(T,ordm1)),
                zeroArray(knots),
                zeroArray(zeros(T,order)))
end
function BSplineBasis(boundary_knots :: Tuple{T,T},
                      interior_knots :: Union{Array{T,1}, Nothing} = nothing,
                      order :: Int32 = 4,
                      intercept :: Bool = false) where T<:Real
    l_interior_knots = interior_knots == nothing ? 0 : length(interior_knots)
    df = Int32(intercept) + order - 1 + l_interior_knots
    nknots = l_interior_knots + 2*order
    ncoef = nknots - order
    knots = zeros(T,nknots)
    for i=1:order
        knots[i] = boundary_knots[1]
        knots[nknots-i+1] = boundary_knots[2]
    end
    if (l_interior_knots>0)
        for i=1:l_interior_knots
            knots[i+order]=interior_knots[i]
        end
    end
    BSplineBasis(SplineBasis(knots,order), boundary_knots, interior_knots, intercept, df)
end
function NSplineBasis(boundary_knots :: Tuple{T,T},
                      interior_knots :: Union{Array{T,1}, Nothing} = nothing,
                      order :: Int32 = 4,
                      intercept :: Bool = false) :: NSplineBasis{T} where T<:Real
    bs = BSplineBasis(boundary_knots,interior_knots,order,intercept)
    co = basis(bs,[bs.boundary_knots[1],bs.boundary_knots[2]],2)
    qmat_ = LinearAlgebra.qr(transpose(co)).Q
    qmat = collect(transpose(qmat_[:, 3:size(qmat_,2)]))
    tl0 = qmat * basis(bs,boundary_knots[1])
    tl1 = qmat * basis(bs,boundary_knots[1],1)
    tr0 = qmat * basis(bs,boundary_knots[2])
    tr1 = qmat * basis(bs,boundary_knots[2],1)
    NSplineBasis(bs,
                 qmat,
                 tl0,tl1,tr0,tr1)
end

function basis(bs :: SplineBasis{T}, x :: T, ders :: Int32 = 0) where T<:Real
    function set_cursor!(x :: T)
        bs.curs = -1
        bs.boundary = T(0)
        for i = 0:(bs.nknots - 1)
	    if (bs.knots[i] >= x)
                bs.curs = i
            end
	    if (bs.knots[i] > x)
                break
            end
        end
        if (bs.curs > bs.nknots - bs.order)
	    lastLegit = bs.nknots - bs.order
	    if (x == bs.knots[lastLegit])
	        bs.boundary = 1
                bs.curs = lastLegit
            end
        end
        return bs.curs
    end
    function diff_table!(x :: T, ndiff :: Int32)
        for i = 0:(ndiff-1)
	    bs.rdel[i] = bs.knots[bs.curs + i] - x
	    bs.ldel[i] = x - bs.knots[bs.curs - (i + 1)]
        end
    end
    function slow_evaluate(x :: T)
        ti    = bs.curs
        outer = bs.ordm1
        nder = ders # copy
        if (bs.boundary>0 && nder == bs.ordm1) # value is arbitrary
	    return T(0)
        end
        while(nder>0)   # FIXME: divides by zero
            nder -= 1
            apt = 0
            lpt = ti - outer
            for inner = (outer-1):-1:0
	        bs.a[apt] = T(outer) * (bs.a[apt + 1] - bs.a[apt])/(bs.knots[lpt + outer] - bs.knots[lpt])
                apt += 1
                lpt += 1
            end
	    outer -= 1
        end
        diff_table!(x, outer)
        while(outer>0)
            outer -= 1
            apt = 0
            lpt = outer
            rpt = 0
            for inner = outer:-1:0
	        # FIXME: divides by zero
	        bs.a[apt] = (bs.a[apt + 1] * bs.ldel[lpt] + bs.a[apt] * bs.rdel[rpt])/(bs.rdel[rpt] + bs.ldel[lpt])
                lpt -= 1
                rpt += 1
                apt += 1
            end
        end
        return(bs.a[0])
    end
    # fast evaluation of basis functions
    function basis_funcs(x :: T)
        b = zeroArray(zeros(T,bs.order))
        diff_table!(x, bs.ordm1)
        b[0] = T(1)
        for j = 1:bs.ordm1
	    saved = T(0)
	    for r = 0:(j-1) # do not divide by zero
	        den = bs.rdel[r] + bs.ldel[j - 1 - r]
	        if(den != T(0))
	            term = b[r]/den
	            b[r] = saved + bs.rdel[r] * term
	            saved = bs.ldel[j - 1 - r] * term
	        else
	            if(r != T(0) || bs.rdel[r] != T(0))
	                b[r] = saved
                    end
	            saved = T(0)
	        end
	    end
	    b[j] = saved
        end
        return(b)
    end
    val = zeroArray(zeros(T,bs.ncoef))
    set_cursor!(x)
    io = bs.curs - bs.order
    if (io < 0 || io > bs.nknots)
	for j=0:(bs.order-1)
	    val[j+io] = T(0) # R_NaN
	end
    elseif (ders > 0) # slow method for derivatives
	for i = 0:(bs.order-1)
	    for j=0:(bs.order-1)
                bs.a[j] = T(0)
            end
	    bs.a[i] = T(1)
	    val[i+io] = slow_evaluate(x)
	end
    else	# fast method for value
	valtmp = basis_funcs(x)
	for i=0:(length(valtmp)-1)
	    val[i+io]=valtmp[i]
        end
    end
    return parent(val)
end

function basis(bs :: BSplineBasis{T}, x :: T, ders :: Int32 = 0) where T<:Real
    if (x < bs.boundary_knots[1] || x > bs.boundary_knots[2])
        if (x < bs.boundary_knots[1])
            k_pivot = T(0.75)*bs.boundary_knots[1]+
            T(0.25)*bs.spline_basis.knots[5-1] # 0-based
        else
            k_pivot = T(0.75)*bs.boundary_knots[2]+
            T(0.25)*bs.spline_basis.knots[length(bs.spline_basis.knots)-4-1] # 0-based
        end
        delta = x - k_pivot
        if (ders==0)
            vec = basis(bs.spline_basis,k_pivot,0) +
                basis(bs.spline_basis,k_pivot,1)*delta +
                basis(bs.spline_basis,k_pivot,2)*delta*delta/T(2.0) +
                basis(bs.spline_basis,k_pivot,3)*delta*delta*delta/T(6.0)
        elseif (ders==1)
            vec = splines.basis(bs.spline_basis,k_pivot,1) +
                splines.basis(bs.spline_basis,k_pivot,2)*delta +
                splines.basis(bs.spline_basis,k_pivot,3)*delta*delta/T(2.0)
        elseif (ders==2)
            vec = splines.basis(bs.spline_basis,k_pivot,2) +
                splines.basis(bs.spline_basis,k_pivot,3)*delta
        elseif (ders==3)
            vec = splines.basis(bs.spline_basis,k_pivot,3)
        else
            vec = k_pivot .* T(0)
        end
    else
        vec = basis(bs.spline_basis, x, ders)
    end
    if (!bs.intercept)
        vec = vec[2:length(vec)]
    end
    vec
end

function basis(ns :: NSplineBasis{T}, x::T, ders :: Int32 = 0) where T<:Real
    if (x < ns.b_spline_basis.boundary_knots[1])
        if (ders==0)
            return(ns.tl0 + (x - ns.b_spline_basis.boundary_knots[1])*ns.tl1)
        elseif (ders==1)
            return(ns.tl1)
        else
            return(ns.tl1.* T(0))
        end
    elseif (x > ns.b_spline_basis.boundary_knots[2])
        if (ders==0)
            return(ns.tr0 + (x - ns.b_spline_basis.boundary_knots[2])*ns.tr1)
        elseif (ders==1)
            return(ns.tr1)
        else
            return(ns.tr1.* T(0))
        end
    else
        return(ns.qmat*basis(ns.b_spline_basis, x, ders))
    end
end

function basis(bs :: AbstractSplineBasis{T}, x :: Array{T,1}, ders :: Int32 = 0) where T<:Real
    f(xi) = basis(bs, xi, ders)
    copy(transpose(reduce(hcat, f.(x))))
end

# utility function for processing the spline arguments
function spline_args(x :: Array{T,1};
                     boundary_knots :: Union{Tuple{T,T},Nothing} = nothing,
                     interior_knots :: Union{Array{T,1},Nothing} = nothing,
                     order :: Int32 = 4,
                     intercept :: Bool = false,
                     df :: Int32 = 3 + Int32(intercept),
                     knots :: Union{Array{T,1}, Nothing} = nothing,
                     knots_offset :: Int32 = 0) where T<:Real
    if (interior_knots != nothing && boundary_knots != nothing)
        # pass
    elseif (knots != nothing)
        boundary_knots = extrema(knots)
        interior_knots = length(knots)==2 ? nothing : knots[2:(length(knots)-1)]
    else
        if (boundary_knots == nothing)
            boundary_knots = extrema(x)
        end
        iKnots = df - order + knots_offset + 1 - Int32(intercept)
        if (iKnots>0)
            p = range(T(0), length=iKnots+2, stop=T(1))[2:(iKnots+1)]
            index = (x .>= boundary_knots[1]) .* (x .<= boundary_knots[2])
            interior_knots = Statistics.quantile(x[index], p)
        end
    end
    (boundary_knots, interior_knots)
end

"""
    bs_(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for B-splines and return a function with signature
`(x:: Array{T,1}; ders :: Int32 = 0)` for evaluation of `ders`
derivative for the splines at `x`.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int32 = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int32 = order - 1 + Int32(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines

# Examples
```jldoctest
julia> Splines2.bs_(collect(0.0:0.2:1.0), df=3)(collect(0.0:0.2:1.0))
6×3 Array{Float64,2}:
 0.0    0.0    0.0
 0.384  0.096  0.008
 0.432  0.288  0.064
 0.288  0.432  0.216
 0.096  0.384  0.512
 0.0    0.0    1.0
```
"""
function bs_(x :: Array{T,1};
            boundary_knots :: Union{Tuple{T,T},Nothing} = nothing,
            interior_knots :: Union{Array{T,1},Nothing} = nothing,
            order :: Int32 = 4,
            intercept :: Bool = false,
            df :: Int32 = order - 1 + Int32(intercept),
            knots :: Union{Array{T,1}, Nothing} = nothing,
            centre :: Union{T,Nothing} = nothing) where T<:Real
    (boundary_knots, interior_knots) =
        spline_args(x, boundary_knots=boundary_knots, interior_knots=interior_knots,
                    order=order, intercept=intercept,
                    df=df, knots=knots)
    spline = BSplineBasis(boundary_knots, interior_knots, order, intercept)
    function eval(x :: Array{T,1}; ders :: Int32 = 0)
        b = basis(spline, x, ders)
        if (centre != nothing && ders==0)
            bc = basis(spline, centre, ders)
            for i=1:size(b,1)
                b[i,:] -= bc
            end
        end
        b
    end
    eval
end

"""
    bs(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for B-splines.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int32 = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int32 = order - 1 + Int32(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines
- `ders :: Int32 = 0`: derivatives of the splines

# Examples
```jldoctest
julia> Splines2.bs(collect(0.0:0.2:1.0), df=3)
6×3 Array{Float64,2}:
 0.0    0.0    0.0
 0.384  0.096  0.008
 0.432  0.288  0.064
 0.288  0.432  0.216
 0.096  0.384  0.512
 0.0    0.0    1.0
```
"""
function bs(x :: Array{T,1}; ders :: Int32 = 0, kwargs...) where T<:Real
    bs_(x; kwargs...)(x, ders=ders)
end

"""
    ns_(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for natural B-splines and return a function with signature
`(x:: Array{T,1}; ders :: Int32 = 0)` for evaluation of `ders`
derivative for the splines at `x`.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int32 = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int32 = order - 1 + Int32(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines

# Examples
```jldoctest
julia> Splines2.ns_(collect(0.0:0.2:1.0), df=3)(collect(0.0:0.2:1.0))
6×3 Array{Float64,2}:
  0.0       0.0        0.0
 -0.100444  0.409332  -0.272888
  0.102383  0.540852  -0.359235
  0.501759  0.386722  -0.172481
  0.418872  0.327383   0.217745
 -0.142857  0.428571   0.714286
```
"""
function ns_(x :: Array{T,1};
            boundary_knots :: Union{Tuple{T,T},Nothing} = nothing,
            interior_knots :: Union{Array{T,1},Nothing} = nothing,
            order :: Int32 = 4,
            intercept :: Bool = false,
            df :: Int32 = order - 3 + Int32(intercept),
            knots :: Union{Array{T,1}, Nothing} = nothing,
            centre :: Union{T,Nothing} = nothing) where T<:Real
    (boundary_knots, interior_knots) =
        spline_args(x, boundary_knots=boundary_knots, interior_knots=interior_knots,
                    order=order, intercept=intercept, df=df, knots=knots, knots_offset=2)
    spline = NSplineBasis(boundary_knots, interior_knots, order, intercept)
    function eval(x :: Array{T,1}; ders :: Int32 = 0)
        b = basis(spline, x, ders)
        if (centre != nothing && ders==0)
            bc = basis(spline, centre, ders)
            for i=1:size(b,1)
                b[i,:] -= bc
            end
        end
        b
    end
    eval
end

"""
    ns(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for natural B-splines.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int32 = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int32 = order - 1 + Int32(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines
- `ders :: Int32 = 0`: derivatives of the splines

# Examples
```jldoctest
julia> Splines2.ns(collect(0.0:0.2:1.0), df=3)
6×3 Array{Float64,2}:
  0.0       0.0        0.0
 -0.100444  0.409332  -0.272888
  0.102383  0.540852  -0.359235
  0.501759  0.386722  -0.172481
  0.418872  0.327383   0.217745
 -0.142857  0.428571   0.714286
```
"""
function ns(x :: Array{T,1}; ders :: Int32 = 0, kwargs...) where T<:Real
    ns_(x; kwargs...)(x, ders=ders)
end

"""
    is_(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for I-splines and return a function with signature
`(x:: Array{T,1}; ders :: Int32 = 0)` for evaluation of `ders`
derivative for the splines at `x`.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int32 = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int32 = order - 1 + Int32(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)

# Examples
```jldoctest
julia> Splines2.is_(collect(0.0:0.2:1.0), df=3)(collect(0.0:0.2:1.0))
6×3 Array{Float64,2}:
 0.0     0.0     0.0
 0.1808  0.0272  0.0016
 0.5248  0.1792  0.0256
 0.8208  0.4752  0.1296
 0.9728  0.8192  0.4096
 1.0     1.0     1.0
```
"""
function is_(x :: Array{T,1};
            boundary_knots :: Union{Tuple{T,T},Nothing} = nothing,
            interior_knots :: Union{Array{T,1},Nothing} = nothing,
            order :: Int32 = 4,
            intercept :: Bool = false,
            df :: Int32 = order - 1 + Int32(intercept),
            knots :: Union{Array{T,1}, Nothing} = nothing) where T<:Real
    (boundary_knots, interior_knots) =
        spline_args(x, boundary_knots=boundary_knots, order=order+1,
                    intercept=intercept, df=df+1, knots=knots)
    spline = BSplineBasis(boundary_knots, interior_knots, order+1, false)
    knots = parent(spline.spline_basis.knots)
    findj(x) = interior_knots==nothing ? (order+1) : searchsortedlast(knots,x)
    function eval(x :: Array{T,1}; ders :: Int32 = 0)
        b = basis(spline,x,ders)
        (nrow,ncol) = size(b)
        for i=1:nrow
            js = findj(x[i])
            for j in ncol:-1:1
                if j>js
                    b[i,j]=T(0)
                elseif j==ncol
                    b[i,j]=b[i,j]
                else
                    b[i,j]=b[i,j]+b[i,j+1]
                end
            end
            for j in ncol:-1:1
                if j < js-(order+1)
                    b[i,j]=T(1)
                end
            end
        end
        intercept ? b : b[:, 2:size(b,2)]
    end
    eval
end

"""
    is(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for I-splines.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int32 = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int32 = order - 1 + Int32(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `ders :: Int32 = 0`: derivatives of the splines

# Examples
```jldoctest
julia> Splines2.is(collect(0.0:0.2:1.0), df=3)
6×3 Array{Float64,2}:
 0.0     0.0     0.0
 0.1808  0.0272  0.0016
 0.5248  0.1792  0.0256
 0.8208  0.4752  0.1296
 0.9728  0.8192  0.4096
 1.0     1.0     1.0
```
"""
function is(x :: Array{T,1}; ders :: Int32 = 0, kwargs...) where T<:Real
    is_(x; kwargs...)(x, ders=ders)
end

"""
    ms_(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for M-splines and return a function with signature
`(x:: Array{T,1}; ders :: Int32 = 0)` for evaluation of `ders`
derivative for the splines at `x`.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int32 = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int32 = order - 1 + Int32(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines

# Examples
```jldoctest
julia> Splines2.ms_(collect(0.0:0.2:1.0), df=3)(collect(0.0:0.2:1.0))
6×3 Array{Float64,2}:
 0.0    0.0    0.0
 1.536  0.384  0.032
 1.728  1.152  0.256
 1.152  1.728  0.864
 0.384  1.536  2.048
 0.0    0.0    4.0
```
"""
function ms_(x :: Array{T,1};
            boundary_knots :: Union{Tuple{T,T},Nothing} = nothing,
            interior_knots :: Union{Array{T,1},Nothing} = nothing,
            order :: Int32 = 4,
            intercept :: Bool = false,
            df :: Int32 = order - 1 + Int32(intercept),
            knots :: Union{Array{T,1}, Nothing} = nothing,
            centre :: Union{T,Nothing} = nothing) where T<:Real
    (boundary_knots, interior_knots) =
        spline_args(x, boundary_knots=boundary_knots, interior_knots=interior_knots,
                    order=order, intercept=intercept, df=df, knots=knots)
    spline = BSplineBasis(boundary_knots, interior_knots, order, intercept)
    function eval(x :: Array{T,1}; ders :: Int32 = 0)
        b = basis(spline, x, ders)
        knots = parent(spline.spline_basis.knots)
        function loop(j)
            denom = knots[j+order]-knots[j]
            denom > T(0) ? order/denom : T(0)
        end
        transCoef = map(loop,1:(length(knots)-order))
        if (!intercept)
            transCoef = transCoef[2:length(transCoef)]
        end
        if (centre != nothing && ders==0)
            bc = basis(spline, centre, ders)
            for i=1:size(b,1)
                b[i,:] -= bc
            end
        end
        for i=1:size(b,1)
            b[i,:] .*= transCoef
        end
        b
    end
    eval
end

"""
    ms(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for M-splines.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int32 = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int32 = order - 1 + Int32(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines
- `ders :: Int32 = 0`: derivatives of the splines

# Examples
```jldoctest
julia> Splines2.ms(collect(0.0:0.2:1.0), df=3)
6×3 Array{Float64,2}:
 0.0    0.0    0.0
 1.536  0.384  0.032
 1.728  1.152  0.256
 1.152  1.728  0.864
 0.384  1.536  2.048
 0.0    0.0    4.0
```
"""
function ms(x :: Array{T,1}; ders :: Int32 = 0, kwargs...) where T<:Real
    ms_(x; kwargs...)(x, ders=ders)
end

end # module


text:
module UserFileIO

export loadCSV, saveCSV!, saveJLD!, loadJLD, printLog!, collectMC
export printFixCostBounds
using TypeDef, Helper
using Printf, Dates, JLD, CSV, DataFrames, DelimitedFiles, Statistics

function loadCSV(
  mcPath,
  dim::dimType)

  fname = ["eq1d","eq2d","eq3d"]
  fullName = string(mcPath,b,"/",fname[1],".csv")
  DF = CSV.File(fullName; header=false) |> DataFrame!
  nBS = length(unique(DF[:,1]))
  eqArray = Array{Array{equilType,1},1}(undef,nBS)
  for n=1:size(DF,1)
    b = DF[n,1]
    q = DF[n,2]
    b0 = minimum(DF[:,1])
    eq = allocateEquil(dim)
    eq[q].csplus .= DF[n,2 .+ (1:dim.qtr)]
    eq[q].fail   .= DF[n,2+dim.qtr .+ (1:dim.qtr)]
    eqArray[b-b0+1] = deepcopy(eq)
  end
  return eqArray

end

function saveCSV!(
  b   ::Int64,
  eq  ::Array{equilType,1},
  dim ::dimType,
  env ::envType)

  io3d = open(string(env.pathOut,"eq3d.csv"),"a")
  for q in eachindex(eq)
    for k=1:2
      for i=1:dim.product
        @printf(io3d,"%3d, %2d, %1d, %2d",b,q,k,i)
        for t=1:dim.qtr
          @printf(io3d,", %15.2f",eq[q].profit[i,t,k])
        end
        for t=1:dim.qtr
          @printf(io3d,", %10.4f",eq[q].markup[i,t,k])
        end
        @printf(io3d,"\n")
      end
    end
  end
  close(io3d)

  io2d = open(string(env.pathOut,"eq2d.csv"),"a")
  for q in eachindex(eq)
    for i=1:dim.product
      @printf(io2d,"%3d, %2d, %2d",b,q,i)
      for t=1:dim.qtr
        @printf(io2d,", %10.0f",eq[q].demand[i,t])
      end
      for t=1:dim.qtr
        @printf(io2d,", %10.4f",eq[q].price[i,t])
      end
      for t=1:dim.qtr
        @printf(io2d,", %10.4f",eq[q].wprice[i,t])
      end
      for t=1:dim.qtr
        @printf(io2d,", %10.4f",eq[q].mcost[i,t])
      end
      @printf(io2d,"\n")
    end
  end
  close(io2d)

  io1d = open(string(env.pathOut,"eq1d.csv"),"a")
  for q in eachindex(eq)
    @printf(io1d,"%3d, %2d",b,q)
    for t=1:dim.qtr
      @printf(io1d,", %21.6f",eq[q].csplus[t])
    end
    for t=1:dim.qtr
      @printf(io1d,", %2d",Int(eq[q].fail[t]))
    end
    @printf(io1d,"\n")
  end
  close(io1d)

end

function saveJLD!(
  eqa::Array{Array{equilType,1},1},
  env::envType,
  filename::String,
  optFull::Bool=false)

  D = Dict([("eqa",eqa),("env",env)])

  if isfile(filename)==true
    rm(filename)
  end
  jldopen(filename,"w") do file
    addrequire(file,TypeDef)
    write(file,"D",D)
  end

end

function loadJLD(
  fileJLD::String,
  optFull::Bool=false)

  D = jldopen(fileJLD,"r") do file
    read(file,"D")
  end

  if optFull==true
    eqa = D["eqa"]
    dat = D["dat"]
    dim = D["dim"]
    env = D["env"]
    mat = D["mat"]
    par = D["par"]
    set = D["set"]
    sim = D["sim"]
    indx = D["indx"]
    mask = D["mask"]
    return eqa,dat,dim,env,mat,par,set,sim,indx,mask
  else
    eqa = D["eqa"]
    env = D["env"]
    return eqa,env
  end

end

function collectMC(
  listBoot::Vector{Int64},
  model   ::Int64,
  # vInt    ::Bool,
  # mc_mvno ::Int64,
  dim     ::dimType,
  env     ::envType,
  mat     ::matType,
  par0    ::parType)

  @printf("\n--------------- Collecting Bootstrap samples for Model %2d ",model)
  @printf("---------------\n")

  δ   = zeros(dim.obs)
  eqb = Array{Array{equilType,1},1}(undef,0)
  pab = Array{parType,1}(undef,0)
  imp    = readdlm(string(env.pathImp,"/paramDraws.csv"),',')

  for b in listBoot # Batch loading
    fileJLD = string(env.pathOut,"sim",env.ver,".jld")
    eqa,env = loadJLD(fileJLD)
    for i in eachindex(eqa) # each bootstrap
      fail = false
      for j in eachindex(eqa[i])
        fail = fail | any(eqa[i][j].fail)
      end

      par = sampleParam!(b+i-1,imp,dim,env,mat,par0)
      if fail==false
        eqb = [eqb; [eqa[i]]]
        pab = [pab; par]
      else
        @printf("Failed mc %3d: α = %.4f\n",b+i-1,par.θ₂[1])
      end
    end
  end
  @printf("\nSummary: %d bootstrap samples collected.\n",length(eqb))
  return eqb, pab
end

function printLog!(
  b,
  tic,
  eq     ::Array{equilType,1},
  fileLog::String)

  fail = false
  for i=1:16
    fail = (fail | any(eq[i].fail))
  end
  fid = open(fileLog,"a")
  if fail==false
    @printf(fid,"Monte Carlo %3d succeded in %.0f sec.\n",b,time()-tic)
  else
    @printf(fid,"Monte Carlo %3d failed in %.0f sec.\n",b,time()-tic)
  end
  close(fid)
end

function printFixCostBounds(
  tabfb,
  t  ::Int64,
  lab::labType)

  @printf("%10s  %10s  %10s  %10s\n","Mean","Lower","Upper","f̄̄")
  for i=1:3
    lb = mean(tabfb[i,1,t,:])
    ub = mean(tabfb[i,2,t,:])
    f̄̄  = mean(tabfb[i,3,t,:])
    ls = std(tabfb[i,1,t,:])
    us = std(tabfb[i,2,t,:])
    f̄̄s = std(tabfb[i,3,t,:])
    d1 = mean(tabfb[i,2,t,:]-tabfb[i,1,t,:])
    d2 = mean(tabfb[i,3,t,:]-tabfb[i,1,t,:])
    s1 = std(tabfb[i,2,t,:]-tabfb[i,1,t,:])
    s2 = std(tabfb[i,3,t,:]-tabfb[i,1,t,:])
    @printf("%10s  %10.0f  %10.0f  %10.0f  %10.0f  %10.0f\n",
            lab.network[i],lb,ub,f̄̄,d1,d2)
    @printf("%10s  %10.0f  %10.0f  %10.0f  %10.0f  %10.0f\n",
            "(se)",ls,us,f̄̄s,s1,s2)
  end

end

end # end of module

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
module SphereSurfaceHistogram

using Requires

include("binning.jl")
include("utils.jl")
export SSHBinner
export partition_sphere1, partition_sphere2
export push!, append!, unsafe_append!

#include("Makie.jl")
# _DIR = @__DIR__
# const PATHS = Dict(
#     :point_gen => _DIR * "/extras/point_generator.jl",
#     :Makie => _DIR * "/extras/Makie.jl"
# )

function __init__()
    @require StaticArrays="90137ffa-7385-5640-81b9-e52037218182" begin
        using .StaticArrays
        include("extras/point_generator.jl")
        export random_unit_vector
    end
    @require GeometryBasics="5c1252a2-5f33-56bf-86c9-59e7332b4326" begin
        using .GeometryBasics
        include("extras/geom.jl")
        export face_mesh, vertex_mesh
        export bin_positions
        export line_segments_minimal, line_segments
    end
    @require AbstractPlotting="537997a7-5e4e-5d89-9595-2241ea00577e" begin
        using .AbstractPlotting
        include("extras/AbstractPlotting.jl")
        export histogram, plot_debug
    end
    @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" begin
        using .Makie
        include("extras/Makie.jl")
        export histogram, histogram!
    end
end

end # module


text:
function GumbelDistribution(arg0::jdouble, arg1::jdouble)
    return GumbelDistribution((jdouble, jdouble), arg0, arg1)
end

function cumulative_probability(obj::GumbelDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::GumbelDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

function get_location(obj::GumbelDistribution)
    return jcall(obj, "getLocation", jdouble, ())
end

function get_numerical_mean(obj::GumbelDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::GumbelDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_scale(obj::GumbelDistribution)
    return jcall(obj, "getScale", jdouble, ())
end

function get_support_lower_bound(obj::GumbelDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::GumbelDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function inverse_cumulative_probability(obj::GumbelDistribution, arg0::jdouble)
    return jcall(obj, "inverseCumulativeProbability", jdouble, (jdouble,), arg0)
end

function is_support_connected(obj::GumbelDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end


text:
#levels2.jl
#Algorithm 5.8 of Random Eigenvalues by Alan Edelman

function levels2(n)
    # Plot exact semicircle formula for GUE
    xfull = (-1:0.001:1) * sqrt(2n) * 1.3

    # Form the T matrix
    T = diagm(sqrt(1:n-1),1)
    T = T + T'

    # Do the eigen-decomposition of T, T = UVU'
    V, U = eig(T)

    # Precompute U'*e_n
    # tmp_en = U' * ((0:n-1) == n-1)'
    tmp_en = U[end, :]'

    y = Array(Float64, length(xfull))
    for i = 1:length(xfull)
        x = xfull[i]
        # generate the v vector as in (2.5)
        v = U * (tmp_en ./ (V - sqrt(2) * x))
        # multiply the normalization term
        y[i] = norm((sqrt(pi)) ^ (-1/2) * exp(-x^2/2) * v/v[1])^2
    end

    # Rescale so that "semicircle" is on [-1, 1] and area is pi/2
    (xfull/sqrt(2n), y * pi / sqrt(2n))
end


text:
# Thompson sampling.
# Based on https://arxiv.org/abs/1111.1797.

# Argument is the type of conjugate distribution.
# Assumption if needed: unit variance. (Only useful for normal distribution/priors, for now.)
mutable struct ThompsonSampling{D <: Distribution} <: Policy end

ThompsonSampling() = ThompsonSampling{Normal}() # Default prior distribution: normal.

struct ThompsonSamplingDetails <: PolicyDetails
  solver_time::Float64
end

# Beta priors: assume initially α = β = 1 (i.e. the prior state is to have observed one failure and one success
# for each arm).
function _ts_sample(::ThompsonSampling{Beta}, state::State{T}) where T
  return Dict(arm => rand(Beta(state.arm_counts[arm] * state.arm_average_reward[arm] + 1,
                               state.arm_counts[arm] * (1 - state.arm_average_reward[arm]) + 1))
                     for arm in keys(state.arm_counts))
end

# Normal priors: assume that the arms have a normal distribution with unknown average and known unit variance.
# Assume that the initial average is zero. This simplifies the formulae a lot!
function _ts_sample(::ThompsonSampling{Normal}, state::State{T}) where T
  n(arm) = state.arm_counts[arm]
  tr(arm) = state.arm_reward[arm]
  return Dict(arm => rand(Normal(tr(arm) / (1 + n(arm)), 1 / (1 + n(arm)))) for arm in keys(state.arm_counts))
end

function choose_action(instance::CombinatorialInstance{T}, algo::ThompsonSampling{D}, state::State{T}; with_trace::Bool=false) where {T, D}
  weights = if any(v == 0 for v in values(state.arm_counts))
    # Initialisation step: for each random variable (i.e. reward source), try arms, preferring those that have never been tested.
    Dict(arm => (state.arm_counts[arm] == 0.0) ? 1.0 : 0.0 for arm in keys(state.arm_counts))
  else
    # Normal step: just sample the prior distribution.
    _ts_sample(algo, state)
  end

  # Avoid having all weights negative, this will yield a zero solution.
  if all(values(weights) .<= 0)
    min_val = maximum(values(weights))
    for k in keys(weights)
      weights[k] -= 2 * min_val
    end
  end

  # Actually find the solution to play.
  t0 = time_ns()
  sol = solve_linear(instance, weights)
  t1 = time_ns()

  if with_trace
    return sol, ThompsonSamplingDetails((t1 - t0) / 1_000_000_000)
  else
    return sol
  end
end


text:
@connector function Pin(;name)
    sts = @variables begin
        v(t)                    # Potential at the pin [V]
        i(t), [connect=Flow]    # Current flowing into the pin [A]
    end
    ODESystem(Equation[], t, sts, [], name=name, defaults=Dict(v=>1.0, i=>1.0))
end
@doc """
    Pin(; name)

A pin in an analog circuit.

# States:
- `v(t)`: [`V`] The voltage at this pin
- `i(t)`: [`A`] The current passing through this pin
""" Pin

"""
    OnePort(; name, v_start=0.0, i_start=0.0)

Component with two electrical pins `p` and `n` and current `i` from `p` to `n`.

# States:
- `v(t)`: [`V`] The voltage across component `p.v - n.v`
- `i(t)`: [`A`] The current passing through positive pin

# Parameters:
- `v_start`: [`V`] Initial voltage across the component
- `i_start`: [`A`] Initial current through the component

# Connectors:
- `p` Positive pin
- `n` Negative pin
"""
function OnePort(;name, v_start=0.0, i_start=0.0)
    @named p = Pin()
    @named n = Pin()
    sts = @variables begin
        v(t)=v_start
        i(t)=i_start
    end
    eqs = [
        v ~ p.v - n.v
        0 ~ p.i + n.i
        i ~ p.i
    ]
    return compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

@connector function DigitalPin(; name)
    @variables val(t) v(t) i(t)
    eqs = [
        val ~ IfElse.ifelse((0.0 <= v) & (v <= 0.8) | (2.0 <= v) & (v <= 5.0),
                                IfElse.ifelse(v > 2.0, 1, 0), X)
    ]
    ODESystem(Equation[], t, [val, v, i], [], defaults=Dict(val=>0, i=>0), name=name)
end
@doc """
    DigitalPin(; name)

A pin in a digital circuit.

# States:
- `v(t)`: [`V`] The voltage at this pin
- `i(t)`: [`A`] The current passing through this pin
- `val(t)`: The binary value of the pin at this point. A voltage from `0V` to `0.8V` is a binary value of `0`.
A voltage in the range `2.0V` to `5.0V` is `1`. Any other value is `X`.
""" DigitalPin


text:
using Distributed
addprocs(16)
@everywhere include("../src/get_unstable_sens.jl")
s = LinRange(0.,1.,25)
get_unstable_sens(s)

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
include("test_analytic_functions.jl")
include("test_constants.jl")
include("test_operators.jl")
include("test_lasers.jl")
include("test_vibrational_modes.jl")
include("test_ions.jl")
include("test_ion_configurations.jl")
include("test_traps.jl")
include("test_hamiltonians.jl")
include("test_dynamics.jl")


text:
# --------------------------------------------------------------

# Methods to automatically generate colorschemes for color selection based on
# background color and a short list of seed colors

# here are some magic constants that could be changed if you really want
const _lightness_darkbg = [80.0]
const _lightness_lightbg = [60.0]
const _lch_c_const = [60]

adjust_lch(color, l, c) = convert(RGBA{Float64}, LCHab(l, c, convert(LCHab, color).h))

function lightness_from_background(bgcolor)
    bglight = convert(LCHab, bgcolor).l
    bglight < 50.0 ? _lightness_darkbg[1] : _lightness_lightbg[1]
end

function generate_colorscheme(
    bgcolor = plot_color(:white);
    color_bases = plot_color([colorant"steelblue", colorant"orangered"]),
    lightness = lightness_from_background(bgcolor),
    chroma = _lch_c_const[1],
    n = 17,
)
    seed_colors = vcat(bgcolor, map(c -> adjust_lch(c, lightness, chroma), color_bases))
    seed_colors = convert(Vector{RGB{Float64}}, seed_colors)
    colors = distinguishable_colors(
        n,
        seed_colors,
        lchoices = Float64[lightness],
        cchoices = Float64[chroma],
        hchoices = range(0; stop = 340, length = 20),
    )[2:end]
    ColorScheme(colors)
end

# ----------------------------------------------------------------------------------

function getpctrange(n::Int)
    n > 0 || error()
    n == 1 && return zeros(1)
    zs = [0.0, 1.0]
    for i in 3:n
        sorted = sort(zs)
        diffs = diff(sorted)
        widestj = 0
        widest = 0.0
        for (j, d) in enumerate(diffs)
            if d > widest
                widest = d
                widestj = j
            end
        end
        push!(zs, sorted[widestj] + 0.5 * diffs[widestj])
    end
    zs
end

function get_zvalues(n::Int)
    offsets = getpctrange(ceil(Int, n / 4) + 1) / 4
    offsets = vcat(offsets[1], offsets[3:end])
    zvalues = Float64[]
    for offset in offsets
        append!(zvalues, offset .+ [0.0, 0.5, 0.25, 0.75])
    end
    vcat(zvalues[1], 1.0, zvalues[2:(n - 1)])
end

# ----------------------------------------------------------------------------------

function darken(c, v = 0.1)
    rgba = convert(RGBA, c)
    r = max(0, min(rgba.r - v, 1))
    g = max(0, min(rgba.g - v, 1))
    b = max(0, min(rgba.b - v, 1))
    RGBA(r, g, b, rgba.alpha)
end
lighten(c, v = 0.3) = darken(c, -v)

# borrowed from http://stackoverflow.com/a/1855903:
lightness_level(c::Colorant) = 0.299 * red(c) + 0.587 * green(c) + 0.114 * blue(c)

isdark(c::Colorant) = lightness_level(c) < 0.5
islight(c::Colorant) = !isdark(c)

Base.convert(::Type{RGB}, h::Unsigned) =
    RGB([(x & 0x0000FF) / 0xFF for x in (h >> 16, h >> 8, h)]...)

make255(x) = round(Int, 255 * x)

rgb_string(c::Colorant) =
    @sprintf("rgb(%d, %d, %d)", make255(red(c)), make255(green(c)), make255(blue(c)))

function rgba_string(c::Colorant)
    @sprintf(
        "rgba(%d, %d, %d, %1.3f)",
        make255(red(c)),
        make255(green(c)),
        make255(blue(c)),
        alpha(c)
    )
end


text:
using SimpleCG

plane_x = Plane(Vector3(0, 1, 0), 0)
plane_y = Plane(Vector3(0, 0, 1), -50)
plane_z = Plane(Vector3(1, 0, 0), -20)
sphere = Sphere(Vector3(0, 10, -10), 10)

# light = DirectionalLight(White, Vector3(-1.75, -2, -1.5))
# light = PointLight(White * 2000, Vector3(30, 40, 20))
light = SpotLight(White * 2000, Vector3(30, 40, 20), Vector3(-1, -1, -1), π/9, π/6, 0.5)

canvas = Canvas(400, 400);
camera = PerspectiveCamera(Vector3(0, 10, 10), Vector3(0, 0, -1), Vector3(0, 1, 0), π/2)
scene = UnionGeometry([plane_x, plane_y, plane_z, sphere])
SimpleCG.inner_render!(canvas, scene, camera, SimpleCG.light_render!, Lights([light]))
imshow(canvas)