text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module SPQR

import Base: \

# ordering options */
const ORDERING_FIXED   = Int32(0)
const ORDERING_NATURAL = Int32(1)
const ORDERING_COLAMD  = Int32(2)
const ORDERING_GIVEN   = Int32(3) # only used for C/C++ interface
const ORDERING_CHOLMOD = Int32(4) # CHOLMOD best-effort (COLAMD, METIS,...)
const ORDERING_AMD     = Int32(5) # AMD(A'*A)
const ORDERING_METIS   = Int32(6) # metis(A'*A)
const ORDERING_DEFAULT = Int32(7) # SuiteSparseQR default ordering
const ORDERING_BEST    = Int32(8) # try COLAMD, AMD, and METIS; pick best
const ORDERING_BESTAMD = Int32(9) # try COLAMD and AMD; pick best#

# Let [m n] = size of the matrix after pruning singletons.  The default
# ordering strategy is to use COLAMD if m <= 2*n.  Otherwise, AMD(A'A) is
# tried.  If there is a high fill-in with AMD then try METIS(A'A) and take
# the best of AMD and METIS.  METIS is not tried if it isn't installed.

# tol options
const DEFAULT_TOL = Int32(-2) # if tol <= -2, the default tol is used
const NO_TOL      = Int32(-1) # if -2 < tol < 0, then no tol is used

# for qmult, method can be 0,1,2,3:
const QTX = Int32(0)
const QX  = Int32(1)
const XQT = Int32(2)
const XQ  = Int32(3)

# system can be 0,1,2,3:  Given Q*R=A*E from SuiteSparseQR_factorize:
const RX_EQUALS_B    = Int32(0) # solve R*X=B      or X = R\B
const RETX_EQUALS_B  = Int32(1) # solve R*E'*X=B   or X = E*(R\B)
const RTX_EQUALS_B   = Int32(2) # solve R'*X=B     or X = R'\B
const RTX_EQUALS_ETB = Int32(3) # solve R'*X=E'*B  or X = R'\(E'*B)

using ..SparseArrays: SparseMatrixCSC
using ..SparseArrays.CHOLMOD: C_Dense, C_Sparse, Dense, ITypes, Sparse, SuiteSparseStruct, VTypes, common

import Base: size
import Base.LinAlg: qrfact
import ..SparseArrays.CHOLMOD: convert, free!

immutable C_Factorization{Tv<:VTypes} <: SuiteSparseStruct
    xtype::Cint
    factors::Ptr{Tv}
end

type Factorization{Tv<:VTypes} <: Base.LinAlg.Factorization{Tv}
    m::Int
    n::Int
    p::Ptr{C_Factorization{Tv}}
    function Factorization(m::Integer, n::Integer, p::Ptr{C_Factorization{Tv}})
        if p == C_NULL
            throw(ArgumentError("factorization failed for unknown reasons. Please submit a bug report."))
        end
        new(m, n, p)
    end
end
Factorization{Tv<:VTypes}(m::Integer, n::Integer, p::Ptr{C_Factorization{Tv}}) = Factorization{Tv}(m, n, p)

size(F::Factorization) = (F.m, F.n)
function size(F::Factorization, i::Integer)
    if i < 1
        throw(ArgumentError("dimension must be positive"))
    end
    if i == 1
        return F.m
    elseif i == 2
        return F.n
    end
    return 1
end

function free!{Tv<:VTypes}(F::Factorization{Tv})
    ccall((:SuiteSparseQR_C_free, :libspqr), Cint,
        (Ptr{Ptr{C_Factorization{Tv}}}, Ptr{Void}),
            &F.p, common()) == 1
end

function backslash{Tv<:VTypes}(ordering::Integer, tol::Real, A::Sparse{Tv}, B::Dense{Tv})
    m, n  = size(A)
    if m != size(B, 1)
        throw(DimensionMismatch("left hand side and right hand side must have same number of rows"))
    end
    d = Dense(ccall((:SuiteSparseQR_C_backslash, :libspqr), Ptr{C_Dense{Tv}},
        (Cint, Cdouble, Ptr{C_Sparse{Tv}}, Ptr{C_Dense{Tv}}, Ptr{Void}),
            ordering, tol, get(A.p), get(B.p), common()))
    finalizer(d, free!)
    d
end

function factorize{Tv<:VTypes}(ordering::Integer, tol::Real, A::Sparse{Tv})
    s = unsafe_load(A.p)
    if s.stype != 0
        throw(ArgumentError("stype must be zero"))
    end
    f = Factorization(size(A)..., ccall((:SuiteSparseQR_C_factorize, :libspqr), Ptr{C_Factorization{Tv}},
        (Cint, Cdouble, Ptr{Sparse{Tv}}, Ptr{Void}),
            ordering, tol, get(A.p), common()))
    finalizer(f, free!)
    f
end

function solve{Tv<:VTypes}(system::Integer, QR::Factorization{Tv}, B::Dense{Tv})
    m, n = size(QR)
    mB = size(B, 1)
    if (system == RX_EQUALS_B || system == RETX_EQUALS_B) && m != mB
        throw(DimensionMismatch("number of rows in factorized matrix must equal number of rows in right hand side"))
    elseif (system == RTX_EQUALS_ETB || system == RTX_EQUALS_B) && n != mB
        throw(DimensionMismatch("number of columns in factorized matrix must equal number of rows in right hand side"))
    end
    d = Dense(ccall((:SuiteSparseQR_C_solve, :libspqr), Ptr{C_Dense{Tv}},
        (Cint, Ptr{C_Factorization{Tv}}, Ptr{C_Dense{Tv}}, Ptr{Void}),
            system, get(QR.p), get(B.p), common()))
    finalizer(d, free!)
    d
end

function qmult{Tv<:VTypes}(method::Integer, QR::Factorization{Tv}, X::Dense{Tv})
    mQR, nQR = size(QR)
    mX, nX = size(X)
    if (method == QTX || method == QX) && mQR != mX
        throw(DimensionMismatch("Q matrix size $mQR and dense matrix has $mX rows"))
    elseif (method == XQT || method == XQ) && mQR != nX
        throw(DimensionMismatch("Q matrix size $mQR and dense matrix has $nX columns"))
    end
    d = Dense(ccall((:SuiteSparseQR_C_qmult, :libspqr), Ptr{C_Dense{Tv}},
            (Cint, Ptr{C_Factorization{Tv}}, Ptr{C_Dense{Tv}}, Ptr{Void}),
                method, get(QR.p), get(X.p), common()))
    finalizer(d, free!)
    d
end

qrfact(A::SparseMatrixCSC, ::Type{Val{true}}) = factorize(ORDERING_DEFAULT, DEFAULT_TOL, Sparse(A, 0))

"""
    qrfact(A) -> SPQR.Factorization

Compute the `QR` factorization of a sparse matrix `A`. A fill-reducing permutation is used.
The main application of this type is to solve least squares problems with `\\`. The function
calls the C library SPQR and a few additional functions from the library are wrapped but not
exported.
"""
qrfact(A::SparseMatrixCSC) = qrfact(A, Val{true})

# With a real lhs and complex rhs with the same precision, we can reinterpret
# the complex rhs as a real rhs with twice the number of columns
#
# This definition is similar to the definition in factorization.jl except that
# here we have to use \ instead of A_ldiv_B! because of limitations in SPQR

## Two helper methods
_ret_size(F::Factorization, b::AbstractVector) = (size(F, 2),)
_ret_size(F::Factorization, B::AbstractMatrix) = (size(F, 2), size(B, 2))

function (\)(F::Factorization{Float64}, B::VecOrMat{Complex{Float64}})
# |z1|z3|  reinterpret  |x1|x2|x3|x4|  transpose  |x1|y1|  reshape  |x1|y1|x3|y3|
# |z2|z4|      ->       |y1|y2|y3|y4|     ->      |x2|y2|     ->    |x2|y2|x4|y4|
#                                                 |x3|y3|
#                                                 |x4|y4|
    c2r = reshape(transpose(reinterpret(Float64, B, (2, length(B)))), size(B, 1), 2*size(B, 2))
    x = F\c2r

# |z1|z3|  reinterpret  |x1|x2|x3|x4|  transpose  |x1|y1|  reshape  |x1|y1|x3|y3|
# |z2|z4|      <-       |y1|y2|y3|y4|     <-      |x2|y2|     <-    |x2|y2|x4|y4|
#                                                 |x3|y3|
#                                                 |x4|y4|
    return reinterpret(Complex{Float64}, transpose(reshape(x, (length(x) >> 1), 2)), _ret_size(F, B))
end

function (\){T<:VTypes}(F::Factorization{T}, B::StridedVecOrMat{T})
    QtB = qmult(QTX, F, Dense(B))
    convert(typeof(B), solve(RETX_EQUALS_B, F, QtB))
end

(\)(F::Factorization, B::StridedVecOrMat) = F\convert(AbstractArray{eltype(F)}, B)

end # module

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
#=
    The File contains Sparsik tests
=#

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

@testset "Sparsik - Base functionality" begin

    RR = Nemo.GF(7)
    BIT = Nemo.GF(2)

    v = Sparsik(7, RR, [1, 3, 5, 7],
        Dict(1 => RR(10), 3 => RR(-2), 5 => RR(8), 7 => RR(4)))
    @test v == from_dense([3, 0, 5, 0, 1, 0, 4], RR)
    v = Sparsik(2, Nemo.QQ, Int[], Dict{Int, elem_type(Nemo.QQ)}())
    @test v == from_dense([0, 0], Nemo.QQ)

    v = from_dense([1, 2, 3], Nemo.QQ)
    @test scale!(v, -2) == from_dense([-2, -4, -6], Nemo.QQ)
    @test scale(v, 0) == from_dense([0, 0, 0], Nemo.QQ)
    @test scale!(v, 2) == from_dense([-4, -8, -12], Nemo.QQ)
    v = from_dense([1, 2, 3, 4], RR)
    @test iszero(scale(v, 14))

    @test iszero(from_dense([3, 6, 9, -3, 0], Nemo.GF(3)))
    @test iszero(from_dense([0, 0, 0], RR))
    @test !iszero(from_dense([0, 1, 0], RR))
    v = from_dense([0, 0, 1, 10], RR)
    @test first_nonzero(v) == 3
    v = from_dense([3, 0, 0, 1, 1], Nemo.GF(3))
    @test first_nonzero(v) == 4
    v = from_dense([0, 0], Nemo.QQ)
    @test first_nonzero(v) == -1

    v, u = from_dense([1, 0, 3, 0, 5], Nemo.QQ), from_dense([0, -1, -4, 3, 3], Nemo.QQ)
    @test inner(v, u) == 3
    v, u = from_dense([1, 2, 3], Nemo.QQ), from_dense([-1, -2, 1], Nemo.QQ)
    @test inner(v, u) == -2
    v, u = from_dense([1, 2, 3], BIT), from_dense([1, 2, 3], BIT)
    @test inner(v, u) == 0

    v, u, c = from_dense([1, 2, 3], Nemo.QQ), from_dense([2, 1, 1], Nemo.QQ), -1
    @test reduce(v, u, c) == from_dense([-1, 1, 2], Nemo.QQ)
    v, u, c = from_dense([2, 0, 6], Nemo.QQ), from_dense([1, 0, 3], Nemo.QQ), -2
    @test iszero(reduce(v, u, c))
    @test iszero(reduce(v, v, -1))
    v = from_dense([2, 0, 7, 14, 0, 5], RR)
    u = from_dense([-2, 1, 1, 7, 1, -5], RR)
    c = 1
    @test length(reduce!(v, u, c)) == 3
    @test scale!(v, -1) == from_dense([0, -1, -1, 0, 6, 0], RR)
    @test scale!(v, 7) == zero_sparsik(size(v)..., RR)
    @test v + zero_sparsik(size(v)..., RR) == zero_sparsik(size(v)..., RR)

    v = from_dense([1, 0, 0, 0, 0, 7, -7, 3, 1, 1, 0, 4, 5, 11], RR)
    u = from_dense([0, 0, 3, 7, 1, -1, 7, 4, 0, 6, 0, 0, 2, 3], RR)
    r = from_dense([1, 0, 3, 0, 1, 6, 0, 0, 1, 0, 0, 4, 0, 0], RR)
    @test reduce(v, u, 1) == r
    @test length(r.data) == length(r.nonzero) == 6

    items = [(i, x) for (i, x) in r]
    @test items == [(1, RR(1)), (3, RR(3)), (5, RR(1)), (6, RR(6)),
                    (9, RR(1)), (12, RR(4))]

end

