text:
#!/usr/bin/env julia
include("amp_mgen.jl")
using Mustache
postpath = joinpath("..", "..", "posts")
csvraw = joinpath("data", "BE0101AW.csv")
csvform = joinpath("data", "popdths15.csv")
run(pipeline(`iconv -f ISO_8859-1 -t UTF8`, stdin=csvraw
	, stdout=pipeline(`./scbpopdths.sed`, csvform)
	))
popdths15 = readtable(csvform)
popdthssplit(sex) = @from i in popdths15 begin
	@where i.sex==sex
        @select {i.age,i.pop,i.dths}
	@collect DataFrame
end
agesum(df) = DataFrame(age=df[:age][1:106]
	, pop = [df[:pop][1:105]; sum(df[:pop][106:end])]
	, dths = [df[:dths][1:105]; sum(df[:dths][106:end])]
	)
chfs = Dict(
	:circ => Dict(:noconv => 0.75
		, :conv => [fill(0.75, 80); (0.75*100:100)./100])
	, :noncirc => Dict(:noconv => 0.8
		, :conv => [fill(0.8, 80); (0.8*100:100)./100; fill(1,5)])
	)
popdths_se15 = map(s -> popdthssplit(s)|>agesum, [1;2])
ch_sexes_se15(ca, chf_ca, chf_nonca, pwr) =
	map(s -> chfac(4290, s, ca, 2015, chf_ca.^pwr, chf_nonca.^pwr, popdths_se15[s]), [1;2])

circse15_conv_noconv(p) =
	map(c -> (c, ch_sexes_se15("circ", chfs[:circ][c], chfs[:noncirc][c], p)),
		[:noconv; :conv]) |> Dict

if !(isdefined(:se3))
	se3 = circse15_conv_noconv(3)
end

percround(v) = replace("$(round(v, 2))", ".", ",")

se3vs = Dict(
	:e50_m => se3[:noconv][1][:lt][:e][51]
	, :e50_m_rc_nc => se3[:noconv][1][:lt_rc][:e][51]
	, :e50_m_rc_c => se3[:conv][1][:lt_rc][:e][51]
	, :e50_f => se3[:noconv][2][:lt][:e][51]
	, :e50_f_rc_nc => se3[:noconv][2][:lt_rc][:e][51]
	, :e50_f_rc_c => se3[:conv][2][:lt_rc][:e][51]
	, :circp_m => se3[:noconv][1][:changes][:ptoages2][end]*100
	, :circp_m_rc_nc => se3[:noconv][1][:changes][:ptoages1][end]*100
	, :circp_m_rc_c => se3[:conv][1][:changes][:ptoages1][end]*100
	, :circp_f => se3[:noconv][2][:changes][:ptoages2][end]*100
	, :circp_f_rc_nc => se3[:noconv][2][:changes][:ptoages1][end]*100
	, :circp_f_rc_c => se3[:conv][2][:changes][:ptoages1][end]*100
	)
se3vs_f = Dict(zip(keys(se3vs), map(percround, values(se3vs))))

tpl = readstring("2018-06-17-extra.mustache")
write(joinpath(postpath, "2018-06-17-extra.md"), render(tpl, se3vs_f))

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
"""
    Debug(i::Int32)
Type for debugging and precompile.
"""
type Debug <: IO_Object
    i::Int32
end
write!(::Debug, val, debug::Bool=false) = nothing
read(::Debug, debug::Bool=false) = -1
teardown(::Debug, debug::Bool=false) = nothing

ext: jl


text:
#=
composite:
- Julia version:
- Author: Dan
- Date: 2021-06-04
=#

export
    Composite

"""
`Composite` implements the minimum set of fields that `System`, `Chain`, `Residue`, `Atom` must also
implement. An object of `Composite` is a Node in a hierarchical tree of many types which subtype
`AbstractComposite`.\n
See also [`AbstractComposite`](@ref),[`System`](@ref), [`Chain`](@ref), [`Residue`](@ref), [`Atom`](@ref)
"""
mutable struct Composite <: AbstractComposite
    name_                                       ::String
    number_of_children_                         ::Int64
    parent_                                     ::Union{AbstractComposite, Nothing}
    previous_                                   ::Union{Nothing, AbstractComposite}
    next_                                       ::Union{AbstractComposite, Nothing}
    first_child_                                ::Union{AbstractComposite, Nothing}
    last_child_                                 ::Union{AbstractComposite, Nothing}
    properties_                                 ::UInt64
    contains_selection_                         ::Bool
    number_of_selected_children_                ::Int64
    number_of_children_containing_selection_    ::Int64
    selection_stamp_                            ::Union{TimeStamp,Nothing}
    modification_stamp_                         ::Union{TimeStamp,Nothing}
    trait_                                      ::Union{Nothing,AbstractComposite}

    #default constructor
    Composite() = Composite("",0,nothing,nothing,nothing,nothing,nothing,UInt64(0),false,0,0,nothing,nothing,nothing)
    Composite(xd::Bool) = new(0,nothing,nothing,nothing,nothing,nothing,nothing,0,false
    ,0,0,nothing,nothing,nothing)

    #full constructor
    Composite(
        name_                                   ::String,
        number_of_children                      ::Int64,
        parent                                  ::Composite,
        previous                                ::Composite,
        next                                    ::Composite,
        first_child                             ::Composite,
        last_child                              ::Composite,
        properties                              ::UInt8,
        contains_selection                      ::Bool,
        number_of_selected_children             ::Int64,
        number_of_children_containing_selection ::Int64,
        selection_stamp                         ::TimeStamp,
        modification_stamp                      ::TimeStamp,
        trait                                   ::AbstractComposite

    ) = new(
            name_,
            number_of_children,
            parent,
            previous,
            next,
            first_child,
            last_child,
            properties,
            contains_selection,
            number_of_selected_children,
            number_of_children_containing_selection,
            selection_stamp,
            modification_stamp,
            trait
        )
end
#partial constructor: Refs to other Composites and the obj itself
Composite(
    number_of_children              ::Int64,
    parent                          ::AbstractComposite,
    prev                            ::AbstractComposite,
    next                            ::AbstractComposite,
    first_child                     ::AbstractComposite,
    last_child                      ::AbstractComposite,
    trait                           ::AbstractComposite
) = begin

    Composite(
        number_of_children,
        parent,
        prev,
        next,
        first_child,
        last_child,
        0,
        false,
        false,
        0,
        0,
        nothing,
        nothing,
        trait,
    )
end


text:
"""
# module InputModel

# Examples

```jldoctest
julia>
```
"""
module InputModel

using Parsso.CoreDataStructure.Namelists
using Parsso.CoreDataStructure.Cards

export NamelistsInput,
    NamelistsCardsInput

abstract type AbstractInput end

struct NamelistsInput <: AbstractInput
    program::String
    namelists::Vector{Namelist}
    function NamelistsInput(program, namelists)
        program in QE_PROGRAMS || throw(DomainError("The program $program is not recognized!"))
        new(program, namelists)
    end
end

struct NamelistsCardsInput <: AbstractInput
    program::String
    namelists::Vector{Namelist}
    cards::Vector{Card}
    function NamelistsCardsInput(program, namelists, cards)
        program in QE_PROGRAMS || throw(DomainError("The program $program is not recognized!"))
        new(program, namelists, cards)
    end
end

end


text:
using PkgBenchmark
include("pprinthelper.jl")
result = PkgBenchmark.readresults(joinpath(@__DIR__, "result.json"))
displayresult(result)


text:
export QuantumCircuit, QulacsQuantumCircuit
export ParametricQuantumCircuit, QulacsParametricQuantumCircuit
export add_parametric_multi_Pauli_rotation_gate!, set_parameter!
export add_CNOT_gate!, add_H_gate!, add_X_gate!

export QuantumState, QulacsQuantumState
export set_computational_basis!, create_hf_state
export FermionOperator, jordan_wigner, get_number_preserving_sparse_operator
export get_vector, get_n_qubit

export QubitOperator, OFQubitOperator
export get_term_count, get_n_qubit, terms_dict, is_hermitian
export get_expectation_value, get_transition_amplitude
export inner_product

export get_expectation_value, create_operator_from_openfermion
export hermitian_conjugated
export create_observable

using PyCall

up_index(i) = 2*(i-1)
down_index(i) = 2*(i-1)+1

@enum PauliID pauli_I=0 pauli_X=1 pauli_Y=2 pauli_Z=3
pauli_id_lookup = Dict("I"=>pauli_I, "X"=>pauli_X, "Y"=>pauli_Y, "Z"=>pauli_Z)
pauli_id_str = ["I", "X", "Y", "Z"]

################################################################################
############################# QUANTUM  CIRCUIT #################################
################################################################################
abstract type QuantumCircuit end
abstract type ParametricQuantumCircuit <: QuantumCircuit end

#  QuantumCircuit
struct QulacsQuantumCircuit <: QuantumCircuit
    pyobj::PyObject
end

function QulacsQuantumCircuit(n_qubit::Int)
    QulacsQuantumCircuit(qulacs.QuantumCircuit(n_qubit))
end

function add_X_gate!(circuit::QulacsQuantumCircuit, idx_qubit::Int)
    circuit.pyobj.add_X_gate(idx_qubit)
end

function add_H_gate!(circuit::QulacsQuantumCircuit, idx_qubit::Int)
    circuit.pyobj.add_H_gate(idx_qubit)
end

function add_CNOT_gate!(circuit::QulacsQuantumCircuit, control::Int, target::Int)
    circuit.pyobj.add_CNOT_gate(control, target)
end

#  ParametricQuantumCircuit
struct QulacsParametricQuantumCircuit <: ParametricQuantumCircuit
    pyobj::PyObject
end

function QulacsParametricQuantumCircuit(n_qubit::Int)
    QulacsParametricQuantumCircuit(qulacs.ParametricQuantumCircuit(n_qubit))
end

function add_parametric_multi_Pauli_rotation_gate!(circuit::QulacsParametricQuantumCircuit,
    pauli_indices::Vector{Int}, pauli_ids::Vector{PauliID}, theta_init::Float64=0.0)
    circuit.pyobj.add_parametric_multi_Pauli_rotation_gate(
        pauli_indices, Int.(pauli_ids), theta_init)
end

function set_parameter!(circuit::QulacsParametricQuantumCircuit, index, theta)
    circuit.pyobj.set_parameter(index, theta)
end

################################################################################
############################# QUANTUM STATE  ###################################
################################################################################
abstract type QuantumState end

# Wrap qulacs.QuantumState
struct QulacsQuantumState <: QuantumState
    pyobj
end

function QulacsQuantumState(n_qubit::Int, int_state=0b0)
    state = QulacsQuantumState(qulacs.QuantumState(n_qubit))
    set_computational_basis!(state, int_state)
    state
end

function Base.:copy(state::QulacsQuantumState)
    QulacsQuantumState(state.pyobj.copy())
end

function set_computational_basis!(state::QulacsQuantumState, int_state)
    state.pyobj.set_computational_basis(int_state)
end

function create_hf_state(n_qubit, n_electron)
    int_state = parse(Int, repeat("0", n_qubit-n_electron) * repeat("1", n_electron), base=2)
    state = QulacsQuantumState(n_qubit)
    set_computational_basis!(state, int_state)
    state
end

function get_n_qubit(state::QulacsQuantumState)
    state.pyobj.get_qubit_count()
end

function get_vector(state::QuantumState)
    state.pyobj.get_vector()
end

################################################################################
################################# OPERATOR #####################################
################################################################################
abstract type SecondQuantOperator end

function Base.:(==)(op1::SecondQuantOperator, op2::SecondQuantOperator)
    op1.pyobj == op2.pyobj
end

function Base.:*(op1::SecondQuantOperator, op2::SecondQuantOperator)
    typeof(op1)(op1.pyobj * op2.pyobj)
end

function Base.:+(op1::SecondQuantOperator, op2::SecondQuantOperator)
    typeof(op1)(op1.pyobj + op2.pyobj)
end

function Base.:-(op1::SecondQuantOperator, op2::SecondQuantOperator)
    typeof(op1)(op1.pyobj - op2.pyobj)
end

function Base.:/(op::SecondQuantOperator, x::Number)
    typeof(op)(op.pyobj/x)
end

################################################################################
############################# FERMION OPERATOR #################################
################################################################################
# Wrap openfermion.ops.operators.FermionOperator
struct FermionOperator <: SecondQuantOperator
    pyobj::PyObject
end

function FermionOperator(ops::Vector{Tuple{Int,Int}}, coeff::Number=1.0)
    FermionOperator(ofermion.ops.operators.FermionOperator(Tuple(ops), coeff))
end

function FermionOperator(op_str::String="", coeff::Number=1.0)
    FermionOperator(ofermion.ops.operators.FermionOperator(op_str, coeff))
end

function get_number_preserving_sparse_operator(ham::FermionOperator, n_qubit::Int, n_electron::Int)::PyObject
    ofermion.linalg.get_number_preserving_sparse_operator(ham.pyobj, n_qubit, n_electron)
end

function jordan_wigner(op::FermionOperator)
    OFQubitOperator(ofermion.transforms.jordan_wigner(op.pyobj))
end

################################################################################
############################# QUBIT OPERATOR ###################################
################################################################################
abstract type QubitOperator <: SecondQuantOperator end

# Wrap openfermion.ops.operators.qubit_operator.QubitOperator
struct OFQubitOperator <: QubitOperator
    pyobj
end

function OFQubitOperator(str::String, coeff::Number=1.0)
    return OFQubitOperator(ofermion.ops.operators.qubit_operator.QubitOperator(str, coeff))
end

function hermitian_conjugated(op::OFQubitOperator)
    OFQubitOperator(ofermion.utils.hermitian_conjugated(op.pyobj))
end

function get_n_qubit(op::OFQubitOperator)
    count_qubit_in_qubit_operator(op.pyobj)
end

function get_term_count(op::OFQubitOperator)
    length(op.pyobj.terms)
end

function terms_dict(op::OFQubitOperator)::Dict{Any,Any}
    op.pyobj.terms
end

function is_hermitian(op::OFQubitOperator)
    ofermion.utils.operator_utils.is_hermitian(op.pyobj)
end

function convert_to_qulacs_op(op::OFQubitOperator, n_qubit::Int)
    convert_openfermion_op(n_qubit, op.pyobj)
end

function get_expectation_value(op::OFQubitOperator, state::QulacsQuantumState)
    if !is_hermitian(op)
        error("op must be hermite for get_expectation_value!")
    end
    convert_to_qulacs_op(op, get_n_qubit(state)).get_expectation_value(state.pyobj)
end

function get_transition_amplitude(op::OFQubitOperator, state_bra::QulacsQuantumState, state_ket::QulacsQuantumState)
    convert_to_qulacs_op(op, get_n_qubit(state_bra)).get_transition_amplitude(state_bra.pyobj, state_ket.pyobj)
end

function inner_product(state_bra::QulacsQuantumState, state_ket::QulacsQuantumState)
    qulacs.state.inner_product(state_bra.pyobj, state_ket.pyobj)
end

"""
Parse a tuple representing a Pauli string
  When x is ((0, "X"), (5, "Y")), returns [0, 5], [PauliID.X, PauliID.Y]
"""
function parse_pauli_str(x)::Tuple{Vector{Int64}, Vector{PauliID}}
    collect(e[1] for e in x), collect(pauli_id_lookup[e[2]] for e in x)
end

"""
Update a state using a circuit
"""
function update_quantum_state!(circuit::QulacsParametricQuantumCircuit, state::QulacsQuantumState)
    circuit.pyobj.update_quantum_state(state.pyobj)
end

function update_quantum_state!(circuit::QulacsQuantumCircuit, state::QulacsQuantumState)
    circuit.pyobj.update_quantum_state(state.pyobj)
end


text:
export proper_divisors, is_abundant

is_abundant(n) = begin
    sum(proper_divisors(n)) > n
end


text:
# create bits types for easy loading of bytes of lengths up to 15
primitive type Bits24 24 end
primitive type Bits40 40 end
primitive type Bits48 48 end
primitive type Bits56 56 end
primitive type Bits72 72 end
primitive type Bits80 80 end
primitive type Bits88 88 end
primitive type Bits96 96 end
primitive type Bits104 104 end
primitive type Bits112 112 end
primitive type Bits120 120 end

# loads `remaining_bytes_to_load` bytes from `ptrs` which is a C-style pointer to a string
# these functions assumes that remaining_bytes_to_load > 0
function load_bits_with_padding(::Type{UInt128}, ptrs::Ptr{UInt8}, remaining_bytes_to_load)::UInt128
    nbits_to_shift_away = 8(sizeof(UInt128) - remaining_bytes_to_load)

    # the below checks if the string is less than 16 bytes away from the page
    # boundary assuming the page size is 4kb
    # see https://discourse.julialang.org/t/is-there-a-way-to-check-how-far-away-a-pointer-is-from-a-page-boundary/8147/11?u=xiaodai
    if (UInt(ptrs) & 0xfff) > 0xff0
        if  remaining_bytes_to_load == 15
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits120}(ptrs))))
        elseif  remaining_bytes_to_load == 14
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits112}(ptrs))))
        elseif  remaining_bytes_to_load == 13
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits104}(ptrs))))
        elseif  remaining_bytes_to_load == 12
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits96}(ptrs))))
        elseif  remaining_bytes_to_load == 11
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits88}(ptrs))))
        elseif  remaining_bytes_to_load == 10
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits80}(ptrs))))
        elseif  remaining_bytes_to_load == 9
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits72}(ptrs))))
        elseif  remaining_bytes_to_load == 8
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt64}(ptrs))))
        elseif  remaining_bytes_to_load == 7
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits56}(ptrs))))
        elseif  remaining_bytes_to_load == 6
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits48}(ptrs))))
        elseif  remaining_bytes_to_load == 5
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits40}(ptrs))))
        elseif  remaining_bytes_to_load == 4
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt32}(ptrs))))
        elseif  remaining_bytes_to_load == 3
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits24}(ptrs))))
        elseif  remaining_bytes_to_load == 2
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt16}(ptrs))))
        else
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt8}(ptrs))))
        end
    else
        return ntoh(unsafe_load(Ptr{UInt128}(ptrs))) >> nbits_to_shift_away << nbits_to_shift_away
    end
end

function load_bits_with_padding(::Type{UInt64}, ptrs::Ptr{UInt8}, remaining_bytes_to_load)::UInt64
    nbits_to_shift_away = 8(sizeof(UInt64) - remaining_bytes_to_load)

    # the below checks if the string is less than 8 bytes away from the page
    # boundary assuming the page size is 4kb
    # see https://discourse.julialang.org/t/is-there-a-way-to-check-how-far-away-a-pointer-is-from-a-page-boundary/8147/11?u=xiaodai
    if (UInt(ptrs) & 0xfff) > 0xff8
        if  remaining_bytes_to_load == 7
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits56}(ptrs))))
        elseif  remaining_bytes_to_load == 6
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits48}(ptrs))))
        elseif  remaining_bytes_to_load == 5
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits40}(ptrs))))
        elseif  remaining_bytes_to_load == 4
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{UInt32}(ptrs))))
        elseif  remaining_bytes_to_load == 3
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits24}(ptrs))))
        elseif  remaining_bytes_to_load == 2
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{UInt16}(ptrs))))
        else
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{UInt8}(ptrs))))
        end
    else
        return ntoh(unsafe_load(Ptr{UInt64}(ptrs))) >> nbits_to_shift_away << nbits_to_shift_away
    end
end

function load_bits_with_padding(::Type{UInt32}, ptrs::Ptr{UInt8}, remaining_bytes_to_load)::UInt32
    nbits_to_shift_away = 8(sizeof(UInt32) - remaining_bytes_to_load)

    # the below checks if the string is less than 4 bytes away from the page
    # boundary assuming the page size is 4kb
    # see https://discourse.julialang.org/t/is-there-a-way-to-check-how-far-away-a-pointer-is-from-a-page-boundary/8147/11?u=xiaodai
    if (UInt(ptrs) & 0xfff) > 0xffc
        if  remaining_bytes_to_load == 3
            return ntoh(Base.zext_int(UInt32, unsafe_load(Ptr{Bits24}(ptrs))))
        elseif  remaining_bytes_to_load == 2
            return ntoh(Base.zext_int(UInt32, unsafe_load(Ptr{UInt16}(ptrs))))
        else
            return ntoh(Base.zext_int(UInt32, unsafe_load(Ptr{UInt8}(ptrs))))
        end
    else
        return ntoh(unsafe_load(Ptr{UInt32}(ptrs))) >> nbits_to_shift_away << nbits_to_shift_away
    end
end

