text:
using Test, AdmittanceModels, LinearAlgebra
import .CircuitExample, .CircuitCascadeExample
const ce = CircuitExample
const cce = CircuitCascadeExample

@testset "Circuit" begin
    @test length(ce.circuit.vertices) == 5
    for m in matrices(ce.circuit)
        @test size(m) == (5,5)
        @test issymmetric(m)
    end
    @test get_inv_inductance(ce.circuit, ce.Q0_SIGNAL, ce.Q0_GROUND) == .1
    @test get_inv_inductance(ce.circuit, ce.Q0_GROUND, ce.Q0_SIGNAL) == .1
    @test get_conductance(ce.circuit, ce.Q0_SIGNAL, ce.Q0_GROUND) == 0
    @test get_conductance(ce.circuit, ce.Q0_GROUND, ce.Q0_SIGNAL) == 0
    @test get_capacitance(ce.circuit, ce.Q0_SIGNAL, ce.Q0_GROUND) == 100
    @test get_capacitance(ce.circuit, ce.Q0_GROUND, ce.Q0_SIGNAL) == 100
end

@testset "coordinate_matrix" begin
    @test coordinate_matrix(3) == [[0,0] [1,0] [0,1]]
    @test coordinate_matrix(ce.circuit) == coordinate_matrix(5)
    coordinate_matrix(ce.circuit, SpanningTree(ce.circuit.vertices)) == coordinate_matrix(5)
    correct_matrix = [[1, 1, 0, 0] [0, 1, 0, 0] [0, 0, 1, 1] [0, 0, 0, 1] [0, 0, 0, 0]]
    @test coordinate_matrix(ce.circuit, ce.tree) == correct_matrix
end

@testset "cascade and unite_vertices" begin
    circ = cascade(cce.circuit0, cce.circuit1, cce.circuit2)
    @test circ.vertices == 0:11
    l = length(circ.vertices)
    @test size(circ.c) == (l, l)
    united_circ = unite_vertices(circ, 0, 4, 8)
    united_circ = unite_vertices(united_circ, 1, 5, 9)
    @test united_circ.vertices == [0, 1, 2, 3, 6, 7, 10, 11]
    l = length(united_circ.vertices)
    @test size(united_circ.c) == (l, l)
    @test get_capacitance(united_circ, 2, 3) == get_capacitance(circ, 2, 3)
    @test get_inv_inductance(united_circ, 6, 1) == get_inv_inductance(circ, 6, 5)
    s = get_conductance(circ, 0, 1) + get_conductance(circ, 4, 5) + get_conductance(circ, 8, 9)
    @test get_conductance(united_circ, 0, 1) == s
    united_circ = unite_vertices(united_circ, 0, 1)
    s = get_capacitance(circ, 3, 0) + get_capacitance(circ, 3, 1)
    @test get_capacitance(united_circ, 3, 0) == s
end


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

"`NIDAQ.jl` provides an interface to NI-DAQmx--- National Instruments' driver for their data acquisition boards"
module NIDAQ

import Base.write, Base.read, Base.start

# tasks
export stop, clear

# channels
export analog_input, analog_output, digital_input, digital_output
export count_edges, quadrature_input, line_to_line, generate_pulses

# properties
export devices, channel_type, getproperties, setproperty!
export analog_input_ranges, analog_output_ranges
export analog_input_channels, analog_output_channels
export digital_input_channels, digital_output_channels
export counter_input_channels, counter_output_channels

const NIDAQmx = "C:\\Windows\\System32\\nicaiu.dll"

bitstype 32 Bool32 <: Integer
export Bool32

try
  global ver
  major = Ref{UInt32}(0)
  ccall((:DAQmxGetSysNIDAQMajorVersion,NIDAQmx),Int32,(Ref{UInt32},),major)
  minor = Ref{UInt32}(0)
  ccall((:DAQmxGetSysNIDAQMinorVersion,NIDAQmx),Int32,(Ref{UInt32},),minor)
  update = Ref{UInt32}(0)
  ccall((:DAQmxGetSysNIDAQUpdateVersion,NIDAQmx),Int32,(Ref{UInt32},),update)
  ver = "$(major[]).$(minor[]).$(update[])"
catch
  error("can not determine NIDAQmx version.")
end

try
  include("constants_V$ver.jl")
  include("functions_V$ver.jl")
catch
  error("NIDAQmx version $ver is not supported.")
end

unsigned_constants = Dict{UInt64,Symbol}()
signed_constants = Dict{Int64,Symbol}()

for sym in names(NIDAQ,true)
    isdefined(sym) || continue
    sym_str = string(sym)
    (length(sym_str)<5 || sym_str[1:5]!="DAQmx") && continue
    sym_str = sym_str[6:end]
    sym_str[1]=='_' && (sym_str = sym_str[2:end])
    if @eval typeof($sym) <: Unsigned
        @eval $(symbol(sym_str)) = UInt32($sym)
        unsigned_constants[eval(:($sym))] = symbol(sym_str)
    elseif @eval typeof($sym) <: Signed
        sym_str[1:min(end,4)]=="Val_" || continue
        @eval $(symbol(sym_str)) = convert(Int32,$sym)
        signed_constants[eval(:($sym))] = symbol(sym_str)
    elseif eval(:(typeof($sym)==Function))
        @eval $(symbol(sym_str)) = $sym
    end
end

function catch_error(code::Int32, extra::ASCIIString=""; err_fcn=error)
    sz = DAQmxGetErrorString(code, convert(Ptr{UInt8},C_NULL), convert(UInt32,0))
    data = zeros(UInt8,sz)
    ret = DAQmxGetErrorString(code, pointer(data), convert(UInt32,sz))
    ret>0 && warn("DAQmxGetErrorString error $ret")
    ret<0 && err_fcn("DAQmxGetErrorString error $ret")
    data = chop(convert(ASCIIString, data))
    code>0 && warn("NIDAQmx: "*extra*data)
    code<0 && err_fcn("NIDAQmx: "*extra*data)
end

include("task.jl")
include("analog.jl")
include("digital.jl")
include("counter.jl")
include("properties.jl")

@doc """`read(task, precision, nsamples) -> Matrix`

receive data from all analog or digital channels in a NIDAQ task
""" read

@doc """`write(task, data)`

send data to all analog or digital channels in a NIDAQ task
""" write

end

ext: jl


text:
# ---
# title: 1431. Kids With the Greatest Number of Candies
# id: problem1431
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/kids-with-the-greatest-number-of-candies/description/>
# hidden: true
# ---
#
# Given the array `candies` and the integer `extraCandies`, where `candies[i]`
# represents the number of candies that the **_ith_** kid has.
#
# For each kid check if there is a way to distribute `extraCandies` among the
# kids such that he or she can have the **greatest** number of candies among
# them. Notice that multiple kids can have the **greatest** number of candies.
#
#
#
# **Example 1:**
#
#
#
#     Input: candies = [2,3,5,1,3], extraCandies = 3
#     Output: [true,true,true,false,true]
#     Explanation:
#     Kid 1 has 2 candies and if he or she receives all extra candies (3) will have 5 candies --- the greatest number of candies among the kids.
#     Kid 2 has 3 candies and if he or she receives at least 2 extra candies will have the greatest number of candies among the kids.
#     Kid 3 has 5 candies and this is already the greatest number of candies among the kids.
#     Kid 4 has 1 candy and even if he or she receives all extra candies will only have 4 candies.
#     Kid 5 has 3 candies and if he or she receives at least 2 extra candies will have the greatest number of candies among the kids.
#
#
# **Example 2:**
#
#
#
#     Input: candies = [4,2,1,1,2], extraCandies = 1
#     Output: [true,false,false,false,false]
#     Explanation: There is only 1 extra candy, therefore only kid 1 will have the greatest number of candies among the kids regardless of who takes the extra candy.
#
#
# **Example 3:**
#
#
#
#     Input: candies = [12,1,12], extraCandies = 10
#     Output: [true,false,true]
#
#
#
#
# **Constraints:**
#
#   * `2 <= candies.length <= 100`
#   * `1 <= candies[i] <= 100`
#   * `1 <= extraCandies <= 50`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
## Find weakly connected components ##

@inline function _root(parent, i)
    while i != parent[i]
        i = parent[i] = parent[parent[i]]
    end
    return i
end

"""
    components(g::Graph)

Enumerate all connected components. Returns `Vector` with component affiliations
of nodes. The number of each component is also an index of one of the nodes
that is included in it.

Implementation based on WQUPC algorithm from Robert Sedgewick's lecture
on Union-Find algorithms
[link](https://www.cs.princeton.edu/~rs/AlgsDS07/01UnionFind.pdf).
"""
function components(g::Graph{N,E}) where {N,E}
    n      = nodecount(g)
    parent = collect(1:n)
    weight = fill(1, n)
    @inbounds for i = 1:n, j = outneighbors(g, i)
        src = _root(parent, i)          # union all
        tgt = _root(parent, Int(j))
        src == tgt && continue
        if weight[src] < weight[tgt]
            parent[src] = tgt
            weight[tgt] += weight[src]
            weight[src] = 0
        else
            parent[tgt] = src
            weight[src] += weight[tgt]
            weight[tgt] = 0
        end
    end
    for i = 1:n                         # find and root all
        while weight[parent[i]] == 0
            parent[i] = parent[parent[i]]
        end
    end
    return relabel(parent)
end

# TODO: implement giantcomponent and componentcount


text:
export UpExpFlow, ExpFlowCircuit, exp_pass_up, ExpectationUpward

#####################
# Expectation Flow circuits
# For use of algorithms depending on pairs of nodes of two circuits
#####################

"A expectation circuit node that has pair of origins of type PC and type LC"
abstract type ExpFlowNode{F} end

const ExpFlowCircuit{O} = Vector{<:ExpFlowNode{<:O}}

struct UpExpFlow{F} <: ExpFlowNode{F}
    p_origin::ProbCircuit
    f_origin::LogisticCircuit
    children::Vector{<:ExpFlowNode{<:F}}
    f::F
    fg::F
end

import LogicCircuits: children
@inline children(x::UpExpFlow) = x.children

"""
Expected Prediction of LC w.r.t PC.
This implementation uses the computation graph approach.
"""
function ExpectationUpward(pc::ProbCircuit, lc::LogisticCircuit, data)
    # 1. Get probability of each observation
    log_likelihoods = marginal(pc, data)
    p_observed = exp.( log_likelihoods )

    # 2. Expectation w.r.t. P(x_m, x_o)
    exps_flow = exp_pass_up(pc, lc, data)
    results_unnormalized = exps_flow[end].fg

    # 3. Expectation w.r.t P(x_m | x_o)
    results = transpose(results_unnormalized) ./ p_observed

    # 4. Add Bias terms
    biases = lc.thetas
    results .+= biases

    results, exps_flow
end

