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

struct IndexMap
    n::Int
    m₁::Int
    m₂::Int
    map::Dict{RowCol, IdxGroup}

    function IndexMap(n::Integer, m₁::Integer, m₂::Integer)
        return new(n, m₁, m₂, Dict{RowCol, IdxGroup}())
    end
end
Base.setindex!(map::IndexMap, val::IdxGroup, key::RowCol) = setindex!(map.map, val, key)
Base.getindex(map::IndexMap, key::RowCol) = getindex(map.map, key)

struct LPData{T <: AbstractFloat, M <: AbstractMatrix}
    c₁::Vector{T}
    c₂::T
    A::M
    b::Vector{T}
    d₁::Vector{T}
    C::M
    d₂::Vector{T}
    lb::Vector{T}
    ub::Vector{T}
    is_binary::Vector{Bool}
    is_integer::Vector{Bool}
    indexmap::IndexMap

    function LPData(::Type{T}, ::Type{<:Matrix}, n::Int, m₁::Int, m₂::Int) where T <: AbstractFloat
        return new{T, Matrix{T}}(
            zeros(T, n),
            zero(T),
            zeros(T, m₁, n),
            zeros(T, m₁),
            fill(convert(T, -Inf), m₂),
            zeros(T, m₂, n),
            fill(convert(T, Inf), m₂),
            fill(convert(T, -Inf), n),
            fill(convert(T, Inf), n),
            fill(false, n),
            fill(false, n),
            IndexMap(n, m₁, m₂))
    end

    function LPData(::Type{T}, ::Type{<:SparseMatrixCSC}, n::Int, m₁::Int, m₂::Int) where T <: AbstractFloat
        return new{T, SparseMatrixCSC{T,Int}}(
            zeros(T, n),
            zero(T),
            spzeros(T, m₁, n),
            zeros(T, m₁),
            fill(convert(T, -Inf), m₂),
            spzeros(T, m₂, n),
            fill(convert(T, Inf), m₂),
            fill(convert(T, -Inf), n),
            fill(convert(T, Inf), n),
            fill(false, n),
            fill(false, n),
            IndexMap(n, m₁, m₂))
    end

    function LPData(::Type{T}, c₁::AbstractVector, c₂::Real,
                    A::AbstractMatrix, b::AbstractVector,
                    d₁::AbstractVector, C::AbstractMatrix, d₂::AbstractVector,
                    lb::AbstractVector, ub::AbstractVector,
                    is_binary::AbstractVector, is_integer::AbstractVector,
                    map::IndexMap) where T <: AbstractFloat
        A_ = convert(Matrix{T}, A)
        C_ = convert(Matrix{T}, C)
        # Sanity check
        _lpcheck(c₁, c₂, A, b, d₁, C, d₂, lb, ub, is_binary, is_integer)
        # Return data
        return new{T, Matrix{T}}(c₁, c₂, A_, b, d₁, C_, d₂, lb, ub, is_binary, is_integer, map)
    end

    function LPData(::Type{T}, c₁::AbstractVector, c₂::Real,
                    A::AbstractSparseMatrix, b::AbstractVector,
                    d₁::AbstractVector, C::AbstractSparseMatrix, d₂::AbstractVector,
                    lb::AbstractVector, ub::AbstractVector,
                    is_binary::AbstractVector, is_integer::AbstractVector,
                    map::IndexMap) where T <: AbstractFloat
        A_ = convert(SparseMatrixCSC{T}, A)
        C_ = convert(SparseMatrixCSC{T}, C)
        # Sanity check
        _lpcheck(c₁, c₂, A, b, d₁, C, d₂, lb, ub, is_binary, is_integer)
        # Return data
        return new{T, SparseMatrixCSC{T,Int}}(c₁, c₂, A_, b, d₁, C_, d₂, lb, ub, is_binary, is_integer, map)
    end
end

LPData(n::Int, m₁::Int, m₂::Int) =
    LPData(Float64, Matrix, n, m₁, m₂)
LPData(t::Type{T}) where T <: AbstractFloat =
    LPData(t, 0, 0, 0)
LPData() =
    LPData(Float64)

SparseLPData(n::Int, m₁::Int, m₂::Int) =
    LPData(Float64, SparseMatrixCSC, n, m₁, m₂)
SparseLPData(t::Type{T}) where T <: AbstractFloat =
    SparseLPData(t, 0, 0, 0)
SparseLPData() =
    SparseLPData(Float64)

function _lpcheck(c₁::AbstractVector, c₂::Real,
                  A::AbstractMatrix, b::AbstractVector,
                  d₁::AbstractVector, C::AbstractMatrix, d₂::AbstractVector,
                  lb::AbstractVector, ub::AbstractVector,
                  is_binary::AbstractVector, is_integer::AbstractVector)
    n₁ = length(c₁)
    m₁, n₂ = size(A)
    m₂, n₃ = size(C)

    if n₁ ≠ n₂
        throw(DomainError(A ,"LPData: `A` must have $(n₁) columns, has $(n₂)"))
    elseif m₁ ≠ length(b)
        throw(DomainError(b ,"LPData: `b` must have $(m₁) elements, , has $(length(b))"))
    elseif n₁ ≠ n₃
        throw(DomainError(D ,"LPData: `D` must have $(n₁) columns, , has $(n₃)"))
    elseif m₂ ≠ length(d₁) || m₂ ≠ length(d₂)
        throw(DomainError(d ,"LPData: `d₁` and `d₂` must have $(m₂) elements, has $(length(d₁)) and $(length(d₂))"))
    elseif n₁ ≠ length(lb) || n₁ ≠ length(ub)
        throw(DomainError((lb, ub) ,"LPData: Both `lb` and `ub` must have $(n₁) elements, , has $(length(lb)) and $(length(ub))"))
    elseif n₁ ≠ length(is_binary) || n₁ ≠ length(is_integer)
        throw(DomainError((lb, ub) ,"LPData: Both `is_binary` and `is_integer` must have $(n₁) elements, , has $(length(is_binary)) and $(length(is_integer))"))
    end
end

function LPData(c₁::AbstractVector, c₂::Real,
                A::AbstractMatrix, b::AbstractVector,
                d₁::AbstractVector, C::AbstractMatrix, d₂::AbstractVector,
                lb::AbstractVector, ub::AbstractVector,
                is_binary::AbstractVector, is_integer::AbstractVector,
                map::IndexMap)
    return LPData(Float64, c₁, c₂, A, b, d₁, C, d₂, lb, ub, is_binary, is_integer, map)
end

function LPData(cor::RawCor{T};
                colrange = collect(keys(cor.cols)),
                rowrange = collect(keys(cor.rows)),
                include_constant = false) where T <: AbstractFloat
    # Compute sizes
    n = length(colrange)
    m₁ = count(x -> x[3] == EQ, [cor.rows[row] for row in rowrange])
    m₂ = length(rowrange) - m₁
    # Prepare index map
    map = IndexMap(n, m₁, m₂)
    # Prepare matrices
    c₁ = zeros(T, n)
    c₂ = zero(T)
    A  = zeros(T, m₁, n)
    b  = zeros(T, m₁)
    d₁ = fill(convert(T, -Inf), m₂)
    C  = zeros(T, m₂, n)
    d₂ = fill(convert(T, Inf), m₂)
    lb = zeros(T, n)
    ub = fill(convert(T, Inf), n)
    is_binary = fill(false, n)
    is_integer = fill(false, n)
    # Aux counters
    ncols = 0
    neqrows = 0
    nineqrows = 0
    for (row, rowval) in cor.rows
        if row in rowrange
            break
        end
        if rowval[3] == EQ
            neqrows += 1
        elseif rowval[3] == LEQ || rowval[3] == GEQ
            nineqrows += 1
        end
    end
    # Fill (in)equalities with ranges
    for (col, val) in cor.cols
        col in colrange || (ncols += 1; continue)
        colidx = cor.vars[col]
        for (row, rowval) in val
            row in rowrange || continue
            idx, rowidx, rowsense = cor.rows[row]
            rowcol = (row, col)
            if rowidx == 0
                j = colidx - ncols
                c₁[j] = rowval
                map[rowcol] = (0,j,OBJ)
            elseif rowsense == EQ
                i = rowidx - neqrows
                j = colidx - ncols
                A[i, j] = rowval
                b[i] = get(cor.rhs, row, zero(T))
                map[rowcol] = (i,j,EQ)
                map[(row,RHS)] = (i,0,EQ)
            elseif rowsense == LEQ
                i = rowidx - nineqrows
                j = colidx - ncols
                C[i, j] = rowval
                d₂[i] = get(cor.rhs, row, zero(T))
                map[rowcol] = (i,j,LEQ)
                map[(row,RHS)] = (i,0,LEQ)
                if haskey(cor.ranges, row)
                    d₁[i] = d₂[i] - abs(cor.ranges[row])
                    map[rowcol] = (i,j,RANGE)
                    map[(row,RHS)] = (i,0,RANGE)
                end
            else
                i = rowidx - nineqrows
                j = colidx - ncols
                C[i, j] = rowval
                d₁[i] = get(cor.rhs, row, zero(T))
                map[rowcol] = (i,j,GEQ)
                map[(row,RHS)] = (i,0,GEQ)
                if haskey(cor.ranges, rowsymbol)
                    d₂[i] = d₁[i] + abs(cor.ranges[row])
                    map[rowcol] = (i,j,RANGE)
                    map[(row,RHS)] = (i,0,RANGE)
                end
            end
        end
    end
    # Get any objective constant
    if cor.objgiven && include_constant
        c₂ = -get(cor.rhs, cor.objname, zero(T))
    end
    # Fill bounds
    fill_bounds!(cor, lb, ub, is_binary, is_integer, colrange, map)
    # Return the problem
    return LPData(T, c₁, c₂, A, b, d₁, C, d₂, lb, ub, is_binary, is_integer, map)
end

function SparseLPData(cor::RawCor{T};
                      colrange = collect(keys(cor.cols)),
                      rowrange = collect(keys(cor.rows)),
                      include_constant = false) where T <: AbstractFloat
    # Compute sizes
    n = length(colrange)
    m₁ = count(x -> x[3] == EQ, [cor.rows[row] for row in rowrange])
    m₂ = length(rowrange) - m₁
    # Prepare index map
    map = IndexMap(n, m₁, m₂)
    # Prepare vectors
    c₁ = zeros(T, n)
    c₂ = zero(T)
    b  = zeros(T, m₁)
    d₁ = fill(convert(T, -Inf), m₂)
    d₂ = fill(convert(T, Inf), m₂)
    lb = zeros(T, n)
    ub = fill(convert(T, Inf), n)
    is_binary = fill(false, n)
    is_integer = fill(false, n)
    # Prepare matrices
    Aᵢ = Vector{Int}()
    Aⱼ = Vector{Int}()
    Aᵥ = Vector{T}()
    Cᵢ = Vector{Int}()
    Cⱼ = Vector{Int}()
    Cᵥ = Vector{T}()
    # Aux counters
    ncols = 0
    neqrows = 0
    nineqrows = 0
    for (row, rowval) in cor.rows
        if row in rowrange
            break
        end
        if rowval[3] == EQ
            neqrows += 1
        elseif rowval[3] == LEQ || rowval[3] == GEQ
            nineqrows += 1
        end
    end
    # Fill (in)equalities with ranges
    for (col, val) in cor.cols
        col in colrange || (ncols += 1; continue)
        colidx = cor.vars[col]
        for (row, rowval) in val
            row in rowrange || row == cor.objname || continue
            idx, rowidx, rowsense = cor.rows[row]
            rowcol = (row, col)
            if rowidx == 0
                j = colidx - ncols
                c₁[j] = rowval
                map[rowcol] = (0,j,OBJ)
            elseif rowsense == EQ
                i = rowidx - neqrows
                j = colidx - ncols
                push!(Aᵢ, i)
                push!(Aⱼ, j)
                push!(Aᵥ, rowval)
                b[i] = get(cor.rhs, row, zero(T))
                map[rowcol] = (i,j,EQ)
                map[(row,RHS)] = (i,0,EQ)
            elseif rowsense == LEQ
                i = rowidx - nineqrows
                j = colidx - ncols
                push!(Cᵢ, i)
                push!(Cⱼ, j)
                push!(Cᵥ, rowval)
                d₂[i] = get(cor.rhs, row, zero(T))
                map[rowcol] = (i,j,LEQ)
                map[(row,RHS)] = (i,0,LEQ)
                if haskey(cor.ranges, row)
                    d₁[i] = d₂[i] - abs(cor.ranges[row])
                    map[rowcol] = (i,j,RANGE)
                    map[(row,RHS)] = (i,0,RANGE)
                end
            else
                i = rowidx - nineqrows
                j = colidx - ncols
                push!(Cᵢ, i)
                push!(Cⱼ, j)
                push!(Cᵥ, rowval)
                d₁[i] = get(cor.rhs, row, zero(T))
                map[rowcol] = (i,j,GEQ)
                map[(row,RHS)] = (i,0,GEQ)
                if haskey(cor.ranges, row)
                    d₂[i] = d₁[i] - abs(cor.ranges[row])
                    map[rowcol] = (i,j,RANGE)
                    map[(row,RHS)] = (i,0,RANGE)
                end
            end
        end
    end
    # Create sparse structures
    A  = sparse(Aᵢ, Aⱼ, Aᵥ, m₁, n)
    C  = sparse(Cᵢ, Cⱼ, Cᵥ, m₂, n)
    # Get any objective constant
    if cor.objgiven && include_constant
        c₂ = -get(cor.rhs, cor.objname, zero(T))
    end
    # Fill bounds
    fill_bounds!(cor, lb, ub, is_binary, is_integer, colrange, map)
    # Return the problem
    return LPData(T, c₁, c₂, A, b, d₁, C, d₂, lb, ub, is_binary, is_integer, map)
end

function fill_bounds!(cor::RawCor{T}, lb, ub, is_binary, is_integer, colrange, map) where T <: AbstractFloat
    ncols = 0
    for (col, val) in cor.bounds
        col in colrange || (ncols += 1; continue)
        colidx = cor.vars[col]
        j = colidx - ncols
        for (bndtype, bndval) in val
            if bndtype == LOWER
                lb[j] = bndval
                map[(LOWER,col)] = (0,j,GEQ)
            elseif bndtype == UPPER
                ub[j] = bndval
                map[(UPPER,col)] = (0,j,LEQ)
            elseif bndtype == FIXED
                lb[j] = ub[j] = bndval
                map[(FIXED,col)] = (0,j,LEQ)
                map[(FIXED,col)] = (0,j,GEQ)
            elseif bndtype == BINARY
                is_binary[j] = true
                lb[j] = convert(T, -Inf)
                map[(BINARY,col)] = (0,j,BINARY)
            elseif bndtype == INTEGER
                is_integer[j] = true
                lb[j] = convert(T, -Inf)
                map[(INTEGER,col)] = (0,j,INTEGER)
            elseif bndtype == INTEGER_LOWER
                is_integer[j] = true
                lb[j] = bndval
                map[(INTEGER,col)] = (0,j,INTEGER_LOWER)
            elseif bndtype == INTEGER_UPPER
                is_integer[j] = true
                ub[j] = bndval
                map[(INTEGER,col)] = (0,j,INTEGER_UPPER)
            else
                lb[j] = convert(T, -Inf)
                map[(FREE,col)] = (0,j,FREE)
            end
        end
    end
    return nothing
end


text: include("models.jl")

ext: jl


text:
to_one(x) = x == true ? 1 : 0
to_neg(x) = x == 0 ? 1 : -1

magnetization(x) = abs(sum(x) / length(x))

function load_generate(T; L=7, n_gibbs=1_000, seed=nothing, n_samples=1)
    # * Create a RNG
    if isnothing(seed)
        rng = Xoroshiro128Plus()
    else
        rng = Xoroshiro128Plus(seed)
    end
    # * Some constants
    L2 = L^2

    # * Allocate to save samples
    samples = Array{Float64}(undef, L, L, n_samples)

    # * Load the respective model for the temperature
    rbm = BernoulliRBM(L2, 32)
    model_name = @sprintf "model_%.2f.jld2" T
    JLD2.@load joinpath("models", model_name) rbm

    # * Create a random state
    confs = Array{Float64}(undef, L, L, n_samples)
    fileising = @sprintf "ising_%.2f.jld2" T
    @load joinpath("data", fileising) confs
    random_sample = confs[:, :, rand(rng, 1:n_samples)]
    random_sample = reshape(random_sample, L2)

    @showprogress "Sampling..." for i = 1:n_samples
        # * Sample until equilibrium
        x_new = generate(rbm, random_sample; n_gibbs=n_gibbs)
        x_new = reshape(x_new, (L, L))

        # * Probabilistiaclly turn visible units on
        x_states = x_new .> rand(rng, Float64, (L, L))
        x_states = to_one.(x_states)

        # * Return to [-1, 1] set of values from the original Ising model
        x_states = to_neg.(x_states)

        # * Save the sample
        samples[:, :, i] = x_states
    end

    return samples
end

function compute_values(;nconfs=10_000, L=8)
    Ts = 1.2:0.1:3.4
    n = size(Ts, 1)
    full_magnetization = Array{Vector{Float64}}(undef, 2, n)
    confs = Array{Float64}(undef, L, L, nconfs)
    ising_magn = Vector{Float64}(undef, n)
    rbm_magn = Vector{Float64}(undef, n)

    for (j, T) in enumerate(Ts)
        confs = Array{Float64}(undef, L, L, nconfs)
        fileising = @sprintf "ising_%.2f.jld2" T
        JLD2.@load joinpath("data", fileising) confs

        rbm_samples = load_generate(T; L=L, n_samples=nconfs)

        for i in axes(confs, 3)
            ising_magn[j] += magnetization(confs[:, :, i])
            rbm_magn[j] += magnetization(rbm_samples[:, :, i])
        end
    end

    ising_magn ./= size(confs, 3)
    display(ising_magn)
    JLD2.@save joinpath("results", "ising_magnetization.jld2") ising_magn
    rbm_magn ./= size(confs, 3)
    display(rbm_magn)
    JLD2.@save joinpath("results", "rbm_magnetization.jld2") rbm_magn
end

# compute_values()

max_stars_repo_name: edwinb-ai/IsingRBM

max_issues_repo_name: edwinb-ai/IsingRBM

max_forks_repo_name: edwinb-ai/IsingRBM


text:
__precompile__(false)
module TicTacToe

export TicTacToeEnv

using ReinforcementLearningBase

const REWARD_LOSE = 0.0
const REWARD_UNFINISH = 0.0
const REWARD_WIN = 1.0
const REWARD_TIE = 0.5
const IDLE_ACTION = 10

struct Offensive end
const offensive = Offensive()
struct Defensive end
const defensive = Defensive()

const ROLES = (offensive, defensive)
const Board = Array{Union{Offensive,Defensive,Nothing}}
const Roles = Union{Offensive,Defensive}

const CHECK_INDS = vcat(
    [(i, i + 3, i + 6) for i = 1:3],     # rows
    [(3i - 2, 3i - 1, 3i) for i = 1:3],  # cols
    [(1, 5, 9), (3, 5, 7)],
)             # diag

function get_winner(board)
    for inds in CHECK_INDS
        w = get_winner(board, inds)
        w ≢ nothing && return w
    end
    nothing
end

get_winner(board, inds) =
    board[inds[1]] == board[inds[2]] == board[inds[3]] ? board[inds[1]] : nothing

RLBase.get_legal_actions_mask(state::Board) = state .=== nothing

function get_next_state(state::Board, role::Roles, action)
    s = copy(state)
    if action != IDLE_ACTION
        s[action] = role
    end
    s
end

get_next_states(state::Board, role::Roles, actions = findall(get_legal_actions_mask(state))) =
    (get_next_state(state, role, action) for action in actions)

get_next_state_id(id::Int, role::Roles, action) = STATE2ID[get_next_state(ID2STATE[id], role, action)]

get_next_role(::Offensive) = defensive
get_next_role(::Defensive) = offensive
get_next_role(::Nothing) = offensive  # offensive role starts the game

function get_states_info()
    init_state, states_info, unfinished_states = Board(nothing, 3, 3), Dict(), Set()
    states_info[init_state] = (isdone = false, winner = nothing)
    push!(unfinished_states, (init_state, offensive))
    while length(unfinished_states) > 0
        state, role = pop!(unfinished_states)
        for s in get_next_states(state, role)
            if !haskey(states_info, s)
                winner = get_winner(s)
                isdone = winner != nothing || sum(get_legal_actions_mask(s)) == 0
                states_info[s] = (isdone = isdone, winner = winner)
                if !isdone
                    push!(unfinished_states, (s, get_next_role(role)))
                end
            end
        end
    end
    states_info
end

const STATES_INFO = get_states_info()
const STATE2ID = Dict(s => i for (i, s) in enumerate(keys(STATES_INFO)))
const ID2STATE = Dict(i => s for (i, s) in enumerate(keys(STATES_INFO)))

#####
# TicTacToeEnv
#####

"""
    TicTacToeEnv()

Using a 3 * 3 Array to simulate the [tic-tac-toe](https://en.wikipedia.org/wiki/Tic-tac-toe) game.
"""
mutable struct TicTacToeEnv <: AbstractEnv
    role::Union{Nothing,Roles}
    board::Board
    observation_space::DiscreteSpace
    action_space::DiscreteSpace
    function TicTacToeEnv()
        init_board = Board(nothing, 3, 3)
        new(
            nothing,
            init_board,
            DiscreteSpace(length(STATES_INFO)),
            DiscreteSpace(length(init_board) + 1),
        )
    end
end

RLBase.get_observation_space(env::TicTacToeEnv) = env.observation_space
RLBase.get_action_space(env::TicTacToeEnv) = env.action_space
RLBase.get_current_player(env::TicTacToeEnv) = get_next_role(env)
RLBase.get_num_players(env::TicTacToeEnv) = 2

get_winner(env::TicTacToeEnv) = STATES_INFO[env.board].winner
is_done(env::TicTacToeEnv) = STATES_INFO[env.board].isdone

