text:
"""
TT version of the restarted GMRES algorithm
"""

#takes a hessenberg matrix and returns q,H resp. orthogonal and upper triangular such that q^T H_out = H_in
function qr_hessenberg(H)
    m = size(H,2)
    q = Matrix{Float64}(I,m+1,m+1)
    T = H
    for i in 1:m
        R = Matrix{Float64}(I,m+1,m+1)
        s = T[i+1,i]/sqrt(T[i,i]^2+T[i+1,i]^2)
        c = T[i,i]/sqrt(T[i,i]^2+T[i+1,i]^2)
        R[i:(i+1),i:(i+1)] = [c s; -s c]
        q = R*q
        T = R*T
    end
    return q,T
end

function tt_gmres(A::ttoperator,b::ttvector,x0::ttvector;Imax=500,tol=1e-8,m=30,hist=false,γ_list=Float64[])
    V = Array{ttvector}(undef,m)
    W = Array{ttvector}(undef,m)
    H = zeros(m+1,m)
    r0 = tt_compression_par(tt_add(b,mult_a_tt(-1.0,tt_compression_par(mult(A,x0)))))
    β = sqrt(tt_dot(r0,r0))
    V[1] = mult_a_tt(1/β,r0)
    W[1] = tt_compression_par(mult(A,V[1]))
    H[1,1] = tt_dot(W[1],V[1])
    W[1] = tt_compression_par(tt_add(W[1],mult_a_tt(-H[1,1],V[1])))
    H[2,1] = sqrt(tt_dot(W[1],W[1]))
    q,r = qr_hessenberg(H[1:2,1])
    γ = abs(β*q[2,1]) #γ = \| Ax_j-b \|_2
    if hist
        γ_list = vcat(γ_list,γ)
    end
    j = 1
    if Imax <=0 || isapprox(H[j+1,j],0.,atol=tol) || isapprox(γ,0,atol=tol)
        if hist
            return x0,γ_list,H[2,1]
        else
            return x0,[γ],H[2,1]
        end
    else
        while j <= min(m-1,Imax) && !isapprox(H[j+1,j],0.,atol=tol) && !isapprox(γ,0,atol=tol)
            V[j+1] = mult_a_tt(1/H[j+1,j],W[j])
            j+=1
            W[j] = mult(A,V[j])
            W[j] = tt_compression_par(W[j])
            for i in 1:j
                H[i,j] = tt_dot(W[j],V[i])
                W[j] = tt_add(W[j],mult_a_tt(-H[i,j],V[i]))
                W[j] = tt_compression_par(W[j])
            end
            H[j+1,j] = sqrt(tt_dot(W[j],W[j]))
            q,r = qr_hessenberg(@view H[1:(j+1),1:j])
            γ = abs(β*q[j+1,1])
            if hist
                γ_list = vcat(γ_list,γ)
            end
        end
        z = r[1:j,1:j]\q[1:j,1]
        for i in 1:j
            x0 = tt_add(x0,mult_a_tt(β*z[i],V[i]))
        end
        x0 = tt_compression_par(x0)
        return tt_gmres(A,b,x0,tol=tol,Imax=Imax-j,m=m,hist=hist,γ_list=γ_list)
    end
end

function tt_cg(A::ttoperator,b::ttvector,x0::ttvector;Imax=500,tol=1e-8)
    p = tt_compression_par(tt_add(b,mult_a_tt(-1.0,tt_compression_par(mult(A,x0)))))
    r = p
    j=1
    res= zeros(Imax)
    res[1] = sqrt(tt_dot(p,p))
    while j < Imax && res[j]>tol
        Ap = mult(A,p)
        Ap = tt_compression_par(Ap)
        a = res[j]^2/tt_dot(p,Ap)
        x0 = tt_add(x0,mult_a_tt(a,p))
        x0 = tt_compression_par(x0)
        r = tt_add(r,mult_a_tt(-a,Ap))
        r = tt_compression_par(r)
        res[j+1] = sqrt(tt_dot(r,r))
        p = tt_add(r,mult_a_tt(res[j+1]^2/res[j]^2,p))
        p = tt_compression_par(p)
        j+=1
    end
    return x0, res[1:j]
end

