text:
using Test, NiLang

@testset "complex" begin
    a = 1.0+ 2im
    @instr (a |> real) += 2
    @instr (a |> imag) += 2
    @test a == 3.0 + 4im

    a = 1.0+ 2im
    @instr a += complex(2.0, 2.0)
    @test a == 3.0 + 4.0im
    @i function f(loss, a::Complex{T}, b) where T
        @routine begin
            c ← zero(a)
            sq ← zero(T)
            c += a * b
            sq += (c |> real) ^ 2
            sq += (c |> imag) ^ 2
        end
        loss += sq ^ 0.5
        ~@routine
    end
    a = 1.0 + 2.0im
    b = 2.0 + 1.0im
    loss = 0.0
    @instr f(loss, a, b)
    @test loss ≈ abs(a*b)
end

@testset "complex arithmetics" begin
    for op in [exp, log, identity]
        x, y = 2.0+1.0im, 0.5+0.2im
        @instr x += op(y)
        @test x ≈ 2.0+1.0im + op(0.5+0.2im)
    end
    for op in [SWAP, HADAMARD]
        x, y = 2.0+1.0im, 0.5+0.2im
        @instr op(x, y)
        @test x ≈ op(2.0+1.0im, 0.5+0.2im)[1]
        @test y ≈ op(2.0+1.0im, 0.5+0.2im)[2]
    end
    for op in [NEG, INC, DEC]
        x = 2.0+1.0im
        @instr op(x)
        @test x ≈ op(2.0+1.0im)
    end
    for op in [^, /, +, -]
        x, y, z = 2.0+1.0im, 0.5+0.2im, 0.8-2.0im
        @instr PlusEq(op)(x, y, z)
        @test x ≈ 2.0+1.0im + op(0.5+0.2im, 0.8-2.0im)
    end
end


text:
using ZMQ
include("zmqext.jl")

REQUEST_TIMEOUT = 2500
REQUEST_RETRIES = 3
SERVER_ENDPOINT = "tcp://127.0.0.1:5555"

context = ZMQ.Context()
client2 = ZMQ.Socket(context, REQ)

function startServer(context::ZMQ.Context)
    @info "Connecting to server..."
    client = ZMQ.Socket(context, ZMQ.REQ)
    ZMQ.connect(client, SERVER_ENDPOINT)

    counter = 1
    while true
        request = counter
        @info "Sending ($(request))"
        ZMQ.send(client, request)

        retries_left = REQUEST_RETRIES
        while true
            items, err = ZMQExt.poll([client], ZMQ.POLLIN, REQUEST_TIMEOUT)
            if err == -1
                println("interrupted, break")
                break
            end

            if items[1] & ZMQ.POLLIN != 0
                reply = ZMQ.recv(client, Int)
                @info "got reply ($(reply)) - type: $(typeof(reply))"
                @info "Server replied OK ($(reply))"
                retries_left = REQUEST_RETRIES
                break
            end

            retries_left -= 1
            @warn "No response from server"
            client.linger = 0
            ZMQ.close(client)
            if retries_left == 0
                @error "Server seems to be offline, abandoning"
                exit(1)
            end

            @info "Reconnecting to server"
            client = ZMQ.Socket(context, ZMQ.REQ)
            ZMQ.connect(client, SERVER_ENDPOINT)
            @info "Resending ($(request))"
            ZMQ.send(client, request)
        end

        counter += 1
    end
end

function test()
    println("testing..")

end

if abspath(PROGRAM_FILE) == @__FILE__
    startServer(context)
end

ext: jl

max_stars_repo_name: jwdevantier/julia-zmqtest

max_issues_repo_name: jwdevantier/julia-zmqtest

max_forks_repo_name: jwdevantier/julia-zmqtest


text:
"""
```
MaterialParameters(;
    crystalStruct::AbstractCrystalStruct,
    μ = 1.0,
    μMag = 1.0,
    ν = 0.2,
    σPN = 0.0,
)
```
Creates [`MaterialParameters`](@ref).
"""
function MaterialParameters(;
    crystalStruct::AbstractCrystalStruct,
    μ = 1.0,
    μMag = 1.0,
    ν = 0.2,
    σPN = 0.0,
)
    omνInv = 1 / (1 - ν)
    opνInv = 1 / (1 + ν)
    νomνInv = ν * omνInv        # ν / (1 - ν)
    νopνInv = ν * opνInv        # ν / (1 + ν)
    μ4π = μ / (4π)              # μ / (4π)
    μ8π = μ4π / 2               # μ / (8π)
    μ4πν = μ4π * omνInv         # μ / (1 - ν)
    νμ4πν = ν * μ4πν            # νμ / (1 - ν)
    omνInv8π = omνInv / (8π)    # 1 / (8π (1 - ν))
    om2νomνInv8π = (1 - 2ν) * omνInv8π # (1 - 2ν) / (8π (1 - ν))

    return MaterialParameters(
        crystalStruct,
        μ,
        μMag,
        ν,
        omνInv,
        opνInv,
        νomνInv,
        νopνInv,
        μ4π,
        μ8π,
        μ4πν,
        νμ4πν,
        omνInv8π,
        om2νomνInv8π,
        σPN,
    )
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
# code by IainNZ
function histogram(x, bins::Int=sturges(length(x)); args...)
  edges, counts = hist(x, bins)
  labels = (@compat AbstractString)[]
  binwidth = edges.step / edges.divisor
  for i in 1:length(counts)
		val = floatRoundLog10(edges[i])
    push!(labels, string("(", val, ",", floatRoundLog10(val+binwidth), "]"))
  end
  barplot(labels, counts; symb="▇", args...)
end

function histogram(x; bins::Int=sturges(length(x)), args...)
  histogram(x, bins; args...)
end

# Sturges' rule
sturges(n) = n == 0 ? 1 : ceil(Int, log2(n) + 1)


text:
"""
Object with a *term* and an applicable *BasisFunction* and a eventfield that are later passed to the basisfunction.

$(TYPEDEF)
$(TYPEDSIGNATURES)
$(FIELDS)
# Examples
```julia-repl
julia>  b = TimeExpandedTerm(term,kernel,[:latencyTR,:durationTR])
```
"""
struct TimeExpandedTerm{T<:AbstractTerm} <: AbstractTerm
    "Term that the basis function is applied to. This is regularly called in other functions to get e.g. term-coefnames and timeexpand those"
    term::T
    "Kernel that determines what should happen to the designmatrix of the term"
    basisfunction::BasisFunction
    "Which fields of the event-table should be passed to the basisfunction.Important: The first entry has to be the event-latency in samples!"
    eventfields::Array{Symbol}
end

function TimeExpandedTerm(
    term::NTuple{N,Union{<:AbstractTerm,<:MixedModels.RandomEffectsTerm}},
    basisfunction,
    eventfields::Array{Symbol,1},
) where {N}
    # for mixed models, apply it to each Term
    TimeExpandedTerm.(term, Ref(basisfunction), Ref(eventfields))
end
function TimeExpandedTerm(term, basisfunction, eventfields::Nothing)
    # if the eventfield is nothing, call default
    TimeExpandedTerm(term, basisfunction)
end

function TimeExpandedTerm(term, basisfunction, eventfields::Symbol)
    # call with symbol, convert to array of symbol
    TimeExpandedTerm(term, basisfunction, [eventfields])
end

function TimeExpandedTerm(term, basisfunction; eventfields = [:latency])
    # default call, use latency
    TimeExpandedTerm(term, basisfunction, eventfields)
end

collabel(term::TimeExpandedTerm) = collabel(term.basisfunction)
StatsModels.width(term::TimeExpandedTerm) = width(term.basisfunction)
StatsModels.terms(t::TimeExpandedTerm) = terms(t.term)

function Base.show(io::IO, p::TimeExpandedTerm)
    #print(io, "timeexpand($(p.term), $(p.basisfunction.type),$(p.basisfunction.times))")
    println(io, "$(coefnames(p))")
end

"""
$(SIGNATURES)
Combine two UnfoldDesignmatrices. This allows combination of multiple events.

This also allows to define events with different lengths.

Not supported for models without timebasis, as it is not needed there (one can simply run multiple models)
# Examples
```julia-repl
julia>  basisfunction1 = firbasis(τ=(0,1),sfreq = 10,name="basis1")
julia>  basisfunction2 = firbasis(τ=(0,0.5),sfreq = 10,name="basis2")
julia>  Xdc1          = designmatrix(UnfoldLinearModelContinuousTime(Dict(Any=>(@formula 0~1,basisfunction1)),tbl_1)
julia>  Xdc2          = designmatrix(UnfoldLinearModelContinuousTime(Dict(Any=>(@formula 0~1,basisfunction2)),tbl_2)
julia>  Xdc = Xdc1+Xdc2
```

"""
function combineDesignmatrices(X1::DesignMatrix, X2::DesignMatrix)

    X1 = deepcopy(X1)
    X2 = deepcopy(X2)
    Xs1 = get_Xs(X1.Xs)
    Xs2 = get_Xs(X2.Xs)

    sX1 = size(Xs1, 1)
    sX2 = size(Xs2, 1)

    # append 0 to the shorter designmat
    if sX1 < sX2
        Xs1 = SparseMatrixCSC(sX2, Xs1.n, Xs1.colptr, Xs1.rowval, Xs1.nzval)
    elseif sX2 < sX1
        Xs2 = SparseMatrixCSC(sX1, Xs2.n, Xs2.colptr, Xs2.rowval, Xs2.nzval)
    end

    Xcomb = hcat(Xs1, Xs2)

    #if !(Float64(X1.formulas.rhs.basisfunction.times.step) ≈ Float64(X2.formulas.rhs.basisfunction.times.step))
    #        @warn("Concatenating formulas with different sampling rates. Be sure that this is what you want.")
    #end
    if typeof(X1.Xs) <: Tuple
        # we have random effects
        # combine REMats in single-eventtpe formulas ala y ~ (1|x) + (a|x)
        Xs1 = MixedModels._amalgamate([X1.Xs[2:end]...], Float64)
        Xs2 = MixedModels._amalgamate([X2.Xs[2:end]...], Float64)

        Xcomb = (Xcomb, Xs1..., Xs2...)

        # Next we make the ranefs all equal size
        equalizeReMatLengths!(Xcomb[2:end])

        # check if ranefs can be amalgamated. If this fails, then MixedModels tried to amalgamate over different eventtypes and we should throw the warning
        # if it success, we have to check if the size before and after is identical. If it is not, it tried to amalgamize over different eventtypes which were of the same length

        try
            reterms = MixedModels._amalgamate([Xcomb[2:end]...], Float64)
            if length(reterms) != length(Xcomb[2:end])
                throw("IncompatibleRandomGroupings")
            end
        catch e
            @error "Error, you seem to have two different eventtypes with the same random-effect grouping variable. \n
            This is not allowed, you have to rename one. Example:\n
            eventA: y~1+(1|item) \n
            eventB: y~1+(1|item)  \n
            This leads to this error. Rename the later one\n
            eventB: y~1+(1|itemB) "
            throw("IncompatibleRandomGroupings")
        end

    end

    if typeof(X1.formulas) <: FormulaTerm
        DesignMatrix([X1.formulas X2.formulas], Xcomb, [X1.events, X2.events])
    else

        DesignMatrix([X1.formulas... X2.formulas], Xcomb, [X1.events..., X2.events])
    end
end

function changeMatSize!(m, fe::AbstractSparseMatrix, remats)
    changeReMatSize!.(remats, Ref(m))
    fe = SparseMatrixCSC(m, fe.n, fe.colptr, fe.rowval, fe.nzval)
    return (fe, remats...)
end
function changeMatSize!(m, fe::AbstractMatrix, remats)
    changeReMatSize!.(remats, Ref(m))
    fe = fe[1:m,:]
    return (fe, remats...)
end
function changeReMatSize!(remat::MixedModels.AbstractReMat, m::Integer)

    n = m - length(remat.refs) # missing elements
    if n < 0
        deleteat!(remat.refs, range(m + 1, stop = length(remat.refs)))
        remat.adjA = remat.adjA[:, 1:(m+1)]
        remat.wtz = remat.wtz[:, 1:(m+1)]
        remat.z = remat.z[:, 1:(m+1)]
    elseif n > 0
        append!(remat.refs, repeat([remat.refs[end]], n))
        remat.adjA = [remat.adjA sparse(repeat([0], size(remat.adjA)[1], n))]
        remat.wtz = [remat.wtz zeros((size(remat.wtz)[1], n))]
        remat.z = [remat.z zeros((size(remat.z)[1], n))]
    end

    #ReMat{T,S}(remat.trm, refs, levels, cnames, z, wtz, λ, inds, adjA, scratch)
end
function equalizeReMatLengths!(remats::NTuple{A,MixedModels.AbstractReMat}) where {A}
    # find max length
    m = maximum([x[1] for x in size.(remats)])
    @debug print("combining lengths: $m")
    # for each reMat
    for k in range(1, length = length(remats))
        remat = remats[k]
        if size(remat)[1] == m
            continue
        end
        # prolong if necessary
        changeReMatSize!(remat, m)

    end
end

Base.:+(X1::DesignMatrix, X2::DesignMatrix) = combineDesignmatrices(X1, X2)

Base.:+(X1::Nothing, X2::DesignMatrix) = X2

function get_Xs(Xs::Tuple)
    return Xs[1]
end
function get_Xs(Xs::SparseMatrixCSC)
    return Xs
end
function get_Xs(Xs::Matrix)
    # mass univariate case
    return Xs
end
function get_Xs(X::DesignMatrix)
    return get_Xs(X.Xs)
end

"""
$(SIGNATURES)
designmatrix(type, f, tbl; kwargs...)
Return a *DesignMatrix* used to fit the models.
# Arguments
- type::UnfoldModel
- f::FormulaTerm: Formula to be used in this designmatrix
- tbl: Events (usually a data frame) to be modelled
- basisfunction::BasisFunction: basisfunction to be used in modeling (if specified)
- contrasts::Dict: (optional) contrast to be applied to formula
- eventfields::Array: (optional) Array of symbols which are passed to basisfunction event-wise.
First field of array always defines eventonset in samples. Default is [:latency]

# Examples
```julia-repl
julia>  designmatrix(UnfoldLinearModelContinuousTime,Dict(Any=>(f,basisfunction1),tbl)
```

"""
function designmatrix(
    type,
    f::Union{Tuple,FormulaTerm},
    tbl,
    basisfunction;
    contrasts = Dict{Symbol,Any}(),
    kwargs...,
)
    @debug("generating DesignMatrix")
    form = apply_schema(f, schema(f, tbl, contrasts), MixedModels.LinearMixedModel)
    form =
        apply_basisfunction(form, basisfunction, get(Dict(kwargs), :eventfields, nothing))

    # Evaluate the designmatrix
    if (!isnothing(basisfunction)) & (type <: UnfoldLinearMixedModel)
        X = modelcols.(form.rhs, Ref(tbl))
    else
        X = modelcols(form.rhs, tbl)
    end
    return DesignMatrix(form, X, tbl)
end
function designmatrix(type, f, tbl; kwargs...)
    return designmatrix(type, f, tbl, nothing; kwargs...)
end

# specify for abstract interface
designmatrix(uf::UnfoldModel) = uf.designmatrix

function designmatrix(
    uf::UnfoldModel,
    tbl;
    eventcolumn = :event,
    contrasts = Dict{Symbol,Any}(),
    kwargs...,
)
    X = nothing
    fDict = design(uf)
    for (eventname, f) in pairs(fDict)
        if eventname == Any
            eventTbl = tbl
        else
            if !(eventcolumn ∈ names(tbl))
                error(
                    "Couldnt find columnName: " *
                    string(eventcolumn) *
                    " in event-table.  Maybe need to specify eventcolumn=:correctColumnName (default is ':event') \n names(tbl) = " *
                    join(names(tbl), ","),
                )
            end
            eventTbl = tbl[tbl[:, eventcolumn].==eventname, :]
        end
        if isempty(eventTbl)
            error(
                "eventTable empty after subsetting. Couldnt find event '" *
                string(eventname) *
                "'{" *
                string(typeof(eventname)) *
                "}, in field tbl[:,:" *
                string(eventcolumn) *
                "].? - maybe you need to specify it as a string instead of a symbol?",
            )
        end

        fIx = collect(typeof.(f) .<: FormulaTerm)
        bIx = collect(typeof.(f) .<: BasisFunction)

        if any(bIx)
            # timeContinuos way
            # TODO there should be a julian way to do this distinction

            X =
                X + designmatrix(
                    typeof(uf),
                    f[fIx],
                    eventTbl,
                    collect(f[bIx])[1];
                    contrasts = contrasts,
                    kwargs...,
                )
        else
            # normal way
            X =
                X +
                designmatrix(typeof(uf), f[fIx], eventTbl; contrasts = contrasts, kwargs...)
        end
    end
    return X
end

import Base.isempty
Base.isempty(d::DesignMatrix) = isempty(d.Xs)

"""
$(SIGNATURES)
timeexpand the rhs-term of the formula with the basisfunction

"""
function apply_basisfunction(form, basisfunction::BasisFunction, eventfields)
    @debug("apply_basisfunction")
    return FormulaTerm(form.lhs, TimeExpandedTerm(form.rhs, basisfunction, eventfields))
end

function apply_basisfunction(form, basisfunction::Nothing, eventfields)
    # in case of no basisfunctin, do nothing
    return form
end

function designmatrix!(uf::UnfoldModel, evts; kwargs...)
    X = designmatrix(uf, evts; kwargs...)
    uf.designmatrix = X
end

function StatsModels.modelmatrix(uf::UnfoldLinearModel,basisfunction)
    if basisfunction
        @warn("basisfunction not defined for this kind of model")
    else
        return modelmatrix(uf)
    end
end

function StatsModels.modelmatrix(uf::UnfoldLinearModelContinuousTime,basisfunction = true)
    if basisfunction
        return modelmatrix(designmatrix(uf))
    else
        # replace basisfunction with non-timeexpanded one
        f = formula(uf)

        # probably a more julian way to do this...
        if isa(f,AbstractArray)
            return modelcols_nobasis.(f[1,:],events(uf))
        else
            return modelcols_nobasis(f,events(uf))
        end

    end
end

modelcols_nobasis(f::FormulaTerm,tbl::DataFrame) = modelcols(f.rhs.term,tbl)
StatsModels.modelmatrix(uf::UnfoldModel) = modelmatrix(designmatrix(uf))#modelmatrix(uf.design,uf.designmatrix.events)
StatsModels.modelmatrix(d::DesignMatrix) = d.Xs

StatsModels.modelmatrix(d::Dict, events) = modelcols(formula(d).rhs, events)

formula(uf::UnfoldModel) = formula(designmatrix(uf))
formula(d::DesignMatrix) = d.formulas

events(uf::UnfoldModel) = events(designmatrix(uf))
events(d::DesignMatrix) = d.events

design(uf::UnfoldModel) = uf.design

function formula(d::Dict) #TODO Specify Dict better
    if length(values(d)) == 1

        return [c[1] for c in collect(values(d))][1]
    else
    return [c[1] for c in collect(values(d))]
    end
    #return collect(values(d))[1][1] # give back first formula for now
end

"""
$(SIGNATURES)
calculates the actual designmatrix for a timeexpandedterm. Multiple dispatch on StatsModels.modelcols
"""
function StatsModels.modelcols(term::TimeExpandedTerm, tbl)

    X = modelcols(term.term, tbl)

    time_expand(X, term, tbl)
end

mutable struct SparseReMat{T,S} <: MixedModels.AbstractReMat{T}
    trm::Any
    refs::Vector{Int32}
    levels::Any
    cnames::Vector{String}
    z::SparseMatrixCSC{T}
    wtz::SparseMatrixCSC{T}
    λ::LowerTriangular{T,Matrix{T}}
    inds::Vector{Int}
    adjA::SparseMatrixCSC{T,Int32}
    scratch::Matrix{T}
end

"""
$(SIGNATURES)
This function timeexpands the random effects and generates a ReMat object
"""
function StatsModels.modelcols(
    term::TimeExpandedTerm{<:Union{<:RandomEffectsTerm,<:MixedModels.AbstractReTerm}},
    tbl,
)
    # exchange this to get ZeroCorr to work

    tbl = DataFrame(tbl)
    # get the non-timeexpanded reMat
    reMat = modelcols(term.term, tbl)

    # Timeexpand the designmatrix
    z = transpose(time_expand(transpose(reMat.z), term, tbl))

    # First we check if there is overlap in the timeexpanded term. If so, we cannot continue. Later implementations will remedy that
    #println(dump(term,))
    if hasfield(typeof(term.term), :rhs)
        rhs = term.term.rhs

    elseif hasfield(typeof(term.term.term), :rhs)
        # we probably have something like zerocorr, which does not need to show a .rhs necessarily
        rhs = term.term.term.rhs
    else

        println("term.term: $(dump(term.term))")
        error("unknown RE structure, has no field .rhs:$(typeof(term.term))")

    end

    group = tbl[!, rhs.sym]
    time = tbl[!, term.eventfields[1]]

    # get the from-to onsets of the grouping varibales
    onsets = time_expand_getRandomGrouping(group, time, term.basisfunction)
    #print(size(reMat.z))
    refs = zeros(size(z)[2]) .+ 1
    for (i, o) in enumerate(onsets[2:end])
        # check for overlap
        if (minimum(o) <= maximum(onsets[i+1])) & (maximum(o) <= minimum(onsets[i+1]))
            error("overlap in random effects structure detected, not currently supported")
        end
    end

    # From now we can assume no overlap
    # We want to fnd which subject is active when
    refs = zeros(size(z)[2]) .+ 1
    uGroup = unique(group)

    for (i, g) in enumerate(uGroup[1:end])

        ix_start = findfirst(g .== group)
        ix_end = findlast(g .== group)
        if i == 1
            time_start = 1
        else
            time_start = time[ix_start]
            # XXX Replace this functionality by shiftOnset?
            time_start = time_start - sum(times(term.basisfunction) .<= 0)
        end
        if i == length(uGroup)
            time_stop = size(refs, 1)
        else
            time_stop = time[ix_end]
            time_stop = time_stop + sum(times(term.basisfunction) .> 0)
        end
        if time_start < 0
            time_start = 1
        end

        if time_stop > size(refs, 1)
            time_stop = size(refs, 1)
        end

        #println("$g,$time_start,$time_stop")
        refs[Int64(time_start):Int64(time_stop)] .= i
    end

    # Other variables with implementaions taken from the LinerMixedModel function
    wtz = z
    trm = term

    S = size(z, 1)
    T = eltype(z)
    λ = LowerTriangular(Matrix{T}(I, S, S))

    inds = MixedModels.sizehint!(Int[], (S * (S + 1)) >> 1) # double trouble? can this line go?
    m = reshape(1:abs2(S), (S, S))
    inds = sizehint!(Int[], (S * (S + 1)) >> 1)
    for j = 1:S
        for i = j:S
            # We currently restrict to diagonal entries
            # Once mixedmodels#293 is pushed, we can relax this and use zerocorr()
            if !(typeof(term.term) <: MixedModels.ZeroCorr) || (i == j) # for diagonal
                push!(inds, m[i, j])
            end
        end
    end

    levels = reMat.levels
    refs = refs

    # reMat.levels doesnt change
    cnames = coefnames(term)
    #print(refs)
    adjA = MixedModels.adjA(refs, z)
    scratch = Matrix{T}(undef, (S, length(uGroup)))

    # Once MixedModels.jl supports it, can be replaced with: SparseReMat
    ReMat{T,S}(rhs, refs, levels, cnames, z, wtz, λ, inds, adjA, scratch)
end

"""
$(SIGNATURES)
Get the timeranges where the random grouping variable was applied
"""
function time_expand_getRandomGrouping(tblGroup, tblLatencies, basisfunction)
    ranges = time_expand_getTimeRange.(tblLatencies, Ref(basisfunction))
end

# helper function to get the ranges from where to where the basisfunction is added
function time_expand_getTimeRange(onset, basisfunction)
    npos = sum(times(basisfunction) .>= 0)
    nneg = sum(times(basisfunction) .< 0)

    #basis = kernel(basisfunction)(onset)

    fromRowIx = floor(onset) - nneg
    toRowIx = floor(onset) + npos

    range(fromRowIx, stop = toRowIx)
end

"""
$(SIGNATURES)
performs the actual time-expansion in a sparse way.

 - Get the non-timeexpanded designmatrix X from StatsModels.
 - evaluate the basisfunction kernel at each event
 - calculate the necessary rows, cols and values for the sparse matrix
 Returns SparseMatrixCSC
"""
function time_expand(X, term, tbl)
    ncolsBasis = size(kernel(term.basisfunction)(0), 2)
    X = reshape(X, size(X, 1), :)

    ncolsX = size(X)[2]
    nrowsX = size(X)[1]
    ncolsXdc = ncolsBasis * ncolsX

    # this is the predefined eventfield, usually "latency"
    #println(term.eventfields)
    #println(tbl[term.eventfields[1]][1:10])
    tbl = DataFrame(tbl)
    onsets = tbl[!, term.eventfields[1]]

    if typeof(term.eventfields) <: Array && length(term.eventfields) == 1
        bases = kernel(term.basisfunction).(tbl[!, term.eventfields[1]])
    else
        bases = kernel(term.basisfunction).(eachrow(tbl[!, term.eventfields]))
    end

    # generate rowindices
    rows = copy(rowvals.(bases))
    # this shift is necessary as some basisfunction time-points can be negative. But a matrix is always from 1:τ. Thus we have to shift it backwards in time.
    # The onsets are onsets-1 XXX not sure why.
    for r = 1:length(rows)
        rows[r] .+= floor(onsets[r] - 1) + shiftOnset(term.basisfunction)
    end

    rows = vcat(rows...)
    rows = repeat(rows, ncolsX)

    # generate column indices
    cols = []
    for Xcol = 1:ncolsX
        for b = 1:length(bases)
            for c = 1:ncolsBasis
                push!(
                    cols,
                    repeat([c + (Xcol - 1) * ncolsBasis], length(nzrange(bases[b], c))),
                )
            end
        end
    end
    cols = vcat(cols...)

    # generate values
    vals = []
    for Xcol = 1:ncolsX
        push!(vals, vcat(nonzeros.(bases) .* X[:, Xcol]...))
    end

    vals = vcat(vals...)
    ix = rows .> 0
    A = sparse(rows[ix], cols[ix], vals[ix])

    return A
end

"""
$(SIGNATURES)
coefnames of a TimeExpandedTerm concatenates the basis-function name with the kronecker product of the term name and the basis-function colnames. Separator is ' : '
Some examples for a firbasis:
        basis_313 : (Intercept) : 0.1
        basis_313 : (Intercept) : 0.2
        basis_313 : (Intercept) : 0.3
        ...
"""
function StatsModels.coefnames(term::TimeExpandedTerm)
    terms = coefnames(term.term)
    colnames = Unfold.colnames(term.basisfunction)
    name = Unfold.name(term.basisfunction)
    if typeof(terms) == String
        terms = [terms]
    end
    return name .* " : " .* kron(terms .* " : ", string.(colnames))
end

function termnames(term::TimeExpandedTerm)
    terms = coefnames(term.term)
    colnames = colnames(term.basisfunction)
    if typeof(terms) == String
        terms = [terms]
    end
    return vcat(repeat.([[t] for t in terms], length(colnames))...)
end

function colname_basis(term::TimeExpandedTerm)
    terms = coefnames(term.term)
    colnames = colnames(term.basisfunction)
    if typeof(terms) == String
        terms = [terms]
    end
    return repeat(colnames, length(terms))
end

function StatsModels.coefnames(terms::AbstractArray{<:FormulaTerm})
    return coefnames.(Base.getproperty.(terms, :rhs))
end
function StatsModels.coefnames(term::MixedModels.ZeroCorr)
    coefnames(term.term)
end

function StatsModels.coefnames(term::RandomEffectsTerm)
    coefnames(term.lhs)
end

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
"""
    function CombineStencils(varstencil::StencilDefinition...)

Adds 2 or more stencils together, combines their coefficients.
Result has the radius of the bigger stencil.
Can only add stencils, substruction should be embedded in the stencil.

Returns:
    symbolic expression (Basic)
    max_radius          (Int)
"""
function CombineStencils(varstencil::StencilDefinition...)

    total_max_radius = -1
    for s in varstencil
        if s.max_radius > total_max_radius
            total_max_radius = s.max_radius
        end

    end

    Dsym = [symbols("w$(i)_$(j)_$(k)") for i = 1:(2total_max_radius+1),
                                        j = 1:(2total_max_radius+1),
                                        k = 1:(2total_max_radius+1)]
    global D = Dsym
    tot_s = 0
    for s in varstencil
        resss = convert_to_sym(s, total_max_radius)
        tot_s += resss

    end
    return expand(tot_s),total_max_radius
end

"""
function convert_to_sym(st::StencilDefinition, radius::Integer)

Converts a stencil expression, either array notation or coeff matrix, to
symbolic math. The extra radius parameter can be bigger than the defined
stencil's radius.
"""
function convert_to_sym(st::StencilDefinition, radius::Integer)
    if st.expr isa Nothing
        ## For coef stencil
        coefs = copy(st.coefs)
        if st.max_radius < radius
            pad = radius - st.max_radius
            coefs = zeros(Float32, 2*radius+1, 2*radius+1, 2*radius+1)

            for i = 1:(2st.max_radius+1), j = 1:(2st.max_radius+1), k = 1:(2st.max_radius+1)
                coefs[pad+i, pad+j,pad+k] = st.coefs[i,j,k]
            end
        end
        #println(size(coefs))
        global D
        st_sym = sum(coefs.*D)
        if st.uses_vsq
            global v
            st_sym *= v
        end
        return st_sym
    else
        ## For expression stencil
        coefs = nothing #copy(st.coefs[:])
        expr = copy(st.expr)
        for i in [:x :y :z]
            expr = replace_symbol(i, 0, expr)
        end
        if !(st.coefs isa Nothing)
            coefs = copy(st.coefs[:])
            if length(coefs) < radius
                pad = radius - length(coefs)
                for i = 1:pad
                    push!(coefs, 0)
                end
            end
            global c = [reverse(coefs[2:end]); coefs[:]]
        end

        expr = fix_indices(expr, radius+1)
        #println(expr)
        st_sym = eval(expr)
        return st_sym
    end
end

# function CheckTimestepCombinationLimit(limit_radius, radius, m_steps)
#     cr = radius
#     cm = 2
#     for i = 2:m_steps
#         cr = 2*cr
#         if cr > 16
#             cm = i-1
#             break
#         end
#     end
#
#     if cm < 2
#         return false
#     else
#         return cm
#     end
#
# end
function CreateSymOffMatrices(radius)
    Dsym = [symbols("w$(i)_$(j)_$(k)") for i = 1:(2radius+1),
                                        j = 1:(2radius+1),
                                        k = 1:(2radius+1)]

    Doff = [[i,j,k] for i = -radius:radius, j=-radius:radius, k=-radius:radius]

    return Dsym, Doff
end
function CreateOffset(radius)
    Doff = [[i,j,k] for i = -radius:radius, j=-radius:radius, k=-radius:radius]

    return Doff
end

function CombineTimeSteps(base_expr::Basic, base_radius::Int,
                            expr::Basic, radius::Int)

    b_sym, b_off = CreateSymOffMatrices(base_radius)
    sym,off = CreateSymOffMatrices(radius)
    new_sym, new_off = CreateSymOffMatrices(radius+base_radius)

    new_eq = 0

    for i = 1:(2*radius+1), j = 1:(2*radius+1),k = 1:(2*radius+1)
        c1 = coeff(expr, sym[i,j,k])
        (c1 == 0) && continue
        for ii = 1:(2*base_radius+1), jj = 1:(2*base_radius+1),kk = 1:(2*base_radius+1)
            c2 = coeff(base_expr, b_sym[ii,jj,kk])
            (c2 == 0) && continue
            off2 = b_off[ii,jj,kk] + [i;j;k] .+ base_radius
            new_eq += c1*c2*new_sym[off2...]
        end
    end
    return new_eq
end

function combine_more_time_steps(base_expr, base_radius, steps)
    steps = steps - 1
    final_radius = (1+steps)*base_radius
    all_syms,temp = CreateSymOffMatrices(final_radius)
    new_eq = 0
    b_sym = @view all_syms[1:2base_radius+1,1:2base_radius+1,1:2base_radius+1]
    b_off = CreateOffset(base_radius)
    expr = base_expr
    for t = 1:steps
        new_eq = 0
        radius = (t)*base_radius
        next_rad = (t+1)*base_radius
        off = CreateOffset(radius)
        sym = @view all_syms[1:2radius+1,1:2radius+1,1:2radius+1]
        new_off = CreateOffset(next_rad)
        new_sym = @view all_syms[1:2next_rad+1,1:2next_rad+1,1:2next_rad+1]
        for i = 1:(2*radius+1), j = 1:(2*radius+1),k = 1:(2*radius+1)
            c1 = coeff(expr, sym[i,j,k])
            (c1 == 0) && continue
            for ii = 1:(2*base_radius+1), jj = 1:(2*base_radius+1),kk = 1:(2*base_radius+1)
                c2 = coeff(base_expr, b_sym[ii,jj,kk])
                (c2 == 0) && continue
                off2 = b_off[ii,jj,kk] + [i;j;k] .+ base_radius
                new_eq += c1*c2*new_sym[off2...]
            end
        end
        expr = new_eq
    end
    return expr

end

function matrix_from_expression(expr, radius)
    sym,temp = CreateSymOffMatrices(radius)
    template = zeros(Float32, 2radius+1,2radius+1,2radius+1)
    for i = 1:(2*radius+1), j = 1:(2*radius+1),k = 1:(2*radius+1)
        c1 = coeff(expr, sym[i,j,k])
        (c1 == 0) && continue
        template[i,j,k] = c1
    end
    return template
end
"""
    function CombineTimeSteps(expr, radius)

Transform a stencil so it applies an additional timestep each execution.
"""
# function CombineTimeSteps(expr::Basic, radius::Integer)
#
#     Dsym = [symbols("w$(i)_$(j)_$(k)") for i = 1:(2radius+1),
#                                         j = 1:(2radius+1),
#                                         k = 1:(2radius+1)]
#
#     Doff = [[i,j,k] for i = -radius:radius, j=-radius:radius, k=-radius:radius]
#     new_Dsym = [symbols("w$(i)_$(j)_$(k)") for i = 1:(4radius+1),
#                                         j = 1:(4radius+1),
#                                         k = 1:(4radius+1)]
#     new_Doff =  [[i,j,k] for i = -2radius:2radius, j=-2radius:2radius, k=-2radius:2radius]
#
#     new_eq = 0
#     #@show 2*radius + 1
#     for i = 1:(2*radius+1), j = 1:(2*radius+1),k = 1:(2*radius+1)
#         #global d, Doff, Dsym, new_Dsym, new_Doff,new_eq,radius
#         #@show (i,j,k)
#         if coeff(expr,Dsym[i,j,k]) != 0
#             Dsym[i,j,k]
#             c1 = coeff(expr,Dsym[i,j,k])
#             for ii = 1:(2*radius+1), jj = 1:(2*radius+1),kk = 1:(2*radius+1)
#                 if coeff(expr,Dsym[ii,jj,kk]) != 0
#                     c = coeff(expr,Dsym[ii,jj,kk])
#                     off2 = Doff[ii,jj,kk] + [i;j;k] .+ radius
#                     new_eq += c1*c*new_Dsym[off2...]
#                 end
#             end
#         end
#
#     end
#     return new_eq
# end


text:
function finalrelease!(node::Int64, error, params::Parameters)
  params.del[4, 1, node] += params.si[4, 2, node] * params.f[2, 1, node]
  params.del[5, 1, node] += params.si[5, 1, node] * params.f[1, 1, node]
  params.del[1, 1, node] += params.si[1, 1, node] * params.f[1, 1, node]
  params.del[2, 1, node] += params.si[2, 2, node] * params.f[2, 1, node]

  if abs(params.nodes[node].fc) > eps() &&
      (abs(params.bf[1, 1, node]) > eps() || abs(params.bf[2, 1, node]) > eps())
    params.friction[node] =
      params.nodes[node].fc * sqrt(params.bf[1, 1, node]^2 + params.bf[2, 1, node]^2)
  end

  for j in 1:6
    error += abs(params.del[j, 1, node] - params.rdel[j, 1])
  end
  error
end


text:
"Augmented Lagrangian solve"
function solve!(prob::Problem{T,Discrete}, solver::AugmentedLagrangianSolver{T}) where T<:AbstractFloat
    reset!(solver)
    t_start = time()

    solver_uncon = solver.solver_uncon

    prob_al = AugmentedLagrangianProblem(prob, solver)
    logger = default_logger(solver)

    rollout!(prob)

    with_logger(logger) do
        record_iteration!(prob_al, solver, cost(prob_al), solver_uncon)
        println(logger,OuterLoop)
        for i = 1:solver.opts.iterations
            set_tolerances!(solver,solver_uncon,i)
            J = step!(prob_al, solver, solver_uncon)

            record_iteration!(prob, solver, J, solver_uncon)

            converged = evaluate_convergence(solver)
            println(logger,OuterLoop)
            converged ? break : nothing

            reset!(solver_uncon)
        end
    end
    solver.stats[:time] = time() - t_start
    return solver
end

function solve!(prob::Problem{T,Discrete},opts::AugmentedLagrangianSolverOptions{T}) where T<:AbstractFloat
    !is_constrained(prob) ? solver = AbstractSolver(prob,opts.opts_uncon) : solver = AbstractSolver(prob,opts)
    solve!(prob,solver)
end

"Set intermediate convergence tolerances for unconstrained solve"
function set_tolerances!(solver::AugmentedLagrangianSolver{T},
        solver_uncon::AbstractSolver{T},i::Int) where T
    if i != solver.opts.iterations
        solver_uncon.opts.cost_tolerance = solver.opts.cost_tolerance_intermediate
        solver_uncon.opts.gradient_norm_tolerance = solver.opts.gradient_norm_tolerance_intermediate
    else
        solver_uncon.opts.cost_tolerance = solver.opts.cost_tolerance
        solver_uncon.opts.gradient_norm_tolerance = solver.opts.gradient_norm_tolerance
    end

    return nothing
end

"Augmented Lagrangian step"
function step!(prob::Problem{T}, solver::AugmentedLagrangianSolver{T},
        unconstrained_solver::AbstractSolver) where T
    to = solver.stats[:timer]

    # Solve the unconstrained problem
    @timeit to "uncon solve" solve!(prob, unconstrained_solver)
    @timeit to "cost" J = cost(prob)

    # Outer loop update
    @timeit to "dual update" dual_update!(prob, solver)
    @timeit to "penalty update" penalty_update!(prob, solver)
    copyto!(solver.C_prev,solver.C)

    return J
end

"Evaluate maximum constraint violation as metric for Augmented Lagrangian solve convergence"
function evaluate_convergence(solver::AugmentedLagrangianSolver{T}) where T
    converged = false
    if solver.opts.kickout_max_penalty
        converged = (max_penalty(solver) == solver.opts.penalty_max ?
            begin @logmsg OuterLoop "Max Penalty"; true end : false)
    end
    converged || (solver.stats[:c_max][end] < solver.opts.constraint_tolerance ? true : false)
end

function record_iteration!(prob::Problem{T}, solver::AugmentedLagrangianSolver{T}, J::T,
        unconstrained_solver::AbstractSolver) where T
    c_max = max_violation(solver)

    solver.stats[:iterations] += 1
    solver.stats[:iterations_total] += unconstrained_solver.stats[:iterations]
    push!(solver.stats[:iterations_inner], unconstrained_solver.stats[:iterations])
    push!(solver.stats[:cost],J)
    push!(solver.stats[:c_max],c_max)
    push!(solver.stats[:penalty_max],max_penalty(solver))
    push!(solver.stats_uncon, copy(unconstrained_solver.stats))
    push!(solver.stats[:S_cond], maximum(map(S->cond(S.xx), solver.solver_uncon.S)))
    push!(solver.stats[:Q_cond], maximum(map(Q->cond(Q.uu), solver.solver_uncon.Q[1:end-1])))

    @logmsg OuterLoop :iter value=solver.stats[:iterations]
    @logmsg OuterLoop :total value=solver.stats[:iterations_total]
    @logmsg OuterLoop :cost value=J
    @logmsg OuterLoop :c_max value=c_max
end

function max_penalty(solver::AugmentedLagrangianSolver)
    maximum(maximum(solver.μ))
end

"Saturate a vector element-wise with upper and lower bounds"
saturate(input::AbstractVector{T}, max_value::T, min_value::T) where T = max.(min_value, min.(max_value, input))

"Dual update (first-order)"
function dual_update!(prob::Problem{T}, solver::AugmentedLagrangianSolver{T}) where T
    c = solver.C; λ = solver.λ; μ = solver.μ

    for k = 1:prob.N
        copyto!(λ[k],saturate(λ[k] + μ[k].*c[k], solver.opts.dual_max,
            solver.opts.dual_min))
        copyto!(λ[k].inequality,max.(0.0, λ[k].inequality))
    end

    # Update active set after updating multipliers (need to calculate c_max)
    update_active_set!(solver.active_set, solver.C, solver.λ)
end

"Penalty update (default) - update all penalty parameters"
function penalty_update!(prob::Problem{T}, solver::AugmentedLagrangianSolver{T}) where T
    μ = solver.μ
    for k = 1:prob.N
        copyto!(μ[k], saturate(solver.opts.penalty_scaling * μ[k], solver.opts.penalty_max, 0.0))
    end
end

"$(TYPEDEF) Augmented Lagrangian Objective: stores stage cost(s) and terminal cost functions"
struct AugmentedLagrangianObjective{T} <: AbstractObjective where T
    cost::CostTrajectory
    constraints::Constraints
    C::PartedVecTrajectory{T}  # Constraint values
    ∇C::PartedMatTrajectory{T} # Constraint jacobians
    λ::PartedVecTrajectory{T}  # Lagrange multipliers
    μ::PartedVecTrajectory{T}  # Penalty Term
    active_set::PartedVecTrajectory{Bool}  # Active set
end

function AugmentedLagrangianObjective(cost::CostTrajectory,constraints::Constraints,N::Int;
        μ_init::T=1.,λ_init::T=0.) where T
    # Get sizes
    n,m = get_sizes(cost)
    C,∇C,λ,μ,active_set = init_constraint_trajectories(constraints,n,m,N)
    AugmentedLagrangianObjective{T}(cost,constraint,C,∇C,λ,μ,active_set)
end

function AugmentedLagrangianObjective(cost::CostTrajectory,constraints::Constraints,
        λ::PartedVecTrajectory{T}; μ_init::T=1.) where T
    # Get sizes
    n,m = get_sizes(cost)
    N = length(λ)
    C,∇C,_,μ,active_set = init_constraint_trajectories(constraints,n,m,N)
    AugmentedLagrangianObjective{T}(cost,constraint,C,∇C,λ,μ,active_set)
end

getindex(obj::AugmentedLagrangianObjective,i::Int) = obj.cost[i]

"Generate augmented Lagrangian cost from unconstrained cost"
function AugmentedLagrangianObjective(prob::Problem{T},
        solver::AugmentedLagrangianSolver{T}) where T
    AugmentedLagrangianObjective{T}(prob.obj.cost,prob.constraints,solver.C,solver.∇C,solver.λ,solver.μ,solver.active_set)
end

"Generate augmented Lagrangian problem from constrained problem"
function AugmentedLagrangianProblem(prob::Problem{T},solver::AugmentedLagrangianSolver{T}) where T
    obj_al = AugmentedLagrangianObjective(prob,solver)
    prob_al = update_problem(prob,obj=obj_al,constraints=Constraints(prob.N),newProb=false)
end

"Evaluate maximum constraint violation"
function max_violation(solver::AugmentedLagrangianSolver{T}) where T
    c_max = 0.0
    C = solver.C
    N = length(C)
    for k = 1:N
        if length(C[k]) > 0
            c_max = max(norm(C[k].equality,Inf), c_max)
            if length(C[k].inequality) > 0
                c_max = max(pos(maximum(C[k].inequality)), c_max)
            end
        end
    end
    return c_max
end

function cost_expansion!(Q::ExpansionTrajectory{T},obj::AugmentedLagrangianObjective{T},
        X::VectorTrajectory{T},U::VectorTrajectory{T},dt::Vector{T}) where T
    N = length(X)

    cost_expansion!(Q, obj.cost, X, U, dt)

    for k = 1:N-1
        c = obj.C[k]
        λ = obj.λ[k]
        μ = obj.μ[k]
        a = active_set(c,λ)
        Iμ = Diagonal(a .* μ)
        jacobian!(obj.∇C[k],obj.constraints[k],X[k],U[k])
        cx = obj.∇C[k].x
        cu = obj.∇C[k].u

        # Second Order pieces
        Q[k].xx .+= cx'Iμ*cx
        Q[k].uu .+= cu'Iμ*cu
        Q[k].ux .+= cu'Iμ*cx

        # First order pieces
        g = (Iμ*c + λ)
        Q[k].x .+= cx'g
        Q[k].u .+= cu'g
    end

    c = obj.C[N]
    λ = obj.λ[N]
    μ = obj.μ[N]
    a = active_set(c,λ)
    Iμ = Diagonal(a .* μ)
    cx = obj.∇C[N]

    jacobian!(cx,obj.constraints[N],X[N])

    # Second Order pieces
    Q[N].xx .+= cx'Iμ*cx

    # First order pieces
    Q[N].x .+= cx'*(Iμ*c + λ)

    return nothing
end

function cost_expansion_sqrt!(Q::ExpansionTrajectory{T},obj::AugmentedLagrangianObjective{T},
        X::VectorTrajectory{T},U::VectorTrajectory{T},dt::Vector{T}) where T
    N = length(X)

    cost_expansion_sqrt!(Q, obj.cost, X, U, dt)

    for k = 1:N-1
        c = obj.C[k]
        λ = obj.λ[k]
        μ = obj.μ[k]
        a = active_set(c,λ)
        Iμ = Diagonal(a .* μ)
        Iμ_sqrt = Diagonal(a.* sqrt.(μ))
        jacobian!(obj.∇C[k],obj.constraints[k],X[k],U[k])
        cx = obj.∇C[k].x
        cu = obj.∇C[k].u

        # Second Order pieces
        chol_plus!(Q[k].xx,Iμ_sqrt*cx)
        chol_plus!(Q[k].uu,Iμ_sqrt*cu)

        # First order pieces
        g = (Iμ*c + λ)
        Q[k].x .+= cx'g
        Q[k].u .+= cu'g
    end

    c = obj.C[N]
    λ = obj.λ[N]
    μ = obj.μ[N]
    a = active_set(c,λ)
    Iμ = Diagonal(a .* μ)
    Iμ_sqrt = Diagonal(a.* sqrt.(μ))

    cx = obj.∇C[N]

    jacobian!(cx,obj.constraints[N],X[N])

    # Second Order pieces
    chol_plus!(Q[N].xx,Iμ_sqrt*cx)

    # First order pieces
    Q[N].x .+= cx'*(Iμ*c + λ)

    return nothing
end

function update_active_set!(obj::AugmentedLagrangianObjective{T},tol::T=0.0) where T
    update_active_set!(obj.active_set,obj.C,obj.λ,tol)
end

"Cost function terms for Lagrangian and quadratic penalty"
function aula_cost(a::AbstractVector{Bool},c::AbstractVector{T},λ::AbstractVector{T},μ::AbstractVector{T}) where T
    λ'c + 1/2*c'Diagonal(a .* μ)*c
end

function stage_constraint_cost(c,λ,μ,
        a,x::AbstractVector{T},u::AbstractVector{T}) where T
    aula_cost(a,c,λ,μ)
end

function stage_constraint_cost(c,λ,μ,
        a,x::AbstractVector{T}) where T
    aula_cost(a,c,λ,μ)
end

function cost(obj::AugmentedLagrangianObjective{T}, X::VectorTrajectory{T}, U::VectorTrajectory{T}, dt::Vector{T}) where T <: AbstractFloat
    N = length(X)
    J = cost(obj.cost,X,U,dt)

    update_constraints!(obj.C,obj.constraints, X, U)
    update_active_set!(obj)

    Jc = 0.0
    for k = 1:N-1
        Jc += stage_constraint_cost(obj.C[k], obj.λ[k], obj.μ[k],obj.active_set[k],X[k],U[k])
    end

    Jc += stage_constraint_cost(obj.C[N], obj.λ[N], obj.μ[N],obj.active_set[N],X[N])

    return J + Jc
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
include("try.jl")

 true_list = ["../PDFs/Data/XOR-0.0-0.dens","../PDFs/Data/AND-0.0-0.dens","../PDFs/Data/UNQ-0.0-0.dens","../PDFs/Data/RDN-0.0-0.dens","../PDFs/Data/XORAND-0.0-0.dens","../PDFs/Data/RDNXOR-0.0-0.dens","../PDFs/Data/RDNUNQXOR-0.0-0.dens"]
# true_list = ["../PDFs/Data/XOR-0.0-0.dens"]
# true_list = ["../PDFs/Data/AND-0.0-0.dens"]
# true_list = ["../PDFs/Data/UNQ-0.0-0.dens"]
# true_list = ["../PDFs/Data/RDN-0.0-0.dens"]
# true_list = ["../PDFs/Data/XORAND-0.0-0.dens"]
# true_list = ["../PDFs/Data/RDNXOR-0.0-0.dens"]
# true_list = ["../PDFs/Data/RDNUNQXOR-0.0-0.dens"]

 list = ["../PDFs/Data/XOR-0.1-10000.dens","../PDFs/Data/AND-0.1-10000.dens","../PDFs/Data/UNQ-0.1-10000.dens","../PDFs/Data/RDN-0.1-10000.dens","../PDFs/Data/XORAND-0.1-10000.dens","../PDFs/Data/RDNXOR-0.1-10000.dens","../PDFs/Data/RDNUNQXOR-0.1-10000.dens"]
# list = ["../PDFs/Data/XOR-0.1-10000.dens"]
# list = ["../PDFs/Data/AND-0.1-10000.dens"]
# list = ["../PDFs/Data/UNQ-0.1-10000.dens"]
# list = ["../PDFs/Data/RDN-0.1-10000.dens"]
# list = ["../PDFs/Data/XORAND-0.1-10000.dens"]
# list = ["../PDFs/Data/RDNXOR-0.1-10000.dens"]
# list = ["../PDFs/Data/RDNUNQXOR-0.1-10000.dens"]
list_new=String[]
for s in 1:length(list)
    for n in 1:101
        if n == 1
            push!(list_new, list[s])
        elseif n == 101
            continue
        end
        r = n - 1
        push!(list_new, string(list[s],"-$r"))
    end#^n
end#^s

# for s in 1:length(list)
#     for n in 1:10
#         push!(list_new, string(list[s],"-$n"))
#     end#^n
# end#^s

# solvers_list = [:Mosek, :Ipopt, :Knitro_Ip, :Knitro_IpCG, :Knitro_AS, :Knitro_SQP, :Cvxopt]
solvers_list = [:Mosek,:Ipopt]
tmpFloatDatatype = Float64
for s in 1:length(list_new)
    for solver in solvers_list
        println("///////////////////////////")
        println("Starting computations of the distribution: ", list_new[s], " with ", solver, " and prc ", tmpFloatDatatype)
        println("///////////////////////////")
        test(solver,list_new[s], true_list[convert(Int32,ceil(s/101))], tmpFloatDatatype)
        println("///////////////////////////")
        println("Finishing computations of the distribution: ", list_new[s]," with ", solver, " and prc ", tmpFloatDatatype)
        println("///////////////////////////")
    end
end

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
# ---
# title: 864. Shortest Path to Get All Keys
# id: problem864
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Heap, Breadth-first Search
# link: <https://leetcode.com/problems/shortest-path-to-get-all-keys/description/>
# hidden: true
# ---
#
# We are given a 2-dimensional `grid`. `"."` is an empty cell, `"#"` is a wall,
# `"@"` is the starting point, (`"a"`, `"b"`, ...) are keys, and (`"A"`, `"B"`,
# ...) are locks.
#
# We start at the starting point, and one move consists of walking one space in
# one of the 4 cardinal directions.  We cannot walk outside the grid, or walk
# into a wall.  If we walk over a key, we pick it up.  We can't walk over a lock
# unless we have the corresponding key.
#
# For some 1 <= K <= 6, there is exactly one lowercase and one uppercase letter
# of the first `K` letters of the English alphabet in the grid.  This means that
# there is exactly one key for each lock, and one lock for each key; and also
# that the letters used to represent the keys and locks were chosen in the same
# order as the English alphabet.
#
# Return the lowest number of moves to acquire all keys.  If it's impossible,
# return `-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: ["@.a.#","###.#","b.A.B"]
#     Output: 8
#
#
# **Example 2:**
#
#
#
#     Input: ["@..aA","..B#.","....b"]
#     Output: 6
#
#
#
#
# **Note:**
#
#   1. `1 <= grid.length <= 30`
#   2. `1 <= grid[0].length <= 30`
#   3. `grid[i][j]` contains only` '.'`, `'#'`, `'@'`, `'a'-``'f``'` and `'A'-'F'`
#   4. The number of keys is in `[1, 6]`.  Each key has a different letter and opens exactly one lock.
#
#
## @lc code=start
using LeetCode

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


text:
mutable struct LassoSCOP <: LassoRoot
    ## It's just a collection of data.

    A::Matrix # Feature matrix
    y::Matrix # label vector
    μ::Matrix # feature mean
    u::Number # label mean
    Z::Matrix # Starndardized Matrix
    l::Matrix # Zero mean label

    OptModel::Model  # The JuMP model for getting it right.
    λ::Float64       # the regularization parameter.

    LassoPath::Matrix # going along a fixed row is the fixing
    # the feature while varying the lambda quantity.
    λs::Vector # the lambda values.

    # Paramters for the iterator
    Tol::Float64
    λMin::Float64

    function LassoSCOP(A::Matrix, y::Union{Matrix, Vector}, λ::Float64=0.0)
        A = copy(A)
        m, _ = size(A)
        @assert size(y, 1) == m ""*
        "The rows of X should match of the columns of y, but the size of"*
        string("X, Y is: ", size(A), " ", size(y))

        y = copy(y)
        if ndims(y) == 1
            y = reshape(y, (length(y), 1))
        end
        μ = mean(A, dims=1)::Matrix
        u = mean(y)
        Z = A .- μ
        l = y .- u
        OptModel = MakeLassoOptimizationProblem(Z, l, λ)
        this = new(A, y, μ, u, Z, l, OptModel, λ)
        this.Tol = 1e-4
        this.λMin = 1e-8
        return this
    end
end

# ==============================================================================
# Type specific functions
# ==============================================================================

function Changeλ!(this::LassoSCOP, λ)
    """
        Change the Lasso regularizer of the current model, this mutates
        the state of the object.
    """
    model = this.OptModel
    x = model[:x]  # objects can be indexed with symbols!
    η = model[:η]
    y = this.l
    A = this.Z
    @objective(model, Min, λ*sum(η) + sum((A*x - y).^2))

end

function SolveLasso(
        this::LassoSCOP;
        x0::T=nothing
    )::Vector{Float64} where {T <: Union{Vector{Float64}, Nothing}}
    """
        Solve for the weights of the current model,
        given the current configuration of the model.
    """
    if !(x0 === nothing)
        setvalue.(this.OptModel[:x], x0) # warm start!
    end
    optimize!(this.OptModel)
    TernimationStatus = termination_status(this.OptModel)
    # @assert TernimationStatus == MOI.OPTIMAL "Terminated with non-optimal value when solving for x. "*
    # string("The status is: ", TernimationStatus)*"\n this is the results \n $(OptResults)"

    if !(TernimationStatus == MOI.OPTIMAL)
        Warn("\nWarning: convergence status for solver: $(TernimationStatus)")
        Warn("Current Value λ: $(this.λ)")
        PrintTitle("Here is the summary for the solution: ")
        display(solution_summary(this.OptModel))
    end
    return value.(this.OptModel[:x])
end

function Base.iterate(this::LassoSCOP)
    """
        Iterate through the parameters for the Lasso program:
            * λ
            * solutions
    """
    u = this.u
    A = this.Z
    y = this.l

    function λMax(A, y)
        ToMax = A'*(y .- u)
        ToMax *= 2
        return maximum(abs.(ToMax))
    end
    λ = λMax(A, y)
    Changeλ!(this, λ)
    x = SolveLasso(this)

    return (x, λ), (x, λ)
end

function Base.iterate(this::LassoSCOP, state::Tuple{Vector{Float64}, Float64})
    x, λ = state
    λ /= 2
    setvalue.(this.OptModel[:x], x)
    Changeλ!(this, λ)
    y = SolveLasso(this)

    if norm(x - y, Inf) <= this.Tol || λ <= this.λMin
        return nothing
    end

    return (y, λ), (y, λ)

end


text:
#Notes:
#CovarianceOperator is a custom hermitian operator that explicitly defines the multiplication operation
#so that you never need the whole covariance matrix in memory
#Why is there a tmp here?
function covmul!(y, Xc, n, x)
    y = Xc'*(Xc*x)./(n-1)
end

mutable struct CovarianceOperator{T<:Real} <: AbstractLinearOperator{T}
    p::Int
    n::Int
    mul!::Function
    _tmp::Nullable{Array{T}}
end

function CovarianceOperator(X)
    T = eltype(X)
    n, p = size(X)
    Xm = mean(X, dims=1)
    ml! = (y, _, x) -> covmul!(y, X.-Xm, n, x)
    CovarianceOperator{T}(p, n, ml!, nothing)
end

Base.convert(::Type{LinearOperator}, A::CovarianceOperator) = convert(LinearOperator{eltype(A)}, A)
Base.convert(::Type{LinearOperator{T}}, A::CovarianceOperator{T}) where T = LinearOperator{T}(A.p, A.p, A.mul!, A.mul!, A._tmp)
Base.transpose(A::CovarianceOperator) = A
LinearAlgebra.adjoint(A::CovarianceOperator) = A
LinearAlgebra.ishermitian(A::CovarianceOperator) = true
LinearAlgebra.issymmetric(A::CovarianceOperator) = isreal(A)
Base.size(A::CovarianceOperator) = (A.p, A.p)
Base.size(A::CovarianceOperator, dim::Integer) = (dim == 1 || dim == 2) ? A.p : 1
function Base.show(io::IO, C::CovarianceOperator)
    print(io, "CovarianceOperator(size = ", C.p)
    print(io, "x", C.p)
    print(io, ", samples = ", C.n)
    print(io, ")")
end

function crosscovmul!(y, Xc, Yc, n, x)
    y = Xc'*(Yc*x)./(n-1)
end

mutable struct CrossCovarianceOperator{T<:Real} <: AbstractLinearOperator{T}
    p::Int
    q::Int
    n::Int
    mul!::Function
    mulc!::Function
    _tmp::Nullable{Array{T}}
end

function CrossCovarianceOperator(X, Y)
    if (size(X) == size(Y)) && isapprox(X, Y)
        return CovarianceOperator(X)
    end
    Tx = eltype(X)
    Ty = eltype(Y)
    @assert Tx == Ty
    nx, p = size(X)
    ny, q = size(Y)
    @assert nx == ny
    Xm = mean(X,dims=1)
    Ym = mean(Y,dims=1)
    ml! = (y, _, x) -> crosscovmul!(y, X.-Xm, Y.-Ym, nx, x)
    mulc! = (y, _, x) -> crosscovmul!(y, Y.-Ym, X.-Xm, nx, x)
    CrossCovarianceOperator{Tx}(p, q, nx, ml!, mulc!, nothing)
end

Base.convert(::Type{LinearOperator}, A::CrossCovarianceOperator) = convert(LinearOperator{eltype(A)}, A)
Base.convert(::Type{LinearOperator{T}}, A::CrossCovarianceOperator{T}) where T = LinearOperator{T}(A.p, A.q, A.mul!, A.mulc!, A._tmp)
Base.transpose(A::CrossCovarianceOperator{T}) where {T} = CrossCovarianceOperator{T}(A.q, A.p, A.mulc!, A.mul!, nothing)
LinearAlgebra.adjoint(A::CrossCovarianceOperator{T}) where {T} = CrossCovarianceOperator{T}(A.q, A.p, A.mulc!, A.mul!, nothing)
LinearAlgebra.ishermitian(A::CrossCovarianceOperator) = false
LinearAlgebra.issymmetric(A::CrossCovarianceOperator) = false
Base.size(A::CrossCovarianceOperator) = (A.p, A.q)
#this is a ``read only'' shortcircuiting way to match the dim - see julia shortcircuiting rules to understand
Base.size(A::CrossCovarianceOperator, dim::Integer) = (dim == 1 && return A.p) || (dim == 2) ? A.q : 1
function Base.show(io::IO, C::CrossCovarianceOperator)
    print(io, "CrossCovarianceOperator(size = ", C.p)
    print(io, "x", C.q)
    print(io, ", samples = ", C.n)
    print(io, ")")
end


text:
function toh06c()
#
#  generalization of K.C. Toh and L. N. Trefethen
#
    p = ones(6)
    p = conv(p,p)
    p = conv(p,p);
    z = [0.50000000000000 + 0.86602540378444im
         0.50000000000000 - 0.86602540378444im
        -1.00000000000000
        -0.50000000000000 + 0.86602540378444im
        -0.50000000000000 - 0.86602540378444im]
    z = [z 4*ones(5)]

    p, PolyZeros(z)
end


text:
@testset "Composition" begin
  c₀ = Composition(1,1,1)
  c₁ = Composition(1,2,3)
  c₂ = Composition(10,20,30)
  @test c₁ == c₂
  @test c₁ - c₂ == c₀
  @test norm(c₀) ≈ 0.
  @test isapprox(distance(c₁, c₁), 0., atol=1e-6)
end


text:
# Stacks

mutable struct Stack{T}
  idx::Int
  data::Vector{T}
end

Stack(data::Vector{T}) where T =
  Stack{T}(length(data), data)

function Base.pop!(stk::Stack)
  i = stk.idx
  stk.idx = i == 1 ? length(stk.data) : i-1
  @inbounds return stk.data[i]
end

xstack(T) = (Vector{T}, Expr(:call, Vector{T}))

function _push!(a::Vector{T}, x::T) where T
  Base._growend!(a, 1)
  @inbounds a[end] = x
  return
end

# Emit

function alphauses(ir, bi)
  us = []
  for i = range(ir.cfg.blocks[bi]), u in userefs(ir.stmts[i])
    u[] isa Alpha && push!(us, SSAValue(u[].id))
  end
  return unique(us)
end

xtuple(xs...) = xcall(:tuple, xs...)

concrete(T::DataType) = T
concrete(::Type{Type{T}}) where T = typeof(T)

# TODO: combine stacks by type
function forward_stacks!(adj, F)
  stks, recs = [], []
  for fb = 1:length(adj.perm)
    for α in alphauses(adj.back, invperm(adj.perm)[fb])
      if fb == 1
        push!(recs, α)
      else
        T = exprtype(adj.forw, α)
        stk = insert_node!(adj.forw, 1, xstack(T)...)
        push!(recs, stk)
        insert_blockend!(adj.forw, blockidx(adj.forw, α.id), Any, xcall(Zygote, :_push!, stk, α))
      end
      push!(stks, (invperm(adj.perm)[fb], alpha(α)))
    end
  end
  args = [Argument(i) for i = 3:length(adj.forw.argtypes)]
  T = Tuple{concrete.(exprtype.(Ref(adj.forw), (args..., recs...)))...}
  isconcretetype(T) || (T = Any)
  rec = insert_node!(adj.forw, length(adj.forw.stmts), T,
                     xtuple(args..., recs...))
  if usetyped
    rec = insert_node!(adj.forw, length(adj.forw.stmts), J{F,T},
                       Expr(:call, J{F,T}, rec))
  else
    rec = insert_node!(adj.forw, length(adj.forw.stmts), Any,
                       Expr(:call, J{F}, rec))
  end
  ret = xtuple(adj.forw.stmts[end].val, rec)
  R = exprtype(adj.forw, adj.forw.stmts[end].val)
  ret = insert_node!(adj.forw, length(adj.forw.stmts), Tuple{R,J{F,T}}, ret)
  adj.forw.stmts[end] = ReturnNode(ret)
  forw = compact!(adj.forw)
  return forw, stks
end

function reverse_stacks!(adj, stks, nargs)
  ir = adj.back
  t = insert_node!(ir, 1, Any, xcall(Base, :getfield, Argument(1), QuoteNode(:t)))
  for b = 1:length(ir.cfg.blocks)
    repl = Dict()
    for (i, (b′, α)) in enumerate(stks)
      b == b′ || continue
      loc = max(2,afterphi(ir, range(ir.cfg.blocks[b])[1]))
      if adj.perm[b′] == 1
        val = insert_node!(ir, loc, Any, xcall(:getindex, t, i+nargs))
      else
        stk = insert_node!(ir, 1, Any, xcall(:getindex, t, i+nargs))
        stk = insert_node!(ir, 1, Any, xcall(Zygote, :Stack, stk))
        val = insert_node!(ir, loc, Any, xcall(:pop!, stk))
      end
      repl[α] = val
    end
    for i in range(ir.cfg.blocks[b]), u in userefs(ir.stmts[i])
      if u.stmt == Expr(:call, :Δ)
        u.stmt = Argument(2)
      elseif u[] isa Argument
        x = insert_node!(ir, i, Any, xcall(:getindex, t, u[].n-2))
        u[] = x
      elseif haskey(repl, u[])
        u[] = repl[u[]]
      else continue
      end
      ir.stmts[i] = u.stmt
    end
  end
  return compact!(ir)
end

function stacks!(adj, T)
  forw, stks = forward_stacks!(adj, T)
  back = reverse_stacks!(adj, stks, length(forw.argtypes)-2)
  return forw, back
end

varargs(m::Method, n) = m.isva ? n - m.nargs + 1 : nothing

function _lookup_grad(T)
  (m = meta(T)) == nothing && return
  usetyped && m.ret == Union{} && return
  va = varargs(m.method, length(T.parameters))
  forw, back = stacks!(Adjoint(IRCode(m), varargs = va), T)
  # verify_ir(forw)
  # verify_ir(back)
  m, forw, back
end


text:
"""

Generates a unique ID for a proper DNA sequence (consisting of A,C,T,G
bases).  If there are other letters in the sequence it returns
UInt(0).

"""
function gen_id(seq::Vector{UInt8})
    if length(seq) >= 32
        error("Sequence is too long ($(length(seq)) bases when only 31 are allowed) to assign an ID")
    end
    i = UInt(0)
    for nt in seq
        i <<= 2
        if nt == UInt8('A')
            i += 0
        elseif nt == UInt8('C')
            i += 1
        elseif nt == UInt8('T')
            i += 2
        elseif nt == UInt8('G')
            i += 3
        else
            return UInt(0)
        end
    end
    # We add 1 at the end in case there were only As in the sequence,
    # which would result in UInt(0) and be indistinguishable from an
    # inproper sequence.
    return i+1
end

gen_id(seq::String) = gen_id(Vector{UInt8}(seq))

"""
the inverse of gen_id, given `i` it returns an associated sequence
"""
function idtoname!(name::Vector{UInt8},i::UInt)
    i-=1
    for n in length(name):-1:1
        if xor(i&0b11,0b00) == 0
            name[n]=UInt8('A')
        elseif xor(i&0b11,0b01) == 0
            name[n]=UInt8('C')
        elseif xor(i&0b11,0b10) == 0
            name[n]=UInt8('T')
        elseif xor(i&0b11,0b11) ==0
            name[n]=UInt8('G')
        end
        i>>=2
    end
    return name
end

function idtoname(i::UInt,len::Int)
    name = zeros(UInt8,len)
    idtoname!(name,i)
    return String(name)
end

"""
merge the contents of `files` into `output`
"""
function catfiles(output,files)
    fout = open(output,"w")
    for f in files
        write(fout,read(f))
    end
    close(fout)
end

max_stars_repo_name: pwl/FASTQDemultiplexer

max_issues_repo_name: pwl/FASTQDemultiplexer

max_forks_repo_name: pwl/FASTQDemultiplexer


text:
# script to parse oneAPI headers and generate Julia wrappers

#
# Parsing
#

using Clang.Generators
using JuliaFormatter

function wrap(name, headers...; library="lib$name", defines=[], include_dirs=[])
    options =  load_options(joinpath(@__DIR__, "wrap.toml"))
    options["general"]["library_name"] = library
    options["general"]["output_file_path"] = "lib$(name).jl"

    args = get_default_args()
    ctx = create_context([headers...], args, options)

    build!(ctx, BUILDSTAGE_NO_PRINTING)

    # tweak exprs
    for node in get_nodes(ctx.dag)
        exprs = get_exprs(node)
        for (i, expr) in enumerate(exprs)
            exprs[i] = add_check_pass(expr)
        end
    end

    build!(ctx, BUILDSTAGE_PRINTING_ONLY)

    return options["general"]["output_file_path"]
end

## rewrite passes

# insert `@checked` before each function with a `ccall` returning a checked type`
checked_types = [
    "ze_result_t",
]
function add_check_pass(x::Expr)
    Meta.isexpr(x, :function) || return x
    body = x.args[2].args[1]
    if Meta.isexpr(body, :macrocall)  # `@ccall`
        ret_type = string(body.args[3].args[2])
        if ret_type in checked_types
            return Expr(:macrocall, Symbol("@checked"), nothing, x)
        end
    else
        @assert Meta.isexpr(body, :call)  # `ccall`
        ret_type = string(body.args[3])
        if ret_type in checked_types
            return Expr(:macrocall, Symbol("@checked"), nothing, x)
        end
    end
    return x
end

#
# Main application
#

using oneAPI_Level_Zero_Headers_jll

function process(name, headers...; modname=name, kwargs...)
    output_file = wrap(name, headers...; kwargs...)

    let file = output_file
        text = read(file, String)

        ## header

        squeezed = replace(text, "\n\n\n"=>"\n\n")
        while length(text) != length(squeezed)
            text = squeezed
            squeezed = replace(text, "\n\n\n"=>"\n\n")
        end
        text = squeezed

        write(file, text)
    end

    ## manual patches

    patchdir = joinpath(@__DIR__, "patches", name)
    if isdir(patchdir)
        for entry in readdir(patchdir)
            if endswith(entry, ".patch")
                path = joinpath(patchdir, entry)
                run(`patch -p1 -i $path`)
            end
        end
    end

    ## move to destination and format

    let src = output_file
        dst = joinpath(dirname(@__DIR__), "lib", modname, src)
        cp(src, dst; force=true)
        format(dst, YASStyle(), always_use_return=false)
    end

    return
end

function main()
    process("ze", oneAPI_Level_Zero_Headers_jll.ze_api; library="libze_loader", modname="level-zero")
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Deprecated functions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                              Introduced in v1.0
################################################################################

@deprecate create_rotation_matrix angle_to_dcm

@deprecate zeros(T::Type{Quaternion{P}}) where P zero(T)
@deprecate zeros(T::Type{Quaternion}) zero(T)
@deprecate zeros(q::Quaternion{T}) where T zero(q)


text:
# Copyright 2020 Google LLC

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
@enum Op nop acc jmp

mutable struct Instruction
    op::Op
    arg::Int
end

mutable struct ComputerState
    accumulator::Int
    pc::Int
    memory::AbstractVector{Instruction}
end

function parseline(line::AbstractString)::Instruction
    if startswith(line, "nop")
        return Instruction(nop, parse(Int, line[4:end]))
    elseif startswith(line, "acc")
        return Instruction(acc, parse(Int, line[4:end]))
    elseif startswith(line, "jmp")
        return Instruction(jmp, parse(Int, line[4:end]))
    else
        throw(ErrorException("Unrecognized op"))
    end
end

function initializecomputer(lines::AbstractVector{String})::ComputerState
    memory = map(parseline, lines)
    return ComputerState(0, 1, memory)
end

function execute(state::ComputerState, visited::AbstractSet{Int}, termination_predicate)
    while true
        instr = state.memory[state.pc]
        push!(visited, state.pc)
        if instr.op == nop
            state.pc += 1
        elseif instr.op == acc
            state.pc += 1
            state.accumulator += instr.arg
        elseif instr.op == jmp
            state.pc += instr.arg
        end
        if termination_predicate(state)
            return
        end
    end
end


text:
function Explore(A,b,C,d,XLB,XUB,VariableDef,YR;TimeLim = 1e10,Cuts = [])
        const EPS = 1e-5;
        t0 = time_ns();
        nObj = length(YR)

        MMModel= MinMaxModel(A,b,C,d,XLB,XUB,VariableDef,YR,Cuts = Cuts)

        # Time commands ---------------------------------------------------
        if (TimeLim-(time_ns() - t0)/1e9) < 0.001;
            return [],[],[],false
        end
        # -----------------------------------------------------------------

        setsolver(MMModel,CplexSolver(CPXPARAM_Threads = 1,CPXPARAM_ScreenOutput = 0,CPXPARAM_TimeLimit = TimeLim-(time_ns() - t0)/1e9));
        status = solve(MMModel);

        # Time commands ---------------------------------------------------
        if (status != :Optimal) | (TimeLim-(time_ns() - t0)/1e9 < 0.001)

            return [],[],[],false
        end
        # -----------------------------------------------------------------

        Var = getvalue(getVar(MMModel,:X))
        Y_bar1 = zeros(size(YR))
        for j = 1:nObj
            Y_bar1[j] = sum(Var.*C[j,:]) + d[j];
        end

        Y = YR + maximum(abs.(Y_bar1-YR))

        return Y_bar1,Var,Y,true
end


text:
using WaterLily
using LinearAlgebra: norm2
include("ThreeD_Plots.jl")

function TGV(p=6,Re=1e5)
    # Define vortex size, velocity, viscosity
    L = 2^p; U = 1; ν = U*L/Re

    # Taylor-Green-Vortex initial velocity field
    function uλ(i,vx)
        x,y,z = @. (vx-1.5)*π/L                # scaled coordinates
        i==1 && return -U*sin(x)*cos(y)*cos(z) # u_x
        i==2 && return  U*cos(x)*sin(y)*cos(z) # u_y
        return 0.                              # u_z
    end

    # Initialize simulation
    return Simulation((L+2,L+2,L+2),zeros(3),L;U,uλ,ν)
end

function ω_mag_data(sim)
    # plot the vorticity modulus
    @inside sim.flow.σ[I] = WaterLily.ω_mag(I,sim.flow.u)*sim.L/sim.U
    return @view sim.flow.σ[2:end-1,2:end-1,2:end-1]
end

sim,fig = volume_video!(TGV(),ω_mag_data,name="TGV.mp4",duration=10)


text:
using Weave
files = readdir(pwd())
filename(path) = splitext(basename(path))[1]

for f in files
    if occursin(".ipynb", f) && f != ".ipynb_checkpoints"
        convert_doc(f, filename(f) * ".jmd")
    end
end


text:
# -----------------------------------------------------------------------------
# Scaled-Diagonally-dominant bridge [MOI]
# -----------------------------------------------------------------------------

"""
Converts a Semidefinite Program (SDP) into a Scaled=Diagonally-Dominant Program
(SDDP). Positive-semidefiniteness constraints are relaxed into Scaled-diagonal-
dominance constraints.
"""
struct SDDPBridge{T,F,G} <: MOIB.Constraint.AbstractBridge
        side_dimension
        constraint
end

function MOIB.Constraint.bridge_constraint(
        ::Type{SDDPBridge{T,F,G}},
        model::MOI.ModelLike,
        f::MOI.AbstractVectorFunction,
        s::MOI.PositiveSemidefiniteConeTriangle,
) where {T,F,G}
        n = s.side_dimension
        matrix = square_a_triangle(f, n)

        for j in 1:n, i in 1:j-1
                zs = MOI.add_variables(model, 3)
                vs = MOI.SingleVariable.(zs)
                MOIU.operate!(-, T, matrix[i, i], vs[1])
                MOIU.operate!(-, T, matrix[i, j], vs[2])
                MOIU.operate!(-, T, matrix[j, j], vs[3])

                g = MOIU.operate(vcat, T, vs[1], vs[3], T(√2) * vs[2])
                l = MOI.add_constraint(model, g, MOI.RotatedSecondOrderCone(3))
        end

        # Handle inadvertent PSD constraints on a single variable
        e = undef
        if n == 1
                g = matrix[1]
                e = MOI.add_constraint(model, g, MOI.GreaterThan(zero(T)))
        else
                g = MOIU.operate(vcat, T, matrix...)
                e = MOI.add_constraint(model, g, MOI.Zeros(n * n))
        end

        return SDDPBridge{T,F,G}(n, e)
end

function MOI.supports_constraint(
        ::Type{<:SDDPBridge},
        ::Type{<:MOI.AbstractVectorFunction},
        ::Type{<:MOI.PositiveSemidefiniteConeTriangle},
)
        return true
end

function MOIB.added_constrained_variable_types(::Type{<:SDDPBridge})
        return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:SDDPBridge{T,F,G}}) where {T,F,G}
        return [
                (G, MOI.RotatedSecondOrderCone)
                (G, MOI.Zeros)
                (F, MOI.GreaterThan{T})
        ]
end
function MOIB.Constraint.concrete_bridge_type(
        ::Type{<:SDDPBridge{T}},
        G::Type{<:MOI.AbstractVectorFunction},
        ::Type{MOI.PositiveSemidefiniteConeTriangle},
) where {T}
        S = MOIU.scalar_type(G)
        F = MOIU.promote_operation(-, T, S, MOI.SingleVariable)
        return SDDPBridge{T,F,G}
end

function MOI.get(ml::MOI.ModelLike, cd::MOI.ConstraintDual, bridge::SDDPBridge)
        dual = MOI.get(ml, cd, bridge.constraint)
        side = bridge.side_dimension
        mat = reshape(dual, side, side)
        [(i == j) ? mat[i, j] : mat[i, j] / 2 for j in 1:side for i in 1:j]
end


text:
function bubble_sort(x) #function declaration
	for i in range(1,length(x)-1) #outer for loop
		for j in range(1,length(x)-1) #inner for loop
            		if x[j]>x[j+1]  #comparison to test acesending
                		t=x[j]       #exchange of variable
                		x[j]=x[j+1]  # ...
                		x[j+1]=t     #using a second one
            		end
    		end
    end
    for i in range(1,length(x)) #loop for printing
        println(x[i])
    end
end

x = [1,4,6,8,7,9,5,6]
bubble_sort(x)

max_stars_repo_name: rajatenzyme/Coding-Journey-

max_issues_repo_name: rajatenzyme/Coding-Journey-

max_forks_repo_name: rajatenzyme/Coding-Journey-


text:
import Avalon: convNd

@testset "conv" begin
    kwargs = (padding=1, stride=2)
    for (x_sz, w_sz) in [
            ((8, 3, 10), (3, 3, 1)),              # 1D
            ((8, 8, 3, 10), (3, 3, 3, 1)),        # 2D
        ]
        x = rand(Float32, x_sz); w = rand(Float32, w_sz)
        test_rrule(convNd, x, w; rtol=1e-3, atol=1e-3, check_inferred=false)
        test_rrule(convNd, x, w; fkwargs=kwargs, rtol=1e-3, atol=1e-3, check_inferred=false)
    end

    # 3D is very noisy for numeric gradient, using larger atol & rtol for it
    x = rand(Float32, 8, 8, 8, 3, 1); w = rand(Float32, 3, 3, 3, 3, 1)
    test_rrule(convNd, x, w; rtol=1e-1, atol=1e-1, check_inferred=false)
    test_rrule(convNd, x, w; fkwargs=kwargs, rtol=1e-1, atol=1e-1, check_inferred=false)
end

@testset "pooling" begin
    x = rand(7, 7, 3, 10);
    test_rrule(maxpool2d, x, 2; check_inferred=false)
    test_rrule(maxpool2d, x, 2; fkwargs=(stride=2,), check_inferred=false)
end


text:
module HomotopyContinuation

export ModelKit

using Arblib: Arblib
import DelimitedFiles
using DynamicPolynomials: @polyvar
import ElasticArrays: ElasticArray
import FiniteDiff
import LinearAlgebra
import LoopVectorization
using LRUCache: LRU
import MixedSubdivisions
import MultivariatePolynomials
const MP = MultivariatePolynomials
using Parameters: @unpack
import ProgressMeter
import Random
import Printf
import PrettyTables
using Reexport: @reexport
import SemialgebraicSets
import StructArrays
import TreeViews

const LA = LinearAlgebra

include("DoubleDouble.jl")
using .DoubleDouble

include("interval_arithmetic.jl")
using .IntervalArithmetic

include("ModelKit.jl")
export @polyvar

using ProjectiveVectors: PVector, dims, dimension_indices

const COMPILE_DEFAULT = Ref(:mixed)

export set_default_compile

"""
    set_default_compile(mode::Symbol)

Set the default value for the `compile` flag in [`solve`](@ref) and other functions.
Possible values are `:mixed` (default), `:all` and `:none`.
"""
function set_default_compile(mode::Symbol)
    mode ∈ [:mixed, :all, :none] ||
        error("Invalid value `:$mode`, valid values are `:mixed`, `:all`, `:none`.")
    COMPILE_DEFAULT[] = mode
end

include("utils.jl")
include("norm.jl")
include("voronoi_tree.jl")
include("unique_points.jl")
include("linear_algebra.jl")
include("linear_subspaces.jl")
include("systems.jl")
include("homotopies.jl")
include("predictor.jl")
include("newton_corrector.jl")
include("newton.jl")
include("tracker.jl")
include("valuation.jl")
include("path_result.jl")
include("endgame_tracker.jl")
include("path_info.jl")
include("total_degree.jl")
include("binomial_system.jl")
include("polyhedral.jl")
include("overdetermined.jl")
include("result.jl")
include("solve.jl")
include("monodromy.jl")
include("witness_set.jl")
include("certification.jl")
include("semialgebraic_sets.jl")

function __init__()
    # Overwrite default IJulia behaviour. See
    # https://github.com/timholy/ProgressMeter.jl/issues/162
    if isdefined(ProgressMeter, :ijulia_behavior)
        ProgressMeter.ijulia_behavior(:clear)
    end
end

include("precompile.jl")
_precompile_()

end #module


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(BinaryProvider.safe_isfile), String})
    precompile(Tuple{typeof(BinaryProvider.info_onchange), String, String, Int64})
    precompile(Tuple{typeof(BinaryProvider.detect_libgfortran_abi), String, BinaryProvider.Windows})
    precompile(Tuple{typeof(BinaryProvider.detect_libgfortran_abi), String, BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.parse_7z_list), String})
    precompile(Tuple{typeof(BinaryProvider.platform_key_abi), Base.SubString{String}})
    precompile(Tuple{typeof(BinaryProvider.libdir), BinaryProvider.Prefix})
    precompile(Tuple{typeof(BinaryProvider.detect_libgfortran_abi), String, BinaryProvider.Linux})
    precompile(Tuple{typeof(BinaryProvider.readuntil_many), Base.Pipe, Array{Char, 1}})
    precompile(Tuple{typeof(BinaryProvider.parse_dl_name_version), String, BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.parse_dl_name_version), String, BinaryProvider.Windows})
    precompile(Tuple{typeof(BinaryProvider.probe_symlink_creation), String})
    precompile(Tuple{typeof(BinaryProvider.locate), BinaryProvider.ExecutableProduct})
    precompile(Tuple{typeof(BinaryProvider.extract_platform_key), String})
    precompile(Tuple{typeof(BinaryProvider.triplet), BinaryProvider.Windows})
    precompile(Tuple{typeof(BinaryProvider.triplet), BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.triplet), BinaryProvider.Linux})
    precompile(Tuple{typeof(BinaryProvider.triplet), BinaryProvider.FreeBSD})
    precompile(Tuple{typeof(BinaryProvider.platforms_match), BinaryProvider.MacOS, BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.platform_key_abi), String})
    precompile(Tuple{typeof(BinaryProvider.platforms_match), BinaryProvider.Linux, BinaryProvider.Linux})
    precompile(Tuple{typeof(BinaryProvider.package), BinaryProvider.Prefix, String, Base.VersionNumber})
    precompile(Tuple{typeof(BinaryProvider.parse_tar_list), String})
    precompile(Tuple{typeof(BinaryProvider.choose_download), Base.Dict{BinaryProvider.Platform, String}, BinaryProvider.Windows})
    precompile(Tuple{typeof(BinaryProvider.choose_download), Base.Dict{BinaryProvider.Platform, Tuple{String, String}}, BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.choose_download), Base.Dict{BinaryProvider.Platform, String}, BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.__init__)})
    precompile(Tuple{typeof(BinaryProvider.choose_download), Base.Dict{BinaryProvider.Platform, String}, BinaryProvider.Linux})
    precompile(Tuple{typeof(BinaryProvider.platforms_match), BinaryProvider.Linux, BinaryProvider.UnknownPlatform})
    precompile(Tuple{typeof(BinaryProvider.platforms_match), BinaryProvider.Windows, BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.wordsize), BinaryProvider.Linux})
    precompile(Tuple{typeof(BinaryProvider.wordsize), BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.wordsize), BinaryProvider.FreeBSD})
    precompile(Tuple{typeof(BinaryProvider.parse_dl_name_version), String, BinaryProvider.Linux})
    precompile(Tuple{typeof(BinaryProvider.libdir), BinaryProvider.Prefix, BinaryProvider.Linux})
    precompile(Tuple{typeof(BinaryProvider.platforms_match), BinaryProvider.UnknownPlatform, BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.platform_key_abi)})
    precompile(Tuple{typeof(BinaryProvider.platforms_match), BinaryProvider.MacOS, BinaryProvider.Linux})
    precompile(Tuple{typeof(BinaryProvider.libdir), BinaryProvider.Prefix, BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.platforms_match), BinaryProvider.Linux, BinaryProvider.Windows})
    precompile(Tuple{typeof(BinaryProvider.list_tarball_files), String})
    precompile(Tuple{typeof(BinaryProvider.arch), BinaryProvider.Linux})
    precompile(Tuple{typeof(BinaryProvider.platforms_match), BinaryProvider.Linux, BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.arch), BinaryProvider.Windows})
    precompile(Tuple{typeof(BinaryProvider.platform_dlext)})
    precompile(Tuple{typeof(BinaryProvider.locate), BinaryProvider.LibraryProduct})
    precompile(Tuple{typeof(BinaryProvider.triplet), BinaryProvider.UnknownPlatform})
    precompile(Tuple{typeof(BinaryProvider.platform_dlext), BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.arch), BinaryProvider.FreeBSD})
    precompile(Tuple{typeof(BinaryProvider.collect_stdout), BinaryProvider.OutputCollector})
    precompile(Tuple{typeof(BinaryProvider.extract_name_version_platform_key), String})
    precompile(Tuple{typeof(BinaryProvider.platform_dlext), BinaryProvider.FreeBSD})
    precompile(Tuple{typeof(BinaryProvider.platform_dlext), BinaryProvider.Linux})
    precompile(Tuple{typeof(BinaryProvider.collect_stderr), BinaryProvider.OutputCollector})
    precompile(Tuple{typeof(BinaryProvider.platform_dlext), BinaryProvider.Windows})
    precompile(Tuple{typeof(BinaryProvider.platform_dlext), BinaryProvider.UnknownPlatform})
    precompile(Tuple{typeof(BinaryProvider.arch), BinaryProvider.UnknownPlatform})
    precompile(Tuple{typeof(BinaryProvider.valid_dl_path), String, BinaryProvider.Windows})
    precompile(Tuple{typeof(BinaryProvider.valid_dl_path), String, BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.valid_dl_path), String, BinaryProvider.Linux})
end


text:
using Knet
using ArgParse
using JLD
using Combinatorics

include("env.jl")
include("model.jl")
include("data.jl")
include("vocab.jl")

const CAPACITY = 20000
const EPS_INIT = 0.9
const EPS_FINAL = 0.005
const EPS_DECAY = 200

function main(args)
    s = ArgParseSettings()
    s.description = ""

    @add_arg_table s begin
        # load/save files
        ("--loadfile"; default=nothing)
        ("--savefile"; default=nothing)
        ("--task"; default="copy")
        ("--seed"; default=-1; arg_type=Int64; help="random seed")
        ("--optim"; default="Rmsprop()")
        ("--units"; default=200)
        ("--controller"; default="feedforward"; help="feedforward or lstm")
        ("--discount"; default=-1.; arg_type=Float64)
        ("--start"; default=6; arg_type=Int64)
        ("--end"; default=50; arg_type=Int64)
        ("--step"; default=4; arg_type=Int64)
        ("--batchsize"; default=20; arg_type=Int64)
        ("--nvalid"; default=60; arg_type=Int64)
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}":"Array{Float32}"))
        ("--period"; default=100; arg_type=Int64)
        ("--dist"; arg_type=String; default="randn"; help="[randn|xavier]")
        ("--supervised"; action=:store_true; help="if not, q-learning")
        ("--capacity"; default=CAPACITY; arg_type=Int64)
        ("--nsteps"; default=20; arg_type=Int64)
        ("--update"; default=5000; arg_type=Int64)
        ("--nepisodes"; default=5000; arg_type=Int64)
        ("--threshold"; default=0.98; arg_type=Float64)
        ("--noreward"; default=10; arg_type=Int64)
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)
    sr = o[:seed] > 0 ? srand(o[:seed]) : srand()
    o[:atype] = eval(parse(o[:atype]))
    data_generator = get_data_generator(o[:task])

    # init model, params etc.
    w = opts = s2i = i2s = nothing
    a2i, i2a = initvocab(get_actions(o[:task]))
    s2i, i2s = initvocab(get_symbols(o[:task]))
    if o[:loadfile] == nothing
        w = initweights(
            o[:atype],o[:units],length(s2i),length(a2i),o[:controller],o[:dist])
    else
        o[:task] = load(o[:loadfile], "task")
        w = load(o[:loadfile], "w")
        w = Dict(k=>convert(o[:atype],v) for (k,v) in w)
        # opts - not yet!
        # opts = load(o[:loadfile])
        # o[:start] = load(o[:loadfile], "complexity")
    end
    opts = initopts(w,o[:optim])

    # C => complexity
    # c => controller cell
    steps_done = 0
    for C = o[:start]:o[:step]:o[:end]
        # prepare validation data
        seqlen = div(C, complexities[o[:task]])
        data = map(xi->data_generator(seqlen), [1:o[:nvalid]...])
        valid = []; actions = []
        for (input,output,action) in data
            push!(valid, Game(input, output, o[:task]))
            push!(actions, action)
        end
        empty!(data)

        iter = 1
        lossval = 0
        while true
            # get examples for training
            trn = map(xi->data_generator(seqlen), 1:o[:nepisodes])

            # build environments
            games = []; actions = []
            for (input,output,action) in trn
                push!(games, Game(input, output, o[:task]))
                push!(actions, [action..., STOP_ACTION])
            end

            # train network while running episodes
            this_loss = train!(
                w,games,s2i,i2s,a2i,i2a,actions,opts,o[:supervised]; o=o)
            lossval = update_lossval(lossval, this_loss, iter)

            lossval = update_lossval(lossval,this_loss,iter)
            if iter % o[:period] == 0
                accuracy, average_reward = validate(w,valid,s2i,i2s,a2i,i2a;o=o)
                println("(iter:$iter,loss:$lossval,accuracy:$accuracy,reward:$average_reward)")
                if accuracy >= o[:threshold] && o[:savefile] != nothing
                    save(o[:savefile],
                         "w", Dict(k=>Array(v) for (k,v) in w),
                         # need something like above for opts
                         # "opts", opts,
                         "task", o[:task],
                         "complexity", C)
                end

                if accuracy >= o[:threshold]
                    println("$C converged in $iter iterations")
                    break
                end
            end

            iter += 1
        end # while true
    end # one complexity step
end

