text:
# developed with Julia 1.1.1
#
# struct for EMS simulation

abstract type AbstractController end

mutable struct DummyController <: AbstractController
end

mutable struct Result
    cost::Array{Float64,1}
    soc::Array{Float64,1}
end

Result(h::Int64) = Result(zeros(h), zeros(h))

struct Id
    site_id::String
    period_id::String
    price_id::String
    model_type::String
end

struct Simulation
    result::Result
    timer::Array{Float64,1}
    id::Id
end

struct Battery
    capacity::Float64
    power::Float64
    charge_efficiency::Float64
    discharge_efficiency::Float64
end

struct Site
    id::String
    battery::Battery
    path_to_test_data_csv::Union{String, Nothing}
    path_to_train_data_csv::Union{String, Nothing}
    path_to_save_folder::String
end

function Site(data::DataFrame, row::Int64, path_to_test_data_folder::Union{String, Nothing},
    path_to_train_data_folder::Union{String, Nothing}, path_to_save_jld_file::String)

    id = string(data[row, :site_id])
    battery = Battery([float(x) for x in data[row, 3:end]]...)
    if path_to_test_data_folder !== nothing
        path_to_test_data_csv = joinpath(path_to_test_data_folder, id*".csv.gz")
    else
        path_to_test_data_csv = nothing
    end
    if path_to_train_data_folder !== nothing
        path_to_train_data_csv = joinpath(path_to_train_data_folder, id*".csv.gz")
    else
        path_to_train_data_csv = nothing
    end

    return Site(id, battery, path_to_test_data_csv, path_to_train_data_csv,
        path_to_save_jld_file)

end

mutable struct Period
    id::String
    data::DataFrame
    site::Site
end

struct Prices
    name::String
    buy::Array{Float64,1}
    sell::Array{Float64,1}
end

struct Information
    t::Int64
    soc::Float64
    pv::Array{Float64,1}
    forecast_pv::Array{Float64,1}
    load::Array{Float64,1}
    forecast_load::Array{Float64,1}
    prices::Prices
    battery::Battery
    site_id::String
end

function Information(t::Int64, prices::Prices, period::Period, soc::Float64)

    data = sort!(period.data[t+1:t+96, :], rev=true)
    pv = data[!, :actual_pv]
    forecast_pv = Vector{Float64}(data[end, 102:197])
    load = data[!, :actual_consumption]
    forecast_load = Vector{Float64}(data[end, 6:101])

    return Information(t, soc, pv, forecast_pv, load, forecast_load, prices, period.site.battery,
        period.site.id)
end


text:
export MtlResource, storage, options, cpuCacheMode, hazardTrackingMode

const MTLResource = Ptr{MtResource}

abstract type MtlResource end

Base.unsafe_convert(::Type{MTLResource}, res::MtlResource) = Base.bitcast(MTLResource, res.handle)

device(res::MtlResource) = convert(MtlDevice, mtResourceDevice(res))

function label(l::MtlResource)
    ptr = mtResourceLabel(l)
    return ptr == C_NULL ? "" : unsafe_string(ptr)
end

cpuCacheMode(res::MtlResource) = mtResourceCPUCacheMode(res)
storage(res::MtlResource) = mtResourceStorageMode(res)
hazardTrackingMode(res::MtlResource) = mtResourceHazardTrackingMode(res)
options(res::MtlResource) = mtResourceOptions(res)

ext: jl


text:
def flatten(arr) {
	result = []
	for i=0, i < |arr|, i=i+1
        if arr[i] is Array
            for j=0, j < |arr[i]|, j=j+1
                result :<= arr[i][j]
        else result :<= arr[i]
	result
}

def join(arr, seq) {
	result = []
    if (0 < |arr|) {
        if arr[0] is Array
            for a : arr[0]
                result :<= a
        else result :<= arr[0]

        for i=1, i < |arr|, i=i+1 {
            for c : seq
                result :<= c
            if arr[i] is Array
                for a : arr[i]
                    result :<= a
            else result :<= arr[i]
        }
    }
    result
}

def map(arr, f) {
    res = []
    for a : arr
        res :<= f(a)
    res
}

def filter(arr, f) {
    res = []
    for a : arr
        if f(a)
            res :<= a
    res
}

; Future idea: Allow type checking in function declaration
; Capitalized are classes, lowercase are consistency checking.
; Allow _ or ?, decide on one for types we don't care about.
; def collect(Array arr, @(a,a)->a f) -> a
def collect(arr, f) {
    res = arr[0]
    for i=1, i < |arr|, i=i+1
        res = f(res, arr[i])
    res
}

def collecti(arr, f, i) {
    res = i
    for a : arr
        res = f(res, a)

    res
}

def foreach(arr, f)
    for a : arr
        f(a)

def max(array, cmp) {
    max_elt = array[0]
    for i=1, i<|array|, i=i+1
        if cmp(array[i],max_elt) > 0
            max_elt = array[i]
    max_elt
}

def min(array, cmp) {
    max_elt = array[0]
    for i=1, i<|array|, i=i+1
        if cmp(array[i],max_elt) < 0
            max_elt = array[i]
    max_elt
}

def contains(array, elt) {
    for a : array
        if a == elt
            return True
    False
}

; A bad impl of unique
def unique(array) {
    res = []
    for i=0, i < |array|, i=i+1
        if ~contains(res, array[i])
            res :<= array[i]
    res
}

max_stars_repo_name: jeffreymanzione/jeff-vm-lang

max_issues_repo_name: jeffreymanzione/jeff-vm-lang

max_forks_repo_name: jeffreymanzione/jeff-vm-lang


text:
## PART 4

println("\nPart 4");
include("../q2/jacknife.jl");

# Set bin size b to the largest autocorrelation time for each temp
numBodies = 1000;
println("\n\tFinding number of steps between independent data recordings\n")
@show bp_1069 = 2maximum([τtemp_1069, τpe_1069, τvir_1069])
@show bp_1304 = 2maximum([τtemp_1304, τpe_1304, τvir_1304])

b=10
println("\n\tPick a bin size of $b\n")

# Function for calculating pressure at various temperatures
function pressure(numBodies, temperatures, virials)
    # Take the time average of the virial
    virialavg = mean(virials);

    # Set Boltzmann constant
    kB = 1.3806e-23;

    kBTs = temperatures * kB;

    return 1 - virialavg * ( 1 ./ (6 * numBodies * kBTs ) );
end

println("\n\tEstimating σpressure:");

mdjack_1069 = jacknife(md_1069, lengthmd, b);
mdjack_1304 = jacknife(md_1304, lengthmd, b);

# Define a function that we can pass to jacknifeσ
pressureN(md) = pressure(numBodies, md[:,1], md[:,3]);

# Estimate σpressure
σpressure_1069 = jacknifeσ(pressureN, mdjack_1069);
σpressure_1304 = jacknifeσ(pressureN, mdjack_1304);

println("\n\tAverage pressure for T = 1.069 is ", mean(pressureN(md_1069)))
println("\n\tAverage pressure for T = 1.304 is ", mean(pressureN(md_1304)))

println("\n\tStandard deviation of pressure for T = 1.069: $σpressure_1069")
println("\n\tStandard deviation of pressure for T = 1.304: $σpressure_1304")

# Should also plot pressure


text:
module Radixsort

include("sort.jl")
#export BInt, BUInt, bitsizeof,  AsciiChar, Latin1Char

end


text:
# Autogenerated wrapper script for P4est_MPICH_jll for x86_64-linux-gnu-libgfortran4
export libp4est, libsc

using Zlib_jll
using MPICH_jll
using CompilerSupportLibraries_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `libp4est`
const libp4est_splitpath = ["lib", "libp4est.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libp4est_path = ""

# libp4est-specific global declaration
# This will be filled out by __init__()
libp4est_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libp4est = "libp4est-2.2.so"

# Relative path to `libsc`
const libsc_splitpath = ["lib", "libsc.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libsc_path = ""

# libsc-specific global declaration
# This will be filled out by __init__()
libsc_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libsc = "libsc-2.2.so"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"P4est_MPICH")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Zlib_jll.PATH_list, MPICH_jll.PATH_list, CompilerSupportLibraries_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Zlib_jll.LIBPATH_list, MPICH_jll.LIBPATH_list, CompilerSupportLibraries_jll.LIBPATH_list,))

    global libp4est_path = normpath(joinpath(artifact_dir, libp4est_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libp4est_handle = dlopen(libp4est_path)
    push!(LIBPATH_list, dirname(libp4est_path))

    global libsc_path = normpath(joinpath(artifact_dir, libsc_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libsc_handle = dlopen(libsc_path)
    push!(LIBPATH_list, dirname(libsc_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
# `info` returns a dictionary of model trait values suitable, after
# encoding, to serializing to TOML file. Not intended to be exposed to
# user. The `traits` function, defined in MLJ, returns the trait
# values as a named-tuple, more friendly for user-interaction. One can
# similarly call `traits` on performance measures.

info_dict(M::Type{<:Supervised}) = info_dict(M, SUPERVISED_TRAITS)
info_dict(M::Type{<:Unsupervised}) = info_dict(M, UNSUPERVISED_TRAITS)
info_dict(model::Model) = info_dict(typeof(model))

function info_dict(M::Type{<:Model}, traits)

    message = "$M has a bad trait declaration.\n"

    ismissing(is_pure_julia(M)) || is_pure_julia(M) isa Bool ||
        error(message*"is_pure_julia($M) must return true, false or missing. ")

    :supports_weights in traits &&
        !(supports_weights(M) isa Bool) &&
        error(message*"supports_weights($M) must return true, false or missing. ")

    :is_wrapper in traits &&
        !(is_wrapper(M) isa Bool) &&
        error(message*"is_wrapper($M) must return true, false or missing. ")

    return LittleDict{Symbol,Any}(trait => eval(:($trait))(M)
                                  for trait in traits)
end


text:
abstract abstract_linear_system_solver;

include("julia.jl")
#include("matlab.jl")
include("mumps.jl")

function initialize!(solver::abstract_linear_system_solver, SparseMatrix::SparseMatrixCSC{Float64,Int64}, sym::Symbol)
		solver._SparseMatrix = SparseMatrix;
		solver.sym = sym
end

function finalize!(solver::abstract_linear_system_solver)

end

function inertia_status(pos_eigs::Int64, neg_eigs::Int64, zero_eigs::Int64, num_vars::Int64, num_constraints::Int64)
	try
		@assert(pos_eigs + neg_eigs + zero_eigs == num_vars + num_constraints)

		#println(pos_eigs, " ", neg_eigs);

		if pos_eigs == num_vars && neg_eigs == num_constraints
			return true
		elseif pos_eigs > num_vars || neg_eigs + pos_eigs != num_vars + num_constraints
			println("Warning: numerical instability in LDL factorization")
			if pos_eigs > num_vars
				println("more positive eigenvalues than variables")
				return false
			else
				println("zero eigenvalues")
				return false
			end
		else
			return false
		end
	catch e
		println("ERROR in compute_inertia")
		throw(e)
	end
end

max_forks_repo_name: ohinder/Everest


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CSV
using DataFrames
using Dates
using Plots
using BSON
using StatsBase
using Underscores
using ColorSchemes
using PrettyTables
colors = ColorSchemes.tableau_10.colors

## ===========================================================================

function consistent!(raw; nonnegative = false, positive = false)
    df = @_ raw |>
        disallowmissing(__) |>
        sort(__, :date) |>
        continuous_time!(__)

    if nonnegative || positive
        dat = @_ df |>
            select(__, Not(:date)) |>
            Array(__)
        nonnegative && @assert all( dat .>= 0 )
        positive && @assert all( dat .> 0 )
    end
    df
end

function Base.filter!(conditions::Vector{Pair{Symbol, Function}}, df::AbstractDataFrame)
    for condition in conditions
        filter!(condition, df)
    end
    df
end
# Base.filter!(condition::Pair{Symbol, Function}, df::AbstractDataFrame) =
#     filter!(condition, df)
Base.filter!(nothing, df::AbstractDataFrame) = df

function read_survey(region; fnames_survey, shift = 0, kwargs...)
    @_ fnames_survey[region] |>
        CSV.read(__, DataFrame) |>
        consistent!(__, positive = true) |>
        transform(__, :date =>( x -> x + Day(shift) )=> :date)
end

function read_mobil(region; fnames_mobil, kwargs...)
    @_ fnames_mobil[region] |>
        CSV.read(__, DataFrame) |>
        consistent!(__)
end

# read_predictors(region; predictors = All(), predictors_kwargs...)
# read_mobil(region; predictors_kwargs...)

function read_predictors(region; predictors = All(), kwargs...)
    @_ outerjoin(
            read_survey(region; kwargs...),
            read_mobil(region; kwargs...);
            on=:date) |>
        select(__, :date, predictors) |>
        sort(__, :date)
end

function read_epidata(region; conditions = nothing, kwargs...)
    @_ projectdir("data/epidata/epidata_region=$(region).csv") |>
        CSV.read(__, DataFrame) |>
        consistent!(__, nonnegative = true) |>
        filter!(conditions, __)
end

function relativechange!(xs::AbstractVector; standardize=false, kwargs...)
    xs[1] == 0. && (@warn "predictor data stream starts with 0, so we assume it is already normalized."; return nothing)
    xs ./= abs(xs[1])
    xs .-= xs[1]
    standardize && (xs ./= std(xs))
    return nothing
end

function relativechange!(df::AbstractDataFrame; kwargs...)
    for n in names(df, Not(:date))
        relativechange!(df[!,n], kwargs...)
    end
    df
end

function continuous_time!(df)
    @assert all( df.date[1]:Day(1):df.date[end] .== df.date )
    df
end

function clean_predictors(pred, epid; conditions, kwargs...)
    df = @_ pred |>
        leftjoin(epid, __; on = :date) |>
        sort(__, :date) |>
        filter!(conditions, __) |>
        select(__, names(pred)) |>
        disallowmissing(__) |>
        relativechange!(__) |>
        continuous_time!(__)
end

function join_predictors!(epid, pred)
    @_ epid |>
        continuous_time!(__) |>
        disallowmissing!(__) |>
        leftjoin(__,
            continuous_time!(pred);
            on=:date) |>
        sort(__, :date) |>
        coalesce.(__, 0.) |>
        disallowmissing!(__)
end

function stepindex(n, j0 = 1; step = 1, kwargs...)
    index = Vector{Int64}(undef, n)
    j = j0
    for i in 1:n
        index[i] = j
        i % step == 0 && ( j += 1 )
    end
    index
end

function earlyphase!(df, startdate)
    j = findlast( <(startdate), df.date)
    !isnothing(j) && (df[1:j,:rw] .= 1)
end

function latephase!(df)
    i = findlast(x->!ismissing(x), df[!, :rw])
    df[i+1:end, :rw] .= df[i, :rw]
end

function last_step(df)
    i = findlast(x->!ismissing(x), df[!, :rw])
    isnothing(i) && return 1
    return df[i, :rw]
end

function step!(df, sdate, edate; kwargs...)
    i0 = last_step(df) + 1
    is = findfirst(==(sdate), df.date)
    ie = findfirst(==(edate), df.date)
    n  = ie - is + 1
    df[is:ie, :rw] = stepindex(n, i0; kwargs...)
end

function wait!(df, sdate, edate)
    i0 = last_step(df)
    is = findfirst(==(sdate), df.date)
    ie = findfirst(==(edate), df.date)
    df[is:ie, :rw] .= i0
end

function randomwalk!(data, pred; semiparam = true, startdate = nothing, enddate = nothing, kwargs...)
    rw = DataFrame(:date => data.date, :rw => Vector{Union{Missing, Int64}}(missing, size(data,1)))
    startdate = isnothing(startdate) ? rw.date[1] : startdate
    enddate   = isnothing(enddate) ? data.date[end] : enddate

    earlyphase!(rw, startdate)
    if semiparam
        s = pred.date[1]
        e = pred.date[end]
        step!(rw, startdate, s; kwargs...)
        wait!(rw, s + Day(1), e)
        step!(rw, e + Day(1), enddate; kwargs...)
    else
        step!(rw, startdate, enddate; kwargs...)
    end
    latephase!(rw)

    @_ leftjoin(data, rw; on=:date) |>
        sort(__, :date)# |>
        # disallowmissing!(__)
end

function load_data(region, epidata_kwargs, predictor_kwargs, rw_kwargs)
    pred = read_predictors(region ; predictor_kwargs... )
    epid = read_epidata(region; epidata_kwargs...)
    pred = clean_predictors(pred, epid; predictor_kwargs... )
    data = @_ join_predictors!(epid, pred) |>
                randomwalk!(__, pred; rw_kwargs... ) |>
                continuous_time!(__) |>
                disallowmissing!(__)
    data[!,:region] .= region
    data
end
## ===========================================================================
regions = ["France", "Hungary", "UK", "Sweden", "Denmark", "Germany", "Italy", "USA"]
params = (
      epidata_kwargs = Dict(
        :conditions => [:date => x -> x >= Date("2020-08-01")]
        )

    , predictors_kwargs = Dict(
        :fnames_survey => Dict(region => projectdir("data/contacts/int/smoothed_contact_rates_region=$(region).csv") for region in regions),
        :fnames_mobil  => Dict(region => projectdir("data/mobility/mobility_country=$(region).csv") for region in regions),
        :shift         => -1,
        :predictors    => [:strangers],#[:strangers, :google, :apple],
        :conditions    => [
            :cases_per_100000    => x -> x >
            20,
            :fraction_vaccinated => x -> x < 5/100
            ]
        )

    , rw_kwargs = Dict(
        :step      => 7,
        :startdate => nothing,
        :enddate   => nothing,
        :semiparam => true
        )
    )

function load_data(regions::AbstractVector, args...)
    data = load_data(regions[1], params...)
    for region in regions[2:end]
        @show region
        data = vcat(data,
            load_data(region, args...)
            )
    end
    data
end

load_data(regions, params...)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
include("common.jl")
jllversion=v"1.7.0"
build_julia(ARGS, v"1.6.3"; jllversion)
build_julia(ARGS, v"1.7.0-rc1"; jllversion)
build_julia(ARGS, v"1.8.0-DEV"; jllversion)

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
module CollegeStratPlots

using DocStringExtensions, StatsBase
using Plots, Plots.PlotMeasures, StatsPlots
using FilesLH
using CollegeStratBase

export plot_settings, plot_setting, plot_defaults, subplot_layout
export blank_plot, sub_plots, histogram_plot, scatter_plot, contour_plot
export bar_graph, grouped_bar_graph
export line_plot, add_line!, add_title!, plot_text
export figsave

include("plot_settings.jl");
include("saving.jl");
include("bar_graphs.jl");

"""
	$(SIGNATURES)

Make a legend (which is a Matrix) from a vector of strings.
"""
make_legend(v) = reshape(v, 1, length(v));

"""
	$(SIGNATURES)

Make layout for subplots of `nPlots` equally sized plots.
"""
function subplot_layout(nPlots)
    if iseven(nPlots)
        nRows = round(Int, nPlots / 2);
        nCols = round(Int, 2);
    elseif nPlots == 3
        nRows = 1;
        nCols = nPlots;
    else
        error("Not implemented for $nPlots");
    end
    @assert (nRows * nCols) == nPlots  "Adding up fails: $nRows, $nCols"
    return nRows, nCols
end

"""
	$(SIGNATURES)

Render a set of plots as subplots.
"""
function sub_plots(pV; kwargs...)
    p = plot(pV...; kwargs...);
    return p
end

add_title!(p, titleStr; kwargs...) = title!(p, titleStr; kwargs...);

plot_text(lbl, sz, pos) = Plots.text(lbl, sz, pos);

# """
# 	$(SIGNATURES)

# Write table by [gpa, parental].
# """
# function table_by_gpa_parental(io,  data_gpM :: Matrix;  digits = 3)
#     ng, np = size(data_gpM);
#     if eltype(data_gpM) <: Number
#         data_gpM = string.(round.(data_gpM, digits = digits));
#     end
#     pretty_table(io,  hcat(gpa_labels(ng), data_gpM),  vcat(" ", parental_labels(np)));
#     return nothing;
# end

# # Grouped bar graph. Groups are dropouts and graduates.
# function plot_by_qual_grad(x_qgM, yStr :: AbstractString, fPath :: AbstractString)
#     # nq = size(x_qgM, 1);
#     p = grouped_bar_graph(["Dropouts", "Graduates"], x_qgM';
#         leg = :none, ylabel = yStr);

#     !isempty(fPath)  &&  figsave(p, fPath; dataM = x_qgM);
#     return p
# end

"""
	$(SIGNATURES)

Line graph. Simple wrapper around `plot`.
"""
function line_plot(xV, yM; kwargs ...)
    p = plot(xV, yM; kwargs...);
    return p
end

function line_plot(yM; kwargs ...)
    p = plot(yM; kwargs...);
    return p
end

"""
	$(SIGNATURES)

Add line to a plot.
"""
function add_line!(p, x, y; kwargs...)
    plot!(p, x, y; kwargs...);
end

"""
	$(SIGNATURES)

Scatter plot.
"""
function scatter_plot(x, y; kwargs...)
    p = scatter(x, y; kwargs...);
    return p
end

function histogram_plot(dataV; kwargs...)
    h = fit(Histogram, dataV);  # ++++++ where?
    p = plot(h; kwargs...)
    return p
end

function contour_plot(x, y, z; kwargs...)
    p = contour(x, y, z; kwargs...);
    return p
end

"""
	$(SIGNATURES)

Blank plot with options.
"""
blank_plot(; kwargs...) = plot(; kwargs...);

end


text:
using JMlem
using Test

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

max_stars_repo_name: jjgomezcadenas/JMlem

max_issues_repo_name: jjgomezcadenas/JMlem

max_forks_repo_name: jjgomezcadenas/JMlem


text:
# This file was generated, do not modify it. # hide
W = X |> Standardizer()
z = y |> UnivariateBoxCoxTransformer()

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
using sem, Feather, ModelingToolkit, Statistics, LinearAlgebra,
    Optim, SparseArrays, Test, Zygote, LineSearches, ForwardDiff

using Feather, Statistics

abstract type SemObs end

abstract type LossFunction end

miss20_dat = Feather.read("test/comparisons/dat_miss20_dat.feather")
miss30_dat = Feather.read("test/comparisons/dat_miss30_dat.feather")
miss50_dat = Feather.read("test/comparisons/dat_miss50_dat.feather")

miss20_par = Feather.read("test/comparisons/dat_miss20_par.feather")
miss30_par = Feather.read("test/comparisons/dat_miss30_par.feather")
miss50_par = Feather.read("test/comparisons/dat_miss50_par.feather")

miss20_mat = Matrix(miss20_dat)

struct SemObsMissing{
        A <: AbstractArray,
        C <: Union{AbstractArray, Nothing},
        D <: AbstractFloat,
        O <: AbstractFloat,
        P <: Vector,
        R <: Vector,
        PD <: AbstractArray,
        PS <: AbstractArray,
        PO <: AbstractArray} <: SemObs
    data::A
    obs_mean::C
    n_man::D
    n_obs::O
    patterns::P # missing patterns
    rows::R # coresponding rows in the data or matrices
    pattern_data::PD # list of data per missing pattern
    pattern_S::PS
    pattern_n_obs::PO #
end

function SemObsMissing(data; meanstructure = false)

    n_obs = Float64(size(data, 1))
    n_man = Float64(size(data, 2))

    # compute and store the different missing patterns with their rowindices
    missings = ismissing.(data)
    patterns = [missings[i, :] for i = 1:size(missings, 1)]
    remember = Vector{BitArray{1}}()
    rows = [Vector{Int64}(undef, 0) for i = 1:size(patterns, 1)]
    for i = 1:size(patterns, 1)
        unknown = true
        for j = 1:size(remember, 1)
            if patterns[i] == remember[j]
                push!(rows[j], i)
                unknown = false
            end
        end
        if unknown
            push!(remember, patterns[i])
            push!(rows[size(remember, 1)], i)
        end
    end
    rows = rows[1:length(remember)]
    n_patterns = size(rows, 1)

    # sort by number of missings
    sort_n_miss = sortperm(sum.(remember))
    remember = remember[sort_n_miss]
    remember_cart = findall.(!, remember)
    rows = rows[sort_n_miss]

    # store the data belonging to the missing patterns
    pattern_data = Vector{Array{Float64}}(undef, n_patterns)
    for i = 1:n_patterns
        pattern_data[i] = miss20_mat[rows[i], remember_cart[i]]
    end
    pattern_data = convert.(Array{Float64}, pattern_data)

    pattern_n_obs = length.(remember_cart)

    # if a meanstructure is needed, don't compute observed means
    if meanstructure
        obs_mean = nothing
        pattern_S = nothing
    else
        pattern_S = Array{Array{Float64, 2}}(undef, length(pattern_data))
        obs_mean = skipmissing_mean(data)

        for i in 1:length(pattern_data)
            S = zeros(pattern_n_obs[i], pattern_n_obs[i])
            for j in 1:size(pattern_data[i], 1)
                diff = pattern_data[i][j, :] - obs_mean[remember_cart[i]]
                S += diff*diff'
            end
            pattern_S[i] = S
        end
    end

    return SemObsMissing(data, obs_mean, n_man, n_obs, remember_cart,
    rows, pattern_data, pattern_S, Float64.(pattern_n_obs))
end

function skipmissing_mean(mat)
    means = Vector{Float64}(undef, size(mat, 2))
    for i = 1:size(mat, 2)
        @views means[i] = mean(skipmissing(mat[:,i]))
    end
    return means
end

@benchmark myobs = SemObsMissing(miss20_mat)

struct SemFIML_Loss{
        INV <: AbstractArray,
        C <: AbstractArray,
        L <: AbstractArray,
        M <: AbstractArray,
        I <: AbstractArray,
        T <: AbstractArray,
        U,
        V} <: LossFunction
    inverses::INV #preallocated inverses of imp_cov
    choleskys::C #preallocated choleskys
    logdets::L #logdets of implied covmats
    meandiff::M
    imp_inv::I
    mult::T # is this type known?
    objective::U
    grad::V
end

function SemFIML_Loss(observed::O where {O <: SemObs}, objective, grad)

    inverses = broadcast(x -> zeros(x, x), Int64.(observed.pattern_n_obs))
    choleskys = copy(inverses)

    n_patterns = size(observed.rows, 1)
    logdets = zeros(n_patterns)

    meandiff = zeros.(Int64.(observed.pattern_n_obs))

    imp_inv = zeros(Int64(observed.n_obs), Int64(observed.n_obs))
    mult = copy(inverses)

    return SemFIML_Loss(
    inverses,
    choleskys,
    logdets,
    meandiff,
    imp_inv,
    mult,
    copy(objective),
    copy(grad)
    )
end

SemFIML_Loss(myobs, 5.0, 5.0)

## Matrix inversion

LinearAlgebra.lowrankupdate()

## lavaan style matrix updates
S = rand(20,20)

S = S*transpose(S)

save = copy(S)

S_inv = inv(S)

rm_idx = [3,4]
not_idx = vcat(1,2,collect(5:20))

function lv_inv(S_inv, rm_idx, not_idx)
    A = S_inv[rm_idx, not_idx]
    H = S_inv[rm_idx,  rm_idx]
    out = S_inv[not_idx, not_idx] - transpose(A)*(H\A)
    return out
end

## structures

# struct SemFIML{O <: SemObs, I <: Imply, L <: Loss, D <: SemDiff} <: AbstractSem
#     observed::O
#     imply::I # former ram
#     loss::L # list of loss functions
#     diff::D
# end
#
# function (model::SemFIML)(par)
#     model.imply(par)
#     F = model.loss(par, model)
#     return F
# end

function (semfiml::SemFIML_Loss)(par, model::Sem{O, I, L, D}) where
            {O <: SemObs, L <: Loss, I <: Imply, D <: SemFiniteDiff}

    copyto!(semfiml.inverses[1], model.imply.imp_cov)
    semfiml.choleskys[1] = cholesky!(Hermitian(semfiml.inverses[1]); check = false)

    if !isposdef(semfiml.choleskys[1])
        F = Inf
    else
        @views for i = 2:size(inverses, 1)
            semfiml.inverses[i] .= model.imply.imp_cov[patterns[i], patterns[i]]
        end

        for i = 2:size(inverses, 1)
            semfilm.choleskys[i] = cholesky!(semfiml.inverses[i]; check = false)
        end

        #ld = logdet(a)
        logdets .= logdet.(semfiml.choleskys)

        #semml.imp_inv .= LinearAlgebra.inv!(a)
        for i = 1:size(inverses, 1)
            semfiml.inverses[i] .= LinearAlgebra.inv!(semfiml.choleskys[i])
        end

        F = zero(eltype(par))
        for i = 1:size(semfiml.rows, 1)

            F_missingpattern(
                model.imply.imp_mean,
                model.observed.obs_mean,
                semfiml.meandiff,
                model.observed.patterns,
                semfiml.inverses,
                model.observed.pattern_S,
                model.observed.pattern_data,
                semfiml.logdets,
                semfiml.mult
                model.observed.pattern_n_obs,
                F,
                i
                )
        end
    end
    return F
end

function myf(a)
    x = size(a, 1)
    x = x^2
end

using BenchmarkTools

@benchmark myf([5.0, 6.0])

# function F_missingpattern(rows, pattern, inverse, S, logdet, n_obs)
#     if size(rows, 1) == n_obs # no missings
#         F =
#     else if size(rows, 1) > 5 # find a good treshold
#
#     else if size(rows, 1) == 1 #no one else has the specific missing pattern
#
#     else #few persons, so better compute individually
#
#     end
#     return F
# end

function F_missingpattern(imp_mean::Nothing, obs_mean, meandiff,
    pattern, inverse, S, data, ld, mult, n_obs, F, i)

    # if n_obs == 1
    #     mult = inverse*S
    # else
    #     mul!(mult, inverse, S)
    # end

    mul!(mult, inverse, S)
    F += n_obs*(ld + tr(mult))

end

function F_missingpattern(imp_mean, obs_mean, meandiff,
    pattern, inverse, S, data, ld, mult, n_obs, F, i)

    F += n_obs*logdet

    @views for i = 1:n_obs
        @. meandiff = data[i, :] - imp_mean[pattern]
        F += meandiff'*inverse*meandiff
    end

end

v1 = rand(10)

v2 = rand(10)

ind = [1, 3, 5, 7, 9]

@benchmark @. getindex([v1], ind)

function f1(v1, v2, ind)
    res = v1[ind] - v2[ind]
    return res
end

function f2(v1, v2, ind)
    @views res = v1[ind] - v2[ind]
    return res
end

function f3(v1, v2, ind, res)
    res .= v1[ind] - v2[ind]
    return res
end

function f3(v1, v2, ind, res)
    @inbounds @views @. res = v1[ind] - v2[ind]
    return res
end

function f4(v1, v2, ind, res)
    for i = 1:length(ind)
        res[i] = v1[ind[i]] - v2[ind[i]]
    end
    return res
end

v1_arr = [v1]
v2_arr = [v2]

function f5(v1, v2, ind, res)
    @inbounds res .= view(v1, ind) .- view(v2, ind)
    return res
end

using BenchmarkTools

@benchmark f1(v1, v2, ind)

@benchmark f2(v1, v2, ind)

@benchmark f3($v1, $v2, $ind, $res)

@benchmark f4($v1, $v2, $ind, $res)

@benchmark f5($v1, $v2, $ind, $res)

res = zeros(5)

@code_lowered f3(v1, v2, ind, res)

@code_lowered f4(v1, v2, ind, res)


text:
import Lumberjack
dbglog(x::Any) = Lumberjack.debug(string(x))
dbglog(x::Any...) = Lumberjack.debug(string(x))
infolog(x::Any) = Lumberjack.info(string(x))
infolog(x::Any...) = Lumberjack.info(string(x))

formatnum(x::Any; prec=4) = map((x) -> round(x, prec), x)

max_stars_repo_name: zxteloiv/julia-ex

max_issues_repo_name: zxteloiv/julia-ex

max_forks_repo_name: zxteloiv/julia-ex


text:
#=
second-order cone
cuts are (1, -y/||y||_2) for (z, y) outside the cone
TODO extended formulation
=#
function get_init_cuts(cone::MOI.SecondOrderCone)
    cuts = Vector{Float64}[]
    dim = MOI.dimension(cone)

    # z is nonnegative
    cut = zeros(dim)
    cut[1] = 1.0
    push!(cuts, cut)

    # z is at least absolute value of each y_i
    for i in 2:dim
        cut = zeros(dim)
        cut[1] = 1.0
        cut[i] = 1.0
        push!(cuts, cut)
        cut = copy(cut)
        cut[i] = -1.0
        push!(cuts, cut)
    end

    return cuts
end

function get_sep_cuts(x::Vector{Float64}, cone::MOI.SecondOrderCone, tol::Float64)
    dim = MOI.dimension(cone)
    y_norm = norm(x[i] for i in 2:dim)

    if y_norm - x[1] > 0.1 * tol
        cut = similar(x)
        cut[1] = 1
        for i in 2:dim
            cut[i] = -x[i] / y_norm
        end

        if dot(x, cut) < -tol # TODO tolerance option
            return [cut]
        end
    end

    return Vector{Float64}[]
end


text:
"""
    AbstractSchedule

Supertype for objects that schedule `OutputWriter`s and `Diagnostics`.
Schedule must define the functor `Schedule(model)` that returns true or
false.
"""
abstract type AbstractSchedule end

# Default behavior is no alignment.
aligned_time_step(schedule, clock, Δt) = Δt

#####
##### TimeInterval
#####

"""
    struct TimeInterval <: AbstractSchedule

Callable `TimeInterval` schedule for periodic output or diagnostic evaluation
according to `model.clock.time`.
"""
mutable struct TimeInterval <: AbstractSchedule
    interval :: Float64
    previous_actuation_time :: Float64
end

"""
    TimeInterval(interval)

Returns a callable `TimeInterval` that schedules periodic output or diagnostic evaluation
on a `interval` of simulation time, as kept by `model.clock`.
"""
TimeInterval(interval) = TimeInterval(Float64(interval), 0.0)

function (schedule::TimeInterval)(model)
    time = model.clock.time

    if time == schedule.previous_actuation_time + schedule.interval
        schedule.previous_actuation_time = time
        return true
    elseif time > schedule.previous_actuation_time + schedule.interval
        # Shave overshoot off previous_actuation_time to prevent overshoot from accumulating
        schedule.previous_actuation_time = time - rem(time, schedule.interval)
        return true
    else
        return false
    end
end

function aligned_time_step(schedule::TimeInterval, clock, Δt)
    next_actuation_time = schedule.previous_actuation_time + schedule.interval
    return min(Δt, next_actuation_time - clock.time)
end

#####
##### IterationInterval
#####

struct IterationInterval <: AbstractSchedule
    interval :: Int
    offset :: Int
end

"""
    IterationInterval(interval; offset=0)

Returns a callable `IterationInterval` that "actuates" (schedules output or callback execution)
whenever the model iteration (modified by `offset`) is a multiple of `interval`.

For example,

* `IterationInterval(100)` actuates at iterations `[100, 200, 300, ...]`.
* `IterationInterval(100, offset=-1)` actuates at iterations `[99, 199, 299, ...]`.
"""
IterationInterval(interval; offset=0) = IterationInterval(interval, offset)

(schedule::IterationInterval)(model) = (model.clock.iteration - schedule.offset) % schedule.interval == 0

#####
##### WallTimeInterval
#####

mutable struct WallTimeInterval <: AbstractSchedule
    interval :: Float64
    previous_actuation_time :: Float64
end

"""
    WallTimeInterval(interval; start_time = time_ns() * 1e-9)

Returns a callable `WallTimeInterval` that schedules periodic output or diagnostic evaluation
on a `interval` of "wall time" while a simulation runs, in units of seconds.

The "wall time" is the actual real world time in seconds, as kept by an actual
or hypothetical clock hanging on your wall.

The keyword argument `start_time` can be used to specify a starting wall time
other than the moment `WallTimeInterval` is constructed.
"""
WallTimeInterval(interval; start_time = time_ns() * 1e-9) = WallTimeInterval(Float64(interval), Float64(start_time))

function (schedule::WallTimeInterval)(model)
    wall_time = time_ns() * 1e-9

    if wall_time >= schedule.previous_actuation_time + schedule.interval
        # Shave overshoot off previous_actuation_time to prevent overshoot from accumulating
        schedule.previous_actuation_time = wall_time - rem(wall_time, schedule.interval)
        return true
    else
        return false
    end
end

#####
##### SpecifiedTimes
#####

"""
    struct SpecifiedTimes <: AbstractSchedule

Callable `TimeInterval` schedule for periodic output or diagnostic evaluation
according to `model.clock.time`.
"""
mutable struct SpecifiedTimes <: AbstractSchedule
    times :: Vector{Float64}
    previous_actuation :: Int
end

SpecifiedTimes(times::Vararg{<:Number}) = SpecifiedTimes(sort([Float64(t) for t in times]), 0)
SpecifiedTimes(times) = SpecifiedTimes(times...)

function next_appointment_time(st::SpecifiedTimes)
    if st.previous_actuation >= length(st.times)
        return Inf
    else
        return st.times[st.previous_actuation+1]
    end
end

function (st::SpecifiedTimes)(model)
    current_time = model.clock.time

    if current_time >= next_appointment_time(st)
        st.previous_actuation += 1
        return true
    end

    return false
end

align_time_step(schedule::SpecifiedTimes, clock, Δt) = min(Δt, next_appointment_time(schedule) - clock.time)

function specified_times_str(st)
    str_elems = ["$(prettytime(t)), " for t in st.times]
    str_elems = str_elems[1:end-2]
    return string("[", str_elems..., "]")
end

#####
##### ConsecutiveIterations
#####

mutable struct ConsecutiveIterations{S} <: AbstractSchedule
    parent :: S
    consecutive_iterations :: Int
    previous_parent_actuation_iteration :: Int
end

"""
    ConsecutiveIterations(parent_schedule)

Return a `schedule::ConsecutiveIterations` that actuates both when `parent_schedule`
actuates, and at iterations immediately following the actuation of `parent_schedule`.
This can be used, for example, when one wants to use output to reproduce a first-order approximation
of the time derivative of a quantity.
"""
ConsecutiveIterations(parent_schedule, N=1) = ConsecutiveIterations(parent_schedule, N, 0)

function (schedule::ConsecutiveIterations)(model)
    if schedule.parent(model)
        schedule.previous_parent_actuation_iteration = model.clock.iteration
        return true
    elseif model.clock.iteration - schedule.consecutive_iterations <= schedule.previous_parent_actuation_iteration
        return true # The iteration _after_ schedule.parent actuated!
    else
        return false
    end
end

aligned_time_step(schedule::ConsecutiveIterations, clock, Δt) =
    aligned_time_step(schedule.parent, clock. Δt)

#####
##### Any and AndSchedule
#####

struct AndSchedule{S} <: AbstractSchedule
    schedules :: S
    AndSchedule(schedules::S) where S <: Tuple = new{S}(schedules)
end

"""
    AndSchedule(child_schedule_1, child_schedule_2, other_child_schedules...)

Return a schedule that actuates when all `child_schedule`s actuate.
"""
AndSchedule(schedules...) = AndSchedule(Tuple(schedules))

# Note that multiple schedules that have a "state" (like TimeInterval and WallTimeInterval)
# could cause the logic of AndSchedule to fail, due to the short-circuiting nature of `all`.
(as::AndSchedule)(model) = all(schedule(model) for schedule in as.schedules)

struct OrSchedule{S} <: AbstractSchedule
    schedules :: S
    OrSchedule(schedules::S) where S <: Tuple = new{S}(schedules)
end

"""
    OrSchedule(child_schedule_1, child_schedule_2, other_child_schedules...)

Return a schedule that actuates when any of the `child_schedule`s actuates.
"""
OrSchedule(schedules...) = OrSchedule(Tuple(schedules))

function (as::OrSchedule)(model)
    # Ensure that all `schedules` get queried
    actuations = Tuple(schedule(model) for schedule in as.schedules)
    return any(actuations)
end

align_time_step(any_or_all_schedule::Union{OrSchedule, AndSchedule}, clock, Δt) =
    minimum(align_time_step(clock, Δt, schedule) for schedule in any_or_all_schedule.schedules)

#####
##### Show methods
#####

Base.summary(schedule::IterationInterval) = string("IterationInterval(", schedule.interval, ")")
Base.summary(schedule::TimeInterval) = string("TimeInterval(", prettytime(schedule.interval), ")")
Base.summary(schedule::SpecifiedTimes) = string("SpecifiedTimes(", specified_times_str(schedule), ")")
Base.summary(schedule::ConsecutiveIterations) = string("ConsecutiveIterations(",
                                                       summary(schedule.parent), ", ",
                                                       schedule.consecutive_iterations, ")")


text:
export active_blocks, occupancy, launch_configuration

"""
    active_blocks(fun::CuFunction, threads; shmem=0)

Calculate the maximum number of active blocks per multiprocessor when running `threads`
threads of a kernel `fun` requiring `shmem` bytes of dynamic shared memory.
"""
function active_blocks(fun::CuFunction, threads::Integer; shmem::Integer=0)
    blocks_ref = Ref{Cint}()
    cuOccupancyMaxActiveBlocksPerMultiprocessor(blocks_ref, fun, threads, shmem)
    return blocks_ref[]
end

"""
    occupancy(fun::CuFunction, threads; shmem=0)

Calculate the theoretical occupancy of launching `threads` threads of a kernel `fun`
requiring `shmem` bytes of dynamic shared memory.

"""
function occupancy(fun::CuFunction, threads::Integer; shmem::Integer=0)
    # https://devblogs.nvidia.com/cuda-pro-tip-occupancy-api-simplifies-launch-configuration/
    blocks = active_blocks(fun, threads; shmem=shmem)

    mod = fun.mod
    ctx = mod.ctx
    dev = device(ctx)

    threads_per_sm = attribute(dev, DEVICE_ATTRIBUTE_MAX_THREADS_PER_MULTIPROCESSOR)
    warp_size = attribute(dev, DEVICE_ATTRIBUTE_WARP_SIZE)

    return (blocks * threads ÷ warp_size) / (threads_per_sm ÷ warp_size)
end

# HACK: callback function for `launch_configuration` on platforms without support for
#       trampolines as used by `@cfunction` (JuliaLang/julia#27174, JuliaLang/julia#32154)
_shmem_cb = nothing
_shmem_cint_cb(x::Cint) = Cint(something(_shmem_cb)(x))
_shmem_cb_lock = Threads.ReentrantLock()

"""
    launch_configuration(fun::CuFunction; shmem=0, max_threads=0)

Calculate a suggested launch configuration for kernel `fun` requiring `shmem` bytes of
dynamic shared memory. Returns a tuple with a suggested amount of threads, and the minimal
amount of blocks to reach maximal occupancy. Optionally, the maximum amount of threads can
be constrained using `max_threads`.

In the case of a variable amount of shared memory, pass a callable object for `shmem`
instead, taking a single integer representing the block size and returning the amount of
dynamic shared memory for that configuration.
"""
function launch_configuration(fun::CuFunction; shmem::Union{Integer,Base.Callable}=0,
                              max_threads::Integer=0)
    blocks_ref = Ref{Cint}()
    threads_ref = Ref{Cint}()
    if isa(shmem, Integer)
        cuOccupancyMaxPotentialBlockSize(blocks_ref, threads_ref, fun, C_NULL, shmem, max_threads)
    elseif Sys.ARCH == :x86 || Sys.ARCH == :x86_64
        shmem_cint = threads -> Cint(shmem(threads))
        cb = @cfunction($shmem_cint, Cint, (Cint,))
        cuOccupancyMaxPotentialBlockSize(blocks_ref, threads_ref, fun, cb, 0, max_threads)
    else
        lock(_shmem_cb_lock) do
            global _shmem_cb
            _shmem_cb = shmem
            cb = @cfunction(_shmem_cint_cb, Cint, (Cint,))
            cuOccupancyMaxPotentialBlockSize(blocks_ref, threads_ref, fun, cb, 0, max_threads)
            _shmem_cb = nothing
        end
    end
    return (blocks=Int(blocks_ref[]), threads=Int(threads_ref[]))
end


text:
using DDD, Test, SparseArrays, LinearAlgebra, StaticArrays
cd(@__DIR__)

@testset "Integral tests" begin
    dlnParams = DislocationParameters(; mobility = mobBCC())
    matParams = MaterialParameters(; crystalStruct = BCC())
    femParams = FEMParameters(;
        type = DispatchRegularCuboidMesh(),
        order = LinearElement(),
        model = CantileverLoad(),
        dx = 1013.0,
        dy = 1987.0,
        dz = 2999.0,
        mx = 19,
        my = 23,
        mz = 29,
    )
    intParams = IntegrationParameters(; method = AdaptiveEulerTrapezoid())
    slipSystem = SlipSystem(;
        crystalStruct = BCC(),
        slipPlane = Float64[-1; 1; 0],
        bVec = Float64[1; 1; 1],
    )
    dx, dy, dz = femParams.dx, femParams.dy, femParams.dz
    segLen = (dx + dy + dz) / 30
    prismLoop = DislocationLoop(;
        loopType = loopPrism(),    # Prismatic loop, all segments are edge segments.
        numSides = 8,   # 8-sided loop.
        nodeSide = 1,   # One node per side, if 1 nodes will be in the corners.
        numLoops = 1,   # Number of loops of this type to generate when making a network.
        segLen = segLen * SVector{8}(ones(8)),  # Length of each segment between nodes, equal to the number of nodes.
        slipSystemIdx = 1, # Slip System (assuming slip systems are stored in a file, this is the index).
        slipSystem = slipSystem,  # Slip plane of the segments.
        label = SVector{8, nodeTypeDln}(1, 1, 1, 1, 1, 1, 1, 1),    # Node labels, has to be equal to the number of nodes.
        buffer = 0,   # Buffer to increase the dislocation spread.
        range = SMatrix{3, 2, Float64}(dx / 2, dy / 2, dz / 2, dx / 2, dy / 2, dz / 2),  # Distribution range
        dist = Zeros(),  # Loop distribution.
    )
    shearLoop = DislocationLoop(;
        loopType = loopShear(),    # Prismatic loop, all segments are edge segments.
        numSides = 8,   # 8-sided loop.
        nodeSide = 1,   # One node per side, if 1 nodes will be in the corners.
        numLoops = 1,   # Number of loops of this type to generate when making a network.
        segLen = segLen * SVector{8}(ones(8)),  # Length of each segment between nodes, equal to the number of nodes.
        slipSystemIdx = 1, # Slip System (assuming slip systems are stored in a file, this is the index).
        slipSystem = slipSystem,  # Slip plane of the segments.
        label = SVector{8, nodeTypeDln}(1, 1, 1, 1, 1, 1, 1, 1),    # Node labels, has to be equal to the number of nodes.
        buffer = 0,   # Buffer to increase the dislocation spread.
        range = SMatrix{3, 2, Float64}(dx / 2, dy / 2, dz / 2, dx / 2, dy / 2, dz / 2),  # Distribution range
        dist = Zeros(),  # Loop distribution.
    )
    network = DislocationNetwork([prismLoop, shearLoop])
    regularCuboidMesh = buildMesh(matParams, femParams)
    cantileverBC, forceDisplacement = Boundaries(femParams, regularCuboidMesh)

    coordFE = regularCuboidMesh.coord
    uGamma = cantileverBC.uGamma.node
    uDofsDln = cantileverBC.uDofsDln
    numSeg = network.numSeg[1]

    segForce =
        calcSegForce(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network)
    calcSegForce!(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network)
    @test isapprox(segForce, network.segForce[:, :, 1:numSeg])

    idx = 1:2:16
    segForceIdx = calcSegForce(
        dlnParams,
        matParams,
        regularCuboidMesh,
        forceDisplacement,
        network,
        idx,
    )
    calcSegForce!(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network, idx)
    @test isapprox(segForceIdx, network.segForce[:, :, idx])

    σTilde = calc_σTilde(coordFE[:, uGamma], dlnParams, matParams, network)
    σTilde2 = zeros(size(σTilde))
    calc_σTilde!(σTilde2, coordFE[:, uGamma], dlnParams, matParams, network)
    @test isapprox(σTilde, σTilde2)

    uTilde = calc_uTilde(regularCuboidMesh, cantileverBC, matParams, network)
    calc_uTilde!(forceDisplacement, regularCuboidMesh, cantileverBC, matParams, network)
    @test isapprox(uTilde, forceDisplacement.uTilde[uDofsDln])

    idx = [
        5771
        6852
        7419
        7436
        6903
        5822
        5255
        5238
        7442
        6852
        6298
        5213
        5232
        5822
        6376
        7461
    ]
    uHatIdx = unique(regularCuboidMesh.connectivity[:, idx]) * 3
    uHatDofs = [uHatIdx .- 2; uHatIdx .- 1; uHatIdx]
    randUHat = sprand(length(uHatDofs), 0.5)
    randIdx = findall(!iszero, randUHat)
    randUHat[randIdx] .= 1:(length(randIdx) * 0.01)
    forceDisplacement.uHat[uHatDofs] = randUHat

    segForce =
        calcSegForce(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network)
    calcSegForce!(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network)
    @test isapprox(segForce, network.segForce[:, :, 1:numSeg])

    idx = 1:2:16
    segForceIdx = calcSegForce(
        dlnParams,
        matParams,
        regularCuboidMesh,
        forceDisplacement,
        network,
        idx,
    )
    calcSegForce!(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network, idx)
    @test isapprox(segForceIdx, network.segForce[:, :, idx])

    network2 = deepcopy(network)
    network2 = remeshSurfaceNetwork!(regularCuboidMesh, cantileverBC, network2)
    getSegmentIdx!(network2)
    @test compStruct(network2, network)

    prismLoopIntersect = DislocationLoop(;
        loopType = loopShear(),    # Prismatic loop, all segments are edge segments.
        numSides = 8,   # 8-sided loop.
        nodeSide = 1,   # One node per side, if 1 nodes will be in the corners.
        numLoops = 1,   # Number of loops of this type to generate when making a network.
        segLen = segLen * SVector{8}(ones(8)),  # Length of each segment between nodes, equal to the number of nodes.
        slipSystemIdx = 1, # Slip System (assuming slip systems are stored in a file, this is the index).
        slipSystem = slipSystem,  # Slip plane of the segments.
        label = SVector{8, nodeTypeDln}(1, 1, 1, 1, 1, 1, 1, 1),    # Node labels, has to be equal to the number of nodes.
        buffer = 0,   # Buffer to increase the dislocation spread.
        range = SMatrix{3, 2, Float64}(
            segLen / 2,
            segLen / 2,
            segLen / 2,
            segLen / 2,
            segLen / 2,
            segLen / 2,
        ),  # Distribution range
        dist = Zeros(),  # Loop distribution.
    )
    DislocationNetwork!(network2, prismLoopIntersect)
    numNode = network2.numNode[1]

    network2.nodeVel[1, 1:numNode] .= [
        0.848409575718539,
        0.6652242160082342,
        0.25113206533518206,
        0.17655640207336298,
        0.33256805603630246,
        0.9969648349174243,
        0.12851456965799923,
        0.05591976789414588,
        0.29757306575603804,
        0.2723507312786162,
        0.8874071844807814,
        0.4429058105879311,
        0.8868742664817155,
        0.4598309743323694,
        0.03931969223629794,
        0.12373984299000429,
        0.03861761480269377,
        0.7603101804979915,
        0.4543587651493177,
        0.27460243590637345,
        0.3,
        0.5434785594380296,
        0.5395078480115256,
        0.8654239757909008,
    ]
    network2.nodeVel[2, 1:numNode] .= [
        0.8636080299591222,
        0.37959594965548926,
        0.9920015465135401,
        0.702809412405236,
        0.14450700182749543,
        0.4816990411158426,
        0.2905765122000159,
        0.319012344824527,
        0.05952526906892097,
        0.4392216093355392,
        0.8752415268798401,
        0.3592278150943291,
        0.8529140139462188,
        0.0993328662955082,
        0.14584336379153573,
        0.43739168117035354,
        0.3338215613113067,
        0.08518695747920568,
        0.7462534959446294,
        0.5272909447878065,
        0.5,
        0.7878515928217127,
        0.5106406121544798,
        0.22360673480968973,
    ]
    network2.nodeVel[3, 1:numNode] .= [
        0.38025510658061124,
        0.028912486176383645,
        0.24395899994533154,
        0.8387909328083565,
        0.041107268905112626,
        0.6138145136425359,
        0.5325226796794968,
        0.21895716122452558,
        0.27194053775758875,
        0.6585257914167739,
        0.6004564313623306,
        0.10944691655246075,
        0.7382692437252569,
        0.9260924108803184,
        0.46106334186849685,
        0.5657795344686234,
        0.09638787333766241,
        0.2854907820375585,
        0.04887167777640866,
        0.850652937369534,
        0.2,
        0.7285145955076195,
        0.39460521446188235,
        0.0044684000305827976,
    ]

    network2 = remeshSurfaceNetwork!(regularCuboidMesh, cantileverBC, network2)
    getSegmentIdx!(network2)

    network2.numNode[1] == numNode + 2
    label = network2.label
    ext = findall(x -> x == 5, label)
    surf = findall(x -> x == 3, label)
    @test length(ext) == 3
    @test length(surf) == 1
    @test compStruct(network2, network) == false

    numSeg = network2.numSeg[1]
    segForce =
        calcSegForce(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network2)
    calcSegForce!(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network2)
    @test isapprox(segForce, network2.segForce[:, :, 1:numSeg])

    idx = 1:2:16
    segForceIdx = calcSegForce(
        dlnParams,
        matParams,
        regularCuboidMesh,
        forceDisplacement,
        network2,
        idx,
    )
    calcSegForce!(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network2, idx)
    @test isapprox(segForceIdx, network2.segForce[:, :, idx])

    σTilde = calc_σTilde(coordFE[:, uGamma], dlnParams, matParams, network2)
    σTilde2 = zeros(size(σTilde))
    calc_σTilde!(σTilde2, coordFE[:, uGamma], dlnParams, matParams, network2)
    @test isapprox(σTilde, σTilde2)

    uTilde = calc_uTilde(regularCuboidMesh, cantileverBC, matParams, network2)
    calc_uTilde!(forceDisplacement, regularCuboidMesh, cantileverBC, matParams, network2)
    @test isapprox(uTilde, forceDisplacement.uTilde[uDofsDln])

    idx = [
        5771
        6852
        7419
        7436
        6903
        5822
        5255
        5238
        7442
        6852
        6298
        5213
        5232
        5822
        6376
        7461
    ]
    uHatIdx = unique(regularCuboidMesh.connectivity[:, idx]) * 3
    uHatDofs = [uHatIdx .- 2; uHatIdx .- 1; uHatIdx]
    randUHat = sprand(length(uHatDofs), 0.5)
    randIdx = findall(!iszero, randUHat)
    randUHat[randIdx] .= 1:(length(randIdx) * 0.01)
    forceDisplacement.uHat[uHatDofs] = randUHat

    segForce =
        calcSegForce(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network2)
    calcSegForce!(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network2)
    @test isapprox(segForce, network2.segForce[:, :, 1:numSeg])

    idx = 1:2:16
    segForceIdx = calcSegForce(
        dlnParams,
        matParams,
        regularCuboidMesh,
        forceDisplacement,
        network2,
        idx,
    )
    calcSegForce!(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network2, idx)
    @test isapprox(segForceIdx, network2.segForce[:, :, idx])

    network3 = deepcopy(network2)
    network3 = remeshSurfaceNetwork!(regularCuboidMesh, cantileverBC, network3)
    getSegmentIdx!(network3)

    compStruct(network3, network2)
    coarsenVirtualNetwork!(dlnParams, network3)
    @test compStruct(network3, network2)
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
using Pkg, Scratch, Preferences

# Use `cc` to build `sandbox.c` into a scratch space owned by `Sandbox`
sandbox_uuid = Base.UUID("9307e30f-c43e-9ca7-d17c-c2dc59df670d")
sdir = get_scratch!(sandbox_uuid, "local_sandbox")
sandbox_path = joinpath(sdir, "sandbox")
sandbox_src = joinpath(@__DIR__, "userns_sandbox.c")
run(`cc -std=c99 -O2 -static -static-libgcc -g -o $(sandbox_path) $(sandbox_src)`)

# Tell UserNSSandbox_jll to load our `sandbox` instead of the default artifact one
set_preferences!(
    joinpath(dirname(@__DIR__), "LocalPreferences.toml"),
    "UserNSSandbox_jll",
    "sandbox_path" => sandbox_path;
    force=true,
)


text:
# create a new "build_series_args" which converts all inputs into xs = Any[xitems], ys = Any[yitems].
# Special handling for: no args, xmin/xmax, parametric, dataframes
# Then once inputs have been converted, build the series args, map functions, etc.
# This should cut down on boilerplate code and allow more focused dispatch on type
# note: returns meta information... mainly for use with automatic labeling from DataFrames for now

typealias FuncOrFuncs Union{Function, AVec{Function}}

all3D(d::KW) = trueOrAllTrue(st -> st in (:contour, :contourf, :heatmap, :surface, :wireframe, :contour3d, :image), get(d, :seriestype, :none))

# missing
convertToAnyVector(v::Void, d::KW) = Any[nothing], nothing

# fixed number of blank series
convertToAnyVector(n::Integer, d::KW) = Any[zeros(0) for i in 1:n], nothing

# numeric vector
convertToAnyVector{T<:Number}(v::AVec{T}, d::KW) = Any[v], nothing

# string vector
convertToAnyVector{T<:AbstractString}(v::AVec{T}, d::KW) = Any[v], nothing

function convertToAnyVector(v::AMat, d::KW)
    if all3D(d)
        Any[Surface(v)]
    else
        Any[v[:,i] for i in 1:size(v,2)]
    end, nothing
end

# function
convertToAnyVector(f::Function, d::KW) = Any[f], nothing

# surface
convertToAnyVector(s::Surface, d::KW) = Any[s], nothing

# volume
convertToAnyVector(v::Volume, d::KW) = Any[v], nothing

# # vector of OHLC
# convertToAnyVector(v::AVec{OHLC}, d::KW) = Any[v], nothing

# # dates
# convertToAnyVector{D<:Union{Date,DateTime}}(dts::AVec{D}, d::KW) = Any[dts], nothing

# list of things (maybe other vectors, functions, or something else)
function convertToAnyVector(v::AVec, d::KW)
    if all(x -> typeof(x) <: Number, v)
        # all real numbers wrap the whole vector as one item
        Any[convert(Vector{Float64}, v)], nothing
    else
        # something else... treat each element as an item
        vcat(Any[convertToAnyVector(vi, d)[1] for vi in v]...), nothing
        # Any[vi for vi in v], nothing
    end
end

convertToAnyVector(t::Tuple, d::KW) = Any[t], nothing

function convertToAnyVector(args...)
    error("In convertToAnyVector, could not handle the argument types: $(map(typeof, args[1:end-1]))")
end

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

# TODO: can we avoid the copy here?  one error that crops up is that mapping functions over the same array
#       result in that array being shared.  push!, etc will add too many items to that array

compute_x(x::Void, y::Void, z)      = 1:size(z,1)
compute_x(x::Void, y, z)            = 1:size(y,1)
compute_x(x::Function, y, z)        = map(x, y)
compute_x(x, y, z)                  = copy(x)

# compute_y(x::Void, y::Function, z)  = error()
compute_y(x::Void, y::Void, z)      = 1:size(z,2)
compute_y(x, y::Function, z)        = map(y, x)
compute_y(x, y, z)                  = copy(y)

compute_z(x, y, z::Function)        = map(z, x, y)
compute_z(x, y, z::AbstractMatrix)  = Surface(z)
compute_z(x, y, z::Void)            = nothing
compute_z(x, y, z)                  = copy(z)

nobigs(v::AVec{BigFloat}) = map(Float64, v)
nobigs(v::AVec{BigInt}) = map(Int64, v)
nobigs(v) = v

@noinline function compute_xyz(x, y, z)
    x = compute_x(x,y,z)
    y = compute_y(x,y,z)
    z = compute_z(x,y,z)
    nobigs(x), nobigs(y), nobigs(z)
end

# not allowed
compute_xyz(x::Void, y::FuncOrFuncs, z)       = error("If you want to plot the function `$y`, you need to define the x values!")
compute_xyz(x::Void, y::Void, z::FuncOrFuncs) = error("If you want to plot the function `$z`, you need to define x and y values!")
compute_xyz(x::Void, y::Void, z::Void)        = error("x/y/z are all nothing!")

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

# we are going to build recipes to do the processing and splitting of the args

# ensure we dispatch to the slicer
immutable SliceIt end

# the catch-all recipes
@recipe function f(::Type{SliceIt}, x, y, z)

    # handle data with formatting attached
    if typeof(x) <: Formatted
        xformatter := x.formatter
        x = x.data
    end
    if typeof(y) <: Formatted
        yformatter := y.formatter
        y = y.data
    end
    if typeof(z) <: Formatted
        zformatter := z.formatter
        z = z.data
    end

    xs, _ = convertToAnyVector(x, d)
    ys, _ = convertToAnyVector(y, d)
    zs, _ = convertToAnyVector(z, d)

    fr = pop!(d, :fillrange, nothing)
    fillranges, _ = if typeof(fr) <: Number
        ([fr],nothing)
    else
        convertToAnyVector(fr, d)
    end
    mf = length(fillranges)

    # @show zs

    mx = length(xs)
    my = length(ys)
    mz = length(zs)
    if mx > 0 && my > 0 && mz > 0
        for i in 1:max(mx, my, mz)
            # add a new series
            di = copy(d)
            xi, yi, zi = xs[mod1(i,mx)], ys[mod1(i,my)], zs[mod1(i,mz)]
            di[:x], di[:y], di[:z] = compute_xyz(xi, yi, zi)

            # handle fillrange
            fr = fillranges[mod1(i,mf)]
            di[:fillrange] = isa(fr, Function) ? map(fr, di[:x]) : fr

            push!(series_list, RecipeData(di, ()))
        end
    end
    nothing  # don't add a series for the main block
end

# this is the default "type recipe"... just pass the object through
@recipe f{T<:Any}(::Type{T}, v::T) = v

# this should catch unhandled "series recipes" and error with a nice message
@recipe f{V<:Val}(::Type{V}, x, y, z) = error("The backend must not support the series type $V, and there isn't a series recipe defined.")

_apply_type_recipe(d, v) = RecipesBase.apply_recipe(d, typeof(v), v)[1].args[1]

# Handle type recipes when the recipe is defined on the elements.
# This sort of recipe should return a pair of functions... one to convert to number,
# and one to format tick values.
function _apply_type_recipe(d, v::AbstractArray)
    args = RecipesBase.apply_recipe(d, typeof(v[1]), v[1])[1].args
    if length(args) == 2 && typeof(args[1]) <: Function && typeof(args[2]) <: Function
        numfunc, formatter = args
        Formatted(map(numfunc, v), formatter)
    else
        v
    end
end

# # special handling for Surface... need to properly unwrap and re-wrap
# function _apply_type_recipe(d, v::Surface)
#     T = eltype(v.surf)
#     @show T
#     if T <: Integer || T <: AbstractFloat
#         v
#     else
#         ret = _apply_type_recipe(d, v.surf)
#         if typeof(ret) <: Formatted
#             Formatted(Surface(ret.data), ret.formatter)
#         else
#             v
#         end
#     end
# end

# don't do anything for ints or floats
_apply_type_recipe{T<:Union{Integer,AbstractFloat}}(d, v::AbstractArray{T}) = v

# handle "type recipes" by converting inputs, and then either re-calling or slicing
@recipe function f(x, y, z)
    did_replace = false
    newx = _apply_type_recipe(d, x)
    x === newx || (did_replace = true)
    newy = _apply_type_recipe(d, y)
    y === newy || (did_replace = true)
    newz = _apply_type_recipe(d, z)
    z === newz || (did_replace = true)
    if did_replace
        newx, newy, newz
    else
        SliceIt, x, y, z
    end
end
@recipe function f(x, y)
    did_replace = false
    newx = _apply_type_recipe(d, x)
    x === newx || (did_replace = true)
    newy = _apply_type_recipe(d, y)
    y === newy || (did_replace = true)
    if did_replace
        newx, newy
    else
        SliceIt, x, y, nothing
    end
end
@recipe function f(y)
    newy = _apply_type_recipe(d, y)
    if y !== newy
        newy
    else
        SliceIt, nothing, y, nothing
    end
end

# if there's more than 3 inputs, it can't be passed directly to SliceIt
# so we'll apply_type_recipe to all of them
@recipe function f(v1, v2, v3, v4, vrest...)
    did_replace = false
    newargs = map(v -> begin
        newv = _apply_type_recipe(d, v)
        if newv !== v
            did_replace = true
        end
        newv
    end, (v1, v2, v3, v4, vrest...))
    if !did_replace
        error("Couldn't process recipe args: $(map(typeof, (v1, v2, v3, v4, vrest...)))")
    end
    newargs
end

# # --------------------------------------------------------------------
# # 1 argument
# # --------------------------------------------------------------------

@recipe f(n::Integer) = is3d(get(d,:seriestype,:path)) ? (SliceIt, n, n, n) : (SliceIt, n, n, nothing)

# return a surface if this is a 3d plot, otherwise let it be sliced up
@recipe function f{T<:Union{Integer,AbstractFloat}}(mat::AMat{T})
    if all3D(d)
        n,m = size(mat)
        SliceIt, 1:m, 1:n, Surface(mat)
    else
        SliceIt, nothing, mat, nothing
    end
end

# if a matrix is wrapped by Formatted, do similar logic, but wrap data with Surface
@recipe function f{T<:AbstractMatrix}(fmt::Formatted{T})
    if all3D(d)
        mat = fmt.data
        n,m = size(mat)
        SliceIt, 1:m, 1:n, Formatted(Surface(mat), fmt.formatter)
    else
        SliceIt, nothing, fmt, nothing
    end
end

# assume this is a Volume, so construct one
@recipe function f{T<:Number}(vol::AbstractArray{T,3}, args...)
    seriestype := :volume
    SliceIt, nothing, Volume(vol, args...), nothing
end

# # images - grays

@recipe function f{T<:Gray}(mat::AMat{T})
    if is_seriestype_supported(:image)
        seriestype := :image
        n, m = size(mat)
        SliceIt, 1:m, 1:n, Surface(mat)
    else
        seriestype := :heatmap
        yflip --> true
        fillcolor --> ColorGradient([:black, :white])
        SliceIt, 1:m, 1:n, Surface(convert(Matrix{Float64}, mat))
    end
end

# # images - colors

@recipe function f{T<:Colorant}(mat::AMat{T})
    if is_seriestype_supported(:image)
        seriestype := :image
        n, m = size(mat)
        SliceIt, 1:m, 1:n, Surface(mat)
    else
        seriestype := :heatmap
        yflip --> true
        z, d[:fillcolor] = replace_image_with_heatmap(mat)
        SliceIt, 1:m, 1:n, Surface(z)
    end
end

#
# # plotting arbitrary shapes/polygons

@recipe function f(shape::Shape)
    seriestype --> :shape
    shape_coords(shape)
end

@recipe function f(shapes::AVec{Shape})
    seriestype --> :shape
    shape_coords(shapes)
end

@recipe function f(shapes::AMat{Shape})
    seriestype --> :shape
    for j in 1:size(shapes,2)
        @series shape_coords(vec(shapes[:,j]))
    end
end

# function without range... use the current range of the x-axis

@recipe function f(f::FuncOrFuncs)
    plt = d[:plot_object]
    xmin,xmax = axis_limits(plt[1][:xaxis])
    f, xmin, xmax
end

#
# # --------------------------------------------------------------------
# # 2 arguments
# # --------------------------------------------------------------------
#
#
# # if functions come first, just swap the order (not to be confused with parametric functions...
# # as there would be more than one function passed in)

@recipe function f(f::FuncOrFuncs, x)
    @assert !(typeof(x) <: FuncOrFuncs)  # otherwise we'd hit infinite recursion here
    x, f
end

#
# # --------------------------------------------------------------------
# # 3 arguments
# # --------------------------------------------------------------------
#
#
# # 3d line or scatter

@recipe function f(x::AVec, y::AVec, z::AVec)
    # st = get(d, :seriestype, :none)
    # if st == :scatter
    #     d[:seriestype] = :scatter3d
    # elseif !is3d(st)
    #     d[:seriestype] = :path3d
    # end
    SliceIt, x, y, z
end

@recipe function f(x::AMat, y::AMat, z::AMat)
    # st = get(d, :seriestype, :none)
    # if size(x) == size(y) == size(z)
    #     if !is3d(st)
    #         seriestype := :path3d
    #     end
    # end
    SliceIt, x, y, z
end

#
# # surface-like... function

@recipe function f(x::AVec, y::AVec, zf::Function)
    # x = X <: Number ? sort(x) : x
    # y = Y <: Number ? sort(y) : y
    SliceIt, x, y, Surface(zf, x, y)  # TODO: replace with SurfaceFunction when supported
end

#
# # surface-like... matrix grid

@recipe function f(x::AVec, y::AVec, z::AMat)
    if !like_surface(get(d, :seriestype, :none))
        d[:seriestype] = :contour
    end
    SliceIt, x, y, Surface(z)
end

#
#
# # --------------------------------------------------------------------
# # Parametric functions
# # --------------------------------------------------------------------

#
# # special handling... xmin/xmax with parametric function(s)
@recipe f(f::FuncOrFuncs, xmin::Number, xmax::Number) = linspace(xmin, xmax, 100), f
@recipe f(fx::FuncOrFuncs, fy::FuncOrFuncs, u::AVec)  = mapFuncOrFuncs(fx, u), mapFuncOrFuncs(fy, u)
@recipe f(fx::FuncOrFuncs, fy::FuncOrFuncs, umin::Number, umax::Number, n = 200) = fx, fy, linspace(umin, umax, n)

#
# # special handling... 3D parametric function(s)
@recipe function f(fx::FuncOrFuncs, fy::FuncOrFuncs, fz::FuncOrFuncs, u::AVec)
    mapFuncOrFuncs(fx, u), mapFuncOrFuncs(fy, u), mapFuncOrFuncs(fz, u)
end
@recipe function f(fx::FuncOrFuncs, fy::FuncOrFuncs, fz::FuncOrFuncs, umin::Number, umax::Number, numPoints = 200)
    fx, fy, fz, linspace(umin, umax, numPoints)
end

#
#
# # --------------------------------------------------------------------
# # Lists of tuples and FixedSizeArrays
# # --------------------------------------------------------------------
#
# # if we get an unhandled tuple, just splat it in
@recipe f(tup::Tuple) = tup

#
# # (x,y) tuples
@recipe f{R1<:Number,R2<:Number}(xy::AVec{Tuple{R1,R2}}) = unzip(xy)
@recipe f{R1<:Number,R2<:Number}(xy::Tuple{R1,R2})       = [xy[1]], [xy[2]]

#
# # (x,y,z) tuples
@recipe f{R1<:Number,R2<:Number,R3<:Number}(xyz::AVec{Tuple{R1,R2,R3}}) = unzip(xyz)
@recipe f{R1<:Number,R2<:Number,R3<:Number}(xyz::Tuple{R1,R2,R3})       = [xyz[1]], [xyz[2]], [xyz[3]]

# these might be points+velocity, or OHLC or something else
@recipe f{R1<:Number,R2<:Number,R3<:Number,R4<:Number}(xyuv::AVec{Tuple{R1,R2,R3,R4}}) = get(d,:seriestype,:path)==:ohlc ? OHLC[OHLC(t...) for t in xyuv] : unzip(xyuv)
@recipe f{R1<:Number,R2<:Number,R3<:Number,R4<:Number}(xyuv::Tuple{R1,R2,R3,R4})       = [xyuv[1]], [xyuv[2]], [xyuv[3]], [xyuv[4]]

#
# # 2D FixedSizeArrays
@recipe f{T<:Number}(xy::AVec{FixedSizeArrays.Vec{2,T}}) = unzip(xy)
@recipe f{T<:Number}(xy::FixedSizeArrays.Vec{2,T})       = [xy[1]], [xy[2]]

#
# # 3D FixedSizeArrays
@recipe f{T<:Number}(xyz::AVec{FixedSizeArrays.Vec{3,T}}) = unzip(xyz)
@recipe f{T<:Number}(xyz::FixedSizeArrays.Vec{3,T})       = [xyz[1]], [xyz[2]], [xyz[3]]

#
# # --------------------------------------------------------------------
# # handle grouping
# # --------------------------------------------------------------------

# @recipe function f(groupby::GroupBy, args...)
#     for (i,glab) in enumerate(groupby.groupLabels)
#         # create a new series, with the label of the group, and an idxfilter (to be applied in slice_and_dice)
#         # TODO: use @series instead
#         @show i, glab, groupby.groupIds[i]
#         di = copy(d)
#         get!(di, :label, string(glab))
#         get!(di, :idxfilter, groupby.groupIds[i])
#         push!(series_list, RecipeData(di, args))
#     end
#     nothing
# end

# split the group into 1 series per group, and set the label and idxfilter for each
@recipe function f(groupby::GroupBy, args...)
    for (i,glab) in enumerate(groupby.groupLabels)
        @series begin
            label     --> string(glab)
            idxfilter --> groupby.groupIds[i]
            args
        end
    end
end


text:
# General helper functions, That implementation used in julia6 with Knet 0.8.4

# KnetArray cleaning
function free_KnetArray()
    gc(); Knet.knetgc(); gc()
end

# To fill the constant word embeddings
function fillwvecs!(sentences, isents, wembed; GPUFEATS=false)
    for (s, isents) in zip(sentences, isents)
        empty!(s.wvec)
        for w in isents
            if GPUFEATS
                push!(s.wvec, wembed[:, w])
            else
                push!(s.wvec, Array(wembed[:, w]))
            end
        end
    end
end

# To fill the context embeddings of words in given sentences one by one
function fillcvecs!(sentences, forw, back; GPUFEATS=false)
    T = length(forw)
    for i in 1:length(sentences)
        s = sentences[i]
        empty!(s.fvec)
        empty!(s.bvec)
        N = length(s)
        for n in 1:N
            t = T-N+n
            if GPUFEATS #GPU
                push!(s.fvec, forw[t][:,i])
                push!(s.bvec, back[t][:,i])
            else #CPU
                push!(s.fvec, Array(forw[t][:,i]))
                push!(s.bvec, Array(back[t][:,i]))
            end
        end
    end
end

# fill both word and context embeddings
function fillvecs!(wmodel, sentences, vocab; batchsize=128)

    words, sents, maxwordlen, maxsentlen = maptoint(sentences, vocab)
    sow = vocab.cdict[vocab.sowchar]
    eow = vocab.cdict[vocab.eowchar]
    paw = vocab.cdict[vocab.unkchar]

    # word-embeddings calcutation
    wembed = Any[]
    free_KnetArray();
    for i=1:batchsize:length(words)
        j = min(i+batchsize-1,length(words))
        wij = view(words,i:j)
        maxij = maximum(map(length, wij))
        cdata, cmask = tokenbatch(wij, maxij, sow, eow)
        push!(wembed, charlstm(wmodel, cdata, cmask))
    end
    wembed =hcat(wembed...) # Here I applied the changes from hcatn-> hcat in newer version
    fillwvecs!(sentences, sents, wembed)

    sos,eos,unk = vocab.idict[vocab.sosword], vocab.idict[vocab.eosword], vocab.odict[vocab.unkword]
    result = zeros(2)
    free_KnetArray()
    for i=1:batchsize:length(sents)
        j = min(i+batchsize-1, length(sents))
        isentij = view(sents, i:j)
        maxij = maximum(map(length, isentij))
        wdata, wmask = tokenbatch(isentij, maxij, sos, eos)
        forw, back = wordlstm(wmodel, wdata, wmask, wembed)
        sentij = view(sentences, i:j)
        fillcvecs!(sentij, forw, back)
        odata, omask = goldbatch(sentij, maxij, vocab.odict, unk)
        lmloss(wmodel,odata,omask,forw,back; result=result)
    end
    return exp(-result[1]/result[2])
end

# calculate the coverage
function unkrate(sentences)
    total = known = 0
    for s in sentences
        for w in s.word
            total += 1
            known += haskey(s.vocab.odict, w)
        end
    end
    return (total-known,total)
end

# convert lm-trained file to new col-major parser style file
function convertfile(infile, outfile)
    atr(x)=transpose(Array(x)) # LM stores in row-major Array, we use col-major Array
    d = load(infile); m = d["model"]
    save(outfile, "cembed", atr(m[:cembed]), "forw", map(atr,m[:forw]),
         "back",map(atr,m[:back]), "soft",map(atr,m[:soft]), "char",map(atr,m[:char]),
         "char_vocab",d["char_vocab"], "word_vocab",d["word_vocab"],
         "sosword","<s>","eosword","</s>","unkword","<unk>",
         "sowchar",'\x12',"eowchar",'\x13',"unkchar",'\x11')
end

# convert KnetArrays to Arrays
function convert2cpu(infile, outfile)
    d = load(infile)
    jldopen(outfile, "w") do file
        for (k,v) in d
            write(file,k,map2cpu(v))
        end
    end
end

# save file in new format
function savefile(file, vocab, wmodel, pmodel, optim, arctype, feats)
    save(file,  "cembed", map2cpu(cembed(wmodel)),
         "char", map2cpu([wchar(wmodel),bchar(wmodel)]),
         "forw", map2cpu([wforw(wmodel),bforw(wmodel)]),
         "back", map2cpu([wback(wmodel),bback(wmodel)]),
         "soft", map2cpu([wsoft(wmodel),bsoft(wmodel)]),

         "char_vocab",vocab.cdict, "word_vocab",vocab.odict,
         "sosword",vocab.sosword,"eosword",vocab.eosword,"unkword",vocab.unkword,
         "sowchar",vocab.sowchar,"eowchar",vocab.eowchar,"unkchar",vocab.unkchar,
         "postags",vocab.postags,"deprels",vocab.deprels,"xpostags",vocab.xpostags,

         "postagv",map2cpu(postagv(pmodel)),"deprelv",map2cpu(deprelv(pmodel)),
         "lcountv",map2cpu(lcountv(pmodel)),"rcountv",map2cpu(rcountv(pmodel)),
         "distancev",map2cpu(distancev(pmodel)),"parserv",map2cpu(parserv(pmodel)),
         "xpostagv",map2cpu(xpostagv(pmodel)),

         "postago",map2cpu(postagv(optim)),"deprelo",map2cpu(deprelv(optim)),
         "lcounto",map2cpu(lcountv(optim)),"rcounto",map2cpu(rcountv(optim)),
         "distanceo",map2cpu(distancev(optim)),"parsero",map2cpu(parserv(optim)),"xpostago",map2cpu(xpostagv(optim)),

         "arctype",arctype,"feats",feats,
    )
end

function writeconllu(sentences, inputfile, outputfile)
    # We only replace the head and deprel fields of the input file
    out = open(outputfile,"w")
    v = sentences[1].vocab
    deprels = Array(String, length(v.deprels))
    for (k,v) in v.deprels; deprels[v]=k; end
    s = p = nothing
    ns = nw = nl = 0
    for line in eachline(inputfile)
        nl += 1
        if ismatch(r"^\d+\t", line)
            # info("$nl word")
            if s == nothing
                s = sentences[ns+1]
                p = s.parse
            end
            f = split(line, '\t');
            nw += 1
            if f[1] != "$nw"; error(); end;
            if f[2] != s.word[nw]; error(); end
            f[7] = string(p.head[nw])
            f[8] = deprels[p.deprel[nw]]
            print(out, join(f, "\t"))
            print(out, "\n")
        else
            if line == "" # "\n" in Julia 6 introduce that bug!
                # info("$nl blank")
                if s == nothing; error(); end
                if nw != length(s.word); error(); end
                ns += 1; nw = 0
                s = p = nothing
            else
                # info("$nl non-word")
            end
            print(out, line)
            print(out, "\n")
        end
    end
    if ns != length(sentences); error(); end
    close(out)
end

map2cpu(x)=(if isbits(x); x; else; map2cpu2(x); end)
map2cpu(x::KnetArray)=Array(x)
map2cpu(x::Tuple)=map(map2cpu,x)
map2cpu(x::AbstractString)=x
map2cpu(x::DataType)=x
map2cpu(x::Array)=map(map2cpu,x)
map2cpu{T<:Number}(x::Array{T})=x
map2cpu(x::Associative)=(y=Dict();for (k,v) in x; y[k] = map2cpu(x[k]); end; y)
map2cpu2(x)=(y=deepcopy(x); for f in fieldnames(x); setfield!(y,f,map2cpu(getfield(x,f))); end; y)

map2gpu(x)=(if isbits(x); x; else; map2gpu2(x); end)
map2gpu(x::KnetArray)=x
map2gpu(x::AbstractString)=x
map2gpu(x::DataType)=x
map2gpu(x::Tuple)=map(map2gpu,x)
map2gpu(x::Array)=map(map2gpu,x)
map2gpu{T<:AbstractFloat}(x::Array{T})=KnetArray(x)
map2gpu(x::Associative)=(y=Dict();for (k,v) in x; y[k] = map2gpu(x[k]); end; y)
map2gpu2(x)=(y=deepcopy(x); for f in fieldnames(x); setfield!(y,f,map2gpu(getfield(x,f))); end; y)

# Optimization parameters initialization
# initoptim creates optimization parameters for each numeric weight
# array in the model.  This should work for a model consisting of any
# combination of tuple/array/dict.
initoptim{T<:Number}(::KnetArray{T},otype)=eval(parse(otype))
initoptim{T<:Number}(::Array{T},otype)=eval(parse(otype))
initoptim(a::Associative,otype)=Dict(k=>initoptim(v,otype) for (k,v) in a)
initoptim(a,otype)=map(x->initoptim(x,otype), a)

macro msg(_x)
    :(join(STDOUT,[Dates.format(now(),"HH:MM:SS"), $(esc(_x)),'\n'],' '); flush(STDOUT))
end

date(x)=(join(STDOUT,[Dates.format(now(),"HH:MM:SS"), x,'\n'],' '); flush(STDOUT))

type StopWatch;
    tstart; nstart; ncurr; nnext;
    StopWatch()=new(time(),0,0,1000)
end

function inc(s::StopWatch, n, step=1000)
    s.ncurr += n
    if s.ncurr >= s.nnext
        tcurr = time()
        dt = tcurr - s.tstart
        dn = s.ncurr - s.nstart
        s.tstart = tcurr
        s.nstart = s.ncurr
        s.nnext += step
        return dn/dt
    end
end

# (:epoch,0,:beam,1,:las,0.007938020871520395,0.013241609670749164) for randseed 4
# (:epoch,0,:beam,1,:las,0.007938020871520395,0.013241609670749164)


text:
# Structs and functions related to the C file pointer type as used in eccodes

export GribFile, destroy

mutable struct File
end

"""
Represents a grib file. Implements many of the Julia I/O functions.
"""
mutable struct GribFile
    ptr::Ptr{File}
    filename::String
    mode::String
    pos::Int
    nmessages::Int32
end

"""
    `GribFile(filename::AbstractString, mode="r")`
Open a grib file. `mode` is a mode as described by `Base.open`.
"""
function GribFile(filename::AbstractString; mode="r")
    if !isfile(filename)
        throw(SystemError("opening file $filename: no such file"))
    end
    f = ccall(:fopen, Ptr{File}, (Cstring, Cstring), filename, mode)

    nref = Ref(Int32(0))
    err = ccall((:codes_count_in_file, eccodes), Cint, (Ptr{codes_context}, Ptr{File}, Ref{Cint}),
                C_NULL, f, nref)
    errorcheck(err)

    return GribFile(f, filename, mode, 0, nref[])
end

"""
    `GribFile(f::Function, filename::AbstractString, mode="r")`
Open a grib file and automatically close after exucuting `f`.
`mode` is a mode as described by `Base.open`.

# Example
```julia
GribFile(filename) do f
    # do things in read mode
end
```
"""
function GribFile(f::Function, filename::AbstractString; mode="r")
    file = GribFile(filename, mode=mode)
    try
        f(file)
    finally
        destroy(file)
    end
end

function Base.iterate(f::GribFile, state=())
    next = Message(f)
    if next == nothing
        return nothing
    else
        return (next, ())
    end
end

Base.eltype(f::GribFile) = Message

Base.IteratorSize(f::GribFile) = Base.SizeUnknown()

function Base.seekstart(f::GribFile)
    newptr = ccall(:fopen, Ptr{File}, (Cstring, Cstring), f.filename, f.mode)
    close(f.ptr)
    f.ptr = newptr
    f.pos = 0
end

""" Read without allocating a return vector """
function readnoreturn(f::GribFile, nm::Integer)
    total = f.nmessages
    nm = nm + f.pos > total ? total - f.pos : nm
    nm <= 0 && return
    for i in 1:nm
        Message(f)
    end
end

"""
    `read(f::GribFile[, nm::Integer])`

Read `nm` messages from `f` and return as vector. Default is 1.
"""
function Base.read(f::GribFile, nm::Integer)
    total = f.nmessages
    nm = nm + f.pos > total ? total - f.pos : nm
    nm <= 0 && return nothing
    ret = Vector{Message}(undef, nm)
    for i in 1:nm
        ret[i] = Message(f)
    end
    return ret
end
Base.read(f::GribFile) = Message(f)

"""
    `position(f::GribFile)`

Get the current position of the file.
"""
Base.position(f::GribFile) = f.pos

"""
    `seek(f::GribFile, n::Integer)`

Seek the file to the given position `n`
"""
function Base.seek(f::GribFile, n::Integer)
    if n < 0 || n > f.nmessages
        throw(DomainError("n is out of range for file length $(f.nmessages)"))
    end
    if n < f.pos
        seekstart(f)
        readnoreturn(f, n)
    elseif n > f.pos
        readnoreturn(f, n-f.pos)
    end
end

"""
    `skip(f::GribFile, offset::Integer)`

Seek the file relative to the current position
"""
function Base.skip(f::GribFile, offset::Integer)
    if f.pos + offset > f.nmessages
        throw(DomainError("offset is out of range for file length $(f.nmessages)"))
    end
    if offset < 0
        oldpos = f.pos
        seekstart(f)
        readnoreturn(f, oldpos+offset)
    else
        readnoreturn(f, offset)
    end
end

"""
    `destroy(f::GribFile)`
Safely close the file.
"""
function destroy(f::GribFile)
    err = ccall(:fclose, Cint, (Ptr{File},), f.ptr)
    if err !=0
        throw(ErrorException("GribFile closed with errorcode $err"))
    end
end

function close(f::Ptr{File})
    err = ccall(:fclose, Cint, (Ptr{File},), f)
    if err !=0
        throw(ErrorException("GribFile closed with errorcode $err"))
    end
end

# Functions related to printing
function Base.show(io::IO, mime::MIME"text/plain", f::GribFile)
    str = "GribFile $(f.filename) at position $(f.pos) in mode $(f.mode)"
    print(io, str)
end


text:
# I don't like the current implementation for timestepping!
"Explicit time advance."
function advance!(param, state_GV)
	(;verbose) = param

	faceState, faceGradient = init_face_value(param)

	faceFluxLR, faceFlux, speedFlux = init_flux(param)

	source_GV, U, div_G = init_source(param)

	time_G = init_timestep(param)

	t, it = 0.0, 0

	verbose && @info "Finished initialization..."

	if param.TimeAccurate # Advance with time
		if param.Order == 1
			while t < param.tEnd
				# Set boundary conditions
				@timeit timer() "BC" set_cell_boundary!(
					param, state_GV)

				# Calculate face value
				@timeit timer() "face value" calc_face_value!(
					param, state_GV, faceState, faceGradient)

				# Calculate face flux
				@timeit timer() "face flux" calc_face_flux!(
					param, faceState, faceFlux, speedFlux, faceFluxLR)

				# Calculate source
				@timeit timer() "source" calc_source!(
					param, state_GV, source_GV, U, div_G)

				# Calculate time step
				@timeit timer() "update dt" dt = calc_timestep!(
					param, speedFlux, time_G)

				if t + dt > param.tEnd dt = param.tEnd - t end

				# Update state
				@timeit timer() "update state" update_state!(
					param, state_GV, dt, faceFlux, source_GV)

				t  += dt
				it += 1

				if param.DoPlot && mod(it, param.PlotInterval) == 0
					@timeit timer() "plotvar" plotvar(param, it, state_GV)
				end

				verbose && @info "it, t = $it, $t"
			end
		elseif param.Order == 2 # 2nd order method
			state1_GV = similar(state_GV)

			while t < param.tEnd
				# 1st stage of modified timestepping

				# Set boundary conditions
				@timeit timer() "BC" set_cell_boundary!(
					param, state_GV)

				# Calculate face value
				@timeit timer() "face value" calc_face_value!(
					param, state_GV, faceState, faceGradient)

				# Calculate face flux
				@timeit timer() "face flux" calc_face_flux!(
					param, faceState, faceFlux, speedFlux, faceFluxLR)

				# Calculate source
				@timeit timer() "source" calc_source!(
					param, state_GV, source_GV, U, div_G)

				# Calculate time step
				@timeit timer() "update dt" dt = calc_timestep!(
					param, speedFlux, time_G)

				# Update state in the 1st stage
				dt *= 0.5
				state1_GV .= state_GV
				@timeit timer() "update state" update_state!(
					param, state1_GV, dt, faceFlux, source_GV)

				# 2nd stage of modified timestepping

				# Set boundary conditions
				@timeit timer() "BC" set_cell_boundary!(
					param, state1_GV)

				# Calculate face value
				@timeit timer() "face value" calc_face_value!(
					param, state1_GV, faceState, faceGradient)

				# Calculate face flux
				@timeit timer() "face flux" calc_face_flux!(
					param, faceState, faceFlux, speedFlux, faceFluxLR)

				# Calculate source
				@timeit timer() "source" calc_source!(
					param, state1_GV, source_GV, U, div_G)

				# Update state
				dt *= 2.0
				@timeit timer() "update state" update_state!(
					param, state_GV, dt, faceFlux, source_GV)

				if t + dt > param.tEnd dt = param.tEnd - t end

				t  += dt
				it += 1

				if param.DoPlot && mod(it, param.PlotInterval) == 0
					@timeit timer() "plotvar" plotvar(param, it, state_GV)
				end

				verbose && @info "it, t = $it, $t"
			end
		else
			@error "Higher Order schemes not yet implemented!"
		end
	else # Advance with steps
		state1_GV = similar(state_GV)

		for iStep = 1:param.nStep
			if param.Order == 1
				# Set boundary conditions
				@timeit timer() "BC" set_cell_boundary!(
					param, state_GV)

				# Calculate face value
				@timeit timer() "face value" calc_face_value!(
					param, state_GV, faceState, faceGradient)

				# Calculate face flux
				@timeit timer() "face flux" calc_face_flux!(
					param, faceState, faceFlux, speedFlux, faceFluxLR)

				# Calculate source
				@timeit timer() "source" calc_source!(
					param, state_GV, source_GV, U, div_G)

				# Calculate time step
				@timeit timer() "update dt" calc_timestep!(
					param, speedFlux, time_G)

				# Update state
				@timeit timer() "update state" update_state!(
					param,state_GV,time_G,faceFlux,source_GV)

				it += 1

				if param.DoPlot && mod(it, param.PlotInterval) == 0
					@timeit timer() "plotvar" plotvar(param, it, state_GV)
				end

				verbose && @info "it = $it"
			elseif param.Order == 2
				# 1st stage of modified timestepping

				# Set boundary conditions
				set_cell_boundary!(param, state_GV)

				# Calculate face value
				calc_face_value!(param, state_GV, faceState, faceGradient)

				# Calculate face flux
				calc_face_flux!(param, faceState, faceFlux, speedFlux, faceFluxLR)

				# Calculate source
				calc_source!(param, state_GV, source_GV, U, div_G)

				# Calculate time step
				calc_timestep!(param, speedFlux, time_G)

				# Update state in the 1st stage
				time_G .*= 0.5
				state1_GV .= state_GV
				update_state!(param, state1_GV, time_G, faceFlux, source_GV)

				# 2nd stage of modified timestepping

				# Set boundary conditions
				set_cell_boundary!(param, state1_GV)

				# Calculate face value
				calc_face_value!(param, state1_GV, faceState, faceGradient)

				# Calculate face flux
				calc_face_flux!(param, faceState, faceFlux, speedFlux, faceFluxLR)

				# Calculate source
				calc_source!(param, state1_GV, source_GV, U, div_G)

				# Update state
				time_G *= 2.0
				update_state!(param, state_GV, time_G, faceFlux, source_GV)

				it += 1

				if param.DoPlot && mod(it, param.PlotInterval) == 0
					plotvar(param, it, state_GV)
				end

				verbose && @info "it = $it"
			end
		end
	end
	return
end


text:
#For a Colab of this example, goto https://colab.research.google.com/drive/1xfUsBn9GEqbRjBF-UX_jnGjHZNtNsMae
using TensorBoardLogger
using Flux
using Logging
using MLDatasets
using Statistics
using CuArrays

#create tensorboard logger
logdir = "content/log"
logger = TBLogger(logdir, tb_overwrite)

#Load data
traindata, trainlabels = FashionMNIST.traindata();
testdata, testlabels = FashionMNIST.testdata();
trainsize = size(traindata, 3);
testsize = size(testdata, 3);

#Log some images
images = TBImage(traindata[:, :, 1:10], WHN)
with_logger(logger) do #log some samples
    @info "fmnist/samples" pics = images log_step_increment=0
end

#Create model
model = Chain(
    x -> reshape(x, :, size(x, 4)),
    Dense(28^2, 32, relu),
    Dense(32, 10),
    softmax
)

loss(x, y) = Flux.crossentropy(model(x), y)

accuracy(x, y) = mean(Flux.onecold(model(x) |> cpu) .== Flux.onecold(y |> cpu))

opt = ADAM()

traindata = permutedims(reshape(traindata, (28, 28, 60000, 1)), (1, 2, 4, 3));
testdata = permutedims(reshape(testdata, (28, 28, 10000, 1)), (1, 2, 4, 3));
trainlabels = Flux.onehotbatch(trainlabels, collect(0:9));
testlabels = Flux.onehotbatch(testlabels, collect(0:9));

#function to get dictionary of model parameters
function fill_param_dict!(dict, m, prefix)
    if m isa Chain
        for (i, layer) in enumerate(m.layers)
            fill_param_dict!(dict, layer, prefix*"layer_"*string(i)*"/"*string(layer)*"/")
        end
    else
        for fieldname in fieldnames(typeof(m))
            val = getfield(m, fieldname)
            if val isa AbstractArray
                val = vec(val)
            end
            dict[prefix*string(fieldname)] = val
        end
    end
end

#function to log information after every epoch
function TBCallback()
  param_dict = Dict{String, Any}()
  fill_param_dict!(param_dict, model, "")
  with_logger(logger) do
    @info "model" params=param_dict log_step_increment=0
    @info "train" loss=loss(traindata, trainlabels) acc=accuracy(traindata, trainlabels) log_step_increment=0
    @info "test" loss=loss(testdata, testlabels) acc=accuracy(testdata, testlabels)
  end
end

minibatches = []
batchsize = 100
for i in range(1, stop = trainsize÷batchsize)
  lbound = (i-1)*batchsize+1
  ubound = min(trainsize, i*batchsize)
  push!(minibatches, (traindata[:, :, :, lbound:ubound], trainlabels[:, lbound:ubound]))
end

Move data and model to gpu
traindata = traindata |> gpu
testdata = testdata |> gpu
trainlabels = trainlabels |> gpu
testlabels = testlabels |> gpu
model = model |> gpu
minibatches = minibatches |> gpu

#Train
@Flux.epochs 15 Flux.train!(loss, params(model), minibatches, opt, cb = Flux.throttle(TBCallback, 5))


text:
using IntroAlgoCh2
using Test
using Random, Plots

@testset "IntroAlgoCh2.jl" begin
    r = 1:1000 |> collect |> shuffle!
    @test issorted(insertionsort!(r)) == true
    @test issorted(mergesort!(r)) == true
end


text:
"""
grabstates(filename)

# Description
Constructs states for visualisation automatically from jld2 file

# Arguments
-  `filename`: a string that corresponds to a filename

# Return states, statenames
- `states`: array of arrays of data
- `statenames`: string names for states

"""
function grabstates(filename)
    file = jldopen(filename)
    ogrid = file["grid"]
    xS, yS, zS = size(ogrid)
    btmp = file["tracers"]["b"]["data"]
    ghost = Int((length(btmp[:,1,1])-xS)/2)
    b = file["tracers"]["b"]["data"][ghost:(xS + ghost - 1), ghost:(yS + ghost - 1), ghost:(zS + ghost - 1)]
    u = file["velocities"]["u"]["data"][ghost:(xS + ghost - 1), ghost:(yS + ghost - 1), ghost:(zS + ghost - 1)]
    v = file["velocities"]["v"]["data"][ghost:(xS + ghost - 1), ghost:(yS + 1 + ghost - 1), ghost:(zS + ghost - 1)]
    w = file["velocities"]["w"]["data"][ghost:(xS + ghost - 1), ghost:(yS + ghost - 1), ghost:(zS + 1 + ghost - 1)]
    x = ogrid.xC[1:xS]
    y = ogrid.yC[1:yS]
    z = ogrid.zC[1:zS]
    f = file["coriolis"].f₀
    β = file["coriolis"].β
    # Cell-Centered Values
    v = (v[:, 1:end-1, :] + v[:, 2:end, :]) .* 0.5
    w = (w[:, :, 1:end-1] + w[:,:, 2:end]) .* 0.5
    close(file)
    ## Create Gradient
    include(pwd() * "/analysis_scripts/" * "post_analysis.jl") # Gradients etc. here
    ∂x = PartialDerivative((x, 1))
    ∂y = PartialDerivative((y, 2))
    ∂z = PartialDerivative((z, 3))
    ∇ = [∂x, ∂y, ∂z]
    ω = ∇ × [u, v, w]
    ∇b = [∂x(b), ∂y(b), ∂z(b)]
    coriolis =  f .+ β * reshape(y, (1, length(y) ,1))
    fpβ = coriolis[:, 1:end-1, :] + coriolis[:, 2:end, :]  # f plus beta
    ω∇b = ω⋅∇b
    pv = ω∇b + fpβ .* ∇b[3]
    alignment = ω∇b ./ ( sqrt.(ω⋅ω) .* sqrt.(∇b⋅∇b) )
    # Hydrostatic Hydrostatic hydrostatic_pressure
    Δz = z[2] - z[1]
    p = cumsum(-b .* Δz, dims = 3) # since evenly spaced
    p = p .- mean(p)
    # Deformation Radius
    relu(x) = x ≥ 0 ? x : 0.0
    Lᵣ = -sqrt.(relu.(∇b[3])) ./ fpβ * abs(z[1]-z[end]) # (NH/f)
    # check Thermal Wind
    ∂xp = ∂x(p)
    fv = v .* coriolis
    ∂yp = ∂y(p)
    fu = u .* coriolis
    deviationx = ∂xp .+ (fv[1:end-1,1:end-1,1:end-1] + fv[2:end,2:end,2:end]) * 0.5 # off by a gridcell (should avg)
    deviationy = ∂yp .- (fu[1:end-1,1:end-1,1:end-1] + fu[2:end,2:end,2:end]) * 0.5
    velocity = [u, v, w]
    speed = sqrt.(velocity ⋅ velocity)
    states = [u, abs.(u), v, w, b, -p,
              speed, ω[3], ω[1], ω[2], abs.(ω∇b), abs.(pv),
              ∇b[3], Lᵣ, -∂xp, fv, -∂yp, -fu,
              deviationx, deviationy]
    statenames =  ["u", "|u|", "v", "w", "b", "-p/ρ₀ (hydrostatic)",
                  "speed", "ω₃", "ω₁", "ω₂", "|ω⋅∇b|",  "|Ertel PV|",
                  "∂z(b)", "Lᵣ", "-∂x(p/ρ₀)", "fv", "-∂y(p/ρ₀)", "-fu",
                  "deviationx", "deviationy"]
    units = [" [m/s]", " [m/s]", " [m/s]", " [m/s]", " [m/s²]", " [m²/s²]",
              " [m/s]", " [1/s]", " [1/s]", " [1/s]", " [1/s³]", " [1/s³]",
              " [m/s²]", " [m]", " [m/s²]", " [m/s²]", " [m/s²]", " [m/s²]",
              " [m/s²]", " [m/s²]"]
    return states, statenames, units
end

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
immutable Node{T}
    state::T
    score::Float64
    prev::Node{T}

    (::Type{Node}){T}(state::T) = new{T}(state, 0.0)
    (::Type{Node}){T}(state::T, score, prev) = new{T}(state, score, prev)
end

lessthan{T}(x::Node{T}, y::Node{T}) = x.score > y.score

"""
* next(state::T) -> T[], Int
* scorefun: (state::T) -> Float64
"""
function aaa{T}(beamsize::Int, initstate::T, target::Vector{Int})
    k = 1
    while k <= length(chart)
        beam = chart[k]
        length(beam) > beamsize && sort!(beam, lt=lessthan)
        for i = 1:beamsize
            prev = beam[i]
            nexts = next(prev)
            t = nexts[k]
            for s in nexts
                push!(chart[l], Node(st,0.0,prev))
            end
        end
    end
end

function beamsearch{T}(beamsize::Int, initstate::T, scorefun)
    chart = Vector{Node{T}}[]
    push!(chart, [Node(initstate,0.0)])
    k = 1
    while k <= length(chart)
        beam = chart[k]
        length(beam) > beamsize && sort!(beam, lt=lessthan)
        for i = 1:beamsize
            i > length(beam) && break
            prev = beam[i]
            for (st::T,l::Int) in next(prev,k)
                while l > length(chart)
                    push!(chart, Node{T}[])
                end
                push!(chart[l], Node(st,0.0,prev))
            end
        end
        k += 1
    end
    sort!(chart[end], lt=lessthan)
    chart
end

function seq{T}(node::Node{T})
    nodes = Node{T}[]
    n = node
    while true
        unshift!(nodes, n)
        isdefined(n, :prev) || break
        n = n.prev
    end
    unshift!(nodes, n)
    nodes
end

"""
* y: correct
* z: predicted
"""
function early(y::Node{T}, z::Node{T})

end

function maxviolation!{T}(y::Node{T}, z::Node{T})
    maxyz, maxv = (y,z), -1.0
    while true
        v = z.score - y.score
        if v > maxv
            maxv = v
            maxyz = (y,z)
        end
        isdefined(y, :prev) || break
        y = y.prev
        z = z.prev
    end
end

function max_violation!{T}(gold::T, pred::T, train_gold, train_pred)
    goldseq, predseq = to_seq(gold), to_seq(pred)
    maxk, maxv = 1, 0.0
    for k = 1:length(goldseq)
        v = predseq[k].score - goldseq[k].score
        if k == 1 || v >= maxv
            maxk = k
            maxv = v
        end
    end
    for k = 2:maxk
        train_gold(goldseq[k])
        train_pred(predseq[k])
    end
end


text:
# find and interact with NEWDATA in neo4j
using Caesar

using RoME, IncrementalInference
using CloudGraphs, Neo4j
# using Distributions
using JSON

## Uncomment out for command line operation
# cloudGraph, addrdict = standardcloudgraphsetup(nparticles=true)
# session = addrdict["session"]
# Nparticles = parse(Int,addrdict["num particles"])

# interactive operation
session = "SESSROX"
Nparticles = 100
include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))

# fieldnames(fg.cg.neo4j)

N=Nparticles
fg = initfg(sessionname=session, cloudgraph=cloudGraph)

# setBackendWorkingSet!(conn, session)
fullLocalGraphCopy!(fg)

# @show length(fg.IDs)
# showcurrentdlapi()

updatenewverts!(fg, N=N)

# getVal(fg, :x6)

@show length(fg.IDs)

writeGraphPdf(fg)
@async run(`evince fg.pdf`)

# Okay now solve locally, or start MM-iSAMCloudSolve.jl
fg = initfg(sessionname=session, cloudgraph=cloudGraph)

setDBAllReady!(fg)
setBackendWorkingSet!(conn, session)
fullLocalGraphCopy!(fg)

tree = wipeBuildNewTree!(fg,drawpdf=true)
@async run(`evince bt.pdf`)

# single core recursive inference, allowing better stacktrace
inferOverTreeR!(fg, tree, N=N)

# totally reset to frontend built state in DB
resetentireremotesession(fg.cg.neo4j.connection, session)

#


text:
using Combinatorics

# Part 1
data = parse.(Int64, readlines("day09_input.txt"))
valid = [d ∈ Set(sum.(combinations(data[i-25:i-1], 2))) for (i, d) in enumerate(data) if i > 25]
target = data[argmin(valid)+25]
println("Part 1: first invalid = $target")

# Part 2
function find_weakness(data, target)
    for start in 1:lastindex(data)
        total = data[start]
        for finish in (start+1):lastindex(data)
            total += data[finish]
            if total == target return minimum(data[start:finish]) + maximum(data[start:finish]) end
            if total > target break end
        end
    end
end
println("Part 2: weakness = $(find_weakness(data, target))")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode