text:
function find_max_subarray(A, start, finish)
    max_sum = -Inf
    sum = 0
    max_index = -1

    start > finish ? step = -1 : step = 1

    for i in range(start, finish, step = step)
        sum += A[i]
        if sum > max_sum
            max_sum = sum
            max_index = i
        end
    end

    return (max_index, max_sum)
end

function find_max_crossing_subarray(A, low, mid, high)
    max_left, left_sum = find_max_subarray(A, mid, low)
    max_right, right_sum = find_max_subarray(A, mid+1, high)

    return (max_left, max_right, left_sum + right_sum)
end

function find_maximum_subarray(A, low, high)
    # base case only one elem
    if high == low
        return (low, high, A[low])
    end

    mid = floor(Int, (low + high) / 2)

    left_low, left_high, left_sum = find_maximum_subarray(A, low, mid)
    right_low, right_high, right_sum = find_maximum_subarray(A, mid+1, high)
    x_low, x_high, x_sum = find_max_crossing_subarray(A, low, mid, high)

    if left_sum ≥ right_sum && left_sum ≥ x_sum
        return (left_low, left_high, left_sum)
    elseif right_sum ≥ left_sum && right_sum ≥ x_sum
        return (right_low, right_high, right_sum)
    else
        return (x_low, x_high, x_sum)
    end
end

using Test

@testset "max sub array" begin
    input = [13 -3 -25 20 -3 -16 -23 18 20 -7 12 -5 -22 15 -4 7]
    output = (8, 11, 43)
    @test find_maximum_subarray(input, 1, length(input)) == output
end

max_stars_repo_name: blallen/CodePractice

max_issues_repo_name: blallen/CodePractice

max_forks_repo_name: blallen/CodePractice


text:
"""
    logreg(X, y)

    Objective function of Logistic Regression on data matrix X and label vector y.
"""
logreg(X::Array{Float64,2}, y::Array{Float64,1}, λ::Float64=0.0) = Objective(
    w::Array{Float64,1} -> _logreg_f(X, y, w, λ),  # f
    w::Array{Float64,1} -> _logreg_g(X, y, w, λ),  # g
    f_b = (w::Array{Float64,1}, b::Array{Int,1}) -> _logreg_f(X, y, w, λ, b),
    g_b = (w::Array{Float64,1}, b::Array{Int,1}) -> _logreg_g(X, y, w, λ, b),
    dim = length(y)
)

function _logreg_f(X::Array{Float64,2}, y::Array{Float64,1}, w::Array{Float64,1},
                   λ::Float64, b::Array{Int,1}=Int[])
    reg = 0.0
    if λ > 0.0
        reg = λ * vecnorm(w, 1)
    end
    if length(b) == 0
        return sum( log(1 + exp(-y .* X*w)) )[1] / length(y) + reg
    else
        return sum( log(1 + exp(-y[b] .* X[b,:]*w)) )[1] / length(b) + reg
    end
end

function _logreg_g(X::Array{Float64,2}, y::Array{Float64,1}, w::Array{Float64,1},
                   λ::Float64, b::Array{Int,1}=Int[])
    reg = 0.0
    if λ > 0.0
        reg = λ * vecnorm(w, 1)
    end
    if length(b) == 0
        return - vec(sum( (y ./ (1 + exp(y .* X*w))) .* X , 1)) ./ length(y) .+ reg
    else
        return - vec(sum( (y[b] ./ (1 + exp(y[b] .* X[b,:]*w))) .* X[b,:] , 1)) ./ length(b) .+ reg
    end
end

"""
    logreg_predict(w, X)

    Predict the examples in X with Logistic Regression under parameter vector w.
"""
logreg_predict(w::Array{Float64,1}, X::Array{Float64,2}) = sign(vec(1 ./ (1 + exp(-X*w))) -.5)

ext: jl

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
include("../LoadGraph/LoadGraph.jl")
using .LoadGraph

function main(dir_path="../../cleaned_data", dtype::Type=UInt32)
    for file in filter(x->endswith(x, ".txt"), readdir(dir_path)) if !startswith(file, "friend")
        println(file[1:end-4])
        try
            @time load_adjacency_list("$dir_path/$file", dtype)
        catch OutOfMemoryError
            println("OOM")
        end
        try
            @time load_edge_list("$dir_path/$file", dtype)
        catch OutOfMemoryError
            println("OOM")
        end
        try
            @time load_adjacency_matrix("$dir_path/$file")
        catch OutOfMemoryError
            println("OOM")
        end
        println()
    end end
end

main()

max_stars_repo_name: quentin-burthier/GrandsGraphes

max_issues_repo_name: quentin-burthier/GrandsGraphes

max_forks_repo_name: quentin-burthier/GrandsGraphes


text:
using Test
using Revise
using BondGraphs
using Graphs
using ModelingToolkit
using Catalyst
using RecipesBase

@testset "BondGraphs.jl" begin
    @testset "Graph functions" begin include("graphfunctions_tests.jl") end
    @testset "Construction" begin include("construction_tests.jl") end
    @testset "Equations" begin include("equation_tests.jl") end
    @testset "Simulations" begin include("simulation_tests.jl") end
    @testset "Conversions" begin include("conversion_tests.jl") end
    @testset "Miscellaneous" begin include("misc_tests.jl") end
    @testset "Modules" begin include("module_tests.jl") end
end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
# Weighted statistics

function wmean{T<:Number,W<:Real}(v::AbstractArray{T}, w::AbstractArray{W})
    sv = 0.0
    sw = 0.0
    for i in 1:length(v)
    	wi = w[i]
        sv += v[i] * wi
        sw += wi
    end
    return sv / sw
end


text:
function newtondist!(yout::Vector{Float64},iters::Vector{Int64},Ws::Float64,W2est::Float64,
    Pn::Float64,C::Float64,R0::Float64,R1::Float64,beta::Float64,c0::Float64,A0::Float64,
    W1::Float64,rho::Float64,f::Function,J::Function,maxiter::Int16,epsJ::Float64,
    epsN::Float64,maxval::Float64,h::Float64)

    # initial guess
    x0 = zeros(1);
    xx = zeros(1);
    x0[1] = W2est;

    # non-dimensionalize
    x0[1] /= Ws;

    # setup for iterations
    xx .= x0;
    x = zeros(1);
    N = 1;

    # allocation
    JJ = zeros(1);
    D = zeros(1);
    d = zeros(1);
    fvec = zeros(1);
    g = zeros(1);

    # max step size for line searches
    stpmax = 100*max(sqrt(norm(xx)),1);
    # println("Maximum step size: $stpmax")

    # useful parameter combinations
    a1 = 2*rho*(1+R1*C/h);
    a2 = (R0+R1+R0*R1*C/h)*(2*rho/beta)^2;

    # Newton iterations
    while N <= maxiter
        # determine Jacobian, check invertibility
        # println("Distal Jacobian:")
        J(JJ,xx,C,R0,R1,Ws,W1,rho,beta,c0,a1,a2,h);
        maximum!(d,abs.(JJ))
        D[1] = 1/d[1];
        # println(JJ)
        # println(D)
        # println(D*JJ)
        # println(cond(D*JJ))
        # println(det(D*JJ))
        if abs.((D.*JJ)[1]) < epsJ
            println(JJ)
            println(D.*JJ)
            error("Distal Newton Jacobian is singular.");
        end
        # compute gradient of line search objective function
        # println("Objective function:")
        f(fvec,xx,D,C,R0,R1,Ws,W1,rho,beta,c0,a1,a2,Pn,h);
        # println(f(xx,system,n,state))
        # println(fvec)
        fold = 0.5*dot(fvec,fvec);
        g = [(D.*JJ)'*fvec];
        # compute newton step
        # println(inv(D*JJ))
        # println("Newton step:")
        s = -inv((D.*JJ)[1])*fvec;
        # line search to update state vector
        # println("Line search:")
        fn, xn, check = LegHemodynamics.linedist(xx,fold,g,s,stpmax,f,J,
            maxiter,Ws,W1,rho,beta,C,R0,R1,c0,a1,a2,Pn,h,JJ,d,fvec,D);
        # println(xn[1]*vs)
        # check if sufficiently close to root
        # println("Root checking:")
        maximum!(d,abs.(JJ))
        D[1] = 1/d[1];
        f(fvec,xn,D,C,R0,R1,Ws,W1,rho,beta,c0,a1,a2,Pn,h);
        # println(fvec)
        if norm(fvec) <= 100*epsN
            x .= xn;
            x[1] = x[1]*Ws;
            # println("Converged estimate for W2: $(x[1])")
            # println(inv(D*JJ))
            # println(fvec)
            iters[1] += N;
            break
        end
        if check
            test = 0.
            den = max(fn,0.5*length(xn));
            for i = 1:length(xn)
                temp = abs(g[i])*max(abs(xn[i]),1.)/den;
                if temp > test
                    test = temp;
                end
            end
            if test < 1e-6
                check = true;
                println("Warning: gradient of objective function close to zero.")
            else
                check = false;
            end
        end
        # check for divergence
        if norm(fvec,Inf) >= maxval
            iters[1] += N;
            println(xn)
            println(fvec)
            error("Newton iteration diverged.");
        end
        N+=1;
        xx = xn;
        if N == maxiter
            println(JJ)
            # println(D)
            # println(D*JJ)
            println(xn)
            println(fvec)
            println(norm(fvec))
            error("Newton iteration failed to converge.");
        end
    end

    # update based on converged solution
    # println("Update:")
    yout[1] = (2*rho/beta)^2*(0.125*(W1-x[1])+c0)^4;
    yout[2] = yout[1]*0.5*(W1+x[1]);
    yout[3] = beta*(sqrt(yout[1]) - sqrt(A0)) - R0*yout[2];
    yout[4] = yout[3]/R1;
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
@inline function always_assert(condition::Bool, message::String)
    if !condition
        throw(AlwaysAssertionError(message))
    end
    return nothing
end


text:
import Core: SSAValue, GotoNode, Compiler
import Core: Typeof
import Core.Compiler: CodeInfo, IRCode, CFG, BasicBlock, Argument, ReturnNode,
  NullLineInfo, just_construct_ssa, compact!, NewNode, InferenceState, OptimizationState,
  GotoIfNot, PhiNode, PiNode, StmtRange, IncrementalCompact, insert_node!, insert_node_here!,
  compact!, finish, DomTree, construct_domtree, dominates, userefs, widenconst, types, verify_ir
using InteractiveUtils: typesof

afterphi(ir, loc) = ir.stmts[loc] isa PhiNode ? afterphi(ir, loc+1) : loc

insert_blockstart!(ir::IRCode, pos, typ, val) =
  insert_node!(ir, afterphi(ir, ir.cfg.blocks[pos].stmts[1]), typ, val)

function insert_blockend!(ir::IRCode, pos, typ, val)
  i = first(ir.cfg.blocks[pos].stmts)
  j = last(ir.cfg.blocks[pos].stmts)
  if !(ir.stmts[j] isa Union{GotoNode,GotoIfNot,ReturnNode})
    return insert_node!(ir, j, typ, val, true)
  end
  while j > i && ir.stmts[j-1] isa Union{GotoNode,GotoIfNot,ReturnNode}
    j -= 1
  end
  insert_node!(ir, j, typ, val)
end

function finish_dc(ic::IncrementalCompact)
  Compiler.non_dce_finish!(ic)
  return Compiler.complete(ic)
end

function _compact!(code::IRCode)
    compact = IncrementalCompact(code)
    foreach(x -> nothing, compact)
    return finish_dc(compact), compact.ssa_rename
end

function argmap(f, @nospecialize(stmt))
    urs = userefs(stmt)
    for op in urs
        val = op[]
        if isa(val, Argument)
            op[] = f(val)
        end
    end
    return urs[]
end

exprtype(ir::IRCode, x::Argument) = widenconst(ir.argtypes[x.n])
exprtype(ir::IRCode, x::SSAValue) = widenconst(types(ir)[x])
exprtype(ir::IRCode, x::GlobalRef) = isconst(x.mod, x.name) ? Typeof(getfield(x.mod, x.name)) : Any
exprtype(ir::IRCode, x::QuoteNode) = Typeof(x.value)
# probably can fall back to any here
exprtype(ir::IRCode, x::Union{Type,Number,Nothing,Tuple,Function,Val,String}) = Typeof(x)
exprtype(ir::IRCode, x::Expr) = error(x)

rename(x, m) = x
rename(x::SSAValue, m) = m[x.id]
rename(xs::AbstractVector, m) = map(x -> rename(x, m), xs)
rename(xs::AbstractSet, m) = Set(rename(x, m) for x in xs)
rename(d::AbstractDict, m) = Dict(k => rename(v, m) for (k, v) in d)

function usages(ir)
  us = Dict()
  for i = 1:length(ir.stmts), u in userefs(ir.stmts[i])
    push!(get!(us, u[], []), SSAValue(i))
  end
  return us
end

function blockidx(ir, i::Integer)
  i = findlast(x -> x <= i, ir.cfg.index)
  i == nothing ? 1 : i+1
end

blockidx(ir, i::SSAValue) = blockidx(ir, i.id)

if VERSION > v"1.1.0-DEV.560"
  Base.range(b::BasicBlock) = b.stmts.start:b.stmts.stop
else
  Base.range(b::BasicBlock) = b.stmts.first:b.stmts.last
end

xcall(mod::Module, f::Symbol, args...) = Expr(:call, GlobalRef(mod, f), args...)
xcall(f::Symbol, args...) = xcall(Base, f, args...)

const unreachable = ReturnNode()

# Dominance frontiers

function domfront(cfg, dt = construct_domtree(cfg))
  fronts = [Set{Int}() for _ in cfg.blocks]
  for b = 1:length(cfg.blocks)
    length(cfg.blocks[b].preds) >= 2 || continue
    for p in cfg.blocks[b].preds
      runner = p
      while runner != dt.idoms[b]
        runner == b && break
        push!(fronts[runner], b)
        runner = dt.idoms[runner]
      end
    end
  end
  return fronts
end


text:
using LarSurf

@testset "Convert 2D array to list and back" begin
    hexV, hexFVarr = LarSurf.hexagon()
    hexFVll = LarSurf.array2ll(hexFVarr)
    hexFVarr2 = LarSurf.ll2array(hexFVll)
    @test hexFVarr[2,1] == hexFVarr2[2,1]
    @test hexFVarr[4,2] == hexFVarr2[4,2]
end

@testset "Test characteristicMatrix implementations" begin
    hexV, hexFVarr = LarSurf.hexagon()
    hexFVll = LarSurf.array2ll(hexFVarr)
    mat0 = LarSurf.characteristicMatrix_for_loop(hexFVarr, size(hexV, 1))
    mat1 = LarSurf.characteristicMatrix_push(hexFVll )
    mat2 = LarSurf.characteristicMatrix_set( hexFVll )
    mat3 = LarSurf.characteristicMatrix_push(hexFVarr )
    mat4 = LarSurf.characteristicMatrix_set( hexFVarr)

    @test LarSurf.match_arr(mat0, mat1)
    @test LarSurf.match_arr(mat0, mat2)
    @test LarSurf.match_arr(mat0, mat3)
    @test LarSurf.match_arr(mat0, mat4)
end


text:
# https://www.chessprogramming.org/Encoding_Moves
# Idea is to use a 16-bit Int to encode the move.
# 6 bits for "from", 6 bits for "to", 4 bits for "flags"
const FLAGS = [:__NORMAL_MOVE, :__DOUBLE_PAWN, :__KING_CASTLE, :__QUEEN_CASTLE, :__ENPASS,
    :__KNIGHT_PROMO, :__BISHOP_PROMO, :__ROOK_PROMO, :__QUEEN_PROMO]

for (num, flag) in enumerate(FLAGS)
    @eval const $flag = UInt16($num - 1) << 12
end

"""
    updatecastling!(board::Board, sqr_from::Integer, sqr_to::Integer)

Update the castling rights of the `board`, given a move is played from `sqr_from` to `sqr_to`.
"""
function updatecastling!(board::Board, sqr_from::Integer, sqr_to::Integer)
    board.hash ⊻= zobookey(board.castling)
    @inbounds board.castling &= CASTLING_RIGHT[sqr_from]
    @inbounds board.castling &= CASTLING_RIGHT[sqr_to]
    board.hash ⊻= zobookey(board.castling)
end

# Used for internally changing the castling rights when a move is played.
const CASTLING_RIGHT = @SVector [~0x01, ~0x00, ~0x00, ~0x05, ~0x00, ~0x00, ~0x00, ~0x04,
                                ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00,
                                ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00,
                                ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00,
                                ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00,
                                ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00,
                                ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00, ~0x00,
                                ~0x02, ~0x00, ~0x00, ~0x0a, ~0x00, ~0x00, ~0x00, ~0x08]

Move() = Move(zero(UInt16))

const MOVE_NONE = Move()
const NULL_MOVE = Move(0xffff)

"""
    Move(move_from::Integer, move_to::Integer, move_flag::Integer)

Encode a move, giving the from & to squares, alongside the move flag.
"""
function Move(move_from::Integer, move_to::Integer, move_flag::Integer)
    Move((move_from - one(move_from)) | ((move_to - one(move_from)) << 6) | move_flag)
end

"""
    Move(move_from::Integer, move_to::Integer, move_flag::Integer)

Encode a move, giving the from & to squares, assuming no move flags.
"""
function Move(move_from::Integer, move_to::Integer)
    Move((move_from - one(move_from)) | ((move_to - one(move_from)) << 6))
end
# the above function saves one bitwise 'or' operation per normal move generation.
# Move(move_from::Integer, move_to::Integer) = Move(move_from, move_to, __NORMAL_MOVE)

"""
    from(move::Move)

Given a move, retrieve the "move from" square, as an `Integer`.
"""
from(move::Move) = (move.val & 0x003f) + 0x0001

"""
    to(move::Move)

Given a move, retrieve the "move to" square, as an `Integer`.
"""
to(move::Move) = ((move.val >> 6) & 0x003f) + 0x0001

"""
    flag(move::Move)

Given a move, return any special flags, as an `Integer`.
"""
flag(move::Move) = move.val & 0xf000 #>> 12

"""
    istactical (board::Board, move::Move)

Given the board position and a move, return true if a move is tactical.
A tactical move involves a capture, promotion, or castling.
"""
function istactical(board::Board, move::Move)
    if flag(move) > __DOUBLE_PAWN
        # All flag values above __DOUBLE_PAWN are tactical.
        return true
    elseif board[to(move)] !== BLANK
        # Capture case
        return true
    else
        return false
    end
end

# Allows a preallocation for MoveStack
# MoveStack(size::Int) = MoveStack(Vector{Move}(undef, size), 0)
MoveStack(size::Int) = MoveStack(repeat([MOVE_NONE], size), 0)

# define useful array methods for MoveStack
Base.iterate(m::MoveStack, state = 1) = (state > m.idx) ? nothing : (m.list[state], state + 1)
Base.length(m::MoveStack) = m.idx
Base.eltype(::Type{MoveStack}) = Move
Base.size(m::MoveStack) = (m.idx, )
Base.IndexStyle(::Type{<:MoveStack}) = IndexLinear()
Base.getindex(m::MoveStack, idx::Int) = m.list[idx]
Base.setindex!(m::MoveStack, val::Move, idx::Int) = @inbounds m.list[idx] = val

# add moves to the MoveStack
function push!(m::MoveStack, move::Move)
    m.idx += 1
    @inbounds m.list[m.idx] = move
end

# pseudo-clear the MoveStack
clear!(m::MoveStack) = m.idx = 0

# Allows a preallocation for MoveStack
UndoStack(size::Int) = UndoStack(Vector{Undo}(undef, size), 0)

# define useful array methods for UndoStack
Base.iterate(u::UndoStack, state = 1) = (state > u.idx) ? nothing : (u.list[state], state + 1)
Base.length(u::UndoStack) = u.idx
Base.eltype(::Type{UndoStack}) = Undo
Base.size(u::UndoStack) = (u.idx, )
Base.IndexStyle(::Type{<:UndoStack}) = IndexLinear()
Base.getindex(u::UndoStack, idx::Int) = u.list[idx]

# add Undos to the UndoStack
function push!(u::UndoStack, undo::Undo)
    u.idx += 1
    @inbounds u.list[u.idx] = undo
end

# pseudo-clear the UndoStack
clear!(u::UndoStack) = u.idx = 0

"""
    apply_move!(board::Board, move::Move, undo::Undo)

Apply the given `move` to the `board`, adding changes to 'undo'.
"""
function apply_move!(board::Board, move::Move)
    undo_checkers = board.checkers
    undo_pinned = board.pinned
    undo_castling = board.castling
    undo_enpass = board.enpass
    undo_halfmovecount = board.halfmovecount
    undo_hash = board.hash
    undo_psqteval = board.psqteval
    undo_pkhash = board.pkhash

    # before we update the enpass
    if board.enpass !== zero(UInt8)
        board.hash ⊻= zobepkey(board.enpass)
    end

    # Apply the moves according to the appropriate flag
    if flag(move) <= __DOUBLE_PAWN
        # If the move is normal or a double pawn push
        undo_captured = apply_normal!(board, move)
    elseif flag(move) === __ENPASS
        # If the move is enpassant
        undo_captured = apply_enpass!(board, move)
    elseif flag(move) <= __QUEEN_CASTLE
        # If the move is castling fallthrough.
        undo_captured = apply_castle!(board, move)
    else
        # Otherwise it's a promotion fallthrough.
        undo_captured = apply_promo!(board, move)
    end

    # Finishing calculations, for the next turn
    board.hash ⊻= zobturnkey()
    switchturn!(board)

    board.checkers = kingAttackers(board)
    board.pinned = findpins(board)
    board.movecount += one(board.movecount)
    @inbounds board.history[board.movecount] = board.hash
    return Undo(undo_checkers, undo_pinned, undo_castling, undo_enpass, undo_captured, undo_halfmovecount, undo_hash, undo_psqteval, undo_pkhash)
end

# This function is the same as above
# ... but includes the option to update the thread movestack and piecestack for history heuristics.
function apply_move!(thread::Thread, move::Move)
    board = thread.board
    push!(thread.movestack, move)
    push!(thread.piecestack, type(board[from(move)]))
    u = apply_move!(board, move)
    return u
end

"""
    apply_normal!(board::Board, move::Move)

Apply the given `move` to the `board`. Assumes the move has either the `__NORMAL_MOVE` or `__DOUBLE_PAWN` flags.
"""
function apply_normal!(board::Board, move::Move)
    sqr_from = from(move)
    sqr_to = to(move)

    # Check for double pawn advance and set enpass square
    if flag(move) === __DOUBLE_PAWN
        if board.turn === WHITE
            board.enpass = UInt8(sqr_from + 8)
        else
            board.enpass = UInt8(sqr_from - 8)
        end
        board.hash ⊻= zobepkey(sqr_from)
    else
        board.enpass = zero(UInt8)
    end

    board.castling > 0x00 && updatecastling!(board, sqr_from, sqr_to)

    p_from = piece(board, sqr_from)
    p_to = piece(board, sqr_to)

    bb_from = Bitboard(sqr_from)
    bb_to = Bitboard(sqr_to)

    ptype_from = type(p_from)

    @inbounds board[ptype_from] ⊻= bb_from ⊻ bb_to
    @inbounds board[board.turn] ⊻= bb_from ⊻ bb_to
    @inbounds board[sqr_from] = BLANK
    @inbounds board[sqr_to] = p_from

    # Update PSQT
    board.psqteval -= psqt(p_from, sqr_from)
    board.psqteval += psqt(p_from, sqr_to)

    if p_to !== BLANK
        ptype_to = type(p_to)
        if ptype_to == PAWN
            board.pkhash ⊻= zobkey(p_to, sqr_to)
        end
        @inbounds board[ptype_to] ⊻= bb_to
        @inbounds board[!board.turn] ⊻= bb_to
        board.hash ⊻= zobkey(p_to, sqr_to)

        # Update PSQT
        board.psqteval -= psqt(p_to, sqr_to)
    end

    if (ptype_from === PAWN) || (p_to !== BLANK)
        board.halfmovecount = 0
    else
        board.halfmovecount += 1
    end

    board.hash ⊻= zobkey(p_from, sqr_from)
    board.hash ⊻= zobkey(p_from, sqr_to)

    if ptype_from === PAWN || ptype_from == KING
        board.pkhash ⊻= zobkey(p_from, sqr_from)
        board.pkhash ⊻= zobkey(p_from, sqr_to)
    end

    return p_to
end

"""
    apply_enpass!(board::Board, move::Move)

Apply the given `move` to the `board`. Assumes the move has the `__ENPASS` flag.
"""
function apply_enpass!(board::Board, move::Move)
    board.enpass = zero(UInt8)

    sqr_from = from(move)
    sqr_to = to(move)

    cap_sqr = sqr_to - 24 + (board.turn.val << 4)
    cap_bb = Bitboard(cap_sqr)

    bb_from = Bitboard(sqr_from)
    bb_to = Bitboard(sqr_to)

    @inbounds board[PAWN] ⊻= bb_from ⊻ bb_to ⊻ cap_bb
    @inbounds board[board.turn] ⊻= bb_from ⊻ bb_to
    @inbounds board[!board.turn] ⊻= cap_bb

    # First set the square to
    @inbounds board[sqr_to] = piece(board, sqr_from)
    # Then clear the square from
    @inbounds board[sqr_from] = BLANK
    # Then clear the captured square
    @inbounds board[cap_sqr] = BLANK

    p_from = makepiece(PAWN, board.turn)
    p_capt = makepiece(PAWN, !board.turn)
    board.hash ⊻= zobkey(p_from, sqr_from)
    board.hash ⊻= zobkey(p_from, sqr_to)
    board.hash ⊻= zobkey(p_capt, cap_sqr)

    # Update PSQT
    board.psqteval -= psqt(p_from, sqr_from)
    board.psqteval += psqt(p_from, sqr_to)
    board.psqteval -= psqt(p_capt, cap_sqr)

    board.halfmovecount = 0

    board.pkhash ⊻= zobkey(p_capt, cap_sqr)
    board.pkhash ⊻= zobkey(p_from, sqr_from)
    board.pkhash ⊻= zobkey(p_from, sqr_to)

    return p_capt
end

"""
    apply_castle!(board::Board, move::Move)

Apply the given `move` to the `board`. Assumes the move has the `__KING_CASTLE` or `__QUEEN_CASTLE` flag.
"""
function apply_castle!(board::Board, move::Move)
    board.enpass = zero(UInt8)

    k_from = from(move)
    k_to = to(move)

    updatecastling!(board, k_from, k_to)

    if flag(move) === __KING_CASTLE
        r_from = k_from - 3
        r_to = k_from - 1
    else
        r_from = k_from + 4
        r_to = k_from + 1
    end

    r_from_bb = Bitboard(r_from)
    r_to_bb = Bitboard(r_to)
    k_from_bb = Bitboard(k_from)
    k_to_bb = Bitboard(k_to)

    @inbounds board[KING] ⊻= k_from_bb ⊻ k_to_bb
    @inbounds board[ROOK] ⊻= r_from_bb ⊻ r_to_bb

    @inbounds board[board.turn] ⊻= k_from_bb ⊻ k_to_bb ⊻ r_from_bb ⊻ r_to_bb

    @inbounds board[k_from] = BLANK
    @inbounds board[r_from] = BLANK

    _king = makepiece(KING, board.turn)
    _rook = makepiece(ROOK, board.turn)
    @inbounds board[k_to] = _king
    @inbounds board[r_to] = _rook

    board.hash ⊻= zobkey(_king, k_from)
    board.hash ⊻= zobkey(_king, k_to)
    board.hash ⊻= zobkey(_rook, r_from)
    board.hash ⊻= zobkey(_rook, r_to)

    # Update PSQT
    board.psqteval -= psqt(_king, k_from)
    board.psqteval -= psqt(_rook, r_from)
    board.psqteval += psqt(_king, k_to)
    board.psqteval += psqt(_rook, r_to)

    board.pkhash ⊻= zobkey(_king, k_from)
    board.pkhash ⊻= zobkey(_king, k_to)

    board.halfmovecount += 1

    return BLANK
end

"""
    apply_promo!(board::Board, move::Move)

Apply the given `move` to the `board`. Assumes the move has the `__<PIECE>_PROMO` flag, where `<PIECE>` is either a `KNIGHT`, `BISHOP`, `ROOK`, or `QUEEN`.
"""
function apply_promo!(board::Board, move::Move)
    board.enpass = zero(UInt8)

    sqr_from = from(move)
    sqr_to = to(move)

    board.castling > 0x00 && updatecastling!(board, sqr_from, sqr_to)

    bb_from = Bitboard(sqr_from)
    bb_to = Bitboard(sqr_to)

    p_to = piece(board, sqr_to)
    ptype_promo = PieceType((flag(move)>>12) - 3)

    @inbounds board[PAWN] ⊻= bb_from
    @inbounds board[ptype_promo] ⊻= bb_to
    @inbounds board[board.turn] ⊻= bb_from ⊻ bb_to

    @inbounds board[sqr_from] = BLANK

    p_promo = makepiece(ptype_promo, board.turn)
    @inbounds board[sqr_to] = p_promo

    # Update PSQT
    board.psqteval -= psqt(makepiece(PAWN, board.turn), sqr_from)
    board.psqteval += psqt(p_promo, sqr_to)

    if p_to !== BLANK
        @inbounds board[type(p_to)] ⊻= bb_to
        @inbounds board[!board.turn] ⊻= bb_to
        board.hash ⊻= zobkey(p_to, sqr_to)

        # Update PSQT
        board.psqteval -= psqt(p_to, sqr_to)
    end

    p_from = makepiece(PAWN, board.turn)
    board.hash ⊻= zobkey(p_from, sqr_from)
    board.hash ⊻= zobkey(p_promo, sqr_to)

    board.pkhash ⊻= zobkey(p_from, sqr_from)

    board.halfmovecount = 0

    return p_to
end

function undo_move!(thread::Thread, move::Move, undo::Undo)
    board = thread.board
    thread.movestack.idx -= 1
    thread.piecestack.idx -= 1
    undo_move!(board, move, undo)
    return
end

function undo_move!(board::Board, move::Move, undo::Undo)
    board.checkers = undo.checkers
    board.pinned = undo.pinned
    board.enpass = undo.enpass
    board.castling = undo.castling
    board.halfmovecount = undo.halfmovecount
    board.movecount -= one(UInt16)
    board.hash = undo.hash
    board.psqteval = undo.psqteval
    board.pkhash = undo.pkhash
    switchturn!(board)
    if flag(move) <= __DOUBLE_PAWN
        # If the move was a normal move or double pawn push
        undo_normal!(board, move, undo)
    elseif flag(move) === __ENPASS
        # If the move was enpassant
        undo_enpass!(board, move, undo)
    elseif flag(move) <= __QUEEN_CASTLE
        # If the move was a castling move, fallthrough.
        undo_castle!(board, move, undo)
    else
        # Otherwise, it was a promotion fallthrough.
        undo_promo!(board, move, undo)
    end
    return
end

function undo_normal!(board::Board, move::Move, undo::Undo)
    sqr_from = to(move)
    sqr_to = from(move)

    p_from = piece(board, sqr_from)
    p_to = undo.captured #piece(board, sqr_to)

    bb_from = Bitboard(sqr_from)
    bb_to = Bitboard(sqr_to)

    @inbounds board[type(p_from)] ⊻= bb_from ⊻ bb_to
    @inbounds board[board.turn] ⊻= bb_from ⊻ bb_to
    @inbounds board[sqr_from] = p_to
    @inbounds board[sqr_to] = p_from

    if p_to !== BLANK
        @inbounds board[type(p_to)] ⊻= bb_from
        @inbounds board[!board.turn] ⊻= bb_from
    end
    return
end

function undo_enpass!(board::Board, move::Move, undo::Undo)
    sqr_from = from(move)
    sqr_to = to(move)

    cap_sqr = sqr_to - 24 + (board.turn.val << 4)
    cap_bb = Bitboard(cap_sqr)

    bb_from = Bitboard(sqr_from)
    bb_to = Bitboard(sqr_to)

    @inbounds board[PAWN] ⊻= bb_from ⊻ bb_to ⊻ cap_bb
    @inbounds board[board.turn] ⊻= bb_from ⊻ bb_to
    @inbounds board[!board.turn] ⊻= cap_bb

    @inbounds board[sqr_from] = piece(board, sqr_to)
    @inbounds board[sqr_to] = BLANK
    @inbounds board[cap_sqr] = undo.captured
    return
end

function undo_castle!(board::Board, move::Move, undo::Undo)
    k_from = from(move)
    k_to = to(move)

    if flag(move) === __KING_CASTLE
        r_from = k_from - 3
        r_to = k_from - 1
    else
        r_from = k_from + 4
        r_to = k_from + 1
    end

    r_from_bb = Bitboard(r_from)
    r_to_bb = Bitboard(r_to)
    k_from_bb = Bitboard(k_from)
    k_to_bb = Bitboard(k_to)

    @inbounds board[KING] ⊻= k_from_bb ⊻ k_to_bb
    @inbounds board[ROOK] ⊻= r_from_bb ⊻ r_to_bb

    @inbounds board[board.turn] ⊻= k_from_bb ⊻ k_to_bb ⊻ r_from_bb ⊻ r_to_bb

    @inbounds board[k_from] = makepiece(KING, board.turn)
    @inbounds board[r_from] = makepiece(ROOK, board.turn)
    @inbounds board[k_to] = BLANK
    @inbounds board[r_to] = BLANK
    return
end

function undo_promo!(board::Board, move::Move, undo::Undo)
    sqr_from = from(move)
    sqr_to = to(move)

    bb_from = Bitboard(sqr_from)
    bb_to = Bitboard(sqr_to)

    p_to = undo.captured
    ptype_promo = PieceType((flag(move) >> 12) - 3)

    @inbounds board[PAWN] ⊻= bb_from
    @inbounds board[ptype_promo] ⊻= bb_to
    @inbounds board[board.turn] ⊻= bb_from ⊻ bb_to

    @inbounds board[sqr_from] = makepiece(PAWN, board.turn)
    @inbounds board[sqr_to] = p_to

    if p_to !== BLANK
        @inbounds board[type(p_to)] ⊻= bb_to
        @inbounds board[!board.turn] ⊻= bb_to
    end
    return
end

# "pass" the go, and let out opponent have another move.
function apply_null!(thread::Thread)
    board = thread.board
    push!(thread.movestack, NULL_MOVE)
    push!(thread.piecestack, VOID)
    undo_checkers = board.checkers
    undo_pinned = board.pinned
    undo_castling = board.castling
    undo_enpass = board.enpass
    undo_halfmovecount = board.halfmovecount
    undo_hash = board.hash
    undo_psqteval = board.psqteval
    undo_pkhash = board.pkhash

    # If the position had an enpassant square, set the key as needed, and turn off the enpass square flag.
    if board.enpass !== zero(UInt8)
        board.hash ⊻= zobepkey(board.enpass)
        board.enpass = zero(UInt8)
    end

    # Finishing calculations, for the next turn
    board.hash ⊻= zobturnkey()
    switchturn!(board)

    #board.checkers = kingAttackers(board)
    board.pinned = findpins(board)
    board.movecount += one(board.movecount)
    @inbounds board.history[board.movecount] = board.hash
    return Undo(undo_checkers, undo_pinned, undo_castling, undo_enpass, BLANK, undo_halfmovecount, undo_hash, undo_psqteval, undo_pkhash)
end

function undo_null!(thread::Thread, undo::Undo)
    board = thread.board
    thread.movestack.idx -= 1
    thread.piecestack.idx -= 1
    board.checkers = undo.checkers
    board.pinned = undo.pinned
    board.enpass = undo.enpass
    board.castling = undo.castling
    board.halfmovecount = undo.halfmovecount
    board.movecount -= one(UInt16)
    board.hash = undo.hash
    board.psqteval = undo.psqteval
    board.pkhash = undo.pkhash
    switchturn!(board)
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
macro return_exception(ex)
    quote
        ans = $(esc(ex))
        ans isa Exception && return ans
        ans
    end
end

macro require_keyword_arguments(input)
    quote
        if !($(esc(input)) isa AbstractDict)
            return ArgumentError("Require keyword argument.")
        end
    end
end

macro ensure_exception(ex)
    quote
        err = $(esc(ex))
        if err isa Exception
            err
        else
            @error "Caught a non-exception type." exception = (err, catch_backtrace())
            ErrorException("Non-exception type `$(typeof(err))` is thrown.")
        end
    end
end

function try_impl(ex)
    quote
        try
            $(esc(ex))
        catch err
            return @ensure_exception err
        end
    end
end

@eval macro $(:try)(ex)
    try_impl(ex)
end

refine_values(x) = x
refine_values(x::AbstractVector) = refine_values.(x)
refine_values(x::AbstractDict) = Dict(zip(keys(x), refine_values.(values(x))))

function getstrings(dict, key)
    try
        return append!(String[], get(dict, key, []))
    catch err
        @debug "Converting `$key` failed." exception = (err, catch_backtrace())
        return ArgumentError("Keyword argument `$key` must be of a vector of strings.")
    end
end

const ReviseError = -3857439168074868810


text:
"""
VshTrack(args...)  -> PyPlot Plot

Return a plot of Vshale logs and optionally Effective porosity in a single track.

 The next table show the list of variables allowed:

 The next table show the list of variables allowed:
|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|Depth|Mandatory| -- | Array{Number,1}| 1D Array of Depths|
|Vsh|Mandatory| -- | Array{Number,1}|1D Array of Vshale Values|
|PhiePlot|Optional| -- | Array{Number,1}|1D Array of Phie Values|
|DepthFrom|Optional| minimum(Depth) | Number |Set the Depth upper limit to plot|
|DepthTo|Optional| maximum(Depth) | Number |  Set the Depth lower limit to plot|
|Dtick|Optional|true|Dtick=Bool|Set the Depth ticks On/Off|
"""
function VshPhieTrack(Depth,Vsh;DepthFrom=minimum(Depth),DepthTo=maximum(Depth),Dtick=false,PhiePlot=false,Phiemax=0.30)
vshax=gca()
vshax.plot(Vsh,Depth,linewidth=1,linestyle="--",color="black")
vshax.set_xlim([0,1])
vshax.set_ylim([DepthTo,DepthFrom])
vshax.set_yticks(range(DepthFrom, stop=DepthTo, length=51),minor=true)
vshax.set_yticks(range(DepthFrom, stop=DepthTo, length=11))
vshax.tick_params("both",labelsize=6)
vshax.set_xlabel("VShale")
setp(vshax.get_yticklabels(),visible=Dtick)
vshax.fill_betweenx(Depth,Vsh,1,where=(Vsh .<= 1),color="gold")

if PhiePlot!=false
phieax=vshax.twiny()
phieax.plot(PhiePlot,Depth,linewidth=1,color="black")
phieax.set_xlim([Phiemax,0])
setp(phieax.get_yticklabels(),visible=Dtick)
phieax.set_xlabel("Phie")
phieax.tick_params("both",labelsize=6)

phieax.fill_betweenx(Depth,PhiePlot,0,where=(Vsh .<= 1),color="red")
else
    vshax.xaxis.set_label_position("top")
    vshax.xaxis.tick_top()

end

end
###########################################################
#VshTrack Function
###########################################################
function VshTrack(Depth,Vsh;DepthFrom=minimum(Depth),DepthTo=maximum(Depth),Dtick=false,CutOff=0.7,fill=true)
vshax=gca()
vshax.plot(Vsh,Depth,linewidth=1,linestyle="--",color="black")
vshax.set_xlim([1,0])
vshax.set_ylim([DepthTo,DepthFrom])
vshax.set_yticks(range(DepthFrom, stop=DepthTo, length=51),minor=true)
vshax.set_yticks(range(DepthFrom, stop=DepthTo, length=11))
vshax.tick_params("both",labelsize=6)
vshax.set_xlabel("VShale")
setp(vshax.get_yticklabels(),visible=Dtick)
if fill==true
    vshax.fill_betweenx(Depth,CutOff,Vsh,where=(CutOff .>= Vsh),color="gold")
end
vshax.xaxis.set_label_position("top")
vshax.xaxis.tick_top()
end

function PhieTrack(Depth,Phie;DepthFrom=minimum(Depth),DepthTo=maximum(Depth),Dtick=false,CutOff=0.08,Phiemax=0.30,fill=false)
    phieax=gca()
    phieax.plot(Phie,Depth,linewidth=1,color="black")
    phieax.set_ylim([DepthTo,DepthFrom])
    phieax.set_xlim([0,Phiemax])
    setp(phieax.get_yticklabels(),visible=Dtick)
    phieax.set_xlabel("Phie")
    phieax.tick_params("both",labelsize=6)
    if fill==true
        phieax.fill_betweenx(Depth,CutOff,Phie,where=(CutOff .<= Phie),color="red")
    end
    phieax.xaxis.set_label_position("top")
    phieax.xaxis.tick_top()
end


text:
@info "Caesar is loading tools related to AprilTags.jl."

using .AprilTags

include("Pose2AprilTag4Corners.jl")

#


text:
wiener_randn!() = nothing

function __init__()
    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
        wiener_randn!(rng::AbstractRNG,rand_vec::CuArrays.CuArray) = randn!(rand_vec)
    end
end


text:
P = 23800
b = P * 0.05 + 300
i = 0.1
N = 1

F = P * (1 + i)^N
ic = (F / (P - b))^(1 / N) - 1


text:
include("held_suarez_resolution.jl")
##
field_function_list = (
    (; comparison_field = "u", func = grab_state),
    (; comparison_field = "T", func = grab_state),
    (; comparison_field = "TT", func = eddy_variance),
    (; comparison_field = "uv", func = eddy_variance),
    (; comparison_field = "vT", func = eddy_variance),
    (; comparison_field = "uu", func = eddy_variance),
    (; comparison_field = "v", func = grab_state),
    (; comparison_field = "vv", func = eddy_variance),
    (; comparison_field = "wT", func = eddy_variance),
    (; comparison_field = "ρ", func = grab_state),
    (; comparison_field = "ρe", func = grab_state),
    (; comparison_field = "w", func = grab_state),
)

field_function_list = field_function_list[1:6]
error_matrix = zeros(5, length(field_function_list))
for (i, field_function) in enumerate(field_function_list)
    comparison_field = field_function.comparison_field
    grab = field_function.func

    jl_file = jldopen(filename[1], "r+")
    u_truth = grab(comparison_field, jl_file)
    u_truth isa Tuple ? u_truth = u_truth[1] : nothing

    jl_file = jldopen(filename[2], "r+")
    u_2 = grab(comparison_field, jl_file)
    u_2 isa Tuple ? u_2 = u_2[1] : nothing

    jl_file = jldopen(filename[3], "r+")
    u_3 = grab(comparison_field, jl_file)
    u_3 isa Tuple ? u_3 = u_3[1] : nothing

    jl_file = jldopen(filename[4], "r+")
    u_4 = grab(comparison_field, jl_file)
    u_4 isa Tuple ? u_4 = u_4[1] : nothing

    jl_file = jldopen(filename[5], "r+")
    u_5 = grab(comparison_field, jl_file)
    u_5 isa Tuple ? u_5 = u_5[1] : nothing

    jl_file = jldopen(filename[6], "r+")
    u_6 = grab(comparison_field, jl_file)
    u_6 isa Tuple ? u_6 = u_6[1] : nothing

    error_list = Float64[]
    er1 = abs.(u_truth - u_2) / maximum(abs.(u_truth))
    er2 = abs.(u_truth - u_3) / maximum(abs.(u_truth))
    er3 = abs.(u_truth - u_4) / maximum(abs.(u_truth))
    er4 = abs.(u_truth - u_5) / maximum(abs.(u_truth))
    er5 = abs.(u_truth - u_6) / maximum(abs.(u_truth))
    qn = 1.0 # quantile number
    push!(error_list, quantile(er1[:], qn))
    push!(error_list, quantile(er2[:], qn))
    push!(error_list, quantile(er3[:], qn))
    push!(error_list, quantile(er4[:], qn))
    push!(error_list, quantile(er5[:], qn))
    # push!(error_list, maximum(abs.(u_truth - u_7)))

    error_matrix[:, i] .= error_list

    minerror = title_names[argmin(error_list)+1] # since u_truth is 1
    println("The minimum error for " * comparison_field, " is ", minerror)
    println("All the errors are ", error_list)
    println("-------------------------")

end
# defined in utils.jl
##
include("utils.jl")
latex_format(error_matrix)

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
using Flo: julia_main

julia_main()


text:
module ProtoBufTestApis
using ProtoBuf
using ..Test
import ProtoBuf.meta

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

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

function test_apis()
    t = TestType()

    @test [:a, :b] == propertynames(t)
    @test !hasproperty(t, :a)
    @test !hasproperty(t, :b)

    @test_throws KeyError getproperty(t, :a)

    t.b = true
    @test hasproperty(t, :b)
    @test (getproperty(t, :b) == true)

    @test !isinitialized(t)
    t.a = "hello world"
    @test isinitialized(t)
    @test (t.a ==  "hello world")

    clear(t, :b)
    @test isinitialized(t)
    clear(t)
    @test !isinitialized(t)

    t = TestType(; a="hello", b=false)
    @test t.a == "hello"
    @test t.b == false
end

function test_deepcopy()
    ts = ProtoBuf.google.protobuf.Timestamp()
    ts.seconds = 123
    @test !hasproperty(ts, :nanos)
    ts2 = deepcopy(ts)
    @test !hasproperty(ts2, :nanos)
    @test hasproperty(ts2, :seconds)
end

end # module ProtoBufTestApis

@testset "Utility APIs" begin
    ProtoBufTestApis.test_apis()
    ProtoBufTestApis.test_deepcopy()
end


text:
module BB
using AA

export likes

likes = AA.something * " banana cards"

end


text:
intlog2(x::T) where {T<:Integer} =  T(8*sizeof(T) - leading_zeros(x-one(T)))
intlog2(x) = x > typemax(UInt) ? ceil(Int, log2(x)) : intlog2(ceil(UInt,x)) % Int

naivemul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) where {T<:LinearAlgebra.BlasFloat} = mul!(C,A,B)
function naivemul!(C, A, B)
    Maxis, Naxis = axes(C)
    # TODO: heuristically pick `Nunroll` and `Munroll` using `sizeof(T)`, and maybe based on size of register file as well.
    # Nunroll = 2
    # Munroll = 4
    # For now, the `4` and `2` are hardcoded to use `Base.Cartesian.@nexprs` without generated functions.
    # But a minor code reorgaanization, e.g. loading and storing tuples while moving the `Base.Cartesian.@nexprs`
    # into `tload` and `tstore!` functions could let us switch APIs.
    nstep = step(Naxis)
    mstep = step(Maxis)
    # I don't want to deal with axes having non-unit step
    if nstep == mstep == 1
        if sizeof(eltype(C)) > 256
            naivemul!(C, A, B, Maxis, Naxis, Val(1), Val(1))
        elseif sizeof(eltype(C)) > 128
            naivemul!(C, A, B, Maxis, Naxis, Val(2), Val(1))
        elseif sizeof(eltype(C)) > 96
            naivemul!(C, A, B, Maxis, Naxis, Val(4), Val(1))
        elseif sizeof(eltype(C)) > 64
            naivemul!(C, A, B, Maxis, Naxis, Val(4), Val(2))
        else
            naivemul!(C, A, B, Maxis, Naxis, Val(4), Val(3))
        end
    else
        mul!(C,A,B)
    end
end
_const(A) = A
_const(A::Array) = Base.Experimental.Const(A)
# Separated to make it easier to test.
@generated function naivemul!(C::AbstractMatrix{T}, A, B, Maxis, Naxis, ::Val{MU}, ::Val{NU}) where {T,MU,NU}
    nrem_body = quote
        m = first(Maxis)-1
        while m < M - $(MU-1)
            Base.Cartesian.@nexprs $MU i -> Cmn_i = zero(T)
            for k ∈ Kaxis
                Base.Cartesian.@nexprs $MU i -> Cmn_i = muladd(_const(A)[m+i,k],_const(B)[k,nn],Cmn_i)
            end
            Base.Cartesian.@nexprs $MU i -> C[m+i,nn] = Cmn_i
            m += $MU
        end
        for mm ∈ 1+m:M
            Cmn = zero(T)
            for k ∈ Kaxis
                Cmn = muladd(_const(A)[mm,k], _const(B)[k,nn], Cmn)
            end
            C[mm,nn] = Cmn
        end
    end
    nrem_quote = if NU > 2
        :(for nn ∈ 1+n:N; $nrem_body; end)
    else
        :(let nn = N; $nrem_body; end)
    end
    quote
        N = last(Naxis)
        M = last(Maxis)
        Kaxis = axes(B,1)
        Base.Experimental.@aliasscope begin
            n = first(Naxis)-1
            @inbounds begin
                while n < N - $(NU-1)
                    m = first(Maxis)-1
                    while m < M - $(MU-1)
                        Base.Cartesian.@nexprs $NU j -> Base.Cartesian.@nexprs $MU i -> Cmn_i_j = zero(T)
                        for k ∈ Kaxis
                            Base.Cartesian.@nexprs $MU i -> Ak_i = _const(A)[m+i,k]
                            Base.Cartesian.@nexprs $NU j -> begin
                                Bk_j = _const(B)[k,n+j]
                                Base.Cartesian.@nexprs $MU i -> Cmn_i_j = muladd(Ak_i, Bk_j, Cmn_i_j)
                            end
                        end
                        Base.Cartesian.@nexprs $NU j -> Base.Cartesian.@nexprs $MU i -> C[m+i,n+j] = Cmn_i_j
                        m += $MU
                    end
                    for mm ∈ 1+m:M
                        Base.Cartesian.@nexprs $NU j -> Cmn_j = zero(T)
                        for k ∈ Kaxis
                            Base.Cartesian.@nexprs $NU j -> Cmn_j = muladd(_const(A)[mm,k],_const(B)[k,n+j],Cmn_j)
                        end
                        Base.Cartesian.@nexprs $NU j -> C[mm,n+j] = Cmn_j
                    end
                    n += $NU
                end
                $(NU > 1 ? nrem_quote : nothing)
            end
        end
        C
    end
end

"""
    struct ExpMethodGeneric{T}
    ExpMethodGeneric()=ExpMethodGeneric{Val{13}}();

Generic exponential implementation of the method `ExpMethodHigham2005`,
for any exp argument `x`  for which the functions
`LinearAlgebra.opnorm`, `+`, `*`, `^`, and `/` (including addition with
UniformScaling objects) are defined. The type `T` is used to adjust the
number of terms used in the Pade approximants at compile time.

See "The Scaling and Squaring Method for the Matrix Exponential Revisited"
by Higham, Nicholas J. in 2005 for algorithm details.

"""
struct ExpMethodGeneric{T}
end
ExpMethodGeneric()=ExpMethodGeneric{Val{13}}();

function _exp!(x,method::ExpMethodGeneric{Vk},cache=alloc_mem(A,method))}()) where Vk
    nx = opnorm(x, 1)
    if !isfinite(nx)
        # This should (hopefully) ensure that the result is Inf or NaN depending on
        # which values are produced by a power series. I.e. the result shouldn't
        # be all Infs when there are both Infs and zero since Inf*0===NaN
        return x*sum(x*nx)
    end
    s = iszero(nx) ? 0 : intlog2(nx)
    (Vk === Val{13}() && x isa AbstractMatrix && ismutable(x)) && return exp_generic_mutable(x, s, Val{13}())
    if s >= 1
        exp_generic(x/(2^s), Vk)^(2^s)
    else
        exp_pade_p(x, Vk, Vk) / exp_pade_q(x, Vk, Vk)
    end
end

function exp_pade_p(x, ::Val{13}, ::Val{13})
    @evalpoly(x,LinearAlgebra.UniformScaling{Float64}(1.0),
                LinearAlgebra.UniformScaling{Float64}(0.5),
                LinearAlgebra.UniformScaling{Float64}(0.12),
                LinearAlgebra.UniformScaling{Float64}(0.018333333333333333),
                LinearAlgebra.UniformScaling{Float64}(0.0019927536231884057),
                LinearAlgebra.UniformScaling{Float64}(0.00016304347826086958),
                LinearAlgebra.UniformScaling{Float64}(1.0351966873706003e-5),
                LinearAlgebra.UniformScaling{Float64}(5.175983436853002e-7),
                LinearAlgebra.UniformScaling{Float64}(2.0431513566525008e-8),
                LinearAlgebra.UniformScaling{Float64}(6.306022705717595e-10),
                LinearAlgebra.UniformScaling{Float64}(1.48377004840414e-11),
                LinearAlgebra.UniformScaling{Float64}(2.529153491597966e-13),
                LinearAlgebra.UniformScaling{Float64}(2.8101705462199623e-15),
                LinearAlgebra.UniformScaling{Float64}(1.5440497506703088e-17))
end

function exp_generic_mutable(x::AbstractMatrix{T}, s, ::Val{13}) where {T}
    y1 = similar(x, promote_type(T,Float64))
    y2 = similar(y1)
    y3 = similar(y1)
    exp_generic!(y1, y2, y3, x, s, Val{13}())
end
function exp_generic!(y1, y2, y3, x, s, ::Val{13})
    if s > 0
        _y3 = exp_generic_core!(y1, y2, y3, x .* (1/(1 << s)), Val{13}())
        if typeof(_y3) === (y3)
            y3 = _y3
        else
            y3 .= _y3
        end
        for _ ∈ 1:s
            naivemul!(y1, y3, y3)
            y3, y1 = y1, y3
        end
        return y3
    else
        return exp_generic_core!(y1, y2, y3, x, Val{13}())
    end
end
# `lu!` is only defined for `StridedMatrix`, and `lu(::StaticArray)` (note `MArray<:StaticArray`) returns `::StaticArrays.LU !== LinearAlgebra.LU`
_rdiv!(A, B::StridedMatrix) = rdiv!(A, lu!(B))
_rdiv!(A, B) = A / B

function exp_generic_core!(y1, y2, y3, x, ::Val{13})
    @inbounds for i ∈ eachindex(y3)
        y3[i] = -x[i]
    end
    exp_pade_p!(y1, y2, y3, Val{13}(), Val{13}())
    exp_pade_p!(y3, y2, x, Val{13}(), Val{13}())
    return _rdiv!(y3, y1)
end

function exp_pade_p!(y1::AbstractMatrix{T}, y2::AbstractMatrix{T}, x::AbstractMatrix, ::Val{13}, ::Val{13}) where {T}
    N = size(x,1) # check square is in `exp_generic`
    y1 .= x .* 1.5440497506703088e-17
    for c ∈ (
        2.8101705462199623e-15, 2.529153491597966e-13, 1.48377004840414e-11, 6.306022705717595e-10,
        2.0431513566525008e-8, 5.175983436853002e-7, 1.0351966873706003e-5, 0.00016304347826086958,
        0.0019927536231884057, 0.018333333333333333, 0.12, 0.5)

        @inbounds for n ∈ 1:N
            y1[n,n] += c
        end
        naivemul!(y2, y1, x)
        y1, y2 = y2, y1
    end
    @inbounds for n ∈ 1:N
        y1[n,n] += one(T)
    end
    return y1
end

function exp_pade_p(x::Number, ::Val{13}, ::Val{13})
    @evalpoly(x,1.0,0.5,0.12,0.018333333333333333,0.0019927536231884057,
              0.00016304347826086958,1.0351966873706003e-5,5.175983436853002e-7,
              2.0431513566525008e-8,6.306022705717595e-10,1.48377004840414e-11,
              2.529153491597966e-13,2.8101705462199623e-15,1.5440497506703088e-17)
end

