text:
"""
    certificationfinder(keyword::AbstractString,
                        userId::AbstractString = ENV["API_CAREERONESTOP_USER"],
                        token::AbstractString = ENV["API_CAREERONESTOP"],
                        directFlag::AbstractString = "0",
                        industry::AbstractString = "0",
                        certType::AbstractString = "0",
                        organization::AbstractString = "0",
                        occupation::AbstractString = "0",
                        agency::AbstractString = "0",
                        sortColumn::AbstractString = "0",
                        sortDirections::AbstractString = "0",
                        startRecord::Integer = 0,
                        limitRecord::Integer = 10,
                        )

Consumes the `List Certifications` endpoint of the CareerOneStop API.

# Example

```jldoctest
julia> nurse = certificationfinder("nurse")
JSON3.Object{Vector{UInt8}, Vector{UInt64}} with 13 entries:
  :CertList                  => JSON3.Object[{…
  :CertAccredAgencyLegend    => JSON3.Object[{…
  :RecordCount               => 561
  :DidYouMean                => nothing
  :AutoCorrection            => nothing
  :IsDirectMatch             => "K"
  :OccupationCountList       => {…
  :IndustryCountList         => {…
  :TypeCountList             => {…
  :OrgCountList              => {…
  :RelatedCountList          => {…
  :CertDetail                => nothing
  :CertAccredAgencyCountList => {…

julia> @views nurse.CertList[1:2]
2-element view(::JSON3.Array{JSON3.Object, Vector{UInt8}, SubArray{UInt64, 1, Vector{UInt64}, Tuple{UnitRange{Int64}}, true}}, 1:2) with eltype JSON3.Object:
 {
                     "Id": "6126-C",
                   "Name": "Gerontological Nursing",
                    "Url": "https://www.nursingworld.org/our-certifications/gerontological-nurse/",
           "Organization": "ANA Enterprise",
                   "Type": "Specialty",
              "TypeOrder": "C",
                "Related": "",
   "CertAccredAgencyList": [
                             {
                                 "Name": "In-Demand",
                                "Order": "1"
                             },
                             {
                                 "Name": "NCCA",
                                "Order": "2"
                             },
                             {
                                 "Name": "MILITARY",
                                "Order": "4"
                             },
                             {
                                 "Name": "ABNS",
                                "Order": "7"
                             }
                           ],
    "OrganizationAddress": "8515 Georgia Avenue, Suite 400, Silver Spring, MD  20910-3492",
        "OrganizationUrl": "https://www.nursingworld.org/",
                "Acronym": "RN-BC",
            "Description": "The ANCC Gerontological Nursing board certification examination is a competency based examination that provides a valid and reliable assessment of the entry-level clinical knowledge and skills of registered nurses in the gerontological specialty after initial RN licensure. Once you complete eligibility requirements to take the certification examination and successfully pass the exam, you are awarded the credential: Registered Nurse-Board Certified (RN-BC).",
         "CertDetailList": [
                             {
                                 "Name": "More than two years of education or training after high school required?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "More than two years of work experience required?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "Oral or Written Exam Required?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "Renewal Required?",
                                "Value": "Every 5  Year(s)"
                             },
                             {
                                 "Name": "Renew through Continuing Educational Units(CEU)?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "Renew through Re-Examination?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "Renew through Continuing Professional Development(CPD)?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "Does applicant have choice of at least two options from above for renewal (CEU, CPD, or exam)?",
                                "Value": "Yes"
                             }
                           ],
         "CertTestOption": [],
             "ExamDetail": nothing
}
 {
                     "Id": "6132-C",
                   "Name": "Pediatric Nursing",
                    "Url": "https://www.nursingworld.org/our-certifications/pediatric-nurse/",
           "Organization": "ANA Enterprise",
                   "Type": "Specialty",
              "TypeOrder": "C",
                "Related": "",
   "CertAccredAgencyList": [
                             {
                                 "Name": "In-Demand",
                                "Order": "1"
                             },
                             {
                                 "Name": "NCCA",
                                "Order": "2"
                             },
                             {
                                 "Name": "ANSI",
                                "Order": "3"
                             },
                             {
                                 "Name": "ABNS",
                                "Order": "7"
                             }
                           ],
    "OrganizationAddress": "8515 Georgia Avenue, Suite 400, Silver Spring, MD  20910-3492",
        "OrganizationUrl": "https://www.nursingworld.org/",
                "Acronym": "RN-BC",
            "Description": "The ANCC Pediatric Nursing board certification examination is a competency based examination that provides a valid and reliable assessment of the entry-level clinical knowledge and skills of registered nurses in the pediatric specialty after initial RN licensure. Once you complete eligibility requirements to take the certification examination and successfully pass the exam, you are awarded the credential: Registered Nurse-Board Certified (RN-BC).",
         "CertDetailList": [
                             {
                                 "Name": "More than two years of education or training after high school required?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "More than two years of work experience required?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "Oral or Written Exam Required?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "Renewal Required?",
                                "Value": "Every 5  Year(s)"
                             },
                             {
                                 "Name": "Renew through Continuing Educational Units(CEU)?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "Renew through Re-Examination?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "Renew through Continuing Professional Development(CPD)?",
                                "Value": "Yes"
                             },
                             {
                                 "Name": "Does applicant have choice of at least two options from above for renewal (CEU, CPD, or exam)?",
                                "Value": "Yes"
                             }
                           ],
         "CertTestOption": [],
             "ExamDetail": nothing
}

```
"""
function certificationfinder(keyword::AbstractString,
                             userID::AbstractString = ENV["API_CAREERONESTOP_USER"],
                             token::AbstractString = ENV["API_CAREERONESTOP"],
                             directFlag::AbstractString = "0",
                             industry::AbstractString = "0",
                             certType::AbstractString = "0",
                             organization::AbstractString = "0",
                             occupation::AbstractString = "0",
                             agency::AbstractString = "0",
                             sortColumn::AbstractString = "0",
                             sortDirections::AbstractString = "0",
                             startRecord::Integer = 0,
                             limitRecord = 10,
                             )
    response = request("GET",
                       URI(scheme = "https",
                           host = "api.careeronestop.org",
                           path = "/v1/certificationfinder/$userID/$keyword/$directFlag/$industry/$certType/$organization/$occupation/$agency/$sortColumn/$sortDirections/$startRecord/$limitRecord"),
                       ["Authorization" => "Bearer $token"],
                       )
    @assert response.status == 200
    json = JSON3.read(response.body)
end


text:
#!/usr/bin/env julia

using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--file", "-f"
        help = "col file to run"
        arg_type = String
        required = true
        "--time_limit", "-t"
        help = "Target commit id or branch"
        arg_type = Int
        default = 1800
    end

    return parse_args(s)
end

if isinteractive() == false
    args = parse_commandline()

    include("cs.jl")
    # for compiling
    main("/home/ole/Julia/ConstraintSolver/instances/eternity/pieces_05x05.txt")
    println("")
    println("======= ACTUAL RUN =======")
    println("")
    flush(stdout)
    # actual run
    main(args["file"]; time_limit = args["time_limit"])
end

ext: jl


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

