text:
using AlphaGo
using AlphaGo: from_kgs, from_sgf, BLACK, WHITE, EMPTY, MISSING_GROUP_ID,
              GoPosition, is_koish, is_eyeish, PlayerMove, is_move_legal,
              is_move_suicidal, score, from_board, flip_playerturn!, add_stone!,
              all_legal_moves, IllegalMove, from_flat, pass_move!, play_move!
using Base.Test

include("test_utils.jl")

pc_set(string, env::GoEnv) = Set(map(from_kgs, split(string),
                                        repmat([env], length(split(string)))))

@testset "go" begin
  env = GoEnv(9)
  EMPTY_ROW = repeat(".", env.N) * "\n"
  TEST_BOARD = load_board("""
                          .X.....OO
                          X........
                          """ * repeat(EMPTY_ROW, 7), env)
    @testset "load_board" begin
    @test assertEqualArray(env.EMPTY_BOARD, zeros(env.N, env.N))
    @test assertEqualArray(env.EMPTY_BOARD, load_board(repeat(". \n", env.N ^ 2), env))
  end

  @testset "parsing" begin
    @test from_kgs("A9", env) == (1, 1)
    @test from_sgf("aa") == (1, 1)
    @test from_kgs("A3", env) == (7, 1)
    @test from_sgf("ac") == (3, 1)
    @test from_kgs("D4", env) == from_sgf("df")
  end

  @testset "neighbors" begin
    corner = from_kgs("A1", env)
    neighbors = [env.EMPTY_BOARD[c...] for c in env.NEIGHBORS[corner...]]
    @test length(neighbors) == 2

    side = from_kgs("A2", env)
    side_neighbors = [env.EMPTY_BOARD[c...] for c in env.NEIGHBORS[side...]]
    @test length(side_neighbors) == 3
  end

  @testset "is_koish" begin
    @test is_koish(TEST_BOARD, from_kgs("A9", env), env) == BLACK
    @test is_koish(TEST_BOARD, from_kgs("B8", env), env) == nothing
    @test is_koish(TEST_BOARD, from_kgs("B9", env), env) == nothing
    @test is_koish(TEST_BOARD, from_kgs("E5", env), env) == nothing
  end

  @testset "is_eyeish" begin
    board = load_board("""
              .XX...XXX
              X.X...X.X
              XX.....X.
              ........X
              XXXX.....
              OOOX....O
              X.OXX.OO.
              .XO.X.O.O
              XXO.X.OO.
          """, env)
    B_eyes = pc_set("A2 A9 B8 J7 H8", env)
    W_eyes = pc_set("H2 J1 J3", env)
    not_eyes = pc_set("B3 E5", env)
    for be in B_eyes
      @test is_eyeish(board, be, env) == BLACK ? true : throw(AssertionError(string(be)))
    end
    for we in W_eyes
      @test is_eyeish(board, we, env) == WHITE ? true : throw(AssertionError(string(we)))
    end
    for ne in not_eyes
      @test is_eyeish(board, ne, env) == nothing ? true : throw(AssertionError(string(ne)))
    end
  end
  @testset "lib_tracker_init" begin
    board = load_board("X........" * repeat(EMPTY_ROW, 8), env)

    lib_tracker = from_board(board, env)
    @test length(lib_tracker.groups) == 1
    @test lib_tracker.group_index[from_kgs("A9", env)...] != MISSING_GROUP_ID
    @test lib_tracker.liberty_cache[from_kgs("A9", env)...] == 2
    sole_group = lib_tracker.groups[lib_tracker.group_index[from_kgs("A9", env)...]]
    @test sole_group.stones == pc_set("A9", env)
    @test sole_group.liberties == pc_set("B9 A8", env)
    @test sole_group.color == BLACK
  end

  @testset "place_stone" begin
    board = load_board("X........" * repeat(EMPTY_ROW, 8), env)
    lib_tracker = from_board(board, env)
    add_stone!(lib_tracker, BLACK, from_kgs("B9", env), env)
    @test length(lib_tracker.groups) == 1
    @test lib_tracker.group_index[from_kgs("A9", env)...] != MISSING_GROUP_ID
    @test lib_tracker.liberty_cache[from_kgs("A9", env)...] == 3
    @test lib_tracker.liberty_cache[from_kgs("B9", env)...] == 3
    sole_group = lib_tracker.groups[lib_tracker.group_index[from_kgs("A9", env)...]]
    @test sole_group.stones == pc_set("A9 B9", env)
    @test sole_group.liberties == pc_set("C9 A8 B8", env)
    @test sole_group.color == BLACK
  end

  @testset "place_stone_opposite_color" begin
    board = load_board("X........" * repeat(EMPTY_ROW, 8), env)
    lib_tracker = from_board(board, env)
    add_stone!(lib_tracker, WHITE, from_kgs("B9", env), env)
    @test length(lib_tracker.groups) == 2
    @test lib_tracker.group_index[from_kgs("A9", env)...] != MISSING_GROUP_ID
    @test lib_tracker.group_index[from_kgs("B9", env)...] != MISSING_GROUP_ID
    @test lib_tracker.liberty_cache[from_kgs("A9", env)...] == 1
    @test lib_tracker.liberty_cache[from_kgs("B9", env)...] == 2
    black_group = lib_tracker.groups[lib_tracker.group_index[from_kgs("A9", env)...]]
    white_group = lib_tracker.groups[lib_tracker.group_index[from_kgs("B9", env)...]]
    @test black_group.stones == pc_set("A9", env)
    @test black_group.liberties == pc_set("A8", env)
    @test black_group.color == BLACK
    @test white_group.stones == pc_set("B9", env)
    @test white_group.liberties == pc_set("C9 B8", env)
    @test white_group.color == WHITE
  end

  @testset "merge_multiple_groups" begin
    board = load_board("""
        .X.......
        X.X......
        .X.......
    """ * repeat(EMPTY_ROW, 6), env)
    lib_tracker = from_board(board,env)
    add_stone!(lib_tracker, BLACK, from_kgs("B8", env), env)
    @test length(lib_tracker.groups) == 1
    @test lib_tracker.group_index[from_kgs("B8", env)...] != MISSING_GROUP_ID
    sole_group = lib_tracker.groups[lib_tracker.group_index[from_kgs("B8", env)...]]
    @test sole_group.stones == pc_set("B9 A8 B8 C8 B7", env)
    @test sole_group.liberties == pc_set("A9 C9 D8 A7 C7 B6", env)
    @test sole_group.color == BLACK

    liberty_cache = lib_tracker.liberty_cache
    for stone in sole_group.stones
      @test liberty_cache[stone...] == 6 ? true : throw(AssertionError(string(stone)))
    end
  end

  @testset "capture_stone" begin
    board = load_board("""
        .X.......
        XO.......
        .X.......
    """ * repeat(EMPTY_ROW, 6), env)
    lib_tracker = from_board(board, env)
    captured = add_stone!(lib_tracker, BLACK, from_kgs("C8", env), env)
    @test length(lib_tracker.groups) == 4
    @test lib_tracker.group_index[from_kgs("B8", env)...] == MISSING_GROUP_ID
    @test captured == pc_set("B8", env)
  end

  @testset "capture_many" begin
    board = load_board("""
        .XX......
        XOO......
        .XX......
    """ * repeat(EMPTY_ROW, 6), env)
    lib_tracker = from_board(board, env)
    captured = add_stone!(lib_tracker, BLACK, from_kgs("D8", env), env)
    @test length(lib_tracker.groups) == 4
    @test lib_tracker.group_index[from_kgs("B8", env)...] == MISSING_GROUP_ID
    @test captured == pc_set("B8 C8", env)

    left_group = lib_tracker.groups[lib_tracker.group_index[from_kgs("A8", env)...]]
    @test left_group.stones == pc_set("A8", env)
    @test left_group.liberties == pc_set("A9 B8 A7", env)

    right_group = lib_tracker.groups[lib_tracker.group_index[from_kgs("D8", env)...]]
    @test right_group.stones == pc_set("D8", env)
    @test right_group.liberties == pc_set("D9 C8 E8 D7", env)

    top_group = lib_tracker.groups[lib_tracker.group_index[from_kgs("B9", env)...]]
    @test top_group.stones == pc_set("B9 C9", env)
    @test top_group.liberties == pc_set("A9 D9 B8 C8", env)

    bottom_group = lib_tracker.groups[lib_tracker.group_index[from_kgs("B7", env)...]]
    @test bottom_group.stones == pc_set("B7 C7", env)
    @test bottom_group.liberties == pc_set("B8 C8 A7 D7 B6 C6", env)

    liberty_cache = lib_tracker.liberty_cache
    for stone in top_group.stones
      @test liberty_cache[stone...] == 4 ? true : throw(AssertionError(string(stone)))
    end
    for stone in left_group.stones
      @test liberty_cache[stone...] == 3 ? true: throw(AssertionError(string(stone)))
    end
    for stone in right_group.stones
      @test liberty_cache[stone...] == 4 ? true : throw(AssertionError(string(stone)))
    end
    for stone in bottom_group.stones
      @test liberty_cache[stone...] == 6 ? true : throw(AssertionError(string(stone)))
    end
    for stone in captured
      @test liberty_cache[stone...] == 0 ? true : throw(AssertionError(string(stone)))
    end
  end

  @testset "capture_multiple_groups" begin
    board = load_board("""
        .OX......
        OXX......
        XX.......
    """ * repeat(EMPTY_ROW, 6), env)
    lib_tracker = from_board(board, env)
    captured = add_stone!(lib_tracker, BLACK, from_kgs("A9", env), env)
    @test length(lib_tracker.groups) == 2
    @test captured == pc_set("B9 A8", env)

    corner_stone = lib_tracker.groups[lib_tracker.group_index[from_kgs("A9", env)...]]
    @test corner_stone.stones == pc_set("A9", env)
    @test corner_stone.liberties == pc_set("B9 A8", env)

    surrounding_stones = lib_tracker.groups[lib_tracker.group_index[from_kgs("C9", env)...]]
    @test surrounding_stones.stones == pc_set("C9 B8 C8 A7 B7", env)
    @test surrounding_stones.liberties == pc_set("B9 D9 A8 D8 C7 A6 B6", env)

    liberty_cache = lib_tracker.liberty_cache
    for stone in corner_stone.stones
      @test liberty_cache[stone...] == 2 ? true : throw(AssertionError(string(stone)))
    end
    for stone in surrounding_stones.stones
      @test liberty_cache[stone...] == 7 ? true : throw(AssertionError(string(stone)))
    end
  end

  @testset "same_friendly_group_neighboring_twice" begin
    board = load_board("""
        XX.......
        X........
    """ * repeat(EMPTY_ROW, 7), env)

    lib_tracker = from_board(board, env)
    captured = add_stone!(lib_tracker, BLACK, from_kgs("B8", env), env)
    @test length(lib_tracker.groups) == 1
    sole_group_id = lib_tracker.group_index[from_kgs("A9", env)...]
    sole_group = lib_tracker.groups[sole_group_id]
    @test sole_group.stones == pc_set("A9 B9 A8 B8", env)
    @test sole_group.liberties == pc_set("C9 C8 A7 B7", env)
    @test captured == Set()
  end

  @testset "same_opponent_group_neighboring_twice" begin
    board = load_board("""
        XX.......
        X........
    """ * repeat(EMPTY_ROW, 7), env)

    lib_tracker = from_board(board, env)
    captured = add_stone!(lib_tracker, WHITE, from_kgs("B8", env), env)
    @test length(lib_tracker.groups) == 2
    black_group = lib_tracker.groups[lib_tracker.group_index[from_kgs("A9", env)...]]
    @test black_group.stones == pc_set("A9 B9 A8", env)
    @test black_group.liberties == pc_set("C9 A7", env)

    white_group = lib_tracker.groups[lib_tracker.group_index[from_kgs("B8", env)...]]
    @test white_group.stones == pc_set("B8", env)
    @test white_group.liberties == pc_set("C8 B7", env)

    @test captured == Set()
  end

  @testset "passing" begin
    start_position = GoPosition(env,
        board = TEST_BOARD,
        n = 0,
        komi = 6.5,
        caps = (1, 2),
        ko = PlayerMove(from_kgs("A1", env)...),
        recent = Vector{PlayerMove}(),
        to_play = BLACK
    )
    expected_position = GoPosition(env,
        board = TEST_BOARD,
        n = 1,
        komi = 6.5,
        caps = (1, 2),
        ko = nothing,
        recent = [PlayerMove(BLACK, nothing)],
        to_play = WHITE
    )
    pass_position = pass_move!(start_position)
    @test assertEqualPositions(pass_position, expected_position)
  end

  @testset "flipturn" begin
    start_position = GoPosition(env,
        board = TEST_BOARD,
        n = 0,
        komi = 6.5,
        caps = (1, 2),
        ko = PlayerMove(from_kgs("A1", env)...),
        recent = Vector{PlayerMove}(),
        to_play = BLACK
    )
    expected_position = GoPosition(env,
        board = TEST_BOARD,
        n = 0,
        komi = 6.5,
        caps = (1, 2),
        ko = nothing,
        recent = Vector{PlayerMove}(),
        to_play = WHITE
    )
    flip_position = flip_playerturn!(start_position)
    @test assertEqualPositions(flip_position, expected_position)
  end

  @testset "is_move_suicidal" begin
    board = load_board("""
        ...O.O...
        ....O....
        XO.....O.
        OXO...OXO
        O.XO.OX.O
        OXO...OOX
        XO.......
        ......XXO
        .....XOO.
    """, env)
    position = GoPosition(env,
        board = board,
        to_play = BLACK
    )
    suicidal_moves = pc_set("E9 H5", env)
    nonsuicidal_moves = pc_set("B5 J1 A9", env)
    for move in suicidal_moves
      @test position.board[move...] == EMPTY #sanity check my coordinate input
      @test is_move_suicidal(position, move) ? true : throw(AssertionError(string(move)))
    end
    for move in nonsuicidal_moves
      @test position.board[move...] == EMPTY # sanity check my coordinate input
      @test !is_move_suicidal(position, move) ? true : throw(AssertionError(string(move)))
    end
  end

  @testset "legal_moves" begin
    board = load_board("""
        .O.O.XOX.
        O..OOOOOX
        ......O.O
        OO.....OX
        XO.....X.
        .O.......
        OX.....OO
        XX...OOOX
        .....O.X.
    """, env)
    position = GoPosition(env,
        board = board,
        to_play = BLACK)
    illegal_moves = pc_set("A9 E9 J9", env)
    legal_moves = pc_set("A4 G1 J1 H7", env)
    for move in illegal_moves
      @test !is_move_legal(position, move)
    end
    for move in legal_moves
      @test is_move_legal(position, move)
    end
    # check that the bulk legal test agrees with move-by-move illegal test.
    bulk_legality = all_legal_moves(position)
    for (i, bulk_legal) in enumerate(bulk_legality)
      @test is_move_legal(position, from_flat(i, position)) == bulk_legal
    end
    # flip the colors and check that everything is still (il)legal
    position = GoPosition(env,board = -board, to_play = WHITE)
    for move in illegal_moves
      @test !is_move_legal(position, move)
    end
    for move in legal_moves
      @test is_move_legal(position, move)
    end
    bulk_legality = all_legal_moves(position)
    for (i, bulk_legal) in enumerate(bulk_legality)
      @test is_move_legal(position, from_flat(i, position)) == bulk_legal
    end
  end

  @testset "move" begin
    start_position = GoPosition(env,
        board=TEST_BOARD,
        n = 0,
        komi = 6.5,
        caps =(1, 2),
        ko = nothing,
        recent = Vector{PlayerMove}(),
        to_play = BLACK
    )
    expected_board = load_board("""
        .XX....OO
        X........
    """ * repeat(EMPTY_ROW, 7), env)
    expected_position = GoPosition(env,
        board = expected_board,
        n = 1,
        komi = 6.5,
        caps = (1, 2),
        ko = nothing,
        recent=[PlayerMove(BLACK, from_kgs("C9", env))],
        to_play = WHITE
    )
    actual_position = play_move!(start_position, from_kgs("C9", env))
    @test assertEqualPositions(actual_position, expected_position)

    expected_board2 = load_board("""
        .XX....OO
        X.......O
    """ * repeat(EMPTY_ROW, 7), env)
    expected_position2 = GoPosition(env,
        board = expected_board2,
        n = 2,
        komi = 6.5,
        caps = (1, 2),
        ko = nothing,
        recent = [PlayerMove(BLACK, from_kgs("C9", env)), PlayerMove(WHITE, from_kgs("J8", env))],
        to_play = BLACK
    )
    actual_position2 = play_move!(actual_position, from_kgs("J8", env))
    @test assertEqualPositions(actual_position2, expected_position2)
  end

  @testset "move_with_capture" begin
    start_board = load_board(repeat(EMPTY_ROW, 5) * """
        XXXX.....
        XOOX.....
        O.OX.....
        OOXX.....
    """, env)
    start_position = GoPosition(env,
        board = start_board,
        n = 0,
        komi = 6.5,
        caps = (1, 2),
        ko = nothing,
        recent = Vector{PlayerMove}(),
        to_play = BLACK
    )
    expected_board = load_board(repeat(EMPTY_ROW, 5) * """
        XXXX.....
        X..X.....
        .X.X.....
        ..XX.....
    """, env)
    expected_position = GoPosition(env,
        board = expected_board,
        n = 1,
        komi = 6.5,
        caps = (7, 2),
        ko = nothing,
        recent = Vector{PlayerMove}([PlayerMove(BLACK, from_kgs("B2", env))]),
        to_play = WHITE
    )
    actual_position = play_move!(start_position, from_kgs("B2", env))
    @test assertEqualPositions(actual_position, expected_position)
  end

  @testset "ko_move" begin
    start_board = load_board("""
        .OX......
        OX.......
    """ * repeat(EMPTY_ROW, 7), env)
    start_position = GoPosition(env,
        board = start_board,
        n = 0,
        komi = 6.5,
        caps = (1, 2),
        ko = nothing,
        recent = Vector{PlayerMove}(),
        to_play = BLACK
    )
    expected_board = load_board("""
        X.X......
        OX.......
    """ * repeat(EMPTY_ROW, 7), env)
    expected_position = GoPosition(env,
        board = expected_board,
        n = 1,
        komi = 6.5,
        caps = (2, 2),
        ko = from_kgs("B9", env),
        recent = Vector{PlayerMove}([PlayerMove(BLACK, from_kgs("A9", env))]),
        to_play = WHITE
    )
    actual_position = play_move!(start_position, from_kgs("A9", env))
    @test assertEqualPositions(actual_position, expected_position)

    # Check that retaking ko is illegal until two intervening moves
    @test_throws IllegalMove play_move!(actual_position, from_kgs("B9", env))
    pass_twice = pass_move!(pass_move!(actual_position))
    ko_delayed_retake = play_move!(pass_twice, from_kgs("B9", env))
    expected_position = GoPosition(env,
        board = start_board,
        n = 4,
        komi = 6.5,
        caps = (2, 3),
        ko = from_kgs("A9", env),
        recent = Vector{PlayerMove}([
            PlayerMove(BLACK, from_kgs("A9", env)),
            PlayerMove(WHITE, nothing),
            PlayerMove(BLACK, nothing),
            PlayerMove(WHITE, from_kgs("B9", env))]),
        to_play = BLACK
        )

    @test assertEqualPositions(ko_delayed_retake, expected_position)
  end

  @testset "is_game_over" begin
    root = GoPosition(env)
    @test !root.done
    first_pass = play_move!(root, nothing)
    @test !first_pass.done
    second_pass = play_move!(first_pass, nothing)
    @test second_pass.done
  end

  @testset "scoring" begin
    board = load_board("""
        .XX......
        OOXX.....
        OOOX...X.
        OXX......
        OOXXXXXX.
        OOOXOXOXX
        .O.OOXOOX
        .O.O.OOXX
        ......OOO
    """, env)
    position = GoPosition(env,
        board = board,
        n = 54,
        komi = 6.5,
        caps = (2, 5),
        ko = nothing,
        recent = Vector{PlayerMove}(),
        to_play = BLACK
    )
    expected_score = 1.5
    @test score(position) == expected_score

    board = load_board("""
        XXX......
        OOXX.....
        OOOX...X.
        OXX......
        OOXXXXXX.
        OOOXOXOXX
        .O.OOXOOX
        .O.O.OOXX
        ......OOO
      """, env)
    position = GoPosition(env,
        board = board,
        n = 55,
        komi = 6.5,
        caps = (2, 5),
        ko = nothing,
        recent = Vector{PlayerMove}(),
        to_play = WHITE
    )
    expected_score = 2.5
    @test score(position) == expected_score
  end
end


text:
const DATEFORMATSTR = "d u y H:M:S"
const RHYTHM_PROC = "Sources/Rhythm FPGA"

### Types ###
"""
    OEChannel{T<:AbstractString}
Type for continuous recording channel metadata

# Fields

**`name`** `T` of channel name

**`number`** `Int` of channel number in GUI

**`bitvolts`** `Float64` of volts per ADC bit

**`position`** `Vector{Int}` vector of start position in the file for each recording's data.

**`filename`** `T` name of associated `.continuous` file
"""
struct OEChannel{T<:AbstractString}
    name::T
    number::Int
    bitvolts::Float64
    position::Vector{Int}
    filename::T
end

"""
    OEProcessor{T<:AbstractString}
Abstract type for recording Open Ephys processors.
"""
abstract type OEProcessor{T<:AbstractString} end

"""
    OERhythmProcessor{T<:AbstractString}(proc_e::LightXML.XMLElement)
Type for Rhythm processor metadata, subtype of [`OEProcessor`](@ref).

Construct with XML element for processor.

# Fields

**`id`** `Int` of processor ID in GUI

**`lowcut`** `Float64` of low pass filter cutoff on headstages

**`highcut`** `Float64` of high pass filter cutoff on headstages

**`adcs_on`** `Bool` `true` if ADCs on

**`noiseslicer`** `Bool` `true` if noiseslicer used for ADC

**`ttl_fastsettle`** `Bool` `true` if TTL fast settle used

**`dac_ttl`** `Bool` `true` if dac ttl is on

**`dac_hpf`** `Bool` `true` if dac hpf is on

**`dsp_offset`** `Bool` `true` if headstage DSP offset removal is used

**`dsp_cutoff`** `Float64` of DSP high pass filter cutoff

**`channels`** `Vector{OEChannel{T}}` list of [`OEChannel`](@ref) in
Rhythm processor
"""
struct OERhythmProcessor{T<:AbstractString} <: OEProcessor{T}
    id::Int
    lowcut::Float64
    highcut::Float64
    adcs_on::Bool
    noiseslicer::Bool
    ttl_fastsettle::Bool
    dac_ttl::Bool
    dac_hpf::Bool
    dsp_offset::Bool
    dsp_cutoff::Float64
    channels::Vector{OEChannel{T}}
end
function OERhythmProcessor(proc_e::LightXML.XMLElement)
    id_attr = attribute(proc_e, "NodeId", required = true)
    id = parse(Int, id_attr)

    channels = channel_arr(proc_e)

    # Editor
    editor_e = required_find_element(proc_e, "EDITOR")
    ## Get attribute strings
    lowcut_attr = attribute(editor_e, "LowCut", required=true)
    highcut_attr = attribute(editor_e, "HighCut", required=true)
    adcs_on_attr = attribute(editor_e, "ADCsOn", required=true)
    noiseslicer_attr = attribute(editor_e, "NoiseSlicer", required=true)
    ttl_fastsettle_attr = attribute(editor_e, "TTLFastSettle", required=true)
    dac_ttl_attr = attribute(editor_e, "DAC_TTL", required=true)
    dac_hpf_attr = attribute(editor_e, "DAC_HPF", required=true)
    dsp_offset_attr = attribute(editor_e, "DSPOffset", required=true)
    dsp_cutoff_attr = attribute(editor_e, "DSPCutoffFreq", required=true)
    ## parse attribute strings
    lowcut = parse(Float64, lowcut_attr)
    highcut = parse(Float64, highcut_attr)
    adcs_on = parse(Int, adcs_on_attr) == 1
    noiseslicer = parse(Int, noiseslicer_attr) == 1
    ttl_fastsettle = parse(Int, ttl_fastsettle_attr) == 1
    dac_ttl = parse(Int, dac_ttl_attr) == 1
    dac_hpf = parse(Int, dac_hpf_attr) == 1
    dsp_offset = parse(Int, dsp_offset_attr) == 1
    dsp_cutoff = parse(Float64, dsp_cutoff_attr)

    return OERhythmProcessor(
        id,
        lowcut,
        highcut,
        adcs_on,
        noiseslicer,
        ttl_fastsettle,
        dac_ttl,
        dac_hpf,
        dsp_offset,
        dsp_cutoff,
        channels
    )
end

"""
    TreeNode{T}
Abstract node type for tree structure, with type `T` content.

Subtypes must have the following fields:

# Required Fields

**`content`** `T` content of node.

**`parent`** `Int` ID of parent node

**`children`** `Vector{Int}` IDs of children node
"""
abstract type TreeNode{T} end

"""
    SignalNode{T<:OEProcessor}
Node type for OEProcessor signal chain, subtype of [`TreeNode`](@ref).

See [`TreeNode`](@ref) for information on fields.
"""
mutable struct SignalNode{T<:OEProcessor} <: TreeNode{T}
    content::T
    parent::Int
    children::Vector{Int}
end

"""
    Tree{T}
Abstract type for tree structure, with type `T` content.

Contains a group of [`TreeNode`](@ref) in the single required field:

# Required Fields

**`nodes`** Indexable list of [`TreeNode`](@ref) elements.
"""
abstract type Tree{T} end

"""
    OESignalTree{T<:OEProcessor}(chain_e::LightXML.XMLElement, [recording_names::Set])

Signal tree for recording processors. Since [`OpenEphysLoader`](@ref)
currently on works on `.continuous` file types, this will search for the
first [`OERhythmProcessor`](@ref) and make a signal tree up to that point.

Construct with a XML signalchain element, and a set of processor names that
are valid recording nodes.

See [`Tree`](@ref) for field information.
"""
struct OESignalTree{T<:OEProcessor} <: Tree{T}
    nodes::Vector{SignalNode{T}}
end
function OESignalTree(
    chain_e::LightXML.XMLElement;
    recording_names = Set([RHYTHM_PROC])
    )
    childs = child_elements(chain_e)
    for ch_e in childs # Break at first recording processor!
        if name(ch_e) == "PROCESSOR"
            procname = attribute(ch_e, "name", required = true)
            if procname in recording_names
                if procname == RHYTHM_PROC
                    proc = OERhythmProcessor(ch_e)
                    signode = SignalNode(proc, 0, Vector{Int}())
                    return OESignalTree([signode])
                end
            end
        end
    end
end

"""
    OEInfo{T<:AbstractString}(info_e::LigthXML.XMLElement)

Type to represent the info element in `settings.xml` made by Open Ephys.

Construct with the XML info element.

# Fields

**`gui_version`** `VersionNumber` GUI version

**`plugin_api_version`** `VersionNumber` plugin API version. If `gui_version`
is less than `0.4.0` then this will be `0`

**`datetime`** `DateTime` date and time that `settings.xml` was made

**`os`** `T` Operating system of computer running GUI

**`machine`** `T` hostname of computer running GUI
"""
struct OEInfo{T<:AbstractString}
    gui_version::VersionNumber
    plugin_api_version::VersionNumber
    datetime::DateTime
    os::T
    machine::T
