text:
# Compensating variations
using CollegeStratWorker, Random, Test

csw = CollegeStratWorker;

function cv_test()
    @testset "Compensating variations" begin
        wk = csw.make_test_worker(2);
        ltIncome = 23.0;
        workStartAge = 3;
        ltuBase = csw.lifetime_utility(wk, workStartAge, ltIncome);
        @test !isinf(ltuBase)

        lb = ltIncome / 10.0;
        ub = ltIncome * 10.0;
        ltIncomeNew = csw.ltincome_from_utility(wk, workStartAge, ltuBase, lb, ub);
        @test ltIncomeNew > 0.0
        @test isapprox(ltIncome, ltIncomeNew, rtol = 0.001)

        # Solve many
        rng = MersenneTwister(43);
        n = 110;
        ltIncomeV = 1000.0 .+ rand(rng, n) .* 3000.0;
        ltUtilV = [csw.lifetime_utility(wk, workStartAge, ltIncomeV[j])
            for j = 1 : n];
        ltIncomeNewV = csw.ltincome_from_utility(wk, workStartAge, ltUtilV,
            100.0, 10_000.0);
        @test isapprox(ltIncomeV, ltIncomeNewV, rtol = 0.001)
	end
end

@testset "Compensating variations" begin
    cv_test()
end

# --------


text:
export inner

"""
    inner(x::AbstractVector{N}, A::AbstractMatrix{N}, y::AbstractVector{N}
         ) where {N}

Compute the inner product ``xᵀ A y``.

### Input

- `x` -- vector on the left
- `A` -- matrix
- `y` -- vector on the right

### Output

The (scalar) result of the multiplication.
"""
function inner(x::AbstractVector{N}, A::AbstractMatrix{N}, y::AbstractVector{N}
              ) where {N}
    dot(x, A * y)
end

ext: jl


text:
## PARTE 1
cd(@__DIR__)
using DelimitedFiles
input = readdlm("input.txt");
get_move(::Val{:forward}, x) = CartesianIndex(x, 0)
get_move(::Val{:up}, x) = CartesianIndex(0, -x)
get_move(::Val{:down}, x) = CartesianIndex(0, x)
get_move(::Val, x) = missing

##
origin = CartesianIndex(0, 0)
res = foldl(eachrow(input); init=origin) do pos, (s, Δ)
    move = get_move(Val(Symbol(s)), Δ)
    ismissing(move) && return
    pos += move
end
prod(res.I)


text:
using BinaryBuilder, Pkg