# train while running episodes
function train!(w,games,s2i,i2s,a2i,i2a,actions,opts,supervised; o=Dict())
    train = true
    run_episodes!(w,games,s2i,i2s,a2i,i2a,train,supervised;
                  o=o,opts=opts,actions=actions)
end

# validate games
function validate(w,games,s2i,i2s,a2i,i2a; o=Dict())
    train = supervised = false
    acc, reward = run_episodes!(w,games,s2i,i2s,a2i,i2a,train,supervised; o=o)
    return acc, reward
end

# this is skeleton for all methods
function run_episodes!(w,games,s2i,i2s,a2i,i2a,train,supervised;
                       o=Dict(), opts=Dict(), actions=[])

    # init state parameters
    atype = get(o, :atype, typeof(w[:wcont]))
    hidden = size(w[:wcont], 1)
    controller = get(o, :controller, "lstm")

    # unrolled
    nsteps = get(o, :nsteps, 20)

    # supervised learning
    if train && supervised
        isempty(actions) && error("actions must not be empty for supervision")

        # (1) prepare input data
        inputs, outputs, masks = make_data(games,s2i,a2i,actions)

        # (2) init controller states
        h, c = initstates(atype, hidden, length(games), controller)

        # (3) train network
        batchloss = sltrain!(w,h,c,inputs,outputs,masks,opts; o=o)

        # (4) return batchloss
        return batchloss
    end

    # needed by qwatkins
    histories = []
    for game in games; push!(histories, []); end

    # run episodes - for both validation and q-learning
    ncorrect = 0
    iter = 1

    cumulative_reward = 0
    for (i,game) in enumerate(games)
        episode_reward = 0
        h, c = initstates(atype, hidden, 1, controller)
        while !game.is_done
            # (1) prepare input data
            input_symbol = read_symbol(game)
            input_action = game.prev_actions[end]
            input = zeros(Cuchar, length(s2i)+length(a2i), 1)
            input[s2i[input_symbol],1] = 1
            input[length(s2i)+a2i[input_action],1] = 1

            # (2) propage controller
            prev_h = prev_c = nothing
            if train && h != nothing && c != nothing
                prev_h = Array(h)
                prev_c = Array(c)
            end
            cout, h, c = propagate(w[:wcont],w[:bcont],convert(atype,input),h,c)

            # (3) take action
            eps = train ? 0.0 : 0.20
            next_action = first(take_action(w[:wact],w[:bact],cout; eps=eps))
            next_action = i2a[next_action]
            move_action, write_action = next_action

            # (4) predict symbol
            y = NO_SYMBOL
            if write_action == WRITE
                y = predict(w[:wsymb],w[:bsymb],cout)
                # symbol = mapslices(indmax, Array(symbol), 1)[1]
                y = indmax(Array(y))
                y = i2s[y]
            end

            move_timestep!(game, y, move_action)
            reward = get_reward(game)
            episode_reward += reward

            # (5) add transition to history if phase is RL train
            if train
                remaining_steps = get_remaining_steps(game)

                # new transition
                transition = Transition(
                    reward,
                    input_symbol,
                    input_action,
                    remaining_steps,
                    prev_h, prev_c,
                    next_action,
                    y,
                    game.is_done
                )

                add_history = true
                if game.is_done
                    add_history = game.gold_tape == game.output_tape
                end

                # push it do history
                add_history && push!(histories[i], transition)
            end
        end

        cumulative_reward += episode_reward
        if game.output_tape == game.gold_tape
            ncorrect += 1
        end
    end

    # reset games
    for (i,game) in enumerate(games)
        reset!(game)
    end

    # q-watkins training if phase is RL train
    if train && !supervised
        discount = get(o, :discount, -1)
        batchsize = get(o, :batchsize, length(games))

        # (1) prepare input data
        batches = make_batches(
            w,histories,s2i,a2i,discount,nsteps,batchsize)

        # (2) train network
        batchloss = rltrain!(w,batches,opts; o=o)

        # (3) return batchloss
        return batchloss
    end

    return ncorrect/length(games), cumulative_reward/length(games)
end

function sltrain!(w,h,c,inputs,outputs,masks,opt; o=Dict())
    dw = similar(w)
    for k in keys(w); dw[k] = similar(w[k]); fill!(dw[k], 0); end

    total = num_samples = 0
    maxlen = get(o, :maxlen, 50)
    for i = 1:maxlen:length(inputs)
        lower = i; upper = min(i+maxlen-1,length(inputs))
        x = inputs[lower:upper]
        y = outputs[lower:upper]
        m = masks[lower:upper]

        values = []
        gloss = slgrad(w,x,y,m,h,c; values=values)

        # track loss and training sample values
        total += values[1]; num_samples += values[2]

        # accumulate gradients
        for k in keys(dw); dw[k] += gloss[k]; end

        # unbox states
        h = AutoGrad.getval(h)
        c = AutoGrad.getval(c)
    end

    # clip with num_samples
    for k in keys(dw); dw[k] = dw[k]/num_samples; end
    # FIXME: try a different thing in here

    # update weights
    update!(w,dw,opt)

    return total/num_samples
end

function rltrain!(w,batches,opt; o=Dict())
    dw = similar(w)
    for k in keys(w); dw[k] = similar(w[k]); fill!(dw[k], 0); end
    atype = get(o, :atype, AutoGrad.getval(typeof(w[:wcont])))

    total = num_samples = 0
    mapconvert(args...) = map(i->convert(atype,i), args)
    asarray(x) = typeof(x) <: Array ? x : typeof(x)[x]
    asarrays(xs...) = map(x->asarray(x), xs)
    for (ts,xs,ys,as,ms,hs,cs) in batches
        ts,ys,as,ms = asarrays(ts,ys,as,ms)
        # @show ts,ys,as,ms
        ts,xs = mapconvert(ts,xs)
        if hs != nothing && cs != nothing
            hs,cs = mapconvert(hs,cs)
        end

        values = []
        gloss = rlgrad(w,ts,xs,ys,as,ms,hs,cs; values=values)
        total += values[1]; num_samples += values[2]
        for k in keys(gloss); dw[k] += gloss[k]; end
    end

    for k in keys(w); dw[k] = dw[k]/num_samples; end
    update!(w,dw,opt)
    return total/num_samples
end

function update_lossval(lossval,batchloss,iter,alpha=0.01)
    if iter < 100
        lossval = (iter-1)*lossval + batchloss
        lossval = lossval / iter
    else
        lossval = alpha * batchloss + (1-alpha) * lossval
    end
    return lossval
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ereday/algorithm-learning

max_issues_repo_name: ereday/algorithm-learning

max_forks_repo_name: ereday/algorithm-learning


text:
using PedCar
using Test
using POMDPs
using AutomotiveDrivingModels
using AutomotivePOMDPs
using ProgressMeter

@testset begin "PedCar basics"
    params = UrbanParams(nlanes_main=1,
                     crosswalk_pos =[VecSE2(6, 0., pi/2), VecSE2(-6, 0., pi/2), VecSE2(0., -5., 0.)],
                     crosswalk_length =  [14.0, 14., 14.0],
                     crosswalk_width = [4.0, 4.0, 3.1],
                     stop_line = 22.0)
    env = UrbanEnv(params=params)
    mdp = PedCarMDP(env=env, pos_res=4.0, vel_res=3., ped_birth=0.3, car_birth=0.3)
    state_space = states(mdp)
    @test length(state_space) == n_states(mdp)
end

include("test_convert.jl")


text:
#workspace()

using EAGOSmoothMcCormickGrad
using IntervalArithmetic
using StaticArrays
#using DataFrames
#using CSV
using Plots
using PyPlot

# Plots the differentiable McCormick relaxations and
# original function x*y over [0,200] by [0,400], respectively
EAGOSmoothMcCormickGrad.set_diff_relax(1)

xl = 100
xu = 200
yl = -400
yu = -100
intv = 10
nx = (xu-xl)/intv+1
ny = (yu-yl)/intv+1
Vals = zeros(Float64,nx*ny,9)
x_grid = zeros(Float64,nx,ny)
y_grid = zeros(Float64,nx,ny)
values = zeros(Float64,nx,ny)
ccval = zeros(Float64,nx,ny)
cvval = zeros(Float64,nx,ny)
ccgrad1 = zeros(Float64,nx,ny)
ccgrad2 = zeros(Float64,nx,ny)
cvgrad1 = zeros(Float64,nx,ny)
cvgrad2 = zeros(Float64,nx,ny)

count1 = 1
for i=yl:intv:yu
    count2 = 1
    for j=xl:intv:xu
        temp1 = SMCg{2,Float64}(Float64(j),Float64(j),seed_g(Float64,1,2),
                                 seed_g(Float64,1,2),Interval(xl,xu),
                                 false, 1.0*[Interval(xl,xu),Interval(yl,yu)],[(xu-xl)/2.0,(yu-yl)/2.0])
        temp2 = SMCg{2,Float64}(Float64(i),Float64(i),seed_g(Float64,2,2),
                                 seed_g(Float64,2,2),Interval(yl,yu),
                                 false, 1.0*[Interval(xl,xu),Interval(yl,yu)],[(xu-xl)/2.0,(yu-yl)/2.0])
        x_grid[count2,count1] = j # x
        y_grid[count2,count1] = i # y
        values[count2,count1] = (j)*(i) #x*y
        temp_MC = temp1*temp2 # relax x*y

        ccval[count2,count1] = temp_MC.cc
        cvval[count2,count1] = temp_MC.cv
        ccgrad1[count2,count1] = temp_MC.cc_grad[1]
        ccgrad2[count2,count1] = temp_MC.cc_grad[2]
        cvgrad1[count2,count1] = temp_MC.cv_grad[1]
        cvgrad2[count2,count1] = temp_MC.cv_grad[2]
        count2 += 1
    end
    count1 += 1
end

#plot(layout = 4)
surface(x_grid,y_grid,cvval, xlab = "x", ylab = "y", zlab = "z")
p1 = surface(x_grid,y_grid,cvgrad1, xlab = "x", ylab = "y", zlab = "z")
#p2 = surface(x_grid,y_grid,cvgrad2, xlab = "x", ylab = "y", zlab = "z")
#pout = plot(p1,p2,layout = 2)
#pout = plot(p1,layout = 1)
#show(pout)
gui()
#=
subplot(221)
surface(x_grid,y_grid,ccval)
subplot(222)
surface(x_grid,y_grid,ccval)
subplot(223)
surface(x_grid,y_grid,ccval)
subplot(224)
surface(x_grid,y_grid,ccval)
suptitle("2x2 Subplot")
=#
#gui()
#DF = convert(DataFrame,Vals)
#CSV.write("/home/mewilhel/Desktop/MultiplicationPlot3.csv", DF)

seed1 = seed_g(Float64,1,2)
seed2 = seed_g(Float64,2,2)
x1 = SMCg{2,Float64}(0.0,0.0,seed1,seed1,Interval(-200.0,200.0),false, [Interval(-200.0,200.0),Interval(0.0,400.0)],[0.0,200.0])
y1 = SMCg{2,Float64}(200.0,200.0,seed2,seed2,Interval(0.0,400.0),false, [Interval(-200.0,200.0),Interval(0.0,400.0)],[0.0,200.0])
z1 = x1*y1
println("z1.cc: $(z1.cc)")
println("z1.cv: $(z1.cv)")

x2 = SMCg{2,Float64}(170.0,170.0,seed1,seed1,Interval(100.0,240.0),false, [Interval(100.0,240.0),Interval(100.0,400.0)],[170.0,250.0])
y2 = SMCg{2,Float64}(250.0,250.0,seed2,seed2,Interval(100.0,400.0),false, [Interval(100.0,240.0),Interval(100.0,400.0)],[170.0,250.0])
z2 = x2*y2
println("z2.cc: $(z2.cc)")
println("z2.cv: $(z2.cv)")

x3 = SMCg{2,Float64}(-200.0,-200.0,seed1,seed1,Interval(-300.0,-100.0),false, [Interval(-300.0,-100.0),Interval(-400.0,-200.0)],[-200.0,-300.0])
y3 = SMCg{2,Float64}(-300.0,-300.0,seed2,seed2,Interval(-400.0,-200.0),false, [Interval(-300.0,-100.0),Interval(-400.0,-200.0)],[-200.0,-300.0])
z3 = x3*y3
println("z3.cc: $(z3.cc)")
println("z3.cv: $(z3.cv)")

x4 = SMCg{2,Float64}(150.0,150.0,seed1,seed1,Interval(100.0,200.0),false, [Interval(100.0,200.0),Interval(-500.0,-100.0)],[150.0,-300.0])
y4 = SMCg{2,Float64}(-250.0,-250.0,seed2,seed2,Interval(-500.0,-100.0),false, [Interval(100.0,200.0),Interval(-500.0,-100.0)],[150.0,-300.0])
z4 = x4*y4
println("z4.cc: $(z4.cc)")
println("z4.cv: $(z4.cv)")

x5 = SMCg{2,Float64}(-150.0,-150.0,seed1,seed1,Interval(-200.0,-100.0),false, [Interval(-200.0,-100.0),Interval(200.0,400.0)],[-150.0,300.0])
y5 = SMCg{2,Float64}(300.0,300.0,seed2,seed2,Interval(200.0,400.0),false, [Interval(-200.0,-100.0),Interval(200.0,400.0)],[-150.0,300.0])
z5 = x5*y5
println("z5.cc: $(z5.cc)")
println("z5.cv: $(z5.cv)")