text:
# unit vectors and identity matrix
const e1 = SVector(1, 0, 0)
const e2 = SVector(0, 1, 0)
const e3 = SVector(0, 0, 1)
const I3 = @SMatrix [1 0 0; 0 1 0; 0 0 1]
const Z3 = @SMatrix zeros(3,3)

"""
    tilde(x)

Construct the cross product operator matrix
"""
@inline tilde(x) = @SMatrix [0 -x[3] x[2]; x[3] 0 -x[1]; -x[2] x[1] 0]

"""
    wiener_milenkovic(c)

Construct a Wiener-Milenkovic transformation matrix, given the three Wiener-Milenkovic
parameters in `c`.

Note that the corresponding rotation matrix is the transpose of this
transformation matrix.
"""
@inline function wiener_milenkovic(c)

    c0 = 2 - c'*c/8

    C = 1/(4-c0)^2*(@SMatrix [
        c0^2 + c[1]^2 - c[2]^2 - c[3]^2  2*(c[1]*c[2] + c0*c[3])          2*(c[1]*c[3]-c0*c[2]);
        2*(c[1]*c[2] - c0*c[3])          c0^2 - c[1]^2 + c[2]^2 - c[3]^2  2*(c[2]*c[3] + c0*c[1]);
        2*(c[1]*c[3] + c0*c[2])          2*(c[2]*c[3] - c0*c[1])          c0^2 - c[1]^2 - c[2]^2 + c[3]^2
        ]
    )

    return C
end

"""
    transform_properties(K, T)

Applies the transformation `T` to the stiffness or mass matrix `K`
"""
transform_properties(X, T) = [T' Z3; Z3 T'] * X * [T Z3; Z3 T]

