text:
######################################################################
# old_cache.jl
#
# like infotheoretic.jl, but allows caching
#
# how this used to be implemented
######################################################################

export make_cache
function make_cache(df::DF)
    cache = zeros(df.n, df.n, df.n, df.n, length(df.obs_list))
    for xi = 1:df.n, yi = 1:df.n
        x = (xi-0.5) * df.cell_size
        y = (yi-0.5) * df.cell_size
        xp = (x,y,0.0)
        for txi = 1:df.n, tyi = 1:df.n
            tx = (txi-0.5) * df.cell_size
            ty = (tyi-0.5) * df.cell_size
            theta = (tx,ty,0.0,0.0)
            for (oi,o) in enumerate(df.obs_list)
                cache[xi,yi,txi,tyi,oi] = O(df.sensor, theta, xp, o)
            end
        end
    end
    return cache
end

function mutual_information(df::DF, cache::Array{Float64,5})
    mut_info = zeros(df.n, df.n)
    for xi = 1:df.n, yi = 1:df.n
        # TODO: really, should have some other heading
        xp = (xi, yi, 0)
        mut_info[xi,yi] = mutual_information(df, xp, cache)
    end
    return mut_info
end

function mutual_information(df::DF, xp::NTuple{3,Int}, cache::Array{Float64,5})
	H_o = 0.0
	H_o_t = 0.0

    for (oi,o) in enumerate(df.obs_list)
		po = p_obs(df, xp, oi, cache)
		if po > 0.0
			H_o -= po * log(po)
		end

		# sum over possible jammer locations
		for txi = 1:df.n, tyi = 1:df.n
            if df.b[txi, tyi] > 0.0
                pot = cache[ xp[1], xp[2], txi, tyi, oi ]
                if pot > 0.0
                    H_o_t -= pot * df.b[txi, tyi] * log(pot)
                end
            end
		end
	end
	return H_o - H_o_t
end

function p_obs(df::DF, xp::NTuple{3,Int}, oi::Int, cache::Array{Float64,5})
    prob = 0.0
    for txi = 1:df.n, tyi = 1:df.n
        if df.b[txi,tyi] > 0.0
            prob += df.b[txi, tyi] * cache[ xp[1], xp[2], txi, tyi, oi ]
        end
    end
    return prob
end


text:
import EnergyExpressions: @above_diagonal_loop, cofactor, nonzero_minors, find_diagonal_blocks, NBodyMatrixElement
using SparseArrays
using Combinatorics

# Used to generate reference data to test the linear complexity
# algorithm against; the naïve algorithm is correct, but is of
# factorial complexity and cannot be used for more than ~20 orbitals.
@generated function naive_nonzero_minors(op::NBodyOperator{N}, overlap, numorbitals) where N
    quote
        ks = Vector{Vector{Int}}()
        ls = Vector{Vector{Int}}()

        @above_diagonal_loop $N k numorbitals begin
            kN = zeros(Int, $N)
            Base.Cartesian.@nexprs $N d -> kN[d] = k_d

            @above_diagonal_loop $N l numorbitals begin
                lN = zeros(Int, $N)
                Base.Cartesian.@nexprs $N d -> lN[d] = l_d

                Dkl = cofactor(kN, lN, overlap) # Determinantal overlap of all other orbitals
                iszero(Dkl) && continue
                push!(ks, kN)
                push!(ls, lN)
            end
        end
        ks, ls
    end
end

struct ConcreteNBodyOperator{N} <: NBodyOperator{N} end

function test_find_minors(cfga, cfgb, overlaps=OrbitalOverlap[])
    a = SlaterDeterminant(cfga)
    b = SlaterDeterminant(cfgb)

    S = overlap_matrix(a, b, overlaps)

    h = ConcreteNBodyOperator{1}()
    g = ConcreteNBodyOperator{2}()

    errors = 0
    for (op,label) in [(h,"γ"), (g,"Γ")]
        ks,ls = naive_nonzero_minors(op, S, length(a))
        N = numbodies(op)
        ks′,ls′ = nonzero_minors(N, S)
        kls′ = collect(zip(sort.(ks′),sort.(ls′)))
        for (k,l) in zip(sort.(ks),sort.(ls))
            if (k,l) ∉ kls′
                errors += 1
                kl = nice_string(k)
                ll = nice_string(l)
                println("$label($kl|$ll) missing from new algorithm")
            end
        end
    end
    @test errors == 0
end

function find_diagonal_blocks_test_case(A, blocks)
    @test find_diagonal_blocks(A) == blocks
    @test find_diagonal_blocks(SparseMatrixCSC(A')) == blocks
end

function test_cofactors(N, ii, jj)
    orbitals = 1:N
    excite = i -> i+100
    virtuals = map(excite, orbitals)

    configurations = [SlaterDeterminant(vcat(1:i-1,excite(i),i+1:N))
                      for i=1:N]
    overlaps = map(((a,b),) -> OrbitalOverlap(a,b),
                   reduce(vcat, [collect(combinations(virtuals,2)),
                                 [[a,a] for a in virtuals]]))

    results = Vector{EnergyExpressions.NBodyMatrixElement}()

    for (i,j) = [(ii,jj),(jj,ii)]
        a = configurations[i]
        b = configurations[j]
        S = overlap_matrix(a, b, overlaps)
        ks,ls = nonzero_minors(1, S)
        for (k,l) in zip(ks,ls)
            push!(results, cofactor(k, l, S))
        end
    end
    results
end

@testset "Minor determinants" begin
    @testset "Find diagonal blocks" begin
        find_diagonal_blocks_test_case((sparse([3,1,2], 1:3, ones(Int, 3))), [1:3])
        find_diagonal_blocks_test_case((sparse([1,4,2,3], 1:4, ones(Int, 4))), [1:1, 2:4])
        find_diagonal_blocks_test_case((sparse([1,2,5,3,4], 1:5, ones(Int, 5))), [1:1, 2:2, 3:5])
        find_diagonal_blocks_test_case((sparse([1,2,3,6,4,5], 1:6, ones(Int, 6))), [1:1, 2:2, 3:3, 4:6])
    end

    @testset "Cofactors" begin
        @test test_cofactors(7,4,7) == [NBodyMatrixElement([-OrbitalOverlap(104,107)]),
                                          NBodyMatrixElement([-OrbitalOverlap(107,104)])]
        @test test_cofactors(8,5,8) == [NBodyMatrixElement([-OrbitalOverlap(105,108)]),
                                          NBodyMatrixElement([-OrbitalOverlap(108,105)])]
        @test test_cofactors(9,6,9) == [NBodyMatrixElement([-OrbitalOverlap(106,109)]),
                                          NBodyMatrixElement([-OrbitalOverlap(109,106)])]
        @test test_cofactors(9,7,9) == [NBodyMatrixElement([-OrbitalOverlap(107,109)]),
                                          NBodyMatrixElement([-OrbitalOverlap(109,107)])]
        @test test_cofactors(54,45,54) == [NBodyMatrixElement([-OrbitalOverlap(145,154)]),
                                             NBodyMatrixElement([-OrbitalOverlap(154,145)])]
    end

    @testset "Base cases" begin
        test_find_minors([:a], [:a])
        test_find_minors([:a,:b], [:a,:b])
        test_find_minors([:a,:b], [:b,:a])
        test_find_minors([:a,:b], [:a,:c])
        test_find_minors([:a,:b], [:c,:a])
        test_find_minors([:a,:b,:c], [:a,:b,:c])
        test_find_minors([:a,:b,:c], [:a,:b,:d])
        test_find_minors([:a,:b,:c], [:a,:b,:d], [OrbitalOverlap(:a,:b)])
        test_find_minors([:a,:b,:c], [:a,:b,:d], [OrbitalOverlap(:a,:c)])
        test_find_minors([:a,:b,:c], [:a,:b,:d], [OrbitalOverlap(:b,:d),OrbitalOverlap(:c,:d)])
        test_find_minors([:a,:b,:c], [:a,:b,:d], [OrbitalOverlap(:a,:c),OrbitalOverlap(:c,:d)])
        test_find_minors([:a,:b,:c], [:a,:b,:d], [OrbitalOverlap(:a,:d),OrbitalOverlap(:c,:d)])
        test_find_minors([:a,:b,:c], [:b,:c,:a])
        test_find_minors([:a,:b,:c], [:c,:b,:a])
        test_find_minors([:a,:b,:c], [:b,:c,:a], [OrbitalOverlap(:a,:c)])
        test_find_minors([:a,:b,:c], [:a,:d,:e])
        test_find_minors([:a,:b,:c,:d], [:a,:b,:c,:d])
        test_find_minors([:a,:b,:c,:d], [:a,:b,:c,:d], [OrbitalOverlap(:c,:d)])
        test_find_minors([:a,:b,:c,:d], [:a,:b,:c,:e], [OrbitalOverlap(:c,:d)])
        test_find_minors([:a,:b,:c,:d], [:a,:b,:c,:e], [OrbitalOverlap(:c,:e)])
        test_find_minors([:a,:b,:c,:d], [:a,:b,:e,:f], [OrbitalOverlap(:c,:b),OrbitalOverlap(:d,:b)])
        test_find_minors([:a,:b,:c,:d], [:a,:b,:e,:f], [OrbitalOverlap(:a,:e),OrbitalOverlap(:c,:f),OrbitalOverlap(:d,:b)])
        test_find_minors([:b,:c,:d,:k_a], [:a,:c,:d,:k_b], [OrbitalOverlap(:k_a,:k_b)])
        test_find_minors([:k_a,:b,:c,:d], [:a,:k_b,:c,:d], [OrbitalOverlap(:k_a,:k_b)])
        test_find_minors([:a, :b, :c, :d], [:a, :e, :f, :g])
        test_find_minors([:a,:b,:c,:d], [:d,:c,:b,:a])
    end

    @testset "Atomic cases" begin
        h = ConcreteNBodyOperator{1}()
        g = ConcreteNBodyOperator{2}()

        for (cfga,cfgb,exp1,exp2) in [
            (c"1s2",c"1s2",2,1),
            (c"1s2",c"1s 2p",1,1),
            (c"1s2 2s2",c"1s2 2s2",4,binomial(4,2)),
            (c"1s2 2s2 2p6",c"1s2 2s2 2p6",10,binomial(10,2)),
            (c"1s2 2s2 2p6 3s2",c"1s2 2s2 2p6 3s2",12,binomial(12,2)),
            (c"1s2 2s2 2p6 3s2 3p6",c"1s2 2s2 2p6 3s2 3p6",18,binomial(18,2)),
            (c"[Ar]",c"[Ar]",18,binomial(18,2)),
            (c"[Kr]",c"[Kr]",36,binomial(36,2)),
            (c"[Xe]",c"[Xe]",54,binomial(54,2))
        ]
            a = SlaterDeterminant(spin_configurations(cfga)[1])
            b = SlaterDeterminant(spin_configurations(cfgb)[1])

            S = overlap_matrix(a, b)

            for (op,ex) in [(h,exp1),(g,exp2)]
                ks,ls = nonzero_minors(numbodies(op), S)
                @test length(ks) == length(ls) == ex
            end
        end
    end
end

ext: jl


text:
#  auto_pull_matlabdatasets.jl
#
#   Puts together the tools needed to parse a collation folder for its datasets
#
#. collatedPath = the path to the collated CSV folder for all the sessions to run
#.     This folder has a folder for each session, within
#			within this folder is singletrial, baseline, LOI so that we can import data for each file
#
using GLM

function run_collated_model(collatedPath::String, modelpackagefunction::Function;
	runID=0, pathIDx = [], runFails=false, failDirs=[],
	postprocessingfunction::Function=template_postprocessingfunction,
	compositesavepath::String="", suppressFigures=false)
	start_dir = pwd()
	if isempty(compositesavepath)
	#
	# Make a results folder ABOVE our collated path
	#
		cd(collatedPath)
		cd("..")
		try
			mkdir(join(["Collated_results_", modelpackagefunction(""; runID=runID, getpackagename=true)]))
		catch
		end
		cd(join(["Collated_results_", modelpackagefunction(""; runID=runID, getpackagename=true)]))
		compositesavepath = pwd()
		cd(start_dir)
	end
	# Can rerun an ID with runID !=0
	if runID==0
		runID = rand(1:10000)
	end
	# Find all the sessions in the folder
	# pathIDx allows you to select only some of the sessions in the path. only use this for forward mode (NOT FAIL MODE)
	println("-----------------------------------------------")
    println("	")
    println("Initializing run_collated_model for ", modelpackagefunction(""; runID=runID, getpackagename=true))
    println("	")

    try
    	root = collatedPath
    	dirs = String[]
    	files = String[]
    	a = readdir(collatedPath)
    	for i in a
		    if i == ".DS_Store"
		        #skip
		    elseif isdir(joinpath(root, i))
		        push!(dirs, i)
		    elseif isfile(joinpath(root, i))
		        push!(files, i)
		    end
		end
    	# for (rooti, dirsi, filesi) in walkdir(collatedPath)
    	# 	root = rooti
    	# 	push!(dirs, dirsi)
    	# 	push!(files, filesi)
    	# end
        if !runFails
	        println("Found ", length(dirs), " sessions in $root:")
	        pretty_print_list(dirs, orient="vertical", enum=true)
	        if isempty(pathIDx)
	        	println("	Using all these sessions.")
	        	dirs = dirs
        	else
				println("	Using only a subset of these sessions:")
				dirs = dirs[pathIDx]
				pretty_print_list(dirs, orient="vertical", enum=false)
        	end
        else
        	println("Re-running failures in $root:")
        	pretty_print_list(failDirs, orient="vertical", enum=true)
        	dirs = failDirs
        end
        sessionIDs = String[]
        sessionPaths = String[]
        sessionIdx = Vector{Int}()
        results = Vector{DataFrame}()
        failDirs = String[]
        result_df = DataFrame([Int,String, DataFrame, String], [:sessionIdx, :sessionIDs, :results, :sessionPaths])
        println("	")

	    for i = 1:length(dirs)
	    	dir = dirs[i]
	        # println(joinpath(root, dir)) # path to directories
	        try
	        	println("-----------------------------------------------")
	        	println("Processing ", dir, "...(", timestamp_now(), ")")
	        	print("   ")
	        	progressbar(i, length(dirs))
	        	println("	")
	        	push!(sessionIdx, i)
		        push!(sessionIDs, dir)
		        push!(sessionPaths, joinpath(root, dir))

		        result = modelpackagefunction(joinpath(root, dir), sessionID = dir, runID = runID, suppressFigures=suppressFigures)
		        push!(results, result)

	        	result_df = DataFrame(sessionIdx=sessionIdx, sessionIDs=sessionIDs, results = results, sessionPaths=sessionPaths)
	        	if suppressFigures
	        		close() # closes all the open figures
        		end
	        catch
	        	println("	!********* Encountered error! Skipping this directory")
	        	println("	")
	        	push!(failDirs, dir)
	        	# rethrow()
	        end
	    end
	    println("-----------------------------------------------")
	    println("	")
	    println("Completed procesing of ", length(dirs) - length(failDirs), " sessions. (", timestamp_now(), ") ~")
	    println("	")
	    println("	Find results in each sessions' folder in: ", modelpackagefunction(""; runID=runID, getpackagename=true))
	    println(collatedPath)
	    println("	")
	    println("Initiating post-modeling collation of results...")
	    postprocessingfunction(result_df, compositesavepath, modelpackagefunction; runID=runID)
	    println("Post-modeling collation of results complete and variables saved to:")
	    println(compositesavepath)
	    println("	")
	    println("-----------------------------------------------")
	    println("	")
	    cd(start_dir)
	    return failDirs, result_df
    catch
        cd(start_dir)
        try
        	typeof(result_df)
    	catch
    		result_df = []
		end
        # return failDirs, result_df
        rethrow()
    end
end

function template_modelpackage(path; sessionID ="", getpackagename=false, runID=0, suppressFigures=false)
	#
	# Use this to build a new analysis
	#
	# name the package and runID
	packagename = join(["templatepkg_",runID])
	if getpackagename
		return packagename
	end

	# do the business of the package on this session
	println("	good!")

	# Try to enter the results folder
	savepath = joinpath(path, join(["results_", packagename]))
	try
		cd(savepath)
	catch
		mkdir(savepath)
		cd(savepath)
	end

	# Save each variable to our results folder
	CSV.write("result_good.csv", DataFrame(good = ["good", "good"]))
	CSV.write("result_othervar.csv", DataFrame(othervar = ["var1", "var2"]))

	# make a working result df with all the results to keep in workspace
	result = DataFrame(
		good = ["good", "good"],
		othervar = ["var1", "var2"]
		)
	return result
end
function template_postprocessingfunction(results::DataFrame, compositesavepath, modelpackagefunction; runID=0)
	#
	# Use this to compile the analysis
	#
	# name the package and runID
	packagename = modelpackagefunction(""; sessionID ="", getpackagename=true, runID=runID)

	# do the business of the package on this session
	println("	This is a template postprocessingfunction. To do work, you need to
		implement what you want to collate from the results here!") #"

	# Try to enter the composite results folder
	cd(compositesavepath)

	# Save each variable to our results folder
	CSV.write(join([packagename, "_result_good_",".csv"]), DataFrame(good = ["good", "good"]))
	CSV.write(join([packagename, "_result_othervar_",".csv"]), DataFrame(othervar = ["var1", "var2"]))
end

function bootlogit_modelpackage1(path; sessionID ="", getpackagename=false, runID=0, suppressFigures=false)
# name the package and runID
	packagename = join(["bootlogit_modelpackage1_",runID])
	if getpackagename
		return packagename
	end
# Try to enter the results folder
	savepath = joinpath(path, join(["results_", packagename]))
	figurePath = joinpath(path, join(["figures_", packagename]))
	try
		cd(savepath)
	catch
		mkdir(savepath)
		mkdir(figurePath)
		cd(savepath)
	end
# do the business of the package on this session
	#
	# first, we extract the relevant data: the singletrial, baseline and LOI sets and make a df
	# (expecting singletrial, baseline, and LOI folders for each dataset with CSV files from matlab)
	#
	ndf = extract_data_with_baselineandLOI(path; normalize=true)
	#
	# next, we want to specify our model formulae, including the nested models
	#
	formulas = [
	    @formula(LickState ~ Y),
	    @formula(LickState ~ Mean_Baseline),
	    @formula(LickState ~ Median_Baseline),
	    @formula(LickState ~ Mean_LOI),
	    @formula(LickState ~ Median_LOI),
	    @formula(LickState ~ Mean_Baseline + Mean_LOI + Y),
	    @formula(LickState ~ Median_Baseline + Median_LOI + Y),
	    @formula(LickState ~ LickTime_1back),
	    @formula(LickState ~ LickTime_2back),
	    @formula(LickState ~ LickTime_2back + LickTime_1back),
	    @formula(LickState ~ Rxn_2back + Early_2back + Reward_2back + ITI_2back),
	    @formula(LickState ~ Rxn_1back + Early_1back + Reward_1back + ITI_1back),
	    @formula(LickState ~ Rxn_2back + Early_2back + Reward_2back + ITI_2back +
	    	Rxn_1back + Early_1back + Reward_1back + ITI_1back),

	    @formula(LickState ~ Rxn_2back + Early_2back + Reward_2back + ITI_2back +
	    	Rxn_1back + Early_1back + Reward_1back + ITI_1back +
	    	Y),
	    @formula(LickState ~ Rxn_2back + Early_2back + Reward_2back + ITI_2back +
	    	Rxn_1back + Early_1back + Reward_1back + ITI_1back +
	        Mean_Baseline + Mean_LOI + Y),
	    @formula(LickState ~ Rxn_2back + Early_2back + Reward_2back + ITI_2back +
	    	Rxn_1back + Early_1back + Reward_1back + ITI_1back +
	        Median_Baseline + Median_LOI + Y),

	    @formula(LickState ~ LickTime_2back + LickTime_1back +
	    	Rxn_2back + Early_2back + Reward_2back + ITI_2back +
	    	Rxn_1back + Early_1back + Reward_1back + ITI_1back +
	        Y),
	    @formula(LickState ~ LickTime_2back + LickTime_1back +
	    	Rxn_2back + Early_2back + Reward_2back + ITI_2back +
	    	Rxn_1back + Early_1back + Reward_1back + ITI_1back +
	        Mean_Baseline + Mean_LOI + Y),
	    @formula(LickState ~ LickTime_2back + LickTime_1back +
	    	Rxn_2back + Early_2back + Reward_2back + ITI_2back +
	    	Rxn_1back + Early_1back + Reward_1back + ITI_1back +
	        Median_Baseline + Median_LOI + Y),
		]

	modelNames = [
	    "DA-only",
	    "μBl-only",
	    "medBl-only",
	    "μLOI-only",
	    "medLOI-only",
	    "DA-μBl-μLOI",
	    "DA-medBl-medLOI",
	    "Lt1b-only",
	    "Lt2b-only",
	    "Lt1b-Lt2b",
	    "oc1b-only",
	    "oc1b-oc2b",
	    "DA-oc1b-oc2b",
	    "DA-μBl-μLOI-oc1b-oc2b",
	    "DA-medBl-medLOI-oc1b-oc2b",
	    "DA-Lt1b-Lt2b-oc1b-oc2b",
	    "DA-μBl-μLOI-Lt1b-Lt2b-oc1b-oc2b",
	    "DA-medBl-medLOI-Lt1b-Lt2b-oc1b-oc2b",
	]

	results = modelSelectionByAICBICxval(ndf, :LickState, formulas, modelNames, "logit";
    		n_iters=100,updownsampleYID=true, figurePath=figurePath, savePath = savepath, suppressFigures=suppressFigures)
# Save each variable to our results folder
	# this is already handled by the modelSelectionByAICBICxval function

# make a working result df with all the results to keep in workspace
	result = results
	return result#, ndf
end

function extract_data_with_baselineandLOI(path; normalize=true)
	data_single_trial = extract_data(joinpath(path, "singletrial"), blmode=false, LOImode=false)
	# data_single_trial.sessionCode = [sessionID for _=1:length(data_single_trial.sessionCode)]
	data_baseline_trial = extract_data(joinpath(path, "baseline"), blmode=true, LOImode=false)
	data_LOI_trial = extract_data(joinpath(path, "LOI"), blmode=false, LOImode=true)
	df = makeSessionDataFrame(data_single_trial; normalize=normalize, includeBL_LOI=true, baseline_data=data_baseline_trial, LOI_data=data_LOI_trial)
	return df
end

function bootlogit_postprocessingfunction1(results::DataFrame, compositesavepath, modelpackagefunction; runID=0)
	#
	# Use this to compile the analysis
	#
	# name the package and runID
	packagename = modelpackagefunction(""; sessionID ="", getpackagename=true, runID=runID)

	# do the business of the package on this session
	combine_th_across_sessions(results, compositesavepath, runID, packagename)
	combine_AICBIC_across_sessions(results, compositesavepath, runID, packagename)

	# Try to enter the composite results folder
	# cd(compositesavepath)

	# Save each variable to our results folder
	# 		This is handled by the combine functions ABOVE
	# CSV.write(join([packagename, "_result_good_",".csv"]), DataFrame(good = ["good", "good"]))
	# CSV.write(join([packagename, "_result_othervar_",".csv"]), DataFrame(othervar = ["var1", "var2"]))
end

#
# 	Collation
#
function combine_th_across_sessions(results, compositesavepath, runID,packagename)
    ret_dir = pwd()
    cd(compositesavepath)
    modelNames = unique(results.results[1].th_summary[1].modelName)
    n_sesh = length(results.results)
    i_sesh=1
    i_model=1
    i_iter=1
    n_models = length(results.results[i_sesh].AICs)
    n_iters = length(results.results[i_sesh].AICs[i_model])
    axs = []
    fs = []
    th_summary = DataFrame(modelName=[], composite_th=[], composite_se=[], composite_CImin=[], composite_CImax=[])
    for model = 1:n_models
        result_df = DataFrame(train_dof = [results.results[1].dofs[model][1]], train_ths = [[results.results[1].ths[model][1]]], train_se_ths=[[results.results[1].se_ths[model][1]]])
        for i_sesh = 1:n_sesh
            if i_sesh==1
                iterstart=2
            else
                iterstart=1
            end
            for i=iterstart:n_iters
                append!(result_df,
                DataFrame(train_dof = [results.results[i_sesh].dofs[model][i]],
                        train_ths = [[results.results[i_sesh].ths[model][i]]],
                        train_se_ths=[[results.results[i_sesh].se_ths[model][i]]],
                    ))
            end
        end

        (composite_th, composite_se, composite_CImin, composite_CImax, ax, f) = theta_summary(result_df; Mode = "sparseFit", result_df=result_df)
        title(modelNames[model])
        push!(axs, ax)
        push!(fs, f)

        append!(th_summary,
            DataFrame(
                modelName = modelNames[model],
                composite_th=composite_th,
                composite_se=composite_se,
                composite_CImin=composite_CImin,
                composite_CImax=composite_CImax,
                ))

        #
        # Save the composite variables
        #
        CSV.write(join(["MODELno",model, "_composite_ths_nboot", n_iters, "_nsesh", n_sesh,
                    "_", packagename, ".csv"]),DataFrame(train_ths = result_df.train_ths))
        CSV.write(join(["MODELno",model, "_composite_se_ths_nboot", n_iters, "_nsesh", n_sesh,
                    "_", packagename, ".csv"]),DataFrame(train_se_ths = result_df.train_se_ths))
        CSV.write(join(["MODELno",model, "_composite_dofs_nboot", n_iters, "_nsesh", n_sesh,
                    "_", packagename, ".csv"]),DataFrame(train_dof = result_df.train_dof))
        CSV.write(join(["MODELno",model, "_composite_th_summary_nboot", n_iters, "_nsesh", n_sesh,
                    "_", packagename, ".csv"]),DataFrame(th_summary = th_summary))
    end

    set_xaxes_same_scale(axs)
    set_yaxes_same_scale(axs)
    for i=1:length(fs)
        # println(i)
        printFigure(join(["composite_", modelNames[i], "_theta_summary_nboot", n_iters, "_nsesh", n_sesh, "_", packagename]); fig=fs[i],figurePath=compositesavepath)
    end
    cd(ret_dir)
end

function combine_AICBIC_across_sessions(results, compositesavepath, runID, packagename)
    n_sesh = length(results.results)
    i_sesh=1
    i_model=1
    i_iter=1
    n_models = length(results.results[i_sesh].AICs)
    n_iters = length(results.results[i_sesh].AICs[i_model])
    allAICs = [[] for _=1:n_models]
    allAICcs = [[] for _=1:n_models]
    allBICs = [[] for _=1:n_models]
    all_Sn_accuracy = [[] for _=1:n_models]
    all_test_accuracy = [[] for _=1:n_models]

    for i_sesh = 1:n_sesh
        for i_model = 1:n_models
            AICs = [results.results[i_sesh].AICs[i_model][i_iter] for i_iter = 1:n_iters]
            append!(allAICs[i_model], AICs)
            AICcs = [results.results[i_sesh].AICcs[i_model][i_iter] for i_iter = 1:n_iters]
            append!(allAICcs[i_model], AICcs)
            BICs = [results.results[i_sesh].BICs[i_model][i_iter] for i_iter = 1:n_iters]
            append!(allBICs[i_model], BICs)
            Sn_accuracies = [results.results[i_sesh].Sn_accuracy[i_model][i_iter] for i_iter = 1:n_iters]
            append!(all_Sn_accuracy[i_model], Sn_accuracies)
            test_accuracies = [results.results[i_sesh].test_accuracy[i_model][i_iter] for i_iter = 1:n_iters]
            append!(all_test_accuracy[i_model], test_accuracies)
        end
    end
    mean_all_AICs = [mean(allAICs[i]) for i=1:n_models]
    mean_all_AICcs = [mean(allAICcs[i]) for i=1:n_models]
    mean_all_BICs = [mean(allBICs[i]) for i=1:n_models]
    mean_all_Sn_accuracy = [mean(all_Sn_accuracy[i]) for i=1:n_models]
    mean_all_test_accuracy = [mean(all_test_accuracy[i]) for i=1:n_models]
    f = figure(figsize=(20,3))
    ax1=subplot(1,3,1)
    compare_AICBIC(mean_all_AICs, allAICs; yl=join(["AIC nsesh=",n_sesh]), iters=n_iters, ax=ax1, minmax="min")
    ax2=subplot(1,3,2)
    compare_AICBIC(mean_all_AICcs, allAICcs; yl=join(["AICc nsesh=",n_sesh]), iters=n_iters, ax=ax2, minmax="min")
    ax3=subplot(1,3,3)
    compare_AICBIC(mean_all_BICs, allBICs; yl=join(["BIC nsesh=",n_sesh]), iters=n_iters, ax=ax3, minmax="min")
    printFigure(join(["compositeAICBIC_summary_nboot", n_iters, "_nsesh", n_sesh, "_", packagename]); fig=f, figurePath=compositesavepath)

    f = figure(figsize=(20,3))
    ax1=subplot(1,3,1)
    compare_AICBIC(mean_all_Sn_accuracy, all_Sn_accuracy; yl=join(["Train Accuracy nsesh=",n_sesh]), iters=n_iters, ax=ax1, minmax="max")
    ax1.set_ylim([0., 1.])
    ax2=subplot(1,3,2)
    compare_AICBIC(mean_all_test_accuracy, all_test_accuracy; yl=join(["Test Accuracy nsesh=",n_sesh]), iters=n_iters, ax=ax2, minmax="max")
    ax2.set_ylim([0., 1.])
    ax3=subplot(1,3,3)
    printFigure(join(["composite_Accuracy_summary_nboot", n_iters, "_nsesh", n_sesh, "_", packagename]); fig=f, figurePath=compositesavepath)
    #
    # Save the variables to the composite folder
    #
    ret_dir = pwd()
    cd(compositesavepath)

    CSV.write(join(["composite_AICs_nboot", n_iters, "_nsesh", n_sesh,
                "_", packagename, ".csv"]),DataFrame(allAICs = allAICs))
    CSV.write(join(["composite_meanAIC_nboot", n_iters, "_nsesh", n_sesh,
                "_", packagename, ".csv"]),DataFrame(mean_all_AICs = mean_all_AICs))
    CSV.write(join(["composite_AICcs_nboot", n_iters, "_nsesh", n_sesh,
                "_", packagename, ".csv"]),DataFrame(allAICcs = allAICcs))
    CSV.write(join(["composite_meanAICc_nboot", n_iters, "_nsesh", n_sesh,
                "_", packagename, ".csv"]),DataFrame(mean_all_AICcs = mean_all_AICcs))
    CSV.write(join(["composite_BICs_nboot", n_iters, "_nsesh", n_sesh,
                "_", packagename, ".csv"]),DataFrame(allBICs = allBICs))
    CSV.write(join(["composite_meanBIC_nboot", n_iters, "_nsesh", n_sesh,
                "_", packagename, ".csv"]),DataFrame(mean_all_BICs = mean_all_BICs))

    CSV.write(join(["composite_Sn_accuracy_nboot", n_iters, "_nsesh", n_sesh,
                "_", packagename, ".csv"]),DataFrame(all_Sn_accuracy = all_Sn_accuracy))
    CSV.write(join(["composite_meanSn_accuracy_nboot", n_iters, "_nsesh", n_sesh,
                "_", packagename, ".csv"]),DataFrame(mean_all_Sn_accuracy = mean_all_Sn_accuracy))
    CSV.write(join(["composite_test_accuracy_nboot", n_iters, "_nsesh", n_sesh,
                "_", packagename, ".csv"]),DataFrame(all_test_accuracy = all_test_accuracy))
    CSV.write(join(["composite_meantest_accuracy_nboot", n_iters, "_nsesh", n_sesh,
                "_", packagename, ".csv"]),DataFrame(mean_all_test_accuracy = mean_all_test_accuracy))

    cd(ret_dir)
end


text:
function pixelshufflehead(k_in, k_out; n_upscale=3, k_mid=64, σ=relu)
    layers = []

    push!(layers, Chain(
        PixelShuffle(2, k_in, k_mid),
        BatchNorm(k_mid, σ),
    ))
    for i = 2:n_upscale
        push!(layers, Chain(
            PixelShuffle(2, k_mid),
            BatchNorm(k_mid, σ),

        ))
    end

    push!(layers, Conv((1, 1), k_mid => k_out))

    return Chain(layers...)
end

function visionhead(
        k_in,
        outsz::Tuple;
        ks_dense=[512],
        p=0.,
        concat_pool=true,
        bn_first=true,
        bn_final=true,
        y_range=nothing,
        act=relu)

    k_out = prod(outsz)
    hs = vcat([concat_pool ? 2k_in : k_in], ks_dense, [k_out])
    n = length(hs)
    bns = trues(n)
    bns[1] = bn_first
    acts = vcat([relu for _ ∈ 1:n-2], [identity])
    pool = concat_pool ? AdaptiveConcatPool((1, 1)) : AdaptiveMeanPool((1, 1))

    layers = [pool, flatten]

    for (h_in, h_out, act) in zip(hs, hs[2:end], acts)
        push!(layers, linbndrop(h_in, h_out, act=act, p=p))
    end

    if !isnothing(y_range)
        min, max = y_range
        push!(layers, x -> Flux.σ.(x) .* (max - min) .- max)
    end
    # reshape to multi-dimensional output if wanted
    if length(outsz) > 1
        push!(layers, xs -> reshape(xs, (outsz..., last(size(xs)))))
    end
    return Chain(layers...)
end
visionhead(k_in, k_out::Int; kwargs...) = visionhead(k_in, (k_out,); kwargs...)

function linbndrop(h_in, h_out; use_bn=true, p=0., act=identity, lin_first=false)
    bn = BatchNorm(lin_first ? h_out : h_in)
    dropout = p == 0 ? identity : Dropout(p)
    dense = Dense(h_in, h_out, act; bias=!use_bn)
    if lin_first
        return Chain(dense, bn, dropout)
    else
        return Chain(bn, dropout, dense)
    end
end


text:
# This file was generated, do not modify it. # hide
bs = MLJ.range(nnregressor, :batch_size, lower=1, upper=5)

tm = MLJ.TunedModel(model=nnregressor, ranges=[bs, ], measure=MLJ.l2)


text: println("No to war")


text:
using Test
using OpenTrick

@testset "open failed" begin
    @debug "open failed test"
    @test_throws SystemError opentrick(open, "nosuchfile.txt")
    @debug "before test length"
    @test length(OpenTrick.tasks_pending) == 0
end


text:
#"Read" data for parametric sinusoidal "simulation" & parameter extraction
#-------------------------------------------------------------------------------

using CMDimData
using CMDimData.MDDatasets

const π = MathConstants.π

#==Emulate reading in simulated data file
===============================================================================#

#=COMMENT
The code below emulates a parametric "simulation" of a sinusoidal response where
the 𝜙, A, and 𝑓 parameters of `signal = A * sin(𝜔*t + 𝜙); 𝜔 = 2π*𝑓` are varied.

The parametric signal can therefore be fully represented as:
	signal(𝜙, A, 𝑓, t)
=#

#But really construct multidimensional DataRS dataset manually:
signal = fill(DataRS, PSweep("phi", [0, 0.5, 1] .* (π/4))) do 𝜙
	fill(DataRS, PSweep("A", [1, 2, 4] .* 1e-3)) do A
	#Inner-most sweep: need to specify element type (DataF1):
	#(Other (scalar) element types: DataInt/DataFloat/DataComplex)
	fill(DataRS{DataF1}, PSweep("freq", [1, 4, 16] .* 1e3)) do 𝑓
		𝜔 = 2π*𝑓
		T = 1/𝑓
		Δt = T/100 #Define resolution from # of samples per period
		Tsim = 4T #Simulated time

		t = DataF1(0:Δt:Tsim) #DataF1 creates a t:{y, x} container with y == x
		sig = A * sin(𝜔*t + 𝜙) #Still a DataF1 sig:{y, x=t} container
		return sig
end; end; end

return signal


text:
"""
Test the flexibility of the scanner instantiation

The scanner should be composable not only of devices defined within
MPIMeasurements, but also from the outside.
"""
@testset "Flexible scanner" begin
  ### This section can be in a different package, which is not coupled to MPIMeasurements.jl as a dependency
  export FlexibleDAQ, FlexibleDAQParams

  Base.@kwdef struct FlexibleDAQParams <: MPIMeasurements.DeviceParams
    samplesPerPeriod::Int
    sendFrequency::typeof(1u"kHz")
  end

  FlexibleDAQParams(dict::Dict) = params_from_dict(FlexibleDAQParams, dict)

  Base.@kwdef mutable struct FlexibleDAQ <: MPIMeasurements.AbstractDAQ
    "Unique device ID for this device as defined in the configuration."
    deviceID::String
    "Parameter struct for this devices read from the configuration."
    params::FlexibleDAQParams
    "Flag if the device is optional."
	  optional::Bool = false
    "Flag if the device is present."
	  present::Bool = false
    "Vector of dependencies for this device."
    dependencies::Dict{String, Union{Device, Missing}}
  end

  function MPIMeasurements.init(daq::FlexibleDAQ)
    @debug "Initializing flexible DAQ with ID `$(daq.deviceID)`."

    daq.present = true
  end

  MPIMeasurements.checkDependencies(daq::FlexibleDAQ) = true

  Base.close(daq::FlexibleDAQ) = nothing

  function MPIMeasurements.startTx(daq::FlexibleDAQ)
  end

  function MPIMeasurements.stopTx(daq::FlexibleDAQ)
  end

  function MPIMeasurements.setTxParams(daq::FlexibleDAQ, amplitude, phase; postpone=false)
  end

  function MPIMeasurements.currentFrame(daq::FlexibleDAQ)
      return 1;
  end

  function MPIMeasurements.currentPeriod(daq::FlexibleDAQ)
      return 1;
  end

  function MPIMeasurements.readData(daq::FlexibleDAQ, startFrame, numFrames)
    uMeas=zeros(2,2,2,2)
    uRef=zeros(2,2,2,2)
    return uMeas, uRef
  end

  function MPIMeasurements.readDataPeriods(daq::FlexibleDAQ, startPeriod, numPeriods)
    uMeas=zeros(2,2,2,2)
    uRef=zeros(2,2,2,2)
    return uMeas, uRef
  end

  ### / External section

  scannerName_ = "TestFlexibleScanner"
  scanner = MPIScanner(scannerName_)

  @testset "Meta" begin
    @test name(scanner) == scannerName_
    @test configDir(scanner) == joinpath(testConfigDir, scannerName_)
  end

  @testset "General" begin
    generalParams_ = generalParams(scanner)
    @test generalParams_ isa MPIScannerGeneral
    @test generalParams_.boreSize == 1337u"mm"
    @test generalParams_.facility == "My awesome institute"
    @test generalParams_.manufacturer == "Me, Myself and I"
    @test generalParams_.name == scannerName_
    @test generalParams_.topology == "FFL"
    @test generalParams_.gradient == 42u"T/m"
    @test scannerBoreSize(scanner) == 1337u"mm"
    @test scannerFacility(scanner) == "My awesome institute"
    @test scannerManufacturer(scanner) == "Me, Myself and I"
    @test scannerName(scanner) == scannerName_
    @test scannerTopology(scanner) == "FFL"
    @test scannerGradient(scanner) == 42u"T/m"
  end

  @testset "Devices" begin

    @testset "DAQ" begin
      daq = getDevice(scanner, "my_daq_id")
      @test daq isa FlexibleDAQ
      @test daq.params.samplesPerPeriod == 1000
      @test daq.params.sendFrequency == 25u"kHz"
    end
  end

  close(scanner)
end


text:
# export NodeNil, is_nil

# import .Nodes.is_nil

# This node is used mostly for tests
mutable struct TrivialNode <: Ranger.AbstractNode
    id::UInt32
    name::String

    # NOTE A test field
    transform::Nodes.TransformProperties{Float64}

    function TrivialNode()
        o = new()
        o.id = 0
        o.name = "TrivialNode"
        o.transform = Nodes.TransformProperties{Float64}()
        o
    end
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
module JuliaAPI

using Dates

using PyCall
using PyCall: pyjlwrap_new
using Requires

using ..PyBase

const NpyNumber = Union{
    (t for t in values(PyCall.npy_typestrs) if t <: Number)...
}

_wrap(obj::Any) = pyjlwrap_new(obj)
# Wrap the object if it is desirable to _not_ invoke PyCall's
# automatic conversion.

_wrap(obj::Union{
    Nothing,
    Integer,
    NpyNumber,            # to exclude ForwardDiff.Dual etc.
    Array{<: NpyNumber},  # ditto
    AbstractString,  # should it be just String?
    Dates.AbstractTime,
    IO,
}) = obj
# It's OK to include some types that are not supported PyCall.  In
# that case, those objects are passed through pyjlwrap_new anyway.
# What should be avoided here is to include some types that would be
# converted by PyCall in an irreversible manner (e.g., Symbol,
# BitArray, etc.).

"""
    wrapcall(f, args...; kwargs...)

Wrap what `f(args...; kwargs...)` returns.
"""
function wrapcall(f, args...; kwargs...)
    return _wrap(f(args...; kwargs...))
    #=
    _f = _unwrap(f)
    _args = _unwrap.(args)
    _kwargs = (k => _unwrap(v) for (k, v) in kwargs)
    return _wrap(_f(_args...; _kwargs...))
    =#
end

struct WrappingCallback
    o::PyObject
    force::Bool
end

function (f::WrappingCallback)(args...; kwargs...)
    wrap = f.force ? pyjlwrap_new : _wrap
    return f.o(wrap.(args)...; (k => wrap(v) for (k, v) in kwargs)...)
end

struct _jlwrap_type end  # a type that would be wrapped as jlwrap by PyCall

get_jlwrap_prototype() = _jlwrap_type()

pybroadcast(op::String, args...) =
    _pybroadcast(eval(Meta.parse(strip(op))), args...)
pybroadcast(op, args...) = _pybroadcast(op, args...)
_pybroadcast(op, args::AbstractArray...) = op.(args...)
_pybroadcast(op, args...) = op(args...)

struct PySlice
    start
    stop
    step
end

struct PyEllipsis end

PyBase.convert_itemkey(self, key) = convert_pyindices(self, key)

function convert_pyindices(x, indices)
    if applicable(firstindex, x, 1)
        # arrays
        return convert_pyindex.((x,),
                                process_pyindices(ndims(x), indices),
                                1:length(indices))
    elseif indices isa Tuple{Integer} && applicable(firstindex, x)
        # tuples, etc.
        return (firstindex(x) + indices[1],)
    else
        # dictionaries, etc.
        return indices
    end
end

function process_pyindices(nd, indices)
    colons = repeat([:], inner=nd - length(indices))
    i = findfirst(i -> i isa PyEllipsis, indices)
    if i === nothing
        return (indices..., colons...)
    else
        return (indices[1:i - 1]..., colons..., ndices[i + 1:end]...)
    end
end

convert_pyindex(x, i, ::Int) = i

convert_pyindex(x, i::Integer, d::Int) = firstindex(x, d) + i

function convert_pyindex(x, slice::PySlice, d::Int)
    start = (slice.start === nothing ? 0 : slice.start) + firstindex(x, d)
    stop = slice.stop === nothing ? lastindex(x, d) : slice.stop + 1
    step = slice.step === nothing ? 1 : slice.step
    if step == 1
        return start:stop
    else
        return start:step:stop
    end
end

function eval_str(code::AbstractString;
                  scope::Module = Main,
                  filename::AbstractString = "string",
                  auto_jlwrap = true,
                  force_jlwrap = false)
    result = include_string(scope, code, filename)
    if force_jlwrap
        return pyjlwrap_new(result)
    elseif auto_jlwrap
        return _wrap(result)
    end
    return result
end

"""
    wrap_pycallable(x)

Wrap a Python callable `x` as a Julia callable which tries to box
arguments whenever automatic-conversion _could_ be lossy.  If `x` is
not a `PyObject`, return `x` as-is.
"""
wrap_pycallable(x::PyObject) = jlfunction(x)
wrap_pycallable(x::Any) = x

const JuliaObject = PyNULL()
const jlfunction = PyNULL()

function __init__()
    # Setup Python package `pybase_jl`:
    pushfirst!(PyVector(pyimport("sys")["path"]), @__DIR__)
    pybase_jl = pyimport("pybase_jl")
    pybase_jl[:setup](eval_str, JuliaAPI)
    copy!(JuliaObject, pybase_jl[:JuliaObject])
    copy!(jlfunction, pybase_jl[:jlfunction])

    # Don't wrap JuliaPy wrappers
    @require Pandas="eadc2687-ae89-51f9-a5d9-86b5a6373a9c" @eval begin
        _wrap(obj::Pandas.PandasWrapped) = obj
    end
    @require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" @eval begin
        _wrap(obj::SymPy.SymbolicObject) = obj
    end
    @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" @eval begin
        _wrap(obj::PyPlot.Figure) = obj
    end
end

"""
    pytest([args::Cmd]; inprocess=true)
"""
function pytest(args=``; inprocess=true,
                cwd = @__DIR__,
                kwargs...)
    cd(cwd) do
        if inprocess
            pytest_inprocess(args; kwargs...)
        else
            pytest_cli(args)
        end
    end
end

function pytest_inprocess(args; check=true)
    @info `pytest $args`
    code = pyimport("pytest")[:main](collect(args))
    if !check
        return code
    end
    if code != 0
        error("$(`pytest $args`) failed with code $code")
    end
end

function pytest_cli(args)
    command = `$(PyCall.pyprogramname) -m pytest $args`
    @info command
    run(command)
end

end  # module

using .JuliaAPI: pytest, _wrap, wrap_pycallable


text:
using TestCodeCov
using Test

@testset "TestCodeCov.jl" begin
    for it in 1:100
        @test TestCodeCov._identity(it) == it
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1FlexPersistentVolumeSource <: SwaggerModel
    driver::Any # spec type: Union{ Nothing, String } # spec name: driver
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    options::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: options
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    secretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: secretRef

    function IoK8sApiCoreV1FlexPersistentVolumeSource(;driver=nothing, fsType=nothing, options=nothing, readOnly=nothing, secretRef=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1FlexPersistentVolumeSource, Symbol("driver"), driver)
        setfield!(o, Symbol("driver"), driver)
        validate_property(IoK8sApiCoreV1FlexPersistentVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1FlexPersistentVolumeSource, Symbol("options"), options)
        setfield!(o, Symbol("options"), options)
        validate_property(IoK8sApiCoreV1FlexPersistentVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1FlexPersistentVolumeSource, Symbol("secretRef"), secretRef)
        setfield!(o, Symbol("secretRef"), secretRef)
        o
    end
end # type IoK8sApiCoreV1FlexPersistentVolumeSource

const _property_map_IoK8sApiCoreV1FlexPersistentVolumeSource = Dict{Symbol,Symbol}(Symbol("driver")=>Symbol("driver"), Symbol("fsType")=>Symbol("fsType"), Symbol("options")=>Symbol("options"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("secretRef")=>Symbol("secretRef"))
const _property_types_IoK8sApiCoreV1FlexPersistentVolumeSource = Dict{Symbol,String}(Symbol("driver")=>"String", Symbol("fsType")=>"String", Symbol("options")=>"Dict{String, String}", Symbol("readOnly")=>"Bool", Symbol("secretRef")=>"IoK8sApiCoreV1SecretReference")
Base.propertynames(::Type{ IoK8sApiCoreV1FlexPersistentVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1FlexPersistentVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1FlexPersistentVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1FlexPersistentVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1FlexPersistentVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1FlexPersistentVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1FlexPersistentVolumeSource)
    (getproperty(o, Symbol("driver")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1FlexPersistentVolumeSource }, name::Symbol, val)
end


text:
module Structures

include("mdl.jl")

mutable struct GenericStructure
    structure_type::String
    structure_role::String
    n_edges_total::Int64
    function GenericStructure(structure_role::String, n_edges_total::Int64)
        return new("generic", structure_role, n_edges_total)
    end
end

mutable struct Clique
    structure_type::String
    n_nodes_total::Int64
    n_edges_total::Int64
    nodes::Array{Int64,1}
    function Clique(n_edges_total::Int64, nodes::Array{Int64,1})
        length_nodes = length(nodes)
        @assert n_edges_total <= (length_nodes * (length_nodes - 1)) / 2 "Impossible number of edges! nv=$(length_nodes), ne=$(n_edges_total)"
        @assert length_nodes == length(unique(nodes)) "Input nodes must be unique!"
        return new("clique", length_nodes, n_edges_total, nodes)
    end
end

mutable struct Star
    structure_type::String
    n_nodes_total::Int64
    n_edges_total::Int64
    n_nodes_in_spokes::Int64
    n_edges_in_spokes::Int64
    hub::Int64
    spokes::Array{Int64,1}
    function Star(n_edges_total::Int64, hub::Int64, spokes::Array{Int64,1})
        length_spokes = length(spokes)
        @assert n_edges_total >= length_spokes "Too few edges!"
        @assert n_edges_total <= length_spokes + length_spokes * (length_spokes - 1) / 2 "Too many edges!"
        @assert length_spokes + 1 == length(unique([spokes...,hub])) "Input nodes must be unique!"
        return new("star", length_spokes + 1, n_edges_total, length_spokes, n_edges_total - length_spokes, hub, spokes)
    end
end

mutable struct Biclique
    structure_type::String
    n_nodes_total::Int64
    n_edges_total::Int64
    n_nodes_in_left::Int64
    n_nodes_in_right::Int64
    n_edges_in_left::Int64
    n_edges_in_right::Int64
    n_edges_across::Int64
    left_nodes::Array{Int64,1}
    right_nodes::Array{Int64,1}
    function Biclique(n_edges_total::Int64, left_nodes::Array{Int64,1}, right_nodes::Array{Int64,1}, n_edges_in_left, n_edges_in_right)
        left_length = length(left_nodes)
        right_length = length(right_nodes)
        n_edges_across = n_edges_total - n_edges_in_left - n_edges_in_right
        max_edges_in_left = left_length * (left_length - 1) / 2
        max_edges_in_right = right_length * (right_length - 1) / 2
        max_edges_across = left_length * right_length
        @assert n_edges_in_left <= max_edges_in_left "Too many edges in left!"
        @assert n_edges_in_right <= max_edges_in_right "Too many edges in right!"
        @assert n_edges_across <= max_edges_across "Too many edges in across!"
        @assert left_length + right_length == length(unique([left_nodes..., right_nodes...])) "Input nodes must be unique!"
        return new("biclique", left_length + right_length, n_edges_total, left_length, right_length,
            n_edges_in_left, n_edges_in_right, n_edges_across, left_nodes, right_nodes)
    end
end

mutable struct Starclique
    structure_type::String
    n_nodes_total::Int64
    n_edges_total::Int64
    n_nodes_in_left::Int64
    n_nodes_in_right::Int64
    n_edges_in_left::Int64
    n_edges_in_right::Int64
    n_edges_across::Int64
    left_nodes::Array{Int64,1}
    right_nodes::Array{Int64,1}
    function Starclique(n_edges_total::Int64, left_nodes::Array{Int64,1}, right_nodes::Array{Int64,1}, n_edges_in_left, n_edges_in_right)
        left_length = length(left_nodes)
        right_length = length(right_nodes)
        n_edges_across = n_edges_total - n_edges_in_left - n_edges_in_right
        max_edges_in_left = left_length * (left_length - 1) / 2
        max_edges_in_right = right_length * (right_length - 1) / 2
        max_edges_across = left_length * right_length
        @assert n_edges_in_left <= max_edges_in_left "Too many edges in left!"
        @assert n_edges_in_right <= max_edges_in_right "Too many edges in right!"
        @assert n_edges_across <= max_edges_across "Too many edges in across!"
        @assert left_length + right_length == length(unique([left_nodes..., right_nodes...])) "Input nodes must be unique!"
        return new("starclique", left_length + right_length, n_edges_total, left_length, right_length,
            n_edges_in_left, n_edges_in_right, n_edges_across, left_nodes, right_nodes)
    end
end

Structure = Union{GenericStructure, Clique, Star, Biclique, Starclique}
GenericBipartite = Union{Biclique, Starclique}

Base.show(io::IO, GS::GenericStructure) = print(io, "Generic Structure: $(GS.structure_role) with $(GS.n_edges_total) edges")
Base.show(io::IO, C::Clique) = print(io, "Clique with $(C.n_nodes_total) nodes and $(C.n_edges_total) edges")
Base.show(io::IO, S::Star) = print(io, "Star with $(S.n_nodes_in_spokes) spokes and $(S.n_edges_in_spokes) edges between spokes")
Base.show(io::IO, B::Biclique) = print(io, "Biclique with $(B.n_nodes_in_left) left nodes, $(B.n_nodes_in_right) right nodes, $(B.n_edges_in_left) edges in left, $(B.n_edges_in_right) edges in right, and $(B.n_edges_across) edges across")
Base.show(io::IO, B::Starclique) = print(io, "Starclique with $(B.n_nodes_in_left) left nodes, $(B.n_nodes_in_right) right nodes, $(B.n_edges_in_left) edges in left, $(B.n_edges_in_right) edges in right, and $(B.n_edges_across) edges across")

"""
generic structure description length
"""
function compute_description_length(structure::GenericStructure, maximum_n_loopy::Int64)
    @assert structure.structure_role in ["Overall number of edges", "Loop-freeness"] "Unexpected generic structure with no proper description length computation!"
    if structure.structure_role == "Loop-freeness"
        return 1.
    else
        return MDL.log2_zero(maximum_n_loopy)
    end
end

function compute_rest(clique::Clique)
    max_n_edges = convert(Int64, MDL.choose(clique.n_nodes_total, 2))
    cost_sparse_or_dense = 1
    cost_n_edges = MDL.log2_zero(MDL.log2_zero(floor(max_n_edges/2))) + MDL.log2_zero(min(clique.n_edges_total, max_n_edges - clique.n_edges_total))
    #cost_n_edges = MDL.universal_integer(min(clique.n_edges_total, max_n_edges - clique.n_edges_total) + 1)
    return cost_sparse_or_dense + cost_n_edges
end

### CLIQUES #######################################################################################################################################

"""
clique description length without connector or with path connector
(note that total_n_remaining is computed differently depending on the connector)
"""
function compute_description_length(clique::Clique, total_n_remaining::Int64, path::Bool; with_node_ids_cost::Bool=true)
    cost_n_nodes = path ? MDL.universal_integer(clique.n_nodes_total - 1) : MDL.universal_integer(clique.n_nodes_total)
    if with_node_ids_cost
        cost_node_ids = path ? MDL.log2_choose(total_n_remaining, clique.n_nodes_total - 1) : MDL.log2_choose(total_n_remaining, clique.n_nodes_total)
    else
        cost_node_ids = 0
    end
    rest = compute_rest(clique)
    return cost_n_nodes + cost_node_ids + rest
end

"""
clique description length with overlap connector
(note that total_n_remaining is computed differently depending on the connector)
"""
function compute_description_length(clique::Clique, parent_size::Int64, overlap_size::Int64, total_n_remaining::Int64)
    cost_n_nodes = MDL.universal_integer(clique.n_nodes_total - overlap_size)
    cost_node_ids = MDL.log2_choose(parent_size, overlap_size) + MDL.log2_choose(total_n_remaining, clique.n_nodes_total - overlap_size)
    rest = compute_rest(clique)
    return cost_n_nodes + cost_node_ids + rest
end

### STARS #########################################################################################################################################

"""
helper to remove redundancy in star description length computation
"""
function compute_rest(star::Star)
    cost_of_hub = MDL.log2_zero(star.n_nodes_total)
    cost_n_edges = MDL.log2_zero(MDL.log2_zero((star.n_nodes_in_spokes*(star.n_nodes_in_spokes-1)/2))) + MDL.log2_zero(star.n_edges_in_spokes)
    #cost_n_edges = MDL.universal_integer(star.n_edges_in_spokes + 1)
    return cost_of_hub + cost_n_edges
end

"""
star description length without connector or with path connector
(note that total_n_remaining is computed differently depending on the connector)
"""
function compute_description_length(star::Star, total_n_remaining::Int64, path::Bool; with_node_ids_cost::Bool=true)
    cost_n_spokes = path ? MDL.universal_integer(star.n_nodes_in_spokes - 1) : MDL.universal_integer(star.n_nodes_in_spokes)
    if with_node_ids_cost
        cost_node_ids = path ? MDL.log2_choose(total_n_remaining, star.n_nodes_total - 1) : MDL.log2_choose(total_n_remaining, star.n_nodes_total)
    else
        cost_node_ids = 0
    end
    rest = compute_rest(star)
    return cost_n_spokes + cost_node_ids + rest
end

"""
star description length with overlap connector
(note that total_n_remaining is computed differently depending on the connector)
"""
function compute_description_length(star::Star, parent_size::Int64, overlap_size::Int64, total_n_remaining::Int64)
    cost_n_spokes = MDL.universal_integer(star.n_nodes_in_spokes)
    cost_node_ids = MDL.log2_choose(parent_size, overlap_size) + MDL.log2_choose(total_n_remaining, star.n_nodes_total - overlap_size)
    rest = compute_rest(star)
    return cost_n_spokes + cost_node_ids + rest
end

### BICLIQUES AND STARCLIQUES #####################################################################################################################

"""
helper to remove redundancy in biclique description length computation
"""
function compute_rest(biclique::GenericBipartite)
    cost_n_left = MDL.log2_zero(biclique.n_nodes_total)
    cost_left_ids = MDL.log2_choose(biclique.n_nodes_total, biclique.n_nodes_in_left)
    # across, we definitely want density
    max_across_edges = biclique.n_nodes_in_left * biclique.n_nodes_in_right
    cost_n_across = MDL.log2_zero(MDL.log2_zero(max_across_edges)) + MDL.log2_zero(biclique.n_edges_across)
    max_left_edges = convert(Int64, MDL.choose(biclique.n_nodes_in_left, 2))
    max_right_edges = convert(Int64, MDL.choose(biclique.n_nodes_in_right, 2))
    # left we can have sparsity or density
    cost_n_left_edges = MDL.log2_zero(MDL.log2_zero(floor(max_left_edges/2))) + MDL.log2_zero(min(biclique.n_edges_in_left, max_left_edges - biclique.n_edges_in_left))
    # right we want sparsity
    cost_n_right_edges = MDL.log2_zero(MDL.log2_zero(floor(max_right_edges/2))) + MDL.log2_zero(biclique.n_edges_in_right)
    return (cost_n_left + cost_left_ids
            + cost_n_across + cost_n_left_edges + cost_n_right_edges)
end

"""
biclique description length without connector or with path connector
(note that total_n_remaining is computed differently depending on the connector)
"""
function compute_description_length(biclique::GenericBipartite, total_n_remaining::Int64, path::Bool; with_node_ids_cost::Bool=true)
    cost_n_nodes = path ? MDL.universal_integer(biclique.n_nodes_total - 1) : MDL.universal_integer(biclique.n_nodes_total)
    if with_node_ids_cost
        cost_node_ids = path ? MDL.log2_choose(total_n_remaining, biclique.n_nodes_total - 1) : MDL.log2_choose(total_n_remaining, biclique.n_nodes_total)
    else
        cost_node_ids = 0
    end
    rest = compute_rest(biclique)
    return (cost_n_nodes + cost_node_ids + rest)
end

"""
biclique description length with overlap connector
(note that total_n_remaining is computed differently depending on the connector)
"""
function compute_description_length(biclique::GenericBipartite, parent_size::Int64, overlap_size::Int64, total_n_remaining::Int64)
    cost_n_nodes = MDL.universal_integer(biclique.n_nodes_total - overlap_size)
    cost_node_ids = MDL.log2_choose(parent_size, overlap_size) + MDL.log2_choose(total_n_remaining, biclique.n_nodes_total - overlap_size)
    rest = compute_rest(biclique)
    return (cost_n_nodes + cost_node_ids + rest)
end

"""
get all nodes from a clique
"""
function get_node_set(clique::Clique)
    return clique.nodes
end

"""
get all nodes from a star
"""
function get_node_set(star::Star)
    return [star.hub, star.spokes...]
end

"""
get all nodes from a biclique
"""
function get_node_set(biclique::GenericBipartite)
    return [biclique.left_nodes..., biclique.right_nodes...]
end

end

max_stars_repo_name: dataspider/momo

max_issues_repo_name: dataspider/momo

max_forks_repo_name: dataspider/momo


text:
using HomotopyContinuation

@time begin
    @var x y
    affine_sqr =
        System([2.3 * x^2 + 1.2 * y^2 + 3x - 2y + 3, 2.3 * x^2 + 1.2 * y^2 + 5x + 2y - 5])
    res = solve(affine_sqr; compile = false)
end
res = solve(affine_sqr; compile = false, start_system = :total_degree)
@var x y z
proj_square = System([
    2.3 * x^2 + 1.2 * y^2 + 3x * z - 2y * z + (3 + 0im) * z^2,
    2.3 * x^2 + 1.2 * y^2 + 5x * z + 2y * z - (5 + 0im) * z^2,
])
solve(proj_square; compile = false)
solve(affine_sqr; compile = false, start_system = :total_degree)

@var x y
affine_ov = System([
    (x^2 + y^2 + x * y - 3) * (x + 3),
    (x^2 + y^2 + x * y - 3) * (y - x + 2.1),
    2x + 5y - 3,
])
solve(affine_ov; compile = false)

function toric_ed(A)
    d, n = size(A)
    @var t[1:d] y[1:n] u[1:n]

    φ = map(j -> prod(i -> t[i]^A[i, j], 1:d), 1:n)
    Dφ = [differentiate(φ[i], t[j]) for i = 1:n, j = 1:d]

    System([φ + y - u; Dφ' * y], parameters = u)
end

F = toric_ed([3 2 1 0; 0 1 2 3])
result = monodromy_solve(
    F,
    target_solutions_count = 21,
    max_loops_no_progress = 20,
    compile = false,
)

solve(
    F,
    solutions(result),
    start_parameters = parameters(result),
    target_parameters = randn(nparameters(F)),
    compile = false,
)


text:
using Test

using Crystallography
using QuantumESPRESSOBase.Inputs.PWscf
using QuantumESPRESSOParser.Inputs

@testset "Parse empty strings" begin
    @test_throws Meta.ParseError parse(ControlNamelist, " ")
    @test tryparse(ControlNamelist, " ") === nothing
    @test tryparse(ControlNamelist, "&control/") == ControlNamelist()
    @test parse(ControlNamelist, "&control\n/\n") == ControlNamelist()
    @test_throws Meta.ParseError parse(ControlNamelist, " ")
end

@testset "Test parsing a `PWInput`: silicon" begin
    # From https://github.com/QEF/q-e/blob/121edf5/PW/examples/example01/run_example#L81-L116
    str = read("data/si.scf.cg.in", String)
    pw = parse(PWInput, str)
    @test pw.control == ControlNamelist(;
        prefix = "silicon",
        tstress = true,
        tprnfor = true,
        pseudo_dir = "pseudo/",
        outdir = "./",
    )
    @test pw.system ==
          SystemNamelist(; ibrav = 2, celldm = [10.20], nat = 2, ntyp = 1, ecutwfc = 18)
    @test pw.electrons == ElectronsNamelist(; diagonalization = "cg", conv_thr = 1e-8)
    @test pw.ions == IonsNamelist()
    @test pw.cell == CellNamelist()
    @test pw.atomic_species ==
          AtomicSpeciesCard([AtomicSpecies("Si", 28.086, "Si.pz-vbc.UPF")])
    @test pw.atomic_positions == AtomicPositionsCard([
        AtomicPosition("Si", [0.0, 0.0, 0.0]),
        AtomicPosition("Si", [0.25, 0.25, 0.25]),
    ])
    @test pw.k_points == SpecialPointsCard(
        [
            0.125 0.125 0.125 1.0
            0.125 0.125 0.375 3.0
            0.125 0.125 0.625 3.0
            0.125 0.125 0.875 3.0
            0.125 0.375 0.375 3.0
            0.125 0.375 0.625 6.0
            0.125 0.375 0.875 6.0
            0.125 0.625 0.625 3.0
            0.375 0.375 0.375 1.0
            0.375 0.375 0.625 3.0
        ],
    )
    @test pw.cell_parameters === nothing
end

@testset "Test parsing a `PWInput`: silicon bands" begin
    # From https://github.com/QEF/q-e/blob/121edf5/PW/examples/example01/run_example#L125-L172
    str = read("data/si.band.cg.in", String)
    pw = parse(PWInput, str)
    @test pw.control == ControlNamelist(;
        calculation = "bands",
        pseudo_dir = "pseudo/",
        outdir = "./",
        prefix = "silicon",
    )
    @test pw.system == SystemNamelist(;
        ibrav = 2,
        celldm = [10.2],
        nat = 2,
        ntyp = 1,
        ecutwfc = 18.0,
        nbnd = 8,
    )
    @test pw.electrons == ElectronsNamelist(; diagonalization = "cg")
    @test pw.atomic_species ==
          AtomicSpeciesCard([AtomicSpecies("Si", 28.086, "Si.pz-vbc.UPF")])
    @test pw.atomic_positions == AtomicPositionsCard([
        AtomicPosition("Si", [0.0, 0.0, 0.0]),
        AtomicPosition("Si", [0.25, 0.25, 0.25]),
    ])
    @test pw.k_points == SpecialPointsCard(
        [
            0.0 0.0 0.0 1.0
            0.0 0.0 0.1 1.0
            0.0 0.0 0.2 1.0
            0.0 0.0 0.3 1.0
            0.0 0.0 0.4 1.0
            0.0 0.0 0.5 1.0
            0.0 0.0 0.6 1.0
            0.0 0.0 0.7 1.0
            0.0 0.0 0.8 1.0
            0.0 0.0 0.9 1.0
            0.0 0.0 1.0 1.0
            0.0 0.0 0.0 1.0
            0.0 0.1 0.1 1.0
            0.0 0.2 0.2 1.0
            0.0 0.3 0.3 1.0
            0.0 0.4 0.4 1.0
            0.0 0.5 0.5 1.0
            0.0 0.6 0.6 1.0
            0.0 0.7 0.7 1.0
            0.0 0.8 0.8 1.0
            0.0 0.9 0.9 1.0
            0.0 1.0 1.0 1.0
            0.0 0.0 0.0 1.0
            0.1 0.1 0.1 1.0
            0.2 0.2 0.2 1.0
            0.3 0.3 0.3 1.0
            0.4 0.4 0.4 1.0
            0.5 0.5 0.5 1.0
        ],
    )
end

@testset "Test parsing a `PWInput`: aluminium" begin
    # From https://github.com/QEF/q-e/blob/121edf5/PW/examples/example01/run_example#L186-L268
    str = read("data/al.scf.david.in", String)
    pw = parse(PWInput, str)
    @test pw.control == ControlNamelist(;
        calculation = "scf",
        restart_mode = "from_scratch",
        pseudo_dir = "pseudo/",
        outdir = "./",
        prefix = "al",
        tprnfor = true,
        tstress = true,
    )
    @test pw.system == SystemNamelist(;
        ibrav = 2,
        celldm = [7.50],
        nat = 1,
        ntyp = 1,
        ecutwfc = 15.0,
        occupations = "smearing",
        smearing = "marzari-vanderbilt",
        degauss = 0.05,
    )
    @test pw.electrons == ElectronsNamelist()
    @test pw.atomic_species ==
          AtomicSpeciesCard([AtomicSpecies("Al", 26.98, "Al.pz-vbc.UPF")])
    @test pw.atomic_positions ==
          AtomicPositionsCard([AtomicPosition("Al", [0.0, 0.0, 0.0])])
    @test pw.k_points == SpecialPointsCard(
        [
            0.0625000 0.0625000 0.0625000 1.00
            0.0625000 0.0625000 0.1875000 3.00
            0.0625000 0.0625000 0.3125000 3.00
            0.0625000 0.0625000 0.4375000 3.00
            0.0625000 0.0625000 0.5625000 3.00
            0.0625000 0.0625000 0.6875000 3.00
            0.0625000 0.0625000 0.8125000 3.00
            0.0625000 0.0625000 0.9375000 3.00
            0.0625000 0.1875000 0.1875000 3.00
            0.0625000 0.1875000 0.3125000 6.00
            0.0625000 0.1875000 0.4375000 6.00
            0.0625000 0.1875000 0.5625000 6.00
            0.0625000 0.1875000 0.6875000 6.00
            0.0625000 0.1875000 0.8125000 6.00
            0.0625000 0.1875000 0.9375000 6.00
            0.0625000 0.3125000 0.3125000 3.00
            0.0625000 0.3125000 0.4375000 6.00
            0.0625000 0.3125000 0.5625000 6.00
            0.0625000 0.3125000 0.6875000 6.00
            0.0625000 0.3125000 0.8125000 6.00
            0.0625000 0.3125000 0.9375000 6.00
            0.0625000 0.4375000 0.4375000 3.00
            0.0625000 0.4375000 0.5625000 6.00
            0.0625000 0.4375000 0.6875000 6.00
            0.0625000 0.4375000 0.8125000 6.00
            0.0625000 0.4375000 0.9375000 6.00
            0.0625000 0.5625000 0.5625000 3.00
            0.0625000 0.5625000 0.6875000 6.00
            0.0625000 0.5625000 0.8125000 6.00
            0.0625000 0.6875000 0.6875000 3.00
            0.0625000 0.6875000 0.8125000 6.00
            0.0625000 0.8125000 0.8125000 3.00
            0.1875000 0.1875000 0.1875000 1.00
            0.1875000 0.1875000 0.3125000 3.00
            0.1875000 0.1875000 0.4375000 3.00
            0.1875000 0.1875000 0.5625000 3.00
            0.1875000 0.1875000 0.6875000 3.00
            0.1875000 0.1875000 0.8125000 3.00
            0.1875000 0.3125000 0.3125000 3.00
            0.1875000 0.3125000 0.4375000 6.00
            0.1875000 0.3125000 0.5625000 6.00
            0.1875000 0.3125000 0.6875000 6.00
            0.1875000 0.3125000 0.8125000 6.00
            0.1875000 0.4375000 0.4375000 3.00
            0.1875000 0.4375000 0.5625000 6.00
            0.1875000 0.4375000 0.6875000 6.00
            0.1875000 0.4375000 0.8125000 6.00
            0.1875000 0.5625000 0.5625000 3.00
            0.1875000 0.5625000 0.6875000 6.00
            0.1875000 0.6875000 0.6875000 3.00
            0.3125000 0.3125000 0.3125000 1.00
            0.3125000 0.3125000 0.4375000 3.00
            0.3125000 0.3125000 0.5625000 3.00
            0.3125000 0.3125000 0.6875000 3.00
            0.3125000 0.4375000 0.4375000 3.00
            0.3125000 0.4375000 0.5625000 6.00
            0.3125000 0.4375000 0.6875000 6.00
            0.3125000 0.5625000 0.5625000 3.00
            0.4375000 0.4375000 0.4375000 1.00
            0.4375000 0.4375000 0.5625000 3.00
        ],
    )
end

@testset "Test parsing a `PWInput`: cobalt relaxation" begin
    # From https://github.com/QEF/q-e/blob/121edf5/PW/examples/example02/run_example#L83-L112
    str = read("data/co.rx.in", String)
    pw = parse(PWInput, str)
    @test pw.control == ControlNamelist(;
        calculation = "relax",
        pseudo_dir = "pseudo/",
        outdir = "./",
        prefix = "CO",
    )
    @test pw.system ==
          SystemNamelist(; ibrav = 0, nat = 2, ntyp = 2, ecutwfc = 24, ecutrho = 144)
    @test pw.electrons == ElectronsNamelist(; conv_thr = 1e-7)
    @test pw.ions == IonsNamelist()
    @test pw.atomic_species == AtomicSpeciesCard([
        AtomicSpecies("O", 1, "O.pz-rrkjus.UPF"),
        AtomicSpecies("C", 1.0, "C.pz-rrkjus.UPF"),
    ])
    @test pw.atomic_positions == AtomicPositionsCard(
        [AtomicPosition("C", [2.256, 0.0, 0.0]), AtomicPosition("O", [0, 0, 0], [0, 0, 0])],
        "bohr",
    )
    @test pw.k_points == GammaPointCard()
end

@testset "Test parsing a `PWInput`: FeO" begin
    # From https://github.com/QEF/q-e/blob/121edf5/PW/examples/example08/run_example#L84-L124
    str = read("data/feo_LDA.in", String)
    pw = parse(PWInput, str)
    @test pw.control == ControlNamelist(;
        calculation = "scf",
        restart_mode = "from_scratch",
        prefix = "feo_af",
        pseudo_dir = "pseudo/",
        outdir = "./",
        tstress = true,
        tprnfor = true,
    )
    @test pw.system == SystemNamelist(;
        ibrav = 0,
        celldm = [8.19],
        nat = 4,
        ntyp = 3,
        ecutwfc = 30.0,
        ecutrho = 240,
        nbnd = 20,
        starting_magnetization = [0, 0.5, -0.5],
        occupations = "smearing",
        smearing = "mv",
        degauss = 0.01,
        nspin = 2,
        # lda_plus_u = true,
        # Hubbard_U = [nothing, 1e-8, 1e-8],
    )
    @test pw.electrons ==
          ElectronsNamelist(; conv_thr = 1e-6, mixing_beta = 0.3, mixing_fixed_ns = 0)
    @test pw.ions == IonsNamelist()
    @test pw.cell_parameters == CellParametersCard([
        0.50 0.50 1.00
        0.50 1.00 0.50
        1.00 0.50 0.50
    ])
    @test pw.atomic_species == AtomicSpeciesCard([
        AtomicSpecies("O1", 1, "O.pz-rrkjus.UPF"),
        AtomicSpecies("Fe1", 1.0, "Fe.pz-nd-rrkjus.UPF"),
        AtomicSpecies("Fe2", 1.0, "Fe.pz-nd-rrkjus.UPF"),
    ])
    @test pw.atomic_positions == AtomicPositionsCard(
        [
            AtomicPosition("O1", [1, 1, 1] / 4),
            AtomicPosition("O1", [3, 3, 3] / 4),
            AtomicPosition("Fe1", [0, 0, 0]),
            AtomicPosition("Fe2", [0.5, 0.5, 0.5]),
        ],
        "crystal",
    )
    @test pw.k_points == KMeshCard(MonkhorstPackGrid([2, 2, 2], [0, 0, 0]))
end


text:
# # Polynomial Regression

# The following example features the LMO for polynomial regression on the ``\ell_1`` norm ball. Given input/output pairs ``\{x_i,y_i\}_{i=1}^N`` and sparse coefficients ``c_j``, where
# ```math
# y_i=\sum_{j=1}^m c_j f_j(x_i)
# ```
# and ``f_j: \mathbb{R}^n\to\mathbb{R}``, the task is to recover those ``c_j`` that are non-zero alongside their corresponding values. Under certain assumptions,
# this problem can be convexified into
# ```math
# \min_{c\in\mathcal{C}}||y-Ac||^2
# ```
# for a convex set ``\mathcal{C}``. It can also be found as example 4.1 [in the paper](https://arxiv.org/pdf/2104.06675.pdf).
# In order to evaluate the polynomial, we generate a total of 1000 data points
# ``\{x_i\}_{i=1}^N`` from the standard multivariate Gaussian, with which we will compute the output variables ``\{y_i\}_{i=1}^N``. Before
# evaluating the polynomial, these points will be contaminated with noise drawn from a standard multivariate Gaussian.
# We run the [`away_frank_wolfe`](@ref) and [`blended_conditional_gradient`](@ref) algorithms, and compare them to Projected Gradient Descent using a
# smoothness estimate. We will evaluate the output solution on test points drawn in a similar manner as the training points.

using FrankWolfe

using LinearAlgebra
import Random

using MultivariatePolynomials
using DynamicPolynomials

using Plots

using LaTeXStrings

const N = 10

DynamicPolynomials.@polyvar X[1:15]

const max_degree = 4
coefficient_magnitude = 10
noise_magnitude = 1

const var_monomials = MultivariatePolynomials.monomials(X, 0:max_degree)

Random.seed!(42)
const all_coeffs = map(var_monomials) do m
    d = MultivariatePolynomials.degree(m)
    return coefficient_magnitude * rand() .* (rand() .> 0.95 * d / max_degree)
end

const true_poly = dot(all_coeffs, var_monomials)

const training_data = map(1:500) do _
    x = 0.1 * randn(N)
    y = MultivariatePolynomials.subs(true_poly, Pair(X, x)) + noise_magnitude * randn()
    return (x, y.a[1])
end

const extended_training_data = map(training_data) do (x, y)
    x_ext = getproperty.(MultivariatePolynomials.subs.(var_monomials, X => x), :α)
    return (x_ext, y)
end

const test_data = map(1:1000) do _
    x = 0.4 * randn(N)
    y = MultivariatePolynomials.subs(true_poly, Pair(X, x)) + noise_magnitude * randn()
    return (x, y.a[1])
end

const extended_test_data = map(test_data) do (x, y)
    x_ext = getproperty.(MultivariatePolynomials.subs.(var_monomials, X => x), :α)
    return (x_ext, y)
end

function f(coefficients)
    return 0.5 / length(extended_training_data) * sum(extended_training_data) do (x, y)
        return (dot(coefficients, x) - y)^2
    end
end

function f_test(coefficients)
    return 0.5 / length(extended_test_data) * sum(extended_test_data) do (x, y)
        return (dot(coefficients, x) - y)^2
    end
end

function coefficient_errors(coeffs)
    return 0.5 * sum(eachindex(all_coeffs)) do idx
        return (all_coeffs[idx] - coeffs[idx])^2
    end
end

function grad!(storage, coefficients)
    storage .= 0
    for (x, y) in extended_training_data
        p_i = dot(coefficients, x) - y
        @. storage += x * p_i
    end
    storage ./= length(training_data)
    return nothing
end

function build_callback(trajectory_arr)
    return function callback(state)
        return push!(
            trajectory_arr,
            (Tuple(state)[1:5]..., f_test(state.x), coefficient_errors(state.x)),
        )
    end
end

gradient = similar(all_coeffs)

max_iter = 10000
random_initialization_vector = rand(length(all_coeffs))

lmo = FrankWolfe.LpNormLMO{1}(0.95 * norm(all_coeffs, 1))

## Estimating smoothness parameter
num_pairs = 1000
L_estimate = -Inf
gradient_aux = similar(gradient)

for i in 1:num_pairs # hide
    global L_estimate # hide
    x = compute_extreme_point(lmo, randn(size(all_coeffs))) # hide
    y = compute_extreme_point(lmo, randn(size(all_coeffs))) # hide
    grad!(gradient, x) # hide
    grad!(gradient_aux, y) # hide
    new_L = norm(gradient - gradient_aux) / norm(x - y) # hide
    if new_L > L_estimate # hide
        L_estimate = new_L # hide
    end # hide
end # hide

function projnorm1(x, τ)
    n = length(x)
    if norm(x, 1) ≤ τ
        return x
    end
    u = abs.(x)
    ## simplex projection
    bget = false
    s_indices = sortperm(u, rev=true)
    tsum = zero(τ)

    @inbounds for i in 1:n-1
        tsum += u[s_indices[i]]
        tmax = (tsum - τ) / i
        if tmax ≥ u[s_indices[i+1]]
            bget = true
            break
        end
    end
    if !bget
        tmax = (tsum + u[s_indices[n]] - τ) / n
    end

    @inbounds for i in 1:n
        u[i] = max(u[i] - tmax, 0)
        u[i] *= sign(x[i])
    end
    return u
end
xgd = FrankWolfe.compute_extreme_point(lmo, random_initialization_vector) # hide
training_gd = Float64[] # hide
test_gd = Float64[] # hide
coeff_error = Float64[] # hide
time_start = time_ns() # hide
gd_times = Float64[] # hide
for iter in 1:max_iter # hide
    global xgd # hide
    grad!(gradient, xgd) # hide
    xgd = projnorm1(xgd - gradient / L_estimate, lmo.right_hand_side) # hide
    push!(training_gd, f(xgd)) # hide
    push!(test_gd, f_test(xgd)) # hide
    push!(coeff_error, coefficient_errors(xgd)) # hide
    push!(gd_times, (time_ns() - time_start) * 1e-9) # hide
end # hide

x00 = FrankWolfe.compute_extreme_point(lmo, random_initialization_vector) # hide
x0 = deepcopy(x00) # hide

trajectory_lafw = [] # hide
callback = build_callback(trajectory_lafw) # hide
x_lafw, v, primal, dual_gap, _ = FrankWolfe.away_frank_wolfe( # hide
    f, # hide
    grad!, # hide
    lmo, # hide
    x0, # hide
    max_iteration=max_iter, # hide
    line_search=FrankWolfe.Adaptive(), # hide
    print_iter=max_iter ÷ 10, # hide
    emphasis=FrankWolfe.memory, # hide
    verbose=false, # hide
    lazy=true, # hide
    gradient=gradient, # hide
    callback=callback, # hide
    L=L_estimate, # hide
) # hide

trajectory_bcg = [] # hide
callback = build_callback(trajectory_bcg) # hide
x0 = deepcopy(x00) # hide
x_bcg, v, primal, dual_gap, _ = FrankWolfe.blended_conditional_gradient( # hide
    f, # hide
    grad!, # hide
    lmo, # hide
    x0, # hide
    max_iteration=max_iter, # hide
    line_search=FrankWolfe.Adaptive(), # hide
    print_iter=max_iter ÷ 10, # hide
    emphasis=FrankWolfe.memory, # hide
    verbose=false, # hide
    weight_purge_threshold=1e-10, # hide
    callback=callback, # hide
    L=L_estimate, # hide
) # hide
x0 = deepcopy(x00) # hide
trajectory_lafw_ref = [] # hide
callback = build_callback(trajectory_lafw_ref) # hide
_, _, primal_ref, _, _ = FrankWolfe.away_frank_wolfe( # hide
    f, # hide
    grad!, # hide
    lmo, # hide
    x0, # hide
    max_iteration=2 * max_iter, # hide
    line_search=FrankWolfe.Adaptive(), # hide
    print_iter=max_iter ÷ 10, # hide
    emphasis=FrankWolfe.memory, # hide
    verbose=false, # hide
    lazy=true, # hide
    gradient=gradient, # hide
    callback=callback, # hide
    L=L_estimate, # hide
) # hide

for i in 1:num_pairs
    global L_estimate
    x = compute_extreme_point(lmo, randn(size(all_coeffs)))
    y = compute_extreme_point(lmo, randn(size(all_coeffs)))
    grad!(gradient, x)
    grad!(gradient_aux, y)
    new_L = norm(gradient - gradient_aux) / norm(x - y)
    if new_L > L_estimate
        L_estimate = new_L
    end
end

# We can now perform projected gradient descent:

xgd = FrankWolfe.compute_extreme_point(lmo, random_initialization_vector)
training_gd = Float64[]
test_gd = Float64[]
coeff_error = Float64[]
time_start = time_ns()
gd_times = Float64[]
for iter in 1:max_iter
    global xgd
    grad!(gradient, xgd)
    xgd = projnorm1(xgd - gradient / L_estimate, lmo.right_hand_side)
    push!(training_gd, f(xgd))
    push!(test_gd, f_test(xgd))
    push!(coeff_error, coefficient_errors(xgd))
    push!(gd_times, (time_ns() - time_start) * 1e-9)
end

x00 = FrankWolfe.compute_extreme_point(lmo, random_initialization_vector)
x0 = deepcopy(x00)

trajectory_lafw = []
callback = build_callback(trajectory_lafw)
x_lafw, v, primal, dual_gap, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=max_iter,
    line_search=FrankWolfe.Adaptive(),
    print_iter=max_iter ÷ 10,
    emphasis=FrankWolfe.memory,
    verbose=false,
    lazy=true,
    gradient=gradient,
    callback=callback,
    L=L_estimate,
)

trajectory_bcg = []
callback = build_callback(trajectory_bcg)

x0 = deepcopy(x00)
x_bcg, v, primal, dual_gap, _ = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=max_iter,
    line_search=FrankWolfe.Adaptive(),
    print_iter=max_iter ÷ 10,
    emphasis=FrankWolfe.memory,
    verbose=false,
    weight_purge_threshold=1e-10,
    callback=callback,
    L=L_estimate,
)

x0 = deepcopy(x00)

trajectory_lafw_ref = []
callback = build_callback(trajectory_lafw_ref)
_, _, primal_ref, _, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=2 * max_iter,
    line_search=FrankWolfe.Adaptive(),
    print_iter=max_iter ÷ 10,
    emphasis=FrankWolfe.memory,
    verbose=false,
    lazy=true,
    gradient=gradient,
    callback=callback,
    L=L_estimate,
)

iteration_list = [
    [x[1] + 1 for x in trajectory_lafw],
    [x[1] + 1 for x in trajectory_bcg],
    collect(eachindex(training_gd)),
]
time_list = [
    [x[5] for x in trajectory_lafw],
    [x[5] for x in trajectory_bcg],
    gd_times,
]
primal_list = [
    [x[2] - primal_ref for x in trajectory_lafw],
    [x[2] - primal_ref for x in trajectory_bcg],
    [x - primal_ref for x in training_gd],
]
test_list = [
    [x[6] for x in trajectory_lafw],
    [x[6] for x in trajectory_bcg],
    test_gd,
]
label = [L"\textrm{L-AFW}", L"\textrm{BCG}", L"\textrm{GD}"]
coefficient_error_values = [
    [x[7] for x in trajectory_lafw],
    [x[7] for x in trajectory_bcg],
    coeff_error,
]

FrankWolfe.plot_results(
    [primal_list, primal_list, test_list, test_list],
    [iteration_list, time_list, iteration_list, time_list],
    label,
    [L"\textrm{Iteration}", L"\textrm{Time}", L"\textrm{Iteration}", L"\textrm{Time}"],
    [L"\textrm{Primal Gap}", L"\textrm{Primal Gap}", L"\textrm{Test loss}", L"\textrm{Test loss}"],
    xscalelog=[:log, :identity, :log, :identity],
    legend_position=[:bottomleft, nothing, nothing, nothing],
)


text:
"""
    read_fasta(filename; seqtype=String)

Read a fasta file, returning sequences as the specified type.

# Examples
```julia-repl
seqnames, seqs = read_fasta(filename)
````
"""
function read_fasta(filename; seqtype=String)
    stream = open(FASTA.Reader, filename)
    seqnames, seqs = String[], seqtype[]
    for entry in stream
        push!(seqs, FASTA.sequence(seqtype, entry))
        push!(seqnames, FASTA.identifier(entry))
    end
    return seqnames, seqs
end

"""
    read_fastq_records(filename; seqtype=String)
Read .fastq file contents.

# Examples
```julia-repl
seqs, phreds, seqnames = read_fastq(filename)
```
"""
function read_fastq(filename; seqtype=String)
    stream = open(FASTQ.Reader, filename)
    seqs, phreds, seqnames = seqtype[], Array{Int8, 1}[], String[]
    for entry in stream
        push!(seqs, FASTQ.sequence(seqtype, entry))
        push!(phreds, FASTQ.quality(entry, :sanger))
        push!(seqnames, FASTQ.identifier(entry))
    end
    return seqs, phreds, seqnames
end

"""
    clean_nt(seq::String; to_strip = ['.', '-'])
Enforce case, remove gaps and convert RNA to DNA in a single step. Alternatively use Biosequences types for more control.

# Examples
```julia-repl
degapped_dna = clean_nt("AUG---GUGUAG")
"ATGGTGTAG"
```
"""
function clean_nt(seq::String; to_strip = ['.', '-'])
    dna_seq = replace(uppercase(seq), 'U' => 'T')
    for c in to_strip
        dna_seq = replace(dna_seq, c => "")
    end
    return dna_seq
end

"""
    clean_nt(seq::String; to_strip = ['.', '-'])
Enforce case and remove gaps in a single step. Alternatively use Biosequences types for more control.

# Examples
```julia-repl
degapped_aa = clean_aa("mvk--lvss*")
"MVKLVSS*"
```
"""
function clean_aa(seq::String; to_strip = ['.', '-'])
    aa_seq = uppercase(seq)
    for c in to_strip
        aa_seq = replace(aa_seq, c => "")
    end
    return aa_seq
end


text:
using Test, GraphUtilities, GenericTensorNetworks
using GraphUtilities: parsegraph, parseproblem, RegularGraphConfig, to_toml, from_toml, parseproperty

@testset "parse config" begin
    gc = RegularGraphConfig(; size=10, degree=3)
    gp = GraphProblemConfig(; problem="IndependentSet", graph=gc)
    to_toml("_test.toml", gp)
    @test gp isa GraphProblemConfig
    @test parsegraph(gp.graph).ne == 15
    @test parseproblem(gp.problem) === IndependentSet
    gp2 = from_toml(GraphProblemConfig{typeof(gc)}, "_test.toml")
    @test gp == gp2
    @test hash(gp) == hash(gp2)
end

@testset "parse property" begin
    @test parseproperty("SizeMax2") == SizeMax(2)
    @test parseproperty("SizeMax") == SizeMax()
    @test parseproperty("ConfigsMax2") == ConfigsMax(2)
    @test parseproperty("ConfigsMin2") == ConfigsMin(2)
    @test parseproperty("SingleConfigMax3") == SingleConfigMax(3)
    @test parseproperty("SingleConfigMin3") == SingleConfigMin(3)
    @test parseproperty("ConfigsAll") == ConfigsAll()
    @test parseproperty("GraphPolynomial") == GraphPolynomial()
    @test parseproperty("CountingMax") == CountingMax()
    @test parseproperty("CountingMin4") == CountingMin(4)
    @test parseproperty("CountingAll") == CountingAll()
    @test parseproperty("ConfigsMaxTree3") == ConfigsMax(3; tree_storage=true)
    @test parseproperty("ConfigsMinTree3") == ConfigsMin(3; tree_storage=true)
    @test parseproperty("ConfigsAllTree") == ConfigsAll(; tree_storage=true)
end


text:
using NBodyIPs, StaticArrays, LinearAlgebra
Sob = NBodyIPs.Sobol

r0 = 1.0
r1 = 2.0
x0 = r0 * (@SVector ones(3))
x1 = r1 * (@SVector ones(3))

Xbl = Sob.filtered_sobol(x0, x1, Sob.bl_is_simplex,
                      npoints=1_000, nfailed=100_000)

const J3 = SVector(2, 3)
cart2r(R) = NBodyIPs.edge_lengths(R, J3)
filt(r) = (r0 <= minimum(r)) && (maximum(r) <= r1)
Xcart = Sob.filtered_cart_sobol(r1, 3, cart2r, filt,
                                npoints = 1_000, nfailed = 100_000
                       )

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
"""
    Hashpipe Module

C Hashpipe functions for the High-Availablity Shared Pipeline Engine ported for Julia usability with pure-Julian extensions.
Written by Max Hawkins
Hashpipe C code written by Dave MacMahon: https://github.com/david-macmahon/hashpipe
"""
module Hashpipe

# Hashpipe error Codes
"""
    HASHPIPE_OK

Status code signifying successful operation.
"""
const HASHPIPE_OK         =  0
"""
    HASHPIPE_TIMEOUT

Status code signifying call timeout.
"""
const HASHPIPE_TIMEOUT    =  1
"""
    HASHPIPE_ERR_GEN

Status code signifying an unknown error?
"""
const HASHPIPE_ERR_GEN    = -1 # Super non-informative???
"""
    HASHPIPE_ERR_SYS

Status code signifying a failed system call.
"""
const HASHPIPE_ERR_SYS    = -2
"""
    HASHPIPE_ERR_PARAM

Status code signifying parameter out of valid range.
"""
const HASHPIPE_ERR_PARAM  = -3
"""
    HASHPIPE_ERR_KEY

Status code signifying requested key doesn't exist.
"""
const HASHPIPE_ERR_KEY    = -4
"""
    HASHPIPE_ERR_PACKET

Status code signifying unexpected packet size.
"""
const HASHPIPE_ERR_PACKET = -5

# Inlcude libhashpipe and libhashpipestatus paths for C-calls
const deps_file = joinpath(dirname(@__FILE__),"..","deps","deps.jl")
if isfile(deps_file)
    include(deps_file)
else
    error("Hashpipe library not properly installed. Have you built Hashpipe? Try Pkg.build(\"Hashpipe\")")
end

include("status.jl")
include("databuf.jl")
include("fitshead.jl")
include("thread.jl")

end # Module Hashpipe


text:
# Dataset files

# This function generates the file names for lookup table files,
# atmospheric state files and reference data file for comparing
# results for test cases.
# The following options are currently supported:
# ftype ∈ (:lookup_tables, :atmos_state, :comparison)
# optics_type ∈ (:clearsky, :cloudy)
# opc ∈ (:OneScalar, :TwoStream, nothing)
# λ ∈ (:lw, :sw, nothing)
# flux_up_dn ∈ (:flux_up, :flux_dn, nothing)

function get_ref_filename(
    ftype,
    optics_type;
    λ = nothing,
    opc = nothing,
    flux_up_dn = nothing,
    lut_pade = nothing,
)
    @assert ftype ∈ (:lookup_tables, :atmos_state, :comparison) &&
            optics_type ∈ (:clearsky, :cloudysky)

    fname = String(optics_type)

    if ftype == :lookup_tables
        @assert λ ∈ (:lw, :sw)
        fname *= "_" * String(λ)
    elseif ftype == :atmos_state
        fname *= "_as"
    else #ftype == :comparison
        if λ isa Symbol
            @assert λ ∈ (:lw, :sw)
            fname *= "_" * String(λ)
        end

        if flux_up_dn isa Symbol
            @assert flux_up_dn ∈ (:flux_up, :flux_dn)
            fname *= "_" * String(flux_up_dn)
        end
        if opc isa Symbol
            @assert opc ∈ (:OneScalar, :TwoStream)
            fname *= "_" * String(opc)
        end
        if lut_pade isa Symbol
            @assert lut_pade ∈ (:lut, :pade)
            fname *= "_" * String(lut_pade)
        end
    end

    return joinpath(
        artifact"RRTMGPReferenceData",
        "RRTMGPReferenceData",
        String(ftype),
        fname * ".nc",
    )
end


text:
struct LoopMulFunc{P,TC,TA,TB,Α,Β,Md,Kd,Nd} <: Function end
function (::LoopMulFunc{P,TC,TA,TB,Α,Β,Md,Kd,Nd})(p::Ptr{UInt}) where {P,TC,TA,TB,Α,Β,Md,Kd,Nd}
  offset, C = load(p, TC, 2*sizeof(UInt))
  offset, A = load(p, TA, offset)
  offset, B = load(p, TB, offset)
  offset, α = load(p, Α, offset)
  offset, β = load(p, Β, offset)
  offset, M = load(p, Md, offset)
  offset, K = load(p, Kd, offset)
  offset, N = load(p, Nd, offset)
  _call_loopmul!(C, A, B, α, β, M, K, N, Val{P}())
  _atomic_store!(p, SPIN)
  nothing
end
@inline _call_loopmul!(C, A, B, α, β, M, K, N, ::Val{false}) = loopmul!(C, A, B, α, β, M, K, N)
@inline function _call_loopmul!(C::StridedPointer{T}, A, B, α, β, M, K, N, ::Val{true}) where {T}
  if M*K < first_cache_size(Val(T)) * R₂Default()
    packaloopmul!(C, A, B, α, β, M, K, N)
    return
  else
    matmul_st_only_pack_A!(C, A, B, α, β, M, K, N, W₁Default(), W₂Default(), R₁Default(), R₂Default())
    return
  end
end
call_loopmul!(C, A, B, α, β, M, K, N, ::Val{P}) where {P} = _call_loopmul!(C, A, B, α, β, M, K, N, Val{P}())

struct SyncMulFunc{TC,TA,TB,Α,Β,Md,Kd,Nd,BCP,ID,TT,W₁,W₂,R₁,R₂} <: Function end
function (::SyncMulFunc{TC,TA,TB,Α,Β,Md,Kd,Nd,BCP,ID,TT,W₁,W₂,R₁,R₂})(p::Ptr{UInt}) where {TC,TA,TB,Α,Β,Md,Kd,Nd,BCP,ID,TT,W₁,W₂,R₁,R₂}
  offset, C = load(p, TC, 2*sizeof(UInt))
  offset, A = load(p, TA, offset)
  offset, B = load(p, TB, offset)
  offset, α = load(p, Α, offset)
  offset, β = load(p, Β, offset)
  offset, M = load(p, Md, offset)
  offset, K = load(p, Kd, offset)
  offset, N = load(p, Nd, offset)
  offset, atomicp = load(p, Ptr{UInt32}, offset)
  offset, bcachep = load(p, BCP, offset)
  offset, id = load(p, ID, offset)
  offset, total_ids = load(p, TT, offset)
  sync_mul!(C, A, B, α, β, M, K, N, atomicp, bcachep, id, total_ids, StaticFloat64{W₁}(), StaticFloat64{W₂}(), StaticFloat64{R₁}(), StaticFloat64{R₂}())
  _atomic_store!(p, SPIN)
  nothing
end

@generated function cfuncpointer(::T) where {T}
  precompile(T(), (Ptr{UInt},))
  quote
    $(Expr(:meta,:inline))
    @cfunction($(T()), Cvoid, (Ptr{UInt},))
  end
end

@inline function setup_matmul!(p::Ptr{UInt}, C::TC, A::TA, B::TB, α::Α, β::Β, M::Md, K::Kd, N::Nd, ::Val{P}) where {P,TC,TA,TB,Α,Β,Md,Kd,Nd}
  offset = store!(p, cfuncpointer(LoopMulFunc{P,TC,TA,TB,Α,Β,Md,Kd,Nd}()), sizeof(UInt))
  offset = store!(p, C, offset)
  offset = store!(p, A, offset)
  offset = store!(p, B, offset)
  offset = store!(p, α, offset)
  offset = store!(p, β, offset)
  offset = store!(p, M, offset)
  offset = store!(p, K, offset)
  offset = store!(p, N, offset)
  nothing
end

@inline function launch_thread_mul!(C, A, B, α, β, M, K, N, tid::UInt32, ::Val{P}) where {P}
  launch(setup_matmul!, tid, C, A, B, α, β, M, K, N, Val{P}())
end

struct SyncMulLauncher{W₁, W₂, R₁, R₂} end
@inline function (::SyncMulLauncher{W₁, W₂, R₁, R₂})(
  p::Ptr{UInt}, C::TC, A::TA, B::TB, α::Α, β::Β, M::Md, K::Kd, N::Nd,
  ap::Ptr{UInt32},bcp::BCP,id::ID,tt::TT
) where {TC,TA,TB,Α,Β,Md,Kd,Nd,BCP,ID,TT,W₁,W₂,R₁,R₂}
  fptr = cfuncpointer(SyncMulFunc{TC,TA,TB,Α,Β,Md,Kd,Nd,BCP,ID,TT,W₁,W₂,R₁,R₂}())
  offset = store!(p, fptr, sizeof(UInt))
  offset = store!(p, C, offset)
  offset = store!(p, A, offset)
  offset = store!(p, B, offset)
  offset = store!(p, α, offset)
  offset = store!(p, β, offset)
  offset = store!(p, M, offset)
  offset = store!(p, K, offset)
  offset = store!(p, N, offset)
  offset = store!(p, ap,  offset)
  offset = store!(p, bcp, offset)
  offset = store!(p, id,  offset)
  offset = store!(p, tt,  offset)
  nothing
end
@inline function launch_thread_mul!(
  C, A, B, α, β, M, K, N, ap, bcp, tid, id, tt,
  ::StaticFloat64{W₁},::StaticFloat64{W₂},::StaticFloat64{R₁},::StaticFloat64{R₂}
) where {W₁,W₂,R₁,R₂}
  launch(SyncMulLauncher{W₁, W₂, R₁, R₂}(), tid, C, A, B, α, β, M, K, N, ap, bcp, id, tt)
end


text:
@testset "Mvn Model" begin
    @testset "1-var" begin
        @test @timev !isnothing(one_var_test((n_vars, dist) -> MvnModel()))
    end
    @testset "2-var" begin
        @test @timev !isnothing(multi_var_test((n_vars, dist) -> MvnModel()))
    end
    @testset "3-var" begin
        @test @timev !isnothing(multi_var_test((n_vars, dist) -> MvnModel(), n_vars=3))
    end
end


text:
function dump_state_and_aux_init(dgngrp, currtime)
    @assert dgngrp.interpol !== nothing
end

function get_dims(dgngrp)
    if dgngrp.interpol !== nothing
        if dgngrp.interpol isa InterpolationBrick
            if Array ∈ typeof(dgngrp.interpol.x1g).parameters
                h_x1g = dgngrp.interpol.x1g
                h_x2g = dgngrp.interpol.x2g
                h_x3g = dgngrp.interpol.x3g
            else
                h_x1g = Array(dgngrp.interpol.x1g)
                h_x2g = Array(dgngrp.interpol.x2g)
                h_x3g = Array(dgngrp.interpol.x3g)
            end
            dims = OrderedDict("x" => h_x1g, "y" => h_x2g, "z" => h_x3g)
        elseif dgngrp.interpol isa InterpolationCubedSphere
            if Array ∈ typeof(dgngrp.interpol.rad_grd).parameters
                h_rad_grd = dgngrp.interpol.rad_grd
                h_lat_grd = dgngrp.interpol.lat_grd
                h_long_grd = dgngrp.interpol.long_grd
            else
                h_rad_grd = Array(dgngrp.interpol.rad_grd)
                h_lat_grd = Array(dgngrp.interpol.lat_grd)
                h_long_grd = Array(dgngrp.interpol.long_grd)
            end
            dims = OrderedDict(
                "rad" => h_rad_grd,
                "lat" => h_lat_grd,
                "long" => h_long_grd,
            )
        else
            error("Unsupported interpolation topology $(dgngrp.interpol)")
        end
    else
        error("Dump of non-interpolated data currently unsupported")
    end

    return dims
end

function dump_state_and_aux_collect(dgngrp, currtime)
    DA = CLIMA.array_type()
    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    FT = eltype(Q.data)
    bl = dg.balancelaw
    mpirank = MPI.Comm_rank(mpicomm)

    # filename (may also want to take out)
    dprefix = @sprintf(
        "%s_%s-%s-num%04d",
        dgngrp.out_prefix,
        dgngrp.name,
        Settings.starttime,
        dgngrp.num
    )
    statefilename = joinpath(Settings.output_dir, dprefix)
    auxfilename = joinpath(Settings.output_dir, "$(dprefix)_aux")

    statenames = flattenednames(vars_state(bl, FT))
    #statenames = ("rho", "rhou", "rhov", "rhow", "rhoe")
    auxnames = flattenednames(vars_aux(bl, FT))

    all_state_data = nothing
    all_aux_data = nothing
    if dgngrp.interpol !== nothing
        istate = DA(Array{FT}(undef, dgngrp.interpol.Npl, num_state(bl, FT)))
        interpolate_local!(dgngrp.interpol, Q.data, istate)

        if dgngrp.project
            if dgngrp.interpol isa InterpolationCubedSphere
                # TODO: get indices here without hard-coding them
                _ρu, _ρv, _ρw = 2, 3, 4
                project_cubed_sphere!(dgngrp.interpol, istate, (_ρu, _ρv, _ρw))
            else
                error("Can only project for InterpolationCubedSphere")
            end
        end

        all_state_data =
            accumulate_interpolated_data(mpicomm, dgngrp.interpol, istate)

        iaux = DA(Array{FT}(undef, dgngrp.interpol.Npl, num_aux(bl, FT)))
        interpolate_local!(dgngrp.interpol, dg.auxstate.data, iaux)

        all_aux_data =
            accumulate_interpolated_data(mpicomm, dgngrp.interpol, iaux)
    else
        error("Dump of non-interpolated data currently unsupported")
    end

    dims = get_dims(dgngrp)

    statevarvals = OrderedDict()
    for i in 1:num_state(bl, FT)
        statevarvals[statenames[i]] = all_state_data[:, :, :, i]
    end
    write_data(dgngrp.writer, statefilename, dims, statevarvals, currtime)

    auxvarvals = OrderedDict()
    for i in 1:num_aux(bl, FT)
        auxvarvals[auxnames[i]] = all_aux_data[:, :, :, i]
    end
    write_data(dgngrp.writer, auxfilename, dims, auxvarvals, currtime)

    return nothing
end

function dump_state_and_aux_fini(dgngrp, currtime) end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
"""
Allows the user to set a starting point for a grid in a PowerDynamics simulation.
Inputs:
        pd: Powergrid to be simulated
        v_vec: Vector of voltages at the nodes
        ω_vec: Vector of frequencies at the nodes
        θ_vec: Vector of phases at the nodes
Outputs:
        Startingpoint: Customized State of the grid
"""
using PowerDynamics: PowerGrid, State

function my_startingpoint(pd::PowerGrid,v_vec, ω_vec, θ_vec)
    numnodes = size(pd.nodes)[1]
    @assert size(v_vec)[1] == numnodes "Size of v_vec not the same as number of nodes."
    @assert size(ω_vec)[1] == numnodes "Size of ω_vec not the same as number of nodes."
    @assert size(θ_vec)[1] == numnodes "Size of θ_vec not the same as number of nodes."

    println("Warning! Manually setting the startingpoint might not result
    in a stable fixpoint.")

    startingpoint = State(pd, ones(systemsize(pd))) # just put in a bunch of ones for now
    for i in 1:numnodes
        startingpoint[i,:ω] = ω_vec[i]
        startingpoint[i,:v] = v_vec[i]
        startingpoint[i,:θ] = θ_vec[i]
    end
    return startingpoint
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
"""
    FixedValueArray{T, N} <: AbstractArray{T, N}

Implicit and arbitrarily-sized array with a common value of type `T` for all elements.

# Constructors
    FixedValueArray(val, dims)
    FixedValueArray(val, dims...)

Create a new matrix of the given dimensions `dims`, where each element is represented
by the given value `val`.

# Examples

```jldoctest; setup = :(using ImplicitArrays)
julia> FixedValueArray(0, (2, 5))
2×5 FixedValueArray{Int64,2}:
 0  0  0  0  0
 0  0  0  0  0

julia> FixedValueArray("(^_^) Hi!", 1, 2)
1×2 FixedValueArray{String,2}:
 "(^_^) Hi!"  "(^_^) Hi!"

julia> FixedValueArray(FixedValueArray(5.0, (1, 2)), (1, 3))
1×3 FixedValueArray{FixedValueArray{Float64,2},2}:
 [5.0 5.0]  [5.0 5.0]  [5.0 5.0]

julia> length(FixedValueArray(1, 1_000_000_000, 1_000_000_000))
1000000000000000000
```
"""
struct FixedValueArray{T, N} <: AbstractArray{T, N}
    "Value representing all matrix elements"
    val::T
    "Array dimensions"
    dims::NTuple{N, Int}
end

FixedValueArray(val::T, dims::Int...) where T = FixedValueArray(val, dims)

@inline Base.size(A::FixedValueArray{T, N}) where {T, N} = A.dims

@inline function Base.getindex(A::FixedValueArray{T, N}, i::Int) where {T, N}
    @boundscheck checkbounds(A, i)
    A.val
end

@inline function Base.getindex(A::FixedValueArray{T, N}, I::Vararg{Int, N}) where {T, N}
    @boundscheck checkbounds(A, I...)
    A.val
end

@inline function Base.setindex!(A::FixedValueArray{T, N}, ::Any, ::Int) where {T, N}
    _throw_canonical_error("setindex! not defined for ", typeof(A))
end

# This one is required to have setindex! at index [] throw an CanonicalIndexError instead
# of a BoundsError (which is used for all other indices).
@inline function Base.setindex!(A::FixedValueArray, ::Any, ::Vararg{Int, N}) where N
    _throw_canonical_error("setindex! not defined for ", typeof(A))
end


text:
#4-echelon supply chain with market demand at levels 1 and 2 of the chain

using Distributions
using InventoryManagement

##define network connectivity
adj_matrix = [0 1 0 0 0;
              0 0 1 0 0;
              0 0 0 1 0;
              0 0 0 0 1;
              0 0 0 0 0]
net = MetaDiGraph(adj_matrix) # 1 (raw) -> 2 (plant) -> 3 (plant) -> 4 (storage; with demand) -> 5 (retail)
materials = [:A, :B, :C]
bom = [0  0  0; # C -> B -> A
      -1  0  0;
       0 -1  0]
set_prop!(net, :materials, materials)
set_prop!(net, :bill_of_materials, bom)

##specify parameters, holding costs and capacity, market demands and penalty for unfilfilled demand
set_props!(net, 1, Dict(:initial_inventory => Dict(:C => Inf),
                        :inventory_capacity => Dict(:A => 0, :B => 0, :C => Inf)))
set_props!(net, 2, Dict(:initial_inventory => Dict(:C => 215), #initial inventory at plant
                        :inventory_capacity => Dict(:A => 0, :B => 0, :C => Inf),
                        :production_time => Dict(:B => 0)))
set_props!(net, 3, Dict(:initial_inventory => Dict(:B => 190), #initial inventory at plant
                        :inventory_capacity => Dict(:A => 0, :B => Inf, :C => 0),
                        :production_time => Dict(:A => 0)))

set_props!(net, 4, Dict(:initial_inventory => Dict(:A => 105), #initial inventory at storage
                        :inventory_capacity => Dict(:A => Inf, :B => 0, :C => 0),
                        :demand_distribution => Dict(:A => [5]), #distribution given
                        :demand_frequency => Dict(:A => 1))) #order every other day on average

set_props!(net, 5, Dict(:initial_inventory => Dict(:A => 60), #initial inventory at retail
                        :inventory_capacity => Dict(:A => Inf, :B => 0, :C => 0),
                        :demand_distribution => Dict(:A => [10]), #Normal Distribution double that of storage
                        :demand_frequency => Dict(:A => 1))) #order daily

##specify lead times
set_props!(net, 1, 2, Dict(:lead_time => Dict(:C => [0])))
set_props!(net, 2, 3, Dict(:lead_time => Dict(:B => [0])))
set_props!(net, 3, 4, Dict(:lead_time => Dict(:A => [0])))
set_props!(net, 4, 5, Dict(:lead_time => Dict(:A => [0])))

##define reorder policy parameters
policy_type = :sS #(s, S) policy
review_period = 1 #continuous review
policy_variable = :ech_position #variable tracked by policy
s = Dict((2,:C) => 215, (3,:B) => 190, (4,:A) => 165, (5,:A) => 60) #lower bound on inventory
S = Dict((2,:C) => 215, (3,:B) => 190, (4,:A) => 165, (5,:A) => 60) #base stock level

##create environment and run simulation with reorder policy
num_periods = 100
env = SupplyChainEnv(net, num_periods, backlog = true, reallocate = true)
simulate_policy!(env, s, S; policy_type, review_period, policy_variable)

##make plots
using DataFrames, StatsPlots
#profit
node_profit = groupby(env.profit, :node)
profit = transform(node_profit, :value => cumsum)
fig1 = @df profit plot(:period, :value_cumsum, group={Node = :node}, legend = :topleft,
                    xlabel="period", ylabel="cumulative profit")

#inventory position
inv_position = filter(i -> i.node > 1, env.ech_position)
fig2 = @df inv_position plot(:period, :level, group={Node = :node, Mat = :material}, linetype=:steppost, legend = :topright,
                    xlabel="period", ylabel="echelon position")