text:
"""
    predict_kernel!
"""
function predict_kernel!(::L, pred::AbstractMatrix{T}, split, feat, cond_float, leaf_pred::AbstractMatrix{T}, X::CuDeviceMatrix, K) where {L,T}
    idx = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    nid = 1
    @inbounds if idx <= size(pred, 2)
        while split[nid]
            X[idx, feat[nid]] < cond_float[nid] ? nid = nid << 1 : nid = nid << 1 + 1
        end
        @inbounds for k in 1:K
            pred[k, idx] += leaf_pred[k, nid]
        end
    end
    sync_threads()
    return nothing
end

"""
    predict_kernel!
        GradientRegression
"""
function predict_kernel!(::L, pred::AbstractMatrix{T}, split, feat, cond_float, leaf_pred::AbstractMatrix{T}, X::CuDeviceMatrix, K) where {L<:GradientRegression,T}
    idx = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    nid = 1
    @inbounds if idx <= size(pred, 2)
        while split[nid]
            X[idx, feat[nid]] < cond_float[nid] ? nid = nid << 1 : nid = nid << 1 + 1
        end
        pred[1, idx] += leaf_pred[1, nid]
    end
    sync_threads()
    return nothing
end

"""
    predict_kernel!
        Logistic
"""
function predict_kernel!(::L, pred::AbstractMatrix{T}, split, feat, cond_float, leaf_pred::AbstractMatrix{T}, X::CuDeviceMatrix, K) where {L<:Logistic,T}
    idx = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    nid = 1
    @inbounds if idx <= size(pred, 2)
        while split[nid]
            X[idx, feat[nid]] < cond_float[nid] ? nid = nid << 1 : nid = nid << 1 + 1
        end
        pred[1, idx] = min(15, max(-15, pred[1,idx] + leaf_pred[1, nid]))
    end
    sync_threads()
    return nothing
end

"""
    predict_kernel!
        GaussianRegression
"""
function predict_kernel!(::L, pred::AbstractMatrix{T}, split, feat, cond_float, leaf_pred::AbstractMatrix{T}, X::CuDeviceMatrix, K) where {L<:GaussianRegression,T}
    idx = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    nid = 1
    @inbounds if idx <= size(pred, 2)
        while split[nid]
            X[idx, feat[nid]] < cond_float[nid] ? nid = nid << 1 : nid = nid << 1 + 1
        end
        pred[1, idx] += leaf_pred[1, nid]
        pred[2, idx] = max(-15, pred[2,idx] + leaf_pred[2, nid])
    end
    sync_threads()
    return nothing
end

# prediction from single tree - assign each observation to its final leaf
function predict!(loss::L, pred::AbstractMatrix{T}, tree::TreeGPU{T}, X::AbstractMatrix, K; MAX_THREADS=1024) where {L,T}
    K = size(pred, 1)
    n = size(pred, 2)
    threads = min(MAX_THREADS, n)
    blocks = ceil(Int, n / threads)
    @cuda blocks = blocks threads = threads predict_kernel!(loss, pred, tree.split, tree.feat, tree.cond_float, tree.pred, X, K)
    CUDA.synchronize()
end

# prediction from single tree - assign each observation to its final leaf
function predict(loss::L, tree::TreeGPU{T}, X::AbstractMatrix, K) where {L,T}
    pred = CUDA.zeros(T, K, size(X, 1))
    predict!(loss, pred, tree, X, K)
    return pred
end