RLBase.observe(env::TicTacToeEnv) = RLBase.observe(env, get_next_role(env))

function RLBase.observe(env::TicTacToeEnv, role::Roles)
    isdone, winner = STATES_INFO[env.board]
    if isdone
        if winner == nothing
            reward = REWARD_TIE
        elseif winner == role
            reward = REWARD_WIN
        else
            reward = REWARD_LOSE
        end
    else
        reward = REWARD_UNFINISH
    end

    (
        reward = reward,
        terminal = isdone,
        state = STATE2ID[env.board],
        legal_actions_mask = get_legal_actions_mask(env, role),
    )
end

function RLBase.reset!(env::TicTacToeEnv)
    fill!(env.board, nothing)
    env.role = nothing
    nothing
end

get_next_role(env::TicTacToeEnv) = is_done(env) ? nothing : get_next_role(env.role)

function (env::TicTacToeEnv)(action::Int)
    next_role = get_next_role(env)
    another_role = get_next_role(next_role)
    env([next_role, another_role] => [action, IDLE_ACTION])
end

function (env::TicTacToeEnv)(act_info::Pair{<:Vector,<:Vector})
    is_done(env) && throw(ArgumentError("env is already done!"))
    roles, actions = act_info
    nextrole = get_next_role(env)

    for (role, action) in zip(roles, actions)
        if role === nextrole
            env.board[action] === nothing || throw(ArgumentError("invalid action: $action"))
            env.role = role
            env.board[action] = role
        else
            action != IDLE_ACTION && throw(ArgumentError("invalid action [$action] for role [$role]"))
        end
    end

    nothing
end

get_roles(::TicTacToeEnv) = ROLES

function RLBase.get_legal_actions_mask(env::TicTacToeEnv, role)
    legal_actions = fill(false, IDLE_ACTION)
    if role === get_next_role(env)
        for i = 1:9
            legal_actions[i] = env.board[i] === nothing
        end
    else
        legal_actions[end] = true
    end
    legal_actions
end

function Base.show(io::IO, env::TicTacToeEnv)
    for r = 1:3
        for c = 1:3
            s = env.board[r, c]
            print(io, s === nothing ? "_" : s)
        end
        println(io)
    end
    println(io, "isdone = [$(is_done(env))], winner = [$(get_winner(env))]")
end

Base.show(io::IO, ::Offensive) = print(io, "X")
Base.show(io::IO, ::Defensive) = print(io, "O")

end


text:
# MIT License

# Copyright (c) 2020 Makoto Suwama

# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

# Code used to compute cycles coming from the sum of squares of digits in base b for the paper Integer Dynamics with Dino Lorenzini, Mentzelos Melistas, Arvind Suresh and Haiyang Wang
# For example, to check in base b, all cycles starting at x+b*y with low<= y <= high, run
# julia cycles.jl low high b
# May need to change the data type of b from UInt32 to other type depending on how large the b is
# Also see correctness.pdf for the correctness of the code.

const numCycles=3;

# Parameters parsed from the command line
low = parse(Int, ARGS[1]);
high = parse(Int, ARGS[2]);
const b = parse(Int, ARGS[3]);

# Directory to save the files
directory = "";
output = directory * "spTestJ" * string(b) * "-" * string(low) * "-" * string(high) * ".txt";

# check a number d has been seen already
# seen is an array of arrays of sets and a number n = d1 + d2*b + d3*b^2 is stored as d1 in seen[d3+1][d2+1]
function isSeen(seen,d)
    return d[1] in seen[d[3]+1][d[2]+1];
end;
# add number d to seen
# seen is an array of arrays of sets and a number n = d1 + d2*b + d3*b^2 is stored as d1 in seen[d3+1][d2+1]
function addSeen!(seen,d)
    push!(seen[d[3]+1][d[2]+1], d[1])
end;
# replaces d with [x,y,0]
function createD!(d,x,y)
    d[1] = x;
    d[2] = y;
    d[3] = 0;
end
# Checks if max(d[1],d[2]) + min(d[1],d[2])*b + d[2]*b^2 is smaller than x+b*y
# Uses separate comparison for speed purpose
function isSmaller(d,x,y)
    return d[3] == 0 && (min(d[1],d[2]) < y || (min(d[1],d[2])==y && max(d[1],d[2]) < x))
end

# start the sum of the square procedure until it hits a cycle or a number already seen
# Returns True if it finds a new cycle
function traverse(x,y,d,seen,cycles,L)
    # if the number already appeared then skip it
    if x in seen[1][y+1]
        return false
    end;
    n = x+b*y;

    # pops L until it is empty and adds n
    # this is faster than L = [n]
    while length(L) > 0
        pop!(L);
    end;
    push!(L,n);
    while !(isSeen(seen,d))
        addSeen!(seen,d)
        n = sum(i^2 for i in d);
        push!(L,n)
        digits!(d, n,base = b)
        # if d is less than x+b*y then we have seen it already
        if isSmaller(d,x,y)
            return false
        end
    end
    # find the position of n in L
    # if n is in L and not the last then there is a cycle starting at the first appearance of n
    p = findfirst(isequal(n),L);
    if p != nothing && p < length(L)
        push!(cycles, L[p:length(L)-1])
        return true
    end
    return false
end

# Check all the numbers x+b*y with low<=y<=high
function checkYInterval(low, high, seen, cycles)
    d=Array{UInt64}(undef,3)
    L=[];
    for y=low:high
        # Check if the progress needs to be printed
        if y % interval == 0
            println(y, " ", time()-t0)
            flush(stdout)
        end
        # See correctness.pdf for the upper bound on x
        for x=max(y,ceil(Int, sqrt(y*b-y^2+1/4)+1/2)):b-1
            createD!(d,x,y)
            if traverse(x,y,d,seen,cycles,L)
                #if length(cycles) >= numCycles
                #    return nothing
                #end
            end
        end
        # delete unncessary seen to save memory. Boolean is chosen to save space
        seen[1][y+1] = Set{Bool}()
    end
end

# Adds a cycle starting at ini
function addKnownCycle(seen, cycles, ini, d)
    c = [];
    while true
        digits!(d,ini,base = b)
        if isSeen(seen, d)
            push!(cycles, c)
            return
        end
        addSeen!(seen,d)
        push!(c, ini)
        ini = sum(i^2 for i in d)
    end
end

# Initialise the variables
# interval is used to determine how often should we print progress
const t0 = time();
const interval =  10^ceil(Int,log(10,(high-low)/100));

# b is UInt32
# indexing is off by 1
# if n = [a,b,c] then check if a is in seen[c+1][b+1]
# keeps track of [a,b,c] where a>=b and c=0 or 1
seen = [[ Set{UInt32}() for i=1:b], [Set{UInt32}() for i=1:b]];
cycles = [[1]];
push!(seen[1][1],1);
d=Array{UInt64}(undef,3);

addKnownCycle(seen, cycles, 91610640434,d);

checkYInterval(low, high, seen, cycles)

io = open(output,"w");
println(io, b);
println(io, [length(c) for c in cycles]);
println(io, [c[1] for c in cycles]);
flush(io);
close(io);

max_stars_repo_name: msuwama/IntegerDynamics

max_issues_repo_name: msuwama/IntegerDynamics

max_forks_repo_name: msuwama/IntegerDynamics


text:
module DeconvolutionTests

import DeconvolveDistribution
import DeconvolveDistribution:
        decon, DeconvolutionMethod,
        FourierDeconv, EfronDeconv
using HypothesisTests
using Distributions
using StatsBase: sample, midpoints

# include("deconvolve.jl")
include("deconv_boot_test.jl")
include("test_statistics.jl")
include("simulate.jl")

export deconv_boot_test, distance_test_statistic, KS_test_statistic

end # module


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

# # Representing data in a computer
#
# The core of data science and machine learning is **data**: we are interested in extracting knowledge from data.
#
# But how exactly do computers represent data? Let's find out exactly what an "artificial intelligence" has at its disposal to learn from.

#-

# ## Data is represented as arrays
#
# Let's take a look at some fruit. Using the `Images.jl` library, we can load in some images:

using Images

apple = load(datapath("data/10_100.jpg"))

#-

banana = load(datapath("data/104_100.jpg"))

# Here we have images of apples and bananas. We would eventually like to build a program that can automatically distinguish between the two. However, the computer doesn't "see" an apple or a banana; instead, it just sees numbers.
#
# An image is encoded in something called an **array**, which is like a container that has boxes or slots for individual pieces of data:
#
# An array is a bunch of numbers in connected boxes; the figure above shows a 1-dimensional array. Our images are instead 2-dimensional arrays, or matrices, of numbers, arranged something like this:

#-

#

#-

# For example, `apple` is an image, consisting of a 100x100 array of numbers:

typeof(apple)

#-

size(apple)

#-

a = [ 1 2 3;4 5 6]

#-

typeof(a)

#-

size(a)

# We can grab the datum stored in the box at row `i` and column `j` by *indexing* using square brackets: `[i, j]`. For example, let's get the pixel (piece of the image) in box $(40, 60)$, i.e. in the 40th row and 60th column of the image:

apple

#-

dump(apple[40, 60])

#-

apple[18:20,29:31]

# We see that Julia displays a coloured box! Julia, via the `Colors.jl` package, is clever enough to display colours in a way that is useful to us humans!
#
# So, in fact, an image is a 2D array, in which each element of the array is an object (a collection of numbers) describing a coloured pixel.

#-

# ## Colors as numbers
#
# How, then, are these colors actually stored? Computers store colors in RGB format, that is they store a value between 0 and 1 for each of three "channels": red, green, and blue. Here, 0 means none of that color and 1 means the brightest form of that color. The overall color is a combination of those three colors.
#
# For example, we can pull out the `red` value using the function `red` applied to the color. Since internally the actual value is stored in a special format, we choose to convert it to a standard floating-point number using the `Float64` function:

Float64(red(apple[40, 60]))

#-

using Statistics

#-

[ mean(float.(c.(img))) for c = [red,green,blue], img = [apple,banana] ]

#-

using Plots
histogram(float.(green.(apple[:])),color="red",label="apple", normalize=true, nbins=25)
histogram!(float.(green.(banana[:])),color="yellow",label="banana",normalize=true, nbins=25)

#-

apple

#-

float(red(banana[50,20]))

#-

banana[50,20]

#-

pixel = apple[40, 60]

red_value   = Float64( red(pixel) )
green_value = Float64( green(pixel) )
blue_value  = Float64( blue(pixel) )

print("The RGB values are ($red_value, $green_value, $blue_value)")

# Since the red value is high while the others are low, this means that at pixel `(40, 60)`, the picture of the apple is very red. If we do the same at one of the corners of the picture, we get the following:

pixel = apple[1, 1]

red_value   = Float64( red(pixel) )
green_value = Float64( green(pixel) )
blue_value  = Float64( blue(pixel) )

print("The RGB values are ($red_value, $green_value, $blue_value)")

#-

apple

# We see that every color is bright, which corresponds to white.

#-

# ## Working on an image as a whole
#
# In Julia, to apply a function to the whole of an array, we place a `.` between the function name and the left parenthesis (`(`), so the following gives us the `red` value of every pixel in the image:

redpartofapple = Float64.(red.(apple))
mean(redpartofapple)

#-

using Plots

#-

gr()

#-

histogram(redpartofapple[:],color=:red,label="redness in the apple")

# Note that we get a 2D array (matrix) back.

#-

# Julia's [mathematical standard library](https://docs.julialang.org/en/stable/stdlib/math/#Mathematics-1) has many mathematical functions built in. One of them is the `mean` function, which computes the average value. If we apply this to our apple:

mean(Float64.(red.(apple)))

# we see that the value indicates that the average amount of red in the image is a value between the amount of red in the apple and the amount of red in the white background.
#
# *Somehow we need to teach a computer to use this information about a picture to recognize that there's an apple there!*

#-

# ## A quick riddle
#
# Here's a quick riddle. Let's check the average value of red in the image of the banana.

mean(Float64.(red.(banana)))

# Oh no, that's more red than our apple? This isn't a mistake and is actually true! Before you move onto the next exercise, examine the images of the apple and the banana very carefully and see if you can explain why this is expected.

#-

# #### Exercise 1
#
# What is the average value of blue in the banana?
#
# (To open a new box use <ESC>+b (b is for "below", what do you think a does?))

#-

# #### Solution
#
# We can calculate the average value of blue in the banana via

mean(Float64.(blue.(banana)))

# The result is approximately `0.8`.

#-

# #### Exercise 2
#
# Does the banana have more blue or more green?

#-

# #### Solution
#
# The average value of green in the banana is

mean(Float64.(green.(banana)))

# which gives approximately `0.88`. The banana has more green on average.

#-

# # Modeling data 1
#
# Machine learning and data science is about modeling data. **Modeling** is the representation of an idea with some parameters and  a mathematical representation which we will encode in software. All machine learning methods are about training a computer to fit a model to some data. Even the fanciest neural networks are simply choices for models. In this notebook, we will begin to start building our first computational model of data.

#-

# ## Modeling data is hard!
#
# Let's pick up where we left off in notebook 1 with fruit. We were left with a riddle: when we load images of apples and bananas,

## using Pkg; Pkg.add("Images")
using Images, Statistics

apple = load(datafile("data/10_100.jpg"))

#-

banana = load(datafile("data/104_100.jpg"))

# and then compare their average value for the color red, we end up with something that is perhaps surprising:

apple_red_amount = mean(Float64.(red.(apple)))
banana_red_amount = mean(Float64.(red.(banana)));

"The average value of red in the apple is $apple_red_amount, " *
"while the average value of red in the banana is $banana_red_amount."

# We see that the banana's mean red value is higher than the apple's, even though the apple looks much redder. Can you guess why?
#
# There are actually two reasons. One of the reasons is the background: the image of the banana has a lot more white background than the apple, and that white background has a red value of 1! In our minds we ignore the background and say "the banana is bright yellow, the apple is dark red", but a computer just has a bundle of numbers and does not know where it should be looking.
#
# The other issue is that "bright yellow" isn't a color that exists in a computer. The computer has three colors: red, green, and blue. "Bright yellow" in a computer is a mixture of red and green, and it just so happens that to get this color yellow, it needs more red than the apple!

"The amount of red in the apple at (60, 60) is $(Float64(red(apple[60, 60]))), " *
"while the amount of red in the banana at (60, 60) is $(Float64(red(banana[60, 60])))."

#-

apple[60,60]

#-

banana[60,60]

# This is a clear example that modeling data is hard!

#-

# ### A note on string interpolation
#
# In the last two input cells, we *interpolated a string*. This means that when we write the string using quotation marks (`"  "`), we insert a placeholder for some **value** we want the string to include. When the string is evaluated, the value we want the string to include replaces the placeholder. For example, in the following string,
#
# ```julia
# mystring = "The average value of red in the apple is $apple_red_amount"
# ```
#
# `$apple_red_amount` is a placeholder for the value stored in the variable `apple_red_amount`. Julia knows that we want to use the value bound to the variable `apple_red_amount` and *not* the word "apple_red_amount" because of the dollar sign, `$`, that comes before `apple_red_amount`.

#-

# #### Exercise 1
#
# Execute the following code to see what the dollar sign does:
#
# ```julia
# mypi = 3.14159
# println("I have a variable called mypi that has a value of $mypi.")
# ```

#-

# #### Solution

mypi = 3.14159
println("I have a variable called mypi that has a value of $mypi.")

# #### Exercise 2
#
# Alter and execute the code that creates `mystring` below
#
# ```julia
# apple_blue_amount = mean(Float64.(blue.(apple)))
# mystring = "The average amount of blue in the apple is apple_blue_amount"
# ```
#
# so that `println(mystring)` prints a string that reports the mean value of blue coloration in our image of an apple.

#-

# #### Solution

#-

# Add a `$` in front of `apple_blue_amount`:

apple_blue_amount = mean(Float64.(blue.(apple)))
mystring = "The average amount of blue in the apple is $apple_blue_amount."

# ## Take some time to think about the data
#
# Apples and bananas are very different, but how could we use the array of RGB values (which is how the images are represented in the computer, as we saw in notebook 1) to tell the difference between the two? Here are some quick ideas:
#
# - We could use the shape of the object in the image. But how can we encode ideas about shape from an array?
# - We could use the size of the object in the image. But how do we calculate that size?
# - We could use another color, or combinations of colors, from the image. Which colors?
#
# Let's go with the last route. The banana is yellow, which is a combination of red and green, while the apple is red. This means that the color that clearly differentiates between the two is not red, but green!

apple_green_amount = mean(Float64.(green.(apple)))
banana_green_amount = mean(Float64.(green.(banana)));

"The average value of green in the apple is $apple_green_amount, " *
"while the average value of green in the banana is $banana_green_amount."

# The processes that we just went through are assigned fancy names: feature selection and data munging.
#
# **Feature selection** is the process of subsetting the data to a more relevant and informative set. We took the full image data and decided to select out the green channel.
#
# **Data munging** is transforming the data into a format more suitable for modeling. Here, instead of keeping the full green channel, we transformed it down to a single data point: the average amount of green.

#-

# ## Building a model
#
# We want to model the connection between "the average amount of green" and "is an apple or banana".
#
#
#
# This model is a mathematical function which takes in our data and spits out a number that we will interpret as "is an apple" or "is a banana".
#
#
#
#
# We will interpret the output of the function as "is an apple" if the output is close to 0, and "is a banana" if it's close to 1. Anything in the middle is something we are unsure about. Here we're using a mathematical function to perform a **classification**.
#
# Knowing how to declare and work with functions will allow us to model our data in the coming sections, so this is the subject of the next notebook!

max_stars_repo_path: Courses/Foundations of machine learning/0300.Representing-data-with-models.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
module TransientFETests

using Gridap
using Test
using GridapODEs.ODETools
using GridapODEs.TransientFETools
using Gridap.FESpaces: get_algebraic_operator

u(x,t) = (x[1] + x[2])*t
u(t::Real) = x -> u(x,t)
∇u(x,t) = VectorValue(1,1)*t
∇u(t::Real) = x -> ∇u(x,t)
import Gridap: ∇
∇(::typeof(u)) = ∇u
∇(u) === ∇u

θ = 1.0

∂tu(t) = x -> x[1]+x[2]
import GridapODEs.TransientFETools: ∂t
∂t(::typeof(u)) = ∂tu
@test ∂t(u) === ∂tu

f(t) = x -> (x[1]+x[2])

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V0 = TestFESpace(
  model,
  reffe,
  conformity=:H1,
  dirichlet_tags="boundary"
)

U = TransientTrialFESpace(V0,u)
U0 = TrialFESpace(V0,u(0.0))
@test test_transient_trial_fe_space(U)

U0 = U(1.0)
ud0 = copy(get_dirichlet_dof_values(U0))
_ud0 = get_dirichlet_dof_values(U0)
U1 = U(2.0)
ud1 = copy(get_dirichlet_dof_values(U1))
_ud1 = get_dirichlet_dof_values(U1)
@test all(ud0 .≈ 0.5ud1)
all(_ud0 .≈ _ud1)

Ut = ∂t(U)
Ut.dirichlet_t
Ut0 = Ut(0.0)
Ut0.dirichlet_values

Ut1 = Ut(1.0)
utd0 = copy(get_dirichlet_dof_values(Ut0))
utd1 = copy(get_dirichlet_dof_values(Ut1))
@test all(utd0 .== utd1)
@test all(utd1 .== ud0)

Ω = Triangulation(model)
degree = 2
dΩ = Measure(Ω,degree)

a(u,v) = ∫(∇(v)⋅∇(u))dΩ
b(v,t) = ∫(v*f(t))dΩ

res(t,u,v) = a(u,v) + ∫(∂t(u)*v)dΩ - b(v,t)
jac(t,u,du,v) = a(du,v)
jac_t(t,u,dut,v) = ∫(dut*v)dΩ

U0 = U(0.0)
_res(u,v) = a(u,v) + 10.0*∫(u*v)dΩ - b(v,0.0)
_jac(u,du,v) = a(du,v) + 10.0*∫(du*v)dΩ
_op = FEOperator(_res,_jac,U0,V0)

uh = interpolate_everywhere(0.0,U0)#1.0)
using Gridap.FESpaces: allocate_residual, allocate_jacobian
_r = allocate_residual(_op,uh)
_J = allocate_jacobian(_op,uh)
using Gridap.FESpaces: residual!, jacobian!
residual!(_r,_op,uh)
jacobian!(_J,_op,uh)

op = TransientFEOperator(res,jac,jac_t,U,V0)
odeop = get_algebraic_operator(op)
cache = allocate_cache(odeop)

r = allocate_residual(op,uh,cache)
J = allocate_jacobian(op,uh,cache)
uh10 = interpolate_everywhere(0.0,U0)#10.0)
xh = TransientCellField(uh,(uh10,))
residual!(r,op,0.0,xh,cache)
jacobian!(J,op,1.0,xh,1,1.0,cache)
jacobian!(J,op,1.0,xh,2,10.0,cache)
@test all(r.≈_r)
@test all(J.≈_J)

U0 = U(0.0)
uh0 = interpolate_everywhere(0.0,U0)
@test test_transient_fe_operator(op,uh0)

u0 = u(0.0)
t0 = 0.0
tF = 1.0
dt = 0.1

ls = LUSolver()
# using LineSearches: BackTracking
tol = 1.0
maxiters = 20
using Gridap.Algebra: NewtonRaphsonSolver
nls = NLSolver(ls;show_trace=true,method=:newton) #linesearch=BackTracking())
ode_solver = ThetaMethod(nls,dt,1.0)
@test test_transient_fe_solver(ode_solver,op,uh0,t0,tF)

xh = TransientCellField(uh,(uh,))
residual!(r,op,0.1,xh,cache)
jacobian!(J,op,1.0,xh,1,1.0,cache)
jacobian!(J,op,1.0,xh,2,10.0,cache)

