text:
module Moves

import Magic
import Util

global bishop_attacks = Magic.initSlidersAttack(true)
global rook_attacks = Magic.initSlidersAttack(false)

function rookMoves(occ, sq)
    #could use whole occupancy board,but would map to too many possibilties, following line gives us a reduced occupancy that is much more manageable
    #use same magic number while instantiating rook_moves array as for rook_magic_numbers
    if(sq == 0) return 0 end
    occ &= UInt128(Magic.rMask(sq))
    occ *= Magic.rook_magic_numbers[sq]
    occ >>= (64 - Magic.n_r_bits[sq])
    return rook_attacks[sq][occ]
    end

function bishopMoves(occ,sq)
    if(sq == 0) return 0 end
    sq = Int(sq) #i dont know why this is necessary, when sq = 41, some error comes up?
    occ &= UInt128(Magic.bMask(sq))
    occ *= Magic.bishop_magic_numbers[sq]
    occ >>= (64 - Magic.n_b_bits[sq])
    return bishop_attacks[sq][occ]
end

function generatePossibleMoves(moves,f1,type)
    possible_moves = []
    for index in 1:64
        el = Util.getBit(moves,index)
        if(el == 1)
            start = f1(index)
            if(type ==='P')
                push!(possible_moves, Util.BBToUCI(start,index) * 'Q')
            else
                push!(possible_moves, Util.BBToUCI(start,index))
            end
        end
    end
    return possible_moves
end

function getPawnMoves(board,colour::Char,history)
    possible_moves = []
    opp_colour = (colour === 'W') ? 'B' : 'W'
    opposing_pieces = board[opp_colour * 'P'] | board[opp_colour * 'N'] | board[opp_colour * 'B'] | board[opp_colour * 'R'] | board[opp_colour * 'Q'] # omit K to avoid illegal capture
    empty = ~(board[colour * 'P'] | board[colour * 'N'] | board[colour * 'B'] | board[colour * 'R'] | board[colour * 'Q'] | board[colour * 'K'] | opposing_pieces | board[opp_colour * 'K'])
    shift_function = (colour === 'W') ? (x,y) -> x << y  : (x,y) -> x >> y
    last_rank = (colour === 'W') ? Util.RANK_MASKS[8] : Util.RANK_MASKS[1]
    left_side = (colour === 'W') ? Util.FILE_MASKS[1] : Util.FILE_MASKS[8]
    right_side = (colour === 'W') ? Util.FILE_MASKS[8] : Util.FILE_MASKS[1]
    two_forward = (colour === 'W') ? Util.RANK_MASKS[4] : Util.RANK_MASKS[5]
    operator = (colour === 'W') ? (a,b) -> a - b : (a,b) -> a + b
    en_passant = (colour === 'W') ? (a,b) -> a + b : (a,b) -> a - b
    pieces_to_move = board[colour * 'P']
    #check for en passant
    if(length(history) > 0)
        last_move = history[length(history)]
        start_index, end_index = Util.UCIToBB(last_move)
    else
        last_move = missing
        start_index, end_index = 0,0
    end
        #if last move was a pawn moving two sqs

    if (~ismissing(last_move) && (Util.getBit(board[opp_colour * 'P'],end_index) == 1) && abs(end_index - start_index) == 16)
        if (Util.getBit(pieces_to_move & ~Util.FILE_MASKS[1], end_index + 1) == 1)  #to the right if white, to the left if black
            push!(possible_moves,Util.BBToUCI(end_index + 1, en_passant(end_index, 8)))
        elseif (Util.getBit(pieces_to_move & ~Util.FILE_MASKS[8], end_index - 1) == 1)
            push!(possible_moves,Util.BBToUCI(end_index- 1, en_passant(end_index, 8)))
        end
    end

    #capture right
    moves_right = shift_function(pieces_to_move, 7) & opposing_pieces & (~last_rank) & (~right_side)
    f1 = num -> operator(num,7)
    possible_moves = append!(possible_moves, generatePossibleMoves(moves_right,f1,'N'))
    #capture left
    cap_left = shift_function(pieces_to_move, 9) & opposing_pieces & ~last_rank & ~left_side
    f1 = num -> operator(num,9)
    possible_moves = append!(possible_moves,generatePossibleMoves(cap_left, f1, 'N'))
    #1 forward
    one_forward = shift_function(pieces_to_move, 8) & empty & ~(last_rank)
    f1 = num -> operator(num, 8)
    possible_moves = append!(possible_moves,generatePossibleMoves(one_forward,f1,'N'))
    #2 forward
    two_forward = shift_function(pieces_to_move, 16) & (empty) & shift_function(empty,8) & two_forward
    f1 = num -> operator(num,16)
    possible_moves = append!(possible_moves,generatePossibleMoves(two_forward, f1,'N'))
    #pawn promotion by capture right
    promo_right = shift_function(pieces_to_move, 7) & (opposing_pieces) & last_rank & ~right_side
    f1 = num -> operator(num, 7)
    possible_moves = append!(possible_moves,generatePossibleMoves(promo_right, f1,  'P'))
    #pawn promotion by 1 forward
    promo_forward = shift_function(pieces_to_move, 8) & empty & last_rank
    f1 = num -> operator(num, 8)
    possible_moves = append!(possible_moves,generatePossibleMoves(promo_forward, f1,  'P'))
    #pawn promotion by capture left
    promo_cap_left = shift_function(pieces_to_move, 9) & opposing_pieces & last_rank & ~left_side
    f1 = num -> operator(num, 9)
    possible_moves = append!(possible_moves,generatePossibleMoves(promo_cap_left, f1,  'P'))
    return possible_moves
end

function pawnAttacks(pawn_board,colour)
    if(pawn_board == 0) return 0 end
    attacks = 0
    if (colour === 'W')
        attacks |= (pawn_board << 7) & ~Util.FILE_MASKS[8]
        attacks |= (pawn_board << 9) & ~Util.FILE_MASKS[1]
    else
        attacks |= (pawn_board >> 7) & ~Util.FILE_MASKS[1]
        attacks |= (pawn_board >> 9) & ~Util.FILE_MASKS[8]
    end
    return [[attacks, 0, "P"]]
end

function pawnAttacksFast(pawn_board,colour)
    if(pawn_board == 0) return 0 end
    attacks = 0
    if (colour === 'W')
        attacks |= (pawn_board << 7) & ~Util.FILE_MASKS[8]
        attacks |= (pawn_board << 9) & ~Util.FILE_MASKS[1]
    else
        attacks |= (pawn_board >> 7) & ~Util.FILE_MASKS[1]
        attacks |= (pawn_board >> 9) & ~Util.FILE_MASKS[8]
    end
    return attacks
end

function knightMoves(sq)
    if(sq == 0) return 0 end

    attacks = 0
    board = Util.setBit(0,sq,1)
    if ((board >> 17) & ~Util.FILE_MASKS[8] > 0) attacks |= (board >> 17) end
    if ((board >> 15) & ~Util.FILE_MASKS[1] > 0) attacks |= (board >> 15) end
    if ((board >> 10) & ~Util.FILE_MASKS[8] & ~Util.FILE_MASKS[7] > 0) attacks |= (board >> 10) end
    if ((board >> 6) & ~Util.FILE_MASKS[1] & ~Util.FILE_MASKS[2] > 0 ) attacks |= (board >> 6) end
    if ((board << 17) & ~Util.FILE_MASKS[1] > 0 ) attacks |= (board << 17) end
    if ((board << 15) & ~Util.FILE_MASKS[8] > 0 ) attacks |= (board << 15) end
    if ((board << 10) & ~Util.FILE_MASKS[1] & ~Util.FILE_MASKS[2] > 0) attacks |= (board << 10) end
    if ((board << 6) & ~Util.FILE_MASKS[8] & ~Util.FILE_MASKS[7] > 0) attacks |= (board << 6) end
    return attacks
end

function queenMoves(occ,sq)
    if(sq == 0) return 0 end

    return bishopMoves(occ,sq) | rookMoves(occ,sq)
end

function kingMoves(sq)
    if(sq == 0) return 0 end
    attacks = 0
    board = 2^(sq-1)
    if ((board >> 1) & ~Util.FILE_MASKS[8] > 0 ) attacks |= (board >> 1) end
    if ((board >> 7) & ~Util.FILE_MASKS[1] > 0 ) attacks |= (board >> 7) end
    attacks |= (board >> 8)
    if ((board >> 9) & ~Util.FILE_MASKS[8] > 0 ) attacks |= (board >> 9) end
    if ((board << 1) & ~Util.FILE_MASKS[1] > 0 ) attacks |= (board << 1) end
    if ((board << 7) & ~Util.FILE_MASKS[8] > 0) attacks |= (board << 7) end
    attacks |= (board << 8)
    if ((board << 9) & ~Util.FILE_MASKS[1] > 0 ) attacks |= (board << 9) end
    return attacks
end

function kingMovesActual(sq,attack_board,pieces)
    return kingMoves(sq) & ~(attack_board) & ~pieces
end

# function getPieceMoves(board,occ,pieces,f,piece_name)
#     ar = []
#     if(board == 0) #none of this type of piece exist
#         return []
#     end
#     piece_1 = mostSignificantBit(board)
#     piece_2 = leastSignificantBit(board) + 1

#     if((piece_1 == piece_2) )
#         # moves = f(occ,piece_1) & ~pieces
#         # while(moves > 0)
#         #     index = leastSignificantBit(moves)
#         #     moves = popBit(moves, index)
#         #     push!(ar,BBToUCI(piece_1,index + 1))
#         # end
#          return [[f(occ,piece_1) & ~pieces,piece_1,piece_name]]
#     else
#         # moves = f(occ,piece_1) & ~pieces
#         # while(moves > 0)
#         #     index = leastSignificantBit(moves)
#         #     moves = popBit(moves, index )
#         #     push!(ar,BBToUCI(piece_1,index +1))
#         # end
#         # moves = f(occ,piece_2) & ~pieces
#         # while(moves > 0)
#         #     index = leastSignificantBit(moves)
#         #     moves = popBit(moves, index )
#         #     push!(ar,BBToUCI(piece_2,index + 1))
#         # end
#          return [[f(occ,piece_1) & ~pieces, piece_1, piece_name],[f(occ,piece_2) & ~pieces, piece_2, piece_name]]
#     end
#     return ar
# end

function getPieceMovesFast(board,occ,pieces,f)
    ar = []
    if(board == 0) #none of this type of piece exist
        return []
    end
    piece_1 = Util.mostSignificantBit(board)
    piece_2 = Util.leastSignificantBit(board) + 1

    if((piece_1 == piece_2) )
        moves = f(occ,piece_1) & ~pieces
        while(moves > 0)
            index = Util.leastSignificantBit(moves)
            moves = Util.popBit(moves, index)
            push!(ar,Util.BBToUCI(piece_1,index + 1))
        end
    else
        moves = f(occ,piece_1) & ~pieces
        while(moves > 0)
            index = Util.leastSignificantBit(moves)
            moves = Util.popBit(moves, index )
            push!(ar,Util.BBToUCI(piece_1,index +1))
        end
        moves = f(occ,piece_2) & ~pieces
        while(moves > 0)
            index = Util.leastSignificantBit(moves)
            moves = Util.popBit(moves, index )
            push!(ar,Util.BBToUCI(piece_2,index + 1))
        end
    end
    return ar
end

function getMoves(board,colour,history)
    occ = reduce((x,y) -> x | y, values(board))
    opp_colour = (colour === 'W') ? 'B' : 'W'
    pieces =  board[colour * 'P'] | board[colour * 'N'] | board[colour * 'B'] | board[colour * 'R'] | board[colour * 'Q'] | board[colour * 'K']
    moves = append!(
        getPieceMovesFast(board[colour * 'K'],occ,pieces,(x,y) -> kingMovesActual(y, getAttackBoard(board,opp_colour, true), pieces)),
        getPieceMovesFast(board[colour * 'Q'],occ,pieces,queenMoves),
        getPieceMovesFast(board[colour *'N'],occ,pieces,(x,y) -> knightMoves(y)),
        getPieceMovesFast(board[colour * 'B'],occ,pieces,bishopMoves),
        getPieceMovesFast(board[colour * 'R'],occ,pieces,rookMoves)
        )
    moves_fin = append!(moves, getPawnMoves(board,colour,history))
    moves_fin = findNoCheckMoves(moves_fin,board,colour)
    moves_with_castle = append!(castleMoves(board,colour,history,occ),moves_fin)
    return moves_with_castle
end

function findNoCheckMoves(moves,board::Dict,colour::Char)::Array{String}
    moves_no_check = []
  #  moves_no_check = Array{String}(undef,length(moves))
    for index in 1:length(moves)
        #just look for pinned pieces
        start_index,end_index = Util.UCIToBB(moves[index])
        key1 = ""
        key2 = ""
        temp1 = 0
        temp2 = 0
        for key in keys(board)
            if(board[key] & (1 << (start_index - 1)) > 0)
                temp1 = board[key]
                key1 = key
                board[key] = Util.setBit(board[key], start_index,0)
                board[key] = Util.setBit(board[key],end_index,1)
            elseif(board[key] & (1 << (end_index - 1)) > 0 && ~('K' in key))
                temp2 = board[key]
                key2 = key
                board[key] = Util.setBit(board[key], end_index,0)
            end
        end
        if(~isCheck(board,colour))
            push!(moves_no_check,moves[index])
        end
        #undo moves
        if(key1 !== "")
        board[key1] = temp1
        end
        if(key2 !== "")
        board[key2] = temp2
        end
    end
    return moves_no_check