# Collection of sources required to build XGBoost
name = "XGBoost"
version = v"1.1.0"
sources = [
    GitSource("https://github.com/dmlc/xgboost.git","f5d4fddafec7ff0e908b776e0cebf554df6b89d0"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/xgboost
git submodule init
git submodule update

# Patch dmlc-core to use case-sensitive windows.h includes
(cd dmlc-core; atomic_patch -p1 "${WORKSPACE}/srcdir/patches/dmlc_windows.patch")

# For Linux, build using CMake
if [[ ${target} == *linux* ]] ||  [[ ${target} == *mingw* ]]; then
    (mkdir build; cd build; cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}")
    make -C build -j ${nproc}
# For Mac and FreeBSD, build without openmp
elif [[ ${target} == *darwin* ]] ||  [[ ${target} == *freebsd* ]]; then
    (mkdir build; cd build; cmake .. -DUSE_OPENMP=OFF -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}")
    make -C build -j ${nproc}
fi

# Install
mkdir -p ${prefix}/{bin,include,lib}
cp -ra include/xgboost ${prefix}/include/
cp -a xgboost${exeext} ${prefix}/bin/xgboost${exeext}

if [[ ${target} == *mingw* ]]; then
    cp -a lib/xgboost.dll ${prefix}/bin
else
    cp -a lib/libxgboost.${dlext} ${prefix}/lib
fi
"""

# 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())

# Disable powerpc for now
platforms = [p for p in platforms if !(arch(p) == :powerpc64le)]
# The products that we will ensure are always built
products = [
    LibraryProduct(["libxgboost", "xgboost"], :libxgboost),
    ExecutableProduct("xgboost", :xgboost)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
]

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


text:
# This file contains tests for result type

# A module which contains the code to
# perform tests on .Result module from Scats
module TestResult

using Scats.Internal.Prec         # Precision module from Scats
using Scats: Internal.Result, API # API and .Result module from Scats
using Test                        # A package to perform tests

# Print the header
println("\e[1;32mRUNNING\e[0m: result_test.jl")

# Creating instance of the API
s = API()

# Test values resetting
@testset "Check resetting" begin

    # Set non-default values
    s.Result.Δt = s.Result.q = s.Result.threshold = 1
    s.Result.t = s.Result.x = s.Result.X_FFT_ABS = [1]
    s.Result.ν = s.Result.D = s.Result.c = s.Result.cw = s.Result.Dw = [1]

    # Reset values
    s.Result.reset!()

    # Test values
    @test s.Result.Δt == 0
    @test s.Result.q == 0
    @test s.Result.threshold == 0
    @test s.Result.t == []
    @test s.Result.x == []
    @test s.Result.X_FFT_ABS == []
    @test s.Result.ν == []
    @test s.Result.D == []
    @test s.Result.c == []
    @test s.Result.cw == []
    @test s.Result.Dw == []

end

end


text:
module CUDArt

using Compat

export
    # pointer symbols
    CudaPtr, rawpointer, CUDA_NULL,
    # array symbols
    AbstractCudaArray, AbstractCudaVector, AbstractCudaMatrix,
    CudaArray, CudaVector, CudaMatrix, CudaVecOrMat,
    CudaPitchedArray, HostArray, to_host,
    # other symbols
    device, devices, device_reset, attribute, capability,
    driver_version, runtime_version,
    CuModule, CuFunction, unload,
    pitchel, pitchbytes,
    launch, device_synchronize, synchronize,
    Stream, null_stream, cudasleep,
    destroy, free, cudafinalizer

import Base: length, size, ndims, eltype, similar, pointer, stride,
    copy, convert, reinterpret, show, summary,
    copy!, get!, fill!, wait

# Prepare the CUDA runtime API bindings
include("libcudart-6.5.jl")
import .CUDArt_gen
const rt = CUDArt_gen

# To load PTX code, we also need access to the driver API module utilities
@windows? (
begin
    const libcuda = find_library(["nvcuda"], [""])
end
: # linux or mac
begin
    const libcuda = find_library(["libcuda"], ["/usr/lib/"])
end)

if isempty(libcuda)
    error("CUDA driver API library cannot be found")
end

include("version.jl")
include("types.jl")
include("module.jl")
include("device.jl")
include("stream.jl")
#include("event.jl")
include("pointer.jl")
include("arrays.jl")
include("execute.jl")

end


text:
module LengthUniqueTest

using HPAT

#HPAT.CaptureAPI.set_debug_level(3)
#HPAT.DomainPass.set_debug_level(3)
#HPAT.DataTablePass.set_debug_level(3)

@acc hpat function q25(d_date, file_name)
    store_sales = DataSource(DataTable{:ss_customer_sk=Int64, :ss_ticket_number=Int64, :ss_sold_date_sk=Int64, :ss_net_paid=Float64}, HDF5, file_name)

    store_agg = aggregate(store_sales, :cid = :ss_customer_sk,
                                           :frequency = length(unique(:ss_ticket_number)),
                                           :most_recent_date = maximum(:ss_sold_date_sk),
                                           :amount = sum(:ss_net_paid))

    return store_agg[:cid], store_agg[:frequency], store_agg[:most_recent_date], store_agg[:amount]
end

using HDF5
using MPI
using Base.Test

function main()

    if MPI.Comm_rank(MPI.COMM_WORLD)==0
        h5write("test_q25.hdf5", "/ss_customer_sk", [1,1,2,2,1])
        h5write("test_q25.hdf5", "/ss_ticket_number", [0,3,5,6,0])
        h5write("test_q25.hdf5", "/ss_sold_date_sk", [37580,34000,35000,36000,37600])
        h5write("test_q25.hdf5", "/ss_net_paid", [101.0,24.0,3.5,50.0,3.2])
    end

    c, f, m, a = q25("33000", "test_q25.hdf5")
    println(c,f,m,a)

    @test c==[1,2]
    @test f==[2,2]
    @test m==[37600,36000]
    @test_approx_eq a [128.2,53.5]

end

end

LengthUniqueTest.main()


text:
using CSV, CairoMakie, DataFrames, LinearAlgebra, Statistics, wgregseq
wgregseq.plotting_style.default_makie!()

# Import data. Change path if necessary
#using Pandas
#Pandas.to_csv(Pandas.read_excel("data/heinemann_data.xlsx", sheet_name="Table S6", header=[2]), "data/heinemann_data_cleaned.csv")
df = CSV.read("data/heinemann_data_cleaned.csv", DataFrames.DataFrame)
df = df[.~ismissing.(df.Gene), :]
df[!, "Annotated functional COG group (description)"] |> unique
##
growth_conditions = ["Glucose",
 "LB",
 "Glycerol + AA",
 "Acetate",
 "Fumarate",
 "Glucosamine",
 "Glycerol",
 "Pyruvate",
 "Chemostat µ=0.5",
 "Chemostat µ=0.35",
 "Chemostat µ=0.20",
 "Chemostat µ=0.12",
 "Stationary phase 1 day",
 "Stationary phase 3 days",
 "Osmotic-stress glucose",
 "42°C glucose",
 "pH6 glucose",
 "Xylose",
 "Mannose",
 "Galactose ",
 "Succinate",
 "Fructose"]

 # Compute mean copy number per growth condition
tot_prot = [sum(df[.~ismissing.(df[!, cond]), cond]) for cond in growth_conditions]

# Plot total protein number
fig = Figure(resolution = (600, 400))
ax1 = Axis(
    fig[1, 1],
    ylabel="Total number of proteins",
    xlabel="Growth Conditions"
)
# Sort growth conditions by total protein count
ind = sortperm(tot_prot, rev=true)
growth_conditions = growth_conditions[ind]

# Plot total protein number
lines!(ax1, 1:22, tot_prot[ind])
scatter!(ax1, 1:22, tot_prot[ind], )

ax1.xticklabelrotation = pi/8
ax1.xticks = 1:length(growth_conditions)
ax1.xtickformat = x -> growth_conditions
ax1.xminorgridwidth = 0

CairoMakie.save("figures/heinemann_total_protein.pdf", fig)
fig

## Function to plot protein number per gene
"""
    copy_number_plot(genes, normalized=False)

Plot the copy number of a list of genes in all
growth conditions.
"""
function copy_number_plot(genes, normalized=false)

    fig = Figure(resolution = (600, 400))
    ax1 = Axis(
        fig[1, 1],
        ylabel="Total number of proteins",
        xlabel="Growth Conditions"
    )

    # Iterate through genes
    for gene in genes
        # Extract row for gene
        if ~(gene in df[!, "Gene"])
            throw(ArgumentError("Gene $gene is not in the data set."))
        end
        row = df[df[!, "Gene"] .== gene, growth_conditions] |> Matrix |> vec
        mask = map(x -> ~ismissing(x), row)

        # Normalize
        if normalized
            row[mask] .= row[mask] ./ tot_prot[mask]
        end

        # Plot and make label
        lines!(ax1, collect(1:22)[mask], row[mask], label=gene)
        scatter!(ax1, collect(1:22)[mask], row[mask])
    end

    # Annotate plot
    if normalized
        ax1.ylabel = "Normalized Protein Copy Number"
    else
        ax1.ylabel = "Protein Copy Number"
    end

    ax1.xticklabelrotation = pi/8
    ax1.xticks = 1:length(growth_conditions)
    ax1.xtickformat = x -> growth_conditions
    ax1.xminorgridwidth = 0
    axislegend(ax1, position= :lt)

    return fig
end

fig = copy_number_plot(["galK", "galM"], true)
fig

##
function compute_variation(df)
    gene_vals = Matrix(df[!, growth_conditions])
    gene_vals = coalesce.(gene_vals, 0)
    return var(gene_vals ./ vec(tot_prot)', dims=2) ./ mean(gene_vals ./ vec(tot_prot)', dims=2)
end

df.variation = vec(compute_variation(df))
sort!(df, :variation, rev=true)
copy_number_plot(df.Gene[1:5])
unknown_ind = coalesce.(df[!, "Annotated functional COG group (description)"] .== "Function unknown", false)
df[unknown_ind, "Gene"][1:20]
copy_number_plot(df[unknown_ind, "Gene"][1:20], true)
## Compare two genes individually

function compare_genes(gene, growth_condition1, growth_condition2)
    if growth_condition1 ∉ growth_conditions
        throw(ArgumentError("Growth conditions $growth_condition1 not found!"))
    elseif growth_condition2 ∉ growth_conditions
        throw(ArgumentError("Growth conditions $growth_condition2 not found!"))
    end
    if (length(df[df[!, "Gene"] .== gene, growth_condition1]) > 0) && (length(df[df[!, "Gene"] .== gene, growth_condition2]) > 0)
        protein1 = df[df[!, "Gene"] .== gene, growth_condition1][1]
        protein2 = df[df[!, "Gene"] .== gene, growth_condition2][1]
        return protein1 / protein2 * sum(df[!, growth_condition2]) / sum(df[!, growth_condition1])
    else
        return NaN
    end
end
g1 = "LB"
g2 = "Stationary phase 1 day"
x = filter(x -> x > 0, compare_genes.(df[!, "Gene"], g1, g2))
fig = Figure(resolution=(400, 400))
ax = Axis(fig[1, 1], xscale = log10)
lines!(ax, sort(x), range(0, stop=1, length=length(x)))
ax.xlabel = "Fold Change"
ax.ylabel = "ECDF"
ax.title = "Fold Change of normalized protein copy number \n $g1/$g2"

fig

## Fold change of gene expression
df = df[.~ismissing.(df.Gene), :]

# Compute fold change
function get_all_fold_changes(df, gene, growth_conditions=growth_conditions)
    total_per_cond = [sum(df[.~ismissing.(df[!, cond]), cond]) for cond in growth_conditions]
    mask =  .~ ismissing.(collect(eachrow(df[df.Gene .== gene, growth_conditions])[1]))
    rel_abundance = collect(eachrow(df[df.Gene .== gene, growth_conditions])[1])[mask] ./ total_per_cond[mask]
    return rel_abundance' ./ rel_abundance, mask
end

# Choose gene and fold change limit
gene = "tolC"
lim = 2

mat, mask = get_all_fold_changes(df, gene)
fig = Figure(resolution=(32*sum(mask), 23*sum(mask)))
ax = Axis(fig)

heatmap!(ax, log.(mat), colorrange=(-log2(lim), log2(lim)))

ax.xticks = 1:sum(mask)
ax.xtickformat = x -> growth_conditions[mask]
ax.xminorgridwidth = 0
ax.xticklabelrotation = pi/8
ax.yticks = 1:sum(mask)
ax.ytickformat = x -> growth_conditions[mask]
ax.yminorgridwidth = 0
ax.yreversed = true
fig[1, 1] = ax
Colorbar(fig[1, 2], limits = (-log2(lim), log2(lim)), colormap = :viridis, label = "Log2 Fold Change")
ax.title = gene
fig
CairoMakie.save("figures/heineman_$(gene)_comp_to_all_fc$lim.pdf", fig)

##

inds_max = findall(x -> x .> lim, mat)
y_unique = unique([x[2] for x in inds_max])
x_unique = unique([x[1] for x in inds_max])
mat_reduced = mat[x_unique, y_unique]

fig = Figure(resolution=(700,500))
ax = Axis(fig, aspect = length(x_unique)/length(y_unique))

heatmap!(ax, log.(mat_reduced), colorrange=(-log2(lim), log2(lim)))

ax.xticks = 1:length(x_unique)
ax.xtickformat = x -> growth_conditions[mask][x_unique]
ax.xminorgridwidth = 0
ax.xticklabelrotation = pi/8
ax.yticks = 1:length(y_unique)
ax.ytickformat = x -> growth_conditions[mask][y_unique]
ax.yminorgridwidth = 0
fig[1, 1] = ax
Colorbar(fig[1, 2], limits = (-log2(lim), log2(lim)), colormap = :viridis, label = "Log2 Fold Change")
ax.title = gene
fig
CairoMakie.save("figures/heineman_$(gene)_comp_to_all_fc$(lim)_reduced.pdf", fig)

## Now we choose subsets of the growth conditions and repeat

sub_growth_conditions = ["Glucose",
 "LB",
 "Glycerol + AA",
 "Chemostat µ=0.5",
 "Chemostat µ=0.35",
 "Chemostat µ=0.20",
 "Chemostat µ=0.12",
 "Stationary phase 1 day",
 "Stationary phase 3 days",
 "Osmotic-stress glucose",
 "42°C glucose",
 "pH6 glucose"]

 carbon_growth_conditions = [
    "Glucose",
    "Acetate",
    "Fumarate",
    "Glucosamine",
    "Glycerol",
    "Pyruvate",
    "Xylose",
    "Mannose",
    "Galactose ",
    "Succinate",
    "Fructose"
    ]

function fold_change_to_any(df, growth_conditions)
    cutoff = 2:10
    rel_identified = Float64[]

    for i in cutoff
        identified = 0

        for gene in df.Gene
            mat = get_all_fold_changes(df, gene, growth_conditions)[1]
            mat = coalesce.(mat, 0)
            identified += any(mat .> i)
        end

        push!(rel_identified, identified / length(df.Gene))
    end
    return rel_identified
end

fig = Figure(resolution=(350, 350))
ax = Axis(fig)

rel_identified_all = fold_change_to_any(df, growth_conditions)
rel_identified_sub = fold_change_to_any(df, sub_growth_conditions)
rel_identified_carbon = fold_change_to_any(df, carbon_growth_conditions)

lines!(ax, cutoff, rel_identified_all, linewidth=2, label="All conditions")
lines!(ax, cutoff, rel_identified_sub, linewidth=2, label="Only Glucose as Carbon Source")
lines!(ax, cutoff, rel_identified_carbon, linewidth=2, label="Only varying Carbon")

axislegend()
ax.xlabel = "Fold Change"
ax.ylabel = "Identified Genes"
ax.yticks = 0:0.2:1
ax.title = "Compared to any \nother growth condition"
ylims!(ax, (-0.02, 1.02))
fig[1, 1] = ax
fig
CairoMakie.save("figures/heineman_identified_genes_comp_to_all.pdf", fig)

## Fold change of gene expression relative to median across conditions
df = df[.~ismissing.(df.Gene), :]
function get_fold_changes(df, gene, growth_conditions=growth_conditions)
    total = [sum(df[.~ismissing.(df[!, cond]), cond]) for cond in growth_conditions]
    mask = .~ ismissing.(collect(eachrow(df[df.Gene .== gene, growth_conditions])[1]))
    rel_abundance = collect(eachrow(df[df.Gene .== gene, growth_conditions])[1])[mask] ./ total[mask]
    return rel_abundance ./ median(rel_abundance), mask
end
gene = "galK"
get_fold_changes(df, gene)
fig = Figure(resolution=(600, 300))
ax = Axis(fig)

y, mask = get_fold_changes(df, gene)
lines!(ax, y, xticklabels=growth_conditions)

ax.xticks = 1:sum(mask)
ax.xtickformat = x -> growth_conditions[mask]
ax.xminorgridwidth = 0
ax.xticklabelrotation = pi/8
fig[1, 1] = ax
ax.title = gene
fig
CairoMakie.save("figures/heineman_$(gene)_rel_to_median.pdf", fig)
##

function fold_change_to_median(df, growth_conditions)
    cutoff = 2:10
    rel_identified = Float64[]
    for i in cutoff
        identified = 0
        for gene in df.Gene
            mat = get_fold_changes(df, gene, growth_conditions)[1]
            mat = coalesce.(mat, 0)
            identified += any((mat .> i) .| (1/i .> mat))
        end
        push!(rel_identified, identified / length(df.Gene))
    end
    return rel_identified
end

fig = Figure(resolution=(350, 350))
ax = Axis(fig)

rel_identified_all = fold_change_to_median(df, growth_conditions)
rel_identified_sub = fold_change_to_median(df, sub_growth_conditions)
rel_identified_carbon = fold_change_to_median(df, carbon_growth_conditions)

lines!(ax, cutoff, rel_identified_all, linewidth=2, label="All conditions")
lines!(ax, cutoff, rel_identified_sub, linewidth=2, label="Only Glucose as Carbon Source")
lines!(ax, cutoff, rel_identified_carbon, linewidth=2, label="Only varying Carbon")
axislegend()

ax.xlabel = "Fold Change"
ax.ylabel = "Identified Genes"
ax.title = "Compared to the median \nacross growth conditions"
ax.yticks = 0:0.2:1
ylims!(ax, (-0.02, 1.02))

fig[1, 1] = ax
fig
CairoMakie.save("figures/heineman_identified_genes_rel_to_median.pdf", fig)

## Find


text:
using Revise
using Pkg
if !isdefined(@__MODULE__, :CAS)
    pkg"""
    activate codegen
    """
end
module M
include("codegen/cas.jl")
using .CAS: sym_multivector, expr_fundef, SubAll, SubDeg, SubEven, wedge, expr_structdef
using .CAS: expr_structdef, exprs_structdef_batteries
using .CAS: exprs_everything, print_everything

#x = sym_multivector(:xx => SubDeg{1}())
#y = sym_multivector(:yy => SubDeg{1}())
#x*y
#y*x
#x*x

#print_everything(stdout)
open("src/auto.jl", "w") do io
    print_everything(io)
end

end


text:
using TopologyPreprocessing
using Test
using LinearAlgebra

# using MatrixOrganization

@testset "MatrixProcessing.jl -> basics" begin
    # Ints
    positive_matrix = [1 2 3;
                       4 5 6]
    negative_matrix = [1 2 -3;
                       -4 5 6]
    @test shift_to_non_negative(positive_matrix) == positive_matrix
    @test isempty(findall(x->x<0, shift_to_non_negative(negative_matrix)))

    # Floats
    positive_matrix2 = [1. 2 3; 4 5 6]
    negative_matrix2 = [1 2 -3; -4 5 6]
    @test shift_to_non_negative(positive_matrix2) == positive_matrix2
    @test isempty(findall(x->x<0, shift_to_non_negative(negative_matrix2)))

    positive_3d_matrix = rand(2,3,4)
    negative_3d_matrix = rand(2,3,4).*2 .-1
    @test shift_to_non_negative(positive_3d_matrix) == positive_3d_matrix
    @test isempty(findall(x->x<0, shift_to_non_negative(negative_3d_matrix)))

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

    @test findmin(normalize_to_01(negative_matrix))[1] == 0
    @test findmax(normalize_to_01(negative_matrix))[1] == 1

    powers_of_2 = [0 ; [2^k for k in 0:2:8]]
    powers_of_3 = [0 ; [3^k for k in 0:2:8]]
    @test findmin(normalize_to_01(powers_of_2, use_factor=true))[1] == 0
    @test findmax(normalize_to_01(powers_of_2, use_factor=true))[1] == 1

    @test findmin(normalize_to_01(powers_of_3, use_factor=true))[1] == 0
    @test findmax(normalize_to_01(powers_of_3, use_factor=true))[1] > 1

    @test findmin(normalize_to_01(powers_of_3, use_factor=true, norm_factor=3^8))[1] == 0
    @test findmax(normalize_to_01(powers_of_3, use_factor=true, norm_factor=3^8))[1] == 1

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

    square_matrix = [1 2 3;
                     4 5 6;
                     7 8 9]
    @test issymmetric(diagonal_symmetrize(square_matrix))
    @test LinearAlgebra.checksquare(diagonal_symmetrize(square_matrix)) == 3

    @test issymmetric(diagonal_symmetrize(positive_matrix))
    @test LinearAlgebra.checksquare(diagonal_symmetrize(positive_matrix)) == 2

    @test diagonal_symmetrize(square_matrix)[end,1] == square_matrix[1,end]
    @test diagonal_symmetrize(square_matrix)[2,1] == square_matrix[1,2]

    @test diagonal_symmetrize(square_matrix, below_over_upper=true)[1,end] == square_matrix[end,1]
    @test diagonal_symmetrize(square_matrix, below_over_upper=true)[1,2] == square_matrix[2,1]
end

@testset "MatrixProcessing.jl -> distances and indices" begin
    # Ints
    positive_matrix = [1 2 3;4 5 6]
    positive_matrix2 = [1. 2 3; 4 5 6]
    positive_3d_matrix = rand(2,3,4)

    negative_matrix = [1 2 -3;-4 5 6]
    negative_matrix2 = [1 2 -3; -4 5 6]
    negative_3d_matrix = rand(2,3,4).*2 .-1

    negative_6d_matrix = rand(2,3,4,5,6,7).*2 .-1

    powers_of_2 = [0 ; [2^k for k in 0:2:8]]
    powers_of_3 = [0 ; [3^k for k in 0:2:8]]

    square_matrix = [1 2 3;
                     4 5 6;
                     7 8 9]

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

    @test length(unique(group_distances(square_matrix,1))) == 1
    @test length(unique(group_distances(square_matrix,2))) == 2
    @test length(unique(group_distances(square_matrix,9))) == 9
    @test_throws DomainError group_distances(square_matrix,20)

    @test length(unique(group_distances(positive_matrix,1))) == 1
    @test length(unique(group_distances(positive_matrix,2))) == 2
    @test length(unique(group_distances(positive_matrix,6))) == 6
    @test_throws DomainError group_distances(positive_matrix,20)

    @test length(unique(group_distances(positive_3d_matrix,2))) == 2
    @test length(unique(group_distances(negative_3d_matrix,2))) == 2
    @test length(unique(group_distances(negative_6d_matrix,2))) == 2

    group_distances(square_matrix,2)

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

    @test length(generate_indices(3))==3^2
    @test length(generate_indices(223))==223^2

    n=3
    @test length(generate_indices(n, symmetry_order=true, include_diagonal=false)) == ((n^2 - n) ÷ 2)
    n=9
    @test length(generate_indices(n, symmetry_order=true, include_diagonal=false)) == ((n^2 - n) ÷ 2)
    index_matrix1 = generate_indices(n, symmetry_order=true, include_diagonal=false)
    @test length(findall(x-> x == CartesianIndex(n,n),index_matrix1)) == 0
    @test length(findall(x-> x == CartesianIndex(n-2,n-1),index_matrix1)) == 1
    @test length(findall(x-> x == CartesianIndex(n-1,n-2),index_matrix1)) == 0

    n=223
    @test length(generate_indices(n, symmetry_order=true, include_diagonal=false)) == ((n^2 - n) ÷ 2)
    @test length(generate_indices(n, symmetry_order=true, include_diagonal=true)) == ((n^2 - n) ÷ 2 + n)

    n=3
    index_matrix2 = generate_indices(n, symmetry_order=true, include_diagonal=true)
    @test length(index_matrix2) == (((n-1)*n)÷2 +n)
    @test findmax(index_matrix2)[1] == CartesianIndex(n,n)
    @test length(findall(x-> x == CartesianIndex(1,n),index_matrix2)) == 1
    @test length(findall(x-> x == CartesianIndex(1,1),index_matrix2)) == 1
    @test length(findall(x-> x == CartesianIndex(n,n),index_matrix2)) == 1

    n=9
    @test length(generate_indices(n, symmetry_order=true, include_diagonal=true)) == (((n-1)*n)÷2 +n)
    n=223
    @test length(generate_indices(n, symmetry_order=true, include_diagonal=true)) == (((n-1)*n)÷2 +n)

    n=9
    @test length(generate_indices((n,n), symmetry_order=true, include_diagonal=true)) == (((n-1)*n)÷2 +n)
    n=223
    @test length(generate_indices((n,n), symmetry_order=true, include_diagonal=true)) == (((n-1)*n)÷2 +n)
end

# @testset "MatrixProcessing.jl -> arrays of arrays" begin
#     # Need to be taken from Bettis
#     # arr_of_arrs
#
#     # reduce_arrs_to_min_len(arr_of_arrs)
#
# end

@testset "MatrixProcessing.jl -> matrix ordering helping functions" begin
    let testing_matrix0 = Array{Float64,2}(undef, 2, 3),
        testing_matrix1 = [1 2 3; 4 5 6; 7 8 9],
        testing_matrix2 = ones((2,3,4))
        testing_matrix3 = [1 2 3; 4 5 6]
        testing_matrix4 = [1, 4, 7, 2, 5, 8, 3, 6, 9]

        @test arr_to_vec(testing_matrix0) isa Vector
        @test length(testing_matrix0) == length(arr_to_vec(testing_matrix0))

        @test arr_to_vec(testing_matrix1) isa Vector
        @test length(testing_matrix1) == length(arr_to_vec(testing_matrix1))
        @test arr_to_vec(testing_matrix1) == [1, 4, 7, 2, 5, 8, 3, 6, 9]

        @test arr_to_vec(testing_matrix2) isa Vector
        @test length(testing_matrix2) == length(arr_to_vec(testing_matrix2))

        @test arr_to_vec(testing_matrix3) isa Vector
        @test length(testing_matrix3) == length(arr_to_vec(testing_matrix3))
        @test arr_to_vec(testing_matrix3) == [1, 4, 2, 5, 3, 6]

        @test arr_to_vec(testing_matrix4) isa Vector
        @test length(testing_matrix4) == length(arr_to_vec(testing_matrix4))
        @test arr_to_vec(testing_matrix4) == [1, 4, 7, 2, 5, 8, 3, 6, 9]
    end

    let testing_matrix1 = CartesianIndices((2,2)),
        testing_matrix2 = CartesianIndices((9,1))

        @test cartesianInd_to_vec(testing_matrix1) isa Vector
        @test length(cartesianInd_to_vec(testing_matrix1)) == length(testing_matrix1)

        @test cartesianInd_to_vec(testing_matrix2) isa Vector
        @test length(cartesianInd_to_vec(testing_matrix2)) == length(testing_matrix2)
    end

    let vals_matrix1a = [1 2 3; 4 5 6],
        vals_matrix1b = [6 5 4; 3 2 1],
        vals_matrix1c = [4 5 6; 1 2 3]
        let index_matrix1a = [CartesianIndex(1,1), CartesianIndex(1,2), CartesianIndex(1,3), CartesianIndex(2,1), CartesianIndex(2,2), CartesianIndex(2,3)]
            @test length(sort_indices_by_values(vals_matrix1a, index_matrix1a)) == length(vals_matrix1a)
            @test sort_indices_by_values(vals_matrix1a, index_matrix1a) == collect(1:6)

            @test length(sort_indices_by_values(vals_matrix1b, index_matrix1a)) == length(vals_matrix1b)
            @test sort_indices_by_values(vals_matrix1b, index_matrix1a) == collect(6:-1:1)

            @test length(sort_indices_by_values(vals_matrix1c, index_matrix1a)) == length(vals_matrix1c)
            @test sort_indices_by_values(vals_matrix1c, index_matrix1a) == [4, 5, 6, 1, 2, 3]
        end

        let index_matrix1b = CartesianIndices((2,3))
            @test_throws TypeError sort_indices_by_values(vals_matrix1a, index_matrix1b)
        end
    end

    let vals_matrix2 = [1 2 3; 4 5 6; 7 8 9],
        index_matrix2a = CartesianIndices((3,3)),
        index_matrix2b =   [CartesianIndex(1,1) CartesianIndex(1,2) CartesianIndex(1,3);
                            CartesianIndex(2,1) CartesianIndex(2,2) CartesianIndex(2,3);
                            CartesianIndex(3,1) CartesianIndex(3,2) CartesianIndex(3,3)],
        index_matrix2c =   [CartesianIndex(1,1), CartesianIndex(1,2), CartesianIndex(1,3),
                            CartesianIndex(2,1), CartesianIndex(2,2), CartesianIndex(2,3),
                            CartesianIndex(3,1), CartesianIndex(3,2), CartesianIndex(3,3)]

        @test_throws TypeError sort_indices_by_values(vals_matrix2, index_matrix2a)
        @test_throws TypeError sort_indices_by_values(vals_matrix2, index_matrix2b)

        @test sort_indices_by_values(vals_matrix2, index_matrix2c) isa Vector
        @test sort_indices_by_values(vals_matrix2, index_matrix2c) == 1:9
        @test length(sort_indices_by_values(vals_matrix2, index_matrix2c)) == length(vals_matrix2)
    end

    let vals_matrix3 = [1, 4, 7, 2, 5, 8, 3, 6, 9],
        index_matrix3a = CartesianIndices((9,1)),
        index_matrix3b = CartesianIndices((9,)),
        index_matrix3c = [1, 4, 7, 2, 5, 8, 3, 6, 9]

        @test_throws TypeError sort_indices_by_values(vals_matrix3, index_matrix3a)
        @test_throws TypeError sort_indices_by_values(vals_matrix3, index_matrix3b)

        @test sort_indices_by_values(vals_matrix3, index_matrix3c) isa Vector
        @test sort_indices_by_values(vals_matrix3, index_matrix3c) == 1:9
        @test length(sort_indices_by_values(vals_matrix3, index_matrix3c)) == length(vals_matrix3)
    end

    let target_coords1 = CartesianIndex(2,3),
        target_value = -20

        let some_matrix = [1 2 3; 4 5 6; 7 8 9]

            set_values!(some_matrix, target_coords1, target_value; do_symmetry=false)
            @test some_matrix[target_coords1] == target_value
        end

        let some_matrix = [1 2 3; 4 5 6; 7 8 9]

            another_matrix = set_values!(some_matrix, target_coords1, target_value; do_symmetry=false)
            @test some_matrix[target_coords1] == target_value
            @test another_matrix[target_coords1] == target_value
            @test another_matrix === some_matrix
        end

        let some_matrix = [1 2 3; 4 5 6; 7 8 9]

            another_matrix = set_values!(some_matrix, target_coords1, target_value; do_symmetry=true)
            @test some_matrix[target_coords1] == target_value
            @test some_matrix[target_coords1[1],target_coords1[2]] == target_value
            @test some_matrix[target_coords1[2],target_coords1[1]] == target_value
            @test another_matrix === some_matrix
        end

        let some_matrix = [1 2 3; 4 5 6; 7 8 9],
            some_matrix2 = [1 2 3; 4 5 6; 7 8 9],
            target_coords2 = CartesianIndex(8,9)

            @test_throws BoundsError set_values!(some_matrix, target_coords2, target_value; do_symmetry=false)
            @test some_matrix == some_matrix2
        end

        let some_matrix = [1 2 3; 4 5 6; 7 8 9],
            target_coords3 = CartesianIndices((2,2))

            @test_throws MethodError set_values!(some_matrix, target_coords3, target_value; do_symmetry=false)
        end
    end
end

@testset "MatrixProcessing.jl -> matrix ordering" begin

    """
        check_for_min_val_position(input_matrix::Matrix; force_symmetry=false, assign_same_values=false)

    Takes an 'input_matrix' and checks if its ordered form has the minimum value
    in the same position as the minimum value of 'input_matrix'.

    """
    function check_for_min_val_position(input_matrix::Matrix; force_symmetry=false, assign_same_values=false)
        # check if min val in uniquely value matrix is in the same position
        ordered_matrix = get_ordered_matrix(input_matrix, force_symmetry=force_symmetry, assign_same_values=assign_same_values)
        if issymmetric(input_matrix) || force_symmetry
            off_diag_ind = generate_indices(size(input_matrix),include_diagonal=false)
        else
            off_diag_ind = generate_indices(size(input_matrix),include_diagonal=true)
        end
        min_val = findmin(input_matrix[off_diag_ind])[1]
        # min_orig = findmin(input_matrix[off_diag_ind])[2]
        all_min_input = findall(x->x==min_val,input_matrix[off_diag_ind])
        all_min_ordered = findall(x->x==1,ordered_matrix[off_diag_ind])
        return off_diag_ind[all_min_input] == off_diag_ind[all_min_ordered]
    end

    # ==
    let power_matrix = zeros(Int, (2,3,4))
        for k = 1:4
            power_matrix[:,:,k] = reshape([(k+1)^n for n =1:6], (2,3))
        end

        ordered_power_matrix = copy(power_matrix)
        ordered_power_matrix[:, :, 1] =[1  6  12;
                                        3  8  13]
        ordered_power_matrix[:, :, 2] =[2  11  17;
                                            7  15  20]
        ordered_power_matrix[:, :, 3] = [4  14  21;
                                            9  18  23]
        ordered_power_matrix[:, :, 4] =[5  16  22;
                                        10  19  24]

        @test !isempty(get_ordered_matrix(power_matrix) == ordered_power_matrix)
        @test get_ordered_matrix(power_matrix) == ordered_power_matrix

        # @test_throws MethodError get_ordered_matrix(power_matrix)
    end

    # ==
    let square_matrix1 = [1 2 3; 4 5 6; 7 8 9],
        square_matrix2 = [1 4 7; 2 5 8; 3 6 9]
        @test get_ordered_matrix(10 .*square_matrix1) == (square_matrix1 )
        @test get_ordered_matrix(10 .*square_matrix2) == (square_matrix2 )
        @test sum(get_ordered_matrix(square_matrix1) .== square_matrix1 ) == 9
        @test sum(get_ordered_matrix(square_matrix2) .== square_matrix2 ) == 9
        @test get_ordered_matrix(10 .*square_matrix1, force_symmetry=true) == get_ordered_matrix(square_matrix1, force_symmetry=true)
        @test get_ordered_matrix(10 .*square_matrix2, force_symmetry=true) == get_ordered_matrix(square_matrix2, force_symmetry=true)

        # check if min val in uniquely value matrix is in the same position
        let input_matrix = 10square_matrix1
            @test check_for_min_val_position(input_matrix)
        end
    end

    # ==
    let square_matrix_same_vals1 = [1 2 3; 3 4 5; 6 7 8],
        square_matrix_same_vals2 = [1 3 6; 2 4 7; 3 5 8]

        @test get_ordered_matrix(10 .*square_matrix_same_vals1, assign_same_values=true) == (square_matrix_same_vals1 )
        @test get_ordered_matrix(10 .*square_matrix_same_vals2, assign_same_values=true) == (square_matrix_same_vals2 )

        # forcing symmetry test
        some_ord_mat = get_ordered_matrix(10 .*square_matrix_same_vals1, force_symmetry=true, assign_same_values=false)
        # remove 1, because 0 is not off diagonal
        @test length(unique(some_ord_mat))-1 == (size(square_matrix_same_vals1,1)*(size(square_matrix_same_vals1,1)-1))/2
    end

    # ==
    let square_matrix_same_vals3 = [1 2 3; 3 4 3; 5 6 7],
        square_matrix_same_vals4 = [1 3 3; 2 4 6; 3 3 7]
        @test get_ordered_matrix(10 .*square_matrix_same_vals3, force_symmetry=true, assign_same_values=true) ==
                            get_ordered_matrix(square_matrix_same_vals3, force_symmetry=true, assign_same_values=true)
        @test get_ordered_matrix(10 .*square_matrix_same_vals4, force_symmetry=true, assign_same_values=true) ==
                            get_ordered_matrix(square_matrix_same_vals4, force_symmetry=true, assign_same_values=true)
    end

    # ==================
    # Tests on symmetric matrices
    let test_mat_b1 = [ 1  1  1  4  5  9 ;
                        1  1  2  3  6  10;
                        1  2  1  3  7  11;
                        4  3  3  1  8  12;
                        5  6  7  8  1  13;
                        9  10 11 12 13 1 ;]

        # no same values
        let test_mat_b1_ord1 = [  0   1   2   6   7  11;
                                  1   0   3   4   8  12;
                                  2   3   0   5   9  13;
                                  6   4   5   0  10  14;
                                  7   8   9  10   0  15;
                                 11  12  13  14  15   0],
            test_mat_b1_indices = generate_indices(size(test_mat_b1))

            ord_mat_b1_1 = get_ordered_matrix(10 .*test_mat_b1, force_symmetry=false, assign_same_values=false)
            @test issymmetric(ord_mat_b1_1)
            @test ord_mat_b1_1[test_mat_b1_indices] == test_mat_b1_ord1[test_mat_b1_indices]

            ord_mat_b1_2 = get_ordered_matrix(10 .*test_mat_b1, force_symmetry=true, assign_same_values=false)
            @test issymmetric(ord_mat_b1_2)
            @test ord_mat_b1_2[test_mat_b1_indices] == test_mat_b1_ord1[test_mat_b1_indices]
        end

        # assign same values
        let test_mat_b1_ord2 = [ 0  1  1  4  5  9 ;
                                 1  0  2  3  6  10;
                                 1  2  0  3  7  11;
                                 4  3  3  0  8  12;
                                 5  6  7  8  0  13;
                                 9 10  11 12 13 0 ;]

            let ord_mat_b1_3 = get_ordered_matrix(10 .*test_mat_b1, force_symmetry=false, assign_same_values=true)
                @test issymmetric(ord_mat_b1_3)
                # Removed, because ordered diagonal is all 1: @test ord_mat_b1_3 == (test_mat_b1 )
                @test ord_mat_b1_3 == test_mat_b1_ord2
            end

            let ord_mat_b1_4 = get_ordered_matrix(10 .*test_mat_b1, force_symmetry=true, assign_same_values=true)
                @test issymmetric(ord_mat_b1_4)
                # Removed, because ordered diagonal is all 1: @test ord_mat_b1_4 == (test_mat_b1 )
                @test ord_mat_b1_4 == test_mat_b1_ord2
            end
        end

        let input_matrix = 10 .*test_mat_b1
            @test check_for_min_val_position(input_matrix; force_symmetry=false, assign_same_values=true)
        end
    end

    # ==
    # Non-symmetric matrix test
    let test_mat_b2 = [ 1  1  3  4  5  9 ;
                        1  1  2  3  6  10;
                        14 2  1  3  7  11;
                        4  15 3  1  8  12;
                        5  5  7  8  1  13;
                        9  10 11 12 13 1 ;]

        let ord_mat_b2_1 = get_ordered_matrix(10 .*test_mat_b2, force_symmetry=false, assign_same_values=false)
            @test !issymmetric(ord_mat_b2_1)
            @test findall(x->x<=8,ord_mat_b2_1) == findall(x->x==1,test_mat_b2) # all values with one are first 8 values used for ordering
            @test length(unique(ord_mat_b2_1)) == length(test_mat_b2) #check if all values are unique
        end

        let test_mat_b2_ord1 = [0  1  6  4  7  11;
                                1  0  2  5  8  12;
                                6  2  0  3  9  13;
                                4  5  3  0 10  14;
                                7  8  9 10  0  15;
                                11 12 13 14 15 0 ;]

            test_mat_b2_indices = generate_indices(size(test_mat_b2))
            filter!(x->x!=CartesianIndex(1,3), test_mat_b2_indices)
            filter!(x->x!=CartesianIndex(1,4), test_mat_b2_indices)
            filter!(x->x!=CartesianIndex(2,4), test_mat_b2_indices)
            filter!(x->x!=CartesianIndex(3,4), test_mat_b2_indices)
            filter!(x->x!=CartesianIndex(3,1), test_mat_b2_indices)
            filter!(x->x!=CartesianIndex(4,1), test_mat_b2_indices)
            filter!(x->x!=CartesianIndex(4,2), test_mat_b2_indices)
            filter!(x->x!=CartesianIndex(4,3), test_mat_b2_indices)

            ord_mat_b2_2 = get_ordered_matrix(10 .*test_mat_b2, force_symmetry=true, assign_same_values=false)
            @test issymmetric(ord_mat_b2_2)
            @test ord_mat_b2_2[test_mat_b2_indices] == test_mat_b2_ord1[test_mat_b2_indices]
            # forcing symmetry test:
            @test length(unique(ord_mat_b2_2))-1 == (size(test_mat_b2,1)*(size(test_mat_b2,1)-1))/2
        end

        # TODO to fix tests, add 1 to all values off diagonal for the input matrix
        let test_mat_b2_ord2 = [0  0  2  3  4  8;
                                0  0  1  2  5  9;
                                13 1  0  2  6  10;
                                3  14 2  0  7  11;
                                4  4  6  7  0  12;
                                8  9  10 11 12 0 ;]
            ord_mat_b2_3 = get_ordered_matrix(10 .*test_mat_b2, force_symmetry=false, assign_same_values=true)
            @test_skip !issymmetric(ord_mat_b2_3)
            @test_skip  ord_mat_b2_3 == test_mat_b2_ord2
        end

        # TODO to fix tests, add 1 to all values off diagonal for the input matrix
        let test_mat_b2_ord3 = [0  0  2   3   4  8
                                0  0  1   2   5  9
                                2  1  0   2   6  10
                                3  2  2   0   7  11
                                4  5  6   7   0  12
                                8  9  10  11  12 0 ;]
            ord_mat_b2_4 = get_ordered_matrix(10 .*test_mat_b2, force_symmetry=true, assign_same_values=true)
            @test_skip issymmetric(ord_mat_b2_4)
            @test_skip ord_mat_b2_4 == test_mat_b2_ord3
        end
    end

    # ==
    let test_mat_b3 = [ 1  1  3  4  5  9 ;
                        1  1  2  3  6  10;
                        3  2  1  3  7  11;
                        4  3  3  1  8  12;
                        5  6  7  8  1  13;
                        9  10 11 12 13 1 ;]

        let test_mat_b3_ord1 = [ 0  0  5  3  6  10;
                                    0  0  1  4  7  11;
                                    5  1  0  2  8  12;
                                    3  4  2  0  9  13;
                                    6  7  8  9  0  14;
                                    10 11 12 13 14 0 ;],
          test_mat_b3_indices = generate_indices(size(test_mat_b3))
            filter!(x->x!=CartesianIndex(1,3), test_mat_b3_indices)
            filter!(x->x!=CartesianIndex(1,4), test_mat_b3_indices)
            filter!(x->x!=CartesianIndex(2,4), test_mat_b3_indices)
            filter!(x->x!=CartesianIndex(3,4), test_mat_b3_indices)
            filter!(x->x!=CartesianIndex(3,1), test_mat_b3_indices)
            filter!(x->x!=CartesianIndex(4,1), test_mat_b3_indices)
            filter!(x->x!=CartesianIndex(4,2), test_mat_b3_indices)
            filter!(x->x!=CartesianIndex(4,3), test_mat_b3_indices)

            ord_mat_b3_1 = get_ordered_matrix(10 .*test_mat_b3, force_symmetry=false, assign_same_values=false)
            @test issymmetric(ord_mat_b3_1)
            @test_skip ord_mat_b3_1[test_mat_b3_indices] == test_mat_b3_ord1[test_mat_b3_indices]

            ord_mat_b3_2 = get_ordered_matrix(10 .*test_mat_b3, force_symmetry=true, assign_same_values=false)
            @test issymmetric(ord_mat_b3_2)
            @test_skip ord_mat_b3_2[test_mat_b3_indices] == test_mat_b3_ord1[test_mat_b3_indices]
        end

        # TODO remove tests that do not add anything new for testing and are just another similar case
        let test_mat_b3_ord2 = [ 0  0  2  3  4  8 ;
                                 0  0  1  2  5  9 ;
                                 2  1  0  2  6  10;
                                 3  2  2  0  7  11;
                                 4  5  6  7  0  12;
                                 8  9  10 11 12 0 ;]

            ord_mat_b3_3 = get_ordered_matrix(10 .*test_mat_b3, force_symmetry=false, assign_same_values=true)
            @test issymmetric(ord_mat_b3_3)
            @test_skip ord_mat_b3_3 == (test_mat_b3 .-1)
            @test_skip ord_mat_b3_3 == test_mat_b3_ord2

            ord_mat_b3_4 = get_ordered_matrix(10 .*test_mat_b3, force_symmetry=true, assign_same_values=true)
            @test issymmetric(ord_mat_b3_4)
            @test_skip ord_mat_b3_4 == (test_mat_b3 .-1)
            @test_skip ord_mat_b3_4 == test_mat_b3_ord2
        end

        let input_matrix = 10 .*test_mat_b3
            @test check_for_min_val_position(input_matrix; force_symmetry=false, assign_same_values=true)
        end
    end

    # ==
    let test_mat_b4 = [  1  1  41 4  5  9  13 17 25 33;
                        1  1  2  42 6  10 14 18 26 34;
                        41 2  1  3  7  11 15 19 27 35;
                        4  42 3  1  8  12 16 20 28 36;
                        5  6  7  8  1  21 43 24 29 37;
                        9  10 11 12 21 1  22 44 30 38;
                        13 14 15 16 43 22 1  23 31 39;
                        17 18 19 20 24 44 23 1  32 40;
                        25 26 27 28 29 30 31 32 1  45;
                        33 34 35 36 37 38 39 40 45 1;]

        @test_skip get_ordered_matrix(10 .*test_mat_b4, force_symmetry=false, assign_same_values=false) == (test_mat_b4 .-1)
        @test_skip get_ordered_matrix(10 .*test_mat_b4, force_symmetry=false, assign_same_values=true) == (test_mat_b4 .-1)

        let ord_mat = get_ordered_matrix(10 .*test_mat_b4, force_symmetry=true, assign_same_values=false)
            @test issymmetric(ord_mat)
            @test_skip ord_mat == (test_mat_b4 .-1)
        end

        let ord_mat = get_ordered_matrix(10 .*test_mat_b4, force_symmetry=true, assign_same_values=true)
            @test issymmetric(ord_mat)
            @test_skip  ord_mat == (test_mat_b4 .-1)
        end

        let input_matrix = 10test_mat_b4
            @test check_for_min_val_position(input_matrix)
        end
    end

    # ==
    let test_mat_b5 = -[1  1  3  4  5  9 ;
                        1  1  2  3  6  10;
                        14 2  1  3  7  11;
                        4  15 3  1  8  12;
                        5  5  7  8  1  13;
                        9  10 11 12 13 1 ;]
        let ord_mat_b5_1 = get_ordered_matrix(10 .*test_mat_b5, force_symmetry=false, assign_same_values=false)
            @test !issymmetric(ord_mat_b5_1)
            @test_skip findall(x->x>=28,ord_mat_b5_1) == findall(x->x==-1,test_mat_b5) # all values with one are first 8 values used for ordering
            @test length(unique(ord_mat_b5_1)) == length(test_mat_b5) #check if all values are unique
        end

        let test_mat_b5_ord1 = [ 0  14  9 11  8  4
                                 14 0  13 10  7  3
                                 9  13 0  12  6  2
                                 11 10 12 0   5  1
                                 8  7  6  5   0  0
                                 4  3  2  1   0  0 ;]

            test_mat_b5_indices = generate_indices(size(test_mat_b5))
            filter!(x->x!=CartesianIndex(1,3), test_mat_b5_indices)
            filter!(x->x!=CartesianIndex(1,4), test_mat_b5_indices)
            filter!(x->x!=CartesianIndex(2,4), test_mat_b5_indices)
            filter!(x->x!=CartesianIndex(3,4), test_mat_b5_indices)
            filter!(x->x!=CartesianIndex(3,1), test_mat_b5_indices)
            filter!(x->x!=CartesianIndex(4,1), test_mat_b5_indices)
            filter!(x->x!=CartesianIndex(4,2), test_mat_b5_indices)
            filter!(x->x!=CartesianIndex(4,3), test_mat_b5_indices)

            ord_mat_b5_2 = get_ordered_matrix(10 .*test_mat_b5, force_symmetry=true, assign_same_values=false)

            @test issymmetric(ord_mat_b5_2)
            @test_skip ord_mat_b5_2[test_mat_b5_indices] == test_mat_b5_ord1[test_mat_b5_indices]
            # forcing symmetry test:
            @test length(unique(ord_mat_b5_2))-1 == (size(test_mat_b5,1)*(size(test_mat_b5,1)-1))/2
        end

        let test_mat_b5_ord2 = [14 14 12 11 10 6;
                                14 14 13 12 9  5;
                                1  13 14 12 8  4;
                                11 0  12 14 7  3;
                                10 10 8  7  14 2;
                                6  5  4  3  2  14],
            ord_mat_b5_3 = get_ordered_matrix(10 .*test_mat_b5, force_symmetry=false, assign_same_values=true)

            @test !issymmetric(ord_mat_b5_3)
            @test_skip ord_mat_b5_3 == test_mat_b5_ord2
        end

        let test_mat_b5_ord3 = [0  12 10 9  8  4;
                                12 0  11 10 7  3;
                                10 11 0  10 6  2;
                                9  10 10 0  5  1;
                                8  7  6  5  0  0;
                                4  3  2  1  0  0]
            ord_mat_b5_4 = get_ordered_matrix(10 .*test_mat_b5, force_symmetry=true, assign_same_values=true)
            @test issymmetric(ord_mat_b5_4)
            @test_skip ord_mat_b5_4 == test_mat_b5_ord3
        end
    end
    # ==================
end


text:
if !isdefined(:__EXPRESSION_HASHES__)
    __EXPRESSION_HASHES__ = Set{UInt64}()
end

macro runonce(expr)
    h = hash(expr)
    return esc(quote
        if !in($h, __EXPRESSION_HASHES__)
            push!(__EXPRESSION_HASHES__, $h)
            $expr
        end
    end)
end

function flatten!(b::Vector, a::Vector)
    for x in a
        if isa(x, Array)
            flatten!(b, x)
        else
            push!(b, x)
        end
    end
    return b
end

flatten(a::Vector) = flatten!(eltype(a)[], a)
flatten{T}(::Type{T}, a::Vector) = convert(Vector{T}, flatten(a))


text:
module Dispatch

export DFTArray, DFTAxisArray
using AxisArrays
using Unitful

macro lintpragma(s) end
@lintpragma("Ignore use of undeclared variable T")
@lintpragma("Ignore use of undeclared variable N")

""" An axis array with dense memory allocation """
const DFTAxisArray = AxisArray{T, N, <: DenseArray{T, N}, A} where A where N where T
""" Set of array types used in LibXC and AtomicDFT """
const DFTArray =
    Union{DenseArray{T, N}, AxisArray{T, N, <: DenseArray{T, N}}} where N where T

""" Dimensions and units of interest for DFT """
const Dimensioned = (
    :Density                         => :ρ,
    :ContractedDensityGradient       => :σ,
    :FirstDensityDerivative          => :∂ϵ_∂ρ,
    :FirstGradientDerivative         => :∂ϵ_∂σ,
    :SecondDensityDerivative         => :∂²ϵ_∂ρ²,
    :SecondGradientDerivative        => :∂²ϵ_∂σ²,
    :SecondDensityGradientDerivative => :∂²ϵ_∂ρ∂σ,
    :ThirdDensityDerivative          => :∂³ϵ_∂ρ³,
    :ThirdGradientDerivative         => :∂³ϵ_∂σ³,
    :ThirdDensity2GradientDerivative => :∂³ϵ_∂ρ²∂σ,
    :ThirdDensityGradient2Derivative => :∂³ϵ_∂ρ∂σ²,
    :Potential                       => :Eₕ
)
include("DispatchDimensions.jl")
include("DispatchHartree.jl")

""" Convenient access to scalar physical dimension types """
const Scalars = Dimensions.Scalars
end

max_stars_repo_name: mdavezac/DFTShims

max_issues_repo_name: mdavezac/DFTShims

max_forks_repo_name: mdavezac/DFTShims


text:
module Gaius

using LoopVectorization: @avx, VectorizationBase.AbstractStridedPointer, VectorizationBase.gesp, VectorizationBase.vload, VectorizationBase.vstore!, VectorizationBase.AVX512F
import LoopVectorization: @avx, VectorizationBase.stridedpointer
using StructArrays: StructArray
using LinearAlgebra: LinearAlgebra

export blocked_mul, blocked_mul!

const DEFAULT_BLOCK_SIZE = AVX512F ? 96 : 64

const Eltypes  = Union{Float64, Float32, Int64, Int32, Int16}
const MatTypesC{T <: Eltypes} = Union{Matrix{T}, SubArray{T, 2, <: Array}} # C for Column Major
const MatTypesR{T <: Eltypes} = Union{LinearAlgebra.Adjoint{T,<:MatTypesC{T}}, LinearAlgebra.Transpose{T,<:MatTypesC{T}}} # R for Row Major
const MatTypes{T <: Eltypes} = Union{MatTypesC{T}, MatTypesR{T}}

# Note this does not support changing the number of threads at runtime
macro _spawn(ex)
    if Threads.nthreads() > 1
        esc(Expr(:macrocall, Expr(:(.), :Threads, QuoteNode(Symbol("@spawn"))), __source__, ex))
    else
        esc(ex)
    end
end
macro _sync(ex)
    if Threads.nthreads() > 1
        esc(Expr(:macrocall, Symbol("@sync"), __source__, ex))
    else
        esc(ex)
    end
end

include("pointermatrix.jl")
include("matmul.jl")
include("block_operations.jl")
include("kernels.jl")

end # module


text:
"""
# Seis.Synth

Create synthetic array data for testing.
"""
module Synth

using Statistics: mean

import Seis
import Geodesics
import ..Beamforming

export synthetic_arrival

"""
    synthetic_arrival(sx, sy; delta=0.01, nsta=20, lon, lat, distance, noise) -> ::Vector{Seis.Trace}

Calculate synthetic traces with an arrival at a set slowness (`sx`, `sy`) s/°.

Optionally specify the number of stations `nsta`, the trace sampling interval
`delta` s and the station coordinates `lon`° and `lat`°.

If `distance` is not `nothing`, then populate the traces' event info with
the coordinates of the event.  Specify the 2σ amplitude of Gaussian `noise`
to add to the trace as a value greater than 0.
"""
function synthetic_arrival(sx, sy;
        delta = 0.01,
        nsamples = 10_000,
        nsta = 20,
        lon = 5rand(nsta) .- 2.5,
        lat = 5rand(nsta) .- 2.5,
        distance = nothing,
        noise = 0.0
        )
    # Synthetic arrival
    tmax = nsamples*0.01/2
    data = exp.(-(-tmax:0.01:tmax).^2 ./ 0.01)
    maxamp = maximum(data)
    # Add noise
    if noise > 0
        data .+= 2maxamp*noise.*randn.()
    end
    # Grid of stations
    s = [Seis.Trace(0, delta, data) for _ in 1:nsta]
    s.sta.lon = lon
    s.sta.lat = lat
    s.sta.sta = string.(1:nsta)
    mlon, mlat = mean(lon), mean(lat)

    # Add event if requested
    evt = Seis.Event{eltype(s[1])}()
    if distance !== nothing
        β = atand(sx, sy)
        evt.lon, evt.lat, _ = Geodesics.angular_step(mlon, mlat, β,
            distance, true)
        evt.dep = 0
        evt.id = "Imaginary event"
    end
    s.evt = evt

    # Shift arrivals to correct time
    for (i, ss) in enumerate(s)
        dist = Beamforming.delta(ss.sta.lon, ss.sta.lat, mlon, mlat, true)
        az = Beamforming.azimuth(mlon, mlat, ss.sta.lon, ss.sta.lat, true)
        dx = dist*sind(az) # °
        dy = dist*cosd(az)
        dt = sx*dx + sy*dy # s
        ss.t = circshift(ss.t, -round(Int, dt/delta))
    end
    s
end

end # module

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
push!(LOAD_PATH,"../")

include("testing_dipole_response.jl")
using Plots

"""
Atomic units:
-------------------------
kb = 1
hbar = 1
4*pi*e0 = 1
a_0 = 1
e = 1
-------------------------
distance = bohr (a_0)
energy = hartree (27.211385 eV)
time = hbar/E_h (2.418884326505(16)×10−17 s)
dipole moment = ea_0 (8.47835326(19)×10−30 C·m = 2.541746 Debye)
polarisability = e2 a_02/E_h (1.65E-41 C^2 m^2 /J)
temperature = E_h/kB (3.1577464E5 K)
"""

#'''SetUp'''
const kB=1 #8.6173324E-5 #eV
const hbar=1.0

const N=50 # Number of sites in model
             # --> Size of tridiagonal Hamiltonian constructed; larger value -> better statistics, slower runtime
const Edisorder=0.0 # Energetic disorder eV, Gaussian form
const Jdisorder=0.0 # Transfer integral disorder, eV.
const r = 20

# Model setup
const J0=0.00001 #(~0.1eV)
modelJ(θ) = J0*cos(θ*π/180.0).^2

const T=0.0032 #(~300K)
const B=1/(T*kB) #300K * k_B in eV

# This effectively reduces down to the 'alpha' parameter in the Frohlich polaron Hamiltonian
const dipolestrength=0.2

function TestNewDipole(dipolefield::Bool = false)
    #'''prepare'''
    maxS=[]
    maxS2 = []

    S,E,H,psi,density = prepare_model2()
    dipole = zeros(N)
    Field = FieldFromDensity(density)
    if dipolefield
        FieldDI = FieldFromDipole(dipole)
        Field = Field+FieldDI
    end
    dipole = UpdateDipole(Field,density,dipole)
    S,E,H = UpdateEnergy2(dipole,Field,S,E)

    for i in 1:100
        S,E,H = UpdateEnergy2(dipole,Field,S,E)
        psi = eigvecs(H)[:,1]
        density = psi.^2
        Field = FieldFromDensity(density)
        if dipolefield
            FieldDI = FieldFromDipole(dipole)
            Field = Field+FieldDI
            norm_FieldDI = norm(FieldDI)
        end
        dipole = UpdateDipole(Field,density,dipole)
        norm_dipole = norm(dipole)
        norm_Field = norm(Field)
        norm_S = norm(S)
        norm_density=norm(density)
        xlims!(1,N)
        ylims!(-1,1)
        plot(real(dipole/norm_dipole), label = "dipole")
        plot!(Field, label="Field")
        plot!(real(density)/norm_density, label = "density")
        plot!(real(S/norm_S), label = "energy")
        gui()
    end

    psi=psi+nondispersive_wavepacket(2,4.0)
    normalise = norm(psi)
    psi = psi/normalise
    density = conj(psi).*psi
    #plot(real(density))
    #gui()
    Field = FieldFromDensity(density)
    if dipolefield
        FieldDI = FieldFromDipole(dipole)
        Field = Field+FieldDI
        norm_FieldDI = norm(FieldDI)
    end
    dipole = UpdateDipole(Field,density,dipole)
    S,E,H = UpdateEnergy2(dipole,Field,S,E)

    #'''propagate'''
    for i in 1:500
        norm_dipole = norm(dipole)
        plot(real(dipole/norm_dipole), label = "dipole")
        norm_Field = norm(Field)
        norm_S = norm(S)
        #plot!(imag(psi))
        #plot!(real(psi))
        xlims!(1,N)
        #ylims!(-1,1)
        norm_density=norm(density)
        plot!(real(density)/norm_density, label = "density")
        #plot!(FieldDI/norm_FieldDI, label = "FieldDI")
        #plot!(Field/norm_Field, label = "Field")
        #plot!((Field+FieldDI)/(norm_Field+norm_FieldDI), label = "Fields")
        plot!(real(S/norm_S), label = "energy")
        gui()
        psi, density = Propagate(H,psi,1000)
        Field = FieldFromDensity(real(density))
        if dipolefield
            FieldDI = FieldFromDipole(dipole)
            Field = Field+FieldDI
            norm_FieldDI = norm(FieldDI)
        end
        dipole = UpdateDipole(Field, density,dipole)
        S,E,H = UpdateEnergy2(dipole,Field,S,E)
        append!(maxS, indmin(S))
        Snew = S
        filter!(x->x≠minimum(Snew),Snew)
        append!(maxS2, indmin(Snew))

    end
    return maxS,maxS2
end

maxNoDipole, maxNoDipole2 = TestNewDipole()
maxDipole, maxDipole2 = TestNewDipole(true)

plot(maxNoDipole, linestyle = :dot, label = "NoDipole")
plot!(maxDipole, linestyle = :solid, label = "Dipole")
plot!(maxNoDipole2, linestyle = :dot, label = "NoDipole2ndPeak")
plot!(maxDipole2, linestyle = :dot, label = "Dipole2ndPeak")


text:
module ComputExp

using JuliaDB
using Missings
using Query
using NamedTuples

include("utils.jl")
include("types.jl")
include("api.jl")
include("run_exp.jl")

export  ComputExpDB,
        New!,
        Language,
        Algorithm,
        Implementation,
        Problem,
        Experiment,
        Run,
        @run_exps

end # module


text:
# YAML Tokens.
# Each token must include at minimum member "span::Span".
abstract Token

# The '%YAML' directive.
immutable DirectiveToken <: Token
    span::Span
    name::AbstractString
    value::(@compat Union{Tuple, Void})
end

# '---'
immutable DocumentStartToken <: Token
    span::Span
end

# '...'
immutable DocumentEndToken <: Token
    span::Span
end

# The stream start
immutable StreamStartToken <: Token
    span::Span
    encoding::AbstractString
end

# The stream end
immutable StreamEndToken <: Token
    span::Span
end

#
immutable BlockSequenceStartToken <: Token
    span::Span
end

#
immutable BlockMappingStartToken <: Token
    span::Span
end

#
immutable BlockEndToken <: Token
    span::Span
end

# '['
immutable FlowSequenceStartToken <: Token
    span::Span
end

# '{'
immutable FlowMappingStartToken <: Token
    span::Span
end

# ']'
immutable FlowSequenceEndToken <: Token
    span::Span
end

# '}'
immutable FlowMappingEndToken <: Token
    span::Span
end

# '?' or nothing (simple keys).
immutable KeyToken <: Token
    span::Span
end

# ':'
immutable ValueToken <: Token
    span::Span
end

# '-'
immutable BlockEntryToken <: Token
    span::Span
end

# ','
immutable FlowEntryToken <: Token
    span::Span
end

# '*anchor'
immutable AliasToken <: Token
    span::Span
    value::AbstractString
end

# '&anchor'
immutable AnchorToken <: Token
    span::Span
    value::AbstractString
end

# '!handle!suffix'
immutable TagToken <: Token
    span::Span
    value
end

# A scalar.
immutable ScalarToken <: Token
    span::Span
    value::AbstractString
    plain::Bool
    style::(@compat Union{Char, Void})
end


text:
using Test
using SymbolicUtils

# == / != syntax is nice, let's use it in tests
macro eqtest(expr)
    @assert expr.head == :call && expr.args[1] in [:(==), :(!=)]
    if expr.args[1] == :(==)
        :(@test isequal($(expr.args[2]), $(expr.args[3])))
    else
        :(@test !isequal($(expr.args[2]), $(expr.args[3])))
    end |> esc
end
SymbolicUtils.show_simplified[] = false

#using SymbolicUtils: Rule

include("basics.jl")
include("order.jl")
include("rewrite.jl")
include("rulesets.jl")
include("interface.jl")
include("fuzz.jl")
if haskey(ENV, "TRAVIS")
    include("benchmark.jl")
end


text:
module LogisticRegression
using LinearAlgebra
using Distributions
using Random
# using UnicodePlots
using SparseArrays
BLAS.set_num_threads(2)

"""
This structure implements the Sklearn.jl api for
the Logistic Regression algorithm
"""
mutable struct LogisticRegressor
    num_steps::Int64
    cutoff::Float64
    learning_rate::Float64
    add_intercept::Bool
    weights::Array{Float64, 1}
end
LogisticRegressor(;
    num_steps=300000,
    cutoff=1e-4,
    learning_rate=5e-5,
    add_intercept=false,
    weights=Float64[]
    ) = LogisticRegressor(  num_steps,
                            cutoff,
                            learning_rate,
                            add_intercept,
                            weights )

function fit!(self::LogisticRegressor, data, labels::Array{T, 1}) where T <: Number
    weights = logistic_regression(  data,
                                    labels,
                                    num_steps = self.num_steps,
                                    cutoff = self.cutoff,
                                    learning_rate = self.learning_rate,
                                    add_intercept = self.add_intercept)
    self.weights = weights
end
function predict(self::LogisticRegressor, data)
    if self.add_intercept
        data = [ones(size(data)[1]) data]
    end
    final_scores = data * self.weights
    preds = round.(par_sigmoid(final_scores))
    return preds
end
function sigmoid(score::T) where T <: Number
    return 1/(1 + exp(-score))
end
function par_sigmoid(scores::Array{T}) where T <: Number
    res = zeros(length(scores))
    Threads.@threads for i in eachindex(scores)
        @inbounds res[i] = sigmoid(scores[i])
    end
    return res
end
function log_likelihood(feats, label, weights)
    scores = feats * weights
    log_likely = sum(label .* scores - log.(1 .+ exp.(scores)))
    return log_likely
end
function logistic_regression(feats, labels; num_steps=20000, cutoff=1e-1000, learning_rate=1e-2, add_intercept=false)
    if add_intercept == true
        intercept = ones(size(feats)[1])
        feats = [intercept feats]
    end
    # should have a weight for each column
    weights = zeros(size(feats)[2])
    scores = zeros(size(feats)[2])
    error_signal = zeros(size(feats)[2])
    last_log_likelihood_delta = 100.0
    current_log_likelihood_delta = 100.0
    last_log_likelihood = 10000.0
    current_log_likelihood = 10000.0
    for step=1:num_steps
        scores = feats * weights
        predictions = par_sigmoid(scores)

        # update using the gradient
        error_signal = labels - predictions
        gradient = transpose(feats) * error_signal
        weights += learning_rate * gradient
        if step % 10000 == 0
            current_log_likelihood = log_likelihood( feats,
                                                     labels,
                                                     weights)
            current_log_likelihood_delta = abs(current_log_likelihood - last_log_likelihood)
            println("Change in log likelihood: " * string(current_log_likelihood_delta))
            if current_log_likelihood_delta < cutoff && last_log_likelihood_delta < cutoff
                break # stop training if we haven't exceeded cutoff in 20000 iterations
            end
            last_log_likelihood_delta = current_log_likelihood_delta
            last_log_likelihood = current_log_likelihood
        end
    end
    return weights
end
function test_log_reg()
    num_instances = 5000
    # generate dense training data
    #x1_v1 = rand(Normal(0, 2), num_instances)
    #x1_v2 = rand(Normal(0, .3), num_instances)
    #x2_v1 = rand(Normal(4, 1), num_instances)
    #x2_v2 = rand(Normal(2, .3), num_instances)

    # generate sparse training data
    x1_v1 = sparse(rand(Binomial(4, .2), num_instances))
    x1_v2 = sparse(rand(Binomial(8, .1),
    num_instances))
    x2_v1 = sparse(rand(Binomial(4, .8), num_instances))
    x2_v2 = sparse(rand(Binomial(8, .7), num_instances))

    x1 = [x1_v1 x1_v2]
    x2 = [x2_v1 x2_v2]
    # generate sparse test data
    x3_v1 = sparse(rand(Binomial(4, .2), num_instances))
    x3_v2 = sparse(rand(Binomial(8, .1),
    num_instances))
    x4_v1 = sparse(rand(Binomial(4, .8), num_instances))
    x4_v2 = sparse(rand(Binomial(8, .7), num_instances))
    x3 = [x3_v1 x3_v2]
    x4 = [x4_v1 x4_v2]
    test_data = [x3 ; x4]
    test_labels = [zeros(num_instances) ; ones(num_instances)]
    # display the dataset we just made
    # plt = scatterplot(x1[:, 1], x2[:, 2], ylim=[-2,3])
    # scatterplot!(plt, x2[:, 1], x2[:, 2], color=:blue)
    # combine these two separate
    # normal distributions
    data = [x1 ; x2]
    # data_with_intercept = [ones(size(data)[1]) data]
    labels = [zeros(num_instances) ; ones(num_instances)]
    # weights = logistic_regression(data, labels, num_steps=300000, learning_rate=5e-5, add_intercept=true)
    # final_scores = data_with_intercept * weights
    # preds = round.(par_sigmoid(final_scores))
    model = LogisticRegressor(add_intercept=true)
    fit!(model, data, labels)
    test_preds = predict(model, test_data)
    # println("Accuracy on training data: ")
    # acc = sum(preds .== labels)/length(labels)
    # println(acc)
    # test_data_with_intercept = [ones(size(test_data)[1]) test_data]
    # test_scores = test_data_with_intercept * weights
    # test_preds = round.(par_sigmoid(test_scores))
    println("Accuracy on test data: ")
    acc = sum(test_preds .== test_labels)/length(test_labels)
    println(acc)
    return acc
end
export LogisticRegressor, fit!, predict
end # module


text:
# This file was generated, do not modify it. # hide
mutable struct MyNetworkBuilder <: MLJFlux.Builder
    n1::Int #Number of cells in the first hidden layer
    n2::Int #Number of cells in the second hidden layer
end


text:
@testset "parser" begin
    function test_str_to_formula_to_str(formula::String)
        @customtest String(Formula(formula)) == formula
    end
    test_str_to_formula_to_str("(P(f(x), g(y)) & G(f(x)))")
    test_str_to_formula_to_str("~*{x}Foo(x)")
    test_str_to_formula_to_str("~~(*{y}Foo(x) | Bar(f(x), g(f(y), f(h, u))))")
    test_str_to_formula_to_str("(~?{y}*{g}P(f(y)) & (?{x}P(x) | *{y}P(y)))")
end


text:
# test the StataCall() functions

using DataFrames, Test
using Random
using Missings

myeps = 1e-5

Random.seed!(1)
df = DataFrame(myint = Int64.(floor.(100 .*rand(Float64, 10))), myfloat = rand(Float64, 10))
instructions = ["gen newvar1 = myint + myfloat";
"gen newvar2 = floor(_n/2)";
"bysort newvar2: egen newvar3 = mean(newvar1)"
]
dfOut = StataCall.stataCall(instructions, df)
testOut = [55.236034f0
42.829613
42.829613
52.74826
52.74826
23.081442
23.081442
15.625783
15.625783
85.986664]
for i in 1:length(testOut)
    @test isapprox(dfOut[!,:newvar3][i], testOut[i], atol = myeps)
end

# With missing values
Random.seed!(1)
df = DataFrame(myint = Int64.(floor.(100 .*rand(Float64, 10))), myfloat = rand(Float64, 10))
allowmissing!(df)
df[4,:myint] = missing
df[2,:myfloat] = missing

instructions = ["gen newvar1 = myint + myfloat";
"gen newvar2 = floor(_n/2)";
"bysort newvar2: egen newvar3 = mean(newvar1)"
]
dfOut = StataCall.stataCall(instructions, df)
testOut = [55.236034
42.312706
42.312706
28.488613
28.488613
23.081442
23.081442
15.625783
15.625783
85.986664]
for i in 1:length(testOut)
    @test isapprox(dfOut[!,:newvar3][i], testOut[i], atol = myeps)
end


text:
using Gadfly

set_default_plot_size(6.6inch, 3.3inch)

p1 = plot(x=[], y=[], Geom.point)
p2 = plot(x=[], y=[], Geom.line,
    Guide.xticks(ticks=[0,10]), Guide.yticks(ticks=[0,10])
)

hstack(p1,p2)


text:
module aoc_03
    using Test
    export read_input, part1, part2
    export test_part1, test_part2, test
    export input, demo

    const input_r = r""

    function read_input(filename = "input.txt")
        lines = readlines(filename)
        # for line in eachline(filename)
        input = [parse_line(line) for line in eachline(filename)]
        permutedims(hcat(input...))
    end

    function parse_line(line)
        Vector{Char}(line) .== '#'
    end

    const input = read_input("input.txt")
    const demo = read_input("demo.txt")

    function part1(input = input, slope = CartesianIndex(1,3))
        pos = CartesianIndex(1,1)
        N = size(input,1)
        M = size(input,2)
        trees = 0
        while pos[1] <= N
            trees += input[ pos ]
            pos += slope
            pos = CartesianIndex( pos[1], mod1( pos[2] , M))
        end
        trees
    end

    function test_part1()
        @test part1( demo ) == 7
        @test part1() == 228
    end

    function part2(input = input)
# Right 1, down 1.
# Right 3, down 1. (This is the slope you already checked.)
# Right 5, down 1.
# Right 7, down 1.
# Right 1, down 2.
        slopes = [
            CartesianIndex(1,1),
            CartesianIndex(1,3),
            CartesianIndex(1,5),
            CartesianIndex(1,7),
            CartesianIndex(2,1)
        ]
        prod( part1.( (input,), slopes ) )
    end

    function test_part2()
        @test part2( demo ) == 336
    end

    function test()
        test_part1()
        test_part2()
    end
end


text:
export RecipeInflation, filter_state!, vortexassim

struct RecipeInflation <: InflationType
    "Parameters"
    p::Array{Float64,1}
end

# Filtering function to bound the strength of the vortices and sources
function filter_state!(x, config::VortexConfig)
	@inbounds for j=1:config.Nv
		# Ensure that vortices stay above the x axis
		x[3*j-1] = clamp(x[3*j-1], 1e-2, Inf)
		# Ensure that the circulation remains positive
    	x[3*j] = clamp(x[3*j], 0.0, Inf)
	end
    return x
end

# This function apply additive inflation to the state components only,
# not the measurements, X is an Array{Float64,2} or a view of it
function (ϵ::RecipeInflation)(X, Ny, Nx, config::VortexConfig)
	ϵX, ϵΓ = ϵ.p
	Nv = config.Nv
	@assert Nx == 3*Nv
	for col in eachcol(X)
		if config.Nv > 0
			for i in 1:Nv
				col[Ny + 3*(i-1) + 1: Ny + 3*(i-1) + 2] .+= ϵX*randn(2)
				col[Ny + 3*i] += ϵΓ*randn()
			end
		end
	end
end

function (ϵ::RecipeInflation)(x::AbstractVector{Float64}, config::VortexConfig)
	ϵX, ϵΓ = ϵ.p
	Nv = config.Nv
	Nx = size(x, 1)

	@assert Nx == 3*Nv
	for i=1:Nv
		x[3*(i-1) + 1:3*(i-1) + 2] .+= ϵX*randn(2)
		x[3*(i-1) + 3] += ϵΓ*randn()
	end
end

# Create a function to perform the sequential assimilation for any sequential filter SeqFilter
function vortexassim(algo::SeqFilter, X, tspan::Tuple{S,S}, config::VortexConfig, data::SyntheticData; withfreestream::Bool = false, P::Parallel = serial) where {S<:Real}

	# Define the additive Inflation
	ϵX = config.ϵX
	ϵΓ = config.ϵΓ
	β = config.β
	ϵY = config.ϵY

	Ny = size(config.ss,1)

	ϵx = RecipeInflation([ϵX; ϵΓ])
	ϵmul = MultiplicativeInflation(β)
	# Set different times
	Δtobs = algo.Δtobs
	Δtdyn = algo.Δtdyn
	t0, tf = tspan
	step = ceil(Int, Δtobs/Δtdyn)

	n0 = ceil(Int64, t0/Δtobs) + 1
	J = (tf-t0)/Δtobs
	Acycle = n0:n0+J-1

	# Array dimensions
	Nypx, Ne = size(X)
	Nx = Nypx - Ny
	ystar = zeros(Ny)

	cachevels = allocate_velocity(state_to_lagrange(X[Ny+1:Ny+Nx,1], config))

	h(x, t) = measure_state(x, t, config; withfreestream = withfreestream)
	press_itp = CubicSplineInterpolation((LinRange(real(config.ss[1]), real(config.ss[end]), length(config.ss)),
	                                   t0:data.Δt:tf), data.yt, extrapolation_bc =  Line())

	yt(t) = press_itp(real.(config.ss), t)
	Xf = Array{Float64,2}[]
	push!(Xf, copy(state(X, Ny, Nx)))

	Xa = Array{Float64,2}[]
	push!(Xa, copy(state(X, Ny, Nx)))

	# Run particle filter
	@showprogress for i=1:length(Acycle)
	   # Forecast step
	   @inbounds for j=1:step
		   tj = t0+(i-1)*Δtobs+(j-1)*Δtdyn
		   X, _ = vortex(X, tj, Ny, Nx, cachevels, config, withfreestream = withfreestream)
	   end

	   push!(Xf, deepcopy(state(X, Ny, Nx)))

	   # Get real measurement
	   ystar .= yt(t0+i*Δtobs)

	   # Perform state inflation
	   ϵmul(X, Ny+1, Ny+Nx)
	   ϵx(X, Ny, Nx, config)

	   # Filter state
	   if algo.isfiltered == true
		   @inbounds for i=1:Ne
			   x = view(X, Ny+1:Ny+Nx, i)
			   x .= filter_state!(x, config)
		   end
	   end

	   observe(h, X, t0+i*Δtobs, Ny, Nx; P = P)

	   ϵ = algo.ϵy.σ*randn(Ny, Ne) .+ algo.ϵy.m

	   Xpert = (1/sqrt(Ne-1))*(X[Ny+1:Ny+Nx,:] .- mean(X[Ny+1:Ny+Nx,:]; dims = 2)[:,1])
	   HXpert = (1/sqrt(Ne-1))*(X[1:Ny,:] .- mean(X[1:Ny,:]; dims = 2)[:,1])
	   ϵpert = (1/sqrt(Ne-1))*(ϵ .- mean(ϵ; dims = 2)[:,1])
	   # Kenkf = Xpert*HXpert'*inv(HXpert*HXpert'+ϵpert*ϵpert')

	   b = (HXpert*HXpert' + ϵpert*ϵpert')\(ystar .- (X[1:Ny,:] + ϵ))
	   view(X,Ny+1:Ny+Nx,:) .+= (Xpert*HXpert')*b

	   # @show cumsum(svd((Xpert*HXpert')*inv((HXpert*HXpert' + ϵpert*ϵpert'))).S)./sum(svd((Xpert*HXpert')*inv((HXpert*HXpert' + ϵpert*ϵpert'))).S)

	   # X = algo(X, ystar, t0+i*Δtobs)

	   # Filter state
	   if algo.isfiltered == true
		   @inbounds for i=1:Ne
			   x = view(X, Ny+1:Ny+Nx, i)
			   x .= filter_state!(x, config)
		   end
	   end

	   push!(Xa, deepcopy(state(X, Ny, Nx)))
	end

	return Xf, Xa
end


text:
@with_kw struct IZParameter_more
    a::SNNFloat = 0.01
    b::SNNFloat = 0.2
    c::SNNFloat = -65
    d::SNNFloat = 2
    C::SNNFloat = 0.9
    vr::SNNFloat = -65
    k::SNNFloat = 1.6
    vPeak::SNNFloat = 20
    vt::SNNFloat = -50
    cellnumber::SNNInt = 1
end

@with_kw mutable struct IZ_more
    param::IZParameter_more = IZParameter_more(a, b, c, d, C, vr, k, vPeak, vt,cellnumber)
    N::SNNInt = 1
    v::Vector{SNNFloat} = fill(param.vr, N)
    u::Vector{SNNFloat} = param.b * v
    fire::Vector{Bool} = zeros(Bool, N)
    I::Vector{SNNFloat} = zeros(N)
    records::Dict = Dict()
end

function integrate!(p::IZ_more, param::IZParameter_more, dt::SNNFloat)
    @unpack a, b, c, d, C, vr, k, vPeak, vt,cellnumber = param
    if cellnumber>=1 && cellnumber<=3
        integrate_one_three!(p, param, SNNFloat(dt))
    elseif cellnumber == 4
        integrate_four!(p,param, SNNFloat(dt))
    elseif cellnumber == 5
        integrate_five!(p, p.param, SNNFloat(dt))
    elseif cellnumber == 6
        integrate_six!(p, p.param, SNNFloat(dt))
    elseif cellnumber == 7
        integrate_seven!(p, p.param, SNNFloat(dt))
    end
end
function integrate_one_three!(p::IZ_more, param::IZParameter_more, dt::SNNFloat)
    @unpack N, v, u, fire, I = p
    @unpack a, b, c, d, C, vr, k, vPeak, vt, cellnumber = param
    @inbounds for i = 1:N
        v[i] = v[i] + 0.25 * (k * (v[i] - vr) * (v[i] - vt) - u[i] + I[i]) / C
        u[i] = u[i] + 0.25*a*(b*(v[i]-vr)-u[i]) # Calculate recovery variable
        fire[i] = v[i] > vPeak

        if v[i]>=vPeak
            v[i]=c
            u[i]=u[i]+d  # reset u, except for FS cells
        end
        v[i] = ifelse(fire[i], c, v[i])
        u[i] = ifelse(fire[i],d, u[i])
    end
end

function integrate_four!(p::IZ_more, param::IZParameter_more, dt::SNNFloat)
    @unpack N, v, u, fire, I = p
    @unpack a, b, c, d, C, vr, k, vPeak, vt, cellnumber = param
    @inbounds for i = 1:N

        v[i] = v[i] + 0.25 * (k * (v[i] - vr) * (v[i] - vt) - u[i] + I[i]) / C
        u[i] = u[i] + 0.25*a*(b*(v[i]-vr)-u[i]) # Calculate recovery variable
        fire[i] = v[i] > vPeak

        if v[i]>=vPeak- 0.1*u[i]
            v[i] = vPeak - 0.1*u[i]
            v[i] = c + 0.04*u[i]; # Reset voltage
            if (u[i]+d)<670
                u[i] = u[i]+d; # Reset recovery variable
            else
                u[i] = 670;
            end
        end
    end
    #v[i] = ifelse(fire[i], c, v[i])
    #u[i] = ifelse(fire[i],d, u[i])
end

function integrate_five!(p::IZ_more, param::IZParameter_more, dt::SNNFloat)
    @unpack N, v, u, fire, I = p
    @unpack a, b, c, d, C, vr, k, vPeak, vt, cellnumber = param
    tau::SNNFloat = 0.25
    @inbounds for i = 1:N
        v[i] = v[i] + 0.25 * (k * (v[i] - vr) * (v[i] - vt) - u[i] + I[i]) / C

        #v[i] = v[i] + tau * (k * (v[i] - vr) * (v[i] - vt) - u[i] + I) / C

        #u[i+1]=u[i]+tau*a*(b*(v[i]-vr)-u[i]); # Calculate recovery variable
        if v[i] < d
            u[i] = u[i] + tau*a*(0-u[i])
        else
            u[i] = u[i] + tau*a*((0.025*(v[i]-d)^3)-u[i])
        end
        if v[i]>=vPeak
            v[i]=vPeak;
            v[i]=c;
        end
        fire[i] = v[i] > vPeak
    end
end
function integrate_six!(p::IZ_more, param::IZParameter_more, dt::SNNFloat)
    @unpack N, v, u, fire, I = p
    @unpack a, b, c, d, C, vr, k, vPeak, vt, cellnumber = param
    tau::SNNFloat = 0.25
    @inbounds for i = 1:N

        v[i] = v[i] + tau * (k * (v[i] - vr) * (v[i] - vt) - u[i] + I[i]) / C
        if v[i] > -65
            b=0;
        else
            b=15;
        end
        u[i] = u[i]+tau*a*(b*(v[i]-vr)-u[i]);

        if v[i] > (vPeak + 0.1*u[i])
            v[i] = vPeak + 0.1*u[i];
            v[i] = c-0.1*u[i]; # Reset voltage
            u[i]=u[i]+d;
        end
        v[i] = ifelse(fire[i], c, v[i])
        u[i] = ifelse(fire[i],d, u[i])
    end
end
function integrate_seven!(p::IZ_more, param::IZParameter_more, dt::SNNFloat)
    @unpack N, v, u, fire, I = p
    @unpack a, b, c, d, C, vr, k, vPeak, vt, cellnumber = param
    tau::SNNFloat = 0.25
    @inbounds for i = 1:N

        v[i] = v[i] + 0.25 * (k * (v[i] - vr) * (v[i] - vt) - u[i] + I[i]) / C

        if v[i] > -65
            b=2;
        else
            b=10;
        end
        u[i]=u[i]+0.25*a*(b*(v[i]-vr)-u[i]);
        if v[i]>=vPeak
            v[i]=vPeak;
            v[i]=c;
            u[i]=u[i]+d;  # reset u, except for FS cells
        end
        v[i] = ifelse(fire[i], c, v[i])
        u[i] = ifelse(fire[i],d, u[i])
    end
end

#=
function integrate!(p::IZ_more, param::IZParameter_more, dt::SNNFloat)
    @unpack N, v, u, fire, I = p
    @unpack a, b, c, d, C, vr, k, vt, vPeak = param
    #tau = 0.25f0; #dt
    @inbounds for i = 1:N
        println(i," i is 1?")
        #v[i] += 0.5f0dt * (0.04f0v[i]^2 + 5f0v[i] + 140f0 - u[i] + I[i])
        v[i] += v[i] + 0.25f0 * (k * (v[i] - vr) * (v[i] - vt) - u[i] + I[i]) / C
        u[i] += u[i] + 0.25f0*a*(b*(v[i]-vr)-u[i]) # Calculate recovery variable
        #v[i+1] = v[i] + 0.25f0 * (k * (v[i] - vr) * (v[i] - vt) - u[i] + I[i]) / C
        #u[i+1] = u[i] + 0.25f0*a*(b*(v[i]-vr)-u[i]); # Calculate recovery variable
        println(v[i],u[i],C," ",I[i]," ",i)

        #fire[i] = v[i] > vPeak
        if v[i+1]>=vPeak
            v[i]=vPeak
            v[i+1]=c
            u[i+1]=u[i+1]+d  # reset u, except for FS cells
        end
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > vPeak
        v[i] = ifelse(fire[i], c, v[i])
        u[i] += ifelse(fire[i], d, 0f0)
    end
end
=#

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


text:
module SweepOperator
export sweep!

using LinearAlgebra
import LinearAlgebra: BlasFloat, checksquare

const AMat = AbstractMatrix
const AVec = AbstractVector

"""
    sweep!(A, k ; inv=false)
    sweep!(A, ks; inv=false)

Perform the sweep operation (or inverse sweep if `inv=true`) on matrix `A` on element `k`
(or each element in `ks`).  Only the upper triangle is read/swept.

# Example:

    x = randn(100, 10)
    xtx = x'x
    sweep!(xtx, 1)
    sweep!(xtx, 1, true)
"""
function sweep!(A::AMat, k::Integer, inv::Bool = false)
    sweep_with_buffer!(Vector{eltype(A)}(undef, size(A, 2)), A, k, inv)
end

function sweep!(A::AMat{T}, ks::AVec{I}, inv::Bool = false) where {T<:BlasFloat, I<:Integer}
    akk = Vector{T}(undef, size(A,1))
    for k in ks
        sweep_with_buffer!(akk, A, k, inv)
    end
    A
end

function sweep_with_buffer!(akk::AVec{T}, A::AMat{T}, k::Integer, inv::Bool = false) where {T}
    # ensure @inbounds is safe
    p = checksquare(A)
    1 ≤ k ≤ p || throw(BoundsError(A, k))
    p == length(akk) || throw(DimensionError("Incorrect buffer size."))
    @inbounds @views begin
        d = one(T) / A[k, k]                # pivot
        copy!(akk, Symmetric(A, :U)[:, k])  # akk = A[:, k]
        syrk!(A, -d, akk)                   # everything not in row/col k
        rmul!(akk, d * (-one(T)) ^ inv)     # akk .* d (negated if inv=true)
        setrowcol!(A, k, akk)
        A[k, k] = -d                        # pivot element
    end
    return A
end

#-----------------------------------------------------------------------------# setrowcol!
# Set upper triangle of: (A[k, :] = x; A[:, k] = x)
function setrowcol!(A::StridedArray, k, x)
    @views copy!(A[1:k-1,k], x[1:k-1])         # col k
    @views copy!(A[k, k+1:end], x[k+1:end])    # row k
end

setrowcol!(A::Union{Hermitian,Symmetric,UpperTriangular}, k, x) = setrowcol!(A.data, k, x)

#-----------------------------------------------------------------------------# syrk!
const BlasNumber = Union{LinearAlgebra.BlasFloat, LinearAlgebra.BlasComplex}

# In-place update of (the upper triangle of) A + α * x * x'
function syrk!(A::StridedMatrix{T}, α::T, x::AbstractArray{<:T}) where {T<:BlasNumber}
    BLAS.syrk!('U', 'N', α, x, one(T), A)
end

function syrk!(A::Hermitian{T, S}, α::T, x::AbstractArray{<:T}) where {T<:BlasNumber, S<:StridedMatrix{T}}
    Hermitian(BLAS.syrk!('U', 'N', α, x, one(T), A.data))
end

function syrk!(A::Symmetric{T, S}, α::T, x::AbstractArray{<:T}) where {T<:BlasNumber, S<:StridedMatrix{T}}
    Symmetric(BLAS.syrk!('U', 'N', α, x, one(T), A.data))
end

function syrk!(A, α, x) where {T}
    p = checksquare(A)
    for i in 1:p, j in i:p
        @inbounds A[i,j] += α * x[i] * x[j]
    end
end

end # module


text:
function op_K( pw_grid::PWGrid, psi::Array{ComplexF64,2} )
    out = zeros(ComplexF64,size(psi))
    Ncol = size(psi,2)
    Ω = pw_grid.Ω
    G2 = pw_grid.G2
    Npoints = pw_grid.Npoints
    for is = 1:Ncol
        for ip = 1:Npoints
            out[ip,is] = psi[ip,is]*G2[ip]
        end
    end
    return 0.5*out # two minus signs
end

function op_K( pw_grid::PWGrid, psi::Array{ComplexF64,1} )
    G2 = pw_grid.G2
    Npoints = pw_grid.Npoints
    out = zeros( size(psi) )
    for ip = 1:Npoints
        out[ip] = psi[ip]*G2[ip]
    end
    return 0.5*out
end