end
function OEInfo(info_e::LightXML.XMLElement)
    ver_e = required_find_element(info_e, "VERSION")
    ver_str = content(ver_e)
    ver_sections = split(ver_str, '.')
    if length(ver_sections) > 3
        gui_version = VersionNumber(join(ver_sections[1:3], '.'))
    else
        gui_version = VersionNumber(ver_str)
    end
    if gui_version > v"0.4.0"
        plugin_e = required_find_element(info_e, "PLUGIN_API_VERSION")
        plugin_api = VersionNumber(content(plugin_e))
    else
        plugin_api = v"0"
    end
    date_e = required_find_element(info_e, "DATE")
    datetime = DateTime(content(date_e), DATEFORMATSTR)
    os_e = required_find_element(info_e, "OS")
    os = content(os_e)
    machine_e = required_find_element(info_e, "MACHINE")
    machine = content(machine_e)
    return OEInfo(gui_version,
                  plugin_api,
                  datetime,
                  os,
                  machine)
end

"""
    OESettings{S<:AbstractString, T<:OEProcessor}(xdoc::LightXML.XMLDocument)
Type to represent information in the `settings.xml` file made by the Open Ephys GUI.

Construct with the XML document for `settings.xml`

# Fields

**`info`** [`OEInfo`](@ref) GUI info.

**`recording_chain`** [`OESignalTree`](@ref) Signal tree that leads to recording
processors.
"""
struct OESettings{S<:AbstractString, T<:OEProcessor}
    info::OEInfo{S}
    # only processors that lead to recording nodes
    recording_chain::OESignalTree{T}
end
function OESettings(xdoc::LightXML.XMLDocument)
    settings_e = root(xdoc)
    @assert name(settings_e) == "SETTINGS" "Not a settings xml"
    info_e = required_find_element(settings_e, "INFO")
    oe_info = OEInfo(info_e)
    chain_e = required_find_element(settings_e, "SIGNALCHAIN")
    signaltree = OESignalTree(chain_e)
    return OESettings(oe_info, signaltree)
end

"""
    OERecordingMeta{T<:OEProcessor}(settings::OESettings, rec_e::LightXML.XMLElement)
Type that represents recording metadata in `Continuous_Data.openephys` file
made by the Open Ephys GUI.

Construct with a [`OESettings`](@ref) from the `settings.xml` file, and the XML
recording element of the `Continuous_Data.openephys` file.

# Fields

**`number`** `Int` Recording number

**`samplerate`** `Float64` Sampling rate

**`recording_processors`** `Vector{T}` list of recording processors
"""
struct OERecordingMeta{T<:OEProcessor}
    number::Int
    samplerate::Float64
    recording_processors::Vector{T}
end
function OERecordingMeta(settings::OESettings{S, T}, rec_e::LightXML.XMLElement) where {S, T}
    no_attr = attribute(rec_e, "number", required=true)
    number = parse(Int, no_attr)
    samprate_attr = attribute(rec_e, "samplerate", required=true)
    samplerate = parse(Float64, samprate_attr)

    proc_es = get_elements_by_tagname(rec_e, "PROCESSOR")
    if isempty(proc_es)
        throw(CorruptedException("Could not find PROCESSOR elements"))
    end
    nproc = length(proc_es)
    @compat rec_procs = Vector{T}(undef, nproc)

    for (i, proc_e) in enumerate(proc_es)
        id = find_matching_proc(settings.recording_chain, proc_e)
        rec_procs[i] = settings.recording_chain.nodes[id].content
    end
    return OERecordingMeta(number, samplerate, rec_procs)
end

"""
    OEExperMeta{S<:AbstractString, T<:OEProcessor}(s::OESettings, exper::LightXML.XMLElement)
Type to represent the Experiment metadata in `Continuous_Data.openephys`.

Construct with the [`OESettings`](@ref) from `settings.xml` and XML experiment element.

# Fields

**`file_version`** `VersionNumber` continuous file format version

**`experiment_number`** `Int` experiment number

**`separate_files`** `Bool` `true` if files are separate

**`recordings`** `Vector{OERecordingMeta{T}}` `Vector` of each [`OERecordingMeta`](@ref) within the experiment

**`settings`** [`OESettings`](@ref) of the `settings.xml` file
"""
struct OEExperMeta{S<:AbstractString, T<:OEProcessor}
    file_version::VersionNumber
    experiment_number::Int
    separate_files::Bool
    recordings::Vector{OERecordingMeta{T}}
    settings::OESettings{S, T}
end
function OEExperMeta(settings::OESettings{S, T}, exper_e::LightXML.XMLElement) where {S, T}
    # Get version
    ver_attr = attribute(exper_e, "version", required=true)
    ## Version numbers in OE files are floats, often see 0.400000000000002
    ## Need to extract the "meaningful" part
    ver_reg = r"^(\d+)\.([1-9]\d*?)0*?\d?$"
    m = match(ver_reg, ver_attr)
    ver_str = string(m.captures[1], '.', m.captures[2])
    ver = VersionNumber(ver_str)

    # Get number
    no_attr = attribute(exper_e, "number", required=true)
    exper_no = parse(Int, no_attr)

    # Get separatefiles
    separate_attr = attribute(exper_e, "separatefiles", required=true)
    separate_files = parse(Int, separate_attr) == 1

    # get recordings
    rec_es = get_elements_by_tagname(exper_e, "RECORDING")
    if isempty(rec_es)
        throw(CorruptedException("Could not find RECORDING elements"))
    end
    nrec = length(rec_es)
    @compat recordings = Vector{OERecordingMeta{T}}(undef, nrec)
    for (i, rec_e) in enumerate(rec_es)
        recordings[i] = OERecordingMeta(settings, rec_e)
    end

    return OEExperMeta(ver,
                       exper_no,
                       separate_files,
                       recordings,
                       settings)
end

### Top level functions ###
"""
    metadata([dirpath::AbstractString = pwd()]; settingsfile = "settings.xml", continuousmeta="Continuous_Data.openephys")

Top-level function to read a directory and parse the `settings.xml` and `Continuous_data.openeephys` files.

returns a [`OEExperMeta`](@ref).
"""
function metadata(
    dirpath::AbstractString = pwd();
    settingsfile::AbstractString = "settings.xml",
    continuousmeta::AbstractString = "Continuous_Data.openephys"
)
    settingspath = joinpath(dirpath, settingsfile)
    continuouspath = joinpath(dirpath, continuousmeta)
    return metadata(settingspath, continuouspath)
end
function metadata(settingspath::AbstractString, continuouspath::AbstractString)
    isfile(settingspath) || error("$settingspath does not exist")
    isfile(continuouspath) || error("$continuouspath does not exist")
    local exper_meta, settings
    parse_file(settingspath) do settingsdoc
        settings = OESettings(settingsdoc)
        parse_file(continuouspath) do contdoc
            exper_e = root(contdoc)
            add_continuous_meta!(settings, exper_e)
            exper_meta = OEExperMeta(settings, exper_e)
        end
    end
    return exper_meta
end

### Open Ephys XML parsing functions ###
"""
Parse XML Element PROCESSOR and recover channel metadata.
"""
function channel_arr(proc_e::LightXML.XMLElement, ::Type{T} = String) where {T<:AbstractString}
    # Assumes that channels are sorted!
    # Channels
    channel_vec = get_elements_by_tagname(proc_e, "CHANNEL")
    if isempty(channel_vec)
        throw(CorruptedException("Could not find CHANNEL elements"))
    end
    nchan = length(channel_vec)
    chan_rec = fill(false, nchan)
    @compat chnos = Array{Int}(undef, nchan)
    for (i, chan_e) in enumerate(channel_vec)
        sel_e = required_find_element(chan_e, "SELECTIONSTATE")
        record_attr = attribute(sel_e, "record", required=true)
        record = parse(Int, record_attr)
        if record == 1
            chan_rec[i] = true
            chno_attr = attribute(chan_e, "number", required=true)
            chno = parse(Int, chno_attr)
            chnos[i] = chno
        end
    end
    nrec = sum(chan_rec)

    # Channel info
    ch_info_e = required_find_element(proc_e, "CHANNEL_INFO")
    chinfo_children = collect(child_elements(ch_info_e))
    @compat channels = Array{OEChannel{T}}(undef, nrec)
    recno = 1
    for (i, chan_e) in enumerate(chinfo_children)
        if chan_rec[i]
            info_chno_attr = attribute(chan_e, "number", required = true)
            info_chno = parse(Int, info_chno_attr)
            @assert info_chno == chnos[i] "Channels not in same order"
            chname = attribute(chan_e, "name", required = true)
            bitvolt_attr = attribute(chan_e, "gain", required = true)
            bitvolts = parse(Float64, bitvolt_attr)
            @compat channels[recno] = OEChannel{String}(chname,
                                                info_chno,
                                                bitvolts,
                                                Vector{Int}(),
                                                "")
            recno += 1
        end
    end

    return channels
end

"""
Add data from `Continuous_Data.openephys` to [`OESettings`](@ref) from `settings.xml`
"""
function add_continuous_meta!(settings::OESettings, exper_e::LightXML.XMLElement)
    rec_es = get_elements_by_tagname(exper_e, "RECORDING")
    if ! recordings_are_consistent(rec_es)
        throw(CorruptedException("Channel information not consistent across recordings"))
    end
    for rec_e in rec_es
        proc_es = get_elements_by_tagname(rec_e, "PROCESSOR")
        if isempty(proc_es)
            throw(CorruptedException("Could not find PROCESSOR elements"))
        end
        for proc_e in proc_es
            id = find_matching_proc(settings.recording_chain, proc_e)
            add_continuous_meta!(settings.recording_chain.nodes[id].content, proc_e)
        end
    end
end
function add_continuous_meta!(proc::OERhythmProcessor, proc_e::LightXML.XMLElement)
    chan_es = get_elements_by_tagname(proc_e, "CHANNEL")
    isempty(chan_es) && throw(CorruptedException("Could not find CHANNEL elements"))
    for (i, chan_e) in enumerate(chan_es)
        name = attribute(chan_e, "name", required=true)
        if name != proc.channels[i].name
            throw(CorruptedException("Channel names don't match!"))
        end
        filename = attribute(chan_e, "filename", required=true)
        position_attr = attribute(chan_e, "position", required=true)
        file_position = parse(Int, position_attr)
        if isempty(proc.channels[i].filename)
            proc.channels[i] = OEChannel(
                proc.channels[i].name,
                proc.channels[i].number,
                proc.channels[i].bitvolts,
                [file_position],
                filename
            )
        else
            push!(proc.channels[i].position, file_position)
        end
    end
end

mutable struct XmlNode
    node_tag::String
    daughter_el::Vector{XmlNode}
    attr::Vector{String}
end

function recordings_are_consistent(rec_es::Vector{LightXML.XMLElement})
    chan_xml = XmlNode("CHANNEL", XmlNode[], ["name", "filename"])
    proc_xml = XmlNode("PROCESSOR", [chan_xml], ["id"])
    nrec = length(rec_es)
    @compat attr_sets = Vector{Set{String}}(undef, nrec)
    for (i, rec_e) in enumerate(rec_es)
        attr_sets[i] = Set(recurse_xml_attr(rec_e, proc_xml))
    end
    return all(x -> x == attr_sets[1], view(attr_sets, 2:nrec))
end

function recurse_xml_attr(e::LightXML.XMLElement, node_tree::XmlNode)
    outs = Vector{String}()
    recurse_xml_attr!(e, node_tree, outs, "")
    return outs
end

function recurse_xml_attr!(e::LightXML.XMLElement, node_tree::XmlNode, outs::Vector{String}, prepend::String)
    target_es = get_elements_by_tagname(e, node_tree.node_tag)
    for target_e in target_es
        node_attrs = map(x -> attribute(target_e, x, required=true), node_tree.attr)
        curr_path = string(prepend, join(node_attrs, "/"), "/")
        if isempty(node_tree.daughter_el)
            push!(outs, curr_path)
        else
            for daughter in node_tree.daughter_el
                recurse_xml_attr!(target_e, daughter, outs, curr_path)
            end
        end
    end
end

"""
Find id of processor in [`OESignalTree`](@ref) that matches id of XML processor element
"""
function find_matching_proc(
    sig_chain::OESignalTree,
    proc_e::LightXML.XMLElement
)
    proc_id_attr = attribute(proc_e, "id", required=true)
    proc_id = parse(Int, proc_id_attr)
    pred = x::OERhythmProcessor -> x.id == proc_id
    maybe_match = find_by(pred, sig_chain)
    if maybe_match < 0
        throw(CorruptedException("Could not find matching processor"))
    end
    return maybe_match
end

### Tree functions ###
# Taken from Scott Jones' StackOverflow answer on 4/13/2016
children(tree::Tree, id::Int) = tree.nodes[id].children
function find_by(pred::Function, tree::Tree)
    return find_by(pred, tree, 1)
end
function find_by(pred::Function, tree::Tree, id::Int)
    if pred(tree.nodes[id].content)
        return id
    else
        for child in children(tree, id)
            maybe_match = find_by(pred, tree, child)
            if maybe_match > 0
                return maybe_match
            end
        end
    end
    return -1 # Didn't find a match
end

### Helper Functions ###
function required_find_element(e::LightXML.XMLElement, name::AbstractString)
    maybe_e = find_element(e, name)
    @compat isa(maybe_e, Nothing) && throw(CorruptedException("Could not find $name element"))
    return maybe_e
end

function parse_file(f::Function, args...)
    xdoc = parse_file(args...)
    try
        f(xdoc)
    finally
        free(xdoc)
    end
