text:
# SPDX-License-Identifier: MIT

"""
    AbstractNormal{D} <: Distribution

Abstract type for Normal distribution implementations.
"""
abstract type AbstractNormal{D} <: Distribution end

function DefaultNormalParameter(μ::AbstractVector{T},
                                Σ::AbstractMatrix{T}) where T
    Λ = pdmat_inverse(Σ)
    ξ = vcat(Λ*μ, -T(.5)*diag(Λ), -vec_tril(Λ))
    DefaultParameter(ξ)
end

"""
    mutable struct Normal{D} <: AbstractNormal{D}
        param::P where P <: AbstractParameter
    end

Normal distribution with full covariance matrix.

# Constructors

Normal(μ, Σ)

where `μ` is the mean and `Σ` is the covariance matrix.

# Examples
```jldoctest
julia> Normal([1.0, 1.0], [2.0 0.5; 0.5 1.0])
Normal{2}:
  μ = [0.9999999999999998, 1.0]
  Σ = [2.0 0.5; 0.5 1.0]
```
"""
mutable struct Normal{D} <: AbstractNormal{D}
    param::P where P <: AbstractParameter
end

function Normal(μ, Σ)
    param = DefaultNormalParameter(μ, Σ)
    P = typeof(param)
    D = length(μ)
    Normal{D}(param)
end

_unpack(D, v) = v[1:D], v[D+1:2*D], v[2*D+1:end]

#######################################################################
# Distribution interface.

function basemeasure(::AbstractNormal, x::AbstractVector{T}) where T
    -T(.5) * length(x) * log(T(2π))
end

function lognorm(n::AbstractNormal{D},
                 η::AbstractVector{T} = naturalform(n.param)) where {T,D}
    η₁, η₂, η₃ = _unpack(D, η)
    H₂ = matrix(η₂, T(.5)*η₃)
    -T(.5)*pdmat_logdet(-T(2)*H₂) + T(.25)*dot(η₁, pdmat_inverse(-H₂) * η₁)
end

function gradlognorm(n::AbstractNormal)
    μ, Σ = stdparam(n, naturalform(n.param))
    M = Σ + μ*μ'
    vcat(μ, diag(M), vec_tril(M))
end

function sample(n::AbstractNormal{D}, size) where D
    μ, Σ = stdparam(n)
    L = cholesky(Σ).L
    ϵ = randn!(similar(μ, D, size))
    μ .+ L*ϵ
end

splitgrad(n::AbstractNormal{D}, m) where D = _unpack(D, m)

function stats(::AbstractNormal{D}, x::AbstractVector) where D
    xxᵀ = x * x'
    vcat(x, diag(xxᵀ), vec_tril(xxᵀ))
end

function stdparam(n::AbstractNormal{D},
                  η::AbstractVector{T} = naturalform(n.param)) where {T,D}
    η₁, η₂, η₃ = _unpack(D, η)
    H₂ = matrix(η₂, T(.5)*η₃)
    Σ = pdmat_inverse(-T(2)*H₂)
    (μ = Σ*η₁, Σ = Σ)
end

max_stars_repo_name: lucasondel/ExpFamilyDistributions-jl

max_issues_repo_name: BUTSpeechFIT/ExpFamilyDistributions

max_forks_repo_name: lucasondel/ExpFamilyDistributions


text:
# initialize DE based on priors, nDays, etc.
function initialize_process_model(;times, PAR, OSat, day, GPP_daily, ER_daily, kGas = 1.05, Q = 1.794, V = 27.0, dayStart = 4, u0)

    if convert(Int64, floor(maximum(times)/24.0 - dayStart/24.0) + 1) > length(GPP_daily)
        error("The simulation timespan contains more days than are provided in GPP_daily")
    end

    p = Dict(
        "PAR" => PAR,
        "OSat" => OSat,
        "day" => day,
        "GPP_daily"=> GPP_daily,
        "ER_daily" => ER_daily,
        "kGas" => kGas,
        "Q" => Q,
        "V" => V)

    tspan = (minimum(times), maximum(times))

    # Define and solve ODE problem
    prob = ODEProblem(mesocosm_metabolism!, u0, tspan, p)
    return prob
end

ext: jl

max_stars_repo_name: elizabethjmohr/mesocosm-metabolism

max_issues_repo_name: elizabethjmohr/mesocosm-metabolism

max_forks_repo_name: elizabethjmohr/mesocosm-metabolism


text:
using Trilinos
using Test
using MPI

if Test.get_testset_depth() == 0
  MPI.Init()
end

lows_factory = Thyra.BelosLinearOpWithSolveFactory(Float64)
pl = Teuchos.ParameterList()
Thyra.setParameterList(lows_factory, pl)

Thyra.setVerbLevel(lows_factory, Teuchos.VERB_NONE)
lows = Thyra.createOp(lows_factory)

comm = Teuchos.MpiComm(MPI.CComm(MPI.COMM_WORLD))
const n = 20
rowmap = Tpetra.Map(n, 0, comm)
A = Tpetra.CrsMatrix(spdiagm(0 => ones(n)),rowmap)

rangespace = Thyra.tpetraVectorSpace(Tpetra.getRangeMap(A))
domainspace = Thyra.tpetraVectorSpace(Tpetra.getDomainMap(A))
A_thyra = Thyra.tpetraLinearOp(rangespace, domainspace, A)
Thyra.initializeOp(lows_factory, A_thyra, lows)

x = Tpetra.Vector(Tpetra.getDomainMap(A))
b = Tpetra.Vector(Tpetra.getRangeMap(A))
Tpetra.randomize(x)

# Test apply
Tpetra.apply(A,x,b)

x_th = Thyra.tpetraVector(domainspace, x)
b_th = Thyra.tpetraVector(rangespace, b)

status = Thyra.solve(lows, Thyra.NOTRANS, b_th, x_th)

# high-level interface
lows2 = Thyra.LinearOpWithSolve(A)
x2 = lows2 \ b

# Check result
bv = Tpetra.device_view(b)
xv = Tpetra.device_view(x2)
@test length(LinearIndices(bv)) == length(LinearIndices(xv))
for (bi,xi) in zip(bv,xv)
  @test bi ≈ xi
end

if Test.get_testset_depth() == 0
  MPI.Finalize()
end


text:
# # How to plot multiple results in subplots for comparisons

# See [How to set up plots](3.0_set_up_plots.md) to get started

# ### Run the simulation and get various results

# ### Plot the results as an array

# ```julia
# plots = stack_plot([results_one, results_two])
# ```

# ![this one](plots-6/P__ThermalStandard_Stack.png)
# ![this one](plots-6/Comparison.png)

# ### Plot by Fuel Type

# ```julia
# plots = fuel_plot([results_one, results_two], c_sys5_re)
# ```

# ![this one](plots-6/Comparison_Stack.png)
# ![this one](plots-6/Comparison_Bar.png)

# ## Multiple results can be compared while also plotting fewer variables

# ```julia
# variables = [Symbol("P__ThermalStandard")]
# plots = stack_plot([results_one, results_two], variables)
# ```

# ![this one](plots-6/P__ThermalStandard_Stack.png)
# ![this one](plots-6/Comparison_with_fewer_variables.png)


text:
using gazel
using Test

@testset "gazel.jl" begin
    # Write your own tests here.
end


text:
using Test
using Example2018

# Polynomials form a ring
# 1. addition is associative
# 2. additive identity
# 3. additive inverse
# 4. addition commutes
#HW add stuff

@testset "Properties polynomial additions are all satisfied." begin

	z = Polynomial([0])
	zz = Polynomial([0,0])

	for n in 1:100
		deg = rand(0:100)
		p = Polynomial(rand(-10000000:10000000,deg+1))
		q = Polynomial(rand(-10000000:10000000,deg+1))
		n = Polynomial(rand(-10000000:10000000,deg+1))

		@test p + z == p #Additive identity element
		@test z + p == p #Additive identity element commutes
		@test p + (-p) == z #Additive inverse
		@test p + (q + n) == (p + q) + n #Addition is associative
		@test q + zz == q #Additive identity element of different length
		@test zz + q == q #Additive identity element commutes
		@test n + (-n) == zz #Additive inverse (of different length indentity)
	end
end

1 == 1
@test 2 == 2


text:
"""
    densityplot(x, y; kw...)

# Description

Draws a density plot for the given points.

The first vector `x` should contain the horizontal positions for all the points.
The second vector `y` should contain the corresponding vertical positions respectively.
The two vectors must thus be of the same length and ordering.

# Usage

    densityplot(x, y; $(keywords(; remove = (:grid,))))

# Arguments

$(arguments(; add = (:x, :y), remove = (:grid,)))

# Returns

A plot object of type `Plot{DensityCanvas}`.

# Author(s)

- Christof Stocker (github.com/Evizero)

# Examples

```julia-repl
julia> densityplot(randn(1000), randn(1000), title = "Density Plot")
                      Density Plot
        ┌────────────────────────────────────────┐
    2.9 │                    ░                   │
        │                ░    ░  ░               │
        │                ░░░  ░ ░ ░      ░       │
        │             ░░░ ░░▒░▒░░░ ▒ ░           │
        │          ░░ ░░ ░▒░░▒▒▒▒░░░░░░░         │
        │         ░░ ░▒░░░░░▓▒▒░▒░▒░▓░▒░░░       │
        │        ░ ░░░░▓▓▒▓▒▒█▒▓▒▒▓▒▒░▒ ░ ░      │
        │           ░ ▒▒▒▒▓▓▒▓▓▓▓▒█▒▒░░░░░░ ░    │
        │     ░   ░░░░░░░▒▒▓░░▓▒█▒▒▓▓▒▒ ░    ░   │
        │           ░ ░░▒▒░░▒░░▒░░░░ ░░  ░ ░     │
        │         ░   ░ ░ ░░░░░░ ░░ ░░ ░         │
        │              ░  ░ ░ ░  ░ ░   ░░░       │
        │                ░░  ░                   │
        │                                        │
   -3.3 │                                        │
        └────────────────────────────────────────┘
        -3.4                                   2.9
```

# See also

[`Plot`](@ref), [`scatterplot`](@ref), [`DensityCanvas`](@ref)
"""
function densityplot(
    x::AbstractVector,
    y::AbstractVector;
    color::UserColorType = KEYWORDS.color,
    grid = false,
    name = KEYWORDS.name,
    kw...,
)
    plot = Plot(x, y, DensityCanvas; grid = grid, kw...)
    scatterplot!(plot, x, y; color = color, name = name)
end

densityplot!(plot::Plot{<:DensityCanvas}, x::AbstractVector, y::AbstractVector; kw...) =
    scatterplot!(plot, x, y; kw...)


text:
# This code is derived from `julia-config.jl` (part of Julia) and should be kept aligned with it.

threadingOn() = ccall(:jl_threading_enabled, Cint, ()) != 0

function shell_escape(str)
    str = replace(str, "'" => "'\''")
    return "'$str'"
end

function libDir()
    return if ccall(:jl_is_debugbuild, Cint, ()) != 0
        dirname(abspath(Libdl.dlpath("libjulia-debug")))
    else
        dirname(abspath(Libdl.dlpath("libjulia")))
    end
end

private_libDir() = abspath(Sys.BINDIR, Base.PRIVATE_LIBDIR)

function includeDir()
    return abspath(Sys.BINDIR, Base.INCLUDEDIR, "julia")
end

function ldflags()
    fl = "-L$(shell_escape(libDir()))"
    if Sys.iswindows()
        fl = fl * " -Wl,--stack,8388608"
    elseif Sys.islinux()
        fl = fl * " -Wl,--export-dynamic"
    end
    return fl
end

function ldlibs()
    libname = if ccall(:jl_is_debugbuild, Cint, ()) != 0
        "julia-debug"
    else
        "julia"
    end
    if Sys.isunix()
        return "-Wl,-rpath,$(shell_escape(libDir())) -Wl,-rpath,$(shell_escape(private_libDir())) -l$libname"
    else
        return "-l$libname -lopenlibm"
    end
end

function cflags()
    flags = IOBuffer()
    print(flags, "-std=gnu99")
    include = shell_escape(includeDir())
    print(flags, " -I", include)
    if threadingOn()
        print(flags, " -DJULIA_ENABLE_THREADING=1")
    end
    if Sys.isunix()
        print(flags, " -fPIC")
    end
    return String(take!(flags))
end

function allflags()
    return "$(cflags()) $(ldflags()) $(ldlibs())"
end

const short_flag_to_jloptions = Dict(
    "C" => :cpu_target,
    "J" => :image_file,
    "O" => :opt_level,
    "g" => :debug_level
)

const flag_to_jloptions = Dict(
    "check-bounds" => :check_bounds,
    "code-coverage" => :code_coverage,
    "compile" => :compile_enabled,
    "compiled-modules" => :use_compiled_modules,
    "cpu-target" => :cpu_target,
    "depwarn" => :depwarn,
    "handle-signals" => :handle_signals,
    "history-file" => :historyfile,
    "machine-file" => :machine_file,
    "math-mode" => :fast_math,
    "optimize" => :opt_level,
    "output-bc" => :outputbc,
    "output-ji" => :outputji,
    "output-jitbc" => :outputjitbc,
    "output-o" => :outputo,
    "output-unoptbc" => :outputunoptbc,
    "project" => :project,
    "startup-file" => :startupfile,
    "sysimage" => :image_file,
    "sysimage-native-code" => :use_sysimage_native_code,
    "trace-compile" => :trace_compile,
    "track-allocation" => :malloc_log,
    "warn-overwrite" => :warn_overwrite,
    "inline" => :can_inline
)

const jl_options_to_flag = Dict(
    :can_inline => "inline",
    :handle_signals => "handle-signals",
    :opt_level => "optimize",
    :depwarn => "depwarn",
    :malloc_log => "track-allocation",
    :outputo => "output-o",
    :startupfile => "startup-file",
    :compile_enabled => "compile",
    :trace_compile => "trace-compile",
    :check_bounds => "check-bounds",
    :outputji => "output-ji",
    :use_sysimage_native_code => "sysimage-native-code",
    :outputunoptbc => "output-unoptbc",
    :historyfile => "history-file",
    :outputbc => "output-bc",
    :warn_overwrite => "warn-overwrite",
    :machine_file => "machine-file",
    :code_coverage => "code-coverage",
    :image_file => "sysimage",
    :cpu_target => "cpu-target",
    :outputjitbc => "output-jitbc",
    :project => "project",
    :fast_math => "math-mode",
    :use_compiled_modules => "compiled-modules",
    :debug_level => "g"
)

const flags_with_cmdval = Set([
    :handle_signals,
    :use_sysimage_native_code,
    :depwarn,
    :can_inline,
    :historyfile,
    :startupfile,
    :use_compiled_modules,
    :warn_overwrite,
    :check_bounds,
    :fast_math,
    :compile_enabled,
    :malloc_log,
    :code_coverage
])

function to_cmd_val(key::Symbol, val)
    # undocumented auto!? well we can only skip it I guess
    if key in (:depwarn, :warn_overwrite)
        val == 0  && return "no"
        val == 1 && return "yes"
        val == 2 && return "error"
    end
    if key in (:code_coverage, :malloc_log)
        val == 0 && return "none"
        val == 1 && return "user"
        val == 2 && return "all"
    end
    if key == :compile_enabled
        val == 0 && return "no"
        val == 1 && return "yes"
        val == 2 && return "all"
    end
    if key == :fast_math
        val == 0 && return "ieee"
        val == 1 && return "fast"
    end
    val in (0, -1) && return ""
    val == 1 && return "yes"
    val == 2 && return "no"
end

function jl_option_value(opts, key)
    value = getfield(opts, key)
    if value isa Ptr{UInt8}
        return value == C_NULL ? "" : unsafe_string(value)
    end
    if key in flags_with_cmdval
        return to_cmd_val(key, value)
    end
    return value
end

is_short_flag(flag) = haskey(short_flag_to_jloptions, flag)

function jl_option_key(flag::Symbol)
    # if symbol is used, we can also check the fields directly.
    # TODO should we also do this for strings?
    flag in fieldnames(Base.JLOptions) && return flag
    jl_option_key(string(flag))
end

function jl_option_key(flag::String)
    haskey(short_flag_to_jloptions, flag) && return short_flag_to_jloptions[flag]
    haskey(flag_to_jloptions, flag) && return flag_to_jloptions[flag]
    m_flag = replace(flag, "_" => "-")
    haskey(flag_to_jloptions, m_flag) && return flag_to_jloptions[m_flag]
    flags = replace.([keys(flag_to_jloptions)..., keys(short_flag_to_jloptions)...], ("-" => "_",))
    error("Flag $flag not a valid Julia Options. Valid flags are:\n$(join(flags, " "))")
end

"""
    extract_flag(flag, jl_cmd = Base.julia_cmd())

Extracts the value for `flag` from `jl_cmd`.
"""
function jl_flag_value(flag, jl_options = Base.JLOptions())
    jl_option_value(jl_options, jl_option_key(flag))
end

current_systemimage() = jl_flag_value("J")

"""
    run_julia(
        code::String;
        g = nothing, O = nothing, output_o = nothing, J = nothing,
        startup_file = "no"
    )

Run `code` in a julia command.
You can overwrite any julia command line flag by setting it to a value.
If the flag has the value `nothing`, the value of the flag of the current julia process is used.
"""
function run_julia(code::String; kw...)
    run(julia_code_cmd(code; kw...))
end

function jl_command(flag, value)
    (value === nothing || isempty(value)) && return ""
    if is_short_flag(flag)
        string("-", flag, value)
    else
        string("--", flag, "=", value)
    end
end

function push_command!(cmd, flag, value)
    command = jl_command(flag, value)
    isempty(command) || push!(cmd.exec, command)
end

function julia_code_cmd(
        code::String, jl_options = Base.JLOptions();
        kw...
    )
    jl_cmd = Base.julia_cmd()
    cmd = `$(jl_cmd.exec[1])` # extract julia exe
    # Add the commands from the keyword arguments
    for (k, v) in kw
        jl_key = jl_option_key(k)
        flag = jl_options_to_flag[jl_key]
        push_command!(cmd, flag, v)
    end
    # add remaining commands from JLOptions
    for key in setdiff(keys(jl_options_to_flag), keys(kw))
        flag = jl_options_to_flag[key]
        push_command!(cmd, flag, jl_option_value(jl_options, key))
    end
    # for better debug, let's not make a tmp file which would get lost!
    file = sysimg_folder("run_julia_code.jl")
    open(io-> println(io, code), file, "w")
    push!(cmd.exec, file)
    cmd
end


text:
using DataFrames, TableOperations, Tables, Random
using TreeParzen
using Recommenders:
    Movielens1M,
    load_dataset,
    ratio_split,
    SLIM,
    evaluate_u2i,
    MeanPrecision,
    MeanRecall,
    MeanNDCG

function main()
    ml1M = Movielens1M()
    download(ml1M)
    rating, _, _ = load_dataset(ml1M)

    rating = rating |> TableOperations.transform(Dict(:rating => x -> 1.0))

    Random.seed!(1234)
    train_valid_table, test_table = ratio_split(rating, 0.8)
    train_table, valid_table = ratio_split(train_valid_table, 0.8)

    prec10 = MeanPrecision(10)
    recall10 = MeanRecall(10)
    ndcg10 = MeanNDCG(10)
    metrics = [prec10, recall10, ndcg10]

    space = Dict(
        :λminratio => HP.LogUniform(:λminratio, log(1e-6), log(0.5)),
        :l1_ratio => HP.LogUniform(:l1_ratio, log(1e-6), log(0.9999)),
        :k => HP.QuantUniform(:k, 1.0, 1000.0, 1.0),
    )

    function invert_output(params)
        @info params
        λminratio = params[:λminratio]
        l1_ratio = params[:l1_ratio]
        k = convert(Int, params[:k])

        model = SLIM(l1_ratio, λminratio, k)

        result = evaluate_u2i(
            model,
            train_table,
            valid_table,
            metrics,
            10,
            col_item = :movieid,
            drop_history = true,
            cd_tol = 1e-4,
            nλ = 10,
        )
        @info result
        return -result[:ndcg10]
    end

    @info "Tuning start."
    best = fmin(invert_output, space, 100, logging_interval = -1)
    @info best

    @info "Evaluate best model."

    best_model = SLIM(best[:l1_ratio], best[:λminratio], convert(Int, best[:k]))
    result = evaluate_u2i(
        best_model,
        train_valid_table,
        test_table,
        metrics,
        10,
        col_item = :movieid,
        drop_history = true,
        cd_tol = 1e-4,
        nλ = 10,
    )
    @info result

end

main()


text:
module Alice

##############################################################################
##
## Dependencies and Reexports
##
##############################################################################

using Reexport
@reexport using Gadfly
@reexport using Distributions
@reexport using Images
using NLopt

using Base.Threads
import Base.show
import Base.size

##############################################################################
##
## Exported methods and types
##
##############################################################################

export Data, InputLayer, FullyConnectedLayer, WordEmbeddingLayer,
       SparseEncoderLayer, ConvolutionLayer, MeanPoolLayer, MaxPoolLayer,
       LinearOutputLayer, MultiLinearOutputLayer, LogisticOutputLayer,
       SoftmaxOutputLayer, NeuralNet, tinynet, check_gradients,
       load_ngrams, load_patches, load_train_subset, load_test_subset,
       load_features,
       fwd_prop!, bwd_prop!, train, train_nlopt,
       loss, val_loss, loss_and_accuracy, val_loss_and_accuracy, accuracy,
       predictions,
       size, viewbatch, display_nearest_words, predict_next_word,
       plot_loss_history, display_rgb_cols, display_rgb_weights

##############################################################################
##
## Load source files
##
##############################################################################

include("layers.jl")
include("init_params.jl")
include("activations.jl")
include("convolutions.jl")
include("pooling.jl")
include("cost_and_delta.jl")
include("regularisation.jl")
include("backprop.jl")
include("train.jl")
include("train_nlopt.jl")
include("performance.jl")
include("utils.jl")
include("check_gradients.jl")
include("plotting.jl");

end # module


text:
#
#  The following functions allow us to parse files, write files, etc
#
#	using DelimitedFiles, DataFrames
#
# 	Classes:
# 		TrialData
#
#
#.  extract_data(path)
#	timestamp_now()
# 	printFigure(name; fig=gcf())
#.  check_imported_data(data::TrialData; idx=nothing)
# 	truncate_at_cue_and_lickp250(td::TrialData; cbuffer_s=0., lbuffer_s=0)

using DelimitedFiles
using DataFrames
using Dates
using CSV

function refresh_tools(path; exact=false)
    # Extracts dependencies from folder
    start_dir = pwd()
    try
        cd(path)
        a = readdir()
        print("Refreshed ")
        for i = 1:length(a)
            if occursin(".jl", a[i])
            	if exact
	                include(joinpath(path, a[i]))
                else
                	include(a[i])
                end
                print(a[i], " ")
            end
        end
        cd(start_dir)
    catch
        cd(start_dir)
        rethrow()
    end
end;

struct TrialData
    xdata::Array{Vector{Float64},1}
    ydata::Array{Vector{Float64},1}
    trialNo::Vector{Int}
    lickTime_s::Vector{Float64}
    path::String
    sessionCode::String
end

