text:
"""
    error_se(target,prediction)

Compute the squared error between the `prediction` and `target` values.

"""
function error_se(target,prediction)
    error_se = zeros(size(target,1),size(target,2));
    for j in 1:size(target,1)
        for i in 1:size(target,2)
            error_se[j,i] = (prediction[j,i] - target[j,i])^2;
        end
    end
    return error_se
end

"""
    error_rel(target,prediction)

Compute the relative error between the `prediction` and `target` values.

"""
function error_rel(target,prediction)
    error_rel = zeros(size(target,1),size(target,2));
    for j in 1:size(target,1)
        for i in 1:size(target,2)
            error_test_rel[j,i] = (prediction[j,i] - target[j,i])/target[j,i];
        end
    end
    return error_rel
end

"""
    mse_error(target,prediction)

Compute the mean squared error between the `prediction` and `target` values.

"""
function mse_error(target,prediction)
    error_test_mse = 0.0;
    for j in 1:size(target,1)
        for i in 1:size(target,2)
            error_test_mse += (prediction[j,i] - target[j,i])^2;
        end
    end
    # return (1/size(target[:],1))*sum((prediction[:].-target[:]).^2);
    return (1/(size(target,1)*size(target,2)))*error_test_mse
end

"""
    norm_rel_error(target,prediction)

Compute the Euclidean two-norm relative error between the `prediction` and `target` values.

"""
function norm_rel_error(target,prediction)
    if size(target) != size(prediction)
        error("Target size and prediction size do not match!")
    end
    error = zeros(size(target,1))
    for i in 1:size(target,1)
        error[i] = norm(prediction[i,:]-target[i,:])/norm(target[i,:]);
    end
    return error
end

"""
    function norm_rel_error_continuous(target,prediction,weights)

Compute the continuous two-norm relative error between the `prediction` and `target` values.

"""
function norm_rel_error_continuous(target,prediction,weights)
  error = zeros(size(target,1))
  W = diagm(0=>weights);
  for i in 1:size(target,1)
    diff = prediction[i,:]-target[i,:];
    error[i] = sqrt(diff'*W*diff)/sqrt(target[i,:]'*W*target[i,:]);
  end
  return error
end

"""
    norm_infinity_error(target,prediction)

Compute the two-norm relative error between the `prediction` and `target` values.

"""
function norm_infinity_error(target,prediction)
  error = zeros(size(target,1))
  for i in 1:size(target,1)
    diff = prediction[i,:]-target[i,:];
    error[i] = maximum(abs.(diff),dims=1)[1];
  end
  return error
end

"""
    average_error(domain,error)

Compute the average error using the trapezoid rule \$\\frac{1}{T} \\int_0^T error(t) dt\$.

"""
function average_error(domain,error)
    return 1/(domain[end]-domain[1])*trapz(domain,error)
end

"""
    ic_error(target,prediction)

Compute the relative error between the `prediction` and `target` values for an initial condition. If the `target` = 0, the `target` in denomenator is augmented by \$\\epsilon_{machine}\$.

"""
function ic_error(target,prediction)
    error = zeros(size(target,1));
    for i in 1:size(target,1)
        error[i] = (prediction[i]-target[i])/(target[i] + eps());
    end
    return error
end

"""
    average_ic_error(target,prediction)

Compute the two-norm relative error between the `prediction` and `target` values for an initial condition.

"""
function average_ic_error(target,prediction)
    return norm(prediction-target)/norm(target)
end

"""
    periodic_fill_domain(x_locations)

Output the full domain from periodic domain specified for `x_locations`.

"""
function periodic_fill_domain(x_locations)
    x_locations_full = vcat(x_locations,(x_locations[end]+(x_locations[end]-x_locations[end-1])))
    return x_locations_full
end

"""
    periodic_fill_solution(solution)

Output the full ``u(t,x)`` solution from periodic `solution`.

"""
function periodic_fill_solution(solution)
    solution_full = zeros(size(solution,1),(size(solution,2)+1));
    for i in 1:size(solution,1)
        solution_full[i,1:end-1] = solution[i,:];
        solution_full[i,end] = solution[i,1];
    end
    return solution_full
end

"""
    solution_interpolation(t_span_original,x_locations_original,t_span_interpolate,x_locations_interpolate,solution)

Compute the ``u(t,x)`` solution at intermediate ``(t,x)`` locations using linear interpolation.

"""
function solution_interpolation(t_span_original,x_locations_original,t_span_interpolate,x_locations_interpolate,solution)
    solution_interpolate_fit = interpolate((t_span_original,x_locations_original),solution,Gridded(Linear()));
    solution_interpolate = zeros(size(t_span_interpolate,1),size(x_locations_interpolate,1));
    for i in 1:size(t_span_interpolate,1)
        for j in 1:size(x_locations_interpolate,1)
            solution_interpolate[i,j] = solution_interpolate_fit(t_span_interpolate[i],x_locations_interpolate[j]);
        end
    end
    return solution_interpolate
end

"""
    reduced_initial_condition(L1,L2,x_reduced,x_locations,initial_condition)

Extract the ``x`` locations and intial condition values ``u(x)`` at a reduced number of equally spaced spatial locations.

"""
function reduced_initial_condition(L1,L2,x_reduced,x_locations,initial_condition)
    # dL = abs(L2-L1);
    # j = reduce(vcat,[0:1:N-1]);
    # x_reduced = (dL.*j)./N;
    ind = [findall(x_reduced[i].==x_locations)[1] for i in 1:length(x_reduced)];
    ic_reduced = zeros(size(ind,1));
    for i in 1:size(ind,1)
        ic_reduced[i] = initial_condition[ind[i]];
    end
    return ic_reduced
end

"""
    fourier_diff(sol,N,dL;format="matrix")

Compute the derivative using a Fourier differentiation matrix (default) or the spectral derivative for periodic functions for domain length `dL` and with `N` discretization points.

"""
function fourier_diff(sol,N,dL;format="matrix")
    if format == "matrix"
        h = 2*pi/N;
        col = vcat(0, 0.5*(-1).^(1:N-1).*cot.((1:N-1)*h/2));
        row = vcat(col[1], col[N:-1:2]);
        diff_matrix = Toeplitz(col,row);
        diff_sol = (2*pi/dL)*diff_matrix*sol; # Make dx calc abs...
    elseif format == "spectral"
        k = reduce(vcat,(2*π/dL)*[0:N/2-1 -N/2:-1]); # Wavenumbers
        sol_k = fft(sol);
        sol_k[Int(N/2)+1] = 0;
        diff_sol = real.(ifft(im*k.*sol_k));
    end
    return diff_sol
end

"""
    trapz(x_range,integrand)

Numerical integration using the multi-application trapezoidal rule.

"""
function trapz(x_range,integrand)
    int = 0.0;
    for i = 2:length(x_range)
        int += ((integrand[i-1] + integrand[i])/2)*(x_range[i] - x_range[i-1]);
    end
    return int
end

"""
    solution_spatial_sampling(x_prediction,x_target,solution)

Extract the solution values ``u(t,x)`` at a reduced number of equally spaced spatial locations.

"""
function solution_spatial_sampling(x_prediction,x_target,solution)
    ind = [findall(x_prediction[i].==x_target)[1] for i in 1:length(x_prediction)];
    solution_reduced = zeros(size(solution,1),size(ind,1));
    for i in 1:size(ind,1)
        solution_reduced[:,i] = solution[:,ind[i]]
    end
    return solution_reduced
end

"""
    solution_temporal_sampling(t_prediction,t_target,solution)

Extract the solution values ``u(t,x)`` at a reduced number of equally spaced temporal locations.

"""
function solution_temporal_sampling(t_prediction,t_target,solution)
    ind = [findall(t_prediction[i].==t_target)[1] for i in 1:length(t_prediction)];
    solution_reduced = zeros(size(ind,1),size(solution,2));
    for i in 1:size(ind,1)
        solution_reduced[i,:] = solution[ind[i],:]
    end
    return solution_reduced
end

"""
    fft_norm(solution)

Compute the FFT normalized by \$\\frac{1}{N}\$.

"""
function fft_norm(solution)
    N = size(solution,1);
    return (1/N)*fft(solution)
end

"""
    ifft_norm(solution)

Compute the IFFT normalized by \$N\$.

"""
function ifft_norm(solution)
    N = size(solution,1);
    return N*ifft(solution)
end

"""
    gauss_legendre(N,L1,L2)

Compute the nodes and weights for Gauss-Legendre quadrature with `N` discretization points and on the interval ``[L1,L2]``.

"""
function gauss_legendre(N,L1,L2)
    nodes, weights = gausslegendre(N);
    return ((L1+L2)/2).+((L2-L1)/2)*nodes, ((L2-L1)/2)*weights;
end

"""
    clenshaw_curtis(N,L1,L2)

Compute the nodes and weights for Clenshaw-Curtis quadrature with `N` discretization points and on the interval ``[L1,L2]``.

"""
function clenshaw_curtis(N,L1,L2)
    theta = pi*(0:N)/N;
    x = -cos.(theta);
    w = zeros(size(x,1));
    ii = 2:N;
    v = ones(N-1);
    if mod(N,2) == 0
        w[1] = 1/(N^2-1);
        w[N+1] = w[1];
        for k in 1:Int(N/2)-1
            v = v - 2*cos.(2*k*theta[ii])/(4*k^2-1);
        end
        v = v - cos.(N*theta[ii])/(N^2-1);
    else
        w[1] = 1/(N^2);
        w[N+1] = w[1];
        for k in 1:Int((N-1)/2)
            v = v - 2*cos.(2*k*theta[ii])/(4*k^2-1);
        end
    end
    w[ii] = 2*v/N;
    return ((L1+L2)/2).+((L2-L1)/2)*x, (L2-L1)/2*w
end

"""
    cheby_grid(N,L1,L2)

Generate the grid of Chebyshev points on the interval ``[L1,L2]`` with `N` discretization points.

"""
function cheby_grid(N,L1,L2)
    x = ((L1+L2)/2).+((L2-L1)/2)*(-cos.(pi*(0:N)/N));
    return x
end

"""
    cheby_diff_matrix(N,L1,L2)

Generate the Chebyshev differentiation matrix for the interval ``[L1,L2]`` with `N` discretization points.

"""
function cheby_diff_matrix(N,L1,L2)
    if N == 0
        D = 0;
        x = 1;
        return D, x
    else
        x = ((L1+L2)/2).+((L2-L1)/2)*(-cos.(pi*(0:N)/N));
        c = vcat(2, ones(N-1,1), 2).*(-1).^(0:N);
        X = repeat(x,1,N+1);
        dx = X-X';
        D = (c*(1 ./c)')./(dx+I);
        D = D - diagm(0 => sum(D,dims = 2)[:]);
        return D, x
    end
end

"""
    cheby_diff(sol,N,L1,L2)

Compute the derivative of using a Chebyshev differentiation matrix on the interval ``[L1,L2]`` with `N` discretization points.

"""
function cheby_diff(sol,N,L1,L2)
    D, x = cheby_diff_matrix(N,L1,L2);
    diff_sol = D*sol;
    return diff_sol
end

"""
    trapezoid(N,L1,L2)

Compute the nodes and weights for trapezoid rule with `N` discretization points and on the interval ``[L1,L2]``.

"""
function trapezoid(N,L1,L2)
    dL = abs(L2-L1);
    j = reduce(vcat,[0:1:N-1]);
    x = (dL.*j)./N;
    wi = ((L2-L1)/N);
    w = repeat([wi],N);
    return x, w
end

"""
    orthonormal_check(basis,weights;tol = 1e-12)

"""
function orthonormal_check(basis,weights;tol = 1e-12)
    W = diagm(0 => weights);
    for i = 1:size(basis,2)
        for j = 1:size(basis,2)
            if i == j
                if abs(sqrt(basis[:,i]'*W*basis[:,i]) - 1) > tol
                    error("Not orthonormal to $tol... $i")
                end
            elseif i != j
                if abs(basis[:,j]'*W*basis[:,i]) > tol
                    error("Not orthogonal to $tol... $i vs $j \n Review the singular value spectrum...")
                end
            end
        end
    end
end

"""
    legendre_norm(x,L1,L2,n)

"""
function legendre_norm(x,L1,L2,n)
    return ((2*n+1)/(L2-L1))^(1/2)*Pl(((2*x-L1-L2)/(L2-L1)),n)
end

"""
    legendre_norm_collect(x,L1,L2,nmax)

"""
function legendre_norm_collect(x,L1,L2,nmax)
    n_list = (0:1:nmax);
    return ((2*n_list.+1)/(L2-L1)).^(1/2).*parent(collectPl(((2*x-L1-L2)/(L2-L1)),lmax=nmax));
end


text:
using LeftChildRightSiblingTrees, AbstractTrees
using Test

@testset "LeftChildRightSiblingTrees" begin
    root = Node(0)
    @test isroot(root)
    @test isleaf(root)
    nchildren = 0
    for c in root
        nchildren += 1
    end
    @test nchildren == 0
    c1 = addchild(root, 1)
    c2 = addchild(root, 2)
    c3 = addsibling(c2, 3)
    @test lastsibling(c1) == c3
    c21 = addchild(c2, 4)
    c22 = addchild(c2, 5)
    @test isroot(root)
    @test !isleaf(root)
    nchildren = 0
    for c in root
        @test !isroot(c)
        nchildren += 1
    end
    @test nchildren == 3
    @test isleaf(c1)
    @test !isleaf(c2)
    @test isleaf(c3)
    for c in c2
        @test !isroot(c)
        @test isleaf(c)
    end
    children2 = [c21,c22]
    i = 0
    for c in c2
        @test c == children2[i+=1]
    end
    io = IOBuffer()
    show(io, c2)
    str = String(take!(io))
    @test str == "Node(2)"
    LeftChildRightSiblingTrees.showedges(io, root)
    str = String(take!(io))
    @test occursin("2 has the following children", str)

    @test depth(root) == 3
    @test depth(c3) == 1

    root1 = deepcopy(root)
    node = collect(root1)[2]
    graftchildren!(root1, node)
    @test isleaf(node)
    @test [c.data for c in root1] == [1,2,3,4,5]
    for c in root1
        @test c.parent == root1
    end
    prunebranch!(node)
    @test [c.data for c in root1] == [1,3,4,5]

    root1 = deepcopy(root)
    chlds = collect(root1)
    p, node = chlds[1], chlds[2]
    @test isleaf(p)
    graftchildren!(p, node)
    @test isleaf(node)
    @test [c.data for c in root1] == [1,2,3]
    @test [c.data for c in p] == [4,5]
    for c in p
        @test c.parent == p
    end

    root1 = deepcopy(root)
    chlds = collect(root1)
    prunebranch!(chlds[end])
    @test [c.data for c in root1] == [1,2]

    root1 = deepcopy(root)
    chlds = collect(root1)
    prunebranch!(chlds[1])
    @test [c.data for c in root1] == [2,3]
    @test_throws ErrorException("cannot prune the root") prunebranch!(root1)
end

@testset "AbstractTrees" begin
    root = Node(0)
    c1 = addchild(root, 1)
    c2 = addchild(root, 2)
    c3 = addsibling(c2, 3)
    c21 = addchild(c2, 4)
    c22 = addchild(c2, 5)
    io = IOBuffer()
    print_tree(io, root)
    @test strip(String(take!(io))) == "0\n├─ 1\n├─ 2\n│  ├─ 4\n│  └─ 5\n└─ 3"

    @test map(x->x.data, @inferred(collect(PostOrderDFS(root)))) == [1,4,5,2,3,0]
    @test map(x->x.data, @inferred(collect(PreOrderDFS(root)))) == [0,1,2,4,5,3]
    @test map(x->x.data, @inferred(collect(Leaves(root)))) == [1,4,5,3]
end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ b99805ba-3931-11eb-0e18-6ffb5497630d
begin
	using DrWatson

	function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
	end
end;

# ╔═╡ d840fcce-3931-11eb-0d08-59ed4b60b2f9
begin
	@quickactivate

	using GameZero
	using Random
    using Serialization
	using PlaymatSimulator
	import PlaymatSimulator.Actors.Image
	import PlaymatSimulator.Actors.GIF

	pd = projectdir()
	US = deserialize("$pd/tmp/user_selection.jls")

	md"""
	## EDH GAME STATE
	Provide an initial game state function for the game engine to load upon launch by defining a `gs` Dict{Symbol,Any}(...) object below with all of the game assets.
	"""
end

# ╔═╡ ebd61216-4da3-11eb-2e77-5900d3980b8e
STAGE = game_include("$pd/Base/notebooks/game_stage.jl")

# ╔═╡ f1e8642e-4da3-11eb-26bc-234df79ed8a6
GS = game_include("$pd/Base/notebooks/game_settings.jl")

# ╔═╡ 2e97f2d2-4c1d-11eb-23d5-3fc9bcd4bc47
#const SCREEN_WIDTH = Int32(1080)

# ╔═╡ 93af0788-4c12-11eb-2c17-8ddd05dfc94d
#const SCREEN_BORDER = Int32(10)

# ╔═╡ 988ca2cc-4c12-11eb-284d-f51814cbf3a3
#const SCREEN_HEIGHT = Int32(1920)

# ╔═╡ 22163bfc-3932-11eb-37cc-0701ba61483e
gs = Dict{Symbol,Any}(
	:GAME_NAME => "EDH",
	:MOUSE_POS => Int32[0,0],
    :MOUSE_OFFSETS => [ Int32[0,0] ],
	:music => readdir("$pd/Base/music"),
	:ui => OrderedDict(
		:cursor_icon => Image("mouse_cursor_icon",
			load("$pd/Base/ui/icons/RavenmoreIconPack/64/swordWood.png")),
        :cursor => Image("mouse_cursor",
			load("$pd/Base/ui/zones/area_wht.png"), alpha=0, w=1, h=1),
		),
    :stage => STAGE,
	:resources => Dict{Symbol,Any}(
        :mana => Dict{Symbol,Int32}(
            :white=>0,
            :blue=>0,
            :black=>0,
            :red=>0,
            :green=>0,
            :colorless=>0,
            ),
        :energy=>0,
        :life=>40,  # EDH life adjustment
        ),
	:resource_spinners => OrderedDict{Symbol,Actor}(
        :life => Text("40:L ",
            "$pd/Base/fonts/OpenSans-Semibold.ttf",
            ),
        :white_mana => Text(" 0:W ",
            "$pd/Base/fonts/OpenSans-Semibold.ttf", font_color = [255,255,255,255]
            ),
        :blue_mana => Text(" 0:U ",
            "$pd/Base/fonts/OpenSans-Semibold.ttf", font_color = [0,0,255,255]
            ),
        :black_mana => Text(" 0:B ",
            "$pd/Base/fonts/OpenSans-Semibold.ttf", font_color = [0,0,0,255]
            ),
        :red_mana => Text(" 0:R ",
            "$pd/Base/fonts/OpenSans-Semibold.ttf", font_color = [255,0,0,255]
            ),
        :green_mana => Text(" 0:G ",
            "$pd/Base/fonts/OpenSans-Semibold.ttf", font_color = [0,255,0,255]
            ),
        :colorless_mana => Text(" 0:C ",
            "$pd/Base/fonts/OpenSans-Semibold.ttf", font_color = [125,125,125,255]
            ),
        :energy => Text(" 0:E ",
            "$pd/Base/fonts/OpenSans-Semibold.ttf",
            ),
        :poison => Text(" 0:P ",
            "$pd/Base/fonts/OpenSans-Semibold.ttf", font_color = [255,0,255,255]
            ),
        ),
    :sfx => OrderedDict(
        :shade_wht=>Image("shade_wht",
			load("$pd/Base/ui/zones/area_wht.png")),
        :shade_blk=>Image("shade_blk",
			load("$pd/Base/ui/zones/area_blk.png")),
		:sel_box => Image("sel_box",
			load("$pd/Base/ui/zones/area_wht.png"), alpha=35, w=0, h=0),
        ),
    :overlay => OrderedDict(
        :texts => Actor[],
        :effects => Actor[],
		:dice => Dice[],
		:counters => Counter[],
		:shades => Actor[],
        ),
    :zone => OrderedDict(
        "Command" => Card[],  # EDH-specific zone
        "Graveyard" => Card[],
        "Battlefield" => Card[],
        "Library" => Card[],
        "Hand" => Card[],
		#"Exile" => Card[],  # replaced with Command
		#"Stack" => Union{Card,Spell}[]  # not implemented
        ),
    :group => OrderedDict(
        :clickables => Actor[],
        :selected => Actor[],
        ),
	:stack => OrderedDict(),
	)

# ╔═╡ 4339c7e0-4c1f-11eb-1372-5b42b11fa88e
begin
	merge!(gs, GS)
	gs
end

# ╔═╡ Cell order:
# ╠═b99805ba-3931-11eb-0e18-6ffb5497630d
# ╠═d840fcce-3931-11eb-0d08-59ed4b60b2f9
# ╠═ebd61216-4da3-11eb-2e77-5900d3980b8e
# ╠═f1e8642e-4da3-11eb-26bc-234df79ed8a6
# ╠═2e97f2d2-4c1d-11eb-23d5-3fc9bcd4bc47
# ╠═93af0788-4c12-11eb-2c17-8ddd05dfc94d
# ╟─988ca2cc-4c12-11eb-284d-f51814cbf3a3
# ╠═22163bfc-3932-11eb-37cc-0701ba61483e
# ╠═4339c7e0-4c1f-11eb-1372-5b42b11fa88e


text:
using ARules
using StatsBase
using BenchmarkTools

const itemlist = randstr(100);

n = 100_000
const m = 25             # number of items in transactions
const mx_depth = 7
t = [sample(itemlist, m, replace = false) for _ in 1:n];

# @code_warntype _frequent(t, 1)
const unq2 = unique_items(t)
# @benchmark occ2 = occurrence(t, unq2)
@time itree2 = frequent_item_tree(t, unq2, round(Int, 0.01*n), mx_depth);

Profile.clear_malloc_data()
itree2 = frequent_item_tree(t, unq2, round(Int, 0.01*n), mx_depth)

# ARules.shownodes(itree2)


text:
__precompile__()

module GLMCopula

using Convex, LinearAlgebra, MathProgBase, Reexport
using LinearAlgebra: BlasReal, copytri!
@reexport using Ipopt
@reexport using NLopt

export GaussianCopulaVCObs, GaussianCopulaVCModel
export fit!, fitted, init_β!, loglikelihood!, standardize_res!
export update_res!, update_Σ!

export GaussianCopulaLMMObs, GaussianCopulaLMMModel

"""
GaussianCopulaVCObs
GaussianCopulaVCObs(y, X, V)

A realization of Gaussian copula variance component data.
"""
struct GaussianCopulaVCObs{T <: BlasReal}
    # data
    y::Vector{T}
    X::Matrix{T}
    V::Vector{Matrix{T}}
    # working arrays
    ∇β::Vector{T}   # gradient wrt β
    ∇τ::Vector{T}   # gradient wrt τ
    ∇Σ::Vector{T}   # gradient wrt σ2
    Hβ::Matrix{T}   # Hessian wrt β
    Hτ::Matrix{T}   # Hessian wrt τ
    res::Vector{T}  # residual vector res_i
    xtx::Matrix{T}  # Xi'Xi
    t::Vector{T}    # t[k] = tr(V_i[k]) / 2
    q::Vector{T}    # q[k] = res_i' * V_i[k] * res_i / 2
    storage_n::Vector{T}
    storage_p::Vector{T}
end

function GaussianCopulaVCObs(
    y::Vector{T},
    X::Matrix{T},
    V::Vector{Matrix{T}}
    ) where T <: BlasReal
    n, p, m = size(X, 1), size(X, 2), length(V)
    @assert length(y) == n "length(y) should be equal to size(X, 1)"
    # working arrays
    ∇β  = Vector{T}(undef, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇Σ  = Vector{T}(undef, m)
    Hβ  = Matrix{T}(undef, p, p)
    Hτ  = Matrix{T}(undef, 1, 1)
    res = Vector{T}(undef, n)
    xtx = transpose(X) * X
    t   = [tr(V[k])/2 for k in 1:m]
    q   = Vector{T}(undef, m)
    storage_n = Vector{T}(undef, n)
    storage_p = Vector{T}(undef, p)
    # constructor
    GaussianCopulaVCObs{T}(y, X, V, ∇β, ∇τ, ∇Σ, Hβ,
        Hτ, res, xtx, t, q, storage_n, storage_p)
end

"""
GaussianCopulaVCModel
GaussianCopulaVCModel(gcs)

Gaussian copula variance component model, which contains a vector of
`GaussianCopulaVCObs` as data, model parameters, and working arrays.
"""
struct GaussianCopulaVCModel{T <: BlasReal} <: MathProgBase.AbstractNLPEvaluator
    # data
    data::Vector{GaussianCopulaVCObs{T}}
    ntotal::Int     # total number of singleton observations
    p::Int          # number of mean parameters in linear regression
    m::Int          # number of variance components
    # parameters
    β::Vector{T}    # p-vector of mean regression coefficients
    τ::Vector{T}    # inverse of linear regression variance parameter
    Σ::Vector{T}    # m-vector: [σ12, ..., σm2]
    # working arrays
    ∇β::Vector{T}   # gradient from all observations
    ∇τ::Vector{T}
    ∇Σ::Vector{T}
    Hβ::Matrix{T}    # Hessian from all observations
    Hτ::Matrix{T}
    XtX::Matrix{T}  # X'X = sum_i Xi'Xi
    TR::Matrix{T}   # n-by-m matrix with tik = tr(Vi[k]) / 2
    QF::Matrix{T}   # n-by-m matrix with qik = res_i' Vi[k] res_i
    storage_n::Vector{T}
    storage_m::Vector{T}
    storage_Σ::Vector{T}
end

function GaussianCopulaVCModel(gcs::Vector{GaussianCopulaVCObs{T}}) where T <: BlasReal
    n, p, m = length(gcs), size(gcs[1].X, 2), length(gcs[1].V)
    β   = Vector{T}(undef, p)
    τ   = Vector{T}(undef, 1)
    Σ   = Vector{T}(undef, m)
    ∇β  = Vector{T}(undef, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇Σ  = Vector{T}(undef, m)
    Hβ  = Matrix{T}(undef, p, p)
    Hτ  = Matrix{T}(undef, 1, 1)
    XtX = zeros(T, p, p) # sum_i xi'xi
    TR  = Matrix{T}(undef, n, m) # collect trace terms
    ntotal = 0
    for i in eachindex(gcs)
        ntotal  += length(gcs[i].y)
        XtX    .+= gcs[i].xtx
        TR[i, :] = gcs[i].t
    end
    QF        = Matrix{T}(undef, n, m)
    storage_n = Vector{T}(undef, n)
    storage_m = Vector{T}(undef, m)
    storage_Σ = Vector{T}(undef, m)
    GaussianCopulaVCModel{T}(gcs, ntotal, p, m, β, τ, Σ,
        ∇β, ∇τ, ∇Σ, Hβ, Hτ, XtX, TR, QF,
        storage_n, storage_m, storage_Σ)
end

"""
GaussianCopulaLMMObs
GaussianCopulaLMMObs(y, X, Z)

A realization of Gaussian copula linear mixed model data instance.
"""
struct GaussianCopulaLMMObs{T <: LinearAlgebra.BlasReal}
    # data
    y::Vector{T}
    X::Matrix{T}
    Z::Matrix{T}
    # working arrays
    ∇β::Vector{T}   # gradient wrt β
    ∇τ::Vector{T}   # gradient wrt τ
    ∇Σ::Matrix{T}   # gradient wrt Σ
    Hβ::Matrix{T}   # Hessian wrt β
    Hτ::Matrix{T}   # Hessian wrt τ
    HΣ::Matrix{T}   # Hessian wrt Σ
    res::Vector{T}  # residual vector
    xtx::Matrix{T}  # Xi'Xi (p-by-p)
    ztz::Matrix{T}  # Zi'Zi (q-by-q)
    xtz::Matrix{T}  # Xi'Zi (p-by-q)
    storage_q1::Vector{T}
    storage_q2::Vector{T}
end

function GaussianCopulaLMMObs(
    y::Vector{T},
    X::Matrix{T},
    Z::Matrix{T}
    ) where T <: BlasReal
    n, p, q = size(X, 1), size(X, 2), size(Z, 2)
    @assert length(y) == n "length(y) should be equal to size(X, 1)"
    # working arrays
    ∇β  = Vector{T}(undef, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇Σ  = Matrix{T}(undef, q, q)
    Hβ  = Matrix{T}(undef, p, p)
    Hτ  = Matrix{T}(undef, 1, 1)
    HΣ  = Matrix{T}(undef, abs2(q), abs2(q))
    res = Vector{T}(undef, n)
    xtx = transpose(X) * X
    ztz = transpose(Z) * Z
    xtz = transpose(X) * Z
    storage_q1 = Vector{T}(undef, q)
    storage_q2 = Vector{T}(undef, q)
    # constructor
    GaussianCopulaLMMObs{T}(y, X, Z,
        ∇β, ∇τ, ∇Σ, Hβ, Hτ, HΣ,
        res, xtx, ztz, xtz,
        storage_q1, storage_q2)
end

"""
GaussianCopulaLMMModel
GaussianCopulaLMMModel(gcs)

Gaussian copula linear mixed model, which contains a vector of
`GaussianCopulaLMMObs` as data, model parameters, and working arrays.
"""
struct GaussianCopulaLMMModel{T <: BlasReal} <: MathProgBase.AbstractNLPEvaluator
    # data
    data::Vector{GaussianCopulaLMMObs{T}}
    ntotal::Int     # total number of singleton observations
    p::Int          # number of mean parameters in linear regression
    q::Int          # number of random effects
    # parameters
    β::Vector{T}    # p-vector of mean regression coefficients
    τ::Vector{T}    # inverse of linear regression variance parameter
    Σ::Matrix{T}    # q-by-q (psd) matrix
    # working arrays
    ΣL::Matrix{T}
    ∇β::Vector{T}   # gradient from all observations
    ∇τ::Vector{T}
    ∇Σ::Matrix{T}
    Hβ::Matrix{T}   # Hessian from all observations
    Hτ::Matrix{T}
    HΣ::Matrix{T}
    XtX::Matrix{T}      # X'X = sum_i Xi'Xi
    storage_qq::Matrix{T}
    storage_nq::Matrix{T}
end

function GaussianCopulaLMMModel(gcs::Vector{GaussianCopulaLMMObs{T}}) where T <: BlasReal
    n, p, q = length(gcs), size(gcs[1].X, 2), size(gcs[1].Z, 2)
    npar = p + 1 + (q * (q + 1)) >> 1
    β   = Vector{T}(undef, p)
    τ   = Vector{T}(undef, 1)
    Σ   = Matrix{T}(undef, q, q)
    ΣL  = similar(Σ)
    ∇β  = Vector{T}(undef, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇Σ  = Matrix{T}(undef, q, q)
    Hβ  = Matrix{T}(undef, p, p)
    Hτ  = Matrix{T}(undef, 1, 1)
    HΣ  = Matrix{T}(undef, abs2(q), abs2(q))
    XtX = zeros(T, p, p) # sum_i xi'xi
    ntotal = 0
    for i in eachindex(gcs)
        ntotal  += length(gcs[i].y)
        XtX    .+= gcs[i].xtx
    end
    storage_qq = Matrix{T}(undef, q, q)
    storage_nq = Matrix{T}(undef, n, q)
    GaussianCopulaLMMModel{T}(gcs, ntotal, p, q,
        β, τ, Σ, ΣL,
        ∇β, ∇τ, ∇Σ, Hβ, Hτ, HΣ,
        XtX, storage_qq, storage_nq)
end

include("gaussian_vc.jl")
include("gaussian_lmm.jl")
include("rand.jl")

end#module


text:
# Based in DrWatson

const DATA_KEY = :dat
const GIT_COMMIT_KEY = :gitcommit
const GIT_PATCH_KEY = :gitpatch
const DIRTY_SUFFIX = "_dirty"
const SHORT_HASH_LENGTH = 7
const ERROR_LOGGER = SimpleLogger(stdout, Logging.Error)

# ------------------------------------------------------------------
function ldat(dfargs...;
        print_fun::Function = global_conf(:PRINT_FUN),
        load_fun::Function = global_conf(:LOAD_FUN),
        mkdir::Bool = global_conf(:MK_DIR),
        verbose::Bool = global_conf(:VERBOSE),
        msg::AbstractString = "",
    )

    datfile = dfname(dfargs...)
    mkdir && mkpath(dirname(datfile))
    file_dat = load_fun(datfile)
    dat = file_dat[DATA_KEY]
    if verbose
        commit_hash = get(file_dat, GIT_COMMIT_KEY, "none")
        _io_print(
            print_fun, "DATA LOADED", msg, dat, datfile,
            "\ncommit: ", _cut_hash(commit_hash)
        )
    end
    return dat
end

function ldat(f::Function, dfargs...;
        load_fun::Function = global_conf(:LOAD_FUN),
        addtag::Bool = global_conf(:ADD_TAG),
        save_fun::Function = global_conf(:SAVE_FUN),
        kwargs...
    )
    datfile = dfname(dfargs...)
    !isfile(datfile) && sdat(f(), datfile; save_fun, addtag, kwargs...)
    return ldat(datfile; load_fun, kwargs...)
end

# ------------------------------------------------------------------
function sdat(f::Function, dfargs...;
        print_fun::Function = global_conf(:PRINT_FUN),
        addtag::Bool = global_conf(:ADD_TAG),
        save_fun::Function = global_conf(:SAVE_FUN),
        mkdir::Bool = global_conf(:MK_DIR),
        verbose::Bool = global_conf(:VERBOSE),
        msg::AbstractString = "",
    )
    dat = f()
    datfile = dfname(dfargs...)
    mkdir && mkpath(dirname(datfile))

    L = verbose ? global_logger() : ERROR_LOGGER
    with_logger(L) do
        dict = Dict(DATA_KEY => dat)
        tagdat = addtag ? DrWatson.tag!(dict) : dict
        save_fun(datfile, tagdat)
        commit_hash = get(tagdat, GIT_COMMIT_KEY, "none")
        verbose && verbose && _io_print(
            print_fun, "DATA SAVED", msg, dat, datfile,
            "\ncommit: ", _cut_hash(commit_hash)
        )
        return datfile
    end
end
sdat(dat, dfargs...; kwargs...) = sdat(() -> dat, dfargs...; kwargs...)

# ------------------------------------------------------------------
function dhash(datfile, l = SHORT_HASH_LENGTH)
    hash = get(_load(datfile), GIT_COMMIT_KEY, "")
    _cut_hash(hash, l)
end

function _cut_hash(commit_hash, l = SHORT_HASH_LENGTH)
    short_hash = first(commit_hash, l)
    endswith(commit_hash, DIRTY_SUFFIX) ? string(short_hash, DIRTY_SUFFIX) : short_hash
end

lgitpatth(args...) = get(ldat(args...), GIT_PATCH_KEY, "")


text:
using BinDeps
using Compat

@BinDeps.setup

deps = [
	omp = library_dependency("libgomp", os= :Unix)
]

run(`make`)

max_stars_repo_name: camruggles/testGenerate

max_issues_repo_name: camruggles/testGenerate

max_forks_repo_name: camruggles/testGenerate


text:
module ShuHua

project(xs...) = normpath(joinpath(dirname(dirname(pathof(ShuHua))), xs...))

include("types.jl")
include("parse.jl")
include("write.jl")
include("build.jl")
include("serve.jl")

using .Types
using .Parse
using .Write

export parse_slide

end


text:
struct Bond

   i :: Int64
   j :: Int64

   kb :: Float64
   b0 :: Float64

end


text:
using JLD, Knet
using AutoGrad: cat1d
include("types.jl")
include("lm_utils.jl")
include("preprocess.jl")
include("parser.jl")
include("helper.jl")
include("inits.jl")
#include("old_features.jl")
# Here you have seen the developments through the model
#include("old_model.jl")
#include("new_model.jl")
#include("new_model2.jl")
#include("new_model3.jl") # include model file in main train-file

const MAXSENT=64       # skip longer sentences during training
const MINSENT=2        # skip shorter sentences during training


text:
struct LazyNode{Name,D} <: AbstractNode
	data::D
end

const LazyNode{Name} = LazyNode{Name, D} where {D}
LazyNode(Name::Symbol, values::T) where {T} = LazyNode{Name,T}(values)
LazyNode{Name}(values::T) where {Name, T} = LazyNode{Name,T}(values)

Base.ndims(x::LazyNode) = Colon()
LearnBase.nobs(a::LazyNode) = length(a.data)
LearnBase.nobs(a::LazyNode, ::Type{ObsDim.Last}) = nobs(a.data)

const VectorOfLazy{N} = Vector{T} where {T<:LazyNode{N,M} where {M}}
function Base.reduce(::typeof(Mill.catobs), as::VectorOfLazy{N}) where {N}
    data = reduce(vcat, [x.data for x in as])
    LazyNode{N}(data)
end

Base.getindex(x::LazyNode{N,T}, i::VecOrRange) where {N,T}= LazyNode{N}(subset(x.data, i))

Base.hash(e::LazyNode{N,D}, h::UInt) where {N,D} = hash((N, D, e.data), h)
Base.:(==)(e1::LazyNode{N,D}, e2::LazyNode{N,D}) where {N,D} = e1.data == e2.data


text:
using SLICOTMath
const SLICOT=SLICOTMath
using Test
using LinearAlgebra
using LinearAlgebra: BlasInt
const BlasBool = BlasInt

LSAME(x,y) = uppercase(x) == uppercase(y)

# helpers for getting data from the examples
parsex(T,s) = parse(T,s)
# CAUTION: this relies on absence of spaces in complex constants
function parsex(::Type{ComplexF64}, s)
  if s == "0"
    return zero(ComplexF64)
  end
  s1 = split(s,",")
  if length(s1) == 1
    return parse(Float64,s1[1])+0.0im
  end
  if s1[2][1] in ['+','-']
    impart = s1[2][1:end-1] * "im"
  else
    impart = '+' * s1[2][1:end-1] * "im"
  end
  s2 = s1[1][2:end] * impart
  parse(ComplexF64, s2)
end
function parsex(::Type{Bool}, s)
    return s[2] in ('T','t')
end

const quiet = Ref(true)
const tdir = joinpath(@__DIR__,"SLICOT-Reference")

if !isdir(joinpath(tdir,"examples"))
    @error """Package tests require that 'test/SLICOT-Reference' hold the input data.
This can be prepared by installing a copy of the SLICOT distribution there,
or making a link to one.
"""
    # force an early exit
    @testset "installation" begin
        @test false
    end
end

files = readdir(@__DIR__)
funclist = [replace(s, "_test.jl" => "") for s in files if (length(s)>7) && (s[end-7:end] == "_test.jl")]

for func2test in funclist
    testcode = lowercase(func2test) * "_test.jl"
    include(testcode)

    rname=Symbol("run_"*func2test)
    eval(:(runme=$rname))

    @testset "$func2test" begin
        datfile = joinpath(tdir, "examples", uppercase(func2test) * ".dat")

        if quiet[]
            runme(datfile, devnull)
        else
            runme(datfile)
        end
    end # testset
end # for


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Parameter <: CGE
        parameter::Symbol
        name::String
        index::Array{Symbol,1}
        units::String
    end

Information about CGE parameters used in the model.

# Arguments
- `parameter::Symbol`: parameter variable abbreviation
- `name::String`: parameter description
- `index::Array{Symbol,1}`: Sets on which the parameter depends
- `units::String`: parameter units
"""
mutable struct Parameter <: CGE
    "parameter variable abbreviation"
    parameter::Symbol
    "parameter description"
    name::String
    "Sets on which the parameter depends"
    index::Array{Symbol,1}
    "parameter units"
    units::String
end

function Parameter(; parameter, name, index, units, )
    Parameter(parameter, name, index, units, )
end

"""Get [`Parameter`](@ref) `parameter`."""
get_parameter(value::Parameter) = value.parameter
"""Get [`Parameter`](@ref) `name`."""
get_name(value::Parameter) = value.name
"""Get [`Parameter`](@ref) `index`."""
get_index(value::Parameter) = value.index
"""Get [`Parameter`](@ref) `units`."""
get_units(value::Parameter) = value.units

"""Set [`Parameter`](@ref) `parameter`."""
set_parameter!(value::Parameter, val) = value.parameter = val
"""Set [`Parameter`](@ref) `name`."""
set_name!(value::Parameter, val) = value.name = val
"""Set [`Parameter`](@ref) `index`."""
set_index!(value::Parameter, val) = value.index = val
"""Set [`Parameter`](@ref) `units`."""
set_units!(value::Parameter, val) = value.units = val

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
export VlasovProblem

abstract type AbstractProblem end

"""
    VlasovProblem( f, method, dev)

"""
struct VlasovProblem{Method<:AbstractMethod} <: AbstractProblem

    f :: DistributionFunction
    method :: Method
    dev :: AbstractDevice

    function VlasovProblem( f, method:: BSLSpline, dev)

        new{BSLSpline}( f, method, dev)

    end

    function VlasovProblem( f, method::Fourier, dev)

        new{Fourier}( f, method, dev)

    end

end

export solve!

"""
    solve!( problem, stepper, dt, nsteps)

"""
function solve!( problem::VlasovProblem{BSLSpline}, stepper::StrangSplitting, dt, nsteps)

  nrj = Vector{Float64}(undef, nsteps)

  for it in 1:nsteps

     advection_x!( problem.f, 0.5dt, problem.method.p)
     sol = advection_v!( problem.f, dt, problem.method.p)
     nrj[it] = sol
     advection_x!( problem.f, 0.5dt, problem.method.p)

  end

  nrj

end

function solve!( problem::VlasovProblem{Fourier}, stepper::StrangSplitting, dt, nsteps)

    # Initialize distribution function
    x = problem.f.xgrid.points |> collect
    v = problem.f.vgrid.points |> collect
    nx = problem.f.xgrid.len
    nv = problem.f.vgrid.len

    f = zeros(Complex{Float64},(nx,nv))
    fᵀ= zeros(Complex{Float64},(nv,nx))

    f .= problem.f.values

    transpose!(fᵀ,f)

    dx  = problem.f.xgrid.step
    dv  = problem.f.vgrid.step
    ρ   = dv .* vec(sum(real(fᵀ), dims=1))
    ρ  .-= mean(ρ)
    e   = zeros(ComplexF64, nx)
    nx  = problem.f.xgrid.len
    Lx  = problem.f.xgrid.stop - problem.f.xgrid.start
    modes  = zeros(Float64, nx)
    modes .= 2π / Lx .* vcat(0:nx÷2-1,-nx÷2:-1)
    modes[1] = 1.0
    ρ̂ = fft(ρ)./modes
    e .= vec(real(ifft(-1im .* ρ̂)))

    nrj = Vector{Float64}(undef, nsteps)

    for it in 1:nsteps

        advection_v!(fᵀ, problem.method, e, 0.5dt)
        transpose!(f,fᵀ)
        advection_x!( f, problem.method, e, v, dt)
        nrj[it] = log(sqrt((sum(e.^2)) * dx))
        transpose!(fᵀ,f)
        advection_v!(fᵀ, problem.method, e, 0.5dt)

    end

    nrj

end


text:
<extend src="layout.jl">
    <block name="content">
        <include src="post-item.jl"/>
        <post-item class="post-info"/>

            {{ post.text }}

        <form method="post" action="/post/" + post.id + "/comment" style="padding-left: 1em;">
            <textarea id="text" name="text" rows="6" cols="60"></textarea>

            <input type="submit" value="add comment">
        </form>

    </block>
</extend>


text:
# Instantiates CR3BPIndirectParams struct for a given scenario
function initCR3BPIndirectParams(scenario::String)
    if scenario == "Low Thrust 10 CR3BP"
        sp = SimpleSpacecraft(1500.0, 1000.0, 10.0, 3000.0)
        cr3 = CR3BPParams(5.9742e24, 7.3483e22,
                          6378.137, 1738.0,
                          1500.0, 384400)
        ps = CR3BPIndirectParams(sp, cr3, 1.0)
    elseif scenario ==  "Low Thrust 1 CR3BP"
        sp = SimpleSpacecraft(1500.0, 1000.0, 1.0, 3000.0)
        cr3 = CR3BPParams(5.9742e24, 7.3483e22,
                          6378.137, 1738.0,
                          1500.0, 384400)
        ps = CR3BPIndirectParams(sp, cr3, 1.0)
    else
        error("Scenario is not implemented!")
    end

    return ps
end

# Instantiates CR3BPIndirectWithSTM struct for a given scenario
function initCR3BPIndirectWithSTMParams(scenario::String)
    if scenario == "Low Thrust 10 CR3BP"
        sp = SimpleSpacecraft(1500.0, 1000.0, 10.0, 3000.0)
        cr3 = CR3BPParams(5.9742e24, 7.3483e22,
                          6378.137, 1738.0,
                          1500.0, 384400)
        ps = CR3BPIndirectWithSTMParams(sp, cr3, 1.0)
    elseif scenario == "Low Thrust 1 CR3BP"
        sp = SimpleSpacecraft(1500.0, 1000.0, 1.0, 3000.0)
        cr3 = CR3BPParams(5.9742e24, 7.3483e22,
                          6378.137, 1738.0,
                          1500.0, 384400)
        ps = CR3BPIndirectWithSTMParams(sp, cr3, 1.0)
    else
        error("Scenario is not implemented!")
    end

    return ps
end


text:
using Test
using Random
using LinearAlgebra
using SparseArrays
using JLD2, FileIO
using TraceEstimation

# Diagonal Approximation method as described by P. Fika works only for SPD matrix with low condition number (< 500)
@testset "Diagonal Approximation (P.Fika)" begin
    Random.seed!(1234323)
    @testset "Dense SPD Hermitian Matrices" begin
        @testset "a[i,j] = exp(-2 * abS(i - j)) (small size)" begin
            println("Executing Test 01: Dense SPD Small Size")
            A = rand(610, 610)
            for i in 1:610
                for j in 1:610
                    A[i, j] = exp(-2 * abs(i - j))
                end
            end
            obv = diagapp(A)
            acv = tr(inv(A))
            @test isapprox(obv, acv, rtol=10)
        end
        @testset "a[i,j] = exp(-2 * abS(i - j)) (large size)" begin
            println("Executing Test 02: Dense SPD Large Size")
            A = rand(5100, 5100)
            for i in 1:5100
                for j in 1:5100
                    A[i, j] = exp(-2 * abs(i - j))
                end
            end
            obv = diagapp(A)
            acv = tr(inv(A))
            @test isapprox(obv, acv, rtol=10)
        end
        @testset "Random generated SPD matrix (small size) (N=30)" begin
            println("Executing Test 03: Random SPD Small Size (N = 30)")
            A = rand(810, 810)
            A = A + A' + 30I
            while isposdef(A) == false
                A = rand(810, 810)
                A = A + A' + 30I
            end
            obv = diagapp(A)
            acv = tr(inv(A))
            @test isapprox(obv, acv, rtol=10)
        end
        @testset "Random generated SPD matrix (small size) (N=60)" begin
            println("Executing Test 04: Random SPD Small Size (N=60)")
            A = rand(810, 810)
            A = A + A' + 30I
            while isposdef(A) == false
                A = rand(810, 810)
                A = A + A' + 30I
            end
            obv = diagapp(A)
            acv = tr(inv(A))
            @test isapprox(obv, acv, rtol=10)
        end
        @testset "Random generated SPD matrix (large size) (N=30)" begin
            println("Executing Test 05: Random SPD Large Size")
            A = rand(5100, 5100)
            A = A + A' + 300I
            while isposdef(A) == false
                A = rand(5100, 5100)
                A = A + A' + 300I
            end
            obv = diagapp(A)
            acv = tr(inv(A))
            @test isapprox(obv, acv, rtol=10)
        end
    end
    @testset "Random generated SPD matrix (large size, large condition number)" begin
            println("Executing Test 05.5: Random SPD Large Size Larger Condition")
            A = Symmetric(rand(5005, 5005))
            A = A + 10I
            while isposdef(A) == false
                A = A + 10I
            end
            obv = diagapp(A)
            acv = tr(inv(A))
            @test isapprox(obv, acv, rtol=10)
        end
    @testset "Sparse SPD Hermitian Matrices" begin
        @testset "Random generated Sparse SPD (small size)" begin
            println("Executing Test 06: Sparse Random SPD Small Size")
            A = Symmetric(sprand(1000, 1000, 0.7))
            A = A+50*I
            while isposdef(A) == false
                A = Symmetric(sprand(1000, 1000, 0.7))
                A = A+50*I
            end
            obv = diagapp(A)
            A = Matrix(A)
            acv = tr(inv(A))
            @test isapprox(obv, acv, rtol=10)
        end
        #= Condition number too large (> 1000)
        @testset "TopOpt Test (Large Condition Number)" begin
            println("Executing Test 07: TopOpt Large Condition Number")
            s = (40, 10) # increase to increase the matrix size
            xmin = 0.9 # decrease to increase the condition number
            problem = HalfMBB(Val{:Linear}, s, (1.0, 1.0), 1.0, 0.3, 1.0)
            solver = FEASolver(Displacement, Direct, problem, xmin = xmin)
            n = length(solver.vars)
            solver.vars[rand(1:n, n÷2)] .= 0
            solver()
            K = solver.globalinfo.K
            obv = diagapp(K)
            M = Matrix(K)
            acv = tr(inv(M))
            @test isapprox(obv, acv, rtol=10)
        end
        =#
        @testset "TopOpt Test (Modified Condition Number)" begin
            println("Executing Test 08: TopOpt Modified Condition Number")
            println("Loading JLD2 file")
            #=
            s = (40, 10) # increase to increase the matrix size
            xmin = 0.9 # decrease to increase the condition number
            problem = HalfMBB(Val{:Linear}, s, (1.0, 1.0), 1.0, 0.3, 1.0)
            solver = FEASolver(Displacement, Direct, problem, xmin = xmin)
            n = length(solver.vars)
            solver.vars[rand(1:n, n÷2)] .= 0
            solver()
            K = solver.globalinfo.K
            K = K+1*I
            =#
            @load "topoptfile.jld2" K
            obv = diagapp(K)
            M = Matrix(K)
            acv = tr(inv(M))
            @test isapprox(obv, acv, rtol=10)
        end
    end
    #=
    @testset "CuArrays Test" begin
        @testset "Dense Random CuArray" begin
            println("Executing Test 09: CuArray Small Size")
            A = cu(rand(400,400))
            A = A+A'+40*I
            function mycurand(type::Type, size::Int64)
                cu(rand(type, size))
            end
            function mycurand(range::StepRange{<:Any, <:Any}, size::Int64)
                cu(rand(range, size))
            end
            obv = diagapp(A, randfunc=mycurand)
            M = Matrix(A)
            acv = tr(inv(M))
            @test isapprox(obv, acv, rtol=10)
        end
    end =#
end


text:
const DevicesModelContainer = Dict{Symbol, DeviceModel}
const BranchModelContainer = Dict{Symbol, DeviceModelForBranches}
const ServicesModelContainer = Dict{Tuple{String, Symbol}, ServiceModel}

"""
    ProblemTemplate(::Type{T}) where {T<:PM.AbstractPowerFormulation}

Creates a model reference of the PowerSimulations Optimization Problem.

# Arguments

  - `model::Type{T<:PM.AbstractPowerFormulation}`:

# Example

template = ProblemTemplate(CopperPlatePowerModel)
"""
mutable struct ProblemTemplate
    network_model::NetworkModel{<:PM.AbstractPowerModel}
    devices::DevicesModelContainer
    branches::BranchModelContainer
    services::ServicesModelContainer
    function ProblemTemplate(network::NetworkModel{T}) where {T <: PM.AbstractPowerModel}
        new(
            network,
            DevicesModelContainer(),
            BranchModelContainer(),
            ServicesModelContainer(),
        )
    end
end

ProblemTemplate(::Type{T}) where {T <: PM.AbstractPowerModel} =
    ProblemTemplate(NetworkModel(T))
ProblemTemplate() = ProblemTemplate(CopperPlatePowerModel)

get_device_models(template::ProblemTemplate) = template.devices
get_branch_models(template::ProblemTemplate) = template.branches
get_service_models(template::ProblemTemplate) = template.services
get_network_model(template::ProblemTemplate) = template.network_model
get_network_formulation(template::ProblemTemplate) =
    get_network_formulation(get_network_model(template))

function get_model(template::ProblemTemplate, device_type)
    if device_type <: PSY.Device
        return get(template.devices, Symbol(device_type), nothing)
    elseif device_type <: PSY.Branch
        return get(template.branches, Symbol(device_type), nothing)
    elseif device_type <: PSY.Service
        return get(template.services, Symbol(device_type), nothing)
    else
        error("not supported: $device_type")
    end
end

# Note to devs. PSY exports set_model! these names are chosen to avoid name clashes

"""
Sets the network model in a template.
"""
function set_network_model!(
    template::ProblemTemplate,
    model::NetworkModel{<:PM.AbstractPowerModel},
)
    template.network_model = model
    return
end

"""
Sets the device model in a template using the component type and formulation.
Builds a default DeviceModel
"""
function set_device_model!(
    template::ProblemTemplate,
    component_type::Type{<:PSY.Device},
    formulation::Type{<:AbstractDeviceFormulation},
)
    set_device_model!(template, DeviceModel(component_type, formulation))
    return
end

"""
Sets the device model in a template using a DeviceModel instance
"""
function set_device_model!(
    template::ProblemTemplate,
    model::DeviceModel{<:PSY.Device, <:AbstractDeviceFormulation},
)
    _set_model!(template.devices, model)
    return
end

function set_device_model!(
    template::ProblemTemplate,
    model::DeviceModel{<:PSY.Branch, <:AbstractDeviceFormulation},
)
    _set_model!(template.branches, model)
    return
end

"""
Sets the service model in a template using a name and the service type and formulation.
Builds a default ServiceModel with use_service_name set to true.
"""
function set_service_model!(
    template::ProblemTemplate,
    service_name::String,
    service_type::Type{<:PSY.Service},
    formulation::Type{<:AbstractServiceFormulation},
)
    set_service_model!(
        template,
        service_name,
        ServiceModel(service_type, formulation, use_service_name=true),
    )
    return
end

"""
Sets the service model in a template using a ServiceModel instance.
"""
function set_service_model!(
    template::ProblemTemplate,
    service_type::Type{<:PSY.Service},
    formulation::Type{<:AbstractServiceFormulation},
)
    set_service_model!(template, ServiceModel(service_type, formulation))
    return
end

function set_service_model!(
    template::ProblemTemplate,
    service_name::String,
    model::ServiceModel{<:PSY.Service, <:AbstractServiceFormulation},
)
    _set_model!(template.services, service_name, model)
    return
end

function set_service_model!(
    template::ProblemTemplate,
    model::ServiceModel{<:PSY.Service, <:AbstractServiceFormulation},
)
    _set_model!(template.services, model)
    return
end


text:
module TestDeprecated

using Test, DataFrames

@testset "by and aggregate" begin
    @test_throws ArgumentError by()
    @test_throws ArgumentError aggregate()
end

@testset "All indexing" begin
    df = DataFrame(a=1, b=2, c=3)

    @test select(df, All(1, 2)) == df[:, 1:2]
    @test select(df, All(1, :b)) == df[:, 1:2]
    @test select(df, All(:a, 2)) == df[:, 1:2]
    @test select(df, All(:a, :b)) == df[:, 1:2]
    @test select(df, All(2, 1)) == df[:, [2, 1]]
    @test select(df, All(:b, 1)) == df[:, [2, 1]]
    @test select(df, All(2, :a)) == df[:, [2, 1]]
    @test select(df, All(:b, :a)) == df[:, [2, 1]]

    @test df[:, All(1, 2)] == df[:, 1:2]
    @test df[:, All(1, :b)] == df[:, 1:2]
    @test df[:, All(:a, 2)] == df[:, 1:2]
    @test df[:, All(:a, :b)] == df[:, 1:2]
    @test df[:, All(2, 1)] == df[:, [2, 1]]
    @test df[:, All(:b, 1)] == df[:, [2, 1]]
    @test df[:, All(2, :a)] == df[:, [2, 1]]
    @test df[:, All(:b, :a)] == df[:, [2, 1]]

    @test df[:, All(1, 1, 2)] == df[:, 1:2]
    @test df[:, All(:a, 1, :b)] == df[:, 1:2]
    @test df[:, All(:a, 2, :b)] == df[:, 1:2]
    @test df[:, All(:a, :b, 2)] == df[:, 1:2]
    @test df[:, All(2, 1, :a)] == df[:, [2, 1]]

    @test select(df, All(1, "b")) == df[:, 1:2]
    @test select(df, All("a", 2)) == df[:, 1:2]
    @test select(df, All("a", "b")) == df[:, 1:2]
    @test select(df, All("b", 1)) == df[:, [2, 1]]
    @test select(df, All(2, "a")) == df[:, [2, 1]]
    @test select(df, All("b", "a")) == df[:, [2, 1]]

    @test df[:, All(1, "b")] == df[:, 1:2]
    @test df[:, All("a", 2)] == df[:, 1:2]
    @test df[:, All("a", "b")] == df[:, 1:2]
    @test df[:, All("b", 1)] == df[:, [2, 1]]
    @test df[:, All(2, "a")] == df[:, [2, 1]]
    @test df[:, All("b", "a")] == df[:, [2, 1]]

    @test df[:, All("a", 1, "b")] == df[:, 1:2]
    @test df[:, All("a", 2, "b")] == df[:, 1:2]
    @test df[:, All("a", "b", 2)] == df[:, 1:2]
    @test df[:, All(2, 1, "a")] == df[:, [2, 1]]

    df = DataFrame(a1=1, a2=2, b1=3, b2=4)
    @test df[:, All(r"a", Not(r"1"))] == df[:, [1, 2, 4]]
    @test df[:, All(Not(r"1"), r"a")] == df[:, [2, 4, 1]]
end

end # module


text:
#################################################################
#
#    Main file of MCMC.jl module
#
#################################################################

module MCMC

using DataFrames, DataArrays
using Distributions
using Stats

import Base.*, Base.show
export show, *
export MCMCTask, MCMCChain, MCMCLikModel, MCMCSampler, MCMCTuner, MCMCRunner

# Abstract types
abstract Model
abstract MCMCModel <: Model

abstract MCMCSampler

abstract MCMCTuner

abstract MCMCRunner

###########  Models  ##########
include("modellers/mcmcmodels.jl")
include("modellers/likmodel.jl")

###########  Autodiff  ##########
include("autodiff/parsing.jl")     #  include model expression parsing function
include("autodiff/diff.jl")        #  include derivatives definitions
include("autodiff/distributions.jl")    #  include distributions definitions

### MCMCTask type, generated by combining a MCMCModel, a MCMCSampler and a MCMCRunner
type MCMCTask
  task::Task
  model::MCMCModel
  sampler::MCMCSampler
  runner::MCMCRunner
end
reset(t::MCMCTask, x) = t.task.storage[:reset](x)

#############  samplers  ########################
include("samplers/samplers.jl")  # Common definitions for samplers
include("samplers/IMH.jl")    # Independent Metropolis-Hastings sampler
include("samplers/RWM.jl")    # Random-walk Metropolis sampler
include("samplers/RAM.jl")    # Robust adaptive Metropolis sampler
include("samplers/MALA.jl")   # Metropolis adjusted Langevin algorithm sampler
include("samplers/HMC.jl")    # Hamiltonian Monte-Carlo sampler
include("samplers/HMCDA.jl")    # Adaptive Hamiltonian Monte-Carlo sampler with dual averaging
include("samplers/NUTS.jl")   # No U-Turn Hamiltonian Monte-Carlo sampler
include("samplers/SMMALA.jl") # Simplified manifold Metropolis adjusted Langevin algorithm sampler
# include("samplers/MMALA.jl")  # Manifold Metropolis adjusted Langevin algorithm sampler (deprecated)
include("samplers/PMALA.jl")  # Position-dependent Metropolis adjusted Langevin algorithm sampler
include("samplers/RMHMC.jl")  # Riemannian manifold Hamiltonian Monte Carlo sampler
include("samplers/ERMLMC.jl") # Explicit Riemannian manifold Lagrangian Monte Carlo (ERMLMC)
include("samplers/RMLMC.jl") # Semi-explicit Riemannian manifold Lagrangian Monte Carlo (RMLMC)

### MCMCChain, the result of running a MCMCTask
type MCMCChain
  range::Range{Int}
  samples::DataFrame
  gradients::DataFrame
  diagnostics::Dict
  task::Union(MCMCTask, Array{MCMCTask})
  runTime::Float64

  function MCMCChain(r::Range{Int}, s::DataFrame, g::DataFrame, d::Dict, t::Union(MCMCTask, Array{MCMCTask}),
    rt::Float64)
    if !isempty(g); @assert size(s) == size(g) "samples and gradients must have the same number of rows and columns"; end
    new(r, s, g, d, t, rt)
  end
end

MCMCChain(r::Range{Int}, s::DataFrame, d::Dict, t::Union(MCMCTask, Array{MCMCTask}), rt::Float64) =
	MCMCChain(r, s, DataFrame(), d, t, rt)
MCMCChain(r::Range{Int}, s::DataFrame, t::Union(MCMCTask, Array{MCMCTask}), rt::Float64) =
	MCMCChain(r, s, DataFrame(), Dict(), t, rt)
MCMCChain(r::Range{Int}, s::DataFrame, d::Dict, t::Union(MCMCTask, Array{MCMCTask})) =
	MCMCChain(r, s, DataFrame(), d, t, NaN)
MCMCChain(r::Range{Int}, s::DataFrame, t::Union(MCMCTask, Array{MCMCTask})) =
	MCMCChain(r, s, DataFrame(), Dict(), t, NaN)

function show(io::IO, res::MCMCChain)
  println("$(ncol(res.samples)) parameters, $(nrow(res.samples)) samples (per parameter), $(round(res.runTime, 1)) sec.")
end

#  Definition of * as a shortcut operator for (model, sampler, runner) combination
*{M<:MCMCModel, S<:MCMCSampler, R<:MCMCRunner}(m::M, s::S, r::R) = spinTask(m, s, r)
*{M<:MCMCModel, S<:MCMCSampler, R<:MCMCRunner}(m::Array{M}, s::S, r::R) = map((me) -> spinTask(me, s, r), m)
*{M<:MCMCModel, S<:MCMCSampler, R<:MCMCRunner}(m::M, s::Array{S}, r::R) = map((se) -> spinTask(m, se, r), s)
*{M<:MCMCModel, S<:MCMCSampler, R<:MCMCRunner}(m::M, s::S, r::Array{R}) = map((re) -> spinTask(m, s, re), r)
*{M<:MCMCModel, S<:MCMCSampler, R<:MCMCRunner}(m::Array{M}, s::Array{S}, r::R) =
  map((me, se) -> spinTask(me, se, r), m, s)
*{M<:MCMCModel, S<:MCMCSampler, R<:MCMCRunner}(m::Array{M}, s::S, r::Array{R}) =
  map((me, re) -> spinTask(me, s, re), m, r)
*{M<:MCMCModel, S<:MCMCSampler, R<:MCMCRunner}(m::M, s::Array{S}, r::Array{R}) =
  map((se, re) -> spinTask(m, se, re), s, r)
*{M<:MCMCModel, S<:MCMCSampler, R<:MCMCRunner}(m::Array{M}, s::Array{S}, r::Array{R}) =
  map((me, se, re) -> spinTask(me, se, re), m, s, r)

#############  runners    ########################
include("runners/runners.jl")
include("runners/SerialMC.jl") # Ordinary serial MCMC runner
include("runners/SerialTempMC.jl") # Serial Tempering Monte-Carlo runner
include("runners/SeqMC.jl") # Sequential Monte-Carlo runner

#############  MCMC output analysis and diagnostics    ########################
include("stats/mean.jl") # MCMC mean estimators
include("stats/var.jl") # MCMC variance estimators
include("stats/ess.jl") # Effective sample size and integrated autocorrelation time functions
include("stats/summary.jl") # Summary statistics for MCMCChain
include("stats/zv.jl")  # ZV-MCMC estimators
end


text:
const BARS = collect("▁▂▃▄▅▆▇█")

function unicode_histogram(data, nbins = 12)
    # @show data
    f = fit(Histogram, data, nbins = nbins)  # nbins: more like a guideline than a rule, really
    # scale weights between 1 and 8 (length(BARS)) to fit the indices in BARS
    # eps is needed so indices are in the interval [0, 8) instead of [0, 8] which could
    # result in indices 0:8 which breaks things
    scaled = f.weights .* (length(BARS) / maximum(f.weights) - eps())
    indices = floor.(Int, scaled) .+ 1
    return join((BARS[i] for i in indices))
end

function precise(df::DataFrame; digits = 2, depth = Inf, alpha = 0.11)
    d = DataFrame()
    d.param = names(df)
    d.mean = mean.(eachcol(df))
    d.std = std.(eachcol(df))
    d[:, "5.5%"] = quantile.(eachcol(df), alpha/2)
    d[:, "50%"] = quantile.(eachcol(df), 0.5)
    d[:, "94.5%"] = quantile.(eachcol(df), 1 - alpha/2)
    u = unicode_histogram.(eachcol(df))
    d.histogram = unicode_histogram.(eachcol(df))

    for col in ["mean", "std", "5.5%", "50%", "94.5%"]
        d[:, col] .= round.(d[:, col], digits = digits)
    end

    pretty_table(d, nosubheader = true, vlines = [0, 1, 7])
end


text:
module ModuleSyntheticData
using MultipleViewGeometry.ModuleOperators, MultipleViewGeometry.ModuleMathAliases
using LinearAlgebra, StaticArrays
export generate_planar_points
include("synthetic_data.jl")
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct PublicIPAddressesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a static or dynamic public IP address.
Param: resourceGroupName::String (required)
Param: publicIpAddressName::String (required)
Param: parameters::PublicIPAddress (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PublicIPAddress
"""
function publicIPAddressesCreateOrUpdate(_api::PublicIPAddressesApi, resourceGroupName::String, publicIpAddressName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", PublicIPAddress, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/publicIPAddresses/{publicIpAddressName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "publicIpAddressName", publicIpAddressName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes the specified public IP address.
Param: resourceGroupName::String (required)
Param: publicIpAddressName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function publicIPAddressesDelete(_api::PublicIPAddressesApi, resourceGroupName::String, publicIpAddressName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/publicIPAddresses/{publicIpAddressName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "publicIpAddressName", publicIpAddressName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the specified public IP address in a specified resource group.
Param: resourceGroupName::String (required)
Param: publicIpAddressName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: PublicIPAddress
"""
function publicIPAddressesGet(_api::PublicIPAddressesApi, resourceGroupName::String, publicIpAddressName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PublicIPAddress, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/publicIPAddresses/{publicIpAddressName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "publicIpAddressName", publicIpAddressName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets all public IP addresses in a resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PublicIPAddressListResult
"""
function publicIPAddressesList(_api::PublicIPAddressesApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PublicIPAddressListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/publicIPAddresses", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets all the public IP addresses in a subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PublicIPAddressListResult
"""
function publicIPAddressesListAll(_api::PublicIPAddressesApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PublicIPAddressListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Network/publicIPAddresses", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export publicIPAddressesCreateOrUpdate, publicIPAddressesDelete, publicIPAddressesGet, publicIPAddressesList, publicIPAddressesListAll


text:
@testset "4.1.13 (d+e x)^m sin(a+b x+c x^2)^n" begin
    (a, b, c, d, e, x, ) = @variables a b c d e x

    #= ::Package:: =#

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

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

    @test_int [x^2*sin(a + b*x + c*x^2), x, 8, (b*cos(a + b*x + c*x^2))/(4*c^2) - (x*cos(a + b*x + c*x^2))/(2*c) + (sqrt(pi/2)*cos(a - b^2/(4*c))*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(2*pi))))/(2*c^(3/2)) + (b^2*sqrt(pi/2)*cos(a - b^2/(4*c))*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(2*pi))))/(4*c^(5/2)) + (b^2*sqrt(pi/2)*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a - b^2/(4*c)))/(4*c^(5/2)) - (sqrt(pi/2)*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a - b^2/(4*c)))/(2*c^(3/2))]
    @test_int [x*sin(a + b*x + c*x^2), x, 4, -(cos(a + b*x + c*x^2)/(2*c)) - (b*sqrt(pi/2)*cos(a - b^2/(4*c))*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(2*pi))))/(2*c^(3/2)) - (b*sqrt(pi/2)*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a - b^2/(4*c)))/(2*c^(3/2))]
    @test_int [sin(a + b*x + c*x^2), x, 3, (sqrt(pi/2)*cos(a - b^2/(4*c))*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(2*pi))))/sqrt(c) + (sqrt(pi/2)*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a - b^2/(4*c)))/sqrt(c)]
    @test_int [sin(a + b*x + c*x^2)/x, x, 0, Unintegrable(sin(a + b*x + c*x^2)/x, x)]
    @test_int [sin(a + b*x + c*x^2)/x^2 - b*cos(a + b*x + c*x^2)/x, x, 5, sqrt(c)*sqrt(2*pi)*cos(a - b^2/(4*c))*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(2*pi))) - sqrt(c)*sqrt(2*pi)*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a - b^2/(4*c)) - sin(a + b*x + c*x^2)/x]

    @test_int [x^2*sin(a + b*x - c*x^2), x, 8, (b*cos(a + b*x - c*x^2))/(4*c^2) + (x*cos(a + b*x - c*x^2))/(2*c) + (sqrt(pi/2)*cos(a + b^2/(4*c))*FresnelC((b - 2*c*x)/(sqrt(c)*sqrt(2*pi))))/(2*c^(3/2)) + (b^2*sqrt(pi/2)*cos(a + b^2/(4*c))*FresnelS((b - 2*c*x)/(sqrt(c)*sqrt(2*pi))))/(4*c^(5/2)) - (b^2*sqrt(pi/2)*FresnelC((b - 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a + b^2/(4*c)))/(4*c^(5/2)) + (sqrt(pi/2)*FresnelS((b - 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a + b^2/(4*c)))/(2*c^(3/2))]
    @test_int [x*sin(a + b*x - c*x^2), x, 4, cos(a + b*x - c*x^2)/(2*c) + (b*sqrt(pi/2)*cos(a + b^2/(4*c))*FresnelS((b - 2*c*x)/(sqrt(c)*sqrt(2*pi))))/(2*c^(3/2)) - (b*sqrt(pi/2)*FresnelC((b - 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a + b^2/(4*c)))/(2*c^(3/2))]
    @test_int [sin(a + b*x - c*x^2), x, 3, (sqrt(pi/2)*cos(a + b^2/(4*c))*FresnelS((b - 2*c*x)/(sqrt(c)*sqrt(2*pi))))/sqrt(c) - (sqrt(pi/2)*FresnelC((b - 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a + b^2/(4*c)))/sqrt(c)]
    @test_int [sin(a + b*x - c*x^2)/x, x, 0, Unintegrable(sin(a + b*x - c*x^2)/x, x)]
    @test_int [sin(a + b*x - c*x^2)/x^2 - b*cos(a + b*x - c*x^2)/x, x, 5, sqrt(c)*sqrt(2*pi)*cos(a + b^2/(4*c))*FresnelC((b - 2*c*x)/(sqrt(c)*sqrt(2*pi))) + sqrt(c)*sqrt(2*pi)*FresnelS((b - 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a + b^2/(4*c)) - sin(a + b*x - c*x^2)/x]

    @test_int [x^2*sin(1/4 + x + x^2), x, 6, (1/4)*cos(1/4 + x + x^2) - (1/2)*x*cos(1/4 + x + x^2) + (1/2)*sqrt(pi/2)*FresnelC((1 + 2*x)/sqrt(2*pi)) + (1/4)*sqrt(pi/2)*FresnelS((1 + 2*x)/sqrt(2*pi))]
    @test_int [x*sin(1/4 + x + x^2), x, 3, (-(1/2))*cos(1/4 + x + x^2) - (1/2)*sqrt(pi/2)*FresnelS((1 + 2*x)/sqrt(2*pi))]
    @test_int [sin(1/4 + x + x^2), x, 2, sqrt(pi/2)*FresnelS((1 + 2*x)/sqrt(2*pi))]
    @test_int [sin(1/4 + x + x^2)/x, x, 0, Unintegrable(sin(1/4 + x + x^2)/x, x)]
    @test_int [sin(1/4 + x + x^2)/x^2, x, 3, sqrt(2*pi)*FresnelC((1 + 2*x)/sqrt(2*pi)) + Unintegrable(cos(1/4 + x + x^2)/x, x) - sin(1/4 + x + x^2)/x]

    @test_int [x^2*sin(a + b*x + c*x^2)^2, x, 10, x^3/6 - (b^2*sqrt(pi)*cos(2*a - b^2/(2*c))*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(pi))))/(16*c^(5/2)) + (sqrt(pi)*cos(2*a - b^2/(2*c))*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(pi))))/(16*c^(3/2)) + (sqrt(pi)*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a - b^2/(2*c)))/(16*c^(3/2)) + (b^2*sqrt(pi)*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a - b^2/(2*c)))/(16*c^(5/2)) + (b*sin(2*a + 2*b*x + 2*c*x^2))/(16*c^2) - (x*sin(2*a + 2*b*x + 2*c*x^2))/(8*c)]
    @test_int [x*sin(a + b*x + c*x^2)^2, x, 6, x^2/4 + (b*sqrt(pi)*cos(2*a - b^2/(2*c))*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(pi))))/(8*c^(3/2)) - (b*sqrt(pi)*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a - b^2/(2*c)))/(8*c^(3/2)) - sin(2*a + 2*b*x + 2*c*x^2)/(8*c)]
    @test_int [sin(a + b*x + c*x^2)^2, x, 5, x/2 - (sqrt(pi)*cos(2*a - b^2/(2*c))*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(pi))))/(4*sqrt(c)) + (sqrt(pi)*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a - b^2/(2*c)))/(4*sqrt(c))]
    @test_int [sin(a + b*x + c*x^2)^2/x, x, 2, (-(1/2))*Unintegrable(cos(2*a + 2*b*x + 2*c*x^2)/x, x) + log(x)/2]

    @test_int [x^2*sin(a + b*x - c*x^2)^2, x, 10, x^3/6 + (b^2*sqrt(pi)*cos(2*a + b^2/(2*c))*FresnelC((b - 2*c*x)/(sqrt(c)*sqrt(pi))))/(16*c^(5/2)) - (sqrt(pi)*cos(2*a + b^2/(2*c))*FresnelS((b - 2*c*x)/(sqrt(c)*sqrt(pi))))/(16*c^(3/2)) + (sqrt(pi)*FresnelC((b - 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a + b^2/(2*c)))/(16*c^(3/2)) + (b^2*sqrt(pi)*FresnelS((b - 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a + b^2/(2*c)))/(16*c^(5/2)) + (b*sin(2*a + 2*b*x - 2*c*x^2))/(16*c^2) + (x*sin(2*a + 2*b*x - 2*c*x^2))/(8*c)]
    @test_int [x*sin(a + b*x - c*x^2)^2, x, 6, x^2/4 + (b*sqrt(pi)*cos(2*a + b^2/(2*c))*FresnelC((b - 2*c*x)/(sqrt(c)*sqrt(pi))))/(8*c^(3/2)) + (b*sqrt(pi)*FresnelS((b - 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a + b^2/(2*c)))/(8*c^(3/2)) + sin(2*a + 2*b*x - 2*c*x^2)/(8*c)]
    @test_int [sin(a + b*x - c*x^2)^2, x, 5, x/2 + (sqrt(pi)*cos(2*a + b^2/(2*c))*FresnelC((b - 2*c*x)/(sqrt(c)*sqrt(pi))))/(4*sqrt(c)) + (sqrt(pi)*FresnelS((b - 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a + b^2/(2*c)))/(4*sqrt(c))]
    @test_int [sin(a + b*x - c*x^2)^2/x, x, 2, (-(1/2))*Unintegrable(cos(2*a + 2*b*x - 2*c*x^2)/x, x) + log(x)/2]

    @test_int [x^2*sin(1/4 + x + x^2)^2, x, 8, x^3/6 - (1/16)*sqrt(pi)*FresnelC((1 + 2*x)/sqrt(pi)) + (1/16)*sqrt(pi)*FresnelS((1 + 2*x)/sqrt(pi)) + (1/16)*sin(1/2 + 2*x + 2*x^2) - (1/8)*x*sin(1/2 + 2*x + 2*x^2)]
    @test_int [x*sin(1/4 + x + x^2)^2, x, 5, x^2/4 + (1/8)*sqrt(pi)*FresnelC((1 + 2*x)/sqrt(pi)) - (1/8)*sin(1/2 + 2*x + 2*x^2)]
    @test_int [sin(1/4 + x + x^2)^2, x, 4, x/2 - (1/4)*sqrt(pi)*FresnelC((1 + 2*x)/sqrt(pi))]
    @test_int [sin(1/4 + x + x^2)^2/x, x, 2, (-(1/2))*Unintegrable(cos(1/2 + 2*x + 2*x^2)/x, x) + log(x)/2]
    @test_int [sin(1/4 + x + x^2)^2/x^2, x, 5, -(1/(2*x)) + cos(1/2 + 2*x + 2*x^2)/(2*x) + sqrt(pi)*FresnelS((1 + 2*x)/sqrt(pi)) + Unintegrable(sin(1/2 + 2*x + 2*x^2)/x, x)]

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

    @test_int [(d + e*x)^2*sin(a + b*x + c*x^2), x, 8, -((e*(2*c*d - b*e)*cos(a + b*x + c*x^2))/(4*c^2)) - (e*(d + e*x)*cos(a + b*x + c*x^2))/(2*c) + (e^2*sqrt(pi/2)*cos(a - b^2/(4*c))*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(2*pi))))/(2*c^(3/2)) + ((2*c*d - b*e)^2*sqrt(pi/2)*cos(a - b^2/(4*c))*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(2*pi))))/(4*c^(5/2)) + ((2*c*d - b*e)^2*sqrt(pi/2)*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a - b^2/(4*c)))/(4*c^(5/2)) - (e^2*sqrt(pi/2)*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a - b^2/(4*c)))/(2*c^(3/2))]
    @test_int [(d + e*x)^1*sin(a + b*x + c*x^2), x, 4, -((e*cos(a + b*x + c*x^2))/(2*c)) + ((2*c*d - b*e)*sqrt(pi/2)*cos(a - b^2/(4*c))*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(2*pi))))/(2*c^(3/2)) + ((2*c*d - b*e)*sqrt(pi/2)*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a - b^2/(4*c)))/(2*c^(3/2))]
    @test_int [(d + e*x)^0*sin(a + b*x + c*x^2), x, 3, (sqrt(pi/2)*cos(a - b^2/(4*c))*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(2*pi))))/sqrt(c) + (sqrt(pi/2)*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(2*pi)))*sin(a - b^2/(4*c)))/sqrt(c)]
    @test_int [sin(a + b*x + c*x^2)/(d + e*x)^1, x, 0, Unintegrable(sin(a + b*x + c*x^2)/(d + e*x), x)]

    @test_int [(d + e*x)^2*sin(a + b*x + c*x^2)^2, x, 10, (d + e*x)^3/(6*e) - ((2*c*d - b*e)^2*sqrt(pi)*cos(2*a - b^2/(2*c))*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(pi))))/(16*c^(5/2)) + (e^2*sqrt(pi)*cos(2*a - b^2/(2*c))*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(pi))))/(16*c^(3/2)) + (e^2*sqrt(pi)*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a - b^2/(2*c)))/(16*c^(3/2)) + ((2*c*d - b*e)^2*sqrt(pi)*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a - b^2/(2*c)))/(16*c^(5/2)) - (e*(2*c*d - b*e)*sin(2*a + 2*b*x + 2*c*x^2))/(16*c^2) - (e*(d + e*x)*sin(2*a + 2*b*x + 2*c*x^2))/(8*c)]
    @test_int [(d + e*x)^1*sin(a + b*x + c*x^2)^2, x, 6, (d + e*x)^2/(4*e) - ((2*c*d - b*e)*sqrt(pi)*cos(2*a - b^2/(2*c))*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(pi))))/(8*c^(3/2)) + ((2*c*d - b*e)*sqrt(pi)*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a - b^2/(2*c)))/(8*c^(3/2)) - (e*sin(2*a + 2*b*x + 2*c*x^2))/(8*c)]
    @test_int [(d + e*x)^0*sin(a + b*x + c*x^2)^2, x, 5, x/2 - (sqrt(pi)*cos(2*a - b^2/(2*c))*FresnelC((b + 2*c*x)/(sqrt(c)*sqrt(pi))))/(4*sqrt(c)) + (sqrt(pi)*FresnelS((b + 2*c*x)/(sqrt(c)*sqrt(pi)))*sin(2*a - b^2/(2*c)))/(4*sqrt(c))]
    @test_int [sin(a + b*x + c*x^2)^2/(d + e*x)^1, x, 2, (-(1/2))*Unintegrable(cos(2*a + 2*b*x + 2*c*x^2)/(d + e*x), x) + log(d + e*x)/(2*e)]
    end

max_stars_repo_path: test/4 Trig functions/4.1 Sine/4.1.13 (d+e x)^m sin(a+b x+c x^2)^n.jl


text:
"""

# cmdline

Recursively parse the model to construct command line.

### Method
```julia
cmdline(m)
```

### Required arguments
```julia
* `m::CmdStanSampleModel`                : CmdStanSampleModel
```

### Related help
```julia
?CmdStanSampleModel                      : Create a CmdStanSampleModel
```
"""
function cmdline(m::HelpModel, id)

  #=
  `.../help help`
  =#

  cmd = ``
  if isa(m, HelpModel)
    # Inserts the executable for unix and windows
    cmd = `$(m.exec_path)`

    # `help` specific portion of the model
    #cmd = `$cmd $(split(lowercase(string(typeof(m))), '.')[end])`
    cmd = `$cmd $(getfield(m.method, :help)) help`
  end

  cmd

end


text:
__precompile__()

"""
A Julia package for accessing the Bitly API. See `BitlyToken`.
"""
module Bitly

#imports
using HTTP
using JSON
using Base64: base64encode
using Printf

#exports
export BitlyToken, requestBitlyToken, shorten, expand, clicks

const BITLY_API_BASE = "https://api-ssl.bitly.com"
const BITLY_SHORTEN = joinpath(BITLY_API_BASE,"v4/shorten")
const BITLY_EXPAND = joinpath(BITLY_API_BASE,"v4/expand")
const BITLY_BITLINKS = joinpath(BITLY_API_BASE,"v4/bitlinks")

const TOKEN_LENGTH = 40
const TOKEN_ENV_NAME = "BITLY_ACCESS_TOKEN"
const TOKEN_FILE_NAME = ".bitlyrc"

# Bitly connection type
"""
A connection to the Bitly API.

Constructors
------------
- `BitlyToken()`: Token detected automatically. First, looks for the environment variable
`BITLY_ACCESS_TOKEN`, then looks for the file `~/.bitlyrc`.
- `BitlyToken(token::AbstractString)`: User specifies token directly

See [`requestBitlyToken`](@ref).
"""
mutable struct BitlyToken
    token::AbstractString
    url::AbstractString
    function BitlyToken(token, url)
        # token validation
        if length(token) > TOKEN_LENGTH
            token = token[1:TOKEN_LENGTH]
            @warn("Bitly API token too long. First $(TOKEN_LENGTH) chars used.")
        elseif length(token) < TOKEN_LENGTH
            error("Invalid Bitly API token -- token too short: $(token)")
        end
        if !all(isxdigit, token)
            error("Invalid Bitly API token -- invalid characters: $(token)")
        end
        return new(token, url)
    end
end

BitlyToken(token::AbstractString) = BitlyToken(token, BITLY_API_BASE)

token_file() = joinpath(homedir(), TOKEN_FILE_NAME)

function load_bitly_token()
    if TOKEN_ENV_NAME in keys(ENV)
        ENV[TOKEN_ENV_NAME]
    elseif isfile(token_file())
        open(token_file(), "r") do file
            rstrip(read(file, String))
        end
    else
        error("Bitly API Key not detected.")
    end
end

has_bitly_token() = TOKEN_ENV_NAME in keys(ENV) || isfile(token_file())

function BitlyToken(;verbose=false)
    token = load_bitly_token()
    verbose && println("API token loaded.")
    return BitlyToken(token)
end

function Base.show(io::IO, b::BitlyToken)
    @printf io "Bitly API Connection\n"
    @printf io "\turl: %s\n" b.url
    @printf io "\ttoken: %s\n" b.token
end

"""
`requestBitlyToken(uname,pw)`

Exchanges your username and password for an access token. Your username and password are *not stored*!

Returns an initialized [`BitlyToken`](@ref) object.

Note: you can also log on to your account on Bitly and request an access token directly.
"""
function requestBitlyToken(uname::AbstractString,pw::AbstractString)
    cred = base64encode(string(uname,":",pw))
    header = ("Content-Type" => "application/x-www-form-urlencoded","Authorization" => "Basic $cred")
    r = HTTP.post(joinpath(BITLY_API_BASE,"oauth/access_token"),header)
    tok = String(r.body)
    return BitlyToken(tok)
end

"""
`shorten(b::Bitly,url::AbstractString)`

Shorten `url` to a `bit.ly/xXxXx` shortened link.

Returns a `NamedTuple` with `link` and `response` items. `link` is the `bit.ly` link,
while `response` is the complete `JSON` response from the Bitly API.

```julia
b = BitlyToken()
S = shorten(b,"https://docs.julialang.org/en/v1/")
S.link
S.response
```
"""
function shorten(b::BitlyToken,url::AbstractString)
    header = ("Content-Type" => "application/json","Authorization" => "Bearer $(b.token)")
    payload = Dict("long_url" => url)
    r = HTTP.post(BITLY_SHORTEN,header,JSON.json(payload))

    js = JSON.parse(String(r.body))
    link = js["link"]

    return (link=link,response=js)
end

"""
`expand(b::BitlyToken,link::AbstractString)`

Recovers the original `url` from a `bit.ly/xXxXx` shortened link.

Returns a `NamedTuple` with `long_url` and `response` items. `long_url` is the original link,
while `response` is the complete `JSON` response from the Bitly API.

```julia
b = BitlyToken()
S = expand(b,"http://bit.ly/2Us5Vl7")
S.long_url
S.response
```
"""
function expand(b::BitlyToken,link::AbstractString)
    s = HTTP.splitpath(link)
    id = joinpath(s[end-1],s[end])

    header = ("Content-Type" => "application/json","Authorization" => "Bearer $(b.token)")
    payload = Dict("bitlink_id" => id)
    r = HTTP.post(BITLY_EXPAND,header,JSON.json(payload))

    js = JSON.parse(String(r.body))
    long_url = js["long_url"]

    return (long_url=long_url,response=js)
end

"""
`clicks(b::BitlyToken,link; summary=false,unit="day",units=-1,size=50)`

Get click information for shortened link `link`.

`summary=true` provides a single count, otherwise returns a time series.

`unit` can be `"minute"`,`"hour"`,`"day"`,`"week"`,`"month"`.

`units` is an integer representing the number of time units to query data for. Pass `-1` to return all units available.

`size` is the quantity of items to be returned.

Returns a `Dict`.
"""
function clicks(b::BitlyToken,link::AbstractString; summary::Bool=false,unit::String="day",units::Int=-1,size::Int=50)
    s = HTTP.splitpath(link)
    id = joinpath(s[end-1],s[end])

    header = ("Content-Type" => "application/json","Authorization" => "Bearer $(b.token)")
    p = summary ? joinpath(BITLY_BITLINKS,id,"clicks","summary") : joinpath(BITLY_BITLINKS,id,"clicks")
    r = HTTP.get(p,header)

    js = JSON.parse(String(r.body))

    return js
end

end # module