"""
Construct a upward expectation flow circuit from a given pair of PC and LC circuits
Note that its assuming the two circuits share the same vtree
"""
function ExpFlowCircuit(pc::ProbCircuit, lc::LogisticCircuit, batch_size::Int, ::Type{El}) where El
    F = Array{El, 2}
    fmem = () -> zeros(El, 1, batch_size) #Vector{El}(undef, batch_size)  #init_array(El, batch_size) # note: fmem's return type will determine type of all UpFlows in the circuit (should be El)
    fgmem = () -> zeros(El, num_classes(lc), batch_size)

    root_pc = pc
    root_lc = children(lc)[1]

    cache = Dict{Pair{DAG, DAG}, ExpFlowNode}()
    sizehint!(cache, (num_nodes(pc) + num_nodes(lc))*10)
    expFlowCircuit = Vector{ExpFlowNode}()

    function ExpflowTraverse(n::PlainSumNode, m::Logistic⋁Node)
        get!(cache, Pair(n, m)) do
            ch = [ ExpflowTraverse(i, j) for i in children(n) for j in children(m)]
            node = UpExpFlow{F}(n, m, ch, fmem(), fgmem())
            push!(expFlowCircuit, node)
            return node
        end
    end
    function ExpflowTraverse(n::PlainMulNode, m::Logistic⋀Node)
        get!(cache, Pair(n, m)) do
            ch = [ ExpflowTraverse(z[1], z[2]) for z in zip(children(n), children(m)) ]
            node = UpExpFlow{F}(n, m, ch, fmem(), fgmem())
            push!(expFlowCircuit, node)
            return node
        end
    end
    function ExpflowTraverse(n::PlainProbLiteralNode, m::Logistic⋁Node)
        get!(cache, Pair(n, m)) do
            ch = Vector{ExpFlowNode{F}}() # TODO
            node = UpExpFlow{F}(n, m, ch, fmem(), fgmem())
            push!(expFlowCircuit, node)
            return node
        end
    end
    function ExpflowTraverse(n::PlainProbLiteralNode, m::LogisticLiteralNode)
        get!(cache, Pair(n, m)) do
            ch = Vector{ExpFlowNode{F}}() # TODO
            node = UpExpFlow{F}(n, m, ch, fmem(), fgmem())
            push!(expFlowCircuit, node)
            return node
        end
    end

    ExpflowTraverse(root_pc, root_lc)
    expFlowCircuit
end

function exp_pass_up(pc::ProbCircuit, lc::LogisticCircuit, data)
    expFlowCircuit = ExpFlowCircuit(pc, lc, num_examples(data), Float32);
    for n in expFlowCircuit
        exp_pass_up_node(n, data)
    end
    expFlowCircuit
end

function exp_pass_up(fc::ExpFlowCircuit, data)
    #TODO write resize_flows similar to flow_circuits
    #     and give as input the expFlowCircuit instead
    #expFlowCircuit = ExpFlowCircuit(pc, lc, num_examples(data), Float64);
    for n in fc
        exp_pass_up_node(n, data)
    end
end

function exp_pass_up_node(node::ExpFlowNode{E}, data) where E
    pType = typeof(node.p_origin)
    fType = typeof(node.f_origin)

    if node.p_origin isa PlainSumNode && node.f_origin isa Logistic⋁Node
        # pthetas = [exp(node.p_origin.log_probs[i])
        #             for i in 1:length(children(node.p_origin)) for j in 1:length(children(node.f_origin))]
        # fthetas = [node.f_origin.thetas[j,:]
        #     for i in 1:length(node.p_origin.children) for j in 1:length(node.f_origin.children)]

        # node.f .= 0.0
        # node.fg .= 0.0
        # for z = 1:length(children(node))
        #     node.f  .+= pthetas[z] .* children(node)[z].f
        #     node.fg .+= (pthetas[z] .* fthetas[z]) .* children(node)[z].f
        #     node.fg .+= pthetas[z] .* children(node)[z].fg
        # end

        node.f .= 0.0
        node.fg .= 0.0
        N = length(children(node.p_origin))
        M = length(children(node.f_origin))
        @inline pthetas(i) = exp(node.p_origin.log_probs[i])
        @inline fthetas(j) = node.f_origin.thetas[j,:]
        @fastmath @inbounds @simd for i=1:N
            @simd for j=1:M
                z = (i-1)*M + j
                node.f  .+= pthetas(i) .* children(node)[z].f
                node.fg .+= (pthetas(i) .* fthetas(j)) .* children(node)[z].f
                node.fg .+= pthetas(i) .* children(node)[z].fg
            end
        end

    elseif node.p_origin isa PlainMulNode && node.f_origin isa Logistic⋀Node

        @fastmath @inbounds  begin
            # assume 2 children
            node.f .= children(node)[1].f .* children(node)[2].f
            node.fg .= (children(node)[1].f .* children(node)[2].fg) .+
                       (children(node)[2].f .* children(node)[1].fg)
        end

    elseif node.p_origin isa PlainProbLiteralNode
        if node.f_origin isa Logistic⋁Node
            m = children(node.f_origin)[1]
        elseif node.f_origin isa LogisticLiteralNode
            m = node.f_origin
        else
            error("Invalid Types of pairs {$pType} - {$fType}")
        end

        var = variable(m)
        X = data
        if ispositive(node.p_origin) && ispositive(m)
            node.f[:, .!isequal.(X[:, var], 0)] .= 1.0 # positive and missing observations
            node.f[:, isequal.(X[:, var], 0)] .= 0.0
        elseif isnegative(node.p_origin) && isnegative(m)
            node.f[:, .!isequal.(X[:, var], 1)] .= 1.0 # negative and missing observations
            node.f[:, isequal.(X[:, var], 1)] .= 0.0
        else
            node.f .= 0.0
        end

        if node.f_origin isa Logistic⋁Node
            node.fg .= node.f .* transpose(node.f_origin.thetas)
        else
            node.fg .= 0.0
        end

    else
        error("Invalid Types of pairs {$pType} - {$fType}")
    end

end

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
using Images, ImageView
using Video4Linux

#create a display function that consumes the video on channel vidchan
function displayliveYonly!(im1, vidchan)

    im1[:] = reinterpret(Gray{N0f8},take!(vidchan))
    ImageView.imshow!(canvas["gui"]["canvas"], im1, canvas["annotations"])

end

##
#Creaete Y component only 4:2:2 decoder type
# yonly = Video4Linux.UYVYonlyY(640,480)
yonly = Video4Linux.YUYVonlyY(640,480)
# create a videoproducer type channel (will only run 100 frames)
vidchan = Channel((c::Channel) -> videoproducer(c, yonly, devicename = "/dev/video0", N = 0; iomethod = Video4Linux.IO_METHOD_MMAP))

##
#capture one frame to create im1 and canvas needed to diplay
im1 = reinterpret(Gray{N0f8},take!(vidchan))
canvas = imshow(im1)

##

# run until channel is closed (With N = 0 it will run untul stopVideoProducer is called)
@async while isopen(vidchan)
    displayliveYonly!(im1, vidchan)
end

## run stopVideoProducer() to stop the infine loop and close the video channel and device
stopVideoProducer()

################################################################################
#full color 4:2:2

function displaylive!(im1, vidchan)

    A =  take!(vidchan)
    im1 = RGB.(YCbCr.(A[:,:,1], A[:,:,2], A[:,:,3]))
    ImageView.imshow!(canvas["gui"]["canvas"], im1, canvas["annotations"])
end

# ycrcb = Video4Linux.UYVY(640,480)
ycrcb = Video4Linux.YUYV(640,480)

vidchan = Channel((c::Channel) -> videoproducer(c, ycrcb, devicename = "/dev/video0", iomethod = Video4Linux.IO_METHOD_MMAP))

##
#capture one frame to create im1 and canvas needed to diplay
A =  take!(vidchan)
im1 = RGB.(YCbCr.(A[:,:,1], A[:,:,2], A[:,:,3]))
canvas = imshow(im1)

##

# run untul channel is closed (100 frames of vidchan)
while isopen(vidchan)
    @time displaylive!(im1, vidchan)
end

################################################################################
#10bit packed

function displaylive10bit!(im1, vidchan)
    im1 = take!(vidchan)
    ImageView.imshow!(canvas["gui"]["canvas"], im1, canvas["annotations"])
end

y10b = Video4Linux.Y10B(640,480)
vidchan = Channel((c::Channel) -> videoproducer(c, y10b), devicename = "/dev/video0")

##
#capture one frame to create im1 and canvas needed to diplay
im1 = take!(vidchan)
canvas = imshow(im1)

##

# run untul channel is closed (100 frames of vidchan)
while isopen(vidchan)
    @time displaylive10bit!(im1, vidchan)
end


text:
export PermMat
export randpermmat

import Base: kron

## object and constructors  ##

immutable PermMat{T<:Real} <: AbstractPerm{T}
    data::Vector{T}
end

length(p::PermMat) = (l=length(p.data); l * l)

# Single dimensional index means different things in different
# contexts. For matrix, it acts like a matrix.
function getindex(m::PermMat, k::Integer)
    i,j = divrem(k-1,plength(m))
    return m.data[j+1] == i+1 ? one(eltype(m)) : zero(eltype(m))
end

# should not be calling this, anyway
setindex!(m::PermMat, i::Int, k::Integer) = m.data[k] = i

# This is needed for eigvals()
similar(p::PermMat, T, dims...) = Array(T,dims...)
# Allows linear indexing, eg a[1:3]
similar(a::PermMat, dims::Dims) = Array(eltype(a), dims)

## operators ##

*{T}(p::PermMat, a::AbstractVector{T}) = PermPlain.permapply(p.data,a)
/(m1::PermMat, m2::PermMat) = PermMat(PermPlain.permcompose(m2.data,invperm(m1.data)))
*(m::PermMat, k::Integer) = full(m) * k

function +(p::PermMat, m::AbstractArray)
    size(p) == size(m) || error("dimensions must match")
    mout = copy(m)
    d = p.data
    mo = one(eltype(m))
    for i in 1:length(d)
        mout[i,d[i]] += mo
    end
    mout
end

function +(m::PermMat, p::PermMat)
    size(p) == size(m) || error("dimensions must match")
    dp = p.data
    dm = m.data
    mo = one(eltype(m))
    mout = full(m)
    for i in 1:length(dp)
        mout[i,dp[i]] += mo
    end
    mout
end

function +(m::AbstractArray, p::PermMat)
    size(p) == size(m) || error("dimensions must match")
    mout = copy(m)
    d = p.data
    mo = one(eltype(m))
    for i in 1:length(d)
        mout[i,d[i]] += mo
    end
    mout
end

function -(p::PermMat, m::AbstractArray)
    size(p) == size(m) || error("dimensions must match")
    mout = copy(m)
    d = p.data
    mo = one(eltype(m))
    for i in 1:length(d)
        j = d[i]
        mout[i,j] = -mout[i,j] + mo
    end
    mout
end

# # This a bit more efficient than full matrix multiplication. Cuts out one loop.
# function kron{T,S}(a::PermMat{T}, b::AbstractMatrix{S})
#     (nrowa, ncola) = size(a)
#     (nrowb, ncolb) = size(a)
#     R = zeros(promote_type(T,S), nrowa * nrowb, ncola * ncolb)
#     d = inv(a).data  # prbly better to not break abstraction everywhere
#     for j = 1:ncola, l = 1:ncolb
#         soff = ncola * (j-1)
#         i = d[j]
#         roff = ncola * (i-1)
#         for k = 1:nrowb
#             R[roff+k,soff+l] = b[k,l]
#         end
#     end
#     R
# end

# # For testing. If there is a difference, it is small
# function kron2{T,S}(a::PermMat{T}, b::AbstractMatrix{S})
#     (nrowa, ncola) = size(a)
#     (nrowb, ncolb) = size(a)
#     R = zeros(promote_type(T,S), nrowa * nrowb, ncola * ncolb)
#     d = a.data
#     for j = 1:ncola, l = 1:ncolb
#         soff = ncola * (j-1)
#         i = d[j]
#         roff = ncola * (i-1)
#         for k = 1:nrowb
#             R[soff+l,roff+k] = b[l,k]
#         end
#     end
#     R
# end

include("matlistcommon.jl")


text:
# # Beam under on/off loading: transient response

# Source code: [`beam_load_on_off_tut.jl`](beam_load_on_off_tut.jl)

# ## Description

# A cantilever beam is loaded by a trapezoidal-pulse traction load at its free
# cross-section. The load is applied within 0.015 seconds and taken off after
# 0.37 seconds. The beam oscillates about its equilibrium configuration.

