text:
println("\n\nload_s_storm: start")

# number of files to load as sample (0 disables)
# FILE_LIMIT = 2
FILE_LIMIT = 0

using Distributed
if FILE_LIMIT == 0 && nprocs() < length(Sys.cpu_info())
    n = length(Sys.cpu_info()) - nprocs()
    addprocs(n, restrict=true)
end
println("load_s_storm: nprocs=$(nprocs())")

@everywhere prepend!(LOAD_PATH, ["Project.toml"])
@everywhere import CSV, LibPQ
@everywhere using Dates, DataFrames, Glob

@everywhere DISASTER_PASSWORD = read("./etc/disaster-pass", String)

FILE_PATHS = glob("StormEvents_details-*.csv", "./data/noaa_storm")
if FILE_LIMIT > 0
    FILE_PATHS = FILE_PATHS[1:FILE_LIMIT]
end
println("load_s_storm: FILE_PATHS=$(FILE_PATHS)")

postgres = LibPQ.Connection("host=127.0.0.1 port=5432 dbname=disaster user=disaster password=$DISASTER_PASSWORD")

# create empty storm staging table
sql =
"""
drop table if exists public.s_storm;
create table if not exists public.s_storm (
    id                   integer primary key,
    state                text,
    state_fips           text,
    begin_yearmonth      integer,
    begin_day            integer,
    begin_location       text,
    begin_latitude       decimal,
    begin_longitude      decimal,
    end_yearmonth        integer,
    end_day              integer,
    end_location         text,
    end_latitude         decimal,
    end_longitude        decimal,
    event_type           text,
    source               text,
    magnitude            decimal,
    magnitude_type       text,
    tor_f_scale          text,
    category             text,
    ogc_fid              integer,
    datestamp            date,
    severity             integer
);

create or replace function public.trigger_s_storm_udf()
returns trigger as \$body\$
declare
    _ogc_fid  integer;
begin
    select ogc_fid into _ogc_fid
        from public.tl_2015_us_county
        where st_contains(
            tl_2015_us_county.wkb_geometry,
            ST_SetSRID(ST_MakePoint(new.begin_longitude, new.begin_latitude), 4326));

    -- raise notice 'update_s_storm_ogc_fid: id=%, ogc_fid=%', new.id, _ogc_fid;
    new.ogc_fid = _ogc_fid;
    new.datestamp = (left(new.begin_yearmonth::text, 4) || '-' || right(new.begin_yearmonth::text, 2) || '-' || new.begin_day::text)::date;

    -- use hail size, wind speed or tor scale to create 1-10 severity value
    if new.event_type = 'Hail' then
        -- magnitude is hail size (in inches)
        if new.magnitude < 0.25 then new.severity = 1;
        elsif new.magnitude < 0.5 then new.severity = 2;
        elsif new.magnitude < 0.75 then new.severity = 3;
        elsif new.magnitude < 1.0 then new.severity = 4;
        elsif new.magnitude < 1.5 then new.severity = 5;
        elsif new.magnitude < 2.0 then new.severity = 6;
        elsif new.magnitude < 2.5 then new.severity = 7;
        elsif new.magnitude < 3.0 then new.severity = 8;
        elsif new.magnitude < 3.5 then new.severity = 9;
        else new.severity = 10;
        end if;
    elsif new.event_type = 'Tornado' then
        -- tornado uses tor scale
        if new.tor_f_scale = 'F0' then new.severity = 1;
        elsif new.tor_f_scale = 'F1' then new.severity = 2;
        elsif new.tor_f_scale = 'F2' then new.severity = 4;
        elsif new.tor_f_scale = 'F3' then new.severity = 6;
        elsif new.tor_f_scale = 'F4' then new.severity = 8;
        elsif new.tor_f_scale = 'F5' then new.severity = 10;
        else new.severity = 1;
        end if;
    elsif new.event_type = 'Thunderstorm Wind' or new.event_type = 'Winter Storm' then
        -- wind speed
        if new.magnitude < 12 then new.severity = 1;
        elsif new.magnitude < 25 then new.severity = 2;
        elsif new.magnitude < 50 then new.severity = 3;
        elsif new.magnitude < 100 then new.severity = 4;
        elsif new.magnitude < 150 then new.severity = 5;
        elsif new.magnitude < 200 then new.severity = 6;
        elsif new.magnitude < 250 then new.severity = 7;
        elsif new.magnitude < 300 then new.severity = 8;
        elsif new.magnitude < 350 then new.severity = 9;
        else new.severity = 10;
        end if;
    else
        new.severity = 1;
    end if;

    return new;
end; \$body\$ language plpgsql;

drop trigger if exists trigger_s_storm on public.s_storm;
create trigger trigger_s_storm
  before insert on public.s_storm
  for each row execute function public.trigger_s_storm_udf();
"""
LibPQ.execute(postgres, sql)
LibPQ.close(postgres)

# BEGIN_YEARMONTH,BEGIN_DAY,BEGIN_TIME,END_YEARMONTH,END_DAY,END_TIME,EPISODE_ID,EVENT_ID,STATE,STATE_FIPS,YEAR,
# MONTH_NAME,EVENT_TYPE,CZ_TYPE,CZ_FIPS,CZ_NAME,WFO,BEGIN_DATE_TIME,CZ_TIMEZONE,END_DATE_TIME,INJURIES_DIRECT,
# INJURIES_INDIRECT,DEATHS_DIRECT,DEATHS_INDIRECT,DAMAGE_PROPERTY,DAMAGE_CROPS,SOURCE,MAGNITUDE,MAGNITUDE_TYPE,
# FLOOD_CAUSE,CATEGORY,TOR_F_SCALE,TOR_LENGTH,TOR_WIDTH,TOR_OTHER_WFO,TOR_OTHER_CZ_STATE,TOR_OTHER_CZ_FIPS,
# TOR_OTHER_CZ_NAME,BEGIN_RANGE,BEGIN_AZIMUTH,BEGIN_LOCATION,END_RANGE,END_AZIMUTH,END_LOCATION,BEGIN_LAT,
# BEGIN_LON,END_LAT,END_LON,EPISODE_NARRATIVE,EVENT_NARRATIVE,DATA_SOURCE

@everywhere function load_s_storm(filepath)
    println("\tfilepath=$filepath")
    postgres = LibPQ.Connection("host=127.0.0.1 port=5432 dbname=disaster user=disaster password=$DISASTER_PASSWORD")

    df = CSV.read(filepath)
    select!(df, [
        :EVENT_ID, :STATE, :STATE_FIPS,
        :BEGIN_YEARMONTH, :BEGIN_DAY, :BEGIN_LOCATION, :BEGIN_LAT, :BEGIN_LON,
        :END_YEARMONTH, :END_DAY, :END_LOCATION, :END_LAT, :END_LON,
        :EVENT_TYPE, :SOURCE, :MAGNITUDE, :MAGNITUDE_TYPE, :TOR_F_SCALE, :CATEGORY
    ])
    println("\tsize=$(size(df))")
    # println("load_s_storm: columns=$(names(df))")

    LibPQ.load!(
        df,
        postgres,
        """
        insert into public.s_storm (
            id,
            state,
            state_fips,
            begin_yearmonth,
            begin_day,
            begin_location,
            begin_latitude,
            begin_longitude,
            end_yearmonth,
            end_day,
            end_location,
            end_latitude,
            end_longitude,
            event_type,
            source,
            magnitude,
            magnitude_type,
            tor_f_scale,
            category
        ) values (
            \$1,
            \$2,
            \$3,
            \$4,
            \$5,
            \$6,
            \$7,
            \$8,
            \$9,
            \$10,
            \$11,
            \$12,
            \$13,
            \$14,
            \$15,
            \$16,
            \$17,
            \$18,
            \$19
        )
        """
    )

    LibPQ.close(postgres)
end

starttime = Dates.now()
jobs = []
for filepath = FILE_PATHS
    job = @spawn load_s_storm(filepath)
    push!(jobs, job)
end
for job in jobs
    fetch(job)
end

endtime = Dates.now()
println("load_s_storm: elapsed time is $(Dates.canonicalize(Dates.CompoundPeriod(endtime - starttime)))")
println("load_s_storm: done")


text:
x = range(-2, stop=2, length=20)
y1 = @. sin(exp(-x)) + 0.5
y2 = @. sin(exp(-x)) - 0.5
plot(x, y1; label="unsmoothed")
plot!(x, y2; smooth=true, label="smoothed")
legend()

ext: jl


text:
#############################################
# Duffing Oscillator controlled by VSSM #
#############################################

K=500
w=1

#########
# Time  #
#########
δt=1e-3
LONG=Int(2e4)

########################################
# Kinematic Block Parameter (2nd order)#
########################################
Λ=12

####################################
# Parameters for Nominal Trajectory#
####################################

ω=0.5
Ampl=2

##########################
# Exact Model Parameters #
##########################

αₑ=1
δₑ=0.2
βₑ=1
#
################################
# Approximate Model Parameters #
################################
#
αₐ=0.8
δₐ=0.1
βₐ=0.9

###############
# Exact Model #
###############
function Exact(q,q_p,u)
  q_pp=αₑ*q+δₑ*q_p+βₑ*q^3+u
  return q_pp
end

################
# Approx Model #
################

function Approx(q,q_p,q_pp)
  u=q_pp-αₐ*q-δₐ*q_p-βₐ*q^3
  return u
end

#ErrorMetrics(hint,h,h_p) = Λ^2*hint+2*Λ*h+h_p # "S" block for 2nd order

#KinBlock(S,h,h_p,qN_pp) = K*tanh(S/w)+Λ^2*h+2*Λ*h_p+qN_pp #q_ppDes -t adja vissza

ErrorMetrics(hint,h,h_p) = Λ^2*hint+2*Λ*h+h_p## # "S" block for 2nd order

KinBlock(S,h,h_p,qN_pp) = K*tanh(S/w)+Λ^2*h+2*Λ*h_p+qN_pp #q_ppDes -t adja vissza

function nominalTraj(t)
  qN=Ampl*sin(ω*t*δt)
  q_pN=ω*Ampl*cos(ω*t*δt)
  q_ppN=-Ampl*ω^2*sin(ω*t*δt)
  return [qN,q_pN,q_ppN]
end

q_mem=zeros(LONG)
q_p_mem=zeros(LONG)
q_pp_mem=zeros(LONG)
q_ppDes_mem=zeros(LONG)
q_N_mem=zeros(LONG)
q_pN_mem=zeros(LONG)
q_ppN_mem=zeros(LONG)
time_mem=zeros(LONG)
S_mem=zeros(LONG)
u_mem=zeros(LONG)

hint=0
##############
# Simulation #
##############

for t=1:LONG-1
  #create time
  time_mem[t]=δt*t
  #Nominal trajectory for the actual time frame
  q_temp=nominalTraj(t)
  q_N_mem[t] = q_temp[1]
  q_pN_mem[t] = q_temp[2]
  q_ppN_mem[t] = q_temp[3]
  #errors
  h=q_N_mem[t]-q_mem[t]
  #println(h)
  #println(h)
  #Error Integral (Euler integral)
  hint=hint+δt*h
  #The time derivative of the error
  h_p=q_pN_mem[t]-q_p_mem[t]
#  h_p=q_pN_mem[t]-q_p_mem[t]
  #the error metric
  S_mem[t]=ErrorMetrics(hint,h,h_p)
#  S_mem[t]=ErrorMetrics(hint,h,h_p)
  #the kinematic block
  q_ppDes_mem[t]=KinBlock(S_mem[t],h,h_p,q_ppN_mem[t])
  u_mem[t]=Approx(q_mem[t],q_p_mem[t],q_ppDes_mem[t])
  #the exact model
  q_pp_mem[t]=Exact(q_mem[t],q_p_mem[t],u_mem[t])
  #println(q_ppDes_mem[t])
  #println(q_ppDes_mem[t])
  #the aprox mq_ppDes_memExact(q_mem[t],q_p_mem[t],u_mem[t])
  #	Integrate back with Euler integral
  q_p_mem[t+1]=q_p_mem[t]+δt*q_pp_mem[t]
  #println(q_p_mem[t+1])
  q_mem[t+1]=q_mem[t]+q_p_mem[t]*δt
end #for

###############################
# Load matplotlib for Figures #
###############################

using PyPlot
#println(time_mem)
figure(1)
plot(time_mem[3:LONG-1],q_N_mem[3:LONG-1],color="red")
plot(time_mem[3:LONG-1],q_mem[3:LONG-1],color="green")
grid("on")
title("Nominal and Realized Trajectory\n Nominal Trajectory components: bluish\n Realized Trajectory components:orangish ")
xlabel("Time [s]")
ylabel("q [m]")

figure(2)
plot(time_mem[3:LONG-1],q_N_mem[3:LONG-1]-q_mem[3:LONG-1])
grid("on")
title("Tracking Error")
xlabel("Time [s]")
ylabel("Error [m]")

figure(3)
plot(q_N_mem[3:LONG-1],q_pN_mem[3:LONG-1])
plot(q_mem[3:LONG-1],q_p_mem[3:LONG-1])
grid("on")
title("Phase Space")

figure(4)
plot(time_mem[3:LONG-1],u_mem[3:LONG-1])
title("Exerted force")
xlabel("Time [s]")
ylabel("u [N]")

figure(5)
plot(time_mem[3:LONG-1],q_ppN_mem[3:LONG-1],color="red")
plot(time_mem[3:LONG-1],q_ppDes_mem[3:LONG-1],color="green")
plot(time_mem[3:LONG-1],q_pp_mem[3:LONG-1],color="blue")
grid("on")
title("Nominal,Des, Rel.")
xlabel("Time [s]")
ylabel(L"\ddot{q^{N}}")

max_stars_repo_name: krisztiankosi/Control

max_issues_repo_name: krisztiankosi/Control

max_forks_repo_name: krisztiankosi/Control


text:
module TabularTDLearning

using POMDPs
using POMDPLinter
using POMDPPolicies
using POMDPSimulators
using Parameters
using Random

import POMDPs: Solver, solve, Policy

export
    QLearningSolver,
    SARSASolver,
    SARSALambdaSolver,
    solve

include("q_learn.jl")
include("sarsa.jl")
include("sarsa_lambda.jl")

end # module


text:
using Pkg

repl_install_klio() = begin
    Pkg.resolve()

    reducedir = joinpath(dirname(Base.find_package(@__MODULE__, "Reduce")), "..")
    if !isfile(joinpath(reducedir, "deps/reduce.tar.gz"))
        run(`curl -g -L -f -o $reducedir/deps/reduce.tar.gz https://fs.quyo.net/repl.it/klio/reduce-csl_4567_amd64.tgz`)
        run(`chmod u+w $reducedir/deps/build.jl`)
        run(`sed -i -e 's|[^#]download(|#download(|g' $reducedir/deps/build.jl`)
        Pkg.build("Reduce")
    end

    if !isdir(joinpath(homedir(), "opt/maxima"))
        run(`$repldir/install-maxima`)
    end

    pkg"add SQLite#master"
    pkg"precompile"
end

repl_install_klio()


text:
function empty_scalar_deviation(F1 = Float64)
    return ScalarDeviation{F1}(name = :empty)
end

get_weights(d :: ScalarDeviation{F1}) where F1 = one(F1);

function scalar_dev(d :: ScalarDeviation; inclScalarWt :: Bool = true)
    # scalarDev = abs(d.modelV - d.dataV) * d.wtV;
    scalarDev = scalar_deviation(d.modelV, d.dataV, d.wtV; p = norm_p(d));
    if inclScalarWt
        scalarDev *= d.scalarWt;
    end
    scalarStr = sprintf1(d.fmtStr, scalarDev);
    return scalarDev, scalarStr
end

data_se(d :: ScalarDeviation{F1}) where F1 = d.stdV;

## -------------  Show

Base.show(io :: IO, d :: ScalarDeviation{F1}) where F1 =
    Base.print(io, "$(name(d)):  ", short_description(d));

"""
    $(SIGNATURES)

Show a scalar deviation. Fallback if not user-defined function is provided.
Appends to the provided file (if any).
"""
function scalar_show_fct(d :: ScalarDeviation; showModel :: Bool = true, fPath :: String = "")
    io = open_show_path(d, fPath = fPath, writeMode = "a");
    write(io, scalar_show_string(d) * "\n");
    close_show_path(d, io);
    return nothing
end

function scalar_show_string(d :: ScalarDeviation{F1}; showModel :: Bool = true) where F1
    if showModel
        mStr = " model: " * sprintf1(d.fmtStr, d.modelV);
    else
        mStr = "";
    end
    dStr = sprintf1(d.fmtStr, d.dataV);
    if data_se(d) > zero(F1)
        seStr = "(" * sprintf1(d.fmtStr, d.stdV) * ")";
    else
        seStr = "";
    end
    return "$(d.name): $mStr  data: $dStr $seStr"
end

# ---------------


text:
using StatsBase, LinearAlgebra , Distributions # for the density/histogram plot

function fitplot!(all_plots, results, lm, plot_args)

    lhs = terms(lm.updformula.lhs)
    rhs_noint = filter(isnotintercept, terms(lm.updformula.rhs))

    length(rhs_noint) == 0 && return nothing
    length(rhs_noint) > 1 && begin
        println("LinearRegressionKit: Fit plot was requested but not appropriate for regression with multiple independent variables")
        return nothing
    end
    length(lhs) == 0 && return nothing

    actx = first(rhs_noint).sym
    acty = first(lhs).sym
    acttitle = "Fit Plot: " * string(actx) * " vs " * string(acty)

    pred_interval = @vlplot(
            mark = {:errorband, color = "lightgrey"},
            y = { field = :ucli, type = :quantitative, title = acty} ,
            y2 = { field = :lcli, type = :quantitative },
            x = {actx, type = :quantitative})
    if lm.weighted
        pred_interval = @vlplot(
                mark = {:errorbar, ticks=true, color = "dimgrey"},
                y = { field = :ucli, type = :quantitative, title = acty} ,
                y2 = { field = :lcli, type = :quantitative },
                x = {actx, type = :quantitative})
    end

    fp = select(results, [actx, acty, :ucli, :lcli, :uclp, :lclp, :predicted]) |> @vlplot() +
        pred_interval+
        @vlplot(
            mark = {:errorband, color = "darkgrey" },
            y = { field = :uclp, type = :quantitative, title = acty},
            y2 = { field = :lclp, type = :quantitative },
            x = {actx, type = :quantitative} ) +
        @vlplot(
            mark = { :line, color = "darkorange" },
            x = {actx, type = :quantitative},
            y = {:predicted, type = :quantitative}) +
        @vlplot(
            :point,
            x = { actx, type = :quantitative, axis = {grid = false}, scale = {zero = false}},
            y = { acty, type = :quantitative, axis = {grid = false}, scale = {zero = false}},
            title = acttitle, width = plot_args["plot_width"], height = plot_args["plot_width"]
        )

    all_plots["fit"] = fp
end

function simple_residuals_plot(results, dep_var=nothing, show_density=false, st_res=false; plot_width=400, loess_bandwidth::Union{Nothing,Float64}=0.99)
    if isnothing(dep_var)
        dep_var = :predicted
    end

    yaxis = :residuals
    if st_res == true
        yaxis = :student
    end

    loess_p = @vlplot()
    if !isnothing(loess_bandwidth)
        loess_p = @vlplot(
            transform = [ { loess = yaxis, on = dep_var, bandwidth = loess_bandwidth } ],
            mark = {:line, color = "firebrick"},
            x = {dep_var, type = :quantitative},
            y = {yaxis, type = :quantitative} )
    end

    title = "Residuals plot: $(string(dep_var))"
    if st_res == true
        title = "st " * title
    end
    sresults = select(results, [yaxis, dep_var])

    p = sresults |>
    @vlplot(title = title, width = plot_width, height = plot_width,
        x = {type = :quantitative, axis = {grid = false}, scale = {zero = false, padding = 5}},
        y = {type = :quantitative, axis = {grid = false}}) +
    @vlplot(:point,
       x = {dep_var, type = :quantitative},
       y = {yaxis, type = :quantitative}) +
    loess_p +
    @vlplot(mark = {:rule, color = :darkgrey}, y = {type = :quantitative, datum = 0})

    if show_density == false
        return p
    end

    mp = sresults |> @vlplot(
        width = 100, height = plot_width,
        mark = {:area, orient = "horizontal"},
        transform = [{density = yaxis, bandwidth = 0.5}],
        x = {"density:q", title = nothing, axis = nothing},
        y = {"value:q", title = nothing, axis = nothing } )

    tp = @vlplot(bounds = :flush, spacing = 5, config = {view = {stroke = :transparent}}) + [p mp]

    return tp
end

function residuals_plots!(all_plots, results, lm, plot_args)
    rhs_noint = filter(isnotintercept, terms(lm.updformula.rhs))
    plots = Vector{VegaLite.VLSpec}()
    length(rhs_noint) == 0 && return nothing

    plot_width = get(plot_args, "plot_width", 400)
    loess_bw = get(plot_args, "loess_bw", 0.6)
    density_requested = get(plot_args, "residuals_with_density", false)

    # main residual plot
    all_plots["residuals"] = simple_residuals_plot(results, plot_width=plot_width, loess_bandwidth=loess_bw)
    all_plots["st residuals"] = simple_residuals_plot(results,nothing, false, true, plot_width=plot_width, loess_bandwidth=loess_bw)

    # additional plot per dependent variable
    for c_dependent_var in rhs_noint
        c_sym = c_dependent_var.sym
        show_density = density_requested && iscontinuousterm(c_dependent_var)
        all_plots[string("residuals ", string(c_sym))] = simple_residuals_plot(results, c_sym, show_density, plot_width=plot_width, loess_bandwidth=loess_bw)
        all_plots[string("st residuals ", string(c_sym))] = simple_residuals_plot(results, c_sym, show_density, true, plot_width=plot_width, loess_bandwidth=loess_bw)
    end

end

function qqplot(results, fitted_residuals, plot_width)
    n = length(results.residuals)
    grid = [1 / (n - 1):1 / (n - 1):1 - (1 / (n - 1));]
    qu_theo = quantile.(fitted_residuals, grid)
    qu_empi = quantile(results.residuals, grid)
    qqdf = DataFrame(x=qu_theo, y=qu_empi)
    qqline = [first(qu_theo), last(qu_theo)]
    qqldf = DataFrame(x=qqline, y=qqline)

    qqplot = qqdf |> @vlplot() + @vlplot(
        title = "Residuals QQ-Plot",
        width = plot_width, height = plot_width,
        :point,
        x = { :x, type = :quantitative, title = "Theoritical quantiles", axis = {grid = false}, scale = {zero = false, padding = 5} },
        y = { :y, type = :quantitative, title = "Empirical quantiles", axis = {grid = false}, scale = {zero = false, padding = 5} }
        ) + @vlplot(
        {:line, color = "darkgrey"},
        data = qqldf, x = {:x, type = :quantitative}, y = {:y, type = :quantitative} )

    return qqplot
end

function default_range(dist::Distribution, alpha=0.0001)
    minval = isfinite(minimum(dist)) ? minimum(dist) : quantile(dist, alpha)
    maxval = isfinite(maximum(dist)) ? maximum(dist) : quantile(dist, 1 - alpha)
    minval, maxval
end

rice_rule(obs) = round(Int, 2. * obs^(1 / 3))

function histogram_density(results, fitted_residuals, plot_width)
    # data for the density curve
    frmin, frmax = default_range(fitted_residuals)
    rangetest = (frmax - frmin) / plot_width
    qpdf = [pdf(fitted_residuals, x) for x in frmin:rangetest:frmax]
    xs = [x for x in frmin:rangetest:frmax]
    tdf = DataFrame(x=xs, y=qpdf)

    # data for the histogram
    hhh = fit(Histogram, results.residuals)
    nhh = normalize(hhh)
    all_edges = collect(first(nhh.edges))
    bin_starts = [x for x in all_edges[1:end - 1]]
    bin_ends = [x for x in all_edges[2:end]]
    counts = nhh.weights

    hdf = DataFrame(bs=bin_starts, be=bin_ends, y=counts)
    step_size = bin_starts[2] - bin_starts[1]

    hdplot = hdf |> @vlplot(width = plot_width, height = plot_width, title = "Residuals: histogram and PDF") +
    @vlplot(
        :bar,
        x = {:bs, type = :quantitative, title = "residuals", bin = {binned = true, step = step_size}, axis = {grid = false}},
        x2 = {:be, type = :quantitative},
        y = {:y, type = :quantitative, stack = "zero",  axis = {grid = false} } ) +
    @vlplot(
        data = tdf,
        {:line, color = "darkorange"},
        x = {:x, type = :quantitative, scale = {zero = false}, axis = {grid = false}},
        y = {:y, type = :quantitative, scale = {zero = false}, axis = {grid = false}} )

    return hdplot
end

function normality_plots!(all_plots, results, lm, plot_args)
    # plots = Vector{VegaLite.VLSpec}()
    plot_width = get(plot_args, "plot_width", nothing)
    fitted_residuals = fit(Normal, results.residuals)

    all_plots["qq plot"] = qqplot(results, fitted_residuals, plot_width)
    all_plots["histogram density"] = histogram_density(results, fitted_residuals, plot_width)
end

function cooksd_plot!(all_plots, results, lm, plot_args)
    plot_width = get(plot_args, "plot_width", nothing)
    plot_height = plot_width / 2

    sdf = select(results, [:cooksd])
    sdf.Observations = rownumber.(eachrow(sdf))

    threshold_cooksd = 4 / lm.observations
    p = sdf |>
        @vlplot(title = "Cook's Distance", width = plot_width, height = plot_height) +
        @vlplot(mark = {:rule, color = :darkgrey}, y = {type = :quantitative, datum = threshold_cooksd}) +
        @vlplot(
            mark = {:rule, color = :steelblue},
            x = {:Observations, type = :quantitative, axis = {grid = false}},
            y = {type = :quantitative, datum = 0}, y2 = {:cooksd, type = :quantitative} )

    all_plots["cooksd"] = p
    end

function scalelocation_plot!(all_plots, results, lm, plot_args)
    plot_width = get(plot_args, "plot_width", nothing)
    sdf = select(results, [:predicted, :student])
    sdf.sqrtstudent = sqrt.(abs.(sdf.student))

    p = sdf |>
        @vlplot() +
        @vlplot(
                title = "Scale and location plot" ,
                width = plot_width , height = plot_width,
                :point,
                x = {:predicted, type = :quantitative, scale = {zero = false}, axis = { grid = false} },
                y = {:sqrtstudent, type = :quantitative, title = "√student",
                    scale = {zero = false}, axis = { grid = false} }) + @vlplot(
                transform = [ { loess = :sqrtstudent, on = :predicted, bandwidth = 0.6 } ],
                mark = {:line, color = "firebrick"},
                x = {:predicted, type = :quantitative}, y = {:sqrtstudent, type = :quantitative} )
    all_plots["scale location"] = p
end

function leverage_plot!(all_plots, results, lm, plot_args)
    threshold_leverage = 2 * lm.p / lm.observations
    plot_width = plot_args["plot_width"]
    p = select(results, [:leverage, :rstudent]) |>
        @vlplot(title = "Leverage vs Rstudent", width = plot_width, height = plot_width,
            x = {axis = {grid = false}}, y = {axis = {grid = false}}   ) +
        @vlplot(:point,
            x = {:leverage, type = :quantitative},
            y = {:rstudent, type = :quantitative}) +
        @vlplot(mark = {:rule, color = :darkgrey}, y = {datum = -2}) +
        @vlplot(mark = {:rule, color = :darkgrey}, x = {datum = threshold_leverage}) +
        @vlplot(mark = {:rule, color = :darkgrey}, y = {datum = 2})
    all_plots["leverage"] = p
end

function escape_javascript!(vec)
    for i in 1:length(vec)
        vec[i] = replace(vec[i], "&"=>" ampersand ")
    end
end

function ridge_traceplots(rdf::AbstractDataFrame; width=400)
    nb_coefs = round(Int, (ncol(rdf) - 5) / 2)
    dfstats = names(rdf)[2:5]
    dfcoefs = nothing
    dfvifs = nothing

    all_names = names(rdf)
    escape_javascript!(all_names)

    if all_names[6] == "(Intercept)"
        dfcoefs = all_names[7:6 - 1 + nb_coefs]
        push!(dfcoefs, "(Intercept)")
        dfvifs = all_names[end - nb_coefs + 2:end]
    else
        dfcoefs = all_names[6:6 - 1 + nb_coefs]
        dfvifs = all_names[end - nb_coefs + 2:end]
    end

    traceplot = rdf |> @vlplot(title = "Ridge trace plot", width = width, height = width,
        transform = [
            {fold = dfcoefs},
        ],
        mark = :line,
        x = {field = :k , type = :quantitative, scale = {zero = false}, axis = { grid = false}} ,
        y = {field = :value, type = :quantitative, scale = {zero = false}, axis = { grid = false}},
        color = {field = :key, type = :nominal},
        )

    viftraceplot = rdf |> @vlplot(title = "VIF trace plot", width = width, height = width,
        transform = [
            {fold = dfvifs},
        ],
        mark = :line,
        x = {field = :k, type = :quantitative, scale = {zero = false}, axis = { grid = false}} ,
        y = {field = :value, type = :quantitative, scale = {type = :linear, zero = false}, axis = { grid = false} },
        color = {field = :key, type = :nominal},
        )

    viftraceplotlog = rdf |> @vlplot(title = "VIF trace plot", width = width, height = width,
        transform = [
            {fold = dfvifs},
        ],
        mark = :line,
        x = {field = :k, type = :quantitative, scale = {zero = false}, axis = { grid = false}} ,
        y = {field = :value, type = :quantitative, scale = {type = :log, zero = false}, axis = { grid = false} },
        color = {field = :key, type = :nominal},
        )

    rmse_r2_plot = select(rdf, ["k", "RMSE", "R2"]) |> @vlplot(title = "Trace plot", width = width, height = width,
        transform = [
            {fold = ["RMSE", "R2"]},
        ],
        mark = :line,
        x = {field = :k, type = :quantitative, scale = {zero = false}, axis = { grid = false}} ,
        y = {field = :value, type = :quantitative, scale = {zero = false}, axis = { grid = false}},
        color = {field = :key, type = :nominal},
    )

    return Dict([("coefs traceplot", traceplot), ("vifs traceplot", viftraceplot), ("vifs traceplot log", viftraceplotlog), ("rmse r2 plot", rmse_r2_plot)])
end


text:
@testset "Interpolate" begin
    @testset "LineSegment" begin
        for T in Floats
            P₁ = Point{3,T}(1, 1, 1)
            P₂ = Point{3,T}(3, 3, 3)
            l = LineSegment(P₁, P₂)
            @test l(0) ≈ P₁
            @test l(1) ≈ P₂
            @test l(1//2) ≈ Point{3,T}(2, 2, 2)
        end
    end

    @testset "QuadraticSegment" begin
        for T in Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(2, 0, 0)
            P₃ = Point{3,T}(1, 1, 0)
            q = QuadraticSegment(P₁, P₂, P₃)
            for r = LinRange{T}(0, 1, 5)
                @test q(r) ≈ Point{3,T}(2r, -(2r)^2 + 4r, 0)
            end
        end
    end

    @testset "Triangle" begin
        for T in Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(0, 1, 0)
            P₃ = Point{3,T}(0, 0, 1)
            tri = Triangle(P₁, P₂, P₃)
            @test tri(0, 0) ≈ P₁
            @test tri(1, 0) ≈ P₂
            @test tri(0, 1) ≈ P₃
            @test tri(1//2, 1//2) ≈ Point{3,T}(0, 1//2, 1//2)
        end
    end

    @testset "Quadrilateral" begin
        for T in Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(0, 1, 0)
            P₃ = Point{3,T}(0, 1, 1)
            P₄ = Point{3,T}(0, 0, 1)
            quad = Quadrilateral(P₁, P₂, P₃, P₄)
            @test quad(0, 0) ≈ P₁
            @test quad(1, 0) ≈ P₂
            @test quad(1, 1) ≈ P₃
            @test quad(0, 1) ≈ P₄
            @test quad(1//2, 1//2) ≈ Point{3,T}(0, 1//2, 1//2)
        end
    end

    @testset "QuadraticTriangle" begin
        for T in Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(0, 1, 0)
            P₃ = Point{3,T}(0, 1, 1)
            P₄ = Point{3,T}(0, 1//2, 0)
            P₅ = Point{3,T}(0, 1, 1//2)
            P₆ = Point{3,T}(0, 1//2, 1//2)
            tri6 = QuadraticTriangle(P₁, P₂, P₃, P₄, P₅, P₆)
            @test tri6(0, 0) ≈ P₁
            @test tri6(1, 0) ≈ P₂
            @test tri6(0, 1) ≈ P₃
            @test tri6(1//2, 0) ≈ P₄
            @test tri6(1//2, 1//2) ≈ P₅
            @test tri6(0, 1//2) ≈ P₆
        end
    end

    @testset "QuadraticQuadrilateral" begin
        for T in Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(0, 1, 0)
            P₃ = Point{3,T}(0, 1, 1)
            P₄ = Point{3,T}(0, 0, 1)
            P₅ = Point{3,T}(0, 1//2,    0)
            P₆ = Point{3,T}(0,    1, 1//2)
            P₇ = Point{3,T}(0, 1//2,    1)
            P₈ = Point{3,T}(0,    0, 1//2)
            quad8 = QuadraticQuadrilateral(P₁, P₂, P₃, P₄, P₅, P₆, P₇, P₈)
            @test quad8(0, 0) ≈ P₁
            @test quad8(1, 0) ≈ P₂
            @test quad8(1, 1) ≈ P₃
            @test quad8(0, 1) ≈ P₄
            @test quad8(1//2,    0) ≈ P₅
            @test quad8(   1, 1//2) ≈ P₆
            @test quad8(1//2,    1) ≈ P₇
            @test quad8(   0, 1//2) ≈ P₈
            @test quad8(1//2, 1//2) ≈ Point{3,T}(0, 1//2, 1//2)
        end
    end

    @testset "Tetrahedron" begin
        for T in Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(1, 0, 0)
            P₃ = Point{3,T}(0, 1, 0)
            P₄ = Point{3,T}(0, 0, 1)
            tet = Tetrahedron(P₁, P₂, P₃, P₄)
            @test tet(0, 0, 0) ≈ P₁
            @test tet(1, 0, 0) ≈ P₂
            @test tet(0, 1, 0) ≈ P₃
            @test tet(0, 0, 1) ≈ P₄
        end
    end

    @testset "Hexahedron" begin
        for T in Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(1, 0, 0)
            P₃ = Point{3,T}(0, 1, 0)
            P₄ = Point{3,T}(1, 1, 0)
            P₅ = Point{3,T}(0, 0, 1)
            P₆ = Point{3,T}(1, 0, 1)
            P₇ = Point{3,T}(0, 1, 1)
            P₈ = Point{3,T}(1, 1, 1)
            hex = Hexahedron(P₁, P₂, P₃, P₄, P₅, P₆, P₇, P₈)
            @test hex(0, 0, 0) ≈ P₁
            @test hex(1, 0, 0) ≈ P₂
            @test hex(1, 1, 0) ≈ P₃
            @test hex(0, 1, 0) ≈ P₄
            @test hex(0, 0, 1) ≈ P₅
            @test hex(1, 0, 1) ≈ P₆
            @test hex(1, 1, 1) ≈ P₇
            @test hex(0, 1, 1) ≈ P₈
        end
    end

    @testset "QuadraticTetrahedron" begin
        for T in Floats
            P₁  = Point{3,T}(0, 0, 0)
            P₂  = Point{3,T}(1, 0, 0)
            P₃  = Point{3,T}(0, 1, 0)
            P₄  = Point{3,T}(0, 0, 1)
            P₅  = Point{3,T}(1//2,    0,    0)
            P₆  = Point{3,T}(1//2, 1//2,    0)
            P₇  = Point{3,T}(0,    1//2,    0)
            P₈  = Point{3,T}(0,       0, 1//2)
            P₉  = Point{3,T}(1//2,    0, 1//2)
            P₁₀ = Point{3,T}(0,    1//2, 1//2)
            tet10 = QuadraticTetrahedron(P₁, P₂, P₃, P₄, P₅, P₆, P₇, P₈, P₉, P₁₀)
            @test tet10(0, 0, 0) ≈ P₁
            @test tet10(1, 0, 0) ≈ P₂
            @test tet10(0, 1, 0) ≈ P₃
            @test tet10(0, 0, 1) ≈ P₄
            @test tet10(1//2,    0,    0) ≈ P₅
            @test tet10(1//2, 1//2,    0) ≈ P₆
            @test tet10(0,    1//2,    0) ≈ P₇
            @test tet10(0,       0, 1//2) ≈ P₈
            @test tet10(1//2,    0, 1//2) ≈ P₉
            @test tet10(0,    1//2, 1//2) ≈ P₁₀
        end
    end

    @testset "QuadraticHexahedon" begin
        for T in Floats
            P₁  = Point{3,T}(0, 0, 0)
            P₂  = Point{3,T}(1, 0, 0)
            P₃  = Point{3,T}(1, 1, 0)
            P₄  = Point{3,T}(0, 1, 0)
            P₅  = Point{3,T}(0, 0, 1)
            P₆  = Point{3,T}(1, 0, 1)
            P₇  = Point{3,T}(1, 1, 1)
            P₈  = Point{3,T}(0, 1, 1)
            P₉  = Point{3,T}(1//2,    0,    0)
            P₁₀ = Point{3,T}(   1, 1//2,    0)
            P₁₁ = Point{3,T}(1//2,    1,    0)
            P₁₂ = Point{3,T}(   0, 1//2,    0)
            P₁₃ = Point{3,T}(1//2,    0,    1)
            P₁₄ = Point{3,T}(   1, 1//2,    1)
            P₁₅ = Point{3,T}(1//2,    1,    1)
            P₁₆ = Point{3,T}(   0, 1//2,    1)
            P₁₇ = Point{3,T}(0, 0, 1//2)
            P₁₈ = Point{3,T}(1, 0, 1//2)
            P₁₉ = Point{3,T}(1, 1, 1//2)
            P₂₀ = Point{3,T}(0, 0, 1//2)
            hex20 = QuadraticHexahedron(P₁,  P₂,  P₃,  P₄,  P₅,  P₆,  P₇,
                                        P₈,  P₉,  P₁₀, P₁₁, P₁₂, P₁₃, P₁₄,
                                        P₁₅, P₁₆, P₁₇, P₁₈, P₁₉, P₂₀)
            @test hex20(0, 0, 0) ≈ P₁
            @test hex20(1, 0, 0) ≈ P₂
            @test hex20(1, 1, 0) ≈ P₃
            @test hex20(0, 1, 0) ≈ P₄
            @test hex20(0, 0, 1) ≈ P₅
            @test hex20(1, 0, 1) ≈ P₆
            @test hex20(1, 1, 1) ≈ P₇
            @test hex20(0, 1, 1) ≈ P₈
            @test hex20(1//2,    0,    0)≈ P₉
            @test hex20(   1, 1//2,    0)≈ P₁₀
            @test hex20(1//2,    1,    0)≈ P₁₁
            @test hex20(   0, 1//2,    0)≈ P₁₂
            @test hex20(1//2,    0,    1)≈ P₁₃
            @test hex20(   1, 1//2,    1)≈ P₁₄
            @test hex20(1//2,    1,    1)≈ P₁₅
            @test hex20(   0, 1//2,    1)≈ P₁₆
            @test hex20(0, 0, 1//2)≈ P₁₇
            @test hex20(1, 0, 1//2)≈ P₁₈
            @test hex20(1, 1, 1//2)≈ P₁₉
            @test hex20(0, 0, 1//2)≈ P₂₀
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
"""
    tinf = SnoopCompile.flatten_demo()

A simple demonstration of [`@snoopi_deep`](@ref). This demo defines a module

```julia
module FlattenDemo
    struct MyType{T} x::T end
    extract(y::MyType) = y.x
    function packintype(x)
        y = MyType{Int}(x)
        return dostuff(y)
    end
    function domath(x)
        y = x + x
        return y*x + 2*x + 5
    end
    dostuff(y) = domath(extract(y))
end
```

It then "warms up" (forces inference on) all of Julia's `Base` methods needed for `domath`,
to ensure that these MethodInstances do not need to be inferred when we collect the data.
It then returns the results of

```julia
@snoopi_deep FlattenDemo.packintypes(1)
```

See [`flatten`](@ref) for an example usage.
"""
function flatten_demo()
    eval(:(
        module FlattenDemo
        struct MyType{T} x::T end
        extract(y::MyType) = y.x
        function packintype(x)
            y = MyType{Int}(x)
            return dostuff(y)
        end
        function domath(x)
            y = x + x
            return y*x + 2*x + 5
        end
        dostuff(y) = domath(extract(y))
        end
    ))
    z = (1 + 1)*1 + 2*1 + 5
    return @snoopi_deep Base.invokelatest(FlattenDemo.packintype, 1)
end

"""
    tinf = SnoopCompile.itrigs_demo()

A simple demonstration of collecting inference triggers. This demo defines a module

```julia
module ItrigDemo
@noinline double(x) = 2x
@inline calldouble1(c) = double(c[1])
calldouble2(cc) = calldouble1(cc[1])
calleach(ccs) = (calldouble2(ccs[1]), calldouble2(ccs[2]))
end
```

It then "warms up" (forces inference on) `calldouble2(::Vector{Vector{Any}})`, `calldouble1(::Vector{Any})`, `double(::Int)`:

```julia
cc = [Any[1]]
ItrigDemo.calleach([cc,cc])
```

Then it collects and returns inference data using

```julia
cc1, cc2 = [Any[0x01]], [Any[1.0]]
@snoopi_deep ItrigDemo.calleach([cc1, cc2])
```

This does not require any new inference for `calldouble2` or `calldouble1`, but it does force inference on `double` with two new types.
See [`inference_triggers`](@ref) to see what gets collected and returned.
"""
function itrigs_demo()
    eval(:(
        module ItrigDemo
        @noinline double(x) = 2x
        @inline calldouble1(c) = double(c[1])
        calldouble2(cc) = calldouble1(cc[1])
        calleach(ccs) = (calldouble2(ccs[1]), calldouble2(ccs[2]))
        end
    ))
    # Call once to infer `calldouble2(::Vector{Vector{Any}})`, `calldouble1(::Vector{Any})`, `double(::Int)`
    cc = [Any[1]]
    Base.invokelatest(ItrigDemo.calleach, [cc,cc])
    # Now use UInt8 & Float64 elements to force inference on double, without forcing new inference on its callers
    cc1, cc2 = [Any[0x01]], [Any[1.0]]
    return @snoopi_deep Base.invokelatest(ItrigDemo.calleach, [cc1, cc2])
end

"""
    tinf = SnoopCompile.itrigs_higherorder_demo()

A simple demonstration of handling higher-order methods with inference triggers. This demo defines a module

```julia
module ItrigHigherOrderDemo
double(x) = 2x
@noinline function mymap!(f, dst, src)
    for i in eachindex(dst, src)
        dst[i] = f(src[i])
    end
    return dst
end
@noinline mymap(f::F, src) where F = mymap!(f, Vector{Any}(undef, length(src)), src)
callmymap(src) = mymap(double, src)
end
```

The key feature of this set of definitions is that the function `double` gets passed as an argument
through `mymap` and `mymap!` (the latter are [higher-order functions](https://en.wikipedia.org/wiki/Higher-order_function)).

It then "warms up" (forces inference on) `callmymap(::Vector{Any})`, `mymap(::typeof(double), ::Vector{Any})`,
`mymap!(::typeof(double), ::Vector{Any}, ::Vector{Any})` and `double(::Int)`:

```julia
ItrigHigherOrderDemo.callmymap(Any[1, 2])
```

Then it collects and returns inference data using

```julia
@snoopi_deep ItrigHigherOrderDemo.callmymap(Any[1.0, 2.0])
```

which forces inference for `double(::Float64)`.

See [`skiphigherorder`](@ref) for an example using this demo.
"""
function itrigs_higherorder_demo()
    eval(:(
        module ItrigHigherOrderDemo
        double(x) = 2x
        @noinline function mymap!(f, dst, src)
            for i in eachindex(dst, src)
                dst[i] = f(src[i])
            end
            return dst
        end
        @noinline mymap(f::F, src) where F = mymap!(f, Vector{Any}(undef, length(src)), src)
        callmymap(src) = mymap(double, src)
        end
    ))
    # Call once to infer `callmymap(::Vector{Any})`, `mymap(::typeof(double), ::Vector{Any})`,
    #    `mymap!(::typeof(double), ::Vector{Any}, ::Vector{Any})` and `double(::Int)`
    Base.invokelatest(ItrigHigherOrderDemo.callmymap, Any[1, 2])
    src = Any[1.0, 2.0]   # double not yet inferred for Float64
    return @snoopi_deep Base.invokelatest(ItrigHigherOrderDemo.callmymap, src)
end


text:
"""
    finetune!(learner, nepochs[, base_lr = 0.002; kwargs...])

Behaves as [`fastai.Learner.fine_tune`](https://github.com/fastai/fastai/blob/master/fastai/callback/schedule.py#L151)
"""
function finetune!(
        learner,
        nepochs,
        base_lr = 0.002;
        freezeepochs = 1,
        trainlayers = 1:1,
        lr_mult = 100,
        div = 5,
        pct_start = 0.3,
        kwargs...)
    FluxTraining.initlearner!(learner, [TrainingPhase()])
    # freeze backbone and train head
    withfields(learner, model = (FluxTraining.model!, freeze(learner.model, trainlayers))) do
        fitonecycle!(learner, freezeepochs, base_lr, pct_start=0.99; kwargs...)
    end

    base_lr /= 2
    # TODO: use discriminative learning rates
    fitonecycle!(
        learner, nepochs, base_lr;
        div = div, pct_start = pct_start, kwargs...)
    return learner
end


text:
using Base.Test
using Distributions
using Klara

fields = Dict{Symbol, Symbol}(
  :pdf=>:pdf,
  :prior=>:prior,
  :spdf=>:setpdf,
  :sprior=>:setprior,
  :ll=>:loglikelihood!,
  :lp=>:logprior!,
  :lt=>:logtarget!,
  :gll=>:gradloglikelihood!,
  :glp=>:gradlogprior!,
  :glt=>:gradlogtarget!,
  :tll=>:tensorloglikelihood!,
  :tlp=>:tensorlogprior!,
  :tlt=>:tensorlogtarget!,
  :dtll=>:dtensorloglikelihood!,
  :dtlp=>:dtensorlogprior!,
  :dtlt=>:dtensorlogtarget!,
  :uptoglt=>:uptogradlogtarget!,
  :uptotlt=>:uptotensorlogtarget!,
  :uptodtlt=>:uptodtensorlogtarget!
)

println("    Testing BasicContMuvParameter constructors:")

println("      Initialization via index and key fields...")

p = BasicContMuvParameter(:p, 1, signature=:low)

for field in values(fields)
  @test getfield(p, field) == nothing
end

println("      Initialization via pdf field...")

pv = [5.18, -7.76]
μv = [6.11, -8.5]
states = VariableState[BasicContMuvParameterState(pv), BasicMuvVariableState(μv)]

p = BasicContMuvParameter(:p, 1, signature=:low, pdf=MvNormal(states[2].value, 1.), states=states)

distribution = MvNormal(μv, 1.)
p.pdf == distribution
lt, glt = logpdf(distribution, pv), gradlogpdf(distribution, pv)
p.logtarget!(states[1])
@test states[1].logtarget == lt
p.gradlogtarget!(states[1])
@test states[1].gradlogtarget == glt

states[1] = BasicContMuvParameterState(pv)

p.uptogradlogtarget!(states[1])
@test (states[1].logtarget, states[1].gradlogtarget) == (lt, glt)

for field in [:prior, :spdf, :sprior, :ll, :lp, :gll, :glp, :tll, :tlp, :tlt, :dtll, :dtlp, :dtlt, :uptotlt, :uptodtlt]
  @test getfield(p, fields[field]) == nothing
end

pv = [-11.87, -13.44]
states[1].value = pv
μv = [-20.2, -18.91]
states[2].value = μv

p.pdf = MvNormal(states[2].value, 1.)

distribution = MvNormal(μv, 1.)
p.pdf == distribution
lt, glt = logpdf(distribution, pv), gradlogpdf(distribution, pv)
p.logtarget!(states[1])
@test states[1].logtarget == lt
p.gradlogtarget!(states[1])
@test states[1].gradlogtarget == glt

states[1] = BasicContMuvParameterState(pv)

p.uptogradlogtarget!(states[1])
@test (states[1].logtarget, states[1].gradlogtarget) == (lt, glt)

for field in [:prior, :spdf, :sprior, :ll, :lp, :gll, :glp, :tll, :tlp, :tlt, :dtll, :dtlp, :dtlt, :uptotlt, :uptodtlt]
  @test getfield(p, fields[field]) == nothing
end

println("      Initialization via prior field...")

pv = [1.25, 1.8]
pvlen = length(pv)
σv = [10., 2.]
states = VariableState[BasicContMuvParameterState(pv), BasicMuvVariableState(σv)]

p = BasicContMuvParameter(:p, 1, signature=:low, prior=MvNormal(zeros(pvlen), states[2].value), states=states)

distribution = MvNormal(zeros(pvlen), σv)
p.prior == distribution
p.logprior!(states[1])
@test states[1].logprior == logpdf(distribution, pv)
p.gradlogprior!(states[1])
@test states[1].gradlogprior == gradlogpdf(distribution, pv)

for field in [
  :pdf, :spdf,
  :sprior,
  :ll, :lt,
  :gll, :glt,
  :tll, :tlp, :tlt,
  :dtll, :dtlp, :dtlt,
  :uptoglt, :uptotlt, :uptodtlt
]
  @test getfield(p, fields[field]) == nothing
end

pv = [-0.21, 0.98]
pvlen = length(pv)
states[1].value = pv
σv = ones(pvlen)
states[2].value = σv

p.prior = MvNormal(zeros(pvlen), states[2].value)

distribution = MvNormal(zeros(pvlen), σv)
p.prior == distribution
p.logprior!(states[1])
@test states[1].logprior == logpdf(distribution, pv)
p.gradlogprior!(states[1])
@test states[1].gradlogprior == gradlogpdf(distribution, pv)

for field in [
  :pdf, :spdf,
  :sprior,
  :ll, :lt,
  :gll, :glt,
  :tll, :tlp, :tlt,
  :dtll, :dtlp, :dtlt,
  :uptoglt, :uptotlt, :uptodtlt
]
  @test getfield(p, fields[field]) == nothing
end

println("      Initialization via setpdf field...")

pv = [3.79, 4.64]
μv = [5.4, 5.3]
states = VariableState[BasicContMuvParameterState(pv), BasicMuvVariableState(μv)]

p = BasicContMuvParameter(:p, 1, signature=:low, setpdf=(state, states) -> MvNormal(states[2].value, 1.), states=states)
setpdf!(p, states[1])

distribution = MvNormal(μv, 1.)
p.pdf == distribution
lt, glt = logpdf(distribution, pv), gradlogpdf(distribution, pv)
p.logtarget!(states[1])
@test states[1].logtarget == lt
p.gradlogtarget!(states[1])
@test states[1].gradlogtarget == glt

states[1] = BasicContMuvParameterState(pv)

p.uptogradlogtarget!(states[1])
@test (states[1].logtarget, states[1].gradlogtarget) == (lt, glt)

for field in [:prior, :sprior, :ll, :lp, :gll, :glp, :tll, :tlp, :tlt, :dtll, :dtlp, :dtlt, :uptotlt, :uptodtlt]
  @test getfield(p, fields[field]) == nothing
end

pv = [-1.91, -0.9]
states[1].value = pv
μv = [0.12, 0.99]
states[2].value = μv

setpdf!(p, states[1])

distribution = MvNormal(μv, 1.)
p.pdf == distribution
lt, glt = logpdf(distribution, pv), gradlogpdf(distribution, pv)
p.logtarget!(states[1])
@test states[1].logtarget == lt
p.gradlogtarget!(states[1])
@test states[1].gradlogtarget == glt

states[1] = BasicContMuvParameterState(pv)

p.uptogradlogtarget!(states[1])
@test (states[1].logtarget, states[1].gradlogtarget) == (lt, glt)

for field in [:prior, :sprior, :ll, :lp, :gll, :glp, :tll, :tlp, :tlt, :dtll, :dtlp, :dtlt, :uptotlt, :uptodtlt]
  @test getfield(p, fields[field]) == nothing
end

println("      Initialization via setprior field...")

pv = [3.55, 9.5]
pvlen = length(pv)
σv = [2., 10.]
states = VariableState[BasicContMuvParameterState(pv), BasicMuvVariableState(σv)]

p = BasicContMuvParameter(
  :p, 1, signature=:low, setprior=(state, states) -> MvNormal(zeros(pvlen), states[2].value), states=states
)
setprior!(p, states[1])

distribution = MvNormal(zeros(pvlen), σv)
p.prior == distribution
p.logprior!(states[1])
@test states[1].logprior == logpdf(distribution, pv)
p.gradlogprior!(states[1])
@test states[1].gradlogprior == gradlogpdf(distribution, pv)

for field in [:pdf, :spdf, :ll, :lt, :gll, :glt, :tll, :tlp, :tlt, :dtll, :dtlp, :dtlt, :uptoglt, :uptotlt, :uptodtlt]
  @test getfield(p, fields[field]) == nothing
end

pv = [-2.67, 7.71]
states[1].value = pv
σv = [5., 4.]
states[2].value = σv

setprior!(p, states[1])

distribution = MvNormal(zeros(pvlen), σv)
p.prior == distribution
p.logprior!(states[1])
@test states[1].logprior == logpdf(distribution, pv)
p.gradlogprior!(states[1])
@test states[1].gradlogprior == gradlogpdf(distribution, pv)

for field in [:pdf, :spdf, :ll, :lt, :gll, :glt, :tll, :tlp, :tlt, :dtll, :dtlp, :dtlt, :uptoglt, :uptotlt, :uptodtlt]
  @test getfield(p, fields[field]) == nothing
end

# Normal-normal: log-likelihood follows N(μ, σ) and log-prior follows MvNormal(μ0, σ0)
println("      Initialization via loglikelihood! and logprior! fields...")

μv = [-2.637, -1.132]
μvlen = length(μv)
xv = [-1.88, 2.23]
Σv = eye(μvlen)
μ0v = zeros(μvlen)
Σ0v = eye(μvlen)
states = VariableState[
  BasicContMuvParameterState(μv),
  BasicMuvVariableState(xv),
  BasicMavVariableState(Σv),
  BasicMuvVariableState(μ0v),
  BasicMavVariableState(Σ0v)
]

llf(state, states) =
  state.loglikelihood =
  -0.5*(
    (states[2].value-state.value)'*inv(states[3].value)*(states[2].value-state.value)+
    μvlen*log(2*pi)+
    logdet(states[3].value)
  )[1]

lpf(state, states) =
  state.logprior =
  -0.5*(
    (state.value-states[4].value)'*inv(states[5].value)*(state.value-states[4].value)+
    μvlen*log(2*pi)+
    logdet(states[5].value)
  )[1]

μ = BasicContMuvParameter(:μ, 1, signature=:low, loglikelihood=llf, logprior=lpf, states=states)

ld = MvNormal(μv, Σv)
pd = MvNormal(μ0v, Σ0v)
ll, lp = logpdf(ld, xv), logpdf(pd, μv)
lt = ll+lp
μ.loglikelihood!(states[1])
@test_approx_eq states[1].loglikelihood ll
μ.logprior!(states[1])
@test_approx_eq states[1].logprior lp

states[1] = BasicContMuvParameterState(μv)

μ.logtarget!(states[1])
@test_approx_eq states[1].loglikelihood ll
@test_approx_eq states[1].logprior lp
@test_approx_eq states[1].logtarget lt

for field in [
  :pdf, :prior,
  :spdf, :sprior,
  :gll, :glp, :glt,
  :tll, :tlp, :tlt,
  :dtll, :dtlp, :dtlt,
  :uptoglt, :uptotlt, :uptodtlt
]
  @test getfield(μ, fields[field]) == nothing
end

# Unnormalized normal target
println("      Initialization via logtarget! field...")

pv = [-1.28, 1.73]
pvlen = length(pv)
μv = [9.4, 3.32]
states = VariableState[BasicContMuvParameterState(pv), BasicMuvVariableState(μv)]

p = BasicContMuvParameter(
  :p,
  1,
  signature=:low,
  logtarget=(state, states) -> state.logtarget = -(state.value-states[2].value)⋅(states[1].value-states[2].value),
  states=states
)

p.logtarget!(states[1])
@test_approx_eq 0.5*(states[1].logtarget-pvlen*log(2*pi))[1] logpdf(MvNormal(μv, 1.), pv)

for field in [
  :pdf, :prior,
  :spdf, :sprior,
  :ll, :lp,
  :gll, :glp, :glt,
  :tll, :tlp, :tlt,
  :dtll, :dtlp, :dtlt,
  :uptoglt, :uptotlt, :uptodtlt
]
  @test getfield(p, fields[field]) == nothing
end

# Normal-normal: log-likelihood follows N(μ, σ) and log-prior follows MvNormal(μ0, σ0)
println("      Initialization via loglikelihood!, gradloglikelihood! and prior fields...")

μv = [5.59, -7.25]
μvlen = length(μv)
xv = [4.11, 8.17]
Σv = eye(μvlen)
μ0v = zeros(μvlen)
Σ0v = eye(μvlen)
states = VariableState[
  BasicContMuvParameterState(μv),
  BasicMuvVariableState(xv),
  BasicMavVariableState(Σv),
  BasicMuvVariableState(μ0v),
  BasicMavVariableState(Σ0v)
]

llf(state, states) =
  state.loglikelihood =
  -0.5*(
    (states[2].value-states[1].value)'*inv(states[3].value)*(states[2].value-states[1].value)+
    μvlen*log(2*pi)+
    logdet(states[3].value)
  )[1]

gllf(state, states) = states[1].gradloglikelihood = (states[3].value\(states[2].value-states[1].value))

μ = BasicContMuvParameter(
  :μ,
  1,
  signature=:low,
  loglikelihood=llf,
  gradloglikelihood=gllf,
  prior=MvNormal(states[4].value, states[5].value),
  states=states
)

ld = MvNormal(μv, Σv)
pd = MvNormal(μ0v, Σ0v)
ll, lp = logpdf(ld, xv), logpdf(pd, μv)
lt = ll+lp
gll, glp = -gradlogpdf(ld, xv), gradlogpdf(pd, μv)
glt = gll+glp
μ.loglikelihood!(states[1])
@test_approx_eq states[1].loglikelihood ll
μ.logprior!(states[1])
@test states[1].logprior == lp
μ.gradloglikelihood!(states[1])
@test_approx_eq states[1].gradloglikelihood gll
μ.gradlogprior!(states[1])
@test states[1].gradlogprior == glp

states[1] = BasicContMuvParameterState(μv)

μ.logtarget!(states[1])
@test_approx_eq states[1].loglikelihood ll
@test states[1].logprior == lp
@test_approx_eq states[1].logtarget lt
μ.gradlogtarget!(states[1])
@test_approx_eq states[1].gradloglikelihood gll
@test states[1].gradlogprior == glp
@test_approx_eq states[1].gradlogtarget glt

states[1] = BasicContMuvParameterState(μv)

μ.uptogradlogtarget!(states[1])
@test_approx_eq states[1].loglikelihood ll
@test states[1].logprior == lp
@test_approx_eq states[1].logtarget lt
@test_approx_eq states[1].gradloglikelihood gll
@test states[1].gradlogprior == glp
@test_approx_eq states[1].gradlogtarget glt

for field in [:pdf, :spdf, :sprior, :tll, :tlp, :tlt, :dtll, :dtlp, :dtlt, :uptotlt, :uptodtlt]
  @test getfield(μ, fields[field]) == nothing
end

μv = [4.21, 7.91]
states[1].value = μv
xv = [-3.1, -2.52]
states[2].value = xv
Σv = diagm([2., 1.])
states[3].value = Σv
μ0v = [1., 2.5]
states[4].value = μ0v
Σ0v = diagm([3., 5.])
states[5].value = Σ0v

μ.prior = MvNormal(states[4].value, states[5].value)

ld = MvNormal(μv, Σv)
pd = MvNormal(μ0v, Σ0v)
ll, lp = logpdf(ld, xv), logpdf(pd, μv)
lt = ll+lp
gll, glp = -gradlogpdf(ld, xv), gradlogpdf(pd, μv)
glt = gll+glp
μ.loglikelihood!(states[1])
@test_approx_eq states[1].loglikelihood ll
μ.logprior!(states[1])
@test states[1].logprior == lp
μ.gradloglikelihood!(states[1])
@test_approx_eq states[1].gradloglikelihood gll
μ.gradlogprior!(states[1])
@test states[1].gradlogprior == glp

states[1] = BasicContMuvParameterState(μv)

μ.logtarget!(states[1])
@test_approx_eq states[1].loglikelihood ll
@test states[1].logprior == lp
@test_approx_eq states[1].logtarget lt
μ.gradlogtarget!(states[1])
@test_approx_eq states[1].gradloglikelihood gll
@test states[1].gradlogprior == glp
@test_approx_eq states[1].gradlogtarget glt

states[1] = BasicContMuvParameterState(μv)

μ.uptogradlogtarget!(states[1])
@test_approx_eq states[1].loglikelihood ll
@test states[1].logprior == lp
@test_approx_eq states[1].logtarget lt
@test_approx_eq states[1].gradloglikelihood gll
@test states[1].gradlogprior == glp
@test_approx_eq states[1].gradlogtarget glt

for field in [:pdf, :spdf, :sprior, :tll, :tlp, :tlt, :dtll, :dtlp, :dtlt, :uptotlt, :uptodtlt]
  @test getfield(μ, fields[field]) == nothing
end

# Normal-normal: log-likelihood follows N(μ, σ) and log-prior follows MvNormal(μ0, σ0)
println("      Initialization via loglikelihood!, logprior!, gradloglikelihood! and gradlogprior! fields...")

μv = [6.69, -3.125]
μvlen = length(μv)
xv = [5.43, 9.783]
Σv = eye(μvlen)
μ0v = zeros(μvlen)
Σ0v = eye(μvlen)
states = VariableState[
  BasicContMuvParameterState(μv),
  BasicMuvVariableState(xv),
  BasicMavVariableState(Σv),
  BasicMuvVariableState(μ0v),
  BasicMavVariableState(Σ0v)
]

llf(state, states) =
  state.loglikelihood =
  -0.5*(
    (states[2].value-state.value)'*inv(states[3].value)*(states[2].value-state.value)+
    μvlen*log(2*pi)+
    logdet(states[3].value)
  )[1]

lpf(state, states) =
  state.logprior =
  -0.5*(
    (state.value-states[4].value)'*inv(states[5].value)*(state.value-states[4].value)+
    μvlen*log(2*pi)+
    logdet(states[5].value)
  )[1]

gllf(state, states) = state.gradloglikelihood = (states[3].value\(states[2].value-state.value))

glpf(state, states) = state.gradlogprior = -(states[5].value\(state.value-states[4].value))

μ = BasicContMuvParameter(
  :μ, 1, signature=:low, loglikelihood=llf, logprior=lpf, gradloglikelihood=gllf, gradlogprior=glpf, states=states
)

ld = MvNormal(μv, Σv)
pd = MvNormal(μ0v, Σ0v)
ll, lp = logpdf(ld, xv), logpdf(pd, μv)
lt = ll+lp
gll, glp = -gradlogpdf(ld, xv), gradlogpdf(pd, μv)
glt = gll+glp
μ.loglikelihood!(states[1])
@test_approx_eq states[1].loglikelihood ll
μ.logprior!(states[1])
@test_approx_eq states[1].logprior lp
μ.gradloglikelihood!(states[1])
@test_approx_eq states[1].gradloglikelihood gll
μ.gradlogprior!(states[1])
@test_approx_eq states[1].gradlogprior glp

states[1] = BasicContMuvParameterState(μv)

μ.logtarget!(states[1])
@test_approx_eq states[1].loglikelihood ll
@test_approx_eq states[1].logprior lp
@test_approx_eq states[1].logtarget lt
μ.gradlogtarget!(states[1])
@test_approx_eq states[1].gradloglikelihood gll
@test_approx_eq states[1].gradlogprior glp
@test_approx_eq states[1].gradlogtarget glt

states[1] = BasicContMuvParameterState(μv)

μ.uptogradlogtarget!(states[1])
@test_approx_eq states[1].loglikelihood ll
@test_approx_eq states[1].logprior lp
@test_approx_eq states[1].logtarget lt
@test_approx_eq states[1].gradloglikelihood gll
@test_approx_eq states[1].gradlogprior glp
@test_approx_eq states[1].gradlogtarget glt

for field in [:pdf, :prior, :spdf, :sprior, :tll, :tlp, :tlt, :dtll, :dtlp, :dtlt, :uptotlt, :uptodtlt]
  @test getfield(μ, fields[field]) == nothing
end

# Unnormalized normal target
println("      Initialization via logtarget! and gradlogtarget! fields...")

pv = [-4.29, 2.91]
μv = [2.2, 2.02]
states = VariableState[BasicContMuvParameterState(pv), BasicMuvVariableState(μv)]

p = BasicContMuvParameter(
  :p,
  1,
  signature=:low,
  logtarget=(state, states) -> state.logtarget = -(state.value-states[2].value)⋅(states[1].value-states[2].value),
  gradlogtarget=(state, states) -> state.gradlogtarget = -2*(state.value-states[2].value),
  states=states
)

distribution = MvNormal(μv, 1.)
lt, glt = logpdf(distribution, pv), gradlogpdf(distribution, pv)
p.logtarget!(states[1])
@test_approx_eq 0.5*(states[1].logtarget-pvlen*log(2*pi))[1] lt
p.gradlogtarget!(states[1])
@test_approx_eq 0.5*states[1].gradlogtarget glt

states[1] = BasicContMuvParameterState(pv)

p.uptogradlogtarget!(states[1])
@test_approx_eq 0.5*(states[1].logtarget-pvlen*log(2*pi))[1] lt
@test_approx_eq 0.5*states[1].gradlogtarget glt

for field in [
  :pdf, :prior,
  :spdf, :sprior,
  :ll, :lp,
  :gll, :glp,
  :tll, :tlp, :tlt,
  :dtll, :dtlp, :dtlt,
  :uptotlt, :uptodtlt
]
  @test getfield(p, fields[field]) == nothing
end


text:
struct SBMEstimator
    exact::Bool                 # Whether the method is exact, i.e. returns a global optimum
    method::String              # Name of the method
    time_limit::Float64         # Time limit in seconds
    verbose::Bool               # Verbose
    seed::Union{Nothing,Int}    # Random seed (for assignments initialization)
    accept_early::Bool          # Specific option for local search 1

    function SBMEstimator(method::String, time_limit::Float64, verbose::Bool, seed::Union{Nothing,Int}, accept_early::Bool)
        if ~(method in ["ls1","ls2","ls3","exact"])
            throw(ArgumentError("Invalid method: $method."))
        end
        if time_limit <= 0
            throw(DomainError("Time limit must be a positive value (in seconds)."))
        end
        if method == "exact"
            exact = true
        else
            exact = false
        end
        return new(exact, method, time_limit, verbose, seed, accept_early)
    end

    function SBMEstimator(method::String, time_limit::Float64, verbose::Bool, accept_early::Bool)
        seed = nothing
        return SBMEstimator(method, time_limit, verbose, seed, accept_early)
    end
end

include("local_search1.jl")
include("local_search2.jl")
include("local_search3.jl")

function run(estimator::SBMEstimator, dataset::Dataset)::Tuple{SBM, Matrix{Int}, OptResults}
    if estimator.method == "ls1"
        return localSearch1(estimator, dataset)
    elseif estimator.method == "ls2"
        return localSearch2(estimator, dataset)
    elseif estimator.method == "ls3"
        return localSearch3(estimator, dataset)
    elseif estimator.method == "exact"
        return MINLP(estimator, dataset)
    end
end

function calculateObjective(dataset::Dataset, w::Matrix{Float64}, x::Matrix{Int})::Float64
    """ Calculates the maximum log-likelihood value

    Parameters
    ----------
    dataset : Dataset
        Dataset representing an observed graph.
    sbm : SBM
        An Stochastic Block Model
    x : Matrix{Int}
        Matrix of assignments of nodes to groups

    Returns
    -------
    L : Float64
        Maximum log-likelihood value
    """
    A = dataset.A
    n = dataset.n
    m = dataset.m
    q = dataset.n_communities
    k = dataset.k
    if ~((n,q) == size(x))
        throw(ArgumentError("The dimensions of the assignment matrix do not match the dimensions of the dataset."))
    end
    # Calculate objective value
    L = 0
    for g=1:q, h=1:q, i=1:n, j=1:n
        if x[i,g] * x[j,h] == 1
            if A[i,j] != 0
                L += 0.5*(- A[i,j]*log(w[g,h]) + ((k[i]*k[j])/(2*m)) * w[g,h])
            else
                L += 0.5*( ((k[i]*k[j])/(2*m)) * w[g,h] )
            end
        end
    end
    return L
end

function optimalProbMatrix(dataset::Dataset, x::Matrix{Int})::Matrix{Float64}
    """ Calculates the optimal probability matrix given the assignments

    Parameters
    ----------
    dataset : Dataset
        Dataset representing an observed graph.
    x : Array{Int64,2}
        Matrix of assignments of nodes to groups

    Returns
    -------
    w : Array{Int64,2}
        Matrix of probabilities
    """
    A = dataset.A
    n = dataset.n
    m = dataset.m
    q = dataset.n_communities
    k = dataset.k
    if ~((n,q) == size(x))
        throw(ArgumentError("The dimensions of the assignment matrix do not match the dimensions of the dataset."))
    end

    w = zeros(q,q)
    for g=1:q, h=1:q
        numerator = 0.5*sum(A[i,j]*x[i,g]*x[j,h] for i=1:n, j=1:n)
        denominator = 0.5*sum( ((k[i]*k[j])/(2*m)) *x[i,g]*x[j,h] for i=1:n, j=1:n)
        if numerator == 0
            w[g,h] = 0
        else
            w[g,h] = numerator / denominator
        end
    end
    return w
end

function estimate(dataset::Dataset; method::String="ls1", time_limit::Float64=400.0, verbose::Bool=false)::Tuple{SBM, Matrix{Int}, OptResults}
    estimator = SBMEstimator(method, time_limit, verbose, true)
    sbm, x, opt_results = run(estimator, dataset)
    return sbm, x, opt_results
end

function estimate(A::Matrix{Int}, q::Int; method::String="ls1", time_limit::Float64=400.0, verbose::Bool=false)::Tuple{SBM, Matrix{Int}, OptResults}
    dataset = Dataset(A, q)
    estimator = SBMEstimator(method, time_limit, verbose, true)
    sbm, x, opt_results = run(estimator, dataset)
    return sbm, x, opt_results
end


text:
module DiffEqsSolver

export EulerDES, RK2DES, RK4DES, EulerCromerDES, EulerMidPiontDES, VerletDES, VelVerDES

function EulerDES(; ẋ::Function, x₀::Vector{T}, t₀::T, t₁::T, h::T) where {T<:AbstractFloat}
    tcoll = collect(t₀:h:t₁)
    xcoll = [x₀]
    x = x₀
    for t ∈ tcoll[2:end]
        x += h * ẋ(t, x)
        push!(xcoll, x)
    end
    return tcoll, xcoll
end

function UnstableDES(; ẋ::Function, x₀::Vector{T}, t₀::T, t₁::T, h::T) where {T<:AbstractFloat}
    tcoll = collect(t₀:h:t₁)
    xcoll = [x₀, x₀ + h * ẋ(t₀, x₀)]
    for t ∈ tcoll[3:end]
        xₙ₊₁ = 2 * h * ẋ(t, xcoll[end]) + xcoll[end-1]
        push!(xcoll, xₙ₊₁)
    end
    return tcoll, xcoll
end

function EulerYaghoubDES(; ẋ::Function, x₀::Vector{T}, t₀::T, t₁::T, h::T) where {T<:AbstractFloat}
    tcoll = collect(t₀:h:t₁)
    xcoll = [x₀]
    x = x₀
    ẋval = ẋ(t₀, x₀)
    for t ∈ tcoll[2:end]
        x += h * ẋval
        ẋval = ẋ(t, x)
        push!(xcoll, x)
    end
    return tcoll, xcoll
end

function RK2DES(; ẋ::Function, x₀::Vector{T}, t₀::T, t₁::T, h::T) where {T<:AbstractFloat}
    tcoll = collect(t₀:h:t₁)
    xcoll = [x₀]
    x = x₀
    for t ∈ tcoll[2:end]
        K₁ = ẋ(t, x) * h
        K₂ = h * ẋ(t + h / 2, x + (K₁ * h) / 2)
        x += K₂
        push!(xcoll, x)
    end
    return tcoll, xcoll
end

function RK4DES(; ẋ::Function, x₀::Vector{T}, t₀::T, t₁::T, h::T) where {T<:AbstractFloat}
    tcoll = collect(t₀:h:t₁)
    xcoll = [x₀]
    x = x₀
    for t ∈ tcoll[2:end]
        K₁ = ẋ(t, x)
        K₂ = ẋ(t + h / 2, x .+ h / 2 * K₁)
        K₃ = ẋ(t + h / 2, x .+ h / 2 * K₂)
        K₄ = ẋ(t + h, x .+ h * K₃)
        x += (K₁ + 2 * K₂ + 2 * K₃ + K₄) * h / 6
        push!(xcoll, x)
    end
    return tcoll, xcoll
end

function EulerCromerDES(; ẍ::Function, ẋ₀::Vector{T}, x₀::Vector{T}, t₀::T, t₁::T, h::T) where {T<:AbstractFloat}
    tcoll = collect(t₀:h:t₁)
    xcoll = [x₀]
    ẋcoll = [ẋ₀]
    x = x₀
    ẋ = ẋ₀
    for t ∈ tcoll[2:end]
        ẋ += h * ẍ(t, x, ẋ)
        x += h * ẋ
        push!(xcoll, x)
        push!(ẋcoll, ẋ)
    end

    return tcoll, xcoll, ẋcoll
end

function EulerMidPiontDES(; ẍ::Function, ẋ₀::Vector{T}, x₀::Vector{T}, t₀::T, t₁::T, h::T) where {T<:AbstractFloat}
    tcoll = collect(t₀:h:t₁)
    xcoll = [x₀]
    ẋcoll = [ẋ₀ - (ẍ(t₀, x₀, ẋ₀) * h / 2)]
    x = x₀
    ẋ = ẋ₀ - (ẍ(t₀, x₀, ẋ₀) * h / 2)

    for t ∈ tcoll[2:end]
        ẋ += h * ẍ(t, x, ẋ)
        x += h * ẋ

        push!(xcoll, x)
        push!(ẋcoll, ẋ)
    end

    return tcoll, xcoll, ẋcoll
end

function VerletDES(; ẍ::Function, ẋ₀::Vector{T}, x₀::Vector{T}, t₀::T, t₁::T, h::T) where {T<:AbstractFloat}
    tcoll = collect(t₀:h:t₁)
    xcoll = [x₀, x₀ + ẋ₀ * h + ẍ(t₀, x₀, ẋ₀) * (h^2) / 2]
    ẋcoll = [ẋ₀, ẋ₀ + ẍ(t₀, x₀, ẋ₀) * h]

    for i ∈ 3:length(tcoll)
        xᵢ = 2 * xcoll[i-1] - xcoll[i-2] + ẍ(tcoll[i-1], xcoll[i-1], ẋcoll[i-1]) * (h^2)
        push!(xcoll, xᵢ)
        ẋᵢ = (xcoll[i] - xcoll[i-1]) / h
        push!(ẋcoll, ẋᵢ)
    end

    return tcoll, xcoll, ẋcoll
end

function VelVerDES(; ẍ::Function, ẋ₀::Vector{T}, x₀::Vector{T}, t₀::T, t₁::T, h::T) where {T<:AbstractFloat}
    tcoll = collect(t₀:h:t₁)
    xcoll = [x₀]
    ẋcoll = [ẋ₀]
    x = x₀
    ẋ = ẋ₀

    for t ∈ tcoll[2:end]
        ẍₙ = ẍ(t, x, ẋ)
        x += ẋ * h + ẍₙ * (h^2) / 2
        ẍₙ₊₁ = ẍ(t + h, x, ẋ)
        ẋ += (ẍₙ + ẍₙ₊₁) * h / 2
        push!(xcoll, x)
        push!(ẋcoll, ẋ)
    end

    return tcoll, xcoll, ẋcoll
end

end


text:
using WavePropBase
using Documenter

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

makedocs(;
    modules=[WavePropBase],
    authors="Luiz M. Faria <maltezfaria@gmail.com> and contributors",
    repo="https://github.com/WaveProp/WavePropBase.jl/blob/{commit}{path}#{line}",
    sitename="WavePropBase.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://WaveProp.github.io/WavePropBase.jl",
        assets=String[],
    ),
    pages=[
        "WavePropBase" => "index.md",
        "Utils module" => "utils.md",
        "Geometry module" => "geometry.md",
        "Interpolation module" => "interpolation.md",
        "Integration module" => "integration.md",
        "Mesh module" => "mesh.md",
        "Trees module" => "trees.md",
        "IO module" => "io.md",
        "References" => "references.md"
    ],
)

deploydocs(;
    repo="github.com/WaveProp/WavePropBase.jl",
    devbranch="main",
)

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
n_values_1 = [200]
n_values_2 = [200]
m_values_1 = [200]
m_values_2 = [200]
sz_values  = [210]

testset_name_suffix = "(coverage)"

include("_public_mul.jl")


text:
using InteractBase, Observables, DataStructures, Colors, WebIO, CSSUtil
using Widgets
import InteractBase: widgettype
import Widgets: components

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

include("test_observables.jl")


text:
function EuropeanCallMC_Antithetic(K, T, S, sig, r, div, N, M)

## setting parameters

dt = T/N
nudt = (r - div - 0.5 * sig^2) * dt
sigsdt = sig * sqrt(dt)
lnS = log(S)

## variables to store the sums of the option values

sum_CT = 0
sum_CT2 = 0

## outer loop is the number of simulations bigger is better

for j = 1:M

    ## setting up two path variables
  lnSt1 = lnS
  lnSt2 = lnS

  for i = 1:N

      # antithetic variance reduction through the perectly negative correlated error term
    error = randn()
    lnSt1 = lnSt1 + nudt + sigsdt * (error)
    lnSt2 = lnSt2 + nudt + sigsdt * (-error)

  end

    ## wraps up and stores all the values of each path

  ST1 = exp(lnSt1)
  ST2 = exp(lnSt2)
  CT = 0.5 * ( max(0, ST1 - K) + max(0, ST2 - K))
  sum_CT = sum_CT + CT
  sum_CT2 = sum_CT2 + CT^2

end

  ## discounting the average call value
call_value = sum_CT / M * exp(-r*T)

SD = (sqrt(sum_CT2 - 1/M * (sum_CT^2)) * exp(-2*r*T)) / (M-1)
SE = SD / sqrt(M)

println("The Call Value is $call_value")
println("The SD is $SD")
println("The SE is $SE")

MC = [call_value, SD, SE]

end


text:
"""
Numeric representation of SAT problem in DIMACS-CNF form

Input:
- `num_variables` is the number of symbolic Boolean variables
- `clauses` are integer vectors; negative numbers indicate negated variables
"""

"""SAT problem represented numerically in DIMACS-CNF form:
clauses are integer vectors with negative numbers indicating negated literals
"""
mutable struct SATProblem
    num_variables::Int
    clauses::Vector{Vector{Int}}
end

max_var(clause) = maximum(abs.(clause))

"Constructor which calculates the number of variables automatically"
function SATProblem(clauses::Vector{Vector{Int}})
    num_variables = maximum(max_var.(clauses))

    return SATProblem(num_variables, clauses)
end

Base.push!(p::SATProblem, clause) = push!(p.clauses, clause)

"Ban a clause, i.e. add the negative of the clause"
ban!(p::SATProblem, clause) = push!(p, .-(clause))

function all_solutions(p::SATProblem)

    solns = []
    status = :sat

    while true
        status, results = solve(p)
        status == :unsat && break

        push!(solns, results)
        ban!(p, results)

    end

    return solns
end


text:
using Knet,AutoGrad
a = map(t->map(a->isa(a,Array)?KnetArray(a):a,t), AutoGrad.alltests());
for ai in a; println(ai); AutoGrad.runtests([ai]); end


text:
using Pkg

pkg"""
add Documenter Literate Plots Yao
"""

using Documenter, Literate, QuAlgorithmZoo

const PATH = (
    tutorial = joinpath(@__DIR__, "src/tutorial"),
    examples = joinpath(@__DIR__, "..", "examples")
)

function process_literate_scripts(;excludes=["make.jl", "README.md"])
    TUTORIALS = []
    for (root, dirs, files) in walkdir(PATH.examples)
        for file in files
            file in excludes && continue
            filepath = joinpath(root, file)
            Literate.markdown(filepath, PATH.tutorial)

            filename, _ = splitext(file)
            mdfile = join([filename, ".md"])
            # TODO: use PATH.tutorial rather then manual path
            push!(TUTORIALS, relpath(joinpath("tutorial", mdfile)))
        end
    end
    TUTORIALS
end

const TUTORIALS = process_literate_scripts()

#-----------------------------------------------

makedocs(
    modules = [QuAlgorithmZoo],
    clean = false,
    format = :html,
    sitename = "Quantum Algorithm Zoo",
    linkcheck = !("skiplinks" in ARGS),
    analytics = "UA-89508993-1",
    pages = [
        "Home" => "index.md",
        "Tutorial" => TUTORIALS,
        "Manual" => Any[
            "man/zoo.md",
        ],
    ],
    html_prettyurls = !("local" in ARGS),
    html_canonical = "https://quantumbfs.github.io/QuAlgorithmZoo.jl/latest/",
)

deploydocs(
    repo = "github.com/QuantumBFS/QuAlgorithmZoo.jl.git",
    target = "build",
    julia = "1.0",
    deps = nothing,
    make = nothing,
)


text:
using Test
using TimeseriesSurrogates
ENV["GKSwstype"] = "100"

ts = cumsum(randn(1000))
ts_nan =cumsum(randn(100))
ts_nan[1] = NaN

@testset "Constrained surrogates" begin
    @testset "Random shuffle" begin
        surrogate = randomshuffle(ts)
        @test length(ts) == length(surrogate)
        @test all(sort(ts) .== sort(surrogate))
        @test !all(ts .== surrogate)
    end

    @testset "Random phases" begin
        surrogate = randomphases(ts)
        @test length(ts) == length(surrogate)
        @test !all(ts .== surrogate)
    end

    @testset "Random amplitudes" begin
        surrogate = randomphases(ts)
        @test length(ts) == length(surrogate)
        @test all(ts .!= surrogate)
    end

    @testset "AAFT" begin
        surrogate = aaft(ts)
        @test length(ts) == length(surrogate)
        #@test all(ts .!= surrogate)
        @test all(sort(ts) .== sort(surrogate))
	@test_throws DomainError aaft(ts_nan)
    end

    @testset "IAAFT" begin
        # Single realization
        surrogate = iaaft(ts)
        @test length(ts) == length(surrogate)
        #@test all(ts .!= surrogate)
        @test all(sort(ts) .== sort(surrogate))

        # Storing all realizations during iterations (the last vector contains the final
        # surrogate).
        surrogates = iaaft_iters(ts)
        @test length(ts) == length(surrogates[1])
        #@test all(ts .!= surrogates[end])
        @test all(sort(ts) .== sort(surrogates[end]))
	@test_throws DomainError iaaft(ts_nan)
    end

    @testset "WIAAFT" begin
        wiaaft(ts)
    end
end


text:
module PkgSearch

const DEBUG = false

using Requests, Sparklines, MetadataTools, Memoize
DEBUG && using Lumberjack

type SearchState
  page::Int
  terms::String
  packages::Dict{String,Any}
  results::Dict{String,Any}
end

const search_api_uri = "http://genieframework.com/api/v1"
const search_endpoint = "/packages/search"
const details_endpoint = "/packages"

const items_per_page = 10
const search = SearchState(1, "", Dict{String,Any}(), Dict{String,Any}())

"""
  lookup(keyword1, keyword2, ...)

Searches for Julia packages matching the provided keywords. Search results limited to $items_per_page per lookup.
Returns the results and displays the info in the REPL.

#Examples
```julia
julia> PkgSearch.lookup("Redis")

4-element Array{Any,1}:
 Dict{String,Any}("search"=>Dict{String,Any} ... )
```
"""
function lookup(keywords...)
  search.terms = join(map(x -> strip(string(x)), collect(keywords)), "+")
  search.page = 1

  _lookup(search.terms);
end
@memoize function _lookup(keywords::String; autorender::Bool = true, page::Int = 1)
  if isempty(strip(keywords))
    Base.error("Search term can not be empty")
  end
  if search.page <= 0
    Base.error("Invalid page number")
    search.page = 1
  end

  try
    page_uri = search_api_uri * search_endpoint * "?q=" * replace(keywords, " ", "+") * "&page[number]=$(search.page)&page[size]=$items_per_page"
    DEBUG && Lumberjack.info(page_uri)
    search.results = Requests.get(page_uri) |> Requests.json
  catch ex
    rethrow(ex) # do something better with the exception
  end

  if ! isempty(search.results["data"])
    return process_results(search.results, autorender = autorender)
  else
    return []
  end
end

"""
  next([page_jump])

Performs a new lookup to bring the next batch of $items_per_page results, if available. If page_jump is provided, it will go directly to that page.
"""
function next(page_jump::Int = 1)
  search.page += page_jump

  _lookup(search.terms, page = search.page)
end

"""
  prev([page_jump])

Performs a new lookup to bring the previous batch of $items_per_page results, if available. If page_jump is provided, it will go directly to that page.
"""
function prev(page_jump::Int = 1)
  next(page_jump * -1)
end

"""
Returns and optionally displays in the REPL more details about the requested package
"""
function details(package_name::AbstractString; autorender::Bool = true)
  _details(package_name, autorender = autorender)
end
@memoize function _details(package_name::AbstractString; autorender::Bool = true)
  if ! haskey(search.packages, package_name)
    _lookup(package_name, autorender = false)
  end

  package_details = Requests.get(search_api_uri * details_endpoint * "/" * string(search.packages[package_name]["id"])) |> Requests.json

  if isempty(package_details) error("Package not found") end
  autorender && render_package(package_details["data"]["attributes"])

  package_details["data"]["attributes"]
end

"""
Installs the requested package using the appropriate technique: Pkg.add if it's an official METADATA package, or Pkg.clone if not.
"""
function install(package_name::AbstractString)
  if is_official_package(package_name)
    Pkg.add(package_name)
  else
    pkg_details = details(package_name, autorender = false)
    Pkg.clone(pkg_details["url"])
  end
end

"""
Returns true if the requested package is official, that is part of METADATA. False otherwise.
"""
function is_official_package(package_name::AbstractString)
  haskey(metadata_packages(), package_name)
end

"""
Returns the current page number
"""
function current_page()
  search.page
end

"""
Sets the current page number
"""
function current_page(page_number::Int)
  page_number > 0 ? search.page = page_number : error("Invalid page number")
end

"""
Returns the official packages from METADATA
"""
function metadata_packages()
  _metadata_packages()
end
@memoize function _metadata_packages()
  MetadataTools.get_all_pkg()
end

function process_results(search_results::Dict; autorender::Bool = false)
  for p in search_results["data"]["packages"]
    p["search"] = Dict{String,String}()
    p["search"]["headline"] = p["attributes"]["readme"][1:searchindex(p["attributes"]["readme"], "\n", 100)] |> String
    search.packages[p["attributes"]["name"]] = p
    autorender && render(p)
  end

  search_results["data"]
end

function render{T<:AbstractString, U<:Any}(p::Dict{T,U})
  println("=====================================================")
  println(mark_official_package(p["attributes"]["name"]) * " " * p["attributes"]["name"] * " " * p["attributes"]["url"])
  println("-----------------------------------------------------")
  display(p["search"]["headline"] |> Markdown.parse)
  println("_____________________________________________________\n\n")
end

function render_package{T<:AbstractString, U<:Any}(p::Dict{T,U})
  println("=====================================================")
  println(mark_official_package(p["name"]) * " " * p["name"] * " " * p["url"])
  println("-----------------------------------------------------")
  println("Weekly activity: ")
  spark(p["participation"])
  println("")
  println("-----------------------------------------------------")
  display(replace(p["readme"], r"(<b>|</b>)", "**") |> Markdown.parse)
  println("_____________________________________________________\n\n")
end

function mark_official_package(package_name::AbstractString)
  is_official_package(package_name) ? "✅ " : "⚠ "
end

function help()
  println("""\n
    PkgSearch quick options: \n
    * PkgSearch.lookup("words...")
      Returns the first page of results matching the search terms \n
    * PkgSearch.next()
      Returns the next page of search results \n
    * PkgSearch.prev()
      Returns the previous page of search results \n
    * PkgSearch.details("package name")
      Shows the details for "package name" \n
    * Legend:
      ✅   official package (metadata)
      ⚠   unofficial package
  """)
end

println("Preloading METADATA packages... please wait...")
metadata_packages()
println("Ready - use `PkgSearch.help()` for help.")

end

max_stars_repo_name: essenciary/PkgSearch

max_issues_repo_name: essenciary/PkgSearch

max_forks_repo_name: essenciary/PkgSearch


text:
using Combinatorics

# Part 1
ratings = sort(parse.(Int32, readlines("day10_input.txt")))
diffs = vcat(ratings, [ratings[end]+3]) - vcat([0], ratings)
println("Part 1: product = $(count(x -> x==1, diffs) * count(x -> x==3, diffs))")

# Part 2
valid(chain) = maximum([chain; chain[end]+3] - [chain[1]-3; chain]) <= 3
innercombinations(chain) = combinations(chain[2:(lastindex(chain)-1)])
swapguts(chain, inner) = [first(chain); inner; last(chain)]
arrangements(chain) = count([[valid(swapguts(chain, c)) for c in innercombinations(chain)]; valid(swapguts(chain, []))])
function subchains(chain)
    results = []
    current = [first(chain)]
    for i in 2:lastindex(chain)
        if last(current) + 3 == chain[i]
            push!(results, current)
            current = []
        end
        push!(current, chain[i])
    end
    push!(results, current)
    results
end
println("Part 2: arrangements = $(prod(arrangements.(subchains([0; ratings; ratings[end]+3]))))")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
# -------- Datasets ---------------------------------------------------------- #

"""
Structure that holds a list of games and labels, i.e., targets for learning.

Dataset are usually created through playing games from start to finish by an
MCTSPlayer. The value label corresponds to the results of the game, and the
policy label to the improved policy in the single MCTS steps. Furthermore,
features can be enabled for a dataset and are stored as part of the label, as
well.
"""
mutable struct DataSet{G <: Game}

  games  :: Vector{G}                 # games saved in the dataset
  label  :: Vector{Vector{Float32}}   # target value/policy labels
  flabel :: Vector{Vector{Float32}}   # target feature values

  features :: Vector{Feature}         # with which features was the ds created

end

"""
      DataSet{G}([; features])

Initialize an empty dataset for concrete game type `G` and `features` enabled.
"""
function DataSet{G}(; features = Feature[]) where {G <: Game}

  DataSet( Vector{G}()
         , Vector{Vector{Float32}}()
         , Vector{Vector{Float32}}()
         , features )

end

function DataSet(games, label, flabel; features = Feature[])

  DataSet(games, label, flabel, features)

end

features(ds :: DataSet) = ds.features
Base.length(d :: DataSet) = length(d.games)

# -------- Dataset Operations ------------------------------------------------ #

function Base.merge(d :: DataSet{G}, ds...) where {G <: Game}

  # Make sure that all datasets have compatible features
  features = d.features

  if !all(x -> x.features == features, ds)

    error("Merging databases with incompatible features")

  end

  # Create and return the merged dataset
  dataset = DataSet{G}()
  dataset.features = features

  dataset.games  = vcat([d.games,  (x.games  for x in ds)...]...)
  dataset.label  = vcat([d.label,  (x.label  for x in ds)...]...)
  dataset.flabel = vcat([d.flabel, (x.flabel for x in ds)...]...)

  dataset

end

function Base.split(d :: DataSet{G}, size :: Int; shuffle = true) where {G}

  n = length(d)
  @assert size <= n "Cannot split dataset of length $n at position $size."

  idx = shuffle ? randperm(n) : 1:n
  idx1, idx2 = idx[1:size], idx[size+1:end]

  d1 = DataSet( d.games[idx1], d.label[idx1]
              , d.flabel[idx1], features = d.features)
  d2 = DataSet( d.games[idx2], d.label[idx2]
              , d.flabel[idx2], features = d.features)

  d1, d2

end

function augment(d :: DataSet{G}) :: Vector{DataSet{G}} where G <: Game

  # NOTE: Augmentation will render flabel information useless.
  # Therefore, one must recalculate them after applying this function!

  # Problem: we must augment d in such a way that playthroughs are still
  # discernible after augmentation. I.e., for each symmetry transformation, one
  # DataSet will be returned.

  gs, ls = unzip([augment(g, l) for (g, l) in zip(d.games, d.label)])

  map(1:length(gs[1])) do j

    games = map(x -> x[j], gs)
    label = map(x -> x[j], ls)

    DataSet(games, label, copy(d.flabel), features = d.features)

  end

end

# -------- Raw DataSet Representation: Caches -------------------------------- #

struct DataCache{G <: Game, GPU}

  data        # game representation data

  vlabel      # value labels
  plabel      # policy labels
  flabel      # feature labels

end

function DataCache( ds :: DataSet{G}
                  ; gpu = false
                  , use_features = false
                  ) where {G <: Game}

  # Preparation
  at = atype(gpu)
  vplabel = hcat(ds.label...)

  # Convert to at
  data = convert(at, representation(ds.games))
  vlabel = convert(at, vplabel[1, :])
  plabel = convert(at, vplabel[2:end, :])
  flabel = use_features ? convert(at, hcat(ds.flabel...)) : nothing

  DataCache{G, gpu}(data, vlabel, plabel, flabel)

end

Base.length(c :: DataCache) = size(c.data)[end]

# -------- Iterating Datasets: Batches --------------------------------------- #

struct Batches{G <: Game, GPU}

  cache :: DataCache{G, GPU}

  batchsize :: Int
  shuffle :: Bool
  partial :: Bool

  indices :: Vector{Int}

end

function Batches( d :: DataSet{G}
                , batchsize
                ; shuffle = false
                , partial = true
                , gpu = false
                , use_features = false
                ) where {G <: Game}

  indices = collect(1:length(d))

  cache = DataCache(d, gpu = gpu, use_features = use_features)

  Batches{G, gpu}(cache, batchsize, shuffle, partial, indices)

end

function Base.length(b :: Batches)
  n = length(b.cache) / b.batchsize
  b.partial ? ceil(Int, n) : floor(Int, n)
end

function Base.iterate(b :: Batches{G, GPU}, start = 1) where {G <: Game, GPU}

  # Preparations
  l = length(b.cache)
  b.shuffle && start == 1 && (b.indices .= randperm(l))

  # start:stop is the range in b.indices that selected
  stop = min(start + b.batchsize - 1, l)

  # Check for end of iteration
  if start > l || !b.partial && stop - start < b.batchsize - 1
    return nothing
  end

  # Build the data cache
  idx = b.indices[start:stop]

  data = b.cache.data[:, :, :, idx]
  vlabel = b.cache.vlabel[idx]
  plabel = b.cache.plabel[:, idx]
  flabel = isnothing(b.cache.flabel) ? nothing : b.cache.flabel[:, idx]

  cache = DataCache{G, GPU}(data, vlabel, plabel, flabel)

  # Return the (cache, new_start) state tuple
  cache, stop + 1

end

# -------- Generating Datasets: Helpers -------------------------------------- #

function _record( play :: Function # maps game to dataset
                , root :: G
                , n :: Int
                ; features
                , ntasks
                , augment = true
                , callback = () -> nothing
                , merge = true
                , prepare = prepare(steps = 0)
                , branch = branch(prob = 0, steps = 1)
                ) where {G <: Game}

  # Extend the provided play function by random branching
  # and the call to callback
  play_with_branching = _ -> begin

    # Play the (prepared) game once without branching inbetween
    dataset = play(prepare(root))

    # Create branchpoints and play them, too
    branchpoints = filter(!isnothing, branch.(dataset.games[1:end-1]))
    branches = play.(branchpoints)

    # Filter branches where we were directly lead to the end of the game
    filter!(x -> length(x) > 0, branches)

    # Collect the original and the branched datasets
    datasets = [dataset; branches]

    # Augment the datasets (if desired)
    augment && (datasets = mapreduce(Jtac.augment, vcat, datasets))

    # Calculate features (if there are any)
    datasets = map(d -> _record_features!(d, features), datasets)

    # Signal that we are done with one iteration
    callback()

    Base.merge(datasets...)

  end

  # Call play_with_branching n times
  if ntasks == 1 || n == 1

    # Asyncmap is assumed to be of disadvantage
    ds = map(play_with_branching, 1:n)

  else

    # Asyncmap is assumed to be of advantage
    ds = asyncmap(play_with_branching, 1:n, ntasks = ntasks)

  end

  # Mark the given features as activated
  for d in ds d.features = features end

  # Return one big dataset or the different playings separately?
  merge ? Base.merge(ds...) : ds

end

function _record_move!( dataset :: DataSet{G}
                      , game
                      , player
                      ) where {G <: Game}

  pl = policy_length(G)

  # Record the current game state
  push!(dataset.games, copy(game))

  # Get the improved policy from the player
  policy = think(player, game)

  # Advance the game randomly according to the policy
  apply_action!(game, choose_index(policy))

  # Prepare the label for this game state with policy filled in
  label = zeros(Float32, 1 + pl)
  label[2:pl+1] .= policy

  push!(dataset.label, label)

end

function _record_final!( dataset :: DataSet{G}
                      , game ) where {G <: Game}

  pl = policy_length(G)
  push!(dataset.games, copy(game))
  push!(dataset.label, zeros(Float32, 1 + pl))

end

function _record_value!(dataset :: DataSet{G}) where {G <: Game}

  result = status(dataset.games[end])

  for i = 1:length(dataset.games)

    dataset.label[i][1] = current_player(dataset.games[i]) * result

  end

end

# Add feature label to an extended dataset (that contains the final game state)
function _record_features!( dataset :: DataSet{G}
                          , features
                          ) where {G <: Game}

  fl = feature_length(features, G)

  dataset.flabel = Array{Float32}[]

  # Go over the unextended dataset
  for i = 1:(length(dataset.games)-1)

    j = 0
    flabel = zeros(Float32, fl)

    # Set all features
    for f in features
      l = feature_length(f, G)
      flabel[j+1:j+l] = f(dataset.games[i], dataset.games)
      j += l
    end

    push!(dataset.flabel, flabel)

  end

  # Remove the final game (and the corresponding zero-label) from the dataset
  # now that the feature labels are calculated.
  pop!(dataset.games)
  pop!(dataset.label)

  dataset

end

# -------- Generating Datasets: API ------------------------------------------ #

"""
    record_self(player [, n]; <keyword arguments>)

Record a dataset by letting `player` play against itself `n` times.

The value label for the dataset is derived from the outcome of the game (-1, 0,
or +1). The policy label is the policy proposed by the player for the respective
game state. The feature labels are calculated by applying the provided features
(see keyword arguments) to each stored game state and the corresponding history.

# Arguments
- `game`: Initial game state that is compatible with `player`. Derived by default.
- `features`: List of features for which labels are created. Derived by default.
- `augment = true`: Whether to apply augmentation on the generated dataset.
- `merge = true`: Whether to return one merged dataset or `n` seperate ones.
- `callback`: Procedure that is called afer each completed playing.
- `prepare`: Function applied on `game` that returns the initial state for
   a single selfplay. Can be used to randomize the initial positions, for
   example by random actions. Default preparation maps can be generated by
   `prepare(; steps).
- `branch`: Function applied to each game state of a first selfplay. It can
   return `nothing` or branched game states, which are then used as roots for
   new playings (without recursive branching). Default branching maps can be
   generated by `branch(; prob, steps)`.
- `distributed = false`: Whether to conduct the self-playings on several
   processes in parallel. If `true`, all available workers are used. Alternatively,
   a list of worker ids can be passed.
- `tickets`: Number of chunks in which the workload is distributed if
   `distributed != false`. By default, it is set to the number of workers.

# Examples
```julia
# Record 20 self-playings of an classical MCTS player with power 250
player = MCTSPlayer(power = 250)
dataset = record_self(player, 20, game = TicTacToe(), branch = branch(0.25))

# Record 10 self-playings of MCTS player with shallow predictor network and
# power 50
model = NeuralModel(TicTacToe, @chain TicTacToe Dense(50))
player = MCTSPlayer(model, power = 50)
dataset = record_self(player, 10, augment = false)
```
"""
function record_self( p :: Player{G}
                    , n :: Int = 1
                    ; game :: T = G()
                    , features = features(p)
                    , distributed = false
                    , tickets = nothing
                    , kwargs...
                    ) where {G, T <: G}

  # If self-playing is to be distributed, get the list of workers and
  # cede to the corresponding distributed function, defined below
  if distributed != false

    workers = distributed == true ? Distributed.workers() : distributed
    tickets = isnothing(tickets) ? length(workers) : tickets

    return record_self_distributed( p, n
                                  ; game = game
                                  , features = features
                                  , workers = workers
                                  , tickets = tickets
                                  , kwargs... )

  end

  # Function that plays a single game of player against itself
  play = game -> begin

    dataset = DataSet{T}()

    while !is_over(game)

      # Play one move, prepare the label, and record the game and player policy
      _record_move!(dataset, game, p)

    end

    # Here, we push the final game state to the dataset.
    # This game is needed for calculating the features later, but it has to be
    # removed again after that. See _record_features!
    _record_final!(dataset, game)

    # Complete the label value information now that we know the total game history
    _record_value!(dataset)

    # Note: the features are calculated later, after augmentation steps took
    # place in _record below

    # Return the value-policy labeled dataset
    dataset

  end

  # Record several games, potentially with branching
  _record(play, game, n; features = features, ntasks = ntasks(p), kwargs...)

end

"""
    record_against(player, enemy [, n]; <keyword arguments>)

Record a dataset by letting `player` play against `enemy` a number of `n` times.
The datasets are recorded similarly to `record_self`.

# Arguments
The function takes the following arguments:
- `game`: Initial game state that is compatible with `player`.
- `features`: List of features for which feature labels are created.
- `start`: Function that determines the starting player (-1: enemy, 1: player).
- `augment`: Whether to apply symmetry augmentation on the generated dataset.
- `merge = false`: Whether to return one merged dataset or seperate playings.
- `callback`: Function that is called afer each finished game.
- `prepare`: Function applied on `game` that returns the initial state for
   a single playing. Can be used to randomize the initial positions, for
   example by random actions. Default preparation maps can be generated by
   `prepare(; steps).
- `branch`: Function applied to each game state of a first playing. It can
   return `nothing` or branched game states, which are then used as roots for
   new playings (without recursive branching). Default branching maps can be
   generated by `branch(; prob, steps)`.
- `distributed = false`: Whether to conduct the playings on several processes in
   parallel. If `true`, all available workers are used. Alternatively, a list of
   worker ids can be passed.
- `tickets`: Number of chunks in which the workload is distributed if
   `distributed != false`. By default, it is set to the number of workers.
"""
function record_against( p :: Player{G}
                       , enemy :: Player{H}
                       , n :: Int = 1
                       ; game :: T = G()
                       , start :: Function = () -> rand([-1, 1])
                       , features = features(p)
                       , distributed = false
                       , tickets = nothing
                       , kwargs...
                       ) where {G, H, T <: typeintersect(G, H)}

  # If playing is to be distributed, get the list of workers and
  # cede to the corresponding function in distributed.jl
  if distributed != false

    workers = distributed == true ? Distributed.workers() : distributed
    tickets = isnothing(tickets) ? length(workers) : tickets

    return record_against_distributed( p, enemy, n
                                     ; game = game
                                     , start = start
                                     , features = features
                                     , workers = workers
                                     , tickets = tickets
                                     , kwargs... )

  end

  play = game -> begin

    dataset = DataSet{T}()

    # Roll the dice to see who starts. Our player (i = 1) or the enemy (i = -1)?
    s = start()

    while !is_over(game)

      # If it is our turn, record the move. Otherwise, let the enemy play
      play = (current_player(game) == s)
      play ? _record_move!(dataset, game, p) : turn!(game, enemy)

    end

    # Finish the dataset (except for features)
    _record_final!(dataset, game)
    _record_value!(dataset)

    dataset

  end

  ntasks = Jtac.ntasks(p) + Jtac.ntasks(enemy) - 1

  _record(play, game, n; features = features, ntasks = ntasks, kwargs...)

end

"""
    record_model(model, games [; use_features])

Record a dataset by applying `model` to `games`. The features enabled for
`model` are recorded if `use_features` is true. `callback` is a function that is
called after each application of `model` to a game state.
"""
function record_model( model :: Model{G}
                     , games :: Vector{T}
                     ; use_features = true
                     , augment = false
                     ) where {G, T <: G}

  games = augment ? mapreduce(Jtac.augment, vcat, games) : games

  vplabel = Vector{Float32}[]
  flabel  = Vector{Float32}[]

  v, p, f = apply_features(model, game)

  for i in 1:length(games)
    push!(vplabel, vcat(v[i], p[:,i]))
    use_features ? push!(flabel, f[:,i]) : push!(flabel, similar(f[:,i], 0))
  end

  features = use_features ? Jtac.features(model) : Feature[]

  DataSet(games, vplabel, flabel, features = features)

end

# -------- Distributed Recording --------------------------------------------- #

function record_self_distributed( p :: Player
                                , n :: Int = 1
                                ; workers = workers()
                                , merge = false
                                , kwargs... )

  # Create the record function
  record = (ps, n; kwargs...) -> begin
    record_self(ps[1], n; merge = merge, kwargs...)
  end

  # Use the with_workers function defined in src/distributed.jl
  ds = with_workers(record, [p], n; workers = workers, kwargs...)
  ds = vcat(ds...)

  merge ? Base.merge(ds) : ds

end

function record_against_distributed( p :: Player
                                   , enemy :: Player
                                   , n :: Int = 1
                                   ; workers = workers()
                                   , merge = false
                                   , kwargs... )

  # Create the record function
  record = (ps, n; kwargs...) -> begin
    record_against(ps[1], ps[2], n; merge = merge, kwargs...)
  end

  # Use the with_workers function defined in src/distributed.jl
  ds = with_workers(record, [p, enemy], n; workers = workers, kwargs...)
  ds = vcat(ds...)

  merge ? Base.merge(ds) : ds

end

max_stars_repo_name: roSievers/jtac


text:
import Pkg
Pkg.activate(".")

using KMarkovGenSeq
println("Using CPU")
for i in 3:9
    @time query_freq("data/test.fa", "data/genomes/", i, nothing, false)
end

println("Using GPU")
for i in 3:9
    @time query_freq("data/test.fa", "data/genomes/", i, nothing, true)
end


text:
"""
Compute the inverse of the lattice. Takes special care of 1D or 2D cases.
"""
function compute_inverse_lattice(lattice::AbstractMatrix{T}) where {T}
    # Note: pinv pretty much does the same, but the implied SVD causes trouble
    #       with interval arithmetic and dual numbers, so we go for this version.
    n_dim = count(!iszero, eachcol(lattice))
    @assert 1 ≤ n_dim ≤ 3
    if n_dim == 3
        inv(lattice)
    else
        Mat3{T}([
            inv(lattice[1:n_dim, 1:n_dim])   zeros(T, n_dim, 3 - n_dim);
            zeros(T, 3 - n_dim, 3)
        ])
    end
end

"""
Compute the reciprocal lattice.
We use the convention that the reciprocal lattice is the set of G vectors such
that G ⋅ R ∈ 2π ℤ for all R in the lattice.
"""
function compute_recip_lattice(lattice::AbstractMatrix{T}) where {T}
    2T(π) * compute_inverse_lattice(lattice')
end

"""
Compute unit cell volume volume. In case of 1D or 2D case, the volume is the length/surface.
"""
function compute_unit_cell_volume(lattice)
    n_dim = count(!iszero, eachcol(lattice))
    abs(det(lattice[1:n_dim, 1:n_dim]))
end

"""Compute the diameter of the unit cell"""
function diameter(lattice::AbstractMatrix)
    # brute force search
    diam = zero(eltype(lattice))
    for vec in Vec3.(Iterators.product(-1:1, -1:1, -1:1))
        diam = max(diam, norm(lattice * vec))
    end
    diam
end

"""
Estimate integer bounds for dense space loops from a given inequality ||Mx|| ≤ δ.
For 1D and 2D systems the limit will be zero in the auxiliary dimensions.
"""
function estimate_integer_lattice_bounds(M::AbstractMatrix{T}, δ, shift=zeros(3); tol=sqrt(eps(T))) where T
    # As a general statement, with M a lattice matrix, then if ||Mx|| <= δ,
    # then xi = <ei, M^-1 Mx> = <M^-T ei, Mx> <= ||M^-T ei|| δ.
    inv_lattice_t = compute_inverse_lattice(M')
    xlims = [norm(inv_lattice_t[:, i]) * δ + shift[i] for i in 1:3]

    # Round up, unless exactly zero (in which case keep it zero in
    # order to just have one x vector for 1D or 2D systems)
    xlims = [xlim == 0 ? 0 : ceil(Int, xlim .- tol) for xlim in xlims]
    xlims
end