@testset "Sparsik - modular" begin

    R7 = Nemo.GF(7)
    R37 = Nemo.GF(37)
    BIT = Nemo.GF(2)

    diff_rationals(x, y) = numerator(x)*denominator(y) - numerator(y)*denominator(x)

    v = from_dense([1, 0, 0, 4, 14, 13], Nemo.QQ)
    u = from_dense([1, 0, 0, 4, 0, 6], R7)
    @test modular_reduction(v, R7) == u
    r = rational_reconstruction(u)
    @test all(i -> iszero(R7(diff_rationals(r[i], v[i]))), 1:dim(v))

    v = from_dense([11//27, 4//13, 1//1, 0, 22//113, 7//2, 4//3], Nemo.QQ)

    u = modular_reduction(v, R7)
    r = rational_reconstruction(u)
    @test all(i -> iszero(R7(diff_rationals(r[i], v[i]))), 1:dim(v))

    u = modular_reduction(v, R37)
    r = rational_reconstruction(u)
    @test all(i -> iszero(R37(diff_rationals(r[i], v[i]))), 1:dim(v))

end

@testset "Sparsik - Advanced functionality" begin

    # to_dense tests
    sz = [1, 2, 3, 5]
    densities = [0.5, 1.0]
    ZZ = Nemo.GF(2^31 - 1)

    # to_dense tests
    for n in sz
        for λ in densities
            v = random_sparsik(n, ZZ, density=λ)
            @test v == from_dense(to_dense(v), ZZ)
        end
    end

end

@info "OK"

ext: jl


text:
module RoundingEmulator

export add_up, add_down, sub_up, sub_down, mul_up, mul_down, div_up, div_down, sqrt_up, sqrt_down

include("rounding.jl")

end # module


text:
# ------------------------------------------------------------------------------------------
# # 绘图
#
# ## 基础
# Julia 有包括调用 `PyPlot` 在内的、多种不同的绘图方式。
#
# 这里将介绍如何使用 `Plots.jl`。
# 如果还没安装 `Plots.jl`，你需要通过包管理器（package manager）来安装它，Julia会在你第一次导入它的时候将它预编译：
# ------------------------------------------------------------------------------------------

# using Pkg
# Pkg.add("Plots")
using Plots

# ------------------------------------------------------------------------------------------
# `Plots.jl` 的优势之一是可以无缝地切换后端（backends）。
# 在这个 notebook 中，我们将尝试 `gr()` 和 `plotlyjs()` 后端。
#
# 以科学调查之名，我们用这个 notebook 调查一下大概1860到2000年的全球温度和海盗数量的关系。
# ------------------------------------------------------------------------------------------

globaltemperatures = [14.4, 14.5, 14.8, 15.2, 15.5, 15.8]
numpirates = [45000, 20000, 15000, 5000, 400, 17];

# ------------------------------------------------------------------------------------------
# Plots 通过同样的接口可以支持多种后端——指的是实际执行绘制的库。
# 刚开始，我们先使用 GR 后端。通过调用 `gr()` 来选择 GR 后端：
# ------------------------------------------------------------------------------------------

gr()

# ------------------------------------------------------------------------------------------
# 现在我们可以调用如 `plot` 和 `scatter` 的函数来绘制图像。
# ------------------------------------------------------------------------------------------

plot(numpirates, globaltemperatures, label="line")
scatter!(numpirates, globaltemperatures, label="points")

# ------------------------------------------------------------------------------------------
# `scatter!` 函数名后面的 `!` 意味着它是一个原地修改传入变量的函数，表示散点图会被添加到已存在的图像上。
#
# 与此对应的，你可以试试换成 `scatter` 看看会发生什么。
#
# 接着，我们通过 `xlabel!`，`ylabel!` 和 `title!` 函数来给图像加上更多信息。
# ------------------------------------------------------------------------------------------

xlabel!("Number of Pirates [Approximate]")
ylabel!("Global Temperature (C)")
title!("Influence of pirate population on global warming")

# ------------------------------------------------------------------------------------------
# 这看上去仍不对劲。自1860以来海盗数量是减少的，而从左往右看其实时间上是倒序的。
# 我们来把 X 轴反过来，可以更清楚地看到时间顺序下海盗人口是如何导致全球温度的变化！
# ------------------------------------------------------------------------------------------

xflip!()

# ------------------------------------------------------------------------------------------
# 绘图就是这样简单！
#
# 注意：这是一个关于人们是如何经常结合相关性和因果性的笑话。
#
# **不需要修改语法，我们可以在 `unicodeplots()` 后端中绘制同样的图像**
# ------------------------------------------------------------------------------------------

using Pkg; Pkg.add("UnicodePlots")
unicodeplots()

plot(numpirates, globaltemperatures, label="line")
scatter!(numpirates, globaltemperatures, label="points")
xlabel!("Number of Pirates [Approximate]")
ylabel!("Global Temperature (C)")
title!("Influence of pirate population on global warming")

# ------------------------------------------------------------------------------------------
# 注意到第二幅图和第一幅的区别！
# 既然 Jupyter notebook 有酷炫的绘图能力，使用文本符号来绘图就显得有些傻了。
# 但其实这在终端中进行快速简陋的可视化是很实用的。
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 练习Exercises
#
# ### 8.1
# 已有
# ```julia
# x = -10:10
# ```
# 要绘制 y 随 x 变化曲线，其中 $y = x^2$。你可能需要先把后端切换回去。
# ------------------------------------------------------------------------------------------

x = -10:10

# 在此作答

# ------------------------------------------------------------------------------------------
# ### 8.2
# 执行以下代码
# ------------------------------------------------------------------------------------------

p1 = plot(x, x)
p2 = plot(x, x.^2)
p3 = plot(x, x.^3)
p4 = plot(x, x.^4)
plot(p1, p2, p3, p4, layout = (2, 2), legend = false)

# ------------------------------------------------------------------------------------------
# 然后使用 `p1`、`p2`、`p3` 和 `p4` 作为子图创建一个四行一列的图。
#
# 提示：你可以参考 Plots 的帮助文档 [GR - Plots](http://docs.juliaplots.org/latest/examples/gr/#magic-
# grid-argument)
# ------------------------------------------------------------------------------------------

# 在此作答

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
using CUDAnative
using CUDAdrv
include("conv_2d.jl")
include("conv_2d_cuda.jl")

img_lena=rand(Float32,32,72,72,32)
fil_lena=rand(Float32,32,3,3,32)
dest_lena=zeros(Float32,32,72,72,32)
conv_2d_combined!(img_lena,fil_lena,dest_lena,(1,1),(1,1))

img_cuda=CuArray(permutedims(img_lena,(4,2,3,1)))
fil_cuda=CuArray(permutedims(fil_lena,(4,2,3,1)))
layer=conv2dWinograd(fil_cuda,img_cuda)
forward(layer,img_cuda)
dest=permutedims(Array(layer.dest),(4,2,3,1))

max_stars_repo_name: JerryLingjieMei/WinogradConvolution

max_issues_repo_name: JerryLingjieMei/WinogradConvolution

max_forks_repo_name: JerryLingjieMei/WinogradConvolution


text:
module TestDevice

include("../device_tests.jl")
using CUDAKernels

test_all(CUDADevice())

end  # module


text:
function gradcheck(m, grad, loss; gcheck=10, _eps=cbrt(eps(eltype(m))), delta=_eps, atol=_eps, rtol=_eps, o...)
    # 6e-6 for Float64, 5e-3 for Float32 works best
    # rnum = 42 # time_ns() #DBG
    # isdefined(data,:rng) && (data_rng_save = data.rng; data.rng=MersenneTwister(0); srand(data.rng,rnum))
    # isdefined(m,:rng) && (m_rng_save = m.rng; m.rng=MersenneTwister(0); srand(m.rng,rnum))
    # l = zeros(2)
    # train(m, data, loss; gcheck=true, losscnt=fill!(l,0), o...)
    # loss0 = l[1]
    loss0 = loss(m)
    grad(m)
    pp = params(m)
    dcopy(p)=(isa(p.dif,Array) ? copy(p.dif) : convert(Array,p.dif))
    ppdif = map(dcopy, pp) # This may get reset in loss, so make a copy now
    maxbad = 0
    for n = 1:length(pp)
        p = pp[n]
        pdif = ppdif[n]
        psave = copy(p.out)     # DBG
        wlen = length(p.out)
        irange = (wlen <= gcheck ? (1:wlen) : sortperm(abs(vec(pdif)),rev=true)[1:gcheck]) # rand(1:wlen, gcheck))
        for i in irange
            wi0 = p.out[i]
            wi1 = (wi0 >= 0 ? wi0 + delta : wi0 - delta)
            p.out[i] = wi1
            # isdefined(data,:rng) && srand(data.rng, rnum)
            # isdefined(m,:rng) && srand(m.rng, rnum)
            # train(m, data, loss; gcheck=true, losscnt=fill!(l,0), o...)
            # loss1 = l[1]
            loss1 = loss(m)
            p.out[i] = wi0
            dw0 = pdif[i]
            dwi = (loss1 - loss0) / (wi1 - wi0)
            z=abs(dw0-dwi)-rtol*max(abs(dw0),abs(dwi)); z>maxbad && (maxbad=z)
            if !isapprox(dw0, dwi; rtol=rtol, atol=atol)
                println("gc:fail@$(findfirst(params(m),p))/$i/$wlen wi:$wi0->$wi1 loss:$loss0->$loss1 dw_bprop=$dw0 dw_numeric=$dwi")
            end
        end
        @assert isequal(p.out, psave) # DBG
    end
    # isdefined(data,:rng) && (data.rng = data_rng_save)
    # isdefined(m,:rng) && (m.rng = m_rng_save)
    @printf("gc:atol=%g for rtol=%g, delta=%g for largest %d x %d gradients: abs(x-y) <= atol+rtol*max(abs(x),abs(y))\n",
            maxbad, rtol, delta, length(params(m)), gcheck)
    # println("gc:atol=$maxbad for rtol=$rtol, delta=$delta for largest $(length(params(m)))x$(gcheck) gradients: abs(x-y) <= atol+rtol*max(abs(x),abs(y))")
end

# Why loss0 = l[1], instead of l[1]/l[2]?
# gcheck=true ensures we process one item/sequence. In case of item
# l[2]=1 so it doesn't matter.  In case of a sequence l[2]=tokencnt,
# but we want the total loss, because that is what determined dw.


text:
"""
    write_figure(f)

Internal function to generate and write the GLE script associated with the
figure object `f`.
"""
function write_figure(
            f::Figure,
            opath::String = ""
        )::String

    g = f.g
    "size $(f.size[1]) $(f.size[2])" |> g

    # >> apply background color if different than nothing
    # Note that if it's white, we apply it as this makes a
    # difference for transparent outputs like SVG.
    if isdef(f.bgcolor)
        # add a box that is slightly larger than the size
        "\namove -0.05 -0.05" |> g
        "\nbox $(f.size[1]+0.1) $(f.size[2]+0.1)" |> g
        "fill $(col2str(f.bgcolor)) nobox" |> g
    end

    # >> apply latex
    # check if has latex
    haslatex = any(isdef, (f.texscale, f.texpreamble))
    isdef(f.texlabels) && (haslatex = f.texlabels)
    # line for texstyle, it may be empty if nothing is given
    apply_textstyle!(g, f.textstyle, addset=true)
    # latex if required
    if haslatex
        if isdef(f.texpreamble)
            "\nbegin texpreamble\n" |> g
            f.texpreamble           |> g
            "\nend texpreamble"     |> g
        end
        "\nset texlabels 1" |> g
        "\nset texscale"    |> g
        ifelse(isdef(f.texscale), f.texscale, "scale") |> g
    end
    write(g.io, "\n")

    # >> apply axes
    # NOTE: this organisation is so that if axes need extra
    # subroutines, these will be generated after applying axes
    # but need to be put before in the GLE script
    gtemp = GS();
    for (i, aᵢ) ∈ enumerate(f.axes)
        apply_axes!(gtemp, aᵢ, f.id, i)
    end
    if !isempty(f.subroutines)
        ks = keys(f.subroutines)
        # subroutines for palettes
        for key ∈ Iterators.filter(k -> startswith(k, "cmap_"), ks)
            f.subroutines[key] |> g
        end
        # subroutines for drawings
        # -- boxplot, heatmap
        for key ∈ Iterators.filter(k -> startswith(k, "bp_") || startswith(k, "hm_"), ks)
            f.subroutines[key] |> g
        end
        # subroutines for markercolors
        for key ∈ Iterators.filter(k->startswith(k, "mk_"), ks) # special markers
            f.subroutines[key] |> g
        end
        # --------------
        # 3D subroutines
        # --------------
        "plot3" ∈ ks && f.subroutines["plot3"] |> g
    end
    gtemp |> g

    # >> either return as string or write to file
    isempty(opath) && return String(take!(g))
    write(opath, take!(g))
    return ""
end

"""
    debug_gle(f)

Print the GLE script associated with figure `f` for debugging.
"""
debug_gle(f::Figure) = println(write_figure(f))
debug_gle()          = debug_gle(gcf())


text:
using RandomBasedArrays
using StableRNGs, Test

@testset "RandomBasedArrays.jl" begin
    rng = StableRNG(42)
    p = reshape(collect(1:6), 2, 3)
    A = RandomBasedArray(p, rng)
    @test parent(A) === p
    @test parent(A) == [1 3 5; 2 4 6]
    @test size(A) === size(p)
    @test A[42] == 2
    @test A[-35] == 1
    @test A[4,0] == 5
    A[-314] = 0
    @test parent(A) == [0 3 5; 2 4 6]
    A[0,0] = 42
    @test parent(A) == [42 3 5; 2 4 6]
    @test occursin(r"2×3 Array{Int[2346]+,2}:\n 42 +6 +3\n +6 +3 +42", repr("text/plain", A))
end

max_stars_repo_name: giordano/RandomBasedArrays

max_issues_repo_name: giordano/RandomBasedArrays

max_forks_repo_name: giordano/RandomBasedArrays


text:
include("../src/fdfd_local_field.jl")
cd("examples/")

using BenchmarkTools
using DelimitedFiles

println("$(Threads.nthreads()) threads!")
number_training_data = 150000 # ideally 1e5 here

r = rand()
δ = 1e-2
fname = "test_data_$(round(Int, r*1e5))_$δ.csv"
number_holes = 10
lb = 0.1 # for a blue wavelength this corresponds to a minimum feature of ≈40 nm
ub = 0.95 - lb # this assumes a period of 0.95 λ

# create parameters
ps_array = create_parameters(number_holes, number_training_data, δ, lb, ub)

# reference simulation
refractive_indexes_ref = ones(3) * 1.45
ref_local_field = get_local_field(Float64[0.5], refractive_indexes=refractive_indexes_ref)[3]

# simulations
refractive_indexes_sim = Float64[1.0, 1.0, 1.45]

# embarassingly parallel loop (used to be serial hence the name of the file)
sim_local_fields = zeros(ComplexF64, (number_holes+1)*number_training_data)
@time @inbounds Base.Threads.@threads for it_ps=1:size(ps_array)[1]
    sim_local_fields[it_ps] = get_local_field(ps_array[it_ps, :], refractive_indexes=refractive_indexes_sim)[3]
end
local_field_data = sim_local_fields./ref_local_field

writedlm("./data/"*fname,hcat(ps_array, real.(local_field_data), imag.(local_field_data)), ',')


text:
# FHIP_mobility.jl

"""
Reproduces figures 1-3 in the 1962 FHIP paper, and extends them to a larger range of frequencies and temperatures.
"""

using PolaronPathIntegrals
using Plots
using SpecialFunctions
using LaTeXStrings
using DSP
using Images

# Interactive Figures
plotly()
Plots.PlotlyBackend()

# Static Figures
# pgfplots()

function plot_figure(
    Ω_range,
    α,
    w,
    v;
    β_range = [
        gamma(i) for i in range(1.461632144968, stop = 5.89251869634377, length = 7)
    ],
    xtickstep = 1,
    ytickstep = 1,
    ylog = false,
    zero = false,
)

    if ylog
        zero = false
    end

    if !zero
        # Imχs = [[PolaronPathIntegrals.ℑχ(Ω, β, α, v, w) for Ω in Ω_range] for β in β_range]
        Imχs = [imag(hilbert(vcat(repeat([Float64(PolaronPathIntegrals.ℑχ(Ω_range[1], β, α, v, w))], 2 * length(Ω_range) + 1), [Float64(PolaronPathIntegrals.ℑχ(Ω, β, α, v, w)) for Ω in Ω_range], repeat([Float64(PolaronPathIntegrals.ℑχ(Ω_range[end], β, α, v, w))], 2 * length(Ω_range) + 1))))[2 * length(Ω_range) + 2:end - 2 * length(Ω_range) - 1] for β in β_range]
        @show(Imχs[1])
    else
        Imχs = [vcat([0.0], [PolaronPathIntegrals.ℑχ(Ω, β, α, v, w) for Ω in Ω_range]) for β in β_range]
        Ω_range = vcat([0.0], Ω_range)
    end

    if ylog
        plot = Plots.plot(
            Ω_range,
            Imχs[1],
            xticks = Ω_range[1]:xtickstep:Ω_range[end],
            yaxis = :log,
            xlabel = "Ω",
            ylabel = "Imχ(Ω)",
            label = "\$\\beta\$ = $(round(β_range[1], digits = 3))",
            legend = :bottomright,
            size = (1800, 1200),
            linewidth = 2,
            xtickfontsize = 18,
            ytickfontsize = 18,
            xguidefontsize = 20,
            yguidefontsize = 20,
            legendfontsize = 18,
            thickness_scaling = 1.2,
        )
        annotate!(plot,  Plots.xlims(plot)[1] + 0.8 * (Plots.xlims(plot)[2] - Plots.xlims(plot)[1]), 0.5e-2, ("\$\\alpha\$ = $α, v = $v, w = $w", :black, :right, 20))
        for i = 2:length(Imχs)
            Plots.plot!(
                plot,
                Ω_range,
                Imχs[i],
                label = "\$\\beta\$ = $(round(β_range[i], digits = 3))",
                linewidth = 2,
            )
        end

    else
        plot = Plots.plot(
            Ω_range,
            Imχs[1],
            xticks = Ω_range[1]:xtickstep:Ω_range[end],
            yticks = floor(min((Imχs...)...)):ytickstep:ceil(max((Imχs...)...)),
            xlabel = "Ω",
            ylabel = "Imχ(Ω)",
            label = "\$\\beta\$ = $(round(β_range[1], digits = 3))",
            legend = :topleft,
            size = (1800, 1200),
            linewidth = 2,
            xtickfontsize = 18,
            ytickfontsize = 18,
            xguidefontsize = 20,
            yguidefontsize = 20,
            legendfontsize = 18,
            thickness_scaling = 1.2,
        )
        annotate!(plot, Plots.xlims(plot)[1] + 0.4 * (Plots.xlims(plot)[2] - Plots.xlims(plot)[1]), Plots.ylims(plot)[1] + 0.95 * (Plots.ylims(plot)[2] - Plots.ylims(plot)[1]), ("\$\\alpha\$ = $α, v = $v, w = $w", :black, :right, 20))
        for i = 2:length(Imχs)
            Plots.plot!(
                plot,
                Ω_range,
                Imχs[i],
                label = "\$\\beta\$ = $(round(β_range[i], digits = 3))",
                linewidth = 2,
            )
        end
    end
    return plot
end

"""
Replicated Plots
"""

setprecision(BigFloat, 64)

# replicated_figure_one = plot_figure(
#     range(3.0, stop = 9.0, length = 1000),
#     3.0,
#     2.5,
#     3.4;
#     β_range = [100.0],
# )
#
# replicated_figure_two = plot_figure(
#     range(2.0, stop = 22.0, length = 1000),
#     5.0,
#     2.1,
#     4.0;
#     β_range = [100.0],
#     xtickstep = 2,
#     ytickstep = 2,
# )
#
# replicated_figure_three = plot_figure(
#     range(6.0, stop = 28.0, length = 1000),
#     7.0,
#     1.6,
#     5.8;
#     β_range = [100.0],
#     xtickstep = 2,
#     ylog = true,
# )

"""
Extended Plots
"""

# extended_figure_one = plot_figure(
#     range(3.0, stop = 9.0, length = 1000),
#     3.0,
#     2.5,
#     3.4;
#     ylog = false,
#     xtickstep = 1,
#     ytickstep = 1,
#     zero = false,
# )

extended_figure_two = plot_figure(
    range(2.0, stop = 28.0, length = 1000),
    5.0,
    2.1,
    4.0;
    ylog = false,
    xtickstep = 2,
    ytickstep = 5,
    zero = false,
)
#
# extended_figure_three = plot_figure(
#     range(6.0, stop = 28.0, length = 1000),
#     7.0,
#     1.6,
#     5.8;
#     ylog = false,
#     xtickstep = 2,
#     ytickstep = 5,
#     zero = false,
# )

# figpath = "C:/Users/neutr/OneDrive - Imperial College London/PhD/Code/PolaronPathIntegrals/examples/plots/"

# for i in ["pdf", "tex", "svg"]
    # savefig(replicated_figure_one, figpath * "FHIP_replicated_fig_1.$i")
    # savefig(replicated_figure_two, figpath * "FHIP_replicated_fig_2.$i")
    # savefig(replicated_figure_three, figpath * "FHIP_replicated_fig_3.$i")
#     savefig(extended_figure_one, figpath * "FHIP_extend_fig_1_short.$i")
#     savefig(extended_figure_two, figpath * "FHIP_extend_fig_2_short.$i")
#     savefig(extended_figure_three, figpath * "FHIP_extend_fig_3_short.$i")
# end


text:
""" Example usage :
@observation MyModel begin
    Categorical(epsilon_argmax(ucb(Q, U, c)))
end
or
@observation MyModel begin
    Categorical(epsilon_greedy(softmax(ucb(Q, U, c))))
end
"""

function epsilon_greedy!(P, ϵ)
    N = length(P)
    P .*= 1 - ϵ
    P .+= ϵ/(N)
    return P
end

function epsilon_argmax(Q::Vector{T}, ϵ) where T
    P = zero(Q)
    P[argmax(Q)] = one(T)
    return epsilon_greedy!(P, ϵ)
end

function ucb!(Q, U, c)
    Q .+= c * sqrt.(U)
end

max_stars_repo_name: sqwayer/AnimalBehavior

max_issues_repo_name: sqwayer/AnimalBehavior


text:
"""
CrystalViscoPlastic(E, ν, τ_y, H_iso, H_kin, q, α_∞, t_star, σ_c, m, slipsystems)

Crystal plasticity with
# Arguments
- `E::Float64`: Young's modulus
- `ν::Float64`: Poisson's ratio
- `τ_y`: critical resolved shear stress
- `H_iso`: hardening modulus for isotropic hardening
- `H_kin`: hardening modulus for kinematic hardening
- `q`: cross-hardening parameter for isotropic hardening
- `α_∞`: saturation stress for kinematic hardening
- `t_star`: visco-plastic characteristic time
- `σ_c`: visco-plastic characteristic stress
- `m`: visco-plastic exponent
- `slipsystems`: slipsystems
"""
struct CrystalViscoPlastic{S} <: AbstractMaterial
    E::Float64
    ν::Float64
    τ_y::Float64
    H_iso::Float64
    H_kin::Float64
    q::Float64
    α_∞::Float64
    t_star::Float64
    σ_c::Float64
    m::Float64

    # precomputed elastic tangent tensor and hardening matrix
    Eᵉ::SymmetricTensor{4,3,Float64,36}
    H::Matrix{Float64}
    MS::SVector{S,Tensor{2,3,Float64,9}}
    # store information about slip sytems
    nslipsystems::Int
    slipsystems::Vector{NTuple{2, Vec{3}}}

    function CrystalViscoPlastic(E, ν, τ_y, H_iso, H_kin, q, α_∞, t_star, σ_c, m, slipsystems)
        Eᵉ = elastic_tangent_3D(E, ν)
        nslipsystems = length(slipsystems)
        H = H_iso*(q*ones(nslipsystems, nslipsystems) + (1-q)*diagm(ones(nslipsystems)))
        MS = SVector{nslipsystems,Tensor{2,3,Float64,9}}([slipsystems[i][1] ⊗ slipsystems[i][2] for i=1:nslipsystems])
        return new{nslipsystems}(E, ν, τ_y, H_iso, H_kin, q, α_∞, t_star, σ_c, m, Eᵉ, H, MS, nslipsystems, slipsystems)
    end
end

# keyword argument constructor
CrystalViscoPlastic(; E, ν, τ_y, H_iso, H_kin, q, α_∞, t_star, σ_c, m, slipsystems) = CrystalViscoPlastic(E, ν, τ_y, H_iso, H_kin, q, α_∞, t_star, σ_c, m, slipsystems)

get_n_slipsystems(::CrystalViscoPlastic{S}) where S = S

struct CrystalViscoPlasticState{dim,T,M,S} <: AbstractMaterialState
    σ::SymmetricTensor{2,dim,T,M}
    κ::SVector{S, Float64}
    α::SVector{S, Float64}
    μ::SVector{S, Float64}
end

Base.zero(::Type{CrystalViscoPlasticState{dim,T,M,S}}) where {dim,T,M,S} = CrystalViscoPlasticState{dim,T,M,S}(zero(SymmetricTensor{2,dim,T,M}), SVector{S,T}(zeros(T, S)), SVector{S,T}(zeros(T, S)), SVector{S,T}(zeros(T, S)))
initial_material_state(::CrystalViscoPlastic{S}) where S = zero(CrystalViscoPlasticState{3,Float64,6,S})

function get_cache(m::CrystalViscoPlastic)
    state = initial_material_state(m)
    # it doesn't actually matter which state and strain step we use here,
    # f is overwritten in the constitutive driver before being used.
    f(r_vector, x_vector) = vector_residual!(((x)->MaterialModels.residuals(x, m, state, zero(SymmetricTensor{2,3}), 1.0)), r_vector, x_vector, m)
    v_cache = Vector{Float64}(undef, get_n_scalar_equations(m))
    cache = NLsolve.OnceDifferentiable(f, v_cache, v_cache; autodiff = :forward)
    return cache
end

# constitutive driver operates in 3D, so these can always be 3D
# TODO: Should Residuals have a Type Parameter N for the number of scalar equations?
struct ResidualsCrystalViscoPlastic{S, T}
    σ::SymmetricTensor{2,3,T,6}
    κ::SVector{S, T}
    α::SVector{S, T}
    μ::SVector{S, T}
    # function ResidualsCrystalViscoPlastic(σ::SymmetricTensor{2,3,T,6}, κ::MVector{S, T}, α::MVector{S, T}, μ::MVector{S, T}) where {S, T}
    #     # S = length(κ)
    #     # length(α) == S || error("$S slipsystems for κ, but $(length(α)) slipsystems for α handed over to Residuals{CrystalViscoPlastic}.")
    #     # length(μ) == S || error("$S slipsystems for κ, but $(length(μ)) slipsystems for μ handed over to Residuals{CrystalViscoPlastic}.")
    #     new{S,T}(σ, κ, α, μ)
    # end
end

get_n_scalar_equations(::CrystalViscoPlastic{S}) where S = 6 + 3S
# get_residual_type(::CrystalViscoPlastic{S}) where S = ResidualsCrystalViscoPlastic{S}
Tensors.get_base(::Type{CrystalViscoPlastic{S}}) where S = ResidualsCrystalViscoPlastic{S} # needed for frommandel

function Tensors.tomandel!(v::Vector{T}, r::ResidualsCrystalViscoPlastic{S, T}) where {S, T}
    M=6
    # TODO check vector length
    tomandel!(view(v, 1:M), r.σ)
    v[M+1:M+S] = view(r.κ, :)
    v[M+S+1:M+2S] = view(r.α, :)
    v[M+2S+1:M+3S] = view(r.μ, :)
    return v
end

function Tensors.frommandel(::Type{ResidualsCrystalViscoPlastic{S}}, v::Vector{T}) where {S, T}
    M = 6
    σ = frommandel(SymmetricTensor{2,3}, view(v, 1:M))
    κ = SVector{S,T}(view(v, M+1:M+S))
    α = SVector{S,T}(view(v, M+S+1:M+2S))
    μ = SVector{S,T}(view(v, M+2S+1:M+3S))
    return ResidualsCrystalViscoPlastic(σ, κ, α, μ)
end

function material_response(m::CrystalViscoPlastic{S}, Δε::SymmetricTensor{2,3,T,6}, state::CrystalViscoPlasticState{3},
    Δt=1.0; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {S, T}

    # set the current residual function that depends only on the variables
    f(r_vector, x_vector) = vector_residual!((x->residuals(x,m,state,Δε,Δt)), r_vector, x_vector, m)
    update_cache!(cache, f)
    # initial guess
    σ_trial = state.σ + m.Eᵉ ⊡ Δε
    x0 = ResidualsCrystalViscoPlastic(σ_trial, state.κ, state.α, state.μ)
    # convert initial guess to vector
    tomandel!(cache.x_f, x0)
    # solve for variables x
    nlsolve_options = get(options, :nlsolve_params, Dict{Symbol, Any}(:method=>:newton))
    haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method
    result = NLsolve.nlsolve(cache, cache.x_f; nlsolve_options...)

    if result.f_converged
        x_mandel = result.zero::Vector{T}
        x = frommandel(ResidualsCrystalViscoPlastic{S}, x_mandel)
        dRdx = cache.DF
        inv_J_σσ = frommandel(SymmetricTensor{4,3}, inv(dRdx))
        ∂σ∂ε = inv_J_σσ ⊡ m.Eᵉ
        return x.σ, ∂σ∂ε, CrystalViscoPlasticState{3,T,6,S}(x.σ, SVector{S,T}(x.κ), SVector{S,T}(x.α), SVector{S,T}(x.μ))
    else
        error("Material model not converged. Could not find material state.")
    end
end

function residuals(vars::ResidualsCrystalViscoPlastic{S,T}, mat::CrystalViscoPlastic{S}, material_state::CrystalViscoPlasticState{3}, Δε, Δt) where {S,T}
    σ=vars.σ; κ=vars.κ; α=vars.α; μ=vars.μ

    temp_sum_σ = zero(Tensor{2,3})
    Rκ = MVector{S, T}(undef)
    Rα = MVector{S, T}(undef)
    Rμ = MVector{S, T}(undef)
    for i=1:S
        # m = mat.slipsystems[i][1]
        # s = mat.slipsystems[i][2]
        # ms_tensor = m ⊗ s
        ms_tensor = mat.MS[i]
        τ = σ ⊡ ms_tensor
        σ_red = τ - α[i]
        stress_sign = sign(σ_red)
        Φ = norm(σ_red) - (mat.τ_y + κ[i])

        # R_σ
        temp_sum_σ += μ[i] * ms_tensor * stress_sign

        # R_κ
        # temp_sum_κ = 0.0
        # for j=1:S
        #     # H is symmetric, thus we can take columns instead of rows
        #     temp_sum_κ += mat.H[j,i]*μ[i]
        # end
        Rκ[i] = κ[i] - material_state.κ[i] - dot(view(mat.H, :,i), μ)

        # R_α
        Rα[i] = α[i] - material_state.α[i] - μ[i]*mat.H_kin*(stress_sign-α[i]/mat.α_∞)

        # R_Φ
        Rμ[i] = Δt*(((Φ+abs(Φ))/2)/mat.σ_c)^mat.m - mat.t_star*μ[i]
    end

    Rσ = σ - material_state.σ - mat.Eᵉ ⊡ Δε + mat.Eᵉ ⊡ temp_sum_σ
    return  ResidualsCrystalViscoPlastic(Rσ, SVector{S,T}(Rκ), SVector{S,T}(Rα), SVector{S,T}(Rμ))
end

# # use StaticVectors everywhere
# function residuals_var(vars::ResidualsCrystalViscoPlastic{S,T}, mat::CrystalViscoPlastic{S}, material_state::CrystalViscoPlasticState{3}, Δε, Δt) where {S,T}
#     σ=vars.σ; κ=vars.κ; α=vars.α; μ=vars.μ

#     temp_sum_σ = zero(Tensor{2,3})
#     for i=1:S
#         m = mat.slipsystems[i][1]
#         s = mat.slipsystems[i][2]
#         ms_tensor = m ⊗ s
#         τ = σ ⊡ ms_tensor

#         stress_sign = sign(τ - α[i])
#         Φ = norm(τ - α[i]) - (mat.τ_y + κ[i])

#         # R_σ
#         temp_sum_σ += μ[i] * ms_tensor * stress_sign

#         # R_κ
#         temp_sum_κ = 0.0
#         for j=1:S
#             # H is symmetric, thus we can take columns instead of rows
#             temp_sum_κ += mat.H[j,i]*μ[i]
#         end
#         Rκ[i] = κ[i] - material_state.κ[i] - temp_sum_κ

#         # R_α
#         Rα[i] = α[i] - material_state.α[i] - μ[i]*mat.H_kin*(stress_sign-α[i]/mat.α_∞)

#         # R_Φ
#         Rμ[i] = Δt*(((Φ+abs(Φ))/2)/mat.σ_c)^mat.m - mat.t_star*μ[i]
#     end

#     Rσ = σ - material_state.σ - mat.Eᵉ ⊡ Δε + mat.Eᵉ ⊡ temp_sum_σ
#     return  ResidualsCrystalViscoPlastic(Rσ, Rκ, Rα, Rμ)
# end


text:
#model url: https://github.com/FluxML/Metalhead.jl/releases/download/v0.1.1/resnet.bson
#reference: https://github.com/FluxML/Metalhead.jl
using Knet, KnetLayers, BSON, Images

struct ResidualBlock
  layers
  shortcut
end

function ResidualBlock(filters, kernels::Array{Tuple{Int,Int}}, pads::Array{Tuple{Int,Int}}, strides::Array{Tuple{Int,Int}}, shortcut = identity)
  layers = []
  for i in 2:length(filters)
    push!(layers, Conv(activation=nothing, height=kernels[i-1][1], width=kernels[i-1][2], inout=filters[i-1]=>filters[i], padding = pads[i-1], stride = strides[i-1], mode=1, binit=nothing))
    if i != length(filters)
      push!(layers, Chain(BatchNorm(filters[i]),ReLU())) # I think we need batchnorm with relu activation
    else
      push!(layers, BatchNorm(filters[i]))
    end
  end
  ResidualBlock(Chain(layers...), shortcut)
end

ResidualBlock(filters, kernels::Array{Int}, pads::Array{Int}, strides::Array{Int}, shortcut = identity) =
  ResidualBlock(filters, [(i,i) for i in kernels], [(i,i) for i in pads], [(i,i) for i in strides], shortcut)

(r::ResidualBlock)(input) = relu.(r.layers(input) + r.shortcut(input))

function BasicBlock(filters::Int, downsample::Bool = false, res_top::Bool = false)
  # NOTE: res_top is set to true if this is the first residual connection of the architecture
  # If the number of channels is to be halved set the downsample argument to true
  if !downsample || res_top
    return ResidualBlock([filters for i in 1:3], [3,3], [1,1], [1,1])
  end
  shortcut = Chain(Conv(activation=nothing, height=3, width=3, inout=filters÷2=>filters, padding = (1,1), stride = (2,2), mode=1, binit=nothing), BatchNorm(filters))
  ResidualBlock([filters÷2, filters, filters], [3,3], [1,1], [1,2], shortcut)
end

function Bottleneck(filters::Int, downsample::Bool = false, res_top::Bool = false)
  # NOTE: res_top is set to true if this is the first residual connection of the architecture
  # If the number of channels is to be halved set the downsample argument to true
  if !downsample && !res_top
    ResidualBlock([4 * filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,1])
  elseif downsample && res_top
    ResidualBlock([filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,1], Chain(Conv(activation=nothing, height=1,width=1, inout=filters=>4 * filters, padding = (0,0), stride = (1,1), mode=1, binit=nothing), BatchNorm(4 * filters)))
  else
    shortcut = Chain(Conv(activation=nothing, height=1, width=1, inout=2 * filters=>4 * filters, padding = (0,0), stride = (2,2), mode=1, binit=nothing), BatchNorm(4 * filters))
    ResidualBlock([2 * filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,2], shortcut)
  end
end

const KA = KnetLayers.arrtype
transfer!(p::Param, x)  = copyto!(p.value,x)
transfer!(p, x) = copyto!(p,x)
to4D(x) = reshape(x,1,1,length(x),1)

function trained_resnet50_layers(model=KnetLayers.dir("examples","resnet.bson"))
    if !isfile(model)
        download("https://github.com/FluxML/Metalhead.jl/releases/download/v0.1.1/resnet.bson",model)
    end
    weight = BSON.load(model)
    weights = Dict{Any ,Any}()
    for ele in keys(weight)
        weights[string(ele)] = weight[ele]
    end
    ls = load_resnet(resnet_configs["resnet50"]...)
    transfer!(ls[1][1].weight, weights["gpu_0/conv1_w_0"])
    ls[1][2].moments.var =  KA(weights["gpu_0/res_conv1_bn_riv_0"] |> to4D)
    ls[1][2].moments.mean = KA(weights["gpu_0/res_conv1_bn_rm_0"] |> to4D)
    ls[1][2].params = KA(vcat(weights["gpu_0/res_conv1_bn_s_0"],weights["gpu_0/res_conv1_bn_b_0"]))
    count = 2
    for j in [3:5, 6:9, 10:15, 16:18]
        for p in j
            transfer!(ls[p].layers[1].weight, weights["gpu_0/res$(count)_$(p-j[1])_branch2a_w_0"] )
            ls[p].layers[2][1].moments.var = KA(weights["gpu_0/res$(count)_$(p-j[1])_branch2a_bn_riv_0"] |> to4D)
            ls[p].layers[2][1].moments.mean = KA(weights["gpu_0/res$(count)_$(p-j[1])_branch2a_bn_rm_0"] |> to4D)
            ls[p].layers[2][1].params =  KA(vcat(weights["gpu_0/res$(count)_$(p-j[1])_branch2a_bn_s_0"],weights["gpu_0/res$(count)_$(p-j[1])_branch2a_bn_b_0"]))
            transfer!(ls[p].layers[3].weight , weights["gpu_0/res$(count)_$(p-j[1])_branch2b_w_0"])
            ls[p].layers[4][1].moments.var = KA(weights["gpu_0/res$(count)_$(p-j[1])_branch2b_bn_riv_0"] |> to4D)
            ls[p].layers[4][1].moments.mean = KA(weights["gpu_0/res$(count)_$(p-j[1])_branch2b_bn_rm_0"] |> to4D)
            ls[p].layers[4][1].params =  KA(vcat(weights["gpu_0/res$(count)_$(p-j[1])_branch2b_bn_s_0"],weights["gpu_0/res$(count)_$(p-j[1])_branch2b_bn_b_0"]))
            transfer!(ls[p].layers[5].weight , weights["gpu_0/res$(count)_$(p-j[1])_branch2c_w_0"])
            ls[p].layers[6].moments.var = KA(weights["gpu_0/res$(count)_$(p-j[1])_branch2c_bn_riv_0"] |> to4D)
            ls[p].layers[6].moments.mean = KA(weights["gpu_0/res$(count)_$(p-j[1])_branch2c_bn_rm_0"] |> to4D)
            ls[p].layers[6].params =  KA(vcat(weights["gpu_0/res$(count)_$(p-j[1])_branch2c_bn_s_0"],weights["gpu_0/res$(count)_$(p-j[1])_branch2c_bn_b_0"]))
        end
        transfer!(ls[j[1]].shortcut[1].weight , weights["gpu_0/res$(count)_0_branch1_w_0"])
        ls[j[1]].shortcut[2].moments.var = KA(weights["gpu_0/res$(count)_0_branch1_bn_riv_0"] |> to4D)
        ls[j[1]].shortcut[2].moments.mean = KA(weights["gpu_0/res$(count)_0_branch1_bn_rm_0"] |> to4D)
        ls[j[1]].shortcut[2].params =  KA(vcat(weights["gpu_0/res$(count)_0_branch1_bn_s_0"], weights["gpu_0/res$(count)_0_branch1_bn_b_0"]))
        count += 1
    end
    transfer!(ls[21].mult.weight, permutedims(weights["gpu_0/pred_w_0"], (2,1)));
    transfer!(ls[21].bias, weights["gpu_0/pred_b_0"]);
    return ls
end

function load_resnet(Block, layers, initial_filters::Int = 64, nclasses::Int = 1000)
  local top = []
  local residual = []
  local bottom = []

  push!(top, Chain(Conv(activation=nothing, width=7, height=7, inout=3=>initial_filters, padding = (3,3), stride = (2,2), mode=1, binit=nothing), BatchNorm(initial_filters)))
  push!(top, Pool(window=(3,3), padding = (1,1), stride = (2,2)))

  for i in 1:length(layers)
    push!(residual, Block(initial_filters, true, i==1))
    for j in 2:layers[i]
      push!(residual, Block(initial_filters))
    end
    initial_filters *= 2
  end

  push!(bottom, Pool(window=(7,7), mode=1))
  push!(bottom, x -> mat(x))
  if Block == Bottleneck
    push!(bottom, (Linear(input=2048,output=nclasses)))
  else
    push!(bottom, (Dense(input=512,output=nclasses)))
  end
  push!(bottom, softmax)

  Chain(top..., residual..., bottom...)
end

resnet_configs =
  Dict("resnet18" => (BasicBlock, [2, 2, 2, 2]),
       "resnet34" => (BasicBlock, [3, 4, 6, 3]),
       "resnet50" => (Bottleneck, [3, 4, 6, 3]),
       "resnet101" => (Bottleneck, [3, 4, 23, 3]),
       "resnet152" => (Bottleneck, [3, 8, 36, 3]))

struct ResNet18
  layers::Chain
end

ResNet18() = ResNet18(load_resnet(resnet_configs["resnet18"]...))

trained(::Type{ResNet18}) = error("Pretrained Weights for ResNet18 are not available")

Base.show(io::IO, ::ResNet18) = print(io, "ResNet18()")

(m::ResNet18)(x) = m.layers(x)

struct ResNet34
  layers::Chain
end

ResNet34() = ResNet34(load_resnet(resnet_configs["resnet34"]...))

trained(::Type{ResNet34}) = error("Pretrained Weights for ResNet34 are not available")

Base.show(io::IO, ::ResNet34) = print(io, "ResNet34()")

(m::ResNet34)(x) = m.layers(x)

struct ResNet50
  layers::Chain
end

ResNet50() = ResNet50(load_resnet(resnet_configs["resnet50"]...))

trained(::Type{ResNet50}) = ResNet50(trained_resnet50_layers())

Base.show(io::IO, ::ResNet50) = print(io, "ResNet50()")

(m::ResNet50)(x) = m.layers(x)

struct ResNet101
  layers::Chain
end

ResNet101() = ResNet101(load_resnet(resnet_configs["resnet101"]...))

trained(::Type{ResNet101}) = error("Pretrained Weights for ResNet101 are not available")

Base.show(io::IO, ::ResNet101) = print(io, "ResNet101()")

(m::ResNet101)(x) = m.layers(x)

struct ResNet152
  layers::Chain
end

ResNet152() = ResNet152(load_resnet(resnet_configs["resnet152"]...))

trained(::Type{ResNet152}) = error("Pretrained Weights for ResNet152 are not available")

Base.show(io::IO, ::ResNet152) = print(io, "ResNet152()")

(m::ResNet152)(x) = m.layers(x)

function preprocess(img::AbstractMatrix{<:AbstractRGB})
  # Resize such that smallest edge is 256 pixels long
    img = resize_smallest_dimension(img, 256)
    im = center_crop(img, 224)
    z = (channelview(im) .* 255 .- 128)./128;
    z1 = Float32.(permutedims(z, (3, 2, 1))[:,:,:,:]);
end

# Resize an image such that its smallest dimension is the given length
function resize_smallest_dimension(im, len)
  reduction_factor = len/minimum(size(im)[1:2])
  new_size = size(im)
  new_size = (
      round(Int, size(im,1)*reduction_factor),
      round(Int, size(im,2)*reduction_factor),
  )
  if reduction_factor < 1.0
    # Images.jl's imresize() needs to first lowpass the image, it won't do it for us
    im = imfilter(im, KernelFactors.gaussian(0.75/reduction_factor), Inner())
  end
  return imresize(im, new_size)
end

# Take the len-by-len square of pixels at the center of image `im`
function center_crop(im, len)
  l2 = div(len,2)
  adjust = len % 2 == 0 ? 1 : 0
  return im[div(end,2)-l2:div(end,2)+l2-adjust,div(end,2)-l2:div(end,2)+l2-adjust]
end


text:
module LinearLeastSquares

using Compat
import Compat.String

include("types/expressions.jl")
include("types/constraints.jl")
include("types/problems.jl")

include("utils/utilities.jl")
include("utils/display.jl")

include("atoms/repmat.jl")
include("atoms/reshape.jl")
include("atoms/add_subtract.jl")
include("atoms/multiply_divide.jl")
include("atoms/transpose.jl")
include("atoms/index.jl")
include("atoms/stack.jl")
include("atoms/diag.jl")
include("atoms/sum.jl")
include("atoms/mean.jl")
include("atoms/variance.jl")

include("solvers/backslash_solver.jl")

include("solve.jl")
end # module


text:
"""
`module DB`
## Overview
This module implements a "database" for storing a precomputed LSQ system. This
is useful for fitting e.g. NBodyIPs since it allows one to precompute the <basis,
data> inner products which are very expensive, and then quickly construct LSQ
systems from them using e.g., many variants of weighting and regularisation.
A db called (e.g.) 'lsqdb' consists of two files:
 * `lsqdb_info.json` : this stores a list of basis functions and a list
 of configurations with attached observations (normally DFT data)
 * `lsqdb_kron.h5` : this stores all the inner products <basis, data>
 as a single large matrix. The observations and basis functions in
 `lsqdb_info` store the corresponding column and row indices.
## More details and remarks
### INFO
* BASIS: The basis is (for the time being) simply represented as a list (Vector) of
JuLIP.AbstractCalculator. These are stored as `INFO["basis"]`.
* DATA: Each "datum" is an atomistic configuration, which must have a
"configtype", stored as a `Dat`.
Informally two pieces of data with the same configtype should represent
"similar" kinds of configurations. `data = INFO["data"]` is a Dictionary
where the keys are the configtypes.
### (De-)Vectorising Data
Forces in `JuLIP` are represented as `Vector{SVector{...}}`, which is the
same memory layout as a 3 x N matrix (`JuLIP.mat`), which can can then be
vectorised (`[:]`), and this vectoriation is readily undone again.
Analogously, any data must be stored in such a vectorised format. This is
achieved (e.g. for forces) via
* `vec_obs(::Val{:F}, F) -> Vector`
* `devec_obs(::Val{:F}, Fvec) -> Vector{SVector{...}}`
or equivalently
* `vec_obs("F", F) -> Vector`
* `devec_obs("F", Fvec) -> Vector{SVector{...}}`
(the `Val` versions are used for performance optimisation)
"""
module DB

using Base.Threads:          SpinLock, nthreads
using StaticArrays:          SVector
using JuLIP:                 AbstractCalculator, AbstractAtoms, Atoms, energy, forces,
                             save_dict, load_dict, read_dict, write_dict
using JuLIP.MLIPs:           IPBasis
using IPFitting:        Dat, LsqDB, basis, eval_obs, observations,
                             observation, vec_obs, devec_obs,
                             tfor_observations
using IPFitting.Data:   configtype
using HDF5:                  h5open, read

using DistributedArrays
using Distributed

import Base: flush, append!, union

export LsqDB, LsqDB_dist, info, configtypes

const VERSION = 2
const KRONFILE = "_kron.h5"
const INFOFILE = "_info.json"

"""
`dbpath(db::LsqDB)` : return the absolute path to the database files, not
including the 'kron.h5' and 'info.json' endings.
"""
dbpath(db::LsqDB) = db.dbpath

infofile(dbpath::AbstractString) = dbpath * INFOFILE

kronfile(dbpath::AbstractString) = dbpath * KRONFILE

# ------------ Save and load the info file

"if the file exists, append a random string to avoid overwriting"
function _backupfile(fname)
   if isfile(fname)
      fnew = fname * "." * String(rand('a':'z', 5))
      @warn("The file $fname already exists. It will be renamed to $fnew to avoid overwriting.")
      run(`mv $fname $fnew`)
   end
   return nothing
end

function load_info(dbpath::String)
   dbinfo = load_dict(infofile(dbpath))
   if !haskey(dbinfo, "version")
      @error("This LsqDB was produced with an older version of IPFitting.")
   end
   version = dbinfo["version"]
   if version < VERSION
      @error("This LsqDB was produced with an older version of IPFitting.")
   end
   basis = read_dict(dbinfo["basis"])
   configs = Dat.(dbinfo["configs"])   # here we already know the type
   return basis, configs
end

function save_info(dbpath::String, db)
   _backupfile(infofile(dbpath))
   save_dict(infofile(dbpath),
             Dict("version" => VERSION,
                  "basis" => write_dict(db.basis),
                  "configs" => write_dict.(db.configs))
            )
   return nothing
end

save_info(db) = save_info(dbpath(db), db)

# ----------- Save and load the KRON file

"save a single matrix to HDF5"
_savemath5(A, fname) =
   h5open(fname, "w") do fid
      fid["A"] = A
      nothing
   end

"load a single matrix from HDF5"
_loadmath5(fname) =
   h5open(fname, "r") do fid
      read(fid["A"])
   end

function save_kron(dbpath, db)
   _backupfile(kronfile(dbpath))
   _savemath5(db.Ψ, kronfile(dbpath))
end

save_kron(db) = save_kron(dbpath(db), db)

load_kron(dbpath::String; mmap=false) = _loadmath5(kronfile(dbpath))

load_kron(db::LsqDB; mmap=false) = load_kron(dbpath(db); mmap=mmap)

function LsqDB(dbpath::AbstractString; mmap=true)
   basis, configs = load_info(dbpath)
   Ψ = load_kron(dbpath; mmap=mmap)
   return LsqDB(basis, configs, Ψ, dbpath)
end

"""
for the time being, this just checks that the db director and db name are
admissible. In the future, this should initialise the DB files.
"""
function initdb(basedir, dbname)
   @assert !('/' in dbname)
   @assert isdir(basedir)
   dbpath = joinpath(basedir, dbname)
   initdb(dbpath)
   return nothing
end

function initdb(dbpath)
   # TODO: seems this fails to detect an existing database?
   @assert !isfile(infofile(dbpath))
   @assert !isfile(kronfile(dbpath))
   # check that we can actually create and delete this file
   save_dict(infofile(dbpath), Dict("version" => VERSION, "basis" => [], "data" => []))
   rm(infofile(dbpath))
   _savemath5(rand(5,5), kronfile(dbpath))
   rm(kronfile(dbpath))
   return nothing
end

function flush(db::LsqDB)
   save_info(db)
   save_kron(db)
   return nothing
end

function LsqDB(dbpath::AbstractString,
               basis::IPBasis,
               configs::AbstractVector{Dat};
               verbose=true,
               maxnthreads=nthreads())
   # assign indices, count observations and allocate a matrix
   Ψ = _alloc_lsq_matrix(configs, basis)
   # create the struct where everything is stored
   db = LsqDB(basis, configs, Ψ, dbpath)
   # parallel assembly of the LSQ matrix
   tfor_observations( configs,
      (n, okey, cfg, lck) -> safe_append!(db, lck, cfg, okey),
      msg = "Assemble LSQ blocks",
      verbose=verbose,
      maxnthreads=maxnthreads )
   # save to file
   if dbpath != ""
      verbose && @info("Writing db to disk...")
      try
         flush(db)
      catch
         @warn("""something went wrong trying to save the db to disk, but the data
               should be ok; if it is crucial to keep it, try to save manually.""")
      end
      verbose && @info("... done")
   else
      verbose && @info("db is not written to disk since `dbpath` is empty.")
   end
   return db
end

# function LsqDB_dist(basis::IPBasis,
#                      configs::AbstractVector{Dat})
#       # assign indices, count observations and allocate a matrix
#       Ψ = _alloc_lsq_matrix_dist(configs, basis, 4)
#       irow = 0
#       @sync @distributed for i = 1:nworkers()
#          # add energy to the lsq system
#          irow += 1
#          #y[irow] = wE * E / length(at)
#          #Ψ_part = localpart(Ψ)
#          #@show size(Ψ_part)
#          #@show Ψ[irow, :] .= repeat([irow], 36)
#          Ψ_part = localpart(Ψ)
#          vec(Ψ_part) .= (1:length(Ψ_part)) .+ 1000*myid()
#          #Ψ_part .= energy(basis, at.at) / length(at.at)

#          # add forces to the lsq system
#          # nf = 3*length(at.at)
#          # #y[(irow+1):(irow+nf)] = wF * mat(F)[:]
#          # Fb = forces(basis, at.at)
#          # for ib = 1:length(basis)
#          #    Ψ[(irow+1):(irow+nf), ib] = mat(Fb[ib])[:]
#          # end
#          # irow += nf
#       end
#       return Ψ
# end

function set_matrows!(d::Dat, okey::String, irows::Vector{Int})
   d.rows[okey] = irows
   return d
end

matrows(d::Dat, okey::String) = d.rows[okey]

function _alloc_lsq_matrix(configs, basis)
   # the indices associated with the basis are simply the indices within
   # the array - there is nothing else to do here
   #
   # loop through all observations and assign indices
   nrows = 0
   for (okey, d, _) in observations(configs)
      len = length(observation(d, okey))
      set_matrows!(d, okey, collect(nrows .+ (1:len)))
      nrows += len
   end
   # allocate and return the matrix
   return zeros(Float64, nrows, length(basis))
end

function _alloc_lsq_matrix_dist(configs, basis, nprocs)
   # the indices associated with the basis are simply the indices within
   # the array - there is nothing else to do here
   #
   # loop through all observations and assign indices
   nrows = 0
   for (okey, d, _) in observations(configs)
      len = length(observation(d, okey))
      set_matrows!(d, okey, collect(nrows .+ (1:len)))
      nrows += len
   end
   # allocate and return the matrix
   return dzeros((nrows,length(basis)), workers()[1:nprocs], [1,nprocs])
end

function safe_append!(db::LsqDB, db_lock, cfg, okey)
   # computing the lsq blocks ("rows") can be done in parallel,
   lsqrow = eval_obs(okey, basis(db), cfg)
   vec_lsqrow = vec_obs(okey, lsqrow)
   irows = matrows(cfg, okey)
   # but writing them into the DB must be done in a threadsafe way
   lock(db_lock)
   db.Ψ[irows, :] = vec_lsqrow
   unlock(db_lock)
   return nothing
end

# ------------------- Evaluating LSQ Blocks -----------------

# TODO: SOMETHING LIKE THIS IS STILL NEEDED...
# """
# evaluate one specific kind of datum, such as energy, forces, etc
# for one Dat across a large section of the basis; implicitly this will
# normally be done via `evaluate_many` or similar. The line
# ```
# vals = eval_obs(vDT, B, at)
# ```
# should likely be the bottleneck of `evallsq`
# """
# function _evallsq(vDT::Val,
#                   B::AbstractVector{<:AbstractCalculator},
#                   at::AbstractAtoms)
#    # vals will be a vector containing multiple evaluations
#    vals = eval_obs(vDT, B, at)
#    # vectorise the first so we know the length of the observation
#    vec1 = vec_obs(vDT, vals[1])
#    # create a multi-D array to reshape these into
#    A = Array{Float64}(undef, length(vec1), length(B))
#    A[:, 1] = vec1
#    for n = 2:length(B)
#       A[:, n] = vec_obs(vDT, vals[n])
#    end
#    return A
# end

# =========================== Convenience =================================

configtypes(db::LsqDB) = unique(configtype.(db.configs))

_nconfigs(db::LsqDB, ct::AbstractString) =
   length(find(configtype.(db.configs) .== ct))

# TODO: rewrite this!
function info(db::LsqDB)
   # config names, how many
   all_cts = configtype.(db.configs)
   unq_cts = unique(all_cts)
   configs_info = Dict{String, Int}()
   for ct in unq_cts
      configs_info[ct] = sum(all_cts .== ct)
   end

   println("======================================================")
   println("       LsqDB Summary ")
   println("------------------------------------------------------")
   println(" Datagroup: configname  =>  number of configs ")
   for (i, (key, val)) in enumerate(configs_info)
      println("        $i : $key  =>  $val ")
   end
   # basis groups, how many, max degree
   B = db.basis
   Bord, Iord = split_basis(B)
   println("------------------------------------------------------")
   println(" Basis Group  | length | description ")
   for (i, B1) in enumerate(Bord)
      lenstr = replace(string(length(B1), pad=5), '0' => ' ')
      idxstr = replace(string(i, pad=2), '0' => ' ')
      desc = string(B1)
      desc = desc[1:min(50, length(desc))]
      println("          $idxstr  | $lenstr  | $desc" )
      #
   end
   println("======================================================")
end

end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct RenewableFix <: RenewableGen
        name::String
        available::Bool
        bus::Bus
        activepower::Float64
        reactivepower::Float64
        rating::Float64
        primemover::PrimeMovers.PrimeMover
        powerfactor::Float64
        basepower::Float64
        services::Vector{Service}
        dynamic_injector::Union{Nothing, DynamicInjection}
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

Data Structure for fixed renewable generation technologies.

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `activepower::Float64`
- `reactivepower::Float64`
- `rating::Float64`: Thermal limited MVA Power Output of the unit. <= Capacity, validation range: (0, nothing), action if invalid: error
- `primemover::PrimeMovers.PrimeMover`: PrimeMover Technology according to EIA 923
- `powerfactor::Float64`, validation range: (0, 1), action if invalid: error
- `basepower::Float64`: Base power of the unit in system base per unit, validation range: (0, nothing), action if invalid: warn
- `services::Vector{Service}`: Services that this device contributes to
- `dynamic_injector::Union{Nothing, DynamicInjection}`: corresponding dynamic injection device
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct RenewableFix <: RenewableGen
    name::String
    available::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    "Thermal limited MVA Power Output of the unit. <= Capacity"
    rating::Float64
    "PrimeMover Technology according to EIA 923"
    primemover::PrimeMovers.PrimeMover
    powerfactor::Float64
    "Base power of the unit in system base per unit"
    basepower::Float64
    "Services that this device contributes to"
    services::Vector{Service}
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, DynamicInjection}
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function RenewableFix(name, available, bus, activepower, reactivepower, rating, primemover, powerfactor, basepower, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    RenewableFix(name, available, bus, activepower, reactivepower, rating, primemover, powerfactor, basepower, services, dynamic_injector, ext, forecasts, InfrastructureSystemsInternal(), )
end

function RenewableFix(; name, available, bus, activepower, reactivepower, rating, primemover, powerfactor, basepower, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    RenewableFix(name, available, bus, activepower, reactivepower, rating, primemover, powerfactor, basepower, services, dynamic_injector, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function RenewableFix(::Nothing)
    RenewableFix(;
        name="init",
        available=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        rating=0.0,
        primemover=PrimeMovers.OT,
        powerfactor=1.0,
        basepower=0.0,
        services=Device[],
        dynamic_injector=nothing,
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::RenewableFix) = value.name
"""Get RenewableFix available."""
get_available(value::RenewableFix) = value.available
"""Get RenewableFix bus."""
get_bus(value::RenewableFix) = value.bus
"""Get RenewableFix activepower."""
get_activepower(value::RenewableFix) = value.activepower
"""Get RenewableFix reactivepower."""
get_reactivepower(value::RenewableFix) = value.reactivepower
"""Get RenewableFix rating."""
get_rating(value::RenewableFix) = value.rating
"""Get RenewableFix primemover."""
get_primemover(value::RenewableFix) = value.primemover
"""Get RenewableFix powerfactor."""
get_powerfactor(value::RenewableFix) = value.powerfactor
"""Get RenewableFix basepower."""
get_basepower(value::RenewableFix) = value.basepower
"""Get RenewableFix services."""
get_services(value::RenewableFix) = value.services
"""Get RenewableFix dynamic_injector."""
get_dynamic_injector(value::RenewableFix) = value.dynamic_injector
"""Get RenewableFix ext."""
get_ext(value::RenewableFix) = value.ext

InfrastructureSystems.get_forecasts(value::RenewableFix) = value.forecasts
"""Get RenewableFix internal."""
get_internal(value::RenewableFix) = value.internal

InfrastructureSystems.set_name!(value::RenewableFix, val::String) = value.name = val
"""Set RenewableFix available."""
set_available!(value::RenewableFix, val::Bool) = value.available = val
"""Set RenewableFix bus."""
set_bus!(value::RenewableFix, val::Bus) = value.bus = val
"""Set RenewableFix activepower."""
set_activepower!(value::RenewableFix, val::Float64) = value.activepower = val
"""Set RenewableFix reactivepower."""
set_reactivepower!(value::RenewableFix, val::Float64) = value.reactivepower = val
"""Set RenewableFix rating."""
set_rating!(value::RenewableFix, val::Float64) = value.rating = val
"""Set RenewableFix primemover."""
set_primemover!(value::RenewableFix, val::PrimeMovers.PrimeMover) = value.primemover = val
"""Set RenewableFix powerfactor."""
set_powerfactor!(value::RenewableFix, val::Float64) = value.powerfactor = val
"""Set RenewableFix basepower."""
set_basepower!(value::RenewableFix, val::Float64) = value.basepower = val
"""Set RenewableFix services."""
set_services!(value::RenewableFix, val::Vector{Service}) = value.services = val
"""Set RenewableFix ext."""
set_ext!(value::RenewableFix, val::Dict{String, Any}) = value.ext = val

InfrastructureSystems.set_forecasts!(value::RenewableFix, val::InfrastructureSystems.Forecasts) = value.forecasts = val
"""Set RenewableFix internal."""
set_internal!(value::RenewableFix, val::InfrastructureSystemsInternal) = value.internal = val


text:
"""
    parity_sign(P)

Given an array or tuple of integers representing
a permutation or a subset of a permutation,
compute the parity sign defined as -1 for a
permutation consisting of an odd number of swaps
and +1 for an even number of swaps. This
implementation uses an O(n^2) algorithm and is
intended for small permutations only.
"""
function parity_sign(P)::Int
  L = length(P)
  s = +1
  for i=1:L, j=i+1:L
    s *= sign(P[j]-P[i])
  end
  return s
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"""SelfSubjectRulesReview enumerates the set of actions the current user can perform within a namespace. The returned list of actions may be incomplete depending on the server's authorization mode, and any errors experienced during the evaluation. SelfSubjectRulesReview should be used by UIs to show/hide actions, or to quickly let an end user reason about their permissions. It should NOT Be used by external systems to drive authorization decisions as this raises confused deputy, cache lifetime/revocation, and correctness concerns. SubjectAccessReview, and LocalAccessReview are the correct way to defer authorization decisions to the API server.

    IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview(;
        apiVersion=nothing,
        kind=nothing,
        metadata=nothing,
        spec=nothing,
        status=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta
    - spec::IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec : Spec holds information about the request being evaluated.
    - status::IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus : Status is filled in by the server and indicates the set of actions a user can perform.
"""
mutable struct IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus } # spec name: status

    function IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview

const _property_map_IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec", Symbol("status")=>"IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus")
Base.propertynames(::Type{ IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview }) = collect(keys(_property_map_IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview))
Swagger.property_type(::Type{ IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview }, property_name::Symbol) =  _property_map_IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview[property_name]

function check_required(o::IoK8sApiAuthorizationV1beta1SelfSubjectRulesReview)
    (getproperty(o, Symbol("spec")) === nothing) && (return false)
    true
end

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


text:
export ForwardEulerODEFunction

"""
    ForwardEulerODEFunction(f)

An ODE function wrapper where `f(un, u, p, t, dt)` provides a forward Euler update
```
un .= u .+ dt * f(u, p, t)
```

"""
struct ForwardEulerODEFunction{F} <: DiffEqBase.AbstractODEFunction{true}
    f::F
end
(f::ForwardEulerODEFunction{F})(un, u, p, t, dt) where {F} = f.f(un, u, p, t, dt)

"""
    OffsetODEFunction(f,α,β,γ,x)

An ODE function wrapper which evaluates `f` with an offset.

Evaluates as
```math
f(u,p,α+β*t) .+ γ .* x
```

It supports the 3, 4, 5, and 6 argument forms.
"""
mutable struct OffsetODEFunction{iip,F,S,A} <: DiffEqBase.AbstractODEFunction{iip}
    f::F
    α::S
    β::S
    γ::S
    x::A
end
function OffsetODEFunction{iip}(f,α,β,γ,x) where {iip}
    α,β,γ = promote(α,β,γ)
    OffsetODEFunction{iip, typeof(f), typeof(γ), typeof(x)}(f,α,β,γ,x)
end
OffsetODEFunction(f::DiffEqBase.AbstractODEFunction{iip},α,β,γ,x) where {iip} =
    OffsetODEFunction{iip}(f,α,β,γ,x)

function (o::OffsetODEFunction)(u,p,t)
    o.f(u,p,o.α+o.β*t) .+ o.γ .* o.x
end
function (o::OffsetODEFunction)(du,u,p,t)
    o.f(du,u,p,o.α+o.β*t)
    du .+= o.γ .* o.x
end
function (o::OffsetODEFunction)(du,u,p,t,α)
    o.f(du,u,p,o.α+o.β*t,α)
    du .+= α .* o.γ .* o.x
end
function (o::OffsetODEFunction)(du,u,p,t,α,β)
    o.f(du,u,p,o.α+o.β*t,α,β)
    du .+= α .* o.γ .* o.x
end


text:
#=
Define a transition model for the car MDP problem
=#

function POMDPs.transition(mdp::CarMDP, s::CarMDPState, a::CarMDPAction)
    # ego transition
    acc = LonAccelDirection(a.acc, 4) #TODO paremeterize
    ego_p = propagate(s.ego, acc, mdp.env.roadway, mdp.ΔT)
    # interpolate ego_p in discrete space
    ego_v_space = get_car_vspace(mdp.env, mdp.vel_res)
    ego_ps, ego_probs = interpolate_state(mdp, ego_p, ego_v_space)
    # car transition
    car_probs = Float64[]
    car_routes = Vector{Lane{Float64}}[]
    if s.car.posG == mdp.off_grid
        car_ps, car_routes = car_starting_states(mdp)
        push!(car_ps, get_off_the_grid(mdp))
        push!(car_routes, SVector{2, LaneTag}(LaneTag(0,0), LaneTag(0,0)))
        car_probs = ones(length(car_ps))
        car_probs[end] = 1-mdp.car_birth
        car_probs[1:end-1] = mdp.car_birth/(length(car_ps)-1)
        @assert sum(car_probs) ≈ 1.0
    else
        # update model
        mdp.car_model.route = [mdp.env.roadway[tag] for tag in find_route(mdp.env, s.route)] # s.route is a static vector but not car_model.route
        lane = get_lane(mdp.env.roadway, s.car)
        set_direction!(mdp.car_model, lane, mdp.env.roadway)
        car_actions, car_actions_probs = get_distribution(mdp.car_model)
        car_ps = VehicleState[]
        for (i, car_a) in enumerate(car_actions)
            p_a = car_actions_probs[i]
            car_p = propagate(s.car, car_a, mdp.env.roadway, mdp.ΔT)
            if car_p.posF.s >= get_end(get_lane(mdp.env.roadway, car_p)) && isempty(get_lane(mdp.env.roadway, car_p).exits)
                car_pp = get_off_the_grid(mdp)
                index_car_pp = findall(x -> x==car_pp, car_ps)
                if isempty(index_car_pp)
                    push!(car_ps, car_pp)
                    push!(car_probs, p_a)
                else
                    car_probs[index_car_pp] += p_a
                end
                car_routes = fill(s.route, length(car_ps))
                continue
            end
            # interpolate car_p in discrete space
            car_v_space = get_car_vspace(mdp.env, mdp.vel_res)
            itp_car_ps, itp_car_weights = interpolate_state(mdp, car_p, car_v_space)
            for (j, car_pss) in enumerate(itp_car_ps)
                index_itp_state = findall(x -> x==car_pss, car_ps)
                if isempty(index_itp_state)
                    push!(car_ps, car_pss)
                    push!(car_probs, itp_car_weights[j]*p_a)
                else
                    car_probs[index_itp_state] += itp_car_weights[j]*p_a
                end
            end
            car_routes = fill(s.route, length(car_ps))
        end
        @assert length(car_probs) == length(car_ps)
        @assert length(car_routes) == length(car_ps)
        normalize!(car_probs, 1)
    end

    # combine car and ego states
    states_p = Vector{CarMDPState}(length(ego_ps)*length(car_ps)) # future states
    probs_p = Vector{Float64}(length(ego_ps)*length(car_ps)) # proba of future states
    k = 1
    for (i, e) in enumerate(ego_ps)
        for (j, p) in enumerate(car_ps)
            crash =  is_colliding(Vehicle(e, mdp.ego_type, 1), Vehicle(p, mdp.car_type, 2))
            states_p[k] = CarMDPState(crash, e, p, SVector(car_routes[j]...))
            probs_p[k] = ego_probs[i]*car_probs[j]
            k += 1
        end
    end
    normalize!(probs_p, 1)
    return SparseCat(states_p, probs_p)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
function FieldHansenZonalLinear(arg0::jint, arg1::jint, arg2::Field)
    return FieldHansenZonalLinear((jint, jint, Field), arg0, arg1, arg2)
end

function compute_init_values(obj::FieldHansenZonalLinear, arg0::RealFieldElement)
    return jcall(obj, "computeInitValues", void, (RealFieldElement,), arg0)
end

function get_derivative(obj::FieldHansenZonalLinear, arg0::jint, arg1::RealFieldElement)
    return jcall(obj, "getDerivative", RealFieldElement, (jint, RealFieldElement), arg0, arg1)
end

function get_value(obj::FieldHansenZonalLinear, arg0::jint, arg1::RealFieldElement)
    return jcall(obj, "getValue", RealFieldElement, (jint, RealFieldElement), arg0, arg1)
end


text:
function mysum_eachindex(a::AbstractArray)
s = zero(eltype(a))
for i in eachindex(a)
s = s + a[i]
end
end

max_stars_repo_path: Module 2/Chapter06/mysum_index.jl

max_issues_repo_path: Module 2/Chapter06/mysum_index.jl

max_forks_repo_path: Module 2/Chapter06/mysum_index.jl


text:
module JacobiDavidson

using LinearMaps
using LinearAlgebra
using IterativeSolvers
using Random

import LinearAlgebra: ldiv!
import LinearAlgebra.BLAS: axpy!, gemv!

include("subspaces.jl")
include("correction_eqn_solvers.jl")
include("orthogonalization.jl")
include("eigenvalue_targets.jl")
include("jdqr.jl")
include("jdqz.jl")

end

max_forks_repo_name: haampie/jacobi-davidson-julia


text:
@doc """
    gen_snapshot!(model, T, [N=1])

generate and return `N` snapshots (spin configuration).

# Arguments
* `model::Model` : model
* `T::Real` : temperature
* `N::Integer=1` : the number of snapshots to be generated
* `MCS::Integer=8192` : the number of Monte Carlo steps followed by snapshot generation
"""
function gen_snapshot!(model::Model, T::Real; MCS::Integer=8192)
    for m in 1:MCS
        SW_update!(model, T)
    end
    return copy(model.spins)
end

function gen_snapshot!(model::Model, T::Real, N::Integer; MCS::Integer=8192)
    nsites = numsites(model.lat)
    X = zeros(nsites, N)
    for n in 1:N
        spins[:,n] = gen_snapshot!(model, T, MCS=MCS)
    end
    return X
end

@doc """
    gensave_snapshot!(io, model, T, [N=1])
    gensave_snapshot!(filename, model, T, [N=1])

generate and write `N` snapshots into `io` or `filename`.

# Arguments
* `io::IO` : output stream where snapshots will be written
* `filename::String` : the name of file where snapshots will be written
* `model::Model` : model to be simulated
* `T::Real` : temperature
* `N::Integer=1` : the number of snapshots to be generated
* `MCS::Integer=1` : the number of Monte Carlo steps followed by snapshot generation
* `sep::String=" "` : field separator of snapshot
"""
function gensave_snapshot!(io::IO, model::Model, T::Real, N::Integer=1; MCS::Integer=8192, sep::String=" ")
    nsites = numsites(model.lat)
    for n in 1:N
        for m in 1:MCS
            SW_update!(model, T)
        end
        spins = model.spins
        print(io, spins[1])
        for s in 2:nsites
            print(io, sep, spins[s])
        end
        println(io)
    end
end

function gensave_snapshot!(filename::String, model::Model, T::Real, N::Integer=1; MCS::Integer=8192, sep::String=" ", append::Bool=false)
    c = ifelse(append, "a", "w")
    open(filename, c) do io
        gensave_snapshot!(io, model, T, N, MCS=MCS, sep=sep)
    end
end

function load_snapshot_impl(io::IO, splitter)
    for line in eachline(io)
        if ismatch(r"^\s*($|#)", line)
            continue
        end
        produce(float(split(strip(line), splitter)))
    end
end
function load_snapshot_impl(filename::String, splitter)
    io = open(filename)
    for line in eachline(io)
        if ismatch(r"^\s*($|#)", line)
            continue
        end
        produce(float(split(strip(line), splitter)))
    end
    close(io)
end
const _default_delims = [' ','\t','\n','\v','\f','\r']

@doc """
    load_snapshot(io, splitter)
    load_snapshot(filename, splitter)

load snapshot from `io` or `filename`.
`splitter` is the set of field splitter and will be passed `Base.split` as the second argument (see the doc of `Base.split`).
"""
load_snapshot(io::IO, splitter = _default_delims) = @task load_snapshot_impl(io, splitter)
load_snapshot(filename::String, splitter = _default_delims) = @task load_snapshot_impl(filename, splitter)


text:
module MazeUndirectedModule

import ..GridWorlds as GW
import Random
import ReinforcementLearningBase as RLBase

#####
##### game logic
#####

const NUM_OBJECTS = 3
const AGENT = 1
const WALL = 2
const GOAL = 3
const NUM_ACTIONS = 4

mutable struct MazeUndirected{R, RNG} <: GW.AbstractGridWorld
    tile_map::BitArray{3}
    agent_position::CartesianIndex{2}
    reward::R
    rng::RNG
    done::Bool
    terminal_reward::R
    goal_position::CartesianIndex{2}
end

function MazeUndirected(; R = Float32, height = 9, width = 9, rng = Random.GLOBAL_RNG)
    @assert isodd(height) && isodd(width) "height and width must be odd numbers"

    tile_map = falses(NUM_OBJECTS, height, width)

    tile_map[WALL, :, :] .= true
    for j in 2 : 2 : width - 1
        for i in 2 : 2 : height - 1
            tile_map[WALL, i, j] = false
        end
    end

    agent_position = CartesianIndex(2, 2)
    tile_map[AGENT, agent_position] = true

    goal_position = CartesianIndex(height - 1, width - 1)
    tile_map[GOAL, goal_position] = true

    reward = zero(R)
    done = false
    terminal_reward = one(R)

    env = MazeUndirected(tile_map, agent_position, reward, rng, done, terminal_reward, goal_position)

    GW.reset!(env)

    return env
end

function GW.reset!(env::MazeUndirected)
    tile_map = env.tile_map
    rng = env.rng
    _, height, width = size(tile_map)

    tile_map[AGENT, env.agent_position] = false
    tile_map[GOAL, env.goal_position] = false

    tile_map[WALL, :, :] .= true
    for j in 2 : 2 : width - 1
        for i in 2 : 2 : height - 1
            tile_map[WALL, i, j] = false
        end
    end

    generate_maze!(env)

    new_agent_position = GW.sample_empty_position(rng, tile_map)
    env.agent_position = new_agent_position
    tile_map[AGENT, new_agent_position] = true

    new_goal_position = GW.sample_empty_position(rng, tile_map)
    env.goal_position = new_goal_position
    tile_map[GOAL, new_goal_position] = true

    env.reward = zero(env.reward)
    env.done = false

    return nothing
end

function GW.act!(env::MazeUndirected, action)
    @assert action in Base.OneTo(NUM_ACTIONS) "Invalid action $(action). Action must be in Base.OneTo($(NUM_ACTIONS))"

    tile_map = env.tile_map
    agent_position = env.agent_position

    if action == 1
        new_agent_position = GW.move_up(agent_position)
    elseif action == 2
        new_agent_position = GW.move_down(agent_position)
    elseif action == 3
        new_agent_position = GW.move_left(agent_position)
    else
        new_agent_position = GW.move_right(agent_position)
    end

    if !tile_map[WALL, new_agent_position]
        tile_map[AGENT, agent_position] = false
        env.agent_position = new_agent_position
        tile_map[AGENT, new_agent_position] = true
    end

    if tile_map[GOAL, env.agent_position]
        env.reward = env.terminal_reward
        env.done = true
    else
        env.reward = zero(env.reward)
        env.done = false
    end

    return nothing
end

function generate_maze!(env::MazeUndirected)
    tile_map = env.tile_map
    rng = env.rng
    _, height, width = size(tile_map)

    vertical_range = 2:2:height-1
    horizontal_range = 2:2:width-1

    visited = falses(height, width)
    stack = CartesianIndex{2}[]

    first_cell = CartesianIndex(rand(rng, vertical_range), rand(rng, horizontal_range))
    visited[first_cell] = true
    push!(stack, first_cell)

    valid_neighbors = CartesianIndex{2}[]

    while !isempty(stack)
        current_cell = pop!(stack)

        i, j = current_cell.I
        all_neighbors = (CartesianIndex(i - 2, j), CartesianIndex(i, j - 2), CartesianIndex(i, j + 2), CartesianIndex(i + 2, j))

        empty!(valid_neighbors)
        for neighbor in all_neighbors
            if (neighbor.I[1] in vertical_range) && (neighbor.I[2] in horizontal_range) && !visited[neighbor]
                push!(valid_neighbors, neighbor)
            end
        end

        if length(valid_neighbors) > 0
            push!(stack, current_cell)

            next_cell = rand(rng, valid_neighbors)

            mid_pos = CartesianIndex((current_cell.I .+ next_cell.I) .÷ 2)
            tile_map[WALL, mid_pos] = false

            visited[next_cell] = true
            push!(stack, next_cell)
        end
    end

    return nothing
end

#####
##### miscellaneous
#####

GW.get_height(env::MazeUndirected) = size(env.tile_map, 2)
GW.get_width(env::MazeUndirected) = size(env.tile_map, 3)

GW.get_action_names(env::MazeUndirected) = (:MOVE_UP, :MOVE_DOWN, :MOVE_LEFT, :MOVE_RIGHT)
GW.get_object_names(env::MazeUndirected) = (:AGENT, :WALL, :GOAL)

function GW.get_pretty_tile_map(env::MazeUndirected, position::CartesianIndex{2})
    characters = ('☻', '█', '♥', '⋅')

    object = findfirst(@view env.tile_map[:, position])
    if isnothing(object)
        return characters[end]
    else
        return characters[object]
    end
end

function GW.get_pretty_sub_tile_map(env::MazeUndirected, window_size, position::CartesianIndex{2})
    tile_map = env.tile_map
    agent_position = env.agent_position

    characters = ('☻', '█', '♥', '⋅')

    sub_tile_map = GW.get_sub_tile_map(tile_map, agent_position, window_size)

    object = findfirst(@view sub_tile_map[:, position])
    if isnothing(object)
        return characters[end]
    else
        return characters[object]
    end
end

function Base.show(io::IO, ::MIME"text/plain", env::MazeUndirected)
    str = "tile_map:\n"
    str = str * GW.get_pretty_tile_map(env)
    str = str * "\nsub_tile_map:\n"
    str = str * GW.get_pretty_sub_tile_map(env, GW.get_window_size(env))
    str = str * "\nreward: $(env.reward)"
    str = str * "\ndone: $(env.done)"
    str = str * "\naction_names: $(GW.get_action_names(env))"
    str = str * "\nobject_names: $(GW.get_object_names(env))"
    print(io, str)
    return nothing
end

GW.get_action_keys(env::MazeUndirected) = ('w', 's', 'a', 'd')

#####
##### RLBase API
#####

RLBase.StateStyle(env::GW.RLBaseEnv{E}) where {E <: MazeUndirected} = RLBase.InternalState{Any}()
RLBase.state_space(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: MazeUndirected} = nothing
RLBase.state(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: MazeUndirected} = env.env.tile_map

RLBase.reset!(env::GW.RLBaseEnv{E}) where {E <: MazeUndirected} = GW.reset!(env.env)

RLBase.action_space(env::GW.RLBaseEnv{E}) where {E <: MazeUndirected} = Base.OneTo(NUM_ACTIONS)
(env::GW.RLBaseEnv{E})(action) where {E <: MazeUndirected} = GW.act!(env.env, action)

RLBase.reward(env::GW.RLBaseEnv{E}) where {E <: MazeUndirected} = env.env.reward
RLBase.is_terminated(env::GW.RLBaseEnv{E}) where {E <: MazeUndirected} = env.env.done

end # module


text:
using Pkg, DrWatson
#@quickactivate "StatisticalRethinkingStan"
using StanSample
using XGBoost
using StatisticalRethinking

include("$(@__DIR__)/rstar.jl")
include("$(@__DIR__)/clip-04-31s.jl");

chns = read_samples(sm, :mcmcchains)

rstar(chns, 1)

rs = rstar(chns, 100)

mean(rs) |> display

histogram(rs)

gui()

closeall()


text:
struct DataFrame
    stream_identifier::UInt32
    is_end_stream::Bool
    data::Array{UInt8, 1}
end

==(a::DataFrame, b::DataFrame) =
    a.stream_identifier == b.stream_identifier &&
    a.is_end_stream == b.is_end_stream &&
    a.data == b.data

function decode_data(header, payload)
    is_end_stream = header.flags & 0x1 == 0x1

    payload = remove_padding(header, payload)

    return DataFrame(header.stream_identifier, is_end_stream, payload)
end

function encode_data(frame)
    typ = DATA
    flags = 0x0 | (frame.is_end_stream ? 0x1 : 0x0)
    stream_identifier = frame.stream_identifier

    return wrap_payload(frame.data, typ, flags, frame.stream_identifier)
end


text:
using Flux, HDF5, Random

include("BayesianSZNet.jl")
using .BayesianSZNet

h5open("data/dr16q_superset.hdf5", "r+") do datafile
    X = gpu(read(datafile, "X"))
    model = SZNet("models/sznet.bson")
    Random.seed!(45)
    ẑs = sample(model, X)
    write_dataset(datafile, "zs_pred", ẑs)
end

max_stars_repo_name: podondra/bdl-redshift

max_issues_repo_name: podondra/bdl-redshift

max_forks_repo_name: podondra/bdl-redshift


text:
using NewtonMethods
include("ChebyPoly.jl")

fSol = open("nleq1.dat","w")
fRest = open("nleq1.out","w")
dimMax = 9
for dim = 2:dimMax
    n1 = dim + 1

    # Initialize the options
    opt = OptionsNLEQ(OPT_MODE              => 1,
                      OPT_JACGEN            => 1,
                      OPT_JACFCN            => chebyQuadJac,
                      OPT_MSTOR             => 0,
                      OPT_NOROWSCAL         => 0,
                      OPT_PRINTWARNING      => 1,
                      OPT_PRINTITERATION    => 3,
                      OPT_PRINTSOLUTION     => 2,
                      OPT_PRINTIOWARN       => fRest,
                      OPT_PRINTIOMON        => fRest,
                      OPT_PRINTIOSOL        => fSol,
                      OPT_NITMAX            => 10,
                      OPT_RTOL              => 1e-5)

    x0    = collect(1:dim)./n1
    xScal = zeros(x0)

    retCode = -1

    i = 1

    println("Calling the while loop for solving the Cheby Quad equation of dimension $dim")
    while retCode == -1
        (x0, _, retCode) = nleq1(chebyQuad, x0, xScal, opt)
        write(fRest, @sprintf("Returned from call %4i of NLEQ1\n",i))
		flush(fSol)
	    flush(fRest)
        i += 1
    end
    println("Solution = $x0")

    println("DONE")
end
close(fSol)
close(fRest)


text:
using ProfileView

function profile_test(n)
    for i = 1:n
        A = randn(100,100,20)
        m = maximum(A)
        Am = mapslices(sum, A, dims = 2)
        B = A[:,:,5]
        Bsort = mapslices(sort, B, dims = 1)
        b = rand(100)
        C = B.*b
    end
end

function profile_simple_test(n)
    for i = 1:n
        A = randn(100,100,20)
        m = max(A)
    end
end

unstable(x) = x > 0.5 ? true : 0.0

function profile_unstable_test(m, n)
    s = s2 = 0
    for i = 1:n
        for k = 1:m
            s += unstable(rand())
        end
        x = collect(1:20)
        s2 += sum(x)
    end
    s, s2
end

profile_test(1)
@profview profile_test(10)

using Profile

Profile.clear()
@profile profile_test(10)
ProfileView.view()
ProfileView.view(C=true)
ProfileView.view(fontsize=18)

Profile.clear()
@profile profile_unstable_test(1, 1)
ProfileView.view(ProfileView.FlameGraphs.StackFrameCategory())

Profile.clear()
@profile profile_unstable_test(10, 10^6)
ProfileView.view()

ProfileView.closeall()


text:
using jInv.Mesh;
using Helmholtz
using Multigrid
using LinearAlgebra
using SparseArrays

const plotting = false;

if plotting
	using PyPlot;
	close("all");
end

# m = readdlm("SEGmodel2Dsalt.dat"); m = Matrix(m'); m = m*1e-3;
m = 1.5*ones(257,129);

Minv = getRegularMesh([0.0,13.5,0.0,4.2],collect(size(m)).-1);

pad = 16;

m = 1.0./m.^2

f = 2.5;

w = 2*pi*f

println("omega*h:");
println(w*Minv.h*sqrt(maximum(m)));
pad = pad*ones(Int64,Minv.dim);

maxOmega = getMaximalFrequency(m,Minv);
ABLamp = maxOmega;

Sommerfeld = true;
NeumannAtFirstDim = true;
H,gamma = GetHelmholtzOperator(Minv,m,w,w*ones(size(m))*0.01,NeumannAtFirstDim,pad,ABLamp,Sommerfeld);
shift = [0.2;0.2;0.2];
SH = H .+ GetHelmholtzShiftOP(m, real(w),shift[1]);

# n = Minv.n .+1 ; n_tup = tuple(n...);
# src = div.(n,2);
# src[end] = 1;
# q = zeros(ComplexF64,n_tup)
# q[loc2cs(n,src)] = 1.0/(Minv.h[1]^2);

levels      = 3;
numCores 	= 2;
maxIter     = 30;
relativeTol = 1e-6;
# relaxType   = "Jac-GMRES";
relaxType   = "Jac";
relaxParam  = 0.75;
relaxPre 	= 2;
relaxPost   = 2;
cycleType   ='W';
coarseSolveType = "NoMUMPS";

MG = getMGparam(ComplexF64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,
				relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0)

## Preparing a point shource RHS ############################
n = Minv.n .+ 1; n_tup = tuple(n...);
src = div.(n,2);
src[end] = 1;
q = zeros(ComplexF64,n_tup)
q[loc2cs(n,src)] = 1.0/(Minv.h[1]^2);
###############################################
b = q[:];

Hparam = HelmholtzParam(Minv,gamma,vec(m),w,NeumannAtFirstDim,Sommerfeld);
Ainv = getShiftedLaplacianMultigridSolver(Hparam, MG,shift,"GMRES",5,true);
Ainv = copySolver(Ainv);

x = solveLinearSystem(sparse(SH'),b,Ainv)[1];
MG.relaxType = "Jac";
MG.cycleType = 'W';
Ainv = getShiftedLaplacianMultigridSolver(Hparam, MG,shift,"BiCGSTAB",0,true);
x = solveLinearSystem(sparse(SH'),b,Ainv)[1];

if plotting
	reX = real(reshape(x,n_tup));
	reX[My_sub2ind(n,src)] = 0.0;
	println(norm(q[:]-H*x)/norm(q[:]));
	println(norm(q[:]-H*y)/norm(q[:]));
	figure();
	imshow(reX');title("Helmholtz Iterative Solution");
end

if plotting
	s = H\q[:];
	s = real(reshape(s,n_tup));
	s[My_sub2ind(n,src)] = 0.0;
	figure();
	imshow(s');title("Helmholtz True Solution");
end

# println("Doing Transpose")
# y = solveLinearSystem(sparse(SH'),b,Ainv,1)[1];

# println(norm(H'*y-b)/norm(b));

################## WITHOUT THE COARSEST GRID SOL ########################
# coarseSolveType = "DDNodal";
# MG = getMGparam(ComplexF64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,
				# relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0)
# Ainv = getShiftedLaplacianMultigridSolver(Hparam, MG,shift,"BiCGSTAB",0,true);

# y = solveLinearSystem(SH',b,Ainv)[1];
###########################################################################

println("SOLVING MULTIPLE RHSs")
coarseSolveType = "Julia";
MG = getMGparam(ComplexF64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,
				relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0)
Ainv = getShiftedLaplacianMultigridSolver(Hparam, MG,shift,"BiCGSTAB",0,true);
nrhs = 2;
b = rand(ComplexF64,length(b),nrhs);
clear!(Ainv);
Ainv.helmParam = Hparam;
x = solveLinearSystem(sparse(SH'),b,Ainv)[1];
println(norm(H*x .- b)/norm(b))

MG = getMGparam(ComplexF64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType);
MG.relaxType = "Jac-GMRES";
MG.cycleType = 'K';
Ainv = getShiftedLaplacianMultigridSolver(Hparam, MG,shift,"GMRES",5,true);
x = solveLinearSystem(sparse(SH'),b,Ainv)[1];
println(norm(H*x .- b)/norm(b))


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# tests that interpreter matches codegen
using Test
using Core: GotoIfNot, ReturnNode

# test that interpreter correctly handles PhiNodes (#29262)
let m = Meta.@lower 1 + 1
    @assert Meta.isexpr(m, :thunk)
    src = m.args[1]::Core.CodeInfo
    src.code = Any[
        # block 1
        QuoteNode(:a),
        QuoteNode(:b),
        GlobalRef(@__MODULE__, :test29262),
        GotoIfNot(Core.SSAValue(3), 6),
        # block 2
        Core.PhiNode(Any[4], Any[Core.SSAValue(1)]),
        Core.PhiNode(Any[4, 5], Any[Core.SSAValue(2), Core.SSAValue(5)]),
        ReturnNode(Core.SSAValue(6)),
    ]
    nstmts = length(src.code)
    src.ssavaluetypes = Any[ Any for _ = 1:nstmts ]
    src.codelocs = fill(Int32(1), nstmts)
    src.inferred = true
    Core.Compiler.verify_ir(Core.Compiler.inflate_ir(src))
    global test29262 = true
    @test :a === @eval $m
    global test29262 = false
    @test :b === @eval $m
end

let m = Meta.@lower 1 + 1
    @assert Meta.isexpr(m, :thunk)
    src = m.args[1]::Core.CodeInfo
    src.code = Any[
        # block 1
        QuoteNode(:a),
        QuoteNode(:b),
        QuoteNode(:c),
        GlobalRef(@__MODULE__, :test29262),
        # block 2
        Core.PhiNode(Any[4, 16], Any[false, true]), # false, true
        Core.PhiNode(Any[4, 16], Any[Core.SSAValue(1), Core.SSAValue(2)]), # :a, :b
        Core.PhiNode(Any[4, 16], Any[Core.SSAValue(3), Core.SSAValue(6)]), # :c, :a
        Core.PhiNode(Any[16], Any[Core.SSAValue(7)]), # NULL, :c
        # block 3
        Core.PhiNode(Any[], Any[]), # NULL, NULL
        Core.PhiNode(Any[17, 8], Any[true, Core.SSAValue(4)]), # test29262, test29262, [true]
        Core.PhiNode(Any[17], Vector{Any}(undef, 1)), # NULL, NULL
        Core.PhiNode(Any[8], Vector{Any}(undef, 1)), # NULL, NULL
        Core.PhiNode(Any[], Any[]), # NULL, NULL
        Core.PhiNode(Any[17, 8], Any[Core.SSAValue(2), Core.SSAValue(8)]), # NULL, :c, [:b]
        Core.PhiNode(Any[], Any[]), # NULL, NULL
        GotoIfNot(Core.SSAValue(5), 5),
        # block 4
        GotoIfNot(Core.SSAValue(10), 9),
        # block 5
        Expr(:call, GlobalRef(Core, :tuple), Core.SSAValue(6), Core.SSAValue(7), Core.SSAValue(8), Core.SSAValue(14)),
        ReturnNode(Core.SSAValue(18)),
    ]
    nstmts = length(src.code)
    src.ssavaluetypes = Any[ Any for _ = 1:nstmts ]
    src.codelocs = fill(Int32(1), nstmts)
    src.inferred = true
    Core.Compiler.verify_ir(Core.Compiler.inflate_ir(src))
    global test29262 = true
    @test (:b, :a, :c, :c) === @eval $m
    global test29262 = false
    @test (:b, :a, :c, :b) === @eval $m
end

let m = Meta.@lower 1 + 1
    @assert Meta.isexpr(m, :thunk)
    src = m.args[1]::Core.CodeInfo
    src.code = Any[
        # block 1
        QuoteNode(:a),
        QuoteNode(:b),
        GlobalRef(@__MODULE__, :test29262),
        # block 2
        Expr(:enter, 11),
        # block 3
        Core.UpsilonNode(),
        Core.UpsilonNode(),
        Core.UpsilonNode(Core.SSAValue(2)),
        GotoIfNot(Core.SSAValue(3), 10),
        # block 4
        Core.UpsilonNode(Core.SSAValue(1)),
        # block 5
        Expr(:throw_undef_if_not, :expected, false),
        # block 6
        Core.PhiCNode(Any[Core.SSAValue(5), Core.SSAValue(7), Core.SSAValue(9)]), # NULL, :a, :b
        Core.PhiCNode(Any[Core.SSAValue(6)]), # NULL
        Expr(:leave, 1),
        # block 7
        ReturnNode(Core.SSAValue(11)),
    ]
    nstmts = length(src.code)
    src.ssavaluetypes = Any[ Any for _ = 1:nstmts ]
    src.codelocs = fill(Int32(1), nstmts)
    src.inferred = true
    Core.Compiler.verify_ir(Core.Compiler.inflate_ir(src))
    global test29262 = true
    @test :a === @eval $m
    global test29262 = false
    @test :b === @eval $m
end


text:
import Libdl

const DEFAULT_PATH_URL = "http://pages.cs.wisc.edu/~ferris/path/julia/"

function _prefix_suffix(str)
    if Sys.islinux() &&  Sys.ARCH == :x86_64
        return "lib$(str).so"
    elseif Sys.isapple()
        return "lib$(str).dylib"
    elseif Sys.iswindows()
        return "$(str).dll"
    end
    error(
        "Unsupported operating system. Only 64-bit linux, OSX, and Windows " *
        "are supported."
    )
end

function download_path()
    libpath = _prefix_suffix("path50")
    ENV["PATH_JL_LOCATION"] = joinpath(@__DIR__, libpath)
    download(joinpath(DEFAULT_PATH_URL, libpath), ENV["PATH_JL_LOCATION"])
    return
end

function install_path()
    if !haskey(ENV, "PATH_JL_LOCATION")
        download_path()
    end
    lusol = _prefix_suffix("lusol")
    lusol_filename = joinpath(@__DIR__, lusol)
    download(joinpath(DEFAULT_PATH_URL, lusol), lusol_filename)
    local_filename = get(ENV, "PATH_JL_LOCATION", nothing)
    if local_filename === nothing
        error("Environment variable `PATH_JL_LOCATION` not found.")
    elseif Libdl.dlopen(local_filename) == C_NULL
        error("Unable to open the path library $(local_filename).")
    end
    open("deps.jl", "w") do io
        write(io, "const PATH_SOLVER = \"$(escape_string(local_filename))\"\n")
        write(io, "const LUSOL_LIBRARY_PATH = \"$(escape_string(lusol_filename))\"\n")
    end
end

install_path()