# The beam is modeled as a solid. Trapezoidal rule is used to integrate the
# equations of motion in time. Rayleigh mass- and stiffness-proportional
# damping is incorporated. The dynamic stiffness is factorized for efficiency.

# ![](beam_load_on_off_tut.png)

# ## Goals

# - Show how to create the discrete model, with implicit dynamics and proportional damping.
# - Apply distributed loading varying in time.
# - Demonstrate  trapezoidal-rule time stepping.

##
# ## Definitions

# Basic imports.
using LinearAlgebra
using Arpack

# This is the finite element toolkit itself.
using FinEtools

# The linear stress analysis application is implemented in this package.
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule

# Input parameters
E = 205000*phun("MPa");# Young's modulus
nu = 0.3;# Poisson ratio
rho = 7850*phun("KG*M^-3");# mass density
loss_tangent = 0.005;
L = 200*phun("mm");
W = 4*phun("mm");
H = 8*phun("mm");
tolerance = W/500;
qmagn = 0.1*phun("MPa");
tend = 0.5*phun("SEC");

##
# ## Create the discrete model

MR = DeforModelRed3D
fens,fes  = H8block(L, W, H, 50, 2, 4)

geom = NodalField(fens.xyz)
u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

nl = selectnode(fens, box=[L L -Inf Inf -Inf Inf], inflate=tolerance)
setebc!(u, nl, true, 1)
setebc!(u, nl, true, 2)
setebc!(u, nl, true, 3)
applyebc!(u)
numberdofs!(u)

corner = selectnode(fens, nearestto=[0 0 0])
cornerzdof = u.dofnums[corner[1], 3]

material = MatDeforElastIso(MR, rho, E, nu, 0.0)

femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3,2)), material)
femm = associategeometry!(femm, geom)
K = stiffness(femm, geom, u)
femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
M = mass(femm, geom, u)

# Find the boundary finite elements at the tip cross-section of the beam. The
# uniform distributed loading will be applied to these elements.
bdryfes = meshboundary(fes)
# Those facing in the positive X direction will be chosen:
tipbfl = selectelem(fens, bdryfes, facing=true, direction=[-1.0 0.0 0.0])
# A base finite element model machine will be created to evaluate the loading.
# The force intensity is created as driven by a function, but the function
# really only just fills the buffer with the constant loading vector.
function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
    forceout .=  [0.0, 0.0, qmagn]
    return forceout
end
fi = ForceIntensity(FFlt, 3, pfun);

# The loading vector is lumped from the distributed uniform loading by
# integrating on the boundary. Hence, the dimension of the integration domain
# is 2.
el1femm =  FEMMBase(IntegDomain(subset(bdryfes,tipbfl), GaussRule(2,2)))
F = distribloads(el1femm, geom, u, fi, 2);

# The loading function is defined as a time -dependent multiplier of the
# constant distribution of the loading on the structure.
function tmult(t)
    if (t <= 0.015)
        t/0.015
    else
        if (t >= 0.4)
            0.0
        else
            if (t <= 0.385)
                1.0
            else
                (t - 0.4)/(0.385 - 0.4)
            end
        end
    end
end

##
# ## Time step determination

# We figure out the fundamental mode frequency, which will determine the time
# step is a fraction of the period.
evals, evecs = eigs(K, M; nev=1, which=:SM);

# The fundamental angular frequency is then:
@show omega_f = sqrt(evals[1]);

# We take  the time step to be a fraction of the period of vibration  in the
# fundamental mode.
@show dt = 0.05 * 1/(omega_f/2/pi);

##
# ## Damping model

# We take the damping to be representative of what's happening at the
# fundamental vibration frequency.

# For a given loss factor at a certain  frequency $\omega_f$, the
# stiffness-proportional damping coefficient may be estimated as
# 2*loss_tangent/$\omega_f$, and the mass-proportional damping coefficient may be
# estimated as 2*loss_tangent*$\omega_f$.
Rayleigh_mass = (loss_tangent/2)*omega_f;
Rayleigh_stiffness = (loss_tangent/2)/omega_f;

# Now we construct the Rayleigh damping matrix as a linear combination of the
# stiffness and mass matrices.
C = Rayleigh_stiffness * K + Rayleigh_mass * M

# The time stepping loop is protected by `let end` to avoid unpleasant surprises
# with variables getting clobbered by globals.
ts, corneruzs = let dt = dt, F = F
    # Initial displacement, velocity, and acceleration.
    U0 = gathersysvec(u)
    v = deepcopy(u)
    V0 = gathersysvec(v)
    U1 = fill(0.0, length(V0))
    V1 = fill(0.0, length(V0))
    F0 = deepcopy(F)
    F1 = fill(0.0, length(V0))
    R  = fill(0.0, length(V0))

    # Factorize the dynamic stiffness
    DSF = cholesky((M + (dt/2)*C + ((dt/2)^2)*K))

    # The times and displacements of the corner will be collected into two vectors
    ts = Float64[]
    corneruzs = Float64[]
    # Let us begin the time integration loop:
    t = 0.0;
    step = 0;
    F0 .= tmult(t) .* F
    while t < tend
        push!(ts, t)
        push!(corneruzs, U0[cornerzdof])
        t = t+dt;
        step = step + 1;
        (mod(step,100)==0) && println("Step $(step): $(t)")
        # Set the time-dependent load
        F1 .= tmult(t) .* F
        # Compute the out of balance force.
        R = (M*V0 - C*(dt/2*V0) - K*((dt/2)^2*V0 + dt*U0) + (dt/2)*(F0+F1));
        # Calculate the new velocities.
        V1 = DSF\R;
        # Update the velocities.
        U1 = U0 + (dt/2)*(V0+V1);
        # Switch the temporary vectors for the next step.
        U0, U1 = U1, U0;
        V0, V1 = V1, V0;
        F0, F1 = F1, F0;
        if (t == tend) # Are we done yet?
            break;
        end
        if (t+dt > tend) # Adjust the last time step so that we exactly reach tend
            dt = tend-t;
        end
    end
    ts, corneruzs # return the collected results
end

##
# ## Plot the results

using Gnuplot

@gp "set terminal windows 0 " :-
@gp  :- ts corneruzs./phun("mm") "lw 2 lc rgb 'red' with lines title 'Displacement of the corner' "
@gp  :- "set xlabel 'Time [s]'"
@gp  :- "set ylabel 'Displacement [mm]'"

# The end.
true


text:
module LMCLUS

using LinearAlgebra
using Logging
import StatsBase: fit, Histogram, sturges
import Statistics: mean
import MultivariateStats: PCA, fit, indim, outdim, projection
import Clustering: ClusteringResult, assignments, counts, nclusters
import Distributed: Future, remotecall, nprocs, myid
import Random: _randjump, MersenneTwister, randperm, DSFMT, GLOBAL_RNG, AbstractRNG

export  lmclus,

        separation,
        distance_to_manifold,

        Manifold,
        indim,
        outdim,
        points,
        threshold,
        projection,
        criteria,

        nclusters,
        counts,
        assignments,
        manifold,
        manifolds,
        LMCLUSResult,
        refine

include("types.jl")
include("params.jl")
include("results.jl")
include("utils.jl")
include("separation.jl")
include("mdl.jl")

const TRACE = LogLevel(Base.CoreLogging.Debug.level-1000)
const DEBUG_SAMPLING = LogLevel(Base.CoreLogging.Debug.level-900)
const DEBUG_SEPARATION = LogLevel(Base.CoreLogging.Debug.level-800)

#
# Linear Manifold Clustering
#
function lmclus(X::AbstractMatrix{T}, params::Parameters, np::Int=nprocs()) where T<:Real
    # Setup RNG
    seed = getseed(params)
    mts = if np == 1
        [MersenneTwister(seed)]
    else
        #randjump(MersenneTwister(seed), np)
        mt = MersenneTwister(seed)
        jumps = fill(2^30, np-1)
        rngs = accumulate((mt, j)->_randjump(mt, DSFMT.calc_jump(j)), jumps, init=mt)
        pushfirst!(rngs, mt)
    end
    return LMCLUSResult(lmclus(X, params, mts)...)
end

function lmclus(X::AbstractMatrix{T}, params::Parameters, prngs::Vector{MersenneTwister}) where T<:Real

    @assert length(prngs) >= nprocs() "Number of PRNGS cannot be less then processes."

    N, n = size(X)
    index = collect(1:n)
    number_of_clusters = 0
    manifolds = Manifold[]
    separations = Separation[]

    # Check if manifold maximum dimension is less then full dimension
    if N <= params.max_dim
        params.max_dim = N - 1
        @debug "Adjusting maximum manifold dimension to $(params.max_dim)"
    end

    # Main loop through dataset
    while length(index) > params.min_cluster_size

        # Find one manifold
        best_manifold, best_separation, remains = find_manifold(X, index, params, prngs, length(manifolds))

        # Add a new manifold cluster to collection
        push!(manifolds, best_manifold)
        push!(separations, best_separation)

        number_of_clusters += 1
        @debug "Found cluster" number_of_clusters size=size(best_manifold) dimension=outdim(best_manifold)

        # Stop clustering if found specified number of clusters
        length(manifolds) == params.stop_after_cluster && break

        # Look at the rest of the dataset
        index = remains
    end

    # Rest of the points considered as noise
    if length(index) > 0
        @debug "Outliers" number=length(index)
        outliers = Manifold(0, zeros(T, N), zeros(T, N, 0), index)
        if params.basis_alignment
            adjustbasis!(outliers, X, adjust_dim=params.dim_adjustment, adjust_dim_ratio=params.dim_adjustment_ratio)
        end
        push!(manifolds, outliers)
        push!(separations, Separation())
    end

    return manifolds, separations
end

# Find manifold in multiple dimensions
# The algorithm contains following state machine:
# 0. Set inspected subspace dimension to 1.
# 1. Find linear subspace (translation & basis) with best possible separation
# 2. Evaluate the separation criteria
# 3. If separation criterion is met then
# 3.1. Filter out points within separation threshold
# 3.2. Form cluster from the points and remove them from the dataset, go to 5.
# 4. If separation criterion not is met then
# 4.1. If the option is set and the manifold basis wasn't adjusted before then
#      adjust the linear manifold basis of the cluster subspace and go to 2.
# 4.2. If the option is set and the separation search within the manifold subspace wasn't performed then
#      search for the separation within the manifold subspace and go to 2.
# 5. If the number of points in the cluster is less then manimal acceptable size
#    then go to 9.
# 6. If the option is set and the separation was found then
#    perform MDL evaluation of the former cluster.
# 6.1. If compression ration is not met, return cluster points to the dataset
#      then go to 8.
# 7. Return the formed cluster.
# 8. Increase the dimension of the inspected subspace
# 8.1. If the subspace dimension equals to N-1, go to 9.
# 8.2. Otherwise, go to 1.
# 9. Form cluster from the rest of the dataset points.
function find_manifold(X::AbstractMatrix{T}, index::Vector{Int},
                       params::Parameters,
                       prngs::Vector{MersenneTwister},
                       found::Int=0) where T<:AbstractFloat
    filtered = Int[]
    selected = copy(index)
    N = size(X,1) # full space dimension
    best_manifold = Manifold{T}(params.min_dim, index)
    best_separation = Separation()

    sep_dim = params.min_dim
    while sep_dim <= params.max_dim
        separations = 0
        state = :SEPARATION

        # get thresholds
        θ, σ = T(Inf), T(Inf)

        # search appropriate linear manifold subspace for best distance separation
        while true
            @debug "Algorithm State" state best_manifold
            if state == :SEPARATION
                # get separation manifold
                @debug("manifold: find suitable subspace...",
                    data_size=length(selected),
                    linear_manifold_dimension=sep_dim,
                    space_dim=N
                )
                sep, origin, basis = find_separation_basis(X[:, selected], sep_dim, params, prngs, found)
            elseif state == :ALIGNMENT && params.basis_alignment && size(best_manifold) > 0
                # check if the adjusted basis provides better separation
                previous_outdim = outdim(best_manifold)
                adjustbasis!(best_manifold, X, adjust_dim=params.dim_adjustment, adjust_dim_ratio=params.dim_adjustment_ratio)
                origin, basis = mean(best_manifold), projection(best_manifold)
                @debug "manifold: perform basis adjustment..." outdim=outdim(best_manifold) previous_outdim origin basis
                idxs = points(best_manifold)
                sep = find_separation(view(X, :, idxs), origin, basis, params, prngs[1])
            elseif state == :BOUND && params.bounded_cluster && size(best_manifold) > 0
                # check if the bounded cluster provides better separation
                origin, basis = mean(best_manifold), projection(best_manifold)
                @debug "manifold: separating within manifold subspace..." origin basis
                idxs = points(best_manifold)
                sep = find_separation(view(X, :, idxs), origin, basis, params, prngs[1], ocss = true)
            else
                break
            end
            @debug("Separation",
                width=sep.discriminability,
                depth=sep.depth,
                threshold=threshold(sep),
                extrema=extrema(sep),
                global_minimum=sep.globalmin,
                total_bins=sep.bins,
                criteria=criteria(sep),
                best_bound=params.best_bound
            )

            # No good separation found
            if criteria(sep) < params.best_bound
                currthr  = state == :BOUND ? σ : θ
                thr = extrema(sep)[2]
                if currthr > 0.0 && thr > 0.0
                    currthr = min(currthr, thr)
                end
                if state == :BOUND
                    σ = best_manifold.σ = currthr
                else
                    θ = best_manifold.θ = currthr
                end

                # swith to next state
                state = (state == :SEPARATION) ? :ALIGNMENT :
                        (state == :ALIGNMENT) ?  :BOUND     :
                        (state == :BOUND) ?      :FINISHED  : :NONE
            else
                thr = threshold(sep)
                θ, σ  = state == :BOUND ? (θ, thr) : (thr, σ)

                # filter separated points
                separated, removed = filter_separated(selected, X, origin, basis, sep,
                                                    ocss=(state == :BOUND))

                # small amount of points is considered noise, try to bump dimension
                if length(separated) <= params.min_cluster_size
                    @debug "Separated points cannot form cluster. Its size is too small ($(length(separated)))."
                    state = :FINISHED
                else
                    # create manifold cluster from good separation found
                    best_manifold = Manifold(sep_dim, origin, basis, separated, θ, σ)
                    best_separation = sep

                    # partition cluster from the dataset
                    append!(filtered, removed)

                    # refine dataset
                    selected = separated
                    @debug "Separated points" size=length(selected)
                    state = :SEPARATION
                    separations += 1
                end
            end
            # @debug("Best manifold",
            #     dim=outdim(best_manifold),
            #     origin=mean(best_manifold),
            #     basis=projection(best_manifold),
            #     threshold=threshold(best_manifold),
            # )
        end

        if length(selected) <= params.min_cluster_size # no more points left - finish clustering
            @debug "Noise" minimum_cluster_size=params.min_cluster_size current_cluster_size=length(selected)
            break
        end

        # check compression ratio
        if params.mdl && outdim(best_manifold) > 0 && separations > 0
            BM = best_manifold
            BMdata = X[:, selected]
            Pm = params.mdl_model_precision
            Pd = params.mdl_data_precision
            adjustbasis!(BM, X) # generate full basis - needed to calculate MDL
            mmdl = MDL.calculate(params.mdl_algo, BM, BMdata, Pm, Pd, ɛ = T(params.mdl_quant_error))
            mraw = MDL.calculate(MDL.Raw, BM, BMdata, Pm, Pd)

            cratio = mraw/mmdl
            @debug("Minimal Description Length",
                mdl_value=mmdl,
                raw_value=mraw,
                compression=cratio,
                compression_threshold=params.mdl_compres_ratio,
                action=(cratio < params.mdl_compres_ratio ? :reject : :accept)
            )
            if cratio < params.mdl_compres_ratio
                # reset dataset to original state
                append!(selected, filtered)
                filtered = Int[]
                separations = 0
            end
        end

        !params.force_max_dim && separations > 0 && break

        sep_dim += 1
        @debug "Separation not found" action=(sep_dim == params.max_dim ? "forming cluster" : "increasing dimension") dim=sep_dim
    end

    # Cannot find any manifold in data then form last cluster
    require_alignment = false
    if outdim(best_manifold) == N
        if size(best_manifold) == 0
            best_manifold.points = selected
        end
        @debug "Forming a cluster from the rest of $(size(best_manifold)) points"
        best_manifold.d = 1

        # calculate bounds
        best_manifold.θ = maximum(distance_to_manifold(view(X, :,selected), best_manifold))
        best_manifold.σ = if params.bounded_cluster
            maximum(distance_to_manifold(view(X, :,selected), best_manifold, ocss=true))
        else
            Inf
        end
        require_alignment = true
    end

    # Adjust cluster basis & dimension if it is the last clusters or a corresponding setting is set
    if params.basis_alignment || require_alignment
        adjustbasis!(best_manifold, X, adjust_dim=params.dim_adjustment, adjust_dim_ratio=params.dim_adjustment_ratio)
    end

    best_manifold, best_separation, filtered
end

# LMCLUS main function:
# 1- sample trial linear manifolds by sampling points from the data
# 2- create distance histograms of the data points to each trial linear manifold
# 3- of all the linear manifolds sampled select the one whose associated distance histogram
#    shows the best separation between to modes.
function find_separation_basis(X::AbstractMatrix{T}, lm_dim::Int,
                               params::Parameters,
                               prngs::Vector{MersenneTwister},
                               found::Int=0) where {T<:Real}
    full_space_dim, data_size = size(X)

    # determine number of samples of lm_dim+1 points
    Q = sample_quantity( lm_dim, full_space_dim, data_size, params, found)

    # divide samples between PRNGs
    samples_proc = round(Int, Q / length(prngs))

    # enable parallel sampling
    bases = Array{Future}(undef, length(prngs))
    np = nprocs()
    nodeid = myid()
    for i in 1:length(prngs)
        pid = nodeid == 1 ? (i%np)+1 : nodeid #  if running from node 1
        bases[i] = remotecall(sample_manifold, pid, X, lm_dim+1, params, prngs[i], samples_proc)
    end

    # reduce values of manifolds from remote sources
    best_sep, best_origin, best_basis = Separation(), zeros(T, 0), zeros(T, 0, 0)
    for (i, rr) in enumerate(bases)
        res = fetch(rr)
        # isa(res, Exception) && throw(res)
        sep, origin, basis, mt = res
        if criteria(sep) > criteria(best_sep) # get largest separation
            best_sep = sep
            best_origin = origin
            best_basis = basis
        end
        prngs[i] = mt # save state of RNG
    end

    # bad sampling
    criteria(best_sep) <= 0. &&  @debug "no good histograms to separate data !!!"
    return best_sep, best_origin, best_basis
end

function sample_manifold(X::Matrix{T}, lm_dim::Int,
                         params::Parameters, prng::MersenneTwister,
                         num_samples::Int) where {T<:Real}
    best_sep = Separation()
    best_origin = zeros(T, 0)
    best_basis = zeros(T, 0, 0)

    for i in 1:num_samples
        sample = sample_points(X, lm_dim, rng=prng)
        length(sample) == 0 && continue

        origin, basis = form_basis(X, sample)
        sep = find_separation(X, origin, basis, params, prng)
        if criteria(sep) > criteria(best_sep)
            best_sep = sep
            best_origin = origin
            best_basis = basis
            @logmsg TRACE "Found new separation" separation=best_sep origin=best_origin basis=best_basis
        end
    end

    return best_sep, best_origin, best_basis, prng
end

"""Find separation criteria

Given a dataset `X`, an tanslation vector `origin` and set of basis vectors `basis` of the linear manifold.

`ocss` parameter enables separation for the orthogonal complement subspace of the linera manifold
"""
function find_separation(X::AbstractMatrix, origin::AbstractVector,
                         basis::AbstractMatrix, params::Parameters,
                         prng::MersenneTwister;
                         ocss::Bool = false)
    # Define sample for distance calculation
    if params.histogram_sampling
        Z_01=2.576  # Z random variable, confidence interval 0.99
        delta_p=0.2
        delta_mu=0.1
        P=1.0/params.number_of_clusters
        Q=1-P
        n1=(Q/P)*((Z_01*Z_01)/(delta_p*delta_p))
        p=( P<=Q ? P : Q )
        n2=(Z_01*Z_01)/(delta_mu*delta_mu*p)
        n3= ( n1 >= n2 ? n1 : n2 )
        n4= round(Int, n3)
        n= ( size(X, 2) <= n4 ? size(X, 2)-1 : n4 )

        @logmsg TRACE "sample distances for histogram" samples=n
        sampleIndex = sample_points(X, n, rng=prng)
    end

    # calculate distances to the manifold
    distances = distance_to_manifold(params.histogram_sampling ? view(X, :,sampleIndex) : X ,
                                     origin, basis, ocss = ocss)
    # define histogram size
    bins = gethistogrambins(distances, params.max_bin_portion, params.hist_bin_size, params.min_bin_num)
    # perform separation
    return separation(params.sep_algo, distances, bins=bins)
end

# Determine the number of times to sample the data in order to guaranty
# that the points sampled are from the same cluster with probability
# of error that does not exceed an error bound.
# Three different types of heuristics may be used depending on LMCLUS's input parameters.
function sample_quantity(k::Int, full_space_dim::Int, data_size::Int,
                         params::Parameters, S_found::Int)

    S_max = params.number_of_clusters
    if S_max <= 1
        return 1 # case where there is only one cluster
    end

    #p = 1.0 / S_max    # p = probability that 1 point comes from a certain cluster
    p = 1.0 / max(2, S_max-S_found)
    P = p^k            # P = probability that "k+1" points are from the same cluster
    N = abs(log10(params.error_bound)/log10(1-P))
    num_samples = 0

    if params.sampling_heuristic == 1
        num_samples = isinf(N) ? typemax(Int) : round(Int, N)
    elseif params.sampling_heuristic == 2
        NN = data_size*params.sampling_factor
        num_samples = isinf(NN) ? typemax(Int) : round(Int, NN)
    elseif params.sampling_heuristic == 3
        NN = min(N, data_size*params.sampling_factor)
        num_samples = isinf(NN) ? typemax(Int) : round(Int, NN)
    end
    num_samples = num_samples > 1 ? num_samples : 1

    @debug "number of samples" first_heuristic=N second_heuristic=data_size*params.sampling_factor

    return num_samples
end

"""Calculate histogram size

Given a `x` value vector, the number of bins for the histogram is `hist_bin_size`.
If `max_bin_portion` value is in closed unit range, then the number of bins estimated from the bin size determined as `x_n+i - x_i`,
i.e. find the smallest difference between i successive points, where `x_n` is a point with index `n`,
and `i` is the max number of points we allow to be stored in a single bin.

*Note: the number of bins cannot be less then `min_bin_num`*
"""
function gethistogrambins(x::Vector, max_bin_portion::Float64, hist_bin_size::Int, min_bin_num::Int)
    bns = hist_bin_size
    if max_bin_portion > 0.0
        l = length(x)
        sort!(x)
        xmin, xmax = x[1], x[end]
        xrng = xmax - xmin
        mbp = round(Int, l * max_bin_portion)

        binwidth = xmax
        for i in mbp:mbp:l
            diff = x[i] - xmin
            if binwidth > diff && diff > 0.0
                binwidth = diff
            end
            xmin = x[i]
        end
        bns = round(Int, xrng/binwidth)
        bns = bns > l ? sturges(l) : bns # deal with special cases
    end
    return max(min_bin_num, bns) # special cases: use max bin size
end

"""Calculate the distances from points, `X`, to a linear manifold with `basis` vectors, translated to `origin` point.

`ocss` parameter turns on the distance calculatation to orthogonal complement subspace of the given manifold.
"""
function distance_to_manifold(X::AbstractMatrix{T}, origin::AbstractVector{T}, basis::AbstractMatrix{T};
                              ocss::Bool = false) where T<:Real
    N, n = size(X)
    M = size(basis,2)
    # vector to hold distances of points from basis
    distances = zeros(T, n)
    tran = similar(X)
    @fastmath @inbounds for i in 1:n
        @simd for j in 1:N
            tran[j,i] = X[j,i] - origin[j]
            if !ocss
                distances[i] += tran[j,i]*tran[j,i]
            end
        end
    end
    proj = transpose(basis) * tran
    @fastmath @inbounds for i in 1:n
        b = 0.0
        @simd for j in 1:M
            b += proj[j,i]*proj[j,i]
        end
        distances[i] = sqrt(abs(distances[i]-b))
    end
    return distances
end
distance_to_manifold(X::AbstractMatrix{<:Real}, M::Manifold; ocss::Bool = false) =
    distance_to_manifold(X, mean(M), projection(M), ocss=ocss)

"""Calculate the distance from the `point` to the manifold, described by a `basis` matrix, and its orthoganal complement."""
function distance_to_manifold(point::AbstractVector{T}, basis::AbstractMatrix{T}) where T<:Real
    dpnt = sum(abs2, point)
    dprj = sum(abs2, transpose(basis) * point)
    return sqrt(abs(dpnt-dprj)), sqrt(dprj)
end
distance_to_manifold(point::AbstractVector, origin::AbstractVector, basis::AbstractMatrix) = distance_to_manifold(point - origin, basis)
distance_to_manifold(point::AbstractVector, M::Manifold) = distance_to_manifold(point - mean(M), projection(M))

end # module


text:
# # Faster convergence for linear dynamcs

#=
In this example, I show how the SSA can be used as objective to guarantee
convergence or faster convergence for a linear dynamical system.
=#

# ## Optimization of all matrix elements

# ### Initialization
using Plots,NamedColors
using LinearAlgebra
using Random
using SmoothedSpectralAbscissa ; const SSA=SmoothedSpectralAbscissa
Random.seed!(0);

# ### Linear Dynamics

#=
Consider continuous linear dynamics, regulated by
```math
\frac{\text{d} \mathbf{x}}{\text{d}t} = A \mathbf{x}
```
The soluton is analytic and takes the form:
```math
\mathbf{x}(t) = \exp( A\,t )\;\mathbf{x}_0
```
Where $\mathbf{x}_0$ are the initial conditions. The function below computes the dynamical
evolution of the system.
=#

function run_linear_dyn(A::Matrix{R},x0::Vector{R},tmax::Real,dt::Real=0.01) where R
  ts = range(0,tmax;step=dt)
  ret = Matrix{R}(undef,length(x0),length(ts))
  for (k,t) in enumerate(ts)
    ret[:,k]=exp(A.*t)*x0
  end
  retnrm = mapslices(norm,ret;dims=1)[:]
  return ts,ret,retnrm
end;

# ### The optimization of the objective is done through Optim.jl and BFGS
using Optim

function objective_and_grad_simple(x::Vector{R},grad::Union{Nothing,Vector{R}},
    n::Integer,ssa_eps::R,alloc::SSA.SSAAlloc) where R
  mat=reshape(x,(n,n))
  gradmat = isnothing(grad) ? nothing : similar(mat)
  obj = SSA.ssa!(mat,gradmat,alloc,ssa_eps)
  if !isnothing(grad)
    for i in eachindex(gradmat)
      grad[i]=gradmat[i]
    end
  end
  return obj
end;
## #src
# ### Start with an unstable matrix

n = 50
A = randn(n,n) ./ sqrt(n) + 0.2I
x0=randn(n)
times,_,dyn_norms = run_linear_dyn(A,x0,3.,0.1)

plot(times,dyn_norms; leg=false,linewidth=3,color=:black,xlabel="time",ylabel="norm(x(t))")

# as expected, the norm grows exponentially with time.

# ### Now do the gradient-based optimization
const ssa_eps=0.001
const alloc = SSA.SSAAlloc(n)
const y0 = A[:];

#=
The objective function to be minimized is
```math
\text{obj}(A) = \text{SSA}(A) + \lambda \frac12 \left\| A - A_0 \right\|^2
```
Where $\lambda$ sets the relative weight.
We are reducing the SSA while keeping the matrix elements close to their initial value.
Adding some form of regularization is **always** necessary when otpimizing. If not,
the SSA would run to $-\infty$.
=#

function objfun!(F,G,y)
  λ = 50.0/length(y) # regularizer weight
  obj=objective_and_grad_simple(y,G,n,ssa_eps,alloc) # SSA and gradient
  ydiffs = y.-y0
  obj += 0.5*λ*mapreduce(x->x^2,+,ydiffs) # add the regularizer
  if !isnothing(G)
    @. G += λ*ydiffs # add gradient of regularizer
  end
  return obj
end;

# ### Optimize and show the results

opt_out = optimize(Optim.only_fg!(objfun!),A[:],BFGS(),Optim.Options(iterations=50))
y_opt=Optim.minimizer(opt_out)
A_opt = reshape(y_opt,(n,n))
times,_,dyn_norms_opt = run_linear_dyn(A_opt,x0,3.,0.1)
plot(times,dyn_norms_opt;
  leg=false,linewidth=3,color=:blue,xlabel="time",ylabel="norm(x(t)) optimized")

#=
The optimized matrix produces stable dynamics, as shown in the plot above.

We can also take a look at the matrices before and after optimization
=#

heatmap(hcat(A,fill(NaN,n,10),A_opt);ratio=1,axis=nothing,ticks=nothing,border=:none,
  colorbar=nothing)

#=
The optimized version simply has negative diagonal terms.

This may appear a bit trivial. In the next part, I optimize a system *excluding* the diagonal.
=#

# ## Optimization that excludes the diagonal

#=
Here I consider a matrix that is stable, but produces a large nonlinear amplification.
It is generated by the function:
=#
function rand_nonnormal(n::Integer,(ud::Real)=1.01)
  mat = randn(n,n) ./ sqrt(n)
  @show SSA.spectral_abscissa(mat)
  mat = mat -  (1.1*SSA.spectral_abscissa(mat))*I
  sh = schur(mat)
  upd = diagm(0=>fill(1.0,n),1=>fill(ud,n-1))
  return sh.vectors*upd*sh.Schur*inv(upd)*sh.vectors'
end
# Let's make one and see how it looks like
A = rand_nonnormal(n,1.0)
x0=randn(n)
times,dyn_t,dyn_norms = run_linear_dyn(A,x0,30.,0.5)
plot(times,dyn_norms;
  leg=false,linewidth=3,color=:black,xlabel="time",ylabel="norm(x(t))")

#=
The norm is initally amplified, and decreases slowly. This is due to the
non-normality of matrix $A$.
=#

# ### Objective function that excludes diagonal

function objective_and_grad_nodiag(x::Vector{R},grad::Union{Nothing,Vector{R}},
    n::Integer,ssa_eps::R,alloc::SSA.SSAAlloc,A0::Matrix{R}) where R
  mat=reshape(x,(n,n))
  for i in 1:n
    mat[i,i]=A0[i,i] # diagonal copied from original matrix
  end
  gradmat = isnothing(grad) ? nothing : similar(mat)
  obj = SSA.ssa!(mat,gradmat,alloc,ssa_eps)
  if !isnothing(grad)
    for i in 1:n
      gradmat[i,i] = 0.0 # diagonal has zero gradient
    end
    for i in eachindex(gradmat)
      grad[i]=gradmat[i] # copy the gradient
    end
  end
  return obj
end;

# ### Optimizer and optimization

#=
The only difference from before is using
`objective_and_grad_nodiag` rather than `objective_and_grad_simple`
=#

const ssa_eps=0.001
const alloc = SSA.SSAAlloc(n)
const y0 = A[:];

function objfun!(F,G,y)
  λ = 1.0/length(y) # regularizer weight
  obj=objective_and_grad_nodiag(y,G,n,ssa_eps,alloc,A) # add the regularizer
  ydiffs = y.-y0
  obj += 0.5*λ*mapreduce(x->x^2,+,ydiffs)
  if !isnothing(G)
    @. G += λ*ydiffs # gradient of regularizer
  end
  return obj
end

opt_out = optimize(Optim.only_fg!(objfun!),A[:],BFGS(),Optim.Options(iterations=50))
y_opt=Optim.minimizer(opt_out)
A_opt = reshape(y_opt,(n,n));

# Now the optimized matrix looks remarkably similar to ro the original one,
# as shown below.
heatmap(hcat(A,fill(NaN,n,10),A_opt);ratio=1,axis=nothing,ticks=nothing,border=:none)

# Here I show the differences between $A$ and $A_{\text{opt}}$.
# (in case you don't believe they are different)
heatmap(A-A_opt;ratio=1,axis=nothing,ticks=nothing,border=:none)

# Let's compare the time evolution of the norms
times,dyn_t_opt,dyn_norms_opt = run_linear_dyn(A_opt,x0,30.,0.5)
plot(times,[dyn_norms dyn_norms_opt];
      leg=:topright,linewidth=3,color=[:black :blue],
      xlabel="time",ylabel="norm(x(t))", label=["before otpimization" "after optimization"])
#=
![So much stability!](./meme1.png)
=#

# ## Extras

#=
In gradient based optimization with no automatic differentiation, it is
always necessary to test the gradient of the objective function.
The procedure is illustrated below.
=#

using Calculus
function test_gradient(myobjfun,y0)
  grad_an = similar(y0)
  _ = myobjfun(1.0,grad_an,y0) # compute gradient analytically
  grad_num = Calculus.gradient(y->myobjfun(1.0,nothing,y),y0) # compute it numerically
  return (grad_an,grad_num)
end

_ = let do_the_test = false
  if do_the_test
    (x1,x2)=test_gradient(objfun!,randn(n^2))
    scatter(x1,x2;ratio=1)
    plot!(identity)
  end
end


text:
# read the output ilimit
export psp_readilimit

const READ_ILIMIT = [0x49, 0x0d]

"""
    psp_readilimit(io)

Read the output current limit.
"""
function psp_readilimit(io_psp::IO)
  write(io_psp, READ_ILIMIT)
  ilimitstring = bytestring(read(io_psp,UInt8, 8))
  if  uppercase(ilimitstring[1]) != 'I' ||
      ilimitstring[3] != '.' ||
      ilimitstring[6:8] != "\r\r\n"
        error("PSP Invalid Response: ",ilimitstring)
  end
  parse(Float64,ilimitstring[2:5])
end


text:
include("constants.jl")

"""
Parse input string s as a float, otherwise return original string.
"""
function parsestring(s)
    x = tryparse(Float64, s)
    if isnothing(x)
        x = s
    end
    return x
end

"""
Return non-comment (#) lines from input file parsed as 2D array
"""
function readfile(filename)
    lines = readlines(filename)
    dat = [parsestring.(split(strip(line))) for line in lines if line[1] ≠ '#']
    return hcat(dat...)'
end

function loadmethyls(filename)
    println("load methyls: $filename")
    return readlines(filename)
end

function loadR2(filename)
    println("load R2: $filename")

    tmp = readfile(filename)
    R2_zq = tmp[:,1:8:end] .± tmp[:,2:8:end]
    R2_dq = tmp[:,3:8:end] .± tmp[:,4:8:end]
    R2_dqq = tmp[:,5:8:end] .± tmp[:,6:8:end]
    R2_qq = tmp[:,7:8:end] .± tmp[:,8:8:end]

    return R2_zq, R2_dq, R2_dqq, R2_qq
end

function loadS2tc(filename)
    println("load S2tc: $filename")

    S2tc = readfile(filename)
    return S2tc[:,1] .± S2tc[:,2]
end

function load1Hcsa(ηHfilename, S2tc, B0_MHz)
    println("load and calculate 1H CSA:")
    println("   ηH filename: $ηHfilename")
    println("   B0 (1H, in MHz) = $B0_MHz")

    ηH = readfile(ηHfilename)
    ηH = ηH[:,1] .± ηH[:,2]

    B0 = B0_MHz * 1e6 * 2π / γH;

    ΔσH = @. 1e6*ηH/(-(4/15)*(μ₀/(4π))*ħ*γH^2*γC*P2cosβ*rCH^-3*B0*S2tc*1e-9)
    return ΔσH
end

function load13Ccsa(filename)
    println("load 13C csa: $filename")

    ΔσC = readfile(filename)
    return ΔσC[:,1] .± ΔσC[:,2]
end

max_stars_repo_name: chriswaudby/qq-hahnecho

max_issues_repo_name: chriswaudby/qq-hahnecho

max_forks_repo_name: chriswaudby/qq-hahnecho


text:
# a Node contains the minimum needed to hold the tree structure
struct Node{ACTIONC}
    action_ix::Int32 # action index leading to this state
    parent::Nullable{Node{ACTIONC}}
    children::Dict{Int32, Node{ACTIONC}}
    id::Int32
    function Node(
                  action_ix::Int32,
                  parent::Nullable{Node{ACTIONC}},
                  id::Int32) where ACTIONC
        children = Dict{Int32, Node{ACTIONC}}()
        return new{ACTIONC}(action_ix, parent, children, id)
    end
end

struct HeapCell
    is_used::Bool # because we cannot efficiently pop a random element from a heap
    score::Float32
    action_ix::Int32
    parent_id::Int32
end

mutable struct Heap
    cells::Vector{HeapCell}
    total_used::Int64
    Heap() = new(Vector{HeapCell}(), 0)
end

# The tree contains the nodes (describing the structure)
# and any additional data per node
struct Tree{STATE, SENSOR, ENV, ACTIONC}
    env::ENV
    heap::Heap
    # these 6 vectors are indexed by node.id (SOA style)
    nodes::Vector{Node{ACTIONC}}
    sensors::Vector{SENSOR}
    children_qs::Vector{Vector{Float32}}
    states::Vector{STATE}
    accumulated_rewards::Vector{Float32}
    dones::Vector{Bool}
end


text:
export PowerDivergenceTest, ChisqTest, MultinomialLRTest

const Levels{T} = Tuple{UnitRange{T},UnitRange{T}}

function boundproportion(x::T) where T<:Real
    max(min(convert(Float64,x),1.0),0.0)
end

struct PowerDivergenceTest <: HypothesisTest
    lambda::Float64
    theta0::Vector{Float64}
    stat::Float64
    df::Int64
    observed::Matrix{Int64}
    n::Int64
    thetahat::Vector{Float64}

    expected::Matrix{Float64}
    residuals::Matrix{Float64}
    stdresiduals::Matrix{Float64}
end

#test name
function testname(x::PowerDivergenceTest)
    if x.lambda == 1
        return "Pearson's Chi-square Test"
    elseif x.lambda == 0
        return "Multinomial Likelihood Ratio Test"
    elseif x.lambda == -1
        return "Minimum Discrimination Information Test"
    elseif x.lambda == -2
        return "Neyman Modified Chi-square Test"
    elseif x.lambda == -0.5
        return "Freeman-Tukey Test"
    else
        return "Power Divergence Test"
    end
end

# parameter of interest: name, value under h0, point estimate
population_param_of_interest(x::PowerDivergenceTest) = ("Multinomial Probabilities", x.theta0, x.thetahat)
default_tail(test::PowerDivergenceTest) = :right

pvalue(x::PowerDivergenceTest; tail=:right) = pvalue(Chisq(x.df),x.stat; tail=tail)

"""
    confint(test::PowerDivergenceTest, alpha = 0.05; tail = :both, method = :sison_glaz)

Compute a confidence interval with coverage 1-`alpha` for multinomial proportions using
one of the following methods. Possible values for `method` are:

  - `:sison_glaz` (default): Sison-Glaz intervals
  - `:bootstrap`: Bootstrap intervals
  - `:quesenberry_hurst`: Quesenberry-Hurst intervals
  - `:gold`: Gold intervals (asymptotic simultaneous intervals)

# References

  * Agresti, Alan. Categorical Data Analysis, 3rd Edition. Wiley, 2013.
  * Sison, C.P and Glaz, J. Simultaneous confidence intervals and sample size determination
    for multinomial proportions. Journal of the American Statistical Association,
    90:366-369, 1995.
  * Quesensberry, C.P. and Hurst, D.C. Large Sample Simultaneous Confidence Intervals for
    Multinational Proportions. Technometrics, 6:191-195, 1964.
  * Gold, R. Z. Tests Auxiliary to ``χ^2`` Tests in a Markov Chain. Annals of
    Mathematical Statistics, 30:56-74, 1963.
"""
function StatsBase.confint(x::PowerDivergenceTest, alpha::Float64=0.05;
                           tail::Symbol=:both, method::Symbol=:sison_glaz, correct::Bool=true,
                           bootstrap_iters::Int64=10000, GC::Bool=true)
    check_alpha(alpha)

    m  = length(x.thetahat)

    if tail == :left
        i = StatsBase.confint(x, alpha*2,method=method, GC=GC)
        Tuple{Float64,Float64}[(0.0, i[j][2]) for j in 1:m]
    elseif tail == :right
        i = StatsBase.confint(x, alpha*2,method=method, GC=GC)
        Tuple{Float64,Float64}[(i[j][1], 1.0) for j in 1:m]
    elseif tail == :both
        if method == :gold
            ci_gold(x,alpha,correct=correct,GC=GC)
        elseif method == :sison_glaz
            ci_sison_glaz(x,alpha, skew_correct=correct)
        elseif method == :quesenberry_hurst
            ci_quesenberry_hurst(x,alpha,GC=GC)
        elseif method == :bootstrap
            ci_bootstrap(x,alpha,bootstrap_iters)
        else
            throw(ArgumentError("method=$(method) is invalid or not implemented yet"))
        end
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

# Bootstrap
function ci_bootstrap(x::PowerDivergenceTest,alpha::Float64, iters::Int64)
    m = mapslices(x -> quantile(x, [alpha / 2, 1 - alpha / 2]), rand(Multinomial(x.n, convert(Vector{Float64}, x.thetahat)),iters) / x.n, dims=2)
    Tuple{Float64,Float64}[(boundproportion(m[i,1]), boundproportion(m[i,2])) for i in 1:length(x.thetahat)]
end

# Quesenberry, Hurst (1964)
function ci_quesenberry_hurst(x::PowerDivergenceTest, alpha::Float64; GC::Bool=true)
    m  = length(x.thetahat)
    cv = GC ? quantile(Chisq(1), 1 - alpha / m) : quantile(Chisq(m - 1), 1 - alpha)

    u = (cv .+ 2 .* x.observed .+ sqrt.(cv .* (cv .+ 4 .* x.observed .* (x.n .- x.observed) ./ x.n))) ./ (2 .* (x.n .+ cv))
    l = (cv .+ 2 .* x.observed .- sqrt.(cv .* (cv .+ 4 .* x.observed .* (x.n .- x.observed) ./ x.n))) ./ (2 .* (x.n .+ cv))
    Tuple{Float64,Float64}[(boundproportion(l[j]), boundproportion(u[j])) for j in 1:m]
end

# asymptotic simultaneous confidence interval
# Gold (1963)
function ci_gold(x::PowerDivergenceTest, alpha::Float64; correct::Bool=true, GC::Bool=true)
    m  = length(x.thetahat)
    cv = GC ? quantile(Chisq(1), 1 - alpha / 2m) : quantile(Chisq(m - 1), 1 - alpha / 2)

    u = [x.thetahat[j] + sqrt.(cv * x.thetahat[j] * (1 - x.thetahat[j]) / x.n) + ifelse(correct, inv(2x.n), 0) for j in 1:m]
    l = [x.thetahat[j] - sqrt.(cv * x.thetahat[j] * (1 - x.thetahat[j]) / x.n) - ifelse(correct, inv(2x.n), 0) for j in 1:m]
    Tuple{Float64,Float64}[ (boundproportion(l[j]), boundproportion(u[j])) for j in 1:m]
end

# Simultaneous confidence interval
# Sison, Glaz (1995)
# adapted from SAS macro by May, Johnson (2000)
function ci_sison_glaz(x::PowerDivergenceTest, alpha::Float64; skew_correct::Bool=true)
    k = length(x.thetahat)
    probn = inv(pdf(Poisson(x.n), x.n))

    c = 0
    p = 0.0
    p_old = 0.0
    m1, m2, m3, m4, m5 = zeros(k), zeros(k), zeros(k), zeros(k), zeros(k)
    mu = zeros(4)

    for _c in 1:x.n
        #run truncpoi
        for i in 1:k
            lambda = x.observed[i]
            #run moments
            a = lambda + _c
            b = max(lambda - _c, 0)
            poislama = cdf(Poisson(lambda), a)
            poislamb = cdf(Poisson(lambda), b - 1)
            den = b > 0.0 ? poislama-poislamb : poislama

            for r in 1:4
                plar = cdf(Poisson(lambda), a - r)
                plbr = cdf(Poisson(lambda), b - r - 1)

                poisA = ifelse( (a - r) >= 0, poislama - plar, poislama)
                poisB = 0.0
                if  (b - r - 1) >= 0
                    poisB = poislamb - plbr
                end
                if (b - r - 1) < 0 && b - 1 >= 0
                    poisB = poislamb
                end
                if (b - r - 1) < 0 && (b - 1) < 0
                    poisB = 0.0
                end
                mu[r] = lambda^r * (1 - (poisA - poisB) / den)
            end
            # end of moments
            m1[i] = mu[1]
            m2[i] = mu[2] + mu[1] - mu[1]^2
            m3[i] = mu[3] + mu[2] * (3 - 3mu[1]) + (mu[1] - 3mu[1]^2 + 2mu[1]^3)
            m4[i] = mu[4] + mu[3] * (6 - 4mu[1]) + mu[2] * (7 - 12mu[1] + 6mu[1]^2) + mu[1] - 4mu[1]^2 + 6mu[1]^3 - 3mu[1]^4
            m5[i] = den
        end
        for i in 1:k
            m4[i] -= 3m2[i]^2
        end
        s1, s2, s3, s4 = sum(m1), sum(m2), sum(m3), sum(m4)
        z  = (x.n - s1) / sqrt(s2)
        g1 = s3 / s2^(3/2)
        g2 = s4 / s2^2

        poly = 1 + g1 * (z^3 - 3z) / 6 + g2 * (z^4 - 6z^2 + 3) / 24 + g1^2 * (z^6 - 15z^4 + 45z^2 - 15) / 72
        f = poly * exp(-z^2 / 2) / sqrt(2π)
        probx = prod(m5)

        p = probn * probx * f / sqrt(s2)
        # end of truncpoi

        if p > 1 - alpha && p_old < 1 - alpha
            c = _c
            break
        else
            c = x.n
        end
        p_old = p
    end

    delta = (1 - alpha - p_old) / (p - p_old)
    out = zeros(k, 5)
    num = zeros(k, 1)

    c -= 1

    vol1 = 1
    vol2 = 1
    for i in 1:k
        num[i,1] = i
        cn = c / x.n
        onen = 1 / x.n
        out[i,1] = x.thetahat[i]

        out[i,2] = max(x.thetahat[i] - cn, 0)
        out[i,3] = min(x.thetahat[i] + cn + 2delta / x.n, 1)

        out[i,4] = max(x.thetahat[i] - cn - onen, 0)
        out[i,5] = min(x.thetahat[i] + cn + onen, 1)
    end
    if skew_correct
        return Tuple{Float64,Float64}[(boundproportion(out[i,2]), boundproportion(out[i,3])) for i in 1:k]
    else
        return Tuple{Float64,Float64}[(boundproportion(out[i,4]), boundproportion(out[i,5])) for i in 1:k]
    end
