text:
export judiProjection, judiWavelet, judiLRWF, judiRHS

############################################################################################################################
# Abstract types
abstract type judiNoopOperator{D} <: joAbstractLinearOperator{D, D} end

############################################################################################################################
# Concrete types
# Lazy adjoint like in LinearAlgebra
struct jAdjoint{T}
    op::T
end

# Projection operator
struct judiProjection{D} <: judiNoopOperator{D}
    m::AbstractSize
    n::AbstractSize
    geometry::Geometry
end

# Wavelet "projection"
struct judiWavelet{D} <: judiNoopOperator{D}
    m::AbstractSize
    n::AbstractSize
    wavelet::Vector{Array{D, N}} where N
    dt::Vector{D}
end

# Poorly named backward compat
const judiLRWF{T} = judiWavelet{T}
const Projection{D} = Union{judiProjection{D}, judiWavelet{D}}
const AdjointProjection{D} = jAdjoint{<:Projection{D}}

"""
    judiRHS
        dt::Vector{T}
        geometry::Geometry
        data
Abstract sparse vector for right-hand-sides of the modeling operators. The `judiRHS` vector has the\\
dimensions of the full time history of the wavefields, but contains only the data defined at the \\
source or receiver positions (i.e. wavelets or shot records).

Constructor
==========
    judiRHS(dt, geometry, data)

Examples
========
Assuming `Pr` and `Ps` are projection operators of type `judiProjection` and `dobs` and `q` are\\
seismic vectors of type `judiVector`, then a `judiRHS` vector can be created as follows:
    rhs = Pr'*dobs    # right-hand-side with injected observed data
    rhs = Ps'*q    # right-hand-side with injected wavelet
"""
struct judiRHS{D} <: judiMultiSourceVector{D}
    nsrc::Integer
    P::jAdjoint{judiProjection{D}}
    d::judiVector
end

"""
    LazyAdd
        nsrc
        A
        B
        sign

Lazy addition of two RHS (currently only judiVector). The addition isn't evaluated to avoid
large memory allocation but instead evaluates the addition (with sign `sign`) `A + sign * B`
for a single source at propagation time.
"""
struct LazyAdd{D} <: judiMultiSourceVector{D}
    nsrc::Integer
    A
    B
    sign
end

############################################################################################################################
# Constructors
"""
    judiProjection(geometry)

Projection operator for sources/receivers to restrict or inject data at specified locations.

Examples
========
`F` is a modeling operator of type `judiModeling` and `q` is a seismic source of type `judiVector`:
    Pr = judiProjection(rec_geometry)
    Ps = judiProjection(q.geometry)
    dobs = Pr*F*Ps'*q
    qad = Ps*F'*Pr'*dobs
"""
judiProjection(G::Geometry) = judiProjection{Float32}(rec_space(G), time_space_src(get_nsrc(G), G.nt, 3), G)

"""
    judiWavelet(dt, wavelet)

Low-rank wavefield operator which injects a wavelet q at every point of the subsurface.

Examples
========
`F` is a modeling operator of type `judiModeling` and `w` is a weighting matrix of type `judiWeights`:
    Pr = judiProjection(rec_geometry)
    Pw = judiWavelet(rec_geometry.dt, q.data) # or judiLRWF(rec_geometry.dt, q.data)
    dobs = Pr*F*Pw'*w
    dw = Pw*F'*Pr'*dobs
"""
judiWavelet(nsrc::Integer, dt::T, w::Array{T, N}) where {T<:AbstractFloat, N} = judiWavelet{Float32}(space_src(nsrc), time_space_src(nsrc, length.(w)), [w for i=1:nsrc], [dt for i=1:nsrc])
judiWavelet(dt::T, w::Array{T, N}) where {T<:AbstractFloat, N} = judiWavelet(1, dt, w)
judiWavelet(dt::Vector{<:Number}, w::Vector{T}) where T<:Array = judiWavelet{Float32}(space_src(length(dt)), time_space_src(length(dt), length.(w)), w, dt)
judiWavelet(dt::Vector{<:Number}, w::Vector{<:Number}) = judiWavelet{Float32}(space_src(length(dt)), time_space_src(length(dt), [length(w) for i=1:length(dt)]), [w for i=1:length(dt)], dt)
judiWavelet(dt::dtT, w::Array{T, N}) where {dtT<:Number, T<:Array, N} = judiWavelet([dt for i=1:length(w)], w)

# Deprecation error
judiWavelet(::Integer, ::Array{T, N}) where {T<:AbstractFloat, N} = throw(ArgumentError("Time sampling of the wavelet need to be suplied `judiWavelet(nsrc, dt, wavelet)`"))

############################################################################################################################
# Base overload
getproperty(P::judiProjection, s::Symbol) = s == :data ? P.geometry : getfield(P, s)
getproperty(P::judiWavelet, s::Symbol) = s == :data ? P.wavelet : getfield(P, s)

function getproperty(jA::jAdjoint, s::Symbol)
    s == :m && (return jA.op.n)
    s == :n && (return jA.op.m)
    s == :op && (return getfield(jA, s))
    return getproperty(jA.op, s)
end

size(jA::jAdjoint) = (jA.op.n, jA.op.m)
display(P::jAdjoint) = println("Adjoint($(P.op))")
display(P::judiProjection{D}) where D = println("JUDI projection operator $(repr(P.n)) -> $(repr(P.m))")

############################################################################################################################
# Indexing
getindex(jA::jAdjoint{T}, i) where T = jAdjoint{T}(jA.op[i])
getindex(P::judiProjection{D}, i) where D = judiProjection{D}(P.m[i], P.n[i], P.geometry[i])
getindex(P::judiProjection{D}, i::Integer) where D = judiProjection{D}(P.m[i], P.n[i], P.geometry[i:i])
getindex(P::judiWavelet{D}, i) where D = judiWavelet{D}(P.m[i], P.n[i], P.wavelet[i], P.dt[i])
getindex(P::judiWavelet{D}, i::Integer) where D = judiWavelet{D}(P.m[i], P.n[i], P.wavelet[i:i], P.dt[i:i])
getindex(rhs::judiRHS{D}, i::Integer) where D = judiRHS{D}(length(i), rhs.P[i], rhs.d[i])
getindex(rhs::judiRHS{D}, i::RangeOrVec) where D = judiRHS{D}(length(i), rhs.P[i], rhs.d[i])
getindex(la::LazyAdd{D}, i::RangeOrVec) where D = LazyAdd{D}(length(i), la.A[i], la.B[i], la.sign)

# Backward compatible subsample
subsample(P::judiNoopOperator{D}, i) where D = getindex(P, i)

############################################################################################################################
# Linear algebra
adjoint(jA::jAdjoint) = jA.op
adjoint(P::judiNoopOperator{D}) where D = jAdjoint(P)

transpose(jA::jAdjoint) = jA.op
transpose(P::judiNoopOperator{D}) where D = jAdjoint(P)

conj(jA::jAdjoint) = jA
conj(P::judiNoopOperator{D}) where D = P

*(P::jAdjoint{judiProjection{D}}, d::judiVector{D, AT}) where {D, AT} = judiRHS{D}(d.nsrc, P, d)

# Combinations of rhs
for (op, s) in zip([:+, :-], (1, -1))
    for T1 in [judiRHS, LazyAdd]
        for T2 in [judiRHS, LazyAdd]
            @eval function $(op)(r1::$(T1){D}, r2::$(T2){D}) where D
                r1.nsrc == r2.nsrc || throw(ArgumentError("Incompatible number of source experiment"))
                LazyAdd{D}(r1.nsrc, r1, r2, $s)
            end
        end
    end
end

############################################################################################################################
# Comparison
==(jA1::jAdjoint, jA2::jAdjoint) = jA1.op  == jA2.op
==(P1::judiProjection, P2::judiProjection) = P1.data == P2.data
==(P1::judiWavelet, P2::judiWavelet) = (P1.wavelet == P2.wavelet && P1.dt == P2.dt)

############################################################################################################################
###### Inpout processing for propagation
make_src(q, P::jAdjoint{judiProjection{D}}) where D = (P.data[1], make_input(q))
make_src(q, P::jAdjoint{judiLRWF{D}}) where D = (make_input(q), P.data[1])
make_src(q, P::judiProjection{D}) where D = (P.data[1], make_input(q))
make_src(q, P::judiLRWF{D}) where D = (make_input(q), P.data[1])
make_src(rhs::judiRHS) = make_src(rhs.d, rhs.P)

get_nsrc(P::Projection) = P.m[:src]
get_nsrc(P::jAdjoint{<:Projection}) = P.op.m[:src]
get_nt(P::Projection) = P.n[:time]
get_nt(P::jAdjoint{<:Projection}) = P.op.n[:time]

############################################################################################################################
###### Evaluate lazy operation
eval(rhs::judiRHS) = rhs.d

function eval(ls::LazyAdd{D}) where D
    aloc = eval(ls.A)
    bloc = eval(ls.B)
    ga = aloc.geometry
    gb = bloc.geometry
    @assert (ga.nt == gb.nt && ga.dt == gb.dt && ga.t == gb.t)
    xloc = [vcat(ga.xloc[1], gb.xloc[1])]
    yloc = [vcat(ga.yloc[1], gb.yloc[1])]
    zloc = [vcat(ga.zloc[1], gb.zloc[1])]
    geom = GeometryIC{D}(xloc, yloc, zloc, ga.dt, ga.nt, ga.t)
    data = hcat(aloc.data[1], ls.sign*bloc.data[1])
    judiVector{D, Matrix{D}}(1, geom, [data])
end

function make_src(ls::LazyAdd{D}) where D
    q = eval(ls)
    return q.geometry[1], q.data[1]
end


text:
using Plots

gr()
include("PRM.jl")

plot([1,2],[2,3])

start = rrt.Point(0,0)
goal = rrt.Point(18,18)

obstacleList = Vector{rrt.Obstacle}()
obs1 = rrt.Obstacle(rrt.Point(8,3),rrt.Point(10,18)) #Todo
push!(obstacleList,obs1)
#obs2 = rrt.Obstacle(rrt.Point(15,13),rrt.Point(17,15)) #Todo
#push!(obstacleList,obs2)

print("HIIIIIIIIIIIIII")
isCollidingRects(obs1, obs2)

#iterList = [40 80 100 120 150 180 200 300]
nTrials = 30
cost = 0
listCosts = Vector{Float32}()
listpSucc= Vector{Float32}()

print("\n\n")
connectDist = 10
iterList = [20 30 50 100 150 200 300]

for testIter in iterList
    totalcost =0
    idx = 0
    nSuccess = 0
    cost = 0
    while idx < nTrials
        nodeslist, edgeslist = preprocessPRM(testIter, connectDist, obstacleList)
        cost, isPathFound, winPath = queryPRM(start, goal, nodeslist, edgeslist)
        if isPathFound
            nSuccess += 1
        end
        totalcost += cost
        idx += 1
    end
    avgCost = totalcost /  nSuccess
    pSucc = nSuccess / nTrials
    #@printf("For the iter of %d the avg cost was %d across %d trials", maxIter, avgCost, nTrials)
    push!(listCosts, avgCost)
    push!(listpSucc, pSucc)
    @show listCosts

end

sizeplot = (400, 600)
pPRMcost = plot(iterList, listCosts', show=true, size=sizeplot, seriestype=:scatter,
     legend=false, yaxis=((0,100), 0:20:100), xaxis=((0,320), 0:50:300), color=:black)

title!("PRM nIter vs pathcost")
xlabel!("iterations requested of PRM alg")
ylabel!("euclidean path cost")

