text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    SynchronousExecution

Functor object for using synchronous execution in an L-shaped algorithm (assuming multiple Julia cores are available). Create by supplying a [`Synchronous`](@ref) object through `execution` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
struct SynchronousExecution{T <: AbstractFloat,
                            A <: AbstractVector,
                            F <: AbstractFeasibilityAlgorithm,
                            I <: AbstractIntegerAlgorithm} <: AbstractLShapedExecution
    subworkers::Vector{SubWorker{T,F,I}}
    decisions::Vector{DecisionChannel}
    subobjectives::A
    model_objectives::A
    metadata::MetaDataChannel
    remote_metadata::Vector{MetaDataChannel}
    cutqueue::CutQueue{T}

    function SynchronousExecution(structure::StageDecompositionStructure{2, 1, <:Tuple{DistributedScenarioProblems}},
                                  feasibility_strategy::AbstractFeasibilityStrategy,
                                  integer_strategy::AbstractIntegerStrategy,
                                  ::Type{T},
                                  ::Type{A}) where {T <: AbstractFloat,
                                                    A <: AbstractVector}
        F = worker_type(feasibility_strategy)
        I = worker_type(integer_strategy)
        execution =  new{T,A,F,I}(Vector{SubWorker{T,F,I}}(undef, nworkers()),
                                  scenarioproblems(structure).decisions,
                                  A(),
                                  A(),
                                  RemoteChannel(() -> MetaChannel()),
                                  Vector{MetaDataChannel}(undef, nworkers()),
                                  RemoteChannel(() -> Channel{QCut{T}}(4 * nworkers() * num_scenarios(structure))))
        # Start loading subproblems
        load_subproblems!(execution.subworkers,
                          scenarioproblems(structure, 2),
                          execution.decisions,
                          feasibility_strategy,
                          integer_strategy)
        return execution
    end
end

function finish_initilization!(lshaped::AbstractLShaped, execution::SynchronousExecution)
    append!(execution.subobjectives, fill(1e10, num_thetas(lshaped)))
    append!(execution.model_objectives, fill(-1e10, num_thetas(lshaped)))
    for w in workers()
        execution.remote_metadata[w-1] = RemoteChannel(() -> MetaChannel(), w)
    end
    return lshaped
end

function mutate_subproblems!(mutator::Function, execution::SynchronousExecution)
    mutate_subproblems!(mutator, execution.subworkers)
    return nothing
end

function restore_subproblems!(::AbstractLShaped, execution::SynchronousExecution)
    restore_subproblems!(execution.subworkers)
    return nothing
end

function resolve_subproblems!(lshaped::AbstractLShaped, execution::SynchronousExecution{T}) where T <: AbstractFloat
    # Update metadata
    for w in workers()
        put!(execution.remote_metadata[w-1], timestamp(lshaped), :gap, gap(lshaped))
    end
    @sync begin
        for (i,w) in enumerate(workers())
            worker_aggregator = remote_aggregator(lshaped.aggregation, scenarioproblems(lshaped.structure), w)
            @async remotecall_fetch(resolve_subproblems!,
                                    w,
                                    execution.subworkers[w-1],
                                    execution.decisions[w-1],
                                    lshaped.x,
                                    execution.cutqueue,
                                    worker_aggregator,
                                    timestamp(lshaped),
                                    execution.metadata,
                                    execution.remote_metadata[w-1])
        end
    end
    # Assume no cuts are added
    added = false
    # Collect incoming cuts
    while isready(execution.cutqueue)
        _, cut::SparseHyperPlane{T} = take!(execution.cutqueue)
        added |= add_cut!(lshaped, cut)
    end
    # Return current objective value and cut_added flag
    return current_objective_value(lshaped), added
end

# API
# ------------------------------------------------------------
function (execution::Synchronous)(structure::StageDecompositionStructure{2, 1, <:Tuple{DistributedScenarioProblems}},
                                  feasibility_strategy::AbstractFeasibilityStrategy,
                                  integer_strategy::AbstractIntegerStrategy,
                                  ::Type{T},
                                  ::Type{A}) where {T <: AbstractFloat,
                                                    A <: AbstractVector}
    return SynchronousExecution(structure,
                                feasibility_strategy,
                                integer_strategy,
                                T,
                                A)
end

function str(::Synchronous)
    return "Synchronous "
end


text:
using Cxx

include("llvmincludes.jl")

RequireCompleteType(d::pcpp"clang::Type") = ccall(:RequireCompleteType,Cint,(Ptr{Void},),d.ptr) > 0
function cxxsizeof(d::pcpp"clang::CXXRecordDecl")
    executionEngine = pcpp"llvm::ExecutionEngine"(ccall(:jl_get_llvm_ee,Ptr{Void},()))
    cgt = pcpp"clang::CodeGen::CodeGenTypes"(ccall(:clang_get_cgt,Ptr{Void},()))
    dl = @cxx executionEngine->getDataLayout()
    RequireCompleteType(@cxx d->getTypeForDecl())
    t = @cxx cgt->ConvertRecordDeclType(d)
    @assert @cxx t->isSized()
    div((@cxx dl->getTypeSizeInBits(t)),8)
end
@assert cxxsizeof(pcpp"clang::CXXRecordDecl"(lookup_name(["llvm","ExecutionEngine"]).ptr)) >= 152

code_llvmf(f,t) = pcpp"llvm::Function"(ccall(:jl_get_llvmf, Ptr{Void}, (Any,Any,Bool), f, t, false))
function code_graph(f,args)
    v = @cxx std::string()
    os = @cxx llvm::raw_string_ostream(v)
    graphf = code_llvmf(f,args)
    @cxx llvm::WriteGraph(os,graphf)
    @cxx os->flush()
    bytestring((@cxx v->data()), (@cxx v->length()))
end

gt = code_graph(factorize,(typeof(rand(4,4)),))

@assert sizeof(gt) > 0

cxx"""
void f() {}
"""
clangmod = pcpp"llvm::Module"(unsafe_load(cglobal(:clang_shadow_module,Ptr{Void})))
ptr = @cxx clangmod->getFunction(pointer("_Z1fv"))
@assert ptr != C_NULL

jns = cglobal((:julia_namespace,libcxxffi),Ptr{Void})
ns = createNamespace("julia")

# This is basically the manual expansion of the cxx_str macro
unsafe_store!(jns,ns.ptr)
ctx = toctx(pcpp"clang::Decl"(ns.ptr))
d = CreateVarDecl(ctx,"var1",cpptype(Int64))
AddDeclToDeclCtx(ctx,pcpp"clang::Decl"(d.ptr))
cxxparse("""
extern llvm::Module *clang_shadow_module;
extern llvm::LLVMContext &jl_LLVMContext;
uint64_t foo() {
    return __julia::var1;
}
""")
unsafe_store!(jns,C_NULL)
#GV = @cxx dyn_cast{llvm::GlobalVariable}(@cxx (@cxx clang_shadow_module)->getNamedValue(pointer("_ZN5julia4var1E")))
GV = pcpp"llvm::GlobalVariable"((@cxx (@cxx clang_shadow_module)->getNamedValue(pointer("_ZN5julia4var1E"))).ptr)
@assert GV != C_NULL
@cxx (@cxx GV->getType())->dump()
@cxx GV->setInitializer(@cxx llvm::ConstantInt::get((@cxx llvm::Type::getInt64Ty(*(@cxx &jl_LLVMContext))),uint64(0)))
@cxx GV->setConstant(true)
@assert (@cxx foo()) == uint64(0)

# LLDB test

addHeaderDir(joinpath(basepath,"deps/llvm-svn/tools/lldb/include"))

# Because LLDB includes private headers from public ones! For shame.
addHeaderDir(joinpath(basepath,"deps/llvm-svn/tools/lldb/include/lldb/Target"))
defineMacro("LLDB_DISABLE_PYTHON") # NO!
cxxinclude("lldb/Interpreter/CommandInterpreter.h")
cxxinclude("lldb/Interpreter/CommandReturnObject.h")
cxxinclude("string", isAngled = true)

initd() = @cxx lldb_private::Debugger::Initialize(cast(C_NULL,vcpp"lldb_private::Debugger::LoadPluginCallbackType"))
initd()

debugger() = @cxx lldb_private::Debugger::CreateInstance()
const dbg = @cxx (debugger())->get()
ci(dbg::pcpp"lldb_private::Debugger") = @cxx dbg->GetCommandInterpreter()
function hc(ci::rcpp"lldb_private::CommandInterpreter",cmd)
    cro = @cxx lldb_private::CommandReturnObject()
    @cxx ci->HandleCommand(pointer(cmd),0,cro)
    @cxx cro->GetOutputData()
end

stdout = pcpp"FILE"(ccall(:fdopen,Ptr{Void},(Int32,Ptr{Uint8}),1,"a"))

function setSTDOUT(dbg,stdout)
    @cxx dbg->SetOutputFileHandle(stdout,false)
end

function setSTDERR(dbg,stdout)
    @cxx dbg->SetErrorFileHandle(stdout,false)
end

setSTDOUT(dbg,stdout)
setSTDERR(dbg,stdout)

ext: jl


text:
@testset "Player" begin

    # Test Player
    T = Float64
    p = 3
    model = UnicycleGame(p=p)
    x0 = @SVector [1.0, 2.0, 3.0, 4.0]
    id = 2
    lane_id = 1
    Q = Diagonal(rand(SVector{model.ni[1],T}))
    R = Diagonal(rand(SVector{model.mi[1],T}))
    xf = rand(SVector{model.ni[1],T})
    uf = rand(SVector{model.mi[1],T})
    u_min = rand(SVector{model.mi[1],T})
    u_max = rand(SVector{model.mi[1],T})
    v_min = rand()
    v_max = rand()
    r_col = 1.0
    r_cost = 1.0
    μ = 1.0
    player = Player(x0, id, lane_id, Q, R, xf, uf, u_min, u_max, v_min, v_max, r_col, r_cost, μ)
    @test typeof(player) <: Player

    lane = Lane()
    player1 = Player(model, lane)
    @test typeof(player1) <: Player

end


text:
include("../module/CustomStencil.jl")
include("../misc/misc.jl")
include("../misc/cpu_stencils.jl")
## Star Stencil definition with radius = 4
coefs = [0.75;0.5]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1,1, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst1 = NewStencilInstance(st_def, m_step=5)
st_inst2 = NewStencilInstance(st_def, m_step=6)
st_inst3 = NewStencilInstance(st_def, m_step=7)
st_inst4 = NewStencilInstance(st_def, m_step=8)

## Input Data size Definition
radius = 4
nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

data = CreateData(dx,dy,dz)

t_steps = 56
## Compile Functions
gpu_out = ApplyStencil(st_inst1, data, t_steps)
gpu_out = ApplyStencil(st_inst2, data, t_steps)
gpu_out = ApplyStencil(st_inst3, data, t_steps)
gpu_out = ApplyStencil(st_inst4, data, t_steps)

CUDA.cuProfilerStart()
gpu_out = ApplyStencil(st_inst1, data, t_steps)
gpu_out = ApplyStencil(st_inst2, data, t_steps)
gpu_out = ApplyStencil(st_inst3, data, t_steps)
gpu_out = ApplyStencil(st_inst4, data, t_steps)

exit()


text:
@non_differentiable Core.print(::Any...)
@non_differentiable Core.println(::Any...)

@non_differentiable Core.show(::Any)
@non_differentiable Core.show(::Any, ::Any)

@non_differentiable Core.isdefined(::Any, ::Any)
@non_differentiable Core.:(<:)(::Any, ::Any)

@non_differentiable Core.apply_type(::Any, ::Any...)
@non_differentiable Core.typeof(::Any)

if isdefined(Core, :_typevar)
    @non_differentiable Core._typevar(::Any...)
end
@non_differentiable TypeVar(::Any...)
@non_differentiable UnionAll(::Any, ::Any)

frule((_, ẋ, _), ::typeof(typeassert), x, T) = (typeassert(x, T), ẋ)
function rrule(::typeof(typeassert), x, T)
    typeassert_pullback(Δ) = (NoTangent(), Δ, NoTangent())
    return typeassert(x, T), typeassert_pullback
end

frule((_, _, ȧ, ḃ), ::typeof(ifelse), c, a, b) = (ifelse(c, a, b), ifelse(c, ȧ, ḃ))
function rrule(::typeof(ifelse), c, a, b)
    ifelse_pullback(Δ) = (NoTangent(), NoTangent(), ifelse(c, Δ, ZeroTangent()), ifelse(c, ZeroTangent(), Δ))
    return ifelse(c, a, b), ifelse_pullback
end
# ensure type stability for numbers
function rrule(::typeof(ifelse), c, a::Number, b::Number)
    ifelse_pullback(Δ) = (NoTangent(), NoTangent(), ifelse(c, Δ, zero(Δ)), ifelse(c, zero(Δ), Δ))
    return ifelse(c, a, b), ifelse_pullback
end


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

@doc raw"""
    ObjsEnterpriseUser(;
        enterprise_id=nothing,
        enterprise_name=nothing,
        id=nothing,
        is_admin=nothing,
        is_owner=nothing,
        teams=nothing,
    )

    - enterprise_id::DefsEnterpriseId
    - enterprise_name::DefsEnterpriseName
    - id::DefsEnterpriseUserId
    - is_admin::Bool
    - is_owner::Bool
    - teams::Vector{DefsTeam}
"""
mutable struct ObjsEnterpriseUser <: SwaggerModel
    enterprise_id::Any # spec type: Union{ Nothing, DefsEnterpriseId } # spec name: enterprise_id
    enterprise_name::Any # spec type: Union{ Nothing, DefsEnterpriseName } # spec name: enterprise_name
    id::Any # spec type: Union{ Nothing, DefsEnterpriseUserId } # spec name: id
    is_admin::Any # spec type: Union{ Nothing, Bool } # spec name: is_admin
    is_owner::Any # spec type: Union{ Nothing, Bool } # spec name: is_owner
    teams::Any # spec type: Union{ Nothing, Vector{DefsTeam} } # spec name: teams

    function ObjsEnterpriseUser(;enterprise_id=nothing, enterprise_name=nothing, id=nothing, is_admin=nothing, is_owner=nothing, teams=nothing)
        o = new()
        validate_property(ObjsEnterpriseUser, Symbol("enterprise_id"), enterprise_id)
        setfield!(o, Symbol("enterprise_id"), enterprise_id)
        validate_property(ObjsEnterpriseUser, Symbol("enterprise_name"), enterprise_name)
        setfield!(o, Symbol("enterprise_name"), enterprise_name)
        validate_property(ObjsEnterpriseUser, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ObjsEnterpriseUser, Symbol("is_admin"), is_admin)
        setfield!(o, Symbol("is_admin"), is_admin)
        validate_property(ObjsEnterpriseUser, Symbol("is_owner"), is_owner)
        setfield!(o, Symbol("is_owner"), is_owner)
        validate_property(ObjsEnterpriseUser, Symbol("teams"), teams)
        setfield!(o, Symbol("teams"), teams)
        o
    end