end

function getAttackBoard(board,colour,for_king)
    occ = reduce((x,y) -> x | y,values(board))
    opp_colour = (colour === 'W') ? 'B' : 'W'

    occ = (for_king) ? occ &= ~(board[opp_colour * 'K']) : occ
    pieces =  board[colour * 'P'] | board[colour * 'N'] | board[colour * 'B'] | board[colour * 'R'] | board[colour * 'Q'] | board[colour * 'K']
    attacks = 0
    for key in keys(board)
        temp = board[key]
        while(temp > UInt64(0))
            piece_1 = Util.leastSignificantBit(temp) + 1
            piece_2 = Util.mostSignificantBit(temp)
            temp = Util.popBit(temp,piece_1 - UInt64(1))
            temp = Util.popBit(temp,piece_2 - UInt64(1))
            #this is separated out, as using * for string concat is slow....
            if(colour == 'W')
                if("WQ" ===  key)
                    attacks |= queenMoves(occ,piece_1) & ~pieces
                elseif("WK" === key)
                    attacks |= kingMoves(piece_1) & ~pieces
                elseif("WB" === key)
                    attacks |= bishopMoves(occ,piece_1) & ~pieces
                    attacks |= bishopMoves(occ,piece_2) & ~pieces
                elseif("WN" === key)
                    attacks |= knightMoves(piece_1) & ~pieces
                    attacks |= knightMoves(piece_2) & ~pieces
                elseif("WR" === key)
                    attacks |= rookMoves(occ,piece_1) & ~pieces
                    attacks |= rookMoves(occ,piece_2) & ~pieces
                end
            elseif (colour == 'B')
                if("BQ" === key)
                    attacks |= queenMoves(occ,piece_1) & ~pieces
                elseif("BK" === key)
                    attacks |= kingMoves(piece_1) & ~pieces
                elseif("BB" === key)
                    attacks |= bishopMoves(occ,piece_1) & ~pieces
                    attacks |= bishopMoves(occ,piece_2) & ~pieces
                elseif( "BN" === key)
                    attacks |= knightMoves(piece_1) & ~pieces
                    attacks |= knightMoves(piece_2) & ~pieces
                elseif("BR" === key)
                    attacks |= rookMoves(occ,piece_1) & ~pieces
                    attacks |= rookMoves(occ,piece_2) & ~pieces
                end
            end
        end
    end
    attacks |= pawnAttacksFast(board[colour * 'P'],colour)
    return attacks
end

function castleMoves(board,colour,history,occ)
    moves = []
    rook_pos = (colour === 'W') ? 8 : 64
    opp_colour = (colour === 'W') ? 'B' : 'W'
    to_move = (colour === 'W') ? Util.setBitRange(0,6,7) : Util.setBitRange(0,62,63)
    king_pos = (colour === 'W') ? 5 : 61
    piece_moved = false
    for move in history
        start_index, end_index = Util.UCIToBB(move)
        if(Util.getBit(board[colour * 'K'],end_index) == 1 || Util.getBit(board[colour * 'R'],end_index) == 1) #make sure it is correct piece
            piece_moved = true
        end
    end

    #kingside
    if(Util.getBit(occ,king_pos) > 0  && Util.getBit(occ,rook_pos) > 0 && ~piece_moved && (occ & to_move) == 0 && (getAttackBoard(board,opp_colour,false) & to_move) == 0)
         #if both king and rook exist, no pieces between, no challenged squares, no piece has moved
         (colour === 'W') ? push!(moves,"e1g1") : push!(moves,"e8g8")
    end

    #queenside
    rook_pos = (colour === 'W') ? 1 : 57
    to_move = (colour === 'W') ? Util.setBitRange(0,2,4) : Util.setBitRange(0,58,60)
    if(Util.getBit(occ,king_pos) > 0 && Util.getBit(occ,rook_pos) > 0 && ~piece_moved && (occ & to_move) == 0 && (getAttackBoard(board,opp_colour,false) & to_move) == 0)
        colour === 'W' ? push!(moves,"e1c1") : push!(moves,"e8c8")
    end
    return moves
end

########### CHECK LOGIC ###############

function isCheck(board,colour)
    opp_colour = (colour === 'W') ? 'B' : 'W'
    attack = getAttackBoard(board,opp_colour,true)
    return ((board[colour * 'K'] & attack) > 0) ? true : false
end

function isCheckMate(board,colour,history)
    check = isCheck(board,colour)
    legal_moves = getMoves(board,colour,history)
    while(check && (length(legal_moves) > 0))
        move = legal_moves[1]
        board_state = makeMoveUCI(move,board,colour)
        if(!isCheck(board_state,colour))
            check = false
        else
            popfirst!(legal_moves)
        end
    end
    if (length(legal_moves) == 0 && check)
        return 1
    elseif (length(legal_moves) == 0 && ~check) #stalemate
        return 2
    else
        return 0
    end
end

function makeMoveUCI(move,board,colour)
    output = copy(board)
    start_index, end_index = Util.UCIToBB(move)
    #see if cap
    is_capture = false
    for key in keys(board)
        if(board[key] & (1 << (end_index - 1)) > 0 && (board[key] & (1 << (start_index - 1))) == 0)
            output[key] = Util.setBit(board[key],end_index,0)
            is_capture = true
            break
        end
    end
    for (key,value) in board
        if (Util.getBit(value,start_index) === UInt64(1))
            #if en passant
            if ('P' in key && ~is_capture && abs(end_index - start_index) !== UInt64(8) && abs(end_index - start_index) !== UInt64(16))
                if(colour === 'W')
                     output["BP"] = Util.setBit(board["BP"], end_index - 8, 0)
                else
                     output["WP"] = Util.setBit(board["WP"], end_index + 8,0)
                end
            #if castle
            elseif ('K' in key)
                if (move === "e1g1")
                     output["WR"] = Util.setBit(board["WR"],8,0)
                     output["WR"] = Util.setBit(output["WR"],6,1)
                elseif (move === "e1c1")
                    output["WR"] = Util.Util.setBit(board["WR"],1,0)
                    output["WR"] = Util.setBit(output["WR"],4,1)
                elseif (move === "e8g8")
                    output["BR"] = Util.setBit(board["BR"],64,0)
                    output["BR"] = Util.setBit(output["BR"],62,1)
                elseif (move === "e8c8")
                    output["BR"] = Util.setBit(board["BR"],57,0)
                    output["BR"] = Util.setBit(output["BR"],60,1)
                end
            end
            output[key] = Util.setBit(output[key],start_index,0)
            output[key] = Util.setBit(output[key],end_index,1)
        end

    end
    if(length(move) === 5) #pawn promotion
        output[colour * 'P'] = Util.setBit(board[colour * 'P'], start_index,0)
        output[colour * move[5]] = Util.setBit(board[colour * move[5]],end_index,1)
        return output
    end
    return output
end

end

max_stars_repo_name: grahamclyne/chessEngine

max_issues_repo_name: grahamclyne/chessEngine

max_forks_repo_name: grahamclyne/chessEngine


text:
using Optimize
using NLPModels

using Plots
# pyplot()
plotly()

n_min = 1
n_max = 1

# Math Prog Base collection
include("MPBProblems-1D.jl")
test_probs = mpb_probs_1D

using ScalarSolvers
solvers = [TR_Sec, TR_SecA, bissect_sec, ARC_SecA, bissect_nwt, ARC_Sec]
labels = ["TR_Sec"; "TR_SecA"; "bissect_sec"; "ARC_SecA"; "bissect_nwt"; "ARC_Sec"]

options = [Dict{Symbol,Any}(:verbose=>false)
           Dict{Symbol,Any}(:verbose=>false)
           Dict{Symbol,Any}(:verbose=>false)
           Dict{Symbol,Any}(:verbose=>false)
           Dict{Symbol,Any}(:verbose=>false)
           Dict{Symbol,Any}(:verbose=>false)
           ]

include("compare_solvers.jl")

# s1, P1, t1, Pt1 = compare_solvers_with_options2(solvers, options, labels, test_probs, n_min, n_max, printskip = false)
P2 = compare_solvers_with_options2(solvers, options, labels, test_probs, n_min, n_max, printskip = false)[2]

ext: jl


text:
@info "starting two stream instability model"

e = species.e

function TwoStream(vth2, vs1, vs2)

    function P(x,v)
        if !(v isa Array)
            v = [v]
        end

        if !(x isa Array)
            x = [x]
        end

        v = sqrt(sum(v .^2))
        x = sqrt(sum(x .^2))

        0.5/sqrt(vth2 * π) * exp(-(v-vs1)*(v-vs1)/vth2) + 0.5/sqrt(vth2 * π) * exp(-(v-vs2)*(v-vs2)/vth2) * (1+0.02*cos(3*π*x))
    end
end

D_e = Distribution(TwoStream(0.02, 1.6, -1.4), e)
G = Geometry()

plasma = ElectrostaticPlasma([D_e], G)

sol = Plasma.solve(plasma, dim=1, GPU=false, time_ub = 100.0, ub=2.0)

@save "two_stream.bson" sol

Plasma.plot(sol, 0.01)


text:
# take a sequence of matrices and multiply them in an optimal order
function fastmatmul(args::AbstractMatrix...)
    length(args) ≤ 1 && return *(args...)
    sizes = size.(args)
    if !all(sizes[i][2] == sizes[i+1][1] for i in 1:length(sizes)-1)
         throw(ArgumentError("matrix dimensions mismatch"))
    end
    partcost = Dict{Tuple{Int,Int}, Tuple{Int, Int}}()
    from, to = 1, length(sizes)
    solvemul(sizes, partcost, from, to)
    domul(args, partcost, from, to)
end

# find the optimal sequence of multiplications
function solvemul(sizes, partcost, from, to)
    if from == to
        partcost[(from, to)] = (0, from)
        return
    end
    mincost = typemax(Int)
    minj = -1
    for j in from:to-1
        haskey(partcost, (from, j)) || solvemul(sizes, partcost, from, j)
        haskey(partcost, (j+1, to)) || solvemul(sizes, partcost, j+1, to)
        curcost = sizes[from][1]*sizes[j][2]*sizes[to][2] +
                  partcost[(from, j)][1] + partcost[(j+1, to)][1]
        if curcost < mincost
            minj = j
            mincost = curcost
        end
    end
    partcost[(from, to)] = (mincost, minj)
end

# perform the multiplication given precomputed information
# about the optimal multiplication order
function domul(args, partcost, from, to)
    from == to && return args[from]
    from+1 == to && return args[from]*args[to]
    j = partcost[(from, to)][2]
    domul(args, partcost, from, j) * domul(args, partcost, j+1, to)
end

max_stars_repo_path: Chapter02/02. Fast matrix multiplication/fastmatmul.jl


text:
"""
    struct ChebyshevNodes{T} <: AbstractIntervalGrid{T}

A grid with chebyshev nodes on [-1,1].

# Example
```jldocs
julia> ChebyshevExtremae(4)
4-element ChebyshevExtremae{Float64}:
  1.0
  0.5000000000000001
 -0.4999999999999998
 -1.0
```
"""
struct ChebyshevTNodes{T} <: AbstractIntervalGrid{T}
    n   ::  Int
end

const ChebyshevNodes = ChebyshevTNodes
const ChebyshevPoints = ChebyshevTNodes

size(g::ChebyshevTNodes) = (g.n,)

unsafe_grid_getindex(g::ChebyshevTNodes{T}, i::Int) where {T} = T(-1)*cos((i-T(1)/2) * T(pi) / (g.n) )

"""
    struct ChebyshevExtremae{T} <: AbstractIntervalGrid{T}

A grid with chebyshev extrema on [-1,1].

# Example
```jldocs
julia> ChebyshevExtremae(4)
4-element ChebyshevExtremae{Float64}:
  1.0
  0.5000000000000001
 -0.4999999999999998
 -1.0
```
"""
struct ChebyshevExtremae{T} <: AbstractIntervalGrid{T}
    n   ::  Int
end

const ChebyshevPointsOfTheSecondKind = ChebyshevExtremae

size(g::ChebyshevExtremae) = (g.n,)

# TODO: flip the values so that they are sorted
unsafe_grid_getindex(g::ChebyshevExtremae{T}, i::Int) where {T} = i == 0 ? T(0) : cos((i-1)*T(pi) / (g.n-1) )

struct ChebyshevUNodes{T} <: AbstractIntervalGrid{T}
    n   :: Int
end

size(g::ChebyshevUNodes) = (g.n,)

unsafe_grid_getindex(nodes::ChebyshevUNodes{T}, i::Int) where T = cos((nodes.n + 1 - i) * convert(T,π) / (nodes.n + 1))

struct LegendreNodes{T} <: AbstractIntervalGrid{T}
    nodes   :: Vector{T}