"""
    rotation_parameter_scaling(θ)

Extract a scaling parameter which may be multiplied by the angular parameters
to yield the Wiener-Milenkovic rotation parameters.  Use of this scaling
parameter allows deflections greater than 360 degrees.
"""
@inline function rotation_parameter_scaling(θ)

    θ = real(θ)

    θnorm = sqrt(θ'*θ)

    scaling = 1

    m = div(div(θnorm, 4) + 1, 2) # number of times to add/subtract

    if θnorm > 4
        scaling -= 8*m/θnorm
    end

    return scaling
end

@inline function rotation_parameter_scaling_θ(θ)

    return rotation_parameter_scaling_θ(rotation_parameter_scaling(θ), θ)
end

@inline function rotation_parameter_scaling_θ(scaling, θ)

    return ifelse(isone(scaling), zero(θ), (1 - scaling)*θ/(θ'*θ))
end

@inline function rotation_parameter_scaling_θ_θ(θ)

    return rotation_parameter_scaling_θ_θ(rotation_parameter_scaling(θ), θ)
end

@inline function rotation_parameter_scaling_θ_θ(scaling, θ)

    return ifelse(isone(scaling), (@SMatrix zeros(eltype(θ),3,3)), (1 - scaling)/(θ'*θ)*(I - 3*θ*θ'/(θ'*θ)))
end

"""
    get_C(θ)

Returns the transformation matrix `C` given the three angular parameters in `θ`.
"""
@inline function get_C(θ)

    scaling = rotation_parameter_scaling(θ)

    c = scaling*θ

    return wiener_milenkovic(c)
end

"""
    get_C_θ([C, ] θ)

Calculate the derivative of the Wiener-Milenkovic transformation matrix `C` with
respect to each of the rotation parameters in `θ`.
"""
@inline get_C_θ(c) = get_C_θ(get_C(c), c)

@inline function get_C_θ(C, θ)

    scaling = rotation_parameter_scaling(θ)
    scaling_θ = rotation_parameter_scaling_θ(scaling, θ)

    c = scaling*θ
    c_θ = scaling_θ*θ' + scaling*I3

    c0 = 2 - c'*c/8
    c0_c = -c/4

    tmp = 1/(4-c0)^2
    tmp_c = 2/(4-c0)^3*c0_c

    C_c1 = tmp_c[1]*C/tmp + tmp*(@SMatrix [
        2*c0*c0_c[1] + 2*c[1]    2*(c[2] + c0_c[1]*c[3]) 2*(c[3] - c0_c[1]*c[2]);
         2*(c[2] - c0_c[1]*c[3])  2*c0*c0_c[1] - 2*c[1]   2*(c0_c[1]*c[1] + c0);
         2*(c[3] + c0_c[1]*c[2]) -2*(c0_c[1]*c[1] + c0)   2*c0*c0_c[1] - 2*c[1]
        ]
    )

    C_c2 = tmp_c[2]*C/tmp + tmp*(@SMatrix [
        2*c0*c0_c[2] - 2*c[2]   2*(c[1] + c0_c[2]*c[3]) -2*(c0_c[2]*c[2] + c0);
        2*(c[1] - c0_c[2]*c[3]) 2*c0*c0_c[2] + 2*c[2]    2*(c[3] + c0_c[2]*c[1]);
         2*(c0_c[2]*c[2] + c0)   2*(c[3] - c0_c[2]*c[1])  2*c0*c0_c[2] - 2*c[2]
        ]
    )

    C_c3 = tmp_c[3]*C/tmp + tmp*(@SMatrix [
         2*c0*c0_c[3] - 2*c[3]   2*(c0_c[3]*c[3] + c0)    2*(c[1] - c0_c[3]*c[2]);
        -2*(c0_c[3]*c[3] + c0)   2*c0*c0_c[3] - 2*c[3]    2*(c[2] + c0_c[3]*c[1]);
         2*(c[1] + c0_c[3]*c[2]) 2*(c[2] - c0_c[3]*c[1])  2*c0*c0_c[3] + 2*c[3]
        ]
    )

    C_θ1 = C_c1*c_θ[1,1] + C_c2*c_θ[2,1] + C_c3*c_θ[3,1]
    C_θ2 = C_c1*c_θ[1,2] + C_c2*c_θ[2,2] + C_c3*c_θ[3,2]
    C_θ3 = C_c1*c_θ[1,3] + C_c2*c_θ[2,3] + C_c3*c_θ[3,3]

    return C_θ1, C_θ2, C_θ3
end

"""
    get_C_t([C, ] θ, θ_t)

Calculate the derivative of the Wiener-Milenkovic transformation matrix `C` with
respect to the scalar parameter `t`. `θ_t` is the derivative of the angular parameter
`θ` with respect to `t`.
"""
@inline get_C_t(θ, θ_t) = get_C_t(get_C(θ), θ, θ_t)

@inline function get_C_t(C, θ, θ_t)

    scaling = rotation_parameter_scaling(θ)
    scaling_θ = rotation_parameter_scaling_θ(scaling, θ)

    c = scaling*θ
    c_t = scaling_θ'*θ_t*θ + scaling*θ_t

    c0 = 2 - c'*c/8
    c0dot = -c'*c_t/4
    tmp = 1/(4-c0)^2
    tmpdot = 2/(4-c0)^3*c0dot

    Cdot = tmpdot*C/tmp + tmp*(@SMatrix [
        2*c0*c0dot + 2*c[1]*c_t[1] - 2*c[2]*c_t[2] - 2*c[3]*c_t[3] (
        2*(c_t[1]*c[2] + c[1]*c_t[2] + c0dot*c[3] + c0*c_t[3])) (
        2*(c_t[1]*c[3] + c[1]*c_t[3] - c0dot*c[2] - c0*c_t[2]));
        2*(c_t[1]*c[2] + c[1]*c_t[2] - c0dot*c[3] - c0*c_t[3]) (
        2*c0*c0dot - 2*c[1]*c_t[1] + 2*c[2]*c_t[2] - 2*c[3]*c_t[3]) (
        2*(c_t[2]*c[3] + c[2]*c_t[3] + c0dot*c[1] + c0*c_t[1]));
        2*(c_t[1]*c[3] + c[1]*c_t[3] + c0dot*c[2] + c0*c_t[2]) (
        2*(c_t[2]*c[3] + c[2]*c_t[3] - c0dot*c[1] - c0*c_t[1])) (
        2*c0*c0dot - 2*c[1]*c_t[1] - 2*c[2]*c_t[2] + 2*c[3]*c_t[3])]
    )

    return Cdot
end

"""
    get_C_t_θ(θ, θ_t)

Calculate the derivative of the time derivative of the Wiener-Milenkovic
transformation matrix `C` with respect to `θ`.
"""
@inline function get_C_t_θ(θ, θ_t)

    θ1, θ2, θ3 = θ
    C_t_θ1 = ForwardDiff.derivative(θ1 -> get_C_t(SVector(θ1, θ2, θ3), θ_t), θ1)
    C_t_θ2 = ForwardDiff.derivative(θ2 -> get_C_t(SVector(θ1, θ2, θ3), θ_t), θ2)
    C_t_θ3 = ForwardDiff.derivative(θ3 -> get_C_t(SVector(θ1, θ2, θ3), θ_t), θ3)

    return C_t_θ1, C_t_θ2, C_t_θ3
end

"""
    get_C_t_θdot(θ)

Calculate the derivative of the time derivative of the Wiener-Milenkovic
transformation matrix `C` with respect to each of the time derivatives of `θ`.
"""
@inline function get_C_t_θdot(θ)

    Cdot_θdot1 = ForwardDiff.derivative(cdot1 -> GXBeam.get_C_t(θ, SVector(cdot1, 0, 0)), 0)
    Cdot_θdot2 = ForwardDiff.derivative(cdot2 -> GXBeam.get_C_t(θ, SVector(0, cdot2, 0)), 0)
    Cdot_θdot3 = ForwardDiff.derivative(cdot3 -> GXBeam.get_C_t(θ, SVector(0, 0, cdot3)), 0)

    return Cdot_θdot1, Cdot_θdot2, Cdot_θdot3
end

"""
    get_Q(θ)

Calculate the matrix Q as defined in the paper "Geometrically nonlinear analysis
of composite beams using Wiener-Milenković parameters" by Qi Wang and Wenbin Yu
given the rotational parameters in `θ`.
"""
@inline function get_Q(θ)

    scaling = rotation_parameter_scaling(θ)

    c = scaling*θ

    c0 = 2 - c'*c/8

    return 1/(4-c0)^2*((4 - 1/4*c'*c)*I - 2*tilde(c) + 1/2*c*c')
end

"""
    get_Q_θ(θ)
    get_Q_θ(Q, θ)

Calculate the derivative of the matrix `Q` with respect to each of the rotation
parameters in `θ`.
"""
@inline get_Q_θ(θ) = get_Q_θ(get_Q(θ), θ)

@inline function get_Q_θ(Q, θ)

    scaling = rotation_parameter_scaling(θ)
    scaling_θ = rotation_parameter_scaling_θ(scaling, θ)

    c = scaling*θ
    c_θ = scaling_θ*θ' + scaling*I3

    c0 = 2 - c'*c/8
    c0_c = -c/4

    tmp = 1/(4-c0)^2
    tmp_c = 2/(4-c0)^3*c0_c

    Q_c1 = tmp_c[1]*Q/tmp + tmp*(-c[1]/2*I - 2*tilde(e1) + (e1*c' + c*e1')/2)

    Q_c2 = tmp_c[2]*Q/tmp + tmp*(-c[2]/2*I - 2*tilde(e2) + (e2*c' + c*e2')/2)

    Q_c3 = tmp_c[3]*Q/tmp + tmp*(-c[3]/2*I - 2*tilde(e3) + (e3*c' + c*e3')/2)

    Q_θ1 = Q_c1*c_θ[1,1] + Q_c2*c_θ[2,1] + Q_c3*c_θ[3,1]
    Q_θ2 = Q_c1*c_θ[1,2] + Q_c2*c_θ[2,2] + Q_c3*c_θ[3,2]
    Q_θ3 = Q_c1*c_θ[1,3] + Q_c2*c_θ[2,3] + Q_c3*c_θ[3,3]

    return Q_θ1, Q_θ2, Q_θ3
end

"""
    get_Qinv(θ)

Calculate the matrix inverse `Qinv` as defined in the paper "Geometrically
nonlinear analysis of composite beams using Wiener-Milenković parameters" by
Qi Wang and Wenbin Yu given the rotational parameters in `θ`.
"""
@inline function get_Qinv(θ)

    scaling = rotation_parameter_scaling(θ)

    c = scaling*θ

    Qinv = (1 - 1/16*c'*c)*I3 + 1/2*tilde(c) + 1/8*c*c'

    return Qinv
end

"""
    get_Qinv_θ(θ)

Calculate the derivative of the matrix inverse `Qinv` with respect to each of the
rotation parameters in `θ`.
"""
@inline function get_Qinv_θ(θ)

    scaling = rotation_parameter_scaling(θ)
    scaling_θ = rotation_parameter_scaling_θ(scaling, θ)

    c = scaling*θ
    c_θ = scaling_θ*θ' + scaling*I3

    Qinv_c1 = -c[1]/8*I3 + 1/2*tilde(e1) + 1/8*(e1*c' + c*e1')
    Qinv_c2 = -c[2]/8*I3 + 1/2*tilde(e2) + 1/8*(e2*c' + c*e2')
    Qinv_c3 = -c[3]/8*I3 + 1/2*tilde(e3) + 1/8*(e3*c' + c*e3')

    Qinv_θ1 = Qinv_c1*c_θ[1,1] + Qinv_c2*c_θ[2,1] + Qinv_c3*c_θ[3,1]
    Qinv_θ2 = Qinv_c1*c_θ[1,2] + Qinv_c2*c_θ[2,2] + Qinv_c3*c_θ[3,2]
    Qinv_θ3 = Qinv_c1*c_θ[1,3] + Qinv_c2*c_θ[2,3] + Qinv_c3*c_θ[3,3]

    return Qinv_θ1, Qinv_θ2, Qinv_θ3
end

"""
    mul3(A_1, A_2, A_3, b)

Return the product of a 3x3x3 tensor represented by `A_1`, `A_2`, and `A_3` with
the vector `b`.
"""
@inline mul3(A_1, A_2, A_3, b) = hcat(A_1*b, A_2*b, A_3*b)

"""
    gauss_quadrature(f, a, b)

Default gauss-quadrature function used for integrating distributed loads.
"""
@inline function gauss_quadrature(f, a, b)
    h = b - a
    c = (a + b)/2
    x = h/2*GAUSS_NODES .+ c
    return h/2*GAUSS_WEIGHTS'*f.(x)
end

# this function is not used
# @inline function linf_norm_scaling(A)
#
#     droptol!(A, eps(eltype(A)))
#
#     m, n = size(A)
#     rows = rowvals(A)
#     vals = nonzeros(A)
#
#     # initialize outputs
#     rowmax = fill(typemin(eltype(A)), n)
#     rowmin = fill(typemax(eltype(A)), n)
#     @inbounds for j = 1:n
#         @inbounds for i in nzrange(A, j)
#            row = rows[i]
#            val = abs(vals[i])
#            # perform sparse wizardry...
#            if val > rowmax[row]
#                rowmax[row] = val
#            end
#            if val < rowmin[row]
#                rowmin[row] = val
#            end
#         end
#     end
#
#     # use storage already allocated for rowmax for r
#     r = rowmax
#     @inbounds for i = 1:m
#         if i in rows
#             r[i] = nextpow(2, 1/sqrt(rowmax[i]*rowmin[i]))
#         else
#             r[i] = 1
#         end
#     end
#
#     # use storage already allocated for rowmin for s
#     s = rowmin
#     @inbounds for j = 1:n
#         colrange = nzrange(A, j)
#         if isempty(colrange)
#             s[j] = 1
#         else
#             colmax = maximum(i -> r[rows[i]]*abs(vals[i]), colrange)
#             colmin = minimum(i -> r[rows[i]]*abs(vals[i]), colrange)
#             s[j] = nextpow(2, 1/sqrt(colmax*colmin))
#         end
#     end
#
#     return r, s
# end


text:
@testset "Some performance testing with memory allocations" begin

    @testset "idx performance" begin
        function f(s, T)
            res = @benchmark idx($T, $s) samples = 2 evals = 2;
            @test res.memory < 500

            res = @benchmark exp.(sum.(idx($T, $s))) samples = 2 evals = 2;
            res2 = @benchmark randn($T, $s) samples=2 evals = 2

            # maximum 10% overhead in memory
            @test res.memory < res2.memory * 1.1
        end

        f((100, 10), Float32)
        f((256, 256), ComplexF64)
    end

    @testset "rr2 performance" begin
        function f(s, T)
            res = @benchmark rr2($T, $s) samples = 2 evals = 2;
            @test res.memory < 500

            res = @benchmark exp.((rr2($T, $s))) samples = 2 evals = 2;
            res2 = @benchmark randn($T, $s) samples=2 evals = 2

            # maximum 10% overhead in memory
            @test res.memory < res2.memory * 1.1
        end

        f((256, 256), ComplexF64)
    end

    @testset "exp_ikx performance" begin
        function f(s, T)
            res = @benchmark exp_ikx($T, $s) samples = 2 evals = 2;
            @test res.memory < 1000

            res = @benchmark exp.((exp_ikx($T, $s))) samples = 2 evals = 2;
            res2 = @benchmark randn($T, $s) samples=2 evals = 2

            # maximum 10% overhead in memory
            @test res.memory < res2.memory * 1.1
        end

        f((256, 256), ComplexF64)
        f((256, 256), ComplexF32)
    end

    @testset "cpx performance" begin
        function f(s, T)
            res = @benchmark cpx($T, $s) samples = 2 evals = 2;
            @test res.memory < 1000

            res = @benchmark exp.((cpx($T, $s))) samples = 2 evals = 2;
            res2 = @benchmark randn(Complex{$T}, $s) samples=2 evals = 2

            # maximum 10% overhead in memory
            @test res.memory < res2.memory * 1.1
        end

        f((256, 256), Float32)
    end

end

ext: jl


text:
abstract type Proposal end

struct Substitution <: Proposal
    pos::Int
    base::DNA
end

struct Insertion <: Proposal
    pos::Int  # insert after this position
    base::DNA
end

struct Deletion <: Proposal
    pos::Int
end

function update_pos(p::Substitution, pos)
    return Substitution(pos, p.base)
end

function update_pos(p::Insertion, pos)
    return Insertion(pos, p.base)
end

function update_pos(p::Deletion, pos)
    return Deletion(pos)
end

struct ScoredProposal
    proposal::Proposal
    score::Score
end

"""
Whether proposals can be applied in any order.

- ensure only one substitution or deletion at every position
- ensure only one insertion at every position

"""
function are_ambiguous(ms::Vector{Proposal})
    ins_positions = []
    other_positions = []
    for i in 1:length(ms)
        m = ms[i]
        t = typeof(m)
        if t == Insertion
            push!(ins_positions, m.pos)
        else
            push!(other_positions, m.pos)
        end
    end
    ins_good = length(Set(ins_positions)) == length(ins_positions)
    others_good = length(Set(other_positions)) == length(other_positions)
    return !ins_good || !others_good
end

function get_proposal_base(seq::DNASeq, proposal::Substitution,
                           last_del_pos::Int)
    return DNASeq([proposal.base])
end

function get_proposal_base(seq::DNASeq, proposal::Insertion,
                           last_del_pos::Int)
    if proposal.pos > 0 && last_del_pos != proposal.pos
        return DNASeq([seq[proposal.pos], proposal.base])
    end
    return DNASeq([proposal.base])
end

function get_proposal_base(seq::DNASeq, proposal::Deletion,
                           last_del_pos)
    return DNASeq([])
end

struct AmbiguousProposalsError <: Exception end

function apply_proposals(seq::DNASeq,
                         proposals::Vector{Proposal})
    if are_ambiguous(proposals)
        throw(AmbiguousProposalsError())
    end
    result = DNASeq[]
    next = 1
    # keep track of deletions, so an insertion following a deletion
    # knows not to include the deleted base
    last_del_pos = 0
    # sort by position, making sure deletions go before insertions
    proposals = sort(proposals, by = p -> (p.pos, typeof(p) == Deletion ? 0 : 1))
    for p in proposals
        push!(result, seq[next:p.pos - 1])
        push!(result, get_proposal_base(seq, p, last_del_pos))
        next = p.pos + 1
        if typeof(p) == Deletion
            last_del_pos = p.pos
        end
    end
    push!(result, seq[next:end])
    return DNASeq(result...)
end

function choose_candidates(candidates::Vector{ScoredProposal}, min_dist::Int)
    final_cands = ScoredProposal[]
    posns = Set{Int}()
    for c in sort(candidates, by=(c) -> c.score, rev=true)
        if any(Bool[(abs(c.proposal.pos - p) < min_dist) for p in posns])
            continue
        end
        push!(posns, c.proposal.pos)
        push!(final_cands, c)
    end
    return final_cands
end


text:
using Triangulate
using WriteVTK
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCell
include("../../../test/useful_routines.jl")

function update_product_stress!(
    qpstress,
    basis,
    stiffness,
    transfstress,
    theta0,
    celldisp,
    points,
    jac,
    vectosymmconverter,
)
    dim = CutCell.dimension(basis)
    lambda, mu = CutCell.lame_coefficients(stiffness, +1)
    nump = size(points)[2]

    for qpidx = 1:nump
        p = points[:, qpidx]
        grad = CutCell.transform_gradient(gradient(basis, p), jac)
        NK = sum([CutCell.make_row_matrix(vectosymmconverter[k], grad[:, k]) for k = 1:dim])
        symmdispgrad = NK * celldisp

        inplanestress = (stiffness[+1] * symmdispgrad) - transfstress
        s33 = lambda * (symmdispgrad[1] + symmdispgrad[2]) - (lambda + 2mu / 3) * theta0

        numericalstress = vcat(inplanestress, s33)

        append!(qpstress, numericalstress)
    end
end

function update_parent_stress!(
    qpstress,
    basis,
    stiffness,
    celldisp,
    points,
    jac,
    vectosymmconverter,
)

    dim = CutCell.dimension(basis)
    lambda, mu = CutCell.lame_coefficients(stiffness, -1)
    nump = size(points)[2]
    for qpidx = 1:nump
        p = points[:, qpidx]
        grad = CutCell.transform_gradient(gradient(basis, p), jac)
        NK = sum([CutCell.make_row_matrix(vectosymmconverter[k], grad[:, k]) for k = 1:dim])
        symmdispgrad = NK * celldisp

        inplanestress = stiffness[-1] * symmdispgrad
        s33 = lambda * (symmdispgrad[1] + symmdispgrad[2])

        numericalstress = vcat(inplanestress, s33)

        append!(qpstress, numericalstress)
    end
end

function compute_stress_at_quadrature_points(
    nodaldisplacement,
    basis,
    stiffness,
    transfstress,
    theta0,
    cellquads,
    cutmesh,
)

    dim = CutCell.dimension(basis)
    ncells = CutCell.number_of_cells(cutmesh)
    qpstress = zeros(0)
    jac = CutCell.jacobian(cutmesh)
    vectosymmconverter = CutCell.vector_to_symmetric_matrix_converter()

    for cellid = 1:ncells
        s = CutCell.cell_sign(cutmesh, cellid)
        @assert s == -1 || s == 0 || s == +1
        if s == +1 || s == 0
            nodeids = CutCell.nodal_connectivity(cutmesh, +1, cellid)
            celldofs = CutCell.element_dofs(nodeids, dim)
            celldisp = nodaldisplacement[celldofs]
            points = cellquads[+1, cellid].points
            update_product_stress!(
                qpstress,
                basis,
                stiffness,
                transfstress,
                theta0,
                celldisp,
                points,
                jac,
                vectosymmconverter,
            )
        end
        if s == -1 || s == 0
            nodeids = CutCell.nodal_connectivity(cutmesh, -1, cellid)
            celldofs = CutCell.element_dofs(nodeids, dim)
            celldisp = nodaldisplacement[celldofs]
            points = cellquads[-1, cellid].points
            update_parent_stress!(
                qpstress,
                basis,
                stiffness,
                celldisp,
                points,
                jac,
                vectosymmconverter,
            )
        end
    end
    return reshape(qpstress, 4, :)
end

function compute_quadrature_points(cellquads, cutmesh)
    ncells = CutCell.number_of_cells(cutmesh)
    coords = zeros(2, 0)
    for cellid = 1:ncells
        s = CutCell.cell_sign(cutmesh, cellid)
        cellmap = CutCell.cell_map(cutmesh, cellid)
        @assert s == -1 || s == 0 || s == +1
        if s == +1 || s == 0
            coords = hcat(coords, cellmap(cellquads[+1, cellid].points))
        end
        if s == -1 || s == 0
            coords = hcat(coords, cellmap(cellquads[-1, cellid].points))
        end
    end
    return coords
end

function onboundary(x, L, W)
    return x[2] ≈ 0.0 || x[1] ≈ L || x[2] ≈ W || x[1] ≈ 0.0
end

function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
    a = zeros(3, 3)
    a[1, 1] = inradius
    a[1, 2] = -inradius
    a[1, 3] = -1.0 / inradius
    a[2, 1] = 2 * (lc + mc)
    a[2, 2] = -2 * (ls + ms)
    a[2, 3] = 2ms / inradius^2
    a[3, 2] = 2(ls + ms)
    a[3, 3] = -2ms / outradius^2
    return a
end

function analytical_coefficient_rhs(ls, ms, theta0)
    r = zeros(3)
    Ks = bulk_modulus(ls, ms)
    r[2] = -Ks * theta0
    r[3] = Ks * theta0
    return r
end

struct AnalyticalSolution
    inradius::Any
    outradius::Any
    center::Any
    A1c::Any
    A1s::Any
    A2s::Any
    ls::Any
    ms::Any
    lc::Any
    mc::Any
    theta0::Any
    function AnalyticalSolution(inradius, outradius, center, ls, ms, lc, mc, theta0)
        a = analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
        r = analytical_coefficient_rhs(ls, ms, theta0)
        coeffs = a \ r
        new(
            inradius,
            outradius,
            center,
            coeffs[1],
            coeffs[2],
            coeffs[3],
            ls,
            ms,
            lc,
            mc,
            theta0,
        )
    end
end

function radial_displacement(A::AnalyticalSolution, r)
    if r <= A.inradius
        return A.A1c * r
    else
        return A.A1s * r + A.A2s / r
    end
end

function (A::AnalyticalSolution)(x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    ur = radial_displacement(A, r)
    if ur ≈ 0.0
        [0.0, 0.0]
    else
        costheta = (x[1] - A.center[1]) / r
        sintheta = (x[2] - A.center[2]) / r
        u1 = ur * costheta
        u2 = ur * sintheta
        return [u1, u2]
    end
end

function solve_cylindrical_bc_elasticity(
    width,
    center,
    inradius,
    outradius,
    stiffness,
    transfstress,
    theta0,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor,
)

    lambda1, mu1 = CutCell.lame_coefficients(stiffness, +1)
    lambda2, mu2 = CutCell.lame_coefficients(stiffness, -1)

    analyticalsolution =
        AnalyticalSolution(inradius, outradius, center, lambda1, mu1, lambda2, mu2, theta0)

    dx = width / nelmts
    meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
    penalty = penaltyfactor / dx * meanmoduli

    basis = TensorProductBasis(2, polyorder)
    mesh = CutCell.Mesh([0.0, 0.0], [width, width], [nelmts, nelmts], basis)
    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCell.levelset_coefficients(
        x -> -circle_distance_function(x, center, inradius),
        mesh,
    )

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    facequads = CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition =
        CutCell.incoherent_interface_condition(basis, interfacequads, stiffness, cutmesh, penalty)

    displacementbc = CutCell.DisplacementCondition(
        analyticalsolution,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> onboundary(x, width, width),
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, cutmesh)
    CutCell.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )
    CutCell.assemble_incoherent_interface_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_bc!(sysmatrix, sysrhs, displacementbc, cutmesh)
    CutCell.assemble_penalty_displacement_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onboundary(x, width, width),
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    sol = matrix \ rhs

    return sol, basis, cellquads, cutmesh
end

function solve_zero_traction_bc_elasticity(
    width,
    center,
    inradius,
    outradius,
    stiffness,
    transfstress,
    theta0,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor,
)

    lambda1, mu1 = CutCell.lame_coefficients(stiffness, +1)
    lambda2, mu2 = CutCell.lame_coefficients(stiffness, -1)

    dx = width / nelmts
    meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
    penalty = penaltyfactor / dx * meanmoduli

    basis = TensorProductBasis(2, polyorder)
    mesh = CutCell.Mesh([0.0, 0.0], [width, width], [nelmts, nelmts], basis)
    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCell.levelset_coefficients(
        x -> -circle_distance_function(x, center, inradius),
        mesh,
    )

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    facequads = CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition =
        CutCell.incoherent_interface_condition(basis, interfacequads, stiffness, cutmesh, penalty)

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, cutmesh)
    CutCell.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )
    CutCell.assemble_incoherent_interface_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        cutmesh,
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    topleftnodeid = CutCell.nodes_per_mesh_side(mesh)[2]
    CutCell.apply_dirichlet_bc!(matrix, rhs, [1, topleftnodeid], 1, 0.0, 2)
    CutCell.apply_dirichlet_bc!(matrix, rhs, [1], 2, 0.0, 2)

    sol = matrix \ rhs

    return sol, basis, cellquads, cutmesh
end

K1, K2 = 247.0, 192.0
mu1, mu2 = 126.0, 87.0
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)

theta0 = -0.067
stiffness = CutCell.HookeStiffness(lambda1, mu1, lambda2, mu2)

width = 1.0
penaltyfactor = 1e2

nelmts = 11
polyorder = 2
numqp = required_quadrature_order(polyorder) + 4
center = [width / 2, width / 2]
inradius = width / 4
outradius = width

transfstress = CutCell.plane_strain_transformation_stress(lambda1, mu1, theta0)

# sol, basis, cellquads, cutmesh = solve_cylindrical_bc_elasticity(
#     width,
#     center,
#     inradius,
#     outradius,
#     stiffness,
#     transfstress,
#     theta0,
#     nelmts,
#     polyorder,
#     numqp,
#     penaltyfactor,
# )

sol, basis, cellquads, cutmesh = solve_zero_traction_bc_elasticity(
    width,
    center,
    inradius,
    outradius,
    stiffness,
    transfstress,
    theta0,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor,
)

qpstress = compute_stress_at_quadrature_points(
    sol,
    basis,
    stiffness,
    transfstress,
    theta0,
    cellquads,
    cutmesh,
)

numqp = size(qpstress)[2]
pressure = -[(qpstress[1, i] + qpstress[2, i] + qpstress[4, i])/3 for i = 1:numqp]

quadcoords = compute_quadrature_points(cellquads, cutmesh)

triin = Triangulate.TriangulateIO()
triin.pointlist = quadcoords
(triout, vorout) = triangulate("", triin)
connectivity = triout.trianglelist
cells = [
    MeshCell(VTKCellTypes.VTK_TRIANGLE, connectivity[:, i]) for i = 1:size(connectivity)[2]
]
vtkfile = vtk_grid(
    "examples/transformation-strain/cylindrical-particle/zero-traction-pressure",
    quadcoords[1, :],
    quadcoords[2, :],
    cells,
)
vtkfile["pressure"] = pressure
vtkfile["s11"] = qpstress[1,:]
vtkfile["s22"] = qpstress[2,:]
vtkfile["s12"] = qpstress[3,:]
vtkfile["s33"] = qpstress[4,:]
outfiles = vtk_save(vtkfile)


text:
# Cross-version compatibility
if VERSION < v"0.5.0-dev" # Julia v0.4
    readstring(cmd) = readall(cmd)
end

"""
This function handles automatic conversion of OS X/Linux directories
by converting / into \ on Windows platfors.
"""
@inline function parsedir(dirname)
    @static if is_windows()
        return replace(dirname, "/", "\\")
    else
        return dirname
    end
end

"""
This function selects a single implementation from `fnames` and associates it
with `fname`. It benchmarks every function in `fnames` against the single
type `Vector{T}` of length `N` and writes the resulting association into the
IOstream `file`.

    fname: name to be associated with Vectorize.(fname)
    fnames: a list of strings "Vectorize.VML.sin" to be benchmarked
    T: the type of the Vector{T} to pass to the function
    file: the IOStream to write the final string onto
    N: the length of the vector to benchmark over.
"""
function benchmarkSingleArgFunction(fname, fnames,
                                    T::DataType, file::IOStream, N::Integer)

    # dictionary to store elpased values
    val = Dict()
    println("TESTING: $(fname)(Vector{$T})")

    # Iterate over all functions
    for fstr in fnames
        f = eval(parse(fstr)) # convert string to function
        X = convert(Vector{T}, randn(N)) # function arguments
        f(X) # force compilation
        time = 0
        for i in 1:10 # benchmark ten times
            time += @elapsed f(X)
        end
        val[fstr] = time/10.0
    end

    # default to first function
    t = val[fnames[1]]
    fbest = fnames[1]

    # find fastest function
    for f in fnames
        if val[f] < t
            t = val[f]
            fbest = f
        end
    end

    # write chosen file into IOStream
    write(file, "\n$(fname)(X::Vector{$T}) = $(fbest)(X)\n")
    # docs = @doc fbest
    # write(file, "\n\"\"\"\n", docs, "\n\"\"\"\n")
    # write(file, "$(fname)(X::Vector{$T}) = $(fbest)(X)\n")
    println("BENCHMARK: $(fname)(Vector{$T}) mapped to $(fbest)()\n")
end

"""
This function selects a single implementation from `fnames` and associates it
with `fname`. It benchmarks every function in `fnames` against two `Vector{T}`'s
of length `N` and writes the resulting association into the
IOstream `file`.

    fname: name to be associated with Vectorize.(fname)
    fnames: a list of strings "Vectorize.VML.sin" to be benchmarked
    T: tuple of types of Vector{T} to pass to the function
    file: the IOStream to write the final string onto
    N: the length of the vector to benchmark over.
"""
function benchmarkTwoArgFunction(fname, fnames,
                                 T::Tuple{DataType, DataType}, file::IOStream, N::Integer)

    # dictionary of elapsed times for each implementation
    val = Dict()
    println("TESTING: $(fname)(Vector{$(T[1])}, Vector{$(T[2])})")

    # Iterate over all functions
    for fstr in fnames
        f = eval(parse(fstr)) # convert string to function
        X = convert(Vector{T[1]}, randn(N)) # function arguments
        Y = convert(Vector{T[2]}, randn(N)) # function arguments
        f(X, Y) # force compilation
        time = 0
        for i in 1:10 # benchmark ten times
            time += @elapsed f(X, Y)
        end
        val[fstr] = time/10.0
    end

    # default to first function
    t = val[fnames[1]]
    fbest = fnames[1]

    # find fastest function
    for f in fnames
        if val[f] < t
            t = val[f]
            fbest = f
        end
    end

    # write result into IOStream
    write(file, "\n$(fname)(X::Vector{$(T[1])}, Y::Vector{$(T[2])}) = $(fbest)(X, Y)\n")
    println("BENCHMARK: $(fname)(Vector{$(T[1])}, Vector{$(T[2])}) mapped to $(fbest)()\n")
end

"""
This function selects a single implementation from `fnames` and associates it
with `fname`. It benchmarks every function in `fnames` against three `Vector{T}`'s
of length `N` and writes the resulting association into the
IOstream `file`.

    fname: name to be associated with Vectorize.(fname)
    fnames: a list of strings "Vectorize.VML.sin" to be benchmarked
    T: tuple of types of Vector{T} to pass to the function
    file: the IOStream to write the final string onto
    N: the length of the vector to benchmark over.
"""
function benchmarkThreeArgFunction(fname, fnames,
                                   T::Tuple{DataType, DataType, DataType}, file::IOStream, N::Integer)

    # dictionary to store elapsed values
    val = Dict()
    println("TESTING: $(fname)(Vector{$(T[1])}, Vector{$(T[2])}, Vector{$(T[3])})")

    # Iterate over all functions
    for fstr in fnames
        f = eval(parse(fstr)) # convert string to function
        X = convert(Vector{T[1]}, randn(N)) # function arguments
        Y = convert(Vector{T[2]}, randn(N)) # function arguments
        Z = convert(Vector{T[3]}, randn(N)) # function arguments
        f(X, Y, Z) # force compilation
        time = 0
        for i in 1:10 # benchmark 10 times
            time += @elapsed f(X, Y, Z)
        end
        val[fstr] = time/10.0
    end

    # default to first function
    t = val[fnames[1]]
    fbest = fnames[1]

    # find fastest function
    for f in fnames
        if val[f] < t
            t = val[f]
            fbest = f
        end
    end

        # write result into IOStream
    write(file, "\n$(fname)(X::Vector{$(T[1])}, Y::Vector{$(T[2])}, Z::Vector{$(T[3])}) = $(fbest)(X, Y, Z)\n")
    println("BENCHMARK: $(fname)(Vector{$(T[1])}, Vector{$(T[2])}, Vector{$(T[3])}) mapped to $(fbest)()\n")
end

"""
`trycmd(cmd::Cmd, msg::String="", err::String="")::String`
This function attemps to run the shell command specified by `cmd` using `run`.
If `cmd` returns successfully, it will print `msg` using info().
If `run(cmd)` throws an exception (the shell command returns non-zero),
then the function throws an error and prints the `err` argument.

Note: This function does not return the value of the command to the caller. However,
prints to STD_OUT by the `run` command will print automatically from the shell process.
"""
function trycmd(cmd::Cmd; msg::String="", err::String="")
    try
        run(cmd)
        if msg == "" return else info(msg) end
    catch
        if err != ""
            error(err)
        end
    end
end

"""
`trycmd_read(cmd::Cmd, msg::String="", err::String="")::String`
This function attemps to run the shell command specified by `cmd` using `readall`.
If `cmd` returns successfully, it will return its result to the caller and print `msg`.
If `readall(cmd)`/`readstring(cmd)` throws an exception (the shell command returns non-zero),
then the function throws an error and prints the `err` argument.
"""
function trycmd_read(cmd::Cmd; msg::String="", err::String="")
    try
        result = readstring(cmd) # Julia v0.5 and above
        if msg == "" return result else info(msg) end
        return result
    catch
        if err != ""
            error(err)
        end
    end
end

"""
`prompt_yn(prompt::String)::Bool`

Prompts the user for a y or n response (y/n). Continually prompts
until user enters a valid response. Prompt string should be provided without
colon and whitespace.
"""
function prompt_yn(prompt::String)
    print(prompt, " [y/n]:  ")
    input = chomp(readline())
    while input != "y" && input != "n"
        print("You did not enter y or n. Please enter a valid response: [y/n]  ")
        input = chomp(readline())
    end

    return (input == "y") ? true : false
end


text:
#!/usr/bin/env julia

type MdlParams
    L::Float64
    Ts::Float64
    MdlParams(L=1.0,Ts=0.1) = new(L,Ts)
end

type MpcParams          # parameters for MPC solver
    N::Int64
    Q::Array{Float64,2}
    R::Array{Float64,2}
    MpcParams(N=0,Q=Float64[],R=Float64[]) = new(N,Q,R)
end

type MpcSol             # MPC solution output
    acc::Float64
    df::Float64
    uOL::Array{Float64}
    zOL::Array{Float64}
    J::Array{Float64}
    solverStatus::Symbol
    MpcSol(acc=0.0,df=0.0,uOL=Float64[],zOL=Float64[],J=Float64[],solverStatus=Symbol()) = new(acc,df,uOL,zOL,J,solverStatus)
end

type MpcModel
    z0::Array{JuMP.NonlinearParameter,1}
    zRef::Array{JuMP.NonlinearParameter,2}
    uRef::Array{JuMP.NonlinearParameter,2}
    zOL::Array{JuMP.Variable,2}
    uOL::Array{JuMP.Variable,2}
    costZ::JuMP.NonlinearExpression
    costU::JuMP.NonlinearExpression
    mdl::JuMP.Model

    function MpcModel(mpcParams::MpcParams, mdlParams::MdlParams)
        println("Building mpc controller .....")
        m = new()

        # get model parameters
        Ts          = mdlParams.Ts
        L           = mdlParams.L

        # get mpc parameters
        N           = mpcParams.N
        Q           = mpcParams.Q
        R           = mpcParams.R

        # Create Model
        mdl = Model(solver = IpoptSolver(print_level=0,max_cpu_time=5.10))

        # Create variables (these are going to be optimized)
        @variable( mdl, zOL[1:4, 1:(N+1)],  start = 0)       # z = x, y, psi, v
        @variable( mdl, uOL[1:2, 1:N],      start = 0)           # u = df, acc

        # Set bounds on inputs
        u_lb = [-0.6  -0.25]' * ones(1,mpcParams.N)       # lower bounds [steering, acceleration]
        u_ub = [0.6    0.25]' * ones(1,mpcParams.N)       # upper bounds

        for i=1:N
            for j=1:2
                setlowerbound(uOL[j,i], u_lb[j,i])
                setupperbound(uOL[j,i], u_ub[j,i])
            end
        end

        @NLparameter(mdl, z0[i=1:4] == 0)
        @NLparameter(mdl, zRef[j=1:4,i=1:N+1] == 0)
        @NLparameter(mdl, uRef[j=1:2,i=1:N] == 0)

        # System dynamics
        @NLconstraint(mdl, [i=1:4], zOL[i,1] == z0[i])         # initial condition

        for i=1:N
            @NLconstraint(mdl, zOL[1,i+1] == zOL[1,i] + Ts*(zOL[4,i]*cos( zOL[3,i] ))  )                # x
            @NLconstraint(mdl, zOL[2,i+1] == zOL[2,i] + Ts*(zOL[4,i]*sin( zOL[3,i] ))  )                # y
            @NLconstraint(mdl, zOL[3,i+1] == zOL[3,i] + Ts*(zOL[4,i]*tan( uOL[1,i] )/L)  )              # psi
            @NLconstraint(mdl, zOL[4,i+1] == zOL[4,i] + Ts*(uOL[2,i]))                                  # v
        end

        # Cost definitions
        # Input cost
        # ---------------------------------
        @NLexpression(mdl, costU, 0.5*sum(R[j,j] * sum((uOL[j,i] - uRef[j,i])^2 for i=1:N) for j = 1:2))

        # State cost
        # ---------------------------------
        @NLexpression(mdl, costZ, 0.5*sum(Q[j,j] * sum((zOL[j,i] - zRef[j,i])^2 for i=2:N+1) for j = 1:4))

        # Objective function
        @NLobjective(mdl, Min, costZ + costU)
        #@NLobjective(mdl, Min, costZ)

        # create first artificial solution (for warm start)
        for i=1:N+1
            setvalue(zOL[:,i],[0.0, 0.0, 0.0, 0.0])
        end
        for i=1:N
            setvalue(uOL[:,i],[0.0, 0.0])
        end

        # First solve
        println("Attemping first solve  .....")
        sol_stat = solve(mdl)
        println("Finished solve 1: $sol_stat")

        m.mdl   = mdl
        m.z0    = z0
        m.zOL   = zOL
        m.uOL   = uOL
        m.costZ = costZ
        m.costU = costU
        m.zRef  = zRef
        m.uRef  = uRef
        return m
    end
end

function SolveMpcProblem(mdl::MpcModel,mpcSol::MpcSol,zCurr::Array{Float64}, zRef::Array{Float64,2}, uRef::Array{Float64,2})
    # update current initial condition
    setvalue(mdl.z0,   zCurr)
    setvalue(mdl.zRef, zRef)
    setvalue(mdl.uRef, uRef)

    # solve mpc
    sol_status  = solve(mdl.mdl)

    # get solution
    sol_u       = getvalue(mdl.uOL)
    sol_z       = getvalue(mdl.zOL)

    mpcSol.df   = sol_u[1,1]
    mpcSol.acc  = sol_u[2,1]
    mpcSol.uOL  = sol_u
    mpcSol.zOL  = sol_z
    mpcSol.solverStatus = sol_status
end

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc


text:
begin
    using JLD2, FileIO
    using LinearAlgebra
    using Interpolations
    using Roots
    using Plots
    using Optim
    pyplot()

    CPROB=[0.0262397497796233 0.152923483594817 0.361483063911416 0.328567584707170 0.114741751017859 0.0152660804766739 0.000778286512441828; 0.0160443669891157 0.114741751017859 0.328567584707170 0.361483063911415 0.152923483594817 0.0247005562212874 0.00153919355833587; 0.00945177150556068 0.0828345049143809 0.287445156270025 0.382789297365280 0.196113758787202 0.0384369616149092 0.00292854954264210; 0.00536221880153007 0.0575309935179400 0.242023809517260 0.390165956326541 0.242023809517260 0.0575309935179400 0.00536221880153009; 0.00292854954264206 0.0384369616149092 0.196113758787202 0.382789297365280 0.287445156270025 0.0828345049143808 0.00945177150556065; 0.00153919355833585 0.0247005562212875 0.152923483594817 0.361483063911416 0.328567584707170 0.114741751017859 0.0160443669891157; 0.000778286512441806 0.0152660804766738 0.114741751017859 0.328567584707170 0.361483063911416 0.152923483594817 0.0262397497796233];

    nos=size(CPROB)[1];
    [sum(CPROB[i,:]) for i=1:nos];

    mM=CPROB'-I;
    mM[end,:]=[1 for i=1:nos];
    rhs=[floor(i/nos) for i=1:nos];
    BSP=mM\rhs; #the list of steady-state probabilities

    hours = exp.([-1.20000000000000, -0.800000000000000, -0.400000000000000, 0, 0.400000000000000, 0.800000000000000, 1.20000000000000]);
    wage=0.2;
    inc=wage*hours;
    ι=inc'*BSP; #payoff from the bond (= the aggregate output in the economy)

    #ρ=0.05
    ρ=-log(0.96)
    R=3.0;
end

function U(ξ::Float64,R::Float64)
    if ξ <= 0.
        -9e10
        else ((1/ξ^(R-1.0)-1)/(1.0 - R))
    end
end

function II(aa::Int64,k::Int64,a::Int64,tbl2::Array{Int64,2})
    if aa==tbl2[k,a]
        return 1
    else
        return 0
    end
end

function LS_method(intr::Float64,accu_pow1::Int64,accu_pow2::Int64,nos::Int64,hours::Array{Float64,1},wage::Float64,β::Float64,grda::Array{Float64,1},Lg::Int64,CPROB::Array{Float64,2},R::Float64)
    local TBL,TBL1,TBL2,tbl,tbl1,tbl2,Λ0,Λ1
    currentU=[U((1+intr)*grda[i]+wage*hours[k]-grda[j],R) for k=1:nos, i=1:Lg, j=1:Lg];
    nextV=[U((1+intr)*grda[j]+wage*hours[k],R) for k=1:nos, j=1:Lg];

    TBL=[findmax(currentU[k,i,:]+β*(CPROB[k,:]'*nextV)') for k=1:nos, i=1:Lg];
    TBL1=first.(TBL);TBL2=last.(TBL);

    tbl=[findmax(currentU[k,i,:]+β*(CPROB[k,:]'*TBL1)') for k=1:nos, i=1:Lg];
    tbl1=first.(tbl);tbl2=last.(tbl);

    while maximum(abs.(TBL2-tbl2))>0||maximum(abs.(TBL1-tbl1))>1/(10^accu_pow1)
        TBL1=tbl1;
        TBL2=tbl2;
        tbl=[findmax(currentU[k,i,:]+β*(CPROB[k,:]'*TBL1)') for k=1:nos, i=1:Lg];
        tbl1=first.(tbl);tbl2=last.(tbl);
    end

    Λ0=[Float64(1/(nos*Lg)) for k=1:nos, i=1:Lg];

    Λ1=[sum([Λ0[k,i]*CPROB[k,kk]*II(ii,k,i,tbl2) for k=1:nos, i=1:Lg]) for kk=1:nos, ii=1:Lg];

    Λ0=Λ1;
    Λ1=[sum([Λ0[k,i]*CPROB[k,kk]*II(ii,k,i,tbl2) for k=1:nos, i=1:Lg]) for kk=1:nos, ii=1:Lg];
    while maximum(abs.(Λ1-Λ0))>1/10^accu_pow2
        Λ0=Λ1;
        Λ1=[sum([Λ0[k,i]*CPROB[k,kk]*II(ii,k,i,tbl2) for k=1:nos, i=1:Lg]) for kk=1:nos, ii=1:Lg];
    end
    return Λ1,tbl2
end


text:
# This file contains the functions to execute the normal episode:

#Functions focusing on the episode execution: Value and Actions diverse
#-------------------------------------------------------------------------------
@everywhere function run_episode!(params::Params,θa::Any, θv::Any, traj::Trajectory, rms::RunMeanStd)
    dt          = params.dt
    t_horizon   = params.t_horizon
    s_init      = params.s_init + 0.1*sum(params.σ)*randn(length(params.s_init))
    l_traj      = params.l_traj
    σ           = params.σ
    γ           = params.γ_disc
    u_range     = params.u_range
    l_out       = length(u_range)
    s_init[1:2] = 0.2 .* randn(2)
    s_init[3]   = 1.2 .* randn()

    #Initialize :
    s_t1 = copy(s_init)
    s_t2 = similar(s_t1)
    totalr = 0
    rand_vec = zeros(l_out)

    for i_t = 1:l_traj

        #select for unnormalized output
        s_norm = s_t1

        #Only for Segway
        s_mod = s_norm
        #=---------------------------------
        s_mod        = zeros(8)
        s_mod[1:2]   = s_norm[1:2]
        s_mod[3]     = sin(s_norm[3])
        s_mod[4]     = cos(s_norm[3])
        s_mod[5:end] = s_norm[4:end]
        #-------------------------------=#

        #Evaluate the Network Output
        out_t = mlp_nlayer_policy_norm(θa,s_mod,params)
        val_t = mlp_nlayer_value_norm(θv,s_mod,params)[1]

        #Sample:
        rand_vec = randn(l_out)

        #Get Action with stochastic policy:
        loc_a_t = sample_action(out_t, params.σ, rand_vec)
        p_a_t   = get_π_prob(out_t,σ,loc_a_t)

        #Interact with the Environment:
        s_t2, loc_r_t, terminal = environment_call(s_t1,loc_a_t,params)

        #Save data in trajectory:
        traj.s_t[i_t,:]    .= s_norm
        traj.s_t2[i_t,:]   .= s_t2
        traj.a_t[i_t,:]    .= loc_a_t
        traj.p_a_t[i_t,:]   = p_a_t
        traj.r_t[i_t]       = loc_r_t*0.01
        traj.val_t[i_t]     = val_t
        traj.out_t[i_t,:]   = out_t
        traj.len            = i_t
        totalr             += loc_r_t

        #Check if the terminal condition is fullfilled and break
        if terminal == 1
            traj.s_t[i_t,:] .= s_t1
            totalr          += (1/(1-γ)) * traj.r_t[traj.len]
            break
        end

        #Override the state:
        s_t1 .= s_t2;

    end

    #Return the average cost of the trajectory
    return totalr / (traj.len*2)
end

#Get the required action
@everywhere function sample_action(μ, σ, rand_vec)
    μ = convert(Array{Float64}, μ)
    return  μ .+ rand_vec .* σ
end

#Interact with environment
@everywhere function environment_call(s_t1::Array{Float64,1},a_t::Array{Float64,1}, params::Params)
    #Ineract with environment:
    s_t2 = rk4_func(s_t1, a_t, params.dt)

    #Claim the reward:
    r_t, terminal = cost_func(s_t2,a_t)

    return s_t2, r_t, terminal
end

#Functions focusing on the episode execution: Value and Actions combined
#-------------------------------------------------------------------------------
@everywhere function run_episode!(params::Params,θ::Any, traj::Trajectory, rms::RunMeanStd)
    dt         = params.dt;
    t_horizon  = params.t_horizon;
    s_init     = params.s_init + 0.1*sum(params.σ)*rand(length(params.s_init))
    l_traj     = params.l_traj;
    σ          = params.σ;
    u_range    = params.u_range;
    s_init[3]  = randn()

    #Initialize :
    s_t1 = copy(s_init)
    s_t2 = similar(s_t1)
    totalr = 0

    for i_t = 1:l_traj
        #Normalize the Input
        s_norm = normalize_rms(rms,s_t1)
        #s_norm = s_t1

        #Evaluate the Network Output
        out_t, val_t = mlp_nlayer(θ,s_norm,params)

        #Get Action with stochastic policy:
        loc_a_t = sample_action(out_t, params.σ)
        p_a_t   = get_π_prob(out_t,σ,loc_a_t)

        #Interact with the Environment:
        s_t2, loc_r_t, terminal = environment_call(s_t1,loc_a_t,params)

        #Save data in trajectory:
        traj.s_t[i_t,:]    .= s_norm
        traj.s_t2[i_t,:]   .= s_t2
        traj.a_t[i_t,:]    .= loc_a_t
        traj.p_a_t[i_t,:]   = p_a_t
        traj.r_t[i_t]       = loc_r_t*0.01
        traj.val_t[i_t]     = val_t
        traj.out_t[i_t,:]   = out_t
        traj.len            = i_t
        totalr             += loc_r_t

        #Override the state:
        s_t1 .= s_t2;

        #Check if the terminal condition is fullfilled and break
        if terminal == 1
                break
        end

    end

    return totalr / traj.len
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using Test
using Turing: BinomialLogit
using Distributions: Binomial, logpdf
using StatsFuns: logistic

n = 10
logitp = randn()
d1 = BinomialLogit(n, logitp)
d2 = Binomial(n, logistic(logitp))
k = 3
@test logpdf(d1, k) ≈ logpdf(d2, k)


text:
using BDF, Lexicon
include("extract_docstrings.jl")
#Lexicon.save("../docs/API.md", BDF)
extract_docstrings(["../src/BDF.jl"], "../docs/API.md")
cd("../")
run(`mkdocs build`)
cd("prep-release")


text:
using PyLCM
using Base.Test

# We need to add the current folder to the python path so it will pick up our message type
using PyCall
unshift!(PyVector(pyimport("sys")["path"]), ".")

test_data = rand(Int32, 2, 3, 4)
test_strings =  AbstractString["foo$(rand(Int32))" for i in 1:2, j in 1:size(test_data, 3)]

global got_message = false

@pyimport pylcm_test_message
msg = pylcm_test_message.multidimensional_array_t()
msg[:size_a] = size(test_data, 1)
msg[:size_b] = size(test_data, 2)
msg[:size_c] = size(test_data, 3)
msg[:data] = test_data
msg[:strarray] = test_strings

lc = LCM()

function handle_msg(channel, msg)
	global got_message = true
	for i = 1:size(test_data, 1)
		for j = 1:size(test_data, 2)
			for k = 1:size(test_data, 3)
				@test msg[:data][i,j][k] == test_data[i,j,k]
			end
		end
	end
	@test msg[:strarray] == test_strings
end

subscribe(lc, "TEST", handle_msg, pylcm_test_message.multidimensional_array_t)
publish(lc, "TEST", msg)
@test handle(lc)
@test got_message

global decoded_message = false

function decode_and_handle_msg(channel, msg_data)
	global decoded_message = true
	msg = pylcm_test_message.multidimensional_array_t[:decode](msg_data)
	for i = 1:size(test_data, 1)
		for j = 1:size(test_data, 2)
			for k = 1:size(test_data, 3)
				@test msg[:data][i,j][k] == test_data[i,j,k]
			end
		end
	end
	@test msg[:strarray] == test_strings
end

subscribe(lc, "TEST_DECODE", decode_and_handle_msg)
publish(lc, "TEST_DECODE", msg)
@test handle(lc, Dates.Second(1))
@test decoded_message


text:
"""
    taylor(df,a,b,y0,M)

berisi program untuk mencari solusi persamaan differensial `y'=f(t,y)`
dengan masalah nilai awal `y(a) = y0` pada interval `[a, b]`. Program ini secara default
berisi 5 masukan, yaitu fungsi `df` yang berisi nilai dari fungsi y', y'', y''', y'''',
titik ujung interval penyelesaian `[a,b]`, nilai awal
`y0`, dan jumlah sub-interval `M`.
# Examples
```jldoctest
julia> function df(t,y)
         y1 = (t-y)/2;
         y2 = (2-t+y)/4;
         y3 = (-2+t-y)/8;
         y4 = (2-t+y)/16;
         z = [y1,y2,y3,y4];
       end

julia> sol1 = taylor(df,0,3,1,6)
7×2 Array{Float64,2}:
 0.0  1.0
 0.5  0.836426
 1.0  0.819628
 1.5  0.917142
 2.0  1.10368
 2.5  1.35956
 3.0  1.66943
```
return solusi masalah nilai awal `sol`.
"""
function taylor(df,a,b,y0,M)
  M = Int(M)
  h = (b-a)/M;
  T = a:h:b;
  Y = Array{Float64}(undef,length(T),1)
  Y[1]=y0;
  for j = 1:M
    D = df(T[j],Y[j]);
    Y[j+1]=Y[j]+h*(D[1]+h*(D[2]/2+h*(D[3]/6+h*D[4]/24)));
  end
  sol = [T Y];
end


text:
using LinearAlgebra, SparseArrays, DiffEqOperators, Random, Test
################################################################################
# Test 2d extension
################################################################################

Random.seed!(7373)

#Create Array
n = 8
m = 15
A = rand(n,m)

#Create atomic BC
q1 = RobinBC((1.0, 2.0, 3.0), (0.0, -1.0, 2.0), 0.1, 4.0)
q2 = PeriodicBC(Float64)

BCx = vcat(fill(q1, div(m,2)), fill(q2, m-div(m,2)))  #The size of BCx has to be all size components *except* for x
BCy = vcat(fill(q1, div(n,2)), fill(q2, n-div(n,2)))

Qx = MultiDimBC{1}(BCx)
Qy = MultiDimBC{2}(BCy)

Ax = Qx*A
Ay = Qy*A

@test size(Ax)[1] == size(A)[1]+2
@test size(Ay)[2] == size(A)[2]+2

for j in 1:m
    @test Ax[:, j]  == Array(BCx[j]*A[:, j])
end
for i in 1:n
    @test Ay[i,:] == Array(BCy[i]*A[i,:])
end

################################################################################
# Test 3d extension
################################################################################

#Create Array
n = 8
m = 11
o = 12
A = rand(n,m, o)

#Create atomic BC
q1 = RobinBC((1.0, 2.0, 3.0), (0.0, -1.0, 2.0), 0.1, 4.0)
q2 = PeriodicBC(Float64)

BCx = vcat(fill(q1, (div(m,2), o)), fill(q2, (m-div(m,2), o)))  #The size of BCx has to be all size components *except* for x
BCy = vcat(fill(q1, (div(n,2), o)), fill(q2, (n-div(n,2), o)))
BCz = fill(Dirichlet0BC(Float64), (n,m))

Qx = MultiDimBC{1}(BCx)
Qy = MultiDimBC{2}(BCy)
Qz = MultiDimBC{3}(Dirichlet0BC(Float64), size(A)) #Test the other constructor

Ax = Qx*A
Ay = Qy*A
Az = Qz*A

Q = compose(Qx,Qy,Qz)
QL, Qb = Array(Q, size(A))
QLs, Qbs = sparse(Q, size(A))

A_conc = QL*reshape(A, prod(size(A))) .+ Qb
A_conc_sp = QLs*reshape(A,prod(size(A))) .+ Qbs

#test BC concretization
A_arr = Array(Q*A)
@test reshape(A_arr, prod(size(A_arr))) ≈ A_conc_sp ≈ A_conc

@test size(Ax)[1] == size(A)[1]+2
@test size(Ay)[2] == size(A)[2]+2
@test size(Az)[3] == size(A)[3]+2
for j in 1:m, k in 1:o
    @test Ax[:, j, k] == Array(BCx[j, k]*A[:, j, k])
end
for i in 1:n, k in 1:o
    @test Ay[i, :, k] == Array(BCy[i, k]*A[i, :, k])
end
for i in 1:n, j in 1:m
    @test Az[i, j, :] == Array(BCz[i, j]*A[i, j, :])
end

#test compositions to higher dimension
for N in 2:6
    sizes = rand(4:7, N)
    local A = rand(sizes...)

    Q1_N = RobinBC(Tuple(rand(3)), Tuple(rand(3)), fill(0.1, N), 4.0, size(A))

    local Q = compose(Q1_N...)

    A1_N = Q1_N.*fill(A, N)

    local A_arr = Array(Q*A)
    Q_l, Q_b = sparse(Q, size(A))

    @test A_arr ≈ Array(compose(A1_N...))
    @test A_arr ≈ reshape(Q_l*reshape(A, length(A)) .+ Q_b, size(A_arr)) #Test concretization
end


text:
struct Model
    fw
    fc
    fs
end

function Model(wordembeds::Matrix{T}, charembeds::Matrix{T}, ntags::Int) where {T}
    fw = Lookup(wordembeds)
    fc = @graph x begin
        x = Lookup(T,100,10)(x)
        x = Conv1D(T,50,50,20,10)(x)
        max(x, 2)
    end
    fs = @graph x begin
        x = Conv1D(T,750,300,300,150)(x)
        x = relu(x)
        Linear(T,300,ntags)(x)
    end
    Model(fw, fc, fs)
end

function (m::Model)(word::Var, char::Var)
    w = m.fw(word)
    c = m.fc(char)
    c.batchdims = w.batchdims
    s = cat(1, w, c)
    m.fs(s)
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
module TypedTables

using Unicode
using Tables
using SplitApplyCombine

using Base: @propagate_inbounds, @pure, OneTo, Fix2
import Tables.columns, Tables.rows

export @Compute, @Select
export Table, FlexTable, columns, rows, columnnames, showtable

# Resultant element type of given column arrays
@generated function _eltypes(a::NamedTuple{names, T}) where {names, T <: Tuple{Vararg{AbstractArray}}}
    Ts = []
    for V in T.parameters
        push!(Ts, eltype(V))
    end
    return NamedTuple{names, Tuple{Ts...}}
end

_ndims(a::NamedTuple{<:Any, T}) where {T} = _ndims(T)
_ndims(::Type{<:Tuple{Vararg{AbstractArray{<:Any, n}}}}) where {n} = n

# The following code causes newer versions of Julia to hang in precompilation
# and the workaround should not be needed after JuliaLang/julia#30577 was
# merged.
if VERSION < v"1.2.0-DEV.291"
    # Workaround for JuliaLang/julia#29970
    let
        for n in 1:32
            Ts = [Symbol("T$i") for i in 1:n]
            xs = [:(x[$i]) for i in 1:n]
            NT = :(Core.NamedTuple{names, Tuple{$(Ts...)}})
            eval(quote
                $NT(x::Tuple{$(Ts...)}) where {names, $(Ts...)} = $(Expr(:new, NT, xs...))
            end)
        end
    end
end

include("properties.jl")
include("Table.jl")
include("FlexTable.jl")
include("columnops.jl")
include("show.jl")

end # module


text:
module ECharts

	using JSON, Parameters, NoveltyColors, ColorBrewer, StatsBase, Blink
	using Dates, Random, REPL, DataFrames, WebIO, JSExpr
	using KernelDensity, Distributions
	import Base: print, show
	import LinearAlgebra: triu!, triu, tril, tril!

	# import JSON
	import JSON: StructuralContext, Serializations.CommonSerialization

	export print

	# Create base color library
	const colorpalettes = merge(ColorBrewer.colorSchemes, NoveltyColors.ColorDict)
	export colorpalettes

	export EChart
	export BoxPlotSeries, CandleStickSeries, EffectScatterSeries, FunnelSeries, GaugeSeries, GraphSeries,
	HeatmapSeries, LinesSeries, MapSeries, ParallelSeries, PictorialBarSeries, PieSeries, RadarSeries,
	SankeySeries, SunburstSeries, ThemeRiverSeries, TreeSeries, TreemapSeries, XYSeries
	export Title, Axis, Series, Toolbox, DataZoom
	export Tooltip, Legend, Grid, Timeline
	export LineStyle, AreaStyle, ItemStyle, TextStyle
	export AxisLine, AxisTick, AxisLabel, SplitLine, SplitArea, MarkLine, MarkArea, MarkPoint
	export JSFunction, Theme, VisualMap

	export xy_plot, bar, radialbar, line, scatter, area, waterfall
	export box, candlestick, sankey
	export radar, funnel
	export pie, donut, gauge, polar
	export histogram, heatmap, streamgraph, bubble
	export corrplot

	export title!, yaxis!, xaxis!, toolbox!, colorscheme!, flip!, seriesnames!, legend!, slider!, smooth!
	export yline!, xline!, lineargradient, radialgradient, text!, xarea!, yarea!, xgridlines!, ygridlines!
	export radial!, jitter!, labels!, theme!, tooltip!, aria!

	# This is a package local function, it is NOT overloading JSON.json
	# Define custom JSON serialization rule
	struct JSSerialization <: CommonSerialization end
	struct JSFunction
		data::String
	end

	function JSON.show_json(io::StructuralContext, ::JSSerialization, f::JSFunction)
		for line in split(f.data, '\n')
			Base.print(io, line)
		end
	end

	json(x) = sprint(JSON.show_json, JSSerialization(), x)
	# end custom JSON serialization rule

	# Primitives - in order of descending dependency within files
	include("theme.jl")
	include("definetypes.jl")
	include("render.jl")

	# Visualization mutating functions
	include("chartopts/utilities.jl")
	include("chartopts/seriesnames.jl")
	include("chartopts/xyaxis.jl")
	include("chartopts/gridlines.jl")
	include("chartopts/colorscheme.jl")
	include("chartopts/titletext.jl")
	include("chartopts/smooth.jl")
	include("chartopts/toolbox.jl")
	include("chartopts/tooltip.jl")
	include("chartopts/legend.jl")
	include("chartopts/xyline.jl")
	include("chartopts/xyarea.jl")
	include("chartopts/flip.jl")
	include("chartopts/slider.jl")
	include("chartopts/labels.jl")
	include("chartopts/jitter.jl")
	include("chartopts/theme.jl")
	include("chartopts/radial.jl")
	include("chartopts/gradients.jl")
	include("chartopts/aria.jl")

	# Plots
	include("plots/xy_plot.jl")
	include("plots/bar.jl")
	include("plots/radialbar.jl")
	include("plots/line.jl")
	include("plots/scatter.jl")
	include("plots/area.jl")
	include("plots/radar.jl")
	include("plots/circular_plot.jl")
	include("plots/pie.jl")
	include("plots/donut.jl")
	include("plots/gauge.jl")
	include("plots/polar.jl")
	include("plots/funnel.jl")
	include("plots/waterfall.jl")
	include("plots/box.jl")
	include("plots/candlestick.jl")
	include("plots/sankey.jl")
	include("plots/heatmap.jl")
	include("plots/histogram.jl")
	include("plots/streamgraph.jl")
	include("plots/bubble.jl")
	include("plots/corrplot.jl")

	# From JMW originally
	makevalidjson(s::Symbol) = string(s)
	makevalidjson(v::Vector) = [makevalidjson(v_i) for v_i in v]

	function makevalidjson(v::Dict)
	    res = Dict()
	    for (k, v) in v
	        res[k] = makevalidjson(v)
	    end
	    return res
	end
	makevalidjson(x::Missing) = "-" # this is a backstop to a value of missing getting to render step. "-" represents missing in ECharts
	makevalidjson(x::Any) = x

	# By convention, using single underscore at beginning to get around reserved words
	# Check for this, strip off _ before creating js
	function makevalidjson(x::AbstractEChartType)
	    res = Dict()

	    for f in fieldnames(typeof(x))
	        if getfield(x, f) != nothing
	        	startswith(string(f), "_") ? res[string(f)[2:end]] = makevalidjson(getfield(x, f)) : res[string(f)] = makevalidjson(getfield(x, f))
	        end
	    end
	    return res
	end

end # module


text:
#####
##### generic code
#####

"""
Abstract type for convenient dispatch.

**Not part of the API**, takes care of the following defaults:

1. `hypercube_dimension` falls back to `dimension`,

2. `logdensity` works through `logdensity_and_gradient` (inefficient, but should not matter)
"""
abstract type SamplingLogDensity end

"""
$(SIGNATURES)

Dimension for the second argument argument [`hypercube_transform`].

Falls back to `dimension` for `SamplingLogDensity`, since only mixtures need to override it.

Mostly for internal use, see [`samples`](@ref).
"""
hypercube_dimension(ℓ::SamplingLogDensity) = dimension(ℓ)

logdensity(ℓ::SamplingLogDensity, x) = first(logdensity_and_gradient(ℓ, x))

"""
$(FUNCTIONNAME)(ℓ, x)

Transform ``x ∈ [0,1]ⁿ`` into a random variable that has the given distribution when `x` is
uniform on the given hypercube, where `n` is `hypercube_dimension(ℓ)`.

The result needs to be distinct from `x`, as that is a buffer that may be reused.

Mostly for internal use, see [`samples`](@ref).
"""
function hypercube_transform end

"""
$(SIGNATURES)

`N` samples from `ℓ`, as columns of a matrix.
"""
function samples(ℓ::SamplingLogDensity, N::Integer)
    K = hypercube_dimension(ℓ)
    Z = Matrix{Float64}(undef, dimension(ℓ), N)
    s = SobolSeq(K)
    x = Vector{Float64}(undef, K)
    for i in 1:N
        next!(s, x)
        Z[:, i] = hypercube_transform(ℓ, x)
    end
    Z
end

capabilities(::Type{<:SamplingLogDensity}) = LogDensityOrder{1}()


text:
using ThreeBodyProblem
using Test

sys = System(ThreeBodyProblem.JUPITER,ThreeBodyProblem.CALLISTO)
@test typeof(sys.prim) == Body
@test typeof(sys.sec) == Body
@test typeof(sys) == System

# Create the ice planet Hoth
m = 1e24 # {kg} mass
R = 1e6 # {km} mean radius
a = 1.5*AU # {km} mean semimajor axis
T = 0.6*JY # {sec} mean orbital period
name = "Hoth"
newplanet = Body(m, R, a, T, name)
@test typeof(newplanet) == Body
@test newplanet.m == m
@test newplanet.R == R
@test newplanet.a == a
@test newplanet.T == T
@test newplanet.name == name

# Create the forest moon Endor
m = 1e22 # {kg} mass
R = 1e5 # {km} mean radius
a = 0.2*AU # {km} mean semimajor axis
T = 2π*sqrt(a^3/(G*newplanet.m)) # {sec} mean orbital period
newmoon = Body(m, R, a, T)
@test typeof(newmoon) == Body
@test newmoon.m == m
@test newmoon.R == R
@test newmoon.a == a
@test newmoon.T == T
@test newmoon.name == "NewPlanet"
@test newmoon.color == :blue

sys = System(newplanet, newmoon)
@test typeof(sys.prim) == Body
@test typeof(sys.sec) == Body
@test typeof(sys) == System

sys = sun_mercury()
@test sys.prim.name == "Sun"
@test sys.sec.name == "Mercury"
@test sys.name == "Sun/Mercury"

sys = sun_venus()
@test sys.name == "Sun/Venus"
sys = sun_earth()
@test sys.name == "Sun/Earth"
sys = earth_moon()
@test sys.name == "Earth/Moon"
sys = sun_mars()
@test sys.name == "Sun/Mars"
sys = sun_jupiter()
@test sys.name == "Sun/Jupiter"
sys = jupiter_europa()
@test sys.name == "Jupiter/Europa"
sys = sun_saturn()
@test sys.name == "Sun/Saturn"
sys = saturn_enceladus()
@test sys.name == "Saturn/Enceladus"
sys = sun_uranus()
@test sys.name == "Sun/Uranus"
sys = sun_neptune()
@test sys.name == "Sun/Neptune"


text:
using SoftRobots
using Base.Test

include("test_snake.jl")
# include("test_two_snakes.jl")


text:
# Usefull functions for subcarrier and carrier mappign in BF-OFDM

"""
---
Get the index of the oversampled carrier (i.e subcarriers index in the precoding field) based on the index of the allocated carriers. It is assumed that for each carriers, all subcarriers are enable (with orthogonality, i.e Nδ allocated subcarriers).
# --- Syntax
subcarriers	= getBFOFDM_subCarrierFromFBMCCarriers(nFBMC,nOFDM,δ,fbmcCarriers)
# --- Input parameters
- nFBMC	  : FBMC carrier size (number of carriers) [Int]
- nOFDM	  : OFDM subcarrier size [Int]
- δ		  : Rate factor [Float64]
- fbmcCarriers : Vector of allocated carrier [Array{Int,L}], L < nFBMC
# --- Output parameters
- subcarriers : Vector of allocated subcarriers [Array{Int, L δ ]
# v 1.0 - Robin Gerzaguet.
"""
function getBFOFDM_subCarrierFromFBMCCarriers(nFBMC, nOFDM, δ, fbmcCarriers)
	# ----------------------------------------------------
	# --- Getting parameters
	# ----------------------------------------------------
	nbDataCarrier	= length(fbmcCarriers);
	subcarrierPerCarrier  = Int(nOFDM*δ);
	# --- Index of occupied carriers in the oversampled plan
	carrierOversamp	      = zeros(Int,nbDataCarrier*subcarrierPerCarrier);
	for idx = 1:nbDataCarrier
		carrierOversamp[ collect(1+(idx-1)*Int(nOFDM*δ):idx*Int(nOFDM*δ)) ]  = 1 .+ mod.((fbmcCarriers[idx]-1)*Int(nOFDM*δ)-Int(nOFDM*δ/2)  .+  collect(1:Int(nOFDM*δ)) .-1,Int(nOFDM*nFBMC*δ));
		end
	return carrierOversamp
end

"""
---
Get the index of the allocated carriers based on the index of the subcarriers. We consider here that a FBMC carriers bear data (i.e is allocated) is at least one of its subcarrier is enable.
# --- Syntax
fbmcCarriers	= getBFOFDM_carrierFromSubcarriers(nFBMC,nOFDM,δ,subcarrier)
# --- Input parameters
- nFBMC	  : FBMC carrier size (number of carriers) [Int]
- nOFDM	  : OFDM subcarrier size [Int]
- δ		  : Rate factor [Float64]
- fbmcCarriers : Vector of allocated carrier [Array{Int,L}], L < nFBMC*nOFDM*δ
# --- Output parameters
- fbmcCarriers : Vector of allocated carriers [Array{Int, P ] P < nFBMC
# ---
# v 1.0 - Robin Gerzaguet.
"""
function getBFOFDM_carrierFromSubcarriers(nFBMC,nOFDM,δ,subcarriers)
	# ----------------------------------------------------
	# --- Getting parameters
	# ----------------------------------------------------
	# --- Allocate worst case scenario
	subcarrierPerCarrier =  nOFDM * δ;
	allCarriers = zeros(Int,nFBMC);		# --- Carrier buffer
	nbC         = 1;					# --- Carrier counter
	# --- Iterative loop
	for iN = 1 : 1 : nFBMC
		# --- For each carrier, we calculate associated "subcarrier index length
		# If at least 1 index of the subcarrier index is in subcarrier, we decide to allocate the index.
		# --- Calculation of subcarrier index for current FBMC carrier
		tmp	 = 1 .+ mod.( (iN-1)*Int(nOFDM*δ) .- Int(nOFDM*δ/2) .+ collect(1:nOFDM*δ) .-1,nOFDM*nFBMC*δ );
		# --- If at least one index belongs to subcarrier and tmp, store as an allocated carrier
		tmpInt	= intersect(tmp,subcarriers);
		if length(tmpInt) != 0
			# --- Savign currennt carrier index
			allCarriers[nbC] = iN;
			# --- Setting next element for saving next carrier
			nbC = nbC + 1;
		end
	end
	# --- Returnns onnly allocated vector
	return allCarriers[1:nbC-1];
end

"""
---
Returns the oversampled grid at the subcarrier level for BF-OFDM. In BF-OFDM as we use a nOFDM*nFBMC grid with nOFDM*δ allocated subcarrier per carriers (for a precodinng stage of size nOFDM), it can be usefull to get all the index of the allocated subcarrier per carrier. The output is a vector of the allocated subcarrier index in the oversampled frequency grid.
# --- Syntax
subcarrierTx	= getCarrierFromSubcarriers(nOFDM,δ,fbmcCarriers)
# --- Input parameters
- nOFDM	  : OFDM subcarrier size [Int]
- δ		  : Rate factor [Float64]
- fbmcCarriers : Vector of allocated carrier [Array{Int,L}], L < nFBMC*nOFDM*δ
# --- Output parameters
- subcarrierTx : Vector of allocated carriers [Array{Int, L*nOFDM*δ  ]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function getBFOFDM_oversampledGridSubcarriers(nFBMC,nOFDM,δ,fbmcCarriers)
	# --- Getting alloc. parameters
	nbDataCarrier	= length(fbmcCarriers);
	subcarrierPerCarrier  = Int(nOFDM*δ);
	# --- Index of occupied carriers in the oversampled plan
	subcarrierTx          = zeros(Int,nbDataCarrier*subcarrierPerCarrier);
	for idx = 1:nbDataCarrier
		# --- Calculating index of subcarriers index into carriers
		# It depends on delta (delta mapping pattern exists)
		# We also have a shift which is unknown
		# --- Specific pattern index for given carrier
		tmp  = mod.((-Int(nOFDM*δ/2)+(fbmcCarriers[idx]-1)*subcarrierPerCarrier: Int(nOFDM*δ/2) +(fbmcCarriers[idx]-1)*subcarrierPerCarrier-1), nOFDM) .+1;
		# --- Store this pattern as function of carrier index
		# No need to have a absolute positionning as after that we will use
		# only matrix with data
		subcarrierTx[collect(1 + (idx-1)*Int(nOFDM*δ):idx*Int(nOFDM*δ)) ] = tmp .+ nOFDM*(idx-1);
	end
	return subcarrierTx
end


text:
#
# tools.jl --
#
# Tools and utilities for scientific cameras.
#
#------------------------------------------------------------------------------
#
# This file is part of the `ScientificCameras.jl` package which is licensed
# under the MIT "Expat" License.
#
# Copyright (C) 2017, Éric Thiébaut.
#

"""
    mean([T = Float64,] cam, num; kwds...) -> img

yields the mean image `img` from the sample averaging of `num` images acquired
with camera `cam`.  Optional argument `T` is the floating-point type of the
elements of the result.

See also: [`stat`](@ref).

"""
mean(cam::ScientificCamera, args...; kwds...) =
    mean(Float64, cam, args...; kwds...)

function mean(::Type{T},
              cam::ScientificCamera,
              num::Integer = 100;
              kwds...) where {T <: AbstractFloat}
    roi = getroi(cam)
    dims = (roi.width, roi.height)
    sum, cnt = processimages(cam, num,
                             (sum, img, ticks, count) -> (sum .+= img; sum),
                             zeros(T, dims); kwds...)

    return convert(T, 1/cnt)*sum
end

"""
    stat([T = Float64,] cam, num; kwds...) -> avg, std, cnt

yields the mean image `avg` and standard deviation `std` of images from the
sample averaging of `num` images acquired with camera `cam`.  Optional argument
`T` is the floating-point type of the elements of the result.  The returned
value `cnt` is the actual number of samples and may be less than `num` if a
timeout occured and keyword `truncate` is true.

See also: [`stat`](@ref).

"""
stat(cam::ScientificCamera, args...; kwds...) =
    stat(Float64, cam, args...; kwds...)

function stat(::Type{T},
              cam::ScientificCamera,
              num::Integer = 100;
              kwds...) where {T <: AbstractFloat}
    # Integrate information.
    @assert num ≥ 2 "number of images must be at least 2"
    roi = getroi(cam)
    dims = (roi.width, roi.height)
    tup, cnt = processimages(cam, num,
                             (stat, img, ticks, count) -> updatestatistics!(stat, img),
                             (zeros(T, dims), zeros(T, dims)); kwds...)
    @assert cnt ≥ 2 "number of images must be at least 2"

    # Post-processing.
    s1, s2 = tup
    c1 = convert(T, 1/cnt)
    c2 = convert(T, 1/(cnt - 1))
    @inbounds @simd for i in 1:length(s1)
        v1 = convert(T, s1[i])
        v2 = c1*v1
        v3 = s2[i] - v2*v1
        s1[i] = v2
        s2[i] = (v3 > zero(T) ? sqrt(c2*v3) : zero(T))
    end

    return s1, s2, convert(Int, cnt)
end

function updatestatistics!(state::DenseArray{Ts,Ns},
                           arr::DenseArray{Ti,Ni}) where {Ts <: AbstractFloat,
                                                          Ns, Ti <: Real, Ni}
    @assert Ns == Ni + 1
    @assert size(state,1) == 2
    @assert size(state)[2:end] == size(arr)
    @inbounds @simd for i in 1:length(arr)
        val = convert(Ts, arr[i])
        state[1,i] += val
        state[2,i] += val*val
    end
    return state
end

function updatestatistics!(state::NTuple{2,DenseArray{Ts,N}},
                           arr::DenseArray{Ti,N}) where {Ts <: AbstractFloat,
                                                         Ti <: Real, N}
    sum1, sum2 = state
    @assert size(sum1) == size(sum2) == size(arr)
    @inbounds @simd for i in 1:length(arr)
        val = convert(Ts, arr[i])
        sum1[i] += val
        sum2[i] += val*val
    end
    return state
end

"""
    processimages(cam, num, proc, state;
                  skip=0, timeout=sec, truncate=false) -> state, count

processes `num` images from camera `cam` by calling the function `proc` to
process each image and update `state`:

    state = proc(state, img, ticks, count)

the final state and the actual number of processed images are returned.

"""
function processimages(cam::ScientificCamera,
                       number::Integer,
                       process::Function,
                       state;
                       skip::Integer = 0,
                       timeout::Real = defaulttimeout(cam),
                       truncate::Bool = false)
    # Check arguments.
    number ≥ 1 || throw(ArgumentError("invalid number of images to process"))
    skip ≥ 0 || throw(ArgumentError("invalid number of images to skip"))
    timeout > zero(timeout) || error("invalid timeout")

    # Acquire and process a sequence of images.
    count = zero(number)
    start(cam, 4)
    while count < number
        try
            img, ticks = wait(cam, timeout)
            if skip > zero(skip)
                # Skip this image.
                skip -= one(skip)
            else
                # Process this image.
                count += one(count)
                state = process(state, img, ticks, count)
            end
            release(cam)
        catch err
            if truncate && isa(err, TimeoutError)
                quiet || warn("Acquisition timeout after $cnt image(s)")
                number = count
            else
                abort(cam)
                rethrow(err)
            end
        end
    end
    abort(cam)
    return state, count
end


text:
using DataDeps

for (file, checksum) in [("train", "bcbbe17c487d0939d48c2d694622303edb3637ca9c4944776628cd1815c5cb34"),
                        ("test", "2695b931a505b65b1d3fe2a5aa0d749edcd53d786278da6d7726714e1426fbe8")]
    register(DataDep("CoNLL 2000 $file",
        """
        Website: https://www.clips.uantwerpen.be/conll2000/chunking/$file.txt.gz
        Data Website:https://www.clips.uantwerpen.be/conll2000/chunking/

        The Fourth Conference on Natural Language Learning (CoNLL-2003) was held on September 13 and September 14, 2000 in conjunction with ICGI-2000 and LLL-2000 at the Instituto Superior Técnico in Lisbon, Portugal.

        The details of the task are in the following paper, please cite the following publication if you are using the corpora.

        Erik F. Tjong Kim Sang, Sabine Buchholz. "Introduction to the CoNLL-2000 Shared Task: Chunking." In: Proceedings of CoNLL-2000 and LLL-2000, pages 127-132, Lisbon, Portugal, 2000.
        https://www.clips.uantwerpen.be/conll2000/pdf/12732tjo.pdf
        """,
        "https://www.clips.uantwerpen.be/conll2000/chunking/$file.txt.gz",
        checksum;
        post_fetch_method = unpack
    ))
end


text:
data = parse.(Int32, readlines("day01.txt"))

function part1(data::Array{Int32,1})
    l = length(data)
    for i in 1:(l - 1), j in (i + 1):l
        if data[i] + data[j] == 2020
            return data[i] * data[j]
        end
    end
end

function part2(data::Array{Int32,1})
    l = length(data)
    for i in 1:(l - 2), j in (i + 1):(l - 1), k in (i + 2):l
        if data[i] + data[j] + data[k] == 2020
            return data[i] * data[j] * data[k]
        end
    end
end

println("Part 1 solution $(part1(data))")
println("Part 2 solution $(part2(data))")


text:
# # Block Basics
# ## Table of Contents
# * Construction and Matrix Representation
# * Block Tree Architecture
# * Tagging System
# * Parameter System
# * Differentiable Blocks
# * Time Evolution and Hamiltonian

using Yao, Yao.Blocks
using LinearAlgebra

# ## Construction and Matrix Representation
# Blocks are operations on registers, we call those with matrix representation (linear) `MatrixBlock`.

# A `MatrixBlock` can be
# * isunitary, $O^\dagger O=I$
# * ishermitian, $O^\dagger = O$
# * isreflexive, $O^2 = 1$

@show X
@show X |> typeof
@show isunitary(X)
@show ishermitian(X)
@show isreflexive(X);

# **matrix representation**
mat(X)

# **composite gates**
# Embeding an X gate into larger Hilbert space, the first parameter of most non-primitive constructors are always qubit numbers

@show X2 = put(3, 2=>X)
@show isunitary(X2)
@show ishermitian(X2)
@show isreflexive(X2);

#---------------
mat(X2)
#---------------
@show cx = control(3, 3, 1=>X)
@show isunitary(cx)
@show ishermitian(cx)
@show isreflexive(cx);
#---------------
mat(cx)
# hermitian and reflexive blocks can be used to construct rotation gates
@show rx = rot(X, π/4)
@show isunitary(rx)
@show ishermitian(rx)
@show isreflexive(rx);

#-----------------
mat(rx)

# now let's build a random circuit for following demos
using Yao.Intrinsics: rand_unitary
circuit = chain(5, control(5, 3=>Rx(0.25π)), put(5, (2,3)=>matrixgate(rand_unitary(4))), swap(5, 3, 4), repeat(5, H, 2:5), put(5, 2=>Ry(0.6)))
# to apply it on some register, we can use
reg = zero_state(10)
focus!(reg, 1:5) do reg_focused
    apply!(reg_focused, circuit)
end
@show reg ≈ zero_state(10);   # reg is changed!
# then we reverse the process and check the correctness
focus!(reg, 1:5) do reg_focused
    reg_focused |> circuit'
end
@show reg ≈ zero_state(10);   # reg is restored!
# Here, we have used the pip "eye candy" `reg |> block` to represent applying a block on register, which is equivalent to `apply!(reg, block)`
#---------------

# **Type Tree**
# To see a full list of block types
using InteractiveUtils: subtypes
function subtypetree(t, level=1, indent=4)
   level == 1 && println(t)
   for s in subtypes(t)
     println(join(fill(" ", level * indent)) * string(s))
     subtypetree(s, level+1, indent)
   end
end

subtypetree(Yao.Blocks.AbstractBlock);

# In the top level, we have
# * `MatrixBlock`, linear operators
# * `AbstractMeasure`, measurement operations
# * `FunctionBlock`, a wrapper for register function that take register as input, change the register inplace and return the register.
# * `Sequential`, a container for **block tree**, which is similar to `ChainBlock`, but has less constraints.

# ## Block Tree Architecture
# A block tree is specified the following two APIs
# * subblocks(block), siblings of a block.
# * chsubblocks, change siblings of a node.
crx = circuit[1]
@show crx
@show subblocks(crx)
@show chsubblocks(crx, (Y,));

# if we want to define a function that travals over the tree in depth first order, we can write something like
function print_block_tree(root, depth=0)
    println("  "^depth * "- $(typeof(root).name)")
    print_block_tree.(root |> subblocks, depth+1)
end
print_block_tree(circuit);

# there are some functions defined using this strategy, like `collect(circuit, block_type)`, it can filter out any type of blocks
rg = collect(circuit, RotationGate)

# ## Tagging System
# We proudly introduced our tag system here.
# In previous sections, we have introduced the magic operation `circuit'` to get the dagger a circuit, its realization is closely related to the tagging mechanism of `Yao`.
@show X'    # hermitian gate
@show Pu'   # special gate
@show Rx(0.5)';   # rotation gate

# The dagger of above gates can be translated to other gates easily.
# but some blocks has no predefined dagger operations, then we put a tag for it as a default behavior, e.g.
daggered_gate = matrixgate(randn(4, 4))'
@show daggered_gate |> typeof
daggered_gate
# Here, `Daggered` is a subtype of `TagBlock`.
#--------------------
# Other tag blocks include
#-----------------
# **Scale**, static scaling
2X
# **CachedBlock**, get the matrix representation of a block when applying it on registers, and cache it in memory (or `CacheServer` more precisely).
# This matrix can be useful in future calculation, like boosting time evolution.
put(5, 2=>X) |> cache
# **AbstactDiff**, marks a block as differentiable, either in classical `back propagation` mode (with extra memory cost to store intermediate data)
put(5, 2=>Rx(0.3)) |> autodiff(:BP)
# or non-cheating quantum circuit simulation
put(5, 2=>Rx(0.3)) |> autodiff(:QC)

# ## Parameter System
# using the depth first searching strategy, we can find all parameters in a tree or subtree. Two relevant APIs are
# * parameters(block), get all parameters in a (sub)tree rooted on block
# * dispatch!([func], block, params), dispatch `params` into (sub)tree rooted on `block`, optional parameter `func` can be used to custom parameter update rule.
@show parameters(circuit)
dispatch!(circuit, [0.1, 0.9])
@show parameters(circuit)
dispatch!(+, circuit, [0.1, 0.1])
@show parameters(circuit)
dispatch!(circuit, :zero)
@show parameters(circuit)
dispatch!(circuit, :random)
@show parameters(circuit);

# ### Intrinsic parameters
# Intrinsic parameters are block's net contribution to total paramters, normally, we define these two APIs for subtyping blocks
# * `iparameters(block)`,
# * `setiparameters!(block, params...)`,
@show iparameters(Rx(0.3))
@show setiparameters!(Rx(0.3), 1.2)
@show chain(Rx(0.3), Ry(0.5)) |> iparameters;

# ## Differentiable Blocks
# see the independant chapter [Automatic Differentiation](@ref autodiff)
#-----------------
# ## Time Evolution and Hamiltonian
# docs are under preparation


text:
module GFIlogisticRegression

import Polyhedra
using CDDLib
import LinearAlgebra
import Distributions
import Optim
import StatsBase
import StatsModels
import DataFrames
import DataFramesMeta

export fidSummary
export fidConfInt
export fidQuantile
export fidProb
export fidSampleLR

function logit(u) # = qlogis
  return log(u / (1.0 - u))
end

function expit(x) # = plogis
  return 1.0 / (1.0 + exp(-x))
end

function rtlogis1(x)
  b = expit(x)
  if b == 0.0
    return x
  end
  return logit(b * rand())
end

function rtlogis2(x)
  a = expit(x)
  if a == 1.0
    return x
  end
  return logit(a + (1.0-a) * rand())
end

function from01(u)
  return map(log, u ./ (1.0 .- u))
end

function from01scalar(u)
  return log(u / (1.0 - u))
end

function to01(x)
  return 1.0 / (1.0 + exp(-x))
end

function dfrom01(u)
  return 1.0 / (u * (1.0 - u))
end

function dlogis(x)
  expminusx = map(exp, -x)
  oneplusexpminusx = 1.0 .+ expminusx
  return expminusx ./ (oneplusexpminusx .* oneplusexpminusx)
end

function ldlogis(x)
  return -x - 2.0 * map(log1p, map(exp, -x))
end

function dldlogis(x)
  return 1.0 .- 2.0 ./ (1.0 .+ map(exp, -x))
end

function forig(x, P, b)
  return prod(dlogis(P * x + b))
end

function f(u, P, b)
  return prod(dlogis(P * from01(u) + b))
end

function logf(u, P, b)
  return sum(ldlogis(P * from01(u) + b))
end

function df(ui, Pi, y1, y2)
  return y1 * dfrom01(ui) * sum(Pi .* y2)
end

function dlogf(ui, Pi, y2)
  return dfrom01(ui) * sum(Pi .* y2)
end

function get_umax(P, b)
  d = size(P, 2)
  fn = function(u)
    return -logf(u, P, b)
  end
  grfn! = function(storage, u)
    y2 = dldlogis(P * from01(u) + b)
    for i in 1:d
      storage[i] = -dlogf(u[i], P[:, i], y2)
    end
  end
  lower = sqrt(eps()) * ones(d)
  upper = 1.0 .- lower
  init = 0.5 * ones(d)
  opt = Optim.optimize(
    fn, grfn!, lower, upper, init, Optim.Fminbox(Optim.LBFGS())
  )
  return (
    mu = from01(opt.minimizer),
    umax = exp(-opt.minimum)^(2.0 / (d + 2))
  )
end

function get_vmin_i(P, b, j, mu)
  d = size(P, 2)
  alpha = 1.0 / (d + 2)
  fn = function(u)
    return f(u, P, b)^alpha * (from01scalar(u[j]) - mu[j])
  end
  grfn! = function(storage, u)
    y1alpha = f(u, P, b)^alpha
    y2 = dldlogis(P * from01(u) + b)
    diff = from01scalar(u[j]) - mu[j]
    storage[j] = y1alpha * dfrom01(u[j]) *
      (alpha * sum(P[:, j] .* y2) * diff + 1.0)
    others = deleteat!(collect(1:d), j)
    for i in others
      storage[i] = y1alpha * dfrom01(u[i]) * alpha * sum(P[:, i] .* y2) * diff
    end
  end
  eta = sqrt(eps())
  lower = eta * ones(d)
  upper = (1.0-eta) * ones(d)
  upper[j] = to01(mu[j])
  init = 0.5 * ones(d)
  init[j] = to01(mu[j]) / 2.0
  opt = Optim.optimize(
    fn, grfn!, lower, upper, init, Optim.Fminbox(Optim.LBFGS())
  )
  return opt.minimum
end

function get_vmin(P, b, mu)
  d = size(P, 2)
  vmin = Vector{Float64}(undef, d)
  for i in 1:d
    vmin[i] = get_vmin_i(P, b, i, mu)
  end
  return vmin
end

function get_vmax_i(P, b, j, mu)
  d = size(P, 2)
  alpha = 1.0 / (d + 2)
  fn = function(u)
    return -f(u, P, b)^alpha * (from01scalar(u[j]) - mu[j])
  end
  grfn! = function(storage, u)
    y1alpha = f(u, P, b)^alpha
    y2 = dldlogis(P * from01(u) + b)
    diff = from01scalar(u[j]) - mu[j]
    storage[j] = -y1alpha * dfrom01(u[j]) *
      (alpha * sum(P[:, j] .* y2) * diff + 1.0)
    others = deleteat!(collect(1:d), j)
    for i in others
      storage[i] = -y1alpha * dfrom01(u[i]) * alpha * sum(P[:, i] .* y2) * diff
    end
  end
  eta = sqrt(eps())
  lower = eta * ones(d)
  lower[j] = to01(mu[j])
  upper = (1.0-eta) * ones(d)
  init = 0.5 * ones(d)
  init[j] = (to01(mu[j]) + 1.0) / 2.0
  opt = Optim.optimize(
    fn, grfn!, lower, upper, init, Optim.Fminbox(Optim.LBFGS())
  )
  return -opt.minimum
end

function get_vmax(P, b, mu)
  d = size(P, 2)
  vmax = Vector{Float64}(undef, d)
  for i in 1:d
    vmax[i] = get_vmax_i(P, b, i, mu)
  end
  return vmax
end

function get_bounds(P, b)
  (mu, umax) = get_umax(P, b)
  return (
    mu = mu,
    umax = umax,
    vmin = get_vmin(P, b, mu),
    vmax = get_vmax(P, b, mu)
  )
end

function rcd(n, P, b)
  d = size(P, 2)
  sims = Vector{Vector{Float64}}(undef, n)
  (mu, umax, vmin, vmax) = get_bounds(P, b)
  if any(vmin .>= vmax)
    error("vmin .>= vmax")
  end
  alpha = 2.0 / (d+2)
  k = 0
  while k < n
    u = umax * rand()
    v = vmin + (vmax - vmin) .* rand(d)
    x = v / sqrt(u) + mu
    if u < forig(x, P, b)^alpha
      k += 1
      sims[k] = x
    end
  end
  return sims
end

"""
    fidSampleLR(formula, data, N[, gmp][, thresh])

Fiducial sampling of the parameters of the logistic regression model.

# Arguments
- `formula`: a formula describing the model
- `data`: data frame in which the variables of the model can be found
- `N`: number of simulations
- `gmp`: whether to use exact arithmetic in the algorithm
- `thresh`: the threshold used in the sequential sampler; the default `N/2` should not be changed

# Example

    using GFIlogisticRegression, DataFrames, StatsModels
    data = DataFrame(
      y = [0, 0, 1, 1, 1],
      x = [-2, -1, 0, 1, 2]
    )
    fidsamples = fidSampleLR(@formula(y ~ x), data, 3000)
"""
function fidSampleLR(formula, data, N, gmp = false, thresh = N/2)
  y = StatsModels.response(formula, data)
  mf = StatsModels.ModelFrame(formula, data)
  mm = StatsModels.ModelMatrix(mf)
  X = mm.m
  (_, colnames) = StatsModels.coefnames(mf.f)
  (n, p) = size(X)
  weight = ones(n, N)
  local WTnorm
  ESS = N .* ones(n)
  T = gmp ? Rational{BigInt} : Float64
  cdd = gmp ? :exact : :float
  CC = Vector{Array{T,2}}(undef, N)
  cc = Vector{Vector{T}}(undef, N)
  # Kstart ####
  Kstart = [1]
  i = 1
  rk = 1
  while rk < p
    i += 1
    Kstart_plus_i = vcat(Kstart, i)
    if LinearAlgebra.rank(X[Kstart_plus_i, :]) == rk + 1
      Kstart = Kstart_plus_i
      rk += 1
    end
  end
  @inbounds Xstart = X[Kstart, :]
  qXstart = convert(Array{T,2}, Xstart)
  ystart = y[Kstart]
  K = setdiff(1:n, Kstart)
  @inbounds XK = X[K, :]
  qXK = convert(Array{T,2}, XK)
  @inbounds yK = y[K]
  # t = 1 to p ####
  At = Array{Float64}(undef, p, N)
  for i in 1:N
    a = map(logit, rand(p))
    @inbounds At[:, i] = a
    C = Array{T,2}(undef, p, p)
    c = Vector{T}(undef, p)
    for j in 1:p
      @inbounds if ystart[j] == 0
        @inbounds C[j, :] = qXstart[j, :]
        @inbounds c[j] = convert(T, a[j])
      else
        @inbounds C[j, :] = -qXstart[j, :]
        @inbounds c[j] = convert(T, -a[j])
      end
    end
    @inbounds CC[i] = C
    @inbounds cc[i] = c
  end
  # t from p+1 to n ####
  for t in 1:(n-p)
    #println(t)
    At = vcat(At, Array{Float64,2}(undef, 1, N))
    @inbounds qXt = qXK[t, :]
    qXt_row = reshape(qXt, 1, :)
    qXtt = LinearAlgebra.transpose(qXt)
    @inbounds if yK[t] == 0
      for i in 1:N
        @inbounds H = Polyhedra.hrep(CC[i], cc[i])
        plyhdrn = Polyhedra.polyhedron(H, CDDLib.Library(cdd))
        pts = collect(Polyhedra.points(plyhdrn))
        MIN = convert(
          Float64, minimum(qXtt * hcat(pts...))
        )
        atilde = rtlogis2(MIN)
        @inbounds weight[t, i] = 1 - expit(MIN)
        @inbounds CC[i] = vcat(CC[i], qXt_row)
        @inbounds cc[i] = vcat(cc[i], convert(T, atilde))
        @inbounds At[p+t, i] = atilde
      end
    else
      for i in 1:N
        @inbounds H = Polyhedra.hrep(CC[i], cc[i])
        plyhdrn = Polyhedra.polyhedron(H, CDDLib.Library(cdd))
        pts = collect(Polyhedra.points(plyhdrn))
        MAX = convert(
          Float64, maximum(qXtt * hcat(pts...))
        )
        atilde = rtlogis1(MAX)
        @inbounds weight[t, i] = expit(MAX)
        @inbounds CC[i] = vcat(CC[i], -qXt_row)
        @inbounds cc[i] = vcat(cc[i], convert(T, -atilde))
        @inbounds At[p+t, i] = atilde
      end
    end
    WT = prod(weight; dims = 1)[1, :]
    WTnorm = WT ./ sum(WT)
    @inbounds ESS[p+t] = 1.0 / sum(WTnorm .* WTnorm)
    @inbounds if ESS[p+t] < thresh || t == n-p
      #println("alteration...")
      Nsons = rand(Distributions.Multinomial(N, WTnorm))
      counter = 1
      At_new = Array{Float64, 2}(undef, p+t, 0)
      @inbounds D = vcat(Xstart, XK[1:t, :])
      QR = LinearAlgebra.qr(D)
      @inbounds P = QR.Q[:, 1:p]#orth(D)
      Pt = LinearAlgebra.transpose(P)
      QQt = LinearAlgebra.I - P*Pt
      M = inv(LinearAlgebra.transpose(D)*D) * LinearAlgebra.transpose(D) * P
      CCtemp = Vector{Array{T,2}}(undef, N)
      cctemp = Vector{Vector{T}}(undef, N)
      for i in 1:N
        @inbounds ncopies = Nsons[i]
        if ncopies >= 1
          @inbounds CCtemp[counter] = CC[i]
          @inbounds cctemp[counter] = cc[i]
          @inbounds At_new = hcat(At_new, At[:, i])
          if ncopies > 1
            @inbounds H = Polyhedra.hrep(CC[i], cc[i])
            plyhdrn = Polyhedra.polyhedron(H, CDDLib.Library(cdd))
            pts = collect(Polyhedra.points(plyhdrn))
            lns = collect(Polyhedra.lines(plyhdrn))
            rys = collect(Polyhedra.rays(plyhdrn))
            @inbounds b = QQt * At[:, i]
            @inbounds B = Pt * At[:, i]
            BTILDES = rcd(ncopies-1, P, b)
            for j in 2:ncopies
              VT_new = Vector{Vector{T}}(undef, length(pts))
              @inbounds Btilde = BTILDES[j-1]
              At_tilde = P * Btilde .+ b
              At_new = hcat(At_new, At_tilde)
              for k in 1:length(pts)
                pt = convert(Vector{Float64}, pts[k])
                VT_new[k] =
                  convert(Vector{T}, pt .- M * (B .- Btilde))
              end
              V = Polyhedra.vrep(VT_new, lns, rys)
              plyhdrn = Polyhedra.polyhedron(V, CDDLib.Library(cdd))
              #HlfSpcs = Polyhedra.hrep(plyhdrn).halfspaces
              #A = map(x -> x.a, HlfSpcs)
              H = Polyhedra.hrep(plyhdrn)
              Hmat = convert(Array{T,2}, CDDLib.extractA(unsafe_load(H.matrix)))
              @inbounds CCtemp[counter + j - 1] = -Hmat[:, 2:end] #LinearAlgebra.transpose(hcat(A...))
              @inbounds cctemp[counter + j - 1] = Hmat[:, 1]#map(x -> x.β, HlfSpcs)
            end
          end
          counter += ncopies
        end
      end
      CC = CCtemp
      cc = cctemp
      At = At_new
      if t < n-p
        weight = ones(n, N)
      end
    end
  end
  Beta = Array{Float64, 2}(undef, N, p)
  for i in 1:N
    @inbounds H = Polyhedra.hrep(CC[i], cc[i])
    plyhdr = Polyhedra.polyhedron(H, CDDLib.Library(cdd))
    vertices = collect(Polyhedra.points(plyhdr))
    pts = hcat(vertices...)
    for j in 1:p
      pt = convert(Vector{Float64}, pts[j, :])
      if rand() < 0.5
        @inbounds Beta[i, j] = minimum(pt)
      else
        @inbounds Beta[i, j] = maximum(pt)
      end
    end
  end
  dfBeta = DataFrames.DataFrame(Beta, Symbol.(colnames))
  return (Beta = dfBeta, Weights = WTnorm)
end # fidSampleLR

"""
    fidSummary(fidsamples)

Summary of the fiducial simulations.

# Argument
- `fidsamples`: an output of `fidSampleLR`

# Example

    using GFIlogisticRegression, DataFrames, StatsModels
    data = DataFrame(
      y = [0, 0, 1, 1, 1],
      x = [-2, -1, 0, 1, 2]
    )
    fidsamples = fidSampleLR(@formula(y ~ x), data, 3000)
    fidSummary(fidsamples)
"""
function fidSummary(fidsamples)
  (_, p) = size(fidsamples.Beta)
  gdf = DataFrames.groupby(DataFrames.stack(fidsamples.Beta, 1:p), :variable)
  wghts = fidsamples.Weights
  function fsummary(x)
    return (
      mean = sum(wghts .* x),
      median = StatsBase.quantile(x, StatsBase.weights(wghts), 0.5),
      lwr = StatsBase.quantile(x, StatsBase.weights(wghts), 0.025),
      upr = StatsBase.quantile(x, StatsBase.weights(wghts), 0.975)
    )
  end
  return DataFrames.combine(gdf, DataFrames.valuecols(gdf) .=> fsummary => DataFrames.AsTable)
end

"""
    fidConfInt(parameter, fidsamples, conf)

Fiducial confidence interval of a parameter of interest.

# Arguments
- `parameter`: an expression of the parameter of interest given as a string; see the example
- `fidsamples`: an output of `fidSampleLR`
- `conf`: confidence level

# Example

    using GFIlogisticRegression, DataFrames, StatsModels
    data = DataFrame(
      y = [0, 0, 1, 1, 1, 1],
      group = ["A", "A", "A", "B", "B", "B"]
    )
    fidsamples = fidSampleLR(@formula(y ~ 0 + group), data, 3000)
    fidConfInt(":\\\"group: A\\\" - :\\\"group: B\\\"", fidsamples, 0.95)
"""
function fidConfInt(parameter, fidsamples, conf = 0.95)
  x = eval(:(DataFramesMeta.@with($(fidsamples.Beta), $(Meta.parse(parameter)))))
  wghts = fidsamples.Weights
  halpha = (1.0 - conf) / 2.0
  qntls = StatsBase.quantile(x, StatsBase.weights(wghts), [halpha; 1.0-halpha])
  return (lower = qntls[1], upper = qntls[2])
end

"""
    fidQuantile(parameter, fidsamples, p)

Fiducial quantile of a parameter of interest.

# Arguments
- `parameter`: an expression of the parameter of interest given as a string; see the example
- `fidsamples`: an output of `fidSampleLR`
- `p`: quantile level, between 0 and 1

# Example

    using GFIlogisticRegression, DataFrames, StatsModels
    data = DataFrame(
      y = [0, 0, 1, 1, 1, 1],
      group = ["A", "A", "A", "B", "B", "B"]
    )
    fidsamples = fidSampleLR(@formula(y ~ 0 + group), data, 3000)
    fidQuantile(":\\\"group: A\\\" ./ :\\\"group: B\\\"", fidsamples, 0.5)
"""
function fidQuantile(parameter, fidsamples, p)
  x = eval(:(DataFramesMeta.@with($(fidsamples.Beta), $(Meta.parse(parameter)))))
  wghts = fidsamples.Weights
  return StatsBase.quantile(x, StatsBase.weights(wghts), p)
end

"""
    fidProb(parameter, fidsamples, q)

Fiducial non-exceedance probability of a parameter of interest.

# Arguments
- `parameter`: an expression of the parameter of interest given as a string; see the example
- `fidsamples`: fiducial simulations, an output of `fidSampleLR`
- `q`: the non-exceedance threshold

# Example

    using GFIlogisticRegression, DataFrames, StatsModels
    data = DataFrame(
      y = [0, 0, 1, 1, 1],
      x = [-2, -1, 0, 1, 2]
    )
    fidsamples = fidSampleLR(@formula(y ~ x), data, 3000)
    fidProb("map(exp, :x)", fidsamples, 1) # this is Pr(exp(x) <= 1)
"""
function fidProb(parameter, fidsamples, q)
  x = eval(:(DataFramesMeta.@with($(fidsamples.Beta), $(Meta.parse(parameter)))))
  wghts = fidsamples.Weights
  return sum(wghts[findall(x .<= q)])
end

end # module

#=
P = [
  -0.82 -0.18;
  -0.41 -0.37;
  0.0   -0.55;
  0.41  -0.73
]
b = [-0.50; 1.58; -1.66; 0.58]
=#

#=
y = [0, 0, 1, 1, 1]
X = [1 -2; 1 -1; 1 0; 1 1; 1 2]
fidsamples = fidSampleLR(y, X, 5)
=#

#=
using DataFrames, StatsModels
data = DataFrame(
  y = [0, 0, 1, 1, 1],
  x = [-2, -1, 0, 1, 2]
)
fidsamples = fidSampleLR(@formula(y ~ x), data, 5)

fidSummary(fidsamples)
fidConfInt("map(exp, :x)", fidsamples)
fidConfInt(":x ./ :\"(Intercept)\"", fidsamples)
=#

#=
data = DataFrame(
  y = [0, 0, 1, 1, 1, 1],
  group = ["A", "A", "A", "B", "B", "B"]
)
fidsamples = fidSampleLR(@formula(y ~ 0 + group), data, 3000)
fidConfInt(":\"group: A\" - :\"group: B\"", fidsamples, 0.95)
=#

#=
    for i in 1:N
      @inbounds H = Polyhedra.hrep(CC[i], cc[i])
      plyhdrn = Polyhedra.polyhedron(H, CDDLib.Library(cdd))
      pts = collect(Polyhedra.points(plyhdrn))
      @inbounds if yK[t] == 0
        MIN = convert(
          Float64, minimum(qXtt * hcat(pts...))
        )
        atilde = rtlogis2(MIN)
        @inbounds weight[t, i] = 1 - expit(MIN)
        @inbounds CC[i] = vcat(CC[i], qXt_row)
        @inbounds cc[i] = vcat(cc[i], convert(T, atilde))
      else
        MAX = convert(
          Float64, maximum(qXtt * hcat(pts...))
        )
        atilde = rtlogis1(MAX)
        @inbounds weight[t, i] = expit(MAX)
        @inbounds CC[i] = vcat(CC[i], -qXt_row)
        @inbounds cc[i] = vcat(cc[i], convert(T, -atilde))
      end
      @inbounds At[p+t, i] = atilde
    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 ApplicationGatewayBackendAddressPoolPropertiesFormat2 <: SwaggerModel
    backendIPConfigurations::Any # spec type: Union{ Nothing, Vector{NetworkInterfaceIPConfiguration} } # spec name: backendIPConfigurations
    backendAddresses::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayBackendAddress2} } # spec name: backendAddresses
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function ApplicationGatewayBackendAddressPoolPropertiesFormat2(;backendIPConfigurations=nothing, backendAddresses=nothing, provisioningState=nothing)
        o = new()
        validate_property(ApplicationGatewayBackendAddressPoolPropertiesFormat2, Symbol("backendIPConfigurations"), backendIPConfigurations)
        setfield!(o, Symbol("backendIPConfigurations"), backendIPConfigurations)
        validate_property(ApplicationGatewayBackendAddressPoolPropertiesFormat2, Symbol("backendAddresses"), backendAddresses)
        setfield!(o, Symbol("backendAddresses"), backendAddresses)
        validate_property(ApplicationGatewayBackendAddressPoolPropertiesFormat2, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ApplicationGatewayBackendAddressPoolPropertiesFormat2

const _property_map_ApplicationGatewayBackendAddressPoolPropertiesFormat2 = Dict{Symbol,Symbol}(Symbol("backendIPConfigurations")=>Symbol("backendIPConfigurations"), Symbol("backendAddresses")=>Symbol("backendAddresses"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ApplicationGatewayBackendAddressPoolPropertiesFormat2 = Dict{Symbol,String}(Symbol("backendIPConfigurations")=>"Vector{NetworkInterfaceIPConfiguration}", Symbol("backendAddresses")=>"Vector{ApplicationGatewayBackendAddress2}", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ ApplicationGatewayBackendAddressPoolPropertiesFormat2 }) = collect(keys(_property_map_ApplicationGatewayBackendAddressPoolPropertiesFormat2))
Swagger.property_type(::Type{ ApplicationGatewayBackendAddressPoolPropertiesFormat2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayBackendAddressPoolPropertiesFormat2[name]))}
Swagger.field_name(::Type{ ApplicationGatewayBackendAddressPoolPropertiesFormat2 }, property_name::Symbol) =  _property_map_ApplicationGatewayBackendAddressPoolPropertiesFormat2[property_name]

function check_required(o::ApplicationGatewayBackendAddressPoolPropertiesFormat2)
    true
end

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