text:
CategoricalValue(pool::CategoricalPool{T, R}, level::Integer) where {T, R} =
    CategoricalValue(pool, convert(R, level))

"""
    CategoricalValue(value, source::Union{CategoricalValue, CategoricalArray})

Return a `CategoricalValue` object wrapping `value` and attached to
the [`CategoricalPool`](@ref) of `source`.
"""
function CategoricalValue(value, source::Union{CategoricalValue, CatArrOrSub})
    p = pool(source)
    i = get(p, value, nothing)
    i === nothing && throw(ArgumentError("level $value not found in source pool"))
    return CategoricalValue(p, i)
end

leveltype(::Type{<:CategoricalValue{T}}) where {T} = T
leveltype(::Type{T}) where {T} = T
leveltype(x::Any) = leveltype(typeof(x))
# to fix ambiguity
leveltype(x::Type{Union{}}) = Union{}

# integer type of category reference codes used by categorical value
reftype(::Type{<:CategoricalValue{<:Any, R}}) where {R} = R
reftype(x::Any) = reftype(typeof(x))

pool(x::CategoricalValue) = x.pool
refcode(x::CategoricalValue) = x.ref
isordered(x::CategoricalValue) = isordered(x.pool)

# extract the type of the original value from array eltype `T`
unwrap_catvaluetype(::Type{T}) where {T} = T
unwrap_catvaluetype(::Type{T}) where {T >: Missing} =
    Union{unwrap_catvaluetype(nonmissingtype(T)), Missing}
unwrap_catvaluetype(::Type{Union{}}) = Union{} # prevent incorrect dispatch to T<:CategoricalValue method
unwrap_catvaluetype(::Type{Any}) = Any # prevent recursion in T>:Missing method
unwrap_catvaluetype(::Type{T}) where {T <: CategoricalValue} = leveltype(T)

"""
    unwrap(x::CategoricalValue)
    unwrap(x::Missing)

Get the value wrapped by categorical value `x`. If `x` is `Missing` return `missing`.
"""
DataAPI.unwrap(x::CategoricalValue) = levels(x)[refcode(x)]

"""
    levelcode(x::CategoricalValue)

Get the code of categorical value `x`, i.e. its index in the set
of possible values returned by [`levels(x)`](@ref DataAPI.levels).
"""
levelcode(x::CategoricalValue) = Signed(widen(refcode(x)))

"""
    levelcode(x::Missing)

Return `missing`.
"""
levelcode(x::Missing) = missing

DataAPI.levels(x::CategoricalValue) = levels(pool(x))

function cat_promote_type(::Type{S}, ::Type{T}) where {S, T}
    U = promote_type(S, T)
    U <: Union{SupportedTypes, Missing} ?
        U : typeintersect(Union{SupportedTypes, Missing}, Union{S, T})
end

cat_promote_eltype() = Union{}
cat_promote_eltype(v1, vs...) = cat_promote_type(eltype(v1), cat_promote_eltype(vs...))

Base.promote_rule(::Type{C}, ::Type{T}) where {C <: CategoricalValue, T} =
    promote_type(leveltype(C), T)

Base.promote_rule(::Type{C}, ::Type{T}) where {C <: CategoricalValue, T >: Missing} =
    Union{promote_rule(C, nonmissingtype(T)), Missing}

# To fix ambiguities with definitions from Base
Base.promote_rule(::Type{C}, ::Type{Missing}) where {C <: CategoricalValue} = Union{C, Missing}
Base.promote_rule(::Type{C}, ::Type{Any}) where {C <: CategoricalValue} = Any

Base.promote_rule(::Type{C1}, ::Type{C2}) where
    {R1<:Integer, R2<:Integer,
     C1<:CategoricalValue{<: SupportedTypes, R1},
     C2<:CategoricalValue{<: SupportedTypes, R2}} =
    CategoricalValue{cat_promote_type(leveltype(C1), leveltype(C2)), promote_type(R1, R2)}
Base.promote_rule(::Type{C1}, ::Type{C2}) where {C1<:CategoricalValue, C2<:CategoricalValue} =
    CategoricalValue{cat_promote_type(leveltype(C1), leveltype(C2))}

# General fallbacks
Base.convert(::Type{S}, x::CategoricalValue) where {S <: SupportedTypes} =
    convert(S, unwrap(x))
Base.convert(::Type{Union{S, Missing}}, x::CategoricalValue) where {S <: SupportedTypes} =
    convert(Union{S, Missing}, unwrap(x))
Base.convert(::Type{Union{S, Nothing}}, x::CategoricalValue) where {S <: SupportedTypes} =
    convert(Union{S, Nothing}, unwrap(x))

(::Type{T})(x::T) where {T <: CategoricalValue} = x

Base.Broadcast.broadcastable(x::CategoricalValue) = Ref(x)

function Base.show(io::IO, x::CategoricalValue)
    if nonmissingtype(get(io, :typeinfo, Any)) === nonmissingtype(typeof(x))
        show(io, unwrap(x))
    else
        print(io, typeof(x))
        print(io, ' ')
        show(io, unwrap(x))
        if isordered(pool(x))
            @printf(io, " (%i/%i)", levelcode(x), length(pool(x)))
        end
    end
end

Base.print(io::IO, x::CategoricalValue) = print(io, unwrap(x))
Base.string(x::CategoricalValue) = string(unwrap(x))
Base.write(io::IO, x::CategoricalValue) = write(io, unwrap(x))
Base.String(x::CategoricalValue{<:AbstractString}) = String(unwrap(x))

@inline function Base.:(==)(x::CategoricalValue, y::CategoricalValue)
    if pool(x) === pool(y) || pool(x) == pool(y)
        return refcode(x) == refcode(y)
    else
        return unwrap(x) == unwrap(y)
    end
end

Base.:(==)(x::CategoricalValue, y::SupportedTypes) = unwrap(x) == y
Base.:(==)(x::SupportedTypes, y::CategoricalValue) = x == unwrap(y)

@inline function Base.isequal(x::CategoricalValue, y::CategoricalValue)
    if pool(x) === pool(y) || pool(x) == pool(y)
        return refcode(x) == refcode(y)
    else
        return isequal(unwrap(x), unwrap(y))
    end
end

Base.isequal(x::CategoricalValue, y::SupportedTypes) = isequal(unwrap(x), y)
Base.isequal(x::SupportedTypes, y::CategoricalValue) = isequal(x, unwrap(y))

Base.in(x::CategoricalValue, y::AbstractRange{T}) where {T<:Integer} = unwrap(x) in y

Base.hash(x::CategoricalValue, h::UInt) = hash(unwrap(x), h)

# Method defined even on unordered values so that sort() works
function Base.isless(x::CategoricalValue, y::CategoricalValue)
    if pool(x) === pool(y) || pool(x) == pool(y)
        return levelcode(x) < levelcode(y)
    else
        throw(ArgumentError("CategoricalValue objects from unequal pools cannot be tested for order"))
    end
end

Base.isless(x::CategoricalValue, y::SupportedTypes) =
    throw(ArgumentError("cannot compare a `CategoricalValue` to value `v` of type " *
                        "`$(typeof(x))`: wrap `v` using `CategoricalValue(v, catvalue)` " *
                        "or `CategoricalValue(v, catarray)` first"))
Base.isless(y::SupportedTypes, x::CategoricalValue) = isless(x, y)

function Base.:<(x::CategoricalValue, y::CategoricalValue)
    poolx = pool(x)
    pooly = pool(y)
    if poolx === pooly || poolx == pooly
        if isordered(poolx) && isordered(pooly)
            return levelcode(x) < levelcode(y)
        else
            throw(ArgumentError("Unordered CategoricalValue objects cannot be tested for order using <. Use isless instead, or call the ordered! function on the parent array to change this"))
        end
    else
        throw(ArgumentError("CategoricalValue objects from unequal pools cannot be tested for order"))
    end
end

Base.:<(x::CategoricalValue, y::SupportedTypes) =
    throw(ArgumentError("cannot compare a `CategoricalValue` to value `v` of type " *
                        "`$(typeof(x))`: wrap `v` using `CategoricalValue(v, catvalue)` " *
                        "or `CategoricalValue(v, catarray)` first"))
Base.:<(y::SupportedTypes, x::CategoricalValue) = x < y

# JSON of CategoricalValue is JSON of the value it refers to
JSON.lower(x::CategoricalValue) = JSON.lower(unwrap(x))
DataAPI.defaultarray(::Type{CategoricalValue{T, R}}, N) where {T, R} =
  CategoricalArray{T, N, R}
DataAPI.defaultarray(::Type{Union{CategoricalValue{T, R}, Missing}}, N) where {T, R} =
  CategoricalArray{Union{T, Missing}, N, R}

# define appropriate handlers for JSON3 interface
StructTypes.StructType(x::CategoricalValue) = StructTypes.StructType(unwrap(x))
StructTypes.StructType(::Type{<:CategoricalValue{T}}) where {T} = StructTypes.StructType(T)
StructTypes.numbertype(::Type{<:CategoricalValue{T}}) where {T <: Number} = T
StructTypes.construct(::Type{T}, x::CategoricalValue{T}) where {T} = T(unwrap(x))

RecipesBase.@recipe function f(::Type{T}, v::T) where T <: CategoricalValue
    level_strings = [map(string, levels(v)); missing]
    ticks --> eachindex(level_strings)
    v -> ismissing(v) ? length(level_strings) : Int(refcode(v)),
    i -> level_strings[Int(i)]
end


text:
function _result_dataframe_vars(variable::JuMP.Containers.DenseAxisArray, get_function)

    if length(axes(variable)) == 1
        result = Vector{Float64}(undef, length(first(variable.axes)))

        for t in variable.axes[1]
            result[t] = get_function(variable[t])
        end

        return DataFrames.DataFrame(var = result)

    elseif length(axes(variable)) == 2

        result = Array{Float64, 2}(undef, length(last(variable.axes)), length(first(variable.axes)))
        names = Array{Symbol, 1}(undef, length(variable.axes[1]))

        for t in variable.axes[2], (ix, name) in enumerate(variable.axes[1])
            result[t, ix] = get_function(variable[name, t])
            names[ix] = Symbol(name)
        end

        return DataFrames.DataFrame(result, names)

    elseif length(axes(variable)) == 3
        extra_dims = sum(length(axes(variable)[2:end-1]))
        extra_vars = [Symbol("S$(s)") for s in 1:extra_dims]
        result_df = DataFrames.DataFrame()
        names = vcat(extra_vars, Symbol.(axes(variable)[1]))

        for i in variable.axes[2]
            third_dim = collect(fill(i,size(variable)[end]))
            result = Array{Float64 ,2}(undef, length(last(variable.axes)),
                                              length(first(variable.axes)))
            for t in last(variable.axes), (ix, name) in enumerate(first(variable.axes))
                result[t, ix] = get_function(variable[name, i, t])
            end
            res = DataFrames.DataFrame(hcat(third_dim, result))
            result_df = vcat(result_df, res)
        end

        return DataFrames.names!(result_df, names)

    else
        error("Dimension Number $(length(axes(variable))) not Supported")

    end

end

function _result_dataframe_vars(variable::JuMP.Containers.SparseAxisArray, get_function)
    return Vector()
end

function get_model_result(m::JuMP.Model)

    results_dict = Dict{Symbol, Any}()

    for (k, v) in m.obj_dict
        results_dict[k] = _result_dataframe_vars(v, JuMP.value)
    end

    results_dict[:pw_ub] = _result_dataframe_vars(m.obj_dict[:pw], JuMP.upper_bound)
    results_dict[:il_ub] = _result_dataframe_vars(m.obj_dict[:pl], JuMP.upper_bound)

    return results_dict

end

function get_optimizer_log(JuMPmodel::JuMP.Model)

    optimizer_log = Dict{Symbol, Any}()

    optimizer_log[:obj_value] = JuMP.objective_value(JuMPmodel)
    optimizer_log[:termination_status] = JuMP.termination_status(JuMPmodel)
    optimizer_log[:primal_status] = JuMP.primal_status(JuMPmodel)
    optimizer_log[:dual_status] = JuMP.dual_status(JuMPmodel)
    optimizer_log[:solver] =  JuMP.solver_name(JuMPmodel)
    try
        optimizer_log[:solve_time] = MOI.get(JuMPmodel, MOI.SolveTime())
    catch
        @warn("SolveTime() property not supported by $(optimizer_log[:solver])")
        optimizer_log[:solve_time] = "Not Supported by $(optimizer_log[:solver])"
    end
    return optimizer_log
end

ext: jl


text:
function Nk(J, Bk)
    return exp(-h*c*Bk*J*(J+1)/(KB*Tk))
end

function Hrot(τ, nJ, gas; pre=-1.0536477760095002e12)
    qt = zeros(Float64, size(τ))
    if gas == :N2
        Bk = 1.99*100 #(m^-1) in QSM paper the constant s in cm^-1
        τrot_k = 61e-12
        for j in 1:nJ
            if j % 2 == true
                qk = 1
            else
                qk = 2
            end
            A_k =(Nk(j+2, Bk) - Nk(j, Bk))*qk*(j + 2)*(j + 1)/(2*j + 3)
            qt += @. (exp(-τ/τrot_k)*A_k*sin(4pi*Bk*c*(2*j+3)*τ))
        end
    end
    if gas == :O2
        Bk = 1.44*100 #(m^-1) in QSM paper the constant s in cm^-1
        τrot_k = 68e-12
        for j in 1:nJ
            if j % 2 == true
                qk = 1
            else
                qk = 0
            end
            A_k =(Nk(j+2, Bk) - Nk(j, Bk))*qk*(j + 2)*(j + 1)/(2*j + 3)
            qt += @. (exp(-τ/τrot_k)*A_k*sin(4pi*Bk*c*(2*j+3)*τ))
        end
    end
     # ret[τ<0] = 0
    @. qt[τ<0] = 0
    return qt .* pre