"""
    load_bits([type,] s, skipbytes)

Load the underlying bits of a string `s` into a `type` of the user's choosing.
The default is `UInt`, so on a 64 bit machine it loads 64 bits (8 bytes) at a time.
If the `String` is shorter than 8 bytes then it's padded with 0.

- `type`:       any bits type that has `>>`, `<<`, and `&` operations defined
- `s`:          a `String`
- `skipbytes`:  how many bytes to skip e.g. load_bits("abc", 1) will load "bc" as bits
"""
# Some part of the return result should be padded with 0s.
# To prevent any possibility of segfault we load the bits using
# successively smaller types
# it is assumed that the type you are trying to load into needs padding
# i.e. `remaining_bytes_to_load > 0`
function load_bits(::Type{T}, s::String, skipbytes = 0)::T where T
    n = sizeof(s)
    load_bits(T, pointer(s), n, skipbytes)
    # if n < skipbytes
    #     res = zero(T)
    # elseif n - skipbytes >= sizeof(T)
    #     res = ntoh(unsafe_load(Ptr{T}(pointer(s, skipbytes+1))))
    # else
    #     ptrs  = pointer(s) + skipbytes
    #     remaining_bytes_to_load = n - skipbytes
    #     res = load_bits_with_padding(T, ptrs, remaining_bytes_to_load)
    # end
    # return res
end

function load_bits(::Type{T}, s::Ptr{UInt8}, n, skipbytes = 0)::T where T
    if n < skipbytes
        res = zero(T)
    elseif n - skipbytes >= sizeof(T)
        res = ntoh(unsafe_load(Ptr{T}(s + skipbytes)))
    else
        ptrs  = s + skipbytes
        remaining_bytes_to_load = n - skipbytes
        res = load_bits_with_padding(T, ptrs, remaining_bytes_to_load)
    end
    return res
end


text:
using ConstrainedDynamics

# Parameters
joint_axis = [1.0;0.0;0.0]

length1 = 0.5
width, depth = 0.5, 0.5
box1 = Box(width, depth, length1, 1., color = RGBA(1., 1., 0.))

# Links
origin = Origin{Float64}()
link1 = Body(box1)

# Constraints
joint1 = InequalityConstraint(Impact(link1, [0;-0.1;1.0]), Impact(link1, [0;0.1;1.0]))
joint2 = InequalityConstraint(Impact(link1, [0.1;0;1.0]))
joint3 = InequalityConstraint(Impact(link1, [-0.1;0;1.0]))

links = [link1]
ineqs = [joint1;joint2;joint3]
shapes = [box1]

mech = Mechanism(origin, links, ineqs, shapes = shapes)
setVelocity!(link1,v = [1;0.5;5])


text:
# Write a function that combines two lists by alternatingly taking elements,
# e.g. [a,b,c], [1,2,3] → [a,1,b,2,c,3].
using BenchmarkTools

a = ["a","b","c","d","e","F","g","h"]
b = [1,2,3,4,5,6,7,8]

function inter_concat(x::Vector{String},y::Vector{Int64})
    c = Vector{Any}(undef, 2*length(a))
    i = 0
    for (x,y) in zip(a,b)
        c[i += 1] = x
        c[i += 1] = y
    end
    return c
end

@benchmark inter_concat(a, b)

@benchmark collect(Iterators.flatten(zip(a, b)))

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
using Dates
using GaussQuadrature
using JLD2
using LinearAlgebra
using Logging
using MPI
using Printf
using StaticArrays
using Test
using Unitful

#add https://github.com/CliMA/ClimateMachine.jl#tb/refactoring_ans_sphere
using ClimateMachine
using ClimateMachine.Atmos: NoReferenceState
using ClimateMachine.MPIStateArrays
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Geometry
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.BalanceLaws
using ClimateMachine.ODESolvers
using ClimateMachine.SystemSolvers
using ClimateMachine.Orientations
using ClimateMachine.VTK

# to be removed
using ClimateMachine.GenericCallbacks: EveryXWallTimeSeconds, EveryXSimulationSteps
using ClimateMachine.Thermodynamics: soundspeed_air
using ClimateMachine.TemperatureProfiles
using ClimateMachine.VariableTemplates: flattenednames

# to be removed
using CLIMAParameters: AbstractEarthParameterSet
using CLIMAParameters.Planet: grav, R_d, cp_d, cv_d, planet_radius, MSLP, Omega, T_0

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import CLIMAParameters
import CLIMAParameters.Planet.T_0
CLIMAParameters.Planet.T_0(::EarthParameterSet) = 0.0

const total_energy = true
const fluctuation_gravity = true

# utils
include("../utilities/operations.jl")
include("../utilities/sphere_utils.jl")
# grid
include("../grid/domains.jl")
include("../grid/grids.jl")
# numerics
include("../models.jl")
include("../numerics/filters.jl")
# physics
include("../physics/physics.jl")
# include("../grid/orientations.jl")
include("../physics/advection.jl")
include("../physics/thermodynamics.jl")
include("../physics/coriolis.jl")
include("../physics/gravity.jl")
include("../physics/pressure_force.jl")
include("../physics/diffusion.jl")
# boundary conditions
include("../boundary_conditions/boundary_conditions.jl")
include("../boundary_conditions/bc_first_order.jl")
include("../boundary_conditions/bc_second_order.jl")
include("../boundary_conditions/bc_gradient.jl")
# parameters and initial conditions
include("../parameters_initialconditions.jl")
# coupler
include("../coupler/surface_fluxes.jl")
include("../coupler/put_get_calls.jl")
# interface
include("../balance_law_interface_atmos.jl")
include("../balance_law_interface_slabland.jl")
include("../numerics/numerical_volume_fluxes.jl")
include("../numerics/numerical_interface_fluxes.jl")
include("../simulations.jl")
include("../utilities/callbacks.jl")

ClimateMachine.init()

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using LevyArea
using Test, BenchmarkTools
import Random
import LinearAlgebra: diag

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1
Random.seed!(638278)

@testset "Iterated Integrals" begin
    @testset "m=1, h=$h" for h in rand(5)
        W = √h * randn()
        Ints = iterated_integrals(W, h, h^(3/2))
        @test Ints == 0.5W^2 - 0.5h
    end

    #trigger compilation
    for alg ∈ LevyArea.ITER_INT_ALGS
        iterated_integrals(randn(2),0.1,0.1, alg=alg)
    end

    h = 0.0001
    ϵ = h^(3/2)
    println("h = $h, ϵ = $ϵ")
    @testset "Alg $alg: m=$m" for alg ∈ LevyArea.ITER_INT_ALGS, m in [2;10;50;100;500]
        W = √h * randn(m)
        print("Alg $alg: m=$m n=$(terms_needed(m,h,ϵ,alg,MaxL2())) cost=$(LevyArea.effective_cost(m,h,ϵ,alg,MaxL2()))")
        @btime iterated_integrals($W, $h, $ϵ, alg=$alg) # stepsize h
        Ints = iterated_integrals(W, h, ϵ, alg=alg) # stepsize h
        @test diag(Ints) ≈ 0.5*W.^2 .- 0.5h
    end
end

max_stars_repo_name: mrbauff/MasterThesis

max_issues_repo_name: mrbauff/MasterThesis


text:
"""
Returns a very (very) rough estimate of the time per SCF step (in seconds),
assuming FFTs are the limiting operation.
"""
function estimate_time_of_scf_step(basis::PlaneWaveBasis)
    # TODO pretty print this instead of just seconds.
    # There's nothing in julia stdlib for this apparently...
    # Super rough figure from various tests on cluster, laptops etc.
    # on a 128^3 FFT grid.
    time_per_FFT_per_grid_point = 30 #= ms =# / 1000 / 128^3

    (time_per_FFT_per_grid_point
     * prod(basis.fft_size)
     * length(basis.kpoints)
     * div(basis.model.n_electrons, filled_occupation(basis.model), RoundUp)
     * 8  # mean number of FFT steps per state per k-point per iteration
     )
end


text:
# Julia1_3.4.jl Methods
#
# restricting cylinder function parameters to Float64
# "multiple dispatch" means that all parameters are used
#   to decide which method to use.
# So for efficiency we restrict the call parameters using
#  type assertion operators so that the operations
#  in the expression π*h*r^2 won't have to decide.
cylinder(r::Float64,h::Float64) = π*h*r^2  # for cylinder volume
println(cylinder(5.0,20.0))
println(cylinder(5,20.0))
println(typeof(5))
#
# Generalizing but still efficient for Float64 since it
# will use method 1.
# The most specific specifications take dispatch priority
cylinder(r::Real,h::Real) = π*h*r^2   # for cylinder volume
println(cylinder(5,20.0))

methods(cylinder) # we defined 2 methods for cylinder
#
# What if we just defined cylinder(r,h) = π*h*r^2
#  r and h would be of type 'any' which includes
#  Number, String, Char or any user defined type, often
#    undesirable
cuboid(l,h,d)=l*h*d   # will even take strings without complaining
cuboid("This"," is ","wrong!")

# perhaps want functions to check if types match
sametype(x::T, y::T) where {T}=true
sametype(x,y)=false
println(sametype(3,3.0))
println(sametype(3,2))
println(sametype(1+2im,3.1+4im))
println(sametype(1.0+2.0im,3.1+4.0im))

# Making sure myappend appends data of the same type
myappend(v::Vector{T}, x::T) where {T} = [v..., x]
y=myappend([1,2,3],4)

w=myappend(["stop","forword","slow"],"reverse")

w[3]


text:
include("start.jl")
include("move_user_move.jl")

#=client.jl - Auth: Kino Roy for Shogi, an assignment for CMPT276:
connects to a server to play a game =#
module networking
using start
using move_user_move

global connection
function getAt(arr, index,backup)
  local res
  try
    res = arr[index]
  catch err
    res = backup

  return res
end
end

function clientStart(idatabase,gameType,cheating,timeLimit,limitAdd,sente_time,gote_time)
  global connection = connect(8080)
  global database = idatabase
  startShogi(database,gameType,cheating,timeLimit,limitAdd,sente_time,gote_time)
  play(gameType,cheating,timeLimit,limitAdd)
end

function clientClose()
  global connection
  close(connection)
end

function play(gameType,cheating,timeLimit,limitAdd)
  global connection
  wcheating = cheating == "F" ? 0 : 1
  write(connection,"""0:$gameType:$wcheating:$timeLimit:$limitAdd\n""")

  while true
    message = chomp(readline(connection))
    sL = split(message,r":",keep=false)
    wincode = getAt(sL,1,-1)
    println(message)
   if wincode == "0"
      #I"m player 1
        global authString = getAt(sL,2,-1)[1]
    elseif wincode == "1"
      #I"m player 2
        global authString = getAt(sL,2,-1)[1]
    elseif wincode == "2"
      #Server quits
      println("gameover:server quit")
    elseif wincode == "3"
      #There"s a draw
      println("gameover:draw")
    elseif wincode == "8"
      #Not your turn
      println("wait")
    elseif wincode == "9"
      #This is your opponents move, it's now your turn
      #Prompt user to make move (GUI)
      #!!!!!! SAMPLE DATA!!!!!!!!!!!!!!!:
      database="game1.db"
      sourcex=1
      soucey=1
      targetx=1
      targety=1
      promote = false

      move_user_move.moveUserMove(database,sourcex,sourcey,targetx,targety,promote::Bool)
      write(connection,"2:$sourcex:$sourcey:$targetx:$targety")
    elseif wincode == "10"
      #bad payload
    elseif wincode == "e"
  end
end
end

end

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


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

"""
    GitRemote(repo::GitRepo, rmt_name::AbstractString, rmt_url::AbstractString) -> GitRemote

Look up a remote git repository using its name and URL. Uses the default fetch refspec.

# Example

```julia
repo = LibGit2.init(repo_path)
remote = LibGit2.GitRemote(repo, "upstream", repo_url)
```
"""
function GitRemote(repo::GitRepo, rmt_name::AbstractString, rmt_url::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_create, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring, Cstring),
                rmt_ptr_ptr, repo.ptr, rmt_name, rmt_url)
    return GitRemote(repo, rmt_ptr_ptr[])
end

"""
    GitRemote(repo::GitRepo, rmt_name::AbstractString, rmt_url::AbstractString, fetch_spec::AbstractString) -> GitRemote

Look up a remote git repository using the repository's name and URL,
as well as specifications for how to fetch from the remote
(e.g. which remote branch to fetch from).

# Example

```julia
repo = LibGit2.init(repo_path)
refspec = "+refs/heads/mybranch:refs/remotes/origin/mybranch"
remote = LibGit2.GitRemote(repo, "upstream", repo_url, refspec)
```
"""
function GitRemote(repo::GitRepo, rmt_name::AbstractString, rmt_url::AbstractString, fetch_spec::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_create_with_fetchspec, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring, Cstring, Cstring),
                rmt_ptr_ptr, repo.ptr, rmt_name, rmt_url, fetch_spec)
    return GitRemote(repo, rmt_ptr_ptr[])
end

"""
    GitRemoteAnon(repo::GitRepo, url::AbstractString) -> GitRemote

Look up a remote git repository using only its URL, not its name.

# Example

```julia
repo = LibGit2.init(repo_path)
remote = LibGit2.GitRemoteAnon(repo, repo_url)
```
"""
function GitRemoteAnon(repo::GitRepo, url::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_create_anonymous, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring),
                rmt_ptr_ptr, repo.ptr, url)
    return GitRemote(repo, rmt_ptr_ptr[])
end

function get(::Type{GitRemote}, repo::GitRepo, rmt_name::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_lookup, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring),
                rmt_ptr_ptr, repo.ptr, rmt_name)
    return GitRemote(repo, rmt_ptr_ptr[])
end

"""
    url(rmt::GitRemote)

Get the URL of a remote git repository.

# Example

```julia-repl
julia> repo_url = "https://github.com/JuliaLang/Example.jl";

julia> repo = LibGit2.clone(cache_repo, "test_directory");

julia> remote = LibGit2.GitRemote(repo, "origin", repo_url);

julia> url(remote)
"https://github.com/JuliaLang/Example.jl"
```
"""
function url(rmt::GitRemote)
    url_ptr = ccall((:git_remote_url, :libgit2), Cstring, (Ptr{Void},), rmt.ptr)
    url_ptr == C_NULL && return ""
    return unsafe_string(url_ptr)
end

"""
    name(rmt::GitRemote)

Get the name of a remote repository, for instance `"origin"`.
If the remote is anonymous (see [`GitRemoteAnon`](@ref))
the name will be an empty string `""`.

# Example

```julia-repl
julia> repo_url = "https://github.com/JuliaLang/Example.jl";

julia> repo = LibGit2.clone(cache_repo, "test_directory");

julia> remote = LibGit2.GitRemote(repo, "origin", repo_url);

julia> name(remote)
"origin"
```
"""
function name(rmt::GitRemote)
    name_ptr = ccall((:git_remote_name, :libgit2), Cstring, (Ptr{Void},), rmt.ptr)
    name_ptr == C_NULL && return ""
    return unsafe_string(name_ptr)
end

"""
    fetch_refspecs(rmt::GitRemote) -> Vector{String}

Get the *fetch* refspecs for the specified `rmt`. These refspecs contain
information about which branch(es) to fetch from.
"""
function fetch_refspecs(rmt::GitRemote)
    sa_ref = Ref(StrArrayStruct())
    @check ccall((:git_remote_get_fetch_refspecs, :libgit2), Cint,
                 (Ptr{StrArrayStruct}, Ptr{Void}), sa_ref, rmt.ptr)
    res = convert(Vector{String}, sa_ref[])
    free(sa_ref)
    res
end

"""
    push_refspecs(rmt::GitRemote) -> Vector{String}

Get the *push* refspecs for the specified `rmt`. These refspecs contain
information about which branch(es) to push to.
"""
function push_refspecs(rmt::GitRemote)
    sa_ref = Ref(StrArrayStruct())
    @check ccall((:git_remote_get_push_refspecs, :libgit2), Cint,
                 (Ptr{StrArrayStruct}, Ptr{Void}), sa_ref, rmt.ptr)
    res = convert(Vector{String}, sa_ref[])
    free(sa_ref)
    res
end

"""
    add_fetch!(repo::GitRepo, rmt::GitRemote, fetch_spec::String)

Add a *fetch* refspec for the specified `rmt`. This refspec will contain
information about which branch(es) to fetch from.

# Example
```julia-repl
julia> LibGit2.add_fetch!(repo, remote, "upstream");

julia> LibGit2.fetch_refspecs(remote)
String["+refs/heads/*:refs/remotes/upstream/*"]
```
"""
function add_fetch!(repo::GitRepo, rmt::GitRemote, fetch_spec::String)
    @check ccall((:git_remote_add_fetch, :libgit2), Cint,
                 (Ptr{Void}, Cstring, Cstring), repo.ptr,
                 name(rmt), fetch_spec)
end

"""
    add_push!(repo::GitRepo, rmt::GitRemote, push_spec::String)

Add a *push* refspec for the specified `rmt`. This refspec will contain
information about which branch(es) to push to.

# Example
```julia-repl
julia> LibGit2.add_push!(repo, remote, "refs/heads/master");

julia> remote = LibGit2.get(LibGit2.GitRemote, repo, branch);

julia> LibGit2.push_refspecs(remote)
String["refs/heads/master"]
```

!!! note
    You may need to [`close`](@ref) and reopen the `GitRemote`
    in question after updating its push refspecs in order for
    the change to take effect and for calls to [`push`](@ref)
    to work.
"""
function add_push!(repo::GitRepo, rmt::GitRemote, push_spec::String)
    @check ccall((:git_remote_add_push, :libgit2), Cint,
                 (Ptr{Void}, Cstring, Cstring), repo.ptr,
                 name(rmt), push_spec)
end

"""
    fetch(rmt::GitRemote, refspecs; options::FetchOptions=FetchOptions(), msg="")

Fetch from the specified `rmt` remote git repository, using `refspecs` to
determine which remote branch(es) to fetch.
The keyword arguments are:
  * `options`: determines the options for the fetch, e.g. whether to prune afterwards.
  * `msg`: a message to insert into the reflogs.
"""
function fetch(rmt::GitRemote, refspecs::Vector{<:AbstractString};
               options::FetchOptions = FetchOptions(),
               msg::AbstractString="")
    msg = "libgit2.fetch: $msg"
    @check ccall((:git_remote_fetch, :libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{FetchOptions}, Cstring),
                 rmt.ptr, isempty(refspecs) ? C_NULL : refspecs, Ref(options), msg)
end

"""
    push(rmt::GitRemote, refspecs; force::Bool=false, options::PushOptions=PushOptions())

Push to the specified `rmt` remote git repository, using `refspecs` to
determine which remote branch(es) to push to.
The keyword arguments are:
  * `force`: if `true`, a force-push will occur, disregarding conflicts.
  * `options`: determines the options for the push, e.g. which proxy headers to use.

!!! note
    You can add information about the push refspecs in two other ways: by setting
    an option in the repository's `GitConfig` (with `push.default` as the key) or
    by calling [`add_push!`](@ref). Otherwise you will need to explicitly specify
    a push refspec in the call to `push` for it to have any effect, like so:
    `LibGit2.push(repo, refspecs=["refs/heads/master"])`.
"""
function push(rmt::GitRemote, refspecs::Vector{<:AbstractString};
              force::Bool = false, options::PushOptions = PushOptions())
    @check ccall((:git_remote_push, :libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{PushOptions}),
                 rmt.ptr, isempty(refspecs) ? C_NULL : refspecs, Ref(options))
end

Base.show(io::IO, rmt::GitRemote) = print(io, "GitRemote:\nRemote name: ", name(rmt), " url: ", url(rmt))

max_stars_repo_name: fcard/julia

max_issues_repo_name: fcard/julia

max_forks_repo_name: fcard/julia


text:
# # Using training an MLJFlux regression model on the Boston house
# # price dataset

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

# **Julia version** is assumed to be 1.6.*

using MLJ
using MLJFlux
using Flux

MLJ.color_off()

using Plots

# This tutorial uses some MLJ's `IteratedModel` wrapper wrappers, to
# transform the MLJFLux `NeuralNetworkRegressor` into a model that
# **automatically selects the number of epochs** required to optimize
# an out-of-sample loss.

# We also show how to include the model in a **pipeline** to carry out
# standardization of the features and target.

# ## Loading data

data = OpenML.load(531); # Loads from https://www.openml.org/d/531

# The target `y` is `:MEDV` and everything else except `:CHAS` goes
# into the features `X`:

y, X = unpack(data, ==(:MEDV), !=(:CHAS));

# (The Charles River dummy variable, `:CHAS`, is not deemed to be
# relevant.)

# Inspect the scientific types:

scitype(y)

#-

schema(X)

# We'll regard `:RAD` (index of accessibility to radial highways) as
# `Continuous` as Flux models don't handle ordered factors:

X = coerce(X, :RAD => Continuous);

# Let's split off a test set for final testing:

X, Xtest = partition(X, 0.7);
y, ytest = partition(y, 0.7);

# ## Defining a builder

# In the macro call below, `n_in` is expected to represent the number
# of inputs features and `rng` a RNG (builders are generic, ie can be
# applied to data with any number of input features):

builder = MLJFlux.@builder begin
    init=Flux.glorot_uniform(rng)
    Chain(Dense(n_in, 64, relu, init=init),
          Dense(64, 32, relu, init=init),
          Dense(32, 1, init=init))
end

# ## Defining a MLJFlux model:

NeuralNetworkRegressor = @load NeuralNetworkRegressor
    model = NeuralNetworkRegressor(builder=builder,
                                   rng=123,
                                   epochs=20)

# ## Standardization

# The following wraps our regressor in feature and target standardizations:

pipe = @pipeline Standardizer model target=Standardizer

# Notice that our original neural network model is now a
# hyper-parameter of the composite `pipe`, with the automatically
# generated name, `:neural_network_regressor`.

# ## Choosing an initial learning rate

# Let's see how the training losses look for the default optimiser. For
# MLJFlux models, `fit!` will print these losses if we bump the
# verbosity level (default is always 1):

mach = machine(pipe, X, y)
fit!(mach, verbosity=2)

# They are also extractable from the training report (which includes
# the pre-train loss):

report(mach).neural_network_regressor.training_losses

# Next, let's visually compare a few learning rates:

plt = plot()
rates = [5e-4, 0.001, 0.005, 0.01]

# By default, changing only the optimiser will not trigger a
# cold-restart when we `fit!` (to allow for adaptive learning rate
# control). So we call `fit!` with the `force=true` option. We'll skip
# the first few losses to get a better vertical scale in our plot.

foreach(rates) do η
    pipe.neural_network_regressor.optimiser.eta = η
    fit!(mach, force=true)
    losses = report(mach).neural_network_regressor.training_losses[3:end]
    plot!(:length(losses), losses, label=η)
end
plt #!md

#-

savefig("learning_rate.png")

# ![](learing_rate.png) #md

# We'll use a relatively conservative rate here:

pipe.neural_network_regressor.optimiser.eta = 0.001

# ## Wrapping in iteration control

# We want a model that trains until an out-of-sample loss satisfies
# the `NumberSinceBest(4)` stopping criterion. We'll add some fallback
# stopping criterion `InvalidValue` and `TimeLimit(1/60)`, and
# controls to print traces of the losses.

# For intializing or clearing the traces:

clear() = begin
    global losses = []
    global training_losses = []
    global epochs = []
    return nothing
end

# And to update the traces:

update_loss(loss) = push!(losses, loss)
update_training_loss(report) =
    push!(training_losses, report.neural_network_regressor.training_losses[end])
update_epochs(epoch) = push!(epochs, epoch)

# The controls to apply:

controls=[Step(1),
          NumberSinceBest(4),
          InvalidValue(),
          TimeLimit(1/60),
          WithLossDo(update_loss),
          WithReportDo(update_training_loss),
          WithIterationsDo(update_epochs)]

# Next we create a "self-iterating" version of the pipeline. Note
# that the iteration parameter is a nested hyperparameter:

iterated_pipe =
    IteratedModel(model=pipe,
                  controls=controls,
                  resampling=Holdout(fraction_train=0.8),
                  iteration_parameter=:(neural_network_regressor.epochs),
                  measure = l2)

# Training the wrapped model on all the train/validation data:

clear()
mach = machine(iterated_pipe, X, y)
fit!(mach)

# And plotting the traces:

plot(epochs, losses,
     xlab = "epoch",
     ylab = "mean sum of sqaures error",
     label="out-of-sample",
     legend = :topleft);
scatter!(twinx(), epochs, training_losses, label="training", color=:red) #!md

#-

savefig("loss.png")

# ![](losses.png) #md

# **How `IteratedModel` works.** Training an `IteratedModel` means
# holding out some data (80% in this case) so an out-of-sample loss
# can be tracked and used in the specified stopping criterion,
# NumberSinceBest(4)`. However, once the stop is triggered, the model
# wrapped by `IteratedModel` (our pipeline model) is retrained on all
# data for the same number of iterations. Calling `predict(mach,
# Xnew)` on new data uses the updated learned parameters.

# In other words, `iterated_model` is a "self-iterating" version of
# the original model, where `epochs` has been transformed from
# hyper-parameter to *learned* parameter.

# ## An evaluation the self-iterating model

# Here's an estimate of performance of our "self-iterating"
# model:

e = evaluate!(mach,
              resampling=CV(nfolds=5),
              measures=[l2, l1])

using Measurements
err = e.measurement[1] ± std(e.per_fold[1])/sqrt(4)
@show err

# which we can see has substantial uncertainty.

# ## Comparing the model with other models

# Although we cannot assign them statistical significance, here are
# comparisons of our fully automated model with a couple other models
# (trained using default hyperparameters):

function performance(model)
    mach = machine(model, X, y) |> fit!
    yhat = predict(mach, Xtest)
    l2(yhat, ytest) |> mean
end
performance(iterated_pipe)

three_models = [(@load EvoTreeRegressor)(), # tree boosting model
                (@load LinearRegressor pkg=MLJLinearModels)(),
                tuned_iterated_pipe]

errs = performance.(three_models)

(models=typeof.(three_models), mean_square_errors=errs) |> pretty

# So, apparently better than linear, but no match for the
# tree-booster. And our cv-estimate of performance on the
# train/validate dataset is wildly optimistic.


text:
abstract type InitializationStrategyGamma <: InitializationStrategy end

struct FixedGammaStrategy <: InitializationStrategyGamma
    kernel
end

calculate_gamma(model, strategy::FixedGammaStrategy) = MLKernels.getvalue(strategy.kernel.alpha)
function calculate_gamma(model::SubSVDD, strategy::FixedGammaStrategy, subspace_idx)
    MLKernels.getvalue(strategy.kernel[subspace_idx].alpha)
end

"""
Original publication:
Silverman, Bernard W. Density estimation for statistics and data analysis. Routledge, 2018.
"""
struct RuleOfThumbSilverman <: InitializationStrategyGamma end

function rule_of_thumb_silverman(data::Array{T,2}) where T <: Real
    return (size(data, 2) * (size(data, 1) + 2) / 4.0)^(-1.0 / (size(data,1) + 4.0))
end

calculate_gamma(model, strategy::RuleOfThumbSilverman) = rule_of_thumb_silverman(model.data)

function calculate_gamma(model::SubSVDD, strategy::RuleOfThumbSilverman, subspace_idx)
    return rule_of_thumb_silverman(model.data[model.subspaces[subspace_idx], :])
end

"""
Original publication:
Scott, David W. Multivariate density estimation: theory, practice, and visualization. John Wiley & Sons, 2015.
"""
struct RuleOfThumbScott <: InitializationStrategyGamma end

function rule_of_scott(data::Array{T,2}) where T <: Real
    return size(data, 2)^(-1.0/(size(data, 1) + 4))
end

calculate_gamma(model, strategy::RuleOfThumbScott) = rule_of_scott(model.data)

function calculate_gamma(model::SubSVDD, strategy::RuleOfThumbScott, subspace_idx)
    return rule_of_scott(model.data[model.subspaces[subspace_idx], :])
end

"""
Generate binary data to tune a one class classifier according to the following paper:
Wang, S. et al. 2018. Hyperparameter selection of one-class support vector machine by self-adaptive data
shifting. Pattern Recognition. 74, 2018.

Changes compared to the paper:
1. We only use the strategy here to fit gamma and not a combination of gamma and C.
2. We also fit on the real data. This is because the algorithm sometimes generates no pseudo inliers.
3. We test on pseudo-inliers, pseudo-outliers and real data. We assume that the real data are all inliers when testing.
4. We use the F1 score instead of an averaged FNR and FPR.
"""
struct WangGammaStrategy <: InitializationStrategyGamma
    solver
    gamma_search_range
    C
    scoring_function
end

WangGammaStrategy(solver) = WangGammaStrategy(solver, 1.0)
WangGammaStrategy(solver, C) = WangGammaStrategy(solver, 10.0.^range(-2, stop=2, length=50), C, f1_scoring)
WangGammaStrategy(solver, gamma_search_range, C) = WangGammaStrategy(solver, gamma_search_range, C, f1_scoring)

function generate_binary_data_for_tuning(data, k=nothing, threshold=0.1)
    if k === nothing
        k = round(Int, ceil(5 * log10(size(data, 2))))
    end
    tree = KDTree(data)
    edge_idx = Int[]
    norm_vec = []
    data_target = []
    l_ns = 0

    for i in 1:size(data, 2)
        idx, dist = knn(tree, data[:, i], k + 1, true)
        v_ij = mapslices(normalize, data[:, i] .- data[:, idx[2:end]], dims=1)
        n_i = sum(v_ij, dims=2)
        θ_ij = sum(v_ij .* n_i, dims=1)
        l_i = 1 / k * sum(θ_ij .>= 0)
        if l_i >= 1 - threshold
            # add new edge
            push!(edge_idx, i)
            push!(norm_vec, n_i)
            l_ns += 1 / k * sum(dist)
        end

        # generate pseudo inlier (shift direction is ∇p(x) = -n)
        n_i = -normalize(vec(n_i))
        Λ_i_positive = sum(n_i .* (data[:, idx[2:end]] .- data[:, i]), dims=1)
        if length(Λ_i_positive[Λ_i_positive .> 0]) > 0
            # shift along positive direction of data density gradient
            x_ij_min_positive = minimum(Λ_i_positive[Λ_i_positive .> 0])
            push!(data_target, data[:, i] + x_ij_min_positive * n_i)
        end
    end

    # compute negative shift amount
    l_ns *= 1 / length(edge_idx)
    # generate pseudo outliers by shifting along negative direction
    data_outliers = data[:, edge_idx] + mapslices(normalize, hcat(norm_vec...), dims=1) * l_ns

    data_target = hcat(data_target...)

    return data_target, data_outliers
end

function f1_scoring(predictions, ground_truth)
    return MLBase.f1score(MLBase.roc(ground_truth .== :outlier, predictions .== :outlier))
end

function calculate_gamma(model::SubSVDD, strategy::WangGammaStrategy, subspace_idx)
    info(LOGGER, "[Gamma Search] Using VanillaSVDD to search for gamma in subspace $subspace_idx")
    calculate_gamma(VanillaSVDD(model.data[model.subspaces[subspace_idx], :]), strategy)
end

function calculate_gamma(model::SubSVDD, strategy::WangGammaStrategy)
    info(LOGGER, "[Gamma Search] Using VanillaSVDD to estimate a global gamma.")
    calculate_gamma(VanillaSVDD(model.data), strategy)
end
function calculate_gamma(model, strategy::WangGammaStrategy)
    m = deepcopy(model)
    data_target, data_outliers = generate_binary_data_for_tuning(m.data)
    if isempty(data_target)
        data_predictions = hcat(m.data, data_outliers)
        ground_truth = vcat(fill(:inlier, size(m.data, 2)),
                            fill(:outlier, size(data_outliers, 2)))
    else
        data_predictions = hcat(m.data, data_target, data_outliers)
        ground_truth = vcat(fill(:inlier, size(m.data, 2) + size(data_target, 2)),
                            fill(:outlier, size(data_outliers, 2)))
    end

    debug(LOGGER, "[Gamma Search] Searching for parameter gamma.")
    best_gamma = 1.0
    best_score = -Inf
    for gamma in strategy.gamma_search_range
        debug(LOGGER, "[Gamma Search] Testing gamma = $gamma")
        init_strategy = FixedParameterInitialization(GaussianKernel(gamma), strategy.C)
        initialize!(m, init_strategy)
        set_adjust_K!(m, true)
        try
            fit!(m, strategy.solver)
        catch e
            debug(LOGGER, "[Gamma Search] Fitting failed for gamma $gamma.")
            debug(LOGGER, e)
            continue
        end
        predictions = classify.(predict(m, data_predictions))
        score = strategy.scoring_function(ground_truth, predictions)
        if score > best_score
            debug(LOGGER, "[Gamma Search] New best found with gamma = $gamma and score = $score.")
            best_gamma = gamma
            best_score = score
        end
    end
    info(LOGGER, "[Gamma Search] Finished.")
    return best_gamma
end


text:
# Load X and y variable
using JLD
data = load("groupData.jld")
(X,y,Xtest,ytest) = (data["X"],data["y"],data["Xtest"],data["ytest"])

# Fit multi-class logistic regression classifer
include("logReg.jl")
model = logRegSoftmax(X,y)

# Compute training and validation error
yhat = model.predict(X)
trainError = mean(yhat .!= y)
yhat = model.predict(Xtest)
validError = mean(yhat .!= ytest)

# Count number of parameters in model and number of features used
nModelParams = sum(model.w .!= 0)
nFeaturesUsed = sum(sum(abs.(model.w),2) .!= 0)
@show(trainError)
@show(validError)
@show(nModelParams)
@show(nFeaturesUsed)

# Show the image as a matrix
using PyPlot
imshow(model.w);


text:
module Optimisers

using Functors: functor, fmap, isleaf
using LinearAlgebra

include("interface.jl")

include("destructure.jl")
export destructure, total, total2

include("rules.jl")
export Descent, ADAM, Momentum, Nesterov, RMSProp,
       ADAGrad, AdaMax, ADADelta, AMSGrad, NADAM, ADAMW, RADAM, OADAM, AdaBelief,
       WeightDecay, ClipGrad, ClipNorm, OptimiserChain

"""
    Optimisers.apply!(rule::RuleType, state, parameters, gradient) -> (state, gradient)

This defines the action of any optimisation rule. It should return the modified gradient
which will be subtracted from the parameters, and the updated state (if any) for use at
the next iteration, as a tuple `(state, gradient)`.

For efficiency it is free to mutate the old state, but only what is returned will be used.
Ideally this should check `iswriteable(x)`, which the built-in rules do via [`@..`](@ref).

The initial state is `init(rule::RuleType, parameters)`.

# Example
```jldoctest
julia> Optimisers.init(Descent(0.1), Float32[1,2,3]) === nothing
true

julia> Optimisers.apply!(Descent(0.1), nothing, Float32[1,2,3], [4,5,6])
(nothing, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}}(*, ([4, 5, 6], 0.1f0)))
```
"""
apply!

"""
    Optimisers.init(rule::RuleType, parameters) -> state

Sets up the initial state for a given optimisation rule, and an array of parameters.
This and [`apply!`](@ref) are the two functions which any new optimisation rule must define.

# Examples
```jldoctest
julia> Optimisers.init(Descent(), Float32[1,2,3])  # is `nothing`

julia> Optimisers.init(Momentum(), [1.0, 2.0])
2-element Vector{Float64}:
 0.0
 0.0
```
"""
init

"""
    Optimisers.setup(rule, model) -> tree

Initialises the given optimiser for every trainable parameter within the model.
Returns a tree of the relevant states, which must be passed to [`update`](@ref)
or [`update!`](@ref).

# Example
```jldoctest
julia> Optimisers.setup(Descent(0.1f0), (x = rand(3), y = (true, false), z = tanh))
(x = Leaf(Descent{Float32}(0.1), nothing), y = (nothing, nothing), z = nothing)
```
"""
setup

"""
    Optimisers.update(tree, model, gradient) -> (tree, model)

Uses the optimiser and the gradient to change the trainable parameters in the model.
Returns the improved model, and the optimiser states needed for the next update.
The initial tree of states comes from [`setup`](@ref).

See also [`update!`](@ref), which will be faster for models of ordinary `Array`s or `CuArray`s.

# Example
```jldoctest
julia> m = (x = Float32[1,2,3], y = tanh);

julia> t = Optimisers.setup(Descent(0.1f0), m)
(x = Leaf(Descent{Float32}(0.1), nothing), y = nothing)

julia> g = (x = [1,1,1], y = nothing);  # fake gradient

julia> Optimisers.update(t, m, g)
((x = Leaf(Descent{Float32}(0.1), nothing), y = nothing), (x = Float32[0.9, 1.9, 2.9], y = tanh))
```
"""
update

"""
    Optimisers.update!(tree, model, gradient) -> (tree, model)

Uses the optimiser and the gradient to change the trainable parameters in the model.
Returns the improved model, and the optimiser states needed for the next update.
The initial tree of states comes from [`setup`](@ref).

This is used in exactly the same manner as [`update`](@ref), but because it may mutate
arrays within the old model (and the old state), it will be faster for models of ordinary
`Array`s or `CuArray`s. However, you should not rely on the old model being fully updated.
"""
update!

end # module


text:
import Base.Math: +, -, *, /
import Base: abs2

#generic case to apply a binary operation to two splines
function apply_binary_op{T}(f::Spline{T}, g::Spline{T}, op::Function)
	fdomain = domain(f)
	gdomain = domain(g)

	newdomain_low = max(fdomain[1], gdomain[1])
	newdomain_high = min(fdomain[end], gdomain[end])
	newdomain = (newdomain_low, newdomain_high)
	newknots = Array{T, 1}(0)

	for knot in union(f.knots, g.knots)
		if knot >= newdomain_low && knot <= newdomain_high
			push!(newknots, knot)
		end
	end

	newknots = sort(unique(newknots))
	@assert !any(isnan(newknots))

	newvalues = similar(newknots)

	for (i, knot) in enumerate(newknots)
		newvalues[i] = op(f(knot), g(knot))
	end

	return Spline(newknots, newvalues)
end

#Arithmetic for splines is defined just as for functions, since they are really approximations that can live in a finite computer to functions that cannot. Therefore, addition and scalar multiplication are defined just as you would expect

#f(x)+g(x)
function +{T}(f::Spline{T}, g::Spline{T})
	return apply_binary_op(f, g, (x, y)->x+y)
end

function -{T}(f::Spline{T}, g::Spline{T})
	return apply_binary_op(f, g, (x, y)->x-y)
end

#f(x)*g(x)
function *{T}(f::Spline{T}, g::Spline{T})
	return apply_binary_op(f, g, (x, y)->x*y)
end

#f(x)/g(x)
function /{T}(f::Spline{T}, g::Spline{T})
	return apply_binary_op(f,  g, (x, y)->x/y)
end

#f(x)+c
function +{T}(f::Spline{T}, c::Number)
	f_new = deepcopy(f)
	f_new.values += c
	f_new.a += c

	return f_new
end

#c+f(x)
function +{T}(c::Number, f::Spline{T})
	return f+c
end

#c*f(x)
function *{T}(c::Number, f::Spline{T})
	f_new = deepcopy(f)
	f_new.values *= c
	f_new.a *= c
	f_new.b *= c
	f_new.c *= c
	f_new.d *= c

	return f_new
end

#f(x)*c
function *{T}(f::Spline{T}, c::Number)
	return c*f
end

#f(x)/c
function /{T}(f::Spline{T}, c::Number)
	f_new = deepcopy(f)
	f_new.values /= c
	f_new.a /= c
	f_new.b /= c
	f_new.c /= c
	f_new.d /= c

	return f_new
end

function abs2{T<:Real}(f::Spline{T})
	return f^2
end

###   Now we define function composition for splines and functions, introducing the ∘ operator that does this   ###
function compose{T}(f::Function, g::Spline{T})
	newknots = adaptive_discretize_mesh(g)
	newvalues = similar(newknots)

	for (i, knot) in enumerate(newknots)
		newvalues[i] = f(g(knot))
	end

	return Spline(newknots, newvalues)
end

function ∘{T}(f::Function, g::Spline{T})
	return compose(f, g)
end

export compose
export ∘


text:
module Cntr_Test

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOSmoothMcCormickGrad

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type FlashPHSteady
	FlashPHSteady()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(),
			DanaReal(Dict{Symbol,Any}(
				:Default=>1000,
				:Brief=>"Regularization Factor"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Feed Stream",
				:PosX=>0,
				:PosY=>0.5421,
				:Symbol=>"_{in}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream",
				:PosX=>0.4790,
				:PosY=>1,
				:Symbol=>"_{outL}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream",
				:PosX=>0.4877,
				:PosY=>0,
				:Symbol=>"_{outV}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Rate of heat supply",
				:PosX=>1,
				:PosY=>0.7559,
				:Symbol=>"_{in}"
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Vaporization fraction",
				:Symbol=>"\\phi"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Lower=>-0.1,
				:Upper=>1.1,
				:Brief=>"Vaporization fraction if saturated",
				:Symbol=>"\\phi_{sat}"
			)),
			temperature(Dict{Symbol,Any}(
				:Lower=>173,
				:Upper=>1473,
				:Brief=>"Temperature if saturated"
			)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Lower=>0,
				:Upper=>1,
				:Brief=>"Liquid composition if saturated"
			)),(NComp)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Lower=>0,
				:Upper=>1,
				:Brief=>"Vapour composition if saturated"
			)),(NComp)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Pressure Ratio",
				:Symbol=>"P_{ratio}"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P"
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Regularization Variable"
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Regularization Variable"
			)),
			[
				:(PP.LiquidFugacityCoefficient(Tsat, OutletL.P, xsat)*xsat = PP.VapourFugacityCoefficient(Tsat, OutletV.P, ysat)*ysat),
				:(Inlet.F = OutletV.F + OutletL.F),
				:(OutletV.F = Inlet.F * vfrac),
				:(Inlet.F*Inlet.z = OutletL.F*xsat + OutletV.F*ysat),
				:(sum(xsat) = sum(ysat)),
				:(Inlet.F*Inlet.h + InletQ = Inlet.F*(1-vsat)*PP.LiquidEnthalpy(Tsat, OutletL.P, xsat) + Inlet.F*vsat*PP.VapourEnthalpy(Tsat, OutletV.P, ysat)),
				:(Inlet.F*Inlet.h + InletQ = Inlet.F*(1-vfrac)*OutletL.h + Inlet.F*vfrac*OutletV.h),
				:(OutletV.T = OutletL.T),
				:(OutletV.P = OutletL.P),
				:(OutletL.P = Inlet.P - Pdrop),
				:(OutletL.P = Inlet.P * Pratio),
				:(zero_one = (1 + tanh(B * vsat))/2),
				:(one_zero = (1 - tanh(B * (vsat - 1)))/2),
				:(vfrac = zero_one * one_zero * vsat + 1 - one_zero),
				:(OutletL.z = zero_one*one_zero*xsat + (1-zero_one*one_zero)*Inlet.z),
				:(OutletV.z = zero_one*one_zero*ysat + (1-zero_one*one_zero)*Inlet.z),
			],
			[
				"Chemical equilibrium","Global Molar Balance","","Component Molar Balance","","Energy Balance if saturated","Real Energy Balance","Thermal Equilibrium","Mechanical Equilibrium","Pressure Drop","Pressure Ratio","","","","","",
			],
			[:PP,:NComp,:B,],
			[:Inlet,:OutletL,:OutletV,:InletQ,:vfrac,:vsat,:Tsat,:xsat,:ysat,:Pratio,:Pdrop,:zero_one,:one_zero,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	B::DanaReal
	Inlet::stream
	OutletL::liquid_stream
	OutletV::vapour_stream
	InletQ::power
	vfrac::fraction
	vsat::DanaReal
	Tsat::temperature
	xsat::Array{DanaReal}
	ysat::Array{DanaReal}
	Pratio::positive
	Pdrop::press_delta
	zero_one::fraction
	one_zero::fraction
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export FlashPHSteady
function setEquationFlow(in::FlashPHSteady)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	# regularization functions
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
end
function atributes(in::FlashPHSteady,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/Flash"
	fields[:Brief]="Another model of a static PH flash."
	fields[:Info]="This model shows how to model a pressure enthalpy flash
directly with the EMSO modeling language.

This model is for demonstration purposes only, the flashPH
routine available on VRTherm is much more robust.

== Assumptions ==
* perfect mixing of both phases;

== Specify ==
* the feed stream;
* the heat duty;
* the outlet pressure.
"
	drive!(fields,_)
	return fields
end
FlashPHSteady(_::Dict{Symbol,Any})=begin
	newModel=FlashPHSteady()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
function transform_date_fields(data_frame::DataFrame, col_key::Symbol)

    # creates a col of string dates -
    date_col = string.(data_frame[!, col_key])
    tmp_array = Array{Date,1}()
    for (_, date_value) in enumerate(date_col)

        # split around "/"
        date_components = split(date_value, "/")

        # process the month -
        month_value = date_components[1]
        if (length(month_value) < 2)
            month_value = lpad(month_value, 1, "0")
        end

        # process the day -
        day_value = date_components[2]
        if (length(day_value) < 2)
            day_value = lpad(day_value, 1, "0")
        end

        # process the year -
        year_value = date_components[3]
        if (length(year_value) < 3)
            year_value = "20$(year_value)"
        end

        # new date -
        new_date_string = "$(year_value)-$(month_value)-$(day_value)"
        new_date_object = Dates.Date(new_date_string)
        push!(tmp_array, new_date_object)
    end

    # remove the col_key -
    select!(data_frame, Not(col_key))
    insertcols!(data_frame, col_key => tmp_array)
end

function filter(data::DataFrame, start_date::Date, end_date::Date;
    key::Symbol = :timestamp)::DataFrame

    return filter(key => (x) -> (x>=start_date && x<=end_date), data);
end

function extract_data_block_for_date_range(data_table::DataFrame, start_date::Date, end_date::Date;
    key::Symbol = :timestamp)::DataFrame

    # find the range -
    idx_range = findall(x->(x>=start_date && x<=end_date), data_table[!, key]) |> collect

    # extract the range -
    data_block = data_table[idx_range,:]

    # return -
    return data_block
end

function reverse_row_order_in_table(table::DataFrame)

    # how many rows do we have?
    (number_of_rows, _) = size(table);

    # create a data array -
    tmp_date_array = Array{Date,1}()
    for row_index ∈ 1:number_of_rows
    end

    df_tmp = DataFrame()
    for row_index ∈ number_of_rows:-1:1
        df_row = table[row_index,:]
        push!(df_tmp, df_row)
    end

    return df_tmp
end

max_stars_repo_name: Pooksoft/Serenity

max_issues_repo_name: Alba-Intelligence/Serenity

max_forks_repo_name: Alba-Intelligence/Serenity


text:
using ParticlePlots
using Test

@testset "ParticlePlots.jl" begin
    # Write your tests here.
end


text:
GradientDescent(execution::ExecutionPolicy; γ = 1.) = ProxGradient(execution,Boosting.None(),Step.Constant(; γ = γ),Smoothing.None(),Prox.None())
GradientDescent(; kw...) = GradientDescent(Execution.Serial(); kw...)

MomentumDescent(execution::ExecutionPolicy; μ = 0.9, ϵ = 1e-1, γ = 1.) = ProxGradient(execution,Boosting.Momentum(; μ = μ, ϵ = ϵ),Step.Constant(; γ = γ),Smoothing.None(),Prox.None())
MomentumDescent(; kw...) = MomentumDescent(Execution.Serial(); kw...)

NesterovDescent(execution::ExecutionPolicy; μ = 0.9, ϵ = 1e-1, γ = 1.) = ProxGradient(execution,Boosting.Nesterov(; μ = μ, ϵ = ϵ),Step.Constant(; γ = γ),Smoothing.None(),Prox.None())
NesterovDescent(; kw...) = NesterovDescent(Execution.Serial(); kw...)

AdaGrad(execution::ExecutionPolicy; γ = 1., ϵ = 1e-6) = ProxGradient(execution,Boosting.None(),Step.Constant(; γ = γ),Smoothing.Adagrad(; ϵ = ϵ),Prox.None())
AdaGrad(; kw...) = AdaGrad(Execution.Serial(); kw...)

AdaDelta(execution::ExecutionPolicy; γ = 1., ρ = 0.95, ϵ = 1e-6) = ProxGradient(execution,Boosting.None(),Step.Constant(; γ = γ),Smoothing.Adadelta(; ρ = ρ, ϵ = ϵ),Prox.None())
AdaDelta(; kw...) = AdaDelta(Execution.Serial(); kw...)

Adam(execution::ExecutionPolicy; μ = 0.99, ϵ = 1e-2, γ = 1., ρ = 0.999, ϵ_rms = 1e-8) = ProxGradient(execution,Boosting.Momentum(; μ = μ, ϵ = ϵ),Step.Constant(; γ = γ),Smoothing.RMSprop(; ρ = ρ, ϵ = ϵ_rms),Prox.None())
Adam(; kw...) = Adam(Execution.Serial(); kw...)

Nadam(execution::ExecutionPolicy; μ = 0.99, ϵ = 1e-2, γ = 1., ρ = 0.999, ϵ_rms = 1e-8) = ProxGradient(execution,Boosting.Nesterov(; μ = μ, ϵ = ϵ),Step.Constant(; γ = γ),Smoothing.RMSprop(; ρ = ρ, ϵ = ϵ_rms),Prox.None())
Nadam(; kw...) = Nadam(Execution.Serial(); kw...)

BB(execution::ExecutionPolicy; γ₀ = 1.) = ProxGradient(execution,Boosting.None(),Step.BB(; γ₀ = γ₀),Smoothing.None(),Prox.None())
BB(; kw...) = BB(Execution.Serial(); kw...)


text:
module TestEnsembleNest

using Test: @test, @testset

using Ensemble
using .EnsembleNest

using Statistics

function test5DGaussian()
    nd = 5

    function logp(x)
        if any(x .> 1.0) || any(x .< 0.0)
            -Inf
        else
            0.0
        end
    end

    mu = 0.5.*ones(nd) .+ 0.1.*rand(Float64, (nd,))
    sigma = 0.01.*ones(nd)

    function logl(x)
        -0.5*nd*log(2.0*pi) - sum(log.(sigma)) - 0.5*sum((x-mu).*(x-mu)./(sigma.*sigma))
    end

    nmcmc = 100
    nl = 1024
    xs = rand(Float64, (nd, nl))
    ns = NestState(logl, logp, xs, nmcmc)

    run!(ns, 0.1)

    lZ, dlZ = logZ(ns)

    @testset "5D Gaussian tests" begin
        @test isapprox(lZ, 0.0, atol=1.0)

        post, lnprobs = postsample(ns)
        npost = size(post, 2)

        for i in 1:nd
            smu = mean(post[i,:])
            @test isapprox(smu, mu[i], atol=0.1*sigma[i])
        end
    end
end

function testall()
    @testset "EnsembleNest tests" begin
        test5DGaussian()
    end
end

end


text:
#Code to do type I error simulation for TrajGWAS.jl
# after reworking SPA

using Interpolations
import Optim: minimizer, optimize, LBFGS, NelderMead

function ecgf(z::AbstractVector; knots::Integer=10000)

    range = (-200.0, 200.0)
    idx0 = quantile.(Cauchy(), (1:knots)/(knots+1))
    idx1 = idx0 .* maximum(range) / maximum(idx0)

    K0_int, K1_int, K2_int = ecgf(idx1, z)

    return K0_int, K1_int, K2_int
end

function ecgf(vs::AbstractVector, zs::AbstractVector)
    exp_zv = Vector{Float64}(undef, length(zs))
    zexp_zv = Vector{Float64}(undef, length(zs))
    z2exp_zv = Vector{Float64}(undef, length(zs))
    K0 = Vector{Float64}(undef, length(vs))
    K1 = similar(K0)
    K2 = similar(K0)
    for (i, v) in enumerate(vs)
        exp_zv .= zs .* v
        alpha = maximum(exp_zv)
        exp_zv .= exp.(exp_zv .- alpha)

        zexp_zv .= zs .* exp_zv
        z2exp_zv .= zs .^ 2 .* exp_zv
        M0 = mean(exp_zv)
        M1 = mean(zexp_zv)
        M2 = mean(z2exp_zv)
        K0[i] = log(M0) + alpha
        K1[i] = M1 / M0
        K2[i] = M2 / M0 - (M1 / M0) ^ 2
    end
    K0_int = LinearInterpolation(vs, convert(Vector{Float64}, K0);
        extrapolation_bc=Line())
    K1_int = LinearInterpolation(vs, convert(Vector{Float64}, K1);
        extrapolation_bc=Line())
    K2_int = LinearInterpolation(vs, convert(Vector{Float64}, K2);
        extrapolation_bc=Line())
    return K0_int, K1_int, K2_int
end

struct TrajGWASEcgfCollection
    K0_β
    K1_β
    K2_β
    K0_τ
    K1_τ
    K2_τ
end

function ecgf(st::WSVarScoreTestInvariant)
    K0_β, K1_β, K2_β = ecgf(st.ψ_β1_pre)
    K0_τ, K1_τ, K2_τ = ecgf(st.ψ_τ1_pre)
    TrajGWASEcgfCollection(K0_β, K1_β, K2_β, K0_τ, K1_τ, K2_τ)
end

"""
    spa(g, st, p_alt, Ks; mode, g_norm, tmp_ecgf, r)

Perform saddlepoint approximation for a single variant.

## arguments
- `g`::AbstractVector:  genotype vector
- `st::WSVarScoreTestInvariant`
- `p_alt` : alternative p-value if saddlepoint approximation is skipped due to cutoff `r`
- `Ks::TrajGWASEcgfCollection`
- `mode`: `:bed` for hard calls of 0, 1, 2; `:ukbbgen` for bgen format provided by UK Biobank; anything else for others
"""
function spa(g::AbstractVector,
    st::WSVarScoreTestInvariant,
    p_alt,
    dir_alt,
    Ks::TrajGWASEcgfCollection;
    g_norm=similar(g),
    cnts=nothing,
    ref_vals= cnts !== nothing ? similar(g, length(cnts)) : nothing,
    vals_norm= cnts !== nothing ? similar(g, length(cnts)) : g_norm,
    tmp_ecgf = cnts !== nothing ? similar(g, length(cnts)) : similar(g),
    r=2.0,
    adjustor::Union{Adjustor, Nothing}=nothing,
    adj_cutoff::Real=5e-5
    )
    # involves internal normalization
    m = mean(skipmissing(g))
    #s = std(skipmissing(g))

    # screening for β
    g_norm .= g .- m
    @inbounds for i in 1:length(g_norm)
        g_norm[i] = isnan(g_norm[i]) ? 0.0 : g_norm[i]
    end
    cutoff_factor = sum(x -> x^2, g_norm)
    std_β = sqrt(st.var_β1_pre * cutoff_factor)
    score_β = dot(g_norm, st.ψ_β1_pre)
    z_β = score_β/std_β

    if abs(z_β) < r
        p_β = p_alt[1]
        dir_β = dir_alt[1]
    else
        g_norm .= (g .- m) ./ std_β
        if ref_vals !== nothing
            vals_norm .= (ref_vals .- m) ./ std_β
            @inbounds for i in 1:length(vals_norm)
                if isnan(vals_norm[i])
                    vals_norm[i] = zero(eltype(vals_norm))
                end
            end
        end
        p_β, dir_β = spa(g_norm, st.ψ_β1_pre, vals_norm,
            cutoff_factor, 0.0, p_alt[1], dir_alt[1], Ks.K0_β, Ks.K1_β, Ks.K2_β;
            cnts=cnts, tmp_ecgf=tmp_ecgf, pre_vec_var=st.var_β1_pre,
            adjustor=adjustor, adj_cutoff=adj_cutoff)
    end

    # screening for τ

    g_norm .= g .- m
    @inbounds for i in 1:length(g_norm)
        g_norm[i] = isnan(g_norm[i]) ? 0.0 : g_norm[i]
    end
    # cutoff_factor = sum(x -> x^2, g_norm)
    std_τ = sqrt(st.var_τ1_pre * cutoff_factor)
    score_τ = dot(g_norm, st.ψ_τ1_pre)
    z_τ = score_τ/std_τ

    if abs(z_τ) < r
        p_τ = p_alt[2]
        dir_τ = dir_alt[2]
    else
        g_norm .= (g .- m) ./ std_τ
        if ref_vals !== nothing
            vals_norm .= (ref_vals .- m) ./ std_τ
            @inbounds for i in 1:length(vals_norm)
                if isnan(vals_norm[i])
                    vals_norm[i] = zero(eltype(vals_norm))
                end
            end
        end
        p_τ, dir_τ = spa(g_norm, st.ψ_τ1_pre, vals_norm,
            cutoff_factor, 0.0, p_alt[2], dir_alt[2], Ks.K0_τ, Ks.K1_τ, Ks.K2_τ;
            cnts=cnts, tmp_ecgf=tmp_ecgf, pre_vec_var=st.var_τ1_pre,
            adjustor=adjustor, adj_cutoff=adj_cutoff)
    end
    # p_β, dir_β = spa(g_norm, st.ψ_β1_pre, vals_norm,
    #     cutoff_factor, r, p_alt[1], dir_alt[1], Ks.K0_β, Ks.K1_β, Ks.K2_β;
    #     cnts=cnts, tmp_ecgf=tmp_ecgf, pre_vec_var=st.var_β1_pre,
    #     adjustor=adjustor, adj_cutoff=adj_cutoff)
    # p_τ, dir_τ = spa(g_norm, st.ψ_τ1_pre, vals_norm,
    #     cutoff_factor, r, p_alt[2], dir_alt[2], Ks.K0_τ, Ks.K1_τ, Ks.K2_τ;
    #     cnts=cnts, tmp_ecgf=tmp_ecgf, pre_vec_var=st.var_τ1_pre,
    #     adjustor=adjustor, adj_cutoff=adj_cutoff)

    p_β, p_τ, dir_β, -dir_τ
end

function spa(g_norm::AbstractVector, pre_vec::AbstractVector,
    vals_norm::AbstractVector,
    cutoff_factor::Real, r::Real, p_alt::Real, dir_alt, K0, K1, K2;
    cnts=nothing,
    tmp_ecgf = g_norm === nothing ? similar(g_norm) : similar(g_norm, length(cnts)),
    pre_vec_var::Real=var(pre_vec),
    adjustor::Union{Adjustor, Nothing}=nothing,
    adj_cutoff::Real=1e-3
    )
    s = dot(g_norm, pre_vec)
    cutoff = r * sqrt(pre_vec_var * cutoff_factor)
    pval, dir = _get_pval(s, cutoff, p_alt, dir_alt, vals_norm, K0, K1, K2, tmp_ecgf; cnts=cnts)
    if pval < adj_cutoff && adjustor !== nothing
        g_norm .= g_norm .- adjustor.X_XtXinv * (transpose(adjustor.X) * g_norm)
        g_norm /= sqrt(sum(x -> x^2, g_norm) * pre_vec_var) # renormalization
        s = dot(g_norm, pre_vec)
        pval, dir = _get_pval(s, Inf, 0.0, 0, g_norm, K0, K1, K2, tmp_ecgf; cnts=nothing)
    end
    return pval, dir
end

const normal = Normal()
const lbfgs = LBFGS()
const singleton_zero = [0.0]
function _get_pval(s, cutoff, p_alt, dir_alt, vals_norm, K0, K1, K2, tmp_ecgf; cnts=nothing)
    function K0_(z)
        if cnts === nothing
            @inbounds for i in 1:length(vals_norm)
                tmp_ecgf[i] = K0(vals_norm[i] * z)
            end
            return sum(tmp_ecgf)
        else
            fill!(tmp_ecgf, 0.0)
            @inbounds for i in 1:length(tmp_ecgf)
                if cnts[i] != 0
                    tmp_ecgf[i] = K0(vals_norm[i] * z)
                end
            end
            return dot(cnts, tmp_ecgf)
        end
    end
    function K1_(z)
        if cnts === nothing
            @inbounds for i in 1:length(vals_norm)
                tmp_ecgf[i] = vals_norm[i] * K1(vals_norm[i] * z)
            end
            return sum(tmp_ecgf)
        else
            fill!(tmp_ecgf, 0.0)
            @inbounds for i in 1:length(tmp_ecgf)
                if cnts[i] != 0
                    tmp_ecgf[i] = vals_norm[i] * K1(vals_norm[i] * z)
                end
            end
            return dot(cnts, tmp_ecgf)
        end
    end
    function K2_(z)
        if cnts === nothing
            @inbounds for i in 1:length(vals_norm)
                tmp_ecgf[i] = vals_norm[i] ^ 2 * K2(vals_norm[i] * z)
            end
            return sum(tmp_ecgf)
        else
            fill!(tmp_ecgf, 0.0)
            @inbounds for i in 1:length(tmp_ecgf)
                if cnts[i] != 0
                    tmp_ecgf[i] = vals_norm[i] ^ 2 * K2(vals_norm[i] * z)
                end
            end
            return dot(cnts, tmp_ecgf)
        end
    end
    function f(x)
        r = K0_(x[1]) - x[1] * s
        return r
    end
    function g!(storage, x)
        storage[1] = K1_(x[1]) - s
    end
    function h!(storage, x)
        storage[1] = K2_(x[1])
    end
    if abs(s) < cutoff
        return p_alt, dir_alt
    else
        r = optimize(f, g!, h!, singleton_zero, lbfgs)
        zeta = minimizer(r)[1]
        omega = sign(zeta) * sqrt(max(0.0, 2 * (zeta * s - K0_(zeta))))
        nu = zeta * sqrt(max(0, K2_(zeta)))
        z2 = omega + 1.0/omega * log(nu / omega)

        return ccdf(normal, abs(z2)) + cdf(normal, -abs(z2)), isnan(z2) ? 0 : convert(Int, sign(z2))
    end
end