text:
# Profiles transit_poly.jl
#using ForwardDiff
#using DiffResults

include("../../../src/transit_poly_struct.jl")
#include("../../../src/transit_structure.jl")

# Define function for autodiff:

#function tp_auto(r::T,b::T,u_n::Array{T,1}) where {T <: Real}
#
## Set up call:
#  x= zeros(T,2+length(u_n))
#  x[1] = r; x[2] = b
#  for i=1:length(u_n)
#    x[2+i] = u_n[i]
#  end
#
#  # Now, define a wrapper of s2 for use with ForwardDiff:
#  function diff_tp(x::Array{T,1}) where {T <: Real}
#  # x should be a two-element vector with values [r,b]
#  rc = x[1]; bc = x[2]; uc = x[3:length(x)]
#  return transit_poly(rc,bc,uc)
#  end
#
## Now apply autodiff:
#  # Set up a type to store s_n and it's Jacobian with respect to x:
#  out = DiffResults.GradientResult(x)
#  # Compute the Jacobian (and value):
#  out = ForwardDiff.gradient!(out,diff_tp,x)
#  # Place the value in the s_2 vector:
#  tp_flux = DiffResults.value(out)
#  # And, place the Jacobian in an array:
#  tp_grad = DiffResults.gradient(out)
#return tp_flux,tp_grad
#end

#include("../../../src/compute_c_n.jl")

#include("../../../src/dJv_seriesdk.jl")
using PyPlot

# Now, run the same with transit_poly:

function profile_transit_poly(trans)
nb = 1000001; b=zeros(nb)
flux = zeros(nb); fgrad = zeros(4)
flux_grad = zeros(nb,4)
tic()
for i=0:1000000
  b[i+1] = sqrt(((i-500000.0)/500000.0*(1.0+2.0*trans.r))^2)
  trans.b = b[i+1]
  flux[i+1]= transit_poly!(trans)
  flux_grad[i+1,:] = trans.dfdrbu
end
toc()
plot(b,flux)
println("b: ",minimum(b)," ",maximum(b))
println("f: ",minimum(flux)," ",maximum(flux))
return
end

r = 0.1  # radius ratio is 0.1
u = [0.2,0.3]
b=0.0
trans = transit_init(r,b,u,true)

profile_transit_poly(trans)

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
const hpf_all_orders = 1:28
const hpf_flux_hdu = 2
const hpf_var_hdu = 5
const hpf_lambda_hdu = 8

const hpf_manifest_format = [
		    ManifestFormatEntry(:UTTimestamp,String, 4:26),
                    ManifestFormatEntry(:UTDate,String,     30:37),
                    ManifestFormatEntry(:ObsNum,Int,        43:46),
                    ManifestFormatEntry(:Frame,String,      50:89),
                    ManifestFormatEntry(:ObsType,String,    97:99),
                    ManifestFormatEntry(:iTime,Float64,    100:117),
                    ManifestFormatEntry(:QProg,String,     118:145),
                    ManifestFormatEntry(:Object,String,    146:190)
                    ]

ext: jl


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

using Documenter, ArtGallery

makedocs(
    sitename = "ArtGallery.jl",
    authors  = "Zenan Huo",
    format   = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    pages = [
        "Home" => "index.md",
        "API Reference" => "api.md"
    ]
)

deploydocs(
    target = "build",
    repo   = "github.com/ZenanH/ArtGallery.jl.git",
    branch = "gh-pages"
)


text:
function bool3(model::CAGD.Model)
    arranged_model = CAGD.spatial_arrangement(model, outerCell = true)
    innerpoints, intersectedfaces = CAGD.internalpoints(arranged_model)
    # number of input complexes
    nInC = size(model, 3, 1)
    infaces = [CAGD.getModelLoCell(model, 3, i) for i = 1 : nInC]
    # Determine in which initial element was h
    span(h) = [j for j = 1 : nInC if h in infaces[j]]
	containmenttest = CAGD.testinternalpoint(model)
	# currently copCF contains the outercell in first column ...
	# TODO remove first row and column, in case (look at file src/)
	boolmatrix = BitArray(undef, length(innerpoints)+1, nInC+1)
	boolmatrix[1,1] = 1
	for (k, point) in enumerate(innerpoints) # k runs on columns
		cells = containmenttest(point) # contents of columns
		#println(k," ",faces)
		rows = [span(h) for h in cells]
		for l in cat(rows)
			boolmatrix[k+1, l+1] = 1
		end
	end
	return arranged_model, boolmatrix
end

function internalpoints(model::CAGD.Model)
    return Lar.internalpoints(model.G, model.T[1], model.T[2], model.T[3][2:end, :])
end

"""
	testinternalpoint(V::Lar.Points, EV::Lar.Cells, FV::Lar.Cells)

"""
function testinternalpoint(model)
	copEV = abs.(model.T[1])
    copFE = abs.(model.T[2])
    V = model.G
    EV = Lar.cop2lar(model.T[1])
    FV = Lar.cop2lar(map(x -> floor(Int8, x / 2), copFE * copEV))

	I,J,Val = findnz(copFE)
	Val = convert(Array{Int8,1},Val)
	copFE = sparse(I,J,Val)
	function testinternalpoint0(testpoint)
		intersectedfaces = Int64[]
		# spatial index for possible intersections with ray
		faces = Lar.spaceindex(testpoint)((V, FV))
		depot = []
		# face in faces :  indices of faces of possible intersection with ray
		for face in faces
			value = Lar.rayintersection(testpoint)(V, FV, face)
			if typeof(value) == Array{Float64,1} push!(depot, (face,value)) end
		end
		# actual containment test of ray point in faces within depot
		for (face,point3d) in depot
			vs, edges, point2d = Lar.planemap(V, copEV, copFE, face)(point3d)
			classify = Lar.pointInPolygonClassification(vs,edges)
			inOut = classify(point2d)
			if inOut != "p_out"
				push!(intersectedfaces,face)
			end
		end
		return intersectedfaces
	end
	return testinternalpoint0
end


text:
function reserve_slacks(
    container::OptimizationContainer,
    service::T,
) where {T <: PSY.Reserve}
    time_steps = get_time_steps(container)
    variable = add_variable_container!(
        container,
        ReserveRequirementSlack(),
        T,
        PSY.get_name(service),
        time_steps,
    )

    for t in time_steps
        variable[t] = JuMP.@variable(
            container.JuMPmodel,
            base_name = "slack_{$(PSY.get_name(service)), $(t)}",
            lower_bound = 0.0
        )
        add_to_objective_invariant_expression!(container, variable[t] * SERVICES_SLACK_COST)
    end
    return variable
end


text:
# Resource constrained shortest path problem

using SDDiP, JuMP, Ipopt, Base.Test, GLPKMathProgInterface

const JuMPversion = Pkg.installed("JuMP")

struct RCSPData
    A::Array{Int,2}
    costs::Vector{Float64}
    resources::Vector{Float64}
    capacity::Float64
end

function getdata(arcs::Vector{Tuple{Int, Int}}, costs::Vector{Float64}, resources::Vector{Float64}, capacity::Float64)
    narcs = length(arcs)
    @assert narcs == length(costs) == length(resources)
    nnodes = maximum(maximum.(arcs))
    A = zeros(nnodes, narcs)
    for (i, a) in enumerate(arcs)
        A[a[1], i] = 1
        A[a[2], i] = -1
    end
    RCSPData(A, costs, resources, capacity)
end

function getmodel(arcs::Vector{Tuple{Int, Int}}, costs::Vector{Float64}, resources::Vector{Float64}, capacity::Float64)
    data = getdata(arcs, costs, resources, capacity)

    nnodes, narcs = size(data.A)

    m = Model(solver=GLPKSolverMIP())
    # Include arc or don't include arc
    @variable(m, 0 <= x[i=1:narcs] <= 1, Int)
    # Flow in = flow out
    @constraint(m, data.A * x .== 0)
    # Something has to flow
    @constraint(m, x[end] == 1)
    # Resource constraint
    complicating = @constraint(m, dot(data.resources, x) <= data.capacity)
    # Minimise costs
    @objective(m, :Min, dot(data.costs, x))

    return m, complicating
end

arcs = [
(1, 2)
(1, 3)
(2, 5)
(1, 4)
(1, 6)
(2, 3)
(2, 4)
(3, 5)
(3, 6)
(3, 7)
(4, 5)
(4, 7)
(5, 7)
(6, 7)
(7, 1) # last arc is a dummy arc
]

costs = [12; 12; 20; 18; 11; 12; 8; 1; 6; 10; 3; 9; 4; 5; 0]
resources = [5; 18; 3; 12; 24; 3; 4; 6; 1; 10; 7; 16; 14; 16; 0]
capacity = 32.0

# Sove MIP
model, complicating = getmodel(arcs, float(costs), float(resources), capacity)
@assert solve(model, relaxation=false) == :Optimal
println("The MIP objective is: ", getobjectivevalue(model))

# Solve using a linear relaxation
@assert solve(model, relaxation=true) == :Optimal
lpbound = getobjectivevalue(model)
println("The LP relaxation objective is: ", lpbound)

# ========================================================
# Solve using a Lagrangian relaxation

# It doesn't matter whether we relax integrality as well
if JuMPversion < v"0.17"
    for xx in getvariable(model, :x)
        setcategory(xx, :Cont)
    end
else
    for xx in model[:x]
        setcategory(xx, :Cont)
    end
end
setsolver(model, GLPKSolverLP())

println("\n**** Solved using the subgradient method. **** ")

π0 = [0.2]
dualbound = 50.0
relaxed_bound = 100.0;

subgradient = SubgradientMethod(initialbound=dualbound, wait=20)
RCSPdata    = LinearProgramData(model.obj,
                                [complicating],
                                [relaxed_bound],
                                method=subgradient)

status, bestvalue = lagrangiansolve!(RCSPdata, model, π0)
# (Although we could have solved the relaxed problem much more efficiently, it's a shortest path)
println("The optimal value is $bestvalue")
println("and the best multiplier is $π0.")

# In this problem the Lagrangian and LP relaxation bounds are the same
@test isapprox(bestvalue, lpbound, atol=1e-6)
@test isapprox(π0[1], 0.5, atol=1e-4)
@test isapprox(RCSPdata.old_bound[1], 32.0, atol=1e-9)
constr = model.linconstr[complicating.idx]
@test isapprox(getdual(complicating), -0.5, atol=1e-4)
@test isapprox(constr.ub, 32.0, atol=1e-9)
@test constr.lb == -Inf

println("\n**** Solved using the level method. **** ")

levelmethod = LevelMethod(initialbound=dualbound, quadsolver=IpoptSolver(print_level=0))
RCSPdata    = LinearProgramData(model.obj,
                                [complicating],
                                [relaxed_bound],
                                method=levelmethod)

status, bestvalue = lagrangiansolve!(RCSPdata, model, π0)
println("The optimal value is $bestvalue")
println("and the best multiplier is $π0.")

# In this problem the Lagrangian and LP relaxation bounds are the same
@test isapprox(bestvalue, lpbound, atol=1e-9)
@test isapprox(π0[1], 0.5, atol=1e-4)


text:
__precompile__(true)

module PSO
    function make_constraints(::Type{Val{nothing}}, args, kwargs, verbose)
        verbose && println("No constraints given.")
        return x -> [0.0]
    end

    function make_constraints(eqs::Vector, args, kwargs, verbose)
        verbose && println("Converting ieqcons to a single constraint function.")
        return x -> [f(x, args...; kwargs...) for f in eqs]
    end

    function make_constraints(eqs::Function, args, kwargs, verbose)
        verbose && println("Single constraint function given in f_ieqcons.")
        return x -> eqs(x, args...; kwargs...)
    end

    make_constraints(eqs, args, kwargs, verbose) = make_constraints(Val{eqs}, args, kwargs, verbose)

    function update_position!(x, p, fx, fp, fs)
        i_update = (fx .< fp) .& fs
        p[i_update, :] = copy(x[i_update, :])
        fp[i_update, :] = fx[i_update, :]
    end

    function pso(func::Function, lb::Vector, ub::Vector, constraints, args, kwargs,
                 swarmsize, ω, ϕp, ϕg, maxiter, minstep, minfunc, verbose)
        assert(length(ub) == length(lb))
        assert(all(ub .> lb))

        obj = x -> func(x, args...; kwargs...)
        cons = make_constraints(constraints, args, kwargs, verbose)
        is_feasible = x -> all(cons(x) .>= 0)

        # Initialize the particle swarm
        vhigh = abs.(ub .- lb)
        vlow = -vhigh
        S = swarmsize
        D = length(lb)  # the number of dimensions each particle has

        x = lb' .+ rand(S, D) .* (ub .- lb)'  # particle positions
        v = vlow' .+ rand(S, D) .* (vhigh .- vlow)'  # particle velocities
        p = zeros(S, D)  # best particle positions

        fx = [obj(x[i, :]) for i = 1:S]  # current particle function values
        fs = [is_feasible(x[i, :]) for i = 1:S]  # feasibility of each particle
        fp = ones(S) * Inf  # best particle function values

        g = copy(x[1, :])  # best swarm position
        fg = Inf  # best swarm position starting value

        # Store particle's best position (if constraints are satisfied)
        update_position!(x, p, fx, fp, fs)

        # Update swarm's best position
        i_min = indmin(fp)
        if fp[i_min] < fg
            g = copy(p[i_min, :])
            fg = fp[i_min]
        end

        # Iterate until termination criterion met
        it = 1
        while it <= maxiter
            rp = rand(S, D)
            rg = rand(S, D)

            # Update the particles' velocities and positions
            v = ω*v .+ ϕp*rp.*(p .- x) .+ ϕg*rg.*(g' .- x)
            x += v
            # Correct for bound violations
            maskl = x .< lb'
            masku = x .> ub'
            x = x.*(.~(maskl .| masku)) .+ lb'.*maskl .+ ub'.*masku

            # Update objectives and constraints
            for i = 1:S
                fx[i] = obj(x[i, :])
                fs[i] = is_feasible(x[i, :])
            end

            # Store particle's best position (if constraints are satisfied)
            update_position!(x, p, fx, fp, fs)

            # Compare swarm's best position with global best position
            i_min = indmin(fp)
            if fp[i_min] < fg
                verbose && println("New best for swarm at iteration $(it): $(p[i_min, :]) $(fp[i_min])")

                p_min = copy(p[i_min, :])
                stepsize = √(sum((g .- p_min).^2))

                if abs.(fg .- fp[i_min]) <= minfunc
                    verbose && println("Stopping search: Swarm best objective change less than $(minfunc)")
                    return (g, fg, p, fp)
                end
                if stepsize <= minstep
                    verbose && println("Stopping search: Swarm best position change less than $(minstep)")
                    return (g, fg, p, fp)
                end

                g = copy(p_min)
                fg = fp[i_min]
            end

            verbose && println("Best after iteration $(it): $(g) $(fg)")
            it += 1
        end

        println("Stopping search: maximum iterations reached --> $(maxiter)")
        is_feasible(g) || print("However, the optimization couldn't find a feasible design. Sorry")
        return (g, fg, p, fp)
    end

    function pso(func, lb, ub; constraints=nothing, args=(), kwargs=Dict(), swarmsize=100,
                 omega=0.5, phip=0.5, phig=0.5, maxiter=100, minstep=1e-8, minfunc=1e-8,
                 verbose=false, particle_output=false)
        g, fg, p, fp = pso(func, lb, ub, constraints, args, kwargs,
            swarmsize, omega, phip, phig, maxiter, minstep, minfunc, verbose)
        return particle_output? (g, fg, p, fp) : (g, fg)
    end

    export pso;

end # module


text:
# This file was generated, do not modify it. # hide
df.price = df.price ./ 1000;


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

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

# lattice visualizations are drawn with Luxor because it is easier to do 2D drawings with Luxor than with Makie.

function draw(tilebasis::AbstractBasis, tilesize, i, j, color, name)
    vertices = Repeat.tilevertices(tilebasis)
    tile = tilesize * [Luxor.Point(vertices[:,i]...) for i in 1:size(vertices)[2]]
    pt = tilesize * tilebasis[i,j]
    offset = Luxor.Point(pt[1], -pt[2]) # flip y so indices show up correctly

    Luxor.translate(offset)

    Luxor.sethue(color)

    Luxor.poly(tile, :fill, close=true)
    Luxor.sethue("lightgrey")
    Luxor.poly(tile, :stroke, close=true)
    Luxor.sethue("black")

    # scale and offset text so coordinates are readable
    Luxor.fontsize(tilesize / 3)
    Luxor.text("$i, $j", Luxor.Point(-tilesize / 3, tilesize / 2.5))
    if name !== nothing
        Luxor.text(name, Luxor.Point(-tilesize / 4, -tilesize / 4))
    end

    Luxor.translate(-offset)
end

"""Draws a list of hexagonal cells, represented by their lattice coordinates, which are represented as a 2D matrix, with each column being one lattice coordinate."""
function drawcells(tilebasis::AbstractBasis, tilesize, cells::AbstractMatrix; color::Union{AbstractArray,Nothing}=nothing, name::Union{AbstractArray{String},Nothing}=nothing, format=:png, resolution=(500, 500))

    numcells = size(cells)[2]
    Luxor.Drawing(resolution[1], resolution[2], format)
    Luxor.origin()
    Luxor.background(Colors.RGBA(0, 1, 1, 0.0))
    if color === nothing
        color = Colors.distinguishable_colors(numcells, lchoices=30:250) # type unstable but not performance critical code
    end

    for i in 1:numcells
        cell = cells[:,i]
        cellname = name === nothing ? nothing : name[i]
        draw(tilebasis, tilesize, cell[1], cell[2], color[i], cellname)
    end

    Luxor.finish()
    if (format == :svg)
        res = Luxor.svgstring()
        return res
    end
    if (format == :png)
        Luxor.preview()
    end
end

""" draw the LatticeCluster offset to (0,0) """
draw(clstr::LatticeCluster, cluster_coordinate_offset=[0;0;;], scale=50.0) = drawcells(clstr.elementbasis, scale, clustercoordinates(clstr, cluster_coordinate_offset[1,1],cluster_coordinate_offset[2,1]))

""" draw the ClusterWithProperties at coordinates specified by lattice_coordinate_offset """
function draw(clstr::Repeat.ClusterWithProperties, cluster_coordinate_offset=[0;0;;], scale=50.0) # this is how you have to make a 2x1 matrix. One would expect [0;0] to work but it doesn't.
    dims = size(cluster_coordinate_offset)
    @assert dims[1] == 2 "dims[1] must be 2 instead was $(dims[1])"
    clstrsize = clustersize(clstr)
    points = Matrix(undef, dims[1], dims[2] * clstrsize)
    for i in 1:dims[2]
        points[:,(i - 1) * clstrsize + 1:i * clstrsize] = clustercoordinates(clstr, cluster_coordinate_offset[1,i], cluster_coordinate_offset[2,i])
    end
    props = repeat(Repeat.properties(clstr), dims[2])
    drawcells(elementbasis(cluster(clstr)), scale, points, color=props[:,:Color], name=props[:,:Name])
end


text:
function signature(ir::YaoIR)
    defs = Dict(:name => ir.name, :args => ir.args)
    if !isempty(ir.whereparams)
        defs[:whereparams] = ir.whereparams
    end
    return defs
end

function build_codeinfo(m::Module, defs::Dict, ir::IR)
    defs[:body] = :(return)
    @timeit_debug to "make empty CI" ci = Meta.lower(m, combinedef(defs))
    mt = ci.args[].code[end-1]
    mt_ci = mt.args[end]

    # update method CodeInfo
    @timeit_debug to "copy IR" ir = copy(ir)
    @timeit_debug to "insert self" Inner.argument!(ir, at = 1)
    @timeit_debug to "update CI" Inner.update!(mt_ci, ir)
    @timeit_debug to "validate CI" Core.Compiler.validate_code(mt_ci)
    return ci
end

function build_codeinfo(ir::YaoIR)
    build_codeinfo(ir.mod, signature(ir), ir.body)
end

function arguements(ir::YaoIR)
    args = map(rm_annotations, ir.args)
    if (ir.name isa Expr) && (ir.name.head === :(::))
        insert!(args, 1, ir.name.args[1])
    end
    return args
end

"""
    rm_annotations(x)

Remove type annotation of given expression.
"""
function rm_annotations(x)
    x isa Expr || return x
    if x.head == :(::)
        return x.args[1]
    elseif x.head in [:(=), :kw] # default values
        return rm_annotations(x.args[1])
    else
        return x
    end
end

function annotations(x)
    x isa Expr || return x
    if x.head == :(::)
        return x.args[2]
    elseif x.head in [:(=), :kw]
        return annotations(x.args[1])
    else
        return x
    end
end

function splatting_variables(variables, free)
    Expr(:(=), Expr(:tuple, variables...), free)
end

# function argtypes(def::Dict)
#     ex = Expr(:curly, :Tuple)
#     if haskey(def, :args)
#         for each in def[:args]
#             if each isa Symbol
#                 push!(ex.args, :Any)
#             elseif (each isa Expr) && (each.head === :(::))
#                 push!(ex.args, each.args[end])
#             end
#         end
#     end

#     return ex
# end

generic_circuit(name::Symbol) = Expr(:curly, GlobalRef(YaoLang, :GenericCircuit), QuoteNode(name))
circuit(name::Symbol) = Expr(:curly, GlobalRef(YaoLang, :Circuit), QuoteNode(name))
# custom struct
generic_circuit(name) = annotations(name)
circuit(name) =
    Expr(:curly, GlobalRef(YaoLang, :Circuit), QuoteNode(gensym(Symbol(annotations(name)))))

to_locations(x) = :(Locations($x))
to_locations(x::Int) = Locations(x)


text:
module ShootoutBenchmarks

# This module contains the Julia implementations of the benchmarks used by the
# Computer Language Benchmarks Game (http://benchmarksgame.alioth.debian.org/),
# previously known as the Great Computer Language Shootout.
#
# See https://github.com/JuliaLang/julia/issues/660 for details.

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools
using Compat

if v"0.7.0-DEV.2915" <= VERSION < v"0.7.0-DEV.3393"
    using Unicode
end
if VERSION >= v"0.7.0-DEV.3052"
    using Printf
end

const SUITE = BenchmarkGroup(["example", "regex", "pi", "π", "tree"])
const SHOOTOUT_DATA_PATH = joinpath(dirname(@__FILE__), "data")

##################################################
# Allocate and deallocate many many binary trees #
##################################################

include("binary_trees.jl")
SUITE["binary_trees"] = @benchmarkable perf_binary_trees(10)

###########################################
# Indexed-access to tiny integer-sequence #
###########################################

include("fannkuch.jl")
SUITE["fannkuch"] = @benchmarkable perf_fannkuch(7)

###########################################
# Generate and write random DNA sequences #
###########################################

include("fasta.jl")
SUITE["fasta"] = @benchmarkable perf_fasta(100)

#############################################
# Hashtable update and k-nucleotide strings #
#############################################

include("k_nucleotide.jl")
SUITE["k_nucleotide"] = @benchmarkable perf_k_nucleotide()

################################################
# Generate Mandelbrot set portable bitmap file #
################################################

include("mandelbrot.jl")
SUITE["mandelbrot"] = @benchmarkable perf_mandelbrot(200)

################################################
# Search for solutions to shape packing puzzle #
################################################

include("meteor_contest.jl")
SUITE["meteor_contest"] = @benchmarkable perf_meteor_contest()

######################################
# Double-precision N-body simulation #
######################################

include("nbody.jl")
include("nbody_vec.jl")
SUITE["nbody"] = @benchmarkable NBody.perf_nbody()
SUITE["nbody_vec"] = @benchmarkable NBodyVec.perf_nbody_vec()

############################################
# Streaming arbitrary-precision arithmetic #
############################################

include("pidigits.jl")
SUITE["pidigits"] = @benchmarkable perf_pidigits(1000)

#############################################################
# Match DNA 8-mers and substitute nucleotides for IUB codes #
#############################################################

include("regex_dna.jl")
SUITE["regex_dna"] = @benchmarkable perf_regex_dna()

#######################################################
# Read DNA sequences - write their reverse-complement #
#######################################################

include("revcomp.jl")
SUITE["revcomp"] = @benchmarkable perf_revcomp() time_tolerance=0.25

#####################################
# Eigenvalue using the power method #
#####################################

include("spectralnorm.jl")
SUITE["spectralnorm"] = @benchmarkable perf_spectralnorm()

end # module


text:
wavelengths(spectrum::Spectrum_29) = 1 ./ spectrum.cyclefrequency

function amplitude_interpolated(spectrum::Spectrum_29, wavelength::Length)
    λs = wavelengths(spectrum)
    as = spectrum.spectralamplitude
    iabove = findlast(x -> isless(wavelength, x), λs)
    ibelow = min(iabove + 1, length(λs))
    Δx = λs[iabove] - λs[ibelow]
    Δy = as[iabove] - as[ibelow]
    frac = (wavelength - λs[ibelow] ) / Δx
    as[ibelow] + frac * Δy
end


text:
module GMM

#main module

end

max_stars_repo_name: Michael-Abbott-W/GMM

max_issues_repo_name: Michael-Abbott-W/GMM

max_forks_repo_name: Michael-Abbott-W/GMM


text:
function nonstandard1(f)
  symbol_f = Symbol("@" * string(f) )
  macro_f = Expr(:quote, Expr(:macrocall, symbol_f))
  string_f = "See documentation of [`$f`](@ref)"
  quote
      macro $f(args...)
          esc($f(args...) )
      end
      @doc $string_f $macro_f
  end
end

export nonstandard
"""
    @nonstandard(fs...)

Will create a nonstandard evaluation macro for each of the `fs` functions.

Each function should be a function that takes and returns expressions. The
nonstandard macro will have the same name but will take in code, not
expressions, and will evaluate the result locally when the macro is called. Will
write a docstring for the nonstandard version pointing to the documentation of
the standard version.

# Examples
```julia
binaryfun(a, b, c) = Expr(:call, b, a, c)
chainback(a, b, c) = Expr(:call, c, b, a)

@nonstandard binaryfun chainback

@test vcat(1, 2) == @binaryfun 1 vcat 2
@test vcat(3, 2) == @chainback 2 3 vcat

@test "See documentation of [`binaryfun`](@ref)" ==
    @chain (@doc @binaryfun) string chomp

```
"""
nonstandard(fs...) = Expr(:block, map(nonstandard1, fs)...)

export @nonstandard
eval(nonstandard(:nonstandard))


text:
# permute the columns until the principle subblock is nonsingular
# assume B approximates K x ∞
function nonsingular_permute(B)
    K = size(B,1)

    kr=[1]
    while rank(B[:,kr])<K
        if rank(B[:,kr])==length(kr)
            kr=[kr;kr[end]+1]  # add another col
        else
            kr[end]+=1  # integrate last col by 1
        end
    end

    kr=[kr;setdiff(1:size(B,2),kr)]
    eye(size(B,2))[:,kr]
end

## rearrange the rows and columns of bcs so that the principle block
# is identity
# assume we have the condition
#    BX=G
# where B is  K x n, X is n x m and G is K x *
#
# We also need to update L & M in
#    LX* + MX* = *
# to reflect the new ordering of rows of X
function regularize_bcs(B::Array)

    # permute rows of X/columns of B so principle block of B is nonsingular
    P = nonsingular_permute(B)

    B = B*P

    # we apply Q' to upper triangularize B,
    # avoiding any need to permute rows
    Q,R = qr(B)
    Q=Q[:,1:size(B,1)]

    K = size(B,1)

    # we invert the principle block of R
    # so that the BC leads with the identity
    Q = inv(R[:,1:K])*Q'
    R = inv(R[:,1:K])*R

    R,Q,P
end

function regularize_bcs{T<:Number}(B::Array, G::Matrix{T})
    R,Q,P=regularize_bcs(B)
    # we invert the principle block of R
    # so that the BC leads with the identity
    G = Q*G

    R,G,P
end

function regularize_bcs(B::Array, G::Array, L::Array, M::Array)
    if length(B) == 0
        R = B
        P = eye(size(L,2))
    else
        R,G,P=regularize_bcs(B,G)
        L = L*P
        M = M*P
    end

    R,G, L, M, P
end

# Result satisfies
#       inv(Q)*R*P' == B
#       L*P' == L (input)
#       M*P' == M (input)
function regularize_bcs(B::Array, L::Array, M::Array)
    if length(B) == 0
        R = B
        P = eye(size(L,2))
        Q= eye(0)
    else
        R,Q,P=regularize_bcs(B)

        L = L*P
        M = M*P
    end

    R,Q,L,M,P
end

function regularize_bcs(B::Array,Ls::Vector)
   if length(B) == 0
        R = B
        P = eye(size(B,2))
        Q= eye(0)
    else
        R,Q,P=regularize_bcs(B)
        Ls=[L*P for L in Ls]
    end
    R,Q,Ls,P
end

function reduce_dofs( R,G, Mx, My, F )
    if length(R) > 0
        GM = G*My.'
        for k = 1:size(R,1)
            F = F - Mx[:,k]*GM[k:k,:]
            Mx = Mx - Mx[:,k]*R[k:k,:]
        end
    end

    Mx, F
end

function reduce_constrained(Bx,Gx,By,Gy,Lx,Ly,Mx,My,F)
    Rx,Gx,Lx,Mx,Px=regularize_bcs(Bx,Gx,Lx,Mx)
    Ry,Gy,Ly,My,Py=regularize_bcs(By,Gy,Ly,My)

    Lx,F = reduce_dofs(Rx,Gx,Lx,Ly,F)
    Mx,F = reduce_dofs(Rx,Gx,Mx,My,F)
    Ly,F = reduce_dofs(Ry,Gy,Ly,Lx,F.');    F = F.';
    My,F = reduce_dofs(Ry,Gy,My,Mx,F.');    F = F.';

    Kx = size(Bx,1); Ky = size(By,1);
    A=Lx[:,Kx+1:end];B=Ly[:,Ky+1:end];
    C=Mx[:,Kx+1:end];D=My[:,Ky+1:end];

    A,B,C,D,F,Rx,Gx,Px,Ry,Gy,Py
end

constrained_lyap(B,L,M,F)=constrained_lyap(B[1,1],B[1,2],B[2,1],B[2,2],L[1],L[2],M[1],M[2],F)

function constrained_lyap(Bx,Gx,By,Gy,Lx,Ly,Mx,My,F)
    A,B,C,D,F,Rx,Gx,Px,Ry,Gy,Py = reduce_constrained(Bx,Gx,By,Gy,Lx,Ly,Mx,My,F)

    X22=lyap(A,B,C,D,F)

    Kx = size(Bx,1); Ky = size(By,1);
    X12 = Gx[:,Ky+1:end] - Rx[:,Kx+1:end]*X22
    X21 = Gy[:,Kx+1:end].' - X22*Ry[:,Ky+1:end].'
    X11 = Gx[:,1:Ky] - Rx[:,Kx+1:end]*X21
    X11a= Gy[:,1:Kx].' - X12*Ry[:,Ky+1:end].'

    tol = 1e-13;
    bcerr = norm(X11 - X11a);

    if bcerr>tol
       warn("Boundary conditions differ by " * string(bcerr));
    end

    X = [X11 X12; X21 X22];

    X = Px*X*Py.';

end


text:
callerror() = error("ObjectiveC call: use [obj method] or [obj method:param ...]")

function flatvcat(ex::Expr)
  any(ex->isexpr(ex, :row), ex.args) || return ex
  flat = Expr(:hcat)
  for row in ex.args
    isexpr(row, :row) ?
      push!(flat.args, row.args...) :
      push!(flat.args, row)
  end
  return calltransform(flat)
end

function calltransform(ex::Expr)
  obj = objcm(ex.args[1])
  args = ex.args[2:end]
  isempty(args) && callerror()
  if isexpr(args[1], Symbol)
    length(args) > 1 && callerror()
    return :($message($obj, $(Selector(args[1]))))
  end

  all(arg->isexpr(arg, :call) && isexpr(arg.args[1], Symbol), args) || callerror()
  msg = join(vcat([arg.args[2] for arg in args], ""), ":") |> Selector
  args = [objcm(arg.args[3]) for arg in args]
  :($message($obj, $msg, $(args...)))
end

objcm(ex::Expr) =
  isexpr(ex, :hcat) ? calltransform(ex) :
  isexpr(ex, :vcat) ? flatvcat(ex) :
  isexpr(ex, :block, :let) ? Expr(:block, map(objcm, ex.args)...) :
  esc(ex)

objcm(ex) = ex

macro objc(ex)
  esc(objcm(ex))
end

# Import Classes

macro classes(names)
  isexpr(names, Symbol) ? (names = [names]) : (names = names.args)
  Expr(:block, [:(const $(esc(name)) = Class($(Expr(:quote, name))))
                for name in names]..., nothing)
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# iis.jl
#
# Extract an Irreducible Inconsistent Subsystem (IIS)
# from an infeasible model.
# This code is specific to Gurobi.
#############################################################################

using JuMP, Gurobi
import MathProgBase

function print_iis_gurobi(m::Model)

    grb = getrawsolver(m)
    Gurobi.computeIIS(grb)
    numconstr = Gurobi.num_constrs(grb)
    numvar = Gurobi.num_vars(grb)

    iisconstr = Gurobi.get_intattrarray(grb, "IISConstr", 1, numconstr)
    iislb = Gurobi.get_intattrarray(grb, "IISLB", 1, numvar)
    iisub = Gurobi.get_intattrarray(grb, "IISUB", 1, numvar)

    println("Irreducible Inconsistent Subsystem (IIS)")
    println("Variable bounds:")
    for i in 1:numvar
        v = Variable(m, i)
        if iislb[i] != 0 && iisub[i] != 0
            println(getlowerbound(v), " <= ", getname(v), " <= ", getupperbound(v))
        elseif iislb[i] != 0
            println(getname(v), " >= ", getlowerbound(v))
        elseif iisub[i] != 0
            println(getname(v), " <= ", getupperbound(v))
        end
    end

    println("Constraints:")
    for i in 1:numconstr
        if iisconstr[i] != 0
            println(ConstraintRef{Model,LinearConstraint}(m, i))
        end
    end
    println("End of IIS")
end

# The function below produces the following output:
#=
Irreducible Inconsistent Subsystem (IIS)
Variable bounds:
x <= 1.0
y <= 1.0
Constraints:
x + y >= 3
End of IIS
=#

function iis_example()
    m = Model(solver=GurobiSolver())

    @variable(m, 0 <= x <= 1)
    @variable(m, 0 <= y <= 1)

    @constraint(m, x+y >= 3)

    status = solve(m)
    @assert status == :Infeasible

    print_iis_gurobi(m)
end


text:
#############################################################################
# customized_submod.jl
# allow users to customize submodular functions
# no check is done for the submodularity or cardinality of the function
#############################################################################

export submod
export sign, monotonicity, modularity, evaluate

mutable struct CustomizedSubmodAtom <: SubmodFunc
  head::Symbol
  id_hash::UInt64
  children::Tuple{Function, SetVariable}
  size::Tuple{Int, Int}
  func::Function
  setvariables::Array{CombiSet}

  function CustomizedSubmodAtom(F::Function, S::SetVariable)
    info("Please make sure the function is submodular, the function returns 0 at the empty set, and the cardinality of the base set is equal to that of the set variable.")
    children = (F, S)
    setvariables = get_sv(S)
    return new(:customized_submod, hash(children), children, (1, 1), F, setvariables)
  end
end

submod(F::Function, S::SetVariable) = CustomizedSubmodAtom(F, S)

function sign(F::CustomizedSubmodAtom)
  return NoSign()
end

function monotonicity(F::CustomizedSubmodAtom)
  return NoMonotinicity()
end

function modularity(F::CustomizedSubmodAtom)
  return SubModularity()
end

function evaluate(F::CustomizedSubmodAtom)
  return F.func(F.children[2].elements)
end


text:
"""
    Implements Thompson Sampling for Bandits with r ∈ [ 0,1]
"""

type TS <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    α0::Vector{Int64}
    β0::Vector{Int64}
    cummSuccess::Vector{Int64}
    cummFailure::Vector{Int64}

    samplingDist::Vector{Distributions.Beta}

    function TS( noOfArms::Integer )
        new( noOfArms,
             0,
             0,
             ones(Int64,noOfArms),
             ones(Int64,noOfArms),
             zeros(Float64,noOfArms),
             zeros(Float64,noOfArms),
             fill(Distributions.Beta(1,1),noOfArms)
        )
    end

    function TS( armParams::Array{Tuple{Int64,Int64},1} )
        _noOfArms   = length( armParams )
        _priorDist  = [ Distributions.Beta(armParams[idx][1],armParams[idx][2]) for idx=1:_noOfArms ]
        new( _noOfArms,
             0,
             0,
             [ armParams[idx][1] for idx=1:_noOfArms ],
             [ armParams[idx][2] for idx=1:_noOfArms ],
             zeros(Float64,_noOfArms),
             zeros(Float64,_noOfArms),
             _priorDist
        )
    end
end

function get_arm_index( agent::TS )
    agent.lastPlayedArm = findmax( [rand(armSample) for armSample in agent.samplingDist] )[2]
    return agent.lastPlayedArm
end

function update_reward!( agent::TS, r::Int64 )
    # Update S and F
    agent.cummSuccess[agent.lastPlayedArm] += (r==0?0:1)
    agent.cummFailure[agent.lastPlayedArm] += (r==0?1:0)

    # Update Distributions
    agent.samplingDist[agent.lastPlayedArm] = Distributions.Beta(
                                                agent.cummSuccess[agent.lastPlayedArm]+agent.α0[agent.lastPlayedArm],
                                                agent.cummFailure[agent.lastPlayedArm]+agent.β0[agent.lastPlayedArm]
                                            )

    # Update time steps
    agent.noOfSteps += 1

    nothing
end

function update_reward!( agent::TS, r::Float64 )
    rTilde = rand( Distributions.Bernoulli(r) )
    update_reward!( agent, rTilde )

    nothing
end

function reset!( agent::TS )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0

    agent.cummSuccess   = zeros( Float64, agent.noOfArms )
    agent.cummFailure   = zeros( Float64, agent.noOfArms )
    agent.samplingDist  = [ Distributions.Beta(agent.α0[idx],agent.β0[idx]) for idx=1:agent.noOfArms ]

    nothing
end

function info_str( agent::TS, latex::Bool )
    return @sprintf( "Thompson Sampling" )
end

"""-----------------------------------------------------------------------------
    Dynamic Thompson Sampling
    Based on: Gupta, N., Granmo, O. C., & Agrawala, A. (2011). Thompson sampling for dynamic multi-armed bandits. Proceedings - 10th International Conference on Machine Learning and Applications, ICMLA 2011, 1, 484–489. http://doi.org/10.1109/ICMLA.2011.144
"""

type DynamicTS <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    α::Vector{Float64}
    β::Vector{Float64}
    C::Float64

    samplingDist::Vector{Distributions.Beta}

    function DynamicTS( noOfArms, _C )
        new( noOfArms,
             0,
             0,
             2*ones(Float64,noOfArms),
             2*ones(Float64,noOfArms),
             _C,
             fill(Distributions.Beta(2,2),noOfArms)
        )
    end
end

function get_arm_index( agent::DynamicTS )
    agent.lastPlayedArm = findmax( [rand(armSample) for armSample in agent.samplingDist] )[2]
    return agent.lastPlayedArm
end

function update_reward!( agent::DynamicTS, r::Integer )
    # Update reward to arm played
    if agent.α[agent.lastPlayedArm]+agent.β[agent.lastPlayedArm] < agent.C
        agent.α[agent.lastPlayedArm]    += (r==0?0:1)
        agent.β[agent.lastPlayedArm]    += (r==0?1:0)
    else
        agent.α[agent.lastPlayedArm]    = (agent.α[agent.lastPlayedArm]+(r==0?0:1)) * agent.C/(agent.C+1)
        agent.β[agent.lastPlayedArm]    = (agent.β[agent.lastPlayedArm]+(r==0?1:0)) * agent.C/(agent.C+1)
    end

    # Update arm's distribution
    agent.samplingDist[agent.lastPlayedArm]  = Distributions.Beta(
                                                agent.α[agent.lastPlayedArm],
                                                agent.β[agent.lastPlayedArm]
                                            )

    # Update time steps
    agent.noOfSteps += 1

    nothing
end

function update_reward!( agent::DynamicTS, r::AbstractFloat )
    update_reward!( agent, rand(Distributions.Bernoulli(r)) )   # Do a Bernoulli Trial to update the posterior

    nothing
end

function reset!( agent::DynamicTS )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0

    agent.α             = 2 * ones( Float64, agent.noOfArms )
    agent.β             = 2 * ones( Float64, agent.noOfArms )
    agent.samplingDist  = fill( Distributions.Beta(2,2), agent.noOfArms )

    nothing
end

function info_str( agent::DynamicTS, latex::Bool )
    if latex
        return @sprintf( "Dynamic Thompson Sampling (\$C = %4.3f\$)", agent.C )
    else
        return @sprintf( "Dynamic Thompson Sampling (C = %4.3f)", agent.C )
    end
end

"""-----------------------------------------------------------------------------
    Optimistic Thompson Sampling
    Based on: May, B. C., Korda, N., Lee, A., & Leslie, D. S. (2012). Optimistic bayesian sampling in contextual-bandit problems. Journal of Machine Learning Research, 13, 2069–2106.
"""

type OTS <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    α0::Vector{Int64}
    β0::Vector{Int64}
    cummSuccess::Vector{Int64}
    cummFailure::Vector{Int64}

    samplingDist::Vector{Distributions.Beta}

    function OTS( noOfArms::Integer )
        new( noOfArms,
             0,
             0,
             ones(Int64,noOfArms),
             ones(Int64,noOfArms),
             zeros(Float64,noOfArms),
             zeros(Float64,noOfArms),
             fill(Distributions.Beta(1,1),noOfArms)
        )
    end

    function OTS( armParams::Array{Tuple{Int64,Int64},1} )
        _noOfArms   = length( armParams )
        _priorDist  = [ Distributions.Beta(armParams[idx][1],armParams[idx][2]) for idx=1:_noOfArms ]
        new( _noOfArms,
             0,
             0,
             [ armParams[idx][1] for idx=1:_noOfArms ],
             [ armParams[idx][2] for idx=1:_noOfArms ],
             zeros(Float64,_noOfArms),
             zeros(Float64,_noOfArms),
             _priorDist
        )
    end
end

function get_arm_index( agent::OTS )
    agent.lastPlayedArm = findmax( [ max(mean(armSample),rand(armSample)) for armSample in agent.samplingDist] )[2]
    return agent.lastPlayedArm
end

function update_reward!( agent::OTS, r::Int64 )
    # Update S and F
    agent.cummSuccess[agent.lastPlayedArm] += (r==1?1:0)
    agent.cummFailure[agent.lastPlayedArm] += (r==0?1:0)

    # Update Distributions
    agent.samplingDist[agent.lastPlayedArm] = Distributions.Beta(
                                                agent.cummSuccess[agent.lastPlayedArm]+agent.α0[agent.lastPlayedArm],
                                                agent.cummFailure[agent.lastPlayedArm]+agent.β0[agent.lastPlayedArm]
                                            )
    # Update time steps
    agent.noOfSteps += 1

    nothing
end

function update_reward!( agent::OTS, r::Float64 )
    rTilde = rand( Distributions.Bernoulli(r) )
    update_reward!( agent, rTilde )

    nothing
end

function reset!( agent::OTS )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0

    agent.cummSuccess   = zeros( Float64, agent.noOfArms )
    agent.cummFailure   = zeros( Float64, agent.noOfArms )
    agent.samplingDist  = fill( Distributions.Beta(1,1), agent.noOfArms )

    nothing
end

function info_str( agent::OTS, latex::Bool )
    return @sprintf( "Optimistic Thompson Sampling" )
end

"""-----------------------------------------------------------------------------
    Thompson Sampling for Normal Rewards
    Based on: Sec 6.3, Russo, D., & Van Roy, B. (2014). Learning to Optimize Via Posterior Sampling. Mathematics of Operations Research. http://doi.org/10.1287/xxxx.0000.0000
"""

type TSNormal <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    η::Float64
    μ::Vector{Float64}
    σ::Vector{Float64}

    samplingDist::Vector{Distributions.Normal}

    function TSNormal( noOfArms::Int64, η::Real = 1 )
        new(
            noOfArms,
            0,
            0,
            η,
            zeros(Float64,noOfArms),
            100*ones(Float64,noOfArms),    # σ - Some high value
            fill(Distributions.Normal(0,100),noOfArms)
        )
    end

    function TSNormal( armParams::Array{Tuple{Real,Real},1}, η::Real = 1 )
        _K = length(armParams)
        new(
            _K,
            0,
            0,
            η,
            [ armParams[idx][1] for i = 1:_K ],
            [ armParams[idx][2] for i = 1:_K ],
            [ Distributions.Normal(armParams[idx][1],armParams[idx][2]) for i = 1:_K ]
        )
    end
end

function get_arm_index( agent::TSNormal )
    agent.lastPlayedArm = findmax( [ max(mean(armSample),rand(armSample)) for armSample in agent.samplingDist] )[2]
    return agent.lastPlayedArm
end

function update_reward!( agent::TSNormal, r::Float64 )

    agent.μ[agent.lastPlayedArm] = ( agent.μ[agent.lastPlayedArm]/agent.σ[agent.lastPlayedArm]^2 + r/agent.η^2 ) /
                                        ( 1/agent.σ[agent.lastPlayedArm]^2 + 1/agent.η^2 )
    agent.σ[agent.lastPlayedArm] = ( 1/agent.σ[agent.lastPlayedArm]^2 + 1/agent.η^2 )^(-1/2)

    # Update Distributions
    agent.samplingDist[agent.lastPlayedArm] = Distributions.Normal(
                                                agent.μ[agent.lastPlayedArm],
                                                agent.σ[agent.lastPlayedArm]
                                            )

    agent.noOfSteps += 1

    nothing
end

function reset!( agent::TSNormal )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0

    agent.μ   = zeros( Float64, agent.noOfArms )
    agent.σ   = 100*ones(Float64,agent.noOfArms)
    agent.samplingDist  = [ Distributions.Normal(agent.μ[i],agent.σ[i]) for i = 1:agent.noOfArms ]

    nothing
end

function info_str( agent::TSNormal, latex::Bool )
    if latex
        return @sprintf( "TS-Normal(\$\\eta=%3.2f\$)", agent.η )
    else
        return @sprintf( "TS-Normal(η=%3.2f)", agent.η )
    end
end

"""-----------------------------------------------------------------------------
    Discounted Thompson Sampling
    Based on: Vishnu Raj, Sheetal Kalyani (2017). Taming Non-stationary Bandits: A Bayesian Approach. https://arxiv.org/abs/1707.09727
"""

type dTS <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    γ::Float64                      # Discounting factor
    αₒ::Vector{Float64}             # Initial α value
    βₒ::Vector{Float64}             # Initial β value
    cummSuccess::Vector{Float64}    # Cummulative success
    cummFailure::Vector{Float64}    # Cummulative failure

    samplingDist::Vector{Distributions.Beta}    # Internal Distributions

    function dTS( noOfArms, γ::Float64 = 1.00,
                    αₒ::Float64 = 1.00, βₒ::Float64 = 1.00 )
        new( noOfArms,
             0,
             0,
             γ,
             αₒ * ones(Float64,noOfArms),
             βₒ * ones(Float64,noOfArms),
             αₒ * ones(Float64,noOfArms),
             βₒ * ones(Float64,noOfArms),
             fill(Distributions.Beta(αₒ,βₒ),noOfArms)
        )
    end

    function dTS( noOfArms, γ::Float64 ,
                    α₀::Vector{Float64}, β₀::Vector{Float64} )
        new( noOfArms,
             0,
             0,
             γ,
             α₀,
             β₀,
             α₀,
             β₀,
             Distributions.Beta.(α₀,β₀)
        )
    end
end

function get_arm_index( agent::dTS )
    agent.lastPlayedArm = findmax( [rand(armSample) for armSample in agent.samplingDist] )[2]
    return agent.lastPlayedArm
end

function update_reward!( agent::dTS, r::Int64 )
    # Update S and F
    agent.cummSuccess *= agent.γ
    agent.cummFailure *= agent.γ
    agent.cummSuccess[agent.lastPlayedArm] += (r==0?0:1)
    agent.cummFailure[agent.lastPlayedArm] += (r==0?1:0)

    # Update Distributions
    agent.samplingDist  = Distributions.Beta.(
                                                agent.cummSuccess + agent.αₒ,
                                                agent.cummFailure + agent.βₒ
                                            )

    # Update time steps
    agent.noOfSteps += 1
end

function update_reward!( agent::dTS, r::Float64 )
    rTilde = rand( Distributions.Bernoulli(r) )
    update_reward!( agent, rTilde )
end

function reset!( agent::dTS )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0

    agent.cummSuccess   = copy( agent.αₒ )
    agent.cummFailure   = copy( agent.βₒ )
    agent.samplingDist  = Distributions.Beta.(agent.αₒ,agent.βₒ)
end

function info_str( agent::dTS, latex::Bool )
    if latex
        return @sprintf( "dTS(\$\\gamma = %4.3f\$)", agent.γ )
    else
        return @sprintf( "dTS(γ = %4.3f)", agent.γ )
    end
end

""" ----------------------------------------------------------------------------
    Discounted Optimistic Thompson Sampling
    Based on: Vishnu Raj, Sheetal Kalyani (2017). Taming Non-stationary Bandits: A Bayesian Approach. https://arxiv.org/abs/1707.09727
"""
type dOTS <: BanditAlgorithmBase
    _dTS::dTS

    function dOTS( noOfArms, γ::Float64 = 1.00,
                    αₒ::Float64 = 1.00, βₒ::Float64 = 1.00 )
        new(
             dTS( noOfArms, γ, αₒ, βₒ )
        )
    end
end

function get_arm_index( agent::dOTS )
    agent._dTS.lastPlayedArm = findmax( [max(rand(armSample),mean(armSample)) for armSample in agent._dTS.samplingDist] )[2]
    return agent._dTS.lastPlayedArm
end

function update_reward!( agent::dOTS, r::Int64 )
    update_reward!( agent._dTS, r )
    nothing
end

function update_reward!( agent::dOTS, r::Float64 )
    update_reward!( agent._dTS, r )
    nothing
end

function reset!( agent::dOTS )
    reset!( agent._dTS )
    nothing
end

function info_str( agent::dOTS, latex::Bool )
    if latex
        return @sprintf( "dOTS(\$\\gamma = %4.3f\$)", agent._dTS.γ )
    else
        return @sprintf( "dOTS(γ = %4.3f)", agent._dTS.γ )
    end
end

""" ----------------------------------------------------------------------------
    Restarting Thompson Sampling
    Based on:
"""

type RestartTS <: BanditAlgorithmBase
    _TS::TS
    Δ::Int64

    function RestartTS( noOfArms::Int, Δ )
        new( TS(noOfArms), Δ )
    end
end

function get_arm_index( agent::RestartTS )
    return get_arm_index( agent._TS )
end

function update_reward!( agent::RestartTS, r::Int64 )
    update_reward!( agent._TS, r )
    if agent._TS.noOfSteps >= agent.Δ
        reset!( agent._TS )
    end
    nothing
end

function update_reward!( agent::RestartTS, r::Float64 )
    update_reward!( agent._TS, r )
    if agent._TS.noOfSteps >= agent.Δ
        reset!( agent._TS )
    end
    nothing
end

function reset!( agent::RestartTS )
    reset!( agent._TS )
    nothing
end

function info_str( agent::RestartTS, latex::Bool )
    if latex
        return @sprintf( "Restarting TS(\$\\Delta = %d\$)", agent.Δ )
    else
        return @sprintf( "Restarting TS(Δ = %d)", agent.Δ )
    end
end

""" ----------------------------------------------------------------------------------------
    Thompson Sampling with Gaussian Priors
    Based on: Algorithm 2, S. Agrawal and N. Goyal, “Near-Optimal Regret Bounds for Thompson Sampling,” J. ACM, vol. 64, no. 5, 2017.
"""
type TSGaussPrior <: BanditAlgorithmBase
    noOfArms::Integer
    noOfSteps::Integer
    lastPlayedArm::Integer

    # Algorithm specific parameters
    μ::Vector{Float64}
    noOfPulls::Vector{Float64} # Even though an integer, keeping this as Float for performance

    samplingDist::Vector{Distributions.Normal}

    function TSGaussPrior( noOfArms::Integer )
        new(
            noOfArms,
            0,
            0,
            zeros( Float64, noOfArms),
            zeros( Float64, noOfArms ),
            fill( Distributions.Normal(0,1), noOfArms )
        )
    end
end

function get_arm_index( agent::TSGaussPrior )
    agent.lastPlayedArm = findmax( [rand(armSample) for armSample in agent.samplingDist] )[2]
    return agent.lastPlayedArm
end

function  update_reward!( agent::TSGaussPrior, r::Real )
    # Update observation
    agent.μ[agent.lastPlayedArm] = (agent.μ[agent.lastPlayedArm]*agent.noOfPulls[agent.lastPlayedArm]+r) /
                                        (agent.noOfPulls[agent.lastPlayedArm]+1)
    agent.noOfPulls[agent.lastPlayedArm] += 1

    # Update distribution
    agent.samplingDist[agent.lastPlayedArm] = Distributions.Normal(
                                                        agent.μ[agent.lastPlayedArm],
                                                        1.0/(agent.noOfPulls[agent.lastPlayedArm]+1)
                                                )

    # Update time steps
    agent.noOfSteps += 1

    nothing
end

function  reset!( agent::TSGaussPrior )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0

    agent.μ             = zeros( Float64, agent.noOfArms )
    agent.noOfPulls     = zeros( Float64, agent.noOfArms )

    agent.samplingDist  = fill( Distributions.Normal(0,1), agent.noOfArms )
end

function  info_str( agent::TSGaussPrior, latex::Bool )
    return @sprintf( "TS Gauss Prior" )
end


text:
# Suitable for any linear transformation on any type of multi-dimensional arrays.
# Implementation details and advantages:
# * has the option to keep all previous residual vectors for orthogonalizing the newest
#   residual vectors obtained.
#
# * The re-orthogonalization process is computed in parallel because it uses a matrix of zero that
#   resizes to double the number of columns each time.
#
# * (TO ADD) Always store the Tridiaognalizations of the Linear operator while running the conjugate gradient.

mutable struct ConjGradModified{T <: Number}
    A::Function              # Linear opeartor
    b::AbstractArray{T}      # RHS vector
    tensor_size::Tuple       # The size of the tensor the linear operator is acting on.

    x::Vector{T}            # current guess, started with initial guess.
    r::Vector{T}            # previous computed residual
    rnew::Vector{T}         # Current residual
    d::Vector{T}            # Conjugate Direction
    Ad::Vector{T}           # for reducing garbage collector time.

    itr::UInt64              # Iteration count.
    Q::Union{AbstractMatrix{T}, Nothing}        # Re-Orthogonalization basis.
    Q_size::UInt64           # Number of orthogonalization vectors.
    over_write::UInt64       # position for overwrite once storage is limit is reached.

    storage_limit::UInt64    # storage limit for the Q vector.
    reorthogonalize::Bool    # Whether to perform reorthogonalization.

    function ConjGradModified(
        A::Function,
        b::AbstractArray,
        x0::Union{AbstractArray,Nothing}=nothing
    )
        x0 = x0===nothing ? b .+ 0.1 : x0
        r = reshape(b - A(x0), :)
        this = new{typeof(r[1])}()
        this.r = r
        this.A = A
        this.b = b
        this.tensor_size = size(b)
        this.x = x0
        # this.r = ComputeResidualVec(this, this.x)
        this.rnew = similar(this.r)
        this.d = this.r
        this.Ad = similar(this.r)
        this.x = x0
        this.itr = 0

        this.Q = zeros(typeof(this.r[1]), length(this.r), 4)
        this.Q[:, 1] = this.r/norm(this.r)
        this.Q_size = 1
        this.over_write = 0

        this.storage_limit = length(this.r) - 1       # Maximal limit, if not, we have a problem.
        this.reorthogonalize = true
        return this
    end

    function ConjGradModified(A::AbstractMatrix, b::AbstractVector, x0::Union{AbstractArray,Nothing}=nothing)
        return ConjGradModified((x)-> A*x, b, x0)
    end

end

"""
    Turn on the reorthogonalizations using the residual vectors,
    and add current residual to the list of residuals.
"""
function TurnOnReorthgonalize(this::ConjGradModified)
    this.reorthogonalize = true
    this.Q = zeros(typeof(this.r[1]), length(this.r), storage_limit)
    this.Q[:, 1] = this.r/norm(this.r)
    this.Q_size = 1
    this.over_write = 0
    return
end

"""
    Turn off the reorthogonalization on the residual vectors. And then
    clear all the stored residual vectors.
"""
function TurnOffReorthgonalize(this::ConjGradModified)
    this.reorthogonalize = false
    this.Q = nothing
    return
end

"""
    Change the storage limit for the number of vectors used for
    re-orthogonalization.
"""
function ChangeStorageLimit(this::ConjGradModified, storage_limit)
    this.storage_limit = min(length(this.r) - 1, storage_limit)
    return
end

"""
    Set the number of vectors in the Q matrix by force. This will clear
    everything that is already in Q.
"""
function SetStorageLimit(this::ConjGradModified, storage_limit)
    this.Q = zeros(typeof(this.r[1]), length(this.r), min(length(this.r) - 1, storage_limit))
    this.Q[:, 1] = this.r/norm(this.r)
    this.storage_limit = storage_limit
    this.Q_size = 1
    this.over_write = 0
end

"""
    Apply the linear transformation on the tensor and return the
    residual vector after transformation.
"""
function ComputeResidualVec(this::ConjGradModified, x::AbstractArray)
    resVec = this.b - this.A(reshape(x, this.tensor_size))
    return reshape(resVec, length(resVec))
end

"""
    Get the norm of the current residual of Modified Conjugate
    Gradient.
"""
function GetResidualNorm(this::ConjGradModified)
    return norm(this.r)
end

"""
    Compute the vector using the linear transformation.
"""
function ComputeVec(this::ConjGradModified, x::AbstractArray)
    vec = this.A(reshape(x, this.tensor_size))
    return reshape(vec, length(vec))
end

"""
    Pre-allocates a vector and get the values to the given
    pre-allocated vector.
"""
function ComputeVec!(this::ConjGradModified, x::AbstractArray, vec::AbstractArray)
    vec .= this.A(reshape(x, this.tensor_size))
    return reshape(vec, :)
end

"""
    Performs One step of conjugate Gradient.
"""
function (this::ConjGradModified)()
    r = this.r
    if norm(r) == 0
        return 0 # The problem is solved already.
    end
    d = this.d
    Ad = this.Ad
    ComputeVec!(this, d, Ad)
    a = dot(r, r)/dot(d, Ad)

    if a < 0
        error("CG got a non-definite matrix")
    end

    this.x += a*d
    this.rnew .= r - a*Ad                    # update rnew
    rnewNorm = norm(this.rnew)

    if this.reorthogonalize
        this.rnew .-= this.Q*(this.Q'*this.rnew)

        if this.Q_size == this.storage_limit  # starts overwrite
            this.Q[:, this.over_write + 1] = this.rnew/rnewNorm
            this.over_write = (this.over_write + 1)%this.storage_limit
        elseif this.Q_size == size(this.Q, 2) # Resize
            newQ = zeros(
                typeof(this.r[1]),
                size(this.Q, 1),
                min(2*size(this.Q, 2), this.storage_limit)
            )
            newQ[:, 1:this.Q_size] = this.Q
            newQ[:, this.Q_size + 1] = this.rnew/rnewNorm
            this.Q = newQ
            this.Q_size += 1
        else                                  # add one more.
            this.Q[:, this.Q_size + 1] = this.rnew/rnewNorm
            this.Q_size += 1
        end
    end

    b = dot(this.rnew, this.rnew)/dot(r, r)
    # @assert abs(dot(rnew + β*d, Ad)) < 1e-8 "Not conjugate"
    this.d = this.rnew + b*d
    this.r = copy(this.rnew)                      # Override
    this.itr += 1
    return convert(Float64, rnewNorm)
end

"""
    Performs j iterations of conjugate gradient, collect the
    2norm of the residual and return all of them in a vector.
"""
function (this::ConjGradModified)(j::Int64)
    results = Vector()
    for _ in 1: j
        push!(results, this())
    end
    return results
end


text:
# LLVM IR generation

## method compilation tracer

# this functionality is used to detect recursion, and functions that shouldn't be called.
# it is a hack, and should disappear over time. don't add new features to it.

# generate a pseudo-backtrace from a stack of methods being emitted
function backtrace(job::CompilerJob, call_stack::Vector{Core.MethodInstance})
    bt = StackTraces.StackFrame[]
    for method_instance in call_stack
        method = method_instance.def
        if method.name === :overdub && isdefined(method, :generator)
            # The inline frames are maintained by the dwarf based backtrace, but here we only have the
            # calls to overdub directly, the backtrace therefore is collapsed and we have to
            # lookup the overdubbed function, but only if we likely are using the generated variant.
            actual_sig = Tuple{method_instance.specTypes.parameters[3:end]...}
            m = ccall(:jl_gf_invoke_lookup, Any, (Any, UInt), actual_sig, typemax(UInt))
            method = m.func::Method
        end
        frame = StackTraces.StackFrame(method.name, method.file, method.line)
        pushfirst!(bt, frame)
    end
    bt
end

# NOTE: we use an exception to be able to display a stack trace using the logging framework
struct MethodSubstitutionWarning <: Exception
    original::Method
    substitute::Method
end
Base.showerror(io::IO, err::MethodSubstitutionWarning) =
    print(io, "You called $(err.original), maybe you intended to call $(err.substitute) instead?")
const method_substitution_whitelist = [:hypot, :exp]

mutable struct MethodCompileTracer
    job::CompilerJob
    call_stack::Vector{Core.MethodInstance}
    last_method_instance::Union{Nothing,Core.MethodInstance}

    MethodCompileTracer(job, start) = new(job, Core.MethodInstance[start])
    MethodCompileTracer(job) = new(job, Core.MethodInstance[])
end

function Base.push!(tracer::MethodCompileTracer, method_instance)
    push!(tracer.call_stack, method_instance)

    if VERSION < v"1.5.0-DEV.393"
        # check for recursion
        if method_instance in tracer.call_stack[1:end-1]
            throw(KernelError(tracer.job, "recursion is currently not supported";
                              bt=backtrace(tracer.job, tracer.call_stack)))
        end
    end

    # check for Base functions that exist in the GPU package
    # FIXME: this might be too coarse
    method = method_instance.def
    if Base.moduleroot(method.module) == Base &&
        isdefined(runtime_module(tracer.job), method_instance.def.name) &&
        !in(method_instance.def.name, method_substitution_whitelist)
        substitute_function = getfield(runtime_module(tracer.job), method.name)
        tt = Tuple{method_instance.specTypes.parameters[2:end]...}
        if hasmethod(substitute_function, tt)
            method′ = which(substitute_function, tt)
            if method′.module == runtime_module(tracer.job)
                @warn "calls to Base intrinsics might be GPU incompatible" exception=(MethodSubstitutionWarning(method, method′), backtrace(tracer.job, tracer.call_stack))
            end
        end
    end
end

function Base.pop!(tracer::MethodCompileTracer, method_instance)
    @compiler_assert last(tracer.call_stack) == method_instance tracer.job
    tracer.last_method_instance = pop!(tracer.call_stack)
end

Base.last(tracer::MethodCompileTracer) = tracer.last_method_instance

## Julia compiler integration

if VERSION >= v"1.5.0-DEV.393"

# JuliaLang/julia#25984 significantly restructured the compiler

function compile_method_instance(job::CompilerJob, method_instance::Core.MethodInstance, world)
    # set-up the compiler interface
    tracer = MethodCompileTracer(job, method_instance)
    hook_emit_function(method_instance, code) = push!(tracer, method_instance)
    hook_emitted_function(method_instance, code) = pop!(tracer, method_instance)
    param_kwargs = [:track_allocations  => false,
                    :code_coverage      => false,
                    :static_alloc       => false,
                    :prefer_specsig     => true,
                    :emit_function      => hook_emit_function,
                    :emitted_function   => hook_emitted_function]
    if LLVM.version() >= v"8.0" && VERSION >= v"1.3.0-DEV.547"
        push!(param_kwargs, :gnu_pubnames => false)

        debug_info_kind = if Base.JLOptions().debug_level == 0
            LLVM.API.LLVMDebugEmissionKindNoDebug
        elseif Base.JLOptions().debug_level == 1
            LLVM.API.LLVMDebugEmissionKindLineTablesOnly
        elseif Base.JLOptions().debug_level >= 2
            LLVM.API.LLVMDebugEmissionKindFullDebug
        end

        # LLVM's debug info crashes older CUDA assemblers
        if job.target isa PTXCompilerTarget # && driver_version(job.target) < v"10.2"
            # FIXME: this was supposed to be fixed on 10.2
            @debug "Incompatibility detected between CUDA and LLVM 8.0+; disabling debug info emission" maxlog=1
            debug_info_kind = LLVM.API.LLVMDebugEmissionKindNoDebug
        end

        push!(param_kwargs, :debug_info_kind => Cint(debug_info_kind))
    end
    params = Base.CodegenParams(;param_kwargs...)

    # generate IR
    if VERSION >= v"1.5.0-DEV.851"
        native_code = ccall(:jl_create_native, Ptr{Cvoid},
                            (Vector{Core.MethodInstance}, Base.CodegenParams, Cint),
                            [method_instance], params, #=extern policy=# 1)
    else
        native_code = ccall(:jl_create_native, Ptr{Cvoid},
                            (Vector{Core.MethodInstance}, Base.CodegenParams),
                            [method_instance], params)
    end
    @assert native_code != C_NULL
    llvm_mod_ref = ccall(:jl_get_llvm_module, LLVM.API.LLVMModuleRef,
                         (Ptr{Cvoid},), native_code)
    @assert llvm_mod_ref != C_NULL
    llvm_mod = LLVM.Module(llvm_mod_ref)

    # get the top-level code
    code = if VERSION >= v"1.6.0-DEV.12"
        # TODO: use our own interpreter
        interpreter = Core.Compiler.NativeInterpreter(world)
        Core.Compiler.inf_for_methodinstance(interpreter, method_instance, world, world)
    else
        Core.Compiler.inf_for_methodinstance(method_instance, world, world)
    end

    # get the top-level function index
    llvm_func_idx = Ref{Int32}(-1)
    llvm_specfunc_idx = Ref{Int32}(-1)
    ccall(:jl_breakpoint, Nothing, ())
    ccall(:jl_get_function_id, Nothing,
          (Ptr{Cvoid}, Any, Ptr{Int32}, Ptr{Int32}),
          native_code, code, llvm_func_idx, llvm_specfunc_idx)
    @assert llvm_func_idx[] != -1
    @assert llvm_specfunc_idx[] != -1

    # get the top-level function)
    llvm_func_ref = ccall(:jl_get_llvm_function, LLVM.API.LLVMValueRef,
                     (Ptr{Cvoid}, UInt32), native_code, llvm_func_idx[]-1)
    @assert llvm_func_ref != C_NULL
    llvm_func = LLVM.Function(llvm_func_ref)
    llvm_specfunc_ref = ccall(:jl_get_llvm_function, LLVM.API.LLVMValueRef,
                         (Ptr{Cvoid}, UInt32), native_code, llvm_specfunc_idx[]-1)
    @assert llvm_specfunc_ref != C_NULL
    llvm_specfunc = LLVM.Function(llvm_specfunc_ref)

    # configure the module
    triple!(llvm_mod, llvm_triple(job.target))
    if llvm_datalayout(job.target) !== nothing
        datalayout!(llvm_mod, llvm_datalayout(job.target))
    end

    return llvm_specfunc, llvm_mod
end

else

function module_setup(job::CompilerJob, mod::LLVM.Module)
    ctx = context(mod)

    # configure the module
    triple!(mod, llvm_triple(job.target))
    datalayout!(mod, llvm_datalayout(job.target))

    # add debug info metadata
    if LLVM.version() >= v"8.0"
        # Set Dwarf Version to 2, the DI printer will downgrade to v2 automatically,
        # but this is technically correct and the only version supported by NVPTX
        LLVM.flags(mod)["Dwarf Version", LLVM.API.LLVMModuleFlagBehaviorWarning] =
            Metadata(ConstantInt(Int32(2), ctx))
        LLVM.flags(mod)["Debug Info Version", LLVM.API.LLVMModuleFlagBehaviorError] =
            Metadata(ConstantInt(DEBUG_METADATA_VERSION(), ctx))
    else
        push!(metadata(mod), "llvm.module.flags",
             MDNode([ConstantInt(Int32(1), ctx),    # llvm::Module::Error
                     MDString("Debug Info Version"),
                     ConstantInt(DEBUG_METADATA_VERSION(), ctx)]))
    end
end

function compile_method_instance(job::CompilerJob, method_instance::Core.MethodInstance, world)
    function postprocess(ir)
        # get rid of jfptr wrappers
        for llvmf in functions(ir)
            startswith(LLVM.name(llvmf), "jfptr_") && unsafe_delete!(ir, llvmf)
        end

        return
    end

    # set-up the compiler interface
    tracer = MethodCompileTracer(job)
    hook_emit_function(method_instance, code, world) = push!(tracer, method_instance)
    hook_emitted_function(method_instance, code, world) = pop!(tracer, method_instance)
    dependencies = MultiDict{Core.MethodInstance,LLVM.Function}()
    function hook_module_setup(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        ir = LLVM.Module(ref)
        module_setup(job, ir)
    end
    function hook_module_activation(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        ir = LLVM.Module(ref)
        postprocess(ir)

        # find the function that this module defines
        llvmfs = filter(llvmf -> !isdeclaration(llvmf) &&
                                 linkage(llvmf) == LLVM.API.LLVMExternalLinkage,
                        collect(functions(ir)))

        llvmf = nothing
        if length(llvmfs) == 1
            llvmf = first(llvmfs)
        elseif length(llvmfs) > 1
            llvmfs = filter!(llvmf -> startswith(LLVM.name(llvmf), "julia_"), llvmfs)
            if length(llvmfs) == 1
                llvmf = first(llvmfs)
            end
        end

        @compiler_assert llvmf !== nothing job

        insert!(dependencies, last(tracer), llvmf)
    end
    param_kwargs = [:cached             => false,
                    :track_allocations  => false,
                    :code_coverage      => false,
                    :static_alloc       => false,
                    :prefer_specsig     => true,
                    :module_setup       => hook_module_setup,
                    :module_activation  => hook_module_activation,
                    :emit_function      => hook_emit_function,
                    :emitted_function   => hook_emitted_function]
    if LLVM.version() >= v"8.0" && VERSION >= v"1.3.0-DEV.547"
        push!(param_kwargs, :gnu_pubnames => false)

        debug_info_kind = if Base.JLOptions().debug_level == 0
            LLVM.API.LLVMDebugEmissionKindNoDebug
        elseif Base.JLOptions().debug_level == 1
            LLVM.API.LLVMDebugEmissionKindLineTablesOnly
        elseif Base.JLOptions().debug_level >= 2
            LLVM.API.LLVMDebugEmissionKindFullDebug
        end

        # LLVM's debug info crashes older CUDA assemblers
        if job.target isa PTXCompilerTarget # && driver_version(job.target) < v"10.2"
            # FIXME: this was supposed to be fixed on 10.2
            @debug "Incompatibility detected between CUDA and LLVM 8.0+; disabling debug info emission" maxlog=1
            debug_info_kind = LLVM.API.LLVMDebugEmissionKindNoDebug
        end

        push!(param_kwargs, :debug_info_kind => Cint(debug_info_kind))
    end
    params = Base.CodegenParams(;param_kwargs...)

    # get the code
    ref = ccall(:jl_get_llvmf_defn, LLVM.API.LLVMValueRef,
                (Any, UInt, Bool, Bool, Base.CodegenParams),
                method_instance, world, #=wrapper=#false, #=optimize=#false, params)
    if ref == C_NULL
        throw(InternalCompilerError(job, "the Julia compiler could not generate LLVM IR"))
    end
    llvmf = LLVM.Function(ref)
    ir = LLVM.parent(llvmf)
    postprocess(ir)

    # link in dependent modules
    entry = llvmf
    mod = LLVM.parent(entry)
    @timeit_debug to "linking" begin
        # we disable Julia's compilation cache not to poison it with GPU-specific code.
        # as a result, we might get multiple modules for a single method instance.
        cache = Dict{String,String}()

        for called_method_instance in keys(dependencies)
            llvmfs = dependencies[called_method_instance]

            # link the first module
            llvmf = popfirst!(llvmfs)
            llvmfn = LLVM.name(llvmf)
            link!(mod, LLVM.parent(llvmf))

            # process subsequent duplicate modules
            for dup_llvmf in llvmfs
                if Base.JLOptions().debug_level >= 2
                    # link them too, to ensure accurate backtrace reconstruction
                    link!(mod, LLVM.parent(dup_llvmf))
                else
                    # don't link them, but note the called function name in a cache
                    dup_llvmfn = LLVM.name(dup_llvmf)
                    cache[dup_llvmfn] = llvmfn
                end
            end
        end

        # resolve function declarations with cached entries
        for llvmf in filter(isdeclaration, collect(functions(mod)))
            llvmfn = LLVM.name(llvmf)
            if haskey(cache, llvmfn)
                def_llvmfn = cache[llvmfn]
                replace_uses!(llvmf, functions(mod)[def_llvmfn])

                @compiler_assert isempty(uses(llvmf)) job
                unsafe_delete!(LLVM.parent(llvmf), llvmf)
            end
        end
    end

    return entry, mod
end

end

function irgen(job::CompilerJob, method_instance::Core.MethodInstance, world)
    entry, mod = @timeit_debug to "emission" compile_method_instance(job, method_instance, world)
    ctx = context(mod)

    # clean up incompatibilities
    @timeit_debug to "clean-up" begin
        for llvmf in functions(mod)
            # only occurs in debug builds
            delete!(function_attributes(llvmf), EnumAttribute("sspstrong", 0, ctx))

            if VERSION < v"1.5.0-DEV.393"
                # make function names safe for ptxas
                llvmfn = LLVM.name(llvmf)
                if !isdeclaration(llvmf)
                    llvmfn′ = safe_name(llvmfn)
                    if llvmfn != llvmfn′
                        LLVM.name!(llvmf, llvmfn′)
                        llvmfn = llvmfn′
                    end
                end
            end

            if Sys.iswindows()
                personality!(llvmf, nothing)
            end
        end

        # remove the exception-handling personality function
        if Sys.iswindows() && "__julia_personality" in functions(mod)
            llvmf = functions(mod)["__julia_personality"]
            @compiler_assert isempty(uses(llvmf)) job
            unsafe_delete!(mod, llvmf)
        end
    end

    # target-specific processing
    process_module!(job, mod)

    # sanitize function names
    # FIXME: Julia should do this, but apparently fails (see maleadt/LLVM.jl#201)
    for f in functions(mod)
        LLVM.isintrinsic(f) && continue
        llvmfn = LLVM.name(f)
        startswith(llvmfn, "julia.") && continue # Julia intrinsics
        startswith(llvmfn, "llvm.") && continue # unofficial LLVM intrinsics
        llvmfn′ = safe_name(llvmfn)
        if llvmfn != llvmfn′
            @assert !haskey(functions(mod), llvmfn′)
            LLVM.name!(f, llvmfn′)
        end
    end

    # rename the entry point
    if job.source.name !== nothing
        LLVM.name!(entry, safe_name(string("julia_", job.source.name)))
    end

    # promote entry-points to kernels and mangle its name
    if job.source.kernel
        entry = promote_kernel!(job, mod, entry)
        LLVM.name!(entry, mangle_call(entry, job.source.tt))
    end

    # minimal required optimization
    @timeit_debug to "rewrite" ModulePassManager() do pm
        global current_job
        current_job = job

        linkage!(entry, LLVM.API.LLVMExternalLinkage)
        internalize!(pm, [LLVM.name(entry)])

        can_throw(job) || add!(pm, ModulePass("LowerThrow", lower_throw!))

        add_lowering_passes!(job, pm)

        run!(pm, mod)

        # NOTE: if an optimization is missing, try scheduling an entirely new optimization
        # to see which passes need to be added to the target-specific list
        #     LLVM.clopts("-print-after-all", "-filter-print-funcs=$(LLVM.name(entry))")
        #     ModulePassManager() do pm
        #         add_library_info!(pm, triple(mod))
        #         add_transform_info!(pm, tm)
        #         PassManagerBuilder() do pmb
        #             populate!(pm, pmb)
        #         end
        #         run!(pm, mod)
        #     end
    end

    return mod, entry
end

## name mangling

# we generate function names that look like C++ functions, because many NVIDIA tools
# support them, e.g., grouping different instantiations of the same kernel together.

function mangle_param(t, substitutions)
    t == Nothing && return "v"

    if isa(t, DataType) || isa(t, Core.Function)
        tn = safe_name(t)

        # handle substitutions
        sub = findfirst(isequal(tn), substitutions)
        if sub === nothing
            str = "$(length(tn))$tn"
            push!(substitutions, tn)
        elseif sub == 1
            str = "S_"
        else
            str = "S$(sub-2)_"
        end

        # encode typevars as template parameters
        if !isempty(t.parameters)
            str *= "I"
            for t in t.parameters
                str *= mangle_param(t, substitutions)
            end
            str *= "E"
        end

        str
    elseif isa(t, Integer)
        "Li$(t)E"
    else
        tn = safe_name(t)
        "$(length(tn))$tn"
    end
end

function mangle_call(f, tt)
    fn = safe_name(f)
    str = "_Z$(length(fn))$fn"

    substitutions = String[]
    for t in tt.parameters
        str *= mangle_param(t, substitutions)
    end

    return str
end

# make names safe for ptxas
safe_name(fn::String) = replace(fn, r"[^A-Za-z0-9_]"=>"_")
safe_name(f::Union{Core.Function,DataType}) = safe_name(String(nameof(f)))
safe_name(f::LLVM.Function) = safe_name(LLVM.name(f))
safe_name(x) = safe_name(repr(x))

## exception handling

# this pass lowers `jl_throw` and friends to GPU-compatible exceptions.
# this isn't strictly necessary, but has a couple of advantages:
# - we can kill off unused exception arguments that otherwise would allocate or invoke
# - we can fake debug information (lacking a stack unwinder)
#
# once we have thorough inference (ie. discarding `@nospecialize` and thus supporting
# exception arguments) and proper debug info to unwind the stack, this pass can go.
function lower_throw!(mod::LLVM.Module)
    job = current_job::CompilerJob
    ctx = context(mod)
    changed = false
    @timeit_debug to "lower throw" begin

    throw_functions = [
        # unsupported runtime functions that are used to throw specific exceptions
        "jl_throw"                      => "exception",
        "jl_error"                      => "error",
        "jl_too_few_args"               => "too few arguments exception",
        "jl_too_many_args"              => "too many arguments exception",
        "jl_type_error"                 => "type error",
        "jl_type_error_rt"              => "type error",
        "jl_undefined_var_error"        => "undefined variable error",
        "jl_bounds_error"               => "bounds error",
        "jl_bounds_error_v"             => "bounds error",
        "jl_bounds_error_int"           => "bounds error",
        "jl_bounds_error_tuple_int"     => "bounds error",
        "jl_bounds_error_unboxed_int"   => "bounds error",
        "jl_bounds_error_ints"          => "bounds error",
        "jl_eof_error"                  => "EOF error",
        # Julia-level exceptions that use unsupported inputs like interpolated strings
        r"julia_throw_exp_domainerror_\d+"      => "DomainError",
        r"julia_throw_complex_domainerror_\d+"  => "DomainError"
    ]

    for f in functions(mod)
        fn = LLVM.name(f)
        for (throw_fn, name) in throw_functions
            occursin(throw_fn, fn) || continue

            for use in uses(f)
                call = user(use)::LLVM.CallInst

                # replace the throw with a PTX-compatible exception
                let builder = Builder(ctx)
                    position!(builder, call)
                    emit_exception!(builder, name, call)
                    dispose(builder)
                end

                # remove the call
                call_args = collect(operands(call))[1:end-1] # last arg is function itself
                unsafe_delete!(LLVM.parent(call), call)

                # HACK: kill the exceptions' unused arguments
                for arg in call_args
                    # peek through casts
                    if isa(arg, LLVM.AddrSpaceCastInst)
                        cast = arg
                        arg = first(operands(cast))
                        isempty(uses(cast)) && unsafe_delete!(LLVM.parent(cast), cast)
                    end

                    if isa(arg, LLVM.Instruction) && isempty(uses(arg))
                        unsafe_delete!(LLVM.parent(arg), arg)
                    end
                end

                changed = true
            end

            @compiler_assert isempty(uses(f)) job
            break
         end
     end

    end
    return changed
end

# report an exception in a GPU-compatible manner
#
# the exact behavior depends on the debug level. in all cases, a `trap` will be emitted, On
# debug level 1, the exception name will be printed, and on debug level 2 the individual
# stack frames (as recovered from the LLVM debug information) will be printed as well.
function emit_exception!(builder, name, inst)
    job = current_job::CompilerJob
    bb = position(builder)
    fun = LLVM.parent(bb)
    mod = LLVM.parent(fun)
    ctx = context(mod)

    # report the exception
    if Base.JLOptions().debug_level >= 1
        name = globalstring_ptr!(builder, name, "exception")
        if Base.JLOptions().debug_level == 1
            call!(builder, Runtime.get(:report_exception), [name])
        else
            call!(builder, Runtime.get(:report_exception_name), [name])
        end
    end

    # report each frame
    if Base.JLOptions().debug_level >= 2
        rt = Runtime.get(:report_exception_frame)
        ft = convert(LLVM.FunctionType, rt, ctx)
        bt = backtrace(inst)
        for (i,frame) in enumerate(bt)
            idx = ConstantInt(parameters(ft)[1], i)
            func = globalstring_ptr!(builder, String(frame.func), "di_func")
            file = globalstring_ptr!(builder, String(frame.file), "di_file")
            line = ConstantInt(parameters(ft)[4], frame.line)
            call!(builder, rt, [idx, func, file, line])
        end
    end

    # signal the exception
    call!(builder, Runtime.get(:signal_exception))

    emit_trap!(job, builder, mod, inst)
end

function emit_trap!(job::CompilerJob, builder, mod, inst)
    ctx = context(mod)
    trap = if haskey(functions(mod), "llvm.trap")
        functions(mod)["llvm.trap"]
    else
        LLVM.Function(mod, "llvm.trap", LLVM.FunctionType(LLVM.VoidType(ctx)))
    end
    call!(builder, trap)
end

## kernel promotion

# promote a function to a kernel
function promote_kernel!(job::CompilerJob, mod::LLVM.Module, kernel::LLVM.Function)
    # pass non-opaque pointer arguments by value (this improves performance,
    # and is mandated by certain back-ends like SPIR-V). only do so for values
    # that aren't a Julia pointer, so we ca still pass those directly.
    kernel_ft = eltype(llvmtype(kernel)::LLVM.PointerType)::LLVM.FunctionType
    kernel_sig = Base.signature_type(job.source.f, job.source.tt)::Type
    kernel_types = filter(dt->!isghosttype(dt) &&
                              (VERSION < v"1.5.0-DEV.581" || !Core.Compiler.isconstType(dt)),
                          [kernel_sig.parameters...])
    @compiler_assert length(kernel_types) == length(parameters(kernel_ft)) job
    for (i, (param_ft,arg_typ)) in enumerate(zip(parameters(kernel_ft), kernel_types))
        if param_ft isa LLVM.PointerType && issized(eltype(param_ft)) &&
           !(arg_typ <: Ptr) && !(VERSION >= v"1.5-" && arg_typ <: Core.LLVMPtr)
            push!(parameter_attributes(kernel, i), EnumAttribute("byval"))
        end
    end

    # target-specific processing
    kernel = process_kernel!(job, mod, kernel)

    return kernel
end


text:
# using MKL
using Test
using DataFlowTasks
using DataFlowTasks.TiledFactorization
using LinearAlgebra
using BenchmarkTools
using VectorizationBase

nc = min(Int(VectorizationBase.num_cores()), Threads.nthreads())

capacity = 50
sch = DataFlowTasks.JuliaScheduler(capacity)
DataFlowTasks.setscheduler!(sch)

tilesize = 256
TiledFactorization.TILESIZE[] = tilesize

PLOT    = true
SAVEFIG = true

nn    = 1000:100:5000 |> collect
t_openblas = Float64[]
t_forkjoin = Float64[]
t_dataflow = Float64[]

for m in nn
    # create an SPD matrix
    A = rand(m,m)
    A = (A + adjoint(A))/2
    A = A + m*I
    nt = Threads.nthreads()

    # test open blas
    @info "OpenBLAS cholesky"
    @info BLAS.get_config()
    BLAS.set_num_threads(nc)
    @info "NBLAS = $(BLAS.get_num_threads())"
    b = @benchmark cholesky!(B) setup=(B=copy($A)) evals=1
    push!(t_openblas,median(b).time)

    # test DataFlowTasks parallelism
    @info "pseudotiled cholesky"
    @info "NBLAS = $(BLAS.get_num_threads())"
    b = @benchmark TiledFactorization.cholesky!(B) setup=(B=copy($A)) evals=1
    push!(t_dataflow,median(b).time)

    # test fork-join parallelism
    b = @benchmark TiledFactorization._cholesky_forkjoin!(B) setup=(B=PseudoTiledMatrix(copy($A),TiledFactorization.TILESIZE[]))
    push!(t_forkjoin,median(b).time)

    # compute the error
    DataFlowTasks.TASKCOUNTER[] = 1 # reset task counter to display how many tasks were created for
    F = TiledFactorization.cholesky(A)

    BLAS.set_num_threads(8)
    er = norm(F.L*F.U-A,Inf)/max(norm(A),norm(F.L*F.U))

    F = cholesky!(copy(A))
    er_blas = norm(F.L*F.U-A,Inf)/max(norm(A),norm(F.L*F.U))

    # for m×m tiled matrix, there should be 1/6*(m^3-m) + 1/2*(m^2-m) + m tasks
    # created
    println("="^80)
    @info "m                 = $(m)"
    @info "Number of tasks   = $(DataFlowTasks.TASKCOUNTER[])"
    @info "Number of threads = $(nt)"
    @info "Number of cores   = $(nc)"
    @info "er                = $er"
    @info "er_blas           = $er_blas"
    @info "t_blas            = $(t_openblas[end])"
    @info "t_dataflow        = $(t_dataflow[end])"
    @info "t_forkjoin        = $(t_forkjoin[end])"
    println("="^80)
end

if PLOT
    using Plots
    flops = @. 1/3*nn^3 + 1/2*nn^2 # I think this is correct (up to O(n)), but double check
    plot(nn,flops./(t_openblas),label="OpenBLAS",xlabel="n",ylabel="GFlops/second",m=:x,title="Cholesky factorization",legend=:bottomright)
    plot!(nn,flops./(t_forkjoin),label="TiledFactorization (forkjoin)",m=:x)
    plot!(nn,flops./(t_dataflow),label="TiledFactorization (dataflow)",m=:x)
    SAVEFIG && savefig(joinpath(DataFlowTasks.PROJECT_ROOT,"benchmarks/choleskyperf_capacity_$(capacity)_tilesize_$(tilesize).png"))
    # peakflops vary, not sure how to measure it. Maybe use cpuinfo?
    # peak = LinearAlgebra.peakflops()
    # plot!(nn,peak/1e9*ones(length(nn)))
    nothing
end


text:
using MicrobiomeAnalysis
using DataFrames, DataStructures
using Plots, MultivariateStats

##### SIMULATION #####

# evaluate numerical solution for experiment 1
t1, Xapp1 = LVmodel(8, [0, 12], h = 3)
# evaluate numerical solution for experiment 2
t2, Xapp2 = LVmodel(12, [0, 10], h = 2)

##### SUMMARIZED EXPERIMENTS #####

# convert transposed time series into Dictionary and store it into assays
assays1 = OrderedDict{String, AbstractArray}("foo" => Xapp1)
assays2 = OrderedDict{String, AbstractArray}("bar" => Xapp2)

# produce feature data including feature name (because it's required by the
# SummarizedExperiment function) and information on genus and species
rowdata1 = DataFrame(
    name = ["strain$i" for i in 1:size(Xapp1, 1)],
    genus = ["g$i" for i in 1:size(Xapp1, 1)],
    species = ["s$i" for i in 1:size(Xapp1, 1)]
)
# produce feature data including feature name (because it's required by the
# SummarizedExperiment function) and information on expression and enzyme
rowdata2 = DataFrame(
    name = ["metabolite$i" for i in 1:size(Xapp2, 1)],
    expression = rand(["intracellular", "extracellular"], size(Xapp2, 1)),
    enzyme = ["enzyme$i" for i in 1:size(Xapp2, 1)]
)

# produce sample data including sample name (because it's required by the
# SummarizedExperiment function) and sampling site
coldata1 = DataFrame(
    name = ["day$i" for i in 1:size(Xapp1, 2)],
    origin = rand(["lake", "ocean", "river"], size(Xapp1, 2)),
    time = t1
)
# produce sample data including sample name (because it's required by the
# SummarizedExperiment function) and wheather conditions
coldata2 = DataFrame(
    name = ["month$i" for i in 1:size(Xapp2, 2)],
    wheather = rand(["sunny", "cloudy", "rainy"], size(Xapp2, 2)),
    time = t2
)

##### MULTIASSAY EXPERIMENT ######

# create a dictionary of multiple experiments
expo = OrderedDict{String, SummarizedExperiment}()
expo["microbiome"] = SummarizedExperiment(assays1, rowdata1, coldata1)
expo["metabolome"] = SummarizedExperiment(assays2, rowdata2, coldata2)

# define general sample data that is not inherently present in experiment columns
sample_data = DataFrame(
     name = ["sample$i" for i in 1:6],
     state = rand(["pure", "polluted", "contaminated", "remediated"], 6)
)

# define mapping between samples and experiment columns
sample_map = DataFrame(
    sample = ["sample1", "sample3", "sample3", "sample3", "sample6", "sample5", "sample2", "sample1", "sample4", "sample4", "sample2"],
    experiment = ["microbiome", "metabolome", "microbiome", "microbiome", "metabolome", "metabolome", "microbiome", "metabolome", "metabolome", "microbiome", "metabolome"],
    colname = ["day1", "month3", "day5", "day3", "month6", "month5", "day2", "month1", "month4", "day4", "month2"]
)

# assemble ingredients into a MultiAssayExperiment object
mae = MultiAssayExperiment(expo, sample_data, sample_map)

##### MANIPULATION #####

# extract microbiome experiment from mae, without sample data
se1 = experiment(mae, "microbiome", sampledata = false)
# extract metabolome experiment from mae, with sample data
se2 = experiment(mae, "metabolome", sampledata = true)

# take a 3-sample subset of mae, keeping all experiments
sample_subset = multifilter(mae, samples = ["sample1", "sample3", "sample3"])
# take a 1-experiment subset of mae, keeping all samples
exp_subset = multifilter(mae, experiments = "metabolome")
# take a 3-sample 1-experiment subset of mae
sample_exp_subset = multifilter(mae, experiments = "metabolome", samples = ["sample1", "sample3", "sample3"])

# extract experiments object
experiments(mae)
# extract sample data object
sampledata(mae)
# extract sample map object
samplemap(mae)


text:
using DataFrames
using CSV
using ArgParse
using MultiChannelFactorAnalysis

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--datafile"
            help = "Path to single_vehicle.csv"
        "--name", "-n"
            help = "Prefix for output files"
            default = "output"
        "--maxmissing"
            help = "radius around focus vehicle to include"
            arg_type = Float64
            default = 0.1
        "--window"
            help = "How many timesteps of history to use"
            default = 50
            arg_type = Int
        "--step"
            help="How many timesteps to step ahead by"
            default=1
            arg_type=Int
        "--tmin"
            help="Which timestep to start at"
            default = -1
            arg_type=Int
        "--tmax"
            help="Which timestep to stop at"
            default = -1
            arg_type=Int
        "--adaptivefac"
            help="Whether to adaptively identify factor structure"
            action= :store_true
        "--nshared"
            help = "Number of common factors"
            arg_type=Int
            default = 1
        "--outfolder"
            help = "Path to write outputs"
            default = "."
        "nspecific"
            nargs = '*'
            arg_type = Int
            help = "Numbers of channel-specific factors"
            default = [1, 1]
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in parsed_args
        println("  $arg  =>  $val")
    end

    println("Reading dataset...")
    dataset = DataFrame(CSV.File(parsed_args["datafile"]))

    tmin = parsed_args["tmin"] == -1 ? minimum(dataset.sample_id) : parsed_args["tmin"]
    tmax = parsed_args["tmin"] == -1 ? maximum(dataset.sample_id) : parsed_args["tmax"]
    windowsize = parsed_args["window"]
    if tmax - tmin < windowsize
        error("$tmax less than $windowsize from $tmin")
    end

    tstart = tmin + windowsize
    params_longs = []
    factors_longs = []
    fitted_longs = []
    for t in tstart:windowsize:tmax
        print("Fitting $t...")
        windowdf = dataset[t - windowsize .< dataset.sample_id .<= t, :]
        windowdata = data_fromlong(windowdf)
        sample_ids = windowdata.sample_ids
        obs_ids = windowdata.obs_ids
        if parsed_args["adaptivefac"]
            error("Adaptive fac not implemented.")
        else
            #nshared = parse(Int, parsed_args[""])
            #nspecific = [parse(Int, p) for p in parsed_args["nspecific"]]
            factorlayout = MultiChannelFactorLayout(parsed_args["nshared"], parsed_args["nspecific"])
        end
        windowfit = fit(windowdata, factorlayout; maxiter=20)
        windowparams = windowfit.finalparams
        windowfactors = extract_factors(windowdata, windowparams)
        windowfitted = predict(windowfactors, windowparams)

        windowparams_long = to_long(windowparams; obs_ids)
        windowfactors_long = to_long(windowfactors; sample_ids)
        windowfitted_long = to_long(MultiChannelFactorAnalysis.unstack(windowfitted, windowdata.channellayout);
            obs_ids=obs_ids, sample_ids=sample_ids)

        windowparams_long[!, :windowtime] .= t
        windowfactors_long[!, :windowtime] .= t
        windowfitted_long[!, :windowtime] .= t
        push!(params_longs, windowparams_long)
        push!(factors_longs, windowfactors_long)
        push!(fitted_longs, windowfitted_long)
    end
    ## combine
    all_params = vcat(params_longs...)
    all_factors = vcat(factors_longs...)
    all_fitted = vcat(fitted_longs...)

    ## output
    CSV.write("$(parsed_args["outfolder"])/$(parsed_args["name"])_params.csv", all_params)
    CSV.write("$(parsed_args["outfolder"])/$(parsed_args["name"])_factors.csv", all_factors)
    CSV.write("$(parsed_args["outfolder"])/$(parsed_args["name"])_fitted.csv", all_fitted)
end

main()


text:
using ProgressMeter

function unbiasedEstimate(model::SMCModel, lM::F1, h::F2,
  ccsmcio::CCSMCIO{Particle}, b::Int64, algorithm::Symbol = :BS,
  independentInitialization::Bool = false, rngCouple::Bool = false,
  maxit::Int64 = typemax(Int64)) where {F1<:Function, F2<:Function, Particle}

  initializeCCSMC(model, lM, ccsmcio, algorithm, independentInitialization)

  ref1 = ccsmcio.ref1
  ref2 = ccsmcio.ref2

  v = h(ref1) # just to get the type of v

  for i in 1:maxit
    ccXpf!(model, lM, ccsmcio, algorithm, rngCouple)
    if i >= b
      if i == b
        v = h(ref1)
      else
        v += h(ref1) - h(ref2)
      end
      if checkEqual(ref1, ref2)
        return (i, v)
      end
    end
  end
  # error("maxit iterations exceeded")
  return (maxit, v)
end

function unbiasedEstimate(model::SMCModel, h::F, ccsmcio::CCSMCIO{Particle},
  b::Int64, independentInitialization::Bool = false, rngCouple::Bool = false,
  maxit::Int64 = typemax(Int64)) where {F<:Function, Particle}
  return unbiasedEstimate(model, error, h, ccsmcio, b, :AT,
    independentInitialization, rngCouple, maxit)
end

function unbiasedEstimates(model::SMCModel, lM::F1, h::F2,
  N::Int64, n::Int64, b::Int64, m::Int64, algorithm::Symbol = :BS,
  independentInitialization::Bool = false, rngCouple::Bool = false,
  maxit::Int64 = typemax(Int64)) where {F1<:Function, F2<:Function, Particle}

  nt = Threads.nthreads()
  ccsmcios::Vector{CCSMCIO{model.particle, model.pScratch}} =
    Vector{CCSMCIO{model.particle, model.pScratch}}(undef, nt)
  Threads.@threads for i in 1:nt
    ccsmcios[i] = CCSMCIO{model.particle, model.pScratch}(N, n)
  end

  # just to get the type of v
  initializeCCSMC(model, ccsmcios[1], independentInitialization)
  v = h(ccsmcios[1].ref1)
  T = typeof(v)

  iterations = Vector{Int64}(undef, m)
  values = Vector{T}(undef, m)

  p = Progress(div(m, nt), 10)
  Threads.@threads for i in 1:m
    tid = Threads.threadid()
    v = unbiasedEstimate(model, lM, h, ccsmcios[tid], b, algorithm,
      independentInitialization, rngCouple, maxit)
    iterations[i] = v[1]
    values[i] = v[2]
    tid == 1 && update!(p, i)
  end

  return iterations, values
end

function unbiasedEstimates(model::SMCModel, h::F, N::Int64, n::Int64,
  b::Int64, m::Int64, independentInitialization::Bool = false,
  rngCouple::Bool = false, maxit::Int64 = typemax(Int64)) where
  {F<:Function, Particle}
  return unbiasedEstimates(model, error, h, N, n, b, m, :AT,
    independentInitialization, rngCouple, maxit)
end