u0 = get_free_dof_values(uh0)
solver = ode_solver
t0 = 0.0
ode_cache = allocate_cache(odeop)
cache = nothing
uf = copy(u0)
dt = solver.dt
tf = t0+dt
update_cache!(ode_cache,odeop,tf)
using GridapODEs.ODETools: ThetaMethodNonlinearOperator
vf = copy(u0)
nlop = ThetaMethodNonlinearOperator(odeop,tf,dt,u0,ode_cache,vf)

x = copy(nlop.u0)

b1 = allocate_residual(nlop,x)
residual!(b1,nlop,x)
b2 = allocate_residual(nlop,x)
residual!(b2,nlop.odeop,nlop.tθ,(x,10.0*x),nlop.ode_cache)
@test all(b1 .≈ b2)
J1 = allocate_jacobian(nlop,x)
jacobian!(J1,nlop,x)
J2 = allocate_jacobian(nlop,x)
jacobian!(J2,nlop.odeop,nlop.tθ,(x,10.0*x),1,1.0,nlop.ode_cache)
jacobian!(J2,nlop.odeop,nlop.tθ,(x,10.0*x),2,10.0,nlop.ode_cache)
@test all(J1 .≈ J2)
using Gridap.Algebra: test_nonlinear_operator
test_nonlinear_operator(nlop,x,b1,jac=J1)

x .= 0.0
r = allocate_residual(nlop,x)
residual!(r,nlop,x)
J = allocate_jacobian(nlop,x)
jacobian!(J,nlop,x)

cache = solve!(uf,solver.nls,nlop)
df = cache.df
ns = cache.ns

function linsolve!(x,A,b)
  numerical_setup!(ns,A)
  solve!(x,ns,b)
end

p = copy(x)
p .= 0.0
l_sol = linsolve!(p,J,-r)
J*l_sol .≈ -r
x = x + l_sol
@test all(abs.(residual!(r,nlop,x)) .< 1e-6)

residual!(r,nlop,x)
jacobian!(J,nlop,x)
p .= 0.0
l_sol = linsolve!(p,J,-r)

cache = solve!(uf,solver.nls,nlop)
@test all(uf .≈ x)
solve!(uf,solver.nls,nlop,cache)
@test all(uf .≈ x)

uf .= 0.0
x = copy(nlop.u0)
cache = Gridap.Algebra._new_nlsolve_cache(x,nls,nlop)
df = cache.df
ns = cache.ns
x .= 0.0
l_sol = linsolve!(x,df.DF,df.F)
@test all(df.DF*l_sol.≈df.F)
x .= 0
Gridap.Algebra.nlsolve(df,x;linsolve=linsolve!,nls.kwargs...)

using Gridap.FESpaces: get_algebraic_operator
odeop = get_algebraic_operator(op)
sol_ode_t = solve(ode_solver,odeop,u0,t0,tF)

test_ode_solution(sol_ode_t)
_t_n = t0
for (u_n, t_n) in sol_ode_t
  global _t_n
  _t_n += dt
  @test t_n≈_t_n
  @test all(u_n .≈ t_n)
end

ode_solver = ThetaMethod(nls,dt,θ)
sol_ode_t = solve(ode_solver,odeop,u0,t0,tF)
test_ode_solution(sol_ode_t)
_t_n = t0
un, tn = Base.iterate(sol_ode_t)
for (u_n, t_n) in sol_ode_t
  global _t_n
  _t_n += dt
  @test t_n≈_t_n
  @test all(u_n .≈ t_n)
end

sol_t = solve(ode_solver,op,uh0,t0,tF)
@test test_transient_fe_solution(sol_t)

_t_n = 0.0
for (u_n, t_n) in sol_t
  global _t_n
  _t_n += dt
  @test t_n≈_t_n
  @test all(u_n.free_values .≈ t_n)
end

l2(w) = w*w

# h1(w) = a(w,w) + l2(w)

_t_n = t0
for (uh_tn, tn) in sol_t
  global _t_n
  _t_n += dt
  @test tn≈_t_n
  e = u(tn) - uh_tn
  el2 = sqrt(sum( ∫(l2(e))dΩ ))
  @test el2 < tol
  # writevtk(trian,"sol at time: $tn",cellfields=["u" => uh_tn])
end

end #module


text:
"""
    update!(
        KS::X,
        ctr::Y,
        face::Z,
        dt,
        residual; # 1D / 2D
        coll = :bgk::Symbol,
        bc = :fix::Symbol,
    ) where {
        X<:AbstractSolverSet,
        Y<:AbstractArray{ControlVolume1D1F,1},
        Z<:AbstractArray{Interface1D1F,1},
    }

    update!(
        KS::X,
        ctr::Y,
        face::Z,
        dt,
        residual; # 1D / 2D
        coll = :bgk::Symbol,
        bc = :extra::Symbol,
    ) where {
        X<:AbstractSolverSet,
        Y<:AbstractArray{ControlVolume1D2F,1},
        Z<:AbstractArray{Interface1D2F,1},
    }

    update!(
        KS::X,
        ctr::Y,
        face::Z,
        dt,
        residual; # 1D / 2D
        coll = :bgk::Symbol,
        bc = :extra::Symbol,
        isMHD = true::Bool,
    ) where {
        X<:AbstractSolverSet,
        Y<:AbstractArray{ControlVolume1D3F,1},
        Z<:AbstractArray{Interface1D3F,1},
    }

    update!(
        KS::X,
        ctr::Y,
        face::Z,
        dt,
        residual; # 1D / 2D
        coll = :bgk::Symbol,
        bc = :extra::Symbol,
        isMHD = true::Bool,
    ) where {
        X<:AbstractSolverSet,
        Y<:AbstractArray{ControlVolume1D4F,1},
        Z<:AbstractArray{Interface1D4F,1},
    }

    update!(
        KS::X,
        ctr::Y,
        a1face::Z,
        a2face::Z,
        dt,
        residual; # 1D / 2D
        coll = :bgk::Symbol,
        bc = :fix::Symbol,
    ) where {
        X<:AbstractSolverSet,
        Y<:AbstractArray{ControlVolume2D2F,2},
        Z<:AbstractArray{Interface2D2F,2},
    }

Update flow variables over control volumes

"""
function update!(
    KS::X,
    ctr::Y,
    face::Z,
    dt,
    residual::T; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolume1D,1},
    Z<:AbstractArray{Interface1D,1},
    T<:AbstractFloat,
}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds Threads.@threads for i = 1:KS.pSpace.nx
        ctr[i].w, sumRes, sumAvg = step!(
            face[i].fw,
            ctr[i].w,
            ctr[i].prim,
            face[i+1].fw,
            KS.gas.a,
            KS.ps.dx[i],
            sumRes,
            sumAvg,
        )
    end

    residual = sqrt(sumRes * KS.pSpace.nx) / (sumAvg + 1.e-7)

    if bc[1] == :period
        ctr[0].w = ctr[KS.pSpace.nx].w
        ctr[0].sw = ctr[KS.pSpace.nx].sw
        ctr[KS.pSpace.nx+1].w = ctr[1].w
        ctr[KS.pSpace.nx+1].sw = ctr[1].sw
    elseif bc[1] == :extra
        ctr[0].w = ctr[1].w
    end

    if bc[2] == :extra
        ctr[KS.pSpace.nx+1].w = ctr[KS.pSpace.nx].w
    end

    return residual

end

function update!(
    KS::X,
    ctr::Y,
    face::Z,
    dt,
    residual::AbstractArray{T}; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolume1D,1},
    Z<:AbstractArray{Interface1D,1},
    T<:AbstractFloat,
}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    if ndims(sumRes) == 1
        @inbounds Threads.@threads for i = 2:KS.pSpace.nx-1
            step!(
                face[i].fw,
                ctr[i].w,
                ctr[i].prim,
                face[i+1].fw,
                KS.gas.γ,
                KS.ps.dx[i],
                sumRes,
                sumAvg,
            )
        end
    elseif ndims(sumRes) == 2
        @inbounds Threads.@threads for i = 2:KS.pSpace.nx-1
            step!(
                face[i].fw,
                ctr[i].w,
                ctr[i].prim,
                face[i+1].fw,
                KS.gas.γ,
                KS.gas.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
                KS.ps.dx[i],
                dt,
                sumRes,
                sumAvg,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * KS.pSpace.nx) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(KS, ctr, face, dt, residual; bc = [bc, bc])
    else
        update_boundary!(KS, ctr, face, dt, residual; bc = bc)
    end

    return nothing

end

function update!(
    KS::X,
    ctr::Y,
    face::Z,
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolume1D1F,1},
    Z<:AbstractArray{Interface1D1F,1},
}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    phase = KS.set.space[3:end]
    if phase == "1f1v"
        @inbounds Threads.@threads for i = 2:KS.pSpace.nx-1
            step!(
                face[i].fw,
                face[i].ff,
                ctr[i].w,
                ctr[i].prim,
                ctr[i].f,
                face[i+1].fw,
                face[i+1].ff,
                KS.vSpace.u,
                KS.vSpace.weights,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.ps.dx[i],
                dt,
                sumRes,
                sumAvg,
                coll,
            )
        end
    elseif phase == "1f3v"
        @inbounds Threads.@threads for i = 2:KS.pSpace.nx-1
            step!(
                face[i].fw,
                face[i].ff,
                ctr[i].w,
                ctr[i].prim,
                ctr[i].f,
                face[i+1].fw,
                face[i+1].ff,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.w,
                KS.vSpace.weights,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.ps.dx[i],
                dt,
                sumRes,
                sumAvg,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * KS.pSpace.nx) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = [bc, bc])
    else
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc)
    end

    return nothing

end

function update!(
    KS::X,
    ctr::Y,
    face::Z,
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolume1D2F,1},
    Z<:AbstractArray{Interface1D2F,1},
}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    if ndims(sumRes) == 1
        @inbounds Threads.@threads for i = 2:KS.pSpace.nx-1
            step!(
                face[i].fw,
                face[i].fh,
                face[i].fb,
                ctr[i].w,
                ctr[i].prim,
                ctr[i].h,
                ctr[i].b,
                face[i+1].fw,
                face[i+1].fh,
                face[i+1].fb,
                KS.vSpace.u,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.ps.dx[i],
                dt,
                sumRes,
                sumAvg,
                coll,
            )
        end
    elseif ndims(sumRes) == 2
        @inbounds Threads.@threads for i = 2:KS.pSpace.nx-1
            step!(
                face[i].fw,
                face[i].fh,
                face[i].fb,
                ctr[i].w,
                ctr[i].prim,
                ctr[i].h,
                ctr[i].b,
                face[i+1].fw,
                face[i+1].fh,
                face[i+1].fb,
                KS.vSpace.u,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
                KS.gas.Pr,
                KS.ps.dx[i],
                dt,
                sumRes,
                sumAvg,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * KS.pSpace.nx) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = [bc, bc])
    else
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc)
    end

    return nothing

end

function update!(
    KS::X,
    ctr::Y,
    face::Z,
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
    isMHD = true::Bool,
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolume1D3F,1},
    Z<:AbstractArray{Interface1D3F,1},
}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds Threads.@threads for i = 2:KS.pSpace.nx-1
        step!(KS, face[i], ctr[i], face[i+1], KS.ps.dx[i], dt, sumRes, sumAvg, coll, isMHD)
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * KS.pSpace.nx) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(
            KS,
            ctr,
            face,
            dt,
            residual;
            coll = coll,
            bc = [bc, bc],
            isMHD = isMHD,
        )
    else
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc, isMHD = isMHD)
    end

    #=
    ng = 1 - first(eachindex(KS.pSpace.x))
    if bc == :extra
        for i in 1:ng
            ctr[1-i].w .= ctr[1].w
            ctr[1-i].prim .= ctr[1].prim
            ctr[1-i].h0 .= ctr[1].h0
            ctr[1-i].h1 .= ctr[1].h1
            ctr[1-i].h2 .= ctr[1].h2
            ctr[1-i].E .= ctr[1].E
            ctr[1-i].B .= ctr[1].B
            ctr[1-i].ϕ = deepcopy(ctr[1].ϕ)
            ctr[1-i].ψ = deepcopy(ctr[1].ψ)
            ctr[1-i].lorenz .= ctr[1].lorenz

            ctr[KS.pSpace.nx+i].w .= ctr[KS.pSpace.nx].w
            ctr[KS.pSpace.nx+i].prim .= ctr[KS.pSpace.nx].prim
            ctr[KS.pSpace.nx+i].h0 .= ctr[KS.pSpace.nx].h0
            ctr[KS.pSpace.nx+i].h1 .= ctr[KS.pSpace.nx].h1
            ctr[KS.pSpace.nx+i].h2 .= ctr[KS.pSpace.nx].h2
            ctr[KS.pSpace.nx+i].E .= ctr[KS.pSpace.nx].E
            ctr[KS.pSpace.nx+i].B .= ctr[KS.pSpace.nx].B
            ctr[KS.pSpace.nx+i].ϕ = deepcopy(ctr[KS.pSpace.nx].ϕ)
            ctr[KS.pSpace.nx+i].ψ = deepcopy(ctr[KS.pSpace.nx].ψ)
            ctr[KS.pSpace.nx+i].lorenz .= ctr[KS.pSpace.nx].lorenz
        end
    elseif bc == :period
        for i in 1:ng
            ctr[1-i].w .= ctr[KS.pSpace.nx+1-i].w
            ctr[1-i].prim .= ctr[KS.pSpace.nx+1-i].prim
            ctr[1-i].h0 .= ctr[KS.pSpace.nx+1-i].h0
            ctr[1-i].h1 .= ctr[KS.pSpace.nx+1-i].h1
            ctr[1-i].h2 .= ctr[KS.pSpace.nx+1-i].h2
            ctr[1-i].E .= ctr[KS.pSpace.nx+1-i].E
            ctr[1-i].B .= ctr[KS.pSpace.nx+1-i].B
            ctr[1-i].ϕ = deepcopy(ctr[KS.pSpace.nx+1-i].ϕ)
            ctr[1-i].ψ = deepcopy(ctr[KS.pSpace.nx+1-i].ψ)
            ctr[1-i].lorenz .= ctr[KS.pSpace.nx+1-i].lorenz

            ctr[KS.pSpace.nx+i].w .= ctr[i].w
            ctr[KS.pSpace.nx+i].prim .= ctr[i].prim
            ctr[KS.pSpace.nx+i].h0 .= ctr[i].h0
            ctr[KS.pSpace.nx+i].h1 .= ctr[i].h1
            ctr[KS.pSpace.nx+i].h2 .= ctr[i].h2
            ctr[KS.pSpace.nx+i].E .= ctr[i].E
            ctr[KS.pSpace.nx+i].B .= ctr[i].B
            ctr[KS.pSpace.nx+i].ϕ = deepcopy(ctr[i].ϕ)
            ctr[KS.pSpace.nx+i].ψ = deepcopy(ctr[i].ψ)
            ctr[KS.pSpace.nx+i].lorenz .= ctr[i].lorenz
        end
    elseif bc == :balance
        @. ctr[0].w = 0.5 * (ctr[-1].w + ctr[1].w)
        @. ctr[0].prim = 0.5 * (ctr[-1].prim + ctr[1].prim)
        @. ctr[0].h0 = 0.5 * (ctr[-1].h0 + ctr[1].h0)
        @. ctr[0].h1 = 0.5 * (ctr[-1].h1 + ctr[1].h1)
        @. ctr[0].h2 = 0.5 * (ctr[-1].h2 + ctr[1].h2)
        @. ctr[0].E = 0.5 * (ctr[-1].E + ctr[1].E)
        @. ctr[0].B = 0.5 * (ctr[-1].B + ctr[1].B)
        ctr[0].ϕ = 0.5 * (ctr[-1].ϕ + ctr[1].ϕ)
        ctr[0].ψ = 0.5 * (ctr[-1].ψ + ctr[1].ψ)
        @. ctr[0].lorenz = 0.5 * (ctr[-1].lorenz + ctr[1].lorenz)

        @. ctr[KS.pSpace.nx+1].w = 0.5 * (ctr[KS.pSpace.nx].w + ctr[KS.pSpace.nx+2].w)
        @. ctr[KS.pSpace.nx+1].prim = 0.5 * (ctr[KS.pSpace.nx].prim + ctr[KS.pSpace.nx+2].prim)
        @. ctr[KS.pSpace.nx+1].h0 = 0.5 * (ctr[KS.pSpace.nx].h0 + ctr[KS.pSpace.nx+2].h0)
        @. ctr[KS.pSpace.nx+1].h1 = 0.5 * (ctr[KS.pSpace.nx].h1 + ctr[KS.pSpace.nx+2].h1)
        @. ctr[KS.pSpace.nx+1].h2 = 0.5 * (ctr[KS.pSpace.nx].h2 + ctr[KS.pSpace.nx+2].h2)
        @. ctr[KS.pSpace.nx+1].E = 0.5 * (ctr[KS.pSpace.nx].E + ctr[KS.pSpace.nx+2].E)
        @. ctr[KS.pSpace.nx+1].B = 0.5 * (ctr[KS.pSpace.nx].B + ctr[KS.pSpace.nx+2].B)
        ctr[KS.pSpace.nx+1].ϕ = 0.5 * (ctr[KS.pSpace.nx].ϕ + ctr[KS.pSpace.nx+2].ϕ)
        ctr[KS.pSpace.nx+1].ψ = 0.5 * (ctr[KS.pSpace.nx].ψ + ctr[KS.pSpace.nx+2].ψ)
        @. ctr[KS.pSpace.nx+1].lorenz = 0.5 * (ctr[KS.pSpace.nx].lorenz + ctr[KS.pSpace.nx+2].lorenz)
    else
    end
    =#
    return nothing

end

function update!(
    KS::X,
    ctr::Y,
    face::Z,
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
    isMHD = true::Bool,
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolume1D4F,1},
    Z<:AbstractArray{Interface1D4F,1},
}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds Threads.@threads for i = 2:KS.pSpace.nx-1
        step!(KS, face[i], ctr[i], face[i+1], KS.ps.dx[i], dt, sumRes, sumAvg, coll, isMHD)
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * KS.pSpace.nx) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(
            KS,
            ctr,
            face,
            dt,
            residual;
            coll = coll,
            bc = [bc, bc],
            isMHD = isMHD,
        )
    else
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc, isMHD = isMHD)
    end

    return nothing

end

function update!(
    KS::X,
    ctr::Y,
    a1face::Z,
    a2face::Z,
    dt,
    residual;
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolume2D,2},
    Z<:AbstractArray{Interface2D,2},
}

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    if ndims(sumRes) == 1

        @inbounds for j = 2:ny-1
            for i = 2:nx-1
                step!(
                    ctr[i, j].w,
                    ctr[i, j].prim,
                    a1face[i, j].fw,
                    a1face[i+1, j].fw,
                    a2face[i, j].fw,
                    a2face[i, j+1].fw,
                    KS.gas.γ,
                    dx[i, j] * dy[i, j],
                    sumRes,
                    sumAvg,
                    coll,
                )
            end
        end

    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * nx * ny) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(
            KS,
            ctr,
            a1face,
            a2face,
            dt,
            residual;
            coll = coll,
            bc = [bc, bc, bc, bc],
        )
    else
        update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = bc)
    end

    return nothing

end

function update!(
    KS::X,
    ctr::Y,
    a1face::Z,
    a2face::Z,
    dt,
    residual;
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolume2D1F,2},
    Z<:AbstractArray{Interface2D1F,2},
}

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    if ndims(sumRes) == 1

        @inbounds for j = 2:ny-1
            for i = 2:nx-1
                step!(
                    ctr[i, j].w,
                    ctr[i, j].prim,
                    ctr[i, j].f,
                    a1face[i, j].fw,
                    a1face[i, j].ff,
                    a1face[i+1, j].fw,
                    a1face[i+1, j].ff,
                    a2face[i, j].fw,
                    a2face[i, j].ff,
                    a2face[i, j+1].fw,
                    a2face[i, j+1].ff,
                    KS.vSpace.u,
                    KS.vSpace.v,
                    KS.vSpace.weights,
                    KS.gas.γ,
                    KS.gas.μᵣ,
                    KS.gas.ω,
                    KS.gas.Pr,
                    dx[i, j] * dy[i, j],
                    dt,
                    sumRes,
                    sumAvg,
                    coll,
                )
            end
        end

    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * nx * ny) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(
            KS,
            ctr,
            a1face,
            a2face,
            dt,
            residual;
            coll = coll,
            bc = [bc, bc, bc, bc],
        )
    else
        update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = bc)
    end

    return nothing

end

#--- 2d2f ---#
function update!(
    KS::X,
    ctr::Y,
    a1face::Z,
    a2face::Z,
    dt,
    residual;
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolume2D2F,2},
    Z<:AbstractArray{Interface2D2F,2},
}

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    if ndims(sumRes) == 1

        @inbounds for j = 2:ny-1
            for i = 2:nx-1
                step!(
                    ctr[i, j].w,
                    ctr[i, j].prim,
                    ctr[i, j].h,
                    ctr[i, j].b,
                    a1face[i, j].fw,
                    a1face[i, j].fh,
                    a1face[i, j].fb,
                    a1face[i+1, j].fw,
                    a1face[i+1, j].fh,
                    a1face[i+1, j].fb,
                    a2face[i, j].fw,
                    a2face[i, j].fh,
                    a2face[i, j].fb,
                    a2face[i, j+1].fw,
                    a2face[i, j+1].fh,
                    a2face[i, j+1].fb,
                    KS.vSpace.u,
                    KS.vSpace.v,
                    KS.vSpace.weights,
                    KS.gas.K,
                    KS.gas.γ,
                    KS.gas.μᵣ,
                    KS.gas.ω,
                    KS.gas.Pr,
                    dx[i, j] * dy[i, j],
                    dt,
                    sumRes,
                    sumAvg,
                    coll,
                )
            end
        end

    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * nx * ny) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(
            KS,
            ctr,
            a1face,
            a2face,
            dt,
            residual;
            coll = coll,
            bc = [bc, bc, bc, bc],
        )
    else
        update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = bc)
    end

    return nothing

end

function update!(
    KS::X,
    ctr::Y,
    face::Z,
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolumeUS,1},
    Z<:AbstractArray{Interface2D,1},
}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds Threads.@threads for i in eachindex(ctr)
        if KS.ps.cellType[i] in (0, 2)
            dirc = [sign(dot(ctr[i].n[j], face[KS.ps.cellFaces[i, j]].n)) for j = 1:3]

            step!(
                ctr[i].w,
                ctr[i].prim,
                face[KS.ps.cellFaces[i, 1]].fw,
                face[KS.ps.cellFaces[i, 2]].fw,
                face[KS.ps.cellFaces[i, 3]].fw,
                KS.gas.γ,
                KS.pSpace.cellArea[i],
                dirc,
                sumRes,
                sumAvg,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * size(KS.ps.cellid, 1)) / (sumAvg[i] + 1.e-7)
    end

    update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc)

    return nothing

end

function update!(
    KS::X,
    ctr::Y,
    face::Z,
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolumeUS1F,1},
    Z<:AbstractArray{Interface2D1F,1},
}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds Threads.@threads for i in eachindex(ctr)
        if KS.ps.cellType[i] in (0, 2)
            dirc = [sign(dot(ctr[i].n[j], face[KS.ps.cellFaces[i, j]].n)) for j = 1:3]

            step!(
                ctr[i].w,
                ctr[i].prim,
                ctr[i].f,
                face[KS.ps.cellFaces[i, 1]].fw,
                face[KS.ps.cellFaces[i, 1]].ff,
                face[KS.ps.cellFaces[i, 2]].fw,
                face[KS.ps.cellFaces[i, 2]].ff,
                face[KS.ps.cellFaces[i, 3]].fw,
                face[KS.ps.cellFaces[i, 3]].ff,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.pSpace.cellArea[i],
                dirc,
                dt,
                sumRes,
                sumAvg,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * size(KS.ps.cellid, 1)) / (sumAvg[i] + 1.e-7)
    end

    update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc)

    return nothing

end

function update!(
    KS::X,
    ctr::Y,
    face::Z,
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{ControlVolumeUS2F,1},
    Z<:AbstractArray{Interface2D2F,1},
}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds Threads.@threads for i in eachindex(ctr)
        if KS.ps.cellType[i] in (0, 2)
            dirc = [sign(dot(ctr[i].n[j], face[KS.ps.cellFaces[i, j]].n)) for j = 1:3]

            step!(
                ctr[i].w,
                ctr[i].prim,
                ctr[i].h,
                ctr[i].b,
                face[KS.ps.cellFaces[i, 1]].fw,
                face[KS.ps.cellFaces[i, 1]].fh,
                face[KS.ps.cellFaces[i, 1]].fb,
                face[KS.ps.cellFaces[i, 2]].fw,
                face[KS.ps.cellFaces[i, 2]].fh,
                face[KS.ps.cellFaces[i, 2]].fb,
                face[KS.ps.cellFaces[i, 3]].fw,
                face[KS.ps.cellFaces[i, 3]].fh,
                face[KS.ps.cellFaces[i, 3]].fb,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.pSpace.cellArea[i],
                dirc,
                dt,
                sumRes,
                sumAvg,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * size(KS.ps.cellid, 1)) / (sumAvg[i] + 1.e-7)
    end

    update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc)

    return nothing

end


text:
struct Solver{T} <: MOI.AbstractNLPEvaluator
    nlp::NLPData{T}
    data::SolverData
end

function Solver(dynamics::Vector{Dynamics{T}}, objective::Objective{T}, constraints::Constraints{T}, bounds::Bounds{T};
    evaluate_hessian=false,
    general_constraint=GeneralConstraint(),
    options=Options(),
    parameters=[[zeros(num_parameter) for num_parameter in dimensions(dynamics)[3]]..., zeros(0)]) where T

    trajopt = TrajectoryOptimizationData(objective, dynamics, constraints, bounds,
        parameters=parameters)
    nlp = NLPData(trajopt,
        general_constraint=general_constraint,
        evaluate_hessian=evaluate_hessian)
    data = SolverData(nlp,
        options=options)

    Solver(nlp, data)
end

function initialize_states!(solver::Solver, states)
    for (t, xt) in enumerate(states)
        n = length(xt)
        for i = 1:n
            MOI.set(solver.data.optimizer, MOI.VariablePrimalStart(), solver.data.variables[solver.nlp.indices.states[t][i]], xt[i])
        end
    end
end

function initialize_controls!(solver::Solver, actions)
    for (t, ut) in enumerate(actions)
        m = length(ut)
        for j = 1:m
            MOI.set(solver.data.optimizer, MOI.VariablePrimalStart(), solver.data.variables[solver.nlp.indices.actions[t][j]], ut[j])
        end
    end
end

function get_trajectory(solver::Solver)
    return solver.nlp.trajopt.states, solver.nlp.trajopt.actions[1:end-1]
end

function solve!(solver::Solver)
    MOI.optimize!(solver.data.optimizer)
end


text:
@model function linear_regression(x, y)
    σ ~ Exponential(1)
    α ~ Normal(50, 10)
    βₗ ~ LogNormal(-5, 2)

    N = length(y)
    for n ∈ 1:N
        μ[n] ~ α + βₗ * x
    end

    y ~ Normal(μ, σ)
end


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

## initdefs.jl - initialization and runtime management definitions

"""
    PROGRAM_FILE

A string containing the script name passed to Julia from the command line. Note that the
script name remains unchanged from within included files. Alternatively see
[`@__FILE__`](@ref).
"""
global PROGRAM_FILE = ""

"""
    ARGS

An array of the command line arguments passed to Julia, as strings.
"""
const ARGS = String[]

"""
    exit(code=0)

Stop the program with an exit code. The default exit code is zero, indicating that the
program completed successfully. In an interactive session, `exit()` can be called with
the keyboard shortcut `^D`.
"""
exit(n) = ccall(:jl_exit, Cvoid, (Int32,), n)
exit() = exit(0)

const roottask = current_task()

is_interactive = false

"""
    isinteractive() -> Bool

Determine whether Julia is running an interactive session.
"""
isinteractive() = (is_interactive::Bool)

## package depots (registries, packages, environments) ##

const DEPOT_PATH = String[]

function init_depot_path()
    if haskey(ENV, "JULIA_DEPOT_PATH")
        depots = split(ENV["JULIA_DEPOT_PATH"], Sys.iswindows() ? ';' : ':')
        append!(empty!(DEPOT_PATH), map(expanduser, depots))
    else
        push!(empty!(DEPOT_PATH), joinpath(homedir(), ".julia"))
        push!(DEPOT_PATH, abspath(Sys.BINDIR, "..", "local", "share", "julia"))
        push!(DEPOT_PATH, abspath(Sys.BINDIR, "..", "share", "julia"))
    end
end

## LOAD_PATH, HOME_PROJECT & ACTIVE_PROJECT ##

# JULIA_LOAD_PATH: split on `:` (or `;` on Windows)
# first empty entry is replaced with DEFAULT_LOAD_PATH, the rest are skipped
# entries starting with `@` are named environments:
#  - the first three `#`s in a named environment are replaced with version numbers
#  - `@stdlib` is a special name for the standard library and expands to its path

# if you want a current env setup, use direnv and
# have your .envrc do something like this:
#
#   export JULIA_LOAD_PATH="$(pwd):$JULIA_LOAD_PATH"
#
# this will inherit an existing JULIA_LOAD_PATH value or if there is none, leave
# a trailing empty entry in JULIA_LOAD_PATH which will be replaced with defaults.

const DEFAULT_LOAD_PATH = ["@", "@v#.#", "@stdlib"]

"""
    LOAD_PATH

An array of paths for `using` and `import` statements to consider as project
environments or package directories when loading code. It is populated based on
the [`JULIA_LOAD_PATH`](@ref JULIA_LOAD_PATH) environment variable if set;
otherwise it defaults to `["@", "@v#.#", "@stdlib"]`. Entries starting with `@`
have special meanings:

- `@` refers to the "current active environment", the initial value of which is
  initially determined by the [`JULIA_PROJECT`](@ref JULIA_PROJECT) environment
  variable or the `--project` command-line option.

- `@stdlib` expands to the absolute path of the current Julia installation's
  standard library directory.

- `@name` refers to a named environment, which are stored in depots (see
  [`JULIA_DEPOT_PATH`](@ref JULIA_DEPOT_PATH)) under the `environments`
  subdirectory. The user's named environments are stored in
  `~/.julia/environments` so `@name` would refer to the environment in
  `~/.julia/environments/name` if it exists and contains a `Project.toml` file.
  If `name` contains `#` characters, then they are replaced with the major, minor
  and patch components of the Julia version number. For example, if you are
  running Julia 1.2 then `@v#.#` expands to `@v1.2` and will look for an
  environment by that name, typically at `~/.julia/environments/v1.2`.

The fully expanded value of `LOAD_PATH` that is searched for projects and packages
can be seen by calling the `Base.load_path()` function.

See also:
[`JULIA_LOAD_PATH`](@ref JULIA_LOAD_PATH),
[`JULIA_PROJECT`](@ref JULIA_PROJECT),
[`JULIA_DEPOT_PATH`](@ref JULIA_DEPOT_PATH), and
[Code Loading](@ref Code-Loading).
"""
const LOAD_PATH = copy(DEFAULT_LOAD_PATH)
const HOME_PROJECT = Ref{Union{String,Nothing}}(nothing)
const ACTIVE_PROJECT = Ref{Union{String,Nothing}}(nothing)

function current_project(dir::AbstractString)
    # look for project file in current dir and parents
    home = homedir()
    while true
        for proj in project_names
            file = joinpath(dir, proj)
            isfile_casesensitive(file) && return file
        end
        # bail at home directory
        dir == home && break
        old, dir = dir, dirname(dir)
        dir == old && break
    end
end

function current_project()
    dir = try pwd()
    catch err
        err isa IOError || rethrow()
        return nothing
    end
    return current_project(dir)
end

function parse_load_path(str::String)
    envs = String[]
    isempty(str) && return envs
    for env in split(str, Sys.iswindows() ? ';' : ':')
        if isempty(env)
            for env′ in DEFAULT_LOAD_PATH
                env′ in envs || push!(envs, env′)
            end
        else
            if env == "@."
                env = current_project()
                env === nothing && continue
            end
            env in envs || push!(envs, env)
        end
    end
    return envs
end

function init_load_path()
    if Base.creating_sysimg
        paths = ["@stdlib"]
    elseif haskey(ENV, "JULIA_LOAD_PATH")
        paths = parse_load_path(ENV["JULIA_LOAD_PATH"])
    else
        paths = filter!(env -> env !== nothing,
            [env == "@." ? current_project() : env for env in DEFAULT_LOAD_PATH])
    end
    project = (JLOptions().project != C_NULL ?
        unsafe_string(Base.JLOptions().project) :
        get(ENV, "JULIA_PROJECT", nothing))
    HOME_PROJECT[] =
        project == nothing ? nothing :
        project == "" ? nothing :
        project == "@." ? current_project() : abspath(project)
    append!(empty!(LOAD_PATH), paths)