end # type ObjsEnterpriseUser

const _property_map_ObjsEnterpriseUser = Dict{Symbol,Symbol}(Symbol("enterprise_id")=>Symbol("enterprise_id"), Symbol("enterprise_name")=>Symbol("enterprise_name"), Symbol("id")=>Symbol("id"), Symbol("is_admin")=>Symbol("is_admin"), Symbol("is_owner")=>Symbol("is_owner"), Symbol("teams")=>Symbol("teams"))
const _property_types_ObjsEnterpriseUser = Dict{Symbol,String}(Symbol("enterprise_id")=>"DefsEnterpriseId", Symbol("enterprise_name")=>"DefsEnterpriseName", Symbol("id")=>"DefsEnterpriseUserId", Symbol("is_admin")=>"Bool", Symbol("is_owner")=>"Bool", Symbol("teams")=>"Vector{DefsTeam}")
Base.propertynames(::Type{ ObjsEnterpriseUser }) = collect(keys(_property_map_ObjsEnterpriseUser))
Swagger.property_type(::Type{ ObjsEnterpriseUser }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsEnterpriseUser[name]))}
Swagger.field_name(::Type{ ObjsEnterpriseUser }, property_name::Symbol) =  _property_map_ObjsEnterpriseUser[property_name]

function check_required(o::ObjsEnterpriseUser)
    (getproperty(o, Symbol("enterprise_id")) === nothing) && (return false)
    (getproperty(o, Symbol("enterprise_name")) === nothing) && (return false)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    (getproperty(o, Symbol("is_admin")) === nothing) && (return false)
    (getproperty(o, Symbol("is_owner")) === nothing) && (return false)
    (getproperty(o, Symbol("teams")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsEnterpriseUser }, name::Symbol, val)
end


text:
"""
    simulate!(rng::AbstractRNG, m::LinearMixedModel{T}; β=m.β, σ=m.σ, θ=T[])
    simulate!(m::LinearMixedModel; β=m.β, σ=m.σ, θ=m.θ)

Overwrite the response (i.e. `m.trms[end]`) with a simulated response vector from model `m`.
"""
function simulate!(
    rng::AbstractRNG,
    m::LinearMixedModel{T};
    β = coef(m),
    σ = m.σ,
    θ = T[],
) where {T}
    length(β) == length(fixef(m)) ||
        length(β) == length(coef(m)) ||
            throw(ArgumentError("You must specify all (non-singular) βs"))

    β = convert(Vector{T},β)
    σ = T(σ)
    θ = convert(Vector{T},θ)
    isempty(θ) || setθ!(m, θ)

    if length(β) ≠ length(coef(m))
        padding = length(coef(m)) - length(β)
        for ii in 1:padding
            push!(β, -0.0)
        end
    end

    y = randn!(rng, response(m))      # initialize y to standard normal

    for trm in m.reterms              # add the unscaled random effects
        unscaledre!(rng, y, trm)
    end
                    # scale by σ and add fixed-effects contribution
    BLAS.gemv!('N', one(T), m.X, β, σ, y)
    m
end

function simulate!(m::LinearMixedModel{T}; β = coef(m), σ = m.σ, θ = T[]) where {T}
    simulate!(Random.GLOBAL_RNG, m, β = β, σ = σ, θ = θ)
end

"""
    unscaledre!(y::AbstractVector{T}, M::ReMat{T}, b) where {T}
    unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, M::ReMat{T}) where {T}

Add unscaled random effects defined by `M` and `b` to `y`.  When `rng` is present the `b`
vector is generated as `randn(rng, size(M, 2))`
"""
function unscaledre! end

function unscaledre!(y::AbstractVector{T}, A::ReMat{T,1}, b::AbstractVector{T}) where {T}
    m, n = size(A)
    length(y) == m && length(b) == n || throw(DimensionMismatch(""))
    z = A.z
    @inbounds for (i, r) in enumerate(A.refs)
        y[i] += b[r] * z[i]
    end
    y
end

unscaledre!(y::AbstractVector{T}, A::ReMat{T,1}, B::AbstractMatrix{T}) where {T} =
    unscaledre!(y, A, vec(B))

function unscaledre!(y::AbstractVector{T}, A::ReMat{T,S}, b::AbstractMatrix{T}) where {T,S}
    Z = A.z
    k, n = size(Z)
    l = nlevs(A)
    length(y) == n && size(b) == (k, l) || throw(DimensionMismatch(""))
    @inbounds for (i, ii) in enumerate(A.refs)
        for j = 1:k
            y[i] += Z[j, i] * b[j, ii]
        end
    end
    y
end

function unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, A::ReMat{T,S}) where {T,S}
    rng_nums = randn(rng, S, nlevs(A))

    unscaledre!(y, A, lmul!(A.λ, rng_nums))
end

function unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, A::ReMat{T,1}) where {T}
    rng_nums = randn(rng, 1, nlevs(A))

    unscaledre!(y, A, lmul!(first(A.λ), rng_nums))
end

unscaledre!(y::AbstractVector, A::ReMat) = unscaledre!(Random.GLOBAL_RNG, y, A)


text:
using Millboard
using Base.Test

board = [(+) (*); (-) (/)]
@test """
+===+===+===+
|   | 1 | 2 |
+===+===+===+
| 1 | + | * |
+---+---+---+
| 2 | - | / |
+---+---+---+""" == table(board) |> string


text:
using Distributions
using LinearAlgebra
using StatsBase

mutable struct CategoricalSmoothed{CT <: Real} <: Distributions.Distribution{Distributions.Multivariate,Distributions.Discrete}
    counts::Vector{CT};
    smooth::Float64;
    sum_counts::CT;

    CategoricalSmoothed(counts::Vector{IT}; smooth::Real=1.0, sum_counts::IT=sum(counts)) where IT <: Real = new{IT}(counts, Float64(smooth), sum_counts)
end

function pdf(dist::CategoricalSmoothed, x::Int; use_smoothing::Bool=true)::Float64
    cnt = dist.counts[x]
    if !use_smoothing
        return cnt / dist.sum_counts
    end

    # Can't simply add smoothing to each vector component because of sparsity
    return fmax(Float64(cnt), dist.smooth) / (dist.sum_counts + dist.smooth)
end

function maximize!(dist::CategoricalSmoothed, x::T where T <: AbstractArray{Int, 1}, confidences::T2 where T2 <: AbstractVector{Float64})
    count_array!(dist.counts, x, confidences)
    dist.sum_counts = sum(confidences)
    return dist
end

max_stars_repo_name: YinuoJin/Baysor

max_issues_repo_name: YinuoJin/Baysor

max_forks_repo_name: YinuoJin/Baysor


text:
#=
             12  6
              │  ▲
              ▼  │
            ┌────┴──┐
      1◄────┤       │◄────11
      2◄────┤       │◄────10
            │   M   │
      7────►│       ├────►5
      8────►│       ├────►4
            └─┬─────┘
              │  ▲
              ▼  │
              3  9

Index convention for position in index array, is codomain counterclockwise from most
west-northerly (west first) positon, then domain in same order.

                           ┌───┐
     v ^               1◄──┤ a │◄──4
     2 3                   └┬──┘
--------------              │ ▲
     4 1                    ▼ │
     v ^                    2 3

3 ->[ a']-> 2               4 1
                            │ ▲
      3                     ▼ │
      v                    ┌──┴┐
2 ->[ a']->1           3──►│ a'│──►2
                           └───┘
=#

function applyH1(
    AC::TensorMap{S,2,2}, FL::TensorMap{S,3,3}, FR::TensorMap{S,3,3}, M::TensorMap{S,6,6}
) where {S<:IndexSpace}
    @tensor HAC[l ld; rd r] :=
        FL[l 8 7; 2 1 a] *
        AC[a 12; 6 b] *
        M[1 2 ld 4 5 6; 7 8 rd 10 11 12] *
        FR[11 10 b; 5 4 r]
    return HAC
end
function applyH1(
    AC::TensorMap{S,2,1}, FL::TensorMap{S,1,2}, FR::TensorMap{S,2,1}, M::TensorMap{S,2,2}
) where {S<:IndexSpace}
    @tensor HAC[l s; r] := FL[l; 1 a] * AC[a 4; b] * M[1 s; 3 4] * FR[3 b; r]
    return HAC
end
function applyH0(
    C::TensorMap{S,1,1}, FL::TensorMap{S,3,3}, FR::TensorMap{S,3,3}
) where {S<:IndexSpace}
    @tensor H0[L; R] := FL[L 2 3; 4 5 a] * C[a; b] * FR[5 4 b; 3 2 R]
    return H0
end
function applyH0(
    C::TensorMap{S,1,1}, FL::TensorMap{S,1,2}, FR::TensorMap{S,2,1}
) where {S<:IndexSpace}
    @tensor H0[L; R] := FL[L; m a] * C[a; b] * FR[m b; R]
    return H0
end

function _randFL(AL::TensorMap{S,2,2}, M::TensorMap{S,6,6}) where {S<:IndexSpace}
    spM = hcat(codomain(M)..., domain(M)...)
    spA = hcat(codomain(AL)..., domain(AL)...)
    return TensorMap(randn, ComplexF64, spA[1] * spM[8] * spM[7], spM[2] * spM[1] * spA[1])
end
function _randFL(AL::TensorMap{S,2,1}, M::TensorMap{S,2,2}) where {S<:IndexSpace}
    spM = hcat(codomain(M)..., domain(M)...)
    spA = hcat(codomain(AL)..., domain(AL)...)
    return TensorMap(randn, ComplexF64, spA[1], spM[1] * spA[1])
end

function _randFR(AR::TensorMap{S,2,2}, M::TensorMap{S,6,6}) where {S<:IndexSpace}
    spM = hcat(codomain(M)..., domain(M)...)
    spA = hcat(codomain(AR)..., domain(AR)...)
    return TensorMap(
        randn, ComplexF64, spM[11] * spM[10] * spA[4], spM[5] * spM[4] * spA[1]
    )
end
function _randFR(AR::TensorMap{S,2,1}, M::TensorMap{S,2,2}) where {S<:IndexSpace}
    spM = hcat(codomain(M)..., domain(M)...)
    spA = hcat(codomain(AR)..., domain(AR)...)
    return TensorMap(randn, ComplexF64, spM[3] * spA[3], spA[3])
end

function leftenv(
    AL::TensorMap{S}, M::TensorMap{S}; FL::TensorMap{S}=_randFL(AL, M)
) where {S<:IndexSpace}
    λs, FLs, info = eigsolve(
        x -> _flsolve(AL, M, x), FL, 1, :LM; eager=true, ishermitian=false
    )
    return FLs[1], λs[1], info
end

function _flsolve(
    AL::TensorMap{S,2,2}, M::TensorMap{S,6,6}, FL::TensorMap{S,3,3}
) where {S<:IndexSpace}
    @tensor FLu[-1 -2 -3; -4 -5 -6] :=
        AL[u 12; 6 -6] *
        FL[d 8 7; 2 1 u] *
        M[1 2 3 -2 -3 6; 7 8 9 -4 -5 12] *
        (AL')[9 -1; d 3]
    return FLu
end
function _flsolve(
    AL::TensorMap{S,2,1}, M::TensorMap{S,2,2}, FL::TensorMap{S,1,2}
) where {S<:IndexSpace}
    @tensor FLu[-1; -2 -3] := AL[u 4; -3] * FL[d; 1 u] * M[1 2; -2 4] * (AL')[-1; d 2]
    return FLu
end

function rightenv(
    AR::TensorMap{S}, M::TensorMap{S}; FR::TensorMap{S}=_randFR(AR, M)
) where {S<:IndexSpace}
    λs, FRs, info = eigsolve(
        x -> _frsolve(AR, M, x), FR, 1, :LM; eager=true, ishermitian=false
    )
    return FRs[1], λs[1], info
end

function _frsolve(
    AR::TensorMap{S,2,2}, M::TensorMap{S,6,6}, FR::TensorMap{S,3,3}
) where {S<:IndexSpace}
    @tensor FRu[-1 -2 -3; -4 -5 -6] :=
        AR[-3 12; 6 u] *
        M[-1 -2 3 4 5 6; -4 -5 9 10 11 12] *
        FR[11 10 u; 5 4 d] *
        (AR')[9 d; -6 3]
    return FRu
end
function _frsolve(
    AR::TensorMap{S,2,1}, M::TensorMap{S,2,2}, FR::TensorMap{S,2,1}
) where {S<:IndexSpace}
    @tensor FRu[-1 -2; -3] := AR[-2 4; u] * M[-1 2; 3 4] * FR[3 u; d] * (AR')[d; -3 2]
    return FRu
end

@doc raw"""
    leftgauge(A::TensorMap, [C]; kwargs...)

Given a MPS tensor `A`, return a left-canonical MPS tensor `AL`, a gauge transform `C` and
a scalar factor `λ` such that ``λ AL^s C = C A^s``, where an initial guess for `C` can be
provided.

                   ┌────┐   ┌───┐          ┌───┐   ┌───┐
            λ * ◄──┤ AL │◄──┤ C │◄─── = ◄──┤ C │◄──┤ A │◄──
                   └─┬──┘   └─┬─┘          └─┬─┘   └─┬─┘
                     ▼        ▼              ▼       ▼
                     s        s              s       s

"""
function leftgauge(
    A,
    C=TensorMap(randn, ComplexF64, codomain(A)[1], codomain(A)[1]);
    tol=1e-12,
    maxiter=100,
)
    AL, C, R, λ = _leftgaugeinitial(A, C; tol=tol)
    numiter = 1
    while norm(C - R) > tol && numiter < maxiter
        λs, Cs, info = eigsolve(
            x -> _lgsolve(A, AL, x),
            R,
            1,
            :LM;
            ishermitian=false,
            tol=tol,
            eager=true,
            maxiter=1,
        )
        AL, C, R, λ = _leftgaugefactor(Cs[1], A)
        numiter += 1
    end
    C = R
    return AL, C, λ
end

function _leftgaugefactor(C, A)
    _, C = leftorth(C)
    CA = _CA(C, A)

    n1 = numout(CA)
    inds = Tuple(1:numind(CA))

    Q, R = leftorth(CA, inds[1:(end - 1)], inds[end:end])
    # Q, R = leftorth(CA, (1, 2, 3), (4,))

    AL = permute(Q, inds[1:n1], inds[(n1 + 1):end])
    # AL = permute(Q, (1, 2), (3, 4))

    λ = norm(R)
    R = R * 1 / λ
    return AL, C, R, λ
end

@doc raw"""

    _CA(C::TensorMap{S,1,1}, A::TensorMap{S,2,2}) where {S<:IndexSpace}

Apply a matrix `C` to the left bond of a MPS tensor `A`

    ┌─┐   ┌───┐
 ◄──│C│◄──┤ A │◄──
    └─┘   └─┬─┘
            ▼

"""
function _CA(C::TensorMap{S,1,1}, A::TensorMap{S,2,2}) where {S<:IndexSpace}
    return @tensor CA[1 2; 3 4] := C[1; a] * A[a 2; 3 4]
end
function _CA(C::TensorMap{S,1,1}, A::TensorMap{S,2,1}) where {S<:IndexSpace}
    return @tensor CA[1 2; 3] := C[1; a] * A[a 2; 3]
end

# Apply a matrix C to the right bond of a MPS Tensor
# assumes left bond is first index, and right bond is last index, with an arb number
# of physical indicies.
function _AC(A::TensorMap{S}, C::TensorMap{S,1,1}) where {S<:IndexSpace}
    N = numind(A)
    n = numout(A)
    pind = (N, Tuple(2:(N - 1))..., 1)
    CA = _CA(permute(C, (2,), (1,)), permute(A, pind[1:n], pind[(n + 1):N]))
    return permute(CA, pind[1:n], pind[(n + 1):N])
end

function _leftgaugeinitial(A, C; tol=1e-12)
    λs, Ls, info = eigsolve(
        x -> _lgsolveinit(A, x),
        C' * C,
        1,
        :LM;
        ishermitian=false,
        tol=tol,
        maxiter=1,
        eager=true,
    )

    L = Ls[1]
    L = Ls[1] + Ls[1]'
    L = L / tr(L)
    _, Σ, Vt = tsvd(L)
    l = sqrt(Σ) * Vt

    AL, l, R, λ = _leftgaugefactor(l, A)
    return AL, l, R, λ
end

function _lgsolveinit(A::TensorMap{S,2,2}, L::TensorMap{S,1,1}) where {S<:IndexSpace}
    @tensor Lu[d; u] := A[a 2; 3 u] * (A')[3 d; b 2] * L[b; a]
    return Lu
end
function _lgsolveinit(A::TensorMap{S,2,1}, L::TensorMap{S,1,1}) where {S<:IndexSpace}
    @tensor Lu[d; u] := A[a 2; u] * (A')[d; b 2] * L[b; a]
    return Lu
end

function _lgsolve(
    A::TensorMap{S,2,2}, AL::TensorMap{S,2,2}, X::TensorMap{S,1,1}
) where {S<:IndexSpace}
    @tensor Y[d; u] := A[a 2; 3 u] * (AL')[3 d; b 2] * X[b; a]
    return Y
end
function _lgsolve(
    A::TensorMap{S,2,1}, AL::TensorMap{S,2,1}, X::TensorMap{S,1,1}
) where {S<:IndexSpace}
    @tensor Y[d; u] := A[a 2; u] * (AL')[d; b 2] * X[b; a]
    return Y
end

function rightgauge(
    A;
    C=TensorMap(randn, ComplexF64, codomain(A)[1], codomain(A)[1]),
    tol=1e-12,
    maxiter=100,
)
    N = numind(A)
    n = numout(A)
    pind = (N, Tuple(2:(N - 1))..., 1)
    AR, C, λ = leftgauge(permute(A, pind[1:n], pind[(n + 1):N]), permute(C, (2,), (1,)))
    return permute(AR, pind[1:n], pind[(n + 1):N]), permute(C, (2,), (1,)), λ
end

function _vumpsstep(AL, C, AR, M, FL, FR)
    # @tensor AC[-1 -2; -3 -4] := AL[-1 -2; -3 4] * C[4, -4]
    AC = _AC(AL, C)
    μ1s, ACs, info1 = eigsolve(
        x -> applyH1(x, FL, FR, M), AC, 1, :LM; ishermitian=false, eager=true, maxiter=1
    )
    μ0s, Cs, info0 = eigsolve(
        x -> applyH0(x, FL, FR), C, 1, :LM; ishermitian=false, eager=true, maxiter=1
    )

    λ = real(μ1s[1] / μ0s[1])
    AC = ACs[1]
    C = Cs[1]

    N = numind(AC)
    n = numout(AC)

    QAC, RAC = leftorth(AC, Tuple(1:(N - 1)), (N,))
    QC, RC = leftorth(C)
    AL = permute(QAC * QC', Tuple(1:n), Tuple((n + 1):N))
    errL = norm(RAC - RC)

    LAC, QAC = rightorth(AC, (1,), Tuple(2:N))
    LC, QC = rightorth(C)
    AR = permute(QC' * QAC, Tuple(1:n), Tuple((1 + n):N))
    errR = norm(LAC - LC)

    return λ, AL, C, AR, errL, errR
end

function _vumpsupdate(AL, C, AR, M, FL, FR)
    λ, AL, C, AR, errL, errR = _vumpsstep(AL, C, AR, M, FL, FR)
    AL, C = leftgauge(AR, C) # Regauge MPS: not really required
    FL, λL = leftenv(AL, M; FL=FL)
    FR, λR = rightenv(AR, M; FR=FR)

    # Convergence measure: norm of the projection of the residual onto the tangent space
    # @tensor AC[1 2; 3 4] := AL[1 2; 3 a] * C[a; 4]
    AC = _AC(AL, C)
    MAC = applyH1(AC, FL, FR, M)
    MC = applyH0(C, FL, FR)
    # @tensor MAC[1 2; 3 4] -= λL * AL[1 2; 3 a] * MC[a; 4]
    MAC = MAC - _AC(AL, MC) * λL
    err = norm(MAC)

    λs = (λ, λL, λR)
    mps = (AL, C, AR)
    env = (FL, FR)
    info = (err, errL, errR)

    return λs, mps, env, info
end

function vumps(
    A::TensorMap{S}, M::TensorMap{S}; verbose=false, tol=1e-8, maxiter=100
) where {S<:IndexSpace}
    AL, C, _ = leftgauge(A)
    AR, C, _ = rightgauge(AL)

    FL, λL, _ = leftenv(AL, M)
    FR, λR, _ = rightenv(AR, M)

    verbose && @info "Starting point has λ ≈ $λL ≈ $λR"

    λs, mps, env, info = _vumpsupdate(AL, C, AR, M, FL, FR)
    λ = λs[1]
    λL = real(λs[2])
    λR = real(λs[3])
    err = info[1]
    numiter = 1
    verbose && @info "Step $numiter: λ ≈ $(λ) ≈ $(λL) ≈ $(λR), err ≈ $(err)"

    while err > tol && numiter < maxiter
        λs, mps, env, info = _vumpsupdate(mps..., M, env...)
        λ = λs[1]
        λL = real(λs[2])
        λR = real(λs[3])
        err = info[1]
        numiter += 1
        verbose && @info "Step $numiter: λ ≈ $(λ) ≈ $(λL) ≈ $(λR), err ≈ $(err)"
    end
    return λ, AL, C, AR, FL, FR
end

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


text:
# [18](https://projecteuler.net/problem=18) / [67](https://projecteuler.net/problem=67)

#In a triangle like this:
#
#      a
#     b c
#    d e f
#
#the best way to find the anwser is not to get the maximum from the summaries of
#all possible branches from top to bottom.
#
#There is such a recursive pattern
#
#a + max( b+max(b,c), c+max(e,f)
#
# tri.txt
#75
#95 64
#17 47 82
#18 35 87 10
#20 04 82 47 65
#19 01 23 75 03 34
#88 02 77 73 07 63 67
#99 65 04 28 06 16 70 92
#41 41 26 56 83 40 80 70 33
#41 48 72 33 47 32 37 16 94 29
#53 71 44 65 25 43 91 52 97 51 14
#70 11 33 28 77 73 17 78 39 68 17 57
#91 71 52 38 17 14 91 43 58 50 27 29 48
#63 66 04 68 89 53 67 30 73 16 69 87 40 31
#04 62 98 27 23 09 70 98 73 93 38 53 60 04 23

import Base.zero
zero(::SubString{String}) = 0 # Julia 0.5

ZeroString(::SubString{String}) = 0
ZeroString(x::Int64) = x

A = readdlm("tri.txt")
A = ZeroString.(A)

function myreduction(m)
	if size(m)[1] == 1
		return m[1,1]
	else
		mprime = m[1:(end-1), :]
		for k in 1:(size(m, 1)-1)
			mprime[size(mprime, 1), k] += max(
									 m[size(m, 1),k], m[size(m, 1),k+1])
		end
		return myreduction(mprime)
	end
end

println(myreduction(A))

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
"""

`unique(arr, dims...)`

Return unique elements of an array `arr` of type `LabeledArray`/`DictArray`/`Nullable AbstractArrayWrapper`.

##### Arguments

* `arr` : an array
* `dims...` : either an integer or, if an array is a DictArray or a LabeledArray, a list of integers. It specifies the directions along which to traverse. Any duplicate elements will be replaced by Nullable{T}(). If all components along some direction are missing, those components will be removed and the whole array size will shrink.
If `dims...` is missing, unique elements along the whole directions will be found. It is equivalent to `unique(arr, 1, 2, ..., ndims(arr))`.
Note that it compares each slice spanned by directions orthogonal to `dims...`.

##### Examples

```julia
julia> unique(nalift([1 2 3;3 4 1]))
2x2 MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)  Nullable(2)
 Nullable(3)  Nullable(4)

julia> unique(nalift([1 2 3;3 4 1]), 1)
2x3 MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)  Nullable(2)  Nullable(3)
 Nullable(3)  Nullable(4)  Nullable(1)

julia> unique(nalift([1 2 3;3 4 1]), 2)
2x3 MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)  Nullable(2)  Nullable(3)
 Nullable(3)  Nullable(4)  Nullable(1)

julia> unique(nalift([1 2 3;1 2 3;4 5 6]), 1)
2x3 MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)  Nullable(2)  Nullable(3)
 Nullable(4)  Nullable(5)  Nullable(6)

julia> t = darr(a=[1 2 1;1 5 1;1 2 1], b=[:a :b :a;:a :c :a;:a :b :a])
3 x 3 DictArray

a b |a b |a b
----+----+----
1 a |2 b |1 a
1 a |5 c |1 a
1 a |2 b |1 a

julia> unique(t, 1)
2 x 3 DictArray

a b |a b |a b
----+----+----
1 a |2 b |1 a
1 a |5 c |1 a

julia> unique(t, 2)
3 x 2 DictArray

a b |a b
----+----
1 a |2 b
1 a |5 c
1 a |2 b

julia> m = larr(a=[1 2 1;1 5 1;1 2 1], b=[:a :b :a;:a :c :a;:a :b :a], axis1=["X","Y","Z"])
3 x 3 LabeledArray

  |1   |2   |3
--+----+----+----
  |a b |a b |a b
--+----+----+----
X |1 a |2 b |1 a
Y |1 a |5 c |1 a
Z |1 a |2 b |1 a

julia> unique(m, 1)
2 x 3 LabeledArray

  |1   |2   |3
--+----+----+----
  |a b |a b |a b
--+----+----+----
X |1 a |2 b |1 a
Y |1 a |5 c |1 a

julia> unique(m, 2)
3 x 2 LabeledArray

  |1   |2
--+----+----
  |a b |a b
--+----+----
X |1 a |2 b
Y |1 a |5 c
Z |1 a |2 b
```

"""
Base.unique{T}(arr::AbstractArrayWrapper{Nullable{T}}) = unique(arr, 1:ndims(arr)...)
Base.unique(arr::DictArray, dim::Int, dims::Int...) = unique(arr, 1:ndims(arr)...)
Base.unique(arr::LabeledArray, dim::Int, dims::Int...) = unique(arr, 1:ndims(arr)...)
Base.unique{T}(arr::AbstractArrayWrapper{Nullable{T}}, dim::Int, dims::Int...) = unique_inner(arr, (dim, dims...))[2]
Base.unique(arr::DictArray, dim::Int, dims::Int...) = unique_inner(arr, (dim, dims...))[2]
Base.unique(arr::LabeledArray, dim::Int, dims::Int...) = begin
  dimall = (dim, dims...)
  (coords_toshow, newdata) = unique_inner(arr.data, dimall)
  newaxes = ntuple(length(arr.axes)) do d
    axis = arr.axes[d]
    if d in dimall
      coords = coords_toshow[findfirst(dimall, d)]
      if isa(axis, DefaultAxis)
        DefaultAxis(length(coords))
      else
        arr.axes[d][coords]
      end
    else
      arr.axes[d]
    end
  end
  LabeledArray(newdata, newaxes)
end

@generated unique_inner{N}(arr::DictArray, dims::NTuple{N,Int}) = quote
  # first rearrange the axes so that the calculatioin can be more convenient.
  # this might be later optimized.
  ndimsarr = ndims(arr)
  sizearr = size(arr)
  remainingdirs = filter(x->!(x in dims), 1:ndimsarr)
  remainingcoords = ntuple(d->Colon(), length(remainingdirs))
  newdirs = (remainingdirs...,dims...)
  permuted_arr = permutedims_if_necessary(arr, newdirs)
  elems_sofar = Set()
  axiselems_toshow = [falses(sizearr[dims[i]]) for i in 1:$N]
  result = similar(permuted_arr, eltype(permuted_arr))
  @nloops $N i d->1:sizearr[dims[d]] begin
    coords = @ntuple($N,i)
    full_coords = (remainingcoords...,coords...)
    oneslice = getindexvalue(permuted_arr, full_coords...)
    if oneslice in elems_sofar
      setna!(result, full_coords...)
    else
      push!(elems_sofar, oneslice)
      for (index,coord) in enumerate(coords)
        axiselems_toshow[index][coord] = true
      end
      result[full_coords...] = oneslice
    end
  end
  coords_toshow = [find(axiselems_toshow[i]) for i in 1:$N]
  (coords_toshow, ipermutedims_if_necessary(result[remainingcoords...,coords_toshow...], newdirs))
end

@generated unique_inner{T,N}(arr::AbstractArray{Nullable{T}}, dims::NTuple{N,Int}) = quote
  # first rearrange the axes so that the calculatioin can be more convenient.
  # this might be later optimized.
  ndimsarr = ndims(arr)
  sizearr = size(arr)
  remainingdirs = filter(x->!(x in dims), 1:ndimsarr)
  remainingcoords = ntuple(d->Colon(), length(remainingdirs))
  newdirs = (remainingdirs...,dims...)
  permuted_arr = permutedims_if_necessary(arr, newdirs)
  elems_sofar = Set()
  axiselems_toshow = [falses(sizearr[dims[i]]) for i in 1:$N]
  result = similar(permuted_arr, eltype(permuted_arr))
  @nloops $N i d->1:sizearr[dims[d]] begin
    coords = @ntuple($N,i)
    full_coords = (remainingcoords...,coords...)
    oneslice = permuted_arr[full_coords...]
    if oneslice in elems_sofar
      setna!(result, full_coords...)
    else
      push!(elems_sofar, oneslice)
      for (index,coord) in enumerate(coords)
        axiselems_toshow[index][coord] = true
      end
      result[full_coords...] = oneslice
    end
  end
  coords_toshow = [find(axiselems_toshow[i]) for i in 1:$N]
  (coords_toshow, ipermutedims_if_necessary(result[remainingcoords...,coords_toshow...], newdirs))
end


text:
# A hacky function to return an empty Int64[] for the first conditioning set.
@inline cond_ixs(j, r) = j == 1 ? Int64[] : collect(max(1,j-r):max(1,j-1))

# number of elements in the lower triangle of an n x n matrix.
ltrisz(n) = div(n*(n+1), 2)

# Update kernel matrix buffer, exploiting redundancy for symmetric case. Not
# necessarily faster unless the kernel is very expensive to evaluate, but not
# slower in any case in my experimentation.
#
# Note that this does _NOT_ use threads, since I am already assuming that the
# nll function itself will be using threads, and in my benchmarking putting
# threaded constructors here slows things down a bit and increases allocations.
function updatebuf!(buf, pts1, pts2, kfun::F, params; skipltri=false) where{F}
  (F <: MemoizedKernel) && (@assert hash(params) === kfun.phash "params for memoized kernel don't agree with provided params. This shouldn't happen and is a bug.")
  if pts1 == pts2 && skipltri
    for k in eachindex(pts2)
      ptk = pts2[k]
      @inbounds buf[k,k] = kfun(ptk, ptk, params)
      @inbounds for j in 1:(k-1)
        buf[j,k] = kfun(pts1[j], ptk, params)
      end
    end
  elseif pts1 == pts2 && !skipltri
    for k in eachindex(pts2)
      ptk = pts2[k]
      @inbounds buf[k,k] = kfun(ptk, ptk, params)
      @inbounds for j in 1:(k-1)
        buf[j,k] = kfun(pts1[j], ptk, params)
        buf[k,j] = kfun(pts1[j], ptk, params)
      end
    end
  else
    @inbounds for k in eachindex(pts2), j in eachindex(pts1)
      buf[j,k] = kfun(pts1[j], pts2[k], params)
    end
  end
  nothing
end

# This function works pretty differently: now we assume that the points have all
# been catted together and that the kernel function takes entirely scalar
# inputs. With this formatting, we can then actually use the SIMD tools of
# LoopVectorization.jl and get some serious speedup.
#
# Very grateful to Chris Elrod (@elrod on discourse, @chriselrod on Github) for
# the help in making this work.
@generated function updatebuf_avx!(buf, ::Val{D}, pts1, pts2,
                                   kfun, params; skipltri=false) where{D}
  quote
    if pts1 == pts2 && skipltri
      for _k in 0:div(length(pts2)-1,$D)
        @turbo for _j in 0:_k # @turbo
        #@inbounds for _j in 0:_k
          val = kfun($([:(pts1[_j*$D+$d]) for d in 1:D]...),
                     $([:(pts2[_k*$D+$d]) for d in 1:D]...),
                     params)
          buf[_j+1,_k+1] = val
        end
      end
    elseif pts1 == pts2 && !skipltri
      for _k in 0:div(length(pts2)-1,$D)
        @turbo for _j in 0:_k # @turbo
        #@inbounds for _j in 0:_k
          val = kfun($([:(pts1[_j*$D+$d]) for d in 1:D]...),
                     $([:(pts2[_k*$D+$d]) for d in 1:D]...),
                      params)
          buf[_j+1,_k+1] = val
          buf[_k+1,_j+1] = val
        end
      end
    else
      @turbo for  _k in 0:div(length(pts2)-1,$D),  _j in 0:div(length(pts1)-1,$D) # @turbo
      #@inbounds for  _k in 0:div(length(pts2)-1,$D),  _j in 0:div(length(pts1)-1,$D)
        val = kfun($([:(pts1[_j*$D+$d]) for d in 1:D]...),
                   $([:(pts2[_k*$D+$d]) for d in 1:D]...),
                    params)
        buf[_j+1,_k+1] = val
      end
    end
    nothing
  end
end

# Primarily kept because it's readable and helps to transition from math
# notation in a paper to code.
function sunsteinchunk_naive(T, n, cpbuf, solve, fccov,
                             pts_ixs::AbstractVector{Int64},
                             cnd_ixs::AbstractVector{Int64})
  # allocate the sparse matrix for b_j^T. Eventually, a possible
  # micro-optimization would be to avoid this entirely and just build the final
  # one from the nzindices. But for now let's just do this.
  bt = spzeros(T, length(pts_ixs), n)
  # fill in the identity matrix part, doing the solves along the way.
  # about 50% of the time.
  for (j, ptixj) in enumerate(pts_ixs)
    ej    = zeros(T, length(pts_ixs))
    ej[j] = one(T)
    ldiv!(fccov.U', ej)
    bt[:,ptixj] .= ej
  end
  # solve the linear system, overwriting input "solve", and update the
  # corresponding columns of bt:
  ldiv!(fccov.U', Adjoint(solve))
  for (j, rowj) in enumerate(eachrow(solve))
    bt[:,cnd_ixs[j]] .= -rowj
  end
  # return the conjugation.
  return bt'bt # almost _all_ the allocs.
end

function sparseIJvecs_ltri(ixs)
  len = ltrisz(length(ixs))
  (Iv, Jv) = (Vector{Int64}(undef, len), Vector{Int64}(undef, len))
  current  = 1
  for j in eachindex(ixs)
    view(Iv, current:(current+j-1)) .= ixs[j]
    view(Jv, current:(current+j-1)) .= view(ixs, 1:j)
    current += j
  end
  (Iv, Jv)
end

function sparseIJvecs_full(ixs)
  ilen = length(ixs)
  len  = ilen^2
  (Iv, Jv) = (Vector{Int64}(undef, len), Vector{Int64}(undef, len))
  # fill in the I vector. Slightly more involved this time.
  current = 1
  for j in eachindex(ixs)
    view(Iv, current:(current+ilen-1)) .= ixs[j]
    current += ilen
  end
  # fill in the J vector. Again, slightly more involved.
  current = 1
  for j in eachindex(ixs)
    view(Jv, current:(current+ilen-1)) .= ixs
    current += ilen
  end
  (Iv, Jv)
end

function update_IJ!(ixs, I, J)
  @inbounds for l in eachindex(I,J)
    I[l] = ixs[I[l]]
    J[l] = ixs[J[l]]
  end
  nothing
end

function prepare_columns(T, solve, fccov)
  sz  = size(solve, 1) + size(fccov,1)
  out = Array{T}(undef, size(fccov,1), sz)
  for j in 1:size(solve,1)
    view(out, :, j) .= -view(solve, j, :)
  end
  for (row_counter, j) in enumerate((size(solve, 1)+1):sz)
    colj = view(out, :, j)
    fill!(colj, zero(T))
    colj[row_counter] = one(T)
  end
  lastchunk = view(out, :, (size(solve, 1)+1):sz)
  ldiv!(fccov.U', lastchunk)
  out
end

# TODO (cg 2021/04/24 18:26): this could always be optimized better.
# I suppose I could get rid of the allocation in combined by making the cross
# buffer slightly larger than it needs to be and updating that object. One of
# the dimensions is right already. But it's hard to imagine that boosting
# performance.
function sunsteinchunk(T, n, solve, fccov, mulbuf,
                       pts_ixs::AbstractVector{Int64},
                       cnd_ixs::AbstractVector{Int64})::Tuple{Vector{Int64}, Vector{Int64}, Vector{T}}
  # Incredibly, this hcat code seems to be more efficient than prepare_columns
  # above. I'm really confused about how it has fewer allocations.
  ldiv!(fccov.U', Adjoint(solve))
  combined     = hcat(-permutedims(solve), inv(fccov.U'))
  combined_ixs = vcat(cnd_ixs, pts_ixs)
  # index set:
  (Iv, Jv) = sparseIJvecs_ltri(1:length(combined_ixs))
  # fill in the matrix entries. The extra allocation isn't ideal, but BLAS-3!
  # having experimented with dropping the ltri-only approach so that I can just
  # return Vv = vec(mulbuf), this approach oddly has slightly more allocations
  # but was significantly faster---by about a factor of 2.
  mul!(mulbuf, Adjoint(combined), combined)
  Vv  = Vector{T}(undef, ltrisz(length(combined_ixs)))
  @turbo for l in eachindex(Iv, Jv)
  #@inbounds for l in eachindex(Iv, Jv)
    j     = Iv[l]
    k     = Jv[l]
    Vv[l] = mulbuf[j,k]
  end
  update_IJ!(combined_ixs, Iv, Jv)
  return (Iv, Jv, Vv)
end

# Primarily kept because it's readable and helps to transition from math
# notation in a paper to code.
function sunsteinchunk1_naive(T, n, buf)
  out = spzeros(T, n, n)
  len = size(buf, 1)
  out[1:len, 1:len] .= inv(buf)
  out
end

# something that instead returns nzindices instead of constructing the matrix.
function sunsteinchunk1(T, n, buf)::Tuple{Vector{Int64}, Vector{Int64}, Vector{T}}
  ibuf    = inv(buf)
  Vv      = Vector{T}(undef, ltrisz(size(buf, 1)))
  counter = 1
  @inbounds for j in 1:size(buf, 1)
    @inbounds for k in 1:j
      Vv[counter] = ibuf[j,k]
      counter += 1
    end
  end
  (Iv, Jv) = sparseIJvecs_ltri(1:size(buf, 1))
  return (Iv, Jv, Vv)
end

function globalidxs(datavv)
  (out, start) = (Vector{UnitRange{Int64}}(undef, length(datavv)), 1)
  for (j, datvj) in enumerate(datavv)
    len = size(datvj,1)
    out[j] = start:(start+len-1)
    start += len
  end
  out
end

function checksorted(V::VecchiaConfig{D,F}) where{D,F}
  all(issorted, V.condix) || throw(error("This function requires that every conditioning vector be sorted."))
  nothing
end

# TODO (cg 2022/04/21 16:16): This is totally not good.
function vec_of_vecs_to_matrows(vv)
  Matrix(reduce(hcat, vv)')
end

# For debugging. This gives M = U*U'.
function rchol(M)
  tmp = cholesky(Symmetric(reverse(Matrix(M), dims=(1,2)))).L
  UpperTriangular(reverse(Matrix(tmp), dims=(1,2)))
end

# Again, for debugging. Don't use this.
irchol(M) = inv(cholesky(M).U)

function prepare_v_buf!(buf, v::Matrix, idxv)
  _ix = 1
  for ixs in idxv
    for ix in ixs
      view(buf, _ix, :) .= view(v, ix, :)
      _ix += 1
    end
  end
  view(buf, 1:(_ix-1), :)
end

function updateptsbuf!(ptbuf, ptvv, idxs)
  ix = 1
  for idx in idxs
    for pt in ptvv[idx]
      ptbuf[ix] = pt
      ix += 1
    end
  end
  view(ptbuf, 1:(ix-1))
end

# Not a clever function at all,
function rchol_nnz(U::RCholesky{T}) where{T}
  # diagonal elements:
  out = sum(U.idxs) do ix
    n = length(ix)
    div(n*(n+1), 2)
  end
  # off-diagonal elements:
  out += sum(enumerate(U.condix)) do (j,ix_c)
    isempty(ix_c) && return 0
    tmp = 0
    len = length(U.idxs[j])
    for ix in ix_c
      tmp += len*length(U.idxs[ix])
    end
    tmp
  end
  out
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 57f367d6-025a-11eb-182a-773094da4307
using Yao, YaoPlots #calling the Yao and YaoPlots package

# ╔═╡ 3c9cf124-0257-11eb-16ce-0d857dce310f
md"# Using Yao - The basics of quantum computing in Julia using Yao.jl"

# ╔═╡ 5d5644b2-0257-11eb-36e9-7f4803cb3558
md"At the current moment, we don't have have quantum computers. How do we make quantum circuits then? Well, two things we can do right now are, simulate a few qubits or use the qubits created by corporates like IBM and D'Wave. Using Yao, we can simulate the qubits, without having a quantum computer\(based on the known mathematical and physics rules\), although the support to run your circuits on Yao using real qubits is coming to Yao soon."

# ╔═╡ 0f7d5bdc-0258-11eb-05bb-074d1b38c551
md"We can make a cicuit in Yao using *chain* function. For parameters we define the number of qubits and the operations we've to perform on them. Lets say we want to pass two qubits through two X gates. We do this by chain(number of qubits, operations). To use the X gate, we use the put() parameter. Run the cell below to see what happens."

# ╔═╡ ee15c64e-025a-11eb-2158-bd9cbea5f262
md" **\#\#** or single **#** sign mean comments in julia. It means that anything written after # or ## won't be read as a part of the program, in the line you used them."

# ╔═╡ dcf63604-0258-11eb-1466-fb21170c2e1a
let
	circuit = chain(2, put(1=>X), put(2=>X)); #define a variable "circuit" and "put" an X gate on the first qubit, and then put an X gate on the second qubit
	plot(circuit) #plot function, which takes a circuit for a parameter and prints the circuit diagram.
end

# ╔═╡ 841aa488-0259-11eb-351a-c1f1770b3a51
md"Assume we have 5 qubits and we have to pass each through an X gate. We can use the *repeat\(\)* parameter to pass the given number of qubits through the same gate."

# ╔═╡ af1aa2b4-0259-11eb-29fe-19ba9462deb9
plot(chain(5, repeat(X,1:5))) #plot function takes a circuit, which repeats the X gate on the qubits 1:5 or from 1st qubit to 5th qubit

# ╔═╡ 1efc3566-025a-11eb-2aae-fd953f22edfb
md"What about the Y, Z and H gate? Figure that out yourself. The next cell is left for you to play around."

# ╔═╡ 443afbfa-025a-11eb-03c6-e74bb5344e36
let
	##circuit =
	##plot(circuit)
end

# ╔═╡ 686ec0a0-025b-11eb-339e-1158d9b25529
md"What about multiqubit gates? We can use the control gate in Yao using the control\( \) parameter."

# ╔═╡ 22f24af0-025c-11eb-23d3-b9945bc05a36
plot(chain(2, control(1, 2=>X))) #Which translates to if the state of the 1st qubit is |1>, perform X gate to the 2nd qubit or "put" the 2nd qubit through the X gate.

# ╔═╡ Cell order:
# ╟─3c9cf124-0257-11eb-16ce-0d857dce310f
# ╟─5d5644b2-0257-11eb-36e9-7f4803cb3558
# ╟─0f7d5bdc-0258-11eb-05bb-074d1b38c551
# ╟─ee15c64e-025a-11eb-2158-bd9cbea5f262
# ╠═57f367d6-025a-11eb-182a-773094da4307
# ╠═dcf63604-0258-11eb-1466-fb21170c2e1a
# ╟─841aa488-0259-11eb-351a-c1f1770b3a51
# ╠═af1aa2b4-0259-11eb-29fe-19ba9462deb9
# ╟─1efc3566-025a-11eb-2aae-fd953f22edfb
# ╠═443afbfa-025a-11eb-03c6-e74bb5344e36
# ╟─686ec0a0-025b-11eb-339e-1158d9b25529
# ╠═22f24af0-025c-11eb-23d3-b9945bc05a36

max_stars_repo_name: logankilpatrick/tutorials

max_issues_repo_name: logankilpatrick/tutorials

max_forks_repo_name: logankilpatrick/tutorials


text:
struct StaticIndividualState <: AbstractArray{StaticPokemon, 1}
    mon1::StaticPokemon
    mon2::StaticPokemon
end

Base.size(s::StaticIndividualState) = (2,)
Base.IndexStyle(::Type{<:StaticIndividualState}) = IndexLinear()
Base.getindex(s::StaticIndividualState, i::UInt8) = i == 0x01 ? s.mon1 : s.mon2

StaticIndividualState(teams::Array{Int64}; league = "great", cup = "all") =
    StaticIndividualState(
        StaticPokemon(teams[1], league = league, cup = cup),
        StaticPokemon(teams[2], league = league, cup = cup)
    )

StaticIndividualState(teams::Array{String}; league = "great", cup = "all") =
    StaticIndividualState(
        StaticPokemon(teams[1], league = league, cup = cup),
        StaticPokemon(teams[2], league = league, cup = cup)
    )

struct DynamicIndividualState <: AbstractArray{DynamicPokemon, 1}
    mon1::DynamicPokemon
    mon2::DynamicPokemon
    data::UInt32
end

Base.size(d::DynamicIndividualState) = (2,)
Base.IndexStyle(::Type{<:DynamicIndividualState}) = IndexLinear()
Base.getindex(d::DynamicIndividualState, i::UInt8) = i == 0x01 ? d.mon1 : d.mon2

has_shield(state::DynamicIndividualState, agent::UInt8) = agent == 0x01 ?
  !iszero(state.data % UInt32(3)) :
  !iszero((state.data ÷ UInt32(3)) % UInt32(3))

get_fast_moves_pending(state::DynamicIndividualState) =
    (state.data ÷ UInt32(9)) % UInt32(7), (state.data ÷ UInt32(63)) % UInt32(7)

get_cmp(state::DynamicIndividualState) = (state.data ÷ UInt32(441)) % UInt32(5)

get_chance(state::DynamicIndividualState) = (state.data ÷ UInt32(2205)) %
    UInt32(6)

get_buffs(data::UInt32, attacker::Int64) = attacker == 1 ?
    ((data ÷ UInt32(13230))   % 9, (data ÷ UInt32(119070))  % 9) :
    ((data ÷ UInt32(1071630)) % 9, (data ÷ UInt32(9644670)) % 9)

DynamicIndividualState(state::StaticIndividualState; shields::Int8 = Int8(2)) =
    DynamicIndividualState(
        DynamicPokemon(state[0x01]),
        DynamicPokemon(state[0x02]),
        0x04*shields + UInt32(43394400)
    )


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "liblsl"
version = v"1.13.0-b13"

# Collection of sources required to build liblsl
sources = [
    "https://github.com/sccn/liblsl/archive/1.13.0-b13.tar.gz" =>
    "b7a1050cf4705d3d9917f7570e27e0c90c189634ecfc76906ecbfed45435b565",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd $WORKSPACE/srcdir
cd liblsl-1.13.0-b13/
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_BUILD_TYPE=Release -DLSL_UNIXFOLDERS=1 -DLSL_LSLBOOST_PATH="lslboost" ../
cmake --build . --config Release --target lsl

if [[ ${target} == *-linux-* ]]; then
    cp liblsl*.so.1.13.0 $prefix/liblsl.so
fi

if [[ ${target} == *-freebsd* ]]; then
    cp liblsl*.so.1.13.0 $prefix/liblsl.so
fi

if [[ ${target} == *-apple-* ]]; then
    cp liblsl64.dylib $prefix/liblsl.dylib
fi

if [[ ${target} == *-mingw* ]]; then
    cp liblsl*.dll $prefix/liblsl.dll
fi

exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "liblsl", :liblsl)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
struct Entity
    name::String
    address::String
end


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

mutable struct IssueContext <: SwaggerModel

    function IssueContext(;)
        o = new()
        o
    end
end # type IssueContext

const _property_map_IssueContext = Dict{Symbol,Symbol}()
const _property_types_IssueContext = Dict{Symbol,String}()
Base.propertynames(::Type{ IssueContext }) = collect(keys(_property_map_IssueContext))
Swagger.property_type(::Type{ IssueContext }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IssueContext[name]))}
Swagger.field_name(::Type{ IssueContext }, property_name::Symbol) =  _property_map_IssueContext[property_name]

function check_required(o::IssueContext)
    true
end

function validate_property(::Type{ IssueContext }, name::Symbol, val)
end


text:
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

include("constants.jl")

function mreplace(x)
    return replace(x, missing=>NaN)
end

function integrate(x, dydx)
    y = copy(dydx) * 0.0

    for i = 2:length(x)
        y[i] = y[i-1] + (dydx[i-1] + dydx[i]) * (x[i] - x[i-1]) / 2.0
    end

    return y
end

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file"
            help = "Atm data file."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["domain-file"], "r") do ds
    global mask = ds["mask"][:] |> mreplace
    global lat = ds["yc"][1, :] |> mreplace

    global lat_weight = Re * cos.(deg2rad.(lat)) * 2π
    global y          = Re * deg2rad.(lat)

    mask[mask.!=0] .= 1
end

Dataset(parsed["data-file"], "r") do ds

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12
    rng = (:,:,beg_t:end_t)

    global FFLX_TOA = ( - mean(ds["FSNT"][rng...], dims=(1, ))  + mean(ds["FLNT"][rng...],  dims=(1, )) )[1, :, :]
    global FFLX_SFC = ( - mean(ds["FSNS"][rng...], dims=(1, ))  + mean(ds["FLNS"][rng...],  dims=(1, )) )[1, :, :]
    global HFLX_SFC = (   mean(ds["SHFLX"][rng...], dims=(1, )) + mean(ds["LHFLX"][rng...], dims=(1, )) )[1, :, :]

    global (Ny, Nt) = size(FFLX_TOA)

end

IAET       = zeros(Float64, Ny, Nt)
EFLX_CONV  = zeros(Float64, Ny, Nt)

for t = 1:Nt

    FFLX_TOA[:, t] .*= lat_weight
    FFLX_SFC[:, t] .*= lat_weight
    HFLX_SFC[:, t] .*= lat_weight
    EFLX_CONV[:, t] = - ( FFLX_TOA[:, t] - FFLX_SFC[:, t] - HFLX_SFC[:, t] )
    IAET[:, t] = integrate(y, EFLX_CONV[:, t])

end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "Ny", Ny)
    defDim(ds, "time", Inf)

    for (varname, vardata, vardim, attrib) in [
        ("FFLX_TOA",  FFLX_TOA,  ("Ny", "time"), Dict()),
        ("FFLX_SFC",  FFLX_SFC,  ("Ny", "time"), Dict()),
        ("HFLX_SFC",  HFLX_SFC,  ("Ny", "time"), Dict()),
        ("EFLX_CONV", EFLX_CONV, ("Ny", "time"), Dict()),
        ("IAET",      IAET,      ("Ny", "time"), Dict()),
        ("FFLX_TOA_mean",  mean(FFLX_TOA, dims=(2,))[:, 1],  ("Ny", ), Dict()),
        ("FFLX_SFC_mean",  mean(FFLX_SFC, dims=(2,))[:, 1],  ("Ny", ), Dict()),
        ("HFLX_SFC_mean",  mean(HFLX_SFC, dims=(2,))[:, 1],  ("Ny", ), Dict()),
        ("EFLX_CONV_mean", mean(EFLX_CONV, dims=(2,))[:, 1], ("Ny", ), Dict()),
        ("IAET_mean",      mean(IAET, dims=(2,))[:, 1],      ("Ny", ), Dict()),
    ]

        println("Doing var: ", varname)

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
function value_iter(prop,reward,gamma,theta)

    a = size(prop,1)
    old = zeros(a,1,1)
    new = zeros(a)
    limit = 2*theta

    while limit > theta
        oyy = reward .+ gamma * permutedims(old,[2,1,3])
        #oyy = reward #.+ gamma * old
        apu = prop.*(oyy)
        #print("apu: \n", apu, "\n")
        summia = sum(apu, dims=2)
        #print(size(summia),"\n")
        #print("summia: ", summia, "\n")
        #new = reshape(maximum(summia,dims=3),s,1)
        new = maximum(summia,dims=3)
        #print("new: ", new, "\n")
        limit = maximum(abs.(old .- new), dims=1)
        limit = limit[1]
        #print("old: ", old, "\n")
        old = new
        print("Tolr: ", limit, "\n")
        #print("old: ", old, "\n")
        print("-----\n")
    end
    return old
end

function policya(prop,reward,gamma,old)
    summia = sum(prop.*(reward .+ gamma * permutedims(old,[2,1,3])), dims = 2)
    #apustus = reshape(summia,3,2)
    #print("\n summia polixcy \n", summia, "\n")
    #print("\n summia polixcy \n", size(summia), "\n")
    #print("\n apustus \n", apustus, "\n")
    #print("\n summia polixcy \n", size(apustus), "\n")
    apu_2 = Array{Any,1}(undef,size(summia,1))
    for i = 1:size(summia,1)
        helpp = summia[i,:,:]
        apu, apu_2[i] = findmax(helpp,dims=2)
        #print(summia[i,:,:], "\n")
        #print(apu_2, "\n")
    end
    #apu, apu_2 = findmax(summia,dims=3)
    return apu_2, summia
end

#i,j,k   state i -> j, action k

prop = zeros(6,6,3)
p_0_0 = 0.5
p_1_0 = 0.5

p_0_1 = 0.5
p_1_1 = 0.25
p_2_1 = 0.25

p_0_2 = 0.5
p_1_2 = 0.5

p_0_3 = 0.5
p_1_3 = 0.25
p_2_3 = 0.25

p_0_4 = 0.5
p_1_4 = 0.5

#state i -> j, action k
#S_0
prop[1,1,1]=1

prop[1,1,2]=p_0_0
prop[1,2,2]=p_1_0

#S_1
prop[2,2,1]=1

prop[2,2,2]=p_0_1
prop[2,5,2]=p_1_1
prop[2,3,2]=p_2_1

prop[2,2,3]=p_0_2
prop[2,3,3]=p_1_2

#S_2
prop[3,3,1]=1

prop[3,3,2]=p_0_3
prop[3,6,2]=p_1_3
prop[3,4,2]=p_2_3

prop[3,3,3]=p_0_4
prop[3,4,3]=p_1_4

# state i -> j, action k
reward = zeros(6,6,3)
reward[:,1,:] .= 1
reward[:,2,:] .= 1
reward[:,3,:] .= 1
reward[:,4,:] .= 100

tulos = value_iter(prop,reward,0.3,0.0001)
print("Values: ", tulos, "\n")

policy, summia = policya(prop,reward,0.3,tulos)

print("Best policy: ", "\n")
for i = 1:size(policy,1)
    print(policy[i][1][2], "\n")
end


text:
module CombineLeeching
    export CombineLeeching, Leeching, Classname, Playername, Gender, Race, Leecher1, Leechstatus
    struct LEECHING end
    struct CLASSNAME end
    struct PLAYERNAME end
    struct GENDER end
    struct LEECHER1 end
    struct LEECHSTATUS end
    struct RACE end
        const Leeching = LEECHING()
        const Classname = CLASSNAME()
        const Playername = PLAYERNAME()
        const Leecher1 = LEECHER1()
        const Leechstatus = LEECHSTATUS()
        const Gender = GENDER()
        const Race = RACE()
        CombineLeeching(Leeching,Playername,Gender,Race,Leechstatus,Leecher1)
end;

max_stars_repo_name: Physolia/-THE-Magic-Button

max_issues_repo_name: Physolia/-The-System-

max_forks_repo_name: Physolia/-The-System-


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BlockWeighter(sides)
    BlockWeighter(side₁, side₂, ...)

A weighting method that assigns weights to points in spatial object
based on blocks of given `sides`. The number `n` of points inside a
block determines the weights `1/n` of these points.
"""
struct BlockWeighter{T,N} <: AbstractWeighter
  sides::SVector{N,T}
end

BlockWeighter(sides::NTuple{N,T}) where {N,T} =
  BlockWeighter{T,N}(sides)

BlockWeighter(sides::Vararg) = BlockWeighter(sides)

weight(object, weighter::BlockWeighter) =
  weight(domain(object), weighter)

function weight(object::AbstractDomain, weighter::BlockWeighter)
  p = partition(object, BlockPartitioner(weighter.sides))

  weights = Vector{Float64}(undef, nelms(object))
  for s in subsets(p)
    n = length(s)
    weights[s] .= 1/n
  end

  SpatialWeights(object, weights)
end


text:
module TargetedLearning

using Reexport

include("Common.jl")
include("LReg.jl")
include("Qmodels.jl")
include("Parameters.jl")
include("TMLEs.jl")
include("CTMLEs.jl")

@reexport using .CTMLEs, .TMLEs, .LReg
import .Parameters: Mean, ATE
export Mean, ATE
import .Common: name!
export name!

end


text:
module lolPackage

greet() = print("Hello World!")
greet2()=print("WHY!")

using Distributions, Expectations, ForwardDiff, LinearAlgebra

function foo(μ = 1.0, σ = 2.0)
    d = Normal(μ, σ)
    E = expectation(d)
    return E(x -> sin(x))
end

function Newtonsroot1(f,f_prime;xiv,tol=1e-7,maxiter=1000)
    x_old = xiv
    error = Inf
    iter = 1
    while (error > tol || error < -1*tol) && iter <= maxiter
        x_new = x_old-(f(x_old)/f_prime(x_old)) # use the passed in map
        error = x_new - x_old
        x_old = x_new
        iter = iter + 1
    end
     if iter < maxiter+1
        return (root = x_old, normdiff = error, iter = iter) # A named tuple
    else
        return nothing
    end
end
D(f)= x ->ForwardDiff.derivative(f,x)

Newtonsroot1(f;xiv,tol=1e-7,maxiter=1000)= Newtonsroot1(f,D(f);xiv=xiv,tol=tol,maxiter=maxiter)

export foo, Newtonsroot1

end

max_stars_repo_name: anubhavpcjha/lolPackage

max_issues_repo_name: anubhavpcjha/lolPackage

max_forks_repo_name: anubhavpcjha/lolPackage


text:
#####################################
#Basis###############################
#####################################

abstract AbstractBasis{S<:Single} <: Dirac

immutable Basis{S<:Single} <: AbstractBasis{S}
	states::Vector{S}
	statemap::Dict
	adjhash::Uint #hash statemap for quick isdual() functionality
end

statemapper{S<:State}(sv::Vector{S}) = (labeltype(S)=>Int)[label(sv[i])=>i for i=1:length(sv)]
hashbasis{S<:Single}(sv::Vector{S}, statemap::Dict) = Basis(sv, statemap, hash(statemap))
consbasis{S<:Single}(sv::Vector{S}) = hashbasis(sv, statemapper(sv))
consbasis(bsym::Symbol,labelvec::Array) = consbasis(svec(bsym,labelvec))

#ensures basis states share a concrete type
makebasis{b,T}(sv::Vector{Ket{b,T}}) = consbasis(sv)
makebasis{b,T}(sv::Vector{Bra{b,T}}) = consbasis(sv)
makebasis{S<:Single}(sv::Vector{S}) = error("cannot construct basis: all states must have same label type and same basis identifier")

basis{S<:Single}(states::Vector{S}) = makebasis(unique(states))
basis{S<:Single}(states::Array{S}) = basis(vec(states))

basis(labelvec::Array, bsym::Symbol) = basis(svec(bsym, labelvec))

basis{S<:Single}(states::Array{S}) = basis(vec(states))
basis(bsym::Symbol, labelvec::Array) = basis(svec(bsym, labelvec))

#####################################
#TensorBasis#########################
#####################################

immutable TensorBasis{S<:Single,B<:Basis} <: AbstractBasis{S}
	bases::Vector{B}
	states::Vector{Tensor{S}}
	statemap::Dict
end

consbasis{B<:Basis,T<:Tensor}(bv::Vector{B}, sv::Vector{T}) = TensorBasis(bv, sv, statemapper(sv))

separate(b::Basis) = [b]
separate(b::TensorBasis) = b.bases

kind{K<:Ket}(b::AbstractBasis{K}) = Ket
kind{B<:Bra}(b::AbstractBasis{B}) = Bra

sepstates(sv) = hcat(map(separate, sv)...)

function consbvec!(bases, seps)
	for i=1:size(seps, 1)
	 	bases[i] = basis([seps[i, :]...]) #forcing correct typing with vcat
	end
	bases
end

function consbvec!{B<:Basis}(bases, seps::Vector{B})
	for i=1:length(seps)
	 	bases[i] = seps[i]
	end
	bases
end

prepbvec{S<:Single}(sv::Vector{Tensor{S}}) = Array(Basis{S}, length(sv[1]))
prepbvec{S<:Tensor}(sv::Vector{S}) = Array(Basis, length(sv[1]))

maketensorbasis{S<:Tensor}(sv::Vector{S}) = consbasis(consbvec!(prepbvec(sv),sepstates(sv)), sv)

basis{S<:Tensor}(states::Vector{S}) = maketensorbasis(unique(states))
basis{S<:Tensor}(states::Array{S}) = basis(vec(states))

#improved type inferencing for statejoin;
#based on the idea that each column of sarr
#has its own type due to the behavior of
#crossjoin
sjointype(sarr::Array)=typejoin(map(eltype, sarr[1,:])...)

statejoin{S<:Single}(sarr::Array{S,2}) = Tensor{S}[tensor(sarr[i, :]) for i=1:size(sarr, 1)]
statejoin{S<:State}(sarr::Array{S,2}) = Tensor{sjointype(sarr)}[tensor(sarr[i, :]) for i=1:size(sarr, 1)]

statecross(v) = statejoin(reduce(crossjoin,v))

tensor() = error("tensor needs arguments of type T<:Basis or T<:State")
tensor(b::AbstractBasis) = b

tensor(b::AbstractBasis...) = consbasis([[separate(i) for i in b]...], statecross([i.states for i in b]))

for t=(:Ket,:Bra)
	@eval begin
	tensor{S1<:($t), S2<:($t)}(s::State{S1}, b::AbstractBasis{S2}) = tensor(basis(s), b)
	tensor{S1<:($t), S2<:($t)}(b::AbstractBasis{S1}, s::State{S2}) = tensor(b, basis(s))
	end
end

#####################################
#Misc Functions######################
#####################################
basis{S<:State}(s::S...) = basis(collect(s))
basis(s::State...) = basis(kind(s[1])[s...])
basis(arr::Array) = basis([arr...]) #vcat forces retyping

copy{B<:Basis}(b::B) = B(copy(b.states), copy(b.statemap), copy(b.adjhash))
copy{B<:TensorBasis}(b::B)= B(copy(b.bases), copy(b.states), copy(b.statemap))
for t=(:Ket,:Bra)
	@eval bsym{b,T}(bs::Basis{($t){b,T}}) = b
end
bsym(b::TensorBasis) = map(bsym, b.bases)
adjhash(b::Basis) = b.adjhash
adjhash(b::TensorBasis) = map(adjhash, b.bases)

for op=(:isequal, :(==))
	@eval ($op){S<:Single}(a::Basis{S}, b::Basis{S}) = adjhash(a)==adjhash(b)
	@eval ($op){b,T}(x::TensorBasis{Bra{b,T}}, y::TensorBasis{Bra{b,T}}) = adjhash(x)==adjhash(y)
	@eval ($op){b,T}(x::TensorBasis{Ket{b,T}}, y::TensorBasis{Ket{b,T}}) = adjhash(x)==adjhash(y)
	@eval ($op){S1<:Ket,S2<:Ket}(a::TensorBasis{S1}, b::TensorBasis{S2}) = adjhash(a)==adjhash(b) && samebasis(a,b)
	@eval ($op){S1<:Bra,S2<:Bra}(a::TensorBasis{S1}, b::TensorBasis{S2}) = adjhash(a)==adjhash(b) && samebasis(a,b)
end

for op=(:length, :endof)
	@eval ($op)(b::AbstractBasis) = ($op)(b.states)
end

dual{S<:Single}(t::Type{Basis{S}}) = Basis{dual(S)}
dual{S<:Single}(t::Type{TensorBasis{S}}) = TensorBasis{dual(S)}

ctranspose(b::Basis) = Basis([ctranspose(i) for i in b.states], b.statemap, b.adjhash)
ctranspose(b::TensorBasis) = TensorBasis([ctranspose(i) for i in b.bases], [ctranspose(i) for i in b.states], b.statemap)

isdual{b,T}(x::Basis{Ket{b,T}}, y::Basis{Bra{b,T}}) = adjhash(x)==adjhash(y)
isdual{b,T}(x::Basis{Bra{b,T}}, y::Basis{Ket{b,T}}) = isdual(y,x)
isdual{b,T}(x::TensorBasis{Ket{b,T}}, y::TensorBasis{Bra{b,T}}) = adjhash(x)==adjhash(y)
isdual{b,T}(x::TensorBasis{Bra{b,T}}, y::TensorBasis{Ket{b,T}}) = isdual(y,x)

isdual{S1<:Ket,S2<:Bra}(a::TensorBasis{S1}, b::TensorBasis{S2}) = adjhash(a)==adjhash(b) && samebasis(a,b)
isdual{S1<:Bra,S2<:Ket}(a::TensorBasis{S1}, b::TensorBasis{S2}) = isdual(b,a)

isdual(a::AbstractBasis,b::AbstractBasis)=false

getindex(b::AbstractBasis, x) = b.states[x]

for t=(:Ket,:Bra)
	@eval begin
	get{b,T}(bs::Basis{($t){b,T}}, s::($t){b,T}, notfound) = get(bs.statemap, label(s), notfound)
	get{b,T}(bs::Basis{($t){b,T}}, s::($t){b,T}) = bs.statemap[label(s)]
	get{b,T}(bs::TensorBasis{($t){b,T}}, s::Tensor{($t){b,T}}, notfound) = get(bs.statemap, label(s), notfound)
	get{b,T}(bs::TensorBasis{($t){b,T}}, s::Tensor{($t){b,T}}) = bs.statemap[label(s)]

	get{S1<:($t),S2<:($t)}(b::TensorBasis{S1}, s::Tensor{S2}, notfound) = samebasis(b,s) ? get(b.statemap, label(s), notfound) : notfound
	get{S1<:($t),S2<:($t)}(b::TensorBasis{S1}, s::Tensor{S2}) = samebasis(b,s) ? b.statemap[label(s)] : throw(KeyError(s))

	in{b,T}(s::($t){b,T}, bs::Basis{($t){b,T}})= get(bs,s,"FALSE")=="FALSE" ? false : true
	in{S1<:($t),S2<:($t)}(s::Tensor{S1}, b::TensorBasis{S2})= get(b,s,"FALSE")=="FALSE" ? false : true
	end
end

get(b::AbstractBasis, s::State, notfound) = notfound
get(b::AbstractBasis, s::State) = throw(KeyError(s))

getpos(b::AbstractBasis, args...) = get(b, ars...)

in(s::State, b::AbstractBasis)= false
in(s, b::AbstractBasis)= error("in(s,b::AbstractBasis) is defined for s of type State, not $(typeof(s))")

######################################
##Show Functions######################
######################################

reprlabel(b::Basis) = bsym(b)
function reprlabel(b::TensorBasis)
	labels = bsym(b)
	#terrible way to grow a string
	str = "$(labels[1])"
	for i=2:length(labels)
		str = "$str$otimes$(labels[i])"
	end
	return str
end

summary(io::IO, b::AbstractBasis)= print(io, "$(typeof(b))")

showcompact(io::IO, b::AbstractBasis) = (summary(io, b); show(io, b.states))

function show(io::IO, b::AbstractBasis)
	summary(io, b)
	println(", $(length(b)) states:")
	if length(b)>20
		for i=1:10
			println(io, b.states[i])
		end
		println(vdots)
		show(io, b.states[length(b)-10])
		for i=length(b)-9:length(b)
			print('\n')
			show(io, b.states[i])
		end
	else
		show(io, b.states[1])
		for i=2:length(b.states)
			print('\n')
			show(io, b.states[i])
		end
	end
end

######################################
##Function-Passing Functions##########
######################################

find(f::Function, b::Basis) = find(f, b.states)
filter(f::Function, b::Basis) = consbasis(filter(f, b.states))
filter(f::Function, b::TensorBasis) = maketensorbasis(filter(f, b.states))
map(f::Function, b::AbstractBasis) = basis(map(f, b.states))

######################################
##Joining/Separating Functions########
######################################

function stateappend(d::Dict,s::State)
	resmap = copy(d)
	resmap[label(s)] = length(d)+1
	return resmap
end

function bcat{S}(a::Basis{S}, b::Basis{S})
	consbasis(unique(vcat(a.states, b.states)))
end

function bcat{S}(b::Basis{S}, s::S)
	if in(s, b)
		return b
	else
		return hashbasis(vcat(b.states, s), stateappend(b.statemap,s))
	end
end

function bcat{S}(s::S, b::Basis{S})
	if in(s, b)
		return b
	else
		return consbasis(vcat(s,b.states))
	end
end

function bcat{S}(b::TensorBasis{S}, s::Tensor{S})
	if in(s, b)
		return b
	else
		@assert samebasis(b, s) "BasisMismatch"
		return TensorBasis(b.bases, vcat(b.states, s), stateappend(b.statemap,s))
	end
end

function bcat{S}(s::Tensor{S},b::TensorBasis{S})
	if in(s, b)
		return b
	else
		@assert samebasis(b, s) "BasisMismatch"
		return consbasis(b.bases, vcat(s,b.states))
	end
end

bcat{S}(a::TensorBasis{S}, b::TensorBasis{S}) = basis(vcat(a.states, b.states))

kron(a::AbstractBasis, b::AbstractBasis) = tensor(a,b)
kron(a::AbstractBasis, b::State) = tensor(a,b)
kron(a::State, b::AbstractBasis) = tensor(a,b)

setdiff(a::AbstractBasis,b::AbstractBasis) = setdiff(a.states, b.states)


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Unit test file for materials.jl

@testset "pigments" begin
    @testset "UniformPigment" begin
        color = RGB(1f0, 2f0, 3f0)
        pigment = UniformPigment(color=color)

        @test pigment(0f0, 0f0) ≈ color
        @test pigment(1f0, 0f0) ≈ color
        @test pigment(0f0, 1f0) ≈ color
        @test pigment(1f0, 1f0) ≈ color
    end

    @testset "ImagePigment" begin
        image = HdrImage(2, 2)
        image[1, 1] = RGB(1f0, 2f0, 3f0)
        image[2, 1] = RGB(2f0, 3f0, 1f0)
        image[1, 2] = RGB(2f0, 1f0, 3f0)
        image[2, 2] = RGB(3f0, 2f0, 1f0)

        pigment = ImagePigment(image)
        @test pigment(0f0, 0f0) ≈ RGB(1f0, 2f0, 3f0)
        @test pigment(1f0, 0f0) ≈ RGB(2f0, 3f0, 1f0)
        @test pigment(0f0, 1f0) ≈ RGB(2f0, 1f0, 3f0)
        @test pigment(1f0, 1f0) ≈ RGB(3f0, 2f0, 1f0)
    end

    @testset "CheckeredPigment" begin
        color1 = RGB(1f0, 2f0, 3f0)
        color2 = RGB(10f0, 20f0, 30f0)

        pigment = CheckeredPigment{2}(color_on=color1, color_off=color2)

        # With num_of_steps == 2, the pattern should be the following:
        #
        #              (0.5, 0)
        #   (0, 0) +------+------+ (1, 0)
        #          |      |      |
        #          | col1 | col2 |
        #          |      |      |
        # (0, 0.5) +------+------+ (1, 0.5)
        #          |      |      |
        #          | col2 | col1 |
        #          |      |      |
        #   (0, 1) +------+------+ (1, 1)
        #              (0.5, 1)
        @test pigment(0.25f0, 0.25f0) ≈ color1
        @test pigment(0.75f0, 0.25f0) ≈ color2
        @test pigment(0.25f0, 0.75f0) ≈ color2
        @test pigment(0.75f0, 0.75f0) ≈ color1
    end
end


text:
using PairAsPipe
using Test

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

using DataFrames, DataFramesMeta, Pipe

data = DataFrame(a = rand(1:8, 100), b = rand(1:8, 100), c = rand(100))

@pipe data |>
    groupby(_, :a) |>
    @orderby(_, :b)


text:
using Knet
using ArgParse
using JLD
using HDF5
using MAT
using JSON
using AutoGrad

include("lib/vocab.jl")
include("lib/base.jl")
include("lib/init.jl")
include("lib/model.jl")
include("lib/batch.jl")
include("lib/convnet.jl")
include("lib/train.jl")
include("lib/eval.jl")
include("lib/util.jl")
include("lib/data.jl")

function main(args)
    s = ArgParseSettings()
    s.description = string(
        "Show and Tell: A Neural Image Caption Generator",
        " Knet implementation by Ilker Kesen [ikesen16_at_ku.edu.tr], 2016.")

    @add_arg_table s begin
        # load/save files
        ("--images"; help="images JLD file")
        ("--captions"; help="captions zip file (shared by Karpathy)")
        ("--vocabfile"; help="vocabulary JLD file")
        ("--loadfile"; default=nothing; help="pretrained model file if any")
        ("--savefile"; default=nothing; help="model save file after train")
        ("--cnnfile"; help="pre-trained CNN MAT file")
        ("--extradata"; action=:store_true;
         help="use restval split for training")

        # model options
        ("--winit"; arg_type=Float32; default=Float32(0.08))
        ("--hidden"; arg_type=Int; default=512)
        ("--embed"; arg_type=Int; default=512)
        ("--convnet"; default="vgg16")
        ("--lastlayer"; default="relu7")
        ("--visual"; arg_type=Int; nargs='+'; default=[4096])
        ("--cnnmode"; arg_type=Int; default=1)

        # training options
        ("--nogpu"; action=:store_true)
        ("--epochs"; arg_type=Int; default=1)
        ("--batchsize"; arg_type=Int; default=200)
        ("--lr"; arg_type=Float32; default=Float32(0.001))
        ("--gclip"; arg_type=Float32; default=Float32(5.0))
        ("--seed"; arg_type=Int; default=-1; help="random seed")
        ("--gcheck"; arg_type=Int; default=0; help="gradient checking")
        ("--finetune"; action=:store_true; help="fine tune convnet")
        ("--optim"; default="Adam"; help="optimizer (Sgd|Adam|Adagrad)")
        ("--eps"; default=1e-6; help="epsilon for Adagrad optimizer")
        ("--adam"; action=:store_true; help="use adam optimizer")
        ("--decay"; arg_type=Float32; default=Float32(1.0); help="lr decay")
        ("--decayperiod"; arg_type=Int64; default=0)
        ("--fast"; action=:store_true; help="do not compute train loss")
        ("--saveperiod"; arg_type=Int; default=0)
        ("--newoptimizer"; action=:store_true)
        ("--evalmetric"; default="bleu")
        ("--beamsize"; arg_type=Int; default=1)
        ("--checkpoints"; arg_type=Int; default=1)
        ("--sortbylen"; action=:store_true)

        # dropout values
        ("--fc6drop"; arg_type=Float32; default=Float32(0.5))
        ("--fc7drop"; arg_type=Float32; default=Float32(0.5))
        ("--softdrop"; arg_type=Float32; default=Float32(0.5))
        ("--wembdrop"; arg_type=Float32; default=Float32(0.5))
        ("--vembdrop"; arg_type=Float32; default=Float32(0.5))
        ("--membdrop"; arg_type=Float32; default=Float32(0.5))
        ("--attdrop"; arg_type=Float32; default=Float32(0.5))
    end

    # parse args
    @printf("\nScript started. [%s]\n", now()); flush(STDOUT)
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); display(o); flush(STDOUT)

    # random seed
    s = o[:seed] > 0 ? srand(o[:seed]) : srand()
    @printf("\nrandom seed:\n")
    display(s.seed); println(); flush(STDOUT)

    # load vocabulary
    vocab = load(o[:vocabfile], "vocab")
    o[:vocabsize] = vocab.size
    println("Vocabulary loaded."); flush(STDOUT)

    # initialize state and weights
    o[:atype] = !o[:nogpu] ? KnetArray{Float32} : Array{Float32}
    prevscore = bestscore =
        o[:loadfile] == nothing ? Inf : load(o[:loadfile], "score")
    prevloss = bestloss =
        o[:loadfile] == nothing ? Inf : load(o[:loadfile], "lossval")
    w, srnn = get_weights(o)
    wcnn =  get_wcnn(o)
    if wcnn != nothing && o[:finetune]
        w["wcnn"] = wcnn
    end
    opts = get_opts(o, w)
    println("Model loaded."); flush(STDOUT)

    # get samples used during training process
    train, restval, val = get_entries(o[:captions], ["train", "restval", "val"])
    if o[:extradata]
        train = [train; restval]
    else
        val = [val; restval]
    end
    restval = 0

    # split samples into image/sentence pairs
    train = get_pairs(train)
    valid = get_pairs(val) # keep val for validation
    gc()
    const nsamples = length(train)
    const nbatches = div(nsamples, o[:batchsize])
    println("Data loaded."); flush(STDOUT)

    # gradient check
    # if o[:gcheck] > 0
    #     ids = shuffle([1:nsamples...])[1:o[:batchsize]]
    #     images, captions, masks = make_batch(o, train[ids], vocab)
    #     gradcheck(loss, w, copy(s), images, captions, masks;
    #               gcheck=o[:gcheck], verbose=true, atol=0.01)
    #     images
    #     gc()
    # end

    # checkpoints
    checkpoints = []

    # sort sequences
    if o[:sortbylen]
        sort!(train, by=i->length(i[2]), rev=true)
    end
    sort!(valid, by=i->length(i[2]), rev=true)
    offsets = collect(1:o[:batchsize]:nsamples+1)

    # training
    const saveperiod = o[:saveperiod] > 0 ? o[:saveperiod] : length(offsets)-1
    @printf("Training started (nsamples=%d, nbatches=%d, loss=%g, score=%g). [%s]\n",
            nsamples, nbatches, prevloss, prevscore, now()); flush(STDOUT)
    for epoch = 1:o[:epochs]
        t0 = now()

        # data split training
        losstrn = 0
        nwords  = 0
        if !o[:sortbylen]
            shuffle!(train)
            orders = [1:(length(offsets)-1)...]
        else
            orders = randperm(length(offsets)-1)
        end
        for (i,k) in enumerate(orders)
            iter = (epoch-1)*nbatches+i
            lower, upper = offsets[k:k+1]
            samples = train[lower:upper-1]
            images, x, y, batchsizes = make_batch(o, samples, vocab)
            this_loss, this_words = train!(
                w, srnn, images, x, y, batchsizes, opts, o)
            flush(STDOUT)
            images = 0; captions = 0; ans = 0; gc()
            losstrn += this_loss
            nwords  += this_words

            if iter % saveperiod == 0
                lossval = bulkloss(w,srnn,o,valid,vocab)
                @printf("\n(epoch/iter): %d/%d, loss: %g/%g [%s] ",
                        epoch, iter, losstrn/nwords, lossval, now())
                flush(STDOUT)
                # scores, bp, hlen, rlen =
                #     validate(w, val, vocab, o)
                # @printf("\nBLEU = %.1f/%.1f/%.1f/%.1f ",
                #         map(i->i*100,scores)...)
                # @printf("(BP=%g, ratio=%g, hyp_len=%d, ref_len=%d) [%s]\n",
                #         bp, hlen/rlen, hlen, rlen, now())
                # flush(STDOUT)
                # score = scores[end]
                score = lossval

                # learning rate decay
                decay!(o, opts, lossval, prevloss)
                prevscore = score
                prevloss  = lossval
                gc()

                # check and save best model
                score <= bestscore || continue

                path, ext = splitext(abspath(o[:savefile]))
                filename  = abspath(string(path, "-iter-", iter, ext))
                savemodel(o, w, opts, filename, score, lossval)
                @printf("Model saved to %s.\n", filename); flush(STDOUT)

                # keep track of checkpoints
                push!(checkpoints, (score, -lossval, filename))
                sort!(checkpoints, rev=true)
                if length(checkpoints) > o[:checkpoints]
                    _, _, worst = pop!(checkpoints)
                    rm(worst)
                end
                bestscore = checkpoints[end][1]
            end
        end # batches end

        t1 = now()
        elapsed = Int64(round(Float64(t1-t0)*0.001))
        @printf("\nepoch #%d finished. (time elapsed: %s)\n",
                epoch, pretty_time(elapsed))
        flush(STDOUT)
        gc()
    end # epoch end
end

function decay!(o, opts, lossval, prevloss)
    if o[:decay] < 1.0 && lossval > prevloss
        o[:lr] *= o[:decay]
        @printf("\nlr decay. new lr=%g\n", o[:lr]); flush(STDOUT)
        decay!(o[:lr], opts)
    end
end

function decay!(lr, opts::Dict)
    for k in keys(opts)
        decay!(lr, opts[k])
    end
end

function decay!(lr, opts::Array)
    for k in 1:length(opts)
        decay!(lr, opts[k])
    end
end

function decay!(lr, opt::Union{Knet.Adam,Knet.Sgd,Knet.Adagrad})
    opt.lr = lr
end

function savemodel(o, w, opts, filename, score, lossval)
    if o[:savefile] != nothing
        save(filename,
             "w", copy_weights(w),
             "opts", copy_opts(opts, w, true),
             "score", score,
             "lossval", lossval)
    end
end

function validate(w, data, vocab, o; metric=bleu, split="val")
    wcnn = o[:finetune] ? w[:wcnn] : nothing
    wdec = w
    hyp, ref = [], []
    s = initstate(o[:atype], o[:hidden], 1)

    hyp, ref = h5open(o[:images], "r") do f
        hyp, ref = [], []
        for entry in data
            image = read(f, entry["filename"])
            caption = generate(
                wdec, copy(s), image, vocab; beamsize=o[:beamsize])
            push!(hyp, caption)
            push!(ref, map(s->s["raw"], entry["sentences"]))
        end
        hyp, ref
    end

    return metric(hyp, ref)
end

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

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
module PandorasBoxMarioShell

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/shell" Shell(x::Integer, y::Integer, texture::String="koopa", lights::Bool=false, color::String="Green", colorSpeed::Number=0.8, direction::Integer=0, dangerous::Bool=false, grabbable::Bool=true)

const textures = Dict{String, String}(
    "Koopa" => "koopa",
    "Beetle" => "beetle",
    "Spiny" => "spiny",
    "Bowser Jr." => "bowserjr"
)

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

const placements = Ahorn.PlacementDict(
    "Shell (Pandora's Box)" => Ahorn.EntityPlacement(
        Shell
    )
)

Ahorn.editingOptions(entity::Shell) = Dict{String, Any}(
    "texture" => textures,
    "color" => colors,
    "direction" => Dict{String, Any}(
        "Left" => -1,
        "Still" => 0,
        "Right" => 1
    )
)

function getColor(color)
    if haskey(Ahorn.XNAColors.colors, color)
        return Ahorn.XNAColors.colors[color]

    else
        try
            return ((Ahorn.argb32ToRGBATuple(parse(Int, replace(color, "#" => ""), base=16))[1:3] ./ 255)..., 1.0)

        catch

        end
    end

    return (1.0, 1.0, 1.0, 1.0)
end

function getTextures(entity::Shell)
    texture = get(entity.data, "texture", "koopa")

    return "objects/pandorasBox/shells/$texture/shell_idle00.png", "objects/pandorasBox/shells/$texture/deco_idle00.png"
end

function Ahorn.selection(entity::Shell)
    x, y = Ahorn.position(entity)

    shellSprite, decoSprite = getTextures(entity)

    return Ahorn.coverRectangles([
        Ahorn.getSpriteRectangle(shellSprite, x, y),
        Ahorn.getSpriteRectangle(decoSprite, x, y),
    ])
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Shell, room::Maple.Room)
    shellSprite, decoSprite = getTextures(entity)

    colorRaw = entity.color
    color = getColor(colorRaw)

    Ahorn.drawSprite(ctx, shellSprite, 0, 0, tint=color)
    Ahorn.drawSprite(ctx, decoSprite, 0, 0)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
using ReverseDiff
using IterativeSolvers
using Tools
using Baseline
using Printf
using Statistics
using Random

struct NPGStrategy{dtype<:AbstractFloat} <: LearningStrategy
   vpg::Vector{dtype}
   npg::Vector{dtype}
   grad::Vector{dtype}

   fullFIM::Bool
   gradll::Matrix{dtype}
   FIM::Matrix{dtype}

   adv::Vector{dtype}
   ret::Vector{dtype}

   baseline::Baseline.AbstractBaseline

   norm_step_size::Float64
   gamma::Float64
   gae::Float64

   cg_iter::Int
   cg_reg::Float64
   cg_tol::Float64

   obsmat::Matrix{dtype}
   ctrlmat::Matrix{dtype}

   # For rollout
   f::FunctionSet
   policyaction::Function
   policymean::Function

   function NPGStrategy{dtype}(model,
                               baseline,
                               fullFIM,
                               norm_step_size,
                               gamma, gae,
                               cg_iter, cg_reg, cg_tol,
                               funcset) where dtype<:AbstractFloat
      #dtype = eltype(model.policy.theta) #Float64
      T, K = model.samples.T, model.samples.K

      nparam  = length(model.theta)
      vpg     = zeros(dtype, nparam)
      npg     = zeros(dtype, nparam)
      grad    = zeros(dtype, nparam)

      if fullFIM
         gradll  = zeros(dtype, nparam, T*K)
         FIM     = zeros(dtype, nparam, nparam)
      else
         gradll  = zeros(dtype, 1, 1)
         FIM     = zeros(dtype, 1, 1)
      end

      adv     = zeros(dtype, T*K)
      ret     = zeros(dtype, T*K)

      tmpobs  = rand(dtype, model.mjsys.ns, T*K)
      tmpctrl = rand(dtype, model.mjsys.nu, T*K)

      new(vpg, npg, grad,
          fullFIM, gradll, FIM, adv, ret,
          baseline,
          norm_step_size,
          gamma, gae,
          cg_iter, cg_reg, cg_tol,
          tmpobs, tmpctrl,
          funcset,
          (x,c,o)->Policy.getaction!(c, model.policy, o),
          (x,c,o)->Policy.getmean!(c, model.policy, o)
         )
   end
end
Base.show(io::IO, s::NPGStrategy) = print(io, "NPG: $(s.fullFIM ? "FIM" : "HVP"), δ=$(s.norm_step_size)\n")

function stochasticrollout(mjsys::mjWrap.mjSet,
                           samples::mjWrap.TrajSamples,
                           ctrlfunc!::Function, f::FunctionSet, iter::Int)
   # stochastic rollout
   f.setmodel!(mjsys, iter) # change model for this iteration
   f.initstate!(mjsys, samples)
   f.setcontrol!(mjsys, samples.ctrl) # sets policy noise in ctrl
   @time mjWrap.rollout2(mjsys, samples,
                        ctrlfunc!, f.observe!, f.reward)
end

function update!(model::PolicyGradModel, s::NPGStrategy, iter, null)
   println("policy rollout, stoc & mean")

   mjsys       = model.mjsys
   samples     = model.samples
   meansamples = model.meansamples

   @time begin
   stochasticrollout(mjsys, samples, s.policyaction, s.f, iter)

   # policy's mean rollout
   s.f.initstate!(mjsys, meansamples)
   mjWrap.rollout2(mjsys, meansamples,
                  s.policymean,
                  s.f.observe!, s.f.reward)
   end
end

function grad_builder(f::Function, input)
   ctp = ReverseDiff.GradientTape(f, input)
   if length(ctp.tape) <= 10000
      @info("Compiling tape")
      ctp = ReverseDiff.compile(ctp)
   end
   return (out,in)->ReverseDiff.gradient!(out, ctp, in) # taped
end

# feature of learning strats; need to have iter and null to get iteration num
function hook(s::NPGStrategy, model::PolicyGradModel, iter)

   samples  = model.samples
   mjsys    = model.mjsys
   pol      = model.policy

   b        = s.baseline

   nsamples = samples.T*samples.K

   # Parallize here......................................................
   rsum = sum(samples.reward, dims=1) # rollouts done in different strat

   if eltype(pol.theta) == Float64
      ctrlmat = reshape(samples.ctrl, mjsys.nu, nsamples)
      obsmat  = reshape(samples.obs, mjsys.ns, nsamples)
   else
      copyto!(s.ctrlmat, convert(Matrix{Float32}, reshape(samples.ctrl, mjsys.nu, nsamples)))
      copyto!(s.obsmat, convert(Matrix{Float32}, reshape(samples.obs, mjsys.ns, nsamples)))
      ctrlmat = s.ctrlmat
      obsmat  = s.obsmat
   end

   println("returns")
   @time Tools.compute_returns!(s.ret, samples.reward, s.gamma)

   println("predict")
   @time Baseline.predict!(b, reshape(samples.obs, mjsys.ns, nsamples))

   println("advantages")
   @time Tools.computeGAEadvantages!(s.adv, b.bline,
                                     samples.reward,
                                     reshape(samples.obs, mjsys.ns, nsamples), s.gae, s.gamma)

   println("fit")
   @time Baseline.fit!(b, s.ret) # update fit

   if s.fullFIM
      DT = eltype(s.vpg)

      @info("FULL FIM")
      println("gradll")
      @time Policy.gradLL!(s.gradll, pol, obsmat, ctrlmat)

      println("vpg")
      # vpg[:] = gradll*adv / nsamples
      @time BLAS.gemv!('N', 1/DT(nsamples), s.gradll, s.adv, DT(0.0), s.vpg)

      println("fim")
      #FIM[:] = gradll*gradll' / nsamples
      @time begin
         #const α = 0.9
         #FIM = copy(s.FIM)
         mul!(s.FIM, s.gradll, transpose(s.gradll))
         rmul!(s.FIM, 1/DT(nsamples))
         #@. s.FIM = FIM + α*(s.FIM - FIM) # keep older FIM data around; works for higher α
      end

      println("cg")
      @time s.npg[:] = Tools.cpu_cg_solve(s.FIM, s.vpg, s.cg_iter, s.cg_reg, s.cg_tol)
      #@time IterativeSolvers.cg!(s.npg, s.FIM, s.vpg; tol = s.cg_tol, maxiter = s.cg_iter)
   else
      @info("HVP algorithm")
      @time my∇LLxA = grad_builder( (theta, a)->Policy.loglikelihood(theta, pol, obsmat, ctrlmat)*a,
                                   (pol.theta, s.adv) )
      println("vpg")
      @time begin
         my∇LLxA((s.vpg, s.adv), (pol.theta, s.adv))
         rmul!(s.vpg, 1/nsamples)
      end

      # TODO don't need in-place Loglikelihood for output; need inplace during comp
      minus = similar(s.adv')
      alpha = similar(s.adv')

      function hvp!(z, vec)
         ϵ = sqrt(eps(eltype(vec))) * norm(vec)
         eps_vec = ϵ * vec
         Policy.loglikelihood!(alpha, pol.theta+eps_vec, pol, obsmat, ctrlmat)
         Policy.loglikelihood!(minus, pol.theta-eps_vec, pol, obsmat, ctrlmat)
         alpha .-= minus
         rmul!(alpha, 1.0/2ϵ)
         s.adv .= alpha'
         my∇LLxA((z, s.adv), (pol.theta, s.adv))
         #rmul!(z, 1/nsamples)
      end
      println("cg")
      @time s.npg[:] = Tools.hvp_cg_solve(hvp!, s.vpg, s.cg_iter, s.cg_reg, s.cg_tol)
      rmul!(s.npg, nsamples)
   end

   alpha = sqrt(s.norm_step_size / dot(s.npg, s.vpg))
   println("grad")
   s.grad .= alpha.*s.npg

   println("update")
   Policy.updatetheta(pol, pol.theta + s.grad) # update parameters
   model.theta[:] = pol.theta

   ########################################################## mean rollout
   println("no noise")
   meanreward = sum(model.meansamples.reward, 1)

   ############################################################### cleanup
   iterscore = mean(rsum)
   meanscore = mean(meanreward)
   if iter > 1
      @printf "\n  Stoc: %f, Mean: %f" iterscore meanscore
   else
      @printf "\n  Mean of rewards: %f" iterscore
   end
   @printf "\n         NPG Norm: %f\n\n" norm(s.npg) #normgrad[iter]
   @printf "  Reward of Mean Rollout: %f\n" meanscore
   #@printf "  Total time: %f\n" itertime
   push!(model.trace[:stocR], iterscore)
   push!(model.trace[:meanR], meanscore)
end


text:
# This file was generated, do not modify it. # hide
no_miss = .!ismissing.(y);


text:
# Ref: https://github.com/hofmannmartin/SphericalHarmonicExpansions.jl
using MultivariatePolynomials, TypedPolynomials
import KitBase

#--- Yₗᵐ in Cartesian---#
@polyvar x y z

KitBase.ylm(0, 0, x, y, z) # 1 / 2√π = 0.282
KitBase.ylm(1, 0, x, y, z) # √(3/π) / 2 = 0.4886
KitBase.ylm(1, 1, x, y, z)
KitBase.ylm(2, 2, x, y, z)
f = KitBase.ylm(3, 3, x, y, z)

# evaluate at quadrature point
f(x => 0.5, y => -1.0, z => 0.25) # MultivariatePolynomials.jl


text:
using OpenQuantumBase, Test

H = build_example_hamiltonian(1)

@test size(H) == (2, 2)
@test H(0.5) == π * (σx + σz)
@test evaluate(H, 0.5) == (σx + σz) / 2
@test get_cache(H) ≈ π * (σx + σz)

# update_cache method
C = similar(σz)
update_cache!(C, H, 10, 0.5)
@test C == -1im * π * (σx + σz)

# update_vectorized_cache method
C = get_cache(H)
C = C⊗C
update_vectorized_cache!(C, H, 10, 0.5)
temp = -1im * π * (σx + σz)
@test C == σi ⊗ temp - transpose(temp) ⊗ σi

# in-place update for matrices
du = [1.0 + 0.0im 0; 0 0]
ρ  = PauliVec[1][1] * PauliVec[1][1]'
H(du, ρ, 2, 0.5)
@test du ≈ -1.0im * π * ((σx + σz) * ρ - ρ * (σx + σz))

# eigen-decomposition
w, v = eigen_decomp(H, 0.5)
@test w ≈ [-1, 1] / √2
@test abs(v[:, 1]'*[1-sqrt(2), 1] / sqrt(4-2*sqrt(2))) ≈ 1
@test abs(v[:, 2]'*[1+sqrt(2), 1] / sqrt(4+2*sqrt(2))) ≈ 1

Hrot= rotate(H, v)
@test evaluate(Hrot, 0.5) ≈ [-1 0; 0 1] / sqrt(2)

# error message test
@test_throws ArgumentError DenseHamiltonian([(s)->1-s, (s)->s], [σx, σz], unit=:hh)


text:
"""
    forward_filter(Z, δ, c₀)

Performs the forward filtering algorithm for a Dirichlet Evolutional Process, as
described in Fonseca & Ferreira (2017), and returns the online Dirichlet
parameters `c[1], …, c[T]`.
"""
function forward_filter(Z::Matrix{<: Integer},
                        δ::RT,
                        c₀::Vector{RT} = fill(RT(0.1), size(Z,2))) where RT <: Real

    T, k = size(Z)
    c = Matrix{RT}(undef, T, k)

    c[1,:] = δ * c₀ + Z[1,:]
    for t = 2:T
        c[t,:] = δ * c[t-1,:] + Z[t,:]
    end

    return c
end

"""
    backwards_sampler(c, δ)

Performs the backwards sampling algorithm for a Dirichlet Evolutional Process,
as described in Fonseca & Ferreira (2017), and returns samples `η[1], …, η[T]`
from the posterior distribution.
"""
function backwards_sampler(c::Matrix{RT}, δ::RT) where RT <: Real

    T = size(c, 1)
    η = similar(c)

    η[T,:] = rand(Dirichlet(c[T,:]))
    for t = T-1:-1:1
        cₛ = sum(c[t,:])
        S = rand(Beta(δ * cₛ, (1 - δ) * cₛ))
        u = rand(Dirichlet((1 - δ) * c[t,:]))
        η[t,:] = S * η[t+1,:] + (1 - S) * u
    end

    return η
end

"""
    moddirichlet_mean(c, a, b)

Internal function computing the mean for a Mod-Dirichlet.
"""
@inline function moddirichlet_mean(c::Vector{RT},
                                   a::Vector{RT},
                                   b::Vector{RT}) where RT <: Real

    return (b - a) .* (c / sum(c)) + a
end

"""
    moddirichlet_params(c, δ, η)

Internal function to return the parameters for the Mod-Dirichlet distribution.
"""
function moddirichlet_params(c::Vector{RT},
                             δ::RT,
                             η::Vector{RT}) where RT <: Real

    k = size(η, 1)

    # Compute mode for S

    cₛ = sum(c)
    α = δ * cₛ
    β = (1 - δ) * cₛ

    if (α < 1.) & (β < 1.)
        throw(ArgumentError("Invalid mode for S"))
    elseif α <= 1.
        S = 0
    elseif β <= 1.
        S = 1
    else
        S = (α - 1) / (α + β - 2)
    end

    # Compute parameters as function of the mode of S

    c = (1 - δ) * c
    a = S * η
    b = (1 - S) * ones(k) + S * η

    return c, a, b
end

"""
    backwards_estimator(c, δ)

Similar to `backwards_sampler`, but obtains modes in a backwards step instead
of samples.
"""
function backwards_estimator(c::Matrix{RT}, δ::RT) where RT <: Real

    T, k = size(c)
    η = similar(c)

    # For the last eta's distribution we can do it directly since it's a known
    # Dirichlet(c[T]).

    η[T,:] = moddirichlet_mean(c[T,:], zeros(k), ones(k))

    # For the other time instants we need to find the Mod-Dirichlet parameters
    # conditional to the previous mode being a known value.

    for t = T-1:-1:1
        cₜ, aₜ, bₜ = moddirichlet_params(c[t,:], δ, η[t+1,:])
        η[t,:] = moddirichlet_mean(cₜ, aₜ, bₜ)
    end

    return η
end