text:
"""
MCTS solver with Importance Sampling

Fields:
    depth::Int64
        Maximum rollout horizon and tree depth.
        default: 10
    exploration_constant::Float64
        Specified how much the solver should explore.
        In the UCB equation, Q + c*sqrt(log(t/N)), c is the exploration constant.
        default: 1.0
    n_iterations::Int64
        Number of iterations during each action() call.
        default: 100
    max_time::Float64
        Maximum amount of CPU time spent iterating through simulations.
        default: Inf
    k_action::Float64
    alpha_action::Float64
    k_state::Float64
    alpha_state::Float64
        These constants control the double progressive widening. A new state
        or action will be added if the number of children is less than or equal to kN^alpha.
        defaults: k:10, alpha:0.5
    virtual_loss::Float64
        A temporary loss added to the UCB score for nodes that are currently being
        evaluated by some threads. This can be used to encourage threads to explore
        broader parts of the search space. Relevant when running MCTS with multiple threads.
        default: 0.0
    keep_tree::Bool
        If true, store the tree in the planner for reuse at the next timestep (and every time it is used in the future). There is a computational cost for maintaining the state dictionary necessary for this.
        default: false
    enable_action_pw::Bool
        If true, enable progressive widening on the action space; if false just use the whole action space.
        default: true
    enable_state_pw::Bool
        If true, enable progressive widening on the state space; if false just use the single next state (for deterministic problems).
        default: true
    tree_in_info::Bool
        If true, return the tree in the info dict when action_info is called. False by default because it can use a lot of memory if histories are being saved.
        default: false
    rng::AbstractRNG
        Random number generator
    init_Q::Any
        Function, object, or number used to set the initial Q(s,a) value at a new node.
        If this is a function `f`, `f(mdp, s, a)` will be called to set the value.
        If this is an object `o`, `init_Q(o, mdp, s, a)` will be called.
        If this is a number, Q will always be set to that number.
        default: 0.0
    init_N::Any
        Function, object, or number used to set the initial N(s,a) value at a new node.
        If this is a function `f`, `f(mdp, s, a)` will be called to set the value.
        If this is an object `o`, `init_N(o, mdp, s, a)` will be called.
        If this is a number, N will always be set to that number.
        default: 0
    next_action::Any
        Function or object used to choose the next action to be considered for progressive widening.
        The next action is determined based on the MDP, the state, `s`, and the current `DPWStateNode`, `snode`.
        If this is a function `f`, `f(mdp, s, snode)` will be called to set the value.
        If this is an object `o`, `next_action(o, mdp, s, snode)` will be called.
        default: RandomActionGenerator(rng)
    default_action::Any
        Function, action, or Policy used to determine the action if POMCP fails with exception `ex`.
        If this is a Function `f`, `f(pomdp, belief, ex)` will be called.
        If this is a Policy `p`, `action(p, belief)` will be called.
        If it is an object `a`, `default_action(a, pomdp, belief, ex)` will be called, and if this method is not implemented, `a` will be returned directly.
        default: `ExceptionRethrow()`
    reset_callback::Function
        Function used to reset/reinitialize the MDP to a given state `s`.
        Useful when the simulator state is not truly separate from the MDP state.
        `f(mdp, s)` will be called.
        default: `(mdp, s)->false` (optimized out)
    show_progress::Bool
        Show progress bar during simulation.
        default: false
"""

mutable struct PISSolver
    depth::Int
    exploration_constant::Float64
    n_iterations::Int
    max_time::Float64
    k_action::Float64
    alpha_action::Float64
    k_state::Float64
    alpha_state::Float64
    virtual_loss::Float64
    nominal_distrib_fn::Function
    cost_reduction::Symbol
    action_selection::Symbol
    keep_tree::Bool
    enable_action_pw::Bool
    enable_state_pw::Bool
    tree_in_info::Bool
    rng::AbstractRNG
    init_Q::Any
    init_N::Any
    next_action::Any
    default_action::Any
    reset_callback::Function
    show_progress::Bool
    timer::Function
    # TODO(kykim): Either put all params in the solver or make them method args.
    α::Float64
end

mutable struct UniformActionGenerator{RNG<:AbstractRNG}
    rng::RNG
end
UniformActionGenerator() = UniformActionGenerator(Random.GLOBAL_RNG)

"""
Use keyword arguments to specify values for the fields.
"""
function PISSolver(;depth::Int=10,
                   exploration_constant::Float64=1.0,
                   n_iterations::Int=100,
                   max_time::Float64=Inf,
                   k_action::Float64=10.0,
                   alpha_action::Float64=0.5,
                   k_state::Float64=10.0,
                   alpha_state::Float64=0.5,
                   virtual_loss::Float64=0.0,
                   nominal_distrib_fn=(mdp, s)->Normal(0, 1),
                   cost_reduction::Symbol=:sum,
                   action_selection::Symbol=:adaptive,
                   keep_tree::Bool=false,
                   enable_action_pw::Bool=true,
                   enable_state_pw::Bool=true,
                   tree_in_info::Bool=false,
                   rng::AbstractRNG=Random.GLOBAL_RNG,
                   init_Q::Any=0.0,
                   init_N::Any=1,
                   next_action::Any=UniformActionGenerator(rng),
                   default_action::Any=nothing,
                   reset_callback::Function=(mdp, s)->false,
                   show_progress::Bool=false,
                   timer=()->1e-9*time_ns(),
                   α::Float64=0.1)
    PISSolver(depth, exploration_constant, n_iterations, max_time, k_action,
              alpha_action, k_state, alpha_state, virtual_loss,
              nominal_distrib_fn, cost_reduction, action_selection, keep_tree,
              enable_action_pw, enable_state_pw, tree_in_info, rng, init_Q,
              init_N, next_action, default_action, reset_callback,
              show_progress, timer, α)
end

mutable struct PISPlanner{P<:Union{MDP,POMDP}, S, A, NA, RCB, RNG}
    solver::PISSolver
    mdp::P
    tree::Union{Nothing, PISTree{S,A}}
    next_action::NA
    reset_callback::RCB
    rng::RNG
end

function PISPlanner(solver::PISSolver, mdp::P) where P<:Union{POMDP,MDP}
    return PISPlanner{P,
                      statetype(P),
                      # TODO(kykim): Temporary hack to get around the RMDP case.
                      # Should be actiontype(P).
                      Any,
                      typeof(solver.next_action),
                      typeof(solver.reset_callback),
                      typeof(solver.rng)}(
                          solver,
                          mdp,
                          nothing,
                          solver.next_action,
                          solver.reset_callback,
                          solver.rng)
end

Random.seed!(p::PISPlanner, seed) = Random.seed!(p.rng, seed)

POMDPs.solve(solver::PISSolver, mdp::Union{POMDP,MDP}) = PISPlanner(solver, mdp)


text:
"""
Run an intcode program: 1 adds, 2 multiplies, 99 halts
prog: the program to run as an array of integers, will be modified in place
"""
function run_intcode!(prog)
    pc = 1
    while prog[pc] != 99
        # Remember to add 1 for Julia 1-based indices!
        arg1_index = prog[pc+1] + 1
        arg2_index = prog[pc+2] + 1
        res_index = prog[pc+3] + 1
        if prog[pc] == 1
            # 1 means add
            prog[res_index] = prog[arg1_index] + prog[arg2_index]
        elseif prog[pc] == 2
            # 2 means mul
            prog[res_index] = prog[arg1_index] * prog[arg2_index]
        end
        pc += 4
    end
end

# Part 1: load input and patch it, then say what was in first element after running
program = open(f -> parse.(Int32, split(readline(f), ',')), "day02_input.txt")
patched = copy(program)
patched[2] = 12
patched[3] = 2
run_intcode!(patched)
println("First element: $(patched[1])")

# Part 2: do a stupid grid search to find "noun" and "verb"
for noun in 0:99, verb in 0:99
    test = copy(program)
    test[2] = noun
    test[3] = verb
    run_intcode!(test)
    if test[1] == 19690720
        # This is the one!
        println("Noun=$noun, verb=$verb, 100*noun+verb=$(100*noun + verb)")
    end
end

ext: jl

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
module Human1

    using ProjAssistant
    @gen_sub_proj

    include("HumanGEM/HumanGEM.jl")
    include("HartGEMs/HartGEMs.jl")

    function __init__()
        @create_proj_dirs
    end

end


text:
using CSV, Feather
#using JLD2
#using JLD#, JLSO
using JDF, FileIO, Blosc, StatsPlots, RCall
using DataFrames, WeakRefStrings # required for JLD2, JDF
Blosc.set_num_threads(6)

gen_benchmark(dirpath, largest_file, outpath, data_label; delim = ',', header=true) = begin
    if !isdir(outpath)
        mkpath(outpath)
    end

    csv_read1 = @elapsed df = CSV.read(joinpath(dirpath, largest_file), delim = delim, header = header, threaded=false);
    csv_read2 = @elapsed      CSV.read(joinpath(dirpath, largest_file), delim = delim, header = header, threaded=false);

    csv_write1 = 0
    csv_write2 = 0
    try
        csv_write1 = @elapsed CSV.write(joinpath(outpath, largest_file*".csv"), df);
        csv_write2 = @elapsed CSV.write(joinpath(outpath, largest_file*".csv"), df);
    catch err
    end

    R"""
    library(data.table)
    library(fst)
    # memory.limit(Inf)

    data_table_read1 = system.time(a <- data.table::fread($(joinpath(dirpath, largest_file))))[3]
    data_table_read2 = system.time(data.table::fread($(joinpath(dirpath, largest_file))))[3]

    data_table_write1 = system.time(data.table::fwrite(a, $(joinpath(outpath, largest_file*".data.table.csv"))))[3]
    data_table_write2 = system.time(data.table::fwrite(a, $(joinpath(outpath, largest_file*".data.table.csv"))))[3]

    fst_write1 = system.time(fst::write_fst(a, $(joinpath(outpath, largest_file*".fst"))))[3]
    fst_write2 = system.time(fst::write_fst(a, $(joinpath(outpath, largest_file*".fst"))))[3]

    fst_read1 = system.time(fst::read_fst($(joinpath(outpath, largest_file*".fst"))))[3]
    fst_read2 = system.time(fst::read_fst($(joinpath(outpath, largest_file*".fst"))))[3]

    parquet_r_write1 = system.time(arrow::write_parquet(a, $(joinpath(outpath, largest_file*".parquet"))))[3]
    parquet_r_write2 = system.time(arrow::write_parquet(a, $(joinpath(outpath, largest_file*".parquet"))))[3]

    parquet_r_read1 = system.time(arrow::read_parquet($(joinpath(outpath, largest_file*".parquet"))))[3]
    parquet_r_read2 = system.time(arrow::read_parquet($(joinpath(outpath, largest_file*".parquet"))))[3]
    rm(a)
    gc()
    """

    @rget data_table_read1
    @rget data_table_read2
    @rget data_table_write1
    @rget data_table_write2
    @rget fst_write1
    @rget fst_write2
    @rget fst_read1
    @rget fst_read2
    @rget parquet_r_read1
    @rget parquet_r_read2
    @rget parquet_r_write1
    @rget parquet_r_write2

    # jlso_write1 = 0
    # jlso_write2 = 0
    # try
    #     jlso_write1 = @elapsed JLSO.save(joinpath(outpath, largest_file*".jlso"), df);
    #     jlso_write2 = @elapsed JLSO.save(joinpath(outpath, largest_file*".jlso"), df);
    # catch err
    # end

    # jld2_write1 = 0
    # jld2_write2 = 0
    # try
    #     jld2_write1 = @elapsed save(joinpath(outpath, largest_file*".jld2"), Dict("df" => df));
    #     jld2_write2 = @elapsed save(joinpath(outpath, largest_file*".jld2"), Dict("df" => df));
    # catch err
    # end

    # jld_write1 = 0
    # jld_write2 = 0
    # try
    #     jld_write1 = @elapsed save(joinpath(outpath, largest_file*".jld"), Dict("df" => df));
    #     jld_write2 = @elapsed save(joinpath(outpath, largest_file*".jld"), Dict("df" => df));
    # catch err
    # end

    jdf_write1 = 0
    jdf_write2 = 0
    try
        jdf_write1 = @elapsed savejdf(joinpath(outpath, largest_file*".jdf"), df);
        jdf_write2 = @elapsed savejdf(joinpath(outpath, largest_file*".jdf"), df);
    catch err

    end

    # # feather can't handle all missing
    for n in names(df)
        if eltype(df[!,n]) == Missing
            println("Removed $n for Feather.jl")
            select!(df, Not(n))
            df[!,n] = Vector{Union{Missing, Bool}}(missing, size(df, 1))
        end
    end

    feather_write1 = 0
    feather_write2 = 0
    try
        feather_write1 = @elapsed Feather.write(joinpath(outpath, largest_file*".feather"), df);
        feather_write2 = @elapsed Feather.write(joinpath(outpath, largest_file*".feather"), df);
    catch err
    end

    #########################################  loading
    feather_read1 = 0
    feather_read2 = 0
    try
        feather_read1 = @elapsed Feather.read(joinpath(outpath, largest_file*".feather"));
        feather_read2 = @elapsed Feather.read(joinpath(outpath, largest_file*".feather"));
    catch err
    end

    # jld_read1 = 0
    # jld_read2 = 0
    # try
    #     jld_read1 = @elapsed load(joinpath(outpath, largest_file*".jld"))["df"];
    #     jld_read2 = @elapsed load(joinpath(outpath, largest_file*".jld"))["df"];
    # catch err
    # end

    # jld2_read1 = 0
    # jld2_read2 = 0
    # try
    #     jld2_read1 = @elapsed load(joinpath(outpath, largest_file*".jld2"))["df"];
    #     jld2_read2 = @elapsed load(joinpath(outpath, largest_file*".jld2"))["df"];
    # catch err
    # end
    #
    # jlso_read1 = 0
    # jlso_read2 = 0
    # try
    #     jlso_read1 = @elapsed JLSO.load(joinpath(outpath, largest_file*".jlso"))["data"];
    #     jlso_read2 = @elapsed JLSO.load(joinpath(outpath, largest_file*".jlso"))["data"];
    # catch err
    # end

    jdf_read1 = 0
    jdf_read2 = 0
    try
        jdf_read1 = @elapsed loadjdf(joinpath(outpath, largest_file*".jdf"));
        jdf_read2 = @elapsed loadjdf(joinpath(outpath, largest_file*".jdf"));
    catch err
    end

    # write_perf = [jdf_write1, jdf_write2, csv_write1, csv_write2, feather_write1, feather_write2, jld2_write1, jld2_write2, jlso_write1, jlso_write2]
    # read_perf = [jdf_read1, jdf_read2, csv_read1, csv_read2, feather_read1, feather_read2, jld2_read1, jld2_read2, jlso_read1, jlso_read2]

    write_perf = [jdf_write1, jdf_write2, csv_write1, csv_write2, feather_write1, feather_write2, data_table_write1, data_table_write2, fst_write1, fst_write2, parquet_r_write1, parquet_r_write2]
    read_perf  = [jdf_read1,  jdf_read2,  csv_read1, csv_read2, feather_read1, feather_read2, data_table_read1, data_table_read2, fst_read1, fst_read2, parquet_r_read1, parquet_r_read2]

    pkgs = repeat(["JDF.jl", "CSV.jl", "Feather.jl", "data.table", "fst", "parquet R"], inner = 2)
    run_group = repeat(["1st", "2nd"], outer = 6)

    plot_write = groupedbar(
        pkgs,
        write_perf,
        group = run_group,
        ylab = "Seconds",
        title = "Disk-format Write performance comparison \n Data: $data_label data \n Size: $(size(df)) filesize:$(round(filesize(joinpath(dirpath, largest_file))/1024^3, digits=1))GB \n Julia $(VERSION)"
    )
    savefig(plot_write, joinpath(outpath, largest_file*"plot_write.png"))

    # plot_write_wo_jlso = groupedbar(
    #     repeat(["JDF.jl", "CSV.jl", "Feather.jl", "JLD2.jl"], inner = 2),
    #     [jdf_write1, jdf_write2, csv_write1, csv_write2, feather_write1, feather_write2, jld2_write1, jld2_write2],
    #     group = repeat(["1st", "2nd"], outer = 4),
    #     ylab = "Seconds",
    #     title = "Disk-format Write performance comparison \n Data: $data_label data \n Size: $(size(df)) filesize:$(round(filesize(joinpath(dirpath, largest_file))/1024^3, digits=1))GB \n Julia $(VERSION)"
    # )
    # savefig(plot_write_wo_jlso, joinpath(outpath, largest_file*"plot_write_wo_jlso.png"))

    plot_read = groupedbar(
        pkgs,
        read_perf,
        group = run_group,
        ylab = "Seconds",
        title = "Disk-format Read performance comparison \n Data: $data_label data \n Size: $(size(df)) filesize:$(round(filesize(joinpath(dirpath, largest_file))/1024^3, digits=1))GB \n Julia $(VERSION)"
    )
    savefig(plot_read, joinpath(outpath, largest_file*"plot_read.png"))

    # plot_read_wo_csv_jlso = groupedbar(
    #     repeat(["JDF.jl", "Feather.jl", "JLD2.jl"], inner = 2),
    #     [jdf_read1, jdf_read2, feather_read1, feather_read2, jld2_read1, jld2_read2],
    #     group = repeat(["1st", "2nd"], outer = 3),
    #     ylab = "Seconds",
    #     title = "Disk-format Read performance comparison \n Data: $data_label data \n Size: $(size(df)) filesize:$(round(filesize(joinpath(dirpath, largest_file))/1024^3, digits=1))GB \n Julia $(VERSION)"
    # )
    # savefig(plot_read_wo_csv_jlso, joinpath(outpath, largest_file*"plot_read_wo_csv_jlso.png"))
    #
    # plot_read_wo_jlso = groupedbar(
    #     repeat(["JDF.jl", "CSV.jl", "Feather.jl", "JLD2.jl"], inner = 2),
    #     [jdf_read1, jdf_read2, csv_read1, csv_read2, feather_read1, feather_read2, jld2_read1, jld2_read2],
    #     group = repeat(["1st", "2nd"], outer = 4),
    #     ylab = "Seconds",
    #     title = "Disk-format Read performance comparison \n Data: $data_label data \n Size: $(size(df)) filesize:$(round(filesize(joinpath(dirpath, largest_file))/1024^3, digits=1))GB \n Julia $(VERSION)"
    # )
    # savefig(plot_read_wo_jlso, joinpath(outpath, largest_file*"plot_read_wo_jlso.png"))

    (write_perf, read_perf, dirpath, outpath, largest_file, df)
end

sum_file_size(dir) = begin
    res = Int[]
    for f in joinpath.(dir, readdir(dir))
        if isfile(f)
            push!(res, filesize(f))
        elseif isdir(f)
            println(f)
            push!(res, sum(filesize.(joinpath.(f,readdir(f)))))
        end
    end
    df = DataFrame(file = readdir(dir), fs = res)

    tmpfn = (x->x[end-3:end])
    df[!,:ext] = tmpfn.(df[!, :file])
    filter!(r -> r.ext != ".png", df)
    df = sort!(df, :ext)[[1;3:4;6:7], :]
    df[:pkg] = ["CSV", "fst", "JDF.jl", "parquet R", "Feather.jl"]
    sort!(df, :pkg)
    df
end


text:
using LinearAlgebra, Plots

Rx(θ::Real) = [1 0 0 0;
        0 cos(θ) -sin(θ) 0;
        0 sin(θ) cos(θ) 0;
        0 0 0 1]

Ry(θ::Real) = [cos(θ) 0 -sin(θ) 0;
        0 1 0 0;
        sin(θ) 0 cos(θ) 0;
        0 0 0 1]

Rz(θ::Real) = [cos(θ) -sin(θ) 0 0;
        sin(θ) cos(θ) 0 0;
        0 0 1 0;
        0 0 0 1]

Tx(d::Real) = [1 0 0 d;
        0 1 0 0;
        0 0 1 0;
        0 0 0 1]

π = pi
d = 2
angToRad(ang::Real) = (π/180)*ang
θ(ang::Real = 150) = angToRad(ang) #70
ω(ang::Real = 160) = angToRad(ang)

B1 = I(4) #nada a fazer

#B2 =  Ry(π)*Tx(d)

B2 = zeros(4,4)
B2[1,1] = -1.0
B2[2,2] = 1.0
B2[3,3] = -1.0
B2[4,4] = 1.0
B2[1,4] = -d

#B3 = Rx(π)*Rz(θ())*Ry(π)*Tx(d)

# tird atom
cθ,sθ = (cos(θ()),sin(θ()))
B3 = zeros(4,4)
B3[1,1] = -cθ
B3[1,2] = -sθ
B3[1,4] = -d*cθ
B3[2,1] = sθ
B3[2,2] = -cθ
B3[2,4] = d*sθ
B3[3,3] = 1.0
B3[4,4] = 1.0

function torsionmatrix(cosθ,sinθ,cosω,sinω,d34,sign::Bool)
	if sign == true

		B=zeros(4,4)
		B[1,1] = -cosθ
		B[1,2] = -sinθ
		B[1,4] = -d34*cosθ
		B[2,1] = sinθ*cosω
		B[2,2] = -cosθ*cosω
		B[2,3] = -sinω
		B[2,4] = d34*sinθ*cosω
		B[3,1] = sinθ*sinω
		B[3,2] = -cosθ*sinω
		B[3,3] = cosω
		B[3,4] = d34*sinθ*sinω
		B[4,4] = 1
	else

		B=zeros(4,4)
		B[1,1] = -cosθ
		B[1,2] = -sinθ
		B[1,4] = -d34*cosθ
		B[2,1] = sinθ*cosω
		B[2,2] = -cosθ*cosω
		B[2,3] = sinω
		B[2,4] = d34*sinθ*cosω
		B[3,1] = -sinθ*sinω
		B[3,2] = cosθ*sinω
		B[3,3] = cosω
		B[3,4] = -d34*sinθ*sinω
		B[4,4] = 1

	end
	return B
end

Bi(ω = ω(), θ = θ(), d = d) = torsionmatrix(cos(angToRad(θ)), sin(angToRad(θ)), cos(angToRad(ω)), sin(angToRad(ω)), d, false) #Bi(ω = ω(), θ = θ(), d = d) = Rx(ω)*Rx(π)*Rz(θ)*Ry(π)*Tx(d)

Bi′(B) = B*inv(Tx(d))

unk = [1 0 0 0;
        0 1 0 0;
        0 0 1 0;
        1 1 1 1]

function sysGenerate(B, α = 1; title = true)
    if title
		return ["x/|x|" "y/|y|" "z/|z|" "i"; B*unk*α]
	else
		return B*unk*α
	end
end

function flatProjection(A, x̂ = [cos(π/4) -cos(π/4)]'*0.5; x :: Bool = true)
	if x
		return [A[2:3,1]+x̂*A[1,1] A[2:3,2]+x̂*A[1,2] A[2:3,3]+x̂*A[1,3] A[2:3,4]+x̂*A[1,4]]
	else
		return [A[1:2,1]+x̂*A[3,1] A[1:2,2]+x̂*A[3,2] A[1:2,3]+x̂*A[3,3] A[1:2,4]+x̂*A[3,4]]
	end
end

function plotFlatProjection(B; sprain=1.01, static = false)
	param = [B[1,4]+B[2,4]im, B[1,4]+B[2,4]im+B[1,1]+B[2,1]im, Inf, B[1,4]+B[2,4]im, B[1,4]+B[2,4]im+B[1,2]+B[2,2]im, Inf, B[1,4]+B[2,4]im, B[1,4]+B[2,4]im+B[1,3]+B[2,3]im]
	if static
		plot!(param, arrow = 2)
	else
		plot(param, arrow = 2)
	end
	annotate!(B[1,4]*sprain, B[2,4]*sprain, text("v", :black, :right, 10))
	annotate!(B[1,4]+B[1,1]*sprain, B[2,4]+B[2,1]*sprain, text("x", :black, :right, 10))
	annotate!(B[1,4]+B[1,2]*sprain, B[2,4]+B[2,2]*sprain, text("y", :black, :right, 10))
	annotate!(B[1,4]+B[1,3]*sprain, B[2,4]+B[2,3]*sprain, text("z", :black, :right, 10))
	# ylims!((-1.5,2)); xlims!((-1.5,1.5))
end

print("DmdgpRealization included\n")

#rotationFactor = Rz(π/2)*Rx(π/2)*Rz(-θ(30))*Rx(π)
#sysGenerate(flatProjection(rotationFactor*B2*B3*Bi(ω(230),θ(120))), 20)
#Bo = B2*B3*Bi(ω(230),θ(120))
#B = Bi(ω(230),θ(120))
#sysGenerate(flatProjection(rotationFactor*Bo*B*B*B), 20)

# a good initial image

# plotFlatProjection(flatProjection(B1, [-cos(π/4) -cos(π/4)]'*0.707))
# plotFlatProjection(flatProjection(B1*Ry(π)*Tx(2.5), [-cos(π/4) -cos(π/4)]'*0.707), static = true)
# plotFlatProjection(flatProjection(B1*B2*Rx(π)*Rz(θ(75))*Ry(π)*Tx(2.5), [-cos(π/4) -cos(π/4)]'*0.707), static = true)
# ylims!((-1,2)); xlims!((-1,4))

max_stars_repo_name: caomem/DmdgpRealization

max_issues_repo_name: caomem/DmdgpRealization

max_forks_repo_name: caomem/DmdgpRealization


text:
y = zeros(100_001,1)
tic()
for i=2:100_000
    y[i+1] = y[i] + randn()
end
toc()


text:
"""
    transform(::BoxCox, 𝐱)

Transform an array using Box-Cox method.  The power parameter `λ` is derived
from maximizing a log-likelihood estimator.
If the array contains any non-positive values then a `DomainError` is thrown.
This can be avoided by providing the shift parameter `α` to make all values
positive.

# Keyword arguments
- `α`: added to all values in 𝐱 before transformation. Default = 0.
- `scaled`: scale transformation results.  Default = false.
"""
function transform(m::BoxCox, 𝐱; kwargs...)
    λ, details = lambda(m, 𝐱; kwargs...)
    @debug "BoxCox: estimated lambda = $λ"
    return transform(m, 𝐱, λ; kwargs...)
end

"""
    transform(::BoxCox, 𝐱, λ; α = 0)

Transform an array using Box-Cox method with the provided power parameter `λ`.
If the array contains any non-positive values then a `DomainError` is thrown.

# Keyword arguments
- `α`: added to all values in 𝐱 before transformation. Default = 0.
- `scaled`: scale transformation results.  Default = false.
"""
function transform(::BoxCox, 𝐱, λ; α = 0, scaled = false, kwargs...)
    if α != 0
        𝐱 .+= α
    end
    any(𝐱 .<= 0) && throw(DomainError(
        "Data must be positive and ideally greater than 1. " *
        "You may specify α argument(shift). "))
    if scaled
        gm = geomean(𝐱)
        return @. λ ≈ 0 ? gm * log(𝐱) : (𝐱 ^ λ - 1) / (λ * gm ^ (λ - 1))
    else
        return @. λ ≈ 0 ? log(𝐱) : (𝐱 ^ λ - 1) / λ
    end
end

"""
    lambda(::BoxCox, 𝐱; interval = (-2.0, 2.0), method = :geomean)

Calculate lambda from an array using a log-likelihood estimator.

# Keyword arguments
- `method`: either :geomean or :normal
- you can also pass any other keyword arguments accepted by `Optim.optimize` function e.g. `abs_tol`

See also: [`log_likelihood`](@ref)
"""
function lambda(m::BoxCox, 𝐱; interval = (-2.0, 2.0), kwargs...)
    i1, i2 = interval
    res = optimize(λ -> -log_likelihood(m, 𝐱, λ; kwargs...), i1, i2)
    return (value=minimizer(res), details=res)
end

"""
    log_likelihood(::BoxCox, 𝐱, λ; method = :geomean)

Return log-likelihood for the given array and lambda.

For method `:geomean`:
```
    -N / 2.0 * log(2 * π * σ² / gm ^ (2 * (λ - 1)) + 1)
```

For method `:normal`:
```
    -N / 2.0 * log(σ²) + (λ - 1) * sum(log.(𝐱))
```
"""
function log_likelihood(m::BoxCox, 𝐱, λ; method = :geomean, kwargs...)
    N = length(𝐱)
    𝐲 = transform(m, float.(𝐱), λ)
    σ² = var(𝐲, corrected = false)
    gm = geomean(𝐱)
    if method == :geomean
        -N / 2.0 * log(2 * π * σ² / gm ^ (2 * (λ - 1)) + 1)
    elseif method == :normal
        -N / 2.0 * log(σ²) + (λ - 1) * sum(log.(𝐱))
    else
        throw(ArgumentError("Incorrect method. Please specify :geomean or :normal."))
    end
end


text:
using CVRPSEP
using Test

include("testdata.jl")

@testset "CVRPSEP.jl" begin
    # Write your tests here.
    demand, capacity, edge_head, edge_tail, edge_x = testdata_provider()
    cut_manager = CutManager()

    idx = findall(v -> v > 0.0, edge_x)
    edge_tail = edge_tail[idx]
    edge_head = edge_head[idx]
    edge_x = edge_x[idx]

    S, RHS = rounded_capacity_cuts!(
        cut_manager,
        demand,
        capacity,
        edge_tail,
        edge_head,
        edge_x,
        integrality_tolerance = 1e-6,
        max_n_cuts = 10
    )

    @show S
    @show RHS

    @test RHS == [13.0, 5.0, 7.0]
end


text:
# Plot Gr4j

function plot_sim{T<:Gr4j}(hydro_out::Array{T,1}; q_obs = [], file_name = [])

  time  = [hydro_out[i].time for i in 1:length(hydro_out)]
  st1   = [hydro_out[i].st[1] for i in 1:length(hydro_out)]
  st2   = [hydro_out[i].st[2] for i in 1:length(hydro_out)]
  q_sim = [hydro_out[i].q_sim for i in 1:length(hydro_out)]

  fig = plt[:figure](figsize = (12,7))

  plt[:style][:use]("ggplot")

  ax = plt[:subplot](211)
  plt[:plot](time, q_sim, linewidth = 1.2, color = "r", label = "Sim")
  plt[:title]("GR4J")
  plt[:ylabel]("Runoff")
  if ~isempty(q_obs)
    nse_res = round(nse(q_sim, q_obs), 2)
    kge_res = round(kge(q_sim, q_obs), 2)
    plt[:plot](time, q_obs, linewidth = 1.2, color = "b", label = "Obs")
    plt[:legend]()
    plt[:title]("GR4J | KGE = $(kge_res) | NSE = $(nse_res)")
  end

  plt[:subplot](212, sharex=ax)
  plt[:plot](time, st1, linewidth = 1.2, color = "k", label = "St1")
  plt[:plot](time, st2, linewidth = 1.2, color = "g", label = "St2")
  plt[:ylabel]("States (mm)")
  plt[:legend]()

  if ~isempty(file_name)
    savefig(file_name)
    close(fig)
  end

end

# Plot Hbv

function plot_sim{T<:Hbv}(hydro_out::Array{T,1}; q_obs = [], file_name = [])

  time  = [hydro_out[i].time for i in 1:length(hydro_out)]
  sm    = [hydro_out[i].sm for i in 1:length(hydro_out)]
  suz   = [hydro_out[i].suz for i in 1:length(hydro_out)]
  slz   = [hydro_out[i].slz for i in 1:length(hydro_out)]
  q_sim = [hydro_out[i].q_sim for i in 1:length(hydro_out)]

  fig = plt[:figure](figsize = (12,7))

  plt[:style][:use]("ggplot")

  ax = plt[:subplot](211)
  plt[:plot](time, q_sim, linewidth = 1.2, color = "r", label = "Sim")
  plt[:title]("HBV")
  plt[:ylabel]("Runoff")
  if ~isempty(q_obs)
    nse_res = round(nse(q_sim, q_obs), 2)
    kge_res = round(kge(q_sim, q_obs), 2)
    plt[:plot](time, q_obs, linewidth = 1.2, color = "b", label = "Obs")
    plt[:legend]()
    plt[:title]("HBV | KGE = $(kge_res) | NSE = $(nse_res)")
  end

  plt[:subplot](212, sharex = ax)
  plt[:plot](time, sm, linewidth = 1.2, color = "k", label = "SM")
  plt[:plot](time, suz, linewidth = 1.2, color = "b", label = "SUZ")
  plt[:plot](time, slz, linewidth = 1.2, color = "g", label = "SLZ")
  plt[:ylabel]("States (mm)")
  plt[:legend]()

  if ~isempty(file_name)
    savefig(file_name)
    close(fig)
  end

end

# Plot TinBasic

function plot_sim{T<:TinBasic}(snow_out::Array{T,1}; file_name = [])

  time = [snow_out[i].time for i in 1:length(snow_out)]
  swe  = [mean(snow_out[i].swe) for i in 1:length(snow_out)]

  fig = plt[:figure](figsize = (12,7))

  plt[:style][:use]("ggplot")

  plt[:plot](time, swe, linewidth = 1.2, color = "b")
  plt[:title]("TinBasic")
  plt[:ylabel]("SWE (mm)")

  if ~isempty(file_name)
    savefig(file_name)
    close(fig)
  end

end

# Plot TinStandard

function plot_sim{T<:TinStandard}(snow_out::Array{T,1}; file_name = [])

  time = [snow_out[i].time for i in 1:length(snow_out)]
  swe  = [mean(snow_out[i].swe) for i in 1:length(snow_out)]
  lw   = [mean(snow_out[i].lw) for i in 1:length(snow_out)]

  fig = plt[:figure](figsize = (12,7))

  plt[:style][:use]("ggplot")

  ax = plt[:subplot](211)
  plt[:plot](time, swe, linewidth = 1.2, color = "b")
  plt[:ylabel]("SWE (mm)")
  plt[:title]("TinStandard")

  plt[:subplot](212, sharex = ax)
  plt[:plot](time, lw, linewidth = 1.2, color = "r")
  plt[:ylabel]("LW (mm)")

  if ~isempty(file_name)
    savefig(file_name)
    close(fig)
  end

end

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
import Distributions: truncated

function validate_bounds(lower, upper, uv, constraint)
    if lower >= upper
        msg = "lower_bound < upper_bound required, got $lower <= $upper"
        msg2 = "cannot truncate $uv with $constraint\n"
        throw(ArgumentError(msg2*msg))
    end
end
################################################################
# Truncating uncertain values based on theoretical distributions
# Operating on the union of both TheoreticalFittedUncertainScalar
# and TheoreticalDistributionScalarValue as the type of the
# uncertain value is ok, because Truncated is defined for
# FittedDistribution, which is the .distribution fild for
# fitted scalars.
################################################################
"""
    truncate(uv::TheoreticalDistributionScalarValue, constraint::SamplingConstraint)

Truncate an uncertain value `uv` represented by a theoretical distribution
according to the sampling `constraint`.
"""
Base.truncate(uv::TheoreticalDistributionScalarValue, constraint::SamplingConstraint)

"""
    truncate(uv::TheoreticalDistributionScalarValue,
        constraint::NoConstraint)

Truncate the theoretical distribution furnishing `uv` using a
`NoConstraint` sampling constraint.
"""
function Base.truncate(uv::TheoreticalDistributionScalarValue,
        constraint::NoConstraint)
    s = support(uv.distribution)
    lower_bound, upper_bound = s.lb, s.ub

    return truncated(uv.distribution, lower_bound, upper_bound)
end

"""
    truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateQuantiles)

Truncate the theoretical distribution furnishing `uv` using a
`TruncateQuantiles` sampling constraint.
"""
function Base.truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateQuantiles)
    lower_bound = quantile(uv.distribution, constraint.lower_quantile)
    upper_bound = quantile(uv.distribution, constraint.upper_quantile)

    validate_bounds(lower_bound, upper_bound, uv, constraint)

    truncated(uv.distribution, lower_bound, upper_bound)
end

"""
    truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateLowerQuantile)

Truncate the theoretical distribution furnishing `uv` using a
`TruncateLowerQuantile` sampling constraint.
"""
function Base.truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateLowerQuantile)
    lower_bound = quantile(uv.distribution, constraint.lower_quantile)
    upper_bound = maximum(uv)

    validate_bounds(lower_bound, upper_bound, uv, constraint)

    truncated(uv.distribution, lower_bound, upper_bound)
end

"""
    truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateUpperQuantile)

Truncate the theoretical distribution furnishing `uv` using a
`TruncateUpperQuantile` sampling constraint.
"""
function Base.truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateUpperQuantile)
    lower_bound = minimum(uv)
    upper_bound = quantile(uv.distribution, constraint.upper_quantile)

    validate_bounds(lower_bound, upper_bound, uv, constraint)

    truncated(uv.distribution, lower_bound, upper_bound)
end

"""
    truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateMinumum)

Truncate the theoretical distribution furnishing `uv` using a
`TruncateMinimum` sampling constraint.
"""
function Base.truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateMinimum)

    lower_bound = constraint.min
    upper_bound = maximum(uv)

    validate_bounds(lower_bound, upper_bound, uv, constraint)

    truncated(uv.distribution, lower_bound, upper_bound)
end

"""
    truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateMaximum)

Truncate the theoretical distribution furnishing `uv` using a
`TruncateMaximum` sampling constraint.
"""
function Base.truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateMaximum)

    lower_bound = minimum(uv)
    upper_bound = constraint.max

    validate_bounds(lower_bound, upper_bound, uv, constraint)

    truncated(uv.distribution, lower_bound, upper_bound)
end

"""
    truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateRange)

Truncate the theoretical distribution furnishing `uv` using a
`TruncateRange` sampling constraint.
"""
function Base.truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateRange)

    lower_bound = constraint.min
    upper_bound = constraint.max

    validate_bounds(lower_bound, upper_bound, uv, constraint)

    truncated(uv.distribution, lower_bound, upper_bound)
end

"""
    truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateStd, n::Int = 10000)

Truncate the theoretical distribution furnishing `uv` using a
`TruncateStd` sampling constraint.

This functions needs to compute the mean and standard deviation
of a truncated distribution, so takes an extra optional argument `n_draws` to allow
this.
"""
function Base.truncate(uv::TheoreticalDistributionScalarValue,
        constraint::TruncateStd)

    m = mean(uv.distribution)
    s = std(uv.distribution)
    lower_bound = m - s
    upper_bound = m + s

    validate_bounds(lower_bound, upper_bound, uv, constraint)

    truncated(uv.distribution, lower_bound, upper_bound)
end

export truncate


text:
# Multithreaded multiplication using Polyester.jl @batch
using SparseMatricesCSR: nzrange

# sparse mat-vec multiplication
"""
    bmul!(y::AbstractVector, A::SparseMatrixCSR, x::AbstractVector, alpha::Number, beta::Number)
    bmul!(y::AbstractVector, A::SparseMatrixCSR, x::AbstractVector)

Evaluates `y = alpha*A*x + beta*y` (`y = A*x`)

In-place multithreaded version of sparse csr matrix - vector multiplication, using the threading provided by Polyester.jl
"""
function bmul!(y::AbstractVector, A::SparseMatrixCSR, x::AbstractVector, alpha::Number, beta::Number)

    A.n == size(x, 1) || throw(DimensionMismatch())
    A.m == size(y, 1) || throw(DimensionMismatch())

    o = getoffset(A)

    @batch minbatch = size(y, 1) ÷ matmul_num_threads[] for row in 1:size(y, 1)
        @inbounds begin
            accu = zero(eltype(y))
            for nz in nzrange(A, row)
                col = A.colval[nz] + o
                accu += A.nzval[nz]*x[col]
            end
            y[row] = alpha*accu + beta*y[row]

        end
    end

    return y

end

function bmul!(y::AbstractVector, A::SparseMatrixCSR, x::AbstractVector)

    return bmul!(y, A, x, true, false)

end

"""
    bmul(A::SparseMatrixCSR, x::AbstractVector)

Evaluates `A*x`.

Multithreaded version of sparse csr matrix - vector multiplication, using the threading provided by Polyester.jl
"""
function bmul(A::SparseMatrixCSR, x::AbstractVector)

    y = similar(x, promote_type(eltype(A), eltype(x)), size(A, 1))

    return bmul!(y, A, x, true, false)

end


text:
# Import modules
include("GModel.jl") # Contains Lorenz 96 source code

# Import modules
using Distributions  # probability distributions and associated functions
using LinearAlgebra
using StatsPlots
using Plots
using Random
using JLD2

# CES
using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.Observations
using EnsembleKalmanProcesses.ParameterDistributionStorage
using EnsembleKalmanProcesses.DataStorage

rng_seed = 4137
Random.seed!(rng_seed)

# Output figure save directory
homedir = pwd()
println(homedir)
figure_save_directory = homedir*"/output/"
data_save_directory = homedir*"/output/"
if ~isdir(figure_save_directory)
    mkdir(figure_save_directory)
end
if ~isdir(data_save_directory)
    mkdir(data_save_directory)
end

# Governing settings
# Characteristic time scale
τc = 5. # days, prescribed by the L96 problem
# Stationary or transient dynamics
dynamics = 2 # Transient is 2
# Statistics integration length
# This has to be less than 360 and 360 must be divisible by Ts_days
Ts_days = 90. # Integration length in days
# Stats type, which statistics to construct from the L96 system
# 4 is a linear fit over a batch of length Ts_days
# 5 is the mean over a batch of length Ts_days
stats_type = 5

###
###  Define the (true) parameters
###
# Define the parameters that we want to learn
F_true = 8. # Mean F
A_true = 2.5 # Transient F amplitude
ω_true = 2. * π / (360. / τc) # Frequency of the transient F

if dynamics==2
    params_true = [F_true, A_true]
    param_names = ["F", "A"]
else
    params_true = [F_true]
    param_names = ["F"]
end
n_param = length(param_names)
params_true = reshape(params_true, (n_param,1))

println(n_param)
println(params_true)

###
###  Define the parameter priors
###
# Lognormal prior or normal prior?
log_normal = false # THIS ISN't CURRENTLY IMPLEMENTED

function logmean_and_logstd(μ, σ)
    σ_log = sqrt(log(1.0 + σ^2/μ^2))
    μ_log = log(μ / (sqrt(1.0 + σ^2/μ^2)))
    return μ_log, σ_log
end
#logmean_F, logstd_F = logmean_and_logstd(F_true, 5)
#logmean_A, logstd_A = logmean_and_logstd(A_true, 0.2*A_true)

if dynamics == 2
	prior_means = [F_true+1.0, A_true+0.5]
	prior_stds = [2.0, 0.5*A_true]
	d1 = Parameterized(Normal(prior_means[1], prior_stds[1]))
	d2 = Parameterized(Normal(prior_means[2], prior_stds[2]))
	prior_distns = [d1, d2]
	c1 = no_constraint()
	c2 = no_constraint()
	constraints = [[c1], [c2]]
	prior_names = param_names
	priors = ParameterDistribution(prior_distns, constraints, prior_names)
else
	prior_distns = [Parameterized(Normal(F_true, 1))]
	constraints = [[no_constraint()]]
	prior_names = ["F"]
	priors = ParameterDistribution(prior_distns, constraints, prior_names)
end

###
###  Define the data from which we want to learn the parameters
###
data_names = ["y0", "y1"]

###
###  L96 model settings
###

# Lorenz 96 model parameters
# Behavior changes depending on the size of F, N
N = 36
dt = 1/64.
# Start of integration
t_start = 800.
# Data collection length
if dynamics==1
    T = 2.
else
    T = 360. / τc
end
# Batch length
Ts = 5. / τc # Nondimensionalize by L96 timescale
# Integration length
Tfit = Ts_days / τc
# Initial perturbation
Fp = rand(Normal(0.0, 0.01), N);
kmax = 1
# Prescribe variance or use a number of forward passes to define true interval variability
var_prescribe = false
# Use CES to learn ω?
ω_fixed = true

# Settings
# Constructs an LSettings structure, see GModel.jl for the descriptions
lorenz_settings = GModel.LSettings(dynamics, stats_type, t_start,
				   T, Ts, Tfit, Fp, N, dt, t_start+T, kmax,
				   ω_fixed, ω_true);
lorenz_params = GModel.LParams(F_true, ω_true, A_true)

###
###  Generate (artificial) truth samples
###  Note: The observables y are related to the parameters θ by:
###        y = G(θ) + η
###

# Lorenz forward
# Input: params: [N_params, N_ens]
# Output: gt: [N_data, N_ens]
# Dropdims of the output since the forward model is only being run with N_ens=1
# corresponding to the truth construction
gt = dropdims(GModel.run_G_ensemble(params_true, lorenz_settings), dims=2)

# Compute internal variability covariance
if var_prescribe==true
    n_samples = 100
    yt = zeros(length(gt),n_samples)
    noise_level = 0.05
    Γy = noise_level * convert(Array, Diagonal(gt))
    μ = zeros(length(gt))
    # Add noise
    for i in 1:n_samples
        yt[:,i] = gt .+ rand(MvNormal(μ, Γy))
    end
else
    println("Using truth values to compute covariance")
    n_samples = 20
    yt = zeros(length(gt), n_samples)
    for i in 1:n_samples
	    lorenz_settings_local = GModel.LSettings(dynamics, stats_type, t_start+T*(i-1),
					      T, Ts, Tfit, Fp, N, dt,
					      t_start+T*(i-1)+T, kmax, ω_fixed, ω_true);
	    yt[:, i] = GModel.run_G_ensemble(params_true, lorenz_settings_local)
    end
    # Variance of truth data
    #Γy = convert(Array, Diagonal(dropdims(mean((yt.-mean(yt,dims=1)).^2,dims=1),dims=1)))
    # Covariance of truth data
    Γy = cov(yt, dims=2)

    println(Γy)
end

# Construct observation object
truth = Observations.Obs(yt, Γy, data_names)
truth_sample = truth.mean
###
###  Calibrate: Ensemble Kalman Inversion
###

# L96 settings for the forward model in the EKP
# Here, the forward model for the EKP settings can be set distinctly from the truth runs
lorenz_settings_G = lorenz_settings; # initialize to truth settings

# EKP parameters
log_transform(a::AbstractArray) = log.(a)
exp_transform(a::AbstractArray) = exp.(a)

N_ens = 20 # number of ensemble members
N_iter = 5 # number of EKI iterations
# initial parameters: N_params x N_ens
initial_params = construct_initial_ensemble(priors, N_ens; rng_seed=rng_seed)

ekiobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_params,
						       truth_sample,
						       truth.obs_noise_cov,
						       Inversion())

# EKI iterations
println("EKP inversion error:")
err = zeros(N_iter)
for i in 1:N_iter
    if log_normal==false
        params_i = get_u_final(ekiobj)
    else
        params_i = exp_transform(get_u_final(ekiobj))
    end
    g_ens = GModel.run_G_ensemble(params_i, lorenz_settings_G)
    EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens)
    err[i] = get_error(ekiobj)[end] #mean((params_true - mean(params_i,dims=2)).^2)
    println("Iteration: "*string(i)*", Error: "*string(err[i]))
end

# EKI results: Has the ensemble collapsed toward the truth?
println("True parameters: ")
println(params_true)

println("\nEKI results:")
if log_normal==false
    println(mean(get_u_final(ekiobj), dims=2))
else
    println(mean(exp_transform(get_u_final(ekiobj)), dims=2))
end

u_stored= get_u(ekiobj,return_array=false)
g_stored= get_g(ekiobj,return_array=false)

@save data_save_directory*"parameter_storage.jld2" u_stored
@save data_save_directory*"data_storage.jld2" g_stored

#plots
gr(size=(600,600))
u_init = get_u_prior(ekiobj)
for i in 1:N_iter
    u_i = get_u(ekiobj,i)
    p = plot(u_i[1,:], u_i[2,:], seriestype=:scatter, xlims = extrema(u_init[1,:]), ylims = extrema(u_init[2,:]))
    plot!([params_true[1]], xaxis="u1", yaxis="u2", seriestype="vline",
        linestyle=:dash, linecolor=:red, label = false,
        title = "EKI iteration = " * string(i)
        )
    plot!([params_true[2]], seriestype="hline", linestyle=:dash, linecolor=:red, label = "optimum")
    #display(p)
    figpath = joinpath(figure_save_directory, "posterior_EKP_it_$(i).png")
    savefig(figpath)
    #linkfig(figpath)
    sleep(0.5)
end


text:
int(x) = (Int ∘ round)(x)
fint(x) = (Int ∘ floor)(x)
cint(x) = (Int ∘ ceil)(x)

"""
  or(x, y)

Return `y` if `x` is `Nothing`
"""
or(x, y) = isnothing(x) ? y : x

"""
  loop_last(v)

  Returns an iterable yielding tuples (is_last, value).
"""
function loop_last(v)
    is_last = [i == length(v) for i in 1:length(v)]
    return zip(is_last, v)
end

"""
  get_lr_widths(width::Int)

To split something with `width` in 2, get the lengths
of the left/right widths.

When width is even that's easy, when it's odd we need to
be careful.
"""
function get_lr_widths(width::Int)::Tuple{Int, Int}
  iseven(width) && return (int(width/2), int(width/2))
  return (fint(width/2), cint(width/2))
end


text:
#!/usr/bin/env julia

using Seq

natural_numbers = cons(1, @lazyseq map(n -> n + 1, natural_numbers))

for k in take(10, natural_numbers)
    println(k)
end


text:
using Distributed
addprocs(8)

@everywhere begin
    include("MCTS.jl")
    include("./utils.jl")
    using ProgressMeter
end

@everywhere function playOnce(curr_nn::NeuralNet, firstHand::Bool=true)
    curr_nn = curr_nn[1]
    mcts = MCTS()
    game = Game()
    randomPlace(game, 2)
    if firstHand
        while isnothing(game.winner)
            move, _ = predict(game, mcts, curr_nn)
            place(game, move)
            randomPlace(game, 1)
        end
    else
        while isnothing(game.winner)
            randomPlace(game, 1)
            move, _ = predict(game, mcts, curr_nn)
            place(game, move)
        end
    end
    save_log("evaluation", ("first hand ? $firstHand", game))
    return firstHand ? game.winner == -1 : game.winner == 1
end

@everywhere function playOnce(players::Array{NeuralNet}, firstHand::Bool=true)
    playerX, playerO =  firstHand ? players : (players[2], players[1])
    mctsX = MCTS()
    mctsO = MCTS()
    game = Game()
    randomPlace(game, 2)
    while isnothing(game.winner)
        move, _ = predict(game, mctsX, playerX)
        place(game, move)
        move, _ = predict(game, mctsO, playerO)
        place(game, move)
    end
    save_log("evaluation", ("first hand ? $firstHand", game))
    return firstHand ? game.winner == -1 : game.winner == 1
end

function evaluate(players::Union{Array{NeuralNet}, NeuralNet}, num_games::Int = 20)
    desc = players isa NeuralNet ? "random evaluating: " : "normal evaluating"
    p = Progress(num_games, 0.5, desc, 50)
    win_counts = progress_pmap(1:num_games, progress=p) do i
        playOnce(players, isodd(i))
    end
    players isa NeuralNet ? checkQuality(win_counts, players) : checkQuality(win_counts, players[1])
    nothing
end


text:
# ~~ Type to set\get data along a band
struct Band
    i::Int
end

show(io::IO, r::Band) = print(io, "Band(", r.i, ")")

"""
    band(i)

Represents the `i`-th band of a banded matrix.

```jldoctest
julia> using BandedMatrices

julia> A = BandedMatrix(Ones(5,5),(1,1))
5×5 BandedMatrix{Float64,Array{Float64,2}}:
 1.0  1.0   ⋅    ⋅    ⋅
 1.0  1.0  1.0   ⋅    ⋅
  ⋅   1.0  1.0  1.0   ⋅
  ⋅    ⋅   1.0  1.0  1.0
  ⋅    ⋅    ⋅   1.0  1.0

julia> A[band(1)]
4-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0
```
"""
band(i::Int) = Band(i)

"""
    BandRange

Represents the entries in a row/column inside the bands.

```jldoctest
julia> using BandedMatrices

julia> A = BandedMatrix(Ones(5,5),(1,1))
5×5 BandedMatrix{Float64,Array{Float64,2}}:
 1.0  1.0   ⋅    ⋅    ⋅
 1.0  1.0  1.0   ⋅    ⋅
  ⋅   1.0  1.0  1.0   ⋅
  ⋅    ⋅   1.0  1.0  1.0
  ⋅    ⋅    ⋅   1.0  1.0

julia> A[2,BandRange]
3-element Array{Float64,1}:
 1.0
 1.0
 1.0
```
"""
struct BandRangeType end
const BandRange = BandRangeType()

to_indices(A::AbstractArray, (_, j)::Tuple{BandRangeType,Integer}) = (colrange(A, j), j)
to_indices(A::AbstractArray, (k, _)::Tuple{Integer,BandRangeType}) = (k, rowrange(A, k))

# ~~ Out of band error
struct BandError <: Exception
    A::AbstractMatrix
    i::Int
end

# shorthand to specify k and j without calculating band
BandError(A::AbstractMatrix, (k, j)::Tuple{Int,Int}) = BandError(A, j-k)
BandError(A::AbstractMatrix) = BandError(A, max(size(A)...)-1)

function showerror(io::IO, e::BandError)
    A, i = e.A, e.i
    print(io, "attempt to access $(typeof(A)) with bandwidths " *
              "($(bandwidth(A, 1)), $(bandwidth(A, 2))) at band $i")
end

# length of diagonal
@inline diaglength(A::AbstractMatrix, b::Band) = diaglength(A, b.i)
@inline function diaglength(A::AbstractMatrix, i::Integer)
    max(min(size(A, 2), size(A, 1)+i) - max(0, i), 0)
end

# check indices fall in the band
checkband(A::AbstractMatrix, i::Integer) =
    (-bandwidth(A, 1) ≤ i ≤ bandwidth(A, 2) || throw(BandError(A, i)))

checkband(A::AbstractMatrix, b::Band) = checkband(A, b.i)

checkband(A::AbstractMatrix, k::Integer, j::Integer) = checkband(A, j-k)

checkband(A::AbstractMatrix, kr::AbstractRange, j::Integer) =
    (checkband(A, first(kr), j); checkband(A,  last(kr), j))

checkband(A::AbstractMatrix, k::Integer, jr::AbstractRange) =
    (checkband(A, k, first(jr)); checkband(A, k,  last(jr)))

checkband(A::AbstractMatrix, kr::AbstractRange, jr::AbstractRange) =
    (checkband(A, kr, first(jr)); checkband(A, kr,  last(jr)))

# checks if the bands match A
function checkbandmatch(A::AbstractMatrix{T}, V::AbstractVector, ::Colon, j::Integer) where {T}
    for k = 1:colstart(A,j)-1
        if V[k] ≠ zero(T)
            throw(BandError(A, (k,j)))
        end
    end
    for k = colstop(A,j)+1:size(A,1)
        if V[k] ≠ zero(T)
            throw(BandError(A, (k,j)))
        end
    end
end

function checkbandmatch(A::AbstractMatrix{T}, V::AbstractVector, kr::AbstractRange, j::Integer) where {T}
    a = colstart(A, j)
    b = colstop(A, j)
    i = 0
    for v in V
        k = kr[i+=1]
        if (k < a || k > b) && v ≠ zero(T)
            throw(BandError(A, (k,j)))
        end
    end
end

function checkbandmatch(A::AbstractMatrix{T}, V::AbstractVector, k::Integer, ::Colon) where {T}
    for j = 1:rowstart(A,k)-1
        if V[j] ≠ zero(T)
            throw(BandError(A, (k,j)))
        end
    end
    for j = rowstop(A,j)+1:size(A,2)
        if V[j] ≠ zero(T)
            throw(BandError(A, (k,j)))
        end
    end
end

function checkbandmatch(A::AbstractMatrix{T}, V::AbstractVector, k::Integer, jr::AbstractRange) where {T}
    a = rowstart(A, k)
    b = rowstop(A, k)
    i = 0
    for v in V
        j = jr[i+=1]
        if (j < a || j > b) && v ≠ zero(T)
            throw(BandError(A, (k,j)))
        end
    end
end

function checkbandmatch(A::AbstractMatrix{T}, V::AbstractMatrix, kr::AbstractRange, jr::AbstractRange) where {T}
    u, l = A.u, A.l
    jj = 1
    for j in jr
        kk = 1
        for k in kr
            if !(-l ≤ j - k ≤ u) && V[kk, jj] ≠ zero(T)
                # we index V manually in column-major order
                throw(BandError(A, (k,j)))
            end
            kk += 1
        end
        jj += 1
    end
end

checkbandmatch(A::AbstractMatrix, V::AbstractMatrix, ::Colon, ::Colon) =
    checkbandmatch(A, V, 1:size(A,1), 1:size(A,2))

"""
    BandSlice(band, indices)

Represent an AbstractUnitRange of indices corresponding to a band.

Upon calling `to_indices()`, Bands are converted to BandSlice objects to represent
the indices over which the Band spans.

This mimics the relationship between `Colon` and `Base.Slice`.
"""
struct BandSlice <: OrdinalRange{Int,Int}
    band::Band
    indices::StepRange{Int,Int}
end

for f in (:indices, :unsafe_indices, :axes1, :first, :last, :size, :length,
          :unsafe_length, :start, :step)
    @eval $f(S::BandSlice) = $f(S.indices)
end

getindex(S::BandSlice, i::Int) = getindex(S.indices, i)
show(io::IO, r::BandSlice) = print(io, "BandSlice(", r.band, ",", r.indices, ")")
next(S::BandSlice, s) = next(S.indices, s)
done(S::BandSlice, s) = done(S.indices, s)

to_index(::Band) = throw(ArgumentError("Block must be converted by to_indices(...)"))

@inline to_indices(A, I::Tuple{Band}) = (BandSlice(I[1], diagind(A, I[1].i)),)


text:
function get_order(obj::EmbeddedRungeKuttaIntegrator)
    return jcall(obj, "getOrder", jint, ())
end


text:
immutable PINningSynapseParameter
end

@withkw type PINningSynapse
    param::PINningSynapseParameter = PINningSynapseParameter()
    W::Matrix{Float}  # synaptic weight
    rI::Vector{Float} # postsynaptic rate
    rJ::Vector{Float} # presynaptic rate
    g::Vector{Float}  # postsynaptic conductance
    P::Matrix{Float}  # <rᵢrⱼ>⁻¹
    q::Vector{Float}  # P * r
    f::Vector{Float}  # postsynaptic traget
    records::Dict = Dict()
end

function PINningSynapse(pre, post; σ = 1.5, p = 0.0, α = 1)
    rI, rJ, g = post.r, pre.r, post.g
    W = σ * 1 / √pre.N * randn(post.N, pre.N) # normalized recurrent weight
    P = α * eye(post.N) # initial inverse of C = <rr'>
    f, q = zeros(post.N), zeros(post.N)
    PINningSynapse(;@symdict(W, rI, rJ, g, P, q, f)...)
end

@replace function forward!(c::PINningSynapse, param::PINningSynapseParameter)
    BLAS.A_mul_B!(q, P, rJ)
    BLAS.A_mul_B!(g, W, rJ)
end

@replace function plasticity!(c::PINningSynapse, param::PINningSynapseParameter, dt::Float, t::Float)
    C = 1 / (1 + dot(q, rI))
    BLAS.ger!(C, f - g, q, W)
    BLAS.ger!(-C, q, q, P)
end


text:
input = joinpath(@__DIR__, "input")
eqs = readlines(input)

EQ = r"(\d+)|([+])|([*])|([(])|([)])"

function parse_eqs(eqs)
    acc = []
    for eq in eqs
        matches = eachmatch(EQ, eq)
        r = map(m -> only(filter(!isnothing, m.captures)), matches)
        push!(acc, r)
    end
    return acc
end

function convert(eqs, prec)::Vector{Vector{String}}
    acc = []
    for eq in eqs
        res = postfix(eq, prec)
        push!(acc, res)
    end
    return acc
end

function postfix(eq, prec)
    postfix = []
    operators = []
    for symbol in eq
        if symbol in ["+", "*"]
            if !isempty(operators) && last(operators) in ["+", "*"] && prec[symbol] <= prec[last(operators)]
                push!(postfix, pop!(operators))
                push!(operators, symbol)
            else
                push!(operators, symbol)
            end
        elseif symbol == "("
            push!(operators, symbol)
        elseif symbol == ")"
            while (op = pop!(operators); op != "(")
                push!(postfix, op)
            end
        else
            push!(postfix, symbol)
        end
    end
    while !isempty(operators)
        op = pop!(operators)
        push!(postfix, op)
    end
    return postfix
end

function evaluate(eqs::Vector{Vector{String}})
    return map(eq -> evaluate(eq), eqs)
end

function evaluate(eq::Vector{String})
    stack = []
    for symbol in eq
        if symbol == "+"
            left = pop!(stack)
            right = pop!(stack)
            push!(stack, left + right)
        elseif symbol == "*"
            left = pop!(stack)
            right = pop!(stack)
            push!(stack, left * right)
        else
            push!(stack, parse(Int, symbol))
        end
    end
    return only(stack)
end

eqs = filter(!isempty, eqs)
eqs = parse_eqs(eqs)

p1 = sum(evaluate(convert(eqs, Dict("+" => 1, "*" => 1))))
p2 = sum(evaluate(convert(eqs, Dict("+" => 2, "*" => 1))))

println("-----------------------------------------------------------------------")
println("operation order -- part one :: $p1")
println("operation order -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 21993583522852)
@assert(p2 == 122438593522757)


text:
# ------------------------------------
# time integration with explicit methods
# ------------------------------------
function time_integration_explicit(dt, Qcon_hat, RHS, cellxmax, cellymax, nval)
    Threads.@threads for l in 1:nval
        for j in 2:cellymax-1
            for i in 2:cellxmax-1
                Qcon_hat[i,j,l] = Qcon_hat[i,j,l] + dt*RHS[i,j,l]
            end
        end
    end
    return Qcon_hat
end

# ------------------------------------
# Calculation of the advection Jacobian matrix by one-wave approximation for lusgs
# ------------------------------------
function one_wave(A_adv_hat_m, A_adv_hat_p, B_adv_hat_m, B_adv_hat_p, A_beta_shig, B_beta_shig, I,
                    Qbase, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, volume, nval)
    beta = 1.1

    for j in 2:cellymax
        for i in 2:cellxmax
            for k in 1:2 #A,B
                jacob_temp = zeros(nval,nval)
                if k == 1
                    kx_av = 0.5*(vecAx[i,j,1]+vecAx[i+1,j,1]) / volume[i,j]
                    ky_av = 0.5*(vecAx[i,j,2]+vecAx[i+1,j,2]) / volume[i,j]
                elseif k == 2
                    kx_av = 0.5*(vecAy[i,j,1]+vecAy[i,j+1,1]) / volume[i,j]
                    ky_av = 0.5*(vecAy[i,j,2]+vecAy[i,j+1,2]) / volume[i,j]
                end

                rho = Qbase[i,j,1]
                u = Qbase[i,j,2]
                v = Qbase[i,j,3]
                e = Qcon[i,j,4]
                p = Qbase[i,j,4]

                g = specific_heat_ratio

                Z = kx_av*u+ky_av*v
                q2   = 0.5*(u^2+v^2)
                b1c2 = 0.5*q2*(g-1)
                gebyrho = g*e/rho

                # advection Jacobian matrix in the general coordinate system
                jacob_temp[1,1] = 0.0
                jacob_temp[1,2] = kx_av
                jacob_temp[1,3] = ky_av
                jacob_temp[1,4] = 0.0

                jacob_temp[2,1] = -u*Z + kx_av*b1c2
                jacob_temp[2,2] = Z - (g-2)*kx_av*u
                jacob_temp[2,3] = ky_av*u - kx_av*(g-1)*v
                jacob_temp[2,4] = (g-1)*kx_av

                jacob_temp[3,1] = -v*Z + ky_av*b1c2
                jacob_temp[3,2] = kx_av*v - ky_av*(g-1)*u
                jacob_temp[3,3] = Z - ky_av*(g-2)*v
                jacob_temp[3,4] = (g-1)*ky_av

                jacob_temp[4,1] = Z*(-gebyrho + 2*b1c2)
                jacob_temp[4,2] = kx_av*(gebyrho-b1c2) - (g-1)*u*Z
                jacob_temp[4,3] = ky_av*(gebyrho-b1c2) - (g-1)*v*Z
                jacob_temp[4,4] = g*Z

                c = (g*rho/p)^0.5
                shigma = abs(Z) + c*(kx_av^2+ky_av^2)^0.5

                for l in 1:nval
                    I[l,l] = beta * shigma
                end

                if k == 1
                    for m in 1:nval
                        for l in 1:nval
                            A_adv_hat_p[i,j,l,m] = 0.5*(jacob_temp[l,m] + I[l,m])
                            A_adv_hat_m[i,j,l,m] = 0.5*(jacob_temp[l,m] - I[l,m])
                        end
                    end
                    A_beta_shig[i,j] = beta * shigma
                elseif k ==2
                    for m in 1:nval
                        for l in 1:nval
                            B_adv_hat_p[i,j,l,m] = 0.5*(jacob_temp[l,m] + I[l,m])
                            B_adv_hat_m[i,j,l,m] = 0.5*(jacob_temp[l,m] - I[l,m])
                        end
                    end
                    B_beta_shig[i,j] = beta * shigma
                end

            end
        end
    end

    # reset
    for l in 1:nval
        I[l,l] = 1.0
    end

    return A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, A_beta_shig, B_beta_shig
end

# ------------------------------------
# Calculation of the viscosity Jacobian matrix by center differential for lusgs
# ------------------------------------
function central_diff_jacobian(jalphaP, jbetaP, Qbase, Qcon, cellxmax, cellymax, mu, lambda,
                            vecAx, vecAy, specific_heat_ratio, volume, nval)

    # substitution to Corner cell
    Qbase[1,1,1] = 1.0
    Qbase[cellxmax,1,1] = 1.0
    Qbase[1,cellymax,1] = 1.0
    Qbase[cellxmax,cellymax,1] = 1.0
    volume[1,1,1] = 1.0
    volume[cellxmax,1,1] = 1.0
    volume[1,cellymax,1] = 1.0
    volume[cellxmax,cellymax,1] = 1.0

    for j in 1:cellymax
        for i in 1:cellxmax
            rho = Qbase[i,j,1]
            xi_x_av = 0.5*(vecAx[i,j,1] + vecAx[i+1,j,1]) / volume[i,j]
            xi_y_av = 0.5*(vecAx[i,j,2] + vecAx[i+1,j,2]) / volume[i,j]
            alpha = (xi_x_av^2+xi_y_av^2)^0.5
            jalphaP[i,j] = alpha / volume[i,j] * (2*mu[i,j]/rho)

            eta_x_av = 0.5*(vecAy[i,j,1] + vecAy[i,j+1,1]) / volume[i,j]
            eta_y_av = 0.5*(vecAy[i,j,2] + vecAy[i,j+1,2]) / volume[i,j]
            beta = (xi_x_av^2+xi_y_av^2)^0.5
            jbetaP[i,j] = beta / volume[i,j] * (2*mu[i,j]/rho)
        end
    end
    return jalphaP, jbetaP
end

# ------------------------------------
# lusgs
# ------------------------------------
function lusgs(D, Lx, Ly, Ux, Uy, LdQ, UdQ, RHS_temp, I, dt, dtau, Qcon_hat, Qconn_hat, delta_Q,
                A_adv_hat_p,  A_adv_hat_m,  B_adv_hat_p,  B_adv_hat_m,  A_beta_shig,  B_beta_shig,
                jalphaP,  jbetaP, RHS, cellxmax, cellymax, volume, nval)

    # calculate L and U
    for m in 1:nval
        for l in 1:nval
            for j in 1:cellymax
                for i in 1:cellxmax
                    Lx[i,j,l,m] = dt*(A_adv_hat_p[i,j,l,m] + jalphaP[i,j]*I[l,m])
                    Ly[i,j,l,m] = dt*(B_adv_hat_p[i,j,l,m] + jbetaP[i,j]*I[l,m])
                    Ux[i,j,l,m] = dt*(A_adv_hat_m[i,j,l,m] - jalphaP[i,j]*I[l,m])
                    Uy[i,j,l,m] = dt*(B_adv_hat_m[i,j,l,m] - jbetaP[i,j]*I[l,m])
                end
            end
        end
    end

    for j in 2:cellymax-1
        for i in 2:cellxmax-1
            for l in 1:nval
                for m in 1:nval
                    LdQ[i,j,l] += Lx[i-1,j,l,m]*delta_Q[i-1,j,m] + Ly[i,j-1,l,m]*delta_Q[i,j-1,m]
                    UdQ[i,j,l] += Ux[i+1,j,l,m]*delta_Q[i+1,j,m] + Uy[i,j+1,l,m]*delta_Q[i,j+1,m]
                end
            end
        end
    end

    # diagonal
    for j in 2:cellymax-1
        for i in 2:cellxmax-1
            D[i,j] = dt/dtau[i,j] + 1.0 + dt*(A_beta_shig[i,j]+2*jalphaP[i,j] + B_beta_shig[i,j]+2*jbetaP[i,j])
        end
    end

    # RHS
    for l in 1:nval
        for j in 2:cellymax-1
            for i in 2:cellxmax-1
                RHS_temp[i,j,l] = - (Qcon_hat[i,j,l]-Qconn_hat[i,j,l])*1.0 + dt*RHS[i,j,l]
            end
        end
    end

    # lower sweep
    for l in 1:nval
        for j in 2:cellymax-1
            for i in 2:cellxmax-1
                delta_Q[i,j,l] = D[i,j]^(-1) * (RHS_temp[i,j,l]+LdQ[i,j,l])
            end
        end
    end

    # upepr sweep
    for l in 1:nval
        for j in 2:cellymax-1
            for i in 2:cellxmax-1
                delta_Q[i,j,l] = delta_Q[i,j,l] - D[i,j]^(-1) * UdQ[i,j,l]
            end
        end
    end

    # reset
    for l in 1:nval
        for j in 2:cellymax-1
            for i in 2:cellxmax-1
                LdQ[i,j,l] = 0.0
                UdQ[i,j,l] = 0.0
            end
        end
    end

    return delta_Q
end


text:
using LinkedLists
using Test

@time @testset "Singly-Linked Lists" begin include("testslist.jl") end
@time @testset "Doubly-Linked Lists" begin include("testlist.jl") end


text:
module JuliADAPT

include("fci.jl")

end

max_stars_repo_name: hrgrimsl/JuliADAPT

max_issues_repo_name: hrgrimsl/JuliADAPT

max_forks_repo_name: hrgrimsl/JuliADAPT


text:
"""
    modular_inverse(base::Int64, modulus::Int64)

Compute the multiplicative inverse of the `base` under the `modulus`.

# Arguments
- `base::Int64`: base
- `modulus ::Int64`: modulus

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> ClassicAlgorithmsCollections.modular_inverse(3, 11)
4
```
"""
function modular_inverse(base::Int64, modulus::Int64)
    base = base % modulus
    for i in 1:modulus
        if (base * i) % modulus == 1
            return i
        end
    end
    return 1
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
module cd

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
using LinearAlgebra
import LinearAlgebra:checksquare

"""
    function sweep_op_internal!(A::AbstractMatrix{Float64}, k::Int64)

    Implement the naive sweep operator described by Goodnight, J. (1979). "A Tutorial on the SWEEP Operator." The American Statistician.
"""
function sweep_op_internal!(A::AbstractMatrix{Float64}, k::Int64)
    p = checksquare(A)
    if k >= p
        throw(ArgumentError("Incorrect k value"))
    end

    @inbounds D = A[k,k] # step 1 D = Aₖₖ

    if D == zero(eltype(A))
        throw(error("the element $k,$k is zero"))
    end

    @inbounds for i in 1:p
        A[k,i] = A[k,i] / D   # step 2 divide row k by D
    end

    @inbounds for i in 1:(k - 1) # step 3: for every other row  i != k
        B = A[i,k]   # let B = Aᵢₖ
        for j in 1:p
            A[i,j] = A[i,j] - B * A[k, j]  # Subtract B * row k from row i
        end
        A[i,k] = - B / D # set Aᵢₖ = -B/D
    end

    @inbounds for i in k + 1:p # step 3: for every other row  i != k
        B = A[i,k]   # let B = Aᵢₖ
        for j in 1:p
            A[i,j] = A[i,j] - B * A[k, j]  # Subtract B * row k from row i
        end
        A[i,k] = - B / D # set Aᵢₖ = -B/D
    end

    @inbounds A[k,k] = 1 / D # step 4 Set Aₖₖ = 1/D

    return A
end

function sweep_op!(A::AbstractMatrix{Float64}, k::Int64)
    sweep_op_internal!(A, k)
end

function sweep_op!(A::AbstractMatrix{Float64}, ks::AbstractVector{Int64})
    for k in ks
        sweep_op_internal!(A, k)
    end
    return A
end

"""
    function sweep_op_full!(A::AbstractMatrix{Float64})

    (internal) Get SSE, error sum of squares, for all p-1 models. The last one is the SSE for the full model.
"""
function sweep_op_full!(A::AbstractMatrix{Float64})
    n , p = size(A)

    SSEv = Vector{Float64}(undef, p-1)

    for k in 1:p-1
        sweep_op_internal!(A, k)
        SSEv[k] = A[p,p]
    end
    return SSEv
end


text:
using ExpFamily
using Base.Test

dim = 4
srand(123)

m1  = randn(dim)
m2  = randn(dim)
C1  = randn(dim,dim)
C1 *= C1'
C2  = randn(dim,dim)
C2 *= C2'

gNPa = GaussianNatParam(mean=m1, cov=C1)
gNPb = GaussianNatParam(mean=m2, cov=C2)
gMPa = GaussianMeanParam(mean=m1, cov=C1)
gMPb = GaussianMeanParam(mean=m2, cov=C2)

x   = randn(dim)
ns  = 5
xb  = randn(dim,ns)

Ca = cov(gNPa)
ma = mean(gNPa)
l  = exp( -dot(x-ma, Ca\(x-ma)) / 2 )/sqrt( (2pi)^dim * det(Ca) )
ll = log(l)

lb  = [exp( -dot(xb[:,i]-ma, Ca\(xb[:,i]-ma))/2) for i in 1:ns]
lb /= sqrt( (2pi)^dim * det(Ca) )
llb = log.(lb)

ul  = -0.5(dot(x, Ca\x)-2dot(x, Ca\ma))
ulb = [-0.5(dot(xb[:,i], Ca\xb[:,i])-2dot(xb[:,i], Ca\ma))
        for i in 1:ns]

Cda   = diag(Ca)
gdNPa = DiagGaussianNatParam(mean=ma, cov=Cda)
gdMPa = DiagGaussianMeanParam(mean=ma, cov=Cda)

ldb  = [exp( -dot(xb[:,i]-ma, (xb[:,i]-ma)./Cda)/2) for i in 1:ns]
ldb /= sqrt( (2pi)^dim * prod(Cda) )
lldb = log.(ldb)

uldb  = [-0.5(dot(xb[:,i], xb[:,i]./Cda)-2dot(xb[:,i], ma./Cda))
            for i in 1:ns]

llNPa   = loglik(gNPa,  x)
llNPab  = loglik(gNPa,  xb)
lldNPab = loglik(gdNPa, xb)
llMPab  = loglik(gMPa, xb)

@test isapprox( llNPa,   ll  )
@test isapprox( llNPab,  llb )
@test isapprox( llNPab,  loglik(gNPa, xb', axis=1))
@test isapprox( llNPab,  llMPab )
@test isapprox( llMPab,  loglik(gMPa, xb', axis=1))
@test isapprox( lldNPab, lldb )
@test isapprox( lldNPab, loglik(gdNPa, xb', axis=1) )
@test isapprox( loglik(gdMPa, xb), lldNPab)

@test isapprox( uloglik(gNPa, x), ul)
@test isapprox( uloglik(gNPa, xb), ulb )
@test isapprox( uloglik(gNPa, xb', axis=1), ulb)
@test isapprox( uloglik(gMPa, xb), ulb)
@test isapprox( uloglik(gMPa, xb', axis=1), ulb)

@test isapprox( uloglik(gdNPa, xb), uldb)
@test isapprox( uloglik(gdNPa, xb', axis=1), uldb)
@test isapprox( uloglik(gdMPa, xb), uldb)
@test isapprox( uloglik(gdMPa, xb', axis=1), uldb)

@test isapprox( gradloglik(gNPa, x), Ca\(ma-x) )
@test isapprox( gradloglik(gMPa, x), gradloglik(gNPa, x) )


text:
using Distributions
using ConjugatePriors

using Random: seed!

using ConjugatePriors: NormalGamma, NormalInverseGamma, NormalInverseChisq
using ConjugatePriors: posterior, posterior_rand, posterior_mode, posterior_randmodel, fit_map

n = 100
w = rand(100)

@testset "Conjugates for normal distribution" begin

    @testset "Νormal - Νormal (known sigma)" begin

        pri = Normal(1.0, 5.0)

        x = rand(Normal(2.0, 3.0), n)
        p = posterior((pri, 3.0), Normal, x)
        @test isa(p, Normal)
        @test mean(p) ≈ (mean(pri) / var(pri) + sum(x) / 9.0) / (1.0 / var(pri) + n / 9.0)
        @test var(p) ≈ inv(1.0 / var(pri) + n / 9.0)

        r = posterior_mode((pri, 3.0), Normal, x)
        @test r ≈ mode(p)

        f = fit_map((pri, 3.0), Normal, x)
        @test isa(f, Normal)
        @test f.μ == r
        @test f.σ == 3.0

        p = posterior((pri, 3.0), Normal, x, w)
        @test isa(p, Normal)
        @test mean(p) ≈  (mean(pri) / var(pri) + dot(x, w) / 9.0) / (1.0 / var(pri) + sum(w) / 9.0)
        @test var(p) ≈ inv(1.0 / var(pri) + sum(w) / 9.0)

        r = posterior_mode((pri, 3.0), Normal, x, w)
        @test r ≈ mode(p)

        f = fit_map((pri, 3.0), Normal, x, w)
        @test isa(f, Normal)
        @test f.μ == r
        @test f.σ == 3.0

    end

    @testset "ΙnverseGamma - Νormal (known mu)" begin

        pri = InverseGamma(1.5, 0.5)

        x = rand(Normal(2.0, 3.0), n)
        p = posterior((2.0, pri), Normal, x)
        @test isa(p, InverseGamma)
        @test shape(p) ≈ shape(pri) + n / 2
        @test scale(p) ≈ scale(pri) + sum(abs2.(x .- 2.0)) / 2

        r = posterior_mode((2.0, pri), Normal, x)
        @test r ≈ mode(p)

        f = fit_map((2.0, pri), Normal, x)
        @test isa(f, Normal)
        @test f.μ == 2.0
        @test abs2(f.σ) ≈ r

        p = posterior((2.0, pri), Normal, x, w)
        @test isa(p, InverseGamma)
        @test shape(p) ≈ shape(pri) + sum(w) / 2
        @test scale(p) ≈ scale(pri) + dot(w, abs2.(x .- 2.0)) / 2

        r = posterior_mode((2.0, pri), Normal, x, w)
        @test r ≈ mode(p)

        f = fit_map((2.0, pri), Normal, x, w)
        @test isa(f, Normal)
        @test f.μ == 2.0
        @test abs2(f.σ) ≈ r

    end

    @testset "Gamma - Νormal (known mu)" begin

        pri = Gamma(1.5, 2.0)

        x = rand(Normal(2.0, 3.0), n)
        p = posterior((2.0, pri), Normal, x)
        @test isa(p, Gamma)
        @test shape(p) ≈ shape(pri) + n / 2
        @test scale(p) ≈ scale(pri) + sum(abs2.(x .- 2.0)) / 2

        r = posterior_mode((2.0, pri), Normal, x)
        @test r ≈ mode(p)

        f = fit_map((2.0, pri), Normal, x)
        @test isa(f, Normal)
        @test f.μ == 2.0
        @test abs2(f.σ) ≈ inv(r)

        p = posterior((2.0, pri), Normal, x, w)
        @test isa(p, Gamma)
        @test shape(p) ≈ shape(pri) + sum(w) / 2
        @test scale(p) ≈ scale(pri) + dot(w, abs2.(x .- 2.0)) / 2

        r = posterior_mode((2.0, pri), Normal, x, w)
        @test r ≈ mode(p)

        f = fit_map((2.0, pri), Normal, x, w)
        @test isa(f, Normal)
        @test f.μ == 2.0
        @test abs2(f.σ) ≈ inv(r)

    end

    @testset "NormalInverseGamma - Normal" begin

        mu_true = 2.
        sig2_true = 3.
        x = rand(Normal(mu_true, sig2_true), n)

        mu0 = 2.
        v0 = 3.
        shape0 = 5.
        scale0 = 2.
        pri = NormalInverseGamma(mu0, v0, shape0, scale0)

        post = posterior(pri, Normal, x)
        @test isa(post, NormalInverseGamma)

        @test post.mu ≈ (mu0/v0 + n*mean(x))/(1.0/v0 + n)
        @test post.v0 ≈ 1.0/(1.0/v0 + n)
        @test post.shape ≈ shape0 + 0.5*n
        @test post.scale ≈ scale0 + 0.5*(n-1)*var(x) + n./v0./(n + 1.0/v0)*0.5*(mean(x)-mu0).^2

        ps = posterior_randmodel(pri, Normal, x)

        @test isa(ps, Normal)
        @test insupport(ps,ps.μ) && ps.σ > zero(ps.σ)

        for shape in [0.1, 1., 10.]
            d = NormalInverseGamma(mu0, v0, scale0, shape)
            μ, σ2 = mode(d)
            @test pdf(d, μ, σ2) > pdf(d, μ, σ2 + 0.001)
            @test pdf(d, μ, σ2) > pdf(d, μ, σ2 - 0.001)
        end

    end

    @testset "NormalInverseChisq - Normal" begin

        mu_true = 2.
        sig2_true = 3.
        x = rand(Normal(mu_true, sig2_true), n)

        μ0 = 2.0
        σ20 = 2.0/5.0
        κ0 = 1.0/3.0
        ν0 = 10.0

        pri = NormalInverseChisq(μ0, σ20, κ0, ν0)
        pri2 = convert(NormalInverseGamma, pri)

        @test convert(NormalInverseChisq, pri2) == pri

        @test mode(pri2) == mode(pri)
        @test mean(pri2) == mean(pri)
        @test pdf(pri, mu_true, sig2_true) == pdf(pri2, mu_true, sig2_true)

        @test (seed!(1); rand(pri)) == (seed!(1); rand(pri2))

        # check that updating is consistent between NIχ2 and NIG
        post = posterior(pri, Normal, x)
        post2 = posterior(pri2, Normal, x)
        @test isa(post, NormalInverseChisq)
        @test convert(NormalInverseChisq, post2) == post

        for _ in 1:10
            x = rand(post)
            @test pdf(post, x...) ≈ pdf(post2, x...)
            @test logpdf(post, x...) ≈ logpdf(post2, x...)
        end

        for ν in [0.1, 1., 10.]
            nix = NormalInverseChisq(0., 2., 3., ν)
            μ, σ2 = mode(nix)
            @test pdf(nix, μ, σ2) > pdf(nix, μ, σ2 + 0.001)
            @test pdf(nix, μ, σ2) > pdf(nix, μ, σ2 - 0.001)
        end

    end

    @testset "NormalGamma - Normal" begin

        mu_true = 2.
        tau2_true = 3.
        x = rand(Normal(mu_true, 1.0/tau2_true), n)

        mu0 = 2.
        nu0 = 3.
        shape0 = 5.
        rate0 = 2.
        pri = NormalGamma(mu0, nu0, shape0, rate0)

        post = posterior(pri, Normal, x)
        @test isa(post, NormalGamma)

        @test post.mu ≈ (nu0*mu0 + n*mean(x))./(nu0 + n)
        @test post.nu ≈ nu0 + n
        @test post.shape ≈ shape0 + 0.5*n
        @test post.rate ≈ rate0 + 0.5*(n-1)*var(x) + n*nu0/(n + nu0)*0.5*(mean(x)-mu0).^2

        ps = posterior_randmodel(pri, Normal, x)

        @test isa(ps, Normal)
        @test insupport(ps, ps.μ) && ps.σ > zero(ps.σ)
    end

end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

# The order in which files are included is important,
# since certain things need to be defined before others use them

## Interval type

abstract type AbstractInterval{T <: AbstractFloat} end

struct EmptyInterval{T} <: AbstractInterval{T} end

struct NotAnInterval{T} <: AbstractInterval{T} end
nai(::Type{T}) where T = NotAnInterval{T}()
nai() = nai(Float64)

struct Interval{T} <: AbstractInterval{T}
    lo :: T
    hi :: T

    function Interval(lo::T, hi::T, ::Val{CHECK}) where {T <: AbstractFloat, CHECK}
        if CHECK
            @assert !isnan(lo)
            @assert !isnan(hi)
            @assert lo <= hi
            @assert lo < Inf
            @assert hi > -Inf
        end
        return new{T}(lo, hi)
    end
end

# all intervals except NotAnInterval
const ValidInterval{T} = Union{Interval{T}, EmptyInterval{T}}

emptyinterval(::Type{T}) where T = EmptyInterval{T}()
emptyinterval(x::ValidInterval{T}) where T = emptyinterval(T)
emptyinterval() = emptyinterval(Float64)
const ∅ = emptyinterval()

# check boundaries by default
Interval(lo::T, hi::T) where T = Interval(lo, hi, Val(true))

# methods which call unsafe_interval must guarantee validity of inputs!
unsafe_interval(lo::T, hi::T) where T = Interval(lo, hi, Val(false))
unsafe_interval(x::T) where T = unsafe_interval(x, x)

# flexible constructor, methods calling this constructor may produce nai()!
function interval(lo::T, hi::T) where T <: AbstractFloat
    if lo > hi || lo == Inf || hi == -Inf
        return emptyinterval(T)
    end
    if isnan(lo) || isnan(hi)
        return nai(T)
    end
    return unsafe_interval(lo, hi)
end

# thin interval constructor
interval(a::T) where T <: AbstractFloat =
    isfinite(a) ? unsafe_interval(a, a) : nai(T)

interval(a::T, b::T) where T <: Real = interval(float(T), a, b)
interval(a::T) where T <: Real = interval(float(T), a, a)

# type conversion
interval(::Type{T}, lo, hi) where T =
    interval(T(lo, RoundDown), T(hi, RoundUp))
interval(::Type{T}, x::Real) where T = isnan(x) ? nai(T) : unsafe_atomic(T, x)

# automatic type promotion
interval(a::Tuple{T, T}) where T <: Real = interval(a[1], a[2])
interval(a::T, b::S) where {T <: Real, S <: Real} = interval(promote(a, b))

# for b < 0 this produces an empty interval
±(a::Real, b::Real) = interval(a) ± b # type-piracy?
±(a::Interval{T}, b::AbstractFloat) where T = a + interval(T, -b, b)

# common intervals
unit_ball(::Type{T}) where T = unsafe_interval(-one(T), one(T))
nonnegatives(::Type{T}) where T = unsafe_interval(zero(T), T(Inf))

"""`entireinterval`s represent the whole Real line: [-∞, ∞]."""
entireinterval(::Type{T}) where T = Interval(T(-Inf), T(Inf))
entireinterval(::ValidInterval{T}) where T = entireinterval(T)
entireinterval() = entireinterval(Float64)

show(io::IO, x::Interval) = print(io, string("Interval(", x.lo, ",", x.hi, ")"))
show(io::IO, x::EmptyInterval{T}) where T =
    print(io, string("emptyinterval(", T, ")"))
show(io::IO, x::NotAnInterval{T}) where T = print(io, string("nai(", T, ")"))


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct AVRSimple <: AVR
        Kv::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of a simple proportional AVR in the derivative of EMF
i.e. an integrator controller on EMF

# Arguments
- `Kv::Float64`: Proportional Gain, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`: Fixed AVR has 1 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct AVRSimple <: AVR
    "Proportional Gain"
    Kv::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    "Fixed AVR has 1 states"
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function AVRSimple(Kv, ext=Dict{String, Any}(), )
    AVRSimple(Kv, ext, [:Vf], 1, InfrastructureSystemsInternal(), )
end

function AVRSimple(; Kv, ext=Dict{String, Any}(), )
    AVRSimple(Kv, ext, )
end

# Constructor for demo purposes; non-functional.
function AVRSimple(::Nothing)
    AVRSimple(;
        Kv=0,
        ext=Dict{String, Any}(),
    )
end

"""Get AVRSimple Kv."""
get_Kv(value::AVRSimple) = value.Kv
"""Get AVRSimple ext."""
get_ext(value::AVRSimple) = value.ext
"""Get AVRSimple states."""
get_states(value::AVRSimple) = value.states
"""Get AVRSimple n_states."""
get_n_states(value::AVRSimple) = value.n_states
"""Get AVRSimple internal."""
get_internal(value::AVRSimple) = value.internal


text: println("Testing...")


text:
using Example
using Test

tests = [
    "foo",
]

@testset "Example.jl" begin
    for t in tests
        include("$(t).jl")
    end
end


text:
using JSON

mutable struct Randomizer
    randomization_config::Dict
    default_config::Dict
    key_s::Set
end

function Randomizer(randomization_config_fp="default_dr.json", default_config_fp="default.json")
    randomization_config, default_config = nothing, nothing

    try
        open(get_file_path("src/randomization/config", randomization_config_fp, "json"), "r") do f
            randomization_config = json.load(f)
        end
    catch
        #logger.warning("Couldn't find {} in randomization/config subdirectory".format(randomization_config_fp))
        randomization_config = Dict()
    end

    open(get_file_path("src/randomization/config", default_config_fp, "json"), "r") do f
        default_config = JSON.parse(f)
    end

    key_s = Set(append!(collect.(keys.((randomization_config, default_config)))...))

    Randomizer(randomization_config, default_config, key_s)
end

function randomize(randomizer::Randomizer)
    ##Returns a dictionary of randomized parameters, with key: parameter name and value: randomized
    #value
    ##
    randomization_settings = Dict()

    for k in randomizer.key_s
        setting = nothing
        if k in keys(randomizer.randomization_config)
            randomization_definition = randomizer.randomization_config[k]

            if randomization_definition["type"] == "int"
                try
                    low = randomization_definition["low"]
                    high = randomization_definition["high"]
                    size = get(randomization_definition, "size", 1)
                catch
                    throw(KeyError("Please check your randomization definition for: $k"))
                end
                setting = rand(low:high, size)

            elseif randomization_definition["type"] == "uniform"
                try
                    low = randomization_definition["low"]
                    high = randomization_definition["high"]
                    size = get(randomization_definition, "size", 1)
                catch
                    throw(KeyError("Please check your randomization definition for: $k"))
                end

                setting = rand(low:high, size)

            elseif randomization_definition["type"] == "normal"
                try
                    loc = randomization_definition["loc"]
                    scale = randomization_definition["scale"]
                    size = get(randomization_definition, "size", 1)
                catch
                    throw(KeyError("Please check your randomization definition for: $k"))
                end

                setting = rand(Normal(loc, scale), size)

            else
                throw(NotImplementedError("You've specified an unsupported distribution type"))
            end

        elseif k in keys(randomizer.default_config)
            randomization_definition = randomizer.default_config[k]
            setting = randomization_definition["default"]
        end

        randomization_settings[k] = setting
    end
    return randomization_settings
end


text:
using LegendrePolynomials
using Documenter

DocMeta.setdocmeta!(LegendrePolynomials, :DocTestSetup, :(using LegendrePolynomials); recursive=true)

makedocs(;
    modules=[LegendrePolynomials],
    authors="Jishnu Bhattacharya",
    repo="https://github.com/jishnub/LegendrePolynomials.jl/blob/{commit}{path}#L{line}",
    sitename="LegendrePolynomials",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jishnub.github.io/LegendrePolynomials.jl",
        assets=String[],
    ),
    pages=[
        "Introduction" => "index.md",
        "Tutorial" => "tutorial.md",
        "Derivatives" => "derivatives.md",
    ],
)

deploydocs(;
    repo="github.com/jishnub/LegendrePolynomials.jl",
)


text:
@testset "4.7.4 x^m (a+b trig^n)^p" begin
    (a, b, c, d, x, ) = @variables a b c d x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*x^m*(a+b*trig(c+d*x)^n)^p=#

    #= ::Section:: =#
    #=Integrands*of*the*form*x^m*(a+b*sin(c+d*x)^n)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(a+b*sin(c+d*x)^2)^p=#

    @test_int [x/(a + b*sin(x)^2), x, 9, -((im*x*log(1 - (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b))) + (im*x*log(1 - (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b)) - PolyLog(2, (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b)))/(4*sqrt(a)*sqrt(a + b)) + PolyLog(2, (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b)))/(4*sqrt(a)*sqrt(a + b))]
    @test_int [x^2/(a + b*sin(x)^2), x, 11, -((im*x^2*log(1 - (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b))) + (im*x^2*log(1 - (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b)) - (x*PolyLog(2, (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b)) + (x*PolyLog(2, (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b)) - (im*PolyLog(3, (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(4*sqrt(a)*sqrt(a + b)) + (im*PolyLog(3, (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(4*sqrt(a)*sqrt(a + b))]
    @test_int [x^3/(a + b*sin(x)^2), x, 13, -((im*x^3*log(1 - (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b))) + (im*x^3*log(1 - (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b)) - (3*x^2*PolyLog(2, (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(4*sqrt(a)*sqrt(a + b)) + (3*x^2*PolyLog(2, (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(4*sqrt(a)*sqrt(a + b)) - (3*im*x*PolyLog(3, (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(4*sqrt(a)*sqrt(a + b)) + (3*im*x*PolyLog(3, (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(4*sqrt(a)*sqrt(a + b)) + (3*PolyLog(4, (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(8*sqrt(a)*sqrt(a + b)) - (3*PolyLog(4, (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(8*sqrt(a)*sqrt(a + b))]

    @test_int [x/(a + b*sin(c + d*x)^2)^2, x, 12, -((im*(2*a + b)*x*log(1 - (b*ℯ^(2*im*(c + d*x)))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(4*a^(3/2)*(a + b)^(3/2)*d)) + (im*(2*a + b)*x*log(1 - (b*ℯ^(2*im*(c + d*x)))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(4*a^(3/2)*(a + b)^(3/2)*d) - log(2*a + b - b*cos(2*c + 2*d*x))/(4*a*(a + b)*d^2) - ((2*a + b)*PolyLog(2, (b*ℯ^(2*im*(c + d*x)))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(8*a^(3/2)*(a + b)^(3/2)*d^2) + ((2*a + b)*PolyLog(2, (b*ℯ^(2*im*(c + d*x)))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(8*a^(3/2)*(a + b)^(3/2)*d^2) + (b*x*sin(2*c + 2*d*x))/(2*a*(a + b)*d*(2*a + b - b*cos(2*c + 2*d*x)))]

    @test_int [x*sqrt(sin(x)^2), x, 3, sqrt(sin(x)^2) - x*cot(x)*sqrt(sin(x)^2)]

    #= ::Section:: =#
    #=Integrands*of*the*form*x^m*(a+b*cos(c+d*x)^n)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(a+b*cos(c+d*x)^2)^p=#

    @test_int [x/(a + b*cos(x)^2), x, 9, -((im*x*log(1 + (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b))) + (im*x*log(1 + (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b)) - PolyLog(2, -((b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(4*sqrt(a)*sqrt(a + b)) + PolyLog(2, -((b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(4*sqrt(a)*sqrt(a + b))]
    @test_int [x^2/(a + b*cos(x)^2), x, 11, -((im*x^2*log(1 + (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b))) + (im*x^2*log(1 + (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b)) - (x*PolyLog(2, -((b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b)))))/(2*sqrt(a)*sqrt(a + b)) + (x*PolyLog(2, -((b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b)))))/(2*sqrt(a)*sqrt(a + b)) - (im*PolyLog(3, -((b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b)))))/(4*sqrt(a)*sqrt(a + b)) + (im*PolyLog(3, -((b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b)))))/(4*sqrt(a)*sqrt(a + b))]
    @test_int [x^3/(a + b*cos(x)^2), x, 13, -((im*x^3*log(1 + (b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b))) + (im*x^3*log(1 + (b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(2*sqrt(a)*sqrt(a + b)) - (3*x^2*PolyLog(2, -((b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b)))))/(4*sqrt(a)*sqrt(a + b)) + (3*x^2*PolyLog(2, -((b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b)))))/(4*sqrt(a)*sqrt(a + b)) - (3*im*x*PolyLog(3, -((b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b)))))/(4*sqrt(a)*sqrt(a + b)) + (3*im*x*PolyLog(3, -((b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b)))))/(4*sqrt(a)*sqrt(a + b)) + (3*PolyLog(4, -((b*ℯ^(2*im*x))/(2*a + b - 2*sqrt(a)*sqrt(a + b)))))/(8*sqrt(a)*sqrt(a + b)) - (3*PolyLog(4, -((b*ℯ^(2*im*x))/(2*a + b + 2*sqrt(a)*sqrt(a + b)))))/(8*sqrt(a)*sqrt(a + b))]

    @test_int [x/(a + b*cos(c + d*x)^2)^2, x, 12, -((im*(2*a + b)*x*log(1 + (b*ℯ^(2*im*(c + d*x)))/(2*a + b - 2*sqrt(a)*sqrt(a + b))))/(4*a^(3/2)*(a + b)^(3/2)*d)) + (im*(2*a + b)*x*log(1 + (b*ℯ^(2*im*(c + d*x)))/(2*a + b + 2*sqrt(a)*sqrt(a + b))))/(4*a^(3/2)*(a + b)^(3/2)*d) - log(2*a + b + b*cos(2*c + 2*d*x))/(4*a*(a + b)*d^2) - ((2*a + b)*PolyLog(2, -((b*ℯ^(2*im*(c + d*x)))/(2*a + b - 2*sqrt(a)*sqrt(a + b)))))/(8*a^(3/2)*(a + b)^(3/2)*d^2) + ((2*a + b)*PolyLog(2, -((b*ℯ^(2*im*(c + d*x)))/(2*a + b + 2*sqrt(a)*sqrt(a + b)))))/(8*a^(3/2)*(a + b)^(3/2)*d^2) - (b*x*sin(2*c + 2*d*x))/(2*a*(a + b)*d*(2*a + b + b*cos(2*c + 2*d*x)))]
    end

max_stars_repo_path: test/4 Trig functions/4.7 Miscellaneous/4.7.4 x^m (a+b trig^n)^p.jl


text:
function _numeric_f(integrand, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, use78=false, ϕm=π)
    tgα = -tgα
    ϕl = mylinspace(θ, ϕm, N)
    u0 = 1/rs
    du0 = dudθ(θ, tgα, rs)
    start = [u0; du0];

    ϕ, y = use78?
        ode78(integrand, start, ϕl, abstol=abstol, reltol=reltol) :
        ode45(integrand, start, ϕl, abstol=abstol, reltol=reltol);
    u = map(y -> y[1], y);
    du = map(y -> y[2], y);
    ixs = find(u.>0)
    u=u[ixs]
    du=du[ixs]
    ϕ=ϕ[ixs]

    x = [v for v in -cos(ϕ)./u]
    y = [v for v in sin(ϕ)./u]

    out_angle =  ray_angle(ϕ[end], 1/u[end], du[end])
    in_angle =  atan(tgα)

    x, y, u, du, -in_angle, out_angle, -out_angle-in_angle, ϕ
end

function numeric(massfunc, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π)
    function integrand(ϕ, y)
        (u, v) = y
        u_prime = v
        v_prime = 3*rg(massfunc(1./u)).*u.*u-u
        [u_prime; v_prime]
    end
    _numeric_f(integrand, θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)
end

function numeric_el(massfunc, rhofunc, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π)
    M200 = massfunc(1.0e30) # infinity, ha!
    du = dudθ(θ, tgα, rs)
    u = 1/rs
    _rg = rg(M200)
    _e = 1/(1-2*_rg*u)
    e2l2 = du*du + u*u*_e
    function integrand(ϕ, y)
        (u, v) = y
        _rg = rg(massfunc(1./u))
        u_prime = v
        v_prime = 3*_rg.*u.*u-u+(e2l2-u.*u)*4*π*G*rhofunc(1./u)/C/C/u/u/u
        [u_prime; v_prime]
    end
    _numeric_f(integrand, θ, rs, tgα, N; abstol=abstol, reltol=reltol, use78=true, ϕm=ϕm)
end

# some convenience functions:

numeric_tiso200b(z, M200b, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric(r->tiso_m(M200b, r, R200b(z, M200b)), θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)
numeric_tiso200c(z, M200c, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric(r->tiso_m(M200c, r, R200c(z, M200c)), θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)
numeric_tiso_el200b(z, M200b, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric_el(r->tiso_m(M200b, r, R200b(z, M200b)), r->tiso_ρ(M200b, r, R200b(z, M200b)),
                θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)
numeric_tiso_el200c(z, M200c, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric_el(r->tiso_m(M200c, r, R200c(z, M200c)), r->tiso_ρ(M200c, r, R200c(z, M200c)),
                θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)
numeric_tisoVir(z, Mvir, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric(r->tiso_m(Mvir, r, RVIR(z, MVir)), θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)
numeric_tiso_elVir(z, Mvir, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric_el(r->tiso_m(Mvir, r, RVIR(z, Mvir)), r->tiso_ρ(Mvir, r, RVIR(z, Mvir)),
                θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)

numeric_tnfw200b(z, M200b, Rscale, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric(r->tnfw_m(M200b, Rcalse, r, R200b(z, M200b)), θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)
numeric_tnfw200c(z, M200c, Rscale, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric(r->tnfw_m(M200c, Rscale, r, R200c(z, M200c)), θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)
numeric_tnfw_el200b(z, M200b, Rscale, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric_el(r->tnfw_m(M200b, Rscale, r, R200b(z, M200b)), r->tnfw_ρ(M200b, Rscale, r, R200b(z, M200b)),
                θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)
numeric_tnfw_el200c(z, M200c, Rscale, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric_el(r->tnfw_m(M200c, Rscale, r, R200c(z, M200c)), r->tnfw_ρ(M200c, Rscale, r, R200c(z, M200c)),
                θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)
numeric_tnfwVir(z, Mvir, Rscale, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric(r->tnfw_m(Mvir, Rscale, r, RVIR(z, MVir)), θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)
numeric_tnfw_elVir(z, Mvir, Rscale, θ, rs, tgα, N=100; abstol=1.0e-16, reltol=1.0e-17, ϕm=π) =
    numeric_el(r->tnfw_m(Mvir, Rscale, r, RVIR(z, Mvir)), r->tnfw_ρ(Mvir, Rscale, r, RVIR(z, Mvir)),
                θ, rs, tgα, N; abstol=abstol, reltol=reltol, ϕm=ϕm)

max_stars_repo_name: skariel/Lensing

max_issues_repo_name: skariel/Lensing

max_forks_repo_name: skariel/Lensing