end

size(grid::LegendreNodes) = size(grid.nodes)
unsafe_grid_getindex(grid::LegendreNodes, i::Int) = @inbounds getindex(grid.nodes, i)

struct LaguerreNodes{T} <: AbstractIntervalGrid{T}
    α       :: T
    nodes   :: Vector{T}
end

size(grid::LaguerreNodes) = size(grid.nodes)
unsafe_grid_getindex(grid::LaguerreNodes, i::Int) = @inbounds getindex(grid.nodes, i)

struct HermiteNodes{T} <: AbstractIntervalGrid{T}
    nodes   :: Vector{T}
end

size(grid::HermiteNodes) = size(grid.nodes)
unsafe_grid_getindex(grid::HermiteNodes, i::Int) = @inbounds getindex(grid.nodes, i)

struct JacobiNodes{T} <: AbstractIntervalGrid{T}
    α   ::  T
    β   ::  T
    nodes   ::Vector
end

size(grid::JacobiNodes) = size(grid.nodes)
unsafe_grid_getindex(grid::JacobiNodes, i::Int) = @inbounds getindex(grid.nodes, i)

"A vector that aliases another vector but has its own type."
abstract type VectorAlias{T} <: AbstractVector{T} end
@propagate_inbounds getindex(vector::VectorAlias, i::Int) = getindex(vector.vector, i)
size(vector::VectorAlias) = size(vector.vector)

"A function vector is a lazy vector with entries defined by a function."
abstract type FunctionVector{T} <: AbstractVector{T} end
size(vector::FunctionVector) = (vector.n,)

@propagate_inbounds function getindex(vector::FunctionVector, i::Int)
    @boundscheck (1 <= i <= length(vector)) || throw(BoundsError())
    funvector_getindex(vector, i)
end

struct ChebyshevTWeights{T,F} <: VectorAlias{T}
    vector  ::  F
end

# The Chebyshev weights are constant
function ChebyshevTWeights{T}(n) where {T}
    vector = Fill(convert(T,π) / n, n)
    ChebyshevTWeights{T,typeof(vector)}(vector)
end

struct ChebyshevUWeights{T} <: FunctionVector{T}
    n   :: Int
end

funvector_getindex(weights::ChebyshevUWeights{T}, i::Int) where {T} =
    convert(T,π)/(weights.n + 1) * sin(convert(T,weights.n + 1 -i) / (weights.n + 1) * convert(T,π))^2

gausschebyshev(::Type{T}, n::Int) where T =
    ChebyshevTNodes{T}(n), ChebyshevTWeights{T}(n)
gausschebyshev(n::Int) = gausschebyshev(Float64, n)

gausschebyshevu(::Type{T}, n::Int) where T =
    ChebyshevUNodes{T}(n), ChebyshevUWeights{T}(n)
gausschebyshevu(n::Int) = gausschebyshevu(Float64,n)

struct LegendreWeights{T} <: VectorAlias{T}
    vector  ::  Vector{T}
end

function gausslegendre(::Type{Float64}, n::Int)
    x,w = FastGaussQuadrature.gausslegendre(n)
    LegendreNodes(x), LegendreWeights(w)
end
function gausslegendre(::Type{T}, n::Int) where T
    x,w = GaussQuadrature.legendre(T, n)
    LegendreNodes(x), LegendreWeights(w)
end
gausslegendre(n::Int) = gausslegendre(Float64, n)

struct LaguerreWeights{T} <: VectorAlias{T}
    α       ::  T
    vector  ::  Vector{T}
end

function gausslaguerre(::Type{Float64}, n::Int, α::Float64)
    x,w = FastGaussQuadrature.gausslaguerre(n, α)
    LaguerreNodes(α, x), LaguerreWeights(α, w)
end
function gausslaguerre(::Type{T}, n::Int, α::T) where T
    x,w = GaussQuadrature.laguerre(n, α)
    LaguerreNodes(α, x), LaguerreWeights(α, w)
end
gausslaguerre(n::Int, α::T) where T = gausslaguerre(T, n, α)

struct HermiteWeights{T} <: VectorAlias{T}
    vector  ::  Vector{T}
end
function gausshermite(::Type{Float64}, n::Int)
    x,w = FastGaussQuadrature.gausshermite(n)
    HermiteNodes(x), HermiteWeights(w)
end
function gausshermite(::Type{T}, n::Int) where T
    x,w = GaussQuadrature.hermite(T, n)
    HermiteNodes(x), HermiteWeights(w)
end
gausshermite(n::Int) = gausshermite(Float64, n)

struct JacobiWeights{T} <: VectorAlias{T}
    α       ::  T
    β       ::  T
    vector  ::  Vector{T}
end
function gaussjacobi(::Type{Float64}, n::Int, α::Float64, β::Float64)
    x,w = FastGaussQuadrature.gaussjacobi(n, α, β)
    JacobiNodes(α, β, x), JacobiWeights(α, β, w)
end
function gaussjacobi(::Type{T}, n::Int, α::T, β::T) where T
    x,w = GaussQuadrature.jacobi(n, α, β)
    JacobiNodes(α, β, x), JacobiWeights(α, β, w)
end
gaussjacobi(n::Int, α::T, β::T) where T = gaussjacobi(T, n, α, β)

covering(::ChebyshevTNodes{T}) where T = ChebyshevInterval{T}()
name(g::ChebyshevTNodes) = "ChebyshevT nodes"

name(g::ChebyshevExtremae) = "Chebyshev extremae"
covering(::ChebyshevExtremae{T}) where T = ChebyshevInterval{T}()

name(g::ChebyshevUNodes) = "ChebyshevU nodes"
covering(::ChebyshevUNodes{T}) where T = ChebyshevInterval{T}()

name(g::LegendreNodes) = "Legendre nodes"
covering(::LegendreNodes{T}) where T = ChebyshevInterval{T}()
LegendreNodes{T}(n::Int) where T = gausslegendre(T, n)[1]

name(g::LaguerreNodes) = "Laguerre nodes α=$(g.α)"
covering(::LaguerreNodes{T}) where T = HalfLine{T}()
similargrid(grid::LaguerreNodes, T, n::Int) = LaguerreNodes{T}(n, T(grid.α))
LaguerreNodes(n::Int, α::T) where T = gausslaguerre(T, n, α)[1]
LaguerreNodes{T}(n::Int, α::T) where T = gausslaguerre(T, n, α)[1]

name(g::HermiteNodes) = "Hermite nodes"
covering(::HermiteNodes{T}) where T = DomainSets.FullSpace{T}()
HermiteNodes{T}(n::Int) where T = gausshermite(T, n)[1]

name(g::JacobiNodes) = "Jacobi nodes α=$(g.α), β=$(g.β)"
covering(::JacobiNodes{T}) where T = ChebyshevInterval{T}()
similargrid(grid::JacobiNodes, T, n::Int) = JacobiNodes{T}(n, T(grid.α), T(grid.β))
JacobiNodes{T}(n::Int, α::T, β::T) where T = gaussjacobi(T, n, α, β)[1]
JacobiNodes(n::Int, α::T, β::T) where T = gaussjacobi(T, n, α, β)[1]

# Grids with fixed support and one variable
for GRID in (:ChebyshevNodes, :ChebyshevExtremae, :ChebyshevUNodes, :LegendreNodes, :HermiteNodes)
    @eval similargrid(g::$GRID, ::Type{T}, n::Int) where {T} = $GRID{T}(n)
    @eval $GRID(n::Int) = $GRID{Float64}(n)
    # @eval $GRID(n::Int, d::AbstractInterval) =
        # $GRID(n, endpoints(d)...)
    # @eval $GRID(n::Int, a, b) = rescale($GRID{typeof((b-a)/n)}(n), a, b)
end


text:
immutable RectangularBinGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement

    function RectangularBinGeometry(
            default_statistic::Gadfly.StatisticElement=Gadfly.Stat.identity())
        new(default_statistic)
    end
end

const rectbin = RectangularBinGeometry

function histogram2d(; xbincount=nothing, xminbincount=3, xmaxbincount=150,
                       ybincount=nothing, yminbincount=3, ymaxbincount=150)
    RectangularBinGeometry(
        Gadfly.Stat.histogram2d(xbincount=xbincount,
                                xminbincount=xminbincount,
                                xmaxbincount=xmaxbincount,
                                ybincount=ybincount,
                                yminbincount=yminbincount,
                                ymaxbincount=ymaxbincount))
end

function default_statistic(geom::RectangularBinGeometry)
    geom.default_statistic
end

function element_aesthetics(::RectangularBinGeometry)
    [:x, :y, :xmin, :xmax, :ymin, :ymax, :color]
end

# Render rectangular bin (e.g., heatmap) geometry.
#
# Args:
#   geom: rectbin geometry
#   theme: the plot's theme
#   aes: some aesthetics
#
# Returns
#   A compose form.
#
function render(geom::RectangularBinGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(RGB{Float32}[theme.default_color])
    aes = inherit(aes, default_aes)

    if aes.x === nothing && (aes.xmin === nothing || aes.xmax === nothing)
        error("Geom.rectbin requires either x or both xmin and xmax be defined.")
    end

    if aes.y === nothing && (aes.ymin === nothing || aes.ymax === nothing)
        error("Geom.rectbin requires either y or both ymin and ymax be defined.")
    end

    if aes.xmin != nothing && length(aes.xmin) != length(aes.xmax)
        error("Geom.rectbin requires that xmin and xmax be of equal length.")
    end

    if aes.ymin != nothing && length(aes.ymin) != length(aes.ymax)
        error("Geom.rectbin requires that ymin and ymax be of equal length.")
    end

    nx = aes.xmin === nothing ? length(aes.x) : length(aes.xmin)
    ny = aes.ymin === nothing ? length(aes.y) : length(aes.ymin)

    if nx != ny
        error("""Geom.rectbin requires an equal number of x (or xmin/xmax) and
                 y (or ymin/ymax) values.""")
    end
    n = nx

    if aes.xmin === nothing
        xmin = [x - 0.5 for x in aes.x]
        xwidths = [1.0*cx - theme.bar_spacing]
    else
        xmin = aes.xmin
        xwidths = [(x1 - x0)*cx - theme.bar_spacing
                   for (x0, x1) in zip(aes.xmin, aes.xmax)]
    end

    if aes.ymin === nothing
        ymin = [y - 0.5 for y in aes.y]
        ywidths = [1.0*cy - theme.bar_spacing]
    else
        ymin = aes.ymin
        ywidths = [(y1 - y0)*cy - theme.bar_spacing
                   for (y0, y1) in zip(aes.ymin, aes.ymax)]
    end

    if length(aes.color) == n
        cs = aes.color
    else
        cs = Array(RGB{Float32}, n)
        for i in 1:n
            cs[i] = aes.color[((i - 1) % length(aes.color)) + 1]
        end
    end

    allvisible = true
    for c in cs
        if c == nothing
            allvisible = false
            break
        end
    end

    if !allvisible
        visibility = cs .!= nothing
        cs = cs[visibility]
        xmin = xmin[visibility]
        xmax = xmax[visibility]
        ymin = ymin[visibility]
        ymax = ymax[visibility]
    end

    return compose!(
        context(),
        rectangle(xmin, ymin, xwidths, ywidths),
        fill(cs),
        stroke(nothing),
        svgclass("geometry"),
        svgattribute("shape-rendering", "crispEdges"))
end


text:
export MPO,Matrix
using OMEinsum
using ITensors:siteinds
import Base.Matrix, ITensors.MPO

function MPO2Matrix(mpo::MPO)
    N = length(mpo)
    if N > 14
        @warn "N>14,construct Matrix may be very slow."
    end

    op = Array(mpo[1],inds(mpo[1])...)
    for j = 2:N-1
        opa = Array(mpo[j],inds(mpo[j])...)
        op = reshape(ein"aij,abpq->bipjq"(op,opa),:,2^j,2^j)
    end
    opn = Array(mpo[N],inds(mpo[N])...)
    op = ein"aij,apq->ipjq"(op,opn)
    n1,n2 = size(op)[[1;2]]
    return reshape(op,n1*n2,n1*n2)
end

Matrix(h::Hamiltonian) = MPO2Matrix(MPO(h))
MPO(h::Hamiltonian{Electron,N}) where N = MPO(h.ex.ex,siteinds("Electron",N))


text:
using Test
using QuantumOptics
using CollectiveSpins

@testset "meanfield-equations" begin

edipole = [0, 0, 1]
γ = 1
N = 2
T = [0:0.1:1;]

spinbasis = SpinBasis(1//2)
sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
sp = sigmap(spinbasis)
sm = sigmam(spinbasis)

systemgeometry = CollectiveSpins.geometry.chain(0.3, 2)
system = CollectiveSpins.SpinCollection(systemgeometry, edipole; gammas=γ)
basis = CollectiveSpins.quantum.basis(system)
I = identityoperator(spinbasis)

H = CollectiveSpins.quantum.Hamiltonian(system)
Γ, J = CollectiveSpins.quantum.JumpOperators(system)
Jdagger = [dagger(j) for j=J]
Ω = CollectiveSpins.interaction.OmegaMatrix(system)

Ψ₀ = CollectiveSpins.quantum.blochstate(0., pi/2., N)
ρ₀ = Ψ₀ ⊗ dagger(Ψ₀)

phi1 = 1.3
phi2 = 0.8

function test_meanfield(t, ρ)
    dρ_h = -1im*(H*ρ - ρ*H)
    dρ_l = DenseOperator(basis)
    for m=1:length(J), n=1:length(J)
       dρ_l += Γ[m,n]*(J[m]*ρ*Jdagger[n] - Jdagger[n]*(J[m]*ρ)/Complex(2) - ρ*Jdagger[n]*J[m]/Complex(2))
    end
    # ========== Test Hamiltonian part ==========
    dsx_master = expect(sx⊗I, dρ_h)
    dsy_master = expect(sy⊗I, dρ_h)
    dsz_master = expect(sz⊗I, dρ_h)
    dsp_master = expect(sp⊗I, dρ_h)
    dsm_master = expect(sm⊗I, dρ_h)

    # <sx>, <sy>, <sz>
    dsx_mf = Ω[1,2]*expect(sz⊗sy, ρ)
    dsy_mf = -Ω[1,2]*expect(sz⊗sx, ρ)
    dsz_mf = Ω[1,2]*(expect(sy⊗sx, ρ) - expect(sx⊗sy, ρ))
    # <s+>, <s->, <sz>
    dsp_mf = -1im*Ω[1,2]*expect(sz⊗sp, ρ)
    dsm_mf = 1im*Ω[1,2]*expect(sz⊗sm, ρ)
    dsz_mf2 = 1im*2*Ω[1,2]*expect((sm⊗sp)-(sp⊗sm), ρ)
    # <sx>rot, <sy>rot, <sz>rot
    dsx_master_rot = cos(phi1)*dsx_master - sin(phi1)*dsy_master
    dsy_master_rot = sin(phi1)*dsx_master + cos(phi1)*dsy_master
    dsz_master_rot = dsz_master
    dsx_rot = (Ω[1,2]*sin(phi1-phi2)*expect(sz⊗(cos(phi2)*sx-sin(phi2)*sy), ρ)
                + Ω[1,2]*cos(phi1-phi2)*expect(sz⊗(sin(phi2)*sx+cos(phi2)*sy), ρ))
    dsy_rot = (-Ω[1,2]*cos(phi1-phi2)*expect(sz⊗(cos(phi2)*sx-sin(phi2)*sy), ρ)
                + Ω[1,2]*sin(phi1-phi2)*expect(sz⊗(sin(phi2)*sx+cos(phi2)*sy), ρ))
    dsz_rot = (-Ω[1,2]*sin(phi1-phi2)*(expect((cos(phi1)*sx-sin(phi1)*sy)⊗(cos(phi2)*sx-sin(phi2)*sy), ρ)
                                    + expect((sin(phi1)*sx+cos(phi1)*sy)⊗(sin(phi2)*sx+cos(phi2)*sy), ρ))
                +Ω[1,2]*cos(phi1-phi2)*(expect((sin(phi1)*sx+cos(phi1)*sy)⊗(cos(phi2)*sx-sin(phi2)*sy), ρ)
                                    - expect((cos(phi1)*sx-sin(phi1)*sy)⊗(sin(phi2)*sx+cos(phi2)*sy), ρ))
        )
    @test abs(dsx_master - dsx_mf) < 1e-12
    @test abs(dsy_master - dsy_mf) < 1e-12
    @test abs(dsz_master - dsz_mf) < 1e-12
    @test abs(dsp_master - dsp_mf) < 1e-12
    @test abs(dsm_master - dsm_mf) < 1e-12
    @test abs(dsz_master - dsz_mf2) < 1e-12
    @test abs(dsx_master_rot - dsx_rot) <1e-12
    @test abs(dsy_master_rot - dsy_rot) <1e-12
    @test abs(dsz_master_rot - dsz_rot) <1e-12

    # ========== Test Lindblad part ==========
    dsx_master = expect(sx⊗I, dρ_l)
    dsy_master = expect(sy⊗I, dρ_l)
    dsz_master = expect(sz⊗I, dρ_l)
    dsp_master = expect(sp⊗I, dρ_l)
    dsm_master = expect(sm⊗I, dρ_l)

    # <sx>, <sy>, <sz>
    dsx_mf = -0.5*expect(sx⊗I, ρ) + 0.5*Γ[1,2]*expect(sz⊗sx, ρ)
    dsy_mf = -0.5*expect(sy⊗I, ρ) + 0.5*Γ[1,2]*expect(sz⊗sy, ρ)
    dsz_mf = -(1+expect(sz⊗I, ρ)) - 0.5*Γ[1,2]*(expect(sx⊗sx, ρ) + expect(sy⊗sy, ρ))
    # <s+>, <s->, <sz>
    dsp_mf = -0.5*expect(sp⊗I, ρ) + 0.5*Γ[1,2]*expect(sz⊗sp, ρ)
    dsm_mf = -0.5*expect(sm⊗I, ρ) + 0.5*Γ[1,2]*expect(sz⊗sm, ρ)
    dsz_mf2 = -(1+expect(sz⊗I, ρ)) - Γ[1,2]*(expect(sp⊗sm, ρ) + expect(sm⊗sp, ρ))
    # <sx>rot, <sy>rot, <sz>rot
    dsx_master_rot = cos(phi1)*dsx_master - sin(phi1)*dsy_master
    dsy_master_rot = sin(phi1)*dsx_master + cos(phi1)*dsy_master
    dsz_master_rot = dsz_master
    dsx_rot = (-0.5*expect((cos(phi1)*sx-sin(phi1)*sy)⊗I, ρ)
            +0.5*Γ[1,2]*cos(phi1-phi2)*expect(sz⊗(cos(phi2)*sx-sin(phi2)*sy), ρ)
            -0.5*Γ[1,2]*sin(phi1-phi2)*expect(sz⊗(sin(phi2)*sx+cos(phi2)*sy), ρ))
    dsy_rot = (-0.5*expect((sin(phi1)*sx+cos(phi1)*sy)⊗I, ρ)
            +0.5*Γ[1,2]*sin(phi1-phi2)*expect(sz⊗(cos(phi2)*sx-sin(phi2)*sy), ρ)
            +0.5*Γ[1,2]*cos(phi1-phi2)*expect(sz⊗(sin(phi2)*sx+cos(phi2)*sy), ρ))
    dsz_rot = (-(1+expect(sz⊗I, ρ))
            -0.5*Γ[1,2]*cos(phi1-phi2)*(expect((cos(phi1)*sx-sin(phi1)*sy)⊗(cos(phi2)*sx-sin(phi2)*sy), ρ)
                                    + expect((sin(phi1)*sx+cos(phi1)*sy)⊗(sin(phi2)*sx+cos(phi2)*sy), ρ))
            -0.5*Γ[1,2]*sin(phi1-phi2)*(expect((sin(phi1)*sx+cos(phi1)*sy)⊗(cos(phi2)*sx-sin(phi2)*sy), ρ)
                                    - expect((cos(phi1)*sx-sin(phi1)*sy)⊗(sin(phi2)*sx+cos(phi2)*sy), ρ)))

    @test abs(dsx_master - dsx_mf) < 1e-12
    @test abs(dsy_master - dsy_mf) < 1e-12
    @test abs(dsz_master - dsz_mf) < 1e-12
    @test abs(dsp_master - dsp_mf) < 1e-12
    @test abs(dsm_master - dsm_mf) < 1e-12
    @test abs(dsz_master - dsz_mf2) < 1e-12
    @test abs(dsx_master_rot - dsx_rot) <1e-12
    @test abs(dsy_master_rot - dsy_rot) <1e-12
    @test abs(dsz_master_rot - dsz_rot) <1e-12
end

CollectiveSpins.quantum.timeevolution(T, system, ρ₀; fout=test_meanfield)

end # testset


text:
"""

```
remove_drift!(v)
```

Removes possible drift from velocities.

"""
function remove_drift!(v)
  vmean = mean(v)
  for i in 1:length(v)
    v[i] -= vmean
  end
end


text:
# ------------------------------------------------------------------------------------------
# # Data structures
#
# Once we start working with many pieces of data at once, it will be convenient for us to
# store data in structures like arrays or dictionaries (rather than just relying on
# variables).

#
# Types of data structures covered:
# 1. Tuples
# 2. Dictionaries
# 3. Arrays
#
#

# As an overview, tuples and arrays are both ordered sequences of elements (so we can index
# into them). Dictionaries and arrays are both mutable.
# We'll explain this more below!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Tuples
#
# We can create a tuple by enclosing an ordered collection of elements in `( )`.
#
# Syntax:

# ```julia
# (item1, item2, ...)```
# ------------------------------------------------------------------------------------------

myfavoriteanimals = ("penguins", "cats", "sugargliders")

# ------------------------------------------------------------------------------------------
# We can index into this tuple,
# ------------------------------------------------------------------------------------------

myfavoriteanimals[1]

# ------------------------------------------------------------------------------------------
# but since tuples are immutable, we can't update it
# ------------------------------------------------------------------------------------------

myfavoriteanimals[1] = "otters"

# ------------------------------------------------------------------------------------------
# ## Now in 1.0: NamedTuples
#
# As you might guess, `NamedTuple`s are just like `Tuple`s except that each element
# additionally has a name! They have a special syntax using `=` inside a tuple:
#
# ```julia
# (name1 = item1, name2 = item2, ...)
# ```
# ------------------------------------------------------------------------------------------

myfavoriteanimals = (bird = "penguins", mammal = "cats", marsupial = "sugargliders")

# ------------------------------------------------------------------------------------------
# Like regular `Tuples`, `NamedTuples` are ordered, so we can retrieve their elements via
# indexing:
# ------------------------------------------------------------------------------------------

myfavoriteanimals[1]

# ------------------------------------------------------------------------------------------
# They also add the special ability to access values by their name:
# ------------------------------------------------------------------------------------------

myfavoriteanimals.bird

# ------------------------------------------------------------------------------------------
# ## Dictionaries
#
# If we have sets of data related to one another, we may choose to store that data in a
# dictionary. We can create a dictionary using the `Dict()` function, which we can
# initialize as an empty dictionary or one storing key, value pairs.
#
# Syntax:
# ```julia
# Dict(key1 => value1, key2 => value2, ...)```
#
# A good example is a contacts list, where we associate names with phone numbers.
# ------------------------------------------------------------------------------------------

myphonebook = Dict("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

# ------------------------------------------------------------------------------------------
# In this example, each name and number is a "key" and "value" pair. We can grab Jenny's
# number (a value) using the associated key
# ------------------------------------------------------------------------------------------

myphonebook["Jenny"]

# ------------------------------------------------------------------------------------------
# We can add another entry to this dictionary as follows
# ------------------------------------------------------------------------------------------

myphonebook["Kramer"] = "555-FILK"

# ------------------------------------------------------------------------------------------
# Let's check what our phonebook looks like now...
# ------------------------------------------------------------------------------------------

myphonebook

# ------------------------------------------------------------------------------------------
# We can delete Kramer from our contact list - and simultaneously grab his number - by using
# `pop!`
# ------------------------------------------------------------------------------------------

pop!(myphonebook, "Kramer")

myphonebook

# ------------------------------------------------------------------------------------------
# Unlike tuples and arrays, dictionaries are not ordered. So, we can't index into them.
# ------------------------------------------------------------------------------------------

myphonebook[1]

# ------------------------------------------------------------------------------------------
# In the example above, `julia` thinks you're trying to access a value associated with the
# key `1`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Arrays
#
# Unlike tuples, arrays are mutable. Unlike dictionaries, arrays contain ordered
# collections.

# We can create an array by enclosing this collection in `[ ]`.
#
# Syntax:

# ```julia
# [item1, item2, ...]```
#
#
# For example, we might create an array to keep track of my friends
# ------------------------------------------------------------------------------------------

myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

# ------------------------------------------------------------------------------------------
# The `1` in `Array{String,1}` means this is a one dimensional vector.  An `Array{String,2}`
# would be a 2d matrix, etc.  The `String` is the type of each element.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# or to store a sequence of numbers
# ------------------------------------------------------------------------------------------

fibonacci = [1, 1, 2, 3, 5, 8, 13]

mixture = [1, 1, 2, 3, "Ted", "Robyn"]

# ------------------------------------------------------------------------------------------
# Once we have an array, we can grab individual pieces of data from inside that array by
# indexing into the array. For example, if we want the third friend listed in `myfriends`,
# we write
# ------------------------------------------------------------------------------------------

myfriends[3]

# ------------------------------------------------------------------------------------------
# We can use indexing to edit an existing element of an array
# ------------------------------------------------------------------------------------------

myfriends[3] = "Baby Bop"

# ------------------------------------------------------------------------------------------
# Yes, Julia is 1-based indexing, not 0-based like Python.  Wars are fought over lesser
# issues. I have a friend with the wisdom of Solomon who proposes settling this once and for
# all with ½ 😃
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can also edit the array by using the `push!` and `pop!` functions. `push!` adds an
# element to the end of an array and `pop!` removes the last element of an array.
#
# We can add another number to our fibonnaci sequence
# ------------------------------------------------------------------------------------------

push!(fibonacci, 21)

# ------------------------------------------------------------------------------------------
# and then remove it
# ------------------------------------------------------------------------------------------

pop!(fibonacci)

fibonacci

# ------------------------------------------------------------------------------------------
# So far I've given examples of only 1D arrays of scalars, but arrays can have an arbitrary
# number of dimensions and can also store other arrays.
#

# For example, the following are arrays of arrays:
# ------------------------------------------------------------------------------------------

favorites = [["koobideh", "chocolate", "eggs"],["penguins", "cats", "sugargliders"]]

numbers = [[1, 2, 3], [4, 5], [6, 7, 8, 9]]

# ------------------------------------------------------------------------------------------
# Below are examples of 2D and 3D arrays populated with random values.
# ------------------------------------------------------------------------------------------

rand(4, 3)

rand(4, 3, 2)

# ------------------------------------------------------------------------------------------
# Be careful when you want to copy arrays!
# ------------------------------------------------------------------------------------------

fibonacci

somenumbers = fibonacci

somenumbers[1] = 404

fibonacci

# ------------------------------------------------------------------------------------------
# Editing `somenumbers` caused `fibonacci` to get updated as well!
#
# In the above example, we didn't actually make a copy of `fibonacci`. We just created a new
# way to access the entries in the array bound to `fibonacci`.
#
# If we'd like to make a copy of the array bound to `fibonacci`, we can use the `copy`
# function.
# ------------------------------------------------------------------------------------------

# First, restore fibonacci
fibonacci[1] = 1
fibonacci

somemorenumbers = copy(fibonacci)

somemorenumbers[1] = 404

fibonacci

# ------------------------------------------------------------------------------------------
# In this last example, fibonacci was not updated. Therefore we see that the arrays bound to
# `somemorenumbers` and `fibonacci` are distinct.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Exercises
#
# #### 3.1
# Create an array, `a_ray`, with the following code:
#
# ```julia
# a_ray = [1, 2, 3]
# ```
#
# Add the number `4` to the end of this array and then remove it.
# ------------------------------------------------------------------------------------------

@assert a_ray == [1, 2, 3]

# ------------------------------------------------------------------------------------------
# #### 3.2
# Try to add "Emergency" as key to `myphonebook` with the value `string(911)` with the
# following code
# ```julia
# myphonebook["Emergency"] = 911
# ```
#
# Why doesn't this work?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### 3.3
# Create a new dictionary called `flexible_phonebook` that has Jenny's number stored as an
# integer and Ghostbusters' number stored as a string with the following code
#
# ```julia
# flexible_phonebook = Dict("Jenny" => 8675309, "Ghostbusters" => "555-2368")
# ```
# ------------------------------------------------------------------------------------------

@assert flexible_phonebook == Dict("Jenny" => 8675309, "Ghostbusters" => "555-2368")

# ------------------------------------------------------------------------------------------
# #### 3.4
# Add the key "Emergency" with the value `911` (an integer) to `flexible_phonebook`.
# ------------------------------------------------------------------------------------------

@assert haskey(flexible_phonebook, "Emergency")

@assert flexible_phonebook["Emergency"] == 911

# ------------------------------------------------------------------------------------------
# #### 3.5
# Why can we add an integer as a value to `flexible_phonebook` but not `myphonebook`? How
# could we have initialized `myphonebook` so that it would accept integers as values?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Please click on `Validate` button at the top, upon completion of the exercise
# ------------------------------------------------------------------------------------------


text:
# new Sandshark example
# add more julia processes
using Distributed
addprocs(4)

using Caesar, RoME
@everywhere using Caesar, RoME
using Interpolations
using Distributions

using RoMEPlotting
using Gadfly, DataFrames
using ProgressMeter
using DelimitedFiles

# const TU = TransformUtils

Gadfly.set_default_plot_size(35cm,25cm)

include(joinpath(@__DIR__,"Plotting.jl"))
include(joinpath(@__DIR__,"SandsharkUtils.jl"))

# Step: Selecting a subset for processing and build up a cache of the factors.
epochs = timestamps[50:3:500]
lastepoch = 0
for ep in epochs
  global lastepoch
  if lastepoch != 0
    # @show interp_yaw(ep)
    deltaAng = interp_yaw(ep) - interp_yaw(lastepoch)

    wXi = TU.SE2([interp_x(lastepoch);interp_y(lastepoch);interp_yaw(lastepoch)])
    wXj = TU.SE2([interp_x(ep);interp_y(ep);interp_yaw(ep)])
    iDXj = se2vee(wXi\wXj)
    NAV[ep] = iDXj
    # println("$(iDXj[1]), $(iDXj[2]), $(iDXj[3])")

    odoDict[ep] = Pose2Pose2(MvNormal(NAV[ep], Matrix(Diagonal([0.1;0.1;0.005].^2))))
  end
  rangepts = rangedata[ep][:]
  rangeprob = kde!(rangepts)
  azipts = azidata[ep][:,1]
  aziprob = kde!(azipts)

  # prep the factor functions
  ppbrDict[ep] = Pose2Point2BearingRange(aziprob, rangeprob)
  lastepoch = ep
end

## Step: Building the factor graph
fg = initfg()
# Add a central beacon with a prior
addVariable!(fg, :l1, Point2)
# Pinger location is (0.6; -16)
addFactor!(fg, [:l1], PriorPose2( MvNormal([0.6; -16], Matrix(Diagonal([0.1; 0.1].^2)) ) ), autoinit=true)

# init tree for simpler code later down
tree = wipeBuildNewTree!(fg)

## and some settings

# enable fixed lag operation
getSolverParams(fg).isfixedlag = true
getSolverParams(fg).qfl = 20

# first solve and initialization
getSolverParams(fg).drawtree = true
# getSolverParams(fg).showtree = false

index = 0

for STEP in 40:10:50
    global tree
    for ep in epochs[(STEP+1):(STEP+10)]
      global index
      curvar = Symbol("x$index")
      addVariable!(fg, curvar, Pose2)

      # xi -> l1 - nonparametric factor
      if index % 5 == 0
          addFactor!(fg, [curvar; :l1], ppbrDict[ep], autoinit=true)
      end

      if ep != epochs[1]
        # Odo factor x(i-1) -> xi
        addFactor!(fg, [Symbol("x$(index-1)"); curvar], odoDict[ep], autoinit=true)
      else
        # Prior to the first pose location (a "GPS" prior)
        initLoc = [interp_x(ep);interp_y(ep);interp_yaw(ep)]
        println("Adding a prior at $curvar, $initLoc")
        addFactor!(fg, [curvar], PriorPose2( MvNormal(initLoc, Matrix(Diagonal([0.1;0.1;0.05].^2))) ), autoinit=true)
      end
      # Heading partial prior
      addFactor!(fg, [curvar], RoME.PartialPriorYawPose2(Normal(interp_yaw(ep), deg2rad(3))), autoinit=true)
      index+=1
    end
    tree, smt, hist = solveTree!(fg, tree)
    pla = drawPosesLandmarksAndOdo(fg, ppbrDict, navkeys, X, Y, lblkeys, lblX, lblY)
    pla |> PDF(joinpath(getSolverParams(fg).logpath, "sandshark-beacon_$STEP.pdf"))
end

# drawGraph(fg)

drawPosesLandms(fg)

# pl = drawPoses(fg, spscale=2.75) # Just for odo plot
# Roll again for inspiration check
## PLOT BEAM PATTERNS
# Gadfly.push_theme(:default)

Gadfly.draw(PDF("sandshark-beacon_$t.pdf", 12cm, 15cm), pla)
Gadfly.draw(PNG("sandshark-beacon_$t.png", 12cm, 15cm), pla)

####  DEBUG PPBRDict ====================

LL = [0.6; -16; 0]
function expectedLocalBearing(curPose, LL)
    v = [LL[1] - curPose[1]; LL[2] - curPose[2]]
    world = atan2(v[2], v[1])
    loc = rad2deg(TU.wrapRad(world - curPose[3]))
    while loc < 0
        loc += 360
    end
    return loc
end
# println("$(interp_x[ep]), $(interp_y[ep]), world yaw = $(interp_yaw[ep] * 180 / pi)")
# XX = [interp_x(ep); interp_y[ep]; interp_yaw[ep]]
# expectedLocalBearing(XX, LL)
# println("Expected local yaw = $(expectedLocalBearing(XX, LL)). Sonar local yaw = $(rad2deg(getKDEMax(ppbr.bearing)[1]))")
# xyt = se2vee(SE2(XX[1:3]) \ SE2(LL))
# bear = rad2deg(TU.wrapRad(atan2(xyt[2],xyt[1]) -XX[3]))

bearExp = []
bearAcoustics = []
for ep in epochs
    XX = [interp_x(ep); interp_y[ep]; interp_yaw[ep]]
    ppbr = ppbrDict[ep]
    push!(bearExp, expectedLocalBearing(XX, LL))
    push!(bearAcoustics, rad2deg(getKDEMax(ppbr.bearing)[1]))
    println("local beacon heading: Expected = $(expectedLocalYaw(XX, LL)). Sonar = $(rad2deg(getKDEMax(ppbr.bearing)[1]))")
end

layers = []
push!(layers, Gadfly.layer(x=1:length(bearExp), y=bearExp, Geom.path())[1], Theme(default_color=color("red")))
push!(layers, Gadfly.layer(x=1:length(bearAcoustics), y=bearAcoustics, Geom.path(), Theme(default_color=color("blue"))))
Gadfly.plot(layers...)
# end

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

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

# pl = drawPoses(fg, spscale=2.75)
# drawPosesLandms(fg, spscale=0.75) #Means so we don't run into MM == Union() || Dict{} in
# You rolled 20!

# Roll again for debug check
## Debugging strange headings
# posePts = get2DPoseMeans(fg)
# landPts = get2DLandmMeans(fg)
## Debugging landmark bearing range

# ppbrDict[epoch_slice[1]]
# getSample(ppbrDict[epoch_slice[1]],100)
#
#
# addVariable!(fg, :l1, Point2)
# addFactor!(fg, [:x0; :l1], ppbrDict[epoch_slice[1]], autoinit=true)
#
# ls(fg, :l1)
# pts = IIF.approxConv(fg, :x0l1f1, :l1)
#
# fct = getData(getVert(fg, :x0l1f1, nt=:fnc))
#
# fct.fnc.zDim

# dev
# ep = epochs[1]
# # azidata[ep][:,1]
# pll = layerBeamPatternRose(ppbrDict[ep].bearing)
# Gadfly.plot(pll...)
#
# pll = layerBeamPatternRose(ppbrDict[ep].bearing, wRr=TU.R(pi/2))
# Gadfly.plot(pll...)

#
# diff(epochs[[end;1]])

## GADFLY EXAMPLE
# plotLocalProduct(fg, :x49, dims=[1;2])
# stuff = IncrementalInference.localProduct(fg, :x49)
# plotKDE(stuff[2], dims=[1;2], levels=10)
#
# # STUFF
# fsym = :x49l1f1
# const TU = TransformUtils
# XX, LL = (KDE.getKDEMax.(IIF.getBelief.(fg, IIF.lsf(fg, fsym)))...)
# @show xyt = se2vee(SE2(XX[1:3]) \ SE2([LL[1:2];0.0]))
# bear= rad2deg(TU.wrapRad(atan2(-xyt[2],xyt[1]) -XX[3]))
# b = IncrementalInference.getData(fg, fsym, nt=:fnc).fnc.usrfnc!
# b.bearing
# p = plotKDE(b.range)
# Gadfly.draw(PDF("multimodal_bearing_range.pdf", 12cm, 15cm), p)
#
# # SAVE THIS PLOT
# # epochs = timestamps[11:2:200]
# g = getBelief.(fg, [:l1, :x49])
# m1 = KDE.marginal(g[1], [1;2])
# m2 = KDE.marginal(g[2], [1;2])
# norm(diff(KDE.getKDEMax.([m1; m2])))
#
# kde!(rand(Normal(4.0, 0.1),100))
#
# ###############
#
# PL = []
#
# push!(PL, Gadfly.layer(x=collect(1:10), y=randn(10), Geom.path, Theme(default_color=colorant"red"))[1])
# push!(PL, Gadfly.layer(x=collect(1:10), y=randn(10), Geom.path, Theme(default_color=colorant"green"))[1])
# push!(PL, Coord.Cartesian(xmin=-1.0,xmax=11, ymin=-5.0, ymax=5.0))
# push!(PL, Guide.ylabel("test"))
# push!(PL, Guide.ylabel("test"))
# pl = drawPoses(fg)
# push!(pl.layers, )
#
# pl = drawPoses(fg)
# pl = Gadfly.plot(pl.layers...)
#
#
#
# pl = Gadfly.plot(PL...)
#
# @show fieldnames(pl)
# push!(pl.layers, Gadfly.layer(x=collect(1:10), y=randn(10), Geom.path, Theme(default_color=colorant"magenta"))[1] )
#
# pl
# # vstasck, hstack
#
#
# Gadfly.push_theme(:default)
# # SVG, PDF, PNG
# Gadfly.draw(PDF("/tmp/testfig.pdf", 12cm, 8cm), pl)
# run(`evince /tmp/testfig.pdf`)
#
#
# ### PATH + Poses
#
# pl = drawPosesLandms(fg, spscale=2.5)
# # Add odo
# navdf = DataFrame(
#   ts = navkeys,
#   x = X,
#   y = Y
# )
# # pl = Gadfly.layer(navdf, x=:x, y=:y, Geom.path())
# push!(pl.layers, Gadfly.layer(navdf, x=:x, y=:y, Geom.path())[1])
# Gadfly.plot(pl.layers)

fsym = :x49l1f1
const TU = TransformUtils
XX, LL = (KDE.getKDEMax.(IIF.getBelief.(fg, IIF.lsf(fg, fsym)))...)
@show xyt = se2vee(SE2(XX[1:3]) \ SE2([LL[1:2];0.0]))
bear= rad2deg(TU.wrapRad(atan2(-xyt[2],xyt[1]) -XX[3]))

#


text:
using RationalRoots, Test

bitstypes  = (Int8, Int16, Int32, Int64, Int128)
inttypes = (bitstypes..., BigInt)
@testset "Constructor" begin
    r = rand(Float64)
    x = RationalRoot{Int}(Val{:inner}(), -3//5)
    for T in bitstypes
        @test (@inferred RationalRoot{T}(1)).signedsquare === Rational{T}(1)
        @test (@inferred RationalRoot{T}(-1)).signedsquare === Rational{T}(-1)
        @test (@inferred RationalRoot{T}(2)).signedsquare === Rational{T}(4)
        @test (@inferred RationalRoot{T}(-2)).signedsquare === Rational{T}(-4)
        @test (@inferred RationalRoot{T}(1//2)).signedsquare === Rational{T}(1//4)
        @test (@inferred RationalRoot{T}(-1//2)).signedsquare === Rational{T}(-1//4)
        @test (@inferred RationalRoot{T}(r)).signedsquare ===
                                                rationalize(T, signedsquare(r))
        @test (@inferred RationalRoot{T}(x)).signedsquare === Rational{T}(x.signedsquare)

        y = RationalRoot{T}(Val{:inner}(), 0.4)
        @test RationalRoot(y) === y
    end
    T = BigInt
    s = (@inferred RationalRoot{T}(1)).signedsquare
    @test s == convert(Rational{T}, 1) && typeof(s) == Rational{T}
    s = (@inferred RationalRoot{T}(-1)).signedsquare
    @test s == convert(Rational{T}, -1) && typeof(s) == Rational{T}
    s = (@inferred RationalRoot{T}(2)).signedsquare
    @test s == convert(Rational{T}, 4) && typeof(s) == Rational{T}
    s = (@inferred RationalRoot{T}(-2)).signedsquare
    @test s == convert(Rational{T}, -4) && typeof(s) == Rational{T}
    s = (@inferred RationalRoot{T}(1//2)).signedsquare
    @test s == convert(Rational{T}, 1//4) && typeof(s) == Rational{T}
    s = (@inferred RationalRoot{T}(-1//2)).signedsquare
    @test s == convert(Rational{T}, -1//4) && typeof(s) == Rational{T}
    s = (@inferred RationalRoot{T}(r)).signedsquare
    @test s == rationalize(T, signedsquare(r)) && typeof(s) == Rational{T}
    s = (@inferred RationalRoot{T}(x)).signedsquare
    @test s == x.signedsquare && typeof(s) == Rational{T}

    y = RationalRoot{T}(Val{:inner}(), 0.4)
    @test RationalRoot(y) === y

    for T in (Float32, Float64, BigFloat)
        @test (@inferred RationalRoot(T(-1//4))).signedsquare === -1//16
    end

    for T in inttypes
        @test typeof(@inferred Rational(one(T))) == Rational{T}
    end
end

@testset "signedroot and signedsquare" begin
    for T in bitstypes
        @test (@inferred signedsquare(T(2))) === T(4)
        @test (@inferred signedroot(T(4))) === RationalRoot{T}(2)
        @test (@inferred signedsquare(Rational{T}(-2//3))) === Rational{T}(-4//9)
        @test (@inferred signedroot(Rational{T}(-4//9))) === RationalRoot{T}(-2//3)
    end
    @test (@inferred signedsquare(2.5)) === 6.25
    @test (@inferred signedroot(6.25)) === 2.5
    @test (@inferred signedroot(RationalRoot, 6.25)).signedsquare ===
        Rational{Int}(625//100)
    @test (@inferred signedsquare(-1.5f0)) === -2.25f0
    @test (@inferred signedroot(RationalRoot, -1.5f0)).signedsquare ===
        Rational{Int}(-15//10)
    @test (@inferred signedsquare(1/big(3))) ≈ 1/big(9)
end

@testset "algebra" begin
    for T in bitstypes
        a = -rand(T(1):T(5)) // rand(T(1):T(5))
        b = rand(T(1):T(5)) // rand(T(1):T(5))

        for op in (+, -, inv)
            @test (@inferred op(signedroot(a))) === signedroot(op(a))
        end

        for op in (*, /, \, //)
            @test (@inferred op(signedroot(a),signedroot(b))) === signedroot(op(a,b))
            @test (@inferred op(signedroot(a),b)) === signedroot(op(a,signedsquare(b)))
            @test (@inferred op(a,signedroot(b))) === signedroot(op(signedsquare(a),b))
        end

        @test (@inferred sign(signedroot(a))) === sign(a)*one(RationalRoot{T})

        x = signedroot(a)
        @test one(x) * x === x
        @test zero(x) === zero(typeof(x))
        @test one(x) === one(typeof(x))
        @test isone(one(x))
        @test !iszero(one(x))
        @test iszero(zero(x))
        @test !isone(zero(x))
        @test zero(x) * x === zero(x)
    end
end

@testset "conversion and promotion" begin
    for T in inttypes
        x = signedroot(T(1)//T(3))
        for T2 in inttypes
            @test promote_type(RationalRoot{T}, T2) == RationalRoot{promote_type(T,T2)}
            @test promote_type(RationalRoot{T}, Rational{T2}) ==
                                                        RationalRoot{promote_type(T,T2)}

            @test typeof(x*one(T2)) == RationalRoot{promote_type(T,T2)}
            @test typeof(x*one(Rational{T2})) == RationalRoot{promote_type(T,T2)}

            @test typeof(x+one(T2)) == Float64
            @test typeof(x-one(Rational{T2})) == Float64
        end
        for T2 in (Float32, Float64, BigFloat)
            @test promote_type(RationalRoot{T}, T2) == T2

            @test typeof(x+one(T2)) == T2
            @test typeof(x-one(T2)) == T2
            @test typeof(x*one(T2)) == T2
            @test typeof(x/one(T2)) == T2
        end
    end

    for T in inttypes
        @test isinteger(signedroot(T(-4)))
        @test !isinteger(signedroot(T(3)))
        @test signedroot(RationalRoot, T(4)) == 2
        @test signedroot(RationalRoot{T}, 4) == 2
        @test signedroot(RationalRoot{T}, -4//9) == -2//3

        h = rand(UInt)
        @test hash(signedroot(RationalRoot{T}, 4), h) == hash(2, h)
        @test hash(signedroot(RationalRoot{T}, -4//9), h) == hash(-2//3, h)
        @test hash(signedroot(1//3)) != hash(sqrt(1/3))

        @test typemax(RationalRoot{T}) == typemax(Rational{T})
        @test typemin(RationalRoot{T}) == typemin(Rational{T})
    end
    for T in (Float32, Float64, BigFloat)
        @test signedroot(RationalRoot, T(4.0)) == T(2.0)
        @test signedroot(RationalRoot, T(2.0)) ≈ sqrt(T(2.0))
        @test signedroot(RationalRoot, T(2.0)) != sqrt(T(2.0))
        @test signedroot(RationalRoot{BigInt}, one(T)/T(4.0)) == one(T)/T(2.0)
        @test signedroot(RationalRoot{BigInt}, -one(T)/T(3.0)) ≈ -sqrt(one(T)/T(3.0))
        @test signedroot(RationalRoot{BigInt}, one(T)/T(3.0)) != sqrt(one(T)/T(3.0))
    end

    for T in bitstypes
        @test convert(T, signedroot(4)) === T(2)
        @test convert(T, signedroot(T(4))) === T(2)
        @test_throws InexactError convert(T, signedroot(3))

        @test convert(Rational{T}, signedroot(4//9)) === Rational{T}(2//3)
        @test_throws InexactError convert(Rational{T}, signedroot(3//2))

        for T2 in (Float32, Float64, BigFloat)
            @test T2(signedroot(T(-2))) ≈ -sqrt(T2(T(2)))
        end
    end
    T = BigInt
    a = big(typemax(Int64))^2
    b = a + 1
    @test convert(T, signedroot(a)) == typemax(Int64)
    @test convert(Rational{T}, signedroot(a)//b) == typemax(Int64)//b
    for T in (Float32, Float64, BigFloat)
        @test T( signedroot(a)//b ) ≈ T( typemax(Int64)/b )
    end

    x = big(signedroot(-4//9))
    @test x == -2//3 && typeof(signedsquare(x)) == Rational{BigInt}

    for T in (Int8, Int16, Int32, Int64)
        T2 = widen(T)
        @test widen(RationalRoot{T}) == RationalRoot{T2}
        @test widen(signedroot(T(3)//T(5))) === signedroot(T2(3)//T2(5))
    end
end

@testset "comparison" begin
    for T in inttypes
        @test RationalRoot{T}(3//4) <= 3//4
        @test RationalRoot{T}(3//4) <= 3/4
        @test 3//4 <= RationalRoot{T}(3//4)
        @test 3/4 <= RationalRoot{T}(3//4)
        @test 2//5 <= RationalRoot{T}(3//4) <= 4//5
        @test 4//5 >= RationalRoot{T}(3//4) >= 2//5
        @test 0 <= RationalRoot{T}(3//5) <= 1
        @test 1 >= RationalRoot{T}(3//5) >= 0
        @test 2//5 < RationalRoot{T}(3//5) < 4//5
        @test 4//5 > RationalRoot{T}(3//5) > 2//5
        @test 0 < RationalRoot{T}(3//5) < 1
        @test 1 > RationalRoot{T}(3//5) > 0
        @test -Inf < RationalRoot{T}(3//5) < Inf
        @test typemin(RationalRoot{T}) < RationalRoot{T}(3//5) < typemax(RationalRoot{T})
        @test typemin(Rational{T}) < RationalRoot{T}(3//5) < typemax(Rational{T})
        if isbitstype(T)
            @test typemin(T) < RationalRoot{T}(3//5) < typemax(T)
        end
    end
end

@testset "show" begin
    @test sprint(show, signedroot(+1//3)) == "+√(1//3)"
    @test sprint(show, signedroot(-2//5)) == "-√(2//5)"
end


text:
# SCSG implementation

using Random

export SCSG,
        minimize,
        BatchSizeInfo,
        get_batchsizeinfo

include("evaluator.jl")

"""
# SCSG

This is the structure describing parameters used in Stochastic Controlled Stochastic Gradient

## Fields

- eta_zero : initial step size
- m_zero : fraction of terms to consider in initialisation of the number of mini batch to run for each large batch (mⱼ)
- mini\\_batch\\_size\\_init : governs size of mini batches (1 is a good initialization)
- large\\_batch\\_size\\_init : fraction of terms to consider in initialisation of large batch size

"""
mutable struct SCSG
    # initial step size
    eta_zero :: Float64
    # fraction of nbterms to consider in initialization of mⱼ governing evolution of nb small iterations
    m_zero :: Float64
    # b_0 in the paper
    mini_batch_size_init :: Int64
    # related to B_0 in Paper. Fraction of terms to consider in initialisation of B_0
    large_batch_size_init :: Float64
end

"""
# struct BatchSizeInfo

stores all parameters describing one batch characteristics

"""
mutable struct BatchSizeInfo
    iteration :: Int64
    large_batchsize:: Int64
    mini_batchsize:: Int64
    # parameter used in determining number of mini batch
    nbminibatchparam :: Float64
    # step used in position update
    stepsize :: Float64
end

"""
returns a struct BatchSizeInfo for a given iteration
"""
function get_batchsizeinfo(scsg_pb :: SCSG, batch_growing_factor :: Float64 , nbterms::Int64, iteration :: Int64)
    alfa_j = batch_growing_factor^iteration
    # max size of large batch is 100 or 0.1 * the number of terms
    max_large_batch_size = nbterms > 100 ? ceil(Int64, nbterms/10.) : nbterms
    # ensure max_mini_batch_size is at least 1.
    max_mini_batch_size = ceil(Int64, nbterms/100.)
    # B_j
    large_batch_size = min(ceil(Int64, scsg_pb.large_batch_size_init * nbterms * alfa_j * alfa_j), max_large_batch_size)
    # b_j  grow slowly
    mini_batch_size = min(floor(Int64, scsg_pb.mini_batch_size_init * alfa_j), max_mini_batch_size)
    # m_j  computed to ensure mean number of mini batch < large_batch_size as mini_batch_size_init < large_batch_size_init is enfored
    nb_mini_batch_parameter = scsg_pb.m_zero * nbterms * alfa_j^1.5
    step_size = scsg_pb.eta_zero / sqrt(alfa_j)
    #
    BatchSizeInfo(iteration, large_batch_size, mini_batch_size, nb_mini_batch_parameter, step_size)
end

# select sizeasked terms selected without replacement
# returns the selected terms in a Vector{Int64}
function samplewithoutreplacement(size_asked::Int64, terms:: Vector{Int64}, rng::MersenneTwister)
    size_in = length(terms)
    out_terms = Vector{Int64}()
    for i in 1:length(terms)
        if rand(rng) * (size_in - i) < size_asked - length(out_terms)
            push!(out_terms, terms[i])
        end
    end
    @assert size_asked == length(out_terms)
    out_terms
end

# used to sample mini batches, Really faster.
# See:
#  1. Faster methods for Random Sampling J.S Vitter Comm ACM 1984
#  2. Kim-Hung Li Reservoir Sampling Algorithms : Comm ACM Vol 20, 4 December 1994
#  3. https://en.wikipedia.org/wiki/Reservoir_sampling
#
# For passing the rng as arg see also:
# https://bkamins.github.io/julialang/2020/11/20/rand.html or julia bloggers : "The need for rand speed"

function samplewithoutreplacement_reservoir(size_asked::Int64, in_terms::Vector{Int64}, rng::MersenneTwister)
    out_terms = in_terms[1:size_asked]
    w  = exp(log(rand(rng)/(size_asked)))
    s = size_asked
    while s <= length(in_terms)
        s = s + floor(Int64, log(rand(rng))/log1p(-w)) + 1
        if s <= length(in_terms)
            k = 1 + floor(Int64, size_asked * rand(rng))
            out_terms[k] = in_terms[s]
            w = w * exp(log(rand(rng))/size_asked)
        end
    end
    out_terms
end

"""
# function estimate_batch_growing_factor(batchinfo:: BatchSizeInfo, nbiterations::Int64, nbterms::Int64)

computes batch growing factor used to control large and mini batch sizes.
It cannot be too large, so it must be adjusted according to nbterms value, and nb_max_iterations
The growing factor alfa should be such that :
       B_0 * alfa^(2nbiterations) < nbterms

"""
function get_batchgrowingfactor(scsg_pb::SCSG, nbiterations::Int64, nbterms::Int64)
    if scsg_pb.m_zero > nbterms
        error("m_zero > nbterms in function to minimize, exiting")
    end
    log_alfa = -log(scsg_pb.large_batch_size_init) / (2. * nbiterations)
    batch_growing_factor = exp(log_alfa)
    if batch_growing_factor <= 1.
        @warn("batch growing factor shoud be greater than 1. , possibly you can reduce number of iterations ")
    end
    @debug "batch growing factor : " batch_growing_factor
    batch_growing_factor
end

"""
# function get_nbminibatch(batch_info::BatchSizeInfo)
"""
function get_nbminibatch(batch_info::BatchSizeInfo)
    m_j =  batch_info.nbminibatchparam
    b_j = batch_info.mini_batchsize
    # we return mean of geometric. Sampling too much instable due to large variance of geometric distribution.
    # ensure that it is at least 1 with ceil
    n_j = ceil(Int64, m_j/b_j)
    n_j = min(n_j, batch_info.large_batchsize)
    @debug " nb mini batch = " n_j
    n_j
end

"""
# function minimize(scsgpb::SCSG, evaluation::Evaluator{F,G}, max_iter::Int64, initialposition::Array{Float64,N})

Generic function used in minimisations for all algorithms.

The function has 3 types parameters F, G and N.

- F and G corresponds to the types of function used in instantiating the structure TermFunction{F} and TermmGradient{G}
(see in test and examples directories the cases: TestLinearRegression or TestLogisticRegression).

- N is the dimension of array for searched parameters (position and gradients) which depends upon the problem and the way we modelize it.
(See Logistic regression where we used N=2)

## Args
- scsgpb : the structure describing main parameters of the batch strategy
- Evaluator{F,G}
- max_iter : maximum number of iterations
- initialposition : initial position of the iterations

"""
function minimize(scsg_pb::SCSG, evaluation::Evaluator{F,G}, max_iterations, initial_position::Array{Float64,N}) where {F,G,N}
    #
    direction = zeros(Float64, size(initial_position))
    large_batch_gradient = zeros(Float64, size(initial_position))
    mini_batch_gradient_current = zeros(Float64, size(initial_position))
    mini_batch_gradient_origin = zeros(Float64, size(initial_position))
    nbterms = get_nbterms(evaluation)
    all_index = collect(1:nbterms)
    batch_growing_factor = get_batchgrowingfactor(scsg_pb, max_iterations, nbterms)
    #
    mt = MersenneTwister(117);
    position = Array{Float64}(initial_position)
    iteration = 0
    initial_value = 0
    more = true
    while more
        iteration += 1
        # get iteration parameters
        batch_info = get_batchsizeinfo(scsg_pb, batch_growing_factor, nbterms, iteration)
        @debug "batch_info" batch_info
        # batch sampling
        batch_indexes = samplewithoutreplacement(batch_info.large_batchsize, all_index, mt)
        # compute gradient on large batch index set and store initial position
        compute_gradient!(evaluation, position , batch_indexes, large_batch_gradient)
        # get mean number  binomial law for number Nj of small batch iterations
        nb_mini_batch = get_nbminibatch(batch_info)
        position_before_mini_batch = Array{Float64}(position)
        # loop on small batch iterations
        for i in 1:nb_mini_batch
            # sample mini batch terms
            batch_indexes = samplewithoutreplacement_reservoir(batch_info.mini_batchsize, all_index, mt)
            compute_gradient!(evaluation, position , batch_indexes, mini_batch_gradient_current)
            compute_gradient!(evaluation, position_before_mini_batch , batch_indexes, mini_batch_gradient_origin)
            direction = mini_batch_gradient_current - mini_batch_gradient_origin + large_batch_gradient;
            position = position - batch_info.stepsize * direction;
        end
        @debug "norm L2 direction" norm(direction)
        value = compute_value(evaluation, position)
        if iteration == 1
            initial_value = value
            @debug "initial value " initial_value
        end
        if iteration >= max_iterations
            @info "Reached maximal number of iterations required , initial_value value " initial_value value
            return position, value
        end
    end
end


text:
## --- Test sorting functions directly

    A = rand(100)

    # SortNaNs
    B, iₗ, iᵤ = NaNStatistics.sortnans!(copy(A))
    @test B == A
    @test (iₗ, iᵤ) == (1, 100)

    # Quicksort
    NaNStatistics.quicksort!(A)
    sort!(B)
    @test A == B

    A = rand(1_000)
    B = sort(A)
    NaNStatistics.quicksort!(A)
    @test A == B

    A = rand(1_000_000)
    B = sort(A)
    NaNStatistics.quicksort!(A)
    @test A == B

    # Multithreaded quicksort
    A = rand(100)
    B = sort(A)
    NaNStatistics.quicksortt!(A)
    @test A == B

    A = rand(1_000)
    B = sort(A)
    NaNStatistics.quicksortt!(A)
    @test A == B

    A = rand(1_000_000)
    B = sort(A)
    NaNStatistics.quicksortt!(A)
    @test A == B

    # Partialsort
    A = rand(101)
    m = median(A)
    NaNStatistics.quickselect!(A, 1, 101, 51)
    @test A[51] == m

    A = rand(1_000_001)
    m = median(A)
    NaNStatistics.quickselect!(A, 1, 1_000_001, 500_001)
    @test A[500_001] == m

    # Quicksort of already-sorted arrays
    @test NaNStatistics.quicksort!(collect(1:100)) == 1:100
    @test NaNStatistics.quicksort!(collect(100:-1:1)) == 1:100
    @test NaNStatistics.quicksortt!(collect(1:100)) == 1:100
    @test NaNStatistics.quicksortt!(collect(100:-1:1)) == 1:100

    # # Vsort, Float64
    # A = rand(100)
    # B = NaNStatistics.vsort(A, multithreaded=false)
    # @test issorted(B)
    # A = rand(100)
    # B = NaNStatistics.vsort(A, multithreaded=true)
    # @test issorted(B)
    #
    # # Vsort, Int64
    # A = rand(Int, 100)
    # B = NaNStatistics.vsort(A, multithreaded=false)
    # @test issorted(B)
    # A = rand(Int, 100)
    # B = NaNStatistics.vsort(A, multithreaded=true)
    # @test issorted(B)

## --- Test nanmedian!

    @test nanmedian!(0:10) == 5
    @test nanmedian!(1:10) == 5.5

    A = rand(100)
    @test nanmedian!(copy(A)) == median(A)
    A[rand(1:100, 20)] .= NaN; B = A[.!isnan.(A)];
    @test nanmedian!(A) == median(B)
    A = rand(10_000)
    @test nanmedian!(copy(A)) == median(A)
    A[rand(1:10_000, 100)] .= NaN; B = A[.!isnan.(A)]
    @test nanmedian!(A) == median(B)
    A = rand(100_000)
    @test nanmedian!(copy(A)) == median(A)
    A[rand(1:100_000, 100)] .= NaN; B = A[.!isnan.(A)]
    @test nanmedian!(A) == median(B)

    A = rand(55,82)
    @test nanmedian!(copy(A)) == median(A)
    @test nanmedian!(copy(A), dims=1) == median(A, dims=1)
    @test nanmedian!(copy(A), dims=2) == median(A, dims=2)

    A = rand(10,11,12)
    @test nanmedian!(copy(A)) == median(A)
    @test nanmedian!(copy(A), dims=1) == median(A, dims=1)
    @test nanmedian!(copy(A), dims=2) == median(A, dims=2)
    @test nanmedian!(copy(A), dims=3) == median(A, dims=3)
    @test nanmedian!(copy(A), dims=(1,2)) == median(A, dims=(1,2))
    @test nanmedian!(copy(A), dims=(2,3)) == median(A, dims=(2,3))

## --- Test nanpctile! / nanquantile!

    @test nanpctile!(0:10, 0) == 0
    @test nanpctile!(0:10, 1) ≈ 0.1
    @test nanpctile!(0:10, 100) == 10
    @test nanpctile!(0:10, 13.582) ≈ 1.3582
    @test nanpctile!(collect(1:10), 50) == 5.5

    A = rand(100)
    @test nanpctile!(copy(A), 50) == median(A)
    A[rand(1:100, 20)] .= NaN; B = A[.!isnan.(A)];
    @test nanpctile!(A, 50) == median(B)
    A = rand(10_000)
    @test nanpctile!(copy(A), 50) == median(A)
    A[rand(1:10_000, 100)] .= NaN; B = A[.!isnan.(A)];
    @test nanpctile!(A, 50) == median(B)
    A = rand(100_000)
    @test nanpctile!(copy(A), 50) == median(A)
    A[rand(1:100_000, 100)] .= NaN; B = A[.!isnan.(A)];
    @test nanpctile!(A, 50) == median(B)

    A = rand(55,82)
    @test nanpctile!(copy(A), 50) == median(A)
    @test nanpctile!(copy(A), 50, dims=1) == median(A, dims=1)
    @test nanpctile!(copy(A), 50, dims=2) == median(A, dims=2)

    A = rand(10,11,12)
    @test nanpctile!(copy(A), 50) == median(A)
    @test nanpctile!(copy(A), 50, dims=1) == median(A, dims=1)
    @test nanpctile!(copy(A), 50, dims=2) == median(A, dims=2)
    @test nanpctile!(copy(A), 50, dims=3) == median(A, dims=3)
    @test nanpctile!(copy(A), 50, dims=(1,2)) == median(A, dims=(1,2))
    @test nanpctile!(copy(A), 50, dims=(2,3)) == median(A, dims=(2,3))

## ---


text:
@testset "Testing functions" begin
    product_year = 2018
    tracer = "p"
    period = 0 # Annual, 1 month and 1 season — no need to test every month and season
    resolution = "5°"
    field = "mn"
    ds = WOA.WOA_Dataset(product_year, tracer, period, resolution)
    @test ds isa Dataset
    @testset "get_3D_field" begin
        field3D = WOA.get_3D_field(product_year, tracer, period, resolution, field)
        @test field3D isa Array{Union{Float32, Missing},3}
        @test size(field3D) == (36, 72, 102)
    end
    @testset "get_gridded_3D_field" begin
        field3D, lat, lon, depth = WOA.get_gridded_3D_field(ds, tracer, field)
        @test field3D isa Array{Union{Float32, Missing},3}
        @test size(field3D) == (36, 72, 102)
        @test lat isa Vector
        @test length(lat) == 36
        @test lon isa Vector
        @test length(lon) == 72
        @test depth isa Vector
        @test length(depth) == 102
    end
    @testset "filter_gridded_3D_field" begin
        field3D, lat, lon, depth = WOA.get_gridded_3D_field(ds, tracer, field)
        fieldvec, latvec, lonvec, depthvec, CI = WOA.filter_gridded_3D_field(field3D, lat, lon, depth)
        @test fieldvec isa Vector
        @test latvec isa Vector
        @test lonvec isa Vector
        @test depthvec isa Vector
    end
    @testset "fit_to_grid" begin
        grid = OceanGrid(90,180,24)
        a, b = WOA.fit_to_grid(grid, product_year, tracer, period, resolution, field)
        println(typeof(a))
        println(typeof(b))
        @test a isa Array{Float64, 3}
        @test size(a) == size(grid)
        @test b isa Array{Float64, 3}
        @test size(b) == size(grid)
        I = findall(a .≠ 0)
        @test 1e-3 < sum(a[I]) / length(I) < 3e-3 # mean PO₄ ≈ 2e-3 mol m⁻³
    end
end


text:
module FranklinParser

import REPL.REPLCompletions: emoji_symbols
import Base.isempty

const SS = SubString{String}
const SubVector{T} = SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int64}}, true}

subv(v::Vector{T}) where T = @view v[1:length(v)]
subv(v::SubVector{T}) where T = v

include("utils/strings.jl")
include("utils/types.jl")
include("utils/regex.jl")
include("utils/errors.jl")

include("tokens/utils.jl")
include("tokens/find_tokens.jl")
include("tokens/md_utils.jl")
include("tokens/_md_tokens.jl")
include("tokens/_html_tokens.jl")
include("tokens/_args_tokens.jl")

# see partition, we put this here because it's also used in 'form.jl'
const INLINE_BLOCKS = [
    :TEXT,
    :COMMENT,                                 # <!-- ... -->
    :RAW, :RAW_HTML, :RAW_LATEX,              # ???...???, ~~~...~~~, %%%...%%%
    :EMPH_EM, :EMPH_STRONG, :EMPH_EM_STRONG,  # * ** ***, _ __ ____
    :LINEBREAK,                               # \\
    #
    :CODE_INLINE,                             # `...`
    :MATH_INLINE,                             # $...$
    #
    :AUTOLINK,                                # <...>
    :LINK_A, :LINK_AB, :LINK_AR,              # ![...] ![...](...) ![...][...]
    :IMG_A, :IMG_AB, :IMG_AR,                 # [...] [...](...) [...][...]
    #
    :CU_BRACKETS, :LX_COMMAND,
    :LX_NEWENVIRONMENT, :LX_NEWCOMMAND,
    :DBB,
    # derived by reconstructing commands (Franklin)
    :RAW_INLINE
]

# these are blocks which, if present in a paragraph, necessarily make that paragraph
# a paragraph. Here's an example:
#
# [:TEXT][:RAW_HTML][:TEXT] --> if the [:TEXT] are non-empty, it's necessarily a paragraph
#
# [:RAW_HTML] --> not a paragraph because it's on it's own
#
# [:TEXT][:RAW_HTML] -> not a paragraph if the stripped text is empty
#
const INLINE_BLOCKS_CHECKP = [
    :COMMENT,
    :RAW, :RAW_HTML, :RAW_LATEX,
    :LINEBREAK,
    :CU_BRACKETS, :LX_COMMAND,
    :LX_NEWENVIRONMENT, :LX_NEWCOMMAND,
    :DBB
]

include("blocks/form.jl")
include("blocks/find_blocks.jl")
include("blocks/_md_blocks.jl")
include("blocks/_html_blocks.jl")
include("blocks/_args_blocks.jl")
include("blocks/utils.jl")

include("partition.jl")

end


text:
Ux = rand()
Uy = rand()
α₀ = rand()
α̇₀ = rand()
ax = rand()
ay = rand()
Ω = 2π
Ax = rand()
ϕx = rand()
Ay = rand()
ϕy = rand()
Δα = rand()
ϕα = rand()

@testset "Motions" begin

  oscil = RigidBodyTools.Oscillation(Ux,Uy,α̇₀,ax,ay,Ω,Ax,Ay,ϕx,ϕy,α₀,Δα,ϕα)

  t = rand()
  c, ċ, c̈, α,α̇,α̈ = oscil(t)

  @test α ≈ α₀ + α̇₀*t + Δα*sin(Ω*t-ϕα)
  @test α̇ ≈ α̇₀ + Δα*Ω*cos(Ω*t-ϕα)
  @test α̈ ≈ -Δα*Ω^2*sin(Ω*t-ϕα)

  a = complex(ax,ay)
  @test real(c) ≈ Ux*t + Ax*sin(Ω*t-ϕx) - ax*cos(α) + ay*sin(α)
  @test imag(c) ≈ Uy*t + Ay*sin(Ω*t-ϕy) - ay*cos(α) - ax*sin(α)
  @test real(ċ) ≈ Ux + Ax*Ω*cos(Ω*t-ϕx) + α̇*(ax*sin(α) + ay*cos(α))
  @test imag(ċ) ≈ Uy + Ay*Ω*cos(Ω*t-ϕy) + α̇*(ay*sin(α) - ax*cos(α))
  @test real(c̈) ≈ -Ax*Ω^2*sin(Ω*t-ϕx) + α̈*(ax*sin(α) + ay*cos(α)) + α̇^2*(ax*cos(α) - ay*sin(α))
  @test imag(c̈) ≈ -Ay*Ω^2*sin(Ω*t-ϕy) + α̈*(ay*sin(α) - ax*cos(α)) + α̇^2*(ay*cos(α) + ax*sin(α))

  ph = RigidBodyTools.PitchHeave(Ux,ax,Ω,α₀,Δα,ϕα,Ay,ϕy)

  t = rand()
  c, ċ, c̈, α,α̇,α̈ = ph(t)

  @test α ≈ α₀ + Δα*sin(Ω*t-ϕα)
  @test α̇ ≈ Δα*Ω*cos(Ω*t-ϕα)
  @test α̈ ≈ -Δα*Ω^2*sin(Ω*t-ϕα)

  @test real(c) ≈ Ux*t - ax*cos(α)
  @test imag(c) ≈ Ay*sin(Ω*t-ϕy) - ax*sin(α)
  @test real(ċ) ≈ Ux + α̇*ax*sin(α)
  @test imag(ċ) ≈ Ay*Ω*cos(Ω*t-ϕy) - α̇*ax*cos(α)
  @test real(c̈) ≈ α̈*ax*sin(α) + α̇^2*ax*cos(α)
  @test imag(c̈) ≈ -Ay*Ω^2*sin(Ω*t-ϕy) - α̈*ax*cos(α) + α̇^2*ax*sin(α)

  ox = RigidBodyTools.OscillationX(Ux,Ω,Ax,ϕx)

  t = rand()
  c, ċ, c̈, α,α̇,α̈ = ox(t)
  @test α ≈ 0.0
  @test α̇ ≈ 0.0
  @test α̈ ≈ 0.0
  @test real(c) ≈ Ux*t + Ax*sin(Ω*t-ϕx)
  @test imag(c) ≈ 0.0
  @test real(ċ) ≈ Ux + Ax*Ω*cos(Ω*t-ϕx)
  @test imag(ċ) ≈ 0.0
  @test real(c̈) ≈ -Ax*Ω^2*sin(Ω*t-ϕx)
  @test imag(c̈) ≈ 0.0

  oy = RigidBodyTools.OscillationY(Uy,Ω,Ay,ϕy)

  t = rand()
  c, ċ, c̈, α,α̇,α̈ = oy(t)
  @test α ≈ 0.0
  @test α̇ ≈ 0.0
  @test α̈ ≈ 0.0
  @test real(c) ≈ 0.0
  @test imag(c) ≈ Uy*t + Ay*sin(Ω*t-ϕy)
  @test real(ċ) ≈ 0.0
  @test imag(ċ) ≈ Uy + Ay*Ω*cos(Ω*t-ϕy)
  @test real(c̈) ≈ 0.0
  @test imag(c̈) ≈ -Ay*Ω^2*sin(Ω*t-ϕy)

  ro = RigidBodyTools.RotationalOscillation(Ω,Δα,ϕα)

  t = rand()
  c, ċ, c̈, α,α̇,α̈ = ro(t)

  @test α ≈ Δα*sin(Ω*t-ϕα)
  @test α̇ ≈ Δα*Ω*cos(Ω*t-ϕα)
  @test α̈ ≈ -Δα*Ω^2*sin(Ω*t-ϕα)
  @test c ≈ 0.0
  @test ċ ≈ 0.0
  @test c̈ ≈ 0.0

end

@testset "Direct motions" begin

  m = RigidBodyTools.BasicDirectMotion()

  ml = RigidBodyTools.DirectlySpecifiedMotionList([m])

  push!(ml,m)

end


text:
# https://turing.ml/dev/tutorials/07-poisson-regression/

# Import Turing, Distributions and DataFrames
using Turing, Distributions, DataFrames, Distributed

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

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

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

# The true Poisson parameter (mean) of toy data
theta_noalcohol_meds = 1    # no alcohol, took medicine
theta_alcohol_meds = 3      # alcohol, took medicine
theta_noalcohol_nomeds = 6  # no alcohol, no medicine
theta_alcohol_nomeds = 36   # alcohol, no medicine

# no of samples for each of the above cases
q = 100

# Generate data from different Poisson distributions
noalcohol_meds = Poisson(theta_noalcohol_meds)
alcohol_meds = Poisson(theta_alcohol_meds)
noalcohol_nomeds = Poisson(theta_noalcohol_nomeds)
alcohol_nomeds = Poisson(theta_alcohol_nomeds)

nsneeze_data = vcat(rand(noalcohol_meds, q), rand(alcohol_meds, q), rand(noalcohol_nomeds, q), rand(alcohol_nomeds, q))
alcohol_data = vcat(zeros(q), ones(q), zeros(q), ones(q))
meds_data = vcat(zeros(q), zeros(q), ones(q), ones(q))

df = DataFrame(nsneeze = nsneeze_data, alcohol_taken = alchol_data, nomeds_taken = meds_data, product_alcohol_meds = meds_data.*alcohol_data)
df[sample(1:nrow(df), 5, replace = false), :]

# Data Plotting
p1 = Plots.histogram(df[(df[:,:alcohol_taken] .== 0) .& (df[:,:nomeds_taken] .== 0), 1], title = "no_alcohol+meds")
p2 = Plots.histogram((df[(df[:,:alcohol_taken] .== 1) .& (df[:,:nomeds_taken] .== 0), 1]), title = "alcohol+meds")
p3 = Plots.histogram((df[(df[:,:alcohol_taken] .== 0) .& (df[:,:nomeds_taken] .== 1), 1]), title = "no_alcohol+no_meds")
p4 = Plots.histogram((df[(df[:,:alcohol_taken] .== 1) .& (df[:,:nomeds_taken] .== 1), 1]), title = "alcohol+no_meds")
plot(p1, p2, p3, p4, layout = (2, 2), legend = false)

# Convert the DataFrame object to matrices.
data = Matrix(df[:, [:alcohol_taken, :nomeds_taken, :product_alcohol_meds]])
data_labels = df[:, :nsneeze]
data

# Rescale our matrices to help Turing sampler
# in initializing the parameter estimates(normalizing).
data = (data .- mean(data, dims=1)) ./ std(data, dims=1)

@model poisson_regression(x, y, n, σ²) = begin
    b0 ~ Normal(0, σ²)
    b1 ~ Normal(0, σ²)
    b2 ~ Normal(0, σ²)
    b3 ~ Normal(0, σ²)
    for i = 1:n
        theta = b0 + b1*x[i, 1] + b2*x[i,2] + b3*x[i,3]
        y[i] ~ Poisson(exp(theta))
    end
end;

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

# Sample using NUTS.
num_chains = 4
m = poisson_regression(data, data_labels, n, 10)
chain = sample(m, NUTS(200, 0.65), MCMCThreads(), 2_500, num_chains; discard_adapt=false)

# check whether our chains have converged
gelmandiag(chain)

# Taking the first chain
c1 = chain[:,:,1]

# Calculating the exponentiated means
b0_exp = exp(mean(c1[:b0]))
b1_exp = exp(mean(c1[:b1]))
b2_exp = exp(mean(c1[:b2]))
b3_exp = exp(mean(c1[:b3]))

print("The exponent of the meaned values of the weights (or coefficients are): \n")
print("b0: ", b0_exp, " \n", "b1: ", b1_exp, " \n", "b2: ", b2_exp, " \n", "b3: ", b3_exp, " \n")
print("The posterior distributions obtained after sampling can be visualised as :\n")

plot(chain)

# removing warmup samples

# Note the standard deviation before removing the warmup samples
describe(chain)

# Removing the first 200 values of the chains.
chains_new = chain[201:2500,:,:]
describe(chains_new)

plot(chains_new)


text:
push!(LOAD_PATH,"../src/")
using Documenter
using IntervalUnions

makedocs(
        modules  = [IntervalUnions],
        format   = :html,
        sitename = "IntervalUnions.jl",
        doctest  = false,
        pages    = Any[
                        "Getting Started" => "index.md",
                        "Intervals" => "intervals.md",
                        "IntervalUnions" => "intervalunions.md"
                        ]
        )

deploydocs(
            deps   = nothing,
            make   = nothing,
            repo   = "github.com/juliastreamgraphs/IntervalUnions.jl.git",
            branch = "gh-pages",
            target = "build",
          )