text:
using Test
using PlayingCards
using TexasHoldem: Player, Bot5050, TransactionManager, dealer_id, Table
TH = TexasHoldem

@testset "TransactionManagers - Lowest bank roll goes all-in and wins it all" begin
    table_cards = (A♢, K♢, Q♢, 2♠, 3♠)
    players = (
        Player(Bot5050(), 1, (A♠, A♣); bank_roll = 1*100),
        Player(Bot5050(), 2, (K♠, K♣); bank_roll = 2*100),
        Player(Bot5050(), 3, (Q♠, Q♣); bank_roll = 3*100),
    )
    tm = TH.TransactionManager(players)
    table = Table(;players=players,cards=table_cards,transactions=tm)
    @test TH.seat_number.(tm.side_pots) == [1,2,3]

    TH.raise_to!(table, players[1], 100) # raise all-in
    TH.call!(table, players[2]) # call
    TH.call!(table, players[3]) # call

    @test TH.amount.(tm.side_pots) == [300.0, 0.0, 0.0]
    TH.distribute_winnings!(players, tm, table_cards)
    @test TH.amount.(tm.side_pots) == [0.0, 0.0, 0.0]

    @test bank_roll(players[1]) == 300
    @test bank_roll(players[2]) == 100
    @test bank_roll(players[3]) == 200
end

@testset "TransactionManagers - Highest bank roll goes all-in and wins it all" begin
    table_cards = (A♢, K♢, Q♢, 2♠, 3♠)
    players = (
        Player(Bot5050(), 1, (A♠, A♣); bank_roll = 3*100),
        Player(Bot5050(), 2, (K♠, K♣); bank_roll = 2*100),
        Player(Bot5050(), 3, (Q♠, Q♣); bank_roll = 1*100),
    )
    tm = TH.TransactionManager(players)
    table = Table(;players=players,cards=table_cards,transactions=tm)

    TH.raise_to!(table, players[1], 100) # Raise
    TH.call!(table, players[2]) # call
    TH.call!(table, players[3]) # all-in

    @test TH.amount.(tm.side_pots) == [300.0, 0.0, 0.0]
    TH.distribute_winnings!(players, tm, table_cards)
    @test TH.amount.(tm.side_pots) == [0.0, 0.0, 0.0]

    @test bank_roll(players[1]) == 500
    @test bank_roll(players[2]) == 100
    @test bank_roll(players[3]) == 0
end

@testset "TransactionManagers - Lowest bank roll goes all-in and wins a split pot" begin
    table_cards = (A♢, K♢, Q♢, 2♠, 3♠)
    players = (
        Player(Bot5050(), 1, (A♠, A♣); bank_roll = 1*100),
        Player(Bot5050(), 2, (K♠, K♣); bank_roll = 2*100),
        Player(Bot5050(), 3, (Q♠, Q♣); bank_roll = 3*100),
    )
    tm = TH.TransactionManager(players)
    table = Table(;players=players,cards=table_cards,transactions=tm)

    TH.raise_to!(table, players[1], 100) # Raise all-in
    TH.call!(table, players[2]) # call
    TH.call!(table, players[3]) # call

    @test TH.amount.(tm.side_pots) == [300.0, 0.0, 0.0]

    @test_throws AssertionError TH.call!(table, players[1]) # already all-in!

    TH.reset_round!(table)

    TH.raise_to!(table, players[2], 100) # Raise all-in
    TH.call!(table, players[3]) # call

    @test TH.amount.(tm.side_pots) == [300.0, 200.0, 0.0]
    TH.distribute_winnings!(players, tm, table_cards)
    @test TH.amount.(tm.side_pots) == [0.0, 0.0, 0.0]

    @test bank_roll(players[1]) == 300
    @test bank_roll(players[2]) == 200
    @test bank_roll(players[3]) == 100
end

@testset "TransactionManagers - Highest bank roll goes all-in and wins it all" begin
    table_cards = (A♢, K♢, Q♢, 2♠, 3♠)
    players = (
        Player(Bot5050(), 1, (A♠, A♣); bank_roll = 3*100),
        Player(Bot5050(), 2, (K♠, K♣); bank_roll = 2*100),
        Player(Bot5050(), 3, (Q♠, Q♣); bank_roll = 1*100),
    )
    tm = TH.TransactionManager(players)
    table = Table(;players=players,cards=table_cards,transactions=tm)

    TH.raise_to!(table, players[1], 100) # Raise
    TH.call!(table, players[2]) # call
    TH.call!(table, players[3]) # all-in

    @test TH.amount.(tm.side_pots) == [300.0, 0.0, 0.0]

    TH.reset_round!(table)

    TH.raise_to!(table, players[1], 100) # call
    TH.call!(table, players[2]) # all-in
    @test_throws AssertionError TH.call!(table, players[3]) # already all-in!

    @test TH.amount.(tm.side_pots) == [300.0, 200.0, 0.0]

    TH.distribute_winnings!(players, tm, table_cards)
    @test TH.amount.(tm.side_pots) == [0.0, 0.0, 0.0]

    @test bank_roll(players[1]) == 600
    @test bank_roll(players[2]) == 0
    @test bank_roll(players[3]) == 0
end

@testset "TransactionManagers - Semi-complicated split pot (shared winners)" begin
    table_cards = (T♢, Q♢, A♠, 8♠, 9♠)
    players = (
        Player(Bot5050(), 1, (4♠, 5♣); bank_roll = 1*100), # bust
        Player(Bot5050(), 2, (K♠, K♣); bank_roll = 2*100), # win, split with player 3
        Player(Bot5050(), 3, (K♡,K♢); bank_roll = 3*100), # win, split with player 2
        Player(Bot5050(), 4, (2♡, 3♢); bank_roll = 4*100), # bust
        Player(Bot5050(), 5, (7♠, 7♣); bank_roll = 5*100), # 2nd to players 2 and 3, win remaining pot
        Player(Bot5050(), 6, (2♠, 3♣); bank_roll = 6*100), # lose, but not bust
    )
    tm = TH.TransactionManager(players)
    table = Table(;players=players,cards=table_cards,transactions=tm)

    TH.raise_to!(table, players[1], 100) # raise all-in
    TH.call!(table, players[2]) # call
    TH.call!(table, players[3]) # call
    TH.call!(table, players[4]) # call
    TH.call!(table, players[5]) # call
    TH.call!(table, players[6]) # call
    @test TH.amount.(tm.side_pots) == [600.0, 0.0, 0.0, 0.0, 0.0, 0.0]

    TH.reset_round!(table)

    TH.raise_to!(table, players[2], 100) # raise all-in
    TH.call!(table, players[3]) # call
    TH.call!(table, players[4]) # call
    TH.call!(table, players[5]) # call
    TH.call!(table, players[6]) # call
    @test TH.amount.(tm.side_pots) == [600.0, 500.0, 0.0, 0.0, 0.0, 0.0]

    TH.reset_round!(table)

    TH.raise_to!(table, players[3], 100) # raise all-in
    TH.call!(table, players[4]) # call
    TH.call!(table, players[5]) # call
    TH.call!(table, players[6]) # call
    @test TH.amount.(tm.side_pots) == [600.0, 500.0, 400.0, 0.0, 0.0, 0.0]

    TH.reset_round!(table)

    TH.raise_to!(table, players[4], 100) # raise all-in
    TH.call!(table, players[5]) # call
    TH.call!(table, players[6]) # call
    @test TH.amount.(tm.side_pots) == [600.0, 500.0, 400.0, 300.0, 0.0, 0.0]

    TH.reset_round!(table)

    TH.raise_to!(table, players[5], 100) # raise all-in
    TH.call!(table, players[6]) # call
    @test TH.amount.(tm.side_pots) ≈ [600.0, 500.0, 400.0, 300.0, 200.0, 0.0]

    TH.distribute_winnings!(players, tm, table_cards)
    @test TH.amount.(tm.side_pots) == [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

    @test bank_roll(players[1]) == 0 # bust
    @test bank_roll(players[2]) == 550 # = 600/2+500/2
    @test bank_roll(players[3]) == 950 # = 600/2+500/2+400
    @test bank_roll(players[4]) == 0 # bust
    @test bank_roll(players[5]) == 500 # all contributions after 3rd all-in (3*100+2*100)
    @test bank_roll(players[6]) == 100 # lost (but not all-in)
    @test sum(bank_roll.(players)) == 2100
end

@testset "TransactionManagers - Single round split pot (shared winners), with simple re-raises" begin
    table_cards = (T♢, Q♢, A♠, 8♠, 9♠)
    players = (
        Player(Bot5050(), 1, (4♠, 5♣); bank_roll = 1*100), # bust
        Player(Bot5050(), 2, (K♠, K♣); bank_roll = 2*100), # win, split with player 3
        Player(Bot5050(), 3, (K♡,K♢); bank_roll = 3*100), # win, split with player 2
        Player(Bot5050(), 4, (2♡, 3♢); bank_roll = 4*100), # bust
        Player(Bot5050(), 5, (7♠, 7♣); bank_roll = 5*100), # 2nd to players 2 and 3, win remaining pot
        Player(Bot5050(), 6, (2♠, 3♣); bank_roll = 6*100), # lose, but not bust
    )
    tm = TH.TransactionManager(players)
    table = Table(;players=players,cards=table_cards,transactions=tm)
    @test TH.amount.(tm.side_pots) == [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

    TH.raise_to!(table, players[1], 100) # raise all-in
    @test TH.amount.(tm.side_pots) == [100.0, 0.0, 0.0, 0.0, 0.0, 0.0]

    TH.raise_to!(table, players[2], 200) # raise all-in
    @test TH.amount.(tm.side_pots) == [200.0, 100.0, 0.0, 0.0, 0.0, 0.0]

    TH.raise_to!(table, players[3], 300) # raise all-in
    @test TH.amount.(tm.side_pots) == [300.0, 200.0, 100.0, 0.0, 0.0, 0.0]

    TH.raise_to!(table, players[4], 400) # raise all-in
    @test TH.amount.(tm.side_pots) == [400.0, 300.0, 200.0, 100.0, 0.0, 0.0]

    TH.raise_to!(table, players[5], 500) # raise all-in
    @test TH.amount.(tm.side_pots) == [500.0, 400.0, 300.0, 200.0, 100.0, 0.0]

    TH.call!(table, players[6]) # call
    @test TH.amount.(tm.side_pots) ≈ [600.0, 500.0, 400.0, 300.0, 200.0, 0.0]

    TH.distribute_winnings!(players, tm, table_cards)
    @test TH.amount.(tm.side_pots) == [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

    @test bank_roll(players[1]) == 0 # bust
    @test bank_roll(players[2]) == 550 # = 600/2+500/2
    @test bank_roll(players[3]) == 950 # = 600/2+500/2+400
    @test bank_roll(players[4]) == 0 # bust
    @test bank_roll(players[5]) == 500 # all contributions after 3rd all-in (3*100+2*100)
    @test bank_roll(players[6]) == 100 # lost (but not all-in)
    @test sum(bank_roll.(players)) == 2100
end

@testset "TransactionManagers - Single round split pot (shared winners), with simple re-raises, reversed bank roll order" begin
    table_cards = (T♢, Q♢, A♠, 8♠, 9♠)
    players = (
        Player(Bot5050(), 1, (2♠, 3♣); bank_roll = 6*100), # lose, but not bust
        Player(Bot5050(), 2, (7♠, 7♣); bank_roll = 5*100), # 2nd to players 2 and 3, win remaining pot
        Player(Bot5050(), 3, (2♡, 3♢); bank_roll = 4*100), # bust
        Player(Bot5050(), 4, (K♡,K♢); bank_roll = 3*100), # win, split with player 2
        Player(Bot5050(), 5, (K♠, K♣); bank_roll = 2*100), # win, split with player 3
        Player(Bot5050(), 6, (4♠, 5♣); bank_roll = 1*100), # bust
    )
    tm = TH.TransactionManager(players)
    table = Table(;players=players,cards=table_cards,transactions=tm)
    @test TH.amount.(tm.side_pots) == [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

    TH.raise_to!(table, players[1], 500) # raise to 500
    @test TH.amount.(tm.side_pots) == [100.0, 100.0, 100.0, 100.0, 100.0, 0.0]

    TH.call!(table, players[2]) # call
    @test TH.amount.(tm.side_pots) == [200.0, 200.0, 200.0, 200.0, 200.0, 0.0]

    TH.call!(table, players[3]) # call
    @test TH.amount.(tm.side_pots) == [300.0, 300.0, 300.0, 300.0, 200.0, 0.0]

    TH.call!(table, players[4]) # call
    @test TH.amount.(tm.side_pots) == [400.0, 400.0, 400.0, 300.0, 200.0, 0.0]

    TH.call!(table, players[5]) # call
    @test TH.amount.(tm.side_pots) == [500.0, 500.0, 400.0, 300.0, 200.0, 0.0]

    TH.call!(table, players[6]) # call
    @test TH.amount.(tm.side_pots) == [600.0, 500.0, 400.0, 300.0, 200.0, 0.0]

    TH.distribute_winnings!(players, tm, table_cards)
    @test TH.amount.(tm.side_pots) == [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

    @test bank_roll(players[1]) == 100 # lost (but not all-in)
    @test bank_roll(players[2]) == 500 # all contributions after 3rd all-in (3*100+2*100)
    @test bank_roll(players[3]) == 0 # bust
    @test bank_roll(players[4]) == 950 # = 600/2+500/2+400
    @test bank_roll(players[5]) == 550 # = 600/2+500/2
    @test bank_roll(players[6]) == 0 # bust
    @test sum(bank_roll.(players)) == 2100
end


text:
module QASM

export @qasm_str

using RBNF
using ExprTools
using OpenQASM
using OpenQASM.Types
using OpenQASM.Types: Gate
using ..YaoCompiler

mutable struct VirtualRegister
    type::Symbol
    address::UnitRange{Int}
end

mutable struct RegisterRecord
    map::Dict{String,VirtualRegister}
    nqubits::Int
    ncbits::Int
end

mutable struct GateRegisterRecord
    map::Dict
    total::Int
end

Base.getindex(x::RegisterRecord, key) = x.map[key]
Base.getindex(x::VirtualRegister, xs...) = x.address[xs...]

RegisterRecord() = RegisterRecord(Dict{String,UnitRange{Int}}(), 0, 0)
GateRegisterRecord() = GateRegisterRecord(Dict(), 0)

struct Ctx
    m::Module
    source::LineNumberNode
    record::Any
end

# tokens don't need context
transpile(::Ctx, x::RBNF.Token) = transpile(x)

transpile(x::RBNF.Token{:unnamed}) = Symbol(x.str)
transpile(x::RBNF.Token{:reserved}) = Symbol(x.str)

function transpile(x::RBNF.Token{:id})
    x.str == "pi" && return Base.pi
    return Symbol(x.str)
end

function transpile(x::RBNF.Token{:float64})
    return Base.parse(Float64, x.str)
end

function transpile(x::RBNF.Token{:int})
    return Base.parse(Int, x.str)
end

transpile_list(::Ctx, ::Nothing) = Any[]
transpile_list(ctx::Ctx, x) = Any[transpile(ctx, x)]

function transpile_list(ctx::Ctx, xs::Vector)
    [transpile(ctx, each) for each in xs]
end

transpile(m::Module, ast::MainProgram) = transpile(m, LineNumberNode(0), ast)

function transpile(m::Module, l::LineNumberNode, ast::MainProgram)
    # check sure minimum compatibility
    @assert v"2.0.0" <= ast.version < v"3.0.0"

    code = Expr(:block)
    body = Expr(:block)
    routines = []
    record = scan_registers(ast)
    ctx = Ctx(m, l, record)

    for stmt in ast.prog
        if stmt isa RegDecl
            continue
        elseif stmt isa Types.Gate
            push!(routines, transpile(ctx, stmt))
        elseif stmt isa Include
            file = stmt.file.str[2:end-1]
            # use relative path to current file if not in REPL
            # isinteractive can be true in IDEs
            if !isnothing(l.file) && !(isinteractive() && isempty(PROGRAM_FILE))
                file = joinpath(dirname(string(l.file)), file)
            end
            source = read(file, String)
            ast = OpenQASM.parse(source)
            push!(code.args, transpile(m, l, ast))
        else
            ex = transpile(ctx, stmt)
            if !isnothing(ex)
                push!(body.args, ex)
            end
        end
    end

    # if there are classical registers
    # return them in a NamedTuple
    ret = Expr(:tuple)
    for (k, r) in record.map
        if r.type === :classical
            name = Symbol(k)
            push!(ret.args, Expr(:(=), name, name))
        end
    end

    if !isempty(body.args)
        if isempty(ret.args)
            push!(body.args, :(return))
        else
            push!(body.args, :(return $ret))
        end
    end

    # routines
    for each in routines
        push!(code.args, each)
    end

    # create an anoymous routine
    # if there are global statements
    if !isempty(body.args)
        def = Dict{Symbol,Any}(:name => gensym(:qasm), :body => body)
        push!(code.args, YaoCompiler.device_def(def))
    end
    return code
end

function transpile(ctx::Ctx, stmt::Types.Gate)
    name = transpile(stmt.decl.name)
    args = transpile_list(ctx, stmt.decl.cargs)
    record = transpile_gate_registers(stmt.decl.qargs)
    body = Expr(:block)
    new_ctx = Ctx(ctx.m, ctx.source, record)

    for each in stmt.body
        push!(body.args, transpile(new_ctx, each))
    end

    def = Dict(:name => name, :args => args, :body => body)
    return YaoCompiler.device_def(def)
end

function transpile_gate_registers(stmt::Vector)
    record = GateRegisterRecord()
    for each in stmt
        haskey(record.map, each.str) && throw(Meta.ParseError("duplicated register name $(each.str)"))
        record.total += 1
        record.map[each.str] = record.total
    end
    return record
end

semantic_gate(gate, locs) = Expr(:call, GlobalRef(YaoCompiler.Semantic, :gate), gate, locs)
semantic_ctrl(gate, locs, ctrl) =
    Expr(:call, GlobalRef(YaoCompiler.Semantic, :ctrl), gate, locs, ctrl)

function transpile(ctx::Ctx, stmt::UGate)
    code = Expr(:block)
    locs = transpile(ctx, stmt.qarg)
    push!(
        code.args,
        semantic_gate(Expr(:call, GlobalRef(Intrinsics, :Rz), transpile(ctx, stmt.z1)), locs),
    )
    push!(
        code.args,
        semantic_gate(Expr(:call, GlobalRef(Intrinsics, :Ry), transpile(ctx, stmt.y)), locs),
    )
    push!(
        code.args,
        semantic_gate(Expr(:call, GlobalRef(Intrinsics, :Rz), transpile(ctx, stmt.z2)), locs),
    )
    return code
end

function transpile(ctx::Ctx, stmt::CXGate)
    return semantic_ctrl(
        GlobalRef(Intrinsics, :X),
        transpile(ctx, stmt.qarg),
        CtrlLocations(transpile(ctx, stmt.ctrl)),
    )
end

function transpile(ctx::Ctx, stmt::IfStmt)
    return :(
        if $(transpile(ctx, stmt.left)) == $(transpile(ctx, stmt.right))
            $(transpile(ctx, stmt.body))
        end
    )
end

function transpile(ctx::Ctx, stmt::Measure)
    locs = transpile(ctx, stmt.qarg)
    name = transpile(ctx, stmt.carg)
    return Expr(:(=), name, Expr(:call, GlobalRef(YaoCompiler.Semantic, :measure), locs))
end

function transpile(ctx::Ctx, stmt::Barrier)
    return Expr(
        :call,
        GlobalRef(YaoCompiler.Semantic, :barrier),
        transpile_locations(ctx, stmt.qargs),
    )
end

function transpile(ctx::Ctx, stmt::Instruction)
    op = stmt.name
    # NOTE: these are not intrinsic function in QASM
    # users need qelib1.inc to get the definition
    # but for convenience we treat them as intrinsic
    # function here in YaoCompiler, since they are predefined
    # as stdlib in YaoCompiler.

    # isnothing(stmt.lst1) || throw(Meta.ParseError("$op gate should not have classical parameters"))
    locs = transpile_locations(ctx, stmt.qargs)

    if op == "x"
        semantic_gate(Intrinsics.X, locs)
    elseif op == "y"
        semantic_gate(Intrinsics.Y, locs)
    elseif op == "z"
        semantic_gate(Intrinsics.Z, locs)
    elseif op == "h"
        semantic_gate(Intrinsics.H, locs)
    elseif op == "s"
        semantic_gate(Intrinsics.S, locs)
    elseif op == "ccx"
        semantic_ctrl(Intrinsics.X, locs[3], CtrlLocations(locs[1:2]))
    else # some user defined routine
        gate = Expr(:call, GlobalRef(ctx.m, Symbol(op)), transpile_list(ctx, stmt.cargs)...)
        semantic_gate(gate, locs)
    end
end

function transpile(ctx::Ctx, stmt::FnExp)
    return Expr(:call, stmt.fn, transpile(ctx, stmt.arg))
end

function transpile_locations(ctx, stmts::Vector)
    locs = map(stmts) do stmt
        transpile(ctx, stmt)
    end
    return merge_locations(locs...)
end

function transpile(ctx::Ctx, stmt::Bit)
    record = ctx.record
    if record isa RegisterRecord
        r = record[stmt.name.str]
        r.type === :classical && return Symbol(stmt.name.str)

        if isnothing(stmt.address)
            return Locations(r[:])
        else
            address = transpile(stmt.address)
            return Locations(r[address+1])
        end
    else
        return Locations(record.map[stmt.name.str])
    end
end

transpile(ctx::Ctx, stmt::Negative) = Expr(:call, -, transpile(ctx, stmt.value))

function transpile(ctx::Ctx, stmt::Tuple)
    length(stmt) == 3 || throw(Meta.ParseError("unrecognized expression: $stmt"))
    stmt[2]::RBNF.Token
    if stmt[2].str in ("+", "-", "*", "/")
        return Expr(:call, Symbol(stmt[2].str), transpile(ctx, stmt[1]), transpile(ctx, stmt[3]))
    else
        throw(Meta.ParseError("unrecognized expression: $stmt"))
    end
end

function scan_registers(ast::MainProgram)
    return scan_registers!(RegisterRecord(), ast)
end

function scan_registers!(record::RegisterRecord, ast::MainProgram)
    for stmt in ast.prog
        scan_registers!(record, stmt)
    end
    return record
end

function scan_registers!(record::RegisterRecord, ast::RegDecl)
    nbits = transpile(ast.size)
    nqubits = record.nqubits
    ncbits = record.ncbits

    if ast.type.str == "qreg"
        record.map[ast.name.str] = VirtualRegister(:quantum, (nqubits+1):(nqubits+nbits))
        record.nqubits += nbits
    else # classical
        record.map[ast.name.str] = VirtualRegister(:classical, (ncbits+1):(ncbits+nbits))
        record.ncbits += nbits
    end
    return record
end

scan_registers!(record::RegisterRecord, ast) = record

function qasm_str_m(m::Module, l::LineNumberNode, source::String)
    ast = OpenQASM.parse(source)
    return transpile(m, l, ast)
end

macro qasm_str(source::String)
    return esc(qasm_str_m(__module__, __source__, source))
end

macro qasm_str(source::Expr)
    source.head === :string || error("expect a String")

    args = map(source.args) do x
        x isa String && return x
        return Base.eval(__module__, x)
    end

    return esc(qasm_str_m(__module__, __source__, join(args)))
end

macro include_str(path)
    if path isa Expr
        path.head === :string || error("expect a String")
        file = map(path.args) do x
            x isa String && return x
            return Base.eval(__module__, x)
        end |> join
    elseif path isa String
        file = path
    else
        error("expect a String")
    end

    return esc(qasm_str_m(__module__, __source__, read(file, String)))
end

end # end module

ext: jl


text:
module Planets

using ..Ahorn, Maple

placements = Maple.Planets

function Ahorn.editingOptions(effect::Maple.Planets)
    return Dict{String, Any}(
        "size" => Maple.planet_effect_sizes
    )
end

function Ahorn.canFgBg(effect::Maple.Planets)
    return true, true
end

end


text:
function _degree_centrality(g::AGraphOrDiGraph, gtype::Integer; normalize=true)
   n_v = nv(g)
   c = zeros(n_v)
   for v in 1:n_v
       if gtype == 0    # count both in and out degree if appropriate
           deg = out_degree(g, v) + (is_directed(g) ? in_degree(g, v) : 0.0)
       elseif gtype == 1    # count only in degree
           deg = in_degree(g, v)
       else                 # count only out degree
           deg = out_degree(g, v)
       end
       s = normalize ? (1.0 / (n_v - 1.0)) : 1.0
       c[v] = deg*s
   end
   return c
end

# TODO avoid repetition of this docstring
"""Calculates the [degree centrality](https://en.wikipedia.org/wiki/Centrality#Degree_centrality)
of the graph `g`, with optional (default) normalization."""
degree_centrality(g::AGraphOrDiGraph; kws...) = _degree_centrality(g, 0; kws...)
"""Calculates the [degree centrality](https://en.wikipedia.org/wiki/Centrality#Degree_centrality)
of the graph `g`, with optional (default) normalization."""
in_degree_centrality(g::AGraphOrDiGraph; kws...) = _degree_centrality(g, 1; kws...)
"""Calculates the [degree centrality](https://en.wikipedia.org/wiki/Centrality#Degree_centrality)
of the graph `g`, with optional (default) normalization."""
out_degree_centrality(g::AGraphOrDiGraph; kws...) = _degree_centrality(g, 2; kws...)


text:
module ReferenceMicrocalFiles

const ljhdir = normpath(joinpath( dirname(@__FILE__), "../ljh"))

mutable struct ReferenceMicrocalFile
	filename::String
	noise_filename::AbstractString
	line_energies_ev::Vector{Float64}
	note::String
end
const dict = Dict()
function add_ref_file(key, fname, noisename, line_energies_ev, note)
	fname_full = joinpath(ljhdir, fname)
	noisename_full = joinpath(ljhdir, noisename)
	noisename_full  = isfile(noisename_full) ? noisename_full : ""
	dict[key] = ReferenceMicrocalFile(fname_full, noisename_full, line_energies_ev, note)
end
add_ref_file("tupac_calibronium", "20140516_tupac_calibronium_4_row_4hz_chan1.ljh", "", [4952.22, 5898.801, 6490.59, 6404.01, 6930.378, 7649.45, 8047.823, 7058.18, 5426.962, 8905.41, 4090.74, 7478.252],
"calibronium in microxsi meausrement in tupac, roughly 5 eV FWHM")
add_ref_file("tupac_fe_emission", "20150730_50mM_irontris_100ps_xes_chan1.ljh", "20150730_50mM_irontris_100ps_xes_noise_chan1.ljh", [6404.1, 7058.18],
"fe emission measurement from iron-tris-bypyridine excited by water-jet-laser-source in tupac, roughly 5 eV FWHM, large background, few counts")
add_ref_file("good_mnka_mystery", "20150707_D_chan13.ljh", "20150707_C_chan13.noi", [5898.801, 6490.59],
"mn emission measurement with xsi in mystery(?), measured by Kelsey, 2.5 eV predicted, 2.44 eV achieved")

end # module


text:
new_pardiso_handle() = zeros(Int,64)

new_iparm() = zeros(Int32,64)

new_iparm_64() = zeros(Int64,64)


text:
include("iar.jl");
include("iar_lr.jl");
include("tiar.jl");

function infgmres(prob,m;method=:iar)
    n=size(prob.b,1);
    σ=0.0;
    start_time=time();
    println("Infgmres version: $method");
    println("Precomputing");
    Afact=factorize(compute_Mder(prob,σ));

    btilde=-(Afact\prob.b);
    precomp_time=time()-start_time;

    if (method == :iar)
        (Q,H,hist)=iar(prob,btilde,m,Afact);
    elseif (method == :iar_lr)
        (Q,H,hist)=iar_lr(prob,btilde,m,Afact);
    elseif (method == :tiar)
        (Z,A_tens,H,hist)=tiar(Float64,prob,btilde,m,Afact)
        Q = zeros(n,m);
        for j=1:m
            for l=1:m
                Q[1:n,j] += A_tens[1,j,l]*Z[:,l];
            end
        end
    else
        error("Unknown method $method");
    end

    hist[:precomp_time]=[precomp_time];

    e = zeros(m+1,1); e[1] = 1;

    solver=(mu,l) -> Q[1:n,1:l]*((mu*H[1:(l+1),1:l]-Matrix{Float64}(I,l+1,l))\e[1:(l+1)]*norm(btilde));
    return (solver,hist)

end

max_stars_repo_name: siobhanie/InfGMRES

max_issues_repo_name: siobhanie/InfGMRES

max_forks_repo_name: siobhanie/InfGMRES


text:
@testset "Pauli matrices" begin
    @test σˣ^2 == σ⁰
    @test σʸ^2 == σ⁰
    @test σᶻ^2 == σ⁰
    @test σˣ*σʸ - σʸ*σˣ == 2im*σᶻ
    @test σʸ*σᶻ - σᶻ*σʸ == 2im*σˣ
    @test σᶻ*σˣ - σˣ*σᶻ == 2im*σʸ
end


text:
import Dagger: debug_compute, get_logs!, LocalEventLog

export start_tracking_time, stop_tracking_time, tracktime, fetch_timings!

function time_table(log; profile=false)

    idx = Columns(proc=map(x->x.pid, getfields(log, :timeline)),
                  event_type=getfields(log, :category),
                  event_id=getfields(log, :id))

    if profile
        data = Columns(start=getfields(log, :start),   finish=getfields(log, :finish),
                       gc_diff=getfields(log, :gc_diff), profile=getfields(log, :profiler_samples))
    else
        data = Columns(start=getfields(log, :start),   finish=getfields(log, :finish),
                       gc_diff=getfields(log, :gc_diff))
    end

    NDSparse(idx, data)
end

function add_gc_diff(x,y)
    Base.GC_Diff(
        x.allocd     + y.allocd,
        x.malloc     + y.malloc,
        x.realloc    + y.realloc,
        x.poolalloc  + y.poolalloc,
        x.bigalloc   + y.bigalloc,
        x.freecall   + y.freecall,
        x.total_time + y.total_time,
        x.pause      + y.pause,
        x.full_sweep + y.full_sweep
    )
end

function aggregate_profile(xs)
    treereduce(Dagger.mix_samples, xs)
end

function aggregate_events(xs)
    sort!(xs, by=x->x.start)
    gc_diff = reduce(add_gc_diff, map(x -> x.gc_diff, xs))
    time_spent = sum(map(x -> x.finish - x.start, xs))
    if isdefined(xs[1], :profile)
        time_spent, gc_diff, aggregate_profile(map(x->x.profile, xs))
    end
    time_spent, gc_diff
end

function show_timings(t; maxdepth=5)
    # first aggregate by type of event
    t1 = reducedim_vec(aggregate_events, t, [:proc, :event_id])

    foreach(t1.index, t1.data) do i, x
        time_spent, gc_diff = x
        print(string(i[1]),": ")
        Base.time_print(time_spent, gc_diff.allocd, gc_diff.total_time, Base.gc_alloc_count(gc_diff))
    end

    t2 = reducedim_vec(aggregate_events, t, :event_id)
    println("Breakdown:")
    println(map(x->first(x)/1e9, t2))
    if isdefined(t.data.columns, :profile)
        p = aggregate_profile(t.data.columns.profile)
        if !isempty(p.samples)
            println("\nProfile output:")
            Profile.print(p.samples, p.lineinfo, maxdepth=maxdepth)
        end
    end
end

function getfields(log, fieldname)
    map(x->getfield(x, fieldname), log)
end

function start_tracking_time(;profile=false)
    ctx = get_context()
    dbgctx = Context(procs(ctx), LocalEventLog(), profile)
    compute_context[] = dbgctx
end

function stop_tracking_time()
    compute_context[] = nothing
end

"""
`tracktime(f)`

Track the time spent on different processes in different
categories in running `f`.
"""
function tracktime(f; profile=false, maxdepth=5)
    start_tracking_time(profile=profile)
    res = f()
    ctx = compute_context[]
    stop_tracking_time()
    t = fetch_timings!(ctx, profile=profile)
    show_timings(t, maxdepth=maxdepth)
    t, res
end

function fetch_timings!(ctx=get_context(); profile=true)
    time_table(Dagger.get_logs!(ctx.log_sink), profile=profile)
end


text:
"Module for tracing a ray."
module trace

export trace_pixel

using ..model
using ..constants
using ..geometry

using DifferentialEquations
using SpecialFunctions

# Some of setups for the integrator.
# const dt = 1e-3 * AU
alg = Tsit5()
# const rt=1e-6
# const at=1e-8

const rt=Float64[1e-8, 1e-12]
const at=Float64[1e-8, 1e-12]
# const dtmax = 5 * AU

const tau_thresh = 8.0

# Set up the terminator callback, stop ODE integration when this condition is satisfied.
function condition(u, t, integrator)
    u[1] > tau_thresh # when tau = 8
end
affect!(integrator) = terminate!(integrator)

cb = DiscreteCallback(condition, affect!)

# The integrand. p encapsulates the passed parameters.
function f(du, u, args, t)

    zstart, v0, pars, interp, gpre = args

    zprime = zstart - t

    # Write in coordinate conversions for querying alpha and sfunc
    rcyl, z, vlos = get_coords(zprime, gpre)

    # Calculate Delta v at this position
    Deltav = v0 * 1.e5 - vlos # get_vlos(xprime, yprime, zprime, pars)

    # Look up RT quantities from nearest-neighbor interp.
    DV2, sfunc, Upsilon = interp(rcyl, z)

    # Calculate alpha
    alpha = Upsilon * exp(-Deltav^2/DV2)

    # Update dtau/ds and dI/ds
    du[1] = alpha
    du[2] = exp(-u[1]) * sfunc * alpha
end

"v0 is the central velocity of the channel, relative to the disk systemic velocity."
function trace_pixel(xprime, yprime, v0, mol::Molecule, pars::AbstractParameters, interp; dt=0.5 * AU, dtmax=5 * AU)

    rmax = 700.0 * AU

    # Calculate DeltaVmax for this ray
    rcyl_min = abs(xprime)
    DeltaVmax = sqrt(2 * kB * temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2) * 1e-5 # km/s

    # Calculate the bounding positions for emission
    # If it's two, just trace it.
    # If it's four, break it up into two integrals.

    bounds = get_bounding_zps(xprime, yprime, pars, v0, DeltaVmax, rmax)
    if length(bounds) == 4
        z1start, z1end, z2start, z2end = bounds
    elseif length(bounds) == 2
        z1start, z1end = bounds
    end

    # Start the ode to integrate the first path.
    u0 = Float64[0.0; 0.0]

    gpre = precalc_geo(xprime, yprime, pars)

    args = (z1start, v0, pars, interp, gpre)

    tspan = (0.0, z1start - z1end)

    prob = ODEProblem(f, u0, tspan, args)

    sol = solve(prob, alg, callback=cb, reltol=rt, abstol=at, dense=true, save_everystep=true, dtmax=dtmax, dt=dt)
    # println(sol)

    tau_final, int_final = sol.u[end]dtmax

    zps = z1start - sol.t

    # println(sol.t ./ AU)
    # println("Number of evaluations ", length(sol.t))

    if (tau_final < tau_thresh) & (length(bounds) == 4)
        #We didn't reach the threshold in the first bounding region, so continue in the second.
        tspan = (0, z2start - z2end)

        # Start the ode to integrate the first path from the end of the second.
        u0 = Float64[tau_final; int_final]
        args = (z2start, v0, pars, interp, gpre)
        prob = ODEProblem(f, u0, tspan, args)
        sol = solve(prob, alg, callback=cb, reltol=rt, abstol=at, dense=true, save_everystep=true, dtmax=dtmax, dt=dt)
        tau_final, int_final = sol.u[end]
        zps = z2start - sol.t
        # return int_final
        return (zps, sol)

    else
        # return int_final
        return (zps, sol)
    end
end

end # module


text:
function lll{Td}(H::Array{Td,2},δ=3/4)
# (B,T,Q,R) = LLL(H,δ=3/4)
#  Do Lenstra–Lenstra–Lovász lattice reduction of matrix H using optional
#  parameter δ.  The output is B, an LLL-reduced basis; T, a unimodular
#  (det(T)=+/-1) transformation matrix such that B= H*T; Q and R such that
#  B=Q*R and Q is orthonormal, and R is upper triangular.  So H = Q*R*inv(T)
#
#  H can be of Integer, FloatingPoint, BigInt, or BigFloat types. The core
#  algorithm is designed for floating-point.
#
# Example with large real matrix:
#   N=500;H = randn(N,N); (B,T) = lll(H)
# Example with 2x2 complex matrix:
#   N=2;H = randn(N,N)+im*randn(N,N); (B,T) = lll(H)
# Example with 2x2 BigInt matrix:
#   N=2;H = ones(BigInt,N,N); rand!(H,-1e10:1e10); (B,T) = lll(H)

# Follows LLL in D. Wueben, et al, MMSE-Based Lattice-Reduction for Near-ML
# Detection of MIMO Systems International IEEE Workshop on Smart Antennas,
# Munich, March 2004.

# A few cycles can be saved by skipping updates of the Q matrix.

if !(0.25 < δ < 1.0)
    error("δ must be between 1/4 and 1.");
end

B = copy(H);
L = size(B,2);
(Q,R) = qr(B);
if Td<:Complex
    T = eye(Complex{Int},L)
    roundf(r) = round(real(r)) + im*round(imag(r));
else
    T = eye(Int,L);
    roundf(r) = round(r);
end

lx  = 2;
while lx <= L

    # reduce lx-th column of B
    for k=lx-1:-1:1
        rk = R[k,lx]/R[k,k]
        mu = roundf(rk)
        if abs(mu)>0
            #B[:,lx]   = B[:,lx]   - mu * B[:,k]
            axpy!(-mu, B[:,k], view(B,:,lx))
            #R[1:k,lx] = R[1:k,lx] - mu * R[1:k,k]
            axpy!(-mu, R[1:k,k], view(R,1:k,lx))
            #T[:,lx]   = T[:,lx]   - mu * T[:,k]
            axpy!(-mu, T[:,k],view(T,:,lx))
        end
    end

    nrm = norm(R[lx-1:lx,lx])
    if δ*abs(R[lx-1,lx-1])^2 > nrm^2

        # swap columns lx-1 and lx in B, T and R
        B[:,[lx-1,lx]]    = B[:,[lx,lx-1]];
        T[:,[lx-1,lx]]    = T[:,[lx,lx-1]];
        R[1:lx,[lx-1,lx]] = R[1:lx,[lx,lx-1]];

        # upper triangular by Givens rotation
        # mult with matrix Θ achieves R[lx,lx-1] = 0
        cc = R[lx-1,lx-1] / nrm # nrm = ||R[lx-1:lx,lx-1]|| after swapping
        ss = R[lx,lx-1]   / nrm
        Θ = [cc' ss; -ss cc]

        # Don't have BLAS working here yet
        R[lx-1:lx,lx-1:end] = Θ * R[lx-1:lx,lx-1:end]
        #gemm!('N','N',1.0,Θ,R[lx-1:lx,lx-1:end],0.0,R[lx-1:lx,lx-1:end])
        Q[:,lx-1:lx] = Q[:,lx-1:lx] * Θ'
        #gemm!('N','C', 1.0, Q[:,lx-1:lx], Θ, 0.0, Q[:,lx-1:lx])
        lx = max(lx-1,2)
    else
        lx = lx+1;
    end
end

return (B,T,Q,R)
end


text:
###############
# Remote Data #
###############

DATA_DIR = abspath(joinpath(abspath(string(@__DIR__)), "../data"))

function download_file(url, file)
    filepath = abspath(joinpath(DATA_DIR, file))
    @debug("Downloading datafile. URL: $url DESTINATION: $filepath")

    # Remove any temp files
    run(`rm -f $(filepath)_tmp`)

    # Attempt to download data
    run(`curl -S -L $url -o $(filepath)_tmp`)

    # Move temporary file into permanent location
    run(`rm -f $filepath`)
    run(`mv $(filepath)_tmp $filepath`)
end

export download_kp
"""
Download geomagnetic indices.

Arguments:
- `year_start::Int` First year to download data for. Default: 2000
- `year_end::Int` Last year to download data for. Default: 2019

Notes:
1. Data source is GFZ Potsdam Geomagnetic WDC tables: https://www.gfz-potsdam.de/en/kp-index/
"""
function download_kp(year_start::Int=2000, year_end::Int=Dates.year(Dates.today()))
    for year in year_start:year_end
        download_file("ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/wdc/yearly/kp$year.wdc", "kp$year.wdc")
    end

    # Merge Data and clean up
    open(abspath(joinpath(DATA_DIR, "kpall.wdc")), "w") do fp_kpall
        # For each file on disk: 1. Open and read it, 2. Write it to all 3. Delete it
        for year in year_start:year_end
            kp_filepath = abspath(joinpath(DATA_DIR, "kp$year.wdc"))
            open(kp_filepath) do fp_kpfile
                # Read and write all in one operation
                write(fp_kpall, read(fp_kpfile, String))

                # Delete File from disk
                rm(kp_filepath)
            end
        end
    end
end

export download_solar_flux
"""
Download F10.7cm Solar Flux data.

10.7cm solar flux is the standard measure of solar activity in space weather models.

Notes:
1. Data source is NRC Canada solar flux tables: ftp://ftp.seismo.nrcan.gc.ca/spaceweather/solar_flux/daily_flux_values/fluxtable.txt
"""
function download_solar_flux()
    download_file("ftp://ftp.seismo.nrcan.gc.ca/spaceweather/solar_flux/daily_flux_values/fluxtable.txt", "fluxtable.txt")
end

export download_all_data
"""
Downloads package datafiles into folders `\$PACKAGE_ROOT/DIR`

Downloads the following files:
- IERS C04 IAU2000A Earth Orientation Data
- IERS C04 IAU1980 Earth Orientation Data
- IERS Bulletin A/B IAU2000 Earth Orientation Data
"""
function download_all_data()
    # Earth Orientation Data
    download_file("https://datacenter.iers.org/data/latestVersion/224_EOP_C04_14.62-NOW.IAU2000A224.txt", "EOP_C04_14.62-NOW.IAU2000A.txt")
    download_file("https://datacenter.iers.org/data/latestVersion/223_EOP_C04_14.62-NOW.IAU1980223.txt", "EOP_C04_80.62-NOW.IAU2000A.txt")
    download_file("https://datacenter.iers.org/data/latestVersion/9_FINALS.ALL_IAU2000_V2013_019.txt", "FINALS.ALL_IAU2000.txt")

    # Geomagnetic inidies
    download_kp()

    # Solar Flux
    download_solar_flux()
end

##########################
# Earth Orientation Data #
##########################

# Define product dictionary
eop_products = Dict(
    "C04_14" => ("https://datacenter.iers.org/data/latestVersion/224_EOP_C04_14.62-NOW.IAU2000A224.txt", abspath(string(@__DIR__), "../data/EOP_C04_14.62-NOW.IAU2000A.txt")),
    "C04_80" => ("https://datacenter.iers.org/data/latestVersion/223_EOP_C04_14.62-NOW.IAU1980223.txt", abspath(string(@__DIR__), "../data/EOP_C04_80.62-NOW.IAU2000A.txt")),
    "FINALS_2000" => ("https://datacenter.iers.org/data/latestVersion/9_FINALS.ALL_IAU2000_V2013_019.txt", abspath(string(@__DIR__), "../data/FINALS.ALL_IAU2000.txt"))
)

export EarthOrientationData
"""
The EarthOrientationData constains a single data member of type
`Dict{Int, Tuple{Float64, Float64, Float64}}` that stores the Earth
Orientation parameters `UT1-UTC`, `xp`, and `yp` whose units are _meters_,
_radians_, and _radians_, respectively. `xp` and `yp` are the x- and
y-components of Earth's polar motion. The dictionary key is the Epoch the
parameters are for as a Modified Julian Day at 0h UTC.

Arguments:
- `product::Symbol` The IERS product type can be `"C04_14"`, `"C04_80"`, or `"FINALS_2000"`
"""
struct EarthOrientationData
    data::Dict{Int, Tuple{Float64, Float64, Float64}}
end

function EarthOrientationData(product::String)
    # Initialize Data Array
    eop_data = Dict{Int, Tuple{Float64, Float64, Float64}}()

    # Load in Data from filepath
    if product == "FINALS_2000"
        for line in readlines(eop_products[product][2])
            if line[17] == 'P' || line[17] == 'I'
                mjd_utc = parse(Int, line[8:12])            # MJD (UTC)
                ut1_utc = parse(Float64, line[59:68])         # UT1-UTC [s]
                xp      = parse(Float64, line[19:27])*AS2RAD  # xp [rad]
                yp      = parse(Float64, line[38:46])*AS2RAD  # yp [rad]

                eop_data[mjd_utc] = (ut1_utc, xp, yp)
            end
        end
    elseif product == "C04_14" || product == "C04_80"
        open(eop_products[product][2], "r") do product_file
            for i in 1:14
                # Read first 14 lines to skip to data
                readline(product_file)
            end

            for line in readlines(product_file)
                split_line = split(line)
                mjd_utc = parse(Int, split_line[4])
                ut1_utc = parse(Float64, split_line[7])
                xp      = parse(Float64, split_line[5])*AS2RAD
                yp      = parse(Float64, split_line[6])*AS2RAD

                eop_data[mjd_utc] = (ut1_utc, xp, yp)
            end
        end
    else
        error("Unknown symbol $(String(:product))")
    end

    return EarthOrientationData(eop_data)
end

# Declare global Earth Orientation Data Object used by Reference System Calls
export EOP
"""
Module-wide global EarthOrientationData object. This data object is used as the
default source of Earth Orientation Data by reference system transformations if
no explicit EarthOrientationData file is provided to those transformations.

This value can be overridden in your own code as follows:

```julia
SatelliteDynamics.EOP = EarthOrientationData(:EOP_PRODUCT_CHOICE)
```

This global variable defaults to use the module's internal version of `"FINALS_2000"`
if it is not otherwise set/provided.
"""
global EOP = EarthOrientationData("FINALS_2000")

# Access Methods
export UT1_UTC
"""
Compute the offset between the UT1 and UTC time systems in seconds. If the EarthOrientationData argument is ommitted the function will use the default module-global value.

Arguments:
- `eop::EarthOrientationData` EarthOrientationData object to use to compute the offset
- `mjd::Real` Modified Julian Date in UTC of the Epoch for which the UT1-UTC offset is desired.
- `interp::Bool` Whether to linearly interpolate the parameter data to the input MJD.

Returns:
- `ut1_utc::Float` UT1 - UTC offset. [s]
"""
function UT1_UTC(eop::EarthOrientationData, mjd::Real; interp::Bool=false)
    if interp
        x1 = floor(mjd)
        x2 = floor(mjd) + 1
        y1 = eop.data[convert(Int, floor(mjd))][1]
        y2 = eop.data[convert(Int, floor(mjd)+1)][1]
        x  = (y2 - y1)/(x2 - x1) * (mjd - x1) + y1
        return x
    else
        return eop.data[convert(Int, floor(mjd))][1]
    end
end

UT1_UTC(mjd::Real; interp::Bool=false) = UT1_UTC(EOP, mjd, interp=interp)

export POLE_LOCATOR
"""
Compute the location of the pole. Returns x- and y- components as a tuple with the units of [radians].  If the EarthOrientationData argument is ommitted the function will use the default module-global value.

Arguments:
- `eop::EarthOrientationData` EarthOrientationData object to use to compute the offset
- `mjd::Real` Modified Julian Date in UTC of the Epoch for which the pole locator is desired.
- `interp::Bool` Whether to linearly interpolate the parameter data to the input MJD.

Returns:
- `pole_locator::Tuple{ -Float, Float}` (x, y) pole location in radians.
"""
function POLE_LOCATOR(eop::EarthOrientationData, mjd::Real; interp::Bool=false)
    if interp
        x1 = floor(mjd)
        x2 = floor(mjd) + 1

        # Get values converted to array for interpolation
        y1 = [v for v in eop.data[convert(Int, floor(mjd))][2:3]]
        y2 = [v for v in eop.data[convert(Int, floor(mjd)+1)][2:3]]
        x  = (y2 - y1)/(x2 - x1) * (mjd - x1) + y1
        return x
    else
        return eop.data[convert(Int, floor(mjd))][2:3]
    end
end

POLE_LOCATOR(mjd::Real; interp::Bool=false) = POLE_LOCATOR(EOP, mjd, interp=interp)

export XP
"""
Compute the x-component of the pole locator in [radians]. If the first EarthOrientationData argument is ommitted the function will use the default module-global value.

Arguments:
- `eop::EarthOrientationData` EarthOrientationData object to use to compute the offset
- `mjd::Real` Modified Julian Date in UTC of the Epoch for which the xp value is desired.
- `interp::Bool` Whether to linearly interpolate the parameter data to the input MJD.

Returns:
- `xp::Float` x-component of pole locator in radians.
"""
function XP(eop::EarthOrientationData, mjd::Real; interp=false)
    if interp
        x1 = floor(mjd)
        x2 = floor(mjd) + 1
        y1 = eop.data[convert(Int, floor(mjd))][2]
        y2 = eop.data[convert(Int, floor(mjd)+1)][2]
        x  = (y2 - y1)/(x2 - x1) * (mjd - x1) + y1
        return x
    else
        return eop.data[convert(Int, floor(mjd))][2]
    end
end

XP(mjd::Real; interp::Bool=false) = XP(EOP, mjd, interp=interp)

export YP
"""
Compute the y-component of the pole locator in [radians]. If the first EarthOrientationData argument is ommitted the function will use the default module-global value.

Arguments:
- `eop::EarthOrientationData` EarthOrientationData object to use to compute the offset
- `mjd::Real` Modified Julian Date in UTC of the Epoch for which the yp value is desired.
- `interp::Bool` Whether to linearly interpolate the parameter data to the input MJD.

Returns:
- `yp::Float` y-component of pole locator in radians.
"""
function YP(eop::EarthOrientationData, mjd::Real; interp::Bool=false)
    if interp
        x1 = floor(mjd)
        x2 = floor(mjd) + 1
        y1 = eop.data[convert(Int, floor(mjd))][3]
        y2 = eop.data[convert(Int, floor(mjd)+1)][3]
        x  = (y2 - y1)/(x2 - x1) * (mjd - x1) + y1
        return x
    else
        return eop.data[convert(Int, floor(mjd))][3]
    end
end

YP(mjd::Real; interp::Bool=false) = YP(EOP, mjd, interp=interp)

export set_eop
"""
Set Earth orientation data values for a specific date in the module global EarthOrientationData object.

Arguments:
- `mjd::Real` Modified Julian Date in UTC of the Epoch for which the Earth orientation data is aligned to.
- `ut1_utc::Real` Offset between UT1 and UTC in seconds.
- `xp::Real` x-component of the pole locator in radians.
- `yp::Real` y-component of the pole locator in radians.
"""
function set_eop(mjd::Real, ut1_utc::Real, xp::Real, yp::Real)
    EOP.data[convert(Int, floor(mjd))] = (ut1_utc, xp*AS2RAD, yp*AS2RAD)
end

export load_eop
"""
Load new Earth orientation data into the module global EarthOrientationData object. The product can be one of the symbols: `"C04_14"`, `"C04_80"`, or `"FINALS_2000"`.

Arguments:
- `product::String` Loads a different set of EarthOrientationData values into the module-wide global EarthOrientationData parameters.
"""
function load_eop(product::String)
    global EOP = EarthOrientationData(product::String)
end

#################
# Gravity Model #
#################

grav_products = Dict(
    :EGM2008_20 => abspath(string(@__DIR__), "../data/EGM2008_20.gfc"),
    :EGM2008_90 => abspath(string(@__DIR__), "../data/EGM2008_90.gfc"),
    :GGM01S => abspath(string(@__DIR__), "../data/GGM01S.gfc"),
    :GGM01S => abspath(string(@__DIR__), "../data/GGM05S.gfc"),
)

function line_starts_with(line::String, str::String)
    if length(line) > length(str) && line[1:length(str)] == str
        return true
    else
        return false
    end
end

export GravModel
"""
GravModel stores a spherical harmonic gravity field in memory. Can store normalized or denomalized coefficients. Package contains EGM2008, GGM01S, and GGM0S gravity models, as well as the default gravity model of EGM2008 truncated to degree and order 90.

Additional gravity field models can be downloaded from: <http://icgem.gfz-potsdam.de/home>

Arguments:
- `filepath::string` Path to spherical harmonic gravity model file.
"""
struct GravModel
    name::String
    normalized::Bool
    R::Float64
    GM::Float64
    n_max::Int64
    m_max::Int64
    data::Array{Float64, 2}
end

function GravModel(filepath::String)
    model_name = ""
    normalized = false
    R          = 0.0
    GM         = 0.0
    n_max      = 0.0
    m_max      = 0.0
    data       = Array{Float64, 2}(undef, 1, 1)

    # Parse File
    for line in readlines(filepath)
        # Replace non-standard float formatting in GFC files
        line = replace(line, "D+" => "e+")
        line = replace(line, "D-" => "e-")

        if line_starts_with(line, "modelname")
            model_name = split(line)[2]
        elseif line_starts_with(line, "max_degree")
            n_max = parse(Int16, split(line)[2])
            m_max = n_max
            data = zeros(Float64, n_max+1, m_max+1)
        elseif line_starts_with(line, "earth_gravity_constant")
            GM = parse(Float64, split(line)[2])
        elseif line_starts_with(line, "radius")
            R = parse(Float64, split(line)[2])
        elseif line_starts_with(line, "norm")
            if split(line)[2] == "fully_normalized"
                normalized = true
            else
                normalized = false
            end
        elseif line_starts_with(line, "gfc")
            line_split = split(line)
            n = parse(Int16, line_split[2])
            m = parse(Int16, line_split[3])
            C = parse(Float64, line_split[4])
            S = parse(Float64, line_split[5])

            data[n+1, m+1] = C
            if m != 0
                data[m+1-1, n+1] = S
            end
        end
    end

    # Pre-Allocate
    return GravModel(model_name, normalized, R, GM, n_max, m_max, data)
end

function GravModel(product_name::Symbol)
    return GravModel(grav_products[product_name])
end

# Declare glrobal Gravity Model used by dynamics model calls
export GRAVITY_MODEL
"""
Module-wide global GravityModel object. This data object is used as the
default spherical harmonic gravity field unless one is otherwise provided.

This value can be overridden in your own code as follows:

```julia
SatelliteDynamics.GravityModel = GravityModel(PATH_TO_YOUR_GRAVITY_MODEL)
```

This global variable defaults to use the module's internal version of the EGM2008 model truncated to order and degree 90, if it is not otherwise set.
"""
global GRAVITY_MODEL = GravModel(abspath(@__DIR__, "../data/EGM2008_90.gfc"))

export load_gravity_model
"""
Load new gravity model into module global EarthOrientationData object. The product can be one of the symbols: `:EGM2008_20`, `:EGM2008_90`, `:GGM01S`, `:GGM05S`, or the filepath to a text-encoded gravity model file.

Arguments:
- `gfc_file::String` File path of gravity field model
- `product_name::Symbol` _OR_ a symbol of a known gravity field product. Valid ones are: `:EGM2008_20`, `:EGM2008_90`, `:GGM01S`, `:GGM05S`
"""
function load_gravity_model_file(gfc_file::String)
    global GRAVITY_MODEL = GravModel(gfc_file::String)
end

function load_gravity_model(product_name::Symbol)
    global GRAVITY_MODEL = GravModel(product_name::Symbol)
end

export GRAV_COEF
function GRAV_COEF(i::Int, j::Int)
    # Offset into matrix to deal with julia indexing
    return GRAVITY_MODEL.data[i+1, j+1]
end

max_stars_repo_name: sisl/SatelliteDynamics

max_issues_repo_name: sisl/SatelliteDynamics

max_forks_repo_name: sisl/SatelliteDynamics


text:
N = 10
NG = N#5

fooErgmDirBin1
deltaN = 3

tmpdegsIO = round.([ Vector(linspace(deltaN,N-1-deltaN,N));Vector(linspace(deltaN,N-1-deltaN,N))])
tmpdegsIO  =   DegSeq2graphDegSeq(fooErgmDirBin1,tmpdegsIO)
linSpacedPar(fooErgmDirBin1,N; Ngroups = NG, deltaN=deltaN)
## test graphicability functions

##
groupsInds = distributeAinVecN(Vector(1:NG),N)
Mod = ErgmDirBin1(deg,zeros(N),groupsInds)
groupsNames = unique(groupsInds)
NG = length(groupsNames)
tmpDic = countmap(groupsInds)
#number of nodes per group
NGMemb = sortrows([[i for i in keys(tmpDic)] [  tmpDic[i] for i in keys(tmpDic) ] ])
# agrregate degree of each group
degGroups = zeros(2NG) ; for i =1:N degGroups[groupsInds[i]] += deg[i]; degGroups[NG + groupsInds[i]] += deg[N + i];  end
L = sum(degGroups,1)[1]
L/(N^2-N)
groupsPar = ones(2NG)
firstOrderCond(Mod;degIO = deg,parGroupsIO = groupsPar, groupsInds = groupsInds )
expMatrix(Mod,groupsPar)

targetErrValStaticNets = 0.005
estPar,estIt,estMod = estimate(Mod)

##
sam = sampl(estMod,1000)

meanMatSam = squeeze(mean(sam,3),3)
degsSam = expValStatsFromMat(estMod,meanMatSam)
degGroupsSam = zeros(2NG) ; for i =1:N degGroupsSam[groupsInds[i]] += degsSam[i]; degGroupsSam[NG + groupsInds[i]] += degsSam[ N + i] end
relErrSam = maximum(abs.(( degGroupsSam - degGroups )./degGroups )[degGroups.!=0])

##

#


text:
export dimension, transform, transform_and_logjac, transform_logdensity, inverse, inverse!,
    inverse_eltype, as, random_arg, random_value

###
### log absolute Jacobian determinant
###

"""
$(TYPEDEF)

Flag used internally by the implementation of transformations, as explained below.

When calculating the log jacobian determinant for a matrix, initialize with
```julia
logjac_zero(flag, x)
```
and then accumulate with log jacobians as needed with `+`.

When `flag` is `LogJac`, methods should return the log Jacobian as the second
argument, otherwise `NoLogJac`, which simply combines to itself with `+`,
serving as an empty placeholder. This allows methods to share code of the two
implementations.
"""
abstract type LogJacFlag end

"""
Calculate log Jacobian as the second value.
"""
struct LogJac <: LogJacFlag end

const LOGJAC = LogJac()

"""
Don't calculate log Jacobian, return `NOLOGJAC` as the second value.
"""
struct NoLogJac <: LogJacFlag end

const NOLOGJAC = NoLogJac()

Base.:+(::NoLogJac, ::NoLogJac) = NOLOGJAC

"""
$(SIGNATURES)

Initial value for log Jacobian calculations.
"""
logjac_zero(::LogJac, T::Type{<:Real}) = log(one(T))

logjac_zero(::NoLogJac, _) = NOLOGJAC

###
### internal methods that implement transformations
###

"""
    transform_with(flag::LogJacFlag, transformation, x::AbstractVector, index)

Transform elements of `x` from `index`, using `transformation`.

Return `(y, logjac), index′`, where

- `y` is the result of the transformation,

- `logjac` is the the log Jacobian determinant or a placeholder, depending on `flag`,

- `index′` is the next index in `x` after the elements used for the transformation

**Internal function**. Implementations

1. can assume that `x` has enough elements for `transformation` (ie `@inbounds` can be
used),

2. should work with generalized indexing on `x`.
"""
function transform_with end

"""
    inverse_at!(x, index, transformation, y)

Invert transformation at `y` and put the result in `x` starting at `index`.

**Internal function**. Implementations

1. can assume that `x` has enough elements for the result (ie `@inbounds` can be used),

2. should work with generalized indexing on `x`.
"""
function inverse_at! end

####
#### API
####

"""
$(TYPEDEF)

Supertype for all transformations in this package.

# Interface

The user interface consists of

- [`dimension`](@ref)
- [`transform`](@ref)
- [`transform_and_logjac`](@ref)
- [`inverse`]@(ref), [`inverse!`](@ref)
- [`inverse_eltype`](@ref).
"""
abstract type AbstractTransform end

Base.broadcastable(t::AbstractTransform) = Ref(t)

"""
$(TYPEDEF)

Inverse of the wrapped transform. Use the 1-argument version of
[`inverse`](@ref) to construct.
"""
struct InverseTransform{T}
    transform::T
end

"""
    inverse(t::AbstractTransform, y)

Return `x` so that `transform(t, x) ≈ y`.

    $SIGNATURES

Return a callable equivalen to `y -> inverse(t, y)`.
"""
inverse(t::AbstractTransform) = InverseTransform(t)

(ι::InverseTransform)(y) = inverse(ι.transform, y)

"""
    inverse_eltype(t::AbstractTransform, y)

The element type for vector `x` so that `inverse!(x, t, y)` works.
"""
function inverse_eltype end

"""
$(SIGNATURES)

Put `inverse(t, y)` into a preallocated vector `x`, returning `x`.

Generalized indexing should be assumed on `x`.

See [`inverse_eltype`](@ref) for determining the type of `x`.
"""
function inverse!(x::AbstractVector, transformation::AbstractTransform, y)
    @argcheck dimension(transformation) == length(x)
    inverse_at!(x, firstindex(x), transformation, y)
    x
end

"""
$(SIGNATURES)

Let ``y = t(x)``, and ``f(y)`` a log density at `y`. This function evaluates `f
∘ t` as a log density, taking care of the log Jacobian correction.
"""
function transform_logdensity(t::AbstractTransform, f, x)
    y, ℓ = transform_and_logjac(t, x)
    ℓ + f(y)
end

"""
    dimension(t::AbstractTransform)

The dimension (number of elements) that `t` transforms.

Types should implement this method.
"""
function dimension end

"""
    as(T, args...)

Shorthand for constructing transformations with image in `T`. `args` determines
or modifies behavior, details depend on `T`.

Not all transformations have an `as` method, some just have direct constructors.
See `methods(as)` for a list.

# Examples

```julia
as(Real, -∞, 1)          # transform a real number to (-∞, 1)
as(Array, 10, 2)         # reshape 20 real numbers to a 10x2 matrix
as(Array, as𝕀, 10)       # transform 10 real numbers to (0, 1)
as((a = asℝ₊, b = as𝕀)) # transform 2 real numbers a NamedTuple, with a > 0, 0 < b < 1
```
"""
function as end

####
#### vector transformations
####

"""
$(TYPEDEF)

Transformation that transforms `<: AbstractVector`s to other values.

# Implementation

Implements [`transform`](@ref) and [`transform_and_logjac`](@ref) via
[`transform_with`](@ref), and [`inverse`](@ref) via [`inverse!`](@ref).
"""
abstract type VectorTransform <: AbstractTransform end

"""
$(SIGNATURES)

Transform `x` using `t`.
"""
function transform(t::VectorTransform, x::AbstractVector)
    @argcheck dimension(t) == length(x)
    first(transform_with(NOLOGJAC, t, x, firstindex(x)))
end

"""
$(SIGNATURES)

Transform `x` using `t`; calculating the log Jacobian determinant, returned as
the second value.
"""
function transform_and_logjac(t::VectorTransform, x::AbstractVector)
    @argcheck dimension(t) == length(x)
    y, ℓ, _ = transform_with(LOGJAC, t, x, firstindex(x))
    y, ℓ
end

function inverse(t::VectorTransform, y)
    inverse!(Vector{inverse_eltype(t, y)}(undef, dimension(t)), t, y)
end

"""
$(SIGNATURES)

A random argument for a transformation.

# Keyword arguments

$(_RANDOM_REALS_KWARGS_DOC)
"""
random_arg(x::VectorTransform; kwargs...) = random_reals(dimension(x); kwargs...)

"""
$(SIGNATURES)

Random value from a transformation.

# Keyword arguments

$(_RANDOM_REALS_KWARGS_DOC)
"""
random_value(t::AbstractTransform; kwargs...) = transform(t, random_arg(t; kwargs...))


text:
 # dropout layer
function dropout(x,d)
    if d > 0
        return x .* (rand!(similar(AutoGrad.getval(x))) .> d) * (1/(1-d))
    else
        return x
    end
end

# LSTM model - input * weight, concatenated weights
function lstm(weight, bias, hidden, cell, input)
    gates   = hcat(input,hidden) * weight .+ bias
    hsize   = size(hidden,2)
    forget  = sigm(gates[:,1:hsize])
    ingate  = sigm(gates[:,1+hsize:2hsize])
    outgate = sigm(gates[:,1+2hsize:3hsize])
    change  = tanh(gates[:,1+3hsize:end])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh(cell)
    return (hidden,cell)
end

function logprob(output, ypred, mask=nothing)
    nrows,ncols = size(ypred)
    index = similar(output)
    @inbounds for i=1:length(output)
        index[i] = i + (output[i]-1)*nrows
    end
    o1 = logp(ypred,2)
    index = mask == nothing ? index : index[mask]
    o2 = o1[index]
    o3 = sum(o2)
    return o3
end

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
# Testing of simple and reference expressions on scalar context

import Devectorize
import Devectorize.@devec
import Devectorize.@inspect_devec
import Devectorize.dump_devec

using Base.Test

#data

a = [1., 2., 3., 4., 5., 6., 7., 8.]
b = [3., 4., 5., 6., 8., 7., 6., 5.]
c = [9., 8., 7., 6., 4., 2., 3., 1.]
cv = 12.0
cv_a = [12., 12., 12., 12., 12., 12., 12., 12.]

abc = [a b c]

at = a'
bt = b'
ct = c'
abct = [at; bt; ct]

type X
    x
end

ax = X(a)

###########################################################
#
#   trivial/simple assignments
#
###########################################################

@devec r = a
@test r === a

@devec r = ax.x
@test r === ax.x

r = X(0)
@devec r.x = a
@test r.x === a

r = X(0)
@devec r.x = ax.x
@test r.x === ax.x

r = similar(a)
r0 = r
#dump_devec(:(r[:] = b))
@devec r[:] = b
@test r === r0
@test isequal(r, b)

###########################################################
#
#   scalar or scalar references
#
###########################################################

@devec r[:] = 12.
@test isequal(r, cv_a)

@devec r[:] = cv
@test isequal(r, cv_a)

r = similar(a)
@devec r[:] = a[3]
@test isequal(r, fill(a[3], size(a)))

i = 2
@devec r[:] = a[i]
@test isequal(r, fill(a[i], size(a)))

@devec r[:] = abc[1,3]
@test isequal(r, fill(c[1], size(a)))

i = 2
@devec r[:] = abc[i,3]
@test isequal(r, fill(c[i], size(a)))

@devec r[:] = abc[2,i]
@test isequal(r, fill(b[2], size(a)))

j = 1
@devec r[:] = abc[i,j]
@test isequal(r, fill(a[i], size(a)))

###########################################################
#
#   1D references
#
###########################################################

# 1D ref on RHS

@devec r = a[:]
@test isequal(r, a)

@devec r = a[1:end]
@test isequal(r, a)

@devec r = a[1:5]
@test isequal(r, a[1:5])

v = 4
@devec r = a[1:v]
@test isequal(r, a[1:v])

@devec r = a[3:end]
@test isequal(r, a[3:end])

@devec r = a[3:6]
@test isequal(r, a[3:6])

v = 7
@devec r = a[3:v]
@test isequal(r, a[3:v])

u = 2
@devec r = a[u:end]
@test isequal(r, a[u:end])

@devec r = a[u:6]
@test isequal(r, a[u:6])

v = 7
@devec r = a[u:v]
@test isequal(r, a[u:v])

# 1D ref on LHS

r = zeros(size(a))
@devec r[:] = b[:]
rr = zeros(size(a))
rr[:] = b[:]
@test isequal(r, rr)

r = zeros(size(a))
@devec r[3:6] = cv
rr = zeros(size(a))
rr[3:6] = cv
@test isequal(r, rr)

r = zeros(size(a))
@devec r[3:6] = b[2:5]
rr = zeros(size(a))
rr[3:6] = b[2:5]
@test isequal(r, rr)

u = 1

r = zeros(size(a))
@devec r[1:end] = b[u:end]
rr = zeros(size(a))
rr[1:end] = b[u:end]
@test isequal(r, rr)

r = zeros(size(a))
@devec r[u:end] = b[:]
rr = zeros(size(a))
rr[u:end] = b[:]
@test isequal(r, rr)

u = 2

r = zeros(size(a))
@devec r[2:end] = b[u:end]
rr = zeros(size(a))
rr[2:end] = b[u:end]
@test isequal(r, rr)

r = zeros(size(a))
@devec r[u:end] = b[2:end]
rr = zeros(size(a))
rr[u:end] = b[2:end]
@test isequal(r, rr)

u, v = 4, 7

r = zeros(size(a))
@devec r[3:6] = b[u:v]
rr = zeros(size(a))
rr[3:6] = b[u:v]
@test isequal(r, rr)

r = zeros(size(a))
@devec r[u:v] = b[2:5]
rr = zeros(size(a))
rr[u:v] = b[2:5]
@test isequal(r, rr)

###########################################################
#
#   vector references
#
###########################################################

# integer vector
r[:] = a
rr[:] = a
u = 1:2:length(r)
@devec r[u] = 1
rr[u] = 1
@test isequal(r, rr)

@devec r[u] = b[u]
rr[u] = b[u]
@test isequal(r, rr)

# logical vector
r[:] = a
rr[:] = a
u = Bool[isodd(x) for x=1:length(r)]
@devec r[u] = 1
rr[u] = 1
@test isequal(r, rr)

@devec r[u] = b[u]
rr[u] = b[u]
@test isequal(r, rr)

###########################################################
#
#   column references
#
###########################################################

# col-ref on RHS

j = 2

@devec r = abc[:,1]
@test isequal(r, a)

@devec r = abc[1:end,1]
@test isequal(r, a)

@devec r = abc[1:5,1]
@test isequal(r, a[1:5])

@devec r = abc[:,j]
@test isequal(r, b)

@devec r = abc[1:end,j]
@test isequal(r, b)

@devec r = abc[1:5,j]
@test isequal(r, b[1:5])

v = 4
@devec r = abc[1:v,1]
@test isequal(r, a[1:v])

@devec r = abc[3:end,1]
@test isequal(r, a[3:end])

@devec r = abc[3:6,1]
@test isequal(r, a[3:6])

@devec r = abc[1:v,j]
@test isequal(r, b[1:v])

@devec r = abc[3:end,j]
@test isequal(r, b[3:end])

@devec r = abc[3:6,j]
@test isequal(r, b[3:6])

v = 7
@devec r = abc[3:v,1]
@test isequal(r, a[3:v])

@devec r = abc[3:v,j]
@test isequal(r, b[3:v])

u = 2
@devec r = abc[u:end,1]
@test isequal(r, a[u:end])

@devec r = abc[u:6,1]
@test isequal(r, a[u:6])

@devec r = abc[u:end,j]
@test isequal(r, b[u:end])

@devec r = abc[u:6,j]
@test isequal(r, b[u:6])

v = 7
@devec r = abc[u:v,1]
@test isequal(r, a[u:v])

@devec r = abc[u:v,j]
@test isequal(r, b[u:v])

# col-ref on LHS

r = zeros(size(abc))
@devec r[:,1] = b[:]
rr = zeros(size(abc))
rr[:,1] = b[:]
@test isequal(r, rr)

r = zeros(size(abc))
@devec r[3:6,1] = cv
rr = zeros(size(abc))
rr[3:6,1] = cv
@test isequal(r, rr)

r = zeros(size(abc))
@devec r[3:6,1] = b[2:5]
rr = zeros(size(abc))
rr[3:6,1] = b[2:5]
@test isequal(r, rr)

u = 1

r = zeros(size(abc))
@devec r[1:end,1] = b[u:end]
rr = zeros(size(abc))
rr[1:end,1] = b[u:end]
@test isequal(r, rr)

r = zeros(size(abc))
@devec r[u:end,1] = b[:]
rr = zeros(size(abc))
rr[u:end,1] = b[:]
@test isequal(r, rr)

j = 3
u = 2

r = zeros(size(abc))
@devec r[2:end,j] = b[u:end]
rr = zeros(size(abc))
rr[2:end,j] = b[u:end]
@test isequal(r, rr)

r = zeros(size(abc))
@devec r[u:end,j] = b[2:end]
rr = zeros(size(abc))
rr[u:end,j] = b[2:end]
@test isequal(r, rr)

u, v = 4, 7

r = zeros(size(abc))
@devec r[3:6,j] = b[u:v]
rr = zeros(size(abc))
rr[3:6,j] = b[u:v]
@test isequal(r, rr)

r = zeros(size(abc))
@devec r[u:v,j] = b[2:5]
rr = zeros(size(abc))
rr[u:v,j] = b[2:5]
@test isequal(r, rr)

###########################################################
#
#   row references
#
###########################################################

# row-ref on RHS

i = 2

@devec r = abct[1,:]
@test isequal(r, at)

@devec r = abct[1,1:end]
@test isequal(r, abct[1,1:end])

@devec r = abct[1,1:5]
@test isequal(r, abct[1,1:5])

@devec r = abct[i,:]
@test isequal(r, abct[i,:])

@devec r = abct[i,1:end]
@test isequal(r, abct[i,1:end])

@devec r = abct[i,1:5]
@test isequal(r, abct[i,1:5])

v = 4
@devec r = abct[1,1:v]
@test isequal(r, abct[1,1:v])

@devec r = abct[1,3:end]
@test isequal(r, abct[1,3:end])

@devec r = abct[1,3:6]
@test isequal(r, abct[1,3:6])

@devec r = abct[i,1:v]
@test isequal(r, abct[i,1:v])

@devec r = abct[i,3:end]
@test isequal(r, abct[i,3:end])

@devec r = abct[i,3:6]
@test isequal(r, abct[i,3:6])

v = 7
@devec r = abct[1,3:v]
@test isequal(r, abct[1,3:v])

@devec r = abct[i,3:v]
@test isequal(r, abct[i,3:v])

u = 2
@devec r = abct[1,u:end]
@test isequal(r, abct[1,u:end])

@devec r = abct[1,u:6]
@test isequal(r, abct[1,u:6])

@devec r = abct[i,u:end]
@test isequal(r, abct[i,u:end])

@devec r = abct[i,u:6]
@test isequal(r, abct[i,u:6])

v = 7
@devec r = abct[1,u:v]
@test isequal(r, abct[1,u:v])

@devec r = abct[i,u:v]
@test isequal(r, abct[i,u:v])

# row-ref on LHS

r = zeros(size(abct))
@devec r[1,:] = bt[:]
rr = zeros(size(abct))
rr[1,:] = bt[:]
@test isequal(r, rr)

r = zeros(size(abct))
@devec r[1,3:6] = cv
rr = zeros(size(abct))
rr[1,3:6] = cv
@test isequal(r, rr)

r = zeros(size(abct))
@devec r[1,3:6] = bt[2:5]
rr = zeros(size(abct))
rr[1,3:6] = bt[2:5]
@test isequal(r, rr)

u = 1

r = zeros(size(abct))
@devec r[1,1:end] = bt[u:end]
rr = zeros(size(abct))
rr[1,1:end] = bt[u:end]
@test isequal(r, rr)

r = zeros(size(abct))
@devec r[1,u:end] = bt[:]
rr = zeros(size(abct))
rr[1,u:end] = bt[:]
@test isequal(r, rr)

i = 2
u = 2

r = zeros(size(abct))
@devec r[i,2:end] = bt[u:end]
rr = zeros(size(abct))
rr[i,2:end] = bt[u:end]
@test isequal(r, rr)

r = zeros(size(abct))
@devec r[i,u:end] = bt[2:end]
rr = zeros(size(abct))
rr[i,u:end] = bt[2:end]
@test isequal(r, rr)

u, v = 4, 7

r = zeros(size(abct))
@devec r[i,3:6] = bt[u:v]
rr = zeros(size(abct))
rr[i,3:6] = bt[u:v]
@test isequal(r, rr)

r = zeros(size(abct))
@devec r[i,u:v] = bt[2:5]
rr = zeros(size(abct))
rr[i,u:v] = bt[2:5]
@test isequal(r, rr)

###########################################################
#
#   2D references
#
###########################################################

B = reshape(linspace(1., 64., 64), 8, 8)

# 2D ref on RHS

u = 3
v = 7

@devec r = B[:,:]
@test isequal(r, B[:,:])

@devec r = B[:,2:end]
@test isequal(r, B[:,2:end])

@devec r = B[:,2:5]
@test isequal(r, B[:,2:5])

@devec r = B[u:end,:]
@test isequal(r, B[u:end,:])

@devec r = B[u:end,2:end]
@test isequal(r, B[u:end,2:end])

@devec r = B[u:end,2:5]
@test isequal(r, B[u:end,2:5])

@devec r = B[u:v,:]
@test isequal(r, B[u:v,:])

@devec r = B[u:v,2:end]
@test isequal(r, B[u:v,2:end])

@devec r = B[u:v,2:5]
@test isequal(r, B[u:v,2:5])

@devec r = B[2:end,:]
@test isequal(r, B[2:end,:])

@devec r = B[2:5,:]
@test isequal(r, B[2:5,:])

@devec r = B[:,u:end]
@test isequal(r, B[:,u:end])

@devec r = B[2:end,u:end]
@test isequal(r, B[2:end,u:end])

@devec r = B[2:5,u:end]
@test isequal(r, B[2:5,u:end])

@devec r = B[:,u:v]
@test isequal(r, B[:,u:v])

@devec r = B[2:end,u:v]
@test isequal(r, B[2:end,u:v])

@devec r = B[2:5,u:v]
@test isequal(r, B[2:5,u:v])

# 2D ref on LHS

r = zeros(size(B))
rr = zeros(size(B))
@devec r[:,:] = B[:,:]
rr[:,:] = B[:,:]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[:,2:end] = B[:,2:end]
rr[:,2:end] = B[:,2:end]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[:,2:5] = B[:,2:5]
rr[:,2:5] = B[:,2:5]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[u:end,:] = B[u:end,:]
rr[u:end,:] = B[u:end,:]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[u:end,2:end] = B[u:end,2:end]
rr[u:end,2:end] = B[u:end,2:end]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[u:end,2:5] = B[u:end,2:5]
rr[u:end,2:5] = B[u:end,2:5]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[u:v,:] = B[u:v,:]
rr[u:v,:] = B[u:v,:]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[u:v,2:end] = B[u:v,2:end]
rr[u:v,2:end] = B[u:v,2:end]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[u:v,2:5] = B[u:v,2:5]
rr[u:v,2:5] = B[u:v,2:5]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[2:end,:] = B[2:end,:]
rr[2:end,:] = B[2:end,:]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[2:5,:] = B[2:5,:]
rr[2:5,:] = B[2:5,:]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[:,u:end] = B[:,u:end]
rr[:,u:end] = B[:,u:end]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[2:end,u:end] = B[2:end,u:end]
rr[2:end,u:end] = B[2:end,u:end]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[2:5,u:end] = B[2:5,u:end]
rr[2:5,u:end] = B[2:5,u:end]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[:,u:v] = B[:,u:v]
rr[:,u:v] = B[:,u:v]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[2:end,u:v] = B[2:end,u:v]
rr[2:end,u:v] = B[2:end,u:v]
@test isequal(r, rr)

r = zeros(size(B))
rr = zeros(size(B))
@devec r[2:5,u:v] = B[2:5,u:v]
rr[2:5,u:v] = B[2:5,u:v]
@test isequal(r, rr)

###########################################################
#
#   Interactions of different kinds
#
#   LHS: Ref1D, RefCol, RelRow, Ref2D
#   RHS: Var, Scalar, Ref1D, RefCol, RefRow, Ref2D
#
###########################################################

# LHS is Ref1D

r = zeros(size(a))
@devec r[:] = a
@test isequal(r, a)

r = zeros(size(a))
@devec r[:] = cv
@test isequal(r, fill(cv, size(a)))

r = zeros(size(a))
@devec r[:] = a[:]
@test isequal(r, a)

r = zeros(size(a))
@devec r[:] = abc[:,2]
@test isequal(r, b)

r = zeros(size(a))
@devec r[:] = abct[2,:]
@test isequal(r, b)

# LHS is sub-range of Ref1D

m = 10

r = zeros(m)
rr = zeros(m)
@devec r[2:9] = a
rr[2:9] = a
@test isequal(r, rr)

r = zeros(m)
rr = zeros(m)
@devec r[2:9] = cv
rr[2:9] = cv
@test isequal(r, rr)

r = zeros(m)
rr = zeros(m)
@devec r[2:9] = a[:]
rr[2:9] = a[:]
@test isequal(r, rr)

r = zeros(m)
rr = zeros(m)
@devec r[2:9] = abc[:,2]
rr[2:9] = abc[:,2]
@test isequal(r, rr)

r = zeros(m)
rr = zeros(m)
@devec r[2:9] = abct[2,:]
rr[2:9] = abct[2,:]
@test isequal(r, rr)

# LHS is RefCol

r = zeros(size(abc))
rr = zeros(size(abc))
@devec r[:,2] = a
rr[:,2] = a
@test isequal(r, rr)

r = zeros(size(abc))
rr = zeros(size(abc))
@devec r[:,2] = cv
rr[:,2] = cv
@test isequal(r, rr)

r = zeros(size(abc))
rr = zeros(size(abc))
@devec r[:,2] = a[:]
rr[:,2] = a[:]
@test isequal(r, rr)

r = zeros(size(abc))
rr = zeros(size(abc))
@devec r[:,2] = abc[:,2]
rr[:,2] = abc[:,2]
@test isequal(r, rr)

r = zeros(size(abc))
rr = zeros(size(abc))
@devec r[:,2] = abct[2,:]
rr[:,2] = abct[2,:]
@test isequal(r, rr)

# LHS is RefRow

r = zeros(size(abct))
rr = zeros(size(abct))
@devec r[2,:] = a
rr[2,:] = a
@test isequal(r, rr)

r = zeros(size(abct))
rr = zeros(size(abct))
@devec r[2,:] = cv
rr[2,:] = cv
@test isequal(r, rr)

r = zeros(size(abct))
rr = zeros(size(abct))
@devec r[2,:] = a[:]
rr[2,:] = a[:]
@test isequal(r, rr)

r = zeros(size(abct))
rr = zeros(size(abct))
@devec r[2,:] = abc[:,2]
rr[2,:] = abc[:,2]
@test isequal(r, rr)

r = zeros(size(abct))
rr = zeros(size(abct))
@devec r[2,:] = abct[2,:]
rr[2,:] = abct[2,:]
@test isequal(r, rr)

# LHS is Ref2D

r = zeros(size(abc))
@devec r[:,:] = abc
@test isequal(r, abc)

r = zeros(size(abc))
@devec r[:,:] = cv
@test isequal(r, fill(cv, size(abc)))

r = zeros(size(abc))
@devec r[:,:] = abc[:,:]
@test isequal(r, abc)

# LHS is a subrange of Ref2D

m = 12
n = 5

r = zeros(m, n)
rr = zeros(m, n)
@devec r[3:10, 2:4] = abc
rr[3:10, 2:4] = abc
@test isequal(r, rr)

r = zeros(m, n)
rr = zeros(m, n)
@devec r[3:10, 2:4] = cv
rr[3:10, 2:4] = cv
@test isequal(r, rr)

r = zeros(m, n)
rr = zeros(m, n)
@devec r[3:10, 2:4] = abc[:,:]
rr[3:10, 2:4] = abc[:,:]
@test isequal(r, rr)

# using qualified names

ax = X(a)
@devec r = ax.x[:]
@test isequal(r, a)

r = X(0)
@devec r.x = abc[:,:]
@test isequal(r.x, abc)

r = X(0)
abcx = X(abc)
@devec r.x = abcx.x[:,2]
@test isequal(r.x, b)


text:
# ------------------------------------------------------------------------------------------
# ## Getting started with Jupyter notebooks
#
# Jupyter notebooks are a convenient way to run, display, and present interactive code. The
# main concept is a cell — a single chunk of text. Cells may be markdown (like this one) or
# code (like the next).
#
# ### Running a cell
# To execute code within a cell, select that cell and either (1) hit `Shift` and `Enter` or
# (2) hit the run button (the right pointing arrow) above.
# ------------------------------------------------------------------------------------------

1 + 1
2 + 2

# ------------------------------------------------------------------------------------------
# If you're new to jupyter notebooks, note that only the last line of a cell prints by
# default when you execute that cell and that you can suppress this output with a semicolon
# ------------------------------------------------------------------------------------------

1 + 1
2 + 2;

# ------------------------------------------------------------------------------------------
# ## Changing modes and inserting new cells
#
# You can use the menu above or key combinations — ESC will drop out of editing a cell and
# into a command mode where there are special keyboard shortcuts to insert/cut/modify cells
# themselves. Try `ESC` and then `a` or `b` for "above" and "below".
#
# See the Help menu for all available shortcuts
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Special modes for cells
#
# Cells can change their behavior depending on their very first character
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### How to get docs for Julia functions
#
# To get docs for a function you're not familiar with, precede it with a question mark.
# (This works at the REPL too!)
# ------------------------------------------------------------------------------------------

?println

# ------------------------------------------------------------------------------------------
# ### How to use shell commands
#
# Type `;` and then you can use shell commands. For example,
# ------------------------------------------------------------------------------------------

;ls

;pwd

# ------------------------------------------------------------------------------------------
# ### Interacting with the package manager
#
# Julia's package manager has a special "command" syntax mode — you can enter it with a `]`
# character.
# ------------------------------------------------------------------------------------------

]status

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
export SNP, Insertion, Deletion

abstract Mut

immutable SNP <: Mut
    pos::i32
    ref::Byte
    alt::Byte
end

immutable Insertion <: Mut
    pos::i32
    bases::Bytes
end

immutable Deletion <: Mut
    pos::i32
    bases::Bytes
end

==(x::SNP, y::SNP)             = x.pos==y.pos && x.ref==y.ref && x.alt==y.alt
==(x::Insertion, y::Insertion) = x.pos==y.pos && x.bases==y.bases
==(x::Deletion, y::Deletion)   = x.pos==y.pos && x.bases==y.bases
hash(x::SNP, y::u64)           = hash(x.pos, hash(x.ref, hash(x.alt, y)))
hash(x::Insertion, y::u64)     = hash(x.pos, hash(x.bases, y))
hash(x::Deletion, y::u64)      = hash(x.pos, hash(x.bases, y))
show(io::IO, snp::SNP)         = io << "SNP(" << snp.pos << ":" << snp.ref << "->" << snp.alt << ')'
show(io::IO, indel::Insertion) = io << "Insertion(" << indel.pos << ":" << indel.bases << ')'
show(io::IO, indel::Deletion)  = io << "Deletion(" << indel.pos << ":" << indel.bases << ')'

macro advance_cigar()
    esc(quote
        i += 1
        i > length(cigar) && break
        op, len = cigar[i] & 0x0f, cigar[i] >> 4
    end)
end

substring2byte(s::SubString{String}) = s.string.data[s.offset+1]

function reconstruct_mut_by_md(cigar, md, seq)
    md = matchall(r"\d+|[ATCG]|\^[ATCG]+", md)
    i, j, r, p = 1, 1, 0, 1
    op, len = cigar[i] & 0x0f, cigar[i] >> 4
    r = parse(Int, md[j])
    muts = Mut[]

    while true
        # ‘MIDNSHP=X’→‘012345678’
        if op == 0
            if j % 2 == 0 # SNP
                push!(muts, SNP(p, substring2byte(md[j]), seq[p]))
                len -= 1
                j += 1
                r = parse(Int, md[j])
                p += 1
            else # match
                l = min(len, r)
                len -= l
                r -= l
                p += l
                if r == 0
                    j += 1
                end
            end
            len == 0 && @advance_cigar
        elseif op == 1
            push!(muts, Insertion(p, seq[p:p+len-1]))
            p += len
            @advance_cigar
        elseif op == 2
            s = md[j][2:end] |> String
            push!(muts, Deletion(p-1, s.data))
            @advance_cigar
            j += 1
            r = parse(Int, md[j])
        elseif op == 4 # NOTE: `md` doesn't contains info about softcliped bases
            p += len
            @advance_cigar
        elseif op == 5
            @advance_cigar
        else
            error("TODO: cigar op: $op")
        end
    end

    muts
end

function reconstruct_mut_by_ref()

end