pPRMsuccess = scatter(iterList, listpSucc', size=sizeplot, legend=false, yaxis=((0,1), 0:0.1:1), xaxis=((0,1), 0:50:300), color=:orange)
title!("RRM nIter vs pSuccess")
ylabel!("# P(success) = numSucc / 30 trials")

plot(pPRMcost, pPRMsuccess, layout=(2,1),legend=false)

############################
############################

obstacleList = Vector{rrt.Obstacle}()
#obs1 = rrt.Obstacle(rrt.Point(8,3),rrt.Point(10,18)) #Todo
push!(obstacleList, obs1)

obs2 = rrt.Obstacle(rrt.Point(15,13),rrt.Point(17,15)) #Todo
push!(obs2)

ext: jl


text:
using Random,StatsBase,Distributions,Polynomials
Random.seed!(150000)

#LSM原始算法，针对美式PUT
function LSM_PRIMAL(S0::Float64,K::Float64,T::Float64,r::Float64,σ::Float64,I::Int,M::Int)
    """

    # Arguments
    - `S0`: 初始标的值.
    - `K`: 行权值
    - `T`: 到期时间，年为单位
    - `r`: 无风险利率
    - `σ`: 波动率
    - `I`: 原始算法中路径数
    - `M`: 时间步数

    """

    dt=T/M
    df=exp(-r*dt)

    p1=(r-0.5*(σ^2))*dt
    p2=rand(Normal(0,1),M+1,I)
    p3=p1 .+σ*√(dt)*p2
    p4=cumsum(p3,dims=1)

    S=S0*(p4.|>exp)
    S[1,:]=repeat([S0],I)
    h=map(x->max(K-x,0),S)

    V=copy(h[end,:])

    @simd for t=M:-1:1
        rg=polyfit(S[t,:],V*df,5)
        C=map(x->polyval(rg,x),S[t,:])
        res=h[t,:].>C
        for i=1:I
            if res[i]==true
                V[i]=h[t,i]
            else
                V[i]=V[i]*df
            end
        end
    end

    v0=df*sum(V)/I

end

#LSM对偶算法
function LSM_DUAL(S0::Float64,K::Float64,T::Float64,r::Float64,σ::Float64,J::Int,I1::Int,I2::Int,M::Int,reg::Int)
    """

    # Arguments
    - `S0`: 初始::Int标的值::Int.
    ::Int- ::Int`K`: 行权值
    - `T`: 到期时间，年为单位
    - `r`: 无风险利率
    - `σ`: 波动率
    - `J`: 对偶算法中子循环模拟路径
    - `I1`: 原始算法中路径数
    - `I2`: 对偶算法中路径数
    - `M`: 时间步数
    - `reg`: 多项是回归阶数

    """
    dt=T/M
    df=exp(-r*dt)

    #路径I1,产生回归系数
    S1=genrate_path(S0,r,σ,dt,M,I1)
    h1=map(x->max(K-x,0),S1)
    V1=map(x->max(K-x,0),S1)

    #路径I1,产生回归系数
    rg0=[]
    a=Poly(0)
    push!(rg0,a)
    @simd for t=M:-1:2
        rgr=polyfit(S1[t,:],V1[t+1,:]*df,5)
        push!(rg0,rgr)
    end
    push!(rg0,a)
    rg=reverse(rg0)

    #路径I2，对偶算法
    Q=zeros(M+1,I2)
    U=zeros(M+1,I2)
    S=genrate_path(S0,r,σ,dt,M,I2)
    h=map(x->max(K-x,0),S)
    V=map(x->max(K-x,0),S)

    @simd for t=2:M+1
        for i=1:I2
            tmp=polyval(rg[t],S[t,i])
            Vt=max(h[t,i],tmp)
            St=generate_nest_mc(r,σ,dt,S[t-1,i],J)
            Ct=map(x->polyval(rg[t],x),St)
            ht=map(x->max(K-x),St)
            res=ht.>Ct
            Vtj=zeros(J)
            for k=1:length(res)
                if res[k]==true
                    Vtj[k]=ht[k]
                else
                    Vtj[k]=Ct[k]
                end
            end
            sum_VT=sum(Vtj)/length(St)

            Q[t,i]=Q[t-1,i]/df+(Vt-sum_VT)
            U[t,i]=max(U[t-1,i]/df,(h[t,i]-Q[t,i]))

            if t==M+1
                U[t,i]=max(U[t-1,i]/df,mean(ht)-Q[t,i])
            end
        end
    end
    k=df^M
    U0=sum(U[M,:])/I2*(df^M)

end

#产生标的仿真路径-单路径
function generate_nest_mc(r::Float64,σ::Float64,dt::Float64,st::Float64,J::Int)

    p1=(r-0.5*(σ^2))*dt
    p2=rand(Normal(0,1),J)
    p3=p1 .+σ*√(dt)*p2
    S=st*(p3.|>exp)

    return S

end

#产生标的仿真路径-多路径
function genrate_path(S0::Float64,r::Float64,σ::Float64,dt::Float64,M::Int,I::Int)
    p1=(r-0.5*(σ^2))*dt
    S=zeros(M+1,I)
    S[1,:]=map(x->x=S0,S[1,:])
    @simd for t=2:M+1
        p2=rand(Normal(0,1),I)
        p3=p1 .+σ*√(dt)*p2
        S[t,:]=S[t-1,:].*(p3.|>exp)
    end

    return S
end

S0=36.0
K=40.0
T=1.0
r=0.06
σ=0.2
I1=16384
I2=1024

M=10
J=50
reg=5

res=LSM_DUAL(S0,K,T,r,σ,J,I1,I2,M,reg)
println(res)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct ConstValueDistribution{VF<:VariateForm,T} <: Distribution{VF,Continuous}
    value::T
end

ConstValueDistribution(x::T) where {T<:Real} = ConstValueDistribution{Univariate,T}(x)
ConstValueDistribution(x::T) where {T<:AbstractVector{<:Real}} = ConstValueDistribution{Multivariate,T}(x)
ConstValueDistribution(x::T) where {T<:AbstractMatrix{<:Real}} = ConstValueDistribution{Matrixvariate,T}(x)

Distributions.pdf(d::ConstValueDistribution{Univariate}, x::Real) = d.value == x ? float(eltype(d))(Inf) : float(eltype(d))(0)
Distributions._logpdf(d::ConstValueDistribution, x::AbstractArray{<:Real}) = d.value == x ? float(eltype(d))(Inf) : float(eltype(d))(0)

Distributions.cdf(d::ConstValueDistribution{Univariate}, x::Real) = d.value <= x ? Float32(1) : Float32(0)
Distributions.quantile(d::ConstValueDistribution{Univariate}, q::Real) = d.value # Sensible?
Distributions.minimum(d::ConstValueDistribution{Univariate}) = x.value
Distributions.maximum(d::ConstValueDistribution{Univariate}) = x.value
Distributions.insupport(d::ConstValueDistribution{Univariate}, x::Real) = x == d.value

Base.size(d::ConstValueDistribution) = size(d.value)
Base.length(d::ConstValueDistribution) = prod(size(d))
Base.eltype(d::ConstValueDistribution) = eltype(d.value)

Random.rand(rng::AbstractRNG, d::ConstValueDistribution) = d.value

function Random._rand!(rng::AbstractRNG, d::ConstValueDistribution, params::AbstractArray{<:Real})
    copyto!(params, d.value)
end

_np_valshape(d::Distribution{Univariate,Continuous}) = ScalarShape{Real}()
_np_bounds(d::Distribution{Univariate,Continuous}) = HyperRectBounds([quantile(d, 0f0)], [quantile(d, 1f0)], reflective_bounds)

_np_valshape(d::Distribution{Multivariate,Continuous}) = ArrayShape{Real}(size(d)...)
_np_bounds(d::Distribution{Multivariate,Continuous}) = HyperRectBounds(fill(Float32(-Inf), length(d)), fill(Float32(+Inf), length(d)), hard_bounds)

_np_bounds(d::ConstValueDistribution) = HyperRectBounds(Vector{eltype(d)}(), Vector{eltype(d)}(), hard_bounds)
_np_bounds(d::Product{Continuous}) = HyperRectBounds(quantile.(d.v, 0f0), quantile.(d.v, 1f0), reflective_bounds)

_np_valshape(d::Distribution{Matrixvariate,Continuous}) = ArrayShape{Real}(size(d)...)

_np_dist_shape_bounds(d::Distribution) = (d, _np_valshape(d), _np_bounds(d))

_np_distribution(s::ConstValueShape) = ConstValueDistribution(s.value)
_np_bounds(s::ConstValueShape) = HyperRectBounds(Float64[], Float64[], hard_bounds)
_np_dist_shape_bounds(s::ConstValueShape) = (_np_distribution(s), s, _np_bounds(s))

_np_dist_shape_bounds(s::AbstractInterval) = _np_dist_shape_bounds(Uniform(minimum(s), maximum(s)))
_np_dist_shape_bounds(xs::AbstractVector{<:AbstractInterval}) = _np_dist_shape_bounds(Product((s -> Uniform(minimum(s), maximum(s))).(xs)))
_np_dist_shape_bounds(xs::AbstractVector{<:Distribution}) = _np_dist_shape_bounds(Product(xs))
_np_dist_shape_bounds(x::Number) = _np_dist_shape_bounds(ConstValueShape(x))
_np_dist_shape_bounds(x::AbstractArray{<:Number}) = _np_dist_shape_bounds(ConstValueShape(x))

"""
    NamedTupleDist <: MultivariateDistribution
    NamedTupleDist <: MultivariateDistribution

A distribution with `NamedTuple`-typed variates.

Can be used to describe the distribution of each parameter in a set of
named parameters. If the distribution is used as a Bayesian prior,
the `NamedTupleDist` then specifies the prior on each named parameter.

A `NamedTupleDist` implies a `NamedTupleShape`:

    valshape(d::NamedTupleDist) isa NamedTupleShape
"""
struct NamedTupleDist{
    names,
    DT <: (NTuple{N,Distribution} where N),
    AT <: (NTuple{N,ValueShapes.ValueAccessor} where N),
    BT <: (NTuple{N,AbstractParamBounds} where N)
} <: Distribution{Multivariate,Continuous}
    _distributions::NamedTuple{names,DT}
    _shapes::NamedTupleShape{names,AT}
    _bounds::NamedTuple{names,BT}
end

export NamedTupleDist

function NamedTupleDist(param_priors::NamedTuple{names}) where {names}
    dsb = map(_np_dist_shape_bounds, param_priors)
    NamedTupleDist(
        map(x -> x[1], dsb),
        NamedTupleShape(map(x -> x[2], dsb)),
        map(x -> x[3], dsb)
    )
end

@inline NamedTupleDist(;named_priors...) = NamedTupleDist(values(named_priors))

@inline _distributions(p::NamedTupleDist) = getfield(p, :_distributions)
@inline _shapes(p::NamedTupleDist) = getfield(p, :_shapes)
@inline _bounds(p::NamedTupleDist) = getfield(p, :_bounds)

@inline Base.keys(p::NamedTupleDist) = keys(_distributions(p))

@inline Base.values(p::NamedTupleDist) = values(_distributions(p))

@inline function Base.getproperty(p::NamedTupleDist, s::Symbol)
    # Need to include internal fields of NamedTupleShape to make Zygote happy:
    if s == :_distributions
        getfield(p, :_distributions)
    elseif s == :_shapes
        getfield(p, :_shapes)
    elseif s == :_bounds
        getfield(p, :_bounds)
    else
        getproperty(_distributions(p), s)
    end
end

@inline Base.propertynames(p::NamedTupleDist) = propertynames(_distributions(p))

ValueShapes.valshape(p::NamedTupleDist) = _shapes(p)

param_bounds(p::NamedTupleDist) = vcat(map(_np_bounds, values(p))...)

_np_length(d::Distribution) = length(d)
_np_length(d::ConstValueDistribution) = 0

Base.length(p::NamedTupleDist) = sum(_np_length, values(p))

function _np_logpdf(
    dist::ConstValueDistribution,
    acc::ValueShapes.ValueAccessor{<:ConstValueShape},
    params::AbstractVector{<:Real}
)
    float(zero(eltype(params)))
end

function _np_logpdf(
    dist::Distribution,
    acc::ValueShapes.ValueAccessor,
    params::AbstractVector{<:Real}
)
    logpdf(dist, float(params[acc]))
end

function Distributions.logpdf(p::NamedTupleDist, params::AbstractVector{<:Real})
    distributions = values(p)
    accessors = values(valshape(p))
    sum(map((dist, acc) -> _np_logpdf(dist, acc, params), distributions, accessors))
end

# ConstValueDistribution has no dof, so NamedTupleDist logpdf contribution must be zero:
_np_logpdf(dist::ConstValueDistribution, params::Any) = zero(Float32)

_np_logpdf(dist::Distribution, params::Any) = logpdf(dist, params)

function Distributions.logpdf(p::NamedTupleDist{names}, params::NamedTuple{names}) where names
    distributions = values(p)
    parvalues = values(params)
    sum(map((dist, p) -> _np_logpdf(dist, p), distributions, parvalues))
end

function _np_rand!(
    rng::AbstractRNG, dist::ConstValueDistribution,
    acc::ValueShapes.ValueAccessor{<:ConstValueShape},
    params::AbstractVector{<:Real}
)
    nothing
end

function _np_rand!(
    rng::AbstractRNG, dist::Distribution,
    acc::ValueShapes.ValueAccessor,
    params::AbstractVector{<:Real}
)
    rand!(rng, dist, view(params, acc))
    nothing
end

function Distributions._rand!(rng::AbstractRNG, p::NamedTupleDist, params::AbstractVector{<:Real})
    distributions = values(p)
    accessors = values(valshape(p))
    map((dist, acc) -> _np_rand!(rng, dist, acc, params), distributions, accessors)
    params
end

#Random.rand(rng::AbstractRNG, p::NamedTupleDist) = rand!(rng, p, Vector{Float64}(undef, length(p)))

function _np_var_or_cov!(A_cov::AbstractArray{<:Real,0}, dist::Distribution{Univariate})
    A_cov[] = var(dist)
    nothing
end

function _np_var_or_cov!(A_cov::AbstractArray{<:Real,2}, dist::Distribution{Multivariate})
    A_cov[:, :] = cov(dist)
    nothing
end

function _np_cov!(
    dist::ConstValueDistribution,
    acc::ValueShapes.ValueAccessor{<:ConstValueShape},
    A_cov::AbstractMatrix{<:Real}
)
    nothing
end

function _np_cov!(
    dist::Distribution,
    acc::ValueShapes.ValueAccessor,
    A_cov::AbstractMatrix{<:Real}
)
    _np_var_or_cov!(view(A_cov, acc, acc), dist)
    nothing
end

function Statistics.cov(p::NamedTupleDist)
    let n = length(p), A_cov = zeros(n, n)
        distributions = values(p)
        accessors = values(valshape(p))
        map((dist, acc) -> _np_cov!(dist, acc, A_cov), distributions, accessors)
        A_cov
    end
end


text:
using SparseInverseProblems
import SparseInverseProblems: lmo, phi, solveFiniteDimProblem, localDescent
using SparseInverseProblems.Util

struct MatrixCompletion <: ForwardModel
  m :: Int64
  k :: Int64
  nnz :: Int64
  nz_i :: Vector{Int64}
  nz_j :: Vector{Int64}
end

function solveFiniteDimProblem(model :: MatrixCompletion, loss :: Loss, parameters :: Matrix{Float64}, y :: Vector{Float64}, tau :: Float64)
  A = runForwardSimulator(model, parameters)
  return nnlasso_kernel(A'*A,A'y,tau)
end

function localDescent(sim :: MatrixCompletion, loss :: LSLoss, thetas :: Matrix{Float64}, w :: Vector{Float64}, y :: Vector{Float64})
  rank = size(thetas,2)
  m,k = sim.m,sim.k
  u = thetas[1:m,:]
  v = thetas[(m+1):end,:]
  #M = uv'
  # min ||Phi(u*v') -y ||^2
  # ||u_i|| \le sqrt(w[i])
  # ||v_i|| \le sqrt(w[i])
  #GD w/ backtracking linesearch on the sphere.
  alpha = 0.05
  beta = 0.5
  c_1 = 1E-4
    for iter = 1:1
      (f,g_u,g_v) = fg_c(u,v,y,w,sim) #50% of the time is spent here..
      #project gradient onto tangent plane
      h_u = g_u - scale(u,vec(sum(g_u.*u,1)))
      h_v = g_v - scale(v,vec(sum(g_v.*v,1)))
      #normalize
      norm_u = vec(sqrt(sumabs2(h_u,1)))
      norm_v = vec(sqrt(sumabs2(h_v,1)))
      scale!(h_u,1.0./norm_u)
      scale!(h_v,1.0./ norm_v)
      #BTLS
      while(true)
        if alpha < 1E-6
          break;
        end
        sw = sum(w)
        ca = cos(alpha)
        sa = sin(alpha)
        u_n = scale(u,ca) - scale(h_u,sa)
        v_n = scale(v,ca) - scale(h_v,sa)
        f_n = f_only_c(u_n,v_n,y,w,sim)
        #check for sufficient decrease
        if f_n <= f + c_1*dot(vec(u_n - u),vec(g_u)) + c_1*dot(vec(v_n -v),vec(g_v))
          break;
        end
        alpha = beta*alpha
      end
      if (alpha < 1E-6)
        println("small stepsize, breaking!")
        break;
      end
      u = u*cos(alpha) - h_u*sin(alpha)
      v = v*cos(alpha) - h_v*sin(alpha)
    end
  return [u;v]
end

function runForwardSimulator(s :: MatrixCompletion, parameters)
  u = parameters[1:s.m, :]
  v = parameters[s.m+1:end, :]
  rank = size(parameters,2)
  result = zeros(s.nnz,rank)
  # out_v = zeros(s.nnz)
  n = length(s.nz_i)
  for idx = 1:n
    @fastmath @inbounds @simd  for r = 1:rank
      result[idx,r] = u[s.nz_i[idx],r]*v[s.nz_j[idx],r]
    end
  end
  return result
end

function phi(s :: MatrixCompletion, parameters :: Matrix{Float64}, w :: Vector{Float64})
  u = parameters[1:s.m, :]
  v = parameters[s.m+1:end, :]
  rank = size(parameters,2)
  result = zeros(s.nnz)
  n = length(s.nz_i)
  for idx = 1:n
    @fastmath @inbounds  @simd for r = 1:rank
       result[idx] += w[r]*u[s.nz_i[idx],r]*v[s.nz_j[idx],r]
    end
  end
  return result
end

function lmo(s :: MatrixCompletion,vector :: Vector{Float64})
  M = sparse(s.nz_i,s.nz_j, vector)
  (u,sv,v,t1,t2,t3,t4) = svds(M; nsv = 5, tol = 1E-7)
  p = vec([-u[:,1];v[:,1]])
  o = runForwardSimulator(s,p)
  return (p, dot(vec(vector),vec(o)) )
end

function fg_c(u,v,y,w,simulator)
  rank = length(w)
  u=u'
  v=v'
  u_aa = Array{Array{Float64}}(simulator.m)
  v_aa = Array{Array{Float64}}(simulator.k)
  for i = 1:simulator.m
    u_aa[i] = copy(vec(u[:,i]))
  end
  for i = 1:simulator.k
    v_aa[i] = copy(vec(v[:,i]))
  end
  u = u_aa
  v = v_aa
  up = deepcopy(u)
  vp = deepcopy(v)
  f_val = ccall((:fg,"c/libmc.so.1.0"),Float64,(Ptr{Ptr{Float64}},Ptr{Ptr{Float64}},Ptr{Ptr{Float64}},Ptr{Ptr{Float64}},Ptr{Int64},Ptr{Int64},Ptr{Float64},Ptr{Float64},Int64,Int64),up,vp,u,v,sim.nz_i,sim.nz_j,w,y,length(sim.nz_j),rank)
  u = zeros(rank,simulator.m)
  v = zeros(rank,simulator.k)
  for i = 1:simulator.m
    u[:,i] = up[i]
  end
  for i = 1:simulator.k
    v[:,i] = vp[i]
  end
  return f_val, u',v'
end

function f_only_c(u, v,y,w ::Vector{Float64},simulator)
  rank = length(w)
  u=u'
  v=v'
  u_aa = Array{Array{Float64}}(simulator.m)
  v_aa = Array{Array{Float64}}(simulator.k)
  for i = 1:simulator.m
    u_aa[i] = copy(vec(u[:,i]))
  end
  for i = 1:simulator.k
    v_aa[i] = copy(vec(v[:,i]))
  end
  u = u_aa
  v = v_aa
  return ccall((:f_only,"c/libmc.so.1.0"),Float64,(Ptr{Ptr{Float64}},Ptr{Ptr{Float64}},Ptr{Int64},Ptr{Int64},Ptr{Float64},Ptr{Float64},Int64,Int64),u',v',sim.nz_i,sim.nz_j,w,y,length(sim.nz_j),rank)
end


text:
"""
    fgsm_s(R, chars)

Compute the Sᵢⱼ's in FGSM, without taking the product over the choices of the signs of the square roots. Return the smallest absolute value of the FGSM relations.
"""
function fgsm_s(R, chars)
    z = zeros(5);
    theta_constants = [[theta(z, R, char=c) for c in chars_s] for chars_s in chars];
    p = [sqrt(prod(theta_constant)) for theta_constant in theta_constants];
    S = collect(Iterators.flatten([p[1] + i2*p[2] + i3*p[3] + i4*p[4] + i5*p[5] + i6*p[6] for i2 in [1,-1], i3 in [1,-1], i4 in [1,-1], i5 in [1,-1], i6 in [1,-1]]));
    min_fgsm = minimum(abs.(S));
    return min_fgsm;
end

"""
    fgsm(τ, chars=fgsm_chars())

Compute the FGSM relations in genus 5. Return the largest absolute value of the relation.
"""
function fgsm(τ::Array{<:Number}, chars=fgsm_chars())
    R = RiemannMatrix(τ, nderivs=0);
    return fgsm(R, chars);
end

"""
    fgsm(R, chars=fgsm_chars())

Compute the FGSM relations in genus 5. Return the largest absolute value of the relation.
"""
function fgsm(R::RiemannMatrix, chars=fgsm_chars())
    max_fgsm = maximum([fgsm_s(R, c) for c in chars]);
    return max_fgsm;
end

"""
    random_nonfgsm(tol=0.1, trials=100)

Find a random genus 5 matrix in the Siegel upper half space which is not in the FGSM locus, such that the largest FGSM relation has absolute value at least tol, using input number of trials.
"""
function random_nonfgsm(tol::Real=0.1, trials::Integer=100)
    chars = fgsm_chars();
    t = 0; # stores value of largest FGSM relation encountered
    i = 0; # counter for number of trials
    max_matrix = rand(5,5); # stores candidate matrix
    while t < tol && i < trials
        τ = random_siegel(5);
        s = fgsm(τ, chars);
        if s > t
            t = s;
            max_matrix = τ;
        end
        i += 1;
    end
    return [max_matrix, t];
end

"""
    fgsm_chars()

Compute the characteristics used in the FGSM relations in genus 5.
"""
function fgsm_chars()
    chars_34 = [[[[0,0,0,0,0], [0,0,0,0,0]],
                 [[0,0,0,0,0], [1,1,1,1,1]],
                 [[0,0,1,1,0], [0,1,0,0,0]],
                 [[0,0,1,1,0], [1,0,1,1,1]],
                 [[1,1,0,0,0], [0,0,0,1,0]],
                 [[1,1,0,0,0], [1,1,1,0,1]],
                 [[1,1,1,1,0], [0,1,0,1,0]],
                 [[1,1,1,1,0], [1,0,1,0,1]]],
                [[[1,0,1,0,0], [0,0,0,0,0]],
                 [[1,0,1,0,0], [1,1,1,1,1]],
                 [[1,0,0,1,0], [0,1,0,0,0]],
                 [[1,0,0,1,0], [1,0,1,1,1]],
                 [[0,1,1,0,0], [0,0,0,1,0]],
                 [[0,1,1,0,0], [1,1,1,0,1]],
                 [[0,1,0,1,0], [0,1,0,1,0]],
                 [[0,1,0,1,0], [1,0,1,0,1]]],
                [[[0,0,0,0,0], [0,0,1,1,0]],
                 [[0,0,0,0,0], [1,1,0,0,1]],
                 [[0,0,1,1,0], [0,1,1,1,0]],
                 [[0,0,1,1,0], [1,0,0,0,1]],
                 [[1,1,0,0,0], [0,0,1,0,0]],
                 [[1,1,0,0,0], [1,1,0,1,1]],
                 [[1,1,1,1,0], [0,1,1,0,0]],
                 [[1,1,1,1,0], [1,0,0,1,1]]],
                [[[1,0,0,0,1], [0,0,0,0,0]],
                 [[1,0,0,0,1], [1,1,1,1,1]],
                 [[1,0,1,1,1], [0,1,0,0,0]],
                 [[1,0,1,1,1], [1,0,1,1,1]],
                 [[0,1,0,0,1], [0,0,0,1,0]],
                 [[0,1,0,0,1], [1,1,1,0,1]],
                 [[0,1,1,1,1], [0,1,0,1,0]],
                 [[0,1,1,1,1], [1,0,1,0,1]]],
                [[[0,0,1,0,1], [0,0,0,0,0]],
                 [[0,0,1,0,1], [1,1,1,1,1]],
                 [[0,0,0,1,1], [0,1,0,0,0]],
                 [[0,0,0,1,1], [1,0,1,1,1]],
                 [[1,1,1,0,1], [0,0,0,1,0]],
                 [[1,1,1,0,1], [1,1,1,0,1]],
                 [[1,1,0,1,1], [0,1,0,1,0]],
                 [[1,1,0,1,1], [1,0,1,0,1]]],
                [[[1,0,0,0,1], [0,0,1,1,0]],
                 [[1,0,0,0,1], [1,1,0,0,1]],
                 [[1,0,1,1,1], [0,1,1,1,0]],
                 [[1,0,1,1,1], [1,0,0,0,1]],
                 [[0,1,0,0,1], [0,0,1,0,0]],
                 [[0,1,0,0,1], [1,1,0,1,1]],
                 [[0,1,1,1,1], [0,1,1,0,0]],
                 [[0,1,1,1,1], [1,0,0,1,1]]]];
    chars_35 = [[[[0,0,0,0,0], [0,0,0,0,0]],
                 [[0,0,0,0,0], [1,1,1,1,1]],
                 [[0,0,1,0,1], [0,1,0,0,0]],
                 [[0,0,1,0,1], [1,0,1,1,1]],
                 [[1,1,0,0,0], [0,0,0,0,1]],
                 [[1,1,0,0,0], [1,1,1,1,0]],
                 [[1,1,1,0,1], [0,1,0,0,1]],
                 [[1,1,1,0,1], [1,0,1,1,0]]],
                [[[1,0,1,0,0], [0,0,0,0,0]],
                 [[1,0,1,0,0], [1,1,1,1,1]],
                 [[1,0,0,0,1], [0,1,0,0,0]],
                 [[1,0,0,0,1], [1,0,1,1,1]],
                 [[0,1,1,0,0], [0,0,0,0,1]],
                 [[0,1,1,0,0], [1,1,1,1,0]],
                 [[0,1,0,0,1], [0,1,0,0,1]],
                 [[0,1,0,0,1], [1,0,1,1,0]]],
                [[[0,0,0,0,0], [0,0,1,0,1]],
                 [[0,0,0,0,0], [1,1,0,1,0]],
                 [[0,0,1,0,1], [0,1,1,0,1]],
                 [[0,0,1,0,1], [1,0,0,1,0]],
                 [[1,1,0,0,0], [0,0,1,0,0]],
                 [[1,1,0,0,0], [1,1,0,1,1]],
                 [[1,1,1,0,1], [0,1,1,0,0]],
                 [[1,1,1,0,1], [1,0,0,1,1]]],
                [[[1,0,0,1,0], [0,0,0,0,0]],
                 [[1,0,0,1,0], [1,1,1,1,1]],
                 [[1,0,1,1,1], [0,1,0,0,0]],
                 [[1,0,1,1,1], [1,0,1,1,1]],
                 [[0,1,0,1,0], [0,0,0,0,1]],
                 [[0,1,0,1,0], [1,1,1,1,0]],
                 [[0,1,1,1,1], [0,1,0,0,1]],
                 [[0,1,1,1,1], [1,0,1,1,0]]],
                [[[0,0,1,1,0], [0,0,0,0,0]],
                 [[0,0,1,1,0], [1,1,1,1,1]],
                 [[0,0,0,1,1], [0,1,0,0,0]],
                 [[0,0,0,1,1], [1,0,1,1,1]],
                 [[1,1,1,1,0], [0,0,0,0,1]],
                 [[1,1,1,1,0], [1,1,1,1,0]],
                 [[1,1,0,1,1], [0,1,0,0,1]],
                 [[1,1,0,1,1], [1,0,1,1,0]]],
                [[[1,0,0,1,0], [0,0,1,0,1]],
                 [[1,0,0,1,0], [1,1,0,1,0]],
                 [[1,0,1,1,1], [0,1,1,0,1]],
                 [[1,0,1,1,1], [1,0,0,1,0]],
                 [[0,1,0,1,0], [0,0,1,0,0]],
                 [[0,1,0,1,0], [1,1,0,1,1]],
                 [[0,1,1,1,1], [0,1,1,0,0]],
                 [[0,1,1,1,1], [1,0,0,1,1]]]];
    chars_45 = [[[[0,0,0,0,0], [0,0,0,0,0]],
                 [[0,0,0,0,0], [1,1,1,1,1]],
                 [[0,0,0,1,1], [0,1,0,0,0]],
                 [[0,0,0,1,1], [1,0,1,1,1]],
                 [[1,1,0,0,0], [0,0,0,0,1]],
                 [[1,1,0,0,0], [1,1,1,1,0]],
                 [[1,1,0,1,1], [0,1,0,0,1]],
                 [[1,1,0,1,1], [1,0,1,1,0]]],
                [[[1,0,0,1,0], [0,0,0,0,0]],
                 [[1,0,0,1,0], [1,1,1,1,1]],
                 [[1,0,0,0,1], [0,1,0,0,0]],
                 [[1,0,0,0,1], [1,0,1,1,1]],
                 [[0,1,0,1,0], [0,0,0,0,1]],
                 [[0,1,0,1,0], [1,1,1,1,0]],
                 [[0,1,0,0,1], [0,1,0,0,1]],
                 [[0,1,0,0,1], [1,0,1,1,0]]],
                [[[0,0,0,0,0], [0,0,0,1,1]],
                 [[0,0,0,0,0], [1,1,1,0,0]],
                 [[0,0,0,1,1], [0,1,0,1,1]],
                 [[0,0,0,1,1], [1,0,1,0,0]],
                 [[1,1,0,0,0], [0,0,0,1,0]],
                 [[1,1,0,0,0], [1,1,1,0,1]],
                 [[1,1,0,1,1], [0,1,0,1,0]],
                 [[1,1,0,1,1], [1,0,1,0,1]]],
                [[[1,0,1,0,0], [0,0,0,0,0]],
                 [[1,0,1,0,0], [1,1,1,1,1]],
                 [[1,0,1,1,1], [0,1,0,0,0]],
                 [[1,0,1,1,1], [1,0,1,1,1]],
                 [[0,1,1,0,0], [0,0,0,0,1]],
                 [[0,1,1,0,0], [1,1,1,1,0]],
                 [[0,1,1,1,1], [0,1,0,0,1]],
                 [[0,1,1,1,1], [1,0,1,1,0]]],
                [[[0,0,1,1,0], [0,0,0,0,0]],
                 [[0,0,1,1,0], [1,1,1,1,1]],
                 [[0,0,1,0,1], [0,1,0,0,0]],
                 [[0,0,1,0,1], [1,0,1,1,1]],
                 [[1,1,1,1,0], [0,0,0,0,1]],
                 [[1,1,1,1,0], [1,1,1,1,0]],
                 [[1,1,1,0,1], [0,1,0,0,1]],
                 [[1,1,1,0,1], [1,0,1,1,0]]],
                [[[1,0,1,0,0], [0,0,0,1,1]],
                 [[1,0,1,0,0], [1,1,1,0,0]],
                 [[1,0,1,1,1], [0,1,0,1,1]],
                 [[1,0,1,1,1], [1,0,1,0,0]],
                 [[0,1,1,0,0], [0,0,0,1,0]],
                 [[0,1,1,0,0], [1,1,1,0,1]],
                 [[0,1,1,1,1], [0,1,0,1,0]],
                 [[0,1,1,1,1], [1,0,1,0,1]]]];
    return [chars_34, chars_35, chars_45];
end


text:
Pkg.add("Clustering")

Pkg.update()

using Clustering
using RDatasets

xclara = dataset("cluster", "xclara")

names!(xclara, [symbol(i) for i in ["x", "y"]])

using Clustering
using Gadfly
iris = dataset("datasets", "iris")
features = array(iris[:, 1:4])'
# group the data onto 3 clusters
result = kmeans( features, 3 )
plot(iris, x = "PetalLength", y = "PetalWidth",
        color = result.assignments, Geom.point)

max_stars_repo_path: Chapter 7/B05321_CH07_Code01.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 7/B05321_CH07_Code01.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 7/B05321_CH07_Code01.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
__precompile__(true)

module AsyPlots

import Random,
       ColorTypes,
       Colors,
       Requires,
       ProgressMeter,
       Contour,
       Interpolations,
       DelimitedFiles,
       Statistics

import DataStructures: OrderedDict

import Base: show, print, string, convert,
             complex, abs, iterate, zero,
             +, *, -, /, ==

import LinearAlgebra: dot

export Point, Point2D, Point3D, PointShape, PointCloud,
       Path, Path2D, Path3D,
       Circle2D, Circle,
       Polygon2D, Polygon3D, Polygon, box,
       Label2D, Label3D, Label,
       PixelMap, heatmap,
       RawString,
       Surface,
       Plot, Plot2D, Plot3D,
       GraphicElement, GraphicElement2D, GraphicElement3D,
       NamedColor, Pen, NoPen, Arrow, Arrow3, NoArrow, Ticks,
       Shift, Shift2D, Scale, Scale2D, layout,
       defaultcolor, isdefault,
       isinside, iswellinside,
       asy, cairo,
       plot,
       interactive,
       isolines,
       save,
       junosvg, junopng, junodefault,
       fitwidth, fitheight,
       animate,
       piechart, histogram

abstract type Plot end
abstract type GraphicElement end
abstract type GraphicString <: GraphicElement end

struct AsyString
    str::String
    data::Dict
end

_BACKEND = "asy"
_DEFAULT_WIDTH = 256

function defaultwidth(n::Integer)
    global _DEFAULT_WIDTH
    _DEFAULT_WIDTH = n
    "Default image width set to $n"
end

"""
    asy()

Set the backend to Asymptote
"""
function asy()
    global _BACKEND
    _BACKEND = "asy"
    "Backend set to Asymptote"
end

"""
    cairo()

Set the backend to Cairo
"""
function cairo()
    global _BACKEND
    _BACKEND = "cairo"
    "Backend set to Cairo"
end

_HAVE_ASY = try success(`asy -version`); catch; false; end

if !_HAVE_ASY
    println("asy not found, using Cairo backend")
    cairo()
end

AsyString(s::AbstractString) = AsyString(s,Dict())

RealOrComplex = Union{Real,Complex}

include("util.jl")
include("pens_and_arrows.jl")
include("asy2d.jl")
include("asy3d.jl")
include("cairo.jl")
include("output.jl")
include("transformations.jl")
include("tools.jl")
include("animate.jl")
include("stats.jl")
include("contours.jl")
include("init.jl")

end # module


text:
#!/home/jburdick/bin/julia
# Writes out the coefficients of the bound matrix, both
# using exact binomials, and using various approximations
# (and log-transformed numbers), as a check of the math.

push!(LOAD_PATH, ".")

using CountingBound, LogCountingBound

"""
  Writes counts for some values of k, r, and n.
  k: the value of k
  maxN: the maximum value of n (this will use
		powers of two	<= maxN)
	logADropOff: this will start r at k, and continue as the log bound
		increases, until the bound drops below (max - logADropOff)
  outputDir: directory in which to write output files
  Side effects: writes files A and b.
    (The choices of r are a bit arbitrary).
"""
function writeCounts(k, maxN, logADropOff, outputDir)
  # values of n to use, from k to maxN
  # (for now, spaced at powers of 2)
  nList = [BigInt(2^i) for i in 1:trunc(log(2,maxN))]
  # print(nList)
  # print("\n")
  # print(typeof(nList[1]))
  # print("\n")
  # ??? for some reason, this was giving a syntax error
  # nList = [n for n in nList if n>=k && n<=maxN]
  # therefore, using filter()...
  nList = filter(n -> n>=k && n<=maxN, nList)

  # FIXME should create output directory
  # first, write the bound b
  of = open(outputDir * "/b_k=" * string(k) * "_maxN=" * string(maxN) * ".csv", "w")
  write(of, "k,n,logBound\n")
  for n in nList
    bound = countingBound(binomial(n, 2), binomial(n, k))
		logBound = NaN
		if bound > 0
			logBound = log(bound)
		end
    write(of, string(k) * "," * string(n) * ","
      * string(logBound) * "\n")
  end
  close(of)

  # then, write the coefficients
  of = open(outputDir * "/A_k=" * string(k) * "_maxN=" * string(maxN) * ".csv", "w")
  write(of, "k,r,n,logA\n")
  for n = nList
		# tracks the maximum (initially undefined)
		maxLogA = NaN
		# loop through possible values of r (this will usually stop early)
    for r = k:n
      print("k=" * string(k) * " r=" * string(r) * " n=" * string(n) * "\n")
			# compute log of number of cliques, and print it
      logA = logNumMaximalCliques(k, r, n)
      write(of, string(k) * "," * string(r) * "," * string(n)
        * "," * string(logA) * "\n")
			# keep track of maximum
			if isnan(maxLogA) || logA > maxLogA
				maxLogA = logA
			end
			# if we've fallen far enough below the maximum, stop the loop
			if logA <= maxLogA - logADropOff
				break
			end
    end
  end
  close(of)

end

outputDir = "logCoef1"
mkpath(outputDir)

for k in [6,8,10,12,14,16,18,20,50,100,200,500,1000]
	writeCounts(k, 1e60, 100 * log(10), outputDir)
end

max_stars_repo_name: joshtburdick/misc

max_issues_repo_name: joshtburdick/misc

max_forks_repo_name: joshtburdick/misc


text:
"The most basic constructor"
mutable struct SingleCellExpr
    data::SparseMatrixCSC
    samples::DataFrame
    annotation::DataFrame
    name::String
    raw::Any
    zscored::Any
    outpath::String
end

"Read the data from a cellexalvrR database"
function BioData( file::String )
      if ! isfile(file)
        throw(ArgumentError("database file does not exists"))
      end

      db = SQLite.DB(file)
      so = SQLite.query(db, "SELECT  gene_id, cell_id, value from datavalues")
      #println(so)
      if typeof(so[2]) == DataType(Array{Union{Int64, Missings.Missing},1})
        A = sparse( so[1], so[2], so[3] );
      else
        A = sparse( map( Int64,collect(so[1]) ),  map(Int64,collect(so[2])) , map( Int64, collect(so[3]))  )
      end
      samples = SQLite.query(db, "SELECT * from cells")
      #println(samples)
      annotation = SQLite.query(db, "SELECT * from genes")
      #println(annotation)
      return ( SingleCellExpr( A, samples, annotation, "Unnamed SingleCellData", 0 , 0, pwd() ) )
end


text:
# benchmarks inner loop acquisition functions
using LinearAlgebra
using GaussianDistributions
using SARA
using SARA: random_sampling, uncertainty_sampling, ucb_sampling,
            ucb_acquisition,
            stripe_uncertainty_sampling, stripe_ucb_sampling,
            TemperatureProfile, stripe_to_global, get_relevant_T
using ForwardDiff
using Base.Threads: @threads
using SARA: oSARA, likelihood_ratio

############### generate sampling paths with different policies
function sampling_path(policy, G::Gaussian, x::AbstractMatrix, y::AbstractMatrix,
                       σ²::AbstractMatrix, stripe_conditions::AbstractVector, nsample::Int)
    C = G
    rmse, r2 = zeros(nsample), zeros(nsample) # rmse for data and estimated gradients
    indices = zeros(Int, nsample)
    helper(x, i) = @views vec(x[:, i]) # selects the columns indexed by i, and vectorizes the result
    for i in 1:nsample
        println("sample $i")
        xj = policy(C, stripe_conditions)
        j = findfirst(==(xj), stripe_conditions)
        indices[i] = j
        indices_i = @view indices[1:i]
        sort!(indices_i)
        xi, yi, σ²i = helper(x, indices_i), helper(y, indices_i), helper(σ², indices_i)
        C = conditional(G, xi, Gaussian(yi, σ²i), tol = 1e-6)
        μ = mean(C, vec(x))
        δ = vec(y) - μ
        rmse[i] = sqrt(mean(δ.^2 ./ vec(σ²))) # weighted root mean squared deviation
        r2[i] = likelihood_ratio(vec(y), μ, vec(σ²))
    end
    return rmse, r2
end

function outer_acquisition_benchmark(x::AbstractMatrix, y::AbstractMatrix,
                                σ²::AbstractMatrix, T_offset, stripe_conditions)
    subsampling = 2 # subsample observations per stripe (for quick testing only)
    stride = 2
    samplesperstripe = 16
    subsample(x::AbstractMatrix) = @view x[end+1-stride*samplesperstripe:stride:end, :]
    subsample(x::AbstractVector) = @view x[end+1-stride*samplesperstripe:stride:end]
    x = subsample(x)
    y = subsample(y)
    σ² = subsample(σ²)

    # standardize global data
    y = y .- mean(y) # otherwise r2 scores are initially negative
    y_normalization = std(y)
    y /= y_normalization
    σ² /= y_normalization^2 # noise variance

    tl = 30
    dl = .4
    k_scale = 1/4
    # kernel = k_scale^2 * oSARA([tl, dl]) # performance bug, sidestep with the following
    kernel = oSARA([tl, dl])
    σ² ./= k_scale^2

    # experimental parameters
    nstripes = 48 # how many stripes to acquire using active learning

    # stripe sampling acquisition functions
    constant_offset_bool = true
    nout = 32 # number of points on temperature curve that are observed per stripe
    constant_offset = Val(constant_offset_bool)
    # T_offset = passed as input # number of degrees from T_max we are generating data for
    T_range = range(T_offset[1], T_offset[2], length = nout)
    T_range = subsample(T_range) # [end-stride*samplesperstripe+1 : stride : end]
    T_offset = [maximum(T_range), minimum(T_range)]
    c_min, c_max = T_offset
    nout = length(T_range)
    relevant_T = get_relevant_T(constant_offset, c_min, c_max, nout)

    # random acquisition
    nrandom = 2
    random_funs = [random_sampling for _ in 1:nrandom]

    # stripe uncertainty sampling
    stripe_uncertainty = stripe_uncertainty_sampling(relevant_T)
    # stripe upper confidence bound sampling
    α_range = [1., 5.] # 10., 25., 50.] # coefficient of standard deviation in acquisition function
    ucb_funs = [ucb_sampling(α) for α in α_range]
    stripe_ucb_funs = [stripe_ucb_sampling(relevant_T, α) for α in α_range]

    policies = vcat(random_funs,
                    uncertainty_sampling, stripe_uncertainty,
                    ucb_funs, stripe_ucb_funs)

    npolicies = length(policies)
    rmse = zeros(nstripes, npolicies)
    r2 = zeros(nstripes, npolicies)
    G = Gaussian(kernel)
    @threads for i in eachindex(policies)
        println("policy $i of $(length(policies))")
        policy = policies[i]
        rmse[:, i], r2[:, i] = sampling_path(policy, G, x, y, σ², stripe_conditions, nstripes)
    end

    rmse_random_mean = vec(mean(rmse[:, 1:nrandom], dims = 2)) # averaging
    rmse_random_std = vec(std(rmse[:, 1:nrandom], dims = 2))
    rmse = hcat(rmse_random_mean, rmse[:, nrandom+1:end])

    r2_random_mean = vec(mean(r2[:, 1:nrandom], dims = 2))
    r2_random_std = vec(std(r2[:, 1:nrandom], dims = 2))
    r2 = hcat(r2_random_mean, r2[:, nrandom+1:end])

    if true
        f = h5open("outer_acquisition_benchmark.h5", "w")
        f["rmse"] = rmse
        f["r2"] = r2
        f["nstripes"] = nstripes
        f["kernel"] = ["oSARA"]
        f["temperature lengthscale"] = tl
        f["dwelltime lengthscale"] = dl
        f["npolicies"] = npolicies
        f["nrandom"] = nrandom
        f["ucb alpha range"] = α_range
        f["y_normalization"] = y_normalization
        f["metrics"] = ["rmse", "r2"]
        f["stride"] = stride
        f["samples per stripe"] = samplesperstripe
        f["subsampling"] = subsampling
        f["policy names"] = ["random", "uncertainty", "stripe uncertainty", "ucb", "stripe ucb"]
        close(f)
    end
    return rmse, r2
end

# load outer loop data
using HDF5
path = "SARA/ScienceAdvances2021/outer/data/" # put path to outer gradient data
# path = "ScienceAdvances2021/outer/"
# file = "Bi2O3_19F44_01_outer_loop_data.h5" # without input noise correction
file = "Bi2O3_19F44_01_outer_optical_gradients_l_2_input_noise_50.h5" # with input noise correction

f = h5open(path * file, "r")
dwelltimes = read(f, "dwelltimes")
temperatures = read(f, "temperatures")
gradients = read(f, "gradients")
var_gradients = read(f, "var(gradients)")
T_offset = read(f, "T_offset")
close(f)

# normalize global data
y_normalization = std(gradients)
gradients ./= y_normalization
var_gradients ./= y_normalization^2 # noise variance

x = tuple.(temperatures, dwelltimes)
y = gradients
σ² = var_gradients

# since the last condition of each column in the x array is offset from
# T_max by T_offset[2], need to add T_offset[2] back in to get the
# stripe condition
temperature_correction(x::NTuple{2}) = (x[1] + T_offset[2], x[2])
stripe_conditions = temperature_correction.(x[end, :])

rmse, r2 = outer_acquisition_benchmark(x, y, σ², T_offset, stripe_conditions)

r2_random = r2[:, 1]
r2_sigma = r2[:, 2]
r2_stripe_sigma = r2[:, 3]
nucb = (size(r2, 2)-3) ÷ 2
r2_ucb = r2[:, 4:3+nucb]
r2_stripe_ucb = r2[:, 4+nucb:end]

using Plots
plotly()
# names = ["random", "sigma", "ucb", "stripe sigma", "stripe ucb"]
plot(legend = :bottomright, xlabel = "number of stripes", ylabel = "R2")
plot!(r2_random, label = "random")
plot!(r2_sigma, label = "sigma")
plot!(r2_stripe_sigma, label = "stripe sigma")
plot!(r2_ucb, label = "ucb")
plot!(r2_stripe_ucb, label = "stripe ucb")
gui()

using HDF5
using Statistics
f = h5open("SARA/ScienceAdvances2021/outer/results/outer_acquisition_benchmark.h5")
r2 = read(f, "r2")
nrandom = read(f, "nrandom")
random = r2[:, 1:nrandom]

m = vec(mean(random, dims = 2))
s = vec(std(random, dims = 2))
threshold(p) = findfirst(>(p), m)
println(threshold(.7))
close(f)

f = h5open("SARA/ScienceAdvances2021/outer/results/outer_acquisition_benchmark_input_noise_0.h5")
r2 = read(f, "r2")
nrandom = read(f, "nrandom")
random = r2[:, 1:nrandom]

m = vec(mean(random, dims = 2))
s = vec(std(random, dims = 2))
threshold(p) = findfirst(>(p), m)
println(threshold(.7))
close(f)


text:
import Distributions: _rand!

function jac_inv_transform(dist::Distribution, x::T where T<:Real)
    ForwardDiff.derivative(x -> invlink(dist, x), x)
end

function jac_inv_transform(dist::Distribution, x::AbstractArray{T} where T <: Real)
    ForwardDiff.jacobian(x -> invlink(dist, x), x)
end

function jac_inv_transform(dist::Distribution, x::TrackedArray{T} where T <: Real)
    Tracker.jacobian(x -> invlink(dist, x), x)
end

# instead of creating a diagonal matrix, we just do elementwise multiplication
center_diag_gaussian(x, μ, σ) = (x .- μ) ./ σ
center_diag_gaussian_inv(η, μ, σ) = (η .* σ) .+ μ

# Mean-field approximation used by ADVI
struct MeanField{TDists, V} <: VariationalPosterior where {V <: AbstractVector{<: Real}, TDists <: AbstractVector{<: Distribution}}
    μ::V
    ω::V
    dists::TDists
    ranges::Vector{UnitRange{Int}}
end

Base.length(advi::MeanField) = length(advi.μ)

function _rand!(
    rng::AbstractRNG,
    q::MeanField{TDists},
    x::AbstractVector{T}
) where {T<:Real, TDists <: AbstractVector{<: Distribution}}
    # extract parameters for convenience
    μ, ω = q.μ, q.ω
    num_params = length(q)

    for i = 1:size(q.dists, 1)
        prior = q.dists[i]
        r = q.ranges[i]

        # initials
        μ_i = μ[r]
        ω_i = ω[r]

        # # sample from VI posterior
        η = randn(rng, length(μ_i))
        ζ = center_diag_gaussian_inv(η, μ_i, exp.(ω_i))
        θ = invlink(prior, ζ)

        x[r] = θ
    end

    return x
end


text:
# Unlike normal Signals, RaggedSignals do *not* have a common timebase;
# their timebase may vary from vector to vector. They are predominantly here to
# support windowed repetitions of irregular signals (e.g., spike rasters).
type RaggedSignal{T, S, N, M, D} <: AbstractSignal{T, S, N}
    time::Array{T, M}
    data::D

    dims::NTuple{N, Symbol}

    meta::Dict{Symbol, Any}
    RaggedSignal(t,d,n,m) = new(t,d,n,m)
end


text:
module EGI

using Polyhedra, QHull, Makie
export ExtendedGaussImage, Surface, Volume, Polyhedron, Normal, randPolyhedron, RenderPolyhedron, Little, SavePolyhedronArray

struct ExtendedGaussImage
    Normal::Array{Float64,2} where T <: Union{Int,Float64}
    Area::Union{Array{T,2},Array{T,1}} where T <: Union{Int,Float64}
    function ExtendedGaussImage((Normal, Area))
        if size(Normal)[2] != length(Area)
            error("Number of surface normal does not match number of surface area. ")
        end
        if sum([norm(Normal[:,i]) for i in 1:length(Area)] - ones(length(Area))) > 1e-10
            error("Normal not unit vector")
        end
        if sum(Normal*Area) > 1e-10
            error("It is not a valid convex polyhedron (Does not satisfy sum of surfaces normal times area equals 0")
        end
        new(Normal,Area)
    end
end

function randEGI(faces::Int64)
    for i = 1:(10*faces)
        global Normed = zeros(3,1)
        for j = 1:faces
            N = rand(3)*2 .-1
            N = N ./ norm(N)
            Normed = [Normed N]
        end
        global Normed = Normed[:,2:end]

        v2 = rand(faces - 3, 1)
        v1 = inv(Normed[:, 1:3])*(-Normed[:, 4:end]*v2)
        global v = [v1;v2]
        if (sum(v .> zeros(faces,1)) == faces)
            break
        end
    end
    return Normed,v
end

function Surface(p::QHull.Polyhedron)
    S = Float64[]
    for pi in eachindex(halfspaces(p))
        global n = get(p, pi).a
        n = n/norm(n)
        get(p, pi).β
        global pts = incidentpoints(p, pi)
        global b = pts[2] - pts[1]
        b = b/norm(b)
        global a = cross(n, b)
        norm(a)
        append!(S, Float64(Polyhedra.volume(polyhedron(vrep([[a b]'*pt for pt in pts]), QHull.Library()))))
    end
    return S
end

function Volume(H::Array{Float64,1}, S::Array{Float64,1})
    return 1/3*(H'*S)
end

function Volume(p::QHull.Polyhedron)
    faces = length(hrep(p).b)
    H = MixedMatHRep(hrep(p)).b ./[norm(MixedMatHRep(hrep(p)).A[i,:]) for i in 1:faces]
    S = Float64[]
    for pi in eachindex(halfspaces(p))
        global n = get(p, pi).a
        n = n/norm(n)
        get(p, pi).β
        global pts = incidentpoints(p, pi)
        global b = pts[2] - pts[1]
        b = b/norm(b)
        global a = cross(n, b)
        norm(a)
        append!(S, Float64(Polyhedra.volume(polyhedron(vrep([[a b]'*pt for pt in pts]), QHull.Library()))))
    end
    return 1/3*(H'*S)
end

function Polyhedron(H,N)
    h = hrep(N', [norm(N[:,i]) for i in 1:length(H)] .* H)
    p = polyhedron(h, QHull.Library())
    return p
end

function Normal(p::QHull.Polyhedron)
    faces = length(hrep(p).b)
    return [MixedMatHRep(hrep(p)).A[i,:]./norm(MixedMatHRep(hrep(p)).A[i,:]) for i in 1:faces]
end

function randPolyhedron(faces)
    egi = randEGI(faces)
    img = ExtendedGaussImage(egi)
    for i in 1:(10*faces)
        H = rand(faces)
        global p = Polyhedron(H, img.Normal)
        if length(Surface(p)) == faces
            break
        end
    end
    return p, ExtendedGaussImage((img.Normal, Surface(p)[:,:]))
end

function RenderPolyhedron(p::QHull.Polyhedron)
    m = Polyhedra.Mesh(p)
    scene = mesh(m, color=:white)
end

Base.length(q::QHull.Polyhedron) = 1
Base.iterate(q::QHull.Polyhedron) = (q, nothing)
function Little(k, ϵ, img::ExtendedGaussImage)
    faces = length(img.Area)
    H = ones(faces)
    p1 = img.Area / sum(img.Area)
    KLArray = []
    KL = 0.05
    qArray = QHull.Polyhedron[]
    q = Polyhedron(H, img.Normal)
    while KL >= ϵ
        H = H / (Volume(q))^(1/3)
        q = Polyhedron(H, img.Normal)
        ∇ = 1/3*Surface(q)
        Δ = ∇ - (img.Area * norm(∇))/ norm(img.Area)
        H = vec(H + k .* Δ)
        H = H / (Volume(q))^(1/3)
        q = Polyhedron(H, img.Normal)
        p2 = Surface(q) / sum(Surface(q))
        KL = sum([p2[i]*log(p2[i]/p1[i]) for i in 1:faces])
        @show p2 ./ p1
        @show KL
        append!(KLArray, KL)
        append!(qArray, q)
    end
    return qArray, KLArray
end

function SavePolyhedronArray(qArray::Array{QHull.Polyhedron,1})
    for i in 1:length(qArray)
        scene = Scene()
        m = Polyhedra.Mesh(qArray[i])
        mesh!(scene, m, color=:white)
        Makie.save(string(i)*".png", scene)
        sleep(1/24)
    end
end

end # module


text:
#Duffig Oscillator
#Load the toolbox
using FPTControl

# if  Adaptive parameter value is 1:it will be use Fixed Point Transformation Method, else (example 0) it will be just a simple PID Controller
Adaptive=1
#Simulation Lenght
LONG=Int(1e4)
# time step
δt=1e-3

#PID tune parameter
Λ=3

#Nominal Trajectory for a simple sin() type Nominal Trajectory
ω=5
Ampl=0.4

#Controller Parameters
K=1e5
B=-1
A=1e-5

# Exact model parameters
αₑ=1
δₑ=0.2
βₑ=1

# Approx. model parameters
αₐ=0.8
δₐ=0.1
βₐ=0.9

###############
# Exact Model #
###############

function exact(q,q_p,u)
  out=αₑ*q+δₑ*q_p+βₑ*q^3+u
  return out #q_pp
end

##########################
# Approx Model (inverse model)    #
##########################

function approx(q,q_p,q_pp)
  out=q_pp-αₐ*q-δₐ*q_p-βₐ*q^3
  return out #u
end

######################################
# Arrays for Plotting and simulation #
######################################

qN=zeros(LONG)
qN_p=zeros(LONG)
qN_pp=zeros(LONG)

qDes_pp=zeros(LONG)

qDef_pp=zeros(LONG)
u=zeros(LONG)

q=zeros(LONG)
q_p=zeros(LONG)
q_pp=zeros(LONG)

t=zeros(LONG)

#Initial condtions
error_int=0
past_input=0
past_response=0

#simulation

for i=1:LONG-1
 #store time
 t[i]=δt*i
 #Nominal Trajectory and the derivatives
 qN[i]=Ampl*sin(ω*t[i])
 qN_p[i]=Ampl*ω*cos(ω*t[i])
 qN_pp[i]=-Ampl*ω^2*sin(ω*t[i])
 #error (by definition)
 h=qN[i]-q[i]
 #first time derivative of the error
 h_p=qN_p[i]-q_p[i]
 #the error vector, for the Kinematic Block
 errors=[error_int,h,h_p]
 #the kinematic block, (Simple PID type controller for 2nd order system) produce the Desired Acceleration.
 qDes_pp[i]=KB(3,Λ,errors,qN_pp[i])
 # If we turn RFPT on, and from the 10th cycle (make sure, the buffers is full)
 if Adaptive==1 && i>10
   #Deformate the system
   qDef_pp[i]=G_SISO(past_input,past_response,qDes_pp[i],K,B,A)
 else
   #if RFPT is off, pass the desired value (thats why it is a simple PID controller)
   qDef_pp[i]=qDes_pp[i]
 end #if
 #store the output of G function, for the next cycle
 past_input=qDef_pp[i]
 # create the control signal
 u[i]=approx(q[i],q_p[i],qDef_pp[i])
 # get the realized response from our system
 q_pp[i]=exact(q[i],q_p[i],u[i])
 # store it for the next cycle (G function needs it.)
 past_response=q_pp[i]
 # Integral of the errror  for the next cycle
 error_int=Integ(error_int,δt,h)
 # Integrals to produce realized q_p and q
 q_p[i+1]=Integ(q_p[i],δt,q_pp[i])
 q[i+1]=Integ(q[i],δt,q_p[i])
end #for
#the last step
t[LONG]=δt*LONG
qN[LONG]=Ampl*sin(ω*t[LONG])
qN_p[LONG]=Ampl*ω*cos(ω*t[LONG])
qN_pp[LONG]=-Ampl*ω^2*sin(ω*t[LONG])
h=qN[LONG]-q[LONG]
h_p=qN_p[LONG]-q_p[LONG]
errors=[error_int,h,h_p]
qDes_pp[LONG]=KB(3,Λ,errors,qN_pp[LONG])
if Adaptive==1
  qDef_pp[LONG]=G_SISO(past_input,past_response,qDes_pp[LONG],K,B,A)
else
  qDef_pp[LONG]=qDes_pp[LONG]
end #if
u[LONG]=approx(q[LONG],q_p[LONG],qDef_pp[LONG])
q_pp[LONG]=exact(q[LONG],q_p[LONG],u[LONG])

#figures

using PyPlot
figure(1)
grid("on")
title("Nominal and realized traj. Nominal: red Realized: blue")
plot(t,qN,color="red")
plot(t,q,color="blue","r--")

figure(2)
grid("on")
title("the control signal")
plot(t,u,color="orange")

figure(3)
grid("on")
title("Tracking Error")
plot(t,qN-q,color="blue")

figure(4)
grid("on")
title("Trajectory Phase Space")
plot(qN,qN_p,color="red")
plot(q,q_p,color="blue","r--")

figure(5)
grid("on")
title("Accelerations: Nominal: red, Desired: orange, Realized: Blue")
plot(t,qN_pp,color="red")
plot(t,qDes_pp,color="orange","r-.")
plot(t,q_pp,color="blue","r--")

max_stars_repo_path: examples/RFPT_SISO/Duffing Oscillator with FPTControl Toolbox/Duffing_RFPT_with_FPTControl_v2.jl

max_stars_repo_name: krisztiankosi/Control

max_issues_repo_name: krisztiankosi/Control

max_forks_repo_name: krisztiankosi/Control


text:
function approximate_Floquet_multipliers(filename::AbstractString)
    x = load(filename)["x"]
    return _approximate_Floquet_multipliers(x)
end

function _approximate_Floquet_multipliers(x)
    x_₆₄ = complex.(Float64.(mid.(real.(x))), Float64.(mid.(imag.(x))))
    return eigvals(monodromy(x_₆₄); sortby = abs2)
end

function monodromy(p)
    M = zeros(5, 5)
    tspan = (0.0, 2π)
    for i ∈ 1:5
        u0 = I[1:5,i]
        prob = ODEProblem(vector_field!, u0, tspan, p)
        sol = solve(prob, Tsit5())
        M[:,i] = sol(2π)
    end
    return M
end

function vector_field!(dx, x, p, t)
    u = real.(coefficients(component(p, 1)(t)))
    γ = real(component(p, 2)[1])
    sqrt3 = sqrt(3.0)
    dx[1] = 2γ*(-(4u[3] - u[4] - u[5]) * x[3] + (u[3] + 2(u[4] - u[5])) * x[4] + (u[3] - 2(u[4] - u[5])) * x[5])
    dx[2] = 2γ*sqrt3*((u[5] - u[4]) * x[3] + (2u[4] - u[3]) * x[4] + (u[3] - 2u[5]) * x[5])
    dx[3] = 4γ*(u[3] * x[1] + u[1] * x[3])
    dx[4] = -2γ*(u[4] * x[1] + sqrt3 * u[4] * x[2] + (u[1] + sqrt3 * u[2]) * x[4])
    dx[5] = -2γ*(u[5] * x[1] - sqrt3 * u[5] * x[2] + (u[1] - sqrt3 * u[2]) * x[5])
end


text:
#!/usr/bin/env julia

module PE14

# collatz(n::Int) = iseven(n) ? collatz(n÷2)+1 : collatz(3n+1)+1

function collatz(n::Integer)
    local len=0;
    while n>1
        n = iseven(n) ? n÷2 : 3n+1;
        len += 1;
    end
    len
end

end

export pe14

pe14(n::Integer=10^6)=PE14.collatz.(1:n)|>findmax|>last

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe14()|>print
end

max_stars_repo_path: 14. Longest Collatz sequence/14.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 14. Longest Collatz sequence/14.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 14. Longest Collatz sequence/14.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
@testset "Sellmeier" begin
    @testset "String representation" begin
        Bstr = string(BK7.B)
        Cstr = join(string.(BK7.C), ", ")
        @test string(BK7) == "Medium($(Bstr), [$(Cstr)])"
    end

    @test BK7(0.5876u"μm") ≈ 1.5168 atol=1e-3
    @test SiO₂(0.5876u"μm") ≈ 1.4585 atol=1e-3

    @testset "Dispersion" begin
        # Dispersing a pulse through a positive amount of glass should
        # /delay/ the pulse, i.e. its maximum should arrive /later/,
        # and vice versa for a negative amount of glass (achievable
        # through precompensation, common in pulse compressors).

        λ = 500.0u"nm"
        f₀ = u"c"/λ |> u"THz"
        ω₀ = 2π*f₀
        τ = 6.2u"fs" # Pulse duration, intensity FWHM
        γ = τ^2/8log(2)

        f = range(0,stop=30,length=2000)*f₀
        ω = 2π*f

        Ê = exp.(-(ω .- ω₀).^2*γ)
        Ê′ = Ê.*dispersion(BK7, 6u"μm", f)
        Ê′′ = Ê.*dispersion(BK7, -6u"μm", f)
        Ê′′′ = Ê.*dispersion(BK7, -6u"μm", f, f₀)

        time_domain_envelope(spectrum) = abs.(fftshift(ifft(spectrum)*√(length(spectrum))))

        @test argmax(time_domain_envelope(Ê′)) > argmax(time_domain_envelope(Ê))
        @test argmax(time_domain_envelope(Ê′′)) < argmax(time_domain_envelope(Ê))
        @test argmax(time_domain_envelope(Ê′′′)) == argmax(time_domain_envelope(Ê))
    end
end


text:
using Images
using Clustering

function img_to_arr(img)
    return convert(Array{Float64}, channelview(img))
end

"""
    pixel(img, [, n_color, w, h]) -> Array{RGB{Float64}}

# Examples

```
julia> using PixelArt

julia> using Images

julia> img = load("img.jpg");

julia> img_pixel = pixel(img);

julia> save("img_pixel.jpg", img_pixel)
```
"""
function pixel(img::AbstractArray; n_color=5, w=64, h=64)
    img = imresize(img, (w, h))
    img = reshape(img_to_arr(img), (3, :))
    color_class = kmeans(img, n_color)
    img = hcat((x -> color_class.centers[:, x]).(color_class.assignments)...)
    return convert(Array{RGB{Float64}}, colorview(RGB, reshape(img, (3, w, h))))
end


text:
__precompile__()

module MultiScaleArrays

import Base: length, push!, deleteat!, getindex, setindex!, eachindex, ndims,
       size, print_matrix, similar, hcat, vcat,
       ==, *, +, /, -, show, vec, reshape

import RecursiveArrayTools: recursivecopy,recursivecopy!

import RecursiveArrayTools: chain

abstract type AbstractMultiScaleArray{B}     <: AbstractVector{B} end
abstract type AbstractMultiScaleArrayLeaf{B} <: AbstractMultiScaleArray{B} end
abstract type AbstractMultiScaleArrayHead{B} <: AbstractMultiScaleArray{B} end

using DiffEqBase, Statistics, LinearAlgebra, FiniteDiff
import OrdinaryDiffEq, StochasticDiffEq, ForwardDiff, SparseDiffTools

Base.show(io::IO, x::AbstractMultiScaleArray) = invoke(show, Tuple{IO, Any}, io, x)
Base.show(io::IO, ::MIME"text/plain", x::AbstractMultiScaleArray) = show(io, x)

include("shape_construction.jl")
include("addition_deletion.jl")
include("indexing.jl")
include("math.jl")
include("level_iterations.jl")
include("diffeq.jl")
include("print_human_readable.jl")

using TreeViews
TreeViews.hastreeview(x::AbstractMultiScaleArray) = true
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::AbstractMultiScaleArray) = x

# Types
export AbstractMultiScaleArray, AbstractMultiScaleArrayLeaf,
       AbstractMultiScaleArrayHead

# Constructors
export construct, recursivecopy!

# Addition Deletion
export add_node!, remove_node!

# Indexing
export num_nodes, getindices

# Misc
export LevelIterIdx, LevelIter, level_iter

# print_human_readable
export print_human_readable

end # module


text:
"""
Problem 1:

Given an input array of integers, count the number of values which increase relatively
to their precedent.

Problem 2:
Same, but using the convolution of the signal with a 3-position rectangular window.

"""

using DSP

# Read input
f = open("./inputs/day1.txt", "r")
data = map( x -> parse(Int,x), eachline(f))

# Problem
function countIncreasing(v::Vector{Int}) :: Int
  counter = 0
  for i in 2:length(v)
    if v[i] > v[i-1]
      counter += 1
    end
  end
  return counter
end

countIncreasingSums(v) = countIncreasing(conv(v,[1,1,1])[3:end-2])

println("Problem 1: ")
println(countIncreasing(data))
println("Problem 2: ")
println(countIncreasingSums(data))


text:
## Permutation Population Generators

# Kmedian epsilon with K minpts
function PGen_Kmedian{T <: Real}(A::Matrix{T},k::Int)
  ϵ = eps_Kmedian(A,k)
  return PGen_clust(A,ϵ,k)
end

# Kmean epsilon with K minpts
function PGen_Kmean{T <: Real}(A::Matrix{T},k::Int)
  ϵ = eps_Kmean(A,k)
  return PGen_clust(A,ϵ,k)
end

# Grid estimate (uniform) epsilon with K minpts
function PGen_Kgrid{T <: Real}(A::Matrix{T},k::Int)
  ϵ = eps_grid(A)
  return PGen_clust(A,ϵ,k)
end

# Sample random from K [median,85%] epsilon with K minpts
function PGen_Km85{T <: Real}(A::Matrix{T},k::Int)
  eps = eps_randm85(A,k)
  return PGen_clust(A,eps,k)
end

# Generic cluster generator function
function PGen_clust{T<:Real}(D::DenseMatrix{T}, ϵ::Real, minpts::Int)
  rslt = dbscan(D,ϵ,minpts)
  fillZero!(rslt.assignments)
  cmap = clustmap(rslt.assignments)
  k = length(keys(cmap))
  #solver for the high level clusters of GA
  #the function will eventually return a generator
  function solveTop(pfx::Function)
    tops = zeros(Int,k)
    bottoms = Vector{Vector{Int}}(k)
    #Extract a random cluster member to represent cluster
    for i in 1:k
      tops[i] = rand(cmap[i])
    end
    cord = pcGA(pfx,tops,false)
    #shortest path for each cluster
    for i in 1:k
      bottoms[i] = pcGA(pfx,cmap[i],true)
    end
    function clustPop(N::Int)
      tgt = zeros(Int,N)
      p = 1
      for i in cord
        l = length(bottoms[i])
        tgt[p:p+l-1] = circshift(bottoms[i],rand(1:l))
        p += l
      end
      return tgt
    end
    return clustPop
  end
  return solveTop
end


text:
using Quadrature
using Base.Test

xsamp = linspace(0.0, 10.0, 11)
ysamp = map(x->x^3, xsamp)
@test simps(ysamp, xsamp) == 2500.0

xsamp = linspace(0.0, 10.0, 1001)
ysamp = map(x->x^4, xsamp)
@test_approx_eq simps(ysamp, xsamp) 20000.0

# make sure we give errors for wrong input shape
# input ans samples are not the same length
@test_throws ArgumentError simps([0.5, 0.5, 0.3], [0.1, 0.2])
# simpsons rule needs odd input length
@test_throws ArgumentError simps([0.5, 0.5, 0.3, 0.2], [0.1, 0.2, 0.3, 0.4])

max_stars_repo_name: gabrielgellner/Quadrature

max_issues_repo_name: gabrielgellner/Quadrature

max_forks_repo_name: gabrielgellner/Quadrature


text:
using Distributions
using StaticArrays

struct DCHSBM_sampler{K, I <: Integer} <: Random.Sampler{NTuple{K, I}} # TODO make fully parameterized
    groups::Vector{AliasTable{I, UnitRange{I}}}
    ms::AliasTable{NTuple{K, I}, Vector{NTuple{K, I}}}#TODO try inline groups
end

"""
    DCHSBM_sampler(Z, θ, kmax)

Create a hypergraph sampler whose probability distribution follows the specified
degree corrected hyperstochastic block model.

Draw a hypergraph from the sampler with `rand(sampler)`.

# Arguments
- `Z::AbstractVector{<:Integer}` the group assignment of each node
- `θ::AbstractVector{<:Real}` the degree correction weight of each node.
- `kmax::Integer` is the largest possible hyperedge size the graph will generate
`Z` and `θ` must have the same length.

# Runtime

This method:
O(`tensor_cells * kmax + nodes`) ≈
400ns * tensor_cells * kmax + 30ns * nodes

Sampling:
O(`edges * kmax`) ≈
150ns * edges * kmax + 1300ns * edges
"""
DCHSBM_sampler(intensity::Function, Z::AbstractVector{<:Integer}, θ::AbstractVector{<:Real}, kmax::Integer) =
    DCHSBM_sampler(intensity, Z, θ, Val{kmax}())
function DCHSBM_sampler(intensity::Function, Z::AbstractVector{<:Integer}, θ::AbstractVector{<:Real}, kmax::Val{K}) where K
    #Note that sorting is O(n) here with a moderate constant factor, and sort checking is practically free
    @assert axes(Z) == axes(θ)
    if !issorted(Z)
        perm = sortperm(Z)
        Z = Z[perm]
        θ = θ[perm]
    end

    I = eltype(Z)
    groups = AliasTable{I, UnitRange{I}}[]
    group_sizes = I[]
    start = 1
    while start <= length(Z)
        next = start
        while next <= length(Z) && Z[next] == Z[start]; next += 1 end
        weights = θ[start:next-1]
        push!(groups, AliasTable(weights, start:next-1))
        push!(group_sizes, next-start)
        start = next
    end

    ms = with_replacement_combinations(length(groups), kmax)
    ms_weights = similar(ms, Float64)
    for (i, m) in enumerate(ms)
        a = intensity(m)
        b = multiply_by_cell_count(a, group_sizes, m)
        ms_weights[i] = b
    end
    expected_edges = sum(ms_weights)
    distribution_of_ms = AliasTable!(ms_weights, ms)

    DCHSBM_sampler(groups, distribution_of_ms)
end

function with_replacement_combinations(n::Integer, k::Val{K}) where K
    m = @MVector fill(1, K)
    out = Vector{NTuple{K, typeof(n)}}(undef, binomial(n+K-1, K))
    for i in eachindex(out)
        out[i] = Tuple(m)
        for j in eachindex(m)
            if j == K || m[j] < m[j+1]
                m[j] += 1
                break
            else
                m[j] = 1
            end
        end
    end
    out
end

# multiply_by_cell_count could be integrated into with_replacement_combinations to get
# asymtotic speedups for sparse graphs
function multiply_by_cell_count(x, group_sizes, m)
    i0 = firstindex(m)
    m0 = first(m)
    i = i0 + 1
    while true
        if i > lastindex(m) || m[i] != m0
            elements = i-i0
            x *= binomial(group_sizes[m0]+elements-1, elements)
            if i > lastindex(m)
                return x
            end
            i0 = i
            m0 = m[i]
        end
        i += 1
    end
end

"""
    rand(s::DCHSBM_sampler[, edges::Integer])

Draw a hypergraph from the sampler `s` with `edges` edges
"""
function Base.rand(rng::AbstractRNG, s::DCHSBM_sampler)
    map(group -> rand(rng, s.groups[group]), rand(rng, s.ms))
end

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
################################################################################
#
#  FakeFmpqMat.jl : Model fmpq_mat's as fmpz_mat's with denominator
#
################################################################################

export num, den

num(x::FakeFmpqMat) = x.num

den(x::FakeFmpqMat) = x.den

rows(x::FakeFmpqMat) = rows(x.num)

cols(x::FakeFmpqMat) = cols(x.num)

function simplify_content!(x::FakeFmpqMat)
  c = content(x.num)
  c = gcd(c, x.den)
  if c != 1
    divexact!(x.num, x.num, c)
    x.den = div(x.den, c)
  end
  y = canonical_unit(x.den)
  if !isone(y)
    x.den *= y
    mul!(x.num, x.num, y)
  end
end

################################################################################
#
#  Hashing
#
################################################################################

function hash(a::FakeFmpqMat, b::UInt)
  h = Base.hash(num(a)) $ Base.hash(den(a))
  h = h $ Base.hash(b)
  h = (h << 1) | (h >> (sizeof(Int)*8 - 1))
  return h
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, x::FakeFmpqMat)
  print(io, "FakeFmpqMat with numerator\n", num(x), "\nand denominator ", den(x))
end

################################################################################
#
#  Unary operations
#
################################################################################

function -(x::FakeFmpqMat)
  z = parent(x)()
  z.num = -x.num
  return z
end

function inv(x::FakeFmpqMat)
  i, d_i = pseudo_inv(num(x))
  i *= den(x)
  z = FakeFmpqMat(i,d_i)
  simplify_content!(z)
  return z
end

################################################################################
#
#  Binary operations
#
################################################################################

function +(x::FakeFmpqMat, y::FakeFmpqMat)
  t = den(y)*num(x) + den(x)*num(y)
  d = den(x)*den(y)
  z = FakeFmpqMat(t,d)
  simplify_content!(z)
  return z
end

function *(x::FakeFmpqMat, y::FakeFmpqMat)
  t = num(x)*num(y)
  d = den(x)*den(y)
  z = FakeFmpqMat(t,d)
  simplify_content!(z)
  return z
end

################################################################################
#
#  Adhoc binary operations
#
################################################################################

function *(x::FakeFmpqMat, y::fmpz_mat)
  t = x.num*y
  z = FakeFmpqMat(t, den(x))
  simplify_content!(z)
  return z
end

function *(x::fmpz_mat, y::FakeFmpqMat)
  t = x*y.num
  z = FakeFmpqMat(t, den(y))
  simplify_content!(z)
  return z
end

################################################################################
#
#  Comparison
#
################################################################################

==(x::FakeFmpqMat, y::FakeFmpqMat) = (num(x) == num(y)) && (den(x) == den(y))

################################################################################
#
#  Conversion
#
################################################################################

to_array(x::FakeFmpqMat) = (x.num, x.den)

function to_fmpz_mat(x::FakeFmpqMat)
  den(x) != 1 && error("Denominator has to be 1")
  return num(x)
end

################################################################################
#
#  Hermite normal form for numerator
#
################################################################################

function hnf!(x::FakeFmpqMat)
  h = hnf(num(x))
  x.num = h
end

function hnf(x::FakeFmpqMat, shape = :lowerleft)
  h = _hnf(num(x), shape)
  return FakeFmpqMat(h,den(x))
end

################################################################################
#
#  Sub
#
################################################################################

function sub(x::FakeFmpqMat, r::UnitRange{Int}, c::UnitRange{Int})
  z = FakeFmpqMat(sub(num(x),r,c),den(x))
  return z
end

function deepcopy(x::FakeFmpqMat)
  z = FakeFmpqMat(deepcopy(num(x)), den(x))
  return z
end

################################################################################
#
#  Determinant
#
################################################################################

function det(x::FakeFmpqMat)
  rows(x) != cols(x) && error("Matrix must be square")

  return det(x.num)//(x.den)^rows(x)
end


text:
using NPZ
using LinearAlgebra
using .Threads
using Arpack
using ProgressMeter
import Base: isless

struct TestResult
    i::Int64
    f::Float64
end
isless(a::TestResult, b::TestResult) = isless(a.f, b.f)

function main()
    X1 = npzread("Kernels/G1.npy")
    X2 = npzread("Kernels/G2.npy")
    X3 = npzread("Kernels/G3.npy")
    nfactor = npzread("Testing/nfactor.npy")

    n = size(X3, 1)
    ndas = size(X1, 1)

    T = Float64
    ϵ = 1.0
    λ = 100.0
    α = 100.0

    X = [X3 zeros(T, size(X3));
        zeros(T, size(X3)) X3;
        nfactor*X1 nfactor*X2]
    Xn = maximum(sqrt.(sum(X.^2, dims=2)))
    X ./= Xn

    n_eigs = 5

    v, V = eigs(X'*X + ϵ*I; nev=n_eigs, which=:LM)

    Γ = V[:, (end-n_eigs+1):end]

    function Ml(μ, X, ϵ)
        M = ϵ*I + X'*Diagonal(μ)*X
        Mc = cholesky(Hermitian(M))
        return (M, Mc)
    end

    Ml(μ) = Ml(μ, X, ϵ)

    function f(μ, Γ, α, λ, n_eigs)
        M, Mc = Ml(μ)
        vu, Vu = eigs(M; nev=5, which=:LM)
        U = Vu[:, (end-n_eigs+1):end]
        return -logdet(M) + λ*tr(U'*(Mc\U)) + α*sum((U.-Γ).^2)
    end

    f(μ) = f(μ, Γ, α, λ, n_eigs)

    function testf(u, i)
        u_test = copy(u)
        u_test[i] = 1.0
        u_test[i+256] = 1.0
        return TestResult(i, f(u_test))
    end

    u = vcat(zeros(256+256), ones(size(X1,1)))
    to_test = Set(1:256)
    selection_order = Int[]
    f_vals = T[]
    p = Progress(128*257, dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)

    while length(to_test) > 0
        test_vec = TestResult[]
        @threads for i in collect(to_test)
            push!(test_vec, testf(u, i))
            next!(p)
        end
        best = minimum(test_vec)
        delete!(to_test, best.i)
        push!(selection_order, best.i)
        push!(f_vals, best.f)
        u[best.i] = 1.0
        u[best.i+256] = 1.0
    end

    npzwrite("seq_order.npy", selection_order)
    npzwrite("seq_vals.npy", f_vals)

end

main()

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
#julia4 --machinefile=$PBS_NODEFILE -- run100.jl 1 2 110
#run on all available cores using:
#addprocs(readcsv(ENV["PBS_NODEFILE"], ASCIIString)[:, 1])

if nprocs() < 2
    if ARGS[4] == "slurm"
        using ClusterManagers
        np=parse(Int, ENV["SLURM_NTASKS"])
        addprocs(SlurmManager(np))
    elseif ARGS[4] == "pbs"
        addprocs(readcsv(ENV["PBS_NODEFILE"], ASCIIString)[:, 1])
    end
end

@everywhere args = @fetchfrom 1 ARGS

include(joinpath(Pkg.dir("LatentGaussianMixtureModel"), "teststat200/TestAsymptoticDistribution.jl"))
@everywhere Ctrue = parse(Int, args[1])
@everywhere Calternative = parse(Int, args[2])
teststat = pmap((b) -> Brun(b, Ctrue, Calternative, showpower=true), 1:parse(Int, args[3]))
writecsv("seq$(Ctrue)$(Calternative).csv", teststat)
run(`sed -i "s/(//g" seq$(Ctrue)$(Calternative).csv`)
run(`sed -i "s/)//g" seq$(Ctrue)$(Calternative).csv`)


text:
using SparseArrays
using LinearAlgebra
using Plots

function invvekit(A, y, m, e)
    l=[0];
    # Berechnung gemäß Algorithmus (Arbeitsblatt)
    Q,R,p = lu(A-m*Diagonal(ones(size(A))))
    while true
        y_h=R\(Q' * y) #TODO: geht nicht immer, warum
        l_h = 1/(dot(y_h,y)) + m
        # Närung Eigenwert zu Liste Eigenwert hinzufügen
        l=[l;l_h]
        # neuen Wert Eigenvektor setzen
        y=y_h/norm(y_h)
        # Abbruchkriterium prüfen
        if (norm(A*y-l_h*y)<e)
            # Liste Eigenwert und Eigenvektor zurückgeben
            return l,y;
        end
    end
end

max_stars_repo_name: SoulKathi/vekit

max_issues_repo_name: SoulKathi/vekit

max_forks_repo_name: SoulKathi/vekit


text:
include_dpo()

# Model
include_model("double_integrator")

# Horizon
T = 51

# Bounds

# ul <= u <= uu
ul, uu = control_bounds(model, T, 0.0, 0.0)

# Initial and final states
xl, xu = state_bounds(model, T, zeros(model.n), zeros(model.n))

# Problem
prob_nom = trajectory_optimization(
			model,
			EmptyObjective(),
			T,
			xl = xl,
			xu = xu,
			ul = ul,
			uu = uu)

# DPO
β = 1.0
δ = 1.0

x1 = resample(ones(model.n), Diagonal(ones(model.n)), β)

# Mean problem
prob_mean = trajectory_optimization(
				model,
				EmptyObjective(),
				T,
				dynamics = false)

# Sample problems
prob_sample = [trajectory_optimization(
				model,
				EmptyObjective(),
				T,
				dynamics = false,
				xl = state_bounds(model, T, x1 = x1[i])[1],
				xu = state_bounds(model, T, x1 = x1[i])[2]) for i = 1:2 * model.n]

# Sample objective
Q = [Diagonal(ones(model.n)) for t = 1:T]
R = [Diagonal(ones(model.m)) for t = 1:T-1]

obj_sample = sample_objective(Q, R)
policy = linear_feedback(model.n, model.m)
dist = disturbances([Diagonal(δ * ones(model.d)) for t = 1:T-1])
sample = sample_params(β, T)

prob_dpo = dpo_problem(
	prob_nom, prob_mean, prob_sample,
	obj_sample,
	policy,
	dist,
	sample)

z0 = zeros(prob_dpo.num_var)
include_snopt()
z, info = solve(prob_dpo, copy(z0),
	tol = 1.0e-7, c_tol = 1.0e-7,
	nlp = :SNOPT7)

# Solve
if true
	include_snopt()
	z, info = solve(prob_dpo, copy(z0),
		tol = 1.0e-7, c_tol = 1.0e-7,
		nlp = :SNOPT7)
	@save joinpath(@__DIR__, "sol_dpo.jld2") z
else
	println("Loading solution...")
	@load joinpath(@__DIR__, "sol_dpo.jld2") z
end

# DPO policy
θ = get_policy(z, prob_dpo)

# TVLQR policy
A, B = get_dynamics(model)
K, P = tvlqr(
	[A for t = 1:T-1],
	[B for t = 1:T-1],
	Q, R)

# Policy difference
policy_diff = [norm(vec(θ[t] - K[t])) / norm(vec(K[t])) for t = 1:T-1]
println("policy difference (inf. norm): $(norm(policy_diff, Inf))")

# Monte Carlo
using Random, Distributions
Random.seed!(1)

function monte_carlo(;M = 5)
	pd = []
	uni_dist = Distributions.Uniform(-1.0, 1.0)
	z_failure = []
	for i = 1:M
		z0 = rand(uni_dist, prob_dpo.num_var)
		z, info = solve(prob_dpo, copy(z0),
			tol = 1.0e-7, c_tol = 1.0e-7,
			nlp = :SNOPT7,
			time_limit = 60)
			# mapl = 5, mipl = 0)

		θ = get_policy(z, prob_dpo)
		policy_diff = [norm(vec(θ[t] - K[t])) / norm(vec(K[t])) for t = 1:T-1]
		pdi = norm(policy_diff, Inf)

		# if SNOPT fails, try Ipopt
		if pdi > 1.0e-4
			z, info = solve(prob_dpo, copy(z0),
				tol = 1.0e-9, c_tol = 1.0e-9,
				nlp = :ipopt,
				max_iter = 250)
			θ = get_policy(z, prob_dpo)
			policy_diff = [norm(vec(θ[t] - K[t])) / norm(vec(K[t])) for t = 1:T-1]
			pdi = norm(policy_diff, Inf)

			if pdi > 1.0e-4
				push!(z_failure, (z0, pdi))
			end
		end

		push!(pd, pdi)

	end

	return pd, z_failure
end

M = 1000
pd, z_failure = monte_carlo(M = M)
println("Monte Carlo: M = $M")
println("failures: $(count(pd .> 1.0e-4))")
println("maximum diff.: $(maximum(pd[pd .<= 1.0e-4]))")
println("mean diff.: $(mean(pd[pd .<= 1.0e-4]))")
println("std diff.: $(std(pd[pd .<= 1.0e-4]))")
#
# norm(z_failure[1][1])
# z, info = solve(prob_dpo, copy(z_failure[1][1]),
# 	tol = 1.0e-9, c_tol = 1.0e-9,
# 	nlp = :ipopt,
# 	time_limit = 20)
# θ = get_policy(z, prob_dpo)
# policy_diff = [norm(vec(θ[t] - K[t])) / norm(vec(K[t])) for t = 1:T-1]
# pdi = norm(policy_diff, Inf)


text:
## Core CAOL routines

using OffsetArrays, ImageFiltering, LinearAlgebra

# Allocation and initialization
function _initvars(x,H0,R)
    K = size(H0,2)

    # Form padded images
    xpad = [padarray(xl,Pad(:circular,ntuple(_->0,ndims(xl)),R)) for xl in x]
    # xpad = [view(xl,ntuple(i->vcat(1:size(xl,i),1:R[i]),ndims(xl))...) for xl in x]
    # more memory efficient but seems slower, make a copy for now
    # may be relevant: https://docs.julialang.org/en/v1/manual/performance-tips/#Copying-data-is-not-always-bad-1

    # Initialize filters
    H = copy(H0)
    h = _filterlist(H,R)
    Hprev = similar(H)  # for stopping condition

    # Initialize temporary variables
    zlk = similar(first(x),map(n->0:n-1,size(first(x))))
    ΨZ = similar(H)
    ψz = _filterlist(ΨZ,R)
    ψztemp = similar(first(ψz))
    HΨZ = similar(H,K,K)
    UVt = HΨZ

    return xpad, H, h, Hprev, zlk, ΨZ, ψz, ψztemp, HΨZ, UVt
end

# Update steps
@inline function _updateΨZ!(ΨZ,ψz,xpad,h,λ,zlk,ψztemp)
    obj = 0.0
    fill!(ΨZ,zero(eltype(ΨZ)))
    for xpadl in xpad, k in 1:length(h)
        imfilter!(zlk,xpadl,(h[k],),NoPad(),Algorithm.FIR())
        obj += _obj(zlk,λ)
        zlk .= hard.(zlk,sqrt(2λ))
        imfilter!(ψztemp,xpadl,(zlk,),NoPad(),Algorithm.FIR())
        ψz[k] .+= ψztemp
    end
    return obj
end
function _updateH!(H,ΨZ,H0,HΨZ,UVt)
    mul!(HΨZ,H0',ΨZ)
    F = svd!(HΨZ)
    mul!(UVt,F.U,F.Vt)
    mul!(H,H0,UVt)
end

# Iterative algorithms
function _CAOL(x,λ,H0,R,maxiters,tol)
    @assert H0'H0 ≈ (1/prod(R))*I
    xpad, H, h, Hprev, zlk, ΨZ, ψz, ψztemp, HΨZ, UVt = _initvars(x,H0,R)

    for t in 1:maxiters
        copyto!(Hprev,H)                      # Copy previous filters
        _updateΨZ!(ΨZ,ψz,xpad,h,λ,zlk,ψztemp) # Update ΨZ
        _updateH!(H,ΨZ,H0,HΨZ,UVt)            # Update filters as polar factor

        sqrt(sosdiff(Hprev,H) / (size(H,2)/size(H,1))) <= tol && break
    end

    return H
end
function _CAOLtrace(x,λ,H0,R,maxiters,tol)
    @assert H0'H0 ≈ (1/prod(R))*I
    xpad, H, h, Hprev, zlk, ΨZ, ψz, ψztemp, HΨZ, UVt = _initvars(x,H0,R)

    Htrace     = fill(H0,0)
    objtrace   = fill(NaN,0)
    Hdifftrace = fill(NaN,0)

    for t in 1:maxiters
        copyto!(Hprev,H)                            # Copy previous filters
        obj = _updateΨZ!(ΨZ,ψz,xpad,h,λ,zlk,ψztemp) # Compute objective, update ΨZ
        _updateH!(H,ΨZ,H0,HΨZ,UVt)                  # Update filters as polar factor

        push!(Htrace,copy(H))
        push!(objtrace,obj)
        push!(Hdifftrace, sqrt(sosdiff(Hprev,H) / (size(H,2)/size(H,1))))

        Hdifftrace[end] <= tol && break
    end

    return H, Htrace, objtrace, Hdifftrace
end


text:
module UtilitiesTests

import Compat
using Compat.Test
import Compat: occursin
import Compat.Base64: stringmime

import Documenter
import Documenter: IdDict

module UnitTests
    module SubModule end

    # Does `submodules` collect *all* the submodules?
    module A
        module B
            module C
                module D end
            end
        end
    end

    mutable struct T end
    mutable struct S{T} end

    "Documenter unit tests."
    Base.length(::T) = 1

    f(x) = x

    const pi = 3.0
end

module OuterModule
module InnerModule
import ..OuterModule
export OuterModule
end
end

@testset "Utilities" begin
    let doc = @doc(length)
        a = Documenter.Utilities.filterdocs(doc, Set{Module}())
        b = Documenter.Utilities.filterdocs(doc, Set{Module}([UnitTests]))
        c = Documenter.Utilities.filterdocs(doc, Set{Module}([Base]))
        d = Documenter.Utilities.filterdocs(doc, Set{Module}([UtilitiesTests]))

        @test a !== nothing
        @test a === doc
        @test b !== nothing
        @test occursin("Documenter unit tests.", stringmime("text/plain", b))
        @test c !== nothing
        @test !occursin("Documenter unit tests.", stringmime("text/plain", c))
        @test d === nothing
    end

    # Documenter.Utilities.issubmodule
    @test Documenter.Utilities.issubmodule(Main, Main) === true
    @test Documenter.Utilities.issubmodule(UnitTests, UnitTests) === true
    @test Documenter.Utilities.issubmodule(UnitTests.SubModule, Main) === true
    @test Documenter.Utilities.issubmodule(UnitTests.SubModule, UnitTests) === true
    @test Documenter.Utilities.issubmodule(UnitTests.SubModule, Base) === false
    @test Documenter.Utilities.issubmodule(UnitTests, UnitTests.SubModule) === false

    @test UnitTests.A in Documenter.Utilities.submodules(UnitTests.A)
    @test UnitTests.A.B in Documenter.Utilities.submodules(UnitTests.A)
    @test UnitTests.A.B.C in Documenter.Utilities.submodules(UnitTests.A)
    @test UnitTests.A.B.C.D in Documenter.Utilities.submodules(UnitTests.A)
    @test OuterModule in Documenter.Utilities.submodules(OuterModule)
    @test OuterModule.InnerModule in Documenter.Utilities.submodules(OuterModule)
    @test length(Documenter.Utilities.submodules(OuterModule)) == 2

    @test Documenter.Utilities.isabsurl("file.md") === false
    @test Documenter.Utilities.isabsurl("../file.md") === false
    @test Documenter.Utilities.isabsurl(".") === false
    @test Documenter.Utilities.isabsurl("https://example.org/file.md") === true
    @test Documenter.Utilities.isabsurl("http://example.org") === true
    @test Documenter.Utilities.isabsurl("ftp://user:pw@example.org") === true
    @test Documenter.Utilities.isabsurl("/fs/absolute/path") === false

    @test Documenter.Utilities.doccat(UnitTests) == "Module"
    @test Documenter.Utilities.doccat(UnitTests.T) == "Type"
    @test Documenter.Utilities.doccat(UnitTests.S) == "Type"
    @test Documenter.Utilities.doccat(UnitTests.f) == "Function"
    @test Documenter.Utilities.doccat(UnitTests.pi) == "Constant"

    # repo type
    @test Documenter.Utilities.repo_host_from_url("https://bitbucket.org/somerepo") == Documenter.Utilities.RepoBitbucket
    @test Documenter.Utilities.repo_host_from_url("https://www.bitbucket.org/somerepo") == Documenter.Utilities.RepoBitbucket
    @test Documenter.Utilities.repo_host_from_url("http://bitbucket.org/somethingelse") == Documenter.Utilities.RepoBitbucket
    @test Documenter.Utilities.repo_host_from_url("http://github.com/Whatever") == Documenter.Utilities.RepoGithub
    @test Documenter.Utilities.repo_host_from_url("https://github.com/Whatever") == Documenter.Utilities.RepoGithub
    @test Documenter.Utilities.repo_host_from_url("https://www.github.com/Whatever") == Documenter.Utilities.RepoGithub
    @test Documenter.Utilities.repo_host_from_url("https://gitlab.com/Whatever") == Documenter.Utilities.RepoGitlab

    # line range
    let formatting = Documenter.Utilities.LineRangeFormatting(Documenter.Utilities.RepoGithub)
        @test Documenter.Utilities.format_line(1:1, formatting) == "L1"
        @test Documenter.Utilities.format_line(123:123, formatting) == "L123"
        @test Documenter.Utilities.format_line(2:5, formatting) == "L2-L5"
        @test Documenter.Utilities.format_line(100:9999, formatting) == "L100-L9999"
    end

    let formatting = Documenter.Utilities.LineRangeFormatting(Documenter.Utilities.RepoGitlab)
        @test Documenter.Utilities.format_line(1:1, formatting) == "L1"
        @test Documenter.Utilities.format_line(123:123, formatting) == "L123"
        @test Documenter.Utilities.format_line(2:5, formatting) == "L2-5"
        @test Documenter.Utilities.format_line(100:9999, formatting) == "L100-9999"
    end

    let formatting = Documenter.Utilities.LineRangeFormatting(Documenter.Utilities.RepoBitbucket)
        @test Documenter.Utilities.format_line(1:1, formatting) == "1"
        @test Documenter.Utilities.format_line(123:123, formatting) == "123"
        @test Documenter.Utilities.format_line(2:5, formatting) == "2:5"
        @test Documenter.Utilities.format_line(100:9999, formatting) == "100:9999"
    end

    # URL building
    filepath = string(first(methods(Documenter.Utilities.url)).file)
    Compat.Sys.iswindows() && (filepath = replace(filepath, "/" => "\\")) # work around JuliaLang/julia#26424
    let expected_filepath = "/src/Utilities/Utilities.jl"
        Compat.Sys.iswindows() && (expected_filepath = replace(expected_filepath, "/" => "\\"))
        @test endswith(filepath, expected_filepath)
        @show filepath expected_filepath
    end
    commit = Documenter.Utilities.repo_commit(filepath)

    @test Documenter.Utilities.url("//blob/{commit}{path}#{line}", filepath) == "//blob/$(commit)/src/Utilities/Utilities.jl#"
    @test Documenter.Utilities.url(nothing, "//blob/{commit}{path}#{line}", Documenter.Utilities, filepath, 10:20) == "//blob/$(commit)/src/Utilities/Utilities.jl#L10-L20"

    # repo_root & relpath_from_repo_root
    @test Documenter.Utilities.repo_root(@__FILE__) == dirname(abspath(joinpath(@__DIR__, ".."))) # abspath() keeps trailing /, hence dirname()
    @test Documenter.Utilities.repo_root(@__FILE__; dbdir=".svn") == nothing
    @test Documenter.Utilities.relpath_from_repo_root(@__FILE__) == joinpath("test", "utilities.jl")
    # We assume that a temporary file is not in a repo
    @test Documenter.Utilities.repo_root(tempname()) == nothing
    @test Documenter.Utilities.relpath_from_repo_root(tempname()) == nothing

    import Documenter.Documents: Document, Page, Globals
    let page = Page("source", "build", [], IdDict(), Globals()), doc = Document()
        code = """
        x += 3
        γγγ_γγγ
        γγγ
        """
        exprs = Documenter.Utilities.parseblock(code, doc, page)

        @test isa(exprs, Vector)
        @test length(exprs) === 3

        @test isa(exprs[1][1], Expr)
        @test exprs[1][1].head === :+=
        @test exprs[1][2] == "x += 3\n"

        @test exprs[2][2] == "γγγ_γγγ\n"

        @test exprs[3][1] === :γγγ
        @test exprs[3][2] == "γγγ\n"
    end

    @testset "TextDiff" begin
        import Documenter.Utilities.TextDiff: splitby
        @test splitby(r"\s+", "X Y  Z") == ["X ", "Y  ", "Z"]
        @test splitby(r"[~]", "X~Y~Z") == ["X~", "Y~", "Z"]
        @test splitby(r"[▶]", "X▶Y▶Z") == ["X▶", "Y▶", "Z"]
        @test splitby(r"[▶]+", "X▶▶Y▶Z▶") == ["X▶▶", "Y▶", "Z▶"]
        @test splitby(r"[▶]+", "▶▶Y▶Z▶") == ["▶▶", "Y▶", "Z▶"]
        @test splitby(r"[▶]+", "Ω▶▶Y▶Z▶") == ["Ω▶▶", "Y▶", "Z▶"]
        @test splitby(r"[▶]+", "Ω▶▶Y▶Z▶κ") == ["Ω▶▶", "Y▶", "Z▶", "κ"]
    end

    @static if isdefined(Base, :with_logger)
        @testset "withoutput" begin
            _, _, _, output = Documenter.Utilities.withoutput() do
                println("println")
                @info "@info"
                f() = (Base.depwarn("depwarn", :f); nothing)
                f()
            end
            @test startswith(output, "println\n[ Info: @info\n┌ Warning: depwarn\n")
        end
    end
end

end


text:
#
# rUK taxes 2022/3, where these are different.
#
sys.it.non_savings_rates = [20.0,40.0,45.0]
sys.it.non_savings_thresholds = [37_700, 150_000.0]
sys.it.non_savings_basic_rate = 1 # above this counts as higher rate

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
@testset "Access and Iterations" begin
    dna_kmer = mer"ACTG"dna
    rna_kmer = mer"ACUG"rna
    big_dna_kmer = bigmer"ACTG"dna
    big_rna_kmer = bigmer"ACUG"rna

    @testset "Access DNA Kmer" begin
        @test dna_kmer[1] == DNA_A
        @test dna_kmer[2] == DNA_C
        @test dna_kmer[3] == DNA_T
        @test dna_kmer[4] == DNA_G

        @test big_dna_kmer[1] == DNA_A
        @test big_dna_kmer[2] == DNA_C
        @test big_dna_kmer[3] == DNA_T
        @test big_dna_kmer[4] == DNA_G

        # Access indexes out of bounds
        @test_throws Exception dna_kmer[-1]
        @test_throws Exception dna_kmer[0]
        @test_throws Exception dna_kmer[5]
        @test_throws Exception getindex(dna_kmer,-1)
        @test_throws Exception getindex(dna_kmer, 0)
        @test_throws Exception getindex(dna_kmer, 5)

        @test_throws Exception big_dna_kmer[-1]
        @test_throws Exception big_dna_kmer[0]
        @test_throws Exception big_dna_kmer[5]
        @test_throws Exception getindex(big_dna_kmer,-1)
        @test_throws Exception getindex(big_dna_kmer, 0)
        @test_throws Exception getindex(big_dna_kmer, 5)
    end

    @testset "Iteration through DNA Kmer" begin
        @test iterate(DNAMer("ACTG")) == (DNA_A, 2)
        @test iterate(BigDNAMer("ACTG")) == (DNA_A, 2)

        @test iterate(DNAMer("ACTG"), 1) == (DNA_A, 2)
        @test iterate(DNAMer("ACTG"), 4) == (DNA_G, 5)
        @test iterate(BigDNAMer("ACTG"), 1) == (DNA_A, 2)
        @test iterate(BigDNAMer("ACTG"), 4) == (DNA_G, 5)

        @test iterate(DNAMer("ACTG"), 1)  !== nothing
        @test iterate(DNAMer("ACTG"), 4)  !== nothing
        @test iterate(DNAMer("ACTG"), 5)  === nothing
        @test iterate(DNAMer("ACTG"), -1) !== nothing

        @test iterate(BigDNAMer("ACTG"), 1)  !== nothing
        @test iterate(BigDNAMer("ACTG"), 4)  !== nothing
        @test iterate(BigDNAMer("ACTG"), 5)  === nothing
        @test iterate(BigDNAMer("ACTG"), -1) !== nothing

        dna_vec = [DNA_A, DNA_C, DNA_T, DNA_G]
        @test all([nt === dna_vec[i] for (i, nt) in enumerate(dna_kmer)])
        @test all([nt === dna_vec[i] for (i, nt) in enumerate(big_dna_kmer)])
    end

    @testset "Access RNA Kmer" begin
        @test rna_kmer[1] == RNA_A
        @test rna_kmer[2] == RNA_C
        @test rna_kmer[3] == RNA_U
        @test rna_kmer[4] == RNA_G

        @test big_rna_kmer[1] == RNA_A
        @test big_rna_kmer[2] == RNA_C
        @test big_rna_kmer[3] == RNA_U
        @test big_rna_kmer[4] == RNA_G

        # Access indexes out of bounds
        @test_throws Exception rna_kmer[-1]
        @test_throws Exception rna_kmer[0]
        @test_throws Exception rna_kmer[5]
        @test_throws Exception getindex(rna_kmer, -1)
        @test_throws Exception getindex(rna_kmer, 0)
        @test_throws Exception getindex(rna_kmer, 5)

        @test_throws Exception big_rna_kmer[-1]
        @test_throws Exception big_rna_kmer[0]
        @test_throws Exception big_rna_kmer[5]
        @test_throws Exception getindex(big_rna_kmer, -1)
        @test_throws Exception getindex(big_rna_kmer, 0)
        @test_throws Exception getindex(big_rna_kmer, 5)
    end

    @testset "Iteration through RNA Kmer" begin
        @test iterate(RNAMer("ACUG")) == (RNA_A, 2)
        @test iterate(BigRNAMer("ACUG")) == (RNA_A, 2)

        @test iterate(RNAMer("ACUG"), 1) == (RNA_A, 2)
        @test iterate(RNAMer("ACUG"), 4) == (RNA_G, 5)
        @test iterate(BigRNAMer("ACUG"), 1) == (RNA_A, 2)
        @test iterate(BigRNAMer("ACUG"), 4) == (RNA_G, 5)

        @test iterate(RNAMer("ACUG"), 1)  !== nothing
        @test iterate(RNAMer("ACUG"), 4)  !== nothing
        @test iterate(RNAMer("ACUG"), 5)  === nothing
        @test iterate(RNAMer("ACUG"), -1) !== nothing

        @test iterate(BigRNAMer("ACUG"), 1)  !== nothing
        @test iterate(BigRNAMer("ACUG"), 4)  !== nothing
        @test iterate(BigRNAMer("ACUG"), 5)  === nothing
        @test iterate(BigRNAMer("ACUG"), -1) !== nothing

        rna_vec = [RNA_A, RNA_C, RNA_U, RNA_G]
        @test all([nt === rna_vec[i] for (i, nt) in enumerate(rna_kmer)])
        @test all([nt === rna_vec[i] for (i, nt) in enumerate(big_rna_kmer)])
    end

    @testset "Take slices of Mers" begin
	@test DNAMer("ACGT")[2:3] == DNAMer("CG")
	@test BigDNAMer("ACGT")[2:3] == BigDNAMer("CG")
	@test RNAMer("ACGU")[2:3] == RNAMer("CG")
	@test BigRNAMer("ACGU")[2:3] == BigRNAMer("CG")

	@test DNAMer("ACGT")[2:1:3] == DNAMer("CG")
	@test DNAMer("ACGT")[1:2:end] == DNAMer("AG")
	@test DNAMer("ACGT")[3:-1:2] == DNAMer("GC")

	@test_throws ArgumentError DNAMer("ACGT")[3:2]
	@test_throws BoundsError DNAMer("ACGT")[3:5]
    end
end