function extract_data(path; blmode=false, LOImode=false)
	#
	# Extracts single trial CSV data from folder
	#
    start_dir = pwd()
    cd(path)
    a = readdir()
    nx2 = length(a)
    xstart::Int=0
    ystart::Int=floor(Int, nx2/2)
    xs = []
    ys = []
    trialNo = Vector{Int}(undef, ystart)
    lickTime_s = Vector{Float64}(undef, ystart)
    for i = 1:ystart
        xx = vec(readdlm(a[xstart+i], ',', Float64, '\n'));
        yy = vec(readdlm(a[ystart+i], ',', Float64, '\n'));
        tt = split(a[xstart+i], "_")
        if blmode # if baseline, we have extra baseline text on the filename
        	tNo = parse(Int, tt[3][2:end])
	        lt = parse(Float64, tt[4][3:end-1])
        elseif LOImode
        	tNo = parse(Int, tt[2][2:end])
	        lt = parse(Float64, tt[3][3:end-1])
        else
	        tNo = parse(Int, tt[2][2:end])
	        lt = parse(Float64, tt[3][3:end-1])
        end
        push!(xs,xx)
        push!(ys,yy)
        trialNo[i] = tNo
        lickTime_s[i] = lt
    end
    if blmode
    	tt = split(a[xstart+1], "_")
    	sessionCode = join(tt[5:end])
	    sessionCode = sessionCode[1:end-4]

    else
    	tt = split(a[xstart+1], "_")
	    sessionCode = join(tt[4:end])
	    sessionCode = sessionCode[1:end-4]
    end
    data = TrialData(xs,ys,trialNo,lickTime_s,path, sessionCode)
    cd(start_dir)
    return data
end;

function timestamp_now()
    dform = Dates.DateFormat("yyyy-mm-dd_HHMM")
    return Dates.format(now(), dform)
end

function printFigure(name; fig=gcf(), figurePath=figurePath, verbose=false)
	if Sys.iswindows()
		dlim="\" #"
	else
		dlim="/"
	end

    # gcf() = fig
    # fig
    timestamp = timestamp_now()
    figname = join([figurePath,dlim, name, "_", timestamp, ".eps"])
    if verbose
	    println(figname)
    end
    fig.savefig(figname, transparent=true, format="eps")
end

function check_imported_data(data::TrialData; idx=nothing)
    if length(data.xdata[end]) != length(data.ydata[end])
        error("data fields not matching")
    end
    if isnothing(idx)
        idx = rand(1:length(data.xdata))
    else
        idx = findall(x->x==idx, data.trialNo)[1]
        println(idx)
    end
    xs = data.xdata[idx]
    ys = data.ydata[idx]
    tNo = data.trialNo[idx]
    lt = data.lickTime_s[idx]
    path = data.path
    seshcode = data.sessionCode
    figure(figsize=(5,3))
    plot(xs, ys, color="red", linewidth=2)
    xlabel("time (s)")
    ylabel("dF/F")
    title(join([seshcode," t=", tNo, " lt=", lt, "s"]));
    println("mean: ", mean(ys))
    println("median: ", median(ys))
    println(path)
    return idx
end;

function makeSessionDataFrame(data::TrialData; normalize=false, includeBL_LOI=false, baseline_data=[], LOI_data=[])
    # useful data from the TrialData struct to add to dataframe
    xs = data.xdata
    ys = data.ydata
    tNo = data.trialNo
    lt = data.lickTime_s
    path = data.path
    seshcode = data.sessionCode
    if includeBL_LOI
    	meanbl = [mean(x) for x in baseline_data.ydata]
    	medianbl = [median(x) for x in baseline_data.ydata]
    	meanloi = [mean(x) for x in LOI_data.ydata]
    	medianloi = [median(x) for x in LOI_data.ydata]
	end

    # start by sorting the trial numbers so we can get prev trial info...
    sorted_trial_idxs = sortperm(tNo)
    tNo_sorted = tNo[sorted_trial_idxs]
    lt_sorted_nonorm = lt[sorted_trial_idxs]
    xs_sorted = xs[sorted_trial_idxs]

    if normalize
        v = reduce(vcat, ys)
        ys = [normalize_vector_0_1(x, mx=maximum(v), mn=minimum(v)) for x in ys]
        v = reduce(vcat, lt)
        lt = [normalize_vector_0_1(x, mx=maximum(v), mn=minimum(v)) for x in lt]
        if includeBL_LOI
        	v = reduce(vcat, meanbl)
        	meanbl = [normalize_vector_0_1(x, mx=maximum(v), mn=minimum(v)) for x in meanbl]
        	v = reduce(vcat, medianbl)
        	medianbl = [normalize_vector_0_1(x, mx=maximum(v), mn=minimum(v)) for x in medianbl]
        	v = reduce(vcat, meanloi)
        	meanloi = [normalize_vector_0_1(x, mx=maximum(v), mn=minimum(v)) for x in meanloi]
        	v = reduce(vcat, medianloi)
        	medianloi = [normalize_vector_0_1(x, mx=maximum(v), mn=minimum(v)) for x in medianloi]
		end
    end

    lt_sorted = lt[sorted_trial_idxs]
    ys_sorted = ys[sorted_trial_idxs]
    if includeBL_LOI
		meanbl_sorted = meanbl[sorted_trial_idxs]
		medianbl_sorted = medianbl[sorted_trial_idxs]
		meanloi_sorted = meanloi[sorted_trial_idxs]
		medianloi_sorted = medianloi[sorted_trial_idxs]
    end

    yss = Array{Float64}(undef, 0)
    xss = Array{Float64}(undef, 0)
    time2lick = Array{Float64}(undef, 0)
    pcInterval = Array{Float64}(undef, 0)
    tNos = Array{Int64}(undef, 0)
    lts = Array{Float64}(undef, 0)
    # baseline info
    # if includeBL_LOI
    	meanbls = Array{Float64}(undef, 0)
    	medianbls = Array{Float64}(undef, 0)
    	meanLOIs = Array{Float64}(undef, 0)
    	medianLOIs = Array{Float64}(undef, 0)
    # end
    #previous trial info
    ltm1s = Array{Float64}(undef, 0)
    rxnm1s = Array{Bool}(undef, 0)
    earlym1s = Array{Bool}(undef, 0)
    rewm1s = Array{Bool}(undef, 0)
    itim1s = Array{Bool}(undef, 0)
    # 2 trials back...
    ltm2s = Array{Float64}(undef, 0)
    rxnm2s = Array{Bool}(undef, 0)
    earlym2s = Array{Bool}(undef, 0)
    rewm2s = Array{Bool}(undef, 0)
    itim2s = Array{Bool}(undef, 0)
    #
    lickstate=Array{Bool}(undef, 0)
    #
    paths = Array{String}(undef, 0)
    sessioncodes = Array{String}(undef, 0)

    for i = 1:length(tNo_sorted)
        # get previous trial info
        if tNo_sorted[i] != 1 && tNo_sorted[1]==1
            if tNo_sorted[i] != 2
                # 2 trials back data...
                if tNo_sorted[i] - tNo_sorted[i-2] == 2
                    # we have 2 trials back!
                    ltm2 = lt_sorted[i-2]
                    outcomes2 = [false, false, false, false]
                    if lt_sorted_nonorm[i-2] < 0.7
                        outcomes2[1] = true
                    elseif lt_sorted_nonorm[i-2] >= 0.7 && lt_sorted_nonorm[i-2] < 3.333
                        outcomes2[2] = true
                    elseif lt_sorted_nonorm[i-2] >= 3.333 && lt_sorted_nonorm[i-2] < 7
                        outcomes2[3] = true
                    else
                        outcomes2[4] = true
                    end
                else
                    # no 2 trials back! we have no outcome info...
                    if normalize
                        ltm2 = 1.0
                    else
                        ltm2 = 17.1
                    end
                    outcomes2 = [false, false, false, false]
                end
            else
                # no 2 trials back! we have no outcome info...
                if normalize
                    ltm2 = 1.0
                else
                    ltm2 = 17.1
                end
                outcomes2 = [false, false, false, false]
            end

            # previous trial data...
            if tNo_sorted[i] - tNo_sorted[i-1] == 1
                # we have a prev trial!
                ltm1 = lt_sorted[i-1]
                outcomes = [false, false, false, false]
                if lt_sorted_nonorm[i-1] < 0.7
                    outcomes[1] = true
                elseif lt_sorted_nonorm[i-1] >= 0.7 && lt_sorted_nonorm[i-1] < 3.333
                    outcomes[2] = true
                elseif lt_sorted_nonorm[i-1] >= 3.333 && lt_sorted_nonorm[i-1] < 7
                    outcomes[3] = true
                else
                    outcomes[4] = true
                end
            elseif tNo_sorted[i] - tNo_sorted[i-1] == 2
                # no previous trial! But we do have 2 trials back...
                if normalize
                    ltm1 = 1.0
                else
                    ltm1 = 17.1
                end
                outcomes = [false, false, false, false]
                # we have 2 trials back!
                ltm2 = lt_sorted[i-1]
                outcomes2 = [false, false, false, false]
                if lt_sorted_nonorm[i-1] < 0.7
                    outcomes2[1] = true
                elseif lt_sorted_nonorm[i-1] >= 0.7 && lt_sorted_nonorm[i-1] < 3.333
                    outcomes2[2] = true
                elseif lt_sorted_nonorm[i-1] >= 3.333 && lt_sorted_nonorm[i-1] < 7
                    outcomes2[3] = true
                else
                    outcomes2[4] = true
                end
            else
                # no previous trial! we have no outcome info...
                if normalize
                    ltm1 = 1.0
                else
                    ltm1 = 17.1
                end
                outcomes = [false, false, false, false]
            end
        else
            # no previous trial! we have no outcome info...
            if normalize
                ltm1 = 1.0
            else
                ltm1 = 17.1
            end
            outcomes = [false, false, false, false]
            if normalize
                ltm2 = 1.0
            else
                ltm2 = 17.1
            end
            outcomes2 = [false, false, false, false]
        end

        for j = 1:length(ys_sorted[i])
            # extract each timepoint and signal value and relevant data
            push!(yss, ys_sorted[i][j])
            push!(xss, xs_sorted[i][j])
            if !normalize
                push!(time2lick, lt_sorted[i]-xs_sorted[i][j])
                if !isnan(xs_sorted[i][j]/lt_sorted[i])
                    push!(pcInterval, xs_sorted[i][j]/lt_sorted[i])
                else
                    push!(pcInterval, 1.)
                end
            else
                push!(time2lick, lt_sorted_nonorm[i]-xs_sorted[i][j])
                if !isnan(xs_sorted[i][j]/lt_sorted_nonorm[i])
                    push!(pcInterval, xs_sorted[i][j]/lt_sorted_nonorm[i])
                else
                    push!(pcInterval, 1.)
                end

            end
            push!(tNos, tNo_sorted[i])
            push!(lts, lt_sorted[i])
            push!(paths, path)
            push!(sessioncodes, seshcode)

            if includeBL_LOI
		    	push!(meanbls, meanbl_sorted[i])
		    	push!(medianbls, medianbl_sorted[i])
		    	push!(meanLOIs, meanloi_sorted[i])
		    	push!(medianLOIs, medianloi_sorted[i])
	    	else
	    		push!(meanbls, NaN)
		    	push!(medianbls, NaN)
		    	push!(meanLOIs, NaN)
		    	push!(medianLOIs, NaN)
	    	end

            push!(ltm1s, ltm1)
            push!(rxnm1s, outcomes[1])
            push!(earlym1s, outcomes[2])
            push!(rewm1s, outcomes[3])
            push!(itim1s, outcomes[4])

            push!(ltm2s, ltm2)
            push!(rxnm2s, outcomes2[1])
            push!(earlym2s, outcomes2[2])
            push!(rewm2s, outcomes2[3])
            push!(itim2s, outcomes2[4])
            if j==length(ys_sorted[i])
                push!(lickstate, true)
            else
                push!(lickstate, false)
            end
        end

    end

    # normalize the time2lick variable
    if normalize
	    time2lick = zero_to_one(time2lick)
    end

    y_shuffle = Random.shuffle(yss)

    df = DataFrame(DataID=1:length(xss),
        TrialNo = tNos,
        X=xss,
        Y=yss,
        LickTime=lts,
        Time2Lick=time2lick,
        PcInterval=pcInterval,
        LickTime_1back=ltm1s,
        LickTime_2back=ltm2s,
        Mean_Baseline = meanbls,
		Median_Baseline = medianbls,
		Mean_LOI = meanLOIs,
		Median_LOI = medianLOIs,
        Rxn_1back=rxnm1s,
        Early_1back=earlym1s,
        Reward_1back=rewm1s,
        ITI_1back=itim1s,
        Rxn_2back=rxnm2s,
        Early_2back=earlym2s,
        Reward_2back=rewm2s,
        ITI_2back=itim2s,
        LickState = lickstate,
        Yshuffle = y_shuffle,
        SessionCode=sessioncodes,
        Path=paths)
    return df
end;

function truncate_at_cue_and_lickp250(td::TrialData; cbuffer_s=0., lbuffer_s=0)
    newxs = []
    newys = []
    trialNo = []
    lickTime_s = []

    for i = 1:length(td.xdata)
        xx = td.xdata[i]
        cidx = findfirst(x->x>=0. +cbuffer_s, xx)
        yy = td.ydata[i]
        lidx = findlast(x->x<=(td.lickTime_s[i]+lbuffer_s), xx)
        if !isnothing(lidx)
	        xx = xx[1:lidx]
	        xx = xx[cidx:end]
	        yy = yy[1:lidx]
	        yy = yy[cidx:end]
	        push!(newxs, xx)
	        push!(newys, yy)
	        push!(trialNo, td.trialNo[i])
	        push!(lickTime_s, td.lickTime_s[i])
        end
    end
    trial_data = TrialData(newxs, newys, trialNo, lickTime_s, td.path, td.sessionCode)
    return trial_data
end;

function saveDataFrame(df, filename; path=nothing)
	if isnothing(path)
		path = pwd()
	end
	if Sys.iswindows()
		dlim="\" #"
	else
		dlim="/"
	end
	file = join([path, dlim, filename, timestamp_now(), ".csv"])
	println("saving... \n", file)
	CSV.write(file, df)
	return file
end

function importDataFrame(filename; suppresswarn=false)
	# NB! If you have nested DFs or arrays in the df, you need to instead use parse_imported_modelresult(filename)
	if !suppresswarn
		warning("NB! If you have nested DFs or arrays in the df, you need to instead use parse_imported_modelresult(filename)")
		println("Loading... \n", filename)
	end
	df = CSV.read(filename)
	return df
end

function unwrap_CSV_modelresult(filename)
	println("Loading... \n", filename)
	# empty_similar_df is a function to construct an empty dataframe

	# [eval(Meta.parse((f.train_ths[i][4:end]))) for i = 1:length(f.train_ths)]
	f = CSV.File(filename);
	colnames = propertynames(f[1,:][1])
	df_new = DataFrame()
	for i = 1:length(colnames)
	    # get the column data from CSV
	    coldata = [f[ii][i] for ii=1:length(f)]
	    if typeof(coldata[i]) == String
	        try
	            df_new[colnames[i]] = [eval(Meta.parse((coldata[j][4:end]))) for j = 1:length(coldata)] #[1,2]
	        catch
	#             println([eval(Meta.parse((coldata[j][4:end]))) for j = 1:length(coldata)])
	            df_new[colnames[i]] = coldata
	        end
	    else
       	 	df_new[colnames[i]] = coldata
	    end
    end
    return df_new
end


text:
using TerminalRegressionTests

const thisdir = dirname(@__FILE__)
TerminalRegressionTests.automated_test(
                joinpath(thisdir,"TRT.multiout"),
                ["Julia\n","Yes!!\n"]) do emuterm
    print(emuterm, "Please enter your name: ")
    name = strip(readline(emuterm))
    print(emuterm, "\nHello $name. Do you like tests? ")
    resp = strip(readline(emuterm))
    @assert resp == "Yes!!"
end


text:
using Plots, MAT

nprocs = (2, 2) # nprocs (x, y) dim

@views function vizme2D_mpi(nprocs)
    T  = []
    ip = 1
    for ipx = 1:nprocs[1]
        for ipy = 1:nprocs[2]
            file = matopen("T_$(ip-1).mat"); T_loc = read(file, "T"); close(file)
            nx_i, ny_i = size(T_loc,1)-2, size(T_loc,2)-2
            ix1, iy1   = 1+(ipx-1)*nx_i, 1+(ipy-1)*ny_i
            if (ip==1)  T = zeros(nprocs[1]*nx_i, nprocs[2]*ny_i)  end
            T[ix1:ix1+nx_i-1,iy1:iy1+ny_i-1] .= T_loc[2:end-1,2:end-1]
            ip += 1
        end
    end
    display(heatmap(T', framestyle=:box, aspect_ratio=1, xlims=(1, size(T,1)), ylims=(1, size(T,2)), c=:hot, title="heat diffusion"))
    return
end

vizme2D_mpi(nprocs)

max_stars_repo_name: luraess/geo-hpc-course

max_issues_repo_name: luraess/geo-hpc-course

max_forks_repo_name: luraess/geo-hpc-course


text:
"""
    Embedding{D, T}

## Fields
Assume `r` is a embedding, then

- **`r.points`**: The embedded state vectors, represented as an array of column
    vectors.
- **`r.embeddingdata`**: Information on what data and lags were used
    to construct the embedding.
"""
mutable struct Embedding{D, T} <: AbstractEmbedding{D, T}
    points::AbstractArray{T, 2}
    embeddingdata::EmbeddingData{D, T}
end

export Embedding


text:
using HDF5, JLD, Compat, Compat.Test, Compat.LinearAlgebra, LegacyStrings
using Compat: @warn

# Define variables of different types
x = 3.7
A = reshape(collect(1:15), 3, 5)
Aarray = Vector{Float64}[[1.2,1.3],[2.2,2.3,2.4]]
str = "Hello"
stringsA = String["It", "was", "a", "dark", "and", "stormy", "night"]
stringsU = String["It", "was", "a", "dark", "and", "stormy", "night"]
strings16 = convert(Array{UTF16String}, stringsA)
strings16_2d = reshape(strings16[1:6], (2,3))
empty_string = ""
empty_string_array = String[]
empty_array_of_strings = String[""]
tf = true
TF = A .> 10
B = [-1.5 sqrt(2) NaN 6;
     0.0  Inf eps() -Inf]
AB = Any[A, B]
t = (3, "cat")
c = ComplexF32(3,7)
cint = 1+im  # issue 108
C = Array(reinterpret(ComplexF64, vec(B)))
emptyA = zeros(0,2)
emptyB = zeros(2,0)
try
    global MyStruct
    mutable struct MyStruct
        len::Int
        data::Array{Float64}
        MyStruct(len::Int) = new(len)
        MyStruct(len::Int, data::Array{Float64}) = new(len, data)
    end
catch
end
ms = MyStruct(2, [3.2, -1.7])
msempty = MyStruct(5, Float64[])
sym = :TestSymbol
syms = [:a, :b]
d = Dict([(syms[1],"aardvark"), (syms[2], "banana")])
ex = quote
    function incrementby1(x::Int)
        x+1
    end
end
T = UInt8
char = 'x'
unicode_char = '\U10ffff'
α = 22
β = Any[[1, 2], [3, 4]]  # issue #93
vv = Vector{Int}[[1,2,3]]  # issue #123
typevar = Array{Int}[[1]]
typevar_lb = Vector{<:Integer}[[1]]
typevar_ub = (Vector{U} where U>:Int)[[1]]
typevar_lb_ub = (Vector{U} where Int<:U<:Real)[[1]]
# Immutable type:
rng = 1:5
# Type with a pointer field (#84)
struct ObjWithPointer
    a::Ptr{Cvoid}
end
objwithpointer = ObjWithPointer(0)
# Custom BitsType (#99)
primitive type MyBT 64 end
bt = reinterpret(MyBT, Int64(55))
# Symbol arrays (#100)
sa_asc = [:a, :b]
sa_utf8 = [:α, :β]
# SubArray (to test tuple type params)
subarray = view([1:5;], 1:5)
# Array of empty tuples (to test tuple type params)
arr_empty_tuple = Tuple{}[]
struct EmptyImmutable end
emptyimmutable = EmptyImmutable()
arr_emptyimmutable = [emptyimmutable]
mutable struct EmptyType end
emptytype = EmptyType()
arr_emptytype = [emptytype]
struct EmptyII
    x::EmptyImmutable
end
emptyii = EmptyII(EmptyImmutable())
struct EmptyIT
    x::EmptyType
end
emptyit = EmptyIT(EmptyType())
mutable struct EmptyTI
    x::EmptyImmutable
end
emptyti = EmptyTI(EmptyImmutable())
mutable struct EmptyTT
    x::EmptyType
end
emptytt = EmptyTT(EmptyType())
struct EmptyIIOtherField
    x::EmptyImmutable
    y::Float64
end
emptyiiotherfield = EmptyIIOtherField(EmptyImmutable(), 5.0)

# Unicode type field names (#118)
mutable struct MyUnicodeStruct☺{τ}
    α::τ
    ∂ₓα::τ
    MyUnicodeStruct☺{τ}(α::τ, ∂ₓα::τ) where {τ} = new{τ}(α, ∂ₓα)
end
unicodestruct☺ = MyUnicodeStruct☺{Float64}(1.0, -1.0)
# Arrays of matrices (#131)
array_of_matrices = Matrix{Int}[[1 2; 3 4], [5 6; 7 8]]
# Tuple of arrays and bitstype
tup = (1, 2, [1, 2], [1 2; 3 4], bt)
# Empty tuple
empty_tup = ()
# Non-pointer-free immutable
struct MyImmutable{T}
    x::Int
    y::Vector{T}
    z::Bool
end
nonpointerfree_immutable_1 = MyImmutable(1, [1., 2., 3.], false)
nonpointerfree_immutable_2 = MyImmutable(2, Any[3., 4., 5.], true)
struct MyImmutable2
    x::Vector{Int}
    MyImmutable2() = new()
end
nonpointerfree_immutable_3 = MyImmutable2()
# Immutable with a non-concrete datatype (issue #143)
struct Vague
    name::String
end
vague = Vague("foo")
# Immutable with a union of BitsTypes
struct BitsUnion
    x::Union{Int64, Float64}
end
bitsunion = BitsUnion(5.0)
# Immutable with a union of Types
struct TypeUnionField
    x::Union{Type{Int64}, Type{Float64}}
end
typeunionfield = TypeUnionField(Int64)
# Generic union type field
struct GenericUnionField
    x::Union{Vector{Int},Int}
end
genericunionfield = GenericUnionField(1)
# Array references
arr_contained = [1, 2, 3]
arr_ref = typeof(arr_contained)[]
push!(arr_ref, arr_contained, arr_contained)
# Object references
mutable struct ObjRefType
    x::ObjRefType
    y::ObjRefType
    ObjRefType() = new()
    ObjRefType(x, y) = new(x, y)
end
ref1 = ObjRefType()
obj_ref = ObjRefType(ObjRefType(ref1, ref1), ObjRefType(ref1, ref1))
# Immutable that requires padding between elements in array
struct PaddingTest
    x::Int64
    y::Int8
end
padding_test = PaddingTest[PaddingTest(i, i) for i = 1:8]
# Empty arrays of various types and sizes
empty_arr_1 = Int[]
empty_arr_2 = Array{Int}(undef, 56, 0)
empty_arr_3 = Any[]
empty_arr_4 = Array{Any}(undef, 0, 97)
# Moderately big dataset (which will be mmapped)
bigdata = [1:10000;]
# BigFloats and BigInts
bigints = big(3).^(1:100)
bigfloats = big(3.2).^(1:100)
# None
none = Union{}
nonearr = Array{Union{}}(undef, 5)
# nothing
scalar_nothing = nothing
vector_nothing = Union{Int,Nothing}[1,nothing]

# some data big enough to ensure that compression is used:
Abig = kron(Matrix(1.0I, 10, 10), reshape(1:400, 20, 20))
Bbig = Any[i for i=1:3000]
Sbig = "A test string "^1000

# Bitstype type parameters
mutable struct BitsParams{x}; end
bitsparamfloat  = BitsParams{1.0}()
bitsparambool   = BitsParams{true}()
bitsparamsymbol = BitsParams{:x}()
bitsparamint    = BitsParams{1}()
bitsparamuint   = BitsParams{0x01}()
bitsparamint16  = BitsParams{Int16(1)}()

# Tuple of tuples
tuple_of_tuples = (1, 2, (3, 4, [5, 6]), [7, 8])

iseq(x,y) = isequal(x,y)
iseq(x::MyStruct, y::MyStruct) = (x.len == y.len && x.data == y.data)
iseq(x::MyImmutable, y::MyImmutable) = (isequal(x.x, y.x) && isequal(x.y, y.y) && isequal(x.z, y.z))
@static if VERSION ≥ v"0.7.0-DEV.3693" # empty mutable structs are no longer singletons
    iseq(x::EmptyType, y::EmptyType) = true
    iseq(x::EmptyIT, y::EmptyIT) = true
    iseq(x::Array{EmptyType}, y::Array{EmptyType}) = size(x) == size(y)
    iseq(x::BitsParams{T}, y::BitsParams{T}) where {T} = true
    iseq(x::BitsParams, y::BitsParams) = false
end
iseq(x::Union{EmptyTI, EmptyTT}, y::Union{EmptyTI, EmptyTT}) = iseq(x.x, y.x)
iseq(c1::Array{Base.Sys.CPUinfo}, c2::Array{Base.Sys.CPUinfo}) = length(c1) == length(c2) && all([iseq(c1[i], c2[i]) for i = 1:length(c1)])
function iseq(c1::Base.Sys.CPUinfo, c2::Base.Sys.CPUinfo)
    for n in fieldnames(Base.Sys.CPUinfo)
        if getfield(c1, n) != getfield(c2, n)
            return false
        end
    end
    true
end
iseq(x::MyUnicodeStruct☺, y::MyUnicodeStruct☺) = (x.α == y.α && x.∂ₓα == y.∂ₓα)
iseq(x::Array{Union{}}, y::Array{Union{}}) = size(x) == size(y)
macro check(fid, sym)
    ex = quote
        let tmp
            try
                tmp = read($fid, $(string(sym)))
            catch e
                @warn("Error reading ", $(string(sym)))
                rethrow(e)
            end
            if !iseq(tmp, $sym)
                written = $sym
                error("For ", $(string(sym)), ", read value $tmp does not agree with written value $written")
            end
            written_type = typeof($sym)
            if typeof(tmp) != written_type
                error("For ", $(string(sym)), ", read type $(typeof(tmp)) does not agree with written type $(written_type)")
            end
        end
    end
    esc(ex)
end

# Test for equality of expressions, skipping line numbers
checkexpr(a, b) = @assert a == b
function checkexpr(a::Expr, b::Expr)
    @assert a.head == b.head
    i = 1
    j = 1
    while i <= length(a.args) && j <= length(b.args)
        if (isa(a.args[i], Expr) && a.args[i].head == :line) || isa(a.args[i], LineNumberNode)
            i += 1
            continue
        end
        if (isa(b.args[j], Expr) && b.args[j].head == :line) || isa(b.args[j], LineNumberNode)
            j += 1
            continue
        end
        checkexpr(a.args[i], b.args[j])
        i += 1
        j += 1
    end
    @assert i >= length(a.args) && j >= length(b.args)
end

for fname in ("v0.4.13-julia-0.3.jld", "v0.4.14-julia-0.4.0-dev+4483.jld")
    undefv = Array{Any}(undef, 1)
    undefm = Array{Any}(undef, 2, 2)
    ms_undef = MyStruct(0)

    fidr = jldopen(joinpath(dirname(@__FILE__), fname), "r")
    @check fidr x
    @check fidr A
    dsetA = fidr["A"]
    @test ndims(dsetA) == ndims(A)
    @test size(dsetA) == size(A)
    @test size(dsetA, 1) == size(A, 1)
    @test size(dsetA, 2) == size(A, 2)
    @test size(dsetA, 3) == size(A, 3)
    @check fidr Aarray
    @check fidr str
    @check fidr stringsA
    @check fidr stringsU
    @check fidr strings16
    @check fidr strings16_2d
    @check fidr empty_string
    @check fidr empty_string_array
    @check fidr empty_array_of_strings
    @check fidr tf
    @check fidr TF
    @check fidr AB
    @check fidr t
    @check fidr c
    @check fidr cint
    @check fidr C
    @check fidr emptyA
    @check fidr emptyB
    @check fidr ms
    @check fidr msempty
    @check fidr sym
    @check fidr syms
    @check fidr d
    exr = read(fidr, "ex")   # line numbers are stripped, don't expect equality
    checkexpr(ex, exr)
    @check fidr T
    @check fidr char
    @check fidr unicode_char
    @check fidr α
    @check fidr β
    @check fidr vv
    @check fidr rng
    @check fidr typevar
    @check fidr typevar_lb
    @check fidr typevar_ub
    @check fidr typevar_lb_ub

    # Special cases for reading undefs
    undefv = read(fidr, "undef")
    if !isa(undefv, Array{Any, 1}) || length(undefv) != 1 || isassigned(undefv, 1)
        error("For undef, read value does not agree with written value")
    end
    undefm = read(fidr, "undefs")
    if !isa(undefm, Array{Any, 2}) || length(undefm) != 4 || any(map(i->isassigned(undefm, i), 1:4))
        error("For undefs, read value does not agree with written value")
    end
    ms_undef = read(fidr, "ms_undef")
    if !isa(ms_undef, MyStruct) || ms_undef.len != 0 || isdefined(ms_undef, :data)
        error("For ms_undef, read value does not agree with written value")
    end

    @check fidr bt
    @check fidr sa_asc
    @check fidr sa_utf8
    # @check fidr subarray
    @check fidr arr_empty_tuple
    @check fidr emptyimmutable
    @check fidr emptytype
    @check fidr arr_emptyimmutable
    @check fidr arr_emptytype
    @check fidr emptyii
    @check fidr emptyit
    @check fidr emptyti
    @check fidr emptytt
    @check fidr emptyiiotherfield
    @check fidr unicodestruct☺
    @check fidr array_of_matrices
    @check fidr tup
    @check fidr empty_tup
    @check fidr nonpointerfree_immutable_1
    @check fidr nonpointerfree_immutable_2
    @check fidr nonpointerfree_immutable_3
    vaguer = read(fidr, "vague")
    @test typeof(vaguer) == typeof(vague) && vaguer.name == vague.name
    @check fidr bitsunion
    @check fidr typeunionfield
    @check fidr genericunionfield

    arr = read(fidr, "arr_ref")
    @test arr == arr_ref
    @test arr[1] === arr[2]

    obj = read(fidr, "obj_ref")
    @test obj.x.x === obj.x.y == obj.y.x === obj.y.y
    @test obj.x !== obj.y

    @check fidr padding_test
    @check fidr empty_arr_1
    @check fidr empty_arr_2
    @check fidr empty_arr_3
    @check fidr empty_arr_4
    @check fidr bigdata
    @check fidr bigfloats
    @check fidr bigints
    @check fidr none
    @check fidr nonearr
    @check fidr scalar_nothing
    @check fidr vector_nothing
    @check fidr Abig
    @check fidr Bbig
    @check fidr Sbig
    @check fidr bitsparamfloat
    @check fidr bitsparambool
    @check fidr bitsparamsymbol
    @check fidr bitsparamint
    @check fidr bitsparamuint
    @check fidr tuple_of_tuples

    x1 = read(fidr, "group1/x")
    @assert x1 == Any[1]
    x2 = read(fidr, "group2/x")
    @assert x2 == Any[2]

    close(fidr)
end


text:
using FermiCG
using Printf
using Test
using LinearAlgebra
using TensorOperations

function test1()
    g1 = rand(4,10) .- .5
    g2 = rand(4,30) .- .5

    c = .1
    @tensor V[I,J] := g1[i,I] * g2[i,J] *c

    b = FermiCG.upper_bound(g1,g2,c=c)
    maximum(abs.(V)) <= b || @printf(" %12.8f %12.8f\n", maximum(abs.(V)), b)
    return maximum(abs.(V)) <= b
end

function test2()
    v = rand(2,4)
    g1 = rand(2,10)
    g2 = rand(4,30)

    c = .1
    @tensor V[I,J] := v[i,j] * g1[i,I] * g2[j,J] * c

    b = FermiCG.upper_bound(v,g1,g2,c=c)
    maximum(abs.(V)) <= b || @printf(" %12.8f %12.8f\n", maximum(abs.(V)), b)
    return maximum(abs.(V)) <= b
end

function test3()
    v = rand(2,4,2)
    g1 = rand(2,10)
    g2 = rand(4,30)
    g3 = rand(2,20)

    c = .1
    @tensor V[I,J,K] := v[i,j,k] * g1[i,I] * g2[j,J] * g3[k,K] * c

    b = FermiCG.upper_bound(v,g1,g2,g3,c=c)
    maximum(abs.(V)) <= b || @printf(" %12.8f %12.8f\n", maximum(abs.(V)), b)
    return maximum(abs.(V)) <= b
end

function test4()
    v = rand(2,4,2,3)
    g1 = rand(2,10).- .5
    g2 = rand(4,30).- .5
    g3 = rand(2,20).- .5
    g4 = rand(3,9) .- .5

    c = .1
    @tensor V[I,J,K,L] := v[i,j,k,l] * g1[i,I] * g2[j,J] * g3[k,K] * g4[l,L]* c

    b = FermiCG.upper_bound(v,g1,g2,g3,g4,c=c)
    maximum(abs.(V)) <= b || @printf(" %12.8f %12.8f\n", maximum(abs.(V)), b)
    return maximum(abs.(V)) <= b
end

function test4a()
    v = rand(2,4,2,3) .- .5
    g1 = rand(2,10) .- .5
    g2 = rand(4,30) .- .5
    g3 = rand(2,20) .- .5
    g4 = rand(3,9)  .- .5

    c = .1
    @tensor V[I,J,K,L] := v[i,j,k,l] * g1[i,I] * g2[j,J] * g3[k,K] * g4[l,L] * c

    thresh = 1.0
    new1, new2, new3, new4 = FermiCG.upper_bound2(v,g1,g2,g3,g4,thresh, c=c)
    g1a = g1[:,new1]
    g2a = g2[:,new2]
    g3a = g3[:,new3]
    g4a = g4[:,new4]

    l1 = count(V.>thresh)
    l2 = 0
    if minimum(length.([new1,new2,new3,new4])) > 0
        @tensor V2[I,J,K,L] := v[i,j,k,l] * g1a[i,I] * g2a[j,J] * g3a[k,K] * g4a[l,L] * c
        l2 = count(V2.>thresh)
    end
    return l1 == l2
end

function test3a()
    v = rand(2,4,2) .- .5
    g1 = rand(2,10) .- .5
    g2 = rand(4,30) .- .5
    g3 = rand(2,20) .- .5

    c = .1
    @tensor V[I,J,K] := v[i,j,k] * g1[i,I] * g2[j,J] * g3[k,K] * c

    thresh = 1.0
    new1, new2, new3 = FermiCG.upper_bound2(v,g1,g2,g3,thresh, c=c)
    g1a = g1[:,new1]
    g2a = g2[:,new2]
    g3a = g3[:,new3]

    l1 = count(V.>thresh)
    l2 = 0
    if minimum(length.([new1,new2,new3])) > 0
        @tensor V2[I,J,K] := v[i,j,k] * g1a[i,I] * g2a[j,J] * g3a[k,K] * c
        l2 = count(V2.>thresh)
    end
    return l1 == l2
end

function tuck_bound()
    v = FermiCG.recompose(FermiCG.Tucker(rand(6,7,8,9).-.5,max_number=4)) .+ .01*(rand(6,7,8,9) .- .5)
    v = v./norm(v)
    tuck = FermiCG.Tucker_tot(v, thresh=.1,verbose=1)
    err1 = norm(v-FermiCG.recompose(tuck)) <= .1

    v = (rand(6,7,8,9) .- .5)
    v = v./norm(v)
    tuck = FermiCG.Tucker_tot(v, thresh=.4,verbose=1)
    err2 = norm(v-FermiCG.recompose(tuck)) <= .4
    return err1 && err2
end

@testset "bounds" begin
    for i in 1:10
        @test tuck_bound()
        @test test1()
        @test test2()
        @test test3()
        @test test4()
        @test test3a()
        @test test4a()
    end
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
struct Cluster <: AbstractProcessingStep
    table::Observable{SimpleTable}
    card::ProcessingCard
end

# TODO: add indirection layer with distance support, as in dimensionality reduction?
const clusterings = (
    kmeans=kmeans,
    mcl=mcl,
    affinityprop=affinityprop,
)

function Cluster(table::Observable{SimpleTable})
    analysis_names = collect(map(string, keys(clusterings)))

    analysis_options = vecmap(analysis_names) do name
        if name == "kmeans"
            name * " classes" => [string(i) for i in 1:100]
        else
            name => String[]
        end
    end

    default_names = ":cluster"

    wdgs = (
        inputs=RichTextField("Inputs", data_options(table, keywords=["", "weights"]), ""),
        outputs=RichTextField("Outputs", ["" => ["cluster"]], default_names),
        method=RichTextField("Method", analysis_options, ""),
    )

    card = ProcessingCard(:Cluster; wdgs...)
    return Cluster(table, card)
end

function (cluster::Cluster)(data)
    card = cluster.card
    args, kwargs = extract_all_arguments(card.inputs)
    name = extract_positional_argument(card.outputs)

    dist = Euclidean() # TODO: make configurable
    cols = Tables.getcolumn.(Ref(data), Symbol.(args))
    X = reduce(vcat, transpose.(cols))
    D = pairwise(dist, X, dims=2)

    f = extract_function(card.method)
    an = clusterings[Symbol(f)]
    input = an === kmeans ? X : D
    positional, named = [], Pair{Symbol, Any}[Symbol(k) => Tables.getcolumn(data, Symbol(v)) for (k, v) in kwargs]
    for (k, v) in extract_named_arguments(card.method)
        if an === kmeans && k == "classes"
            push!(positional, parse(Int, v))
        else
            push!(named, Symbol(k) => v)
        end
    end
    anres = an(input, positional...; named...)
    return SimpleTable(Symbol(name) => map(nonnumeric, Clustering.assignments(anres)))
end


text:
"""
    Molecule

Photophysical properties of a molecule.

This is the most general type of luminecent or scattering single molecule.
Inherited types will defines the properties of a class of molecules.

"""
abstract type Molecule end

"""
    GenericFluor

Defines a fluorophore

# Fields
- γ: photon emission rate in Hz, Default: 1e3
- q: state transision matrix. Default: q=[1.0]
"""
mutable struct GenericFluor <: Molecule
    γ::AbstractFloat
    q::Array{<:AbstractFloat}
end
function GenericFluor(;
    γ::AbstractFloat=1e5,
    q::Array{<:AbstractFloat}=[1]
    )
    return GenericFluor(γ,q)
end


text:
# Submodule with some operations that will help to pad multidimensional arrays and Upsample in the Shearlab3D.m fashion

###################################################
# Padding array function to resize an array padding
# zeros
"""
...
padarray(array, newSize,gpu)
pads array with zeros to have a new size if the input array
is bigger than the targeted size, it will centered it at zero
and cut it off to fit
...
"""
function padarray(array::AbstractArray{T},newSize,gpu = false) where {T<:Number}
    # If you want to use gpu via ArrayFire array needs to be ArrayFire.AFArray{Float32,2}
    # Small patch if the array is one dimensional
    if size([size(array)...])[1]==1
        array = transpose(array)
    end
    padSizes = zeros(Integer,1,length(newSize))
    n = length(newSize)
    sizes = size(array)
    @fastmath @inbounds @simd for k = 1:n
        sizeDiff = newSize[k]-sizes[k];
        if mod(sizeDiff,2) == 0
            padSizes[k] = Int(sizeDiff/2);
        else
            padSizes[k] = Int(ceil.(sizeDiff/2))
        end
    end
    # Correct the matrix that is bigger than the targeted size
    idbig = [1:sizes[i] for i in 1:n]
    @fastmath @inbounds @simd for i in 1:length(newSize)
        if padSizes[i] < 0
            idbig[i] = (round(Int64,(size(array,i)-newSize[i])/2)+1):(round(Int64,(size(array,i)+newSize[i])/2))
            padSizes[i] = 0
        end
    end
    # We need to check if some of the padsizes are negative and cut the array from the center
    # Initialize the padded array with zeros
    if gpu
        paddedArray = AFArray(zeros(typeof(array[1]), newSize...))
    else
        paddedArray = zeros(typeof(array[1]),newSize...)
    end
    # lets create the indices array
    if gpu
        paddedArray[[idbig[1]+padSizes[1],idbig[2]+padSizes[2]]...] = array[idbig...]
    else
        view(paddedArray,[idbig[1].+padSizes[1],idbig[2].+padSizes[2]]...) .= view(array,idbig...)
    end
    return paddedArray
end #padarray

################################################################
# Function that flips from left to right an array in the second dimension
"""
...
fliplr(array,gpu) flips an array from left to right in the second dimension
...
"""
function fliplr(array::AbstractArray{T},gpu = false) where {T<:Number}
    return reverse(array,dims=2)
end #fliplr

##################################################################
# function that upsample an multidimensional array based on the same
# function at the matlab version
"""
..
upsample(array,dims,nZeros,gpu) upsample an array, in the dimensions dims
with nZeros
...
"""
function upsample(array::AbstractArray{T},dims::Integer,nZeros::Integer, gpu = false) where {T<:Number}
    sz = [size(array)...]
    szUpsampled = sz
    szUpsampled[dims] = (szUpsampled[dims]-1).*(nZeros+1)+1
    if gpu
        arrayUpsampled = AFArray(zeros(typeof(array[1]),szUpsampled...))
    else
        arrayUpsampled = zeros(typeof(array[1]),szUpsampled...);
    end
    # Generate the indices per dimension
    ids = [1:1:size(array,i) for i in 1:length(size(array))]
    @fastmath @inbounds @simd  for k in 1:length(dims)
        ids[dims[k]] = 1:(nZeros[k]+1):szUpsampled[dims[k]]
    end
    arrayUpsampled[ids...] = array
    return arrayUpsampled
end #upsample

#####################################################################
# Function that rounds a number to the nearest integer towards zero
"""
...
fix(x) rounds a number x to the nearest integer towards zero
...
"""
function fix(x::T) where {T<:Number}
    if x < 0
        fixed = ceil(x)
    else
        fixed = floor(x)
    end
    return convert(Int64,fixed)
end #fix

#######################################################################
# Function that shears an array in order k in the direction of axis
# based on the same function on the Matlab version
"""
...
dshear(array,k,axis,gpu) shears and array in order k in the direction of
axis
...
"""
function dshear(array::AbstractArray{T},k::Integer,axis::Integer, gpu = false) where {T<:Number}
    if gpu
        array = Array(array)
    end
    if k == 0
        sheared = array;
    else
        rows = size(array,1);
        cols = size(array,2);
        sheared = zeros(typeof(array[1]),size(array)...)
        if axis == 1
            for col = 1:cols
                sheared[:,col] = reshape(circshift(reshape(array[:,col],rows,1),[-k*(floor(cols/2)+1-col) 0]),(rows,));
            end
        else
            for row = 1:rows
                sheared[row,:] = reshape(circshift(reshape(array[row,:],1,cols),[0 -k*(floor(rows/2)+1-row)]),(cols,));
            end
        end
    end
    if gpu
        sheared = AFArray(sheared)
    end
    return sheared
end #dshear

######################################################################
## Type of filter configurations
struct Filterconfigs
    directionalFilter
    scalingFilter
    waveletFilter
    scalingFilter2
end #Filterconfigs

#######################################################################
# Function that check the sizes of the filters to know if it is possible
"""
...
checkfiltersizes(rows,cols,shearLevels,directionalFilter,scalingFilter,waveletFilter,scalingFilter2) function that check
wBand
the size of the filters and set new possible configurations
...
"""
function checkfiltersizes(rows,cols,shearLevels,directionalFilter,scalingFilter,waveletFilter,scalingFilter2)
    # println("input is rows=$(rows)")
    # println("cols=$(cols)")
    # println("shearLevels=$(shearLevels)")
    # println("directionalFilter=$(directionalFilter)")
    # println("scalingFilter=$(scalingFilter)")
    # println("waveletFilter=$(waveletFilter)")
    # println("scalingFilter2 = $(scalingFilter2)")
    # Lets initialize the FilterConfig array,
    filterSetup = []
    # Set all configurations

    # Configuration 1
    push!(filterSetup,Filterconfigs(directionalFilter,scalingFilter,
                                    waveletFilter,scalingFilter2));

    # Configuration 2
    # Check the default configuration
    scalingFilter = filt_gen("scaling_shearlet");
    directionalFilter = filt_gen("directional_shearlet");
    waveletFilter = mirror(scalingFilter);
    scalingFilter2 =  scalingFilter;
    push!(filterSetup,Filterconfigs(directionalFilter,scalingFilter,
                                    waveletFilter,scalingFilter2));

    # Configuration 3
    # Just change the directionalFilter
    directionalFilter = filt_gen("directional_shearlet2");
    push!(filterSetup,Filterconfigs(directionalFilter,scalingFilter,
                                    waveletFilter,scalingFilter2));

    # Configuration 4
    # The same as 3
    push!(filterSetup,Filterconfigs(directionalFilter,scalingFilter,
                                    waveletFilter,scalingFilter2));

    # Configuration 5
    # It changes the scalingFilter, the waveletFilter and the scalingFilter2
    scalingFilter = filt_gen("Coiflet1");
    waveletFilter = mirror(scalingFilter);
    scalingFilter2 =  scalingFilter;
    push!(filterSetup,Filterconfigs(directionalFilter,scalingFilter,
                                    waveletFilter,scalingFilter2));

    # Configuration 6
    # It changes the scalingFilter, the waveletFilter and the scalingFilter2
    scalingFilter = filt_gen(WT.db2)[2:5];
    waveletFilter = mirror(scalingFilter);
    scalingFilter2 =  scalingFilter;
    push!(filterSetup,Filterconfigs(directionalFilter,scalingFilter,
                                    waveletFilter,scalingFilter2));

    # Configuration 7
    # Just change the directionalFilter
    directionalFilter = filt_gen("directional_shearlet3");
    push!(filterSetup,Filterconfigs(directionalFilter,scalingFilter,
                                    waveletFilter,scalingFilter2));

    # Configuration 8
    # It changes the scalingFilter, the waveletFilter and the scalingFilter2
    scalingFilter = filt_gen(WT.haar)[2:3]
    waveletFilter = mirror(scalingFilter);
    scalingFilter2 =  scalingFilter;
    push!(filterSetup,Filterconfigs(directionalFilter,scalingFilter,
                                    waveletFilter,scalingFilter2));

    # Check the sizes of the filters in comparison with the rows and cols
    kk = 0;
    success1 = 0;
    # println("lets find out why this breaks")
    for k = 1:8
        # println("k=$(k)")
        ## check 1
        lwfilter = length(filterSetup[k].waveletFilter);
        lsfilter = length(filterSetup[k].scalingFilter);
        lcheck1 = lwfilter;
        for j = 1:(length(shearLevels)-1)
            lcheck1 = lsfilter + 2*lcheck1 - 2;
        end
        # println("lcheck1 = $(lcheck1)")
        if lcheck1 > cols || lcheck1 > rows
            continue;
        end

        ## check  2
        rowsdirfilter = size(filterSetup[k].directionalFilter,1);
        colsdirfilter = size(filterSetup[k].directionalFilter,2);
        lcheck2 = (rowsdirfilter-1)*2^(maximum(shearLevels)+1) + 1;

        lsfilter2 = length(filterSetup[k].scalingFilter2);
        lcheck2help = lsfilter2;
        for j = 1:maximum(shearLevels)
            lcheck2help = lsfilter2 + 2*lcheck2help - 2;
        end
        lcheck2 = lcheck2help + lcheck2 - 1;
        # println("lcheck2 = $(lcheck2)")
        if lcheck2 > cols || lcheck2 > rows || colsdirfilter > cols || colsdirfilter > rows
            continue;
        end
        success1 = 1;
	kk = k;
        break;
    end
    if success1 == 0
        @error "The specified Shearlet system was not available for data of size "* string(rows) *"x",string(cols)* ". Filters were not set (see operations.jl)."
    end
    if success1 == 1 && kk > 1
        @warn "The specified Shearlet system was not available for data of size "*string(rows)*"x"*string(cols)*". Filters were automatically set to configuration "*string(kk)*" (see operations.jl)."
    end
    filterSetup[kk]
end #checkfiltersizes

function describeConfig(kk::Int)
    if kk==1
        # the requested settings are available
        return "you shouldn't see this"
    elseif kk==2
        # Configuration 2
        return "the scaling shearlet as both scaling function and  and and directional shearlet"
    elseif kk==3
        # Configuration 3
        return "the scaling shearlet and the 2nd directional shearlet"
    elseif kk==4
        # Configuration 4
        return "the scaling shearlet and the 2nd directional shearlet"
    elseif kk==5
        # Configuration 5
        return "the Coiflet1 scaling function, the mirrored scaling shearlet as the wavelet, and the 2nd directional shearlet"
    elseif kk==6
        # Configuration 6
        return "the Daubechies2 scaling function, the mirrored version as the wavelet, and the 2nd directional shearlet"
    elseif kk==7
        # Configuration 7
        return "the Daubechies2 scaling function, the mirrored version as the wavelet, and the 3nd directional shearlet"
    elseif kk==8
        # Configuration 8
        return "the Haar scaling function, the mirrored version as the wavelet, and the 3nd directional shearlet"
    end
end

"""
    findSupport(shears::Array{<:Number, 3}; tolerance::Float64=1e-10)

The default tolerance keeps ~99% of the coefficients (and a significantly high fraction of the mass)
TODO: implement this in an efficient way (i.e. use the fact that these originate from sheared versions of the same functions)
"""
function findSupport(shears::Array{<:Number, 3}; tolerance::Float64=1e-10)
  nShears = size(shears,3)
  effSupport = Array{Tuple{Tuple{Int,Int}, Tuple{Int,Int}}}(undef, nShears)
  spaceShears = real.(fftshift(ifft(ifftshift(shears), (1, 2))))
  actuallySupported = abs.(spaceShears).>=tolerance
  doesThisRowHaveTolEls= maximum(actuallySupported,dims=1)
  doesThisColHaveTolEls = maximum(actuallySupported,dims=2)
  size(doesThisRowHaveTolEls)
  size(doesThisColHaveTolEls)
  for i=1:nShears
     effSupport[i] = ((findfirst(doesThisColHaveTolEls[:,1,i]), findlast(doesThisColHaveTolEls[:,1,i])), (findfirst(doesThisRowHaveTolEls[1,:,i]), findlast(doesThisRowHaveTolEls[1,:,i])))
  end
  return effSupport
end

"""
    padBy =

getPadBy(supports; tolerance::Float64=1e-10)

Find the amount to pad by to guarantee that no shearlet's support wraps around
the boundary. Support here being the space window that captures `1-tolerance`
of the mass

"""
function getPadBy(shears; tolerance::Float64=1e-10)
    supports = findSupport(shears; tolerance = tolerance)
    padBy = (0, 0)

    for sysPad in supports
        padBy = (max(sysPad[1][2]-sysPad[1][1], padBy[1]), max(padBy[2],
                                                               sysPad[2][2]-sysPad[2][1]))
    end
    return padBy
end

"""
    pad(x, padBy)
create a padded version of X so that it has +/-padBy[1] extra zeros in the
first dimension and +/-padBy[2] in the second dimension. All other dimensions
are untouched.
"""
function pad(x, padBy)
    T = eltype(x)
    szx = size(x)
    corner = zeros(T, padBy...,  szx[3:end]...)
    firstRow = cat(corner,
                   zeros(T, padBy[1], szx[2:end]...),
                   corner, dims = 2)
    secondRow = cat(zeros(T, szx[1] , padBy[2], szx[3:end]...),
                    x,
                    zeros(T, szx[1] , padBy[2], szx[3:end]...),
                    dims=2)
    thirdRow = cat(corner,
                   zeros(T, padBy[1], szx[2:end]...),
                   corner,
                   dims = 2)
    return cat(firstRow, secondRow, thirdRow, dims = 1)
end


text:
# This file is a part of EasyPkg, licensed under the MIT License (MIT).

using Base.Test

@testset "This should fail" begin
	@test 1 == 1
	@test 2 == 0
    @test 3 == 3
    error("Some bogus error")
end


text:
# https://github.com/vtjnash/Glob.jl/issues/11
rglob(pat, topdir) = Base.Iterators.flatten(map(d -> glob(pat, d[1]), walkdir(topdir)))

function parse_obsxlsx(obs_path::String, input_dir::String)
    ws = readxlsheet(string(obs_path), "2017년")

    re_opened = r"신규"
    re_closed = r"폐쇄"
    re_moved = r"이전"
    re_hangul = r"[가-힣]+"
    re_opened_date = r"(\d+)\.\W*(\d+)\.?"
    #re_closed_date = r"(\d+)\.\W*(\d+)\.?\W*(\d+)\."
    re_closed_date = r"(\d+)\.\W*(\d+)\.?(\W*(\d+)\.?)*"
    KST = tz"Asia/Seoul"
    df = DataFrame(stationCode=[], stationName=[], oDate=[], cDate=[], lon=[], lat=[])
    cnt = 0

    #=
        row[1] = state (시, 도)
        row[2] = city (도시)
        row[3] = stationCode (측정소코드)
        row[4] = stationName (측정소명)
        row[5] = address (주소)
        row[6] = langitude (경도)
        row[7] = latitude (위도)
        row[8] = remakrs (비고)
    =#

    for i in 1:size(ws, 1)
        # just pick old date
        opened_date = ZonedDateTime(1970, 1, 1, tz"Asia/Seoul")
        # close value to maximum value in Int64
        closed_date = ZonedDateTime(2037, 12, 31, tz"Asia/Seoul")

        row = ws[i, :]
        # prev_row = ws[i - 1, :]
        # if stationCode is blank, skip, so no move
        if isa(row[3], DataValue) && isna(row[3])
            continue
        end

        # if there is not stationCode (i.e. wrong column), skip loop
        try
            if isa(row[3], String)
                station_code = parse(Int, row[3])
            elseif isa(row[3], AbstractFloat)
                station_code = Int(row[3])
            end
        catch e
            if isa(e, ArgumentError)
                continue
            end
        end

        # check lat / long
        if isa(row[6], Real) && isa(row[7], Real)
            station_code = 0
            station_name = ""
            long_X = float(row[6])
            lat_Y = float(row[7])

            cnt += 1

            # parse station code and name, doesn't support moved case
            if isa(row[3], String)
                station_code = parse(Int, row[3])
            elseif isa(row[3], AbstractFloat)
                station_code = Int(row[3])
            end

            if occursin(re_hangul, string(row[4]))
                station_name = string(row[4])
            end

            # parse date of station opened
            # check there is a stationCode and "신규"(opened new station) in remarks
            if occursin(re_opened, string(row[8]))
                # parse date (year & month only)
                m = match(re_opened_date, string(row[8]))
                opened_year = parse(Int, m.captures[1])
                opened_month = parse(Int, m.captures[2])

                opened_date = ZonedDateTime(opened_year, opened_month, 1, tz"Asia/Seoul")
            end

            # parse date of station closed
            if occursin(re_closed, string(row[8]))
                m = match(re_closed_date, string(row[8]))
                closed_year = parse(Int, m.captures[1])
                closed_month = parse(Int, m.captures[2])
                closed_day = m.captures[4]
                if closed_day == nothing
                    closed_day = day(Dates.lastdayofmonth(DateTime(closed_year, closed_month)))
                else
                    closed_day = parse(Int, m.captures[4])
                end

                closed_date = ZonedDateTime(closed_year, closed_month, closed_day, tz"Asia/Seoul")
            end

            # parse date of station moved
            #=
            if hasvalue(row[3]) == true && occursin(re_moved, string(row[8]))
                if row[3] !=
                    station_code = int(row[3])
                    station_name = string(row[4])
                else
                    station_code = int(prev_row[3])
                    station_name = string(prev_row[4])
                end

                m = match(re_closed_date, string(row[8])
                moved_year = m.captures[1]
                moved_month = m.captures[2]
                opened_date = ZonedDateTime(int(moved_year), int(moved_month), 1, 0, 0, tz"Asia/Seoul")

                # compute closed date of modified one
                closed_date_prev = opened_date - Minute(1)

                # find exising one and modify closed_date
                # TODO: not sure how to do if there are multiple movement of station
                prev = @from i in stn_df begin
                    @where i.statioCode == station_code
                    @select i
                    @collect DataFrame
                end

                station_code = Msrstn_list[prev_msrstn_idx][stationCode]
                station_name = Msrstn_list[prev_msrstn_idx][stationName]

                Msrstn_list[prev_Msrstn_idx][closedDate] = closed_date_ex.strftime(strftime_fmt)
            end
            =#

            # check type
            @assert isa(station_code, Int)
            @assert isa(long_X, Real)
            @assert isa(lat_Y, Real)
            push!(df, (station_code, station_name, opened_date, closed_date, long_X, lat_Y))
        end
    end

    # check path
    mkpath(input_dir)
    filename = joinpath(input_dir, "msrstn_korea.csv")
    CSV.write(filename, df)

    df
end

function parse_aerosols(aes_dir::String, input_dir::String)
    re_aes_fn = r"([0-9]+)년\W*([0-9]+)(분기|월).csv"
    re_ext_fn = r".*.csv"
    date_str = "yyyymmddHH"

    aes_paths = joinpath(input_dir, aes_dir)
    aes_globs = rglob([re_aes_fn], aes_paths)

    # check there is a input file
    @assert isempty(aes_globs) == false

    df = DataFrame(stationCode = [], date = [],
                   SO2 = [], CO = [], O3 = [], NO2 = [],
                   PM10 = [], PM25 = [])

    p = Progress(length(collect(aes_globs)), dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)
    for _aes in aes_globs
        #=
        if match(re_ext_fn, _aes)[1] == "csv"
            df_raw = CSV.read(_aes)
        elseif match(re_ext_fn, _aes)[1] == "xlsx"
            df_raw = ExcelReaders.readxlsheet("Filename.xlsx", "Sheet1")
        else
            continue
        end
        =#
        df_raw = CSV.read(_aes, copycols=true)

        DataFrames.rename!(df_raw, [:지역 => :region, :측정소코드 => :stationCode, :측정소명 => :stationName, :측정일시 => :date,
                   :주소 => :addr])
        dropmissing!(df_raw, :date)

        # convert to ZonedDateTime to parse time correctly (i.e. 20181231T24:00 and 20190101T00:00)
        dates = [ZonedDateTime(DateTime(string(d), date_str), tz"Asia/Seoul") for d in df_raw[!, :date]]

        df_tmp = DataFrame(stationCode = df_raw[!, :stationCode],
            date = dates,
            SO2 = df_raw[!, :SO2],
            CO = df_raw[!, :CO],
            O3 = df_raw[!, :O3],
            NO2 = df_raw[!, :NO2],
            PM10 = df_raw[!, :PM10],
            PM25 = df_raw[!, :PM25])

        df = vcat(df, df_tmp)

        ProgressMeter.next!(p; showvalues = [(:filename, _aes)])
        flush(stdout); flush(stderr)
    end

    sort!(df, :date)

    df
end

function parse_weathers(wea_dir::String, input_dir::String, wea_stn_code::Integer)
    re_wea_fn = Regex("SURFACE_ASOS_$(string(wea_stn_code))_HR_([0-9]+)_([0-9]+)_([0-9]+).csv")
    date_str = "yyyy-mm-dd HH:MM"

    wea_paths = joinpath(input_dir, wea_dir)
    wea_globs = rglob([re_wea_fn], wea_paths)
    @assert isempty(wea_globs) == false

    df = DataFrame(date = [], temp = [],
        u = [], v = [], wind_spd = [], wind_dir = [], wind_sdir = [], wind_cdir = [],
        pres = [], humid = [], prep = [], snow = [])

    p = Progress(length(collect(wea_globs)), dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)
    for _wea in wea_globs

        df_raw = CSV.read(_wea, copycols=true, normalizenames=true)

        DataFrames.rename!(df_raw, [Symbol("일시") => :date,
            Symbol("기온(°C)") => :temp,
            Symbol("강수량(mm)") => :prep,
            Symbol("풍속(m/s)") => :wind_spd,
            Symbol("풍향(16방위)") => :wind_dir,
            Symbol("습도(%)") => :humid,
            Symbol("현지기압(hPa)") => :pres,
            Symbol("적설(cm)") => :snow])
        dropmissing!(df_raw, :date)

        dates = [ZonedDateTime(DateTime(string(d), date_str), tz"Asia/Seoul") for d in df_raw[!, :date]]

        # missings to zero except temp, humid, and pres (pressure)
        df_raw[!, :wind_spd] = coalesce.(df_raw[!, :wind_spd], 0.0)
        df_raw[!, :wind_dir] = coalesce.(df_raw[!, :wind_dir], 0.0)
        df_raw[!, :prep] = coalesce.(df_raw[!, :prep], 0.0)
        df_raw[!, :snow] = coalesce.(df_raw[!, :snow], 0.0)

        # http://colaweb.gmu.edu/dev/clim301/lectures/wind/wind-uv

        dict_wind_dir = Dict([
            ("0", 90),
            ("360", 90 ), ("20" , 67.5 ), ("50" , 45 ), ("70", 22.5),
            ("90" , 0  ), ("110", 337.5), ("140", 315), ("160", 292.5),
            ("180", 270), ("200", 247.5), ("230", 225), ("250", 202.5),
            ("270", 180), ("290", 157.5), ("320", 135), ("340", 112.5)])

        _wind_cdir = [Base.Math.cos(Base.Math.deg2rad(-dict_wind_dir[string(Int(w))] + 90)) for w in df_raw[!, :wind_dir]]
        _wind_sdir = [Base.Math.sin(Base.Math.deg2rad(-dict_wind_dir[string(Int(w))] + 90)) for w in df_raw[!, :wind_dir]]

        _u = [w[1] * Base.Math.cos(Base.Math.deg2rad(-dict_wind_dir[string(Int(w[2]))] + 90)) for w in zip(df_raw[!, :wind_spd], df_raw[!, :wind_dir])]
        _v = [w[1] * Base.Math.sin(Base.Math.deg2rad(-dict_wind_dir[string(Int(w[2]))] + 90)) for w in zip(df_raw[!, :wind_spd], df_raw[!, :wind_dir])]

        df_tmp = DataFrame(
            date = dates,
            temp = df_raw[!, :temp],
            u = _u,
            v = _v,
            wind_spd = df_raw[!, :wind_spd],
            wind_dir = df_raw[!, :wind_dir],
            wind_sdir = _wind_sdir,
            wind_cdir = _wind_cdir,
            pres = df_raw[!, :pres],
            humid = df_raw[!, :humid],
            prep = df_raw[!, :prep],
            snow = df_raw[!, :snow])

        df = vcat(df, df_tmp)

        ProgressMeter.next!(p; showvalues = [(:filename, _wea)])
        flush(stdout); flush(stderr)
    end

    sort!(df, :date)

    df
end

"""
    join_data(input_dir, obs_path, aes_dir, wea_dir, start_date, end_date)
join dataframe from 3 input
input_dir: input directory  (i.e. "/input")
obs_path : station xlsx path (full path)
aes_dir : subdirectory name of aerosol data (i.e. "aerosol")
wea_dir : subdirectory name of weather data (i.e. joinpath("weather", "seoul"))
"""
function join_data(input_dir::String, obs_path::String, aes_dir::String, wea_dir::String, start_date::ZonedDateTime, end_date::ZonedDateTime)
    seoul_stn_code = 108

    @info "Parsing Station dataset..."
    flush(stdout); flush(stderr)
    df_obs = parse_obsxlsx(obs_path, input_dir)
    @info "Parsing Weather dataset..."
    flush(stdout); flush(stderr)
    df_wea = parse_weathers(wea_dir, input_dir, seoul_stn_code)
    @info "Parsing Aerosol dataset..."
    flush(stdout); flush(stderr)
    df_aes = parse_aerosols(aes_dir, input_dir)
    # @info "Parsing Vertical Wind dataset..."
    # flush(stdout); flush(stderr)
    # df_vwn = parse_vwind(vwind_dir, input_dir

    # filter in date range
    df_obs = df_obs[(df_obs.oDate .<= start_date) .& (df_obs.cDate .>= end_date), :]
    df_aes = df_aes[(df_aes.date .>= start_date) .& (df_aes.date .<= end_date), :]
    df_wea = df_wea[(df_wea.date .>= start_date) .& (df_wea.date .<= end_date), :]
    # df_vwn = df_vwn[(df_wea.date .>= start_date) .& (df_wea.date .<= end_date), :]

    df1 = join(df_obs, df_aes, on = :stationCode)
    df2 = join(df1, df_wea, on = :date, kind = :inner)
    # df3 = join(df2, df_vwn, on = :date, kind = :inner)

    # sort by stationCode and date
    sort!(df2, (:stationCode, :date))

    # scaling except :PM10, :PM25
    features = [:SO2, :CO, :O3, :NO2, :PM10, :PM25,
        :temp, :u, :v, :wind_sdir, :wind_cdir, :wind_spd, :wind_dir, :pres, :humid, :prep, :snow]

    df = df2[!, [:stationCode, :date, :lat, :lon,
                   :SO2, :CO, :O3, :NO2, :PM10, :PM25,
                   :temp, :pres,
                   :u, :v, :wind_spd, :wind_dir, :wind_sdir, :wind_cdir,
                   :humid, :prep, :snow]]

    # check path
    mkpath(input_dir)
    filename = joinpath(input_dir, "input_2021.csv")
    CSV.write(filename, df)

    df
end


text:
mutable struct MinimaxModel <: AbstractModel
    settings::Settings
    constraints::Vector{Constraint}
    obj::MinimaxObjective
    output::Output
    MinimaxModel(settings, constraints, obj, output) =
        new(settings, constraints, obj, output)
    MinimaxModel() = new(Settings(), Constraint[], MinimaxObjective(), Output())
end


text:
@testset "Renewable data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    warn_message = "The data doesn't include devices of type RenewableDispatch, consider changing the device models"
    model = DeviceModel(PSY.RenewableDispatch, PSI.RenewableFullDispatch)
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5)
    @test_logs (:warn, warn_message) construct_device!(op_model, :Renewable, model)
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys14)
    @test_logs (:warn, warn_message) construct_device!(op_model, :Renewable, model)
end

@testset "Renewable DCPLossLess FullDispatch" begin
    model = DeviceModel(PSY.RenewableDispatch, PSI.RenewableFullDispatch)

    #5 Bus testing case
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_re)
    construct_device!(op_model, :Renewable, model)
    moi_tests(op_model, false, 72, 0, 72, 0, 0, false)

    psi_checkobjfun_test(op_model, GAEVF)

    # Using Parameters Testing
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_re; parameters = true)
    construct_device!(op_model, :Renewable, model)
    moi_tests(op_model, true, 72, 0, 72, 0, 0, false)
    psi_checkobjfun_test(op_model, GAEVF)

    # No Forecast - No Parameters Testing
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_re; forecast = false)
    construct_device!(op_model, :Renewable, model)
    moi_tests(op_model, false, 3, 0, 3, 0, 0, false)
    psi_checkobjfun_test(op_model, GAEVF)

    # No Forecast Testing
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_re; parameters = true, forecast = false)
    construct_device!(op_model, :Renewable, model)
    moi_tests(op_model, true, 3, 0, 3, 0, 0, false)
    psi_checkobjfun_test(op_model, GAEVF)
end

@testset "Renewable ACPPower Full Dispatch" begin
    model = DeviceModel(PSY.RenewableDispatch, PSI.RenewableFullDispatch)
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5_re; parameters = p)
        construct_device!(op_model, :Renewable, model)
        if p
            moi_tests(op_model, p, 144, 24, 72, 0, 48, false)
            psi_checkobjfun_test(op_model, GAEVF)
        else
            moi_tests(op_model, p, 144, 24, 72, 0, 48, false)

            psi_checkobjfun_test(op_model, GAEVF)
        end
    end
    # No Forecast Test
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5_re; forecast = false, parameters = p)
        construct_device!(op_model, :Renewable, model)
        if p
            moi_tests(op_model, p, 6, 1, 3, 0, 2, false)
            psi_checkobjfun_test(op_model, GAEVF)
        else
            moi_tests(op_model, p, 6, 1, 3, 0, 2, false)

            psi_checkobjfun_test(op_model, GAEVF)
        end
    end
end

@testset "Renewable DCPLossLess ConstantPowerFactor" begin
    model = DeviceModel(PSY.RenewableDispatch, PSI.RenewableConstantPowerFactor)

    #5 Bus testing case
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_re)
    construct_device!(op_model, :Renewable, model)
    moi_tests(op_model, false, 72, 0, 72, 0, 0, false)

    psi_checkobjfun_test(op_model, GAEVF)

    # Using Parameters Testing
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_re; parameters = true)
    construct_device!(op_model, :Renewable, model)
    moi_tests(op_model, true, 72, 0, 72, 0, 0, false)
    psi_checkobjfun_test(op_model, GAEVF)

    # No Forecast - No Parameters Testing
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_re; forecast = false)
    construct_device!(op_model, :Renewable, model)
    moi_tests(op_model, false, 3, 0, 3, 0, 0, false)
    psi_checkobjfun_test(op_model, GAEVF)

    # No Forecast Testing
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_re; parameters = true, forecast = false)
    construct_device!(op_model, :Renewable, model)
    moi_tests(op_model, true, 3, 0, 3, 0, 0, false)
    psi_checkobjfun_test(op_model, GAEVF)
end

@testset "Renewable ACPPower ConstantPowerFactor" begin
    model = DeviceModel(PSY.RenewableDispatch, PSI.RenewableConstantPowerFactor)
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5_re; parameters = p)
        construct_device!(op_model, :Renewable, model)
        if p
            moi_tests(op_model, p, 144, 0, 72, 0, 72, false)
            psi_checkobjfun_test(op_model, GAEVF)
        else
            moi_tests(op_model, p, 144, 0, 72, 0, 72, false)

            psi_checkobjfun_test(op_model, GAEVF)
        end
    end
    # No Forecast Test
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5_re; forecast = false, parameters = p)
        construct_device!(op_model, :Renewable, model)
        if p
            moi_tests(op_model, p, 6, 0, 3, 0, 3, false)
            psi_checkobjfun_test(op_model, GAEVF)
        else
            moi_tests(op_model, p, 6, 0, 3, 0, 3, false)

            psi_checkobjfun_test(op_model, GAEVF)
        end
    end
end

@testset "Renewable DCPLossLess FixedOutput" begin
    model = DeviceModel(PSY.RenewableDispatch, PSI.RenewableFixed)
    for p in [true, false]
        op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_re; parameters = p)
        construct_device!(op_model, :Renewable, model)
        if p
            moi_tests(op_model, p, 0, 0, 0, 0, 0, false)
            psi_checkobjfun_test(op_model, GAEVF)
        else
            moi_tests(op_model, p, 0, 0, 0, 0, 0, false)
            psi_checkobjfun_test(op_model, GAEVF)
        end
    end
end

@testset "Renewable ACPPowerModel FixedOutput" begin
    model = DeviceModel(PSY.RenewableDispatch, PSI.RenewableFixed)
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5_re; parameters = p)
        construct_device!(op_model, :Renewable, model)
        if p
            moi_tests(op_model, p, 0, 0, 0, 0, 0, false)
            psi_checkobjfun_test(op_model, GAEVF)
        else
            moi_tests(op_model, p, 0, 0, 0, 0, 0, false)
            psi_checkobjfun_test(op_model, GAEVF)
        end
    end
end