text:
#include("move.jl")
#include("move_helpers.jl")
#include("validate_alone.jl")

# simulates a game starting at current board state and runs until any king is killed or until max_moves reached
# returns king_killed, move_number of kill (1 for enemy king killed; -1 for our king killed)
# returns 0, move_number if simuation exceeds 100 simulated moves
# WARNING: this function is recursive
function simulate_normal(board, move_number, currentMove)

    #println("entering simulate_normal()") # for debugging

    const max_moves = 20

    if move_number - currentMove > max_moves
        return 0, move_number
    end

    ourPieces = []
    enemyPieces = []
    ourKing = piece("", "", 0, 0, -1, false)
    enemyKing = piece("", "", 0, 0, -1, false)

    for y = 1:length(board)
        for x = 1:length(board[1])
            if board[y][x].side == move_number % 2
                if board[y][x].name == "king"
                    ourKing = board[y][x]
                end
                push!(ourPieces, board[y][x])
            elseif board[y][x].side != -1
                if board[y][x].name == "king"
                    enemyKing = board[y][x]
                end
                push!(enemyPieces, board[y][x])
            end
        end
    end

    if ourKing == piece("", "", 0, 0, -1, false)
        if move_number%2 == currentMove%2
            return -1, move_number-1
        else
            return 1, move_number-1
        end
    end

    # check if enemy king wants to be killed
    for apiece in ourPieces
        if inRange(board, move(move_number, "move", apiece.x, apiece.y, enemyKing.x, enemyKing.y, "", false, -1, -1, -1, -1))
            if move_number%2 == currentMove%2
                return 1, move_number
            else
                return -1, move_number
            end
        end
    end

    safeMoves = [] # where any pieces cannot move to; i.e. safe moves for our pieces
    for y = 1:length(board)
        for x = 1:length(board[1])
            if board[y][x].side != currentMove%2
                for apiece in enemyPieces
                    if !inRange(board, move(currentMove, "move", apiece.x, apiece.y, x, y, "", false, -1, -1, -1, -1))
                        push!(safeMoves, (x,y))
                    end
                end
            end
        end
    end

    possibleMoves = []

    for apiece in ourPieces
        for amove in safeMoves
            # moving apiece to move is valid
            if inRange(board, move(move_number+1, "move", apiece.x, apiece.y, amove[1], amove[2], "", false, -1, -1, -1, -1))
                push!(possibleMoves, (apiece, amove))
            end
        end
    end

    srand()
    sourcePiece, target = possibleMoves[rand(1:length(possibleMoves))]
    sourceX = sourcePiece.x
    sourceY = sourcePiece.y
    targetX = target[1]
    targetY = target[2]

    board, move_number = simulateMove(board, move_number, sourceX, sourceY, targetX, targetY)

    return simulate_normal(board, move_number + 1, currentMove)

end

function AI_normal(board, currentMove)

    #println("entering AI_normal()") # for debugging

    gameActive = true

    #= start timing =#
    tic()

    currentSide = currentMove % 2 #getting current side for convenience

    #AI check 0: check if won

    ourPieces = Array(piece, 0) #0 element piece array
    enemyPieces = Array(piece, 0)
    enemyKing = empty
    ourKing = empty

    # go through board
    for i in 1:length(board)
        for j in 1:length(board[1])

            #get all our pieces in board into ourPieces array
            if currentSide == board[i][j].side
                push!(ourPieces, board[i][j])
            elseif board[i][j].side != -1
                push!(enemyPieces, board[i][j])
            end

            #get enemy king
            if board[i][j].name == "king" && board[i][j].side != currentSide
                enemyKing = board[i][j]
            #get our king
            elseif board[i][j].name == "king" && board[i][j].side == currentSide
                ourKing = board[i][j]
            end

        end
    end

    #= for debugging =#
    println("ourPieces: ", ourPieces)
    println("enemyPieces: ", enemyPieces)
    println("ourKing: ", ourKing)
    println("enemyKing: ", enemyKing)

    killerPieces = [] # array of pieces that can kill our king

    # identify pieces that can kill our king
    for apiece in enemyPieces
        if inRange(board, move(currentMove, "move", apiece.x, apiece.y, ourKing.x, ourKing.y, "", false, -1, -1, -1, -1))
            ##println  ("found a killer: $apiece")
            push!(killerPieces, apiece)
        end
    end

    println("killerPieces: ", killerPieces) # for debugging

    sX, sY, tX, tY = checkForWin(board, ourPieces,enemyKing)

    if sX != 0 # we can win this game
        moveMade = move(currentMove, "move", sX, sY, tX, tY, "", false, -1, -1, -1, -1)
        gameActive = false
    end

    println("result of checkForWin: sX $sX, sY $sY, tX $tX, tY $tY") # for debugging

    inCheck = checkForCheck(board, enemyPieces, ourKing)

    println("inCheck: $inCheck") # for debugging

    if inCheck && gameActive

        if length(killerPieces) == 1

            println("1 killer piece") # for debugging

            for apiece in ourPieces
                if inRange(board, move(currentMove, "move", apiece.x, apiece.y, killerPieces[1].x, killerPieces[1].y, "", false, -1, -1, -1, -1))
                    moveMade = move(currentMove,  "move", apiece.x, apiece.y, killerPieces[1].x, killerPieces[1].y, "", false, -1, -1, -1, -1)
                    board, currentMove = simulateMove(board, currentMove, apiece.x, apiece.y, killerPieces[1].x, killerPieces[1].y)
                    ##println  ("move made on line 916")
                    gameActive = false
                    break
                end
            end

            println("gameActive: $gameActive") # for debugging

            if gameActive

                # check if we can block the move
                sX, sY, tX, tY = checkForBlock(board, currentMove)

                println("result of checkForBlock: sX $sX, sY $sY, tX $tX, tY $tY") # for debugging

                if sX != 0 # block the move
                    moveMade = move(currentMove,  "move", sX, sY, tX, tY, "", false, -1, -1, -1, -1)
                    board, currentMove = simulateMove(board, currentMove, sX, sY, tX, tY)
                    ##println  ("move made on line 923")
                    gameActive = false
                else
                    tX, tY = checkForRun(board, currentMove)
                    ##println  ("checkForRun: $tX, $tY")
                    if tX != 0         # make king run
                        moveMade = move(currentMove,  "move", ourKing.x, ourKing.y, tX, tY, "", false, -1, -1, -1, -1)
                        board, currentMove = simulateMove(board, currentMove, ourKing.x, ourKing.y, tX, tY)
                        ##println  ("move made on line 930")
                        gameActive = false
                    else #king cannot run
                        ##println  ("move made on line 934")
                        moveMade = move(currentMove,  "resign", -1, -1, -1, -1, "", false, -1, -1, -1, -1)
                        gameActive = false
                    end
                end

            end
        else

            println("more than 1 killerPiece") # for debugging

            # check if king can run

            tX, tY = checkForRun(board, currentMove)
            ##println  ("checkForRun: $tX, $tY")
            if tX != 0         # make king run
                moveMade = move(currentMove,  "move", ourKing.x, ourKing.y, tX, tY, "", false, -1, -1, -1, -1)
                board, currentMove = simulateMove(board, currentMove, 0,0, tX, tY)
                #println  ("move made on line 947")
                gameActive = false
            else #king cannot run
                #println  ("move made on line 951")
                moveMade = move(currentMove,  "resign", -1, -1, -1, -1, "", false, -1, -1, -1, -1)
                gameActive = false
            end
        end

    end

    if gameActive

        println("!inCheck && gameActive")

        loopCount = 1
        validMoves = [] # valid moves for move_piece
        move_piece = piece("", "", 0, 0, -1, false)
        while gameActive
            srand()
            move_piece = ourPieces[rand(1:length(ourPieces))]
            for y = 1:length(board)
                for x = 1:length(board[1])
                    if board[y][x].side != move_piece.side # can't move onto friendly piece
                        if inRange(board, move(currentMove+1, "move", move_piece.x, move_piece.y, x, y, "", false, -1, -1, -1, -1))
                            push!(validMoves, (x,y))
                        end
                    end
                end
            end
            loopCount += 1
            if length(validMoves) >= 1
                break
            end
            if loopCount > 100
                println("resigning because loopCount exceeded 100")
                moveMade = move(currentMove,  "resign", -1, -1, -1, -1,"", false, -1, -1, -1, -1)
                gameActive = false
            end
        end

    end
    #println("piece selected: ",move_piece)
    #println("valid moves: \n",validMoves)

    if gameActive

        simulateResults = []

        myRes = 0
        gameLength = 0

        for move in validMoves
            #global myRes
            #global gameLength
           # println("piece selected 2: ",move_piece)
            newBoard, move_number = simulateMove(board, currentMove, move_piece.x, move_piece.y, move[1], move[2])
           # println("piece selected 3: ",move_piece)
            const numOfTrials = 10

            wins = 0
            losses = 0
            draws = 0
            winGameLenTotal = 0
            lossGameLenTotal = 0
            drawGameLenTotal = 0
            trialResults = []
            for i = 1:numOfTrials
                res, gameLen = simulate_normal(board, move_number, currentMove)
                if res == 1
                    wins += 1
                    winGameLenTotal += gameLen
                elseif res == -1
                    losses += 1
                    lossGameLenTotal += gameLen
                else
                    draws += 1
                    drawGameLenTotal += gameLen
                end
            end
            if wins >= losses && wins >= draws
                myRes = 1
                gameLength = winGameLenTotal / wins
            elseif losses >= draws
                myRes = -1
                gameLength = lossGameLenTotal / losses
            else
                myRes = 0
                gameLength = drawGameLenTotal / draws
            end
            #myRes, gameLength = simulate_normal(board, move_number, currentMove)
            push!(simulateResults, (myRes, gameLength))
        end

        #=
        insert 50 trials of simulate for each valid move
        find best move to make
        place final results into simulateResults
        =#

        winningMoves = []
        losingMoves = []
        miscMoves = []
        for i = 1:length(simulateResults)
            if simulateResults[i][1] == 1
                #pushing gameLength, sourcex, sourcey, targetx,targety
                push!(winningMoves, (simulateResults[i][2], move_piece.x, move_piece.y, validMoves[i][1], validMoves[i][2]))
            elseif simulateResults[i][1] == -1
                push!(losingMoves, (simulateResults[i][2], move_piece.x, move_piece.y, validMoves[i][1], validMoves[i][2]))
            else
                push!(miscMoves, (simulateResults[i][2], move_piece.x, move_piece.y, validMoves[i][1], validMoves[i][2]))
            end
        end

        if length(winningMoves) > 0
            # make the possibly winning move
            bestMove = winningMoves[1]
            for i = 2:length(winningMoves)
                if winningMoves[i][1] < bestMove[1] # smaller game length for winning
                    bestMove = winningMoves[i]
                end
            end
            sourceX = bestMove[2]
            sourceY = bestMove[3]
            targetX = bestMove[4]
            targetY = bestMove[5]
            #println("bestMove is: ", bestMove)
            deservesPromotion = promotion_check(board, sourceX, sourceY, targetX, targetY)
            option = deservesPromotion == true? "!":""

            moveMade = move(currentMove,  "move", sourceX, sourceY, targetX, targetY, option, false, -1, -1, -1, -1)
            board, currentMove = simulateMove(board, currentMove, sourceX, sourceY, targetX, targetY)
            #moveMade = true

        elseif length(losingMoves) > 0

            bestMove = losingMoves[1]
            for i = 2:length(losingMoves)
                if losingMoves[i][1] < bestMove[1] # smaller game length for winning
                    bestMove = losingMoves[i]
                end
            end
            sourceX = bestMove[2]
            sourceY = bestMove[3]
            targetX = bestMove[4]
            targetY = bestMove[5]
            #println("bestMove is: ", bestMove)
            deservesPromotion = promotion_check(board, sourceX, sourceY, targetX, targetY)
            option = deservesPromotion == true? "!":""

            moveMade = move(currentMove,  "move", sourceX, sourceY, targetX, targetY, option, false, -1, -1, -1, -1)
            board, currentMove = simulateMove(board, currentMove, sourceX, sourceY, targetX, targetY)
            #moveMade = true

        elseif length(miscMoves) > 0

            bestMove = miscMoves[1]
            for i = 2:length(miscMoves)
                if miscMoves[i][1] < bestMove[1] # smaller game length for winning
                    bestMove = miscMoves[i]
                end
            end
            sourceX = bestMove[2]
            sourceY = bestMove[3]
            targetX = bestMove[4]
            targetY = bestMove[5]
            #println("bestMove is: ", bestMove)
            deservesPromotion = promotion_check(board, sourceX, sourceY, targetX, targetY)
            option = deservesPromotion == true? "!":""

            moveMade = move(currentMove,  "move", sourceX, sourceY, targetX, targetY, option, false, -1, -1, -1, -1)
            board, currentMove = simulateMove(board, currentMove, sourceX, sourceY, targetX, targetY)
            #moveMade = true

        else # were screwed, so resign with dignity!

            moveMade = move(currentMove,  "resign", -1, -1, -1, -1, "", false, -1, -1, -1, -1)
            #moveMade = true

        end

    end # if moveMade != true

    #= end timing =#
    time_taken = toq()

    if length(board) == 16
        board = demon_burning(board)
    end

    return board, time_taken, moveMade
end

max_stars_repo_name: DarienFlamont/Shogi

max_issues_repo_name: DarienFlamont/Shogi

max_forks_repo_name: DarienFlamont/Shogi


text:
# NOTE: Simple test BCs to start with
function apply_BCs!(cache)
    u = cache.u
    u[1,:]   .= 0.0
    u[end,:] .= 0.0
    u[:,1]   .= 0.0
    u[:,end] .= 0.0
    return nothing
end

ext: jl


text:
using Wavelets

wf = waveletfilter("db2")
wl = waveletls("db2")
N = 1024*1024;
x0 = rand(N);
L = nscales(N)  # int(log2(N)-2)
tn = 10


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra
using Quaternions
using BenchmarkTools

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.open(vis)

@show Threads.nthreads()

# include("continuous_policy.jl")
include("continuous_policy_v2.jl")

# ## Simulation
s = get_simulation("point_foot_quadruped", "flat_3D_lc", "flat")
model = s.model
env = s.env

# ref_traj.z
# ref_traj.θ
# ref_traj.q
# ref_traj.u
# ref_traj.H
# κ = 2e-4
# N = 200
# z = [[deepcopy(ref_traj.z) for i = 1:N]...;]
# θ = [[deepcopy(ref_traj.θ) for i = 1:N]...;]
# rz = [LinearizedStep(s, z[t], θ[t], κ).rz for t = 1:N*H]
# rθ = [LinearizedStep(s, z[t], θ[t], κ).rθ for t = 1:N*H]
# idx_z = indices_z(s)
# idx_θ = indices_θ(model, nf=1)
# ∂z∂θ = [-rz[t] \ rθ[t] for t = 1:N*H]
# ∂q3∂q1 = [∂z∂θ[t][idx_z.q, idx_θ.q1] for t = 1:N*H]
# ∂q3∂q2 = [∂z∂θ[t][idx_z.q, idx_θ.q2] for t = 1:N*H]
# ∂q2∂q1 = [zeros(model.nq, model.nq) for t = 1:N*H]
# ∂q2∂q2 = [I(model.nq) for t = 1:N*H]
# ∂q2∂u = [zeros(model.nq, model.nu) for t = 1:N*H]
# ∂q3∂u = [∂z∂θ[t][idx_z.q, idx_θ.u] for t = 1:N*H]
#
# A = [[∂q2∂q1[t] ∂q3∂q1[t]; ∂q2∂q2[t] ∂q3∂q2[t]] for t = 1:N*H]
# B = [[∂q2∂u[t]; ∂q3∂u[t]] for t = 1:N*H]
# Q = [0.01 * I(2model.nq) for t = 1:N*H+1]
# R = [1000 * I(model.nu) for t = 1:N*H]
#
# A = [[1 1; 0 sin(2π*t/1000)] for t = 1:N*H]
# B = [[0 1.]' for t = 1:N*H]
# Q = [1 * I(2) for t = 1:N*H+1]
# R = [1 * I(1) for t = 1:N*H]
#
# K, P = tvlqr(A, B, Q, R)

# plot(hcat([reshape(K[t], 2) for t = 1:N*H]...)')
# function tvlqr(A, B, Q, R)
#     T = length(Q)
#
#     P = [zero(A[1]) for t = 1:T]
#     K = [zero(B[1]') for t = 1:T-1]
#     P[T] = Q[T]
#
#     for t = T-1:-1:1
#         K[t] = (R[t] + B[t]' * P[t+1] *  B[t]) \ (B[t]' * P[t+1] * A[t])
#         P[t] = (Q[t] + K[t]' * R[t] * K[t]
#                 + (A[t] - B[t] * K[t])' * P[t+1] * (A[t] - B[t] * K[t]))
#     end
#
#     return K, P
# end

# ## Reference Trajectory
ref_traj = deepcopy(get_trajectory(s.model, s.env,
# joinpath(module_dir(), "src/dynamics/point_foot_quadruped/gaits/inplace_trot_v0.jld2"),
	joinpath(module_dir(), "src/dynamics/point_foot_quadruped/gaits/inplace_trot_v1.jld2"),
    load_type = :split_traj_alt));

H = ref_traj.H
h = ref_traj.h

# ## MPC setup
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 4000
κ_mpc = 2.0e-4

v0 = 0.2
obj = TrackingVelocityObjective(model, env, H_mpc,
    v = [Diagonal(1e-3 * [[1,1,1]; 1e+3*[1,1,1]; fill([1,1,1], 4)...]) for t = 1:H_mpc],
	q = [relative_state_cost(1e-0*[1e-2,1e-2,1], 3e-1*[1,1,1], 1e-0*[0.2,0.2,1]) for t = 1:H_mpc],
	u = [Diagonal(3e-3 * vcat(fill([1,1,1], 4)...)) for t = 1:H_mpc],
	v_target = [1/ref_traj.h * [v0;0;0; 0;0;0; v0;0;0; v0;0;0; v0;0;0; v0;0;0] for t = 1:H_mpc],)

p = ci_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
	mode = :configuration,
	ip_opts = InteriorPointOptions(
					undercut = 5.0,
					κ_tol = κ_mpc,
					r_tol = 1.0e-4, # TODO maybe relax this parameter
					diff_sol = true,
					solver = :empty_solver,
					max_time = 1e5),
    n_opts = NewtonOptions(
        r_tol = 3e-5,
        max_time=10.0e-1,
		solver=:ldl_solver,
        threads=false,
        verbose=false,
        max_iter = 5),
    mpc_opts = CIMPCOptions(
		# live_plotting=true
		));

# ## Disturbances
w = [[0.0,0.0,0.0] for i=1:H_sim/N_sample]
d = open_loop_disturbances(w, N_sample)

# ## Initial conditions
q1_sim, v1_sim = initial_conditions(ref_traj);

# ## Simulator
sim = simulator(s, H_sim, h=h_sim, policy=p, dist=d);

# ## Simulate
q1_sim0 = deepcopy(q1_sim)
RoboDojo.simulate!(sim, q1_sim0, v1_sim)

# ## Visualize
set_light!(vis)
set_floor!(vis, grid=true)
set_background!(vis)
anim = visualize!(vis, model, sim.traj.q; Δt=h_sim)
# anim = visualize!(vis, model, ref_traj.q; Δt=h)

# # ## Timing result
# # Julia is [JIT-ed](https://en.wikipedia.org/wiki/Just-in-time_compilation) so re-run the MPC setup through Simulate for correct timing results.
process!(sim.stats, N_sample) # Time budget
H_sim * h_sim / sum(sim.stats.policy_time) # Speed ratio
plot(sim.stats.policy_time, xlabel="timestep", ylabel="mpc time (s)",
	ylims=[-0.001, 0.1],
	label="", linetype=:steppost)

sim
sim.traj

plt = plot()
plot!(plt, hcat(Vector.([(sim.traj.q[i+1][1:1] - sim.traj.q[i][1:1]) / sim.h for i=1:H_sim])...)')
plot!(plt, hcat(Vector.([(sim.traj.q[i+1][7:7] - sim.traj.q[i][7:7]) / sim.h for i=1:H_sim])...)')
plot!(plt, hcat(Vector.([(sim.traj.q[i+1][10:10] - sim.traj.q[i][10:10]) / sim.h for i=1:H_sim])...)')
plot!(plt, hcat(Vector.([(sim.traj.q[i+1][13:13] - sim.traj.q[i][13:13]) / sim.h for i=1:H_sim])...)')
plot!(plt, hcat(Vector.([(sim.traj.q[i+1][16:16] - sim.traj.q[i][16:16]) / sim.h for i=1:H_sim])...)')


text:
include("Include.jl")

# load the data dictionary -
data_dictionary = maximize_product_dictionary(0,10,1)

# solve the lp problem -
(objective_value, flux_array, dual_array, uptake_array, exit_flag) = FluxDriver(data_dictionary)
(calculated_flux_array, dual_value_array) = calculate_flux_variabilty(data_dictionary)

# FBA
println("objective: $objective_value")
println("flux array: ", flux_array)
println("uptake array: ", uptake_array)

# FVA
println("FVA results: ")
for i = 1:size(calculated_flux_array)[1]
  println("    ", calculated_flux_array[i,:])
end


text:
lines = readlines(ARGS[1])

function solution1(lines)
    val_matrix = []
    for line in lines
        line_vals = [parse(Int, v) for v in line]
        append!(val_matrix, [line_vals])
    end
    num_rows = length(val_matrix)
    row_width = length(val_matrix[1])
    minima = []
    for (i, row) in enumerate(val_matrix)
        for (j, val) in enumerate(row)
            is_minimal = true
            # check above
            if i > 1
                if val >= val_matrix[i-1][j]
                    is_minimal = false
                end
            end
            # check below
            if i < num_rows
                if val >= val_matrix[i+1][j]
                    is_minimal = false
                end
            end
            # check left
            if j > 1
                if val >= val_matrix[i][j-1]
                    is_minimal = false
                end
            end
            # check right
            if j < row_width
                if val >= val_matrix[i][j+1]
                    is_minimal = false
                end
            end
            if is_minimal
                append!(minima, val + 1)
            end
        end
    end
    return sum(minima)
end

function crawl_basin(basin, val_matrix, num_rows, row_width)
    to_search = Set(basin)
    found = Set(basin)
    while length(to_search) > 0
        check_point = pop!(to_search)
        check_x = check_point[1]
        check_y = check_point[2]
        # check above
        if check_y - 1 >= 1 && val_matrix[check_x][check_y-1] != 9 && !in([check_x, check_y - 1], found)
            push!(to_search, [check_x, check_y - 1])
            push!(found, [check_x, check_y - 1])
        end
        # check below
        if check_y + 1 <= num_rows && val_matrix[check_x][check_y+1] != 9 && !in([check_x, check_y + 1], found)
            push!(to_search, [check_x, check_y + 1])
            push!(found, [check_x, check_y + 1])
        end
        # check left
        if check_x - 1 >= 1 && val_matrix[check_x-1][check_y] != 9 && !in([check_x - 1, check_y], found)
            push!(to_search, [check_x - 1, check_y])
            push!(found, [check_x - 1, check_y])
        end
        # check right
        if check_x + 1 <= row_width && val_matrix[check_x+1][check_y] != 9 && !in([check_x + 1, check_y], found)
            push!(to_search, [check_x + 1, check_y])
            push!(found, [check_x + 1, check_y])
        end
    end
    return [x for x in found]
end

function solution2(lines)
    val_matrix = []
    for line in lines
        line_vals = [parse(Int, v) for v in line]
        append!(val_matrix, [line_vals])
    end
    num_rows = length(val_matrix)
    row_width = length(val_matrix[1])
    minima = []
    basins = []
    for (i, row) in enumerate(val_matrix)
        for (j, val) in enumerate(row)
            is_minimal = true
            basin = [[i, j]]
            # check above
            if i > 1
                if val >= val_matrix[i-1][j]
                    is_minimal = false
                else
                    if val_matrix[i-1][j] != 9
                        append!(basin, [[i - 1, j]])
                    end

                end
            end
            # check below
            if i < num_rows
                if val >= val_matrix[i+1][j]
                    is_minimal = false
                else
                    if val_matrix[i+1][j] != 9
                        append!(basin, [[i + 1, j]])
                    end
                end
            end
            # check left
            if j > 1
                if val >= val_matrix[i][j-1]
                    is_minimal = false
                else
                    if val_matrix[i][j-1] != 9
                        append!(basin, [[i, j - 1]])
                    end
                end
            end
            # check right
            if j < row_width
                if val >= val_matrix[i][j+1]
                    is_minimal = false
                else
                    if val_matrix[i][j+1] != 9
                        append!(basin, [[i, j + 1]])
                    end
                end
            end
            if is_minimal
                append!(minima, val + 1)
                append!(basins, [basin])
            end
        end
    end

    crawled_basins = []
    for basin in basins
        append!(crawled_basins, length(crawl_basin(basin, val_matrix, num_rows, row_width)))
    end

    rev_lengths = reverse(sort(crawled_basins))

    return rev_lengths[1] * rev_lengths[2] * rev_lengths[3]
end

println(solution1(lines))

println(solution2(lines))


text:
@testset "innerjoin" begin
    @testset "UnitRange" begin
	    l = 2:5
	    r = 3:6
	    @test innerjoin(identity, identity, tuple, l, r) == [(3,3), (4,4), (5,5)]
	    @test innerjoin(identity, identity, tuple, isequal, l, r) == [(3,3), (4,4), (5,5)]
	    @test innerjoin(identity, identity, tuple, ==, l, r) == [(3,3), (4,4), (5,5)]

	    @test innerjoin(identity, identity, tuple, isless, l, r) == [(2,3), (2,4), (2,5), (2,6), (3,4), (3,5), (3,6), (4,5), (4,6), (5,6)]
	end

	@testset "Arrays of NamedTuple" begin
	    l = [(a=1, b=2.0), (a=2, b=4.0), (a=3, b=6.0)]
	    r = [(a=1, c=:a), (a=2, c=:b), (a=4,c=:d)]
	    @test innerjoin(x->x.a, x->x.a, l, r) == [(a=1, b=2.0, c=:a), (a=2, b=4.0, c=:b)]
	end
end


text:
## Evaluation

function Evaluation(S::Jacobi,x::Union{typeof(leftendpoint),typeof(leftendpoint)},order)
    if order ≤ 2
        ConcreteEvaluation(S,x,order)
    else
        # assume Derivative is available
        D = Derivative(S,order)
        EvaluationWrapper(S,x,order,Evaluation(rangespace(D),x)*D)
    end
end

function Evaluation(S::Jacobi,x,order)
    if order == 0
        ConcreteEvaluation(S,x,order)
    else
        # assume Derivative is available
        D = Derivative(S,order)
        EvaluationWrapper(S,x,order,Evaluation(rangespace(D),x)*D)
    end
end

# avoid ambiguity
for OP in (:leftendpoint,:rightendpoint)
    @eval getindex(op::ConcreteEvaluation{<:Jacobi,typeof($OP)},k::Integer) =
        op[k:k][1]
end

getindex(op::ConcreteEvaluation{<:Jacobi},k::Integer) = op[k:k][1]

function getindex(op::ConcreteEvaluation{<:Jacobi,typeof(leftendpoint)},kr::AbstractRange)
    @assert op.order <= 2
    sp=op.space
    T=eltype(op)
    RT=real(T)
    a=convert(RT,sp.a);b=convert(RT,sp.b)

    if op.order == 0
        jacobip(T,kr.-1,a,b,-one(T))
    elseif op.order == 1&&  b==0
        d=domain(op)
        @assert isa(d,IntervalOrSegment)
        T[tocanonicalD(d,leftendpoint(d))/2*(a+k)*(k-1)*(-1)^k for k=kr]
    elseif op.order == 1
        d=domain(op)
        @assert isa(d,IntervalOrSegment)
        if kr[1]==1 && kr[end] ≥ 2
            tocanonicalD(d,leftendpoint(d)).*(a .+ b .+ kr).*T[zero(T);jacobip(T,0:kr[end]-2,1+a,1+b,-one(T))]/2
        elseif kr[1]==1  # kr[end] ≤ 1
            zeros(T,length(kr))
        else
            tocanonicalD(d,leftendpoint(d)) .* (a.+b.+kr).*jacobip(T,kr.-1,1+a,1+b,-one(T))/2
        end
    elseif op.order == 2
        @assert b==0
        @assert domain(op) == ChebyshevInterval()
        T[-0.125*(a+k)*(a+k+1)*(k-2)*(k-1)*(-1)^k for k=kr]
    else
        error("Not implemented")
    end
end

function getindex(op::ConcreteEvaluation{<:Jacobi,typeof(rightendpoint)},kr::AbstractRange)
    @assert op.order <= 2
    sp=op.space
    T=eltype(op)
    RT=real(T)
    a=convert(RT,sp.a);b=convert(RT,sp.b)

    if op.order == 0
        jacobip(T,kr.-1,a,b,one(T))
    elseif op.order == 1
        d=domain(op)
        @assert isa(d,IntervalOrSegment)
        if kr[1]==1 && kr[end] ≥ 2
            tocanonicalD(d,leftendpoint(d))*((a+b).+kr).*T[zero(T);jacobip(T,0:kr[end]-2,1+a,1+b,one(T))]/2
        elseif kr[1]==1  # kr[end] ≤ 1
            zeros(T,length(kr))
        else
            tocanonicalD(d,leftendpoint(d))*((a+b).+kr).*jacobip(T,kr.-1,1+a,1+b,one(T))/2
        end
    else
        error("Not implemented")
    end
end

function getindex(op::ConcreteEvaluation{<:Jacobi,<:Number},kr::AbstractRange)
    @assert op.order == 0
    jacobip(eltype(op),kr-1,op.space.a,op.space.b,tocanonical(domain(op),op.x))
end

## Derivative

Derivative(J::Jacobi,k::Integer)=k==1 ? ConcreteDerivative(J,1) : DerivativeWrapper(TimesOperator(Derivative(Jacobi(J.b+1,J.a+1,J.domain),k-1),ConcreteDerivative(J,1)),k)

rangespace(D::ConcreteDerivative{J}) where {J<:Jacobi}=Jacobi(D.space.b+D.order,D.space.a+D.order,domain(D))
bandwidths(D::ConcreteDerivative{J}) where {J<:Jacobi}=0,D.order

getindex(T::ConcreteDerivative{J},k::Integer,j::Integer) where {J<:Jacobi} =
    j==k+1 ? eltype(T)((k+1+T.space.a+T.space.b)/complexlength(domain(T))) : zero(eltype(T))

## Integral

function Integral(J::Jacobi,k::Integer)
    if k > 1
        Q=Integral(J,1)
        IntegralWrapper(TimesOperator(Integral(rangespace(Q),k-1),Q),k)
    elseif J.a > 0 && J.b > 0   # we have a simple definition
        ConcreteIntegral(J,1)
    else   # convert and then integrate
        sp=Jacobi(J.b+1,J.a+1,domain(J))
        C=Conversion(J,sp)
        Q=Integral(sp,1)
        IntegralWrapper(TimesOperator(Q,C),1)
    end
end

rangespace(D::ConcreteIntegral{J}) where {J<:Jacobi}=Jacobi(D.space.b-D.order,D.space.a-D.order,domain(D))
bandwidths(D::ConcreteIntegral{J}) where {J<:Jacobi}=D.order,0

function getindex(T::ConcreteIntegral{J},k::Integer,j::Integer) where J<:Jacobi
    @assert T.order==1
    if k≥2 && j==k-1
        complexlength(domain(T))./(k+T.space.a+T.space.b-2)
    else
        zero(eltype(T))
    end
end

## Volterra Integral operator

Volterra(d::IntervalOrSegment) = Volterra(Legendre(d))
function Volterra(S::Jacobi, order::Integer)
    @assert S.a == S.b == 0.0
    @assert order==1
    ConcreteVolterra(S,order)
end

rangespace(V::ConcreteVolterra{J}) where {J<:Jacobi}=Jacobi(-1.0,0.0,domain(V))
bandwidths(V::ConcreteVolterra{J}) where {J<:Jacobi}=1,0

function getindex(V::ConcreteVolterra{J},k::Integer,j::Integer) where J<:Jacobi
    d=domain(V)
    C = complexlength(d)/2
    if k≥2
        if j==k-1
            C/(k-1.5)
        elseif j==k
            -C/(k-0.5)
        else
            zero(eltype(V))
        end
    else
        zero(eltype(V))
    end
end

for (Func,Len,Sum) in ((:DefiniteIntegral,:complexlength,:sum),(:DefiniteLineIntegral,:arclength,:linesum))
    ConcFunc = Meta.parse("Concrete"*string(Func))

    @eval begin
        $Func(S::Jacobi{<:IntervalOrSegment}) = $ConcFunc(S)

        function getindex(Σ::$ConcFunc{Jacobi{D,R},T},k::Integer) where {D<:IntervalOrSegment,R,T}
            dsp = domainspace(Σ)

            if dsp.b == dsp.a == 0
                # TODO: copy and paste
                k == 1 ? convert(T,$Sum(Fun(dsp,[one(T)]))) : zero(T)
            else
                convert(T,$Sum(Fun(dsp,[zeros(T,k-1);1])))
            end
        end

        function bandwidths(Σ::$ConcFunc{Jacobi{D,R}}) where {D<:IntervalOrSegment,R}
            if domainspace(Σ).b == domainspace(Σ).a == 0
                0,0  # first entry
            else
                0,∞
            end
        end
    end
end

## Conversion
# We can only increment by a or b by one, so the following
# multiplies conversion operators to handle otherwise

function Conversion(L::Jacobi,M::Jacobi)
    domain(L) == reverseorientation(domain(M)) &&
        return ConversionWrapper(Conversion(reverseorientation(L), M)*ReverseOrientation(L))

    domain(L) == domain(M) || domain(L) == reverseorientation(domain(M)) ||
        throw(ArgumentError("Domains must be the same"))

    if isapproxinteger(L.a-M.a) && isapproxinteger(L.b-M.b)
        dm=domain(M)
        D=typeof(dm)
        if isapprox(M.a,L.a) && isapprox(M.b,L.b)
            ConversionWrapper(Operator(I,L))
        elseif (isapprox(M.b,L.b+1) && isapprox(M.a,L.a)) || (isapprox(M.b,L.b) && isapprox(M.a,L.a+1))
            ConcreteConversion(L,M)
        elseif M.b > L.b+1
            ConversionWrapper(TimesOperator(Conversion(Jacobi(M.b-1,M.a,dm),M),Conversion(L,Jacobi(M.b-1,M.a,dm))))
        else  #if M.a >= L.a+1
            ConversionWrapper(TimesOperator(Conversion(Jacobi(M.b,M.a-1,dm),M),Conversion(L,Jacobi(M.b,M.a-1,dm))))
        end
    elseif L.a ≈ L.b ≈ 0. && M.a ≈ M.b ≈ 0.5
        Conversion(L,Ultraspherical(L),Ultraspherical(M),M)
    elseif L.a ≈ L.b ≈ 0. && M.a ≈ M.b ≈ -0.5
        Conversion(L,Ultraspherical(L),Chebyshev(M),M)
    elseif L.a ≈ L.b ≈ -0.5 && M.a ≈ M.b ≈ 0.5
        Conversion(L,Chebyshev(L),Ultraspherical(M),M)
    else # L.a - M.a ≈ L.b - M.b
        error("Implement for $L → $M")
    end
end

bandwidths(C::ConcreteConversion{J1,J2}) where {J1<:Jacobi,J2<:Jacobi}=(0,1)

function Base.getindex(C::ConcreteConversion{J1,J2,T},k::Integer,j::Integer) where {J1<:Jacobi,J2<:Jacobi,T}
    L=C.domainspace
    if L.b+1==C.rangespace.b
        if j==k
            k==1 ? convert(T,1) : convert(T,(L.a+L.b+k)/(L.a+L.b+2k-1))
        elseif j==k+1
            convert(T,(L.a+k)./(L.a+L.b+2k+1))
        else
            zero(T)
        end
    elseif L.a+1==C.rangespace.a
        if j==k
            k==1 ? convert(T,1) : convert(T,(L.a+L.b+k)/(L.a+L.b+2k-1))
        elseif j==k+1
            convert(T,-(L.b+k)./(L.a+L.b+2k+1))
        else
            zero(T)
        end
    else
        error("Not implemented")
    end
end

# return the space that has banded Conversion to the other
function conversion_rule(A::Jacobi,B::Jacobi)
    if !isapproxinteger(A.a-B.a) || !isapproxinteger(A.b-B.b)
        NoSpace()
    else
        Jacobi(min(A.b,B.b),min(A.a,B.a),domain(A))
    end
end

## Ultraspherical Conversion

# Assume m is compatible

function Conversion(A::PolynomialSpace,B::Jacobi)
    J = Jacobi(A)
    J == B ? ConcreteConversion(A,B) :
             ConversionWrapper(TimesOperator(Conversion(J,B),Conversion(A,J)))
end

function Conversion(A::Jacobi,B::PolynomialSpace)
    J = Jacobi(B)
    J == A ? ConcreteConversion(A,B) :
             ConversionWrapper(TimesOperator(Conversion(J,B),Conversion(A,J)))
end

function Conversion(A::Jacobi,B::Chebyshev)
    if A.a == A.b == -0.5
        ConcreteConversion(A,B)
    elseif A.a == A.b == 0
        ConversionWrapper(
            SpaceOperator(
                Conversion(Ultraspherical(1//2),B),
                A,B))
    elseif A.a == A.b
        US = Ultraspherical(A)
        ConversionWrapper(Conversion(US,B)*Conversion(A,US))
    else
        J = Jacobi(B)
        Conversion(J,B)*Conversion(A,J)
    end
end

function Conversion(A::Chebyshev,B::Jacobi)
    if B.a == B.b == -0.5
        ConcreteConversion(A,B)
    elseif B.a == B.b == 0
        ConversionWrapper(
            SpaceOperator(
                Conversion(A,Ultraspherical(1//2,domain(B))),
                A,B))
    elseif B.a == B.b
        US = Ultraspherical(B)
        ConversionWrapper(Conversion(US,B)*Conversion(A,US))
    else
        J = Jacobi(A)
        Conversion(J,B)*Conversion(A,J)
    end
end

function Conversion(A::Jacobi,B::Ultraspherical)
    if A.a == A.b == -0.5
        ConversionWrapper(Conversion(Chebyshev(domain(A)),B)*Conversion(A,Chebyshev(domain(A))))
    elseif A.a == A.b == order(B)-0.5
        ConcreteConversion(A,B)
    elseif A.a == A.b == 0
        ConversionWrapper(
            SpaceOperator(
                Conversion(Ultraspherical(1//2),B),
                A,B))
    elseif A.a == A.b
        US = Ultraspherical(A)
        ConversionWrapper(Conversion(US,B)*Conversion(A,US))
    else
        J = Jacobi(B)
        Conversion(J,B)*Conversion(A,J)
    end
end

function Conversion(A::Ultraspherical,B::Jacobi)
    if B.a == B.b == -0.5
        ConversionWrapper(Conversion(Chebyshev(domain(A)),B)*Conversion(A,Chebyshev(domain(A))))
    elseif B.a == B.b == order(A)-0.5
        ConcreteConversion(A,B)
    elseif B.a == B.b == 0
        ConversionWrapper(
            SpaceOperator(
                Conversion(A,Ultraspherical(1//2,domain(B))),
                A,B))
    elseif B.a == B.b
        US = Ultraspherical(B)
        ConversionWrapper(Conversion(US,B)*Conversion(A,US))
    else
        J = Jacobi(A)
        Conversion(J,B)*Conversion(A,J)
    end
end

bandwidths(C::ConcreteConversion{US,J}) where {US<:Chebyshev,J<:Jacobi} = 0,0
bandwidths(C::ConcreteConversion{J,US}) where {US<:Chebyshev,J<:Jacobi} = 0,0

bandwidths(C::ConcreteConversion{US,J}) where {US<:Ultraspherical,J<:Jacobi} = 0,0
bandwidths(C::ConcreteConversion{J,US}) where {US<:Ultraspherical,J<:Jacobi} = 0,0

#TODO: Figure out how to unify these definitions
function getindex(C::ConcreteConversion{CC,J,T},k::Integer,j::Integer) where {J<:Jacobi,CC<:Chebyshev,T}
    if j==k
        one(T)/jacobip(T,k-1,-one(T)/2,-one(T)/2,one(T))
    else
        zero(T)
    end
end

function BandedMatrix(S::SubOperator{T,ConcreteConversion{CC,J,T},Tuple{UnitRange{Int},UnitRange{Int}}}) where {J<:Jacobi,CC<:Chebyshev,T}
    ret=BandedMatrix(Zeros, S)
    kr,jr = parentindices(S)
    k=(kr ∩ jr)

    vals = one(T)./jacobip(T,k .- 1,-one(T)/2,-one(T)/2,one(T))

    ret[band(bandshift(S))] = vals
    ret
end

function getindex(C::ConcreteConversion{J,CC,T},k::Integer,j::Integer) where {J<:Jacobi,CC<:Chebyshev,T}
    if j==k
        jacobip(T,k-1,-one(T)/2,-one(T)/2,one(T))
    else
        zero(T)
    end
end

function BandedMatrix(S::SubOperator{T,ConcreteConversion{J,CC,T},Tuple{UnitRange{Int},UnitRange{Int}}}) where {J<:Jacobi,CC<:Chebyshev,T}
    ret=BandedMatrix(Zeros, S)
    kr,jr = parentindices(S)
    k=(kr ∩ jr)

    vals = jacobip(T,k.-1,-one(T)/2,-one(T)/2,one(T))

    ret[band(bandshift(S))] = vals
    ret
end

function getindex(C::ConcreteConversion{US,J,T},k::Integer,j::Integer) where {US<:Ultraspherical,J<:Jacobi,T}
    if j==k
        S=rangespace(C)
        jp=jacobip(T,k-1,S.a,S.b,one(T))
        um=convert(Operator{T}, Evaluation(setcanonicaldomain(domainspace(C)),rightendpoint,0))[k]::T
        (um/jp)::T
    else
        zero(T)
    end
end

function BandedMatrix(S::SubOperator{T,ConcreteConversion{US,J,T},Tuple{UnitRange{Int},UnitRange{Int}}}) where {US<:Ultraspherical,J<:Jacobi,T}
    ret=BandedMatrix(Zeros, S)
    kr,jr = parentindices(S)
    k=(kr ∩ jr)

    sp=rangespace(parent(S))
    jp=jacobip(T,k.-1,sp.a,sp.b,one(T))
    um=Evaluation(T,setcanonicaldomain(domainspace(parent(S))),rightendpoint,0)[k]
    vals = um./jp

    ret[band(bandshift(S))] = vals
    ret
end

function getindex(C::ConcreteConversion{J,US,T},k::Integer,j::Integer) where {US<:Ultraspherical,J<:Jacobi,T}
    if j==k
        S=domainspace(C)
        jp=jacobip(T,k-1,S.a,S.b,one(T))
        um=Evaluation(T,setcanonicaldomain(rangespace(C)),rightendpoint,0)[k]
        jp/um::T
    else
        zero(T)
    end
end

function BandedMatrix(S::SubOperator{T,ConcreteConversion{J,US,T},Tuple{UnitRange{Int},UnitRange{Int}}}) where {US<:Ultraspherical,J<:Jacobi,T}
    ret=BandedMatrix(Zeros, S)
    kr,jr = parentindices(S)
    k=(kr ∩ jr)

    sp=domainspace(parent(S))
    jp=jacobip(T,k.-1,sp.a,sp.b,one(T))
    um=Evaluation(T,setcanonicaldomain(rangespace(parent(S))),rightendpoint,0)[k]
    vals = jp./um

    ret[band(bandshift(S))] = vals
    ret
end

function union_rule(A::Jacobi,B::Jacobi)
    if domainscompatible(A,B)
        Jacobi(min(A.b,B.b),min(A.a,B.a),domain(A))
    else
        NoSpace()
    end
end

function maxspace_rule(A::Jacobi,B::Jacobi)
    if isapproxinteger(A.a-B.a) && isapproxinteger(A.b-B.b)
        Jacobi(max(A.b,B.b),max(A.a,B.a),domain(A))
    else
        NoSpace()
    end
end

function union_rule(A::Chebyshev,B::Jacobi)
    if isapprox(B.a,-0.5) && isapprox(B.b,-0.5)
        # the spaces are the same
        A
    else
        union(Jacobi(A),B)
    end
end
function union_rule(A::Ultraspherical,B::Jacobi)
    m=order(A)
    if isapprox(B.a,m-0.5) && isapprox(B.b,m-0.5)
        # the spaces are the same
        A
    else
        union(Jacobi(A),B)
    end
end

for (OPrule,OP) in ((:conversion_rule,:conversion_type), (:maxspace_rule,:maxspace))
    @eval begin
        function $OPrule(A::Chebyshev,B::Jacobi)
            if B.a ≈ -0.5 && B.b ≈ -0.5
                # the spaces are the same
                A
            elseif isapproxinteger(B.a+0.5) && isapproxinteger(B.b+0.5)
                $OP(Jacobi(A),B)
            else
                NoSpace()
            end
        end
        function $OPrule(A::Ultraspherical,B::Jacobi)
            m = order(A)
            if B.a ≈ m-0.5 && B.b ≈ m-0.5
                # the spaces are the same
                A
            elseif isapproxinteger(B.a+0.5) && isapproxinteger(B.b+0.5)
                $OP(Jacobi(A),B)
            else
                NoSpace()
            end
        end
    end
end

hasconversion(a::Jacobi,b::Chebyshev) = hasconversion(a,Jacobi(b))
hasconversion(a::Chebyshev,b::Jacobi) = hasconversion(Jacobi(a),b)

hasconversion(a::Jacobi,b::Ultraspherical) = hasconversion(a,Jacobi(b))
hasconversion(a::Ultraspherical,b::Jacobi) = hasconversion(Jacobi(a),b)

## Special Multiplication
# special multiplication operators exist when multiplying by
# (1+x) or (1-x) by _decreasing_ the parameter.  Thus the

# represents [b+(1+z)*d/dz] (false) or [a-(1-z)*d/dz] (true)
struct JacobiSD{T} <:Operator{T}
    lr::Bool
    S::Jacobi
end

JacobiSD(lr,S)=JacobiSD{Float64}(lr,S)

convert(::Type{Operator{T}},SD::JacobiSD) where {T}=JacobiSD{T}(SD.lr,SD.S)

domain(op::JacobiSD)=domain(op.S)
domainspace(op::JacobiSD)=op.S
rangespace(op::JacobiSD)=op.lr ? Jacobi(op.S.b+1,op.S.a-1,domain(op.S)) : Jacobi(op.S.b-1,op.S.a+1,domain(op.S))
bandwidths(::JacobiSD)=0,0

function getindex(op::JacobiSD,A,k::Integer,j::Integer)
    m=op.lr ? op.S.a : op.S.b
    if k==j
        k+m-1
    else
        zero(eltype(op))
    end
end


text:
__precompile__(true)

module RationalExtensions

import Base.convert,
       Base.promote_rule,
       Base.show,
       Base.norm,
       Base.conj,
       Base.*,
       Base.//,
       Base./,
       Base.+,
       Base.-,
       Base.==

import Primes

export Rad,
       Sqrt,
       AbstractRational

AbstractRational = Union{Rational,Integer}

immutable Rad{S<:AbstractRational,T<:Integer} <: Real
    a::S
    b::S
    n::T

    function Rad{S,T}(a::S,b::S,n::T) where {S,T}
        if b == 0 || n == 0
            return new(a,convert(S,0),convert(T,1))
        end
        sgn = sign(n)
        d = Primes.factor(abs(n))
        for p in keys(d)
            b *= p^div(d[p],2)
            d[p] = mod(d[p],2)
        end
        n = isempty(d) ? 1 : prod(p^(d[p]) for p in keys(d))
        if sgn * n != 1
            return new(a,b,sgn * n)
        else
            return new(a+b,0,1)
        end
    end
end

Rad{S<:AbstractRational,T<:Integer}(a::S,b::S,n::T) = Rad{S,T}(a,b,n)
Rad{S<:AbstractRational,U<:AbstractRational,T<:Integer}(a::S,b::U,n::T) = Rad(promote(a,b)...,n)
Rad{S<:AbstractRational,T<:Integer}(a::S,b::S,r::Rational{T}) =
    Rad(a,b//denominator(r),numerator(r)*denominator(r))

Sqrt{T<:Integer}(n::T) = Rad{T,T}(zero(n),one(n),n)
Sqrt(r::Rational) = Rad(0,1,r)
Sqrt(x::Rad) = (x.n == 1 ? Sqrt(b(x)) : throw(InexactError()))
Sqrt(x::AbstractFloat) = sqrt(x)

==(x::Rad{S,T},y::Rad{S,T}) where {S,T} = (x.a,x.b,x.n) == (y.a,y.b,y.n)

function *(x::Rad,y::Rad)
    if x.n == y.n
        return Rad(x.a * y.a + x.b * y.b * x.n, y.a * x.b + y.b * x.a, x.n)
    elseif x.a == y.a == 0
        return Rad(zero(x.a), x.b * y.b, x.n * y.n)
    else
        error("Error multiplying radicals: radicands must mutiply to give a perfect square")
    end
end

*(x::Bool,y::Rad) = Rad( x * y.a, x * y.b, y.n)
*{T<:AbstractRational}(x::T,y::Rad) = Rad(x * y.a, x * y.b, y.n)
*(x::Rad,y::Bool) = y * x
*{T<:AbstractRational}(x::Rad,y::T) = y * x

function //(x::Rad,y::Rad)
    if x.n == y.n || x.b == 0 || y.b == 0
        return (x.a * y.a + y.a * x.b * Sqrt(x.n) - x.a * y.b * Sqrt(y.n)
                - x.b * y.b * Sqrt(x.n * y.n))//(y.a^2 - y.n * y.b^2)
    else
        error("Division with different radicands not supported")
    end
end

//{T<:AbstractRational}(x::T,y::Rad) = Rad(x,zero(x),one(x)) // y
//{T<:AbstractRational}(x::Rad,y::T) = Rad(x.a//y, x.b//y,x.n)

/(x::Rad,y::Rad) = float(x)/float(y)
/(x::Rad,y::AbstractRational) = convert(typeof(float(y)),x)/float(y)
/{T<:Real}(x::Rad,y::T) = /(promote(x,y)...)
/{T<:Real}(x::T,y::Rad) = /(promote(x,y)...)

function +(x::Rad,y::Rad)
    if x.n == y.n || x.b == 0
        return Rad(x.a + y.a, x.b+y.b, y.n)
    elseif y.b == 0
        return Rad(x.a + y.a, x.b, x.n)
    else
        error("Addition with different radicands not supported")
    end
end

-(x::Rad) = Rad(-x.a,-x.b,x.n)
-(x::Rad,y::Rad) = x + (-y)
+{T<:AbstractRational}(x::Rad,y::T) = Rad(x.a + y, x.b, x.n)
+{T<:AbstractRational}(x::T,y::Rad) = y + x
-{T<:AbstractRational}(x::Rad,y::T) = x + (-y)
-{T<:AbstractRational}(x::T,y::Rad) = x + (-y)

convert{S<:AbstractRational,T<:Integer}(::Type{Rad{S,T}},n::Integer) = Rad(n,zero(n),one(n))
convert{S<:AbstractRational,T<:Integer}(::Type{Rad{S,T}},r::Rational) = Rad(r,zero(r),one(r))
convert(::Type{AbstractFloat}, x::Rad) = float(x.a) + float(x.b) * (sqrt∘float)(x.n)

convert{S,T}(::Type{BigFloat},x::Rad{S,T}) = big(x.a) + big(x.b) * (sqrt∘float)(big(x.n))

convert{S<:AbstractFloat,T,U}(::Type{S},x::Rad{T,U}) =
    ((a,b,n) -> a + b*(sqrt∘float)(n))(promote(x.a,x.b,x.n)...)

convert(::Type{Bool},x::Rad) = (x.b != 0)
convert{T<:Integer}(::Type{T}, x::Rad) = (x.n == 1 ? convert(T,x.a + x.b) : throw(InexactError()))
convert{S<:AbstractRational,T<:Integer}(::Type{Rad{S,T}},x::Rad) =
     Rad(convert(S,x.a),convert(S,x.b),convert(T,x.n))

promote_rule{S<:AbstractRational,T<:Integer,U<:Integer}(::Type{Rad{S,T}}, ::Type{U}) =
     Rad{S,promote_type(T,U)}
promote_rule{S<:AbstractRational,T<:Integer,U<:Integer}(::Type{Rad{S,T}}, ::Type{Rational{U}}) =
     Rad{promote_type(S,U),T}
promote_rule{S<:AbstractRational,T<:Integer,U<:AbstractRational,V<:Integer}(::Type{Rad{S,T}},
     ::Type{Rad{U,V}}) = Rad{promote_type(S,U),promote_type(T,V)}
promote_rule{T<:Rad}(::Type{BigFloat},::Type{T}) = BigFloat
promote_rule{S<:AbstractFloat,T<:Rad}(::Type{S},::Type{T}) = S

convert(::Type{Rad}, x::Rad) = x
convert(::Type{Rad}, r::AbstractRational) = Rad(r,zero(r),one(r))

conj(x::Rad) = Rad(x.a,-x.b,x.n)
norm(x::Rad) = x * conj(x)

function show(io::IO,x::Rad)
    if x == 0
        show(io,0)
    elseif isa(x.b,Rational)
        if x.a != 0
            if x.a != 1
                show(io,x.a)
            else
                show(io,numerator(x.a))
            end
            if sign(x.b) == 1
                print(io," + ")
            elseif sign(x.b) == -1
                print(io," - ")
            end
        end
        if ~in(numerator(x.b),[-1,0,1])
            show(io,abs(numerator(x.b)))
        end
        if x.n != 1
            print(io,"√(")
            show(io,x.n)
            print(io,")")
        end
        if denominator(x.b) != 1
            print(io,"//")
            show(io,denominator(x.b))
        end
    else
        if x.a != 0
            show(io,x.a)
            if sign(x.b) == 1
                print(io," + ")
            elseif sign(x.b) == -1
                print(io," - ")
            end
        elseif sign(x.b) == -1
            print(io,"-")
        end
        if x.b != 0
            if abs(x.b) != 1
                show(io,abs(x.b))
            end
            if x.n != 1
                print(io,"√(")
                show(io,x.n)
                print(io,")")
            end
            if x.n == 1 && x.b == 1
                show(io,1)
            end
        end
    end
end

end # module


text:
function Kprec( pw::PWGrid, psi::Array{ComplexF64,2} )

    Ngwx  = size(psi)[1]
    Ncols = size(psi)[2]
    G2    = pw.gvec.G2[pw.gvecw.idx_gw2r]
    Kpsi  = zeros( ComplexF64, Ngwx, Ncols )

    for ic = 1:Ncols
        for ig = 1:Ngwx
            Kpsi[ig,ic] = psi[ig,ic] / ( 1.0 + G2[ig] )
        end
    end
    return Kpsi
end

function Kprec(psi)
    return psi
end


text:
using Test
using QueryOperators

@testset "NamedTupleUtilities" begin

@test QueryOperators.NamedTupleUtilities.select((a = 1, b = 2, c = 3), Val(:a)) == (a = 1,)
@test QueryOperators.NamedTupleUtilities.select((a = 1, b = 2, c = 3), Val(:d)) == NamedTuple()
@inferred QueryOperators.NamedTupleUtilities.remove((a = 1, b = 2), Val(:b))

@test QueryOperators.NamedTupleUtilities.remove((a = 1, b = 2, c = 3),Val(:c)) == (a = 1, b = 2)
@test QueryOperators.NamedTupleUtilities.remove((a = 1, b = 2),Val(:c)) == (a = 1, b = 2)
@inferred QueryOperators.NamedTupleUtilities.remove((a = 1, b = 2),Val(:c))

@test QueryOperators.NamedTupleUtilities.range((a = 1, b = 2, c = 3),Val(:a),Val(:b)) == (a = 1, b = 2)
@test QueryOperators.NamedTupleUtilities.range((a = 1, b = 2, c = 3),Val(:b),Val(:d)) == (b = 2, c = 3)
@test QueryOperators.NamedTupleUtilities.range((a = 1, b = 2, c = 3),Val(:d),Val(:c)) == NamedTuple()
@inferred QueryOperators.NamedTupleUtilities.range((a = 1, b = 2, c = 3),Val(:a),Val(:b))

@test QueryOperators.NamedTupleUtilities.rename((a = 1, b = 2, c = 3),Val(:a),Val(:d)) == (d = 1, b = 2, c = 3)
@test QueryOperators.NamedTupleUtilities.rename((a = 1, b = 2, c = 3),Val(:m),Val(:d)) == (a = 1, b = 2, c = 3)
@test_throws ErrorException QueryOperators.NamedTupleUtilities.rename((a = 1, b = 2, c = 3),Val(:a),Val(:c))
@inferred QueryOperators.NamedTupleUtilities.rename((a = 1, b = 2, c = 3),Val(:a),Val(:d))

@test QueryOperators.NamedTupleUtilities.startswith((abc=1,bcd=2,cde=3),Val(:a)) == (abc = 1,)
@test QueryOperators.NamedTupleUtilities.not_startswith((abc=1,bcd=2,cde=3),Val(:a)) == (bcd = 2, cde = 3)
@test QueryOperators.NamedTupleUtilities.endswith((abc=1,bcd=2,cde=3),Val(:d)) == (bcd = 2,)
@test QueryOperators.NamedTupleUtilities.not_endswith((abc=1,bcd=2,cde=3),Val(:d)) == (abc = 1, cde = 3)
@test QueryOperators.NamedTupleUtilities.occursin((abc=1,bcd=2,cde=3),Val(:d)) == (bcd = 2, cde = 3)
@test QueryOperators.NamedTupleUtilities.not_occursin((abc=1,bcd=2,cde=3),Val(:d)) == (abc = 1,)
@inferred QueryOperators.NamedTupleUtilities.startswith((abc=1,bcd=2,cde=3),Val(:a))
@inferred QueryOperators.NamedTupleUtilities.not_startswith((abc=1,bcd=2,cde=3),Val(:a))
@inferred QueryOperators.NamedTupleUtilities.endswith((abc=1,bcd=2,cde=3),Val(:d))
@inferred QueryOperators.NamedTupleUtilities.not_endswith((abc=1,bcd=2,cde=3),Val(:d))
@inferred QueryOperators.NamedTupleUtilities.occursin((abc=1,bcd=2,cde=3),Val(:d))
@inferred QueryOperators.NamedTupleUtilities.not_occursin((abc=1,bcd=2,cde=3),Val(:d))

nt = (a=4,b=true,c="Named")
@test QueryOperators.NamedTupleUtilities.oftype(nt, Val(Int)) == (a=4,)
@test QueryOperators.NamedTupleUtilities.oftype(nt, Val(Any)) == nt
@test QueryOperators.NamedTupleUtilities.oftype(nt, Val(Float64)) == NamedTuple()
@inferred QueryOperators.NamedTupleUtilities.oftype(nt, Val(Int))

end


text:
using Pkg
import Pkg; Pkg.add("Test")
using Test, PandaMPC

using ECOS
using Convex
using Plots

@testset "SimulationB" begin
    include("SimulationBtest.jl")
end


text:
"""
```
usual_model_settings!(m, vint; cdvt = vint, dsid = data_id(m), cdid = cond_id(m),
    fcast_date = Dates.lastdayofquarter(Dates.today()),
    altpolicy = AltPolicy(:historical, eqcond, solve))
```

Apply usual defaults for the following settings:

- `data_vintage` and `cond_vintage`: given by input argument `vint`
- `date_forecast_start` and `date_conditional_end`: given by kwarg `fcast_date`
- `use_population_forecast`: `true`
- `alternative_policy`: given by input argument `altpolicy`. If this argument is
  specified, then `altpolicy_settings!` and `altpolicy.setup` are also called.
"""
function usual_model_settings!(m::AbstractDSGEModel, vint::String;
                               cdvt::String = vint,
                               dsid::Int = data_id(m),
                               cdid::Int = cond_id(m),
                               fcast_date::Dates.Date = Dates.lastdayofquarter(Dates.today()),
                               altpolicy::AltPolicy = AltPolicy(:historical, eqcond, solve))
    m <= Setting(:data_vintage, vint)
    m <= Setting(:cond_vintage, cdvt)
    m <= Setting(:data_id, dsid)
    m <= Setting(:cond_id, cdid)
    m <= Setting(:date_forecast_start,  fcast_date)
    m <= Setting(:date_conditional_end, fcast_date)
    m <= Setting(:use_population_forecast, true)

    if altpolicy.key != :historical
        altpolicy_settings!(m, altpolicy)
    end
end

"""
```
usual_model_forecast(m, input_type, cond_type,
    output_vars = [:histobs, :histpseudo, :forecastobs, :forecastpseudo];
    est_override = "", forecast_string = "",
    density_bands = [0.5, 0.6, 0.7, 0.8, 0.9],
    mb_matrix = false, check_empty_columns = true)
```

Forecast, compute means and bands, and optionally (if `mb_matrix`) convert
`MeansBands` to matrices. If the path `est_override` is provided, it will be
added to `forecast_input_file_overrides(m)`.
"""
function usual_model_forecast(m::AbstractDSGEModel, input_type::Symbol, cond_type::Symbol,
                              output_vars::Vector{Symbol} = [:histobs, :histpseudo, :forecastobs, :forecastpseudo,
                                                             :shockdecobs, :shockdecpseudo];
                              est_override::String = "",
                              forecast_string::String = "",
                              density_bands::Vector{Float64} = [0.5, 0.6, 0.7, 0.8, 0.9],
                              mb_matrix::Bool = false,
                              shock_name::Symbol = :none,
                              shock_var_name::Symbol = :none,
                              shock_var_value::Float64 = 0.0,
                              check_empty_columns = true)

    # Override estimation file if necessary
    if !isempty(est_override)
        overrides = forecast_input_file_overrides(m)
        overrides[input_type] = est_override
    end

    # Forecast
    forecast_one(m, input_type, cond_type, output_vars; forecast_string = forecast_string,
                 verbose = :high,
                 shock_name = shock_name,
                 shock_var_name = shock_var_name,
                 shock_var_value = shock_var_value,
                 check_empty_columns = check_empty_columns)

    # Compute means and bands
    compute_meansbands(m, input_type, cond_type, output_vars; forecast_string = forecast_string,
                       density_bands = density_bands, check_empty_columns = check_empty_columns,
                       verbose = :high)

    # Convert MeansBands to matrices if desired
    if mb_matrix
        meansbands_to_matrix(m, input_type, cond_type, output_vars; forecast_string = forecast_string,
                             verbose = :high)
    end
end


text:
include(joinpath(@__DIR__, "branch_variable.jl"))
include(joinpath(@__DIR__, "branch_special.jl"))

_branchvar_lower_solution(m,i) = m._lower_solution[m._branch_to_sol_map[i]]
function _fractional_integer_branch!(t::ExtensionType, m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    if !m._current_node.cont
        integer_count = 0
        for i in m.branch_variables
            if _is_unfixed_integer(m, i)
                isinteger(_lower_solution(BranchVar, m, i))
                integer_count += 1
            end
        end
    end
    return false
end

#=
5.3.1 of Belotti, Pietro, et al. "Branching and bounds tighteningtechniques for
non-convex MINLP." Optimization Methods & Software 24.4-5 (2009): 597-634.
=#
function _select_branch_point!(m::GlobalOptimizer{N,T,S}, i) where {N,T<:Real,S}
    l = _lower_bound(BranchVar, m, i)
    u = _upper_bound(BranchVar, m, i)
    s = _lower_solution(BranchVar, m, i)
    α = _branch_cvx_α(m)
    b = _branch_offset_β(m)*(u - l)
    return max(l + b, min(u - b, α*s + (one(T) - α)*_mid(BranchVar, m, i)))
end

function _continuous_branch!(t::ExtensionType, m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}

    bi = _select_branch_variable!(m)
    bp = _select_branch_point!(m)

    # rounds into branch points, which in turn prevents the
    # solution at the branch point from being discarded
    lower_bound = max(n.lower_bound, m._lower_objective_value)
    upper_bound = min(n.upper_bound, m._upper_objective_value)
    new_depth = n.depth + 1

    m._maximum_node_id += 1
    push!(m._stack, NodeBB(_upper_bound(:branch, m),
                           ntuple(i -> (i == bi ? nextfloat(bp) : _upper_bound(BranchVar, m, i)), N),
                           lower_bound, upper_bound, new_depth, m._maximum_node_id,
                           bi, BD_NEG))

    m._maximum_node_id += 1
    push!(m._stack, NodeBB(ntuple(i -> (i == bi ? prevfloat(bp) : _lower_bound(BranchVar, m, i)), N),
                           _upper_bound(:branch, m),
                           lower_bound, upper_bound, new_depth, m._maximum_node_id,
                           bi, BD_POS))

    m._node_repetitions = 1
    m._node_count += 2
    return true
end

"""
$(SIGNATURES)

Creates two nodes from `current_node` using information available the `x`
and stores them to the stack. By default, relative width bisection is perfomed
at a point `branch_pnt` which is a convex combination
(parameter: `branch_cvx_factor`) of the solution to the relaxation and
the midpoint of the node. If this solution lies within `branch_offset/width` of
a bound then the branch point is moved to a distance of `branch_offset/width`
from the bound.
"""
function branch_node!(t::ExtensionType, m::GlobalOptimizer)
    was_branched = _conic_branch!(t, m)
    !was_branched && (was_branched &= _fractional_integer_branch!(t, m))
    !was_branched && (was_branched &= _continuous_branch!(t, m))
    return nothing
end


text:
function update_beta!(state::State, data::Data, prior::Prior,
                      flags::Vector{Flag})
  if UpdateBetaWithSkewT() in flags
    llC = marginal_loglike_beta('C', state, data)
    llT = marginal_loglike_beta('T', state, data)
  else
    llC = marginal_loglike_beta_latent_var('C', state, data)
    llT = marginal_loglike_beta_latent_var('T', state, data)
  end

  p1 = logistic(logit(state.p) + llT - llC)
  state.beta = p1 > rand()
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
"""
        function find_initial_MAP(model,trans;serowaningrate = 0.0)

This finds a MAP estimator for the `model` object with variable transformation `trans`. Step 1 is BlackBox optimisation
        using a differential evolution optimizer, which is then refined in step 2 by a local Nelder-Mead optimisation.
"""
function find_initial_MAP(model::CoVAreaModel,trans;serowaningrate = 0.0,VE_acquisition=zeros(3),VE_infectiousness=zeros(3),VE_severe_disease_risk=zeros(3))
	println("Searching for a good initial condition for $(model.areaname)")
	l_area(x) = model.log_likelihood(x,model,serowaningrate,VE_acquisition,VE_infectiousness,VE_severe_disease_risk) + model.log_priors(x)
	f(x) = -transform_logdensity(trans, l_area,x)
	searchrange = fill((-3.,3.),TransformVariables.dimension(trans))
	res = bboptimize(f; SearchRange = searchrange,PopulationSize=100,MaxSteps=10000)
	q₀ = best_candidate(res)
	res_nm = optimize(f,q₀,NelderMead(),Optim.Options(iterations = 1000,show_trace = true))
	return res_nm.minimizer
end

"""
        function PSGLD!(model,samples,trans,stepsize_params,q₀;serowaningrate = 0.0,beta = 0.99,λ =1e-5)

This draws a chain of parameter samples of length `samples` from the log-posterior defined by the `model` object using pre-conditioned Stochastic
Gradient Langevin Dynamics (PSGLD). The step size at time step `t` is `a*(b + t)^(-γ)`. The `a`, `b` and `γ` values must be given
in the `NamedTuple` object `stepsize_params`. The starting point is at `q₀`. The method is based on _C. Li, et al, Preconditioned stochastic gradient Langevin dynamics for deep neural networks. Proceedings of the AAAI Conference on Artificial Intelligence. 30, 1788–1794 (2016)_
As in the paper, the variable `beta` is the learning rate for the pre-conditioner and `λ` sets the maximum implied curvature of the parameter space. These are set to defaults from the paper.

The MCMC results are stored in the model object as per usual.
"""
function PSGLD!(model::CoVAreaModel,samples,trans,stepsize_params,q₀;serowaningrate = 0.0,VE_acquisition=zeros(3),VE_infectiousness=zeros(3),VE_severe_disease_risk=zeros(3),beta = 0.99,λ =1e-5)
	θ = copy(q₀)
	l_area(x) = model.log_likelihood(x,model,serowaningrate,VE_acquisition,VE_infectiousness,VE_severe_disease_risk) + model.log_priors(x)
	f(x) = transform_logdensity(trans, l_area,x)
	#ℓ = TransformedLogDensity(trans, l_area)#transformed log-likelihood
	#∇ℓ = LogDensityProblems.ADgradient(:ForwardDiff, ℓ)#transformed log-likelihood gradient wrt the parameters

	log_post_den_vals = zeros(samples)
	param_draws = zeros(samples/10,length(θ)) # save every 10th draw
	precond = zeros(length(θ))
	G_diag = zeros(length(θ))

	for t in 1:samples
		################GRADIENT CALCULATIONS
		#lpd,∇L = LogDensityProblems.logdensity_and_gradient(∇ℓ,θ)
		lpd = f(θ)
		∇L =  ForwardDiff.gradient(f,θ)
		#Step size
		ϵ = stepsize_params.a*(stepsize_params.b + t)^(-stepsize_params.γ)

		###############PRECONDITIONING#####################
		#Learns the curvature of the space
		precond *= beta
		precond += (1-beta)*(∇L .*∇L)
		G_diag .= 1.0./(λ .+ sqrt.(precond))

		###############LOG-POSTERIOR ASCENT WITH RANDOM WALK COMPONENT#############
		θ .+= 0.5.*ϵ.*G_diag.*∇L  .+ sqrt.(ϵ.*G_diag).*randn(length(∇L))

		#################BOOKKEEPING############################
		if t % 10 == 0  # edited to save every 10the draw
			param_draws[t,:] = θ
			log_post_den_vals[t] = lpd
			println(θ[1])
			println("On step $(t) log-posterior density is $(lpd). Step size is $(ϵ*1e4)e-4. Max component of G is $(maximum(G_diag)).")
		end
	end
	transformed_samples = [TransformVariables.transform(trans,param_draws[t,:]) for t in 1:size(param_draws,1)]
	val = zeros(length(transformed_samples),length(transformed_samples[1]),1)
	for i = 1:size(val,1),j = 1:size(val,2)
		val[i,j,1] = transformed_samples[i][j]
	end
	chn = Chains(val,[String(k) for k in keys(transformed_samples[1])])
	model.MCMC_results = KenyaCoVaccines.MCMCResults(chn,
                                        [l_area(transformed_samples[i]) - model.log_priors(transformed_samples[i]) for i = 1:length(transformed_samples)],
                                        DynamicHMC.TreeStatisticsNUTS[])

	return param_draws,log_post_den_vals
end

"""
        function inferparameters!(model,samples,trans,psgld_stepsize::NamedTuple;serowaningrate = 0.0,beta = 0.99,λ = 1e-5)

                Infer the unknown transmission and observation parameters for the `model` county. Inference is performed by drawing `samples` number of replicates from the posterior distribution with log-density
    given in the model with sero-waning rate fixed. The parameters and their transformation are encoded in `trans`. Initial parameter search has a pre-optimisation using the default adaptive differential evolution algorithm implemented by the `BlackBoxOptim` package,
MCMC draws are done via PSGLD. Inference output is stored in-place.
"""
function inferparameters!(model::CoVAreaModel,samples,trans,psgld_stepsize::NamedTuple;serowaningrate = 0.0,VE_acquisition=zeros(3),VE_infectiousness=zeros(3),VE_severe_disease_risk=zeros(3),beta = 0.99,λ = 1e-5)
	q₀ = find_initial_MAP(model,trans;serowaningrate =serowaningrate,VE_acquisition=VE_acquisition,VE_infectiousness=VE_infectiousness,VE_severe_disease_risk=VE_severe_disease_risk)
	PSGLD!(model,samples,trans,psgld_stepsize,q₀;serowaningrate =serowaningrate,VE_acquisition=VE_acquisition,VE_infectiousness=VE_infectiousness,VE_severe_disease_risk=VE_severe_disease_risk,beta = beta,λ = λ)
	return nothing
end

"""
        function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal,q₀;serowaningrate = 1/180,num_chains = 1)
Infer the unknown transmission and observation parameters for the `k`th county in by `country_model`. Inference is performed by drawing `samples` number of replicates from the posterior distribution with log-density
    given in the model with sero-waning rate fixed. The parameters and their transformation are encoded in `trans`. Initial parameter search starts at `q₀` in the transformed domain defined by `trans`.
    whilst the initial kinetic energy search begins at `D` and the HMC step size is fixed to be `stepsize`. Inference output is stored in-place.
    Serowaning rate can be added, as well as additional chains for cross-comparison and validation.
"""
function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal,q₀;serowaningrate = 1/365,num_chains = 1)
        println("Starting MCMC parameter inference for $(model.areaname)")
        n = length(q₀)
        l_area(x) = model.log_likelihood(x,model,serowaningrate) + model.log_priors(x)
		ℓ = TransformedLogDensity(trans, l_area)#transformed log-likelihood
        ∇ℓ = LogDensityProblems.ADgradient(:ForwardDiff, ℓ)#transformed log-likelihood gradient wrt the parameters
		#θ = copy(q₀)
		#f(x) = - transform_logdensity(trans, l_area,x)
		#∇ℓ =  ForwardDiff.gradient(f,θ)

        results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, samples,
                                initialization = (q = q₀,κ=GaussianKineticEnergy(D),ϵ = stepsize),
                                warmup_stages = fixed_stepsize_warmup_stages(M=Symmetric),
                                reporter = NoProgressReport())

        transformed_results = TransformVariables.transform.(trans,results.chain)
        val = zeros(length(transformed_results),length(transformed_results[1]),1)

        for i = 1:size(val,1),j = 1:size(val,2)
                val[i,j,1] = transformed_results[i][j]
        end

        for chain_reps = 2:num_chains
                println("Running additional chain $(chain_reps) for $(model.areaname)")
                results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, samples,
                                        initialization = (q = q₀,κ=GaussianKineticEnergy(D),ϵ = stepsize),
                                        warmup_stages = fixed_stepsize_warmup_stages(M=Symmetric),
                                        reporter = NoProgressReport())
                transformed_results = TransformVariables.transform.(trans,results.chain)
                for i = 1:size(val,1),j = 1:size(val,2)
                        val[i,j,chain_reps] = transformed_results[i][j]
                end
        end

        chn = Chains(val,[String(k) for k in keys(transformed_results[1])])

        model.MCMC_results = MCMCResults(chn,
                                        [l_area(transformed_results[i]) - model.log_priors(transformed_results[i]) for i = 1:length(transformed_results)],
                                        results.tree_statistics)

        return nothing
end

"""
        function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal,searchrange;serowaningrate = 1/180,num_chains = 1)
Infer the unknown transmission and observation parameters for the `k`th county in by `country_model`. Inference is performed by drawing `samples` number of replicates from the posterior distribution with log-density
    given in the model with sero-waning rate fixed. The parameters and their transformation are encoded in `trans`. Initial parameter search has a pre-optimisation using the default adaptive differential evolution algorithm implemented by the `BlackBoxOptim` package,
    whilst the initial kinetic energy search begins at `D` and the HMC step size is fixed to be `stepsize`. Inference output is stored in-place.
    serowaning rate can be added, as well as additional chains for cross-comparison and validation.
"""
function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal;serowaningrate = 1/365,num_chains = 1)
        println("Searching for a good initial condition for $(model.areaname)")
        l_area(x) = model.log_likelihood(x,model,serowaningrate) + model.log_priors(x)
        f(x) = -transform_logdensity(trans, l_area,x)
        searchrange = fill((-3.,3.),TransformVariables.dimension(trans))
        res = bboptimize(f; SearchRange = searchrange,PopulationSize=100,MaxSteps=30000,TraceMode = :silent)
        q₀ = best_candidate(res)
        inferparameters!(model,samples,trans,stepsize,D,q₀;serowaningrate=serowaningrate,num_chains=num_chains)
        return nothing
end


text:
struct TuneDPWSim
    mdp
    seed::Int
    n_iters::Int
    max_steps::Int
    k_action::Float64
    alpha_action::Float64
    ec::Float64
    k_state::Float64
    alpha_state::Float64
    metadata::Dict
end
function TuneDPWSim()
end
function POMDPs.simulate(sim::TuneDPWSim)
    rng = MersenneTwister(sim.seed)
    mdp = sim.mdp
    b = DPWBandit(mdp,
                     k_action=sim.k_action,
                     alpha_action=sim.alpha_action,
                     exploration_constant=sim.ec,
                    )
    solver = ModularSolver(mdp, b, sim.seed; n_iterations=sim.n_iters,
                     k_state=sim.k_state,
                     alpha_state=sim.alpha_state)
    policy = solve(solver, mdp)
    s0 = initial_state(mdp, rng)
    hr = HistoryRecorder(; max_steps=sim.max_steps, rng=rng)
    h = simulate(hr, mdp, policy, s0)
    r_total = sum(h.reward_hist)
    r_total
end

@with_kw struct TuneDPW
    mdp=Walk1D()
    n_seeds=5
    n_iters=[10,25,50]
    max_steps=20
    k_actions=[0.0,1.0]
    alpha_actions=[0.0,1.0]
    ecs=[0.0,1.0]
    k_states=[0.0,1.0]
    alpha_states=[0.0,1.0]
end
struct TuneDPWResult
    title::String
    data::DataFrame
end
TuneDPW(::Type{Walk1D}) = TuneDPW(Walk1D(), 200,
                                       [300],  #n_iters
                                       20, #max_steps
                                       [0.1,0.5,1.0,2.0,5.0,10.0,20.0], #k_actions
                                       [0.1,0.25,0.5,0.75,0.9], #alpha_actions
                                       [0.25,0.5,0.75,1.0], #ec
                                       [0.1,0.5,1.0,2.0,5.0,10.0,20.0], #k_states
                                       [0.0,0.1,0.25,0.5,0.75,0.9]) #alpha_states
TuneDPW(::Type{VDPTagMDP}) = TuneDPW(VDPTagMDP(), 500,
                                       [500],  #n_iters
                                       10, #max_steps
                                       [0.1,0.5,1.0,2.0,5.0,10.0], #k_actions
                                       [0.1,0.25,0.5,0.75,0.9], #alpha_actions
                                       [0.25,0.5,0.75,1.0], #ec
                                       [0.1,0.5,1.0,2.0,5.0,10.0], #k_states
                                       [0.1,0.25,0.5,0.75,0.9]) #alpha_states
TuneDPW(::Type{MountainCar}) = TuneDPW(MountainCar(), 20,
                                       [1000],  #n_iters
                                       20, #max_steps
                                       [0.1,0.5,1.0,2.0,5.0,10.0], #k_actions
                                       [0.1,0.25,0.5,0.75,0.9], #alpha_actions
                                       [0.25,0.5,0.75,1.0], #ec
                                       [0.1,0.5,1.0,2.0,5.0,10.0], #k_states
                                       [0.1,0.25,0.5,0.75,0.9]) #alpha_states
TuneDPW(::Type{LightDark2D}) = TuneDPW(LightDark2D(), 200,
                                       [500],  #n_iters
                                       10, #max_steps
                                       [0.1,0.5,1.0,2.0,5.0,10.0], #k_actions
                                       [0.1,0.25,0.5,0.75,0.9], #alpha_actions
                                       [0.25,0.5,0.75,1.0], #ec
                                       [0.1,0.5,1.0,2.0,5.0,10.0], #k_states
                                       [0.1,0.25,0.5,0.75,0.9]) #alpha_states

function generate_sim_q(study::TuneDPW)
    q = []
    for seed in 1:study.n_seeds, n_iters in study.n_iters, k_action in study.k_actions, alpha_action in study.alpha_actions, ec in study.ecs, k_state in study.k_states, alpha_state in study.alpha_states
        push!(q, TuneDPWSim(study.mdp, seed, n_iters, study.max_steps, k_action, alpha_action, ec, k_state, alpha_state, Dict()))
    end
    q
end
function run_study(study::TuneDPW)
    df = DataFrame([String, Int, Int, Int, Float64, Float64, Float64, Float64, Float64, Float64], [:alg, :seed, :n_iters, :max_steps, :k_action, :alpha_action, :ec, :k_state, :alpha_state, :r_total], 0)
    q = generate_sim_q(study)
    result = pmap(POMDPs.simulate, q)
    for (r, sim) in zip(result, q)
        push!(df, ["DPWBandit", sim.seed, sim.n_iters, sim.max_steps, sim.k_action, sim.alpha_action, sim.ec, sim.k_state, sim.alpha_state, r])
    end
    TuneDPWResult(string(study.mdp), df)
end
function report(result::TuneDPWResult)
    df = aggregate(result.data, [:alg, :max_steps, :k_action, :alpha_action, :ec, :k_state, :alpha_state], [mean, std])[[:alg, :k_action, :alpha_action, :ec, :k_state, :alpha_state, :r_total_mean, :r_total_std]]
    dd = sort(df[:], :r_total_mean, rev=true)
    dd
end


text:
using GitHub, JSON, HTTP, MbedTLS
using Dates, Test, Base64
using GitHub: Branch, name
using GitHub.Checks

include("ghtype_tests.jl")
include("event_tests.jl")
include("read_only_api_tests.jl")
include("auth_tests.jl")

@testset "SSH keygen" begin
    pubkey, privkey = GitHub.genkeys(keycomment="GitHub.jl")
    @test endswith(pubkey, "GitHub.jl")
    @test isa(privkey, String)
end


text:
abstract type AbstractContainerView{T} end

struct ContainerView{T} <: AbstractContainerView{T}
    model::EnergyModel
    devices::Dict{Symbol,T}
end

struct SubContainerView{T} <: AbstractContainerView{T}
    model::EnergyModel
    devices::Dict{Symbol,T}
    buses::Axis
end

tupleaxisvals(c::Bus, A) = A.val
tupleaxisvals(c) = tupleaxisvals(c, axis(c))
tupleaxisvals(c, A) = AxisArrays.CategoricalVector(tuple.(c.class, A.val))

_maybesubset(cv::ContainerView, c, A) = A
_maybesubset(cv::SubContainerView, c::Component, A) = A[findall(in(cv.buses), c), ntuple(i->:,ndims(A)-1)...]

_wrapaxis(cv::AbstractContainerView{T}, a) where T = Axis{Symbol(naming(T))}(a)
function _wrapaxisarray(cv::AbstractContainerView, c::Component, a::AxisArray)
    a = _maybesubset(cv, c, a)
    ax = first(AxisArrays.axes(a))
    AxisArray(a, _wrapaxis(cv, tupleaxisvals(c, ax)))
end
function _wrapaxisarray(cv::AbstractContainerView, c::Component, a::Union{Float64,Symbol})
    ax = _maybesubset(cv, c, tupleaxisvals(c))
    AxisArray(fill(a, length(ax)), _wrapaxis(cv, ax))
end

devices(cv::AbstractContainerView) = values(cv.devices)
devices(cv::AbstractContainerView, T::Type{<:Device}) = (c for c = devices(cv) if isa(c, T))

axis(cv::AbstractContainerView) = _wrapaxis(cv, vcat((_maybesubset(cv, c, tupleaxisvals(c)) for c = devices(cv))...))
axis(cv::AbstractContainerView, T::Type{<:Component}) = axis(m[T])

# Base.cat(cv::AbstractContainerView, cs::Vector{<:Component}, as::Vector) = @consense(as, "Single values do not agree")
"""
    cat(containerview, components, attrvalues)

Combine the `attrvalues` for the `components` in a `containerview` to a single AxisArray.
"""
Base.cat(cv::AbstractContainerView, components::Vector{<:Component}, attrvalues::Vector) =
    cat((_wrapaxisarray(cv, c, a) for (c, a) = zip(components, attrvalues))...; dims=1)

mapcat(f::Function, cv::AbstractContainerView) = cat(cv, collect(devices(cv)), map(f, devices(cv)))

Base.get(cv::AbstractContainerView, i) = mapcat(c->get(c, i), cv)
Base.get(cv::AbstractContainerView, i, axes...) = WrappedArray(get(cv, i), axes...)

Base.getindex(cv::AbstractContainerView, i) = get(cv, i)

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
#Benchmarking of SGU variants
using DifferentiableStateSpaceModels, BenchmarkTools
using DifferentiableStateSpaceModels.Examples
const SGU = BenchmarkGroup()
const SGU["sgu"] = BenchmarkGroup()

# Prepare the model for the sgu observables
const m_sgu = @include_example_module(Examples.sgu)
# Basic Steady State
const p_f_sgu = nothing
const p_d_sgu = (γ=2.0, ω=1.455, ρ=0.42, σe=0.0129, δ=0.1, ψ=0.000742, α=0.32, ϕ=0.028,
β=1.0 / (1.0 + 0.04), r_w=0.04, d_bar=0.7442)
const c_sgu_1 = SolverCache(m_sgu, Val(1), p_d_sgu)
const c_sgu_2 = SolverCache(m_sgu, Val(2), p_d_sgu)

# The derivatives caches need to be already filled in case they are done out of order
const c_sgu_1_p = SolverCache(m_sgu, Val(1), p_d_sgu)
const c_sgu_2_p = SolverCache(m_sgu, Val(2), p_d_sgu)
generate_perturbation(m_sgu, p_d_sgu, p_f_sgu, Val(1); cache = c_sgu_1_p)
generate_perturbation(m_sgu, p_d_sgu, p_f_sgu, Val(2); cache = c_sgu_2_p)

# Can't hurt to  call the derivatives to precompile as well, eventhough they will be written over
generate_perturbation_derivatives!(m_sgu, p_d_sgu, p_f_sgu, c_sgu_1_p)
generate_perturbation_derivatives!(m_sgu, p_d_sgu, p_f_sgu, c_sgu_2_p)

# sgu
const SGU["sgu"]["first_order"] =
    @benchmarkable generate_perturbation($m_sgu, $p_d_sgu, $p_f_sgu, Val(1); cache = $c_sgu_1)
const SGU["sgu"]["first_order_p"] =
    @benchmarkable generate_perturbation_derivatives!($m_sgu, $p_d_sgu, $p_f_sgu, $c_sgu_1_p)
const SGU["sgu"]["second_order"] =
    @benchmarkable generate_perturbation($m_sgu, $p_d_sgu, $p_f_sgu, Val(2); cache = $c_sgu_2)
const SGU["sgu"]["second_order_p"] =
    @benchmarkable generate_perturbation_derivatives!($m_sgu, $p_d_sgu, $p_f_sgu, $c_sgu_2_p)
const SGU["sgu"]["second_order_no_cache"] =
    @benchmarkable generate_perturbation($m_sgu, $p_d_sgu, $p_f_sgu, Val(2)) # just to see how long the cache takes
const SGU["sgu"]["first_order_solver_cache"] = @benchmarkable SolverCache($m_sgu, Val(1), $p_d_sgu)
const SGU["sgu"]["second_order_solver_cache"] = @benchmarkable SolverCache($m_sgu, Val(2), $p_d_sgu)

# return for the suite
SGU

# Uncomment to run on its own, otherwise nests into full benchmarking
# results = run(SGU; verbose = true)
# results["sgu"] # etc.


text:
using CompScienceMeshes

function mesh_box_width_lid(size, h1, h2)
    width, height, depth = size, size, size

    m = -meshrectangle(width, height, h1, 3)
    m_front = rotate(m, 0.5*pi*[1,0,0])
    m_left  = rotate(m, -0.5*pi*[0,1,0])
    m_back  = -translate(m_front, [0,height,0])
    m_right = -translate(m_left, [width,0,0])
    m_bottom = -m
    m_open = weld(m_front, m_left, m_back, m_right, m_bottom)

    m_fine = -meshrectangle(width, height, h2)
    m_top  = translate(m_fine, [0,0,height])

    # m_coarse = -meshrectangle(width, height, h1)
    # m_coarse = translate(m_coarse, [0,0,height])
    # m_conf = weld(m_open, m_coarse)

    return m_open, m_top
end

"""
    meshcube(size, h) -> mesh

Mesh a conforming cube.
"""
meshcube(size, h) =  weld(mesh_box_width_lid(size,h,h)...)


text:
import OpenNI
using OpenNI

function distance(pnt1, pnt2, norm=2)
    sum = 0
    for i=1:length(pnt1)
        sum += (pnt1[i] - pnt2[i]) ^ norm
    end
    return sum ^ (1 / norm)
end

###Record only once skeleton is detected and then the user puts their right hand by torso

g_elapseSinceLastEvent = 0
g_addedNodes = false
g_userFound = false
g_startedRecording = false
g_stoppedRecording = false

#Init without any preloaded nodes
print("Creating context...")

context = NIContext()
InitContext(context)

println("Done!")
print("Attaching generators...")

#Create an ImageGenerator node
imagegen = NIImageGenerator()
CreateImageGenerator(imagegen, context)

#Create a DepthGenerator node
depthgen = NIDepthGenerator()
CreateDepthGenerator(depthgen, context)

#Now create a UserGenerator node
usergen = NIUserGenerator()
CreateUserGenerator(usergen, context)

#Function for callbacks (new and lost user)
function userFound()
    println("User recognized registered!")
end

function userLost()
    g_stoppedRecording = true
    println("User lost. Stopping recording")
end

RegisterUserCallbacks(usergen, userFound, userLost)

println("Done!")
print("Creating Recorder...")

#Create recorder
recorder = NIRecorder()
CreateRecorder(recorder, context)

#specify the output of the recorder
SetRecorderDestination(recorder,  joinpath(@__DIR__, "../../output.oni"))

println("Done!")

#Begin generators
StartGeneratingAll(context)

print("Starting device and recording...")

#init array of user IDs
userIDs = Vector{UInt32}(2)

tic()
#start program loop
while !g_stoppedRecording
    #each frame update, but wait for usergen to have new data
    WaitOneUpdateAll(context, usergen)

    userIDs = GetUsers(usergen, userIDs)

    for i=1:length(userIDs)
        #If tracking, print head coordinates
        if (IsTrackingUser(usergen, userIDs[i]))

            #attach generators to recorder, do this now so that the skeleton is detected
            if (!g_addedNodes)
                print("Adding nodes to record...")
                AddNodeToRecording(recorder, imagegen)
                AddNodeToRecording(recorder, depthgen)
                g_addedNodes = true
                println("Done!")
            end

            skeleton = NISkeleton(usergen, userIDs[i])
            jntHandID = Int(OpenNI.NI_SKEL_LEFT_HAND) #mirrored
            jntTorsoID = Int(OpenNI.NI_SKEL_TORSO)
            dist = distance(skeleton.jointPos[jntHandID, :], skeleton.jointPos[jntTorsoID, :])

            g_elapseSinceLastEvent += toq()
            tic()

            if (!g_userFound && dist < 250)
                g_userFound = true
                print("User struck pose, begin motion in 2 seconds...")
                g_elapseSinceLastEvent = 0
                tic()
            end

            if (g_userFound && !g_startedRecording && g_elapseSinceLastEvent > 2)
                println("Starting recording...")
                g_startedRecording = true
                g_elapseSinceLastEvent = 0
                tic()
            end

            if (g_startedRecording && g_elapseSinceLastEvent > 5)
                println("Stopping recording and saving to file")
                g_stoppedRecording = true
            end
        end #IsTrackingUser
    end #for

    #record the frame
    Record(recorder)
end

#cleanup recorder
RemoveNodeFromRecording(recorder, imagegen)
RemoveNodeFromRecording(recorder, depthgen)

#force cleanup to print to file?
finalize(recorder)

println("Finished recording")


text:
using Pkg

dependencies = [
    "PyCall",
    "ColorSchemes",
    "ColorTypes",
    "Images",
    "Luxor",
    "CSV",
    "DataFrames",
    "Optim",
    "JuMP",
    "Clp",
    "Cbc",
    "GLPK",
    "PackageCompiler",
    "MultivariateStats",
    "Glob",
    "ImageIO",
    "Ipopt"]

Pkg.add(dependencies)

max_stars_repo_name: nmaterise/optmap

max_issues_repo_name: nmaterise/optmap

max_forks_repo_name: nmaterise/optmap


text:
# # Simulation Chain: Inverted Coax Detector

using Plots
using SolidStateDetectors
using Unitful

T = Float32
sim = Simulation{T}(SSD_examples[:InvertedCoax])

plot(sim.detector, size = (700, 700))
#jl savefig("tutorial_det.pdf") # hide
#md savefig("tutorial_det.pdf") # hide
#md savefig("tutorial_det.svg"); nothing # hide
#md # [![tutorial_det](tutorial_det.svg)](tutorial_det.pdf)

# One can also have a look at how the initial conditions look like on the grid (its starts with a very coarse grid):

apply_initial_state!(sim, ElectricPotential) # optional
plot(
    plot(sim.electric_potential), # initial electric potential (boundary conditions)
    plot(sim.point_types), # map of different point types: fixed point / inside or outside detector volume / depleted/undepleted
    plot(sim.q_eff_imp), # charge density distribution
    plot(sim.ϵ_r), # dielectric distribution
    layout = (1, 4), size = (1600, 500)
)
#jl savefig("tutorial_initial_condition.pdf") # hide
#md savefig("tutorial_initial_condition.pdf") # hide
#md savefig("tutorial_initial_condition.svg"); nothing # hide
#md # [![tutorial_initial_condition](tutorial_initial_condition.svg)](tutorial_initial_condition.pdf)

# Next, calculate the electric potential:

calculate_electric_potential!( sim,
                               refinement_limits = [0.2, 0.1, 0.05, 0.01])

plot(
    plot(sim.electric_potential, φ = 20), # initial electric potential (boundary conditions)
    plot(sim.point_types), # map of different point types: fixed point / inside or outside detector volume / depleted/undepleted
    plot(sim.q_eff_imp), # charge density distribution
    plot(sim.ϵ_r), # dielectric distribution
    layout = (1, 4), size = (1600, 500)
)
#jl savefig("tutorial_calculated_potential.pdf") # hide
#md savefig("tutorial_calculated_potential.pdf") # hide
#md savefig("tutorial_calculated_potential.svg"); nothing # hide
#md # [![tutorial_calculated_potential](tutorial_calculated_potential.svg)](tutorial_calculated_potential.pdf)

# SolidStateDetectors.jl supports active (i.e. depleted) volume calculation:

is_depleted(sim.point_types)

#

get_active_volume(sim.point_types) # approximation (sum of the volume of cells marked as depleted)

# ## Partially depleted detectors

# SolidStateDetectors.jl can also calculate the electric potential of a partially depleted detector:

sim_undep = deepcopy(sim)
sim_undep.detector = SolidStateDetector(sim_undep.detector, contact_id = 2, contact_potential = 500); # V  <-- Bias Voltage of Mantle

calculate_electric_potential!( sim_undep,
                               depletion_handling = true,
                               convergence_limit = 1e-6,
                               refinement_limits = [0.2, 0.1, 0.05, 0.01],
                               verbose = false)

plot(
    plot(sim_undep.electric_potential),
    plot(sim_undep.point_types),
    plot(sim_undep.imp_scale),
    layout = (1, 3), size = (1200, 600)
)
#jl savefig("tutorial_calculated_potential_undep.pdf") # hide
#md savefig("tutorial_calculated_potential_undep.pdf") # hide
#md savefig("tutorial_calculated_potential_undep.svg"); nothing # hide
#md # [![tutorial_calculated_potential_undep](tutorial_calculated_potential_undep.svg)](tutorial_calculated_potential_undep.pdf)

#

is_depleted(sim_undep.point_types)

# Compare both volumes:

println("Depleted:   ", get_active_volume(sim.point_types))
println("Undepleted: ", get_active_volume(sim_undep.point_types));

# ## Electric field calculation

# Calculate the electric field of the fully depleted detector, given the already calculated electric potential:

calculate_electric_field!(sim, n_points_in_φ = 72)

plot(sim.electric_field, full_det = true, φ = 0.0, size = (700, 700))
plot_electric_fieldlines!(sim, full_det = true, φ = 0.0)
#jl savefig("tutorial_electric_field.pdf") # hide
#md savefig("tutorial_electric_field.pdf") # hide
#md savefig("tutorial_electric_field.svg"); nothing # hide
#md # [![tutorial_electric_field](tutorial_electric_field.svg)](tutorial_electric_field.pdf)

# ## Simulation of charge drifts

# Any charge drift model can be used for the calculation of the electric field. If no model is explicitely given, the `ElectricFieldChargeDriftModel` is used. Other configurations are saved in their configuration files and can be found under:

# `<package_directory>/examples/example_config_files/ADLChargeDriftModel/<config_filename>.yaml.`

# Set the charge drift model of the simulation:

charge_drift_model = ADLChargeDriftModel()
sim.detector = SolidStateDetector(sim.detector, charge_drift_model)

# Now, let's create an "random" multi-site event:

starting_positions = [ CylindricalPoint{T}( 0.020, deg2rad(10), 0.015 ),
                       CylindricalPoint{T}( 0.015, deg2rad(20), 0.045 ),
                       CylindricalPoint{T}( 0.022, deg2rad(35), 0.025 ) ]
energy_depos = T[1460, 609, 1000] * u"keV" # are needed later in the signal generation

evt = Event(starting_positions, energy_depos);

time_step = 5u"ns"
drift_charges!(evt, sim, Δt = time_step)

plot(sim.detector, size = (700, 700))
plot!(evt.drift_paths)
#jl savefig("tutorial_drift_paths.pdf") # hide
#md savefig("tutorial_drift_paths.pdf") # hide
#md savefig("tutorial_drift_paths.svg"); nothing # hide
#md # [![tutorial_drift_paths](tutorial_drift_paths.svg)](tutorial_drift_paths.pdf)

# ## Weighting potential calculation

# We need weighting potentials to simulate the detector charge signal induced by drifting charges. We'll calculate the weighting potential for the point contact and the outer shell of the detector:

for contact in sim.detector.contacts
    calculate_weighting_potential!(sim, contact.id, refinement_limits = [0.2, 0.1, 0.05, 0.01], n_points_in_φ = 2, verbose = false)
end

plot(
    plot(sim.weighting_potentials[1]),
    plot(sim.weighting_potentials[2]),
    size = (900, 700)
)
#jl savefig("tutorial_weighting_potentials.pdf") # hide
#md savefig("tutorial_weighting_potentials.pdf") # hide
#md savefig("tutorial_weighting_potentials.svg"); nothing # hide
#md # [![tutorial_weighting_potentials](tutorial_weighting_potentials.svg)](tutorial_weighting_potentials.pdf)

# ## Detector Capacitance Matrix

# After the weighting potentials are calculated, the detector capacitance matrix can be calculated in the **Maxwell Capacitance Matrix Notation**:

calculate_capacitance_matrix(sim)

# See [Capacitances](@ref) for more information.

# ## Detector waveform generation

# Given an interaction at an arbitrary point in the detector, we can now simulate the charge drift and the resulting detector charge signals (e.g. at the point contact):

simulate!(evt, sim) # drift_charges + signal generation of all channels

p_pc_signal = plot( evt.waveforms[1], lw = 1.5, xlims = (0, 1100), xlabel = "Time", unitformat = :slash,
                    legend = false, tickfontsize = 12, ylabel = "Charge", guidefontsize = 14)
#jl savefig("tutorial_waveforms.pdf") # hide
#md savefig("tutorial_waveforms.pdf") # hide
#md savefig("tutorial_waveforms.svg"); nothing # hide
#md # [![tutorial_waveforms](tutorial_waveforms.svg)](tutorial_waveforms.pdf)

# SolidStateDetectors.jl also allows to separate the waveform into the two contributions from electrons and holes

contact_id = 1
plot_electron_and_hole_contribution(evt, sim, contact_id, xlims = (0, 1100), xlabel = "Time",
                    legend = :topleft, tickfontsize = 12, ylabel = "Charge", guidefontsize = 14)
#jl savefig("tutorial_waveform_contributions.pdf") # hide
#md savefig("tutorial_waveform_contributions.pdf") # hide
#md savefig("tutorial_waveform_contributions.svg"); nothing # hide
#md # [![tutorial_waveform_contributions](tutorial_waveform_contributions.svg)](tutorial_waveform_contributions.pdf)


text:
module MathJaxRenderer

using Artifacts, Librsvg_jll, NodeJS, JSON3

export Math

"""
    Math(src::AbstractString)

Wrap a string `src` containing LaTeX math syntax for conversion to one of the
following formats: `svg`, `png`, `pdf`, or `eps`. The MathJax NodeJS library
is used to typeset `src`.

When writing to either a file (via `write(filename, math)`) or an `IO` object
(via `show(io, mime, math)`) the following keyword options are available for
either `show` or `write`:

  - The following apply to the conversion from `tex` to `svg`:

      - `inline = false`, whether to typeset as inline or display math.
      - `em = 16`, a number giving the number of pixels in an `em` for the surrounding font.
      - `ex = 8`, a number giving the number of pixels in an `ex` for the surrounding font.
      - `width = 80em`, a number giving the width of the container, in pixels.

  - and the following apply to conversion from `tex` to `svg`, `png`, `pdf`, and `eps`:

      - `dpi_x = 90`, pixels per inch (horizontal).
      - `dpi_y = 90`, pixels per inch (vertical).
      - `x_zoom = 1.0`, zoom factor (horizontal).
      - `y_zoom = 1.0`, zoom factor (vertical).
      - `zoom = 1.0`, zoom factor (both horizontal and vertical).
      - `keep_aspect_ratio = false`, whether to preserve the aspect ratio.
      - `background_color = :white`, set the background color, CSS color names and syntax.
      - `unlimited = false`, allow for huge SVG files.

When writing to file with `write` the follow extensions are supported:

```julia-repl
julia> m = Math("\\sqrt{x^2 + y^2}");

julia> write("math.svg", m); # MIME("image/svg+xml")

julia> write("math.png", m); # MIME("image/png")

julia> write("math.pdf", m); # MIME("application/pdf")

julia> write("math.eps", m); # MIME("application/postscript")

```

and when using `show` use the noted `MIME` types for each extension.
"""
struct Math
    src::String
end

Base.read(fn::String, ::Type{Math}) = Math(read(fn, String))
Base.write(fn::String, d::Math; kws...) = open(io -> show(io, mimetype(fn)(), d; kws...), fn, "w")

Base.show(io::IO, d::Math) = print(io, "$Math(...)")

### Keep next section in sync.

const SVG = MIME"image/svg+xml"
const PNG = MIME"image/png"
const PDF = MIME"application/pdf"
const EPS = MIME"application/postscript"

const SUPPORTED_MIMES = Union{SVG,PNG,PDF,EPS}

extension(::SVG) = "svg"
extension(::PNG) = "png"
extension(::PDF) = "pdf"
extension(::EPS) = "eps"
extension(other) = error("unknown extension '$other'.")

function mimetype(filename::AbstractString)
    _, ext = splitext(filename)
    return mimetype(Val{Symbol(ext[2:end])}())
end
mimetype(::Val{:svg}) = SVG
mimetype(::Val{:png}) = PNG
mimetype(::Val{:pdf}) = PDF
mimetype(::Val{:eps}) = EPS
mimetype(::Val{s}) where s = error("unknown format '$s'.")

###

Base.show(io::IO, mime::SUPPORTED_MIMES, d::Math; kws...) = write(io, converter(d, mime; kws...))

mathjax_bin() = joinpath(artifact"mathjax", "index.js")
mathjax() = `$(nodejs_cmd()) $(mathjax_bin())`

function exec(cmd::Cmd, input::IOBuffer)
    output, errors = IOBuffer(), IOBuffer()
    yes = success(pipeline(cmd; stdout = output, stdin = input, stderr = errors))
    return (yes = yes, io = yes ? output : errors)
end
exec(cmd::Cmd, input="") = exec(cmd, IOBuffer(input))

function extract_svg(s::String)
    fn(str, a::UnitRange, b::UnitRange) = String(str[a[begin]:b[end]])
    fn(str, ::Any, ::Any) = str
    return fn(s, findfirst("<svg", s), findlast("</svg>", s))
end
extract_svg(vec::Vector{UInt8}) = extract_svg(String(vec))

function svg_converter(
    f::Math;
    inline = false,
    em = 16,
    ex = 8,
    width = 80em,
    kws...,
)
    conf = (
        src = f.src,
        config = (
            display = !inline,
            em = em,
            ex = ex,
            containerWidth = width,
        ),
    )
    str = take!(check(exec(mathjax(), JSON3.write(conf))).io)
    return IOBuffer(extract_svg(str))
end

function converter(
    f::Math,
    m::SUPPORTED_MIMES;
    dpi_x = 90,
    dpi_y = 90,
    x_zoom = 1.0,
    y_zoom = 1.0,
    zoom = 1.0,
    keep_aspect_ratio = false,
    background_color = :white,
    unlimited = false,
    kws...,
)
    Librsvg_jll.rsvg_convert() do bin
        ext = extension(m)
        cmd = [
            bin,
            "--format=$(extension(m))",
            "--dpi-x=$dpi_x",
            "--dpi-y=$dpi_y",
            "--x-zoom=$x_zoom",
            "--y-zoom=$y_zoom",
            "--zoom=$zoom",
            "--background-color=$background_color",
        ]
        keep_aspect_ratio && push!(cmd, "--keep-aspect-ratio")
        unlimited && push!(cmd, "--unlimited")
        take!(check(exec(Cmd(cmd), svg_converter(f; kws...))).io)
    end
end
check(result) = result.yes ? result : error(String(take!(result.io)))

end # module