# prediction from single tree - assign each observation to its final leaf
function predict(model::GBTreeGPU{T}, X::AbstractMatrix) where {T}
    K = size(model.trees[1].pred, 1)
    pred = CUDA.zeros(T, K, size(X, 1))
    X_gpu = CuArray(X)
    for tree in model.trees
        predict!(model.params.loss, pred, tree, X_gpu, model.K)
    end
    if typeof(model.params.loss) == Logistic
        @. pred = sigmoid(pred)
    elseif typeof(model.params.loss) == Poisson
        @. pred = exp(pred)
    elseif typeof(model.params.loss) == Gaussian
        pred[2,:] = exp.(pred[2,:])
    elseif typeof(model.params.loss) == Softmax
        pred = transpose(reshape(pred, model.K, :))
        for i in 1:size(pred, 1)
            pred[i,:] .= softmax(pred[i,:])
        end
    end
    return Array(pred')
end

# prediction in Leaf - GradientRegression
function pred_leaf_gpu!(::S, p::AbstractMatrix{T}, n, ∑::AbstractVector{T}, params::EvoTypes) where {S <: GradientRegression,T}
    p[1,n] = - params.η * ∑[1] / (∑[2] + params.λ * ∑[3])
    return nothing
end

# prediction in Leaf - Gaussian
function pred_leaf_gpu!(::S, p::AbstractMatrix{T}, n, ∑::AbstractVector{T}, params::EvoTypes) where {S <: GaussianRegression,T}
    p[1,n] = - params.η * ∑[1] / (∑[3] + params.λ * ∑[5])
    p[2,n] = - params.η * ∑[2] / (∑[4] + params.λ * ∑[5])
    return nothing
end

ext: jl


text:
#name: t-test
#description: Welch's t-test
#help-url: https://en.wikipedia.org/wiki/Welch%27s_t-test
#language: julia
#tags: demo, hide-suggestions
#sample: TSLA.csv
#input: dataframe data [Input data table]
#input: column x {type:numerical, allowNulls:false} [X axis column name]
#input: column y {type:numerical, allowNulls:false} [Y axis column name]
#output: double pValue [P-value of t-statistics]

using HypothesisTests

pValue = pvalue(UnequalVarianceTTest(disallowmissing(data[x]), disallowmissing(data[y])));

max_stars_repo_name: cesarlpb/public

max_issues_repo_name: cesarlpb/public

max_forks_repo_name: cesarlpb/public


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

const Session_Type = (;[
    Symbol("DATA") => Int32(0),
    Symbol("SCHEMA") => Int32(1),
]...)

mutable struct Session_Open_Req <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Session_Open_Req(; kwargs...)
        obj = new(meta(Session_Open_Req), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Session_Open_Req
const __meta_Session_Open_Req = Ref{ProtoMeta}()
function meta(::Type{Session_Open_Req})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Session_Open_Req)
            __meta_Session_Open_Req[] = target = ProtoMeta(Session_Open_Req)
            allflds = Pair{Symbol,Union{Type,String}}[:database => AbstractString, :_type => Int32, :options => Options]
            meta(target, Session_Open_Req, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Session_Open_Req[]
    end
end
function Base.getproperty(obj::Session_Open_Req, name::Symbol)
    if name === :database
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :_type
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :options
        return (obj.__protobuf_jl_internal_values[name])::Options
    else
        getfield(obj, name)
    end
end

mutable struct Session_Open_Res <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Session_Open_Res(; kwargs...)
        obj = new(meta(Session_Open_Res), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Session_Open_Res
const __meta_Session_Open_Res = Ref{ProtoMeta}()
function meta(::Type{Session_Open_Res})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Session_Open_Res)
            __meta_Session_Open_Res[] = target = ProtoMeta(Session_Open_Res)
            allflds = Pair{Symbol,Union{Type,String}}[:session_id => Vector{UInt8}, :server_duration_millis => Int32]
            meta(target, Session_Open_Res, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Session_Open_Res[]
    end
end
function Base.getproperty(obj::Session_Open_Res, name::Symbol)
    if name === :session_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :server_duration_millis
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct Session_Open <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Session_Open(; kwargs...)
        obj = new(meta(Session_Open), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Session_Open
const __meta_Session_Open = Ref{ProtoMeta}()
function meta(::Type{Session_Open})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Session_Open)
            __meta_Session_Open[] = target = ProtoMeta(Session_Open)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Session_Open, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Session_Open[]
    end
end

mutable struct Session_Close_Req <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Session_Close_Req(; kwargs...)
        obj = new(meta(Session_Close_Req), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Session_Close_Req
const __meta_Session_Close_Req = Ref{ProtoMeta}()
function meta(::Type{Session_Close_Req})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Session_Close_Req)
            __meta_Session_Close_Req[] = target = ProtoMeta(Session_Close_Req)
            allflds = Pair{Symbol,Union{Type,String}}[:session_id => Vector{UInt8}]
            meta(target, Session_Close_Req, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Session_Close_Req[]
    end
end
function Base.getproperty(obj::Session_Close_Req, name::Symbol)
    if name === :session_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    else
        getfield(obj, name)
    end
end

mutable struct Session_Close_Res <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Session_Close_Res(; kwargs...)
        obj = new(meta(Session_Close_Res), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Session_Close_Res
const __meta_Session_Close_Res = Ref{ProtoMeta}()
function meta(::Type{Session_Close_Res})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Session_Close_Res)
            __meta_Session_Close_Res[] = target = ProtoMeta(Session_Close_Res)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Session_Close_Res, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Session_Close_Res[]
    end
end

mutable struct Session_Close <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Session_Close(; kwargs...)
        obj = new(meta(Session_Close), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Session_Close
const __meta_Session_Close = Ref{ProtoMeta}()
function meta(::Type{Session_Close})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Session_Close)
            __meta_Session_Close[] = target = ProtoMeta(Session_Close)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Session_Close, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Session_Close[]
    end
end

mutable struct Session_Pulse_Req <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Session_Pulse_Req(; kwargs...)
        obj = new(meta(Session_Pulse_Req), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Session_Pulse_Req
const __meta_Session_Pulse_Req = Ref{ProtoMeta}()
function meta(::Type{Session_Pulse_Req})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Session_Pulse_Req)
            __meta_Session_Pulse_Req[] = target = ProtoMeta(Session_Pulse_Req)
            allflds = Pair{Symbol,Union{Type,String}}[:session_id => Vector{UInt8}]
            meta(target, Session_Pulse_Req, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Session_Pulse_Req[]
    end
end
function Base.getproperty(obj::Session_Pulse_Req, name::Symbol)
    if name === :session_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    else
        getfield(obj, name)
    end
end

mutable struct Session_Pulse_Res <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Session_Pulse_Res(; kwargs...)
        obj = new(meta(Session_Pulse_Res), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Session_Pulse_Res
const __meta_Session_Pulse_Res = Ref{ProtoMeta}()
function meta(::Type{Session_Pulse_Res})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Session_Pulse_Res)
            __meta_Session_Pulse_Res[] = target = ProtoMeta(Session_Pulse_Res)
            allflds = Pair{Symbol,Union{Type,String}}[:alive => Bool]
            meta(target, Session_Pulse_Res, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Session_Pulse_Res[]
    end
end
function Base.getproperty(obj::Session_Pulse_Res, name::Symbol)
    if name === :alive
        return (obj.__protobuf_jl_internal_values[name])::Bool
    else
        getfield(obj, name)
    end
end

mutable struct Session_Pulse <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Session_Pulse(; kwargs...)
        obj = new(meta(Session_Pulse), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Session_Pulse
const __meta_Session_Pulse = Ref{ProtoMeta}()
function meta(::Type{Session_Pulse})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Session_Pulse)
            __meta_Session_Pulse[] = target = ProtoMeta(Session_Pulse)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Session_Pulse, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Session_Pulse[]
    end
end

mutable struct Session <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Session(; kwargs...)
        obj = new(meta(Session), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Session
const __meta_Session = Ref{ProtoMeta}()
function meta(::Type{Session})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Session)
            __meta_Session[] = target = ProtoMeta(Session)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Session, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Session[]
    end
end

export Session_Type, Session_Open_Req, Session_Open_Res, Session_Open, Session_Close_Req, Session_Close_Res, Session_Close, Session_Pulse_Req, Session_Pulse_Res, Session_Pulse, Session


text:
module meanfield

export ProductState, densityoperator

import ..integrate

using QuantumOpticsBase, LinearAlgebra
using ..interaction, ..CollectiveSpins

# Define Spin 1/2 operators
const spinbasis = SpinBasis(1//2)
const id = dense(identityoperator(spinbasis))
const sigmax_ = dense(sigmax(spinbasis))
const sigmay_ = dense(sigmay(spinbasis))
const sigmaz_ = dense(sigmaz(spinbasis))
const sigmap_ = dense(sigmap(spinbasis))
const sigmam_ = dense(sigmam(spinbasis))

"""
Class describing a Meanfield state (Product state).

The data layout is [sx1 sx2 ... sy1 sy2 ... sz1 sz2 ...]

# Arguments
* `N`: Number of spins.
* `data`: Vector of length 3*N.
"""
mutable struct ProductState{T1<:Int,T2<:Real}
    N::T1
    data::Vector{T2}
end

"""
    meanfield.ProductState(N)

Meanfield state with all Pauli expectation values equal to zero.
"""
ProductState(N::T) where T<:Int = ProductState(N, zeros(float(T), 3*N))

"""
    meanfield.ProductState(data)

Meanfield state created from real valued vector of length 3*spinnumber.
"""
ProductState(data::Vector{<:Real}) = ProductState(dim(data), data)

"""
    meafield.ProductState(rho)

Meanfield state from density operator.
"""
function ProductState(rho::DenseOpType)
    N = quantum.dim(rho)
    basis = quantum.basis(N)
    state = ProductState(N)
    sx, sy, sz = splitstate(s)
    f(ind, op) = real(expect(embed(basis, ind, op), rho))
    for k=1:N
        sx[k] = f(k, sigmax_)
        sy[k] = f(k, sigmay_)
        sz[k] = f(k, sigmaz_)
    end
    return state
end

"""
    meanfield.blochstate(phi, theta[, N=1])

Product state of `N` single spin Bloch states.

All spins have the same azimuthal angle `phi` and polar angle `theta`.
"""
function blochstate(phi::Vector{T1}, theta::Vector{T2}) where {T1<:Real, T2<:Real}
    N = length(phi)
    @assert length(theta)==N
    state = ProductState(N)
    sx, sy, sz = splitstate(state)
    for k=1:N
        sx[k] = cos(phi[k])*sin(theta[k])
        sy[k] = sin(phi[k])*sin(theta[k])
        sz[k] = cos(theta[k])
    end
    return state
end

function blochstate(phi::Real, theta::Real, N::Int=1)
    state = ProductState(N)
    sx, sy, sz = splitstate(state)
    for k=1:N
        sx[k] = cos(phi)*sin(theta)
        sy[k] = sin(phi)*sin(theta)
        sz[k] = cos(theta)
    end
    return state
end

"""
    meanfield.dim(state)

Number of spins described by this state.
"""
function dim(state::Vector{T}) where T<:Real
    N, rem = divrem(length(state), 3)
    @assert rem==0
    return N
end

"""
    meanfield.splitstate(N, data)
    meanfield.splitstate(state)

Split state into sx, sy and sz parts.
"""
splitstate(N::Int, data::Vector{<:Real}) = view(data, 1:1*N), view(data, 1*N+1:2*N), view(data, 2*N+1:3*N)
splitstate(state::ProductState) = splitstate(state.N, state.data)

"""
    meanfield.densityoperator(sx, sy, sz)
    meanfield.densityoperator(state)

Create density operator from independent sigma expectation values.
"""
function densityoperator(sx::Real, sy::Real, sz::Real)
    return 0.5*(id + sx*sigmax_ + sy*sigmay_ + sz*sigmaz_)
end
function densityoperator(state::ProductState)
    sx, sy, sz = splitstate(state)
    rho = densityoperator(sx[1], sy[1], sz[1])
    for i=2:state.N
        rho = tensor(rho, densityoperator(sx[i], sy[i], sz[i]))
    end
    return rho
end

"""
    meanfield.sx(state)

Sigma x expectation values of state.
"""
sx(x::ProductState) = view(x.data, 1:x.N)

"""
    meanfield.sy(state)

Sigma y expectation values of state.
"""
sy(x::ProductState) = view(x.data, x.N+1:2*x.N)

"""
    meanfield.sz(state)

Sigma z expectation values of state.
"""
sz(x::ProductState) = view(x.data, 2*x.N+1:3*x.N)

"""
    meanfield.timeevolution(T, S::SpinCollection, state0[; fout])

Meanfield time evolution.

# Arguments
* `T`: Points of time for which output will be generated.
* `S`: [`SpinCollection`](@ref) describing the system.
* `state0`: Initial ProductState.
* `fout` (optional): Function with signature `fout(t, state)` that is called whenever output
    should be generated.
"""
function timeevolution(T, S::SpinCollection, state0::ProductState; fout=nothing, kwargs...)
    N = length(S.spins)
    @assert N==state0.N
    Ω = interaction.OmegaMatrix(S)
    Γ = interaction.GammaMatrix(S)

    function f(dy, y, p, t)
        sx, sy, sz = splitstate(N, y)
        dsx, dsy, dsz = splitstate(N, dy)
        @inbounds for k=1:N
            dsx[k] = -S.spins[k].delta*sy[k] - 0.5*S.gammas[k]*sx[k]
            dsy[k] = S.spins[k].delta*sx[k] - 0.5*S.gammas[k]*sy[k]
            dsz[k] = -S.gammas[k]*(1+sz[k])
            for j=1:N
                if j==k
                    continue
                end
                dsx[k] += Ω[k,j]*sy[j]*sz[k] + 0.5*Γ[k,j]*sx[j]*sz[k]
                dsy[k] += -Ω[k,j]*sx[j]*sz[k] + 0.5*Γ[k,j]*sy[j]*sz[k]
                dsz[k] += Ω[k,j]*(sx[j]*sy[k] - sy[j]*sx[k]) - 0.5*Γ[k,j]*(sx[j]*sx[k] + sy[j]*sy[k])
            end
        end
    end

    if isa(fout, Nothing)
        fout_(t, state::ProductState) = deepcopy(state)
    else
        fout_ = fout
    end

    return integrate(T, f, state0, fout_; kwargs...)
end

"""
    meanfield.timeevolution_symmetric(T, state0, Ωeff, Γeff[; γ, δ0, fout])

Symmetric meanfield time evolution.

# Arguments
* `T`: Points of time for which output will be generated.
* `state0`: Initial ProductState.
* `Ωeff`: Effective dipole-dipole interaction.
* `Γeff`: Effective collective decay rate.
* `γ=1`: Single spin decay rate.
* `δ0=0`: Phase shift for rotated symmetric meanfield time evolution.
* `fout` (optional): Function with signature `fout(t, state)` that is called whenever output
    should be generated.
"""
function timeevolution_symmetric(T, state0::ProductState, Ωeff::Real, Γeff::Real; γ::Real=1.0, δ0::Real=0., fout=nothing, kwargs...)
    N = 1
    @assert state0.N==N
    function f(dy, y, p, t)
        sx, sy, sz = splitstate(N, y)
        dsx, dsy, dsz = splitstate(N, dy)
        dsx[1] = -δ0*sy[1] + Ωeff*sy[1]*sz[1] - 0.5*γ*sx[1] + 0.5*Γeff*sx[1]*sz[1]
        dsy[1] = δ0*sx[1] - Ωeff*sx[1]*sz[1] - 0.5*γ*sy[1] + 0.5*Γeff*sy[1]*sz[1]
        dsz[1] = -γ*(1+sz[1]) - 0.5*Γeff*(sx[1]^2+sy[1]^2)
    end

    if isa(fout, Nothing)
        fout_ = (t, state) -> deepcopy(state)
    else
        fout_ = fout
    end

    return integrate(T, f, state0, fout_; kwargs...)

end

"""
    meanfield.rotate(axis, angles, state)

Rotations on the Bloch sphere for the given [`ProductState`](@ref).

# Arguments
* `axis`: Rotation axis.
* `angles`: Rotation angle(s).
* `state`: [`ProductState`](@ref) that should be rotated.
"""
function rotate(axis::Vector{T1}, angles::Vector{T2}, state::ProductState) where {T1<:Real, T2<:Real}
    @assert length(axis)==3
    @assert length(angles)==state.N
    w = axis/norm(axis)
    sx, sy, sz = splitstate(state)
    state_rot = ProductState(state.N)
    sx_rot, sy_rot, sz_rot = splitstate(state_rot)
    v = zeros(T1, 3)
    for i=1:state.N
        v[1], v[2], v[3] = sx[i], sy[i], sz[i]
        θ = angles[i]
        sx_rot[i], sy_rot[i], sz_rot[i] = cos(θ)*v + sin(θ)*(w × v) + (1-cos(θ))*(w ⋅ v)*w
    end
    return state_rot
end

rotate(axis::Vector{T}, angle::Real, state::ProductState) where {T<:Real} = rotate(axis, ones(T, state.N)*angle, state)

end # module


text:
"EOS values from the literature to verify our implementations against"
module function_testvalues

const vinetpars = Dict(
    :ρ₀ => 8300.,
    :K₀ => 156.2e9,
    :dK₀ => 6.08,
)

const vinetdata = [
    # P           rho
    1.22307e+10   8848.49;
    4.96384e+10   9826.75;
    2.48563e+11   12772.1;
    1.33498e+12   19427.9;
    6.92316e+12   31970.1;
    2.9099e+13    54006.7;
    1.50906e+14   118578;
    6.12452e+14   260351;
    2.01457e+15   602402;
    5.76045e+15   1.32264e+06;
    1.43184e+16   2.98115e+06;
    3.55902e+16   6.89779e+06;
    1.05393e+17   2.24411e+07;
    2.27727e+17   5.76648e+07;
    3.34746e+17   9.4892e+07;
]

const bme3pars = Dict(
    :ρ₀ => 1460.,
    :K₀ => 23.7e9,
    :dK₀ => 4.15,
)

const bme3data = [
    # P           rho
    1e+09         1515.37;
    1.78526e+09   1554.68;
    3.34048e+09   1629.43;
    6.15341e+09   1752.08;
    1.1335e+10    1908.25;
    2.05554e+10   2141.38;
    3.61263e+10   2454.84;
    6.15341e+10   2802.2;
    1.03182e+11   3267.71;
    1.81344e+11   3876.19;
    3.04084e+11   4558.88;
    4.94171e+11   5407.78;
    8.96151e+11   6553.16;
    1.78526e+12   8358.54;
    3.90694e+12   11078.9;
    7.085e+12     13773.8;
    9.69158e+12   15390.7;
]

const bme4pars = Dict(
    :ρ₀ => 4100.,
    :K₀ => 247e9,
    :dK₀ => 3.97,
    :d2K₀ => -0.016e-9
)

const bme4data = [
    # P           rho
    1e+07         4271.43;
    2.6777e+08    4133.51;
    4.66182e+09   4120.99;
    2.82114e+10   4491.47;
    1.91993e+11   6010.25;
    1.35877e+12   10169.3;
    6.25222e+12   17726.5;
    2.66029e+13   34749.6;
    # unfortunately the BME doesn't appear to work properly beyond this point
]

const cgdata = [
    # P      V         T
    208.902  1.05615   251.15;
    212.166  0.889231  251.15;
    349.258  0.868462  256.43;
    352.522  0.818462  256.43;
    616.914  0.801538  272.23;
    620.178  0.748462  272.73;
]

const mgdpars = Dict(
    :V₀ => 12.3e-6,   # m3 /mol
    :K₀ => 23.9e9,    # Pa
    :dK₀ => 4.2,
    :α₀ => 150e-5,
    :T₀ => 300.,      # K
    :γ₀ => 1.2,
    :θD₀ => 1470.,    # K
    :q => -2,
    :n => 3,
)

const mgddata = Dict(
    300 => [
        # P       V
        0.89486	  11.9024;
        3.79145	  10.8593;
        6.44658	  10.2259;
        9.78022	  9.60449;
        12.6811	  9.19058;
        15.7994	  8.82874;
        18.9452	  8.52464;
    ],
    500 => [
        1.68223	  11.9606;
        4.06415	  11.0557;
        7.72163	  10.169;
        10.2963	  9.67982;
        13.7123	  9.19118;
        16.5594	  8.85804;
        19.5425	  8.57116;
    ],
    700 => [
        3.58171	  11.9443;
        6.61391	  10.8898;
        9.48582	  10.2104;
        12.4943	  9.65801;
        15.368	  9.23831;
        18.9201	  8.82477;
    ]
)

end # module function_testvalues


text:
el(
  todos = [
    include(joinpath("partials", "item.json.jl"))(t) for t = @vars(:todos)
  ]
)

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

__precompile__(true)

module ValidatedNumerics

using CRlibm
using Compat
using FixedSizeArrays
using ForwardDiff

import Base:
    +, -, *, /, //, fma,
    <, >, ==, !=, ⊆, ^, <=,
    in, zero, one, abs, real, min, max,
    sqrt, exp, log, sin, cos, tan, inv,
    exp2, exp10, log2, log10,
    asin, acos, atan, atan2,
    sinh, cosh, tanh, asinh, acosh, atanh,
    union, intersect, isempty,
    convert, promote_rule, eltype,
    BigFloat, float, widen, big,
    ∩, ⊆, eps,
    floor, ceil, trunc, sign, round,
    expm1, log1p,
    precision,
    isfinite, isnan,
    show, showall,
    isinteger, setdiff

export
    Interval, AbstractInterval,
    @interval, @biginterval, @floatinterval, @make_interval,
    diam, radius, mid, mag, mig, hull,
    emptyinterval, ∅, ∞, isempty, interior, isdisjoint, ⪽,
    precedes, strictprecedes, ≺,
    entireinterval, isentire, nai, isnai, isthin, iscommon,
    widen, infimum, supremum,
    parameters, eps, dist, roughly,
    pi_interval,
    midpoint_radius, interval_from_midpoint_radius,
    RoundTiesToEven, RoundTiesToAway,
    cancelminus, cancelplus, isunbounded,
    .., @I_str, ±

export
    displaymode

export RootFinding

if VERSION >= v"0.5.0-dev+1182"
    import Base: rounding, setrounding, setprecision
else
    import Compat:
        rounding, setrounding, setprecision

    export rounding, setrounding, setprecision  # reexport
end

## Multidimensional
export
    IntervalBox, @intervalbox

## Decorations
export
    @decorated,
    interval_part, decoration, DecoratedInterval,
    com, dac, def, trv, ill

function __init__()
    setrounding(BigFloat, RoundNearest)
    setrounding(Float64, RoundNearest)

    setprecision(Interval, 256)  # set up pi
    setprecision(Interval, Float64)
end

## Includes

include("intervals/intervals.jl")
include("multidim/multidim.jl")
include("decorations/decorations.jl")

include("display.jl")

include("root_finding/root_finding.jl")

end # module ValidatedNumerics


text:
using SingularityAppBuilder
using Test

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


text:
using WeakRefStrings, Missings
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end
@testset "WeakRefString{UInt8}" begin
    data = Vector{UInt8}("hey there sailor")

    str = WeakRefStrings.WeakRefString(pointer(data), 3)

    @test str.len == 3
    @test string(str) == "hey"
    @test String(str) == "hey"

    io = IOBuffer()
    show(io, str)
    @test String(take!(io)) == "\"hey\""

    show(io, typeof(str))
    @test String(take!(io)) == "WeakRefString{UInt8}"
end

@testset "WeakRefString{UInt16}" begin
    # simulate UTF16 data
    data = [0x0068, 0x0065, 0x0079]

    str = WeakRefStrings.WeakRefString(pointer(data), 3)
    @test typeof(str) == WeakRefStrings.WeakRefString{UInt16}
    @test String(str) == "hey"
    @test str.len == 3
end

@testset "WeakRefString{UInt32}" begin
    # simulate UTF32 data with trailing null byte
    data = [0x00000068, 0x00000065, 0x00000079, 0x00000000]

    str = WeakRefStrings.WeakRefString(pointer(data), 4)
    @test typeof(str) == WeakRefStrings.WeakRefString{UInt32}
    @test String(str) == "hey"
    @test str.len == 4
end

@testset "WeakRefArray" begin
    data = Vector{UInt8}("hey there sailor")
    str = WeakRefStrings.WeakRefString(pointer(data), 3)
    C = WeakRefStringArray(data, [str])
    @test length(C) == 1
    @test eltype(C) === String

    A = WeakRefStringArray(UInt8[], WeakRefStrings.WeakRefString{UInt8}, 10)
    @test size(A) == (10,)
    @test A[1] == ""
    @test A[1, 1] == ""
    A[1] = "hey"
    A[1, 1] = "hey"
    @test A[1] == "hey"
    resize!(A, 5)
    @test size(A) == (5,)
    push!(A, str)
    @test size(A) == (6,)
    @test A[end] == str
    push!(A, "hi")
    @test length(A) == 7
    @test A[end] == convert(WeakRefStrings.WeakRefString{UInt8}, "hi")

    B = WeakRefStringArray(UInt8[], WeakRefStrings.WeakRefString{UInt8}, 10)
    B[1] = "ho"
    append!(A, B)
    @test size(A) == (17,)

    D = WeakRefStringArray(UInt8[], Union{Missing, WeakRefStrings.WeakRefString{UInt8}}, 0)
    push!(D, "hey")
    push!(D, str)
    push!(D, missing)
    @test length(D) == 3
    @test eltype(D) == Union{Missing, String}
    @test D[2] == str
    @test D[3] === missing
    D[2] = missing
    @test D[2] === missing

    E = WeakRefStringArray(data, [str missing])
    @test size(E) == (1, 2)
    @test eltype(E) === Union{String, Missing}

end


text:
using SpikesCorrelations
using Test

@testset "SpikesCorrelations.jl" begin
    # Write your tests here.
    include("test_CorrelationIndex.jl")
    include("test_SpikeTrain.jl")
    include("test_TilingCoefficient.jl")
    include("test_VPDistance.jl")
end


text:
using AllenCahn
using Test
using LinearAlgebra

@testset "Neumann" begin
    @test include("neumanntest1.jl")
    @test include("neumanntest2.jl")
end

@testset "Periodic" begin
    @test include("periodictest1.jl")
    @test include("periodictest2.jl")
end


text:
# From HyperDualNumbers
const DIFF_RULES = [
    (  :sqrt,    :(1/2/sqrt(x)),                       :(-1/4/(sqrt(x) * x)))
    (  :cbrt,    :(1/3/x^(2/3)),                       :(-2/9/x^(5/3)))
    (  :abs2,    :(2*x),                               :(2))
    (  :abs,     :(signbit(x) ? -one(x) : one(x)),     :(0))
    (  :inv,     :(-1/x^2),                            :(2/x^3))
    (  :log,     :(1/x),                               :(-1/x^2))
    (  :log10,   :(1/(log(10)*x)),                     :(-1/(log(10)*x^2)))
    (  :log2,    :(1/(log(2)*x)),                      :(-1/(log(2)*x^2)))
    (  :log1p,   :(1/(x + 1)),                         :(-1/(x + 1)^2))
    (  :exp,     :(exp(x)),                            :(exp(x)))
    (  :exp2,    :(2^x*log(2)),                        :(2^x*log(2)^2))
    (  :exp10,   :(10^x*log(10)),                      :(10^x*log(10)^2))
    (  :expm1,   :(exp(x)),                            :(exp(x)))
   #( :sin,     :(cos(x)),                             :(-sin(x)))
   #( :cos,     :(-sin(x)),                            :(-cos(x)))
    (  :tan,     :(tan(x)^2 + 1),                      :(2*(tan(x)^2 + 1)*tan(x)))
    (  :sec,     :(sec(x)*tan(x)),                     :(sec(x)*tan(x)^2 + (tan(x)^2 + 1)*sec(x)))
    (  :csc,     :(-cot(x)*csc(x)),                    :(cot(x)^2*csc(x) + (cot(x)^2 + 1)*csc(x)))
    (  :cot,     :(-cot(x)^2 - 1),                     :(2*(cot(x)^2 + 1)*cot(x)))
    (  :sind,    :(π*cos(π*x/180)/180),                :(-π^2*sin(π*x/180)/180^2))
    (  :cosd,    :(-π*sin(π*x/180)/180),               :(-π^2*cos(π*x/180)/180^2))
    (  :tand,    :(π*(tan(π*x/180)^2 + 1)/180),        :(2*π^2*(tan(π*x/180)^2 + 1)*tan(π*x/180)/180^2))
    (  :secd,    :(π*sec(π*x/180)*tan(π*x/180)/180),   :(π^2*sec(π*x/180)*tan(π*x/180)^2/180^2 + π^2*(tan(π*x/180)^2 + 1)*sec(π*x/180)/180^2))
    (  :cscd,    :(-π*cot(π*x/180)*csc(π*x/180)/180),  :(π^2*cot(π*x/180)^2*csc(π*x/180)/180^2 + π^2*(cot(π*x/180)^2 + 1)*csc(π*x/180)/180^2))
    (  :cotd,    :(-π*(cot(π*x/180)^2 + 1)/180),       :(2*π^2*(cot(π*x/180)^2 + 1)*cot(π*x/180)/180^2))
    (  :asin,    :(1/sqrt(-x^2 + 1)),                  :(x/(-x^2 + 1)^(3/2)))
    (  :acos,    :(-1/sqrt(-x^2 + 1)),                 :(-x/(-x^2 + 1)^(3/2)))
    (  :atan,    :(1/(x^2 + 1)),                       :(-2*x/(x^2 + 1)^2))
    (  :asec,    :(1/(sqrt(x^2 - 1)*x)),               :(-1/(x^2 - 1)^(3/2) - 1/(sqrt(x^2 - 1)*x^2)))
    (  :acsc,    :(-1/(sqrt(x^2 - 1)*x)),              :(1/(x^2 - 1)^(3/2) + 1/(sqrt(x^2 - 1)*x^2)))
    (  :acot,    :(-1/(x^2 + 1)),                      :(2*x/(x^2 + 1)^2))
    (  :asind,   :(π/(sqrt(-π^2*x^2/180^2 + 1)*180)),  :(π^3*x/((-π^2*x^2/180^2 + 1)^(3/2)*180^3)))
    (  :acosd,   :(-π/(sqrt(-π^2*x^2/180^2 + 1)*180)), :(-π^3*x/((-π^2*x^2/180^2 + 1)^(3/2)*180^3)))
    (  :atand,   :(π/((π^2*x^2/180^2 + 1)*180)),       :(-2*π^3*x/((π^2*x^2/180^2 + 1)^2*180^3)))
    (  :asecd,   :(1/(sqrt(π^2*x^2/180^2 - 1)*x)),     :(-π^2/((π^2*x^2/180^2 - 1)^(3/2)*180^2) - 1/(sqrt(π^2*x^2/180^2 - 1)*x^2)))
    (  :acscd,   :(-1/(sqrt(π^2*x^2/180^2 - 1)*x)),    :(π^2/((π^2*x^2/180^2 - 1)^(3/2)*180^2) + 1/(sqrt(π^2*x^2/180^2 - 1)*x^2)))
    (  :acotd,   :(-π/((π^2*x^2/180^2 + 1)*180)),      :(2*π^3*x/((π^2*x^2/180^2 + 1)^2*180^3)))
    (  :sinh,    :(cosh(x)),                           :(sinh(x)))
    (  :cosh,    :(sinh(x)),                           :(cosh(x)))
    (  :tanh,    :(-tanh(x)^2 + 1),                    :(2*(tanh(x)^2 - 1)*tanh(x)))
    (  :sech,    :(-sech(x)*tanh(x)),                  :(sech(x)*tanh(x)^2 + (tanh(x)^2 - 1)*sech(x)))
    (  :csch,    :(-coth(x)*csch(x)),                  :(coth(x)^2*csch(x) + csch(x)/sinh(x)^2))
    (  :coth,    :(-1/sinh(x)^2),                      :(2*cosh(x)/sinh(x)^3))
    (  :asinh,   :(1/sqrt(x^2 + 1)),                   :(-x/(x^2 + 1)^(3/2)))
    (  :acosh,   :(1/(sqrt(x + 1)*sqrt(x - 1))),       :(-1/2/(sqrt(x + 1)*(x - 1)^(3/2)) - 1/2/((x + 1)^(3/2)*sqrt(x - 1))))
    (  :atanh,   :(-1/(x^2 - 1)),                      :(2*x/(x^2 - 1)^2))
    (  :asech,   :(-1/(sqrt(-x^2 + 1)*x)),             :(-1/(-x^2 + 1)^(3/2) + 1/(sqrt(-x^2 + 1)*x^2)))
    (  :acsch,   :(-1/(sqrt(x^2 + 1)*x)),              :(1/(x^2 + 1)^(3/2) + 1/(sqrt(x^2 + 1)*x^2)))
    (  :acoth,   :(-1/(x^2 - 1)),                      :(2*x/(x^2 - 1)^2))
    (  :deg2rad, :(π/180),                             :(0))
    (  :rad2deg, :(180/π),                             :(0))
  # ( :erf,     :(2*exp(-x^2)/sqrt(π)),               :(-4*x*exp(-x^2)/sqrt(π)))
  # ( :erfinv,  :(1/2*sqrt(π)*exp(erfinv(x)^2)),      :(1/2*π*erfinv(x)*exp(2*erfinv(x)^2)))
  # ( :erfc,    :(-2*exp(-x^2)/sqrt(π)),              :(4*x*exp(-x^2)/sqrt(π)))
  # ( :erfi,    :(2*exp(x^2)/sqrt(π)),                :(4*x*exp(x^2)/sqrt(π)))
    (  :cospi,   :(-π*sinpi(x)),                       :(-π^2*cospi(x)))
    (  :sinpi,   :(π*cospi(x)),                        :(-π^2*sinpi(x)))
]

changeprecision(x) = x
changeprecision(x::Integer) = :(T($x))
function changeprecision(x::Symbol)
    x == :π  ?  :(T($x)) : x
end
function changeprecision(ex::Expr)
    if Meta.isexpr(ex, :call, 3) && ex.args[1] == :^ && ex.args[3] isa Int
        return Expr(:call, :^, changeprecision(ex.args[2]), ex.args[3])
    else
        return Expr(ex.head, changeprecision.(ex.args)...)
    end
end

for (i, rule) in enumerate(DIFF_RULES)
    prec_f′ = changeprecision(rule[2])
    prec_f′′ = changeprecision(rule[3])
    DIFF_RULES[i] = (rule[1], prec_f′, prec_f′′)
end


text:
#=
  Tug of war
  From Hakaru example documentation/tugofwar_rejection.hk

  This variant include laziness

  cf ~/webppl/tug_of_war.wppl
     ~/blog/tug_of_war.blog
=#
using Turing, StatsPlots, DataFrames

@model function tug_of_war(n)

    # strength = Vector{Real}(undef, n)
    # lazy = Vector{Bool}(undef, n)
    # pulls = Vector{Real}(undef, n)
    strength = TArray{Real}(undef, n)
    lazy = TArray{Real}(undef, n)
    pulls = TArray{Real}(undef, n)

    for p in 1:n
        strength[p] ~ truncated(Normal(0,1),0,10) # latent
        lazy[p]     ~ Bernoulli(0.1)
        # Note: This works, but since the pulls[p] is a derived
        #       value, then the pulls data are not in the chns.
        #       So we have to restore this in post process.
        if lazy[p]
            pulls[p] = strength[p] / 2.0
        else
            pulls[p] = strength[p]
        end

    end

    function winner(a,b)
        return pulls[a] > pulls[b] ? a : b
    end

    for p1 in 1:n-1, p2 in p1+1:n
        true ~ Dirac(winner(p1,p2) == p1)
    end
    return pulls
end

num_people = 3
model = tug_of_war(num_people)

num_chns = 4
num_samples = 10_000 # 10_000

# 2.1s Good rhat
chns = sample(model,  MH(), MCMCThreads(), num_samples, num_chns)

# PG(20) doesn't differentiate the people's pull/strength!
# Very slow: 24.6s
# chns = sample(model,  PG(20), MCMCThreads(), num_samples, num_chns)

# IS() doesn't differentiate the people's pull/strength
# chns = sample(model,  IS(), MCMCThreads(), num_samples, num_chns)

# Error
# chns = sample(model,  HMC(0.1, 5), MCMCThreads(), num_samples, num_chns)

# chns = sample(model,  HMCDA(200, 0.65, 0.3), MCMCThreads(), num_samples, num_chns)
# chns = sample(model,  NUTS(1000, 0.65), MCMCThreads(), num_samples, num_chns)

# Too slow
# chns = sample(model,  SMC(), MCMCThreads(), num_samples, num_chns)

# chns = sample(model,  NUTS(1000,0.65), MCMCThreads(), num_samples, num_chns) # Error
# Not correct: lazy is all 0 and neither strength nor pulls are "ordered"
# (and lots of rejected proposals)
# chns = sample(model,  Gibbs(NUTS(1000,0.65,:pulls,:strength),MH(:lazy)), MCMCThreads(), num_samples, num_chns)

# chns = sample(model,  HMC(0.1, 5), MCMCThreads(), num_samples, num_chns) # error
# This works: 6.6s. Good rhats.
# chns = sample(model,  Gibbs(HMC(0.1,5,:pulls,:strength),MH(:lazy)), MCMCThreads(), num_samples, num_chns)

display(chns)
# display(gelmandiag(chns))
# display(plot(chns))
df = DataFrame(chns)

# Show the latent strength
println("strength (latent):")
for p1 in 1:num_people-1, p2 in p1+1:num_people
    println("strength[$p1] > strength[$p2]: ", mean(df[!,"strength[$p1]"] .>df[!,"strength[$p2]"]))
end

if "pulls[1]" in names(df)
    # This is for the now commented version of pulls.
    # But that don't work well.
    println("pulls (observed):")
    for p1 in 1:num_people-1, p2 in p1+1:num_people
        println("pulls[$p1] > pulls[$p2]: ", mean(df[!,"pulls[$p1]"] .>df[!,"pulls[$p2]"]))
    end
else
   println("Restoring pulls")
   #
   # Since the pulls are not available from the chns we have to restore it
   # (Note: This is one of the disadvantages that Turing has compared to e.g. WebPPL!)
   println("\npulls")
   pulls = []
   for row in eachrow(df)
       r = vcat([ row["lazy[$p]"] == 1.0 ? row["strength[$p]"] / 2 : row["strength[$p]"]  for p in 1:num_people]...)
       push!(pulls, r)
   end

   println("\npulls (restored):")
   println("pulls overall mean: ", mean(pulls))

   println("\ncomparison of pulls")
   for p1 in 1:num_people-1, p2 in p1+1:num_people
       # println("pulls[$p1] > pulls[$p2]: ", mean(pulls[:,p1] .> pulls[:,p2]))
       println("pulls[$p1] > pulls[$p2]: ", mean([p[p1] > p[p2] for p in pulls]))
   end

end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
struct SubTabular{N, I<:Tuple, D} <: AbstractTabular{N}
    indices::I
    data::D
end

const SubSeries{I <: Tuple, D} = SubTabular{1,I,D}
const SubTable{I <: Tuple, D} = SubTabular{2,I,D}

SubTabular{N}(indices::I, data::D) where {N,I,D} = SubTabular{N,I,D}(indices, data)

# function aresubindices(new_inds::NTuple{N,Any}, data) where {N}
#     old_inds = _indices(Dims(N), data)
#     shape = index_shape(Dims(N), new_inds, data)
#
#     return all(_map(issubindex, shape, old_inds, new_inds))
# end

@inline function indices(t::SubTabular{N}) where {N}
    full_inds = _indices(Dims(length(t.indices)), t.data)
    subshape = index_subshape(t)
    if subshape === (NonscalarIndex(),)
        return (replace_colon(t.indices[1], full_inds[1]),)
    elseif subshape === (NonscalarIndex(), ScalarIndex())
        return (replace_colon(t.indices[1], full_inds[1]),)
    elseif subshape === (ScalarIndex(), NonscalarIndex())
        return (replace_colon(t.indices[2], full_inds[2]),)
    elseif subshape === (NonscalarIndex(), NonscalarIndex())
        return (replace_colon(t.indices[1], full_inds[1]), replace_colon(t.indices[2], full_inds[2]))
    else
        error("Index subshape = $subshape")
    end
end

replace_colon(::Colon, fullindex) = fullindex
replace_colon(subindex, fullindex) = subindex

issubindex(::ScalarIndex, i, inds) = i ∈ inds
issubindex(::NonscalarIndex, i, inds) = i ⊆ inds
issubindex(::NonscalarIndex, ::Colon, inds) = true
issubindex(::NonscalarIndex, i, ::Colon) = true # TODO error("Unsure")
issubindex(::NonscalarIndex, ::Colon, ::Colon) = true

# inject_indices(::Tuple{ScalarIndex}, ::Tuple{}, old_inds) = old_inds
@inline function inject_indices(::Tuple{NonscalarIndex}, new_inds::Tuple{Any}, old_inds)
    @boundscheck if !issubindex(NonscalarIndex(), new_inds[1], old_inds[1])
        throw(IndexError("Indices $(new_inds[1]) not found in Tabular"))
    end
    if new_inds isa Tuple{Colon}
        return old_inds
    else
        return new_inds
    end
end

# inject_indices(::Tuple{ScalarIndex,ScalarIndex}, ::Tuple{}, old_inds) = old_inds
@inline function inject_indices(::Tuple{NonscalarIndex,ScalarIndex}, new_inds::Tuple{Any}, old_inds)
    @boundscheck if !issubindex(NonscalarIndex(), new_inds[1], old_inds[1])
        throw(IndexError("Indices $(new_inds[1]) not found in Tabular"))
    end
    if new_inds isa Tuple{Colon}
        return old_inds
    else
        return (new_inds[1], old_inds[2])
    end
end
@inline function inject_indices(::Tuple{ScalarIndex,NonscalarIndex}, new_inds::Tuple{Any}, old_inds)
    @boundscheck if !issubindex(NonscalarIndex(), new_inds[1], old_inds[2])
        throw(IndexError("Indices $(new_inds[1]) not found in Tabular"))
    end
    if new_inds isa Tuple{Colon}
        return old_inds
    else
        return (old_inds[1], new_inds[1])
    end
end
@inline function inject_indices(::Tuple{NonscalarIndex,NonscalarIndex}, new_inds::Tuple{Any,Any}, old_inds)
    @boundscheck if !issubindex(NonscalarIndex(), new_inds[1], old_inds[1])
        throw(IndexError("Indices $(new_inds[1]) not found in Tabular"))
    end
    @boundscheck if !issubindex(NonscalarIndex(), new_inds[2], old_inds[2])
        throw(IndexError("Indices $(new_inds[2]) not found in Tabular"))
    end
    if new_inds isa Tuple{Colon, Any}
        if new_inds isa Tuple{Colon, Colon}
            return old_inds
        else
            return (old_inds[1], new_inds[2])
        end
    else
        if new_inds isa Tuple{Any, Colon}
            return (new_inds[1], old_inds[2])
        else
            return new_inds
        end
    end
end

# Helper function to get shape of subindices with respect to the full structure
function index_subshape(s::SubSeries)
    if length(s.indices) === 1
        return _index_shape(Dims(1), typeof(s.data), typeof(s.indices[1]))
    else
        return _index_shape(Dims(2), typeof(s.data), typeof(s.indices[1]), typeof(s.indices[2]))
    end
end
function index_subshape(s::SubTable)
    _index_shape(Dims(2), typeof(s.data), typeof(s.indices[1]), typeof(s.indices[2]))
end

index_shape(s::SubSeries, ::Any) = (ScalarIndex(),)
index_shape(s::SubSeries, ::AbstractVector) = (NonscalarIndex(),)
index_shape(s::SubSeries, ::Colon) = (NonscalarIndex(),)

index_shape(s::SubTable, ::Any, ::Any) = (ScalarIndex(), ScalarIndex())
index_shape(s::SubTable, ::AbstractVector, ::Any) = (NonscalarIndex(), ScalarIndex())
index_shape(s::SubTable, ::Colon, ::Any) = (NonscalarIndex(), ScalarIndex())
index_shape(s::SubTable, ::Any, ::AbstractVector) = (ScalarIndex(), NonscalarIndex())
index_shape(s::SubTable, ::AbstractVector, ::AbstractVector) = (NonscalarIndex(), NonscalarIndex())
index_shape(s::SubTable, ::Colon, ::AbstractVector) = (NonscalarIndex(), NonscalarIndex())
index_shape(s::SubTable, ::Any, ::Colon) = (ScalarIndex(), NonscalarIndex())
index_shape(s::SubTable, ::AbstractVector, ::Colon) = (NonscalarIndex(), NonscalarIndex())
index_shape(s::SubTable, ::Colon, ::Colon) = (NonscalarIndex(), NonscalarIndex())

#@inline function index_shape(s::SubSeries, i1)
#    inds = indices[i1]

#    subshape = index_subshape(s)

#    fullinds = _indices(Dims(N), t.data)
#    subinds = inject_indices((i1,) fullinds)
#    _index_shape(Dims(length(fullinds)), typeof(t.data), typeof(subinds[1]))
#end

# function index_shape(t::SubSeries, i1, i2)
#     fullinds = _indices(Dims(N), t.data)
#     subinds = inject_indices((i1,i2), fullinds)
#     _index_shape(Dims(length(fullinds)), typeof(t.data), typeof(subinds[1]), typeof(subinds[2]))
# end

@propagate_inbounds function getindex(s::SubSeries, i)
    inds = inject_indices(index_subshape(s), (i,), s.indices)#, _indices(Dims(length(s.indices)), s.data))
    if length(inds) === 1
        shape = _index_shape(Dims(length(inds)), typeof(s.data), typeof(inds[1]))
    elseif length(inds) === 2
        shape = _index_shape(Dims(length(inds)), typeof(s.data), typeof(inds[1]), typeof(inds[2]))
    end
    data = _getindex(Dims(length(inds)), s.data, inds...)
    if shape isa Tuple{Vararg{ScalarIndex}}
        return data
    else
        return Series(data)
    end
end
@propagate_inbounds function getindex(t::SubTable, i1, i2)
    inds = inject_indices(index_subshape(t), (i1,i2), t.indices)#, _indices(Dims(length(t.indices)), t.data))
    shape = _index_shape(Dims(length(inds)), typeof(t.data), typeof(inds[1]), typeof(inds[2]))
    data = _getindex(Dims(length(inds)), t.data, inds[1], inds[2])
    if shape[1] isa ScalarIndex
        if shape[2] isa ScalarIndex
            return data
        else
            return Series(data)
        end
    else
        if shape[2] isa ScalarIndex
            return Series(data)
        else
            return Table(data)
        end
    end
end

@propagate_inbounds function setindex!(t::SubTabular{N}, val, i::Vararg{Any, N}) where N
    inds = inject_indices(index_subshape(t), i, t.indices)
    if length(inds) === 1
        shape = _index_shape(Dims(length(inds)), typeof(t.data), typeof(inds[1]))
    else
        shape = _index_shape(Dims(length(inds)), typeof(t.data), typeof(inds[1]), typeof(inds[2]))
    end

    if shape isa Tuple{Vararg{ScalarIndex}}
        _setindex!(Dims(length(inds)), t.data, val, inds...)
    else
        error("Non-scalar setindex! has not yet been implemented.")
    end
    return val
end

@inline function view(s::Series, i)
    shape = index_shape(s, i)

    if shape[1] == ScalarIndex()
        throw("Cannot take the view of a scalar index")
    end

    # Check if all indices in `i` are a part of the original indices
    @boundscheck begin
        if !(issubindex(NonscalarIndex(), i, indices(s)[1]))
            throw(IndexError("Indices $i are not in the input Series"))
        end
    end

    return SubSeries((i,), s.data)
end
@inline function view(s::SubSeries, i)
    subshape = index_subshape(s)
    shape = index_shape(s, i)

    if shape[1] == ScalarIndex()
        throw("Cannot take the view of a scalar index")
    end

    # Check if all indices in `i` were a part of the original indices
    @boundscheck begin
        if !(issubindex(NonscalarIndex(), i, indices(s)[1]))
            throw(IndexError("Indices $i are not in the input Series"))
        end
    end

    if subshape === (NonscalarIndex(),)
        if i isa Colon
            return SubSeries((indices(s)[1],), s.data)
        else
            return SubSeries((i,), s.data)
        end
    elseif subshape === (NonscalarIndex(), ScalarIndex())
        if i isa Colon
            return SubSeries((indices(s)[1], s.indices[2]), s.data)
        else
            return SubSeries((i, s.indices[2]), s.data)
        end
    elseif subshape === (ScalarIndex(), NonscalarIndex())
        if i isa Colon
            return SubSeries((s.indices[1], indices(s)[1]), s.data)
        else
            return SubSeries((s.indices[1], i), s.data)
        end
    end
end

@inline function view(t::Table, i1, i2)
    shape = index_shape(t, i1, i2)

    @boundscheck begin
        inds = indices(t)
        if !issubindex(shape[1], replace_colon(i1, inds[1]), inds[1])
            throw(IndexError("Indices $i1 are not in the input Table columns"))
        end
        if !issubindex(shape[2], i2, inds[2])
            throw(IndexError("Indices $i2 are not in the input Table rows"))
        end
    end

    if shape[1] == ScalarIndex()
        if shape[2] == ScalarIndex()
            throw("Cannot take the view of a scalar index")
        else
            SubSeries((i1, i2), t.data)
        end
    else
        if shape[2] == ScalarIndex()
            SubSeries((i1, i2), t.data)
        else
            SubTable((i1, i2), t.data)
        end
    end
end

@inline function view(t::SubTable, i1, i2)
    shape = index_shape(t, i1, i2)

    @boundscheck begin
        # Don't need the output, but this will check for inclusion
        inject_indices(index_subshape(t), (i1,i2), indices(t))
    end

    if shape[1] == ScalarIndex()
        if shape[2] == ScalarIndex()
            throw("Cannot take the view of a scalar index")
        else
            if i2 isa Colon
                SubSeries((i1, t.indices[2]), t.data)
            else
                SubSeries((i1, i2), t.data)
            end
        end
    else
        if shape[2] == ScalarIndex()
            if i1 isa Colon
                SubSeries((indices(t)[1], i2), t.data)
            else
                SubSeries((i1, i2), t.data)
            end
        else
            if i1 isa Colon
                if i2 isa Colon
                    SubTable(indices(t), t.data)
                else
                    SubTable((indices(t)[1], i2), t.data)
                end
            else
                if i2 isa Colon
                    SubTable((i1, indices(t)[2]), t.data)
                else
                    SubTable((i1, i2), t.data)
                end
            end
        end
    end
end


text:
int main() {
    int x;

    {
        int x = 10;
        while (x > 0) {
            printInt(x);
            x--;
        }
    }

    printInt(x);

    return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
using NonUniformRandomVariateGeneration
using Distributions
using BenchmarkTools

function benchNative(α::Float64, β::Float64, N::Int64)
  m::Float64 = 0.0
  for i in 1:N
    m += sampleGamma(α, β)
  end
  return m / N
end

function benchDistributions(α::Float64, β::Float64, N::Int64)
  gam = Distributions.Gamma(α, 1/β)
  m::Float64 = 0.0
  for i in 1:N
    m += rand(gam)
  end
  return m / N
end

@btime benchNative(0.1, 1.0, 1024*1024*8)
@btime benchDistributions(0.1, 1.0, 1024*1024*8)

@btime benchNative(15.0, 1.0, 1024*1024*8)
@btime benchDistributions(15.0, 1.0, 1024*1024*8)


text:
# This file is part of Kpax3. License is MIT.

include("mcmc.jl")
include("write.jl")


text:
@testset "Comparison" begin
    @testset "Graded Lex Order" begin
        Mod.@polyvar x y z
        @test x > y > z
        @test x^2*y > y^3 > z
        @test y^2 >= x
        # Examples from p. 58, 59 of the 4th edition of "Ideals, Varieties, and Algorithms" of Cox, Little and O'Shea
        @test x^1*y^2*z^3 > x^3*y^2
        @test !(x^1*y^2*z^3 < x^3*y^2)
        @test x^1*y^2*z^4 > x^1*y^1*z^5
        @test !(x^1*y^2*z^4 < x^1*y^1*z^5)
        @test x^4*y^7*z > x^4*y^2*z^3
        @test !(x^4*y^7*z < x^4*y^2*z^3)
        @test x*y^5*z^2 < x^4*y*z^3
        @test !(x*y^5*z^2 > x^4*y*z^3)
        @test x^5*y*z > x^4*y*z^2
        @test !(x^5*y*z < x^4*y*z^2)
    end
    @testset "Equality" begin
        @testset "Monomial equality" begin
            Mod.@polyvar x y
            @test 1 == constantmonomial(x)
            @test 2 != constantmonomial(x)
            @test 2 != x
            @test 2 != x*y
            @test x*y != x
            @test 1x*y == x*y
            @test 2x*y != x*y
            @test x == monomial(x)
            @test x*y^0 == x
            @test x != x^0*y
            @test x^2 != x*y
        end
        @testset "Polynomial equality" begin
            Mod.@polyvar x y
            @test isapproxzero(0.0x, ztol=0.0)
            @test polynomial(CustomTerms(x + 1 - x)) isa AbstractPolynomial
            @test MultivariatePolynomials.eqconstant(polynomial(CustomTerms(x + 1 - x)), 1)
            @test MultivariatePolynomials.eqconstant(CustomTerms(x + 1 - x), 1)
            @test 2 * CustomPoly(x + 1 - x) == 2
            @test 2 != CustomTerms(x + 1 - x)
            @test 3x^2 == CustomTerms(x - x + x^2) * 3
            @test CustomPoly(-x + x^2) != x^2
            @test 2*x*y + 3*y^2 == 3*y^2 + 2*y*x
            @test 3*x*y + 2*y^2 != 3*y^2 + 2*y*x
            @test x + y != x * (1 + y)
            @test x*y == 3x + 2x*y - x - x*y - 2x
            @test isapproxzero((1+1e-8)x - x, ztol=1e-7)
            @test !isapproxzero((1+1e-6)x - x, ztol=1e-7)
            @test isapprox((2-1e-3)*x*y + (3+1e-3)*y^2, 3*y^2 + 2*y*x, rtol=1e-2)
            @test !isapprox((2-1e-3)*x*y + (3+1e-1)*y^2, 3*y^2 + 2*y*x, rtol=1e-2)
            @test isapprox(1e-3*x*y + 3*y^2 + x^2, x^2 + 3*y^2, rtol=1e-2, ztol=1e-2)
            @test isapprox(1e-3*x*y + 3*y^2 + x^2, x^2 + 3*y^2, rtol=1e-2, atol=1e-2)
            @test isapprox(3*y^2 + x^2, x^2 + 1e-3*x*y + 3*y^2, rtol=1e-2, ztol=1e-2)
            @test isapprox(3*y^2 + x^2, x^2 + 1e-3*x*y + 3*y^2, rtol=1e-2, atol=1e-2)
            @test !isapprox(3*y^2 + x^2, x^2 + 1e-1*x*y + 3*y^2, rtol=1e-2, ztol=1e-2)
            @test !isapprox(3.0*y^2 + x + x^2, x + 3*y^2, rtol=1e-2, ztol=1e-2)
            @test isapprox(x+1-x, 1+1e-8, rtol=1e-7)
            @test !isapprox(1+1e-8, x+1-x, rtol=1e-9)
        end
        @testset "RationalPoly equality" begin
            Mod.@polyvar x y
            @test (x^2 - x - 6) / (x + 2) != x + 3
            @test x - 3 == (x^2 - x - 6) / (x + 2)
            @test (x^2 - x - 6) / (x - 3) == x + 2
            @test 3 != 4x / 2x
            @test 4x / 2x == 2
            @test 3 != 4x / 2x
            @test 1 - 1/x == (x - 1) / x
            @test 1 - 1/x != (1 - x) / x
            @test x + x/x == 1 + x^2/x
            @test x - x/x == -(1 - x^2/x)
            @test (1+x)/x - 1 == 1/x
            @test isapprox((1+1e-8)x, (x*y)/y, rtol=1e-7)
            @test isapproxzero(((1+1e-8)x - x)/y, ztol=1e-7)
            @test !isapproxzero(((1+1e-8)x - y)/y, ztol=1e-9)
            @test isapprox(((1+1e-8)x*y) / y^2, x / y, rtol=1e-7)
            @test isapprox((2x) / x, 2.001, rtol=1e-2)
            @test !isapprox(2.001,  (2x) / x, rtol=1e-4)
        end
    end
    @testset "Equality between a Polynomial and a type not defining zero #22" begin
        Mod.@polyvar x
        # Polynomial of multiple terms
        p = x + x^2
        @test nothing != p
        @test p != nothing
        @test Dict{Int,Int}() != p
        @test p != Dict{Int,Int}()
        # Polynomial of one term
        p = x + x^2 - x
        @test p != nothing
        @test p != Dict{Int,Int}()
        # Polynomial of no term
        # See https://github.com/blegat/MultivariatePolynomials.jl/issues/22
        p = x - x
        @test p != nothing
        @test p != Dict{Int,Int}()
    end
end


text:
#
# Provides out-of-the-box Two-body systems in our solar system.
#


text:
using Test
using DataStructures

include("../optimizers/cma_es.jl")
include("optimizers/cma_es_deap.jl")
include("optimizers/compare_optimizer_states.jl")

number_generations = 100
population_size = 200
sigma = 1.5
free_parameters = 1000

tolerance = 0.00001

@testset "Optimizers" begin

    optimizer_configuration = OrderedDict("sigma" => sigma, "population_size" => population_size)

    # Initialize Optimizers
    # OptimizerCmaEsDeap: Original Deap CMA-ES optimizer implemented in Python using PyCall
    # OptimizerCmaEs: CMA-ES optimizer implemented in Julia that we use in our framework
    optimizer1, eigenvectors1, indx1 = OptimizerCmaEsDeap(free_parameters, optimizer_configuration)
    optimizer2, eigenvectors2, indx2, diagD = OptimizerCmaEs(free_parameters, optimizer_configuration, test = true, eigenvectors1 = eigenvectors1, indx1 = indx1)

    # Compare return values
    @test size(indx1) == size(indx2)
    @test diagD[indx1] ≈ diagD[indx2] atol = 0.00001
    @test size(eigenvectors1) == size(eigenvectors2)

    # Compare internal states of both optimizers
    compare_optimizer_states(optimizer1, optimizer2, tolerance)

    for generation = 1:number_generations

        # Ask optimizers for new population
        genomes1, randoms1 = ask(optimizer1)
        genomes2, randoms2 = ask(optimizer2, test = true, randoms1 = randoms1)

        # Compare return values
        @test genomes1 ≈ genomes2 atol = tolerance
        @test size(randoms1) == size(randoms2)
        @test mean(randoms1) ≈ mean(randoms2) atol = 0.1
        @test std(randoms1) ≈ std(randoms2) atol = 0.01

        # Generate random rewards
        rewards_training = rand(population_size)

        # Tell optimizers new rewards
        eigenvectors1, indx1 = tell(optimizer1, rewards_training)
        eigenvectors2, indx2, diagD = tell(optimizer2, rewards_training, test = true, eigenvectors1 = eigenvectors1, indx1 = indx1)

        # Compare return values
        @test size(indx1) == size(indx2)
        @test diagD[indx1] ≈ diagD[indx2] atol = 0.00001
        @test size(eigenvectors1) == size(eigenvectors2)

        # Compare internal states of both optimizers
        compare_optimizer_states(optimizer1, optimizer2, tolerance)

        # Test if C is a Hermitian matrix
        @test optimizer2.C ≈ Hermitian(optimizer2.C) atol = tolerance

    end
end


text:
using FunctionalCollections
using Test

@testset "Persistent Lists" begin

    @testset "length" begin
        @test length(PersistentList([1])) == 1
        @test length(PersistentList([1, 2, 3])) == 3
        @test length(EmptyList()) == 0
    end

    @testset "equality" begin
        @test PersistentList(1:100) == PersistentList(1:100)
        @test PersistentList(1:100) != PersistentList(1:99)
        @test PersistentList(1:100) == collect(1:100)
    end

    @testset "head" begin
        @test head(PersistentList(1:100)) == 1
        @test head(PersistentList([1]))   == 1
        @test try head(EmptyList()); false catch e true end
    end

    @testset "tail" begin
        @test tail(PersistentList(1:100)) == PersistentList(2:100)
        @test tail(PersistentList([1]))     == EmptyList()
        @test try tail(EmptyList()); false catch e true end
    end

    @testset "cons" begin
        @test cons(1, cons(2, cons(3, EmptyList()))) == PersistentList([1, 2, 3])
        @test 1..(2..(3..EmptyList())) == PersistentList([1, 2, 3])
    end

    @testset "sharing" begin
        l = PersistentList(1:100)
        l2 = 0..l
        @test l === tail(l2)
    end

    @testset "iteration" begin
        arr2 = Int[]
        for i in PersistentList(1:1000)
            push!(arr2, i)
        end
        @test collect(1:1000) == arr2
    end

    @testset "map" begin
        @test map(x->x+1, PersistentList([1,2,3,4,5])) == PersistentList([2,3,4,5,6])
    end

    @testset "filter" begin
        @test filter(x->x%2==0, PersistentList([1,2,3,4,5,6])) == PersistentList([2,4,6])
    end

    @testset "reverse" begin
        @test reverse(PersistentList(1:10)) == 10:-1:1
    end

    @testset "hash" begin
        @test hash(PersistentList(1:1000)) == hash(PersistentList(1:1000))
    end

    @testset "isempty" begin
        @test isempty(EmptyList())
        @test !isempty(PersistentList([1]))
    end

end


text:
function demoSumProdStandardiseVariables()
    #DEMOSUMPROD Sum-Product algorithm test :
    # Variable order is arbitary
    r=randperm(10)
    variables=r[1:5] # In this code, the variables do not need to be numbered 1:end
    a, b, c, d, e=variables
    nstates=round(Int64,3*rand(1,5)+2) # random number of states for each variable

    pot=Array(PotArray,5)
    pot[1]=PotArray([a b],rand(nstates[memberinds([a b],variables)]...))
    pot[2]=PotArray([b c d],rand(nstates[memberinds([b c d],variables)]...));
    pot[3]=PotArray([c],rand(nstates[memberinds([c],variables)]...));
    pot[4]=PotArray([e d],rand(nstates[memberinds([e d],variables)]...));
    pot[5]=PotArray([d],rand(nstates[memberinds([d],variables)]...));

    newpot,originalvariables=standardisevariables(pot) # translate the variables into 1:end form
    A = FactorGraph(newpot);
    newmarg, mess, normconst=sumprod(newpot,A)

    marg=returnvariables(newmarg,originalvariables) # translate back to orginal variables
    # check if this is correct:
    jointpot =prod(pot); V=length(pot);
    for i=variables
        mtable=[condpot(marg[whichpot(marg,i)]).content condpot(jointpot,i).content]
        println("marginal of variable $i:\n Factor Graph\t\t Raw summation");
        println(mtable)
    end
    tmp=sum(prod(pot))
    println("Normalisation constant:\nFactorGraph\t=$(normconst.content)\nRaw summation\t=$(tmp.content)")

end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the -(x, y) function
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- File inclusions

# Test utilities
include("test_utils.jl")

# --- Tests

# ------ M::Multivector

@testset "-(M::Multivector, N::Multivector)" begin
    @test_skip 1
end

@testset "-(M::Multivector, N::Blade)" begin
    @test_skip 1
end

@testset "-(M::Multivector, N::Pseudoscalar)" begin
    @test_skip 1
end

@testset "-(M::Multivector, N::Scalar)" begin
    @test_skip 1
end

@testset "-(M::Multivector, N::One)" begin
    @test_skip 1
end

@testset "-(M::Multivector, N::Zero)" begin
    @test_skip 1
end

@testset "-(M::Multivector, N::Real)" begin
    @test_skip 1
end

@testset "-(M::Multivector, N::Vector)" begin
    @test_skip 1
end

# ------ B::Blade

@testset "-(B::Blade, M::Multivector)" begin
    @test_skip 1
end

@testset "-(B::Blade, C::Blade)" begin
    @test_skip 1
end

@testset "-(B::Blade, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "-(B::Blade, C::Scalar)" begin
    @test_skip 1
end

@testset "-(B::Blade, C::One)" begin
    @test_skip 1
end

@testset "-(B::Blade, C::Zero)" begin
    test_dim = 10
    B = Blade(rand(test_dim, 3))
    C = Zero()
    @test B - C === B
end

@testset "-(B::Blade, C::Real)" begin
    @test_skip 1
end

@testset "-(B::Blade, C::Vector)" begin
    @test_skip 1
end

# ------ B::Pseudoscalar

@testset "-(B::Pseudoscalar, M::Multivector)" begin
    @test_skip 1
end

@testset "-(B::Pseudoscalar, C::Blade)" begin
    @test_skip 1
end

@testset "-(B::Pseudoscalar, C::Pseudoscalar)" begin
    test_dim = 10

    # B != C
    test_value_1 = get_random_value(1)  # add 1 to keep value away from 0
    B = Pseudoscalar(test_dim, test_value_1)

    test_value_2 = abs(test_value_1) + rand() + 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Pseudoscalar(test_dim, test_value_2)

    B_minus_C = B - C
    @test B_minus_C isa Pseudoscalar
    @test B_minus_C.dim == test_dim
    @test B_minus_C.value == test_value_1 - test_value_2

    # B == C
    B = Pseudoscalar(test_dim, test_value_1)
    C = Pseudoscalar(test_dim, test_value_1)
    @test iszero(B - C)
end

@testset "-(B::Pseudoscalar, C::Scalar)" begin
    @test_skip 1
end

@testset "-(B::Pseudoscalar, C::One)" begin
    @test_skip 1
end

@testset "-(B::Pseudoscalar, C::Zero)" begin
    test_dim = 10
    test_value = get_random_value(1)  # add 1 to keep value away from 0

    B = Pseudoscalar(test_dim, test_value)
    C = Zero()
    @test B - C === B
end

@testset "-(B::Pseudoscalar, C::Real)" begin
    @test_skip 1
end

@testset "-(B::Pseudoscalar, C::Vector)" begin
    @test_skip 1
end

# ------ B::Scalar

@testset "-(B::Scalar, M::Multivector)" begin
    @test_skip 1
end

@testset "-(B::Scalar, C::Blade)" begin
    @test_skip 1
end

@testset "-(B::Scalar, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "-(B::Scalar, C::Scalar)" begin
    # --- Preparations

    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value_1)

    # --- Tests

    # B != C, C != B - 1
    test_value_2 = abs(test_value_1) + rand() + 2
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Scalar(test_value_2)

    B_minus_C = B - C
    @test B_minus_C isa Scalar
    @test B_minus_C.value == test_value_1 - test_value_2

    # B == C
    C = Scalar(test_value_1)
    @test iszero(B - C)

    # C == B + 1
    C = Scalar(test_value_1 - 1)

    B_minus_C = B - C
    @test B_minus_C isa One
end

@testset "-(B::Scalar, C::One)" begin
    # --- Preparations

    C = One()

    # --- Tests

    # B != 2
    test_value = get_random_value(3)  # add 3 to keep value away from 0, 1, and 2
    B = Scalar(test_value)

    B_minus_C = B - C
    @test B_minus_C isa Scalar
    @test B_minus_C.value == test_value - 1

    # B == 2
    B = Scalar(2)

    B_minus_C = B - C
    @test B_minus_C isa One
end

@testset "-(B::Scalar, C::Zero)" begin
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value)

    C = Zero()

    @test B - C === B
end

@testset "-(B::Scalar, C::Real)" begin
    # --- Preparations

    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value_1)

    # --- Tests

    # B != C, C != B - 1
    test_value_2 = abs(test_value_1) + rand() + 2
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = test_value_2

    B_minus_C = B - C
    @test B_minus_C isa Scalar
    @test B_minus_C.value == test_value_1 - test_value_2

    # B == C
    C = test_value_1
    @test iszero(B - C)

    # C == B - 1
    C= test_value_1 - 1

    B_minus_C = B - C
    @test B_minus_C isa One
end

@testset "-(B::Scalar, C::Vector)" begin
    @test_skip 1
end

# ------ B::One

@testset "-(B::One, M::Multivector)" begin
    @test_skip 1
end

@testset "-(B::One, C::Blade)" begin
    @test_skip 1
end

@testset "-(B::One, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "-(B::One, C::Scalar)" begin
    B = One()

    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value)

    B_minus_C = B - C
    @test B_minus_C isa Scalar
    @test B_minus_C.value == 1 - test_value
end

@testset "-(B::One, C::One)" begin
    B = One()
    C = One()
    @test iszero(B - C)
end

@testset "-(B::One, C::Zero)" begin
    B = One()
    C = Zero()

    B_minus_C = B - C
    @test B_minus_C === B
end

@testset "-(B::One, C::Real)" begin
    # B != C, C != B - 1
    B = One()

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = test_value

    B_minus_C = B - C
    @test B_minus_C isa Scalar
    @test B_minus_C.value == 1 - test_value

    # B == C
    B = One()
    C = 1
    @test iszero(B - C)

    # C == B - 1
    C = 0

    B_minus_C = B - C
    @test B_minus_C isa One
end

@testset "-(B::One, C::Vector)" begin
    @test_skip 1
end

# ------ B::Zero

@testset "-(B::Zero, M::Multivector)" begin
    @test_skip 1
end

@testset "-(B::Zero, C::Blade)" begin
    test_dim = 10
    B = Zero()
    C = Blade(rand(test_dim, 3))
    @test B - C == -C
end

@testset "-(B::Zero, C::Pseudoscalar)" begin
    test_dim = 10
    test_value = get_random_value(1)  # add 1 to keep value away from 0

    B = Zero()
    C = Pseudoscalar(test_dim, test_value)
    @test B - C == -C
end

@testset "-(B::Zero, C::Scalar)" begin
    # --- Preparations

    B = Zero()

    # --- Tests

    # C != -1
    test_value = get_random_value(2)  # add 2 to keep value away from 0, 1, and -1
    C = Scalar(test_value)

    B_minus_C = B - C
    @test B_minus_C isa Scalar
    @test B_minus_C.value == -test_value

    # C == -1
    C = Scalar(-1)

    B_minus_C = B - C
    @test B_minus_C isa One
end

@testset "-(B::Zero, C::One)" begin
    B = Zero()
    C = One()
    B_minus_C = B - C
    @test B_minus_C isa Scalar
    @test B_minus_C.value == -1
end

@testset "-(B::Zero, C::Zero)" begin
    B = Zero()
    C = Zero()

    B_minus_C = B - C
    @test B_minus_C === B
end

@testset "-(B::Zero, C::Real)" begin
    # --- Preparations

    B = Zero()

    # --- Tests

    # B != C, C != -1
    test_value = get_random_value(2)  # add 2 to keep value away from 0, 1, and -1
    C = test_value

    B_minus_C = B - C
    @test B_minus_C isa Scalar
    @test B_minus_C.value == -test_value

    # B == C
    B = Zero()
    C = 0
    @test iszero(B - C)

    # C == -1
    C = -1

    B_minus_C = B - C
    @test B_minus_C isa One
end

@testset "-(B::Zero, C::Vector)" begin
    B = Zero()
    C = rand(5)
    @test B - C == -C
end

# ------ B::Real

@testset "-(B::Real, M::Multivector)" begin
    @test_skip 1
end

@testset "-(B::Real, C::Blade)" begin
    @test_skip 1
end

@testset "-(B::Real, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "-(B::Real, C::Scalar)" begin
    # --- Preparation

    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = test_value_1

    # --- Tests

    # B != C, C != B - 1
    test_value_2 = abs(test_value_1) + rand() + 2
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Scalar(test_value_2)

    B_minus_C = B - C
    @test B_minus_C isa Scalar
    @test B_minus_C.value == test_value_1 - test_value_2

    # B == C
    C = Scalar(test_value_1)
    @test iszero(B - C)

    # C == B - 1
    C = Scalar(test_value_1 - 1)

    B_minus_C = B - C
    @test B_minus_C isa One
end

@testset "-(B::Real, C::One)" begin
    # --- Preparations

    C = One()

    # --- Tests

    # B != C, B != C + 1
    test_value = get_random_value(3)  # add 2 to keep value away from 0, 1, and 2
    B = test_value

    B_minus_C = B - C
    @test B_minus_C isa Scalar
    @test B_minus_C.value == test_value - 1

    # B == C
    B = 1
    @test iszero(B - C)

    # B = C + 1
    B = 2

    B_minus_C = B - C
    @test B_minus_C isa One
end

@testset "-(B::Real, C::Zero)" begin
    # --- Preparations

    C = Zero()

    # --- Tests

    # B != C, B != C + 1
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = test_value

    B_minus_C = B - C
    @test B_minus_C isa Scalar
    @test B_minus_C.value == test_value

    # B == C
    B = 0
    C = Zero()
    @test iszero(B - C)

    # B == C + 1
    B = 1

    B_minus_C = B - C
    @test B_minus_C isa One
end

# ------ B::Vector

@testset "-(B::Vector, M::Multivector)" begin
    @test_skip 1
end

@testset "-(B::Vector, C::Blade)" begin
    @test_skip 1
end

@testset "-(B::Vector, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "-(B::Vector, C::Scalar)" begin
    @test_skip 1
end

@testset "-(B::Vector, C::One)" begin
    @test_skip 1
end

@testset "-(B::Vector, C::Zero)" begin
    B = rand(5)
    C = Zero()
    @test B - C === B
end


text:
# return (ipython, version) tuple, where ipython is the string of the
# IPython executable, and version is the VersionNumber.
function find_ipython()
    try
        "ipython",convert(VersionNumber, chomp(readall(`ipython --version`)))
    catch e1
        try
            "ipython2",convert(VersionNumber, chomp(readall(`ipython2 --version`)))
        catch e2
            try
                "ipython3",convert(VersionNumber, chomp(readall(`ipython3 --version`)))
            catch e3
                try
                    "ipython.bat",convert(VersionNumber, chomp(readall(`ipython.bat --version`)))
                catch e4
                    error("IPython is required for IJulia, got errors\n",
                          "   $e1\n   $e2\n   $e3" * (@windows ? "\n$e4\n" : "") )
                end
            end
        end
    end
end


text:
"""
    pyiter(x)

Equivalent to `iter(x)` in Python.
"""
pyiter(x) = pynew(errcheck(@autopy x C.PyObject_GetIter(getptr(x_))))
export pyiter

pynext(x::Py) = pynew(errcheck_ambig(C.PyIter_Next(getptr(x))))
export pyiter


text:
"""

```
coor(atoms; selection, options)
```

Returns the coordinates of the atoms. The input may be one atom (type `Atom`), a vector of atoms, or a `Residue`. For vectors or residues, selection strings or functions may be applied. If `xyz_in_cols=true`, the 3D coordinates will be the three columns of the matrix, thus likely the highest dimension will be the first, what is usually optimal for speed since Julia is column-major.

### Examples

```julia-repl
julia> protein = wget("1LBD");

julia> coor(protein[1])
(45.228, 84.358, 70.638)

julia> coor(protein,"index <= 2")
2×3 Matrix{Float64}:
 45.228  84.358  70.638
 46.08   83.165  70.327

julia> coor(protein,"residue = 1",xyz_in_cols=false)
3×6 Matrix{Float64}:
 45.228  46.08   45.257  45.823  47.147  46.541
 84.358  83.165  81.872  80.796  82.98   82.639
 70.638  70.327  70.236  69.974  71.413  72.662

julia> coor(protein, only = at -> at.resname == "ALA", xyz_in_cols=false)
3×110 Matrix{Float64}:
 43.94   43.02   41.996  41.228  42.293  28.854  27.821  26.415  …  17.282  -16.085  -16.377  -17.866  -18.496  -15.888
 81.982  80.825  80.878  79.937  80.676  82.285  81.232  81.715     81.631   84.599   84.019   84.088   83.942   82.583
 70.474  70.455  69.34   69.157  71.812  60.096  60.093  60.408     80.842   53.001   51.738   51.741   52.777   51.706

julia> residues = collect(eachresidue(protein));

julia> coor(residues[1])
6×3 Matrix{Float64}:
 45.228  84.358  70.638
 46.08   83.165  70.327
 45.257  81.872  70.236
 45.823  80.796  69.974
 47.147  82.98   71.413
 46.541  82.639  72.662

```

"""
coor(atom::Atom) = atom.x, atom.y, atom.z

function coor(atoms::AbstractVector{Atom}, selection::String; xyz_in_cols::Bool=true)
  query = parse_query(selection)
  return coor(atoms,only = atom -> apply_query(query,atom), xyz_in_cols=xyz_in_cols)
end

function coor(atoms::AbstractVector{Atom}; only=all, xyz_in_cols::Bool=true)
  n = 0
  for atom in atoms
    if only(atom)
      n = n + 1
    end
  end
  if xyz_in_cols
    x = Matrix{Float64}(undef,n,3)
  else
    x = Matrix{Float64}(undef,3,n)
  end
  i = 0
  for atom in atoms
    if only(atom)
      i = i + 1
      if xyz_in_cols
        x[i,1:3] .= coor(atom)
      else
        x[1:3,i] .= coor(atom)
      end
    end
  end
  return x
end

#
# Coordinates of the atoms of a residue/molecule
#
coor(residue::Residue; only=all, xyz_in_cols::Bool=true) =
  coor(residue.atoms[residue.range], only=only, xyz_in_cols=xyz_in_cols)
coor(residue::Residue, selection::String; xyz_in_cols::Bool=true) =
  coor(residue.atoms[residue.range], selection, xyz_in_cols=xyz_in_cols)


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
function visit_depth_first(a::ANode, previsit::Function, postvisit::Function, revisit::Function, visited::Set{ANode}=Set{ANode}())
  if sym(a) == :constant return end
  if a ∉ visited
    push!(visited, a)
    previsit(a)
    for b in a.input
      visit_depth_first(b, previsit, postvisit, revisit, visited)
    end
    postvisit(a)
  else
    revisit(a)
  end
end

function find_subgraph(f::ANode, x::ANode...)
  subgraph = Set{ANode}(x)
  path = ANode[]

  add_path() = for i in path push!(subgraph, i) end
  previsit(n::ANode) = push!(path, n)
  revisit(n::ANode) = if n ∈ subgraph add_path() end
  function postvisit(n::ANode)
    @assert pop!(path) == n
    if n ∈ subgraph add_path() end
  end
  visit_depth_first(f, previsit, postvisit, revisit)
  subgraph
end

immutable GradientCache
  objective::ANode
  gradients::Dict{ANode, ANode}
end

gradients(f::ANode) = GradientCache(f, Dict{ANode, ANode}(f=>ones(ANode, size(f))))

constant(a::ANode) = ANode(:constant, (a,), size(a))
do_forward!(a::ANode{:constant}, y, x) = copy!(y, x)

function Base.getindex(c::GradientCache, x::ANode)
  subgraph = find_subgraph(c.objective, x)
  outputs = Dict{ANode, Vector{ANode}}()

  for node in subgraph
    for input in node.input
      if input ∈ subgraph
        push!(get!(outputs, input, ANode[]), node)
      end
    end
  end
  @assert x ∈ subgraph "Invalid gradient.  There is no path between the variable being differentiated and the independent variable.  Technically the gradient is zero, but this usually indicates that you're computing a gradient w.r.t. the wrong node."

  visited = Set{ANode}([c.objective])
  ready = ANode[c.objective]

  while !haskey(c.gradients, x)
    if isempty(ready)
     return zeros(ANode, size(x))
    end
    node = pop!(ready)

    if node ∉ subgraph continue end
    if !haskey(c.gradients, node)
      c.gradients[node] = reduce(+, map(i->gradient_node(i, node, c.gradients[i]), outputs[node]))
    end

    for a in node.input
      if a in visited continue end
      if a ∉ subgraph continue end
      if sym(a) == :constant continue end

      isready = true
      for b in outputs[a]
        if !haskey(c.gradients, b)
          isready = false
          break
        end
      end

      if isready
        push!(ready, a)
        push!(visited, a)
      end
    end
  end
  c.gradients[x]
end


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Foundations of machine learning")

# # Tools - Using arrays to store data

#-

# ### Introduction to arrays
#
# **Arrays are collections of boxes that we can use to store data.** In the last notebook, we saw an image that can help us to picture a 1-dimensional array:
#
#
#
# *Why do we want an object like this to store our data?*

#-

# An alternative to using an array in some contexts would be to name every individual piece of data, as follows:
#
# ```julia
# a = 1.1
# b = 2.2
# c = 3.3
# ```
#
# We can visualize how this data is stored:
#
#
#
#
#
# This is like having a separate box for every piece of data, rather than a series of connected boxes for all our data.
#
# The more data we have, the more annoying it becomes to keep track of all these boxes and their names. Furthermore, if we want to do the same thing with many pieces of data, it's much easier to put all of these pieces of data in one place to work with them at once.
#
# For example, we may want to multiply `a`, `b`, and `c` by `2`. We could multiply three times:
#
# ```julia
# a * 2
# b * 2
# c * 2
# ```
#
# Or, instead, we could create one array (let's call it `numbers`) and multiply that array by `2`:
# ```julia
# numbers * 2
# ```
#
# The syntax for creating this array, `numbers`, is
#
# ```julia
# numbers = [a, b, c]
# ```
#
# Or, we could have just written
#
# ```julia
# numbers = [1.1, 2.2, 3.3]
# ```
#
# It's worth noting that in Julia, 1-dimensional arrays are also called "vectors".

#-

# ### Creating arrays
#
# In the last section, we saw that we could create the array `numbers` by typing our elements, `a`, `b`, and `c` (or `1.1`, `2.2`, and `3.3`), inside square brackets.

#-

# #### Exercise 1
#
#
#
# Create an array called `first_array` that stores the numbers 10 through 20.
#
# Reminder: ESC+b to open a box below.

#-

# #### Solution

first_array = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]

#-

[√i for i=1:10] # \sqrt + <tab> gives √

#-

π

# ### Array comprehensions

#-

# Alternatively we can create arrays via *"array comprehensions"*. This is a nice way to automate creating an array, if you don't want to have to type long lists of numbers inside square brackets.
#
# In an array comprehension, you write code inside square brackets that will generate the array as that code is executed. Let's see an example.
#
# ```julia
# counting = [2 * i for i in 1:10]
# ```
#
# The above line of code is saying that we want to create an array called `counting` that stores two times each of the integers between 1 and 10. This array comprehension has a few different parts, so let's dissect them:
#
#
#
# If we wanted to create another array, `roots`, that stored the square roots of all integers between 1 and 10, we could modify the code above to
#
# ```julia
# roots = [sqrt(i) for i in 1:10]
# ```
#
# #### Exercise
#
# Create an array, `squares`, that stores the squares of all integers between 1 and 100.

#-

# #### Solution

squares = [x^2 for x in 1:100]

# ### Looking inside an array
#
# We can "index" into an array to grab contents inside the array at a particular position. To index into our `counting` array, we place square brackets after the name of the array and put the position number of the element/data we want inside those brackets.
#
# For example, we might say
#
# ```julia
# counting[3]
# ```
#
# to grab the third element in the array called `counting`.

#-

# #### Exercise 2
#
# Execute the code in the next cell. It will generate an array called `myprimes` that stores the first 168 prime numbers. Index into `myprimes` to grab the 89th smallest prime. What is this prime?

using Primes

myprimes = primes(1000); # The semicolon suppresses the output, try removing it

# #### Solution
#
# Via `myprimes[89]`, we see the 89th smallest prime is `461`:

myprimes[89]

# ### Slicing

#-

# Instead of grabbing a single number in an array, we can also take a **slice** of an array, i.e. a subset of the array, which can include multiple values. To take a slice of `counting` that includes the 3rd, 4th, and 5th entries, we use the following syntax with a colon (`:`):
#
# ```julia
# counting[3:5]
# ```

#-

# #### Exercise 3
#
# Index into `myprimes` to grab the 89th through the 99th smallest primes (inclusive).

#-

# #### Solution

myprimes[89:99]

# ### Modifying an array
#
# In this section, we'll see how to edit an array we've already created. We'll see how to
#
# 1) Update an item at an existing position in an array;
#
# 2) Add an item to the end of an array;
#
# 3) Remove an item from the end of an array.

#-

# #### Update an item at an existing position in an array
#
# First off, we can update an item at an existing position in an array by indexing into the array and changing the value at the desired index. For example, let's say we have an array called `myfriends`:
#
# ```julia
# myfriends = ["Ted", "Robin", "Barney", "Lily", "Marshall"]
# ```
#
# We can grab Barney's name from `myfriends` via
#
# ```julia
# myfriends[3]
# ```
#
# and we can change "Barney" to something else by reassigning `myfriends[3]`:
#
# ```julia
# myfriends[3] = "Baby Bop"
# ```
#
# Note that a single `=` assigns a new variable to the value on the left of the `=` sign.

#-

# #### Exercise 4
#
# Use an array comprehension to create an array, `mysequence`, that stores the numbers 4 through 10. Index into `mysequence` and update it to replace the last element, `10`, with `11`.

#-

# #### Solution
#
# Create `mysequence` via

mysequence = [x for x in 4:10]

# then update the last element using the `end` keyword

mysequence[end] = 11

# or using the length of `mysequence`, such as

mysequence[length(mysequence)] = 11

# #### Add an item to the end of an array
#
# We can add an item to the end of an array with the `push!` function. Don't worry about the exclamation mark at the end of `push!`; we'll talk about why that's there later, when we discuss functions.
#
# For now, note that when you call `push!` on an array and a value, you add that value to the end of the array. For example,
#
# ```julia
# push!(counting, 1000)
# ```
#
# will turn `counting`, which was formerly
#
# ```julia
# [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
# ```
# into
# ```julia
# [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1000]
# ```
#
# Prove to yourself this works!

#-

# #### Exercise 5
#
# Copy the following code to declare the array `fibonacci`:
#
# ```julia
# fibonacci = [1, 1, 2, 3, 5, 8, 13]
# ```
#
# Use `push!` to add `21` to `fibonacci` after the number `13`.

#-

# #### Solution

fibonacci = [1, 1, 2, 3, 5, 8, 13]
push!(fibonacci, 21)

# #### Remove an item from the end of an array
#
# To remove an item from the end of an array, use the `pop!` function. When using `pop!` on arrays, `pop!` only needs one input argument, namely the array you want to change. For example, if
#
# ```julia
# counting == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1000]
# ```
#
# then
# ```julia
# pop(counting)
# ```
# will remove `1000` from the end of `counting`.

#-

# #### Exercise 6
#
# Use `pop!` to remove `21` from `fibonacci`. What does this function call return?

#-

# **Solution**
#
# We can remove `21` via the following code, which also **returns** `21`:

pop!(fibonacci)

# #### Exercise 7
#
# What is the last element of `fibonacci`? Try
#
# ```julia
# fibonacci[end]
# ```

#-

# **Solution**
#
# The last element is **`13`**. *Hint*: If you're getting a different answer, you may have run the exercises in this notebook out of order, or called `pop!` multiple times. Start with the array
#
# ```julia
# fibonacci = [1, 1, 2, 3, 5, 8, 13, 21]
# ```
#
# and call
#
# ```julia
# pop!(fibonacci)
# ```
#
# *Now* what is the last element of `fibonacci`?

#-

# ### Arrays of arrays and multidimensional arrays
#
# So far we've only seen 1D arrays of scalars, but arrays can also store other arrays and can have an arbitrary number of dimensions.
#

# For example, the following are arrays of arrays:

favorites = [ ["koobideh", "chocolate", "eggs"], ["penguins", "cats", "sugargliders"] ]

#-

numbers = [ [1, 2, 3], [4, 5], [6, 7, 8, 9] ]

# One way to create a multidimensional array is to use the `rand` function. `rand` takes an arbitrary number of arguments where the number of arguments determines the number of dimensions the array created by `rand` will have!
#
# For example, if we pass `rand` two integer inputs, it will generate a 2D array, while three integer inputs will generate a 3D array:

rand(4, 3)

#-

rand(4, 3, 2)

# If we want to grab a value from a 2D array, we index into the array, specifying the row and column of the element of interest! For example, if we have the array `A` given by
#
# ```julia
# A = [ 1  2  3
#       4  5  6
#       7  8  9 ]
#  ```
#  we could grab the number 6 by saying
#
#  ```
#  A[2, 3]
#  ```
#
#  since 6 is in the 2nd row and 3rd column of `A`.

#-

# #### Exercise 8
#
#  Copy and execute the following code to get an array, `myprimematrix`, which stores the first 100 primes
#
#  ```
#  myprimematrix = reshape(primes(541), (10, 10))
#  ```
#
#  Grab the prime in the 8th row and 5th column via indexing. What is the value of this prime?

#-

# #### Solution
#
# Grabbing the prime in the 8th row and 5th column of `myprimematrix` yields the prime `223`:

myprimematrix = reshape(primes(541), (10, 10))
myprimematrix[8, 5]

# ### Copying arrays
#
# Be careful when you want to copy arrays!
#
# Let's say you want to copy an existing array, like `fibonacci`, to an array called `somenumbers`. Remember how we defined `fibonacci`:
#
# ```julia
# fibonacci = [1, 1, 2, 3, 5, 8, 13]
# ```
#
# What if we try to copy fibonacci by saying
#
# ```julia
# somenumbers = fibonacci
# ```
# ?
# Execute this code to try to copy fibonacci. Look at `somenumbers` to see if it stores what you want.

fibonacci

#-

somenumbers = fibonacci

#-

somenumbers

# Now let's update `somenumbers` by changing its first element.

somenumbers[1] = 404

# Now let's look inside fibonacci.

fibonacci

# #### Exercise 9
#
# What is the first element in `fibonacci`?

#-

# #### Solution
#
# After running

fibonacci = [1, 1, 2, 3, 5, 8, 13]
somenumbers = fibonacci
somenumbers[1] = 404

# we can check the value of the first element of `fibonacci` via

fibonacci[1]

# to see that the first element of `fibonacci` is `404`.

#-

# ### Copying or not?

#-

# Did copying `fibonacci` like this work?
#
# No, unfortunately not. When we tried to copy, all we did was give `fibonacci` a new name, `somenumbers`. Now when we update `somenumbers`, we're also updating `fibonacci` because they are the same object!

#-

# If we actually want to make a *copy* of the array bound to `fibonacci`, we can use the `copy` function:

## First, restore fibonacci

fibonacci[1] = 1
fibonacci

#-

somemorenumbers = copy(fibonacci)

#-

somemorenumbers[1] = 404

#-

fibonacci

# In this last example, fibonacci was not updated. Therefore we see that the arrays bound to `somemorenumbers` and `fibonacci` are distinct.

#-

# #### Exercise 10
#
# Copy `myprimematrix` to `mynewprimematrix`. Update `mynewprimematrix[3,3]` to `1234`.

#-

# #### Solution

mynewprimematrix = copy(myprimematrix)
mynewprimematrix[3, 3] = 1234

# **Tests**

#-

# ```julia
# mynewprimematrix[3, 3] == 1234
# myprimematrix[3, 3] == 83
# ```

max_stars_repo_path: Courses/Foundations of machine learning/0100.Representing-data-in-a-computer-quiz.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials