text:
module Crazyflie

using PyCall

export scan, connect, disconnect, play
export bootloader, crtp, drivers, positioning, utils
export crazyflie, synccrazyflie, logger, synclogger, mem
export record, log_posn, log_quat, log_config
#  export reset_estimator, upload_traj, enable_commander
#  export land

export motor_ramp_test, log_posn, log_quat

const bootloader = PyNULL()
const crtp = PyNULL()
const drivers = PyNULL()
const positioning = PyNULL()
const utils = PyNULL()
const crazyflie = PyNULL()
const synccrazyflie = PyNULL()
const logger = PyNULL()
const synclogger = PyNULL()
const mem = PyNULL()

function __init__()
    copy!(bootloader, pyimport("cflib.bootloader"))
    copy!(crtp, pyimport("cflib.crtp"))
    copy!(drivers, pyimport("cflib.drivers"))
    copy!(positioning, pyimport("cflib.positioning"))
    copy!(utils, pyimport("cflib.utils"))
    copy!(crazyflie, pyimport("cflib.crazyflie"))
    copy!(synccrazyflie, pyimport("cflib.crazyflie.syncCrazyflie"))
    copy!(logger, pyimport("cflib.crazyflie.log"))
    copy!(synclogger, pyimport("cflib.crazyflie.syncLogger"))
    copy!(mem, pyimport("cflib.crazyflie.mem"))

    crtp.init_drivers(enable_debug_driver=false)
    nothing
end

function scan()
    available = crtp.scan_interfaces()
    if isempty(available)
        print("No crazyflies found.\n")
    end
    print("Found $(size(available, 1)) crazyflies:\n")
    for device in eachrow(available)
        print("\t" * device[1] * "\n")
    end
    nothing
end

function _first_available()
    available = crtp.scan_interfaces()
    if isempty(available)
        error("No crazyflies found. " *
              "Please specifiy a URI or try again. \n")
    end
    return available[1,1]
end

function connect(uri=_first_available())
    scf = synccrazyflie.SyncCrazyflie(uri, cf=crazyflie.Crazyflie(rw_cache="./cache"))
    scf.open_link()
    return scf
end

function play(alg, uri=_first_available())
    scf = connect(uri)
    try
        alg(scf.cf)
    catch e
        @show e
        # Catch keyboard interrupts if desired
    end
    disconnect(scf)
end

function disconnect(scf)
    scf.close_link()
    nothing
end
#
#  function upload_traj(cf, id, traj)
#      trajmem = cf.mem.get_mems(mem.MemoryElement.TYPE_TRAJ)[1]
#      total_duration = 0
#      for row in traj
#          duration = row[1]
#          x = mem.Poly4D.Poly(row[2:9])
#          y = mem.Poly4D.Poly(row[10:17])
#          z = mem.Poly4D.Poly(row[18:25])
#          yaw = mem.Poly4D.Poly(row[26:33])
#          push!(trajmem.poly4Ds, mem.Poly4D(duration, x, y, z, yaw))
#          total_duration += duration
#      end
#
#      py"""
#      import time
#      class Uploader:
#          def __init__(self):
#              self._is_done = False
#
#          def upload(self, trajectory_mem):
#              print('Uploading data')
#              trajectory_mem.write_data(self._upload_done)
#
#              while not self._is_done:
#                  time.sleep(0.2)
#
#          def _upload_done(self, mem, addr):
#              print('Data uploaded')
#              self._is_done = True
#      """
#      py"Uploader().upload($trajmem)"
#      cf.high_level_commander.define_trajectory(id, 0, length(trajmem.poly4Ds))
#      return total_duration
#  end
#
#  function reset_estimator(cf)
#      cf.param.set_value("kalman.resetEstimation", "1")
#      sleep(0.1)
#      cf.param.set_value("kalman.resetEstimation", "0")
#      lconfig = logger.LogConfig(name="Kalman Variance", period_in_ms=500)
#      lconfig.add_variable("kalman.varPX", "float")
#      lconfig.add_variable("kalman.varPY", "float")
#      lconfig.add_variable("kalman.varPZ", "float")
#      varx = ones(10)*1000
#      vary = ones(10)*1000
#      varz = ones(10)*1000
#      record(cf, lconfig) do logs
#          for entry in logs
#              data = entry[2]
#              push!(varx, data["kalman.varPX"])
#              push!(vary, data["kalman.varPY"])
#              push!(varz, data["kalman.varPZ"])
#              popfirst!(varx); popfirst!(vary); popfirst!(varz)
#              threshx = (maximum(varx) - minimum(varx)) < 0.001
#              threshy = (maximum(vary) - minimum(vary)) < 0.001
#              threshz = (maximum(varz) - minimum(varz)) < 0.001
#              if threshx && threshy && threshz
#                  break
#              end
#          end
#      end
#  end
#
#  function enable_commander(cf)
#      cf.param.set_value("commander.enHighLevel", "1")
#      cf.high_level_commander
#  end
#
#  # commander take off because commander.land is not reliable
#  function land(commander, duration)
#      print("Landing\n")
#      for i = 1:10
#          # keep sending land messages because
#          # it doesn't recieve it sometimes
#          commander.land(0.0, duration)
#          sleep(0.1)
#      end
#      sleep(duration)
#      commander.stop()
#  end

include("log.jl")
include("examples.jl")

end # module


text:
using DReal

ctx = Context(qf_nra)
X = Var(ctx, Float64,"Xleaks2", 0.0,1.0)
add!(ctx,(<=)(ctx,X,0.3))
for i = 1:1000
  push_ctx!(ctx)
  a = rand()
  b = rand()
  add!(ctx,(>=)(ctx,X,0.1))
  is_satisfiable(ctx)
  pop_ctx!(ctx)
end

ext: jl


text:
function show_me(df::DataFrame, var=:peakArea; logy=true)
  figure()
  if logy
    semilogy(df[:date], df[var], "ok")
  else
    plot(df[:date], df[var], "ok")
  end
  grid(true)
  ylabel(string(var))
end

function show_me(dfs::Vector{DataFrame}, var=:peakArea; logy=true)
  figure()
  for df in dfs
    if logy
      semilogy(df[:date], df[var], "o", lw=2)
    else
      plot(df[:date], df[var], "o", lw=2)
    end
  end
  ylabel(string(var))
  grid(true)
end

function show_me(s::Spectrum; logy=true)
  if logy
    figure()
    semilogy(s.baseline_A, "ok", markerfacecolor="white")
    semilogy(s.baseline_B, "or", markerfacecolor="white")
    semilogy(s.row_A, "-k", lw=2, label="Row A")
    semilogy(s.row_B, "-r", lw=2, label="Row B")
    legend()

    figure()
    semilogy(s.ionsPerSpectrum_A, lw=2, label="ions per spectrum A")
    semilogy(s.ionsPerSpectrum_B, lw=2, label="ions per spectrum B")
    semilogy(s.peakIndices_A, s.peakAmplitudes_A, "or")
    semilogy(s.peakIndices_B, s.peakAmplitudes_B, "sr")
  else
    figure()
    plot(s.baseline_A, "ok", markerfacecolor="white")
    plot(s.baseline_B, "or", markerfacecolor="white")
    plot(s.row_A, "-k", lw=2, label="Row A")
    plot(s.row_B, "-r", lw=2, label="Row B")
    legend()

    figure()
    plot(s.ionsPerSpectrum_A, lw=2, label="ions per spectrum A")
    plot(s.ionsPerSpectrum_B, lw=2, label="ions per spectrum B")
    plot(s.peakIndices_A, s.peakAmplitudes_A, "or")
    plot(s.peakIndices_B, s.peakAmplitudes_B, "sr")

  end
  grid(true)
  legend()
end

function show_me(fileName::ASCIIString; logy=false)
  y, gainStep, t, m0, p_pt = parseDataFileBothRows(fileName)
  @show(gainStep)
  @show(t)
  @show(m0)
  @show(p_pt)
  figure()
  if logy
    semilogy(y[:,1], "-ok", label="Row A")
    semilogy(y[:,2], "-sr", label="Row B")
  else
    plot(y[:,1], "-ok", label="Row A")
    plot(y[:,2], "-sr", label="Row B")
  end
  grid(true)
  legend()
end

function show_me(y::Vector{Real}; logy=true)
  figure()
  if logy
    semilogy(y, "-k", lw=2)
  else
    plot(y, "-k", lw=2)
  end
  grid(true)
  nothing
end


text:
struct DurationResolver{T₁ <: Number} <: Function
    intersample_duration::T₁
    nested_interval::NestedInterval
end

function (resolve_duration::DurationResolver)(n::Nothing)
    # All times are in milliseconds
    nested_interval = resolve_duration.nested_interval
    tₛ = (get_stop(nested_interval) - get_start(nested_interval)) * resolve_duration.intersample_duration
    # Display timestamp.
    #timestamp = Dates.epochms2datetime(Dates.value(Millisecond(tₛ)))
    #timestr = Dates.format(timestamp , dateformat"HH:MM:SS")
    timestr = Dates.canonicalize(Dates.CompoundPeriod(Dates.Millisecond(tₛ)))
    string(timestr)
end


text:
using HTTP: HTTP, request, URI
using JSON3: JSON3
using XLSX: XLSX, readxlsx, readtable, readdata, sheetnames
using CSV: CSV
using DataFrames: DataFrames, DataFrame, select!, unstack,
                  subset!, ByRow
const BEA_API_KEY = get(ENV, "api_bea_token", "")

#=
Find the mapping to of industry to NAICS2 (sector) through the APIs
=#
url = URI(scheme = "https", host = "apps.bea.gov", path = "/api/data",
          query = ["UserID" => BEA_API_KEY,
                   "method" => "GetData",
                   "DatasetName" => "InputOutput",
                   "TableID" => 259,
                   "Year" => join(2008:2019, ','),
                   "ResultFormat" => "JSON"
                  ])
response = request("GET", url)
@assert response.status == 200
json = JSON3.read(response.body)
x = DataFrame((year = row.Year, descr = row.RowDescr, value = parse(Int, row.DataValue))
                for row in json.BEAAPI.Results.Data
                  if (row.RowDescr ∈ ["Total Intermediate",
                                      "Compensation of employees",
                                      "Gross operating surplus",
                                      "Other taxes on production",
                                      ""]) &
                     isequal("5415", row.ColCode))
y = unstack(x, :descr, :value)
transform!(y,
           ["Total Intermediate", "Compensation of employees"] => ByRow(/) => :ti_coe,
           ["Gross operating surplus", "Compensation of employees"] => ByRow(/) => :gos_coe,
           ["Other taxes on production", "Compensation of employees"] => ByRow(/) => :top_coe)
select!(y, [:year, :ti_coe, :gos_coe, :top_coe])
CSV.write(joinpath("data", "usetbl5415.csv"), y)

max_stars_repo_name: Nosferican/Estimating-Resource-Costs

max_issues_repo_name: Nosferican/Estimating-Resource-Costs

max_forks_repo_name: Nosferican/Estimating-Resource-Costs


text:
using Oceananigans.Architectures
using Oceananigans.BoundaryConditions
using Oceananigans.TurbulenceClosures: calculate_diffusivities!
using Oceananigans.ImmersedBoundaries: mask_immersed_field!

import Oceananigans.TimeSteppers: update_state!