end

# power divergence statistic for testing goodness of fit
# Cressie and Read 1984; Read and Cressie 1988
# lambda    =  1: Pearson's Chi-square statistic
# lambda ->  0: Converges to Likelihood Ratio test stat
# lambda -> -1: Converges to minimum discrimination information statistic (Gokhale, Kullback 1978)
# lambda =  -2: Neyman Modified chi-squared statistic (Neyman 1949)
# lambda = -.5: Freeman-Tukey statistic (Freeman, Tukey 1950)

# Under regularity conditions, their asymptotic distributions are all the same (Drost 1989)
# Chi-squared null approximation works best for lambda near 2/3
"""
    PowerDivergenceTest(x[, y]; lambda = 1.0, theta0 = ones(length(x))/length(x))

Perform a Power Divergence test.

If `y` is not given and `x` is a matrix with one row or column, or `x` is a vector, then
a goodness-of-fit test is performed (`x` is treated as a one-dimensional contingency
table). In this case, the hypothesis tested is whether the population probabilities equal
those in `theta0`, or are all equal if `theta0` is not given.

If `x` is a matrix with at least two rows and columns, it is taken as a two-dimensional
contingency table. Otherwise, `x` and `y` must be vectors of the same length. The contingency
table is calculated using the `counts` function from the `StatsBase` package. Then the power
divergence test is conducted under the null hypothesis that the joint distribution of the
cell counts in a 2-dimensional contingency table is the product of the row and column
marginals.

Note that the entries of `x` (and `y` if provided) must be non-negative integers.

The power divergence test is given by
```math
    \\dfrac{2}{λ(λ+1)}\\sum_{i=1}^I \\sum_{j=1}^J n_{ij} \\left[(n_{ij}
    /\\hat{n}_{ij})^λ -1\\right]
```
where ``n_{ij}`` is the cell count in the ``i`` th row and ``j`` th column and ``λ`` is a
real number determining the nature of the test to be performed:

  * ``λ = 1``: equal to Pearson's chi-squared statistic
  * ``λ \\to 0``: converges to the likelihood ratio test statistic
  * ``λ \\to -1``: converges to the minimum discrimination information statistic
    (Gokhale and Kullback, 1978)
  * ``λ = -2``: equals Neyman modified chi-squared (Neyman, 1949)
  * ``λ = -1/2``: equals the Freeman-Tukey statistic (Freeman and Tukey, 1950).

Under regularity conditions, the asymptotic distributions are identical (see Drost et. al.
1989). The ``χ^2`` null approximation works best for ``λ`` near ``2/3``.

Implements: [`pvalue`](@ref), [`confint(::PowerDivergenceTest)`](@ref)

# References

  * Agresti, Alan. Categorical Data Analysis, 3rd Edition. Wiley, 2013.
"""
function PowerDivergenceTest(x::AbstractMatrix{T}; lambda::U=1.0, theta0::Vector{U} = ones(length(x))/length(x)) where {T<:Integer,U<:AbstractFloat}

    nrows, ncols = size(x)
    n = sum(x)

    #validate date
    (any(x .< 0) || any(!isfinite, x)) && throw(ArgumentError("all entries must be nonnegative and finite"))
    n == 0 && throw(ArgumentError("at least one entry must be positive"))
    (!isfinite(nrows) || !isfinite(ncols) || !isfinite(nrows*ncols)) && throw(ArgumentError("invalid number of rows or columns"))

    if nrows > 1 && ncols > 1
        rowsums = sum(x, dims=2)
        colsums = sum(x, dims=1)
        df = (nrows - 1) * (ncols - 1)
        thetahat = x ./ n
        xhat = rowsums * colsums / n
        theta0 = xhat / n
        V = Float64[(colsums[j]/n) * (rowsums[i]/n) * (1 - rowsums[i]/n) * (n - colsums[j]) for i in 1:nrows, j in 1:ncols]
    elseif nrows == 1 || ncols == 1
        df = length(x) - 1
        xhat = reshape(n * theta0, size(x))
        thetahat = x / n
        V = reshape(n .* theta0 .* (1 .- theta0), size(x))

        abs(1 - sum(theta0)) <= sqrt(eps()) || throw(ArgumentError("Probabilities must sum to one"))
    else
        throw(ArgumentError("Number of rows or columns must be at least 1"))
    end

    stat = 0
    if lambda == 0
        for i in 1:length(x)
            stat += x[i] * (log(x[i]) - log(xhat[i]))
        end
        stat *= 2
    elseif lambda == -1
        for i in 1:length(x)
            stat += xhat[i] * (log(xhat[i]) - log(x[i]))
        end
        stat *= 2
    else
        for i in 1:length(x)
            stat += x[i] * ((x[i] / xhat[i])^lambda - 1)
        end
        stat *= 2 / (lambda * (lambda + 1))
    end

    PowerDivergenceTest(lambda, vec(theta0), stat, df, x, n, vec(thetahat), xhat, (x - xhat) ./ sqrt.(xhat), (x - xhat) ./ sqrt.(V))
end

#convenience functions

#PDT
function PowerDivergenceTest(x::AbstractVector{T}, y::AbstractVector{T}, levels::Levels{T}; lambda::U=1.0) where {T<:Integer,U<:AbstractFloat}
    d = counts(x, y, levels)
    PowerDivergenceTest(d, lambda=lambda)
end

function PowerDivergenceTest(x::AbstractVector{T}, y::AbstractVector{T}, k::T; lambda::U=1.0) where {T<:Integer,U<:AbstractFloat}
    d = counts(x, y, k)
    PowerDivergenceTest(d, lambda=lambda)
end

PowerDivergenceTest(x::AbstractVector{T}; lambda::U=1.0, theta0::Vector{U} = ones(length(x))/length(x)) where {T<:Integer,U<:AbstractFloat} =
    PowerDivergenceTest(reshape(x, length(x), 1), lambda=lambda, theta0=theta0)

#ChisqTest
"""
    ChisqTest(x[, y][, theta0 = ones(length(x))/length(x)])

Perform a Pearson chi-squared test (equivalent to a [`PowerDivergenceTest`](@ref)
with ``λ = 1``).

If `y` is not given and `x` is a matrix with one row or column, or `x` is a vector, then
a goodness-of-fit test is performed (`x` is treated as a one-dimensional contingency
table). In this case, the hypothesis tested is whether the population probabilities equal
those in `theta0`, or are all equal if `theta0` is not given.

If `x` is a matrix with at least two rows and columns, it is taken as a two-dimensional
contingency table. Otherwise, `x` and `y` must be vectors of the same length. The contingency
table is calculated using `counts` function from the `StatsBase` package. Then the power
divergence test is conducted under the null hypothesis that the joint distribution of the
cell counts in a 2-dimensional contingency table is the product of the row and column
marginals.

Note that the entries of `x` (and `y` if provided) must be non-negative integers.

Implements: [`pvalue`](@ref), [`confint`](@ref)
"""
function ChisqTest(x::AbstractMatrix{T}) where T<:Integer
    PowerDivergenceTest(x, lambda=1.0)
end

function ChisqTest(x::AbstractVector{T}, y::AbstractVector{T}, levels::Levels{T}) where T<:Integer
    d = counts(x, y, levels)
    PowerDivergenceTest(d, lambda=1.0)
end

function ChisqTest(x::AbstractVector{T}, y::AbstractVector{T}, k::T) where T<:Integer
    d = counts(x, y, k)
    PowerDivergenceTest(d, lambda=1.0)
end

ChisqTest(x::AbstractVector{T}, theta0::Vector{U} = ones(length(x))/length(x)) where {T<:Integer,U<:AbstractFloat} =
    PowerDivergenceTest(reshape(x, length(x), 1), lambda=1.0, theta0=theta0)

#MultinomialLRTest
"""
    MultinomialLRTest(x[, y][, theta0 = ones(length(x))/length(x)])

Perform a multinomial likelihood ratio test (equivalent to a [`PowerDivergenceTest`](@ref)
with ``λ = 0``).

If `y` is not given and `x` is a matrix with one row or column, or `x` is a vector, then
a goodness-of-fit test is performed (`x` is treated as a one-dimensional contingency
table). In this case, the hypothesis tested is whether the population probabilities equal
those in `theta0`, or are all equal if `theta0` is not given.

If `x` is a matrix with at least two rows and columns, it is taken as a two-dimensional
contingency table. Otherwise, `x` and `y` must be vectors of the same length. The contingency
table is calculated using `counts` function from the `StatsBase` package. Then the power
divergence test is conducted under the null hypothesis that the joint distribution of the
cell counts in a 2-dimensional contingency table is the product of the row and column
marginals.

Note that the entries of `x` (and `y` if provided) must be non-negative integers.

Implements: [`pvalue`](@ref), [`confint`](@ref)
"""
function MultinomialLRTest(x::AbstractMatrix{T}) where T<:Integer
    PowerDivergenceTest(x, lambda=0.0)
end

function MultinomialLRTest(x::AbstractVector{T}, y::AbstractVector{T}, levels::Levels{T}) where T<:Integer
    d = counts(x, y, levels)
    PowerDivergenceTest(d, lambda=0.0)
end

function MultinomialLRTest(x::AbstractVector{T}, y::AbstractVector{T}, k::T) where T<:Integer
    d = counts(x, y, k)
    PowerDivergenceTest(d, lambda=0.0)
end

MultinomialLRTest(x::AbstractVector{T}, theta0::Vector{U} = ones(length(x))/length(x)) where {T<:Integer,U<:AbstractFloat} =
    PowerDivergenceTest(reshape(x, length(x), 1), lambda=0.0, theta0=theta0)

function show_params(io::IO, x::PowerDivergenceTest, ident="")
    println(io, ident, "Sample size:        $(x.n)")
    println(io, ident, "statistic:          $(x.stat)")
    println(io, ident, "degrees of freedom: $(x.df)")
    print(io, ident, "residuals:          ")
    show(io, vec(x.residuals))
    println(io)
    print(io, ident, "std. residuals:     ")
    show(io, vec(x.stdresiduals))
    println(io)
end


text:
using DifferentialEquations, Plots

a=3.0
f(y,a) =   sign(y)  + a -  y
f(y,a,t) = f(y,a) # just for the difeq solver

sol = solve(  ODEProblem( f, 3, (0,10.0), 0.2 ));
plot( sol , legend=false, background_color_inside=:black , ylims=(-7,7), lw=3, c=:red)

max_stars_repo_name: matutenun/projectbinder

max_issues_repo_name: matutenun/projectbinder

max_forks_repo_name: matutenun/projectbinder