end

function Hvib(τ, nJ, gas; pre=2.8945247923481406e12)
    qt = zeros(Float64, size(τ))

    if gas == :N2
        Bk = 1.99*100
        τvib_k = 10e-12
        Ω_k = 2330.0*100.0 #(m^-1)
        η_k = 0.0173*100.0
        for j in 1:nJ
            if j % 2 == true
                qk = 1
            else
                qk = 2
            end
            M_k = qk*(2*j+1)*exp(-h*c*Bk*j*(j+1)/(KB*Tk))
            ω_k = 2*pi*c*(Ω_k - η_k*j*(j+1))
            @. qt += (exp(-τ/τvib_k)*(M_k*sin(ω_k*τ)))
        end
    end
    if gas == :O2
        Bk = 1.44*100
        τvib_k = 10e-12
        Ω_k = 1556.0*100.0 #(m^-1)
        η_k = 0.0159*100.0
        for j in 1:nJ
            if j % 2 == true
                qk = 1
            else
                qk = 0
            end
            M_k = qk*(2*j+1)*exp(-h*c*Bk*j*(j+1)/(KB*Tk))
            ω_k = 2*pi*c*(Ω_k - η_k*j*(j+1))
            @. qt += (exp(-τ/τvib_k)*(M_k*sin(ω_k*τ)))
        end
    end
    @. qt[τ<0] = 0
    return qt .* pre
end

function Htot(τ, gas, p, Tk; nJ=30)
    Ht = zeros(Float64, size(τ))
    if gas == :Air
        μrot1 = 0.986
        μrot2 = 1
        H1 = μrot1.*Hrot(τ, nJ, :N2) + (1-μrot1).*Hvib(τ, nJ, :O2)
        H2 = μrot2.*Hrot(τ, nJ, :N2) + (1-μrot2).*Hvib(τ, nJ, :O2)
        nrN2 = n_2(:N2,  p, Tk)
        nrO2 = n_2(:O2,  p, Tk)
        nrAir = n_2(:Air,  p, Tk)
        Ht = 0.8*nrN2/nrAir*(μrot1.*Hrot(τ, nJ, :N2) + (1-μrot1).*Hvib(τ, nJ, :N2)) +
             0.2*nrO2/nrAir*(μrot2.*Hrot(τ, nJ, :N2) + (1-μrot2).*Hvib(τ, nJ, :N2))
        return Ht
    elseif gas == :N2
        μrot = 0.986
        μrot.*Hrot(τ, nJ, gas) + (1-μrot).*Hvib(τ, nJ, gas)
    end
end

function getpre(nJ, gas)
    τ = LinRange(-1.0, 1000.0,3*10^6)*1e-12
    Hr = Hrot(τ, nJ, gas, pre=1.0)
    totr = trapz(τ, Hr)
    Hv = Hvib(τ, nJ, gas, pre=1.0)
    totv = trapz(τ, Hv)
    return 1.0/totr, 1.0/totv
end

function sdo(τ; τ1 = 62.5e-15, τ2=77e-15)
    RT = (τ1^2+τ2^2)/(τ1*τ2^2).*exp.(-τ/τ2).*sin.(τ/τ1)
    RT[τ.<0] .= 0.0
    return RT
end

function raman(fr, γ, convt)
    if fr != 0.0
        γ = γ/(1-fr)
        RW = Htot(τ)
        @. qt[τ<0] = 0
        Ht = f2(Ht, τ)
        @. RW[τ<0] = 0
        RW /= trapz(τ, RW)
        RW = shifted_ifft(RW)*convt
    else
        RW = 0.0
    end
    return RW
end

function covRaman(U, RW, convt, convw, fr)
    # println(fr)
    Uw = shifted_ifft(abs2.(U)).*convt
    Uw = shifted_fft(Uw.*RW) .* U .* convw
    Uw = shifted_ifft(Uw) .* convt
    @. Uw =  (1.0 - fr) * Uw + fr* Uw
    return Uw
end

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
using Test, ScientificTypesBase, Tables
const ST = ScientificTypesBase

include("convention.jl")

include("scitype.jl")


text:
begin
	SIZE=16
	A = reshape(collect(0:SIZE*SIZE-1) % SIZE + 1,SIZE,SIZE)
	function F(x,y)
			A[y,x]
	end
	B = [F(x,y) for x=collect(1:SIZE),y=collect(1:SIZE)]

	C = Array(Array,6)
	C[1] = ((A.*B)&(2^0)).==2^0
	C[2] = ((A.*B)&(2^1)).==2^1
	C[3] = ((A.*B)&(2^2)).==2^2
	C[4] = ((A.*B)&(2^3)).==2^3
	C[5] = ((A.*B)&(2^4)).==2^4
	C[6] = ((A.*B)&(2^5)).==2^5

	K = Array(Array,6)
	K[1] = (A.*B).<=2^1
	K[2] = (A.*B).<=2^2
	K[3] = (A.*B).<=2^3
	K[4] = (A.*B).<=2^4
	K[5] = (A.*B).<=2^5
	K[6] = (A.*B).<=2^6

	function custom_factorize(p)
		L = 1
		W = trues(SIZE,SIZE)
		while 2^L <= SIZE
			W = K[L]
			L += 1
		end

		N = [2^x for x=[0:L]]
		I = 0
		for s in N
			I += 1
			if (p & s) != 0
				W &= C[I]
			else
				W &= !C[I]
			end
		end
		#*1 for smaller output
		return W*1
	end
	custom_factorize(12)
end

max_stars_repo_name: TheCrazyT/factorisation-maps

max_issues_repo_name: TheCrazyT/factorisation-maps

max_forks_repo_name: TheCrazyT/factorization-maps


text:
#=
   Credit-card fraud.

   From David Heckerman
   "A Tutorial on Learning With Bayesian Networks"
   page 12ff
   """
   To illustrate the process of building a Bayesian network, consider the problem of de-
   tecting credit-card fraud. We begin by determining the variables to model. One possible
   choice of variables for our problem is Fraud(F), Gas(G), Jewelry(J) Age(A) and Sex(S)
   representing whether or not the current purchase is fraudulent, whether or not there
   was a gas purchase in the last 24 hours, whether or not there was a jewelry purchase in
   the last 24 hours, and the age and sex of the card holder, respectively.
   """

   Distributions of variable fraud (num:0)
   1.00000 =>   20244  (0.506100)
   0.00000 =>   19756  (0.493900)

   Distributions of variable sex
   male       =>   21678  (0.541950)
   female     =>   18322  (0.458050)

   Distributions of variable age
   age_30_to_50 =>   17315  (0.432875)
   age_above_50 =>   15237  (0.380925)
   age_below_30 =>    7448  (0.186200)

   Distributions of variable gas (num:0)
   0.00000 =>   39821  (0.995525)
   1.00000 =>     179  (0.004475)

   Distributions of variable jewlery (num:0)
   0.00000 =>   39980  (0.999500)
   1.00000 =>      20  (0.000500)

   Cf ~/webppl/credit_card_fraud.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function credit_card_fraud()
    male = 1
    female = 2

    #  type Age;
    age_below_30 = 1
    age_30_to_50 = 2
    age_above_50 = 3

    fraud ~   flip(0.00001)
    age ~ Categorical(simplex([0.25, 0.40, 0.35]))  # [age_below_30, age_30_to_50, age_above_50]
    sex ~ Categorical(simplex([0.5,0.5]))           # [male, female]

    gas ~ fraud ? flip(0.2) : flip(0.01);

    jewlery ~  fraud ? flip(0.05) :
        (age == age_below_30 && sex==male)   ? flip(0.0001) :
        (age == age_30_to_50 && sex==male)   ? flip(0.0004) :
        (age == age_above_50 && sex==male)   ? flip(0.0002) :
        (age == age_below_30 && sex==female) ? flip(0.0005) :
        (age == age_30_to_50 && sex==female) ? flip(0.002)  :
        (age == age_above_50 && sex==female) ? flip(0.001)  : Dirac(0)

    # true ~ Dirac(age == age_below_30)
    # true ~ Dirac(age == age_30_to_50)
    # true ~ Dirac(age == age_above_50)
    # true ~ Dirac(gas == true)
    # true ~ Dirac(sex == female)
    # true ~ Dirac(sex == male)
    # true ~ Dirac(jewlery == true)
    true ~ Dirac(fraud == true)
    # true ~ Dirac(fraud == false)

end

model = credit_card_fraud()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
chns = sample(model, PG(5), 40_000)
# chns = sample(model, SMC(), 40_000)
# chns = sample(model, IS(), 100_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :fraud)
show_var_dist_pct(chns, :sex,["male","female"])
show_var_dist_pct(chns, :age,["age_below_30","age_30_to_50","age_above_50"])
show_var_dist_pct(chns, :gas)
show_var_dist_pct(chns, :jewlery)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
@testset "HeapKNNGraphEdge tests" begin

    @testset "constructor" begin
        @test_throws ErrorException HeapKNNGraphEdge(1, 1, rand())
        @test_throws ErrorException HeapKNNGraphEdge(2, 2, rand(), true)
        @test HeapKNNGraphEdge(1, 2, rand()) isa HeapKNNGraphEdge{Int, Float64}
        @test HeapKNNGraphEdge(2, 1, rand(), false) isa HeapKNNGraphEdge
    end
    @testset "interface" begin
        @testset "comparison" begin
            e1 = HeapKNNGraphEdge(1, 2, .5, true)
            e2 = HeapKNNGraphEdge(1, 2, .5, true)
            e3 = HeapKNNGraphEdge(1, 2, .5, false)
            @test e1 == e2
            @test e1 === e2
            @test e1 == e3 # `flag` does not participate in normal equality
            @test e1 !== e3 # but does for `===`
            # floating point equality has a tolerance
            e4 = HeapKNNGraphEdge(1, 2, .5 + .25*sqrt(eps(Float64)), false)
            @test e3 == e4 # within tol
            @test !(e3 < e4)
            e5 = HeapKNNGraphEdge(1, 2, .5 + sqrt(eps(Float64)), false)
            @test e3 != e5 # outside of tol
            @test e3 < e5
            # rational weights
            e1 = HeapKNNGraphEdge(1, 2, 1//3, true)
            e2 = HeapKNNGraphEdge(1, 2, 1//3, true)
            e3 = HeapKNNGraphEdge(1, 2, 1//2, true)
            @test e1 == e2
            @test e1 === e2
            @test e1 != e3
            @test e1 < e3
        end
        e1 = HeapKNNGraphEdge(1, 2, .5, true)
        @test eltype(e1) == Int

        @testset "property methods" begin
            e1 = HeapKNNGraphEdge(3, 4, 0., true)
            @test src(e1) == 3
            @test dst(e1) == 4
            @test weight(e1) == 0.
            @test flag(e1)
            rev_e1 = reverse(e1)
            e2 = HeapKNNGraphEdge(4, 3, 0., true)
            @test rev_e1 === e2
        end

    end
end


text:
function preprocessSpectral(H ,g, params::Tparams,n1,n2)
    #global V, D, rho, ncomp

    # Δ, V = eig(H)
    X = eigen(H)
    Δ = X.values; V = X.vectors

    l_m, = findmin(Δ)
    g̃ = V'*g
    n_g = norm(g)
    ɛ = 1.0e-10
    λ = max(-l_m,0.0)

    return  PDataSpectral(V,Δ,g̃,λ,true,true)
end

function AInv(X :: PDataSpectral, d̃ ::  Array{Float64,1})
    return X.V * d̃
end

function reconstructH(X :: PDataSpectral)
    return X.V*diagm(X.Δ)*X.V'
end


text:
function build_geometric_brownian_motion_model(μ::Float64 = 0.0,σ::Float64 = 0.0)::PSResult

    # TODO - checks: ate the parameters passed in legit?
    if (σ < 0.0)
        return PSResult{ArgumentError}(ArgumentError("Postive values are expected for σ"))
    end

    # build -
    type_object = PSGeometricBrownianMotionModelParameters(μ,σ)

    # return wrapped -
    return PSResult{PSGeometricBrownianMotionModelParameters}(type_object)
end

function build_geometric_brownian_motion_model(dict::Dict{String,Any})::PSResult

    # Check: does the dictionary have the correct key?
    if (haskey(dict,"underlying_model_parameters") == false)
        return PSResult{ArgumentError}(ArgumentError("Dictionary does not contain the underlying_model_parameters key"))
    end

    # get data -
    underlying_dictionary = dict["underlying_model_parameters"]

    # Check: is the data correct in the underlying_dictionary?
    μ = underlying_dictionary["price_growth_rate"]
    σ = underlying_dictionary["price_volatility"]
    if (μ < 0.0 || σ < 0.0)
        return PSResult{ArgumentError}(ArgumentError("Postive values are expected for both μ and σ"))
    end

    # call -> return PSResult wrapper model -
    return build_geometric_brownian_motion_model(μ,σ)
end


text:
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()
using Revise
using Zygote, Optim
using OrdinaryDiffEq, Flux, DiffEqFlux, LinearAlgebra, Plots
using DiffEqSensitivity
const EIGEN_EST = Ref(0.0f0)
const USE_GPU = Ref(false)

USE_GPU[] = false # the network is small enough such that CPU works just great

USE_GPU[] && (using CUDA)

_gpu(arg) = USE_GPU[] ? gpu(arg) : cpu(arg)
_cu(arg) = USE_GPU[] ? cu(arg) : identity(arg)

function getops(grid, T=Float32)
    N, dz = length(grid), step(grid)
    d = ones(N-2)
    dl = ones(N-3) # super/lower diagonal
    zv = zeros(N-2) # zero diagonal used to extend D* for boundary conditions

    # D1 first order discretization of ∂_z
    D1= diagm(-1 => -dl, 0 => d)
    D1_B = hcat(zv, D1, zv)
    D1_B[1,1] = -1
    D1_B = _cu((1/dz)*D1_B)

    # D2 discretization of ∂_zz
    D2 = diagm(-1 => dl, 0 => -2*d, 1 => dl)
    κ = 0.05
    D2_B = hcat(zv, D2, zv) #add space for the boundary conditions space for "ghost nodes"
    #we only solve for the interior space steps
    D2_B[1,1] = D2_B[end, end] = 1

    D2_B = _cu((κ/(dz^2)).*D2_B) #add the constant κ as the equation requires and finish the discretization

    # Boundary conditions matrix QQ
    Q = Matrix{Int}(I, N-2, N-2)
    QQ = _cu(vcat(zeros(1,N-2), Q, zeros(1,N-2)))

    D1 = D1_B * QQ
    D2 = D2_B * QQ
    EIGEN_EST[] = maximum(abs, eigvals(Matrix(D2)))
    return (D1=T.(D1), D2=T.(D2))
end

function getu0(grid, T=Float32)
    z = grid[2:N-1]
    f0 = z -> T(exp(-200*(z-0.75)^2))
    u0 = f0.(z) |> _gpu
end

function ode_i(u, p, t)
    Φ = u -> cos.(sin.(u.^3) .+ sin.(cos.(u.^2)))
    return p.D1*Φ(u) + p.D2*u
end

function ground_truth(grid, tspan)
    prob = ODEProblem(ode_i, u0, (tspan[1], tspan[2]+0.1), ops)
    sol = solve(prob, ROCK4(eigen_est = (integ)->integ.eigen_est = EIGEN_EST[]), abstol = 1e-9, reltol = 1e-9)
    return sol
end

N = 32
grid = range(0, 1, length = N)
tspan = (0.0f0, 1.5f0)
u0 = getu0(grid)
ops = getops(grid)
soldata = ground_truth(grid, tspan)

ann = Chain(Dense(30,8,tanh), Dense(8,30,tanh)) |> _gpu
pp, re = Flux.destructure(ann)
lyrs = Flux.params(pp)
function dudt_(u,p,t)
    Φ = re(p)
    return ops.D1*Φ(u) + ops.D2*u
end

function predict_adjoint(fullp)
    Array(concrete_solve(prob,
                         ROCK4(eigen_est = (integ)->integ.eigen_est = EIGEN_EST[]),
    u0, fullp, saveat = saveat))
end

function loss_adjoint(fullp)
    pre = predict_adjoint(fullp)
    sum(abs2, training_data - pre)
end

function cb(opt_state:: Optim.OptimizationState)
    cur_pred = collect(predict_adjoint(opt_state.metadata["x"]))
    n = size(training_data, 1)
    pl = scatter(1:n,training_data[:,10],label="data", legend =:bottomright,title="Spatial Plot at t=$(saveat[10])")
    scatter!(pl,1:n,cur_pred[:,10],label="prediction")
    pl2 = scatter(saveat,training_data[N÷2,:],label="data", legend =:bottomright, title="Timeseries Plot at Middle X")
    scatter!(pl2,saveat,cur_pred[N÷2,:],label="prediction")
    display(plot(pl, pl2, size=(600, 300)))
    display(opt_state.value)
    false
end

cb(trace::Optim.OptimizationTrace) = cb(last(trace))

saveat = range(tspan..., length = 30) #time range
prob = ODEProblem{false}(dudt_,u0,tspan,pp)
training_data = _cu(soldata(saveat))
concrete_solve(prob, ROCK4(eigen_est = (integ)->integ.eigen_est = EIGEN_EST[]), u0, pp)
loss_adjoint(pp)

function loss_adjoint_gradient!(G, fullp)
    G .= Zygote.gradient(loss_adjoint, fullp)[1]
end

result =  optimize(loss_adjoint, loss_adjoint_gradient!, pp, BFGS(), Optim.Options(extended_trace=true,callback = cb))

prob2 = ODEProblem{false}(dudt_,u0,(0f0,10f0),pp)
@time full_sol = solve(prob2,
                       ROCK2(eigen_est = (integ)->integ.eigen_est = EIGEN_EST[]),
                       saveat = saveat, abstol=1e-4, reltol=1e-2)

cur_pred = collect((predict_adjoint(result.minimizer)))
n = size(training_data, 1)
pl = scatter(1:n,training_data[:,10],label="data", legend =:bottomright)
scatter!(pl,1:n,cur_pred[:,10],label="prediction",title="Spatial Plot at t=$(saveat[10])")
pl2 = scatter(saveat,training_data[N÷2,:],label="data", legend =:bottomright, title="Time Series Plot: Middle X")
scatter!(pl2,saveat,cur_pred[N÷2,:],label="prediction")
plot(pl, pl2, size=(600, 300))
savefig("npde_fit.pdf")


text:
# TODO: should probably be elsewhere
module BipedControlUtil
    export
        Side,
        left,
        right,
        flipsign_if_right

    # Side
    # From https://bitbucket.org/ihmcrobotics/ihmc-open-robotics-software/src/575d33e1ab064f4e5957096414376fc011a98bce/IHMCRoboticsToolkit/src/us/ihmc/robotics/robotSide/RobotSide.java?at=develop&fileviewer=file-view-default
    @enum Side left right
    flipsign_if_right(x::Number, side::Side) = ifelse(side == right, -x, x)
    Base.rand(rng::AbstractRNG, ::Type{Side}) = ifelse(rand(rng, Bool), left, right)
    Base.:(-)(side::Side) = ifelse(side == right, left, right)
end


text:
using Plots, StatsBase, Statistics, LaTeXStrings, JLD, StatsPlots

function FindBorder(dim, Network)
    border0 = []
    for i in 1:dim
        for j in 1:dim
            check = false
            if Network[i,j] == 0
                Neighbors = ((i-1,j),(i,j-1),(i+1,j),(i,j+1))
                for nei in Neighbors
                    if dim+1 ∉ nei && 0 ∉ nei && Network[nei...] == 1
                        check = true
                        break
                    end
                end
                if check == true
                    push!(border0,(i,j))
                end
            end
        end
    end
    return border0
end

function FindCorrelationLength(Network,dim)
    FindOnSpots = findall(x->x==1,Network)
    SpotNums = length(FindOnSpots)
    ΣR² = 0.0; iMC = 0.0; jMC = 0.0
    for indx in FindOnSpots
        iMC += indx[1]/SpotNums
        jMC += indx[2]/SpotNums
    end
    for indx in FindOnSpots
        ΣR² += ((indx[1]-iMC)^2 + (indx[2]-jMC)^2)/SpotNums
    end
    return sqrt(ΣR²)
end

function ClusterGrowth(dim,P)
    TotalNetworks = []
    Network = zeros(Int,dim,dim)
    # Network[rand(1:dim^2)] = 1
    Network[round(Int, dim/2), round(Int, dim/2)] = 1
    Borders = FindBorder(dim,Network)
    while length(Borders) != 0
        push!(TotalNetworks, copy(Network))
        Borders = FindBorder(dim,Network)
        for spot in Borders
            if P > rand()
                Network[spot...] = 1
            else
                Network[spot...] = -1
            end
        end
    end
    return TotalNetworks
end

PList = [0.1,0.3,0.5,0.7,0.9]
for P in PList
    save("../../Data/Q3/Q3-anim-P=$P.jld","Networks",Networks)
    Networks = ClusterGrowth(150,0.6)
    anim = Animation()
    for i in 1:length(Networks)
        plt = heatmap(Networks[i], legend = false, border = :none, title = L"Growing\ Cluster\ for\ L= 150\ , P\ =\ %$P")
        frame(anim, plt)
    end
    gif(anim, "../../Figs/Q3/Q3-anim-P=$P.gif")
end


text:
using PyPlot
using JLD

function main(code::AbstractString,n::Int64)
    log_path_LMPC   = "$(homedir())/simulations/output-LMPC-$(code).jld"
    log_path_record = "$(homedir())/simulations/output-record-$(code).jld"
    #log_path_record = "$(homedir())/simulations/output-record-$(code).jld"
    #d_rec       = load(log_path_record)
    d_lmpc      = load(log_path_LMPC)
    d_rec       = load(log_path_record)

    delay_df    = 5

    #n = 2

    #idx0 = (oldTraj.oldTraj[:,6,n].> 0) & (oldTraj.oldTraj[:,6,n].< 17.76)
    #oldTraj.oldTraj  = oldTraj.oldTraj[idx0,:,:]
    #oldTraj.oldInput = oldTraj.oldInput[idx0,:,:]
    oldTraj     = d_lmpc["oldTraj"]
    pos_info    = d_rec["pos_info"]
    imu_meas    = d_rec["imu_meas"]
    sz = size(oldTraj.oldTraj[oldTraj.oldTraj[:,1,n].<1000,:,n],1)
    oldTraj.oldTraj[1:sz,:,n] = smooth(oldTraj.oldTraj[1:sz,:,n],20)
    #oldTraj.oldInput[1:sz,:,n] = smooth(oldTraj.oldInput[1:sz,:,n],2)
    # plot(oldTraj.oldInput[1:sz,:,n])
    # for i=2:sz
    #     oldTraj.oldInput[i,:,n] = oldTraj.oldInput[i-1,:,n] + (oldTraj.oldInput[i,:,n]-oldTraj.oldInput[i-1,:,n])*0.5
    # end
    # plot(oldTraj.oldInput[1:sz,:,n],"--")
    #sz = 200
    idx_d = 20:5:sz-5
    idx   = 20:5:sz-10
    sz1 = size(idx,1)

    y_psi = zeros(5*sz1)
    A_psi = zeros(5*sz1,3)
    y_vy = zeros(5*sz1)
    A_vy = zeros(5*sz1,4)
    y_vx = zeros(5*sz1)
    A_vx = zeros(5*sz1,4)

    for i=0:4
        y_psi[(1:sz1)+i*sz1] = diff(oldTraj.oldTraj[idx_d+i,3,n])
        A_psi[(1:sz1)+i*sz1,:] = [oldTraj.oldTraj[idx+i,3,n]./oldTraj.oldTraj[idx+i,1,n] oldTraj.oldTraj[idx+i,2,n]./oldTraj.oldTraj[idx+i,1,n] oldTraj.oldInput[idx+i-delay_df,2,n]]
    end
    c_psi = (A_psi'*A_psi)\A_psi'*y_psi

    for i=0:4
        y_vy[(1:sz1)+i*sz1] = diff(oldTraj.oldTraj[idx_d+i,2,n])
        A_vy[(1:sz1)+i*sz1,:] = [oldTraj.oldTraj[idx+i,2,n]./oldTraj.oldTraj[idx+i,1,n] oldTraj.oldTraj[idx+i,1,n].*oldTraj.oldTraj[idx+i,3,n] oldTraj.oldTraj[idx+i,3,n]./oldTraj.oldTraj[idx+i,1,n] oldTraj.oldInput[idx+i-delay_df,2,n]]
    end
    c_vy = (A_vy'*A_vy)\A_vy'*y_vy

    for i=0:4
        y_vx[(1:sz1)+i*sz1] = diff(oldTraj.oldTraj[idx_d+i,1,n])
        A_vx[(1:sz1)+i*sz1,:] = [oldTraj.oldTraj[idx+i,2,n] oldTraj.oldTraj[idx+i,3,n] oldTraj.oldTraj[idx+i,1,n] oldTraj.oldInput[idx+i,1,n]]
    end
    c_vx = (A_vx'*A_vx)\A_vx'*y_vx

    println("psi:")
    println(c_psi)
    println("vy:")
    println(c_vy)
    println("vx:")
    println(c_vx)

    figure(1)
    plot(A_psi*c_psi)
    plot(y_psi,"--")
    #plot(A_psi,"--")
    title("Psi")

    figure(2)
    plot(A_vy*c_vy)
    plot(y_vy,"--")
    #plot(A_vy,"--")
    title("Vy")

    figure(3)
    plot(A_vx*c_vx)
    plot(y_vx,"--")
    #plot(A_vx,"--")
    title("Vx")

    # figure(4)
    # plot(oldTraj.oldTimes[:,1],oldTraj.oldTraj[:,3,1])
    # plot(pos_info.t,pos_info.z[:,11])
    # plot(imu_meas.t,imu_meas.z[:,3])
    nothing
end

function smooth(x,n)
    y = zeros(size(x))
    for i=1:size(x,1)
        start = max(1,i-n)
        fin = min(size(x,1),start + 2*n)
        y[i,:] = mean(x[start:fin,:],1)
    end
    return y
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
module CUDA

using ..CuArrays

if !applicable(CuArray{UInt8}, undef, 1)
  (T::Type{<:CuArray})(::UndefInitializer, sz...) = T(sz...)
end

if CuArrays.libcudnn != nothing
  include("cudnn.jl")
else
  @warn("CUDNN is not installed, some functionality will not be available.")
end

end


text:
## Model
## Two equations NK model on the form (1)-(2):
# (1) A*Y(+1) = B*Y + C*X
# (2) X(+1) = D*X + E*ϵ, ϵ ~ N(0,Σ)
# Solution has the form of (3):
# (3) Y(+1) = inv(A)*B*Y + inv(A)*C*X
# provided A full rank and eigenvalues in inv(A)*B respecting Blanchard and Kahn (1980) conditions

## Example
# NKIS: x(t+1) + σ^(-1)π(t+1) = x(t) + σ^(-1)δπ(t) + σ^(-1)v(t)
# NKPC: βπ(t+1) = -kx(t) +π(t)
using DsgeSolves, Plots
plotly()

# Declare number of variables
nb = 2 # #backwardlooking variables
nf = 2 # #forward looking variables
ns = 2 # #forcing variables
m = DSGE(nb,nf,ns)

# Preallocations and hyperparameters
rep = 1000
burn_in = 100

## Indexing variables for filling A
iv = 1
iu = 2
ix = 3
iπ = 4

# Indexing shock
iϵ_v = 1
iϵ_u = 2

## Calibration
# Structural parameters
β = 0.99
σ = 1
χ = 1.55
η = 0
θ = 2.064
ω = 0.5
α = 3
δ = 1.5
κ = (1-ω)*(1-β*ω)/(α*ω)

# Non-structural parameters
ρ_u = 0.5
ρ_v = 0.8
σ_v = 1
σ_u = 0.5
Σ = diagm([σ_v, σ_u])

## Filling matrices A, B, C
# four variables:
#                  2 forward-looking -- output gap and inflation
#                  2 predetermined -- interest rate shock

# Equation #1 monetary policy shock
# v(t+1) = ρ_vv(t) + ϵ_v(t+1)
nequ = 1
m.A[nequ, iv] = 1
m.B[nequ, iv] = ρ_v
m.C[nequ, iϵ_v] = 1

# Equation #2 cost-push shock
# u(t+1) = ρ_uu(t) + ϵ_u(t+1)
nequ = 2
m.A[nequ, iu] = 1
m.B[nequ, iu] = ρ_u
m.C[nequ, iϵ_u] = 1

# Equation #3 IS curve
# x(t+1) + σ^(-1)π(t+1) = x(t) + σ^(-1)δπ(t) + σ^(-1)v(t)
nequ = 3
m.A[nequ, ix] = 1
m.A[nequ, iπ] = 1/σ
m.B[nequ, iv] = 1/σ
m.B[nequ, ix] = 1
m.B[nequ, iπ] = δ/σ

# Equation #4 NKPC
# βπ(t+1) = -kx(t) +π(t)
nequ = 4
m.A[nequ, iπ] = β
m.B[nequ, iu] = -1
m.B[nequ, ix] = -κ
m.B[nequ, iπ] = 1

## Solve and simulate the model
(A, B, Λ, P, pstar, R) = create_matrices(m)
(w, y) = simulating_dsge(m,pstar,Λ,R,Σ,rep)
# Plot model simulation
pSimulation = plot(layout = grid(2,2));
plot!(pSimulation, w[1,burn_in+1:end], subplot = 1, color = "blue", title = ":ϵ_u");
plot!(pSimulation, w[2,burn_in+1:end], subplot = 2, color = "blue", title = ":ϵ_v");
plot!(pSimulation, y[1,burn_in+1:end], subplot = 3, color = "blue", title = ":x");
plot!(pSimulation, y[2,burn_in+1:end], subplot = 4, color = "blue", title = ":π");
gui(pSimulation)

# Impulse response functions
shk = 1 # choose the shock number
h = 20
(irf_w,irf_y) = irf_dsge(m,pstar,Λ,R,shk,h)
# Plot model IRFs
pIRF = plot(layout = grid(1,2));
plot!(pIRF, irf_y[1,:], subplot = 1, color = "blue", title = ":x");
plot!(pIRF, irf_y[2,:], subplot = 2, color = "blue", title = ":π");
gui(pIRF)

# Computing stylized facts as standard deviation and correlation of series
mean(w,2), mean(y,2)
std(w,2), std(y,2)


text:
local a::Int
local b::Bool

a = 1
b = a


text:
const CMat = Ptr{Cvoid}

abstract type AbstractMat{T} <: AbstractMatrix{T} end
scalartype(::AbstractMat{T}) where {T} = T

# allows us to pass XXMat objects directly into CMat ccall signatures
Base.cconvert(::Type{CMat}, obj::AbstractMat) = obj.ptr
# allows us to pass XXMat objects directly into Ptr{CMat} ccall signatures
Base.unsafe_convert(::Type{Ptr{CMat}}, obj::AbstractMat) =
    convert(Ptr{CMat}, pointer_from_objref(obj))

Base.eltype(::Type{A}) where {A<:AbstractMat{T}} where {T} = T
Base.eltype(A::AbstractMat{T}) where {T} = T

"""
    MatSeqAIJ{T}

PETSc sparse array using AIJ format (also known as a compressed sparse row or CSR format).

Memory allocation is handled by PETSc.
"""
mutable struct MatSeqAIJ{T} <: AbstractMat{T}
    ptr::CMat
    comm::MPI.Comm
end

"""
    MatSeqDense{T}

PETSc dense array. This wraps a Julia `Matrix{T}` object.
"""
mutable struct MatSeqDense{T} <: AbstractMat{T}
    ptr::CMat
    comm::MPI.Comm
    array::Matrix{T}
end

@for_libpetsc begin
    function MatSeqAIJ{$PetscScalar}(m::Integer, n::Integer, nnz::Vector{$PetscInt})
        initialize($PetscScalar)
        comm = MPI.COMM_SELF
        mat = MatSeqAIJ{$PetscScalar}(C_NULL, comm)
        @chk ccall((:MatCreateSeqAIJ, $libpetsc), PetscErrorCode,
                (MPI.MPI_Comm, $PetscInt, $PetscInt, $PetscInt, Ptr{$PetscInt}, Ptr{CMat}),
                comm, m, n, 0, nnz, mat)
        finalizer(destroy, mat)
        return mat
    end
    function MatSeqDense(A::Matrix{$PetscScalar})
        initialize($PetscScalar)
        comm = MPI.COMM_SELF
        mat = MatSeqDense(C_NULL, comm, A)
        @chk ccall((:MatCreateSeqDense, $libpetsc), PetscErrorCode,
            (MPI.MPI_Comm, $PetscInt, $PetscInt, Ptr{$PetscScalar}, Ptr{CMat}),
            comm, size(A,1), size(A,2), A, mat)
        finalizer(destroy, mat)
        return mat
    end

    function destroy(M::AbstractMat{$PetscScalar})
        finalized($PetscScalar) ||
        @chk ccall((:MatDestroy, $libpetsc), PetscErrorCode, (Ptr{CMat},), M)
        return nothing
    end

    function setvalues!(M::AbstractMat{$PetscScalar}, row0idxs::Vector{$PetscInt}, col0idxs::Vector{$PetscInt}, rowvals::Array{$PetscScalar}, insertmode::InsertMode)
        @chk ccall((:MatSetValues, $libpetsc), PetscErrorCode,
            (CMat, $PetscInt, Ptr{$PetscInt}, $PetscInt, Ptr{$PetscInt}, Ptr{$PetscScalar},InsertMode),
            M, length(row0idxs), row0idxs, length(col0idxs), col0idxs, rowvals, insertmode)
        return nothing
    end

    function Base.setindex!(M::AbstractMat{$PetscScalar}, val, i::Integer, j::Integer)
        @chk ccall((:MatSetValues, $libpetsc), PetscErrorCode,
            (CMat, $PetscInt, Ptr{$PetscInt}, $PetscInt, Ptr{$PetscInt}, Ptr{$PetscScalar}, InsertMode),
            M, 1, Ref{$PetscInt}(i-1), 1, Ref{$PetscInt}(j-1), Ref{$PetscScalar}(val), INSERT_VALUES)
        return val
    end

    function assemblybegin(M::AbstractMat{$PetscScalar}, t::MatAssemblyType=MAT_FINAL_ASSEMBLY)
        @chk ccall((:MatAssemblyBegin, $libpetsc), PetscErrorCode, (CMat, MatAssemblyType), M, t)
        return nothing
    end
    function assemblyend(M::AbstractMat{$PetscScalar}, t::MatAssemblyType=MAT_FINAL_ASSEMBLY)
        @chk ccall((:MatAssemblyEnd, $libpetsc), PetscErrorCode, (CMat, MatAssemblyType), M, t)
        return nothing
    end
    function view(mat::AbstractMat{$PetscScalar}, viewer::Viewer{$PetscScalar}=ViewerStdout{$PetscScalar}(mat.comm))
        @chk ccall((:MatView, $libpetsc), PetscErrorCode,
                    (CMat, CPetscViewer),
                mat, viewer);
        return nothing
    end

    function Base.size(A::AbstractMat{$PetscScalar})
        m = Ref{$PetscInt}()
        n = Ref{$PetscInt}()
        @chk ccall((:MatGetSize, libpetsc), PetscErrorCode,
            (CMat, Ptr{$PetscInt}, Ptr{$PetscInt}),
            A, m, n)
        return (m[], n[])
    end
    function Base.:(==)(A::AbstractMat{$PetscScalar}, B::AbstractMat{$PetscScalar})
        fr = Ref{PetscBool}()
        @chk ccall((:MatEqual, $libpetsc), PetscErrorCode,
             (CMat, CMat, Ptr{PetscBool}),
             A, B, fr)
        return fr[] == PETSC_TRUE
    end

    function LinearAlgebra.issymmetric(A::AbstractMat{$PetscScalar}; tol=zero($PetscReal))
        fr = Ref{PetscBool}()
        @chk ccall((:MatIsSymmetric, $libpetsc), PetscErrorCode,
            (CMat, $PetscReal, Ptr{PetscBool}),
            A, tol, fr)
        return fr[] == PETSC_TRUE
    end
    function LinearAlgebra.ishermitian(A::AbstractMat{$PetscScalar}; tol=zero($PetscReal))
        fr = Ref{PetscBool}()
        @chk ccall((:MatIsHermitian, $libpetsc), PetscErrorCode,
            (CMat, $PetscReal, Ptr{PetscBool}),
            A, tol, fr)
        return fr[] == PETSC_TRUE
    end
    function LinearAlgebra.norm(M::AbstractMat{$PetscScalar}, normtype::NormType=NORM_FROBENIUS)
        r_val = Ref{$PetscReal}()
        @chk ccall((:MatNorm, $libpetsc), PetscErrorCode,
            (CMat, NormType, Ptr{$PetscReal}),
            M, normtype, r_val)
        return r_val[]
    end

    function LinearAlgebra.mul!(y::AbstractVec{$PetscScalar}, M::AbstractMat{$PetscScalar}, x::AbstractVec{$PetscScalar})
        @chk ccall((:MatMult, $libpetsc), PetscErrorCode, (CMat, CVec, CVec), M, x, y)
        return y
    end
    function LinearAlgebra.mul!(y::AbstractVec{$PetscScalar}, M::Adjoint{T,A}, x::AbstractVec{$PetscScalar}) where {T,A<:AbstractMat{$PetscScalar}}
        @chk ccall((:MatMultHermitianTranspose, $libpetsc), PetscErrorCode, (CMat, CVec, CVec), parent(M), x, y)
        return y
    end
    function LinearAlgebra.mul!(y::AbstractVec{$PetscScalar}, M::Transpose{T,A}, x::AbstractVec{$PetscScalar}) where {T,A<:AbstractMat{$PetscScalar}}
        @chk ccall((:MatMultTranspose, $libpetsc), PetscErrorCode, (CMat, CVec, CVec), parent(M), x, y)
        return y
    end

end

function assemble(M::AbstractMat, t::MatAssemblyType=MAT_FINAL_ASSEMBLY)
    assemblybegin(M, t)
    assemblyend(M, t)
end

function MatSeqAIJ(S::SparseMatrixCSC{T}) where {T}
    PetscInt = inttype(T)
    m,n = size(S)
    nnz = zeros(PetscInt,m)
    for r in S.rowval
        nnz[r] += 1
    end
    M = MatSeqAIJ{T}(m, n, nnz)
    for j = 1:n
        for ii = S.colptr[j]:S.colptr[j+1]-1
            i = S.rowval[ii]
            M[i,j] = S.nzval[ii]
        end
    end
    assemble(M)
    return M
end

function Base.copyto!(M::PETSc.MatSeqAIJ{T}, S::SparseMatrixCSC{T}) where {T}
    for j = 1:size(S,2)
        for ii = S.colptr[j]:S.colptr[j+1]-1
            i = S.rowval[ii]
            M[i,j] = S.nzval[ii]
        end
    end
    assemble(M);
end

function Base.show(io::IO, ::MIME"text/plain", mat::AbstractMat)
    _show(io, mat)
end
AbstractMat(A::Matrix) = MatSeqDense(A)
AbstractMat(A::SparseMatrixCSC) = MatSeqAIJ(A)

const MatAT{T} = Union{AbstractMat{T}, Transpose{T, <:AbstractMat{T}}, Adjoint{T, <:AbstractMat{T}}}

LinearAlgebra.mul!(y::AbstractVector{T}, M::MatAT{T}, x::AbstractVector{T}) where {T} =
    parent(LinearAlgebra.mul!(AbstractVec(y), M, AbstractVec(x)))


text:
export Diff

"""
    Diff(f)

Create the derivative of the scalar function `f`, based on `ForwardDiff`.
"""
Diff(f) = x -> ForwardDiff.derivative(f, float(x))

"""
    Diff(f,k)

By recurrence, create the k-th derivative of the scalar function `f`, based on `ForwardDiff`.
"""
Diff(f, k) = k == 0 ? f : k == 1 ? Diff(f) : k>1 ? Diff(Diff(f),k-1) : print("Differential operator not defined for k<0")


text:
import PredictMD
import Test

a = PredictMD.predictionsassoctodataframe(Dict(:x => [0.3], :y => [0.7]), Symbol[])
Test.@test a[1, :x] == 0.3
Test.@test a[1, :y] == 0.7


text:
############################################################################################
# Include files
include("flatten/flatten.jl")
include("utility.jl")
include("checks.jl")
include("random.jl")

############################################################################################
# Export


text:
using POMDPs
using LaserTag
using POMDPPolicies
using POMDPSimulators
using ParticleFilters
using ProfileView
using Random
using Profile

p = gen_lasertag(3,4,2)

pol = RandomPolicy(p, rng=MersenneTwister(1))

sim = RolloutSimulator(max_steps=100, rng=MersenneTwister(2))

fltr = SIRParticleFilter(p, 10000)

@time simulate(sim, p, pol, fltr)
@time simulate(sim, p, pol, fltr)
simulate(sim, p, pol, fltr)

Profile.clear()
@profile simulate(sim, p, pol, fltr)

ProfileView.view()


text:
function grow_region_unwrap!(
    wrapped, weights, visited=zeros(UInt8, size(wrapped)), pqueue=PQueue{Int}(NBINS);
    maxseeds=1, merge_regions=false, correct_regions=false, wrap_addition=0, keyargs...
    )
    ## Init
    maxseeds = min(255, maxseeds) # Hard Limit, Stored in UInt8
    dimoffsets = getdimoffsets(wrapped)
    notvisited(i) = checkbounds(Bool, visited, i) && (visited[i] == 0)
    seeds = Int[]
    new_seed_thresh = 256
    if isempty(pqueue) # no seed added yet
        addseed! = getseedfunction(seeds, pqueue, visited, weights, wrapped, keyargs)
        new_seed_thresh = addseed!()
    end
    ## MST loop
    while !isempty(pqueue)
        if length(seeds) < maxseeds && pqueue.min > new_seed_thresh
            new_seed_thresh = addseed!()
        end
        edge = dequeue!(pqueue)
        oldvox, newvox = getvoxelsfromedge(edge, visited, dimoffsets)
        if visited[newvox] == 0
            unwrapedge!(wrapped, oldvox, newvox, visited, wrap_addition)
            visited[newvox] = visited[oldvox]
            for i in 1:6 # 6 directions
                e = getnewedge(newvox, notvisited, dimoffsets, i)
                if e != 0 && weights[e] > 0
                    enqueue!(pqueue, e, weights[e])
                end
            end
        end
    end
    ## Region merging
    if merge_regions regions = merge_regions!(wrapped, visited, length(seeds), weights) end
    if merge_regions && correct_regions correct_regions!(wrapped, visited, regions) end
    return visited
end

function getvoxelsfromedge(edge, visited, stridelist)
    dim = getdimfromedge(edge)
    vox = getfirstvoxfromedge(edge)
    neighbor = vox + stridelist[dim] # direct neigbor in dim
    if visited[neighbor] == 0
        return vox, neighbor
    else
        return neighbor, vox
    end
end

# edge calculations
getdimfromedge(edge) = (edge - 1) % 3 + 1
getfirstvoxfromedge(edge) = div(edge - 1, 3) + 1
getedgeindex(leftvoxel, dim) = dim + 3(leftvoxel-1)

function unwrapedge!(wrapped, oldvox, newvox, visited, x)
    oo = 2oldvox - newvox
    d = 0
    if checkbounds(Bool, wrapped, oo) && visited[oo] != 0 # neighbor behind is visited
        v = wrapped[oldvox] - wrapped[oo]
        d = if v < -x # threshold
            -x
        elseif v > x
            x
        else
            v
        end
    end
    wrapped[newvox] = unwrapvoxel(wrapped[newvox], wrapped[oldvox] + d)
end
unwrapvoxel(new, old) = new - 2pi * round((new - old) / 2pi)

function getnewedge(v, notvisited, stridelist, i)
    iDim = div(i+1,2)
    n = stridelist[iDim] # neigbor-offset in dimension iDim
    if iseven(i)
        if notvisited(v+n) getedgeindex(v, iDim) else 0 end
    else
        if notvisited(v-n) getedgeindex(v-n, iDim) else 0 end
    end
end


text:
using Test
using DataFrames, DBInterface, CSV, Dates
using DuckDB
#include("../src/DuckDB.jl")

@testset "DB Connection" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    @test isa(con, Base.RefValue{Ptr{Nothing}})
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Appender API" begin
    # Open the database
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)

    # Create the table the data is appended to
    DuckDB.execute(con, "CREATE TABLE dtypes(bol BOOLEAN, tint TINYINT, sint SMALLINT, int INTEGER, bint BIGINT, utint UTINYINT, usint USMALLINT, uint UINTEGER, ubint UBIGINT, float FLOAT, double DOUBLE, date DATE, time TIME, vchar VARCHAR, nullval INTEGER)")

    # Create the appender
    appender = DuckDB.appender_create(con, "dtypes")

    # Append the different data types
    DuckDB.duckdb_append_bool(appender, true)
    DuckDB.duckdb_append_int8(appender, 1)
    DuckDB.duckdb_append_int16(appender, 2)
    DuckDB.duckdb_append_int32(appender, 3)
    DuckDB.duckdb_append_int64(appender, 4)
    DuckDB.duckdb_append_uint8(appender, 1)
    DuckDB.duckdb_append_uint16(appender, 2)
    DuckDB.duckdb_append_uint32(appender, 3)
    DuckDB.duckdb_append_uint64(appender, 4)
    DuckDB.duckdb_append_float(appender, 1.0)
    DuckDB.duckdb_append_double(appender, 2.0)
    DuckDB.duckdb_append_date(appender, 100)
    DuckDB.duckdb_append_time(appender, 200)
    DuckDB.duckdb_append_varchar(appender, "Foo")
    DuckDB.duckdb_append_null(appender)
    # End the row of the appender
    DuckDB.duckdb_appender_end_row(appender)
    # Destroy the appender and flush the data
    DuckDB.duckdb_appender_destroy(appender)

    # Retrive the data from the table and store it in  a vector
    df = DuckDB.toDataFrame(con, "select * from dtypes;")
    data = Matrix(df)

    # Test if the correct types have been appended to the table
    @test data[1] === true
    @test data[2] === Int8(1)
    @test data[3] === Int16(2)
    @test data[4] === Int32(3)
    @test data[5] === Int64(4)
    @test data[6] === UInt8(1)
    @test data[7] === UInt16(2)
    @test data[8] === UInt32(3)
    @test data[9] === UInt64(4)
    @test data[10] === Float32(1.0)
    @test data[11] === Float64(2.0)
    @test data[12] === Dates.Date("1970-04-11")
    @test data[13] === Dates.Time(0, 0, 0, 0, 200)
    @test data[14] === "Foo"
    @test data[15] === missing

    # Disconnect and close the database
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Test append DataFrame" begin
    # Open the database
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)

    # Create the table the data is appended to
    DuckDB.execute(con, "CREATE TABLE dtypes(bool BOOLEAN, tint TINYINT, sint SMALLINT, int INTEGER, bint BIGINT, utint UTINYINT, usint USMALLINT, uint UINTEGER, ubint UBIGINT, float FLOAT, double DOUBLE, date DATE, time TIME, vchar VARCHAR, nullval INTEGER)")

    # Create test DataFrame
    input_df = DataFrame(bool=[true, false], tint=Int8.(1:2), sint=Int16.(1:2), int=Int32.(1:2), bint=Int64.(1:2), utint=UInt8.(1:2), usint=UInt16.(1:2), uint=UInt32.(1:2), ubint=UInt64.(1:2), float=Float32.(1:2), double=Float64.(1:2), date=[Dates.Date("1970-04-11"), Dates.Date("1970-04-12")], time=[Dates.Time(0, 0, 0, 100, 0), Dates.Time(0, 0, 0, 200, 0)], vchar=["Foo", "Bar"], nullval=[missing, Int32(2)])

    # append the DataFrame to the table
    DuckDB.appendDataFrame(input_df, con, "dtypes")

    # Output the data from the table
    output_df = DuckDB.toDataFrame(con, "select * from dtypes;")

    # Compare each column of the input and output dataframe with each other
    for (col_pos, input_col) in enumerate(eachcol(input_df))
        @test isequal(input_col, output_df[:, col_pos])
    end

    # Disconnect and close the database
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Test README" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    res = DuckDB.execute(con, "CREATE TABLE integers(date DATE, jcol INTEGER)")
    res = DuckDB.execute(
        con,
        "INSERT INTO integers VALUES ('2021-09-27', 4), ('2021-09-28', 6), ('2021-09-29', 8)",
    )
    res = DuckDB.execute(con, "SELECT * FROM integers")
    @test isa(DuckDB.toDataFrame(res), DataFrame)
    df = DuckDB.toDataFrame(con, "SELECT * FROM integers")
    @test isa(df, DataFrame)
    DuckDB.appendDataFrame(df, con, "integers")
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "HUGE Int test" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    res = DuckDB.execute(con,"CREATE TABLE huge(id INTEGER,data HUGEINT);")
    res = DuckDB.execute(con,"INSERT INTO huge VALUES (1,NULL), (2, 1761718171), (3, 171661889178);")
    res = DuckDB.toDataFrame(con,"SELECT * FROM huge")
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Interval type" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    res = DuckDB.execute(con, "CREATE TABLE interval(interval INTERVAL);")
    res = DuckDB.execute(
        con,
        """
INSERT INTO interval VALUES
(INTERVAL 5 HOUR),
(INTERVAL 12 MONTH),
(INTERVAL 12 MICROSECOND),
(INTERVAL 1 YEAR);
""",
    )
    res = DuckDB.toDataFrame(con, "SELECT * FROM interval;")
    @test isa(res, DataFrame)
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Timestamp" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    res = DuckDB.execute(con, "CREATE TABLE timestamp(timestamp TIMESTAMP , data INTEGER);")
    res = DuckDB.execute(
        con,
        "INSERT INTO timestamp VALUES ('2021-09-27 11:30:00.000', 4), ('2021-09-28 12:30:00.000', 6), ('2021-09-29 13:30:00.000', 8);",
    )
    res = DuckDB.execute(con, "SELECT * FROM timestamp;")
    res = DuckDB.toDataFrame(res)
    @test isa(res, DataFrame)
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Items table" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    res = DuckDB.execute(
        con,
        "CREATE TABLE items(item VARCHAR, value DECIMAL(10,2), count INTEGER);",
    )
    res = DuckDB.execute(
        con,
        "INSERT INTO items VALUES ('jeans', 20.0, 1), ('hammer', 42.2, 2);",
    )
    res = DuckDB.toDataFrame(con, "SELECT * FROM items;")
    @test isa(res, DataFrame)
    DuckDB.disconnect(con)
end

@testset "Integers and dates table" begin
    db = DuckDB.DB()
    res = DBInterface.execute(db, "CREATE TABLE integers(date DATE, data INTEGER);")
    res = DBInterface.execute(
        db,
        "INSERT INTO integers VALUES ('2021-09-27', 4), ('2021-09-28', 6), ('2021-09-29', 8);",
    )
    res = DBInterface.execute(db, "SELECT * FROM integers;")
    res = DuckDB.toDataFrame(res)
    @test isa(res, DataFrame)
    DBInterface.close!(db)
end

@testset "Query CSV and output DataFrame" begin
    df = DataFrame(a=1:100, b=1:100)
    CSV.write("test_dataframe.csv", df)
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    df1 = DuckDB.toDataFrame(con, "SELECT * FROM 'test_dataframe.csv';")
    @test df == df1
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Export and Query Parquet" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    res = DuckDB.execute(con, "CREATE TABLE integers(date DATE, jcol INTEGER)")
    res = DuckDB.execute(con,
        "INSERT INTO integers VALUES ('2021-09-27', 4), ('2021-09-28', 6), ('2021-09-29', 8)")
    res = DuckDB.execute(con, "COPY (SELECT * FROM integers) TO 'test.parquet' (FORMAT 'parquet');")
    df1 = DuckDB.toDataFrame(con, "SELECT * FROM integers;")
    df2 = DuckDB.toDataFrame(con, "SELECT * FROM 'test.parquet';")
    @test df1 == df2
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

max_issues_repo_name: kimmolinna/julia-duckdb


text:
module SEMFitting
using Images
import DataFrames
using DataFrames: DataFrame, sort!
using Clustering

export approxscale!
export sizerange
export colrow
export groundplane

"""
```
approxscale!(img::AbstractImage, scale, scalerow=1758)
```

Given an SEM image `img` with scale bar of physical length `scale`,
returns physical length per pixel and sets the `spacedirections` property
of the image. Assumes that the scale bar is alone in pixel row `scalerow` of the
image, which defaults to 1758 (good for Painter group SEM).

Currently uses the far extrema of the scale bar to determine length.
"""
function approxscale!(img::AbstractImage, scale, scalerow=1758)
    gray = convert(Image{Images.Gray}, img)
    tf = convert(Array{Bool}, data(gray[:,scalerow]))
    last = tf[1]
    flag1 = true
    local firstedge
    local lastedge
    for (i,x) in enumerate(tf[2:end])
        if x == true && last == false && flag1
            firstedge = i
            flag1 = false
        elseif x == false && last == true
            lastedge = i-1
        end
        last = x
    end
    distpp = scale/(lastedge-firstedge)
    img["spacedirections"] = [distpp,0,0,distpp]
    distpp
end

"""
```
sizerange(img::AbstractImage, scale, min, max; thresh=0.5, tone=true)
```

Does grayscale conversion and thresholding of image `img`, with `thresh` for
thresholding and `tone` for optional inversion (default to true).

Next, extract components of the binary image that are within a specified size
range `min` to `max` where these are units of area in terms of the approximate
physical scale of the image. The `scale` parameter is the length of the scale bar.

Returns `im, df` where `im` is a black and white image only showing the
objects in the specified size range, and `df` is a `DataFrame` with the sizes,
image indices, and centroids of the found objects.
"""
function sizerange(img::AbstractImage, scale, min, max; thresh=0.5, tone=true)
    distpp = approxscale!(img,scale)
    a = convert(Image{Images.Gray}, img)
    a[a.<thresh] = 0.0
    a[a.>=thresh] = 1.0
    tone && (a = 1.0 - a)
    b = convert(Array{Bool}, Images.data(a))
    arr = label_components(b)
    df = DataFrames.DataFrame(len=component_lengths(arr).*(distpp^2), ind=component_indices(arr),
        cen=component_centroids(arr))
    df[:cen] = map(x->(x[1]*distpp, x[2]*distpp), df[:cen])
    sort!(df, cols=[:len])
    similar(arr, Bool)
    inds = reduce(vcat, df[min .<= df[:len] .<= max, :ind])
    c = similar(arr, Bool)
    c[:] = false
    c[inds] = true
    im = grayim(convert(Array{Float64}, c))
    im["spacedirections"] = img["spacedirections"]
    im, df[min .<= df[:len] .<= max, :]
end

"""
```
colrow(df::DataFrame, cols, rows)
```

Uses some clustering to determine rows and columns of found objects in `df`.
You should tell it the number of columns and rows for the correct clustering.
"""
function colrow(df::DataFrame, cols, rows)
    xlen = Array{Float64}(size(df)[1])
    ylen = Array{Float64}(size(df)[1])
    for (i,v) in enumerate(df[:cen])
        x,y = v
        xlen[i] = x
        ylen[i] = y
    end
    xcen = kmeans(Array{Float64}(transpose(xlen)), cols).centers
    ycen = kmeans(Array{Float64}(transpose(ylen)), rows).centers
    sort!(reshape(transpose(xcen), (length(xcen),))),
        sort!(reshape(transpose(ycen), (length(ycen),)))
    # transpose(xcen), transpose(ycen)
end

# img has correct scaling set
"""
```
groundplane(img::AbstractImage, cols_or_rows, dim, spacing, noiselen)
```

Given some image without rotation with rectangular holes in a ground plane...

- `img` - a binary image (expects grayscale).
- `cols_or_rows`- typically one of the results from [`colrow`](@ref).
- `dim` - 1 or 2 (should be 1 for col, 2 for row)
- `spacing` - The design spacing between hole centers in the ground plane.
- `noiselen` - Contiguous lengths of 1 or 0 shorter than this are ignored in
length determination.
"""
function groundplane(img::AbstractImage, cols_or_rows, dim, spacing, noiselen)
    sd = img["spacedirections"]
    distpp = sd[1]
    a = convert(Image{Images.Gray}, img)
    b = convert(Array{Bool}, Images.data(a))

    pixlengths = Array{Array{Int,1},1}(length(cols_or_rows))
    # pixlengthsrow = Array{Array{Int,1},1}(length(rows))

    for (i,x) in enumerate(cols_or_rows)
        ctr = 1
        flag = false
        pix = Array{Int}(0)
        pixlengths[i] = Array{Int}(0)
        cut = (dim == 1 ? b[Int(round(x/distpp)),:] : b[:, Int(round(x/distpp))])
        last = cut[1]
        for z in cut[2:end]
            if z != last
                push!(pix, ctr)
                ctr = 0
            end
            ctr += 1
            last = z
        end
        ctr != 1 && push!(pix, ctr)

        minpix = round(noiselen/distpp)
        runninglength = 0
        for v in pix
            runninglength += v
            if runninglength > minpix
                push!(pixlengths[i], runninglength)
                runninglength = 0
            end
        end
    end

    l = 0.0
    for arr in pixlengths
        deleteat!(arr, 1)
        pop!(arr)
        mod(length(arr), 2) == 1 && pop!(arr)
        l += reduce(+, arr) / (length(arr) / 2)
    end
    l /= length(pixlengths)
    # Now l is the average length in pixels
    distpp = spacing/l

    display(pixlengths)
    physlengths = reduce(hcat, [map(x->x*distpp, arr) for arr in pixlengths])
    # distpp, physlengths
end

end # module


text:
module NestedViews

export combine, slices

using ArgCheck: @argcheck
using Base: @propagate_inbounds
using DocStringExtensions: SIGNATURES

"""
T, innersize = $(SIGNATURES)

Use the element type and size from the first array.
"""
function use_first(A)
    f = first(A)
    eltype(f), size(f)
end

struct CombinedView{T,N,S,P} <: AbstractArray{T,N}
    innersize::S
    parent::P
    function CombinedView{T,N}(innersize::S, parent::P
                               ) where {T,N,S<:Tuple{Vararg{Int}},P<:AbstractArray}
        @argcheck N isa Int && N ≥ 0
        @argcheck all(innersize .> 0)
        new{T,N,S,P}(innersize, parent)
    end
end

"""
$(SIGNATURES)

Return an `AbstractArray` which is a view that combines the elements of `A`, in the sense
that `combine(A)[i..., j...] ≡ A[i...][j...]` for indices that are compatible.

`common` is a function of `A` that is used to determine the common element type and size.
"""
function combine(A::AbstractArray; common = use_first)
    T, innersize = common(A)
    CombinedView{T, ndims(A) + length(innersize)}(innersize, A)
end

Base.parent(c::CombinedView) = c.parent

Base.IndexStyle(::Type{<:CombinedView}) = Base.IndexCartesian()

Base.size(c::CombinedView) = (size(parent(c))..., c.innersize...)

@propagate_inbounds function Base.getindex(c::CombinedView{T,N,S},
                                           I::Vararg{Int,N}) where {T,N,S}
    K = length(c.innersize)
    M = N - K
    outer_ix = ntuple(i -> I[i], Val{M}())
    inner_ix = ntuple(i -> I[M + i], Val{K}())
    (parent(c)[outer_ix...]::AbstractArray{T,K})[inner_ix...]
end

struct SlicedView{T,N,P<:AbstractArray} <: AbstractArray{T,N}
    parent::P
end

@propagate_inbounds @inline function _slices_view(A, I, ::Val{N}) where N
    view(A, I..., ntuple(_ -> Colon(), Val{ndims(A) - N}())...)
end

"""
$(SIGNATURES)

Return an `AbstractArray` which is a view that contains slices of `A` following the index
`N`, ie `slices(A, Val(K))[i...][j...] ≡ A[i..., j...]` for indices that are compatible,
particularly `length(i) ≡ N`.
"""
function slices(A::AbstractArray{T,M}, ::Val{N}) where {T,M,N}
    @argcheck 0 ≤ N ≤ M
    # only the type is neeed, let's hope the compiler optimizes this out
    elt = _slices_view(A, ntuple(i -> firstindex(A, i), Val{N}()), Val{N}())
    SlicedView{typeof(elt),N,typeof(A)}(A)
end

slices(A::AbstractArray, N::Integer) = slices(A, Val{N}())

Base.parent(s::SlicedView) = s.parent

Base.IndexStyle(::Type{<:SlicedView}) = Base.IndexCartesian()

Base.size(s::SlicedView{T,N}) where {T,N} = size(parent(s))[1:N]

@propagate_inbounds function Base.getindex(s::SlicedView{T,N},
                                                I::Vararg{Int,N}) where {T,N}
    _slices_view(parent(s), I, Val{N}())
end

end # module


text:
# ---
# title: 167. Two Sum II - Input array is sorted
# id: problem167
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Two Pointers, Binary Search
# link: <https://leetcode.com/problems/two-sum-ii-input-array-is-sorted/description/>
# hidden: true
# ---
#
# Given an array of integers that is already **_sorted in ascending order_** ,
# find two numbers such that they add up to a specific target number.
#
# The function twoSum should return indices of the two numbers such that they
# add up to the target, where index1 must be less than index2.
#
# **Note:**
#
#   * Your returned answers (both index1 and index2) are not zero-based.
#   * You may assume that each input would have _exactly_ one solution and you may not use the _same_ element twice.
#
#
#
# **Example 1:**
#
#
#
#     Input: numbers = [2,7,11,15], target = 9
#     Output: [1,2]
#     Explanation: The sum of 2 and 7 is 9. Therefore index1 = 1, index2 = 2.
#
#
# **Example 2:**
#
#
#
#     Input: numbers = [2,3,4], target = 6
#     Output: [1,3]
#
#
# **Example 3:**
#
#
#
#     Input: numbers = [-1,0], target = -1
#     Output: [1,2]
#
#
#
#
# **Constraints:**
#
#   * `2 <= nums.length <= 3 * 104`
#   * `-1000 <= nums[i] <= 1000`
#   * `nums` is sorted in **increasing order**.
#   * `-1000 <= target <= 1000`
#
#
## @lc code=start
using LeetCode

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


text:
#=

# 107: 1D Nonlinear Storage
([source code](SOURCE_URL))

This equation comes from the transformation of the nonlinear diffuision equation.
```math
\partial_t u^\frac{1}{m} -\Delta u = 0
```
in $\Omega=(-1,1)$ with homogeneous Neumann boundary conditions.
We can derive an exact solution from the Barenblatt solution of the previous
example.

=#

module Example107_NonlinearStorage1D
using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize

function barenblatt(x,t,m)
    tx=t^(-1.0/(m+1.0))
    xx=x*tx
    xx=xx*xx
    xx=1- xx*(m-1)/(2.0*m*(m+1));
    if xx<0.0
        xx=0.0
    end
    return tx*xx^(1.0/(m-1.0))
end

function main(;n=20,m=2.0,Plotter=nothing,verbose=false, unknown_storage=:sparse,tend=0.01,tstep=0.0001)

    ## Create a one-dimensional discretization
    h=1.0/convert(Float64,n/2)
    X=collect(-1:h:1)
    grid=VoronoiFVM.Grid(X)

    ## Flux function which describes the flux
    ## between neigboring control volumes
    function flux!(f,u0,edge)
        u=unknowns(edge,u0)
        f[1]=u[1,1]-u[1,2]
    end

    ϵ=1.0e-10

    ## Storage term
    ## This needs to be regularized as its derivative
    ## at 0 is infinity
    function storage!(f,u,node)
        f[1]=(ϵ+u[1])^(1.0/m)
    end

    ## Create a physics structure
    physics=VoronoiFVM.Physics(
        flux=flux!,
        storage=storage!)

    ## Create a finite volume system - either
    ## in the dense or  the sparse version.
    ## The difference is in the way the solution object
    ## is stored - as dense or as sparse matrix
    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)

    ## Add species 1 to region 1
    enable_species!(sys,1,[1])

    ## Create a solution array
    inival=unknowns(sys)
    solution=unknowns(sys)
    t0=0.001

    ## Broadcast the initial value
    inival[1,:].=map(x->barenblatt(x,t0,m)^m,X)

    ## Create solver control info
    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose
    control.Δu_opt=0.1
    control.force_first_step=true
    tsol=solve(inival,sys,[t0,tend],control=control)

    p=GridVisualizer(Plotter=Plotter,layout=(1,1),fast=true)
    for i=1:length(tsol)
        time=tsol.t[i]
        scalarplot!(p[1,1],grid,tsol[1,:,i],title=@sprintf("t=%.3g",time),color=:red,label="numerical")
        scalarplot!(p[1,1],grid,map(x->barenblatt(x,time,m)^m,grid),clear=false,color=:green,label="exact")
        reveal(p)
        sleep(1.0e-2)
    end

    return sum(tsol[end])
end

function test()
    testval=175.20261258406686
    main(unknown_storage=:sparse) ≈ testval && main(unknown_storage=:dense) ≈ testval
end

end


text:
@doc raw"""
    has_approx_invariant_metric(
        G::AbstractDecoratorManifold,
        p,
        X,
        Y,
        qs::AbstractVector,
        conv::ActionDirection = LeftAction();
        kwargs...,
    ) -> Bool

Check whether the metric on the group $\mathcal{G}$ is (approximately) invariant using a set of predefined
points. Namely, for $p ∈ \mathcal{G}$, $X,Y ∈ T_p \mathcal{G}$, a metric $g$, and a
translation map $τ_q$ in the specified direction, check for each $q ∈ \mathcal{G}$ that the
following condition holds:

````math
g_p(X, Y) ≈ g_{τ_q p}((\mathrm{d}τ_q)_p X, (\mathrm{d}τ_q)_p Y).
````

This is necessary but not sufficient for invariance.

Optionally, `kwargs` passed to `isapprox` may be provided.
"""
has_approx_invariant_metric(::AbstractDecoratorManifold, p, X, Y, qs, ::ActionDirection)
@trait_function has_approx_invariant_metric(
    M::AbstractDecoratorManifold,
    p,
    X,
    Y,
    qs,
    conv::ActionDirection=LeftAction();
    kwargs...,
)
function has_approx_invariant_metric(
    ::TraitList{<:IsGroupManifold},
    M::AbstractDecoratorManifold,
    p,
    X,
    Y,
    qs,
    conv::ActionDirection;
    kwargs...,
)
    gpXY = inner(M, p, X, Y)
    for q in qs
        τq = translate(M, q, p, conv)
        dτqX = translate_diff(M, q, p, X, conv)
        dτqY = translate_diff(M, q, p, Y, conv)
        isapprox(gpXY, inner(M, τq, dτqX, dτqY); kwargs...) || return false
    end
    return true
end

"""
    direction(::AbstractDecoratorManifold) -> AD

Get the direction of the action a certain Lie group with its implicit metric has
"""
direction(::AbstractDecoratorManifold)

@trait_function direction(M::AbstractDecoratorManifold)

direction(::TraitList{HasLeftInvariantMetric}, ::AbstractDecoratorManifold) = LeftAction()

direction(::TraitList{HasRightInvariantMetric}, ::AbstractDecoratorManifold) = RightAction()

function exp(::TraitList{HasLeftInvariantMetric}, M::AbstractDecoratorManifold, p, X)
    return retract(M, p, X, GroupExponentialRetraction(LeftAction()))
end
function exp!(::TraitList{HasLeftInvariantMetric}, M::AbstractDecoratorManifold, q, p, X)
    return retract!(M, q, p, X, GroupExponentialRetraction(LeftAction()))
end
function exp(::TraitList{HasRightInvariantMetric}, M::AbstractDecoratorManifold, p, X)
    return retract(M, p, X, GroupExponentialRetraction(RightAction()))
end
function exp!(::TraitList{HasRightInvariantMetric}, M::AbstractDecoratorManifold, q, p, X)
    return retract!(M, q, p, X, GroupExponentialRetraction(RightAction()))
end
function exp(::TraitList{HasBiinvariantMetric}, M::MetricManifold, p, X)
    return exp(M.manifold, p, X)
end
function exp!(::TraitList{HasBiinvariantMetric}, M::MetricManifold, q, p, X)
    return exp!(M.manifold, q, p, X)
end

@trait_function has_invariant_metric(M::AbstractDecoratorManifold, op::ActionDirection)

# Fallbacks / false
has_invariant_metric(::AbstractManifold, op::ActionDirection) = false
function has_invariant_metric(
    ::TraitList{<:HasLeftInvariantMetric},
    ::AbstractDecoratorManifold,
    ::LeftAction,
)
    return true
end
function has_invariant_metric(
    ::TraitList{<:HasRightInvariantMetric},
    ::AbstractDecoratorManifold,
    ::RightAction,
)
    return true
end

@trait_function has_biinvariant_metric(M::AbstractDecoratorManifold)

# fallbavk / default: false
has_biinvariant_metric(::AbstractManifold) = false
function has_biinvariant_metric(
    ::TraitList{<:HasBiinvariantMetric},
    ::AbstractDecoratorManifold,
)
    return true
end

function inner(
    t::TraitList{IT},
    M::AbstractDecoratorManifold,
    p,
    X,
    Y,
) where {IT<:AbstractInvarianceTrait}
    conv = direction(t, M)
    Xₑ = inverse_translate_diff(M, p, p, X, conv)
    Yₑ = inverse_translate_diff(M, p, p, Y, conv)
    return inner(next_trait(t), M, Identity(M), Xₑ, Yₑ)
end
function inner(
    t::TraitList{<:IsGroupManifold},
    M::AbstractDecoratorManifold,
    ::Identity,
    X,
    Y,
)
    return inner(next_trait(t), M, identity_element(M), X, Y)
end

function inverse_translate_diff(
    ::TraitList{IsMetricManifold},
    M::MetricManifold,
    p,
    q,
    X,
    conv::ActionDirection,
)
    return inverse_translate_diff(M.manifold, p, q, X, conv)
end
function inverse_translate_diff!(
    ::TraitList{IsMetricManifold},
    M::MetricManifold,
    Y,
    p,
    q,
    X,
    conv::ActionDirection,
)
    return inverse_translate_diff!(M.manifold, Y, p, q, X, conv)
end

function log(::TraitList{HasLeftInvariantMetric}, M::AbstractDecoratorManifold, p, q)
    return inverse_retract(M, p, q, GroupLogarithmicInverseRetraction(LeftAction()))
end
function log!(::TraitList{HasLeftInvariantMetric}, M::AbstractDecoratorManifold, X, p, q)
    return inverse_retract!(M, X, p, q, GroupLogarithmicInverseRetraction(LeftAction()))
end
function log(::TraitList{HasRightInvariantMetric}, M::AbstractDecoratorManifold, p, q)
    return inverse_retract(M, p, q, GroupLogarithmicInverseRetraction(RightAction()))
end
function log!(::TraitList{HasRightInvariantMetric}, M::AbstractDecoratorManifold, X, p, q)
    return inverse_retract!(M, X, p, q, GroupLogarithmicInverseRetraction(RightAction()))
end
function log(::TraitList{HasBiinvariantMetric}, M::MetricManifold, p, q)
    return log(M.manifold, p, q)
end
function log!(::TraitList{HasBiinvariantMetric}, M::MetricManifold, X, p, q)
    return log!(M.manifold, X, p, q)
end

function LinearAlgebra.norm(
    t::TraitList{IT},
    M::AbstractDecoratorManifold,
    p,
    X,
) where {IT<:AbstractInvarianceTrait}
    conv = direction(t, M)
    Xₑ = inverse_translate_diff(M, p, p, X, conv)
    return norm(next_trait(t), M, Identity(M), Xₑ)
end
function LinearAlgebra.norm(
    t::TraitList{<:IsGroupManifold},
    M::AbstractDecoratorManifold,
    ::Identity,
    X,
)
    return norm(next_trait(t), M, identity_element(M), X)
end

function translate_diff!(
    ::TraitList{IsMetricManifold},
    M::MetricManifold,
    p,
    q,
    X,
    conv::ActionDirection,
)
    return translate_diff(M.manifold, p, q, X, conv)
end
function translate_diff!(
    ::TraitList{IsMetricManifold},
    M::MetricManifold,
    Y,
    p,
    q,
    X,
    conv::ActionDirection,
)
    return translate_diff!(M.manifold, Y, p, q, X, conv)
end


text:
using Distributed
using exactQmodule
using Arpack
using MPOmodule2
using finiteMPS
#@everywhere using ProgressMeter
@everywhere using DelimitedFiles
@everywhere include("finiteTDVP_SSE.jl")
@everywhere include("MPOcollection.jl")
include("PauliComutator.jl")
include("decomposeJump.jl")
include("dePara.jl")

#
#	sys size
N = 16
Np = 4

V = 1.5
U = 0.5
J = 1.0 #intra
J′ = 0.5 #inter
μ = 0.0

#
#	physical units
#
GHz = 1
MHz = 0.001
THz = 1000
ns  = 1
mus = 1/MHz
ps  = 1/THz
K   = 20.837
Tesla  = 1

#
#	Bath spectral function
#
Temp = 0.1*K
γ = 20*MHz
Λ = 30*GHz
df = 100*MHz
freqvec = -20*Λ:df:20*Λ

function J_CutOff(ω, Λ, T)
	res = ω == 0 ?  T/GHz : (ω/(1*GHz))*exp(-ω^2/(2*Λ^2)) / (1-exp(-ω/T))
	return res
end

function S(x) return  J_CutOff(x, Λ, Temp) end
g(x) = sqrt(S(x)/(2*pi))
G_time, time_vec = get_time_domain_func(g, freqvec, sorted = true) # get g in time

"""**************************** run SSE-TDVP ***********************************"""
N_realize = 30
dτ = 0.01
dmax = 256
order = 4

#load instance
A = loadMPS_SSE("./MPS_SSE_instance_N$(N)_T$(Temp)_gamma$(γ)_order$(order)_U$(U).jld")

println("start evolution")

# creat dir to save runs!
dir_name = "./results/N$(N)_Np$(Np)_Temp$(Temp)_γ$(γ)_dτ$(dτ)_U$(U)_order$(order)_dmax$(dmax)"
try
	mkdir(dir_name)
	mkdir(string(dir_name,"/inter_states"))
catch
	println(string("directory",dir_name," already exists."))
end

@everywhere global f = open("error_message.txt", "w+")
@everywhere redirect_stderr(f)

f = open("datalist.txt", "r")
f = readlines(f)

#=
mps = load(f[1])["MPSvec"]
i = findlast("id", f[1])
x = f[1][i[2]+1:end]
i = findfirst(".", x)
x = x[1:i[1]-1]
@show x
name = (dir_name, x)
x = MPS_SSE(mps, A.H_MPO, A.Heff, A.Decay_MPO[1:1:end], A.Jump_MPO[1:1:end], 0.025, (250,260), 150, A.measureDict)
evo_TDVP_SSE(x, MPOvec = A.Heff, dτ = dτ, name = name)
=#

pmap([1:length(f)...]) do x
	mps = load(f[x])["MPSvec"]

	i = findlast("id", f[x])
	x = f[x][i[2]+1:end]
	i = findfirst(".", x)
	x = x[1:i[1]-1]
	@show x
	name = (dir_name, string(x,"_continued"))
	x = MPS_SSE(mps, A.H_MPO, A.Heff, A.Decay_MPO[1:1:end], A.Jump_MPO[1:1:end], 0.025, (150,170), dmax, A.measureDict)
	evo_TDVP_SSE(x, MPOvec = A.Heff, dτ = dτ, name = name, saveState = true)
end


text:
# Julia script to compute robsutness coefficients -
include("Sample.jl")
using PyCall
@pyimport numpy as np

const data_plasma = readdlm("../data/Average_Plasma_05.csv")
const dilution_factor = 0.79
const mixing_delay_time = 490.0
const tpa_level = 4.0
const experimental_data_array = [data_plasma[:,1] data_plasma[:,round(Int,tpa_level)+2]]
const rank_threshold = 1
const percentage_coverage = 0.60
const SMALL = 1e-6

# main: calculates the ratio of the AUC in the perturbed state, to the base state
function main()

  # Calculate the base AUC -
  base_auc_array = calculate_auc_base()
  idx_small = find(base_auc_array.< SMALL)
  base_auc_array[idx_small] = SMALL

  # how many sets did we sample?
  number_of_samples = length(base_auc_array)

  # Calculate the perturbed case -
  number_of_divisions = 50
  plasminogen_array = collect(linspace(0,3,number_of_divisions))
  tafi_array = collect(linspace(0,100,number_of_divisions))

  # Calculate the perturbed cases -
  mean_robsutness_array = -1.0*ones(number_of_divisions,number_of_divisions)
  std_robustness_array = -1.0*ones(number_of_divisions,number_of_divisions)
  for (plasminogen_array_index,plasminogen_value) in enumerate(plasminogen_array)
    for (tafi_array_index,tafi_value) in enumerate(tafi_array)

      # Solves the model over the ensemble, calculates the AUC and returns back the array -
      pertubed_auc_array = calculate_auc_perturbed(plasminogen_value,tafi_value)

      # calculate the robsutness coeffcients -
      tmp_array = zeros(number_of_samples)
      for sample_index = 1:number_of_samples

        base_value = base_auc_array[sample_index]
        perturbed_value = pertubed_auc_array[sample_index]
        tmp_array[sample_index] = perturbed_value/base_value # calculate robustness
      end

      # calc the mean and std -
      mean_robustness_value = mean(tmp_array)
      std_robustness_value = std(tmp_array)

      # cache -
      mean_robsutness_array[plasminogen_array_index,tafi_array_index] = mean_robustness_value
      std_robustness_array[plasminogen_array_index,tafi_array_index] = std_robustness_value

      println("Plasminogen index: "*string(plasminogen_array_index)*" and TAFI index: "*string(tafi_array_index))
    end

    # save tmp results to disk -
    tmp_mean_file_path = "./robustness/tmp-mean.dat"
    writedlm(tmp_mean_file_path,mean_robsutness_array)
    tmp_std_file_path = "./robustness/tmp-std.dat"
    writedlm(tmp_std_file_path,std_robustness_array)
  end

  return (mean_robsutness_array,std_robustness_array,plasminogen_array,tafi_array)
end

function calculate_auc_base()

  # Calculate the base AUC -

  # setup plasma specific constants -
  plasminogen_mean_level = 1.5
  tafi_mean_level = 50.0

  # Run the model -
  (T,A,mean_array,std_array,sample_idx_vec) = sample_ensemble(dilution_factor::Float64,tpa_level::Float64,
                                                                plasminogen_mean_level::Float64,
                                                                tafi_mean_level::Float64,
                                                                rank_threshold::Int,
                                                                percentage_coverage::Float64,
                                                                mixing_delay_time::Float64,
                                                                experimental_data_array::Array{Float64,2})

  # Calculate the AUC for the simulation archive A -
  (number_of_rows,number_of_cols) = size(A)
  auc_array = zeros(number_of_cols)
  for col_index = 1:number_of_cols
    auc_array[col_index] = trapz(T,A[:,col_index]) # integration
  end

  return (auc_array)
end

function calculate_auc_perturbed(plasminogen_mean_level,tafi_mean_level)

  # Run the model -
  (T,A,mean_array,std_array,sample_idx_vec) = sample_ensemble(dilution_factor::Float64,tpa_level::Float64,
                                                                plasminogen_mean_level::Float64,
                                                                tafi_mean_level::Float64,
                                                                rank_threshold::Int,
                                                                percentage_coverage::Float64,
                                                                mixing_delay_time::Float64,
                                                                experimental_data_array::Array{Float64,2})

  # Calculate the AUC for the simulation archive A -
  (number_of_rows,number_of_cols) = size(A)
  auc_array = zeros(number_of_cols)
  for col_index = 1:number_of_cols
    auc_array[col_index] = trapz(T,A[:,col_index])
  end

  return (auc_array)
end

function trapz{Tx<:Number, Ty<:Number}(x::Vector{Tx}, y::Vector{Ty})
    # Trapezoidal integration rule
    local n = length(x)
    if (length(y) != n)
        error("Vectors 'x', 'y' must be of same length")
    end
    r = zero(zero(Tx) + zero(Ty))
    if n == 1; return r; end
    for i in 2:n
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    #= correction -h^2/12 * (f'(b) - f'(a))
    ha = x[2] - x[1]
    he = x[end] - x[end-1]
    ra = (y[2] - y[1]) / ha
    re = (y[end] - y[end-1]) / he
    r/2 - ha*he/12 * (re - ra)
    =#
    return r/2
end


text:
function predict_svm(θ,x)
	return dot(θ,x) > 0 ? 1 : -1
end

function predict_logreg(theta,x)
	prediction =  1 / ( 1 + exp(dot(-theta,x)))
	if prediction  < 0.5
		return -1
	else
		return 1
	end
end

function evaluate(X,y,theta,mode,verbose=false)
	positiveHits = 0
	for i in 1:size(X,1)
		if mode == "logreg"
			label = predict_logreg(theta,X[i:i,:])
		elseif mode == "svm"
			label = predict_svm(theta,X[i:i,:])
		else
			error("Prediction mode not aviable")
		end
		if label == y[i]
			positiveHits += 1
		end
		if label == 1
			label = "+1"
		end
		labelString = "-1.0"
		if y[i] == 1.0
			labelString = "+1.0"
		end
		if verbose
			@printf("Prediction : %s | Label : %s | Accuracy %f \n",label,labelString,positiveHits / i)
		end
	end
	return positiveHits / size(X,1)
end

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
using SparseArrays, LinearAlgebra

X =  sparse([1,2,3,4,5,6,7,8,1,3,4,5,2,6,7,8], [1,1,2,2,3,3,4,4,5,5,5,5,6,6,6,6],[1,1,1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,-1,-1])
Xprime = sparse([1,2,3,4,5,6,7,8,1,3,4,5], [1,1,2,2,3,3,4,4,5,5,5,5],[1,1,1,1,1,1,1,1,-1,-1,-1,-1])

settings_default = Settings()
settings_direct = Settings(lls_algorithm = DirectLLS())
settings_CMG = Settings(lls_algorithm = CMGPreconditionedLLS())

@testset "solver" begin
    #Direct Method
    @test lss(settings_default.lls_algorithm, X, X[1,:], settings_default)  ≈ [ 0.5166666666666669 ; -0.23333333333333325 ; 0.01666666666666683; 0.2666666666666666; -0.23333333333333314; 0.26666666666666666]

    #AMGPreconditionedLLS
    @test lss(settings_direct.lls_algorithm, Xprime, Xprime[1,:], settings_direct) ≈ [0.25; -0.5; -0.25; 0.0; -0.5]
    @test_throws SingularException(6) lss(settings_direct.lls_algorithm, X, X[1,:], settings_direct)

    #CMGPreconditionedLLS (using MATLAB; don't run them in CI)
    # @test_throws MATLAB.MEngineError("failed to get variable iter from MATLAB session") lss(settings_CMG.lls_algorithm, Xprime, Xprime[1,:], settings_CMG)  # ≈ [0.25; -0.5; -0.25; 0.0; -0.5]
    # @test_throws MATLAB.MEngineError("failed to get variable iter from MATLAB session") lss(settings_CMG.lls_algorithm, X, X[1,:], settings_CMG)  # ≈ [0.25; -0.5; -0.25; 0.0; -0.5; 0.0]
end


text:
module TestDictoVec
using Test
using RData

function show2string(dv::DictoVec)
    iob = IOBuffer()
    show(iob, dv)
    return String(take!(iob))
end

@testset "RData DictoVec type" begin

@testset "Empty DictoVec" begin
    @test_throws MethodError DictoVec(Symbol[], [:a])
    @test_throws DimensionMismatch DictoVec(Symbol[], ["a"])
    dv = DictoVec(Symbol[])

    @test typeof(dv) === DictoVec{Symbol}
    @test eltype(dv) === Symbol
    @test length(dv) == 0
    @test isempty(dv)
    @test collect(keys(dv)) == RData.RString[]
    @test values(dv) == Symbol[]
    @test !haskey(dv, "a")
    @test !haskey(dv, :a)
    @test !haskey(dv, 1)
    @test_throws BoundsError dv[1]
    @test_throws KeyError dv["a"]
    @test_throws KeyError dv[:a]

    @test get(dv, 1, :x) == :x
    @test get(() -> :y, dv, 1) == :y
    @test get(dv, "a", :x) == :x
    @test get(() -> :y, dv, "a") == :y
    @test length(dv) == 0

    @test convert(Dict{RData.RString, Any}, dv) == Dict{RData.RString, Any}()

    @test show2string(dv) == "DictoVec{Symbol}()"

    # add an element
    @test_throws MethodError (dv[1] = "xx") # incompatible value
    @test length(dv) == 0
    @test_throws BoundsError (dv[1] = :xx) # cannot create new elements by integer index
    @test length(dv) == 0
    @test_throws MethodError (dv[:a] = "xx") # incompatible key
    @test length(dv) == 0
    @test_throws MethodError (dv["a"] = "xx") # incompatible value
    @test length(dv) == 0
    dv["a"] = :xx
    @test !isempty(dv)
    @test length(dv) == 1
    @test haskey(dv, "a")
    @test dv["a"] == :xx
    @test dv[1] == :xx
    @test collect(keys(dv)) == ["a"]
    @test values(dv) == [:xx]

    @test show2string(dv) == "DictoVec{Symbol}(\"a\"=>:xx)"

    # reassign element
    dv[1] = :yy
    @test length(dv) == 1
    @test dv[1] == :yy
    @test dv["a"] == :yy

    # deleting an element by key
    @test delete!(dv, 1) === dv
    @test length(dv) == 1
    delete!(dv, "b")
    @test length(dv) == 1
    delete!(dv, "a")
    @test length(dv) == 0
    @test !haskey(dv, "a")
    @test_throws BoundsError dv[1]
end

@testset "Nameless DictoVec" begin
    dv = DictoVec([2.0, 5.0, 4.0])
    @test typeof(dv) === DictoVec{Float64}
    @test eltype(dv) === Float64
    @test length(dv) == 3
    @test !isempty(dv)
    @test !haskey(dv, 1)
    @test !haskey(dv, "a")
    @test !haskey(dv, :a)
    @test show2string(dv) == "DictoVec{Float64}(2.0,5.0,4.0)"
    @test collect(keys(dv)) == RData.RString[]
    @test values(dv) == [2.0, 5.0, 4.0]

    @test_throws BoundsError dv[0]
    @test_throws BoundsError dv[4]
    @test dv[1] == 2.0
    @test dv[[1, 3]] == [2.0, 4.0]
    @test_throws KeyError dv["a"]
    @test_throws KeyError dv[:a]

    # add new element by name
    dv["a"] = 3
    @test haskey(dv, "a")
    @test length(dv) == 4
    @test dv["a"] === 3.0
    @test dv[4] === 3.0
    @test show2string(dv) == "DictoVec{Float64}(2.0,5.0,4.0,\"a\"=>3.0)"

    @test delete!(dv, "a") === dv
    @test length(dv) == 3
    @test !haskey(dv, "a")
end

@testset "DictoVec with names" begin
    dv = DictoVec([2.0, 5.0, 4.0], ["a", "b", "c"])
    @test typeof(dv) === DictoVec{Float64}
    @test eltype(dv) === Float64
    @test length(dv) == 3
    @test !isempty(dv)
    @test !haskey(dv, 1)
    @test haskey(dv, "a")
    @test !haskey(dv, :a)
    @test sort!(collect(keys(dv))) == ["a", "b", "c"]
    @test values(dv) == [2.0, 5.0, 4.0]
    @test show2string(dv) == "DictoVec{Float64}(\"a\"=>2.0,\"b\"=>5.0,\"c\"=>4.0)"

    @test dv[1] === 2.0
    @test dv["a"] === 2.0
    @test dv[[1, 3]] == [2.0, 4.0]
    @test_throws KeyError dv[:a]

    # reassign element by name
    dv["a"] = 6
    @test haskey(dv, "a")
    @test length(dv) == 3
    @test dv["a"] === 6.0
    @test dv[1] === 6.0
    @test show2string(dv) == "DictoVec{Float64}(\"a\"=>6.0,\"b\"=>5.0,\"c\"=>4.0)"

    # deleting an element from the middle
    @test delete!(dv, "b") === dv
    @test length(dv) == 2
    @test !haskey(dv, "b")
    @test dv[2] == 4.0 # indices has updated
    @test show2string(dv) == "DictoVec{Float64}(\"a\"=>6.0,\"c\"=>4.0)"
end

end

end # TestDictoVec


text:
#=---------------------------------------------------------------
Analysis Stage 2: For each RNN type, fit a GAM to predict
the log-likelihood ratio (and effect direction) from the weighted
average surprisals. Create data frames with all values for plotting.
---------------------------------------------------------------=#

using RCall
print("Fitting GAMs")

function do_gam(lmm_rnn)
    # Fit GAM in R, generate a plot (without actually plotting it), and return plot values

    model = R"""m = gam(effect ~ s(logprob)+s(logprob,rep,bs="re"), data=$lmm_rnn)"""
    plot  = R"plot($model, select=0)"

    intercept = model[:coefficients][1]
    fit = convert(Array{Float32,1}, plot[1][:fit])
    se  = convert(Array{Float32,1}, plot[1][:se])
    x   = convert(Array{Float32,1}, plot[1][:x])
    (intercept, fit, se, x)
end

function get_gam_curve(lmm_results)
    lmm_results[!,:effect]  = lmm_results[!,:chi2] .* sign.(lmm_results[!,:b_surp])                  # reverse effect direction if coefficient negative. What if b_surp and b_prevsurp of unequal sign?!
    lmm_results[!,:logprob] = -lmm_results[!,:avsurp]

    lmm_srn  = lmm_results[lmm_results[!,:rnntype].=="SRN",:]
    lmm_gru  = lmm_results[lmm_results[!,:rnntype].=="GRU",:]
    lmm_lstm = lmm_results[lmm_results[!,:rnntype].=="LSTM",:]

    (intercept_srn,  fit_srn,  se_srn,  x_srn)  = do_gam(lmm_srn)
    (intercept_gru,  fit_gru,  se_gru,  x_gru)  = do_gam(lmm_gru)
    (intercept_lstm, fit_lstm, se_lstm, x_lstm) = do_gam(lmm_lstm)

    gam_curve_srn  = DataFrame(rnntype="SRN",  fit=intercept_srn .+ fit_srn,  se=se_srn,  x=x_srn)
    gam_curve_gru  = DataFrame(rnntype="GRU",  fit=intercept_gru .+ fit_gru,  se=se_gru,  x=x_gru)
    gam_curve_lstm = DataFrame(rnntype="LSTM", fit=intercept_lstm .+ fit_lstm, se=se_lstm, x=x_lstm)

    [gam_curve_srn; gam_curve_gru; gam_curve_lstm]
end

R"library(mgcv)"

gam_curve_SPR = get_gam_curve(CSV.read("lmm_SPR.csv", delim='\t'))
CSV.write("gam_SPR.csv", gam_curve_SPR; delim='\t')

gam_curve_ET = get_gam_curve(CSV.read("lmm_ET.csv", delim='\t'))
CSV.write("gam_ET.csv", gam_curve_ET; delim='\t')

gam_curve_EEG = get_gam_curve(CSV.read("lmm_EEG.csv", delim='\t'))
CSV.write("gam_EEG.csv", gam_curve_EEG; delim='\t')

max_stars_repo_name: caurnhammer/gated-rnns

max_issues_repo_name: caurnhammer/gated-rnns

max_forks_repo_name: caurnhammer/gated-rnns


text:
function whichUInt(n::Integer)
    n == 8 && return UInt8
    n == 16 && return UInt16
    n == 32 && return UInt32
    n == 64 && return UInt64
    throw(error("Only n=8,16,32,64 supported."))
end

whichUInt(::Type{T}) where T = whichUInt(sizeof(T)*8)