"""
    update_state!(model::NonhydrostaticModel)

Update peripheral aspects of the model (halo regions, diffusivities, hydrostatic pressure) to the current model state.
"""
function update_state!(model::NonhydrostaticModel)

    # Mask immersed tracers
    tracer_masking_events = Tuple(mask_immersed_field!(c) for c in model.tracers)

    wait(device(model.architecture), MultiEvent(tracer_masking_events))

    # Fill halos for velocities and tracers
    fill_halo_regions!(merge(model.velocities, model.tracers), model.architecture,  model.clock, fields(model))

    # Calculate diffusivities
    calculate_diffusivities!(model.diffusivity_fields, model.closure, model)
    fill_halo_regions!(model.diffusivity_fields, model.architecture, model.clock, fields(model))

    # Calculate hydrostatic pressure
    pressure_calculation = launch!(model.architecture, model.grid, :xy, update_hydrostatic_pressure!,
                                   model.pressures.pHY′, model.grid, model.buoyancy, model.tracers,
                                   dependencies=Event(device(model.architecture)))

    # Fill halo regions for pressure
    wait(device(model.architecture), pressure_calculation)

    fill_halo_regions!(model.pressures.pHY′, model.architecture)

    return nothing
end


text:
# Adds packages that don't require special setup
using Pkg
#Pkg.develop(path = joinpath(@__DIR__, "McCormick.jl-master"))
Pkg.develop(path = joinpath(@__DIR__, "McCormick.jl"))
Pkg.develop(path = joinpath(@__DIR__, "MINLPLib.jl"))
Pkg.develop(path = joinpath(@__DIR__, "EAGO.jl-SIPextension"))
Pkg.add("JSON"); Pkg.add("DataFrames"); Pkg.add("CSV")
Pkg.add("JuMP"); Pkg.add("IntervalArithmetic"); Pkg.add("LaTeXStrings")
Pkg.add("Plots"); Pkg.add("StatsBase"); Pkg.add("BenchmarkProfiles")

# Loads solver_benchmarking module
include(joinpath(@__DIR__, "solver_benchmarking.jl"))

# Loads relevent modules
using JuMP, MINLPLib, SCIP, EAGO, IntervalArithmetic, GAMS
using CSV, DataFrames, LaTeXStrings, Plots, StatsBase, BenchmarkProfiles

function baron_factory()
    m = GAMS.Optimizer(GAMS.GAMSWorkspace("C:\\GAMS\\37"))

    m.gams_options["nlp"] = "BARON"
    m.gams_options["optca"] = 1E-4
    m.gams_options["optcr"] = 1E-4

    return m
end

function scip_factory()
    m = SCIP.Optimizer()
    MOI.set(m, SCIP.Param("limits/gap"), 1E-4)
    MOI.set(m, SCIP.Param("limits/absgap"), 1E-4)
    m
end

function eago_factory()
    m = EAGO.Optimizer()
    MOI.set(m, MOI.RawOptimizerAttribute("absolute_tolerance"), 1E-4)
    MOI.set(m, MOI.RawOptimizerAttribute("relative_tolerance"), 1E-4)
    MOI.set(m, MOI.RawOptimizerAttribute("mul_relax_style"), 0)
    m
end

function write_new_model(model_file_name, f, W, b, n, l, p, flag)

    # Gets interval lower bound
    function ANN(w, b, x)
        y = b[1] + W[1]*x
        for i = 2:length(b)
            y = b[i] + W[2]*y
        end
        return sum(y)
    end
    X = [Interval(-1,1) for i=1:n]
    qIntv = ANN(W, b, X)
    qL = qIntv.lo
    qU = qIntv.hi

    # Builds nn
    layer_val = [String[] for i=1:l]
    for k = 1:l               # compute first layer
        bl = b[k]
        Wl = W[k]
        for i = 1:p           # iterate over number of neurons per layer
            s = "$(bl[i])"
            if k != 1
                llayer = layer_val[k-1]
                for j = 1:p
                    s = s*" + $(Wl[i,j])*$(llayer[j])"
                end
            else
                for j = 1:n
                    s = s*" + $(Wl[i,j])*\$(x[$j])"
                end
            end
            s = f(s; lib = flag)
            push!(layer_val[k], s)
        end
    end

    llast = layer_val[l]
    nn = "$(llast[1])"
    for i = 2:p
        nn = nn*" + $(llast[i])"
    end

    open(model_file_name, "w") do file
        write(file, "using JuMP, EAGO\n
                     m = Model()\n
                     EAGO.register_eago_operators!(m)\n
                     @variable(m, -1 <= x[i=1:$n] <= 1)\n
                     @variable(m, $qL <= q <= $qU)\n
                     add_NL_constraint(m, :("*nn*" - \$q <= 0.0))\n
                     @objective(m, Min, q)\n
                     return m\n
                    "
              )
    end
    println("wrote instance to $model_file_name")
end

relu(x::String; lib = true) =     lib ? "relu("*x*")"     : "max("*x*", 0.0)"
silu(x::String; lib = true) =     lib ? "swish("*x*")"    : "("*x*")/(1 + exp(-("*x*")))"
gelu(x::String; lib = true) =     lib ? "gelu("*x*")"     : "("*x*")*(1 + erf(("*x*")/sqrt(2)))/2"
softsign(x::String; lib = true) = lib ? "softsign("*x*")" : "("*x*")/(1 + "*x*")"
sigmoid(x::String; lib = true) =  lib ? "sigmoid("*x*")"  : "1/(1 + exp(-("*x*")))"
softplus(x::String; lib = true) = lib ? "softplus("*x*")" : "log(1 + exp("*x*"))"
maxsig(x::String; lib = true) =   lib ? "maxsig("*x*")"   : "max("*x*", 1/(1 + exp(-("*x*"))))"

function create_lib()
    instance_number = 100
    variable_range = 2:5
    layer_range = 1:4
    neuron_per_layer_range = 2:5
    act_func = [sigmoid; silu; gelu; softplus]
    minlp_folder_instances = joinpath(@__DIR__, "MINLPLib.jl", "instances")

    for i = 1:instance_number
        n = rand(variable_range)
        layers = rand(layer_range)
        neuron_per_layer = rand(neuron_per_layer_range)
        W = []
        b = []
        for j = 1:layers
            if j == 1
                push!(W, rand(neuron_per_layer, n)*2 .- 1)
            else
                push!(W, rand(neuron_per_layer, neuron_per_layer)*2 .- 1)
            end
            push!(b, rand(neuron_per_layer)*2 .- 1)
        end
        j = i > 9 ? "$i" : "0$i"
        for f in act_func
            file_name_env = joinpath(minlp_folder_instances, "ANN_Env", "$(j)_"*String(Symbol(f))*"_$(n)_$(layers)_$(neuron_per_layer).jl")
            file_name_expr = joinpath(minlp_folder_instances, "ANN_Expr", "$(j)_"*String(Symbol(f))*"_$(n)_$(layers)_$(neuron_per_layer).jl")
            file_name_env_15 = joinpath(minlp_folder_instances, "ANN_Env_15", "$(j)_"*String(Symbol(f))*"_$(n)_$(layers)_$(neuron_per_layer).jl")
            file_name_expr_15 = joinpath(minlp_folder_instances, "ANN_Expr_15", "$(j)_"*String(Symbol(f))*"_$(n)_$(layers)_$(neuron_per_layer).jl")

            write_new_model(file_name_env, f, W, b, n, layers, neuron_per_layer, true)
            write_new_model(file_name_expr, f, W, b, n, layers, neuron_per_layer, false)
            write_new_model(file_name_env_15, f, W, b, n, layers, neuron_per_layer, true)
            write_new_model(file_name_expr_15, f, W, b, n, layers, neuron_per_layer, false)
        end
    end
end

# set create_lib_files to true to generate new ANNs in the problem libary
create_lib_files = false
create_lib_files && create_lib()

expr_solvers = Dict{String,Any}()
expr_solvers["SCIP"]  = scip_factory
expr_solvers["EAGO"]  = eago_factory
expr_solvers["BARON"] = baron_factory

env_solvers = Dict{String,Any}()
env_solvers["EAGO"] = eago_factory

function drop_num_underscore(x)
    y = replace(x, "_" => "")
    for i=0:9
        y = replace(y, "$i" => "")
    end
    return y
end

function print_summary_tables(df_env, df_expr, fig_name, env_lib, upper_limit)

    df_expr.Form = ["Exp" for i=1:size(df_expr,1)]
    df_env.Form = ["Env" for i=1:size(df_env,1)]

    df_comb = vcat(df_env, df_expr)
    df_comb[!,:FullSolverName] = map((x,y) -> string(x)*" "*string(y), df_comb.SolverName, df_comb.Form)
    df_comb.SolveTime = map((x,y) -> ifelse(occursin("INFEASIBLE",x), upper_limit, y), df_comb.TerminationStatus, df_comb.SolveTime)
    df_comb.SolveTime = map(x -> min(upper_limit, x), df_comb.SolveTime)
    df_comb.ShiftedSolveTime = df_comb.SolveTime .+ 1
    df_comb.CorrectlySolved = map(x -> (occursin("OPTIMAL",x) || occursin("LOCALLY_SOLVED",x)), df_comb.TerminationStatus)

    df_comb.IncorrectFeasibility = map(x -> occursin("INFEASIBLE",x), df_comb.TerminationStatus)
    # Presolve infeasibility in SCIP is labelled as OPTIMIZE_NOT_CALLED...
    df_comb.IncorrectFeasibility = map((x,y,z) -> ifelse(occursin("SCIP",x) & occursin("NOT_CALLED",y), true, z), df_comb.SolverName, df_comb.TerminationStatus, df_comb.IncorrectFeasibility)

    df_comb.ActFunc = map(drop_num_underscore, df_comb.InstanceName)

    gdf_comb_cs = groupby(df_comb, Symbol[:FullSolverName, :ActFunc])
    @show combine(gdf_comb_cs, :ShiftedSolveTime => x -> StatsBase.geomean(x) - 1)

    gdf_combt_correct = groupby(df_comb , Symbol[:FullSolverName, :CorrectlySolved, :ActFunc])
    @show combine(gdf_combt_correct, :ShiftedSolveTime => x -> count(fill(true,length(x))))

    df_comb_infeas = df_comb[df_comb.IncorrectFeasibility,:]
    gdf_check_infeasible  = groupby(df_comb_infeas, Symbol[:FullSolverName, :ActFunc])
    @show combine(gdf_check_infeasible, :SolveTime => x -> count(fill(true,length(x))))

    #gdf_status = groupby(df_comb, Symbol[:FullSolverName, :SolvedInTime])
    env_folder = joinpath(@__DIR__, "MINLPLib.jl", "instances", env_lib)
    name_anns = filter(x -> !occursin("gelu",x), readdir(joinpath(env_folder)))

    df_comb.UnsolvedRelativeGap = map((x,y,z) -> (~x & ~y & isfinite(z)), df_comb.CorrectlySolved, df_comb.IncorrectFeasibility, df_comb.RelativeGap)
    df_comb_gap = df_comb[df_comb.UnsolvedRelativeGap,:]
    gdf_comb_gap  = groupby(df_comb_gap, Symbol[:FullSolverName, :ActFunc])
    @show combine(gdf_comb_gap, :RelativeGap => x -> mean(x))

    trunc_solved_time = rand(300, 4)
    for (i,n) in enumerate(name_anns)
        plt_sdf = df_comb[occursin.(n[1:end-3], string.(df_comb.InstanceName)), :]
        trunc_solved_time[i,3] = plt_sdf[plt_sdf.FullSolverName .== "SCIP Exp",:].SolveTime[1]  # SCIP
        trunc_solved_time[i,4] = plt_sdf[plt_sdf.FullSolverName .== "BARON Exp",:].SolveTime[1] # BARON
        trunc_solved_time[i,1] = plt_sdf[plt_sdf.FullSolverName .== "EAGO Env",:].SolveTime[1]  # EAGO Env Entry
        trunc_solved_time[i,2] = plt_sdf[plt_sdf.FullSolverName .== "EAGO Exp",:].SolveTime[1]  # EAGO Exp Entry
    end
    plt = performance_profile(PlotsBackend(), trunc_solved_time, ["EAGO (Envelope)", "EAGO (McCormick)", "SCIP", "BARON"], linewidth = 1.5, linestyles=[:solid, :dash, :dashdot, :dot], legend=:bottomright)
    xlabel!("\$\\tau\$")
    #xlims!(0.0,6.5)
    ylabel!("\$P(r_{p,s} \\leq \\tau : 1 \\leq s \\leq n_s)\$")
    ylims!(0.2,1.0)
    savefig(plt, joinpath(result_path, fig_name*".pdf"))
    show(plt)
end

result_path = joinpath(@__DIR__, "solver_benchmark_result")

params = SolverBenchmarking.BenchmarkParams(time = 100, rerun = false, has_obj_bnd = true)
params_15 = SolverBenchmarking.BenchmarkParams(time = 15*60, rerun = false, has_obj_bnd = true)

SolverBenchmarking.run_solver_benchmark(result_path, expr_solvers, "ANN_Expr", "ANN_Expr"; params = params)
SolverBenchmarking.run_solver_benchmark(result_path, env_solvers, "ANN_Env", "ANN_Env"; params = params)
SolverBenchmarking.run_solver_benchmark(result_path, expr_solvers, "ANN_Expr", "ANN_Expr_15"; params = params_15)
SolverBenchmarking.run_solver_benchmark(result_path, env_solvers, "ANN_Env", "ANN_Env_15"; params = params_15)

SolverBenchmarking.summarize_results("ANN_Expr", result_path)
SolverBenchmarking.summarize_results("ANN_Env", result_path)
SolverBenchmarking.summarize_results("ANN_Expr_15", result_path)
SolverBenchmarking.summarize_results("ANN_Env_15", result_path)

df_expr = DataFrame(CSV.File(joinpath(result_path, "ANN_Expr", "result_summary.csv")))
df_env = DataFrame(CSV.File(joinpath(result_path, "ANN_Env", "result_summary.csv")))
df_expr_15 = DataFrame(CSV.File(joinpath(result_path, "ANN_Expr_15", "result_summary.csv")))
df_env_15 = DataFrame(CSV.File(joinpath(result_path, "ANN_Env_15", "result_summary.csv")))

print_summary_tables(df_env, df_expr, "performance_profile", "ANN_Env", 100.00)
print_summary_tables(df_env_15, df_expr_15, "performance_profile_15", "ANN_Env", 15*60.0)

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function genKnap(num_items::Int64, range::Int64 = 1000; type::Symbol = :Uncorrelated, seed::Int64 = 42, num_tests::Int64 = 1000)::KnapData
    Random.seed!(seed)

    sum = 0
    items = KnapItem[]
    for i in 1:num_items

        weight = rand(1:range)
        sum += weight

        if (type == :Uncorrelated)
            profit = rand(1:range)
        elseif (type == :WeakCorrelated)
            profit = rand(1:(range ÷ 5 + 1)) + weight - range ÷ 10
            if (profit <= 0) profit = 1 end
        elseif (type == :StrongCorrelated)
            profit = weight + 10
        elseif (type == :SubsetSum)
            profit = weight
        else
            error("Wrong profit correlation type! Use :Uncorrelated, :WeakCorrelated, :StrongCorrelated or :SubsetSum.")
        end

        push!(items, KnapItem(weight, profit))
    end

    capacity = (seed * sum) ÷ (num_tests + 1)
    if (capacity <= range) capacity = range + 1 end

    return KnapData(capacity, items)
end


text:
# Autogenerated wrapper script for lttng_ust_jll for aarch64-linux-gnu
export liblttng_ust, liblttng_ust_ctl, liblttng_ust_cyg_profile, liblttng_ust_cyg_profile_fast, liblttng_ust_dl, liblttng_ust_fd, liblttng_ust_fork, liblttng_ust_libc_wrapper, liblttng_ust_pthread_wrapper, liblttng_ust_tracepoint

using NUMA_jll
using URCU_jll
JLLWrappers.@generate_wrapper_header("lttng_ust")
JLLWrappers.@declare_library_product(liblttng_ust, "liblttng-ust.so.0")
JLLWrappers.@declare_library_product(liblttng_ust_ctl, "liblttng-ust-ctl.so.4")
JLLWrappers.@declare_library_product(liblttng_ust_cyg_profile, "liblttng-ust-cyg-profile.so.0")
JLLWrappers.@declare_library_product(liblttng_ust_cyg_profile_fast, "liblttng-ust-cyg-profile-fast.so.0")
JLLWrappers.@declare_library_product(liblttng_ust_dl, "liblttng-ust-dl.so.0")
JLLWrappers.@declare_library_product(liblttng_ust_fd, "liblttng-ust-fd.so.0")
JLLWrappers.@declare_library_product(liblttng_ust_fork, "liblttng-ust-fork.so.0")
JLLWrappers.@declare_library_product(liblttng_ust_libc_wrapper, "liblttng-ust-libc-wrapper.so.0")
JLLWrappers.@declare_library_product(liblttng_ust_pthread_wrapper, "liblttng-ust-pthread-wrapper.so.0")
JLLWrappers.@declare_library_product(liblttng_ust_tracepoint, "liblttng-ust-tracepoint.so.0")
function __init__()
    JLLWrappers.@generate_init_header(NUMA_jll, URCU_jll)
    JLLWrappers.@init_library_product(
        liblttng_ust,
        "lib/liblttng-ust.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblttng_ust_ctl,
        "lib/liblttng-ust-ctl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblttng_ust_cyg_profile,
        "lib/liblttng-ust-cyg-profile.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblttng_ust_cyg_profile_fast,
        "lib/liblttng-ust-cyg-profile-fast.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblttng_ust_dl,
        "lib/liblttng-ust-dl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblttng_ust_fd,
        "lib/liblttng-ust-fd.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblttng_ust_fork,
        "lib/liblttng-ust-fork.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblttng_ust_libc_wrapper,
        "lib/liblttng-ust-libc-wrapper.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblttng_ust_pthread_wrapper,
        "lib/liblttng-ust-pthread-wrapper.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblttng_ust_tracepoint,
        "lib/liblttng-ust-tracepoint.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Hygiene

using Test

import MathOptInterface
const MOI = MathOptInterface

# Dict is used in the @model macro but setting Dict in the outer scope
# should not affect it
Dict = nothing

MathOptInterface.Utilities.@model(
    LPModel,                      # Name of model
    (),                                                         # untyped scalar sets
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval), #   typed scalar sets
    (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives),            # untyped vector sets
    (),                                                         #   typed vector sets
    (),                                                         # untyped scalar functions
    (MOI.ScalarAffineFunction,),                                #   typed scalar functions
    (MOI.VectorOfVariables,),                                   # untyped vector functions
    (MOI.VectorAffineFunction,)
)                                #   typed vector functions

model = LPModel{Float64}()

@test model isa MathOptInterface.ModelLike
@test model isa MathOptInterface.Utilities.AbstractModel{Float64}

end


text:
using Base.Iterators
using DataStructures
using StaticArrays

struct State{N}
    value::SVector{N, Int8}
end

struct Movement{N}
    value::SVector{N, Int8}
end

State(x...) = State(SVector{length(x), Int8}(x...))
Movement(x...) = Movement(SVector{length(x), Int8}(x...))

function Base.ones(::Type{State{N}})::State{N} where {N}
    return State(ones(SVector{N, Int8}))
end

function Base.zeros(::Type{State{N}})::State{N} where {N}
    return State(zeros(SVector{N, Int8}))
end

function is_valid(state::State)::Bool
    return all(0 .<= state.value .<= 1)
end

Base.show(io::IO, state::State) = for el in state.value print(io, el) end
function Base.show(io::IO, movement::Movement)
    for el in movement.value
        char = if el == -1
            '-'
        elseif el == 0
            '.'
        elseif el == 1
            '+'
        end
        print(io, char)
    end
end

function Base.:(+)(state::State{N}, movement::Movement{N})::State{N} where {N}
    return State(state.value + movement.value)
end

"""
Get all movements in increasing order of applicability.

That is, movements near the start of the list contain fewer 'no change' actions, and as such
are usable in fewer situations.
"""
function get_all_movements(n::Integer)
    actions = (-1, 0, 1)
    all_movements = collect(product((actions for _ in 1:n)...))
    all_movements = [Movement(movement...) for movement in @view all_movements[:]]
    sort!(all_movements; by=x -> sum(abs.(x.value)), rev=true)
    result = MutableLinkedList()
    for movement in all_movements
        push!(result, movement)
    end
    return result
end

function get_movement_cycle(initial::State{N})::Vector{Movement{N}} where {N}
    movement_cycle = Movement{N}[]
    all_movements = get_all_movements(N)
    state = initial
    while !isempty(all_movements)
        found = false
        for (i, movement) in enumerate(all_movements)
            new_state = state + movement
            if is_valid(new_state)
                delete!(all_movements, i)
                push!(movement_cycle, movement)
                state = new_state
                found = true
                break
            end
        end
        found || error("Was not able to find lower bound cycle!")
    end
    state == initial || error("Did not end in the initial state!")
    return movement_cycle
end

print_movements(state, cycle) = print_movements(stdout, state, cycle)

function print_movements(io::IO, state::State{N}, cycle::Vector{Movement{N}}) where {N}
    for movement in cycle
        println(io, "$state    $movement")
        state += movement
        is_valid(state) || error("Invalid movement $movement from $state")
    end
    println(io, state)
    return nothing
end

for n in 1:7
    initial = zeros(State{n})
    movement_cycle = get_movement_cycle(initial)
    open(joinpath("output", "movements_$n.txt"); write=true) do io
        print_movements(io, initial, movement_cycle)
    end
end


text:
using DataDeps

register(DataDep(
  "connectionist-bench",
  "https://archive.ics.uci.edu/ml/datasets/Connectionist+Bench+(Vowel+Recognition+-+Deterding+Data)",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/undocumented/connectionist-bench/vowel/vowel-context.data",
  "1e5c40bdd0b6705c75586eb23f5d9c3078e9588f48bc5c2c98da219fe83fff29",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=14,
        feature_indices=4:13,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
#=
   Classical Random Walk 1

   From
   Gunnar Blom, Lars Holst, Dennis Sandell:
   "Problems and Snapshots from the World of Probability"
   Page 6f, Problem 1.5 Classical Random Walk 1,
   cases a) Probability of absorption
         b) Expected number of steps until absorption
         c) Ruin problem

   The probability of moving to the left or right is 0.5 (symmetric
   random walk). The walk stops when reaching either -a or +b (
   a and are both positive integers).

   * a) and b)  a:10 and b:10
     a:10 b:10
     Theoretical prob of [-a,b]:[0.5, 0.5]
     Theoretical length of absorption:100

     Summary Statistics
      parameters       mean       std   naive_se      mcse          ess      rhat   ess_per_sec
           Symbol    Float64   Float64    Float64   Float64      Float64   Float64       Float64

             len   100.6312   82.2539     0.8225    0.7312   10156.6554    0.9999     8164.5140
         lastPos     0.0460   10.0004     0.1000    0.0915    9771.6168    1.0000     7854.9975

     Distributions of variable lastPos (num:0)
     10.00000 =>    5023  (0.502300)
     -10.00000 =>    4977  (0.497700)

   * c) The ruin problem, a=1 and b=100:

     a:1 b:100
     Theoretical prob of [-a,b]:[0.9900990099009901, 0.009900990099009901]
     Theoretical length of absorption:100

     Summary Statistics
     parameters       mean        std   naive_se      mcse          ess      rhat   ess_per_sec
         Symbol    Float64    Float64    Float64   Float64      Float64   Float64       Float64

            len   110.7316   612.2311     6.1223    6.2518    9854.2103    0.9999    62765.6706
        lastPos     0.0504    10.2468     0.1025    0.1034   10218.6345    0.9999    65086.8437

    Distributions of variable lastPos (num:0)
    -1.00000 =>    9896  (0.989600)
    100.00000 =>     104  (0.010400)

  Cf ~/webppl/random_walk_1.wppl

=#
using Turing, StatsPlots
include("jl_utils.jl")

random_walk_1_theoretical_prob = (a,b) -> [b/(a+b),a/(a+b)]
random_walk_1_theoretical_length = (a,b) -> a*b

@model function random_walk_1(a=10,b=10)
    t = [-1,1]
    function walk(arr)
        len = length(arr)
        last = len == 0 ? 0 : arr[end]
        if last == -a || last == b
            return arr
        else
            # d ~ UniformDraw([-1,1])
            d = rand(t)
            return walk(push!(arr,last + d) )
        end
    end

    arr = walk([])
    len ~ Dirac(length(arr))
    lastPos ~ Dirac(arr[end])

    return arr
end

function run_random_walk(a,b)
    println("a:$a b:$b")
    println("Theoretical prob of [-a,b]:", random_walk_1_theoretical_prob(a,b))
    println("Theoretical length of absorption:", random_walk_1_theoretical_length(a,b))
    model = random_walk_1(a,b)

    chns = sample(model, Prior(), 10_000)
    # chns = sample(model, MH(), 10_000)
    # chns = sample(model, PG(5), 10_000)
    # chns = sample(model, SMC(), 10_000)
    # chns = sample(model, IS(), 10_000)

    display(chns)
    # display(plot(chns))
    show_var_dist_pct(chns, :lastPos)

    # chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
    # genq = generated_quantities(model, chains_params)
    # show_var_dist_pct(genq)
    println()
end

run_random_walk(10,10)
run_random_walk(1,100)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module SameRand
	using FromFile
    @from "../types.jl" import MyInt
    @from "utils.jl" import add_same_rand

    add_same_rand1(x::MyInt)::MyInt = add_same_rand(x)
end


text:
using BinaryBuilder

name = "HelloWorldRust"
version = v"1.0.2"

# No sources, we're just building the testsuite
sources = [
]

# Bash recipe for building across all platforms
script = raw"""
mkdir -p ${prefix}/bin
rustc -o ${prefix}/bin/hello_world${exeext} -g /usr/share/testsuite/rust/hello_world/hello_world.rs
"""

# We build for a restricted set of platforms, because our rust toolchain is a little broken
platforms = supported_platforms()

# First, FreeBSD has -fPIC problems when linking in `crt.o`
filter!(!Sys.isfreebsd, platforms)

# Next, musl libcs have a hard time linking
filter!(p -> libc(p) != "musl", platforms)

# Finally, windows seems to be broken
# https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/499
filter!(!Sys.iswindows, platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("hello_world", :hello_world),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers=[:c, :rust])

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#reproducing part of our results is as easy as
#running this script from start to end in Julia (>=1.6)
using Pkg
Pkg.add("SparseArrays")
Pkg.add("Roots")
Pkg.add(name="Lasso", version="0.6.2")
Pkg.add("LogDensityProblems")
Pkg.add("Parameters")
Pkg.add("Random")
Pkg.add("DynamicHMC")
Pkg.add("DataFrames")
Pkg.add("MCMCChains")
Pkg.add("Statistics")
Pkg.add("StatsPlots")

using SparseArrays: SparseMatrixCSC,spdiagm
using Roots: find_zero, Bisection
using Lasso: fit, FusedLasso
import LogDensityProblems.dimension,LogDensityProblems.capabilities
import LogDensityProblems.logdensity_and_gradient,LogDensityProblems.LogDensityOrder
using Parameters: @unpack
using Random: GLOBAL_RNG
using DynamicHMC: mcmc_with_warmup, default_warmup_stages, ProgressMeterReport, NoProgressReport
using DynamicHMC.Diagnostics: summarize_tree_statistics
using DataFrames: DataFrame
using MCMCChains: Chains,summarize,quantile
using Statistics: mean,norm,var
using StatsPlots: plot, scatter!, plot!

function getD(k::Int,n::Int,x::AbstractArray{T}=[1.:1.:n;]) where T<:Real
    D = spdiagm(n-1,n,0=>fill(-1.,n-1),1=>fill(1.,n-1))
    if k>= 1
        for i in 1:k
            D1 = spdiagm(n-i-1,n-i,0=>fill(-1.,n-i-1),1=>fill(1.,n-i-1))
            diffx = x[(i+1):n]-x[1:(n-i)]
            xdiag = spdiagm(0=> i ./ diffx)
            D = D1*xdiag*D
        end
    end
    D
end

function getTM(k::Int,n::Int,x::Vector{T}=[1.:1.:n;]) where T<:Real
    if (k==1) & (n<=200)
        D = getD(k,n,x)
        return [spdiagm(k+1,n,0=>fill(1.,k+1));D]
    else
        D = getD(k-1,n,x)
        diffx = x[(k+1):n]-x[1:(n-k)]
        xdiag = spdiagm(0=> k ./ diffx)
        return [spdiagm(k,n,0=> fill(1.,k));xdiag*D]
    end
end

function forward_solve!(β::AbstractArray{T},bw::Int,TM::SparseMatrixCSC{T,Int},θ::AbstractArray{T}) where T<:Real
    n = length(θ)
    copyto!(β,θ)
    @inbounds for i in 1:n
        for j in 1:min(i-1,bw-1)
            β[i] -= TM[i,i-j]*β[i-j]
        end
        β[i] /= TM[i,i]
    end
    return zero(T)
end

function back_solve!(∇β::AbstractArray{T},bw::Int,TMᵀ::SparseMatrixCSC{T,Int},x::AbstractArray{T}) where T<:Real
    n = length(x)
    copyto!(∇β,x)
    @inbounds for i in n:-1:1
        for j in 1:min(n-i,bw-1)
            ∇β[i] -= TMᵀ[i,i+j]*∇β[i+j]
        end
        ∇β[i] /= TMᵀ[i,i]
    end
    return zero(T)
end

function thin(xgrid::AbstractArray{T},w::AbstractArray{T},ybar::AbstractArray{T},sse::T,nbins::Int) where T<:Real
    xmin = minimum(xgrid)
    range = maximum(xgrid)-xmin
    gap = range/nbins
    w_new = zeros(T,nbins)
    xgrid_new = zeros(T,nbins)
    ybar_new = zeros(T,nbins)
    sse_new = 0.0
    for i in 1:nbins
        subset = nothing
        if i==1
            subset = xmin .<= xgrid .<= xmin+gap
        else
            subset = xmin+gap*(i-1) .< xgrid .<= xmin+i*gap
        end
        w_new[i] = sum(w[subset])
        if w_new[i]>0
            xgrid_new[i] = sum(w[subset].*xgrid[subset])/w_new[i]
            ybar_new[i] = sum(w[subset].*ybar[subset])/w_new[i]
            sse_new += sum(w[subset].* abs2.(ybar[subset].-ybar_new[i]))
        else
            xgrid_new[i] = xmin+gap*(i-0.5)
            ybar_new[i] = mean(y)
        end
    end
    sse_new += sse
    return xgrid_new,w_new,ybar_new,sse_new
end

function F₀(x::AbstractArray{T}, α::T, λ::T) where T<:Real
    out = zero(T)
    @inbounds for i in 1:length(x)
        δᵢ = max(abs(x[i])-λ,0)
        out += δᵢ
    end
    out -= α+λ
    out
end

function F₁(x::AbstractArray{T}, α::T, λ::T) where T<:Real
    if λ==0
        return sum(abs,diff(x))
    else
        x_prox = fit(FusedLasso,x,λ).β
        return sum(abs,diff(x_prox))-α-λ
    end
end

function absepi_projection!(θ_proj::AbstractArray{T}, skip::Int, θ::AbstractArray{T}, α::T) where T<:Real
    copyto!(θ_proj,θ)
    n = length(θ)
    @views if sum(abs, θ[(skip+1):n]) <= α
        return α
    else
        λstar = find_zero(λ -> F₀(θ[(skip+1):n],α,λ),(0,maximum(abs,θ[(skip+1):n])+max(0,-α)),Bisection())
        @inbounds for i in 1:(n-skip)
            θ_proj[skip+i] = sign(θ[skip+i])*max(abs(θ[skip+i])-λstar,0)
        end
        return α + λstar
    end
end

function diffepi_projection!(θ_proj::AbstractArray{T}, skip::Int, θ::AbstractArray{T}, α::T, inv_DDᵀ_D::AbstractArray{T}) where T<:Real
    copyto!(θ_proj,θ)
    n = length(θ)
    if sum(abs, diff(θ[(skip+1):n])) <= α
        return α
    else
        λₘ = maximum(abs, inv_DDᵀ_D*θ[(skip+1):n])
        λstar = find_zero(λ -> F₁(θ[(skip+1):n],α,λ),(0,λₘ+max(0,-α)),Bisection())
        θ_proj[(skip+1):n] = fit(FusedLasso,θ[(skip+1):n],λstar).β
        return α + λstar
    end
end

mutable struct pbtfProblem{T<:Real}
    y       :: Vector{T}               # response
    x       :: Vector{T}               # predictor
    k       :: Int                     # order of desired polynomial, only 1 or 2 is allowed
    λ       :: Union{T,Nothing}        # Moreau-Yosida envelope parameter
    s       :: T                       # IG(s, r) prior for σ²
    r       :: T
    s₂      :: Union{T,Nothing}        # β'(n-k,s₂) prior for α
    m       :: Int                     # total number of observations
    n       :: Int                     # number of grid locations
    xgrid   :: Vector{T}               # strictly increasing grid locations (length n)
    w       :: Vector{T}               # number of observations at each location (length n)
    ybar    :: Vector{T}               # the average of y at each location (length n)
    sse     :: T                       # sum of squared errors
    TM      :: SparseMatrixCSC{T,Int}  # n×n transformation matrix
    TMᵀ     :: SparseMatrixCSC{T,Int}  # n×n transformation matrix transposed
    D       :: SparseMatrixCSC{T,Int}  # (n-k-1)×n difference matrix
    inv_DDᵀ_D :: Matrix{T}             # used to determine the bisection interval
    θ_proj  :: Vector{T}               # the following arrays are pre-created containers to avoid creating
    β       :: Vector{T}               # new arrays during the course of MCMC sampling
    res     :: Vector{T}
end

function pbtfProblem(y::AbstractArray{T}, x::AbstractArray{T}, k::Int;
                     thinning::Bool=false,nbins::Int=100,λ::Union{T,Nothing}=nothing,
                     s::T=1e-3,r::T=1e-3,s₂::Union{T,Nothing}=nothing) where T<: Real
    if (k!=1) & (k!=2)
        throw(ArgumentError("Currenly only k=1 and k=2 are accepted for pbtf."))
    end
    m = length(x)
    xgrid = sort(unique(x))
    n = length(xgrid)
    w = zeros(T,n)
    ybar = zeros(T,n)
    sse = 0.0
    for i in 1:n
        w[i] = length(x[x.==xgrid[i]])
        ybar[i] = mean(y[x.==xgrid[i]])
        sse += sum(abs2, y[x.==xgrid[i]].-ybar[i])
    end
    if thinning
        xgrid,w,ybar,sse = thin(xgrid,w,ybar,sse,nbins)
        n = nbins
    end
    if λ==nothing
        λ = min(1/n^2,1e-4*var(y))
    end
    if s₂==nothing
        s₂ = sqrt(n)
    end
    TM = getTM(k,n,xgrid)
    TMᵀ = SparseMatrixCSC(TM')
    D1 = getD(0,n-k)
    inv_DDᵀ_D = inv(Matrix(D1*D1'))*D1
    D = getD(k,n,xgrid)
    θ_proj = zeros(T,n)
    β = zeros(T,n)
    res = zeros(T,n)
    pbtfProblem{T}(y,x,k,λ,s,r,s₂,m,n,xgrid,w,ybar,sse,TM,TMᵀ,D,inv_DDᵀ_D,θ_proj,β,res)
end

dimension(problem::pbtfProblem) = problem.n + 2
capabilities(::Type{<:pbtfProblem}) = LogDensityOrder{1}()

function logdensity_and_gradient(problem::pbtfProblem, z)
    @unpack y,x,k,λ,s,r,s₂,m,n,xgrid,w,ybar,sse,TM,TMᵀ,D,inv_DDᵀ_D,θ_proj,β,res=problem
    θ         = view(z,1:n)
    logσ²     = z[n+1]
    σ²        = exp(logσ²)
    logα      = z[n+2]
    α         = exp(logα)
    α_proj    = 0.0
    if (k==1) & (n<=200)
        α_proj = absepi_projection!(θ_proj,k+1,θ,α)
        forward_solve!(β,k+2,TM,θ)
    else
        α_proj = diffepi_projection!(θ_proj,k,θ,α,inv_DDᵀ_D)
        forward_solve!(β,k+1,TM,θ)
    end
    res      .= ybar .- β
    ssm       = sum(res .* w .* res)
    qf        = (ssm+sse+2*r)/(2*σ²)
    dist²     = norm(θ.-θ_proj,2)^2+(α-α_proj)^2
    logl      = -dist²/(2*λ)-qf-(m/2+s)*logσ²
    ∇         = similar(z)
    if (k==1) & (n<=200)
        back_solve!(view(∇,1:n),k+2,TMᵀ,w.*res)
    else
        back_solve!(view(∇,1:n),k+1,TMᵀ,w.*res)
    end
    logl += logα-(n-k+s₂)*log(1+α)
    ∇[n+2]  = 1 -α/λ*(α-α_proj)-(n-k+s₂)*α/(1+α)
    ∇[1:n]  ./= σ²
    ∇[1:n]  .-= (θ.-θ_proj)./λ
    ∇[n+1]    = qf-(m/2+s)
    logl, ∇
end

function pbtf(
    y::Vector{T},                 # response
    x::Vector{T},                 # predictor
    k::Int;                       # order of desired polynomial, only 1 or 2 is acceptable
    thinning::Bool=false,         # whether to apply the thinning technique
    nbins::Int=100,               # number of bins if thinning is applied
    λ::Union{T,Nothing}=nothing,  # the Moreau envelope parameter
    s::T=1e-3,                    # IG(s,r) prior for σ²
    r::T=1e-3,
    s₂::Union{T,Nothing}=nothing, # β'(n-k,s₂) prior for α
    nsample::Int=3000,            # number of MCMC sample iterations
    nthin::Int=5,                 # thinning the chain after MCMC sampling
    verbose::Bool=true            # whether to display progress bar
    ) where T<: Real
    problem = pbtfProblem(y,x,k,thinning=thinning,nbins=nbins,λ=λ,s=s,r=r,s₂=s₂)
    X = [ones(T,problem.n) problem.xgrid]
    q₀ = [problem.TM*X*(X\problem.ybar);0;0]
    if verbose
        elapse_time = @elapsed results = mcmc_with_warmup(GLOBAL_RNG, problem, nsample, initialization=(ϵ=0.01,q=q₀),
        warmup_stages=default_warmup_stages(;stepsize_search=nothing),reporter=ProgressMeterReport())
    else
        elapse_time = @elapsed results = mcmc_with_warmup(GLOBAL_RNG, problem, nsample, initialization=(ϵ=0.01,q=q₀),
        warmup_stages=default_warmup_stages(;stepsize_search=nothing),reporter=NoProgressReport())
    end
    println(summarize_tree_statistics(results.tree_statistics))
    results_matrix = reshape([(results.chain...)...], length(results.chain[1]), length(results.chain))
    results_matrix[1:problem.n,:] = problem.TM \ results_matrix[1:problem.n,:]
    results_matrix[end-1:end,:] = exp.(results_matrix[end-1:end,:])
    chn = Chains(transpose(results_matrix),["β[" .* string.(1:problem.n) .*"]"; "σ²";"α"],thin=nthin)
    result_summary = DataFrame(summarize(chn))
    result_quantile = DataFrame(quantile(chn))
    return result_summary,result_quantile,problem,elapse_time
end

function visualize(result_quantile::DataFrame,problem::pbtfProblem;legend_position::Symbol=:topleft)
    plt=plot(size=(800,600))
    scatter!(problem.x,problem.y,markersize=5,markercolor=:green,label="data",legendfontsize=20,
             xtickfontsize=20,ytickfontsize=20,markerstrokewidth=0,legend=legend_position)
    plot!(subplot=1,problem.xgrid,result_quantile[1:problem.n, Symbol("50.0%")],linewidth=5,label="posterior",c=:blue,
          ribbon=(result_quantile[1:problem.n, Symbol("50.0%")] .- result_quantile[1:problem.n, Symbol("2.5%")],
          result_quantile[1:problem.n, Symbol("97.5%")] .- result_quantile[1:problem.n, Symbol("50.0%")]),fillalpha=0.3)
end

#x = [1:1.:100;]
#fx = map(xi -> 0 ≤ xi ≤ 35 ? xi : 35 < xi ≤ 70 ? 70-xi : 0.5xi-35, x)
#σ = 3.0
#y = fx .+  σ.*randn(100)
#result_summary,result_quantile,problem,elapse_time = pbtf(y,x,1)
#visualize(result_quantile,problem)

#x = sort(100*rand(1000))
#σ = 3.0
#fx  = 13 .* sin.((4*π/100).*x)
#y = fx .+  σ.*randn(1000)
#result_summary,result_quantile,problem,elapse_time=pbtf(y,x,2;thinning=true)
#visualize(result_quantile,problem)


text:
function insertSort(lst)
    for i in 2:length(lst)
        j = i - 1
        key = lst[i]
        while j >= 1 && key > lst[j]
            lst[j+1] = lst[j]
            j -= 1
        end
        lst[j+1] = key
    end
end

lst = [1, 2, 4, 5, 6, 8, 9, 11, 13, 14, 19, 22]
println(lst)
insertSort(lst)
println(lst)

max_stars_repo_name: Bingwen-Hu/hackaway

max_issues_repo_name: Bingwen-Hu/hackaway

max_forks_repo_name: Bingwen-Hu/hackaway


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the project(x, y) function
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using LinearAlgebra: norm, ⋅, UniformScaling, I
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- File inclusions

# Test utilities
include("test_utils.jl")

# --- Tests

# ------ M::Multivector

@testset "project(M::Multivector, N::Multivector)" begin
    @test_skip 1
end

@testset "project(M::Multivector, N::Blade)" begin
    @test_skip 1
end

@testset "project(M::Multivector, N::Pseudoscalar)" begin
    @test_skip 1
end

@testset "project(M::Multivector, N::Scalar)" begin
    @test_skip 1
end

@testset "project(M::Multivector, N::One)" begin
    @test_skip 1
end

@testset "project(M::Multivector, N::Zero)" begin
    @test_skip 1
end

@testset "project(M::Multivector, N::Real)" begin
    @test_skip 1
end

@testset "project(M::Multivector, N::Vector)" begin
    @test_skip 1
end

# ------ B::Blade

@testset "project(B::Blade, M::Multivector)" begin
    @test_skip 1
end

@testset "project(B::Blade, C::Blade)" begin
    # --- Preparations

    test_dim = 10

    # --- grade(B) < grade(C)

    B = Blade(rand(test_dim, 5))
    C = Blade(rand(test_dim, 7))

    # ------ return_blade == true

    B_proj_C = project(B, C)
    @test B_proj_C isa Blade
    @test dim(B_proj_C) == dim(B)

    # Check norm
    expected_norm_B_proj_C = norm(B) *
        norm(Blade(basis(C) * transpose(basis(C)) * basis(B)))
    @test norm(B_proj_C) ≈ expected_norm_B_proj_C

    # Check that project(B, C) is contained in C
    projection_coefficients = transpose(basis(C)) * basis(B_proj_C)
    @test norm(projection_coefficients)^2 ≈ grade(B)

    # ------ return_blade == false

    B_proj_C_matrix = project(B, C, return_blade=false)
    @test B_proj_C_matrix isa Matrix
    @test Blade(B_proj_C_matrix) == B_proj_C

    # --- grade(B) > grade(C)

    C = Blade(rand(test_dim, 3))

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C_matrix = project(B, C, return_blade=false)
    @test B_proj_C_matrix isa Real
    @test B_proj_C_matrix == 0

    # --- dim(B) != dim(C)

    C = Blade(rand(test_dim + 1, 4))
    @test_throws DimensionMismatch project(B, C, return_blade=true)
    @test_throws DimensionMismatch project(B, C, return_blade=false)

    # --- Check consistency with project(v::Vector, B::Blade)

    v = Vector(rand(test_dim))
    @test project(v, B) ≈ project(Blade(v), B)
end

@testset "project(B::Blade, C::Pseudoscalar)" begin
    # --- Preparations

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value

    test_dim = 5
    test_vectors = rand(test_dim, 3)

    # --- dim(B) == dim(C)

    B = Blade(test_vectors)
    C = Pseudoscalar(test_dim, test_value)

    # return_blade == true
    @test project(B, C) === B

    # return_blade == false
    B_proj_C_matrix = project(B, C, return_blade=false)
    @test B_proj_C_matrix isa Matrix
    @test Blade(B_proj_C_matrix) ≈ B

    # --- dim(B) != dim(C)

    B = Pseudoscalar(test_dim + 1, test_value)
    C = Blade(test_vectors)
    @test_throws DimensionMismatch project(B, C)
    @test_throws DimensionMismatch project(B, C, return_blade=false)
end

@testset "project(B::Blade, C::Scalar)" begin
    # Preparations
    test_dim = 5
    B = Blade(rand(test_dim, 3))

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Scalar(test_value)

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C_matrix = project(B, C, return_blade=false)
    @test B_proj_C_matrix isa Real
    @test B_proj_C_matrix == 0
end

@testset "project(B::Blade, C::One)" begin
    # Preparations
    test_dim = 5
    B = Blade(rand(test_dim, 3))
    C = One()

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C_matrix = project(B, C, return_blade=false)
    @test B_proj_C_matrix isa Real
    @test B_proj_C_matrix == 0
end

@testset "project(B::Blade, C::Zero)" begin
    # Preparations
    test_dim = 5
    B = Blade(rand(test_dim, 3))
    C = Zero()

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C_matrix = project(B, C, return_blade=false)
    @test B_proj_C_matrix isa Real
    @test B_proj_C_matrix == 0
end

@testset "project(B::Blade, C::Real)" begin
    # Preparations
    test_dim = 5
    B = Blade(rand(test_dim, 3))

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = test_value

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C_matrix = project(B, C, return_blade=false)
    @test B_proj_C_matrix isa Real
    @test B_proj_C_matrix == 0
end

@testset "project(B::Blade, C::Vector)" begin
    # --- Preparations

    test_dim = 15
    C = rand(test_dim)

    # --- Tests

    # grade(B) > 1, return_blade == true
    B = Blade(rand(test_dim, 5))

    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # grade(B) > 1, return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0

    # grade(B) == 1, return_blade == true
    B = Blade(rand(test_dim, 1))
    projection_vectors = (C ⋅ basis(B)) * C
    @test project(B, C) ≈ Blade(projection_vectors)

    # grade(B) == 1, return_blade == false
    @test project(B, C, return_blade=false) ≈ projection_vectors

    # dim(B) != length(C)
    B = Blade(rand(test_dim + 1, 3))
    @test_throws DimensionMismatch project(B, C)
    @test_throws DimensionMismatch project(B, C, return_blade=false)
end

# ------ B::Pseudoscalar

@testset "project(B::Pseudoscalar, M::Multivector)" begin
    @test_skip 1
end

@testset "project(B::Pseudoscalar, C::Blade)" begin
    # --- Preparations

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value

    test_dim = 10
    test_vectors = rand(test_dim, 3)

    # --- Tests

    # ------ dim(B) == dim(C)

    B = Pseudoscalar(test_dim, test_value)
    C = Blade(test_vectors)

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0

    # ------ dim(B) != dim(C)

    B = Blade(test_vectors)
    C = Pseudoscalar(test_dim + 1, test_value)
    @test_throws DimensionMismatch project(B, C)
    @test_throws DimensionMismatch project(B, C, return_blade=false)
end

@testset "project(B::Pseudoscalar, C::Pseudoscalar)" begin
    # --- Preparations

    test_value_1 = get_random_value(1)  # add 1 to keep value away from 0
    test_value_2 = get_random_value(1)  # add 1 to keep value away from 0

    # --- Tests

    for test_dim in 5:8
        # --- dim(B) == dim(C)

        B = Pseudoscalar(test_dim, test_value_1)
        C = Pseudoscalar(test_dim, test_value_2)

        # return_blade == true
        B_proj_C = project(B, C)
        @test B_proj_C === B

        # return_blade == false
        B_proj_C = project(B, C, return_blade=false)
        @test B_proj_C isa UniformScaling
        @test B_proj_C == value(B) * I

        # --- dim(B) != dim(C)

        B = Pseudoscalar(test_dim, test_value_1)
        C = Pseudoscalar(test_dim + 1, test_value_2)
        @test_throws DimensionMismatch project(B, C)
        @test_throws DimensionMismatch project(B, C, return_blade=false)
    end
end

@testset "project(B::Pseudoscalar, C::Scalar)" begin
    test_dim = 10
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Pseudoscalar(test_dim, test_value_1)

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Scalar(test_value_2)

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Pseudoscalar, C::One)" begin
    test_dim = 10
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Pseudoscalar(test_dim, test_value)

    C = One()

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Pseudoscalar, C::Zero)" begin
    test_dim = 10
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Pseudoscalar(test_dim, test_value)

    C = Zero()

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Pseudoscalar, C::Real)" begin
    test_dim = 10
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Pseudoscalar(test_dim, test_value_1)

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = test_value_2

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Pseudoscalar, C::Vector)" begin
    # --- Preparations

    test_dim = 10

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Pseudoscalar(test_dim, test_value)

    C = Vector(rand(test_dim))

    # --- Tests

    # ------ dim(B) == length(C)

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0

    # ------ dim(B) != length(C)

    B = Pseudoscalar(test_dim + 1, test_value)
    @test_throws DimensionMismatch project(B, C)
    @test_throws DimensionMismatch project(B, C, return_blade=false)
end

# ------ B::Scalar

@testset "project(B::Scalar, M::Multivector)" begin
    @test_skip 1
end

@testset "project(B::Scalar, C::Blade)" begin
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Scalar(test_value)

    test_dim = 5
    C = Blade(rand(test_dim, 3))

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C isa AbstractScalar
    @test B_proj_C == test_value

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == test_value
end

@testset "project(B::Scalar, C::Pseudoscalar)" begin
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Scalar(test_value_1)

    test_dim = 10
    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Pseudoscalar(test_dim, test_value_2)

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C isa AbstractScalar
    @test B_proj_C == test_value_1

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == test_value_1
end

@testset "project(B::Scalar, C::Scalar)" begin
    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value_1)

    test_value_2 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value_2)

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C === B

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == test_value_1
end

@testset "project(B::Scalar, C::One)" begin
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value)

    C = One()

    # return_blade == true
    @test project(B, C) === B

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == test_value
end

@testset "project(B::Scalar, C::Zero)" begin
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value)

    C = Zero()

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Scalar, C::Real)" begin
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Scalar(test_value_1)

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = test_value_2

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C isa AbstractScalar
    @test B_proj_C == test_value_1

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == test_value_1
end

@testset "project(B::Scalar, C::Vector)" begin
    # Preparations
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Scalar(test_value)

    test_dim = 15
    C = rand(test_dim)

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C isa AbstractScalar
    @test B_proj_C == test_value

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == test_value
end

# ------ B::One

@testset "project(B::One, M::Multivector)" begin
    @test_skip 1
end

@testset "project(B::One, C::Blade)" begin
    B = One()

    test_dim = 5
    C = Blade(rand(test_dim, 3))

    # return_blade == true
    B_proj_C = project(B, C)
    @test isone(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 1
end

@testset "project(B::One, C::Pseudoscalar)" begin
    B = One()

    test_dim = 10
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Pseudoscalar(test_dim, test_value)

    # return_blade == true
    B_proj_C = project(B, C)
    @test isone(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 1
end

@testset "project(B::One, C::Scalar)" begin
    B = One()

    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value)

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C === B

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 1
end

@testset "project(B::One, C::One)" begin
    B = One()
    C = One()

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C === B

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 1
end

@testset "project(B::One, C::Zero)" begin
    B = One()
    C = Zero()

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::One, C::Real)" begin
    B = One()

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = test_value

    # return_blade == treu
    B_proj_C = project(B, C)
    @test isone(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 1
end

@testset "project(B::One, C::Vector)" begin
    # return_blade == true
    B = One()
    C = Vector(rand(10))

    # return_blade == true
    B_proj_C = project(B, C)
    @test isone(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 1
end

# ------ B::Zero

@testset "project(B::Zero, M::Multivector)" begin
    @test_skip 1
end

@testset "project(B::Zero, C::Blade)" begin
    B = Zero()

    test_dim = 5
    C = Blade(rand(test_dim, 3))

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Zero, C::Pseudoscalar)" begin
    B = Zero()

    test_dim = 10
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Pseudoscalar(test_dim, test_value)

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Zero, C::Scalar)" begin
    B = Zero()

    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value)

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C === B

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Zero, C::One)" begin
    B = Zero()
    C = One()

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C === B

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Zero, C::Zero)" begin
    B = Zero()
    C = Zero()

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C === B

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Zero, C::Real)" begin
    B = Zero()

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = test_value

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Zero, C::Vector)" begin
    B = Zero()
    C = Vector(rand(10))

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

# ------ B::Real

@testset "project(B::Real, M::Multivector)" begin
    @test_skip 1
end

@testset "project(B::Real, C::Blade)" begin
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = test_value

    test_dim = 5
    C = Blade(rand(test_dim, 3))

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C isa AbstractScalar
    @test B_proj_C == test_value

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == test_value
end

@testset "project(B::Real, C::Pseudoscalar)" begin
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = test_value_1

    test_dim = 10
    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Pseudoscalar(test_dim, test_value_2)

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C isa AbstractScalar
    @test B_proj_C == test_value_1

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == test_value_1
end

@testset "project(B::Real, C::Scalar)" begin
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = test_value_1

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Scalar(test_value_2)

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C isa AbstractScalar
    @test B_proj_C == test_value_1

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == test_value_1
end

@testset "project(B::Real, C::One)" begin
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = test_value

    C = One()

    # return_blade == true
    B_proj_C = project(B, C)
    @test B_proj_C isa AbstractScalar
    @test B_proj_C == test_value

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == test_value
end

@testset "project(B::Real, C::Zero)" begin
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = test_value

    C = Zero()

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

# ------ B::Vector

@testset "project(B::Vector, M::Multivector)" begin
    @test_skip 1
end

@testset "project(B::Vector, C::Blade)" begin
    # --- Preparations

    test_dim = 15
    B = rand(test_dim)

    # --- Tests

    # ------ grade(C) > 1

    # return_blade == true
    C = Blade(rand(test_dim, 5))
    projection_vectors = basis(C) * transpose(basis(C)) * B
    @test project(B, C) ≈ Blade(projection_vectors)

    # return_blade == false
    @test project(B, C, return_blade=false) ≈ projection_vectors

    # ------ grade(C) == 1

    # return_blade == true
    C = Blade(rand(test_dim, 1))
    projection_vectors = B ⋅ basis(C) * basis(C)
    @test project(B, C) ≈ Blade(projection_vectors)

    # return_blade == false
    projection_vectors = (B ⋅ basis(C)) * basis(C)
    @test project(B, C, return_blade=false) ≈ projection_vectors

    # ------ length(B) != dim(C)

    C = Blade(rand(test_dim + 1, 3))
    @test_throws DimensionMismatch project(B, C)
    @test_throws DimensionMismatch project(B, C, return_blade=false)
end

@testset "project(B::Vector, C::Pseudoscalar)" begin
    # --- Preparations

    test_dim = 10
    B = Vector(rand(test_dim))

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Pseudoscalar(test_dim, test_value)

    # --- Tests

    # ------ length(B) == dim(C)

    # return_blade == true
    @test project(B, C) == Blade(B)

    # return_blade == false
    @test project(B, C, return_blade=false) == B

    # ------ length(B) != dim(C)

    B = Vector(rand(test_dim + 1))
    @test_throws DimensionMismatch project(B, C)
    @test_throws DimensionMismatch project(B, C, return_blade=false)
end

@testset "project(B::Vector, C::Scalar)" begin
    # --- Preparations

    test_dim = 15
    B = rand(test_dim)

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Scalar(test_value)

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Vector, C::One)" begin
    # Preparations
    B = Vector(rand(10))
    C = One()

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end

@testset "project(B::Vector, C::Zero)" begin
    B = Vector(rand(10))
    C = Zero()

    # return_blade == true
    B_proj_C = project(B, C)
    @test iszero(B_proj_C)

    # return_blade == false
    B_proj_C = project(B, C, return_blade=false)
    @test B_proj_C isa Real
    @test B_proj_C == 0
end


text:
struct ALView{S,E,N} <: AbstractArray{E,N}
    parent::S
    ALView(parent::S) where S= new{S,site_type(S),length(size(parent))}(parent);
end

function Base.getindex(v::ALView{<:FiniteMPS,E},i::Int)::E where E
    ismissing(v.parent.ALs[i]) && v.parent.CR[i] # by getting CL[i+1], we are garantueeing that AL[i] exists
    return v.parent.ALs[i]
end

function Base.getindex(v::ALView{<:MPSComoving,E},i::Int)::E where E
    i <= length(v.parent) || throw(ArgumentError("out of bounds"))
    i<1 && return v.parent.left_gs.AL[i]
    ALView(v.parent.window)[i]
end

Base.getindex(v::ALView{<:Multiline},i::Int,j::Int) = v.parent[i].AL[j]
Base.setindex!(v::ALView{<:Multiline},vec,i::Int,j::Int) = setindex!(v.parent[i].AL,vec,j);

struct ARView{S,E,N} <: AbstractArray{E,N}
    parent::S
    ARView(parent::S) where S= new{S,site_type(S),length(size(parent))}(parent);
end

function Base.getindex(v::ARView{<:FiniteMPS,E},i::Int)::E where E
    ismissing(v.parent.ARs[i]) && v.parent.CR[i-1] # by getting CL[i], we are garantueeing that AR[i] exists
    return v.parent.ARs[i]
end

function Base.getindex(v::ARView{<:MPSComoving,E},i::Int)::E where E
    i >=1 || throw(ArgumentError("out of bounds"))
    i>length(v.parent) && return v.parent.right_gs.AR[i]
    ARView(v.parent.window)[i]
end

Base.getindex(v::ARView{<:Multiline},i::Int,j::Int) = v.parent[i].AR[j]
Base.setindex!(v::ARView{<:Multiline},vec,i::Int,j::Int) = setindex!(v.parent[i].AR,vec,j);

struct CRView{S,E,N} <: AbstractArray{E,N}
    parent::S
    CRView(parent::S) where S= new{S,bond_type(S),length(size(parent))}(parent);
end

function Base.getindex(v::CRView{<:FiniteMPS,E},i::Int)::E where E
    if ismissing(v.parent.CLs[i+1])
        if i == 0 || !ismissing(v.parent.ALs[i])
            (v.parent.CLs[i+1],temp) = rightorth(_transpose_tail(v.parent.AC[i+1]),alg=LQpos())
            v.parent.ARs[i+1] = _transpose_front(temp);
        else
            (v.parent.ALs[i],v.parent.CLs[i+1]) = leftorth(v.parent.AC[i],alg=QRpos())
        end
    end
    return v.parent.CLs[i+1]