end

### Display methods ###
show(io::IO, a::OESettings) = showfields(io, a)
show(io::IO, a::OEInfo) = showfields(io, a)
show(io::IO, a::OEExperMeta) = showfields(io, a)
function show(io::IO, a::OERecordingMeta)
    compact = get(io, :compact, false)
    compact ? print(io, "recording $(a.number)") : showfields(io, a)
end
show(io::IO, a::OEProcessor) = showfields(io, a)
function show(io::IO, a::OEChannel)
    compact = get(io, :compact, false)
    compact ? print(io, a.name) : showfields(io, a)
end
show(io::IO, a::SignalNode) = show(io, a.content)

showfields(io::IO, a::Any) = showfields(IOContext(io, :depth => 0), a)
function showfields(io::IOContext, a::T) where T
    depth = get(io, :depth, 0)
    pad = "  " ^ depth
    fields = fieldnames(T)
    depth > 0 && print(io, '\n')
    next_io = IOContext(IOContext(io, :typeinfo => Any), :depth => depth + 1)
    for field in fields
        print(io, pad, "$field: ")
        show(next_io, getfield(a, field))
        print(io, '\n')
    end
end

ext: jl


text:
module MoistAirProp

import Roots
"""
    psat_water(T)

Calculates the saturation pressure of the water in ``kPa`` at a given temperature in ``Kelvin``.
`T` should range between 173.15 to 473.15.
Ref -> ASHRAE Handbook - Fundamentals - Chapter-1

"""
function psat_water(T)
    C1 = -5.6745359E+03
    C2 = 6.3925247
    C3 = -9.6778430E-03
    C4 = 6.2215701E-07
    C5 = 2.0747825E-09
    C6 = -9.4840240E-13
    C7 = 4.1635019
    C8 = -5.8002206E+03
    C9 = 1.3914993
    C10 = -4.8640239E-02
    C11 = 4.1764768E-05
    C12 = -1.4452093E-08
    C13 = 6.5459673
    if (T > 173.15 && T < 273.15)
        p_ws = exp(C1/T + C2 + C3*T + C4*T^2 + C5*T^3 + C6*T^4 + C7*log(T))
    elseif (T >= 273.15 && T < 473.15)
        p_ws = exp(C8/T + C9 + C10*T + C11*T^2 + C12*T^3 + C13*log(T))
    else
        error("Temperature should be between -100 and 200 deg C")
    end #if
    return p_ws/1000
end

"""
    humidity_ratio(p, t, str, val)

Calculates humidity ratio provided `p` total pressure in ``kPa``, `t` Temperature in kelvin, `str` string could be (``rh`` for relative humidty, ``wbt`` for wet bulb temperature, ``dpt`` for dewpoint temperature, ``w`` for humidity ratio), `val` value of variable specified in the string
"""
function humidity_ratio(p, t, str, val)
    t0 = 273.15
    if str == "rh" || str == "RH"
        rh = val
        p_ws = psat_water(t)
        w = 0.621945*rh*(p_ws)/(p-rh*p_ws)
    elseif str == "WBT" || str == "wbt"
        wbt = val
        p_ws_wbt = psat_water(wbt)
        ws_wbt = 0.621945*(p_ws_wbt)/(p-p_ws_wbt)
        if t >= t0
            w = ((2501-2.326*(wbt-t0))*ws_wbt - 1.006*((t-t0) - (wbt-t0)))/(2501 + 1.86*(t-t0) - 4.186*(wbt-t0))
        else
            w = ((2830-0.24*(wbt-t0))*ws_wbt - 1.006*((t-t0) - (wbt-t0)))/(2830 + 1.86*(t-t0) - 2.1*(wbt-t0))
        end #if
    elseif str == "dpt" || str == "DPT"
        dpt = val
        p_ws = psat_water(dpt)
        w = 0.621945*(p_ws)/(p-p_ws)
    elseif str == "w" || str == "W"
        w = val
    else
        error("use 'rh' for relative humidity,'wbt' for wet bulb temp,'dpt' for dew point temp")
    end #if
    return w
end #function

"""
    density_moistair(p,t,str,val)

Calculates density of moist air in ``kg/m^3`` provided `p` total pressure in ``kPa``, `t` Temperature in kelvin, `str` string could be (``rh`` for relative humidty, ``wbt`` for wet bulb temperature, ``dpt`` for dewpoint temperature, ``w`` for humidity ratio), `val` value of variable specified in the string
"""
function density_moistair(p,t,str,val)
    if str == "w" || str == "w"
        w = val
    else
        w = humidity_ratio(p, t, str, val)
    end
    v = 0.287042*(t)*(1 + 1.607858*w)/p
    return 1/v
end
"""
    dewpoint(p,t,str,val)

Calculates dewpoint temperature of moist air in ``kelvin`` provided `p` total pressure in ``kPa``, `t` Temperature in ``kelvin``, `str` string could be (``rh`` for relative humidty, ``wbt`` for wet bulb temperature, ``w`` for humidity ratio), `val` value of variable specified in the string
"""
function dewpoint(p,t,str,val)
    if str == "w" || str == "W"
        w = val
    else
        w = humidity_ratio(p, t, str, val)
    end
    p_w = p/(0.621945/w +1)
    t_dpt = Roots.find_zero(x -> (psat_water(x)-p_w),(t-100,t))
end
"""
    enthalpy_moistair(p,t,str,val)

Calculates enthalpy of moist air in ``kJ/kg`` provided `p` total pressure in ``kPa``, `t` Temperature in ``kelvin``, `str` string could be (``rh`` for relative humidty, ``wbt`` for wet bulb temperature, ``dpt`` for dewpoint temperature, ``w`` for humidity ratio), `val` value of variable specified in the string
"""
function enthalpy_moistair(p, t, str, val)
    if str == "W" || str == "w"
        w = val
        h = 1.006*(t-273.15)+ w*(2501+1.86*(t-273.15))
    else
        w = humidity_ratio(p, t, str, val)
        h = 1.006*(t-273.15)+ w*(2501+1.86*(t-273.15))
    end
    return h
end
"""
    wetbulb(p,t,str,val)

Calculates wetbulb temperature of moist air in ``Kelvin`` provided `p` total pressure in ``kPa``, `t` Temperature in ``kelvin``, `str` string could be (``rh`` for relative humidty, ``dpt`` for dewpoint temperature, ``w`` for humidity ratio), `val` value of variable specified in the string
"""
function wetbulb(p,t,str,val)
    w = humidity_ratio(p,t,str,val)
    if str == "dpt" || str == "DPT"
        dpt = val
    else
        dpt = dewpoint(p,t,str,val)
    end
    t_wbt = Roots.find_zero(temp -> (humidity_ratio(p,t,"wbt",temp)-w), (dpt,t))
end

export psat_water, humidity_ratio, density_moistair, dewpoint, enthalpy_moistair, wetbulb

end # module


text:
"""
Calculates the exact boundingbox of a Scene/Plot
"""
function boundingbox(x)
    data_limits(x)
end

boundingbox(x::Combined) = boundingbox(x.plots)

boundingbox(scene::Scene) = boundingbox(plots_from_camera(scene))

function boundingbox(plots::Vector)
    isempty(plots) && return FRect3D(Vec3f0(0), Vec3f0(0))
    plot_idx = iterate(plots)
    bb = FRect3D()
    while plot_idx !== nothing
        plot, idx = plot_idx
        plot_idx = iterate(plots, idx)
        bb2 = boundingbox(plot)
        isfinite(bb) || (bb = bb2)
        isfinite(bb2) || continue
        bb = union(bb, bb2)
    end
    bb
end

function boundingbox(x::Text, text::String)
    position = value(x[:position])
    @get_attribute x (textsize, font, align, rotation, model)
    boundingbox(text, position, textsize, font, align, rotation, model)
end
boundingbox(x::Text) = boundingbox(x, value(x[1]))

function boundingbox(
        text::String, position, textsize;
        font = "default", align = (:left, :bottom), rotation = 0.0,
        model::Mat4f0 = Mat4f0(I)
    )
    boundingbox(
        text, position, textsize,
        to_font(font), to_align(align), to_rotation(rotation), model
    )

end

function boundingbox(text::String, position, textsize, font, align, rotation, model)
    atlas = get_texture_atlas()
    N = length(text)
    ctext_state = iterate(text)
    ctext_state === nothing && return AABB(Vec3f0(0), Vec3f0(0))
    pos_per_char = !isa(position, VecTypes)
    start_pos = Vec(pos_per_char ? first(position) : position)
    start_pos2D = to_ndim(Point2f0, start_pos, 0.0)
    last_pos = Point2f0(0, 0)
    aoffsetn = to_ndim(Vec3f0, align, 0f0)
    start_pos3d = to_ndim(Vec3f0, start_pos, 0.0)
    bb = AABB(start_pos3d, Vec3f0(0))
    broadcast_foreach(1:N, rotation, font, textsize) do i, rotation, font, scale
        c, text_state = ctext_state
        ctext_state = iterate(text, text_state)
        if c != '\r'
            pos = if pos_per_char
                to_ndim(Vec3f0, position[i], 0.0)
            else
                last_pos = calc_position(last_pos, Point2f0(0, 0), atlas, c, font, scale)
                rotation * (start_pos3d .+ to_ndim(Vec3f0, last_pos, 0.0))
            end
            s = glyph_scale!(atlas, c, font, scale)
            srot = rotation * to_ndim(Vec3f0, s, 0.0)
            bb = GeometryTypes.update(bb, pos)
            bb = GeometryTypes.update(bb, pos .+ srot)
        end
    end
    FRect3D(minimum(bb) .- (widths(bb) .* aoffsetn), widths(bb))
end


text:
using DelimitedFiles

struct MatchingSphere
    x::String
    y::String
    z::String
    ms_id::Int64
end

function add_sphere!(ms_dict, receptor, coords)
    if !haskey(ms_dict, receptor)
        ms_dict[receptor] = Dict()
    end

    if !haskey(ms_dict[receptor], coords)
        ms_id = length(ms_dict[receptor]) + 1
        ms_dict[receptor][coords] = ms_id
    end

    sphere = MatchingSphere(
        coords[1], coords[2], coords[3],
        ms_dict[receptor][coords]
    )

    return sphere
end

function add_spheres!(ms_dict, values, spheres)
    receptor = values[1]
    coords_1 = Tuple(values[9:11])
    coords_2 = Tuple(values[12:14])
    coords_3 = Tuple(values[15:17])
    coords_4 = Tuple(values[18:20])

    for (index, coords) in enumerate([coords_1, coords_2, coords_3, coords_4])
        spheres[index] = add_sphere!(ms_dict, receptor, coords)
    end
end

function update_usage!(sph_usage, values, s)
    index = Tuple(values[1:3])
    if !haskey(sph_usage, index)
        sph_usage[index] = Dict()
    end

    if !haskey(sph_usage[index], s)
        sph_usage[index][s] = Dict(
            :Usage => 0,
            :Ligand_Usage => 0,
            :Decoy_Usage => 0
        )
    end

    sph_usage[index][s][:Usage] += 1
    if values[end] == "ligand"
        sph_usage[index][s][:Ligand_Usage] += 1
    else
        sph_usage[index][s][:Decoy_Usage] += 1
    end
end

function update_cooccurrence!(co_occurrence, values, spheres)
    index = Tuple(values[1:3])

    if !haskey(co_occurrence, index)
        co_occurrence[index] = zeros((45,45))
    end

    for s_i in spheres
        for s_j in spheres
            if s_i.ms_id != s_j.ms_id
                co_occurrence[index][s_i.ms_id,s_j.ms_id] += 1
            end
        end
    end
end

function write_usage(sph_usage)
    open("sphere_usage.tab", "w") do io
        header = [
            "receptor", "match_type", "cluster_id",
            "x", "y", "z", "ms_id",
            "Usage", "Ligand_Usage", "Decoy_Usage"
        ]
        writedlm(io, [header])
        for idx in keys(sph_usage)
            for s in keys(sph_usage[idx])
                values = [
                    idx[1], idx[2], idx[3],
                    s.x, s.y, s.z, s.ms_id,
                    sph_usage[idx][s][:Usage],
                    sph_usage[idx][s][:Ligand_Usage],
                    sph_usage[idx][s][:Decoy_Usage],
                    ]
                writedlm(io,[values])
            end
        end

    end
end

function write_cooccurrence(co_occurrence)
    open("co-occurrence.tab", "w") do io
        header = vcat(
            ["receptor", "match_type", "cluster_id"],
            ["sph.$i" for i in 1:45]
            )
        writedlm(io, [header])

        for index in keys(co_occurrence)
            for i in eachrow(co_occurrence[index])
                values = vcat([_ for _ in index], [c for c in i])
                writedlm(io, [values])
            end
        end
    end
end

function main(fn)

    ms_dict = Dict()
    sph_usage = Dict()
    co_occurrence = Dict()

    open(fn) do io

        # preallocated
        spheres = Array{MatchingSphere, 1}(undef, 4)

        while !eof(io)

            # read and split line
            line = readline(io)
            values = split(strip(line))

            # create matching spheres
            add_spheres!(ms_dict, values, spheres)

            # update specific sample usage
            for s in spheres
                update_usage!(sph_usage, values, s)
            end

            # updates sphere co-occurrence matrix
            update_cooccurrence!(co_occurrence, values, spheres)

        end
    end

    write_usage(sph_usage)
    write_cooccurrence(co_occurrence)
end

if length(ARGS) != 1
    println("Error : Provide a coordinate dataframe to process")
    exit()
end

filename = ARGS[1]
main(filename)

max_stars_repo_name: noamteyssier/ClusterDock

max_issues_repo_name: noamteyssier/ClusterDock

max_forks_repo_name: noamteyssier/ClusterDock


text:
import DelayEmbeddings: AbstractDataset

function invariantize(pts::AbstractDataset, bc = "circular")
    invariant_pts = copy(pts.data)
    if bc == "circular"
        return push!(invariant_pts, pts.data[1])
    elseif bc == "random"
        return push!(invariant_pts, pts.data[rand(1:length(pts))])
    end
    return invariant_pts
end


text:
# # How Julia works?
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__notebooks/04_HowJuliaWorks.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__notebooks/04_HowJuliaWorks.ipynb)
#
# ## Type system
#
# Until now, we have not used type annotations, but Julia has a rich type system.
# Julia is an optionally- and dynamically-typed programming language. That means
# that you can change the type of a variable...

a = 10
typeof(a)

#-

a = true
typeof(a)

# ... and that type annotations are optional. You can use type annotations to:
#  - make a program more robust (type checking)

a::Bool # type assertion

#-

a::Int

#  - optimize code by giving a hint to the compiler
#  - documment the code and use multiple dispatch

## function_name(arg::ArgumentType)::ReturnType = function_body
first_character(str::AbstractString)::Char = str[1]

#-

character = first_character("ABC")

#-

character = first_character(10_000)

# ### Type hierarchy
#
# Any Julia object has a type that belongs to a fully connected type graph.
# There are abstract and concrete types. Concrete types are final, i.e. they
# cannot have subtypes, while abstract types can have multiple subtypes but
# only one supertype.

using JuliaForBioinformatics
show_type_tree(Number)

# In Julia, all values are instances of the abstract type `Any`.
#
# The functions supertype and subtypes are useful to navegate the type graph.

supertype(Real)

#-

subtypes(Real)

# You can use `isa` to test if an object is of a given type

isa("I'm a string", String)

# And the subtype operator `<:` to test if a type is a subtype of another

String <: AbstractString

# You can also use `Union` of types, for example, if the possible types don't
# share a meaningful supertype

String <: Union{AbstractString, AbstractChar}

# ### Multiple dispatch
#
# We can define multiple [methods](https://docs.julialang.org/en/v1/manual/methods/#Methods-1)
# for a function by using different method signatures by indicating the argument
# types using `::` or `<:`.
#
# For example we are going to define 3 methods for the function `say_my_type`:

say_my_type(x) = println(x, " is a ", typeof(x))
## say_my_type(x) is the same that say_my_type(x::Any)

say_my_type(x::Real) = println(x, " is a Real number of type ",  typeof(x))
say_my_type(x::Float64) = println(x, " is a Float64 number")

# When the function is called, Julia selects the method with the most specific
# method signature.

say_my_type('A') # 'A' is a Char, a subtype of Any
say_my_type(2) # 2 is an Int, a subtype of Real
say_my_type(2.0)

# `say_my_type(x::Real)` can also be written using the `where` keyword as

say_my_type(x::T) where {T <: Real} = println(x, " is a Real number of type ",  T)

#-

methods(say_my_type)

# ### Parametric types
#
# Julia types can have parameters. We have already used parametric types, one
# of them is Array:

three_d_array = zeros(Int, 4, 3, 2)

#-

typeof(three_d_array)

# Julia Arrays take two parameters, the type of the elements stored in the array
# and the array dimensions.
#
# This allows to write specific methods depending on those parameters

say_my_type(x::Array{T, 1}) where {T} = println(x, " is vector with ",  T, " elements")
say_my_type(x::Array{T, 2}) where {T} = println(x, " is matrix with ",  T, " elements")

#-

say_my_type(Rational[0.5, 1, 1.5])
say_my_type(Float64[1 3 5; 2 4 6])

# #### Exercise 1
#
# Add a method to `say_my_type` that prints the number of unique values of an
# array of characters or strings and its dimensions. Hint: You can use the
# `unique` function.

## ... println(x, " is a text array with ... dimensions and ... unique values of type ...

# That means that the function call:
# ```julia
# test_array = ['a', 'b', 'b']
# say_my_type(test_array)
# ```
# Should print something like:
# ```
# ['a', 'b', 'b'] is a text array with 1 dimensions and 2 unique values of type Char
# ```

# ### Which method is being used?
#
# You can use the `@which` macro to ask Julia which method is being used

@which say_my_type(2 + 0im)

#-

@which say_my_type(2.0)

# ### Julia compiler
#
# Julia uses Just-in-time (JIT) compilation to achieve close to C performance.
# After selecting the most specific method, Julia (generally) compiles the
# method for the particular argument types.
#
# For this reason, the first time a function is called, it is compiled (slow).
# If you call the same function a second time with the same argument types, it
# will use the already compiled code (fast).

@time sum(1:10_000)

#-

@time sum(1:10_000)

# While compilation times can be annoying sometimes, this mechanism allows
# Julia generality, composability and its capacity to generate efficient code
# for user-defined types. That's mean that you do not need to use built-in
# types or functions or to code some parts in C/Fortran to get a good
# performance like in other high-level languages.
#
#
# Also, Julia represents its own code as a Julia data structure. This allows a
# program to transform and generate its own code, using **macros** and
# **generated functions**, for example, and powerful reflection capabilities
# to explore the internals of a program. You can read the
# [metaprogramming section of the manual](https://docs.julialang.org/en/v1/manual/metaprogramming/)
# to learn more about this topic.
#

##?@elapsed

#-

macroexpand(Main, :(@elapsed sum(1:10_000)))

#
# ![](JuliaCompiler.png)
#

@which 1.0 + 2.0

#-

@which 1 + 2

#-

@code_lowered 1.0 + 2.0

#-

@code_typed 1.0 + 2.0

#-

@code_llvm 1.0 + 2.0

#-

@code_native 1.0 + 2.0

# Julia can do extra optimizations inside functions, e.g. constant propagation:

f() = 1.0 + 2.0
@code_typed f()

# #### Exercise 2
#
# Modify the `sum_numbers` function to make it
# [type stable](https://docs.julialang.org/en/v1/manual/performance-tips/#Write-%22type-stable%22-functions-1)
# by using `zero` and `eltype`.
#
# To avoid performance issues for using a
# [global variable](https://docs.julialang.org/en/v1/manual/performance-tips/index.html#Avoid-global-variables-1)
# we are going to define it as a constant using the
# [const keyword.](https://docs.julialang.org/en/v1/base/base/#const)

const rand_matrix = rand(4, 4)

#-

function sum_numbers(numbers)
    total = 0
    for value in numbers
        total += value
    end
    total
end

# Check that you don't have any warning in the output of `@code_warntype`:

@code_warntype sum_numbers(rand_matrix)

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
module Nonlinear
import Luna
import Luna.PhysData: ε_0, e_ratio
import Luna: Maths, Utils
import FFTW
import LinearAlgebra: mul!, ldiv!

function KerrScalar!(out, E, fac)
    @. out += fac*E^3
end

function KerrVector!(out, E, fac)
    for i = 1:size(E,1)
        Ex = E[i,1]
        Ey = E[i,2]
        Ex2 = Ex^2
        Ey2 = Ey^2
        out[i,1] += fac*(Ex2 + Ey2)*Ex
        out[i,2] += fac*(Ex2 + Ey2)*Ey
    end
end

"Kerr response for real field"
function Kerr_field(γ3)
    Kerr = let γ3 = γ3
        function Kerr(out, E, ρ)
            if size(E,2) == 1
                KerrScalar!(out, E, ρ*ε_0*γ3)
            else
                KerrVector!(out, E, ρ*ε_0*γ3)
            end
        end
    end
end

"Kerr response for real field but without THG"
function Kerr_field_nothg(γ3, n)
    E = Array{Float64}(undef, n)
    hilbert = Maths.plan_hilbert(E)
    Kerr = let γ3 = γ3, hilbert = hilbert
        function Kerr(out, E, ρ)
            out .+= ρ*3/4*ε_0*γ3.*abs2.(hilbert(E)).*E
        end
    end
end

function KerrScalarEnv!(out, E, fac)
    @. out += 3/4*fac*abs2(E)*E
end

function KerrVectorEnv!(out, E, fac)
    for i = 1:size(E,1)
        Ex = E[i,1]
        Ey = E[i,2]
        Ex2 = abs2(Ex)
        Ey2 = abs2(Ey)
        out[i,1] += 3/4*fac*((Ex2 + 2/3*Ey2)*Ex + 1/3*conj(Ex)*Ey^2)
        out[i,2] += 3/4*fac*((Ey2 + 2/3*Ex2)*Ey + 1/3*conj(Ey)*Ex^2)
    end
end

"Kerr response for envelope"
function Kerr_env(γ3)
    Kerr = let γ3 = γ3
        function Kerr(out, E, ρ)
            if size(E,2) == 1
                KerrScalarEnv!(out, E, ρ*ε_0*γ3)
            else
                KerrVectorEnv!(out, E, ρ*ε_0*γ3)
            end
        end
    end
end

"Kerr response for envelope but with THG"
# see Eq. 4, Genty et al., Opt. Express 15 5382 (2007)
function Kerr_env_thg(γ3, ω0, t)
    C = exp.(2im*ω0.*t)
    Kerr = let γ3 = γ3, C = C
        function Kerr(out, E, ρ)
            @. out += ρ*ε_0*γ3/4*(3*abs2(E) + C*E^2)*E
        end
    end
end

"Response type for cumtrapz-based plasma polarisation, adapted from:
M. Geissler, G. Tempea, A. Scrinzi, M. Schnürer, F. Krausz, and T. Brabec, Physical Review Letters 83, 2930 (1999)."
struct PlasmaCumtrapz{R, EType, tType}
    ratefunc::R # the ionization rate function
    ionpot::Float64 # the ionization potential (for calculation of ionization loss)
    rate::tType # buffer to hold the rate
    fraction::tType # buffer to hold the ionization fraction
    phase::EType # buffer to hold the plasma induced (mostly) phase modulation
    J::EType # buffer to hold the plasma current
    P::EType # buffer to hold the plasma polarisation
    δt::Float64 # the time step
end

"""
    PlasmaCumtrapz(t, E, ratefunc, ionpot)

Construct the Plasma polarisation response for a field on time grid `t`
with example electric field like `E`, an ionization rate callable
`ratefunc` and ionization potential `ionpot`.
"""
function PlasmaCumtrapz(t, E, ratefunc, ionpot)
    rate = similar(t)
    fraction = similar(t)
    phase = similar(E)
    J = similar(E)
    P = similar(E)
    return PlasmaCumtrapz(ratefunc, ionpot, rate, fraction, phase, J, P, t[2]-t[1])
end

"The plasma response for a scalar electric field"
function PlasmaScalar!(Plas::PlasmaCumtrapz, E)
    Plas.ratefunc(Plas.rate, E)
    Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt)
    @. Plas.fraction = 1-exp(-Plas.fraction)
    @. Plas.phase = Plas.fraction * e_ratio * E
    Maths.cumtrapz!(Plas.J, Plas.phase, Plas.δt)
    for ii in eachindex(E)
        if abs(E[ii]) > 0
            Plas.J[ii] += Plas.ionpot * Plas.rate[ii] * (1-Plas.fraction[ii])/E[ii]
        end
    end
    Maths.cumtrapz!(Plas.P, Plas.J, Plas.δt)
end

"""
The plasma response for a vector electric field.

We take the magnitude of the electric field to calculate the ionization
rate and fraction, and then solve the plasma polarisation component-wise
for the vector field.

A similar approach was used in: C Tailliez et al 2020 New J. Phys. 22 103038.
"""
function PlasmaVector!(Plas::PlasmaCumtrapz, E)
    Ex = E[:,1]
    Ey = E[:,2]
    Em = @. hypot.(Ex, Ey)
    Plas.ratefunc(Plas.rate, Em)
    Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt)
    @. Plas.fraction = 1-exp(-Plas.fraction)
    @. Plas.phase = Plas.fraction * e_ratio * E
    Maths.cumtrapz!(Plas.J, Plas.phase, Plas.δt)
    for ii in eachindex(Em)
        if abs(Em[ii]) > 0
            pre = Plas.ionpot * Plas.rate[ii] * (1-Plas.fraction[ii])/Em[ii]^2
            Plas.J[ii,1] += pre*Ex[ii]
            Plas.J[ii,2] += pre*Ey[ii]
        end
    end
    Maths.cumtrapz!(Plas.P, Plas.J, Plas.δt)
end

"Handle plasma polarisation routing to `PlasmaVector` or `PlasmaScalar`."
function (Plas::PlasmaCumtrapz)(out, Et, ρ)
    if ndims(Et) > 1
        if size(Et, 2) == 1 # handle scalar case but within modal simulation
            PlasmaScalar!(Plas, reshape(Et, size(Et,1)))
            out .+= ρ .* reshape(Plas.P, size(Et))
        else
            PlasmaVector!(Plas, Et) # vector case
            out .+= ρ .* Plas.P
        end
    else
        PlasmaScalar!(Plas, Et) # straight scalar case
        out .+= ρ .* Plas.P
    end
end

"Raman polarisation response type"
abstract type RamanPolar end

"Raman polarisation response type for a carrier resolved field"
struct RamanPolarField{TR, Tt, Thv, Tω, Tv, FTt, HTt} <: RamanPolar
    r::TR # Raman response
    h::Tt # doubled buffer to hold response + padding
    ht::Thv # buffer to hold time domain response
    hω::Tω # the frequency domain Raman response function
    Eω2::Tω # buffer to hold the Fourier transform of E^2
    Pω::Tω # buffer to hold the frequency domain polarisation
    E2::Tt # buffer to hold E^2
    E2v::Tv # view into first half of E2
    P::Tt # buffer to hold the time domain polarisation
    Pout::Tt # buffer to hold the output portion of the time domain polarisation
    FT::FTt # Fourier transform plan
    HT::HTt # Hilbert transform
    thg::Bool # do we include third harmonic generation
    dt::Float64 # time step for scaling
