text:
import Base: iterate

mutable struct SubBagIter
    bag::Bag
    k::Int
    child::Union{Nothing,SubBagIter}
end

function SubBagIter(bag, k)
   return SubBagIter(bag, k, nothing)
end

########################
# HELPER FUNCTIONS
########################
function create_child_bag(bag, k, idx)
    child = Int[min(k, count) for count in bag[idx:end]]
    child[1] -= 1
    return child
end

function create_item(bag, idx)
    item = zero(bag)
    item[idx] += 1
    return item
end

function combine_items(bag, idx, child_item)
    item = create_item(bag, idx)
    item[idx:end] .+= child_item
    return item
end

function combine_states(idx, child_state)
    return [[idx]; child_state]
end

###################################
# ITERATOR INTERFACE DEFINITION
###################################
function iterate(sbi::SubBagIter)

    # No such subbags exist -- no items
    # to iterate through!
    if sum(sbi.bag) < sbi.k
        return nothing
    end

    # Otherwise, advance to the first nonzero entry
    # (which is guaranteed to exist)
    idx = 1
    while sbi.bag[idx] <= 0
        idx += 1
    end

    # Base case: k == 1
    if sbi.k <= 1
        state = [idx]
        item = create_item(sbi.bag, idx)
        return item, state
    end

    # Recursive case: k > 1

    # Create the "child" bag of remaining items and
    # construct a child iterator to advance through it.
    child_bag = create_child_bag(sbi.bag, sbi.k, idx)
    sbi.child = SubBagIter(child_bag, sbi.k-1)

    # Call "iterate" on the child in order to initialize it.
    # Guaranteed to return (item, state) tuple, rather than nothing.
    child_item, child_state = iterate(sbi.child)

    # combine the child item with this iterator's item
    item = combine_items(sbi.bag, idx, child_item)
    state = combine_states(idx, child_state)

    return item, state
end

function iterate(sbi::SubBagIter, state)

    # Base case: k = 1
    if sbi.k == 1

        # Advance to the next nonzero entry
        idx = state[1]
        for (i, count) in enumerate(sbi.bag[idx+1:end])
            if count > 0
                idx += i
                break
            end
        end
        # If the idx is unchanged, then we've reached the end! return `nothing`
        if idx == state[1]
             return nothing
        else # Otherwise, return the current item and state
            item = create_item(sbi.bag, idx)
            state = [idx]
            return item, state
        end
    end

    # Recursive case: k > 1

    # Iterate the child
    child_result = iterate(sbi.child, state[2:end])

    # If the child returned (item, state),
    # then we return combined items/states
    if child_result != nothing
        child_item = child_result[1]
        child_state = child_result[2]

        item = combine_items(sbi.bag, state[1], child_item)
        state = combine_states(state[1], child_state)
        return item, state
    else
    # If the child returned `nothing`, then
    # we advance idx to the next nonzero entry
        idx = state[1]
        for (i, count) in enumerate(sbi.bag[idx+1:end])
            if count > 0
                idx += i
                break
            end
        end

        # If the idx was not updated, then we've run out of
        # nonzero entries -- we're done!
        if idx == state[1]
            return nothing
        else
        # Otherwise, re-initialize the child at this new location
            sbi.child.bag = create_child_bag(sbi.bag, sbi.k, idx)
            child_result = iterate(sbi.child)

            # If the re-initialized child returns nothing,
            # then we're done!
            if child_result == nothing
                return nothing
            end

            # Otherwise, we just return the next item/state
            child_item = child_result[1]
            child_state = child_result[2]
            item = combine_items(sbi.bag, idx, child_item)
            state = combine_states(idx, child_state)
            return item, state
        end
    end
end

#########################################
# Iterate through ALL subbags
#########################################

mutable struct AllSubBagIter
    bag::Bag
    subbag_iter::SubBagIter
end

function AllSubBagIter(bag::Bag)
    subbag_iter = SubBagIter(bag, 1)
    return AllSubBagIter(bag,subbag_iter)
end

function iterate(asbi::AllSubBagIter)

    # Construct subbag iterator; initialize it
    asbi.subbag_iter = SubBagIter(asbi.bag, 1)
    sbi_iter_result = iterate(asbi.subbag_iter)

    # If it returns nothing, then we're done!
    if sbi_iter_result == nothing
        return nothing
    end

    # Otherwise: return the split and the state
    sbi_item = sbi_iter_result[1]
    sbi_state = sbi_iter_result[2]

    sbi_complement = bag_complement(asbi.bag, sbi_item)
    return sbi_item, sbi_state
end

function iterate(asbi::AllSubBagIter, state)

    # Try to advance the subbag_iter
    sbi_iter_result = iterate(asbi.subbag_iter, state)

    # If we get "nothing",
    if sbi_iter_result == nothing
        # then increment k
        asbi.subbag_iter.k = asbi.subbag_iter.k + 1

        # then reset and initialize the subbag_iter
        sbi_iter_result = iterate(asbi.subbag_iter)

        if sbi_iter_result == nothing
            return nothing
        end
    end

    sbi_item = sbi_iter_result[1]
    sbi_state = sbi_iter_result[2]

    # return the split and the state (k, subbag_iter state)
    return sbi_item, sbi_state
end

#########################################
# Iterate through all possible "splits"
# of a bag
#########################################

mutable struct BagSplitIter
    bag::Bag
    N::BagType
    subbag_iter::SubBagIter
end

function BagSplitIter(bag)
    subbag_iter = SubBagIter(bag, 1)
    N = sum(bag)
    return BagSplitIter(bag, N, subbag_iter)
end

function count_subbags(bag, k)
    N = sum(bag)
    redundant = binomial(N,k)
    non_redundant = div(redundant, prod([factorial(b) for b in bag]))
end

function iterate(bsi::BagSplitIter)

    # Initialize the subbag_iter
    subbag_iter_result = iterate(bsi.subbag_iter)

    if subbag_iter_result == nothing
        return nothing
    end

    subbag_item = subbag_iter_result[1]
    subbag_state = subbag_iter_result[2]

    complement = bag_complement(bsi.bag, subbag_item)
    return (subbag_item, complement), ((false,-1,-1), subbag_state)
end

function iterate(bsi::BagSplitIter, state)

    subbag_iter_state = state[2]
    middle_state = state[1]

    # Try iterating the subbag_iter
    subbag_iter_result = iterate(bsi.subbag_iter, subbag_iter_state)

    new_middle_state = (false,-1,-1)

    # If it returned nothing:
    if subbag_iter_result == nothing
        # Increment the iterator's `k`
        bsi.subbag_iter.k = bsi.subbag_iter.k + 1

        N_d_2 = div(bsi.N, 2)

        # Check k's size
        if bsi.subbag_iter.k >= N_d_2

            # If k is too large, we're done!
            if bsi.subbag_iter.k > N_d_2
                return nothing
            # Otherwise, we've reached the "middle size".
            # We need to start tracking our state in the
            # middle size, to avoid redundancies.
            else
                if N_d_2 % 2 == 0
                    max_iter = div(count_subbags(bsi.bag, bsi.subbag_iter.k), 2)
                    new_middle_state = (true, 1, max_iter)
                end
            end
        end
        # Re-initialize the iterator
        subbag_iter_result = iterate(bsi.subbag_iter)
    end
    # If it returned (item, state)
    sub_item = subbag_iter_result[1]
    sub_state = subbag_iter_result[2]
    complement = bag_complement(bsi.bag, sub_item)

    # if we're in the "middle size", then
    # we need to update some additional state
    if middle_state[1]
        new_middle_state = (true, middle_state[2]+1, middle_state[3])
    end

    return (sub_item, complement), (new_middle_state, sub_state)
end

######################################
# Iterate through the decompositions
# encoded in a table
######################################

mutable struct DecompIter
    cache::Dict{Bag,Vector{Bag}}
    bag::Bag
    threshold::Bag
    sub_iter::Union{Nothing,DecompIter}
    depth::Int
end

function DecompIter(cache, bag::Bag, threshold::Bag)
    return DecompIter(cache, bag, threshold, nothing, 0)
end

function DecompIter(cache, bag::Bag, threshold::Bag, depth::Int)
    return DecompIter(cache, bag, threshold, nothing, depth)
end

function DecompIter(cache, bag)
    return DecompIter(cache, bag, zeros(BagType, N_CHAR))
end

function iterate(dci::DecompIter)

    # Base case: bag is empty
    if all(dci.bag .== 0)
        return Bag[], nothing
    end

    # Advance the child_idx to the first child with a valid decomposition.
    child_vec = dci.cache[dci.bag]
    child_idx = 1
    while child_idx <= length(child_vec)

        child = child_vec[child_idx]

        # The child must be >= the threshold, though
        if bag_leq(dci.threshold, child)
            complement = bag_complement(dci.bag, child)

            dci.sub_iter = DecompIter(dci.cache, complement, child, dci.depth+1)
            sub_result = iterate(dci.sub_iter)

            if sub_result != nothing
                sub_item = sub_result[1]
                sub_state = sub_result[2]
                return [Bag[child]; sub_item], (child_idx, sub_state)
            end
        end

        # Otherwise, we move on to the next child
        child_idx += 1
    end

    return nothing
end

function iterate(dci::DecompIter, state)

    # Base case: bag is empty
    if all(dci.bag .== 0)
        return nothing
    end

    # Otherwise: search for the next valid decomposition
    child_idx = state[1]
    sub_state = state[2]
    child_vec = dci.cache[dci.bag]

    # Try to iterate the sub_iter
    sub_result = iterate(dci.sub_iter, sub_state)
    if sub_result != nothing
        sub_item = sub_result[1]
        sub_state = sub_result[2]
        return [Bag[child_vec[child_idx]]; sub_item], (child_idx, sub_state)
    end

    # If the sub_iter terminates, advance the child index
    while true
        while true
            child_idx += 1
            if child_idx > length(child_vec)
                return nothing
            end
            if bag_leq(dci.threshold, child_vec[child_idx])
                break
            end
        end

        # ...and update the sub_iter
        complement = bag_complement(dci.bag, child_vec[child_idx])
        dci.sub_iter = DecompIter(dci.cache, complement, child_vec[child_idx], dci.depth+1)
        sub_result = iterate(dci.sub_iter)

        if sub_result != nothing
            sub_item = sub_result[1]
            sub_state = sub_result[2]
            return [Bag[child_vec[child_idx]]; sub_item], (child_idx, sub_state)
        end
    end

end

#################################
# Cartesian Product Iterator
#################################

mutable struct ProductIter
    collections::AbstractVector
    child::Union{Nothing,ProductIter}
end

function ProductIter(collections::AbstractVector)
    return ProductIter(collections, nothing)
end

function iterate(prit::ProductIter)

    # Base cases
    if length(prit.collections[1]) == 0
        return nothing
    end
    if length(prit.collections) == 1
        return [prit.collections[1][1]], (1, nothing)
    end

    # Recursive case
    prit.child = ProductIter(prit.collections[2:end])
    child_result = iterate(prit.child)
    if child_result == nothing
        return nothing
    else
        child_item = child_result[1]
        child_state = child_result[2]
        return [[prit.collections[1][1]] ; child_item], (1, child_state)
    end
end

function iterate(prit::ProductIter, state)

    cur_idx = state[1]
    child_state = state[2]

    # Base case
    if length(prit.collections) == 1
        cur_idx += 1
        if cur_idx > length(prit.collections[1])
            return nothing
        else
            return [prit.collections[1][cur_idx]], (cur_idx, nothing)
        end
    end

    # Recursive case
    child_result = iterate(prit.child, child_state)
    if child_result != nothing
        child_item = child_result[1]
        child_state = child_result[2]
    else
        cur_idx += 1
        if cur_idx > length(prit.collections[1])
            return nothing
        end
        prit.child = ProductIter(prit.collections[2:end])
        child_item, child_state = iterate(prit.child)
    end
    return [prit.collections[1][cur_idx]; child_item], (cur_idx, child_state)
end

##########################################
# Anagram Iterator
##########################################

mutable struct AnagramIter
    bag::Bag
    decomp_cache::Dict{Bag,Vector{Bag}}
    decomp_iter::Union{Nothing,DecompIter}
    prod_iter::Union{Nothing,ProductIter}
end

function AnagramIter(bag::Bag, decomp_cache)
    return AnagramIter(bag, decomp_cache, nothing, nothing)
end

function iterate(ani::AnagramIter)
    ani.decomp_iter = DecompIter(ani.decomp_cache, ani.bag)
    decomp_result = iterate(ani.decomp_iter)
    if decomp_result == nothing
        return nothing
    end
    decomp_item = decomp_result[1]
    decomp_state = decomp_result[2]

    ani.prod_iter = ProductIter([VOCAB_DICT[bag] for bag in decomp_item])
    prod_result = iterate(ani.prod_iter)
    if prod_result == nothing
        return nothing
    end
    prod_item = prod_result[1]
    prod_state = prod_result[2]

    return prod_item, (decomp_state, prod_state)
end

function iterate(ani::AnagramIter, state)

    decomp_state = state[1]
    prod_state = state[2]

    prod_result = iterate(ani.prod_iter, prod_state)
    if prod_result == nothing
        decomp_result = iterate(ani.decomp_iter, decomp_state)
        if decomp_result == nothing
            return nothing
        else
            decomp_item = decomp_result[1]
            decomp_state = decomp_result[2]
            ani.prod_iter = ProductIter([VOCAB_DICT[bag] for bag in decomp_item])

            prod_result = iterate(ani.prod_iter)
            if prod_result == nothing
                return nothing
            end
            prod_item = prod_result[1]
            prod_state = prod_result[2]
            return prod_item, (decomp_state, prod_state)
        end
    else
        prod_item = prod_result[1]
        prod_state = prod_result[2]
        return prod_item, (decomp_state, prod_state)
    end
end

##########################################
# Subword Iterator
##########################################

mutable struct SubwordIter
    bag::Bag
    subbag_iter::Union{Nothing,AllSubBagIter}
end

function SubwordIter(bag::Bag)
    return SubwordIter(bag, nothing)
end

function iterate(swi::SubwordIter)

    swi.subbag_iter = AllSubBagIter(swi.bag)
    sbi_result = iterate(swi.subbag_iter)
    if sbi_result == nothing
        return nothing
    end

    while true
        if sbi_result == nothing
            return nothing
        end

        sbi_item = sbi_result[1]
        sbi_state = sbi_result[2]

        if haskey(VOCAB_DICT, sbi_item)
            word_idx = 1
            return VOCAB_DICT[sbi_item][word_idx], (sbi_state, sbi_item, word_idx)
        else
            sbi_result = iterate(swi.subbag_iter, sbi_state)
        end
    end
end

function iterate(swi::SubwordIter, state)

    sbi_state = state[1]
    cur_bag = state[2]
    word_idx = state[3]

    word_idx += 1
    if word_idx > length(VOCAB_DICT[cur_bag])

        while true
            sbi_result = iterate(swi.subbag_iter, sbi_state)

            if sbi_result == nothing
                return nothing
            else
                sbi_item = sbi_result[1]
                sbi_state = sbi_result[2]

                if haskey(VOCAB_DICT, sbi_item)
                    word_idx = 1
                    return VOCAB_DICT[sbi_item][word_idx], (sbi_state, sbi_item, word_idx)
                else
                    continue
                end
            end
        end

    else
        return VOCAB_DICT[cur_bag][word_idx], (sbi_state, cur_bag, word_idx)
    end

end


text:
import Base.isless
import Base.isequal
import Base.colon
import Base.endof

typealias SDMContainer Union{SortedDict, SortedMultiDict}
typealias SAContainer Union{SDMContainer, SortedSet}

typealias Token Tuple{SAContainer, IntSemiToken}
typealias SDMToken Tuple{SDMContainer, IntSemiToken}
typealias SetToken Tuple{SortedSet, IntSemiToken}

## Function startof returns the semitoken that points
## to the first sorted order of the tree.  It returns
## the past-end token if the tree is empty.

@inline startof(m::SAContainer) = IntSemiToken(beginloc(m.bt))

## Function endof returns the semitoken that points
## to the last item in the sorted order,
## or the before-start marker if the tree is empty.

@inline endof(m::SAContainer) = IntSemiToken(endloc(m.bt))

## Function pastendsemitoken returns the token past the end of the data.

@inline pastendsemitoken(::SAContainer) = IntSemiToken(2)

## Function beforestarttoken returns the token before the start of the data.

@inline beforestartsemitoken(::SAContainer) = IntSemiToken(1)

## delete! deletes an item given a token.

@inline function delete!(ii::Token)
    has_data(ii)
    delete!(ii[1].bt, ii[2].address)
end

## Function advances takes a token and returns the
## next token in the sorted order.

@inline function advance(ii::Token)
    not_pastend(ii)
    IntSemiToken(nextloc0(ii[1].bt, ii[2].address))
end

## Function regresss takes a token and returns the
## previous token in the sorted order.

@inline function regress(ii::Token)
    not_beforestart(ii)
    IntSemiToken(prevloc0(ii[1].bt, ii[2].address))
end

## status of a token is 0 if the token is invalid, 1 if it points to
## ordinary data, 2 if it points to the before-start location and 3 if
## it points to the past-end location.

@inline status(ii::Token) =
       !(ii[2].address in ii[1].bt.useddatacells)? 0 :
         ii[2].address == 1?                       2 :
         ii[2].address == 2?                       3 : 1

"""
    compare(m::SAContainer, s::IntSemiToken, t::IntSemiToken)

Determines the  relative positions of the  data items indexed
by `(m,s)` and  `(m,t)` in the sorted order. The  return value is `-1`
if `(m,s)` precedes `(m,t)`, `0` if they are equal, and `1` if `(m,s)`
succeeds `(m,t)`. `s`  and `t`  are semitokens  for the  same container `m`.
"""
@inline compare(m::SAContainer, s::IntSemiToken, t::IntSemiToken) =
      compareInd(m.bt, s.address, t.address)

@inline function deref(ii::SDMToken)
    has_data(ii)
    return Pair(ii[1].bt.data[ii[2].address].k, ii[1].bt.data[ii[2].address].d)
end

@inline function deref(ii::SetToken)
    has_data(ii)
    return ii[1].bt.data[ii[2].address].k
end

@inline function deref_key(ii::SDMToken)
    has_data(ii)
    return ii[1].bt.data[ii[2].address].k
end

@inline function deref_value(ii::SDMToken)
    has_data(ii)
    return ii[1].bt.data[ii[2].address].d
end

## Functions setindex! and getindex for semitokens.
## Note that we can't use SDMContainer here; we have
## to spell it out otherwise there is an ambiguity.

@inline function getindex(m::SortedDict,
                          i::IntSemiToken)
    has_data((m,i))
    return m.bt.data[i.address].d
end

@inline function getindex(m::SortedMultiDict,
                          i::IntSemiToken)
    has_data((m,i))
    return m.bt.data[i.address].d
end

@inline function setindex!(m::SortedDict,
                           d_,
                           i::IntSemiToken)
    has_data((m,i))
    m.bt.data[i.address] = KDRec{keytype(m),valtype(m)}(m.bt.data[i.address].parent,
                                                         m.bt.data[i.address].k,
                                                         convert(valtype(m),d_))
    m
end

@inline function setindex!(m::SortedMultiDict,
                           d_,
                           i::IntSemiToken)
    has_data((m,i))
    m.bt.data[i.address] = KDRec{keytype(m),valtype(m)}(m.bt.data[i.address].parent,
                                                         m.bt.data[i.address].k,
                                                         convert(valtype(m),d_))
    m
end

## This function takes a key and returns the token
## of the first item in the tree that is >= the given
## key in the sorted order.  It returns the past-end marker
## if there is none.

@inline function searchsortedfirst(m::SAContainer, k_)
    i = findkeyless(m.bt, convert(keytype(m), k_))
    IntSemiToken(nextloc0(m.bt, i))
end

## This function takes a key and returns a token
## to the first item in the tree that is > the given
## key in the sorted order.  It returns the past-end marker
## if there is none.

@inline function searchsortedafter(m::SAContainer, k_)
    i, exactfound = findkey(m.bt, convert(keytype(m), k_))
    IntSemiToken(nextloc0(m.bt, i))
end

## This function takes a key and returns a token
## to the last item in the tree that is <= the given
## key in the sorted order.  It returns the before-start marker
## if there is none.

@inline function searchsortedlast(m::SAContainer, k_)
    i, exactfound = findkey(m.bt, convert(keytype(m),k_))
    IntSemiToken(i)
end

## The next four are correctness-checking routines.  They are
## not exported.

@inline not_beforestart(i::Token) =
    (!(i[2].address in i[1].bt.useddatacells) ||
     i[2].address == 1) && throw(BoundsError())

@inline not_pastend(i::Token) =
    (!(i[2].address in i[1].bt.useddatacells) ||
     i[2].address == 2) && throw(BoundsError())

@inline has_data(i::Token) =
    (!(i[2].address in i[1].bt.useddatacells) ||
     i[2].address < 3) && throw(BoundsError())

ext: jl


text:
using .Networks

using Dates
using JuMP
using Printf

@with_kw mutable struct EnergyMarket <: AbstractMarket
    configs = EnergyMarketConfigs()
end

function run(runnable::Union{EnergyMarket,BalanceMarket},
            ech::Dates.DateTime, firmness, TS::Vector{Dates.DateTime},
            context::AbstractContext)
    problem_name_l = @sprintf("%s_%s", typeof(runnable), ech)

    gratis_starts = Set{Tuple{String,Dates.DateTime}}()
    if runnable.configs.CONSIDER_GRATIS_STARTS
        gratis_starts = init_gratis_start(context, runnable.configs.REF_SCHEDULE_TYPE)
    end
    @debug("gratis_starts : ", gratis_starts)

    runnable.configs.out_path = context.out_dir
    runnable.configs.problem_name = problem_name_l

    return energy_market(get_network(context),
                        TS,
                        get_generators_initial_state(context),
                        get_scenarios(context, ech),
                        get_uncertainties(context, ech),
                        firmness,
                        get_market_schedule(context),
                        get_tso_schedule(context),
                        get_tso_actions(context),
                        gratis_starts,
                        runnable.configs
                        )
end

function update_market_schedule!(context::AbstractContext, ech,
                                result::EnergyMarketModel,
                                firmness,
                                runnable::Union{EnergyMarket,BalanceMarket})
    market_schedule = get_market_schedule(context)
    market_schedule.decider_type = DeciderType(runnable)
    market_schedule.decision_time = ech

    # Pilotables levels
    for ((gen_id, ts, s), p_injected_var) in result.pilotable_model.p_injected
        if get_power_level_firmness(firmness, gen_id, ts) == FREE
            set_prod_value!(market_schedule, gen_id, ts, s, value(p_injected_var))
        elseif get_power_level_firmness(firmness, gen_id, ts) == TO_DECIDE
            set_prod_definitive_value!(market_schedule, gen_id, ts, value(p_injected_var))
        elseif get_power_level_firmness(firmness, gen_id, ts) == DECIDED
            @assert( ismissing(get_prod_value(market_schedule, gen_id, ts))
                    || (value(p_injected_var) ≈ get_prod_value(market_schedule, gen_id, ts)) )
            set_prod_definitive_value!(market_schedule, gen_id, ts, value(p_injected_var))
        end
    end

    # Commitment
    for ((gen_id, ts, s), b_on_var) in result.pilotable_model.b_on
        gen_state_value = parse(GeneratorState, value(b_on_var))
        if get_commitment_firmness(firmness, gen_id, ts) == FREE
            set_commitment_value!(market_schedule, gen_id, ts, s, gen_state_value)
        elseif get_commitment_firmness(firmness, gen_id, ts) in [TO_DECIDE, DECIDED]
            set_commitment_definitive_value!(market_schedule, gen_id, ts, gen_state_value)
        end
    end

    # Capping
    update_schedule_capping!(market_schedule, context, ech, result.limitable_model, runnable.configs.CONSIDER_TSOACTIONS_LIMITATIONS)

    # Limitables levels : needs to be after capping update
    uncertainties_l = get_uncertainties(context, ech)
    for gen in Networks.get_generators_of_type(get_network(context), Networks.LIMITABLE)
        gen_id = Networks.get_id(gen)
        for ts in get_target_timepoints(context)
            for s in get_scenarios(context)
                capped_l = safeget_capping(market_schedule, gen_id, ts, s)
                injected_l = get_uncertainties(uncertainties_l,gen_id,ts,s) - capped_l
                set_prod_value!(market_schedule, gen_id, ts, s, injected_l)
            end
        end
    end

    # loss_of_load (load-shedding)
    update_schedule_loss_of_load!(market_schedule, context, ech, result.lol_model)

    return market_schedule
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
module TestSpatialDiscretization

using Modia

# Desired:
#   using ModiaMath: plot
#
# In order that these packages need not to be defined in the user environment, they are included via Modia:
using Modia.ModiaMath: plot

const n = 5
@model SpatialDiscretization begin
    u = Var(start=zeros(n))
    c = 1
    @equations begin
        # der(u[2:n]) = ( u[2:n] - u[1:n-1] ) / c
        der(u) = [1; ( u[2:n] - u[1:n - 1] ) / c]
    end
end

result = simulate(SpatialDiscretization, 1)
plot(result, "u", figure=1)

splice(uodd, ueven) = [if isodd(i); uodd[div(i,2)+1] else ueven[div(i,2)] end for i in 1:length(ueven)+length(uodd)]

@model SpatialDiscretization2 begin
    u = Var(size=(5,))
    ueven = Var(start=[1,2])
    uodd = Var(size=(3,))
    @equations begin
        u = splice(uodd, ueven)
        der(ueven) = - u[[2,4]] + u[[1,3]]
        uodd = ones(3)*time
    end
end

result = simulate(SpatialDiscretization2, 1, logTranslation=true, storeEliminated=false)
plot(result, "ueven", figure=2)

n = 100
@model SpatialDiscretization3 begin
    u = Var(size=(2*n+1,))
    ueven = Var(start=[i for i in 1:n])
    uodd = Var(size=(n+1,))
    @equations begin
        u = splice(uodd, ueven)
        der(ueven) = - u[2*(1:n)] + u[2(1:n)-1]
        uodd = ones(n+1)*time
    end
end

@static if VERSION < v"0.7.0-DEV.2005"
    @time result = simulate(SpatialDiscretization3, 1, logTranslation=true, storeEliminated=false)
    plot(result, "ueven", figure=3)
end

@model SpatialDiscretization4 begin
    u = Var(size=(2*n+1,))
    ueven = Var(start=[i for i in 1:n])
    uodd = Var(size=(n+1,))
    @equations begin
        u = splice(uodd, ueven)
        der(ueven) = ueven
        uodd = zeros(n+1)
    end
end

@time result = simulate(SpatialDiscretization4, 1, logTranslation=true, storeEliminated=false)
plot(result, "ueven", figure=4)

maskEven(n) = [if isodd(i); 0 else 1 end for i in 1:n]

@model SpatialDiscretization5 begin
    u = Float(start=[if iseven(i); div(i,2) else 0 end for i in 1:2*n+1], state=false)
    @equations begin
        broadcast(*, maskEven(2*n+1), der(u)) = u
    end
end

# @time result = simulate(SpatialDiscretization5, 1, logTranslation=true, storeEliminated=false, removeSingularities=false)
# plot(result, "u", figure=5)

end


text:
#=
In the 20×20 grid below, four numbers along a diagonal line have been marked in red.

08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70
67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21
24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72
21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95
78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48
The product of these numbers is 26 × 63 × 78 × 14 = 1788696.

What is the greatest product of four adjacent numbers in the same direction (up, down, left, right, or diagonally) in the 20×20 grid?
=#

using DelimitedFiles

function calc(x, i, j, adjacents)
  ar = Function[]
  push!(ar, (i,j)->i<size(x,1)-adjacents) # right
  push!(ar, (i,j)->j<size(x,2)-adjacents) # down
  push!(ar, (i,j)->j>=adjacents) # up
  push!(ar, (i,j)->i>=adjacents) # left
  push!(ar, (i,j)->i+adjacents<size(x,1) && j+adjacents<size(x,2)) # right diagonal
  push!(ar, (i,j)->i-adjacents>0 && j+adjacents<size(x,2)) # left diagonal
  vectors = Array[[1, 0], [0, 1], [0, -1], [-1, 0], [1, 1], [-1, 1]]
  maxwalk = 1
  for element in 1:length(ar)
    if ar[element](i,j)
      maxwalk = max(maxwalk, walk(x, i, j, vectors[element], adjacents))
    end
  end
  maxwalk
end

function walk(x, i, j, vector, adjacents)
  product = 1
  for a in 1:adjacents
      product *= x[j, i]
      i += vector[1]
      j += vector[2]
  end
  product
end

function product_adjacent_numbers()
  x = readdlm(dirname(@__FILE__()) * "/../Resources/p11.txt", ' ', Int64) # Load matrix from file
  adjacents = 4
  maxproduct = 1
  for i in 1:size(x,1), j in 1:size(x,2)
    maxproduct = max(calc(x, i, j, adjacents), maxproduct)
  end
  maxproduct
end
@time println(product_adjacent_numbers())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using Test
import ArchGDAL as AG

@testset "test_geotransform.jl" begin
    function f(gt, pixel, line)
        return [
            gt[1] + pixel * gt[2] + line * gt[3],
            gt[4] + pixel * gt[5] + line * gt[6],
        ]
    end
    gt1 = [
        -1111950.519667,
        463.3127165279167,
        0.0,
        6671703.118,
        0.0,
        -463.3127165279165,
    ]
    for (pixel, line) in ((0.0, 0.0), (3.0, 2.0), (15.5, 9.9))
        @test AG.applygeotransform(gt1, pixel, line) ≈ f(gt1, pixel, line)
    end

    gt2 = AG.invgeotransform(gt1)
    @test gt2 ≈ [
        2400.0,
        0.0021583694216166533,
        0.0,
        14400.0,
        0.0,
        -0.002158369421616654,
    ]

    gt3 = AG.composegeotransform(gt1, gt2)
    @test gt3 ≈ [0.0, 1.0, 0.0, 0.0, 0.0, 1.0]

    for (pixel, line) in ((0.0, 0.0), (3.0, 2.0), (15.5, 9.9))
        @test AG.applygeotransform(gt3, pixel, line) ≈ [pixel, line]
    end
end


text:
include("libavutil_h.jl")

#include("aes_ctr.jl")
#include("audio_fifo.jl")
#include("avutil.jl")
#include("buffer.jl")
#include("camellia.jl")
#include("cast5.jl")
#include("channel_layout.jl")
#include("des.jl")
#include("dict.jl")
#include("display.jl")
#include("downmix_info.jl")
#include("fifo.jl")
#include("file.jl")
#include("frame.jl")
#include("hash.jl")
#include("imgutils.jl")
include("log.jl")
#include("mastering_display_metadata.jl")
include("mem.jl")
#include("motion_vector.jl")
#include("opencl.jl")
#include("opt.jl")
include("pixdesc.jl")
#include("pixelutils.jl")
#include("pixfmt.jl")
#include("rational.jl")
#include("rc4.jl")
#include("replaygain.jl")
#include("samplefmt.jl")
#include("stereo3d.jl")
#include("tea.jl")
#include("threadmessage.jl")
#include("timecode.jl")
#include("tree.jl")
#include("twofish.jl")
#include("version.jl")
#include("xtea.jl")

function AVFrame()
    ns = fieldnames(AVFrame)
    ts = AVFrame.types
    parms = [T <: Ptr ? C_NULL : zero(T) for T in ts]

    fmt_pos = findfirst(ns, :format)
    parms[fmt_pos] = -one(ts[fmt_pos])

    pts_pos = findfirst(ns, :pts)
    parms[pts_pos] = AV_NOPTS_VALUE

    AVFrame(parms...)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
using ExBroadcast
using Test
using OffsetArrays, CUDA, StructArrays, Interpolations, Adapt
@testset "@tab.jl" begin
    a = randn(1000)
    b, c = @tab sincos.(a)
    @test b == sin.(a)
    @test c == cos.(a)
    b′ = @view b[1:999]
    @test_throws ArgumentError b′, c .= sincos.(a)
    b″ = OffsetArray(b,-1)
    @test_throws ArgumentError b″, c .= sincos.(a)
end

@testset "@mtb for OffsetVector" begin
    a = randn(1000)
    aᵒ = OffsetArray(a,-1)
    @mtb b = parent(aᵒ .+ 1)
    @test (a .+ 1) == b
end

@testset "@mtb for small array" begin
    a = [1]
    @mtb 2 a .+= 1
    @test a[1] == 2
end

@testset "@mtb for BitArray" begin
    a = randn(4096 * 6 - 2048)
    b = a .> 0
    @test @mtb (a .> 0) == b
end

@testset "CUDA supports" begin
    a = CUDA.randn(1000)
    b, c = @tab sincos.(a)
    @test b ≈ sin.(a)
    @test c ≈ cos.(a)
    d = StructArray{ComplexF32}((b, c))
    @test abs.(d) isa CuVector
    e = OffsetArray(d, -1) .+ 1
    @test e isa StructArray
    @test eachindex(e) == 0:999
end

@testset "Lazy test" begin
    a = ((i,k...) for i in 1:10, k in Iterators.product(1:10,1:2))
    @test identity.(a) == identity.(Lazy(a))
    a = (randn() for i in 1:100)
    b = Lazy(a) .+ [1 2] .- [1 2]
    @test b[:,1] != b[:,2]
    a = (randn() for i in 1:100, j in tuple())
    @test identity.(a) == identity.(Lazy(a))
end

@testset "CUDA Interp test" begin
    a = randn(Float32,101,101)
    itp = CubicSplineInterpolation((-50f0:50f0,-50f0:50f0), a)
    cuitp = adapt(CuArray, itp)
    res = cuitp.(-10:10,-10:10)
    @test res isa CuArray
    @test collect(res) ≈ itp.(-10:10,-10:10)

    itp = itp.itp
    a,b,c = preweight(itp,-1:0.01f0:1,-1:0.01f0:1)
    @test a.(b,c) == itp.(-1:0.01f0:1,-1:0.01f0:1)
end


text:
# https://turing.ml/dev/tutorials/02-logistic-regression/

# Import Turing and Distributions.
using Turing, Distributions

# Import RDatasets.
using RDatasets

# Import MCMCChains, Plots, and StatsPlots for visualizations and diagnostics.
using MCMCChains, Plots, StatsPlots

# We need a logistic function, which is provided by StatsFuns.
using StatsFuns: logistic

# Functionality for splitting and normalizing the data
using MLDataUtils: shuffleobs, stratifiedobs, rescale!

# Set a seed for reproducibility.
using Random
Random.seed!(0);

# Turn off progress monitor.
Turing.setprogress!(false)

# Import the "Default" dataset.
data = RDatasets.dataset("ISLR", "Default");

# Show the first six rows of the dataset.
first(data, 6)

# Convert "Default" and "Student" to numeric values.
# if else --> condition ? true_value : false_value
data[!, :DefaultNum] = [r.Default == "Yes" ? 1.0 : 0.0 for r in eachrow(data)]
data[!, :StudentNum] = [r.Student == "Yes" ? 1.0 : 0.0 for r in eachrow(data)]

# Delete the old columns which say "Yes" and "No".
select!(data, Not([:Default, :Student]))

# Show the first six rows of our edited dataset.
first(data, 6)

# summary
describe(data)

# train and split (; enable to change default argument `at`)
function split_data(df, target; at = 0.70)
    shuffled = shuffleobs(df)
    trainset, testset = stratifiedobs(row -> row[target], shuffled, p = at)
end

features = [:StudentNum, :Balance, :Income]
numerics = [:Balance, :Income]
target = :DefaultNum

trainset, testset = split_data(data, target, at = 0.05)

# Rescale variables for efficient parameter searching.
for feature in numerics
    μ, σ = rescale!(trainset[!, feature], obsdim=1) # rescale!() overwrite the memory
    rescale!(testset[!, feature], μ, σ, obsdim=1)
end

describe(trainset)
describe(testset)

# Turing requires data in matrix form, not dataframe
train = Matrix(trainset[:, features])
test = Matrix(testset[:, features])
train_label = trainset[:, target]
test_label = testset[:, target]

# Bayesian logistic regression (LR)
"""
logistic regression

    x is our set of independent variables;
    y is the element we want to predict;
    n is the number of observations we have; and
    σ is the standard deviation we want to assume for our priors.
"""
@model function logistic_regression(x, y, n, σ)
    intercept ~ Normal(0, σ)

    student ~ Normal(0, σ)
    balance ~ Normal(0, σ)
    income ~ Normal(0, σ)

    for i = 1:n
        v = logistic(intercept + student*x[i,1] + balance*x[i,2] + income*x[i,3])
        y[i] ~ Bernoulli(v)
    end
end

# Retrieve the number of observations.
n, _ = size(train)

# Sample using HMC.
m = logistic_regression(train, train_label, n, 1)
chain = sample(m, HMC(0.05, 10), MCMCThreads(), 1_500, 3)
describe(chain)
# size(chain)
plot(chain)

# the labels to use
l = [:student, :balance, :income]

# Use the corner function. Requires StatsPlots and MCMCChains.
corner(chain, l)

# Fortunately the corner plot appears to demonstrate unimodal distributions for each of our parameters, so it should be straightforward to take the means of each parameter's sampled values to estimate our model to make predictions.

# making prediction
function prediction(x::Matrix, chain, threshold)
    # Pull the means from each parameter's sampled values in the chain.
    intercept = mean(chain[:intercept])
    student = mean(chain[:student])
    balance = mean(chain[:balance])
    income = mean(chain[:income])

    # Retrieve the number of rows.
    n, _ = size(x)

    # Generate a vector to store our predictions.
    v = Vector{Float64}(undef, n)

    # Calculate the logistic function for each element in the test set.
    for i in 1:n
        num = logistic(intercept .+ student * x[i,1] + balance * x[i,2] + income * x[i,3])
        if num >= threshold
            v[i] = 1
        else
            v[i] = 0
        end
    end
    return v
end;

# Set the prediction threshold.
threshold = 0.07 # 0.5

# Make the predictions.
predictions = prediction(test, chain, threshold)

# Calculate MSE for our test set.
loss = sum((predictions - test_label).^2) / length(test_label)

# Calculate percentage of correctly predicted.
defaults = sum(test_label)
not_defaults = length(test_label) - defaults

predicted_defaults = sum(test_label .== predictions .== 1)
predicted_not_defaults = sum(test_label .== predictions .== 0)