end

function Base.setindex!(v::CRView{<:FiniteMPS},vec,i::Int)
    if ismissing(v.parent.CLs[i+1])
        if !ismissing(v.parent.ALs[i])
            (v.parent.CLs[i+1],temp) = rightorth(_transpose_tail(v.parent.AC[i+1]),alg=LQpos())
            v.parent.ARs[i+1] = _transpose_front(temp);
        else
            (v.parent.ALs[i],v.parent.CLs[i+1]) = leftorth(v.parent.AC[i],alg=QRpos())
        end
    end

    v.parent.CLs.=missing;
    v.parent.ACs.=missing;
    v.parent.ALs[i+1:end].=missing;
    v.parent.ARs[1:i].=missing;

    v.parent.CLs[i+1] = vec
end

Base.getindex(v::CRView{<:MPSComoving},i::Int) = CRView(v.parent.window)[i]
Base.setindex!(v::CRView{<:MPSComoving},vec,i::Int) = Base.setindex!(CRView(v.parent.window),vec,i);
Base.getindex(v::CRView{<:Multiline},i::Int,j::Int) = v.parent[i].CR[j]
Base.setindex!(v::CRView{<:Multiline},vec,i::Int,j::Int) = setindex!(v.parent[i].CR,vec,j);

struct ACView{S,E,N} <: AbstractArray{E,N}
    parent::S
    ACView(parent::S) where S= new{S,site_type(S),length(size(parent))}(parent);
end

function Base.getindex(v::ACView{<:FiniteMPS,E},i::Int)::E where E
    if ismissing(v.parent.ACs[i]) && !ismissing(v.parent.ARs[i])
        c = v.parent.CR[i-1];
        ar = v.parent.ARs[i];
        v.parent.ACs[i] = _transpose_front(c*_transpose_tail(ar))
    elseif ismissing(v.parent.ACs[i]) && !ismissing(v.parent.ALs[i])
        c = v.parent.CR[i];
        al = v.parent.ALs[i];
        v.parent.ACs[i] = al*c;
    end
    return v.parent.ACs[i]
end

function Base.setindex!(v::ACView{<:FiniteMPS},vec::GenericMPSTensor,i::Int)
    if ismissing(v.parent.ACs[i])
        i<length(v) && v.parent.AR[i+1]
        i>1 && v.parent.AL[i-1]
    end

    v.parent.ACs.=missing;
    v.parent.CLs.=missing;
    v.parent.ALs[i:end].=missing;
    v.parent.ARs[1:i].=missing;
    v.parent.ACs[i] = vec;
end

function Base.setindex!(v::ACView{<:FiniteMPS},vec::Tuple{<:GenericMPSTensor,<:GenericMPSTensor},i::Int)
    if ismissing(v.parent.ACs[i])
        i<length(v) && v.parent.AR[i+1]
        i>1 && v.parent.AL[i-1]
    end

    v.parent.ACs.=missing;
    v.parent.CLs.=missing;
    v.parent.ALs[i:end].=missing;
    v.parent.ARs[1:i].=missing;

    (a,b) = vec
    if isa(a,MPSBondTensor) #c/ar
        v.parent.CLs[i] = a;
        v.parent.ARs[i] = b;
    else #al/c
        @assert isa(b,MPSBondTensor)

        v.parent.CLs[i+1] = b;
        v.parent.ALs[i] = a;
    end

end

function Base.getindex(v::ACView{<:MPSComoving,E},i::Int)::E where E
    (i >= 1 && i <= length(v.parent)) || throw(ArgumentError("out of bounds"))
    ACView(v.parent.window)[i]
end
Base.setindex!(v::ACView{<:MPSComoving},vec,i::Int) = Base.setindex!(ACView(v.parent.window),vec,i)
Base.getindex(v::ACView{<:Multiline},i::Int,j::Int) = v.parent[i].AC[j]
Base.setindex!(v::ACView{<:Multiline},vec,i::Int,j::Int) = setindex!(v.parent[i].AC,vec,j);

#--- define the rest of the abstractarray interface
Base.size(psi::Union{ACView,ALView,ARView}) = size(psi.parent);

#=
CRView is tricky. It starts at 0 for finitemps/mpscomoving, but for multiline Infinitemps objects, it should start at 1.
=#
Base.size(psi::CRView{S}) where S <: Union{FiniteMPS,MPSComoving} = (length(psi.parent)+1,);
Base.axes(psi::CRView{S}) where S <: Union{FiniteMPS,MPSComoving} = map(n->0:n-1, size(psi));

Base.size(psi::CRView{<:Multiline{<:InfiniteMPS}}) = size(psi.parent);
Base.size(psi::CRView{<:Multiline{S}}) where S <: Union{FiniteMPS,MPSComoving} = (length(psi.parent.data),length(first(psi.parent.data))+1);
Base.axes(psi::CRView{<:Multiline{S}}) where S <: Union{FiniteMPS,MPSComoving} = (Base.OneTo(length(psi.parent.data)),0:length(first(psi.parent.data)));

#the checkbounds for multiline objects needs to be changed, as the first index is periodic
#however if it is a Multiline(Infinitemps), then the second index is also periodic!
Base.checkbounds(::Type{Bool},psi::Union{ACView{<:Multiline},ALView{<:Multiline},ARView{<:Multiline},CRView{<:Multiline}}, a, b) = first(psi.parent.data) isa InfiniteMPS ? true : checkbounds(Bool, CRView(first(psi.parent.data)), b);


text:
"""
    ModelNet

Base ModelNet for MN10/40.

### Fields:

* `root::String`                           - Root directory of dataset
* `path::String`                           - Directory of dataset
* `train::Bool`                            - Specifies the trainset
* `length::Int`                            - Length of dataset
* `datapaths::Array`                       - Array containing the shape and path for each datapoint
* `transform`                              - Transform to be applied to data point
* `categories::Vector{String}`             - Categories to be used in dataset
* `classes_to_idx::Dict{String, UInt8}`    - Dict mapping from shape name to class_idx
* `idx_to_classes::Dict{UInt8, String}`    - Dict mapping from class_idx to shape name
"""
struct ModelNet <: AbstractDataset
    root::String
    path::String
    train::Bool
    length::Int
    datapaths::Array
    transform::Any
    categories::Vector{String}
    classes_to_idx::Dict{String,UInt8}
    idx_to_classes::Dict{UInt8,String}
end

function _mn_dataset(
    root::String,
    train::Bool,
    download::Bool,
    transform,
    categories::Vector{String},
    variant::Int,
)
    root = normpath(root)
    _path = load_dataset(root, download, variant)
    train ? _split = "train" : _split = "test"
    datapaths = []
    for category in categories
        category in MN10_classes ||
            error("given category: $(category) is not a valid ModelNet10 category.")
        datapath = [
            (category, joinpath(_path, category, _split, filename)) for
            filename in readdir(joinpath(_path, category, _split)) if
            split(filename, ".")[end] == "off"
        ]
        append!(datapaths, datapath)
    end

    classes_to_idx = Dict{String,UInt8}([(categories[i], i) for i = 1:length(categories)])
    idx_to_classes = Dict{UInt8,String}([(i, categories[i]) for i = 1:length(categories)])

    _length = length(datapaths)
    return ModelNet(
        root,
        _path,
        train,
        _length,
        datapaths,
        transform,
        categories,
        classes_to_idx,
        idx_to_classes,
    )
end

function load_dataset(root::String, download::Bool, variant::Int)
    ispath(root) || mkpath(root)
    local_dir = joinpath(root, "ModelNet$(variant)")
    local_path = joinpath(root, "ModelNet$(variant).zip")
    hash_mn10 = "9d8679435fc07d1d26f13009878db164a7aa8ea5e7ea3c8880e42794b7307d51"
    hash_mn40 = "42dc3e656932e387f554e25a4eb2cc0e1a1bd3ab54606e2a9eae444c60e536ac"
    if (variant == 10)
        hash = hash_mn10
    elseif (variant == 40)
        hash = hash_mn40
    end

    if (!isdir(local_dir))
        if (!isfile(local_path))
            if (download)
                download_and_verify(
                    "http://3dvision.princeton.edu/projects/2014/3DShapeNets/ModelNet$(variant).zip",
                    local_path,
                    hash,
                )
            else
                error("dataset not found and auto-download option is set false.")
            end
        end
        run(`unzip -qo $local_path -d $root`)
    end

    return local_dir
end

function Base.getindex(v::ModelNet, idx::Int)
    data = Flux3D.load_trimesh(v.datapaths[idx][2])
    cls = v.classes_to_idx[v.datapaths[idx][1]]
    category = v.datapaths[idx][1]
    if v.transform != nothing
        data = v.transform(data)
    end
    return DataPoint(idx, data, cls, category)
end

Base.size(v::ModelNet) = (v.length,)
Base.length(v::ModelNet) = v.length

function Base.show(io::IO, dset::ModelNet)
    print(
        io,
        "ModelNet Dataset:",
        "\n    root: $(dset.root)",
        "\n    train: $(dset.train)",
        "\n    length: $(dset.length)",
        "\n    transform: $(dset.transform)",
        "\n    categories: $(length(dset.categories))",
    )
end


text:
DimensionalData.name(::Coord) = "Coordinates"

# printing of a single clim array
function Base.show(io::IO, ::MIME"text/plain", A::ClimArray)
    summary(io, A)
    print(io, "and")
    printstyled(io, " data: "; color=:green)
    dataA = gnv(A)
    print(io, summary(dataA), "\n")
    x = 2length(dims(A)) + attriblength(A.attrib) + 5
    custom_show(io, dataA, x)
end

attriblength(d::AbstractDict) = length(d)
attriblength(d) = 0

# Define summary
function Base.summary(io::IO, A::ClimArray)
    printstyled(io, "ClimArray"; color=:blue)
    if A.name ≠ Symbol("")
        print(io, " (named ")
        printstyled(io, A.name; color=:blue)
        print(io, ")")
    end

    DimensionalData.Dimensions.print_dims(io, MIME"text/plain"(), dims(A))
    print(io, "\n")

    if !isnothing(A.attrib)
        printstyled(io, "attributes: "; color=:magenta)
        show(io, MIME"text/plain"(), A.attrib)
        print(io, '\n')
    end
end

# Thanks to Michael Abbott for the following function
function custom_show(io::IO, A::AbstractArray{T,0}, x) where T
    Base.show(IOContext(io, :compact => true, :limit => true), A)
end
function custom_show(io::IO, A::AbstractArray{T,1}, x) where T
    Base.show(IOContext(io, :compact => true, :limit => true, :displaysize => displaysize(io) .- (x, 0)), A)
end
function custom_show(io::IO, A::AbstractArray{T,2}, x) where T
    Base.print_matrix(IOContext(io, :compact => true, :limit => true, :displaysize => displaysize(io) .- (x, 0)), A)
end
function custom_show(io::IO, A::AbstractArray{T,N}, x) where {T,N}
    o = ones(Int, N-2)
    frame = A[:, :, o...]
    onestring = join(o, ", ")
    println(io, "[:, :, $(onestring)]")
    Base.print_matrix(
        IOContext(io, :compact => true, :limit=>true, :displaysize => displaysize(io) .- (x, 0)),
        frame)
    print(io, "\n[and ", prod(size(A,d) for d=3:N) - 1," more slices...]")
end

# printing of a tuple of climarrays
function Base.show(io::IO, A::ClimArray)
    printstyled(io, "ClimArray"; color=:blue)
    if A.name ≠ Symbol("")
        print(io, " (named ")
        printstyled(io, A.name; color=:blue)
        print(io, ")")
    end
    s = join(string.(size(A)), "×")
    d = join([replace(string(basetypeof(d)), "Ti"=>"Time") for d in dims(A)], "×")
    print(io, " with $s $d data")
end
# printing of a vector of climarrays
function Base.show(io::IO, ::MIME"text/plain", Xs::Vector{<:ClimArray})
    println(io, "$(length(Xs))-element Vector{<:ClimArray}:")
    for X in Xs
        println(io, " ", X)
    end
end


text:
# Common functions for GBMs.
#
# Copyright (c) 2020: Jeffrey W. Miller.
# This file is released under the MIT "Expat" License.

module Common

export SS, logprior, Prior, validate_inputs, validate_outputs, sum_of_squares, residuals, tsvd

using Distributions
import PROPACK
import LinearAlgebra
Identity = LinearAlgebra.I

include("outcome_nb.jl"); using .Outcome_NB

# Wrapper for tsvd function in PROPACK.  Re-tries multiple times if errors occur, and suppresses standard output.
function tsvd(A; k=1, trynum=1, maxtries=5)
    if trynum == maxtries
        U,D,V,bnd,nprod,ntprod = PROPACK.tsvd(A; k=k)
    else
        try
            U,D,V,bnd,nprod,ntprod = PROPACK.tsvd(A; k=k)
            return U,D,V
        catch
        end
        U,D,V = tsvd(A; k=k, trynum=trynum+1, maxtries=maxtries)
    end
    return U,D,V
end

minsvdval(A) = PROPACK.tsvdvals(A; k=size(A,2))[1][end]

# Sum of squares
SS(A) = sum(A.*A)

logprior(A,lambda,mu=0.0) = sum(logpdf.(Normal(mu,1/lambda),A))

# Prior structure
#    lambda_* = precision of the Normal prior on the entries of A, B, C, D, U, V, S, or T.
#    mu_S, mu_T = means of the Normal priors on the entries of S and T, respectively, e.g., S[i] ~ Normal(mu_S, 1/lambda_S).
#    The prior mean for A, B, C, D, U, and V is zero, e.g., A[j,k] ~ N(0, 1/lambda_A).
mutable struct Prior
    lambda_A::Float64
    lambda_B::Float64
    lambda_C::Float64
    lambda_D::Float64
    lambda_U::Float64
    lambda_V::Float64
    lambda_S::Float64
    lambda_T::Float64
    mu_S::Float64
    mu_T::Float64
    Prior() = (p = new();
        p.lambda_A = p.lambda_B = p.lambda_C = p.lambda_D = p.lambda_U = p.lambda_V = 1.0;
        p.lambda_S = p.lambda_T = 1.0;
        p.mu_S = p.mu_T = 0.0;
        return p
    )
end

# Check constraints on Y, X, Z, and M
function validate_inputs(Y,X,Z,M)
    @assert(isa(Y,Matrix) && isa(X,Matrix) && isa(Z,Matrix), "Y, X, and Z must be matrices.")
    I,K = size(X)
    J,L = size(Z)
    @assert(I == size(Y,1), "Y and X must have the same number of rows.")
    @assert(J == size(Y,2), "Y and Z' must have the same number of columns.")
    @assert(K > 0, "X must have at least one column.")
    @assert(L > 0, "Z must have at least one column.")
    @assert((M>=0) && isa(M,Int), "M must be a nonnegative integer.")
    @assert(M < min(I,J))
    @assert(all(X[:,1].==1), "The first column of X must be all ones.")
    @assert(all(Z[:,1].==1), "The first column of Z must be all ones.")
    if K>1; @assert(maximum(abs.(sum(X[:,2:end],dims=1))) < 1e-10, "Except for column 1, the columns of X must sum to zero."); end
    if L>1; @assert(maximum(abs.(sum(Z[:,2:end],dims=1))) < 1e-10, "Except for column 1, the columns of Z must sum to zero."); end
    @assert(minsvdval(X) > 0, "X must be full rank, that is, X'*X must be invertible.")
    @assert(minsvdval(Z) > 0, "Z must be full rank, that is, Z'*Z must be invertible.")
end

# Check identifiability constraints and interpretation identities
function validate_outputs(Y,X,Z,A,B,C,D,U,V,S,T)
    I,M = size(U)
    @assert(maximum(abs.(X'*B)) < 1e-10, "X'*B is nonzero.")
    @assert(maximum(abs.(Z'*A)) < 1e-10, "Z'*A is nonzero.")
    if M > 0
        @assert(maximum(abs.(X'*U)) < 1e-10, "X'*U is nonzero.")
        @assert(maximum(abs.(Z'*V)) < 1e-10, "Z'*V is nonzero.")
        @assert(maximum(abs.(U'*U - Identity)) < 1e-8, "U is not orthonormal, that is, U'*U-Identity is nonzero.")
        @assert(maximum(abs.(V'*V - Identity)) < 1e-8, "V is not orthonormal, that is, V'*V-Identity is nonzero.")
        if M > 1; @assert(maximum(diff(D)) <= 0, "The entries of D are not in decreasing order."); end
        # (could check sign constraints here)
    end
    @assert(abs(mean(exp.(S)) .- 1) < 1e-10, "The identifiability constraint on S is not satisfied.")
    @assert(abs(mean(exp.(T)) .- 1) < 1e-10, "The identifiability constraint on T is not satisfied.")
    logMu = X*A' + B*Z' + X*C*Z' + U*(D.*V')
    @assert(abs.(mean(logMu) - C[1,1]) < 1e-12)
    @assert(maximum(abs.(vec(mean(logMu, dims=2)) - B[:,1] - X*C[:,1])) < 1e-12, "One or more interpretability identities are not satisfied.")
    @assert(maximum(abs.(vec(mean(logMu, dims=1)) - A[:,1] - Z*C[1,:])) < 1e-12, "One or more interpretability identities are not satisfied.")
end

# Compute sum-of-squares of each model component
function sum_of_squares(X,Z,A,B,C,D,U,V)
    return [SS(X*A'), SS(B*Z'), SS(X*C*Z'), SS(U*(D.*V'))]
end

# Compute residuals (or adjusted residuals)
#    rx = indices of X covariates to retain in residuals (do not adjust them out)
#    rz = indices of Z covariates to retain in residuals (do not adjust them out)
#    ru = indices of latent factors to retain in residuals (do not adjust them out)
function residuals(Y,X,Z,A,B,C,D,U,V,S,T,omega; rx=[], rz=[], ru=[])
    logMu = X*A' + B*Z' + X*C*Z' + U*(D.*V')
    logMu_retain = X[:,rx]*A[:,rx]' + B[:,rz]*Z[:,rz]' + X[:,rx]*C[rx,rz]*Z[:,rz]' + U[:,ru]*(D[ru].*V[:,ru]')
    Eps = log.(Y .+ 0.125) .- logMu
    Eps = logMu_retain .+ Eps
    r = exp.(-(S .+ T' .+ omega))
    Mu,W,E = compute_MuWE(Y,logMu,r)
    Sigma_Eps = 1.0./sqrt.(W)
    return Eps,Sigma_Eps
end

end # module Common


text:
 n=rand(1:10)
 @test diag(Matrix(Toeplitz(collect(-n:n)))) == zeros(n+1)


text:
immutable Animation
  dir::ASCIIString
  frames::Vector{ASCIIString}
end

function Animation()
  tmpdir = convert(ASCIIString, mktempdir())
  Animation(tmpdir, ASCIIString[])
end

function frame{P<:PlottingObject}(anim::Animation, plt::P=current())
  i = length(anim.frames) + 1
  filename = @sprintf("%06d.png", i)
  png(plt, joinpath(anim.dir, filename))
  push!(anim.frames, filename)
end

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

"Wraps the location of an animated gif so that it can be displayed"
immutable AnimatedGif
  filename::ASCIIString
end

function gif(anim::Animation, fn::@compat(AbstractString) = "tmp.gif"; fps::Integer = 20)
  fn = abspath(fn)

  try

    # high quality
    speed = round(Int, 100 / fps)
    run(`convert -delay $speed -loop 0 $(anim.dir)/*.png $fn`)

  catch err
    warn("Tried to create gif using convert (ImageMagick), but got error: $err\nWill try ffmpeg, but it's lower quality...)")

    # low quality
    run(`ffmpeg -v 0 -framerate $fps -i $(anim.dir)/%06d.png -y $fn`)
    # run(`ffmpeg -v warning -i  "fps=$fps,scale=320:-1:flags=lanczos"`)
  end

  info("Saved animation to ", fn)
  AnimatedGif(fn)
end

# write out html to view the gif... note the rand call which is a hack so the image doesn't get cached
function Base.writemime(io::IO, ::MIME"text/html", agif::AnimatedGif)
  write(io, "<img src=\"$(relpath(agif.filename))?$(rand())>\" />")
end


text:
# idea: use a vector of Nodes, where each Node is a operator and arguments
using DataStructures
import Base: convert, promote_rule, +, -, *, /, >, <, log

immutable RDiff <: Real
    val::Float64
    adj::Float64

    RDiff(val, adj) = new(val, adj)
end
typealias Converts Union{Float64, Int, Irrational, Rational}
RDiff(x::Converts) = RDiff(x, 0.)

immutable Node
    inputs::Vector{RDiff}
    outputs::Vector{RDiff}
    op::Function
    container::Vector{Node}
end

convert(::Type{RDiff}, x::Converts) = RDiff(x)
promote_rule(::Type{Bool}, ::Type{RDiff}) = RDiff
promote_rule(::Type{Float64}, ::Type{RDiff}) = RDiff
promote_rule{T<:Integer}(::Type{T}, ::Type{RDiff}) = RDiff

value(x::RDiff) = x.val
adjoint(x::RDiff) = x.adj
>(x::RDiff, y::RDiff) = x.val > y.val
<(x::RDiff, y::RDiff) = x.val < y.val

function grad(f)
    global varstack = Vector{Node}()

    function ∇f(x)
        y = map(x -> RDiff(x, 0., varstack), x)
        f(y)
        varstack[end].outputs[1].adj = 1.
        for n in varstack
            backprop(n)
        end
        map(value, y)
    end

    return ∇f

end

macro diff_rule(sig, code)
    fn = _fn_of_vals(sig)
    argsyms = Symbol[s.args[1] for s in sig.args[2:end]]

    ex = quote
        function $(sig.args[1])($(sig.args[2:end]...))
            vstack = get_container($(argsyms...))
            this = RDiff($fn, 0.0, vstack)
            function backprop()
                $code
                nothing
            end
            push!(vstack, this)
            this
        end
    end
    esc(ex)
end

_fn_of_vals(x) = x
function _fn_of_vals(ex::Expr)
    if ex.head == :(::)
        out = :(value($(ex.args[1])))
    else
        out = copy(ex)
        for i in eachindex(out.args)
            out.args[i] = _fn_of_vals(out.args[i])
        end
    end
    out
end

function get_container(vars...)
    has_stack = filter(x -> !isnull(x.stack), vars)
    get(first(has_stack).stack)
end

@diff_rule log(x::RDiff) begin
    x.adj += this.adj / x.val
end

@diff_rule +(x::RDiff, y::RDiff) begin
    x.adj += this.adj
    y.adj += this.adj
end

@diff_rule -(x::RDiff, y::RDiff) begin
    x.adj += this.adj
    y.adj += -this.adj
end

@diff_rule *(x::RDiff, y::RDiff) begin
    x.adj += this.adj * y.val
    y.adj += this.adj * x.val
end

@diff_rule /(x::RDiff, y::RDiff) begin
    x.adj += this.adj / y.val
    y.adj += -this.adj * x.val / y.val^2
end

@diff_rule -(x::RDiff) begin
    x.adj += -this.adj
end


text:
# Autogenerated wrapper script for libpsl_jll for i686-linux-musl
export libpsl, psl

using Libiconv_jll
JLLWrappers.@generate_wrapper_header("libpsl")
JLLWrappers.@declare_library_product(libpsl, "libpsl.so.5")
JLLWrappers.@declare_executable_product(psl)
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll)
    JLLWrappers.@init_library_product(
        libpsl,
        "lib/libpsl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        psl,
        "bin/psl",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()