"""
A_k : n_k x n_k x R_{k-1} x R_k
X_k : n_k x r^X_{k-1} x r^X_k
B_k : n_k x r^B_{k-1} x r^B_k
Ql : R_{k-1} r^X_{k-1} x r^B_{k-1}
Qr : R_k r^X_k x r^B_k
"""
function init_core(A_k,dim_X,B_k,Ql,Qr)
    X_k = zeros(dim_X...)
    #B right shape
    B = reshape(B_k,:,size(B_k,3)) # n_k r^B_{k-1} x r^B_k
    B = reshape(B*(Qr'),size(B_k,1),size(B_k,2),size(A_k,4),size(X_k,3)) #n_k x r^B_{k-1} x R_k x r^X_k
    B = permutedims(B,[2,1,3,4]) #r^B_{k-1} x n_k x R_k x r^X_k
    B = Ql*reshape(B,size(B,1),:) #R_{k-1} r^X_{k-1} x n_k R_k r^X_k
    B = permutedims(reshape(B,size(A_k,3),size(X_k,2),size(A_k,1),size(A_k,4),size(X_k,3)),[3,1,4,2,5]) #n_k x R_{k-1} x R_k x r^X_{k-1} x r^X_k
    B = reshape(B,size(A_k,1)*size(A_k,3)*size(A_k,4),dim_X[2]*dim_X[3]) #n_k R_{k-1} R_k x r^X_{k-1} r^X_k
    A = reshape(permutedims(A_k,[1,3,4,2]),size(A_k,1)*size(A_k,3)*size(A_k,4),size(A_k,2))
    ua,sa,va = svd(A)
    X = va*inv(Diagonal(sa))*ua'*B
    println(norm(A*X-B))
    return reshape(X,dim_X[1],dim_X[2],dim_X[3])
end

function rand_norm(n,m) #m>=n
    A = randn(n,m)
    for i in 1:m
        A[:,i] = A[:,i]./norm(A[:,i])
    end
    return A
end

function rand_struct_orth(r_A,r_X,r_b)
    A = zeros(r_X,r_A,r_b)
    q1 = rand_norm(r_A,r_b)
    q2 = rand_orthogonal(r_b,r_X)
    for ia in 1:r_A
        for ix in 1:r_X
            for ib in 1:r_b
                A[ix,ia,ib] = q1[ia,ib]*q2[ib,ix]
            end
        end
    end
    return reshape(A,r_A*r_X,r_b)
end

function init(A::ttoperator,b::ttvector,opt_rks)
    @assert(A.tto_dims == b.ttv_dims,DimensionMismatch)
    d = length(A.tto_dims)
    opt_rks = vcat([1],opt_rks)
    Q_list = Array{Array{Float64},1}(undef,d+1)
    ttvec = Array{Array{Float64},1}(undef,d)
    Q_list[1] = [1]
    Q_list[d+1] = [1]
    for k in 1:(d-1)
        Q_list[k+1] = rand_struct_orth(A.tto_rks[k],opt_rks[k+1],b.ttv_rks[k])
    end
    for k in 1:d
        ttvec[k] = init_core(A.tto_vec[k],[A.tto_dims[k],opt_rks[k],opt_rks[k+1]],b.ttv_vec[k],Q_list[k],Q_list[k+1])
    end
    return ttvector(ttvec,A.tto_dims,opt_rks[2:(d+1)],ones(Int64,d))
end

#automatically determines the initial tt ranks
function init_adapt(A::ttoperator,b::ttvector)
    d = length(A.tto_dims)
    opt_rks = ones(Int64,d)
    for k in 1:(d-1)
        opt_rks[k] = lcm(A.tto_rks[k],b.ttv_rks[k])/A.tto_rks[k]
    end
    println(opt_rks)
    return init(A,b,opt_rks)
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

# These functions handle I/O for result files.

# *function* __`openTempFiles`__
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `v`            `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning:
# ----------------------------------------------------------------------------
# `.temp` files for the current vessel are `open`ed in the current directory
# with `w`riting permits.
# ----------------------------------------------------------------------------
#
function openTempFiles(v :: Vessel)

	v.temp_P = open(v.temp_P_name, "w")
	v.temp_Q = open(v.temp_Q_name, "w")
	v.temp_A = open(v.temp_A_name, "w")
	v.temp_c = open(v.temp_c_name, "w")
	v.temp_u = open(v.temp_u_name, "w")

end

# *function* __`openTempFiles`__
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `vessels`      `::Array{Vessel, 1}` collection of vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning:
# ----------------------------------------------------------------------------
# `.temp` files are `open`ed in the current directory with `w`riting permits
# for all the vessels in the array provided.
# ----------------------------------------------------------------------------
#
function openTempFiles(vessels :: Array{Vessel, 1})

	for v in vessels
		openTempFiles(v)
	end

end

# *function* __`closeTempFiles`__
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `v`            `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning:
# ----------------------------------------------------------------------------
# `.temp` files for the current vessel are `close`d in the current directory.
# ----------------------------------------------------------------------------
#
function closeTempFiles(v :: Vessel)

	close(v.temp_P)
	close(v.temp_Q)
	close(v.temp_A)
	close(v.temp_c)
	close(v.temp_u)

end

# *function* __`closeTempFiles`__
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `vessels`      `::Array{Vessel, 1}` collection of vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning:
# ----------------------------------------------------------------------------
# `.temp` files are `close`d in the current directory
# for all the vessels in the array provided.
# ----------------------------------------------------------------------------
#
function closeTempFiles(vessels :: Array{Vessel, 1})

	for v in vessels
		closeTempFiles(v)
	end

end

function openCloseLastFiles(v :: Vessel)

	v.last_P = open(v.last_P_name, "w")
	v.last_Q = open(v.last_Q_name, "w")
	v.last_A = open(v.last_A_name, "w")
	v.last_c = open(v.last_c_name, "w")
	v.last_u = open(v.last_u_name, "w")

	close(v.last_P)
	close(v.last_Q)
	close(v.last_A)
	close(v.last_c)
	close(v.last_u)

end

function openCloseLastFiles(vessels :: Array{Vessel, 1})

	for v in vessels
		openCloseLastFiles(v)
	end

end

# *function* __`saveTempData`__
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `t`            `::Float` current simulation time.
#
# `v`            `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning:
# ----------------------------------------------------------------------------
# `.temp` files are written with informations contained in the vessel data
# structure. For each quantity, 5 nodes are saved: the first, the middle, the
# last, and two intermediate nodes in between those three.
# ----------------------------------------------------------------------------
#
# function saveTempData(t :: Float64, v :: Vessel)

# 	write(v.temp_P, t, " ", v.P[1], " ", v.P[v.node2], " ", v.P[v.node3], " ", v.P[v.node4], " ", v.P[end], "\n")
# 	write(v.temp_Q, t, " ", v.Q[1], " ", v.Q[v.node2], " ", v.Q[v.node3], " ", v.Q[v.node4], " ", v.Q[end], "\n")
# 	write(v.temp_A, t, " ", v.A[1], " ", v.A[v.node2], " ", v.A[v.node3], " ", v.A[v.node4], " ", v.A[end], "\n")
# 	write(v.temp_c, t, " ", v.c[1], " ", v.c[v.node2], " ", v.c[v.node3], " ", v.c[v.node4], " ", v.c[end], "\n")
# 	write(v.temp_u, t, " ", v.u[1], " ", v.u[v.node2], " ", v.u[v.node3], " ", v.u[v.node4], " ", v.u[end], "\n")

# 	# write(v.temp_Q, join((t, v.Q[1], v.Q[v.node2], v.Q[v.node3],
# 	# 												 v.Q[v.node4], v.Q[end]), " "), "\n")

# 	# write(v.temp_A, join((t, v.A[1], v.A[v.node2], v.A[v.node3],
# 	# 												 v.A[v.node4], v.A[end]), " "), "\n")

# 	# write(v.temp_c, join((t, v.c[1], v.c[v.node2], v.c[v.node3],
# 	# 												 v.c[v.node4], v.c[end]), " "), "\n")

# 	# write(v.temp_u, join((t, v.u[1], v.u[v.node2], v.u[v.node3],
# 	# 												 v.u[v.node4], v.u[end]), " "), "\n")

# end

function saveTempData(t :: Float64, v :: Vessel)

	println(v.temp_P, t, " ", v.P[1], " ", v.P[v.node2], " ", v.P[v.node3], " ", v.P[v.node4], " ", v.P[end])
	println(v.temp_A, t, " ", v.A[1], " ", v.A[v.node2], " ", v.A[v.node3], " ", v.A[v.node4], " ", v.A[end])
	println(v.temp_Q, t, " ", v.Q[1], " ", v.Q[v.node2], " ", v.Q[v.node3], " ", v.Q[v.node4], " ", v.Q[end])
	println(v.temp_c, t, " ", v.c[1], " ", v.c[v.node2], " ", v.c[v.node3], " ", v.c[v.node4], " ", v.c[end])
	println(v.temp_u, t, " ", v.u[1], " ", v.u[v.node2], " ", v.u[v.node3], " ", v.u[v.node4], " ", v.u[end])

end

# *function* __`saveTempData`__
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `t`            `::Float` current simulation time.
#
# `vessels`      `::Array{Vessel, 1}` array containing vessel data structures.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning:
# ----------------------------------------------------------------------------
# `.temp` files are written with informations contained in the vessel data
# structure for all the vessels provided in the array `vessels`.
# ----------------------------------------------------------------------------
#
function saveTempData(t :: Float64, vessels :: Array{Vessel, 1})

	for v in vessels
		saveTempData(t, v)
	end

end

# *function* __`transferTempToOut`__
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `v`            `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning:
# ----------------------------------------------------------------------------
# `.temp` files are concatenated to `.out` files by `appender.sh` script.
# ----------------------------------------------------------------------------
#
function transferTempToOut(v :: Vessel)

	tempP = v.temp_P_name
	tempQ = v.temp_Q_name
	tempA = v.temp_A_name
	tempc = v.temp_c_name
	tempu = v.temp_u_name
	temps = [tempP, tempQ, tempA, tempc, tempu]

	outP  = v.out_P_name
	outQ  = v.out_Q_name
	outA  = v.out_A_name
	outc  = v.out_c_name
	outu  = v.out_u_name
	outs = [outP, outQ, outA, outc, outu]

	for (a, b) in zip(outs, temps)
		out_file = open(a, "a")
		temp_file = open(b, "r")
		write(out_file, temp_file)
		close(out_file)
		close(temp_file)
	end
	#
	# run(`sh appender.sh $tempP $outP`)
	# run(`sh appender.sh $tempQ $outQ`)
	# run(`sh appender.sh $tempA $outA`)
	# run(`sh appender.sh $tempc $outc`)
	# run(`sh appender.sh $tempu $outu`)

end

function transferLastToOut(v :: Vessel)

	lastP = v.last_P_name
	lastQ = v.last_Q_name
	lastA = v.last_A_name
	lastc = v.last_c_name
	lastu = v.last_u_name
	lasts = [lastP, lastQ, lastA, lastc, lastu]

	outP  = v.out_P_name
	outQ  = v.out_Q_name
	outA  = v.out_A_name
	outc  = v.out_c_name
	outu  = v.out_u_name
	outs = [outP, outQ, outA, outc, outu]

	for (a, b) in zip(outs, lasts)
		out_file = open(a, "a")
		last_file = open(b, "r")
		write(out_file, last_file)
		close(out_file)
		close(last_file)
	end

	# run(`sh appender.sh $lastP $outP`)
	# run(`sh appender.sh $lastQ $outQ`)
	# run(`sh appender.sh $lastA $outA`)
	# run(`sh appender.sh $lastc $outc`)
	# run(`sh appender.sh $lastu $outu`)

end

# *function* __`transferTempToOut`__
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `vessels`      `::Array{Vessel, 1}` array containing vessel data structures.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning:
# ----------------------------------------------------------------------------
# `.temp` files are concatenated to `.out` files by `appender.sh` script for
# each vessel in `vessels`.
# ----------------------------------------------------------------------------
#
function transferTempToOut(vessels :: Array{Vessel, 1})

	for v in vessels
		transferTempToOut(v)
	end

end

function transferLastToOut(vessels :: Array{Vessel, 1})

	for v in vessels
		transferLastToOut(v)
	end

end

function transferTempToLast(v :: Vessel)

	tempP = v.temp_P_name
	tempQ = v.temp_Q_name
	tempA = v.temp_A_name
	tempc = v.temp_c_name
	tempu = v.temp_u_name
	temps = [tempP, tempQ, tempA, tempc, tempu]

	lastP  = v.last_P_name
	lastQ  = v.last_Q_name
	lastA  = v.last_A_name
	lastc  = v.last_c_name
	lastu  = v.last_u_name
	lasts = [lastP, lastQ, lastA, lastc, lastu]

	# run(`sh appender.sh $tempP $lastP`)
	# run(`sh appender.sh $tempQ $lastQ`)
	# run(`sh appender.sh $tempA $lastA`)
	# run(`sh appender.sh $tempc $lastc`)
	# run(`sh appender.sh $tempu $lastu`)

	for (a, b) in zip(temps, lasts)
		temp_file = open(a, "r")
		last_file = open(b, "w")
		write(last_file, temp_file)
		close(temp_file)
		close(last_file)
	end

end

function transferTempToLast(vessels :: Array{Vessel, 1})

	for v in vessels
		transferTempToLast(v)
	end

end

function cleanTemps(vessels :: Array{Vessel, 1})

	for v in vessels
		cleanTemps(v)
	end

end

function cleanTemps(v :: Vessel)

	tempP = v.temp_P_name
	tempQ = v.temp_Q_name
	tempA = v.temp_A_name
	tempc = v.temp_c_name
	tempu = v.temp_u_name
	temps = [tempP, tempQ, tempA, tempc, tempu]

	for temp in temps
		rm(temp)
	end

end

function cleanOuts(vessels :: Array{Vessel, 1})

	for v in vessels
		cleanOuts(v)
	end

end

function cleanOuts(v :: Vessel)

	outP = v.out_P_name
	outQ = v.out_Q_name
	outA = v.out_A_name
	outc = v.out_c_name
	outu = v.out_u_name
	outs = [outP, outQ, outA, outc, outu]

	for out in outs
		rm(out)
	end

end

ext: jl


text:
using Convex

using Convex: DotMultiplyAtom
using SDPA_GMP
using Random
using Test

import LinearAlgebra.eigen
import LinearAlgebra.I
import LinearAlgebra.opnorm
import Random.shuffle
import Statistics.mean

TOL = 1e-3
eye(n) = Matrix(1.0I, n, n)

# Seed random number stream to improve test reliability
Random.seed!(2)

solvers = [SDPAGMPoptimizer(BigFloat)]

@info "Starting tests with mode PARAMETER_DEFAULT"

@testset "Convex Mode=PARAMETER_DEFAULT" begin
    include(joinpath("test", "test_utilities.jl"))
    include(joinpath("test", "test_const.jl"))
    include(joinpath("test", "test_affine.jl"))
    include(joinpath("test", "test_lp.jl"))
    # include(joinpath("test", "test_socp.jl"))
    include(joinpath("test", "test_sdp.jl"))
    # include(joinpath("test", "test_exp.jl"))
    # include(joinpath("test", "test_sdp_and_exp.jl"))
    # include(joinpath("test", "test_mip.jl"))
end


text:
"""
This function creates parallel workers for making price predictions.
"""
function create_parallel_workers(case::CaseDefinition, hpc::Bool)
    data_dir = get_data_dir(case)
    dir_name = joinpath(data_dir, "investors")
    investor_names = readdir(dir_name)

    investor_scenarios = AxisArrays.AxisArray(zeros(length(investor_names)), investor_names)

    for investor in investor_names
        investor_dir = joinpath(dir_name, investor)
        forecast_type = get_forecast_type(case)

        if forecast_type == "perfect"
            num_scenarios = 1
        elseif forecast_type == "imperfect"
            if get_uncertainty(case)
                file_name = joinpath(investor_dir, "markets_data", "scenario_multiplier_data.csv")
                @assert isfile(file_name)
                scenario_df = read_data(file_name)
                num_scenarios = DataFrames.nrow(scenario_df)
            else
                num_scenarios = 1
            end
        end

        investor_scenarios[investor] = num_scenarios
    end

    parallelize_investors = get_parallel_investors(case)
    parallelize_scenarios = get_parallel_scenarios(case)

    if parallelize_investors
        if parallelize_scenarios
            num_workers_required = sum(investor_scenarios)
        else
            num_workers_required = length(investor_scenarios)
        end
    else
        if parallelize_scenarios
            num_workers_required = maximum(investor_scenarios)
        else
            num_workers_required = 0
        end
    end

    if num_workers_required > 0
        if hpc
          nodes = split(ENV["SLURM_NODELIST"], ",")
          num_procs = min(Int(ceil(num_workers_required / length(nodes))), 3)

          node_pairs = [(n, num_procs) for n in  nodes]
          Distributed.addprocs(node_pairs)
        else
          Distributed.addprocs(min(Int(num_workers_required), 3), lazy=false)
        end
    end

    return
end

"""
This function runs price prediction if investors are parallelized but scenarios are sequential.
"""
function parallelize_only_investors(investor::Investor,
                                    sys_data_dir::String,
                                    expected_portfolio::Vector{<: Project{<: BuildPhase}},
                                    zones::Vector{String},
                                    lines::Vector{ZonalLine},
                                    peak_load::Float64,
                                    average_capital_cost_multiplier::Float64,
                                    iteration_year::Int64,
                                    yearly_horizon::Int64,
                                    solver::JuMP.MOI.OptimizerWithAttributes,
                                    sys_results_dir::String,
                                    investor_name::String)

    investor_name,
    investor_dir,
    market_names,
    rep_hour_weight,
    scenarios = gather_prediction_parameters(investor, sys_data_dir, iteration_year)

    for scenario in scenarios
        create_expected_marketdata(investor_dir,
                                market_names,
                                expected_portfolio,
                                zones,
                                lines,
                                peak_load,
                                rep_hour_weight,
                                average_capital_cost_multiplier,
                                scenario,
                                iteration_year,
                                yearly_horizon,
                                solver,
                                sys_results_dir,
                                investor_name)
    end

    return
end


text:
# # An example of DataLoader using image data

# In this example, we show how to load image data in Flux DataLoader and process it in mini-batches.
# We use the [DataLoader](https://fluxml.ai/Flux.jl/stable/data/dataloader/#Flux.Data.DataLoader) type
# to handle iteration over mini-batches of data.
# Moreover, we load the [MNIST dataset](https://juliaml.github.io/MLDatasets.jl/stable/datasets/MNIST/)
# using the [MLDatasets](https://juliaml.github.io/MLDatasets.jl/stable/) package.

# Load the packages we need:

using MLDatasets: MNIST
using Flux.Data: DataLoader
using Flux: onehotbatch
using Parameters: @with_kw

# We set a default value for the size of the mini-batches:

@with_kw mutable struct Args
    minibath_size::Int = 128  ## Size of mini-batch
end

# ## Data

# We create the function `get_data` to get, preprare and load the data onto a DataLoader object.

function get_data(args)

    ## Load the MNIST train and test data from MLDatasets
    train_x, train_y = MNIST(:train)[:]
    test_x, test_y = MNIST(:test)[:]

    ## Reshape data to 28x28x1 multi-dimensional array
    train_x = reshape(train_x, 28, 28, 1, :)
    test_x = reshape(test_x, 28, 28, 1, :)

    ## Labels must be encoded as a vector with the same dimension
    ## as the number of categories (unique handwritten digits) in the data set
    train_y, test_y = onehotbatch(train_y, 0:9), onehotbatch(test_y, 0:9)

    ## Now, we load the train and test images and their labels onto a DataLoader object
    data_loader_train = DataLoader(train_x, train_y, batchsize=args.minibath_size, shuffle=true)
    data_loader_test = DataLoader(train_x, train_y, batchsize=args.minibath_size, shuffle=true)

    return data_loader_train, data_loader_test

end

# This function performs the following tasks:
# * Loads the MNIST train and test images as Float32 as well as their labels. The dataset `train_x` is a 28×28×60000 multi-dimensional array.
#   It contains 60000 elements and each one of it contains a 28x28 array. Each array represents a 28x28 image (in grayscale) of a handwritten digit.
#   Moreover, each element of the 28x28 arrays is a pixel that represents the amount of light that it contains. On the other hand, `test_y` is a 60000 element vector and each element of this vector represents the label or actual value (0 to 9) of a handwritten digit.
# * Reshapes the train and test data to a 28x28x1 multi-dimensional array.
# * One-hot encodes the train and test labels. It creates a batch of one-hot vectors so we can pass the labels of the data as arguments for the loss function.
# * Creates two DataLoader objects that handle data mini-batches of the size defined above.

# ## Iterate over data

# Now, we can iterate over the train data during the training routine we want to define.

function train(; kws...)
    args = Args(; kws...)

    @info("Loading data...")
    data_loader_train, data_loader_test = get_data(args)

    ## Iterating over train data
    for (x, y) in data_loader_train
        @assert size(x) == (28, 28, 1, 128) || size(x) == (28, 28, 1, 96)
        @assert size(y) == (10, 128) || size(y) == (10, 96)
     end
end

# ## Run the example

# We call the `train` function:

cd(@__DIR__)
train()


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_sastre_m4_opt_rho0_69(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sastre_m4_opt_rho0_69!(A_copy)
end

@inline function exp_sastre_m4_opt_rho0_69!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=7
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 T2k7
    # Computing Bb2 = x*I+x*A
    coeff1=-0.26484760452432327
    coeff2=0.9903047529247974
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=-0.26484760452432327
    coeff2=0.9903047529247974
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.3183560628971706
    coeff2=0.9451759356891595
    coeff3=0.08341693362546203
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.041562405017987775
    coeff2=-0.30378066944691157
    coeff3=0.949941947908575
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=0.0002727973329217783
    coeff2=0.05658254301026135
    coeff3=-0.31616784713314605
    coeff4=0.9736101476608433
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.0045397184361004065
    coeff2=0.0026543192111867
    coeff3=0.0004235929351940947
    coeff4=4.468866390298052e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.10536695630497932
    coeff2=1.3909408033494417
    coeff3=0.3494801360488722
    coeff4=0.029870845474797018
    coeff5=0.4576696798781078
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.14340598505310592
    coeff2=0.18419588627720188
    coeff3=0.06331721212120305
    coeff4=0.016928981566343947
    coeff5=0.6928011352478407
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing T2k7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.9584219043602061
    coeff2=1.0556647430407748
    coeff3=0.4759594481148073
    coeff4=0.11560289591544066
    coeff5=5.152589478327363
    coeff6=0.7070319120970353
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    return memslots1 # Returning T2k7
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using JuMP, Gurobi, OneClassSampling, SVDD, OneClassActiveLearning
include(joinpath(@__DIR__, "config.jl"))

experiment_name="dami_large_outperc"

data_dirs = ["ALOI", "Annthyroid", "InternetAds", "KDDCup99", "PageBlocks", "PenDigits", "SpamBase", "Waveform", "Wilt"]

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)
quality_metrics = Dict(:mcc => matthews_corr, :kappa => cohens_kappa, :f1 => f1_score)

threshold_strategies = [
    OutlierPercentageThresholdStrategy(0.05)
    OutlierPercentageThresholdStrategy(0.25)
]

model_init_strategy = Dict{Symbol, Any}(
                            :svdd => Dict{Symbol, Any}(
                                :ROT_CTax => SVDD.SimpleCombinedStrategy(SVDD.RuleOfThumbScott(), SVDD.BoundedTaxErrorEstimate(0.05, 0.02, 0.98)),
                                :ROT_C1 => SVDD.SimpleCombinedStrategy(SVDD.RuleOfThumbScott(), SVDD.FixedCStrategy(1)),
                            ),
                            :sample_qe => Dict{Symbol, Any}(
                                :threshold_strategies => threshold_strategies
                        ))

KDE_GAMMA = :scott

competitors = [
    BPS(nothing, 0.05),
    DAEDS(30, 0.1, 0.3),
    DBSRSVDD(7, 0.5),
    HSR(20, 0.01),
    IESRSVDD(KDE_GAMMA, 0.5),
    FBPE(360),
    KFNCBD(100, 0.2),
    NDPSR(20, 10)
]

density_strategies = []
for t in threshold_strategies
    push!(density_strategies, RAPID(t, KDE_GAMMA))
end

sampling_strategies = Vector{Dict{Symbol, Any}}()

for c in competitors
    push!(sampling_strategies, Dict{Symbol, Any}(
        :method => c,
        :seed => 0
    ))
    for t in threshold_strategies
        push!(sampling_strategies, Dict{Symbol, Any}(
            :method => PreFilteringWrapper(c, KDE_GAMMA, t),
            :seed => 0
        ))
    end
end

for d in density_strategies
    push!(sampling_strategies, Dict{Symbol, Any}(
        :method => d,
        :seed => 0
    ))
end

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
function depwarn_ex(msg, name)
    return quote
        if VERSION >= v"0.6.0"
            Base.depwarn($msg, Symbol($name))
        end
    end
end

function primarytype(@nospecialize(t))
    tn = t.name
    if isdefined(tn, :primary)
        return tn.primary
    else
        return tn.wrapper
    end
end

export @functorize
macro functorize(f)
    code = f === :scalarmax          ? :(Base.scalarmax) :
           f === :scalarmin          ? :(Base.scalarmin) :
           f === :centralizedabs2fun ? :(primarytype(typeof(Base.centralizedabs2fun(0)))) :
           f
    warning = depwarn_ex("@functorize is deprecated as functor objects are no longer supported in julia", "@functorize")
    return quote
        $warning
        $code
    end
end

@static if VERSION >= v"0.6.0"
    Base.@deprecate_binding KERNEL Sys.KERNEL
    Base.@deprecate_binding UTF8String Core.String
    Base.@deprecate_binding ASCIIString Core.String
    Base.@deprecate_binding unsafe_convert Base.unsafe_convert
    Base.@deprecate_binding remote_do Distributed.remote_do
    Base.@deprecate_binding Filesystem Base.Filesystem
    Base.@deprecate_binding AsyncCondition Base.AsyncCondition
    Base.@deprecate_binding promote_eltype_op Base.promote_eltype_op
    @eval Base.@deprecate_binding $(Symbol("@irrational")) Base.$(Symbol("@irrational"))
    @eval Base.@deprecate_binding $(Symbol("@blasfunc")) Compat.LinearAlgebra.BLAS.$(Symbol("@blasfunc"))
else
    const KERNEL = Sys.KERNEL
    const UTF8String = Core.String
    const ASCIIString = Core.String
    import Base.unsafe_convert
    import Base.remote_do
    import Base.Filesystem
    import Base.AsyncCondition
    import Base.promote_eltype_op
    import Base.@irrational
    import Base.LinAlg.BLAS.@blasfunc
end

if VERSION < v"0.7.0-DEV.2915"
    const textwidth = Compat.Unicode.textwidth
    export textwidth
end


text:
module APE
#header
import SparseArrays, LinearAlgebra, ProgressMeter
using ..Meshutils, ..NLEVP
include("Meshutils_exports.jl")
include("NLEVP_exports.jl")
include("./FEM/FEM.jl")
export compute_potflow_field, discretize

function discretize(mesh::Mesh,dscrp,U;output=true)
    L=LinearOperatorFamily(["s","λ"],complex([0.,Inf]))
    N_points=size(mesh.points)[2]
    N_lines=length(mesh.lines)
    dim=N_points+3*(N_points+N_lines)
    #gas properties (air)
    P=101325 #ambient pressure (one atmosphere)
    ρ=1.225 #density
    γ=1.4 #ratio of specific heats

    triangles,tetrahedra=aggregate_elements(mesh,:quad)

    #initialize progress bar
    if output
        n_task=4*length(mesh.tetrahedra)
        for dom in keys(dscrp)
            n_task+=length(mesh.domains[dom]["simplices"])
        end

        prog = ProgressMeter.Progress(n_task,desc="Discretize... ", dt=.5,
             barglyphs=ProgressMeter.BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
             barlen=20)

    end

    #interior (term I + III)
    MM=ComplexF64[]
    II=Int64[]
    JJ=Int64[]
    #for tet in mesh.tetrahedra
    for tet in tetrahedra
        J=CooTrafo(mesh.points[:,tet[1:4]])
        #velocity components term I
        mm=ρ*s43v2u2(J) #TODO: space-dependent ρ
        ii, jj =create_indices(tet)
        iip, jjp =create_indices(tet[1:4],tet[1:4])
        #u
        for dd=1:3
            append!(MM,mm[:])
            append!(II,ii[:].+N_points.+(dd-1)*(N_points+N_lines))
            append!(JJ,jj[:].+N_points.+(dd-1)*(N_points+N_lines))
        end

        #pressure component Term III
        mm=s43v1u1(J)
        append!(MM,mm[:])
        append!(II,iip[:])
        append!(JJ,jjp[:])
        if output
            ProgressMeter.next!(prog)
        end
    end

    M=SparseArrays.sparse(II,JJ,MM,dim,dim)
    push!(L,Term(M,(pow1,),((:s,),),"s","M"))

    #Boundary matrix
    for (dom,val) in dscrp
        if dom=="Inlet"
            Ysym=:Y_in
        elseif dom =="Outlet"
            Ysym=:Y_out
        end
        L.params[Ysym]=-sqrt(γ*P/ρ)/(val/mesh.domains[dom]["size"])

        MM=ComplexF64[]
        II=Int64[]
        JJ=Int64[]
        smplx_list=mesh.domains[dom]["simplices"]
        for tri in mesh.triangles[smplx_list]
            J=CooTrafo(mesh.points[:,tri])
            ii, jj =create_indices(tri)
            mm=sqrt(γ*P/ρ)*s33v1u1(J)
            append!(MM,mm[:])
            append!(II,ii[:])
            append!(JJ,jj[:])
            if output
                ProgressMeter.next!(prog)
            end
        end
        M=SparseArrays.sparse(II,JJ,MM,dim,dim)
        push!(L,Term(M,(pow1,),((Ysym,),),string(Ysym),"B"))
    end

    #Term II and IV
    MM=ComplexF64[]
    II=Int64[]
    JJ=Int64[]
    for tet in tetrahedra
        J=CooTrafo(mesh.points[:,tet[1:4]])
        ii, jj = create_indices(tet[1:end],tet[1:4])
        iip,jjp= create_indices(tet[1:4],tet[1:end])
        for dd=1:3
            #TERM II
            #u†p
            mm=s43v2du1(J,dd)
            append!(MM,mm[:])
            append!(II,ii[:].+N_points.+(dd-1)*(N_points+N_lines))
            append!(JJ,jj[:])

            #termIV
            mm=-γ*P*s43dv1u2(J,dd) #TODO: space-dependent P
            append!(MM,mm[:])
            append!(II,iip[:])
            append!(JJ,jjp[:].+N_points.+(dd-1)*(N_points+N_lines))
        end
        if output
            ProgressMeter.next!(prog)
        end
    end
    M=SparseArrays.sparse(II,JJ,MM,dim,dim)
    push!(L,Term(M,(),(),"","K"))

    #term V and VI (mean flow)
    MM=ComplexF64[]
    II=Int64[]
    JJ=Int64[]
    for tet in tetrahedra
        J=CooTrafo(mesh.points[:,tet[1:4]])
        ii, jj =create_indices(tet)
        iip, jjp =create_indices(tet[1:4])
        #term V
        for dd=1:3
            for ee=1:3
                u=U[ee,tet[1:4]]
                #mtx= #TODO: speed up by correct ordering
                mm=ρ*(s43diffc1(J,u,dd)*s43v2u2(J)+s43v2du2c1(J,u,dd))
                append!(MM,mm[:])
                append!(II,ii[:].+N_points.+(dd-1)*(N_points+N_lines))
                append!(JJ,jj[:].+N_points.+(ee-1)*(N_points+N_lines))
            end
            #term VI
            u=U[dd,tet[1:4]]
            mm=s43v1du1c1(J,u,dd)
            append!(MM,mm[:])
            append!(II,iip[:])
            append!(JJ,jjp[:])
        end
        if output
            ProgressMeter.next!(prog)
        end
    end
    L.params[:v]=1.0
    M=SparseArrays.sparse(II,JJ,MM,dim,dim)
    push!(L,Term(M,(pow1,),((:v,),),"v","U"))

    #grid mass matrix
    MM=ComplexF64[]
    II=Int64[]
    JJ=Int64[]
    for tet in tetrahedra
        J=CooTrafo(mesh.points[:,tet[1:4]])
        mm=ρ*s43v2u2(J)
        ii, jj =create_indices(tet)
        iip, jjp =create_indices(tet[1:4])

        for dd=1:3
            #u
            append!(MM,mm[:])
            append!(II,ii[:].+N_points.+(dd-1)*(N_points+N_lines))
            append!(JJ,jj[:].+N_points.+(dd-1)*(N_points+N_lines))
        end
        #pressure component Term III
        mm=s43v1u1(J)
        append!(MM,mm[:])
        append!(II,iip[:])
        append!(JJ,jjp[:])
        if output
            ProgressMeter.next!(prog)
        end
    end
    M=SparseArrays.sparse(II,JJ,MM,dim,dim)
    push!(L,Term(-M,(pow1,),((:λ,),),"-λ","__aux__"))
    return L#II,JJ,MM,M
end

"""
    U=compute_potflow_field(mesh::Mesh,dscrp;order=:lin,output=true)

Compute potential flow field `U` on mesh `mesh` using the boundary conditions
specified in `dscrp`.

# Arguments
- `mesh::Mesh`: discretization of the computational domain
- `dscrp::Dict`: Dictionairy mapping domain names to volume flow values. Positive and negative volume flow correspond to inflow and outflow, respectively.
- `order::Symbol=:lin`: (optional) toggle whether computed velocity field is constant on a tetrahedron (`:const`) or linear interpolated from the vertices (`:lin`).
- `output::Bool=false`: (optional) toggle showing progress bar.

# Returns
- `U::Array`: 3×`N`-Array. Each column corresponds to a velocity vector. Depending on whether "order==:const" or "order==:lin" the number of columns will be `N==size(mesh.points,2)` or `N==length(mesh.tetrahedra)`, respectively.

# Notes
The flow is computed by solving the Poisson equation. The specified volume flows
must sum up to 0, otherwise the continuum equation is violated.
"""
function compute_potflow_field(mesh::Mesh,dscrp;order=:lin,output=false)
    #the potential is a stem function of the velocity
    #the element order for the potential must therefore be one order higher
    if order==:const
        triangles, tetrahedra, dim = aggregate_elements(mesh,:lin)
    elseif order==:lin
        triangles, tetrahedra, dim = aggregate_elements(mesh,:herm)
    else
        println("Error: order $order not supported for potential flow!")
        return nothing #force crash
    end
    function s43nvnu(J)
        if order==:const
            return s43nv1nu1(J)
        #elseif order==:quad
        #    return s43nv2nu2(J)
        elseif order==:lin
            return s43nvhnuh(J)
        else
            return nothing #force crash
        end
    end
    function s33v(J)
        if order==:const
            return s33v1(J)
        #elseif order==:quad
        #    return s33v2(J)
    elseif order==:lin
            return s33vh(J)
        else
            return nothing #force crash
        end
    end

    MM=Float64[]
    II=Int64[]
    JJ=Int64[]

    if output
        n_task=length(mesh.tetrahedra)
        for dom in keys(dscrp)
            n_task+=length(mesh.domains[dom]["simplices"])
        end

        p = ProgressMeter.Progress(n_task,desc="Discretize... ", dt=.5,
             barglyphs=ProgressMeter.BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
             barlen=20)

    end

    #discretize stiffness matrix (Laplacian)
    for smplx in tetrahedra
        J=CooTrafo(mesh.points[:,smplx[1:4]])
        ii, jj =create_indices(smplx)
        mm=s43nvnu(J)
        append!(MM,mm[:])
        append!(II,ii[:])
        append!(JJ,jj[:])
        if output
            ProgressMeter.next!(p)
        end
    end
    M=SparseArrays.sparse(II,JJ,MM,dim,dim)

    #discretize source terms from B.C.
    v=zeros(dim)
    for (dom,val) in dscrp
        compute_size!(mesh,dom)
        a=val/mesh.domains[dom]["size"]
        simplices=mesh.domains[dom]["simplices"]
        MM=Float64[]
        II=Int64[]
        for smplx in triangles[simplices]
            J=CooTrafo(mesh.points[:,smplx[1:3]])
            mm=-s33v(J)*a
            v[smplx[:]]+=mm[:]
            if output
                ProgressMeter.next!(p)
            end
        end
    end

    #solve for potential
    phi=M\v
    #compute velocity from grad phi
    if order==:const
        U=Array{Float64}(undef,3,length(mesh.tetrahedra))
        D=[1 0 0;
            0 1 0;
            0 0 1;
            -1 -1 -1]
        for (idx,tet) in enumerate(mesh.tetrahedra)
            J=CooTrafo(mesh.points[:,tet])
            U[:,idx]=transpose(phi[tet])*D*J.inv
        end
    elseif order==:lin
        #potential field was calculated with hermitian elements
        #both the potential and its derivative are part of the solution vector
        N_points=size(mesh.points,2)
        U=Array{Float64}(undef,3,N_points)
        U[1,:]=phi[1*N_points+1:2*N_points]
        U[2,:]=phi[2*N_points+1:3*N_points]
        U[3,:]=phi[3*N_points+1:4*N_points]
    end
    return  U
end

end


text:
mtg = read_mtg("files/simple_plant.mtg", MutableNamedTuple, NodeMTG);
classes = get_classes(mtg)
features = get_features(mtg)

@testset "test classes" begin
    @test typeof(classes) == DataFrame
    @test size(classes) == (5, 5)
    @test String.(classes.SYMBOL) == ["\$", "Individual", "Axis", "Internode", "Leaf"]
    @test classes.SCALE == [0, 1, 2, 3, 3]
    @test classes.DECOMPOSITION == ["FREE" for i = 1:5]
    @test classes.INDEXATION == ["FREE" for i = 1:5]
    @test classes.DEFINITION == ["IMPLICIT" for i = 1:5]
end

@testset "test description" begin
    @test typeof(mtg[:description]) == DataFrame
    @test size(mtg[:description]) == (2, 4)
    @test mtg[:description].LEFT == ["Internode", "Internode"]
    @test mtg[:description].RELTYPE == ["+", "<"]
    @test mtg[:description].MAX == ["?", "?"]
end

@testset "test features" begin
    @test typeof(features) == DataFrame
    @test size(features) == (3, 2)
    @test features.NAME == [:Length, :Width, :XEuler]
    @test features.TYPE == ["REAL", "REAL", "REAL"]
end

@testset "test mtg content" begin
    @test length(mtg) == 7
    @test typeof(mtg) == Node{NodeMTG,MutableNamedTuple}
    @test mtg.name == "node_1"
    @test mtg.attributes[:scales] == [0, 1, 2, 3, 3]
    @test mtg.attributes[:symbols] == ["\$", "Individual", "Axis", "Internode", "Leaf"]
    @test mtg.MTG == NodeMTG("/", "\$", 0, 0)
    @test typeof(mtg.children) == Dict{Int,Node}
    @test typeof(mtg[1]) == Node{NodeMTG,MutableNamedTuple}
    @test mtg[1].name == "node_2"
    @test mtg[1].parent === mtg
end

@testset "test mtg mutation" begin
    @test (mtg.name = "first_node") == "first_node"
    @test (mtg.attributes[:scales] .= [0, 1, 2, 3, 4]) == [0, 1, 2, 3, 4]
    @test (mtg.MTG = MultiScaleTreeGraph.NodeMTG("<", "Leaf", 2, 0)) == MultiScaleTreeGraph.NodeMTG("<", "Leaf", 2, 0)
    @test (mtg[1].parent = nothing) === nothing
    node2 = mtg[1]
    @test (mtg.children = nothing) === nothing
    @test (mtg.children = Dict(2 => node2)) == Dict(2 => node2)
    mtg.attributes = MutableNamedTuple(a = 1)
end

@testset "test mtg with NamedTuples" begin
    mtg = read_mtg("files/simple_plant.mtg", NamedTuple, NodeMTG)

    @test length(mtg) == 7
    @test typeof(mtg) == Node{NodeMTG,NamedTuple}
    @test mtg.name == "node_1"
    @test mtg.MTG == MultiScaleTreeGraph.NodeMTG("/", "\$", 0, 0)
    @test typeof(mtg.children) == Dict{Int,Node}
    @test mtg[1].name == "node_2"
    @test mtg[1].parent === mtg
end

@testset "test mtg with Dict" begin
    mtg = read_mtg("files/simple_plant.mtg", Dict, NodeMTG)
    @test length(mtg) == 7
    @test typeof(mtg) == Node{MultiScaleTreeGraph.NodeMTG,Dict{Symbol,Any}}
    @test mtg.name == "node_1"
    @test mtg.attributes == Dict(:symbols => ["\$", "Individual", "Axis", "Internode", "Leaf"],
        :scales => [0, 1, 2, 3, 3], :description => mtg[:description])
    @test mtg.MTG == MultiScaleTreeGraph.NodeMTG("/", "\$", 0, 0)
    @test typeof(mtg.children) == Dict{Int,Node}
    @test mtg[1].name == "node_2"
    @test mtg[1].parent === mtg
end

@testset "test mtg with Dict: mutation" begin
    mtg = read_mtg("files/simple_plant.mtg", Dict, NodeMTG)
    @test (mtg.name = "first_node") == "first_node"
    @test (mtg.attributes[:scales] = [0, 1, 2, 3, 4]) == [0, 1, 2, 3, 4]
    @test (mtg.MTG = MultiScaleTreeGraph.NodeMTG("<", "Leaf", 2, 0)) == MultiScaleTreeGraph.NodeMTG("<", "Leaf", 2, 0)
    @test (mtg[1].parent = nothing) === nothing
    node2 = mtg[1]
    @test (mtg.children = nothing) === nothing
    @test (mtg.children = Dict(2 => node2)) == Dict(2 => node2)
    mtg.attributes = Dict(:a => 3, :b => "test")
end


text:
#!/usr/bin/env julia

using HarmonieNamelist, JSON

io = stdout
# println(typeof(ARGS))
for arg in ARGS
    schema2namelist(io, JSON.parsefile(arg))
end


text:
#  importall Base.Operators
import Base: +,
             -,
             *,
             /,
             ^,
             %,
             !,
             >,
             <,
             >=,
             <=,
             one,
             ones,
             zero,
             zeros,
             rand,
             randn,
             trues,
             falses,
             isnan,
             sum,
             maximum,
             minimum,
             round,
             prod,
             cumsum,
             cumprod,
             diff,
             all,
             any,
             sign,
             findall,
             findfirst,
             findlast,
             log
             #  broadcast

import Base.Broadcast.broadcast

import Statistics: mean

findall(x::TS) = findall(x.values)
findfirst(x::TS) = findfirst(x.values)
findlast(x::TS) = findlast(x.values)

ones(x::TS{V,T}) where {V<:Real,T<:TimeType} = TS{V,T}(ones(V, size(x)), x.index, x.fields)
ones(::Type{TS{V}}, n::Int) where {V<:Real} = TS{V}(ones(V, n))
ones(::Type{TS{V}}, dims::Tuple{Int,Int}) where {V<:Real} = TS{V}(ones(V, dims))
ones(::Type{TS{V}}, r::Int, c::Int) where {V<:Real} = TS{V}(ones(V, r, c))
one(x::TS{V}) where {V<:Real} = one(V)
one(::Type{TS{V}}) where {V<:Real} = one(V)

zeros(x::TS{V,T}) where {V<:Real,T<:TimeType} = TS{V,T}(zeros(size(x.values)), x.index, x.fields)
zeros(::Type{TS{V}}, n::Int) where {V<:Real} = TS{V}(zeros(V, n))
zeros(::Type{TS{V}}, r::Int, c::Int) where {V<:Real} = TS{V}(zeros(V, r, c))
zeros(::Type{TS{V}}, dims::Tuple{Int,Int}) where {V<:Real} = TS{V}(zeros(V, dims))
zero(x::TS{V}) where {V<:Real} = zero(V)
zero(::Type{TS{V}}) where {V<:Real} = zero(V)

rand(::Type{TS}, n::Int=1) = TS(rand(Float64, n))
rand(::Type{TS}, r::Int, c::Int) = TS(rand(Float64, r, c))
rand(::Type{TS}, dims::Tuple{Int,Int}) = TS(rand(Float64, dims))

randn(::Type{TS}, n::Int=1) = TS(randn(Float64, n))
randn(::Type{TS}, r::Int, c::Int) = TS(randn(Float64, r, c))
randn(::Type{TS}, dims::Tuple{Int,Int}) = TS(randn(Float64, dims))

trues(x::TS) = ts(trues(size(x)), x.index, x.fields)
falses(x::TS) = ts(falses(size(x)), x.index, x.fields)
isnan(x::TS) = ts(isnan.(x.values), x.index, x.fields)
sign(x::TS) = ts(sign.(x.values), x.index, x.fields)
log(x::TS) = ts(log.(x.values), x.index, x.fields)
log(b::Number, x::TS) = ts(log.(b, x.values), x.index, x.fields)

# Number functions
round(x::TS{V}; digits::Int=0) where {V} = ts(round.(x.values, digits=digits), x.index, x.fields)
round(::Type{R}, x::TS) where {R} = ts(round.(R, x.values), x.index, x.fields)
sum(x::TS{V}) where {V} = sum(x.values)
sum(x::TS{V}, dim::Int) where {V} = sum(x.values, dim)
sum(f::Function, x::TS{V}) where {V} = sum(f, x.values)
mean(x::TS{V}) where {V} = mean(x.values)
mean(x::TS{V}, dim::Int) where {V} = mean(x, dim)
mean(f::Function, x::TS{V}) where {V} = mean(f, x.values)
prod(x::TS{V}) where {V} = prod(x.values)
prod(x::TS{V}, dim::Int) where {V} = prod(x.values, dim)
maximum(x::TS{V}) where {V} = maximum(x.values)
maximum(x::TS{V}, dim::Int) where {V} = maximum(x.values, dim)
minimum(x::TS{V}) where {V} = minimum(x.values)
minimum(x::TS{V}, dim::Int) where {V} = minimum(x.values, dim)
cumsum(x::TS; dims::Int=1) = ts(cumsum(x.values, dims=dims), x.index, x.fields)
cummin(x::TS; dims::Int=1) = ts(cummin(x.values, dims=dims), x.index, x.fields)
cummax(x::TS; dims::Int=1) = ts(cummax(x.values, dims=dims), x.index, x.fields)
cumprod(x::TS; dims::Int=1) = ts(cumprod(x.values, dims=dims), x.index, x.fields)

nans(r::Int, c::Int) = fill(NaN, 1, 2)
nans(dims::Tuple{Int,Int}) = fill(NaN, dims)

function rowdx!(dx::AbstractArray{T,N}, x::AbstractArray{T,N}, n::Int, r::Int=size(x,1)) where {T,N}
    idx = n > 0 ? (n+1:r) : (1:r+n)
    @inbounds for i=idx
        dx[i,:] .= x[i,:] - x[i-n,:]
    end
    nothing
end

function coldx!(dx::AbstractArray{T,N}, x::AbstractArray{T,N}, n::Int, c::Int=size(x,2)) where {T,N}
    idx = n > 0 ? (n+1:c) : (1:c+n)
    @inbounds for j=idx
        dx[:,j] .= x[:,j] - x[:,j-n]
    end
    nothing
end

function diffn(x::Array{T,N}, dim::Int=1, n::Int=1) where {T<:Number,N}
    @assert dim == 1 || dim == 2 "Argument `dim` must be 1 (rows) or 2 (columns)."
    @assert abs(n) < size(x,dim) "Argument `n` out of bounds."
    if n == 0
        return x
    end
    dx = zeros(eltype(x), size(x))
    if dim == 1
        rowdx!(dx, x, n)
    else
        coldx!(dx, x, n)
    end
    return dx
end

function diff(x::TS{V,T}, n::Int=1; dim::Int=1, pad::Bool=false, padval::V=zero(V)) where {V,T}
    @assert dim == 1 || dim == 2 "Argument dim must be either 1 (rows) or 2 (columns)."
    r = size(x, 1)
    c = size(x, 2)
    dx = diffn(x.values, dim, n)
    if dim == 1
        if pad
            idx = n>0 ? (1:n) : (r+n+1:r)
            dx[idx,:] .= padval
            return ts(dx, x.index, x.fields)
        else
            idx = n > 0 ? (n+1:r) : (1:r+n)
            return ts(dx[idx,:], x.index[idx], x.fields)
        end
    else
        if pad
            idx = n > 0 ? (1:c) : (c+1+1:c)
            dx[:,idx] .= padval
            return ts(dx, x.index, x.fields[idx])
        else
            idx = n > 0 ? (n+1:c) : (1:c+n)
            return ts(dx[:,idx], x.index, x.fields[idx])
        end
    end
end

function lag(x::TS{V,T}, n::Int=1; pad::Bool=false, padval::V=zero(V)) where {V,T}
    @assert abs(n) < size(x,1) "Argument `n` out of bounds."
    if n == 0
        return x
    elseif n > 0
        out = zeros(eltype(x), size(x))
        out[n+1:end,:] .= x.values[1:end-n,:]
    elseif n < 0
        out = zeros(eltype(x), size(x))
        out[1:end+n,:] .= x.values[1-n:end,:]
    end
    r = size(x, 1)
    c = size(x, 2)
    if pad
        idx = n>0 ? (1:n) : (r+n+1:r)
        out[idx,:] .= padval
        return ts(out, x.index, x.fields)
    else
        idx = n > 0 ? (n+1:r) : (1:r+n)
        return ts(out[idx,:], x.index[idx], x.fields)
    end
end

const shift = lag

function pct_change(x::TS{V}, n::Int=1; continuous::Bool=true, pad::Bool=false, padval::V=zero(V)) where {V}
    if continuous
        return diff(log(x), n; pad=pad, padval=padval)
    else
        if pad
            return diff(x, n, pad=pad, padval=padval) / x
        else
            return diff(x, n) / x[n+1:end,:]
        end
    end
end

# Artithmetic operators
function compare_elementwise(x::TS, y::TS, f::Function)
    x_cols = 1:size(x,2)
    y_cols = (1:size(y,2)) .+ size(x,2)
    merged = [x y]
    result = f.(merged.values[:,x_cols], merged.values[:,y_cols])
    return ts(result, merged.index)
end

+(x::TS) = ts(+x.values, x.index, x.fields)
-(x::TS) = ts(-x.values, x.index, x.fields)

+(x::TS, y::TS) = compare_elementwise(x, y, +)
-(x::TS, y::TS) = compare_elementwise(x, y, -)
*(x::TS, y::TS) = compare_elementwise(x, y, *)
/(x::TS, y::TS) = compare_elementwise(x, y, /)
^(x::TS, y::TS) = compare_elementwise(x, y, ^)
%(x::TS, y::TS) = compare_elementwise(x, y, %)

+(x::TS, y::AbstractArray) = x + TS(y, x.index, x.fields)
-(x::TS, y::AbstractArray) = x - TS(y, x.index, x.fields)
*(x::TS, y::AbstractArray) = x * TS(y, x.index, x.fields)
/(x::TS, y::AbstractArray) = x / TS(y, x.index, x.fields)
%(x::TS, y::AbstractArray) = x % TS(y, x.index, x.fields)
^(x::TS, y::AbstractArray) = x ^ TS(y, x.index, x.fields)

+(y::AbstractArray, x::TS) = x + y
-(y::AbstractArray, x::TS) = x - y
*(y::AbstractArray, x::TS) = y * x
/(y::AbstractArray, x::TS) = y / x
%(y::AbstractArray, x::TS) = y % x
^(y::AbstractArray, x::TS) = y ^ x

+(x::TS, y::Number) = ts(x.values .+ y, x.index, x.fields)
-(x::TS, y::Number) = ts(x.values .- y, x.index, x.fields)
*(x::TS, y::Number) = ts(x.values .* y, x.index, x.fields)
/(x::TS, y::Number) = ts(x.values ./ y, x.index, x.fields)
%(x::TS, y::Number) = ts(x.values .% y, x.index, x.fields)
^(x::TS{B}, y::E) where {B<:Real,E<:Real} = ts(x.values .^ y, x.index, x.fields)

+(y::Number, x::TS) = x + y
-(y::Number, x::TS) = x - y
*(y::Number, x::TS) = x * y
/(y::Number, x::TS) = x / y
%(y::Number, x::TS) = x % y
^(y::Number, x::TS) = x ^ y

# Logical operators
all(x::TS) = all(x.values)
any(x::TS) = any(x.values)

!(x::TS) = TS(.!x.values, x.index, x.fields)

# compared to another ts object
==(x::TS, y::TS) = x.values == y.values && x.index == y.index && x.fields == y.fields
!=(x::TS, y::TS) = !(x == y)
>(x::TS, y::TS) = compare_elementwise(x, y, >)
<(x::TS, y::TS) = compare_elementwise(x, y, <)
>=(x::TS, y::TS) = compare_elementwise(x, y, >=)
<=(x::TS, y::TS) = compare_elementwise(x, y, <=)

#  # Vectorizized operations
#  broadcast(::typeof(+), x::TS, y::TS) = ts(x.values .+ y.values, x.index, x.fields)
#  broadcast(::typeof(-), x::TS, y::TS) = ts(x.values .- y.values, x.index, x.fields)
#  broadcast(::typeof(*), x::TS, y::TS) = ts(x.values .* y.values, x.index, x.fields)
#  broadcast(::typeof(/), x::TS, y::TS) = ts(x.values ./ y.values, x.index, x.fields)
#  broadcast(::typeof(^), x::TS, y::TS) = ts(x.values .^ y.values, x.index, x.fields)
#  broadcast(::typeof(%), x::TS, y::TS) = ts(x.values .% y.values, x.index, x.fields)
#  broadcast(::typeof(+), x::TS, y::AbstractArray) = ts(x.values .+ y, x.index, x.fields)
#  broadcast(::typeof(-), x::TS, y::AbstractArray) = ts(x.values .- y, x.index, x.fields)
#  broadcast(::typeof(*), x::TS, y::AbstractArray) = ts(x.values .* y, x.index, x.fields)
#  broadcast(::typeof(/), x::TS, y::AbstractArray) = ts(x.values ./ y, x.index, x.fields)
#  broadcast(::typeof(%), x::TS, y::AbstractArray) = ts(x.values .% y, x.index, x.fields)
#  broadcast(::typeof(^), x::TS, y::AbstractArray) = ts(x.values .^ y, x.index, x.fields)
#  broadcast(::typeof(+), y::AbstractArray, x::TS) = x .+ y
#  broadcast(::typeof(-), y::AbstractArray, x::TS) = x .- y
#  broadcast(::typeof(*), y::AbstractArray, x::TS) = ts(y .* x.values, x.index, x.fields)
#  broadcast(::typeof(/), y::AbstractArray, x::TS) = ts(y ./ x.values, x.index, x.fields)
#  broadcast(::typeof(%), y::AbstractArray, x::TS) = ts(y .% x.values, x.index, x.fields)
#  broadcast(::typeof(^), y::AbstractArray, x::TS) = ts(y .^ x.values, x.index, x.fields)
#  broadcast(::typeof(+), x::TS, y::Number) = ts(x.values .+ y, x.index, x.fields)
#  broadcast(::typeof(-), x::TS, y::Number) = ts(x.values .- y, x.index, x.fields)
#  broadcast(::typeof(*), x::TS, y::Number) = ts(x.values .* y, x.index, x.fields)
#  broadcast(::typeof(/), x::TS, y::Number) = ts(x.values ./ y, x.index, x.fields)
#  broadcast(::typeof(%), x::TS, y::Number) = ts(x.values .% y, x.index, x.fields)
#  broadcast(::typeof(^), x::TS, y::Number) = ts(x.values .^ y, x.index, x.fields)
#  broadcast(::typeof(+), y::Number, x::TS) = ts(y .+ x.values, x.index, x.fields)
#  broadcast(::typeof(-), y::Number, x::TS) = ts(y .- x.values, x.index, x.fields)
#  broadcast(::typeof(*), y::Number, x::TS) = ts(y .* x.values, x.index, x.fields)
#  broadcast(::typeof(/), y::Number, x::TS) = ts(y ./ x.values, x.index, x.fields)
#  broadcast(::typeof(%), y::Number, x::TS) = ts(y .% x.values, x.index, x.fields)
#  broadcast(::typeof(^), y::Number, x::TS) = ts(y .^ x.values, x.index, x.fields)
#  broadcast(::typeof(!), x::TS) = ts(.!x.values, x.index, x.fields)
#  broadcast(::typeof(!=), x::TS, y::TS) = compare_elementwise(x, y, !=)
#  broadcast(::typeof(==), x::TS, y::TS) = compare_elementwise(x, y, ==)
#  broadcast(::typeof(>), x::TS, y::TS) = compare_elementwise(x, y, >)
#  broadcast(::typeof(<), x::TS, y::TS) = compare_elementwise(x, y, <)
#  broadcast(::typeof(>=), x::TS, y::TS) = compare_elementwise(x, y, >=)
#  broadcast(::typeof(<=), x::TS, y::TS) = compare_elementwise(x, y, <=)


text:
using Test

import FFTW

using Oceananigans
using Oceananigans.Operators

archs = [CPU()]
@hascuda archs = [CPU(), GPU()]

float_types = [Float32, Float64]

@testset "Oceananigans" begin
    println("Testing Oceananigans...")

    @testset "Grid" begin
        println("  Testing grids...")
        include("test_grids.jl")

        @testset "Grid initialization" begin
            println("    Testing grid initialization...")
            for ft in float_types
                @test correct_grid_size(ft)
                @test correct_cell_volume(ft)
                @test faces_start_at_zero(ft)
                @test end_faces_match_grid_length(ft)
            end
        end

        @testset "Grid dimensions" begin
            println("    Testing grid dimensions...")
            L = (100, 100, 100)
            for ft in float_types
                @test isbitstype(typeof(RegularCartesianGrid(ft, (16, 16, 16), (1, 1, 1))))

                @test RegularCartesianGrid(ft, (25, 25, 25), L).dim == 3
                @test RegularCartesianGrid(ft, (5, 25, 125), L).dim == 3
                @test RegularCartesianGrid(ft, (64, 64, 64), L).dim == 3
                @test RegularCartesianGrid(ft, (32, 32,  1), L).dim == 2
                @test RegularCartesianGrid(ft, (32,  1, 32), L).dim == 2
                @test RegularCartesianGrid(ft, (1,  32, 32), L).dim == 2
                @test RegularCartesianGrid(ft, (1,  1,  64), L).dim == 1

                @test_throws ArgumentError RegularCartesianGrid(ft, (32,), L)
                @test_throws ArgumentError RegularCartesianGrid(ft, (32, 64), L)
                @test_throws ArgumentError RegularCartesianGrid(ft, (1, 1, 1), L)
                @test_throws ArgumentError RegularCartesianGrid(ft, (32, 32, 32, 16), L)
                @test_throws ArgumentError RegularCartesianGrid(ft, (32, 32, 32), (100,))
                @test_throws ArgumentError RegularCartesianGrid(ft, (32, 32, 32), (100, 100))
                @test_throws ArgumentError RegularCartesianGrid(ft, (32, 32, 32), (100, 100, 1, 1))
                @test_throws ArgumentError RegularCartesianGrid(ft, (32, 32, 32), (100, 100, -100))

                @test_throws ArgumentError RegularCartesianGrid(ft, (32, 32, 32.0), (1, 1, 1))
                @test_throws ArgumentError RegularCartesianGrid(ft, (20.1, 32, 32), (1, 1, 1))
                @test_throws ArgumentError RegularCartesianGrid(ft, (32, nothing, 32), (1, 1, 1))
                @test_throws ArgumentError RegularCartesianGrid(ft, (32, "32", 32), (1, 1, 1))
                @test_throws ArgumentError RegularCartesianGrid(ft, (32, 32, 32), (1, nothing, 1))
                @test_throws ArgumentError RegularCartesianGrid(ft, (32, 32, 32), (1, "1", 1))
            end
        end
    end

    @testset "Fields" begin
        println("  Testing fields...")
        include("test_fields.jl")

        N = (4, 6, 8)
        L = (2π, 3π, 5π)

        field_types = [CellField, FaceFieldX, FaceFieldY, FaceFieldZ, EdgeField]

        @testset "Field initialization" begin
            println("    Testing field initialization...")
            for arch in archs, ft in float_types
                grid = RegularCartesianGrid(ft, N, L)

                for field_type in field_types
                    @test correct_field_size(arch, grid, field_type)
                end
            end
        end

        int_vals = Any[0, Int8(-1), Int16(2), Int32(-3), Int64(4), Int128(-5)]
        uint_vals = Any[6, UInt8(7), UInt16(8), UInt32(9), UInt64(10), UInt128(11)]
        float_vals = Any[0.0, -0.0, 6e-34, 1.0f10]
        rational_vals = Any[1//11, -23//7]
        other_vals = Any[π]
        vals = vcat(int_vals, uint_vals, float_vals, rational_vals, other_vals)

        @testset "Setting fields" begin
            println("    Testing field setting...")

            for arch in archs, ft in float_types
                grid = RegularCartesianGrid(ft, N, L)

                for field_type in field_types, val in vals
                    @test correct_field_value_was_set(arch, grid, field_type, val)
                end
            end
        end

        # @testset "Field operations" begin
        #     for arch in archs, ft in float_types
        #         grid = RegularCartesianGrid(ft, N, L)
        #
        #         for field_type in field_types, val1 in vals, val2 in vals
        #             @test correct_field_addition(arch, grid, field_type, val1, val2)
        #         end
        #     end
        # end
    end

    @testset "Operators" begin
        println("  Testing operators...")

        @testset "2D operators" begin
            println("    Testing 2D operators...")
            Nx, Ny, Nz = 32, 16, 8
            Lx, Ly, Lz = 100, 100, 100

            grid = RegularCartesianGrid((Nx, Ny, Nz), (Lx, Ly, Lz))
            A3 = rand(Nx, Ny, Nz)

            A2yz = A3[1:1, :, :]  # A yz-slice with Nx==1.
            grid_yz = RegularCartesianGrid((1, Ny, Nz), (Lx, Ly, Lz))

            A2xz = A3[:, 1:1, :]  # An xz-slice with Ny==1.
            grid_xz = RegularCartesianGrid((Nx, 1, Nz), (Lx, Ly, Lz))

            test_indices_3d = [(4, 5, 5), (21, 11, 4), (16, 8, 4),  (30, 12, 3), (11, 3, 6), # Interior
                               (2, 10, 4), (31, 5, 6), (10, 2, 4), (17, 15, 5), (17, 10, 2), (23, 5, 7),  # Borderlands
                               (1, 5, 5), (32, 10, 3), (16, 1, 4), (16, 16, 4), (16, 8, 1), (16, 8, 8),  # Edges
                               (1, 1, 1), (32, 16, 8)] # Corners

            test_indices_2d_yz = [(1, 1, 1), (1, 1, 2), (1, 2, 1), (1, 2, 2),
                                  (1, 1, 5), (1, 5, 1), (1, 5, 5), (1, 11, 4),
                                  (1, 15, 7), (1, 15, 8), (1, 16, 7), (1, 16, 8)]

            test_indices_2d_xz = [(1, 1, 1), (1, 1, 2), (2, 1, 1), (2, 1, 2),
                                  (1, 1, 5), (5, 1, 1), (5, 1, 5), (17, 1, 4),
                                  (31, 1, 7), (31, 1, 8), (32, 1, 7), (32, 1, 8)]

            for idx in test_indices_2d_yz
                @test δx_f2c(grid_yz, A2yz, idx...) ≈ 0
                @test δx_c2f(grid_yz, A2yz, idx...) ≈ 0
                @test δy_f2c(grid_yz, A2yz, idx...) ≈ δy_f2c(grid_yz, A3, idx...)
                @test δy_c2f(grid_yz, A2yz, idx...) ≈ δy_c2f(grid_yz, A3, idx...)
                @test δz_f2c(grid_yz, A2yz, idx...) ≈ δz_f2c(grid_yz, A3, idx...)
                @test δz_c2f(grid_yz, A2yz, idx...) ≈ δz_c2f(grid_yz, A3, idx...)
            end

            for idx in test_indices_2d_xz
                @test δx_f2c(grid_xz, A2xz, idx...) ≈ δx_f2c(grid_xz, A3, idx...)
                @test δx_c2f(grid_xz, A2xz, idx...) ≈ δx_c2f(grid_xz, A3, idx...)
                @test δy_f2c(grid_xz, A2xz, idx...) ≈ 0
                @test δy_c2f(grid_xz, A2xz, idx...) ≈ 0
                @test δz_f2c(grid_xz, A2xz, idx...) ≈ δz_f2c(grid_xz, A3, idx...)
                @test δz_c2f(grid_xz, A2xz, idx...) ≈ δz_c2f(grid_xz, A3, idx...)
            end
        end
    end

    @testset "Poisson solvers" begin
        println("  Testing Poisson solvers...")
        include("test_poisson_solvers.jl")

        @testset "FFTW commutativity" begin
            println("    Testing FFTW commutativity...")

            # Testing for nice powers of two and not-so-nice prime numbers.
            for N in [4, 7, 8, 10, 48, 64, 79, 128]
                @test mixed_fft_commutes(N)
                @test mixed_ifft_commutes(N)
            end
        end

        @testset "FFTW plans" begin
            println("    Testing FFTW planning...")

            for ft in float_types
                @test fftw_planner_works(ft, 32, 32, 32, FFTW.ESTIMATE)
                @test fftw_planner_works(ft, 1,  32, 32, FFTW.ESTIMATE)
                @test fftw_planner_works(ft, 32,  1, 32, FFTW.ESTIMATE)
                @test fftw_planner_works(ft,  1,  1, 32, FFTW.ESTIMATE)
            end
        end

        @testset "Divergence-free solution [CPU]" begin
            println("    Testing divergence-free solution [CPU]...")

            for N in [5, 7, 10, 15, 20, 29, 32]
                for ft in float_types
                    @test poisson_ppn_planned_div_free_cpu(ft, N, N, N, FFTW.ESTIMATE)
                    @test poisson_ppn_planned_div_free_cpu(ft, 1, N, N, FFTW.ESTIMATE)
                    @test poisson_ppn_planned_div_free_cpu(ft, N, 1, N, FFTW.ESTIMATE)
                    @test poisson_ppn_planned_div_free_cpu(ft, 1, 1, N, FFTW.ESTIMATE)

                    # Commented because https://github.com/climate-machine/Oceananigans.jl/issues/99
                    # for planner_flag in [FFTW.ESTIMATE, FFTW.MEASURE]
                    #     @test test_3d_poisson_ppn_planned!_div_free(mm, N, N, N, planner_flag)
                    #     @test test_3d_poisson_ppn_planned!_div_free(mm, 1, N, N, planner_flag)
                    #     @test test_3d_poisson_ppn_planned!_div_free(mm, N, 1, N, planner_flag)
                    # end
                end
            end

            Ns = [5, 7, 10, 15, 20, 29, 32]
            for Nx in Ns, Ny in Ns, Nz in Ns, ft in float_types
                @test poisson_ppn_planned_div_free_cpu(ft, Nx, Ny, Nz, FFTW.ESTIMATE)
            end
        end
    end

    @testset "Model" begin
        println("  Testing model...")

        for arch in archs, ft in float_types
            model = Model(N=(4, 5, 6), L=(1, 2, 3), arch=arch, float_type=ft)

            # Just testing that a Model was constructed with no errors/crashes.
            @test true
        end
    end

    @testset "Time stepping" begin
        println("  Testing time stepping...")
        include("test_time_stepping.jl")

        for arch in archs, ft in float_types
            @test time_stepping_works(arch, ft)
        end

        @testset "Adams-Bashforth 2" begin
            for arch in archs, ft in float_types
                run_first_AB2_time_step_tests(arch, ft)
            end
        end
    end

    @testset "Boundary conditions" begin
        println("  Testing boundary conditions...")
        include("test_boundary_conditions.jl")

        Nx, Ny, Nz = 3, 4, 5 # for simple test
        funbc(args...) = π

        for fld in (:u, :v, :T, :S)
            for bctype in (Gradient, Flux)
                for bc in (0.6, rand(Nx, Ny), funbc)
                    @test test_z_boundary_condition_simple(fld, bctype, bc, Nx, Ny, Nz)
                end
            end

            @test test_diffusion_simple(fld)
            @test test_diffusion_budget(fld)
            @test test_diffusion_cosine(fld)
            @test test_flux_budget(fld)
        end
    end

    @testset "Forcing" begin
        println("  Testing forcings...")
        add_one(args...) = 1.0
        function test_forcing(fld)
            kwarg = Dict(Symbol(:F, fld)=>add_one)
            forcing = Forcing(; kwarg...)
            f = getfield(forcing, fld)
            f() == 1.0
        end

        for fld in (:u, :v, :w, :T, :S)
            @test test_forcing(fld)
        end
    end

    @testset "Output writers" begin
        println("  Testing output writers...")
        include("test_output_writers.jl")

        @testset "Checkpointing" begin
            println("    Testing checkpointing...")
            run_thermal_bubble_checkpointer_tests()
        end

        @testset "NetCDF" begin
            println("    Testing NetCDF output writer...")
            run_thermal_bubble_netcdf_tests()
        end
    end

    @testset "Golden master tests" begin
        include("test_golden_master.jl")

        @testset "Thermal bubble" begin
            for arch in archs
                run_thermal_bubble_golden_master_tests(arch)
            end
        end

        @testset "Deep convection" begin
            run_deep_convection_golden_master_tests()
        end
    end
end # Oceananigans tests


text:
using Test
using SafeTestsets

@testset "DCA" begin
    @safetestset "Size tests" begin include("dca_test.jl") end
    @safetestset "ECA ruleset tests" begin include("eca_ruleset_test.jl") end
end

@testset "TCA" begin
    @safetestset "Size tests" begin include("tca_test.jl") end
end

@testset "CCA" begin
    @safetestset "Size tests" begin include("cca_test.jl") end
end

@testset "Life-like" begin
    @safetestset "Life glider" begin include("glider_test.jl") end
    @safetestset "Life blinker" begin include("blinker_test.jl") end
end


text:
import Base: getindex, lastindex

@doc doc"""
    load_chains(filename; burnin=0, thin=1, join=false)

Load a single chain or multiple parallel chains which were written to a file by
[`sample_joint`](@ref).

Keyword arguments:

* `burnin` — Remove this many samples from the start of each chain.
* `thin` — If `thin` is an integer, thin the chain by this factor. If
  `thin == :hasmaps`, return only samples which have maps saved. If thin is a
  `Function`, filter the chain by this function (e.g. `thin=haskey(:g)` on Julia 1.5+)
* `unbatch` — If true, [unbatch](@ref) the chains if they are batched.
* `join` — If true, concatenate all the chains together.
* `skip_missing_chunks` — Skip missing chunks in the chain instead of
  terminating the chain there.

The object returned by this function is a `Chain` or `Chains` object, which
simply wraps an `Array` of `Dicts` or an `Array` of `Array` of `Dicts`,
respectively (each sample is a `Dict`). The wrapper object has some extra
indexing properties for convenience:

* It can be indexed as if it were a single multidimensional object, e.g.
  `chains[1,:,:accept]` would return the `:accept` key of all samples in the
  first chain.
* Leading colons can be dropped, i.e. `chains[:,:,:accept]` is the same as
  `chains[:accept]`.
* If some samples are missing a particular key, `missing` is returned for those
  samples insted of an error.
* The recursion goes arbitrarily deep into the objects it finds. E.g., since
  sampled parameters are stored in a `NamedTuple` like `(Aϕ=1.3,)` in the `θ`
  key of each sample `Dict`, you can do `chain[:θ,:Aϕ]` to get all `Aϕ` samples
  as a vector.

"""
function load_chains(filename; burnin=0, thin=1, join=false, unbatch=true, dropmaps=false)
    chains = jldopen(filename) do io
        ks = keys(io)
        chunk_ks = [k for k in ks if startswith(k,"chunks_")]
        for (isfirst,k) in flagfirst(sort(chunk_ks, by=k->parse(Int,k[8:end])))
            if isfirst
                chains = read(io,k)
            else
                append!.(chains, read(io,k))
            end
            if dropmaps
                map!(chain -> filter!.(((_,v),)->!(v isa Field), chain), chains, chains)
            end
        end
        chains
    end
    if thin isa Int
        chains = [chain[(1+burnin):thin:end] for chain in chains]
    elseif thin == :hasmaps
        chains = [[samp for samp in chain[(1+burnin):end] if :ϕ in keys(samp)] for chain in chains]
    elseif thin isa Function
        chains = [filter(thin,chain) for chain in chains]
    else
        error("`thin` should be an Int or :hasmaps")
    end
    chains = wrap_chains(chains)
    if unbatch
        chains = CMBLensing.unbatch(chains)
    end
    return join ? Chain(reduce(vcat, chains)) : chains
end

# multiple parallel chains
struct Chains{T} <: AbstractVector{T}
    chains :: Vector{T}
end
Chains(chains::Chains) = chains
getindex(c::Chains, k::Symbol, ks...) = _getindex(c.chains, :, :, k, ks...)
getindex(c::Chains, k,         ks...) = _getindex(c.chains, k, ks...)
size(c::Chains) = size(c.chains)
function Base.print_array(io::IO, cs::Chains; indent="  ")
    for c in cs
        print(io, "  ")
        Base.summary(io, c)
        println(io)
        Base.print_array(io, c, indent="    ")
    end
end
function lastindex(c::Chains, d)
    d==1 ? lastindex(c.chains) :
    d==2 ? lastindex(c.chains[1]) :
    error("`end` only valid in first or second dim of Chains")
end

# single chain
struct Chain{T} <: AbstractVector{T}
    chain :: Vector{T}
end
getindex(c::Chain, k::Symbol,              ks...) = _getindex(c.chain, :, k, ks...)
getindex(c::Chain, k,                      ks...) = _getindex(c.chain, k, ks...)
getindex(c::Chain, k::Union{Colon,AbstractRange}) = Chain(getindex(c.chain, k))
lastindex(c::Chain) = lastindex(c.chain)
lastindex(c::Chain, d) = d==1 ? lastindex(c.chain) : error("`end` only valid in first dim of Chain")
size(c::Chain) = size(c.chain)
function Base.print_array(io::IO, c::Chain; indent="  ")
    _,cols = displaysize(io)
    for k in keys(c[end])
        str = string("$(indent)$(k) => ", repr(c[k]; context=(:limit => true)))
        println(io, Base._truncate_at_width_or_chars(str, cols))
    end
end

# recurse
_getindex(x, k,                             k2, ks...) = _getindex(_getindex(x, k), k2, ks...)
_getindex(x, k::Union{Colon,AbstractRange}, k2, ks...) = broadcast(y -> _getindex(y, k2, ks...),  x[k])

# base cases
_getindex(x::Union{Dict,NamedTuple}, k::Symbol) = haskey(x,k) ? getindex(x, k) : missing
_getindex(x,                         k) = getindex(x, k)

wrap_chains(chains::Vector{<:Vector{<:Dict}}) = Chains(Chain.(chains))
wrap_chains(chain::Vector{<:Dict}) = Chain(chain)

# batching
@doc doc"""
    unbatch(chain::Chain)

Convert a chain of batch-length-`D` fields to `D` chains of unbatched fields.
"""
function unbatch(chain::Chain)
    D = batchsize(chain[end][:lnP])
    (D==1) && return chain
    Chains(map(1:D) do I
        Chain(map(chain) do samp
            Dict(map(collect(samp)) do (k,v)
                if v isa Union{BatchedReal,FlatField}
                    k => batchindex(v, I)
                elseif v isa NamedTuple{<:Any, <:NTuple{<:Any,<:BatchedVals}}
                    k => map(x -> (x isa BatchedVals ? batchindex(x,I) : x), v)
                elseif v isa AbstractArray{<:BatchedVals}
                    k => batchindex.(v, I)
                else
                    k => v
                end
            end)
        end)
    end)
end

@doc doc"""
    unbatch(chains::Chains)

Expand each chain in this `Chains` object by unbatching it.
"""
unbatch(chains::Chains) = Chains(mapreduce(unbatch, vcat, chains))

@doc doc"""
    mean_std_and_errors(samples; N_bootstrap=10000)

Get the mean and standard deviation of a set of correlated `samples` from a
chain where the error on the mean and standard deviation is estimated with
bootstrap resampling using the calculated "effective sample size" of the chain.
"""
function mean_std_and_errors(samples; N_bootstrap=10000)

    Neff = round(Int, length(samples) / @ondemand(PyCall.pyimport)(:emcee).autocorr.integrated_time(samples)[1])

    μ = mean(samples)
    σ = std(samples)

    SEμ = std([mean(samples[rand(1:end, Neff)]) for i=1:N_bootstrap])
    SEσ = std([ std(samples[rand(1:end, Neff)]) for i=1:N_bootstrap])

    "$(paren_errors(μ, SEμ)) ± $(paren_errors(σ, SEσ))"

end

@doc doc"""
    paren_errors(μ, σ; N_in_paren=2)

Get a string represntation of `μ ± σ` in "parenthesis" format, e.g. `1.234 ±
0.012` becomes `1.234(12)`.
"""
function paren_errors(μ, σ; N_in_paren=2)
    N = round(Int, floor(log10(1/σ))) + N_in_paren
    fmt = "%.$(N)f"
    @ondemand(Formatting.sprintf1)(fmt, μ)*"($(round(Int,σ*10^N)))"
end


text:
julia> print(j(41,3))
Prisoner killing order: 2, 5, 8, 11, 14, 17, 20, 23, 26, 29, 32, 35, 38, 0, 4, 9, 13, 18, 22, 27, 31, 36, 40, 6, 12, 19, 25, 33, 39, 7, 16, 28, 37, 10, 24, 1, 21, 3, 34, 15.
Survivor: 30

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode