text:
"""
Dynamic Programming Tree Solver
===============================

Merge the event queues within the tree by sorting arrays.

Compared to previous versions, all weighting (λ and μ) is templated;
thus we only need one implementation and the compiler will do the
rest.

"""
module TreeDP

import GraphIdx
import GraphIdx: Ones, Const, Vec, Weights, is_const
import GraphIdx.Utils: MutRef
import GraphIdx.Tree: ChildrenIndex, reset!, dfs_walk
import ..Pwl: clip, Range, Event, EPS, DEBUG
import ..Pwl: clip_front, clip_back
import ..Pwl.QueueUnion: Queues, merge, reset!

const Tree = GraphIdx.Tree.RootedTree

"""
Contains all memory needed for `tree_dp!`.
"""
struct TreeDPMem{F,I}
    lb::Vector{F}
    queues::Queues
    proc_order::Vector{I}
    kidz::ChildrenIndex
    stack::Vector{Int}
end

function TreeDPMem(n::Integer, F::Type = Float64, I::Type = Int)
    lb = Vector{F}(undef, n)
    proc_order = Vector{I}(undef, n)
    kidz = ChildrenIndex(n)
    stack = Vector{Int}(undef, n)
    TreeDPMem{F,I}(lb, Queues{Event}(n), proc_order, kidz, stack)
end

"""
    reset!(mem::TreeDPMem, t::Tree)

Re-initialize the memory given the new tree, i.e.
0. Assert that `mem` has the required size for `t`.
1. Compute processing order.
2. Initialize queues to fit `t`.
"""
function reset!(mem::TreeDPMem{F,I}, tree::Tree) where {F,I}
    reset!(mem.kidz, tree.parent, tree.root)
    sizehint!(mem.proc_order, length(tree))
    reset!(mem.queues.pq, mem.kidz, mem.stack, mem.proc_order)
    mem
end

"""
    tree_dp(y, t::Tree, λ [, μ = 1.0])

Perform the dynamic programming algorithm on `y` and return the optimized `x`.

The edge weighting λ should be either a constant (e.g. `Float64`) or a
callable such that `λ(i)` returns the weight of the edge `(i, t.parent[i])`;
for `λ(t.root)` it might return anything but not causean error.
"""
tree_dp(y::Array{F,N}, t::Tree, λ::F)  where {F,N} =
    tree_dp!(similar(y), y, t, Const{F}(λ), Ones{F}())

tree_dp(y::Array{F,N}, t::Tree, λ::F, µ::F)  where {F,N} =
    tree_dp!(similar(y), y, t, Const{F}(λ), Const{F}(µ))

tree_dp(y::Array{F,N}, t::Tree, λ::Lam, µ::Mu) where {F,N,Lam,Mu} =
    tree_dp!(similar(y), y, t, λ, µ)

"""
    tree_dp!(x, y, t, λ, μ [,mem])

Like `tree_dp` but do not allocate an output buffer for `x`.

If `mem::TreeDPMem` is provided, no additional allocations will be
needed.
"""
tree_dp!(x::Array{F,N}, y::Array{F,N}, t::Tree, λ::La, µ::Mu) where {F,N,La,Mu} =
    tree_dp!(x, y, t, λ, µ, TreeDPMem(length(y)))

function tree_dp!(
    x::Array{F,N},
    y::Array{F,N},
    t::Tree,
    λ::Lam,
    µ::Mu,
    mem::TreeDPMem{F,I},
)::Array{F,N} where {F, I, N, Lam<:Weights, Mu<:Weights}
    tree_dp!(x, y, t, λ, μ, mem, Val(!is_const(Mu)))
end

function tree_dp!(
    x::Array{F,N},
    y::Array{F,N},
    t::Tree,
    λ::Lam,
    µ::Mu,
    mem::TreeDPMem{F,I},
    C::Val{check},
)::Array{F,N} where {F, I, N, Lam<:Weights, Mu<:Weights, check}
    reset!(mem, t)
    local lb::Vector{F} = mem.lb
    local ub::Vector{F} = vec(x)
    local sig::Vector{F} = mem.lb
    local ev::Vector{Event} = mem.queues.events
    local pq::Vector{Range} = mem.queues.pq
    local ref::MutRef{Range} = MutRef(pq, 1)

    sig .= 0
    for i in @view mem.proc_order[1:end-1]
        local sig_i::F = sig[i]
        ref.i = i
        lb[i] = clip(ev, ref, +µ[i], -µ[i]*y[i] - sig_i + λ[i], Val(true), C)
        ub[i] = clip(ev, ref, -µ[i], +µ[i]*y[i] - sig_i + λ[i], Val(false), C)
        merge(mem.queues.events::Vector{Event},
              mem.queues.pq::Vector{Range},
              t.parent[i]::Int, i::Int)
        sig[t.parent[i]] += μ[i] > EPS ? λ[i] : min(sig_i, λ[i])
    end

    let r = t.root
        ref.i = r
        x[r] = clip(ev, ref,  +µ[r], -µ[r]*y[r] -sig[r] + 0.0, Val(true), C)
    end
    for v in @view mem.proc_order[end-1:-1:1]
        x[v] = clamp(x[t.parent[v]], lb[v], ub[v])
    end

    return x
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
### Now

## TODO take 1 arg as well
@mkapprule Now :nargs => 0

@doap Now() = Dates.now()

### Timing

@sjdoc Timing """
    Timing(expr)

evaluate `expr` and return a `List` of the elapsed CPU time
and the result.
"""
@mkapprule Timing  :nodefault => true

@sjseealso_group(Timing,Allocated,Time,Trace)

@doap function Timing(exprs...)
    local mxnew
    t = @elapsed begin
        reset_meval_count()
        for x in exprs
            mxnew = doeval(x)
        end
        setsymval(:ans,mxnew)
    end
    mxpr(:List,t,mxnew)
end

### Pause

@sjdoc Pause """
    Pause(x)

pauses (i.e.sleeps) for `x` seconds.
"""
@mkapprule Pause  :nargs => 1
@doap Pause(x::T) where {T<:Real} = sleep(x)

## Partial implementations of these...

### UnixTime

@mkapprule UnixTime :nargs => 0

@doap UnixTime() = round(Int,time())

### AbsoluteTime

#mma2datetime(t) =

@mkapprule AbsoluteTime :nargs => 0

## FIXME: works in v0.6. Will break for earlier versions
@doap AbsoluteTime() = Dates.value(Dates.now() - Dates.DateTime(1900))/1000

### DateList

## Note: we can use DateParser.jl at some point.

function _datelist(t)
    MList(Dates.year(t),Dates.month(t),Dates.day(t),Dates.hour(t),Dates.minute(t),Dates.second(t))
end

@mkapprule DateList

@doap DateList(x::ListT) = _datelist(Dates.DateTime(margs(x)...))

@doap DateList() = _datelist(Dates.now())

#@doap DateList(t::AbstractFloat) = _datelist(Dates.DateTime(t))

### Run

@sjdoc Run """
    Run(\`cmd\`)

runs the symstem command `cmd`. Note the backticks around `cmd`.
"""
@mkapprule Run :nargs => 1

# eg Run(`echo hello`)
@doap Run(x::Cmd) = run(x)

### Environment

@mkapprule Environment :nargs => 1

## Mma returns $Failed on failure
@doap Environment(x::String) = get(ENV,x,Null)

### GetEnvironment

@mkapprule GetEnvironment

## Mma returns a rule. We just return the val here.
@doap GetEnvironment(x::String) = get(ENV,x,Null)

## We return a Dict, not a list of Rules
@doap GetEnvironment() = ENV

## Mma gets and receives objects somehow

### CopyToClipboard

@mkapprule CopyToClipboard :nargs => 1
@doap CopyToClipboard(x) = clipboard(wrapout(x))

### Paste

@mkapprule Paste :nargs => 0
@doap function Paste()
    try
        clipboard()
    catch
        Null
    end
end

ext: jl


text:
@testset "UniquePoints" begin
    data = [randn(ComplexF64, 12) for i = 1:10_000]
    up = HC.UniquePoints(first(data), 1)
    HC.add!.(up, data, 1:length(data), 1e-8)
    @test length(up) == 10_000
    HC.add!.(up, data, 1:length(data), 1e-8)

    # Test with group action
    x = randn(ComplexF64, 4)
    permutation1(x) = ([x[2]; x[1]; x[3]; x[4]],)
    permutation2(x) = ([x[1]; x[2]; x[4]; x[3]],)
    group_actions = GroupActions(permutation1, permutation1)
    X = [[x]; group_actions(x)]

    X = [v for v in GroupActions(permutation1, permutation2)(x)]

    # One group action
    data = HC.UniquePoints(first(X), 1, group_action = permutation1)
    HC.add!.(data, X, 1:length(X), 1e-5)
    @test length(data) == 2

    data =
        HC.UniquePoints(X[1], 1, group_actions = GroupActions(permutation1, permutation2))
    HC.add!.(data, X, 1:4, 1e-5)
    @test length(data) == 1

    # Group action and reality check
    x = randn(4)
    X = [im .* x, x, randn(ComplexF64, 4)]
    data = HC.UniquePoints(X[1], 1, group_action = x -> (im .* x, (-1) .* x, (-im) .* x))
    HC.add!.(data, X, 1:length(X), 1e-5)
    @test length(data) == 2
end

@testset "Multiplicities" begin
    # test empty input
    @test isempty(multiplicities(Vector{Vector{ComplexF64}}()))

    V = [randn(3) + im .* randn(3) for i = 1:10]
    W = [
        map(v -> v + [1e-7 * v[1]; 1.0; 1.0], V)
        V
        map(v -> v + [1e-6 * v[1]; 0; 0], V)
    ]
    U = [[cis(rand()) for _ = 1:3] for _ = 1:20]
    X = [
        [
            0.499064 - 0.707804im,
            -0.164223 + 0.0933025im,
            -0.00695332 - 0.139516im,
            0.000716617 - 0.135055im,
            0.009846 - 0.12984im,
            0.020204 - 0.123609im,
            0.34326 + 0.161402im,
        ],
        [
            -0.795428 - 0.0975931im,
            -0.20942 - 0.382162im,
            -0.119855 - 0.0414809im,
            -0.123105 - 0.0293171im,
            -0.128161 - 0.0148951im,
            -0.134641 + 0.00272884im,
            0.258239 + 0.185188im,
        ],
        [
            0.605937 - 0.289401im,
            0.42308 - 0.188036im,
            -0.057496 + 0.200812im,
            0.284444 + 0.0477843im,
            0.191189 - 0.225824im,
            -0.233732 - 0.181827im,
            -0.0939606 + 0.154403im,
        ],
        [
            -0.493372 + 0.512561im,
            -0.385166 + 0.204525im,
            -0.120853 - 0.142237im,
            0.0742939 + 0.163197im,
            -0.298456 + 0.202635im,
            0.100108 - 0.219369im,
            0.193075 + 0.105514im,
        ],
        [
            0.483174 - 0.581023im,
            0.23432 - 0.0528133im,
            -0.111264 - 0.067633im,
            -0.292867 - 0.014453im,
            0.24805 - 0.267034im,
            0.141232 - 0.13631im,
            0.168603 + 0.261794im,
        ],
        [
            -0.455116 + 0.4482im,
            -0.366004 + 0.264082im,
            -0.0578486 + 0.110708im,
            -0.296085 - 0.0619238im,
            0.216907 - 0.327015im,
            0.207934 + 0.242906im,
            -0.121251 + 0.101268im,
        ],
        [
            0.526416 - 0.465825im,
            0.39214 - 0.0908531im,
            -0.260222 + 0.00391701im,
            -0.110348 - 0.1151im,
            0.101057 - 0.190265im,
            0.308379 - 0.137148im,
            0.0488701 + 0.29665im,
        ],
        [
            0.147886 - 0.750864im,
            -0.235168 - 0.257487im,
            0.0394847 + 0.194512im,
            -0.0325228 - 0.227443im,
            -0.111526 - 0.214717im,
            0.0245065 + 0.162303im,
            0.304189 - 0.150105im,
        ],
        [
            -0.372175 - 0.538886im,
            -0.117563 - 0.341994im,
            0.0752294 - 0.0141672im,
            0.143671 - 0.249302im,
            -0.374165 - 0.228281im,
            -0.0971581 + 0.346633im,
            0.168499 + 0.0394346im,
        ],
        [
            -0.503759 + 0.48im,
            -0.447123 + 0.262057im,
            -0.214529 + 0.109917im,
            0.148493 - 0.0330798im,
            -0.0510215 - 0.262778im,
            -0.110165 + 0.198317im,
            0.200704 + 0.0508037im,
        ],
        [
            0.623685 - 0.513102im,
            0.0522679 + 0.0223332im,
            -0.188478 - 0.10101im,
            0.229926 - 0.157889im,
            -0.146424 - 0.0562049im,
            0.209112 - 0.185147im,
            0.205229 + 0.276164im,
        ],
        [
            0.127017 - 0.655061im,
            -0.147053 - 0.569588im,
            -0.134863 + 0.0901295im,
            0.152888 + 0.049255im,
            -0.124224 + 0.0624774im,
            0.152924 + 0.0778815im,
            0.105575 - 0.310906im,
        ],
        [
            0.575711 + 0.384346im,
            0.336129 + 0.0531313im,
            0.155698 - 0.152359im,
            -0.104555 + 0.197167im,
            -0.254208 - 0.215823im,
            0.327345 + 0.0728999im,
            -0.0274605 + 0.288691im,
        ],
        [
            0.479139 - 0.593553im,
            -0.0960412 - 0.315975im,
            0.0438803 - 0.245038im,
            0.000516985 + 0.193173im,
            -0.0828034 + 0.172297im,
            -0.00479366 - 0.21881im,
            0.353615 - 0.0168729im,
        ],
        [
            -0.37103 + 0.609109im,
            -0.201984 + 0.259807im,
            -0.158043 + 0.252162im,
            -0.211405 - 0.106454im,
            0.0274292 - 0.262351im,
            0.333102 + 0.0455649im,
            -0.0847849 + 0.220536im,
        ],
        [
            -0.200697 + 0.682536im,
            -0.247073 + 0.359358im,
            -0.172902 - 0.0702972im,
            0.141297 + 0.110399im,
            -0.171049 + 0.317616im,
            -0.0125869 - 0.240803im,
            0.219982 + 0.00435721im,
        ],
        [
            -0.0896611 - 0.805671im,
            -0.0393637 - 0.183565im,
            0.224508 - 0.167715im,
            0.0918438 - 0.180778im,
            -0.0879493 - 0.157071im,
            -0.250607 - 0.0700599im,
            0.0521971 + 0.291749im,
        ],
        [
            -0.226054 + 0.605127im,
            -0.210133 + 0.55141im,
            0.218162 - 0.0576634im,
            0.08772 - 0.0936im,
            -0.0820928 - 0.101495im,
            -0.22419 - 0.0507194im,
            -0.0146449 + 0.311527im,
        ],
        [
            0.498516 + 0.360158im,
            0.219394 + 0.0657585im,
            -0.0836376 + 0.0106629im,
            0.00702255 - 0.322935im,
            0.481634 + 0.13367im,
            -0.133322 + 0.397338im,
            -0.171875 - 0.0532251im,
        ],
        [
            0.440673 - 0.62568im,
            -0.110019 - 0.330907im,
            -0.043096 + 0.193744im,
            0.0628576 - 0.220991im,
            -0.0144861 - 0.24152im,
            -0.0436694 + 0.158228im,
            0.338948 - 0.0133119im,
        ],
        [
            0.622749 + 0.219955im,
            0.628838 + 0.216896im,
            0.0331851 + 0.145177im,
            -0.123675 - 0.0636893im,
            0.0348129 + 0.117408im,
            -0.133898 - 0.0420959im,
            0.155277 - 0.144857im,
        ],
        [
            0.652186 - 0.0569317im,
            0.552174 - 0.169926im,
            -0.0977967 - 0.0521929im,
            -0.161552 - 0.169677im,
            0.00437006 + 0.237406im,
            -0.0164646 + 0.116656im,
            0.310812 - 0.0599935im,
        ],
        [
            0.448293 + 0.574542im,
            0.247779 + 0.106134im,
            -0.0955677 - 0.215494im,
            0.0935008 + 0.252581im,
            0.336093 - 0.0655203im,
            -0.242972 - 0.0256941im,
            -0.0177608 + 0.301502im,
        ],
        [
            0.445422 + 0.671948im,
            -0.0726023 + 0.313409im,
            -0.044775 + 0.21726im,
            0.154346 + 0.0261763im,
            -0.0305679 + 0.188435im,
            0.20006 + 0.0391233im,
            -0.0810809 - 0.297159im,
        ],
        [
            -0.00979382 - 0.81029im,
            0.106391 - 0.195197im,
            0.108622 - 0.139201im,
            0.234362 - 0.1806im,
            -0.20136 - 0.131455im,
            -0.11823 - 0.118552im,
            -0.0308846 + 0.297292im,
        ],
        [
            0.318611 + 0.584183im,
            0.105547 + 0.221636im,
            0.248174 + 0.0583734im,
            -0.00256072 - 0.290164im,
            -0.302039 - 0.0824955im,
            -0.114841 + 0.409322im,
            0.25213 + 0.0736154im,
        ],
        [
            0.696901 + 0.0684572im,
            0.505426 + 0.266373im,
            -0.152493 + 0.0711982im,
            0.197532 - 0.0412531im,
            0.214441 + 0.0165324im,
            -0.096866 + 0.0847148im,
            -0.103987 - 0.201388im,
        ],
        [
            0.115188 - 0.657244im,
            -0.157297 - 0.566844im,
            -0.133216 + 0.092546im,
            0.153751 + 0.0464909im,
            -0.123077 + 0.0647066im,
            0.154303 + 0.0751121im,
            0.0999536 - 0.312759im,
        ],
        [
            0.653199 - 0.277909im,
            0.556407 - 0.00960311im,
            0.204776 - 0.112518im,
            -0.112199 + 0.0771487im,
            -0.0695128 + 0.163961im,
            0.167741 - 0.0723464im,
            -0.16148 - 0.148729im,
        ],
        [
            0.411538 + 0.550657im,
            0.413916 + 0.420472im,
            -0.0555052 + 0.089203im,
            -0.0609812 + 0.227821im,
            -0.0913546 - 0.187807im,
            -0.0607857 - 0.0922333im,
            0.237526 + 0.0198403im,
        ],
        [
            0.602049 + 0.184366im,
            0.695878 + 0.270759im,
            -0.0420156 + 0.0359236im,
            -0.0431941 + 0.0331572im,
            -0.0446403 + 0.0302374im,
            -0.046326 + 0.027106im,
            0.0979775 - 0.156785im,
        ],
        [
            -0.270799 + 0.567571im,
            -0.113714 + 0.377218im,
            0.0592799 - 0.219174im,
            0.226587 + 0.091218im,
            0.0269123 + 0.326853im,
            -0.402766 + 0.100383im,
            -0.018186 - 0.240642im,
        ],
    ]

    M = multiplicities(V)
    @test length(M) == 0

    N = multiplicities(W, metric = InfNorm(), atol = 1e-5)
    sort!(N, by = first)
    @test length(N) == 10
    @test unique([length(m) for m in N]) == [2]
    @test N == [[i, i + 10] for i = 11:20]

    O = multiplicities([U; U])
    @test length(O) == 20

    P = multiplicities(X, metric = (x, y) -> 1 - abs(LinearAlgebra.dot(x, y)), atol = 1e-5)
    @test length(P) == 3

    # Test with group action
    x = randn(ComplexF64, 4)
    permutation1(x) = ([x[2]; x[1]; x[3]; x[4]],)
    permutation2(x) = ([x[1]; x[2]; x[4]; x[3]],)
    X = GroupActions(permutation1, permutation2)(x)

    # One group action
    m = sort(multiplicities(X, group_action = permutation1))
    @test m == [[1, 2], [3, 4]]

    # Two group actions
    m = multiplicities(X, group_actions = GroupActions(permutation1, permutation2))
    @test m == [[1, 2, 3, 4]]

    # docs
    X = [[1, 2, 3, 4], [2, 1, 3, 4], [1, 2, 4, 3], [2, 1, 4, 3]]
    permutation(x) = [x[2], x[1], x[3], x[4]]
    m = sort(multiplicities(X, group_action = permutation))
    @test m == [[1, 2], [3, 4]]