end

## load path expansion: turn LOAD_PATH entries into concrete paths ##

function load_path_expand(env::AbstractString)::Union{String, Nothing}
    # named environment?
    if startswith(env, '@')
        # `@` in JULIA_LOAD_PATH is expanded early (at startup time)
        # if you put a `@` in LOAD_PATH manually, it's expanded late
        env == "@" && return active_project(false)
        env == "@." && return current_project()
        env == "@stdlib" && return Sys.STDLIB
        env = replace(env, '#' => VERSION.major, count=1)
        env = replace(env, '#' => VERSION.minor, count=1)
        env = replace(env, '#' => VERSION.patch, count=1)
        name = env[2:end]
        # look for named env in each depot
        for depot in DEPOT_PATH
            path = joinpath(depot, "environments", name)
            isdir(path) || continue
            for proj in project_names
                file = abspath(path, proj)
                isfile_casesensitive(file) && return file
            end
        end
        isempty(DEPOT_PATH) && return nothing
        return abspath(DEPOT_PATH[1], "environments", name, project_names[end])
    end
    # otherwise, it's a path
    path = abspath(env)
    if isdir(path)
        # directory with a project file?
        for proj in project_names
            file = joinpath(path, proj)
            isfile_casesensitive(file) && return file
        end
    end
    # package dir or path to project file
    return path
end
load_path_expand(::Nothing) = nothing

function active_project(search_load_path::Bool=true)
    for project in (ACTIVE_PROJECT[], HOME_PROJECT[])
        project == "@" && continue
        project = load_path_expand(project)
        project === nothing && continue
        if !isfile_casesensitive(project) && basename(project) ∉ project_names
            project = abspath(project, "Project.toml")
        end
        return project
    end
    search_load_path || return
    for project in LOAD_PATH
        project == "@" && continue
        project = load_path_expand(project)
        project === nothing && continue
        isfile_casesensitive(project) && return project
        ispath(project) && continue
        basename(project) in project_names && return project
    end
end

function load_path()
    paths = String[]
    for env in LOAD_PATH
        path = load_path_expand(env)
        path !== nothing && path ∉ paths && push!(paths, path)
    end
    return paths
end

## atexit: register exit hooks ##

const atexit_hooks = []

"""
    atexit(f)

Register a zero-argument function `f()` to be called at process exit. `atexit()` hooks are
called in last in first out (LIFO) order and run before object finalizers.
"""
atexit(f::Function) = (pushfirst!(atexit_hooks, f); nothing)

function _atexit()
    for f in atexit_hooks
        try
            f()
        catch err
            show(stderr, err)
            println(stderr)
        end
    end
end

## hook for disabling threaded libraries ##

library_threading_enabled = true
const disable_library_threading_hooks = []

function at_disable_library_threading(f)
    push!(disable_library_threading_hooks, f)
    if !library_threading_enabled
        disable_library_threading()
    end
    return
end

function disable_library_threading()
    global library_threading_enabled = false
    while !isempty(disable_library_threading_hooks)
        f = pop!(disable_library_threading_hooks)
        try
            f()
        catch err
            @warn("a hook from a library to disable threading failed:",
                  exception = (err, catch_backtrace()))
        end
    end
    return
end

max_stars_repo_name: briochemc/julia

max_issues_repo_name: briochemc/julia

max_forks_repo_name: briochemc/julia


text:
module QViM

include("Utils/Utils.jl")
include("QuantInstr.jl")
include("VirtualMachine.jl")

using .VirtualMachine
export QVM, showstate, step!, execute!

using .QuantInstr
export @quip, QInstr

end


text:
"""
$(SIGNATURES)

Convenience function to call the base markdown to html converter on "simple" strings (i.e. strings
that don't need to be further considered and don't contain anything else than markdown tokens).
The boolean `stripp` indicates whether to remove the inserted `<p>` and `</p>` by the base markdown
processor, this is relevant for things that are parsed within latex commands etc.
"""
function md2html(ss::AS; stripp::Bool=false)::AS
    # if there's nothing, return that...
    isempty(ss) && return ss
    # Use Julia's Markdown parser followed by Julia's MD->HTML conversion
    partial = ss |> fix_inserts |> Markdown.parse |> Markdown.html
    # In some cases, base converter adds
...
\n which we might not want
    stripp || return partial
    startswith(partial, "
")    && (partial = chop(partial, head=3))
    endswith(partial,   "
")   && return chop(partial, tail=4)
    endswith(partial,   "
\n") && return chop(partial, tail=5)
    return partial
end

"""
$(SIGNATURES)

Convenience function to check if `idx` is smaller than the length of `v`, if it is, then return the starting point of `v[idx]` (via `from`), otherwise return `BIG_INT`.
"""
from_ifsmaller(v::Vector, idx::Int, len::Int)::Int = (idx > len) ? BIG_INT : from(v[idx])

"""
$(SIGNATURES)

Since divs are recursively processed, once they've been found, everything inside them needs to be
deactivated and left for further re-processing to avoid double inclusion.
"""
function deactivate_divs(blocks::Vector{OCBlock})::Vector{OCBlock}
    active_blocks = ones(Bool, length(blocks))
    for (i, β) ∈ enumerate(blocks)
        fromβ, toβ = from(β), to(β)
        active_blocks[i] || continue
        if β.name == :DIV
            innerblocks = findall(b -> (fromβ < from(b) < toβ), blocks)
            active_blocks[innerblocks] .= false
        end
    end
    return blocks[active_blocks]
end

"""
$(SIGNATURES)

The insertion token have whitespaces around them: ` ##JDINSERT## `, this mostly helps but causes
a problem when combined with italic or bold markdown mode since `_blah_` works but not `_ blah _`.
This function looks for any occurrence of `[\\*_] ##JDINSERT##` or the opposite and removes the
extraneous whitespace.
"""
fix_inserts(s::AS)::String =
    replace(replace(s, r"([\*_]) ##JDINSERT##" => s"\1##JDINSERT##"),
                       r"##JDINSERT## ([\*_])" => s"##JDINSERT##\1")


text:
# Various helper functions to calculate dimensions for operations
include("dim_helpers/ConvDims.jl")
include("dim_helpers/DenseConvDims.jl")
include("dim_helpers/DepthwiseConvDims.jl")
include("dim_helpers/PoolDims.jl")

"""
    transpose_swapbatch(x::AbstractArray)

Given an AbstractArray, swap its batch and channel axes, as we must during transposed
convolution.  We do this to the operands during convolution, and then again to the
output once we're done.
"""
function transpose_swapbatch(x::AbstractArray)
    return permutedims(x, ((1:(ndims(x)-2))..., ndims(x), ndims(x)-1))
end
function transpose_swapbatch(x::Tuple)
    return (x[1:end-2]..., x[end], x[end-1])
end

"""
    transpose_pad(cdims::ConvDims)

Transposed convolution can be calculated in terms of typical convolution with some extra
padding.  This method computes the padding of the convolution that would result in the
transposed convolution of two operands, in essence taking care of that "extra padding".
Note that this method should almost always be accompanied by a call that predilates one
of the operands.
"""
function transpose_pad(cdims::ConvDims)
    I = input_size(cdims)
    K = kernel_size(cdims)
    D = dilation(cdims)
    P = padding(cdims)
    S = stride(cdims)
    return ntuple(length(P)) do i
        hi = ceil(Int, i/2)
        if mod(i, 2) == 1
            return (K[hi] - 1)*D[hi] - P[i]
        else
            return (K[hi] - 1)*D[hi] - P[i] + mod(I[hi] + P[i-1] + P[i] - (K[hi] - 1)*D[hi] - 1, S[hi])
        end
    end
end

"""
    insert_singleton_spatial_dimension(cdims::DenseConvDims)

When converting a 1d convolution to a 2d, or a 2d to a 3d, we need to insert a singleton
spatial dimension at the end of the spatial dimensions.  This does so for a ConvDims.
"""
@inline function insert_singleton_spatial_dimension(cdims::C) where {C <: ConvDims}
    return basetype(C)(cdims;
        N=spatial_dims(cdims) + 1,
        I=(input_size(cdims)..., 1),
        K=(kernel_size(cdims)..., 1),
        S=(stride(cdims)..., 1),
        # Padding is always the problem child....
        P=(padding(cdims)..., 0, 0),
        D=(dilation(cdims)..., 1),
    )
end

# We specialize common cases
@inline function insert_singleton_spatial_dimension(x::AbstractArray{T,3}) where {T}
    return reshape(x, size(x,1), 1, size(x,2), size(x,3))
end
@inline function insert_singleton_spatial_dimension(x::AbstractArray{T,4}) where {T}
    return reshape(x, size(x,1), size(x,2), 1, size(x,3), size(x,4))
end

# Helper to do this as many times as needed
@inline function insert_singleton_spatial_dimension(x, reps::Int)
    for r in 1:reps
        x = insert_singleton_spatial_dimension(x)
    end
    return x
end

"""
    predilated_size(x_size::Tuple, dilation::Tuple)

Calculate the size of a predilated `x` given a particular dilation factor.  This is used
within `predilate()` and `transpose_cdims()`.
"""
function predilated_size(x_size::NTuple{N}, dilation::NTuple{M}) where {N, M}
    @assert (M == N - 2) DimensionMismatch("len(dilation) != number of spatial dims")
    return ntuple(N) do idx
        if idx <= N - 2
            return (x_size[idx] - 1)*dilation[idx] + 1
        else
            x_size[idx]
        end
    end
end

"""
    predilate(x, dilation::Tuple)

Places elements of `x` within a lattice of zeros, used in expressing a transposed
convolution in terms of normal convolution.  Note that while we call this "predilation"
for aesthetic reasons, you are typically passing a "stride" value into here.  Yes,
transposed convolution is confusing.
"""
function predilate(x::AbstractArray{T,N}, dilation::NTuple{M}) where {T, N, M}
    @assert (M == N - 2) DimensionMismatch("len(dilation) != number of spatial dims")

    # If there is no dilation to be done, then ignore it.
    if all(dilation .== 1)
        return x
    end

    # Validate dilation factors
    for idx in 1:length(dilation)
        @assert dilation[idx] >= 1 ArgumentError("dilation cannot be less than 1")
    end

    # Create new x that is bigger and holier
    x_dil = zeros(eltype(x), predilated_size(size(x), dilation))

    # Fill in strategic locations within `x_dil`, such that there are `dilation[idx] - 1`
    # zeros between each element of `x` along each spatial dimension.
    x_dil[(1:dilation[idx]:size(x_dil,idx) for idx in 1:(N-2))..., :, :] .= x
    return x_dil
end

"""
    flipweight(w::AbstractArray)

Reorders the weight tensor for supporting both convolution and cross-correlation operations.
"""

# For any array with ndims <= 3 it makes no sense to flip the weights so simply return the
# original array
@inline flipweight(w::AbstractArray) = w

@inline flipweight(w::AbstractArray{T, 4}) where {T} = w[end:-1:1, end:-1:1, :, :]

@inline flipweight(w::AbstractArray{T, 5}) where {T} = w[end:-1:1, end:-1:1, end:-1:1, :, :]