text:
"""
    function calculate_metric(;
        metric::String,
        players::Array{Player},
        ranking_scheme1::String,
        ranking_scheme2::String
    )::Float64

Calculates the desired metric.

Arguments
=========
- `metric`:
    the metric to be calculated.
- `players`:
    the players whose rankings to consider.
- `ranking_scheme1`:
    ranking scheme 1 to use.
- `ranking_scheme2`:
    ranking scheme 2 to use.
"""

function calculate_metric(;
    metric::String,
    players::Array{Player},
    ranking_scheme1::String,
    ranking_scheme2::String
)::Float64
    if metric == "Spearman Correlation"
        return spearman_corr(
            players = players,
            ranking_scheme1 = ranking_scheme1,
            ranking_scheme2 = ranking_scheme2
        )
    elseif metric == "Kendall Tau Correlation"
        return kendtau_corr(
            players = players,
            ranking_scheme1 = ranking_scheme1,
            ranking_scheme2 = ranking_scheme2
        )
    end
    throw(DomainError(metric, "this metric isn't supported"))
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
# Apportion total learning to different human capital types

abstract type AbstractLearnSwitches end
abstract type AbstractLearningSet <: ModelObject end
abstract type AbstractLearning end

## ----------------  All colleges only teach college h

struct LearnCollegeOnlySet <: AbstractLearningSet
    objId :: ObjectId
end
struct LearnCollegeOnlySwitches <: AbstractLearnSwitches end
struct LearnCollegeOnly <: AbstractLearning end

ModelParams.has_pvector(l :: LearnCollegeOnlySet) = false;

init_learning_set(objId :: ObjectId, switches :: LearnCollegeOnlySwitches) =
    LearnCollegeOnlySet(objId);

init_learning(l :: LearnCollegeOnlySet, isTwoYear :: Bool) =
    LearnCollegeOnly();

frac_h(l :: LearnCollegeOnly) = 1.0;

# Map learning and endowments into human capital stocks at work start
h_stocks(l :: LearnCollegeOnly, dh, h0, hHs) =
    (dh .+ h0, hHs .+ zeros(size(dh)));

StructLH.describe(switches :: LearnCollegeOnlySwitches) =
    ["Learning"  "Colleges only teach college human capital"];

## ---------  Two and four year colleges

Base.@kwdef mutable struct LearnTwoFourSwitches <: AbstractLearnSwitches
    # Fraction learned that becomes hCollege for both types of colleges
    fracHc2 :: Double = 0.2
    fracHc4 :: Double = 0.8
    calFracHc2 :: Bool = false
    calFracHc4 :: Bool = false
    # 4y colleges produce more college h than 2y colleges
    isIncreasing :: Bool = true
end

mutable struct LearnTwoFourSet <: AbstractLearningSet
    objId :: ObjectId
    switches :: LearnTwoFourSwitches
    pvec :: ParamVector
    fracHc2 :: Double
    fracHc4 :: Double
end

struct LearnTwoFour <: AbstractLearning
    fracHc :: Double
end

function StructLH.describe(l :: LearnTwoFourSet)
    frac2y = round(frac_h(l, true); digits = 2);
    frac4y = round(frac_h(l, false); digits = 2);
    return [
        "Learning in college"  "2y/4y produce different combinations of h";
        "Fraction college h 2y"  "$frac2y";
        "Fraction college h 4y"  "$frac4y"
    ]
end

function StructLH.describe(switches :: LearnTwoFourSwitches)
    if switches.calFracHc2
        h2Str = "calibrated";
    else
        h2Str = "$(switches.fracHc2)"
    end
    if switches.calFracHc4
        h4Str = "calibrated";
    else
        h4Str = "$(switches.fracHc4)"
    end
    return ["Learning"  "Two and four year colleges teach different h mixes.";
    "Frac college h"  "Two yr: $h2Str   Four yr: $h4Str"];
end

# Displays parameters in levels, not as intercept and increments.
function ModelParams.param_table(l :: LearnTwoFourSet, isCalibrated :: Bool)
    if l.switches.isIncreasing
        if (l.switches.calFracHc2 == isCalibrated) ||
            (l.switches.calFracHc4 == isCalibrated)
            symbols = chain_strings(param_symbols(l), ", ");
            values = [frac_h(l, true), frac_h(l, false)];
            # This is where we get the description and symbol from
            p = l.pvec[1];
            pt = ModelParams.ParamTable(1);
            ModelParams.set_row!(pt, 1, string(p.name), symbols,
                "Fraction college h produced by quality",
                format_vector("", values, 2));
        else
            pt = nothing;
        end
    else
        # The standard reporting is easily readable.
        pt = ModelParams.report_params(l.pvec, isCalibrated);
    end
    return pt
end

param_symbols(l :: LearnTwoFourSet) = [p.symbol  for p in l.pvec];

function init_learning_set(objId :: ObjectId, switches :: LearnTwoFourSwitches)
    fracHc2 = switches.fracHc2;
    pFrac2 = Param(:fracHc2, ldescription(:fracHcTwo), lsymbol(:fracHcTwo),
        fracHc2, fracHc2, 0.0, 1.0, switches.calFracHc2);

    fracHc4 = switches.fracHc4;
    descrStr = ldescription(:fracHcFour);
    if switches.isIncreasing
        descrStr = descrStr * " (increment)";
    end
    pFrac4 = Param(:fracHc4, descrStr, lsymbol(:fracHcFour),
        fracHc4, fracHc4, 0.0, 1.0, switches.calFracHc4);

    pvec = ParamVector(objId, [pFrac2, pFrac4]);
    return LearnTwoFourSet(objId, switches, pvec, fracHc2, fracHc4);
end

init_learning(l :: LearnTwoFourSet, isTwoYear :: Bool) =
    LearnTwoFour(frac_h(l, isTwoYear));

# Fraction of learning that goes to college h
frac_h(switches :: LearnTwoFourSwitches, isTwoYear :: Bool) =
    frac_h(switches.isIncreasing, switches.fracHc2, switches.fracHc4, isTwoYear);
frac_h(l :: LearnTwoFourSet, isTwoYear :: Bool) =
    frac_h(l.switches.isIncreasing, l.fracHc2, l.fracHc4, isTwoYear);
frac_h(l :: LearnTwoFour) = l.fracHc;

function frac_h(isIncreasing, fracHc2, fracHc4, isTwoYear)
    if isTwoYear
        fracH = fracHc2;
    elseif isIncreasing
        fracH = fracHc2 + fracHc4 * (1.0 - fracHc2);
    else
        fracH = fracHc4;
    end
    return fracH
end

# Map learning and endowments into human capital stocks at work start
h_stocks(l :: LearnTwoFour, dh, h0, hHs) =
    (h0 .+ frac_h(l) .* dh, hHs .+ (1.0 - frac_h(l)) .* dh);

## ---------  Generic

make_test_learning() = LearnCollegeOnly();

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


text:
cd(abs(parent(Path(@__FILE__)))) do
    @testset "Simple Path Usage" begin
        reg = is_windows() ? "..\\src\\FilePaths.jl" : "../src/FilePaths.jl"
        @test ispath(reg)

        p = Path(reg)

        @test p == p"../src/FilePaths.jl"
        @test String(p) == reg
        @test String(cwd()) == pwd()
        @test String(home()) == homedir()

        @test parts(p) == ("..", "src", "FilePaths.jl")
        @test hasparent(p)
        @test parent(p) == p"../src"
        @test parents(p) == [p"..", p"../src"]
        @test_throws ErrorException parents(p".")

        @test basename(p) == "FilePaths.jl"
        @test join(parent(p), Path(basename(p))) == p
        @test filename(p) == "FilePaths"

        @test extension(p) == "jl"
        @test extension(p"../REQUIRE") == ""
        @test extensions(p"foo.tar.gz") == ["tar", "gz"]
        @test length(extensions(p"../REQUIRE")) == 0

        @test exists(p)
        @test !isabs(p)
        @test String(norm(p"../src/../src/FilePaths.jl")) == normpath("../src/../src/FilePaths.jl")
        @test String(abs(p)) == abspath(String(p))
        # This works around an issue with Base.relpath: that function does not take
        # into account the paths on Windows should be compared case insensitive.
        homedir_patched = homedir()
        if is_windows()
            conv_f = isupper(abspath(String(p))[1]) ? uppercase : lowercase
            homedir_patched = conv_f(homedir_patched[1]) * homedir_patched[2:end]
        end
        @test String(relative(p, home())) == relpath(String(p), homedir_patched)

        @test isa(relative(Path(".")), AbstractPath)
        @test relative(Path(".")) == Path(".")

        s = stat(p)
        lstat(p)

        show_str = sprint(show, s)
        #@test "device" in show_str
        #@test "blocks" in show_str

        @test size(p) == stat(p).size
        @test modified(p) == stat(p).mtime
        @test created(p) == stat(p).ctime

        @test isfile(p)
        @test isdir(parent(p))
        @test !islink(p)
        @test !issocket(p)
        @test !isfifo(p)
        @test !ischardev(p)
        @test !isblockdev(p)
    end
end

mktmpdir() do d
    cd(d) do
        @testset "Modifying Path Usage" begin
            new_path = p"foo"

            mkdir(new_path)
            @test_throws ErrorException mkdir(new_path)
            remove(new_path)

            new_path = p"foo/bar"
            @test_throws ErrorException mkdir(new_path)
            mkdir(new_path; recursive=true)
            mkdir(new_path; recursive=true, exist_ok=true)

            other_path = p"car/bar"
            copy(new_path, other_path; recursive=true)
            copy(new_path, other_path; exist_ok=true, overwrite=true)
            @test_throws ErrorException copy(new_path, other_path)
            @test_throws ErrorException copy(p"badpath", other_path; exist_ok=true, overwrite=true)
            remove(p"car"; recursive=true)

            move(new_path, other_path; recursive=true)
            mkdir(new_path; recursive=true)
            move(new_path, other_path; exist_ok=true, overwrite=true)
            @test_throws ErrorException move(new_path, other_path)
            @test_throws ErrorException move(p"badpath", other_path; exist_ok=true, overwrite=true)
            remove(p"car"; recursive=true)

            mkdir(new_path; recursive=true)

            symlink(new_path, p"mysymlink")
            symlink(new_path, p"mysymlink"; exist_ok=true, overwrite=true)
            @test_throws ErrorException symlink(new_path, p"mysymlink")
            @test_throws ErrorException symlink(p"badpath", p"mysymlink"; exist_ok=true, overwrite=true)

            touch(p"newfile")
            mktmp(d) do f, io
                println(f)
                println(io)
            end

            @static if is_unix()
                if haskey(ENV, "USER")
                    if ENV["USER"] == "root"
                        chown(p"newfile", "nobody", "nogroup"; recursive=true)
                    else
                        @test_throws ErrorException chown(p"newfile", "nobody", "nogroup"; recursive=true)
                    end
                end
            else
                @test_throws ErrorException chown(p"newfile", "nobody", "nogroup"; recursive=true)
            end

            @static if is_unix()
                chmod(p"newfile", user=(READ+WRITE+EXEC), group=(READ+EXEC), other=READ)
                @test String(mode(p"newfile")) == "-rwxr-xr--"
                @test isexecutable(p"newfile")
                @test iswritable(p"newfile")
                @test isreadable(p"newfile")

                chmod(p"newfile", "-x")
                @test !isexecutable(p"newfile")

                @test String(mode(p"newfile")) == "-rw-r--r--"
                chmod(p"newfile", "+x")
                write(p"newfile", "foobar")
                @test read(p"newfile") == "foobar"
                chmod(p"newfile", "u=rwx")

                chmod(new_path, mode(p"newfile"); recursive=true)
            end
        end
    end
end