end

"Raman polarisation response type for an envelope"
struct RamanPolarEnv{TR, Tt, Thv, Tω, Tv, FTt} <: RamanPolar
    r::TR # Raman response
    h::Tt # doubled buffer to hold response + padding
    ht::Thv # buffer to hold time domain response
    hω::Tω # the frequency domain Raman response function
    Eω2::Tω # buffer to hold the Fourier transform of E^2
    Pω::Tω # buffer to hold the frequency domain polarisation
    E2::Tω # buffer to hold E^2
    E2v::Tv # view into first half of E2
    P::Tω # buffer to hold the time domain polarisation
    Pout::Tω # buffer to hold the output portion of the time domain polarisation
    FT::FTt # Fourier transform plan
    dt::Float64 # time step for scaling
end

"""
    RamanPolarField(t, ht; thg=true)

Construct Raman polarisation response for a field on time grid `t`
using response function `r`. If `thg=false` then exclude the third
harmonic generation component of the response.
"""
function RamanPolarField(t, r; thg=true)
    h = zeros(length(t)*2) # note double grid size, see explanation below
    ht = view(h, 1:length(t))
    Utils.loadFFTwisdom()
    FT = FFTW.plan_rfft(h, 1, flags=Luna.settings["fftw_flag"])
    inv(FT)
    Utils.saveFFTwisdom()
    hω = FT * h
    Eω2 = similar(hω)
    Pω = similar(hω)
    E2 = similar(h)
    E2v = view(E2, 1:length(t))
    P = similar(h)
    Pout = similar(t)
    HT = Maths.plan_hilbert(Pout)
    fill!(E2, 0.0)
    RamanPolarField(r, h, ht, hω, Eω2, Pω, E2, E2v, P, Pout, FT, HT, thg, t[2] - t[1])
end

"""
    RamanPolarEnv(t, ht)

Construct Raman polarisation response for an envelope on time grid `t`
using response function `r`.
"""
function RamanPolarEnv(t, r)
    h = zeros(length(t)*2) # note double grid size, see explanation below
    ht = view(h, 1:length(t))
    Utils.loadFFTwisdom()
    FT = FFTW.plan_fft(h, 1, flags=Luna.settings["fftw_flag"])
    inv(FT)
    Utils.saveFFTwisdom()
    hω = FT * h
    Eω2 = similar(hω)
    Pω = similar(hω)
    E2 = similar(hω)
    P = similar(hω)
    Pout = Array{ComplexF64,}(undef,size(t))
    E2v = view(E2, 1:length(t))
    fill!(E2, 0.0)
    RamanPolarEnv(r, h, ht, hω, Eω2, Pω, E2, E2v, P, Pout, FT, t[2] - t[1])
end

"Square the field or envelope"
function sqr!(R::RamanPolarField, E)
    if !R.thg
        # see documentation for factor of 1/2 here
        R.E2v .= 1/2 .* abs2.(R.HT(E))
    else
        R.E2v .= E.^2
    end
end

function sqr!(R::RamanPolarEnv, E)
    # see documentation for factor of 1/2 here
    R.E2v .= 1/2 .* abs2.(E)
end

"Calculate Raman polarisation for field/envelope Et"
function (R::RamanPolar)(out, Et, ρ)
    # get the field as a 1D Array
    n = size(Et, 1)
    if ndims(Et) > 1
        if size(Et, 2) == 1 # handle scalar case but within modal simulation
            E = reshape(Et, n)
        else
            # handle vector case
            error("vector Raman not yet implemented")
        end
    else
        E = Et # handle straight scalar case
    end

    # square the field or envelope in first half
    # corresponding to the field/envelope grid size
    sqr!(R, E)

    # update frequency domain response function `hω`.
    # we fill only up to the first half of h (using the view ht)
    # i.e. only the part corresponding to the original time grid
    # note that the response function time 0 is put into the first element of the response array
    # this ensures that causality is maintained, and no artificial delay between the field and
    # the start of the response function occurs, at each convolution point.
    R.r(R.ht, ρ)
    R.hω .= R.FT * R.h

    # convolution by multiplication in frequency domain
    # The double grid gives us accurate full convolution between the full field grid
    # and full response function. It is unnecessary for highly damped responses, like
    # in glass. But for gases with very long decay times it prevents artefacts due to
    # truncation of the response function. There is likely a more efficient way. But
    # this is safe, until we come up with one.
    # we scale to correct for missing dt*dt*df from IFFT(FFT*FFT)
    # the ifft already scales by 1/n = dt*df, so we need an additional dt
    R.Eω2 .= R.FT * R.E2
    @. R.Pω = R.hω * R.Eω2 * R.dt
    R.P .= R.FT \ R.Pω

    # calculate full polarisation, extracting only the valid
    # grid region, which is the first length(E) part.
    for i = 1:length(E)
        R.Pout[i] = ρ*E[i]*R.P[i]
    end

    # copy to output in dimensions requested
    if ndims(Et) > 1
        out .+= reshape(R.Pout, size(Et))
    else
        out .+= R.Pout
    end
end

end


text:
"""
    SeisPlotFK(d; <keyword arguments>)

Plot time-space, frequency-wavenumber or amplitude-frequency 2D seismic data `d`
with color, wiggles or overlay.

# Arguments
- `d::Array{Real,2}`: 2D data to plot.

# Keyword arguments
- `cmap="PuOr"`: colormap for  `"color"` or `"overlay"` style.
- `pclip=99.9`: percentile for determining clip.
- `vmin="NULL"`: minimum value used in colormapping data.
- `vmax="NULL"`: maximum value used in colormapping data.
- `aspect="auto"`: color image aspect ratio.
- `interpolation="Hanning"`: interpolation method for colormapping data.
- `fmax=100`: maximum frequency for `"FK"` or `"Amplitude"` plot.
- `title=" "`: title of plot.
- `titlesize=16`: size of title.
- `xlabel=" "`: label on x-axis.
- `xunits=" "`: units of y-axis.
- `ylabel=" "`: label on x-axis.
- `yunits=" "`: units of y-axis.
- `labelsize=14`: size of labels on axis.
- `ox=0`: first point of x-axis.
- `dx=1`: increment of x-axis.
- `oy=0`: first point of y-axis.
- `dy=1`: increment of y-axis.
- `xticks="NULL"`: ticks on x-axis.
- `yticks="NULL"`: ticks on y-axis.
- `xticklabels="NULL"`: labels on ticks of x-axis.
- `yticklabels="NULL"`: labels on ticks of y-axis.
- `ticksize=11`: size of labels on ticks.
- `fignum="NULL"`: number of figure.
- `wbox=6`: width of figure in inches.
- `hbox=6`: height of figure in inches.
- `dpi=100`: dots-per-inch of figure.
- `name="NULL"`: name of the figure to save (only if `name` is given).

# Example
```julia
julia> d = SeisLinearEvents(); SeisPlotFK(d);
```

Credits: Aaron Stanton, 2015
"""
function SeisPlotFK(d::Array{T,2};
                           cmap="PuOr", pclip=99.9, vmin="NULL", vmax="NULL",
                           aspect="auto", interpolation="Hanning", fmax=100,
                           scal="NULL", title=" ", titlesize=16, xlabel=" ", xunits=" ",
                           ylabel=" ", yunits=" ", labelsize=14, ox=0, dx=1,
                           oy=0, dy=1, xticks="NULL", yticks="NULL",
                           xticklabels="NULL", yticklabels="NULL", ticksize=11,
                           fignum="NULL", wbox=6, hbox=6, dpi=100, name="NULL") where {T<:Real}

    pl.ion()
    if (fignum == "NULL")
	fig = pl.figure(figsize=(wbox, hbox), dpi=dpi, facecolor="w",
                           edgecolor="k")
    else
	fig = pl.figure(num=fignum, figsize=(wbox, hbox), dpi=dpi,
                           facecolor="w", edgecolor="k")
    end

	xlabel = "Wavenumber"
	xunits = "(1/m)"
	ylabel = "Frequency"
	yunits = "Hz"

	dk = 1/dx/size(d, 2)
	kmin = -dk*size(d, 2)/2
	kmax =  dk*size(d, 2)/2

	df = 1/dy/size(d, 1)
	FMAX = df*size(d, 1)/2
	if fmax > FMAX
	    fmax = FMAX
	end
	nf = convert(Int32, floor((size(d, 1)/2)*fmax/FMAX))
	D = abs.(fftshift(fft(d)))
	D = D[round(Int,end/2):round(Int,end/2)+nf, :]
	if (vmin=="NULL" || vmax=="NULL")
	    a = 0.
	    if (pclip<=100)
		b = quantile(abs.(D[:]), (pclip/100))
	    else
		b = quantile(abs.(D[:]), 1)*pclip/100
	    end
	end

    im = pl.imshow(D, cmap=cmap, vmin=a, vmax=b, extent=[kmin,kmax,fmax,0.0],
                            aspect=aspect)
    pl.title(title, fontsize=titlesize)
    pl.xlabel(join([xlabel " " xunits]), fontsize=labelsize)
    pl.ylabel(join([ylabel " " yunits]), fontsize=labelsize)
    xticks == "NULL" ? nothing : pl.xticks(xticks)
    yticks == "NULL" ? nothing : pl.yticks(yticks)
    ax = pl.gca()
    xticklabels == "NULL" ? nothing : ax.set_xticklabels(xticklabels)
    yticklabels == "NULL" ? nothing : ax.set_yticklabels(yticklabels)
    pl.setp(ax.get_xticklabels(), fontsize=ticksize)
    pl.setp(ax.get_yticklabels(), fontsize=ticksize)
    if (name == "NULL")
	pl.show()
    else
	pl.savefig(name, dpi=dpi)
	pl.close()
    end
    return im
end


text:
function remove_undef!(ExpressionArray::T) where T <: JumpExpressionMatrix
    for j in 1:size(ExpressionArray)[2]
        for i in 1:size(ExpressionArray)[1]
            !isassigned(ExpressionArray,i,j) ? ExpressionArray[i,j] = AffExpr(0.0) : continue
        end
    end
end


text:
if VERSION <= v"1.1"
    const TH_EXP = 100
    const TH_INV = 1000
else
    const TH_EXP = 5000
    const TH_INV = 100_000
end

#========== exp! log! inv! ==========#

export exp0, exp_, exp!, exp!!
export log0, log_, log!, log!!
export inv0, inv_, inv!, inv!!

"""
    exp!(A)
    exp_(A) = exp!(similar(A), A)
    exp0(A) ≈ exp.(A)

Element-wise in-place exponential, and friends.
Multi-threaded when `length(A) >= $TH_EXP`.
Will be handled by `Yeppp` or `AppleAccelerate` or `IntelVectorMath` if you load one of them,
note that `Yeppp` may well be slower.
"""
function exp! end

exp0(A) = similar(A) .= exp.(A) # maps Adjoint -> Adjoint etc

@doc @doc(exp!)
exp_(A) = exp!(similar(A), A)
exp!(A) = exp!(A, A)
exp!!(A) = exp!(A) # differs in gradient

function exp!(B, A)
    @assert size(A)==size(B)
    if length(A) < TH_EXP
        for I in eachindex(A)
            @inbounds B[I] = exp1(A[I])
        end
    else
        Threads.@threads for I in eachindex(A)
            @inbounds B[I] = exp1(A[I])
        end
    end
    B
end

"""
    log!(A)
    log_(A) ≈ log!(similar(A), A)
    log0(A) = log.(A)

Element-wise in-place natural logarithm, and friends.
Multi-threaded when `length(A) >= $TH_EXP`.
Will be handled by `Yeppp` or `AppleAccelerate` or `IntelVectorMath` if you load one of them.
"""
function log! end

log0(A) = similar(A) .= log.(A)

@doc @doc(log!)
log_(A) = log!(similar(A), A)
log!(A) = log!(A, A)
log!!(A) = log!(A) # differs in gradient

function log!(B, A)
    @assert size(A)==size(B)
    if length(A) < TH_EXP
        for I in eachindex(A)
            @inbounds B[I] = log1(A[I])
        end
    else
        Threads.@threads for I in eachindex(A)
            @inbounds B[I] = log1(A[I])
        end
    end
    B
end

# These are a little faster than Julia's built-in functions?
exp1(x::Float64) = ccall(:exp, Cdouble, (Cdouble,), x)
exp1(x) = exp(x)
log1(x::Float64) = ccall(:log, Cdouble, (Cdouble,), x)
log1(x) = log(x)

# Versions which use cache
exp_(name::Symbol, A) = exp!(similar_(name, A), A)
log_(name::Symbol, A) = log!(similar_(name, A), A)

"""
    inv!(A) ≈ 1 ./ A
    inv!(A, b::Number) ≈ b ./ A

And `inv_(A)` which copies, and `inv0(A)` simple broadcasting.
Multi-threaded when `length(A) >= $TH_INV`.
Will be handled by `AppleAccelerate` if you load it.
"""
function inv! end

inv0(A::AbstractArray, b::Number=1) = similar(A) .= b ./ A # maps Adjoint -> Adjoint etc

@doc @doc(inv!)
inv_(A::AbstractArray, b::Number=1) = inv!(similar(A), A, b)
inv_(a::Number) = 1/a # for iscale_

inv!(A::AbstractArray, b::Number=1) = inv!(A, A, 1)
inv!(a::Number) = 1/a
function inv!(C::AbstractArray, A::AbstractArray, b::Number=1)
    @assert size(A)==size(C)
    if length(A) < TH_INV
        for I in eachindex(A)
            @inbounds C[I] = b / A[I]
        end
    else
        Threads.@threads for I in eachindex(A)
            @inbounds C[I] = b / A[I]
        end
    end
    C
end

inv_(name::Symbol, A::AbstractArray, b::Number=1) = inv!(similar_(name, A), A, b)
inv_(name::Symbol, a::Number=1) = 1/a # for iscale_

#========== scale! iscale! ==========#

export scale0, scale_, scale!, scale!!
export iscale0, iscale_, iscale!, iscale!!

using LinearAlgebra: Adjoint, Transpose

const ARVector = Union{Adjoint{<:Any, <:AbstractVector}, Transpose{<:Any, <:AbstractVector}}
const RVector = Union{Adjoint{<:Any, <:Vector}, Transpose{<:Any, <:Vector}}

"""
    scale!(A, b::Number) ≈ A .* b
    scale!(M, v::Vector) ≈ A .* v       # M::Matrix
    scale!(M, r::Adjoint) ≈ A .* r      # r::RowVector / Transpose etc.
    scale!(A, B) ≈ A .* B               # A,B same ndims

In-place scaling by a constant or (in the case of a matrix) by a row- or column-vector.
For each of these, there is also also `scale_(A, ...)` non-mutating but perhaps accellerated,
and `scale0(A, ...)` simple broadcasting.
"""
function scale! end

using LinearAlgebra

scale0(A::AbstractArray, b) = similar(A) .= A .* b

scale_(A::Array, b::Number) = rmul!(copy(A), b)
scale!(A::Array, b::Number) = rmul!(A, b)
scale!!(A::Array, b) = scale!(A,b) # differs in gradient

scale_(A::RVector, b::Number) = rmul!(copy(A), b) # scale_(::Abstract...) causes flux ambiguities
scale!(A::RVector, b::Number) = rmul!(A, b)
scale!!(A::RVector, b) = scale!(A,b)

@doc @doc(scale!)
scale_(A::Matrix, v::Vector) = lmul!(Diagonal(v), copy(A))
scale!(A::Matrix, v::Vector) = lmul!(Diagonal(v), A)

scale_(A::Matrix, r::RVector) = rmul!(copy(A), Diagonal(transpose(r)))
scale!(A::Matrix, r::RVector) = rmul!(A, Diagonal(transpose(r)))

# scale_(A::AbstractArray{T,N}, B::AbstractArray{T,N}) where {T,N} = similar(A) .= A .* B
# scale!(A::AbstractArray{T,N}, B::AbstractArray{T,N}) where {T,N} = A .= A .* B

function scale!(C::AbstractArray{T,N}, A::AbstractArray{TA,N}, B::AbstractArray{TB,N}) where {T,N, TA, TB}
    @assert size(A) == size(B) == size(C)
    for i in eachindex(A)
        @inbounds C[i] = A[i] * B[i]
    end
    C
end
scale_(A::AbstractArray{T,N}, B::AbstractArray{T,N}) where {T,N} = scale!(similar(A), A, B)
scale!(A::AbstractArray{T,N}, B::AbstractArray{T,N}) where {T,N} = scale!(A, A, B)

scale_(name::Symbol, A::Array, b::Number) = rmul!(copy_(name, A), b)
scale_(name::Symbol, A::Matrix, v::Vector) = lmul!(Diagonal(v), copy_(name, A))
scale_(name::Symbol, A::Matrix, r::RVector) = rmul!(copy_(name, A), Diagonal(transpose(r)))
scale_(name::Symbol, A::AbstractArray{T,N}, B::AbstractArray{T,N}) where {T,N} =
    scale!(similar_(name, A), A, B)

"""
    iscale!(A, b::Number) ≈ A ./ b
    iscale!(A, v::Vector) ≈ A ./ v      # A::Matrix
    iscale!(A, r::Adjoint) ≈ A ./ r     # r::RowVector / Transpose etc.
    iscale!(A, B) ≈ A ./ B

For each of these, there is also `iscale_(A, ...)` non-mutating but perhaps accellerated,
and `iscale0(A, ...)` simple broadcasting.
Finally there is `iscale!!(A, x)` which mutate both arguments, wihch may be a terrible idea.
"""
function iscale! end

iscale0(A::AbstractArray, b) = similar(A) .= A ./ b

@doc @doc(iscale!)
iscale_(A::AbstractArray, b) = scale_(A, inv_(b))
iscale!(A::AbstractArray, b) = scale!(A, inv_(b))
iscale!!(A::AbstractArray, b) = scale!(A, inv!(b))

function iscale!(C::Array{T,N}, A::Array{TA,N}, B::Array{TB,N}) where {T,N, TA, TB}
    @assert size(A) == size(B) == size(C)
    for i in eachindex(A)
        @inbounds C[i] = A[i] / B[i]
    end
    C
end
iscale_(A::Array{T,N}, B::Array{T,N}) where {T,N} = iscale!(similar(A), A, B)
iscale!(A::Array{T,N}, B::Array{T,N}) where {T,N} = iscale!(A, A, B)

# On square matrices this is a tie, but on (n,N) ./ N' it is faster
# The equivalent for iscale(Matrix, Vector) however is slower, wrong order
function iscale!(C::Matrix, A::Matrix, r::RVector)
    @assert size(A)==size(C)
    axes(A,2) == axes(r,2) || throw(DimensionMismatch("size disagreement in iscale?(Matrix,RowVector)"))
    @inbounds for j in axes(A,2)
        invr = inv(r[j])
        for i in axes(A,1)
            C[i,j] = A[i,j] * invr
        end
    end
    C
end
iscale_(A::Matrix, r::RVector) = iscale!(similar(A), A, r)
iscale!(A::Matrix, r::RVector) = iscale!(A, A, r)

iscale_(name::Symbol, A::AbstractArray, b) = scale_(name, A, inv_(name, b))
iscale_(name::Symbol, A::Matrix, r::RVector) = iscale!(similar_(name, A), A, r)
iscale_(name::Symbol, A::Array{T,N}, B::Array{T,N}) where {T,N} =
    iscale!(similar_(name, A), A, B)

#========== Accelerators ==========#

const CFloat = Union{Float64, Float32}
const CFloatArray{N} = Array{<:CFloat, N}
const CFloatMatrix = Matrix{<:CFloat}

IVERBOSE = false

VEC = ""
function load_note(str)
    global VEC
    if VEC == ""
        @info "ArrayAllez loaded code for $str"
        VEC = str
    else
        @warn "ArrayAllez loaded code for $str, perhaps overwriting $VEC"
        VEC *= " then $str"
    end
end

using Requires

@init @require Yeppp = "6310b701-1812-5374-a82f-9f6f2d54a40a" begin
    using .Yeppp

    exp!(B::CFloatArray, A::CFloatArray) = Yeppp.exp!(B, A)

    log!(B::CFloatArray, A::CFloatArray) = Yeppp.log!(B, A) # log_(A) calls log!(B,A)

    scale_(A::Array{T,N}, B::Array{T,N}) where {T<:CFloat,N} = Yeppp.multiply(A,B)
    scale!(A::Array{T,N}, B::Array{T,N}) where {T<:CFloat,N} = Yeppp.multiply!(A,A,B)

    IVERBOSE && load_note("Yeppp")
end

@init @require AppleAccelerate = "13e28ba4-7ad8-5781-acae-3021b1ed3924" begin
    using .AppleAccelerate

    exp!(B::CFloatArray, A::CFloatArray) = AppleAccelerate.exp!(B, A)

    log!(B::CFloatArray, A::CFloatArray) = AppleAccelerate.log!(B, A)

    inv!(A::CFloatArray) = AppleAccelerate.rec!(A, A)

    scale_(A::Vector{T}, B::Vector{T}) where {T<:CFloat} = AppleAccelerate.vmul(A,B)
    scale!(A::Vector{T}, B::Vector{T}) where {T<:CFloat} = AppleAccelerate.vmul!(A,A,B)

    scale_(A::Array{T,N}, B::Array{T,N}) where {T<:CFloat,N} =
        reshape(AppleAccelerate.vmul(vec(A),vec(B)), size(A))  # vmul is literally only vectors
    scale!(A::Array{T,N}, B::Array{T,N}) where {T<:CFloat,N} =
        begin AppleAccelerate.vmul!(vec(A),vec(A),vec(B)); A end

    iscale_(A::Vector{T}, B::Vector{T}) where {T<:CFloat} = AppleAccelerate.vdiv(A,B)
    iscale!(A::Vector{T}, B::Vector{T}) where {T<:CFloat} = AppleAccelerate.vdiv!(A,A,B)

    iscale_(A::Array{T,N}, B::Array{T,N}) where {T<:CFloat,N} =
        reshape(AppleAccelerate.vdiv(vec(A),vec(B)), size(A))
    iscale!(A::Array{T,N}, B::Array{T,N}) where {T<:CFloat,N} =
        begin AppleAccelerate.vdiv!(vec(A),vec(A),vec(B)); A end

    IVERBOSE && load_note("AppleAccelerate")
end

@init @require IntelVectorMath = "c8ce9da6-5d36-5c03-b118-5a70151be7bc" begin
    using .IntelVectorMath

    exp!(B::CFloatArray, A::CFloatArray) = IVM.exp!(B, A)

    log!(B::CFloatArray, A::CFloatArray) = IVM.log!(B, A) # log_(A) calls log!(B,A)

    iscale_(A::Array{T,N}, B::Array{T,N}) where {T<:CFloat,N} = IVM.divide(A,B)
    iscale!(A::Array{T,N}, B::Array{T,N}) where {T<:CFloat,N} = IVM.divide!(A,A,B)

    IVERBOSE && load_note("IntelVectorMath")
end

#========== The End ==========#


text:
module TalusPlotScript

using SynapseClient
using DISSEQT
using DISSEQT.Plots
using DISSEQT.AnnotatedArrays
using DISSEQT.SynapseTools
using Statistics
using DataFrames
using CSV
using JLD
using Gadfly
using Cairo
using Fontconfig

function cachemetadata(syn,alignmentFolder,metadataID)
    metadataCacheFilename = "metadatacache.csv"
    metadata = DataFrame()
    if isfile(metadataCacheFilename)
        metadata = CSV.read(metadataCacheFilename; missingstrings=["","NA"])
    else
        # Get metadata (metadataID can be the synapse ID of a .csv file or of a folder containing .csv files).
        metadata = getsamplemetadata(syn, metadataID, :Rejected=>"No") # Always remove rejected samples

        # extend metadata
        appendswarmids!(syn,metadata,alignmentFolder)
        appendalignedids!(syn,metadata,alignmentFolder)

        # download swarms and consensus sequences
        downloadswarms!(syn,metadata,cacheDir="MutantSwarms")
        downloadconsensuses!(syn,metadata,cacheDir="MutantSwarms")

        CSV.write(metadataCacheFilename, metadata; missingstring="NA")
    end
    metadata
end

function main()
    syn = SynapseClient.login()

    # projectFolder should point to MyProject folder (containing subfolders Metadata and Analysis)
    projectFolder   = "syn11639899" # VignuzziLabPublic/Projects/FitnessLandscapes
    analysisFolder  = getchildbyname(syn, projectFolder, "Analysis")
    alignmentFolder = getchildbyname(syn, analysisFolder, "Alignment")
    referenceFolder = getchildbyname(syn, alignmentFolder, "ReferenceGenomes")

    isdir("MutantSwarms") || mkdir("MutantSwarms")
    isdir("plots") || mkdir("plots")

    # set to true to upload files
    doUpload = true
    uploadName = "fitness landscapes invitro"

    metadataID = getchildbyname(syn, projectFolder, "Metadata", "fitness landscapes invitro.csv")
    metadata = cachemetadata(syn, alignmentFolder, metadataID)

    # Get Reference Genomes
    referenceNames = ["WT", "Less", "More", "Stop"]
    referenceIDs, referenceGenomes = getreferencegenomes(syn, referenceFolder, referenceNames)

    # get CVB3 feature list
    featureFileID = getchildbyname(syn,referenceFolder,"CVB3_features.csv")
    features = CSV.read(localpath(syn, featureFileID); missingstrings=["","NA"])

    positions = featurepositions(features, "ORF") # coding positions

    # load swarms
    swarms  = loadswarm(metadata[:SampleID],metadata[:SwarmPath],metadata[:ConsensusPath])

    # filter by positions
    swarms = swarms[:,:,:,inset(:position,positions),:]

    coverage = dropdims(swarms[:coverage]; dims=(1,2,3))
    meanCoverage = mean(coverage; dims=1)[:]

    # filter by mean read coverage at coding sites
    sampleMask = meanCoverage .>= 1000
    println("Removed ", count(.~sampleMask), " samples due to low coverage.")
    swarms  = swarms[:,:,:,:,sampleMask]
    metadata = metadata[sampleMask,:]

    # remove mixed samples
    mixes = estimatemix(swarms, referenceGenomes)
    mixMask = maximum(mixes;dims=2)[:] .>= 0.98
    println("Removed ", count(.~mixMask), " mixed samples.")
    swarms  = swarms[:,:,:,:,mixMask]
    metadata = metadata[mixMask,:]

    println("Total number of samples used: ", size(swarms,5))

    # load limit of detection
    limitOfDetectionID = getchildbyname(syn,analysisFolder,"LimitOfDetection","fitness landscapes invitro","limitofdetection.jld")
    limitOfDetection = convert(AnnotatedArray,load(localpath(syn,limitOfDetectionID,downloadLocation="MutantSwarms")))
    limitOfDetection = limitOfDetection[:,:,:,inset(:position,positions)] # filter by positions
    @assert isequal(swarms[:position][:], limitOfDetection[:position][:]) # the positions must be the same

    # get rid of annotations to simplify math stuff
    # arrays are 4x4x4xPxN - nuc3 x nuc2 x nuc1 x positions x samples
    s = convert(Array, swarms)
    α = convert(Array, limitOfDetection)

    # Put a lower limit on the limit of detection
    α = max.(α,1e-3)

    # apply transformation
    X = log2.(s.+α)

    # center variables
    X = X .- mean(X;dims=5) # remove mean over samples

    # reshape to VxN - variables x samples
    X = reshape(X, 64*size(X,4), size(X,5))

    # Non-intrusive prefiltering of variables that will never contribute
    σThreshold = 1e-2 # same as the lower limit we use for projection score
    σ = dropdims(std(X;dims=2);dims=2)
    σ = σ/maximum(σ) # normalize
    X = X[σ.>=σThreshold,:]

    uploads = UploadList()
    dependencies = unique(vcat(featureFileID, referenceIDs, limitOfDetectionID, metadata[:SwarmID], metadata[:MetadataID]))

    nbrDims = 40

    pl = talusplot(X, 1:nbrDims, Guide.title("Talus plot"))
    markforupload!(uploads, "plots", saveplot(pl, [:png,:svg,:pdf], "plots/talus", width=29.7cm/2, height=21cm/2), dependencies)

    doUpload && uploadfiles(syn, createchildfolder(syn, analysisFolder, "Dimensionality", uploadName), uploads, @__FILE__, "Talus plot")

    nothing
end

end
TalusPlotScript.main();


text:
mutable struct OrbitSystem{T, WT<:AbstractArray, ST<:AbstractArray,
                           SiT<:AbstractArray, GT<:AbstractArray,
                           TT<:AbstractArray}
    # sum(((W*f .- d)./err).^2) + (f .- mu)'*inv(alpha*S)*(f .- mu)
    W::WT
    d::Vector{T}
    err::Vector{T}
    S::ST
    S_inv::SiT
    G::GT
    T::TT
    mu::Vector{T}
    alpha::T
end

function Base.show(io::IO, OS::OrbitSystem)
    println(io, "OrbitSystem")
    println(io, " # of measurements: $(length(OS.d))")
    print(io,   " # of orbits: $(size(OS.W,2))")
end

function OrbitSystem(W, d, err, S; S_inv = inv(S), G = cholesky(S_inv).U,
                     T=Diagonal(ones(size(W,2))), mu = zeros(size(W,2)), alpha=1.0)

    OrbitSystem(W,d,err,S,S_inv,G,T,mu,alpha)
end

function lcurve_point(OS::OrbitSystem, f::Vector)
    W = OS.W
    d = OS.d
    err = OS.err
    mu = OS.mu
    S_inv = OS.S_inv

    chisq = sum(((W*f .- d)./err).^2)
    R = (f .- mu)'*S_inv*(f .- mu)
    return chisq, R
end

function lcurve_point(OS::OrbitSystem; kwargs...)
    f = solve(OS; kwargs...)
    return lcurve_point(OS,f)
end

function lcurve(OS; n = 50, log_bounds = (-6,6), kwargs...)
    alpha = 10.0.^range(log_bounds...,length=n)
    chisq = zeros(n)
    R = zeros(n)
    for i=1:n
        OS.alpha = alpha[i]
        chisq[i], R[i] = lcurve_point(OS; kwargs...)
    end
    return alpha, chisq, R
end

function marginal_loglike(OS::OrbitSystem; norm=1e18/size(OS.W,2),nonneg=false,max_iter=30*size(OS.W,1))

    d = OS.d
    err = OS.err
    Σ_d = Diagonal(err.^2)
    Σ_d_inv = inv(Σ_d)

    K = norm*OS.W

    mu_X = OS.mu/norm

    # Scale covariance matrices
    Σ_X = OS.alpha*OS.S
    Σ_X_inv = inv(OS.alpha)*OS.S_inv

    Σ_inv = Σ_X_inv .+ K'*Σ_d_inv*K

    if nonneg
        Γ = sqrt(inv(OS.alpha))*OS.G
        X = vec(nonneg_lsq(vcat(K./err, Γ), vcat(d./err, Γ*mu_X),
                           alg=:fnnls,use_parallel=false,max_iter=max_iter))
    else
        X = Σ_inv \ (K'*(Σ_d_inv*d) + Σ_X_inv*mu_X)
    end

    l = 0.0
    try
        dhat = K*X
        l = -logabsdet(Array(Σ_X))[1] - logabsdet(Array(Σ_inv))[1] -
            ((d .- dhat)'*Σ_d_inv*(d .- dhat))[1] - ((X .- mu_X)'*Σ_X_inv*(X .- mu_X))[1]
    catch er
        if isa(er,InterruptException)
            rethrow(er)
        end
        l = -Inf
    end

    return l
end

function optimize_alpha!(OS; log_bounds = (-6,6), kwargs...)
    f = x -> begin
        OS.alpha = 10.0^x
        return -marginal_loglike(OS; kwargs...)
    end

    op = optimize(f, log_bounds[1], log_bounds[2], Brent(),rel_tol=1e-3)

    OS.alpha = 10.0^Optim.minimizer(op)

    return Optim.minimum(op)
end

function estimate_rtol(S)
    s = filter(x -> x > 0, eigvals(S))
    sn = log10.(reverse(s)/s[end])
    n = length(sn)
    x = collect(range(0,1,length=n))
    f = PolyharmonicSpline(3,x,sn,s=0.1)
    df2(t) = ForwardDiff.derivative(tt->ForwardDiff.derivative(f,tt),t)
    roots = find_zeros(df2,0.2,0.8)
    return 10.0^(f(roots[1]))
end

function optimize_parameters(make_orbit_system::Function, lbounds, ubounds;
                             niter=10, nseed = 15, objective::Function = optimize_alpha!,
                             warmstart=false, checkpoint=true,
                             file="optim_progress.jld2",verbose=true, kwargs...)

    s = Sobol.SobolSeq(lbounds,ubounds)
    points = Vector{Float64}[]
    values = Float64[]

    if warmstart && isfile(file) && (filesize(file) != 0)
        s, points, values = load(file, "s", "points", "values")
    end

    if checkpoint
        touch(file)
    end

    while s.s.n < nseed
        p = Sobol.next!(s)
        OS = make_orbit_system(p)
        v = objective(OS; kwargs...)
        push!(points, p)
        push!(values, v)
        if checkpoint
            @save file s points values
        end
        if verbose
            println("Seed Number: $(s.s.n)")
            println("Seed Point:  ", repr(p))
            println("Seed Value:  ", v)
            println("Best Seed: ", repr(points[argmin(values)]))
            println("")
        end
    end

    for i=1:niter
        spl = PolyharmonicSpline(2, Array(hcat(points...)'), values)
        xstart = s.lb .+ s.ub.*rand(length(s.ub))
        op = optimize(x -> spl(x...), s.lb, s.ub, xstart, Fminbox(NelderMead()))
        p = Optim.minimizer(op)
        pv = Optim.minimum(op)
        if p in points
            p = Sobol.next!(s)
            pv = spl(p...)
        end
        if verbose
            println("Iteration: $i")
            println("Guess Point:  ",repr(p))
            println("Guess Value:  ", pv)
        end
        OS = make_orbit_system(p)
        v = objective(OS; kwargs...)
        if isfinite(v)
            push!(points, p)
            push!(values, v)
        end
        if verbose
            println("Actual Value: ", v)
            println("Best Point: ", repr(points[argmin(values)]))
            println("")
        end
        if checkpoint
            @save file s points values
        end
    end

    return points[argmin(values)]
end

function inv_chol(Σ; rtol=0.0)
    if rtol == 0.0
        Σ_inv = inv(Σ)
        Γ = cholesky(Hermitian(Σ_inv)).U
    else
        S = svd(Σ)
        smax = maximum(S.S)
        svals_inv = Diagonal([s >= smax*rtol ? 1/s : 0.0 for s in S.S])
        Σ_inv = (svals_inv*S.Vt)' * S.U'
        Γ = sqrt.(svals_inv)*S.Vt
    end
    return Σ_inv, Γ
end

function inv_chol(S::LinearAlgebra.SVD; rtol=0.0)
    if rtol == 0.0
        Σ_inv = inv(S)
        Γ = cholesky(Hermitian(Σ_inv)).U
    else
        smax = maximum(S.S)
        svals_inv = Diagonal([s >= smax*rtol ? 1/s : 0.0 for s in S.S])
        Σ_inv = (svals_inv*S.Vt)' * S.U'
        Γ = sqrt.(svals_inv)*S.Vt
    end
    return Σ_inv, Γ
end

function solve(OS::OrbitSystem; norm=1e18/size(OS.W,2), nonneg=true,max_iter=30*size(OS.W,1))

    d = OS.d
    err = OS.err
    Σ_d = Diagonal(err.^2)
    Σ_d_inv = inv(Σ_d)

    K = norm*OS.W

    mu_X = OS.mu/norm

    # Scale covariance matrices
    Σ_X = OS.alpha*OS.S
    Σ_X_inv = inv(OS.alpha)*OS.S_inv
    Γ = sqrt(inv(OS.alpha))*OS.G

    if nonneg
        try
            X = vec(nonneg_lsq(vcat(K./err, Γ), vcat(d./err, Γ*mu_X),
                               alg=:fnnls,use_parallel=false,max_iter=max_iter))
        catch er
            if isa(er,InterruptException)
                rethrow(er)
            end
            @warn "Non-negative Least Squares failed. Using MAP estimate"
            println(err)
            Σ_inv = Σ_X_inv .+ K'*Σ_d_inv*K
            X = Σ_inv \ (K'*(Σ_d_inv*d) + Σ_X_inv*mu_X)
        end
    else
        Σ_inv = Σ_X_inv .+ K'*Σ_d_inv*K
        X = Σ_inv \ (K'*(Σ_d_inv*d) + Σ_X_inv*mu_X)
    end

    return norm*X
end


text:
# CUDA pointer types

export CuPtr, CU_NULL, PtrOrCuPtr

#
# CUDA pointer
#

"""
    CuPtr{T}

A memory address that refers to data of type `T` that is accessible from the GPU. A `CuPtr`
is ABI compatible with regular `Ptr` objects, e.g. it can be used to `ccall` a function that
expects a `Ptr` to GPU memory, but it prevents erroneous conversions between the two.
"""
CuPtr

if sizeof(Ptr{Cvoid}) == 8
    primitive type CuPtr{T} 64 end
else
    primitive type CuPtr{T} 32 end
end

# constructor
CuPtr{T}(x::Union{Int,UInt,CuPtr}) where {T} = Base.bitcast(CuPtr{T}, x)

const CU_NULL = CuPtr{Cvoid}(0)

## getters

Base.eltype(::Type{<:CuPtr{T}}) where {T} = T

## conversions

# to and from integers
## pointer to integer
Base.convert(::Type{T}, x::CuPtr) where {T<:Integer} = T(UInt(x))
## integer to pointer
Base.convert(::Type{CuPtr{T}}, x::Union{Int,UInt}) where {T} = CuPtr{T}(x)
Int(x::CuPtr)  = Base.bitcast(Int, x)
UInt(x::CuPtr) = Base.bitcast(UInt, x)

# between regular and CUDA pointers
Base.convert(::Type{<:Ptr}, p::CuPtr) =
    throw(ArgumentError("cannot convert a GPU pointer to a CPU pointer"))

# between CUDA pointers
Base.convert(::Type{CuPtr{T}}, p::CuPtr) where {T} = Base.bitcast(CuPtr{T}, p)

# defer conversions to unsafe_convert
Base.cconvert(::Type{<:CuPtr}, x) = x

# fallback for unsafe_convert
Base.unsafe_convert(::Type{P}, x::CuPtr) where {P<:CuPtr} = convert(P, x)

## limited pointer arithmetic & comparison

isequal(x::CuPtr, y::CuPtr) = (x === y)
isless(x::CuPtr{T}, y::CuPtr{T}) where {T} = x < y

Base.:(==)(x::CuPtr, y::CuPtr) = UInt(x) == UInt(y)
Base.:(<)(x::CuPtr,  y::CuPtr) = UInt(x) < UInt(y)
Base.:(-)(x::CuPtr,  y::CuPtr) = UInt(x) - UInt(y)

Base.:(+)(x::CuPtr, y::Integer) = oftype(x, Base.add_ptr(UInt(x), (y % UInt) % UInt))
Base.:(-)(x::CuPtr, y::Integer) = oftype(x, Base.sub_ptr(UInt(x), (y % UInt) % UInt))
Base.:(+)(x::Integer, y::CuPtr) = y + x

#
# GPU or CPU pointer
#

"""
    PtrOrCuPtr{T}

A special pointer type, ABI-compatible with both `Ptr` and `CuPtr`, for use in `ccall`
expressions to convert values to either a GPU or a CPU type (in that order). This is
required for CUDA APIs which accept pointers that either point to host or device memory.
"""
PtrOrCuPtr

if sizeof(Ptr{Cvoid}) == 8
    primitive type PtrOrCuPtr{T} 64 end
else
    primitive type PtrOrCuPtr{T} 32 end
end

function Base.cconvert(::Type{PtrOrCuPtr{T}}, val) where {T}
    # `cconvert` is always implemented for both `Ptr` and `CuPtr`, so pick the first result
    # that has done an actual conversion

    gpu_val = Base.cconvert(CuPtr{T}, val)
    if gpu_val !== val
        return gpu_val
    end

    cpu_val = Base.cconvert(Ptr{T}, val)
    if cpu_val !== val
        return cpu_val
    end

    return val
end

function Base.unsafe_convert(::Type{PtrOrCuPtr{T}}, val) where {T}
    # FIXME: this is expensive; optimize using isapplicable?
    ptr = try
        Base.unsafe_convert(Ptr{T}, val)
    catch
        try
            Base.unsafe_convert(CuPtr{T}, val)
        catch
            throw(ArgumentError("cannot convert to either a CPU or GPU pointer"))
        end
    end
    return Base.bitcast(PtrOrCuPtr{T}, ptr)
end