mutable struct VirtualMachineSoftwarePatchProperties <: SwaggerModel
    patchId::Any # spec type: Union{ Nothing, String } # spec name: patchId
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    version::Any # spec type: Union{ Nothing, String } # spec name: version
    kbid::Any # spec type: Union{ Nothing, String } # spec name: kbid
    classifications::Any # spec type: Union{ Nothing, Vector{String} } # spec name: classifications
    rebootBehavior::Any # spec type: Union{ Nothing, String } # spec name: rebootBehavior
    activityId::Any # spec type: Union{ Nothing, String } # spec name: activityId
    publishedDate::Any # spec type: Union{ Nothing, DateTime } # spec name: publishedDate
    lastModifiedDateTime::Any # spec type: Union{ Nothing, DateTime } # spec name: lastModifiedDateTime
    assessmentState::Any # spec type: Union{ Nothing, String } # spec name: assessmentState

    function VirtualMachineSoftwarePatchProperties(;patchId=nothing, name=nothing, version=nothing, kbid=nothing, classifications=nothing, rebootBehavior=nothing, activityId=nothing, publishedDate=nothing, lastModifiedDateTime=nothing, assessmentState=nothing)
        o = new()
        validate_property(VirtualMachineSoftwarePatchProperties, Symbol("patchId"), patchId)
        setfield!(o, Symbol("patchId"), patchId)
        validate_property(VirtualMachineSoftwarePatchProperties, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(VirtualMachineSoftwarePatchProperties, Symbol("version"), version)
        setfield!(o, Symbol("version"), version)
        validate_property(VirtualMachineSoftwarePatchProperties, Symbol("kbid"), kbid)
        setfield!(o, Symbol("kbid"), kbid)
        validate_property(VirtualMachineSoftwarePatchProperties, Symbol("classifications"), classifications)
        setfield!(o, Symbol("classifications"), classifications)
        validate_property(VirtualMachineSoftwarePatchProperties, Symbol("rebootBehavior"), rebootBehavior)
        setfield!(o, Symbol("rebootBehavior"), rebootBehavior)
        validate_property(VirtualMachineSoftwarePatchProperties, Symbol("activityId"), activityId)
        setfield!(o, Symbol("activityId"), activityId)
        validate_property(VirtualMachineSoftwarePatchProperties, Symbol("publishedDate"), publishedDate)
        setfield!(o, Symbol("publishedDate"), publishedDate)
        validate_property(VirtualMachineSoftwarePatchProperties, Symbol("lastModifiedDateTime"), lastModifiedDateTime)
        setfield!(o, Symbol("lastModifiedDateTime"), lastModifiedDateTime)
        validate_property(VirtualMachineSoftwarePatchProperties, Symbol("assessmentState"), assessmentState)
        setfield!(o, Symbol("assessmentState"), assessmentState)
        o
    end
end # type VirtualMachineSoftwarePatchProperties

const _property_map_VirtualMachineSoftwarePatchProperties = Dict{Symbol,Symbol}(Symbol("patchId")=>Symbol("patchId"), Symbol("name")=>Symbol("name"), Symbol("version")=>Symbol("version"), Symbol("kbid")=>Symbol("kbid"), Symbol("classifications")=>Symbol("classifications"), Symbol("rebootBehavior")=>Symbol("rebootBehavior"), Symbol("activityId")=>Symbol("activityId"), Symbol("publishedDate")=>Symbol("publishedDate"), Symbol("lastModifiedDateTime")=>Symbol("lastModifiedDateTime"), Symbol("assessmentState")=>Symbol("assessmentState"))
const _property_types_VirtualMachineSoftwarePatchProperties = Dict{Symbol,String}(Symbol("patchId")=>"String", Symbol("name")=>"String", Symbol("version")=>"String", Symbol("kbid")=>"String", Symbol("classifications")=>"Vector{String}", Symbol("rebootBehavior")=>"String", Symbol("activityId")=>"String", Symbol("publishedDate")=>"DateTime", Symbol("lastModifiedDateTime")=>"DateTime", Symbol("assessmentState")=>"String")
Base.propertynames(::Type{ VirtualMachineSoftwarePatchProperties }) = collect(keys(_property_map_VirtualMachineSoftwarePatchProperties))
Swagger.property_type(::Type{ VirtualMachineSoftwarePatchProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineSoftwarePatchProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineSoftwarePatchProperties }, property_name::Symbol) =  _property_map_VirtualMachineSoftwarePatchProperties[property_name]

const _allowed_VirtualMachineSoftwarePatchProperties_rebootBehavior = ["NeverReboots", "AlwaysRequiresReboot", "CanRequestReboot"]

const _allowed_VirtualMachineSoftwarePatchProperties_assessmentState = ["Installed", "Failed", "Excluded", "NotSelected", "Pending", "Available"]

function check_required(o::VirtualMachineSoftwarePatchProperties)
    true
end

function validate_property(::Type{ VirtualMachineSoftwarePatchProperties }, name::Symbol, val)
    if name === Symbol("rebootBehavior")
        Swagger.validate_param(name, "VirtualMachineSoftwarePatchProperties", :enum, val, _allowed_VirtualMachineSoftwarePatchProperties_rebootBehavior)
    end
    if name === Symbol("assessmentState")
        Swagger.validate_param(name, "VirtualMachineSoftwarePatchProperties", :enum, val, _allowed_VirtualMachineSoftwarePatchProperties_assessmentState)
    end
end


text:
using POMDPs
using POMCPOWRAVE
using ProfileView
using POMDPModels

#=
using Gallium
breakpoint(Pkg.dir("POMCPOWRAVE", "src", "solver2.jl"), 100)
=#

solver = POMCPOWRAVESolver(tree_queries=250_000,
                     eps=0.01,
                     criterion=MaxUCB(10.0),
                     enable_action_pw=true,
                     check_repeat_obs=false,
                     alpha_observation=1/8,
                     rng=MersenneTwister(2))

problem = LightDark1D()
policy = POMCPOWRAVEPlanner(solver, problem)
ib = initial_state_distribution(problem)

a = action(policy, ib)

@time a = action(policy, ib)

Profile.clear()
@profile a = action(policy, ib)
ProfileView.view()

#=
P = typeof(problem)
S = state_type(P)
A = action_type(P)
O = obs_type(P)
tree = POMCPOWRAVETree{POWRAVENodeBelief{S,A,O,P},A,O,typeof(ib)}(ib, 500_000)

@code_warntype POMCPOWRAVE.simulate(policy, POMCPOWRAVE.POWRAVETreeObsNode(tree, 1), rand(Base.GLOBAL_RNG, ib), 10)
=#

#=
solver = POMCPOWRAVESolver(tree_queries=30,
                     eps=0.01,
                     criterion=MaxUCB(10.0),
                     enable_action_pw=true,
                     check_repeat_obs=false,
                     alpha_observation=1/8,
                     rng=MersenneTwister(2))

policy = POMCPOWRAVEPlanner(solver, problem)
a = action(policy, ib)
blink(policy)
=#


text:
#=##############################################################################
# DESCRIPTION
    Simulation postprocessing

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jan 2020
  * License   : MIT
=###############################################################################

"""
    Given the path of a simulation `read_path`, it will read the particle field
in time steps `nums` and it will calculate and output the velocity and vorticity
field in a prescribed fluid domain.
"""
function postprocessing_fluiddomain(read_path::String, nums, save_path::String;
                                    b=1.693, R=0.127, Vvehicle=23.384*[-1, 0, 0],
                                    dt=0.195516/2160, ndx=10,
                                    run_name="blownwing", Vinf=(x,t)->zeros(3),
                                    verbose=true, v_lvl=0, prompt=true)

    gt.create_path(save_path, prompt)

    # ------------------- Generate fluid domain grid ---------------------------
    if verbose; println("\t"^(v_lvl)*"Creating grid..."); end;
    lx = 1.00*b                         # Length in x-direction
    ly = 1.10*b                         # Length in y-direction
    lz = 2.35*R                         # Length in z-direction

    dx = R/ndx                          # Element size in x-direction
    dy = dx                             # Element size in y-direction
    dz = dx                             # Element size in z-direction

    P_min = [0.0, -ly/2, -lz/2*1.1]         # Minimum bounding point of grid
    P_max = P_min + [lx, ly, lz]        # Maximum bounding point of grid

    NDIVS = ceil.(Int, [lx/dx, ly/dy, lz/dz])   # Number of cells in each direction

    grid = gt.Grid(P_min, P_max, NDIVS)

    # ------------------- Evaluate fuild domain --------------------------------
    if verbose; println("\t"^(v_lvl)*"Evaluating fluid domain...\t$(grid.nnodes) nodes"); end;

    prevX = zeros(3)

    for (i, num) in enumerate(nums)
        if verbose; print("\t"^(v_lvl+1)*"Processing step $i of $(length(nums)):"); end;

        curX = (dt*num)*Vvehicle  # Current position of grid

        # Translate grid
        gt.lintransform!(grid, eye(3), curX-prevX)

        # Read particle field
        pfield = uns.vpm.read_vtk(run_name*"_pfield", Vinf, nothing; path=read_path,
                                                num=num, solver_method="ExaFMM")

        # Eliminates particles outside the domain to save computation
        xmax = gt.get_node(grid, NDIVS+1)[1]
        for p in uns.vpm.get_np(pfield):-1:1
            if uns.vpm.get_x(pfield, p)[1]>xmax
                uns.vpm.delparticle(pfield, p)
            end
        end

        if verbose; println("\t$(uns.vpm.get_np(pfield)) particles"); end;

        # Evaluate particle field velocity
        nodes = [gt.get_node(grid, n) for n in 1:grid.nnodes]
        @time _, U, W, dUdX = uns.vpm.conv(pfield, "ExaFMM"; Uprobes=nodes, vorticity=true)

        gt.add_field(grid, "U", "vector", U, "node"; raise_warn=false)
        gt.add_field(grid, "W", "vector", W, "node"; raise_warn=false)
        gt.save(grid, run_name*"_Grid"; path=save_path, num=num)

        prevX .= curX
    end

    return nothing
end


text:
"Holds the tableau of an Specialised Partitioned Additive Runge-Kutta method for Variational systems."
struct TableauVSPARKsecondary{DT <: Number} <: AbstractTableau{DT}
    name::Symbol
    o::Int
    s::Int
    σ::Int

    q::CoefficientsSPARK{DT}
    p::CoefficientsSPARK{DT}

    q̃::CoefficientsSPARK{DT}
    p̃::CoefficientsSPARK{DT}

    ω::Matrix{DT}
    d::Vector{DT}

    function TableauVSPARKsecondary(name::Symbol, o::Int, s::Int, σ::Int,
                        q::CoefficientsSPARK{DT}, p::CoefficientsSPARK{DT},
                        q̃::CoefficientsSPARK{DT}, p̃::CoefficientsSPARK{DT},
                        ω::Matrix{DT}, d=DT[]) where {DT}

        @assert s > 0 "Number of stages s must be > 0"
        @assert σ > 0 "Number of stages σ must be > 0"

        @assert s==q.s==p.s==q̃.σ==p̃.σ
        @assert σ==q.σ==p.σ==q̃.s==p̃.s
        @assert size(ω,1)==σ
        @assert size(ω,2)==σ+1

        @assert length(d)==0 || length(d)==σ

        new{DT}(name, o, s, σ, q, p, q̃, p̃, ω, d)
    end
end

"Parameters for right-hand side function of Specialised Partitioned Additive Runge-Kutta methods for Variational systems."
mutable struct ParametersVSPARKsecondary{DT,TT,D,S,Σ,θT,FT,GT,G̅T,ϕT,ψT,tabType} <: Parameters{DT,TT}
    f_ϑ::θT
    f_f::FT
    f_g::GT
    f_g̅::G̅T
    f_ϕ::ϕT
    f_ψ::ψT

    Δt::TT

    tab::tabType

    t::TT
    q::Vector{DT}
    p::Vector{DT}
    λ::Vector{DT}

    function ParametersVSPARKsecondary{DT,D,S,Σ}(f_ϑ::θT, f_f::FT, f_g::GT, f_g̅::G̅T, f_ϕ::ϕT, f_ψ::ψT, Δt::TT, tab::tabType) where {DT,TT,D,S,Σ,θT,FT,GT,G̅T,ϕT,ψT,tabType}
        # create solution vectors
        q = zeros(DT,D)
        p = zeros(DT,D)
        λ = zeros(DT,D)

        new{DT,TT,D,S,Σ,θT,FT,GT,G̅T,ϕT,ψT,tabType}(f_ϑ, f_f, f_g, f_g̅, f_ϕ, f_ψ, Δt, tab, zero(TT), q, p, λ)
    end
end

function update_params!(params::ParametersVSPARKsecondary, sol::AtomicSolutionPDAE)
    # set time for nonlinear solver and copy previous solution
    params.t  = sol.t
    params.q .= sol.q
    params.p .= sol.p
    params.λ .= sol.λ
end

@doc raw"""
Specialised Partitioned Additive Runge-Kutta integrator for degenerate
variational systems with projection on secondary constraint.

This integrator solves the following system of equations for the internal stages,
```math
\begin{align}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a^1_{ij} V_{n,j} + h \sum \limits_{j=1}^{\sigma} a^2_{ij} \tilde{\Lambda}_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \bar{a}^1_{ij} F^1_{n,j} + h \sum \limits_{j=1}^{s} \bar{a}^2_{ij} F^2_{n,j} + h \sum \limits_{j=1}^{\sigma} \bar{a}^3_{ij} \tilde{F}^3_{n,j} , & i &= 1, ..., s , \\
0 &= \Phi_{n,i} , & i &= 1, ..., s ,
\end{align}
```
the projective stages
```math
\begin{align}
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \alpha^1_{ij} \tilde{V}_{n,j} + h \sum \limits_{j=1}^{\sigma} \alpha^2_{ij} \tilde{\Lambda}_{n,j} , & i &= 1, ..., \sigma , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \bar{\alpha}^1_{ij} F^1_{n,j} + h \sum \limits_{j=1}^{s} \bar{\alpha}^2_{ij} \tilde{F}^2_{n,j} + h \sum \limits_{j=1}^{\sigma} \bar{\alpha}^3_{ij} \tilde{F}^3_{n,j} , & i &= 1, ..., \sigma , \\
0 &= \tilde{\Phi}_{n,i} , & i &= 1, ..., \sigma , \\
0 &= \sum \limits_{j=1}^{\sigma} \omega_{ij} \tilde{\Psi}_{n,i} , & i &= 1, ..., \sigma-1 ,
\end{align}
```
and update rule
```math
\begin{align}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b^1_{i} V_{n,i} + h \sum \limits_{i=1}^{\sigma} b^2_{i} \tilde{\Lambda}_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b^1_{i} F^1_{n,i} + h \sum \limits_{i=1}^{s} b^2_{i} F^2_{n,i} + h \sum \limits_{i=1}^{\sigma} b^3_{i} \tilde{F}^3_{n,i} , \\
0 &= \phi (q_{n+1}, p_{n+1}) ,
\end{align}
```
with definitions
```math
\begin{align}
F^1_{n,i} &= \nabla H (Q_{n,i}) , & i &= 1, ..., s , \\
F^2_{n,i} &= \nabla \vartheta (Q_{n,i}) \cdot V_{n,i} , & i &= 1, ..., s , \\
\tilde{F}^3_{n,i} &= \nabla \phi (\tilde{Q}_{n,i}, \tilde{P}_{n,i}) \cdot \tilde{\Lambda}_{n,i} , & i &= 1, ..., \sigma , \\
\Phi_{n,i} &= \phi (Q_{n,i}, P_{n,i}) = P_{n,i} - \vartheta (Q_{n,i}) , & i &= 1, ..., s , \\
\tilde{\Phi}_{n,i} &= \phi (\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = \tilde{P}_{n,i} - \vartheta (\tilde{Q}_{n,i}) , & i &= 1, ..., \sigma , \\
\tilde{\Psi}_{n,i} &= \psi (\tilde{Q}_{n,i}, \tilde{V}_{n,i}, \tilde{P}_{n,i}, \tilde{F}_{n,i}) = \tilde{F}_{n,i} - \tilde{V}_{n,i} \cdot \nabla \vartheta (\tilde{Q}_{n,i}) , & i &= 1, ..., \sigma ,
\end{align}
```
so that
```math
\begin{align}
F^1_{n,i} + F^2_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
\phi (Q_{n,i}, P_{n,i}) &= P_{n,i} - \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
\psi (\tilde{Q}_{n,i}, \tilde{V}_{n,i}, \tilde{P}_{n,i}, \tilde{F}_{n,i}) &= \tilde{F}_{n,i} - \tilde{V}_{n,i} \cdot \nabla \frac{\partial L}{\partial v} (\tilde{Q}_{n,i}, \tilde{V}_{n,i}) , & i &= 1, ..., \sigma .
\end{align}
```
"""
struct IntegratorVSPARKsecondary{DT, TT, ET <: VDAE{DT,TT},
                                         PT <: ParametersVSPARKsecondary{DT,TT},
                                         ST <: NonlinearSolver{DT},
                                         IT <: InitialGuessPODE{DT,TT}, D, S, Σ} <: AbstractIntegratorVSPARK{DT, TT}
    equation::ET
    tableau::TableauVSPARKsecondary{TT}

    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCacheSPARK{DT,TT,D,S,Σ}
end

function IntegratorVSPARKsecondary(equation::VDAE{DT,TT},
                                   tableau::TableauVSPARKsecondary{TT}, Δt::TT) where {DT,TT}
    D = equation.d
    S = tableau.s
    Σ = tableau.σ

    # println(tableau)

    N = 4*D*Σ

    if isdefined(tableau, :d) && length(tableau.d) > 0
        N += D
    end

    # create params
    params = ParametersVSPARKsecondary{DT,D,S,Σ}(equation.ϑ, equation.f, equation.g, equation.g̅, equation.ϕ, equation.ψ, Δt, tableau)

    # create solver
    solver = create_nonlinear_solver(DT, N, params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create cache
    cache = IntegratorCacheSPARK{DT,TT,D,S,Σ}()

    # create integrator
    IntegratorVSPARKsecondary{DT, TT, typeof(equation), typeof(params), typeof(solver), typeof(iguess), D, S, Σ}(
                              equation, tableau, params, solver, iguess, cache)
end

pstages(int::IntegratorVSPARKsecondary) = int.tableau.σ

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,TT,D,S,Σ},
                                        params::ParametersVSPARKsecondary{DT,TT,D,S,Σ}) where {ST,DT,TT,D,S,Σ}
    local t::TT

    for i in 1:Σ
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[4*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[4*(D*(i-1)+k-1)+2]
            cache.Vp[i][k] = x[4*(D*(i-1)+k-1)+3]
            cache.Λp[i][k] = x[4*(D*(i-1)+k-1)+4]

            # compute Q and P
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        t = params.t + params.Δt * params.tab.p̃.c[i]
        params.f_f(t, cache.Qp[i], cache.Pp[i], cache.Fp[i])
        params.f_g(t, cache.Qp[i], cache.Pp[i], cache.Vp[i], cache.Gp[i])
        params.f_g(t, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.G̅p[i])

        cache.Hp[i] .= cache.Fp[i] .+ cache.Gp[i]

        params.f_ϕ(t, cache.Qp[i], cache.Pp[i], cache.Φp[i])
        params.f_ψ(t, cache.Qp[i], cache.Pp[i], cache.Vp[i], cache.Hp[i], cache.Ψp[i])
    end

    if length(params.tab.d) > 0
        for k in 1:D
            cache.μ[k] = x[4*D*Σ+k]
        end
    end

    for i in 1:S
        # compute Q
        for k in 1:D
            cache.Yi[i][k] = 0
            for j in 1:Σ
                cache.Yi[i][k] += params.tab.q.a[1][i,j] * cache.Vp[j][k]
                cache.Yi[i][k] += params.tab.q.a[2][i,j] * cache.Λp[j][k]
            end
            cache.Qi[i][k] = params.q[k] + params.Δt * cache.Yi[i][k]
        end

        # compute f(X)
        t = params.t + params.Δt * params.tab.p.c[i]
        params.f_f(t, cache.Qi[i], cache.Fi[i])
    end

    # compute q and p
    cache.q̃ .= params.q
    cache.p̃ .= params.p
    for i in 1:S
        cache.p̃ .+= params.Δt .* params.tab.p.b[1][i] .* cache.Fi[i]
    end
    for i in 1:Σ
        cache.q̃ .+= params.Δt .* params.tab.q.b[1][i] .* cache.Vp[i]
        cache.q̃ .+= params.Δt .* params.tab.q.b[2][i] .* cache.Λp[i]
        cache.p̃ .+= params.Δt .* params.tab.p.b[2][i] .* cache.Gp[i]
        cache.p̃ .+= params.Δt .* params.tab.p.b[3][i] .* cache.G̅p[i]
    end

    # compute ϕ(q,p)
    t = params.t + params.Δt
    params.f_ϕ(t, cache.q̃, cache.p̃, cache.ϕ̃)
end

"Compute stages of specialised partitioned additive Runge-Kutta methods for variational systems."
@generated function function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersVSPARKsecondary{DT,TT,D,S,Σ}) where {ST,DT,TT,D,S,Σ}
    cache = IntegratorCacheSPARK{ST,TT,D,S,Σ}()

    quote
        compute_stages!(y, $cache, params)

        # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ, ωΨ]
        for i in 1:Σ
            for k in 1:D
                b[4*(D*(i-1)+k-1)+1] = - $cache.Yp[i][k]
                b[4*(D*(i-1)+k-1)+2] = - $cache.Zp[i][k]
                b[4*(D*(i-1)+k-1)+3] = - $cache.Φp[i][k]
                b[4*(D*(i-1)+k-1)+4] = 0
                for j in 1:S
                    b[4*(D*(i-1)+k-1)+2] += params.tab.p̃.a[1][i,j] * $cache.Fi[j][k]
                end
                for j in 1:Σ
                    b[4*(D*(i-1)+k-1)+1] += params.tab.q̃.a[1][i,j] * $cache.Vp[j][k]
                    b[4*(D*(i-1)+k-1)+1] += params.tab.q̃.a[2][i,j] * $cache.Λp[j][k]
                    b[4*(D*(i-1)+k-1)+2] += params.tab.p̃.a[2][i,j] * $cache.Gp[j][k]
                    b[4*(D*(i-1)+k-1)+2] += params.tab.p̃.a[3][i,j] * $cache.G̅p[j][k]
                end
                for j in 1:Σ
                    b[4*(D*(i-1)+k-1)+4] -= params.tab.ω[i,j] * $cache.Ψp[j][k]
                end
                b[4*(D*(i-1)+k-1)+4] -= params.tab.ω[i,Σ+1] * $cache.ϕ̃[k]
            end
        end

        if length(params.tab.d) > 0
            for i in 1:Σ
                for k in 1:D
                    b[4*(D*(i-1)+k-1)+2] -= $cache.μ[k] * params.tab.d[i] / params.tab.p.b[2][i]
                end
            end

            for k in 1:D
                b[4*D*Σ+k] = 0
                for i in 1:Σ
                    b[4*D*Σ+k] -= $cache.Vp[i][k] * params.tab.d[i]
                end
            end
        end
    end
end

function initial_guess!(int::IntegratorVSPARKsecondary, sol::AtomicSolutionPDAE)
    for i in 1:pstages(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.p̃, int.cache.ṽ, int.cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in 1:ndims(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+1] = (int.cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+2] = (int.cache.p̃[k] - sol.p[k])/timestep(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+3] = int.cache.ṽ[k]
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+4] = 0
        end
    end

    if isdefined(tableau(int), :d) && length(tableau(int).d) > 0
        for k in 1:ndims(int)
            int.solver.x[4*ndims(int)*pstages(int)+k] = 0
        end
    end
end

function update_solution!(int::IntegratorVSPARKsecondary{DT,TT}, sol::AtomicSolutionPDAE{DT,TT}) where {DT,TT}
    # compute final update
    update_solution!(sol.p, sol.p̃, int.cache.Fi, int.params.tab.p.b[1], timestep(int))

    # compute projection
    update_solution!(sol.q, sol.q̃, int.cache.Vp, int.params.tab.q.b[1], timestep(int))
    update_solution!(sol.q, sol.q̃, int.cache.Λp, int.params.tab.q.b[2], timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.Gp, int.params.tab.p.b[2], timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.G̅p, int.params.tab.p.b[3], timestep(int))
end


text:
import Flux.testmode!
import Base.deepcopy

include("resnet.jl")

mutable struct NeuralNet
  base_net::Chain
  value::Chain
  policy::Chain
  opt

  function NeuralNet(env::T; base_net = nothing, value = nothing, policy = nothing,
                          tower_height::Int = 19) where T <: GameEnv
    if base_net == nothing
      res_block() = ResidualBlock([256,256,256], [3,3], [1,1], [1,1])
      # 19 residual blocks
      tower = [res_block() for i = 1:tower_height]
      base_net = Chain(Conv((3,3), env.planes=>256, pad=(1,1)), BatchNorm(256, relu),
                        tower...) |> gpu
    end
    if value == nothing
      value = Chain(Conv((1,1), 256=>1), BatchNorm(1, relu), x->reshape(x, :, size(x, 4)),
                    Dense(env.N*env.N, 256, relu), Dense(256, 1, tanh)) |> gpu
    end
    if policy == nothing
      policy = Chain(Conv((1,1), 256=>2), BatchNorm(2, relu), x->reshape(x, :, size(x, 4)),
                      Dense(2env.N*env.N, env.action_space), x -> softmax(x)) |> gpu
    end

    all_params = vcat(params(base_net), params(value), params(policy))
    opt = Momentum(all_params, 0.02f0)
    new(base_net, value, policy, opt)
  end
end

#=
function deepcopy(nn::NeuralNet)
  base_net = deepcopy(nn.base_net)
  value = deepcopy(nn.value)
  policy = deepcopy(nn.policy)

  return NeuralNet(;base_net = base_net, value = value, policy = policy)
end
=#

function testmode!(nn::NeuralNet, val::Bool=true)
  testmode!(nn.base_net, val)
  testmode!(nn.policy, val)
  testmode!(nn.value, val)
end

function (nn::NeuralNet)(input::Vector{T}, train = false) where T <: Position
  nn_in = cat(4, get_feats.(input)...) |> gpu
  if !train testmode!(nn) end
  common_out = nn.base_net(nn_in |> gpu)
  π, val = nn.policy(common_out), nn.value(common_out)
  if !train testmode!(nn, false) end
  return π, val
end

function (nn::NeuralNet)(input::Position)
  p, v = nn([input])
  return p[:, 1], v[1]
end

loss_π(π, p) = crossentropy(p, π; weight = 0.01f0)

loss_value(z, v) = 0.01f0 * mse(z, v)

function loss_reg(nn::NeuralNet)
  sum_sqr(x) = sum([sum(i.^2) for i in x])
   0.0001f0 * (sum_sqr(params(nn.base_net)) + sum_sqr(params(nn.value)) + sum_sqr(params(nn.policy)))
end

function train!(nn::NeuralNet, input_data::Tuple{Vector{Position}, Matrix{Float32}, Vector{Int}}; epochs = 1)
  positions = input_data[1]
  π, z = input_data[2:3]
  loss_avg = 0.0f0
  data_size = length(z)
  for i = 1:epochs
    for j = 1:32:data_size
      p, v = nn(positions[j:j+31], true)
      loss = loss_π(cu(π[:, j:j+31]), p) + loss_value(cu(z[j:j+31]),v) + loss_reg(nn)
      back!(loss)
      loss_avg += loss.tracker.data
      nn.opt()
    end
  end
  return loss_avg / epochs
end

function evaluate(env::GameEnv, black_net::NeuralNet, white_net::NeuralNet; num_games = 400, ro = 800)
  games_won = 0

  testmode!(black_net)
  testmode!(white_net)

  black = MCTSPlayer(env, black_net, num_readouts = ro, two_player_mode = true)
  white = MCTSPlayer(env, white_net, num_readouts = ro, two_player_mode = true)

  for i = 1:num_games
    num_move = 0  # The move number of the current game

    initialize_game!(black)
    initialize_game!(white)

    while true
      active = num_move % 2 == true ? white : black
      inactive = num_move % 2 == true? black : white

      current_readouts = N(active.root)
      readouts = active.num_readouts

      while N(active.root) < current_readouts + readouts
        tree_search!(active)
      end

      # First, check the roots for hopeless games.
      if should_resign(active)  # Force resign
        set_result!(active, -active.root.position.to_play, true)
        set_result!(inactive, -active.root.position.to_play, true)
	break
      end

      move = pick_move(active)
      play_move!(active, move)
      play_move!(inactive, move)
      num_move += 1

      if is_done(active)
        winner = result(active.root.position)
        set_result!(active,winner,false)
	set_result!(inactive, winner, false)
	break
      end
    end
    games_won += result(black.root.position) == BLACK
  end

  testmode!(black_net, false)
  testmode!(white_net, false)
  print("Won $games_won / $num_games. Win rate $(games_won/num_games). ")
  return games_won / num_games ≥ 0.55
end


text:
module Empirikos

using Reexport

import Base: broadcast, broadcast!, broadcasted, eltype, zero, <=
using DataStructures
@reexport using Distributions
import Distributions:
    ntrials, pdf, support, location, cf, cdf, ccdf, logpdf, logdiffcdf, logccdf, components

import Intervals: Interval, Closed, Open, Unbounded, Bounded, AbstractInterval, isbounded,
    RightEndpoint
export Interval, Closed, Open, Unbounded # instead of @reexport

import JuMP
import JuMP: @constraint, @variable, set_lower_bound, @expression,
    Model, @objective, optimize!, objective_value, set_objective,
    set_normalized_rhs, RotatedSecondOrderCone, SecondOrderCone, set_value

using KernelDensity
using LinearAlgebra
using LinearFractional
using MathOptInterface
using Optim
using ParameterJuMP
using QuadGK
using Random
using RecipesBase
using Setfield
using Statistics
import Statistics: std, var

using StatsBase
import StatsBase: loglikelihood, response, fit, nobs, weights, confint

using UnPack

include("set_defaults.jl")
include("ebayes_samples.jl")
include("compound.jl")
include("interval_discretizer.jl")
include("dict_function.jl")
include("ebayes_methods.jl")
include("ebayes_targets.jl")
include("mixtures.jl")
include("convex_priors.jl")
include("flocalizations.jl")
include("NPMLE.jl")
include("samples/binomial.jl")
include("samples/normal.jl")
include("samples/poisson.jl")
include("samples/truncatedpoisson.jl")
include("example_priors.jl")
include("confidence_interval_tools.jl")
include("flocalization_intervals.jl")
include("flocalization_kde.jl")
include("amari.jl")

include("datasets/LordCressie/LordCressie.jl")
include("datasets/Prostate/Prostate.jl")
include("datasets/Neighborhoods/neighborhoods.jl")
include("datasets/Butterfly/Butterfly.jl")
include("datasets/Surgery/Surgery.jl")

export EBayesSample,
    NormalSample,
    StandardNormalSample,
    BinomialSample,
    PoissonSample,
    TruncatedPoissonSample,
    marginalize,
    compound,
    discretize,
    summarize,
    likelihood,
    likelihood_distribution,
    EBInterval,
    nuisance_parameter,
    skedasticity,
    MethodOfMoments,
    ParametricMLE,
    PosteriorMean,
    PosteriorVariance,
    MarginalDensity,
    PriorDensity,
    DiscretePriorClass,
    MixturePriorClass,
    GaussianScaleMixtureClass,
    NPMLE,
    nominal_alpha,
    integer_discretizer,
    interval_discretizer

export loglikelihood,
    response,
    fit,
    nobs,
    weights,
    confint

# F-Localizations
export DvoretzkyKieferWolfowitz,
    ChiSquaredFLocalization,
    InfinityNormDensityBand

# utilities
export DictFunction

# default

export DataBasedDefault

export FLocalizationInterval,
       AMARI

end


text:
using HDF5

struct Instance
    clients
    stations
    Γ
    ϕ
    Σ
    Π
    k
    c
    T
end

function load_istance(file_name, instance)
    h5open(file_name) do file
        inst = file[instance]
        clients = read(inst["clients"])
        stations = read(inst["stations"])
        Γ = read(inst["gamma"])
        ϕ = read(inst["phi"])
        Σ = read(inst["sigma"])
        Π = read(inst["pi"])
        k = read(inst["k"])
        c = read(inst["c"])
        T = read(inst["T"])
        Instance(clients, stations, Γ, ϕ, Σ, Π, k, c, T)
    end
end


text:
function toh06a()
#
#  generalization of K.C. Toh and L. N. Trefethen
#
    p = ones(11)
    p = conv(p,p)
    z = [0.84125353283118 + 0.54064081745560im
         0.84125353283118 - 0.54064081745560im
         0.41541501300189 + 0.90963199535452im
         0.41541501300189 - 0.90963199535452im
        -0.95949297361449 + 0.28173255684143im
        -0.95949297361449 - 0.28173255684143im
        -0.65486073394529 + 0.75574957435426im
        -0.65486073394529 - 0.75574957435426im
        -0.14231483827329 + 0.98982144188093im
        -0.14231483827329 - 0.98982144188093im];
    z = [z 2*ones(10)]

    p, PolyZeros(z)
end


text:
using Plots, DataFrames, Roots

import RecipesBase.plot

function plot(ht::pwr.htest; backend::Function = gr)

    methodlist = (
        "One-sample t-test power calculation",
        "Two-sample t-test power calculation",
        "Paired t-test power calculation",
        "T-test power calculation",
        "Difference of proportion power calculation for binomial distribution (same sample)",
        "Difference of proportion power calculation for binomial distribution (different sample)",
        "Balanced one-way analysis of variance power calculation",
        "Chi-square test power calculation",
        "Mean power calculation for normal distribution with known variance",
        "Proportion power calculation for binomial distribution (arcsine transformation)",
        "Approximate correlation power calculation (arctangh transformation)"
    )

    if backend in (gr, pyplot, plotly, plotlyjs) == false
        error(backend," is not supported")
    end

    if backend in (gr, pyplot)
        sepstr = "\n"
    else
        sepstr = "<br>"
    end

    ytitle = "Test Power = 1 - β"
    if backend == gr
        ytitle = "Test Power = 1 - \\beta"
    end

    d = ht.d

    if !(ht.title in methodlist)
        warn("The method `",ht.title,"` is not supported. Supported methods include: \n\t",join(methodlist,"\n\t"))
        exit(0)
    end

    breaks = 20

    # case: One-sample, Two-sample or Paired t test
    if ht.title in ("One-sample t-test power calculation",
        "Two-sample t-test power calculation",
        "Paired t-test power calculation")

        n = d["n"]
        n_upper = ceil(Int64,max(n*1.5, n+30))
        n_increment = ceil(Int64,(n_upper - 10)/breaks)
        sample_sizes = collect(10:n_increment:n_upper)

        # obtain power at each sample size
        power = [ss->powerTTest(n=ss,d=d["d"],alpha = d["alpha"],sampletype=d["sampletype"],alternative=d["alternative"]) for ss in sample_sizes]

        # create labels
        title = ht.title
        legend_string = string(
                    "tails = ", d["alternative"],sepstr,
                    "effect size d = ", d["d"],sepstr,
                    "alpha = ", d["alpha"])
        optimal_string = string(
                    "optimal sample size",sepstr,
                    "n = ", ceil(Int64,n),sepstr,d["note"])

    # case: Two-sample t test with n1 and n2
    elseif ht.title == "T-test power calculation"
        n = d["n1"] + d["n2"]
        n_upper = ceil(Int64,max(n*1.5, n+30)) # upper at least 30 above n
        n_rel = d["n1"] / n # relative sample size; will be kept constant in claculations

        # generate data
        n_increment = ceil(Int64,(n_upper - 10)/breaks)
        n_lower = ceil(Int64,fzero(x->(x*n_rel) - 2.1,2,d["n2"]))
        sample_sizes = collect(n_lower:n_increment:n_upper)

        function pwrt2test(ss)
            na = ceil(Int64,ss*n_rel)
            nb = ss - na
            if (na < 2 || nb < 2)
                return 0.0
            end
            return powerT2nTest(n1=na, n2=nb, d=d["d"], alpha = d["alpha"], alternative = d["alternative"])
        end

        power = [pwrt2test(ss) for ss in sample_sizes]

        # create labels
        title = ht.title
        legend_string = string(
                    "tails = ", d["alternative"], sepstr,
                    "effect size d = ", d["d"], sepstr,
                    "alpha = ", d["alpha"], sepstr,
                    "n1/n2 = ",round(n_rel,2))
        optimal_string = string("optimal sample size",sepstr,
                    "n = ",  d["n1"], " + ", d["n2"], " = ", n)

    # case: Difference of proportion (same sample size)
    elseif ht.title == "Difference of proportion power calculation for binomial distribution (same sample)"
        n = d["n"]
        n_upper = ceil(Int64, max(n*1.5, n+30)) # upper at least 30 above n
        n_increment = ceil(Int64,(n_upper - 10)/breaks)

        # generate data
        sample_sizes = collect(10:n_increment:n_upper)
        power = [ ss-> powerTwopTest(n=ss, h=d["h"], alpha = d["alpha"], alternative = d["alternative"]) for ss in sample_sizes]

        # create labels
        title = string("Difference of proportion power calculation",sepstr,"for binomial distribution (same sample)")
        legend_string = string(
                    "tails = ", d["alternative"], sepstr,
                    "effect size h = ", d["h"], sepstr,
                    "alpha = ", d["alpha"])
        optimal_string = string("optimal sample size", sepstr,
                    "n = ", ceil(Int64,n), sepstr,d["note"])

    # case: difference of proportion (different sample size)
    elseif ht.title == "Difference of proportion power calculation for binomial distribution (different sample)"
        n = d["n1"] + d["n2"]
        n_upper = ceil(Int64, max(n*1.5, n+30)) # upper at least 30 above n
        n_rel = d["n1"] / n # relative sample size; will be kept constant in claculations
        n_increment = ceil(Int64,(n_upper - 10)/breaks)
        n_lower = ceil(Int64,fzero(x->(x*n_rel) - 2.1,2,d["n2"]))

        # generate data
        sample_sizes = collect(n_lower:n_increment:n_upper)

        function pwr2p2ntest(ss)
            na = ceil(Int64,ss*n_rel)
            nb = ss - na
            if (na < 2 || nb < 2)
                return (0.0)
            end
            return powerTwop2nTest(n1=na, n2=nb, h=d["h"], alpha = d["alpha"], alternative = d["alternative"])
        end

        power = [ss->pwr2p2ntest(ss) for ss in sample_sizes]

        # create labels
        title = string("Difference of proportion power calculation",sepstr,"for binomial distribution (different sample)")
        legend_string = string(
                    "tails = ", d["alternative"],sepstr,
                    "effect size h = ", d["h"], sepstr,
                    "alpha = ", d["alpha"], sepstr,
                    "n1/n2 = ",round(n_rel,2))
        optimal_string = string("optimal sample size", sepstr,
                    "n = ",  d["n1"], " + ", d["n2"], " = ", n)

    # case: ANOVA
    elseif ht.title== "Balanced one-way analysis of variance power calculation"
        n = d["n"]
        n_upper = ceil(Int64, max(n*1.5, n+30)) # upper at least 30 above n
        n_increment = ceil(Int64,(n_upper - 10)/breaks)

        # generate data
        sample_sizes = collect(10:n_increment:n_upper)
        power = [ss->powerAnovaTest(n=ss, k=d["k"], f=d["f"], alpha = d["alpha"]) for ss in sample_sizes]

        # create labels
        title = ht.title
        legend_string = string(
                    "groups k = ", d["k"], sepstr,
                    "effect size f = ", d["f"], sepstr,
                    "alpha = ", d["alpha"])
        optimal_string = string("optimal sample size",sepstr,
                    "n = ", ceil(Int64,n),sepstr,d["note"])

    # case: Chi Squared
    elseif ht.title == "Chi-square test power calculation"
        n = d["N"]
        n_upper = ceil(Int64, max(n*1.5, n+30)) # upper at least 30 above n
        n_increment = ceil(Int64,(n_upper - 10)/breaks)

        # generate data
        sample_sizes = collect(10:n_increment:n_upper)
        power = [ss->powerChisqTest(N=ss, w=d["w"], alpha = d["alpha"], df=d["df"]) for ss in sample_sizes]

        # create labels
        title = ht.title
        legend_string = string(
                    "df = ", d["df"], sepstr,
                    "effect size w = ", d["w"], sepstr,
                    "alpha = ", d["alpha"])
        optimal_string = string("optimal sample size", sepstr,
                    "N = ", ceil(Int64,n),sepstr,d["note"])

    # case: Normal distribution
    elseif ht.title == "Mean power calculation for normal distribution with known variance"
        n = d["n"]
        n_upper = ceil(Int64, max(n*1.5, n+30)) # upper at least 30 above n
        n_increment = ceil(Int64,(n_upper - 10)/breaks)

        # generate data
        sample_sizes = collect(10:n_increment:n_upper)
        power = [ss->powerNormTest(n=ss, d=d["d"], alpha = d["alpha"], alternative = d["alternative"]) for ss in sample_sizes]

        # create labels
        title = string("Mean power calculation for normal distribution",sepstr,"with known variance")
        legend_string = string(
                    "tails = ", d["alternative"],sepstr,
                    "effect size d = ", d["d"],sepstr,
                    "alpha = ", d["alpha"])
        optimal_string = string("optimal sample size", sepstr,
                    "n = ", ceil(Int64,n))

    # case: proportion
    elseif ht.title == "Proportion power calculation for binomial distribution (arcsine transformation)"
        n = d["n"]
        n_upper = ceil(Int64, max(n*1.5, n+30)) # upper at least 30 above n
        n_increment = ceil(Int64,(n_upper - 10)/breaks)

        # generate data
        sample_sizes = collect(10:n_increment:n_upper)
        power = [ss->powerPTest(n=ss, h=d["h"], alpha = d["alpha"], alternative = d["alternative"]) for ss in sample_sizes]

        # create labels
        title = string("Proportion power calculation",sepstr,"for binomial distribution (arcsine transformation)")
        legend_string = string(
                    "tails = ", d["alternative"],sepstr,
                    "effect size h = ", d["h"],sepstr,
                    "alpha = ", d["alpha"])
        optimal_string = string("optimal sample size",sepstr,
                    "n = ", ceil(Int64,n))

    # case: correlation
    elseif ht.title == "Approximate correlation power calculation (arctangh transformation)"
        n = d["n"]
        n_upper = ceil(Int64, max(n*1.5, n+30)) # upper at least 30 above n
        n_increment = ceil(Int64,(n_upper - 10)/breaks)

        # generate data
        sample_sizes = collect(10:n_increment:n_upper)
        power = [ss->powerRTest(n=ss, r=d["r"], alpha = d["alpha"], alternative = d["alternative"]) for ss in sample_sizes]

        # create labels
        title = string("Approximate correlation power calculation",sepstr,"(arctangh transformation)")
        legend_string = string(
                    "tails = ", d["alternative"],sepstr,
                    "effect size r = ", d["r"],sepstr,
                    "alpha = ", d["alpha"])
        optimal_string = string("optimal sample size",sepstr,
                    "n = ", ceil(Int64,n))
    end

    # use DataFrame
    len = length(sample_sizes)
    if len < 2
        error("cannot display power for this test")
    end

    # select the backend
    backend()

    n_lower = sample_sizes[1]
    xlim_lower = sample_sizes[len] / 80

    # plot with title and x-axis and y-axis labels
    plot(sample_sizes,
        power,
        linecolor = :blue,
        marker=(2,.5,:circle,:blue),
        #left_margin = 8mm,
        top_margin = 8mm,
        right_margin = 4mm,
        title = title,
        yaxis = (ytitle,(-0.03,1.03),0.0:0.2:1.0),
        xaxis = ("Sample Size",(n_lower - xlim_lower,n_upper + xlim_lower)),
        legend=false,
        annotations=([((n_lower+xlim_lower),0.99,text(legend_string,9,:blue,:left,:top)),
            (sample_sizes[floor(Int64,len/2)],0.03,text(optimal_string,9,:red,:left,:bottom))]))
    vline!([n],line=(1,:dot,.8,:red))