println("Defaults: $defaults
    Predictions: $predicted_defaults
    Percentage defaults correct $(predicted_defaults/defaults)")

println("Not defaults: $not_defaults
    Predictions: $predicted_not_defaults
    Percentage non-defaults correct $(predicted_not_defaults/not_defaults)")

# 課題
# 行列でモデルのpriorや尤度を書く
# バーンインを捨てる


text:
using LBFGSB
using Test

include("driver1.jl")
include("driver2.jl")
include("driver3.jl")
include("wrapper.jl")


text:
using Compat

@static if VERSION >= v"0.7"
    using Libdl
end

const DEPS_FILE = joinpath(@__DIR__, "deps.jl")

if isfile(DEPS_FILE)
    rm(DEPS_FILE)
end

function write_depsfile(path)
    if Compat.Sys.iswindows()
        # When `path` gets written out to a file, it will escape any
        # backslashes, so we need to doubly escape them. If your path uses
        # forward slashes, this operation won't do anything.
        path = replace(path, "\\" => "\\\\")
    end
    open(DEPS_FILE, "w") do io
        println(io, "const libgurobi = \"$(path)\"")
    end
end

const ALIASES = [
    "gurobi81", "gurobi80",
    "gurobi75", "gurobi70"
]

paths_to_try = copy(ALIASES)

for a in ALIASES
    if haskey(ENV, "GUROBI_HOME")
        if Compat.Sys.isunix()
            push!(paths_to_try, joinpath(ENV["GUROBI_HOME"], "lib", string("lib", a, ".so")))
        end
        if Compat.Sys.iswindows()
            push!(paths_to_try, joinpath(ENV["GUROBI_HOME"], "bin", string(a, ".", Libdl.dlext)))
        end
        if Compat.Sys.isapple()
            push!(paths_to_try, joinpath(ENV["GUROBI_HOME"], "lib", string("lib", a, ".dylib")))
        end
    end
    # gurobi uses .so on OS X for some reason
    if Compat.Sys.isapple()
        push!(paths_to_try, string("lib$a.so"))
        push!(paths_to_try, string("lib$a.dylib"))
    end
end

found = false
for l in paths_to_try
    d = Libdl.dlopen_e(l)
    if d != C_NULL
        global found = true
        write_depsfile(l)
        break
    end
end

if !found && !haskey(ENV, "GUROBI_JL_SKIP_LIB_CHECK")
    error("Unable to locate Gurobi installation. Note that this must be downloaded separately from gurobi.com")
end


text:
# This file contains all the basic functions needed to run a quantum simulator
# using the stabiliser formalism
# The format used is as described in
# Improved Simulation of Stabilizer Circuits
# Scott Aaronson and Daniel Gottesman
# arXiv:quant-ph/0406196v5

# It is, effectively, a Julia Port of the CHP program by Scott Aaronson
# The original CHP program can be found at http://www.scottaaronson.com/chp

# Primarily, I have tried to write for clarity rather than efficiency.
# (Although thanks to the algorithms by Aaronson/Gottesman its fast enough for hundreds of qubits)

# I am partially through rationalising some things for ease of REPL
# Currently  you

# include("Initial.l")
# using Main.CHP

# Then you can e.g.
# t = Tableau(4)
# hadamard(t,2,4) # Applies hadamard to qubits 2 and 4
# cnot(t,2,1)
# cnot(t,2,3)
# cnot(t,4,3)
# measure(t,1,2)
# print(kets(t))

# Copyright Robin Harper 2015-2020

"""
    Tableau(n::Integer)

 sets up and returns the a CHP Tableau for n qubits.

 This is based on the formalism given by:
 *Improved Simulation of Stabilizer Circuits*,
 Scott Aaronson and Daniel Gottesman,
 arXiv:quant-ph/0406196v5

 The initial tableau represents a |00...0\$\\rangle\$ ket in the stabiliser state
 This stabilises with "Z" and anti-stabilises with "X"
 So the tableau starts off with an Identity in the 'anti commute' top half in the X section
 and identity in the commuting (bottom) half in the Z section.

 For the purposes of this port, the tableau is exactly replicated as per the paper
 i.e. the "state" (Tableau.state) is an Int32 array (used as a bit array)
 containing the following information.

```
 x11   .....  x1n | z11   ...      z1n | r1
  .    \\       .  |  .    \\          . |  .
  .     \\      .  |  .     \\         . |  .     Destabilisers
  .      \\     .  |  .      \\        . |  .
 xn1      \\   xnn | zn1      \\      znn| rn
 ______________________________________________
 x(n+1)1. x(n+1)n | z(n+1) ... z(n+1)n | r(n+1)
  .    \\       .  |  .      \\        . |  .
  .     \\      .  |  .       \\       . |  .     Stabilisers
  .      \\     .  |  .        \\      . |  .
 x(2n)1   \\x(2n)n | z(2n)1     \\ z(2n)n| r(2n)
```
Set Tableau.showRaw to true to see the underlying state as a matrix.

See also: [`cnot`](@ref), [`hadamard`](@ref), [`phase`](@ref), [`measure`](@ref), [`X`](@ref), [`Z`](@ref)

Related: [`cliffordToTableau`](@ref)
"""
function Tableau(x::Integer)
    t  = Tableau(setup(x),x,false,["initialise($x)"],[],true)
    t.executeCommands = [Expr(:call,:initialise,:_t)]
    return t
end

"""
    Tableau(cs::Array{String,1})

    creates a tableau from a string of commands.
    String parsing is not very robust, it assumes they are in exactly
    the form stored internally in a tableau.
    e.g

    >t = Tableau(2)
    >hadamard(t,2)
    >phase(t,1)

    > t.commands
    3-element Array{String,1}:
     "initialise(2)"
     "hadamard(2)"
     "phase(1)"

     >x=Tableau(t.commands)
"""
function Tableau(cs::Array{String,1})
    m = match(r"initialise\((.*)\)",cs[1])
    if (m==nothing)
        @warn("Command string needs to start with an intialise(::Integer) command")
        return nothing
    end
    bits =parse(Int,(m.captures[1]))
    t = Tableau(Juqst.setup(bits),bits,false,["initialise($bits)"],[],true)
    t.executeCommands = [Expr(:call,:initialise,:_t)]
    for c in cs[2:end]
        m=match(r"phase\((.*)\)",c)
        if (m!=nothing)
           bit = parse(Int,(m.captures[1]))
           phase(t,bit)
        else
            m=match(r"hadamard\((.*)\)",c)
            if (m!=nothing)
               bit = parse(Int,(m.captures[1]))
               hadamard(t,bit)
            else
               m=match(r"cnot\((.*),(.*)\)",c)
               if (m!=nothing)
                  cbit = parse(Int,(m.captures[1]))
                  tbit = parse(Int,(m.captures[2]))
                  cnot(t,cbit,tbit)
                else
                    m=match(r"measure\((.*)\)",c)
                    if (m!=nothing)
                        bit = parse(Int,(m.captures[1]))
                        measure(t,bit)
                    end
                end

            end
        end
    end
    return t
end

"""
    reinitialise(t::Tableau)

Reinitialise the Tableau, setting it support the |000000> state.
Resets all commands associated with the Tableau.

See also: [`initialise`](@ref)
"""
function reinitialise(t::Tableau)
    initialise(t)
    t.trackCommands = true
    t.commands = ["initialise($(t.qubits))"]
    t.executeCommands = [Expr(:call,:initialise,:_t)]
    return t
end

function initialize(t::Tableau)
    initialise(t)
end

"""
    initialise(t::Tableau)

Initialises the Tableau, setting it support the |000000> state.
However, unlike reinitialse it does not change the commands associated with the Tableau.

See also: [`reinitialise`](@ref)
"""
function initialise(t::Tableau)
    t.state = setup(t.qubits)
    return t
end

Base.show(io::IO,z::Tableau) = print(io,pformat(z))
Base.show(stream, ::MIME"text/latex", x::Tableau) = write(stream, nicePrintTableauState(x))

"""
Used by `show` to 'text' print out the Tableau based on the state.
"""
function pformat(tab::Tableau)::String
    # prints the Heisnberg representation of the state array supplied
    state = tab.state

    toOutput::String = "\n"
    n=div(size(state,1),2)
    #  1 in a and a+n = Y
    #  else 1 in a = X
    #  else 1 in a+n=Z
    #  otherwise I
    #  last bit is sign. (exponent of -1).
    for row=1:n
      if state[row,2*n+1] == 1
        toOutput *= "-"
      else
        toOutput *= "+"
      end
      for col=1:n
         toOutput *= decode(state[row,col],state[row,col+n])
      end
      toOutput *= "\n"
    end
    for dot=1:n+1
        toOutput *= "-"
    end
    toOutput *= "\n"
    for row=n+1:2*n
      if state[row,2*n+1] == 1
          toOutput *= "-"
      else
          toOutput *= "+"
      end
      for col=1:n
          toOutput *= decode(state[row,col],state[row,col+n])
      end
      toOutput *= "\n"
    end
    return toOutput
end

"""
    cnot(t::Tableau,control::Integer,target::Integer,showOutput=false)

Performs a CNOT (controlled not) on the tableau, with the specified qubits
as control and target.

If showOutput is passed as true the state of the tableau subsequent to the cnot
is printed.

See also: [`hadamard`](@ref), [`phase`](@ref), [`measure`](@ref), [`X`](@ref), [`Z`](@ref)
"""
function cnot(t::Tableau,control::Integer,target::Integer,showOutput=false)
  # From control a to taget b
  # You can supress the output by calling with the fourth parameter false
  # otherwise it will "output" the stabiliser state.
  # note the state is changed, (in Julia this should really be cnot!)
  state = t.state
  n = t.qubits
  endC=2*n+1

  if control < 1 || control > n
    print("Control qubit $(control) out of range of state which has $(n) qubits\n")
    throw(BoundsError())
    return
   end
  if target < 1 || target > n
    print("Target qubit $(target) out of range of state which has $(n) qubits\n")
    throw(BoundsError())
    return
  end
  if t.trackCommands
    append!(t.commands,["cnot($control,$target)"])
    append!(t.executeCommands,[Expr(:call,:cnot,:_t,control,target)])
  end
  for i = 1:(2*n)
      ri=state[i,endC]
      xia = state[i,control]
      xib = state[i,target]
      zia = state[i,control+n]
      zib = state[i,target+n]
      state[i,endC] = j_xor(ri,xia*zib*(j_xor(xib,j_xor(zia,1))))
      state[i,target]=j_xor(xib,xia)
      state[i,control+n]=j_xor(zia,zib)
      #print("a is $(a) b is $(b), xia is $(xia) xib is $(xib), setting $([i,b]) to $(state[i,b])\n")
   end
   if showOutput
       print(t)
   end
end

"""
    hadamard(t::Tableau,qubit::Integer,showOutput::Bool=false)

Performs a hadamard gate on the specified qubit of the tableau.

See also: [`cnot`](@ref), [`phase`](@ref), [`measure`](@ref), [`X`](@ref), [`Z`](@ref)
"""
function hadamard(t::Tableau,qubit::Integer,showOutput::Bool=false)
  # Apply a hadamard gate to quibit a in the state, state.
  # again note it changes the state
  state = t.state
  n = t.qubits
  endC=2*n+1

  if qubit < 0 || qubit > n
    print("Qubit $(qubit) out of range of state (tableau only has $(n) qubits)\n")
    throw(BoundsError())
    return
  end
  if t.trackCommands
    append!(t.commands,["hadamard($qubit)"])
    append!(t.executeCommands,[Expr(:call,:hadamard,:_t,qubit)])
  end
  for i = 1:(2*n)
      ri=state[i,endC]
      xia = state[i,qubit]
      zia = state[i,qubit+n]
      state[i,endC]=j_xor(ri,xia*zia)
      state[i,qubit]=zia
      state[i,qubit+n]=xia
   end
   if showOutput
      print(t)
   end
end

"""
    hadamard(t::Tableau,qubits::Integer...;showOutput::Bool = false)

Peforms hadamards on **each** of the qubits supplied.

# Examples
```julia-repl
julia> hadamard(t,2,3,1)
```

See also: [`cnot`](@ref) [`phase`](@ref), [`measure`](@ref), [`X`](@ref), [`Z`](@ref)
"""
function hadamard(t::Tableau,qubits::Integer...;showOutput::Bool = false)
    for i in qubits
        hadamard(t,i,false)
    end
    if showOutput
       print(t)
    end
end

"""
    phase(t::Tableau,qubit::Integer,showOutput::Bool=false)

Performs a phase gate on the specified qubit of the tableau.

See also: [`cnot`](@ref), [`hadamard`](@ref), [`measure`](@ref), [`X`](@ref), [`Z`](@ref)
"""
function phase(t::Tableau,qubit::Integer,showOutput::Bool=false)
  # Apply the phase gate to qubit a
  # again changes the state (phase!)
  state = t.state
  n = t.qubits
  endC=2*n+1

  if qubit < 0 || qubit > n
    print("Qubit $(qubit) out of range of state (tableau only has $(n) qubits)\n")
    throw(BoundsError())
    return
   end
  # so we need two loops!
  if t.trackCommands
      append!(t.commands,["phase($qubit)"])
      append!(t.executeCommands,[Expr(:call,:phase,:_t,qubit)])
  end
  for i = 1:(2*n)
      ri=state[i,endC]
      xia = state[i,qubit]
      zia = state[i,qubit+n]
      state[i,endC]=j_xor(ri,xia*zia)
   end
   for i = 1:(2*n)
      ri=state[i,endC]
      xia = state[i,qubit]
      zia = state[i,qubit+n]
      state[i,qubit+n]=j_xor(xia,zia)
   end
   if showOutput
      print(t)
   end
end

"""
    phase(t::Tableau,qubits::Integer...;showOutput::Bool = false)

Peforms phase on each of the qubits supplied.

# Examples
```julia-repl
julia> phase(t,2,3,1)
```

See also: [`cnot`](@ref) [`hadamard`](@ref), [`measure`](@ref), [`X`](@ref), [`Z`](@ref)
"""
function phase(t::Tableau,qubits::Integer...;showOutput::Bool = false)
    for i in qubits
        phase(t,i,false)
    end
    if showOutput
       print(t)
    end
end

"""
    Z(t::Tableau,qubit::Integer,showOutput::Bool=false)

Convenience function to perform a 'Z' gate on the tableau. Performs two phase gates.

See also: [`cnot`](@ref), [`hadamard`](@ref), [`measure`](@ref), [`X`](@ref), [`phase`](@ref)
"""
function Z(t::Tableau,qubit::Integer,showOutput::Bool=false)
    phase(t,qubit,false)
    phase(t,qubit,false)
    if showOutput
       print(t)
    end
end

"""
    Z(t::Tableau,qubits::Integer...;showOutput::Bool = false)

Convenience function to perform a 'Z' gate on the tableau on each of the qubits supplied. Performs the gate by doing two phase gates.

# Examples
```julia-repl
julia> Z(t,2,3,1)
```

See also: [`cnot`](@ref) [`hadamard`](@ref), [`measure`](@ref), [`X`](@ref), [`Z`](@ref)
"""
function Z(t::Tableau,qubits::Integer...;showOutput::Bool = false)
    for i in qubits
        Z(t,i)
    end
    if showOutput
       print(t)
    end
end

"""
    X(t::Tableau,qubit::Integer,showOutput::Bool=false)

Convenience function to perform an 'X' gate on the tableau. Performs hadamard, two phase gates and a hadamard.

See also: [`cnot`](@ref), [`hadamard`](@ref), [`measure`](@ref), [`Z`](@ref), [`phase`](@ref)
"""
function X(t::Tableau,qubit::Integer,showOutput::Bool=false)
    hadamard(t,qubit,false)
    phase(t,qubit,false)
    phase(t,qubit,false)
    hadamard(t,qubit,false)
    if showOutput
       print(t)
    end
end

"""
    X(t::Tableau,qubits::Integer...;showOutput::Bool = false)

Convenience function to perform an 'X' gate on the tableau on each of the qubits supplied.
Performs the gate by doing a hadamard, two phase gates and a hadamard on each of the qubits.

# Examples
```julia-repl
julia> X(t,2,3,1)
```

See also: [`cnot`](@ref) [`hadamard`](@ref), [`measure`](@ref), [`phase`](@ref), [`Z`](@ref)
"""
function X(t::Tableau,qubits::Integer...;showOutput::Bool = false)
    for i in qubits
        X(t,i,false)
    end
    if showOutput
       print(t)
    end
end

"""
    measure(t::Tableau,qubit::Integer)

Performs a measurement on the supplied qubit. If the Tableau is not in a Z eigenstate, it
will (pseudo) randomly choose the result and collapse the Tableau approriately.

"""
function measure(t::Tableau,qubit::Integer)
  state = t.state
  n=t.qubits
  endC=2*n+1

  if qubit < 0 || qubit > n
     print("Qubit $(qubit) out of range of state (tableau only has $(n) qubits)\n")
     throw(BoundsError())
     return
  end
  if t.trackCommands
      append!(t.commands,["measure($qubit)"])
      append!(t.executeCommands,[Expr(:call,:measure,:_t,qubit)])
  end
  if sum(state[(n+1):end,qubit]) > 0
     # we have a one on the xia somewhere.
     # This makes it a random measurement
     p = n+argmax(state[(n+1):end,qubit]); # Find the first stabiliser with an X in this qubit.
     # set p-nth to equal p, (because of the measurement the stabiliser
     # has now become a destabiliser. I.e. it will anti-commute with the measurement, but
     # commutes with everything else).
     state[p-n,:]=state[p,:]
     # Now set up the measurment 'stabiliser' this can be just IIs with a Z in the appropriate
     # qubit (and the 'r' set to 1 or 0 depending on the random outcome)
     value = round.(Integer,rand()) # check distribution (did we roll a 1 or a 0)
     state[p,:]=zeros(2*n+1) # All of them zero
     state[p,endC]=value # Set the end approrpiately.
     state[p,n+qubit]=1 # Set the 'z' qubit
     # Finally we need to make everything commute with the measurment.
     # we do this my right multiplying everything that doesn't commute with the de-stabiliser.
     for h=1:2n # h is what we use in rowmult!
       if (h!=p-n) && (state[h,qubit]==1) # We have an X (or Y) in the relevant qubit
         rowmult!(state,h,p-n)
       end
     end
    return value # the result of the measurement.
  else
      # print("Deterministic\n\n")
       #first set the 2n+1 row i.e. scratch space to be zero
       scratch=zeros(Int32,2n+1)' # hey 2n works dont you know
       state = vcat(state,scratch) # add it on.
       for i = 1:n
         if state[i,qubit]==1
           rowmult!(state,2n+1,i+n)
         end
       end
       value = state[2n+1,endC]
       state=state[1:end-1,:]
       return value
   end
end

"""
    measure(t::Tableau,qubits::Integer...;showOutput = false)

Performs a measurement on the supplied qubits. Results are reproted in a dictionary of outcomes.

# Examples
```julia-repl
julia> using Main.CHP

julia> t = Tableau(4)

+XIII
+IXII
+IIXI
+IIIX
-----
+ZIII
+IZII
+IIZI
+IIIZ

julia> measure(t,1,2)
Dict{Any,Any}(2=>0,1=>0)
julia>
```
"""
function measure(t::Tableau,qubits::Integer...;showOutput = false)
    outcomes = Dict()
    for i in qubits
        outcomes[i] = measure(t,i)
    end
    print(outcomes)
    if showOutput
       print(t)
    end
end

"""
    kets(tableau::Tableau)

returns a string showing the kets of the state stabilised by the tableau.
"""
function kets(originalTableau::Tableau)

    #Use gaussian Elimination to provide a minimal set of generators
    # containing X's and Y's in upper triangular form
    ss1 = gaussianElimination(originalTableau.state)

    state=ss1[1]
    n=ss1[2]

    start=seed(state,n)
    scratchState=copy(start) # the first nqubits are the destabilisers
    outputString = printBasisState(vec(scratchState))

    # So the idea is we have the stabilisers with an X in them in upper triangular form
    # Each X represents a "1" in the Ket
    # But we need to generate each of the states stabilised by these Paulis.
    # For the first row, the ones generated are where there is an X in the first row.
    # On the second row, the ones generated are the second row, and the second row multiplied by the first.
    # On the third row, its the third row, third times first, third times second and third times first and second
    # etc.
    # Aaronsen has a really neat algorithm to do this, but its not easily understood.
    # I have used this, slighly less efficient, but (I hope) more comprehenisible algorithm
    # For a particular row (say we are looking at our third row with an X in it (n 3)
    # loop through the numbers 4:7 (being 2^(3-1) to (2^3)-1)
    # In binary, this loops as 100, 101, 110, 111, So we can see This is applying
    # Row 3 with each possible combination of Rows 1 and 2.
    # Then we need to multiply the rows together if 2^(row-1) (i.e. row 2 - equates to bit 2) is not zero when
    # logically anded with the count.

    number=2^(n-1)
    nqubits=div(size(state,1),2)
    for count = 0:(2^(n-1)-2)
        c2=xor(count,(count+1))
        for rows=0:(n-2)
            if (2^(rows) & c2) > 0
               scratchState=rowmult(scratchState,state[nqubits+rows+1,:])
            end
        end
        # Note the use of vec here, just makes sure that the scratch state is a vectdor
        # This was only an issue for the |00> one.
        outputString *= printBasisState(vec(scratchState))
    end
    return outputString
end

"""
    nicePrintTableauState(t::Tableau)

Used by show to print out the latex equation of the current tableau for e.g. Jupyter notebook
    In theory, if you want the latex equation call this and post into tex!

"""
function nicePrintTableauState(t::Tableau)
  state = t.state
  n = t.qubits
  toOuptut = ""
  if !t.showRaw
      start = string("\\begin{equation*} \\left(\\begin{array}{l|c*{$n}c}")
      for row=1:n
          if state[row,2*n+1] == 1
            toOutput = "-&"
          else
            toOutput = "+&"
          end
          start = "$start $toOutput &"
          for col=1:n
             start = "$start $(decode(state[row,col],state[row,col+n])) "
         if (col != n)
            start = "$start & "
          end
          end
          start = "$start \\\\\n"
        end
        start = "$start \\hline\n"
        for row=n+1:2*n
          if state[row,2*n+1] == 1
            toOutput = "-&"
          else
            toOutput = "+&"
          end
          start = "$start $toOutput &"
          for col=1:n
            start = "$start $(decode(state[row,col],state[row,col+n])) "
            if (col != n)
              start = "$start & "
            end

          end
          start = "$start \\\\\n"
        end
        start = "$start \\end{array}\\right)\\\\\\end{equation*}"
      return start
    else
        start = string("\\begin{equation*} \\begin{array}{l$(join(["c" for i =1:n]))|$(join(["c" for i =1:n]))|c}")
        toshift = round(Integer,n/2,RoundUp)
        start = "$start\n$(join(["&" for i =1:toshift]))\\textbf{X}$(join(["&" for i =1:(n-toshift)]))$(join(["&" for i =1:toshift]))\\textbf{Z}$(join(["&" for i =1:(n-toshift)]))\\\\\n"
        for row=1:n
          start = "$start D$(row): &"
          for col=1:(2n)
             start = "$start $(state[row,col]) &"
          end
          start = "$start $(state[row,2n+1]) \\\\\n"
        end
        start = "$start \\hline\n"
        for row=n+1:2*n
            start = "$start S$(row-n): &"
             for col=1:(2n)
             start = "$start $(state[row,col]) &"
          end
          start = "$start $(state[row,2n+1])\\\\\n"
        end
        start = "$start \\end{array}\\\\\\end{equation*}"
        return start
    end
end

# ============== Internal Functions ==========================
function setup(n)
  # sets up and returns the initial state vector for n qubits.
  # Heisenberg representation
  # We are representing a |00...0> ket in the stabiliser state
  # This stabilisises with "Z" and anti-stabilises with "X"
  # So the tableau starts off with an Identity in the Anti commute top X section
  # and identity in the commuting - Z section.

  # For the purposes of this port, the tableau is exactly replicated as drawn
  # i.e. the "state" is an Int32 array (should really be just binary)
  # of the form
  #
  #
  # x11   .....  x1n | z11   ...      z1n | r1
  #  .    \       .  |  .    \          . |  .
  #  .     \      .  |  .     \         . |  .     Destabilisers
  #  .      \     .  |  .      \        . |  .
  # xn1      \   xnn | zn1      \      znn| rn
  # ______________________________________________
  # x(n+1)1. x(n+1)n | z(n+1) ... z(n+1)n | r(n+1)
  #  .    \       .  |  .      \        . |  .
  #  .     \      .  |  .       \       . |  .     Stabilisers
  #  .      \     .  |  .        \      . |  .
  # x(2n)1   \x(2n)n | z(2n)1     \ z(2n)n| r(2n)

  top = hcat(Matrix{Int32}(I,n,n),zeros(Int32,n,n),zeros(Int32,n,1))
  bottom = hcat(zeros(Int32,n,n),Matrix{Int32}(I,n,n),zeros(Int32,n,1))
  state = vcat(top,bottom)
end

# This is just mod 2 arithmetic.
function j_xor(a,b)
  return mod.(a+b,2)
end

function rowswap!(alteredState,i,k)
      temp = alteredState[i,:]
      alteredState[i,:]=alteredState[k,:]
      alteredState[k,:]=temp
end

"""
    cliffordPhase(i,k)

    Pass in two vectors representing rows in a Tableau.

    Returns the exponent i is raised to (0,1,2 or 3) if i is left multiplied by k.
    ---
    Note: Assumes that the 'r' part of i (in position 2n+1) has ALREADY been changed from
    an exponent of -1 (which is how it is in the Tableau) to an exponent of i
    (by multiplying it by 2). It also assumes k is still in the range of 0..1 i.e. an
    exponent of -1
"""
function cliffordPhase(i,k)
  # Returns the phase (0,1,2,3) when row i is left multiplied by row k
  # println("CliffordPhase $i,$k")
  e=0
  n=div(length(i),2)
  for j=1:n
    if k[j] == 1 && k[j+n]==0 # its an X
      if i[j] == 1 && i[j+n] == 1 # its a Y
        e+=1 # XY=iZ
      elseif   i[j] == 0 && i[j+n] == 1 # its a Z
        e+=-1 # XZ=-iY
      end
    elseif  k[j] == 1 && k[j+n]==1 # its a Y
      if i[j] == 0 && i[j+n] == 1       # its a Z
        e+=1 # YZ=iXZ
      elseif i[j] == 1 && i[j+n] == 0 # its a X
        e+=-1 # YX=-iZ
      end
    elseif  k[j] == 0 && k[j+n]==1 # its a Z
      if i[j] == 1 && i[j+n] == 0       # its a X
        e+=1 # ZX=iY
      elseif i[j] == 1 && i[j+n] == 1 # its a Y
        e+=-1 # ZY=-iX
      end
    end
  end

  e = e + i[2*n+1] + 2*k[2*n+1]
  e = e%4
  if (e < 0)
    e+=4
  end

  return e
end

"""
    rowmult!(stateArray,i,k)

    stateArray - the state of a Tableau
    i the (integer) row of the stateArray that will hold the result of i left multiplied by
      row pointed to by k.
    k the (integer) row with which i is to be left multiplied.

    Note: assumes that both 'r' components of the Tableau rows are exponents of -1
    use _rowmult! if i has already been lifted up to be an exponent of i.
"""
function rowmult!(stateArray,i,k)
    n=div(size(stateArray,1),2)
    # Note cliffordPhase assumes the first row has the 'r' component changed to represent
    # a phase of i, we do this here by multiplying it by 2
    # and then after the row multiplication dividing it by 2.
    stateArray[i,2n+1] *=2
    _rowmult!(stateArray,i,k)
    stateArray[i,2n+1] = div(stateArray[i,2n+1],2)

end

function _rowmult!(stateArray,i,k)
    n=div(size(stateArray,1),2)
    # print(stateArray[i,:],stateArray[k,:])
    phaseBack = cliffordPhase(stateArray[i,:],stateArray[k,:])
    if (phaseBack  == 1 || phaseBack  ==3)
        println("phaseBack $phaseBack WAS 1 OR 3")
    end
    stateArray[i,:] = xor.(stateArray[i,:],stateArray[k,:])
    stateArray[i,2*n+1] = phaseBack
end

"""
    rowmult(i,k)

    i,k two vectors, equivalent to rwos in a Tableau.
    Left multiplies i with k and returns the result.
    Note it assumes that the 'r' in i is an exponent of i
    Whereas the 'r' in k is an exponent of -1.
"""
function rowmult(i,k) # supply two vectors and get the multiplied vector out.
  n=div(length(i),2)
  temp = cliffordPhase(i,k)
  tempState = [xor(i[c],k[c]) for c =1:length(i)]
  tempState[length(i)]=temp
  return tempState
 end

#returns a tableau state and a size where gaussian elimination has
#been carried out on the supplied state, in order to
#provide a minimal set of generators containing X's and Y's in
#upper triangular form

function gaussianElimination(state)
  alteredState = copy(state)
  #println(alteredState)
  n=div(size(alteredState,1),2)
  i = n+1 # first row of commuting generators
  k=0
  k2=0
  j=0
  for j=1:n
    found = 0;

    for kidx=i:2*n
      # Find a generator containing X in the jth column
      if alteredState[kidx,j] == 1
        found = 1
        k = kidx
        break
      end
    end
    if found  == 1
      #swap row with the row pointed to by i
      if (i!=k)
        rowswap!(alteredState,i,k)
        rowswap!(alteredState,i-n,k-n) # swap the non-commutators as well
      end
      # then use the row to eliminate X from that bit in the rest of the tableau
      for k2=i+1:2*n
        if alteredState[k2,j] == 1
          #println("Before, $(k2) and j is $j and i is $i")
          #println(alteredState)
          #print("ROWMULTA $k,$i\n")
          rowmult!(alteredState,k2,i)
          #println(alteredState)
          #println("Doing $(i-n) and $(k2-n)")
          #print("ROWMULTB $(i-n) $(k2-n)\n")
          rowmult!(alteredState,i-n,k2-n)
          #println("Row mult $k2, $i $alteredState")
        end
        #println(alteredState)
      end
      #println(alteredState)
      i+=1
    end
  end
  gen = i - n
  # The first gen generators are X/Ys and in quasi upper triangular form.
  for j=1:n
    found = 0
    for kidx=i:2*n
      if alteredState[kidx,n+j] == 1 # we have found a Z in the jth bit
        k = kidx
        found = 1
        break
      end
    end
    if found == 1
      if (i!=k)
        rowswap!(alteredState,i,k)
        rowswap!(alteredState,i-n,k-n)
      end
      for k2=i+1:2*n
        if (alteredState[k2,j+n] ==1 ) # z in this "bit"
          #print("ROWMULTC $k2,$i\n")

          rowmult!(alteredState,k2,i)
          #print("ROWMULTD $(i-n),$(k2-n)\n")

          rowmult!(alteredState,i-n,k2-n)
        end
      end
      i+=1
    end

  end
  return (alteredState,gen)
end

function printBasisState(sS)
    toPrint = ""
    nqubits=div(size(sS,1),2)
    #println("\npBS: $sS and $nqubits and $(size(sS)) e=$(sS[2*nqubits+1])\n")
    e=sS[2*nqubits+1]
    for i=1:nqubits
      if sS[i]==1 && sS[nqubits+i]==1 # its a Y
        e = (e+1) % 4
      end
    end
    if (e==0) toPrint *= "+|"
    elseif e== 1 toPrint *= "+i|"
    elseif e==2 toPrint *= "-|"
    else toPrint *= "-i|"
    end
    for i=1:nqubits
      if sS[i]==1
        toPrint *= "1"
      else
        toPrint *= "0"
      end
    end
    toPrint *= "> \n"
    return toPrint
end

function seed(state,n)
  nqubits=div(size(state,1),2)
  scratchState = zeros(Int32,1,2*nqubits+1)
   for i=2*nqubits:-1:(nqubits+n)
    min = nqubits
    f = 2*state[i,2*nqubits+1]
    # count down through the qubits.
    for j=nqubits:-1:1
      # if the z part of this is 1
      if state[i,nqubits+j] == 1
        min = j
        if scratchState[1,j]==1
          f = (f+2)%4
        end
      end
    end
    if f==2
      #print("F was 2 so we need to flip the x");
      scratchState[1,min] = xor(scratchState[1,min],1)

    end
  end
  return scratchState
end

function decode(z,x)
    if z==1
        if x == 1
            return "Y"
        else
            return "X"
        end
    elseif x == 1
        return "Z"
    else
        return "I"
    end
end


text:
using Dates
using GitHub
using Pkg
using Printf
using RegistryCI
using Test
using TimeZones

const AutoMerge = RegistryCI.AutoMerge

include("automerge-integration-utils.jl")

AUTOMERGE_INTEGRATION_TEST_REPO = ENV["AUTOMERGE_INTEGRATION_TEST_REPO"]::String
TEST_USER_GITHUB_TOKEN = ENV["BCBI_TEST_USER_GITHUB_TOKEN"]::String
GIT = "git"
auth = GitHub.authenticate(TEST_USER_GITHUB_TOKEN)
whoami = RegistryCI.AutoMerge.username(auth)
repo_url_without_auth = "https://github.com/$(AUTOMERGE_INTEGRATION_TEST_REPO)"
repo_url_with_auth = "https://$(whoami):$(TEST_USER_GITHUB_TOKEN)@github.com/$(AUTOMERGE_INTEGRATION_TEST_REPO)"
repo = GitHub.repo(AUTOMERGE_INTEGRATION_TEST_REPO; auth = auth)
@test success(`$(GIT) --version`)
@info("Authenticated to GitHub as \"$(whoami)\"")

close_all_pull_requests(repo; auth = auth, state = "open")
delete_stale_branches(repo_url_with_auth; GIT = GIT)

with_master_branch(templates("master_1"), "master"; GIT = GIT, repo_url = repo_url_with_auth) do master_1
    with_feature_branch(templates("feature_1"), master_1; GIT = GIT, repo_url = repo_url_with_auth) do feature_1
        title = "New package: Requires v1.0.0"
        head = feature_1
        base = master_1
        params = Dict("title" => title,
                      "head" => head,
                      "base" => base)
        pr_1_1 = GitHub.create_pull_request(repo; auth = auth, params = params)
        pr_1_1 = wait_pr_compute_mergeability(repo, pr_1_1; auth = auth)
        @test pr_1_1.mergeable
        with_pr_merge_commit(pr_1_1, repo_url_without_auth; GIT = GIT) do build_dir
            withenv("AUTOMERGE_GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN,
                    "TRAVIS_BRANCH" => master_1,
                    "TRAVIS_BUILD_DIR" => build_dir,
                    "TRAVIS_EVENT_TYPE" => "pull_request",
                    "TRAVIS_PULL_REQUEST" => string(pr_1_1.number),
                    "TRAVIS_PULL_REQUEST_SHA" => string(AutoMerge.pull_request_head_sha(pr_1_1))) do
                AutoMerge.run(;
                              merge_new_packages = true,
                              merge_new_versions = true,
                              new_package_waiting_period = Minute(typemax(Int32)),
                              new_version_waiting_period = Minute(typemax(Int32)),
                              registry = AUTOMERGE_INTEGRATION_TEST_REPO,
                              authorized_authors = String[whoami],
                              master_branch = master_1,
                              master_branch_is_default_branch = false)

            end
            withenv("AUTOMERGE_GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN,
                    "TRAVIS_BRANCH" => master_1,
                    "TRAVIS_BUILD_DIR" => build_dir,
                    "TRAVIS_EVENT_TYPE" => "cron",
                    "TRAVIS_PULL_REQUEST" => "false",
                    "TRAVIS_PULL_REQUEST_SHA" => "") do
                AutoMerge.run(;
                              merge_new_packages = true,
                              merge_new_versions = true,
                              new_package_waiting_period = Minute(typemax(Int32)),
                              new_version_waiting_period = Minute(typemax(Int32)),
                              registry = AUTOMERGE_INTEGRATION_TEST_REPO,
                              authorized_authors = String[whoami],
                              master_branch = master_1,
                              master_branch_is_default_branch = false)

            end
        end
    end
end

with_master_branch(templates("master_2"), "master"; GIT = GIT, repo_url = repo_url_with_auth) do master_2
    with_feature_branch(templates("feature_2"), master_2; GIT = GIT, repo_url = repo_url_with_auth) do feature_2
        title = "New version: Requires v2.0.0"
        head = feature_2
        base = master_2
        params = Dict("title" => title,
                      "head" => head,
                      "base" => base)
        pr_2_2 = GitHub.create_pull_request(repo; auth = auth, params = params)
        pr_2_2 = wait_pr_compute_mergeability(repo, pr_2_2; auth = auth)
        @test pr_2_2.mergeable
        with_pr_merge_commit(pr_2_2, repo_url_without_auth; GIT = GIT) do build_dir
            withenv("AUTOMERGE_GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN,
                    "TRAVIS_BRANCH" => master_2,
                    "TRAVIS_BUILD_DIR" => build_dir,
                    "TRAVIS_EVENT_TYPE" => "pull_request",
                    "TRAVIS_PULL_REQUEST" => string(pr_2_2.number),
                    "TRAVIS_PULL_REQUEST_SHA" => string(AutoMerge.pull_request_head_sha(pr_2_2))) do
                AutoMerge.run(;
                              merge_new_packages = true,
                              merge_new_versions = true,
                              new_package_waiting_period = Minute(typemax(Int32)),
                              new_version_waiting_period = Minute(typemax(Int32)),
                              registry = AUTOMERGE_INTEGRATION_TEST_REPO,
                              authorized_authors = String[whoami],
                              master_branch = master_2,
                              master_branch_is_default_branch = false)

            end
            withenv("AUTOMERGE_GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN,
                    "TRAVIS_BRANCH" => master_2,
                    "TRAVIS_BUILD_DIR" => build_dir,
                    "TRAVIS_EVENT_TYPE" => "cron",
                    "TRAVIS_PULL_REQUEST" => "false",
                    "TRAVIS_PULL_REQUEST_SHA" => "") do
                AutoMerge.run(;
                              merge_new_packages = true,
                              merge_new_versions = true,
                              new_package_waiting_period = Minute(typemax(Int32)),
                              new_version_waiting_period = Minute(typemax(Int32)),
                              registry = AUTOMERGE_INTEGRATION_TEST_REPO,
                              authorized_authors = String[whoami],
                              master_branch = master_2,
                              master_branch_is_default_branch = false)

            end
        end
    end
end


text:
using ADVUMPS
using Test

@testset "ADVUMPS.jl" begin
    @testset "hamiltonianmodels" begin
        println("hamiltonianmodels tests running...")
        include("hamiltonianmodels.jl")
    end

    @testset "example tensors" begin
        println("exampletensors tests running...")
        include("exampletensors.jl")
    end

    @testset "autodiff" begin
        println("autodiff tests running...")
        include("autodiff.jl")
    end

    @testset "fixedpoint" begin
        println("fixedpoint tests running...")
        include("fixedpoint.jl")
    end

    @testset "vumps" begin
        println("vumps tests running...")
        include("vumps.jl")
    end

    @testset "variationalipeps" begin
        println("variationalipeps tests running...")
        include("variationalipeps.jl")
    end
end


text:
# ## Heat

# ## Model

module heat #jl

using ReachabilityAnalysis, MAT, SparseArrays
using ReachabilityModels: @relpath

file = matopen(@relpath "heat.mat")

## system matrix
A = read(file, "A")

## input matrix
B = sparse([67], [1], [1.0], size(A, 1), 1)

## state domain
X = Universe(200)

## input domain
U = BallInf([0.0], 0.5)

function model(X0)
    S = @system(x' = Ax + Bu, x ∈ X, u ∈ U)
    return IVP(S, X0)
end

end  #jl

# ## References

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
# include("/Users/tarikz/Documents/Code/01-julia/05-PlotData.jl")

using CSV

#data = CSV.read("/Users/tarikz/Documents/Code/01-julia/05-NAEXKP01TRQ657S.csv")
data_lux = CSV.read("/Users/tarikz/Documents/Code/01-julia/05-b3019.csv", datarow = 2, delim = ',')
# Use 'showcols(data_lux)' to see the variable names

### Method 1: Using Plots
#using Plots
#plotly()
#plot( data.DATE, data.NAEXKP01TRQ657S, linewidth=2, title="GDP Turkey"  )

### Method 2: Using Gadfly (visit: http://gadflyjl.org)
using Gadfly
#draw(SVG("/Users/tarikz/Documents/Code/01-julia/output.svg", 6inch, 3inch), plot(data.DATE, data.NAEXKP01TRQ657S))
#plot(data.DATE, data.NAEXKP01TRQ657S)

### Reorganise the data (STATEC organises time series from latest to oldest)

#│ Row │ variable                                              │ eltype   │ nmissing │ first          │ last         │
#│     │ Symbol                                                │ DataType │ Int64    │ Any            │ Any          │
#├─────┼───────────────────────────────────────────────────────┼──────────┼──────────┼────────────────┼──────────────┤
#│ 1   │ Spécification                                         │ String   │ 0        │ septembre 2018 │ janvier 1995 │
#│ 2   │ 1. Salariés résidents sortants - données cvs          │ Float64  │ 0        │ 12551.5        │ 8628.11      │
#│ 3   │ 2. Salariés frontaliers entrants - données cvs        │ Float64  │ 0        │ 1.93647e5      │ 53773.7      │
#│ 4   │ 3. Emploi salarié national - données cvs              │ Float64  │ 0        │ 2.43207e5      │ 1.49295e5    │
#│ 5   │ 4. Emploi salarié intérieur - données cvs (3 + 2 - 1) │ Float64  │ 0        │ 4.24302e5      │ 1.9444e5     │
#│ 6   │ 5. Emploi non-salarié national - données cvs          │ Float64  │ 0        │ 22144.2        │ 16912.3      │
#│ 7   │ 6. Emploi non-salarié intérieur - données cvs         │ Float64  │ 0        │ 27091.1        │ 17462.6      │
#│ 8   │ 7. Emploi total national - données cvs (3 + 5)        │ Float64  │ 0        │ 2.65351e5      │ 1.66207e5    │
#│ 9   │ 8. Emploi total intérieur - données cvs (4 + 6)       │ Float64  │ 0        │ 4.51394e5      │ 2.11903e5    │
#│ 10  │ 9. Nombre de chômeurs - données cvs                   │ Float64  │ 0        │ 15027.0        │ 4331.61      │
#│ 11  │ 10. Population active - données cvs (7 + 9)           │ Float64  │ 0        │ 2.80378e5      │ 1.70538e5    │
#│ 12  │ 11. Taux de chômage (en %) - données cvs (9 / 10)     │ Float64  │ 0        │ 5.35954        │ 2.53996      │

nrows, ncols = size(data_lux)
df = zeros(nrows, 1)
for r in 0:nrows-1
     df[r+1, 1] = data_lux[nrows-r, 12]
end
# Working with dates: https://en.wikibooks.org/wiki/Introducing_Julia/Working_with_dates_and_times
dr = Dates.Date(1995,1):Dates.Month(1):Dates.Date(2018,9)

using Gadfly, Compose
# Visit: http://gadflyjl.org
plot(dr, df, label="empl")

max_stars_repo_name: tarikgit/julia-lux

max_issues_repo_name: tarikgit/julia-lux

max_forks_repo_name: tarikgit/julia-lux


text:
using QuantumESPRESSOExpress
using Test

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


text:
# explore warping functions.

include("../src/MonotoneMaps.jl")
import .MonotoneMaps

using FFTW
import PyPlot
import BSON
import Optim
import Random
using LinearAlgebra

import Interpolations

PyPlot.close("all")
fig_num = 1

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

Random.seed!(25)

## gives a random map.
D = 1 #2 # number of regions.
lb = 0.0
ub = 1.0
scale = 1.0 #0.5

# amount of input region used to map to z_lens.
input_range_percentage = 0.9
c = input_range_percentage*(ub-lb)

# generate random boundary points that define the regions.

# # for lb = -1.0
# z_st = [ -0.12; ]
# z_fin = [ 0.76; ]

# for lb = 0.0
# z_st = [ 0.12; ]
# z_fin = [ 0.76; ]

# z_st = [ 0.1; ]
# z_fin = [ 0.3; ]
#
#
# z_st = [ 0.6; ]
# z_fin = [ 0.8; ]

# z_st = [ -0.999999999999; ]
# z_fin = [ -0.9; ]
#
# z_st = [ 0.9; ]
# z_fin = [ 0.999; ]

xs, ys, ms, bs, len_s, len_z = MonotoneMaps.getpiecewiselines(z_st, z_fin, c; lb = lb, ub = ub)
f = xx->MonotoneMaps.evalpiecewise2Dlinearfunc(xx, xs, ys, ms, bs)*scale
finv = xx->MonotoneMaps.evalinversepiecewise2Dlinearfunc(xx/scale, xs, ys, ms, bs)

x_range = LinRange(lb, ub, 5000)
f_x = f.(x_range)
finv_y = finv.(f_x)

sanity_check = norm(sort(f_x)-f_x)
println("f_x: sanity_check = ", sanity_check)
@assert sanity_check < 1e-12

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(x_range, f_x)

PyPlot.xlabel("x")
PyPlot.ylabel("f")
PyPlot.title("target warp func")


text:
function create(c::Client, ::Type{Role}, g::AbstractGuild; kwargs...)
    return create_guild_role(c, g.id; kwargs...)
end

retrieve(c::Client, ::Type{Role}, g::AbstractGuild) = get_guild_roles(c, g.id)

function update(c::Client, r::Role, g::AbstractGuild; kwargs...)
    return modify_guild_role(c, g.id, r.id; kwargs...)
end

delete(c::Client, r::Role, g::AbstractGuild) = delete_guild_role(c, g.id, r.id)


text:
module ModelFitting
using ..Concepts
using ..BetterMGF
using ..Estimator
import LinearAlgebra

export DistributionDeduction
export ObservedOrMissing
export MatrixCompletionModel

struct DistributionDeduction <: Concepts.AbstractModelView end
struct ObservedOrMissing     <: Concepts.AbstractModelView end

struct MatrixCompletionModel end

@overload
function Concepts.check(object::Union{Type{Val{:continuous}},Continuous},
                        data::Array{T}) where T<:Real
    dist_to_integral_cast = LinearAlgebra.norm(data .- round.(data),2)
    return dist_to_integral_cast > 1e-2
end

@overload
function Concepts.check(object::Union{Type{Val{:integral}},Categorical},
                        data::Array{T}) where T<:Real
    dist_to_integral_cast = LinearAlgebra.norm(data .- round.(data),2)
    return dist_to_integral_cast <= 1e-4
end

@overload
function Concepts.check(object::Union{Type{Binary},
                                      Type{Val{:binary}},
                                      Binary},
                        data::Array{T}) where T<:Real
    return length(unique(data)) == 2
end

# const Concepts.check(object::Symbol,data::Array{T}) where T<:Real = Concepts.check(Val{object},data)

# @overload
# const Concepts.check(object::Symbol,arg...;kwargs...) = Concepts.check(Val{object},arg...;kwargs...)

@overload
function Concepts.check(object::Union{Type{Support},
                                      Type{Val{:support}},
                                      Support},
                        data::Array{T};
                        layout::Symbol = :flatten) where T<:Real
    if layout == :flatten
        return (minimum(data[:]), maximum(data[:]))
    end
    # TODO: implement for other data layouts
end

#const Concepts.check(object::Symbol,data::Array{T;layout::Symbol =:flatten) where T<:Real  =
#    Concepts.check(Val{object},data,layout=layout)

# for now we only compare with Gamma. More smooth distributions will be added later.
@overload
function Concepts.check(object::Union{Type{Val{:Gaussian}},
                                      Type{AbstractGaussian},
                                      AbstractGaussian},
                        data::Array{T}) where T<:Real
    if check(:continuous,data) == false
        return false
    end
    # check against gamma when we get positive support
    if check(:support,data) > (0,0) && choose(Gaussian,Gamma) == :Gamma
        return false
    end
    return true
end

@overload
function Concepts.check(object::Union{Type{Val{:Gamma}},
                                      Type{AbstractGamma},
                                      AbstractGamma},
                        data::Array{T}) where T<:Real
    if check(:continuous,data) == false
        return false
    end
    if check(:support,data) < (0,0)
        return false
    end
    if choose(Gaussian,Gamma) == :Gaussian
        return false
    end
    return true
end

# For now we let non binary categorical array be Poisson. We will refine this later.
@overload
function Concepts.check(object::Union{Type{Val{:Poisson}},
                                      Type{AbstractPoisson},
                                      AbstractPoisson},
                        data::Array{T}) where T<:Real
    if check(:integral,data) == false
        return false
    end
    if check(:binary,data) == true
        return false
    end
    return true
end

@overload
function Concepts.check(object::Union{Type{Val{:Bernoulli}},
                                      Type{AbstractBernoulli},
                                      AbstractBernoulli},
                        data::Array{T}) where T<:Real
    if check(:integral,data) == false
        return false
    end
    if check(:binary,data) == false
        return false
    end
    return true
end

@overload
function Concepts.check(object::Union{Type{Val{:NegativeBinomial}},
                                      Type{AbstractNegativeBinomial},
                                      AbstractNegativeBinomial},
                        data::Array{T}) where T<:Real
    return false #TODO: implement this properly
end

@overload
function Concepts.choose(a::Union{AbstractPoisson,
                                  Type{AbstractPoisson},
                                  Type{Val{:Poisson}}},
                         b::Union{AbstractNegativeBinomial,
                                  Type{AbstractNegativeBinomial},
                                  Type{Val{:NegativeBinomial}}};
                         data::AutoboxedArray{Real} = nothing,
                         comp::Comparator = Comparator{MGF}(MGF()))
    # for now we return poisson all the time
    return :Poisson
end

# we choose to return symbol for maximum flexibility
@overload
function Concepts.choose(a::Union{AbstractGaussian,
                                  Type{AbstractGaussian},
                                  Type{Val{:Gaussian}}},
                         b::Union{AbstractGamma,
                                  Type{AbstractGamma},
                                  Type{Val{:Gamma}}};
                         data::AutoboxedArray{Real} = nothing,
                         comp::Comparator = Comparator{MGF}(MGF()))
    if check(:support,data) <= (0,)
        @info "Data has negative support, force cast to Gaussian"
        return :Gaussian
    end
    mle_est_gaussian = estimator(MLE{AbstractGaussian}(),data)
    mle_est_gamma    = estimator(MLE{AbstractGamma}(),data)
    t = nothing
    if isnothing(comp.field[:eval_at])
        if mle_est_gamma[:α] <= 15 && mle_est_gamma[:θ] <= 15
            t = collect(0.01:0.001:0.02)
        else
            t = collect(0:0.0001:0.001)
        end
    else
        t = comp.field[:eval_at]
    end
    empirical_mgf = evaluate(SampleMGF(),t, data=data, order=20)
    gaussian_mgf  = evaluate(MGF(:Gaussian),t,μ = mle_est_gaussian[:μ],σ=mle_est_gaussian[:σ])
    gamma_mgf     = evaluate(MGF(:Gamma),t,α=mle_est_gamma[:α],θ=mle_est_gamma[:θ])
    if check(:l2diff,empirical_mgf, gaussian_mgf) < check(:l2diff,empirical_mgf,gamma_mgf)
        return :Gaussian
    end
    return :Gamma
end

@overload
function Concepts.fit(model::DistributionDeduction, data = Array{MaybeMissing{T}, 2};
                      data_layout::Symbol = :by_col) where T<: Number
    distribution_view = Array{Symbol, 2}(undef, size(data))
    if data_layout == :by_col
        for i = 1:size(data)[2]
            current_column = convert(Array{Float64}, filter(x -> !ismissing(x), data[:, i]))
            # @show(typeof(current_column))
            if check(:continuous, current_column)
                distribution_view[:, i] .= choose(:Gaussian, :Gamma, data = current_column)
                # display(distribution_view[:, i])
            end
            if check(:integral, current_column)
                # @show(i)
                if check(:binary, current_column)
                    distribution_view[:, i] .= :Bernoulli
                end
                if check(:Poisson, current_column)
                    distribution_view[:, i] .= :Poisson
                end
            end
        end
    end

    if data_layout == :by_row
        # TODO
    end
    # @show(distribution_view)
    return distribution_view
end

@overload
function Concepts.fit(model::ObservedOrMissing, data = Array{MaybeMissing{T}, 2})
    data_view = Array{Symbol, 2}(undef, size(data))
    for r in 1:size(data)[1]
        for c in 1:size(data)[2]
            if ismissing(data[r, c])
                data_view[r, c] = :Missing
            else
                data_view[r, c] = :Observed
            end
        end
    end
    return data_view
end

@overload
function Concepts.predict(model::MatrixCompletionModel;
                          completed_matrix,
                          type_tracker,
                          estimators = nothing)
    predicted_matrix = similar(completed_matrix)
    for dist in setdiff(keys(type_tracker.indices), [:Missing, :Observed])
        # idx = type_tracker[convert(Symbol, dist)]
        idx = type_tracker[type_conversion(Symbol, dist)]
        # if convert(Symbol, dist) != :NegativeBinomial
        if type_conversion(Symbol, dist) != :NegativeBinomial
            predicted_matrix[idx] .= predict(dist, forward_map(Val{dist}, completed_matrix[idx]))
        else
            predicted_matrix[idx] .= predict(dist, forward_map(Val{dist}, completed_matrix[idx],
                                                               r_estimate = estimators[:NegativeBinomial][:r]))
        end
    end
    return predicted_matrix
end

@overload
function Base.summary(model::MatrixCompletionModel;
                      predicted_matrix,
                      truth_matrix,
                      type_tracker,
                      tracker)
    ret = Dict()
    for dist in setdiff(keys(type_tracker.indices), [:Missing, :Observed])
        summary_missing_only = provide(Diagnostics{Any}(),
                                       reference  = truth_matrix[tracker[type_conversion(Symbol, dist)][:Missing]],
                                       input_data = predicted_matrix[tracker[type_conversion(Symbol, dist)][:Missing]])
        summary_all  = provide(Diagnostics{Any}(),
                               reference  = truth_matrix[type_tracker[type_conversion(Symbol, dist)]],
                               input_data = predicted_matrix[type_tracker[type_conversion(Symbol, dist)]])
        ret[dist] = Dict()
        ret[dist][:MissingOnly] = summary_missing_only
        ret[dist][:All]         = summary_all
    end
    return ret
end

end


text:
# Input and output routines

"Maximum polynomial degree of Mineos polynomial model files"
MINEOS_MAX_DEGREE = 4

"""
    read_mineos(file) -> m::LinearLayeredModel

Read a 1D seismic model from a Mineos tabular format file.

### References:
- https://geodynamics.org/cig/software/mineos/mineos-manual.pdf
"""
function read_mineos(io::IO)
    name = readline(io)
    ifanis, tref, ifdeck = parse.((Int, Float64, Int), split(readline(io)))
    aniso = ifanis == 1
    ifdeck != 1 && error("File '$file' is not in tabular format (`ifdeck` is `$ifdeck`)")
    n, nic, noc = parse.(Int, split(readline(io)))
    d = readdlm(io)
    r = d[:,1]./1e3
    a = maximum(r)
    rho = d[:,2]./1e3
    if aniso
        vpv = d[:,3]./1e3
        vsv = d[:,4]./1e3
        vph = d[:,7]./1e3
        vsh = d[:,8]./1e3
        eta = d[:,9]
        voigt_vels = voigt_velocities.(vpv, vsv, vph, vsh, eta)
        vp = [v[1] for v in voigt_vels]
        vs = [v[2] for v in voigt_vels]
    else
        vp = d[:,3]./1e3
        vs = d[:,4]./1e3
        vpv = vph = vsv = vsh = eta = []
    end
    Qκ = d[:,5]
    Qμ = d[:,6]
    attenuation = !(all(isequal(0), Qκ) && all(isequal(0), Qμ))
    if !attenuation
        Qμ = []
        Qκ = []
    end
    LinearLayeredModel(a, n, r, vp, vs, rho, aniso, vph, vpv, vsh, vsv, eta,
                       attenuation, Qμ, Qκ)
end

read_mineos(file::AbstractString) = (open(file, "r") do io; read_mineos(io); end)

# TODO: Remove when bumping minor or major version
Base.@deprecate write_mineos(model::LinearLayeredModel, file::AbstractString, args...) write_mineos(file, model, args...)

"""
    write_mineos(file, m::LinearLayeredModel, freq=1.0, title="Model from SeisModels.jl")
    write_mineos(io::IO, m, args...)

Save a `SeisModel1D` as a Mineos 'tabular' format file to `file` on disk or
to an `io`.  Optionally supply the reference
frequency `freq` in Hz (which defaults to 1 Hz) and a `title`.

Note that Mineos has two types of model file; one parameterised by radial knots
with values interpolated linearly between them ('tabular' format),
and one with PREM-style polynomials ('polynomial' format).
Mineos uses the tabular format internally even if a polynomial file is read in,
and testing suggests that polynomial files are not used correctly by
the [CIG](https://geodynamics.org/cig/software/mineos/) version of Mineos,
so it is not recommended to use them at all.  Consequently, SeisModels does
not support writing polynomial format files.

## Reference
- Mineos manual, https://geodynamics.org/cig/software/mineos/mineos-manual.pdf
"""
function write_mineos(io::IO, m::LinearLayeredModel, freq=1.0, title="Model from SeisModels.jl")
    _check_mineos_title(title)
    ifanis = isanisotropic(m) ? 1 : 0
    tref = freq
    ifdeck = 1
    N = _check_mineos_num_layers(m)
    nic, noc = core_interface_layers(m)
    rho, vp, vs = m.density.*1e3, m.vp.*1e3, m.vs.*1e3
    if m.aniso
        vpv, vph, vsv, vsh, eta = m.vpv.*1e3, m.vph.*1e3, m.vsv.*1e3, m.vsh.*1e3, m.eta
    else
        vpv = vph = vp
        vsh = vsv = vs
        eta = ones(m.n)
    end
    if hasattenuation(m)
        Qμ, Qκ = m.Qμ, m.Qκ
    else
        Qμ = Qκ = zeros(m.n)
    end
    println(io, title)
    println(io, ifanis, " ", tref, " ", ifdeck)
    println(io, N, " ", nic, " ", noc)
    writedlm(io, [m.r*1e3 rho vpv vsv Qκ Qμ vph vsh eta])
    nothing
end

write_mineos(file::AbstractString, m::LinearLayeredModel, args...) =
    (open(file, "w") do io; write_mineos(io, m, args...); end)

"""
    _check_mineos_title(title)

Issue a warning if `title` is longer than Mineos's maximum title length
(80 characters).
"""
_check_mineos_title(title) = length(title) > 80 &&
        @warn("Mineos model files can have titles only 80 characters long " *
              "('$title' is $(length(title)) characters)")

"""
    _check_mineos_num_layers(m::SeisModel) -> N

Issue a warning if the model `m` contains more than Mineos's maximum
number of layers (350), and return the number of layers `N`.
"""
_check_mineos_num_layers(m::SeisModel) = (N = m.n;
    N > 350 && @warn("Mineos model files are limited to N ≤ 350 (have $N layers)"); N)

"""
    _check_mineos_max_poly_degree(m::PREMPolyModel)

Throw an error if the model `m` has any fields with polynomial degree greater
than Mineos's maximum degree (4).
"""
function _check_mineos_max_poly_degree(m::PREMPolyModel)
    # Mineos only deals with up to degree-four polynomials with 5 coefficients
    max_degree = 4
    for field in (:vp, :vs, :density, :vph, :vpv, :vsh, :vsv, :eta, :Qμ, :Qκ)
        v = getfield(m, field)
        if ndims(v) > 1 && size(v, 1) > MINEOS_MAX_DEGREE + 1
            throw(ArgumentError("maximum degree polynomial in Mineos polynomial files is 4"))
        end
    end
    nothing
end

"""
    read_tvel(file) -> model::LinearLayeredModel
    read_tvel(io::IO) -> model::LinearLayeredModel

Read a `model` from a file in the ttimes 'tvel' format.  This format
specifies a `LinearLayeredModel` with only isotropic Vp, Vs and
density.

The two header lines are ignored but must be present.
"""
function read_tvel(io::IO; file=nothing)
    lines = readlines(io)
    length(lines) > 2 || error("input .tvel string has too few lines")
    r, Vp, Vs, rho = Float64[], Float64[], Float64[], Float64[]
    for i in 3:length(lines)
        tokens = split(lines[i])
        length(tokens) >= 4 || error("too few columns at line $i of " *
            (file === nothing ? "input" : "file '$file'"))
        rr, vvpp, vvss, ρρ = parse.(Float64, tokens)
        push!(r, rr)
        push!(Vp, vvpp)
        push!(Vs, vvss)
        push!(rho, ρρ)
    end
    R = r[end]
    r .= R .- r
    reverse!(r)
    reverse!(Vp)
    reverse!(Vs)
    reverse!(rho)
    LinearLayeredModel(r=r, vp=Vp, vs=Vs, density=rho)
end

read_tvel(file::AbstractString) = (open(file, "r") do io
    read_tvel(io, file=file); end)

"""
    write_tvel(file, model::LinearLayeredModel; comment1="SeisModels.jl - P", comment2="SeisModels.jl - S")
    write_tvel(io::IO, model; kwargs...)

Write a `model` to `file` or an arbitrary `IO` object
in the `ttimes` 'tvel' format.

Optionally supply the two comment lines, `comment1` and `comment2`
which form the header of the file.
"""
function write_tvel(io::IO, model::LinearLayeredModel; file=nothing,
    comment1="SeisModels.jl - P", comment2="SeisModels.jl - S")
    println(io, comment1)
    println(io, comment2)
    dep = surface_radius(model) .- reverse(model.r)
    writedlm(io, [dep reverse(model.vp) reverse(model.vs) reverse(model.density)])
end

write_tvel(file::AbstractString, model::LinearLayeredModel; kwargs...) =
    (open(file, "w") do io; write_tvel(io, model; file=file, kwargs...); end)

"""
    core_interface_layers(m) -> i_icb, i_cmb

Return the indices of the inner-core boundary, `i_icb`, and the core-mantle boundary,
`i_cmb`, for the model `m` based on the following assumptions:

1. The outer core is a liquid and has Vs = 0.
2. There is only one outer core
3. The inner core is solid and has Vs > 0.
4. There is only one inner core
5. The mantle above the outer core is solid and has Vs > 0.

These numbers refer to the layer which is the **top** of the inner/outer core.
"""
function core_interface_layers(m)
    i = 1
    i_icb = i_cmb = 0
    got_icb = false
    while i <= m.n
        if !got_icb && isapprox(m.vs[i], 0, atol=eps(eltype(m.vs)))
            got_icb = true
            i > 1 || error("Unexpectedly found Vs = 0 at centre of model")
            i_icb = i - 1
        end
        if got_icb && m.vs[i] > eps(eltype(m.vs))
            i_cmb = i - 1
            return i_icb, i_cmb
        end
        i += 1
    end
end

"""
    core_interface_layers(m::PREMPolyModel) -> i_icb, i_cmb

Return the ICB and CMB indices for a `PREMPolyModel`, using the following
heuristic:

1. The ICB and CMB are located at boundaries between layers.
2. The velocity 1 m above the inner core is liquid and has Vs = 0.
3. There is a solid inner core at the centre of the model.
4. The mantle is above a liquid, above a solid, and is itself solid.
"""
function core_interface_layers(m::PREMPolyModel)
    i = 1
    i_icb = i_cmb = 0
    got_icb = false
    while i <= m.n
        # Radii one meter either side of a layer boundary
        r_plus_one = m.r[i] + 1e-3
        if !got_icb && isapprox(vs(m, r_plus_one), 0, atol=eps(eltype(m.vs)))
            got_icb = true
            i_icb = i
        end
        if got_icb && vs(m, r_plus_one) > eps(eltype(m.vs))
            i_cmb = i
            return i_icb, i_cmb
        end
        i += 1
    end
end