end


text:
using DrWatson
using OptimalEstimationProject
using LinearAlgebra
using AstroTime
using DataFrames
using CSV

# Number of trials
nTrials = 100

function main(nTrials)
    # Spacecraft simulator
    t0      = TAIEpoch("2021-08-08T12:00:00")
    scsim   = SpacecraftSim(t0.second)

    # Initial state error variance
    σr      = 0.1   # [km]
    σv      = 0.01  # [km / s]
    σm      = 1e-6  # [kg]
    sqrtP0  = Diagonal([σr, σr, σr, σv, σv, σv, σm])
    P0      = sqrtP0.^2

    # Initial true state
    (y0true, us) = OptimalEstimationProject.GetStateAndControl(scsim, t0.second)

    # Measurement time stamps
    Δt      = 60.0        # [sec]
    gpsΔt   = 15*60     # [sec]
    ts      = range(t0.second + Δt; step = Δt, stop = t0.second + scsim.ts[2]*86400)

    # Measurement statistics
    σρ      = 10.0e-3 # [km]     Pseudorange noise standard deviation
    σr      = 5.0e-3  # [km]     GPS broadcast ephemeris standard deviation
    σa      = 1.0e-6  # [km/s^2] Accelerometer noise standard deviation

    # Process noise covariance
    Q       = Diagonal([1.0e-16, 1.0e-16, 1.0e-16, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-3])

    # GPS Simulation Span
    startWeek   = 2170
    startDay    = 0
    endWeek     = 2174
    endDay      = 6

    # Create GPS simulation object
    gpssim = GPSSim(startWeek, startDay, endWeek, endDay; σρ = σρ, σr = σr)

    # Create IMU simulation object
    imusim = IMUSim(σa)

    # MC Loop
    Threads.@threads for i in 1:nTrials
        # Initial state estimate
        xhat0   = y0true[1:7] + sqrtP0*randn(7)

        # Create UKF
        ukf = UKF(xhat0, P0, Q, (σρ^2 + 3*σr^2), σa^2, ts, gpsΔt, deepcopy(gpssim), deepcopy(imusim), deepcopy(scsim);
            lunaPerts = true, α = 1.0, β = 0.0, κ = 3.0+14.0);

        runFilter!(ukf)

        # Write results to file (Just looking at filter uncertainty quantification so only need time and error)
        df = DataFrame(
            Time    = ukf.txp,
            erx     = ukf.es[:,1],
            ery     = ukf.es[:,2],
            erz     = ukf.es[:,3],
            evx     = ukf.es[:,4],
            evy     = ukf.es[:,5],
            evz     = ukf.es[:,6],
            em      = ukf.es[:,7])

        CSV.write(datadir("UKF_MC", "trial_" * string(i) * ".csv"), df)
    end
end

main(nTrials)

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
"""
    PyObjectArray(undef, dims...)
    PyObjectArray(array)

An array of `PyObject`s which supports the Python buffer protocol.

Internally, the objects are stored as an array of pointers.
"""
mutable struct PyObjectArray{N} <: AbstractArray{PyObject,N}
    ptrs::Array{CPyPtr,N}
    function PyObjectArray{N}(::UndefInitializer, dims::NTuple{N,Integer}) where {N}
        x = new{N}(fill(CPyPtr(C_NULL), dims))
        finalizer(x) do x
            if CONFIG.isinitialized
                with_gil(false) do
                    for ptr in x.ptrs
                        C.Py_DecRef(ptr)
                    end
                end
            end
        end
    end
end
PyObjectArray{N}(::UndefInitializer, dims::Vararg{Integer,N}) where {N} =
    PyObjectArray(undef, dims)
PyObjectArray(::UndefInitializer, dims::NTuple{N,Integer}) where {N} =
    PyObjectArray{N}(undef, dims)
PyObjectArray(::UndefInitializer, dims::Vararg{Integer,N}) where {N} =
    PyObjectArray{N}(undef, dims)
PyObjectArray{N}(x::AbstractArray{T,N}) where {T,N} =
    copyto!(PyObjectArray{N}(undef, size(x)), x)
PyObjectArray(x::AbstractArray{T,N}) where {T,N} = PyObjectArray{N}(x)
export PyObjectArray

Base.IndexStyle(x) = Base.IndexStyle(x.ptrs)

Base.length(x::PyObjectArray) = length(x.ptrs)

Base.size(x::PyObjectArray) = size(x.ptrs)

Base.isassigned(x::PyObjectArray, i::Vararg{Integer}) =
    checkbounds(Bool, x.ptrs, i...) && x.ptrs[i...] != C_NULL

function Base.getindex(x::PyObjectArray, i::Integer...)
    ptr = x.ptrs[i...]
    ptr == C_NULL ? throw(UndefRefError()) : pyborrowedobject(ptr)
end

function Base.setindex!(x::PyObjectArray, _v, i::Integer...)
    vo = C.PyObject_From(_v)
    isnull(vo) && pythrow()
    C.Py_DecRef(x.ptrs[i...])
    C.Py_IncRef(vo)
    x.ptrs[i...] = vo
    x
end

function Base.deleteat!(x::PyObjectArray, i::Integer)
    C.Py_DecRef(x.ptrs[i])
    deleteat!(x.ptrs, i)
    x
end

C._pyjlarray_get_buffer(o, buf, flags, x::PyObjectArray) = C.pyjl_get_buffer_impl(
    o,
    buf,
    flags,
    pointer(x.ptrs),
    sizeof(CPyPtr),
    length(x),
    ndims(x),
    "O",
    size(x),
    strides(x.ptrs) .* sizeof(CPyPtr),
    true,
)

C._pyjlarray_array_interface(x::PyObjectArray) = C.PyDict_From(
    Dict(
        "shape" => size(x),
        "typestr" => "|O",
        "data" => (UInt(pointer(x.ptrs)), false),
        "strides" => strides(x.ptrs) .* sizeof(CPyPtr),
        "version" => 3,
    ),
)


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

struct VirtualNetworkGatewayConnectionsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a virtual network gateway connection in the specified resource group.
Param: resourceGroupName::String (required)
Param: virtualNetworkGatewayConnectionName::String (required)
Param: parameters::VirtualNetworkGatewayConnection (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualNetworkGatewayConnection
"""
function virtualNetworkGatewayConnectionsCreateOrUpdate(_api::VirtualNetworkGatewayConnectionsApi, resourceGroupName::String, virtualNetworkGatewayConnectionName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", VirtualNetworkGatewayConnection, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/connections/{virtualNetworkGatewayConnectionName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkGatewayConnectionName", virtualNetworkGatewayConnectionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes the specified virtual network Gateway connection.
Param: resourceGroupName::String (required)
Param: virtualNetworkGatewayConnectionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function virtualNetworkGatewayConnectionsDelete(_api::VirtualNetworkGatewayConnectionsApi, resourceGroupName::String, virtualNetworkGatewayConnectionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/connections/{virtualNetworkGatewayConnectionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkGatewayConnectionName", virtualNetworkGatewayConnectionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the specified virtual network gateway connection by resource group.
Param: resourceGroupName::String (required)
Param: virtualNetworkGatewayConnectionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualNetworkGatewayConnection
"""
function virtualNetworkGatewayConnectionsGet(_api::VirtualNetworkGatewayConnectionsApi, resourceGroupName::String, virtualNetworkGatewayConnectionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualNetworkGatewayConnection, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/connections/{virtualNetworkGatewayConnectionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkGatewayConnectionName", virtualNetworkGatewayConnectionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

The Get VirtualNetworkGatewayConnectionSharedKey operation retrieves information about the specified virtual network gateway connection shared key through Network resource provider.
Param: resourceGroupName::String (required)
Param: virtualNetworkGatewayConnectionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ConnectionSharedKey
"""
function virtualNetworkGatewayConnectionsGetSharedKey(_api::VirtualNetworkGatewayConnectionsApi, resourceGroupName::String, virtualNetworkGatewayConnectionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ConnectionSharedKey, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/connections/{virtualNetworkGatewayConnectionName}/sharedkey", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkGatewayConnectionName", virtualNetworkGatewayConnectionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

The List VirtualNetworkGatewayConnections operation retrieves all the virtual network gateways connections created.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualNetworkGatewayConnectionListResult
"""
function virtualNetworkGatewayConnectionsList(_api::VirtualNetworkGatewayConnectionsApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualNetworkGatewayConnectionListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/connections", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

The VirtualNetworkGatewayConnectionResetSharedKey operation resets the virtual network gateway connection shared key for passed virtual network gateway connection in the specified resource group through Network resource provider.
Param: resourceGroupName::String (required)
Param: virtualNetworkGatewayConnectionName::String (required)
Param: parameters::ConnectionResetSharedKey (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ConnectionResetSharedKey
"""
function virtualNetworkGatewayConnectionsResetSharedKey(_api::VirtualNetworkGatewayConnectionsApi, resourceGroupName::String, virtualNetworkGatewayConnectionName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", ConnectionResetSharedKey, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/connections/{virtualNetworkGatewayConnectionName}/sharedkey/reset", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkGatewayConnectionName", virtualNetworkGatewayConnectionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

The Put VirtualNetworkGatewayConnectionSharedKey operation sets the virtual network gateway connection shared key for passed virtual network gateway connection in the specified resource group through Network resource provider.
Param: resourceGroupName::String (required)
Param: virtualNetworkGatewayConnectionName::String (required)
Param: parameters::ConnectionSharedKey (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ConnectionSharedKey
"""
function virtualNetworkGatewayConnectionsSetSharedKey(_api::VirtualNetworkGatewayConnectionsApi, resourceGroupName::String, virtualNetworkGatewayConnectionName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", ConnectionSharedKey, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/connections/{virtualNetworkGatewayConnectionName}/sharedkey", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkGatewayConnectionName", virtualNetworkGatewayConnectionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export virtualNetworkGatewayConnectionsCreateOrUpdate, virtualNetworkGatewayConnectionsDelete, virtualNetworkGatewayConnectionsGet, virtualNetworkGatewayConnectionsGetSharedKey, virtualNetworkGatewayConnectionsList, virtualNetworkGatewayConnectionsResetSharedKey, virtualNetworkGatewayConnectionsSetSharedKey


text:
# lstm weights initialization
# w[2k-1], w[2k] : weight and bias for kth layer respectively
function initweights(atype, hiddens, embedding, vocab, init=xavier)
    weights = Array(Any, 2length(hiddens))
    input = embedding
    for k = 1:length(hiddens)
        weights[2k-1] = init(input+hiddens[k], 4hiddens[k])
        weights[2k] = zeros(1, 4hiddens[k])
        weights[2k][1:hiddens[k]] = 1 # forget gate bias
        input = hiddens[k]
    end
    return map(w->convert(atype, w), weights)
end

# state initialization
# s[2k-1], s[2k] : hidden and cell respectively
function initstate(atype, hiddens, batchsize)
    state = Array(Any, 2length(hiddens))
    for k=1:length(hiddens)
        state[2k-1] = atype(zeros(batchsize, hiddens[k]))
        state[2k] = atype(zeros(batchsize, hiddens[k]))
    end
    return state
end

function initmodel(atype, hiddens, embedding, vocabsize, init=xavier)
    model = Dict{Symbol, Any}()
    model[:forw] = initweights(atype, hiddens, embedding, vocabsize, init)
    model[:back] = initweights(atype, hiddens, embedding, vocabsize, init)
    model[:fembed] = atype(init(vocabsize, embedding))
    model[:bembed] = atype(init(vocabsize, embedding))
    model[:soft] = [ atype(init(2hiddens[end], vocabsize)), atype(init(1, vocabsize)) ]
    return model
end

function lstm(weight,bias,hidden,cell,input)
    gates   = hcat(input,hidden) * weight .+ bias
    hsize   = size(hidden,2)
    forget  = sigm(gates[:,1:hsize])
    ingate  = sigm(gates[:,1+hsize:2hsize])
    outgate = sigm(gates[:,1+2hsize:3hsize])
    change  = tanh(gates[:,1+3hsize:end])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh(cell)
    return (hidden,cell)
end

# multilayer lstm forward, returns the final hidden
function forward(weight, states, input)
    x = input
    for i=1:2:length(states)
        (states[i], states[i+1]) = lstm(weight[i], weight[i+1], states[i], states[i+1], x)
        x = states[i]
    end
    return x
end

# dropout layer
function m_dropout(x, p)
    if p > 0
        return x .* (rand!(similar(AutoGrad.getval(x))) .> p) * (1/(1-p))
    else
        return x
    end
end

function logprob(output, ypred)
    nrows,ncols = size(ypred)
    index = similar(output)
    @inbounds for i=1:length(output)
        index[i] = i + (output[i]-1)*nrows
    end
    o1 = logp(ypred,2)
    o2 = o1[index]
    o3 = sum(o2)
    return o3
end

function bilstm(model, states, sequence; pdrop=(0, 0))
    total = 0.0
    count = 0
    atype = typeof(AutoGrad.getval(model[:fembed]))

    # forward lstm
    fhiddens = Array(Any, length(sequence))
    sf = copy(states)
    for i=1:length(sequence)-1
        x = model[:fembed][sequence[i], :]
        h = forward(model[:forw], sf, x)
        fhiddens[i+1] = copy(h)
    end
    fhiddens[1] = convert(atype, zeros(size(fhiddens[2])))

    # backward lstm
    bhiddens = Array(Any, length(sequence))
    sb = copy(states)
    for i=length(sequence):-1:2
        x = model[:bembed][sequence[i], :]
        h  = forward(model[:back], sb, x)
        bhiddens[i-1] = copy(h)
    end
    bhiddens[end] = convert(atype, zeros(size(bhiddens[2])))

    # concatenate layer
    for i=1:length(fhiddens)
        hf = m_dropout(fhiddens[i], pdrop[1])
        hb = m_dropout(bhiddens[i], pdrop[2])
        ypred = hcat(hf, hb) * model[:soft][1] .+ model[:soft][2]
        total += logprob(sequence[i], ypred)
        count += length(sequence[i])
    end
    return - total / count
end

bilstmgrad = grad(bilstm)

oparams{T<:Number}(::KnetArray{T},otype; o...)=otype(;o...)
oparams{T<:Number}(::Array{T},otype; o...)=otype(;o...)
oparams(a::Associative,otype; o...)=Dict(k=>oparams(v,otype;o...) for (k,v) in a)
oparams(a,otype; o...)=map(x->oparams(x,otype;o...), a)

function train(model, state, sequence, opts; pdrop=(0,0))
    gloss = bilstmgrad(model, state, sequence; pdrop=pdrop)
    update!(model, gloss, opts)
end

convertmodel{T<:Number}(x::KnetArray{T}) = convert(Array{T}, x)
convertmodel{T<:Number}(x::Array{T}) = convert(Array{T}, x)
convertmodel(a::Associative)=Dict(k=>convertmodel(v) for (k,v) in a)
convertmodel(a) = map(x->convertmodel(x), a)

function devperp(m, s, dev)
    devloss = 0
    for d in dev
        devloss += bilstm(m, s, d)
    end
    return exp(devloss / length(dev))
end

max_stars_repo_name: kirnap/wordsubs

max_issues_repo_name: kirnap/wordsubs

max_forks_repo_name: kirnap/wordsubs


text:
xdoc = XMLDocument()

xroot = create_root(xdoc, "States")

add_pi(xroot, "State", "Massachusetts")
add_pi(xroot, "State", "New Jersey")
add_pi(xroot, "State", "New York")

rtxt = """
<?xml version="1.0" encoding="utf-8"?>
<States>
  <?State Massachusetts?>
  <?State New Jersey?>
  <?State New York?>
</States>
"""

@test strip(string(xdoc)) == strip(rtxt)

free(xdoc)


text:
module remesh_smoothing_examples
using FinEtools
using FinEtools.MeshExportModule
using FinEtools.TetRemeshingModule

function remesh1()
    L= 0.3;
    W = 0.3;
    a = 0.15;
    nL=46; nW=46; na=36;

    fens,fes = T4block(a,L,W,nL,nW,na,:a);
    t = deepcopy(connasarray(fes));
    v = deepcopy(fens.xyz);
    tmid = ones(Int, size(t,1));

    desired_ts =a;
    bfes = meshboundary(fes);
    f = connectednodes(bfes);
    bv = zeros(Bool, size(v,1));
    bv[f] = true;

    println("Mesh size: initial = $(size(t,1))")
    t0 = time()

    t, v, tmid = TetRemeshingModule.coarsen(t, v, tmid; bv = bv, desired_ts = desired_ts);

    println("Mesh size: final = $(size(t,1)) [$(time() - t0) sec]")

    fens.xyz = deepcopy(v)
    fes = fromarray!(fes, t)
    setlabel!(fes, tmid)
    geom  =  NodalField(fens.xyz)

    femm  =  FEMMBase(IntegData(fes, SimplexRule(3, 1)))
    V = integratefunction(femm, geom, (x) ->  1.0)
    println("V = $(V) compared to $(L * W * a)")

    # File = "test1.vtk"
    # MeshExportModule.vtkexportmesh(File, t, v, MeshExportModule.T4)
    # @async run(`"paraview.exe" $File`)
end # remesh1

function mesh_smoothing()
    println("""
    Meshing, deforming  and smoothing
    """)

    A = 100. # strip width
    N = 16
    tolerance = A / N / 1.0e5
    fens, fes = T3block(A, A, N, N)

    bnl = connectednodes(meshboundary(fes))
    for ixxxx = 1:length(bnl)
        x, y = fens.xyz[bnl[ixxxx], :]
        fens.xyz[bnl[ixxxx], 1] += A / N * sin(2 * pi * y / A)
        fens.xyz[bnl[ixxxx], 2] += -A / N * sin(2 * pi * x / A)
    end

    File =  "mesh_smoothing_before.vtk"
    vtkexportmesh(File, fens, fes);
    @async run(`"paraview.exe" $File`)
    println("$(fens.xyz[Int(N^2 / 2), :] )")

    fixedv = falses(count(fens))
    fixedv[bnl] = true
    fens = meshsmoothing(fens, fes; fixedv=fixedv, method=:taubin, npass=100)

    println("$(fens.xyz[Int(N^2 / 2), :] )")

    geom = NodalField(fens.xyz)

    File =  "mesh_smoothing_after.vtk"
    vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.T3);
    @async run(`"paraview.exe" $File`)

    println("Done")
    true

end # mesh_smoothing

function allrun()
    println("#####################################################")
    println("# remesh1 ")
    remesh1()
    println("#####################################################")
    println("# mesh_smoothing ")
    mesh_smoothing()
    return true
end # function allrun

end # module remesh_smoothing_examples


text:
using DataFrames
using DataArrays
using DifferentialEquations
using ParameterizedFunctions
using Gadfly

# definition of differential equations
g = @ode_def_nohes LotkaVolterra begin
         dx = A*x - B*x*y
         dy = -C*y + D*x*y
end A=>1.0 B =>1.0 C =>1.0 D =>1.0

# solving differential equations
function lotka_volterra(a, b, c, d, t0, t1, x0, y0)
  u0 = [x0;y0]
  tspan = (t0,t1)
  h = LotkaVolterra(A=a,B=b,C=c,D=d)
  prob = ODEProblem(h,u0,tspan)
  sol = solve(prob, RK4(), dt=0.01)
  return sol
end

# extracting results and saving into array
function into_array(sol, expNum)
  results = Array{Any,1}[]
  mintime = first(sol.t)
  maxtime = last(sol.t)
  i = mintime
  push!(results, ["t", "x", "y", "experiment"])
  while i < maxtime
    push!(results, [i, sol(i)[1], sol(i)[2], "exp$expNum"])
    i += 0.1
  end
  return results
end

# saving results (array) into csv file
function into_file(filename, results)
  open(filename,"w")
  writecsv(filename, results)
end

function experiments()
  s1 = lotka_volterra(1.5,1.0,3.0,1.0,0.0,5.0,8.0,4.0)
  s2 = lotka_volterra(1.5,1.0,3.0,1.0,0.0,5.0,16.0,2.0)
  s3 = lotka_volterra(1.5,1.0,3.0,1.0,0.0,5.0,12.0,4.0)
  s4 = lotka_volterra(1.5,1.0,3.0,1.0,0.0,5.0,2.0,10.0)
  into_file("res1.csv",into_array(s1,1))
  into_file("res2.csv",into_array(s2,2))
  into_file("res3.csv",into_array(s3,3))
  into_file("res4.csv",into_array(s4,4))
  return [s1,s2,s3,s4]
end

# read in csv files into DataFrame
function read_in(res...)
  df = vcat([readtable(res[i]) for i in 1:length(res)])
end

# print out min, max and mean value for each experiment
function min_max_mean(data)
  experiments = sort(unique(data[:experiment]))
  i = 1
  for exp in experiments
    println("Experiment $i\n")
    exp_data = data[data[:experiment] .== exp, :]
    println("Prey population\n")
    println("Minimum: $(minimum(exp_data[:x]))")
    println("Maximum: $(maximum(exp_data[:x]))")
    println("Mean:    $(mean(exp_data[:x]))")
    println("\nPredator population\n")
    println("Minimum: $(minimum(exp_data[:y]))")
    println("Maximum: $(maximum(exp_data[:y]))")
    println("Mean:    $(mean(exp_data[:y]))")
    println("\n")
    i += 1
  end
end

# adding new column
function diff_column(data)
  new_data = data
  new_data[:diff] = new_data[:y] - new_data[:x]
  return new_data
end

# plotting functions
function plot_prey_population(sol)
  df=DataFrame(t=sol.t, prey=map(x->x[1],sol.u))
  plot(df,  x="t", y="prey") #prey population
end

function plot_predator_population(sol)
  df2=DataFrame(t=sol.t, predator=map(y->y[2],sol.u))
  plot(df2,  x="t", y="predator") #predator population
end

function plot_prey_predator(data)
  dataframe = DataFrame()
  exp1_data = data[data[:experiment] .== "exp1", :]
  exp2_data = data[data[:experiment] .== "exp2", :]
  exp3_data = data[data[:experiment] .== "exp3", :]
  exp4_data = data[data[:experiment] .== "exp4", :]
  dataframe[:x1] = exp1_data[:x]
  dataframe[:y1] = exp1_data[:y]
  dataframe[:x2] = exp2_data[:x]
  dataframe[:y2] = exp2_data[:y]
  dataframe[:x3] = exp3_data[:x]
  dataframe[:y3] = exp3_data[:y]
  dataframe[:x4] = exp4_data[:x]
  dataframe[:y4] = exp4_data[:y]

  plot(dataframe, layer(x="x1", y="y1", Geom.point, Theme(default_color= colorant"blue")),
                  layer(x="x2", y="y2", Geom.point, Theme(default_color= colorant"red")),
                  layer(x="x3", y="y3", Geom.point, Theme(default_color= colorant"green")),
                  layer(x="x4", y="y4", Geom.point, Theme(default_color= colorant"yellow")),
                  Guide.XLabel("Prey population"),
                  Guide.YLabel("Predator population"),
                  Guide.Title("Lotka - Volterra"))
end

solutions = experiments()

solutions_data = read_in("res1.csv", "res2.csv", "res3.csv", "res4.csv")

min_max_mean(solutions_data)

diff_column(solutions_data)

plot_prey_population(solutions[1])

plot_predator_population(solutions[1])

plot_prey_predator(solutions_data)

max_stars_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_issues_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_forks_repo_name: wchmielarz/NobodyDoesItLikeJuliet


text:
VM = Union{Vector{Float64},Matrix{Float64}}

immutable State_lm
  b::Vector{Float64}
  sig2::Float64
end

immutable LM
  post_params::Vector{State_lm}
  DIC::Float64
end

immutable Summary_lm
  mean_beta::Vector{Float64}
  std_beta::Vector{Float64}
  quantile_beta::Matrix{Float64}
  mean_sig::Float64
  std_sig::Float64
  quantile_sig::Vector{Float64}
  DIC::Float64
  B::Int
end

function show(io::IO, SL::Summary_lm)
  const P = length(SL.mean_beta)
  zeroInCI = [SL.quantile_beta[2,p] <= 0 <= SL.quantile_beta[6,p] ? "" : "*" for p in 1:P]

  @printf "%5s%10s%10s%3s\n" "" "mean" "std" "≠0"
  for k in 1:P
    @printf "%5s%10.4f%10.4f%3s\n" string("β",k-1) SL.mean_beta[k] SL.std_beta[k] zeroInCI[k]
  end
  @printf "%5s%10.4f%10.4f\n" "σ" SL.mean_sig SL.std_sig
  @printf "%5s%10.2f\n" "DIC" SL.DIC
end

function dic{T <: VM}(post_samples::Vector{State_lm},
                      y::Vector{Float64}, X::T;
                      addIntercept=true)
  if addIntercept
    X::Matrix{Float64} = [ones(size(X,1)) X]
  end
  const (N,P) = size(X)

  function loglike(param)
    -N/2 * log(2pi*param.sig2) - sum((y - X*param.b).^2) / (2param.sig2)
  end

  dic(post_samples,loglike)
end

function cpo(model::LM, y::Vector{Float64}, X::Matrix{Float64})
  function den(p::State_lm, yi::Float64, xi::Vector{Float64})
    xb = xi'p.b
    assert(length(xb) == 1)
    # Declaring the types really speeds things up here
    return pdf(Normal(xb[1]::Float64, sqrt(p.sig2)::Float64),
               yi::Float64)
  end
  cpo(model.post_params, den, y, X)
end

function summary(out::LM; alpha::Float64=.05)
  const post_beta= hcat(map(o -> o.b, out.post_params)...)'
  const (B,P) = size(post_beta)
  const mean_beta = vec(mean(post_beta,1))
  const std_beta = vec(std(post_beta,1))
  const quantile_beta = hcat([ quantile(post_beta[:,p], quants) for p in 1:P]...)
  const post_sig = map(o -> sqrt(o.sig2), out.post_params)
  const mean_sig = mean(post_sig)
  const std_sig = std(post_sig)
  const quantile_sig = quantile(post_sig, quants)

  return Summary_lm(mean_beta,std_beta,quantile_beta,
                    mean_sig,std_sig,quantile_sig,
                    out.DIC,B)
end

function lm{T <: VM}(y::Vector{Float64}, X::T; B::Int=10000, burn::Int=10,
                     addIntercept::Bool=true)
  if addIntercept
    X = [ones(size(X,1)) X]
  end
  const (N,K) = size(X)
  const XXi = inv(X'X)
  const a = (N-K) / 2
  const BETA_HAT = XXi * X'y # use QR to do this instead?
  const beta_init = fill(0.0,K)
  const s2_init = 1.0

  function update(state::State_lm)
    const s2_new = rig(a, sum((y-X*state.b).^2)/2)
    const b_new =
      rand(MultivariateNormal(BETA_HAT::Vector{Float64},
                              s2_new*XXi::Matrix{Float64}))
    State_lm(b_new,s2_new)
  end

  const post_params = MCMC.gibbs(State_lm(beta_init,s2_init), update, B, burn)
  const DIC = dic(post_params, y, X, addIntercept=false)

  LM(post_params,DIC)
end

#= For general covariance matrix
function lm(y:: Vector, X::Matrix, V::Matrix)
end
=#

#=
using BayesLM
n = 100
X = randn(n)
b = [0,3]
y = [ones(n) X] * b + randn(n)
@time model = lm(y,X);
summary(model)

using Distributions
@time cpo(model.post_params, cpo_f, [ones(n) X])
=#


text:
__precompile__()

"""
*MultipleTesting* package

* Adjusting p-values for multiple testing
* Estimating the fraction π₀ of tests under the null hypothesis
* Combination of p-values
"""
module MultipleTesting

using StatsBase
import StatsBase: fit

using Distributions
import Distributions: estimate

import SpecialFunctions: digamma

import Random: shuffle!

export
    PValues,
    adjust,
    PValueAdjustment,
    Bonferroni,
    BenjaminiHochberg,
    BenjaminiHochbergAdaptive,
    BenjaminiYekutieli,
    BenjaminiLiu,
    Hochberg,
    Holm,
    Hommel,
    Sidak,
    ForwardStop,
    BarberCandes,
    estimate_pi0,
    Pi0Estimator,
    Storey,
    StoreyBootstrap,
    LeastSlope,
    Oracle,
    TwoStep,
    RightBoundary,
    CensoredBUM,
    CensoredBUMFit,
    BUM,
    BUMFit,
    FlatGrenander,
    ConvexDecreasing,
    ConvexDecreasingFit,
    isin,
    fit,
    BetaUniformMixtureModel,
    PValueCombination,
    combine,
    FisherCombination,
    LogitCombination,
    StoufferCombination,
    TippettCombination,
    SimesCombination,
    WilkinsonCombination,
    MinimumCombination,
    estimate,
    HigherCriticismScores,
    HigherCriticismThreshold

include("types.jl")
include("utils.jl")
include("pval-adjustment.jl")
include("pi0-estimators.jl")
include("model.jl")
include("combinations.jl")
include("higher-criticism.jl")

end


text:
using OrdinaryDiffEq
using Trixi
using Plots

# Vgl. Gauss on an cartesian Grid with coordinates (-1,-1),(1,1) with periodic BC
###############################################################################
CFL = 0.01
tspan = (0.0, 10)
N = 3
c = 8
mu = 0.01

function boundary_condition_constant( u_inner, q1, q2, orientation, direction, x, t,
    surface_flux_function,
    equations::Trixi.AbstractEquations)
    # Far Field Conditions
    u_boundary = initial_condition(x , t, equations)

    # Calculate boundary flux
    if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
        flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
    else # direction == 4 # u_boundary is "left" of boundary, u_inner is "right" of boundary
        flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
    end

    return flux
end

function initial_condition_density_pulse_new(x, t, equations::Trixi.AbstractEquations)
    rho = 1 + exp(-(x[1]^2 + x[2]^2)) / 2
    v1 = 0
    v2 = 0
    rho_v1 = rho * v1
    rho_v2 = rho * v2
    p = 1
    rho_e = p / (1.4 - 1) + 1 / 2 * rho * (v1^2 + v2^2)
    return SVector(rho, rho_v1, rho_v2, rho_e)
end

function initial_condition_point(x, t, equations::CompressibleEulerEquations2D)
    if x[1] <= 0 && x[2] <= 0
        rho = 5
        v1 = 0
        v2 = 0
        rho_v1 = rho * v1
        rho_v2 = rho * v2
        p = 1
        rho_e = p / (equations.gamma - 1) + 1 / 2 * rho * (v1^2 + v2^2)
    else
        rho = 1
        v1 = 0
        v2 = 0
        rho_v1 = rho * v1
        rho_v2 = rho * v2
        p = 1
        rho_e = p / (equations.gamma - 1) + 1 / 2 * rho * (v1^2 + v2^2)
    end
    return SVector(rho, rho_v1, rho_v2, rho_e)
end

equations = CompressibleEulerEquations2D(1.4, viscous=true, mu=mu)

initial_condition = initial_condition_density_pulse_new

boundary_conditions = boundary_condition_constant

# boundary_conditions = (x_neg=boundary_condition_periodic,
#                        x_pos=boundary_condition_periodic,
#                        y_neg=boundary_condition_periodic,
#                        y_pos=boundary_condition_periodic)

cells_per_dimension = (c, c)
mesh = CurvedMesh(cells_per_dimension, (-1.0, -1.0), (1.0, 1.0))

surface_flux = flux_lax_friedrichs # FluxPlusDissipation(flux_chandrashekar, DissipationLocalLaxFriedrichs(max_abs_speed_naive))
basis = LobattoLegendreBasis(N)
volume_flux  = flux_chandrashekar
volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
solver = DGSEM(basis, surface_flux, volume_integral)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_condition_constant)

###############################################################################
# ODE solvers, callbacks etc.

ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 500
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
                                      # extra_analysis_errors=(:conservation_error,),
                                      # extra_analysis_integrals=(entropy, energy_total,
                                      # energy_kinetic, energy_internal)
                                      )

alive_callback = AliveCallback(analysis_interval=analysis_interval)
# visualization = VisualizationCallback(interval=28, plot_creator=Trixi.save_plot)
visualization = VisualizationCallback(interval=500, variable_names=["rho"], plot_creator=Trixi.save_plot, clims=(1.0, 1.5))

save_solution = SaveSolutionCallback(interval=50,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=CFL)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        visualization,
                        #save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

max_stars_repo_path: src/WR2/MS5 Gauss.jl

max_issues_repo_path: src/WR2/MS5 Gauss.jl

max_forks_repo_path: src/WR2/MS5 Gauss.jl


text:
using Test

@assert 1 == 1

@assert 1 == 2

@assert 1 == 2 "Very Bad Error"

@macroexpand @assert 1 == 2 "Very Bad Error"

@macroexpand @assert 1 == 2

@test 1 == 1

@testset "Helpful Tests" begin
    @testset "Tautologies" begin
        @test 1 == 1
        @test "a" == "a"
        @test pi == 3.1415
    end
    @testset "Hopefully True" begin
        @test isapprox(2, 2.0002, atol = 1)
        @test isapprox(2 + 2, 5, atol = 1)
        @test isapprox(2 + 2, 3)
        @test isapprox(1e-10, 2e-9, atol = 1e-8)
    end
end

max_stars_repo_name: phrb/PDP

max_issues_repo_name: phrb/PDP

max_forks_repo_name: phrb/PDP


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

doc"""
    AbstractProposalDist

The following functions must be implemented for subtypes:

* `BAT.distribution_logpdf`
* `BAT.proposal_rand!`
* `BAT.nparams`, returning the number of parameters (i.e. dimensionality).
* `Base.issymmetric`, indicating whether p(a -> b) == p(b -> a) holds true.

In some cases, it may be desirable to override the default implementation
of `BAT.distribution_logpdf!`.
"""
abstract type AbstractProposalDist end
export AbstractProposalDist

doc"""
    distribution_logpdf(
        pdist::AbstractProposalDist,
        params_new::AbstractVector,
        params_old:::AbstractVector
    )

Analog to `distribution_logpdf!`, but for a single parameter vector.
"""
function distribution_logpdf end
export distribution_logpdf

# TODO: Implement distribution_logpdf for included proposal distributions

doc"""
    distribution_logpdf!(
        p::AbstractArray,
        pdist::AbstractProposalDist,
        params_new::AbstractVecOrMat,
        params_old:::AbstractVecOrMat
    )

Returns log(PDF) value of `pdist` for transitioning from old to new parameter
values for multiple parameter sets.

end

Input:

* `params_new`: New parameter values (column vectors)
* `params_old`: Old parameter values (column vectors)

Output is stored in

* `p`: Array of PDF values, length must match, shape is ignored

Array size requirements:

* `size(params_old, 1) == size(params_new, 1) == length(pdist)`
* `size(params_old, 2) == size(params_new, 2)` or `size(params_old, 2) == 1`
* `size(params_new, 2) == length(p)`

Implementations of `distribution_logpdf!` must be thread-safe.
"""
function distribution_logpdf! end
export distribution_logpdf!

# TODO: Default implementation of distribution_logpdf!

doc"""
    function proposal_rand!(
        rng::AbstractRNG,
        pdist::GenericProposalDist,
        params_new::AbstractVecOrMat,
        params_old::AbstractVecOrMat
    )

Generate one or multiple proposed parameter vectors, based on one or multiple
previous parameter vectors.

Input:

* `rng`: Random number generator to use
* `pdist`: Proposal distribution to use
* `params_old`: Old parameter values (vector or column vectors, if a matrix)

Output is stored in

* `params_new`: New parameter values (vector or column vectors, if a matrix)

The caller must guarantee:

* `size(params_old, 1) == size(params_new, 1)`
* `size(params_old, 2) == size(params_new, 2)` or `size(params_old, 2) == 1`
* `params_new !== params_old` (no aliasing)

Implementations of `proposal_rand!` must be thread-safe.
"""
function proposal_rand! end
export proposal_rand!

struct GenericProposalDist{D<:Distribution{Multivariate},SamplerF,S<:Sampleable} <: AbstractProposalDist
    d::D
    sampler_f::SamplerF
    s::S

    function GenericProposalDist{D,SamplerF}(d::D, sampler_f::SamplerF) where {D<:Distribution{Multivariate},SamplerF}
        s = sampler_f(d)
        new{D,SamplerF, typeof(s)}(d, sampler_f, s)
    end

end

export GenericProposalDist

GenericProposalDist(d::D, sampler_f::SamplerF) where {D<:Distribution{Multivariate},SamplerF} =
    GenericProposalDist{D,SamplerF}(d, sampler_f)

GenericProposalDist(d::Distribution{Multivariate}) = GenericProposalDist(d, bat_sampler)

GenericProposalDist(D::Type{<:Distribution{Multivariate}}, n_params::Integer, args...) =
    GenericProposalDist(D, Float64, n_params, args...)

Base.similar(q::GenericProposalDist, d::Distribution{Multivariate}) =
    GenericProposalDist(d, q.sampler_f)

function Base.convert(::Type{AbstractProposalDist}, q::GenericProposalDist, T::Type{<:AbstractFloat}, n_params::Integer)
    n_params != nparams(q) && throw(ArgumentError("q has wrong number of parameters"))
    q
end

get_cov(q::GenericProposalDist) = get_cov(q.d)
set_cov(q::GenericProposalDist, Σ::AbstractMatrix{<:Real}) = similar(q, set_cov(q.d, Σ))

function distribution_logpdf!(
    p::AbstractArray,
    pdist::GenericProposalDist,
    params_new::AbstractMatrix,
    params_old::AbstractVecOrMat
)
    params_diff = params_new .- params_old # TODO: Avoid memory allocation
    Distributions.logpdf!(p, pdist.d, params_diff)
end

function distribution_logpdf!(
    p::AbstractArray,
    pdist::GenericProposalDist,
    params_new::AbstractVector,
    params_old::AbstractVector
)
    p[1] = distribution_logpdf(pdist, params_new, params_old)
    p
end

function distribution_logpdf(
    pdist::GenericProposalDist,
    params_new::AbstractVector,
    params_old::AbstractVector
)
    params_diff = params_new .- params_old # TODO: Avoid memory allocation
    Distributions.logpdf(pdist.d, params_diff)
end

function proposal_rand!(
    rng::AbstractRNG,
    pdist::GenericProposalDist,
    params_new::AbstractVecOrMat,
    params_old::AbstractVecOrMat
)
    rand!(rng, pdist.s, params_new)
    params_new .+= params_old
end

nparams(pdist::GenericProposalDist) = length(pdist.d)

Base.issymmetric(pdist::GenericProposalDist) = issymmetric_around_origin(pdist.d)

struct GenericUvProposalDist{D<:Distribution{Univariate},T<:Real,SamplerF,S<:Sampleable} <: AbstractProposalDist
    d::D
    scale::Vector{T}
    sampler_f::SamplerF
    s::S
end

export GenericUvProposalDist

GenericUvProposalDist(d::Distribution{Univariate}, scale::Vector{<:AbstractFloat}, samplerF) =
    GenericUvProposalDist(d, scale, samplerF, samplerF(d))

GenericUvProposalDist(d::Distribution{Univariate}, scale::Vector{<:AbstractFloat}) =
    GenericUvProposalDist(d, scale, bat_sampler)

BAT.nparams(pdist::GenericUvProposalDist) = size(pdist.scale, 1)

Base.issymmetric(pdist::GenericUvProposalDist) = issymmetric_around_origin(pdist.d)

function BAT.distribution_logpdf(
    pdist::GenericUvProposalDist,
    params_new::AbstractVecOrMat,
    params_old::AbstractVecOrMat
)
    params_diff = (params_new .- params_old) ./ pdist.scale  # TODO: Avoid memory allocation
    sum_first_dim(Distributions.logpdf.(pdist.d, params_diff))  # TODO: Avoid memory allocation
end

function BAT.proposal_rand!(
    rng::AbstractRNG,
    pdist::GenericUvProposalDist,
    params_new::AbstractVector,
    params_old::AbstractVector
)
    params_new .= params_old
    dim = rand(rng, eachindex(pdist.scale))
    params_new[dim] += pdist.scale[dim] * rand(rng, pdist.s)
    params_new
end

abstract type ProposalDistSpec end

export ProposalDistSpec

struct MvTDistProposalSpec <: ProposalDistSpec
    df::Float64
end

export MvTDistProposalSpec

MvTDistProposalSpec() = MvTDistProposalSpec(1.0)

(ps::MvTDistProposalSpec)(T::Type{<:AbstractFloat}, n_params::Integer) =
    GenericProposalDist(MvTDist, T, n_params, convert(T, ps.df))

function GenericProposalDist(::Type{MvTDist}, T::Type{<:AbstractFloat}, n_params::Integer, df = one(T))
    Σ = PDMat(full(ScalMat(n_params, one(T))))
    zeromean = true
    μ = fill(zero(T), n_params)
    M = typeof(Σ)
    d = Distributions.GenericMvTDist{T,M}(convert(T, df), n_params, zeromean, μ, Σ)
    GenericProposalDist(d)
end

struct UvTDistProposalSpec <: ProposalDistSpec
    df::Float64
end

export UvTDistProposalSpec

(ps::UvTDistProposalSpec)(T::Type{<:AbstractFloat}, n_params::Integer) =
    GenericUvProposalDist(TDist(convert(T, ps.df)), fill(one(T), n_params))


text:
path=dirname(@__FILE__)
f = load("$path/modelsolutions.jld")

NKsol=solve(twoEQNK)
FRWZsol=solve(FRWZ)

@testset "Number of solutions" begin
    @test length(NKsol)==1
    @test length(FRWZsol) == 2
end

@testset "MX matrices" begin
    @test isapprox(NKsol[1].MX, f["NKsolMX"])
    @test isapprox(FRWZsol[1].MX, f["FRWZsolMX1"])
    @test isapprox(FRWZsol[2].MX, f["FRWZsolMX2"])
end

@testset "ME matrices" begin
    @test isapprox(NKsol[1].ME, f["NKsolME"])
    @test isapprox(FRWZsol[1].ME, f["FRWZsolME1"])
    @test isapprox(FRWZsol[2].ME, f["FRWZsolME2"])
end

@testset "MC matrices" begin
    @test isapprox(NKsol[1].MC, zeros(Float64,2,1,2))
    @test isapprox(FRWZsol[1].MC, f["FRWZsolMC1"])
    @test isapprox(FRWZsol[2].MC, f["FRWZsolMC2"])
end


text:
using MathProgComplex, ProfileView

function main()

    instance_path = getinstancepath("Matpower", "QCQP", "case9")
    ## Build real problem
    problem_C, point = import_from_dat(instance_path)
    problem = pb_cplx2real(problem_C)

    pbsolved = :MomentRelaxation

    ## Build relaxation context
    relax_ctx = set_relaxation(problem; hierarchykind=hierarchykind,
                                        d=d,
                                        symmetries=symmetries,
                                        params = Dict(:opt_outlev=>1,
                                                      :opt_outmode=>0,
                                                      :opt_outcsv=>0,
                                                    #   :opt_msk_maxtime=>2*3600,
                                                      :opt_pbsolved=>pbsolved,
                                                      :opt_solver=>:None,
                                                      :opt_outname=>joinpath(logpath, "momentsos.log"),
                                                      :opt_outcsvname=>joinpath(logpath, "momentsos_solve.csv")))

    max_cliques = get_maxcliques(relax_ctx, problem)

    momentmat_param, localizingmat_param = build_sparsity(relax_ctx, problem, max_cliques)

    mmtrel_pb = build_momentrelaxation(relax_ctx, problem, momentmat_param, localizingmat_param, max_cliques);

    sdp_pb = build_SDP_Instance_from_SDPDual(mmtrel_pb);

    Profile.clear()
    @profile build_momentrelaxation(relax_ctx, problem, momentmat_param, localizingmat_param, max_cliques)
    ProfileView.view()

end


text:
module FEMAssembler

    # initialization of structures
    include("init.jl")

    # elements
    include("elements/linearTrussStiffnessMatrix.jl")
    include("elements/linearTetraSolid.jl")

    include("mesh/dofsComputations.jl")
    include("assembler.jl")

    # exports
    include("exports.jl")

end # module


text:
using AlgebraOfGraphics, CategoricalArrays, CairoMakie, DataFrames
using LaTeXStrings, HypergeometricFunctions, SpecialFunctions
set_aog_theme!()
update_theme!(font = "Arial", fontsize = 28)
theme_minimal
function weight(j, s, xi)
    exp(
        logabsbinomial(j + s - 2, j - 1)[1] -
        log(j) +
        s * log(xi) +
        (j - 1) * log(1 - xi) +
        log(_₂F₁(j + s - 1, 1, j + 1, 1 - xi))
    )
end

df_j = DataFrame(j = 1:20);
df_ϕ = DataFrame(ϕ = [0.1, 0.3, 0.5]);
df_s = DataFrame(s = [1, 2, 3, 4]);
df = crossjoin(df_j, df_ϕ, df_s);
df[!, :w] = weight.(df[!, :j], df[!, :s], df[!, :ϕ]);
df[!, :s] = string.(df[!, :s]);
df[!, :ϕ] = string.(df[!, :ϕ]);
df[!, :ϕ] = LaTeXString.("φ(x_i ⋅ β) = " .* df[!, :ϕ]);
plt = data(df) *
    visual(Lines) *
    mapping(:j, :w, col = :ϕ, color = :s)
axis = (xlabel = L"j", ylabel = L"w_j(x)", color = L"s");
fg = draw(plt; axis, figure = (resolution = (900, 450), px_per_unit = 8))
fg
save("figures/final/lineplot-dgsb-weights.pdf", fg, px_per_unit = 3)


text:
for (T1, T2) in ((:T, :(float(T))), (:(Union{Missing,T}), :(Union{Missing,float(T)})))
  @eval begin

    # unweighted windowed function application that tapers

    function running(fun::Function, data::V, windowspan::Int) where {T, V<:AbstractVector{$T1}}
        ndata   = length(data)
        nvals   = nrolled(ndata, windowspan)
        ntapers = ndata - nvals

        result = zeros($T2, ndata)

        result[1:ntapers] = tapers(fun, data[1:ntapers])
        ntapers += 1
        result[ntapers:ndata] = rolling(fun, data, windowspan)

        return result
    end

    # weighted windowed function application that tapers

    function running(fun::Function, data::V, windowspan::Int, weighting::F) where
                     {T, V<:AbstractVector{$T1}, N<:Number, F<:Vector{N}}

        ndata   = length(data)
        nvals   = nrolled(ndata, windowspan)
        ntapers = ndata - nvals

        result = zeros($T2, ndata)

        result[1:ntapers] = tapers(fun, data[1:ntapers], weighting[end-(ntapers-1):end])
        ntapers += 1
        result[ntapers:ndata] = rolling(fun, data, windowspan, weighting)

        return result
    end

    running(fun::Function, data::V, windowspan::Int, weighting::W) where
                    {T, V<:AbstractVector{$T1}, N<:Number, W<:AbstractWeights} =
        running(fun, data, windowspan, weighting.values)

  end
end


text:
using Luxor, Test

function polyconvextest()

    Drawing(100, 100, :png)
    origin()

    # polygons are convex
    for i in 3:16
        @test ispolyconvex(ngon(O, 100, i, vertices=true)) == true
    end

    # 4+ pointed stars aren't convex
    # 3 pointed star is effectively a triangle
    for i in 4:16
        @test ispolyconvex(star(O, 100, i, 0.5, 0, vertices=true)) == false
    end

    pts = [Point(x, 0) for x in 0:100]
    @test ispolyconvex(pts) == true

    @test finish() == true
end

function polyclockwisetest()

    Drawing(100, 100, :png)
    origin()

    # polygons are clockwise
    for i in 3:16
        @test ispolyclockwise(ngon(O, 100, i, vertices=true)) == true
    end
    # stars are too
    for i in 3:16
        @test ispolyclockwise(star(O, 100, i, 0.5, 0, vertices=true)) == true
    end

    pts = [Point(x, 0) for x in 0:100]
    @test ispolyclockwise(pts) == false

    pts = [Point(10cos(x), 10sin(x)) for x in 0:100]
    @test ispolyclockwise(pts) == true

    @test finish() == true
end

polyconvextest()
polyclockwisetest()

println("...finished convex and clockwise tests")


text:
using LinearAlgebra, ForwardDiff, Random

# test Mehrotra predictor-corrector algorithm
m = 15
n = 10
p = 5
x0 = Random.randn(n)
Q = Random.randn(n, n)
Q = Q' * Q
c = Random.randn(n)
G = Random.randn(m, n)
h = G * x0 + ones(m)
A = Random.randn(p, n)
b = A * x0

nθ = 0
p0 = zeros(nθ)
a0 = zeros(m)

# unpack
function unpack(w)
    x = view(w, 1:n)
    s = view(w, n .+ (1:m))
    yx = view(w, n + m .+ (1:p))
    ys = view(w, n + m + p .+ (1:m))
    z = view(w, n + m + p + m .+ (1:m))

    return x, s, yx, ys, z
end

function obj(w)
    x, s, yx, ys, z = unpack(w)
    0.5 * x' * Q * x + c' * x
end

function con(w)
    x, s, yx, ys, z = unpack(w)
    A * x - b
end

# barrier check
function check_variables(w)
    x, s, yx, ys, z = unpack(w)
    if any(s .<= 0.0)
        return true
    end
    if any(z .<= 0.0)
        return true
    end
    return false
end

function r(w, θ, a)
    x, s, yx, ys, z = unpack(w)

    lag_x = Q * x + c + transpose(A) * yx + transpose(G) * ys
    lag_s = ys - z
    c_eq = A * x - b
    c_ineq = s - (h - G * x)

    comp = s .* z
    comp -= a

    return [lag_x; lag_s; c_eq; c_ineq; comp]
end

function Rz(w, θ)
    return ForwardDiff.jacobian(q -> r(q, θ, a0), w)
end

function Rθ(w, θ)
    return ForwardDiff.jacobian(q -> r(w, q, a0), θ)
end

function interior_point(x0, θ;
    x_init_tol = 0.1,
    r_tol = 1.0e-8,
    μ_tol = 1.0e-6,
    max_iter_inner = 100,
    max_iter_ls = 25)

    # initialization
    x = max.(x0, x_init_tol)
    s = ones(m)
    yx = zeros(p)
    ys = zeros(m)
    z = ones(m)

    w = [x; s; yx; ys; z]
    μ = 1.0
    e = ones(m)

    flag = false
    total_iter = 0

    while μ > μ_tol
        for i = 1:max_iter_inner

            # compute residual, residual Jacobian
            res = r(w, θ, μ * e)
            res_norm = norm(res)

            # println("iteration: $i, residual: $res_norm")

            if res_norm < r_tol
                flag = true
                continue
            end

            jac = Rz(w, θ)

            # compute step
            Δ = -jac \ res

            # line search the step direction
            α = 1.0

            iter = 0
            while check_variables(w + α * Δ) # backtrack inequalities
                α = 0.5 * α
                # println("   α = $α")
                iter += 1
                if iter > max_iter_ls
                    @error "backtracking line search fail"
                    flag = false
                    return w, zeros(n + me + mi, p), false, total_iter
                end
            end

            while norm(r(w + α * Δ, θ, μ * e))^2.0 >= (1.0 - 0.001 * α) * res_norm^2.0
                α = 0.5 * α
                # println("   α = $α")
                iter += 1
                if iter > max_iter_ls
                    @error "line search fail"
                    flag = false
                    return w, zeros(n + m + p + m + m, ), false, total_iter
                end
            end

            # update
            w .+= α * Δ

            total_iter += 1
        end

        μ = 0.1 * μ
        println("μ: $μ")
    end

    δw = -1.0 * Rz(w, θ) \ Rθ(w, θ)

    return w, δw, flag, total_iter
end

w_soli, δw_soli, flag, total_iter = interior_point(rand(n), p0)
@show total_iter

function mehrotra(x0, θ;
    x_init_tol = 0.1,
    r_tol = 1.0e-8,
    μ_tol = 1.0e-6,
    max_iter_inner = 100,
    max_iter_ls = 25)

    # initialization
    x = max.(x0, x_init_tol)
    s = ones(m)
    yx = zeros(p)
    ys = zeros(m)
    z = ones(m)
    w = [x; s; yx; ys; z]
    μ = 1.0
    e = ones(m)

    flag = false
    total_iter = 0

    while μ > μ_tol
        for i = 1:max_iter_inner

            # unpack
            x, s, yx, ys, z = unpack(w)

            # compute residual, residual Jacobian
            res_aff = r(w, θ, a0)
            res_aff_norm = norm(res_aff, 1)

            # println("iteration: $i, residual: $res_norm")

            if res_aff_norm < r_tol
                flag = true
                continue
            end

            jac = Rz(w, θ)

            # compute affine step
            Δ_aff = -jac \ res_aff

            Δx_aff, Δs_aff, Δyx_aff, Δys_aff, Δz_aff = unpack(Δ_aff)

            α_aff_pr = min(1.0, minimum([Δsi < 0.0 ? -s[i] / Δsi : Inf for (i, Δsi) in enumerate(Δs_aff)]))
            α_aff_du = min(1.0, minimum([Δzi < 0.0 ? -z[i] / Δzi : Inf for (i, Δzi) in enumerate(Δz_aff)]))
            μ_aff = transpose(s + α_aff_pr * Δs_aff) * (z + α_aff_du * Δz_aff) / m

            # compute corrector step
            σ = (μ_aff / μ)^3.0

            res_corr = r(w, θ, -Δs_aff .* Δz_aff + σ * μ * e)

            Δ_corr = -jac \ res_corr

            # line search the step direction
            η = 1.0
            α_pr = α_aff_pr
            α_du = α_aff_du

            α = η .* [α_pr * ones(n + m + p + m); α_du * ones(m)]

            iter = 0
            while check_variables(w + α .* Δ_corr) # backtrack inequalities
                η = 0.99 * η
                α = η .* [α_pr * ones(n + m + p + m); α_du * ones(m)]

                # println("   α = $α")
                iter += 1
                if iter > max_iter_ls
                    @error "backtracking line search fail"
                    flag = false
                    return w, zeros(n + m + p + m + m, nθ), false, total_iter
                end
            end

            while norm(r(w + α .* Δ_corr, θ, a0))^2.0 >= (1.0 - 0.001 * η) * res_aff_norm^2.0
                η = 0.99 * η
                α = η .* [α_pr * ones(n + m + p + m); α_du * ones(m)]
                # println("   α = $α")
                iter += 1
                if iter > max_iter_ls
                    @error "line search fail"
                    flag = false
                    return w, zeros(n + m + p + m + m, nθ), false, total_iter
                end
            end

            # update
            w .+= α .* Δ_corr

            total_iter += 1
        end

        μ = 0.1 * μ
        println("μ: $μ")
    end

    δw = -1.0 * Rz(w, θ) \ Rθ(w, θ)

    return w, δw, flag, total_iter
end

w_solm, δw_solm, flag, total_iter = mehrotra(rand(n), p0)
@show total_iter

obj(w_soli)
obj(w_solm)

norm(con(w_soli), Inf)
norm(con(w_solm), Inf)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
###############
# Discrete maps
###############
@testset "Discrete example systems" begin

    # Check that initialisation happens correctly.
    @test isa(ar1_unidir(), DiscreteDynamicalSystem)
    @test isa(ar1_bidir(), DiscreteDynamicalSystem)
    @test isa(anishchenko1(), DiscreteDynamicalSystem)
    @test isa(henon2(), DiscreteDynamicalSystem)
    @test isa(ikeda(), DiscreteDynamicalSystem)
    @test isa(nonlinear3d(), DiscreteDynamicalSystem)
    @test isa(ulam(2), DiscreteDynamicalSystem)
    @test isa(logistic2_unidir(), DiscreteDynamicalSystem)
    @test isa(logistic2_bidir(), DiscreteDynamicalSystem)
    @test isa(logistic3(), DiscreteDynamicalSystem)
    @test isa(logistic4(), DiscreteDynamicalSystem)
    @test isa(var1(), DiscreteDynamicalSystem)
    @test isa(verdes(), DiscreteDynamicalSystem)

    # Discrete maps with lags directly return a Dataset
    @test isa(nontrivial_pegiun(n = 10), Dataset)
    @test isa(henon_triple(n = 10), Dataset)

    # Initialise all the systems and generate trajectories
    @test isa(trajectory(ar1_unidir(), 10), Dataset)
    @test isa(trajectory(ar1_bidir(), 10), Dataset)

    @test isa(trajectory(anishchenko1(), 10), Dataset)
    @test isa(trajectory(henon2(), 10), Dataset)
    @test isa(trajectory(ikeda(), 10), Dataset)
    @test isa(trajectory(nonlinear3d(), 10), Dataset)
    @test isa(trajectory(ulam(2), 10), Dataset)
    @test isa(trajectory(logistic2_unidir(), 10), Dataset)
    @test isa(trajectory(logistic2_bidir(), 10), Dataset)

    @test isa(trajectory(logistic3(), 10), Dataset)
    @test isa(trajectory(logistic4(), 10), Dataset)
    @test isa(trajectory(var1(), 10), Dataset)
    @test isa(trajectory(verdes(), 10), Dataset)

end


text:
function f()
    a::Float64 = 0
    s = 1

    println("start")
    for i in 1:10000000000
        a += 1 / (2 * i - 1) * s
        s *= -1
    end

    println(a * 4)
end

f()