end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
##############################################################################
##
## LSMR MultiCores
##
##############################################################################

struct FixedEffectSolverLSMRCores{T} <: AbstractFixedEffectSolver{T}
    fes::Vector{<:FixedEffect}
    weights::AbstractWeights
end

function AbstractFixedEffectSolver{T}(fes::Vector{<:FixedEffect}, weights::AbstractWeights,
                ::Type{Val{:lsmr_cores}}) where {T}
    FixedEffectSolverLSMRCores{T}(fes, weights)
end

function AbstractFixedEffectSolver{T}(fes::Vector{<:FixedEffect}, weights::AbstractVector,
                ::Type{Val{:lsmr_parallel}}) where {T}
    info(":lsmr_parallel is deprecated. Use :lsmr_cores")
    AbstractFixedEffectSolver{T}(fes, weights, Val{:lsmr_cores})
end

function solve_residuals!(r::AbstractVector, feM::FixedEffectSolverLSMRCores{T}; kwargs...) where {T}
    solve_residuals!(r, AbstractFixedEffectSolver{T}(feM.fes, feM.weights, Val{:lsmr}); kwargs...)
end

function solve_residuals!(X::AbstractMatrix, feM::FixedEffectSolverLSMRCores; kwargs...)
    iterations = zeros(Int, size(X, 2))
    convergeds = zeros(Bool, size(X, 2))
    result = pmap(x -> solve_residuals!(x, feM; kwargs...), [X[:, j] for j in 1:size(X, 2)])
    for j in 1:size(X, 2)
        X[:, j] = result[j][1]
        iterations[j] =  result[j][2]
        convergeds[j] = result[j][3]
    end
    return X, iterations, convergeds
end

function solve_coefficients!(r::AbstractVector, feM::FixedEffectSolverLSMRCores{T}, ; kwargs...) where {T}
    solve_coefficients!(r, AbstractFixedEffectSolver{T}(feM.fes, feM.weights, Val{:lsmr}); kwargs...)
end

##############################################################################
##
## LSMR MultiThreads
##
## One needs to construct a new fe matrix / fe vectirs for each LHS/RHS
## This is because each thread need to have its own cache / u / w, etc
##
##############################################################################

struct FixedEffectSolverLSMRThreads{T} <: AbstractFixedEffectSolver{T}
    x::Vector{FixedEffectSolverLSMR{T}}
end

function AbstractFixedEffectSolver{T}(fes::Vector{<:FixedEffect}, weights::AbstractWeights, ::Type{Val{:lsmr_threads}}) where {T}
    FixedEffectSolverLSMRThreads([AbstractFixedEffectSolver{T}(fes, weights, Val{:lsmr}) for _ in 1:Threads.nthreads()])
end

function solve_residuals!(r::AbstractVector, feM::FixedEffectSolverLSMRThreads; kwargs...)
    solve_residuals!(r, feM.x[1]; kwargs...)
end

function solve_residuals!(X::AbstractMatrix, feM::FixedEffectSolverLSMRThreads; kwargs...)
   iterations = zeros(Int, size(X, 2))
   convergeds = zeros(Bool, size(X, 2))
   Threads.@threads for j in 1:size(X, 2)
        _, iteration, converged = solve_residuals!(view(X, :, j), feM.x[Threads.threadid()]; kwargs...)
        iterations[j] = iteration
        convergeds[j] = converged
   end
   return X, iterations, convergeds
end

function solve_coefficients!(r::AbstractVector, feM::FixedEffectSolverLSMRThreads; kwargs...)
    solve_coefficients!(r, feM.x[1]; kwargs...)
end


text:
using StaticArrays
using LinearAlgebra
import ElectricFields: compute_rotation, rotation_angle, rotation_axis

@testset "Rotations" begin
    @test compute_rotation(nothing) == I
    @test compute_rotation(π*I) == I

    @testset "From axis–angle" begin
        @test compute_rotation((0,[0,0,1])) == I
        @test compute_rotation((π/4,[0,0,1])) ≈ [1/√2 -1/√2 0
                                                 1/√2 1/√2  0
                                                 0    0     1]
        @test compute_rotation((π/2,[0,0,1])) ≈ [0 -1 0
                                                 1  0 0
                                                 0  0 1]
        @test compute_rotation((π,[0,0,1])) ≈ [-1 0 0
                                               0 -1 0
                                               0  0 1]
        @test compute_rotation((π/2,[0,1,1])) ≈ [0    -1/√2 1/√2
                                                 1/√2  1/2  1/2
                                                 -1/√2 1/2  1/2]
    end

    @testset "From rotation matrix" begin
        @test compute_rotation([1.0 1.0 0.0
                                0.0 1.0 1.0
                                0.0 0.0 1.0]) ≈ I rtol=1e-14
        @test compute_rotation([1.0 1.0 1.0
                                0.0 1.0 1.0
                                0.0 0.0 1.0]) ≈ I rtol=1e-14

        @test_throws DimensionMismatch compute_rotation([1.0 0.0
                                                         0.0 1.0])

        @test_throws ArgumentError compute_rotation([1.0 1.0 0.0
                                                     0.0 0.0 0.0
                                                     0.0 0.0 1.0])
    end

    @testset "Rotation axis and angle from matrix" begin
        @test rotation_angle(Matrix(I, 3, 3)) == 0
        for (angle,axis) in [(π/4,[0,0,1]),
                             (π/3,[0,0,1]),
                             (π/2,[0,0,1]),
                             (π,[0,0,1]),
                             (π/2,[0,1,1])]
            R = compute_rotation((angle, axis))
            @test rotation_angle(R) ≈ angle
            @test rotation_axis(R) ≈ normalize(axis)
        end
    end
end


text:
type Cursor{T<:Tuple} <: IO
    stmt::Stmt
    ncols::Int
    cols::Vector{UTF8String}
    status::Cint
    where::String
end

function Cursor(stmt::Stmt, where::String="")
    status = Nierika.execute!(stmt)
    ncols::Int = sqlite3_column_count(stmt.handle)
    types = coltypes(stmt, ncols)
    cols = colnames(stmt, ncols)
    return Cursor{types}(stmt, ncols, cols, status, where)
end

function colnames(curs::Cursor)
    n = curs.ncols
    colnames = [ utf8(sqlite3_column_name(curs.stmt.handle, i-1)) for i in 1:n ]
end

function colnames(stmt::Stmt, ncols=sqlite3_column_count(stmt.handle))
    colnames = [ utf8(sqlite3_column_name(stmt.handle, i-1)) for i in 1:ncols ]
end

function coltypes(stmt::Stmt, ncols::Int)
    types = Array{DataType, 1}()
    for i in 1:ncols
        typ = sqlite3_types[sqlite3_column_type(stmt.handle, i-1)]
        push!(types, Vector{typ})
    end
    return Tuple{types...}
end

Base.close(curs::Cursor) = close(curs.stmt)

function Base.eof(curs::Cursor)
    (curs.status == SQLITE_DONE || curs.status == SQLITE_ROW) || sqliteerror(curs.stmt.db)
    return curs.status == SQLITE_DONE
end

function Base.readline(curs::Cursor, delim::Union{Char,ByteString}=", ", buf::IOBuffer=IOBuffer())
    eof(curs) && return ""
    for i = 1:curs.ncols
        val = sqlite3_column_text(curs.stmt.handle, i-1)
        val != C_NULL && write(buf, bytestring(val))
        write(buf, ifelse(i == curs.ncols, '\n', delim))
    end
    # curs.status = sqlite3_step(curs.stmt.handle)
    return takebuf_string(buf)
end

function step!(curs::Cursor)
    curs.status = execute!(curs.stmt)
    return curs
end

query(db::DB, sql::String) = Cursor(Stmt(db, sql))
function query(sql::String)
    actuality.actualized ? Cursor(Stmt(actuality.db, sql)) : error()
end

function restof(curs::Cursor; delim=',', buf::IOBuffer=IOBuffer())
    while curs.status == 100
        for i = 1:curs.ncols
            val = sqlite3_column_text(curs.stmt.handle, i-1)
            val != C_NULL && write(buf, bytestring(val))
            write(buf, ifelse(i == curs.ncols, '\n', delim))
        end
        step!(curs)
    end
    return seekstart(buf)
end

function results(curs::Cursor)
    return readdlm(restof(curs), ',', '\n')
end

reset!(stmt::Stmt) = sqlite3_reset(stmt.handle)
function reset!(curs::Cursor)
    reset!(curs.stmt)
    curs.status = execute!(curs.stmt)
    return curs
end

function nrows(curs::Cursor)
    return sqlite3_column_int64(curs.stmt.handle, 0)
end

let safe, buf
    safe = false
    buf = IOBuffer()
    where_buf = IOBuffer()

    global select
    function Base.select(f, db::DB, table::String; sel::String="*")
        safe = true
        write(buf, "select $sel from $table")
        f()
        _where = takebuf_string(where_buf)
        write(buf, _where)
        sql = takebuf_string(buf)
        safe = false
        return Cursor(Stmt(db, sql, table), _where)
        # println(sql)
    end

    global where
    function where(sql::String)
        safe == true || error()
        write(where_buf, " where $sql")
    end
end


text:
module ChanceConstraintTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable
using ForneyLab: SPChanceConstraintOutG, truncatedGaussianMoments

#--------
# Helpers
#--------

@testset "truncatedGaussianMoments" begin
    @test truncatedGaussianMoments(0.0, 1.0, 0.0, Inf) == (0.5, 0.7978845608028654, 0.3633802276324186)
    @test truncatedGaussianMoments(0.0, 1.0, -Inf, 0.0) == (0.5, -0.7978845608028654, 0.3633802276324186)
    @test truncatedGaussianMoments(0.0, 1.0, -Inf, Inf) == (1.0, 0.0, 1.0)
    @test truncatedGaussianMoments(0.0, 1.0, 0.0, 0.0) == (0.0, 0.0, 0.0) # Improper truncation
end

#-------------
# Update rules
#-------------

@testset "SPChanceConstraintOutG" begin
    @test SPChanceConstraintOutG <: SumProductRule{ChanceConstraint}
    @test outboundType(SPChanceConstraintOutG) == Message{GaussianWeightedMeanPrecision}
    @test isApplicable(SPChanceConstraintOutG, [Nothing])

    @test ruleSPChanceConstraintOutG(Message(Univariate, GaussianMeanVariance, m=0.0, v=1.0), (0.0, Inf), 0.05) == Message(Univariate, GaussianWeightedMeanPrecision, xi=2.9470311004635064, w=2.2101933871827293)
    @test ruleSPChanceConstraintOutG(Message(Univariate, GaussianMeanVariance, m=0.0, v=1.0), (-Inf, Inf), 0.05) == Message(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=tiny)
end

@testset "averageEnergy" begin
    @test averageEnergy(ChanceConstraint, ProbabilityDistribution(GaussianMeanVariance, m=0.0, v=1.0)) == 0.0
end

#------------
# Integration
#------------

@testset "ChanceConstraint integration" begin
    fg = FactorGraph()
    @RV x ~ GaussianMeanVariance(0.0, 1.0)
    ChanceConstraint(x; G=(0.0, Inf), epsilon=0.05)

    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(x)
    algo_code = algorithmSourceCode(algo)

    @test occursin("ruleSPChanceConstraintOutG(messages[1], (0.0, Inf), 0.05)", algo_code)
end

end # module


text:
function maketinfoplot(cinfo,x,y)
tinfo = Tinfo()

@unpack en,P,M,Ω,basis = cinfo
ωx,ωy=Ω; Mx,My=M;
Tx = eigmat(Mx,x,ω=ωx/ωy)
Ty = eigmat(My,y,ω=ωy/ωy)
@pack tinfo = Tx,Ty

return tinfo
end


text:
# v0.6

# Task 1
function hofstadterconway(n::Integer)::Vector{Int}
    rst = fill(1, n)
    for i in 3:n
        rst[i] = rst[rst[i - 1]] + rst[i - rst[i - 1]]
    end
    return rst
end

function hcfraction(n::Integer)::Vector{Float64}
    rst = Array{Float64}(n)
    for (i, a) in enumerate(hofstadterconway(n))
        rst[i] = abs(a) / i
    end
    return rst
end

