text:
using Dates
using Statistics, DataFrames, Bootstrap
using OrderedCollections
import ..ME, ..MAE

mutable struct TemporalCrossSection{A<:Dates.Date, B<:Real, D<:OrderedDict}
    dates::Array{A}
    ts::Array{B}
    subperiods::D
    cross_section::Union{DataFrame, Nothing}
end

function TemporalCrossSection(dates::Array, ts::Array, period::Function=monthday)
    df = DataFrame(Date=dates, data=ts)
    sp = sort(unique(period.(dates)))

    # Remove leap days (when using monthday)
    deleteat!(sp, findall(x -> x == (2,29), sp))

    res = OrderedDict(
        sp_i => df[period.(df.Date) .== [sp_i], :].data for sp_i in sp
    )

    return TemporalCrossSection(dates, ts, res, nothing)
end

"""
TemporalCrossSection constructor that handles two time series
"""
function TemporalCrossSection(dates::Array, obs::Array, sim::Array,
                              period::Function=monthday)
    Y = ME.(obs, sim)
    return TemporalCrossSection(dates, Y, period)
end

function cross_section(tr::TemporalCrossSection)
    sp = tr.subperiods

    # Array holding boxplot-like data
    # No. of Dims: (75, 95) * (2) + median + min/max = 7
    n_cols = 7  # -q95, -q75, q50, q75, q95, min, max
    boxframe = fill(0.0, length(keys(sp)), n_cols)

    target_pts = [0.0, 0.025, 0.125, 0.875, 0.975, 1.0]
    for (i, (ki, vi)) in enumerate(sp)
        ab_min, lower_95, lower_75, upper_75, upper_95, ab_max = quantile(vi, target_pts)

        boxframe[i, 1:3] .= ab_min, lower_95, lower_75
        boxframe[i, 4] = median(vi)
        boxframe[i, 5:7] .= upper_75, upper_95, ab_max
    end

    cols = [:abs_min, :lower_95, :lower_75, :median,
            :upper_75, :upper_95, :abs_max]

    x_section = DataFrame(boxframe, cols)
    x_section[:, :subperiod] = collect(keys(sp))

    tr.cross_section = x_section

    return x_section
end

"""
    offsets(tr::TemporalCrossSection, period::Function=monthday)

Determine (median) offsets, assuming mean error-based temporal cross-sectional analysis
was performed. Only supports daily time series.
"""
function offsets(tr::TemporalCrossSection)
    x_section = tr.cross_section
    offset = x_section.median
    dates = tr.dates

    offset_vals = zeros(length(dates))

    doy = Dates.dayofyear.(dates)
    mds = monthday.(dates)
    for (idx, (d, md)) in enumerate(zip(doy, mds))
        if md == (2, 29) || d == 366
            offset_vals[idx] = offset[d-1] * -1.0
            continue
        end

        offset_vals[idx] = offset[d] * -1.0
    end

    return offset_vals
end

function Base.getproperty(a::TemporalCrossSection, v::Symbol)
    if v == :cross_section
        if isnothing(getfield(a, :cross_section))
            boxframe = cross_section(a)
            return boxframe
        end
    end

    if v == :offsets
        return offsets(a)
    end

    if in(v, [:cv_median, :cv_95, :cv_75, :mean_95, :mean_75, :std_75, :std_95])
        boxframe = a.cross_section

        if v == :cv_median
            med = boxframe[:, :median]
            cv_med = std(med) ./ mean(med)
            return cv_med
        elseif in(v, [:cv_95, :mean_95, :std_95])
            upper = boxframe[:, :upper_95]
            lower = boxframe[:, :lower_95]
            wr = upper .- lower

            if v == :cv_95
                return std(wr) ./ mean(wr)
            elseif v == :mean_95
                return mean(wr)
            elseif v == :std_95
                return std(wr, corrected=false)
            end
        elseif in(v, [:cv_75, :mean_75, :std_75])
            upper = boxframe[:, :upper_75]
            lower = boxframe[:, :lower_75]
            wr = upper .- lower

            if v == :cv_75
                return std(wr) ./ mean(wr)
            elseif v == :mean_75
                return mean(wr)
            elseif v == :std_75
                return std(wr, corrected=false)
            end
        end
    else
        return getfield(a, v)
    end
end


text:
using LayeredArrays
using LayeredArrays: LazyLayeredArray
using Test

struct MyType{T} <: AbstractLayeredVector{1, T}
    x::Vector{T}
end
Base.size(m::MyType) = size(m.x)
Base.getindex(m::MyType, i::Int) = getindex(m.x, i)
Base.setindex!(m::MyType, v, i::Int) = setindex!(m.x, v, i)

include("AbstractLayeredArray.jl")
include("LayeredArray.jl")
include("LazyLayeredArray.jl")

ext: jl


text:
XFloat32(x::XFloat16) = reinterpret(XFloat32, Float64(reinterpret(Float32,x)))
XFloat16(x::XFloat32) = reinterpret(XFloat16, Float32(reinterpret(Float64,x)))

XFloat32(x::Float64) = reinterpret(XFloat32, x)
XFloat32(x::Float32) = reinterpret(XFloat32, Float64(x))
XFloat32(x::Float16) = reinterpret(XFloat32, Float64(x))
XFloat32(x::Integer) = reinterpret(XFloat32, Float64(Float32(x)))
XFloat32(x::Real) = XFloat32(Float64(x))

XFloat16(x::Float64) = reinterpret(XFloat16, Float32(x))
XFloat16(x::Float32) = reinterpret(XFloat16, x)
XFloat16(x::Float16) = reinterpret(XFloat16, Float32(x))
XFloat16(x::Integer) = reinterpret(XFloat16, Float32(Float16(x)))
XFloat16(x::Real) = XFloat16(Float32(x))

Float64(x::XFloat32) = reinterpret(Float64, x)
Float32(x::XFloat32) = Float32(Float64(x))
Float16(x::XFloat32) = Float16(Float64(x))

Float64(x::XFloat16) = Float64(Float32(x))
Float32(x::XFloat16) = reinterpret(Float32, x)
Float16(x::XFloat16) = Float16(Float32(x))

for I in (:Int8, :Int16, :Int32, :Int64, :Int128, :BigInt,
          :UInt8, :UInt16, :UInt32, :UInt64, :UInt128,)
    @eval begin
        $I(x::XFloat32) = $I(Float64(x))
        $I(x::XFloat16) = $I(Float32(x))
    end
end

XFloat16(x::BigFloat) = reinterpret(XFloat16, Float32(x))
XFloat32(x::BigFloat) = reinterpret(XFloat32, Float64(x))
BigFloat(x::XFloat16) = BigFloat(reinterpret(Float32, x))
BigFloat(x::XFloat32) = BigFloat(reinterpret(Float64, x))

XFloat16(x::Bool) = reinterpret(XFloat16, Float32(x))
XFloat32(x::Bool) = reinterpret(XFloat32, Float64(x))
Bool(x::XFloat16) = Bool(reinterpret(Float32, x))
Bool(x::XFloat32) = Bool(reinterpret(Float64, x))

XFloat16(x::Irrational{S}) where {S} = XFloat16(Float32(x))
XFloat32(x::Irrational{S}) where {S} = XFloat32(Float64(x))

function Base.frexp(x::XFloat16)
  fr,xp = frexp(Float32(x))
  return XFloat16(fr), xp
end

function Base.frexp(x::XFloat32)
  fr,xp = frexp(Float64(x))
  return XFloat32(fr), xp
end

function Base.ldexp(x::XFloat16, xp::I) where {I<:Integer}
   fr = Float32(x)
   ldxp = ldexp(fr, xp)
   return XFloat16(ldxp)
end

function Base.ldexp(x::XFloat32, xp::I) where {I<:Integer}
   fr = Float64(x)
   ldxp = ldexp(fr, xp)
   return XFloat32(ldxp)
end


text:
"""
     l63!(du,u,p,t)

 Steps forward the Lorenz63 equations for time `t`.
 Prognostic variables x, y, z are defined as
 elements 1, 2 and 3 of vector `u` and with parameters
 σ, ρ and β defined as elements of vector `p`.
 Elements of du are set to dxdt, dydt and dzdt respectively.

```
 julia> using DifferentialEquations, Plots
 julia> u0=[1.,0.,0.]
 julia> p=[10.,28.,8/3]
 julia> tspan=(0.,100.)
 julia> prob=ODEProblem(l63!,u0,tspan,p)
 julia> sol=solve(prob)
 julia> plot(sol,vars=(1,2,3))
```

"""
function l63!(du,u,p,t)
        x, y, z=u
        σ, ρ, β=p

        dxdt=σ*(y-x)
        dydt=ρ*x-x*z-y
        dzdt=x*y-β*z

        du[1], du[2], du[3]=
        dxdt , dydt,  dzdt
        return
end


text:
function loadpng(fname::String)
    png_img = load(fname)
    img_channel_view = channelview(png_img)
    img_re = reinterpret(UInt16, img_channel_view[1,:,:])
    return copy(img_re)
end


text:
using TensorCast
using Test
using LinearAlgebra
using StaticArrays
using OffsetArrays
using Einsum
using Strided
using Compat
if VERSION >= v"1.1"
    using LoopVectorization
end

@testset "ex-@shape" begin include("shape.jl") end
@testset "@reduce" begin include("reduce.jl")  end
@testset "@cast"   begin include("casting.jl") end
@testset "@matmul" begin include("mul.jl")     end

@testset "slice/view"  begin include("cat.jl") end
@testset "old readmes" begin include("old.jl") end
@testset "new in 0.2"  begin include("two.jl") end


text:
using BinDeps

@BinDeps.setup

sqlite3_lib = library_dependency("sqlite3_lib",
                                 aliases = ["sqlite3", "sqlite3-64",
                                            "libsqlite3", "libsqlite3-0"])
# @osx_only begin
#     using Homebrew
#     provides(Homebrew.HB, "sqlite", sqlite3_lib, os = :Darwin)
# end

@BinDeps.install [:sqlite3_lib => :sqlite3_lib ]


text:
module TIVM

using GenericInstruments
using GenericInstruments: viWrite, query, viRead,viRead!, readavailable
using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using ImPlot
using Statistics
using FFTW
import CImGui.LibCImGui: ImGuiCond_Always, ImGuiCond_Once

const GI = GenericInstruments

# VISA
include("visa_utils.jl")

# TYPES
include("types.jl")

# PST3201
include("PST3201_driver.jl")
include("PST3201_utils.jl")

# GDM8246
include("GDM8246_driver.jl")
include("GDM8246_utils.jl")

# GFG3015
include("GFG3015_driver.jl")
include("GFG3015_utils.jl")

# DS1000Z scope
include("DS1000Z_driver.jl")
include("DS1000Z_utils.jl")

# DS1000Z fgen
include("DS1000Z_FGEN_driver.jl")
include("DS1000Z_FGEN_utils.jl")

# TDS2002B
include("TDS2002B_driver.jl")
include("TDS2002B_utils.jl")

# Relays
include("Relays_driver.jl")

# GUI
include("gui_utils.jl")
include("psu_gui.jl")
include("dmm_gui.jl")
include("fgen_gui.jl")
include("scope_gui.jl")
include("relays_gui.jl")
include("main_gui.jl")

# ADALM2000
#const jllibm2k = joinpath(@__DIR__, "jllibm2k.dll")
#include("ADALM2000/jllibm2k_api.jl")
#include("ADALM2000/analog/m2kanalogin.jl")
#include("ADALM2000/analog/m2kanalogout.jl")

# export utils
export connect!, disconnect!, write, read, viWrite, viRead, query, start_gui, visaRead, visaWrite, find_resources

# export dmm functions
export get_sense_func, set_sense_func, set_sense_range, get_sense_range, set_sense_range_auto, get_sense_range_auto,
	get_primary_measurement, get_secondary_measurement

# export psu functions
export set_source_lev, get_source_lev, set_volt_protection, get_volt_protection, set_curr_protection, get_curr_protection,
	set_max_curr, get_max_curr, reset_protections, get_meas, set_outp, get_outp

# export fgen functions
export 	set_wfm, get_wfm, set_amplit_unit, get_amplit_unit, set_amplit, get_amplit, set_offs, get_offs, set_freq,
	get_freq, set_duty, get_duty

# export scope functions
export get_vertical_scale, get_horizontal_scale, get_ch_position, set_ch_position, set_vertical_offset, get_trig_data, get_meas_data,
	set_vertical_scale, set_horizontal_scale, set_trig_ch, set_trig_level, set_trig_mode, set_meas_ch, set_meas_type,
	set_meas, conf_acq_ch, Trigger_Aquistion

# export relays functions
export set_state, get_state, get_idn, all_off

end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
X = collect(Float64, reshape(1:40, 10, 4))
x = rand(10) * 10

D = DataFrame(A=rand(10), B=1:10, C=string.(11:20; base=16))
D_NA = allowmissing!(copy(D))
D_NA[1, :A] = missing

@testset "Array" begin
    XX = deepcopy(X)
    mu = center!(XX, obsdim=1)
    @test sum(abs.(mean(XX; dims=1))) == 0
    @test all(std(XX; dims=1) .== std(X; dims=1))
    @test all(mu .== vec(mean(X; dims=1)))

    XX = deepcopy(X)
    mu = center!(XX, ObsDim.First())
    @test sum(abs.(mean(XX; dims=1))) == 0
    @test all(std(XX; dims=1) .== std(X; dims=1))
    @test all(mu .== vec(mean(X; dims=1)))

    XX = deepcopy(X)
    mu = center!(XX, ObsDim.Last())
    @test sum(abs.(mean(XX; dims=2))) == 0
    @test all(std(XX; dims=2) .== std(X; dims=2))
    @test all(mu .== vec(mean(X; dims=2)))

    XX = deepcopy(X)
    mu = center!(XX)
    @test sum(abs.(mean(XX; dims=2))) == 0
    @test all(std(XX; dims=2) .== std(X; dims=2))
    @test all(mu .== vec(mean(X; dims=2)))

    XX = deepcopy(X)
    mu = vec(mean(X; dims=1))
    center!(XX, mu, obsdim=1)
    @test sum(abs.(mean(XX; dims=1))) == 0
    @test all(std(XX; dims=1) .== std(X; dims=1))

    XX = deepcopy(X)
    mu = vec(mean(X; dims=1))
    center!(XX, mu, ObsDim.First())
    @test sum(abs.(mean(XX; dims=1))) == 0
    @test all(std(XX; dims=1) .== std(X; dims=1))

    XX = deepcopy(X)
    mu = vec(mean(XX; dims=2))
    center!(XX, mu, obsdim=2)
    @test sum(abs.(mean(XX; dims=2))) == 0
    @test all(std(XX; dims=2) .== std(X; dims=2))

    XX = deepcopy(X)
    mu = vec(mean(XX; dims=2))
    center!(XX, mu, ObsDim.Last())
    @test sum(abs.(mean(XX; dims=2))) == 0
    @test all(std(XX; dims=2) .== std(X; dims=2))

    XX = deepcopy(X)
    mu = vec(mean(X[:,[1,3]]; dims=1))
    center!(XX, mu, obsdim=1, operate_on=[1, 3])
    @test sum(abs.(mean(XX[:,[1,3]]; dims=1))) == 0
    @test all(XX[:,2] .== X[:,2])
    @test all(std(XX; dims=1) .== std(X; dims=1))

    XX = deepcopy(X)
    mu = vec(mean(X[[1,3],:]; dims=2))
    center!(XX, mu, obsdim=2, operate_on=[1, 3])
    @test sum(abs.(mean(XX[[1,3],:]; dims=2))) == 0
    @test all(XX[2,:] .== X[2,:])
    @test all(std(XX; dims=2) .== std(X; dims=2))
    println()

    xx = deepcopy(x)
    center!(xx)
    @test mean(xx) <= 10e-10

    xx = deepcopy(x)
    mu = mean(xx)
    center!(xx, mu)
    @test mean(xx) <= 10e-10

    xx = deepcopy(x)
    mu = ones(size(xx))
    center!(xx, mu)
    @test mean(xx) - mean(x) ≈ -1

    xx = deepcopy(x)
    mu = ones(size(xx))
    center!(xx, mu)
    @test mean(xx) - mean(x) ≈ -1
end

@testset "DataFrame" begin
    # Center DataFrame
    DD = deepcopy(D)
    center!(DD)
    @test abs.(mean(DD[:A])) <= 10e-10
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D)
    center!(DD, operate_on=[:B])
    @test all(DD[:A] .== D[:A])
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D)
    mu = center!(DD, operate_on=[:A, :B])
    @test abs.(mean(DD[:A])) <= 10e-10
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[:C] .== D[:C])
    @test all(mu .== [mean(D[:A]), mean(D[:B])])

    DD = deepcopy(D_NA)
    mu = center!(DD, [:A, :B, :C])
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[2:end, :A] .== D[2:end, :A])
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D_NA)
    mu = [0.0, mean(DD[:B]), 0.0]
    mu = center!(DD, mu, [:A, :B, :C])
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[2:end, :A] .== D[2:end, :A])
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D_NA)
    mu = 0.0
    mu = center!(DD, mu, :A)
    @test all(DD[2:end, :A] .== D[2:end, :A])
    @test all(DD[:B] .== D[:B])
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D)
    mu =  [mean(D[:A]), mean(D[:B])]
    @test all(center!(DD, mu, operate_on=[:A, :B]) .== mu)
    @test abs.(mean(DD[:A])) <= 10e-10
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D_NA)
    center!(DD)
    @test all(DD[2:end, :A] .== D[2:end, :A])
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[:C] .== D[:C])
    @test ismissing(DD[1, :A])
end


text:
using DataStructures
using Serialization

const iobs = Stack(IOBuffer)

function get_iob()
    try
        return pop!(iobs)
    catch e
        return IOBuffer()
    end
end

function unget_iob(iob)
    push!(iobs, iob)
end

"""
    byte_array(x)

Serialize Julia type *x* such that it can be written to
a db or file
"""
function byte_array(x)
    iob = IOBuffer()
    serialize(iob, x)
    iob.data
end

struct KVSSerializeException <: Exception
    msg::String
end

"""
    array_to_type{T}(arr, ::Type{T})

Deserialize from Array{UInt8} back to Julia type.
"""
function array_to_type(arr, ::Type{T}) where T
    iob = IOBuffer(arr)
    seek(iob, 0)
    t = deserialize(iob)
    if ! isa(t, T)
        throw(KVSSerializeException("Could not deserialize to type"))
    end
    t
end

function array_to_type(arr)
    iob = IOBuffer(arr)
    seek(iob, 0)
    t = deserialize(iob)
    t
end


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
using GenerativeAD.Models
using BSON
using StatsBase
using StatsBase: fit!, predict, sample
using MLDataPattern
using Distances
using Statistics

s = ArgParseSettings()
@add_arg_table! s begin
	"max_seed"
		default = 1
		arg_type = Int
		help = "max_seed"
	"category"
		default = "wood"
		arg_type = String
		help = "category"
	"contamination"
		arg_type = Float64
		help = "contamination rate of training data"
		default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack category, max_seed, contamination = parsed_args

modelname = "ocsvm"

function median_l2_dist(X)
	dists = pairwise(Euclidean(), X)
	# take only the upper diagonal
	ds = []
	for i in 1:size(dists,1)
		for j in (i+1):size(dists,1)
			push!(ds, dists[i,j])
		end
	end
	ml=median(ds)
end
function set_params(data)
	D, N = size(data[1][1])
	gamma = 1/median_l2_dist(data[1][1])
	return (gamma=gamma, kernel="rbf", nu=0.5)
end

function joint_fit(models, data_splits)
	info = []
	for (model, data) in zip(models, data_splits)
		push!(info, fit!(model, data))
	end
	return info
end

function joint_prediction(models, data)
	joint_pred = Array{Float32}(undef, length(models), size(data,2))
	for (i,model) in enumerate(models)
		joint_pred[i,:] = predict(model, data)
	end
	return vec(mean(joint_pred, dims=1))
end

function StatsBase.fit(data, parameters, n_models=1)
	# construct model - constructor should only accept kwargs
	models = [GenerativeAD.Models.OCSVM(;parameters...) for _ = 1:n_models]

	# sumbsample and fit train data
	tr_data = data[1][1]
	M,N = size(tr_data)
	split_size = floor(N/n_models)
	data_splits = [tr_data[:,Int(split_size*i+1):Int(split_size*(i+1))] for i = 0:n_models-1]

	try
		global info, fit_t, _, _, _ = @timed joint_fit(models, data_splits)
 	catch e
		# return an empty array if fit fails so nothing is computed
		return (fit_t = NaN,), []
	end

	# construct return information - put e.g. the model structure here for generative models
	training_info = (
		fit_t = fit_t,
		model = nothing
		)

	# now return the different scoring functions
	training_info, [(x->joint_prediction(models, x), parameters)]
end

try_counter = 0
max_tries = 10*max_seed
cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
while try_counter < max_tries
	for seed in 1:max_seed
		savepath = datadir("experiments/images_mvtec_clean_val_default$cont_string/$(modelname)/$(category)/ac=1/seed=$(seed)")
		mkpath(savepath)

		# get data
		data = GenerativeAD.load_data("MVTec-AD", seed=seed, category=category,
			contamination=contamination, img_size=128)
		data = GenerativeAD.Datasets.vectorize(data)

		# parameters
		parameters = set_params(data)

		# here, check if a model with the same parameters was already tested
		if GenerativeAD.check_params(savepath, parameters)
			training_info, results = fit(data, parameters)
			# here define what additional info should be saved together with parameters, scores, labels and predict times
			save_entries = merge(training_info, (modelname = modelname, seed = seed,
				category = category, dataset = "MVTec-AD_$category",
				contamination=contamination))
			# now loop over all anomaly score funs
			for result in results
				GenerativeAD.experiment(result..., data, savepath; save_entries...)
			end
			global try_counter = max_tries + 1
		else
			@info "Model already present, sampling new hyperparameters..."
			global try_counter += 1
		end
	end
end
(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing


text:
using Pkg
Pkg.activate(".")
using PlotlyJS
using RDatasets
using Dash, DashHtmlComponents, DashCoreComponents, DashBootstrapComponents

struct Style
    layout::Layout
    global_trace::PlotlyBase.PlotlyAttribute
end

iris = RDatasets.dataset("datasets", "iris")
#
# Creation of Dash App starts here
#

app = dash(external_stylesheets=[dbc_themes.CYBORG])
app.title = "Iris Dataset"

app.layout = html_div(
    [html_center([html_h1("Iris Dataset Scatter Plot")]),html_br(),
    html_center([dcc_graph(id = "graph"),html_div(style=Dict(
        "max-width" => "620px"),[html_br(),html_p("Petal Width Range"),dcc_rangeslider(
            id = "slider",
            min = 0.0,
            max = 2.5,
            step=0.1,
            marks = Dict(0=>"0.0",2.5=>"2.5"),
            value = [0.0,2.5]
        )])]),
    html_br(),

])
callback!(app, Output("graph", "figure"), Input("slider", "value")) do input_value
    low = input_value[1]
    high = input_value[2]
    layout = Layout(;title="",width=600, height=400,xaxis=attr(title="Sepal Width"), yaxis=attr(title="Sepal Length"),template=templates["plotly_dark"])
    mask = (low.<iris."PetalWidth") .& (iris."PetalWidth" .< high)
    scale = 5.0

    dc = Dict("setosa"=>1, "virginica"=>2, "versicolor"=>3)
    cat = [dc[string(f)] for f = iris."Species"]
    sw = ["Petal Width = $w" for w in iris."PetalWidth"]
    fig = scatter(iris[mask,:], x=:SepalWidth, y=:SepalLength, marker_color=cat, mode="markers", marker_size=scale*iris."PetalLength", text=sw)
    plot(fig, layout)
end

run_server(app, "0.0.0.0", debug=true)


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec julia --project=@. --color=yes --startup-file=no -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" "$@"
=#

import CPLEX, Gurobi, GLPK
using JuMP

function test_duals(solver_sym, relaxed_mip)
	solver_lib = getfield(Main, solver_sym)
	m = direct_model(solver_lib.Optimizer())
	set_silent(m)
	if relaxed_mip
		@variable(m, x[1:10] >= 0, Int)
	else
		@variable(m, x[1:10] >= 0)
	end
	@objective(m, Min, sum(x))
	relaxed_mip && unset_integer.(x)
	optimize!(m)
	@show solver_sym
	@show relaxed_mip
	@show has_duals(m)
end

test_duals(:GLPK, false)
test_duals(:GLPK, true)
test_duals(:CPLEX, false)
test_duals(:CPLEX, true)
test_duals(:Gurobi, false)
test_duals(:Gurobi, true)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "LHAPDF"
version = v"6.3.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://lhapdf.hepforge.org/downloads/?f=LHAPDF-$(version).tar.gz",
                  "ed4d8772b7e6be26d1a7682a13c87338d67821847aa1640d78d67d2cef8b9b5d"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/LHAPDF*/
atomic_patch -p1 ../patches/undefined-setenv-mingw.patch
if [[ "${target}" == *-mingw* ]]; then
    # The macro `M_PI` is defined conditionally inside `math.h`.  Use
    # `_GNU_SOURCE` to enable it.
    export CPPFLAGS="-D_GNU_SOURCE"
    # Pass `-no-undefined` to the linker
    atomic_patch -p1 ../patches/ldflags-no-undefined.patch
    autoreconf -fiv
fi
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-python
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

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

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Setup.jl: define structs Sets, Game, and Population
# __precompile__()
# module Setup

export Sets, Game, Population
export compute_total_memberships
export assign_parties
export set_members_and_pairs
export random_sets

# load packages
using Base
using Random
using StatsBase
using Statistics
using Combinatorics

mutable struct Sets
    # structure for storing information about set memberships, opinions, and parties

    N::Int64   # number of individuals
    M::Int64   # number of sets
    K::Int64   # number of sets to which each individual belongs
    P::Int64   # number of parties to which each individual belongs
    α::Float64 # fraction of individuals who are independent
    set_members::Array{Array{Int64, 1}, 1}            # lists of individuals in each set
    set_pairs::Array{Array{Tuple{Int64, Int64}, 1},1} # list of tuples, each tuple is a pair of individuals in a set
    h::Array{Int64, 2}      # N-by-M int array denoting *opinions*
                            # call h[:, k] to get all opinions in set k
                            # or h[i, :] to get all of i's opinions
    h_bar::Array{Int64,2}   # N-by-M int array denoting *set memberships*
                            # call h_bar[:, k] to get all members in set k
                            # or h_bar[i, :] to get all of i's set memberships
    h_base10::Array{Int64, 1}      # N-by-1 vector denoting *opinions* in base 10 (from -(3^M-1)/2 to (3^M-1)/2)
    h_bar_base10::Array{Int64, 1}  # N-by-1 vector denoting *set memberships* in base 10 (from 1 to (3^M-1)/2)
    affils::Array{Int64, 1} # N-by-1 int array denoting *party affils* (0, ..., P-1)
                            # call p[i] to get part membership of i
    n_pairs::Int64          # number of pairs = binomial(N,2) for computing pairwise distances later
    n_party_pairs::Int64    # number of within-party pairs = binomial(N÷2, 2)

    # constructors for pre-specified h
    function Sets( h::Array{Int64, 2}, P::Int64=2, α::Float64=0.)

        N, M   = size(h) # extract numbers of individuals and sets
        h_bar  = [abs(x) for x in h]              # extract set memberships from h
        K      = compute_total_memberships(h_bar) # each member must belong to K sets
        affils = assign_parties(N,P,α)            # each member is affiliated with a party
        n_pairs, n_party_pairs = binomial(N, 2), binomial(N÷2, 2) # number of pairs
        set_members, set_pairs = set_members_and_pairs(h_bar)

        h_base10     = [vectobase3(h[i,:]) for i in 1:N] # track opinions in base10
        h_bar_base10 = [vectobase3(h_bar[i,:]) for i in 1:N] # track set memberships in base10

        return new( N, M, K, P, α, set_members, set_pairs, h, h_bar, h_base10, h_bar_base10, affils, n_pairs, n_party_pairs)
    end

    # alternative constructor with affils
    function Sets( h::Array{Int64, 2}, affils::Array{Int64, 1}, P::Int64=2, α::Float64=0.)

        N, M  = size(h) # extract numbers of individuals and sets
        h_bar = [abs(x) for x in h]              # extract set memberships from h
        K     = compute_total_memberships(h_bar) # each member must belong to K sets
        n_pairs, n_party_pairs = binomial(N, 2), binomial(N÷2, 2) # number of pairs
        set_members, set_pairs = set_members_and_pairs(h_bar)

        h_base10     = [vectobase3(h[i,:]) for i in 1:N] # track opinions in base10
        h_bar_base10 = [vectobase3(h_bar[i,:]) for i in 1:N] # track set memberships in base10

        return new( N, M, K, P, α, set_members, set_pairs, h, h_bar, h_base10, h_bar_base10, affils, n_pairs, n_party_pairs)
    end
end

function compute_total_memberships( h_bar::Array{Int64, 2} )::Int64
    # function to compute K, the number of sets to which each individual belongs
    total_memberships = sum( h_bar, dims=2 ) # total no. of sets indivduals belong to

    # check that all individuals belong to the same number of sets (K)
    length( unique(total_memberships) ) == 1 ? K = total_memberships[1] : @warn("check that each individual belongs to K sets")

    return K
end

function set_members_and_pairs( h::Array{Int64, 2} )
    # return a list of list of set members given a boolean array of memberships

    N, M = size(h)
    set_members = [Int64[] for k in 1:M]
    for indices in findall(x->x in [1,-1], h)
        push!(set_members[indices[2]],indices[1])
    end

    set_pairs = [Tuple{Int64, Int64}[] for k in 1:M]
    for k in 1:M
        for i in set_members[k]
            for j in set_members[k]
                if i < j push!(set_pairs[k], (i, j)) end
            end
        end
    end
    return set_members, set_pairs
end

function assign_parties( N::Int64, P::Int64=2, α::Float64=0. )
    # return party assignments p given N individuals and P parties

    # P = 2 only for now; warn if the number of parties is greater than 2
    if P > 2 @warn("Party assignment for P>2 is undefined"); return; end
    # warn if N*α is not an integer
    if N * α != round(N*α) @warn("Number of independents is not an integer, check value of α"); return; end
    # warn if the number of party members is indivisible by P
    if N * (1-α) % P != 0 @warn("Parties are uneven in size!"); return; end

    if P == 1 # if everyone is either independent (0) or in party 1
        affils = vcat( zeros(Int64, Int(N*α)), ones(Int64, Int(N*(1-α))) )
    elseif P == 2 # without loss of generality, first N*α are independents, then N*(1-α)/2 each for party 1 and party 2
        affils = vcat( zeros(Int64, Int(N*α)), ones(Int64, Int(N*(1-α)/2)), repeat([2], Int(N*(1-α)/2)) )
    end

    return affils
end

function random_sets( N::Int64, M::Int64, K::Int64, P::Int64=2, α::Float64=0. )
    # cosntructor for randomized set membership
    # where each individual belongs to K > 0 sets (no loners)
    h = zeros(Int64, N, M) # initializing
    affils = assign_parties( N, P, α ) # determine affils

    if P == 2
        # updated: use affils explicitly in case party assignments are changed
        for i in 1:N
            if affils[i] == 0  h[i,:] = shuffle( vcat( zeros(Int64,M-K), rand([-1,1],K) ) ) end # independents
            if affils[i] == 1  h[i,:] = shuffle( vcat( zeros(Int64,M-K), -ones(Int64,K) ) ) end # party 1
            if affils[i] == 2  h[i,:] = shuffle( vcat( zeros(Int64,M-K),  ones(Int64,K) ) ) end # party 2
        end
    elseif P == 1
        # for each individual, pick K random sets to belong to
        # for each of the K sets, choose -1 or +1 opinion randomly
        for i in 1:N h[i,:] = shuffle( vcat( zeros(Int64,M-K), rand([-1,1],K) ) ) end
    end

    return Sets(h, affils)
end

function vectobase3( vec::Array{Int64, 1} )::Int64
    # convert vector of opinions to int by interpreting the vector as digits in base 3
    foldl( (x::Int64, y::Int64) -> 3 * x + y, vec )
end

struct Game
    # structure for storing game parameters

    b::Float64    # benefit to cooperating
    c::Float64    # cost to cooperating
    β::Float64    # selection strength
    u::Float64    # strategy mutation rate
    v::Float64    # set mutation rate
    # p::Float64    # party bias (p = 0 no bias, p = 1 strong bias)
    ps::Array{Float64, 1}   # array of party bias p (p = 0 no bias, p = 1 strong bias)
                            # if ps = [p], everyone has the same P
                            # if ps = [p1, p2], party 1 has p1, party 2 has p2
    ϵ::Float64    # damping parameter attenuating party bias (ε = 1 no damping, ε = 0 full damping)
    A::Array{Float64, 2} # the game matrix

    # function Game( b::Float64, c::Float64, β::Float64, u::Float64, v::Float64, p::Float64, ε::Float64 )
    #     return new(b, c, β, u, v, p, ε, [0.0 b; -c b-c])
    # end
    # constructor for uniform p: if p is a Float, assume everyone has the same P
    function Game( b::Float64, c::Float64, β::Float64, u::Float64, v::Float64, p::Float64, ϵ::Float64 )
        return new(b, c, β, u, v, [p, p], ϵ, [0.0 b; -c b-c])
    end
    # constructor for two values of p: if p is an Array with two elements, assume parties have different p's
    function Game( b::Float64, c::Float64, β::Float64, u::Float64, v::Float64, ps::Array{Float64, 1}, ϵ::Float64 )
        if length(ps) != 2 @error("ps should have exactly two entries"); return end
        return new(b, c, β, u, v, ps, ϵ, [0.0 b; -c b-c])
    end
end

mutable struct Simulation
    # structure storing information about the simulation
    # this is temporary storage *within* a generation
    # any info to be carried over for tracking should be in pop

    i_action::Int64                   # temporarily store action of i in a given pairwise game
    j_action::Int64                   # temporarily store action of j in a given pairwise game
    learner::Int64                    # temporarily store learner in a given generation
    non_learners::Array{Int64,1}      # temporarily store list of individuals other than the learner
    role_model::Int64                 # temporarily store role model in a given generation
    strategy_copy::Array{Int64,1}     # temporarily store learner's strategies before imitation
    sets_copy::Array{Int64,1}         # temporarily store learner's opinions before imitation
    fitnesses::Array{Float64,1}       # temporarily store fitnesses of everyone but the learner
    cityblock_ind::Array{Int64,2}     # temporarily store pairwise cityblock distance between individuals
    cityblock_party::Array{Int64,2}   # temporarily store within-party pairwise cityblock distance between individuals
    hamming_ind::Array{Int64,2}       # temporarily store pairwise hamming distance between individuals
    hamming_party::Array{Int64,2}     # temporarily store within-party hamming cityblock distance between individuals
    cooperation::Array{Int64, 3}      # temporarily store within-party or between-party cooperation

    # constructor
    function Simulation( sets::Sets )
        # initialize
        i_action, j_action, learner, role_model = -1, -1, -1, -1  # initializing with a value other than 0(D)/1(C)
        non_learners    = zeros(Int64, sets.N-1)
        strategy_copy   = [-1, -1]
        sets_copy       = -ones(Int64, sets.M)
        fitnesses       = -ones(Float64, sets.N) # OLD: -ones(Float64, sets.N-1)

        cityblock_ind   = zeros(Int64, sets.N, sets.N)
        cityblock_party = zeros(Int64, sets.N÷2, sets.N÷2)
        hamming_ind     = zeros(Int64, sets.N, sets.N)
        hamming_party   = zeros(Int64, sets.N÷2, sets.N÷2)
        cooperation     = zeros(Int64, sets.M, sets.N÷2, sets.N÷2)

        return new(i_action, j_action, learner, non_learners, role_model, strategy_copy, sets_copy, fitnesses,
                   cityblock_ind, cityblock_party, hamming_ind, hamming_party, cooperation)
    end
end

mutable struct Population
    # structure storing information about the population

    sets::Sets                        # the set distribution
    game::Game                        # not mutable: game attributes and parameters
    sim::Simulation                   # temporary info in simulations
    strategies::Array{Int64, 2}       # array of game strategies - 0:DD, 1:DC, 2:CD, 3:CC
    payoffs::Array{Float64, 1}        # array of payoffs
    all_strategies::Array{Int64, 1}   # which strategies are allowed to appear
    generation::Int64                 # current generation
    verbose::Bool                     # turn this on for error tracking
    sets_changed::Bool                # track changes in sets
    strategies_changed::Bool          # track changes in strategies
    prev_actions::Array{Int64, 3}     # the last action each individual took toward each other
    prev_learner::Int64               # track learner in the previous round
    payoffs_mat::Array{Float64, 2}    # track pairwise payoffs
    all_opinions::Array{Int64, 1}     # list of all possible opinions (in base10)
    all_memberships::Array{Int64,1}   # list of all possible set memberships (in base10)
    strategies_base10::Array{Int64,1} # list of strategies in base 10
    prev_interactions::Int64          # track number of interactions
    num_imitating::Int64              # track number of imitation events
    num_not_imitating::Int64          # track number of !(imitation events)
    num_strategies_changed::Int64     # track number of strategy changes
    num_sets_changed::Int64           # track number of opinion changes
    num_both_changed::Int64           # track number of times when both strategy + opinions change
    num_random_opinions::Int64        # track number of times random opinions were adopted
    num_biased_opinions::Int64        # track number of times biased opinions were adopted
    num_random_strategies::Int64      # track number of times random strategies were adopted
    num_random_opinions_flipped::Int64        # track number of times random opinions resulted in opinion change
    num_biased_opinions_flipped::Int64        # track number of times biased opinions resulted in opinion change
    num_random_sa_flipped::Int64      # track number of times random strategies resulted in strategy change (sa)
    num_random_sd_flipped::Int64      # track number of times random strategies resulted in strategy change (sd)

    # constructor for the case where initial strategy distribution is specified
    function Population( sets::Sets, game::Game, initial_strategies::Array{Int64, 2}, verbose::Bool=false, sets_changed::Bool=true,    strategies_changed::Bool=true )

        # initializing the population with specified strategies
        payoffs      = zeros(Float64, sets.N)
        generation   = 0
        payoffs_mat  = zeros(Float64, (sets.N, sets.N))

        # conditional: 00 (DD), 01 (DC), 01 (CD), or 11 (CC) (in binary)
        strategies        = initial_strategies
        all_strategies    = collect(0:3)

        # misc
        prev_actions      = zeros(Int64, sets.M, sets.N, sets.N)
        all_opinions      = (-3^sets.M+1)÷2 : (3^sets.M-1)÷2
        all_memberships   = 1 : (3^sets.M-1)÷2
        strategies_base10 = [2*strategies[i,1] + strategies[i,2] for i in 1:sets.N]
        prev_interactions, prev_learner = 0, 0

        # overall trackers
        num_imitating, num_not_imitating                                = 0, 0
        num_strategies_changed, num_sets_changed, num_both_changed      = 0, 0, 0
        num_random_opinions, num_biased_opinions, num_random_strategies = 0, 0, 0
        num_random_opinions_flipped, num_biased_opinions_flipped        = 0, 0
        num_random_sa_flipped, num_random_sd_flipped                    = 0, 0

        # initialize sim::Simulation
        sim = Simulation(sets)

        return new(sets, game, sim, strategies, payoffs, all_strategies,
            generation, verbose, sets_changed, strategies_changed,
            prev_actions, prev_learner, payoffs_mat,
            all_opinions, all_memberships, strategies_base10, prev_interactions,
            num_imitating, num_not_imitating, num_strategies_changed, num_sets_changed,
            num_both_changed, num_random_opinions, num_biased_opinions, num_random_strategies,
            num_random_opinions_flipped, num_biased_opinions_flipped,
            num_random_sa_flipped, num_random_sd_flipped)
    end

    # same constructor as above, but allows randomized initial strategies
    function Population( sets::Sets, game::Game, verbose::Bool=false, sets_changed::Bool=true, strategies_changed::Bool=true )

        # begin by initializing the population with random strategies
        payoffs      = zeros(Float64, sets.N)
        generation   = 0
        payoffs_mat  = zeros(Float64, (sets.N, sets.N))

        # conditional: 00 (DD), 01 (DC), 01 (CD), or 11 (CC) (in binary)
        strategies        = rand(collect(0:1), sets.N, 2)  # each strat is 0 or 1, total 4
        all_strategies    = collect(0:3)

        # misc
        prev_actions      = zeros(Int64, sets.M, sets.N, sets.N)
        all_opinions      = (-3^sets.M+1)÷2 : (3^sets.M-1)÷2
        all_memberships   = 1 : (3^sets.M-1)÷2
        strategies_base10 = [2*strategies[i,1] + strategies[i,2] for i in 1:sets.N]
        prev_interactions, prev_learner = 0, 0

        # overall trackers
        num_imitating, num_not_imitating                                = 0, 0
        num_strategies_changed, num_sets_changed, num_both_changed      = 0, 0, 0
        num_random_opinions, num_biased_opinions, num_random_strategies = 0, 0, 0
        num_random_opinions_flipped, num_biased_opinions_flipped        = 0, 0
        num_random_sa_flipped, num_random_sd_flipped                    = 0, 0

        # initialize sim::Simulation
        sim = Simulation(sets)

        return new(sets, game, sim, strategies, payoffs, all_strategies,
            generation, verbose, sets_changed, strategies_changed,
            prev_actions, prev_learner, payoffs_mat,
            all_opinions, all_memberships, strategies_base10, prev_interactions,
            num_imitating, num_not_imitating, num_strategies_changed, num_sets_changed,
            num_both_changed, num_random_opinions, num_biased_opinions, num_random_strategies,
            num_random_opinions_flipped, num_biased_opinions_flipped,
            num_random_sa_flipped, num_random_sd_flipped)
    end
end


text:
# Multivariate Normal distribution

import StatsFuns
using LinearAlgebra
export MvNormal
using Random
import Base
using Tullio

struct MvNormal{N, T, I, J} <: ParameterizedMeasure{N}
    par::NamedTuple{N, T}
end

function MvNormal(nt::NamedTuple{N,T}) where {N,T}
    I,J = mvnormaldims(nt)

    cache = Vector{Float64}(undef, max(I,J))
    MvNormal{N,T,I,J}(cache, nt)
end

function Base.size(d::MvNormal{N, T, I, J}) where {N,T,I,J}
    return (I,)
end

mvnormaldims(nt::NamedTuple{(:A, :b)}) = size(nt.A)

function MeasureTheory.basemeasure(μ::MvNormal{N, T, I,I}) where {N, T, I}
    return (1 / sqrt2π)^I * Lebesgue(ℝ)^I
end

sampletype(d::MvNormal{N, T, I, J}) where {N,T,I,J} = Vector{Float64}

MvNormal(; kwargs...) = begin
    MvNormal((; kwargs...))
end

function Random.rand!(rng::AbstractRNG, d::MvNormal{(:A, :b),T,I,J}, x::AbstractArray) where {T,I,J}
    z = getcache(d, J)
    rand!(rng, Normal()^J, z)
    copyto!(x, d.b)
    mul!(x, d.A, z, 1.0, 1.0)
    return x
end

function logdensity(d::MvNormal{(:A,:b)}, x)
    cache = getcache(d, size(d))
    z = d.A \ (x - d.b)
    return logdensity(MvNormal(), z) - logabsdet(d.A)
end

≪(::MvNormal, ::Lebesgue{ℝ}) = true

function logdensity(d::MvNormal{(:Σ⁻¹,)}, x)
    @tullio ℓ = -0.5 * x[i] * d.Σ⁻¹[i,j] * x[j]
    return ℓ
end

mvnormaldims(nt::NamedTuple{(:Σ⁻¹,)}) = size(nt.Σ⁻¹)


text:
export load, save

"""
    load(filename::AbstractString)

Read a BVH file with name `filename` and return a `BVHGraph`.

# Examples
```julia
julia> g = load("Test.bvh")
BVHGraph
Name: Test.bvh
[...]
```
"""
function load(filename::AbstractString)
    list = read(filename, String) |> split((' ', '\t', '\n', '\r'))
    g = BVHGraph(1, name = filename,
                    offset = [parse(Float64, e) for e in list[6:8]],
                    sequence = shorten(list[11:13]))

    name!(g, 1, "ROOT" * ' ' * list[3])
    sequence!(g, 1, shorten(list[14:16]))
    i = 17

    function add_joint(v₋₁::Integer)
        while true
            v = add_vertex!(g, name = list[i] * ' ' * list[i + 1])
            add_edge!(g, v₋₁, v, offset = [parse(Float64, e) for e in list[i + 4:i + 6]])

            if list[i] == "JOINT"
                sequence!(g, v, shorten(list[i + 9:i + 11]))
                i += 12
                add_joint(v)
                i += 1
            else
                i += 8
            end

            list[i] != "JOINT" && list[i] != "End" && break
        end
    end

    add_joint(1)
    frames = parse(Int64, list[i + 3])
    frames!(g, frames)
    frametime!(g, parse(Float64, list[i + 6]))
    positions!(g, zeros(Float64, frames, 3))

    for v in vertices(g)
        outneighbors(g, v) != [] && rotations!(g, v, zeros(Float64, frames, 3))
    end

    i += 7

    function add_frame(v::Integer, f::Integer)
        rotations(g, v)[f, :] = [parse(Float64, e) for e in list[i:i + 2]]
        i += 3

        for n in outneighbors(g, v)
            outneighbors(g, n) != [] && add_frame(n, f)
        end
    end

    for f in 1:frames
        positions(g)[f, :] = [parse(Float64, e) for e in list[i:i + 2]]
        rotations(g, 1)[f, :] = [parse(Float64, e) for e in list[i + 3:i + 5]]
        i += 6

        for n in outneighbors(g, 1)
            outneighbors(g, n) != [] && add_frame(n, f)
        end
    end

    return g
end

"""
    save(g::BVHGraph, filename::AbstractString)

Save a BVHGraph `g` to `filename` as a BVH file.
"""
function save(g::BVHGraph, filename::AbstractString)

    function add_hierarchy(io::IO, v₋₁::Integer, v::Integer, tabs::Integer)
        println(io, '\t'^tabs, name(g, v))
        println(io, '\t'^tabs, "{")
        println(io, '\t'^(tabs + 1), "OFFSET", ['\t' * string(e) for e in offset(g, v₋₁, v)]...)

        if outneighbors(g, v) != []
            println(io, '\t'^(tabs + 1), "CHANNELS 3 ", sequence(g, v) |> extend_rotation)

            for n in outneighbors(g, v)
                add_hierarchy(io, v, n, tabs + 1)
            end
        end

        println(io, '\t'^tabs, "}")
    end

    function add_numbers(io::IO, f::Integer, v::Integer)
        print(io, ['\t' * string(e) for e in rotations(g, v, f)]...)

        for n in outneighbors(g, v)
            outneighbors(g, n) != [] && add_numbers(io, f, n)
        end
    end

    open(filename, "w") do io
        println(io, "HIERARCHY")
        println(io, name(g, 1))
        println(io, "{")
        println(io, "\tOFFSET", ['\t' * string(e) for e in offset(g)]...)
        println(io, "\tCHANNELS 6 ", sequence(g) |> extend_position, " ", sequence(g, 1) |> extend_rotation)

        for n in outneighbors(g, 1)
            add_hierarchy(io, 1, n, 1)
        end

        println(io, "}")
        println(io, "MOTION")
        println(io, "Frames: ", frames(g))
        println(io, "Frame Time: ", frametime(g))

        for f in 1:frames(g)
            pos = positions(g)[f, :]
            print(io, pos[1], "\t", pos[2], "\t", pos[3])
            print(io, ['\t' * string(e) for e in rotations(g, 1, f)]...)

            for n in outneighbors(g, 1)
                add_numbers(io, f, n)
            end

            print(io, "\n")
        end
    end
end


text:
const SKEW_4D = 0.309016994374947f0
const UNSKEW_4D = -0.138196601125011f0

const RSQUARED_4D = 4.0f0 / 5.0f0

const S2VAL = -0.211324865405187 # (1/sqrt(N+1)-1)/N where N=2
const WWMUL = 1.118033988749894
const WWVAL = -0.866025403784439
const YZVAL = 0.28867513459481294226
const INVSQRT3 = -0.5773502691896257 # (-1/sqrt(3))

export ImproveXYZ, ImproveXYZ_XY, ImproveXYZ_XZ, ImproveXY_ZW
struct ImproveXYZ_XY <: NoiseType end
struct ImproveXYZ_XZ <: NoiseType end
struct ImproveXYZ    <: NoiseType end
struct ImproveXY_ZW  <: NoiseType end

"""
4D SuperSimplex noise, with XYZ oriented like `noise(ImproveXY, ...)`
and W for an extra degree of freedom. W repeats eventually.
Recommended for time-varied animations which texture a 3D object (W=time)
in a space where Z is vertical
"""
function noise(::Type{ImproveXYZ_XY}, seed::Int64, x, y, z, w)
    xy = x + y
    s2 = xy * S2VAL
    zz = z * YZVAL
    ww = w * WWMUL
    r = zz + ww + s2
    _noise(seed, x + r, y + r, xy * INVSQRT3 + (zz + ww), z * WWVAL + ww)
end

"""
4D SuperSimplex noise, with XYZ oriented like `noise(ImproveXZ, ...)`
and W for an extra degree of freedom. W repeats eventually.
Recommended for time-varied animations which texture a 3D object (W=time)
in a space where Y is vertical
"""
function noise(::Type{ImproveXYZ_XZ}, seed::Int64, x, y, z, w)
    xz = x + z
    s2 = xz * S2VAL
    yy = y * YZVAL
    ww = w * WWMUL
    r = yy + ww + s2
    _noise(seed, x + r, xz * INVSQRT3 + (yy + ww), z + r, y * WWVAL + ww)
end

"""
4D SuperSimplex noise, with XYZ oriented like 3D `noise`
and W for an extra degree of freedom. W repeats eventually.
Recommended for time-varied animations which texture a 3D object (W=time)
where there isn't a clear distinction between horizontal and vertical
"""
function noise(::Type{ImproveXYZ}, seed::Int64, x, y, z, w)
    xyz = x + y + z
    ww = w * WWMUL
    s2 = xyz * -0.16666666666666666 + ww
    _noise(seed, x + s2, y + s2, z + s2, -0.5 * xyz + ww)
end

"""
4D SuperSimplex noise, with XY and ZW forming orthogonal triangular-based planes.
Recommended for 3D terrain, where X and Y (or Z and W) are horizontal.
Recommended for `noise(seed, x, y, sin(time), cos(time))` trick.
"""
function noise(::Type{ImproveXY_ZW}, seed::Int64, x, y, z, w)
    s2 = (x + y) * -0.28522513987434876941 + (z + w) * 0.83897065470611435718
    t2 = (z + w) * 0.21939749883706435719 + (x + y) * -0.48214856493302476942
    return _noise(seed, x + s2, y + s2, z + t2, w + t2)
end

"""
4D SuperSimplex noise, fallback lattice orientation.
"""
function noise(seed::Int64, x, y, z, w)
    # Get points for A4 lattice
    s = SKEW_4D * (x + y + z + w)
    return _noise(seed, x + s, y + s, z + s, w + s)
end

"""
4D SuperSimplex noise base.
Using ultra-simple 4x4x4x4 lookup partitioning.
This isn't as elegant or SIMD/GPU/etc. portable as other approaches,
but it competes performance-wise with optimized 2014 OpenSimplex.
"""
function _noise(seed::Int64, xs, ys, zs, ws)
    # Get base points and offsets
    xsb, ysb, zsb, wsb = fastFloor(xs), fastFloor(ys), fastFloor(zs), fastFloor(ws)
    xsi, ysi, zsi, wsi = Float32(xs - xsb), Float32(ys - ysb), Float32(zs - zsb), Float32(ws - wsb)

    # Unskewed offsets
    ssi = (xsi + ysi + zsi + wsi) * UNSKEW_4D
    xi, yi, zi, wi = xsi + ssi, ysi + ssi, zsi + ssi, wsi + ssi

    # Prime pre-multiplication for hash.
    xsvp = (xsb * PRIME_X)%Int64
    ysvp = (ysb * PRIME_Y)%Int64
    zsvp = (zsb * PRIME_Z)%Int64
    wsvp = (wsb * PRIME_W)%Int64

    # Index into initial table.
    index = ((fastFloor(xs * 4) & 3) << 0) |
            ((fastFloor(ys * 4) & 3) << 2) |
            ((fastFloor(zs * 4) & 3) << 4) |
            ((fastFloor(ws * 4) & 3) << 6)

    # Point contributions
    value = 0f0
    ind_beg = index == 0 ? 1 : Int(offsets_4D[index])
    ind_end = Int(offsets_4D[index+1])
    for i = ind_beg:ind_end
        c = lookup_4D[vertex_ind[i]]
        dx, dy, dz, dw = xi + c.dx, yi + c.dy, zi + c.dz, wi + c.dw
        a = (dx * dx + dy * dy) + (dz * dz + dw * dw)
        if a < RSQUARED_4D
            value += grad(a - RSQUARED_4D, seed,
                          xsvp + c.xsvp,
                          ysvp + c.ysvp,
                          zsvp + c.zsvp,
                          wsvp + c.wsvp,
                          dx, dy, dz, dw)
        end
    end
    return Float32(value)
end


text:
"""
	xyz2grd(cmd0::String="", arg1=[]; kwargs...)

Convert data table to a grid file.

Full option list at [`xyz2grd`](http://gmt.soest.hawaii.edu/doc/latest/xyz2grd.html)

Parameters
----------

- **I** : **inc** : -- Str or Number --     Flags = xinc[unit][+e|n][/yinc[unit][+e|n]]

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/xyz2grd.html#i)
- $(GMT.opt_R)
- **A** : **multiple_nodes** : -- Str --      Flags = [d|f|l|m|n|r|S|s|u|z]

    By default we will calculate mean values if multiple entries fall on the same node.
    Use A to change this behavior.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/xyz2grd.html#a)
- **D** : **header** : -- Str --  Flags = [+xxname][+yyname][+zzname][+sscale][+ooffset][+ninvalid][+ttitle][+rremark]

    Output edges
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/xyz2grd.html#d)
- **G** : **outgrid** : -- Str --

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdclip(....) form.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/xyz2grd.html#g)
- $(GMT.opt_J)
- **S** : **swap** : -- Str or [] --        Flags = [zfile]

    Swap the byte-order of the input only. No grid file is produced.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/xyz2grd.html#s)
- $(GMT.opt_V)
- **Z** : **flags** : -- Str --

    Read a 1-column table. This assumes that all the nodes are present and sorted according to specified ordering convention contained in. ``flags``.
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/grd2xyz.html#z)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function xyz2grd(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("xyz2grd", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	cmd, = parse_R("", d)
	cmd, = parse_J(cmd, d)
	cmd  = parse_V_params(cmd, d)
	cmd, = parse_bi(cmd, d)
	cmd, = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_r(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'A', d, [:A :multiple_nodes])
	cmd = add_opt(cmd, 'D', d, [:D :header])
    cmd = add_opt(cmd, 'G', d, [:G :outgrid])
	cmd = add_opt(cmd, 'I', d, [:I :inc])
	cmd = add_opt(cmd, 'S', d, [:S :swap])
	cmd = add_opt(cmd, 'Z', d, [:Z :flags])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "xyz2grd", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
xyz2grd(arg1::Array, cmd0::String=""; kw...) = xyz2grd(cmd0, arg1; kw...)


text:
struct Empatica <: AbstractVendor end
struct E4 <: AbstractProduct end
struct SkinConductance <: AbstractData end
struct Tags <: AbstractData end


text:
"""
```
splitNode!(network::DislocationNetwork, splitNode, splitConnect, midCoord, midVel)
```
Splits node `splitNode` along connection `splitConnect`, puts it at coordinate `midCoord` with velocity `midVel`. If it is called from within [`refineNetwork!`](@ref), `midCoord` and `midVel` are the coordinate between `splitNode` and the node connected to it via `splitConnect` and gives it the mean velocity of the two nodes.
"""
function splitNode!(network::DislocationNetwork, splitNode, splitConnect, midCoord, midVel)
    elemT = eltype(network.bVec)
    # newNode gets inserted between splitNode and the node it is connected to via the connection splitConnect. We want to take this connection and remove it from splitNode. We then connect splitNode to newNode. Then we assign splitConnect to newNode so that it can connect to the node splitNode used to be connected to, this way we close the loop and the connections move from splitNode -> other, to splitNode -> newNode -> other. We copy the connection to a temporary variable, guaranteeing the data isn't modified by removeConnection!().
    tmpConnect = network.connectivity[:, splitNode]
    # Remove connection from node in preparation of adding the new node.
    removeConnection!(network, splitNode, splitConnect)

    # New node created at the very end.
    network.numNode[1] += 1
    newNode = network.numNode[1]

    # Allocate memory.
    if newNode > length(network.label)
        # A good heuristic for memory allocation is to allocate an extra N log₂(N) entries.
        numNewEntries = Int(round(newNode * log2(newNode)))
        network = DislocationNetwork(;
            links = network.links,
            slipPlane = network.slipPlane,
            bVec = network.bVec,
            coord = hcat(
                network.coord,
                zeros(elemT, size(network.coord, 1), numNewEntries),
            ),
            label = vcat(network.label, zeros(nodeTypeDln, numNewEntries)),
            nodeVel = hcat(
                network.nodeVel,
                zeros(elemT, size(network.nodeVel, 1), numNewEntries),
            ),
            nodeForce = hcat(
                network.nodeForce,
                zeros(elemT, size(network.nodeForce, 1), numNewEntries),
            ),
            numNode = network.numNode,
            numSeg = network.numSeg,
            maxConnect = network.maxConnect,
            connectivity = hcat(
                network.connectivity,
                zeros(Int, size(network.connectivity, 1), numNewEntries),
            ),
            linksConnect = network.linksConnect,
            extSeg = network.extSeg,
            segIdx = network.segIdx,
            segForce = network.segForce,
        )
    end

    links = network.links
    bVec = network.bVec
    network.coord[:, newNode] = midCoord
    network.label[newNode] = 1
    network.nodeVel[:, newNode] = midVel
    connectivity = network.connectivity
    linksConnect = network.linksConnect

    # Create connectivity for the new node and update links and linksConnect.
    # This only does the first connection. The rest are made after in case they are needed.
    connectivity[1, newNode] = 1
    connectivity[2, newNode] = tmpConnect[(2 * splitConnect)]
    connectivity[3, newNode] = tmpConnect[(2 * splitConnect + 1)]

    link = connectivity[2, newNode]
    colLink = connectivity[3, newNode]
    links[colLink, link] = newNode
    linksConnect[colLink, link] = 1

    # Check if we need a new link between newNode and splitNode for Burgers vector conservation.
    b = SVector{3, elemT}(bVec[1, link], bVec[2, link], bVec[3, link])
    # If burgers vector is conserved return.
    b ⋅ b == 0 && return network

    # New link created at the end.
    network.numSeg[1] += 1
    newSeg = network.numSeg[1]

    # Allocate memory.
    if newSeg > size(network.links, 2)
        # A good heuristic for memory allocation is to allocate an extra N log₂(N) entries.
        numNewEntries = Int(round(newSeg * log2(newSeg)))

        network = DislocationNetwork(;
            links = hcat(network.links, zeros(Int, size(network.links, 1), numNewEntries)),
            slipPlane = hcat(
                network.slipPlane,
                zeros(elemT, size(network.slipPlane, 1), numNewEntries),
            ),
            bVec = hcat(network.bVec, zeros(elemT, size(network.bVec, 1), numNewEntries)),
            coord = network.coord,
            label = network.label,
            nodeVel = network.nodeVel,
            nodeForce = network.nodeForce,
            numNode = network.numNode,
            numSeg = network.numSeg,
            maxConnect = network.maxConnect,
            connectivity = network.connectivity,
            linksConnect = hcat(
                network.linksConnect,
                zeros(Int, size(network.linksConnect, 1), numNewEntries),
            ),
            extSeg = vcat(network.extSeg, zeros(Bool, numNewEntries)),
            segIdx = vcat(network.segIdx, zeros(Int, numNewEntries, 3)),
            segForce = cat(
                network.segForce,
                zeros(elemT, size(network.segForce, 1), 2, numNewEntries),
                dims = 3,
            ),
        )
    end

    bVec = network.bVec
    links = network.links
    slipPlane = network.slipPlane
    connectivity = network.connectivity
    linksConnect = network.linksConnect

    # Add new link to network.
    bVec[:, newSeg] = b
    links[colLink, newSeg] = splitNode  # First node in the link is the node that was split.
    colOppLink = 3 - colLink            # Other column in links.
    links[colOppLink, newSeg] = newNode # Second node of the link is the new node.

    # Update connectivity of splitNode.
    newConnect1 = connectivity[1, splitNode] + 1
    connectivity[1, splitNode] = newConnect1
    connectivity[(2 * newConnect1):(2 * newConnect1 + 1), splitNode] .= (newSeg, colLink)

    # Update connectivity of newNode.
    newConnect2 = connectivity[1, newNode] + 1
    connectivity[1, newNode] = newConnect2
    connectivity[(2 * newConnect2):(2 * newConnect2 + 1), newNode] .= (newSeg, colOppLink)

    # Update linksConnect.
    linksConnect[colLink, newSeg] = newConnect1
    linksConnect[colOppLink, newSeg] = newConnect2

    # Fix slip plane.
    # Line direction and velocity of the resultant dislocation.
    coord = network.coord
    nodeVel = network.nodeVel

    # WARNING This calculation's dodgy. Try using the cross product of the adjacent segments.
    t = SVector{3, elemT}(
        coord[1, splitNode] - coord[1, newNode],
        coord[2, splitNode] - coord[2, newNode],
        coord[3, splitNode] - coord[3, newNode],
    )

    v = SVector{3, elemT}(
        nodeVel[1, splitNode] + nodeVel[1, newNode],
        nodeVel[2, splitNode] + nodeVel[2, newNode],
        nodeVel[3, splitNode] + nodeVel[3, newNode],
    )

    # Potential new slip plane.
    n1 = t × b  # For non-screw segments.
    n2 = t × v  # For screw segments.
    if n1 ⋅ n1 > eps(elemT) # non-screw
        slipPlane[:, newSeg] = n1 / norm(n1)
    elseif n2 ⋅ n2 > eps(elemT) # screw
        slipPlane[:, newSeg] = n2 / norm(n2)
    end

    return network
end

"""
```
refineNetwork!(
    dlnParams::DislocationParameters,
    matParams::MaterialParameters,
    mesh::AbstractMesh,
    forceDisplacement::ForceDisplacement,
    network::DislocationNetwork,
)
```
Refines a dislocation network to ensure all the segments are shorter than the maximum allowable length and so no two links form a triangle with an area over the maximum allowed.
"""
function refineNetwork!(
    dlnParams::DislocationParameters,
    matParams::MaterialParameters,
    mesh::AbstractMesh,
    forceDisplacement::ForceDisplacement,
    network::DislocationNetwork,
)
    maxAreaSq = dlnParams.maxAreaSq
    maxSegLen = dlnParams.maxSegLen
    twoMinSegLen = dlnParams.twoMinSegLen

    links = network.links
    coord = network.coord
    label = network.label
    numNode = network.numNode[1]
    nodeVel = network.nodeVel
    connectivity = network.connectivity

    elemT = eltype(network.coord)

    @inbounds for i in 1:numNode
        if connectivity[1, i] == 2 && label[i] == intMobDln
            link1, missing, link1_nodeOppI = findConnectedNode(network, i, 1) # First connection.
            link2, missing, link2_nodeOppI = findConnectedNode(network, i, 2) # Second connection.

            # Create triangle formed by the node and its two links.
            iCoord = SVector{3, elemT}(coord[1, i], coord[2, i], coord[3, i])
            # Side 1
            coordVec1 =
                SVector{3, elemT}(
                    coord[1, link1_nodeOppI],
                    coord[2, link1_nodeOppI],
                    coord[3, link1_nodeOppI],
                ) - iCoord
            # Side 2
            coordVec2 =
                SVector{3, elemT}(
                    coord[1, link2_nodeOppI],
                    coord[2, link2_nodeOppI],
                    coord[3, link2_nodeOppI],
                ) - iCoord
            # Side 3 (close the triangle)
            coordVec3 = coordVec2 - coordVec1
            # Side lengths.
            r1 = norm(coordVec1)
            r2 = norm(coordVec2)
            r3 = norm(coordVec3)

            # Half the triangle's perimeter.
            r0 = (r1 + r2 + r3) / 2
            # Heron's formula for the area of a triangle when we know its sides.
            areaSq = r0 * (r0 - r1) * (r0 - r2) * (r0 - r3)

            # Check if we have to split the second link.
            # Split node the area is greater than the maximum area and if by splitting we would not create two segments smaller than the minimum allowed and if the node is in the simulation. Or if r2 is longer than the maximum length allowed.
            if (areaSq > maxAreaSq && r2 >= twoMinSegLen && link2_nodeOppI <= numNode) ||
               r2 > maxSegLen
                midCoord =
                    SVector{3, elemT}(
                        coord[1, i] + coord[1, link2_nodeOppI],
                        coord[2, i] + coord[2, link2_nodeOppI],
                        coord[3, i] + coord[3, link2_nodeOppI],
                    ) / 2

                midVel =
                    SVector{3, elemT}(
                        nodeVel[1, i] + nodeVel[1, link2_nodeOppI],
                        nodeVel[2, i] + nodeVel[2, link2_nodeOppI],
                        nodeVel[3, i] + nodeVel[3, link2_nodeOppI],
                    ) / 2

                network = splitNode!(network, i, 2, midCoord, midVel)
                getSegmentIdx!(network)
                links = network.links
                slipPlane = network.slipPlane
                coord = network.coord
                label = network.label
                nodeVel = network.nodeVel
                connectivity = network.connectivity

                newNode = network.numNode[1]
                newLink = network.numSeg[1]

                equalSlipPlane = let
                    flag = true
                    @inbounds @simd for j in 1:3
                        flag = flag && slipPlane[j, link2] ≈ slipPlane[j, link1]
                    end
                    flag
                end
                if equalSlipPlane
                    @inbounds @simd for j in 1:3
                        slipPlane[j, newLink] = slipPlane[j, link2]
                    end
                end

                # Calculate force and mobility for the new node's connectivity.
                j = 1:connectivity[1, newNode]
                link, missing, oldNode = findConnectedNode(network, newNode, j) # j connections.

                # Calculate segment force for segment link.
                calcSegForce!(dlnParams, matParams, mesh, forceDisplacement, network, link)

                # Calculate old node velocity.
                dlnMobility!(dlnParams, matParams, network, oldNode)
                nodeVel = network.nodeVel

                # Calculate new node velocity.
                dlnMobility!(dlnParams, matParams, network, newNode)
                nodeVel = network.nodeVel
            end

            # Check if we have to split the first link.
            # Split node the area is greater than the maximum area and if by splitting we would not create two segments smaller than the minimum allowed and if the node is in the simulation. Or if r1 is longer than the maximum length allowed.
            if (areaSq > maxAreaSq && r1 >= twoMinSegLen && link1_nodeOppI <= numNode) ||
               r1 > maxSegLen
                midCoord =
                    SVector{3, elemT}(
                        coord[1, i] + coord[1, link1_nodeOppI],
                        coord[2, i] + coord[2, link1_nodeOppI],
                        coord[3, i] + coord[3, link1_nodeOppI],
                    ) / 2

                midVel =
                    SVector{3, elemT}(
                        nodeVel[1, i] + nodeVel[1, link1_nodeOppI],
                        nodeVel[2, i] + nodeVel[2, link1_nodeOppI],
                        nodeVel[3, i] + nodeVel[3, link1_nodeOppI],
                    ) / 2

                network = splitNode!(network, i, 1, midCoord, midVel)
                getSegmentIdx!(network)
                links = network.links
                slipPlane = network.slipPlane
                coord = network.coord
                label = network.label
                nodeVel = network.nodeVel
                connectivity = network.connectivity

                newNode = network.numNode[1]
                newLink = network.numSeg[1]

                equalSlipPlane = let
                    flag = true
                    @inbounds @simd for j in 1:3
                        flag = flag && slipPlane[j, link1] ≈ slipPlane[j, link2]
                    end
                    flag
                end
                if equalSlipPlane
                    @inbounds @simd for j in 1:3
                        slipPlane[j, newLink] = slipPlane[j, link1]
                    end
                end

                # Calculate force and mobility for the new node's connectivity.
                j = 1:connectivity[1, newNode]
                link, oldNode = findConnectedNode(network, newNode, j) # j connections.

                # Calculate segment force for segment link.
                calcSegForce!(dlnParams, matParams, mesh, forceDisplacement, network, link)

                # Calculate old node velocity.
                dlnMobility!(dlnParams, matParams, network, oldNode)
                nodeVel = network.nodeVel

                # Calculate new node velocity.
                dlnMobility!(dlnParams, matParams, network, newNode)
                nodeVel = network.nodeVel
            end

        elseif connectivity[1, i] > 2 && label[i] == intMobDln
            # Loop through the connections of node i.
            for j in 1:connectivity[1, i]
                link, missing, link_nodeOpp = findConnectedNode(network, i, j) # j connection.

                t = SVector{3, elemT}(
                    coord[1, link_nodeOpp] - coord[1, i],
                    coord[2, link_nodeOpp] - coord[2, i],
                    coord[3, link_nodeOpp] - coord[3, i],
                )
                r1 = norm(t)

                # If the link is smaller than the maximum segment length we skip to the next iteration.
                r1 < maxSegLen && continue

                midCoord =
                    SVector{3, elemT}(
                        coord[1, i] + coord[1, link_nodeOpp],
                        coord[2, i] + coord[2, link_nodeOpp],
                        coord[3, i] + coord[3, link_nodeOpp],
                    ) / 2

                midVel =
                    SVector{3, elemT}(
                        nodeVel[1, i] + nodeVel[1, link_nodeOpp],
                        nodeVel[2, i] + nodeVel[2, link_nodeOpp],
                        nodeVel[3, i] + nodeVel[3, link_nodeOpp],
                    ) / 2

                network = splitNode!(network, i, j, midCoord, midVel)
                getSegmentIdx!(network)
                links = network.links
                slipPlane = network.slipPlane
                coord = network.coord
                label = network.label
                nodeVel = network.nodeVel
                connectivity = network.connectivity

                newNode = network.numNode[1]
                newLink = network.numSeg[1]

                @inbounds @simd for k in 1:3
                    slipPlane[k, newLink] = slipPlane[k, link]
                end

                # Calculate force and mobility for the new node's connectivity.
                k = 1:connectivity[1, newNode]
                link, missing, oldNode = findConnectedNode(network, newNode, k) # k connections.

                # Calculate segment force for segment link.
                calcSegForce!(dlnParams, matParams, mesh, forceDisplacement, network, link)
                # Calculate old node velocity.
                dlnMobility!(dlnParams, matParams, network, oldNode)
                nodeVel = network.nodeVel

                # Calculate new node velocity.
                dlnMobility!(dlnParams, matParams, network, newNode)
                nodeVel = network.nodeVel
            end
        end
    end
    return network
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
# check sparse matrix construction
@test isequal(full(sparse(complex(ones(5,5),ones(5,5)))), complex(ones(5,5),ones(5,5)))

# check matrix operations
se33 = speye(3)
do33 = ones(3)
@test isequal(se33 * se33, se33)

# check sparse binary op
@test all(full(se33 + convert(SparseMatrixCSC{Float32,Int32}, se33)) == 2*eye(3))
@test all(full(se33 * convert(SparseMatrixCSC{Float32,Int32}, se33)) == eye(3))

# check horiz concatenation
@test all([se33 se33] == sparse([1, 2, 3, 1, 2, 3], [1, 2, 3, 4, 5, 6], ones(6)))

# check vert concatenation
@test all([se33; se33] == sparse([1, 4, 2, 5, 3, 6], [1, 1, 2, 2, 3, 3], ones(6)))

# check h+v concatenation
se44 = speye(4)
sz42 = spzeros(4, 2)
sz41 = spzeros(4, 1)
sz34 = spzeros(3, 4)
se77 = speye(7)
@test all([se44 sz42 sz41; sz34 se33] == se77)

# check blkdiag concatenation
@test all(blkdiag(se33, se33) == sparse([1, 2, 3, 4, 5, 6], [1, 2, 3, 4, 5, 6], ones(6)))

# check concatenation promotion
sz41_f32 = spzeros(Float32, 4, 1)
se33_i32 = speye(Int32, 3, 3)
@test all([se44 sz42 sz41_f32; sz34 se33_i32] == se77)

# check mixed sparse-dense concatenation
sz33 = spzeros(3, 3)
de33 = eye(3)
@test  all([se33 de33; sz33 se33] == full([se33 se33; sz33 se33 ]))

# check splicing + concatenation on
# random instances, with nested vcat
# also side-checks sparse ref
for i = 1 : 10
    a = sprand(5, 4, 0.5)
    @test all([a[1:2,1:2] a[1:2,3:4]; a[3:5,1] [a[3:4,2:4]; a[5,2:4]]] == a)
end

# sparse ref
a116 = reshape(1:16, 4, 4)
s116 = sparse(a116)
p = [4, 1, 2, 3, 2]
@test full(s116[p,:]) == a116[p,:]
@test full(s116[:,p]) == a116[:,p]
@test full(s116[p,p]) == a116[p,p]

# sparse assign
p = [4, 1, 3]
a116[p, p] = -1
s116[p, p] = -1
@test a116 == s116

p = [2, 1, 4]
a116[p, p] = reshape(1:9, 3, 3)
s116[p, p] = reshape(1:9, 3, 3)
@test a116 == s116

# matrix-vector multiplication (non-square)
for i = 1:5
    a = sprand(10, 5, 0.5)
    b = rand(5)
    @test maximum(abs(a*b - full(a)*b)) < 100*eps()
end

# complex matrix-vector multiplication and left-division
for i = 1:5
    a = speye(5) + 0.1*sprandn(5, 5, 0.2)
    b = randn(5,3) + im*randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + 0.1*sprandn(5, 5, 0.2) + 0.1*im*sprandn(5, 5, 0.2)
    b = randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + tril(0.1*sprandn(5, 5, 0.2))
    b = randn(5,3) + im*randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + tril(0.1*sprandn(5, 5, 0.2) + 0.1*im*sprandn(5, 5, 0.2))
    b = randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + triu(0.1*sprandn(5, 5, 0.2))
    b = randn(5,3) + im*randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + triu(0.1*sprandn(5, 5, 0.2) + 0.1*im*sprandn(5, 5, 0.2))
    b = randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + triu(0.1*sprandn(5, 5, 0.2))
    b = randn(5,3) + im*randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = spdiagm(randn(5)) + im*spdiagm(randn(5))
    b = randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    b = randn(5,3) + im*randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())
end

# matrix multiplication and kron
for i = 1:5
    a = sprand(10, 5, 0.7)
    b = sprand(5, 10, 0.3)
    @test maximum(abs(a*b - full(a)*full(b))) < 100*eps()
    @test full(kron(a,b)) == kron(full(a), full(b))
end

# reductions
@test sum(se33)[1] == 3.0
@test sum(se33, 1) == [1.0 1.0 1.0]
@test sum(se33, 2) == [1.0 1.0 1.0]'
@test prod(se33)[1] == 0.0
@test prod(se33, 1) == [0.0 0.0 0.0]
@test prod(se33, 2) == [0.0 0.0 0.0]'

# spdiagm
@test full(spdiagm((ones(2), ones(2)), (0, -1), 3, 3)) ==
                       [1.0  0.0  0.0; 1.0  1.0  0.0;  0.0  1.0  0.0]

# elimination tree
## upper triangle of the pattern test matrix from Figure 4.2 of
## "Direct Methods for Sparse Linear Systems" by Tim Davis, SIAM, 2006
rowval = int32([1,2,2,3,4,5,1,4,6,1,7,2,5,8,6,9,3,4,6,8,10,3,5,7,8,10,11])
colval = int32([1,2,3,3,4,5,6,6,6,7,7,8,8,8,9,9,10,10,10,10,10,11,11,11,11,11,11])
A = sparse(rowval, colval, ones(length(rowval)))
P,post = Base.LinAlg.etree(A, true)
@test P == int32([6,3,8,6,8,7,9,10,10,11,0])
@test post == int32([2,3,5,8,1,4,6,7,9,10,11])

# issue #4986, reinterpret
sfe22 = speye(Float64, 2)
mfe22 = eye(Float64, 2)
@test reinterpret(Int64, sfe22) == reinterpret(Int64, mfe22)

# issue #5190
@test_throws sparsevec([3,5,7],[0.1,0.0,3.2],4)

# issue #5169
@test nfilled(sparse([1,1],[1,2],[0.0,-0.0])) == 0

# issue #5386
K,J,V = findnz(SparseMatrixCSC(2,1,[1,3],[1,2],[1.0,0.0]))
@test length(K) == length(J) == length(V) == 1

# issue #5437
@test nfilled(sparse([1,2,3],[1,2,3],[0.0,1.0,2.0])) == 2

# issue #5824
@test sprand(4,5,0.5).^0 == sparse(ones(4,5))

# issue #5985
@test sprandbool(4, 5, 0.0) == sparse(zeros(Bool, 4, 5))
@test sprandbool(4, 5, 1.00) == sparse(ones(Bool, 4, 5))
sprb45 = sprandbool(4, 5, 0.5)
@test length(sprb45) == 20
@test (sum(sprb45)[1] - 10) <= 5

# issue #5853, sparse diff
for i=1:2, a={[1 2 3], [1 2 3]', speye(3)}
    @test all(diff(sparse(a),i) == diff(a,i))
end

# test for "access to undefined error" types that initially allocate elements as #undef
@test all(sparse(1:2, 1:2, Number[1,2])^2 == sparse(1:2, 1:2, [1,4]))
sd1 = diff(sparse([1,1,1], [1,2,3], Number[1,2,3]), 1)

# issue #6036
P = spzeros(Float64, 3, 3)
for i = 1:3
    P[i,i] = i
end
@test minimum(P) == [0]
@test maximum(P) == [3]

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
using MarketData, TimeSeries

facts("Candlesticks") do

    context("doji") do
        @fact ohlc[findall(doji(ohlc))].timestamp[1] --> Date(2000, 5, 30)
        @fact ohlc[findall(doji(ohlc))].timestamp[2] --> Date(2000, 12, 11)
        @fact length(ohlc[findall(doji(ohlc))])      --> 5
    end
end


text:
include("compute_hermite.jl")
function hermitekernel(n,x)
    #Compute the Hermite kernel matrix at points x
    #including degrees from 0 to n inclusive

    # # Method 1: Direct Computation
       phi=compute_hermite(n,x);
       K=phi'*phi;

    # # Method 2: Christoffel-Darboux
    #   phinp1 = compute_hermite(n+1,x,false)
    #   phin = compute_hermite(n,x,false)
    #   K = phinp1'*phin*sqrt((n+1)/2)
    #   x = x[:]*ones(1,length(x))
    #   K = (K-K')./(x-x')

    K
end


text:
horner2(coefs, x) = foldr((u, v) -> u + x * v, 0, coefs)

@show horner2([-19, 7, -4, 6], 3)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test, SafeTestsets

const GROUP = get(ENV, "GROUP", "All")

if GROUP == "All" || GROUP == "Tests"
    @safetestset "InpParser Tests" begin
        include("inp_parser/parser.jl")
    end
    @safetestset "TopOptProblems Tests" begin
        include("topopt_problems/problems.jl")
        include("topopt_problems/metadata.jl")
    end
    @safetestset "Functions" begin
        include("Functions/test_common_fns.jl")
        include("Functions/test_buckling_fns.jl")
    end
    @safetestset "Solver" begin
        include("fea/solvers.jl")
    end
    @safetestset "Truss Problem Tests" begin
        include("truss_topopt_problems/test_problem.jl")
        include("truss_topopt_problems/test_fea.jl")
        include("truss_topopt_problems/test_buckling.jl")
        include("truss_topopt_problems/test_buckling_optimize.jl")
    end
end

if GROUP == "All" || GROUP == "Examples"
    @safetestset "CSIMP example" begin
        include("examples/csimp.jl")
    end
    @safetestset "Global stress example" begin
        include("examples/global_stress.jl")
    end
    @safetestset "Local stress example" begin
        include("examples/local_stress.jl")
    end
    @safetestset "More examples" begin
        include("examples/test_examples.jl")
    end
    @safetestset "Neural network example" begin
        include("examples/neural.jl")
    end
    @safetestset "Integer nonlinear programming for truss optimization example" begin
        include("examples/mixed_integer_truss/truss_compliance_2d1.jl")
    end
end

if GROUP == "All" || GROUP == "WCSMO14"
    # This was originlly part of https://github.com/JuliaTopOpt/TopOpt.jl_WCSMO21
    @safetestset "Continuum demos" begin
        include("wcsmo14/demos/continuum/cont_compliance1.jl")
        include("wcsmo14/demos/continuum/cont_compliance2.jl")
        include("wcsmo14/demos/continuum/cont_stress.jl")
    end
    @safetestset "Truss 2d demos" begin
        include("wcsmo14/demos/truss/truss_compliance_2d1.jl")
        include("wcsmo14/demos/truss/truss_compliance_2d2.jl")
    end
    @safetestset "Truss 3d demos" begin
        include("wcsmo14/demos/truss/truss_compliance_3d1.jl")
        include("wcsmo14/demos/truss/truss_compliance_3d2.jl")
    end
    @safetestset "WCSMO Benchmarks" begin
        include("wcsmo14/jl_benchmarks/compare_neo99_2D.jl")
        #include("wcsmo14/jl_benchmarks/compare_polytop.jl")
        #include("wcsmo14/jl_benchmarks/compare_top3d.jl")
        #include("wcsmo14/jl_benchmarks/compare_top3d125.jl")
        #include("wcsmo14/jl_benchmarks/new_problems.jl")
    end
end


text:
module BenchIncrements

using BenchmarkTools
using RecordArrays

using ..Utils: clobber

@inline function unsafe_inc!(xs::RecordArray)
    @inbounds xs.a[end] += 1
end

@inline function unsafe_inc!(xs)
    x = @inbounds xs[end]
    x = merge(x, (; a = x.a + 1))
    @inbounds xs[end] = x
end

@noinline function repeat_inc!(xs, n = 2^10)
    isempty(xs) && return
    for _ in 1:n
        unsafe_inc!(xs)
        clobber()
    end
end

function setup()
    suite = BenchmarkGroup()

    suite["base"] = @benchmarkable(
        repeat_inc!(xs),
        setup = begin
            xs = fill((a = 0, b = 0, c = 0, d = 0), 10)
        end,
    )

    suite["rec"] = @benchmarkable(
        repeat_inc!(xs),
        setup = begin
            xs = RecordArrays.fill((a = 0, b = 0, c = 0, d = 0), 10)
        end,
    )

    return suite
end

function clear() end

end  # module


text:
#=
CRC and Hashing functions for Str types (to make compatible with String hashes)

Copyright 2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
=#

using MurmurHash3

# Support for higher performance hashing, while still compatible with hashed UTF8 String

mmhash128(str::Union{String, Str}, seed::UInt32) =
    @preserve str mmhash128_a(sizeof(str), pointer(str), seed)

is_aligned(pnt::Ptr) = (reinterpret(UInt, pnt) & (sizeof(UInt) - 1)%UInt) == 0

# Check alignment of substrings first
function mmhash128(str::SubString, seed::UInt32)
    @preserve str begin
        pnt = pointer(str)
        if is_aligned(pnt)
            mmhash128_a(sizeof(str), pnt, seed)
        elseif sizeof(Int) == 8
            mmhash128_c(str, seed)
        else
            s = string(str)
            @preserve s mmhash128_a(sizeof(s), pointer(s), seed)
        end
    end
end

# This are for debugging purposes, to compare against current C implementation, will be removed
_memhash(siz, ptr, seed) =
    ccall(Base.memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), ptr, siz, seed % UInt32)

# Optimized code for hashing empty string
_hash(seed)          = last(mmhash128_a(seed%UInt32)) + seed
# Optimized for hashing a UTF-8 compatible aligned string
_hash(str, seed)     = last(mmhash128(str, seed%UInt32)) + seed
# For hashing generic abstract strings as if UTF-8 encoded
@static if sizeof(Int) == 8
    _hash_abs(str, seed) = last(mmhash128_c(str, seed%UInt32)) + seed
else
    function _hash_abs(str, seed)
        s = string(str)
        @preserve s last(mmhash128_a(sizeof(s), pointer(s), seed%UInt32)) + seed
    end
end

hash(str::Union{S,SubString{S}}, seed::UInt) where {S<:Str} =
    isempty(str) ? _hash(seed + Base.memhash_seed) : _hash_abs(str, seed + Base.memhash_seed)

# Check for UTF-8 compatible (i.e. only ASCII)
function hash(str::Union{S,SubString{S}}, seed::UInt) where {S<:Str{LatinCSE}}
    seed += Base.memhash_seed
    isempty(str) ? _hash(seed) : (is_ascii(str) ? _hash(str, seed) : _hash_abs(str, seed))
end

# Directly calculate hash for "compatible" types

hash(str::Union{S,SubString{S}},
     seed::UInt) where {S<:Str{<:Union{ASCIICSE,UTF8CSE,Binary_CSEs}}} =
         isempty(str) ? _hash(seed + Base.memhash_seed) : _hash(str, seed + Base.memhash_seed)

# Use crc32c to make CRC32c of UTF8 view of string, for use with hashing
# where different string types are supposed to compare as ==

function utf8crc(str::Union{S,SubString{S}}, seed::UInt) where {S<:Str}
    (len = ncodeunits(str)) == 0 && return utf8crc(str, seed)
    @preserve str begin
        pnt = pointer(str)
        len, flags, num4byte, num3byte, num2byte, latin1 = count_chars(S, pnt, len)
        # could be UCS2, _UCS2, UTF32, _UTF32, Text2, Text4
        utf8crc((flags == 0
                 ? _cvtsize(UInt8, pnt, len)
                 : _encode_utf8(pnt, len += latin1 + num2byte + num3byte*2 + num4byte*3)),
                seed)
    end
end

function utf8crc(str::Union{S,SubString{S}},
                 seed::UInt) where {S<:Str{<:Union{Text1CSE,Latin_CSEs}}}
    (len = ncodeunits(str)) == 0 && return utf8crc(str, seed)
    @preserve str begin
        pnt = pointer(str)
        utf8crc((cnt = count_latin(len, pnt)) == 0 ? str : _latin_to_utf8(pnt, len + cnt), seed)
    end
end

function utf8crc(str::Union{S,SubString{S}}, seed::UInt) where {S<:Str{UTF16CSE}}
    (len = ncodeunits(str)) == 0 && return utf8crc(str, seed)
    @preserve str begin
        pnt = pointer(str)
        len, flags, num4byte, num3byte, num2byte, latin1 = count_chars(S, pnt, len)
        utf8crc((flags == 0
                 ? _cvtsize(UInt8, pnt, len)
                 : _cvt_16_to_utf8(S, pnt, len += latin1 + num2byte + num3byte*2 + num4byte*3)),
                seed)
    end
end

utf8crc(str::Union{S,SubString{S}},
    seed::UInt32=0%UInt32) where {S<:Str{<:Union{ASCIICSE,UTF8CSE,BinaryCSE}}} =
        unsafe_crc32c(pointer(str), sizeof(s) % Csize_t, seed)

# Optimize conversion to ASCII or UTF8 to calculate compatible hash value

function cvthash(str::Union{S,SubString{S}}, seed::UInt) where {S<:Str}
    seed += Base.memhash_seed
    (len = ncodeunits(str)) == 0 && return _hash(seed)
    @preserve str begin
        pnt = pointer(str)
        len, flags, num4byte, num3byte, num2byte, latin1 = count_chars(S, pnt, len)
        # could be UCS2, _UCS2, UTF32, _UTF32, Text2, Text4
        _hash((flags == 0
               ? _cvtsize(UInt8, pnt, len)
               : _encode_utf8(pnt, len + latin1 + num2byte + num3byte*2 + num4byte*3)),
              seed)
    end
end

function cvthash(str::Union{S,SubString{S}}, seed::UInt) where {S<:Str{<:Latin_CSEs}}
    seed += Base.memhash_seed
    (len = ncodeunits(str)) == 0 && return _hash(seed)
    @preserve str begin
        pnt = pointer(str)
        _hash((cnt = count_latin(len, pnt)) == 0 ? str : _latin_to_utf8(pnt, len + cnt), seed)
    end
end

function cvthash(str::Union{S,SubString{S}}, seed::UInt) where {S<:Str{UTF16CSE}}
    seed += Base.memhash_seed
    (len = ncodeunits(str)) == 0 && return _hash(seed)
    @preserve str begin
        pnt = pointer(str)
        len, flags, num4byte, num3byte, num2byte, latin1 = count_chars(S, pnt, len)
        _hash((flags == 0
               ? _cvtsize(UInt8, pnt, len)
               : _cvt_16_to_utf8(S, pnt, len + latin1 + num2byte + num3byte*2 + num4byte*3)),
              seed)
    end
end


text:
export Fill

doc"""
    Fill(x)

Fill initializer.
"""
struct Fill
    x
end

(f::Fill)(::Type{T}, dims::Int...) where T = fill(T(f.x), dims)