text:
using Plots; gr()
using Measures
using FFTW
push!(LOAD_PATH, pwd())
using DHC_2DUtils

function plot_filter_bank_QA(filt, info; fname="filter_bank_QA.png")

    # -------- define plot1 to append plot to a list
    function plot1(ps, image; clim=nothing, bin=1.0, fsz=12, label=nothing)
        # ps    - list of plots to append to
        # image - image to heatmap
        # clim  - tuple of color limits (min, max)
        # bin   - bin down by this factor, centered on nx/2+1
        # fsz   - font size
        marg   = -1mm
        nx, ny = size(image)
        nxb    = nx/round(Integer, 2*bin)

        # -------- center on nx/2+1
        i0 = max(1,round(Integer, (nx/2+2)-nxb-1))
        i1 = min(nx,round(Integer, (nx/2)+nxb+1))
        lims = [i0,i1]
        subim = image[i0:i1,i0:i1]
        push!(ps, heatmap(image,aspect_ratio=:equal,clim=clim,
            xlims=lims, ylims=lims, size=(400,400),
            legend=false,xtickfontsize=fsz,ytickfontsize=fsz,#tick_direction=:out,
            rightmargin=marg,leftmargin=marg,topmargin=marg,bottommargin=marg))
        if label != nothing
            annotate!(lims'*[.96,.04],lims'*[.09,.91],text(label,:left,:white,16))
        end
        return
    end

    # -------- initialize array of plots
    ps   = []
    ind  = info["filter_index"]
    jval = info["j_value"]
    J    = length(jval)
    L    = length(info["theta_value"])
    # plot the first 3 and last 3 j values
    jplt = [1:3;J-2:J]

    # -------- loop over j
    for j_ind in jplt
        # first 3 filters in Fourier space
        j = jval[j_ind]
        bfac = max(1,2.0^(j-1))
        for l = 1:3
            label=string("j=",j_ind," ℓ=",l-1)
            #label=latexstring("\\tilde{\\psi}_{$j_ind,$(l-1)}")
            plot1(ps, fftshift(filt[:,:,ind[j_ind,l]]), clim=(0,1),
                  bin=bfac, label=label)
        end
        fac = max(1,2.0^(5-j))
        # real part of config space of 3rd filter
        plot1(ps, fftshift(real(ifft(filt[:,:,ind[j_ind,3]]))),
              bin=fac, label=string("j=",j_ind," ℓ=",3-1))

        # sum all filters to form (partial) ring
        ring = dropdims(sum(filt[:,:,ind[j_ind,:]].^2, dims=3), dims=3)
        plot1(ps, fftshift(ring[:,:,1]), clim=(0,1),
              bin=bfac, label=string("j=",j_ind," ℓ=",0,":",L-1))
    end

    wavepow = fftshift(dropdims(sum(filt[:,:,ind].^2, dims=(3,4)), dims=(3,4)))
    plot1(ps, wavepow, clim=(-0.1,1), label=string("j=1:",J," ℓ=0:",L-1))

    if info["pc"]==1
        wavepow += circshift(wavepow[end:-1:1,end:-1:1],(1,1))
    end
    plot1(ps, wavepow, clim=(-0.1,1), bin=32, label="missing power")
    phi = filt[:,:,info["phi_index"]]
    phi_shift = fftshift(phi)
    plot1(ps, phi_shift.^2, clim=(-0.1,1), bin=32, label="ϕ")
    disc = wavepow+(phi_shift.^2)
    plot1(ps, fftshift(real(ifft(phi))), label="ϕ")
    plot1(ps, disc, clim=(-0.1,1), bin=32, label="all")

    myplot = plot(ps..., layout=(7,5), size=(1400,2000))
    savefig(myplot, fname)
end

# -------- define a filter bank

filt, info  = fink_filter_bank(1, 8, pc=1, wd=1)
hash = fink_filter_hash(1, 8, pc=1, wd=1)
plot_filter_bank_QA(filt, info, fname="filt-8-pc1-wd1.png")

filt, info  = fink_filter_bank(1, 8, pc=1, wd=2)
plot_filter_bank_QA(filt, info, fname="filt-8-pc1-wd2.png")

filt, info  = fink_filter_bank(1, 8, pc=1, wd=3)
plot_filter_bank_QA(filt, info, fname="filt-8-pc1-wd3.png")

filt, info  = fink_filter_bank(1, 8, pc=2, wd=1)
plot_filter_bank_QA(filt, info, fname="filt-8-pc2-wd1.png")

filt, info  = fink_filter_bank(1, 16, pc=2, wd=1)
plot_filter_bank_QA(filt, info, fname="filt-16-pc2-wd1.png")

filt, info  = fink_filter_bank(1, 16, pc=2, wd=2)
plot_filter_bank_QA(filt, info, fname="filt-16-pc2-wd2.png")

filt, info  = fink_filter_bank(1, 16, pc=1, wd=1)
plot_filter_bank_QA(filt, info, fname="filt-16-pc1-wd1.png")

filt, info  = fink_filter_bank(2, 8, pc=1, wd=1)
plot_filter_bank_QA(filt, info, fname="filt2-8-pc1-wd1.png")

filt, info  = fink_filter_bank(2, 8, pc=2, wd=1)
plot_filter_bank_QA(filt, info, fname="filt2-8-pc2-wd1.png")

filt, info  = fink_filter_bank(2, 16, pc=2, wd=1)
plot_filter_bank_QA(filt, info, fname="filt2-16-pc2-wd1.png")

filt, info  = fink_filter_bank(2, 16, pc=2, wd=2)
plot_filter_bank_QA(filt, info, fname="filt2-16-pc2-wd2.png")

filt, info  = fink_filter_bank(1, 8, pc=1, wd=1, shift=true)
plot_filter_bank_QA(filt, info, fname="filt-8-pc1-wd1-shift.png")

filt, info  = fink_filter_bank(1, 16, pc=2, wd=1, shift=true)
plot_filter_bank_QA(filt, info, fname="filt-16-pc2-wd1-shift.png")

filt, info  = fink_filter_bank(1, 16, pc=2, wd=2, shift=true)
plot_filter_bank_QA(filt, info, fname="filt-16-pc2-wd2-shift.png")
DH

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
# import stuff
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate(".")
@everywhere using Statistics
@everywhere using GRASS
@everywhere using SharedArrays
@everywhere using EchelleCCFs
using CSV
using DataFrames
using LaTeXStrings
using LsqFit

# define rms loop function
include(GRASS.moddir * "figures/fig_functions.jl")

# some global stuff
const N = 132
const Nt = 200
const Nloop = 200

# get command line args and output directories
run, plot = parse_args(ARGS)
grassdir, plotdir, datadir = check_plot_dirs()

function main()
    # set up parameters for lines
    lines = [5434.5]
    depths = range(0.05, stop=0.95, step=0.05)
    resolution=700000.0
    top = NaN
    contiguous_only=false

    # allocate shared arrays
    avg_avg_depth = SharedArray{Float64}(length(depths))
    std_avg_depth = SharedArray{Float64}(length(depths))
    avg_rms_depth = SharedArray{Float64}(length(depths))
    std_rms_depth = SharedArray{Float64}(length(depths))

    # loop over depths
    disk = DiskParams(N=N, Nt=Nt)
    @sync @distributed for i in eachindex(depths)
        println("running depth = " * string(depths[i]))

        # create spec instance
        spec = SpecParams(lines=lines, depths=[depths[i]], resolution=resolution,
                          extrapolate=true, contiguous_only=contiguous_only)

        # synthesize spectra, get velocities and stats
        avg_avg1, std_avg1, avg_rms1, std_rms1 = spec_loop(spec, disk, Nloop, top=top)
        avg_avg_depth[i] = avg_avg1
        std_avg_depth[i] = std_avg1
        avg_rms_depth[i] = avg_rms1
        std_rms_depth[i] = std_rms1
    end

    # make data frame
    df = DataFrame()
    df[!,:depths] = depths
    df[!,:avg_avg_depth] = avg_avg_depth
    df[!,:std_avg_depth] = std_avg_depth
    df[!,:avg_rms_depth] = avg_rms_depth
    df[!,:std_rms_depth] = std_rms_depth

    # write to CSV
    fname = datadir * "rms_vs_depth_" * string(N) * ".csv"
    CSV.write(fname, df)
    return nothing
end

# run the simulation
if run
    main()
end

if plot
    # plotting imports
    import PyPlot; plt = PyPlot; mpl = plt.matplotlib; plt.ioff()
    using PyCall; animation = pyimport("matplotlib.animation")
    mpl.style.use(GRASS.moddir * "figures/fig.mplstyle")

    # read in the data
    fname = datadir * "rms_vs_depth_" * string(N) * ".csv"
    df = CSV.read(fname, DataFrame)

    # assign to variable names
    depths = df.depths
    avg_avg_depth = df.avg_avg_depth
    std_avg_depth = df.std_avg_depth
    avg_rms_depth = df.avg_rms_depth
    std_rms_depth = df.std_rms_depth

    # get the errors
    err_avg_depth = std_avg_depth ./ sqrt(Nloop)
    err_rms_depth = std_rms_depth ./ sqrt(Nloop)

    # plot the results
    fig, ax1 = plt.subplots()
    ax1.errorbar(depths, avg_rms_depth, yerr=err_rms_depth, capsize=3.0, color="black", fmt=".")
    ax1.fill_between(depths, avg_rms_depth .- std_rms_depth, avg_rms_depth .+ std_rms_depth, color="tab:blue", alpha=0.3)
    ax1.fill_betweenx(range(0.0, 1.0, length=5), zeros(5), zeros(5) .+ 0.2, hatch="/", fc="black", ec="white", alpha=0.15, zorder=0)

    # set labels, etc.
    ax1.set_xlabel(L"{\rm Line\ Depth}")
    ax1.set_ylabel(L"{\rm RMS}_{\rm RV}\ {\rm (m s}^{-1})")
    ax1.set_xlim(0.0, 1.0)
    ax1.set_ylim(0.375, 0.675)

    # annotate the axes and save the figure
    arrowprops = Dict("facecolor"=>"black", "shrink"=>0.05, "width"=>2.0,"headwidth"=>8.0)
    ax1.annotate(L"{\rm Shallow}", xy=(0.85,0.39), xytext=(0.05,0.388), arrowprops=arrowprops)
    ax1.annotate(L"{\rm Deep}", xy=(0.86, 0.388))
    fig.savefig(plotdir * "fig5.pdf")
    plt.clf(); plt.close()
    println(">>> Figure written to: " * plotdir * "fig5.pdf")
end

ext: jl

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
#using Distances

""" Euclidian distance between data sets
get the euclidian distance between 2 multidimentional arrays
both arrays need to have the same dimentions (and each array within each array should also have the same dimentions)

Input:
	- d1 - first array
	- d2 - second array

Output:
	- Distance in Float64

"""
function euclidian_distance(d1::Array{Array{Float64,1},1},d2::Array{Array{Float64,1},1})

	#get the distance for each array from the two datasets)

	if ( size(d1)[1] != size(d2)[1] )
		error("euclidian_distance: data sets are of different dimentions")
		return
	end

	temp_matrix = Array{Array{Float64,1},1}(size(d1)[1])

	for i in 1:size(d1)[1]
		#get the distance between the vectors
		vec1 = d1[i]
		vec2 = d2[i]

		if ( size(vec1)[1] != size(vec2)[1] )
			error("euclidian_distance: one of the species has different dimentions in the two data sets")
			return
		end

		temp_matrix[i] = (vec1-vec2) .* (vec1-vec2)
	end

	return sqrt(sum(sum(temp_matrix)))

end

"""
Hellinger distance between 2 matrices (of same size)
"""
function hellinger_distance(m1::Array{Float64,2},m2::Array{Float64,2})

	d_sum = 0

	for i in 1:size(m1)[1]
		for j in 1:size(m1)[2]
			d_sum += (sqrt(m1[i,j])-sqrt(m2[i,j]))^2
		end
	end

	return (1/sqrt(2))*sqrt(d_sum)

end

"""
compute the transition matrix for an array of random walks (all need to have the same number of steps)
"""
function transition_matrix(walks::Array{Random_Walk,1};nrows::Int64=15,ncols::Int64=15)

	#create the matrix
	T = Array{Float64,2}(nrows,ncols)

	nbins=15
	half_nbins=round(Int,nbins/2)

	#fill the matrix
	for i in 1:nrows

		"""
		#calculate i interval from 0 to 2pi and convert
		#TODO: need to convert in the right order!!! from -pi to pi

		alpha_i_min , alpha_i_max = (2*pi*(i-1))/nbins , (2*pi*i)/nbins

		#convert to -pi to pi interval
		alpha_i_min = (alpha_i_min + pi)%(2*pi)-pi
		alpha_i_max = (alpha_i_max + pi)%(2*pi)-pi

		if i==15
			#TODO: find a way to correct in formula
			alpha_i_max = 1*pi
		end

		"""

		"""
		#calculate i interval from -pi to pi
		#TODO: find a way to make it 15 intervals
		if i<=half_nbins
			#the first half of the intervals goes from -pi to 0
			alpha_i_min , alpha_i_max = ((i-1-half_nbins)*pi)/half_nbins , ((i-half_nbins)*pi)/half_nbins

			if alpha_i_max==0.0
				#the middle interval goes from before 0.0 to after 0.0 (beacuse the number of intervals is 15)
				alpha_i_max = -alpha_i_min
			end
		else
			#the second half of the intervals goes from 0 to pi
			alpha_i_min , alpha_i_max = ((i-half_nbins)*pi)/half_nbins , ((i+1-half_nbins)*pi)/half_nbins
		end
		"""

		#16 intervals
		#alpha_i_min , alpha_i_max = ((i-1-half_nbins)*pi)/half_nbins , ((i-half_nbins)*pi)/half_nbins

		#calculate intervals from -pi to pi (15 intervals)
		if i-1==0
			alpha_i_min , alpha_i_max = -pi , -pi + (2pi/nbins)

		else
			alpha_i_min , alpha_i_max = -pi+(2pi/nbins*(i-1)) , -pi+(2pi/nbins*i)
		end

		#println("interval ", i ," ", alpha_i_min, ",", alpha_i_max)

		for j in 1:ncols

			"""
			#get the j interval
			alpha_j_min , alpha_j_max = (2*pi*(j-1))/nbins, (2*pi*j)/nbins
			#convert to -pi to pi interval
			alpha_j_min = (alpha_j_min + pi)%(2*pi)-pi
			alpha_j_max = (alpha_j_max + pi)%(2*pi)-pi

			if j==15
				#TODO: find a way to correct in formula
				alpha_j_max = 1*pi
			end
			"""

			"""
			if j<=half_nbins
				#the first half of the intervals goes from -pi to 0
				alpha_j_min , alpha_j_max = ((j-1-half_nbins)*pi)/half_nbins , ((j-half_nbins)*pi)/half_nbins

				if alpha_j_max==0.0
					#the middle interval goes from before 0.0 to after 0.0 (beacuse the number of intervals is 15)
					alpha_j_max = -alpha_j_min
				end
			else
				#the second half of the intervals goes from 0 to pi
				alpha_j_min , alpha_j_max = ((j-half_nbins)*pi)/half_nbins , ((j+1-half_nbins)*pi)/half_nbins
			end
			"""

			#16 intervals
			#alpha_j_min , alpha_j_max = ((j-1-half_nbins)*pi)/half_nbins , ((j-half_nbins)*pi)/half_nbins

			#calculate intervals from -pi to pi (15 intervals)
			if j-1==0
				alpha_j_min , alpha_j_max = -pi , -pi + (2pi/nbins)

			else
				alpha_j_min , alpha_j_max = -pi+(2pi/nbins*(j-1)) , -pi+(2pi/nbins*j)
			end

			#println("interval j: ", alpha_j_min, ",", alpha_j_max)

			current_sum = 0

			n_walks = length(walks)
			#n_steps = length(walks[1].steps_length)-1 #assuming each walk has the same number of steps
			all_steps=0

			for q in 1:n_walks
				#go over each walk
				angles = walks[q].steps_angles
				I_qr_ij = 0

				#in case number of stepts can change for each walk
				n_steps = length(angles)-1
				all_steps += n_steps

				for r in 1:n_steps #last time step is ignored because there is no next one for it
					#go over each step of the walk
					#get the angle for the current step and the next step
					alpha_q_r_i = angles[r+1]
					alpha_q_r_j = angles[r]

					if ( alpha_q_r_i>=alpha_i_min && alpha_q_r_i<alpha_i_max ) && ( alpha_q_r_j>=alpha_j_min &&  alpha_q_r_j<alpha_j_max )
						#both angles are within the interval
						I_qr_ij += 1
					end
				end #end r loop

				current_sum += I_qr_ij
			end #end q loop

			#current_sum = current_sum / (n_walks*n_steps) #assuming each walk has the same number of steps
			current_sum = current_sum / all_steps

			T[i,j] = current_sum

		end #end j loop
	end #end i loop

	return T

end


text:
module TestNodeCommand

using Comonicon.AST
using Comonicon.JuliaExpr
using Comonicon.JuliaExpr: emit, emit_body, emit_norm_body, emit_dash_body
using Test

function called()
    @test true
end

cmd = Entry(;
    version = v"1.2.0",
    root = NodeCommand(;
        name = "node",
        subcmds = Dict(
            "cmd1" => LeafCommand(; fn = called, name = "cmd1"),
            "cmd2" => LeafCommand(; fn = called, name = "cmd2"),
        ),
    ),
)

eval(emit(cmd))

@testset "test node" begin
    @test command_main(["cmd3"]) == 1
    @test command_main(["cmd1", "foo"]) == 1
    @test command_main(["cmd1", "foo", "-h"]) == 0
    @test command_main(["cmd1", "foo", "-V"]) == 0
    @test command_main(String[]) == 1
end

end


text:
using PBDS
using Test
using NBInclude

@testset "Examples" begin
    "no_plots" in ARGS && (global const no_plots = true)
    PBDS_dir = joinpath(@__DIR__, "..", "examples", "PBDS")

    test_notebook = "R7Arm_DynamicMugGrasping"
    @testset "$test_notebook" begin
        file = string(test_notebook, ".ipynb")
        @nbinclude(joinpath(PBDS_dir, file))
        @test root.children[end-1].children[3].traj_log.x[end][1] .< 5e-3
        robot_coord_rep = ChartRep()
        traj = propagate_tasks(xm, vm, root, CM, Time, dt, robot_coord_rep, state, cache,
            mugparams; time_dep, log_tasks=true)
        @test root.children[end-1].children[3].traj_log.x[end][1] .< 5e-3
        println("Finished example ", test_notebook)
    end

    test_notebooks = ["R1_To_R1PointPositionAttractor",
                      "R2_To_R1PointDistanceAttractor_S1Damping",
                      "R2_To_R1PointDistanceAttractor_S1Damping_R1BoxAvoidance",
                      "R2_To_R1PointDistanceAttractor_S1Damping_R1SphereAvoidance",
                      "R2_To_R2PointPositionAttractor",
                      "R3_To_R1PointDistanceAttractor_S2Damping_R1BoxAvoidance",
                      "R3_To_R1PointDistanceAttractor_S2Damping_R1SphereAvoidance",
                      "S2_To_R1Attractor_S2Damping_R1ObstacleAvoidance"]

    for test_notebook in test_notebooks
        @testset "$test_notebook" begin
            file = string(test_notebook, ".ipynb")
            @nbinclude(joinpath(PBDS_dir, file))
            ε = 5e-2
            @test all([norm(traj.xm[end] - xm_goal) for traj in trajs] .< ε)

            if test_notebook == "R3_To_R1PointDistanceAttractor_S2Damping_R1SphereAvoidance"
                log_tasks = true
                robot_coord_rep = ChartRep()
                traj = propagate_tasks(xm, vm, tasks, CM, CNs, Time, dt, robot_coord_rep; log_tasks)
                @test norm(traj.xm[end] - xm_goal) < ε
                robot_coord_rep = EmbRep()
                traj = propagate_tasks(xm, vm, tasks, CM, CNs, Time, dt, robot_coord_rep; log_tasks)
                @test norm(traj.xm[end] - xm_goal) < ε
            end
            println("Finished example ", test_notebook)
        end
    end

    test_notebook = "S2_To_R1Attractor_S2Damping_ConsistencyTest"
    @testset "$test_notebook" begin
        file = string(test_notebook, ".ipynb")
        @nbinclude(joinpath(PBDS_dir, file))
        Δx = 0.
        Δx += sum(@. norm(traj_north.xm - traj_south.xm))
        Δx += sum(@. norm(traj_south.xm - traj_switching.xm))
        Δx += sum(@. norm(traj_switching.xm - traj_north.xm))
        @test Δx < 1e-4
        ε = 5e-3
        @test norm(traj_north.xm[end] - xm_goal) < ε

        RMP_dir = joinpath(@__DIR__, "..", "examples", "RMPflow")
        @nbinclude(joinpath(RMP_dir, file))
        Δx = 0.
        Δx += sum(@. norm(traj_north.xm - traj_south.xm))
        Δx += sum(@. norm(traj_south.xm - traj_switching.xm))
        Δx += sum(@. norm(traj_switching.xm - traj_north.xm))
        @test Δx > 1e3
        Δx = 0.
        Δx += norm(traj_north.xm[end] - xm_goal)
        Δx += norm(traj_south.xm[end] - xm_goal)
        Δx += norm(traj_switching.xm[end] - xm_goal)
        @test Δx < ε
        println("Finished example ", test_notebook)
    end

end


text:
@testset "172.factorial-trailing-zeroes.jl" begin
    @test all(trailing_zeroes(i) == 0 for i in 0:4)
    @test trailing_zeroes(5) == 1
    @test trailing_zeroes(1999) == 496
    @test trailing_zeroes(2000) == trailing_zeroes(2001) == 499
end


text:
module TestRefCount

using AMGX, Test
using AMGX: Config, Resources, AMGXVector, dDDI, AMGXMatrix, Solver, RefCountError

@testset "ref counts" begin
    c = Config("")
    @test c.ref_count[] == 0
    r = Resources(c)
    @test c.ref_count[] == 1
    @test_throws RefCountError close(c)
    v = AMGXVector(r, dDDI)
    @test_throws RefCountError close(r)
    @test r.ref_count[] == 1
    close(v)
    @test r.ref_count[] == 0
    m = AMGXMatrix(r, dDDI)
    @test r.ref_count[] == 1
    @test_throws RefCountError close(r)
    close(m)
    m = AMGXMatrix(r, dDDI)
    s = Solver(r, dDDI, c)
    @test c.ref_count[] == 2
    @test r.ref_count[] == 2
    close(s)
    close(m)
    close(r)
    close(c)
end

end #module


text:
using Test
using ZipFile

Debug = false

@test Any[] == detect_ambiguities(Base, Core, ZipFile)

function findfile(dir, name)
    for f in dir.files
        if f.name == name
            return f
        end
    end
    nothing
end

function fileequals(f, s)
    read(f, String) == s
end

# test a zip file created using Info-Zip
dir = ZipFile.Reader(joinpath(dirname(@__FILE__), "infozip.zip"))
@test length(dir.files) == 4

f = findfile(dir, "ziptest/")
@test f.method == ZipFile.Store
@test f.uncompressedsize == 0
@test fileequals(f, "")

f = findfile(dir, "ziptest/hello.txt")
@test fileequals(f, "hello world!\n")

f = findfile(dir, "ziptest/info.txt")
@test fileequals(f, "Julia\nfor\ntechnical computing\n")

f = findfile(dir, "ziptest/julia.txt")
@test f.method == ZipFile.Deflate
@test fileequals(f, repeat("Julia\n", 10))

close(dir)

tmp = mktempdir()
if Debug
    println("temporary directory $tmp")
end

# write an empty zip file
dir = ZipFile.Writer("$tmp/empty.zip")
close(dir)
dir = ZipFile.Reader("$tmp/empty.zip")
@test length(dir.files) == 0

# write and then read back a zip file
zipdata = [
    ("hello.txt", "hello world!\n", ZipFile.Store),
    ("info.txt", "Julia\nfor\ntechnical computing\n", ZipFile.Store),
    ("julia.txt", "julia\n"^10, ZipFile.Deflate),
    ("empty1.txt", "", ZipFile.Store),
    ("empty2.txt", "", ZipFile.Deflate),
]
# 2013-08-16	9:42:24
modtime = time(Libc.TmStruct(24, 42, 9, 16, 7, 2013-1900, 0, 0, -1))

dir = ZipFile.Writer("$tmp/hello.zip")
@test length(string(dir)) > 0
for (name, data, meth) in zipdata
    local f = ZipFile.addfile(dir, name; method=meth, mtime=modtime)
    @test length(string(f)) > 0
    write(f, data)
end
close(dir)

dir = ZipFile.Reader("$tmp/hello.zip")
@test length(string(dir)) > 0
for (name, data, meth) in zipdata
    local f = findfile(dir, name)
    @test length(string(f)) > 0
    @test f.method == meth
    @test abs(mtime(f) - modtime) < 2
    @test fileequals(f, data)
end
close(dir)

s1 = "this is an example sentence"
s2 = ". hello world.\n"
filename = "$tmp/multi.zip"
dir = ZipFile.Writer(filename)
f = ZipFile.addfile(dir, "data"; method=ZipFile.Deflate)
write(f, s1)
write(f, s2)
close(dir)
dir = ZipFile.Reader(filename)
@test String(read!(dir.files[1], Array{UInt8}(undef, length(s1)))) == s1
@test String(read!(dir.files[1], Array{UInt8}(undef, length(s2)))) == s2
@test eof(dir.files[1])
close(dir)

data = Any[
    UInt8(20),
    Int(42),
    float(3.14),
    "julia",
    rand(5),
    rand(3, 4),
    view(rand(10,10), 2:8,2:4),
]
filename = "$tmp/multi2.zip"
dir = ZipFile.Writer(filename)
f = ZipFile.addfile(dir, "data"; method=ZipFile.Deflate)
@test_throws ErrorException read!(f, Array{UInt8}(undef, 1))
for x in data
    write(f, x)
end
close(dir)

dir = ZipFile.Reader(filename)
@test_throws ErrorException write(dir.files[1], UInt8(20))
for x in data
    if isa(x, String)
        @test x == String(read!(dir.files[1], Array{UInt8}(undef, length(x))))
    elseif isa(x, Array)
        y = similar(x)
        y[:] .= 0
        @test x == read!(dir.files[1], y)
        @test x == y
    elseif isa(x, SubArray)
        continue # Base knows how to write, but not read
    else
        @test x == read(dir.files[1], typeof(x))
    end
end
close(dir)

if !Debug
    rm(tmp, recursive=true)
end

println("done")


text:
#=
--------------------------------------------------------------------------------------------------------------------------------
Convergence
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    vec_ecm_params(estim::EstimSettings, B_star::SubArray{Float64}, C_star::Nothing, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::Nothing)

Return the vec form of the free parameters in B_star and Q_view.

    vec_ecm_params(estim::EstimSettings, B_star::Nothing, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::Nothing, coordinates_free_params_C::CoordinatesVector)

Return the vec form of the free parameters in C_star and Q_view.

    vec_ecm_params(estim::EstimSettings, B_star::SubArray{Float64}, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::CoordinatesVector)

Return the vec form of the free parameters in B_star, C_star and Q_view.
"""
function vec_ecm_params(estim::EstimSettings, B_star::SubArray{Float64}, C_star::Nothing, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::Nothing)

    # Setup and memory pre-allocation
    counter = 1;
    output = zeros(length(coordinates_free_params_B) + length(Q_view));

    # Add parameters in B_star to output
    for ij in coordinates_free_params_B
        output[counter] = B_star[ij];
        counter += 1;
    end

    # Add parameters in Q_view to output
    for ij in eachindex(Q_view)
        output[counter] = Q_view[ij];
        counter += 1;
    end

    # Return output
    return output;
end

function vec_ecm_params(estim::EstimSettings, B_star::Nothing, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::Nothing, coordinates_free_params_C::CoordinatesVector)

    # Setup and memory pre-allocation
    counter = 1;
    output = zeros(length(coordinates_free_params_C) + length(Q_view));

    # Add parameters in C_star to output
    for ij in coordinates_free_params_C
        output[counter] = C_star[ij];
        counter += 1;
    end

    # Add parameters in Q_view to output
    for ij in eachindex(Q_view)
        output[counter] = Q_view[ij];
        counter += 1;
    end

    # Return output
    return output;
end

function vec_ecm_params(estim::EstimSettings, B_star::SubArray{Float64}, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::CoordinatesVector)

    # Setup and memory pre-allocation
    counter = 1;
    output = zeros(length(coordinates_free_params_B) + length(coordinates_free_params_C) + length(Q_view));

    # Add parameters in B_star to output
    for ij in coordinates_free_params_B
        output[counter] = B_star[ij];
        counter += 1;
    end

    # Add parameters in C_star to output
    for ij in coordinates_free_params_C
        output[counter] = C_star[ij];
        counter += 1;
    end

    # Add parameters in Q_view to output
    for ij in eachindex(Q_view)
        output[counter] = Q_view[ij];
        counter += 1;
    end

    # Return output
    return output;
end

"""
    vec_ecm_params!(estim::EstimSettings, output::FloatVector, B_star::SubArray{Float64}, C_star::Nothing, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::Nothing)

Update output in-place with the vec form of the free parameters in B_star and Q_view.

    vec_ecm_params!(estim::EstimSettings, output::FloatVector, B_star::Nothing, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::Nothing, coordinates_free_params_C::CoordinatesVector)

Update output in-place with the vec form of the free parameters in C_star and Q_view.

    vec_ecm_params!(estim::EstimSettings, output::FloatVector, B_star::SubArray{Float64}, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::CoordinatesVector)

Update output in-place with the vec form of the free parameters in B_star, C_star and Q_view.
"""
function vec_ecm_params!(estim::EstimSettings, output::FloatVector, B_star::SubArray{Float64}, C_star::Nothing, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::Nothing)

    # Setup and memory pre-allocation
    counter = 1;

    # Add parameters in B_star to output
    for ij in coordinates_free_params_B
        output[counter] = B_star[ij];
        counter += 1;
    end

    # Add parameters in Q_view to output
    for ij in eachindex(Q_view)
        output[counter] = Q_view[ij];
        counter += 1;
    end
end

function vec_ecm_params!(estim::EstimSettings, output::FloatVector, B_star::Nothing, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::Nothing, coordinates_free_params_C::CoordinatesVector)

    # Setup and memory pre-allocation
    counter = 1;

    # Add parameters in C_star to output
    for ij in coordinates_free_params_C
        output[counter] = C_star[ij];
        counter += 1;
    end

    # Add parameters in Q_view to output
    for ij in eachindex(Q_view)
        output[counter] = Q_view[ij];
        counter += 1;
    end
end

function vec_ecm_params!(estim::EstimSettings, output::FloatVector, B_star::SubArray{Float64}, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::CoordinatesVector)

    # Setup and memory pre-allocation
    counter = 1;

    # Add parameters in B_star to output
    for ij in coordinates_free_params_B
        output[counter] = B_star[ij];
        counter += 1;
    end

    # Add parameters in C_star to output
    for ij in coordinates_free_params_C
        output[counter] = C_star[ij];
        counter += 1;
    end

    # Add parameters in Q_view to output
    for ij in eachindex(Q_view)
        output[counter] = Q_view[ij];
        counter += 1;
    end
end

"""
    update_relative_change!(estim::EstimSettings, params_rel_change::FloatVector, ecm_params_new::FloatVector, ecm_params_old::FloatVector)

Update `params_rel_change` in-place.
"""
@inline function update_relative_change!(estim::EstimSettings, params_rel_change::FloatVector, ecm_params_new::FloatVector, ecm_params_old::FloatVector)
    @turbo for i in eachindex(params_rel_change)
        params_rel_change[i] = abs(ecm_params_new[i] - ecm_params_old[i])/(abs(ecm_params_old[i]) + estim.ε);
    end
end

"""
    ecm_converged_median(estim::EstimSettings, params_rel_change::FloatVector, iter::Int64)

Check whether the ecm has converged looking at the median of the parameters relative change.
"""
function ecm_converged_median(estim::EstimSettings, params_rel_change::FloatVector, iter::Int64)

    # Median relative change
    median_params_rel_change = median(params_rel_change);
    if iter-estim.prerun-1 < 10
        verb_message(estim.verb, "ecm > iter=0$(iter-estim.prerun-1), relative parameter convergence criterion=$(round(median_params_rel_change, digits=16))");
    else
        verb_message(estim.verb, "ecm > iter=$(iter-estim.prerun-1), relative parameter convergence criterion=$(round(median_params_rel_change, digits=16))");
    end

    # Check convergence
    convergence_criterion = median_params_rel_change <= estim.tol;
    return convergence_criterion;
end

"""
    ecm_converged_complete(estim::EstimSettings, params_rel_change::FloatVector, iter::Int64)

Check whether the ecm has converged looking at both the median and 95th quantile of the parameters relative change.
"""
function ecm_converged_complete(estim::EstimSettings, params_rel_change::FloatVector, iter::Int64)

    # Median relative change
    median_params_rel_change = median(params_rel_change);
    quantile_params_rel_change = quantile(params_rel_change, 0.95);
    if iter-estim.prerun-1 < 10
        verb_message(estim.verb, "ecm > iter=0$(iter-estim.prerun-1), median convergence criterion=$(round(median_params_rel_change, digits=16)), 95th quantile convergence criterion=$(round(quantile_params_rel_change, digits=16))");
    else
        verb_message(estim.verb, "ecm > iter=$(iter-estim.prerun-1), median convergence criterion=$(round(median_params_rel_change, digits=16)), 95th quantile convergence criterion=$(round(quantile_params_rel_change, digits=16))");
    end

    # Check convergence
    convergence_criterion = (median_params_rel_change <= estim.tol) & (quantile_params_rel_change <= 10*estim.tol)
    return convergence_criterion;
end

"""
    ecm_converged(estim::EstimSettings, params_rel_change::FloatVector, iter::Int64)

Check whether the ecm has converged.
"""
function ecm_converged(estim::EstimSettings, params_rel_change::FloatVector, iter::Int64)
    if estim.check_quantile == false
        return ecm_converged_median(estim, params_rel_change, iter)::Bool;
    else
        return ecm_converged_complete(estim, params_rel_change, iter)::Bool;
    end
end

#=
--------------------------------------------------------------------------------------------------------------------------------
ECM stats
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    initialise_ecm_stats_measurement(estim::EstimSettings, coordinates_measurement_states::IntVector)

Initialise M, N and O.

    initialise_ecm_stats_measurement(estim::EstimSettings, coordinates_measurement_states::Nothing)

Initialise M, N and O to nothing.
"""
function initialise_ecm_stats_measurement(estim::EstimSettings, coordinates_measurement_states::IntVector)
    M = zeros(estim.n, estim.m);
    N = Vector{SparseMatrixCSC{Float64, Int64}}(undef, estim.T);
    O = Array{VectorsArray{Float64},1}(undef, estim.T);
    buffer_M = zeros(estim.n);
    buffer_O = Symmetric(zeros(estim.m, estim.m));
    return M, N, O, buffer_M, buffer_O;
end

initialise_ecm_stats_measurement(estim::EstimSettings, coordinates_measurement_states::Nothing) = nothing, nothing, nothing, nothing, nothing, nothing;

"""
    initialise_ecm_stats_transition(coordinates_transition_current::IntVector, coordinates_transition_lagged::IntVector, coordinates_transition_PPs::IntVector)

Initialise F, G and H.

    initialise_ecm_stats_transition(coordinates_transition_current::IntVector, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing)

Initialise F, G and H (G and H to nothing)

    initialise_ecm_stats_transition(coordinates_transition_current::Nothing, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing)

Initialise F, G and H to nothing.
"""
function initialise_ecm_stats_transition(coordinates_transition_current::IntVector, coordinates_transition_lagged::IntVector, coordinates_transition_PPs::IntVector)
    F = zeros(length(coordinates_transition_current), length(coordinates_transition_current));
    G = zeros(length(coordinates_transition_current), length(coordinates_transition_PPs));
    H = zeros(length(coordinates_transition_lagged), length(coordinates_transition_lagged));

    return F, G, H;
end

function initialise_ecm_stats_transition(coordinates_transition_current::IntVector, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing)
    F = zeros(length(coordinates_transition_current), length(coordinates_transition_current));
    return F, nothing, nothing;
end

initialise_ecm_stats_transition(coordinates_transition_current::Nothing, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing) = nothing, nothing, nothing;

"""
    update_ecm_stats_measurement!(barrier_M::FloatMatrix, estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::IntVector, ind_not_missings::IntVector, t::Int64, Xs::FloatVector, Ps::SymMatrix)

Update the ECM statistics for the measurement equation.

    update_ecm_stats_measurement!(barrier_M::FloatMatrix, estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::IntVector, ind_not_missings::Nothing, t::Int64, Xs::FloatVector, Ps::SymMatrix)
    update_ecm_stats_measurement!(barrier_M::Nothing, estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::Nothing, ind_not_missings::IntVector, t::Int64, Xs::FloatVector, Ps::SymMatrix)
    update_ecm_stats_measurement!(barrier_M::Nothing, estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::Nothing, ind_not_missings::Nothing, t::Int64, Xs::FloatVector, Ps::SymMatrix)

Do not update ECM statistics for the measurement equation of models for which these calculations are not required.
"""
function update_ecm_stats_measurement!(barrier_M::FloatMatrix, estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::IntVector, ind_not_missings::IntVector, t::Int64, Xs::FloatVector, Ps::SymMatrix)

    # Views
    Y_obs = @view estim.Y[ind_not_missings, t];
    Xs_view = @view Xs[coordinates_measurement_states];
    Ps_view = @view Ps[coordinates_measurement_states, coordinates_measurement_states];

    # Necessary statistics
    A_transpose = spzeros(estim.n, length(ind_not_missings));
    for (i, j) in enumerate(ind_not_missings)
        A_transpose[j, i] = 1.0;
    end

    # Update ECM statistics: M
    mul!(smoother_arrays.buffer_M, A_transpose, Y_obs);
    mul!(smoother_arrays.M, smoother_arrays.buffer_M, Xs_view', 1.0, 1.0);

    # Update ECM statistics: compute and store N_t
    smoother_arrays.N[t] = A_transpose*A_transpose';

    # Update ECM statistics: compute O_t
    mul!(smoother_arrays.buffer_O.data, Xs_view, Xs_view');
    smoother_arrays.buffer_O.data .+= Ps_view;

    # Store O_t
    smoother_arrays.O[t] = [col[:] for col in eachcol(smoother_arrays.buffer_O)];
end

update_ecm_stats_measurement!(barrier_M::FloatMatrix, estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::IntVector, ind_not_missings::Nothing, t::Int64, Xs::FloatVector, Ps::SymMatrix) = nothing;
update_ecm_stats_measurement!(barrier_M::Nothing, estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::Nothing, ind_not_missings::IntVector, t::Int64, Xs::FloatVector, Ps::SymMatrix) = nothing;
update_ecm_stats_measurement!(barrier_M::Nothing, estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::Nothing, ind_not_missings::Nothing, t::Int64, Xs::FloatVector, Ps::SymMatrix) = nothing;

"""
    call_update_ecm_stats_measurement!(estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::IntVector, ind_not_missings::IntVector, t::Int64, Xs::FloatVector, Ps::SymMatrix)
    call_update_ecm_stats_measurement!(estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::IntVector, ind_not_missings::Nothing, t::Int64, Xs::FloatVector, Ps::SymMatrix)
    call_update_ecm_stats_measurement!(estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::Nothing, ind_not_missings::IntVector, t::Int64, Xs::FloatVector, Ps::SymMatrix)
    call_update_ecm_stats_measurement!(estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::Nothing, ind_not_missings::Nothing, t::Int64, Xs::FloatVector, Ps::SymMatrix)

APIs to call `update_ecm_stats_measurement!` with SmootherArrays.
"""
call_update_ecm_stats_measurement!(estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::IntVector, ind_not_missings::IntVector, t::Int64, Xs::FloatVector, Ps::SymMatrix) = update_ecm_stats_measurement!(smoother_arrays.M, estim, smoother_arrays, coordinates_measurement_states, ind_not_missings, t, Xs, Ps);
call_update_ecm_stats_measurement!(estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::IntVector, ind_not_missings::Nothing, t::Int64, Xs::FloatVector, Ps::SymMatrix) = update_ecm_stats_measurement!(smoother_arrays.M, estim, smoother_arrays, coordinates_measurement_states, ind_not_missings, t, Xs, Ps);
call_update_ecm_stats_measurement!(estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::Nothing, ind_not_missings::IntVector, t::Int64, Xs::FloatVector, Ps::SymMatrix) = update_ecm_stats_measurement!(smoother_arrays.M, estim, smoother_arrays, coordinates_measurement_states, ind_not_missings, t, Xs, Ps);
call_update_ecm_stats_measurement!(estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::Nothing, ind_not_missings::Nothing, t::Int64, Xs::FloatVector, Ps::SymMatrix) = update_ecm_stats_measurement!(smoother_arrays.M, estim, smoother_arrays, coordinates_measurement_states, ind_not_missings, t, Xs, Ps);

"""
    update_ecm_stats_transition!(barrier_F::FloatMatrix, barrier_G::FloatMatrix, estim::EstimSettings, smoother_arrays::SmootherArrays, Xs_lagged::FloatVector, Ps_lagged::SymMatrix, coordinates_transition_current::IntVector, coordinates_transition_lagged::IntVector, coordinates_transition_PPs::IntVector)

Update the ECM statistics for the transition equation.

    update_ecm_stats_transition!(barrier_F::FloatMatrix, barrier_G::Nothing, estim::EstimSettings, smoother_arrays::SmootherArrays, Xs_lagged::FloatVector, Ps_lagged::SymMatrix, coordinates_transition_current::IntVector, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing)

Update F only.

    update_ecm_stats_transition!(barrier_F::Nothing, barrier_G::Nothing, estim::EstimSettings, smoother_arrays::SmootherArrays, Xs_lagged::FloatVector, Ps_lagged::SymMatrix, coordinates_transition_current::Nothing, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing)

Do not update ECM statistics for the transition equation of models for which these calculations are not required.
"""
function update_ecm_stats_transition!(barrier_F::FloatMatrix, barrier_G::FloatMatrix, estim::EstimSettings, smoother_arrays::SmootherArrays, Xs_lagged::FloatVector, Ps_lagged::SymMatrix, coordinates_transition_current::IntVector, coordinates_transition_lagged::IntVector, coordinates_transition_PPs::IntVector)

    # Views
    Xs_view = @view smoother_arrays.Xs_leading[coordinates_transition_current];
    Ps_view = @view smoother_arrays.Ps_leading[coordinates_transition_current, coordinates_transition_current];
    Xs_lagged_view = @view Xs_lagged[coordinates_transition_lagged];
    Ps_lagged_view = @view Ps_lagged[coordinates_transition_lagged, coordinates_transition_lagged];
    PPs_view = @view smoother_arrays.Ps_leading[coordinates_transition_current, coordinates_transition_PPs];

    # Update ECM statistics
    smoother_arrays.F .+= Xs_view*Xs_view' + Ps_view;
    smoother_arrays.G .+= Xs_view*Xs_lagged_view' + PPs_view;
    smoother_arrays.H .+= Xs_lagged_view*Xs_lagged_view' + Ps_lagged_view;
end

function update_ecm_stats_transition!(barrier_F::FloatMatrix, barrier_G::Nothing, estim::EstimSettings, smoother_arrays::SmootherArrays, Xs_lagged::FloatVector, Ps_lagged::SymMatrix, coordinates_transition_current::IntVector, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing)

    # Views
    Xs_view = @view smoother_arrays.Xs_leading[coordinates_transition_current];
    Ps_view = @view smoother_arrays.Ps_leading[coordinates_transition_current, coordinates_transition_current];

    # Update ECM statistics
    smoother_arrays.F .+= Xs_view*Xs_view' + Ps_view;
end

update_ecm_stats_transition!(barrier_F::Nothing, barrier_G::Nothing, estim::EstimSettings, smoother_arrays::SmootherArrays, Xs_lagged::FloatVector, Ps_lagged::SymMatrix, coordinates_transition_current::Nothing, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing) = nothing;

"""
    call_update_ecm_stats_transition!(estim::EstimSettings, smoother_arrays::SmootherArrays, Xs::FloatVector, Ps::SymMatrix, coordinates_transition_current::IntVector, coordinates_transition_lagged::IntVector, coordinates_transition_PPs::IntVector)
    call_update_ecm_stats_transition!(estim::EstimSettings, smoother_arrays::SmootherArrays, Xs::FloatVector, Ps::SymMatrix, coordinates_transition_current::IntVector, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing)
    call_update_ecm_stats_transition!(estim::EstimSettings, smoother_arrays::SmootherArrays, Xs::FloatVector, Ps::SymMatrix, coordinates_transition_current::Nothing, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing)

APIs to call `update_ecm_stats_transition!` with SmootherArrays.
"""
call_update_ecm_stats_transition!(estim::EstimSettings, smoother_arrays::SmootherArrays, Xs::FloatVector, Ps::SymMatrix, coordinates_transition_current::IntVector, coordinates_transition_lagged::IntVector, coordinates_transition_PPs::IntVector) = update_ecm_stats_transition!(smoother_arrays.F, smoother_arrays.G, estim, smoother_arrays, Xs, Ps, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);
call_update_ecm_stats_transition!(estim::EstimSettings, smoother_arrays::SmootherArrays, Xs::FloatVector, Ps::SymMatrix, coordinates_transition_current::IntVector, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing) = update_ecm_stats_transition!(smoother_arrays.F, smoother_arrays.G, estim, smoother_arrays, Xs, Ps, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);
call_update_ecm_stats_transition!(estim::EstimSettings, smoother_arrays::SmootherArrays, Xs::FloatVector, Ps::SymMatrix, coordinates_transition_current::Nothing, coordinates_transition_lagged::Nothing, coordinates_transition_PPs::Nothing) = update_ecm_stats_transition!(smoother_arrays.F, smoother_arrays.G, estim, smoother_arrays, Xs, Ps, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);

"""
    reinitialise_ecm_stats!(A::Nothing)
    reinitialise_ecm_stats!(A::FloatVector)
    reinitialise_ecm_stats!(A::FloatMatrix)
    reinitialise_ecm_stats!(A::Array{VectorsArray{Float64},1})
    reinitialise_ecm_stats!(smoother_arrays::SmootherArrays)

Re-initialise ECM statistics for the following iteration.
"""
reinitialise_ecm_stats!(A::Nothing) = nothing;
reinitialise_ecm_stats!(A::FloatVector) = fill!(A, 0.0);
reinitialise_ecm_stats!(A::FloatMatrix) = fill!(A, 0.0);
reinitialise_ecm_stats!(A::SymMatrix) = fill!(A.data, 0.0);
reinitialise_ecm_stats!(A::Array{VectorsArray{Float64},1}) = nothing; # the items of A are replaced in-place in the ecm iteration, when needed - thus, empty!(A) is not needed.
reinitialise_ecm_stats!(A::Vector{SparseMatrixCSC{Float64, Int64}}) = nothing; # the items of A are replaced in-place in the ecm iteration, when needed - thus, empty!(A) is not needed.

function reinitialise_ecm_stats!(smoother_arrays::SmootherArrays)
    reinitialise_ecm_stats!(smoother_arrays.J1);
    reinitialise_ecm_stats!(smoother_arrays.J2);
    reinitialise_ecm_stats!(smoother_arrays.Xs_leading);
    reinitialise_ecm_stats!(smoother_arrays.Ps_leading);
    reinitialise_ecm_stats!(smoother_arrays.F);
    reinitialise_ecm_stats!(smoother_arrays.G);
    reinitialise_ecm_stats!(smoother_arrays.H);
    reinitialise_ecm_stats!(smoother_arrays.M);
    reinitialise_ecm_stats!(smoother_arrays.N);
    reinitialise_ecm_stats!(smoother_arrays.O);
end

#=
--------------------------------------------------------------------------------------------------------------------------------
Kalman smoother and CM-step for initial conditions (i.e., X0 and P0)
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    call_update_smoothing_factors!(sspace::KalmanSettings, status::SizedKalmanStatus, smoother_arrays::SmootherArrays, ind_not_missings::Union{IntVector, Nothing}, e::Union{FloatVector, Nothing}, inv_F::Union{SymMatrix, Nothing}, L::Union{FloatMatrix, Nothing})
    call_update_smoothing_factors!(sspace::KalmanSettings, status::SizedKalmanStatus, smoother_arrays::SmootherArrays, ind_not_missings::Union{IntVector, Nothing})

APIs to call `MessyTimeSeries.update_smoothing_factors!` with SmootherArrays.
"""
call_update_smoothing_factors!(sspace::KalmanSettings, status::SizedKalmanStatus, smoother_arrays::SmootherArrays, ind_not_missings::Union{IntVector, Nothing}, e::Union{FloatVector, Nothing}, inv_F::Union{SymMatrix, Nothing}, L::Union{FloatMatrix, Nothing}) = update_smoothing_factors!(sspace, status, ind_not_missings, smoother_arrays.J1, smoother_arrays.J2, e, inv_F, L);
call_update_smoothing_factors!(sspace::KalmanSettings, status::SizedKalmanStatus, smoother_arrays::SmootherArrays, ind_not_missings::Union{IntVector, Nothing}) = update_smoothing_factors!(sspace, status, ind_not_missings, smoother_arrays.J1, smoother_arrays.J2);

"""
    call_backwards_pass(Xp::FloatVector, Pp::SymMatrix, smoother_arrays::SmootherArrays)
    call_backwards_pass(Pp::SymMatrix, smoother_arrays::SmootherArrays)

APIs to call `MessyTimeSeries.backwards_pass` with SmootherArrays.
"""
call_backwards_pass(Xp::FloatVector, Pp::SymMatrix, smoother_arrays::SmootherArrays) = backwards_pass(Xp, Pp, smoother_arrays.J1);
call_backwards_pass(Pp::SymMatrix, smoother_arrays::SmootherArrays) = backwards_pass(Pp, smoother_arrays.J2);

"""
    ksmoother_ecm_iteration!(sspace::KalmanSettings, status::SizedKalmanStatus, smoother_arrays::SmootherArrays, t::Int64)

Kalman smoother iteration of `ksmoother_ecm!`
"""
function ksmoother_ecm_iteration!(sspace::KalmanSettings, status::SizedKalmanStatus, smoother_arrays::SmootherArrays, t::Int64)

    # Pointers
    Xp = status.history_X_prior[t];
    Pp = status.history_P_prior[t];
    e = status.history_e[t];
    inv_F = status.history_inv_F[t];
    L = status.history_L[t];

    # Handle missing observations
    ind_not_missings = find_observed_data(sspace, t);

    # Smoothed estimates for t
    call_update_smoothing_factors!(sspace, status, smoother_arrays, ind_not_missings, e, inv_F, L);
    Xs = call_backwards_pass(Xp, Pp, smoother_arrays)::FloatVector;
    Ps = call_backwards_pass(Pp, smoother_arrays)::SymMatrix;

    # Return smoothed estimates for t
    return Xs, Ps, ind_not_missings;
end

"""
    cm_step_X0_P0!(sspace::KalmanSettings, status::SizedKalmanStatus, smoother_arrays::SmootherArrays, coordinates_transition_P0::IntVector)

CM-step for X0 and P0.
"""
function cm_step_X0_P0!(sspace::KalmanSettings, status::SizedKalmanStatus, smoother_arrays::SmootherArrays, coordinates_transition_P0::IntVector)

    # CM-step for X0
    mul!(sspace.X0, sspace.P0, smoother_arrays.J1, 1.0, 1.0);

    # CM-step for sspace.P0
    mul!(status.online_status.buffer_J2, sspace.P0, Array(smoother_arrays.J2));         # Array(...) helps speeding up mul!(...) while keeping J2 symmetric (i.e., good trade-off wrt other options incl. ".data")
    mul!(status.online_status.buffer_m_m, status.online_status.buffer_J2, sspace.P0);
    for ij in coordinates_transition_P0
        sspace.P0.data[ij] -= status.online_status.buffer_m_m[ij];
    end
end

"""
    ksmoother_ecm!(estim::EstimSettings, sspace::KalmanSettings, status::SizedKalmanStatus, coordinates_measurement_states::Union{IntVector, Nothing}, coordinates_transition_current::Union{IntVector, Nothing}, coordinates_transition_lagged::Union{IntVector, Nothing}, coordinates_transition_PPs::Union{IntVector, Nothing}, coordinates_transition_P0::IntVector)

Kalman smoother: RTS smoother from the last evaluated time period in status to t==0.

The smoother is implemented following the approach proposed in Durbin and Koopman (2012).

This instance of the smoother returns the ECM statistics and updates the initial conditions in KalmanSettings.

# Arguments
- `estim`: EstimSettings struct
- `sspace`: KalmanSettings struct
- `status`: SizedKalmanStatus struct
- `smoother_arrays`: SmootherArrays struct
- `coordinates_measurement_states`: indices describing the states needed to construct M, N, O (or nothing)
- `coordinates_transition_current`, `coordinates_transition_lagged` and `coordinates_transition_PPs`: indices describing the elements of interest needed to construct F, G and H (or nothing)
- `coordinates_transition_P0`: index describing which element of P0 should be re-estimated in the cm-step
"""
function ksmoother_ecm!(estim::EstimSettings, sspace::KalmanSettings, status::SizedKalmanStatus, smoother_arrays::SmootherArrays, coordinates_measurement_states::Union{IntVector, Nothing}, coordinates_transition_current::Union{IntVector, Nothing}, coordinates_transition_lagged::Union{IntVector, Nothing}, coordinates_transition_PPs::Union{IntVector, Nothing}, coordinates_transition_P0::IntVector)

    # Loop over t (from status.t-1 to 1)
    for t=status.online_status.t:-1:1

        # Kalman smoother iteration
        Xs, Ps, ind_not_missings = ksmoother_ecm_iteration!(sspace, status, smoother_arrays, t);

        # Update ECM statistics for the measurement equation
        call_update_ecm_stats_measurement!(estim, smoother_arrays, coordinates_measurement_states, ind_not_missings, t, Xs, Ps);

        #=
        Update ECM statistics for the transition equation - please note that:
            1. the following if-statement makes sure to have a valid value for the *_leading estimators;
            2. implementing the update this way is beneficial from a computational standpoint, but `update_ecm_stats_transition!` needs to be called an additional time outside the loop to finalise the calculations.
        =#

        if t < status.online_status.t
            call_update_ecm_stats_transition!(estim, smoother_arrays, Xs, Ps, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);
        end

        # Update Xs_leading and Ps_leading
        copyto!(smoother_arrays.Xs_leading, Xs);
        copyto!(smoother_arrays.Ps_leading.data, Ps.data);
    end

    # Compute smoothed estimates for t==0
    call_update_smoothing_factors!(sspace, status, smoother_arrays, nothing);

    # CM-step for X0 and P0
    cm_step_X0_P0!(sspace, status, smoother_arrays, coordinates_transition_P0);

    # Finalise calculations of the ECM statistics for the transition equation
    call_update_ecm_stats_transition!(estim, smoother_arrays, sspace.X0, sspace.P0, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);
end

#=
--------------------------------------------------------------------------------------------------------------------------------
CM-step
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    turbo_dot(x::AbstractVector{Float64}, y::AbstractVector{Float64})

Dot product.

    turbo_dot(x::AbstractVector{Float64}, A::AbstractMatrix{Float64}, y::AbstractVector{Float64})

Generalised dot product.
"""
@inline function turbo_dot(x::AbstractVector{Float64}, y::AbstractVector{Float64})
    s = 0.0;
    @turbo for i in eachindex(x)
        s += x[i]*y[i];
    end
    return s;
end

@inline function turbo_dot(x::AbstractVector{Float64}, A::AbstractMatrix{Float64}, y::AbstractVector{Float64})
    s = 0.0;
    for j in eachindex(y)
        @inbounds yj = y[j];
        temp = 0.0;
        if !iszero(yj)
            @turbo for i in eachindex(x)
                temp += x[i]*A[i,j];
            end
            s += temp*yj;
        end
    end
    return s;
end

"""
    cm_step_time_loop(sspace::KalmanSettings, B_star::SubArray{Float64}, ij::CartesianIndex{2}, N::Vector{SparseMatrixCSC{Float64, Int64}}, O::Array{VectorsArray{Float64},1})

Compute numerator and denominator to compute the (i,j)-th CM update for the measurement equation coefficients in `cm_step!`.
"""
function cm_step_time_loop(sspace::KalmanSettings, B_star::SubArray{Float64}, ij::CartesianIndex{2}, N::Vector{SparseMatrixCSC{Float64, Int64}}, O::Array{VectorsArray{Float64},1})

    # Coordinates
    i,j = ij.I;

    # Store transpose of B_star in Array format and B_{ij} to speed-up calculations
    B_star_transpose = B_star' |> FloatMatrix;
    B_star_ij = B_star[ij];

    # Memory pre-allocation for ij sums
    sum_numerator = 0.0;
    sum_denominator = 0.0;

    for t=1:sspace.Y.T

        # Proceed if N_t and O_t are not #undef (the points in time when at least one series is observed)
        if isassigned(N, t) && isassigned(O, t)

            # Pointers
            N_t = N[t];
            O_t = O[t];
            O_jt = O_t[j];

            # Convenient view into N_t
            N_it = @view N_t[:,i];

            # Shortcut
            @inbounds NO_ij_t = N_it[i]*O_jt[j];

            # Update numerator
            sum_numerator += turbo_dot(O_jt, B_star_transpose, N_it); # this order gives a faster performance, since N_it includes a significant number of zeros!
            sum_numerator -= NO_ij_t*B_star_ij;

            # Update `sum_denominator`
            sum_denominator += NO_ij_t;
        end
    end

    return sum_numerator, sum_denominator, i, j;
end

"""
    call_cm_step!(estim::EstimSettings, sspace::KalmanSettings, smoother_arrays::SmootherArrays, B_star::SubArray{Float64}, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::CoordinatesVector)
    call_cm_step!(estim::EstimSettings, sspace::KalmanSettings, smoother_arrays::SmootherArrays, B_star::SubArray{Float64}, C_star::Nothing, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::Nothing)
    call_cm_step!(estim::EstimSettings, sspace::KalmanSettings, smoother_arrays::SmootherArrays, B_star::Nothing, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::Nothing, coordinates_free_params_C::CoordinatesVector)

APIs to call `cm_step!` with SmootherArrays.
"""
call_cm_step!(estim::EstimSettings, sspace::KalmanSettings, smoother_arrays::SmootherArrays, B_star::SubArray{Float64}, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::CoordinatesVector) = cm_step!(estim, sspace, B_star, C_star, Q_view, smoother_arrays.F, smoother_arrays.G, smoother_arrays.H, smoother_arrays.M, smoother_arrays.N, smoother_arrays.O, coordinates_free_params_B, coordinates_free_params_C);
call_cm_step!(estim::EstimSettings, sspace::KalmanSettings, smoother_arrays::SmootherArrays, B_star::SubArray{Float64}, C_star::Nothing, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::Nothing) = cm_step!(estim, sspace, B_star, C_star, Q_view, smoother_arrays.F, smoother_arrays.G, smoother_arrays.H, smoother_arrays.M, smoother_arrays.N, smoother_arrays.O, coordinates_free_params_B, coordinates_free_params_C);
call_cm_step!(estim::EstimSettings, sspace::KalmanSettings, smoother_arrays::SmootherArrays, B_star::Nothing, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::Nothing, coordinates_free_params_C::CoordinatesVector) = cm_step!(estim, sspace, B_star, C_star, Q_view, smoother_arrays.F, smoother_arrays.G, smoother_arrays.H, smoother_arrays.M, smoother_arrays.N, smoother_arrays.O, coordinates_free_params_B, coordinates_free_params_C);

"""
    enforce_causality_and_invertibility_write!(companion_new::FloatMatrix, companion_proposal::FloatMatrix, params_new::Nothing)
    enforce_causality_and_invertibility_write!(companion_new::AbstractArray{Float64,2}, companion_proposal::FloatMatrix, params_new::FloatMatrix)

Barrier function for `enforce_causality_and_invertibility!` based on @turbo. Replace inplace either `companion_new` or `params_new` depending on whether the latter is `nothing` or not.

    enforce_causality_and_invertibility_write!(companion_new::AbstractArray{Float64,2}, companion_proposal::AbstractArray{Float64,2}, params_new::Nothing)
    enforce_causality_and_invertibility_write!(companion_new::AbstractArray{Float64,2}, companion_proposal::AbstractArray{Float64,2}, params_new::AbstractArray{Float64,2})

Barrier function for `enforce_causality_and_invertibility!`. Replace inplace either `companion_new` or `params_new` depending on whether the latter is `nothing` or not.
"""
@inline function enforce_causality_and_invertibility_write!(companion_new::FloatMatrix, companion_proposal::FloatMatrix, params_new::Nothing)
    @turbo for ij in eachindex(companion_new)
        companion_new[ij] = companion_proposal[ij];
    end
end

@inline function enforce_causality_and_invertibility_write!(companion_new::AbstractArray{Float64,2}, companion_proposal::FloatMatrix, params_new::FloatMatrix)
    @turbo for j in axes(params_new, 2), i in axes(params_new, 1)
        params_new[i,j] = companion_proposal[i,j];
    end
end

function enforce_causality_and_invertibility_write!(companion_new::AbstractArray{Float64,2}, companion_proposal::AbstractArray{Float64,2}, params_new::Nothing)
    for ij in eachindex(companion_new)
        @inbounds companion_new[ij] = companion_proposal[ij];
    end
end

function enforce_causality_and_invertibility_write!(companion_new::AbstractArray{Float64,2}, companion_proposal::AbstractArray{Float64,2}, params_new::AbstractArray{Float64,2})
    for j in axes(params_new, 2), i in axes(params_new, 1)
        @inbounds params_new[i,j] = companion_proposal[i,j];
    end
end

"""
    enforce_causality_and_invertibility!(companion_new::AbstractArray{Float64,2}, companion_old::AbstractArray{Float64,2}, params_new::Union{Nothing, AbstractArray{Float64,2}}=nothing)

Shrink `companion_new` towards `companion_old` (when needed) to ensure it is causal / invertible (via grid search). If `params_new` is not `nothing` the function output replaces it inplace instead.

    enforce_causality_and_invertibility!(params_new::AbstractArray{Float64,2})

Compute the companion form of `params_new` and shrink it towards zero (when needed) to ensure it is causal / invertible (via grid search).
"""
function enforce_causality_and_invertibility!(companion_new::AbstractArray{Float64,2}, companion_old::AbstractArray{Float64,2}, params_new::Union{Nothing, AbstractArray{Float64,2}}=nothing)

    # Loop over grid of weights
    for weight=1.0:-0.1:0.0

        # Compute weighted cycle
        companion_proposal = (weight*companion_new + (1-weight)*companion_old)::FloatMatrix;

        # Compute corresponding largest eigenvalue
        ith_eigmax = maximum(abs.(eigvals(companion_proposal)));

        # Keep current setting and adjust
        if (ith_eigmax <= 0.98) && (weight != 1)
            enforce_causality_and_invertibility_write!(companion_new, companion_proposal, params_new);
            break;

        # No adjustment needed
        elseif (ith_eigmax <= 0.98) && (weight == 1)
            break;
        end
    end
end

function enforce_causality_and_invertibility!(params_new::AbstractArray{Float64,2})
    companion_old = companion_form(zeros(size(params_new)), extended=false);
    companion_new = companion_form(params_new, extended=false);
    enforce_causality_and_invertibility!(companion_new, companion_old, params_new);
end

#=
--------------------------------------------------------------------------------------------------------------------------------
Optimisation
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    sspace_representation(estim::EstimSettings, B_star::FloatMatrix, C_star::FloatMatrix, Q::SymMatrix; extended=true)

Return the default state-space coefficients. If extended=true, `n` additional states at the end for the PPs calculation proposed in Watson and Engle (1983).
"""
function sspace_representation(estim::EstimSettings, B_star::FloatMatrix, C_star::FloatMatrix, Q::SymMatrix; extended=true)

    # Observation equation (fixed) coefficients
    B = [Matrix(1.0I, estim.n, estim.n) B_star zeros(estim.n*extended, estim.n*extended)];
    R = Symmetric(Matrix(estim.ε * I, estim.n, estim.n));

    # Transition equation
    C = companion_form(C_star, extended=extended);
    D = zeros(estim.m+estim.n*extended, estim.n);
    D[1:estim.n, :] = Matrix(1.0I, estim.n, estim.n);

    # Return output
    return B, R, C, D, Q;
end

"""
    update_sspace_data!(sspace::KalmanSettings, output_sspace_data::Union{FloatMatrix, JMatrix{Float64}})

Update `sspace.Y` with `output_sspace_data`.

    update_sspace_data!(sspace::KalmanSettings, output_sspace_data::Nothing)

Return `nothing`.
"""
function update_sspace_data!(sspace::KalmanSettings, output_sspace_data::Union{FloatMatrix, JMatrix{Float64}})
    sspace.Y.data = output_sspace_data;
    sspace.Y.T = size(output_sspace_data, 2);
end

update_sspace_data!(sspace::KalmanSettings, output_sspace_data::Nothing) = nothing;

"""
    ecm(estim::EstimSettings)

Run the ECM algorithm in Pellegrino (2022).

# Arguments
- `estim`: settings used for the estimation
- `output_sspace_data`: Optional argument. If specified, it is used as the output state space data. Otherwise, estim.Y is used instead.

# References
Pellegrino (2022)
"""
function ecm(estim::EstimSettings; output_sspace_data::Union{FloatMatrix, JMatrix{Float64}, Nothing}=nothing)

    # Check inputs
    check_model_bounds(estim);

    #=
    The state vector includes additional terms with respect to standard companion form representations.
    This is to estimate the lag-one covariance smoother as in Watson and Engle (1983).
    =#

    # Initialise state-space
    verb_message(estim.verb, "ecm > initialisation");
    sspace, B_star, C_star, Q_view, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs, coordinates_transition_P0, coordinates_free_params_B, coordinates_free_params_C = initialise(estim);

    # Pre-allocate memory for Kalman filter and smoother
    status = SizedKalmanStatus(sspace.Y.T);
    smoother_arrays = SmootherArrays(estim, sspace, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);

    # ECM controls
    ecm_params_old = vec_ecm_params(estim, B_star, C_star, Q_view, coordinates_free_params_B, coordinates_free_params_C);
    ecm_params_new = similar(ecm_params_old);
    params_rel_change = similar(ecm_params_old);

    # Run ECM
    for iter=1:estim.max_iter

        # Run Kalman filter
        kfilter_full_sample!(sspace, status);

        # Compute Kalman stats
        ksmoother_ecm!(estim, sspace, status, smoother_arrays, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs, coordinates_transition_P0);

        # CM-step
        call_cm_step!(estim, sspace, smoother_arrays, B_star, C_star, Q_view, coordinates_free_params_B, coordinates_free_params_C);

        if iter > estim.prerun

            # Update ecm_params_new
            vec_ecm_params!(estim, ecm_params_new, B_star, C_star, Q_view, coordinates_free_params_B, coordinates_free_params_C);

            # Relative change criterion
            update_relative_change!(estim, params_rel_change, ecm_params_new, ecm_params_old);

            # Stop when the ECM algorithm converges
            if iter > estim.prerun+1
                if ecm_converged(estim, params_rel_change, iter)
                    verb_message(estim.verb, "ecm > converged!\n");
                    break;
                end
            end

            # Store current run information
            copyto!(ecm_params_old, ecm_params_new);

        else
            verb_message(estim.verb, "ecm > prerun $iter (out of $(estim.prerun))");
        end

        reinitialise_ecm_stats!(smoother_arrays);
    end

    # Final operations
    update_sspace_data!(sspace, output_sspace_data); # TBD: Use compress representation
    sspace.Q.data .= Array(Q_view); # update Q (for internal consistency)

    # Return output
    return sspace;
end


text:
function jacobian(f::Function, x::Vector; δ = 0.001)

    y = f(x)

    J = Matrix{Union{eltype(y), Missing}}(undef, length(y), length(x))

    for i in 1:length(x)
        if ismissing(x[i])
            J[:, i] .= missing
        else
            dx = copy(x)
            dx[i] += δ
            J[:, i] .= (y .- f(dx)) ./ δ
        end
    end

    return J
end


text:
# Object-based wrapper for level-1 BLIS matrix routines.
#

macro blis_object_api_level1m_form1(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           a, xObj,
                           b, xObj)
    end
end

@blis_object_api_level1m_form1 addm
@blis_object_api_level1m_form1 copym
@blis_object_api_level1m_form1 subm

macro blis_object_api_level1m_form2(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           α, xObj,
                           a, xObj,
                           b, xObj)
    end
end

@blis_object_api_level1m_form2 axpym
@blis_object_api_level1m_form2 scal2m

macro blis_object_api_level1m_form3(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           α, xObj,
                           a, xObj)
    end
end

@blis_object_api_level1m_form3 scalm
@blis_object_api_level1m_form3 setm

macro blis_object_api_level1m_form4(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           a, xObj,
                           β, xObj,
                           b, xObj)
    end
end

@blis_object_api_level1m_form4 xpbym
@blis_object_api_level1m_form4 xpbym_md


text:
begin
	using StatsBase
	function coin_toss(n)
		toss = []
		head = 0
		for i in 1:n
			push!(toss, sample(["H", "T"], 1, replace = false))
		end
		for i in toss
			if i == ["H"]
				head += 1
			end
		end
		return head/n
	end
end


text:
struct Material
    id::String
    hid::Int
    mat_type::MaterialType

    E::Float64
    ν::Float64
    ρ::Float64
    α::Float64
    G::Float64
    λ::Float64

    E₂::Float64
    f::Float64
    fᵤ::Float64
end

function IsoElastic(id,hid,E,ν,ρ,α)
    G=E/2/(1+ν)
    λ=E*ν/(1+ν)/(1-2ν)
    mat_type=Enums.ISOELASTIC
    Material(string(id),hid,mat_type,E,ν,ρ,α,G,λ,0,0,0)
end

function UniAxialMetal(id,hid,E,ν,ρ,α,E₂,f,fᵤ)
    G=E/2/(1+ν)
    λ=E*ν/(1+ν)/(1-2ν)
    mat_type=Enums.UNIAXIAL_METAL
    Material(id,hid,mat_type,E,ν,ρ,α,G,λ,E₂,f,fᵤ)
end

# end


text:
# simulate circle multimodal vs outlier example

# using Revise

using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--altmode"
            help = "fraction for alternative mode"
            arg_type = Float64
            default = 0.7
        "--resolve"
            help = "should a repeat solve be performed in the loop"
            action = :store_true
        "--CYCLES"
            help = "How many cycles to drive"
            arg_type = Int
            default = 10
        "--spreadNH"
            help = "Scale adjustment for spreading nullhypo (entropy)"
            arg_type = Float64
            default = 3.0
    end

    return parse_args(s)
end

pargs = parse_commandline()

using Distributed
# addprocs(8) # use -p8 instead

using Random

using RoME, DistributedFactorGraphs
using RoMEPlotting
using KernelDensityEstimatePlotting
@everywhere using RoME

using FileIO
using JLD2
using JSON2
using Gadfly
Gadfly.set_default_plot_size(35cm,25cm)

include(joinpath(@__DIR__, "FixedLagCirclePlotting.jl"))

##  Landmark ground truth

# assume all 10 landmarks changing position from rows 1:10 to 11:20, randomly after each circle
L = [0.0 20.0;      # L1
     30.0 30.0;     # L2
     -10.0 20.0;    # L3
     20.0 30.0;     # L4
     -20.0 30.0;    # L5
     -20.0 10.0;    # L6
     20.0 -0.0;     # L7
     10.0 10.0;     # L8
     0.0 40.0;      # L9
     30.0 10.0;     # L10
     10.0 -10.0;    # L1_0
     10.0 40.0;     # L2_0
     -10.0 -10.0;   # L3_0
     -10.0 30.0;    # L4_0
     30.0 -0.0;     # L5_0
     20.0 20.0;     # L6_0
     0.0 10.0;      # L7_0
     30.0 20.0;     # L8_0
     20.0 -10.0;    # L9_0
     -10.0 10.0]    # L10_0

# working from previous save stored right here in the example folder
# @load string(joinpath(@__DIR__, "LassoBearingRange.jld2")) BR

SIZE = 10
CYCLES = pargs["CYCLES"]
LANDMARKS = 10
altFrac = pargs["altmode"]
mainFrac = 1-pargs["altmode"]

## divise landmark rotation schedule

movePerCycle = rand(Categorical([0.6;0.4]),CYCLES)
idx = sum(cumsum(movePerCycle) .< CYCLES)
movePerCycle = movePerCycle[1:idx]
push!(movePerCycle, CYCLES-cumsum(movePerCycle)[end])
order = shuffle(1:CYCLES)
sequence = Vector{Vector{Int}}()
total = 1
for n in 1:length(movePerCycle)
  global total
  push!(sequence, order[total:total+movePerCycle[n]-1] )
  total += movePerCycle[n]
end
if pargs["CYCLES"] == 10
  @assert sum(union(sequence...)) == 55
end

sequence

## convert sequence to landmark lookups

lookup = Vector{Vector{Symbol}}()
push!(lookup, [Symbol("l$i") for i in 1:10])

for i in 1:9
  lu = deepcopy(lookup[i])
  # before each loop, change location of some landmarks
  seq = i <= length(sequence) ? sequence[i] : Int[]
  for s in seq
    lu[s] = Symbol("l$(s)_0")
  end

  # push new vector into lookup
  push!(lookup, lu)
end

# this is the sequence of landmark transitions
lookup

# BRdistr[:l1]
# BRdistr[:l1_0]

### SETUP first circle----------------------------------------------

# drive first circle
fg = generateCanonicalFG_Circle(SIZE, kappaOdo=0.1, loopClosure=false, landmark=false, cyclePoses=10)

getSolverParams(fg).spreadNH = pargs["spreadNH"]

ensureAllInitialized!(fg)

## store arguments in results log

argstr = JSON2.write(pargs)

Base.mkpath(getLogPath(fg))
fid = open(joinLogPath(fg, "args.json"), "w")
JSON2.write(fid, pargs)
close(fid)
fid = open(joinLogPath(fg, "..", "results.log"), "a")
Base.write(fid, "$(getLogPath(fg)), RoME/examples/OutlierVsMultimodal.jl, ")
JSON2.write(fid, pargs)
Base.write(fid, "\n")
close(fid)

## add initial landmarks

BR = Dict{Symbol,Dict{Symbol, Tuple}}()
# add the landmarks
for i in 1:2*CYCLES
  lmid = Symbol(i <= 10 ? "l$i" : "l$(i-10)_0")
  # lmid = Symbol("l$i")
  addVariable!(fg, lmid, RoME.Point2)
  pts = rand(MvNormal(L[i,:],diagm([0.01;0.01].^2)),100)
  initManual!(fg, lmid, manikde!(pts, RoME.Point2))
  setVariablePosteriorEstimates!(fg,lmid)

  BR[lmid] = Dict{Symbol,Tuple}()
  varNear, varDist = findVariablesNear(fg, L[i,:], r"x\d", number=10)
  for vsym in varNear[varDist .< 25]
    # calculate bearing and range factor
    bear, rang = predictBodyBR(fg, vsym, lmid)
    if abs(bear) < pi/3 && vsym != :x10
      BR[lmid][vsym] = (bear, rang)
    end
  end
end
for i in CYCLES+1:2*CYCLES
  lmid = Symbol(i <= 10 ? "l$i" : "l$(i-10)_0")
  # vsym = Symbol("l$i")
  setSolvable!(fg, lmid, 0)
  # ( x->deleteFactor!(fg, x) ).( ls(fg, vsym) )  # no factors added yet
  # deleteVariable!(fg, lmid)
end
# @save joinpath(@__DIR__, "LassoBearingRange.jld2") BR

## see what is going on

plfl1 = drawPosesLandms(fg, spscale=1.0,title=getLogPath(fg)*"\n$argstr")

## prepare factors to use

BRdistr = Dict{Symbol,Dict{Symbol,Tuple}}()
for i in 1:2*SIZE
  lmid = Symbol(i <= SIZE ? "l$i" : "l$(i-SIZE)_0")
  BRdistr[lmid] = Dict{Symbol,Tuple}()
  for (vsym, br) in BR[lmid]
    @show bear, rang = br
    BRdistr[lmid][vsym] = (Normal(bear,0.05),Normal(rang,0.3))
  end
end

BRdistr

## add factors to the landmarks on first cycle

@assert ls(fg, :l1) |> length == 0

OFFSET=0
for id in 1:CYCLES
  lmid = Symbol("l$(id)")
  lmid_0 = Symbol("l$(id)_0")
  for (psid, ppbr) in BRdistr[lmid]
    # global OFFSET
    adjpsid = psid # add offset
    addFactor!(fg, [adjpsid;lmid], Pose2Point2BearingRange(ppbr...))
  end
end

##

getSolverParams(fg).drawtree = true
tree, smt, hist = solveTree!(fg)

plfl1 = drawPosesLandms(fg, spscale=1.0,title=getLogPath(fg)*"\n$argstr")

# lookup
##  change to alternate location, one at a time

# activate solvable on alternative landmark locations for remainer of test
for l in 1:10
  lmid_0 = Symbol("l$(l)_0")
  setSolvable!(fg, lmid_0, 1)
end

POSEOFFSET=0

## store

@save joinLogPath(fg, "lookup.jld2") lookup

plfl1 = drawPosesLandms(fg, spscale=1.0,title=getLogPath(fg)*"\n$argstr", landmsPPE=:max, contour=true) #, posesPPE=:max)
plfl1 |> PDF(joinLogPath(fg, "plot_x20_before.pdf"), 20cm, 17cm)

##

# Drive CYCLES-1 more loops
i = 2
for i in 2:CYCLES-1

global fg
global lookup
global SIZE
global BRdistr
global pargs
global POSEOFFSET += 10
# drive the new loop without landmark detections (dont solve yet)
fg = generateCanonicalFG_Circle(i*SIZE, fg=fg, kappaOdo=0.1, loopClosure=false, landmark=false, cyclePoses=10)

# add modified landmark sighting measurements
for l in 1:LANDMARKS
  lmidBr = lookup[i][l]
  lmid = Symbol("l$(l)")
  lmid_0 = Symbol("l$(l)_0")
  for (psid, ppbr) in BRdistr[lmidBr]
    # apply pose offset
    psnum = parse(Int, string(psid)[2:end]) + POSEOFFSET
    opsid = Symbol("x$psnum")
    # add the few BR factors from this lmid accordingly
    addFactor!(fg, [opsid;lmid;lmid_0], Pose2Point2BearingRange(ppbr...), multihypo=[1;mainFrac;altFrac])
  end

  # reinit all lm_0 with the latest factor
  if 0 < length(ls(fg, lmid_0))
    pts = approxConv(fg, ls(fg, lmid_0)[1], lmid_0  )
    initManual!(fg,lmid_0,pts)
  end
end

getSolverParams(fg).dbg = true
saveDFG(fg, joinLogPath(fg, "fg_x$(POSEOFFSET+10)before"))
getSolverParams(fg).maxincidence = 1000
tree, smt, hist = solveTree!(fg, recordcliqs=ls(fg))
saveDFG(fg, joinLogPath(fg, "fg_x$(POSEOFFSET+10)_solve"))

plfl1 = drawPosesLandms(fg, spscale=1.0,title=getLogPath(fg)*"\n$argstr", landmsPPE=:max, contour=true)
plfl1 |> PDF(joinLogPath(fg, "plot_x$(POSEOFFSET+10)_solve.pdf"), 20cm, 17cm)

if pargs["resolve"]
  tree, smt, hist = solveTree!(fg, recordcliqs=ls(fg))
  saveDFG(fg, joinLogPath(fg, "fg_x$(POSEOFFSET+10)_resolve"))
  plfl1 = drawPosesLandms(fg, spscale=1.0,title=getLogPath(fg)*"\n$argstr", landmsPPE=:max, contour=true)
  plfl1 |> PDF(joinLogPath(fg, "plot_x$(POSEOFFSET+10)_resolve.pdf"), 20cm, 17cm)
end

end

## test solve

getSolverParams(fg).dbg = true
getSolverParams(fg).maxincidence = 1000
tree, smt, hist = solveTree!(fg, recordcliqs=ls(fg)) #[:l9_0;:x14])
saveDFG(fg, joinLogPath(fg, "fg_x$(POSEOFFSET+10)_final"))

##

plfl1 = drawPosesLandms(fg, spscale=1.0,title=getLogPath(fg)*"\n$argstr", landmsPPE=:max)
plfl1 |> PDF(joinLogPath(fg, "plot_x$(POSEOFFSET+10)_final.pdf"), 20cm, 17cm)

#


text:
# NOTE: All tests of this file can be run with any number of processes.
push!(LOAD_PATH, "../src")
using Test
using ImplicitGlobalGrid; GG = ImplicitGlobalGrid
import MPI
using CUDA
import ImplicitGlobalGrid: @require

test_gpu = CUDA.functional()

## Test setup
MPI.Init();
nprocs = MPI.Comm_size(MPI.COMM_WORLD); # NOTE: these tests can run with any number of processes.

@testset "$(basename(@__FILE__)) (processes: $nprocs)" begin
    @testset "1. select_device" begin
        me, = init_global_grid(3, 4, 5; quiet=true, init_MPI=false);
        @static if test_gpu
            gpu_id = select_device();
            @test gpu_id < length(CUDA.devices())
        else
            @test_throws ErrorException select_device()
        end
        finalize_global_grid(finalize_MPI=false);
    end;
end;

## Test tear down
MPI.Finalize()


text:
if basename(pwd()) == "aoc"
    cd("2018/2")
end
using StatsBase: countmap

function part1(filename::AbstractString)
    hastwo = 0
    hasthree = 0
    for line in eachline(filename)
        counts = values(countmap(line))
        hastwo += any(==(2), values(counts))
        hasthree += any(==(3), values(counts))
    end
    hastwo * hasthree
end
@assert part1("example.txt") == 12
part1("input.txt")

function part2(filename::AbstractString)
    lines = readlines(filename)
    pairs = [(l, r) for (l, r) in Iterators.product(lines, lines) if l != r]
    l, r = pairs[findmin(hammingdistance, pairs)[2]]
    join(a for (a, b) in zip(l, r) if a == b)
end

hammingdistance((l, r)) = count(a != b for (a, b) in zip(l, r))

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
export ThinCategory, FreeThinCategory,
  ThinSymmetricMonoidalCategory, FreeThinSymmetricMonoidalCategory,
  Preorder, FreePreorder, Elt, Leq, ≤, lhs, rhs, reflexive, transitive

# Thin category
###############

""" Theory of *thin categories*

Thin categories have at most one morphism between any two objects.
"""
@theory Category(Ob,Hom) => ThinCategory(Ob,Hom) begin
  f == g ⊣ (A::Ob, B::Ob, f::Hom(A,B), g::Hom(A,B))
end

@syntax FreeThinCategory(ObExpr,HomExpr) ThinCategory begin
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))
end

@theory SymmetricMonoidalCategory(Ob,Hom) => ThinSymmetricMonoidalCategory(Ob,Hom) begin
  f == g ⊣ (A::Ob, B::Ob, f::Hom(A,B), g::Hom(A,B))
end

@syntax FreeThinSymmetricMonoidalCategory(ObExpr,HomExpr) ThinSymmetricMonoidalCategory begin
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))
  otimes(A::Ob, B::Ob) = associate_unit(new(A,B), munit)
  otimes(f::Hom, g::Hom) = associate(new(f,g))
end

# Preorder
##########

""" Theory of *preorders*

Preorders encode the axioms of reflexivity and transitivity as term constructors.
"""
@theory Preorder(Elt,Leq) begin
  Elt::TYPE
  Leq(lhs::Elt, rhs::Elt)::TYPE
  @op (≤) := Leq

  # Preorder axioms are lifted to term constructors in the GAT.
  reflexive(A::Elt)::(A≤A) # ∀ A there is a term reflexive(A) which implies A≤A
  transitive(f::(A≤B), g::(B≤C))::(A≤C) ⊣ (A::Elt, B::Elt, C::Elt)

  # Axioms of the GAT are equivalences on terms or simplification rules in the logic
  f == g ⊣ (A::Elt, B::Elt, f::(A≤B), g::(A≤B))
  # Read as (f⟹ A≤B ∧ g⟹ A≤B) ⟹ f ≡ g
end

@syntax FreePreorder(ObExpr,HomExpr) Preorder begin
  transitive(f::Leq, g::Leq) = associate(new(f,g; strict=true))
end

# TODO: a GAT-homomorphism between the Preorder GAT and the ThinCategory GAT
# this is a morphism is *GAT* the category whose objects are GATs
# and whose morphisms are algebraic structure preserving maps

# @functor F(Preorder(Elt, Leq))::ThinCategory(Ob, Hom) begin
#   Elt ↦ Ob
#   Leq ↦ Hom
#   reflexive ↦ id
#   transitive ↦ compose
# end


text:
module Flase

export InfiniteSimulation, FiniteSimulation, runsim
export BrownianMotion, ConstVelocity
export UnicodePlotter
export World
export DenseSheeps

include("motions/Motion.jl")
include("motions/ConstVelocity.jl")
include("motions/BrownianMotion.jl")

include("collectors/Dog.jl")
include("collectors/Dogs.jl")
include("items/Sheeps.jl")
include("items/DenseSheeps.jl")

include("World.jl")

include("item-collector-interactions.jl")

include("plotter/Plotter.jl")
include("plotter/VoidPlotter.jl")
include("plotter/UnicodePlotter.jl")

include("simulations/Simulation.jl")
include("simulations/InfiniteSimulation.jl")
include("simulations/FiniteSimulation.jl")

include("interactions.jl")

end # module


text:
@testset "Test symmetric_analytical_jacobian_position for point vortices with real freestream" begin

    atol = 1000*eps()

    Nv = 10
    Nx = 3*Nv

    x = rand(Nx)
    x0 = deepcopy(x)

    xsensors = collect(-2.0:0.5:10)
    sensors = complex(xsensors)
    Ny = length(sensors)
#     U = 0.0 + 0.0*im
    U = randn() + 0.0*im

    #     U = randn(ComplexF64)
    freestream = Freestream(U);

    config = let Nv = Nv,
             U = U,
             ss = sensors, Δt = 1e-3, δ = 5e-2,
             ϵX = 1e-3, ϵΓ = 1e-4,
             β = 1.0,
             ϵY = 1e-16
            VortexConfig(Nv, U, ss, Δt, δ, ϵX, ϵΓ, β, ϵY)
    end

    sys = state_to_lagrange(x, config; isblob = false)
    @test typeof(sys[1][1])<:Vortex.Point{Float64, Float64}

    sys = vcat(sys...)

    dpdz, dpdzstar = analytical_jacobian_position(sensors, sys, freestream, 0.0)

    dpdzsym, dpdzstarsym = symmetric_analytical_jacobian_position(xsensors, sys, freestream, 0.0)

    @test isapprox(dpdz, dpdzsym, atol = atol)
    @test isapprox(dpdzstar, dpdzstarsym, atol = atol)
end

@testset "Test symmetric_analytical_jacobian_position for regularized vortices with real freestream" begin

    atol = 1000*eps()

    Nv = 10
    Nx = 3*Nv

    x = rand(Nx)
    x0 = deepcopy(x)

    xsensors = collect(-2.0:0.5:10)
    sensors = complex(xsensors)
    Ny = length(sensors)
#     U = 0.0 + 0.0*im
    U = randn() + 0.0*im

    #     U = randn(ComplexF64)
    freestream = Freestream(U);

    config = let Nv = Nv,
             U = U,
             ss = sensors, Δt = 1e-3, δ = 5e-2,
             ϵX = 1e-3, ϵΓ = 1e-4,
             β = 1.0,
             ϵY = 1e-16
            VortexConfig(Nv, U, ss, Δt, δ, ϵX, ϵΓ, β, ϵY)
    end

    sys = state_to_lagrange(x, config; isblob = true)
    @test typeof(sys[1][1])<:Vortex.Blob{Float64, Float64}

    sys = vcat(sys...)

    dpdz, dpdzstar = analytical_jacobian_position(sensors, sys, freestream, 0.0)

    dpdzsym, dpdzstarsym = symmetric_analytical_jacobian_position(xsensors, sys, freestream, 0.0)

    @test isapprox(dpdz, dpdzsym, atol = atol)
    @test isapprox(dpdzstar, dpdzstarsym, atol = atol)
end

@testset "Test symmetric_analytical_jacobian_strength for point vortices with real freestream" begin

    atol = 1000*eps()

    Nv = 10
    Nx = 3*Nv

    x = rand(Nx)
    x0 = deepcopy(x)

    xsensors = collect(-2.0:0.5:10)
    sensors = complex(xsensors)
    Ny = length(sensors)

    U = randn() + 0.0*im
    freestream = Freestream(U);

    config = let Nv = Nv,
             U = U,
             ss = sensors, Δt = 1e-3, δ = 5e-2,
             ϵX = 1e-3, ϵΓ = 1e-4,
             β = 1.0,
             ϵY = 1e-16
            VortexConfig(Nv, U, ss, Δt, δ, ϵX, ϵΓ, β, ϵY)
    end

    sys = state_to_lagrange(x, config; isblob = false)
    @test typeof(sys[1][1])<:Vortex.Point{Float64, Float64}

    sys = vcat(sys...)

    dpdS, dpdSstar = analytical_jacobian_strength(sensors, sys, freestream, 0.0)

    dpdSsym, dpdSstarsym = symmetric_analytical_jacobian_strength(xsensors, sys, freestream, 0.0)

    @test isapprox(dpdS, dpdSsym, atol = atol)
    @test isapprox(dpdSstar, dpdSstarsym, atol = atol)
end

@testset "Test symmetric_analytical_jacobian_strength for regularized vortices with real freestream" begin

    atol = 1000*eps()

    Nv = 10
    Nx = 3*Nv

    x = rand(Nx)
    x0 = deepcopy(x)

    xsensors = collect(-2.0:0.5:10)
    sensors = complex(xsensors)
    Ny = length(sensors)

    U = randn() + 0.0*im
    freestream = Freestream(U);

    config = let Nv = Nv,
             U = U,
             ss = sensors, Δt = 1e-3, δ = 5e-2,
             ϵX = 1e-3, ϵΓ = 1e-4,
             β = 1.0,
             ϵY = 1e-16
            VortexConfig(Nv, U, ss, Δt, δ, ϵX, ϵΓ, β, ϵY)
    end

    sys = state_to_lagrange(x, config; isblob = true)
    @test typeof(sys[1][1])<:Vortex.Blob{Float64, Float64}

    sys = vcat(sys...)

    dpdS, dpdSstar = analytical_jacobian_strength(sensors, sys, freestream, 0.0)

    dpdSsym, dpdSstarsym = symmetric_analytical_jacobian_strength(xsensors, sys, freestream, 0.0)

    @test isapprox(dpdS, dpdSsym, atol = atol)
    @test isapprox(dpdSstar, dpdSstarsym, atol = atol)
end

@testset "Test routine symmetric_analytical_jacobian_pressure with real freestream" begin
    atol = 1000*eps()

    Nv = 10
    Nx = 3*Nv

    x = rand(Nx)
    x0 = deepcopy(x)

    xsensors = collect(-2.0:0.5:10)
    sensors = complex(xsensors)
    Ny = length(sensors)

    U = randn() + 0.0*im
    freestream = Freestream(U);

    config = let Nv = Nv,
             U = U,
             ss = sensors, Δt = 1e-3, δ = 5e-2,
             ϵX = 1e-3, ϵΓ = 1e-4,
             β = 1.0,
             ϵY = 1e-16
            VortexConfig(Nv, U, ss, Δt, δ, ϵX, ϵΓ, β, ϵY)
    end

    sys = state_to_lagrange(x, config; isblob = true)
    @test typeof(sys[1][1])<:Vortex.Blob{Float64, Float64}

    sys = vcat(sys...)

    Jfull = analytical_jacobian_pressure(sensors, sys, freestream, 0.0)
    Jsymfull = symmetric_analytical_jacobian_pressure(xsensors, sys, freestream, 0.0)

    @test isapprox(Jfull, Jsymfull, atol = atol)

    # Test only on a subset of the vortices

    J = analytical_jacobian_pressure(sensors, sys, freestream, 1:Nv, 0.0)
    Jsym = symmetric_analytical_jacobian_pressure(xsensors, sys, freestream, 1:Nv, 0.0)

    @test isapprox(J[:,1:3*Nv], Jfull[:,1:3*Nv], atol = atol)
    @test isapprox(Jsym[:,1:3*Nv], Jfull[:,1:3*Nv], atol = atol)

    @test norm(J[:,3*Nv+1:end]) < atol
    @test norm(Jsym[:,3*Nv+1:end]) < atol
end


text:
# linear
function update_grads!(::Linear, δ𝑤::Matrix{T}, p::Matrix{T}, y::Vector{T}, α::T) where {T <: AbstractFloat}
    @inbounds for i in eachindex(y)
        δ𝑤[1,i] = 2 * (p[1,i] - y[i]) * δ𝑤[3,i]
        δ𝑤[2,i] = 2 * δ𝑤[3,i]
    end
end

# logistic - on linear predictor
function update_grads!(::Logistic, δ𝑤::Matrix{T}, p::Matrix{T}, y::Vector{T}, α::T) where {T <: AbstractFloat}
    @inbounds for i in eachindex(y)
        pred = sigmoid(p[1,i])
        δ𝑤[1,i] = (pred - y[i]) * δ𝑤[3,i]
        δ𝑤[2,i] = pred * (1 - pred) * δ𝑤[3,i]
    end
end

# Poisson
function update_grads!(::Poisson, δ𝑤::Matrix{T}, p::Matrix{T}, y::Vector{T}, α::T) where {T <: AbstractFloat}
    @inbounds for i in eachindex(y)
        δ𝑤[1,i] = (exp(p[1, i]) .- y[i]) * δ𝑤[3,i]
        δ𝑤[2,i] = exp(p[1, i]) * δ𝑤[3,i]
    end
end

# L1
function update_grads!(::L1, δ𝑤::Matrix{T}, p::Matrix{T}, y::Vector{T}, α::T) where {T <: AbstractFloat}
    @inbounds for i in eachindex(y)
        δ𝑤[1, i] = (α * max(y[i] - p[1,i], 0) - (1 - α) * max(p[1,i] - y[i], 0)) * δ𝑤[3, i]
    end
end

# Softmax
function update_grads!(::Softmax, δ𝑤::Matrix{T}, p::Matrix{T}, y::Vector{S}, α::T) where {T <: AbstractFloat,S}
    p .= p .- maximum(p, dims=1)
    sums = sum(exp.(p), dims=1)
    K = (size(δ𝑤, 1) - 1) ÷ 2
    for i in eachindex(y)
        for k in 1:K
            # δ𝑤[k, i] = (exp(p[k, i]) / sums[i] - (onehot(y[i], 1:K))) * δ𝑤[2 * K + 1, i]
            if k == y[i]
                δ𝑤[k, i] = (exp(p[k, i]) / sums[i] - 1) * δ𝑤[2 * K + 1, i]
            else
                δ𝑤[k, i] = (exp(p[k, i]) / sums[i]) * δ𝑤[2 * K + 1, i]
            end
            δ𝑤[k + K, i] = 1 / sums[i] * (1 - exp(p[k, i]) / sums[i]) * δ𝑤[2 * K + 1, i]
        end
    end
end

# Quantile
function update_grads!(::Quantile, δ𝑤::Matrix{T}, p::Matrix{T}, y::Vector{T}, α::T) where {T <: AbstractFloat}
    @inbounds for i in eachindex(y)
        δ𝑤[1,i] = y[i] > p[1,i] ? α * δ𝑤[3,i] : (α - 1) * δ𝑤[3,i]
        δ𝑤[2,i] = y[i] - p[1,i] # δ² serves to calculate the quantile value - hence no weighting on δ²
    end
end

# Gaussian - http://jrmeyer.github.io/machinelearning/2017/08/18/mle.html
# pred[i][1] = μ
# pred[i][2] = log(σ)
function update_grads!(::Gaussian, δ𝑤::Matrix{T}, p::Matrix{T}, y::Vector{T}, α::T) where {T <: AbstractFloat}
    @inbounds @simd for i in eachindex(y)
        # first order
        δ𝑤[1,i] = (p[1,i] - y[i]) / exp(2 * p[2,i]) * δ𝑤[5,i]
        δ𝑤[2,i] = (1 - (p[1,i] - y[i])^2 / exp(2 * p[2,i])) * δ𝑤[5,i]
        # second order
        δ𝑤[3,i] = δ𝑤[5,i] / exp(2 * p[2,i])
        δ𝑤[4,i] = 2 * δ𝑤[5,i] / exp(2 * p[2,i]) * (p[1,i] - y[i])^2
    end
end

# utility functions
function logit(x::AbstractArray{T,1}) where T <: AbstractFloat
    @. x = log(x / (1 - x))
    return x
end

function logit(x::T) where T <: AbstractFloat
    log(x / (1 - x))
end

function sigmoid(x::AbstractArray{T,1}) where T <: AbstractFloat
    @. x = 1 / (1 + exp(-x))
    return x
end

function sigmoid(x::T) where T <: AbstractFloat
    1 / (1 + exp(-x))
end

function softmax(x::AbstractVector{T}) where T <: AbstractFloat
    x .-= maximum(x)
    x = exp.(x) ./ sum(exp.(x))
    return x
end

##############################
# get the gain metric
##############################
# GradientRegression
function get_gain(::S, ∑::Vector{T}, λ::T, K) where {S <: GradientRegression,T <: AbstractFloat}
    ∑[1]^2 / (∑[2] + λ * ∑[3]) / 2
end

# GaussianRegression
function get_gain(::S, ∑::Vector{T}, λ::T, K) where {S <: GaussianRegression,T <: AbstractFloat}
    (∑[1]^2 / (∑[3] + λ * ∑[5]) + ∑[2]^2 / (∑[4] + λ * ∑[5])) / 2
end

# MultiClassRegression
function get_gain(::S, ∑::Vector{T}, λ::T, K) where {S <: MultiClassRegression,T <: AbstractFloat}
    gain = zero(T)
    @inbounds for k in 1:K
        gain += ∑[k]^2 / (∑[k + K] + λ * ∑[2 * K + 1]) / 2
    end
    return gain
end

# QuantileRegression
function get_gain(::S, ∑::Vector{T}, λ::T, K) where {S <: QuantileRegression,T <: AbstractFloat}
    abs(∑[1])
end

# L1 Regression
function get_gain(::S, ∑::Vector{T}, λ::T, K) where {S <: L1Regression,T <: AbstractFloat}
    abs(∑[1])
end

function update_childs_∑!(::L, nodes, n, bin, feat, K) where {L <: Union{GradientRegression,QuantileRegression,L1Regression}}
    nodes[n << 1].∑ .= nodes[n].hL[feat][(3 * bin - 2):(3 * bin)]
    nodes[n << 1 + 1].∑ .= nodes[n].hR[feat][(3 * bin - 2):(3 * bin)]
    return nothing
end

function update_childs_∑!(::L, nodes, n, bin, feat, K) where {L <: GaussianRegression}
    nodes[n << 1].∑ .= nodes[n].hL[feat][(5 * bin - 4):(5 * bin)]
    nodes[n << 1 + 1].∑ .= nodes[n].hR[feat][(5 * bin - 4):(5 * bin)]
    return nothing
end

function update_childs_∑!(::L, nodes, n, bin, feat, K) where {L <: MultiClassRegression}
    KK = 2 * K + 1
    nodes[n << 1].∑ .= nodes[n].hL[feat][(KK * (bin - 1) + 1):(KK * bin)]
    nodes[n << 1 + 1].∑ .= nodes[n].hR[feat][(KK * (bin - 1) + 1):(KK * bin)]
    return nothing
end


text:
# Stand-alone script to manually check compat when not using GitHub workflow
using CompatHelper
dep_to_current_compat_entry,
    dep_to_current_compat_entry_verbatim,
    dep_to_latest_version,
    deps_with_missing_compat_entry = CompatHelper.get_project_deps("Project.toml")

CompatHelper.get_latest_version_from_registries!(dep_to_latest_version,CompatHelper.default_registries)


text:
using DeformableMirrors
using Test

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


text:
abstract type LinearSolver end

"""
    Empty solver
"""
mutable struct EmptySolver <: LinearSolver
    F::Any
end

function empty_solver(A::Any)
    EmptySolver(A)
end

"""
    LU solver
"""
mutable struct LUSolver{T} <: LinearSolver
    A::Array{T,2}
    ipiv::Vector{Int}
    lda::Int
    info::Base.RefValue{Int}
end

function lu_solver(A)
    m, n = size(A)
    ipiv = similar(A, LinearAlgebra.BlasInt, min(m, n))
    lda  = max(1, stride(A, 2))
    info = Ref{LinearAlgebra.BlasInt}()
    LUSolver(copy(A), ipiv, lda, info)
end

function getrf!(A, ipiv, lda, info)
    Base.require_one_based_indexing(A)
    LinearAlgebra.chkstride1(A)
    m, n = size(A)
    lda  = max(1,stride(A, 2))
    ccall((LinearAlgebra.BLAS.@blasfunc(dgetrf_), Base.liblapack_name), Cvoid,
          (Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{Float64},
           Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}),
          m, n, A, lda, ipiv, info)
    return nothing
end

function factorize!(s::LUSolver{T}, A::AbstractMatrix{T}) where T
    fill!(s.A, 0.0)
    fill!(s.ipiv, 0)
    s.lda = 0
    s.A .= A
    getrf!(s.A, s.ipiv, s.lda, s.info)
end

function linear_solve!(s::LUSolver{T}, x::Vector{T}, A::Matrix{T},
        b::Vector{T}; reg::T = 0.0, fact::Bool = true) where T
    fact && factorize!(s, A)
    x .= b
    LinearAlgebra.LAPACK.getrs!('N', s.A, s.ipiv, x)
end

function linear_solve!(s::LUSolver{T}, x::Matrix{T}, A::Matrix{T},
    b::Matrix{T}; reg::T = 0.0, fact::Bool = true) where T
    fill!(x, 0.0)
    n, m = size(x)
    r_idx = 1:n
    fact && factorize!(s, A)
    x .= b
    for j = 1:m
        xv = @views x[r_idx, j]
        LinearAlgebra.LAPACK.getrs!('N', s.A, s.ipiv, xv)
    end
end


text:
using PkgTemplates
using Base.Test

# write your own tests here
@test 1 == 1


text:
export DifferentialData, average!

struct DifferentialData <: TimeResolvedData
    @def prop Vector{Float64} real δreal imag δimag
    logs::Dict{Symbol,Int} # [dims, avg. num]

    DifferentialData(n::Int) = new(
        Vector{Float64}(undef, n), Vector{Float64}(undef, n),
        Vector{Float64}(undef, n), Vector{Float64}(undef, n),
        Dict(:dims => n, :count => 0)
    )

    function DifferentialData(s::String, c::Int=1; Ireal::Int=2)
        data = readdlm(s)
        dims = size(data, 1)
        real = Vector{Float64}(undef, dims)
        @simd for i in eachindex(1:dims)
            @inbounds real[i] = data[i, Ireal]
        end
        return new(
            real, zeros(Float64, dims), zeros(Float64, dims), zeros(Float64, dims),
            Dict(:dims => dims, :count => c)
        )
    end
end

# Perform a single step of Wolford algorithm
function welford_step(μ::Real, s::Real, v::Real, c::Real)
    isone(c) && return v, zero(v)
    s = s * (c - 1)
    m = μ + (v - μ) / c
    s = s + (v - μ) * (v - m)
    μ = m
    return μ, s / (c - 1)
end

function average!(des::DifferentialData, src::MatI, Ireal::Int, Iimag::Int, Iref::Int)
    @get des real δreal imag δimag logs

    one2n = eachindex(1:logs[:dims])
    c = logs[:count] + 1

    @inbounds for i in one2n
        ref = src[i, Iref]
        real[i], δreal[i] = welford_step(real[i], δreal[i], src[i, Ireal] / ref, c)
        imag[i], δimag[i] = welford_step(imag[i], δimag[i], src[i, Iimag] / ref, c)
    end

    logs[:count] = c
    return nothing
end


text:
using DataFrames
using Dates
using DataStructures
using Ipopt
using MaximumGeneratorProfitTCRDD
using PowerSystems
using PowerSimulations
using Test
using TimeSeries

#  Package Constants
const PSY = PowerSystems
const PSI = PowerSimulations

@testset "MaximumGeneratorProfitTCRDD.jl" begin
    # Settings for testset
    include("../examples/data/example_case4/buildup_case4.jl")
    include("../examples/data/example_case4/c_sys_case4.jl")
    include("max_gen_prof.jl")
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
using DifferentiableStateSpaceModels, Test, LinearAlgebra
using DifferentiableStateSpaceModels: vech, inv_vech
@testset "vech" begin
    A = [1 2; 3 4]
    @test vech(A) == [1, 3, 4]
    @test vech(LowerTriangular(A)) == [1, 3, 4]
    @test_throws DimensionMismatch vech([1 2 4; 3 4 5])

    @test inv_vech([1, 3, 4], 2) == LowerTriangular(A)
    @test inv_vech([1, 3, 4]) == LowerTriangular(A)
    @test_throws ErrorException inv_vech([1, 3, 4, 5])
    @test_throws ErrorException inv_vech([1, 3, 4], 3)

    B = LowerTriangular(rand(10, 10))
    @test inv_vech(vech(B)) ≈ B
end