@generated function exp_pade_p(x, ::Val{k}, ::Val{m}) where {k, m}
    factorial = Base.factorial ∘ big
    p = map(Tuple(0:k)) do j
        num = factorial(k + m - j) * factorial(k)
        den = factorial(k + m) * factorial(k - j)*factorial(j)
        (float ∘ eltype)(x)(num // den) * (x <: Number ? 1 : I)
    end
    :(@evalpoly(x, $(p...)))
end

exp_pade_q(x, k, m) = exp_pade_p(-x, m, k)


text:
module TestIndexers

using Base.Test

using Bio.Indexers

@testset "Indexers" begin

    symbolnames = [:First, :Second, :Third, :Fourth, :Fifth]
    textnames = ["First", "Second", "Third", "Fourth", "Fifth"]
    groups = UnitRange{Int}[1:5, 6:10, 11:15, 16:20, 21:25]
    groups_eight = UnitRange{UInt8}[1:5, 6:10, 11:15, 16:20, 21:25]

    int_idxes = (Indexer(symbolnames),
                 Indexer(textnames, UInt64),
                 Indexer(symbolnames, UInt8))

    g_idxes = (Indexer(symbolnames, groups),
               Indexer(textnames, groups),
               Indexer(textnames, groups_eight))

    empty_index = Indexer(UInt64)

    @testset "Constructor Behaviour" for idxset in (int_idxes, g_idxes)
        for idxer in idxset
            @test idxer.names !== symbolnames && idxer.names == symbolnames
        end
    end

    @testset "Index equality" begin
        @testset "Single Indexes" for idxer in int_idxes
            for idx in int_idxes
                @test idxer == idx
            end
            for idx in g_idxes
                @test idxer != g_idxes
            end
            @test idxer != empty_index
        end
        @testset "Group Indexes" for idxer in g_idxes
            for idx in g_idxes
                @test idxer == idx
            end
            for idx in int_idxes
                @test idxer != idx
            end
            @test idxer != empty_index
        end
    end

    @testset "Basic Operators" begin
        # Tests for single indexes and group indexes.
        for idxset in (int_idxes, g_idxes)
            for idxer in idxset
                @test length(idxer) == 5
                @test names(idxer) == symbolnames
                @test names(idxer) !== symbolnames
                @test names(idxer) == idxer.names
                @test names(idxer) !== idxer.names
                for i in 1:5
                    @test haskey(idxer, symbolnames[i]) == true
                    @test haskey(idxer, textnames[i]) == true
                    @test haskey(idxer, i) == true
                end
                @test keys(idxer) == symbolnames
                @test keys(idxer) !== symbolnames
                idxcopy = copy(idxer)
                @test isequal(idxcopy, idxer)
                @test is(idxcopy, idxer) == false
                @test idxcopy.names == idxer.names
                @test idxcopy.names !== idxer.names
                @test idxcopy.lookup == idxer.lookup
                @test idxcopy.lookup !== idxer.lookup
            end
        end
        # Tests for empty index seperate because they are a bit different.
        @test names(empty_index) == Symbol[]
        @test names(empty_index) == Symbol[]
        @test names(empty_index) !== Symbol[]
        @test names(empty_index) !== empty_index.names
        @test names(empty_index) == empty_index.names
        for i in 1:5
            @test haskey(empty_index, symbolnames[i]) == false
            @test haskey(empty_index, textnames[i]) == false
            @test haskey(empty_index, i) == false
        end
        @test keys(empty_index) == Symbol[]
        @test keys(empty_index) !== Symbol[]
        idxcopy = copy(empty_index)
        @test isequal(idxcopy, empty_index)
        @test is(idxcopy, empty_index) == false
        @test idxcopy.names == empty_index.names
        @test idxcopy.names !== empty_index.names
        @test idxcopy.lookup == empty_index.lookup
        @test idxcopy.lookup !== empty_index.lookup
    end

    @testset "Manipulation methods" begin
        @testset "Setting and resetting names" begin

            # Ok let's give int_index_one a load of new names
            newnames = [:Sixth, :Seventh, :Eighth, :Ninth, :Tenth]
            newnamestext = ["Sixth", "Seventh", "Eighth", "Ninth", "Tenth"]
            @test names!(int_idxes[1], newnames) == Indexer(newnames)
            @test_throws ArgumentError names!(int_idxes[1], newnames[1:4])
            @test names!(int_idxes[2], newnamestext) == Indexer(newnamestext)
            @test_throws ArgumentError names!(int_idxes[2], newnamestext[1:4])

            # Now let's manipulate some of those names by renaming them

            # Make a new index by replacing name :Sixth with :First,
            # and test that only :Sixth is different
            rn = rename(int_idxes[1], :Sixth, :First)
            @test (rn != Indexer(symbolnames)) && (rn != Indexer(newnames))
            @test names(rn)[1] == :First
            @test names(rn)[2] == :Seventh
            @test names(rn)[3] == :Eighth
            @test names(rn)[4] == :Ninth
            @test names(rn)[5] == :Tenth

            # Let's modify int_index_one by replacing the names, Sixth Eights, and
            # Tenth, with First, Third, and Fifth.
            rename!(int_idxes[1], [:Sixth, :Eighth, :Tenth], [:First, :Third, :Fifth])
            @test int_idxes[1] != Indexer(symbolnames)
            @test rn != int_idxes[1]
            @test rn !== int_idxes[1]

            # Now let's test that renaming worked correctly, some names should
            # match names in symbolnames, and some should match names in
            # newnames.
            @test names(int_idxes[1])[1] != newnames[1]
            @test names(int_idxes[1])[1] == symbolnames[1]
            @test names(int_idxes[1])[2] == newnames[2]
            @test names(int_idxes[1])[2] != symbolnames[2]
            @test names(int_idxes[1])[3] != newnames[3]
            @test names(int_idxes[1])[3] == symbolnames[3]
            @test names(int_idxes[1])[4] == newnames[4]
            @test names(int_idxes[1])[4] != symbolnames[4]
            @test names(int_idxes[1])[5] != newnames[5]
            @test names(int_idxes[1])[5] == symbolnames[5]
        end
        @testset "getindex" begin
            idxer_one = int_idxes[1]
            @testset "single names" begin
                @test idxer_one[:First] == 1
                @test idxer_one[:Fifth] == 5
                @test idxer_one[:Third] == 3
                @test idxer_one["First"] == 1
                @test idxer_one["Fifth"] == 5
                @test idxer_one["Third"] == 3
                @test g_idxes[1][:First] == 1:5
                @test g_idxes[3][:First] == 1:5
            end
            @testset "vector of names" begin
                @test idxer_one[[:Seventh, :Third, :Ninth]] == [2, 3, 4]
                @test idxer_one[["Seventh", "Third", "Ninth"]] == [2, 3, 4]
                @test g_idxes[1][[:First, :Second]] == UnitRange{Int}[1:5, 6:10]
            end
            @testset "single integer" begin
                @test idxer_one[1] == 1
                @test g_idxes[1][4] == 16:20
            end
            @testset "multiple integers" begin
                @test idxer_one[1:3] == [1:3;]
                @test idxer_one[[1,3,5]] == [1,3,5]
            end
            @testset "booleans" begin
                @test idxer_one[[true, false, false, true, true]] == UInt[1, 4, 5]
            end
        end
    end
end

end


text:
# -------------------------------------------
# Process global projections and observations
# - Read altimetry observation (GSFC/MEASURES)
# - Read TG reconstruction ensemble (F2020)
# - Estimate trend + acceleration using Hector
# - Extrapolate data from tide gauges
# - Read NCA5 projections
# --------------------------------------------
module GlobalProjections
using NetCDF
using DelimitedFiles
using Statistics
using NCDatasets
using MAT
using Interpolations
using CSV
dir_code = homedir()*"/Projects/2021_NOAA_TR/Code/"
include(dir_code*"Hector.jl")
function RunGlobalProjections(settings)
    println("\nGlobal-mean observations, trajectories, and scenarios...")
    GMSL_20c = read_GMSL_20c(settings)
    GMSL_Trajectory = compute_GMSL_trajectory(GMSL_20c,settings)
    GMSL_Altimetry  = read_GMSL_Altimetry(settings)
    GMSL_NCA5       = read_NCA5(settings)
    correct_baselines!(GMSL_20c,GMSL_Trajectory,GMSL_Altimetry,GMSL_NCA5,settings)
    save_data(GMSL_20c,GMSL_Trajectory,GMSL_Altimetry,GMSL_NCA5,settings)
    println("Global-mean observations, trajectories, and scenarios done\n")
    return nothing
end

function read_GMSL_20c(settings)
    println("  Reading tide-gauge GMSL...")
    # Read F2020 GMSL
    fh = Dataset(settings["fn_gmsl_20c_ensembles"],"r")
    t = fh["time"][:]
    Λ = fh["likelihood"][:]
    η = fh["GMSL"][:]

    # Adjust to baseline
    idx_baseline = findall(in(settings["years_baseline"]),t)
    η .-= mean(η[idx_baseline,:],dims=1)
    # Mean and standard deviation
    η_mean = zeros(length(t),3)
    sidx = zeros(Int,length(Λ))
    𝚲    = zeros(Float64,length(Λ))
    Η    = zeros(Float64,length(Λ))
    for ts ∈ 1:length(t)
        sortperm!(sidx, η[ts,:])
        cumsum!(𝚲,Λ[sidx])
        @. Η = @views η[ts,sidx];
        η_mean[ts,1] = Η[findfirst(>(0.17),𝚲)]
        η_mean[ts,2] = Η[findfirst(>(0.5),𝚲)]
        η_mean[ts,3] = Η[findfirst(>(0.83),𝚲)]
    end
    GMSL_20c = Dict()
    GMSL_20c["years"] = t
    GMSL_20c["rsl"] = η_mean
    GMSL_20c["likelihood"] = Λ
    GMSL_20c["rsl_ensemble"] = η
    return GMSL_20c
end

function compute_GMSL_trajectory(GMSL_20c,settings)
    # -----------------------------------------------
    # Compute the trajectory for GMSL and extrapolate
    # -----------------------------------------------
    println("  Computing GMSL trajectories...")
    t_acc = findall(in(settings["years_trajectory_global"]),GMSL_20c["years"])

    # Compute trend and acceleration in each of the ensemble members of F2020 GMSL
    amat = ones(length(t_acc),3)
    amat[:,2] = GMSL_20c["years"][t_acc] .- mean(GMSL_20c["years"][t_acc])
    @. amat[:,3] = 0.5 * (GMSL_20c["years"][t_acc] - $mean(GMSL_20c["years"][t_acc]))^2
    amat_tr = transpose(amat)
    amat_sq = inv(amat_tr*amat)
    sol_arr = zeros(length(GMSL_20c["likelihood"]),3)
    amat_expand = ones(length(settings["years_trajectory_global"]),3)
    amat_expand[:,2] = settings["years_trajectory_global"] .- mean(GMSL_20c["years"][t_acc])
    @. amat_expand[:,3] = 0.5 * (settings["years_trajectory_global"] - $mean(GMSL_20c["years"][t_acc]))^2
    η_projected = zeros(length(settings["years_trajectory_global"]),length(GMSL_20c["likelihood"]))

    # Generate autocorrelated noise to account for uncertainties in the trajectory due to
    # internal variability and add this to the ensemble members
    noise = Hector.GenerateNoise(GMSL_20c["years"][t_acc],GMSL_20c["rsl"][t_acc,2],length(GMSL_20c["likelihood"]),length(GMSL_20c["years"]);accel=true,model="Powerlaw");
    GMSL_20c["rsl_ensemble"] .+= noise

    # Compute trend and acceleration for each ens member, and compute trajectory
    # time series: design matrix * solution
    for λ ∈ 1:length(GMSL_20c["likelihood"])
        sol_arr[λ,:] = amat_sq * (amat_tr*GMSL_20c["rsl_ensemble"][t_acc,λ])
        η_projected[:,λ] = amat_expand * sol_arr[λ,:]
    end

    # Statistics: compute mean and confidence intervals
    # of the extrapolated trajectory
    projection_mean = zeros(length(settings["years_trajectory_global"]),3)
    sidx = zeros(Int,length(GMSL_20c["likelihood"]))
    𝚲    = zeros(Float64,length(GMSL_20c["likelihood"]))
    Η    = zeros(Float64,length(GMSL_20c["likelihood"]))
    for ts ∈ 1:length(settings["years_trajectory_global"])
        sortperm!(sidx, η_projected[ts,:])
        cumsum!(𝚲,GMSL_20c["likelihood"][sidx])
        @. Η = @views η_projected[ts,sidx];
        projection_mean[ts,1] = Η[findfirst(>(0.17),𝚲)]
        projection_mean[ts,2] = Η[findfirst(>(0.5),𝚲)]
        projection_mean[ts,3] = Η[findfirst(>(0.83),𝚲)]
    end
    GMSL_Trajectory = Dict()
    GMSL_Trajectory["rsl"] = projection_mean
    GMSL_Trajectory["likelihood"] = GMSL_20c["likelihood"]

    # Statistics: get trend and acceleration of trajectory
    GMSL_Trajectory["trend"] = zeros(Float32,3)
    sidx = sortperm(sol_arr[:,2])
    𝚲 = cumsum(GMSL_20c["likelihood"][sidx])
    Η = sol_arr[sidx,2];
    GMSL_Trajectory["trend"][1] = Η[findfirst(>(0.17),𝚲)]
    GMSL_Trajectory["trend"][2] = Η[findfirst(>(0.5),𝚲)]
    GMSL_Trajectory["trend"][3] = Η[findfirst(>(0.83),𝚲)]

    GMSL_Trajectory["accel"] = zeros(Float32,3)
    sidx = sortperm(sol_arr[:,3])
    𝚲 = cumsum(GMSL_20c["likelihood"][sidx])
    Η = sol_arr[sidx,3];
    GMSL_Trajectory["accel"][1] = Η[findfirst(>(0.17),𝚲)]
    GMSL_Trajectory["accel"][2] = Η[findfirst(>(0.5),𝚲)]
    GMSL_Trajectory["accel"][3] = Η[findfirst(>(0.83),𝚲)]
    return GMSL_Trajectory
end

function read_GMSL_Altimetry(settings)
    println("  Reading altimetry GMSL...")
    # Read Altimetry GMSL from Goddard and apply corrections for
    # internal variability due to ENSO/PDO and take annual-means
    GMSL_Altimetry = Dict()
    GMSL_Altimetry["years"] = [1993:2020...]
      # Read data
    GSFC_raw = readdlm(settings["fn_gmsl_altimetry_GSFC"],skipstart=48)[:,[3,12]]
    ENSO_raw = matread(settings["fn_gmsl_ENSO"])
    # To monthly-mean GMSL
    tbounds = zeros(length(ENSO_raw["time"]),2)
    tbounds[:,1] = [1993.0:1/12:2021-1/12...]
    tbounds[:,2] = tbounds[:,1] .+ 1/12
    GSFC_monthly = zeros(Float32,length(ENSO_raw["time"]))
    [GSFC_monthly[t] = mean(GSFC_raw[(GSFC_raw[:,1] .>= tbounds[t,1]) .& (GSFC_raw[:,1] .<= tbounds[t,2]),2]) - ENSO_raw["gmsl"][1,t] for t ∈ 1:length(ENSO_raw["time"])]
     GMSL_Altimetry["alt"] =mean(reshape(GSFC_monthly,(12,:)),dims=1)[1,:] # Monthly-mean to annual-mean
    return GMSL_Altimetry
end

function read_NCA5(settings)
    println("  Reading NCA5 scenarios...")
    # Read the GMSL scenarios
    GMSL_NCA5 = Dict()
    GMSL_NCA5["years"] = ncread(settings["dir_NCA5"] * "NCA5_Low_gmsl.nc","years",start=[1],count=[-1])
    for scenario ∈ settings["NCA5_scenarios"]
        fn = settings["dir_NCA5"] * "NCA5_"*scenario*"_gmsl.nc"
        GMSL_NCA5[scenario] = convert.(Float32,ncread(fn,"total",start=[1,1],count=[-1,-1]))
    end
    return GMSL_NCA5
end

function correct_baselines!(GMSL_20c,GMSL_Trajectory,GMSL_Altimetry,GMSL_NCA5,settings)
    # Correct baseline: baseline for everything is trajectory value in 2000
    traj_baseline = findall(in(settings["years_baseline"]),settings["years_trajectory_global"])
    traj_nca5_start = findfirst(settings["years_trajectory_global"].==2005)

    GMSL_Trajectory["rsl"] .-= mean(GMSL_Trajectory["rsl"][traj_baseline,:])

    traj_bl = findall(in(intersect(settings["years_trajectory_global"],GMSL_20c["years"])),settings["years_trajectory_global"])
    obs_bl  = findall(in(intersect(settings["years_trajectory_global"],GMSL_20c["years"])),GMSL_20c["years"])
    GMSL_20c["rsl"] = GMSL_20c["rsl"] .- mean(GMSL_20c["rsl"][obs_bl,2]) .+ mean(GMSL_Trajectory["rsl"][traj_bl,2])

    traj_bl = findall(in(intersect(settings["years_trajectory_global"],GMSL_Altimetry["years"])),settings["years_trajectory_global"])
    obs_bl  = findall(in(intersect(settings["years_trajectory_global"],GMSL_Altimetry["years"])),GMSL_Altimetry["years"])
    GMSL_Altimetry["alt"] = GMSL_Altimetry["alt"] .- mean(GMSL_Altimetry["alt"][obs_bl]) .+ mean(GMSL_Trajectory["rsl"][traj_bl,2])

    for scenario in settings["NCA5_scenarios"]
        GMSL_NCA5[scenario] .+= GMSL_Trajectory["rsl"][traj_nca5_start,2]
    end
    return nothing
end

function save_data(GMSL_20c,GMSL_Trajectory,GMSL_Altimetry,GMSL_NCA5,settings)
    println("  Saving data...")
    fh = Dataset(settings["fn_proj_glb"],"c")
    fh.attrib["title"] = "Global projections"
    fh.attrib["description"] = "Global projections, trajectory, and observations for the interagency report: Global and Regional Sea Level Rise Scenarios for the United States: Updated Mean Projections and Extreme Water Level Probabilities Along U.S. Coastlines"
    defDim(fh,"years", length(settings["years"]))
    defDim(fh,"percentiles",3)

    defVar(fh,"years",settings["years"],("years",),deflatelevel=5,attrib = Dict("units" => "years"))
    defVar(fh,"percentiles",convert.(Int32,[17,50,83]),("percentiles",),deflatelevel=5, attrib = Dict("description" => "The percentile of the projection. The 50th percentile shows the median projection, and the 17th and 83rd show the upper- and lower bound on the 1 sigma level", "units" => "-"))

    # Write observations and trajectory
    defVar(fh,"GMSL_20c",Float32,("years","percentiles"),deflatelevel=5,attrib = Dict("description" => "Estimated 20th-century GMSL change (adapted from Frederikse et al. (2020)), relative to year 2000.", "units" => "mm"))[:] = LinearInterpolation((convert.(Float32,GMSL_20c["years"]),[1.0f0:3.0f0...]),GMSL_20c["rsl"],extrapolation_bc=NaN32)(settings["years"],[1.0f0:3.0f0...])
    defVar(fh,"GMSL_Altimetry",Float32,("years",),deflatelevel=5, attrib = Dict("description" => "Observed GMSL change over the altimetry era (NASA GSFC, Beckley et al. (2017)), relative to year 2000.", "units" => "mm"))[:] = LinearInterpolation((convert.(Float32,GMSL_Altimetry["years"])),GMSL_Altimetry["alt"],extrapolation_bc=NaN32)(settings["years"])
    defVar(fh,"GMSL_Trajectory",Float32,("years","percentiles"),deflatelevel=5, attrib = Dict("description" => "Estimated GMSL trajectory, relative to year 2000.", "units" => "mm"))[:] = LinearInterpolation((convert.(Float32,settings["years_trajectory_global"]),[1.0f0:3.0f0...]),GMSL_Trajectory["rsl"],extrapolation_bc=NaN32)(settings["years"],[1.0f0:3.0f0...])

    defVar(fh,"GMSL_trend", GMSL_Trajectory["trend"], ("percentiles",), deflatelevel=5, attrib = Dict("description" => "Linear trend of estimated trajectory", "units" => "mm yr-1"))
    defVar(fh,"GMSL_accel",GMSL_Trajectory["accel"],("percentiles",),deflatelevel=5, attrib = Dict("description" => "Acceleration (half quadratic) of estimated trajectory", "units" => "mm yr-2"))

    # Write scenarios
    for scenario ∈ settings["NCA5_scenarios"]
        defVar(fh,"GMSL_"*scenario,Float32,("years","percentiles"),deflatelevel=5, attrib = Dict("description" => "Projected sea level for scenario "*scenario*", relative to year 2000.", "units" => "mm"))[:] = LinearInterpolation((convert.(Float32,GMSL_NCA5["years"]),[1.0f0:3.0f0...]),GMSL_NCA5[scenario],extrapolation_bc=NaN32)(settings["years"],[1.0f0:3.0f0...])
    end
    close(fh)
    return nothing
end
end


text:
# Save images for trials from large_scale_weights.jl

using SubmodularMaximization
using PyPlot
using Base.Iterators
using Base.Threads
using JLD2
using Statistics
using HistogramFilters
using RosDataProcess

close("all")

pygui(false)

experiment_name = "large_scale_weights"
data_folder = "./data"

video_path = "./video"

experiment_file = joinpath(data_folder, experiment_name * ".jld2")

#@load experiment_file results

num_robots = 8:8:100

trials = 1:5
steps = 100
trial_steps = 20:steps

#trial = first(values(results))
#visualize_experiment(trial.trial_data, trial.configs)

for num_robots in num_robots
  println("Saving for ", num_robots, " robots")
  trial = results[num_robots, 1]

  data = trial.trial_data
  configs = trial.configs

  close("all")

  plot_state_space(configs.grid)
  ax = gca()
  ax.axis("off")

  folder = joinpath(video_path, string("num_robots_", num_robots))
  mkpath(folder)

  for (ii, step) in enumerate(data)
    file = joinpath(folder, string(ii, ".jpg"))

    if isfile(file)
      continue
    end

    println("  ", file)

    plots = visualize_time_step(robot_states=step.robot_states,
                                target_states=step.target_states,
                                target_filters=step.target_filters,
                                trajectories=step.trajectories)

    savefig(file, transparent=true, pad_inches=0.01, bbox_inches="tight")

    foreach(x->x.remove(), plots)
  end
end

# Record videos
rate = 2
for num_robots in num_robots
  @show folder = joinpath(video_path, string("num_robots_", num_robots))
  @show file = joinpath(video_path, string("num_robots_", num_robots, ".mp4"))

  video_command = `ffmpeg -hide_banner -y -r $(rate) -i $(folder)/%d.jpg -c:v libx264 -vf "scale=400:-1" $(file)`

  if !isfile(file)
    run(video_command)
  end
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
using MPIMeasurements

scanner = MPIScanner("HeadScanner.toml")
daq = getDAQ(scanner)

params = toDict(daq.params)

params["studyName"]="TestTobi"
params["studyDescription"]="A very cool measurement"
params["scannerOperator"]="Tobi"
params["dfStrength"]=[5e-3]  #[5e-3]
params["acqNumAverages"]=100
#params["calibFieldToVolt"]=[1]
#x = linspace(0,1,6)
#params["acqFFValues"] = [0.0]
#params["acqNumFFChannels"] = 2
#params["acqFFValues"] = repeat( x, inner=10)
#params["acqNumPeriodsPerFrame"]= div(length(params["acqFFValues"]),params["acqNumFFChannels"])

enableSlowDAC(daq,true)
@info "Starting Measurement..."
measurementCont(daq, params, controlPhase=false, showFT=true)


text:
##############
# Inputs

# Constant input current (mostly for testing purposes)
# the scaling is regulated here, and the weight is ignored
struct InputSimpleOffset <: NeuronType
  α::Float64 # scaling constant
end

# when the presynaptic is a simple input, just sum linearly to the input vector
function forward_signal!(tnow::Real,dt::Real,p_post::PopulationState,
    c::Connection,p_pre::PSSimpleInput{InputSimpleOffset})
  for i in eachindex(p_post.input)
    p_post.input[i] += p_pre.neurontype.α
  end
  return nothing
end

# if the receiver is a spiking type, update only when NOT refractory
function forward_signal!(tnow::Real,dt::Real,p_post::PSSpikingType,
    c::Connection,p_pre::PSSimpleInput{InputSimpleOffset})
  for i in eachindex(p_post.input)
    if ! p_post.isrefractory[i]
      p_post.input[i] += p_pre.neurontype.α
    end
  end
  return nothing
end

# independent Gaussian noise for each neuron
# weight is ignored
struct InputIndependentNormal <: NeuronType
  α::Float64 # scaling constant
end

# the connection is ignored: all neurons of the same population receive
# the same noise. Postsynaptic neurons must have a τ
function forward_signal!(tnow::Real,dt::Real,p_post::PopulationState,
    c::Connection,p_pre::PSSimpleInput{InputIndependentNormal})
  # std is α for isolated neuron
  _reg = sqrt(2*p_post.neurontype.τ/dt)
  for i in eachindex(p_post.input)
    p_post.input[i] += p_pre.neurontype.α*_reg*randn()
  end
  return nothing
end

# Arbitrary 1D time-dependent function
# injected directly as input current
# the input is regulated by the weight matrix

struct InputSomeFunction <: NeuronType
  f::Function # f(t) = input
end
function forward_signal!(tnow::Real,dt::Real,p_post::PopulationState,
    c::Connection,p_pre::PSSimpleInput{InputSomeFunction})
  in_t::Float64 = p_pre.neurontype.f(tnow)
  ws = nonzeros(c.weights)
  w_idx = rowvals(c.weights)
  for (i,w) in zip(w_idx,ws)
    p_post.input[i] += w*in_t
  end
  return nothing
end

# Fixed spiketrain as input
struct InputFixedSpiketrain <: NeuronType
  trains::Vector{Vector{Float64}} # one neuron -> one vector of spiketimes
  τ_output_decay::Float64
end
struct PSFixedSpiketrain{NT} <: PSSpikingType{NT where NT <: InputFixedSpiketrain}
  neurontype::NT
  n::Int64 # pop size
  isfiring::BitArray{1}
  counter::Vector{Int64} # keep track of elapsed spikes (instead of doing a full search)
end

function PSFixedSpiketrain(nt::InputFixedSpiketrain)
  n = length(nt.trains)
  isfiring=falses(n)
  counter=ones(n)
  return PSFixedSpiketrain{InputFixedSpiketrain}(nt,n,isfiring,counter)
end
function PSFixedSpiketrain(train::Vector{Vector{Float64}},τ::Real)
  return PSFixedSpiketrain(InputFixedSpiketrain(train,τ))
end
function reset_input!(ps::PSFixedSpiketrain)
  return nothing
end

function local_update!(t_now::Float64,dt::Float64,ps::PSFixedSpiketrain)
  reset_spikes!(ps)
  for neu in (1:ps.n) # WARNING : at most 1 spike per neuron in 1 dt
    # pick the spiketrain of each neuron, at the counter index
    c = ps.counter[neu]
    @inbounds if checkbounds(Bool,ps.neurontype.trains[neu],c)
      t_signal = ps.neurontype.trains[neu][c]
      if t_signal < t_now+dt
        ps.isfiring[neu] = true # neuron is firing
        ps.counter[neu] += 1    # move counter forward
      end
    end
  end
  return nothing
end
function reset!(ps::PSFixedSpiketrain)
  reset_spikes!(ps)
  fill!(ps.counter,1)
  return nothing
end

# Poisson firing
struct NTPoisson <: NeuronType
  rate::Ref{Float64} # rate is in Hz and is mutable
  τ_output_decay::Float64 # decay of postsynaptic conductance
end

struct PSPoisson{NT} <: PSSpikingType{NT}
  neurontype::NT
  n::Int64 # pop size
	isfiring::BitArray{1} # firing will be i.i.d. Poisson
  isfiring_alloc::Vector{Float64} # allocate probabilities
end
function PSPoisson(p::NT,n) where NT<:NeuronType
  return PSPoisson{NT}(p,n,falses(n),zeros(Float64,n))
end
function PSPoisson(rate::Real,τ_decay::Real,n::Integer)
  nt = NTPoisson(rate,τ_decay)
  return PSPoisson{NTPoisson}(nt,n,falses(n),zeros(Float64,n))
end
function reset_input!(::PSPoisson)
  return nothing
end
function reset!(::PSPoisson)
  return nothing
end
function local_update!(t_now::Float64,dt::Float64,ps::PSPoisson)
  reset_spikes!(ps)
  Random.rand!(ps.isfiring_alloc)
  _rate = ps.neurontype.rate[]
  @.  ps.isfiring = ps.isfiring_alloc < _rate*dt
  return nothing
end

# Poisson with frequency folloing some function

struct InputPoissonFt <: NeuronType
  ratefunction::Function # sigature ::Float64 -> Float64
  τ_output_decay::Float64
end
struct PSInputPoissonFt{NT} <: PSSpikingType{NT}
  neurontype::NT
  n::Int64 # pop size
	isfiring::BitArray{1}
  isfiring_alloc::Vector{Float64} # allocate probabilities
end
function PSInputPoissonFt(ratefun,τ,n)
  isfiring = falses(n)
  isfi_alloc = zeros(n)
  return PSInputPoissonFt(InputPoissonFt(ratefun,τ),n,isfiring,isfi_alloc)
end

function reset!(ps::PSInputPoissonFt)
  fill!(ps.isfiring,false)
  return nothing
end
function reset_input!(ps::PSInputPoissonFt)
  return nothing
end

function local_update!(t_now::Float64,dt::Float64,ps::PSInputPoissonFt)
  reset_spikes!(ps)
  _rate = ps.neurontype.ratefunction(t_now)
  Random.rand!(ps.isfiring_alloc)
  @.  ps.isfiring = ps.isfiring_alloc < _rate*dt
  return nothing
end

# Poisson with frequency following a *multi dimensional* function
# useful for input patterns

struct InputPoissonFtMulti <: NeuronType
  ratefunction::Function # signature ::Float64 -> ::Vector{Float64}
  τ_output_decay::Float64
end
struct PSInputPoissonFtMulti{NT} <: PSSpikingType{NT}
  neurontype::NT
  n::Int64 # pop size
	isfiring::BitArray{1}
  isfiring_alloc::Vector{Float64} # allocate probabilities
end
function PSInputPoissonFtMulti(ratefun,τ,n)
  isfiring = falses(n)
  isfi_alloc = zeros(n)
  return PSInputPoissonFtMulti(InputPoissonFtMulti(ratefun,τ),n,isfiring,isfi_alloc)
end
function reset!(ps::PSInputPoissonFtMulti)
  fill!(ps.isfiring,false)
  return nothing
end
function reset_input!(ps::PSInputPoissonFtMulti)
  return nothing
end
function local_update!(t_now::Float64,dt::Float64,ps::PSInputPoissonFtMulti)
  #reset_spikes!(ps)  not needed, I rewrite it fully
  Random.rand!(ps.isfiring_alloc)
  _rates::Vector{Float64} = ps.neurontype.ratefunction(t_now)
  @.  ps.isfiring = ps.isfiring_alloc < _rates*dt
  return nothing
end

# Training assemblies, let's do objects
struct PatternPresentation
  Npatt::Int64
  sequence::Vector{Int64} # sequence of patterns presented
  times::Vector{Float64}  # t_start1,t_start2... , size is length(sequence)+1
  patterns::Matrix{Float64} # columns are patterns, Npatt+1 columns: last column is null pattern
end
struct PatternPresentationContent # :-(  just to save it ...
  Npatt::Int64
  sequence::Vector{Int64}
  times::Vector{Float64}
  patterns::Matrix{Float64}
  function PatternPresentationContent(p::PatternPresentation)
    new(p.Npatt,p.sequence,p.times,p.patterns)
  end
end
get_content(p::PatternPresentation) = PatternPresentationContent(p)

function (patt::PatternPresentation)(t::Real,i::Integer)
  if patt.times[1] < t <= patt.times[end]
    it = searchsortedfirst(patt.times,t)-1
    return patt.patterns[i,patt.sequence[it]]
  else
    return patt.patterns[i,end]
  end
end

# helper function for
# training with patterns with no consecutive repetitions
function _generate_uniform_pattern_sequence(Npatterns::Integer,Ntrials::Integer)
  nreps = ceil(Integer,Ntrials/Npatterns)
  last_patt = 0
  _seq = Vector{Int64}[]
  for _ in 1:nreps
    _ss = shuffle(1:Npatterns)
    while _ss[1] == last_patt
      shuffle!(_ss)
    end
    last_patt = _ss[end]
    push!(_seq,_ss)
  end
  patt_seq = vcat(_seq...)[1:Ntrials]
  return patt_seq
end

function PatternPresentation(Δt::R,Ttot::R,
      low::R,high::R,
      npost::Integer,
      idxs_patternpop::Vector{Vector{Int64}} ;
      Δt_pattern_blank::R=0.0,
      t_pattern_delay::R=0.0) where R<:Real
  # time vector
  ts_pattern_start = collect(range(t_pattern_delay,Ttot;step=Δt+Δt_pattern_blank))
  # make pattern sequence
  Npatt = length(idxs_patternpop)
  patt_seq = _generate_uniform_pattern_sequence(Npatt,length(ts_pattern_start)-1)
  # must add blanks
  if Δt_pattern_blank > 0.0
    blank_seq = fill(Npatt+1,length(patt_seq))
    patt_seq = [transpose(hcat(patt_seq,blank_seq))...]
    patttimes = repeat(ts_pattern_start;inner=2)[2:end]
    for i in 2:2:length(patttimes)
      patttimes[i] -=Δt_pattern_blank
    end
  else
    patttimes = ts_pattern_start
  end
  # Generate full scale patterns
  fullp = fill(low,(npost,Npatt+1))
  for (i,neuidxs) in enumerate(idxs_patternpop)
    fullp[neuidxs,i] .= high
  end
  # generate the function and return it
  return PatternPresentation(Npatt,patt_seq,patttimes,fullp)
end

## generalize a little on the Inputs

abstract type SpikeGenerator end

struct SGPoisson <: SpikeGenerator
  rate::Float64
end
function generate_spikes!(t_now::Float64,dt::Float64,isfiring::BitArray{1},sg::SGPoisson)
  _rat = sg.rate*dt
  @inbounds @simd for i in eachindex(isfiring)
    isfiring[i] = rand() < _rat
  end
  return nothing
end
struct SGPoissonF <: SpikeGenerator
  ratefun::Function # t::Float64 -> rate::Float64
end
function generate_spikes!(t_now::Float64,dt::Float64,isfiring::BitArray{1},sg::SGPoissonF)
  _rat::Float64 = sg.ratefun(t_now)*dt
  if _rat > 0.0
    @inbounds @simd for i in eachindex(isfiring)
      isfiring[i] = rand() < _rat
    end
  else
    fill!(false,isfiring)
  end
  return nothing
end
struct SGPoissonMulti <: SpikeGenerator
  rates::Vector{Float64}
end
function generate_spikes!(t_now::Float64,dt::Float64,isfiring::BitArray{1},sg::SGPoissonMulti)
  @inbounds @simd for i in eachindex(isfiring)
    isfiring[i] = rand() < sg.rates[i]*dt
  end
  return nothing
end
struct SGPoissonMultiF <: SpikeGenerator
  ratefun::Function # t::Float64 -> rates::Vector{Float64}
end
function generate_spikes!(t_now::Float64,dt::Float64,isfiring::BitArray{1},sg::SGPoissonMultiF)
  _rates::Vector{Float64} = sg.ratefun(t_now)
  @inbounds @simd for i in eachindex(isfiring)
    if _rates[i] > 0
      isfiring[i] = rand() < _rates[i]*dt
    else
      isfiring[i] = false
    end
  end
  return nothing
end

# for exact spike generation, I need an upper bound for the instantaneous rate
struct SGPoissonFExact <: SpikeGenerator
  ratefunction::Function # (t::Float64,i::Int64) -> rate_i::Float64
  ratefunction_upper::Function # (t::Float64,i::Int64) -> rate_i::Float64
end

struct SGTrains <: SpikeGenerator
  trains::Vector{Vector{Float64}}
  counter::Vector{Int64}
  function SGTrains(trains)
    counter = fill(1,length(trains))
    new(trains,counter)
  end
end
function generate_spikes!(t_now::Float64,dt::Float64,isfiring::BitArray{1},sg::SGTrains)
  fill!(isfiring,false)
  for i in eachindex(isfiring)
    # pick the spiketrain of each neuron, at the counter index
    # I could have used searchsortedfirst, too
    c = sg.counter[i]
    if checkbounds(Bool,sg.trains[i],c)
      @inbounds t_next_spike = sg.trains[i][c]
      if t_next_spike < t_now+dt
        isfiring[i] = true # neuron is firing
        sg.counter[i] += 1    # move counter forward
      end
    end
  end
  return nothing
end

struct NTInputConductance{SK<:SynapticKernel,SG<:SpikeGenerator} <:NTConductance
  spikegenerator::SG
  synaptic_kernel::SK
  v_reversal::Float64 # reversal potential that affects postsynaptic neurons
end
struct PSInputConductance{NT<:NTInputConductance} <: PSSpikingType{NT}
  n::Int64
  neurontype::NT
  isfiring::BitArray{1}
  function PSInputConductance(nt::NT,n_neu::Int64) where NT<:NeuronType
    isfiring=falses(n_neu)
    new{NT}(n_neu,nt,isfiring)
  end
end
reset_input!(::PSInputConductance) = nothing
function reset!(::PSInputConductance{NT}) where {SK,SG,NT<:NTInputConductance{SK,SG}}
  return nothing
end
function reset!(ps::PSInputConductance{NT}) where {SK,NT<:NTInputConductance{SK,SGTrains}}
  fill!(ps.neurontype.spikegenerator.counter,1)
  return nothing
end
function local_update!(t_now::Float64,dt::Float64,ps::PSInputConductance)
  generate_spikes!(t_now,dt,ps.isfiring,ps.neurontype.spikegenerator)
  return nothing
end

#####
# here I also supplement the Forward signal , for simplicity and efficiency
# works with FakeConnection only
# downside : one needs to include the population weight here (vector of weights)
# other downside : NO PLASTICITY
# note that input could be correlated, if I define the spike generation in
# NTInputConductance accordingly ... but it becomes hacky

struct PSInputPoissonConductance{NT<:NTInputConductance} <: PSSpikingType{NT}
  n::Int64 # Must be the same as npost (n neurons that receive the input)
  neurontype::NT
  isfiring::BitArray{1}
  input_weights::Vector{Float64}
  trace1::Vector{Float64}  # alas! traces are still needed
  trace2::Vector{Float64}
  function PSInputPoissonConductance(nt::NT,weights::Vector{Float64}) where NT<:NeuronType
    n_neu = length(weights)
    isfiring=falses(n_neu)
    tr1 = zeros(n_neu)
    tr2 = zeros(n_neu)
    new{NT}(n_neu,nt,isfiring,weights,tr1,tr2)
  end
  function PSInputPoissonConductance(nt::NT,weight::Float64,n_neu::Integer) where NT<:NeuronType
    weights = fill(weight,n_neu)
    isfiring=falses(n_neu)
    tr1 = zeros(n_neu)
    tr2 = zeros(n_neu)
    new{NT}(n_neu,nt,isfiring,weights,tr1,tr2)
  end
end
function reset!(ps::PSInputPoissonConductance{NT}) where {SK,SG,NT<:NTInputConductance{SK,SG}}
  fill!(ps.trace1,0.0)
  fill!(ps.trace2,0.0)
  fill!(ps.isfiring,false)
  # special treatment for train inputs
  if hasproperty(ps.neurontype.spikegenerator,:counter)
    fill!(ps.neurontype.spikegenerator.counter,1)
  end
  return nothing
end
function reset_input!(::PSInputPoissonConductance)
  error("""
  If you are modeling a neural population explicitly, please
  use PSInputConductance instead of PSInputPoissonConductance
  to define it.
    """)
  return nothing
end
# function reset!(ps::PSInputPoissonConductance{NT}) where {SK,NT<:NTInputConductance{SK,SGTrains}}
#   fill!(ps.trace1,0.0)
#   fill!(ps.trace2,0.0)
#   fill!(ps.isfiring,false)
#   return nothing
# end

@inline function trace_decay!(dt::Real,ps::PSInputPoissonConductance)
  trace_decay!(dt,ps.trace1,ps.trace2,ps.neurontype.synaptic_kernel)
end

function forward_signal!(t_now::Real,dt::Real,
      pspost::PSLIFConductance,::FakeConnection,
      pspre::PSInputPoissonConductance)
  preneu = pspre.neurontype
  pre_v_reversal = preneu.v_reversal
  pre_synker = preneu.synaptic_kernel
  # traces time decay (here or at the end? meh)
  trace_decay!(dt,pspre)
  # update spikes
  generate_spikes!(t_now,dt,pspre.isfiring,preneu.spikegenerator)
  @inbounds for (i,firing) in enumerate(pspre.isfiring)
    # update trace if firing
    if firing
      # increment traces of spiking neurons
      # function defined in lif_conductance.jl
      trace_spike_update!(pspre.input_weights[i],pspre.trace1,pspre.trace2,pre_synker,i)
    end
    # update ALL inputs with synaptic traces, unless refractory
    if !pspost.isrefractory[i]
      # function defined in lif_conductance.jl
      postsynaptic_kernel_update!(pspost.input,pspost.state_now,
          pspre.trace1,pspre.trace2,pre_synker,pre_v_reversal,i)
    end
  end
  return nothing
end

#####
# but what happens when the input rate is >> dt ?
# reducing dt would be extremely unefficient... so
# NEW input implementation wiht EXACT spike timing
# and multiple update of traces when needed

# Thinning algorith, e.g.  Laub,Taimre,Pollet 2015
function _rand_by_thinning(t_start::Real,get_rate::Function,get_rate_upper::Function;
    Tmax=50.0,nowarning::Bool=false)
  t = t_start
  while (t-t_start)<Tmax # Tmax is upper limit, if rate too low
    (rup::Float64) = get_rate_upper(t)
    Δt = rand(Exponential())./rup
    t = t+Δt
    u = rand(Uniform(0.0,rup))
    (_r::Float64) = get_rate(t)
    if u <= _r
      return t
    end
  end
  # if we are above Tmax, just return upper limit
  if !nowarning
    @warn "Upper limit reached, input firing below $(inv(Tmax)) Hz"
  end
  return Tmax + t_start
end

struct PSInputPoissonConductanceExact{NT<:NTInputConductance} <: PSSpikingType{NT}
  n::Int64 # Must be the same as npost (n neurons that receive the input)
  neurontype::NT
  firingtimes::Vector{Float64}
  input_weights::Vector{Float64}
  trace1::Vector{Float64}  # alas! traces are still needed
  trace2::Vector{Float64}
  function PSInputPoissonConductanceExact(nt::NT,weights::Vector{Float64}) where NT<:NeuronType
    n_neu = length(weights)
    tr1,tr2,firingtimes = ntuple(_->Vector{Float64}(undef,n_neu),3)
    ret = new{NT}(n_neu,nt,firingtimes,weights,tr1,tr2)
    reset!(ret)
    return ret
  end
  function PSInputPoissonConductanceExact(nt::NT,weight::Float64,n_neu::Integer) where NT<:NeuronType
    weights = fill(weight,n_neu)
    tr1,tr2,firingtimes = ntuple(_->Vector{Float64}(undef,n_neu),3)
    ret = new{NT}(n_neu,nt,firingtimes,weights,tr1,tr2)
    reset!(ret)
    return ret
  end
end

function _get_spiketime_update(t_current_spike::Real,sg::SGPoisson,::Integer)
  return t_current_spike + rand(Exponential())/sg.rate
end

function _get_spiketime_update(t_current_spike::Real,sg::SGPoissonMulti,i::Integer)
  return t_current_spike + rand(Exponential())/sg.rates[i]
end

function _get_spiketime_update(::Real,sg::SGTrains,i::Integer)
  # note that t_current_spike is expected to be
  # sg.trains[i][counter[i]] (before counter update)
  c = sg.counter[i]
  if checkbounds(Bool,sg.trains[i],c+1)
    sg.counter[i] = c+1    # move counter forward
    return sg.trains[i][c+1] # return next spike time
  else
    return Inf
  end
end
function _get_spiketime_update(::Real,sg::SGPoissonMultiF,i::Integer)
  return error("please define input using `SGPoissonFExact`")
end

function _get_spiketime_update(t_current_spike::Real,sg::SGPoissonFExact,i::Integer)
  return _rand_by_thinning(t_current_spike,
    t->sg.ratefunction(t,i),t->sg.ratefunction_upper(t,i))
end

# in this case reset is an initialization step...
# which means I need to compute all first spike proposals

function reset!(ps::PSInputPoissonConductanceExact{NT}) where {SK,SG,NT<:NTInputConductance{SK,SG}}
  fill!(ps.trace1,0.0)
  fill!(ps.trace2,0.0)
  sgen = ps.neurontype.spikegenerator
  for i in eachindex(ps.firingtimes)
    ps.firingtimes[i] = _get_spiketime_update(0.0,sgen,i)
  end
  # special treatment for train inputs
  if hasproperty(sgen,:counter)
    fill!(sgen.counter,1)
  end
  return nothing
end

@inline function trace_decay!(dt::Real,ps::PSInputPoissonConductanceExact)
  trace_decay!(dt,ps.trace1,ps.trace2,ps.neurontype.synaptic_kernel)
end

function forward_signal!(t_now::Real,dt::Real,
      pspost::PSLIFConductance,::FakeConnection,
      pspre::PSInputPoissonConductanceExact)
  preneu = pspre.neurontype
  pre_v_reversal = preneu.v_reversal
  pre_synker = preneu.synaptic_kernel
  sgen = pspre.neurontype.spikegenerator
  # traces time decay (here or at the end? meh)
  trace_decay!(dt,pspre)
  # if t_now moved past a spiketime...
  @inbounds for i in eachindex(pspre.firingtimes)
    tspike = pspre.firingtimes[i]
    _n_spikes = 0
    # count spike occurrences
    while tspike <= t_now
      _n_spikes += 1
      tspike = _get_spiketime_update(tspike,sgen,i)
    end
    # increment traces of spiking neurons
    # ASSUMING IT IS LINEAR
    if _n_spikes > 0
    # function defined in lif_conductance.jl
      trace_spike_update!(_n_spikes*pspre.input_weights[i],
          pspre.trace1,pspre.trace2,pre_synker,i)
      pspre.firingtimes[i]=tspike
    end
  end
  @inbounds @simd for i in eachindex(pspost.input)
    if ! pspost.isrefractory[i]
    # function defined in lif_conductance.jl
    postsynaptic_kernel_update!(pspost.input,pspost.state_now,
            pspre.trace1,pspre.trace2,pre_synker,pre_v_reversal,i)
    end
  end
  return nothing
end

#################################

# Inhibitory stabilization

# from Fiete et al , 2010 , Neuron
# not a firing neuron type! More like an input type
# all constants here, works with FakeConnection only
struct NTConductanceInputInhibitoryStabilization <:NeuronType
  v_reversal::Float64
  Aglo::Float64
  Aloc::Float64
  τglo::Float64
  τloc::Float64
end
struct PSConductanceInputInhibitionStabilization <: PopulationState{NTConductanceInputInhibitoryStabilization}
  neurontype::NTConductanceInputInhibitoryStabilization
  n::Int64 # size of stabilized E population
  post_loc_traces::Vector{Float64}
  post_glo_trace::Ref{Float64}
end
function PSConductanceInputInhibitionStabilization(
    v_reversal::R,Aglo::R,Aloc::R,τglo::R,τloc::R,n::Integer) where R
  nt = NTConductanceInputInhibitoryStabilization(v_reversal,Aglo,Aloc,τglo,τloc)
  loc_traces = zeros(n)
  glo_trace = Ref(0.0)
  return PSConductanceInputInhibitionStabilization(nt,n,
    loc_traces,glo_trace)
end

function reset!(ps::PSConductanceInputInhibitionStabilization)
  fill!(ps.post_loc_traces,0.0)
  ps.post_glo_trace[]=0.0
  return nothing
end

# traces decay in time
@inline function trace_decay!(dt::Real,ps::PSConductanceInputInhibitionStabilization)
  ps.post_loc_traces .*= exp(-dt/ps.neurontype.τloc)
  ps.post_glo_trace[] *= exp(-dt/ps.neurontype.τglo)
end

# adds 1 to traces
function trace_spike_update!(ps::PSConductanceInputInhibitionStabilization,
    isfiring::BitArray{1})
  @inbounds for (i,firing) in enumerate(isfiring)
    if firing
      ps.post_loc_traces[i] += 1.0
      ps.post_glo_trace[] += 1.0
    end
  end
  return nothing
end

# This function should never be called anyway!
# because this PS does not need to appear as pre in a Population object
function local_update!(::Float64,::Float64,
      ::PSConductanceInputInhibitionStabilization)
  return nothing # nothing to update, traces decay in the forward signal function
end

function forward_signal!(t_now::Real,dt::Real,
      pspost::PSLIFConductance,::FakeConnection,
      pspre::PSConductanceInputInhibitionStabilization)
  preneu = pspre.neurontype
	# add the currents to postsynaptic input
	# ONLY non-refractory neurons
	@inbounds @simd for i in eachindex(pspost.input)
    if ! pspost.isrefractory[i]
      ker_loc = preneu.Aloc*pspre.post_loc_traces[i]/preneu.τloc
      ker_glo = preneu.Aglo*pspre.post_glo_trace[]/(pspre.n*preneu.τglo)
	    pspost.input[i] += (ker_loc+ker_glo)*(preneu.v_reversal - pspost.state_now[i])
		end
	end
  # traces time decay
  trace_decay!(dt,pspre)
  # update traces with spiking neurons
	trace_spike_update!(pspre,pspost.isfiring)
  # all done !
  return nothing
end


text:
module OntoNotes

type Sentence
    plain::UTF8String
    tokens::Vector{UTF8String}
    tree
    leaves::Vector{UTF8String}
end

#import Base.show
#show(io::IO, s::Sentence) = show(s.plain, s.tokens, s.tree, s.leaves

"Read .onf file and convert it to Sentence."
function readonf(path)
    text = open(readall, path)
    reg = r"Plain sentence:\n-+?\n(.+?)Treebanked sentence:\n-+?\n(.+?)Tree:\n-+?\n(.+?)Leaves:\n-+?\n(.+?)\n\n"s
    map(eachmatch(reg, text)) do m
        length(m.captures) == 4 || error("invalid format")
        plain = strip(m.captures[1])
        tokens = split(strip(m.captures[2]))
        tree = format_tree(m.captures[3])
        leaves = split(m.captures[4], '\n', keep=false)
        leaves = map(strip, leaves)
        Sentence(plain, tokens, tree, leaves)
    end
end

function format_tree(s)
    s = replace(s, "\n", "")
    s = replace(s, r"\s+\(", "(")
    s = replace(s, " ", "_") # concat preterminal and terminal
    s = replace(s, ",", "comma")
    s = replace(s, ".", "period")
    s = replace(s, "-", "_")
    s
end

end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
const ContinuousMultivariateSampleable = Sampleable{Multivariate, Continuous}

distribution(d::Sampleable) = d

"""
    CenteredAmbiguitySet

Defines the ambiguity related to the distribution of a random variable, often used in
Robust Optimization (RO) and Distributionally Robust Optimization (DRO) problems.
It represents a bounded infinite set of distributions.
"""
abstract type CenteredAmbiguitySet{T<:Real, D<:Sampleable} end

Base.length(s::CenteredAmbiguitySet) = length(distribution(s))

mean(s::CenteredAmbiguitySet) = mean(distribution(s))
cov(s::CenteredAmbiguitySet) = cov(distribution(s))

"""
    AmbiguitySet

Alias for Union{CenteredAmbiguitySet, ContinuousMultivariateSampleable}
"""
AmbiguitySet = Union{CenteredAmbiguitySet, Sampleable}

max_stars_repo_name: andrewrosemberg/PortfolioOpt

max_issues_repo_name: andrewrosemberg/PortfolioOpt

max_forks_repo_name: andrewrosemberg/PortfolioOpt


text:
using Distributions, NLopt, ChaosCommunications

carrier = LogisticCarrier(4)

noise_var = ebn0db2var(5., system=:UTURUniChannelCSK, carrier=carrier)

println("    Testing correlation decoder of coherent single-user CSK...")

system = UTURUniChannelCSK(true, carrier, noise=Normal(0., sqrt(noise_var)))

sim_ber(system, 1, 100000)
# psim_ber(system, 1, 100000)

println("    Testing correlation decoder of non-coherent single-user CSK...")

system = UTURUniChannelCSK(false, carrier, noise=Normal(0., sqrt(noise_var)))

sim_ber(system, 1, 100000)

println("    Testing MCML decoder of non-coherent single user-CSK...")

decoder = mcml_decoder(system=:UTURUniChannelCSK, carrier=carrier)

system = UTURUniChannelCSK(false, carrier, noise=Normal(0., sqrt(noise_var)), decoder=decoder)

sim_ber(system, 1, 1000)