# Task 2
seq = hcfraction(2 ^ 20)
for i in 1:20
    a, b = 2 ^ (i - 1), 2 ^ i
    @printf("max value of a(n)/n in %6i < n < %7i = %5.3f\n", a, b, maximum(seq[a:b]))
end

# Task 3
function lastindex(val::Float64)
    r, p = 1, 0
    # Find the range of 2 power in which the maximum is < val
    seq = hcfraction(2 ^ 15)
    while maximum(seq[2^p:2^(p+1)]) > val; p += 1 end
    r = 2 ^ (p + 1)
    while seq[r] < val; r -= 1 end
    return r
end

println("You too might have won \$1000 with the mallows number of ", lastindex(0.55))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Distributions
# cdf() is equivalent to p[dist] in R
# quantile() is equivalent to q[dist] in R

for n=6:250
    NCP = (88-81)^2/(22.6^2*(2/n))
    power_calc = 1-cdf(NoncentralF(1,n*4-5,NCP),quantile(FDist(1,n*4-5),0.95))
    println(n*2)
    println(power_calc)
end

max_stars_repo_name: dshehtanian/julia-learn

max_issues_repo_name: dshehtanian/julia-learn

max_forks_repo_name: dshehtanian/julia-learn


text:
module util

export Point, Point3D, extract_lat, extract_lon, latlon, format_date,
format_datetime, format_time, Point_deg

import Base:(convert)
using Dates
using Navigation: Point_deg, Point_rad

"Point type with latitude `lat` [deg] and longitude `lon` [deg]"
struct Point{T<:Float64}
    lat::T
    lon::T
end

"Point type with latitude `lat` [deg], longitude `lon` [deg] and altitude [m]"
struct Point3D{T<:Float64}
    lat::T
    lon::T
    alt::T
end

Base.:convert(type::Type{Point_deg}, x::Point) = Point_deg(x.lat, x.lon)
Base.:convert(type::Type{Point_rad}, x::Point) = Point_rad(deg2rad(x.lat), deg2rad(x.lon))
Base.:convert(type::Type{Point}, x::Point_deg) = Point(x.ϕ, x.λ)
Base.:convert(type::Type{Point}, x::Point_rad) = Point(rad2deg(x.ϕ), rad2deg(x.λ))

function extract_lat(str::AbstractString)
    if length(str) > 2
        if str[1] in ['N', 'S']
            sign_lat = str[1] == 'N' ? 1 : -1
            relative = 1
        elseif str[end] in ['N', 'S']
            sign_lat = str[end] == 'N' ? 1 : -1
            relative = 0
        else
            return NaN
        end
    else
        return NaN
    end

    lat_h = parse(Float64, str[relative+1:relative+2])
    lat_m = 0.0
    lat_s = 0.0
    if length(str) > 3
        lat_m = parse(Float64, str[relative+3:relative+4])
    end
    if length(str) > 5
        lat_s = parse(Float64, str[relative+5:end+relative-1])
    end
    return sign_lat*(lat_h + lat_m/60.0 + lat_s/3600.0)
end

function extract_lon(str::AbstractString)
    if length(str) > 3
        if str[1] in ['E', 'W']
            sign_lon = str[1] == 'E' ? 1 : -1
            relative = 1
        elseif str[end] in ['E', 'W']
            sign_lon = str[end] == 'E' ? 1 : -1
            relative = 0
        else
            return NaN
        end
    else
        return NaN
    end

    lon_h = parse(Float64, str[relative+1:relative+3])
    lon_m = 0.0
    lon_s = 0.0
    if length(str) > 4
        lon_m = parse(Float64, str[relative+4:relative+5])
    end
    if length(str) > 6
        lon_s = parse(Float64, str[relative+6:end+relative-1])
    end
    return sign_lon*(lon_h + lon_m/60.0 + lon_s/3600.0)
end

function latlon(str_lat::AbstractString, str_lon::AbstractString)
    lat = extract_lat(str_lat)
    lon = extract_lon(str_lon)
    if !isnan(lat*lon)
        # return Point(lat, lon)
        return Point_deg(lat, lon)
    else
        return NaN
    end
end

function latlon(str::AbstractString)
    nors = "S"
    if occursin("S", str)
        nors = "S"
    else
        nors = "N"
    end

    split_str = split(str, nors)
    str_lat = split_str[1]*nors
    str_lon = split_str[2]
    return latlon(str_lat, str_lon)
end

function format_date(str::Union{AbstractString, Missing}, format; addyear = Year(0))
    if str === missing
        return missing
    elseif str == "!!!!"
        return missing
    end
    maybedate = tryparse(Date, str, format)
    if maybedate == nothing
        return missing
    elseif maybedate > Date(1, 1, 1)
        return Date(str, format) + addyear
    else
        return missing
    end
end

function format_datetime(str::Union{AbstractString, Missing}, format; addyear = Year(0))
    if str === missing
        return missing
    elseif occursin("!!!!", str)
        return missing
    end
    maybedatetime = tryparse(DateTime, str, format)
    if maybedatetime == nothing
        return missing
    elseif maybedatetime > DateTime(1, 1, 1, 0, 0, 0)
        return DateTime(str, format) + addyear
    else
        return missing
    end
end

function format_time(str::Union{AbstractString, Missing}, format)
    if str === missing
        return missing
    end

    # handling special case of 2400 hours used in ATM to indicate 0000 24 hours later
    if str == "2400"
        str = "235959"
        format = DateFormat("HHMMSS")
    end

    maybetime = tryparse(Time, str, format)
    if maybetime == nothing
        return missing
    else
        if str == ""
            return missing
        else
            return Time(str, format)
        end
    end
end

end  #Module


text:
# Distance dendrogram

using JLD2, DataFrames, CSV, Dates, Query
using LinearAlgebra, Clustering, StatsPlots

include("/home/kristen/Documents/coda_utilities/coda_ops_transforms.jl")

#load in the data:
@load "/home/kristen/Documents/V6V8_analysis/analysis_products/oligo_df_zero_plus.jld2"
#should return dataframe oligo_df_plus

## closure on the dataset:
y=closure(convert(Array,oligo_df_plus[!,Between(:csO1,:csO6)])) #seasonally imputed zeros
# findall(y .== 0) #this should be 0 as a double check!
vyp=ait_var(y) #compute aitchison variation

#other metrics to consider:
# test=exp.(-vyp.^2 / 2)
# test=exp.(-sqrt.(vyp))

#visualize! but must be with Plots.jl!
hc = hclust(Symmetric(ait_var(y)), linkage=:complete, uplo=:U) #:complete
#-> a few other linkage options, check hclust help

##
w=plot(hc,xticks=false,linewidth=2)
node_labels=["O1-I","O2-CB5","O3-XV","O4-III/IV","O5-I","O6-I*"] #oligotypes 1-6
#node_labels[hc.order]
xticks!(w,([1,2,3,4,5,6],node_labels[hc.order]),fontsize=16) #maps from hc and matching up oligos to clades
ylabel!(w,"Distance, Aitchision variation")
##
savefig(w, "/home/kristen/Documents/V6V8_analysis/paper_V6V8/V6V8_Syn_CoDa_paper/figures/dendrogram_complete.pdf")


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

"""

    AbstractMenu

The supertype for all Menu types.
See AbstractMenu.jl for descriptions of functions mentioned in this
doc string.

# Functions

The following functions can be called on all <:AbstractMenu types.
Details can be found in

## Exported

  - `request(m::AbstractMenu)`
  - `request(msg::AbstractString, m::AbstractMenu)`

## Hidden

  - `printMenu(m::AbstractMenu, cursor::Int; init::Bool=false)`

# Subtypes

All subtypes must contain the fields `pagesize::Int` and
`pageoffset::Int`. They must also implement the following functions.

## Necessary Functions

These functions must be implemented for all subtypes of AbstractMenu.

  - `pick(m::AbstractMenu, cursor::Int)`
  - `cancel(m::AbstractMenu)`
  - `options(m::AbstractMenu)`
  - `writeLine(buf::IOBuffer, m::AbstractMenu, idx::Int, cur::Bool)`

## Optional Functions

These functions do not need to be implemented for all AbstractMenu
subtypes.

  - `header(m::AbstractMenu)`
  - `keypress(m::AbstractMenu, i::UInt32)`

"""
abstract type AbstractMenu end

# NECESSARY FUNCTIONS
# These functions must be implemented for all subtypes of AbstractMenu
######################################################################

"""
    pick(m::AbstractMenu, cursor::Int)

Defines what happens when a user presses the Enter key while the menu is open.
If `true` is returned, `request()` will exit.
"""
pick(m::AbstractMenu, cursor::Int) = error("unimplemented")

"""
    cancel(m::AbstractMenu)

Define what happens when a user cancels ('q' or ctrl-c) a menu.
`request()` will always exit after calling this function.
"""
cancel(m::AbstractMenu) = error("unimplemented")

"""
    options(m::AbstractMenu)

Return a list of strings to be displayed as options in the current page.
"""
options(m::AbstractMenu) = error("unimplemented")

"""
    writeLine(buf::IOBuffer, m::AbstractMenu, idx::Int, cur::Bool)

Write the option at index `idx` to the buffer. If cursor is `true` display the cursor.
"""
function writeLine(buf::IOBuffer, m::AbstractMenu, idx::Int, cur::Bool)
    error("unimplemented")
end

# OPTIONAL FUNCTIONS
# These functions do not need to be implemented for all menu types
##################################################################

"""
    header(m::AbstractMenu)

Displays the header above the menu when it is rendered to the screen.
"""
header(m::AbstractMenu) = ""

"""
    keypress(m::AbstractMenu, i::UInt32)

Send any non-standard keypress event to this function.
If `true` is returned, `request()` will exit.
"""
keypress(m::AbstractMenu, i::UInt32) = false

"""
    request(m::AbstractMenu)

Display the menu and enter interactive mode. Returns `m.selected` which
varies based on menu type.
"""
request(m::AbstractMenu) = request(terminal, m)

function request(term::REPL.Terminals.TTYTerminal, m::AbstractMenu)
    cursor = 1

    menu_header = header(m)
    !CONFIG[:supress_output] && menu_header != "" && println(term.out_stream, menu_header)

    printMenu(term.out_stream, m, cursor, init=true)

    raw_mode_enabled = enableRawMode(term)
    raw_mode_enabled && print(term.out_stream, "\x1b[?25l") # hide the cursor

    lastoption = length(options(m))
    try
        while true
            c = readKey(term.in_stream)

            if c == Int(ARROW_UP)
                if cursor > 1
                    cursor -= 1 # move selection up
                    if cursor < (2+m.pageoffset) && m.pageoffset > 0
                        m.pageoffset -= 1 # scroll page up
                    end
                elseif CONFIG[:scroll_wrap] # wrap to bottom
                    cursor = lastoption
                    m.pageoffset = lastoption - m.pagesize
                end

            elseif c == Int(ARROW_DOWN)
                if cursor < lastoption
                    cursor += 1 # move selection down
                    if m.pagesize + m.pageoffset <= cursor < lastoption
                        m.pageoffset += 1 # scroll page down
                    end
                elseif CONFIG[:scroll_wrap] # wrap to top
                    cursor = 1
                    m.pageoffset = 0
                end

            elseif c == Int(PAGE_UP)
                # If we're at the bottom, move the page 1 less to move the cursor up from
                # the bottom entry, since we try to avoid putting the cursor at bounds.
                m.pageoffset -= m.pagesize - (cursor == lastoption ? 1 : 0)
                m.pageoffset = max(m.pageoffset, 0)
                cursor -= m.pagesize
                cursor = max(cursor, 1)

            elseif c == Int(PAGE_DOWN)
                m.pageoffset += m.pagesize - (cursor == 1 ? 1 : 0)
                m.pageoffset = min(m.pageoffset, lastoption - m.pagesize)
                cursor += m.pagesize
                cursor = min(cursor, lastoption)

            elseif c == Int(HOME_KEY)
                cursor = 1
                m.pageoffset = 0

            elseif c == Int(END_KEY)
                cursor = lastoption
                m.pageoffset = lastoption - m.pagesize

            elseif c == 13 # <enter>
                # will break if pick returns true
                pick(m, cursor) && break

            elseif c == UInt32('q')
                cancel(m)
                break

            elseif c == 3 # ctrl-c
                cancel(m)
                CONFIG[:ctrl_c_interrupt] ? throw(InterruptException()) : break

            else
                # will break if keypress returns true
                keypress(m, c) && break
            end

            printMenu(term.out_stream, m, cursor)
        end
    finally # always disable raw mode
        if raw_mode_enabled
            print(term.out_stream, "\x1b[?25h") # unhide cursor
            disableRawMode(term)
        end
    end
    println(term.out_stream)

    return m.selected
end

"""
    request([term,] msg::AbstractString, m::AbstractMenu)
Shorthand for `println(msg); request(m)`.
"""
request(msg::AbstractString, m::AbstractMenu) = request(terminal, msg, m)

function request(term::REPL.Terminals.TTYTerminal, msg::AbstractString, m::AbstractMenu)
    println(term.out_stream, msg)
    request(term, m)
end

"""
    printMenu(out, m::AbstractMenu, cursor::Int; init::Bool=false)

Display the state of a menu.
"""
function printMenu(out, m::AbstractMenu, cursor::Int; init::Bool=false)
    CONFIG[:supress_output] && return

    buf = IOBuffer()

    lines = m.pagesize-1

    if init
        m.pageoffset = 0
    else
        # Move the cursor to the beginning of where it should print
        print(buf, "\x1b[999D\x1b[$(lines)A")
    end

    firstline = m.pageoffset+1
    lastline = m.pagesize+m.pageoffset

    for i in firstline:lastline
        # clearline
        print(buf, "\x1b[2K")

        if i == firstline && m.pageoffset > 0
            print(buf, CONFIG[:up_arrow])
        elseif i == lastline && i != length(options(m))
            print(buf, CONFIG[:down_arrow])
        else
            print(buf, " ")
        end

        writeLine(buf, m, i, i == cursor)

        i != lastline && print(buf, "\r\n")
    end

    print(out, String(take!(buf)))
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
PROJECT_IMG_FOLDER = joinpath("data", "book-covers")
DATASETS = ["test", "train"]

if ~isdir(PROJECT_IMG_FOLDER) mkdir(PROJECT_IMG_FOLDER) end

for dataset_name in DATASETS

    lines = nothing

    current_folder = joinpath(PROJECT_IMG_FOLDER, dataset_name)
    if ~isdir(current_folder) mkdir(current_folder) end

    file_name = joinpath(pwd(), "book30-listing-$dataset_name.csv")
    lines = readlines(file_name);

    for line in lines
        line_split = split(line, "\",\"")
        download(line_split[3], joinpath(PROJECT_IMG_FOLDER, dataset_name, line_split[2]))
    end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export badge

"""
    badge(;kwargs...)
    badge(children::Any;kwargs...)
    badge(children_maker::Function;kwargs...)

A Badge component.
Display badge, pill or tag. For more information, see: https://mantine.dev/core/badge/
Keyword arguments:
- `children` (String; optional): Primary content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `color` (a value equal to: "dark", "gray", "red", "pink", "grape", "violet", "indigo", "blue", "cyan", "teal", "green", "lime", "yellow", "orange"; optional): Badge color from theme
- `fullWidth` (Bool; optional): Sets badge width to 100% of parent element, hides overflow text with text-overflow: ellipsis
- `gradient` (optional): Controls gradient settings in gradient variant only. gradient has the following type: lists containing elements 'from', 'to', 'deg'.
Those elements have the following types:
  - `from` (String; required)
  - `to` (String; required)
  - `deg` (Real; optional)
- `radius` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Predefined border-radius value from theme.radius or number for border-radius in px
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Predefined badge size
- `style` (Dict; optional): Inline style override
- `variant` (a value equal to: "light", "filled", "outline", "dot", "gradient"; optional): Controls badge background, color and border styles
"""
function badge(; kwargs...)
        available_props = Symbol[:children, :id, :className, :color, :fullWidth, :gradient, :radius, :size, :style, :variant]
        wild_props = Symbol[]
        return Component("badge", "Badge", "dash_mantine_components", available_props, wild_props; kwargs...)
end

badge(children::Any; kwargs...) = badge(;kwargs..., children = children)
badge(children_maker::Function; kwargs...) = badge(children_maker(); kwargs...)


text:
"""
    B,T,Q,R = lll(H,δ=3/4)

Do Lenstra–Lenstra–Lovász lattice reduction of matrix `H` using optional
parameter `δ`.  The output is `B`, an LLL-reduced basis; `T`, a unimodular
(meaning `det(T)=+/-1`) transformation matrix such that `B= H*T`; and
finally `Q` and `R` which are a QR decomposition of `B`.  So `H = B*inv(T) =
Q*R*inv(T)`.

Follows D. Wuebben, et al, "Lattice Reduction - A Survey with Applications
in Wireless Communications". IEEE Signal Processing Magazine, 2011. When
comparing with the core lattice reduction technique, we belive it is closest
to the floating-point algorithm of C. P. Schnorr. "A more efficient
algorithm for lattice basis reduction". Journal of Algorithms, Vol 9, 1988.

# Examples
```jldoctest
julia> H= [1 2; 3 4];B,_ = lll(H); B
2×2 Array{Int64,2}:
 1  -1
 1   1

julia> H= BigFloat.([1.5 2; 3 4]) .+ 2im; B,_= lll(H); B
2×2 Array{Complex{BigFloat},2}:
 0.50+0.0im  0.0+1.0im
  1.0+0.0im  0.0+0.0im

julia> N=500;H = randn(N,N); B,T = lll(H);

```
"""
function lll(H::AbstractArray{Td,2},δ::Float64=3/4) where {Td<:Number}

if !(0.25 < δ < 1.0)
    error("δ must be between 1/4 and 1.");
end
if Td<:AbstractFloat
    δ= Td(δ)
end

B = copy(H);
L = size(B,2);
Qt,R = qr(B);
Q = Matrix(Qt); # A few cycles can be saved by skipping updates of the Q matrix.

Ti= getIntType(Td)
T = Matrix{Ti}(I, L, L)
zeroTi = real(zero(Ti))

