text:
#====================== Auxilliary precomputation =============================#
"""
    EvalInfo

EvalInfo is an auxilliary data structure for storing useful computations for the evaluation of the board.
"""
struct EvalInfo
    wrammedpawns::Bitboard
    brammedpawns::Bitboard
    wmobility::Bitboard
    bmobility::Bitboard
    stage::Int
end

"""
    EvalAttackInfo

EvalAttackInfo is an auxilliary data structure for storing useful computations for the evaluation of the board.
"""
mutable struct EvalAttackInfo
    wpawnattacks::Bitboard
    wknightattacks::Bitboard
    wbishopattacks::Bitboard
    wrookattacks::Bitboard
    wqueenattacks::Bitboard
    wkingattacks::Bitboard
    bpawnattacks::Bitboard
    bknightattacks::Bitboard
    bbishopattacks::Bitboard
    brookattacks::Bitboard
    bqueenattacks::Bitboard
    bkingattacks::Bitboard
    passed::Bitboard
end

function initEvalInfo(board::Board)
    # extract pawn positions
    wpawns = white(board) & pawns(board)
    bpawns = black(board) & pawns(board)

    # check for rammed pawns
    wrammedpawns = pawnAdvance(bpawns, wpawns, BLACK)
    brammedpawns = pawnAdvance(wpawns, bpawns, WHITE)

    # generate list of all positions attacked by a pawn
    wpawnattacks = pawnCapturesWhite(wpawns, FULL)
    bpawnattacks = pawnCapturesBlack(bpawns, FULL)

    # generate list of all positions attacked by a knight
    wknightattacks = bknightattacks = EMPTY
    wbishopattacks = bbishopattacks = EMPTY
    wrookattacks = brookattacks = EMPTY
    wqueenattacks = bqueenattacks = EMPTY
    wkingattacks = bkingattacks = EMPTY

    # mobility regions
    wmobility = ~((white(board) & kings(board)) | bpawnattacks)
    bmobility = ~((black(board) & kings(board)) | wpawnattacks)

    passed = EMPTY

    gamestage = stage(board)

    ei = EvalInfo(wrammedpawns, brammedpawns, wmobility, bmobility, gamestage)
    ea = EvalAttackInfo(wpawnattacks, wknightattacks, wbishopattacks, wrookattacks, wqueenattacks, wkingattacks, bpawnattacks, bknightattacks, bbishopattacks, brookattacks, bqueenattacks, bkingattacks, passed)
    ei, ea
end

#============================== Game stage functions ==========================#

function stage(board::Board)
    stage = 6156 - 1024count(queens(board)) - 512count(rooks(board)) - 256count(knights(board) | bishops(board))
    stage = fld(stage, 24)
end

#=============================== Endgame scoring scaling ======================#

function scale_factor(board::Board, eval::Int)
    if isone(board[WHITEBISHOP]) && isone(board[BLACKBISHOP]) && isone(bishops(board) & LIGHT)
        if isempty(knights(board) | rooklike(board))
            return SCALE_OCB_BISHOPS
        end
        if isempty(rooklike(board)) && isone(board[WHITEKNIGHT]) && isone(board[BLACKKNIGHT])
            return SCALE_OCB_ONE_KNIGHT
        end
        if isempty(knights(board) | queens(board)) && isone(board[WHITEROOK]) && isone(board[BLACKROOK])
            return SCALE_OCB_ONE_ROOK
        end
    end
    if (eval > 0) && (count(white(board)) == 2) && ismany(board[WHITEKNIGHT] | board[WHITEBISHOP])
        return SCALE_DRAW
    elseif (eval < 0) && (count(black(board)) == 2) && ismany(board[BLACKKNIGHT] | board[BLACKBISHOP])
        return SCALE_DRAW
    end
    return SCALE_NORMAL
end

#=============================== Main evaluation ==============================#

"""
    evaluate(board)

Naive evaluation function to get the code development going.
"""
function evaluate(board::Board, pktable::PawnKingTable)
    ei, ea = initEvalInfo(board)

    score = 0
    score += board.psqteval

    if (pkt_entry = getPKentry(pktable, board.pkhash)) !== PKT_BLANK
        score += pkt_entry.score
        ea.passed = pkt_entry.passed
    else
        kpscore = 0
        kpscore += evaluate_pawns(board, ea)
        kpscore += evaluate_kingpawns(board, ea)
        storePKentry!(pktable, board.pkhash, ea.passed, kpscore)
        score += kpscore
    end

    v = fld(scoreEG(score) + scoreMG(score), 2)
    v = round(Int, v * (100.0 - Float64(board.halfmovecount)) / 100.0)
    if abs(v) > LAZY_THRESH
        if board.turn == WHITE
            return v
        else
            return -v
        end
    end

    score += evaluate_knights(board, ei, ea)
    score += evaluate_bishops(board, ei, ea)
    score += evaluate_rooks(board, ei, ea)
    score += evaluate_queens(board, ei, ea)
    score += evaluate_kings(board, ea)
    score += evaluate_pins(board)
    score += evaluate_space(board, ea)
    score += evaluate_threats(board, ei, ea)

    score += evaluate_passed(board, ea)

    score += evaluate_initiative(board, ea.passed, score)

    scale_f = scale_factor(board, scoreEG(score))

    eval = (256 - ei.stage) * scoreMG(score) + ei.stage * scoreEG(score) * fld(scale_f, SCALE_NORMAL)
    eval = fld(eval, 256)

    if board.turn == WHITE
        eval += TEMPO_BONUS#
        eval = round(Int, eval * (100.0 - Float64(board.halfmovecount)) / 100.0)
        return eval
    else
        eval -= TEMPO_BONUS
        eval = round(Int, eval * (100.0 - Float64(board.halfmovecount)) / 100.0)
        return -eval
    end
end

function evaluate_pawns(board::Board, ea::EvalAttackInfo)

    w_pawns = white(board) & pawns(board)
    b_pawns = black(board) & pawns(board)

    w_king = square(white(board) & kings(board))
    b_king = square(black(board) & kings(board))

    score = 0

    @inbounds for pawn in w_pawns
        file = fileof(pawn)
        rank = rankof(pawn)
        blocked    = b_pawns & (Bitboard(pawn) << 8)
        stoppers   = b_pawns & PASSED_PAWN_MASKS[1][pawn]
        opposed    = stoppers & FILE[file]
        lever      = b_pawns & PAWN_CAPTURES_WHITE[pawn]
        leverpush  = b_pawns & PAWN_CAPTURES_WHITE[pawn + 8]
        doubled    = !isempty(w_pawns & (Bitboard(pawn) >> 8))
        neighbours = w_pawns & NEIGHBOUR_FILE_MASKS[file]
        phalanx    = neighbours & RANK[rank]
        support    = neighbours & RANK[rank - 1]

        # Passed pawns
        if isempty(stoppers)
            ea.passed |= Bitboard(pawn)
        end

        # Backward pawns
        backward = isempty(neighbours & PASSED_PAWN_MASKS[2][pawn - 8] & FILE[file]) && !isempty(leverpush | blocked)

        # file mapping
        file_psqt = FILE_TO_QSIDE_MAP[file]

        # Connected pawns
        if !isempty(support) || !isempty(phalanx)
            score += CONNECTED_PAWN_PSQT[rank][file_psqt]

        # Isolated pawns
        elseif isempty(neighbours)
            score -= ISOLATED_PAWN_PENALTY + WEAK_UNOPPOSED * isempty(opposed)

        # Backward pawns
        elseif backward
            score -= BACKWARD_PAWN + WEAK_UNOPPOSED * isempty(opposed)
        end

        # Doubled pawns
        if isempty(support)
            score -= DOUBLE_PAWN_PENALTY * doubled + WEAK_LEVER * ismany(lever)
        end
    end

    @inbounds for pawn in b_pawns
        file = fileof(pawn)
        rank = rankof(pawn)
        blocked    = w_pawns & (Bitboard(pawn) >> 8)
        stoppers   = w_pawns & PASSED_PAWN_MASKS[2][pawn]
        opposed    = stoppers & FILE[file]
        lever      = w_pawns & PAWN_CAPTURES_BLACK[pawn]
        leverpush  = w_pawns & PAWN_CAPTURES_BLACK[pawn - 8]
        doubled    = !isempty(b_pawns & (Bitboard(pawn) << 8))
        neighbours = b_pawns & NEIGHBOUR_FILE_MASKS[file]
        phalanx    = neighbours & RANK[rank]
        support    = neighbours & RANK[rank + 1]

        # Passed pawns
        if isempty(stoppers)
            ea.passed |= Bitboard(pawn)
        end

        # Backward pawns
        backward = isempty(neighbours & PASSED_PAWN_MASKS[1][pawn + 8] & FILE[file]) && !isempty(leverpush | blocked)

        # file mapping
        file_psqt = FILE_TO_QSIDE_MAP[file]

        # Connected pawns
        if !isempty(support) || !isempty(phalanx)
            score -= CONNECTED_PAWN_PSQT[9 - rank][file_psqt]

        # Isolated pawns
        elseif isempty(neighbours)
            score += ISOLATED_PAWN_PENALTY + WEAK_UNOPPOSED * isempty(opposed)

        # Backward pawns
        elseif backward
            score += BACKWARD_PAWN + WEAK_UNOPPOSED * isempty(opposed)
        end

        # Doubled pawns
        if isempty(support)
            score += DOUBLE_PAWN_PENALTY * doubled + WEAK_LEVER * ismany(lever)
        end
    end

    if isempty(pawns(board) & KINGFLANK[fileof(w_king)])
        score -= PAWNLESS_FLANK
    end
    if isempty(pawns(board) & KINGFLANK[fileof(b_king)])
        score += PAWNLESS_FLANK
    end

    return score
end

function evaluate_knights(board::Board, ei::EvalInfo, ea::EvalAttackInfo)
    w_knights = white(board) & knights(board)
    b_knights = black(board) & knights(board)

    w_pawns = pawns(board) & white(board)
    b_pawns = pawns(board) & black(board)

    w_king_sqr = square(kings(board) & white(board))
    b_king_sqr = square(kings(board) & black(board))

    score = 0

    w_outposts = (RANK_4 | RANK_5 | RANK_6) & ea.wpawnattacks & ~ea.bpawnattacks
    @inbounds for knight in w_knights
        attacks = knightMoves(knight)
        ea.wknightattacks |= attacks
        score += KNIGHT_MOBILITY[count(attacks & ei.wmobility) + 1]
        score -= MINOR_KING_PROTECTION * DISTANCE_BETWEEN[w_king_sqr, knight]
    end

    # Score knight outposts
    score += count(w_outposts & w_knights) * KNIGHT_OUTPOST_BONUS

    # Score reachable knight outposts
    score += count(w_outposts & ea.wknightattacks & ~white(board)) * KNIGHT_POTENTIAL_OUTPOST_BONUS

    # Score knights behind pawns
    score += count((w_knights << 8) & w_pawns) * PAWN_SHIELD_BONUS

    b_outposts = (RANK_3 | RANK_4 | RANK_5) & ~ea.wpawnattacks & ea.bpawnattacks
    @inbounds for knight in b_knights
        attacks = knightMoves(knight)
        ea.bknightattacks |= attacks
        score -= KNIGHT_MOBILITY[count(attacks & ei.bmobility) + 1]
        score += MINOR_KING_PROTECTION * DISTANCE_BETWEEN[b_king_sqr, knight]
    end

    # Score knight outposts
    score -= count(b_outposts & b_knights) * KNIGHT_OUTPOST_BONUS

    # Score reachable knight outposts
    score -= count(b_outposts & ea.bknightattacks & ~black(board)) * KNIGHT_POTENTIAL_OUTPOST_BONUS

    # Score knights behind pawns
    score -= count((b_knights >> 8) & b_pawns) * PAWN_SHIELD_BONUS

    # bonus for knights in rammed positions
    num2_rammed = count(ei.wrammedpawns)^2
    score += div(count(w_knights) * num2_rammed, 4) * KNIGHT_RAMMED_BONUS
    score -= div(count(b_knights) * num2_rammed, 4) * KNIGHT_RAMMED_BONUS

    score
end

function evaluate_bishops(board::Board, ei::EvalInfo, ea::EvalAttackInfo)
    w_bishops = white(board) & bishops(board)
    b_bishops = black(board) & bishops(board)

    w_pawns = pawns(board) & white(board)
    b_pawns = pawns(board) & black(board)

    w_king_sqr = square(kings(board) & white(board))
    b_king_sqr = square(kings(board) & black(board))

    score = 0

    attacks = EMPTY
    occ = occupied(board)

    w_outposts = (RANK_4 | RANK_5 | RANK_6) & ea.wpawnattacks & ~ea.bpawnattacks
    @inbounds for bishop in w_bishops
        attacks = bishopMoves(bishop, occ)
        ea.wbishopattacks |= attacks
        score += BISHOP_MOBILITY[count(attacks & ei.wmobility) + 1]
        if ismany(attacks & CENTRAL_SQUARES)
            score += BISHOP_CENTRAL_CONTROL
        end
        score -= MINOR_KING_PROTECTION * DISTANCE_BETWEEN[w_king_sqr, bishop]
    end

    # Outpost bonus
    score += count(w_outposts & w_bishops) * BISHOP_OUTPOST_BONUS

    # Add a bonus for being behind a pawn.
    score += count((w_bishops << 8) & w_pawns) * PAWN_SHIELD_BONUS

    b_outposts = (RANK_3 | RANK_4 | RANK_5) & ~ea.wpawnattacks & ea.bpawnattacks
    @inbounds for bishop in b_bishops
        attacks = bishopMoves(bishop, occ)
        ea.bbishopattacks |= attacks
        score -= BISHOP_MOBILITY[count(attacks & ei.bmobility) + 1]
        if ismany(attacks & CENTRAL_SQUARES)
            score -= BISHOP_CENTRAL_CONTROL
        end
        score += MINOR_KING_PROTECTION * DISTANCE_BETWEEN[b_king_sqr, bishop]
    end

    # Outpost bonus
    score -= count(b_outposts & b_bishops) * BISHOP_OUTPOST_BONUS

    # Add a bonus for being behind a pawn
    score -= count((b_bishops >> 8) & b_pawns) * PAWN_SHIELD_BONUS

    # bishop pair
    if count(w_bishops) >= 2
        score += BISHOP_PAIR_BONUS
    end
    if count(b_bishops) >= 2
        score -= BISHOP_PAIR_BONUS
    end

    # penalty for bishops on colour of own pawns
    if !isempty(w_bishops & LIGHT)
        score -= BISHOP_COLOR_PENALTY * count(w_pawns & LIGHT)
        score -= BISHOP_RAMMED_COLOR_PENALTY * count(ei.wrammedpawns & LIGHT)
    end
    if !isempty(w_bishops & DARK)
        score -= BISHOP_COLOR_PENALTY * count(w_pawns & DARK)
        score -= BISHOP_RAMMED_COLOR_PENALTY * count(ei.wrammedpawns & DARK)
    end
    if !isempty(b_bishops & LIGHT)
        score += BISHOP_COLOR_PENALTY * count(b_pawns & LIGHT)
        score += BISHOP_RAMMED_COLOR_PENALTY * count(ei.brammedpawns & LIGHT)
    end
    if !isempty(b_bishops & DARK)
        score += BISHOP_COLOR_PENALTY * count(b_pawns & DARK)
        score += BISHOP_RAMMED_COLOR_PENALTY * count(ei.brammedpawns & DARK)
    end

    score
end

function evaluate_rooks(board::Board, ei::EvalInfo, ea::EvalAttackInfo)
    w_rooks = (white(board) & rooks(board))
    b_rooks = (black(board) & rooks(board))

    score = 0

    occ = occupied(board)
    @inbounds for rook in w_rooks
        rfile = file(rook)
        if isempty(rfile & pawns(board))
            score += ROOK_OPEN_FILE_BONUS
        elseif isempty(rfile & white(board) & pawns(board))
            score += ROOK_SEMIOPEN_FILE_BONUS
        end
        if !isempty(rfile & queens(board))
            score += ROOK_ON_QUEEN_FILE
        end
        attacks = rookMoves(rook, occ)
        ea.wrookattacks |= attacks
        mob_cnt = count(attacks & ei.wmobility) + 1
        score += ROOK_MOBILITY[mob_cnt]
    end

    @inbounds for rook in b_rooks
        rfile = file(rook)
        if isempty(rfile & pawns(board))
            score -= ROOK_OPEN_FILE_BONUS
        elseif isempty(rfile & black(board) & pawns(board))
            score -= ROOK_SEMIOPEN_FILE_BONUS
        end
        if !isempty(rfile & queens(board))
            score -= ROOK_ON_QUEEN_FILE
        end
        attacks = rookMoves(rook, occ)
        ea.brookattacks |= attacks
        mob_cnt = count(attacks & ei.bmobility) + 1
        score -= ROOK_MOBILITY[mob_cnt]
    end

    score
end

function evaluate_queens(board::Board, ei::EvalInfo, ea::EvalAttackInfo)
    w_queens = board[WHITEQUEEN]
    b_queens = board[BLACKQUEEN]

    score = 0

    occ = occupied(board)

    @inbounds for queen in w_queens
        attacks = queenMoves(queen, occ)
        ea.wqueenattacks |= attacks
        score += QUEEN_MOBILITY[count(attacks & ei.wmobility) + 1]
    end

    @inbounds for queen in b_queens
        attacks = queenMoves(queen, occ)
        ea.bqueenattacks |= attacks
        score -= QUEEN_MOBILITY[count(attacks & ei.bmobility) + 1]
    end

    score
end

function evaluate_kings(board::Board, ea::EvalAttackInfo)
    w_king = white(board) & kings(board)
    b_king = black(board) & kings(board)
    w_king_sqr = square(w_king)
    b_king_sqr = square(b_king)
    ea.wkingattacks |= kingMoves(w_king_sqr)
    ea.bkingattacks |= kingMoves(b_king_sqr)

    score = 0

    if cancastlekingside(board, WHITE)
        score += CASTLE_OPTION_BONUS
    end
    if cancastlequeenside(board, WHITE)
        score += CASTLE_OPTION_BONUS
    end
    if cancastlekingside(board, BLACK)
        score -= CASTLE_OPTION_BONUS
    end
    if cancastlequeenside(board, BLACK)
        score -= CASTLE_OPTION_BONUS
    end

    score
end

function evaluate_kingpawns(board::Board, ea::EvalAttackInfo)
    w_king = white(board) & kings(board)
    b_king = black(board) & kings(board)
    w_king_sqr = square(w_king)
    b_king_sqr = square(b_king)

    score = 0

    wkingfile = fileof(w_king_sqr)
    lower = max(1, wkingfile - 1)
    upper = min(8, wkingfile + 1)
    _wpawns = pawns(board) & white(board) & PASSED_PAWN_MASKS[1][w_king_sqr]
    _bpawns = pawns(board) & black(board) & PASSED_PAWN_MASKS[1][w_king_sqr]
    wksqr = rankof(w_king_sqr)
    @inbounds for file in lower:upper
        wpawns = _wpawns & FILE[file]
        if isempty(wpawns)
            w_dist = 8
        else
            w_dist = abs(wksqr - rankof(trailing_zeros(wpawns.val)+1))
        end

        bpawns = _bpawns & FILE[file]
        if isempty(bpawns)
            b_dist = 8
        else
            b_dist = abs(wksqr - rankof(trailing_zeros(bpawns.val)+1))
        end

        if file === wkingfile
            score += KING_SHELTER_ONFILE[file][w_dist]
        else
            score += KING_SHELTER_OFFFILE[file][w_dist]
        end

        blocked = (b_dist !== 8) && (w_dist === b_dist - 1)
        fmap = FILE_TO_QSIDE_MAP[file]
        if blocked
            score += KING_STORM_BLOCKED[fmap][b_dist]
        else
            score += KING_STORM_UNBLOCKED[fmap][b_dist]
        end
    end

    bkingfile = fileof(b_king_sqr)
    lower = max(1, bkingfile - 1)
    upper = min(8, bkingfile + 1)
    _bpawns = pawns(board) & black(board) & PASSED_PAWN_MASKS[2][b_king_sqr]
    _wpawns = pawns(board) & white(board) & PASSED_PAWN_MASKS[2][b_king_sqr]
    bksqr = rankof(b_king_sqr)
    @inbounds for file in lower:upper
        bpawns = _bpawns & FILE[file]
        if isempty(bpawns)
            b_dist = 8
        else
            b_dist = abs(bksqr - rankof(64-leading_zeros(bpawns.val)))
        end

        wpawns = _wpawns & FILE[file]
        if isempty(wpawns)
            w_dist = 8
        else
            w_dist = abs(bksqr - rankof(64-leading_zeros(wpawns.val)))
        end

        if file === bkingfile
            score -= KING_SHELTER_ONFILE[file][b_dist]
        else
            score -= KING_SHELTER_OFFFILE[file][b_dist]
        end

        blocked = (w_dist !== 8) && (b_dist === w_dist - 1)
        fmap = FILE_TO_QSIDE_MAP[file]
        if blocked
            score -= KING_STORM_BLOCKED[fmap][w_dist]
        else
            score -= KING_STORM_UNBLOCKED[fmap][w_dist]
        end
    end

    score
end

function evaluate_pins(board::Board)
    eval = 0

    # switch turn and find all pins
    board.turn = !board.turn
    opp_pinned = findpins(board)
    board.turn = !board.turn

    if board.turn == WHITE
        eval -= count(board.pinned)
        eval += count(opp_pinned)
    else
        eval += count(board.pinned)
        eval -= count(opp_pinned)
    end

    eval *= PIN_BONUS

    # specific additional pin bonus
    eval -= count(pinned(board) & queens(board)) * 30
    eval += count(opp_pinned & queens(board)) * 30
    eval -= count(pinned(board) & rooks(board)) * 10
    eval += count(opp_pinned & rooks(board)) * 10

    makescore(eval, eval)
end

function evaluate_space(board::Board, ea::EvalAttackInfo)
    eval = 0

    w_filter = (RANK_2 | RANK_3 | RANK_4) & CENTERFILES
    b_filter = (RANK_5 | RANK_6 | RANK_7) & CENTERFILES
    w_attacks = ea.wpawnattacks | ea.wknightattacks | ea.wbishopattacks | ea.wrookattacks | ea.wqueenattacks | ea.wkingattacks
    b_attacks = ea.bpawnattacks | ea.bknightattacks | ea.bbishopattacks | ea.brookattacks | ea.bqueenattacks | ea.bkingattacks
    eval += count(w_filter & ~b_attacks)
    eval -= count(b_filter & ~w_attacks)
    eval *= SPACE_BONUS

    makescore(eval, eval)
end

function evaluate_threats(board::Board, ei::EvalInfo, ea::EvalAttackInfo)
    w_attacks = ea.wpawnattacks
    w_double_attacks = w_attacks & ea.wknightattacks
    w_attacks |= ea.wknightattacks
    w_double_attacks |= w_attacks & ea.wbishopattacks
    w_attacks |= ea.wbishopattacks
    w_double_attacks |= w_attacks & ea.wrookattacks
    w_attacks |= ea.wrookattacks
    w_double_attacks |= w_attacks & ea.wqueenattacks
    w_attacks |= ea.wqueenattacks
    w_double_attacks |= w_attacks & ea.wkingattacks
    w_attacks |= ea.wkingattacks

    b_attacks = ea.bpawnattacks
    b_double_attacks = b_attacks & ea.bknightattacks
    b_attacks |= ea.bknightattacks
    b_double_attacks |= b_attacks & ea.bbishopattacks
    b_attacks |= ea.bbishopattacks
    b_double_attacks |= b_attacks & ea.brookattacks
    b_attacks |= ea.brookattacks
    b_double_attacks |= b_attacks & ea.bqueenattacks
    b_attacks |= ea.bqueenattacks
    b_double_attacks |= b_attacks & ea.bkingattacks
    b_attacks |= ea.bkingattacks

    occ = occupied(board)

    score = 0

    #========================= Evaluation w.r.t. white ========================#
    # strongly protected by the enemy.
    strongly_protected = ea.bpawnattacks | (b_double_attacks & ~w_double_attacks)
    # well defended by the enemy
    defended = (black(board) & ~pawns(board)) & strongly_protected
    # not well defended by the enemy
    weak = black(board) & ~strongly_protected & w_attacks

    # Case where our opponent is hanging pieces
    case = ~b_attacks | ((black(board) & ~pawns(board)) & w_double_attacks)
    # Bonus if opponent is hanging pieces
    score += HANGING_BONUS * count(weak & case)

    if !isempty(weak & ea.wkingattacks)
        score += THREAT_BY_KING
    end

    # Case where our opponent is defended or weak, and attacked by a bishop or knight.
    case = (defended | weak) & (ea.wknightattacks | ea.wbishopattacks)
    case_pawns = count(case & pawns(board))
    case_knights = count(case & knights(board))
    case_bishops = count(case & bishops(board))
    case_rooks = count(case & rooks(board))
    case_queens = count(case & queens(board))
    @inbounds score += THREAT_BY_MINOR[1] * case_pawns
    @inbounds score += THREAT_BY_MINOR[2] * case_knights
    @inbounds score += THREAT_BY_MINOR[3] * case_bishops
    @inbounds score += THREAT_BY_MINOR[4] * case_rooks
    @inbounds score += THREAT_BY_MINOR[5] * case_queens

    # Case where our opponent is weak and attacked by our rook
    case = weak & ea.wrookattacks
    case_pawns = count(case & pawns(board))
    case_knights = count(case & knights(board))
    case_bishops = count(case & bishops(board))
    case_rooks = count(case & rooks(board))
    case_queens = count(case & queens(board))
    @inbounds score += THREAT_BY_ROOK[1] * case_pawns
    @inbounds score += THREAT_BY_ROOK[2] * case_knights
    @inbounds score += THREAT_BY_ROOK[3] * case_bishops
    @inbounds score += THREAT_BY_ROOK[4] * case_rooks
    @inbounds score += THREAT_BY_ROOK[5] * case_queens

    # Case where we get a bonus for restricting our opponent
    case = b_attacks & ~strongly_protected & w_attacks
    score += count(case) * RESTRICTION_BONUS

    # Find squares which are safe
    safe = ~b_attacks | w_attacks

    # Case where we have safe pawn attacks
    case = pawns(board) & white(board) & safe
    case = pawnCapturesWhite(case, black(board) & ~pawns(board))
    score += THREAT_BY_PAWN * count(case)

    # Case threat by pawn push
    case = pawnAdvance(white(board) & pawns(board), empty(board), WHITE)
    case |= pawnDoubleAdvance(white(board) & pawns(board), empty(board), WHITE)
    case &= ~ea.bpawnattacks & safe
    case = pawnCapturesWhite(case, black(board) & ~pawns(board))
    score += THREAT_BY_PUSH * count(case)

    # Evaluate if there are pieces that can reach a checking square, while being safe.
    # King Danger evaluations.
    # As WHITE, we are seeing how in danger the BLACK king is.
    b_king_sqr = square(black(board) & kings(board))
    b_king_box = ea.bkingattacks | Bitboard(b_king_sqr)
    safe = ~white(board) & (~b_attacks | (weak & w_double_attacks))
    king_bishop_moves = bishopMoves(b_king_sqr, occ)
    king_rook_moves = rookMoves(b_king_sqr, occ)
    king_knight_moves = knightMoves(b_king_sqr)
    knightcheck_sqrs = king_knight_moves & safe & ea.wknightattacks
    rookcheck_sqrs   = king_rook_moves & safe & ea.wrookattacks
    # Remove rook checks from queen checks, because they are more valuable.
    # Don't count if queens trade.
    queencheck_sqrs  = (king_rook_moves | king_bishop_moves) & safe & ea.wqueenattacks & ~rookcheck_sqrs & ~ea.bqueenattacks
    # Remove queen checks from bishop checks, because they are more valuable.
    bishopcheck_sqrs = king_bishop_moves & safe & ea.wbishopattacks & ~queencheck_sqrs

    king_danger = 0
    unsafechecks = EMPTY
    if !isempty(knightcheck_sqrs)
        king_danger += KNIGHT_SAFE_CHECK
    else
        unsafechecks |= king_knight_moves & ea.wknightattacks
    end
    if !isempty(bishopcheck_sqrs)
        king_danger += BISHOP_SAFE_CHECK
    else
        unsafechecks |= king_bishop_moves & ea.wbishopattacks
    end
    if !isempty(rookcheck_sqrs)
        king_danger += ROOK_SAFE_CHECK
    else
        unsafechecks |= king_rook_moves & ea.wrookattacks
    end
    if !isempty(queencheck_sqrs)
        king_danger += QUEEN_SAFE_CHECK
    end
    if isempty(board[WHITEQUEEN])
        king_danger -= KD_QUEEN
    end
    if isempty(board[WHITEROOK])
        king_danger -= KD_ROOK
    end
    if isempty(board[WHITEBISHOP])
        king_danger -= KD_BISHOP
    end
    if isempty(board[WHITEKNIGHT])
        king_danger -= KD_KNIGHT
    end
    king_danger += 9 * count(w_attacks & ea.bkingattacks)
    if king_danger > 0
        score += makescore(king_danger, fld(king_danger, 2))
    end

    # Score the number of attacks on the king's flank
    camp = FULL ⊻ RANK_1 ⊻ RANK_2 ⊻ RANK_3
    black_king_flank = KINGFLANK[fileof(b_king_sqr)]
    king_flank_attacks = w_attacks & black_king_flank & camp
    king_flank_double_attacks = king_flank_attacks & w_double_attacks
    score += (count(king_flank_attacks) + count(king_flank_double_attacks)) * KING_FLANK_ATTACK
    king_flank_defence = b_attacks & black_king_flank & camp
    score -= count(king_flank_defence) * KING_FLANK_DEFEND

    # King ring attack bonus
    score += count(b_king_box & weak)^2 * KING_BOX_WEAK
    score += count(unsafechecks)^2 * UNSAFE_CHECK

    # Threats on black queen
    black_queens = black(board) & queens(board)
    if isone(black_queens)
        safe = ei.wmobility & ~strongly_protected
        case = ea.wknightattacks & knightMoves(square(black_queens)) & safe
        score += count(case) * KNIGHT_ON_QUEEN
    end

    #========================= Evaluation w.r.t. black ========================#

    # strongly protected by the enemy.
    strongly_protected = ea.wpawnattacks | (w_double_attacks & ~b_double_attacks)
    # well defended by the enemy
    defended = (white(board) & ~pawns(board)) & strongly_protected
    # not well defended by the enemy
    weak = white(board) & ~strongly_protected & b_attacks

    # Case where our opponent is hanging pieces
    case = ~w_attacks | ((white(board) & ~pawns(board)) & b_double_attacks)
    # Bonus if opponent is hanging pieces
    score -= HANGING_BONUS * count(weak & case)

    if !isempty(weak & ea.bkingattacks)
        score -= THREAT_BY_KING
    end

    # Case where our opponent is defended or weak, and attacked by a bishop or knight.
    case = (defended | weak) & (ea.bknightattacks | ea.bbishopattacks)
    case_pawns = count(case & pawns(board))
    case_knights = count(case & knights(board))
    case_bishops = count(case & bishops(board))
    case_rooks = count(case & rooks(board))
    case_queens = count(case & queens(board))
    @inbounds score -= THREAT_BY_MINOR[1] * case_pawns
    @inbounds score -= THREAT_BY_MINOR[2] * case_knights
    @inbounds score -= THREAT_BY_MINOR[3] * case_bishops
    @inbounds score -= THREAT_BY_MINOR[4] * case_rooks
    @inbounds score -= THREAT_BY_MINOR[5] * case_queens

    # Case where our opponent is weak and attacked by our rook
    case = weak & ea.brookattacks
    case_pawns = count(case & pawns(board))
    case_knights = count(case & knights(board))
    case_bishops = count(case & bishops(board))
    case_rooks = count(case & rooks(board))
    case_queens = count(case & queens(board))
    @inbounds score -= THREAT_BY_ROOK[1] * case_pawns
    @inbounds score -= THREAT_BY_ROOK[2] * case_knights
    @inbounds score -= THREAT_BY_ROOK[3] * case_bishops
    @inbounds score -= THREAT_BY_ROOK[4] * case_rooks
    @inbounds score -= THREAT_BY_ROOK[5] * case_queens

    # Case where we get a bonus for restricting our opponent
    case = b_attacks & ~strongly_protected & w_attacks
    score -= count(case) * RESTRICTION_BONUS

    # Case threat by pawn
    safe = ~w_attacks | b_attacks
    case = pawns(board) & black(board) & safe
    case = pawnCapturesBlack(case, white(board) & ~pawns(board))
    score -= THREAT_BY_PAWN * count(case)

    # Case threat by pawn push
    case = pawnAdvance(black(board) & pawns(board), empty(board), BLACK)
    case |= pawnDoubleAdvance(black(board) & pawns(board), empty(board), BLACK)
    case &= ~ea.wpawnattacks & safe
    case = pawnCapturesBlack(case, white(board) & ~pawns(board))
    score -= THREAT_BY_PUSH * count(case)

    # Evaluate if there are pieces that can reach a checking square, while being safe.
    # King Danger evaluations.
    # As BLACK, we are seeing how in danger the WHITE king is.
    w_king_sqr = square(white(board) & kings(board))
    w_king_box = ea.wkingattacks | Bitboard(w_king_sqr)
    safe = ~black(board) & (~w_attacks | (weak & b_double_attacks))
    king_bishop_moves = bishopMoves(w_king_sqr, occ)
    king_rook_moves = rookMoves(w_king_sqr, occ)
    king_knight_moves = knightMoves(w_king_sqr)
    knightcheck_sqrs = king_knight_moves & safe & ea.bknightattacks
    rookcheck_sqrs   = king_rook_moves & safe & ea.brookattacks
    # Remove rook checks from queen checks, because they are more valuable.
    # Don't count if queens trade.
    queencheck_sqrs  = (king_rook_moves | king_bishop_moves) & safe & ea.bqueenattacks & ~rookcheck_sqrs & ~ea.wqueenattacks
    # Remove queen checks from bishop checks, because they are more valuable.
    bishopcheck_sqrs = king_bishop_moves & safe & ea.bbishopattacks & ~queencheck_sqrs

    king_danger = 0
    unsafechecks = EMPTY
    if !isempty(knightcheck_sqrs)
        king_danger += KNIGHT_SAFE_CHECK
    else
        unsafechecks |= king_knight_moves & ea.bknightattacks
    end
    if !isempty(bishopcheck_sqrs)
        king_danger += BISHOP_SAFE_CHECK
    else
        unsafechecks |= king_bishop_moves & ea.bbishopattacks
    end
    if !isempty(rookcheck_sqrs)
        king_danger += ROOK_SAFE_CHECK
    else
        unsafechecks |= king_rook_moves & ea.brookattacks
    end
    if !isempty(queencheck_sqrs)
        king_danger += QUEEN_SAFE_CHECK
    end
    if isempty(board[BLACKQUEEN])
        king_danger -= KD_QUEEN
    end
    if isempty(board[BLACKROOK])
        king_danger -= KD_ROOK
    end
    if isempty(board[BLACKBISHOP])
        king_danger -= KD_BISHOP
    end
    if isempty(board[BLACKKNIGHT])
        king_danger -= KD_KNIGHT
    end
    king_danger += 9 * count(b_attacks & ea.wkingattacks)
    if king_danger > 0
        score -= makescore(king_danger, fld(king_danger, 2))
    end

    # Count king flank attacks
    camp = (FULL ⊻ RANK_8 ⊻ RANK_7 ⊻ RANK_6)
    white_king_flank = KINGFLANK[fileof(w_king_sqr)]
    king_flank_attacks = b_attacks & white_king_flank & camp
    king_flank_double_attacks = king_flank_attacks & b_double_attacks
    score -= (count(king_flank_attacks) + count(king_flank_double_attacks)) * KING_FLANK_ATTACK
    king_flank_defence = w_attacks & white_king_flank & camp
    score += count(king_flank_defence) * KING_FLANK_DEFEND

    # King ring attack bonus
    score -= count(w_king_box & weak)^2 * KING_BOX_WEAK
    score -= count(unsafechecks)^2 * UNSAFE_CHECK

    # Attacks on white queen
    white_queens = white(board) & queens(board)
    if isone(white_queens)
        safe = ei.bmobility & ~strongly_protected
        case = ea.bknightattacks & knightMoves(square(white_queens)) & safe
        score -= count(case) * KNIGHT_ON_QUEEN
    end

    score
end

function evaluate_initiative(board::Board, passed::Bitboard, score::Int)
    eg = scoreEG(score)
    mg = scoreMG(score)
    black_king = square(kings(board) & black(board))
    white_king = square(kings(board) & white(board))
    black_king_rank = rankof(black_king)
    white_king_rank = rankof(white_king)
    black_king_file = fileof(black_king)
    white_king_file = fileof(white_king)

    outflank = abs(white_king_file - black_king_file) - abs(white_king_rank - black_king_rank)

    infiltrating = (white_king_rank > 4) || (black_king_rank < 5)

    pawns_both_flanks = !isempty(pawns(board) & KINGSIDE) || !isempty(pawns(board) & QUEENSIDE)

    slimchance = isempty(passed) && (outflank < 0) && !pawns_both_flanks

    complexity = 8 * count(passed) + 10 * count(pawns(board)) + 8 * outflank +
        (infiltrating ? 11 : 0) + (pawns_both_flanks ? 19 : 0) + (has_non_pawn_material(board) ? 0 : 44) -
        (slimchance ? 39 : 0) - 90

    c_mg = (Int(mg > 0) - Int(mg < 0)) * clamp(complexity+45, -abs(mg), 0)
    c_eg = (Int(eg > 0) - Int(eg < 0)) * max(complexity, -abs(eg))

    makescore(c_mg, c_eg)
end

function evaluate_passed(board::Board, ea::EvalAttackInfo)
    black_king = square(kings(board) & black(board))
    white_king = square(kings(board) & white(board))
    white_pawns = white(board) & ea.passed
    black_pawns = black(board) & ea.passed

    w_attacks = ea.wknightattacks | ea.wpawnattacks | ea.wbishopattacks | ea.wrookattacks | ea.wqueenattacks | ea.wkingattacks
    b_attacks = ea.bknightattacks | ea.bpawnattacks | ea.bbishopattacks | ea.brookattacks | ea.bqueenattacks | ea.bkingattacks

    block_sqr = 0

    score = 0
    @inbounds for pawn in white_pawns
        bonus = 0
        file = fileof(pawn)
        rank = rankof(pawn)
        bonus += PASS_PAWN_THREAT[rank]
        if rank > 3
            factor = muladd(5, rank, -13)
            block_sqr = pawn + 8
            bonus += makescore(0, (fld(min(DISTANCE_BETWEEN[block_sqr, black_king],5)*19, 4) - min(DISTANCE_BETWEEN[block_sqr, white_king],5)*2)*factor)
            if rank !== 7
                block_sqr_2 = block_sqr + 8
                bonus -= makescore(0, DISTANCE_BETWEEN[block_sqr_2, white_king]*factor)
            end
            if board[block_sqr] === BLANK
                sqrs_to_queen = PASSED_PAWN_MASKS[1][pawn] & FILE[file]
                unsafe_sqrs = PASSED_PAWN_MASKS[1][pawn]
                case = (PASSED_PAWN_MASKS[2][pawn] & FILE[file]) & rooklike(board)
                if isempty(case & black(board))
                    unsafe_sqrs &= b_attacks
                end

                if isempty(unsafe_sqrs)
                    k = 32
                elseif isempty(unsafe_sqrs & sqrs_to_queen)
                    k = 18
                elseif isempty(unsafe_sqrs & Bitboard(block_sqr))
                    k = 8
                else
                    k = 0
                end

                if !isempty(white(board) & case) || !isempty(w_attacks & Bitboard(block_sqr))
                    k += 5
                end

                bonus += makescore(k * factor, k * factor)
            end
        end
        if board[block_sqr] === WHITEPAWN
            bonus = makescore(fld(scoreMG(bonus), 2), fld(scoreEG(bonus), 2))
        end
        score += bonus
    end

    @inbounds for pawn in black_pawns
        bonus = 0
        file = fileof(pawn)
        rank = 9 - rankof(pawn)
        bonus += PASS_PAWN_THREAT[rank]
        if rank > 3
            factor = muladd(5, rank, -13)
            block_sqr = pawn - 8
            bonus += makescore(0, (fld(min(DISTANCE_BETWEEN[block_sqr, white_king],5)*14, 3) - min(DISTANCE_BETWEEN[block_sqr, black_king],5)*2)*factor)
            if rank !== 7
                block_sqr_2 = block_sqr - 8
                bonus -= makescore(0, DISTANCE_BETWEEN[block_sqr_2, black_king]*factor)
            end
            if board[block_sqr] === BLANK
                sqrs_to_queen = PASSED_PAWN_MASKS[2][pawn] & FILE[file]
                unsafe_sqrs = PASSED_PAWN_MASKS[2][pawn]
                case = (PASSED_PAWN_MASKS[1][pawn] & FILE[file]) & rooklike(board)
                if isempty(case & white(board))
                    unsafe_sqrs &= w_attacks
                end

                if isempty(unsafe_sqrs)
                    k = 32
                elseif isempty(unsafe_sqrs & sqrs_to_queen)
                    k = 18
                elseif isempty(unsafe_sqrs & Bitboard(block_sqr))
                    k = 8
                else
                    k = 0
                end

                if !isempty(black(board) & case) || !isempty(b_attacks & Bitboard(block_sqr))
                    k += 5
                end

                bonus += makescore(k * factor, k * factor)
            end
        end
        if board[block_sqr] === BLACKPAWN
            bonus = makescore(fld(scoreMG(bonus), 2), fld(scoreEG(bonus), 2))
        end
        score -= bonus
    end

    score
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
module SurrogatesRandomForest

import Surrogates: add_point!, AbstractSurrogate
export RandomForestSurrogate

using XGBoost
mutable struct RandomForestSurrogate{X,Y,B,L,U,N} <: AbstractSurrogate
    x::X
    y::Y
    bst::B
    lb::L
    ub::U
    num_round::N
end

function RandomForestSurrogate(x,y,lb::Number,ub::Number; num_round::Int = 1)
    bst = xgboost(reshape(x,length(x),1), num_round, label = y)
    RandomForestSurrogate(x,y,bst,lb,ub,num_round)
end

function (rndfor::RandomForestSurrogate)(val::Number)
    return XGBoost.predict(rndfor.bst,reshape([val],1,1))[1]
end

"""
RandomForestSurrogate(x,y,lb,ub,num_round)

Build Random forest surrogate. num_round is the number of trees.

"""
function RandomForestSurrogate(x,y,lb,ub;num_round::Int = 1)
    X = Array{Float64,2}(undef,length(x),length(x[1]))
    for j = 1:length(x)
        X[j,:] = vec(collect(x[j]))
    end
    bst = xgboost(X, num_round, label = y)
    RandomForestSurrogate(x,y,bst,lb,ub,num_round)
end

function (rndfor::RandomForestSurrogate)(val)
    return XGBoost.predict(rndfor.bst,reshape(collect(val),1,length(val)))[1]
end

function add_point!(rndfor::RandomForestSurrogate,x_new,y_new)
    if length(rndfor.lb) == 1
        #1D
        rndfor.x = vcat(rndfor.x,x_new)
        rndfor.y = vcat(rndfor.y,y_new)
        rndfor.bst = xgboost(reshape(rndfor.x,length(rndfor.x),1), rndfor.num_round, label = rndfor.y)
    else
        n_previous = length(rndfor.x)
        a = vcat(rndfor.x,x_new)
        n_after = length(a)
        dim_new = n_after - n_previous
        n = length(rndfor.x)
        d = length(rndfor.x[1])
        tot_dim = n + dim_new
        X = Array{Float64,2}(undef,tot_dim,d)
        for j = 1:n
            X[j,:] = vec(collect(rndfor.x[j]))
        end
        if dim_new == 1
            X[n+1,:] = vec(collect(x_new))
        else
            i = 1
            for j = n+1:tot_dim
                X[j,:] = vec(collect(x_new[i]))
                i = i + 1
            end
        end
        rndfor.x = vcat(rndfor.x,x_new)
        rndfor.y = vcat(rndfor.y,y_new)
        rndfor.bst = xgboost(X, rndfor.num_round, label = rndfor.y)
    end
    nothing
end
end # module

ext: jl


text:
using Documenter, XGrad

makedocs(
    format = :html,
    sitename = "XGrad.jl",
    pages = [
        "Main" => "index.md",
        "Tutorial" => "tutorial.md",
        "Code Discovery" => "codediscovery.md"
    ]
)

deploydocs(
    repo   = "github.com/dfdx/XGrad.jl.git",
    target = "build",
    julia = "0.6",
    deps   = nothing,
    make   = nothing
)


text:
"""
Overview: collections of functions required to solve the ε⁰-invariance equation
"""

"""
> solve_homological!(Cp,ndofs,p,sys_mat,sys_rhs,sys_res,ic,style,M,C,K)
It solves al homological equations for a given order of the asymtptotic developmet of the ε⁰-invariance.
- Cp : parametrisation structure
- ndofs : number of degrees of freedom
- p : order of the asymptotic developmet
- sys_mat : matrix used to solve the homological equations
- sys_rhs : right hand side of the homological equations
- sys_res : residual of the homological equations
- ic : index combination
- style : parametrisation style
- M : mass matrix
- C : damping matrix
- K : stiffness matrix
"""
function solve_homological!(Cp,ndofs,p,sys_mat,sys_rhs,sys_res,ic,
                            style,M,C,K)
  #
  resonant_modes = zeros(Bool,ndofs)
  #
  for i = 1:Cp[p+1].nc
    check = Cp[p+1].conj[i]
    if (check>0)
      σ = 0.0+0im
      for j = 1:p
        ic[j] = Cp[p+1].comb[j,i]
        σ += Cp[2].f[ic[j],ic[j]]
      end
      check_resonances!(Cp,σ,ndofs,style,resonant_modes)
      #
      assembly_sys_mat!(Cp,sys_mat.data,M.data,C.data,K.data,σ,resonant_modes,ndofs)
      assembly_sys_rhs!(Cp,sys_rhs,sys_res,M.data,C.data,i,p,ndofs,σ,resonant_modes)
      #
      sys_res = sparse(sys_mat)\sys_rhs
      #
      @inbounds for j = 1:M.data.m
        Cp[p+1].W[j,i] = Cp[p+1].W[j+M.data.m,i] + σ*sys_res[j]
        Cp[p+1].W[j+M.data.m,i] = sys_res[j]
      end
      #
      @inbounds for j = 1:ndofs
        if (resonant_modes[j]>0)
          Cp[p+1].f[j,i] = sys_res[M.data.m+j]
        end
      end
      #
      @inbounds for j = 1:ndofs
        for jj = 1:M.data.m
          Cp[p+1].W[jj,i] += Cp[p+1].f[j,i]*Cp[2].W[jj+M.data.m,j]
        end
      end
      #
    end
  end
  #
  nm = Int(ndofs/2)
  for i = 1:Cp[p+1].nc
    check = Cp[p+1].conj[i]
    tmp = abs(check)
    #
    if (check<0)
      #
      for j = 1:nm
        Cp[p+1].f[j+nm,i] = conj(Cp[p+1].f[j,tmp])
        Cp[p+1].f[j,i] = conj(Cp[p+1].f[j+nm,tmp])
      end
      #
      @inbounds for j = 1:M.data.m*2
        Cp[p+1].W[j,i] = conj(Cp[p+1].W[j,tmp])
      end
      #
    end
    #
  end
  #
  return nothing
  #
end

"""
> assembly_sys_rhs!(Cp,sys_rhs,sys_res,M,C,comb,p,ndofs,σ,resonant_modes)
It assembles the right hand side of a given homological equation of the ε⁰-developmet
- Cp : parametrisation structure
- sys_rhs : right hand side of the homological equations
- sys_res : residual of the homological equations
- M : mass matrix
- C : damping matrix
- comb : reference monomial index
- p : order of the asymptotic development
- ndofs : number of degrees of freedom
- σ : eingevalues λ summation
- resonant_modes : book-keeping of the set R of resonant modes
"""
function assembly_sys_rhs!(Cp,sys_rhs,sys_res,M,C,comb,p,ndofs,σ,resonant_modes)
  #
  fill!(sys_rhs,0.0)
  fill!(sys_res,0.0)
  nm = Int(ndofs/2)
  #
  @inbounds for irow = 1:M.m
    pos1 = M.colptr[irow]
    pos2 = M.colptr[irow+1]-1
    for j = pos1:pos2
      jcol = M.rowval[j]
      sys_rhs[irow] -= M.nzval[j]*Cp[p+1].W[jcol,comb] + (σ*M.nzval[j]+C.nzval[j])*Cp[p+1].W[jcol+M.m,comb]
      sys_res[irow] -= M.nzval[j]*Cp[p+1].W[jcol+M.m,comb]
      if (irow!=jcol)
        sys_rhs[jcol] -= M.nzval[j]*Cp[p+1].W[irow,comb] + (σ*M.nzval[j]+C.nzval[j])*Cp[p+1].W[irow+M.m,comb]
        sys_res[jcol] -= M.nzval[j]*Cp[p+1].W[irow+M.m,comb]
      end
    end
  end
  #
  for i = 1:ndofs
    if (resonant_modes[i]>0)
      if (i<=nm)
        @inbounds for j = 1:M.m
          sys_rhs[M.m+i] += Cp[2].W[M.m+j,i]*sys_res[j]
        end
      else
        @inbounds for j = 1:M.m
          sys_rhs[M.m+i] += Cp[2].W[M.m+j,i-nm]*sys_res[j]
        end
      end
    end
  end
  #
  @inbounds for i = 1:M.m
    sys_rhs[i] -= Cp[p+1].nlr[i,comb]
  end
  #
  return nothing
  #
end

"""
> assembly_sys_mat!(Cp,sys_mat,M,C,K,σ,resonant_modes,ndofs)
it assembles the linear system associated to a homological equation
- Cp : parametrisation structure
- sys_mat : matrix associated to the homological equation
- M : mass matrix
- C : damping matrix
- K : stiffness matrix
- σ : summation of eigenvalues λ
- resonant_modes : book-keeping of the set R of resonant modes
- ndofs : number of degrees of freedom
"""
function assembly_sys_mat!(Cp,sys_mat,M,C,K,σ,resonant_modes,ndofs)
  #
  fill!(sys_mat.nzval,0.0)
  nm = Int(ndofs/2)
  #
  for j = 1:K.m
    pos1 = K.colptr[j]
    pos2 = K.colptr[j+1]-1
    ne = pos2-pos1
    pos1_s = sys_mat.colptr[j]
    @inbounds for i = 0:ne
      sys_mat.nzval[pos1_s+i] = σ^2.0*M.nzval[pos1+i] + σ*C.nzval[pos1+i] + K.nzval[pos1+i]
    end
  end
  #
  for d = 1:ndofs
    if (!resonant_modes[d])
      pos = sys_mat.colptr[K.m+d]
      sys_mat.nzval[pos] = 1.0
    end
  end
  #
  for d = 1:nm
    if (resonant_modes[d])
      if (resonant_modes[d+nm])
        pos = sys_mat.colptr[K.m+d] + 1
        sys_mat.nzval[pos] = 1.0
      end
    end
  end
  #
  for d = 1:ndofs
    if (resonant_modes[d])
      if (d<=nm)
        λ = conj(Cp[2].f[d,d])
      else
        λ = Cp[2].f[d-nm,d-nm]
      end
      vec = (σ-λ)*(Symmetric(M,:L)*Cp[2].W[M.m+1:M.m*2,d])
      #
      @inbounds for j = 1:K.m
        pos = sys_mat.colptr[j+1]-1-ndofs+d
        sys_mat.nzval[pos] = vec[j]
      end
      #
      pos = sys_mat.colptr[K.m+d]
      sys_mat.nzval[pos] = 1.0
      #
    end
  end
  #
  return nothing
  #
end

"""
> check_resonances!(Cp,σ,ndofs,style,resonant_modes)
It checks which mode is resonant. It complies with the chosen style.
- Cp : parametrisation data structure
- σ : summation of eigenvalues λ
- ndofs : number of degrees of freedom
- style : parametrisation style
- resonant_modes : book-keeping of the set R of resonant modes
"""
function check_resonances!(Cp,σ,ndofs,style,resonant_modes)
  #
  fill!(resonant_modes,false)
  #
  if (style=='g')
    fill!(resonant_modes,true)
  else
    for i = 1:ndofs
      λ = Cp[2].f[i,i]
      if (abs(σ-λ)/abs(λ)<=ϵ_tol)
        resonant_modes[i] = true
        if (style=='r')
          if (i <= Int(ndofs/2))
            resonant_modes[i+Int(ndofs/2)] = true
          else
            resonant_modes[i-Int(ndofs/2)] = true
          end
        end
      end
    end
  end
  #
  return nothing
end

"""
> recursive_assembly!(Cp,ndofs,p,mesh,U,ic,cc,nls=0)
It assembles all right hand sides of the homological equations for a given order p of the expansion of the ε⁰ development
- Cp : parametrisation data structure
- ndofs : number of degrees of freedom
- p : order of the asymptotic development
- mesh : mesh data structure
- U : dispalcement field
- ic : index combination of a given monomial
- cc : counter of the asymptotic development
- nls : tag for the nonlinear static analysis. 0 = yes, 1 = yes
"""
function recursive_assembly!(Cp,ndofs,p,mesh,U,ic,cc,nls=0)
  #
  for i = 1:ndofs
    ic[cc] = i
    if (cc<p)
      cc += 1
      recursive_assembly!(Cp,ndofs,p,mesh,U,ic,cc,nls)
      cc -= 1
    else
      #
      pos = find_index_global(ic,ndofs,p)

      entry = Cp[p+1].cmap[pos]
      check = Cp[p+1].conj[abs(entry)]
      #
      if (check>0)
        assembly_quadratic_nl!(Cp,entry,ic,ndofs,p,mesh,U,nls)
        assembly_cubic_nl!(Cp,entry,ic,ndofs,p,mesh,U)
        assembly_μ_ν_resid!(Cp,entry,ic,ndofs,p,U.neq)
      end
      #
    end
  end
  #
  return nothing
  #
end

"""
> assembly_μ_ν_resid!(Cp,entry,ic,ndofs,p,neq)
It assembles μ and ν for a given monomial I of order p for the ε⁰ development.
- Cp : parametrisation data structure
- entry : entry associated to a monomial combination
- ic : monomial index combination
- ndofs : number of degrees of freedom
- p : order of the asymptotic development
- neq : number of equations
"""
function assembly_μ_ν_resid!(Cp,entry,ic,ndofs,p,neq)
  #
  ent = abs(entry)
  #
  for j = 2:p-1
    for k = 0:p-j
      posf = find_index_rdyn(j,k,ndofs,ic)
      cmf = Cp[j+1].cmap[posf]
      if (cmf>0)
        f = @view Cp[j+1].f[:,cmf]
        for s = 1:ndofs
          posw = find_index_map(j,k,s,p,ndofs,ic)
          cmw = Cp[p-j+2].cmap[posw]
          if (cmw>0)
            W = @view Cp[p-j+2].W[:,cmw]
            @inbounds for i = 1:neq*2
              Cp[p+1].W[i,ent] += f[s]*W[i]
            end
          end
        end
      end
    end
  end
  #
  return nothing
  #
end

"""
> assembly_cubic_nl!(Cp,entry,ic,ndofs,p,mesh,U)
It assembles cubic nonlinearities vector for a given monomial I of order p
- Cp : parametrisation data structure
- entry : index entry
- ic : index combination of the monomial
- ndofs : number or degrees of freedom
- p : order of the asymptotic development
- mesh : grid data structure
- U : displacement field
"""
function assembly_cubic_nl!(Cp,entry,ic,ndofs,p,mesh,U)
  #
  for j = 1:p-1
    for k = 1:p-j-1
      #
      pos1,pos2,pos3 = find_position_h(j,k,p,ndofs,ic)
      #
      cm1 = Cp[j+1].cmap[pos1]
      cm2 = Cp[k+1].cmap[pos2]
      cm3 = Cp[p-j-k+1].cmap[pos3]
      #
      if (cm1>0 && cm2>0 && cm3>0)
        #
        Ψ₁ = @view Cp[j+1].W[U.neq+1:2*U.neq,cm1]
        Ψ₂ = @view Cp[k+1].W[U.neq+1:2*U.neq,cm2]
        Ψ₃ = @view Cp[p-j-k+1].W[U.neq+1:2*U.neq,cm3]
        #
        assembly_H_nl!(Cp[p+1], abs(entry), Ψ₁, Ψ₂, Ψ₃, mesh, U)
        #
      end
      #
    end
  end
  #
  return nothing
  #
end

"""
> assembly_quadratic_nl!(Cp,entry,ic,ndofs,p,mesh,U,nls)
It assembles quadratic nonlinearities vector for a given monomial I of order p
- Cp : parametrisation data structure
- entry : index entry
- ic : index combination of the monomial
- ndofs : number or degrees of freedom
- p : order of the asymptotic development
- mesh : grid data structure
- U : displacement field
- nls : tag for nonlinear static analysis. 0 = no, 1 = yes
"""
function assembly_quadratic_nl!(Cp,entry,ic,ndofs,p,mesh,U,nls)
  #
  for j = 1:p-1
    pos1, pos2 = find_position_g(j,p,ndofs,ic)
    #
    cm1 = Cp[j+1].cmap[pos1]
    cm2 = Cp[p-j+1].cmap[pos2]
    #
    if (cm1>0 && cm2>0)
      #
      Ψ₁ = @view Cp[j+1].W[U.neq+1:2*U.neq,cm1]
      Ψ₂ = @view Cp[p-j+1].W[U.neq+1:2*U.neq,cm2]
      #
      assembly_G_nl!(Cp[p+1], abs(entry), Ψ₁, Ψ₂, mesh, U)
      if (nls>0)
        assembly_H_nl!(Cp[p+1], abs(entry), Ψ₁, Ψ₂, Cp[1].W, mesh, U, 3.0)
      end
    end
  end
  #
  return nothing
  #
end


text:
@testset "ModAlgAss" begin
  include("ModAlgAss/Lattices.jl")
end


text:
### A Pluto.jl notebook ###
# v0.18.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ ac99559c-880b-11ec-164b-ef9d96149fb9
begin
	import Pkg
	Pkg.activate("Project.toml")
	Pkg.instantiate()
end

# ╔═╡ e5936aaa-59ec-48bf-91af-46f7bb3987c3
using Plots

# ╔═╡ 492be7c0-c211-4c40-8638-508183ddc615
using FinSight

# ╔═╡ 21f7cc65-851a-42a9-9b6d-ad8f5e188f66
using PlutoUI

# ╔═╡ bfb73c6d-b36d-4d9b-b549-8ce8775b1be9
using Measurements

# ╔═╡ ad43845d-81fe-41ae-b3f9-1eebded9e399

values = [-10000, 3000, 4200, 5800]

# ╔═╡ 38bf894e-fcae-4d20-98c6-6de7e73ebcfc
uncertanties = [1000, 100, 300, 400]

# ╔═╡ f8a0b02e-e913-4b97-92c7-dbac971ed569
values .± uncertanties

# ╔═╡ c6275a35-8093-432c-b5ed-db36eb8e7189
@bind rate Slider(range(0.1, stop=0.9, length=10))

# ╔═╡ dfd301bc-8420-41e7-b2d8-4e570dfdaad8
npv(rate, values .± uncertanties)

# ╔═╡ 594461e0-4da5-48c8-98cc-31232a7770ba
rate

# ╔═╡ 57f66f1e-3cf2-4db8-8ae6-3cb8f5069e61
npv(rate, values)

# ╔═╡ d0f16333-dfb6-437b-9007-9b07da70b64a
rates = range(0.1, stop=0.9, length=10)

# ╔═╡ 5702265d-e0dc-45a6-9e72-8df3c708c718
plot(rates,
	map(x->npv(x, values), rates),

	m=:c, xlab="rate", ylab="npv")

# ╔═╡ 3640b925-9670-4b37-a357-64cd91f6ef98
rates

# ╔═╡ 73fd2d92-1d35-451e-be2e-b9975f7296fb
[npv(i, values) for i in rates]

# ╔═╡ Cell order:
# ╠═ac99559c-880b-11ec-164b-ef9d96149fb9
# ╠═e5936aaa-59ec-48bf-91af-46f7bb3987c3
# ╠═492be7c0-c211-4c40-8638-508183ddc615
# ╠═21f7cc65-851a-42a9-9b6d-ad8f5e188f66
# ╠═bfb73c6d-b36d-4d9b-b549-8ce8775b1be9
# ╠═ad43845d-81fe-41ae-b3f9-1eebded9e399
# ╠═38bf894e-fcae-4d20-98c6-6de7e73ebcfc
# ╠═f8a0b02e-e913-4b97-92c7-dbac971ed569
# ╠═dfd301bc-8420-41e7-b2d8-4e570dfdaad8
# ╠═c6275a35-8093-432c-b5ed-db36eb8e7189
# ╠═594461e0-4da5-48c8-98cc-31232a7770ba
# ╠═57f66f1e-3cf2-4db8-8ae6-3cb8f5069e61
# ╠═d0f16333-dfb6-437b-9007-9b07da70b64a
# ╠═5702265d-e0dc-45a6-9e72-8df3c708c718
# ╠═3640b925-9670-4b37-a357-64cd91f6ef98
# ╠═73fd2d92-1d35-451e-be2e-b9975f7296fb

max_stars_repo_name: davidtaylor-fastlanebi/FinSight

max_issues_repo_name: davidtaylor-fastlanebi/FinSight

max_forks_repo_name: davidtaylor-fastlanebi/FinSight


text:
function pluralize(word::String)
    if (endswith(word,"eau")
     || endswith(word,"eu")
     || endswith(word,"au"))
       return "$(word)x"
    elseif (endswith(word,"ou"))
        return "$(word)s"
    elseif (endswith(word,"al"))
        return replace(word,r"al$" => "aux")
    elseif (endswith(word,"ail"))
        return "$(word)s"
    elseif (endswith(word,r"[abcdefghijklmnopqrtuvwy]"))
        return "$(word)s"
    else
        return word
    end
end


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 6c692a11-cb4c-425d-9acb-9529bd382dd2
begin
    import Pkg
	Pkg.activate(Base.current_project())

	using FringeExplorer, Latexify, RadioInterferometry, Symbolics, Plots, PlotlyJS, LaTeXStrings, PlutoUI
	Latexify.set_default(fmt=FancyNumberFormatter(4, "{\\cdot}"))
	#Latexify.set_default(fmt=FancyNumberFormatter("%.1g",
	#	SubstitutionString{String}("\\g<mantissa> {\\\\cdot} 10^{\\g<exp>}")))
	#LaTeXString(numformat(latexify(FringeExplorer.wdw_m)))
	md"""
	# Fringe Explorer
	This notebook is for exploring the relative geometric delays, and their rates of change, as seen by antennas in a terrestrial radio telescope array as it tracks a sidereal source moving across the sky.

	Of course we now know that the ground based antennas are doing (most of) the moving thanks to Earth's rotation, but as observers who live in the same reference frame as the antennas it can sometimes be convenient to think in terms of the sky moving instead.
	"""
end

# ╔═╡ 655091ad-69a0-4672-abb4-37ee03b27cff
md"""
## Getting Oriented
### XYZ frame
For the purposes of this discussion we will represent antenna positions, and the array's reference position, in a right-handed XYZ frame where the Z axis is parallel to the rotational axis of Earth and the XY plane parallel to the equarorial plane of Earth.  One such frame commonly used for antenna positions is the *International Terrestrial Reference Frame* (ITRF), a geocentric frame where Z is the rotational axis of Earth, X is in the equatorial lane and passes through the prime meridian (i.e. longitude 0), and Y completes the right hand frame (i.e. passing though longitude 90 degrees).  The frame used here is a modified version of this frame where the origin has been translated to the array reference position and the frame has been rotated around the Z axis such that X points to the longitude of the array reference position.
"""

# ╔═╡ 7ff31ce8-190e-4018-ba49-c2f11ed82801
md"""
### Hour Angle and Declination frame
We will also choose to represent sidereal positions as an *hour angle*, represented by the variable ``h``, and *declination* angle, represented by the variable ``\delta``.  Being angles, we will store and calculate with them in radians, but often display them as degrees.  This choice of representation is convenient because for a given sidereal position, ``h`` is a linear function of time and ``\delta`` is a constant.  Positions on the local meridian have ``h = 0``, positions to the west have ``h > 0``.  Using ``\Omega`` to represent rate of rotation, ``h`` for a given sideral source can be expressed as ``h = \Omega t``, with ``t = 0`` corresponding to the time when the source is *transiting*, i.e. on/crossing the local meridian.
"""

# ╔═╡ 423f2ef7-932d-4eb8-8ea6-304e2b7cbf2e
md"""
### UVW frame
With these conventions established, we can now discuss how to compute the relative geometric delays between wavefront arrival times at a given antenna and the array reference position.  One way to do this is to rotate the XYZ frame we are using for antenna positions into a new frame known as the UVW frame (also right handed).  The U axis of the UVW frame points eastward, the V axis points northward, and the W axis points at the source.  Because the W axis points to the source that is moving, the UVW frame is not a fixed frame, but a rotating frame.  When representing the antenna positions in the UVW frame, the W component of this position vector represents to *line of sight* distance between the antenna and the reference position as seen from the source position.  In other words, the RF wavefront arrives at one end of the vector first, then travels a distance W before arriving at the other end of the vector.  For our purposes here we will ignore the amount of Earth rotation that happens during this time.  The XYZ coordinates of the antennas are typically specified in meters, but speed of light can be used to convert W from meters to seconds, or more typically nanoseconds, of delay which can be used along with observing frequency to convert W to units of wavelength.
"""

# ╔═╡ c02e8687-be77-4af0-a0a1-6901ad02f85e
md"""
Rotating antenna position vectors from the XYZ frame to the UVW frame can be accomplished by rotating our XYZ frame about the Z (i.e. third) axis by ``h``, producing a (X',U,Z) frame, then rotating that frame about the U (i.e. second) axis by ``\delta``, producing a (W,U,V) frame which can then permuted to (U,V,W) where U is east, V is north, and W is in the direction of ``h`` and ``\delta``.  This operation can be performed by left multiplying the antenna position vector by a 3x3 rotation matrix.

The `RadioInterferometry.jl` package provides a `xyz2uvw` function that can generate these 3x3 rotation matrices for any given ``h`` and ``\delta``.

!!! tip
    The `xyz2uvw` function supports XYZ frames that, unlike ours, are oriented with X parallel to the great circle containing the prime meridian (e.g. the ITRF XYZ frame).  This involves adjusting ``h`` by the logitude of the array's reference position, but since our X axis is already oriented to the local meridian we can simply pass 0 for longitude.

The rotation from XYZ to UVW for a given hour angle ``h`` and declination ``\delta`` can be expressed as:
"""

# ╔═╡ 0c19d2b3-3869-48f1-b2eb-1ff157e931cf
begin
	@variables x y z u v w h δ
	([u,v,w] ~ :($(xyz2uvw(h, δ, 0)) * $([x,y,z])))|>latexify
end

# ╔═╡ 11567a37-55b2-4d9e-88fb-9efdf7ada485
md"""
## Play time!

Change the sliders here to select between hour angle (in minutes) and declination (in degrees).  The values you selected are shown below followed by the XYZ to UVW transformation for the given parameters.

$(@bind h0 Slider(-360:5:360, default=0)) Hour angle (minutes)

$(@bind δ0 Slider(-90:5:90, default=0)) Declination (degrees)
"""

# ╔═╡ 1ceb4ee8-2bbd-45a7-b5a4-886c949c17e3
(h=h0, δ=δ0)

# ╔═╡ aee119f3-b190-4342-9d48-c0ee9f8e2859
begin
	([u,v,w] ~ :($(xyz2uvw(deg2rad(h0/4), deg2rad(δ0), 0)) * $([x,y,z])))|>latexify
end

# ╔═╡ 979c1962-db9a-4581-95a8-15b183afd990

# ╔═╡ d4145813-5e18-4230-8389-d8935d72ba4b
FringeExplorer.wdw_m

# ╔═╡ Cell order:
# ╟─6c692a11-cb4c-425d-9acb-9529bd382dd2
# ╟─655091ad-69a0-4672-abb4-37ee03b27cff
# ╟─7ff31ce8-190e-4018-ba49-c2f11ed82801
# ╟─423f2ef7-932d-4eb8-8ea6-304e2b7cbf2e
# ╟─c02e8687-be77-4af0-a0a1-6901ad02f85e
# ╟─0c19d2b3-3869-48f1-b2eb-1ff157e931cf
# ╟─11567a37-55b2-4d9e-88fb-9efdf7ada485
# ╟─1ceb4ee8-2bbd-45a7-b5a4-886c949c17e3
# ╟─aee119f3-b190-4342-9d48-c0ee9f8e2859
# ╠═979c1962-db9a-4581-95a8-15b183afd990
# ╠═d4145813-5e18-4230-8389-d8935d72ba4b


text:
__precompile__()

module Common

using UnicodePlots
using JLD

function addmachine(name, count, tunnel::Bool=false)
   #topo=:all_to_all
   topo=:master_slave
   println("ADDING $name to cluster")
   if tunnel
      addprocs([(name, count)]; max_parallel=1, topology=topo, tunnel=true)
   else
      addprocs([(name, count)]; max_parallel=1, topology=topo)
   end
   println("WORKERS: ", nworkers())
end

# ornstein-uhlenbeck process for ctrl noise for stochastic search
# μ  = target value of x[end]; get some random target?
# σ  = noise. sure.
# θ  = how quickly to go to μ; fights against σ
# dt = dt of model * skip; dt of environment
# x0 = initial value of x[1]
function ou_process!(x::AbstractMatrix{Float64},
                     σ::Float64=1.0, θ::Float64=2.0,
                     dt::Float64=1/size(x,2))
   nu = size(x,1)
   T  = size(x,2)-1
   W  = zeros(nu)
   x0 = copy(x[:,1])
   μ  = copy(x[:,end])
   for (i,t) in enumerate(0.0:dt:T*dt)
      ex = exp(-θ*t)
      for j=1:nu
         x[j,i] = x0[j]*ex + μ[j]*(1.0-ex) + σ*ex*W[j]
         W[j] += exp(θ*t)*sqrt(dt)*randn()
      end
   end
end

function ou_process!(x::AbstractVector{Float64}, x0::Float64,
                     μ::Float64=(2rand()-1),
                     σ::Float64=1.0, θ::Float64=2.0,
                     dt::Float64=1/length(x))
   T = length(x)-1
   W = 0.0
   for (i,t) in enumerate(0.0:dt:T*dt)
      ex = exp(-θ*t)
      x[i] = x0*ex + μ*(1.0-ex) + σ*ex*W
      W += exp(θ*t)*sqrt(dt)*randn()
   end
end

function plotlines(iter::Integer,title::String,lines::Matrix)
   yin = minimum(lines)
   yax = maximum(lines)
   p = lineplot(lines[1,:],
                xlim=[0, size(lines,2)],
                ylim=[yin, yax],
                title=title,
                width=60, height=10)
   for l=2:size(lines, 1)
      lineplot!(p, lines[l,:])
   end

   display(p)
end

function plotlines(iter::Integer,title::String,lines::Vector)
   yin = minimum(lines)
   yax = maximum(lines)
   p = lineplot(lines,
                xlim=[0, size(lines,2)],
                ylim=[yin, yax],
                title=title,
                width=60, height=10)

   display(p)
end

function plotlines(iter::Integer,title::String,lines...)
   if length(lines) < 1
      error("plotlines needs data and title")
   end

   p = lineplot(lines[1][1][:],
                title=title,
                width=60, height=10, name=lines[1][2])
   for l in lines[2:end]
      lineplot!(p, l[1][:], name=l[2])
   end

   display(p)
end

function plotexpmt(file)
   d = load(file)
   plotlines(length(d["stocR"]),"Reward",
             (d["stocR"],"Stoc"),
             (d["meanR"],"Mean"))
   return d
end

function makesavedir(dir_name::String, exp_file::String,
                     m_file::String, mesh_file::String="",
                     overwrite::Bool=false)
   if isdir(dir_name) == false
      mkdir(dir_name)
   else
      if overwrite == false
         val = 1
         while isdir("$(dir_name)_$(val)") == true
            val += 1
         end
         dir_name = "$(dir_name)_$(val)"
         mkdir(dir_name)
      end
   end
   # copy file to save directory
   filename = basename(exp_file)
   cp(exp_file, "$(dir_name)/$(filename)", remove_destination=overwrite)
   info("Saving experiment $exp_file")
   open(exp_file) do f
      for l in eachline(f)
         if startswith(l, "include")
            includefile = basename(split(l, "\"")[2])
            info("Including $includefile")
            cp(dirname(exp_file)*"/"*includefile, "$(dir_name)/$(includefile)", remove_destination=overwrite)
         end
      end
   end
   filename = basename(m_file)
   cp(m_file, "$(dir_name)/$(filename)", remove_destination=overwrite)
   info("Model file $m_file")
   if (mesh_file != "")
      if isdir(mesh_file)
         symlink(mesh_file, dir_name*"/meshes")
      else
         symlink(dirname(mesh_file), dir_name*"/meshes")
      end
   end
   return dir_name
end

function save(expresults::String, data::Dict{Symbol, Vector{Float64}})
    jldopen(expresults, "w") do file
        for d in data
           write(file, "$(d[1])", d[2])
        end
    end
end

end


text:
# # Wiring diagrams in Graphviz
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/graphics/graphviz_wiring_diagrams.ipynb)
#
# Catlab can draw wiring diagrams using the `dot` program in
# [Graphviz](https://www.graphviz.org/). This feature requires that Graphviz be
# installed, but does not require any additional Julia packages.

using Catlab.WiringDiagrams, Catlab.Graphics

# ## Examples

# ### Symmetric monoidal category

using Catlab.Doctrines

A, B = Ob(FreeSymmetricMonoidalCategory, :A, :B)
f = Hom(:f, A, B)
g = Hom(:g, B, A)
h = Hom(:h, otimes(A,B), otimes(A,B));

# To start, here are a few very simple examples.

to_graphviz(f)
#-
to_graphviz(compose(f,g))
#-
to_graphviz(otimes(f,g))

# In the next example, notice how Graphviz automatically "untwists" the double
# braiding to minimize edge crossings.

to_graphviz(compose(braid(A,A), otimes(f,f), braid(B,B)))

# Here is a larger composite morphism.

composite = compose(otimes(g,f), h, otimes(f,g))
to_graphviz(composite)

# By default, the wiring diagram is laid out from top to bottom. Other layout
# orientations can be requested, such as left-to-right or bottom-to-top:

to_graphviz(composite, orientation=LeftToRight)
#-
to_graphviz(composite, orientation=BottomToTop)

# When working with very large diagrams (larger than the ones shown here), it is
# sometimes convenient to omit the ports of the outer box and any wires attached
# to them.

to_graphviz(composite, outer_ports=false)

# ### Biproduct category

A, B = Ob(FreeBiproductCategory, :A, :B)
f = Hom(:f, A, B)
g = Hom(:g, B, A);

# By default, copies and merges are drawn the way they are represented
# internally, as multiple wires.

f1 = compose(mcopy(A), otimes(f,f))
to_graphviz(f1)
#-
f2 = compose(mcopy(A), otimes(f,f), mmerge(B))
to_graphviz(f2)

# To draw nodes for copies and merges, we need to add junctions to the wiring
# diagram.

to_graphviz(add_junctions!(to_wiring_diagram(f1)))
#-
to_graphviz(add_junctions!(to_wiring_diagram(f2)))

# ### Traced monoidal category

A, B, X, Y = Ob(FreeTracedMonoidalCategory, :A, :B, :X, :Y)
f = Hom(:f, otimes(X,A), otimes(X,B))

to_graphviz(trace(X, A, B, f))
#-
to_graphviz(trace(X, A, B, f), orientation=LeftToRight)
#-
g, h = Hom(:g, A, A), Hom(:h, B, B)

trace_naturality = trace(X, A, B, compose(otimes(id(X),g), f, otimes(id(X),h)))
to_graphviz(trace_naturality, orientation=LeftToRight)

# ## Custom styles

# The visual appearance of wiring diagrams can be customized by setting Graphviz
# [attributes](https://www.graphviz.org/doc/info/attrs.html) at the graph, node,
# edge, and cell levels. Graph, node, and edge attributes are described in the
# Graphviz documentation. Cell attributes are passed to the primary cell of the
# [HTML-like label](https://www.graphviz.org/doc/info/shapes.html#html) used for
# the boxes.

A, B, C = Ob(FreeSymmetricMonoidalCategory, :A, :B, :C)
f, g = Hom(:f, A, B), Hom(:g, B, C)

to_graphviz(compose(f,g),
  labels = true, label_attr=:headlabel,
  node_attrs = Dict(
    :fontname => "Courier",
  ),
  edge_attrs = Dict(
    :fontname => "Courier",
    :labelangle => "25",
    :labeldistance => "2",
  ),
  cell_attrs = Dict(
    :bgcolor => "lavender",
  )
)

# ## Output formats

# The function `to_graphviz` returns an object of a type `Graphviz.Graph`,
# representing a Graphviz graph as an abstract syntax tree. When displayed
# interactively, this object is automatically run through Graphviz and rendered
# as an SVG image. Sometimes it is convenient to perform this process manually,
# to change the output format or further customize the generated dot file.

# To generate a dot file, use the builtin pretty-printer. This feature does not
# require Graphviz to be installed.

using Catlab.Graphics: Graphviz

graph = to_graphviz(compose(f,g))
Graphviz.pprint(graph)

# Catlab provides a simple wrapper around the Graphviz command-line programs.
# For example, here the JSON output for the graph.

import JSON

JSON.parse(Graphviz.run_graphviz(graph, format="json0"))


text:
using CBindingGen
using Libdl
using MPI
import Pkg.TOML
using Pkg.Artifacts
import P4est_jll

# setup configuration using ideas from MPI.jl
const config_toml = joinpath(first(DEPOT_PATH), "prefs", "P4est.toml")
mkpath(dirname(config_toml))

if !isfile(config_toml)
  touch(config_toml)
end

config = TOML.parsefile(config_toml)

# P4est.toml has the following keys:
#  p4est_generate_bindings = "" (default) | if non-empty, re-generate bindings during build stage
#  p4est_path              = "" (default) | path to p4est containing subdirectories lib and include
#  p4est_library           = "" (default) | library name/path
#  p4est_include           = "" (default) | include name/path
#  p4est_uses_mpi          = "" (default) | "yes" indicates that we need the MPI headers
#  mpi_path                = "" (default) | path to MPI containing subdirectories lib and include
#  mpi_include             = "" (default) | include name/path

# Step 1: Check environment variables and update preferences accordingly
if haskey(ENV, "JULIA_P4EST_GENERATE_BINDINGS")
  config["p4est_generate_bindings"] = ENV["JULIA_P4EST_GENERATE_BINDINGS"]
else
  config["p4est_generate_bindings"] = ""
end

if haskey(ENV, "JULIA_P4EST_PATH")
  config["p4est_path"] = ENV["JULIA_P4EST_PATH"]
else
  config["p4est_path"] = ""
end

if haskey(ENV, "JULIA_P4EST_LIBRARY")
  config["p4est_library"] = ENV["JULIA_P4EST_LIBRARY"]
else
  config["p4est_library"] = ""
end

if haskey(ENV, "JULIA_P4EST_INCLUDE")
  config["p4est_include"] = ENV["JULIA_P4EST_INCLUDE"]
else
  config["p4est_include"] = ""
end

if haskey(ENV, "JULIA_P4EST_USES_MPI")
  config["p4est_uses_mpi"] = ENV["JULIA_P4EST_USES_MPI"]
else
  config["p4est_uses_mpi"] = ""
end

if haskey(ENV, "JULIA_P4EST_MPI_PATH")
  config["mpi_path"] = ENV["JULIA_P4EST_MPI_PATH"]
else
  config["mpi_path"] = ""
end

if haskey(ENV, "JULIA_P4EST_MPI_INCLUDE")
  config["mpi_include"] = ENV["JULIA_P4EST_MPI_INCLUDE"]
else
  config["mpi_include"] = ""
end

open(config_toml, "w") do io
  TOML.print(io, config)
end

if isempty(config["p4est_generate_bindings"])
  println("Use pre-generated bindings for p4est")
  const bindings_filename = joinpath(@__DIR__, "libp4est.jl")
  const pre_generated_bindings_filename = joinpath(artifact"libp4est", "libp4est.jl")
  cp(pre_generated_bindings_filename, bindings_filename, force=true)
else
  # Step 2: Choose p4est library according to the settings
  p4est_library = ""
  if !isempty(config["p4est_library"])
    p4est_library = config["p4est_library"]
    println("Use custom p4est library $p4est_library")
  elseif !isempty(config["p4est_path"])
    p4est_library = joinpath(config["p4est_path"], "lib", "libp4est." * Libdl.dlext)
    if isfile(p4est_library)
      println("Use custom p4est library $p4est_library")
    else
      p4est_library = ""
    end
  end

  if isempty(p4est_library)
    p4est_library = P4est_jll.libp4est_path
    println("Use p4est library provided by P4est_jll")
  end

  # Step 3a: Choose the p4est include path according to the settings
  include_directories = String[]
  p4est_include = ""
  if !isempty(config["p4est_include"])
    p4est_include = config["p4est_include"]
    println("Use custom p4est include path $p4est_include")
  elseif !isempty(config["p4est_path"])
    p4est_include = joinpath(config["p4est_path"], "include")
    if isdir(p4est_include)
      println("Use custom p4est include path $p4est_include")
    else
      p4est_include = ""
    end
  end

  if isempty(p4est_include)
    p4est_include = joinpath(dirname(dirname(P4est_jll.libp4est_path)), "include")
    println("Use p4est include path provided by P4est_jll")
  end

  push!(include_directories, p4est_include)

  # Step 3b: Choose the MPI include path according to the settings
  if config["p4est_uses_mpi"] == "yes"
    mpi_include = ""
    if !isempty(config["mpi_include"])
      mpi_include = config["mpi_include"]
      println("Use custom MPI include path $mpi_include")
    elseif !isempty(config["mpi_path"])
      mpi_include = joinpath(config["mpi_path"], "include")
      if isdir(mpi_include)
        println("Use custom MPI include path $mpi_include")
      else
        mpi_include = ""
      end
    end

    if isempty(mpi_include)
      mpi_include = joinpath(dirname(dirname(normpath(dlpath(MPI.libmpi)))), "include")
      println("Use MPI include path based on path to `MPI.jl`'s `libmpi`")
    end

    if !isfile(joinpath(mpi_include, "mpi.h"))
      error("could not find `mpi.h` in ", mpi_include)
    end

    push!(include_directories, mpi_include)
  end

  # Step 4: Generate binding using the include path according to the settings

  # Manually set header files to consider
  hdrs = ["p4est.h", "p4est_extended.h",
          "p6est.h", "p6est_extended.h",
          "p8est.h", "p8est_extended.h"]

  # Build list of arguments for Clang
  include_args = String[]
  @show include_directories
  for dir in include_directories
    append!(include_args, ("-I", dir))
  end

  # Workaround for MacOS: The some headers required by p4est (such as `math.h`) are only available via
  # Xcode
  if Sys.isapple()
    # These two paths *should* - on any reasonably current MacOS system - contain the relevant headers
    const xcode_include_path_cli = "/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/"
    const xcode_include_path_gui = "/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/"
    if !isdir(xcode_include_path_cli) && !isdir(xcode_include_path_gui)
      error("MacOS SDK include paths ('$xcode_include_path_cli' or '$xcode_include_path_gui') do not exist. Have you installed Xcode?")
    end
    append!(include_args, ("-idirafter", xcode_include_path_cli, "-idirafter", xcode_include_path_gui))
  end

  # Convert symbols in header
  cvts = convert_headers(hdrs, args=include_args) do cursor
    header = CodeLocation(cursor).file
    name   = string(cursor)

    # only wrap the libp4est and libsc headers
    dirname, filename = splitdir(header)
    if !(filename in hdrs ||
        startswith(filename, "p4est_") ||
        startswith(filename, "p6est_") ||
        startswith(filename, "p8est_") ||
        startswith(filename, "sc_") ||
        filename == "sc.h" )
      return false
    end

    # Ignore macro hacks
    startswith(name, "sc_extern_c_hack_") && return false

    return true
  end

  # Write generated C bindings to file
  const bindings_filename = joinpath(@__DIR__, "libp4est.jl")
  open(bindings_filename, "w+") do io
    generate(io, p4est_library => cvts)
  end
end


text:
function update_status(client, config, bot_path)
    files = glob(joinpath(bot_path[2:end], "history/*/*.json"))
    profit = 0
    trade_num = 0
    sold = 0
    for file in files
        tmp = JSON.parsefile(file)
        for (day,trades) in tmp
            for trade in trades
                profit += trade["profit"]
                trade_num += 1
                if trade["type"] == "initial"
                    sold += typeof(trade["order"]["executedQty"]) == Float64 ? trade["order"]["executedQty"] : parse(Float64, trade["order"]["executedQty"])
                end
            end
        end
    end
    if !haskey(config["init_state"], "sold")
        config["init_state"]["sold"] = 0
    end
    if config["init_state"]["sold"] != sold
        config["init_state"]["sold"] = sold
        open(joinpath(bot_path, "config.json"), "w") do f
            JSON.print(f, config, 4)
        end
    end
    trade_grid = JSON.parsefile(joinpath(bot_path, "trade_grid.json"))
    current_price = findfirst(x->x["status"] == "current", trade_grid)
    current_gridpoint = findfirst(x->x == current_price, sort(string.(keys(trade_grid))))

    next_gridpoints = Pinguu.get_next_gridpoints(trade_grid, current_price)

    ticker, current_ticker = Pinguu.get_current_ticker(client, config, use_api=true)
    next_sell = findall(x -> x["status"] == "check" && x["side"] == client.SIDE_SELL, trade_grid)
    if length(next_sell) > 0
        next_sell = next_sell[1]
        to_sell = round(100*(1 - current_ticker / parse(Float64, next_sell)), digits=2)
    else
        to_sell = 0
    end

    #
    ## APR
    files = glob(joinpath(bot_path[2:end], "history/*/*.json"));
    total_profit_percent = 0
    profit_thes = 0
    for file in files
        tmp = JSON.parsefile(file)
        for (day, trades) in tmp
            for trade in trades
                total_profit_percent += trade["profit_percent"]
                profit_thes += trade["profit_thes"]
            end
        end
    end
    timespan = now() - DateTime(config["start"])
    rendite = total_profit_percent / (timespan.value/(1000*3600*24*365))
    trade_num_h = trade_num / (timespan.value/(1000*3600))
    trade_num_d = trade_num / (timespan.value/(1000*3600*24))
    if config["thes_volume"] != round(profit_thes, digits=config["sell_precision"])
        config["thes_volume"] = round(profit_thes, digits=config["sell_precision"])
        open(joinpath(bot_path, "config.json"), "w") do f
            JSON.print(f, config, 4)
        end
    end

    #
    ## Currently invested in coin
    bought = config["init_state"]["cummulativeQuoteQty"] * (1-config["init_state"]["sold"] / config["init_state"]["executedQty"])
    for (k,v) in trade_grid
        if haskey(v, "buy_order")
            bought += typeof(v["buy_order"]["cummulativeQuoteQty"]) == Float64 ? v["buy_order"]["cummulativeQuoteQty"] : parse(Float64, v["buy_order"]["cummulativeQuoteQty"])
        end
    end

    #
    ## Daily average profit
    current_value = Pinguu.get_current_volume(config)
    daily_avg = (current_value * rendite/100)/365

    #
    ## Today profit
    today_profit = 0
    yr, m, d = split(string(today()), "-")
    profit_file = joinpath(bot_path[2:end], "history/" * yr * "/" * m * ".json")
    if isfile(profit_file)
        tmp = JSON.parsefile(profit_file)
        if haskey(tmp,string(today()))
            for trade in tmp[string(today())]
                today_profit += trade["profit"]
            end
        end
    end

    #
    ## Status dictionary
    if haskey(next_gridpoints, "BUY")
        progress_line = round((current_ticker - parse(Float64, next_gridpoints["BUY"]))/((parse(Float64, next_gridpoints["SELL"]) - parse(Float64, next_gridpoints["BUY"]))/100))
    else
        progress_line = 0
    end

    status = Dict(
        "name"              => config["name"],
        "uuid"              => config["uuid"],
        "total_value"       => config["total_value"],
        "bought"            => bought,
        "profit"            => profit,
        "today_profit"      => today_profit,
        "daily_avg"         => daily_avg,
        "trade_num"         => trade_num,
        "trade_num_d"       => round(trade_num_d, digits=1),
        "trade_num_h"       => round(trade_num_h, digits=1),
        "time"              => string(now()),
        "current_gridpoint" => current_gridpoint,
        "gridpoint_num"     => length(trade_grid),
        "to_sell"           => to_sell,
        "progress_line"     => progress_line,
        "rendite"           => rendite,
        "profit_thes"       => profit_thes,
        "thes_factor"       => config["thes_factor"],
        "current_volume"    => current_value
    )

    if haskey(next_gridpoints, "BUY")
        status["next_buy"] = next_gridpoints["BUY"]
        progress_line =
        round((current_ticker - parse(Float64, next_gridpoints["BUY"]))/((parse(Float64, next_gridpoints["SELL"]) - parse(Float64, next_gridpoints["BUY"]))/100))
    else
        progress_line = 0
    end
    if haskey(next_gridpoints, "SELL")
        status["next_sell"] = next_gridpoints["SELL"]
    end

    open(joinpath(bot_path, "status.json"), "w") do f
        JSON.print(f, status, 4)
    end

    return profit, trade_num, sold
end
export update_status

function create_grid_super_status()
    user_paths = glob(joinpath(Pinguu.base_dir , "grid_bots/user/*"));
    super_status = Dict()
    for user_path in user_paths
        status_paths = glob(joinpath(user_path, "*/status.json"))
        for status_path in status_paths
            status = JSON.parsefile(status_path)
            config = JSON.parsefile(joinpath(split(status_path, "status.json")[1], "config.json"))
            status["config"] = config

            user_id = split(user_path, "/")[end]
            !haskey(super_status, user_id) ? super_status[user_id] = Dict("active" => [], "inactive" => []) : ""
            if config["status"] == "active"
                push!(super_status[user_id]["active"], status)
            else
                push!(super_status[user_id]["inactive"], status)
            end
        end
    end
    open(joinpath(Pinguu.base_dir, "stats/grid_super_status.json"), "w") do f
        JSON.print(f, super_status, 4)
    end
end
export create_grid_super_status

function create_bots_super_status()
    user_paths = glob(joinpath(Pinguu.base_dir , "users/*"));
    super_status = Dict()
    for user_path in user_paths
        status_paths = glob(joinpath(user_path, "*/status.json"))
        for status_path in status_paths
            status = JSON.parsefile(status_path)
            config = JSON.parsefile(joinpath(split(status_path, "status.json")[1], "config.json"))
            status["config"] = config
#             delete!(status, "plotting")

            user_id = split(user_path, "/")[end]
            !haskey(super_status, user_id) ? super_status[user_id] = Dict("active" => [], "inactive" => []) : ""
            if config["status"] in ["active", "pause", "stop"]
                push!(super_status[user_id]["active"], status)
            else
                push!(super_status[user_id]["inactive"], status)
            end
        end
    end
    open(joinpath(Pinguu.base_dir, "stats/bots_super_status.json"), "w") do f
        JSON.print(f, super_status, 4)
    end
end
export create_bots_super_status


text:
"""
    simpledownup

Compute the mesh topology structures that correspond to the downward and upward
adjacency between the tetrahedra and the vertices. Memory is used only to store
the incidence relation (3, 0), (0, 3) and the locations of the vertices.
"""
module simpledownup
using StaticArrays
using MeshCore: nvertices, nshapes, manifdim, attribute
using MeshCore: ir_skeleton, ir_bbyfacets, ir_transpose
using MeshSteward: vtkwrite, T4block
using BenchmarkTools
using Test

include("usedbytes.jl")

function test()
    n = 3
    membytes = 0; summembytes = 0
    @info "Initial (3, 0)"
    @time connectivity = T4block(1.0, 2.0, 3.0, n*7, n*9, n*10, :a; intbytes = 4)
    ir30 = connectivity
    @show "($(manifdim(ir30.left)), $(manifdim(ir30.right)))"
    @show (nshapes(ir30.left), nshapes(ir30.right))
    @show membytes = usedbytes(ir30._v)
    summembytes += membytes
    geom = attribute(ir30.right, "geom")
    @show membytes = usedbytes(geom.v)
    summembytes += membytes

    @info "Transpose. (0, 3)"
    @time tr03 = ir_transpose(ir30)
    @show "($(manifdim(tr03.left)), $(manifdim(tr03.right)))"
    @show (nshapes(tr03.left), nshapes(tr03.right))
    @show membytes = usedbytes(tr03._v)
    summembytes += membytes

    # Print the total number of megabytes used to store the database
    @show summembytes/2^20

    # vtkwrite("speedtest1", connectivity)
    true
end
end
using .simpledownup
# using BenchmarkTools
# @btime simpledownup.test()
simpledownup.test()


text:
"""
	Rocket
"""

abstract type Rocket end

struct RocketNominal{T} <: Rocket
	n::Int
	m::Int
	d::Int

	g::T  # gravity
	m1::T # mass
	l1::T # length
	J::T  # inertia
end

struct RocketSlosh{T} <: Rocket
	n::Int
	m::Int
	d::Int

	g::T  # gravity
	m1::T # mass
	l1::T # length
	J::T  # inertia

	m2::T # mass of fuel
	l2::T # length of fuel pendulum
	l3::T # length from rocket COM to fuel pendulum
end

g = 9.81 # gravity
m1 = 1.0 # mass of rocket
l1 = 0.5 # length from COM to thruster
J = 1.0 / 12.0 * m1 * (2.0 * l1)^2 # inertia of rocket

m2 = 0.1 # mass of pendulum
l2 = 0.1 # length from COM to pendulum
l3 = 0.1 # length of pendulum

# models
m = 2
d = 0

# nominal model
n_nominal = 6
model_nominal = RocketNominal(n_nominal, m, d, g, m1, l1, J)

# slosh model
n_slosh = 8
model_slosh = RocketSlosh(n_slosh, m, d, g, m1 - m2, l1, J, m2, l2, l3)

function k_thruster(model::Rocket, x)
	y, z, θ = x[1:3]
	px = y + model.l1 * sin(θ)
	pz = z - model.l1 * cos(θ)

	return [py;
		    pz]
end

function jacobian_thruster(model::RocketNominal, x)
	y, z, θ = x[1:3]

	return [1.0 0.0 model.l1 * cos(θ);
			0.0 1.0 model.l1 * sin(θ)]
end

function jacobian_thruster(model::RocketSlosh, x)
	y, z, θ = x[1:3]

	return [1.0 0.0 model.l1 * cos(θ) 0.0;
			0.0 1.0 model.l1 * sin(θ) 0.0]
end

function k_mass(model::RocketSlosh, x)
	yp = x[1] + l2 * sin(x[3]) + l3 * sin(x[4])
	zp = x[2] - l2 * cos(x[3]) - l3 * cos(x[4])

	return [yp;
			zp]
end

function jacobian_mass(model::RocketSlosh, x)
	return [1.0 0.0 l2 * cos(x[3]) l3 * cos(x[4]);
	        0.0 1.0 l2 * sin(x[3]) l3 * sin(x[4])]
end

function lagrangian(model::RocketNominal, q, q̇)
	return (0.5 * model.m1 * (q̇[1]^2 + q̇[2]^2.0) + 0.5 * model.J * q̇[3]^2.0
			- model.m1 * model.g * q[2])
end

function lagrangian(model::RocketSlosh, q, q̇)
	zp = k_mass(model, q)[2]
	vp = jacobian_mass(model, q) * q̇

	return (0.5 * model.m1 * (q̇[1]^2 + q̇[2]^2) + 0.5 * model.J * q̇[3]^2.0
			- model.m1 * model.g * q[2]
			+ 0.5 * model.m2 * vp' * vp
			- model.m2 * model.g * zp)
end

function dLdq(model::Rocket, q, q̇)
	Lq(x) = lagrangian(model, x, q̇)
	ForwardDiff.gradient(Lq, q)
end

function dLdq̇(model::Rocket, q, q̇)
	Lq̇(x) = lagrangian(model, q, x)
	ForwardDiff.gradient(Lq̇, q̇)
end

function f(model::Rocket, x, u, w)
	nq = convert(Int, floor(model.n / 2.0))
	q = x[1:nq]
	q̇ = x[nq .+ (1:nq)]
	tmp_q(z) = dLdq̇(model, z, q̇)
	tmp_q̇(z) = dLdq̇(model, q, z)
	[q̇;
	 ForwardDiff.jacobian(tmp_q̇,q̇) \ (-1.0 * ForwardDiff.jacobian(tmp_q, q) * q̇
	 	+ dLdq(model, q, q̇)
		+ jacobian_thruster(model, q)' * u[1:2])]
end

function state_output(model::RocketSlosh, x)
	x[collect([1, 2, 3, 5, 6, 7])]
end

function state_output_idx(mode::RocketSlosh, idx)
	idx[collect([1, 2, 3, 5, 6, 7])]
end

function visualize!(vis, model::Rocket, x;
       Δt = 0.1, r_rocket = 0.1, r_pad = 0.25)

	obj_rocket = joinpath(pwd(), "src/models/rocket/space_x_booster.obj")
	mtl_rocket = joinpath(pwd(), "src/models/rocket/space_x_booster.mtl")

	rkt_offset = [4.0, -6.35, 0.2]
	ctm = ModifiedMeshFileObject(obj_rocket, mtl_rocket, scale = 1.0)
	setobject!(vis["rocket"], ctm)
	settransform!(vis["rocket"], compose(Translation((x[T][1:3] + rkt_offset)...),
		LinearMap(RotZ(-pi) * RotX(pi / 2.0))))

	obj_platform = joinpath(pwd(), "src/models/rocket/space_x_platform.obj")
	mtl_platform = joinpath(pwd(), "src/models/rocket/space_x_platform.mtl")

	ctm_platform = ModifiedMeshFileObject(obj_platform, mtl_platform, scale = 1.0)
	setobject!(vis["platform"], ctm_platform)
	settransform!(vis["platform"], compose(Translation(0.0, 0.0, 0.0),
		LinearMap(RotZ(pi) * RotX(pi / 2.0))))

   	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(x)
        MeshCat.atframe(anim,t) do
			settransform!(vis["rocket"],
				compose(Translation(x[t][1] + rkt_offset[1],
					0.0 + rkt_offset[2],
					x[t][2] + rkt_offset[3]),
				LinearMap(RotY(-1.0 * x[t][3]) * RotZ(pi) * RotX(pi / 2.0))))
        end
    end

    # settransform!(vis["/Cameras/default"], compose(Translation(0.0, 0.0, 0.0),
	# 	LinearMap(RotZ(pi / 2.0))))

    MeshCat.setanimation!(vis, anim)
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using Test
using SortingLab

x = [1,2,missing, 3]

@test isequal(fsortperm(x), sortperm(x))


text:
abstract type AbstractRaggedMatrix{T,I} <: AbstractMatrix{T} end

# struct RaggedMatrix{T,I} <: AbstractRaggedMatrix{T,I}
#     data::Vector{T} # length(data) == sum(column_lengths) == column_offsets[end]
#     column_offsets::Vector{I} # length(column_offsets) == size(A,2)
#     column_lengths::Vector{I} # length(column_lengths) == size(A,2)
#     nrow::Int
# end

# Column offsets go from 0 (col 1 offset), col 2 offset,...,2nd last col offset
struct RaggedMatrix{T,I,VI<:AbstractVector{I},VT<:AbstractVector{T}} <: AbstractRaggedMatrix{T,I}
    data::VT # length(data) == sum(column_lengths) == column_offsets[end]
    column_offsets::VI # length(column_offsets) == size(A,2)
    column_lengths::VI # length(column_lengths) == size(A,2)
    nrow::Int
end
struct FixedSizeRaggedMatrix{M,N,P,T,I} <: AbstractRaggedMatrix{T,I}
    data::FixedSizeVector{P,T,P}
    column_offsets::FixedSizeVector{N,I,N}
    column_lengths::FixedSizeVector{N,I,N}
end
function FixedSizeRaggedMatrix{M,N,P}(
    data::AbstractVector{T},
    column_offsets::AbstractVector{I},
    column_lengths::AbstractVector{I}
) where {M,N,P,T,I}
    FixedSizeRaggedMatrix{M,N,P,T,I}(
        FixedSizeVector{P,T,P}(data),
        FixedSizeVector{N,I,N}(column_offsets),
        FixedSizeVector{N,I,N}(column_lengths)
    )
end
# const FixedSizeRaggedMatrix{T,I,NC,ND} = RaggedMatrix{T,I,FixedSizeVector{NC,I,NC},FixedSizeVector{ND,T,ND}}
nrow(A::AbstractRaggedMatrix) = A.nrow
ncol(A::AbstractRaggedMatrix) = length(A.column_lengths)
nrow(::FixedSizeRaggedMatrix{M}) where {M} = M
ncol(::FixedSizeRaggedMatrix{M,N}) where {M,N} = N

Base.size(A::AbstractRaggedMatrix) = (nrow(A),ncol(A))
number_not_structurally_zero(A::AbstractRaggedMatrix) = length(A.data)
number_not_structurally_zero(::FixedSizeRaggedMatrix{M,N,P}) where {M,N,P} = P

@Base.propagate_inbounds Base.getindex(A::AbstractRaggedMatrix, i::Integer) = A.data[i]
@inline function Base.getindex(A::AbstractRaggedMatrix{T,I}, i::Integer, j::Integer) where {T,I}
    @boundscheck begin
        ( (i > nrow(A)) || ( (j > ncol(A)) || (min(i,j) < 1) ) ) && PaddedMatrices.ThrowBoundsError("Tried to index array of size $(size(A)) with index ($i,$j).")
    end
    @inbounds begin
        col_j_offset = A.column_offsets[j] #j == one(j) ? zero(I) : A.column_offsets[j-1]#; col_j_length = A.column_descriptions[j,2]
        col_j_nextoffset = j == length(A.column_offsets) ? Base.unsafe_trunc(I, length(A.data)) : A.column_offsets[j+1]
        ind = col_j_offset + Base.unsafe_trunc(I, i)
        col_j_nextoffset >= ind ? A.data[ind] : zero(T)
    end
end

### Iteration returns (row, column, value)
@inline function Base.iterate(A::AbstractRaggedMatrix)
    (1, 1, 1), (2, 2, 1)
end
@inline function Base.iterate(A::AbstractRaggedMatrix, (k,i,j)::NTuple{3,Int})
    if i > A.column_lengths[j]
        j == ncol(A) && return nothing
        i, j = 1, j+1
    end
    # k = 1,...,P; i = 1,...,M; j = 1,...,N
    (k, i, j), (k+1, i+1, j)
end

function expand_by_x_quote!(q, N, ncol, nrl, T, sp::Bool, def::Bool=true)
    push!(q.args, ind = 1)
    loop_body = quote end
    a_getindex = quote end
    return_expression = Expr(:tuple)
    for n in 1:N
        if N > 1
            asym = Symbol(:a_,n)
            push!(q.args, :($asym = @inbounds a[$n]))
        else
            asym = :a
        end
        asymi = Symbol(:a_i_,n)
        bsym = Symbol(:b_,n)
        def && push!(q.args, pointer_vector_expr( bsym, nrl, T, sp, :sptr) )
        push!(a_getindex.args, :($asymi = $asym[i]))
        push!(loop_body.args, :($bsym[ind] = $asymi))
        push!(return_expression.args, bsym)
    end
    return_expression = sp ? :((sptr,$return_expression)) : return_expression
    loop_quote = quote
        @inbounds for i in 1:$ncol
            $a_getindex
            nr = column_lengths[i]
            for r in one(nr):nr
                $loop_body
                ind += 1
            end
        end
        $return_expression
    end
    push!(q.args, loop_body)
    q
end
function contract_by_x_quote!(q, N, ncol, nrl, T, sp::Bool)
    push!(q.args, ind = 1)
    loop_body = quote end
    bset_quote = quote end
    loop_suffix = quote end
    return_expression = Expr(:tuple)
    for n in 1:N
        asym = Symbol(:a_,n)
        # asymi = Symbol(:a_i_,n)
        bsym = Symbol(:b_,n)
        bsymi = Symbol(:b_i_,n)
        push!(q.args, pointer_vector_expr( bsym, ncol, T, sp, :sptr) )
        push!(bset_quote.args, :($bsymi = zero(T)))
        push!(loop_body.args, :($bsymi += $asym[i]))
        push!(loop_suffix.args, :($bsym[ind] = $bymi))
        push!(return_expression.args, bsym)
    end
    return_expression = sp ? :((sptr,$return_expression)) : return_expression
    loop_quote = quote
        @inbounds for i in 1:$ncol
            $bset_quote
            nr = column_lengths[i]
            for r in one(nr):nr
                $loop_body
                ind += 1
            end
            $loop_suffix
        end
        $return_expression
    end
    push!(q.args, loop_body)
    q
end
function expand_contract_prequote(VI, VT)
    q = quote column_lengths = x.column_lengths end
    if VI <: AbstractFixedSizeVector
        ncol = full_length(VI)
    else
        ncol = :ncol
        push!(q.args, :($ncol = :(ncol(x))))
    end
    if VT <: AbstractFixedSizeVector
        nrl = full_length(VI)
    else
        nrl = :nrl
        push!(q.args, :($nrl = :(number_not_structurally_zero(x))))
    end
    ncol, nrl, q
end

@generated function expand_by_x(x::RaggedMatrix{T,I,VI,VT}, a::Vararg{<:AbstractVector,N}) where {N,T,I,VI,VT}
    ncol, nrl, q = expand_contract_prequote(VI, VT)
    expand_by_x_quote!(q, N, ncol, nrl, T, false)
end
@generated function contract_by_x(x::RaggedMatrix{T,I,VI,VT}, a::Vararg{<:AbstractVector,N}) where {N,T,I,VI,VT}
    ncol, nrl, q = expand_contract_prequote(VI, VT)
    contract_by_x_quote!(q, N, ncol, nrl, T, false)
end
function expand_by_x!(b::AbstractVector, x::RaggedMatrix{T,I,VI,VT}, a::AbstractVector) where {T,I,VI,VT}
    ncol, nrl, q = expand_contract_prequote(VI, VT)
    contract_by_x_quote!(q, 1, ncol, nrl, T, false, false)
end
@generated function expand_by_x(sptr::StackPointer, x::RaggedMatrix{T,I,VI,VT}, a::Vararg{<:AbstractVector,N}) where {N,I,T,VI,VT}
    ncol, nrl, q = expand_contract_prequote(VI, VT)
    expand_by_x_quote!(q, N, ncol, nrl, T, true)
end
@generated function contract_by_x(sptr::StackPointer, x::RaggedMatrix{T,I,VI,VT}, a::Vararg{<:AbstractVector,N}) where {N,I,T,VI,VT}
    ncol, nrl, q = expand_contract_prequote(VI, VT)
    contract_by_x_quote!(q, N, ncol, nrl, T, true)
end


text:
# Naive implementation

# using Simulation_Env_module

mutable struct PrioritizedReplayBuffer
    max_size::Int64
    batch_size::Int64
    rng::AbstractRNG
    α::Float64
    β::Float64
    ϵ::Float64
    _curr_size::Int64
    _idx::Int64
    _priorities::Vector{Float64}
    _experience::Vector{DQExperience}

    _s_batch::Array{Float64}
    _a_batch::Vector{Int64}
    _r_batch::Vector{Float64}
    _sp_batch::Array{Float64}
    _done_batch::Vector{Bool}
    _weights_batch::Vector{Float64}

    function PrioritizedReplayBuffer(env::AbstractEnvironment,
                                    max_size::Int64,
                                    batch_size::Int64;
                                    rng::AbstractRNG = MersenneTwister(0),
                                    α::Float64 = 0.6,
                                    β::Float64 = 0.4,
                                    ϵ::Float64 = 1e-3)
        s_dim = obs_dimensions(env)
        experience = Vector{DQExperience}(max_size)
        priorities = Vector{Float64}(max_size)
        _s_batch = zeros(batch_size, s_dim...)
        _a_batch = zeros(Int64, batch_size)
        _r_batch = zeros(batch_size)
        _sp_batch = zeros(batch_size, s_dim...)
        _done_batch = zeros(Bool, batch_size)
        _weights_batch = zeros(Float64, batch_size)
        return new(max_size, batch_size, rng, α, β, ϵ, 0, 1, priorities, experience,
                   _s_batch, _a_batch, _r_batch, _sp_batch, _done_batch, _weights_batch)
    end
end

is_full(r::PrioritizedReplayBuffer) = r._curr_size == r.max_size

max_size(r::PrioritizedReplayBuffer) = r.max_size

function add_exp!(r::PrioritizedReplayBuffer, expe::DQExperience, td_err::Float64=abs(expe.r))
    @assert td_err + r.ϵ > 0.
    priority = (td_err + r.ϵ)^r.α
    r._experience[r._idx] = expe
    r._priorities[r._idx] = priority
    r._idx = mod1((r._idx + 1),r.max_size)
    if r._curr_size < r.max_size
        r._curr_size += 1
    end
end

function update_priorities!(r::PrioritizedReplayBuffer, indices::Vector{Int64}, td_errors::Float64)
    @assert td_errors + r.ϵ .> 0.
    new_priorities = (td_err + r.ϵ)^r.α
    r._priorities[indices] = new_priorities
end

function StatsBase.sample(r::PrioritizedReplayBuffer)
    @assert r._curr_size >= r.batch_size
    @assert r.max_size >= r.batch_size # could be checked during construction
    sample_indices = sample(r.rng, 1:r._curr_size, Weights(r._priorities[1:r._curr_size]), r.batch_size, replace=false)
    @assert length(sample_indices) == size(r._s_batch)[1]
    for (i, idx) in enumerate(sample_indices)
        r._s_batch[Base.setindex(indices(r._s_batch), i, 1)...] = r._experience[idx].s
        r._a_batch[i] = r._experience[idx].a
        r._r_batch[i] = r._experience[idx].r
        r._sp_batch[Base.setindex(indices(r._sp_batch), i, 1)...] = r._experience[idx].sp
        r._done_batch[i] = r._experience[idx].done
        r._weights_batch[i] = r._priorities[idx]
    end
    r._weights_batch = (r._weights_batch.*r._curr_size).^-r.β./maximum(r._weights_batch)
    return r._s_batch, r._a_batch, r._r_batch, r._sp_batch, r._done_batch, sample_indices, r._weights_batch
end

function populate_replay_buffer!(replay::PrioritizedReplayBuffer, env::AbstractEnvironment;
                                 max_pop::Int64=replay.max_size, max_steps::Int64=2000)
    render_option = env.render
    env.render = false
    println("Start populating PRIORITIZED replay buffer...")
    # println("Using closest intruder detection...")
    if !env.random_populate_replay_buffer
        println("Using multi intruder detection ...")
    else
        println("Using random ...")
    end
    o = reset(env)
    done = false
    step = 0
    for t = 1:(max_pop - replay._curr_size)
        if !env.random_populate_replay_buffer
            q_low_fi = lowfi_values(env, o) # multi-threat
            # q_low_fi = get_qvals(env.low_fi_policy, convert_ac_state(env.ego_agent.state)) # closest
            action = actions(env)[indmax(q_low_fi)]
        else
            action = sample_action(env)
        end
        ai = action_index(env.problem, action)
        op, rew, _, done, info = step!(env, action)
        exp = DQExperience(o, ai, rew, op, done)
        add_exp!(replay, exp, abs(rew)) # assume initial td error is r
        o = op
        step += 1
        if done || step >= max_steps
            o = reset(env)
            done = false
            step = 0
        end
    end
    @assert replay._curr_size >= replay.batch_size
    println("Done populating replay buffer!")
    env.render = render_option
end

max_stars_repo_name: parachutel/DeepCorrectionCAS