lx  = 2;
while lx <= L

    # reduce lx-th column of B
    for k=lx-1:-1:1
        rk = R[k,lx]/R[k,k]
        mu = roundf(rk)
        if abs(mu)>zeroTi
            # B[:,lx]   -= mu * B[:,k]
            # R[1:k,lx] -= mu * R[1:k,k]
            # T[:,lx]   -= mu * T[:,k]
            B[:,lx]   .-= mu .* view(B,:,k)
            R[1:k,lx] .-= mu .* view(R,1:k,k)
            T[:,lx]   .-= mu .* view(T,:,k)
        end
    end

    nrm = norm(R[lx-1:lx,lx])
    if δ*abs(R[lx-1,lx-1])^2 > nrm^2

        # swap columns lx-1 and lx in B, T and R
        # B[:,[lx-1,lx]]    = B[:,[lx,lx-1]]
        # T[:,[lx-1,lx]]    = T[:,[lx,lx-1]]
        # R[1:lx,[lx-1,lx]] = R[1:lx,[lx,lx-1]]
        B[:,[lx-1,lx]]    .= view(B,:,[lx,lx-1])
        T[:,[lx-1,lx]]    .= view(T,:,[lx,lx-1])
        R[1:lx,[lx-1,lx]] .= view(R,1:lx,[lx,lx-1])

        # upper triangular by Givens rotation
        cc = R[lx-1,lx-1] / nrm # nrm = ||R[lx-1:lx,lx-1]|| after swapping
        ss = R[lx,lx-1]   / nrm
        Θ = [cc' ss; -ss cc]

        # the following multiply by Θ results in R[lx,lx-1] = 0
        R[lx-1:lx,lx-1:end] .= Θ * R[lx-1:lx,lx-1:end]
        # Q[:,lx-1:lx] = Q[:,lx-1:lx] * Θ'
        Q[:,lx-1:lx] .= view(Q,:,lx-1:lx) * Θ'
        lx = max(lx-1,2)
    else
        lx = lx+1;
    end
end
return B,T,Q,R
end

"""
    B = gauss(H)

Do Gauss/Lagrange reduction on the lattice defined by the two columns of
H.

Follows Fig 2.3 of "Lattice Basis Reduction: An Introduction to the LLL
Algorithm and Its Applications" by Murray R. Bremner, CRC Press, 2012.

# Examples
```jldoctest
julia> H = [1 2; 3 3]; B = gauss(H)
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  3.0

```
"""
function gauss(H::AbstractArray{Td,2}) where Td
    if size(H,2)!=2
        error("Gauss reduction only works on two columns.")
    end
    x = float.(H[:,1])
    y = float.(H[:,2])
    if norm(x)<=norm(y)
        v1,v2=x,y
    else
        v1,v2=y,x
    end
    finished = false;
    while !finished
        m = round((v2'*v1)/(v1'*v1))
        v2 -= m*v1
        if norm(v1)<=norm(v2)
            finished = true
        else
            v1,v2 = v2,v1
        end
    end
    return [v1 v2]
end

"""
    getIntType(Td)

    Return an integer type that matches Td.  I.e. BigInt for Td=BigFloat, Int64
    for Td=Float64.

    It's possible that the integer type could be picked on a per-matrix
    basis depending on the values in the matrix and the decomposition
    (i.e. [`lll`](@ref)) that is to be done. This is not done below.
# Examples
```julia-repl
julia> Pkg.add("DoubleFloats")
julia> using DoubleFloats
julia> using LLLplus
julia> import LLLplus.getIntType # yes, it's type piracy!
julia> x = randn(3,3); xd = Double64(x); lll(xd)
ERROR: MethodError: no method matching getIntType(::Type{DoubleFloat{Float64}})
...
julia> getIntType(Td::Type{Tr}) where {Tr<:Double64} = Int128
julia> B,T =lll(xd); # succeeds
```
"""
getIntType(Td::Type{Tr}) where {Tr<:Integer} = Tr
getIntType(Td::Type{Tr}) where {Tr<:Float64} = Int64
getIntType(Td::Type{Tr}) where {Tr<:Float32} = Int32
getIntType(Td::Type{Tr}) where {Tr<:Float16} = Int16
getIntType(Td::Type{Tr}) where {Tr<:BigFloat} = BigInt
getIntType(Td::Type{Complex{Tr}}) where Tr = Complex{getIntType(Tr)}

"""
    B,T,Q,R = sizereduction(H)
Do size reduction of matrix H.  The output is B, an size-reduced
basis; T, a unimodular (det(T)=+/-1) transformation matrix such that B=
H*T; Q and R such that B=Q*R and Q is orthonormal, and R is upper
triangular.  So H = Q*R*inv(T)

H can be of Integer, FloatingPoint, BigInt, or BigFloat types. The core
algorithm is designed for floating-point.

Size reduction is the first part of LLL lattice reduction.

# Examples
```jldoctest
julia> H= [1 2; 3 4];B,_ = sizereduction(H); B
2×2 Array{Int64,2}:
 3   1
 1   1

julia> H= BigFloat.([1.5 2; 3 4]) .+ 2im; B,_= sizereduction(H); B
2×2 Array{Complex{BigFloat},2}:
 1.5+2.0im  0.50+0.0im
 3.0+2.0im   1.0+0.0im

julia> N=100;H = randn(N,N); B,T = sizereduction(H);

```
"""
function sizereduction(H::Array{Td,2}) where {Td}
    B = copy(H);
    L = size(B,2);
    Q,R = qr(B);

    Ti= getIntType(Td)
    T = Matrix{Ti}(I, L, L)

    for lx=2:L
        # reduce lx-th column of B
        for k=lx-1:-1:1
            rk = R[k,lx]/R[k,k]
            mu = roundf(rk)
            if abs(mu)>0
                # B[:,lx]   -= mu * B[:,k]
                # R[1:k,lx] -= mu * R[1:k,k]
                # T[:,lx]   -= mu * T[:,k]
                B[:,lx]   .-= mu .* view(B,:,k)
                R[1:k,lx] .-= mu .* view(R,1:k,k)
                T[:,lx]   .-= mu .* view(T,:,k)
            end
        end
    end

    return B,T,Q,R
end


text:
mutable struct Mapping

    key::String
    value::Set{String}

    function Mapping()
		this = new()
	end
end

mutable struct CBMetabolite

    symbol::String
    name::String
    kegg_id::String

    function CBMetabolite()
		this = new()
	end
end

mutable struct CBRule

    reaction_name::String
    rule_text::String

    function CBRule()
		this = new()
	end
end

mutable struct CBReactionGeneMap

    reaction_name::String
    genes::Set{String}

    function CBReactionGeneMap()
        this = new()
    end
end

mutable struct CBMetabolicReaction

    record::String
    reaction_name::String
    ec_number::String
    left_phrase::String
    right_phrase::String
    reverse::String
    forward::String

    function CBMetabolicReaction()
		this = new()
	end
end

max_stars_repo_name: varnerlab/CBModelTools

max_issues_repo_name: varnerlab/CBModelTools

max_forks_repo_name: varnerlab/CBModelTools


text:
module Client
using HTTP

# Explicitly use a client whose chunksize is well above our chunk size so that
# we don't accidentally trigger chunked encoding
client = HTTP.Client(chunksize=2*1024*1024)
res = nothing

function upload_chunk(file::IOStream, target::AbstractString,
                      chunk_idx::Integer, chunk_size::Integer)

function upload_file(file_path::AbstractString, target::AbstractString;
                     chunk_size::Integer = 1*1024*1024,
                     mime_type::MIME = MIME(""))
    global client, res

    filename = basename(file_path)
    total_size = filesize(file_path)
    num_chunks = ceil(Int,total_size/chunk_size)

    f = open(file_path, "r")

    for chunk_idx = 1:num_chunks
        params = Dict(
            "resumableChunkNumber" => chunk_idx,
            "resumableChunkSize" => chunk_size,
            "resumableCurrentChunkSize" => chunk_size,
            "resumableTotalSize" => total_size,
            "resumableType" => string(mime_type),
            "resumableIdentifier" => "$(total_size)-$(filename)",
            "resumableFileName" => filename,
            "resumableRelativePath" => filename,
            "resumableTotalChunks" => num_chunks,
        )

        # Check each chunk before uploading the full chunk
        uri = HTTP.URI(target; query=params)
        res = HTTP.get(client, uri)
        if res.status != 200
            # That chunk is not there yet, let's post it up.
            seek(f, (chunk_idx - 1)*chunk_size)
            data = read(f, chunk_size)
            params["file"] = IOBuffer(data)

            println("Uploading chunk $chunk_idx...")
            res = HTTP.post(client, uri; body=params)
            if res.status != 200
                error("Upload failed with HTTP code $(res.status)")
            end
        end
    end
end

end #module client


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

choose the frequency of experiments to minimize a given convex spectral function
of the information matrix and satisfy an experiment budget constraint and other
side constraints

adapted from Boyd and Vandenberghe, "Convex Optimization", section 7.5

minimize    f(V × Diagonal(x) × V')
subject to  x ≥ 0
            Σ x = k
            A x = b
where variable x ∈ ℝᵖ is the frequency of each experiment, k is the number of
experiments to run (let k = p), the columns of V ∈ ℝ^(q × p) correspond to each
experiment (let q ≈ p/2), and f is a convex spectral function
=#

struct ExperimentDesignJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    p::Int
    ext::MatSpecExt # formulation specifier
end

function build(inst::ExperimentDesignJuMP{T}) where {T <: Float64}
    p = inst.p
    @assert p >= 2
    @assert is_domain_pos(inst.ext)

    q = div(p, 2)
    V = randn(T, q, p)
    A = randn(T, round(Int, sqrt(p - 1)), p)
    b = sum(A, dims = 2)

    model = JuMP.Model()
    JuMP.@variable(model, x[1:p] >= 0)
    JuMP.@constraint(model, sum(x) == p)
    JuMP.@constraint(model, A * x .== b)

    vec_dim = Cones.svec_length(q)
    Q = V * diagm(x) * V' # information matrix
    Q_vec = zeros(JuMP.AffExpr, vec_dim)
    Cones.smat_to_svec!(Q_vec, Q, sqrt(T(2)))

    # convex objective
    JuMP.@variable(model, epi)
    JuMP.@objective(model, Min, epi)
    add_homog_spectral(inst.ext, q, vcat(1.0 * epi, Q_vec), model)

    # save for use in tests
    model.ext[:Q_var] = Q

    return model
end

function test_extra(inst::ExperimentDesignJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    (stat == MOI.OPTIMAL) || return

    # check objective
    tol = eps(T)^0.2
    Q_opt = JuMP.value.(model.ext[:Q_var])
    λ = eigvals(Symmetric(Q_opt, :U))
    @test minimum(λ) >= -tol
    obj_result = get_val(pos_only(λ), inst.ext)
    @test JuMP.objective_value(model) ≈ obj_result atol=tol rtol=tol
    return
end


text:
# import Base: show

"""
The structure used for environment requests.
"""
mutable struct EnvironmentRequest
  id::String
  name::String
  description::String
  status::String
end

"""
The structure returned when any environment requests are made.
"""
mutable struct EnvironmentResponse
  sessions::Vector{Dict{String, Any}}
  id::String
  name::String
  description::String
  createdTimestamp::String
  lastUpdatedTimestamp::String
end

function show(io::IO, r::EnvironmentResponse)
    println(io, "GraffSDK Environment:")
    println(io, " - ID: $(r.id)")
    println(io, "   - Name: $(r.name)")
    println(io, "   - Description: $(r.description)")
    println(io, "   - Sessions Associated: $(length(r.sessions)) sessions")
    println(io, "   - Created: $(r.createdTimestamp)")
    println(io, "   - Last Updated: $(r.lastUpdatedTimestamp)")
end


text:
#! /usr/bin/julia

# Rosetta Code, Random number generator (device)

const rdev = "/dev/random"
rstream = try
    open(rdev, "r")
catch
    false
end

if isa(rstream, IOStream)
    b = readbytes(rstream, 4)
    close(rstream)
    i = reinterpret(Int32, b)[1]
    println("A hardware random number is:  ", i)
else
    println("The hardware random number stream, ", rdev, ", was unavailable.")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using CSV, CategoricalArrays, Plots
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.27"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	df = CSV.read(datadir("companyData.csv"))
	companies = levels(df.Type)

	year2012 = df[df.Year .== 2012, :MarketCap]
	year2016 = df[df.Year .== 2016, :MarketCap]
end

# ╔═╡ 9b9baab2-2522-11eb-0bfb-0bdfde4381ea
begin
	p1 = pie(companies, year2012, title="2012 Market Cap \n by company")
	p2 = pie(companies, year2016, title="2016 Market Cap \n by company")
	plot(p1, p2, size=(800, 400))
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.27"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═9b9baab2-2522-11eb-0bfb-0bdfde4381ea
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
@testset "getindex" begin
    @testset "getindex(::Matrix{<:Number},...)" begin
        x = [1.0 2.0 3.0; 10.0 20.0 30.0]
        x̄ = [1.4 2.5 3.7; 10.5 20.1 30.2]
        full_ȳ = [7.4 5.5 2.7; 8.5 11.1 4.2]

        @testset "single element" begin
            rrule_test(getindex, 2.3, (x, x̄), (2, nothing))
            rrule_test(getindex, 2.3, (x, x̄), (2, nothing), (1, nothing))
            rrule_test(getindex, 2.3, (x, x̄), (2, nothing), (2, nothing))

            rrule_test(getindex, 2.3, (x, x̄), (CartesianIndex(2, 3), nothing))
        end

        @testset "slice/index postions" begin
            rrule_test(getindex, [2.3, 3.1], (x, x̄), (2:3, nothing))
            rrule_test(getindex, [2.3, 3.1], (x, x̄), (3:-1:2, nothing))
            rrule_test(getindex, [2.3, 3.1], (x, x̄), ([3,2], nothing))
            rrule_test(getindex, [2.3, 3.1], (x, x̄), ([2,3], nothing))

            rrule_test(getindex, [2.3 3.1; 4.1 5.1], (x, x̄), (1:2, nothing), (2:3, nothing))
            rrule_test(getindex, [2.3 3.1; 4.1 5.1], (x, x̄), (:, nothing), (2:3, nothing))

            rrule_test(getindex, [2.3, 3.1], (x, x̄), (1:2, nothing), (1, nothing))
            rrule_test(getindex, [2.3, 3.1], (x, x̄), (1, nothing), (1:2, nothing))

            rrule_test(getindex, [2.3 3.1; 4.1 5.1], (x, x̄), (1:2, nothing), (2:3, nothing))
            rrule_test(getindex, [2.3 3.1; 4.1 5.1], (x, x̄), (:, nothing), (2:3, nothing))

            rrule_test(getindex, full_ȳ, (x, x̄), (:, nothing), (:, nothing))
            rrule_test(getindex, full_ȳ[:], (x, x̄), (:, nothing))
        end

        @testset "masking" begin
            rrule_test(getindex, full_ȳ, (x, x̄), (trues(size(x)), nothing))
            rrule_test(getindex, full_ȳ[:], (x, x̄), (trues(length(x)), nothing))

            mask = falses(size(x))
            mask[2,3] = true
            mask[1,2] = true
            rrule_test(getindex, [2.3, 3.1], (x, x̄), (mask, nothing))

            rrule_test(
                getindex, full_ȳ[1,:], (x, x̄), ([true, false], nothing), (:, nothing)
            )
        end

        @testset "By position with repeated elements" begin
            rrule_test(getindex, [2.3, 3.1], (x, x̄), ([2, 2], nothing))
            rrule_test(getindex, [2.3, 3.1, 4.1], (x, x̄), ([2, 2, 2], nothing))
            rrule_test(
                getindex, [2.3 3.1; 4.1 5.1], (x, x̄), ([2,2], nothing), ([3,3], nothing)
            )
        end
    end
end


text:
#  MIT License
#
#  Copyright (c) 2020 Michael J Simms. All rights reserved.
#
#  Permission is hereby granted, free of charge, to any person obtaining a copy
#  of this software and associated documentation files (the "Software"), to deal
#  in the Software without restriction, including without limitation the rights
#  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#  copies of the Software, and to permit persons to whom the Software is
#  furnished to do so, subject to the following conditions:
#
#  The above copyright notice and this permission notice shall be included in all
#  copies or substantial portions of the Software.
#
#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#  SOFTWARE.

module Distance

export hamming_distance
export levenshtein_distance
export euclidian_distance_1D
export euclidian_distance_2D

# Hamming Distance
function hamming_distance(str1, str2)
	len1 = length(str1)

	# The strings must be the same length.
	if len1 != length(str2)
		return -1
    end

	distance = 0

    for i in 1:len1
		if str1[i] != str2[i]
			distance = distance + 1
        end
    end

	distance
end

# Levenshtein Distance
function levenshtein_distance(str1, str2)
	m = length(str1)
	n = length(str2)

	# Test for empty strings.
	if m == 0
		return n
    end
	if n == 0
		return m
    end

	# Temp vectors.
	v0 = zeros(n)
	v1 = zeros(n)
    for i in 1:n
		v0[i] = i - 1
    end

    for i in 1:m - 1
		v1[1] = i + 1

        for j in 1:n - 1
            costs = zeros(3) # deletion, insertion, and substitution costs (in that order)
			costs[1] = v0[j + 1] + 1
			costs[2] = v1[j] + 1
			if str1[i] == str2[j]
				costs[3] = v0[j]
			else
				costs[3] = v0[j] + 1
            end

            tempJ = j + 1
			v1[tempJ] = costs[1]
            if costs[2] < v1[tempJ]
                v1[tempJ] = costs[2]
            end
            if costs[3] < v1[tempJ]
                v1[tempJ] = costs[3]
            end
		end

		# Copy v1 to v0.
        tempV = v0
        v0 = v1
        v1 = tempV
    end

	# Save the final distance calculation.
	distance = v0[n]

	distance
end

# 1 dimensional Euclidian distance
function euclidian_distance_1D(pt1, pt2)
	abs(pt1 - pt2)
end

# 2 dimensional Euclidian distance
function euclidian_distance_2D(pt1X, pt1Y, pt2X, pt2Y)
	tempX = pt2X - pt1X
	tempY = pt2Y - pt1Y
	sqrt((tempX * tempX) + (tempY * tempY))
end

end

max_stars_repo_name: msimms/LibMath

max_issues_repo_name: msimms/LibMath

max_forks_repo_name: msimms/LibMath


text:
using Convex, SCS, SchattenNorms

set_default_solver(SCSSolver(verbose=2))

K = 1000
n = 2

A = randn(n, n);
U,_,_ = svd(randn(n,n))
p = 0.0
ρ = (1-p)*U*[1 0; 0 0]*U'+p*eye(2)/2

nzm = round(Int, trace(ρ*[1 0; 0 0]) * K |> real)
nxm = round(Int, trace(ρ*[1 -1; -1 1]/2) * K |> real)

println("(x,z) = ($nxm,$nzm)")

x = Semidefinite(n)

constraints = trace(x) == 1

β = 0.05
problem = maximize( nzm * log(trace(x*[1 0; 0 0])) +
                    (K-nzm) * log(trace(x*[0 0; 0 1])) +
                    nxm * log(trace(x*[1 -1; -1 1])/2) +
                    (K-nxm) * log(trace(x*[1 1; 1 1])/2) +
                    # β * logdet(x),
                    β * log(lambda_min(x)),
                    constraints )

# Solve the problem by calling solve!
solve!(problem)

# Check the status of the problem
println(problem.status) # :Optimal, :Infeasible, :Unbounded etc.

# Get the optimum value
println(problem.optval)

#println(x.value)

println(round(Int, trace(x.value*[1 0; 0 0]) * K))
println(round(Int, trace(x.value*[1 -1; -1 1]/2) * K))

objective(x_) = nzm * log(trace(x_*[1 0; 0 0])) +
                (K-nzm) * log(trace(x_*[0 0; 0 1])) +
                nxm * log(trace(x_*[1 -1; -1 1])/2) +
                (K-nxm) * log(trace(x_*[1 1; 1 1])/2) +
                β * minimum(eigvals(x_))

println(ρ)
println(x.value)
println(trnorm(ρ-x.value)/2)